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