Rheolef  7.2
an efficient C++ finite element environment
field_evaluate.cc
Go to the documentation of this file.
1 //
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"
27 namespace 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 // -----------------------------------------------------
37 template<class T, class M>
38 void
40  const field_basic<T,M>& uh,
41  const basis_on_pointset<T>& bops,
42  reference_element hat_K,
43  const std::vector<size_t>& dis_idof,
44  Eigen::Matrix<T,Eigen::Dynamic,1>& value)
45 {
46  typedef typename field_basic<T,M>::size_type size_type;
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 // -----------------------------------------------------
62 template<class T, class M>
63 void
65  const field_basic<T,M>& uh,
66  const basis_on_pointset<T>& bops,
67  reference_element hat_K,
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 {
73  typedef typename field_basic<T,M>::size_type size_type;
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 // -----------------------------------------------------
135 template<class T, class M>
136 void
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 // -----------------------------------------------------
179 template<class T, class M>
180 void
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 // ----------------------------------------------------------------------------
215 template<class T, class M, class Value>
216 void
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 }
249 template<class T, class M, class Value>
250 void
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) \
267 template 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); \
273 template 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>) \
287 template \
288 void \
289 field_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); \
295 template \
296 void \
297 vector_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); \
305 template \
306 void \
307 tensor_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); \
313 template \
314 void \
315 field_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 
324 _RHEOLEF_instanciation(Float,sequential)
325 #ifdef _RHEOLEF_HAVE_MPI
327 #endif // _RHEOLEF_HAVE_MPI
328 
329 } // namespace rheolef
field::size_type size_type
Definition: branch.cc:430
see the Float page for the full documentation
const space_type & get_space() const
Definition: field.h:270
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
see the geo_element page for the full documentation
Definition: geo_element.h:102
see the reference_element page for the full documentation
size_t size_type
Definition: basis_get.cc:76
distributed
Definition: asr.cc:228
point_basic< T >
Definition: piola_fem.h:135
rheolef::std value
#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)")
string sys_coord
Definition: mkgeo_grid.sh:171
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)
void dis_idof(const basis_basic< T > &b, const geo_size &gs, const geo_element &K, typename std::vector< size_type >::iterator dis_idof_tab)
size_type dis_ndof(const basis_basic< T > &b, const geo_size &gs, size_type map_dim)
This file is part of Rheolef.
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)
T det_jacobian_piola_transformation(const tensor_basic< T > &DF, size_t d, size_t map_d)
Definition: piola_util.cc:137
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 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 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)
_RHEOLEF_instanciation(Float, sequential, std::allocator< Float >) _RHEOLEF_instanciation(Float