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