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