Rheolef  7.1
an efficient C++ finite element environment
basis_fem_Pk_sherwin.cc
Go to the documentation of this file.
1 #include "basis_fem_Pk_sherwin.h"
22 #include "basis_fem_Pk_lagrange.h"
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 
31 namespace rheolef {
32 using namespace std;
33 
34 // =========================================================================
35 // basis members
36 // =========================================================================
37 template<class T>
39 {
40 }
41 template<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 {
59 warning_macro ("cstor...");
60  _alpha = _beta = 1; // 1 or 2, see SheKar-2005 p. 53
62 
63  // piola FEM transformatiion:
64  typedef piola_fem_lagrange<T> piola_fem_type;
65  base::_piola_fem.piola_fem<T>::base::operator= (new_macro(piola_fem_type));
66 warning_macro ("cstor done");
67 }
68 template<class T>
69 void
71 {
73  degree(),
74  base::is_continuous(),
75  base::_ndof_on_subgeo,
76  base::_nnod_on_subgeo,
77  base::_first_idof_by_dimension,
78  base::_first_inod_by_dimension);
79 }
80 template<class T>
81 void
83 {
84  // initialization is similar to Pk-Lagrange
85  size_type k = degree();
86  size_type variant = hat_K.variant();
87 
88  fatal_macro ("sherwin: not yet");
89 #ifdef TODO // requires a _raw_basis here
90  // nodes: TODO: volume-internal nodes should be quadrature, because of integral dofs
91  switch (base::_sopt.get_node()) {
93  pointset_lagrange_equispaced (hat_K, k, base::_hat_node[variant]);
94  break;
96  pointset_lagrange_warburton (hat_K, k, base::_hat_node[variant]); break;
97  default: error_macro ("unsupported node set: "<<base::_sopt.get_node_name());
98  }
99  // vdm:
100  details::basis_on_pointset_evaluate (_raw_basis, hat_K, base::_hat_node[variant], base::_vdm[variant]);
101  check_macro (invert(base::_vdm[variant], base::_inv_vdm[variant]),
102  "unisolvence failed for " << base::name() <<"(" << hat_K.name() << ") basis");
103 #endif // TODO
104 }
105 // evaluation of all basis functions at hat_x:
106 template<class T>
107 void
109  reference_element hat_K,
110  const point_basic<T>& hat_x,
111  Eigen::Matrix<T,Eigen::Dynamic,1>& value) const
112 {
113  base::_initialize_data_guard (hat_K);
114 
115  Eigen::Matrix<T,Eigen::Dynamic,1> _work0, _work1, _work2; // TODO: allocate one time for all in class ?
116 
117  eval_sherwin_basis (hat_x, hat_K, degree(), _alpha, _beta, _work0, _work1, _work2, value);
118 }
119 // evaluate the gradient:
120 template<class T>
121 void
123  reference_element hat_K,
124  const point_basic<T>& hat_x,
125  Eigen::Matrix<point_basic<T>,Eigen::Dynamic,1>& value) const
126 {
127  base::_initialize_data_guard (hat_K);
128  point_basic<ad3_basic<T> > hat_x_ad = ad3::point (hat_x);
129 
130  std::vector<ad3_basic<T> > _work0_ad, _work1_ad, _work2_ad; // TODO: allocate one time for all in class ?
131  std::array<std::vector<ad3_basic<T> >,
132  reference_element::max_variant> _value_ad; // TODO: idem
133 
134  std::vector<ad3_basic<T> >& value_ad = _value_ad [hat_K.variant()];
135 
136  eval_sherwin_basis (hat_x_ad, hat_K, degree(), _alpha, _beta, _work0_ad, _work1_ad, _work2_ad, value_ad);
137  size_t loc_ndof = value_ad.size();
138  value.resize(loc_ndof);
139  for (size_t loc_idof = 0; loc_idof < loc_ndof; loc_idof++) {
140  value[loc_idof] = value_ad[loc_idof].grad();
141  }
142 }
143 // dofs for a scalar-valued function
144 template<class T>
145 void
147  reference_element hat_K,
148  const Eigen::Matrix<T,Eigen::Dynamic,1>& f_xnod,
149  Eigen::Matrix<T,Eigen::Dynamic,1>& dof) const
150 {
151 #ifdef TODO
152  // TODO: nodal values + integrals inside, instead of Lagrange node
153  base::_initialize_data_guard (hat_K);
154  dof = base::_inv_vdm[hat_K.variant()]*f_xnod;
155 #endif // TODO
156 }
157 // ----------------------------------------------------------------------------
158 // instanciation in library
159 // ----------------------------------------------------------------------------
160 #define _RHEOLEF_instanciation(T) \
161 template class basis_fem_Pk_sherwin<T>;
162 
164 
165 }// namespace rheolef
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, 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)
basis_fem_Pk_sherwin(size_type degree, const basis_option &sopt)
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
Definition: basis_option.h:93
reference_element::size_type size_type
Definition: basis.h:214
piola_fem< T > _piola_fem
Definition: basis.h:394
basis_option _sopt
Definition: basis.h:393
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:392
see the reference_element page for the full documentation
static const variant_type max_variant
variant_type variant() const
rheolef::std value
#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:261
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_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