Rheolef  7.2
an efficient C++ finite element environment
 
Loading...
Searching...
No Matches
geo_seq_upgrade.cc
Go to the documentation of this file.
1
21#include "rheolef/geo.h"
22#include "rheolef/geo_domain.h"
23#include "rheolef/rheostream.h"
24#include "rheolef/iorheo.h"
25#include "rheolef/space_numbering.h"
26
27#include "geo_header.h"
28
29namespace rheolef {
30
31// ----------------------------------------------------------------------------
32// build connectivity: for input files such as gmsh that do not have
33// a global numbering of edges and faces
34// ----------------------------------------------------------------------------
35
36template <class T>
37void
40{
41 typedef geo_element_auto<> element_t;
42 typedef disarray<element_t,sequential> disarray_t;
43 typedef typename disarray_t::const_iterator const_iterator_by_variant_tmp;
44 size_type map_dim = geo_base_rep<T,sequential>::_gs._map_dimension;
45 if (map_dim <= 1) return;
46 build_connectivity_sides (map_dim-1, tmp_geo_element);
47 build_connectivity_sides (map_dim-2, tmp_geo_element);
48}
49template <class T>
50void
52 size_type side_dim,
54{
55 typedef geo_element_auto<> element_t;
56 typedef disarray<element_t,sequential> disarray_t;
57 typedef typename disarray_t::const_iterator const_iterator_by_variant_tmp;
58 if (side_dim == 0) return;
59 // ------------------------------------------------------------------------
60 // 1) count sides
61 // ------------------------------------------------------------------------
62 // difficulte : on ne connait pas la taille de la table des aretes et faces par variante...
63 // first pass: count by variant and use a nsid global counter as pseudo side index
64 // this index is not usable since sides indexes may be ordered by increasing variants
65 // and here, during the first pass, sides appears by unordered variant
67 std::fill (size_by_variant, size_by_variant+reference_element::max_variant, 0);
68 index_set empty_set; // TODO: add a global allocator to empty_set
69 size_type map_dim = geo_base_rep<T,sequential>::_gs._map_dimension;
70 size_type nsid = 0;
71 {
72 disarray<index_set,sequential> ball_S (geo_element_ownership(0), empty_set);
74 variant < reference_element:: last_variant_by_dimension(map_dim); variant++) {
75 for (const_iterator_by_variant_tmp iter = tmp_geo_element [variant].begin(),
76 last = tmp_geo_element [variant].end();
77 iter != last; iter++) {
78
79 const geo_element& K = *iter;
80 for (size_type loc_isid = 0, loc_nsid = K.n_subgeo(side_dim); loc_isid < loc_nsid; loc_isid++) {
81 size_type S_variant = reference_element::variant (K.subgeo_size (side_dim, loc_isid), side_dim);
82 size_type loc_jv0 = K.subgeo_local_vertex (side_dim, loc_isid, 0);
83 size_type jv0 = K[loc_jv0];
84 index_set isid_set = ball_S [jv0]; // copy: will be intersected
85 for (size_type sid_jloc = 1, sid_nloc = K.subgeo_size (side_dim, loc_isid); sid_jloc < sid_nloc; sid_jloc++) {
86 size_type loc_jv = K.subgeo_local_vertex (side_dim, loc_isid, sid_jloc);
87 size_type jv = K[loc_jv];
88 const index_set& ball_jv = ball_S [jv];
89 isid_set.inplace_intersection (ball_jv);
90 }
91 check_macro (isid_set.size() <= 1, "connectivity problem");
92 if (isid_set.size() == 1) {
93 continue; // side already exists
94 }
95 // then insert it in ball_S
96 for (size_type sid_jloc = 0, sid_nloc = K.subgeo_size (side_dim, loc_isid); sid_jloc < sid_nloc; sid_jloc++) {
97 size_type loc_jv = K.subgeo_local_vertex (side_dim, loc_isid, sid_jloc);
98 size_type jv = K[loc_jv];
99 ball_S[jv] += nsid;
100 }
101 size_by_variant [S_variant]++;
102 nsid++;
103 }
104 }
105 }
106 }
107 // ------------------------------------------------------------------------
108 // 2) allocate tmp_geo_element[]
109 // ------------------------------------------------------------------------
110 size_type first_by_variant [reference_element::max_variant+1];
111 std::fill (first_by_variant, first_by_variant+reference_element::max_variant+1, 0);
112 for (size_type side_variant = reference_element::first_variant_by_dimension(side_dim);
113 side_variant < reference_element:: last_variant_by_dimension(side_dim); side_variant++) {
114 geo_element_auto<> init_val (side_variant, base::order());
115 tmp_geo_element [side_variant].resize (size_by_variant [side_variant], init_val);
116 base::_gs.ownership_by_variant [side_variant] = tmp_geo_element [side_variant].ownership();
117 first_by_variant [side_variant+1] = first_by_variant [side_variant] + size_by_variant [side_variant];
118 }
119 // ownership_by_dimension[side_dim]: used by connectivity
120 base::_gs.ownership_by_dimension [side_dim] = distributor (nsid, base::comm(), nsid);
121
122 // ------------------------------------------------------------------------
123 // 3) set numbering (ordered by variant) & create sides
124 // ------------------------------------------------------------------------
126 {
127 size_type count_by_variant [reference_element::max_variant];
128 std::fill (count_by_variant, count_by_variant+reference_element::max_variant, 0);
129 disarray<index_set,sequential> ball_S (geo_element_ownership(0), empty_set);
130 nsid = 0;
132 variant < reference_element:: last_variant_by_dimension(map_dim); variant++) {
133 for (const_iterator_by_variant_tmp iter = tmp_geo_element [variant].begin(),
134 last = tmp_geo_element [variant].end();
135 iter != last; iter++) {
136
137 const geo_element& K = *iter;
138 for (size_type loc_isid = 0, loc_nsid = K.n_subgeo(side_dim); loc_isid < loc_nsid; loc_isid++) {
139 size_type S_variant = reference_element::variant (K.subgeo_size (side_dim, loc_isid), side_dim);
140 size_type loc_jv0 = K.subgeo_local_vertex (side_dim, loc_isid, 0);
141 size_type jv0 = K[loc_jv0];
142 index_set isid_set = ball_S [jv0]; // copy: will be intersected
143 for (size_type sid_jloc = 1, sid_nloc = K.subgeo_size (side_dim, loc_isid); sid_jloc < sid_nloc; sid_jloc++) {
144 size_type loc_jv = K.subgeo_local_vertex (side_dim, loc_isid, sid_jloc);
145 size_type jv = K[loc_jv];
146 const index_set& ball_jv = ball_S [jv];
147 isid_set.inplace_intersection (ball_jv);
148 }
149 check_macro (isid_set.size() <= 1, "connectivity problem");
150 if (isid_set.size() == 1) {
151 continue; // side already exists
152 }
153 // side does not exists yet: then insert it in ball_S
154 for (size_type sid_jloc = 0, sid_nloc = K.subgeo_size (side_dim, loc_isid); sid_jloc < sid_nloc; sid_jloc++) {
155 size_type loc_jv = K.subgeo_local_vertex (side_dim, loc_isid, sid_jloc);
156 size_type jv = K[loc_jv];
157 ball_S[jv] += nsid;
158 }
159 size_type isid_by_variant = count_by_variant [S_variant];
160 count_by_variant [S_variant]++;
161 nsid++;
162 // create S
163 S.reset (S_variant, base::order());
164 // set also indexes
165 size_type ige = first_by_variant[S_variant] + isid_by_variant;
166 S.set_ios_dis_ie (ige);
167 S.set_dis_ie (ige);
168 // copy all nodes from K.subgeo into S that goes to tmp_gheo_element (will be conserved)
169 size_type S_n_node = K.subgeo_n_node (side_dim, loc_isid); // TODO: use ref_elem
170 for (size_type sid_jlocnod = 0, sid_nlocnod = S_n_node; sid_jlocnod < sid_nlocnod; sid_jlocnod++) {
171 size_type loc_jnod = K.subgeo_local_node (side_dim, loc_isid, sid_jlocnod); // TODO: use ref_elem
172 size_type jnod = K.node(loc_jnod);
173 S[sid_jlocnod] = jnod;
174 }
175 // then copy the side S in the tmp mesh:
176 tmp_geo_element [S_variant][isid_by_variant] = S;
177 }
178 }
179 }
180 }
181}
182// =========================================================================
183// get upgrade
184// =========================================================================
186template <class T>
189 idiststream& ips,
190 const geo_header& hdr)
191{
192 using namespace std;
193 check_macro (ips.good(), "bad input stream for geo.");
194 istream& is = ips.is();
195 // ------------------------------------------------------------------------
196 // 1) get nodes
197 // ------------------------------------------------------------------------
198 size_type nnod = hdr.dis_size_by_dimension [0];
200 if (hdr.dimension > 0) {
201 node.get_values (ips, _point_get<T>(hdr.dimension));
202 check_macro (ips.good(), "bad input stream for geo.");
203 }
204 // ------------------------------------------------------------------------
205 // 2) get elements
206 // ------------------------------------------------------------------------
207 typedef geo_element_auto<> element_t;
208 typedef disarray<element_t, sequential> disarray_t;
209 std::array<disarray_t, reference_element::max_variant> tmp_geo_element;
210 if (hdr.map_dimension > 0) {
213 geo_element_auto<> init_val (variant, hdr.order);
214 tmp_geo_element [variant] = disarray_t (hdr.dis_size_by_variant [variant], init_val);
215 tmp_geo_element [variant].get_values (ips);
216 check_macro (ips.good(), "bad input stream for geo.");
217 }
218 }
219 // ------------------------------------------------------------------------
220 // 3) build & upgrade
221 // ------------------------------------------------------------------------
222 bool do_upgrade = iorheo::getupgrade(is);
223 build_from_data (hdr, node, tmp_geo_element, do_upgrade);
224 // ------------------------------------------------------------------------
225 // 4) get domain, until end-of-file
226 // ------------------------------------------------------------------------
227 vector<index_set> ball [4];
229 while (dom.get (ips, *this, ball)) {
230 base::_domains.push_back (dom);
231 }
232 return ips;
233}
234template <class T>
235void
237 const geo_header& hdr,
240 tmp_geo_element,
241 bool do_upgrade)
242{
243 using namespace std;
244 typedef geo_element_auto<> element_t;
245 typedef disarray<element_t, sequential> disarray_t;
246 typedef typename disarray_t::const_iterator const_iterator_by_variant_tmp;
247 typedef typename disarray_t::iterator iterator_by_variant_tmp;
248 // ------------------------------------------------------------------------
249 // 1) get header
250 // ------------------------------------------------------------------------
251 base::_have_connectivity = false;
252 base::_version = 4;
253 base::_name = "unnamed";
254 base::_dimension = hdr.dimension;
255 base::_gs._map_dimension = hdr.map_dimension;
256 base::_sys_coord = hdr.sys_coord;
257 base::_piola_basis.reset_family_index (hdr.order);
258 size_type nnod = hdr.dis_size_by_dimension [0];
259 size_type n_edg = hdr.dis_size_by_dimension [1];
260 size_type n_fac = hdr.dis_size_by_dimension [2];
261 size_type n_elt = hdr.dis_size_by_dimension [base::_gs._map_dimension];
262 // ------------------------------------------------------------------------
263 // 2) get coordinates
264 // ------------------------------------------------------------------------
265 base::_node = node;
266 base::_gs.node_ownership = base::_node.ownership();
267 // ------------------------------------------------------------------------
268 // 3) bounding box: _xmin, _xmax
269 // ------------------------------------------------------------------------
270 base::compute_bbox();
271 // ------------------------------------------------------------------------
272 // 4) get elements
273 // ------------------------------------------------------------------------
274 if (base::_gs._map_dimension > 0) {
275 for (size_type variant = reference_element::first_variant_by_dimension(base::_gs._map_dimension);
276 variant < reference_element:: last_variant_by_dimension(base::_gs._map_dimension); variant++) {
277 base::_gs.ownership_by_variant [variant] = tmp_geo_element [variant].ownership();
278 }
279 base::_gs.ownership_by_dimension [base::_gs._map_dimension] = distributor (n_elt, base::comm(), n_elt);
280 }
281 // check that nodes are numbered by increasing node_subgeo_dim
282 {
283 std::vector<size_type> node_subgeo_dim (nnod, size_type(-1));
284 size_type prev_variant = 0;
285 for (size_type variant = reference_element::first_variant_by_dimension(base::_gs._map_dimension);
286 variant < reference_element:: last_variant_by_dimension(base::_gs._map_dimension); variant++) {
287 for (const_iterator_by_variant_tmp iter = tmp_geo_element [variant].begin(),
288 last = tmp_geo_element [variant].end();
289 iter != last; iter++) {
290
291 const geo_element& K = *iter;
292 check_macro (prev_variant <= K.variant(), "elements should be numbered by increasing variants (petqTPH)");
293 prev_variant = K.variant();
294 for (size_type d = 0; d <= base::_gs._map_dimension; d++) {
295 for (size_type loc_inod = K.first_inod(d), loc_nnod = K.last_inod(d); loc_inod < loc_nnod; loc_inod++) {
296 check_macro (K[loc_inod] < nnod, "inod K["<<loc_inod<<"]="<<K[loc_inod]
297 <<" is out of range [0:"<<nnod<<"[; K.dis_ie="<<K.dis_ie());
298 node_subgeo_dim [K[loc_inod]] = d;
299 }
300 }
301 }
302 }
303 size_type prev_node_dim = 0;
304 size_type i_curr_node = 0;
305 for (typename std::vector<size_type>::const_iterator iter = node_subgeo_dim.begin(), last = node_subgeo_dim.end();
306 iter != last; iter++, i_curr_node++) {
307 check_macro (prev_node_dim <= *iter, "nodes should be numbered by increasing subgeo dimension (prev="
308 << prev_node_dim << " > subgeo_dim(node("<<i_curr_node<<"))=" << *iter << ")");
309 prev_node_dim = *iter;
310 }
311 }
312 // compute n_vert (n_vert < nnod when order > 1) and set element indexes (K.dis_ie & K.ios_dis_ie)
313 size_type n_vert = 0;
314 if (base::_gs._map_dimension == 0) {
315 n_vert = nnod;
316 } else {
317 vector<size_t> is_vertex (nnod, 0);
318 size_type ie = 0;
319 for (size_type variant = reference_element::first_variant_by_dimension(base::_gs._map_dimension);
320 variant < reference_element:: last_variant_by_dimension(base::_gs._map_dimension); variant++) {
321 for (iterator_by_variant_tmp iter = tmp_geo_element [variant].begin(),
322 last = tmp_geo_element [variant].end();
323 iter != last; iter++, ie++) {
324
325 geo_element& K = *iter;
326 K.set_ios_dis_ie (ie);
327 K.set_dis_ie (ie);
328 if (base::order() > 1) {
329 for (size_type iloc = 0, nloc = K.size(); iloc < nloc; iloc++) {
330 is_vertex [K[iloc]] = 1;
331 }
332 }
333 }
334 }
335 if (base::order() == 1) {
336 n_vert = nnod;
337 } else {
338 n_vert = accumulate (is_vertex.begin(), is_vertex.end(), 0);
339 }
340 }
341 // ------------------------------------------------------------------------
342 // 5) create vertex-element (0d elements)
343 // ------------------------------------------------------------------------
344 {
345 geo_element_auto<> init_val (reference_element::p, base::order());
346 tmp_geo_element [reference_element::p].resize (n_vert, init_val);
347 size_type first_iv = base::_node.ownership().first_index();
348 size_type iv = 0;
349 for (iterator_by_variant_tmp iter = tmp_geo_element [reference_element::p].begin(),
350 last = tmp_geo_element [reference_element::p].end();
351 iter != last; iter++, iv++) {
352 geo_element& P = *iter;
353 P[0] = first_iv + iv;
354 P.set_dis_ie (first_iv + iv); // TODO: P[0] & dis_ie redundant for `p'
355 P.set_ios_dis_ie (first_iv + iv);
356 }
357 }
358 // ownership_by_dimension[0]: used by connectivity
359 base::_gs.ownership_by_variant [reference_element::p] = tmp_geo_element [reference_element::p].ownership();
360 base::_gs.ownership_by_dimension [0] = tmp_geo_element [reference_element::p].ownership();
361 // ------------------------------------------------------------------------
362 // 6) build connectivity
363 // ------------------------------------------------------------------------
364 build_connectivity (tmp_geo_element);
365 // ------------------------------------------------------------------------
366 // 7) copy tmp_geo_element into _geo_element
367 // ------------------------------------------------------------------------
368 for (size_type dim = 0; dim <= base::_gs._map_dimension; dim++) {
369
371 variant < reference_element:: last_variant_by_dimension(dim); variant++) {
372
373 geo_element::parameter_type param (variant, 1);
374 base::_geo_element [variant].resize (tmp_geo_element [variant].size(), param);
375 size_type igev = 0;
376 iterator_by_variant iter2 = base::_geo_element [variant].begin();
377 for (const_iterator_by_variant_tmp iter = tmp_geo_element [variant].begin(),
378 last = tmp_geo_element [variant].end();
379 iter != last; iter++, iter2++) {
380
381 *iter2 = *iter;
382 }
383 }
384 }
385 // ------------------------------------------------------------------------
386 // 8) set indexes on faces and edges of elements, for P2 approx
387 // ------------------------------------------------------------------------
388 set_element_side_index (1);
389 set_element_side_index (2);
390 // ------------------------------------------------------------------------
391 // 9) node renumbering, when v1 -> v2 aka !have_connectivity
392 // ------------------------------------------------------------------------
393 if (do_upgrade) {
394 // 9.1. compute the new node numbering
395 std::vector<size_type> inod2new_inod (base::_node.size(), std::numeric_limits<size_type>::max());
396 std::vector<size_t> K_new_inod;
397 for (size_type dim = 0; dim < base::_gs._map_dimension; dim++) {
398 for (size_type variant = reference_element::first_variant_by_dimension(base::_gs._map_dimension);
399 variant < reference_element:: last_variant_by_dimension(base::_gs._map_dimension); variant++) {
400 const_iterator_by_variant_tmp tmp_iter = tmp_geo_element [variant].begin();
401 for (const_iterator_by_variant iter = base::_geo_element [variant].begin(),
402 last = base::_geo_element [variant].end();
403 iter != last; iter++, tmp_iter++) {
404
405 const geo_element& K = *iter;
406 const geo_element& tmp_K = *tmp_iter;
407 space_numbering::dis_idof (base::_piola_basis, base::_gs, K, K_new_inod);
408 for (size_type loc_inod = 0; loc_inod < tmp_K.n_node(); loc_inod++) {
409 size_type inod = tmp_K.node (loc_inod);
410 size_type new_inod = K_new_inod [loc_inod];
411 if (inod2new_inod [inod] != std::numeric_limits<size_type>::max()) continue; // already done
412 inod2new_inod [inod] = new_inod;
413 }
414 }
415 }
416 }
417#ifdef _RHEOLEF_PARANO
418 // 9.1.b check permutation
419 std::vector<size_type> new_inod2inod (base::_node.size(), std::numeric_limits<size_type>::max());
420 for (size_type inod = 0, nnod = base::_node.size(); inod < nnod; inod++) {
421 check_macro (inod2new_inod [inod] < nnod, "invalid renumbering: inod2new_inod ["<<inod<<"] = "
422 << inod2new_inod [inod] << " is out of range [0:"<<nnod<<"[");
423 new_inod2inod [inod2new_inod [inod]] = inod;
424 }
425 for (size_type inod = 0, nnod = base::_node.size(); inod < nnod; inod++) {
426 check_macro (new_inod2inod [inod2new_inod [inod]] == inod, "invalid permutation");
427 }
428#endif // _RHEOLEF_PARANO
429 // 9.2. apply the new node numbering to the node table
430 disarray<point_basic<T>,sequential> new_node (base::_node.ownership());
431 for (size_type inod = 0, nnod = base::_node.size(); inod < nnod; inod++) {
432 size_type new_inod = inod2new_inod [inod];
433 new_node [new_inod] = base::_node [inod];
434 }
435 base::_node = new_node;
436 }
437 // ------------------------------------------------------------------------
438 // 11) epilogue
439 // ------------------------------------------------------------------------
440 if (! base::_have_connectivity && ! do_upgrade) {
441 warning_macro ("mesh without connectivity is obsolete (HINT: see geo -upgrade)");
442 }
443 if (do_upgrade) {
444 base::_have_connectivity = true;
445 }
446}
447// ----------------------------------------------------------------------------
448// instanciation in library
449// ----------------------------------------------------------------------------
450template class geo_rep<Float,sequential>;
451
452} // namespace rheolef
field::size_type size_type
Definition branch.cc:430
see the disarray page for the full documentation
Definition disarray.h:497
see the distributor page for the full documentation
Definition distributor.h:69
the finite element boundary domain
base class for M=sequential or distributed meshes representations
Definition geo.h:528
void reset(variant_type variant, size_type order)
see the geo_element page for the full documentation
size_type subgeo_n_node(size_type subgeo_dim, size_type loc_isid) const
size_type n_node() const
size_type subgeo_local_vertex(size_type subgeo_dim, size_type i_subgeo, size_type i_subgeo_vertex) const
size_type first_inod(size_type subgeo_dim) const
size_type size() const
void set_ios_dis_ie(size_type ios_dis_ie)
size_type subgeo_size(size_type subgeo_dim, size_type loc_isid) const
size_type & node(size_type loc_inod)
variant_type variant() const
size_type n_subgeo(size_type subgeo_dim) const
size_type dis_ie() const
size_type subgeo_local_node(size_type subgeo_dim, size_type loc_isid, size_type loc_jsidnod) const
size_type last_inod(size_type subgeo_dim) const
void set_dis_ie(size_type dis_ie)
base::const_iterator_by_variant const_iterator_by_variant
Definition geo.h:799
base::iterator_by_variant iterator_by_variant
Definition geo.h:798
base::size_type size_type
Definition geo.h:791
sequential mesh representation
Definition geo.h:778
idiststream: see the diststream page for the full documentation
Definition diststream.h:336
std::istream & is()
Definition diststream.h:400
void inplace_intersection(const index_set &b)
static const variant_type max_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
variant_type variant() const
#define warning_macro(message)
Definition dis_macros.h:53
check_macro(expr1.have_homogeneous_space(Xh1), "dual(expr1,expr2); expr1 should have homogeneous space. HINT: use dual(interpolate(Xh, expr1),expr2)")
void dis_idof(const basis_basic< T > &b, const geo_size &gs, const geo_element &K, typename std::vector< size_type >::iterator dis_idof_tab)
This file is part of Rheolef.
STL namespace.
point input helper
Definition geo.h:155
size_type map_dimension
Definition geo_header.h:40
size_type dis_size_by_dimension[4]
Definition geo_header.h:44
coordinate_type sys_coord
Definition geo_header.h:41
size_type dimension
Definition geo_header.h:39
size_type dis_size_by_variant[reference_element::max_variant]
Definition geo_header.h:43