22#include "rheolef/space.h"
23#include "rheolef/space_numbering.h"
24#include "rheolef/space_mult.h"
25#include "rheolef/piola_util.h"
26#include "rheolef/piola_on_pointset.h"
33template <
class T,
class M>
35 const geo_basic<T,M>& omega_in,
39 (approx ==
"" || valued ==
"scalar") ? approx : valued+
"("+approx+
")"),
47 _parent_subgeo_owner_initialized(false),
48 _parent_subgeo_owner_reattribution_required(false),
49 _parent_subgeo_owner()
54 if (approx ==
"")
return;
55 _idof2blk_dis_iub.resize (_constit.ownership());
56 _has_nt_basis.resize (_constit.ownership());
57 _normal.resize (_constit.ownership());
60template <
class T,
class M>
61space_base_rep<T,M>::space_base_rep (
62 const geo_basic<T,M>& omega_in,
63 const basis_basic<T>& b)
64: _constit(omega_in,
b.
name()),
72 _parent_subgeo_owner_initialized(false),
73 _parent_subgeo_owner_reattribution_required(false),
74 _parent_subgeo_owner()
79 if (!
b.is_initialized())
return;
80 _idof2blk_dis_iub.resize (_constit.ownership());
81 _has_nt_basis.resize (_constit.ownership());
82 _normal.resize (_constit.ownership());
85template <
class T,
class M>
86space_base_rep<T,M>::space_base_rep (
const space_constitution<T,M>& constit)
95 _parent_subgeo_owner_initialized(false),
96 _parent_subgeo_owner_reattribution_required(false),
97 _parent_subgeo_owner()
99 distributor idof_ownership = _constit.ownership();
100 _idof2blk_dis_iub.resize (idof_ownership);
101 _has_nt_basis.resize (idof_ownership);
102 _normal.resize (idof_ownership);
105template <
class T,
class M>
106space_base_rep<T,M>::space_base_rep (
const space_mult_list<T,M>& expr)
109 _have_freezed(false),
115 _parent_subgeo_owner_initialized(false),
116 _parent_subgeo_owner_reattribution_required(false),
117 _parent_subgeo_owner()
119 _idof2blk_dis_iub.resize (_constit.ownership());
120 _has_nt_basis.resize (_constit.ownership());
121 _normal.resize (_constit.ownership());
124template <
class T,
class M>
126space_base_rep<T,M>::init_xdof()
128 if (_constit.is_hierarchical()) {
129 trace_macro (
"init_xdof: hierarchical space: xdof undefined");
132 const geo_basic<T,M>& omega =
get_geo();
133 const basis_basic<T>&
b = get_basis();
134 size_type dis_nnod = space_numbering::dis_nnod (b, omega.sizes(), omega.map_dimension());
135 size_type nnod = space_numbering:: nnod (b, omega.sizes(), omega.map_dimension());
137 distributor nod_ownership (dis_nnod, comm, nnod);
138 _xdof.resize (nod_ownership);
139 piola_on_pointset<T> pops;
140 pops.initialize (omega.get_piola_basis(), b, integrate_option());
141 std::vector<size_type> dis_inod_tab;
142 for (
size_type ie = 0, ne = omega.size(); ie < ne; ie++) {
143 const geo_element& K = omega [ie];
144 const Eigen::Matrix<piola<T>,Eigen::Dynamic,1>& piola = pops.get_piola (omega, K);
145 space_numbering::dis_inod (b, omega.sizes(), K, dis_inod_tab);
146 for (
size_type loc_inod = 0, loc_nnod = piola.size(); loc_inod < loc_nnod; ++loc_inod) {
147 _xdof.dis_entry (dis_inod_tab[loc_inod]) = piola [loc_inod].F;
150 _xdof.dis_entry_assembly();
152 if (is_distributed<M>::value && comm.size() > 1) {
153 index_set dis_inod_set;
154 for (
size_type ie = 0, ne = omega.size(); ie < ne; ++ie) {
155 const geo_element& K = omega [ie];
156 space_numbering::dis_inod (b, omega.sizes(), K, dis_inod_tab);
157 for (
size_type loc_inod = 0, loc_nnod = dis_inod_tab.size(); loc_inod < loc_nnod; ++loc_inod) {
159 if (nod_ownership.is_owned(dis_inod))
continue;
160 check_macro (dis_inod < dis_nnod,
"invalid dis_inod="<<dis_inod<<
" out of range [0:"<<dis_nnod<<
"[");
161 dis_inod_set.insert (dis_inod);
164 _xdof.set_dis_indexes (dis_inod_set);
167template <
class T,
class M>
169space_base_rep<T,M>::base_freeze_body ()
const
171 _constit.neighbour_guard();
175 disarray<size_type,M> blocked_flag = _constit.build_blocked_flag();
178 for (
size_type idof = 0, ndof = blocked_flag.size(); idof <
ndof; idof++) {
179 _idof2blk_dis_iub [idof].set_blocked (blocked_flag[idof]);
186 for (
size_type idof = 0, ndof = _idof2blk_dis_iub.size(); idof <
ndof; idof++) {
187 bool blk = _idof2blk_dis_iub [idof].is_blocked();
189 _idof2blk_dis_iub[idof].set_dis_iub (n_unknown);
192 _idof2blk_dis_iub[idof].set_dis_iub (n_blocked);
198#ifdef _RHEOLEF_HAVE_MPI
199 if (is_distributed<M>::value) {
200 dis_n_unknown = mpi::all_reduce (comm(), dis_n_unknown, std::plus<T>());
201 dis_n_blocked = mpi::all_reduce (comm(), dis_n_blocked, std::plus<T>());
204 _iu_ownership = distributor (dis_n_unknown, comm(), n_unknown);
205 _ib_ownership = distributor (dis_n_blocked, comm(), n_blocked);
210template <
class T,
class M>
212space_base_rep<T,M>::dis_idof (
const geo_element& K, std::vector<size_type>& dis_idof)
const
215 _constit.assembly_dis_idof (
get_geo(), K, dis_idof);
220template <
class T,
class M>
222space_base_rep<T,M>::name()
const
224 return _constit.name();
239template <
class T,
class M>
240disarray<typename space_base_rep<T,M>::size_type,
M>
241space_base_rep<T,M>::build_dom_dis_idof2bgd_dis_idof (
242 const space_base_rep<T,M>& Wh,
243 const std::string& dom_name
246 const geo_basic<T,M>& bgd_gamma =
get_geo()[dom_name];
247 return build_dom_dis_idof2bgd_dis_idof (Wh, bgd_gamma);
249template <
class T,
class M>
250disarray<typename space_base_rep<T,M>::size_type,
M>
251space_base_rep<T,M>::build_dom_dis_idof2bgd_dis_idof (
252 const space_base_rep<T,M>& Wh,
253 const geo_basic<T,M>& bgd_gamma2
257 const geo_basic<T,M>& bgd_gamma = bgd_gamma2.get_background_domain();
261 const space_base_rep<T,M>& Vh = *
this;
264 const geo_basic<T,M>& dom_gamma = Wh.get_geo();
266 check_macro (dom_gamma.size() == bgd_gamma.size(),
"incompatible domains");
267 distributor dom_ownership = Wh.ownership();
268 distributor bgd_ownership = Vh.ownership();
269 size_type my_proc = dom_ownership.comm().rank();
270trace_macro(
"Vh="<<Vh.name()<<
", Wh="<<Wh.name()<<
", dom_dis_idof2bgd_dis_idof.size=Wh.size="<<dom_ownership.size()<<
"...");
271 size_type first_dom_dis_idof = dom_ownership.first_index();
272 size_type first_bgd_dis_idof = bgd_ownership.first_index();
273 std::vector<size_type> dom_dis_idofs, bgd_dis_idofs;
274 disarray<size_type, M> dom_dis_idof2bgd_dis_idof (dom_ownership, std::numeric_limits<size_type>::max());
275 std::set<size_type> ext_dom_dis_idof;
276 for (
size_type ige = 0, nge = dom_gamma.size(); ige < nge; ige++) {
277 const geo_element& dom_S = dom_gamma[ige];
278 const geo_element& bgd_S = bgd_gamma[ige];
279 Wh.dis_idof (dom_S, dom_dis_idofs);
280 Vh.dis_idof (bgd_S, bgd_dis_idofs);
281trace_macro(
"ige="<<ige<<
": dom_S="<<dom_S.name()<<dom_S.dis_ie()<<
", bgd_S="<<bgd_S.name()<<bgd_S.dis_ie()
282 <<
": dom_dis_idofs.size="<<dom_dis_idofs.size()
283 <<
", bgd_dis_idofs.size="<<bgd_dis_idofs.size());
284 for (
size_type loc_idof = 0, loc_ndof = dom_dis_idofs.size(); loc_idof < loc_ndof; loc_idof++) {
285 size_type dom_dis_idof = dom_dis_idofs [loc_idof];
286 size_type bgd_dis_idof = bgd_dis_idofs [loc_idof];
287 dom_dis_idof2bgd_dis_idof.dis_entry (dom_dis_idof) = bgd_dis_idof;
288 trace_macro(
"ige="<<ige<<
": dom_dis_idof2bgd_dis_idof["<<dom_dis_idof<<
"] = "<<bgd_dis_idof);
289 if (! dom_ownership.is_owned (dom_dis_idof, my_proc)) {
290 ext_dom_dis_idof.insert (dom_dis_idof);
294 trace_macro (
"ext_dom_dis_idof.size="<<ext_dom_dis_idof.size());
295 dom_dis_idof2bgd_dis_idof.dis_entry_assembly();
296 dom_dis_idof2bgd_dis_idof.set_dis_indexes (ext_dom_dis_idof);
299 for (
size_type dom_idof = 0, dom_ndof = dom_dis_idof2bgd_dis_idof.size(); dom_idof < dom_ndof; dom_idof++) {
300 size_type bgd_dis_idof = dom_dis_idof2bgd_dis_idof [dom_idof];
301 size_type dom_dis_idof = dom_idof + first_dom_dis_idof;
302 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()<<
"[");
303 size_type bgd_idof = bgd_dis_idof - first_bgd_dis_idof;
304 dom_dis_idof2bgd_dis_idof [dom_idof] = bgd_idof;
306trace_macro(
"Vh="<<Vh.name()<<
", Wh="<<Wh.name()<<
", dom_dis_idof2bgd_dis_idof.size=Wh.size="<<dom_ownership.size()<<
" done");
308 return dom_dis_idof2bgd_dis_idof;
310#define _RHEOLEF_space_real(M) \
313space_basic<T,M>::real() \
315 return space_basic<T,M> (geo_basic<T,M>(details::zero_dimension()), "P1"); \
318#ifdef _RHEOLEF_HAVE_MPI
321#undef _RHEOLEF_space_real
373template <
class T,
class M>
377 parent_subgeo_owner_check();
378 if (! _parent_subgeo_owner_reattribution_required) {
379 return ownership().find_owner (dis_idof);
381 return _parent_subgeo_owner.dis_at (dis_idof);
384struct pair_dof_parent_owners_min_op {
385 using t = std::pair<size_t,size_t>;
386 t operator() (
const t& a,
const t& b) {
390 if (a.first == a.second)
return a;
391 if (b.first == b.second)
return b;
392 size_t parent_owner = std::min (a.first, b.first);
393 return std::make_pair (parent_owner, a.second);
398template <
class T,
class M>
400space_base_rep<T,M>::parent_subgeo_owner_guard()
const
402 if (_parent_subgeo_owner_initialized) {
403 return _parent_subgeo_owner_reattribution_required;
405 _parent_subgeo_owner_initialized =
true;
408 check_macro (! get_constitution().is_hierarchical(),
"parent_subgeo_owner_reattribution: invalid product-space");
409 const space_constitution_terminal<T,M>& terminal_constit = get_constitution().get_terminal();
410 size_type nproc = ownership().comm().size();
411 size_type my_proc = ownership().comm().rank();
412 const geo_basic<T,M>& omega = terminal_constit.get_geo();
413 _parent_subgeo_owner_reattribution_required
414 = terminal_constit.get_basis().is_continuous()
416 && omega.variant() == geo_abstract_base_rep<T>::geo_domain;
417 if (! _parent_subgeo_owner_reattribution_required) {
418 return _parent_subgeo_owner_reattribution_required;
422 std::vector<size_type> dis_idx;
423 size_type unset = std::numeric_limits<size_type>::max();
424 std::pair<size_type,size_type> unset_pair (unset, unset);
425 disarray<std::pair<size_type,size_type>> dof_parent_choice (ownership(), unset_pair);
427 std::set<size_type> ext_dis_idof;
428 for (
size_type ie = 0, ne = omega.size(); ie < ne; ++ie) {
429 const geo_element& K = omega [ie];
430 get_constitution().assembly_dis_idof (omega, K, dis_idx);
431 for (
size_type loc_idof = 0, loc_ndof = dis_idx.size(); loc_idof < loc_ndof; loc_idof++) {
433 size_type dof_proc = ownership().find_owner (dis_idof);
434 dof_parent_choice.dis_entry (dis_idof) = std::make_pair (parent_proc, dof_proc);
435 if (dof_proc != my_proc) {
436 ext_dis_idof.insert (dis_idof);
440 dof_parent_choice.dis_entry_assembly (details::pair_dof_parent_owners_min_op());
441 dof_parent_choice.set_dis_indexes (ext_dis_idof);
443 _parent_subgeo_owner.resize (ownership());
444 std::fill (_parent_subgeo_owner.begin(), _parent_subgeo_owner.end(), std::numeric_limits<size_type>::max());
445 for (
size_type ie = 0, ne = omega.size(); ie < ne; ++ie) {
446 const geo_element& K = omega [ie];
447 get_constitution().assembly_dis_idof (omega, K, dis_idx);
448 for (
size_type loc_idof = 0, loc_ndof = dis_idx.size(); loc_idof < loc_ndof; loc_idof++) {
450 size_type chosen_parent_proc = dof_parent_choice.dis_at(dis_idof).first;
453 size_type reattributed_dof_proc = (chosen_parent_proc != unset) ? chosen_parent_proc : my_proc;
454 _parent_subgeo_owner.dis_entry (dis_idof) = reattributed_dof_proc;
456 size_type dof_proc = ownership().find_owner (dis_idof);
457 if (dof_proc != my_proc) {
458 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);
464 _parent_subgeo_owner.dis_entry_assembly (details::generic_set_min_op());
465 _parent_subgeo_owner.set_dis_indexes (ext_dis_idof);
466 return _parent_subgeo_owner_reattribution_required;
471#define _RHEOLEF_instanciation(T,M) \
472template class space_base_rep<T,M>; \
473template space_basic<T,M> space_basic<T,M>::real();
476#ifdef _RHEOLEF_HAVE_MPI
479#undef _RHEOLEF_instanciation
field::size_type size_type
see the Float page for the full documentation
see the communicator page for the full documentation
space_pair_type::size_type size_type
#define trace_macro(message)
#define warning_macro(message)
void get_geo(istream &in, my_geo &omega)
check_macro(expr1.have_homogeneous_space(Xh1), "dual(expr1,expr2); expr1 should have homogeneous space. HINT: use dual(interpolate(Xh, expr1),expr2)")
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)
size_type nnod(const basis_basic< T > &b, const geo_size &gs, size_type map_dim)
size_type dis_nnod(const basis_basic< T > &b, const geo_size &gs, size_type map_dim)
size_type ndof(const basis_basic< T > &b, const geo_size &gs, size_type map_dim)
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.
#define _RHEOLEF_instanciation(T, M)
#define _RHEOLEF_space_real(M)