Loading...
Searching...
No Matches
interpolate.cc
Go to the documentation of this file.
#include "rheolef/geo_domain.h"
#include "rheolef/interpolate.h"
#include "rheolef/space_numbering.h"
#include "rheolef/piola_util.h"
#ifdef _RHEOLEF_HAVE_MPI
#include "rheolef/mpi_scatter_init.h"
#include "rheolef/mpi_scatter_begin.h"
#include "rheolef/mpi_scatter_end.h"
#endif // _RHEOLEF_HAVE_MPI
namespace rheolef {
// ============================================================================
// specialized interpolate algorithm:
// interpolate one field on disarray of nodes with element locations
// and put the result into an disarray of values
// => used by characteristic
// ============================================================================
namespace details {
template<class T, class M>
void
interpolate_pass1_symbolic (
const geo_basic<T,M>& omega1,
const disarray<point_basic<T>,M>& x2,
const disarray<geo_element::size_type,M>& ix2_2_dis_ie1, // x2 has been already localized in K1
disarray<index_set,M>& ie1_2_dis_ix2, // K1 -> list of ix2
disarray<point_basic<T>,M>& hat_x2)
{
typedef typename field_basic<T,M>::size_type size_type;
// -----------------------------------------------------------------------------
// 1) in order to call only one time per K element the dis_inod() function,
// we group the x2(i) per K element of omega1
// -----------------------------------------------------------------------------
size_type first_dis_ix2 = x2.ownership().first_index();
distributor ie1_ownership = omega1.geo_element_ownership (omega1.map_dimension());
size_type first_dis_ie1 = ie1_ownership.first_index();
std::set<size_type> ext_ie1_set;
index_set empty_set;
ie1_2_dis_ix2.resize (ie1_ownership, empty_set);
for (size_type ix2 = 0, nx2 = ix2_2_dis_ie1.size(); ix2 < nx2; ix2++) {
size_type dis_ie1 = ix2_2_dis_ie1 [ix2];
size_type dis_ix2 = first_dis_ix2 + ix2;
index_set dis_ix2_set; dis_ix2_set.insert (dis_ix2);
check_macro (dis_ie1 != std::numeric_limits<size_type>::max(), "node("<<ix2<<")="<<x2[ix2]
<<" cannot be localized in "<<omega1.name() << " (HINT: curved geometry ?)");
ie1_2_dis_ix2.dis_entry (dis_ie1) += dis_ix2_set;
}
ie1_2_dis_ix2.dis_entry_assembly();
// -----------------------------------------------------------------------------
// 2) loop on K and list external x(i)
// -----------------------------------------------------------------------------
// TODO: it would be more efficient to send (dis_ix,x) together : less mpi calls
// => how to do an assembly() with a disarray<pair<ix,x> > ?
// ... it does not support it yet
// also: would use more memory: massive coordinates copy...
communicator comm = ie1_ownership.comm();
distributor ix2_ownership = x2.ownership();
index_set ext_dis_ix2_set;
if (is_distributed<M>::value && comm.size() > 1) {
for (size_type ie1 = 0, ne1 = ie1_2_dis_ix2.size(); ie1 < ne1; ie1++) {
const index_set& dis_ix2_set = ie1_2_dis_ix2 [ie1];
for (typename index_set::const_iterator iter = dis_ix2_set.begin(), last = dis_ix2_set.end(); iter != last; ++iter) {
size_type dis_ix2 = *iter;
if (ix2_ownership.is_owned (dis_ix2)) continue;
ext_dis_ix2_set.insert (dis_ix2);
}
}
// TODO: when x=xdofs lives in space, then it modifies (enlarges)
// definitively the set of externals indexes & values
// => massive permanent storage
// otherwise: manage the scatter in a separate map<ext_ix,x>
// map<ix,x> ext_x_map;
// x.append_dis_indexes (ext_dis_ix_set, ext_x_map);
// bug: see field_reinterpolate2_tst.cc
// x.set_dis_indexes (ext_dis_ix_set);
x2.append_dis_indexes (ext_dis_ix2_set);
}
// -----------------------------------------------------------------------------
// 3) loop on K and compute hat_x(dis_ix)
// -----------------------------------------------------------------------------
T infty = std::numeric_limits<T>::max();
hat_x2.resize (ix2_ownership, point_basic<T>(infty,infty,infty));
std::vector<size_type> dis_inod1;
for (size_type ie1 = 0, ne1 = ie1_2_dis_ix2.size(); ie1 < ne1; ie1++) {
const index_set& dis_ix2_set = ie1_2_dis_ix2 [ie1];
if (dis_ix2_set.size() == 0) continue;
const geo_element& K1 = omega1 [ie1];
omega1.dis_inod (K1, dis_inod1);
for (typename index_set::const_iterator iter = dis_ix2_set.begin(), last = dis_ix2_set.end(); iter != last; ++iter) {
size_type dis_ix2 = *iter;
const point_basic<T>& xi2 = x2.dis_at (dis_ix2);
hat_x2.dis_entry (dis_ix2) = inverse_piola_transformation (omega1, K1, dis_inod1, xi2);
}
}
hat_x2.dis_entry_assembly();
hat_x2.set_dis_indexes (ext_dis_ix2_set);
}
template<class T, class M, class Value>
void
interpolate_pass2_valued (
const field_basic<T,M>& u1h,
const disarray<point_basic<T>,M>& x2,
const disarray<index_set,M>& ie1_2_dis_ix2, // K1 -> list of ix2
const disarray<point_basic<T>,M>& hat_x2, // ix2 -> hat_x2
disarray<Value,M>& ux2)
{
typedef typename field_basic<T,M>::size_type size_type;
u1h.dis_dof_update();
T infty = std::numeric_limits<T>::max();
distributor x2_ownership = x2.ownership();
ux2.resize (x2.ownership(), Value(infty));
// -----------------------------------------------------------------------------
// on locally managed K1, evaluate at hat_x2 that could be external
// -----------------------------------------------------------------------------
const geo_basic<T,M>& omega1 = u1h.get_geo();
const space_basic<T,M>& Xh1 = u1h.get_space();
const basis_basic<T>& b1 = u1h.get_space().get_basis();
check_macro (! b1.get_piola_fem().transform_need_piola(), "unsupported piola-based basis");
std::vector<size_type> dis_idof1;
std::vector<T> dof1;
Eigen::Matrix<Value,Eigen::Dynamic,1> b1_value;
for (size_type ie1 = 0, ne1 = ie1_2_dis_ix2.size(); ie1 < ne1; ie1++) {
const index_set& dis_ix2_set = ie1_2_dis_ix2 [ie1];
if (dis_ix2_set.size() == 0) continue;
// extract dis_idof[] & dof[] one time for all on K1
const geo_element& K1 = omega1 [ie1];
Xh1.dis_idof (K1, dis_idof1);
size_type loc_ndof1 = dis_idof1.size();
dof1.resize(loc_ndof1);
for (size_type loc_idof1 = 0; loc_idof1 < loc_ndof1; loc_idof1++) {
dof1 [loc_idof1] = u1h.dis_dof (dis_idof1 [loc_idof1]);
}
b1_value.resize (loc_ndof1);
// loop on all hat_x in K: eval basis at hat_x and combine with dof[]
for (typename index_set::const_iterator iter = dis_ix2_set.begin(), last = dis_ix2_set.end(); iter != last; ++iter) {
size_type dis_ix2 = *iter;
size_type iproc = x2_ownership.find_owner (dis_ix2);
size_type nx2 = x2_ownership.size (iproc);
size_type first_dis_ix2 = x2_ownership.first_index(iproc);
size_type ix2 = dis_ix2 - first_dis_ix2;
const point_basic<T>& hat_xi = hat_x2.dis_at (dis_ix2);
b1.evaluate (K1, hat_xi, b1_value); // TODO: RTk and Hdiv: apply also Piola transform at hat_xi
Value value = Value();
for (size_type loc_idof1 = 0; loc_idof1 < loc_ndof1; loc_idof1++) {
value += dof1 [loc_idof1] * b1_value[loc_idof1]; // sum_i w_coef(i)*hat_phi(hat_x)
}
ux2.dis_entry (dis_ix2) = value;
}
}
ux2.dis_entry_assembly();
}
template<class T, class M, class Value>
void
interpolate_pass3_dof (
const disarray<Value,M>& ux2,
field_basic<T,M>& u2h)
{
typedef typename field_basic<T,M>::size_type size_type;
std::set<size_type> ext_dis_inod_set;
u2h.get_space().get_xdofs().get_dis_indexes (ext_dis_inod_set);
ux2.set_dis_indexes (ext_dis_inod_set);
const space_basic<T,M>& V2h = u2h.get_space();
const geo_basic<T,M>& omega2 = u2h.get_geo();
const basis_basic<T>& b2 = V2h.get_basis();
Eigen::Matrix<Value,Eigen::Dynamic,1> u_xnod2;
Eigen::Matrix<T,Eigen::Dynamic,1> udof2;
std::vector<size_type> dis_inod2, dis_idof2;
for (size_type ie2 = 0, ne2 = omega2.size(); ie2 < ne2; ++ie2) {
const geo_element& K2 = omega2 [ie2];
space_numbering::dis_inod (V2h.get_basis(), omega2.sizes(), K2, dis_inod2);
u_xnod2.resize (dis_inod2.size());
for (size_type loc_inod2 = 0, loc_nnod2 = dis_inod2.size(); loc_inod2 < loc_nnod2; ++loc_inod2) {
check_macro (dis_inod2[loc_inod2] < ux2.ownership().dis_size(), "invalid size");
u_xnod2 [loc_inod2] = ux2.dis_at (dis_inod2[loc_inod2]);
}
b2.compute_dofs (K2, u_xnod2, udof2);
space_numbering::dis_idof (V2h.get_basis(), omega2.sizes(), K2, dis_idof2);
check_macro (dis_idof2.size() == size_type(udof2.size()), "invalid sizes");
for (size_type loc_idof2 = 0, loc_ndof2 = dis_idof2.size(); loc_idof2 < loc_ndof2; ++loc_idof2) {
u2h.dis_dof_entry (dis_idof2[loc_idof2]) = udof2 [loc_idof2];
}
}
u2h.dis_dof_update();
}
template<class T, class M, class Value>
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)
{
const geo_basic<T,M>& omega1 = u1h.get_geo();
disarray<index_set,M> ie2dis_ix;
disarray<point_basic<T>,M> hat_x;
interpolate_pass1_symbolic (omega1, x2, ix2dis_ie, ie2dis_ix, hat_x);
interpolate_pass2_valued (u1h, x2, ie2dis_ix, hat_x, ux2);
}
} // namespace details
// ============================================================================
// interpolate function:
// re-interpolate one field on another mesh and space:
// field u2h = interpolate (V2h, u1h);
// ============================================================================
template<class T, class M>
field_basic<T,M>
interpolate (const space_basic<T,M>& V2h, const field_basic<T,M>& u1h)
{
u1h.dis_dof_update();
u1h.get_geo().boundary(); // force "boundary" domain creation, for geo locator
typedef typename field_basic<T,M>::size_type size_type;
if (u1h.get_space() == V2h) {
size_type have_same_dofs = (u1h.u().size() == V2h.iu_ownership().size());
#ifdef _RHEOLEF_HAVE_MPI
if (is_distributed<M>::value) {
have_same_dofs = mpi::all_reduce (V2h.ownership().comm(), have_same_dofs, mpi::minimum<size_type>());
}
#endif // _RHEOLEF_HAVE_MPI
if (have_same_dofs) {
// spaces are exactly the same: no need to re-interpolate or copy
return u1h;
}
// spaces differs only by blocked/unblocked dofs: need to copy
field_basic<T,M> u2h (V2h);
for (size_type idof = 0, ndof = V2h.ndof(); idof < ndof; ++idof) {
u2h.dof(idof) = u1h.dof(idof);
}
u2h.dis_dof_update();
return u2h;
}
if (u1h.get_geo().get_background_geo() == V2h.get_geo().get_background_geo()) {
// meshes are compatible:
// => buid a wrapper expression and call back interpolate with the
// general nonlinear specialized version
return interpolate (V2h, details::field_expr_v2_nonlinear_terminal_field<T,M,details::differentiate_option::none>(u1h));
}
// -----------------------------------------------------------------------------
// between different meshes:
// 0) create the xdofs2 set of nodes
// -----------------------------------------------------------------------------
const disarray<point_basic<T>,M>& xdof2 = V2h.get_xdofs();
// -----------------------------------------------------------------------------
// 1) locate each xdof of V2 in omega1
// -----------------------------------------------------------------------------
disarray<point_basic<T>,M> xdof2_nearest (xdof2.ownership());
const geo_basic<T,M>& omega1 = u1h.get_geo();
disarray<size_type,M> dis_ie1_tab;
omega1.nearest (xdof2, xdof2_nearest, dis_ie1_tab);
// -----------------------------------------------------------------------------
// 2) interpolate uh1 at xdof2_nearest: get the values in a disarray
// -----------------------------------------------------------------------------
field_basic<T,M> u2h (V2h);
switch (V2h.valued_tag()) {
case space_constant::scalar: {
disarray<T,M> u2h_dof;
details::interpolate_on_a_different_mesh (u1h, xdof2_nearest, dis_ie1_tab, u2h_dof);
details::interpolate_pass3_dof (u2h_dof, u2h);
return u2h;
}
case space_constant::vector: {
disarray<point_basic<T>,M> u2h_dof;
details::interpolate_on_a_different_mesh (u1h, xdof2_nearest, dis_ie1_tab, u2h_dof);
details::interpolate_pass3_dof (u2h_dof, u2h);
return u2h;
}
case space_constant::tensor:
case space_constant::unsymmetric_tensor: {
disarray<tensor_basic<T>,M> u2h_dof;
details::interpolate_on_a_different_mesh (u1h, xdof2_nearest, dis_ie1_tab, u2h_dof);
details::interpolate_pass3_dof (u2h_dof, u2h);
return u2h;
}
default:
error_macro("interpolate between different meshes: unexpected "<<V2h.valued()<<"-valued field");
}
return field_basic<T,M>();
}
// ----------------------------------------------------------------------------
// instanciation in library
// ----------------------------------------------------------------------------
#define _RHEOLEF_instanciation_value(T,M,Value) \
template void details::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_x, \
disarray<Value,M>& ux); \
#define _RHEOLEF_instanciation(T,M) \
_RHEOLEF_instanciation_value(T,M,T) \
_RHEOLEF_instanciation_value(T,M,point_basic<T>) \
template \
void \
details::interpolate_pass1_symbolic ( \
const geo_basic<T,M>&, \
const disarray<point_basic<T>,M>&, \
const disarray<geo_element::size_type,M>&, \
disarray<index_set,M>&, \
disarray<point_basic<T>,M>&); \
template \
field_basic<T,M> \
interpolate (const space_basic<T,M>&, const field_basic<T,M>&);
_RHEOLEF_instanciation(Float,sequential)
#ifdef _RHEOLEF_HAVE_MPI
_RHEOLEF_instanciation(Float,distributed)
#endif // _RHEOLEF_HAVE_MPI
} // namespace rheolef
_RHEOLEF_instanciation
#define _RHEOLEF_instanciation(T, M, A)
Definition asr.cc:223
Float
see the Float page for the full documentation
communicator
see the communicator page for the full documentation
rheolef::basis_basic
Definition basis.h:566
rheolef::basis_basic::get_piola_fem
const piola_fem< T > & get_piola_fem() const
Definition basis.h:834
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:964
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:919
rheolef::details::field_expr_v2_nonlinear_terminal_field
Definition field_expr_terminal.h:807
rheolef::disarray
see the disarray page for the full documentation
Definition disarray.h:497
rheolef::distributor
see the distributor page for the full documentation
Definition distributor.h:69
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::distributor::size
size_type size(size_type iproc) const
Definition distributor.h:170
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:158
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:220
rheolef::distributor::comm
const communicator_type & comm() const
Definition distributor.h:152
rheolef::field_basic
Definition field.h:219
rheolef::field_basic::dis_dof
const T & dis_dof(size_type dis_idof) const
Definition field.cc:377
rheolef::field_basic::dof
T & dof(size_type idof)
Definition field.h:738
rheolef::field_basic::dis_dof_update
void dis_dof_update(const SetOp &=SetOp()) const
Definition field.h:763
rheolef::field_basic::dis_dof_entry
dis_reference dis_dof_entry(size_type dis_idof)
Definition field.cc:398
rheolef::field_basic::get_geo
const geo_type & get_geo() const
Definition field.h:271
rheolef::field_basic::get_space
const space_type & get_space() const
Definition field.h:270
rheolef::field_basic::u
const vec< T, M > & u() const
Definition field.h:282
rheolef::field_basic::size_type
std::size_t size_type
Definition field.h:225
rheolef::geo_basic
generic mesh with rerefence counting
Definition geo.h:1089
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::index_set::insert
void insert(size_type dis_i)
Definition index_set_header.icc:159
rheolef::point_basic
Definition point.h:86
rheolef::space_basic
the finite element space
Definition space.h:382
error_macro
#define error_macro(message)
Definition dis_macros.h:49
T
Expr1::float_type T
Definition field_expr.h:230
Xh1
space_basic< T, M > Xh1
Definition field_expr.h:232
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
Definition compiler_eigen.h:39
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::details::interpolate_pass3_dof
void interpolate_pass3_dof(const disarray< Value, M > &ux2, field_basic< T, M > &u2h)
Definition interpolate.cc:176
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::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::space_constant::tensor
@ tensor
Definition space_constant.h:136
rheolef::space_constant::scalar
@ scalar
Definition space_constant.h:134
rheolef::space_constant::unsymmetric_tensor
@ unsymmetric_tensor
Definition space_constant.h:137
rheolef::space_constant::vector
@ vector
Definition space_constant.h:135
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::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
rheolef
This file is part of Rheolef.
Definition compiler_eigen.h:39
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::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:473
rheolef::is_distributed::value
static const bool value
Definition communicator.h:71
M
Expr1::memory_type M
Definition vec_expr_v2.h:416