DUNE PDELab (git)

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
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_,
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 auto pos = std::ftell(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 // 2.2 is not representable as float and leads to problems on i386
406 // Hence we use >= 2.00001.
407 if( (version_number < 2.0) || (version_number >= 2.20001) ) // 2.2 is not representable as float and leads to problems on i386
408 DUNE_THROW(Dune::IOError, "can only read Gmsh version 2 files");
409 if (verbose) std::cout << "version " << version_number << " Gmsh file detected" << std::endl;
410 readfile(file,1,"%s\n",buf);
411 if (strcmp(buf,"$EndMeshFormat")!=0)
412 DUNE_THROW(Dune::IOError, "expected $EndMeshFormat");
413
414 // physical name section
415 physical_entity_names.clear();
416 readfile(file,1,"%s\n",buf);
417 if (strcmp(buf,"$PhysicalNames")==0) {
418 int number_of_names;
419 readfile(file,1,"%d\n",&number_of_names);
420 if (verbose) std::cout << "file contains " << number_of_names << " physical entities" << std::endl;
421 physical_entity_names.resize(number_of_names);
422 std::string buf_name;
423 for( int i = 0; i < number_of_names; ++i ) {
424 int dim, id;
425 readfile(file,3, "%d %d %s\n", &dim, &id, buf);
426 buf_name.assign(buf);
427 auto begin = buf_name.find_first_of('\"') + 1;
428 auto end = buf_name.find_last_of('\"') - begin;
429 physical_entity_names[id-1].assign(buf_name.substr(begin, end));
430 }
431 readfile(file,1,"%s\n",buf);
432 if (strcmp(buf,"$EndPhysicalNames")!=0)
433 DUNE_THROW(Dune::IOError, "expected $EndPhysicalNames");
434 readfile(file,1,"%s\n",buf);
435 }
436
437 // node section
438 int number_of_nodes;
439
440 if (strcmp(buf,"$Nodes")!=0)
441 DUNE_THROW(Dune::IOError, "expected $Nodes");
442 readfile(file,1,"%d\n",&number_of_nodes);
443 if (verbose) std::cout << "file contains " << number_of_nodes << " nodes" << std::endl;
444
445 // read nodes
446 // The '+1' is due to the fact that gmsh numbers node starting from 1 rather than from 0
447 std::vector< GlobalVector > nodes( number_of_nodes+1 );
448 {
449 int id;
450 double x[ 3 ];
451 for( int i = 1; i <= number_of_nodes; ++i )
452 {
453 readfile(file,4, "%d %lg %lg %lg\n", &id, &x[ 0 ], &x[ 1 ], &x[ 2 ] );
454
455 if (id > number_of_nodes) {
457 "Only dense sequences of node indices are currently supported (node index "
458 << id << " is invalid).");
459 }
460
461 // just store node position
462 for( int j = 0; j < dimWorld; ++j )
463 nodes[ id ][ j ] = x[ j ];
464 }
465 readfile(file,1,"%s\n",buf);
466 if (strcmp(buf,"$EndNodes")!=0)
467 DUNE_THROW(Dune::IOError, "expected $EndNodes");
468 }
469
470 // element section
471 readfile(file,1,"%s\n",buf);
472 if (strcmp(buf,"$Elements")!=0)
473 DUNE_THROW(Dune::IOError, "expected $Elements");
474 int number_of_elements;
475 readfile(file,1,"%d\n",&number_of_elements);
476 if (verbose) std::cout << "file contains " << number_of_elements << " elements" << std::endl;
477
478 //=========================================
479 // Pass 1: Select and insert those vertices in the file that
480 // actually occur as corners of an element.
481 //=========================================
482
483 auto section_element_offset = std::ftell(file);
484 std::map<int,unsigned int> renumber;
485 for (int i=1; i<=number_of_elements; i++)
486 {
487 int id, elm_type, number_of_tags;
488 readfile(file,3,"%d %d %d ",&id,&elm_type,&number_of_tags);
489 for (int k=1; k<=number_of_tags; k++)
490 {
491 int blub;
492 readfile(file,1,"%d ",&blub);
493 // k == 1: physical entity
494 // k == 2: elementary entity (not used here)
495 // if version_number < 2.2:
496 // k == 3: mesh partition 0
497 // else
498 // k == 3: number of mesh partitions
499 // k => 4: mesh partition k-4
500 }
501 pass1HandleElement(file, elm_type, renumber, nodes);
502 }
503 if (verbose) std::cout << "number of real vertices = " << number_of_real_vertices << std::endl;
504 if (verbose) std::cout << "number of boundary elements = " << boundary_element_count << std::endl;
505 if (verbose) std::cout << "number of elements = " << element_count << std::endl;
506 readfile(file,1,"%s\n",buf);
507 if (strcmp(buf,"$EndElements")!=0)
508 DUNE_THROW(Dune::IOError, "expected $EndElements");
509 boundary_id_to_physical_entity.resize(boundary_element_count);
510 element_index_to_physical_entity.resize(element_count);
511
512 //==============================================
513 // Pass 2: Insert boundary segments and elements
514 //==============================================
515
516 std::fseek(file, section_element_offset, SEEK_SET);
517 boundary_element_count = 0;
518 element_count = 0;
519 for (int i=1; i<=number_of_elements; i++)
520 {
521 int id, elm_type, number_of_tags;
522 readfile(file,3,"%d %d %d ",&id,&elm_type,&number_of_tags);
523 int physical_entity = -1;
524
525 for (int k=1; k<=number_of_tags; k++)
526 {
527 int blub;
528 readfile(file,1,"%d ",&blub);
529 if (k==1) physical_entity = blub;
530 }
531 pass2HandleElement(file, elm_type, renumber, nodes, physical_entity);
532 }
533 readfile(file,1,"%s\n",buf);
534 if (strcmp(buf,"$EndElements")!=0)
535 DUNE_THROW(Dune::IOError, "expected $EndElements");
536
537 fclose(file);
538 }
539
545 void pass1HandleElement(FILE* file, const int elm_type,
546 std::map<int,unsigned int> & renumber,
547 const std::vector< GlobalVector > & nodes)
548 {
549 // some data about gmsh elements
550 const int nDofs[16] = {-1, 2, 3, 4, 4, 8, 6, 5, 3, 6, -1, 10, -1, -1, -1, 1};
551 const int nVertices[16] = {-1, 2, 3, 4, 4, 8, 6, 5, 2, 3, -1, 4, -1, -1, -1, 1};
552 const int elementDim[16] = {-1, 1, 2, 2, 3, 3, 3, 3, 1, 2, -1, 3, -1, -1, -1, 0};
553
554 // test whether we support the element type
555 if ( not (elm_type > 0 && elm_type <= 15 // index in suitable range?
556 && (elementDim[elm_type] == dim || elementDim[elm_type] == (dim-1) ) ) ) // real element or boundary element?
557 {
558 skipline(file); // skip rest of line if element is unknown
559 return;
560 }
561
562 // The format string for parsing is n times '%d' in a row
563 std::string formatString = "%d";
564 for (int i=1; i<nDofs[elm_type]; i++)
565 formatString += " %d";
566 formatString += "\n";
567
568 // '10' is the largest number of dofs we may encounter in a .msh file
569 std::vector<int> elementDofs(10);
570
571 readfile(file,nDofs[elm_type], formatString.c_str(),
572 &(elementDofs[0]),&(elementDofs[1]),&(elementDofs[2]),
573 &(elementDofs[3]),&(elementDofs[4]),&(elementDofs[5]),
574 &(elementDofs[6]),&(elementDofs[7]),&(elementDofs[8]),
575 &(elementDofs[9]));
576
577 // insert each vertex if it hasn't been inserted already
578 for (int i=0; i<nVertices[elm_type]; i++)
579 if (renumber.find(elementDofs[i])==renumber.end())
580 {
581 renumber[elementDofs[i]] = number_of_real_vertices++;
582 factory.insertVertex(nodes[elementDofs[i]]);
583 }
584
585 // count elements and boundary elements
586 if (elementDim[elm_type] == dim)
587 element_count++;
588 else
589 boundary_element_count++;
590
591 }
592
593
594
595 // generic-case: This is not supposed to be used at runtime.
596 template <class E, class V, class V2>
597 void boundarysegment_insert(
598 const V&,
599 const E&,
600 const V2&
601 )
602 {
603 DUNE_THROW(Dune::IOError, "tried to create a 3D boundary segment in a non-3D Grid");
604 }
605
606 // 3d-case:
607 template <class E, class V>
608 void boundarysegment_insert(
609 const std::vector<FieldVector<double, 3> >& nodes,
610 const E& elementDofs,
611 const V& vertices
612 )
613 {
614 std::array<FieldVector<double,dimWorld>, 6> v;
615 for (int i=0; i<6; i++)
616 for (int j=0; j<dimWorld; j++)
617 v[i][j] = nodes[elementDofs[i]][j];
618
619 BoundarySegment<dim,dimWorld>* newBoundarySegment
620 = (BoundarySegment<dim,dimWorld>*) new GmshReaderQuadraticBoundarySegment< 3, 3 >( v[0], v[1], v[2],
621 v[3], v[4], v[5] );
622
623 factory.insertBoundarySegment( vertices,
624 std::shared_ptr<BoundarySegment<dim,dimWorld> >(newBoundarySegment) );
625 }
626
627
628
633 virtual void pass2HandleElement(FILE* file, const int elm_type,
634 std::map<int,unsigned int> & renumber,
635 const std::vector< GlobalVector > & nodes,
636 const int physical_entity)
637 {
638 // some data about gmsh elements
639 const int nDofs[16] = {-1, 2, 3, 4, 4, 8, 6, 5, 3, 6, -1, 10, -1, -1, -1, 1};
640 const int nVertices[16] = {-1, 2, 3, 4, 4, 8, 6, 5, 2, 3, -1, 4, -1, -1, -1, 1};
641 const int elementDim[16] = {-1, 1, 2, 2, 3, 3, 3, 3, 1, 2, -1, 3, -1, -1, -1, 0};
642
643 // test whether we support the element type
644 if ( not (elm_type > 0 && elm_type <= 15 // index in suitable range?
645 && (elementDim[elm_type] == dim || elementDim[elm_type] == (dim-1) ) ) ) // real element or boundary element?
646 {
647 skipline(file); // skip rest of line if element is unknown
648 return;
649 }
650
651 // The format string for parsing is n times '%d' in a row
652 std::string formatString = "%d";
653 for (int i=1; i<nDofs[elm_type]; i++)
654 formatString += " %d";
655 formatString += "\n";
656
657 // '10' is the largest number of dofs we may encounter in a .msh file
658 std::vector<int> elementDofs(10);
659
660 readfile(file,nDofs[elm_type], formatString.c_str(),
661 &(elementDofs[0]),&(elementDofs[1]),&(elementDofs[2]),
662 &(elementDofs[3]),&(elementDofs[4]),&(elementDofs[5]),
663 &(elementDofs[6]),&(elementDofs[7]),&(elementDofs[8]),
664 &(elementDofs[9]));
665
666 // correct differences between gmsh and Dune in the local vertex numbering
667 switch (elm_type)
668 {
669 case 3 : // 4-node quadrilateral
670 std::swap(elementDofs[2],elementDofs[3]);
671 break;
672 case 5 : // 8-node hexahedron
673 std::swap(elementDofs[2],elementDofs[3]);
674 std::swap(elementDofs[6],elementDofs[7]);
675 break;
676 case 7 : // 5-node pyramid
677 std::swap(elementDofs[2],elementDofs[3]);
678 break;
679 }
680
681 // renumber corners to account for the explicitly given vertex
682 // numbering in the file
683 std::vector<unsigned int> vertices(nVertices[elm_type]);
684
685 for (int i=0; i<nVertices[elm_type]; i++)
686 vertices[i] = renumber[elementDofs[i]];
687
688 // If it is an element, insert it as such
689 if (elementDim[elm_type] == dim) {
690
691 switch (elm_type)
692 {
693 case 1 : // 2-node line
695 break;
696 case 2 : // 3-node triangle
698 break;
699 case 3 : // 4-node quadrilateral
701 break;
702 case 4 : // 4-node tetrahedron
704 break;
705 case 5 : // 8-node hexahedron
707 break;
708 case 6 : // 6-node prism
710 break;
711 case 7 : // 5-node pyramid
713 break;
714 case 9 : // 6-node triangle
716 break;
717 case 11 : // 10-node tetrahedron
719 break;
720 }
721
722 } else {
723 // it must be a boundary segment then
724 if (insert_boundary_segments) {
725
726 switch (elm_type)
727 {
728 case 1 : // 2-node line
729 factory.insertBoundarySegment(vertices);
730 break;
731
732 case 2 : // 3-node triangle
733 factory.insertBoundarySegment(vertices);
734 break;
735
736 case 3 : // 4-node quadrilateral
737 factory.insertBoundarySegment(vertices);
738 break;
739
740 case 15 : // 1-node point
741 factory.insertBoundarySegment(vertices);
742 break;
743
744 case 8 : { // 3-node line
745 std::array<FieldVector<double,dimWorld>, 3> v;
746 for (int i=0; i<dimWorld; i++) {
747 v[0][i] = nodes[elementDofs[0]][i];
748 v[1][i] = nodes[elementDofs[2]][i]; // yes, the renumbering is intended!
749 v[2][i] = nodes[elementDofs[1]][i];
750 }
751 BoundarySegment<dim,dimWorld>* newBoundarySegment
752 = (BoundarySegment<dim,dimWorld>*) new GmshReaderQuadraticBoundarySegment< 2, dimWorld >(v[0], v[1], v[2]);
753 factory.insertBoundarySegment(vertices,
754 std::shared_ptr<BoundarySegment<dim,dimWorld> >(newBoundarySegment));
755 break;
756 }
757 case 9 : { // 6-node triangle
758 boundarysegment_insert(nodes, elementDofs, vertices);
759 break;
760 }
761 default: {
762 DUNE_THROW(Dune::IOError, "GmshReader does not support using element-type " << elm_type << " for boundary segments");
763 break;
764 }
765
766 }
767
768 }
769 }
770
771 // count elements and boundary elements
772 if (elementDim[elm_type] == dim) {
773 element_index_to_physical_entity[element_count] = physical_entity;
774 element_count++;
775 } else {
776 boundary_id_to_physical_entity[boundary_element_count] = physical_entity;
777 boundary_element_count++;
778 }
779
780 }
781
782 };
783
784 namespace Gmsh {
790 enum class ReaderOptions
791 {
792 verbose = 1,
793 insertBoundarySegments = 2,
794 readElementData = 4,
795 readBoundaryData = 8
796 };
797
799 constexpr ReaderOptions operator | (ReaderOptions a, ReaderOptions b)
800 {
801 return static_cast<ReaderOptions>(
802 static_cast<int>(a) | static_cast<int>(b)
803 );
804 }
805
807 constexpr bool operator & (ReaderOptions a, ReaderOptions b)
808 {
809 return static_cast<int>(a) & static_cast<int>(b);
810 }
811
812 } // end namespace Gmsh
813
838 template<typename GridType>
840 {
842
861 static void doRead(Dune::GridFactory<GridType> &factory,
862 const std::string &fileName,
863 std::vector<int>& boundarySegmentToPhysicalEntity,
864 std::vector<int>& elementToPhysicalEntity,
865 bool verbose, bool insertBoundarySegments)
866 {
867 // register boundary segment to boundary segment factory for possible load balancing
868 // this needs to be done on all cores since the type might not be known otherwise
869 GmshReaderQuadraticBoundarySegment< Grid::dimension, Grid::dimensionworld >::registerFactory();
870
871#ifndef NDEBUG
872 // check that this method is called on all cores
873 factory.comm().barrier();
874#endif
875
876 // create parse object and read grid on process 0
877 if (factory.comm().rank() == 0)
878 {
879 GmshReaderParser<Grid> parser(factory,verbose,insertBoundarySegments);
880 parser.read(fileName);
881
882 boundarySegmentToPhysicalEntity = std::move(parser.boundaryIdMap());
883 elementToPhysicalEntity = std::move(parser.elementIndexMap());
884 }
885 else
886 {
887 boundarySegmentToPhysicalEntity = {};
888 elementToPhysicalEntity = {};
889 }
890 }
891
893
912 template<class T>
913 static T &discarded(T &&value) { return static_cast<T&>(value); }
914
915 struct DataArg {
916 std::vector<int> *data_ = nullptr;
917 DataArg(std::vector<int> &data) : data_(&data) {}
918 DataArg(const decltype(std::ignore)&) {}
919 DataArg() = default;
920 };
921
922 struct DataFlagArg : DataArg {
923 bool flag_ = false;
924 using DataArg::DataArg;
925 DataFlagArg(bool flag) : flag_(flag) {}
926 };
927
928 public:
929 typedef GridType Grid;
930
937 static std::unique_ptr<Grid> read (const std::string& fileName, bool verbose = true, bool insertBoundarySegments=true)
938 {
939 // make a grid factory
941
942 read(factory, fileName, verbose, insertBoundarySegments);
943
944 return factory.createGrid();
945 }
946
966 static std::unique_ptr<Grid> read (const std::string& fileName,
967 std::vector<int>& boundarySegmentToPhysicalEntity,
968 std::vector<int>& elementToPhysicalEntity,
969 bool verbose = true, bool insertBoundarySegments=true)
970 {
971 // make a grid factory
973
974 doRead(
975 factory, fileName, boundarySegmentToPhysicalEntity,
976 elementToPhysicalEntity, verbose, insertBoundarySegments
977 );
978
979 return factory.createGrid();
980 }
981
983 static void read (Dune::GridFactory<Grid>& factory, const std::string& fileName,
984 bool verbose = true, bool insertBoundarySegments=true)
985 {
986 doRead(
987 factory, fileName, discarded(std::vector<int>{}),
988 discarded(std::vector<int>{}), verbose, insertBoundarySegments
989 );
990 }
991
993
1016 static void read (Dune::GridFactory<Grid> &factory,
1017 const std::string &fileName,
1018 DataFlagArg boundarySegmentData,
1019 DataArg elementData,
1020 bool verbose=true)
1021 {
1022 doRead(
1023 factory, fileName,
1024 boundarySegmentData.data_
1025 ? *boundarySegmentData.data_ : discarded(std::vector<int>{}),
1026 elementData.data_
1027 ? *elementData.data_ : discarded(std::vector<int>{}),
1028 verbose,
1029 boundarySegmentData.flag_ || boundarySegmentData.data_
1030 );
1031 }
1032
1053 static void read (Dune::GridFactory<Grid>& factory,
1054 const std::string& fileName,
1055 std::vector<int>& boundarySegmentToPhysicalEntity,
1056 std::vector<int>& elementToPhysicalEntity,
1057 bool verbose, bool insertBoundarySegments)
1058 {
1059 doRead(
1060 factory, fileName, boundarySegmentToPhysicalEntity,
1061 elementToPhysicalEntity, verbose, insertBoundarySegments
1062 );
1063 }
1064
1066 //\{
1067
1069
1070 static constexpr Opts defaultOpts =
1071 Opts::verbose | Opts::insertBoundarySegments | Opts::readElementData | Opts::readBoundaryData;
1072
1074
1097 GmshReader(const std::string& fileName,
1098 Gmsh::ReaderOptions options = defaultOpts)
1099 {
1100 gridFactory_ = std::make_unique<Dune::GridFactory<Grid>>();
1101 readGridFile(fileName, *gridFactory_, options);
1102 }
1103
1111 GmshReader(const std::string& fileName, GridFactory<Grid>& factory,
1112 Gmsh::ReaderOptions options = defaultOpts)
1113 {
1114 readGridFile(fileName, factory, options);
1115 }
1116
1118 const std::vector<int>& elementData () const
1119 {
1120 checkElementData();
1121 return elementIndexToGmshPhysicalEntity_;
1122 }
1123
1125 const std::vector<int>& boundaryData () const
1126 {
1127 checkBoundaryData();
1128 return boundarySegmentIndexToGmshPhysicalEntity_;
1129 }
1130
1135 bool hasElementData () const
1136 { return hasElementData_ && !extractedElementData_; }
1137
1142 bool hasBoundaryData () const
1143 { return hasBoundaryData_ && !extractedBoundaryData_; }
1144
1146 std::vector<int> extractElementData ()
1147 {
1148 checkElementData();
1149 extractedElementData_ = true;
1150 return std::move(elementIndexToGmshPhysicalEntity_);
1151 }
1152
1154 std::vector<int> extractBoundaryData ()
1155 {
1156 checkBoundaryData();
1157 extractedBoundaryData_ = true;
1158 return std::move(boundarySegmentIndexToGmshPhysicalEntity_);
1159 }
1160
1162 std::unique_ptr<Grid> createGrid ()
1163 {
1164 if (!gridFactory_)
1166 "This GmshReader has been constructed with a Dune::GridFactory. "
1167 << "This grid factory has been filled with all information to create a grid. "
1168 << "Please use this factory to create the grid by calling factory.createGrid(). "
1169 << "Alternatively use the constructor without passing the factory in combination with this member function."
1170 );
1171
1172 return gridFactory_->createGrid();
1173 }
1174
1175 //\}
1176
1177 private:
1178 void checkElementData () const
1179 {
1180 if (!hasElementData_)
1182 "This GmshReader has been constructed without the option 'readElementData'. "
1183 << "Please enable reading element data by passing the option 'Gmsh::ReaderOpts::readElementData' "
1184 << "to the constructor of this class."
1185 );
1186
1187 if (extractedElementData_)
1189 "The element data has already been extracted from this GmshReader "
1190 << "via a function call to reader.extractElementData(). Use the extracted data or "
1191 << "read the grid data from file again by constructing a new reader."
1192 );
1193 }
1194
1195 void checkBoundaryData () const
1196 {
1197 if (!hasBoundaryData_)
1199 "This GmshReader has been constructed without the option 'readBoundaryData'. "
1200 << "Please enable reading boundary data by passing the option 'Gmsh::ReaderOpts::readBoundaryData' "
1201 << "to the constructor of this class."
1202 );
1203
1204 if (extractedBoundaryData_)
1206 "The boundary data has already been extracted from this GmshReader "
1207 << "via a function call to reader.extractBoundaryData(). Use the extracted data or "
1208 << "read the grid data from file again by constructing a new reader."
1209 );
1210 }
1211
1212 void readGridFile (const std::string& fileName, GridFactory<Grid>& factory, Gmsh::ReaderOptions options)
1213 {
1214 const bool verbose = options & Opts::verbose;
1215 const bool insertBoundarySegments = options & Opts::insertBoundarySegments;
1216 const bool readBoundaryData = options & Opts::readBoundaryData;
1217 const bool readElementData = options & Opts::readElementData;
1218
1219 doRead(
1220 factory, fileName, boundarySegmentIndexToGmshPhysicalEntity_,
1221 elementIndexToGmshPhysicalEntity_, verbose,
1222 readBoundaryData || insertBoundarySegments
1223 );
1224
1225 // clear unwanted data
1226 if (!readBoundaryData)
1227 boundarySegmentIndexToGmshPhysicalEntity_ = std::vector<int>{};
1228 if (!readElementData)
1229 elementIndexToGmshPhysicalEntity_ = std::vector<int>{};
1230
1231 hasElementData_ = readElementData;
1232 hasBoundaryData_ = readBoundaryData;
1233 }
1234
1235 std::unique_ptr<Dune::GridFactory<Grid>> gridFactory_;
1236
1237 std::vector<int> elementIndexToGmshPhysicalEntity_;
1238 std::vector<int> boundarySegmentIndexToGmshPhysicalEntity_;
1239
1240 bool hasElementData_;
1241 bool hasBoundaryData_;
1242
1243 // for better error messages, we keep track of these separately
1244 bool extractedElementData_ = false;
1245 bool extractedBoundaryData_ = false;
1246 };
1247
1250} // namespace Dune
1251
1252#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
FieldTraits< value_type >::real_type two_norm() const
two norm sqrt(sum over squared values of entries)
Definition: densevector.hh:641
derived_type & axpy(const field_type &a, const DenseVector< Other > &x)
vector space axpy operation ( *this += a x )
Definition: densevector.hh:575
vector space out of a tensor product of fields.
Definition: fvector.hh:91
dimension independent parts for GmshReaderParser
Definition: gmshreader.hh:302
std::vector< std::string > & physicalEntityNames()
Returns the names of the gmsh physical entities (0-based index)
Definition: gmshreader.hh:373
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:545
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
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:633
Read Gmsh mesh file.
Definition: gmshreader.hh:840
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:966
const std::vector< int > & elementData() const
Access element data (maps element index to Gmsh physical entity)
Definition: gmshreader.hh:1118
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:1016
static std::unique_ptr< Grid > read(const std::string &fileName, bool verbose=true, bool insertBoundarySegments=true)
Definition: gmshreader.hh:937
static void read(Dune::GridFactory< Grid > &factory, const std::string &fileName, bool verbose=true, bool insertBoundarySegments=true)
Definition: gmshreader.hh:983
std::unique_ptr< Grid > createGrid()
Create the grid.
Definition: gmshreader.hh:1162
std::vector< int > extractBoundaryData()
Erase boundary data from reader and return the data.
Definition: gmshreader.hh:1154
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:1053
bool hasElementData() const
If element data is available.
Definition: gmshreader.hh:1135
bool hasBoundaryData() const
If boundary data is available.
Definition: gmshreader.hh:1142
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:1111
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:1097
std::vector< int > extractElementData()
Erase element data from reader and return the data.
Definition: gmshreader.hh:1146
const std::vector< int > & boundaryData() const
Access boundary data (maps boundary segment index to Gmsh physical entity)
Definition: gmshreader.hh:1125
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 insertElement(const GeometryType &type, const std::vector< unsigned int > &vertices)
Insert an element into the coarse grid.
Definition: gridfactory.hh:307
virtual void insertVertex(const FieldVector< ctype, dimworld > &pos)
Insert a vertex into the coarse grid.
Definition: gridfactory.hh:296
virtual void insertBoundarySegment(const std::vector< unsigned int > &vertices)
insert a boundary segment
Definition: gridfactory.hh:325
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:791
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.111.3 (Nov 12, 23:30, 2024)