4 #ifndef DUNE_GMSHREADER_HH
5 #define DUNE_GMSHREADER_HH
18 #include <dune/common/exceptions.hh>
19 #include <dune/common/fvector.hh>
21 #include <dune/geometry/type.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;
74 typedef Dune::FieldVector< double, dimWorld >
GlobalVector;
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()
100 virtual GlobalVector operator() (
const Dune::FieldVector<double,1> &local )
const
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 );
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");
169 class GmshReaderQuadraticBoundarySegment< 3, 3 >
172 typedef GmshReaderQuadraticBoundarySegment< 3, 3 > ThisType;
175 GmshReaderQuadraticBoundarySegment (Dune::FieldVector<double,3> p0_, Dune::FieldVector<double,3> p1_,
176 Dune::FieldVector<double,3> p2_, Dune::FieldVector<double,3> p3_,
177 Dune::FieldVector<double,3> p4_, Dune::FieldVector<double,3> p5_)
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()
203 virtual Dune::FieldVector<double,3> operator() (
const Dune::FieldVector<double,2>& local)
const
205 Dune::FieldVector<double,3> y;
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 );
235 Dune::FieldVector<double,3> d1,d2;
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");
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");
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");
265 double phi0 (
const Dune::FieldVector<double,2>& local)
const
267 return (alpha*beta-beta*local[0]-alpha*local[1])*(1-local[0]-local[1])/(alpha*beta);
269 double phi3 (
const Dune::FieldVector<double,2>& local)
const
271 return local[0]*(1-local[0]-local[1])/(alpha*(1-alpha));
273 double phi1 (
const Dune::FieldVector<double,2>& local)
const
275 return local[0]*(gamma*local[0]-(sqrt2-gamma-sqrt2*alpha)*local[1]-alpha*gamma)/(gamma*(1-alpha));
277 double phi5 (
const Dune::FieldVector<double,2>& local)
const
279 return local[1]*(1-local[0]-local[1])/(beta*(1-beta));
281 double phi4 (
const Dune::FieldVector<double,2>& local)
const
283 return local[0]*local[1]/((1-gamma/sqrt2)*gamma/sqrt2);
285 double phi2 (
const Dune::FieldVector<double,2>& local)
const
287 return local[1]*(beta*(1-gamma/sqrt2)-local[0]*(beta-gamma/sqrt2)-local[1]*(1-gamma/sqrt2))/((1-gamma/sqrt2)*(beta-1));
290 Dune::FieldVector<double,3> p0,p1,p2,p3,p4,p5;
291 double alpha,beta,gamma,sqrt2;
297 template<
typename Gr
idType>
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);
335 DUNE_THROW(Dune::IOError,
"Error parsing " <<
fileName <<
" "
337 <<
": Expected '" << format <<
"', only read " << c <<
" entries instead of " << cnt <<
".");
345 c = std::fgetc(file);
346 }
while(c !=
'\n' && c != EOF);
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");
372 DUNE_THROW(Dune::IOError,
"Could not open " <<
fileName);
384 double version_number;
385 int file_type, data_size;
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;
395 if (strcmp(
buf,
"$EndMeshFormat")!=0)
396 DUNE_THROW(Dune::IOError,
"expected $EndMeshFormat");
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;
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) {
418 DUNE_THROW(Dune::IOError,
419 "Only dense sequences of node indices are currently supported (node index "
420 <<
id <<
" is invalid).");
425 nodes[
id ][ j ] = x[ j ];
428 if (strcmp(
buf,
"$EndNodes")!=0)
429 DUNE_THROW(Dune::IOError,
"expected $EndNodes");
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;
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++)
469 if (strcmp(
buf,
"$EndElements")!=0)
470 DUNE_THROW(Dune::IOError,
"expected $EndElements");
478 fseeko(file, section_element_offset, SEEK_SET);
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++)
491 if (k==1) physical_entity = blub;
496 if (strcmp(
buf,
"$EndElements")!=0)
497 DUNE_THROW(Dune::IOError,
"expected $EndElements");
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) ) ) )
525 std::string formatString =
"%d";
526 for (
int i=1; i<nDofs[elm_type]; i++)
527 formatString +=
" %d";
528 formatString +=
"\n";
531 std::vector<int> elementDofs(10);
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]),
540 for (
int i=0; i<nVertices[elm_type]; i++)
544 factory.insertVertex(nodes[elementDofs[i]]);
548 if (elementDim[elm_type] ==
dim)
558 template <
class E,
class V,
class V2>
565 DUNE_THROW(Dune::IOError,
"tried to create a 3D boundary segment in a non-3D Grid");
569 template <
class E,
class V>
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++)
579 v[i][j] = nodes[elementDofs[i]][j];
585 factory.insertBoundarySegment( vertices,
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) ) ) )
614 std::string formatString =
"%d";
615 for (
int i=1; i<nDofs[elm_type]; i++)
616 formatString +=
" %d";
617 formatString +=
"\n";
620 std::vector<int> elementDofs(10);
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]),
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) {
691 factory.insertBoundarySegment(vertices);
695 factory.insertBoundarySegment(vertices);
699 factory.insertBoundarySegment(vertices);
703 factory.insertBoundarySegment(vertices);
707 std::array<FieldVector<double,dimWorld>, 3> v;
709 v[0][i] = nodes[elementDofs[0]][i];
710 v[1][i] = nodes[elementDofs[2]][i];
711 v[2][i] = nodes[elementDofs[1]][i];
715 factory.insertBoundarySegment(vertices,
724 DUNE_THROW(Dune::IOError,
"GmshReader does not support using element-type " << elm_type <<
" for boundary segments");
734 if (elementDim[elm_type] ==
dim) {
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();
805 factory.
comm().barrier();
809 if (factory.
comm().rank() == 0)
812 parser.
read(fileName);
814 boundarySegmentToPhysicalEntity = std::move(parser.
boundaryIdMap());
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) {}
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.
Include standard header files.
Definition: agrid.hh:58
static const int dimWorld
Definition: misc.hh:44
ALBERTA REAL_D GlobalVector
Definition: misc.hh:48
int renumber(const Dune::GeometryType &t, int i)
renumber VTK <-> Dune
Definition: common.hh:184
@ line
Definition: common.hh:132
@ pyramid
Definition: common.hh:139
@ quadrilateral
Definition: common.hh:135
@ tetrahedron
Definition: common.hh:136
@ prism
Definition: common.hh:138
@ hexahedron
Definition: common.hh:137
@ triangle
Definition: common.hh:133
Base class for classes implementing geometries of boundary segments.
Definition: boundarysegment.hh:92
Communication comm() const
Return the Communication used by the grid factory.
Definition: common/gridfactory.hh:295
Provide a generic factory class for unstructured grids.
Definition: common/gridfactory.hh:312
virtual std::unique_ptr< GridType > createGrid()
Finalize grid creation and hand over the grid.
Definition: common/gridfactory.hh:370
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
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
static const int dimWorld
Definition: gmshreader.hh:317
Dune::GridFactory< GridType > & factory
Definition: gmshreader.hh:302
std::vector< int > & elementIndexMap()
Definition: gmshreader.hh:359
unsigned int number_of_real_vertices
Definition: gmshreader.hh:305
void boundarysegment_insert(const std::vector< FieldVector< double, 3 > > &nodes, const E &elementDofs, const V &vertices)
Definition: gmshreader.hh:570
GmshReaderParser(Dune::GridFactory< GridType > &_factory, bool v, bool i)
Definition: gmshreader.hh:351
int element_count
Definition: gmshreader.hh:307
void read(const std::string &f)
Definition: gmshreader.hh:364
void skipline(FILE *file)
Definition: gmshreader.hh:341
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:327
std::vector< int > element_index_to_physical_entity
Definition: gmshreader.hh:313
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
std::vector< int > & boundaryIdMap()
Definition: gmshreader.hh:354
static const int dim
Definition: gmshreader.hh:316
FieldVector< double, dimWorld > GlobalVector
Definition: gmshreader.hh:318
std::string fileName
Definition: gmshreader.hh:310
int boundary_element_count
Definition: gmshreader.hh:306
void boundarysegment_insert(const V &, const E &, const V2 &)
Definition: gmshreader.hh:559
bool verbose
Definition: gmshreader.hh:303
std::vector< int > boundary_id_to_physical_entity
Definition: gmshreader.hh:312
char buf[512]
Definition: gmshreader.hh:309
bool insert_boundary_segments
Definition: gmshreader.hh:304
Read Gmsh mesh file.
Definition: gmshreader.hh:772
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 void read(Dune::GridFactory< Grid > &factory, const std::string &fileName, bool verbose=true, bool insertBoundarySegments=true)
Definition: gmshreader.hh:913
GridType Grid
Definition: gmshreader.hh:861
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
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 std::unique_ptr< Grid > read(const std::string &fileName, bool verbose=true, bool insertBoundarySegments=true)
Definition: gmshreader.hh:869
Provide a generic factory class for unstructured grids.