Rheolef  7.2
an efficient C++ finite element environment
 
Loading...
Searching...
No Matches
geo_domain_seq.cc
Go to the documentation of this file.
1
21
22#include "rheolef/geo_domain.h"
23
24namespace rheolef {
25
26// ------------------------------------------------------------------------
27// build edge and face connectivity from domain
28// ------------------------------------------------------------------------
29template <class T>
30void
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,
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
72 variant < reference_element:: last_variant_by_dimension(sid_dim); variant++) {
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;
93 variant < reference_element:: last_variant_by_dimension(sid_dim); variant++) {
94
95 distributor dom_igev_ownership (distributor::decide, comm, size_by_variant [variant]);
96 geo_element::parameter_type param (variant, 1);
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}
104template <class T>
105void
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,
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// ------------------------------------------------------------------------
154template <class T>
155void
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{
162trace_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;
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 // ------------------------------------------------------------------------
189 variant < reference_element:: last_variant_by_dimension(map_dim); variant++) {
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;
204 variant < reference_element:: last_variant_by_dimension(map_dim); variant++) {
205 distributor dom_igev_ownership (distributor::decide, base::comm(), size_by_variant [variant]);
206 geo_element::parameter_type param (variant, 1);
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 ;
227 reference_element::init_local_nnode_by_variant (base::order(), 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 }
307 base::compute_bbox();
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);
330trace_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);
346trace_macro ("geo("<<base::name()<<") += domain("<<dom.name()<<") map_d="<<dom.map_dimension());
347 }
348}
349// ----------------------------------------------------------------------------
350// instanciation in library
351// ----------------------------------------------------------------------------
352template 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
size_type size(size_type iproc) const
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
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
variant_type variant() const
size_type n_subgeo(size_type subgeo_dim) const
void set_dis_ie(size_type dis_ie)
base::size_type size_type
Definition geo.h:791
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)")
This file is part of Rheolef.