4#ifndef DUNE_GMSHREADER_HH
5#define DUNE_GMSHREADER_HH
48 template<
int dimension,
int dimWorld = dimension >
49 class GmshReaderQuadraticBoundarySegment
64 template<
int dimWorld >
65 struct GmshReaderQuadraticBoundarySegment< 2, dimWorld >
70 GmshReaderQuadraticBoundarySegment (
const GlobalVector &p0_,
const GlobalVector &p1_,
const GlobalVector &p2_)
71 : p0(p0_), p1(p1_), p2(p2_)
78 alpha=d1.two_norm()/(d1.two_norm()+d2.two_norm());
79 if (alpha<1E-6 || alpha>1-1E-6)
87 y.axpy((local[0]-alpha)*(local[0]-1.0)/alpha,p0);
88 y.axpy(local[0]*(local[0]-1.0)/(alpha*(alpha-1.0)),p1);
89 y.axpy(local[0]*(local[0]-alpha)/(1.0-alpha),p2);
94 GlobalVector p0,p1,p2;
123 class GmshReaderQuadraticBoundarySegment< 3, 3 >
130 : p0(p0_), p1(p1_), p2(p2_), p3(p3_), p4(p4_), p5(p5_)
138 if (alpha<1E-6 || alpha>1-1E-6)
144 if (beta<1E-6 || beta>1-1E-6)
150 if (gamma<1E-6 || gamma>1-1E-6)
158 y.
axpy(phi0(local),p0);
159 y.
axpy(phi1(local),p1);
160 y.
axpy(phi2(local),p2);
161 y.
axpy(phi3(local),p3);
162 y.
axpy(phi4(local),p4);
163 y.
axpy(phi5(local),p5);
173 return (alpha*beta-beta*local[0]-alpha*local[1])*(1-local[0]-local[1])/(alpha*beta);
177 return local[0]*(1-local[0]-local[1])/(alpha*(1-alpha));
181 return local[0]*(gamma*local[0]-(sqrt2-gamma-sqrt2*alpha)*local[1]-alpha*gamma)/(gamma*(1-alpha));
185 return local[1]*(1-local[0]-local[1])/(beta*(1-beta));
189 return local[0]*local[1]/((1-gamma/sqrt2)*gamma/sqrt2);
193 return local[1]*(beta*(1-gamma/sqrt2)-local[0]*(beta-gamma/sqrt2)-local[1]*(1-gamma/sqrt2))/((1-gamma/sqrt2)*(beta-1));
197 double alpha,beta,gamma,sqrt2;
203 template<
typename Gr
idType>
210 bool insert_boundary_segments;
211 unsigned int number_of_real_vertices;
212 int boundary_element_count;
216 std::string fileName;
218 std::vector<int> boundary_id_to_physical_entity;
219 std::vector<int> element_index_to_physical_entity;
222 static const int dim = GridType::dimension;
223 static const int dimWorld = GridType::dimensionworld;
224 static_assert( (dimWorld <= 3),
"GmshReader requires dimWorld <= 3." );
233 void readfile(FILE * file,
int cnt,
const char * format,
234 void* t1,
void* t2 = 0,
void* t3 = 0,
void* t4 = 0,
235 void* t5 = 0,
void* t6 = 0,
void* t7 = 0,
void* t8 = 0,
236 void* t9 = 0,
void* t10 = 0)
238 off_t pos = ftello(file);
239 int c = fscanf(file, format, t1, t2, t3, t4, t5, t6, t7, t8, t9, t10);
243 <<
": Expected '" << format <<
"', only read " << c <<
" entries instead of " << cnt <<
".");
247 void skipline(FILE * file)
251 c = std::fgetc(file);
252 }
while(c !=
'\n' && c != EOF);
258 factory(_factory), verbose(v), insert_boundary_segments(i) {}
260 std::vector<int> & boundaryIdMap()
262 return boundary_id_to_physical_entity;
265 std::vector<int> & elementIndexMap()
267 return element_index_to_physical_entity;
270 void read (
const std::string& f)
272 if (verbose) std::cout <<
"Reading " << dim <<
"d Gmsh grid..." << std::endl;
276 FILE* file = fopen(fileName.c_str(),
"r");
285 number_of_real_vertices = 0;
286 boundary_element_count = 0;
290 double version_number;
291 int file_type, data_size;
293 readfile(file,1,
"%s\n",buf);
294 if (strcmp(buf,
"$MeshFormat")!=0)
296 readfile(file,3,
"%lg %d %d\n",&version_number,&file_type,&data_size);
297 if( (version_number < 2.0) || (version_number > 2.2) )
299 if (verbose) std::cout <<
"version " << version_number <<
" Gmsh file detected" << std::endl;
300 readfile(file,1,
"%s\n",buf);
301 if (strcmp(buf,
"$EndMeshFormat")!=0)
307 readfile(file,1,
"%s\n",buf);
308 if (strcmp(buf,
"$Nodes")!=0)
310 readfile(file,1,
"%d\n",&number_of_nodes);
311 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) {
324 "Only dense sequences of node indices are currently supported (node index "
325 <<
id <<
" is invalid).");
329 for(
int j = 0; j < dimWorld; ++j )
330 nodes[
id ][ j ] = x[ j ];
332 readfile(file,1,
"%s\n",buf);
333 if (strcmp(buf,
"$EndNodes")!=0)
338 readfile(file,1,
"%s\n",buf);
339 if (strcmp(buf,
"$Elements")!=0)
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;
349 long section_element_offset = ftell(file);
350 std::map<int,unsigned int> renumber;
351 for (
int i=1; i<=number_of_elements; i++)
353 int id, elm_type, number_of_tags;
354 readfile(file,3,
"%d %d %d ",&
id,&elm_type,&number_of_tags);
355 for (
int k=1; k<=number_of_tags; k++)
358 readfile(file,1,
"%d ",&blub);
367 pass1HandleElement(file, elm_type, renumber, nodes);
369 if (verbose) std::cout <<
"number of real vertices = " << number_of_real_vertices << std::endl;
370 if (verbose) std::cout <<
"number of boundary elements = " << boundary_element_count << std::endl;
371 if (verbose) std::cout <<
"number of elements = " << element_count << std::endl;
372 readfile(file,1,
"%s\n",buf);
373 if (strcmp(buf,
"$EndElements")!=0)
375 boundary_id_to_physical_entity.resize(boundary_element_count);
376 element_index_to_physical_entity.resize(element_count);
382 fseek(file, section_element_offset, SEEK_SET);
383 boundary_element_count = 0;
385 for (
int i=1; i<=number_of_elements; i++)
387 int id, elm_type, number_of_tags;
388 readfile(file,3,
"%d %d %d ",&
id,&elm_type,&number_of_tags);
389 int physical_entity = -1;
390 std::vector<int> mesh_partitions;
391 if ( version_number < 2.2 )
393 mesh_partitions.resize(1);
395 for (
int k=1; k<=number_of_tags; k++)
398 readfile(file,1,
"%d ",&blub);
399 if (k==1) physical_entity = blub;
401 if ( version_number < 2.2 )
403 if (k==3) mesh_partitions[0] = blub;
408 mesh_partitions[k-4] = blub;
410 mesh_partitions.resize(blub);
413 pass2HandleElement(file, elm_type, renumber, nodes, physical_entity);
415 readfile(file,1,
"%s\n",buf);
416 if (strcmp(buf,
"$EndElements")!=0)
423 void pass1HandleElement(FILE* file,
const int elm_type,
424 std::map<int,unsigned int> & renumber,
425 const std::vector< GlobalVector > & nodes)
428 const int nDofs[12] = {-1, 2, 3, 4, 4, 8, 6, 5, 3, 6, -1, 10};
429 const int nVertices[12] = {-1, 2, 3, 4, 4, 8, 6, 5, 2, 3, -1, 4};
430 const int elementDim[12] = {-1, 1, 2, 2, 3, 3, 3, 3, 1, 2, -1, 3};
433 if ( not (elm_type >= 0 && elm_type < 12
434 && (elementDim[elm_type] == dim || elementDim[elm_type] == (dim-1) ) ) )
442 for (
int i=1; i<nDofs[elm_type]; i++)
447 std::vector<int> elementDofs(10);
450 &(elementDofs[0]),&(elementDofs[1]),&(elementDofs[2]),
451 &(elementDofs[3]),&(elementDofs[4]),&(elementDofs[5]),
452 &(elementDofs[6]),&(elementDofs[7]),&(elementDofs[8]),
456 for (
int i=0; i<nVertices[elm_type]; i++)
457 if (renumber.find(elementDofs[i])==renumber.end())
459 renumber[elementDofs[i]] = number_of_real_vertices++;
464 if (elementDim[elm_type] == dim)
467 boundary_element_count++;
474 template <
class E,
class V,
class V2>
475 void boundarysegment_insert(
477 const E& elementDofs,
485 template <
class E,
class V>
486 void boundarysegment_insert(
488 const E& elementDofs,
492 array<FieldVector<double,dimWorld>, 6> v;
493 for (
int i=0; i<6; i++)
494 for (
int j=0; j<dimWorld; j++)
495 v[i][j] = nodes[elementDofs[i]][j];
509 virtual void pass2HandleElement(FILE* file,
const int elm_type,
510 std::map<int,unsigned int> & renumber,
511 const std::vector< GlobalVector > & nodes,
512 const int physical_entity)
515 const int nDofs[12] = {-1, 2, 3, 4, 4, 8, 6, 5, 3, 6, -1, 10};
516 const int nVertices[12] = {-1, 2, 3, 4, 4, 8, 6, 5, 2, 3, -1, 4};
517 const int elementDim[12] = {-1, 1, 2, 2, 3, 3, 3, 3, 1, 2, -1, 3};
520 if ( not (elm_type >= 0 && elm_type < 12
521 && (elementDim[elm_type] == dim || elementDim[elm_type] == (dim-1) ) ) )
529 for (
int i=1; i<nDofs[elm_type]; i++)
534 std::vector<int> elementDofs(10);
537 &(elementDofs[0]),&(elementDofs[1]),&(elementDofs[2]),
538 &(elementDofs[3]),&(elementDofs[4]),&(elementDofs[5]),
539 &(elementDofs[6]),&(elementDofs[7]),&(elementDofs[8]),
546 std::swap(elementDofs[2],elementDofs[3]);
549 std::swap(elementDofs[2],elementDofs[3]);
550 std::swap(elementDofs[6],elementDofs[7]);
553 std::swap(elementDofs[2],elementDofs[3]);
559 std::vector<unsigned int> vertices(nVertices[elm_type]);
561 for (
int i=0; i<nVertices[elm_type]; i++)
562 vertices[i] = renumber[elementDofs[i]];
565 if (elementDim[elm_type] == dim) {
600 if (insert_boundary_segments) {
613 array<FieldVector<double,dimWorld>, 3> v;
614 for (
int i=0; i<dimWorld; i++) {
615 v[0][i] = nodes[elementDofs[0]][i];
616 v[1][i] = nodes[elementDofs[2]][i];
617 v[2][i] = nodes[elementDofs[1]][i];
626 boundarysegment_insert(nodes, elementDofs, vertices);
636 if (elementDim[elm_type] == dim) {
637 element_index_to_physical_entity[element_count] = physical_entity;
640 boundary_id_to_physical_entity[boundary_element_count] = physical_entity;
641 boundary_element_count++;
664 template<
typename Gr
idType>
668 typedef GridType Grid;
671 static Grid*
read (
const std::string& fileName,
bool verbose =
true,
bool insert_boundary_segments=
true)
678 parser.read(fileName);
684 static Grid*
read (
const std::string& fileName,
685 std::vector<int>& boundary_id_to_physical_entity,
686 std::vector<int>& element_index_to_physical_entity,
687 bool verbose =
true,
bool insert_boundary_segments=
true)
694 parser.read(fileName);
696 boundary_id_to_physical_entity.swap(parser.boundaryIdMap());
697 element_index_to_physical_entity.swap(parser.elementIndexMap());
704 bool verbose =
true,
bool insert_boundary_segments=
true)
708 parser.read(fileName);
713 const std::string& fileName,
714 std::vector<int>& boundary_id_to_physical_entity,
715 std::vector<int>& element_index_to_physical_entity,
716 bool verbose =
true,
bool insert_boundary_segments=
true)
720 parser.read(fileName);
722 boundary_id_to_physical_entity.swap(parser.boundaryIdMap());
723 element_index_to_physical_entity.swap(parser.elementIndexMap());
Base class for grid boundary segments of arbitrary geometry.
FieldTraits< value_type >::real_type two_norm() const
two norm sqrt(sum over squared values of entries)
Definition: densevector.hh:589
derived_type & axpy(const value_type &a, const DenseVector< Other > &y)
vector space axpy operation ( *this += a y )
Definition: densevector.hh:523
vector space out of a tensor product of fields.
Definition: fvector.hh:94
Unique label for each type of entities that can occur in DUNE grids.
Definition: type.hh:25
@ cube
Cube element in any nonnegative dimension.
Definition: type.hh:31
@ simplex
Simplicial element in any nonnegative dimension.
Definition: type.hh:30
@ pyramid
Four sided pyramid in three dimensions.
Definition: type.hh:32
@ prism
Prism element in three dimensions.
Definition: type.hh:33
dimension independent parts for GmshReaderParser
Definition: gmshreader.hh:205
Read Gmsh mesh file.
Definition: gmshreader.hh:666
static Grid * read(const std::string &fileName, bool verbose=true, bool insert_boundary_segments=true)
Definition: gmshreader.hh:671
static void read(Dune::GridFactory< Grid > &factory, const std::string &fileName, std::vector< int > &boundary_id_to_physical_entity, std::vector< int > &element_index_to_physical_entity, bool verbose=true, bool insert_boundary_segments=true)
Definition: gmshreader.hh:712
static void read(Dune::GridFactory< Grid > &factory, const std::string &fileName, bool verbose=true, bool insert_boundary_segments=true)
Definition: gmshreader.hh:703
static Grid * read(const std::string &fileName, std::vector< int > &boundary_id_to_physical_entity, std::vector< int > &element_index_to_physical_entity, bool verbose=true, bool insert_boundary_segments=true)
Definition: gmshreader.hh:684
Provide a generic factory class for unstructured grids.
Definition: gridfactory.hh:263
virtual void insertElement(const GeometryType &type, const std::vector< unsigned int > &vertices)
Insert an element into the coarse grid.
Definition: gridfactory.hh:290
virtual GridType * createGrid()
Finalize grid creation and hand over the grid.
Definition: gridfactory.hh:316
virtual void insertVertex(const FieldVector< ctype, dimworld > &pos)
Insert a vertex into the coarse grid.
Definition: gridfactory.hh:279
virtual void insertBoundarySegment(const std::vector< unsigned int > &vertices)
insert a boundary segment
Definition: gridfactory.hh:308
Default exception class for I/O errors.
Definition: exceptions.hh:256
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.
static std::string formatString(const std::string &s, const T &... args)
Format values according to printf format string.
Definition: stringutility.hh:72
#define DUNE_THROW(E, m)
Definition: exceptions.hh:243
Dune namespace.
Definition: alignment.hh:10
Base class for classes implementing geometries of boundary segments.
Definition: boundarysegment.hh:30
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.