22#include "rheolef/config.h"
24#ifdef _RHEOLEF_HAVE_MPI
25#include "rheolef/geo_domain.h"
43 if (sid_dim != 0 && base::_gs._map_dimension <= sid_dim)
return;
51 for (
size_type ioige = 0, noige = indirect.
size(); ioige < noige; ioige++) {
54 for (
size_type loc_isid = 0, loc_nsid = bgd_K.
n_subgeo(sid_dim); loc_isid < loc_nsid; loc_isid++) {
62 case 1: bgd_dis_isid = bgd_K.
edge(loc_isid);
break;
63 case 2: bgd_dis_isid = bgd_K.
face(loc_isid);
break;
64 default:
error_macro (
"domain: unexpected side dimension " << sid_dim);
67 bgd_isid_is_on_domain.dis_entry (bgd_dis_isid) += 1;
68 bgd_ios_isid_is_on_domain.dis_entry (bgd_ios_dis_isid) += 1;
71 bgd_isid_is_on_domain.dis_entry_assembly();
72 bgd_ios_isid_is_on_domain.dis_entry_assembly();
77 if (bgd_isid_is_on_domain[bgd_isid] != 0) dom_nsid++ ;
82 if (bgd_ios_isid_is_on_domain[bgd_ios_isid] != 0) dom_ios_nsid++ ;
91 bgd_isid2dom_dis_isid.resize (bgd_omega.
geo_element_ownership(sid_dim), std::numeric_limits<size_type>::max());
92 dom_isid2bgd_isid.resize (dom_isid_ownership, std::numeric_limits<size_type>::max());
94 if (bgd_isid_is_on_domain[bgd_isid] == 0)
continue;
95 size_type dom_dis_isid = first_dom_dis_isid + dom_isid;
96 bgd_isid2dom_dis_isid [bgd_isid] = dom_dis_isid;
97 dom_isid2bgd_isid [dom_isid] = bgd_isid;
99 size_by_variant [bgd_S.
variant()]++;
106 disarray<size_type> dom_ios_isid2bgd_ios_isid (dom_ios_isid_ownership, std::numeric_limits<size_type>::max());
108 if (bgd_ios_isid_is_on_domain[bgd_ios_isid] == 0)
continue;
109 size_type dom_ios_dis_isid = first_dom_ios_dis_isid + dom_ios_isid;
110 bgd_ios_isid2dom_ios_dis_isid [bgd_ios_isid] = dom_ios_dis_isid;
111 dom_ios_isid2bgd_ios_isid [dom_ios_isid] = bgd_ios_isid;
116 for (
size_type dom_ios_isid = 0; dom_ios_isid < dom_ios_nsid; dom_ios_isid++) {
117 size_type bgd_ios_isid = dom_ios_isid2bgd_ios_isid [dom_ios_isid];
118 size_type dom_ios_dis_isid = dom_ios_isid + first_dom_ios_dis_isid;
120 bgd_isid2dom_ios_dis_isid.dis_entry (bgd_dis_isid) = dom_ios_dis_isid;
122 bgd_isid2dom_ios_dis_isid.dis_entry_assembly();
125 dom_isid2dom_ios_dis_isid.resize (dom_isid_ownership, std::numeric_limits<size_type>::max());
127 _ios_ige2dis_ige [sid_dim].resize (dom_ios_isid_ownership, std::numeric_limits<size_type>::max());
128 for (
size_type dom_isid = 0; dom_isid < dom_nsid; dom_isid++) {
129 size_type bgd_isid = dom_isid2bgd_isid [dom_isid];
130 size_type dom_dis_isid = dom_isid + first_dom_dis_isid;
131 size_type dom_ios_dis_isid = bgd_isid2dom_ios_dis_isid [bgd_isid];
132 dom_isid2dom_ios_dis_isid [dom_isid] = dom_ios_dis_isid;
133 _ios_ige2dis_ige [sid_dim].dis_entry (dom_ios_dis_isid) = dom_dis_isid;
135 _ios_ige2dis_ige [sid_dim].dis_entry_assembly();
146 base::_geo_element[
variant].resize (dom_igev_ownership, param);
147 base::_gs.ownership_by_variant[
variant] = dom_igev_ownership;
148 dom_nge += dom_igev_ownership.
size();
149 dom_dis_nge += dom_igev_ownership.
dis_size();
151 base::_gs.ownership_by_dimension[sid_dim] =
distributor (dom_dis_nge, comm, dom_nge);
166 if (sid_dim != 0 && base::_gs._map_dimension <= sid_dim)
return;
172 std::set<size_type> ext_bgd_dis_iv_set;
173 distributor dom_isid_ownership = dom_isid2bgd_isid.ownership();
175 for (
size_type dom_isid = 0, dom_nsid = dom_isid_ownership.
size(); dom_isid < dom_nsid; dom_isid++) {
176 size_type bgd_isid = dom_isid2bgd_isid [dom_isid];
177 size_type dom_dis_isid = first_dom_dis_isid + dom_isid;
178 size_type dom_ios_dis_isid = dom_isid2dom_ios_dis_isid [dom_isid];
180 geo_element& dom_S = get_geo_element(sid_dim,dom_isid);
186 for (
size_type iloc = 0, nloc = dom_S.
size(); iloc < nloc; iloc++) {
190 ext_bgd_dis_iv_set.insert (bgd_dis_iv);
198 bgd_iv2dom_dis_iv.append_dis_indexes (ext_bgd_dis_iv_set);
199 for (
size_type dom_isid = 0, dom_nsid = dom_isid_ownership.
size(); dom_isid < dom_nsid; dom_isid++) {
200 geo_element& dom_S = get_geo_element(sid_dim,dom_isid);
201 for (
size_type iloc = 0, nloc = dom_S.
size(); iloc < nloc; iloc++) {
204 size_type dom_dis_iv = bgd_iv2dom_dis_iv.dis_at (bgd_dis_iv);
205 size_type dom_dis_inod = base::dis_iv2dis_inod (dom_dis_iv);
206 dom_S[iloc] = dom_dis_inod;
312 std::map<size_type,size_type>& bgd_ie2dom_ie,
313 std::map<size_type,size_type>& bgd_dis_ie2dom_dis_ie)
315 base::_name = bgd_omega.
name() +
"[" + indirect.
name() +
"]";
318 base::_dimension = bgd_omega.
dimension();
327 std::array<disarray<size_type,distributed>,4> bgd_ige2dom_dis_ige;
328 std::array<disarray<size_type,distributed>,4> dom_ige2bgd_ige;
330 domain_set_side_part1 (indirect, bgd_omega, 0,
331 bgd_ige2dom_dis_ige[0], dom_ige2bgd_ige[0],
332 dom_isid2dom_ios_dis_isid [0], size_by_variant);
333 domain_set_side_part2 (indirect, bgd_omega, bgd_ige2dom_dis_ige[0], 0,
334 bgd_ige2dom_dis_ige[0], dom_ige2bgd_ige[0],
335 dom_isid2dom_ios_dis_isid [0], size_by_variant);
343 std::set<size_type> ext_bgd_dis_iv_set;
344 distributor ioige_ownership = indirect.ownership();
346 for (
size_type ioige = 0, noige = indirect.
size(); ioige < noige; ioige++) {
348 bgd_ie2dom_ie [ige] = ioige;
350 size_by_variant [bgd_K.
variant()]++;
355 ext_bgd_dis_iv_set.insert (bgd_dis_iv);
359 bgd_ige2dom_dis_ige[0].append_dis_indexes (ext_bgd_dis_iv_set);
371 base::_geo_element[
variant].resize (dom_igev_ownership, param);
372 base::_gs.ownership_by_variant [
variant] = dom_igev_ownership;
373 dis_nge += dom_igev_ownership.
dis_size();
374 nge += dom_igev_ownership.
size();
376 base::_gs.ownership_by_dimension [
map_dim] =
distributor (dis_nge, base::comm(), nge);
381 for (
size_type sid_dim = 1; sid_dim < base::_gs._map_dimension; sid_dim++) {
382 domain_set_side_part1 (indirect, bgd_omega, sid_dim,
383 bgd_ige2dom_dis_ige[sid_dim], dom_ige2bgd_ige[sid_dim],
384 dom_isid2dom_ios_dis_isid [sid_dim], size_by_variant);
390 std::array<size_type,reference_element::max_variant> loc_nnod_by_variant ;
396 nnod += base::_gs.ownership_by_variant [
variant].size() * loc_nnod_by_variant [
variant];
399 base::_gs.node_ownership = dom_node_ownership;
403 for (
size_type ioige = 0, noige = indirect.
size(); ioige < noige; ioige++) {
404 size_type dis_ioige = first_dis_ioige + ioige;
413 for (
size_type iloc = 0, nloc = dom_K.
size(); iloc < nloc; iloc++) {
416 size_type dom_dis_iv = bgd_ige2dom_dis_ige[0].dis_at (bgd_dis_iv);
417 size_type dom_dis_inod = base::dis_iv2dis_inod (dom_dis_iv);
418 dom_K[iloc] = dom_dis_inod;
422 if (base::_gs._map_dimension > 0) {
423 dom_ige2bgd_ige [base::_gs._map_dimension].resize(indirect.ownership());
424 bgd_ige2dom_dis_ige[base::_gs._map_dimension].resize(bgd_omega.
geo_element_ownership(base::_gs._map_dimension), std::numeric_limits<size_type>::max());
425 size_type first_dis_ioige = indirect.ownership().first_index();
426 for (
size_type ioige = 0, noige = indirect.
size(); ioige < noige; ioige++) {
428 size_type dis_ioige = first_dis_ioige + ioige;
429 dom_ige2bgd_ige [base::_gs._map_dimension] [ioige] = bgd_ige;
430 bgd_ige2dom_dis_ige[base::_gs._map_dimension] [bgd_ige] = dis_ioige;
437 for (
size_type sid_dim = 1; sid_dim < base::_gs._map_dimension; sid_dim++) {
438 domain_set_side_part2 (indirect, bgd_omega, bgd_ige2dom_dis_ige[0], sid_dim,
439 bgd_ige2dom_dis_ige[sid_dim], dom_ige2bgd_ige[sid_dim],
440 dom_isid2dom_ios_dis_isid [sid_dim], size_by_variant);
480 bool is_mixed = ((
dim == 2) &&
534 const distributor& ios_ige_ownership_dim = _ios_ige2dis_ige[
dim].ownership();
535 for (
size_type ie = 0, ne = base::sizes().ownership_by_dimension[
dim].size(); ie < ne; ie++) {
539 loc_ios_size_by_variant_by_proc [K.
variant()][iproc]++;
543 mpi::all_reduce (base::comm(),
545 ios_size_by_variant_by_proc [
reference_element::t].begin().operator->(), std::plus<size_type>());
546 mpi::all_reduce (base::comm(),
548 ios_size_by_variant_by_proc [
reference_element::q].begin().operator->(), std::plus<size_type>());
555 error_macro (
"3D domain \""<<base::_name<<
"\" with mixed element variants: not yet");
568 first_ios_v += ios_size_by_variant [
variant];
579 _igev2ios_dis_igev [
variant].resize (
580 base::_gs.ownership_by_variant [
variant],
581 std::numeric_limits<size_type>::max());
589 for (
size_type igev = 0, ngev = _igev2ios_dis_igev [
variant].size(); igev < ngev; igev++, ige++) {
592 size_type iproc = _ios_gs.ownership_by_dimension [
dim].find_owner(ios_dis_ige);
593 size_type first_ios_dis_ige = _ios_gs.ownership_by_dimension [
dim].first_index (iproc);
594 assert_macro (ios_dis_ige >= first_ios_dis_ige,
"invalid index");
595 size_type ios_ige = ios_dis_ige - first_ios_dis_ige;
599 size_type first_ios_dis_igev = _ios_gs.ownership_by_variant [
variant].first_index (iproc);
600 size_type ios_dis_igev = first_ios_dis_igev + ios_igev;
601 _igev2ios_dis_igev [
variant][igev] = ios_dis_igev;
604 _ios_igev2dis_igev [
variant].resize (
605 _ios_gs.ownership_by_variant [
variant],
606 std::numeric_limits<size_type>::max());
607 _igev2ios_dis_igev [
variant].reverse_permutation (_ios_igev2dis_igev [
variant]);
617 set_ios_permutation (_inod2ios_dis_inod);
619 _ios_inod2dis_inod.resize (ios_dom_node_ownership);
620 _inod2ios_dis_inod.reverse_permutation (_ios_inod2dis_inod);
623 base::_node.resize (dom_node_ownership);
624 disarray<size_type> dom_inod2bgd_inod (dom_node_ownership, std::numeric_limits<size_type>::max());
632 if (loc_nnod_by_variant [
variant] == 0)
continue;
634 for (
size_type dom_igev = 0, dom_ngev = base::_geo_element [
variant].size(); dom_igev < dom_ngev; dom_igev++, dom_ige++) {
638 size_type bgd_igev = bgd_ige - first_bgd_v;
639 for (
size_type loc_inod = 0, loc_nnod = loc_nnod_by_variant [
variant]; loc_inod < loc_nnod; loc_inod++, dom_inod++) {
640 size_type bgd_inod = first_bgd_inod_v + bgd_igev * loc_nnod_by_variant [
variant] + loc_inod;
641 dom_inod2bgd_inod [dom_inod] = bgd_inod;
642 base::_node [dom_inod] = bgd_omega.
node (bgd_inod);
652 set_element_side_index (
dim);
661 size_type dom_dis_iv = first_dom_dis_iv + dom_iv;
662 size_type dom_dis_inod = base::dis_iv2dis_inod (dom_dis_iv);
670 check_macro (ptr_bgd_omega != 0,
"invalid bgd_omega");
674 size_type map_d = base::_gs._map_dimension;
679 const std::map<size_type,geo_element_auto<>>& ext_bgd_gev = bgd_omega1.
_geo_element [
variant].get_dis_map_entries();
680 for (
auto x: ext_bgd_gev) {
684 ext_bgd_dis_ie_set.
insert (bgd_dis_ie);
688 bgd_ige2dom_dis_ige[map_d].append_dis_indexes (ext_bgd_dis_ie_set);
693 const std::map <size_type, size_type>& bgd_dis_ie2dom_dis_ie_tmp = bgd_ige2dom_dis_ige[map_d].get_dis_map_entries();
694 std::array<index_set,reference_element::max_variant> ext_dom_dis_igev_set;
695 for (
auto x : bgd_dis_ie2dom_dis_ie_tmp) {
698 if (dom_dis_ie == std::numeric_limits<size_type>::max()) {
702 bgd_dis_ie2dom_dis_ie [bgd_dis_ie] = dom_dis_ie;
703 ext_dom_dis_ie_set.
insert (dom_dis_ie);
706 size_type dom_dis_igev = base::sizes().dis_ige2dis_igev_by_dimension (map_d, dom_dis_ie,
variant);
707 ext_dom_dis_igev_set [
variant].insert (dom_dis_igev);
712 base::_geo_element [
variant].append_dis_indexes (ext_dom_dis_igev_set[
variant]);
717 build_external_entities();
722 check_macro (ptr != 0,
"cannot build domains on \""<<base::_name<<
"\"");
726 if (bgd_dom.
name() == indirect.
name())
continue;
728 if (dom_map_dim > base::_gs._map_dimension)
continue;
729 std::vector<size_type> ie_list;
730 size_type first_dom_dis_ige = base::_gs.ownership_by_dimension[dom_map_dim].first_index();
734 size_type dom_dis_ige = bgd_ige2dom_dis_ige [dom_map_dim][bgd_ige];
735 if (dom_dis_ige == std::numeric_limits<size_type>::max())
continue;
736 check_macro (dom_dis_ige >= first_dom_dis_ige,
"invalid index");
737 size_type dom_ige = dom_dis_ige - first_dom_dis_ige;
738 ie_list.push_back(dom_ige);
740 size_type ie_list_dis_size = mpi::all_reduce (base::comm(), ie_list.size(), std::plus<size_type>());
741 if (ie_list_dis_size == 0) {
745 base::_domains.push_back (dom);
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 dis_size() const
global and local sizes
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)[
static const size_type decide
const communicator_type & comm() const
const_iterator_ioige ioige_end() const
const_iterator_ioige ioige_begin() const
size_type map_dimension() const
the finite element boundary domain
size_type ioige2ini_dis_ioige(size_type ioige) const
disarray< size_type, distributed > _ini_ioige2dis_ioige
size_type map_dimension() const
const geo_element_indirect & oige(size_type ioige) const
virtual const node_type & node(size_type inod) const =0
virtual const basis_basic< T > & get_piola_basis() const =0
virtual const distributor & geo_element_ownership(size_type dim) const =0
virtual std::string name() const =0
virtual const_reference get_geo_element(size_type dim, size_type ige) const =0
virtual const geo_size & sizes() const =0
virtual size_type dis_inod2dis_iv(size_type dis_inod) const =0
virtual size_type dimension() const =0
virtual coordinate_type coordinate_system() const =0
virtual size_type ios_ige2dis_ige(size_type dim, size_type ios_ige) const =0
virtual distributor geo_element_ios_ownership(size_type dim) const =0
virtual size_type dis_ige2ios_dis_ige(size_type dim, size_type dis_ige) const =0
const domain_indirect_basic< M > & get_domain_indirect(size_type i) const
size_type n_domain_indirect() const
std::array< hack_array< geo_element_hack, M >, reference_element::max_variant > _geo_element
see the geo_element page for the full documentation
size_type edge(size_type i) const
size_type face(size_type i) const
void set_ios_dis_ie(size_type ios_dis_ie)
size_type dimension() const
size_type ios_dis_ie() const
variant_type variant() const
size_type n_subgeo(size_type subgeo_dim) const
void set_dis_ie(size_type dis_ie)
base::size_type size_type
sequential mesh representation
void insert(size_type dis_i)
static const variant_type H
static const variant_type q
static const variant_type e
static const variant_type max_variant
static void init_local_nnode_by_variant(size_type order, std::array< size_type, reference_element::max_variant > &loc_nnod_by_variant)
static variant_type last_variant_by_dimension(size_type dim)
static variant_type first_variant_by_dimension(size_type dim)
static const variant_type p
static const variant_type T
static const variant_type P
static const variant_type t
#define assert_macro(ok_condition, message)
#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)")
size_type nnod(const basis_basic< T > &b, const geo_size &gs, size_type map_dim)
This file is part of Rheolef.
void compute_bbox(const geo_base_rep< T, M > &omega, const geo_element &K, point_basic< T > &xmin, point_basic< T > &xmax)
distributor ownership_by_variant[reference_element::max_variant]
distributor first_by_variant[reference_element::max_variant]