Rheolef  7.2
an efficient C++ finite element environment
 
Loading...
Searching...
No Matches
interpolate.cc
Go to the documentation of this file.
1
21#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
32namespace 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// ============================================================================
40namespace details {
41
42template<class T, class M>
43void
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{
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;
60 index_set empty_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}
120template<class T, class M, class Value>
121void
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
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}
174template<class T, class M, class Value>
175void
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}
208template<class T, class M, class Value>
209void
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,
215{
216 const geo_basic<T,M>& omega1 = u1h.get_geo();
217 disarray<index_set,M> ie2dis_ix;
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// ============================================================================
231template<class T, class M>
232field_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()) {
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 }
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 }
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) \
308template 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>) \
318template \
319void \
320details::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>&); \
326template \
327field_basic<T,M> \
328interpolate (const space_basic<T,M>&, const field_basic<T,M>&);
329
331#ifdef _RHEOLEF_HAVE_MPI
332_RHEOLEF_instanciation(Float,distributed)
333#endif // _RHEOLEF_HAVE_MPI
334
335} // 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 piola_fem< T > & get_piola_fem() const
Definition basis.h:834
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:964
void evaluate(reference_element hat_K, const point_basic< T > &hat_x, Eigen::Matrix< Value, Eigen::Dynamic, 1 > &value) const
Definition basis.h:919
see the disarray page for the full documentation
Definition disarray.h:497
see the distributor page for the full documentation
Definition distributor.h:69
size_type find_owner(size_type dis_i) const
find iproc associated to a global index dis_i: CPU=log(nproc)
size_type size(size_type iproc) const
size_type first_index(size_type iproc) const
global index range and local size owned by ip-th process
bool is_owned(size_type dis_i, size_type iproc) const
true when dis_i in [first_index(iproc):last_index(iproc)[
const communicator_type & comm() const
const T & dis_dof(size_type dis_idof) const
Definition field.cc:377
T & dof(size_type idof)
Definition field.h:738
void dis_dof_update(const SetOp &=SetOp()) const
Definition field.h:763
dis_reference dis_dof_entry(size_type dis_idof)
Definition field.cc:398
const geo_type & get_geo() const
Definition field.h:271
const space_type & get_space() const
Definition field.h:270
const vec< T, M > & u() const
Definition field.h:282
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
void insert(size_type dis_i)
the finite element space
Definition space.h:382
#define error_macro(message)
Definition dis_macros.h:49
Expr1::float_type T
Definition field_expr.h:230
space_basic< T, M > Xh1
Definition field_expr.h:232
check_macro(expr1.have_homogeneous_space(Xh1), "dual(expr1,expr2); expr1 should have homogeneous space. HINT: use dual(interpolate(Xh, expr1),expr2)")
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)
void interpolate_pass3_dof(const disarray< Value, M > &ux2, field_basic< T, M > &u2h)
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)
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)
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)
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)
This file is part of Rheolef.
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
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)
Expr1::memory_type M