Rheolef  7.2
an efficient C++ finite element environment
geo_seq_upgrade.cc
Go to the documentation of this file.
1 #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 
29 namespace 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 
36 template <class T>
37 void
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;
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 }
49 template <class T>
50 void
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
66  size_type size_by_variant [reference_element::max_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
70  size_type nsid = 0;
71  {
72  disarray<index_set,sequential> ball_S (geo_element_ownership(0), empty_set);
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;
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 // =========================================================================
186 template <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  // ------------------------------------------------------------------------
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 }
234 template <class T>
235 void
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);
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  // ------------------------------------------------------------------------
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  {
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 
372 
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 // ----------------------------------------------------------------------------
450 template class geo_rep<Float,sequential>;
451 
452 } // 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
idiststream & get(idiststream &ips, const geo_rep< T, sequential > &omega, std::vector< index_set > *ball)
geo_element_hack::size_type size_type
Definition: geo.h:260
hack_array< geo_element_hack >::iterator iterator_by_variant
Definition: geo.h:269
hack_array< geo_element_hack >::const_iterator const_iterator_by_variant
Definition: geo.h:278
base class for M=sequential or distributed meshes representations
Definition: geo.h:528
void reset(variant_type variant, size_type order)
Definition: geo_element.h:426
see the geo_element page for the full documentation
Definition: geo_element.h:102
size_type subgeo_n_node(size_type subgeo_dim, size_type loc_isid) const
Definition: geo_element.h:217
size_type n_node() const
Definition: geo_element.h:170
size_type subgeo_local_vertex(size_type subgeo_dim, size_type i_subgeo, size_type i_subgeo_vertex) const
Definition: geo_element.h:223
size_type first_inod(size_type subgeo_dim) const
Definition: geo_element.h:225
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 subgeo_size(size_type subgeo_dim, size_type loc_isid) const
Definition: geo_element.h:221
variant_type variant() const
Definition: geo_element.h:161
size_type & node(size_type loc_inod)
Definition: geo_element.h:185
size_type n_subgeo(size_type subgeo_dim) const
Definition: geo_element.h:212
size_type dis_ie() const
Definition: geo_element.h:163
size_type subgeo_local_node(size_type subgeo_dim, size_type loc_isid, size_type loc_jsidnod) const
Definition: geo_element.h:219
size_type last_inod(size_type subgeo_dim) const
Definition: geo_element.h:227
void set_dis_ie(size_type dis_ie)
Definition: geo_element.h:172
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
bool good() const
Definition: diststream.cc:124
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
static std::set< size_t > empty_set
Definition: space_seq.cc:59
size_t size_type
Definition: basis_get.cc:76
#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)
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
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