Rheolef  7.2
an efficient C++ finite element environment
 
Loading...
Searching...
No Matches
geo_domain_mpi.cc
Go to the documentation of this file.
1
21
22#include "rheolef/config.h"
23
24#ifdef _RHEOLEF_HAVE_MPI
25#include "rheolef/geo_domain.h"
26
27namespace rheolef {
28
29// ------------------------------------------------------------------------
30// build edge and face connectivity from domain
31// ------------------------------------------------------------------------
32template <class T>
33void
35 const domain_indirect_rep<distributed>& indirect,
36 const geo_abstract_rep<T,distributed>& bgd_omega,
37 size_type sid_dim,
38 disarray<size_type>& bgd_isid2dom_dis_isid,
39 disarray<size_type>& dom_isid2bgd_isid,
40 disarray<size_type>& dom_isid2dom_ios_dis_isid,
42{
43 if (sid_dim != 0 && base::_gs._map_dimension <= sid_dim) return;
44 communicator comm = bgd_omega.geo_element_ownership(sid_dim).comm();
45 // ------------------------------------------------------------------------
46 // 1) side compact re-numbering
47 // ------------------------------------------------------------------------
48 // 1.1 loop on elements and mark used sides
49 disarray<size_type> bgd_isid_is_on_domain (bgd_omega.geo_element_ownership(sid_dim), 0); // logical, init to "false"
50 disarray<size_type> bgd_ios_isid_is_on_domain (bgd_omega.geo_element_ios_ownership(sid_dim), 0);
51 for (size_type ioige = 0, noige = indirect.size(); ioige < noige; ioige++) {
52 size_type ige = indirect.oige (ioige).index();
53 const geo_element& bgd_K = bgd_omega.get_geo_element (base::_gs._map_dimension, ige);
54 for (size_type loc_isid = 0, loc_nsid = bgd_K.n_subgeo(sid_dim); loc_isid < loc_nsid; loc_isid++) {
55 size_type bgd_dis_isid = 0; // TODO: = bgd_K.subgeo (sid_dim, loc_isid);
56 switch (sid_dim) {
57 case 0: {
58 size_type bgd_dis_inod = bgd_K[loc_isid];
59 bgd_dis_isid = bgd_omega.dis_inod2dis_iv (bgd_dis_inod);
60 break;
61 }
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);
65 }
66 size_type bgd_ios_dis_isid = bgd_omega.dis_ige2ios_dis_ige (sid_dim, bgd_dis_isid);
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;
69 }
70 }
71 bgd_isid_is_on_domain.dis_entry_assembly(); // logical "or" across processes
72 bgd_ios_isid_is_on_domain.dis_entry_assembly();
73
74 // 1.2 counting & distribution for dom_isid
75 size_type dom_nsid = 0;
76 for (size_type bgd_isid = 0, bgd_nsid = bgd_omega.geo_element_ownership(sid_dim).size(); bgd_isid < bgd_nsid; bgd_isid++) {
77 if (bgd_isid_is_on_domain[bgd_isid] != 0) dom_nsid++ ;
78 }
79 // 1.3 counting & distribution for dom_ios_isid
80 size_type dom_ios_nsid = 0;
81 for (size_type bgd_ios_isid = 0, bgd_ios_nsid = bgd_omega.geo_element_ios_ownership(sid_dim).size(); bgd_ios_isid < bgd_ios_nsid; bgd_ios_isid++) {
82 if (bgd_ios_isid_is_on_domain[bgd_ios_isid] != 0) dom_ios_nsid++ ;
83 }
84 // 1.4 numbering dom_isid & permutation: bgd_isid <--> dom_isid
86 variant < reference_element:: last_variant_by_dimension(sid_dim); variant++) {
87 size_by_variant [variant] = 0;
88 }
89 distributor dom_isid_ownership (distributor::decide, comm, dom_nsid);
90 size_type first_dom_dis_isid = dom_isid_ownership.first_index();
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());
93 for (size_type dom_isid = 0, bgd_isid = 0, bgd_nsid = bgd_omega.geo_element_ownership(sid_dim).size(); bgd_isid < bgd_nsid; bgd_isid++) {
94 if (bgd_isid_is_on_domain[bgd_isid] == 0) continue;
95 size_type dom_dis_isid = first_dom_dis_isid + dom_isid; // distributed case !
96 bgd_isid2dom_dis_isid [bgd_isid] = dom_dis_isid;
97 dom_isid2bgd_isid [dom_isid] = bgd_isid;
98 const geo_element& bgd_S = bgd_omega.get_geo_element(sid_dim,bgd_isid);
99 size_by_variant [bgd_S.variant()]++;
100 dom_isid++;
101 }
102 // 1.5 numbering dom_ios_isid & permutation: bgd_ios_isid <--> dom_ios_isid
103 distributor dom_ios_isid_ownership (distributor::decide, comm, dom_ios_nsid);
104 size_type first_dom_ios_dis_isid = dom_ios_isid_ownership.first_index();
105 disarray<size_type> bgd_ios_isid2dom_ios_dis_isid (bgd_omega.geo_element_ios_ownership(sid_dim), std::numeric_limits<size_type>::max());
106 disarray<size_type> dom_ios_isid2bgd_ios_isid (dom_ios_isid_ownership, std::numeric_limits<size_type>::max());
107 for (size_type dom_ios_isid = 0, bgd_ios_isid = 0, bgd_ios_nsid = bgd_omega.geo_element_ios_ownership(sid_dim).size(); bgd_ios_isid < bgd_ios_nsid; bgd_ios_isid++) {
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;
112 dom_ios_isid++;
113 }
114 // 1.6 permutation: bgd_isid --> dom_ios_isid
115 disarray<size_type> bgd_isid2dom_ios_dis_isid (bgd_omega.geo_element_ownership(sid_dim), std::numeric_limits<size_type>::max());
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;
119 size_type bgd_dis_isid = bgd_omega.ios_ige2dis_ige (sid_dim, bgd_ios_isid);
120 bgd_isid2dom_ios_dis_isid.dis_entry (bgd_dis_isid) = dom_ios_dis_isid;
121 }
122 bgd_isid2dom_ios_dis_isid.dis_entry_assembly();
123
124 // 1.7 permutations: dom_ios_isid <--> dom_isid
125 dom_isid2dom_ios_dis_isid.resize (dom_isid_ownership, std::numeric_limits<size_type>::max());
126 // set reverse permutations for i/o:
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;
134 }
135 _ios_ige2dis_ige [sid_dim].dis_entry_assembly();
136 // ------------------------------------------------------------------------
137 // 2. compute sizes and resize _geo_element[variant]
138 // ------------------------------------------------------------------------
139 size_type dom_nge = 0;
140 size_type dom_dis_nge = 0;
142 variant < reference_element:: last_variant_by_dimension(sid_dim); variant++) {
143
144 distributor dom_igev_ownership (distributor::decide, comm, size_by_variant [variant]);
145 geo_element::parameter_type param (variant, 1);
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();
150 }
151 base::_gs.ownership_by_dimension[sid_dim] = distributor (dom_dis_nge, comm, dom_nge);
152
153}
154template <class T>
155void
157 const domain_indirect_rep<distributed>& indirect,
158 const geo_abstract_rep<T,distributed>& bgd_omega,
159 disarray<size_type>& bgd_iv2dom_dis_iv,
160 size_type sid_dim,
161 disarray<size_type>& bgd_isid2dom_dis_isid,
162 disarray<size_type>& dom_isid2bgd_isid,
163 disarray<size_type>& dom_isid2dom_ios_dis_isid,
165{
166 if (sid_dim != 0 && base::_gs._map_dimension <= sid_dim) return;
167 communicator comm = bgd_omega.geo_element_ownership(sid_dim).comm();
168 // ------------------------------------------------------------------------
169 // 2) set _geo_element[variant] and S.set_ios_dis_ie
170 // also reperate external vertices
171 // ------------------------------------------------------------------------
172 std::set<size_type> ext_bgd_dis_iv_set;
173 distributor dom_isid_ownership = dom_isid2bgd_isid.ownership();
174 size_type first_dom_dis_isid = dom_isid_ownership.first_index();
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];
179 const geo_element& bgd_S = bgd_omega.get_geo_element(sid_dim,bgd_isid);
180 geo_element& dom_S = get_geo_element(sid_dim,dom_isid);
181 dom_S = bgd_S;
182 dom_S.set_dis_ie (dom_dis_isid);
183 dom_S.set_ios_dis_ie (dom_ios_dis_isid);
184 if (dom_S.dimension() == 0) continue;
185 // set S face & edge : index, orientation & rotation will be set by propagate_numbering later
186 for (size_type iloc = 0, nloc = dom_S.size(); iloc < nloc; iloc++) {
187 size_type bgd_dis_inod = bgd_S[iloc];
188 size_type bgd_dis_iv = bgd_omega.dis_inod2dis_iv (bgd_dis_inod);
189 if (! bgd_omega.geo_element_ownership(0).is_owned(bgd_dis_iv)) {
190 ext_bgd_dis_iv_set.insert (bgd_dis_iv);
191 }
192 }
193 }
194 // ------------------------------------------------------------------------
195 // 3) propagate new vertex numbering in all `dom_S' new sides
196 // ------------------------------------------------------------------------
197 if (sid_dim > 0) {
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++) {
202 size_type bgd_dis_inod = dom_S[iloc];
203 size_type bgd_dis_iv = bgd_omega.dis_inod2dis_iv (bgd_dis_inod);
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;
207 }
208 }
209 }
210}
307template <class T>
308void
310 const domain_indirect_rep<distributed>& indirect,
311 const geo_abstract_rep<T,distributed>& bgd_omega,
312 std::map<size_type,size_type>& bgd_ie2dom_ie,
313 std::map<size_type,size_type>& bgd_dis_ie2dom_dis_ie)
314{
315 base::_name = bgd_omega.name() + "[" + indirect.name() + "]";
316 base::_version = 4;
317 base::_sys_coord = bgd_omega.coordinate_system();
318 base::_dimension = bgd_omega.dimension();
319 base::_piola_basis = bgd_omega.get_piola_basis();
320 base::_gs._map_dimension = indirect.map_dimension();
321 size_type map_dim = base::_gs._map_dimension;
323 std::fill (size_by_variant, size_by_variant+reference_element::max_variant, 0);
324 // ------------------------------------------------------------------------
325 // 1) _geo_element[0]: compact renumbering
326 // ------------------------------------------------------------------------
327 std::array<disarray<size_type,distributed>,4> bgd_ige2dom_dis_ige;
328 std::array<disarray<size_type,distributed>,4> dom_ige2bgd_ige;
329 disarray<size_type> dom_isid2dom_ios_dis_isid [4];
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);
336 // ------------------------------------------------------------------------
337 // 2) detect extern vertices & count elements by variants
338 // ------------------------------------------------------------------------
340 variant < reference_element:: last_variant_by_dimension(map_dim); variant++) {
341 size_by_variant [variant] = 0;
342 }
343 std::set<size_type> ext_bgd_dis_iv_set;
344 distributor ioige_ownership = indirect.ownership();
345 size_type first_dis_ioige = ioige_ownership.first_index();
346 for (size_type ioige = 0, noige = indirect.size(); ioige < noige; ioige++) {
347 size_type ige = indirect.oige (ioige).index();
348 bgd_ie2dom_ie [ige] = ioige;
349 const geo_element& bgd_K = bgd_omega.get_geo_element (map_dim, ige);
350 size_by_variant [bgd_K.variant()]++;
351 for (size_type iloc = 0; iloc < bgd_K.size(); iloc++) {
352 size_type bgd_dis_inod = bgd_K[iloc];
353 size_type bgd_dis_iv = bgd_omega.dis_inod2dis_iv (bgd_dis_inod);
354 if (! bgd_omega.geo_element_ownership(0).is_owned(bgd_dis_iv)) {
355 ext_bgd_dis_iv_set.insert (bgd_dis_iv);
356 }
357 }
358 }
359 bgd_ige2dom_dis_ige[0].append_dis_indexes (ext_bgd_dis_iv_set);
360 // ------------------------------------------------------------------------
361 // 3) _geo_element[map_dim]: compact vertex numbering
362 // ------------------------------------------------------------------------
363 // resize _geo_element[]
364 size_type dis_nge = 0;
365 size_type nge = 0;
367 variant < reference_element:: last_variant_by_dimension(map_dim); variant++) {
368
369 distributor dom_igev_ownership (distributor::decide, base::comm(), size_by_variant [variant]);
370 geo_element::parameter_type param (variant, 1);
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();
375 }
376 base::_gs.ownership_by_dimension [map_dim] = distributor (dis_nge, base::comm(), nge);
377 // ------------------------------------------------------------------------
378 // 4) count all geo_element[] and set base::_gs.ownership_by_variant[]
379 // => then can determine the node count (order > 1)
380 // ------------------------------------------------------------------------
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);
385 }
386 // ------------------------------------------------------------------------
387 // 5) count nodes & set base::_gs.node_ownership
388 // => then dis_iv2dis_inod works (used at step 6)
389 // ------------------------------------------------------------------------
390 std::array<size_type,reference_element::max_variant> loc_nnod_by_variant ;
391 reference_element::init_local_nnode_by_variant (base::order(), loc_nnod_by_variant);
392 size_type nnod = 0;
393 for (size_type variant = 0;
394 variant < reference_element::last_variant_by_dimension(base::map_dimension());
395 variant++) {
396 nnod += base::_gs.ownership_by_variant [variant].size() * loc_nnod_by_variant [variant];
397 }
398 distributor dom_node_ownership (distributor::decide, base::comm(), nnod);
399 base::_gs.node_ownership = dom_node_ownership;
400 // ------------------------------------------------------------------------
401 // 6) set _geo_element[] values
402 // ------------------------------------------------------------------------
403 for (size_type ioige = 0, noige = indirect.size(); ioige < noige; ioige++) {
404 size_type dis_ioige = first_dis_ioige + ioige;
405 size_type ige = indirect.oige (ioige).index();
406 geo_element& dom_K = get_geo_element (map_dim, ioige);
407 const geo_element& bgd_K = bgd_omega.get_geo_element (map_dim, ige);
408 dom_K = bgd_K;
409 dom_K.set_dis_ie (dis_ioige);
410 size_type ini_dis_ioige = indirect.ioige2ini_dis_ioige (ioige);
411 dom_K.set_ios_dis_ie (ini_dis_ioige);
412 // set K face & edge : index, orientation & rotation will be set by propagate_numbering later
413 for (size_type iloc = 0, nloc = dom_K.size(); iloc < nloc; iloc++) {
414 size_type bgd_dis_inod = bgd_K[iloc];
415 size_type bgd_dis_iv = bgd_omega.dis_inod2dis_iv (bgd_dis_inod);
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;
419 }
420 }
421 // reset also dom_ige2bgd_ige[map_dim] : used by _node[] for order > 1 geometries
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++) {
427 size_type bgd_ige = indirect.oige (ioige).index();
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;
431 }
432 }
433 // ------------------------------------------------------------------------
434 // 7) _geo_element[1&2]: compact renumbering
435 // _ios_ige2dis_ige[1&2]: idem
436 // ------------------------------------------------------------------------
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);
441 }
442 // ------------------------------------------------------------------------
443 // 8) set reverse permutations for i/o
444 // ------------------------------------------------------------------------
445 _ios_ige2dis_ige [map_dim] = indirect._ini_ioige2dis_ioige;
446 // ------------------------------------------------------------------------
447 // 9) set geo_size
448 // ------------------------------------------------------------------------
449 /* There is still a problem with ios_ownership_by_variant[] :
450 - first, its formaly "ini_ownership_by_variant" for domain_indirect
451 since i/o are associated to "ini" numbering for domain_indirect.
452 - physically, geo_elements of the background geo are never distributed
453 with the "ini" ownership so, the associated number of variants are unavailable
454 Notice that this "ini" domain_indirect numbering becomes the "ios" for the geo_domain.
455 Also, the new numbering is denoted by "dom" for geo_domain.
456 Thus, there are at least two solutions for geo_domain :
457 1) on non-mixed domains, the variant disarray is not necessary, and assembly communications also
458 can be avoided. Just identify the used variant related to map_dimension and store its
459 number from dom_ige_domain.size in _ios_ownership_by_variant[].
460 2) on mixed domains & for sid_dim=2,3 :
461 via the new "dom_ige" numbering, access physically to elements and store variants in:
462 disarray<size_t> dom_variant (dom_ige_ownership)
463 then, permut it using ios_ige2dom_dis_ige and get
464 disarray<size_t> ios_variant (ios_ige_ownership)
465 finaly, count variants and store it in
466 _ios_ownership_by_variant[]
467 */
468 size_type dis_size_by_variant [reference_element::max_variant];
469 std::fill (dis_size_by_variant, dis_size_by_variant+reference_element::max_variant, 0);
470 mpi::all_reduce (base::comm(), size_by_variant, reference_element::max_variant, dis_size_by_variant, std::plus<size_type>());
471
472 size_type ios_size_by_variant [reference_element::max_variant];
473 std::fill (ios_size_by_variant, ios_size_by_variant+reference_element::max_variant, 0);
474 ios_size_by_variant [reference_element::p] = _ios_ige2dis_ige[0].ownership().size();
475 ios_size_by_variant [reference_element::e] = _ios_ige2dis_ige[1].ownership().size();
476
477 std::vector<size_type> loc_ios_size_by_variant_by_proc [reference_element::max_variant];
478 std::vector<size_type> ios_size_by_variant_by_proc [reference_element::max_variant];
479 for (size_type dim = 2; dim <= map_dim; dim++) {
480 bool is_mixed = ((dim == 2) &&
481 (dis_size_by_variant[reference_element::t] != 0 &&
482 dis_size_by_variant[reference_element::q] != 0))
483 ||
484 ((dim == 3) &&
485 ((dis_size_by_variant[reference_element::T] != 0 &&
486 dis_size_by_variant[reference_element::P] != 0) ||
487 (dis_size_by_variant[reference_element::P] != 0 &&
488 dis_size_by_variant[reference_element::H] != 0) ||
489 (dis_size_by_variant[reference_element::H] != 0 &&
490 dis_size_by_variant[reference_element::T] != 0)));
491 if (!is_mixed) {
492 switch (dim) {
493 case 2:
494 if ( dis_size_by_variant[reference_element::t] != 0) {
495 ios_size_by_variant[reference_element::t] = _ios_ige2dis_ige[2].ownership().size();
496 } else {
497 ios_size_by_variant[reference_element::q] = _ios_ige2dis_ige[2].ownership().size();
498 }
499 break;
500 case 3:
501 default:
502 if ( dis_size_by_variant[reference_element::T] != 0) {
503 ios_size_by_variant[reference_element::T] = _ios_ige2dis_ige[3].ownership().size();
504 } else
505 if ( dis_size_by_variant[reference_element::P] != 0) {
506 ios_size_by_variant[reference_element::P] = _ios_ige2dis_ige[3].ownership().size();
507 } else {
508 ios_size_by_variant[reference_element::H] = _ios_ige2dis_ige[3].ownership().size();
509 }
510 break;
511 }
512 continue;
513 }
514 // here we have a mixed domain:
515 size_type nproc = base::comm().size();
516 size_type my_proc = base::comm().rank();
517 switch (dim) {
518 case 2:
519 ios_size_by_variant_by_proc [reference_element::t].resize(nproc, 0);
520 ios_size_by_variant_by_proc [reference_element::q].resize(nproc, 0);
521 loc_ios_size_by_variant_by_proc [reference_element::t].resize(nproc, 0);
522 loc_ios_size_by_variant_by_proc [reference_element::q].resize(nproc, 0);
523 break;
524 case 3:
525 default:
526 ios_size_by_variant_by_proc [reference_element::T].resize(nproc, 0);
527 ios_size_by_variant_by_proc [reference_element::P].resize(nproc, 0);
528 ios_size_by_variant_by_proc [reference_element::H].resize(nproc, 0);
529 loc_ios_size_by_variant_by_proc [reference_element::T].resize(nproc, 0);
530 loc_ios_size_by_variant_by_proc [reference_element::P].resize(nproc, 0);
531 loc_ios_size_by_variant_by_proc [reference_element::H].resize(nproc, 0);
532 break;
533 }
534 const distributor& ios_ige_ownership_dim = _ios_ige2dis_ige[dim].ownership(); // has been set by domain_set_sides
535 for (size_type ie = 0, ne = base::sizes().ownership_by_dimension[dim].size(); ie < ne; ie++) {
536 const geo_element& K = get_geo_element (dim, ie);
537 size_type ios_dis_ie = K.ios_dis_ie();
538 size_type iproc = ios_ige_ownership_dim.find_owner(ios_dis_ie);
539 loc_ios_size_by_variant_by_proc [K.variant()][iproc]++;
540 }
541 switch (dim) {
542 case 2:
543 mpi::all_reduce (base::comm(),
544 loc_ios_size_by_variant_by_proc [reference_element::t].begin().operator->(), nproc,
545 ios_size_by_variant_by_proc [reference_element::t].begin().operator->(), std::plus<size_type>());
546 mpi::all_reduce (base::comm(),
547 loc_ios_size_by_variant_by_proc [reference_element::q].begin().operator->(), nproc,
548 ios_size_by_variant_by_proc [reference_element::q].begin().operator->(), std::plus<size_type>());
549 ios_size_by_variant[reference_element::t] = ios_size_by_variant_by_proc [reference_element::t][my_proc];
550 ios_size_by_variant[reference_element::q] = ios_size_by_variant_by_proc [reference_element::q][my_proc];
551 break;
552 case 3:
553 default:
554 // TODO: TPH
555 error_macro ("3D domain \""<<base::_name<<"\" with mixed element variants: not yet");
556 break;
557 }
558 // end of domain with mixed elements
559 }
560 // then, ios_size_by_variant[] is completed and we set ios_sizes:
561 for (size_type dim = 0; dim <= base::_gs._map_dimension; dim++) {
562 size_type first_ios_v = 0;
564 variant < reference_element:: last_variant_by_dimension(dim); variant++) {
565
566 _ios_gs.ownership_by_variant [variant] = distributor (distributor::decide, base::comm(), ios_size_by_variant [variant]);
567 _ios_gs.first_by_variant [variant] = distributor (distributor::decide, base::comm(), first_ios_v);
568 first_ios_v += ios_size_by_variant [variant];
569 }
570 size_type ios_nge = first_ios_v;
571 _ios_gs.ownership_by_dimension [dim] = distributor (distributor::decide, base::comm(), ios_nge);
572 }
573 // ------------------------------------------------------------------------
574 // 10) resize _igev2ios_dis_igev : used by Pk_numbering
575 // ------------------------------------------------------------------------
576 for (size_type variant = 0;
577 variant < reference_element:: last_variant_by_dimension(base::_gs._map_dimension);
578 variant++) {
579 _igev2ios_dis_igev [variant].resize (
580 base::_gs.ownership_by_variant [variant],
581 std::numeric_limits<size_type>::max());
582 }
583 // for dim=0,1 : no variants and igev2ios_dis_igev[] = ige2ios_dis_ige[]
584 for (size_type dim = 0; dim <= base::_gs._map_dimension; dim++) {
585 size_type ige = 0;
588 variant++) {
589 for (size_type igev = 0, ngev = _igev2ios_dis_igev [variant].size(); igev < ngev; igev++, ige++) {
590 const geo_element& K = base::_geo_element [variant][igev];
591 size_type ios_dis_ige = K.ios_dis_ie();
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;
596 size_type first_v = _ios_gs.first_by_variant [variant].size (iproc);
597 assert_macro (ios_ige >= first_v, "invalid index");
598 size_type ios_igev = ios_ige - first_v;
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;
602 }
603#ifdef TODO
604 _ios_igev2dis_igev [variant].resize (
605 _ios_gs.ownership_by_variant [variant], // this ownership is not yet set
606 std::numeric_limits<size_type>::max());
607 _igev2ios_dis_igev [variant].reverse_permutation (_ios_igev2dis_igev [variant]); // not used ?
608#endif // TODO
609 }
610 }
611 // ------------------------------------------------------------------------
612 // 11) raw copy _node
613 // TODO: use shallow copy & indirection :
614 // node(inod) { return bgd_omega.node(inod2bgd_inod(inod))}
615 // ------------------------------------------------------------------------
616 // 11.a) resize _node[]
617 set_ios_permutation (_inod2ios_dis_inod);
618 distributor ios_dom_node_ownership (dom_node_ownership.dis_size(), base::comm(), distributor::decide);
619 _ios_inod2dis_inod.resize (ios_dom_node_ownership);
620 _inod2ios_dis_inod.reverse_permutation (_ios_inod2dis_inod);
621
622 // 11.b) copy _node[] from bgd_omega.node[]
623 base::_node.resize (dom_node_ownership);
624 disarray<size_type> dom_inod2bgd_inod (dom_node_ownership, std::numeric_limits<size_type>::max());
625 size_type dom_inod = 0;
626 size_type first_bgd_inod_v = 0;
627 for (size_type dim = 0; dim <= base::_gs._map_dimension; dim++) {
628 size_type dom_ige = 0;
631 variant++) {
632 if (loc_nnod_by_variant [variant] == 0) continue;
633 size_type first_bgd_v = bgd_omega.sizes().first_by_variant [variant].size();
634 for (size_type dom_igev = 0, dom_ngev = base::_geo_element [variant].size(); dom_igev < dom_ngev; dom_igev++, dom_ige++) {
635 const geo_element& K = base::_geo_element [variant][dom_igev];
636 size_type bgd_ige = dom_ige2bgd_ige [dim][dom_ige];
637 assert_macro (bgd_ige >= first_bgd_v, "invalid index");
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);
643 }
644 }
645 first_bgd_inod_v += bgd_omega.sizes().ownership_by_variant [variant].size() * loc_nnod_by_variant [variant];
646 }
647 }
648 // ------------------------------------------------------------------------
649 // 12) propagate new edge & face numbering to all `dom_S' elements
650 // ------------------------------------------------------------------------
651 for (size_type dim = 1; dim < base::_gs._map_dimension; dim++) {
652 set_element_side_index (dim);
653 }
654 // ------------------------------------------------------------------------
655 // 13) reset vertex P[0] value (useful for band zero vertex domain)
656 // note: should be before "build_external_entities" as some vertex will be exported
657 // ------------------------------------------------------------------------
658 size_type first_dom_dis_iv = base::_geo_element[reference_element::p].ownership().first_index();
659 for (size_type dom_iv = 0, dom_nv = base::_geo_element[reference_element::p].size(); dom_iv < dom_nv; dom_iv++) {
660 geo_element& P = base::_geo_element[reference_element::p][dom_iv];
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);
663 P[0] = dom_dis_inod;
664 }
665 // ------------------------------------------------------------------------
666 // 17) bgd_dis_ie2dom_dis_ie, used by space("square[sides]",Pk) for HDG lambda multiplier
667 // ------------------------------------------------------------------------
668 // step 0: get direct access to bgd_omega _geo_element table
669 const geo_base_rep<T,distributed>* ptr_bgd_omega = dynamic_cast<const geo_base_rep<T,distributed>*>(&bgd_omega);
670 check_macro (ptr_bgd_omega != 0, "invalid bgd_omega");
671 const geo_base_rep<T,distributed>& bgd_omega1 = *ptr_bgd_omega;
672
673 // step 1: list exported elements from bgd_omega
674 size_type map_d = base::_gs._map_dimension;
675 // size_type first_bgd_igev_by_variant = 0;
676 index_set ext_bgd_dis_ie_set;
678 variant < reference_element:: last_variant_by_dimension (map_d); ++variant) {
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) {
681 size_type bgd_dis_igev = x.first;
682 const geo_element_auto<>& bgd_K = x.second;
683 size_type bgd_dis_ie = bgd_K.dis_ie(); // = bgd_dis_igev + first_bgd_igev_by_variant;
684 ext_bgd_dis_ie_set.insert (bgd_dis_ie);
685 }
686 // first_bgd_igev_by_variant += bgd_omega1._geo_element [variant].dis_size();
687 }
688 bgd_ige2dom_dis_ige[map_d].append_dis_indexes (ext_bgd_dis_ie_set);
689
690 // step 2: copy bgd_ige2dom_dis_ige[map_d].externals since this array is a temporary
691 // also: build ext_dom_dis_igev_set [variant], for each variant for map_d
692 index_set ext_dom_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) {
696 size_type bgd_dis_ie = x.first;
697 size_type dom_dis_ie = x.second;
698 if (dom_dis_ie == std::numeric_limits<size_type>::max()) {
699 // bgd_dis_ie is not part of the present domain
700 continue;
701 }
702 bgd_dis_ie2dom_dis_ie [bgd_dis_ie] = dom_dis_ie;
703 ext_dom_dis_ie_set.insert (dom_dis_ie);
704 // convert dom_dis_ie to dom_dis_igev and get its variant:
705 size_type variant;
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);
708 }
709 // step 3: communicate domain external elements
711 variant < reference_element:: last_variant_by_dimension (map_d); ++variant) {
712 base::_geo_element [variant].append_dis_indexes (ext_dom_dis_igev_set[variant]);
713 }
714 // ------------------------------------------------------------------------
715 // 14) external entities
716 // ------------------------------------------------------------------------
717 build_external_entities();
718 // ------------------------------------------------------------------------
719 // 15) domains : intersection with the current domain
720 // ------------------------------------------------------------------------
721 const geo_base_rep<T,distributed> *ptr = dynamic_cast<const geo_base_rep<T,distributed>*> (&bgd_omega);
722 check_macro (ptr != 0, "cannot build domains on \""<<base::_name<<"\"");
723 const geo_base_rep<T,distributed>& bgd_omega2 = *ptr;
724 for (size_type idom = 0, ndom = bgd_omega2.n_domain_indirect(); idom < ndom; ++idom) {
725 const domain_indirect_basic<distributed>& bgd_dom = bgd_omega2.get_domain_indirect(idom);
726 if (bgd_dom.name() == indirect.name()) continue; // myself
727 size_type dom_map_dim = bgd_dom.map_dimension();
728 if (dom_map_dim > base::_gs._map_dimension) continue; // dom has upper dimension
729 std::vector<size_type> ie_list;
730 size_type first_dom_dis_ige = base::_gs.ownership_by_dimension[dom_map_dim].first_index();
731 for (domain_indirect_basic<distributed>::const_iterator_ioige iter = bgd_dom.ioige_begin(), last = bgd_dom.ioige_end(); iter != last; ++iter) {
732 const geo_element_indirect& ioige = *iter;
733 size_type bgd_ige = ioige.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; // side do not belongs to dom
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);
739 }
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) {
742 continue; // empty intersection
743 }
744 domain_indirect_basic<distributed> dom (*this, bgd_dom.name(), bgd_dom.map_dimension(), base::comm(), ie_list);
745 base::_domains.push_back (dom);
746 }
747 // ------------------------------------------------------------------------
748 // 16) some post treatments (after build_external_entities and so one!)
749 // ------------------------------------------------------------------------
750 base::compute_bbox();
751}
752// ----------------------------------------------------------------------------
753// instanciation in library
754// ----------------------------------------------------------------------------
755template class geo_rep<Float,distributed>;
756
757} // namespace rheolef
758#endif // _RHEOLEF_HAVE_MPI
see the disarray page for the full documentation
Definition disarray.h:497
see the distributor page for the full documentation
Definition distributor.h:69
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
static const size_type decide
Definition distributor.h:83
the finite element boundary domain
abstract interface class
Definition geo.h:401
base class for M=sequential or distributed meshes representations
Definition geo.h:528
const domain_indirect_basic< M > & get_domain_indirect(size_type i) const
Definition geo.h:597
size_type n_domain_indirect() const
Definition geo.h:595
std::array< hack_array< geo_element_hack, M >, reference_element::max_variant > _geo_element
Definition geo.h:678
see the geo_element page for the full documentation
size_type edge(size_type i) const
size_type face(size_type i) const
size_type size() 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
size_type dis_ie() const
void set_dis_ie(size_type dis_ie)
base::size_type size_type
Definition geo.h:934
sequential mesh representation
Definition geo.h:778
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)
Definition dis_macros.h:113
#define error_macro(message)
Definition dis_macros.h:49
check_macro(expr1.have_homogeneous_space(Xh1), "dual(expr1,expr2); expr1 should have homogeneous space. HINT: use dual(interpolate(Xh, expr1),expr2)")
This file is part of Rheolef.