Rheolef  7.1
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,
45  std::array<std::array<size_type,reference_element::max_variant>,4>& nnod_on_subgeo,
46  std::array<std::array<size_type,5>,reference_element::max_variant>& first_idof_by_dimension,
47  std::array<std::array<size_type,5>,reference_element::max_variant>& first_inod_by_dimension)
48 {
49  // 1) ndof_on_subgeo
50  if (k == size_t(-1)) {
51  // P_{-1} == empty polynomial space
52  for (size_type map_d = 0; map_d < 4; ++map_d) {
53  ndof_on_subgeo [map_d].fill (0);
54  }
55  } else if (k == 0) {
56  // P0 initialization
57  for (size_type map_d = 0; map_d < 4; ++map_d) {
58  ndof_on_subgeo [map_d].fill (0);
61  ++variant) {
62  ndof_on_subgeo [map_d] [variant] = 1;
63  }
64  }
65  } else { // k > 0
66  for (size_type map_d = 0; map_d < 4; ++map_d) {
67  reference_element::init_local_nnode_by_variant (k, ndof_on_subgeo [map_d]);
68  // clean upper-dimensional subgeos, since init_local_nnode_by_variant do not consider dimension
71  ndof_on_subgeo [map_d] [variant] = 0;
72  }
73  }
74  // when discontinuous, fix:
75  base::_helper_discontinuous_ndof_on_subgeo_inplace_change (is_continuous, ndof_on_subgeo);
76  }
77  // 2) deduce automatically first_idof_by_dimension:
78  base::_helper_initialize_first_ixxx_by_dimension_from_nxxx_on_subgeo (ndof_on_subgeo, first_idof_by_dimension);
79  // 3) Lagrange nodes follow dofs:
80  nnod_on_subgeo = ndof_on_subgeo;
81  first_inod_by_dimension = first_idof_by_dimension;
82 }
83 // =========================================================================
84 // basis members
85 // =========================================================================
86 template<class T>
88 {
89 }
90 template<class T>
92  : basis_rep<T> (sopt),
93  _raw_basis(),
94  _hat_node(),
95  _vdm(),
96  _inv_vdm()
97 {
99  string R = "?";
100  switch (base::_sopt.get_raw_polynomial()) {
101  case basis_option::monomial: R = "M"; break;
102  case basis_option::bernstein: R = "B"; break;
103  case basis_option::dubiner: R = "D"; break;
104  default: error_macro ("unsupported polynomial: "<<sopt.get_raw_polynomial_name());
105  }
108 
109  // piola FEM transformation:
110  typedef piola_fem_lagrange<T> piola_fem_type;
111  base::_piola_fem.piola_fem<T>::base::operator= (new_macro(piola_fem_type));
112 }
113 template<class T>
114 bool
116 {
117  return true;
118 }
119 template<class T>
120 const Eigen::Matrix<point_basic<T>,Eigen::Dynamic,1>&
122 {
123  base::_initialize_data_guard (hat_K);
124  return _hat_node [hat_K.variant()];
125 }
126 template<class T>
127 const Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>&
129 {
130  base::_initialize_data_guard (hat_K);
131  return _vdm [hat_K.variant()];
132 }
133 template<class T>
134 const Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>&
136 {
137  base::_initialize_data_guard (hat_K);
138  return _inv_vdm [hat_K.variant()];
139 }
140 template<class T>
141 void
143 {
145  degree(),
146  base::is_continuous(),
147  base::_ndof_on_subgeo,
148  base::_nnod_on_subgeo,
149  base::_first_idof_by_dimension,
150  base::_first_inod_by_dimension);
151 }
152 template<class T>
153 void
155 {
156  size_type k = degree();
157  size_type variant = hat_K.variant();
158 
159  // nodes:
160  switch (base::_sopt.get_node()) {
162  pointset_lagrange_equispaced (hat_K, k, _hat_node[variant]); break;
164  pointset_lagrange_warburton (hat_K, k, _hat_node[variant]); break;
166  pointset_lagrange_fekete (hat_K, k, _hat_node[variant]); break;
167  default: error_macro ("unsupported node set: "<<base::_sopt.get_node_name());
168  }
169  // vdm:
170  details::basis_on_pointset_evaluate (_raw_basis, hat_K, _hat_node[variant], _vdm[variant]);
171  check_macro (invert(_vdm[variant], _inv_vdm[variant]),
172  "unisolvence failed for " << base::name() <<"(" << hat_K.name() << ") basis");
173 }
174 // evaluation of all basis functions at hat_x:
175 template<class T>
176 void
178  reference_element hat_K,
179  const point_basic<T>& hat_x,
180  Eigen::Matrix<T,Eigen::Dynamic,1>& value) const
181 {
182  base::_initialize_data_guard (hat_K);
183  Eigen::Matrix<T,Eigen::Dynamic,1> raw_value;
184  _raw_basis.evaluate (hat_K, hat_x, raw_value);
185  value = _inv_vdm[hat_K.variant()].transpose()*raw_value;
186 }
187 // evaluate the gradient:
188 template<class T>
189 void
191  reference_element hat_K,
192  const point_basic<T>& hat_x,
193  Eigen::Matrix<point_basic<T>,Eigen::Dynamic,1>& value) const
194 {
195  base::_initialize_data_guard (hat_K);
196  Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>& inv_vdm = _inv_vdm [hat_K.variant()];
197  Eigen::Matrix<point_basic<T>,Eigen::Dynamic,1> raw_v_value;
198  _raw_basis.grad_evaluate (hat_K, hat_x, raw_v_value);
199  size_type loc_ndof = raw_v_value.size();
200  value.resize (loc_ndof);
201  // TODO: point-valued blas2 : value := trans(inv_vdm)*raw_value
202  for (size_type loc_idof = 0; loc_idof < loc_ndof; ++loc_idof) {
203  value [loc_idof] = point_basic<T>(0,0,0);
204  for (size_type loc_jdof = 0; loc_jdof < loc_ndof; ++loc_jdof) {
205  value [loc_idof] += inv_vdm (loc_jdof,loc_idof)*raw_v_value[loc_jdof];
206  }
207  }
208 }
209 // extract local dof-indexes on a side
210 template<class T>
213  reference_element hat_K,
214  const side_information_type& sid) const
215 {
216  return base::ndof (sid.hat);
217 }
218 template<class T>
219 void
221  reference_element hat_K,
222  const side_information_type& sid,
223  Eigen::Matrix<size_type,Eigen::Dynamic,1>& loc_idof) const
224 {
225  details::Pk_get_local_idof_on_side (hat_K, sid, degree(), loc_idof);
226 }
227 // dofs for a scalar-valued function
228 template<class T>
229 void
231  reference_element hat_K,
232  const Eigen::Matrix<T,Eigen::Dynamic,1>& f_xnod,
233  Eigen::Matrix<T,Eigen::Dynamic,1>& dof) const
234 {
235  dof = f_xnod; // TODO: could avoid a physical copy: check it in compute_dof(f) with is_nodal()
236 }
237 // ----------------------------------------------------------------------------
238 // instanciation in library
239 // ----------------------------------------------------------------------------
240 #define _RHEOLEF_instanciation(T) \
241 template class basis_fem_Pk_lagrange<T>;
242 
244 
245 }// namespace rheolef
rheolef::space_numbering::ndof
size_type ndof(const basis_basic< T > &b, const geo_size &gs, size_type map_dim)
Definition: space_numbering.cc:28
rheolef::reference_element::last_variant_by_dimension
static variant_type last_variant_by_dimension(size_type dim)
Definition: reference_element.h:150
warburton.icc
rheolef::basis_fem_Pk_lagrange::local_ndof_on_side
size_type local_ndof_on_side(reference_element hat_K, const side_information_type &sid) const
Definition: basis_fem_Pk_lagrange.cc:212
eigen_util.h
rheolef::point_basic
Definition: point.h:87
rheolef::reference_element::init_local_nnode_by_variant
static void init_local_nnode_by_variant(size_type order, std::array< size_type, reference_element::max_variant > &loc_nnod_by_variant)
Definition: reference_element.cc:206
check_macro
check_macro(expr1.have_homogeneous_space(Xh1), "dual(expr1,expr2); expr1 should have homogeneous space. HINT: use dual(interpolate(Xh, expr1),expr2)")
rheolef::basis_fem_Pk_lagrange::hat_node
const Eigen::Matrix< point_basic< T >, Eigen::Dynamic, 1 > & hat_node(reference_element hat_K) const
Definition: basis_fem_Pk_lagrange.cc:121
rheolef::piola_fem_lagrange
Definition: piola_fem_lagrange.h:65
rheolef::basis_fem_Pk_lagrange::evaluate
void evaluate(reference_element hat_K, const point_basic< T > &hat_x, Eigen::Matrix< T, Eigen::Dynamic, 1 > &value) const
Definition: basis_fem_Pk_lagrange.cc:177
rheolef::pointset_lagrange_equispaced
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
rheolef::details::basis_on_pointset_evaluate
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)
Definition: basis_on_pointset_evaluate.icc:31
rheolef::basis_fem_Pk_lagrange::_raw_basis
basis_raw_basic< T > _raw_basis
Definition: basis_fem_Pk_lagrange.h:132
basis_fem_Pk_lagrange.h
rheolef::_RHEOLEF_instanciation
_RHEOLEF_instanciation(Float, sequential, std::allocator< Float >) _RHEOLEF_instanciation(Float
rheolef::basis_raw_basic
Definition: basis_raw.h:107
rheolef::basis_fem_Pk_lagrange::grad_evaluate
void grad_evaluate(reference_element hat_K, const point_basic< T > &hat_x, Eigen::Matrix< point_basic< T >, Eigen::Dynamic, 1 > &value) const
Definition: basis_fem_Pk_lagrange.cc:190
rheolef::value
rheolef::std value
rheolef::basis_fem_Pk_lagrange::inv_vdm
const Eigen::Matrix< T, Eigen::Dynamic, Eigen::Dynamic > & inv_vdm(reference_element hat_K) const
Definition: basis_fem_Pk_lagrange.cc:135
rheolef::basis_option::warburton
@ warburton
Definition: basis_option.h:103
rheolef::basis_option::bernstein
@ bernstein
Definition: basis_option.h:110
rheolef::basis_option
see the basis_option page for the full documentation
Definition: basis_option.h:93
rheolef::basis_rep::size_type
reference_element::size_type size_type
Definition: basis.h:214
rheolef::reference_element::first_variant_by_dimension
static variant_type first_variant_by_dimension(size_type dim)
Definition: reference_element.h:148
fekete.icc
mkgeo_ball.variant
int variant
Definition: mkgeo_ball.sh:149
rheolef::side_information_type
Definition: reference_element_face_transformation.h:37
rheolef::reference_element
see the reference_element page for the full documentation
Definition: reference_element.h:66
rheolef::basis_fem_Pk_lagrange::local_idof_on_side
void local_idof_on_side(reference_element hat_K, const side_information_type &sid, Eigen::Matrix< size_type, Eigen::Dynamic, 1 > &loc_idof) const
Definition: basis_fem_Pk_lagrange.cc:220
rheolef::basis_rep::_piola_fem
piola_fem< T > _piola_fem
Definition: basis.h:394
rheolef::basis_fem_Pk_lagrange::degree
size_type degree() const
Definition: basis_fem_Pk_lagrange.h:69
rheolef::reference_element::name
char name() const
Definition: reference_element.h:100
Pk_get_local_idof_on_side.icc
rheolef::basis_fem_Pk_lagrange::vdm
const Eigen::Matrix< T, Eigen::Dynamic, Eigen::Dynamic > & vdm(reference_element hat_K) const
Definition: basis_fem_Pk_lagrange.cc:128
rheolef::reference_element::variant
variant_type variant() const
Definition: reference_element.h:99
rheolef::basis_option::monomial
@ monomial
Definition: basis_option.h:109
rheolef::basis_fem_Pk_lagrange::_compute_dofs
void _compute_dofs(reference_element hat_K, const Eigen::Matrix< T, Eigen::Dynamic, 1 > &f_xnod, Eigen::Matrix< T, Eigen::Dynamic, 1 > &dof) const
Definition: basis_fem_Pk_lagrange.cc:230
rheolef::pointset_lagrange_warburton
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::basis_rep::_sopt
basis_option _sopt
Definition: basis.h:393
rheolef::basis_rep::_name
std::string _name
Definition: basis.h:392
rheolef
This file is part of Rheolef.
Definition: compiler_eigen.h:37
rheolef::basis_fem_Pk_lagrange::initialize_local_first
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, 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, std::array< std::array< size_type, 5 >, reference_element::max_variant > &first_inod_by_dimension)
Definition: basis_fem_Pk_lagrange.cc:41
error_macro
#define error_macro(message)
Definition: dis_macros.h:49
Float
see the Float page for the full documentation
basis_on_pointset_evaluate.icc
rheolef::basis_fem_Pk_lagrange::basis_fem_Pk_lagrange
basis_fem_Pk_lagrange(size_type degree, const basis_option &sopt)
Definition: basis_fem_Pk_lagrange.cc:91
rheolef::reference_element::max_variant
static const variant_type max_variant
Definition: reference_element.h:82
rheolef::basis_option::equispaced
@ equispaced
Definition: basis_option.h:102
rheolef::basis_fem_Pk_lagrange::is_nodal
bool is_nodal() const
Definition: basis_fem_Pk_lagrange.cc:115
rheolef::basis_fem_Pk_lagrange::_initialize_cstor_sizes
void _initialize_cstor_sizes() const
Definition: basis_fem_Pk_lagrange.cc:142
rheolef::point_basic< T >
point_basic< T >
Definition: piola_fem.h:135
equispaced.icc
rheolef::basis_rep::standard_naming
static std::string standard_naming(std::string family_name, size_t degree, const basis_option &sopt)
Definition: basis_rep.cc:44
rheolef::invert
void invert(tiny_matrix< T > &a, tiny_matrix< T > &inv_a)
Definition: tiny_lu.h:127
piola_fem_lagrange.h
rheolef::basis_fem_Pk_lagrange::size_type
reference_element::size_type size_type
Definition: basis_fem_Pk_lagrange.h:58
rheolef::basis_fem_Pk_lagrange::~basis_fem_Pk_lagrange
~basis_fem_Pk_lagrange()
Definition: basis_fem_Pk_lagrange.cc:87
rheolef::basis_option::dubiner
@ dubiner
Definition: basis_option.h:111
rheolef::basis_option::fekete
@ fekete
Definition: basis_option.h:104
rheolef::itos
std::string itos(std::string::size_type i)
itos: see the rheostream page for the full documentation
rheolef::side_information_type::hat
reference_element hat
Definition: reference_element_face_transformation.h:43
rheolef::pointset_lagrange_fekete
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
rheolef::basis_fem_Pk_lagrange::_initialize_data
void _initialize_data(reference_element hat_K) const
Definition: basis_fem_Pk_lagrange.cc:154
mkgeo_contraction.name
string name
Definition: mkgeo_contraction.sh:133
rheolef::basis_option::get_raw_polynomial_name
std::string get_raw_polynomial_name() const
Definition: basis_option.cc:104
rheolef::basis_fem_Pk_lagrange::family_name
virtual std::string family_name() const
Definition: basis_fem_Pk_lagrange.h:68
rheolef::basis_rep
Definition: basis.h:209
T
Expr1::float_type T
Definition: field_expr.h:261