Rheolef  7.2
an efficient C++ finite element environment
basis_fem_Pk_lagrange.cc
Go to the documentation of this file.
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 
32 namespace rheolef {
33 
34 using namespace std;
35 
36 // =========================================================================
37 // utilities
38 // =========================================================================
39 template <class T>
40 void
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
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 // =========================================================================
94 template<class T>
96 {
97 }
98 template<class T>
100  : basis_rep<T> (sopt),
101  _raw_basis(),
102  _hat_node(),
103  _vdm(),
104  _inv_vdm()
105 {
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:
118  typedef piola_fem_lagrange<T> piola_fem_type;
119  base::_piola_fem.piola_fem<T>::base::operator= (new_macro(piola_fem_type));
120 }
121 template<class T>
122 bool
124 {
125  return true;
126 }
127 template<class T>
128 const Eigen::Matrix<point_basic<T>,Eigen::Dynamic,1>&
130 {
131  base::_initialize_data_guard (hat_K);
132  return _hat_node [hat_K.variant()];
133 }
134 template<class T>
135 const Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>&
137 {
138  base::_initialize_data_guard (hat_K);
139  return _vdm [hat_K.variant()];
140 }
141 template<class T>
142 const Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>&
144 {
145  base::_initialize_data_guard (hat_K);
146  return _inv_vdm [hat_K.variant()];
147 }
148 template<class T>
149 void
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 }
164 template<class T>
165 void
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:
187 template<class T>
188 void
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:
200 template<class T>
201 void
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
222 template<class T>
225  reference_element hat_K,
226  const side_information_type& sid) const
227 {
228  return base::ndof (sid.hat);
229 }
230 template<class T>
231 void
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
240 template<class T>
241 void
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) \
253 template class basis_fem_Pk_lagrange<T>;
254 
256 
257 }// namespace rheolef
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
Definition: basis_option.h:93
std::string get_raw_polynomial_name() const
reference_element::size_type size_type
Definition: basis.h:214
piola_fem< T > _piola_fem
Definition: basis.h:424
basis_option _sopt
Definition: basis.h:423
static std::string standard_naming(std::string family_name, size_t degree, const basis_option &sopt)
Definition: basis_rep.cc:44
std::string _name
Definition: basis.h:422
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
point_basic< T >
Definition: piola_fem.h:135
rheolef::std value
#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)
size_type ndof(const basis_basic< T > &b, const geo_size &gs, size_type map_dim)
This file is part of Rheolef.
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 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)
Definition: equispaced.icc:44
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)
Definition: warburton.icc:574
_RHEOLEF_instanciation(Float, sequential, std::allocator< Float >) _RHEOLEF_instanciation(Float
void invert(tiny_matrix< T > &a, tiny_matrix< T > &inv_a)
Definition: tiny_lu.h:127