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);
758 #endif // _RHEOLEF_HAVE_MPI