Rheolef  7.1
an efficient C++ finite element environment
interpolate.cc
Go to the documentation of this file.
1 #include "rheolef/geo_domain.h"
22 #include "rheolef/interpolate.h"
23 #include "rheolef/space_numbering.h"
24 #include "rheolef/piola_util.h"
25 
26 #ifdef _RHEOLEF_HAVE_MPI
27 #include "rheolef/mpi_scatter_init.h"
28 #include "rheolef/mpi_scatter_begin.h"
29 #include "rheolef/mpi_scatter_end.h"
30 #endif // _RHEOLEF_HAVE_MPI
31 
32 namespace rheolef {
33 
34 // ============================================================================
35 // specialized interpolate algorithm:
36 // interpolate one field on disarray of nodes with element locations
37 // and put the result into an disarray of values
38 // => used by characteristic
39 // ============================================================================
40 namespace details {
41 
42 template<class T, class M>
43 void
45  const geo_basic<T,M>& omega1,
46  const disarray<point_basic<T>,M>& x2,
47  const disarray<geo_element::size_type,M>& ix2_2_dis_ie1, // x2 has been already localized in K1
48  disarray<index_set,M>& ie1_2_dis_ix2, // K1 -> list of ix2
49  disarray<point_basic<T>,M>& hat_x2)
50 {
51  typedef typename field_basic<T,M>::size_type size_type;
52  // -----------------------------------------------------------------------------
53  // 1) in order to call only one time per K element the dis_inod() function,
54  // we group the x2(i) per K element of omega1
55  // -----------------------------------------------------------------------------
56  size_type first_dis_ix2 = x2.ownership().first_index();
57  distributor ie1_ownership = omega1.geo_element_ownership (omega1.map_dimension());
58  size_type first_dis_ie1 = ie1_ownership.first_index();
59  std::set<size_type> ext_ie1_set;
61  ie1_2_dis_ix2.resize (ie1_ownership, empty_set);
62  for (size_type ix2 = 0, nx2 = ix2_2_dis_ie1.size(); ix2 < nx2; ix2++) {
63  size_type dis_ie1 = ix2_2_dis_ie1 [ix2];
64  size_type dis_ix2 = first_dis_ix2 + ix2;
65  index_set dis_ix2_set; dis_ix2_set.insert (dis_ix2);
66  check_macro (dis_ie1 != std::numeric_limits<size_type>::max(), "node("<<ix2<<")="<<x2[ix2]
67  <<" cannot be localized in "<<omega1.name() << " (HINT: curved geometry ?)");
68  ie1_2_dis_ix2.dis_entry (dis_ie1) += dis_ix2_set;
69  }
70  ie1_2_dis_ix2.dis_entry_assembly();
71  // -----------------------------------------------------------------------------
72  // 2) loop on K and list external x(i)
73  // -----------------------------------------------------------------------------
74  // TODO: it would be more efficient to send (dis_ix,x) together : less mpi calls
75  // => how to do an assembly() with a disarray<pair<ix,x> > ?
76  // ... it does not support it yet
77  // also: would use more memory: massive coordinates copy...
78  communicator comm = ie1_ownership.comm();
79  distributor ix2_ownership = x2.ownership();
80  index_set ext_dis_ix2_set;
81  if (is_distributed<M>::value && comm.size() > 1) {
82  for (size_type ie1 = 0, ne1 = ie1_2_dis_ix2.size(); ie1 < ne1; ie1++) {
83  const index_set& dis_ix2_set = ie1_2_dis_ix2 [ie1];
84  for (typename index_set::const_iterator iter = dis_ix2_set.begin(), last = dis_ix2_set.end(); iter != last; ++iter) {
85  size_type dis_ix2 = *iter;
86  if (ix2_ownership.is_owned (dis_ix2)) continue;
87  ext_dis_ix2_set.insert (dis_ix2);
88  }
89  }
90  // TODO: when x=xdofs lives in space, then it modifies (enlarges)
91  // definitively the set of externals indexes & values
92  // => massive permanent storage
93  // otherwise: manage the scatter in a separate map<ext_ix,x>
94  // map<ix,x> ext_x_map;
95  // x.append_dis_indexes (ext_dis_ix_set, ext_x_map);
96  // bug: see field_reinterpolate2_tst.cc
97  // x.set_dis_indexes (ext_dis_ix_set);
98  x2.append_dis_indexes (ext_dis_ix2_set);
99  }
100  // -----------------------------------------------------------------------------
101  // 3) loop on K and compute hat_x(dis_ix)
102  // -----------------------------------------------------------------------------
103  T infty = std::numeric_limits<T>::max();
104  hat_x2.resize (ix2_ownership, point_basic<T>(infty,infty,infty));
105  std::vector<size_type> dis_inod1;
106  for (size_type ie1 = 0, ne1 = ie1_2_dis_ix2.size(); ie1 < ne1; ie1++) {
107  const index_set& dis_ix2_set = ie1_2_dis_ix2 [ie1];
108  if (dis_ix2_set.size() == 0) continue;
109  const geo_element& K1 = omega1 [ie1];
110  omega1.dis_inod (K1, dis_inod1);
111  for (typename index_set::const_iterator iter = dis_ix2_set.begin(), last = dis_ix2_set.end(); iter != last; ++iter) {
112  size_type dis_ix2 = *iter;
113  const point_basic<T>& xi2 = x2.dis_at (dis_ix2);
114  hat_x2.dis_entry (dis_ix2) = inverse_piola_transformation (omega1, K1, dis_inod1, xi2);
115  }
116  }
117  hat_x2.dis_entry_assembly();
118  hat_x2.set_dis_indexes (ext_dis_ix2_set);
119 }
120 template<class T, class M, class Value>
121 void
123  const field_basic<T,M>& u1h,
124  const disarray<point_basic<T>,M>& x2,
125  const disarray<index_set,M>& ie1_2_dis_ix2, // K1 -> list of ix2
126  const disarray<point_basic<T>,M>& hat_x2, // ix2 -> hat_x2
127  disarray<Value,M>& ux2)
128 {
129  typedef typename field_basic<T,M>::size_type size_type;
130  u1h.dis_dof_update();
131  T infty = std::numeric_limits<T>::max();
132  distributor x2_ownership = x2.ownership();
133  ux2.resize (x2.ownership(), Value(infty));
134  // -----------------------------------------------------------------------------
135  // on locally managed K1, evaluate at hat_x2 that could be external
136  // -----------------------------------------------------------------------------
137  const geo_basic<T,M>& omega1 = u1h.get_geo();
138  const space_basic<T,M>& Xh1 = u1h.get_space();
139  const basis_basic<T>& b1 = u1h.get_space().get_basis();
140  check_macro (! b1.get_piola_fem().transform_need_piola(), "unsupported piola-based basis");
141  std::vector<size_type> dis_idof1;
142  std::vector<T> dof1;
143  Eigen::Matrix<Value,Eigen::Dynamic,1> b1_value;
144  for (size_type ie1 = 0, ne1 = ie1_2_dis_ix2.size(); ie1 < ne1; ie1++) {
145  const index_set& dis_ix2_set = ie1_2_dis_ix2 [ie1];
146  if (dis_ix2_set.size() == 0) continue;
147  // extract dis_idof[] & dof[] one time for all on K1
148  const geo_element& K1 = omega1 [ie1];
149  Xh1.dis_idof (K1, dis_idof1);
150  size_type loc_ndof1 = dis_idof1.size();
151  dof1.resize(loc_ndof1);
152  for (size_type loc_idof1 = 0; loc_idof1 < loc_ndof1; loc_idof1++) {
153  dof1 [loc_idof1] = u1h.dis_dof (dis_idof1 [loc_idof1]);
154  }
155  b1_value.resize (loc_ndof1);
156  // loop on all hat_x in K: eval basis at hat_x and combine with dof[]
157  for (typename index_set::const_iterator iter = dis_ix2_set.begin(), last = dis_ix2_set.end(); iter != last; ++iter) {
158  size_type dis_ix2 = *iter;
159  size_type iproc = x2_ownership.find_owner (dis_ix2);
160  size_type nx2 = x2_ownership.size (iproc);
161  size_type first_dis_ix2 = x2_ownership.first_index(iproc);
162  size_type ix2 = dis_ix2 - first_dis_ix2;
163  const point_basic<T>& hat_xi = hat_x2.dis_at (dis_ix2);
164  b1.evaluate (K1, hat_xi, b1_value); // TODO: RTk and Hdiv: apply also Piola transform at hat_xi
165  Value value = Value();
166  for (size_type loc_idof1 = 0; loc_idof1 < loc_ndof1; loc_idof1++) {
167  value += dof1 [loc_idof1] * b1_value[loc_idof1]; // sum_i w_coef(i)*hat_phi(hat_x)
168  }
169  ux2.dis_entry (dis_ix2) = value;
170  }
171  }
172  ux2.dis_entry_assembly();
173 }
174 template<class T, class M, class Value>
175 void
177  const disarray<Value,M>& ux2,
178  field_basic<T,M>& u2h)
179 {
180  typedef typename field_basic<T,M>::size_type size_type;
181  std::set<size_type> ext_dis_inod_set;
182  u2h.get_space().get_xdofs().get_dis_indexes (ext_dis_inod_set);
183  ux2.set_dis_indexes (ext_dis_inod_set);
184 
185  const space_basic<T,M>& V2h = u2h.get_space();
186  const geo_basic<T,M>& omega2 = u2h.get_geo();
187  const basis_basic<T>& b2 = V2h.get_basis();
188  Eigen::Matrix<Value,Eigen::Dynamic,1> u_xnod2;
189  Eigen::Matrix<T,Eigen::Dynamic,1> udof2;
190  std::vector<size_type> dis_inod2, dis_idof2;
191  for (size_type ie2 = 0, ne2 = omega2.size(); ie2 < ne2; ++ie2) {
192  const geo_element& K2 = omega2 [ie2];
193  space_numbering::dis_inod (V2h.get_basis(), omega2.sizes(), K2, dis_inod2);
194  u_xnod2.resize (dis_inod2.size());
195  for (size_type loc_inod2 = 0, loc_nnod2 = dis_inod2.size(); loc_inod2 < loc_nnod2; ++loc_inod2) {
196  check_macro (dis_inod2[loc_inod2] < ux2.ownership().dis_size(), "invalid size");
197  u_xnod2 [loc_inod2] = ux2.dis_at (dis_inod2[loc_inod2]);
198  }
199  b2.compute_dofs (K2, u_xnod2, udof2);
200  space_numbering::dis_idof (V2h.get_basis(), omega2.sizes(), K2, dis_idof2);
201  check_macro (dis_idof2.size() == size_type(udof2.size()), "invalid sizes");
202  for (size_type loc_idof2 = 0, loc_ndof2 = dis_idof2.size(); loc_idof2 < loc_ndof2; ++loc_idof2) {
203  u2h.dis_dof_entry (dis_idof2[loc_idof2]) = udof2 [loc_idof2];
204  }
205  }
206  u2h.dis_dof_update();
207 }
208 template<class T, class M, class Value>
209 void
211  const field_basic<T,M>& u1h,
212  const disarray<point_basic<T>,M>& x2,
213  const disarray<geo_element::size_type,M>& ix2dis_ie,
214  disarray<Value,M>& ux2)
215 {
216  const geo_basic<T,M>& omega1 = u1h.get_geo();
217  disarray<index_set,M> ie2dis_ix;
218  disarray<point_basic<T>,M> hat_x;
219  interpolate_pass1_symbolic (omega1, x2, ix2dis_ie, ie2dis_ix, hat_x);
220  interpolate_pass2_valued (u1h, x2, ie2dis_ix, hat_x, ux2);
221 }
222 
223 } // namespace details
224 
225 // ============================================================================
226 // interpolate function:
227 // re-interpolate one field on another mesh and space:
228 // field u2h = interpolate (V2h, u1h);
229 // ============================================================================
231 template<class T, class M>
232 field_basic<T,M>
234 {
235  u1h.dis_dof_update();
236  u1h.get_geo().boundary(); // force "boundary" domain creation, for geo locator
237  typedef typename field_basic<T,M>::size_type size_type;
238  if (u1h.get_space() == V2h) {
239  size_type have_same_dofs = (u1h.u().size() == V2h.iu_ownership().size());
240 #ifdef _RHEOLEF_HAVE_MPI
242  have_same_dofs = mpi::all_reduce (V2h.ownership().comm(), have_same_dofs, mpi::minimum<size_type>());
243  }
244 #endif // _RHEOLEF_HAVE_MPI
245  if (have_same_dofs) {
246  // spaces are exactly the same: no need to re-interpolate or copy
247  return u1h;
248  }
249  // spaces differs only by blocked/unblocked dofs: need to copy
250  field_basic<T,M> u2h (V2h);
251  for (size_type idof = 0, ndof = V2h.ndof(); idof < ndof; ++idof) {
252  u2h.dof(idof) = u1h.dof(idof);
253  }
254  u2h.dis_dof_update();
255  return u2h;
256  }
257  if (u1h.get_geo().get_background_geo() == V2h.get_geo().get_background_geo()) {
258  // meshes are compatible:
259  // => buid a wrapper expression and call back interpolate with the
260  // general nonlinear specialized version
262  }
263  // -----------------------------------------------------------------------------
264  // between different meshes:
265  // 0) create the xdofs2 set of nodes
266  // -----------------------------------------------------------------------------
267  const disarray<point_basic<T>,M>& xdof2 = V2h.get_xdofs();
268  // -----------------------------------------------------------------------------
269  // 1) locate each xdof of V2 in omega1
270  // -----------------------------------------------------------------------------
271  disarray<point_basic<T>,M> xdof2_nearest (xdof2.ownership());
272  const geo_basic<T,M>& omega1 = u1h.get_geo();
273  disarray<size_type,M> dis_ie1_tab;
274  omega1.nearest (xdof2, xdof2_nearest, dis_ie1_tab);
275  // -----------------------------------------------------------------------------
276  // 2) interpolate uh1 at xdof2_nearest: get the values in a disarray
277  // -----------------------------------------------------------------------------
278  field_basic<T,M> u2h (V2h);
279  switch (V2h.valued_tag()) {
280  case space_constant::scalar: {
281  disarray<T,M> u2h_dof;
282  details::interpolate_on_a_different_mesh (u1h, xdof2_nearest, dis_ie1_tab, u2h_dof);
283  details::interpolate_pass3_dof (u2h_dof, u2h);
284  return u2h;
285  }
286  case space_constant::vector: {
287  disarray<point_basic<T>,M> u2h_dof;
288  details::interpolate_on_a_different_mesh (u1h, xdof2_nearest, dis_ie1_tab, u2h_dof);
289  details::interpolate_pass3_dof (u2h_dof, u2h);
290  return u2h;
291  }
294  disarray<tensor_basic<T>,M> u2h_dof;
295  details::interpolate_on_a_different_mesh (u1h, xdof2_nearest, dis_ie1_tab, u2h_dof);
296  details::interpolate_pass3_dof (u2h_dof, u2h);
297  return u2h;
298  }
299  default:
300  error_macro("interpolate between different meshes: unexpected "<<V2h.valued()<<"-valued field");
301  }
302  return field_basic<T,M>();
303 }
304 // ----------------------------------------------------------------------------
305 // instanciation in library
306 // ----------------------------------------------------------------------------
307 #define _RHEOLEF_instanciation_value(T,M,Value) \
308 template void details::interpolate_pass2_valued ( \
309  const field_basic<T,M>& uh, \
310  const disarray<point_basic<T>,M>& x, \
311  const disarray<index_set,M>& ie2dis_ix, \
312  const disarray<point_basic<T>,M>& hat_x, \
313  disarray<Value,M>& ux); \
314 
315 #define _RHEOLEF_instanciation(T,M) \
316 _RHEOLEF_instanciation_value(T,M,T) \
317 _RHEOLEF_instanciation_value(T,M,point_basic<T>) \
318 template \
319 void \
320 details::interpolate_pass1_symbolic ( \
321  const geo_basic<T,M>&, \
322  const disarray<point_basic<T>,M>&, \
323  const disarray<geo_element::size_type,M>&, \
324  disarray<index_set,M>&, \
325  disarray<point_basic<T>,M>&); \
326 template \
327 field_basic<T,M> \
328 interpolate (const space_basic<T,M>&, const field_basic<T,M>&);
329 
330 _RHEOLEF_instanciation(Float,sequential)
331 #ifdef _RHEOLEF_HAVE_MPI
333 #endif // _RHEOLEF_HAVE_MPI
334 
335 } // 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::index_set::insert
void insert(size_type dis_i)
Definition: index_set_header.icc:158
rheolef::geo_basic
generic mesh with rerefence counting
Definition: geo.h:1089
Xh1
space_basic< T, M > Xh1
Definition: field_expr.h:263
rheolef::field_basic::get_space
const space_type & get_space() const
Definition: field.h:300
rheolef::basis_basic::get_piola_fem
const piola_fem< T > & get_piola_fem() const
Definition: basis.h:795
rheolef::details::interpolate_pass3_dof
void interpolate_pass3_dof(const disarray< Value, M > &ux2, field_basic< T, M > &u2h)
Definition: interpolate.cc:176
rheolef::point_basic
Definition: point.h:87
rheolef::basis_basic::evaluate
void evaluate(reference_element hat_K, const point_basic< T > &hat_x, Eigen::Matrix< Value, Eigen::Dynamic, 1 > &value) const
Definition: basis.h:852
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::details::interpolate_pass1_symbolic
void interpolate_pass1_symbolic(const geo_basic< T, M > &omega, const disarray< point_basic< T >, M > &x, const disarray< geo_element::size_type, M > &ix2dis_ie, disarray< index_set, M > &ie2dis_ix, disarray< point_basic< T >, M > &hat_y)
Definition: interpolate.cc:44
rheolef::distributor::comm
const communicator_type & comm() const
Definition: distributor.h:145
rheolef::details::interpolate_pass2_valued
void interpolate_pass2_valued(const field_basic< T, M > &uh, const disarray< point_basic< T >, M > &x, const disarray< index_set, M > &ie2dis_ix, const disarray< point_basic< T >, M > &hat_y, disarray< Value, M > &ux)
Definition: interpolate.cc:122
rheolef::distributor::find_owner
size_type find_owner(size_type dis_i) const
find iproc associated to a global index dis_i: CPU=log(nproc)
Definition: distributor.cc:106
rheolef::_RHEOLEF_instanciation
_RHEOLEF_instanciation(Float, sequential, std::allocator< Float >) _RHEOLEF_instanciation(Float
rheolef::value
rheolef::std value
rheolef::space_basic
the finite element space
Definition: space.h:352
rheolef::distributor
see the distributor page for the full documentation
Definition: distributor.h:62
rheolef::geo_element
see the geo_element page for the full documentation
Definition: geo_element.h:102
rheolef::index_set
Definition: index_set_header.icc:58
rheolef::size_type
size_t size_type
Definition: basis_get.cc:76
rheolef::space_numbering::dis_idof
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)
Definition: space_numbering.cc:187
rheolef::details::interpolate_on_a_different_mesh
void interpolate_on_a_different_mesh(const field_basic< T, M > &u1h, const disarray< point_basic< T >, M > &x2, const disarray< geo_element::size_type, M > &ix2dis_ie, disarray< Value, M > &ux2)
Definition: interpolate.cc:210
rheolef::field_basic::dis_dof
const T & dis_dof(size_type dis_idof) const
Definition: field.cc:376
rheolef::empty_set
static std::set< size_t > empty_set
Definition: space_seq.cc:59
rheolef::field_basic::u
const vec< T, M > & u() const
Definition: field.h:310
rheolef::is_distributed
Definition: distributed.h:36
rheolef::basis_basic::compute_dofs
void compute_dofs(reference_element hat_K, const Eigen::Matrix< Value, Eigen::Dynamic, 1 > &f_xnod, Eigen::Matrix< T, Eigen::Dynamic, 1 > &dof) const
Definition: basis.h:897
rheolef::field_basic::dis_dof_entry
dis_reference dis_dof_entry(size_type dis_idof)
Definition: field.cc:397
rheolef::basis_basic
Definition: basis.h:532
rheolef::field_basic::get_geo
const geo_type & get_geo() const
Definition: field.h:301
rheolef::space_constant::tensor
@ tensor
Definition: space_constant.h:138
rheolef::space_constant::scalar
@ scalar
Definition: space_constant.h:136
rheolef::inverse_piola_transformation
point_basic< T > inverse_piola_transformation(const geo_basic< T, M > &omega, const reference_element &hat_K, const std::vector< size_t > &dis_inod, const point_basic< T > &x)
Definition: piola_util.cc:459
rheolef::field_basic::dis_dof_update
void dis_dof_update() const
Definition: field.cc:410
rheolef::interpolate
field_basic< T, M > interpolate(const space_basic< T, M > &V2h, const field_basic< T, M > &u1h)
see the interpolate page for the full documentation
Definition: interpolate.cc:233
rheolef::field_basic
Definition: field.h:235
rheolef::field_basic::dof
T & dof(size_type idof)
Definition: field.h:658
rheolef::space_constant::unsymmetric_tensor
@ unsymmetric_tensor
Definition: space_constant.h:139
rheolef
This file is part of Rheolef.
Definition: compiler_eigen.h:37
error_macro
#define error_macro(message)
Definition: dis_macros.h:49
rheolef::space_numbering::dis_inod
void dis_inod(const basis_basic< T > &b, const geo_size &gs, const geo_element &K, typename std::vector< size_type >::iterator dis_inod_tab)
Definition: space_numbering.cc:177
Float
see the Float page for the full documentation
rheolef::disarray
see the disarray page for the full documentation
Definition: disarray.h:459
rheolef::space_constant::vector
@ vector
Definition: space_constant.h:137
size_type
field::size_type size_type
Definition: branch.cc:425
rheolef::details::field_expr_v2_nonlinear_terminal_field< T, M, details::differentiate_option::none >
rheolef::distributor::is_owned
bool is_owned(size_type dis_i, size_type iproc) const
true when dis_i in [first_index(iproc):last_index(iproc)[
Definition: distributor.h:213
rheolef::distributor::first_index
size_type first_index(size_type iproc) const
global index range and local size owned by ip-th process
Definition: distributor.h:151
rheolef::distributed
distributed
Definition: asr.cc:228
M
Expr1::memory_type M
Definition: vec_expr_v2.h:416
rheolef::distributor::size
size_type size(size_type iproc) const
Definition: distributor.h:163
rheolef::field_basic::size_type
std::size_t size_type
Definition: field.h:239
T
Expr1::float_type T
Definition: field_expr.h:261