Rheolef  7.2
an efficient C++ finite element environment
 
Loading...
Searching...
No Matches
basis_fem_Pk_sherwin.cc
Go to the documentation of this file.
1
23#include "piola_fem_lagrange.h"
24#include "rheolef/rheostream.h"
25#include "equispaced.icc"
26#include "warburton.icc"
27#include "sherwin.icc"
28#include "eigen_util.h"
30
31namespace rheolef {
32using namespace std;
33
34// =========================================================================
35// basis members
36// =========================================================================
37template<class T>
41template<class T>
43 : base (sopt),
44 _degree(degree),
45 _alpha(1.0),
46 _beta(1.0)
47#ifdef TODO
48 ,
49 _value_ad(),
50 _work0_ad(),
51 _work1_ad(),
52 _work2_ad()
53 _work0(),
54 _work1(),
55 _work2()
56#endif // TODO
57{
58 base::_name = base::standard_naming (family_name(), degree, base::_sopt);
59warning_macro ("cstor...");
60 _alpha = _beta = 1; // 1 or 2, see SheKar-2005 p. 53
62
63 // piola FEM transformatiion:
65 base::_piola_fem.piola_fem<T>::base::operator= (new_macro(piola_fem_type));
66warning_macro ("cstor done");
67}
68template<class T>
69void
71{
73 degree(),
74 base::is_continuous(),
75 base::_ndof_on_subgeo_internal,
76 base::_ndof_on_subgeo,
77 base::_nnod_on_subgeo_internal,
78 base::_nnod_on_subgeo,
79 base::_first_idof_by_dimension_internal,
80 base::_first_idof_by_dimension,
81 base::_first_inod_by_dimension_internal,
82 base::_first_inod_by_dimension);
83}
84template<class T>
85void
87{
88 // initialization is similar to Pk-Lagrange
89 size_type k = degree();
90 size_type variant = hat_K.variant();
91
92 fatal_macro ("sherwin: not yet");
93#ifdef TODO // requires a _raw_basis here
94 // nodes: TODO: volume-internal nodes should be quadrature, because of integral dofs
95 switch (base::_sopt.get_node()) {
97 pointset_lagrange_equispaced (hat_K, k, base::_hat_node[variant]);
98 break;
100 pointset_lagrange_warburton (hat_K, k, base::_hat_node[variant]); break;
101 default: error_macro ("unsupported node set: "<<base::_sopt.get_node_name());
102 }
103 // vdm:
104 details::basis_on_pointset_evaluate (_raw_basis, hat_K, base::_hat_node[variant], base::_vdm[variant]);
105 check_macro (invert(base::_vdm[variant], base::_inv_vdm[variant]),
106 "unisolvence failed for " << base::name() <<"(" << hat_K.name() << ") basis");
107#endif // TODO
108}
109// evaluation of all basis functions at hat_x:
110template<class T>
111void
113 reference_element hat_K,
114 const point_basic<T>& hat_x,
115 Eigen::Matrix<T,Eigen::Dynamic,1>& value) const
116{
117 base::_initialize_data_guard (hat_K);
118
119 Eigen::Matrix<T,Eigen::Dynamic,1> _work0, _work1, _work2; // TODO: allocate one time for all in class ?
120
121 eval_sherwin_basis (hat_x, hat_K, degree(), _alpha, _beta, _work0, _work1, _work2, value);
122}
123// evaluate the gradient:
124template<class T>
125void
127 reference_element hat_K,
128 const point_basic<T>& hat_x,
129 Eigen::Matrix<point_basic<T>,Eigen::Dynamic,1>& value) const
130{
131 base::_initialize_data_guard (hat_K);
132 point_basic<ad3_basic<T> > hat_x_ad = ad3::point (hat_x);
133
134 std::vector<ad3_basic<T> > _work0_ad, _work1_ad, _work2_ad; // TODO: allocate one time for all in class ?
135 std::array<std::vector<ad3_basic<T> >,
136 reference_element::max_variant> _value_ad; // TODO: idem
137
138 std::vector<ad3_basic<T> >& value_ad = _value_ad [hat_K.variant()];
139
140 eval_sherwin_basis (hat_x_ad, hat_K, degree(), _alpha, _beta, _work0_ad, _work1_ad, _work2_ad, value_ad);
141 size_t loc_ndof = value_ad.size();
142 value.resize(loc_ndof);
143 for (size_t loc_idof = 0; loc_idof < loc_ndof; loc_idof++) {
144 value[loc_idof] = value_ad[loc_idof].grad();
145 }
146}
147// dofs for a scalar-valued function
148template<class T>
149void
151 reference_element hat_K,
152 const Eigen::Matrix<T,Eigen::Dynamic,1>& f_xnod,
153 Eigen::Matrix<T,Eigen::Dynamic,1>& dof) const
154{
155#ifdef TODO
156 // TODO: nodal values + integrals inside, instead of Lagrange node
157 base::_initialize_data_guard (hat_K);
158 dof = base::_inv_vdm[hat_K.variant()]*f_xnod;
159#endif // TODO
160}
161// ----------------------------------------------------------------------------
162// instanciation in library
163// ----------------------------------------------------------------------------
164#define _RHEOLEF_instanciation(T) \
165template class basis_fem_Pk_sherwin<T>;
166
168
169}// namespace rheolef
#define _RHEOLEF_instanciation(T, M, A)
Definition asr.cc:223
see the Float page for the full documentation
static void initialize_local_first(size_type k, bool is_continuous, std::array< std::array< size_type, reference_element::max_variant >, 4 > &ndof_on_subgeo_internal, std::array< std::array< size_type, reference_element::max_variant >, 4 > &ndof_on_subgeo, std::array< std::array< size_type, reference_element::max_variant >, 4 > &nnod_on_subgeo_internal, std::array< std::array< size_type, reference_element::max_variant >, 4 > &nnod_on_subgeo, std::array< std::array< size_type, 5 >, reference_element::max_variant > &first_idof_by_dimension_internal, std::array< std::array< size_type, 5 >, reference_element::max_variant > &first_idof_by_dimension, std::array< std::array< size_type, 5 >, reference_element::max_variant > &first_inod_by_dimension_internal, std::array< std::array< size_type, 5 >, reference_element::max_variant > &first_inod_by_dimension)
basis_fem_Pk_sherwin(size_type degree, const basis_option &sopt)
reference_element::size_type size_type
void evaluate(reference_element hat_K, const point_basic< T > &hat_x, Eigen::Matrix< T, Eigen::Dynamic, 1 > &value) const
void grad_evaluate(reference_element hat_K, const point_basic< T > &hat_x, Eigen::Matrix< point_basic< T >, Eigen::Dynamic, 1 > &value) const
void _compute_dofs(reference_element hat_K, const Eigen::Matrix< T, Eigen::Dynamic, 1 > &f_xnod, Eigen::Matrix< T, Eigen::Dynamic, 1 > &dof) const
void _initialize_data(reference_element hat_K) const
see the basis_option page for the full documentation
see the reference_element page for the full documentation
static const variant_type max_variant
variant_type variant() const
#define error_macro(message)
Definition dis_macros.h:49
#define fatal_macro(message)
Definition dis_macros.h:33
#define warning_macro(message)
Definition dis_macros.h:53
Expr1::float_type T
Definition field_expr.h:230
check_macro(expr1.have_homogeneous_space(Xh1), "dual(expr1,expr2); expr1 should have homogeneous space. HINT: use dual(interpolate(Xh, expr1),expr2)")
void basis_on_pointset_evaluate(const Basis &b, const reference_element &hat_K, const Eigen::Matrix< point_basic< T >, Eigen::Dynamic, 1 > &hat_x, Eigen::Matrix< Value, Eigen::Dynamic, Eigen::Dynamic > &vdm)
This file is part of Rheolef.
void pointset_lagrange_warburton(reference_element hat_K, size_t degree, Eigen::Matrix< point_basic< T >, Eigen::Dynamic, 1 > &hat_xnod, bool map_on_reference_element=true)
void invert(tiny_matrix< T > &a, tiny_matrix< T > &inv_a)
Definition tiny_lu.h:127
void pointset_lagrange_equispaced(reference_element hat_K, size_t order_in, Eigen::Matrix< point_basic< T >, Eigen::Dynamic, 1 > &hat_xnod, size_t internal=0)
STL namespace.