4#ifndef DUNE_GMSHREADER_HH
5#define DUNE_GMSHREADER_HH
48 template<
int dimension,
int dimWorld = dimension >
49 class GmshReaderQuadraticBoundarySegment
53 static void registerFactory() {}
68 template<
int dimWorld >
69 struct GmshReaderQuadraticBoundarySegment< 2, dimWorld >
72 typedef GmshReaderQuadraticBoundarySegment< 2, dimWorld > ThisType;
76 GmshReaderQuadraticBoundarySegment (
const GlobalVector &p0_,
const GlobalVector &p1_,
const GlobalVector &p2_)
77 : p0(p0_), p1(p1_), p2(p2_)
82 GmshReaderQuadraticBoundarySegment( ObjectStreamType& in )
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 );
92 static void registerFactory()
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);
110 void backup( ObjectStreamType& out )
const
113 out.write( (
const char *) &key(),
sizeof(
int ) );
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 );
124 GlobalVector d1 = p1;
126 GlobalVector d2 = p2;
129 alpha=d1.two_norm()/(d1.two_norm()+d2.two_norm());
130 if (alpha<1E-6 || alpha>1-1E-6)
140 GlobalVector p0,p1,p2;
169 class GmshReaderQuadraticBoundarySegment< 3, 3 >
172 typedef GmshReaderQuadraticBoundarySegment< 3, 3 > ThisType;
178 : p0(p0_), p1(p1_), p2(p2_), p3(p3_), p4(p4_), p5(p5_)
183 GmshReaderQuadraticBoundarySegment( ObjectStreamType& in )
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 );
195 static void registerFactory()
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);
216 void backup( ObjectStreamType& out )
const
219 out.write( (
const char*) &key(),
sizeof(
int ) );
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 );
240 if (alpha<1E-6 || alpha>1-1E-6)
246 if (beta<1E-6 || beta>1-1E-6)
252 if (gamma<1E-6 || gamma>1-1E-6)
267 return (alpha*beta-beta*local[0]-alpha*local[1])*(1-local[0]-local[1])/(alpha*beta);
271 return local[0]*(1-local[0]-local[1])/(alpha*(1-alpha));
275 return local[0]*(gamma*local[0]-(sqrt2-gamma-sqrt2*alpha)*local[1]-alpha*gamma)/(gamma*(1-alpha));
279 return local[1]*(1-local[0]-local[1])/(beta*(1-beta));
283 return local[0]*local[1]/((1-gamma/sqrt2)*gamma/sqrt2);
287 return local[1]*(beta*(1-gamma/sqrt2)-local[0]*(beta-gamma/sqrt2)-local[1]*(1-gamma/sqrt2))/((1-gamma/sqrt2)*(beta-1));
291 double alpha,beta,gamma,sqrt2;
297 template<
typename Gr
idType>
304 bool insert_boundary_segments;
305 unsigned int number_of_real_vertices;
306 int boundary_element_count;
310 std::string fileName;
312 std::vector<int> boundary_id_to_physical_entity;
313 std::vector<int> element_index_to_physical_entity;
316 static const int dim = GridType::dimension;
317 static const int dimWorld = GridType::dimensionworld;
318 static_assert( (dimWorld <= 3),
"GmshReader requires dimWorld <= 3." );
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)
332 off_t pos = ftello(file);
333 int c = fscanf(file, format, t1, t2, t3, t4, t5, t6, t7, t8, t9, t10);
337 <<
": Expected '" << format <<
"', only read " << c <<
" entries instead of " << cnt <<
".");
341 void skipline(FILE * file)
345 c = std::fgetc(file);
346 }
while(c !=
'\n' && c != EOF);
352 factory(_factory), verbose(v), insert_boundary_segments(i) {}
354 std::vector<int> & boundaryIdMap()
356 return boundary_id_to_physical_entity;
359 std::vector<int> & elementIndexMap()
361 return element_index_to_physical_entity;
364 void read (
const std::string& f)
366 if (verbose) std::cout <<
"Reading " << dim <<
"d Gmsh grid..." << std::endl;
370 FILE* file = fopen(fileName.c_str(),
"rb");
379 number_of_real_vertices = 0;
380 boundary_element_count = 0;
384 double version_number;
385 int file_type, data_size;
387 readfile(file,1,
"%s\n",buf);
388 if (strcmp(buf,
"$MeshFormat")!=0)
390 readfile(file,3,
"%lg %d %d\n",&version_number,&file_type,&data_size);
391 if( (version_number < 2.0) || (version_number > 2.2) )
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)
401 readfile(file,1,
"%s\n",buf);
402 if (strcmp(buf,
"$Nodes")!=0)
404 readfile(file,1,
"%d\n",&number_of_nodes);
405 if (verbose) std::cout <<
"file contains " << number_of_nodes <<
" nodes" << std::endl;
409 std::vector< GlobalVector > nodes( number_of_nodes+1 );
413 for(
int i = 1; i <= number_of_nodes; ++i )
415 readfile(file,4,
"%d %lg %lg %lg\n", &
id, &x[ 0 ], &x[ 1 ], &x[ 2 ] );
417 if (
id > number_of_nodes) {
419 "Only dense sequences of node indices are currently supported (node index "
420 <<
id <<
" is invalid).");
424 for(
int j = 0; j < dimWorld; ++j )
425 nodes[
id ][ j ] = x[ j ];
427 readfile(file,1,
"%s\n",buf);
428 if (strcmp(buf,
"$EndNodes")!=0)
433 readfile(file,1,
"%s\n",buf);
434 if (strcmp(buf,
"$Elements")!=0)
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;
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++)
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++)
454 readfile(file,1,
"%d ",&blub);
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)
471 boundary_id_to_physical_entity.resize(boundary_element_count);
472 element_index_to_physical_entity.resize(element_count);
478 fseeko(file, section_element_offset, SEEK_SET);
479 boundary_element_count = 0;
481 for (
int i=1; i<=number_of_elements; i++)
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;
487 for (
int k=1; k<=number_of_tags; k++)
490 readfile(file,1,
"%d ",&blub);
491 if (k==1) physical_entity = blub;
495 readfile(file,1,
"%s\n",buf);
496 if (strcmp(buf,
"$EndElements")!=0)
508 std::map<int,unsigned int> & renumber,
509 const std::vector< GlobalVector > & nodes)
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};
517 if ( not (elm_type > 0 && elm_type <= 15
518 && (elementDim[elm_type] == dim || elementDim[elm_type] == (dim-1) ) ) )
526 for (
int i=1; i<nDofs[elm_type]; i++)
531 std::vector<int> elementDofs(10);
534 &(elementDofs[0]),&(elementDofs[1]),&(elementDofs[2]),
535 &(elementDofs[3]),&(elementDofs[4]),&(elementDofs[5]),
536 &(elementDofs[6]),&(elementDofs[7]),&(elementDofs[8]),
540 for (
int i=0; i<nVertices[elm_type]; i++)
541 if (renumber.find(elementDofs[i])==renumber.end())
543 renumber[elementDofs[i]] = number_of_real_vertices++;
548 if (elementDim[elm_type] == dim)
551 boundary_element_count++;
558 template <
class E,
class V,
class V2>
559 void boundarysegment_insert(
569 template <
class E,
class V>
570 void boundarysegment_insert(
571 const std::vector<FieldVector<double, 3> >& nodes,
572 const E& elementDofs,
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];
581 BoundarySegment<dim,dimWorld>* newBoundarySegment
582 = (BoundarySegment<dim,dimWorld>*)
new GmshReaderQuadraticBoundarySegment< 3, 3 >( v[0], v[1], v[2],
586 std::shared_ptr<BoundarySegment<dim,dimWorld> >(newBoundarySegment) );
596 std::map<int,unsigned int> & renumber,
597 const std::vector< GlobalVector > & nodes,
598 const int physical_entity)
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};
606 if ( not (elm_type > 0 && elm_type <= 15
607 && (elementDim[elm_type] == dim || elementDim[elm_type] == (dim-1) ) ) )
615 for (
int i=1; i<nDofs[elm_type]; i++)
620 std::vector<int> elementDofs(10);
623 &(elementDofs[0]),&(elementDofs[1]),&(elementDofs[2]),
624 &(elementDofs[3]),&(elementDofs[4]),&(elementDofs[5]),
625 &(elementDofs[6]),&(elementDofs[7]),&(elementDofs[8]),
632 std::swap(elementDofs[2],elementDofs[3]);
635 std::swap(elementDofs[2],elementDofs[3]);
636 std::swap(elementDofs[6],elementDofs[7]);
639 std::swap(elementDofs[2],elementDofs[3]);
645 std::vector<unsigned int> vertices(nVertices[elm_type]);
647 for (
int i=0; i<nVertices[elm_type]; i++)
648 vertices[i] = renumber[elementDofs[i]];
651 if (elementDim[elm_type] == dim) {
686 if (insert_boundary_segments) {
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];
711 v[2][i] = nodes[elementDofs[1]][i];
720 boundarysegment_insert(nodes, elementDofs, vertices);
724 DUNE_THROW(
Dune::IOError,
"GmshReader does not support using element-type " << elm_type <<
" for boundary segments");
734 if (elementDim[elm_type] == dim) {
735 element_index_to_physical_entity[element_count] = physical_entity;
738 boundary_id_to_physical_entity[boundary_element_count] = physical_entity;
739 boundary_element_count++;
770 template<
typename Gr
idType>
794 const std::string &fileName,
795 std::vector<int>& boundarySegmentToPhysicalEntity,
796 std::vector<int>& elementToPhysicalEntity,
797 bool verbose,
bool insertBoundarySegments)
801 GmshReaderQuadraticBoundarySegment< Grid::dimension, Grid::dimensionworld >::registerFactory();
812 parser.read(fileName);
814 boundarySegmentToPhysicalEntity = std::move(parser.boundaryIdMap());
815 elementToPhysicalEntity = std::move(parser.elementIndexMap());
819 boundarySegmentToPhysicalEntity = {};
820 elementToPhysicalEntity = {};
845 static T &discarded(T &&value) {
return value; }
848 std::vector<int> *data_ =
nullptr;
849 DataArg(std::vector<int> &data) : data_(&data) {}
850 DataArg(
const decltype(std::ignore)&) {}
854 struct DataFlagArg : DataArg {
856 using DataArg::DataArg;
857 DataFlagArg(
bool flag) : flag_(flag) {}
861 typedef GridType Grid;
869 static std::unique_ptr<Grid>
read (
const std::string& fileName,
bool verbose =
true,
bool insertBoundarySegments=
true)
874 read(factory, fileName, verbose, insertBoundarySegments);
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)
906 do_read(factory, fileName, boundarySegmentToPhysicalEntity,
907 elementToPhysicalEntity, verbose, insertBoundarySegments);
914 bool verbose =
true,
bool insertBoundarySegments=
true)
916 do_read(factory, fileName, discarded(std::vector<int>{}),
917 discarded(std::vector<int>{}), verbose, insertBoundarySegments);
945 const std::string &fileName,
946 DataFlagArg boundarySegmentData,
950 do_read(factory, fileName,
951 boundarySegmentData.data_
952 ? *boundarySegmentData.data_ : discarded(std::vector<int>{}),
954 ? *elementData.data_ : discarded(std::vector<int>{}),
956 boundarySegmentData.flag_ || boundarySegmentData.data_);
980 const std::string& fileName,
981 std::vector<int>& boundarySegmentToPhysicalEntity,
982 std::vector<int>& elementToPhysicalEntity,
983 bool verbose,
bool insertBoundarySegments)
985 do_read(factory, fileName, boundarySegmentToPhysicalEntity,
986 elementToPhysicalEntity, verbose, insertBoundarySegments);
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.