dune-grid 2.9.0
gmshreader.hh
Go to the documentation of this file.
1// SPDX-FileCopyrightText: Copyright (C) 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
21#include <dune/common/exceptions.hh>
22#include <dune/common/fvector.hh>
23
24#include <dune/geometry/type.hh>
25
28
29namespace 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_,
179 Dune::FieldVector<double,3> p2_, Dune::FieldVector<double,3> p3_,
180 Dune::FieldVector<double,3> p4_, Dune::FieldVector<double,3> p5_)
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 {
208 Dune::FieldVector<double,3> y;
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);
238 Dune::FieldVector<double,3> d1,d2;
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
311 // read buffer
312 char buf[512];
313 std::string fileName;
314 // exported data
317
318 // static data
319 static const int dim = GridType::dimension;
320 static const int dimWorld = GridType::dimensionworld;
321 static_assert( (dimWorld <= 3), "GmshReader requires dimWorld <= 3." );
322
323 // typedefs
324 typedef FieldVector< double, dimWorld > GlobalVector;
325
326 // don't use something like
327 // readfile(file, 1, "%s\n", buf);
328 // to skip the rest of of the line -- that will only skip the next
329 // whitespace-separated word! Use skipline() instead.
330 void readfile(FILE * file, int cnt, const char * format,
331 void* t1, void* t2 = 0, void* t3 = 0, void* t4 = 0,
332 void* t5 = 0, void* t6 = 0, void* t7 = 0, void* t8 = 0,
333 void* t9 = 0, void* t10 = 0)
334 {
335 off_t pos = ftello(file);
336 int c = fscanf(file, format, t1, t2, t3, t4, t5, t6, t7, t8, t9, t10);
337 if (c != cnt)
338 DUNE_THROW(Dune::IOError, "Error parsing " << fileName << " "
339 "file pos " << pos
340 << ": Expected '" << format << "', only read " << c << " entries instead of " << cnt << ".");
341 }
342
343 // skip over the rest of the line, including the terminating newline
344 void skipline(FILE * file)
345 {
346 int c;
347 do {
348 c = std::fgetc(file);
349 } while(c != '\n' && c != EOF);
350 }
351
352 public:
353
355 factory(_factory), verbose(v), insert_boundary_segments(i) {}
356
357 std::vector<int> & boundaryIdMap()
358 {
360 }
361
362 std::vector<int> & elementIndexMap()
363 {
365 }
366
367 void read (const std::string& f)
368 {
369 if (verbose) std::cout << "Reading " << dim << "d Gmsh grid..." << std::endl;
370
371 // open file name, we use C I/O
372 fileName = f;
373 FILE* file = fopen(fileName.c_str(),"rb");
374 if (file==0)
375 DUNE_THROW(Dune::IOError, "Could not open " << fileName);
376
377 //=========================================
378 // Header: Read vertices into vector
379 // Check vertices that are needed
380 //=========================================
381
384 element_count = 0;
385
386 // process header
387 double version_number;
388 int file_type, data_size;
389
390 readfile(file,1,"%s\n",buf);
391 if (strcmp(buf,"$MeshFormat")!=0)
392 DUNE_THROW(Dune::IOError, "expected $MeshFormat in first line");
393 readfile(file,3,"%lg %d %d\n",&version_number,&file_type,&data_size);
394 if( (version_number < 2.0) || (version_number > 2.2) )
395 DUNE_THROW(Dune::IOError, "can only read Gmsh version 2 files");
396 if (verbose) std::cout << "version " << version_number << " Gmsh file detected" << std::endl;
397 readfile(file,1,"%s\n",buf);
398 if (strcmp(buf,"$EndMeshFormat")!=0)
399 DUNE_THROW(Dune::IOError, "expected $EndMeshFormat");
400
401 // node section
402 int number_of_nodes;
403
404 readfile(file,1,"%s\n",buf);
405 if (strcmp(buf,"$Nodes")!=0)
406 DUNE_THROW(Dune::IOError, "expected $Nodes");
407 readfile(file,1,"%d\n",&number_of_nodes);
408 if (verbose) std::cout << "file contains " << number_of_nodes << " nodes" << std::endl;
409
410 // read nodes
411 // The '+1' is due to the fact that gmsh numbers node starting from 1 rather than from 0
412 std::vector< GlobalVector > nodes( number_of_nodes+1 );
413 {
414 int id;
415 double x[ 3 ];
416 for( int i = 1; i <= number_of_nodes; ++i )
417 {
418 readfile(file,4, "%d %lg %lg %lg\n", &id, &x[ 0 ], &x[ 1 ], &x[ 2 ] );
419
420 if (id > number_of_nodes) {
421 DUNE_THROW(Dune::IOError,
422 "Only dense sequences of node indices are currently supported (node index "
423 << id << " is invalid).");
424 }
425
426 // just store node position
427 for( int j = 0; j < dimWorld; ++j )
428 nodes[ id ][ j ] = x[ j ];
429 }
430 readfile(file,1,"%s\n",buf);
431 if (strcmp(buf,"$EndNodes")!=0)
432 DUNE_THROW(Dune::IOError, "expected $EndNodes");
433 }
434
435 // element section
436 readfile(file,1,"%s\n",buf);
437 if (strcmp(buf,"$Elements")!=0)
438 DUNE_THROW(Dune::IOError, "expected $Elements");
439 int number_of_elements;
440 readfile(file,1,"%d\n",&number_of_elements);
441 if (verbose) std::cout << "file contains " << number_of_elements << " elements" << std::endl;
442
443 //=========================================
444 // Pass 1: Select and insert those vertices in the file that
445 // actually occur as corners of an element.
446 //=========================================
447
448 off_t section_element_offset = ftello(file);
449 std::map<int,unsigned int> renumber;
450 for (int i=1; i<=number_of_elements; i++)
451 {
452 int id, elm_type, number_of_tags;
453 readfile(file,3,"%d %d %d ",&id,&elm_type,&number_of_tags);
454 for (int k=1; k<=number_of_tags; k++)
455 {
456 int blub;
457 readfile(file,1,"%d ",&blub);
458 // k == 1: physical entity (not used here)
459 // k == 2: elementary entity (not used here either)
460 // if version_number < 2.2:
461 // k == 3: mesh partition 0
462 // else
463 // k == 3: number of mesh partitions
464 // k => 4: mesh partition k-4
465 }
466 pass1HandleElement(file, elm_type, renumber, nodes);
467 }
468 if (verbose) std::cout << "number of real vertices = " << number_of_real_vertices << std::endl;
469 if (verbose) std::cout << "number of boundary elements = " << boundary_element_count << std::endl;
470 if (verbose) std::cout << "number of elements = " << element_count << std::endl;
471 readfile(file,1,"%s\n",buf);
472 if (strcmp(buf,"$EndElements")!=0)
473 DUNE_THROW(Dune::IOError, "expected $EndElements");
476
477 //==============================================
478 // Pass 2: Insert boundary segments and elements
479 //==============================================
480
481 fseeko(file, section_element_offset, SEEK_SET);
483 element_count = 0;
484 for (int i=1; i<=number_of_elements; i++)
485 {
486 int id, elm_type, number_of_tags;
487 readfile(file,3,"%d %d %d ",&id,&elm_type,&number_of_tags);
488 int physical_entity = -1;
489
490 for (int k=1; k<=number_of_tags; k++)
491 {
492 int blub;
493 readfile(file,1,"%d ",&blub);
494 if (k==1) physical_entity = blub;
495 }
496 pass2HandleElement(file, elm_type, renumber, nodes, physical_entity);
497 }
498 readfile(file,1,"%s\n",buf);
499 if (strcmp(buf,"$EndElements")!=0)
500 DUNE_THROW(Dune::IOError, "expected $EndElements");
501
502 fclose(file);
503 }
504
510 void pass1HandleElement(FILE* file, const int elm_type,
511 std::map<int,unsigned int> & renumber,
512 const std::vector< GlobalVector > & nodes)
513 {
514 // some data about gmsh elements
515 const int nDofs[16] = {-1, 2, 3, 4, 4, 8, 6, 5, 3, 6, -1, 10, -1, -1, -1, 1};
516 const int nVertices[16] = {-1, 2, 3, 4, 4, 8, 6, 5, 2, 3, -1, 4, -1, -1, -1, 1};
517 const int elementDim[16] = {-1, 1, 2, 2, 3, 3, 3, 3, 1, 2, -1, 3, -1, -1, -1, 0};
518
519 // test whether we support the element type
520 if ( not (elm_type > 0 && elm_type <= 15 // index in suitable range?
521 && (elementDim[elm_type] == dim || elementDim[elm_type] == (dim-1) ) ) ) // real element or boundary element?
522 {
523 skipline(file); // skip rest of line if element is unknown
524 return;
525 }
526
527 // The format string for parsing is n times '%d' in a row
528 std::string formatString = "%d";
529 for (int i=1; i<nDofs[elm_type]; i++)
530 formatString += " %d";
531 formatString += "\n";
532
533 // '10' is the largest number of dofs we may encounter in a .msh file
534 std::vector<int> elementDofs(10);
535
536 readfile(file,nDofs[elm_type], formatString.c_str(),
537 &(elementDofs[0]),&(elementDofs[1]),&(elementDofs[2]),
538 &(elementDofs[3]),&(elementDofs[4]),&(elementDofs[5]),
539 &(elementDofs[6]),&(elementDofs[7]),&(elementDofs[8]),
540 &(elementDofs[9]));
541
542 // insert each vertex if it hasn't been inserted already
543 for (int i=0; i<nVertices[elm_type]; i++)
544 if (renumber.find(elementDofs[i])==renumber.end())
545 {
546 renumber[elementDofs[i]] = number_of_real_vertices++;
547 factory.insertVertex(nodes[elementDofs[i]]);
548 }
549
550 // count elements and boundary elements
551 if (elementDim[elm_type] == dim)
553 else
555
556 }
557
558
559
560 // generic-case: This is not supposed to be used at runtime.
561 template <class E, class V, class V2>
563 const V&,
564 const E&,
565 const V2&
566 )
567 {
568 DUNE_THROW(Dune::IOError, "tried to create a 3D boundary segment in a non-3D Grid");
569 }
570
571 // 3d-case:
572 template <class E, class V>
574 const std::vector<FieldVector<double, 3> >& nodes,
575 const E& elementDofs,
576 const V& vertices
577 )
578 {
579 std::array<FieldVector<double,dimWorld>, 6> v;
580 for (int i=0; i<6; i++)
581 for (int j=0; j<dimWorld; j++)
582 v[i][j] = nodes[elementDofs[i]][j];
583
584 BoundarySegment<dim,dimWorld>* newBoundarySegment
585 = (BoundarySegment<dim,dimWorld>*) new GmshReaderQuadraticBoundarySegment< 3, 3 >( v[0], v[1], v[2],
586 v[3], v[4], v[5] );
587
588 factory.insertBoundarySegment( vertices,
589 std::shared_ptr<BoundarySegment<dim,dimWorld> >(newBoundarySegment) );
590 }
591
592
593
598 virtual void pass2HandleElement(FILE* file, const int elm_type,
599 std::map<int,unsigned int> & renumber,
600 const std::vector< GlobalVector > & nodes,
601 const int physical_entity)
602 {
603 // some data about gmsh elements
604 const int nDofs[16] = {-1, 2, 3, 4, 4, 8, 6, 5, 3, 6, -1, 10, -1, -1, -1, 1};
605 const int nVertices[16] = {-1, 2, 3, 4, 4, 8, 6, 5, 2, 3, -1, 4, -1, -1, -1, 1};
606 const int elementDim[16] = {-1, 1, 2, 2, 3, 3, 3, 3, 1, 2, -1, 3, -1, -1, -1, 0};
607
608 // test whether we support the element type
609 if ( not (elm_type > 0 && elm_type <= 15 // index in suitable range?
610 && (elementDim[elm_type] == dim || elementDim[elm_type] == (dim-1) ) ) ) // real element or boundary element?
611 {
612 skipline(file); // skip rest of line if element is unknown
613 return;
614 }
615
616 // The format string for parsing is n times '%d' in a row
617 std::string formatString = "%d";
618 for (int i=1; i<nDofs[elm_type]; i++)
619 formatString += " %d";
620 formatString += "\n";
621
622 // '10' is the largest number of dofs we may encounter in a .msh file
623 std::vector<int> elementDofs(10);
624
625 readfile(file,nDofs[elm_type], formatString.c_str(),
626 &(elementDofs[0]),&(elementDofs[1]),&(elementDofs[2]),
627 &(elementDofs[3]),&(elementDofs[4]),&(elementDofs[5]),
628 &(elementDofs[6]),&(elementDofs[7]),&(elementDofs[8]),
629 &(elementDofs[9]));
630
631 // correct differences between gmsh and Dune in the local vertex numbering
632 switch (elm_type)
633 {
634 case 3 : // 4-node quadrilateral
635 std::swap(elementDofs[2],elementDofs[3]);
636 break;
637 case 5 : // 8-node hexahedron
638 std::swap(elementDofs[2],elementDofs[3]);
639 std::swap(elementDofs[6],elementDofs[7]);
640 break;
641 case 7 : // 5-node pyramid
642 std::swap(elementDofs[2],elementDofs[3]);
643 break;
644 }
645
646 // renumber corners to account for the explicitly given vertex
647 // numbering in the file
648 std::vector<unsigned int> vertices(nVertices[elm_type]);
649
650 for (int i=0; i<nVertices[elm_type]; i++)
651 vertices[i] = renumber[elementDofs[i]];
652
653 // If it is an element, insert it as such
654 if (elementDim[elm_type] == dim) {
655
656 switch (elm_type)
657 {
658 case 1 : // 2-node line
659 factory.insertElement(Dune::GeometryTypes::line,vertices);
660 break;
661 case 2 : // 3-node triangle
662 factory.insertElement(Dune::GeometryTypes::triangle,vertices);
663 break;
664 case 3 : // 4-node quadrilateral
665 factory.insertElement(Dune::GeometryTypes::quadrilateral,vertices);
666 break;
667 case 4 : // 4-node tetrahedron
668 factory.insertElement(Dune::GeometryTypes::tetrahedron,vertices);
669 break;
670 case 5 : // 8-node hexahedron
671 factory.insertElement(Dune::GeometryTypes::hexahedron,vertices);
672 break;
673 case 6 : // 6-node prism
674 factory.insertElement(Dune::GeometryTypes::prism,vertices);
675 break;
676 case 7 : // 5-node pyramid
677 factory.insertElement(Dune::GeometryTypes::pyramid,vertices);
678 break;
679 case 9 : // 6-node triangle
680 factory.insertElement(Dune::GeometryTypes::triangle,vertices);
681 break;
682 case 11 : // 10-node tetrahedron
683 factory.insertElement(Dune::GeometryTypes::tetrahedron,vertices);
684 break;
685 }
686
687 } else {
688 // it must be a boundary segment then
690
691 switch (elm_type)
692 {
693 case 1 : // 2-node line
694 factory.insertBoundarySegment(vertices);
695 break;
696
697 case 2 : // 3-node triangle
698 factory.insertBoundarySegment(vertices);
699 break;
700
701 case 3 : // 4-node quadrilateral
702 factory.insertBoundarySegment(vertices);
703 break;
704
705 case 15 : // 1-node point
706 factory.insertBoundarySegment(vertices);
707 break;
708
709 case 8 : { // 3-node line
710 std::array<FieldVector<double,dimWorld>, 3> v;
711 for (int i=0; i<dimWorld; i++) {
712 v[0][i] = nodes[elementDofs[0]][i];
713 v[1][i] = nodes[elementDofs[2]][i]; // yes, the renumbering is intended!
714 v[2][i] = nodes[elementDofs[1]][i];
715 }
716 BoundarySegment<dim,dimWorld>* newBoundarySegment
717 = (BoundarySegment<dim,dimWorld>*) new GmshReaderQuadraticBoundarySegment< 2, dimWorld >(v[0], v[1], v[2]);
718 factory.insertBoundarySegment(vertices,
719 std::shared_ptr<BoundarySegment<dim,dimWorld> >(newBoundarySegment));
720 break;
721 }
722 case 9 : { // 6-node triangle
723 boundarysegment_insert(nodes, elementDofs, vertices);
724 break;
725 }
726 default: {
727 DUNE_THROW(Dune::IOError, "GmshReader does not support using element-type " << elm_type << " for boundary segments");
728 break;
729 }
730
731 }
732
733 }
734 }
735
736 // count elements and boundary elements
737 if (elementDim[elm_type] == dim) {
740 } else {
743 }
744
745 }
746
747 };
748
749 namespace Gmsh {
755 enum class ReaderOptions
756 {
757 verbose = 1,
759 readElementData = 4,
761 };
762
765 {
766 return static_cast<ReaderOptions>(
767 static_cast<int>(a) | static_cast<int>(b)
768 );
769 }
770
773 {
774 return static_cast<int>(a) & static_cast<int>(b);
775 }
776
777 } // end namespace Gmsh
778
803 template<typename GridType>
805 {
807
826 static void doRead(Dune::GridFactory<GridType> &factory,
827 const std::string &fileName,
828 std::vector<int>& boundarySegmentToPhysicalEntity,
829 std::vector<int>& elementToPhysicalEntity,
830 bool verbose, bool insertBoundarySegments)
831 {
832 // register boundary segment to boundary segment factory for possible load balancing
833 // this needs to be done on all cores since the type might not be known otherwise
834 GmshReaderQuadraticBoundarySegment< Grid::dimension, Grid::dimensionworld >::registerFactory();
835
836#ifndef NDEBUG
837 // check that this method is called on all cores
838 factory.comm().barrier();
839#endif
840
841 // create parse object and read grid on process 0
842 if (factory.comm().rank() == 0)
843 {
844 GmshReaderParser<Grid> parser(factory,verbose,insertBoundarySegments);
845 parser.read(fileName);
846
847 boundarySegmentToPhysicalEntity = std::move(parser.boundaryIdMap());
848 elementToPhysicalEntity = std::move(parser.elementIndexMap());
849 }
850 else
851 {
852 boundarySegmentToPhysicalEntity = {};
853 elementToPhysicalEntity = {};
854 }
855 }
856
858
877 template<class T>
878 static T &discarded(T &&value) { return value; }
879
880 struct DataArg {
881 std::vector<int> *data_ = nullptr;
882 DataArg(std::vector<int> &data) : data_(&data) {}
883 DataArg(const decltype(std::ignore)&) {}
884 DataArg() = default;
885 };
886
887 struct DataFlagArg : DataArg {
888 bool flag_ = false;
889 using DataArg::DataArg;
890 DataFlagArg(bool flag) : flag_(flag) {}
891 };
892
893 public:
894 typedef GridType Grid;
895
902 static std::unique_ptr<Grid> read (const std::string& fileName, bool verbose = true, bool insertBoundarySegments=true)
903 {
904 // make a grid factory
906
907 read(factory, fileName, verbose, insertBoundarySegments);
908
909 return factory.createGrid();
910 }
911
931 static std::unique_ptr<Grid> read (const std::string& fileName,
932 std::vector<int>& boundarySegmentToPhysicalEntity,
933 std::vector<int>& elementToPhysicalEntity,
934 bool verbose = true, bool insertBoundarySegments=true)
935 {
936 // make a grid factory
938
939 doRead(
940 factory, fileName, boundarySegmentToPhysicalEntity,
941 elementToPhysicalEntity, verbose, insertBoundarySegments
942 );
943
944 return factory.createGrid();
945 }
946
948 static void read (Dune::GridFactory<Grid>& factory, const std::string& fileName,
949 bool verbose = true, bool insertBoundarySegments=true)
950 {
951 doRead(
952 factory, fileName, discarded(std::vector<int>{}),
953 discarded(std::vector<int>{}), verbose, insertBoundarySegments
954 );
955 }
956
958
981 static void read (Dune::GridFactory<Grid> &factory,
982 const std::string &fileName,
983 DataFlagArg boundarySegmentData,
984 DataArg elementData,
985 bool verbose=true)
986 {
987 doRead(
988 factory, fileName,
989 boundarySegmentData.data_
990 ? *boundarySegmentData.data_ : discarded(std::vector<int>{}),
991 elementData.data_
992 ? *elementData.data_ : discarded(std::vector<int>{}),
993 verbose,
994 boundarySegmentData.flag_ || boundarySegmentData.data_
995 );
996 }
997
1018 static void read (Dune::GridFactory<Grid>& factory,
1019 const std::string& fileName,
1020 std::vector<int>& boundarySegmentToPhysicalEntity,
1021 std::vector<int>& elementToPhysicalEntity,
1022 bool verbose, bool insertBoundarySegments)
1023 {
1024 doRead(
1025 factory, fileName, boundarySegmentToPhysicalEntity,
1026 elementToPhysicalEntity, verbose, insertBoundarySegments
1027 );
1028 }
1029
1031 //\{
1032
1033 [[deprecated("Will be removed after 2.8. Either use other constructors or use static methods without constructing an object")]]
1034 GmshReader() = default;
1035
1037
1038 static constexpr Opts defaultOpts =
1039 Opts::verbose | Opts::insertBoundarySegments | Opts::readElementData | Opts::readBoundaryData;
1040
1042
1065 GmshReader(const std::string& fileName,
1067 {
1068 gridFactory_ = std::make_unique<Dune::GridFactory<Grid>>();
1069 readGridFile(fileName, *gridFactory_, options);
1070 }
1071
1079 GmshReader(const std::string& fileName, GridFactory<Grid>& factory,
1081 {
1082 readGridFile(fileName, factory, options);
1083 }
1084
1086 const std::vector<int>& elementData () const
1087 {
1088 checkElementData();
1089 return elementIndexToGmshPhysicalEntity_;
1090 }
1091
1093 const std::vector<int>& boundaryData () const
1094 {
1095 checkBoundaryData();
1096 return boundarySegmentIndexToGmshPhysicalEntity_;
1097 }
1098
1103 bool hasElementData () const
1104 { return hasElementData_ && !extractedElementData_; }
1105
1110 bool hasBoundaryData () const
1111 { return hasBoundaryData_ && !extractedBoundaryData_; }
1112
1114 std::vector<int> extractElementData ()
1115 {
1116 checkElementData();
1117 extractedElementData_ = true;
1118 return std::move(elementIndexToGmshPhysicalEntity_);
1119 }
1120
1122 std::vector<int> extractBoundaryData ()
1123 {
1124 checkBoundaryData();
1125 extractedBoundaryData_ = true;
1126 return std::move(boundarySegmentIndexToGmshPhysicalEntity_);
1127 }
1128
1130 std::unique_ptr<Grid> createGrid ()
1131 {
1132 if (!gridFactory_)
1133 DUNE_THROW(Dune::InvalidStateException,
1134 "This GmshReader has been constructed with a Dune::GridFactory. "
1135 << "This grid factory has been filled with all information to create a grid. "
1136 << "Please use this factory to create the grid by calling factory.createGrid(). "
1137 << "Alternatively use the constructor without passing the factory in combination with this member function."
1138 );
1139
1140 return gridFactory_->createGrid();
1141 }
1142
1143 //\}
1144
1145 private:
1146 void checkElementData () const
1147 {
1148 if (!hasElementData_)
1149 DUNE_THROW(Dune::InvalidStateException,
1150 "This GmshReader has been constructed without the option 'readElementData'. "
1151 << "Please enable reading element data by passing the option 'Gmsh::ReaderOpts::readElementData' "
1152 << "to the constructor of this class."
1153 );
1154
1155 if (extractedElementData_)
1156 DUNE_THROW(Dune::InvalidStateException,
1157 "The element data has already been extracted from this GmshReader "
1158 << "via a function call to reader.extractElementData(). Use the extraced data or "
1159 << "read the grid data from file again by constructing a new reader."
1160 );
1161 }
1162
1163 void checkBoundaryData () const
1164 {
1165 if (!hasBoundaryData_)
1166 DUNE_THROW(Dune::InvalidStateException,
1167 "This GmshReader has been constructed without the option 'readBoundaryData'. "
1168 << "Please enable reading boundary data by passing the option 'Gmsh::ReaderOpts::readBoundaryData' "
1169 << "to the constructor of this class."
1170 );
1171
1172 if (extractedBoundaryData_)
1173 DUNE_THROW(Dune::InvalidStateException,
1174 "The boundary data has already been extracted from this GmshReader "
1175 << "via a function call to reader.extractBoundaryData(). Use the extraced data or "
1176 << "read the grid data from file again by constructing a new reader."
1177 );
1178 }
1179
1180 void readGridFile (const std::string& fileName, GridFactory<Grid>& factory, Gmsh::ReaderOptions options)
1181 {
1182 const bool verbose = options & Opts::verbose;
1183 const bool insertBoundarySegments = options & Opts::insertBoundarySegments;
1184 const bool readBoundaryData = options & Opts::readBoundaryData;
1185 const bool readElementData = options & Opts::readElementData;
1186
1187 doRead(
1188 factory, fileName, boundarySegmentIndexToGmshPhysicalEntity_,
1189 elementIndexToGmshPhysicalEntity_, verbose,
1190 readBoundaryData || insertBoundarySegments
1191 );
1192
1193 // clear unwanted data
1194 if (!readBoundaryData)
1195 boundarySegmentIndexToGmshPhysicalEntity_ = std::vector<int>{};
1196 if (!readElementData)
1197 elementIndexToGmshPhysicalEntity_ = std::vector<int>{};
1198
1199 hasElementData_ = readElementData;
1200 hasBoundaryData_ = readBoundaryData;
1201 }
1202
1203 std::unique_ptr<Dune::GridFactory<Grid>> gridFactory_;
1204
1205 std::vector<int> elementIndexToGmshPhysicalEntity_;
1206 std::vector<int> boundarySegmentIndexToGmshPhysicalEntity_;
1207
1208 bool hasElementData_;
1209 bool hasBoundaryData_;
1210
1211 // for better error messages, we keep track of these separately
1212 bool extractedElementData_ = false;
1213 bool extractedBoundaryData_ = false;
1214 };
1215
1218} // namespace Dune
1219
1220#endif
Base class for grid boundary segments of arbitrary geometry.
ReaderOptions
Option for the Gmsh mesh file reader.
Definition: gmshreader.hh:756
void swap(Dune::PersistentContainer< G, T > &a, Dune::PersistentContainer< G, T > &b)
Definition: utility/persistentcontainer.hh:83
Include standard header files.
Definition: agrid.hh:60
static const int dimWorld
Definition: misc.hh:46
ALBERTA REAL_D GlobalVector
Definition: misc.hh:50
constexpr bool operator&(ReaderOptions a, ReaderOptions b)
query operator for reader options (is b set in a)
Definition: gmshreader.hh:772
constexpr ReaderOptions operator|(ReaderOptions a, ReaderOptions b)
composition operator for reader options
Definition: gmshreader.hh:764
int renumber(const Dune::GeometryType &t, int i)
renumber VTK <-> Dune
Definition: common.hh:186
@ line
Definition: common.hh:134
@ pyramid
Definition: common.hh:141
@ quadrilateral
Definition: common.hh:137
@ tetrahedron
Definition: common.hh:138
@ prism
Definition: common.hh:140
@ hexahedron
Definition: common.hh:139
@ triangle
Definition: common.hh:135
Base class for classes implementing geometries of boundary segments.
Definition: boundarysegment.hh:94
Communication comm() const
Return the Communication used by the grid factory.
Definition: common/gridfactory.hh:297
Provide a generic factory class for unstructured grids.
Definition: common/gridfactory.hh:314
virtual std::unique_ptr< GridType > createGrid()
Finalize grid creation and hand over the grid.
Definition: common/gridfactory.hh:372
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
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:510
static const int dimWorld
Definition: gmshreader.hh:320
Dune::GridFactory< GridType > & factory
Definition: gmshreader.hh:305
std::vector< int > & boundaryIdMap()
Definition: gmshreader.hh:357
std::vector< int > & elementIndexMap()
Definition: gmshreader.hh:362
unsigned int number_of_real_vertices
Definition: gmshreader.hh:308
void boundarysegment_insert(const std::vector< FieldVector< double, 3 > > &nodes, const E &elementDofs, const V &vertices)
Definition: gmshreader.hh:573
GmshReaderParser(Dune::GridFactory< GridType > &_factory, bool v, bool i)
Definition: gmshreader.hh:354
int element_count
Definition: gmshreader.hh:310
void read(const std::string &f)
Definition: gmshreader.hh:367
void skipline(FILE *file)
Definition: gmshreader.hh:344
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:330
std::vector< int > element_index_to_physical_entity
Definition: gmshreader.hh:316
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:598
static const int dim
Definition: gmshreader.hh:319
FieldVector< double, dimWorld > GlobalVector
Definition: gmshreader.hh:324
std::string fileName
Definition: gmshreader.hh:313
int boundary_element_count
Definition: gmshreader.hh:309
void boundarysegment_insert(const V &, const E &, const V2 &)
Definition: gmshreader.hh:562
bool verbose
Definition: gmshreader.hh:306
std::vector< int > boundary_id_to_physical_entity
Definition: gmshreader.hh:315
char buf[512]
Definition: gmshreader.hh:312
bool insert_boundary_segments
Definition: gmshreader.hh:307
Read Gmsh mesh file.
Definition: gmshreader.hh:805
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:931
const std::vector< int > & elementData() const
Access element data (maps element index to Gmsh physical entity)
Definition: gmshreader.hh:1086
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:981
static std::unique_ptr< Grid > read(const std::string &fileName, bool verbose=true, bool insertBoundarySegments=true)
Definition: gmshreader.hh:902
static void read(Dune::GridFactory< Grid > &factory, const std::string &fileName, bool verbose=true, bool insertBoundarySegments=true)
Definition: gmshreader.hh:948
GridType Grid
Definition: gmshreader.hh:894
std::unique_ptr< Grid > createGrid()
Create the grid.
Definition: gmshreader.hh:1130
std::vector< int > extractBoundaryData()
Erase boundary data from reader and return the data.
Definition: gmshreader.hh:1122
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:1018
bool hasElementData() const
If element data is available.
Definition: gmshreader.hh:1103
bool hasBoundaryData() const
If boundary data is available.
Definition: gmshreader.hh:1110
static constexpr Opts defaultOpts
Definition: gmshreader.hh:1038
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:1079
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:1065
std::vector< int > extractElementData()
Erase element data from reader and return the data.
Definition: gmshreader.hh:1114
const std::vector< int > & boundaryData() const
Access boundary data (maps boundary segment index to Gmsh physical entity)
Definition: gmshreader.hh:1093
GmshReader()=default
Dynamic Gmsh reader interface.
Provide a generic factory class for unstructured grids.