Dune Core Modules (unstable)

gmshreader.hh
1 // SPDX-FileCopyrightText: Copyright © DUNE Project contributors, see file LICENSE.md in module root
2 // SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception
3 // -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
4 // vi: set et ts=4 sw=2 sts=2:
5 
6 #ifndef DUNE_GMSHREADER_HH
7 #define DUNE_GMSHREADER_HH
8 
9 #include <cstdarg>
10 #include <cstdio>
11 #include <cstring>
12 #include <fstream>
13 #include <iostream>
14 #include <map>
15 #include <memory>
16 #include <string>
17 #include <tuple>
18 #include <vector>
19 #include <utility>
20 
22 #include <dune/common/fvector.hh>
23 
24 #include <dune/geometry/type.hh>
25 
28 
29 namespace Dune
30 {
31 
39  {
45  };
46  };
47 
48  namespace {
49 
50  // arbitrary dimension, implementation is in specialization
51  template< int dimension, int dimWorld = dimension >
52  class GmshReaderQuadraticBoundarySegment
53  {
54  public:
55  // empty function since this class does not implement anything
56  static void registerFactory() {}
57  };
58 
59  // quadratic boundary segments in 1d
60  /*
61  Note the points
62 
63  (0) (alpha) (1)
64 
65  are mapped to the points in global coordinates
66 
67  p0 p2 p1
68 
69  alpha is determined automatically from the given points.
70  */
71  template< int dimWorld >
72  struct GmshReaderQuadraticBoundarySegment< 2, dimWorld >
73  : public Dune::BoundarySegment< 2, dimWorld >
74  {
75  typedef GmshReaderQuadraticBoundarySegment< 2, dimWorld > ThisType;
76  typedef typename Dune::BoundarySegment< 2, dimWorld > :: ObjectStreamType ObjectStreamType;
77  typedef Dune::FieldVector< double, dimWorld > GlobalVector;
78 
79  GmshReaderQuadraticBoundarySegment ( const GlobalVector &p0_, const GlobalVector &p1_, const GlobalVector &p2_)
80  : p0(p0_), p1(p1_), p2(p2_)
81  {
82  init();
83  }
84 
85  GmshReaderQuadraticBoundarySegment( ObjectStreamType& in )
86  {
87  // key is read before by the factory
88  const int bytes = sizeof(double)*dimWorld;
89  in.read( (char *) &p0[ 0 ], bytes );
90  in.read( (char *) &p1[ 0 ], bytes );
91  in.read( (char *) &p2[ 0 ], bytes );
92  init();
93  }
94 
95  static void registerFactory()
96  {
97  if( key() < 0 )
98  {
99  key() = Dune::BoundarySegment< 2, dimWorld >::template registerFactory< ThisType >();
100  }
101  }
102 
103  virtual GlobalVector operator() ( const Dune::FieldVector<double,1> &local ) const
104  {
105  GlobalVector y;
106  y = 0.0;
107  y.axpy((local[0]-alpha)*(local[0]-1.0)/alpha,p0);
108  y.axpy(local[0]*(local[0]-1.0)/(alpha*(alpha-1.0)),p1);
109  y.axpy(local[0]*(local[0]-alpha)/(1.0-alpha),p2);
110  return y;
111  }
112 
113  void backup( ObjectStreamType& out ) const
114  {
115  // backup key to identify object
116  out.write( (const char *) &key(), sizeof( int ) );
117  // backup data
118  const int bytes = sizeof(double)*dimWorld;
119  out.write( (const char*) &p0[ 0 ], bytes );
120  out.write( (const char*) &p1[ 0 ], bytes );
121  out.write( (const char*) &p2[ 0 ], bytes );
122  }
123 
124  protected:
125  void init()
126  {
127  GlobalVector d1 = p1;
128  d1 -= p0;
129  GlobalVector d2 = p2;
130  d2 -= p1;
131 
132  alpha=d1.two_norm()/(d1.two_norm()+d2.two_norm());
133  if (alpha<1E-6 || alpha>1-1E-6)
134  DUNE_THROW(Dune::IOError, "ration in quadratic boundary segment bad");
135  }
136 
137  static int& key() {
138  static int k = -1;
139  return k;
140  }
141 
142  private:
143  GlobalVector p0,p1,p2;
144  double alpha;
145  };
146 
147 
148  // quadratic boundary segments in 2d
149  /* numbering of points corresponding to gmsh:
150 
151  2
152 
153  5 4
154 
155  0 3 1
156 
157  Note: The vertices 3, 4, 5 are not necessarily at the edge midpoints but can
158  be placed with parameters alpha, beta , gamma at the following positions
159  in local coordinates:
160 
161 
162  2 = (0,1)
163 
164  5 = (0,beta) 4 = (1-gamma/sqrt(2),gamma/sqrt(2))
165 
166  0 = (0,0) 3 = (alpha,0) 1 = (1,0)
167 
168  The parameters alpha, beta, gamma are determined from the given vertices in
169  global coordinates.
170  */
171  template<>
172  class GmshReaderQuadraticBoundarySegment< 3, 3 >
173  : public Dune::BoundarySegment< 3 >
174  {
175  typedef GmshReaderQuadraticBoundarySegment< 3, 3 > ThisType;
176  typedef typename Dune::BoundarySegment< 3 > :: ObjectStreamType ObjectStreamType;
177  public:
178  GmshReaderQuadraticBoundarySegment (Dune::FieldVector<double,3> p0_, Dune::FieldVector<double,3> p1_,
181  : p0(p0_), p1(p1_), p2(p2_), p3(p3_), p4(p4_), p5(p5_)
182  {
183  init();
184  }
185 
186  GmshReaderQuadraticBoundarySegment( ObjectStreamType& in )
187  {
188  const int bytes = sizeof(double)*3;
189  in.read( (char *) &p0[ 0 ], bytes );
190  in.read( (char *) &p1[ 0 ], bytes );
191  in.read( (char *) &p2[ 0 ], bytes );
192  in.read( (char *) &p3[ 0 ], bytes );
193  in.read( (char *) &p4[ 0 ], bytes );
194  in.read( (char *) &p5[ 0 ], bytes );
195  init();
196  }
197 
198  static void registerFactory()
199  {
200  if( key() < 0 )
201  {
202  key() = Dune::BoundarySegment< 3 >::template registerFactory< ThisType >();
203  }
204  }
205 
206  virtual Dune::FieldVector<double,3> operator() (const Dune::FieldVector<double,2>& local) const
207  {
209  y = 0.0;
210  y.axpy(phi0(local),p0);
211  y.axpy(phi1(local),p1);
212  y.axpy(phi2(local),p2);
213  y.axpy(phi3(local),p3);
214  y.axpy(phi4(local),p4);
215  y.axpy(phi5(local),p5);
216  return y;
217  }
218 
219  void backup( ObjectStreamType& out ) const
220  {
221  // backup key to identify object in factory
222  out.write( (const char*) &key(), sizeof( int ) );
223  // backup data
224  const int bytes = sizeof(double)*3;
225  out.write( (const char*) &p0[ 0 ], bytes );
226  out.write( (const char*) &p1[ 0 ], bytes );
227  out.write( (const char*) &p2[ 0 ], bytes );
228  out.write( (const char*) &p3[ 0 ], bytes );
229  out.write( (const char*) &p4[ 0 ], bytes );
230  out.write( (const char*) &p5[ 0 ], bytes );
231  }
232 
233  protected:
234  void init()
235  {
236  using std::sqrt;
237  sqrt2 = sqrt(2.0);
239 
240  d1 = p3; d1 -= p0;
241  d2 = p1; d2 -= p3;
242  alpha=d1.two_norm()/(d1.two_norm()+d2.two_norm());
243  if (alpha<1E-6 || alpha>1-1E-6)
244  DUNE_THROW(Dune::IOError, "alpha in quadratic boundary segment bad");
245 
246  d1 = p5; d1 -= p0;
247  d2 = p2; d2 -= p5;
248  beta=d1.two_norm()/(d1.two_norm()+d2.two_norm());
249  if (beta<1E-6 || beta>1-1E-6)
250  DUNE_THROW(Dune::IOError, "beta in quadratic boundary segment bad");
251 
252  d1 = p4; d1 -= p1;
253  d2 = p2; d2 -= p4;
254  gamma=sqrt2*(d1.two_norm()/(d1.two_norm()+d2.two_norm()));
255  if (gamma<1E-6 || gamma>1-1E-6)
256  DUNE_THROW(Dune::IOError, "gamma in quadratic boundary segment bad");
257  }
258 
259  static int& key() {
260  static int k = -1;
261  return k;
262  }
263 
264  private:
265  // The six Lagrange basis function on the reference element
266  // for the points given above
267 
268  double phi0 (const Dune::FieldVector<double,2>& local) const
269  {
270  return (alpha*beta-beta*local[0]-alpha*local[1])*(1-local[0]-local[1])/(alpha*beta);
271  }
272  double phi3 (const Dune::FieldVector<double,2>& local) const
273  {
274  return local[0]*(1-local[0]-local[1])/(alpha*(1-alpha));
275  }
276  double phi1 (const Dune::FieldVector<double,2>& local) const
277  {
278  return local[0]*(gamma*local[0]-(sqrt2-gamma-sqrt2*alpha)*local[1]-alpha*gamma)/(gamma*(1-alpha));
279  }
280  double phi5 (const Dune::FieldVector<double,2>& local) const
281  {
282  return local[1]*(1-local[0]-local[1])/(beta*(1-beta));
283  }
284  double phi4 (const Dune::FieldVector<double,2>& local) const
285  {
286  return local[0]*local[1]/((1-gamma/sqrt2)*gamma/sqrt2);
287  }
288  double phi2 (const Dune::FieldVector<double,2>& local) const
289  {
290  return local[1]*(beta*(1-gamma/sqrt2)-local[0]*(beta-gamma/sqrt2)-local[1]*(1-gamma/sqrt2))/((1-gamma/sqrt2)*(beta-1));
291  }
292 
293  Dune::FieldVector<double,3> p0,p1,p2,p3,p4,p5;
294  double alpha,beta,gamma,sqrt2;
295  };
296 
297  } // end empty namespace
298 
300  template<typename GridType>
302  {
303  protected:
304  // private data
306  bool verbose;
307  bool insert_boundary_segments;
308  unsigned int number_of_real_vertices;
309  int boundary_element_count;
310  int element_count;
311  // read buffer
312  char buf[512];
313  std::string fileName;
314  // exported data
315  std::vector<int> boundary_id_to_physical_entity;
316  std::vector<int> element_index_to_physical_entity;
317  std::vector<std::string> physical_entity_names;
318 
319  // static data
320  static const int dim = GridType::dimension;
321  static const int dimWorld = GridType::dimensionworld;
322  static_assert( (dimWorld <= 3), "GmshReader requires dimWorld <= 3." );
323 
324  // typedefs
326 
327  // don't use something like
328  // readfile(file, 1, "%s\n", buf);
329  // to skip the rest of of the line -- that will only skip the next
330  // whitespace-separated word! Use skipline() instead.
331  void readfile(FILE * file, int cnt, const char * format, ...)
332 #ifdef __GNUC__
333  __attribute__((format(scanf, 4, 5)))
334 #endif
335  {
336  std::va_list ap;
337  va_start(ap, format);
338  off_t pos = ftello(file);
339  int c = std::vfscanf(file, format, ap);
340  va_end(ap);
341  if (c != cnt)
342  DUNE_THROW(Dune::IOError, "Error parsing " << fileName << " "
343  "file pos " << pos
344  << ": Expected '" << format << "', only read " << c << " entries instead of " << cnt << ".");
345  }
346 
347  // skip over the rest of the line, including the terminating newline
348  void skipline(FILE * file)
349  {
350  int c;
351  do {
352  c = std::fgetc(file);
353  } while(c != '\n' && c != EOF);
354  }
355 
356  public:
357 
358  GmshReaderParser(Dune::GridFactory<GridType>& _factory, bool v, bool i) :
359  factory(_factory), verbose(v), insert_boundary_segments(i) {}
360 
361  std::vector<int> & boundaryIdMap()
362  {
363  return boundary_id_to_physical_entity;
364  }
365 
367  std::vector<int> & elementIndexMap()
368  {
369  return element_index_to_physical_entity;
370  }
371 
373  std::vector<std::string> & physicalEntityNames()
374  {
375  return physical_entity_names;
376  }
377 
378  void read (const std::string& f)
379  {
380  if (verbose) std::cout << "Reading " << dim << "d Gmsh grid..." << std::endl;
381 
382  // open file name, we use C I/O
383  fileName = f;
384  FILE* file = fopen(fileName.c_str(),"rb");
385  if (file==0)
386  DUNE_THROW(Dune::IOError, "Could not open " << fileName);
387 
388  //=========================================
389  // Header: Read vertices into vector
390  // Check vertices that are needed
391  //=========================================
392 
393  number_of_real_vertices = 0;
394  boundary_element_count = 0;
395  element_count = 0;
396 
397  // process header
398  double version_number;
399  int file_type, data_size;
400 
401  readfile(file,1,"%s\n",buf);
402  if (strcmp(buf,"$MeshFormat")!=0)
403  DUNE_THROW(Dune::IOError, "expected $MeshFormat in first line");
404  readfile(file,3,"%lg %d %d\n",&version_number,&file_type,&data_size);
405  if( (version_number < 2.0) || (version_number > 2.2) )
406  DUNE_THROW(Dune::IOError, "can only read Gmsh version 2 files");
407  if (verbose) std::cout << "version " << version_number << " Gmsh file detected" << std::endl;
408  readfile(file,1,"%s\n",buf);
409  if (strcmp(buf,"$EndMeshFormat")!=0)
410  DUNE_THROW(Dune::IOError, "expected $EndMeshFormat");
411 
412  // physical name section
413  physical_entity_names.clear();
414  readfile(file,1,"%s\n",buf);
415  if (strcmp(buf,"$PhysicalNames")==0) {
416  int number_of_names;
417  readfile(file,1,"%d\n",&number_of_names);
418  if (verbose) std::cout << "file contains " << number_of_names << " physical entities" << std::endl;
419  physical_entity_names.resize(number_of_names);
420  std::string buf_name;
421  for( int i = 0; i < number_of_names; ++i ) {
422  int dim, id;
423  readfile(file,3, "%d %d %s\n", &dim, &id, buf);
424  buf_name.assign(buf);
425  auto begin = buf_name.find_first_of('\"') + 1;
426  auto end = buf_name.find_last_of('\"') - begin;
427  physical_entity_names[id-1].assign(buf_name.substr(begin, end));
428  }
429  readfile(file,1,"%s\n",buf);
430  if (strcmp(buf,"$EndPhysicalNames")!=0)
431  DUNE_THROW(Dune::IOError, "expected $EndPhysicalNames");
432  readfile(file,1,"%s\n",buf);
433  }
434 
435  // node section
436  int number_of_nodes;
437 
438  if (strcmp(buf,"$Nodes")!=0)
439  DUNE_THROW(Dune::IOError, "expected $Nodes");
440  readfile(file,1,"%d\n",&number_of_nodes);
441  if (verbose) std::cout << "file contains " << number_of_nodes << " nodes" << std::endl;
442 
443  // read nodes
444  // The '+1' is due to the fact that gmsh numbers node starting from 1 rather than from 0
445  std::vector< GlobalVector > nodes( number_of_nodes+1 );
446  {
447  int id;
448  double x[ 3 ];
449  for( int i = 1; i <= number_of_nodes; ++i )
450  {
451  readfile(file,4, "%d %lg %lg %lg\n", &id, &x[ 0 ], &x[ 1 ], &x[ 2 ] );
452 
453  if (id > number_of_nodes) {
455  "Only dense sequences of node indices are currently supported (node index "
456  << id << " is invalid).");
457  }
458 
459  // just store node position
460  for( int j = 0; j < dimWorld; ++j )
461  nodes[ id ][ j ] = x[ j ];
462  }
463  readfile(file,1,"%s\n",buf);
464  if (strcmp(buf,"$EndNodes")!=0)
465  DUNE_THROW(Dune::IOError, "expected $EndNodes");
466  }
467 
468  // element section
469  readfile(file,1,"%s\n",buf);
470  if (strcmp(buf,"$Elements")!=0)
471  DUNE_THROW(Dune::IOError, "expected $Elements");
472  int number_of_elements;
473  readfile(file,1,"%d\n",&number_of_elements);
474  if (verbose) std::cout << "file contains " << number_of_elements << " elements" << std::endl;
475 
476  //=========================================
477  // Pass 1: Select and insert those vertices in the file that
478  // actually occur as corners of an element.
479  //=========================================
480 
481  off_t section_element_offset = ftello(file);
482  std::map<int,unsigned int> renumber;
483  for (int i=1; i<=number_of_elements; i++)
484  {
485  int id, elm_type, number_of_tags;
486  readfile(file,3,"%d %d %d ",&id,&elm_type,&number_of_tags);
487  for (int k=1; k<=number_of_tags; k++)
488  {
489  int blub;
490  readfile(file,1,"%d ",&blub);
491  // k == 1: physical entity
492  // k == 2: elementary entity (not used here)
493  // if version_number < 2.2:
494  // k == 3: mesh partition 0
495  // else
496  // k == 3: number of mesh partitions
497  // k => 4: mesh partition k-4
498  }
499  pass1HandleElement(file, elm_type, renumber, nodes);
500  }
501  if (verbose) std::cout << "number of real vertices = " << number_of_real_vertices << std::endl;
502  if (verbose) std::cout << "number of boundary elements = " << boundary_element_count << std::endl;
503  if (verbose) std::cout << "number of elements = " << element_count << std::endl;
504  readfile(file,1,"%s\n",buf);
505  if (strcmp(buf,"$EndElements")!=0)
506  DUNE_THROW(Dune::IOError, "expected $EndElements");
507  boundary_id_to_physical_entity.resize(boundary_element_count);
508  element_index_to_physical_entity.resize(element_count);
509 
510  //==============================================
511  // Pass 2: Insert boundary segments and elements
512  //==============================================
513 
514  fseeko(file, section_element_offset, SEEK_SET);
515  boundary_element_count = 0;
516  element_count = 0;
517  for (int i=1; i<=number_of_elements; i++)
518  {
519  int id, elm_type, number_of_tags;
520  readfile(file,3,"%d %d %d ",&id,&elm_type,&number_of_tags);
521  int physical_entity = -1;
522 
523  for (int k=1; k<=number_of_tags; k++)
524  {
525  int blub;
526  readfile(file,1,"%d ",&blub);
527  if (k==1) physical_entity = blub;
528  }
529  pass2HandleElement(file, elm_type, renumber, nodes, physical_entity);
530  }
531  readfile(file,1,"%s\n",buf);
532  if (strcmp(buf,"$EndElements")!=0)
533  DUNE_THROW(Dune::IOError, "expected $EndElements");
534 
535  fclose(file);
536  }
537 
543  void pass1HandleElement(FILE* file, const int elm_type,
544  std::map<int,unsigned int> & renumber,
545  const std::vector< GlobalVector > & nodes)
546  {
547  // some data about gmsh elements
548  const int nDofs[16] = {-1, 2, 3, 4, 4, 8, 6, 5, 3, 6, -1, 10, -1, -1, -1, 1};
549  const int nVertices[16] = {-1, 2, 3, 4, 4, 8, 6, 5, 2, 3, -1, 4, -1, -1, -1, 1};
550  const int elementDim[16] = {-1, 1, 2, 2, 3, 3, 3, 3, 1, 2, -1, 3, -1, -1, -1, 0};
551 
552  // test whether we support the element type
553  if ( not (elm_type > 0 && elm_type <= 15 // index in suitable range?
554  && (elementDim[elm_type] == dim || elementDim[elm_type] == (dim-1) ) ) ) // real element or boundary element?
555  {
556  skipline(file); // skip rest of line if element is unknown
557  return;
558  }
559 
560  // The format string for parsing is n times '%d' in a row
561  std::string formatString = "%d";
562  for (int i=1; i<nDofs[elm_type]; i++)
563  formatString += " %d";
564  formatString += "\n";
565 
566  // '10' is the largest number of dofs we may encounter in a .msh file
567  std::vector<int> elementDofs(10);
568 
569  readfile(file,nDofs[elm_type], formatString.c_str(),
570  &(elementDofs[0]),&(elementDofs[1]),&(elementDofs[2]),
571  &(elementDofs[3]),&(elementDofs[4]),&(elementDofs[5]),
572  &(elementDofs[6]),&(elementDofs[7]),&(elementDofs[8]),
573  &(elementDofs[9]));
574 
575  // insert each vertex if it hasn't been inserted already
576  for (int i=0; i<nVertices[elm_type]; i++)
577  if (renumber.find(elementDofs[i])==renumber.end())
578  {
579  renumber[elementDofs[i]] = number_of_real_vertices++;
580  factory.insertVertex(nodes[elementDofs[i]]);
581  }
582 
583  // count elements and boundary elements
584  if (elementDim[elm_type] == dim)
585  element_count++;
586  else
587  boundary_element_count++;
588 
589  }
590 
591 
592 
593  // generic-case: This is not supposed to be used at runtime.
594  template <class E, class V, class V2>
595  void boundarysegment_insert(
596  const V&,
597  const E&,
598  const V2&
599  )
600  {
601  DUNE_THROW(Dune::IOError, "tried to create a 3D boundary segment in a non-3D Grid");
602  }
603 
604  // 3d-case:
605  template <class E, class V>
606  void boundarysegment_insert(
607  const std::vector<FieldVector<double, 3> >& nodes,
608  const E& elementDofs,
609  const V& vertices
610  )
611  {
612  std::array<FieldVector<double,dimWorld>, 6> v;
613  for (int i=0; i<6; i++)
614  for (int j=0; j<dimWorld; j++)
615  v[i][j] = nodes[elementDofs[i]][j];
616 
617  BoundarySegment<dim,dimWorld>* newBoundarySegment
618  = (BoundarySegment<dim,dimWorld>*) new GmshReaderQuadraticBoundarySegment< 3, 3 >( v[0], v[1], v[2],
619  v[3], v[4], v[5] );
620 
621  factory.insertBoundarySegment( vertices,
622  std::shared_ptr<BoundarySegment<dim,dimWorld> >(newBoundarySegment) );
623  }
624 
625 
626 
631  virtual void pass2HandleElement(FILE* file, const int elm_type,
632  std::map<int,unsigned int> & renumber,
633  const std::vector< GlobalVector > & nodes,
634  const int physical_entity)
635  {
636  // some data about gmsh elements
637  const int nDofs[16] = {-1, 2, 3, 4, 4, 8, 6, 5, 3, 6, -1, 10, -1, -1, -1, 1};
638  const int nVertices[16] = {-1, 2, 3, 4, 4, 8, 6, 5, 2, 3, -1, 4, -1, -1, -1, 1};
639  const int elementDim[16] = {-1, 1, 2, 2, 3, 3, 3, 3, 1, 2, -1, 3, -1, -1, -1, 0};
640 
641  // test whether we support the element type
642  if ( not (elm_type > 0 && elm_type <= 15 // index in suitable range?
643  && (elementDim[elm_type] == dim || elementDim[elm_type] == (dim-1) ) ) ) // real element or boundary element?
644  {
645  skipline(file); // skip rest of line if element is unknown
646  return;
647  }
648 
649  // The format string for parsing is n times '%d' in a row
650  std::string formatString = "%d";
651  for (int i=1; i<nDofs[elm_type]; i++)
652  formatString += " %d";
653  formatString += "\n";
654 
655  // '10' is the largest number of dofs we may encounter in a .msh file
656  std::vector<int> elementDofs(10);
657 
658  readfile(file,nDofs[elm_type], formatString.c_str(),
659  &(elementDofs[0]),&(elementDofs[1]),&(elementDofs[2]),
660  &(elementDofs[3]),&(elementDofs[4]),&(elementDofs[5]),
661  &(elementDofs[6]),&(elementDofs[7]),&(elementDofs[8]),
662  &(elementDofs[9]));
663 
664  // correct differences between gmsh and Dune in the local vertex numbering
665  switch (elm_type)
666  {
667  case 3 : // 4-node quadrilateral
668  std::swap(elementDofs[2],elementDofs[3]);
669  break;
670  case 5 : // 8-node hexahedron
671  std::swap(elementDofs[2],elementDofs[3]);
672  std::swap(elementDofs[6],elementDofs[7]);
673  break;
674  case 7 : // 5-node pyramid
675  std::swap(elementDofs[2],elementDofs[3]);
676  break;
677  }
678 
679  // renumber corners to account for the explicitly given vertex
680  // numbering in the file
681  std::vector<unsigned int> vertices(nVertices[elm_type]);
682 
683  for (int i=0; i<nVertices[elm_type]; i++)
684  vertices[i] = renumber[elementDofs[i]];
685 
686  // If it is an element, insert it as such
687  if (elementDim[elm_type] == dim) {
688 
689  switch (elm_type)
690  {
691  case 1 : // 2-node line
692  factory.insertElement(Dune::GeometryTypes::line,vertices);
693  break;
694  case 2 : // 3-node triangle
696  break;
697  case 3 : // 4-node quadrilateral
699  break;
700  case 4 : // 4-node tetrahedron
702  break;
703  case 5 : // 8-node hexahedron
705  break;
706  case 6 : // 6-node prism
707  factory.insertElement(Dune::GeometryTypes::prism,vertices);
708  break;
709  case 7 : // 5-node pyramid
711  break;
712  case 9 : // 6-node triangle
714  break;
715  case 11 : // 10-node tetrahedron
717  break;
718  }
719 
720  } else {
721  // it must be a boundary segment then
722  if (insert_boundary_segments) {
723 
724  switch (elm_type)
725  {
726  case 1 : // 2-node line
727  factory.insertBoundarySegment(vertices);
728  break;
729 
730  case 2 : // 3-node triangle
731  factory.insertBoundarySegment(vertices);
732  break;
733 
734  case 3 : // 4-node quadrilateral
735  factory.insertBoundarySegment(vertices);
736  break;
737 
738  case 15 : // 1-node point
739  factory.insertBoundarySegment(vertices);
740  break;
741 
742  case 8 : { // 3-node line
743  std::array<FieldVector<double,dimWorld>, 3> v;
744  for (int i=0; i<dimWorld; i++) {
745  v[0][i] = nodes[elementDofs[0]][i];
746  v[1][i] = nodes[elementDofs[2]][i]; // yes, the renumbering is intended!
747  v[2][i] = nodes[elementDofs[1]][i];
748  }
749  BoundarySegment<dim,dimWorld>* newBoundarySegment
750  = (BoundarySegment<dim,dimWorld>*) new GmshReaderQuadraticBoundarySegment< 2, dimWorld >(v[0], v[1], v[2]);
751  factory.insertBoundarySegment(vertices,
752  std::shared_ptr<BoundarySegment<dim,dimWorld> >(newBoundarySegment));
753  break;
754  }
755  case 9 : { // 6-node triangle
756  boundarysegment_insert(nodes, elementDofs, vertices);
757  break;
758  }
759  default: {
760  DUNE_THROW(Dune::IOError, "GmshReader does not support using element-type " << elm_type << " for boundary segments");
761  break;
762  }
763 
764  }
765 
766  }
767  }
768 
769  // count elements and boundary elements
770  if (elementDim[elm_type] == dim) {
771  element_index_to_physical_entity[element_count] = physical_entity;
772  element_count++;
773  } else {
774  boundary_id_to_physical_entity[boundary_element_count] = physical_entity;
775  boundary_element_count++;
776  }
777 
778  }
779 
780  };
781 
782  namespace Gmsh {
788  enum class ReaderOptions
789  {
790  verbose = 1,
791  insertBoundarySegments = 2,
792  readElementData = 4,
793  readBoundaryData = 8
794  };
795 
797  constexpr ReaderOptions operator | (ReaderOptions a, ReaderOptions b)
798  {
799  return static_cast<ReaderOptions>(
800  static_cast<int>(a) | static_cast<int>(b)
801  );
802  }
803 
805  constexpr bool operator & (ReaderOptions a, ReaderOptions b)
806  {
807  return static_cast<int>(a) & static_cast<int>(b);
808  }
809 
810  } // end namespace Gmsh
811 
836  template<typename GridType>
838  {
840 
859  static void doRead(Dune::GridFactory<GridType> &factory,
860  const std::string &fileName,
861  std::vector<int>& boundarySegmentToPhysicalEntity,
862  std::vector<int>& elementToPhysicalEntity,
863  bool verbose, bool insertBoundarySegments)
864  {
865  // register boundary segment to boundary segment factory for possible load balancing
866  // this needs to be done on all cores since the type might not be known otherwise
867  GmshReaderQuadraticBoundarySegment< Grid::dimension, Grid::dimensionworld >::registerFactory();
868 
869 #ifndef NDEBUG
870  // check that this method is called on all cores
871  factory.comm().barrier();
872 #endif
873 
874  // create parse object and read grid on process 0
875  if (factory.comm().rank() == 0)
876  {
877  GmshReaderParser<Grid> parser(factory,verbose,insertBoundarySegments);
878  parser.read(fileName);
879 
880  boundarySegmentToPhysicalEntity = std::move(parser.boundaryIdMap());
881  elementToPhysicalEntity = std::move(parser.elementIndexMap());
882  }
883  else
884  {
885  boundarySegmentToPhysicalEntity = {};
886  elementToPhysicalEntity = {};
887  }
888  }
889 
891 
910  template<class T>
911  static T &discarded(T &&value) { return static_cast<T&>(value); }
912 
913  struct DataArg {
914  std::vector<int> *data_ = nullptr;
915  DataArg(std::vector<int> &data) : data_(&data) {}
916  DataArg(const decltype(std::ignore)&) {}
917  DataArg() = default;
918  };
919 
920  struct DataFlagArg : DataArg {
921  bool flag_ = false;
922  using DataArg::DataArg;
923  DataFlagArg(bool flag) : flag_(flag) {}
924  };
925 
926  public:
927  typedef GridType Grid;
928 
935  static std::unique_ptr<Grid> read (const std::string& fileName, bool verbose = true, bool insertBoundarySegments=true)
936  {
937  // make a grid factory
938  Dune::GridFactory<Grid> factory;
939 
940  read(factory, fileName, verbose, insertBoundarySegments);
941 
942  return factory.createGrid();
943  }
944 
964  static std::unique_ptr<Grid> read (const std::string& fileName,
965  std::vector<int>& boundarySegmentToPhysicalEntity,
966  std::vector<int>& elementToPhysicalEntity,
967  bool verbose = true, bool insertBoundarySegments=true)
968  {
969  // make a grid factory
970  Dune::GridFactory<Grid> factory;
971 
972  doRead(
973  factory, fileName, boundarySegmentToPhysicalEntity,
974  elementToPhysicalEntity, verbose, insertBoundarySegments
975  );
976 
977  return factory.createGrid();
978  }
979 
981  static void read (Dune::GridFactory<Grid>& factory, const std::string& fileName,
982  bool verbose = true, bool insertBoundarySegments=true)
983  {
984  doRead(
985  factory, fileName, discarded(std::vector<int>{}),
986  discarded(std::vector<int>{}), verbose, insertBoundarySegments
987  );
988  }
989 
991 
1014  static void read (Dune::GridFactory<Grid> &factory,
1015  const std::string &fileName,
1016  DataFlagArg boundarySegmentData,
1017  DataArg elementData,
1018  bool verbose=true)
1019  {
1020  doRead(
1021  factory, fileName,
1022  boundarySegmentData.data_
1023  ? *boundarySegmentData.data_ : discarded(std::vector<int>{}),
1024  elementData.data_
1025  ? *elementData.data_ : discarded(std::vector<int>{}),
1026  verbose,
1027  boundarySegmentData.flag_ || boundarySegmentData.data_
1028  );
1029  }
1030 
1051  static void read (Dune::GridFactory<Grid>& factory,
1052  const std::string& fileName,
1053  std::vector<int>& boundarySegmentToPhysicalEntity,
1054  std::vector<int>& elementToPhysicalEntity,
1055  bool verbose, bool insertBoundarySegments)
1056  {
1057  doRead(
1058  factory, fileName, boundarySegmentToPhysicalEntity,
1059  elementToPhysicalEntity, verbose, insertBoundarySegments
1060  );
1061  }
1062 
1064  //\{
1065 
1067 
1068  static constexpr Opts defaultOpts =
1069  Opts::verbose | Opts::insertBoundarySegments | Opts::readElementData | Opts::readBoundaryData;
1070 
1072 
1095  GmshReader(const std::string& fileName,
1096  Gmsh::ReaderOptions options = defaultOpts)
1097  {
1098  gridFactory_ = std::make_unique<Dune::GridFactory<Grid>>();
1099  readGridFile(fileName, *gridFactory_, options);
1100  }
1101 
1109  GmshReader(const std::string& fileName, GridFactory<Grid>& factory,
1110  Gmsh::ReaderOptions options = defaultOpts)
1111  {
1112  readGridFile(fileName, factory, options);
1113  }
1114 
1116  const std::vector<int>& elementData () const
1117  {
1118  checkElementData();
1119  return elementIndexToGmshPhysicalEntity_;
1120  }
1121 
1123  const std::vector<int>& boundaryData () const
1124  {
1125  checkBoundaryData();
1126  return boundarySegmentIndexToGmshPhysicalEntity_;
1127  }
1128 
1133  bool hasElementData () const
1134  { return hasElementData_ && !extractedElementData_; }
1135 
1140  bool hasBoundaryData () const
1141  { return hasBoundaryData_ && !extractedBoundaryData_; }
1142 
1144  std::vector<int> extractElementData ()
1145  {
1146  checkElementData();
1147  extractedElementData_ = true;
1148  return std::move(elementIndexToGmshPhysicalEntity_);
1149  }
1150 
1152  std::vector<int> extractBoundaryData ()
1153  {
1154  checkBoundaryData();
1155  extractedBoundaryData_ = true;
1156  return std::move(boundarySegmentIndexToGmshPhysicalEntity_);
1157  }
1158 
1160  std::unique_ptr<Grid> createGrid ()
1161  {
1162  if (!gridFactory_)
1164  "This GmshReader has been constructed with a Dune::GridFactory. "
1165  << "This grid factory has been filled with all information to create a grid. "
1166  << "Please use this factory to create the grid by calling factory.createGrid(). "
1167  << "Alternatively use the constructor without passing the factory in combination with this member function."
1168  );
1169 
1170  return gridFactory_->createGrid();
1171  }
1172 
1173  //\}
1174 
1175  private:
1176  void checkElementData () const
1177  {
1178  if (!hasElementData_)
1180  "This GmshReader has been constructed without the option 'readElementData'. "
1181  << "Please enable reading element data by passing the option 'Gmsh::ReaderOpts::readElementData' "
1182  << "to the constructor of this class."
1183  );
1184 
1185  if (extractedElementData_)
1187  "The element data has already been extracted from this GmshReader "
1188  << "via a function call to reader.extractElementData(). Use the extracted data or "
1189  << "read the grid data from file again by constructing a new reader."
1190  );
1191  }
1192 
1193  void checkBoundaryData () const
1194  {
1195  if (!hasBoundaryData_)
1197  "This GmshReader has been constructed without the option 'readBoundaryData'. "
1198  << "Please enable reading boundary data by passing the option 'Gmsh::ReaderOpts::readBoundaryData' "
1199  << "to the constructor of this class."
1200  );
1201 
1202  if (extractedBoundaryData_)
1204  "The boundary data has already been extracted from this GmshReader "
1205  << "via a function call to reader.extractBoundaryData(). Use the extracted data or "
1206  << "read the grid data from file again by constructing a new reader."
1207  );
1208  }
1209 
1210  void readGridFile (const std::string& fileName, GridFactory<Grid>& factory, Gmsh::ReaderOptions options)
1211  {
1212  const bool verbose = options & Opts::verbose;
1213  const bool insertBoundarySegments = options & Opts::insertBoundarySegments;
1214  const bool readBoundaryData = options & Opts::readBoundaryData;
1215  const bool readElementData = options & Opts::readElementData;
1216 
1217  doRead(
1218  factory, fileName, boundarySegmentIndexToGmshPhysicalEntity_,
1219  elementIndexToGmshPhysicalEntity_, verbose,
1220  readBoundaryData || insertBoundarySegments
1221  );
1222 
1223  // clear unwanted data
1224  if (!readBoundaryData)
1225  boundarySegmentIndexToGmshPhysicalEntity_ = std::vector<int>{};
1226  if (!readElementData)
1227  elementIndexToGmshPhysicalEntity_ = std::vector<int>{};
1228 
1229  hasElementData_ = readElementData;
1230  hasBoundaryData_ = readBoundaryData;
1231  }
1232 
1233  std::unique_ptr<Dune::GridFactory<Grid>> gridFactory_;
1234 
1235  std::vector<int> elementIndexToGmshPhysicalEntity_;
1236  std::vector<int> boundarySegmentIndexToGmshPhysicalEntity_;
1237 
1238  bool hasElementData_;
1239  bool hasBoundaryData_;
1240 
1241  // for better error messages, we keep track of these separately
1242  bool extractedElementData_ = false;
1243  bool extractedBoundaryData_ = false;
1244  };
1245 
1248 } // namespace Dune
1249 
1250 #endif
Base class for grid boundary segments of arbitrary geometry.
int rank() const
Return rank, is between 0 and size()-1.
Definition: communication.hh:114
int barrier() const
Wait until all processes have arrived at this point in the program.
Definition: communication.hh:267
derived_type & axpy(const field_type &a, const DenseVector< Other > &x)
vector space axpy operation ( *this += a x )
Definition: densevector.hh:575
FieldTraits< value_type >::real_type two_norm() const
two norm sqrt(sum over squared values of entries)
Definition: densevector.hh:641
vector space out of a tensor product of fields.
Definition: fvector.hh:95
dimension independent parts for GmshReaderParser
Definition: gmshreader.hh:302
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:543
std::vector< int > & elementIndexMap()
Returns a map for the gmsh physical entity id (1-based index) for each entity of codim 0.
Definition: gmshreader.hh:367
std::vector< std::string > & physicalEntityNames()
Returns the names of the gmsh physical entities (0-based index)
Definition: gmshreader.hh:373
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:631
Read Gmsh mesh file.
Definition: gmshreader.hh:838
std::vector< int > extractBoundaryData()
Erase boundary data from reader and return the data.
Definition: gmshreader.hh:1152
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:1014
const std::vector< int > & elementData() const
Access element data (maps element index to Gmsh physical entity)
Definition: gmshreader.hh:1116
static void read(Dune::GridFactory< Grid > &factory, const std::string &fileName, bool verbose=true, bool insertBoundarySegments=true)
Definition: gmshreader.hh:981
std::unique_ptr< Grid > createGrid()
Create the grid.
Definition: gmshreader.hh:1160
std::vector< int > extractElementData()
Erase element data from reader and return the data.
Definition: gmshreader.hh:1144
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:1051
bool hasElementData() const
If element data is available.
Definition: gmshreader.hh:1133
bool hasBoundaryData() const
If boundary data is available.
Definition: gmshreader.hh:1140
GmshReader(const std::string &fileName, GridFactory< Grid > &factory, Gmsh::ReaderOptions options=defaultOpts)
Construct a Gmsh reader object from a file name and a grid factory.
Definition: gmshreader.hh:1109
GmshReader(const std::string &fileName, Gmsh::ReaderOptions options=defaultOpts)
Construct a Gmsh reader object (alternatively use one of the static member functions)
Definition: gmshreader.hh:1095
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:964
static std::unique_ptr< Grid > read(const std::string &fileName, bool verbose=true, bool insertBoundarySegments=true)
Definition: gmshreader.hh:935
const std::vector< int > & boundaryData() const
Access boundary data (maps boundary segment index to Gmsh physical entity)
Definition: gmshreader.hh:1123
Communication comm() const
Return the Communication used by the grid factory.
Definition: gridfactory.hh:258
Provide a generic factory class for unstructured grids.
Definition: gridfactory.hh:275
virtual void insertVertex([[maybe_unused]] const FieldVector< ctype, dimworld > &pos)
Insert a vertex into the coarse grid.
Definition: gridfactory.hh:296
virtual void insertBoundarySegment([[maybe_unused]] const std::vector< unsigned int > &vertices)
insert a boundary segment
Definition: gridfactory.hh:325
virtual void insertElement([[maybe_unused]] const GeometryType &type, [[maybe_unused]] const std::vector< unsigned int > &vertices)
Insert an element into the coarse grid.
Definition: gridfactory.hh:307
virtual std::unique_ptr< GridType > createGrid()
Finalize grid creation and hand over the grid.
Definition: gridfactory.hh:333
Default exception class for I/O errors.
Definition: exceptions.hh:231
Default exception if a function was called while the object is not in a valid state for that function...
Definition: exceptions.hh:281
Provide a generic factory class for unstructured grids.
A few common exception classes.
int renumber(const Dune::GeometryType &t, int i)
renumber VTK <-> Dune
Definition: common.hh:186
Implements a vector constructed from a given type representing a field and a compile-time given size.
#define DUNE_THROW(E, m)
Definition: exceptions.hh:218
constexpr GeometryType line
GeometryType representing a line.
Definition: type.hh:498
constexpr GeometryType prism
GeometryType representing a 3D prism.
Definition: type.hh:528
constexpr GeometryType triangle
GeometryType representing a triangle.
Definition: type.hh:504
constexpr GeometryType quadrilateral
GeometryType representing a quadrilateral (a square).
Definition: type.hh:510
constexpr GeometryType hexahedron
GeometryType representing a hexahedron.
Definition: type.hh:534
constexpr GeometryType pyramid
GeometryType representing a 3D pyramid.
Definition: type.hh:522
constexpr GeometryType tetrahedron
GeometryType representing a tetrahedron.
Definition: type.hh:516
ReaderOptions
Option for the Gmsh mesh file reader.
Definition: gmshreader.hh:789
static std::string formatString(const std::string &s, const T &... args)
Format values according to printf format string.
Definition: stringutility.hh:73
Dune namespace.
Definition: alignedallocator.hh:13
Base class for classes implementing geometries of boundary segments.
Definition: boundarysegment.hh:94
Options for read operation.
Definition: gmshreader.hh:39
GeometryOrder
Definition: gmshreader.hh:40
@ firstOrder
edges are straight lines.
Definition: gmshreader.hh:42
@ secondOrder
quadratic boundary approximation.
Definition: gmshreader.hh:44
A unique label for each type of element that can occur in a grid.
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.80.0 (Mar 28, 23:30, 2024)