Rheolef  7.2
an efficient C++ finite element environment
 
Loading...
Searching...
No Matches
basis_fem_tensor.cc
Go to the documentation of this file.
1
21#include "basis_fem_tensor.h"
22#include "rheolef/rheostream.h"
23
24namespace rheolef {
25
26using namespace std;
27
28// =========================================================================
29// basis members
30// =========================================================================
31template<class T>
35template<class T>
37 : basis_rep<T>(sopt),
38 _n_comp(0),
39 _scalar_basis(scalar_basis),
40 _scalar_value(),
41 _vector_value()
42{
43 base::_sopt.set_valued_tag (space_constant::tensor);
44 base::_name = base::standard_naming (family_name(), family_index(), base::_sopt);
45 base::_piola_fem = _scalar_basis.get_piola_fem();
46 check_macro (base::option().dimension() != std::numeric_limits<basis_option::size_type>::max(),
47 "tensor(basis): basis.option.map_dimension should be initialized for component number");
48 size_type d = base::option().dimension();
49 _n_comp = d*(d+1)/2;
50 if (d == 2 &&
51 (base::option().coordinate_system() == space_constant::axisymmetric_rz ||
52 base::option().coordinate_system() == space_constant::axisymmetric_zr)) {
53 _n_comp++;
54 }
56}
57#ifdef TO_CLEAN
58template<class T>
61{
62 const size_type unset = std::numeric_limits<basis_option::size_type>::max();
63 size_type d = (base::option().dimension() == unset) ? map_d : base::option().dimension();
64 size_type n_comp = d*(d+1)/2;
65 if (d == 2 &&
66 (base::option().coordinate_system() == space_constant::axisymmetric_rz ||
67 base::option().coordinate_system() == space_constant::axisymmetric_zr)) {
68 n_comp++;
69 }
70 return n_comp;
71}
72#endif // TO_CLEAN
73template<class T>
74void
76{
77 for (size_type map_d = 0; map_d < 4; ++map_d) {
78 for (size_type subgeo_variant = 0; subgeo_variant < reference_element::max_variant; ++subgeo_variant) {
79 base::_ndof_on_subgeo [map_d][subgeo_variant] = _n_comp*_scalar_basis.ndof_on_subgeo (map_d, subgeo_variant);
80 base::_nnod_on_subgeo [map_d][subgeo_variant] = _scalar_basis.nnod_on_subgeo (map_d, subgeo_variant);
81 }
82 }
83 for (size_type variant = 0; variant < reference_element::max_variant; ++variant) {
84 reference_element hat_K (variant);
85 for (size_type subgeo_d = 0; subgeo_d < 5; ++subgeo_d) {
86 base::_first_idof_by_dimension [variant][subgeo_d] = _n_comp*_scalar_basis.first_idof_by_dimension (hat_K, subgeo_d);
87 base::_first_inod_by_dimension [variant][subgeo_d] = _scalar_basis.first_inod_by_dimension (hat_K, subgeo_d);
88 }
89 }
90}
91template<class T>
92void
96template<class T>
97const Eigen::Matrix<point_basic<T>,Eigen::Dynamic,1>&
99{
100 return _scalar_basis.hat_node (hat_K);
101}
102template<class T>
103void
105 reference_element hat_K,
106 const point_basic<T>& hat_x,
107 Eigen::Matrix<tensor_basic<T>,Eigen::Dynamic,1>& value) const
108{
109 base::_initialize_data_guard (hat_K);
110 space_constant::coordinate_type sys_coord = base::option().coordinate_system();
111 _scalar_basis.evaluate (hat_K, hat_x, _scalar_value);
112 size_type loc_comp_ndof = _scalar_value.size();
113 size_type loc_ndof = _n_comp*loc_comp_ndof;
114 value.resize (loc_ndof);
115 value.fill (tensor_basic<T>()); // do not remove !
116 for (size_type loc_comp_idof = 0; loc_comp_idof < loc_comp_ndof; ++loc_comp_idof) {
117 for (size_type ij_comp = 0; ij_comp < _n_comp; ++ij_comp) {
118 size_type loc_idof = _n_comp*loc_comp_idof + ij_comp;
119 std::pair<size_type,size_type> ij = space_constant::tensor_subscript (space_constant::tensor, sys_coord, ij_comp);
120 size_type i_comp = ij.first;
121 size_type j_comp = ij.second;
122 value[loc_idof] (i_comp,j_comp) = _scalar_value[loc_comp_idof];
123 if (ij.first == ij.second) continue;
124 value[loc_idof] (j_comp,i_comp) = _scalar_value[loc_comp_idof];
125 }
126 }
127}
128template<class T>
129void
131 reference_element hat_K,
132 const point_basic<T>& hat_x,
133 Eigen::Matrix<tensor3_basic<T>,Eigen::Dynamic,1>& value) const
134{
135 base::_initialize_data_guard (hat_K);
136 space_constant::coordinate_type sys_coord = base::option().coordinate_system();
137 _scalar_basis.grad_evaluate (hat_K, hat_x, _vector_value);
138 size_type loc_comp_ndof = _vector_value.size();
139 size_type loc_ndof = _n_comp*loc_comp_ndof;
140 value.resize (loc_ndof);
141 value.fill (tensor3_basic<T>()); // do not remove !
142 for (size_type loc_comp_idof = 0; loc_comp_idof < loc_comp_ndof; ++loc_comp_idof) {
143 for (size_type ij_comp = 0; ij_comp < _n_comp; ++ij_comp) {
144 size_type loc_idof = _n_comp*loc_comp_idof + ij_comp;
145 std::pair<size_type,size_type> ij = space_constant::tensor_subscript (space_constant::tensor, sys_coord, ij_comp);
146 size_type i_comp = ij.first;
147 size_type j_comp = ij.second;
148 for (size_type k_comp = 0; k_comp < _n_comp; ++k_comp) {
149 value[loc_idof](i_comp,j_comp,k_comp) = _vector_value[loc_comp_idof][k_comp];
150 if (i_comp == j_comp) continue;
151 value[loc_idof](j_comp,i_comp,k_comp) = _vector_value[loc_comp_idof][k_comp];
152 }
153 }
154 }
155}
156// dofs for a scalar-valued function
157template<class T>
158void
160 reference_element hat_K,
161 const Eigen::Matrix<tensor_basic<T>,Eigen::Dynamic,1>& f_xnod,
162 Eigen::Matrix<T,Eigen::Dynamic,1>& dof) const
163{
164 base::_initialize_data_guard (hat_K);
165 space_constant::coordinate_type sys_coord = base::option().coordinate_system();
166 size_type loc_comp_ndof = _scalar_basis.ndof (hat_K);
167 size_type loc_comp_nnod = _scalar_basis.nnod (hat_K);
168 size_type loc_ndof = _n_comp*loc_comp_ndof;
169 Eigen::Matrix<T,Eigen::Dynamic,1> f_comp_xnod (loc_comp_nnod); // TODO: class working array
170 Eigen::Matrix<T,Eigen::Dynamic,1> comp_dof (loc_comp_ndof); // TODO: class working array
171 dof.resize (loc_ndof);
172 for (size_type ij_comp = 0; ij_comp < _n_comp; ++ij_comp) {
173 std::pair<size_type,size_type> ij = space_constant::tensor_subscript (space_constant::tensor, sys_coord, ij_comp);
174 size_type i_comp = ij.first;
175 size_type j_comp = ij.second;
176 for (size_type loc_comp_inod = 0; loc_comp_inod < loc_comp_nnod; ++loc_comp_inod) {
177 // interpolate the symmetric value of the tensor-valued function
178 f_comp_xnod [loc_comp_inod]
179 = 0.5*( f_xnod [loc_comp_inod] (i_comp,j_comp)
180 + f_xnod [loc_comp_inod] (j_comp,i_comp));
181 }
182 _scalar_basis.compute_dofs (hat_K, f_comp_xnod, comp_dof);
183 for (size_type loc_comp_idof = 0; loc_comp_idof < loc_comp_ndof; ++loc_comp_idof) {
184 size_type loc_idof = _n_comp*loc_comp_idof + ij_comp;
185 dof [loc_idof] = comp_dof [loc_comp_idof];
186 }
187 }
188}
189// ----------------------------------------------------------------------------
190// instanciation in library
191// ----------------------------------------------------------------------------
192#define _RHEOLEF_instanciation(T) \
193template class basis_fem_tensor<T>;
194
196
197}// namespace rheolef
#define _RHEOLEF_instanciation(T, M, A)
Definition asr.cc:223
field::size_type size_type
Definition branch.cc:430
see the Float page for the full documentation
void _compute_dofs(reference_element hat_K, const Eigen::Matrix< tensor_basic< T >, Eigen::Dynamic, 1 > &f_xnod, Eigen::Matrix< T, Eigen::Dynamic, 1 > &dof) const
basis_fem_tensor(const basis_basic< T > &scalar_basis, const basis_option &sopt)
void grad_evaluate(reference_element hat_K, const point_basic< T > &hat_x, Eigen::Matrix< tensor3_basic< T >, Eigen::Dynamic, 1 > &value) const
void evaluate(reference_element hat_K, const point_basic< T > &hat_x, Eigen::Matrix< tensor_basic< T >, Eigen::Dynamic, 1 > &value) const
size_type family_index() const
const Eigen::Matrix< point_basic< T >, Eigen::Dynamic, 1 > & hat_node(reference_element hat_K) const
std::string family_name() 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
const size_t dimension
Definition edge.icc:64
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)")
std::pair< size_type, size_type > tensor_subscript(valued_type valued_tag, coordinate_type sys_coord, size_type i_comp)
This file is part of Rheolef.
STL namespace.