Rheolef  7.2
an efficient C++ finite element environment
geo_domain_seq.cc
Go to the documentation of this file.
1 
22 #include "rheolef/geo_domain.h"
23 
24 namespace rheolef {
25 
26 // ------------------------------------------------------------------------
27 // build edge and face connectivity from domain
28 // ------------------------------------------------------------------------
29 template <class T>
30 void
32  const domain_indirect_rep<sequential>& indirect,
33  const geo_abstract_rep<T,sequential>& bgd_omega,
34  size_type sid_dim,
35  disarray<size_type,sequential>& bgd_isid2dom_dis_isid,
36  disarray<size_type,sequential>& dom_isid2bgd_isid,
37  disarray<size_type,sequential>& dom_isid2dom_ios_dis_isid,
38  size_type size_by_variant [reference_element::max_variant])
39 {
40  if (sid_dim != 0 && base::_gs._map_dimension <= sid_dim) return;
41  communicator comm = bgd_omega.geo_element_ownership(sid_dim).comm();
42  // ------------------------------------------------------------------------
43  // 1) side compact re-numbering
44  // ------------------------------------------------------------------------
45  // 1.1 loop on elements and mark used sides
46  disarray<size_type,sequential> bgd_isid_is_on_domain (bgd_omega.geo_element_ownership(sid_dim), 0); // logical, init to "false"
47  for (size_type ioige = 0, noige = indirect.size(); ioige < noige; ioige++) {
48  size_type ige = indirect.oige (ioige).index();
49  const geo_element& bgd_K = bgd_omega.get_geo_element (base::_gs._map_dimension, ige);
50  for (size_type loc_isid = 0, loc_nsid = bgd_K.n_subgeo(sid_dim); loc_isid < loc_nsid; loc_isid++) {
51  size_type bgd_dis_isid = 0; // TODO: = bgd_K.subgeo (sid_dim, loc_isid);
52  switch (sid_dim) {
53  case 0: {
54  size_type bgd_dis_inod = bgd_K[loc_isid];
55  bgd_dis_isid = bgd_omega.dis_inod2dis_iv (bgd_dis_inod);
56  break;
57  }
58  case 1: bgd_dis_isid = bgd_K.edge(loc_isid); break;
59  case 2: bgd_dis_isid = bgd_K.face(loc_isid); break;
60  default: error_macro ("domain: unexpected side dimension " << sid_dim);
61  }
62  bgd_isid_is_on_domain[bgd_dis_isid] += 1;
63  }
64  }
65  // 1.2 counting & distribution for dom_isid
66  size_type dom_nsid = 0;
67  for (size_type bgd_isid = 0, bgd_nsid = bgd_omega.geo_element_ownership(sid_dim).size(); bgd_isid < bgd_nsid; bgd_isid++) {
68  if (bgd_isid_is_on_domain[bgd_isid] != 0) dom_nsid++ ;
69  }
70  // 1.4 numbering dom_isid & permutation: bgd_isid <--> dom_isid
73  size_by_variant [variant] = 0;
74  }
75  distributor dom_isid_ownership (distributor::decide, comm, dom_nsid);
76  bgd_isid2dom_dis_isid.resize (bgd_omega.geo_element_ownership(sid_dim), std::numeric_limits<size_type>::max());
77  dom_isid2bgd_isid.resize (dom_isid_ownership, std::numeric_limits<size_type>::max());
78  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++) {
79  if (bgd_isid_is_on_domain[bgd_isid] == 0) continue;
80  size_type dom_dis_isid = dom_isid; // sequential case !
81  bgd_isid2dom_dis_isid [bgd_isid] = dom_dis_isid;
82  dom_isid2bgd_isid [dom_isid] = bgd_isid;
83  const geo_element& bgd_S = bgd_omega.get_geo_element(sid_dim,bgd_isid);
84  size_by_variant [bgd_S.variant()]++;
85  dom_isid++;
86  }
87  // ------------------------------------------------------------------------
88  // 2. compute sizes and resize _geo_element[variant]
89  // ------------------------------------------------------------------------
90  size_type dom_nge = 0;
91  size_type dom_dis_nge = 0;
94 
95  distributor dom_igev_ownership (distributor::decide, comm, size_by_variant [variant]);
97  base::_geo_element[variant].resize (dom_igev_ownership, param);
98  base::_gs.ownership_by_variant[variant] = dom_igev_ownership;
99  dom_nge += dom_igev_ownership.size();
100  dom_dis_nge += dom_igev_ownership.dis_size();
101  }
102  base::_gs.ownership_by_dimension[sid_dim] = distributor (dom_dis_nge, comm, dom_nge);
103 }
104 template <class T>
105 void
107  const domain_indirect_rep<sequential>& indirect,
108  const geo_abstract_rep<T,sequential>& bgd_omega,
109  disarray<size_type,sequential>& bgd_iv2dom_dis_iv,
110  size_type sid_dim,
111  disarray<size_type,sequential>& bgd_isid2dom_dis_isid,
112  disarray<size_type,sequential>& dom_isid2bgd_isid,
113  disarray<size_type,sequential>& dom_isid2dom_ios_dis_isid,
114  size_type size_by_variant [reference_element::max_variant])
115 {
116  if (sid_dim != 0 && base::_gs._map_dimension <= sid_dim) return;
117  communicator comm = bgd_omega.geo_element_ownership(sid_dim).comm();
118  // ------------------------------------------------------------------------
119  // 2) set _geo_element[variant] and S.set_ios_dis_ie
120  // also reperate external vertices
121  // ------------------------------------------------------------------------
122  distributor dom_isid_ownership = dom_isid2bgd_isid.ownership();
123  for (size_type dom_isid = 0, dom_nsid = dom_isid_ownership.size(); dom_isid < dom_nsid; dom_isid++) {
124  size_type bgd_isid = dom_isid2bgd_isid [dom_isid];
125  size_type dom_dis_isid = dom_isid;
126  size_type dom_ios_dis_isid = dom_dis_isid;
127  const geo_element& bgd_S = bgd_omega.get_geo_element(sid_dim,bgd_isid);
128  geo_element& dom_S = get_geo_element(sid_dim,dom_isid);
129  dom_S = bgd_S;
130  dom_S.set_dis_ie (dom_dis_isid);
131  dom_S.set_ios_dis_ie (dom_ios_dis_isid);
132  // set S face & edge : index, orientation & rotation will be set by propagate_numbering later
133  }
134  // ------------------------------------------------------------------------
135  // 3) propagate new vertex numbering in all `dom_S' new sides
136  // ------------------------------------------------------------------------
137  if (sid_dim > 0) {
138  for (size_type dom_isid = 0, dom_nsid = dom_isid_ownership.size(); dom_isid < dom_nsid; dom_isid++) {
139  geo_element& dom_S = get_geo_element(sid_dim,dom_isid);
140  if (dom_S.dimension() == 0) continue;
141  for (size_type iloc = 0, nloc = dom_S.size(); iloc < nloc; iloc++) {
142  size_type bgd_dis_inod = dom_S[iloc];
143  size_type bgd_dis_iv = bgd_omega.dis_inod2dis_iv (bgd_dis_inod);
144  size_type dom_dis_iv = bgd_iv2dom_dis_iv.dis_at (bgd_dis_iv);
145  size_type dom_dis_inod = base::dis_iv2dis_inod (dom_dis_iv);
146  dom_S[iloc] = dom_dis_inod;
147  }
148  }
149  }
150 }
151 // ------------------------------------------------------------------------
152 // geo construtor: build from domain
153 // ------------------------------------------------------------------------
154 template <class T>
155 void
157  const domain_indirect_rep<sequential>& indirect,
158  const geo_abstract_rep<T,sequential>& bgd_omega,
159  std::map<size_type,size_type>& bgd_ie2dom_ie,
160  std::map<size_type,size_type>& dis_bgd_ie2dis_dom_ie)
161 {
162 trace_macro ("build("<<bgd_omega.name()<<","<<indirect.name()<<")...");
163  base::_name = bgd_omega.name() + "[" + indirect.name() + "]";
164  base::_version = 4;
165  base::_sys_coord = bgd_omega.coordinate_system();
166  base::_dimension = bgd_omega.dimension();
167  base::_piola_basis = bgd_omega.get_piola_basis();
168  base::_gs._map_dimension = indirect.map_dimension();
169  base::_have_connectivity = 1;
170  size_type map_dim = base::_gs._map_dimension;
171  size_type size_by_variant [reference_element::max_variant];
172  std::fill (size_by_variant, size_by_variant+reference_element::max_variant, 0);
173  // ----------------------------------------------------
174  // 1) _geo_element[0]: compact re-numbering
175  // ----------------------------------------------------
176  disarray<size_type,sequential> dom_isid2dom_ios_dis_isid [4];
177  std::array<disarray<size_type,sequential>,4> bgd_ige2dom_dis_ige;
178  std::array<disarray<size_type,sequential>,4> dom_ige2bgd_ige;
179  domain_set_side_part1 (indirect, bgd_omega, 0,
180  bgd_ige2dom_dis_ige[0], dom_ige2bgd_ige[0],
181  dom_isid2dom_ios_dis_isid [0], size_by_variant);
182  domain_set_side_part2 (indirect, bgd_omega, bgd_ige2dom_dis_ige[0], 0,
183  bgd_ige2dom_dis_ige[0], dom_ige2bgd_ige[0],
184  dom_isid2dom_ios_dis_isid [0], size_by_variant);
185  // ------------------------------------------------------------------------
186  // 2) count elements by variants
187  // ------------------------------------------------------------------------
190  size_by_variant [variant] = 0;
191  }
192  for (size_type ioige = 0, noige = indirect.size(); ioige < noige; ioige++) {
193  size_type ige = indirect.oige (ioige).index();
194  bgd_ie2dom_ie [ige] = ioige;
195  const geo_element& bgd_K = bgd_omega.get_geo_element (map_dim, ige);
196  size_by_variant [bgd_K.variant()]++;
197  }
198  // -----------------------------------------------
199  // 3) _geo_element[map_dim]: compact vertex numbering
200  // -----------------------------------------------
201  size_type dis_nge = 0;
202  size_type nge = 0;
205  distributor dom_igev_ownership (distributor::decide, base::comm(), size_by_variant [variant]);
207  base::_geo_element[variant].resize (dom_igev_ownership, param);
208  base::_gs.ownership_by_variant [variant] = dom_igev_ownership;
209  dis_nge += dom_igev_ownership.dis_size();
210  nge += dom_igev_ownership.size();
211  }
212  base::_gs.ownership_by_dimension [map_dim] = distributor (dis_nge, base::comm(), nge);
213  // ------------------------------------------------------------------------
214  // 4) count all geo_element[] and set base::_gs.ownership_by_variant[]
215  // => then can determine the node count (order > 1)
216  // ------------------------------------------------------------------------
217  for (size_type sid_dim = 1; sid_dim < base::_gs._map_dimension; sid_dim++) {
218  domain_set_side_part1 (indirect, bgd_omega, sid_dim,
219  bgd_ige2dom_dis_ige[sid_dim], dom_ige2bgd_ige[sid_dim],
220  dom_isid2dom_ios_dis_isid [sid_dim], size_by_variant);
221  }
222  // ------------------------------------------------------------------------
223  // 5) count nodes & set base::_gs.node_ownership
224  // => then dis_iv2dis_inod works (used at step 6)
225  // ------------------------------------------------------------------------
226  std::array<size_type,reference_element::max_variant> loc_nnod_by_variant ;
228  size_type nnod = 0;
229  for (size_type variant = 0;
230  variant < reference_element::last_variant_by_dimension(base::map_dimension());
231  variant++) {
232  nnod += base::_gs.ownership_by_variant [variant].size() * loc_nnod_by_variant [variant];
233  }
234  distributor dom_node_ownership (distributor::decide, base::comm(), nnod);
235  base::_gs.node_ownership = dom_node_ownership;
236  // ------------------------------------------------------------------------
237  // 6) set _geo_element[] values
238  // ------------------------------------------------------------------------
239  size_type first_dis_ioige = indirect.ownership().first_index();
240  for (size_type ioige = 0, noige = indirect.size(); ioige < noige; ioige++) {
241  size_type dis_ioige = first_dis_ioige + ioige;
242  size_type ige = indirect.oige (ioige).index();
243  geo_element& dom_K = get_geo_element (map_dim, ioige);
244  const geo_element& bgd_K = bgd_omega.get_geo_element (map_dim, ige);
245  dom_K = bgd_K;
246  dom_K.set_dis_ie (dis_ioige);
247  size_type ini_dis_ioige = ioige;
248  dom_K.set_ios_dis_ie (ini_dis_ioige);
249  // set K face & edge : index, orientation & rotation will be set by propagate_numbering later
250  for (size_type iloc = 0, nloc = dom_K.size(); iloc < nloc; iloc++) {
251  size_type bgd_dis_inod = bgd_K[iloc];
252  size_type bgd_dis_iv = bgd_omega.dis_inod2dis_iv (bgd_dis_inod);
253  size_type dom_dis_iv = bgd_ige2dom_dis_ige[0].dis_at (bgd_dis_iv);
254  size_type dom_dis_inod = base::dis_iv2dis_inod (dom_dis_iv);
255  dom_K[iloc] = dom_dis_inod;
256  }
257  }
258  // reset also dom_ige2bgd_ige[map_dim] : used by _node[] for order > 1 geometries
259  if (base::_gs._map_dimension > 0) {
260  dom_ige2bgd_ige [base::_gs._map_dimension].resize(indirect.ownership());
261  bgd_ige2dom_dis_ige[base::_gs._map_dimension].resize(bgd_omega.geo_element_ownership(base::_gs._map_dimension), std::numeric_limits<size_type>::max());
262  for (size_type ioige = 0, noige = indirect.size(); ioige < noige; ioige++) {
263  size_type bgd_ige = indirect.oige (ioige).index();
264  dom_ige2bgd_ige [base::_gs._map_dimension] [ioige] = bgd_ige;
265  bgd_ige2dom_dis_ige[base::_gs._map_dimension] [bgd_ige] = ioige; // sequential!
266  }
267  }
268  // ------------------------------------------------------------------------
269  // 7) _geo_element[1&2]: compact renumbering
270  // _ios_ige2dis_ige[1&2]: idem
271  // ------------------------------------------------------------------------
272  for (size_type sid_dim = 1; sid_dim < base::_gs._map_dimension; sid_dim++) {
273  domain_set_side_part2 (indirect, bgd_omega, bgd_ige2dom_dis_ige[0], sid_dim,
274  bgd_ige2dom_dis_ige[sid_dim], dom_ige2bgd_ige[sid_dim],
275  dom_isid2dom_ios_dis_isid [sid_dim], size_by_variant);
276  }
277  // ------------------------------------------------------------------------
278  // 11) raw copy _node
279  // TODO: use shallow copy & indirection instead (memory usage reduction):
280  // node(inod) { return bgd_omega.node(inod2bgd_inod(inod))}
281  // ------------------------------------------------------------------------
282  base::_node.resize (dom_node_ownership);
283  disarray<size_type,sequential> dom_inod2bgd_inod (dom_node_ownership, std::numeric_limits<size_type>::max());
284  size_type dom_inod = 0;
285  size_type first_bgd_inod_v = 0;
286  for (size_type dim = 0; dim <= base::_gs._map_dimension; dim++) {
287  size_type dom_ige = 0;
290  variant++) {
291  if (loc_nnod_by_variant [variant] == 0) continue;
292  size_type first_bgd_v = bgd_omega.sizes().first_by_variant [variant].size();
293  for (size_type dom_igev = 0, dom_ngev = base::_geo_element [variant].size(); dom_igev < dom_ngev; dom_igev++, dom_ige++) {
294  const geo_element& K = base::_geo_element [variant][dom_igev];
295  size_type bgd_ige = dom_ige2bgd_ige [dim][dom_ige];
296  assert_macro (bgd_ige >= first_bgd_v, "invalid index");
297  size_type bgd_igev = bgd_ige - first_bgd_v;
298  for (size_type loc_inod = 0, loc_nnod = loc_nnod_by_variant [variant]; loc_inod < loc_nnod; loc_inod++, dom_inod++) {
299  size_type bgd_inod = first_bgd_inod_v + bgd_igev * loc_nnod_by_variant [variant] + loc_inod;
300  dom_inod2bgd_inod [dom_inod] = bgd_inod;
301  base::_node [dom_inod] = bgd_omega.node (bgd_inod);
302  }
303  }
304  first_bgd_inod_v += bgd_omega.sizes().ownership_by_variant [variant].size() * loc_nnod_by_variant [variant];
305  }
306  }
308  // ------------------------------------------------------------------------
309  // 12) propagate new edge & face numbering to all `dom_S' elements
310  // ------------------------------------------------------------------------
311  for (size_type dim = 1; dim < base::_gs._map_dimension; dim++) {
312  set_element_side_index (dim);
313  }
314  // ------------------------------------------------------------------------
315  // 13) reset vertex P[0] value (useful for band zero vertex domain)
316  // ------------------------------------------------------------------------
317  for (size_type dom_iv = 0, dom_nv = base::_geo_element[reference_element::p].size(); dom_iv < dom_nv; dom_iv++) {
318  geo_element& P = base::_geo_element[reference_element::p][dom_iv];
319  size_type dom_inod = base::dis_iv2dis_inod (dom_iv);
320  P[0] = dom_inod;
321  }
322  // ------------------------------------------------------------------------
323  // 14) domains : intersection with the current domain
324  // ------------------------------------------------------------------------
325  const geo_base_rep<T,sequential> *ptr = dynamic_cast<const geo_base_rep<T,sequential>*> (&bgd_omega);
326  check_macro (ptr != 0, "cannot build domains on \""<<base::_name<<"\"");
327  const geo_base_rep<T,sequential>& bgd_omega2 = *ptr;
328  for (size_type idom = 0, ndom = bgd_omega2.n_domain_indirect(); idom < ndom; ++idom) {
329  const domain_indirect_basic<sequential>& bgd_dom = bgd_omega2.get_domain_indirect(idom);
330 trace_macro ("build("<<bgd_omega.name()<<"["<<indirect.name()<<"] with " <<bgd_dom.name()<<"...");
331  if (bgd_dom.name() == indirect.name()) continue; // myself
332  size_type dom_map_dim = bgd_dom.map_dimension();
333  if (dom_map_dim > base::_gs._map_dimension) continue; // dom has upper dimension
334  std::vector<size_type> ie_list;
335  for (domain_indirect_basic<sequential>::const_iterator_ioige iter = bgd_dom.ioige_begin(), last = bgd_dom.ioige_end(); iter != last; ++iter) {
336  const geo_element_indirect& ioige = *iter;
337  size_type bgd_ige = ioige.index();
338  size_type dom_dis_ige = bgd_ige2dom_dis_ige [dom_map_dim][bgd_ige];
339  if (dom_dis_ige == std::numeric_limits<size_type>::max()) continue; // side do not belongs to dom
340  size_type dom_ige = dom_dis_ige; // sequential !
341  ie_list.push_back(dom_ige);
342  }
343  if (ie_list.size() == 0) continue; // empty intersection
344  domain_indirect_basic<sequential> dom (*this, bgd_dom.name(), bgd_dom.map_dimension(), communicator(), ie_list);
345  base::_domains.push_back (dom);
346 trace_macro ("geo("<<base::name()<<") += domain("<<dom.name()<<") map_d="<<dom.map_dimension());
347  }
348 }
349 // ----------------------------------------------------------------------------
350 // instanciation in library
351 // ----------------------------------------------------------------------------
352 template class geo_rep<Float,sequential>;
353 
354 } // namespace rheolef
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 dis_size() const
global and local sizes
Definition: distributor.h:214
size_type size(size_type iproc) const
Definition: distributor.h:170
static const size_type decide
Definition: distributor.h:83
const communicator_type & comm() const
Definition: distributor.h:152
const_iterator_ioige ioige_begin() const
the finite element boundary domain
const geo_element_indirect & oige(size_type ioige) const
virtual const geo_size & sizes() const =0
virtual const node_type & node(size_type inod) const =0
virtual std::string name() const =0
geo_element_hack::size_type size_type
Definition: geo.h:260
virtual const_reference get_geo_element(size_type dim, size_type ige) const =0
virtual const distributor & geo_element_ownership(size_type dim) 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 const basis_basic< T > & get_piola_basis() const =0
size_type n_domain_indirect() const
Definition: geo.h:595
const domain_indirect_basic< M > & get_domain_indirect(size_type i) const
Definition: geo.h:597
see the geo_element page for the full documentation
Definition: geo_element.h:102
size_type edge(size_type i) const
Definition: geo_element.h:209
size_type face(size_type i) const
Definition: geo_element.h:210
size_type size() const
Definition: geo_element.h:168
void set_ios_dis_ie(size_type ios_dis_ie)
Definition: geo_element.h:173
size_type dimension() const
Definition: geo_element.h:167
variant_type variant() const
Definition: geo_element.h:161
size_type n_subgeo(size_type subgeo_dim) const
Definition: geo_element.h:212
void set_dis_ie(size_type dis_ie)
Definition: geo_element.h:172
sequential mesh representation
Definition: geo.h:778
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
#define trace_macro(message)
Definition: dis_macros.h:111
#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)")
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)
Definition: geo_locate.cc:50
distributor ownership_by_variant[reference_element::max_variant]
Definition: geo_size.h:64
distributor first_by_variant[reference_element::max_variant]
Definition: geo_size.h:66