Rheolef  7.2
an efficient C++ finite element environment
 
Loading...
Searching...
No Matches
field_evaluate.cc
Go to the documentation of this file.
1
21//
22// evaluate a field on a predefined point set: hat_x[q], q=0..nq
23// See also piola_transformation.h
24//
25#include "rheolef/field_evaluate.h"
26#include "rheolef/piola_util.h"
27namespace rheolef {
28
29// TODO: DVT_FEM_CLASS : should move in a
30// fem_on_pointset & fem class
31// fem: we known how to evaluate : RTk, etc
32// fem_on_pointset: loop on a quadrature
33//
34// -----------------------------------------------------
35// scalar-valued case:
36// -----------------------------------------------------
37template<class T, class M>
38void
40 const field_basic<T,M>& uh,
41 const basis_on_pointset<T>& bops,
43 const std::vector<size_t>& dis_idof,
44 Eigen::Matrix<T,Eigen::Dynamic,1>& value)
45{
47 const Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>&
48 phij_xi = bops.template evaluate<T> (hat_K);
49 size_type loc_nnod = phij_xi.rows();
50 size_type loc_ndof = phij_xi.cols();
51 value.resize (loc_nnod);
52 for (size_type loc_jdof = 0; loc_jdof < loc_ndof; ++loc_jdof) {
53 const T& u_jdof = uh.dis_dof (dis_idof[loc_jdof]);
54 for (size_type loc_inod = 0; loc_inod < loc_nnod; ++loc_inod) {
55 value[loc_inod] += phij_xi (loc_inod,loc_jdof)*u_jdof;
56 }
57 }
58}
59// -----------------------------------------------------
60// vector-valued case:
61// -----------------------------------------------------
62template<class T, class M>
63void
65 const field_basic<T,M>& uh,
66 const basis_on_pointset<T>& bops,
68 const std::vector<size_t>& dis_idof_tab,
69 const basis_on_pointset<T>& piola_on_geo_basis,
70 std::vector<size_t>& dis_inod_geo,
71 Eigen::Matrix<point_basic<T>,Eigen::Dynamic,1>& value)
72{
74 size_type dis_ndof = uh.dis_ndof();
75 size_type loc_ndof = dis_idof_tab.size();
76 if (uh.get_space().get_constitution().is_hierarchical()) {
77 // vector: as a product of scalar basis, e.g. (Pk)^d
78 const Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>&
79 phij_xi = bops.template evaluate<T> (hat_K);
80 size_type loc_nnod = phij_xi.rows();
81 size_type loc_comp_ndof = phij_xi.cols();
82 value.resize (loc_nnod);
83 for (size_type loc_inod = 0; loc_inod < loc_nnod; ++loc_inod) {
84 value[loc_inod] = point_basic<T>(0,0,0);
85 }
86 size_type n_comp = uh.get_geo().dimension();
87 for (size_type loc_comp_jdof = 0; loc_comp_jdof < loc_comp_ndof; ++loc_comp_jdof) {
88 for (size_type k_comp = 0; k_comp < n_comp; k_comp++) {
89 size_type loc_jdof = loc_comp_jdof + k_comp*loc_comp_ndof;
90 assert_macro (loc_jdof < loc_ndof, "invalid local index "<<loc_jdof<<" out of range [0:"<<loc_ndof<<"[");
91 size_type dis_jdof = dis_idof_tab[loc_jdof];
92 assert_macro (dis_jdof < dis_ndof, "invalid distr index");
93 const T& u_jdof = uh.dis_dof (dis_jdof);
94 for (size_type loc_inod = 0; loc_inod < loc_nnod; ++loc_inod) {
95 value[loc_inod][k_comp] += phij_xi(loc_inod,loc_jdof)*u_jdof;
96 }
97 }
98 }
99 } else { // non-hierarchical vector basis: e.g. Hdiv RTk requires piola vector transformation
100 check_macro (piola_on_geo_basis.is_set(), "un-set piola_on_geo_basis");
101 const Eigen::Matrix<point_basic<T>,Eigen::Dynamic,Eigen::Dynamic>&
102 phij_xi = bops.template evaluate<point_basic<T>> (hat_K);
103 size_type loc_nnod = phij_xi.rows();
104 size_type loc_ndof = phij_xi.cols();
105 value.resize (loc_nnod);
106 for (size_type loc_inod = 0; loc_inod < loc_nnod; ++loc_inod) {
107 value[loc_inod] = point_basic<T>(0,0,0);
108 }
109 for (size_type loc_jdof = 0; loc_jdof < loc_ndof; ++loc_jdof) {
110 size_type dis_jdof = dis_idof_tab[loc_jdof];
111 assert_macro (dis_jdof < dis_ndof, "invalid distributed index");
112 const T& u_jdof = uh.dis_dof (dis_jdof);
113 for (size_type loc_inod = 0; loc_inod < loc_nnod; ++loc_inod) {
114 value[loc_inod] += phij_xi(loc_inod,loc_jdof)*u_jdof;
115 }
116 }
117 // piola_on_geo_basis : on evalue la base de la transf geo aux points xq
118 // c'est different de bops : evalue la base fem_RT aux points xq
119 // TODO: l'initialiser en amont de la boucle en q, dans la classe field_evaluate_base
120 // si on est vecteur et non-hierarchical
121 // basis_on_pointset<T> piola_on_geo_basis;
122 // piola_on_geo_basis.set (b, omega.get_piola_basis());
123 // std::vector<size_type> dis_inod_geo;
124 Eigen::Matrix<tensor_basic<T>,Eigen::Dynamic,1> DF;
125 jacobian_piola_transformation (uh.get_geo(), piola_on_geo_basis, hat_K, dis_inod_geo, DF);
126 for (size_type loc_inod = 0; loc_inod < loc_nnod; ++loc_inod) {
127 T J = det_jacobian_piola_transformation (DF[loc_inod], uh.get_geo().dimension(), hat_K.dimension());
128 value[loc_inod] = (1/J)*(DF[loc_inod]*value[loc_inod]); // see RavTho-1977, eqn (3.17)
129 }
130 }
131}
132// -----------------------------------------------------
133// tensor-valued case:
134// -----------------------------------------------------
135template<class T, class M>
136void
138 const field_basic<T,M>& uh,
139 const basis_on_pointset<T>& bops,
140 reference_element hat_K,
141 const std::vector<size_t>& dis_idof_tab,
142 Eigen::Matrix<tensor_basic<T>,Eigen::Dynamic,1>& value)
143{
144 typedef typename field_basic<T,M>::size_type size_type;
145 const Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>&
146 phij_xi = bops.template evaluate<T> (hat_K);
147 size_type loc_nnod = phij_xi.rows();
148 size_type loc_comp_ndof = phij_xi.cols();
149 value.resize (loc_nnod);
150 for (size_type loc_inod = 0; loc_inod < loc_nnod; ++loc_inod) {
151 value[loc_inod] = tensor_basic<T>();
152 }
153 size_type dis_ndof = uh.dis_ndof();
154 size_type d = uh.get_geo().dimension();
155 space_constant::coordinate_type sys_coord = uh.get_geo().coordinate_system();
157 for (size_type loc_comp_jdof = 0; loc_comp_jdof < loc_comp_ndof; ++loc_comp_jdof) {
158 for (size_type kl_comp = 0; kl_comp < n_comp; kl_comp++) {
159 size_type loc_jdof = loc_comp_jdof + kl_comp*loc_comp_ndof;
160 assert_macro (loc_jdof < loc_ndof, "invalid local index "<<loc_jdof<<" out of range [0:"<<loc_ndof<<"[");
161 size_type dis_jdof = dis_idof_tab[loc_jdof];
162 assert_macro (dis_jdof < dis_ndof, "invalid distr index");
163 const T& u_jdof = uh.dis_dof (dis_jdof);
164 std::pair<size_type,size_type> kl
166 size_type k = kl.first;
167 size_type l = kl.second;
168 for (size_type loc_inod = 0; loc_inod < loc_nnod; ++loc_inod) {
169 T tmp = phij_xi (loc_inod, loc_comp_jdof)*u_jdof;
170 value[loc_inod](k,l) += tmp;
171 if (k != l) value[loc_inod](l,k) += tmp;
172 }
173 }
174 }
175}
176// -----------------------------------------------------
177// homogeneous multi-component case: get the i-th value
178// -----------------------------------------------------
179template<class T, class M>
180void
182 const field_basic<T,M>& uh,
183 const basis_on_pointset<T>& bops,
184 reference_element hat_K,
185 const std::vector<size_t>& dis_idof_tab,
186 size_t k_comp,
187 Eigen::Matrix<T,Eigen::Dynamic,1>& value)
188{
189 typedef typename field_basic<T,M>::size_type size_type;
190 const Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>&
191 phij_xi = bops.template evaluate<T> (hat_K);
192 size_type loc_nnod = phij_xi.rows();
193 size_type loc_comp_ndof = phij_xi.cols();
194 value.resize (loc_nnod);
195 for (size_type loc_inod = 0; loc_inod < loc_nnod; ++loc_inod) {
196 value[loc_inod] = 0;
197 }
198 size_type dis_ndof = uh.dis_ndof();
199 size_type loc_ndof = dis_idof_tab.size();
200 size_type n_comp = uh.get_geo().dimension();
201 for (size_type loc_comp_jdof = 0; loc_comp_jdof < loc_comp_ndof; ++loc_comp_jdof) {
202 size_type loc_jdof = loc_comp_jdof + k_comp*loc_comp_ndof;
203 assert_macro (loc_jdof < loc_ndof, "invalid local index "<<loc_jdof<<" out of range [0:"<<loc_ndof<<"[");
204 size_type dis_jdof = dis_idof_tab[loc_jdof];
205 assert_macro (dis_jdof < dis_ndof, "invalid distr index");
206 const T& u_jdof = uh.dis_dof (dis_jdof);
207 for (size_type loc_inod = 0; loc_inod < loc_nnod; ++loc_inod) {
208 value[loc_inod] += phij_xi (loc_inod,loc_comp_jdof)*u_jdof;
209 }
210 }
211}
212// ----------------------------------------------------------------------------
213// evaluate on a pointset: new version wih fops
214// ----------------------------------------------------------------------------
215template<class T, class M, class Value>
216void
218 const field_basic<T,M>& uh,
219 const geo_basic<T,M>& omega_K,
220 const geo_element& K,
221 const Eigen::Matrix<Value,Eigen::Dynamic,Eigen::Dynamic>& phij_xi,
222 Eigen::Matrix<Value,Eigen::Dynamic,1>& value)
223{
224 using size_type = typename field_basic<T,M>::size_type;
225 // 1) extract uh dofs on K
226 // TODO: merge all cases in assembly_dis_idof instead of here
227 std::vector<size_type> dis_idof;
228 uh.get_space().get_constitution().assembly_dis_idof (omega_K, K, dis_idof);
229 size_type loc_ndof = dis_idof.size();
230 Eigen::Matrix<T,Eigen::Dynamic,1> udof (loc_ndof);
231 check_macro (size_type(phij_xi.cols()) == loc_ndof,
232 "invalid sizes phij_xi("<<phij_xi.rows()<<","<<phij_xi.cols()<<") and udof("<<udof.size()<<")");
233 for (size_type loc_jdof = 0; loc_jdof < loc_ndof; ++loc_jdof) {
234 size_type dis_jdof = dis_idof[loc_jdof];
235 udof [loc_jdof] = uh.dis_dof (dis_idof[loc_jdof]);
236 }
237 // 2) interpolate uh on the prescribed pointset
238 // TODO: DVT_EIGEN_BLAS2: value = phij_xi*udof; but Eigen compil pb when Value=point
239 size_type loc_nnod = phij_xi.rows();
240 value.resize (loc_nnod);
241 for (size_type loc_inod = 0; loc_inod < loc_nnod; ++loc_inod) {
242 Value sum = Value();
243 for (size_type loc_jdof = 0; loc_jdof < loc_ndof; ++loc_jdof) {
244 sum += phij_xi (loc_inod,loc_jdof)*udof[loc_jdof];
245 }
246 value[loc_inod] = sum;
247 }
248}
249template<class T, class M, class Value>
250void
252 const field_basic<T,M>& uh,
253 const fem_on_pointset<T>& fops,
254 const geo_basic<T,M>& omega_K,
255 const geo_element& K,
256 Eigen::Matrix<Value,Eigen::Dynamic,1>& value)
257{
258 Eigen::Matrix<Value,Eigen::Dynamic,Eigen::Dynamic> phij_xi;
259 details::differentiate_option none;
260 fops.template evaluate <M,Value,details::differentiate_option::none> (omega_K, K, none, phij_xi);
261 field_evaluate_continued (uh, omega_K, K, phij_xi, value);
262}
263// ----------------------------------------------------------------------------
264// instanciation in library
265// ----------------------------------------------------------------------------
266#define _RHEOLEF_instanciation_value(T,M,Value) \
267template void field_evaluate_continued ( \
268 const field_basic<T,M>& uh, \
269 const geo_basic<T,M>& omega_K, \
270 const geo_element& K, \
271 const Eigen::Matrix<Value,Eigen::Dynamic,Eigen::Dynamic>& phij_xi, \
272 Eigen::Matrix<Value,Eigen::Dynamic,1>& value); \
273template void field_evaluate ( \
274 const field_basic<T,M>& uh, \
275 const fem_on_pointset<T>& fops, \
276 const geo_basic<T,M>& omega_K, \
277 const geo_element& K, \
278 Eigen::Matrix<Value,Eigen::Dynamic,1>& value); \
279
280
281#define _RHEOLEF_instanciation(T,M) \
282_RHEOLEF_instanciation_value(T,M,T) \
283_RHEOLEF_instanciation_value(T,M,point_basic<T>) \
284_RHEOLEF_instanciation_value(T,M,tensor_basic<T>) \
285_RHEOLEF_instanciation_value(T,M,tensor3_basic<T>) \
286_RHEOLEF_instanciation_value(T,M,tensor4_basic<T>) \
287template \
288void \
289field_evaluate ( \
290 const field_basic<T,M>& uh, \
291 const basis_on_pointset<T>& bops, \
292 reference_element hat_K, \
293 const std::vector<size_t>& dis_idof, \
294 Eigen::Matrix<T,Eigen::Dynamic,1>& value); \
295template \
296void \
297vector_field_evaluate<T,M> ( \
298 const field_basic<T,M>& uh, \
299 const basis_on_pointset<T>& bops, \
300 reference_element hat_K, \
301 const std::vector<size_t>& dis_idof_tab, \
302 const basis_on_pointset<T>& piola_on_geo_basis, \
303 std::vector<size_t>& dis_inod_geo, \
304 Eigen::Matrix<point_basic<T>,Eigen::Dynamic,1>& value); \
305template \
306void \
307tensor_field_evaluate<T,M> ( \
308 const field_basic<T,M>& uh, \
309 const basis_on_pointset<T>& bops, \
310 reference_element hat_K, \
311 const std::vector<size_t>& dis_idof_tab, \
312 Eigen::Matrix<tensor_basic<T>,Eigen::Dynamic,1>& value); \
313template \
314void \
315field_component_evaluate<T,M> ( \
316 const field_basic<T,M>& uh, \
317 const basis_on_pointset<T>& bops, \
318 reference_element hat_K, \
319 const std::vector<size_t>& dis_idof_tab, \
320 size_t i_comp, \
321 Eigen::Matrix<T,Eigen::Dynamic,1>& value); \
322
323
325#ifdef _RHEOLEF_HAVE_MPI
326_RHEOLEF_instanciation(Float,distributed)
327#endif // _RHEOLEF_HAVE_MPI
328
329} // 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
const T & dis_dof(size_type dis_idof) const
Definition field.cc:377
size_type dis_ndof() const
Definition field.h:299
const geo_type & get_geo() const
Definition field.h:271
const space_type & get_space() const
Definition field.h:270
std::size_t size_type
Definition field.h:225
generic mesh with rerefence counting
Definition geo.h:1089
see the geo_element page for the full documentation
see the reference_element page for the full documentation
#define assert_macro(ok_condition, message)
Definition dis_macros.h:113
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)")
size_type n_component(valued_type valued_tag, size_type d, coordinate_type sys_coord)
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.
void field_evaluate(const field_basic< T, M > &uh, const basis_on_pointset< T > &bops, reference_element hat_K, const std::vector< size_t > &dis_idof, Eigen::Matrix< T, Eigen::Dynamic, 1 > &value)
void tensor_field_evaluate(const field_basic< T, M > &uh, const basis_on_pointset< T > &bops, reference_element hat_K, const std::vector< size_t > &dis_idof_tab, Eigen::Matrix< tensor_basic< T >, Eigen::Dynamic, 1 > &value)
void jacobian_piola_transformation(const geo_basic< T, M > &omega, const basis_on_pointset< T > &piola_on_pointset, reference_element hat_K, const std::vector< size_t > &dis_inod, Eigen::Matrix< tensor_basic< T >, Eigen::Dynamic, 1 > &DF)
Definition piola_util.cc:80
void field_evaluate_continued(const field_basic< T, M > &uh, const geo_basic< T, M > &omega_K, const geo_element &K, const Eigen::Matrix< Value, Eigen::Dynamic, Eigen::Dynamic > &phij_xi, Eigen::Matrix< Value, Eigen::Dynamic, 1 > &value)
T det_jacobian_piola_transformation(const tensor_basic< T > &DF, size_t d, size_t map_d)
void vector_field_evaluate(const field_basic< T, M > &uh, const basis_on_pointset< T > &bops, reference_element hat_K, const std::vector< size_t > &dis_idof_tab, const basis_on_pointset< T > &piola_on_geo_basis, std::vector< size_t > &dis_inod_geo, Eigen::Matrix< point_basic< T >, Eigen::Dynamic, 1 > &value)
void field_component_evaluate(const field_basic< T, M > &uh, const basis_on_pointset< T > &bops, reference_element hat_K, const std::vector< size_t > &dis_idof_tab, size_t k_comp, Eigen::Matrix< T, Eigen::Dynamic, 1 > &value)