Loading...
Searching...
No Matches
space.cc
Go to the documentation of this file.
#include "rheolef/space.h"
#include "rheolef/space_numbering.h"
#include "rheolef/space_mult.h"
#include "rheolef/piola_util.h"
#include "rheolef/piola_on_pointset.h"
namespace rheolef {
// ======================================================================
// space_base_rep
// ======================================================================
template <class T, class M>
space_base_rep<T,M>::space_base_rep (
const geo_basic<T,M>& omega_in,
std::string approx,
std::string valued)
: _constit(omega_in,
(approx == "" || valued == "scalar") ? approx : valued+"("+approx+")"),
_xdof(),
_have_freezed(false),
_idof2blk_dis_iub(),
_has_nt_basis(),
_normal(),
_iu_ownership(),
_ib_ownership(),
_parent_subgeo_owner_initialized(false),
_parent_subgeo_owner_reattribution_required(false),
_parent_subgeo_owner()
{
// omega_in is compressed by constitution_terminal() allocator
// when it is a domain: then it becomes a geo_domain, with compressed numbering
// so, omega_in can be different from omega:
if (approx == "") return; // empty element => default cstor
_idof2blk_dis_iub.resize (_constit.ownership());
_has_nt_basis.resize (_constit.ownership()); // TODO: _constit.scalar_ownership()
_normal.resize (_constit.ownership());
init_xdof();
}
template <class T, class M>
space_base_rep<T,M>::space_base_rep (
const geo_basic<T,M>& omega_in,
const basis_basic<T>& b)
: _constit(omega_in, b.name()),
_xdof(),
_have_freezed(false),
_idof2blk_dis_iub(),
_has_nt_basis(),
_normal(),
_iu_ownership(),
_ib_ownership(),
_parent_subgeo_owner_initialized(false),
_parent_subgeo_owner_reattribution_required(false),
_parent_subgeo_owner()
{
// omega_in is compressed by constitution_terminal() allocator
// when it is a domain: then it becomes a geo_domain, with compressed numbering
// so, omega_in can be different from omega:
if (! b.is_initialized()) return; // empty element => default cstor
_idof2blk_dis_iub.resize (_constit.ownership());
_has_nt_basis.resize (_constit.ownership());
_normal.resize (_constit.ownership());
init_xdof();
}
template <class T, class M>
space_base_rep<T,M>::space_base_rep (const space_constitution<T,M>& constit)
: _constit(constit),
_xdof(),
_have_freezed(false),
_idof2blk_dis_iub(),
_has_nt_basis(),
_normal(),
_iu_ownership(),
_ib_ownership(),
_parent_subgeo_owner_initialized(false),
_parent_subgeo_owner_reattribution_required(false),
_parent_subgeo_owner()
{
distributor idof_ownership = _constit.ownership();
_idof2blk_dis_iub.resize (idof_ownership);
_has_nt_basis.resize (idof_ownership);
_normal.resize (idof_ownership);
init_xdof();
}
template <class T, class M>
space_base_rep<T,M>::space_base_rep (const space_mult_list<T,M>& expr)
: _constit(expr),
_xdof(),
_have_freezed(false),
_idof2blk_dis_iub(),
_has_nt_basis(),
_normal(),
_iu_ownership(),
_ib_ownership(),
_parent_subgeo_owner_initialized(false),
_parent_subgeo_owner_reattribution_required(false),
_parent_subgeo_owner()
{
_idof2blk_dis_iub.resize (_constit.ownership());
_has_nt_basis.resize (_constit.ownership());
_normal.resize (_constit.ownership());
init_xdof();
}
template <class T, class M>
void
space_base_rep<T,M>::init_xdof()
{
if (_constit.is_hierarchical()) {
trace_macro ("init_xdof: hierarchical space: xdof undefined");
return;
}
const geo_basic<T,M>& omega = get_geo();
const basis_basic<T>& b = get_basis();
size_type dis_nnod = space_numbering::dis_nnod (b, omega.sizes(), omega.map_dimension());
size_type nnod = space_numbering:: nnod (b, omega.sizes(), omega.map_dimension());
communicator comm = ownership().comm();
distributor nod_ownership (dis_nnod, comm, nnod);
_xdof.resize (nod_ownership);
piola_on_pointset<T> pops;
pops.initialize (omega.get_piola_basis(), b, integrate_option());
std::vector<size_type> dis_inod_tab;
for (size_type ie = 0, ne = omega.size(); ie < ne; ie++) {
const geo_element& K = omega [ie];
const Eigen::Matrix<piola<T>,Eigen::Dynamic,1>& piola = pops.get_piola (omega, K);
space_numbering::dis_inod (b, omega.sizes(), K, dis_inod_tab);
for (size_type loc_inod = 0, loc_nnod = piola.size(); loc_inod < loc_nnod; ++loc_inod) {
_xdof.dis_entry (dis_inod_tab[loc_inod]) = piola [loc_inod].F;
}
}
_xdof.dis_entry_assembly();
// add external nodes on current element:
if (is_distributed<M>::value && comm.size() > 1) {
index_set dis_inod_set;
for (size_type ie = 0, ne = omega.size(); ie < ne; ++ie) {
const geo_element& K = omega [ie];
space_numbering::dis_inod (b, omega.sizes(), K, dis_inod_tab);
for (size_type loc_inod = 0, loc_nnod = dis_inod_tab.size(); loc_inod < loc_nnod; ++loc_inod) {
size_type dis_inod = dis_inod_tab[loc_inod];
if (nod_ownership.is_owned(dis_inod)) continue;
check_macro (dis_inod < dis_nnod, "invalid dis_inod="<<dis_inod<<" out of range [0:"<<dis_nnod<<"[");
dis_inod_set.insert (dis_inod);
}
}
_xdof.set_dis_indexes (dis_inod_set);
}
}
template <class T, class M>
void
space_base_rep<T,M>::base_freeze_body () const
{
_constit.neighbour_guard(); // geo could uses S.master(i), a neighbour connnectivity
// -----------------------------------------------------------------------
// 1) loop on domains: mark blocked dofs
// -----------------------------------------------------------------------
disarray<size_type,M> blocked_flag = _constit.build_blocked_flag();
// copy the blocked_flag into the dis_iub disarray, as the "blocked" bit:
for (size_type idof = 0, ndof = blocked_flag.size(); idof < ndof; idof++) {
_idof2blk_dis_iub [idof].set_blocked (blocked_flag[idof]);
}
// -----------------------------------------------------------------------
// 2) init numbering
// -----------------------------------------------------------------------
size_type n_unknown = 0;
size_type n_blocked = 0;
for (size_type idof = 0, ndof = _idof2blk_dis_iub.size(); idof < ndof; idof++) {
bool blk = _idof2blk_dis_iub [idof].is_blocked();
if (! blk) {
_idof2blk_dis_iub[idof].set_dis_iub (n_unknown);
n_unknown++;
} else {
_idof2blk_dis_iub[idof].set_dis_iub (n_blocked);
n_blocked++;
}
}
size_type dis_n_unknown = n_unknown;
size_type dis_n_blocked = n_blocked;
#ifdef _RHEOLEF_HAVE_MPI
if (is_distributed<M>::value) {
dis_n_unknown = mpi::all_reduce (comm(), dis_n_unknown, std::plus<T>());
dis_n_blocked = mpi::all_reduce (comm(), dis_n_blocked, std::plus<T>());
}
#endif // // _RHEOLEF_HAVE_MPI
_iu_ownership = distributor (dis_n_unknown, comm(), n_unknown);
_ib_ownership = distributor (dis_n_blocked, comm(), n_blocked);
}
// ----------------------------------------------------------------------------
// still used by geo_subdivide.cc and internally in space.cc
// ----------------------------------------------------------------------------
template <class T, class M>
void
space_base_rep<T,M>::dis_idof (const geo_element& K, std::vector<size_type>& dis_idof) const
{
freeze_guard();
_constit.assembly_dis_idof (get_geo(), K, dis_idof);
}
// ----------------------------------------------------------------------------
// name : e.g. "P1(square)", for field_expr<Expr> checks
// ----------------------------------------------------------------------------
template <class T, class M>
std::string
space_base_rep<T,M>::name() const
{
return _constit.name();
}
// ----------------------------------------------------------------------------
// u["left"]
// => build here the requiresd temporary indirect disarray
// ----------------------------------------------------------------------------
template <class T, class M>
disarray<typename space_base_rep<T,M>::size_type, M>
space_base_rep<T,M>::build_dom_dis_idof2bgd_dis_idof (
const space_base_rep<T,M>& Wh, // Wh = space on domain
const std::string& dom_name // redundant: contained in Wh.get_geo()
) const
{
const geo_basic<T,M>& bgd_gamma = get_geo()[dom_name];
return build_dom_dis_idof2bgd_dis_idof (Wh, bgd_gamma);
}
template <class T, class M>
disarray<typename space_base_rep<T,M>::size_type, M>
space_base_rep<T,M>::build_dom_dis_idof2bgd_dis_idof (
const space_base_rep<T,M>& Wh, // Wh = space on domain
const geo_basic<T,M>& bgd_gamma2
) const
{
// TODO: move it to the call ?
const geo_basic<T,M>& bgd_gamma = bgd_gamma2.get_background_domain();
// TODO: verifier que le domaine bgd_gamma est compatible:
// => il doit y avoir un meme maillage de base a bgd_gamma & Wh.get_geo
const space_base_rep<T,M>& Vh = *this; // Vh = space on background mesh
Vh.freeze_guard();
Wh.freeze_guard();
const geo_basic<T,M>& dom_gamma = Wh.get_geo();
size_type map_dim = dom_gamma.map_dimension();
check_macro (dom_gamma.size() == bgd_gamma.size(), "incompatible domains");
distributor dom_ownership = Wh.ownership();
distributor bgd_ownership = Vh.ownership();
size_type my_proc = dom_ownership.comm().rank();
trace_macro("Vh="<<Vh.name()<<", Wh="<<Wh.name()<<", dom_dis_idof2bgd_dis_idof.size=Wh.size="<<dom_ownership.size()<<"...");
size_type first_dom_dis_idof = dom_ownership.first_index();
size_type first_bgd_dis_idof = bgd_ownership.first_index();
std::vector<size_type> dom_dis_idofs, bgd_dis_idofs;
disarray<size_type, M> dom_dis_idof2bgd_dis_idof (dom_ownership, std::numeric_limits<size_type>::max());
std::set<size_type> ext_dom_dis_idof;
for (size_type ige = 0, nge = dom_gamma.size(); ige < nge; ige++) {
const geo_element& dom_S = dom_gamma[ige];
const geo_element& bgd_S = bgd_gamma[ige];
Wh.dis_idof (dom_S, dom_dis_idofs);
Vh.dis_idof (bgd_S, bgd_dis_idofs);
trace_macro("ige="<<ige<<": dom_S="<<dom_S.name()<<dom_S.dis_ie()<<", bgd_S="<<bgd_S.name()<<bgd_S.dis_ie()
<< ": dom_dis_idofs.size="<<dom_dis_idofs.size()
<< ", bgd_dis_idofs.size="<<bgd_dis_idofs.size());
for (size_type loc_idof = 0, loc_ndof = dom_dis_idofs.size(); loc_idof < loc_ndof; loc_idof++) {
size_type dom_dis_idof = dom_dis_idofs [loc_idof];
size_type bgd_dis_idof = bgd_dis_idofs [loc_idof];
dom_dis_idof2bgd_dis_idof.dis_entry (dom_dis_idof) = bgd_dis_idof;
trace_macro("ige="<<ige<<": dom_dis_idof2bgd_dis_idof["<<dom_dis_idof<<"] = "<<bgd_dis_idof);
if (! dom_ownership.is_owned (dom_dis_idof, my_proc)) {
ext_dom_dis_idof.insert (dom_dis_idof);
}
}
}
trace_macro ("ext_dom_dis_idof.size="<<ext_dom_dis_idof.size());
dom_dis_idof2bgd_dis_idof.dis_entry_assembly();
dom_dis_idof2bgd_dis_idof.set_dis_indexes (ext_dom_dis_idof); // doit etre en deuxieme...?
#ifdef TO_CLEAN
// move to local numbering:
for (size_type dom_idof = 0, dom_ndof = dom_dis_idof2bgd_dis_idof.size(); dom_idof < dom_ndof; dom_idof++) {
size_type bgd_dis_idof = dom_dis_idof2bgd_dis_idof [dom_idof];
size_type dom_dis_idof = dom_idof + first_dom_dis_idof;
check_macro (bgd_ownership.is_owned (bgd_dis_idof), "bgd_dis_idof="<<bgd_dis_idof<<" is out of range ["<<first_bgd_dis_idof<<":"<< bgd_ownership.last_index()<<"[");
size_type bgd_idof = bgd_dis_idof - first_bgd_dis_idof;
dom_dis_idof2bgd_dis_idof [dom_idof] = bgd_idof;
}
trace_macro("Vh="<<Vh.name()<<", Wh="<<Wh.name()<<", dom_dis_idof2bgd_dis_idof.size=Wh.size="<<dom_ownership.size()<<" done");
#endif // TO_CLEAN
return dom_dis_idof2bgd_dis_idof;
}
#define _RHEOLEF_space_real(M) \
template <class T> \
space_basic<T,M> \
space_basic<T,M>::real() \
{ \
return space_basic<T,M> (geo_basic<T,M>(details::zero_dimension()), "P1"); \
}
_RHEOLEF_space_real(sequential)
#ifdef _RHEOLEF_HAVE_MPI
_RHEOLEF_space_real(distributed)
#endif // _RHEOLEF_HAVE_MPI
#undef _RHEOLEF_space_real
// ======================================================================
// parent_subgeo_owner
// ======================================================================
/* IMPLEMENTATION NOTES:
Consider the following situation:
P0 S0 top
+-----+------+--
|\ K0 |
S1 | \ |
| K1 \|
P1 +-----+
|
left |
+
vertex P0 is shared by edges S0 and S1 in triangle K0
edge S0 belongs to domain "top"
edge S1 belongs to domain "left"
assume that the mesh partition attributes
K0 on partition iproc=0
K1 on partition iproc=1
then, sides and vertex are also attributed to a proc
via a recursive descent on subgeos:
iproc=0 K0, S0, P0
iproc=1 K1, S1, P1
when explorating S1, its vertex P0 is already attributed to iproc=0
Then, consider:
space Wh (omega["left"], "P1");
The element edge S1 contains an **orphean** dof associated to vertex P0 :
it do not belong to S1 owner iproc=1 and is not referenced by any others sides
of "left" that belongs to its owner iproc=0
So, in order to interpolate in a lazy element-by-element loop, we have to
re-attribute this dof to S0 owner, iproc=1.
notes:
- this idof iproc owner reattribution required only when
* nproc > 1
* continuous elements
* omega is a geo_domain, i.e. some parent elements are missing
otherwise could return usual dis_idof owner
- all dofs are associated to a subgeo
- when continuous approx, this subgeo corresponds to a mesh subgeo
- when discontinuous, this subgeo de-multiplicated from mesh:
* verticies 1D, faces 3D or edges 2D : doubled when internals
* edges 3D and verticies 2D and 3D : de-multiplied
in all cases, dofs conceptually conserve the memory of the parent subgeo
*/
template <class T, class M>
typename space_base_rep<T,M>::size_type
space_base_rep<T,M>::get_parent_subgeo_owner (size_type dis_idof) const
{
parent_subgeo_owner_check();
if (! _parent_subgeo_owner_reattribution_required) {
return ownership().find_owner (dis_idof);
}
return _parent_subgeo_owner.dis_at (dis_idof);
}
namespace details {
struct pair_dof_parent_owners_min_op {
using t = std::pair<size_t,size_t>; // pair: first=parent_owner & second=dof_owner
t operator() (const t& a, const t& b) {
// if "a" or "b" parent owner match the dof owner: select it
// otherwise, choose the parent owner arbitrarily
// avoid unset=max(size_t) by tacking the min
if (a.first == a.second) return a;
if (b.first == b.second) return b;
size_t parent_owner = std::min (a.first, b.first); // avoid unset=max
return std::make_pair (parent_owner, a.second);
}
};
} // namespace details
template <class T, class M>
bool
space_base_rep<T,M>::parent_subgeo_owner_guard() const
{
if (_parent_subgeo_owner_initialized) {
return _parent_subgeo_owner_reattribution_required;
}
_parent_subgeo_owner_initialized = true;
freeze_guard();
check_macro (! get_constitution().is_hierarchical(), "parent_subgeo_owner_reattribution: invalid product-space");
const space_constitution_terminal<T,M>& terminal_constit = get_constitution().get_terminal();
size_type nproc = ownership().comm().size();
size_type my_proc = ownership().comm().rank();
const geo_basic<T,M>& omega = terminal_constit.get_geo();
_parent_subgeo_owner_reattribution_required
= terminal_constit.get_basis().is_continuous()
&& nproc > 1
&& omega.variant() == geo_abstract_base_rep<T>::geo_domain;
if (! _parent_subgeo_owner_reattribution_required) {
return _parent_subgeo_owner_reattribution_required;
}
// try to reduce comms: when a parent owner is the dof owner, prefer it
// --> see pair_dof_parent_owners_min_op
std::vector<size_type> dis_idx;
size_type unset = std::numeric_limits<size_type>::max();
std::pair<size_type,size_type> unset_pair (unset, unset);
disarray<std::pair<size_type,size_type>> dof_parent_choice (ownership(), unset_pair);
size_type parent_proc = my_proc;
std::set<size_type> ext_dis_idof;
for (size_type ie = 0, ne = omega.size(); ie < ne; ++ie) {
const geo_element& K = omega [ie];
get_constitution().assembly_dis_idof (omega, K, dis_idx);
for (size_type loc_idof = 0, loc_ndof = dis_idx.size(); loc_idof < loc_ndof; loc_idof++) {
size_type dis_idof = dis_idx[loc_idof];
size_type dof_proc = ownership().find_owner (dis_idof);
dof_parent_choice.dis_entry (dis_idof) = std::make_pair (parent_proc, dof_proc);
if (dof_proc != my_proc) {
ext_dis_idof.insert (dis_idof);
}
}
}
dof_parent_choice.dis_entry_assembly (details::pair_dof_parent_owners_min_op());
dof_parent_choice.set_dis_indexes (ext_dis_idof);
_parent_subgeo_owner.resize (ownership());
std::fill (_parent_subgeo_owner.begin(), _parent_subgeo_owner.end(), std::numeric_limits<size_type>::max());
for (size_type ie = 0, ne = omega.size(); ie < ne; ++ie) {
const geo_element& K = omega [ie];
get_constitution().assembly_dis_idof (omega, K, dis_idx);
for (size_type loc_idof = 0, loc_ndof = dis_idx.size(); loc_idof < loc_ndof; loc_idof++) {
size_type dis_idof = dis_idx[loc_idof];
size_type chosen_parent_proc = dof_parent_choice.dis_at(dis_idof).first;
// when chosen_parent_proc == unset : it means that this dof is orphean
// ie there are no element K in the geo_domain containing this dof and such that K belongs to the dof_owner
size_type reattributed_dof_proc = (chosen_parent_proc != unset) ? chosen_parent_proc : my_proc;
_parent_subgeo_owner.dis_entry (dis_idof) = reattributed_dof_proc;
#ifdef TO_CLEAN
size_type dof_proc = ownership().find_owner (dis_idof);
if (dof_proc != my_proc) {
warning_macro ("K="<<K.name()<<K.dis_ie()<<", dis_idof="<<dis_idof<<": dof_proc="<<dof_proc<<" != my_proc="<<my_proc<<": chosen_parent_proc="<<chosen_parent_proc<<", reattributed_dof_proc="<<reattributed_dof_proc);
}
#endif // TO_CLEAN
}
}
// choose for each dis_idof the proc owner that has the min index:
_parent_subgeo_owner.dis_entry_assembly (details::generic_set_min_op());
_parent_subgeo_owner.set_dis_indexes (ext_dis_idof);
return _parent_subgeo_owner_reattribution_required;
}
// ----------------------------------------------------------------------------
// instanciation in library
// ----------------------------------------------------------------------------
#define _RHEOLEF_instanciation(T,M) \
template class space_base_rep<T,M>; \
template space_basic<T,M> space_basic<T,M>::real();
_RHEOLEF_instanciation(Float,sequential)
#ifdef _RHEOLEF_HAVE_MPI
_RHEOLEF_instanciation(Float,distributed)
#endif // _RHEOLEF_HAVE_MPI
#undef _RHEOLEF_instanciation
} // namespace rheolef
size_type
field::size_type size_type
Definition branch.cc:430
Float
see the Float page for the full documentation
communicator
see the communicator page for the full documentation
rheolef::space_base_rep
Definition space.h:159
rheolef::space_base_rep::space_base_rep
space_base_rep()
rheolef::space_base_rep::size_type
space_pair_type::size_type size_type
Definition space.h:164
trace_macro
#define trace_macro(message)
Definition dis_macros.h:111
warning_macro
#define warning_macro(message)
Definition dis_macros.h:53
get_geo
void get_geo(istream &in, my_geo &omega)
Definition field2gmsh_pos.cc:126
check_macro
check_macro(expr1.have_homogeneous_space(Xh1), "dual(expr1,expr2); expr1 should have homogeneous space. HINT: use dual(interpolate(Xh, expr1),expr2)")
mkgeo_ball.expr
str expr
Definition mkgeo_ball.sh:361
mkgeo_ball.b
int b
Definition mkgeo_ball.sh:152
mkgeo_ball.map_dim
int map_dim
Definition mkgeo_ball.sh:337
mkgeo_contraction.name
str name
Definition mkgeo_contraction.sh:133
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::nnod
size_type nnod(const basis_basic< T > &b, const geo_size &gs, size_type map_dim)
Definition space_numbering.cc:54
rheolef::space_numbering::dis_nnod
size_type dis_nnod(const basis_basic< T > &b, const geo_size &gs, size_type map_dim)
Definition space_numbering.cc:67
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::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_instanciation
#define _RHEOLEF_instanciation(T, M)
Definition space.cc:471
_RHEOLEF_space_real
#define _RHEOLEF_space_real(M)
Definition space.cc:310
rheolef::distributed
Definition communicator.h:128
rheolef::sequential
Definition communicator.h:62
M
Expr1::memory_type M
Definition vec_expr_v2.h:416