4 #ifndef DUNE_GMSHREADER_HH 5 #define DUNE_GMSHREADER_HH 17 #include <dune/common/exceptions.hh> 18 #include <dune/common/fvector.hh> 20 #include <dune/geometry/type.hh> 33 struct GmshReaderOptions
47 template<
int dimension,
int dimWorld = dimension >
48 class GmshReaderQuadraticBoundarySegment
63 template<
int dimWorld >
64 struct GmshReaderQuadraticBoundarySegment< 2, dimWorld >
67 typedef Dune::FieldVector< double, dimWorld >
GlobalVector;
69 GmshReaderQuadraticBoundarySegment (
const GlobalVector &p0_,
const GlobalVector &p1_,
const GlobalVector &p2_)
70 : p0(p0_), p1(p1_), p2(p2_)
77 alpha=d1.two_norm()/(d1.two_norm()+d2.two_norm());
78 if (alpha<1E-6 || alpha>1-1E-6)
79 DUNE_THROW(Dune::IOError,
"ration in quadratic boundary segment bad");
82 virtual GlobalVector operator() (
const Dune::FieldVector<double,1> &local )
const 86 y.axpy((local[0]-alpha)*(local[0]-1.0)/alpha,p0);
87 y.axpy(local[0]*(local[0]-1.0)/(alpha*(alpha-1.0)),p1);
88 y.axpy(local[0]*(local[0]-alpha)/(1.0-alpha),p2);
93 GlobalVector p0,p1,p2;
122 class GmshReaderQuadraticBoundarySegment< 3, 3 >
126 GmshReaderQuadraticBoundarySegment (Dune::FieldVector<double,3> p0_, Dune::FieldVector<double,3> p1_,
127 Dune::FieldVector<double,3> p2_, Dune::FieldVector<double,3> p3_,
128 Dune::FieldVector<double,3> p4_, Dune::FieldVector<double,3> p5_)
129 : p0(p0_), p1(p1_), p2(p2_), p3(p3_), p4(p4_), p5(p5_)
132 Dune::FieldVector<double,3> d1,d2;
136 alpha=d1.two_norm()/(d1.two_norm()+d2.two_norm());
137 if (alpha<1E-6 || alpha>1-1E-6)
138 DUNE_THROW(Dune::IOError,
"alpha in quadratic boundary segment bad");
142 beta=d1.two_norm()/(d1.two_norm()+d2.two_norm());
143 if (beta<1E-6 || beta>1-1E-6)
144 DUNE_THROW(Dune::IOError,
"beta in quadratic boundary segment bad");
148 gamma=sqrt2*(d1.two_norm()/(d1.two_norm()+d2.two_norm()));
149 if (gamma<1E-6 || gamma>1-1E-6)
150 DUNE_THROW(Dune::IOError,
"gamma in quadratic boundary segment bad");
153 virtual Dune::FieldVector<double,3> operator() (
const Dune::FieldVector<double,2>& local)
const 155 Dune::FieldVector<double,3> y;
157 y.axpy(phi0(local),p0);
158 y.axpy(phi1(local),p1);
159 y.axpy(phi2(local),p2);
160 y.axpy(phi3(local),p3);
161 y.axpy(phi4(local),p4);
162 y.axpy(phi5(local),p5);
170 double phi0 (
const Dune::FieldVector<double,2>& local)
const 172 return (alpha*beta-beta*local[0]-alpha*local[1])*(1-local[0]-local[1])/(alpha*beta);
174 double phi3 (
const Dune::FieldVector<double,2>& local)
const 176 return local[0]*(1-local[0]-local[1])/(alpha*(1-alpha));
178 double phi1 (
const Dune::FieldVector<double,2>& local)
const 180 return local[0]*(gamma*local[0]-(sqrt2-gamma-sqrt2*alpha)*local[1]-alpha*gamma)/(gamma*(1-alpha));
182 double phi5 (
const Dune::FieldVector<double,2>& local)
const 184 return local[1]*(1-local[0]-local[1])/(beta*(1-beta));
186 double phi4 (
const Dune::FieldVector<double,2>& local)
const 188 return local[0]*local[1]/((1-gamma/sqrt2)*gamma/sqrt2);
190 double phi2 (
const Dune::FieldVector<double,2>& local)
const 192 return local[1]*(beta*(1-gamma/sqrt2)-local[0]*(beta-gamma/sqrt2)-local[1]*(1-gamma/sqrt2))/((1-gamma/sqrt2)*(beta-1));
195 Dune::FieldVector<double,3> p0,p1,p2,p3,p4,p5;
196 double alpha,beta,gamma,sqrt2;
202 template<
typename Gr
idType>
221 static const int dim = GridType::dimension;
222 static const int dimWorld = GridType::dimensionworld;
223 static_assert( (dimWorld <= 3),
"GmshReader requires dimWorld <= 3." );
232 void readfile(FILE * file,
int cnt,
const char * format,
233 void* t1,
void* t2 = 0,
void* t3 = 0,
void* t4 = 0,
234 void* t5 = 0,
void* t6 = 0,
void* t7 = 0,
void* t8 = 0,
235 void* t9 = 0,
void* t10 = 0)
237 off_t pos = ftello(file);
238 int c = fscanf(file, format, t1, t2, t3, t4, t5, t6, t7, t8, t9, t10);
240 DUNE_THROW(Dune::IOError,
"Error parsing " << fileName <<
" " 242 <<
": Expected '" << format <<
"', only read " << c <<
" entries instead of " << cnt <<
".");
250 c = std::fgetc(file);
251 }
while(c !=
'\n' && c != EOF);
257 factory(_factory), verbose(v), insert_boundary_segments(i) {}
261 return boundary_id_to_physical_entity;
266 return element_index_to_physical_entity;
269 void read (
const std::string& f)
271 if (verbose) std::cout <<
"Reading " << dim <<
"d Gmsh grid..." << std::endl;
275 FILE* file = fopen(fileName.c_str(),
"r");
277 DUNE_THROW(Dune::IOError,
"Could not open " << fileName);
284 number_of_real_vertices = 0;
285 boundary_element_count = 0;
289 double version_number;
290 int file_type, data_size;
292 readfile(file,1,
"%s\n",buf);
293 if (strcmp(buf,
"$MeshFormat")!=0)
294 DUNE_THROW(Dune::IOError,
"expected $MeshFormat in first line");
295 readfile(file,3,
"%lg %d %d\n",&version_number,&file_type,&data_size);
296 if( (version_number < 2.0) || (version_number > 2.2) )
297 DUNE_THROW(Dune::IOError,
"can only read Gmsh version 2 files");
298 if (verbose) std::cout <<
"version " << version_number <<
" Gmsh file detected" << std::endl;
299 readfile(file,1,
"%s\n",buf);
300 if (strcmp(buf,
"$EndMeshFormat")!=0)
301 DUNE_THROW(Dune::IOError,
"expected $EndMeshFormat");
306 readfile(file,1,
"%s\n",buf);
307 if (strcmp(buf,
"$Nodes")!=0)
308 DUNE_THROW(Dune::IOError,
"expected $Nodes");
309 readfile(file,1,
"%d\n",&number_of_nodes);
310 if (verbose) std::cout <<
"file contains " << number_of_nodes <<
" nodes" << std::endl;
314 std::vector< GlobalVector > nodes( number_of_nodes+1 );
318 for(
int i = 1; i <= number_of_nodes; ++i )
320 readfile(file,4,
"%d %lg %lg %lg\n", &
id, &x[ 0 ], &x[ 1 ], &x[ 2 ] );
322 if (
id > number_of_nodes) {
323 DUNE_THROW(Dune::IOError,
324 "Only dense sequences of node indices are currently supported (node index " 325 <<
id <<
" is invalid).");
330 nodes[
id ][ j ] = x[ j ];
332 readfile(file,1,
"%s\n",buf);
333 if (strcmp(buf,
"$EndNodes")!=0)
334 DUNE_THROW(Dune::IOError,
"expected $EndNodes");
338 readfile(file,1,
"%s\n",buf);
339 if (strcmp(buf,
"$Elements")!=0)
340 DUNE_THROW(Dune::IOError,
"expected $Elements");
341 int number_of_elements;
342 readfile(file,1,
"%d\n",&number_of_elements);
343 if (verbose) std::cout <<
"file contains " << number_of_elements <<
" elements" << std::endl;
350 long section_element_offset = ftell(file);
351 std::map<int,unsigned int>
renumber;
352 for (
int i=1; i<=number_of_elements; i++)
354 int id, elm_type, number_of_tags;
355 readfile(file,3,
"%d %d %d ",&
id,&elm_type,&number_of_tags);
356 for (
int k=1; k<=number_of_tags; k++)
359 readfile(file,1,
"%d ",&blub);
368 pass1HandleElement(file, elm_type, renumber, nodes);
370 if (verbose) std::cout <<
"number of real vertices = " << number_of_real_vertices << std::endl;
371 if (verbose) std::cout <<
"number of boundary elements = " << boundary_element_count << std::endl;
372 if (verbose) std::cout <<
"number of elements = " << element_count << std::endl;
373 readfile(file,1,
"%s\n",buf);
374 if (strcmp(buf,
"$EndElements")!=0)
375 DUNE_THROW(Dune::IOError,
"expected $EndElements");
376 boundary_id_to_physical_entity.resize(boundary_element_count);
377 element_index_to_physical_entity.resize(element_count);
383 fseek(file, section_element_offset, SEEK_SET);
384 boundary_element_count = 0;
386 for (
int i=1; i<=number_of_elements; i++)
388 int id, elm_type, number_of_tags;
389 readfile(file,3,
"%d %d %d ",&
id,&elm_type,&number_of_tags);
390 int physical_entity = -1;
392 for (
int k=1; k<=number_of_tags; k++)
395 readfile(file,1,
"%d ",&blub);
396 if (k==1) physical_entity = blub;
398 pass2HandleElement(file, elm_type, renumber, nodes, physical_entity);
400 readfile(file,1,
"%s\n",buf);
401 if (strcmp(buf,
"$EndElements")!=0)
402 DUNE_THROW(Dune::IOError,
"expected $EndElements");
413 std::map<int,unsigned int> &
renumber,
414 const std::vector< GlobalVector > & nodes)
417 const int nDofs[12] = {-1, 2, 3, 4, 4, 8, 6, 5, 3, 6, -1, 10};
418 const int nVertices[12] = {-1, 2, 3, 4, 4, 8, 6, 5, 2, 3, -1, 4};
419 const int elementDim[12] = {-1, 1, 2, 2, 3, 3, 3, 3, 1, 2, -1, 3};
422 if ( not (elm_type >= 0 && elm_type < 12
423 && (elementDim[elm_type] == dim || elementDim[elm_type] == (dim-1) ) ) )
430 std::string formatString =
"%d";
431 for (
int i=1; i<nDofs[elm_type]; i++)
432 formatString +=
" %d";
433 formatString +=
"\n";
436 std::vector<int> elementDofs(10);
438 readfile(file,nDofs[elm_type], formatString.c_str(),
439 &(elementDofs[0]),&(elementDofs[1]),&(elementDofs[2]),
440 &(elementDofs[3]),&(elementDofs[4]),&(elementDofs[5]),
441 &(elementDofs[6]),&(elementDofs[7]),&(elementDofs[8]),
445 for (
int i=0; i<nVertices[elm_type]; i++)
446 if (renumber.find(elementDofs[i])==renumber.end())
448 renumber[elementDofs[i]] = number_of_real_vertices++;
453 if (elementDim[elm_type] == dim)
456 boundary_element_count++;
463 template <
class E,
class V,
class V2>
466 const E& elementDofs,
470 DUNE_THROW(Dune::IOError,
"tried to create a 3D boundary segment in a non-3D Grid");
474 template <
class E,
class V>
476 const std::vector<FieldVector<double, 3> >& nodes,
477 const E& elementDofs,
481 std::array<FieldVector<double,dimWorld>, 6> v;
482 for (
int i=0; i<6; i++)
484 v[i][j] = nodes[elementDofs[i]][j];
501 std::map<int,unsigned int> &
renumber,
502 const std::vector< GlobalVector > & nodes,
503 const int physical_entity)
506 const int nDofs[12] = {-1, 2, 3, 4, 4, 8, 6, 5, 3, 6, -1, 10};
507 const int nVertices[12] = {-1, 2, 3, 4, 4, 8, 6, 5, 2, 3, -1, 4};
508 const int elementDim[12] = {-1, 1, 2, 2, 3, 3, 3, 3, 1, 2, -1, 3};
511 if ( not (elm_type >= 0 && elm_type < 12
512 && (elementDim[elm_type] == dim || elementDim[elm_type] == (dim-1) ) ) )
519 std::string formatString =
"%d";
520 for (
int i=1; i<nDofs[elm_type]; i++)
521 formatString +=
" %d";
522 formatString +=
"\n";
525 std::vector<int> elementDofs(10);
527 readfile(file,nDofs[elm_type], formatString.c_str(),
528 &(elementDofs[0]),&(elementDofs[1]),&(elementDofs[2]),
529 &(elementDofs[3]),&(elementDofs[4]),&(elementDofs[5]),
530 &(elementDofs[6]),&(elementDofs[7]),&(elementDofs[8]),
537 std::swap(elementDofs[2],elementDofs[3]);
540 std::swap(elementDofs[2],elementDofs[3]);
541 std::swap(elementDofs[6],elementDofs[7]);
544 std::swap(elementDofs[2],elementDofs[3]);
550 std::vector<unsigned int> vertices(nVertices[elm_type]);
552 for (
int i=0; i<nVertices[elm_type]; i++)
553 vertices[i] = renumber[elementDofs[i]];
556 if (elementDim[elm_type] == dim) {
591 if (insert_boundary_segments) {
604 std::array<FieldVector<double,dimWorld>, 3> v;
606 v[0][i] = nodes[elementDofs[0]][i];
607 v[1][i] = nodes[elementDofs[2]][i];
608 v[2][i] = nodes[elementDofs[1]][i];
617 boundarysegment_insert(nodes, elementDofs, vertices);
627 if (elementDim[elm_type] == dim) {
628 element_index_to_physical_entity[element_count] = physical_entity;
631 boundary_id_to_physical_entity[boundary_element_count] = physical_entity;
632 boundary_element_count++;
655 template<
typename Gr
idType>
662 static Grid*
read (
const std::string& fileName,
bool verbose =
true,
bool insertBoundarySegments=
true)
669 parser.
read(fileName);
675 static Grid*
read (
const std::string& fileName,
676 std::vector<int>& boundarySegmentToPhysicalEntity,
677 std::vector<int>& elementToPhysicalEntity,
678 bool verbose =
true,
bool insertBoundarySegments=
true)
685 parser.
read(fileName);
687 boundarySegmentToPhysicalEntity.swap(parser.
boundaryIdMap());
695 bool verbose =
true,
bool insertBoundarySegments=
true)
699 parser.
read(fileName);
704 const std::string& fileName,
705 std::vector<int>& boundarySegmentToPhysicalEntity,
706 std::vector<int>& elementToPhysicalEntity,
707 bool verbose =
true,
bool insertBoundarySegments=
true)
711 parser.
read(fileName);
713 boundarySegmentToPhysicalEntity.swap(parser.
boundaryIdMap());
Include standard header files.
Definition: agrid.hh:58
quadratic boundary approximation.
Definition: gmshreader.hh:40
void skipline(FILE *file)
Definition: gmshreader.hh:246
static Grid * read(const std::string &fileName, std::vector< int > &boundarySegmentToPhysicalEntity, std::vector< int > &elementToPhysicalEntity, bool verbose=true, bool insertBoundarySegments=true)
Definition: gmshreader.hh:675
Base class for grid boundary segments of arbitrary geometry.
Definition: common.hh:186
virtual void insertVertex(const FieldVector< ctype, dimworld > &pos)
Insert a vertex into the coarse grid.
Definition: common/gridfactory.hh:279
std::vector< int > element_index_to_physical_entity
Definition: gmshreader.hh:218
GeometryOrder
Definition: gmshreader.hh:36
static void read(Dune::GridFactory< Grid > &factory, const std::string &fileName, bool verbose=true, bool insertBoundarySegments=true)
Definition: gmshreader.hh:694
virtual GridType * createGrid()
Finalize grid creation and hand over the grid.
Definition: common/gridfactory.hh:316
GmshReaderParser(Dune::GridFactory< GridType > &_factory, bool v, bool i)
Definition: gmshreader.hh:256
std::string fileName
Definition: gmshreader.hh:215
Read Gmsh mesh file.
Definition: gmshreader.hh:656
Provide a generic factory class for unstructured grids.
Definition: common/gridfactory.hh:263
std::vector< int > boundary_id_to_physical_entity
Definition: gmshreader.hh:217
std::vector< int > & boundaryIdMap()
Definition: gmshreader.hh:259
Definition: common.hh:180
FieldVector< double, dimWorld > GlobalVector
Definition: gmshreader.hh:223
dimension independent parts for GmshReaderParser
Definition: gmshreader.hh:203
Definition: common.hh:185
void read(const std::string &f)
Definition: gmshreader.hh:269
bool verbose
Definition: gmshreader.hh:208
static Grid * read(const std::string &fileName, bool verbose=true, bool insertBoundarySegments=true)
Definition: gmshreader.hh:662
std::vector< int > & elementIndexMap()
Definition: gmshreader.hh:264
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:500
Definition: common.hh:183
bool insert_boundary_segments
Definition: gmshreader.hh:209
int element_count
Definition: gmshreader.hh:212
void boundarysegment_insert(const std::vector< FieldVector< double, 3 > > &nodes, const E &elementDofs, const V &vertices)
Definition: gmshreader.hh:475
GridType Grid
Definition: gmshreader.hh:659
ALBERTA REAL_D GlobalVector
Definition: misc.hh:48
static void read(Dune::GridFactory< Grid > &factory, const std::string &fileName, std::vector< int > &boundarySegmentToPhysicalEntity, std::vector< int > &elementToPhysicalEntity, bool verbose=true, bool insertBoundarySegments=true)
Definition: gmshreader.hh:703
static const int dimWorld
Definition: misc.hh:44
unsigned int number_of_real_vertices
Definition: gmshreader.hh:210
Definition: common.hh:187
virtual void insertElement(const GeometryType &type, const std::vector< unsigned int > &vertices)
Insert an element into the coarse grid.
Definition: common/gridfactory.hh:290
Provide a generic factory class for unstructured grids.
Definition: common.hh:184
edges are straight lines.
Definition: gmshreader.hh:38
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:412
int boundary_element_count
Definition: gmshreader.hh:211
int renumber(const Dune::GeometryType &t, int i)
renumber VTK <-> Dune
Definition: common.hh:232
Dune::GridFactory< GridType > & factory
Definition: gmshreader.hh:207
Base class for classes implementing geometries of boundary segments.
Definition: boundarysegment.hh:29
void boundarysegment_insert(const V &nodes, const E &elementDofs, const V2 &vertices)
Definition: gmshreader.hh:464
Definition: common.hh:181
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:232
virtual void insertBoundarySegment(const std::vector< unsigned int > &vertices)
insert a boundary segment
Definition: common/gridfactory.hh:308