Rheolef  7.2
an efficient C++ finite element environment
 
Loading...
Searching...
No Matches
basis_fem_Pk_lagrange.cc
Go to the documentation of this file.
1
22#include "rheolef/rheostream.h"
23
24#include "piola_fem_lagrange.h"
25#include "equispaced.icc"
26#include "warburton.icc"
27#include "fekete.icc"
30#include "eigen_util.h"
31
32namespace rheolef {
33
34using namespace std;
35
36// =========================================================================
37// utilities
38// =========================================================================
39template <class T>
40void
42 size_type k,
43 bool is_continuous,
44 std::array<std::array<size_type,reference_element::max_variant>,4>& ndof_on_subgeo_internal,
45 std::array<std::array<size_type,reference_element::max_variant>,4>& ndof_on_subgeo,
46 std::array<std::array<size_type,reference_element::max_variant>,4>& nnod_on_subgeo_internal,
47 std::array<std::array<size_type,reference_element::max_variant>,4>& nnod_on_subgeo,
48 std::array<std::array<size_type,5>,reference_element::max_variant>& first_idof_by_dimension_internal,
49 std::array<std::array<size_type,5>,reference_element::max_variant>& first_idof_by_dimension,
50 std::array<std::array<size_type,5>,reference_element::max_variant>& first_inod_by_dimension_internal,
51 std::array<std::array<size_type,5>,reference_element::max_variant>& first_inod_by_dimension)
52{
53 // 1) ndof_on_subgeo
54 if (k == size_t(-1)) {
55 // P_{-1} == empty polynomial space
56 for (size_type map_d = 0; map_d < 4; ++map_d) {
57 ndof_on_subgeo_internal [map_d].fill (0);
58 }
59 } else if (k == 0) {
60 // P0 initialization
61 for (size_type map_d = 0; map_d < 4; ++map_d) {
62 ndof_on_subgeo_internal [map_d].fill (0);
65 ++variant) {
66 ndof_on_subgeo_internal [map_d] [variant] = 1;
67 }
68 }
69 } else { // k > 0
70 for (size_type map_d = 0; map_d < 4; ++map_d) {
71 reference_element::init_local_nnode_by_variant (k, ndof_on_subgeo_internal [map_d]);
72 // clean upper-dimensional subgeos, since init_local_nnode_by_variant do not consider dimension
74 variant < reference_element::max_variant; ++variant) {
75 ndof_on_subgeo_internal [map_d] [variant] = 0;
76 }
77 }
78 }
79 // when discontinuous, fix, but conserve the subgeo structure in _internal:
80 base::_helper_make_discontinuous_ndof_on_subgeo (is_continuous, ndof_on_subgeo_internal, ndof_on_subgeo);
81
82 // 2) deduce automatically first_idof_by_dimension:
83 base::_helper_initialize_first_ixxx_by_dimension_from_nxxx_on_subgeo (ndof_on_subgeo_internal, first_idof_by_dimension_internal);
84 base::_helper_initialize_first_ixxx_by_dimension_from_nxxx_on_subgeo (ndof_on_subgeo, first_idof_by_dimension);
85 // 3) Lagrange nodes follow dofs:
86 nnod_on_subgeo_internal = ndof_on_subgeo_internal;
87 nnod_on_subgeo = ndof_on_subgeo;
88 first_inod_by_dimension_internal = first_idof_by_dimension_internal;
89 first_inod_by_dimension = first_idof_by_dimension;
90}
91// =========================================================================
92// basis members
93// =========================================================================
94template<class T>
98template<class T>
100 : basis_rep<T> (sopt),
101 _raw_basis(),
102 _hat_node(),
103 _vdm(),
104 _inv_vdm()
105{
106 base::_name = base::standard_naming (family_name(), degree, base::_sopt);
107 string R = "?";
108 switch (base::_sopt.get_raw_polynomial()) {
109 case basis_option::monomial: R = "M"; break;
110 case basis_option::bernstein: R = "B"; break;
111 case basis_option::dubiner: R = "D"; break;
112 default: error_macro ("unsupported polynomial: "<<sopt.get_raw_polynomial_name());
113 }
114 _raw_basis = basis_raw_basic<T> (R+std::to_string(degree));
116
117 // piola FEM transformation:
119 base::_piola_fem.piola_fem<T>::base::operator= (new_macro(piola_fem_type));
120}
121template<class T>
122bool
124{
125 return true;
126}
127template<class T>
128const Eigen::Matrix<point_basic<T>,Eigen::Dynamic,1>&
130{
131 base::_initialize_data_guard (hat_K);
132 return _hat_node [hat_K.variant()];
133}
134template<class T>
135const Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>&
137{
138 base::_initialize_data_guard (hat_K);
139 return _vdm [hat_K.variant()];
140}
141template<class T>
142const Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>&
144{
145 base::_initialize_data_guard (hat_K);
146 return _inv_vdm [hat_K.variant()];
147}
148template<class T>
149void
151{
153 degree(),
154 base::is_continuous(),
155 base::_ndof_on_subgeo_internal,
156 base::_ndof_on_subgeo,
157 base::_nnod_on_subgeo_internal,
158 base::_nnod_on_subgeo,
159 base::_first_idof_by_dimension_internal,
160 base::_first_idof_by_dimension,
161 base::_first_inod_by_dimension_internal,
162 base::_first_inod_by_dimension);
163}
164template<class T>
165void
167{
168 size_type k = degree();
169 size_type variant = hat_K.variant();
170
171 // nodes:
172 switch (base::_sopt.get_node()) {
174 pointset_lagrange_equispaced (hat_K, k, _hat_node[variant]); break;
176 pointset_lagrange_warburton (hat_K, k, _hat_node[variant]); break;
178 pointset_lagrange_fekete (hat_K, k, _hat_node[variant]); break;
179 default: error_macro ("unsupported node set: "<<base::_sopt.get_node_name());
180 }
181 // vdm:
182 details::basis_on_pointset_evaluate (_raw_basis, hat_K, _hat_node[variant], _vdm[variant]);
183 check_macro (invert(_vdm[variant], _inv_vdm[variant]),
184 "unisolvence failed for " << base::name() <<"(" << hat_K.name() << ") basis");
185}
186// evaluation of all basis functions at hat_x:
187template<class T>
188void
190 reference_element hat_K,
191 const point_basic<T>& hat_x,
192 Eigen::Matrix<T,Eigen::Dynamic,1>& value) const
193{
194 base::_initialize_data_guard (hat_K);
195 Eigen::Matrix<T,Eigen::Dynamic,1> raw_value;
196 _raw_basis.evaluate (hat_K, hat_x, raw_value);
197 value = _inv_vdm[hat_K.variant()].transpose()*raw_value;
198}
199// evaluate the gradient:
200template<class T>
201void
203 reference_element hat_K,
204 const point_basic<T>& hat_x,
205 Eigen::Matrix<point_basic<T>,Eigen::Dynamic,1>& value) const
206{
207 base::_initialize_data_guard (hat_K);
208 Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>& inv_vdm = _inv_vdm [hat_K.variant()];
209 Eigen::Matrix<point_basic<T>,Eigen::Dynamic,1> raw_v_value;
210 _raw_basis.grad_evaluate (hat_K, hat_x, raw_v_value);
211 size_type loc_ndof = raw_v_value.size();
212 value.resize (loc_ndof);
213 // TODO: point-valued blas2 : value := trans(inv_vdm)*raw_value
214 for (size_type loc_idof = 0; loc_idof < loc_ndof; ++loc_idof) {
215 value [loc_idof] = point_basic<T>(0,0,0);
216 for (size_type loc_jdof = 0; loc_jdof < loc_ndof; ++loc_jdof) {
217 value [loc_idof] += inv_vdm (loc_jdof,loc_idof)*raw_v_value[loc_jdof];
218 }
219 }
220}
221// extract local dof-indexes on a side
222template<class T>
225 reference_element hat_K,
226 const side_information_type& sid) const
227{
228 return base::ndof (sid.hat);
229}
230template<class T>
231void
233 reference_element hat_K,
234 const side_information_type& sid,
235 Eigen::Matrix<size_type,Eigen::Dynamic,1>& loc_idof) const
236{
237 details::Pk_get_local_idof_on_side (hat_K, sid, degree(), loc_idof);
238}
239// dofs for a scalar-valued function
240template<class T>
241void
243 reference_element hat_K,
244 const Eigen::Matrix<T,Eigen::Dynamic,1>& f_xnod,
245 Eigen::Matrix<T,Eigen::Dynamic,1>& dof) const
246{
247 dof = f_xnod; // TODO: could avoid a physical copy: check it in compute_dof(f) with is_nodal()
248}
249// ----------------------------------------------------------------------------
250// instanciation in library
251// ----------------------------------------------------------------------------
252#define _RHEOLEF_instanciation(T) \
253template class basis_fem_Pk_lagrange<T>;
254
256
257}// namespace rheolef
#define _RHEOLEF_instanciation(T, M, A)
Definition asr.cc:223
see the Float page for the full documentation
const Eigen::Matrix< T, Eigen::Dynamic, Eigen::Dynamic > & vdm(reference_element hat_K) const
reference_element::size_type size_type
void local_idof_on_side(reference_element hat_K, const side_information_type &sid, Eigen::Matrix< size_type, Eigen::Dynamic, 1 > &loc_idof) const
void evaluate(reference_element hat_K, const point_basic< T > &hat_x, Eigen::Matrix< T, Eigen::Dynamic, 1 > &value) const
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)
size_type local_ndof_on_side(reference_element hat_K, const side_information_type &sid) 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
basis_fem_Pk_lagrange(size_type degree, const basis_option &sopt)
const Eigen::Matrix< point_basic< T >, Eigen::Dynamic, 1 > & hat_node(reference_element hat_K) const
virtual std::string family_name() const
const Eigen::Matrix< T, Eigen::Dynamic, Eigen::Dynamic > & inv_vdm(reference_element hat_K) 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
static void init_local_nnode_by_variant(size_type order, std::array< size_type, reference_element::max_variant > &loc_nnod_by_variant)
static variant_type last_variant_by_dimension(size_type dim)
static variant_type first_variant_by_dimension(size_type dim)
variant_type variant() const
#define error_macro(message)
Definition dis_macros.h:49
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 pointset_lagrange_fekete(reference_element hat_K, size_t degree, Eigen::Matrix< point_basic< T >, Eigen::Dynamic, 1 > &hat_xnod, bool map_on_reference_element=true)
Definition fekete.icc:42
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.