GetFEM++  5.3
getfem_import.cc
1 /*===========================================================================
2 
3  Copyright (C) 2000-2017 Julien Pommier
4 
5  This file is a part of GetFEM++
6 
7  GetFEM++ is free software; you can redistribute it and/or modify it
8  under the terms of the GNU Lesser General Public License as published
9  by the Free Software Foundation; either version 3 of the License, or
10  (at your option) any later version along with the GCC Runtime Library
11  Exception either version 3.1 or (at your option) any later version.
12  This program is distributed in the hope that it will be useful, but
13  WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
14  or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
15  License and GCC Runtime Library Exception for more details.
16  You should have received a copy of the GNU Lesser General Public License
17  along with this program; if not, write to the Free Software Foundation,
18  Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA.
19 
20 ===========================================================================*/
21 
22 #include <iostream>
23 #include <iomanip>
24 #include <fstream>
25 
26 #include "getfem/getfem_mesh.h"
27 #include "getfem/getfem_import.h"
28 #include "getfem/getfem_regular_meshes.h"
29 
30 namespace getfem {
31 
32  /* mesh file from gmsh [http://www.geuz.org/gmsh/]*/
33 
34  struct gmsh_cv_info {
35  unsigned id, type, region;
36  bgeot::pgeometric_trans pgt;
37  std::vector<size_type> nodes;
38  void set_pgt() {
39  switch (type) {
40  case 1: { /* LINE */
41  pgt = bgeot::simplex_geotrans(1,1);
42  } break;
43  case 2: { /* TRIANGLE */
44  pgt = bgeot::simplex_geotrans(2,1);
45  } break;
46  case 3: { /* QUADRANGLE */
47  pgt = bgeot::parallelepiped_geotrans(2,1);
48  } break;
49  case 4: { /* TETRAHEDRON */
50  pgt = bgeot::simplex_geotrans(3,1);
51  } break;
52  case 5: { /* HEXAHEDRON */
53  pgt = bgeot::parallelepiped_geotrans(3,1);
54  } break;
55  case 6: { /* PRISM */
56  pgt = bgeot::prism_geotrans(3,1);
57  } break;
58  case 7: { /* PYRAMID */
59  pgt = bgeot::pyramid_QK_geotrans(1);
60  } break;
61  case 8: { /* 2ND ORDER LINE */
62  pgt = bgeot::simplex_geotrans(1,2);
63  } break;
64  case 9: { /* 2ND ORDER TRIANGLE */
65  pgt = bgeot::simplex_geotrans(2,2);
66  } break;
67  case 10: { /* 2ND ORDER QUADRANGLE */
68  pgt = bgeot::parallelepiped_geotrans(2,2);
69  } break;
70  case 11: { /* 2ND ORDER TETRAHEDRON (10-NODE) */
71  pgt = bgeot::simplex_geotrans(3,2);
72  } break;
73  case 12: { /* 2ND ORDER HEXAHEDRON (27-NODE) */
74  pgt = bgeot::parallelepiped_geotrans(3,2);
75  } break;
76  case 15: { /* POINT */
77  GMM_WARNING2("ignoring point element");
78  } break;
79  case 16: { /* INCOMPLETE 2ND ORDER QUADRANGLE (8-NODE) */
80  pgt = bgeot::Q2_incomplete_geotrans(2);
81  } break;
82  case 17: { /* INCOMPLETE 2ND ORDER HEXAHEDRON (20-NODE) */
83  pgt = bgeot::Q2_incomplete_geotrans(3);
84  } break;
85  case 26: { /* 3RD ORDER LINE */
86  pgt = bgeot::simplex_geotrans(1,3);
87  } break;
88  case 21: { /* 3RD ORDER TRIANGLE */
89  pgt = bgeot::simplex_geotrans(2,3);
90  } break;
91  case 23: { /* 4TH ORDER TRIANGLE */
92  pgt = bgeot::simplex_geotrans(2, 4);
93  } break;
94  case 27: { /* 4TH ORDER LINE */
95  pgt = bgeot::simplex_geotrans(1, 4);
96  } break;
97  default: { /* UNKNOWN .. */
98  /* higher order elements : to be done .. */
99  GMM_ASSERT1(false, "gmsh element type " << type << " is unknown.");
100  } break;
101  }
102  }
103 
104  void set_nb_nodes() {
105  /* Especially for the gmsh file format version 2.*/
106  switch (type) {
107  case 1: { /* LINE */
108  nodes.resize(2);
109  } break;
110  case 2: { /* TRIANGLE */
111  nodes.resize(3);
112  } break;
113  case 3: { /* QUADRANGLE */
114  nodes.resize(4);
115  } break;
116  case 4: { /* TETRAHEDRON */
117  nodes.resize(4);
118  } break;
119  case 5: { /* HEXAHEDRON */
120  nodes.resize(8);
121  } break;
122  case 6: { /* PRISM */
123  nodes.resize(6);
124  } break;
125  case 7: { /* PYRAMID */
126  nodes.resize(5);
127  } break;
128  case 8: { /* 2ND ORDER LINE */
129  nodes.resize(3);
130  } break;
131  case 9: { /* 2ND ORDER TRIANGLE */
132  nodes.resize(6);
133  } break;
134  case 10: { /* 2ND ORDER QUADRANGLE */
135  nodes.resize(9);
136  } break;
137  case 11: { /* 2ND ORDER TETRAHEDRON (10-NODE) */
138  nodes.resize(10);
139  } break;
140  case 12: { /* 2ND ORDER HEXAHEDRON (27-NODE) */
141  nodes.resize(27);
142  } break;
143  case 15: { /* POINT */
144  nodes.resize(1);
145  } break;
146  case 16: { /* INCOMPLETE 2ND ORDER QUADRANGLE (8-NODE) */
147  nodes.resize(8);
148  } break;
149  case 17: { /* INCOMPLETE 2ND ORDER HEXAHEDRON (20-NODE) */
150  nodes.resize(20);
151  } break;
152  case 26: { /* 3RD ORDER LINE */
153  nodes.resize(4);
154  } break;
155  case 21: { /* 3RD ORDER TRIANGLE */
156  nodes.resize(10);
157  } break;
158  case 23: { /* 4TH ORDER TRIANGLE */
159  nodes.resize(15);
160  } break;
161  case 27: { /* 4TH ORDER LINE */
162  nodes.resize(5);
163  } break;
164  default: { /* UNKNOWN .. */
165  /* higher order elements : to be done .. */
166  GMM_ASSERT1(false, "the gmsh element type " << type << " is unknown..");
167  } break;
168  }
169  }
170 
171  bool operator<(const gmsh_cv_info& other) const {
172  unsigned this_dim = (type == 15) ? 0 : pgt->dim();
173  unsigned other_dim = (other.type == 15) ? 0 : other.pgt->dim();
174  if (this_dim == other_dim) return region < other.region;
175  return this_dim > other_dim;
176  }
177  };
178 
179  std::map<std::string, size_type> read_region_names_from_gmsh_mesh_file(std::istream& f)
180  {
181  std::map<std::string, size_type> region_map;
182  bgeot::read_until(f, "$PhysicalNames");
183  size_type nb_regions;
184  f >> nb_regions;
185  size_type rt,ri;
186  std::string region_name;
187  for (size_type region_cnt=0; region_cnt < nb_regions; ++ region_cnt) {
188  f >> rt >> ri;
189  std::getline(f, region_name);
190  /* trim the string to the quote character front and back*/
191  size_t pos = region_name.find_first_of("\"");
192  if (pos != region_name.npos) {
193  region_name.erase(0, pos+1);
194  pos = region_name.find_last_of("\"");
195  region_name.erase(pos);
196  }
197  region_map[region_name] = ri;
198  }
199 
200  return region_map;
201  }
202 
203  /*
204  Format version 1 [for gmsh version < 2.0].
205  structure: $NOD list_of_nodes $ENDNOD $ELT list_of_elt $ENDELT
206 
207  Format version 2 [for gmsh version >= 2.0]. Modification of some keywords
208  and more than one tag is authorized for a specific region.
209 
210  structure: $Nodes list_of_nodes $EndNodes $Elements list_of_elt
211  $EndElements
212 
213  Lower dimensions elements in the regions of lower_dim_convex_rg will
214  be imported as independant convexes.
215 
216  If add_all_element_type is set to true, convexes with less dimension
217  than highest dimension pgt and are not part of other element's face will
218  be imported as independent convexes.
219 
220  for gmsh and gid meshes, the mesh nodes are always 3D, so for a 2D mesh
221  if remove_last_dimension == true the z-component of nodes will be removed
222  */
223  static void import_gmsh_mesh_file(std::istream& f, mesh& m, int deprecate=0,
224  std::map<std::string, size_type> *region_map=NULL,
225  std::set<size_type> *lower_dim_convex_rg=NULL,
226  bool add_all_element_type = false,
227  bool remove_last_dimension = true,
228  std::map<size_type, std::set<size_type>> *nodal_map = NULL,
229  bool remove_duplicated_nodes = true)
230  {
231  gmm::stream_standard_locale sl(f);
232  /* print general warning */
233  GMM_WARNING3(" All regions must have different number!");
234 
235  /* print deprecate warning */
236  if (deprecate!=0){
237  GMM_WARNING4("" << endl
238  << " deprecate: " << endl
239  << " static void" << endl
240  << " import_gmsh_mesh_file(std::istream& f,"
241  << " mesh& , int version)" << endl
242  << " replace with:" << endl
243  << " static void" << endl
244  << " import_gmsh_mesh_file(std::istream& f,"
245  << " mesh&)");
246  }
247 
248  /* read the version */
249  int version;
250  std::string header;
251  f >> header;
252  if (bgeot::casecmp(header,"$MeshFormat")==0)
253  f >> version;
254  else if (bgeot::casecmp(header,"$NOD")==0)
255  version = 1;
256  else
257  GMM_ASSERT1(false, "can't read Gmsh format: " << header);
258 
259  /* read the region names */
260  if (region_map != NULL) {
261  if (version == 2) {
262  *region_map = read_region_names_from_gmsh_mesh_file(f);
263  }
264  }
265  /* read the node list */
266  if (version == 2)
267  bgeot::read_until(f, "$Nodes"); /* Format version 2 */
268 
269  size_type nb_node;
270  f >> nb_node;
271  //cerr << "reading nodes..[nb=" << nb_node << "]\n";
272  std::map<size_type, size_type> msh_node_2_getfem_node;
273  for (size_type node_cnt=0; node_cnt < nb_node; ++node_cnt) {
274  size_type node_id;
275  base_node n(3); n[0]=n[1]=n[2]=0.0;
276  f >> node_id >> n[0] >> n[1] >> n[2];
277  msh_node_2_getfem_node[node_id] = m.add_point(n, 0.0, remove_duplicated_nodes);
278  }
279 
280  if (version == 2)
281  bgeot::read_until(f, "$Endnodes"); /* Format version 2 */
282  else
283  bgeot::read_until(f, "$ENDNOD");
284 
285  /* read the convexes */
286  if (version == 2)
287  bgeot::read_until(f, "$Elements"); /* Format version 2 */
288  else
289  bgeot::read_until(f, "$ELM");
290 
291  size_type nb_cv;
292  f >> nb_cv;
293  std::vector<gmsh_cv_info> cvlst; cvlst.reserve(nb_cv);
294  for (size_type cv=0; cv < nb_cv; ++cv) {
295  unsigned id, type, region;
296  unsigned dummy, cv_nb_nodes;
297 
298  if (version == 2) { /* Format version 2 */
299  unsigned nbtags, mesh_part;
300  f >> id >> type >> nbtags;
301  if (nbtags == 0 || nbtags > 3)
302  GMM_ASSERT1(false, "Number of tags " << nbtags
303  << " is not managed.");
304 
305  f >> region;
306  if (nbtags > 1) f >> dummy;
307  if (nbtags > 2) f >> mesh_part;
308  }
309  else
310  f >> id >> type >> region >> dummy >> cv_nb_nodes;
311 
312  id--; /* gmsh numbering starts at 1 */
313 
314  cvlst.push_back(gmsh_cv_info());
315  gmsh_cv_info &ci = cvlst.back();
316  ci.id = id; ci.type = type; ci.region = region;
317 
318 
319  if (version == 2) { /* For version 2 */
320  ci.set_nb_nodes();
321  cv_nb_nodes = unsigned(ci.nodes.size());
322  }
323  else
324  ci.nodes.resize(cv_nb_nodes);
325 
326  // cout << "cv_nb_nodes = " << cv_nb_nodes << endl;
327 
328  for (size_type i=0; i < cv_nb_nodes; ++i) {
329  size_type j;
330  f >> j;
331  std::map<size_type, size_type>::iterator
332  it = msh_node_2_getfem_node.find(j);
333  GMM_ASSERT1(it != msh_node_2_getfem_node.end(),
334  "Invalid node ID " << j << " in gmsh element "
335  << (ci.id + 1));
336  ci.nodes[i] = it->second;
337  }
338  if(ci.type != 15) ci.set_pgt();
339  // Reordering nodes for certain elements (should be completed ?)
340  // http://www.geuz.org/gmsh/doc/texinfo/gmsh.html#Node-ordering
341  std::vector<size_type> tmp_nodes(ci.nodes);
342  switch(ci.type) {
343  case 3 : {
344  ci.nodes[2] = tmp_nodes[3];
345  ci.nodes[3] = tmp_nodes[2];
346  } break;
347  case 5 : { /* First order hexaedron */
348  //ci.nodes[0] = tmp_nodes[0];
349  //ci.nodes[1] = tmp_nodes[1];
350  ci.nodes[2] = tmp_nodes[3];
351  ci.nodes[3] = tmp_nodes[2];
352  //ci.nodes[4] = tmp_nodes[4];
353  //ci.nodes[5] = tmp_nodes[5];
354  ci.nodes[6] = tmp_nodes[7];
355  ci.nodes[7] = tmp_nodes[6];
356  } break;
357  case 7 : { /* first order pyramid */
358  //ci.nodes[0] = tmp_nodes[0];
359  ci.nodes[1] = tmp_nodes[2];
360  ci.nodes[2] = tmp_nodes[1];
361  // ci.nodes[3] = tmp_nodes[3];
362  // ci.nodes[4] = tmp_nodes[4];
363  } break;
364  case 8 : { /* Second order line */
365  //ci.nodes[0] = tmp_nodes[0];
366  ci.nodes[1] = tmp_nodes[2];
367  ci.nodes[2] = tmp_nodes[1];
368  } break;
369  case 9 : { /* Second order triangle */
370  //ci.nodes[0] = tmp_nodes[0];
371  ci.nodes[1] = tmp_nodes[3];
372  ci.nodes[2] = tmp_nodes[1];
373  ci.nodes[3] = tmp_nodes[5];
374  //ci.nodes[4] = tmp_nodes[4];
375  ci.nodes[5] = tmp_nodes[2];
376  } break;
377  case 10 : { /* Second order quadrangle */
378  //ci.nodes[0] = tmp_nodes[0];
379  ci.nodes[1] = tmp_nodes[4];
380  ci.nodes[2] = tmp_nodes[1];
381  ci.nodes[3] = tmp_nodes[7];
382  ci.nodes[4] = tmp_nodes[8];
383  //ci.nodes[5] = tmp_nodes[5];
384  ci.nodes[6] = tmp_nodes[3];
385  ci.nodes[7] = tmp_nodes[6];
386  ci.nodes[8] = tmp_nodes[2];
387  } break;
388  case 11: { /* Second order tetrahedron */
389  //ci.nodes[0] = tmp_nodes[0];
390  ci.nodes[1] = tmp_nodes[4];
391  ci.nodes[2] = tmp_nodes[1];
392  ci.nodes[3] = tmp_nodes[6];
393  ci.nodes[4] = tmp_nodes[5];
394  ci.nodes[5] = tmp_nodes[2];
395  ci.nodes[6] = tmp_nodes[7];
396  ci.nodes[7] = tmp_nodes[9];
397  //ci.nodes[8] = tmp_nodes[8];
398  ci.nodes[9] = tmp_nodes[3];
399  } break;
400  case 12: { /* Second order hexahedron */
401  //ci.nodes[0] = tmp_nodes[0];
402  ci.nodes[1] = tmp_nodes[8];
403  ci.nodes[2] = tmp_nodes[1];
404  ci.nodes[3] = tmp_nodes[9];
405  ci.nodes[4] = tmp_nodes[20];
406  ci.nodes[5] = tmp_nodes[11];
407  ci.nodes[6] = tmp_nodes[3];
408  ci.nodes[7] = tmp_nodes[13];
409  ci.nodes[8] = tmp_nodes[2];
410  ci.nodes[9] = tmp_nodes[10];
411  ci.nodes[10] = tmp_nodes[21];
412  ci.nodes[11] = tmp_nodes[12];
413  ci.nodes[12] = tmp_nodes[22];
414  ci.nodes[13] = tmp_nodes[26];
415  ci.nodes[14] = tmp_nodes[23];
416  //ci.nodes[15] = tmp_nodes[15];
417  ci.nodes[16] = tmp_nodes[24];
418  ci.nodes[17] = tmp_nodes[14];
419  ci.nodes[18] = tmp_nodes[4];
420  ci.nodes[19] = tmp_nodes[16];
421  ci.nodes[20] = tmp_nodes[5];
422  ci.nodes[21] = tmp_nodes[17];
423  ci.nodes[22] = tmp_nodes[25];
424  ci.nodes[23] = tmp_nodes[18];
425  ci.nodes[24] = tmp_nodes[7];
426  ci.nodes[25] = tmp_nodes[19];
427  ci.nodes[26] = tmp_nodes[6];
428  } break;
429  case 16 : { /* Incomplete second order quadrangle */
430  //ci.nodes[0] = tmp_nodes[0];
431  ci.nodes[1] = tmp_nodes[4];
432  ci.nodes[2] = tmp_nodes[1];
433  ci.nodes[3] = tmp_nodes[7];
434  ci.nodes[4] = tmp_nodes[5];
435  ci.nodes[5] = tmp_nodes[3];
436  ci.nodes[6] = tmp_nodes[6];
437  ci.nodes[7] = tmp_nodes[2];
438  } break;
439  case 17: { /* Incomplete second order hexahedron */
440  //ci.nodes[0] = tmp_nodes[0];
441  ci.nodes[1] = tmp_nodes[8];
442  ci.nodes[2] = tmp_nodes[1];
443  ci.nodes[3] = tmp_nodes[9];
444  ci.nodes[4] = tmp_nodes[11];
445  ci.nodes[5] = tmp_nodes[3];
446  ci.nodes[6] = tmp_nodes[13];
447  ci.nodes[7] = tmp_nodes[2];
448  ci.nodes[8] = tmp_nodes[10];
449  ci.nodes[9] = tmp_nodes[12];
450  ci.nodes[10] = tmp_nodes[15];
451  ci.nodes[11] = tmp_nodes[14];
452  ci.nodes[12] = tmp_nodes[4];
453  ci.nodes[13] = tmp_nodes[16];
454  ci.nodes[14] = tmp_nodes[5];
455  ci.nodes[15] = tmp_nodes[17];
456  ci.nodes[16] = tmp_nodes[18];
457  ci.nodes[17] = tmp_nodes[7];
458  ci.nodes[18] = tmp_nodes[19];
459  ci.nodes[19] = tmp_nodes[6];
460  } break;
461  case 26 : { /* Third order line */
462  //ci.nodes[0] = tmp_nodes[0];
463  ci.nodes[1] = tmp_nodes[2];
464  ci.nodes[2] = tmp_nodes[3];
465  ci.nodes[3] = tmp_nodes[1];
466  } break;
467  case 21 : { /* Third order triangle */
468  //ci.nodes[0] = tmp_nodes[0];
469  ci.nodes[1] = tmp_nodes[3];
470  ci.nodes[2] = tmp_nodes[4];
471  ci.nodes[3] = tmp_nodes[1];
472  ci.nodes[4] = tmp_nodes[8];
473  ci.nodes[5] = tmp_nodes[9];
474  ci.nodes[6] = tmp_nodes[5];
475  //ci.nodes[7] = tmp_nodes[7];
476  ci.nodes[8] = tmp_nodes[6];
477  ci.nodes[9] = tmp_nodes[2];
478  } break;
479  case 23: { /* Fourth order triangle */
480  //ci.nodes[0] = tmp_nodes[0];
481  ci.nodes[1] = tmp_nodes[3];
482  ci.nodes[2] = tmp_nodes[4];
483  ci.nodes[3] = tmp_nodes[5];
484  ci.nodes[4] = tmp_nodes[1];
485  ci.nodes[5] = tmp_nodes[11];
486  ci.nodes[6] = tmp_nodes[12];
487  ci.nodes[7] = tmp_nodes[13];
488  ci.nodes[8] = tmp_nodes[6];
489  ci.nodes[9] = tmp_nodes[10];
490  ci.nodes[10] = tmp_nodes[14];
491  ci.nodes[11] = tmp_nodes[7];
492  ci.nodes[12] = tmp_nodes[9];
493  ci.nodes[13] = tmp_nodes[8];
494  ci.nodes[14] = tmp_nodes[2];
495  } break;
496  case 27: { /* Fourth order line */
497  //ci.nodes[0] = tmp_nodes[0];
498  ci.nodes[1] = tmp_nodes[2];
499  ci.nodes[2] = tmp_nodes[3];
500  ci.nodes[3] = tmp_nodes[4];
501  ci.nodes[4] = tmp_nodes[1];
502  } break;
503  }
504  }
505  nb_cv = cvlst.size();
506  if (cvlst.size()) {
507  std::sort(cvlst.begin(), cvlst.end());
508  if (cvlst.front().type == 15){
509  GMM_WARNING2("Only nodes defined in the mesh! No elements are added.");
510  return;
511  }
512 
513  unsigned N = cvlst.front().pgt->dim();
514  for (size_type cv=0; cv < nb_cv; ++cv) {
515  bool cvok = false;
516  gmsh_cv_info &ci = cvlst[cv];
517  bool is_node = (ci.type == 15);
518  unsigned ci_dim = (is_node) ? 0 : ci.pgt->dim();
519  //cout << "importing cv dim=" << int(ci.pgt->dim()) << " N=" << N
520  // << " region: " << ci.region << "\n";
521 
522  //main convex import
523  if (ci_dim == N) {
524  size_type ic = m.add_convex(ci.pgt, ci.nodes.begin());
525  cvok = true;
526  m.region(ci.region).add(ic);
527 
528  //convexes with lower dimensions
529  }
530  else {
531  //convex that lies within the regions of lower_dim_convex_rg
532  //is imported explicitly as a convex.
533  if (lower_dim_convex_rg != NULL &&
534  lower_dim_convex_rg->find(ci.region) != lower_dim_convex_rg->end() &&
535  !is_node){
536  size_type ic = m.add_convex(ci.pgt, ci.nodes.begin()); cvok = true;
537  m.region(ci.region).add(ic);
538  }
539  //find if the convex is part of a face of higher dimension convex
540  else{
541  bgeot::mesh_structure::ind_cv_ct ct = m.convex_to_point(ci.nodes[0]);
542  for (bgeot::mesh_structure::ind_cv_ct::const_iterator
543  it = ct.begin(); it != ct.end(); ++it) {
544  for (short_type face=0;
545  face < m.structure_of_convex(*it)->nb_faces(); ++face) {
546  if (m.is_convex_face_having_points(*it,face,
547  short_type(ci.nodes.size()),
548  ci.nodes.begin())) {
549  m.region(ci.region).add(*it,face);
550  cvok = true;
551  }
552  }
553  }
554  if (is_node && (nodal_map != NULL)) (*nodal_map)[ci.region].insert(ci.id);
555  //if the convex is not part of the face of others
556  if (!cvok)
557  {
558  if (is_node)
559  {
560  if (nodal_map == NULL){
561  GMM_WARNING2("gmsh import ignored a node id: "
562  << ci.id << " region :" << ci.region <<
563  " point is not added explicitly as an element.");
564  }
565  }
566  else if (add_all_element_type){
567  size_type ic = m.add_convex(ci.pgt, ci.nodes.begin());
568  m.region(ci.region).add(ic);
569  cvok = true;
570  }
571  else{
572  GMM_WARNING2("gmsh import ignored an element of type "
573  << bgeot::name_of_geometric_trans(ci.pgt) <<
574  " as it does not belong to the face of another element");
575  }
576  }
577  }
578  }
579  }
580  }
581  if (remove_last_dimension) maybe_remove_last_dimension(m);
582  }
583 
584  /* mesh file from GiD [http://gid.cimne.upc.es/]
585 
586  supports linear and quadratic elements (quadrilaterals, use 9(or 27)-noded elements)
587  */
588  static void import_gid_mesh_file(std::istream& f, mesh& m) {
589  gmm::stream_standard_locale sl(f);
590  /* read the node list */
591  size_type dim;
592  enum { LIN,TRI,QUAD,TETR, PRISM, HEX,BADELTYPE } eltype=BADELTYPE;
593  size_type nnode = 0;
594  std::map<size_type, size_type> msh_node_2_getfem_node;
595  std::vector<size_type> cv_nodes, getfem_cv_nodes;
596  bool nodes_done = false;
597  do {
598  if (!f.eof()) f >> std::ws;
599  if (f.eof() || !bgeot::read_until(f, "MESH")) break;
600  std::string selemtype;
601  f >> bgeot::skip("DIMENSION") >> dim
602  >> bgeot::skip("ELEMTYPE") >> std::ws
603  >> selemtype
604  >> bgeot::skip("NNODE") >> nnode;
605  if (bgeot::casecmp(selemtype, "linear")==0) { eltype = LIN; }
606  else if (bgeot::casecmp(selemtype, "triangle")==0) { eltype = TRI; }
607  else if (bgeot::casecmp(selemtype, "quadrilateral")==0) { eltype = QUAD; }
608  else if (bgeot::casecmp(selemtype, "tetrahedra")==0) { eltype = TETR; }
609  else if (bgeot::casecmp(selemtype, "prisma")==0) { eltype = PRISM; }
610  else if (bgeot::casecmp(selemtype, "hexahedra")==0) { eltype = HEX; }
611  else GMM_ASSERT1(false, "unknown element type '"<< selemtype << "'");
612  GMM_ASSERT1(!f.eof(), "File ended before coordinates");
613  f >> bgeot::skip("COORDINATES");
614  if (!nodes_done) {
615  dal::dynamic_array<base_node> gid_nodes;
616  dal::bit_vector gid_nodes_used;
617  do {
618  //cerr << "reading coordinates " << std::streamoff(f.tellg()) << "\n";
619  std::string ls;
620  f >> std::ws;
621  std::getline(f,ls);
622  if (bgeot::casecmp(ls, "END COORDINATES", 15)==0) break;
623  std::stringstream s; s << ls;
624  size_type id;
625  s >> id;
626 
627  gid_nodes[id].resize(dim); gid_nodes_used.add(id);
628  for (size_type i=0; i < dim; ++i) s >> gid_nodes[id][i];
629  //cerr << "ppoint " << id << ", " << n << endl;
630  } while (true);
631 
632  GMM_ASSERT1(gid_nodes_used.card() != 0, "no nodes in the mesh!");
633 
634  /* suppression of unused dimensions */
635  std::vector<bool> direction_useless(3,true);
636  base_node first_pt = gid_nodes[gid_nodes_used.first()];
637  for (dal::bv_visitor ip(gid_nodes_used); !ip.finished(); ++ip) {
638  for (size_type j=0; j < first_pt.size(); ++j) {
639  if (direction_useless[j] && (gmm::abs(gid_nodes[ip][j]-first_pt[j]) > 1e-13))
640  direction_useless[j] = false;
641  }
642  }
643  size_type dim2=0;
644  for (size_type j=0; j < dim; ++j) if (!direction_useless[j]) dim2++;
645  for (dal::bv_visitor ip(gid_nodes_used); !ip.finished(); ++ip) {
646  base_node n(dim2);
647  for (size_type j=0, cnt=0; j < dim; ++j) if (!direction_useless[j]) n[cnt++]=gid_nodes[ip][j];
648  msh_node_2_getfem_node[ip] = m.add_point(n);
649  }
650  }
651 
652  bgeot::read_until(f, "ELEMENTS");
653  bgeot::pgeometric_trans pgt = NULL;
654  std::vector<size_type> order(nnode); // ordre de GiD cf http://gid.cimne.upc.es/support/gid_11.subst#SEC160
655  for (size_type i=0; i < nnode; ++i) order[i]=i;
656  //cerr << "reading elements " << std::streamoff(f.tellg()) << "\n";
657  switch (eltype) {
658  case LIN: {
659  if (nnode == 2) pgt = bgeot::simplex_geotrans(1,1);
660  else if (nnode == 3) {
661  pgt = bgeot::simplex_geotrans(1,2); // A VERIFIER TOUT CA
662  std::swap(order[1],order[2]);
663  }
664  } break;
665  case TRI: {
666  if (nnode == 3) pgt = bgeot::simplex_geotrans(2,1);
667  else if (nnode == 6) { // validé
668  static size_type lorder[6] = {0,3,1,5,4,2};
669  pgt = bgeot::simplex_geotrans(2,2);
670  std::copy(lorder,lorder+nnode,order.begin());
671  }
672  } break;
673  case QUAD: {
674  if (nnode == 4) pgt = bgeot::parallelepiped_geotrans(2,1);
675  else if (nnode == 9) {
676  static size_type lorder[9] = {0,4,1, 7,8,5, 3,6,2};
677  pgt = bgeot::parallelepiped_geotrans(2,2);
678  std::copy(lorder,lorder+nnode,order.begin());
679  }
680  } break;
681  case TETR: {
682  if (nnode == 4) pgt = bgeot::simplex_geotrans(3,1);
683  else if (nnode == 10) { // validé
684  static size_type lorder[10] = {0,4,1, 7,8, 3, 6, 5, 9, 2};
685  pgt = bgeot::simplex_geotrans(3,2);
686  std::copy(lorder,lorder+nnode,order.begin());
687  }
688  } break;
689  case PRISM: {
690  if (nnode == 6) pgt = bgeot::prism_geotrans(3,1);
691  } break;
692  case HEX: {
693  if (nnode == 8) pgt = bgeot::parallelepiped_geotrans(3,1);
694  else if (nnode == 27) {
695  static size_type lorder[27] = {0,8,1, 12,21,13, 4,16,5,
696  11,20,9, 22,26,24, 19,25,17,
697  3,10,2, 15,23,14, 7,18,6};
698  pgt = bgeot::parallelepiped_geotrans(3,2);
699  std::copy(lorder,lorder+nnode,order.begin());
700  }
701  } break;
702  default: GMM_ASSERT1(false, ""); break;
703  }
704  GMM_ASSERT1(pgt != NULL, "unknown element type " << selemtype
705  << " with " << nnode << "nodes");
706  do {
707  std::string ls;
708  f >> std::ws;
709  std::getline(f,ls);
710  if (bgeot::casecmp(ls, "END ELEMENTS", 12)==0) break;
711  std::stringstream s; s << ls;
712  size_type cv_id;
713  s >> cv_id;
714  cv_nodes.resize(nnode);
715  for (size_type i=0; i < nnode; ++i) {
716  size_type j;
717  s >> j;
718  std::map<size_type, size_type>::iterator
719  it = msh_node_2_getfem_node.find(j);
720  GMM_ASSERT1(it != msh_node_2_getfem_node.end(),
721  "Invalid node ID " << j << " in GiD element " << cv_id);
722  cv_nodes[i] = it->second;
723  }
724  getfem_cv_nodes.resize(nnode);
725  for (size_type i=0; i < nnode; ++i) {
726  getfem_cv_nodes[i] = cv_nodes[order[i]];
727  }
728  //cerr << "elt " << cv_id << ", " << getfem_cv_nodes << endl;
729 
730  // envisager la "simplification" quand on a une transfo non
731  // lineaire mais que la destination est lineaire
732  m.add_convex(pgt, getfem_cv_nodes.begin());
733  } while (true);
734  } while (!f.eof());
735  }
736 
737  /* mesh file from ANSYS
738 
739  Supports solid structural elements stored with cdwrite in blocked format.
740  Use the following command in ANSYS for exporting the mesh:
741 
742  cdwrite,db,filename,cdb
743  */
744  static void import_cdb_mesh_file(std::istream& f, mesh& m,
745  size_type imat_filt=size_type(-1)) {
746 
747  std::map<size_type, size_type> cdb_node_2_getfem_node;
748  std::vector<size_type> getfem_cv_nodes;
749  std::vector<std::string> elt_types;
750  std::vector<size_type> elt_cnt;
751  std::vector<dal::bit_vector> regions;
752 
753  size_type pos, pos2;
754  std::string line;
755  while (true) {
756  std::getline(f,line);
757  pos = line.find_first_not_of(" ");
758  if (bgeot::casecmp(line.substr(pos,2),"ET") == 0) {
759  size_type itype;
760  std::string type_name;
761  pos = line.find_first_of(",")+1;
762  pos2 = line.find_first_of(",", pos);
763  itype = std::stol(line.substr(pos, pos2-pos));
764  pos = line.find_first_not_of(" ,\n\r\t", pos2);
765  pos2 = line.find_first_of(" ,\n\r\t", pos);
766  type_name = line.substr(pos, pos2-pos);
767  bool only_digits
768  = (type_name.find_first_not_of("0123456789") == std::string::npos);
769  const std::locale loc;
770  for (auto&& c : type_name) c = std::toupper(c, loc);
771 
772  if (elt_types.size() < itype+1)
773  elt_types.resize(itype+1);
774 
775  elt_types[itype] = "";
776  if (only_digits) {
777  size_type type_num = std::stol(type_name);
778  if (type_num == 42 || type_num == 82 ||
779  type_num == 182 || type_num == 183)
780  elt_types[itype] = "PLANE";
781  else if (type_num == 45 || type_num == 73 || type_num == 87 ||
782  type_num == 90 || type_num == 92 || type_num == 95 ||
783  type_num == 162 || type_num == 185 || type_num == 186 ||
784  type_num == 187 || type_num == 191)
785  elt_types[itype] = "SOLID";
786  else if (type_num == 89)
787  elt_types[itype] = "VISCO";
788  }
789  elt_types[itype].append(type_name);
790  }
791  else if (bgeot::casecmp(line.substr(pos,5),"KEYOP") == 0) {
792  size_type itype, knum, keyval;
793  pos = line.find_first_of(",")+1;
794  pos2 = line.find_first_of(",", pos);
795  itype = std::stol(line.substr(pos, pos2-pos));
796  pos = pos2+1;
797  pos2 = line.find_first_of(",", pos);
798  knum = std::stol(line.substr(pos+1, pos2-pos));
799  keyval = std::stol(line.substr(pos2+1));
800  if (knum == 1 && itype < elt_types.size() &&
801  elt_types[itype].size() == 7 &&
802  bgeot::casecmp(elt_types[itype].substr(0,7),"MESH200") == 0) {
803  std::stringstream ss("MESH200");
804  ss << "_" << keyval;
805  elt_types[itype] = ss.str();
806  }
807  }
808  else if (bgeot::casecmp(line.substr(pos,6),"NBLOCK") == 0)
809  break;
810  else if (f.eof())
811  return;
812  }
813  elt_cnt.resize(elt_types.size());
814 
815  //(3i8,6e20.13)
816  size_type fields1, fieldwidth1, fields2, fieldwidth2; // 3,8,6,20
817  { // "%8lu%*8u%*8u%20lf%20lf%20lf"
818  std::string fortran_fmt; // "(%lu%*[i]%lu,%lu%*[e,E]%lu.%*u)"
819  std::getline(f,fortran_fmt);
820  pos = fortran_fmt.find_first_of("(")+1;
821  pos2 = fortran_fmt.find_first_of("iI", pos);
822  fields1 = std::stol(fortran_fmt.substr(pos, pos2-pos));
823  pos = pos2+1;
824  pos2 = fortran_fmt.find_first_of(",", pos);
825  fieldwidth1 = std::stol(fortran_fmt.substr(pos, pos2-pos));
826  pos = pos2+1;
827  pos2 = fortran_fmt.find_first_of("eE", pos);
828  fields2 = std::stol(fortran_fmt.substr(pos, pos2-pos));
829  pos = pos2+1;
830  pos2 = fortran_fmt.find_first_of(".", pos);
831  fieldwidth2 = std::stol(fortran_fmt.substr(pos, pos2-pos));
832  GMM_ASSERT1(fields1 >= 1 && fields2 >= 3 ,
833  "Ansys mesh import routine requires NBLOCK entries with at least "
834  "1 integer field and 3 float number fields");
835  }
836 
837  base_node pt(3);
838  for (size_type i=0; i < size_type(-1); ++i) {
839  size_type nodeid;
840  std::getline(f,line);
841  if (line.compare(0,1,"N") == 0)
842  break;
843  // 1 0 0-3.0000000000000E+00 2.0000000000000E+00 1.0000000000000E+00
844  nodeid = std::stol(line.substr(0, fieldwidth1));
845  pos = fields1*fieldwidth1;
846  for (size_type j=0; j < 3; ++j, pos += fieldwidth2)
847  pt[j] = std::stod(line.substr(pos, fieldwidth2));
848 
849  cdb_node_2_getfem_node[nodeid] = m.add_point(pt, 0., false);
850  }
851 
852  while (bgeot::casecmp(line.substr(0,6),"EBLOCK") != 0) {
853  if (f.eof())
854  return;
855  std::getline(f,line);
856  }
857 
858 
859  //(19i8)
860  size_type fieldsno, fieldwidth; // 19,8
861  { // "%8lu%8lu%8lu%8lu%8lu%8lu%8lu%8lu"
862  std::string fortran_fmt;
863  std::getline(f,fortran_fmt);
864 
865  pos = fortran_fmt.find_first_of("(")+1;
866  pos2 = fortran_fmt.find_first_of("iI", pos);
867  fieldsno = std::stol(fortran_fmt.substr(pos, pos2-pos));
868  pos = pos2+1;
869  pos2 = fortran_fmt.find_first_of(")\n", pos);
870  fieldwidth = std::stol(fortran_fmt.substr(pos, pos2-pos));
871  GMM_ASSERT1(fieldsno == 19, "Ansys mesh import routine requires EBLOCK "
872  "entries with 19 fields");
873  }
874 
875  size_type II,JJ,KK,LL,MM,NN,OO,PP,QQ,RR,SS,TT,UU,VV,WW,XX,YY,ZZ,AA,BB;
876  for (size_type i=0; i < size_type(-1); ++i) {
877  GMM_ASSERT1(!f.eof(), "File ended before all elements could be read");
878  size_type imat, itype, nodesno(0);
879  std::getline(f,line);
880  {
881  long int ii = std::stol(line.substr(0,fieldwidth));
882  if (ii < 0)
883  break;
884  else
885  imat = size_type(ii);
886  }
887  itype = std::stol(line.substr(fieldwidth,fieldwidth));
888  nodesno = std::stol(line.substr(8*fieldwidth,fieldwidth));
889  line = line.substr(11*fieldwidth);
890 
891  if (imat_filt != size_type(-1) && imat != imat_filt) { // skip current element
892  if (nodesno > 8)
893  std::getline(f,line);
894  continue;
895  }
896 
897  if (nodesno >= 1) II = std::stol(line.substr(0,fieldwidth));
898  if (nodesno >= 2) JJ = std::stol(line.substr(1*fieldwidth,fieldwidth));
899  if (nodesno >= 3) KK = std::stol(line.substr(2*fieldwidth,fieldwidth));
900  if (nodesno >= 4) LL = std::stol(line.substr(3*fieldwidth,fieldwidth));
901  if (nodesno >= 5) MM = std::stol(line.substr(4*fieldwidth,fieldwidth));
902  if (nodesno >= 6) NN = std::stol(line.substr(5*fieldwidth,fieldwidth));
903  if (nodesno >= 7) OO = std::stol(line.substr(6*fieldwidth,fieldwidth));
904  if (nodesno >= 8) PP = std::stol(line.substr(7*fieldwidth,fieldwidth));
905  if (nodesno >= 9) {
906  std::getline(f,line);
907  if (nodesno >= 9) QQ = std::stol(line.substr(0,fieldwidth));
908  if (nodesno >= 10) RR = std::stol(line.substr(1*fieldwidth,fieldwidth));
909  if (nodesno >= 11) SS = std::stol(line.substr(2*fieldwidth,fieldwidth));
910  if (nodesno >= 12) TT = std::stol(line.substr(3*fieldwidth,fieldwidth));
911  if (nodesno >= 13) UU = std::stol(line.substr(4*fieldwidth,fieldwidth));
912  if (nodesno >= 14) VV = std::stol(line.substr(5*fieldwidth,fieldwidth));
913  if (nodesno >= 15) WW = std::stol(line.substr(6*fieldwidth,fieldwidth));
914  if (nodesno >= 16) XX = std::stol(line.substr(7*fieldwidth,fieldwidth));
915  if (nodesno >= 17) YY = std::stol(line.substr(8*fieldwidth,fieldwidth));
916  if (nodesno >= 18) ZZ = std::stol(line.substr(9*fieldwidth,fieldwidth));
917  if (nodesno >= 19) AA = std::stol(line.substr(10*fieldwidth,fieldwidth));
918  if (nodesno >= 20) BB = std::stol(line.substr(11*fieldwidth,fieldwidth));
919  }
920 
921  if (imat+1 > regions.size())
922  regions.resize(imat+1);
923 
924  if (nodesno == 3) {
925  // TODO MESH200_2, MESH200_3, MESH200_4
926  }
927  else if (nodesno == 4) {
928 
929  // assume MESH200_6 (4-node quadrilateral)
930  std::string eltname("MESH200_6");
931  if (elt_types.size() > itype && elt_types[itype].size() > 0)
932  eltname = elt_types[itype];
933 
934  if (eltname.compare("MESH200_6") == 0 ||
935  eltname.compare("PLANE42") == 0 ||
936  eltname.compare("PLANE182") == 0) {
937  getfem_cv_nodes.resize(4);
938  getfem_cv_nodes[0] = cdb_node_2_getfem_node[II];
939  getfem_cv_nodes[1] = cdb_node_2_getfem_node[JJ];
940  getfem_cv_nodes[2] = cdb_node_2_getfem_node[LL];
941  getfem_cv_nodes[3] = cdb_node_2_getfem_node[KK];
942  regions[imat].add(m.add_convex(bgeot::parallelepiped_geotrans(2,1),
943  getfem_cv_nodes.begin()));
944  if (itype < elt_cnt.size())
945  elt_cnt[itype] += 1;
946  }
947  else if (eltname.compare("MESH200_8") == 0 ||
948  eltname.compare("SOLID72") == 0) {
949  getfem_cv_nodes.resize(4);
950  getfem_cv_nodes[0] = cdb_node_2_getfem_node[II];
951  getfem_cv_nodes[1] = cdb_node_2_getfem_node[JJ];
952  getfem_cv_nodes[2] = cdb_node_2_getfem_node[KK];
953  getfem_cv_nodes[3] = cdb_node_2_getfem_node[LL];
954  regions[imat].add(m.add_convex(bgeot::simplex_geotrans(3,1),
955  getfem_cv_nodes.begin()));
956  if (itype < elt_cnt.size())
957  elt_cnt[itype] += 1;
958  }
959  }
960  else if (nodesno == 6) {
961  // TODO MESH200_5
962  }
963  else if (nodesno == 8) {
964 
965  // assume MESH200_10
966  std::string eltname("MESH200_10");
967  if (elt_types.size() > itype && elt_types[itype].size() > 0)
968  eltname = elt_types[itype];
969 
970  if (eltname.compare("MESH200_7") == 0 ||
971  eltname.compare("PLANE82") == 0 ||
972  eltname.compare("PLANE183") == 0) {
973  getfem_cv_nodes.resize(8);
974  getfem_cv_nodes[0] = cdb_node_2_getfem_node[II];
975  getfem_cv_nodes[1] = cdb_node_2_getfem_node[MM];
976  getfem_cv_nodes[2] = cdb_node_2_getfem_node[JJ];
977  getfem_cv_nodes[3] = cdb_node_2_getfem_node[PP];
978  getfem_cv_nodes[4] = cdb_node_2_getfem_node[NN];
979  getfem_cv_nodes[5] = cdb_node_2_getfem_node[LL];
980  getfem_cv_nodes[6] = cdb_node_2_getfem_node[OO];
981  getfem_cv_nodes[7] = cdb_node_2_getfem_node[KK];
982  regions[imat].add(m.add_convex(bgeot::Q2_incomplete_geotrans(2),
983  getfem_cv_nodes.begin()));
984  if (itype < elt_cnt.size())
985  elt_cnt[itype] += 1;
986  }
987  else if (eltname.compare("MESH200_10") == 0 ||
988  eltname.compare("SOLID45") == 0 ||
989  eltname.compare("SOLID185") == 0) {
990  if (KK == LL && OO == PP) {
991  if (MM == NN && NN == OO) { // 4-node tetrahedral
992  getfem_cv_nodes.resize(4);
993  getfem_cv_nodes[0] = cdb_node_2_getfem_node[II];
994  getfem_cv_nodes[1] = cdb_node_2_getfem_node[JJ];
995  getfem_cv_nodes[2] = cdb_node_2_getfem_node[KK];
996  getfem_cv_nodes[3] = cdb_node_2_getfem_node[MM];
997  regions[imat].add(m.add_convex(bgeot::simplex_geotrans(3,1),
998  getfem_cv_nodes.begin()));
999  }
1000  else { // 6-node prism
1001  getfem_cv_nodes.resize(6);
1002  getfem_cv_nodes[0] = cdb_node_2_getfem_node[II];
1003  getfem_cv_nodes[1] = cdb_node_2_getfem_node[JJ];
1004  getfem_cv_nodes[2] = cdb_node_2_getfem_node[KK];
1005  getfem_cv_nodes[3] = cdb_node_2_getfem_node[MM];
1006  getfem_cv_nodes[4] = cdb_node_2_getfem_node[NN];
1007  getfem_cv_nodes[5] = cdb_node_2_getfem_node[OO];
1008  regions[imat].add(m.add_convex(bgeot::prism_geotrans(3,1),
1009  getfem_cv_nodes.begin()));
1010  }
1011  }
1012  else { // 8-node hexahedral
1013  getfem_cv_nodes.resize(8);
1014  getfem_cv_nodes[0] = cdb_node_2_getfem_node[II];
1015  getfem_cv_nodes[1] = cdb_node_2_getfem_node[JJ];
1016  getfem_cv_nodes[2] = cdb_node_2_getfem_node[LL];
1017  getfem_cv_nodes[3] = cdb_node_2_getfem_node[KK];
1018  getfem_cv_nodes[4] = cdb_node_2_getfem_node[MM];
1019  getfem_cv_nodes[5] = cdb_node_2_getfem_node[NN];
1020  getfem_cv_nodes[6] = cdb_node_2_getfem_node[PP];
1021  getfem_cv_nodes[7] = cdb_node_2_getfem_node[OO];
1022  regions[imat].add(m.add_convex(bgeot::parallelepiped_geotrans(3,1),
1023  getfem_cv_nodes.begin()));
1024  }
1025  if (itype < elt_cnt.size())
1026  elt_cnt[itype] += 1;
1027  }
1028  }
1029  else if (nodesno == 10) {
1030 
1031  // assume MESH200_9 (10-node tetrahedral)
1032  std::string eltname("MESH200_9");
1033  if (elt_types.size() > itype && elt_types[itype].size() > 0)
1034  eltname = elt_types[itype];
1035 
1036  if (eltname.compare("MESH200_9") == 0 ||
1037  eltname.compare("SOLID87") == 0 ||
1038  eltname.compare("SOLID92") == 0 ||
1039  eltname.compare("SOLID162") == 0 ||
1040  eltname.compare("SOLID187") == 0) {
1041  getfem_cv_nodes.resize(10);
1042  getfem_cv_nodes[0] = cdb_node_2_getfem_node[II];
1043  getfem_cv_nodes[1] = cdb_node_2_getfem_node[MM];
1044  getfem_cv_nodes[2] = cdb_node_2_getfem_node[JJ];
1045  getfem_cv_nodes[3] = cdb_node_2_getfem_node[OO];
1046  getfem_cv_nodes[4] = cdb_node_2_getfem_node[NN];
1047  getfem_cv_nodes[5] = cdb_node_2_getfem_node[KK];
1048  getfem_cv_nodes[6] = cdb_node_2_getfem_node[PP];
1049  getfem_cv_nodes[7] = cdb_node_2_getfem_node[QQ];
1050  getfem_cv_nodes[8] = cdb_node_2_getfem_node[RR];
1051  getfem_cv_nodes[9] = cdb_node_2_getfem_node[LL];
1052  regions[imat].add(m.add_convex(bgeot::simplex_geotrans(3,2),
1053  getfem_cv_nodes.begin()));
1054  if (itype < elt_cnt.size())
1055  elt_cnt[itype] += 1;
1056  }
1057  }
1058  else if (nodesno == 20) { // # assume SOLID186/SOLID95
1059 
1060  // assume MESH200_11 (20-node hexahedral)
1061  std::string eltname("MESH200_11");
1062  if (elt_types.size() > itype && elt_types[itype].size() > 0)
1063  eltname = elt_types[itype];
1064 
1065  if (eltname.compare("MESH200_11") == 0 ||
1066  eltname.compare("VISCO89") == 0 ||
1067  eltname.compare("SOLID90") == 0 ||
1068  eltname.compare("SOLID95") == 0 ||
1069  eltname.compare("SOLID186") == 0 ||
1070  eltname.compare("SOLID191") == 0) {
1071  if (KK == LL && MM == NN && NN == OO && OO == PP) { // assume 10-node tetrahedral
1072  getfem_cv_nodes.resize(10);
1073  getfem_cv_nodes[0] = cdb_node_2_getfem_node[II];
1074  getfem_cv_nodes[1] = cdb_node_2_getfem_node[QQ];
1075  getfem_cv_nodes[2] = cdb_node_2_getfem_node[JJ];
1076  getfem_cv_nodes[3] = cdb_node_2_getfem_node[TT];
1077  getfem_cv_nodes[4] = cdb_node_2_getfem_node[RR];
1078  getfem_cv_nodes[5] = cdb_node_2_getfem_node[KK];
1079  getfem_cv_nodes[6] = cdb_node_2_getfem_node[YY];
1080  getfem_cv_nodes[7] = cdb_node_2_getfem_node[ZZ];
1081  getfem_cv_nodes[8] = cdb_node_2_getfem_node[AA];
1082  getfem_cv_nodes[9] = cdb_node_2_getfem_node[MM];
1083  regions[imat].add(m.add_convex(bgeot::simplex_geotrans(3,2),
1084  getfem_cv_nodes.begin()));
1085  if (itype < elt_cnt.size())
1086  elt_cnt[itype] += 1;
1087  } else if (MM == NN && NN == OO && OO == PP) { // assume 13-node pyramid
1088  getfem_cv_nodes.resize(13);
1089  getfem_cv_nodes[0] = cdb_node_2_getfem_node[II];
1090  getfem_cv_nodes[1] = cdb_node_2_getfem_node[QQ];
1091  getfem_cv_nodes[2] = cdb_node_2_getfem_node[JJ];
1092  getfem_cv_nodes[3] = cdb_node_2_getfem_node[TT];
1093  getfem_cv_nodes[4] = cdb_node_2_getfem_node[RR];
1094  getfem_cv_nodes[5] = cdb_node_2_getfem_node[LL];
1095  getfem_cv_nodes[6] = cdb_node_2_getfem_node[SS];
1096  getfem_cv_nodes[7] = cdb_node_2_getfem_node[KK];
1097  getfem_cv_nodes[8] = cdb_node_2_getfem_node[YY];
1098  getfem_cv_nodes[9] = cdb_node_2_getfem_node[ZZ];
1099  getfem_cv_nodes[10] = cdb_node_2_getfem_node[BB];
1100  getfem_cv_nodes[11] = cdb_node_2_getfem_node[AA];
1101  getfem_cv_nodes[12] = cdb_node_2_getfem_node[MM];
1102  regions[imat].add(m.add_convex(bgeot::pyramid_Q2_incomplete_geotrans(),
1103  getfem_cv_nodes.begin()));
1104  if (itype < elt_cnt.size())
1105  elt_cnt[itype] += 1;
1106 
1107  } else if (KK == LL && OO == PP) { // assume 15-node prism
1108  getfem_cv_nodes.resize(15);
1109  getfem_cv_nodes[0] = cdb_node_2_getfem_node[II];
1110  getfem_cv_nodes[1] = cdb_node_2_getfem_node[QQ];
1111  getfem_cv_nodes[2] = cdb_node_2_getfem_node[JJ];
1112  getfem_cv_nodes[3] = cdb_node_2_getfem_node[TT];
1113  getfem_cv_nodes[4] = cdb_node_2_getfem_node[RR];
1114  getfem_cv_nodes[5] = cdb_node_2_getfem_node[LL];
1115  getfem_cv_nodes[6] = cdb_node_2_getfem_node[YY];
1116  getfem_cv_nodes[7] = cdb_node_2_getfem_node[ZZ];
1117  getfem_cv_nodes[8] = cdb_node_2_getfem_node[AA];
1118  getfem_cv_nodes[9] = cdb_node_2_getfem_node[MM];
1119  getfem_cv_nodes[10] = cdb_node_2_getfem_node[UU];
1120  getfem_cv_nodes[11] = cdb_node_2_getfem_node[NN];
1121  getfem_cv_nodes[12] = cdb_node_2_getfem_node[XX];
1122  getfem_cv_nodes[13] = cdb_node_2_getfem_node[VV];
1123  getfem_cv_nodes[14] = cdb_node_2_getfem_node[OO];
1124  regions[imat].add(m.add_convex
1125  (bgeot::prism_incomplete_P2_geotrans(),
1126  getfem_cv_nodes.begin()));
1127  if (itype < elt_cnt.size())
1128  elt_cnt[itype] += 1;
1129  } else {
1130  getfem_cv_nodes.resize(20);
1131  getfem_cv_nodes[0] = cdb_node_2_getfem_node[II];
1132  getfem_cv_nodes[1] = cdb_node_2_getfem_node[QQ];
1133  getfem_cv_nodes[2] = cdb_node_2_getfem_node[JJ];
1134  getfem_cv_nodes[3] = cdb_node_2_getfem_node[TT];
1135  getfem_cv_nodes[4] = cdb_node_2_getfem_node[RR];
1136  getfem_cv_nodes[5] = cdb_node_2_getfem_node[LL];
1137  getfem_cv_nodes[6] = cdb_node_2_getfem_node[SS];
1138  getfem_cv_nodes[7] = cdb_node_2_getfem_node[KK];
1139  getfem_cv_nodes[8] = cdb_node_2_getfem_node[YY];
1140  getfem_cv_nodes[9] = cdb_node_2_getfem_node[ZZ];
1141  getfem_cv_nodes[10] = cdb_node_2_getfem_node[BB];
1142  getfem_cv_nodes[11] = cdb_node_2_getfem_node[AA];
1143  getfem_cv_nodes[12] = cdb_node_2_getfem_node[MM];
1144  getfem_cv_nodes[13] = cdb_node_2_getfem_node[UU];
1145  getfem_cv_nodes[14] = cdb_node_2_getfem_node[NN];
1146  getfem_cv_nodes[15] = cdb_node_2_getfem_node[XX];
1147  getfem_cv_nodes[16] = cdb_node_2_getfem_node[VV];
1148  getfem_cv_nodes[17] = cdb_node_2_getfem_node[PP];
1149  getfem_cv_nodes[18] = cdb_node_2_getfem_node[WW];
1150  getfem_cv_nodes[19] = cdb_node_2_getfem_node[OO];
1151  regions[imat].add(m.add_convex(bgeot::Q2_incomplete_geotrans(3),
1152  getfem_cv_nodes.begin()));
1153  if (itype < elt_cnt.size())
1154  elt_cnt[itype] += 1;
1155  }
1156  }
1157  }
1158  }
1159 
1160  int nonempty_regions=0;
1161  for (size_type i=0; i < regions.size(); ++i)
1162  if (regions[i].card() > 0)
1163  ++nonempty_regions;
1164 
1165  if (nonempty_regions > 1)
1166  for (size_type i=0; i < regions.size(); ++i)
1167  if (regions[i].card() > 0)
1168  m.region(i).add(regions[i]);
1169 
1170  for (size_type i=1; i < elt_types.size(); ++i)
1171  if (elt_cnt[i] > 0)
1172  cout << "Imported " << elt_cnt[i] << " " << elt_types[i] << " elements." << endl;
1173  cout << "Imported " << m.convex_index().card() << " elements in total." << endl;
1174 
1175  }
1176 
1177 
1178  static double round_to_nth_significant_number(double x, int ndec) {
1179  double p = 1.;
1180  double s = (x < 0 ? -1 : 1);
1181  double pdec = pow(10.,double(ndec));
1182  if (x == 0) return 0.;
1183  x = gmm::abs(x);
1184  while (x > 1) { x /= 10.0; p*=10; }
1185  while (x < 0.1) { x *= 10.0; p/=10; }
1186  //cerr << "x=" << x << ", p=" << p << ", pdec=" << pdec << "\n";
1187  x = s * (floor(x * pdec + 0.5) / pdec) * p;
1188  return x;
1189  }
1190 
1191 
1192  /* mesh file from noboite [http://www.distene.com/fr/corp/newsroom16.html] */
1193  static void import_noboite_mesh_file(std::istream& f, mesh& m) {
1194 
1195  using namespace std;
1196  gmm::stream_standard_locale sl(f);
1197 
1198  ofstream fichier_GiD("noboite_to_GiD.gid", ios::out | ios::trunc ); //déclaration du flux et ouverture du fichier
1199 
1200  fichier_GiD << "MESH dimension 3 ElemType Tetrahedra Nnode 4"<<endl;
1201 
1202 
1203  int i,NE,NP,ligne_debut_NP;
1204 
1205  /*NE: nombre d'elements (premier nombre du fichier .noboite)
1206  NP: nombre de points (deuxieme nombre du fichier .noboite)
1207  ligne_debut_NP: la ligne commence la liste des coordonnees des points dans le
1208  fichier .noboite*/
1209 
1210  f >> NE>>NP;
1211  ligne_debut_NP=NE*4/6+3;
1212 
1213  //passer 3 premiers lignes du fichier .noboite (la liste des elements commence a la
1214  //quatrieme ligne)
1215  string contenu;
1216  for (i=1; i<=ligne_debut_NP; i++){
1217  getline(f, contenu);
1218  }
1219 
1220 
1221  /*-----------------------------------------------------------------------
1222  Lire les coordonnees dans .noboite
1223  -----------------------------------------------------------------------*/
1224  fichier_GiD << "Coordinates" <<endl;
1225 
1226  for (i=1; i<=NP; i++){
1227  float coor_x,coor_y,coor_z;
1228 
1229  fichier_GiD << i<<" ";
1230 
1231  f>>coor_x >>coor_y >>coor_z;
1232  fichier_GiD<<coor_x<<" " <<coor_y <<" "<<coor_z <<endl;
1233 
1234  }
1235 
1236  fichier_GiD << "end coordinates" <<endl<<endl;
1237 
1238  /*-----------------------------------------------------------------------
1239  Lire les elements dans .noboite et ecrire au . gid
1240  ------------------------------------------------------------------------*/
1241 
1242  //revenir au debut du fichier . noboite, puis passer les trois premiere lignes
1243  f.seekg(0, ios::beg);
1244  for (i=1; i<=3; i++){
1245  getline(f, contenu);
1246  }
1247 
1248 
1249  fichier_GiD << "Elements" <<endl;
1250 
1251  for (i=1; i<=NE; i++){
1252  float elem_1,elem_2,elem_3,elem_4;
1253 
1254  fichier_GiD << i<<" ";
1255  f>>elem_1>>elem_2>>elem_3>>elem_4;
1256  fichier_GiD<<elem_1<<" " <<elem_2 <<" "<<elem_3<<" "<<elem_4<<" 1"<<endl;
1257 
1258  }
1259  fichier_GiD << "end elements" <<endl<<endl;
1260 
1261  if(fichier_GiD) // si l'ouverture a réussi
1262  {
1263  // instructions
1264  fichier_GiD.close(); // on referme le fichier
1265  }
1266  else // sinon
1267  cerr << "Erreur à l'ouverture !" << endl;
1268 
1269  if(f) // si l'ouverture a réussi
1270  {
1271  // instructions
1272  //f.close(); // on referme le fichier
1273  }
1274  else // sinon
1275  cerr << "Erreur à l'ouverture !" << endl;
1276 
1277  // appeler sunroutine import_gid_mesh_file
1278  //import_mesh(const std::string& "noboite_to_GiD.gid", mesh& msh)
1279  ifstream fichier1_GiD("noboite_to_GiD.gid", ios::in);
1280  import_gid_mesh_file(fichier1_GiD, m);
1281 
1282  // return 0;
1283  }
1284 
1285  /* mesh file from emc2 [http://pauillac.inria.fr/cdrom/prog/unix/emc2/eng.htm], am_fmt format
1286 
1287  (only triangular 2D meshes)
1288  */
1289  static void import_am_fmt_mesh_file(std::istream& f, mesh& m) {
1290  gmm::stream_standard_locale sl(f);
1291  /* read the node list */
1292  std::vector<size_type> tri;
1293  size_type nbs,nbt;
1294  base_node P(2);
1295  f >> nbs >> nbt; bgeot::read_until(f,"\n");
1296  tri.resize(nbt*3);
1297  for (size_type i=0; i < nbt*3; ++i) f >> tri[i];
1298  for (size_type j=0; j < nbs; ++j) {
1299  f >> P[0] >> P[1];
1300  cerr.precision(16);
1301  P[0]=round_to_nth_significant_number(P[0],6); // force 9.999999E-1 to be converted to 1.0
1302  P[1]=round_to_nth_significant_number(P[1],6);
1303  size_type jj = m.add_point(P);
1304  GMM_ASSERT1(jj == j, "ouch");
1305  }
1306  for (size_type i=0; i < nbt*3; i+=3)
1307  m.add_triangle(tri[i]-1,tri[i+1]-1,tri[i+2]-1);
1308  }
1309 
1310  /* mesh file from emc2 [http://pauillac.inria.fr/cdrom/prog/unix/emc2/eng.htm], am_fmt format
1311 
1312  triangular/quadrangular 2D meshes
1313  */
1314  static void import_emc2_mesh_file(std::istream& f, mesh& m) {
1315  gmm::stream_standard_locale sl(f);
1316  /* read the node list */
1317  std::vector<size_type> tri;
1318  size_type nbs=0,nbt=0,nbq=0,dummy;
1319  base_node P(2);
1320  bgeot::read_until(f,"Vertices");
1321  f >> nbs;
1322  for (size_type j=0; j < nbs; ++j) {
1323  f >> P[0] >> P[1] >> dummy;
1324  size_type jj = m.add_point(P);
1325  GMM_ASSERT1(jj == j, "ouch");
1326  }
1327  while (!f.eof()) {
1328  size_type ip[4];
1329  std::string ls;
1330  std::getline(f,ls);
1331  if (ls.find("Triangles")+1) {
1332  f >> nbt;
1333  for (size_type i=0; i < nbt; ++i) {
1334  f >> ip[0] >> ip[1] >> ip[2] >> dummy; ip[0]--; ip[1]--; ip[2]--;
1335  m.add_triangle(ip[0],ip[1],ip[2]);
1336  }
1337  } else if (ls.find("Quadrangles")+1) {
1338  f >> nbq;
1339  for (size_type i=0; i < nbq; ++i) {
1340  f >> ip[0] >> ip[1] >> ip[2] >> ip[3] >> dummy; ip[0]--; ip[1]--; ip[2]--; ip[3]--;
1341  m.add_parallelepiped(2, &ip[0]);
1342  }
1343  } else if (ls.find("End")+1) break;
1344  }
1345  }
1346 
1347  void import_mesh(const std::string& filename, const std::string& format,
1348  mesh& m) {
1349  m.clear();
1350  try {
1351 
1352  if (bgeot::casecmp(format,"structured")==0)
1353  { regular_mesh(m, filename); return; }
1354  else if (bgeot::casecmp(format,"structured_ball")==0)
1355  { regular_ball_mesh(m, filename); return; }
1356 
1357  std::ifstream f(filename.c_str());
1358  GMM_ASSERT1(f.good(), "can't open file " << filename);
1359  /* throw exceptions when an error occurs */
1360  f.exceptions(std::ifstream::badbit | std::ifstream::failbit);
1361  import_mesh(f, format,m);
1362  f.close();
1363  }
1364  catch (std::logic_error& exc) {
1365  m.clear();
1366  throw exc;
1367  }
1368  catch (std::ios_base::failure& exc) {
1369  m.clear();
1370  GMM_ASSERT1(false, "error while importing " << format
1371  << " mesh file \"" << filename << "\" : " << exc.what());
1372  }
1373  catch (std::runtime_error& exc) {
1374  m.clear();
1375  throw exc;
1376  }
1377  }
1378 
1379  void import_mesh_gmsh(std::istream& f, mesh &m,
1380  std::map<std::string, size_type> &region_map,
1381  bool remove_last_dimension,
1382  std::map<size_type, std::set<size_type>> *nodal_map,
1383  bool remove_duplicated_nodes)
1384  {
1385  import_gmsh_mesh_file(f, m, 0, &region_map, nullptr, false, remove_last_dimension, nodal_map,
1386  remove_duplicated_nodes);
1387  }
1388 
1389  void import_mesh_gmsh(std::istream& f, mesh& m,
1390  bool add_all_element_type,
1391  std::set<size_type> *lower_dim_convex_rg,
1392  std::map<std::string, size_type> *region_map,
1393  bool remove_last_dimension,
1394  std::map<size_type, std::set<size_type>> *nodal_map,
1395  bool remove_duplicated_nodes)
1396  {
1397  import_gmsh_mesh_file(f, m, 0, region_map, lower_dim_convex_rg, add_all_element_type,
1398  remove_last_dimension, nodal_map, remove_duplicated_nodes);
1399  }
1400 
1401  void import_mesh_gmsh(const std::string& filename, mesh& m,
1402  bool add_all_element_type,
1403  std::set<size_type> *lower_dim_convex_rg,
1404  std::map<std::string, size_type> *region_map,
1405  bool remove_last_dimension,
1406  std::map<size_type, std::set<size_type>> *nodal_map,
1407  bool remove_duplicated_nodes)
1408  {
1409  m.clear();
1410  try {
1411  std::ifstream f(filename.c_str());
1412  GMM_ASSERT1(f.good(), "can't open file " << filename);
1413  /* throw exceptions when an error occurs */
1414  f.exceptions(std::ifstream::badbit | std::ifstream::failbit);
1415  import_gmsh_mesh_file(f, m, 0, region_map, lower_dim_convex_rg, add_all_element_type,
1416  remove_last_dimension, nodal_map, remove_duplicated_nodes);
1417  f.close();
1418  }
1419  catch (std::logic_error& exc) {
1420  m.clear();
1421  throw exc;
1422  }
1423  catch (std::ios_base::failure& exc) {
1424  m.clear();
1425  GMM_ASSERT1(false, "error while importing " << "gmsh"
1426  << " mesh file \"" << filename << "\" : " << exc.what());
1427  }
1428  catch (std::runtime_error& exc) {
1429  m.clear();
1430  throw exc;
1431  }
1432  }
1433 
1434  void import_mesh_gmsh(const std::string& filename,
1435  mesh& m, std::map<std::string, size_type> &region_map,
1436  bool remove_last_dimension,
1437  std::map<size_type, std::set<size_type>> *nodal_map,
1438  bool remove_duplicated_nodes) {
1439  import_mesh_gmsh(filename, m, false, NULL, &region_map, remove_last_dimension, nodal_map,
1440  remove_duplicated_nodes);
1441  }
1442 
1443  void import_mesh(std::istream& f, const std::string& format,
1444  mesh& m) {
1445  if (bgeot::casecmp(format,"gmsh")==0)
1446  import_gmsh_mesh_file(f,m);
1447  else if (bgeot::casecmp(format,"gmshv2")==0)/* deprecate */
1448  import_gmsh_mesh_file(f,m,2);
1449  else if (bgeot::casecmp(format,"gid")==0)
1450  import_gid_mesh_file(f,m);
1451  else if (bgeot::casecmp(format,"noboite")==0)
1452  import_noboite_mesh_file(f,m);
1453  else if (bgeot::casecmp(format,"am_fmt")==0)
1454  import_am_fmt_mesh_file(f,m);
1455  else if (bgeot::casecmp(format,"emc2_mesh")==0)
1456  import_emc2_mesh_file(f,m);
1457  else if (bgeot::casecmp(format,"cdb")==0)
1458  import_cdb_mesh_file(f,m);
1459  else if (bgeot::casecmp(format.substr(0,4),"cdb:")==0) {
1460  size_type imat(-1);
1461  bool success(true);
1462  try {
1463  size_t sz;
1464  imat = std::stol(format.substr(4), &sz);
1465  success = (sz == format.substr(4).size() && imat != size_type(-1));
1466  } catch (const std::invalid_argument&) {
1467  success = false;
1468  } catch (const std::out_of_range&) {
1469  success = false;
1470  }
1471  if (success)
1472  import_cdb_mesh_file(f,m,imat);
1473  else GMM_ASSERT1(false, "cannot import "
1474  << format << " mesh type : wrong cdb mesh type input");
1475  }
1476  else GMM_ASSERT1(false, "cannot import "
1477  << format << " mesh type : unknown mesh type");
1478  }
1479 
1480  void import_mesh(const std::string& filename, mesh& msh) {
1481  size_type pos = filename.find_last_of(":");
1482  if (pos != std::string::npos)
1483  getfem::import_mesh(filename.substr(pos+1), filename.substr(0,pos), msh);
1484  else
1485  msh.read_from_file(filename);
1486  }
1487 
1488  void maybe_remove_last_dimension(mesh &m) {
1489  bool is_flat = true;
1490  unsigned N = m.dim(); if (N < 1) return;
1491  for (dal::bv_visitor i(m.points().index()); !i.finished(); ++i)
1492  if (m.points()[i][N-1] != 0) is_flat = 0;
1493  if (is_flat) {
1494  base_matrix M(N-1,N);
1495  for (unsigned i=0; i < N-1; ++i) M(i,i) = 1;
1496  m.transformation(M);
1497  }
1498  }
1499 
1500 }
getfem::mesh::add_parallelepiped
size_type add_parallelepiped(dim_type di, const ITER &ipts)
Add a parallelepiped to the mesh.
Definition: getfem_mesh.h:565
getfem::maybe_remove_last_dimension
void maybe_remove_last_dimension(mesh &msh)
for gmsh and gid meshes, the mesh nodes are always 3D, so for a 2D mesh the z-component of nodes shou...
Definition: getfem_import.cc:1488
getfem_mesh.h
Define a getfem::getfem_mesh object.
getfem::mesh::read_from_file
void read_from_file(const std::string &name)
Load the mesh from a file.
Definition: getfem_mesh.cc:613
getfem::mesh
Describe a mesh (collection of convexes (elements) and points).
Definition: getfem_mesh.h:95
bgeot::mesh_structure::structure_of_convex
pconvex_structure structure_of_convex(size_type ic) const
Return the pconvex_structure of the convex ic.
Definition: bgeot_mesh_structure.h:112
getfem::mesh::add_convex
size_type add_convex(bgeot::pgeometric_trans pgt, ITER ipts)
Add a convex to the mesh.
Definition: getfem_mesh.h:226
getfem::mesh::clear
void clear()
Erase the mesh.
Definition: getfem_mesh.cc:273
getfem_regular_meshes.h
Build regular meshes.
getfem_import.h
Import mesh files from various formats.
bgeot::mesh_structure::convex_to_point
const ind_cv_ct & convex_to_point(size_type ip) const
Return a container of the convexes attached to point ip.
Definition: bgeot_mesh_structure.h:166
bgeot::mesh_structure::is_convex_face_having_points
bool is_convex_face_having_points(size_type ic, short_type face_num, short_type nb, ITER pit) const
Return true if the face of the convex contains the given list of points.
Definition: bgeot_mesh_structure.h:269
getfem::regular_mesh
void regular_mesh(mesh &m, const std::string &st)
Build a regular mesh parametrized by the string st.
Definition: getfem_regular_meshes.cc:289
getfem::import_mesh
void import_mesh(const std::string &filename, const std::string &format, mesh &m)
imports a mesh file.
Definition: getfem_import.cc:1347
bgeot::size_type
size_t size_type
used as the common size type in the library
Definition: bgeot_poly.h:49
dal::dynamic_array
Dynamic Array.
Definition: dal_basic.h:47
getfem
GEneric Tool for Finite Element Methods.
Definition: getfem_assembling.h:48
bgeot::short_type
gmm::uint16_type short_type
used as the common short type integer in the library
Definition: bgeot_config.h:79
bgeot::name_of_geometric_trans
std::string name_of_geometric_trans(pgeometric_trans p)
Get the string name of a geometric transformation.
Definition: bgeot_geometric_trans.cc:1167
getfem::import_mesh_gmsh
void import_mesh_gmsh(const std::string &filename, mesh &m, std::map< std::string, size_type > &region_map, bool remove_last_dimension=true, std::map< size_type, std::set< size_type >> *nodal_map=NULL, bool remove_duplicated_nodes=true)
Import a mesh file in format generated by Gmsh.
Definition: getfem_import.cc:1434
getfem::mesh::add_triangle
size_type add_triangle(size_type a, size_type b, size_type c)
Add a triangle to the mesh, given the point id of its vertices.
Definition: getfem_mesh.cc:288
getfem::mesh::transformation
void transformation(const base_matrix &)
apply the given matrix transformation to each mesh node.
Definition: getfem_mesh.cc:255
getfem::mesh::region
const mesh_region region(size_type id) const
Return the region of index &#39;id&#39;.
Definition: getfem_mesh.h:414
getfem::regular_ball_mesh
void regular_ball_mesh(mesh &m, const std::string &st)
Build a regular mesh on a ball, parametrized by the string st.
Definition: getfem_regular_meshes.cc:356
bgeot::mesh_structure::convex_index
const dal::bit_vector & convex_index() const
Return the list of valid convex IDs.
Definition: bgeot_mesh_structure.h:91
bgeot::pgeometric_trans
std::shared_ptr< const bgeot::geometric_trans > pgeometric_trans
pointer type for a geometric transformation
Definition: bgeot_geometric_trans.h:183
getfem::read_region_names_from_gmsh_mesh_file
std::map< std::string, size_type > read_region_names_from_gmsh_mesh_file(std::istream &f)
for gmsh meshes, create table linking region name to region_id.
Definition: getfem_import.cc:179