1 #include "rheolef/geo_domain.h"
22 #include "rheolef/interpolate.h"
23 #include "rheolef/space_numbering.h"
24 #include "rheolef/piola_util.h"
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
42 template<
class T,
class M>
56 size_type first_dis_ix2 = x2.ownership().first_index();
57 distributor ie1_ownership = omega1.geo_element_ownership (omega1.map_dimension());
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++) {
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;
70 ie1_2_dis_ix2.dis_entry_assembly();
78 communicator comm = ie1_ownership.
comm();
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) {
86 if (ix2_ownership.
is_owned (dis_ix2))
continue;
87 ext_dis_ix2_set.
insert (dis_ix2);
98 x2.append_dis_indexes (ext_dis_ix2_set);
103 T infty = std::numeric_limits<T>::max();
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;
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) {
117 hat_x2.dis_entry_assembly();
118 hat_x2.set_dis_indexes (ext_dis_ix2_set);
120 template<
class T,
class M,
class Value>
131 T infty = std::numeric_limits<T>::max();
133 ux2.resize (x2.ownership(), Value(infty));
141 std::vector<size_type> dis_idof1;
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;
149 Xh1.dis_idof (K1, dis_idof1);
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]);
155 b1_value.resize (loc_ndof1);
157 for (
typename index_set::const_iterator iter = dis_ix2_set.begin(), last = dis_ix2_set.end(); iter != last; ++iter) {
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];
169 ux2.dis_entry (dis_ix2) =
value;
172 ux2.dis_entry_assembly();
174 template<
class T,
class M,
class Value>
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);
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) {
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]);
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];
208 template<
class T,
class M,
class Value>
231 template<
class T,
class M>
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>());
244 #endif // _RHEOLEF_HAVE_MPI
245 if (have_same_dofs) {
252 u2h.
dof(idof) = u1h.
dof(idof);
257 if (u1h.
get_geo().get_background_geo() == V2h.get_geo().get_background_geo()) {
274 omega1.nearest (xdof2, xdof2_nearest, dis_ie1_tab);
279 switch (V2h.valued_tag()) {
300 error_macro(
"interpolate between different meshes: unexpected "<<V2h.valued()<<
"-valued field");
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); \
315 #define _RHEOLEF_instanciation(T,M) \
316 _RHEOLEF_instanciation_value(T,M,T) \
317 _RHEOLEF_instanciation_value(T,M,point_basic<T>) \
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>&); \
328 interpolate (const space_basic<T,M>&, const field_basic<T,M>&);
331 #ifdef _RHEOLEF_HAVE_MPI
333 #endif // _RHEOLEF_HAVE_MPI