dune-grid  2.8.0
gmshreader.hh
Go to the documentation of this file.
1 // -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
2 // vi: set et ts=4 sw=2 sts=2:
3 
4 #ifndef DUNE_GMSHREADER_HH
5 #define DUNE_GMSHREADER_HH
6 
7 #include <cstdarg>
8 #include <cstdio>
9 #include <cstring>
10 #include <fstream>
11 #include <iostream>
12 #include <map>
13 #include <memory>
14 #include <string>
15 #include <tuple>
16 #include <vector>
17 
18 #include <dune/common/exceptions.hh>
19 #include <dune/common/fvector.hh>
20 
21 #include <dune/geometry/type.hh>
22 
25 
26 namespace Dune
27 {
28 
36  {
42  };
43  };
44 
45  namespace {
46 
47  // arbitrary dimension, implementation is in specialization
48  template< int dimension, int dimWorld = dimension >
49  class GmshReaderQuadraticBoundarySegment
50  {
51  public:
52  // empty function since this class does not implement anything
53  static void registerFactory() {}
54  };
55 
56  // quadratic boundary segments in 1d
57  /*
58  Note the points
59 
60  (0) (alpha) (1)
61 
62  are mapped to the points in global coordinates
63 
64  p0 p2 p1
65 
66  alpha is determined automatically from the given points.
67  */
68  template< int dimWorld >
69  struct GmshReaderQuadraticBoundarySegment< 2, dimWorld >
70  : public Dune::BoundarySegment< 2, dimWorld >
71  {
72  typedef GmshReaderQuadraticBoundarySegment< 2, dimWorld > ThisType;
73  typedef typename Dune::BoundarySegment< 2, dimWorld > :: ObjectStreamType ObjectStreamType;
74  typedef Dune::FieldVector< double, dimWorld > GlobalVector;
75 
76  GmshReaderQuadraticBoundarySegment ( const GlobalVector &p0_, const GlobalVector &p1_, const GlobalVector &p2_)
77  : p0(p0_), p1(p1_), p2(p2_)
78  {
79  init();
80  }
81 
82  GmshReaderQuadraticBoundarySegment( ObjectStreamType& in )
83  {
84  // key is read before by the factory
85  const int bytes = sizeof(double)*dimWorld;
86  in.read( (char *) &p0[ 0 ], bytes );
87  in.read( (char *) &p1[ 0 ], bytes );
88  in.read( (char *) &p2[ 0 ], bytes );
89  init();
90  }
91 
92  static void registerFactory()
93  {
94  if( key() < 0 )
95  {
96  key() = Dune::BoundarySegment< 2, dimWorld >::template registerFactory< ThisType >();
97  }
98  }
99 
100  virtual GlobalVector operator() ( const Dune::FieldVector<double,1> &local ) const
101  {
102  GlobalVector y;
103  y = 0.0;
104  y.axpy((local[0]-alpha)*(local[0]-1.0)/alpha,p0);
105  y.axpy(local[0]*(local[0]-1.0)/(alpha*(alpha-1.0)),p1);
106  y.axpy(local[0]*(local[0]-alpha)/(1.0-alpha),p2);
107  return y;
108  }
109 
110  void backup( ObjectStreamType& out ) const
111  {
112  // backup key to identify object
113  out.write( (const char *) &key(), sizeof( int ) );
114  // backup data
115  const int bytes = sizeof(double)*dimWorld;
116  out.write( (const char*) &p0[ 0 ], bytes );
117  out.write( (const char*) &p1[ 0 ], bytes );
118  out.write( (const char*) &p2[ 0 ], bytes );
119  }
120 
121  protected:
122  void init()
123  {
124  GlobalVector d1 = p1;
125  d1 -= p0;
126  GlobalVector d2 = p2;
127  d2 -= p1;
128 
129  alpha=d1.two_norm()/(d1.two_norm()+d2.two_norm());
130  if (alpha<1E-6 || alpha>1-1E-6)
131  DUNE_THROW(Dune::IOError, "ration in quadratic boundary segment bad");
132  }
133 
134  static int& key() {
135  static int k = -1;
136  return k;
137  }
138 
139  private:
140  GlobalVector p0,p1,p2;
141  double alpha;
142  };
143 
144 
145  // quadratic boundary segments in 2d
146  /* numbering of points corresponding to gmsh:
147 
148  2
149 
150  5 4
151 
152  0 3 1
153 
154  Note: The vertices 3, 4, 5 are not necessarily at the edge midpoints but can
155  be placed with parameters alpha, beta , gamma at the following positions
156  in local coordinates:
157 
158 
159  2 = (0,1)
160 
161  5 = (0,beta) 4 = (1-gamma/sqrt(2),gamma/sqrt(2))
162 
163  0 = (0,0) 3 = (alpha,0) 1 = (1,0)
164 
165  The parameters alpha, beta, gamma are determined from the given vertices in
166  global coordinates.
167  */
168  template<>
169  class GmshReaderQuadraticBoundarySegment< 3, 3 >
170  : public Dune::BoundarySegment< 3 >
171  {
172  typedef GmshReaderQuadraticBoundarySegment< 3, 3 > ThisType;
173  typedef typename Dune::BoundarySegment< 3 > :: ObjectStreamType ObjectStreamType;
174  public:
175  GmshReaderQuadraticBoundarySegment (Dune::FieldVector<double,3> p0_, Dune::FieldVector<double,3> p1_,
176  Dune::FieldVector<double,3> p2_, Dune::FieldVector<double,3> p3_,
177  Dune::FieldVector<double,3> p4_, Dune::FieldVector<double,3> p5_)
178  : p0(p0_), p1(p1_), p2(p2_), p3(p3_), p4(p4_), p5(p5_)
179  {
180  init();
181  }
182 
183  GmshReaderQuadraticBoundarySegment( ObjectStreamType& in )
184  {
185  const int bytes = sizeof(double)*3;
186  in.read( (char *) &p0[ 0 ], bytes );
187  in.read( (char *) &p1[ 0 ], bytes );
188  in.read( (char *) &p2[ 0 ], bytes );
189  in.read( (char *) &p3[ 0 ], bytes );
190  in.read( (char *) &p4[ 0 ], bytes );
191  in.read( (char *) &p5[ 0 ], bytes );
192  init();
193  }
194 
195  static void registerFactory()
196  {
197  if( key() < 0 )
198  {
199  key() = Dune::BoundarySegment< 3 >::template registerFactory< ThisType >();
200  }
201  }
202 
203  virtual Dune::FieldVector<double,3> operator() (const Dune::FieldVector<double,2>& local) const
204  {
205  Dune::FieldVector<double,3> y;
206  y = 0.0;
207  y.axpy(phi0(local),p0);
208  y.axpy(phi1(local),p1);
209  y.axpy(phi2(local),p2);
210  y.axpy(phi3(local),p3);
211  y.axpy(phi4(local),p4);
212  y.axpy(phi5(local),p5);
213  return y;
214  }
215 
216  void backup( ObjectStreamType& out ) const
217  {
218  // backup key to identify object in factory
219  out.write( (const char*) &key(), sizeof( int ) );
220  // backup data
221  const int bytes = sizeof(double)*3;
222  out.write( (const char*) &p0[ 0 ], bytes );
223  out.write( (const char*) &p1[ 0 ], bytes );
224  out.write( (const char*) &p2[ 0 ], bytes );
225  out.write( (const char*) &p3[ 0 ], bytes );
226  out.write( (const char*) &p4[ 0 ], bytes );
227  out.write( (const char*) &p5[ 0 ], bytes );
228  }
229 
230  protected:
231  void init()
232  {
233  using std::sqrt;
234  sqrt2 = sqrt(2.0);
235  Dune::FieldVector<double,3> d1,d2;
236 
237  d1 = p3; d1 -= p0;
238  d2 = p1; d2 -= p3;
239  alpha=d1.two_norm()/(d1.two_norm()+d2.two_norm());
240  if (alpha<1E-6 || alpha>1-1E-6)
241  DUNE_THROW(Dune::IOError, "alpha in quadratic boundary segment bad");
242 
243  d1 = p5; d1 -= p0;
244  d2 = p2; d2 -= p5;
245  beta=d1.two_norm()/(d1.two_norm()+d2.two_norm());
246  if (beta<1E-6 || beta>1-1E-6)
247  DUNE_THROW(Dune::IOError, "beta in quadratic boundary segment bad");
248 
249  d1 = p4; d1 -= p1;
250  d2 = p2; d2 -= p4;
251  gamma=sqrt2*(d1.two_norm()/(d1.two_norm()+d2.two_norm()));
252  if (gamma<1E-6 || gamma>1-1E-6)
253  DUNE_THROW(Dune::IOError, "gamma in quadratic boundary segment bad");
254  }
255 
256  static int& key() {
257  static int k = -1;
258  return k;
259  }
260 
261  private:
262  // The six Lagrange basis function on the reference element
263  // for the points given above
264 
265  double phi0 (const Dune::FieldVector<double,2>& local) const
266  {
267  return (alpha*beta-beta*local[0]-alpha*local[1])*(1-local[0]-local[1])/(alpha*beta);
268  }
269  double phi3 (const Dune::FieldVector<double,2>& local) const
270  {
271  return local[0]*(1-local[0]-local[1])/(alpha*(1-alpha));
272  }
273  double phi1 (const Dune::FieldVector<double,2>& local) const
274  {
275  return local[0]*(gamma*local[0]-(sqrt2-gamma-sqrt2*alpha)*local[1]-alpha*gamma)/(gamma*(1-alpha));
276  }
277  double phi5 (const Dune::FieldVector<double,2>& local) const
278  {
279  return local[1]*(1-local[0]-local[1])/(beta*(1-beta));
280  }
281  double phi4 (const Dune::FieldVector<double,2>& local) const
282  {
283  return local[0]*local[1]/((1-gamma/sqrt2)*gamma/sqrt2);
284  }
285  double phi2 (const Dune::FieldVector<double,2>& local) const
286  {
287  return local[1]*(beta*(1-gamma/sqrt2)-local[0]*(beta-gamma/sqrt2)-local[1]*(1-gamma/sqrt2))/((1-gamma/sqrt2)*(beta-1));
288  }
289 
290  Dune::FieldVector<double,3> p0,p1,p2,p3,p4,p5;
291  double alpha,beta,gamma,sqrt2;
292  };
293 
294  } // end empty namespace
295 
297  template<typename GridType>
299  {
300  protected:
301  // private data
303  bool verbose;
308  // read buffer
309  char buf[512];
310  std::string fileName;
311  // exported data
314 
315  // static data
316  static const int dim = GridType::dimension;
317  static const int dimWorld = GridType::dimensionworld;
318  static_assert( (dimWorld <= 3), "GmshReader requires dimWorld <= 3." );
319 
320  // typedefs
321  typedef FieldVector< double, dimWorld > GlobalVector;
322 
323  // don't use something like
324  // readfile(file, 1, "%s\n", buf);
325  // to skip the rest of of the line -- that will only skip the next
326  // whitespace-separated word! Use skipline() instead.
327  void readfile(FILE * file, int cnt, const char * format,
328  void* t1, void* t2 = 0, void* t3 = 0, void* t4 = 0,
329  void* t5 = 0, void* t6 = 0, void* t7 = 0, void* t8 = 0,
330  void* t9 = 0, void* t10 = 0)
331  {
332  off_t pos = ftello(file);
333  int c = fscanf(file, format, t1, t2, t3, t4, t5, t6, t7, t8, t9, t10);
334  if (c != cnt)
335  DUNE_THROW(Dune::IOError, "Error parsing " << fileName << " "
336  "file pos " << pos
337  << ": Expected '" << format << "', only read " << c << " entries instead of " << cnt << ".");
338  }
339 
340  // skip over the rest of the line, including the terminating newline
341  void skipline(FILE * file)
342  {
343  int c;
344  do {
345  c = std::fgetc(file);
346  } while(c != '\n' && c != EOF);
347  }
348 
349  public:
350 
351  GmshReaderParser(Dune::GridFactory<GridType>& _factory, bool v, bool i) :
352  factory(_factory), verbose(v), insert_boundary_segments(i) {}
353 
354  std::vector<int> & boundaryIdMap()
355  {
357  }
358 
359  std::vector<int> & elementIndexMap()
360  {
362  }
363 
364  void read (const std::string& f)
365  {
366  if (verbose) std::cout << "Reading " << dim << "d Gmsh grid..." << std::endl;
367 
368  // open file name, we use C I/O
369  fileName = f;
370  FILE* file = fopen(fileName.c_str(),"rb");
371  if (file==0)
372  DUNE_THROW(Dune::IOError, "Could not open " << fileName);
373 
374  //=========================================
375  // Header: Read vertices into vector
376  // Check vertices that are needed
377  //=========================================
378 
381  element_count = 0;
382 
383  // process header
384  double version_number;
385  int file_type, data_size;
386 
387  readfile(file,1,"%s\n",buf);
388  if (strcmp(buf,"$MeshFormat")!=0)
389  DUNE_THROW(Dune::IOError, "expected $MeshFormat in first line");
390  readfile(file,3,"%lg %d %d\n",&version_number,&file_type,&data_size);
391  if( (version_number < 2.0) || (version_number > 2.2) )
392  DUNE_THROW(Dune::IOError, "can only read Gmsh version 2 files");
393  if (verbose) std::cout << "version " << version_number << " Gmsh file detected" << std::endl;
394  readfile(file,1,"%s\n",buf);
395  if (strcmp(buf,"$EndMeshFormat")!=0)
396  DUNE_THROW(Dune::IOError, "expected $EndMeshFormat");
397 
398  // node section
399  int number_of_nodes;
400 
401  readfile(file,1,"%s\n",buf);
402  if (strcmp(buf,"$Nodes")!=0)
403  DUNE_THROW(Dune::IOError, "expected $Nodes");
404  readfile(file,1,"%d\n",&number_of_nodes);
405  if (verbose) std::cout << "file contains " << number_of_nodes << " nodes" << std::endl;
406 
407  // read nodes
408  // The '+1' is due to the fact that gmsh numbers node starting from 1 rather than from 0
409  std::vector< GlobalVector > nodes( number_of_nodes+1 );
410  {
411  int id;
412  double x[ 3 ];
413  for( int i = 1; i <= number_of_nodes; ++i )
414  {
415  readfile(file,4, "%d %lg %lg %lg\n", &id, &x[ 0 ], &x[ 1 ], &x[ 2 ] );
416 
417  if (id > number_of_nodes) {
418  DUNE_THROW(Dune::IOError,
419  "Only dense sequences of node indices are currently supported (node index "
420  << id << " is invalid).");
421  }
422 
423  // just store node position
424  for( int j = 0; j < dimWorld; ++j )
425  nodes[ id ][ j ] = x[ j ];
426  }
427  readfile(file,1,"%s\n",buf);
428  if (strcmp(buf,"$EndNodes")!=0)
429  DUNE_THROW(Dune::IOError, "expected $EndNodes");
430  }
431 
432  // element section
433  readfile(file,1,"%s\n",buf);
434  if (strcmp(buf,"$Elements")!=0)
435  DUNE_THROW(Dune::IOError, "expected $Elements");
436  int number_of_elements;
437  readfile(file,1,"%d\n",&number_of_elements);
438  if (verbose) std::cout << "file contains " << number_of_elements << " elements" << std::endl;
439 
440  //=========================================
441  // Pass 1: Select and insert those vertices in the file that
442  // actually occur as corners of an element.
443  //=========================================
444 
445  off_t section_element_offset = ftello(file);
446  std::map<int,unsigned int> renumber;
447  for (int i=1; i<=number_of_elements; i++)
448  {
449  int id, elm_type, number_of_tags;
450  readfile(file,3,"%d %d %d ",&id,&elm_type,&number_of_tags);
451  for (int k=1; k<=number_of_tags; k++)
452  {
453  int blub;
454  readfile(file,1,"%d ",&blub);
455  // k == 1: physical entity (not used here)
456  // k == 2: elementary entity (not used here either)
457  // if version_number < 2.2:
458  // k == 3: mesh partition 0
459  // else
460  // k == 3: number of mesh partitions
461  // k => 4: mesh partition k-4
462  }
463  pass1HandleElement(file, elm_type, renumber, nodes);
464  }
465  if (verbose) std::cout << "number of real vertices = " << number_of_real_vertices << std::endl;
466  if (verbose) std::cout << "number of boundary elements = " << boundary_element_count << std::endl;
467  if (verbose) std::cout << "number of elements = " << element_count << std::endl;
468  readfile(file,1,"%s\n",buf);
469  if (strcmp(buf,"$EndElements")!=0)
470  DUNE_THROW(Dune::IOError, "expected $EndElements");
473 
474  //==============================================
475  // Pass 2: Insert boundary segments and elements
476  //==============================================
477 
478  fseeko(file, section_element_offset, SEEK_SET);
480  element_count = 0;
481  for (int i=1; i<=number_of_elements; i++)
482  {
483  int id, elm_type, number_of_tags;
484  readfile(file,3,"%d %d %d ",&id,&elm_type,&number_of_tags);
485  int physical_entity = -1;
486 
487  for (int k=1; k<=number_of_tags; k++)
488  {
489  int blub;
490  readfile(file,1,"%d ",&blub);
491  if (k==1) physical_entity = blub;
492  }
493  pass2HandleElement(file, elm_type, renumber, nodes, physical_entity);
494  }
495  readfile(file,1,"%s\n",buf);
496  if (strcmp(buf,"$EndElements")!=0)
497  DUNE_THROW(Dune::IOError, "expected $EndElements");
498 
499  fclose(file);
500  }
501 
507  void pass1HandleElement(FILE* file, const int elm_type,
508  std::map<int,unsigned int> & renumber,
509  const std::vector< GlobalVector > & nodes)
510  {
511  // some data about gmsh elements
512  const int nDofs[16] = {-1, 2, 3, 4, 4, 8, 6, 5, 3, 6, -1, 10, -1, -1, -1, 1};
513  const int nVertices[16] = {-1, 2, 3, 4, 4, 8, 6, 5, 2, 3, -1, 4, -1, -1, -1, 1};
514  const int elementDim[16] = {-1, 1, 2, 2, 3, 3, 3, 3, 1, 2, -1, 3, -1, -1, -1, 0};
515 
516  // test whether we support the element type
517  if ( not (elm_type > 0 && elm_type <= 15 // index in suitable range?
518  && (elementDim[elm_type] == dim || elementDim[elm_type] == (dim-1) ) ) ) // real element or boundary element?
519  {
520  skipline(file); // skip rest of line if element is unknown
521  return;
522  }
523 
524  // The format string for parsing is n times '%d' in a row
525  std::string formatString = "%d";
526  for (int i=1; i<nDofs[elm_type]; i++)
527  formatString += " %d";
528  formatString += "\n";
529 
530  // '10' is the largest number of dofs we may encounter in a .msh file
531  std::vector<int> elementDofs(10);
532 
533  readfile(file,nDofs[elm_type], formatString.c_str(),
534  &(elementDofs[0]),&(elementDofs[1]),&(elementDofs[2]),
535  &(elementDofs[3]),&(elementDofs[4]),&(elementDofs[5]),
536  &(elementDofs[6]),&(elementDofs[7]),&(elementDofs[8]),
537  &(elementDofs[9]));
538 
539  // insert each vertex if it hasn't been inserted already
540  for (int i=0; i<nVertices[elm_type]; i++)
541  if (renumber.find(elementDofs[i])==renumber.end())
542  {
543  renumber[elementDofs[i]] = number_of_real_vertices++;
544  factory.insertVertex(nodes[elementDofs[i]]);
545  }
546 
547  // count elements and boundary elements
548  if (elementDim[elm_type] == dim)
549  element_count++;
550  else
552 
553  }
554 
555 
556 
557  // generic-case: This is not supposed to be used at runtime.
558  template <class E, class V, class V2>
560  const V&,
561  const E&,
562  const V2&
563  )
564  {
565  DUNE_THROW(Dune::IOError, "tried to create a 3D boundary segment in a non-3D Grid");
566  }
567 
568  // 3d-case:
569  template <class E, class V>
571  const std::vector<FieldVector<double, 3> >& nodes,
572  const E& elementDofs,
573  const V& vertices
574  )
575  {
576  std::array<FieldVector<double,dimWorld>, 6> v;
577  for (int i=0; i<6; i++)
578  for (int j=0; j<dimWorld; j++)
579  v[i][j] = nodes[elementDofs[i]][j];
580 
581  BoundarySegment<dim,dimWorld>* newBoundarySegment
582  = (BoundarySegment<dim,dimWorld>*) new GmshReaderQuadraticBoundarySegment< 3, 3 >( v[0], v[1], v[2],
583  v[3], v[4], v[5] );
584 
585  factory.insertBoundarySegment( vertices,
586  std::shared_ptr<BoundarySegment<dim,dimWorld> >(newBoundarySegment) );
587  }
588 
589 
590 
595  virtual void pass2HandleElement(FILE* file, const int elm_type,
596  std::map<int,unsigned int> & renumber,
597  const std::vector< GlobalVector > & nodes,
598  const int physical_entity)
599  {
600  // some data about gmsh elements
601  const int nDofs[16] = {-1, 2, 3, 4, 4, 8, 6, 5, 3, 6, -1, 10, -1, -1, -1, 1};
602  const int nVertices[16] = {-1, 2, 3, 4, 4, 8, 6, 5, 2, 3, -1, 4, -1, -1, -1, 1};
603  const int elementDim[16] = {-1, 1, 2, 2, 3, 3, 3, 3, 1, 2, -1, 3, -1, -1, -1, 0};
604 
605  // test whether we support the element type
606  if ( not (elm_type > 0 && elm_type <= 15 // index in suitable range?
607  && (elementDim[elm_type] == dim || elementDim[elm_type] == (dim-1) ) ) ) // real element or boundary element?
608  {
609  skipline(file); // skip rest of line if element is unknown
610  return;
611  }
612 
613  // The format string for parsing is n times '%d' in a row
614  std::string formatString = "%d";
615  for (int i=1; i<nDofs[elm_type]; i++)
616  formatString += " %d";
617  formatString += "\n";
618 
619  // '10' is the largest number of dofs we may encounter in a .msh file
620  std::vector<int> elementDofs(10);
621 
622  readfile(file,nDofs[elm_type], formatString.c_str(),
623  &(elementDofs[0]),&(elementDofs[1]),&(elementDofs[2]),
624  &(elementDofs[3]),&(elementDofs[4]),&(elementDofs[5]),
625  &(elementDofs[6]),&(elementDofs[7]),&(elementDofs[8]),
626  &(elementDofs[9]));
627 
628  // correct differences between gmsh and Dune in the local vertex numbering
629  switch (elm_type)
630  {
631  case 3 : // 4-node quadrilateral
632  std::swap(elementDofs[2],elementDofs[3]);
633  break;
634  case 5 : // 8-node hexahedron
635  std::swap(elementDofs[2],elementDofs[3]);
636  std::swap(elementDofs[6],elementDofs[7]);
637  break;
638  case 7 : // 5-node pyramid
639  std::swap(elementDofs[2],elementDofs[3]);
640  break;
641  }
642 
643  // renumber corners to account for the explicitly given vertex
644  // numbering in the file
645  std::vector<unsigned int> vertices(nVertices[elm_type]);
646 
647  for (int i=0; i<nVertices[elm_type]; i++)
648  vertices[i] = renumber[elementDofs[i]];
649 
650  // If it is an element, insert it as such
651  if (elementDim[elm_type] == dim) {
652 
653  switch (elm_type)
654  {
655  case 1 : // 2-node line
656  factory.insertElement(Dune::GeometryTypes::line,vertices);
657  break;
658  case 2 : // 3-node triangle
659  factory.insertElement(Dune::GeometryTypes::triangle,vertices);
660  break;
661  case 3 : // 4-node quadrilateral
662  factory.insertElement(Dune::GeometryTypes::quadrilateral,vertices);
663  break;
664  case 4 : // 4-node tetrahedron
665  factory.insertElement(Dune::GeometryTypes::tetrahedron,vertices);
666  break;
667  case 5 : // 8-node hexahedron
668  factory.insertElement(Dune::GeometryTypes::hexahedron,vertices);
669  break;
670  case 6 : // 6-node prism
671  factory.insertElement(Dune::GeometryTypes::prism,vertices);
672  break;
673  case 7 : // 5-node pyramid
674  factory.insertElement(Dune::GeometryTypes::pyramid,vertices);
675  break;
676  case 9 : // 6-node triangle
677  factory.insertElement(Dune::GeometryTypes::triangle,vertices);
678  break;
679  case 11 : // 10-node tetrahedron
680  factory.insertElement(Dune::GeometryTypes::tetrahedron,vertices);
681  break;
682  }
683 
684  } else {
685  // it must be a boundary segment then
687 
688  switch (elm_type)
689  {
690  case 1 : // 2-node line
691  factory.insertBoundarySegment(vertices);
692  break;
693 
694  case 2 : // 3-node triangle
695  factory.insertBoundarySegment(vertices);
696  break;
697 
698  case 3 : // 4-node quadrilateral
699  factory.insertBoundarySegment(vertices);
700  break;
701 
702  case 15 : // 1-node point
703  factory.insertBoundarySegment(vertices);
704  break;
705 
706  case 8 : { // 3-node line
707  std::array<FieldVector<double,dimWorld>, 3> v;
708  for (int i=0; i<dimWorld; i++) {
709  v[0][i] = nodes[elementDofs[0]][i];
710  v[1][i] = nodes[elementDofs[2]][i]; // yes, the renumbering is intended!
711  v[2][i] = nodes[elementDofs[1]][i];
712  }
713  BoundarySegment<dim,dimWorld>* newBoundarySegment
714  = (BoundarySegment<dim,dimWorld>*) new GmshReaderQuadraticBoundarySegment< 2, dimWorld >(v[0], v[1], v[2]);
715  factory.insertBoundarySegment(vertices,
716  std::shared_ptr<BoundarySegment<dim,dimWorld> >(newBoundarySegment));
717  break;
718  }
719  case 9 : { // 6-node triangle
720  boundarysegment_insert(nodes, elementDofs, vertices);
721  break;
722  }
723  default: {
724  DUNE_THROW(Dune::IOError, "GmshReader does not support using element-type " << elm_type << " for boundary segments");
725  break;
726  }
727 
728  }
729 
730  }
731  }
732 
733  // count elements and boundary elements
734  if (elementDim[elm_type] == dim) {
736  element_count++;
737  } else {
740  }
741 
742  }
743 
744  };
745 
770  template<typename GridType>
772  {
774 
793  static void do_read(Dune::GridFactory<GridType> &factory,
794  const std::string &fileName,
795  std::vector<int>& boundarySegmentToPhysicalEntity,
796  std::vector<int>& elementToPhysicalEntity,
797  bool verbose, bool insertBoundarySegments)
798  {
799  // register boundary segment to boundary segment factory for possible load balancing
800  // this needs to be done on all cores since the type might not be known otherwise
801  GmshReaderQuadraticBoundarySegment< Grid::dimension, Grid::dimensionworld >::registerFactory();
802 
803 #ifndef NDEBUG
804  // check that this method is called on all cores
805  factory.comm().barrier();
806 #endif
807 
808  // create parse object and read grid on process 0
809  if (factory.comm().rank() == 0)
810  {
811  GmshReaderParser<Grid> parser(factory,verbose,insertBoundarySegments);
812  parser.read(fileName);
813 
814  boundarySegmentToPhysicalEntity = std::move(parser.boundaryIdMap());
815  elementToPhysicalEntity = std::move(parser.elementIndexMap());
816  }
817  else
818  {
819  boundarySegmentToPhysicalEntity = {};
820  elementToPhysicalEntity = {};
821  }
822  }
823 
825 
844  template<class T>
845  static T &discarded(T &&value) { return value; }
846 
847  struct DataArg {
848  std::vector<int> *data_ = nullptr;
849  DataArg(std::vector<int> &data) : data_(&data) {}
850  DataArg(const decltype(std::ignore)&) {}
851  DataArg() = default;
852  };
853 
854  struct DataFlagArg : DataArg {
855  bool flag_ = false;
856  using DataArg::DataArg;
857  DataFlagArg(bool flag) : flag_(flag) {}
858  };
859 
860  public:
861  typedef GridType Grid;
862 
869  static std::unique_ptr<Grid> read (const std::string& fileName, bool verbose = true, bool insertBoundarySegments=true)
870  {
871  // make a grid factory
872  Dune::GridFactory<Grid> factory;
873 
874  read(factory, fileName, verbose, insertBoundarySegments);
875 
876  return factory.createGrid();
877  }
878 
898  static std::unique_ptr<Grid> read (const std::string& fileName,
899  std::vector<int>& boundarySegmentToPhysicalEntity,
900  std::vector<int>& elementToPhysicalEntity,
901  bool verbose = true, bool insertBoundarySegments=true)
902  {
903  // make a grid factory
904  Dune::GridFactory<Grid> factory;
905 
906  do_read(factory, fileName, boundarySegmentToPhysicalEntity,
907  elementToPhysicalEntity, verbose, insertBoundarySegments);
908 
909  return factory.createGrid();
910  }
911 
913  static void read (Dune::GridFactory<Grid>& factory, const std::string& fileName,
914  bool verbose = true, bool insertBoundarySegments=true)
915  {
916  do_read(factory, fileName, discarded(std::vector<int>{}),
917  discarded(std::vector<int>{}), verbose, insertBoundarySegments);
918  }
919 
921 
944  static void read (Dune::GridFactory<Grid> &factory,
945  const std::string &fileName,
946  DataFlagArg boundarySegmentData,
947  DataArg elementData,
948  bool verbose=true)
949  {
950  do_read(factory, fileName,
951  boundarySegmentData.data_
952  ? *boundarySegmentData.data_ : discarded(std::vector<int>{}),
953  elementData.data_
954  ? *elementData.data_ : discarded(std::vector<int>{}),
955  verbose,
956  boundarySegmentData.flag_ || boundarySegmentData.data_);
957  }
958 
979  static void read (Dune::GridFactory<Grid>& factory,
980  const std::string& fileName,
981  std::vector<int>& boundarySegmentToPhysicalEntity,
982  std::vector<int>& elementToPhysicalEntity,
983  bool verbose, bool insertBoundarySegments)
984  {
985  do_read(factory, fileName, boundarySegmentToPhysicalEntity,
986  elementToPhysicalEntity, verbose, insertBoundarySegments);
987  }
988  };
989 
992 } // namespace Dune
993 
994 #endif
Base class for grid boundary segments of arbitrary geometry.
Include standard header files.
Definition: agrid.hh:58
static const int dimWorld
Definition: misc.hh:44
ALBERTA REAL_D GlobalVector
Definition: misc.hh:48
int renumber(const Dune::GeometryType &t, int i)
renumber VTK <-> Dune
Definition: common.hh:184
@ line
Definition: common.hh:132
@ pyramid
Definition: common.hh:139
@ quadrilateral
Definition: common.hh:135
@ tetrahedron
Definition: common.hh:136
@ prism
Definition: common.hh:138
@ hexahedron
Definition: common.hh:137
@ triangle
Definition: common.hh:133
Base class for classes implementing geometries of boundary segments.
Definition: boundarysegment.hh:92
Communication comm() const
Return the Communication used by the grid factory.
Definition: common/gridfactory.hh:295
Provide a generic factory class for unstructured grids.
Definition: common/gridfactory.hh:312
virtual std::unique_ptr< GridType > createGrid()
Finalize grid creation and hand over the grid.
Definition: common/gridfactory.hh:370
Options for read operation.
Definition: gmshreader.hh:36
GeometryOrder
Definition: gmshreader.hh:37
@ firstOrder
edges are straight lines.
Definition: gmshreader.hh:39
@ secondOrder
quadratic boundary approximation.
Definition: gmshreader.hh:41
dimension independent parts for GmshReaderParser
Definition: gmshreader.hh:299
void pass1HandleElement(FILE *file, const int elm_type, std::map< int, unsigned int > &renumber, const std::vector< GlobalVector > &nodes)
Process one element during the first pass through the list of all elements.
Definition: gmshreader.hh:507
static const int dimWorld
Definition: gmshreader.hh:317
Dune::GridFactory< GridType > & factory
Definition: gmshreader.hh:302
std::vector< int > & elementIndexMap()
Definition: gmshreader.hh:359
unsigned int number_of_real_vertices
Definition: gmshreader.hh:305
void boundarysegment_insert(const std::vector< FieldVector< double, 3 > > &nodes, const E &elementDofs, const V &vertices)
Definition: gmshreader.hh:570
GmshReaderParser(Dune::GridFactory< GridType > &_factory, bool v, bool i)
Definition: gmshreader.hh:351
int element_count
Definition: gmshreader.hh:307
void read(const std::string &f)
Definition: gmshreader.hh:364
void skipline(FILE *file)
Definition: gmshreader.hh:341
void readfile(FILE *file, int cnt, const char *format, void *t1, void *t2=0, void *t3=0, void *t4=0, void *t5=0, void *t6=0, void *t7=0, void *t8=0, void *t9=0, void *t10=0)
Definition: gmshreader.hh:327
std::vector< int > element_index_to_physical_entity
Definition: gmshreader.hh:313
virtual void pass2HandleElement(FILE *file, const int elm_type, std::map< int, unsigned int > &renumber, const std::vector< GlobalVector > &nodes, const int physical_entity)
Process one element during the second pass through the list of all elements.
Definition: gmshreader.hh:595
std::vector< int > & boundaryIdMap()
Definition: gmshreader.hh:354
static const int dim
Definition: gmshreader.hh:316
FieldVector< double, dimWorld > GlobalVector
Definition: gmshreader.hh:318
std::string fileName
Definition: gmshreader.hh:310
int boundary_element_count
Definition: gmshreader.hh:306
void boundarysegment_insert(const V &, const E &, const V2 &)
Definition: gmshreader.hh:559
bool verbose
Definition: gmshreader.hh:303
std::vector< int > boundary_id_to_physical_entity
Definition: gmshreader.hh:312
char buf[512]
Definition: gmshreader.hh:309
bool insert_boundary_segments
Definition: gmshreader.hh:304
Read Gmsh mesh file.
Definition: gmshreader.hh:772
static void read(Dune::GridFactory< Grid > &factory, const std::string &fileName, DataFlagArg boundarySegmentData, DataArg elementData, bool verbose=true)
read Gmsh file, possibly with data
Definition: gmshreader.hh:944
static void read(Dune::GridFactory< Grid > &factory, const std::string &fileName, bool verbose=true, bool insertBoundarySegments=true)
Definition: gmshreader.hh:913
GridType Grid
Definition: gmshreader.hh:861
static void read(Dune::GridFactory< Grid > &factory, const std::string &fileName, std::vector< int > &boundarySegmentToPhysicalEntity, std::vector< int > &elementToPhysicalEntity, bool verbose, bool insertBoundarySegments)
Read Gmsh file, possibly with data.
Definition: gmshreader.hh:979
static std::unique_ptr< Grid > read(const std::string &fileName, std::vector< int > &boundarySegmentToPhysicalEntity, std::vector< int > &elementToPhysicalEntity, bool verbose=true, bool insertBoundarySegments=true)
Read Gmsh file, possibly with data.
Definition: gmshreader.hh:898
static std::unique_ptr< Grid > read(const std::string &fileName, bool verbose=true, bool insertBoundarySegments=true)
Definition: gmshreader.hh:869
Provide a generic factory class for unstructured grids.