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()
{
// omega_in is compressed by space_scalar_constitution() 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()
{
// omega_in is compressed by space_scalar_constitution() 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()
{
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()
{
_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_indirect_array (
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_indirect_array (Wh, bgd_gamma);
}
template <class T, class M>
disarray<typename space_base_rep<T,M>::size_type, M>
space_base_rep<T,M>::build_indirect_array (
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();
trace_macro("Vh="<<Vh.name()<<", Wh="<<Wh.name()<<", dom_idof2bgd_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_idof2bgd_idof (dom_ownership, std::numeric_limits<size_type>::max());
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_idof2bgd_idof.dis_entry (dom_dis_idof) = bgd_dis_idof;
trace_macro("ige="<<ige<<": dom_idof2bgd_idof["<<dom_dis_idof<<"] = "<<bgd_dis_idof);
}
}
dom_idof2bgd_idof.dis_entry_assembly();
// move to local numbering:
for (size_type dom_idof = 0, dom_ndof = dom_idof2bgd_idof.size(); dom_idof < dom_ndof; dom_idof++) {
size_type bgd_dis_idof = dom_idof2bgd_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_idof2bgd_idof [dom_idof] = bgd_idof;
}
trace_macro("Vh="<<Vh.name()<<", Wh="<<Wh.name()<<", dom_idof2bgd_idof.size=Wh.size="<<dom_ownership.size()<<" done");
return dom_idof2bgd_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
// ----------------------------------------------------------------------------
// 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:425
Float
see the Float page for the full documentation
rheolef::basis_basic
Definition: basis.h:532
rheolef::disarray< size_type, M >
rheolef::distributor
see the distributor page for the full documentation
Definition: distributor.h:62
rheolef::distributor::last_index
size_type last_index(size_type iproc) const
Definition: distributor.h:157
rheolef::distributor::size
size_type size(size_type iproc) const
Definition: distributor.h:163
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:151
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:213
rheolef::geo_basic< T, M >
rheolef::geo_element
see the geo_element page for the full documentation
Definition: geo_element.h:102
rheolef::geo_element::dis_ie
size_type dis_ie() const
Definition: geo_element.h:163
rheolef::geo_element::name
char name() const
Definition: geo_element.h:169
rheolef::index_set
Definition: index_set_header.icc:58
rheolef::index_set::insert
void insert(size_type dis_i)
Definition: index_set_header.icc:158
rheolef::integrate_option
see the integrate_option page for the full documentation
Definition: integrate_option.h:125
rheolef::piola_on_pointset
Definition: piola_on_pointset.h:142
rheolef::piola_on_pointset::initialize
void initialize(const basis_basic< T > &piola_basis, const quadrature< T > &quad, const integrate_option &iopt)
Definition: piola_on_pointset.h:190
rheolef::piola_on_pointset::get_piola
const Eigen::Matrix< piola< T >, Eigen::Dynamic, 1 > & get_piola(const geo_basic< T, M > &omega, const geo_element &K) const
Definition: piola_on_pointset.h:239
rheolef::space_base_rep
Definition: space.h:159
rheolef::space_base_rep::_normal
disarray< point_basic< T >, M > _normal
Definition: space.h:259
rheolef::space_base_rep::dis_idof
void dis_idof(const geo_element &K, std::vector< size_type > &dis_idof) const
Definition: space.cc:200
rheolef::space_base_rep::base_freeze_body
void base_freeze_body() const
Definition: space.cc:157
rheolef::space_base_rep::name
std::string name() const
Definition: space.cc:210
rheolef::space_base_rep::build_indirect_array
disarray< size_type, M > build_indirect_array(const space_base_rep< T, M > &Wh, const std::string &dom_name) const
Definition: space.cc:229
rheolef::space_base_rep::_constit
space_constitution< T, M > _constit
Definition: space.h:254
rheolef::space_base_rep::space_base_rep
space_base_rep()
rheolef::space_base_rep::get_geo
const geo_basic< T, M > & get_geo() const
Definition: space.h:189
rheolef::space_base_rep::_has_nt_basis
disarray< int, M > _has_nt_basis
Definition: space.h:258
rheolef::space_base_rep::ownership
const distributor & ownership() const
Definition: space.h:183
rheolef::space_base_rep::size_type
space_pair_type::size_type size_type
Definition: space.h:164
rheolef::space_base_rep::freeze_guard
void freeze_guard() const
Definition: space.h:243
rheolef::space_base_rep::_idof2blk_dis_iub
disarray< space_pair_type, M > _idof2blk_dis_iub
Definition: space.h:257
rheolef::space_base_rep::init_xdof
void init_xdof()
Definition: space.cc:114
rheolef::space_constitution< T, M >
rheolef::space_constitution::ownership
const distributor & ownership() const
Definition: space_constitution.h:410
rheolef::space_mult_list
Definition: space_mult.h:54
trace_macro
#define trace_macro(message)
Definition: dis_macros.h:111
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
expr
Definition: mkgeo_ball.sh:361
mkgeo_ball.b
int b
Definition: mkgeo_ball.sh:152
mkgeo_ball.map_dim
map_dim
Definition: mkgeo_ball.sh:337
mkgeo_contraction.name
string 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::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::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
This file is part of Rheolef.
Definition: compiler_eigen.h:37
rheolef::_RHEOLEF_space_real
_RHEOLEF_space_real(sequential) _RHEOLEF_space_real(distributed) _RHEOLEF_instanciation(Float
rheolef::_RHEOLEF_instanciation
_RHEOLEF_instanciation(Float, sequential, std::allocator< Float >) _RHEOLEF_instanciation(Float
rheolef::distributed
Definition: distributed.h:88
rheolef::is_distributed
Definition: distributed.h:36
rheolef::is_distributed::value
static const bool value
Definition: distributed.h:37
rheolef::piola
Definition: piola.h:67
rheolef::piola::F
point_basic< T > F
Definition: piola.h:79
rheolef::sequential
Definition: distributed.h:28
M
Expr1::memory_type M
Definition: vec_expr_v2.h:416