21#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"
42template<
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);
120template<
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();
174template<
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];
208template<
class T,
class M,
class Value>
231template<
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>());
245 if (have_same_dofs) {
251 for (
size_type idof = 0, ndof = V2h.ndof(); idof < ndof; ++idof) {
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) \
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); \
315#define _RHEOLEF_instanciation(T,M) \
316_RHEOLEF_instanciation_value(T,M,T) \
317_RHEOLEF_instanciation_value(T,M,point_basic<T>) \
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>&); \
328interpolate (const space_basic<T,M>&, const field_basic<T,M>&);
331#ifdef _RHEOLEF_HAVE_MPI
#define _RHEOLEF_instanciation(T, M, A)
field::size_type size_type
see the Float page for the full documentation
const piola_fem< T > & get_piola_fem() const
void compute_dofs(reference_element hat_K, const Eigen::Matrix< Value, Eigen::Dynamic, 1 > &f_xnod, Eigen::Matrix< T, Eigen::Dynamic, 1 > &dof) const
void evaluate(reference_element hat_K, const point_basic< T > &hat_x, Eigen::Matrix< Value, Eigen::Dynamic, 1 > &value) const
see the disarray page for the full documentation
see the distributor page for the full documentation
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
void dis_dof_update(const SetOp &=SetOp()) const
dis_reference dis_dof_entry(size_type dis_idof)
const geo_type & get_geo() const
const space_type & get_space() const
const vec< T, M > & u() const
generic mesh with rerefence counting
see the geo_element page for the full documentation
void insert(size_type dis_i)
#define error_macro(message)
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)