Rheolef  7.1
an efficient C++ finite element environment
geo_domain_mpi.cc
Go to the documentation of this file.
1 
22 #include "rheolef/config.h"
23 
24 #ifdef _RHEOLEF_HAVE_MPI
25 #include "rheolef/geo_domain.h"
26 
27 namespace rheolef {
28 
29 // ------------------------------------------------------------------------
30 // build edge and face connectivity from domain
31 // ------------------------------------------------------------------------
32 template <class T>
33 void
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,
41  size_type size_by_variant [reference_element::max_variant])
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
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;
143 
144  distributor dom_igev_ownership (distributor::decide, comm, size_by_variant [variant]);
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 }
154 template <class T>
155 void
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,
164  size_type size_by_variant [reference_element::max_variant])
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 }
307 template <class T>
308 void
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;
322  size_type size_by_variant [reference_element::max_variant];
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  // ------------------------------------------------------------------------
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;
368 
369  distributor dom_igev_ownership (distributor::decide, base::comm(), size_by_variant [variant]);
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 ;
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;
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;
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:
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
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  // ------------------------------------------------------------------------
751 }
752 // ----------------------------------------------------------------------------
753 // instanciation in library
754 // ----------------------------------------------------------------------------
755 template class geo_rep<Float,distributed>;
756 
757 } // namespace rheolef
758 #endif // _RHEOLEF_HAVE_MPI
rheolef::reference_element::last_variant_by_dimension
static variant_type last_variant_by_dimension(size_type dim)
Definition: reference_element.h:150
rheolef::domain_indirect_rep< distributed >::size
size_type size() const
Definition: domain_indirect.h:264
rheolef::domain_indirect_rep< distributed >::name
std::string name() const
Definition: domain_indirect.h:272
rheolef::geo_abstract_rep< T, distributed >::ios_ige2dis_ige
virtual size_type ios_ige2dis_ige(size_type dim, size_type ios_ige) const =0
rheolef::reference_element::e
static const variant_type e
Definition: reference_element.h:76
rheolef::compute_bbox
void compute_bbox(const geo_base_rep< T, M > &omega, const geo_element &K, point_basic< T > &xmin, point_basic< T > &xmax)
Definition: geo_locate.cc:50
rheolef::index_set::insert
void insert(size_type dis_i)
Definition: index_set_header.icc:158
rheolef::reference_element::H
static const variant_type H
Definition: reference_element.h:81
rheolef::domain_indirect_basic< distributed >::map_dimension
size_type map_dimension() const
Definition: domain_indirect.h:695
rheolef::geo_base_rep::n_domain_indirect
size_type n_domain_indirect() const
Definition: geo.h:595
rheolef::domain_indirect_rep< distributed >::map_dimension
size_type map_dimension() const
Definition: domain_indirect.h:273
rheolef::reference_element::init_local_nnode_by_variant
static void init_local_nnode_by_variant(size_type order, std::array< size_type, reference_element::max_variant > &loc_nnod_by_variant)
Definition: reference_element.cc:206
check_macro
check_macro(expr1.have_homogeneous_space(Xh1), "dual(expr1,expr2); expr1 should have homogeneous space. HINT: use dual(interpolate(Xh, expr1),expr2)")
rheolef::distributor::comm
const communicator_type & comm() const
Definition: distributor.h:145
rheolef::domain_indirect_basic< distributed >::name
std::string name() const
Definition: domain_indirect.h:689
rheolef::geo_size::first_by_variant
distributor first_by_variant[reference_element::max_variant]
Definition: geo_size.h:66
rheolef::distributor::find_owner
size_type find_owner(size_type dis_i) const
find iproc associated to a global index dis_i: CPU=log(nproc)
Definition: distributor.cc:106
rheolef::geo_abstract_base_rep::dimension
virtual size_type dimension() const =0
rheolef::geo_element_indirect::index
size_type index() const
Definition: geo_element_indirect.h:51
rheolef::geo_abstract_base_rep::name
virtual std::string name() const =0
rheolef::geo_base_rep< T, distributed >
rheolef::geo_base_rep::_geo_element
std::array< hack_array< geo_element_hack, M >, reference_element::max_variant > _geo_element
Definition: geo.h:678
rheolef::reference_element::T
static const variant_type T
Definition: reference_element.h:79
rheolef::distributor
see the distributor page for the full documentation
Definition: distributor.h:62
mkgeo_ball.order
order
Definition: mkgeo_ball.sh:343
rheolef::geo_element::n_subgeo
size_type n_subgeo(size_type subgeo_dim) const
Definition: geo_element.h:212
rheolef::geo_element::size
size_type size() const
Definition: geo_element.h:168
rheolef::geo_element
see the geo_element page for the full documentation
Definition: geo_element.h:102
rheolef::geo_abstract_base_rep::node
virtual const node_type & node(size_type inod) const =0
rheolef::geo_rep
sequential mesh representation
Definition: geo.h:778
rheolef::index_set
Definition: index_set_header.icc:58
rheolef::geo_element_auto
Definition: geo_element.h:376
rheolef::distributor::decide
static const size_type decide
Definition: distributor.h:76
rheolef::geo_abstract_base_rep::geo_element_ownership
virtual const distributor & geo_element_ownership(size_type dim) const =0
rheolef::geo_element::variant
variant_type variant() const
Definition: geo_element.h:161
rheolef::domain_indirect_basic< distributed >::ioige_end
const_iterator_ioige ioige_end() const
Definition: domain_indirect.h:737
rheolef::reference_element::first_variant_by_dimension
static variant_type first_variant_by_dimension(size_type dim)
Definition: reference_element.h:148
rheolef::geo_element_indirect
Definition: geo_element_indirect.h:32
mkgeo_ball.variant
int variant
Definition: mkgeo_ball.sh:149
rheolef::geo_abstract_base_rep::coordinate_system
virtual coordinate_type coordinate_system() const =0
rheolef::geo_abstract_base_rep::get_piola_basis
virtual const basis_basic< T > & get_piola_basis() const =0
rheolef::geo_element::set_dis_ie
void set_dis_ie(size_type dis_ie)
Definition: geo_element.h:172
rheolef::domain_indirect_basic< distributed >
Definition: domain_indirect.h:559
rheolef::geo_base_rep::get_domain_indirect
const domain_indirect_basic< M > & get_domain_indirect(size_type i) const
Definition: geo.h:597
rheolef::geo_abstract_rep< T, distributed >::dis_ige2ios_dis_ige
virtual size_type dis_ige2ios_dis_ige(size_type dim, size_type dis_ige) const =0
rheolef::geo_element::dimension
size_type dimension() const
Definition: geo_element.h:167
assert_macro
#define assert_macro(ok_condition, message)
Definition: dis_macros.h:113
rheolef::domain_indirect_basic< distributed >::ioige_begin
const_iterator_ioige ioige_begin() const
Definition: domain_indirect.h:731
rheolef::geo_abstract_base_rep::dis_inod2dis_iv
virtual size_type dis_inod2dis_iv(size_type dis_inod) const =0
rheolef::geo_abstract_base_rep::get_geo_element
virtual const_reference get_geo_element(size_type dim, size_type ige) const =0
rheolef
This file is part of Rheolef.
Definition: compiler_eigen.h:37
rheolef::distributor::dis_size
size_type dis_size() const
global and local sizes
Definition: distributor.h:207
error_macro
#define error_macro(message)
Definition: dis_macros.h:49
rheolef::geo_element::parameter_type
Definition: geo_element.h:136
rheolef::domain_indirect_rep< distributed >::oige
const geo_element_indirect & oige(size_type ioige) const
Definition: domain_indirect.h:267
rheolef::geo_size::ownership_by_variant
distributor ownership_by_variant[reference_element::max_variant]
Definition: geo_size.h:64
rheolef::reference_element::max_variant
static const variant_type max_variant
Definition: reference_element.h:82
rheolef::disarray< size_type >
rheolef::geo_abstract_base_rep::sizes
virtual const geo_size & sizes() const =0
rheolef::reference_element::p
static const variant_type p
Definition: reference_element.h:75
rheolef::reference_element::q
static const variant_type q
Definition: reference_element.h:78
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::geo_abstract_rep< T, distributed >::geo_element_ios_ownership
virtual distributor geo_element_ios_ownership(size_type dim) const =0
rheolef::domain_indirect_rep< distributed >
Definition: domain_indirect.h:229
rheolef::reference_element::P
static const variant_type P
Definition: reference_element.h:80
rheolef::geo_abstract_rep< T, distributed >
Definition: geo.h:457
rheolef::reference_element::t
static const variant_type t
Definition: reference_element.h:77
rheolef::geo_rep< T, distributed >::size_type
base::size_type size_type
Definition: geo.h:934
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::domain_indirect_rep< distributed >::ioige2ini_dis_ioige
size_type ioige2ini_dis_ioige(size_type ioige) const
Definition: domain_indirect.h:280
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
mkgeo_ball.map_dim
map_dim
Definition: mkgeo_ball.sh:337
rheolef::domain_indirect_rep< distributed >::_ini_ioige2dis_ioige
disarray< size_type, distributed > _ini_ioige2dis_ioige
Definition: domain_indirect.h:294
rheolef::geo_element::face
size_type face(size_type i) const
Definition: geo_element.h:210
mkgeo_ball.dim
int dim
Definition: mkgeo_ball.sh:307
rheolef::domain_indirect_basic
the finite element boundary domain
Definition: domain_indirect.h:335
rheolef::geo_element::edge
size_type edge(size_type i) const
Definition: geo_element.h:209
rheolef::distributor::size
size_type size(size_type iproc) const
Definition: distributor.h:163
rheolef::geo_element::ios_dis_ie
size_type ios_dis_ie() const
Definition: geo_element.h:164
rheolef::geo_element::set_ios_dis_ie
void set_ios_dis_ie(size_type ios_dis_ie)
Definition: geo_element.h:173