Dune Core Modules (2.8.0)

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