Rheolef  7.2
an efficient C++ finite element environment
msh2geo.cc
Go to the documentation of this file.
1 //
4 // Copyright (C) 2000-2009 Pierre Saramito <Pierre.Saramito@imag.fr>
5 //
6 // Rheolef is free software; you can redistribute it and/or modify
7 // it under the terms of the GNU General Public License as published by
8 // the Free Software Foundation; either version 2 of the License, or
9 // (at your option) any later version.
10 //
11 // Rheolef is distributed in the hope that it will be useful,
12 // but WITHOUT ANY WARRANTY; without even the implied warranty of
13 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 // GNU General Public License for more details.
15 //
16 // You should have received a copy of the GNU General Public License
17 // along with Rheolef; if not, write to the Free Software
18 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
19 //
20 // =========================================================================
21 // the msh2geo unix command
22 // author: Pierre.Saramito@imag.fr
23 
24 namespace rheolef {
122 } // namespace rheolef
123 
124 #include "rheolef/compiler.h" // after index_set to avoid boost
125 #include "scatch.icc"
126 #include "index_set_header.icc"
127 #include "index_set_body.icc"
128 #include <array>
129 #include <cstring>
130 
131 namespace rheolef {
132 
133 template<class T>
134 struct point_basic: std::array<T,3> {
135  typedef std::array<T,3> base;
136  point_basic(T x0=T(), T x1=T(), T x2=T())
137  : base() {
138  base::operator[](0) = x0;
139  base::operator[](1) = x1;
140  base::operator[](2) = x2;
141  }
142 };
144 
145 namespace edge {
146 #include "edge.icc"
147 } // namespace edge
148 
149 namespace triangle {
150 #include "triangle.icc"
151 } // namespace triangle
152 
153 namespace quadrangle {
154 #include "quadrangle.icc"
155 } // namespace quadrangle
156 
157 namespace tetrahedron {
158 #include "tetrahedron.icc"
159 } // namespace tetrahedron
160 
161 namespace prism {
162 #include "prism.icc"
163 } // namespace prism
164 
165 namespace hexahedron {
166 #include "hexahedron.icc"
167 } // namespace hexahedron
168 
169 } // namespace rheolef
170 
171 #include "reference_element_aux.icc"
172 
173 namespace rheolef {
174 
175 // TODO: move somewhere in reference_element_xxx.cc
176 const char
178  'p',
179  'e',
180  't',
181  'q',
182  'T',
183  'P',
184  'H'
185 };
186 } // namespace rheolef
187 
188 namespace rheolef {
189 
190 typedef int orientation_type;
191 typedef int shift_type;
192 
193 // side (edge & face): element with reduced node list when order > 1, for sides (edges or faces)
194 struct geo_element_side: std::array<size_t,4> {
195  void setname (char name) { _variant = reference_element_variant(name); }
196  void setindex (size_t index) { _index = index; }
197  char name() const { return _reference_element_name [_variant]; }
198  size_t index() const { return _index; }
199  size_t variant() const { return _variant; }
200  size_t dimension()const { return reference_element_dimension_by_variant(variant()); }
201  size_t n_vertex() const { return reference_element_n_vertex (variant()); }
202  size_t size() const { return n_vertex(); }
203  geo_element_side()
204  : array<size_t,4>(), _variant(-1), _index(-1)
205  {}
206 protected:
207  size_t _variant;
208  size_t _index;
209 };
210 // side index with orient and shift for geo_element
211 struct geo_element_indirect {
212  typedef int orientation_type;
213  typedef int shift_type;
214  void setindex (size_t index) { _index = index; }
215  void setorientation (size_t orient) { _orient = orient; }
216  void setshift (size_t shift) { _shift = shift; }
217  size_t index() const { return _index; }
218  int orientation() const { return _orient; }
219  size_t shift() const { return _shift; }
221  : _index(-1), _orient(1), _shift(0) {}
222 protected:
223  size_t _index;
224  int _orient;
225  int _shift;
226 };
227 // element with full node list when order > 1
228 struct geo_element: vector<size_t> {
229  void setname (char name) { _name = name; }
230  void setorder (size_t order) { _order = order; }
231  void setindex (size_t index) { _index = index; }
232  void setgmshtype (size_t gmshtype) { _gmshtype = gmshtype; }
233  geo_element_indirect& edge_indirect (size_t i) { return _edge[i]; }
234  geo_element_indirect& face_indirect (size_t i) { return _face[i]; }
235  char name() const { return _name; }
236  size_t order() const { return _order; }
237  size_t index() const { return _index; }
238  size_t gmshtype() const { return _gmshtype; }
239  size_t variant() const { return reference_element_variant (name()); }
240  size_t dimension()const { return reference_element_dimension_by_name(name()); }
241  size_t n_vertex() const { return reference_element_n_vertex (variant()); }
242  const geo_element_indirect& edge_indirect (size_t i) const { return _edge[i]; }
243  const geo_element_indirect& face_indirect (size_t i) const { return _face[i]; }
244  size_t n_edge() const { return _reference_element_n_edge_by_variant [variant()]; }
245  size_t edge (size_t i) const { return _edge[i].index(); }
246  size_t edge_local_vertex(size_t iedg, size_t edg_iv) const;
247  size_t face (size_t i) const { return _face[i].index(); }
248  size_t n_face() const { return _reference_element_n_face_by_variant [variant()]; }
249  size_t face_n_vertex(size_t loc_ifac) const;
250  size_t face_local_vertex(size_t iedg, size_t edg_iv) const;
252  : vector<size_t>(), _name('z'), _order(-1), _index(-1), _gmshtype(-1), _edge(), _face()
253  {}
254 protected:
255  char _name; // TODO: store variant instead of name
256  size_t _order;
257  size_t _index;
258  size_t _gmshtype;
259  array<geo_element_indirect,12> _edge;
260  array<geo_element_indirect,6> _face;
261 };
262 size_t
263 geo_element::edge_local_vertex (size_t iedg, size_t edg_iv) const
264 {
265  switch (variant()) {
266  case reference_element__e: return vector<size_t>::operator[] (edg_iv);
267  case reference_element__t: return triangle::edge [iedg][edg_iv%2];
268  case reference_element__q: return quadrangle::edge [iedg][edg_iv%2];
269  case reference_element__T: return tetrahedron::edge [iedg][edg_iv%2];
270  case reference_element__P: return prism::edge [iedg][edg_iv%2];
271  case reference_element__H: return hexahedron::edge [iedg][edg_iv%2];
272  default: error_macro ("invalid variant"); return 0;
273  }
274 }
275 size_t
276 geo_element::face_local_vertex (size_t loc_ifac, size_t fac_iv) const
277 {
278  switch (variant()) {
280  case reference_element__q: return vector<size_t>::operator[] (fac_iv);
281  case reference_element__T: return tetrahedron::face [loc_ifac][fac_iv%3];
282  case reference_element__P: return prism::face [loc_ifac][fac_iv%4];
283  case reference_element__H: return hexahedron::face [loc_ifac][fac_iv%4];
284  default: error_macro ("invalid variant"); return 0;
285  }
286 }
287 size_t
288 geo_element::face_n_vertex(size_t loc_ifac) const
289 {
290  switch (variant()) {
291  case reference_element__t: return 3;
292  case reference_element__q: return 4;
293  case reference_element__T: return 3;
294  case reference_element__H: return 4;
295  case reference_element__P: return (loc_ifac < 2) ? 3 : 4;
296  default: error_macro ("invalid variant"); return 0;
297  }
298 }
299 
300 } // namespace rheolef
301 
302 #include "msh2geo_defs.icc"
303 #include "msh2geo_node_renum.icc"
304 #include "geo_element_aux.icc"
305 
306 using namespace std;
307 using namespace rheolef;
308 // necessite de numeroter les aretes
309 // puis de le stoker dans ts les triangles
310 // le num arete et son orient
311 // => doit etre appeler apres la num des aretes
312 void
314  size_t order,
315  const array<size_t,reference_element__max_variant>& size_by_variant,
316  const array<size_t,reference_element__max_variant+1>& first_by_variant,
317  const array<size_t,reference_element__max_variant>& loc_ndof_by_variant,
318  const geo_element& K,
319  vector<size_t>& dis_idof1)
320 {
321  dis_idof1.resize (reference_element_n_node (K.variant(), order));
322  std::fill (dis_idof1.begin(), dis_idof1.end(), std::numeric_limits<size_t>::max());
323 
324  for (size_t subgeo_variant = 0, variant = K.variant(); subgeo_variant <= variant; subgeo_variant++) {
325  size_t subgeo_dim = reference_element_dimension_by_variant (subgeo_variant);
326  size_t loc_sub_ndof = loc_ndof_by_variant [subgeo_variant];
327  for (size_t first_loc_idof = reference_element_first_inod_by_variant (variant, order, subgeo_variant),
328  last_loc_idof = reference_element_last_inod_by_variant (variant, order, subgeo_variant),
329  loc_idof = first_loc_idof;
330  loc_idof < last_loc_idof; loc_idof++) {
331  // 1) local loc_igev on subgeo
332  size_t loc_igev = (loc_idof - first_loc_idof) / loc_sub_ndof;
333  size_t loc_igev_j = (loc_idof - first_loc_idof) % loc_sub_ndof;
334  // 2) then compute ige; computation depends upon the subgeo dimension:
335  size_t ige;
336  switch (subgeo_dim) {
337  case 0: {
338  // convert node numbering to vertex numbering, for geo order > 1
339  size_t loc_inod = loc_idof; // only one dof per vertex
340  ige = K [loc_inod];
341  break;
342  }
343  case 1: {
344  loc_igev_j = geo_element_fix_edge_indirect (K, loc_igev, loc_igev_j, order);
345  size_t loc_ige = loc_igev;
346  ige = K.edge_indirect(loc_ige).index();
347  break;
348  }
349  case 2: {
350  size_t loc_ige = loc_igev;
351  if (subgeo_variant == reference_element__t) {
352  loc_igev_j = geo_element_fix_triangle_indirect (K, loc_igev, loc_igev_j, order);
353  } else {
354  size_t loc_ntri = (K.variant() == reference_element__P) ? 2 : 0;
355  loc_ige += loc_ntri;
356  loc_igev_j = geo_element_fix_quadrangle_indirect (K, loc_ige, loc_igev_j, order);
357  }
358  ige = K.face(loc_ige);
359  break;
360  }
361  case 3: {
362  ige = K.index();
363  break;
364  }
365  default: {
366  error_macro ("unexpected subgeo_dim="<<subgeo_dim);
367  }
368  }
369  // 3) then compute igev, by variant
370  size_t igev = ige;
371  for (size_t prev_subgeo_variant = reference_element_first_variant_by_dimension(subgeo_dim);
372  prev_subgeo_variant < subgeo_variant;
373  prev_subgeo_variant++) {
374  size_t nprev = size_by_variant [prev_subgeo_variant];
375  assert_macro (ige >= nprev, "invalid index");
376  igev -= nprev;
377  }
378  size_t dis_igev = igev;
379  // 4) finally compute dis_idof
380  size_t dis_idof = 0;
381  for (size_t prev_subgeo_variant = 0;
382  prev_subgeo_variant < subgeo_variant;
383  prev_subgeo_variant++) {
384  dis_idof += size_by_variant [prev_subgeo_variant]
385  *loc_ndof_by_variant [prev_subgeo_variant];
386  }
387  dis_idof += dis_igev
388  *loc_ndof_by_variant [subgeo_variant]
389  + loc_igev_j;
390  assert_macro (loc_idof < dis_idof1.size(), "msh2geo: invalid index");
391  dis_idof1 [loc_idof] = dis_idof;
392  }
393  }
394 }
395 void
397  ostream& out,
398  size_t dim,
399  size_t dom_dim,
400  const map<size_t,list<size_t> >& domain_map,
401  map<size_t, string>& phys,
402  const vector<geo_element>& element)
403 {
404  if (dim == dom_dim && domain_map.size() == 1) return;
405  for (map<size_t,list<size_t> >::const_iterator
406  first = domain_map.begin(),
407  last = domain_map.end(); first != last; first++) {
408  size_t dom_idx = (*first).first;
409  const list<size_t>& dom = (*first).second;
410  string dom_name = phys[dom_idx];
411  if (dom_name == "") dom_name = "unnamed" + std::to_string(dom_idx);
412  out << "domain" << endl
413  << dom_name << endl
414  << "1 " << dom_dim << " " << dom.size() << endl;
415  for (size_t variant = reference_element_first_variant_by_dimension(dom_dim);
416  variant < reference_element_last_variant_by_dimension(dom_dim); variant++) {
417  for (list<size_t>::const_iterator fe = dom.begin(), le = dom.end(); fe != le; fe++) {
418  size_t ie = *fe;
419  if (element[ie].variant() != variant) continue;
420  out << element[ie].name() << "\t";
421  if (element[ie].order() > 1) {
422  out << "p" << element[ie].order() << " ";
423  }
424  for (size_t j = 0; j < element[ie].size(); j++) {
425  out << element[ie][j] << " ";
426  }
427  out << endl;
428  }
429  }
430  out << endl;
431  }
432 }
433 void
435  ostream& out,
436  size_t dim,
437  size_t dom_dim,
438  const map<size_t,list<size_t> >& domain_map,
439  map<size_t, string>& phys,
440  const vector<geo_element>& element,
441  const vector<geo_element_side>& edge,
442  const vector<index_set>& edge_ball,
443  const vector<geo_element_side>& face,
444  const vector<index_set>& face_ball)
445 {
446  if (dim == dom_dim && domain_map.size() == 1) return;
447  for (map<size_t,list<size_t> >::const_iterator
448  first = domain_map.begin(),
449  last = domain_map.end(); first != last; first++) {
450  size_t dom_idx = (*first).first;
451  const list<size_t>& dom = (*first).second;
452  string dom_name = phys[dom_idx];
453  if (dom_name == "") dom_name = "unnamed" + std::to_string(dom_idx);
454  out << "domain" << endl
455  << dom_name << endl
456  << "2 " << dom_dim << " " << dom.size() << endl;
457  for (size_t variant = reference_element_first_variant_by_dimension(dom_dim);
458  variant < reference_element_last_variant_by_dimension(dom_dim); variant++) {
459  for (list<size_t>::const_iterator fe = dom.begin(), le = dom.end(); fe != le; fe++) {
460  size_t ie = *fe;
461  if (element[ie].variant() != variant) continue;
462  size_t isid = 0;
463  orientation_type orient = 1;
464  if (element[ie].name() == 'p') {
465  isid = element[ie][0];
466  orient = 1;
467  } else if (element[ie].name() == 'e') {
468  size_t inod0 = element[ie][0];
469  size_t inod1 = element[ie][1];
470  index_set iedge_set = edge_ball[inod0];
471  iedge_set.inplace_intersection (edge_ball[inod1]);
472  check_macro (iedge_set.size() == 1,
473  "msh2geo: error: connectivity problem (iedge.size="<<iedge_set.size()<<")");
474  isid = *(iedge_set.begin());
475  orient = (inod0 == edge[isid][0]) ? 1 : -1;
476  } else if (element[ie].dimension() == 2) {
477  size_t inod0 = element[ie][0];
478  index_set iface_set = face_ball[inod0];
479  for (size_t j = 1; j < element[ie].n_vertex(); ++j) {
480  size_t inodj = element[ie][j];
481  iface_set.inplace_intersection (face_ball[inodj]);
482  }
483  check_macro (iface_set.size() == 1,
484  "msh2geo: error: connectivity problem (iface.size="<<iface_set.size()<<")");
485  isid = *(iface_set.begin());
486  shift_type shift;
487  if (element[ie].name() == 't') {
488  geo_element_get_orientation_and_shift (
489  face[isid],
490  element[ie][0],
491  element[ie][1],
492  element[ie][2],
493  orient, shift);
494  } else {
495  geo_element_get_orientation_and_shift (
496  face[isid],
497  element[ie][0],
498  element[ie][1],
499  element[ie][2],
500  element[ie][3],
501  orient, shift);
502  }
503  } else { // 3d domain
504  isid = element[ie].index();
505  orient = 1;
506  }
507  out << orient*int(isid) << endl;
508  }
509  }
510  out << endl;
511  }
512 }
513 void msh2geo (istream& in, ostream& out, string sys_coord_name, bool do_upgrade)
514 {
515  // ----------------------------------
516  // 1. input gmsh
517  // ----------------------------------
518  // 1.0. preambule
519  check_macro (scatch(in,"$MeshFormat",true),
520  "msh2geo: input stream does not contains a gmsh mesh file ($MeshFormat not found).");
521  double gmsh_fmt_version;
522  size_t file_type, float_data_size;
523  in >> gmsh_fmt_version >> file_type >> float_data_size;
524  check_macro (gmsh_fmt_version <= 2.2,
525  "gmsh format version " << gmsh_fmt_version << " founded ; expect version 2.2 (HINT: use gmsh -format msh2)");
526  check_macro (file_type == 0, "msh2geo: unsupported gmsh non-ascii format");
527  check_macro (scatch(in,"$EndMeshFormat",true), "msh2geo: gmsh input error: $EndMeshFormat not found.");
528  //
529  // 1.1 optional domain names
530  //
531  bool full_match = true;
532  bool partial_match = !full_match;
533  check_macro (scatch(in,"$",partial_match), "msh2geo: gmsh input error: no more label found.");
534  size_t nphys;
535  map<size_t, string> phys;
536  string label;
537  in >> label;
538  size_t n_names = 0;
539  if (label == "PhysicalNames") {
540  in >> nphys;
541  for (size_t i = 0; i < nphys; i++) {
542  string name;
543  size_t dom_dim, dom_idx;
544  in >> dom_dim >> dom_idx;
545  // get name:
546  char c;
547  in >> std::ws >> c;
548  if (c != '"') name.push_back(c);
549  do {
550  in.get(c);
551  name.push_back(c);
552  } while (c != '"' && c != '\n');
553  // strip '"' in name
554  size_t start = 0, end = name.length();
555  if (name[start] == '"') start++;
556  if (name[end-1] == '"') end--;
557  name = name.substr(start,end-start);
558  // rename spaces and tabs
559  for (size_t i = 0; i < name.size(); i++) {
560  if (name[i] == ' ' || name[i] == '\t') name[i] = '_';
561  }
562  phys[dom_idx] = name.substr(start,end-start);
563  }
564  check_macro (scatch(in,"$EndPhysicalNames",true), "msh2geo: gmsh input error ($EndPhysicalNames not found).");
565  check_macro (scatch(in,"$",partial_match), "msh2geo: gmsh input error: no more label found.");
566  in >> label;
567  }
568  //
569  // 1.2. nodes
570  //
571  if (label != "Nodes") {
572  check_macro (scatch(in,"$Nodes",true), "msh2geo: $Nodes not found in gmsh file");
573  }
574  size_t nnode;
575  in >> nnode;
576  vector<point> node (nnode);
577  double infty = numeric_limits<double>::max();
578  point xmin ( infty, infty, infty);
579  point xmax (-infty, -infty, -infty);
580  for (size_t k = 0; k < node.size(); k++) {
581  size_t dummy;
582  in >> dummy;
583  for (size_t i = 0; i < 3; ++i) {
584  in >> node[k][i];
585  }
586  for (size_t j = 0 ; j < 3; j++) {
587  xmin[j] = min(node[k][j], xmin[j]);
588  xmax[j] = max(node[k][j], xmax[j]);
589  }
590  }
591  //
592  // dimension is deduced from bounding box
593  //
594  size_t dim = 3;
595  if (xmax[2] == xmin[2]) {
596  dim = 2;
597  if (xmax[1] == xmin[1]) {
598  dim = 1;
599  if (xmax[0] == xmin[0]) dim = 0;
600  }
601  }
602  //
603  // 1.3. elements
604  //
605  check_macro (scatch(in,"$Elements",true), "msh2geo: $Elements not found in gmsh file");
606  size_t n_element;
607  in >> n_element;
608  vector<geo_element> element (n_element);
609  vector<size_t> node_subgeo_variant (node.size(), std::numeric_limits<size_t>::max());
610  array<map<size_t, list<size_t> >,4> domain_map; // domain_map[dim][idom][ie]
611  size_t map_dim = 0;
612  size_t order = 0;
613  for (size_t i = 0; i < element.size(); i++) {
614  size_t id, dummy, gmshtype;
615  in >> id >> gmshtype;
616  size_t n_tag_gmsh;
617  in >> n_tag_gmsh;
618  size_t domain_idx = 0;
619  for (size_t j = 0 ; j < n_tag_gmsh; j++) {
620  // the first tag is the physical domain index
621  // the second tag is the object index, defined for all elements
622  // the third is zero (in all examples)
623  size_t tag_dummy;
624  in >> tag_dummy;
625  if (j == 0) {
626  domain_idx = tag_dummy;
627  }
628  }
629  check_macro (gmshtype < gmshtype_max,
630  "msh2geo: element #" << id << ": unexpected gmsh type '" << gmshtype << "'");
631  check_macro (gmsh_table[gmshtype].supported,
632  "msh2geo: element #" << id << ": unsupported gmsh type '" << gmshtype << "'");
633 
634  element[i].setname (gmsh_table[gmshtype].name);
635  element[i].setorder(gmsh_table[gmshtype].order);
636  element[i].setgmshtype(gmshtype);
637  size_t nv = gmsh_table[gmshtype].nv;
638  size_t nn = gmsh_table[gmshtype].nn_tot;
639  size_t dim_i = element[i].dimension();
640  size_t variant = element[i].variant();
641  map_dim = max(map_dim,dim_i);
642  if (order == 0) {
643  order = element[i].order();
644  } else {
645  check_macro (order == element[i].order(), "msh2geo: unexpected order="<<element[i].order());
646  }
647  element[i].resize(nn);
648  for (size_t j = 0; j < nn; j++) {
649  in >> element[i][j];
650  element[i][j]--;
651  }
652  for (size_t subgeo_variant = 0; subgeo_variant <= variant; subgeo_variant++) { // set node dimension
653  for (size_t loc_inod = reference_element_first_inod_by_variant(variant,element[i].order(),subgeo_variant),
654  loc_nnod = reference_element_last_inod_by_variant(variant,element[i].order(),subgeo_variant);
655  loc_inod < loc_nnod; loc_inod++) {
656  node_subgeo_variant [element[i][loc_inod]] = subgeo_variant;
657  }
658  }
659  // from gmsh to rheolef/geo local node renumbering: element[i][j] are modified
660  msh2geo_node_renum (element[i], element[i].name(), element[i].order());
661  if (domain_idx != 0) {
662  domain_map[dim_i][domain_idx].push_back(i);
663  }
664  }
665  array<size_t,reference_element__max_variant> loc_ndof_by_variant;
666  reference_element_init_local_nnode_by_variant (order, loc_ndof_by_variant);
667  // set dis_ie: by increasing variant order, for map_dim element only
668  size_t last_index = 0;
669  for (size_t variant = reference_element_first_variant_by_dimension(map_dim);
670  variant < reference_element_last_variant_by_dimension(map_dim); variant++) {
671  for (size_t ie = 0; ie < element.size(); ++ie) {
672  if (element[ie].variant() != variant) continue;
673  element[ie].setindex (last_index);
674  switch (element[ie].dimension()) {
675  case 0: element[ie].resize(1); element[ie][0] = last_index; break;
676  case 1: element[ie].edge_indirect(0).setindex(last_index); break;
677  case 2: element[ie].face_indirect(0).setindex(last_index); break;
678  default: break;
679  }
680  last_index++;
681  }
682  }
683  // -------------------------------------------------
684  // 2. node reordering, first pass
685  // -------------------------------------------------
686  // permut: first vertex, then edge, faces and inernal nodes, by increasing subgeo_dim
687  vector<size_t> old2new_inode (node.size(), std::numeric_limits<size_t>::max());
688  if (true || !do_upgrade) {
689  // 2.1. when no upgrade
690  size_t new_inode = 0;
691  for (size_t subgeo_variant = 0; subgeo_variant < reference_element_last_variant_by_dimension(map_dim); subgeo_variant++) {
692  for (size_t old_inode = 0; old_inode < node.size(); old_inode++) {
693  if (node_subgeo_variant [old_inode] != subgeo_variant) continue;
694  old2new_inode[old_inode] = new_inode++;
695  }
696  }
697  } else {
698  // 2.2. when no upgrade: will be renumbered after side creation
699  for (size_t inode = 0; inode < node.size(); inode++) {
700  old2new_inode[inode] = inode;
701  }
702  }
703  for (size_t inode = 0; inode < node.size(); inode++) {
704  check_macro (old2new_inode[inode] != std::numeric_limits<size_t>::max(), "msh2geo: invalid permutation");
705  }
706  // inv permut
707  vector<size_t> new2old_inode (node.size(), std::numeric_limits<size_t>::max());
708  for (size_t inode = 0; inode < node.size(); inode++) {
709  new2old_inode[old2new_inode[inode]] = inode;
710  }
711  for (size_t inode = 0; inode < node.size(); inode++) {
712  check_macro (new2old_inode[inode] != std::numeric_limits<size_t>::max(), "msh2geo: invalid inverse permutation for inode="<<inode);
713  }
714  // for convenience, apply the new numbering to nodes and elements
715  {
716  vector<point> new_node (node.size());
717  for (size_t old_inode = 0; old_inode < node.size(); old_inode++) {
718  new_node [old2new_inode[old_inode]] = node [old_inode];
719  }
720  for (size_t inode = 0; inode < node.size(); inode++) {
721  node [inode] = new_node [inode];
722  }
723  }
724  for (size_t ie = 0; ie < element.size(); ++ie) {
725  for (size_t i = 0; i < element[ie].size(); ++i) {
726  element[ie][i] = old2new_inode[element[ie][i]];
727  }
728  }
729  // ------------------------------------------------------------------------
730  // 3.1) compute all edges
731  // ------------------------------------------------------------------------
732  vector<index_set> edge_ball (node.size());
733  vector<geo_element_side> edge;
734  if (do_upgrade && map_dim >= 2) {
735  list <geo_element_side> edge_list;
736  size_t nedg = 0;
737  for (size_t ie = 0; ie < element.size(); ++ie) {
738  if (element[ie].dimension() != map_dim) continue;
739  for (size_t iedg = 0; iedg < element[ie].n_edge(); ++iedg) {
740  size_t iv0 = element[ie].edge_local_vertex(iedg, 0);
741  size_t iv1 = element[ie].edge_local_vertex(iedg, 1);
742  size_t inod0 = element[ie][iv0];
743  size_t inod1 = element[ie][iv1];
744  index_set iedge_set = edge_ball[inod0];
745  iedge_set.inplace_intersection (edge_ball[inod1]);
746  check_macro (iedge_set.size() <= 1,
747  "msh2geo: error: connectivity problem (iedge.size="<<iedge_set.size()<<")");
748  if (iedge_set.size() == 1) continue; // edge already exists
749  edge_ball[inod0].insert (nedg);
750  edge_ball[inod1].insert (nedg);
751  geo_element_side new_edge;
752  new_edge.setname('e');
753  new_edge[0] = inod0;
754  new_edge[1] = inod1;
755  new_edge.setindex(nedg);
756  edge_list.push_back (new_edge);
757  nedg++;
758  }
759  }
760  edge.resize (edge_list.size());
761  std::copy (edge_list.begin(), edge_list.end(), edge.begin());
762  }
763  // propagate edge numbering inside elements
764  if (do_upgrade && map_dim >= 2) {
765  for (size_t ie = 0; ie < element.size(); ++ie) {
766  if (element[ie].dimension() != map_dim) continue;
767  for (size_t iedg = 0; iedg < element[ie].n_edge(); ++iedg) {
768  size_t iv0 = element[ie].edge_local_vertex(iedg, 0);
769  size_t iv1 = element[ie].edge_local_vertex(iedg, 1);
770  size_t inod0 = element[ie][iv0];
771  size_t inod1 = element[ie][iv1];
772  index_set iedge_set = edge_ball[inod0];
773  iedge_set.inplace_intersection (edge_ball[inod1]);
774  check_macro (iedge_set.size() == 1,
775  "msh2geo: error: connectivity problem (iedge.size="<<iedge_set.size()<<")");
776  size_t iedge = *(iedge_set.begin());
777  element[ie].edge_indirect(iedg).setindex (iedge);
778  element[ie].edge_indirect(iedg).setorientation((edge[iedge][0] == inod0) ? 1 : -1);
779  }
780  }
781  }
782  // ------------------------------------------------------------------------
783  // 3.2) compute all faces
784  // ------------------------------------------------------------------------
785  vector<index_set> face_ball (node.size());
786  vector<geo_element_side> face;
787  if (do_upgrade && map_dim >= 3) {
788  list <geo_element_side> face_list;
789  size_t nfac = 0;
790  for (size_t ie = 0; ie < element.size(); ++ie) {
791  if (element[ie].dimension() != map_dim) continue;
792  for (size_t loc_ifac = 0; loc_ifac < element[ie].n_face(); ++loc_ifac) {
793  size_t iv0 = element[ie].face_local_vertex(loc_ifac, 0);
794  size_t inod0 = element[ie][iv0];
795  index_set iface_set = face_ball[inod0];
796  for (size_t j = 1; j < element[ie].face_n_vertex(loc_ifac); ++j) {
797  size_t ivj = element[ie].face_local_vertex(loc_ifac, j);
798  size_t inodj = element[ie][ivj];
799  iface_set.inplace_intersection (face_ball[inodj]);
800  }
801  check_macro (iface_set.size() <= 1,
802  "msh2geo: error: connectivity problem (iface.size="<<iface_set.size()<<")");
803  if (iface_set.size() == 1) continue; // face already exists
804  geo_element_side new_face;
805  char face_name = (element[ie].face_n_vertex(loc_ifac) == 3) ? 't' : 'q';
806  new_face.setname(face_name);
807  new_face.setindex(nfac);
808  for (size_t j = 0; j < element[ie].face_n_vertex(loc_ifac); ++j) {
809  size_t ivj = element[ie].face_local_vertex(loc_ifac, j);
810  size_t inodj = element[ie][ivj];
811  face_ball[inodj].insert (nfac);
812  new_face[j] = inodj;
813  }
814  face_list.push_back (new_face);
815  nfac++;
816  }
817  }
818  face.resize (face_list.size());
819  std::copy (face_list.begin(), face_list.end(), face.begin());
820  }
821  // propagate face numbering inside elements
822  if (do_upgrade && map_dim >= 3) {
823  for (size_t ie = 0; ie < element.size(); ++ie) {
824  if (element[ie].dimension() != map_dim) continue;
825  for (size_t loc_ifac = 0; loc_ifac < element[ie].n_face(); ++loc_ifac) {
826  size_t iv0 = element[ie].face_local_vertex(loc_ifac, 0);
827  size_t inod0 = element[ie][iv0];
828  index_set iface_set = face_ball[inod0];
829  for (size_t j = 1; j < element[ie].face_n_vertex(loc_ifac); ++j) {
830  size_t ivj = element[ie].face_local_vertex(loc_ifac, j);
831  size_t inodj = element[ie][ivj];
832  iface_set.inplace_intersection (face_ball[inodj]);
833  }
834  check_macro (iface_set.size() == 1,
835  "msh2geo: error: connectivity problem (iface.size="<<iface_set.size()<<")");
836  size_t iface = *(iface_set.begin());
837  orientation_type orient;
838  shift_type shift;
839  if (element[ie].face_n_vertex(loc_ifac) == 3) {
840  geo_element_get_orientation_and_shift (
841  face[iface],
842  element[ie][element[ie].face_local_vertex(loc_ifac, 0)],
843  element[ie][element[ie].face_local_vertex(loc_ifac, 1)],
844  element[ie][element[ie].face_local_vertex(loc_ifac, 2)],
845  orient, shift);
846  } else {
847  geo_element_get_orientation_and_shift (
848  face[iface],
849  element[ie][element[ie].face_local_vertex(loc_ifac, 0)],
850  element[ie][element[ie].face_local_vertex(loc_ifac, 1)],
851  element[ie][element[ie].face_local_vertex(loc_ifac, 2)],
852  element[ie][element[ie].face_local_vertex(loc_ifac, 3)],
853  orient, shift);
854  }
855  element[ie].face_indirect(loc_ifac).setindex (iface);
856  element[ie].face_indirect(loc_ifac).setorientation (orient);
857  element[ie].face_indirect(loc_ifac).setshift (shift);
858  }
859  }
860  }
861  // ------------------------------------------------------------------------
862  // 3.4) count all elements
863  // ------------------------------------------------------------------------
864  array<size_t,reference_element__max_variant> size_by_variant;
865  array<size_t,reference_element__max_variant+1> first_by_variant;
866  size_by_variant.fill (0);
867  first_by_variant.fill (0);
868  if (true) {
869  size_t n_vertex = 0;
870  for (size_t inode = 0; inode < node.size(); inode++) {
871  if (node_subgeo_variant [inode] == reference_element__p) n_vertex++;
872  }
873  size_by_variant [reference_element__p] = n_vertex;
874  if (map_dim >= 2 && edge.size() != 0) {
875  size_by_variant [reference_element__e] = edge.size();
876  }
877  if (map_dim >= 3 && face.size() != 0) {
878  for (size_t loc_ifac = 0; loc_ifac < face.size(); ++loc_ifac) {
879  size_by_variant [face[loc_ifac].variant()]++;
880  }
881  }
882  for (size_t ie = 0; ie < element.size(); ++ie) {
883  if (element[ie].dimension() != map_dim) continue;
884  size_t variant = element[ie].variant();
885  size_by_variant [variant]++;
886  }
887  for (size_t dim = 0; dim <= 3; ++dim) {
888  for (size_t variant = reference_element_first_variant_by_dimension(dim);
889  variant < reference_element_last_variant_by_dimension(dim); variant++) {
890  first_by_variant [variant+1] = size_by_variant [variant];
891  }
892  }
893  }
894  // -------------------------------------------------
895  // 4. node reordering, second pass
896  // -------------------------------------------------
897  // need edges & faces with index & orient
898  if (do_upgrade) {
899  std::fill (old2new_inode.begin(), old2new_inode.end(), std::numeric_limits<size_t>::max());
900  std::fill (new2old_inode.begin(), new2old_inode.end(), std::numeric_limits<size_t>::max());
901  for (size_t variant = reference_element_first_variant_by_dimension(map_dim);
902  variant < reference_element_last_variant_by_dimension(map_dim); variant++) {
903  for (size_t ie = 0; ie < element.size(); ++ie) {
904  if (element[ie].variant() != variant) continue;
905  const geo_element& old_K = element[ie];
906  std::vector<size_t> new_K;
908  order, size_by_variant, first_by_variant, loc_ndof_by_variant, old_K, new_K);
909  for (size_t loc_inod = 0; loc_inod < old_K.size(); loc_inod++) {
910  size_t old_inod = old_K [loc_inod];
911  size_t new_inod = new_K [loc_inod];
912  if (old2new_inode [old_inod] != std::numeric_limits<size_t>::max()) continue; // already done
913  old2new_inode [old_inod] = new_inod;
914  }
915  }
916  }
917  } else {
918  // when no upgrade:
919  for (size_t inode = 0; inode < node.size(); inode++) {
920  old2new_inode[inode] = inode;
921  }
922  }
923  for (size_t inode = 0; inode < node.size(); inode++) {
924  check_macro (old2new_inode[inode] != std::numeric_limits<size_t>::max(), "msh2geo: invalid permutation");
925  }
926  // inv permut
927  for (size_t inode = 0; inode < node.size(); inode++) {
928  new2old_inode[old2new_inode[inode]] = inode;
929  }
930  for (size_t inode = 0; inode < node.size(); inode++) {
931  check_macro (new2old_inode[inode] != std::numeric_limits<size_t>::max(), "msh2geo: invalid inverse permutation for inode="<<inode);
932  }
933 #ifdef TODO
934 #endif // TODO
935  // ----------------------------------
936  // 5. output geo
937  // ----------------------------------
938  // 5.1. output the mesh
939  size_t version = 4;
940  static const char* geo_variant_name [reference_element__max_variant] = {
941  "points",
942  "edges",
943  "triangles",
944  "quadrangles",
945  "tetrahedra",
946  "prisms",
947  "hexahedra"
948  };
949  out << setprecision(numeric_limits<double>::digits10)
950  << "#!geo" << endl
951  << endl
952  << "mesh" << endl
953  << version << endl
954  << "header" << endl
955  << " dimension\t" << dim << endl;
956  if (sys_coord_name != "cartesian") {
957  out << " coordinate_system " << sys_coord_name << endl;
958  }
959  if (order != 1) {
960  out << " order\t\t" << order << endl;
961  }
962  out << " nodes\t\t" << node.size() << endl;
963 
964  if (map_dim > 0) {
965  for (size_t variant = reference_element_first_variant_by_dimension(map_dim);
966  variant < reference_element_last_variant_by_dimension(map_dim); variant++) {
967  if (size_by_variant[variant] > 0) {
968  out << " " << geo_variant_name [variant] << "\t" << size_by_variant[variant] << endl;
969  }
970  }
971  }
972  if (map_dim >= 3 && size_by_variant[reference_element__t] != 0) {
973  out << " triangles\t" << size_by_variant[reference_element__t] << endl;
974  }
975  if (map_dim >= 3 && size_by_variant[reference_element__q] != 0) {
976  out << " quadrangles\t" << size_by_variant[reference_element__q] << endl;
977  }
978  if (map_dim >= 2 && size_by_variant[reference_element__e] != 0) {
979  out << " edges\t" << size_by_variant[reference_element__e] << endl;
980  }
981  out << "end header" << endl
982  << endl;
983  // nodes:
984  for (size_t inode = 0; inode < node.size(); inode++) {
985  for (size_t i = 0; i < dim; ++i) {
986  out << node[new2old_inode[inode]][i];
987  if (i+1 != dim) out << " ";
988  }
989  out << endl;
990  }
991  out << endl;
992  // elements:
993  for (size_t variant = reference_element_first_variant_by_dimension(map_dim);
994  variant < reference_element_last_variant_by_dimension(map_dim); variant++) {
995  for (size_t ie = 0; ie < element.size(); ie++) {
996  if (element[ie].variant() != variant) continue;
997  if (!do_upgrade) {
998  if (element[ie].name() != 'e' || element[ie].order() > 1) {
999  out << element[ie].name() << "\t";
1000  }
1001  if (element[ie].order() > 1) {
1002  out << "p" << element[ie].order() << " ";
1003  }
1004  for (size_t iloc = 0, nloc = element[ie].size(); iloc < nloc; iloc++) {
1005  out << element[ie][iloc];
1006  if (iloc+1 != nloc) out << " ";
1007  }
1008  } else {
1009  // upgrade: truncate inodes from pk to p1
1010  out << element[ie].name() << "\t";
1011  for (size_t iloc = 0, nloc = element[ie].n_vertex(); iloc < nloc; iloc++) {
1012  out << element[ie][iloc];
1013  if (iloc+1 != nloc) out << " ";
1014  }
1015  }
1016  out << endl;
1017  }
1018  }
1019  out << endl;
1020  // faces: no-empty when upgrade & map_dim >= 3
1021  for (size_t variant = reference_element_first_variant_by_dimension(2);
1022  variant < reference_element_last_variant_by_dimension(2); variant++) {
1023  for (size_t ie = 0; ie < face.size(); ie++) {
1024  if (face[ie].variant() != variant) continue;
1025  out << face[ie].name() << "\t";
1026  for (size_t iloc = 0, nloc = face[ie].n_vertex(); iloc < nloc; iloc++) {
1027  out << face[ie][iloc];
1028  if (iloc+1 != nloc) out << " ";
1029  }
1030  out << endl;
1031  }
1032  }
1033  // edges: no-empty when upgrade & map_dim >= 2
1034  for (size_t ie = 0, ne = edge.size(); ie < ne; ++ie) {
1035  out << "e\t" << old2new_inode[edge[ie][0]]
1036  << " " << old2new_inode[edge[ie][1]] << endl;
1037  }
1038  out << endl;
1039  //
1040  // 5.2. output domains
1041  //
1042  for (size_t d = 0; d <= 3; ++d) {
1043  if (!do_upgrade) {
1044  put_domain_noupgrade (out, map_dim, d, domain_map[d], phys, element);
1045  } else {
1046  put_domain_upgrade (out, map_dim, d, domain_map[d], phys, element,
1047  edge, edge_ball, face, face_ball);
1048  }
1049  }
1050 }
1051 void usage()
1052 {
1053  cerr << "msh2geo: usage:" << endl
1054  << "msh2geo [-[no]upgrade][-rz|-zr] in.msh > out.geo" << endl
1055  << "msh2geo [-[no]upgrade][-rz|-zr] < in.msh > out.geo" << endl;
1056  exit (1);
1057 }
1058 int main (int argc, char**argv) {
1059  // ----------------------------
1060  // scan the command line
1061  // ----------------------------
1062  string input_filename = "";
1063  std::string sys_coord_name = "cartesian";
1064  bool do_upgrade = true; // upgrade: put sides & truncate pk inodes to p1
1065  for (int i = 1; i < argc; i++) {
1066  if (strcmp (argv[i], "-rz") == 0) sys_coord_name = "rz";
1067  else if (strcmp (argv[i], "-zr") == 0) sys_coord_name = "zr";
1068  else if (strcmp (argv[i], "-cartesian") == 0) sys_coord_name = "cartesian";
1069  else if (strcmp (argv[i], "-upgrade") == 0) do_upgrade = true;
1070  else if (strcmp (argv[i], "-noupgrade") == 0) do_upgrade = false;
1071  else if (argv[i][0] == '-') {
1072  cerr << "field: unknown option `" << argv[i] << "'" << endl;
1073  usage();
1074  } else {
1075  input_filename = argv[i];
1076  }
1077  }
1078  if (input_filename == "") {
1079  msh2geo (cin, cout, sys_coord_name, do_upgrade);
1080  } else {
1081  ifstream fin (input_filename.c_str());
1082  if (!fin.good()) {
1083  cerr << "msh2geo: unable to read file \""<<input_filename<<"\"" << endl; exit (1);
1084  }
1085  msh2geo (fin, cout, sys_coord_name, do_upgrade);
1086  }
1087 }
see the edge page for the full documentation
see the hexahedron page for the full documentation
see the point page for the full documentation
see the prism page for the full documentation
see the quadrangle page for the full documentation
void setorientation(size_t orient)
Definition: msh2geo.cc:215
void setindex(size_t index)
Definition: msh2geo.cc:214
void setshift(size_t shift)
Definition: msh2geo.cc:216
see the geo_element page for the full documentation
Definition: geo_element.h:102
size_t face(size_t i) const
Definition: msh2geo.cc:247
array< geo_element_indirect, 12 > _edge
Definition: msh2geo.cc:259
const geo_element_indirect & edge_indirect(size_t i) const
Definition: msh2geo.cc:242
size_t variant() const
Definition: msh2geo.cc:239
geo_element_indirect & edge_indirect(size_t i)
Definition: msh2geo.cc:233
void setindex(size_t index)
Definition: msh2geo.cc:231
geo_element_indirect & face_indirect(size_t i)
Definition: msh2geo.cc:234
size_t edge(size_t i) const
Definition: msh2geo.cc:245
size_type face(size_type i) const
Definition: geo_element.h:210
size_type size() const
Definition: geo_element.h:168
size_t index() const
Definition: msh2geo.cc:237
void setorder(size_t order)
Definition: msh2geo.cc:230
size_t dimension() const
Definition: msh2geo.cc:240
void setname(char name)
Definition: msh2geo.cc:229
size_t edge_local_vertex(size_t iedg, size_t edg_iv) const
Definition: msh2geo.cc:263
size_t order() const
Definition: msh2geo.cc:236
size_t n_face() const
Definition: msh2geo.cc:248
void setgmshtype(size_t gmshtype)
Definition: msh2geo.cc:232
array< geo_element_indirect, 6 > _face
Definition: msh2geo.cc:260
size_t face_local_vertex(size_t iedg, size_t edg_iv) const
Definition: msh2geo.cc:276
const geo_element_indirect & face_indirect(size_t i) const
Definition: msh2geo.cc:243
size_t face_n_vertex(size_t loc_ifac) const
Definition: msh2geo.cc:288
variant_type variant() const
Definition: geo_element.h:161
size_t n_edge() const
Definition: msh2geo.cc:244
const geo_element_indirect & edge_indirect(size_type i) const
Definition: geo_element.h:193
size_t gmshtype() const
Definition: msh2geo.cc:238
char name() const
Definition: geo_element.h:169
size_t n_vertex() const
Definition: msh2geo.cc:241
size_type order() const
Definition: geo_element.h:162
void inplace_intersection(const index_set &b)
point_basic(T x0=T(), T x1=T(), T x2=T())
Definition: msh2geo.cc:136
std::array< T, 3 > base
Definition: msh2geo.cc:135
constexpr size_t reference_element__H
constexpr size_t reference_element__e
constexpr size_t reference_element__P
static size_t _reference_element_n_face_by_variant[reference_element__max_variant]
constexpr size_t reference_element__T
int shift_type
Definition: msh2geo.cc:191
int orientation_type
Definition: msh2geo.cc:190
constexpr size_t reference_element__t
constexpr size_t reference_element__q
constexpr size_t reference_element__max_variant
static size_t _reference_element_n_edge_by_variant[reference_element__max_variant]
const char _reference_element_name[reference_element__max_variant]
see the tetrahedron page for the full documentation
see the triangle page for the full documentation
#define assert_macro(ok_condition, message)
Definition: dis_macros.h:113
#define error_macro(message)
Definition: dis_macros.h:49
edge - reference element
const size_t dimension
Definition: edge.icc:64
const size_t n_vertex
Definition: edge.icc:66
Expr1::float_type T
Definition: field_expr.h:230
check_macro(expr1.have_homogeneous_space(Xh1), "dual(expr1,expr2); expr1 should have homogeneous space. HINT: use dual(interpolate(Xh, expr1),expr2)")
hexahedron - reference element
const size_t face[n_face][4]
Definition: hexahedron.icc:110
const size_t edge[n_edge][2]
Definition: hexahedron.icc:118
void numbering_Pk_dis_idof(size_t order, const array< size_t, reference_element__max_variant > &size_by_variant, const array< size_t, reference_element__max_variant+1 > &first_by_variant, const array< size_t, reference_element__max_variant > &loc_ndof_by_variant, const geo_element &K, vector< size_t > &dis_idof1)
Definition: msh2geo.cc:313
void usage()
Definition: msh2geo.cc:1051
int main(int argc, char **argv)
Definition: msh2geo.cc:1058
void put_domain_noupgrade(ostream &out, size_t dim, size_t dom_dim, const map< size_t, list< size_t > > &domain_map, map< size_t, string > &phys, const vector< geo_element > &element)
Definition: msh2geo.cc:396
void msh2geo(istream &in, ostream &out, string sys_coord_name, bool do_upgrade)
Definition: msh2geo.cc:513
void put_domain_upgrade(ostream &out, size_t dim, size_t dom_dim, const map< size_t, list< size_t > > &domain_map, map< size_t, string > &phys, const vector< geo_element > &element, const vector< geo_element_side > &edge, const vector< index_set > &edge_ball, const vector< geo_element_side > &face, const vector< index_set > &face_ball)
Definition: msh2geo.cc:434
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.
bool scatch(std::istream &in, const std::string &ch, bool full_match=true)
scatch: see the rheostream page for the full documentation
Definition: scatch.icc:44
void msh2geo_node_renum(vector< size_t > &element, size_t variant, size_t order)
prism - reference element
quadrangle - reference element
tetrahedron - reference element
triangle - reference element