dune-grid 2.9.0
Loading...
Searching...
No Matches
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
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 // 2.2 is not representable as float and leads to problems on i386
395 // Hence we use >= 2.00001.
396 if( (version_number < 2.0) || (version_number >= 2.20001) ) // 2.2 is not representable as float and leads to problems on i386
397 DUNE_THROW(Dune::IOError, "can only read Gmsh version 2 files");
398 if (verbose) std::cout << "version " << version_number << " Gmsh file detected" << std::endl;
399 readfile(file,1,"%s\n",buf);
400 if (strcmp(buf,"$EndMeshFormat")!=0)
401 DUNE_THROW(Dune::IOError, "expected $EndMeshFormat");
402
403 // node section
404 int number_of_nodes;
405
406 readfile(file,1,"%s\n",buf);
407 if (strcmp(buf,"$Nodes")!=0)
408 DUNE_THROW(Dune::IOError, "expected $Nodes");
409 readfile(file,1,"%d\n",&number_of_nodes);
410 if (verbose) std::cout << "file contains " << number_of_nodes << " nodes" << std::endl;
411
412 // read nodes
413 // The '+1' is due to the fact that gmsh numbers node starting from 1 rather than from 0
414 std::vector< GlobalVector > nodes( number_of_nodes+1 );
415 {
416 int id;
417 double x[ 3 ];
418 for( int i = 1; i <= number_of_nodes; ++i )
419 {
420 readfile(file,4, "%d %lg %lg %lg\n", &id, &x[ 0 ], &x[ 1 ], &x[ 2 ] );
421
422 if (id > number_of_nodes) {
423 DUNE_THROW(Dune::IOError,
424 "Only dense sequences of node indices are currently supported (node index "
425 << id << " is invalid).");
426 }
427
428 // just store node position
429 for( int j = 0; j < dimWorld; ++j )
430 nodes[ id ][ j ] = x[ j ];
431 }
432 readfile(file,1,"%s\n",buf);
433 if (strcmp(buf,"$EndNodes")!=0)
434 DUNE_THROW(Dune::IOError, "expected $EndNodes");
435 }
436
437 // element section
438 readfile(file,1,"%s\n",buf);
439 if (strcmp(buf,"$Elements")!=0)
440 DUNE_THROW(Dune::IOError, "expected $Elements");
441 int number_of_elements;
442 readfile(file,1,"%d\n",&number_of_elements);
443 if (verbose) std::cout << "file contains " << number_of_elements << " elements" << std::endl;
444
445 //=========================================
446 // Pass 1: Select and insert those vertices in the file that
447 // actually occur as corners of an element.
448 //=========================================
449
450 off_t section_element_offset = ftello(file);
451 std::map<int,unsigned int> renumber;
452 for (int i=1; i<=number_of_elements; i++)
453 {
454 int id, elm_type, number_of_tags;
455 readfile(file,3,"%d %d %d ",&id,&elm_type,&number_of_tags);
456 for (int k=1; k<=number_of_tags; k++)
457 {
458 int blub;
459 readfile(file,1,"%d ",&blub);
460 // k == 1: physical entity (not used here)
461 // k == 2: elementary entity (not used here either)
462 // if version_number < 2.2:
463 // k == 3: mesh partition 0
464 // else
465 // k == 3: number of mesh partitions
466 // k => 4: mesh partition k-4
467 }
468 pass1HandleElement(file, elm_type, renumber, nodes);
469 }
470 if (verbose) std::cout << "number of real vertices = " << number_of_real_vertices << std::endl;
471 if (verbose) std::cout << "number of boundary elements = " << boundary_element_count << std::endl;
472 if (verbose) std::cout << "number of elements = " << element_count << std::endl;
473 readfile(file,1,"%s\n",buf);
474 if (strcmp(buf,"$EndElements")!=0)
475 DUNE_THROW(Dune::IOError, "expected $EndElements");
478
479 //==============================================
480 // Pass 2: Insert boundary segments and elements
481 //==============================================
482
483 fseeko(file, section_element_offset, SEEK_SET);
485 element_count = 0;
486 for (int i=1; i<=number_of_elements; i++)
487 {
488 int id, elm_type, number_of_tags;
489 readfile(file,3,"%d %d %d ",&id,&elm_type,&number_of_tags);
490 int physical_entity = -1;
491
492 for (int k=1; k<=number_of_tags; k++)
493 {
494 int blub;
495 readfile(file,1,"%d ",&blub);
496 if (k==1) physical_entity = blub;
497 }
498 pass2HandleElement(file, elm_type, renumber, nodes, physical_entity);
499 }
500 readfile(file,1,"%s\n",buf);
501 if (strcmp(buf,"$EndElements")!=0)
502 DUNE_THROW(Dune::IOError, "expected $EndElements");
503
504 fclose(file);
505 }
506
512 void pass1HandleElement(FILE* file, const int elm_type,
513 std::map<int,unsigned int> & renumber,
514 const std::vector< GlobalVector > & nodes)
515 {
516 // some data about gmsh elements
517 const int nDofs[16] = {-1, 2, 3, 4, 4, 8, 6, 5, 3, 6, -1, 10, -1, -1, -1, 1};
518 const int nVertices[16] = {-1, 2, 3, 4, 4, 8, 6, 5, 2, 3, -1, 4, -1, -1, -1, 1};
519 const int elementDim[16] = {-1, 1, 2, 2, 3, 3, 3, 3, 1, 2, -1, 3, -1, -1, -1, 0};
520
521 // test whether we support the element type
522 if ( not (elm_type > 0 && elm_type <= 15 // index in suitable range?
523 && (elementDim[elm_type] == dim || elementDim[elm_type] == (dim-1) ) ) ) // real element or boundary element?
524 {
525 skipline(file); // skip rest of line if element is unknown
526 return;
527 }
528
529 // The format string for parsing is n times '%d' in a row
530 std::string formatString = "%d";
531 for (int i=1; i<nDofs[elm_type]; i++)
532 formatString += " %d";
533 formatString += "\n";
534
535 // '10' is the largest number of dofs we may encounter in a .msh file
536 std::vector<int> elementDofs(10);
537
538 readfile(file,nDofs[elm_type], formatString.c_str(),
539 &(elementDofs[0]),&(elementDofs[1]),&(elementDofs[2]),
540 &(elementDofs[3]),&(elementDofs[4]),&(elementDofs[5]),
541 &(elementDofs[6]),&(elementDofs[7]),&(elementDofs[8]),
542 &(elementDofs[9]));
543
544 // insert each vertex if it hasn't been inserted already
545 for (int i=0; i<nVertices[elm_type]; i++)
546 if (renumber.find(elementDofs[i])==renumber.end())
547 {
548 renumber[elementDofs[i]] = number_of_real_vertices++;
549 factory.insertVertex(nodes[elementDofs[i]]);
550 }
551
552 // count elements and boundary elements
553 if (elementDim[elm_type] == dim)
555 else
557
558 }
559
560
561
562 // generic-case: This is not supposed to be used at runtime.
563 template <class E, class V, class V2>
565 const V&,
566 const E&,
567 const V2&
568 )
569 {
570 DUNE_THROW(Dune::IOError, "tried to create a 3D boundary segment in a non-3D Grid");
571 }
572
573 // 3d-case:
574 template <class E, class V>
576 const std::vector<FieldVector<double, 3> >& nodes,
577 const E& elementDofs,
578 const V& vertices
579 )
580 {
581 std::array<FieldVector<double,dimWorld>, 6> v;
582 for (int i=0; i<6; i++)
583 for (int j=0; j<dimWorld; j++)
584 v[i][j] = nodes[elementDofs[i]][j];
585
586 BoundarySegment<dim,dimWorld>* newBoundarySegment
587 = (BoundarySegment<dim,dimWorld>*) new GmshReaderQuadraticBoundarySegment< 3, 3 >( v[0], v[1], v[2],
588 v[3], v[4], v[5] );
589
590 factory.insertBoundarySegment( vertices,
591 std::shared_ptr<BoundarySegment<dim,dimWorld> >(newBoundarySegment) );
592 }
593
594
595
600 virtual void pass2HandleElement(FILE* file, const int elm_type,
601 std::map<int,unsigned int> & renumber,
602 const std::vector< GlobalVector > & nodes,
603 const int physical_entity)
604 {
605 // some data about gmsh elements
606 const int nDofs[16] = {-1, 2, 3, 4, 4, 8, 6, 5, 3, 6, -1, 10, -1, -1, -1, 1};
607 const int nVertices[16] = {-1, 2, 3, 4, 4, 8, 6, 5, 2, 3, -1, 4, -1, -1, -1, 1};
608 const int elementDim[16] = {-1, 1, 2, 2, 3, 3, 3, 3, 1, 2, -1, 3, -1, -1, -1, 0};
609
610 // test whether we support the element type
611 if ( not (elm_type > 0 && elm_type <= 15 // index in suitable range?
612 && (elementDim[elm_type] == dim || elementDim[elm_type] == (dim-1) ) ) ) // real element or boundary element?
613 {
614 skipline(file); // skip rest of line if element is unknown
615 return;
616 }
617
618 // The format string for parsing is n times '%d' in a row
619 std::string formatString = "%d";
620 for (int i=1; i<nDofs[elm_type]; i++)
621 formatString += " %d";
622 formatString += "\n";
623
624 // '10' is the largest number of dofs we may encounter in a .msh file
625 std::vector<int> elementDofs(10);
626
627 readfile(file,nDofs[elm_type], formatString.c_str(),
628 &(elementDofs[0]),&(elementDofs[1]),&(elementDofs[2]),
629 &(elementDofs[3]),&(elementDofs[4]),&(elementDofs[5]),
630 &(elementDofs[6]),&(elementDofs[7]),&(elementDofs[8]),
631 &(elementDofs[9]));
632
633 // correct differences between gmsh and Dune in the local vertex numbering
634 switch (elm_type)
635 {
636 case 3 : // 4-node quadrilateral
637 std::swap(elementDofs[2],elementDofs[3]);
638 break;
639 case 5 : // 8-node hexahedron
640 std::swap(elementDofs[2],elementDofs[3]);
641 std::swap(elementDofs[6],elementDofs[7]);
642 break;
643 case 7 : // 5-node pyramid
644 std::swap(elementDofs[2],elementDofs[3]);
645 break;
646 }
647
648 // renumber corners to account for the explicitly given vertex
649 // numbering in the file
650 std::vector<unsigned int> vertices(nVertices[elm_type]);
651
652 for (int i=0; i<nVertices[elm_type]; i++)
653 vertices[i] = renumber[elementDofs[i]];
654
655 // If it is an element, insert it as such
656 if (elementDim[elm_type] == dim) {
657
658 switch (elm_type)
659 {
660 case 1 : // 2-node line
661 factory.insertElement(Dune::GeometryTypes::line,vertices);
662 break;
663 case 2 : // 3-node triangle
664 factory.insertElement(Dune::GeometryTypes::triangle,vertices);
665 break;
666 case 3 : // 4-node quadrilateral
667 factory.insertElement(Dune::GeometryTypes::quadrilateral,vertices);
668 break;
669 case 4 : // 4-node tetrahedron
670 factory.insertElement(Dune::GeometryTypes::tetrahedron,vertices);
671 break;
672 case 5 : // 8-node hexahedron
673 factory.insertElement(Dune::GeometryTypes::hexahedron,vertices);
674 break;
675 case 6 : // 6-node prism
676 factory.insertElement(Dune::GeometryTypes::prism,vertices);
677 break;
678 case 7 : // 5-node pyramid
679 factory.insertElement(Dune::GeometryTypes::pyramid,vertices);
680 break;
681 case 9 : // 6-node triangle
682 factory.insertElement(Dune::GeometryTypes::triangle,vertices);
683 break;
684 case 11 : // 10-node tetrahedron
685 factory.insertElement(Dune::GeometryTypes::tetrahedron,vertices);
686 break;
687 }
688
689 } else {
690 // it must be a boundary segment then
692
693 switch (elm_type)
694 {
695 case 1 : // 2-node line
696 factory.insertBoundarySegment(vertices);
697 break;
698
699 case 2 : // 3-node triangle
700 factory.insertBoundarySegment(vertices);
701 break;
702
703 case 3 : // 4-node quadrilateral
704 factory.insertBoundarySegment(vertices);
705 break;
706
707 case 15 : // 1-node point
708 factory.insertBoundarySegment(vertices);
709 break;
710
711 case 8 : { // 3-node line
712 std::array<FieldVector<double,dimWorld>, 3> v;
713 for (int i=0; i<dimWorld; i++) {
714 v[0][i] = nodes[elementDofs[0]][i];
715 v[1][i] = nodes[elementDofs[2]][i]; // yes, the renumbering is intended!
716 v[2][i] = nodes[elementDofs[1]][i];
717 }
718 BoundarySegment<dim,dimWorld>* newBoundarySegment
719 = (BoundarySegment<dim,dimWorld>*) new GmshReaderQuadraticBoundarySegment< 2, dimWorld >(v[0], v[1], v[2]);
720 factory.insertBoundarySegment(vertices,
721 std::shared_ptr<BoundarySegment<dim,dimWorld> >(newBoundarySegment));
722 break;
723 }
724 case 9 : { // 6-node triangle
725 boundarysegment_insert(nodes, elementDofs, vertices);
726 break;
727 }
728 default: {
729 DUNE_THROW(Dune::IOError, "GmshReader does not support using element-type " << elm_type << " for boundary segments");
730 break;
731 }
732
733 }
734
735 }
736 }
737
738 // count elements and boundary elements
739 if (elementDim[elm_type] == dim) {
742 } else {
745 }
746
747 }
748
749 };
750
751 namespace Gmsh {
757 enum class ReaderOptions
758 {
759 verbose = 1,
761 readElementData = 4,
763 };
764
767 {
768 return static_cast<ReaderOptions>(
769 static_cast<int>(a) | static_cast<int>(b)
770 );
771 }
772
775 {
776 return static_cast<int>(a) & static_cast<int>(b);
777 }
778
779 } // end namespace Gmsh
780
805 template<typename GridType>
807 {
809
828 static void doRead(Dune::GridFactory<GridType> &factory,
829 const std::string &fileName,
830 std::vector<int>& boundarySegmentToPhysicalEntity,
831 std::vector<int>& elementToPhysicalEntity,
832 bool verbose, bool insertBoundarySegments)
833 {
834 // register boundary segment to boundary segment factory for possible load balancing
835 // this needs to be done on all cores since the type might not be known otherwise
836 GmshReaderQuadraticBoundarySegment< Grid::dimension, Grid::dimensionworld >::registerFactory();
837
838#ifndef NDEBUG
839 // check that this method is called on all cores
840 factory.comm().barrier();
841#endif
842
843 // create parse object and read grid on process 0
844 if (factory.comm().rank() == 0)
845 {
846 GmshReaderParser<Grid> parser(factory,verbose,insertBoundarySegments);
847 parser.read(fileName);
848
849 boundarySegmentToPhysicalEntity = std::move(parser.boundaryIdMap());
850 elementToPhysicalEntity = std::move(parser.elementIndexMap());
851 }
852 else
853 {
854 boundarySegmentToPhysicalEntity = {};
855 elementToPhysicalEntity = {};
856 }
857 }
858
860
879 template<class T>
880 static T &discarded(T &&value) { return value; }
881
882 struct DataArg {
883 std::vector<int> *data_ = nullptr;
884 DataArg(std::vector<int> &data) : data_(&data) {}
885 DataArg(const decltype(std::ignore)&) {}
886 DataArg() = default;
887 };
888
889 struct DataFlagArg : DataArg {
890 bool flag_ = false;
891 using DataArg::DataArg;
892 DataFlagArg(bool flag) : flag_(flag) {}
893 };
894
895 public:
896 typedef GridType Grid;
897
904 static std::unique_ptr<Grid> read (const std::string& fileName, bool verbose = true, bool insertBoundarySegments=true)
905 {
906 // make a grid factory
908
909 read(factory, fileName, verbose, insertBoundarySegments);
910
911 return factory.createGrid();
912 }
913
933 static std::unique_ptr<Grid> read (const std::string& fileName,
934 std::vector<int>& boundarySegmentToPhysicalEntity,
935 std::vector<int>& elementToPhysicalEntity,
936 bool verbose = true, bool insertBoundarySegments=true)
937 {
938 // make a grid factory
940
941 doRead(
942 factory, fileName, boundarySegmentToPhysicalEntity,
943 elementToPhysicalEntity, verbose, insertBoundarySegments
944 );
945
946 return factory.createGrid();
947 }
948
950 static void read (Dune::GridFactory<Grid>& factory, const std::string& fileName,
951 bool verbose = true, bool insertBoundarySegments=true)
952 {
953 doRead(
954 factory, fileName, discarded(std::vector<int>{}),
955 discarded(std::vector<int>{}), verbose, insertBoundarySegments
956 );
957 }
958
960
983 static void read (Dune::GridFactory<Grid> &factory,
984 const std::string &fileName,
985 DataFlagArg boundarySegmentData,
986 DataArg elementData,
987 bool verbose=true)
988 {
989 doRead(
990 factory, fileName,
991 boundarySegmentData.data_
992 ? *boundarySegmentData.data_ : discarded(std::vector<int>{}),
993 elementData.data_
994 ? *elementData.data_ : discarded(std::vector<int>{}),
995 verbose,
996 boundarySegmentData.flag_ || boundarySegmentData.data_
997 );
998 }
999
1020 static void read (Dune::GridFactory<Grid>& factory,
1021 const std::string& fileName,
1022 std::vector<int>& boundarySegmentToPhysicalEntity,
1023 std::vector<int>& elementToPhysicalEntity,
1024 bool verbose, bool insertBoundarySegments)
1025 {
1026 doRead(
1027 factory, fileName, boundarySegmentToPhysicalEntity,
1028 elementToPhysicalEntity, verbose, insertBoundarySegments
1029 );
1030 }
1031
1033 //\{
1034
1035 [[deprecated("Will be removed after 2.8. Either use other constructors or use static methods without constructing an object")]]
1036 GmshReader() = default;
1037
1039
1040 static constexpr Opts defaultOpts =
1041 Opts::verbose | Opts::insertBoundarySegments | Opts::readElementData | Opts::readBoundaryData;
1042
1044
1067 GmshReader(const std::string& fileName,
1069 {
1070 gridFactory_ = std::make_unique<Dune::GridFactory<Grid>>();
1071 readGridFile(fileName, *gridFactory_, options);
1072 }
1073
1081 GmshReader(const std::string& fileName, GridFactory<Grid>& factory,
1083 {
1084 readGridFile(fileName, factory, options);
1085 }
1086
1088 const std::vector<int>& elementData () const
1089 {
1090 checkElementData();
1091 return elementIndexToGmshPhysicalEntity_;
1092 }
1093
1095 const std::vector<int>& boundaryData () const
1096 {
1097 checkBoundaryData();
1098 return boundarySegmentIndexToGmshPhysicalEntity_;
1099 }
1100
1105 bool hasElementData () const
1106 { return hasElementData_ && !extractedElementData_; }
1107
1112 bool hasBoundaryData () const
1113 { return hasBoundaryData_ && !extractedBoundaryData_; }
1114
1116 std::vector<int> extractElementData ()
1117 {
1118 checkElementData();
1119 extractedElementData_ = true;
1120 return std::move(elementIndexToGmshPhysicalEntity_);
1121 }
1122
1124 std::vector<int> extractBoundaryData ()
1125 {
1126 checkBoundaryData();
1127 extractedBoundaryData_ = true;
1128 return std::move(boundarySegmentIndexToGmshPhysicalEntity_);
1129 }
1130
1132 std::unique_ptr<Grid> createGrid ()
1133 {
1134 if (!gridFactory_)
1135 DUNE_THROW(Dune::InvalidStateException,
1136 "This GmshReader has been constructed with a Dune::GridFactory. "
1137 << "This grid factory has been filled with all information to create a grid. "
1138 << "Please use this factory to create the grid by calling factory.createGrid(). "
1139 << "Alternatively use the constructor without passing the factory in combination with this member function."
1140 );
1141
1142 return gridFactory_->createGrid();
1143 }
1144
1145 //\}
1146
1147 private:
1148 void checkElementData () const
1149 {
1150 if (!hasElementData_)
1151 DUNE_THROW(Dune::InvalidStateException,
1152 "This GmshReader has been constructed without the option 'readElementData'. "
1153 << "Please enable reading element data by passing the option 'Gmsh::ReaderOpts::readElementData' "
1154 << "to the constructor of this class."
1155 );
1156
1157 if (extractedElementData_)
1158 DUNE_THROW(Dune::InvalidStateException,
1159 "The element data has already been extracted from this GmshReader "
1160 << "via a function call to reader.extractElementData(). Use the extraced data or "
1161 << "read the grid data from file again by constructing a new reader."
1162 );
1163 }
1164
1165 void checkBoundaryData () const
1166 {
1167 if (!hasBoundaryData_)
1168 DUNE_THROW(Dune::InvalidStateException,
1169 "This GmshReader has been constructed without the option 'readBoundaryData'. "
1170 << "Please enable reading boundary data by passing the option 'Gmsh::ReaderOpts::readBoundaryData' "
1171 << "to the constructor of this class."
1172 );
1173
1174 if (extractedBoundaryData_)
1175 DUNE_THROW(Dune::InvalidStateException,
1176 "The boundary data has already been extracted from this GmshReader "
1177 << "via a function call to reader.extractBoundaryData(). Use the extraced data or "
1178 << "read the grid data from file again by constructing a new reader."
1179 );
1180 }
1181
1182 void readGridFile (const std::string& fileName, GridFactory<Grid>& factory, Gmsh::ReaderOptions options)
1183 {
1184 const bool verbose = options & Opts::verbose;
1185 const bool insertBoundarySegments = options & Opts::insertBoundarySegments;
1186 const bool readBoundaryData = options & Opts::readBoundaryData;
1187 const bool readElementData = options & Opts::readElementData;
1188
1189 doRead(
1190 factory, fileName, boundarySegmentIndexToGmshPhysicalEntity_,
1191 elementIndexToGmshPhysicalEntity_, verbose,
1192 readBoundaryData || insertBoundarySegments
1193 );
1194
1195 // clear unwanted data
1196 if (!readBoundaryData)
1197 boundarySegmentIndexToGmshPhysicalEntity_ = std::vector<int>{};
1198 if (!readElementData)
1199 elementIndexToGmshPhysicalEntity_ = std::vector<int>{};
1200
1201 hasElementData_ = readElementData;
1202 hasBoundaryData_ = readBoundaryData;
1203 }
1204
1205 std::unique_ptr<Dune::GridFactory<Grid>> gridFactory_;
1206
1207 std::vector<int> elementIndexToGmshPhysicalEntity_;
1208 std::vector<int> boundarySegmentIndexToGmshPhysicalEntity_;
1209
1210 bool hasElementData_;
1211 bool hasBoundaryData_;
1212
1213 // for better error messages, we keep track of these separately
1214 bool extractedElementData_ = false;
1215 bool extractedBoundaryData_ = false;
1216 };
1217
1220} // namespace Dune
1221
1222#endif
Base class for grid boundary segments of arbitrary geometry.
ReaderOptions
Option for the Gmsh mesh file reader.
Definition gmshreader.hh:758
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
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:774
constexpr ReaderOptions operator|(ReaderOptions a, ReaderOptions b)
composition operator for reader options
Definition gmshreader.hh:766
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:512
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:575
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:600
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:564
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:807
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:933
const std::vector< int > & elementData() const
Access element data (maps element index to Gmsh physical entity)
Definition gmshreader.hh:1088
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:983
static std::unique_ptr< Grid > read(const std::string &fileName, bool verbose=true, bool insertBoundarySegments=true)
Definition gmshreader.hh:904
static void read(Dune::GridFactory< Grid > &factory, const std::string &fileName, bool verbose=true, bool insertBoundarySegments=true)
Definition gmshreader.hh:950
GridType Grid
Definition gmshreader.hh:896
std::unique_ptr< Grid > createGrid()
Create the grid.
Definition gmshreader.hh:1132
std::vector< int > extractBoundaryData()
Erase boundary data from reader and return the data.
Definition gmshreader.hh:1124
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:1020
bool hasElementData() const
If element data is available.
Definition gmshreader.hh:1105
bool hasBoundaryData() const
If boundary data is available.
Definition gmshreader.hh:1112
static constexpr Opts defaultOpts
Definition gmshreader.hh:1040
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:1081
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:1067
std::vector< int > extractElementData()
Erase element data from reader and return the data.
Definition gmshreader.hh:1116
const std::vector< int > & boundaryData() const
Access boundary data (maps boundary segment index to Gmsh physical entity)
Definition gmshreader.hh:1095
GmshReader()=default
Dynamic Gmsh reader interface.
Provide a generic factory class for unstructured grids.