[Dune] gmsh vs. gmesh
Oliver Sander
sander at mi.fu-berlin.de
Fri Jun 5 13:32:37 CEST 2009
When I look at what seems to be the official homepage
http://geuz.org/gmsh/
the spelling is gmsh, not gmesh
--
Oliver
-------- Original-Nachricht --------
Betreff: [Dune-Commit] dune-grid r5232 - trunk/grid/io/file
Datum: Thu, 04 Jun 2009 21:20:25 +0200
Von: peter at dune-project.org
An: dune-commit at dune-project.org
Author: peter
Date: 2009-06-04 21:20:25 +0200 (Thu, 04 Jun 2009)
New Revision: 5232
Added:
trunk/grid/io/file/gmeshreader.hh
Removed:
trunk/grid/io/file/gmshreader.hh
Log:
gmsh is actually called gmesh!
Also the reader class has been renamed.
Copied: trunk/grid/io/file/gmeshreader.hh (from rev 5231, trunk/grid/io/file/gmshreader.hh)
===================================================================
--- trunk/grid/io/file/gmeshreader.hh (rev 0)
+++ trunk/grid/io/file/gmeshreader.hh 2009-06-04 19:20:25 UTC (rev 5232)
@@ -0,0 +1,887 @@
+// -*- tab-width: 4; indent-tabs-mode: nil -*-
+#ifndef DUNE_GMESHREADER_HH
+#define DUNE_GMESHREADER_HH
+
+#include<iostream>
+#include<fstream>
+#include<string>
+#include<vector>
+#include<map>
+#include<stdio.h>
+
+#include<dune/common/geometrytype.hh>
+#include<dune/common/exceptions.hh>
+#include<dune/common/fvector.hh>
+#include<dune/grid/common/gridfactory.hh>
+#include<dune/grid/common/boundarysegment.hh>
+
+namespace Dune {
+
+ /**
+ \ingroup Gmesh
+ \{
+ */
+
+ //! Options for read operation
+ struct GmeshReaderOptions
+ {
+ enum GeometryOrder {
+ /** @brief edges are straight lines. */
+ firstOrder,
+ /** @brief quadratic boundary approximation. */
+ secondOrder
+ };
+ };
+
+ // arbitrary dimension, implementation is in specialization
+ template<int dimension>
+ class GmeshReaderLinearBoundarySegment
+ {
+ };
+
+
+ // linear boundary segments in 1d
+ template<>
+ class GmeshReaderLinearBoundarySegment<2> : public Dune::BoundarySegment<2>
+ {
+ public:
+ GmeshReaderLinearBoundarySegment (Dune::FieldVector<double,2> p0_, Dune::FieldVector<double,2> p1_)
+ : p0(p0_), p1(p1_)
+ {}
+
+ virtual Dune::FieldVector<double,2> operator() (const Dune::FieldVector<double,1>& local) const
+ {
+ Dune::FieldVector<double,2> y;
+ y = 0.0;
+ y.axpy(1.0-local[0],p0);
+ y.axpy(local[0],p1);
+ return y;
+ }
+
+ private:
+ Dune::FieldVector<double,2> p0,p1;
+ };
+
+
+
+ // arbitrary dimension, implementation is in specialization
+ template<int dimension>
+ class GmeshReaderQuadraticBoundarySegment
+ {
+ };
+
+ // quadratic boundary segments in 1d
+ /*
+ Note the points
+
+ (0) (alpha) (1)
+
+ are mapped to the points in global coordinates
+
+ p0 p2 p1
+
+ alpha is determined automatically from the given points.
+ */
+ template <>
+ class GmeshReaderQuadraticBoundarySegment<2> : public Dune::BoundarySegment<2>
+ {
+ public:
+ GmeshReaderQuadraticBoundarySegment (Dune::FieldVector<double,2> p0_, Dune::FieldVector<double,2> p1_,
+ Dune::FieldVector<double,2> p2_)
+ : p0(p0_), p1(p1_), p2(p2_)
+ {
+ Dune::FieldVector<double,2> d1(p1);
+ d1 -= p0;
+ Dune::FieldVector<double,2> d2(p2);
+ d2 -= p1;
+
+ alpha=d1.two_norm()/(d1.two_norm()+d2.two_norm());
+ if (alpha<1E-6 || alpha>1-1E-6)
+ DUNE_THROW(Dune::IOError, "ration in quadratic boundary segment bad");
+ }
+
+ virtual Dune::FieldVector<double,2> operator() (const Dune::FieldVector<double,1>& local) const
+ {
+ Dune::FieldVector<double,2> y;
+ y = 0.0;
+ y.axpy((local[0]-alpha)*(local[0]-1.0)/alpha,p0);
+ y.axpy(local[0]*(local[0]-1.0)/(alpha*(alpha-1.0)),p1);
+ y.axpy(local[0]*(local[0]-alpha)/(1.0-alpha),p2);
+ return y;
+ }
+
+ private:
+ Dune::FieldVector<double,2> p0,p1,p2;
+ double alpha;
+ };
+
+
+
+ // arbitrary dimension, implementation is in specialization
+ template<typename GridType, int dimension>
+ class GmeshReaderImp
+ {
+ };
+
+ template<typename GridType>
+ class GmeshReaderImp<GridType,2>
+ {
+ public:
+
+ template<typename T>
+ static GridType* read(Dune::GridFactory<GridType>& factory, const std::string& fileName,
+ bool verbose, std::vector<T>& boundary_id_to_physical_entity,
+ std::vector<T>& element_index_to_physical_entity)
+ {
+ // the grid dimension
+ const int dim = GridType::dimension;
+
+ // open file name, we use C I/O
+ FILE* file = fopen(fileName.c_str(),"r");
+ if (file==0)
+ DUNE_THROW(Dune::IOError, "Could not open " << fileName);
+
+ // a read buffer
+ char buf[512];
+
+ //=========================================
+ // Pass 1: Read vertices into vector
+ // Check vertices that are needed
+ // Renumber needed vertices
+ //=========================================
+
+ int boundary_element_count = 0;
+ int element_count = 0;
+
+ // process header
+ fscanf(file,"%s\n",buf);
+ if (strcmp(buf,"$MeshFormat")!=0)
+ DUNE_THROW(Dune::IOError, "expected $MeshFormat in first line");
+ int version_number, file_type, data_size;
+ fscanf(file,"%d %d %d\n",&version_number,&file_type,&data_size);
+ if (version_number!=2)
+ DUNE_THROW(Dune::IOError, "can only read version_number==2");
+ if (verbose) std::cout << "version 2 Gmesh file detected" << std::endl;
+ fscanf(file,"%s\n",buf);
+ if (strcmp(buf,"$EndMeshFormat")!=0)
+ DUNE_THROW(Dune::IOError, "expected $EndMeshFormat");
+
+ // node section
+ fscanf(file,"%s\n",buf);
+ if (strcmp(buf,"$Nodes")!=0)
+ DUNE_THROW(Dune::IOError, "expected $Nodes");
+ int number_of_nodes;
+ fscanf(file,"%d\n",&number_of_nodes);
+ if (verbose) std::cout << "file contains " << number_of_nodes << " nodes" << std::endl;
+ std::vector<Dune::FieldVector<double,dim> > nodes(number_of_nodes+1); // store positions
+ for (int i=1; i<=number_of_nodes; i++)
+ {
+ int id;
+ double x,y,z;
+ fscanf(file,"%d %lg %lg %lg\n",&id,&x,&y,&z);
+ // if (verbose) std::cout << id << " " << x << " " << y << " " << z << std::endl;
+ if (id!=i)
+ DUNE_THROW(Dune::IOError, "expected id " << i);
+
+ Dune::FieldVector<double,dim> position;
+ position[0] = x; position[1] = y;
+ nodes[i] = position; // just store node position
+ }
+ fscanf(file,"%s\n",buf);
+ if (strcmp(buf,"$EndNodes")!=0)
+ DUNE_THROW(Dune::IOError, "expected $EndNodes");
+
+ // element section
+ fscanf(file,"%s\n",buf);
+ if (strcmp(buf,"$Elements")!=0)
+ DUNE_THROW(Dune::IOError, "expected $Elements");
+ int number_of_elements;
+ fscanf(file,"%d\n",&number_of_elements);
+ if (verbose) std::cout << "file contains " << number_of_elements << " elements" << std::endl;
+ unsigned int number_of_real_vertices=0; // count number of vertices that are really needed
+ std::map<int,unsigned int> renumber;
+ for (int i=1; i<=number_of_elements; i++)
+ {
+ int id, elm_type, number_of_tags;
+ fscanf(file,"%d %d %d",&id,&elm_type,&number_of_tags);
+ int physical_entity, elementary_entity, mesh_partition;
+ for (int k=1; k<=number_of_tags; k++)
+ {
+ int blub;
+ fscanf(file,"%d",&blub);
+ if (k==1) physical_entity = blub;
+ if (k==2) elementary_entity = blub;
+ if (k==3) mesh_partition = blub;
+ }
+ std::vector<int> simplexVertices(6);
+ switch (elm_type)
+ {
+ case 1: // 2-node line
+ simplexVertices.resize(2);
+ fscanf(file,"%d %d\n",&(simplexVertices[0]),&(simplexVertices[1]));
+ for (int i=0; i<2; i++)
+ if (renumber.find(simplexVertices[i])==renumber.end())
+ {
+ // std::cout << "real vertex " << number_of_real_vertices << " at " << nodes[simplexVertices[i]] << " old number was " << simplexVertices[i] << std::endl;
+ renumber[simplexVertices[i]] = number_of_real_vertices++;
+ factory.insertVertex(nodes[simplexVertices[i]]);
+ }
+ boundary_element_count++;
+ break;
+ case 8: // 3-node line
+ simplexVertices.resize(3);
+ fscanf(file,"%d %d %d\n",&(simplexVertices[0]),&(simplexVertices[1]),&(simplexVertices[2]));
+ for (int i=0; i<2; i++)
+ if (renumber.find(simplexVertices[i])==renumber.end())
+ {
+ renumber[simplexVertices[i]] = number_of_real_vertices++;
+ factory.insertVertex(nodes[simplexVertices[i]]);
+ }
+ boundary_element_count++;
+ break;
+ case 2: // 3-node triangle
+ simplexVertices.resize(3);
+ fscanf(file,"%d %d %d\n",&(simplexVertices[0]),&(simplexVertices[1]),&(simplexVertices[2]));
+ for (int i=0; i<3; i++)
+ if (renumber.find(simplexVertices[i])==renumber.end())
+ {
+ renumber[simplexVertices[i]] = number_of_real_vertices++;
+ factory.insertVertex(nodes[simplexVertices[i]]);
+ }
+ element_count++;
+ break;
+ case 9: // 6-node triangle
+ simplexVertices.resize(6);
+ fscanf(file,"%d %d %d %d %d %d\n",&(simplexVertices[0]),&(simplexVertices[1]),&(simplexVertices[2]),
+ &(simplexVertices[3]),&(simplexVertices[4]),&(simplexVertices[5]));
+ for (int i=0; i<3; i++)
+ if (renumber.find(simplexVertices[i])==renumber.end())
+ {
+ renumber[simplexVertices[i]] = number_of_real_vertices++;
+ factory.insertVertex(nodes[simplexVertices[i]]);
+ }
+ element_count++;
+ break;
+ default:
+ fgets(buf,512,file); // skip rest of line if no tetrahedron
+ }
+ }
+ if (verbose) std::cout << "number of real vertices = " << number_of_real_vertices << std::endl;
+ if (verbose) std::cout << "number of boundary elements = " << boundary_element_count << std::endl;
+ if (verbose) std::cout << "number of elements = " << element_count << std::endl;
+ fscanf(file,"%s\n",buf);
+ if (strcmp(buf,"$EndElements")!=0)
+ DUNE_THROW(Dune::IOError, "expected $EndElements");
+ boundary_id_to_physical_entity.resize(boundary_element_count);
+ element_index_to_physical_entity.resize(element_count);
+
+ //==============================================
+ // Pass 2: Insert boundary segments and elements
+ //==============================================
+
+ // go to beginning of file
+ rewind(file);
+ boundary_element_count = 0;
+ element_count = 0;
+
+ // process header
+ fscanf(file,"%s\n",buf);
+ if (strcmp(buf,"$MeshFormat")!=0)
+ DUNE_THROW(Dune::IOError, "expected $MeshFormat in first line");
+ fscanf(file,"%d %d %d\n",&version_number,&file_type,&data_size);
+ if (version_number!=2)
+ DUNE_THROW(Dune::IOError, "can only read version_number==2");
+ fscanf(file,"%s\n",buf);
+ if (strcmp(buf,"$EndMeshFormat")!=0)
+ DUNE_THROW(Dune::IOError, "expected $EndMeshFormat");
+
+ // node section
+ fscanf(file,"%s\n",buf);
+ if (strcmp(buf,"$Nodes")!=0)
+ DUNE_THROW(Dune::IOError, "expected $Nodes");
+ fscanf(file,"%d\n",&number_of_nodes);
+ for (int i=1; i<=number_of_nodes; i++)
+ {
+ int id;
+ double x,y,z;
+ fscanf(file,"%d %lg %lg %lg\n",&id,&x,&y,&z);
+ if (id!=i)
+ DUNE_THROW(Dune::IOError, "expected id " << i);
+ }
+ fscanf(file,"%s\n",buf);
+ if (strcmp(buf,"$EndNodes")!=0)
+ DUNE_THROW(Dune::IOError, "expected $EndNodes");
+
+ // element section
+ fscanf(file,"%s\n",buf);
+ if (strcmp(buf,"$Elements")!=0)
+ DUNE_THROW(Dune::IOError, "expected $Elements");
+ fscanf(file,"%d\n",&number_of_elements);
+ for (int i=1; i<=number_of_elements; i++)
+ {
+ int id, elm_type, number_of_tags;
+ fscanf(file,"%d %d %d",&id,&elm_type,&number_of_tags);
+ int physical_entity, elementary_entity, mesh_partition;
+ for (int k=1; k<=number_of_tags; k++)
+ {
+ int blub;
+ fscanf(file,"%d",&blub);
+ if (k==1) physical_entity = blub;
+ if (k==2) elementary_entity = blub;
+ if (k==3) mesh_partition = blub;
+ }
+ std::vector<int> simplexVertices(6);
+ std::vector<unsigned int> vertices(3);
+ switch (elm_type)
+ {
+ case 1: // 2-node line
+
+ // read the vertices, but don't do anything with them. Linear boundary
+ // segments are the default.
+ simplexVertices.resize(2);
+ fscanf(file,"%d %d\n",&(simplexVertices[0]),&(simplexVertices[1]));
+ vertices.resize(2);
+ for (int i=0; i<2; i++)
+ vertices[i] = renumber[simplexVertices[i]]; // renumber vertices
+ factory.insertBoundarySegment(vertices,new GmeshReaderLinearBoundarySegment<2>(nodes[simplexVertices[0]],nodes[simplexVertices[1]]));
+ boundary_id_to_physical_entity[boundary_element_count] = physical_entity;
+ boundary_element_count++;
+ break;
+ case 8: // 3-node line
+ simplexVertices.resize(3);
+ fscanf(file,"%d %d %d\n",&(simplexVertices[0]),&(simplexVertices[1]),&(simplexVertices[2]));
+ vertices.resize(2);
+ for (int i=0; i<2; i++)
+ vertices[i] = renumber[simplexVertices[i]]; // renumber vertices
+ factory.insertBoundarySegment(vertices,new GmeshReaderQuadraticBoundarySegment<2>(nodes[simplexVertices[0]],
+ nodes[simplexVertices[2]],
+ nodes[simplexVertices[1]]));
+ boundary_id_to_physical_entity[boundary_element_count] = physical_entity;
+ boundary_element_count++;
+ break;
+ case 2: // 3-node triangle
+ simplexVertices.resize(3);
+ fscanf(file,"%d %d %d\n",&(simplexVertices[0]),&(simplexVertices[1]),&(simplexVertices[2]));
+ vertices.resize(3);
+ for (int i=0; i<3; i++)
+ vertices[i] = renumber[simplexVertices[i]]; // renumber vertices
+ factory.insertElement(Dune::GeometryType(Dune::GeometryType::simplex,dim),vertices);
+ element_index_to_physical_entity[element_count] = physical_entity;
+ element_count++;
+ break;
+ case 9: // 6-node triangle
+ simplexVertices.resize(6);
+ fscanf(file,"%d %d %d %d %d %d\n",&(simplexVertices[0]),&(simplexVertices[1]),&(simplexVertices[2]),
+ &(simplexVertices[3]),&(simplexVertices[4]),&(simplexVertices[5]));
+ vertices.resize(3);
+ for (int i=0; i<3; i++)
+ vertices[i] = renumber[simplexVertices[i]]; // renumber vertices
+ factory.insertElement(Dune::GeometryType(Dune::GeometryType::simplex,dim),vertices);
+ element_index_to_physical_entity[element_count] = physical_entity;
+ element_count++;
+ break;
+ default:
+ fgets(buf,512,file); // skip rest of line if no tetrahedron
+ }
+ }
+ fscanf(file,"%s\n",buf);
+ if (strcmp(buf,"$EndElements")!=0)
+ DUNE_THROW(Dune::IOError, "expected $EndElements");
+
+ fclose(file);
+
+ return factory.createGrid();
+ }
+ };
+
+
+ // linear boundary segments in 2d
+ template <>
+ class GmeshReaderLinearBoundarySegment<3> : public Dune::BoundarySegment<3>
+ {
+ public:
+ GmeshReaderLinearBoundarySegment (Dune::FieldVector<double,3> p0_, Dune::FieldVector<double,3> p1_,
+ Dune::FieldVector<double,3> p2_)
+ : p0(p0_), p1(p1_), p2(p2_)
+ {
+ // std::cout << "created boundary segment " << p0 << " | " << p1 << " | " << p2 << std::endl;
+ }
+
+ virtual Dune::FieldVector<double,3> operator() (const Dune::FieldVector<double,2>& local) const
+ {
+ Dune::FieldVector<double,3> y;
+ y = 0.0;
+ y.axpy(1.0-local[0]-local[1],p0);
+ y.axpy(local[0],p1);
+ y.axpy(local[1],p2);
+ // std::cout << "eval boundary segment local=" << local << " y=" << y << std::endl;
+ return y;
+ }
+
+ private:
+ Dune::FieldVector<double,3> p0,p1,p2;
+ };
+
+
+ // quadratic boundary segments in 2d
+ /* numbering of points corresponding to gmsh:
+
+ 2
+
+ 5 4
+
+ 0 3 1
+
+ Note: The vertices 3, 4, 5 are not necessarily at the edge midpoints but can
+ be placed with parameters alpha, beta , gamma at the following positions
+ in local coordinates:
+
+
+ 2 = (0,1)
+
+ 5 = (0,beta) 4 = (1-gamma/sqrt(2),gamma/sqrt(2))
+
+ 0 = (0,0) 3 = (alpha,0) 1 = (1,0)
+
+ The parameters alpha, beta, gamma are determined from the given vertices in
+ global coordinates.
+ */
+ template <>
+ class GmeshReaderQuadraticBoundarySegment<3> : public Dune::BoundarySegment<3>
+ {
+ public:
+ GmeshReaderQuadraticBoundarySegment (Dune::FieldVector<double,3> p0_, Dune::FieldVector<double,3> p1_,
+ Dune::FieldVector<double,3> p2_, Dune::FieldVector<double,3> p3_,
+ Dune::FieldVector<double,3> p4_, Dune::FieldVector<double,3> p5_)
+ : p0(p0_), p1(p1_), p2(p2_), p3(p3_), p4(p4_), p5(p5_)
+ {
+ Dune::FieldVector<double,3> d1,d2;
+
+ d1 = p3; d1 -= p0;
+ d2 = p1; d2 -= p3;
+ alpha=d1.two_norm()/(d1.two_norm()+d2.two_norm());
+ if (alpha<1E-6 || alpha>1-1E-6)
+ DUNE_THROW(Dune::IOError, "alpha in quadratic boundary segment bad");
+
+ d1 = p5; d1 -= p0;
+ d2 = p2; d2 -= p5;
+ beta=d1.two_norm()/(d1.two_norm()+d2.two_norm());
+ if (beta<1E-6 || beta>1-1E-6)
+ DUNE_THROW(Dune::IOError, "beta in quadratic boundary segment bad");
+
+ d1 = p4; d1 -= p1;
+ d2 = p2; d2 -= p4;
+ gamma=d1.two_norm()/(d1.two_norm()+d2.two_norm());
+ if (gamma<1E-6 || gamma>1-1E-6)
+ DUNE_THROW(Dune::IOError, "gamma in quadratic boundary segment bad");
+
+ sqrt2 = sqrt(2.0);
+ }
+
+ virtual Dune::FieldVector<double,3> operator() (const Dune::FieldVector<double,2>& local) const
+ {
+ Dune::FieldVector<double,3> y;
+ y = 0.0;
+ y.axpy(phi0(local),p0);
+ y.axpy(phi1(local),p1);
+ y.axpy(phi2(local),p2);
+ y.axpy(phi3(local),p3);
+ y.axpy(phi4(local),p4);
+ y.axpy(phi5(local),p5);
+ return y;
+ }
+
+ private:
+ // The six Lagrange basis function on the reference element
+ // for the points given above
+
+ double phi0 (const Dune::FieldVector<double,2>& local) const
+ {
+ return (alpha*beta-beta*local[0]-alpha*local[1])*(1-local[0]-local[1])/(alpha*beta);
+ }
+ double phi3 (const Dune::FieldVector<double,2>& local) const
+ {
+ return local[0]*(1-local[0]-local[1])/(alpha*(1-alpha));
+ }
+ double phi1 (const Dune::FieldVector<double,2>& local) const
+ {
+ return local[0]*(gamma*local[0]-(sqrt2-gamma-sqrt2*alpha)*local[1]-alpha*gamma)/(gamma*(1-alpha));
+ }
+ double phi5 (const Dune::FieldVector<double,2>& local) const
+ {
+ return local[1]*(1-local[0]-local[1])/(beta*(1-beta));
+ }
+ double phi4 (const Dune::FieldVector<double,2>& local) const
+ {
+ return local[0]*local[1]/((1-gamma/sqrt2)*gamma/sqrt2);
+ }
+ double phi2 (const Dune::FieldVector<double,2>& local) const
+ {
+ return local[1]*(beta*(1-gamma/sqrt2)-local[0]*(beta-gamma/sqrt2)-local[1]*(1-gamma/sqrt2))/((1-gamma/sqrt2)*(beta-1));
+ }
+
+ Dune::FieldVector<double,3> p0,p1,p2,p3,p4,p5;
+ double alpha,beta,gamma,sqrt2;
+ };
+
+
+
+ template<typename GridType>
+ class GmeshReaderImp<GridType,3>
+ {
+ public:
+
+ template<typename T>
+ static GridType* read(Dune::GridFactory<GridType>& factory, const std::string& fileName,
+ bool verbose, std::vector<T>& boundary_id_to_physical_entity,
+ std::vector<T>& element_index_to_physical_entity)
+ {
+ // the grid dimension
+ const int dim = GridType::dimension;
+
+ // open file name, we use C I/O
+ FILE* file = fopen(fileName.c_str(),"r");
+ if (file==0)
+ DUNE_THROW(Dune::IOError, "Could not open " << fileName);
+
+ // a read buffer
+ char buf[512];
+
+ //=========================================
+ // Pass 1: Read vertices into vector
+ // Check vertices that are needed
+ // Renumber needed vertices
+ //=========================================
+
+ int boundary_element_count = 0;
+ int element_count = 0;
+
+ // process header
+ fscanf(file,"%s\n",buf);
+ if (strcmp(buf,"$MeshFormat")!=0)
+ DUNE_THROW(Dune::IOError, "expected $MeshFormat in first line");
+ int version_number, file_type, data_size;
+ fscanf(file,"%d %d %d\n",&version_number,&file_type,&data_size);
+ if (version_number!=2)
+ DUNE_THROW(Dune::IOError, "can only read version_number==2");
+ if (verbose) std::cout << "version 2 Gmesh file detected" << std::endl;
+ fscanf(file,"%s\n",buf);
+ if (strcmp(buf,"$EndMeshFormat")!=0)
+ DUNE_THROW(Dune::IOError, "expected $EndMeshFormat");
+
+ // node section
+ fscanf(file,"%s\n",buf);
+ if (strcmp(buf,"$Nodes")!=0)
+ DUNE_THROW(Dune::IOError, "expected $Nodes");
+ int number_of_nodes;
+ fscanf(file,"%d\n",&number_of_nodes);
+ if (verbose) std::cout << "file contains " << number_of_nodes << " nodes" << std::endl;
+ std::vector<Dune::FieldVector<double,dim> > nodes(number_of_nodes+1); // store positions
+ for (int i=1; i<=number_of_nodes; i++) // note: ids are starting with 1!
+ {
+ int id;
+ double x,y,z;
+ fscanf(file,"%d %lg %lg %lg\n",&id,&x,&y,&z);
+ // if (verbose) std::cout << id << " " << x << " " << y << " " << z << std::endl;
+ if (id!=i)
+ DUNE_THROW(Dune::IOError, "expected id " << i);
+
+ Dune::FieldVector<double,dim> position;
+ position[0] = x; position[1] = y; position[2] = z;
+ nodes[i] = position; // just store node position
+ }
+ fscanf(file,"%s\n",buf);
+ if (strcmp(buf,"$EndNodes")!=0)
+ DUNE_THROW(Dune::IOError, "expected $EndNodes");
+
+ // element section
+ fscanf(file,"%s\n",buf);
+ if (strcmp(buf,"$Elements")!=0)
+ DUNE_THROW(Dune::IOError, "expected $Elements");
+ int number_of_elements;
+ fscanf(file,"%d\n",&number_of_elements);
+ if (verbose) std::cout << "file contains " << number_of_elements << " elements" << std::endl;
+ unsigned int number_of_real_vertices=0; // count number of vertices that are really needed
+ std::map<int,unsigned int> renumber;
+ for (int i=1; i<=number_of_elements; i++) // note: ids are starting with 1!
+ {
+ int id, elm_type, number_of_tags;
+ fscanf(file,"%d %d %d",&id,&elm_type,&number_of_tags);
+ int physical_entity, elementary_entity, mesh_partition;
+ for (int k=1; k<=number_of_tags; k++)
+ {
+ int blub;
+ fscanf(file,"%d",&blub);
+ if (k==1) physical_entity = blub;
+ if (k==2) elementary_entity = blub;
+ if (k==3) mesh_partition = blub;
+ }
+ std::vector<int> simplexVertices(10);
+ switch (elm_type)
+ {
+ case 2: // 3-node triangle
+ simplexVertices.resize(3);
+ fscanf(file,"%d %d %d\n",&(simplexVertices[0]),&(simplexVertices[1]),&(simplexVertices[2]));
+ for (size_t i=0; i<simplexVertices.size(); i++)
+ if (renumber.find(simplexVertices[i])==renumber.end())
+ {
+ renumber[simplexVertices[i]] = number_of_real_vertices++;
+ factory.insertVertex(nodes[simplexVertices[i]]);
+ }
+ boundary_element_count++;
+ break;
+
+ case 4: // 4-node tetrahedron
+ simplexVertices.resize(4);
+ fscanf(file,"%d %d %d %d\n",&(simplexVertices[0]),&(simplexVertices[1]),&(simplexVertices[2]),&(simplexVertices[3]));
+ for (size_t i=0; i<simplexVertices.size(); i++)
+ if (renumber.find(simplexVertices[i])==renumber.end())
+ {
+ renumber[simplexVertices[i]] = number_of_real_vertices++;
+ factory.insertVertex(nodes[simplexVertices[i]]);
+ }
+ element_count++;
+ break;
+
+ case 9: // 6-node triangle
+ simplexVertices.resize(6);
+ fscanf(file,"%d %d %d %d %d %d\n",&(simplexVertices[0]),&(simplexVertices[1]),&(simplexVertices[2]),
+ &(simplexVertices[3]),&(simplexVertices[4]),&(simplexVertices[5]));
+ for (int i=0; i<3; i++) // insert only the first three !
+ if (renumber.find(simplexVertices[i])==renumber.end())
+ {
+ renumber[simplexVertices[i]] = number_of_real_vertices++;
+ factory.insertVertex(nodes[simplexVertices[i]]);
+ }
+ boundary_element_count++;
+ break;
+
+ case 11: // 10-node tetrahedron
+ simplexVertices.resize(10);
+ fscanf(file,"%d %d %d %d %d %d %d %d %d %d\n",&(simplexVertices[0]),&(simplexVertices[1]),&(simplexVertices[2]),
+ &(simplexVertices[3]),&(simplexVertices[4]),&(simplexVertices[5]),
+ &(simplexVertices[6]),&(simplexVertices[7]),&(simplexVertices[8]),&(simplexVertices[9]));
+ for (int i=0; i<4; i++) // insert only the first four !
+ if (renumber.find(simplexVertices[i])==renumber.end())
+ {
+ renumber[simplexVertices[i]] = number_of_real_vertices++;
+ factory.insertVertex(nodes[simplexVertices[i]]);
+ }
+ element_count++;
+ break;
+
+ default:
+ fgets(buf,512,file); // skip rest of line if no tetrahedron
+ }
+ }
+ if (verbose) std::cout << "number of real vertices = " << number_of_real_vertices << std::endl;
+ if (verbose) std::cout << "number of boundary elements = " << boundary_element_count << std::endl;
+ if (verbose) std::cout << "number of elements = " << element_count << std::endl;
+ fscanf(file,"%s\n",buf);
+ if (strcmp(buf,"$EndElements")!=0)
+ DUNE_THROW(Dune::IOError, "expected $EndElements");
+ boundary_id_to_physical_entity.resize(boundary_element_count);
+ element_index_to_physical_entity.resize(element_count);
+
+ //==============================================
+ // Pass 2: Insert boundary segments and elements
+ //==============================================
+
+ // go to beginning of file
+ rewind(file);
+ boundary_element_count = 0;
+ element_count = 0;
+
+ // process header
+ fscanf(file,"%s\n",buf);
+ if (strcmp(buf,"$MeshFormat")!=0)
+ DUNE_THROW(Dune::IOError, "expected $MeshFormat in first line");
+ fscanf(file,"%d %d %d\n",&version_number,&file_type,&data_size);
+ if (version_number!=2)
+ DUNE_THROW(Dune::IOError, "can only read version_number==2");
+ fscanf(file,"%s\n",buf);
+ if (strcmp(buf,"$EndMeshFormat")!=0)
+ DUNE_THROW(Dune::IOError, "expected $EndMeshFormat");
+
+ // node section
+ fscanf(file,"%s\n",buf);
+ if (strcmp(buf,"$Nodes")!=0)
+ DUNE_THROW(Dune::IOError, "expected $Nodes");
+ fscanf(file,"%d\n",&number_of_nodes);
+ for (int i=1; i<=number_of_nodes; i++)
+ {
+ int id;
+ double x,y,z;
+ fscanf(file,"%d %lg %lg %lg\n",&id,&x,&y,&z);
+ if (id!=i)
+ DUNE_THROW(Dune::IOError, "expected id " << i);
+ }
+ fscanf(file,"%s\n",buf);
+ if (strcmp(buf,"$EndNodes")!=0)
+ DUNE_THROW(Dune::IOError, "expected $EndNodes");
+
+ // element section
+ fscanf(file,"%s\n",buf);
+ if (strcmp(buf,"$Elements")!=0)
+ DUNE_THROW(Dune::IOError, "expected $Elements");
+ fscanf(file,"%d\n",&number_of_elements);
+
+ for (int i=1; i<=number_of_elements; i++)
+ {
+ int id, elm_type, number_of_tags;
+ fscanf(file,"%d %d %d",&id,&elm_type,&number_of_tags);
+ int physical_entity, elementary_entity, mesh_partition;
+ for (int k=1; k<=number_of_tags; k++)
+ {
+ int blub;
+ fscanf(file,"%d",&blub);
+ if (k==1) physical_entity = blub;
+ if (k==2) elementary_entity = blub;
+ if (k==3) mesh_partition = blub;
+ }
+ std::vector<int> simplexVertices(10);
+ std::vector<unsigned int> vertices(10);
+ switch (elm_type)
+ {
+ case 2: // 3-node triangle
+
+ // Read the vertices but don't do anything with them. Linear boundary segments
+ // are the default anyways.
+ simplexVertices.resize(3);
+ fscanf(file,"%d %d %d\n",&(simplexVertices[0]),&(simplexVertices[1]),&(simplexVertices[2]));
+ vertices.resize(3);
+ for (int i=0; i<3; i++)
+ vertices[i] = renumber[simplexVertices[i]]; // renumber vertices
+
+ factory.insertBoundarySegment(vertices,new GmeshReaderLinearBoundarySegment<3>(nodes[simplexVertices[0]],nodes[simplexVertices[1]],nodes[simplexVertices[2]]));
+ boundary_id_to_physical_entity[boundary_element_count] = physical_entity;
+ boundary_element_count++;
+ break;
+
+ case 9: // 6-node triangle
+ simplexVertices.resize(6);
+ fscanf(file,"%d %d %d %d %d %d\n",&(simplexVertices[0]),&(simplexVertices[1]),&(simplexVertices[2]),
+ &(simplexVertices[3]),&(simplexVertices[4]),&(simplexVertices[5]));
+ vertices.resize(3);
+ for (int i=0; i<3; i++)
+ vertices[i] = renumber[simplexVertices[i]]; // renumber vertices first three vertices
+
+ factory.insertBoundarySegment(vertices,new GmeshReaderQuadraticBoundarySegment<3>(nodes[simplexVertices[0]],nodes[simplexVertices[1]],nodes[simplexVertices[2]],
+ nodes[simplexVertices[3]],nodes[simplexVertices[4]],nodes[simplexVertices[5]]));
+ boundary_id_to_physical_entity[boundary_element_count] = physical_entity;
+ boundary_element_count++;
+ break;
+
+ case 4: // 4-node tetrahedron
+ simplexVertices.resize(4);
+ fscanf(file,"%d %d %d %d\n",&(simplexVertices[0]),&(simplexVertices[1]),&(simplexVertices[2]),&(simplexVertices[3]));
+ vertices.resize(4);
+ for (int i=0; i<4; i++)
+ vertices[i] = renumber[simplexVertices[i]]; // renumber vertices
+ factory.insertElement(Dune::GeometryType(Dune::GeometryType::simplex,dim),vertices);
+ element_index_to_physical_entity[element_count] = physical_entity;
+ element_count++;
+ break;
+
+ case 11: // 10-node tetrahedron
+ simplexVertices.resize(10);
+ fscanf(file,"%d %d %d %d %d %d %d %d %d %d\n",&(simplexVertices[0]),&(simplexVertices[1]),&(simplexVertices[2]),
+ &(simplexVertices[3]),&(simplexVertices[4]),&(simplexVertices[5]),
+ &(simplexVertices[6]),&(simplexVertices[7]),&(simplexVertices[8]),&(simplexVertices[9]));
+ vertices.resize(4);
+ for (int i=0; i<4; i++)
+ vertices[i] = renumber[simplexVertices[i]]; // renumber vertices
+ factory.insertElement(Dune::GeometryType(Dune::GeometryType::simplex,dim),vertices);
+ element_index_to_physical_entity[element_count] = physical_entity;
+ element_count++;
+ break;
+
+ default:
+ fgets(buf,512,file); // skip rest of line if no tetrahedron
+ }
+ }
+ fscanf(file,"%s\n",buf);
+ if (strcmp(buf,"$EndElements")!=0)
+ DUNE_THROW(Dune::IOError, "expected $EndElements");
+
+ fclose(file);
+
+ return factory.createGrid();
+ }
+ };
+
+ /**
+ \ingroup Gmesh
+
+ \brief Read Gmesh mesh file
+
+ Read a .msh file generated using Gmesh and construct a grid using the grid factory interface.
+ */
+ template<typename GridType>
+ class GmeshReader
+ {
+ public:
+ /** \todo doc me */
+ static GridType* read (const std::string& fileName, bool verbose = true)
+ {
+ // make a grid factory
+ Dune::GridFactory<GridType> factory;
+
+ std::vector<int> boundary_id_to_physical_entity;
+ std::vector<int> element_index_to_physical_entity;
+
+ return GmeshReaderImp<GridType,GridType::dimension>::read(factory,fileName,verbose,
+ boundary_id_to_physical_entity,
+ element_index_to_physical_entity);
+ }
+
+ /** \todo doc me */
+ template<typename T>
+ static GridType* read (const std::string& fileName, std::vector<T>& boundary_id_to_physical_entity,
+ std::vector<T>& element_index_to_physical_entity, bool verbose = true)
+ {
+ // make a grid factory
+ Dune::GridFactory<GridType> factory;
+
+ return GmeshReaderImp<GridType,GridType::dimension>::read(factory,fileName,verbose,
+ boundary_id_to_physical_entity,
+ element_index_to_physical_entity);
+ }
+
+ /** \todo doc me */
+ static GridType* read (GridType& grid, const std::string& fileName,
+ bool verbose = true)
+ {
+ // make a grid factory
+ Dune::GridFactory<GridType> factory(&grid);
+
+ // store dummy
+ std::vector<int> boundary_id_to_physical_entity;
+ std::vector<int> element_index_to_physical_entity;
+
+ return GmeshReaderImp<GridType,GridType::dimension>::read(factory,fileName,verbose,
+ boundary_id_to_physical_entity,
+ element_index_to_physical_entity);
+ }
+
+ /** \todo doc me */
+ template<typename T>
+ static GridType* read (GridType& grid, const std::string& fileName,
+ std::vector<T>& boundary_id_to_physical_entity,
+ std::vector<T>& element_index_to_physical_entity,
+ bool verbose = true)
+ {
+ // make a grid factory
+ Dune::GridFactory<GridType> factory(&grid);
+
+ return GmeshReaderImp<GridType,GridType::dimension>::read(factory,fileName,verbose,
+ boundary_id_to_physical_entity,
+ element_index_to_physical_entity);
+ }
+ };
+
+ /** \} */
+
+} // namespace Dune
+
+#endif
Deleted: trunk/grid/io/file/gmshreader.hh
===================================================================
--- trunk/grid/io/file/gmshreader.hh 2009-06-04 19:19:30 UTC (rev 5231)
+++ trunk/grid/io/file/gmshreader.hh 2009-06-04 19:20:25 UTC (rev 5232)
@@ -1,887 +0,0 @@
-// -*- tab-width: 4; indent-tabs-mode: nil -*-
-#ifndef DUNE_GMESHREADER_HH
-#define DUNE_GMESHREADER_HH
-
-#include<iostream>
-#include<fstream>
-#include<string>
-#include<vector>
-#include<map>
-#include<stdio.h>
-
-#include<dune/common/geometrytype.hh>
-#include<dune/common/exceptions.hh>
-#include<dune/common/fvector.hh>
-#include<dune/grid/common/gridfactory.hh>
-#include<dune/grid/common/boundarysegment.hh>
-
-namespace Dune {
-
- /**
- \ingroup Gmesh
- \{
- */
-
- //! Options for read operation
- struct GmeshReaderOptions
- {
- enum GeometryOrder {
- /** @brief edges are straight lines. */
- firstOrder,
- /** @brief quadratic boundary approximation. */
- secondOrder
- };
- };
-
- // arbitrary dimension, implementation is in specialization
- template<int dimension>
- class GmeshReaderLinearBoundarySegment
- {
- };
-
-
- // linear boundary segments in 1d
- template<>
- class GmeshReaderLinearBoundarySegment<2> : public Dune::BoundarySegment<2>
- {
- public:
- GmeshReaderLinearBoundarySegment (Dune::FieldVector<double,2> p0_, Dune::FieldVector<double,2> p1_)
- : p0(p0_), p1(p1_)
- {}
-
- virtual Dune::FieldVector<double,2> operator() (const Dune::FieldVector<double,1>& local) const
- {
- Dune::FieldVector<double,2> y;
- y = 0.0;
- y.axpy(1.0-local[0],p0);
- y.axpy(local[0],p1);
- return y;
- }
-
- private:
- Dune::FieldVector<double,2> p0,p1;
- };
-
-
-
- // arbitrary dimension, implementation is in specialization
- template<int dimension>
- class GmeshReaderQuadraticBoundarySegment
- {
- };
-
- // quadratic boundary segments in 1d
- /*
- Note the points
-
- (0) (alpha) (1)
-
- are mapped to the points in global coordinates
-
- p0 p2 p1
-
- alpha is determined automatically from the given points.
- */
- template <>
- class GmeshReaderQuadraticBoundarySegment<2> : public Dune::BoundarySegment<2>
- {
- public:
- GmeshReaderQuadraticBoundarySegment (Dune::FieldVector<double,2> p0_, Dune::FieldVector<double,2> p1_,
- Dune::FieldVector<double,2> p2_)
- : p0(p0_), p1(p1_), p2(p2_)
- {
- Dune::FieldVector<double,2> d1(p1);
- d1 -= p0;
- Dune::FieldVector<double,2> d2(p2);
- d2 -= p1;
-
- alpha=d1.two_norm()/(d1.two_norm()+d2.two_norm());
- if (alpha<1E-6 || alpha>1-1E-6)
- DUNE_THROW(Dune::IOError, "ration in quadratic boundary segment bad");
- }
-
- virtual Dune::FieldVector<double,2> operator() (const Dune::FieldVector<double,1>& local) const
- {
- Dune::FieldVector<double,2> y;
- y = 0.0;
- y.axpy((local[0]-alpha)*(local[0]-1.0)/alpha,p0);
- y.axpy(local[0]*(local[0]-1.0)/(alpha*(alpha-1.0)),p1);
- y.axpy(local[0]*(local[0]-alpha)/(1.0-alpha),p2);
- return y;
- }
-
- private:
- Dune::FieldVector<double,2> p0,p1,p2;
- double alpha;
- };
-
-
-
- // arbitrary dimension, implementation is in specialization
- template<typename GridType, int dimension>
- class GmeshReaderImp
- {
- };
-
- template<typename GridType>
- class GmeshReaderImp<GridType,2>
- {
- public:
-
- template<typename T>
- static GridType* read(Dune::GridFactory<GridType>& factory, const std::string& fileName,
- bool verbose, std::vector<T>& boundary_id_to_physical_entity,
- std::vector<T>& element_index_to_physical_entity)
- {
- // the grid dimension
- const int dim = GridType::dimension;
-
- // open file name, we use C I/O
- FILE* file = fopen(fileName.c_str(),"r");
- if (file==0)
- DUNE_THROW(Dune::IOError, "Could not open " << fileName);
-
- // a read buffer
- char buf[512];
-
- //=========================================
- // Pass 1: Read vertices into vector
- // Check vertices that are needed
- // Renumber needed vertices
- //=========================================
-
- int boundary_element_count = 0;
- int element_count = 0;
-
- // process header
- fscanf(file,"%s\n",buf);
- if (strcmp(buf,"$MeshFormat")!=0)
- DUNE_THROW(Dune::IOError, "expected $MeshFormat in first line");
- int version_number, file_type, data_size;
- fscanf(file,"%d %d %d\n",&version_number,&file_type,&data_size);
- if (version_number!=2)
- DUNE_THROW(Dune::IOError, "can only read version_number==2");
- if (verbose) std::cout << "version 2 Gmesh file detected" << std::endl;
- fscanf(file,"%s\n",buf);
- if (strcmp(buf,"$EndMeshFormat")!=0)
- DUNE_THROW(Dune::IOError, "expected $EndMeshFormat");
-
- // node section
- fscanf(file,"%s\n",buf);
- if (strcmp(buf,"$Nodes")!=0)
- DUNE_THROW(Dune::IOError, "expected $Nodes");
- int number_of_nodes;
- fscanf(file,"%d\n",&number_of_nodes);
- if (verbose) std::cout << "file contains " << number_of_nodes << " nodes" << std::endl;
- std::vector<Dune::FieldVector<double,dim> > nodes(number_of_nodes+1); // store positions
- for (int i=1; i<=number_of_nodes; i++)
- {
- int id;
- double x,y,z;
- fscanf(file,"%d %lg %lg %lg\n",&id,&x,&y,&z);
- // if (verbose) std::cout << id << " " << x << " " << y << " " << z << std::endl;
- if (id!=i)
- DUNE_THROW(Dune::IOError, "expected id " << i);
-
- Dune::FieldVector<double,dim> position;
- position[0] = x; position[1] = y;
- nodes[i] = position; // just store node position
- }
- fscanf(file,"%s\n",buf);
- if (strcmp(buf,"$EndNodes")!=0)
- DUNE_THROW(Dune::IOError, "expected $EndNodes");
-
- // element section
- fscanf(file,"%s\n",buf);
- if (strcmp(buf,"$Elements")!=0)
- DUNE_THROW(Dune::IOError, "expected $Elements");
- int number_of_elements;
- fscanf(file,"%d\n",&number_of_elements);
- if (verbose) std::cout << "file contains " << number_of_elements << " elements" << std::endl;
- unsigned int number_of_real_vertices=0; // count number of vertices that are really needed
- std::map<int,unsigned int> renumber;
- for (int i=1; i<=number_of_elements; i++)
- {
- int id, elm_type, number_of_tags;
- fscanf(file,"%d %d %d",&id,&elm_type,&number_of_tags);
- int physical_entity, elementary_entity, mesh_partition;
- for (int k=1; k<=number_of_tags; k++)
- {
- int blub;
- fscanf(file,"%d",&blub);
- if (k==1) physical_entity = blub;
- if (k==2) elementary_entity = blub;
- if (k==3) mesh_partition = blub;
- }
- std::vector<int> simplexVertices(6);
- switch (elm_type)
- {
- case 1: // 2-node line
- simplexVertices.resize(2);
- fscanf(file,"%d %d\n",&(simplexVertices[0]),&(simplexVertices[1]));
- for (int i=0; i<2; i++)
- if (renumber.find(simplexVertices[i])==renumber.end())
- {
- // std::cout << "real vertex " << number_of_real_vertices << " at " << nodes[simplexVertices[i]] << " old number was " << simplexVertices[i] << std::endl;
- renumber[simplexVertices[i]] = number_of_real_vertices++;
- factory.insertVertex(nodes[simplexVertices[i]]);
- }
- boundary_element_count++;
- break;
- case 8: // 3-node line
- simplexVertices.resize(3);
- fscanf(file,"%d %d %d\n",&(simplexVertices[0]),&(simplexVertices[1]),&(simplexVertices[2]));
- for (int i=0; i<2; i++)
- if (renumber.find(simplexVertices[i])==renumber.end())
- {
- renumber[simplexVertices[i]] = number_of_real_vertices++;
- factory.insertVertex(nodes[simplexVertices[i]]);
- }
- boundary_element_count++;
- break;
- case 2: // 3-node triangle
- simplexVertices.resize(3);
- fscanf(file,"%d %d %d\n",&(simplexVertices[0]),&(simplexVertices[1]),&(simplexVertices[2]));
- for (int i=0; i<3; i++)
- if (renumber.find(simplexVertices[i])==renumber.end())
- {
- renumber[simplexVertices[i]] = number_of_real_vertices++;
- factory.insertVertex(nodes[simplexVertices[i]]);
- }
- element_count++;
- break;
- case 9: // 6-node triangle
- simplexVertices.resize(6);
- fscanf(file,"%d %d %d %d %d %d\n",&(simplexVertices[0]),&(simplexVertices[1]),&(simplexVertices[2]),
- &(simplexVertices[3]),&(simplexVertices[4]),&(simplexVertices[5]));
- for (int i=0; i<3; i++)
- if (renumber.find(simplexVertices[i])==renumber.end())
- {
- renumber[simplexVertices[i]] = number_of_real_vertices++;
- factory.insertVertex(nodes[simplexVertices[i]]);
- }
- element_count++;
- break;
- default:
- fgets(buf,512,file); // skip rest of line if no tetrahedron
- }
- }
- if (verbose) std::cout << "number of real vertices = " << number_of_real_vertices << std::endl;
- if (verbose) std::cout << "number of boundary elements = " << boundary_element_count << std::endl;
- if (verbose) std::cout << "number of elements = " << element_count << std::endl;
- fscanf(file,"%s\n",buf);
- if (strcmp(buf,"$EndElements")!=0)
- DUNE_THROW(Dune::IOError, "expected $EndElements");
- boundary_id_to_physical_entity.resize(boundary_element_count);
- element_index_to_physical_entity.resize(element_count);
-
- //==============================================
- // Pass 2: Insert boundary segments and elements
- //==============================================
-
- // go to beginning of file
- rewind(file);
- boundary_element_count = 0;
- element_count = 0;
-
- // process header
- fscanf(file,"%s\n",buf);
- if (strcmp(buf,"$MeshFormat")!=0)
- DUNE_THROW(Dune::IOError, "expected $MeshFormat in first line");
- fscanf(file,"%d %d %d\n",&version_number,&file_type,&data_size);
- if (version_number!=2)
- DUNE_THROW(Dune::IOError, "can only read version_number==2");
- fscanf(file,"%s\n",buf);
- if (strcmp(buf,"$EndMeshFormat")!=0)
- DUNE_THROW(Dune::IOError, "expected $EndMeshFormat");
-
- // node section
- fscanf(file,"%s\n",buf);
- if (strcmp(buf,"$Nodes")!=0)
- DUNE_THROW(Dune::IOError, "expected $Nodes");
- fscanf(file,"%d\n",&number_of_nodes);
- for (int i=1; i<=number_of_nodes; i++)
- {
- int id;
- double x,y,z;
- fscanf(file,"%d %lg %lg %lg\n",&id,&x,&y,&z);
- if (id!=i)
- DUNE_THROW(Dune::IOError, "expected id " << i);
- }
- fscanf(file,"%s\n",buf);
- if (strcmp(buf,"$EndNodes")!=0)
- DUNE_THROW(Dune::IOError, "expected $EndNodes");
-
- // element section
- fscanf(file,"%s\n",buf);
- if (strcmp(buf,"$Elements")!=0)
- DUNE_THROW(Dune::IOError, "expected $Elements");
- fscanf(file,"%d\n",&number_of_elements);
- for (int i=1; i<=number_of_elements; i++)
- {
- int id, elm_type, number_of_tags;
- fscanf(file,"%d %d %d",&id,&elm_type,&number_of_tags);
- int physical_entity, elementary_entity, mesh_partition;
- for (int k=1; k<=number_of_tags; k++)
- {
- int blub;
- fscanf(file,"%d",&blub);
- if (k==1) physical_entity = blub;
- if (k==2) elementary_entity = blub;
- if (k==3) mesh_partition = blub;
- }
- std::vector<int> simplexVertices(6);
- std::vector<unsigned int> vertices(3);
- switch (elm_type)
- {
- case 1: // 2-node line
-
- // read the vertices, but don't do anything with them. Linear boundary
- // segments are the default.
- simplexVertices.resize(2);
- fscanf(file,"%d %d\n",&(simplexVertices[0]),&(simplexVertices[1]));
- vertices.resize(2);
- for (int i=0; i<2; i++)
- vertices[i] = renumber[simplexVertices[i]]; // renumber vertices
- factory.insertBoundarySegment(vertices,new GmeshReaderLinearBoundarySegment<2>(nodes[simplexVertices[0]],nodes[simplexVertices[1]]));
- boundary_id_to_physical_entity[boundary_element_count] = physical_entity;
- boundary_element_count++;
- break;
- case 8: // 3-node line
- simplexVertices.resize(3);
- fscanf(file,"%d %d %d\n",&(simplexVertices[0]),&(simplexVertices[1]),&(simplexVertices[2]));
- vertices.resize(2);
- for (int i=0; i<2; i++)
- vertices[i] = renumber[simplexVertices[i]]; // renumber vertices
- factory.insertBoundarySegment(vertices,new GmeshReaderQuadraticBoundarySegment<2>(nodes[simplexVertices[0]],
- nodes[simplexVertices[2]],
- nodes[simplexVertices[1]]));
- boundary_id_to_physical_entity[boundary_element_count] = physical_entity;
- boundary_element_count++;
- break;
- case 2: // 3-node triangle
- simplexVertices.resize(3);
- fscanf(file,"%d %d %d\n",&(simplexVertices[0]),&(simplexVertices[1]),&(simplexVertices[2]));
- vertices.resize(3);
- for (int i=0; i<3; i++)
- vertices[i] = renumber[simplexVertices[i]]; // renumber vertices
- factory.insertElement(Dune::GeometryType(Dune::GeometryType::simplex,dim),vertices);
- element_index_to_physical_entity[element_count] = physical_entity;
- element_count++;
- break;
- case 9: // 6-node triangle
- simplexVertices.resize(6);
- fscanf(file,"%d %d %d %d %d %d\n",&(simplexVertices[0]),&(simplexVertices[1]),&(simplexVertices[2]),
- &(simplexVertices[3]),&(simplexVertices[4]),&(simplexVertices[5]));
- vertices.resize(3);
- for (int i=0; i<3; i++)
- vertices[i] = renumber[simplexVertices[i]]; // renumber vertices
- factory.insertElement(Dune::GeometryType(Dune::GeometryType::simplex,dim),vertices);
- element_index_to_physical_entity[element_count] = physical_entity;
- element_count++;
- break;
- default:
- fgets(buf,512,file); // skip rest of line if no tetrahedron
- }
- }
- fscanf(file,"%s\n",buf);
- if (strcmp(buf,"$EndElements")!=0)
- DUNE_THROW(Dune::IOError, "expected $EndElements");
-
- fclose(file);
-
- return factory.createGrid();
- }
- };
-
-
- // linear boundary segments in 2d
- template <>
- class GmeshReaderLinearBoundarySegment<3> : public Dune::BoundarySegment<3>
- {
- public:
- GmeshReaderLinearBoundarySegment (Dune::FieldVector<double,3> p0_, Dune::FieldVector<double,3> p1_,
- Dune::FieldVector<double,3> p2_)
- : p0(p0_), p1(p1_), p2(p2_)
- {
- // std::cout << "created boundary segment " << p0 << " | " << p1 << " | " << p2 << std::endl;
- }
-
- virtual Dune::FieldVector<double,3> operator() (const Dune::FieldVector<double,2>& local) const
- {
- Dune::FieldVector<double,3> y;
- y = 0.0;
- y.axpy(1.0-local[0]-local[1],p0);
- y.axpy(local[0],p1);
- y.axpy(local[1],p2);
- // std::cout << "eval boundary segment local=" << local << " y=" << y << std::endl;
- return y;
- }
-
- private:
- Dune::FieldVector<double,3> p0,p1,p2;
- };
-
-
- // quadratic boundary segments in 2d
- /* numbering of points corresponding to gmsh:
-
- 2
-
- 5 4
-
- 0 3 1
-
- Note: The vertices 3, 4, 5 are not necessarily at the edge midpoints but can
- be placed with parameters alpha, beta , gamma at the following positions
- in local coordinates:
-
-
- 2 = (0,1)
-
- 5 = (0,beta) 4 = (1-gamma/sqrt(2),gamma/sqrt(2))
-
- 0 = (0,0) 3 = (alpha,0) 1 = (1,0)
-
- The parameters alpha, beta, gamma are determined from the given vertices in
- global coordinates.
- */
- template <>
- class GmeshReaderQuadraticBoundarySegment<3> : public Dune::BoundarySegment<3>
- {
- public:
- GmeshReaderQuadraticBoundarySegment (Dune::FieldVector<double,3> p0_, Dune::FieldVector<double,3> p1_,
- Dune::FieldVector<double,3> p2_, Dune::FieldVector<double,3> p3_,
- Dune::FieldVector<double,3> p4_, Dune::FieldVector<double,3> p5_)
- : p0(p0_), p1(p1_), p2(p2_), p3(p3_), p4(p4_), p5(p5_)
- {
- Dune::FieldVector<double,3> d1,d2;
-
- d1 = p3; d1 -= p0;
- d2 = p1; d2 -= p3;
- alpha=d1.two_norm()/(d1.two_norm()+d2.two_norm());
- if (alpha<1E-6 || alpha>1-1E-6)
- DUNE_THROW(Dune::IOError, "alpha in quadratic boundary segment bad");
-
- d1 = p5; d1 -= p0;
- d2 = p2; d2 -= p5;
- beta=d1.two_norm()/(d1.two_norm()+d2.two_norm());
- if (beta<1E-6 || beta>1-1E-6)
- DUNE_THROW(Dune::IOError, "beta in quadratic boundary segment bad");
-
- d1 = p4; d1 -= p1;
- d2 = p2; d2 -= p4;
- gamma=d1.two_norm()/(d1.two_norm()+d2.two_norm());
- if (gamma<1E-6 || gamma>1-1E-6)
- DUNE_THROW(Dune::IOError, "gamma in quadratic boundary segment bad");
-
- sqrt2 = sqrt(2.0);
- }
-
- virtual Dune::FieldVector<double,3> operator() (const Dune::FieldVector<double,2>& local) const
- {
- Dune::FieldVector<double,3> y;
- y = 0.0;
- y.axpy(phi0(local),p0);
- y.axpy(phi1(local),p1);
- y.axpy(phi2(local),p2);
- y.axpy(phi3(local),p3);
- y.axpy(phi4(local),p4);
- y.axpy(phi5(local),p5);
- return y;
- }
-
- private:
- // The six Lagrange basis function on the reference element
- // for the points given above
-
- double phi0 (const Dune::FieldVector<double,2>& local) const
- {
- return (alpha*beta-beta*local[0]-alpha*local[1])*(1-local[0]-local[1])/(alpha*beta);
- }
- double phi3 (const Dune::FieldVector<double,2>& local) const
- {
- return local[0]*(1-local[0]-local[1])/(alpha*(1-alpha));
- }
- double phi1 (const Dune::FieldVector<double,2>& local) const
- {
- return local[0]*(gamma*local[0]-(sqrt2-gamma-sqrt2*alpha)*local[1]-alpha*gamma)/(gamma*(1-alpha));
- }
- double phi5 (const Dune::FieldVector<double,2>& local) const
- {
- return local[1]*(1-local[0]-local[1])/(beta*(1-beta));
- }
- double phi4 (const Dune::FieldVector<double,2>& local) const
- {
- return local[0]*local[1]/((1-gamma/sqrt2)*gamma/sqrt2);
- }
- double phi2 (const Dune::FieldVector<double,2>& local) const
- {
- return local[1]*(beta*(1-gamma/sqrt2)-local[0]*(beta-gamma/sqrt2)-local[1]*(1-gamma/sqrt2))/((1-gamma/sqrt2)*(beta-1));
- }
-
- Dune::FieldVector<double,3> p0,p1,p2,p3,p4,p5;
- double alpha,beta,gamma,sqrt2;
- };
-
-
-
- template<typename GridType>
- class GmeshReaderImp<GridType,3>
- {
- public:
-
- template<typename T>
- static GridType* read(Dune::GridFactory<GridType>& factory, const std::string& fileName,
- bool verbose, std::vector<T>& boundary_id_to_physical_entity,
- std::vector<T>& element_index_to_physical_entity)
- {
- // the grid dimension
- const int dim = GridType::dimension;
-
- // open file name, we use C I/O
- FILE* file = fopen(fileName.c_str(),"r");
- if (file==0)
- DUNE_THROW(Dune::IOError, "Could not open " << fileName);
-
- // a read buffer
- char buf[512];
-
- //=========================================
- // Pass 1: Read vertices into vector
- // Check vertices that are needed
- // Renumber needed vertices
- //=========================================
-
- int boundary_element_count = 0;
- int element_count = 0;
-
- // process header
- fscanf(file,"%s\n",buf);
- if (strcmp(buf,"$MeshFormat")!=0)
- DUNE_THROW(Dune::IOError, "expected $MeshFormat in first line");
- int version_number, file_type, data_size;
- fscanf(file,"%d %d %d\n",&version_number,&file_type,&data_size);
- if (version_number!=2)
- DUNE_THROW(Dune::IOError, "can only read version_number==2");
- if (verbose) std::cout << "version 2 Gmesh file detected" << std::endl;
- fscanf(file,"%s\n",buf);
- if (strcmp(buf,"$EndMeshFormat")!=0)
- DUNE_THROW(Dune::IOError, "expected $EndMeshFormat");
-
- // node section
- fscanf(file,"%s\n",buf);
- if (strcmp(buf,"$Nodes")!=0)
- DUNE_THROW(Dune::IOError, "expected $Nodes");
- int number_of_nodes;
- fscanf(file,"%d\n",&number_of_nodes);
- if (verbose) std::cout << "file contains " << number_of_nodes << " nodes" << std::endl;
- std::vector<Dune::FieldVector<double,dim> > nodes(number_of_nodes+1); // store positions
- for (int i=1; i<=number_of_nodes; i++) // note: ids are starting with 1!
- {
- int id;
- double x,y,z;
- fscanf(file,"%d %lg %lg %lg\n",&id,&x,&y,&z);
- // if (verbose) std::cout << id << " " << x << " " << y << " " << z << std::endl;
- if (id!=i)
- DUNE_THROW(Dune::IOError, "expected id " << i);
-
- Dune::FieldVector<double,dim> position;
- position[0] = x; position[1] = y; position[2] = z;
- nodes[i] = position; // just store node position
- }
- fscanf(file,"%s\n",buf);
- if (strcmp(buf,"$EndNodes")!=0)
- DUNE_THROW(Dune::IOError, "expected $EndNodes");
-
- // element section
- fscanf(file,"%s\n",buf);
- if (strcmp(buf,"$Elements")!=0)
- DUNE_THROW(Dune::IOError, "expected $Elements");
- int number_of_elements;
- fscanf(file,"%d\n",&number_of_elements);
- if (verbose) std::cout << "file contains " << number_of_elements << " elements" << std::endl;
- unsigned int number_of_real_vertices=0; // count number of vertices that are really needed
- std::map<int,unsigned int> renumber;
- for (int i=1; i<=number_of_elements; i++) // note: ids are starting with 1!
- {
- int id, elm_type, number_of_tags;
- fscanf(file,"%d %d %d",&id,&elm_type,&number_of_tags);
- int physical_entity, elementary_entity, mesh_partition;
- for (int k=1; k<=number_of_tags; k++)
- {
- int blub;
- fscanf(file,"%d",&blub);
- if (k==1) physical_entity = blub;
- if (k==2) elementary_entity = blub;
- if (k==3) mesh_partition = blub;
- }
- std::vector<int> simplexVertices(10);
- switch (elm_type)
- {
- case 2: // 3-node triangle
- simplexVertices.resize(3);
- fscanf(file,"%d %d %d\n",&(simplexVertices[0]),&(simplexVertices[1]),&(simplexVertices[2]));
- for (size_t i=0; i<simplexVertices.size(); i++)
- if (renumber.find(simplexVertices[i])==renumber.end())
- {
- renumber[simplexVertices[i]] = number_of_real_vertices++;
- factory.insertVertex(nodes[simplexVertices[i]]);
- }
- boundary_element_count++;
- break;
-
- case 4: // 4-node tetrahedron
- simplexVertices.resize(4);
- fscanf(file,"%d %d %d %d\n",&(simplexVertices[0]),&(simplexVertices[1]),&(simplexVertices[2]),&(simplexVertices[3]));
- for (size_t i=0; i<simplexVertices.size(); i++)
- if (renumber.find(simplexVertices[i])==renumber.end())
- {
- renumber[simplexVertices[i]] = number_of_real_vertices++;
- factory.insertVertex(nodes[simplexVertices[i]]);
- }
- element_count++;
- break;
-
- case 9: // 6-node triangle
- simplexVertices.resize(6);
- fscanf(file,"%d %d %d %d %d %d\n",&(simplexVertices[0]),&(simplexVertices[1]),&(simplexVertices[2]),
- &(simplexVertices[3]),&(simplexVertices[4]),&(simplexVertices[5]));
- for (int i=0; i<3; i++) // insert only the first three !
- if (renumber.find(simplexVertices[i])==renumber.end())
- {
- renumber[simplexVertices[i]] = number_of_real_vertices++;
- factory.insertVertex(nodes[simplexVertices[i]]);
- }
- boundary_element_count++;
- break;
-
- case 11: // 10-node tetrahedron
- simplexVertices.resize(10);
- fscanf(file,"%d %d %d %d %d %d %d %d %d %d\n",&(simplexVertices[0]),&(simplexVertices[1]),&(simplexVertices[2]),
- &(simplexVertices[3]),&(simplexVertices[4]),&(simplexVertices[5]),
- &(simplexVertices[6]),&(simplexVertices[7]),&(simplexVertices[8]),&(simplexVertices[9]));
- for (int i=0; i<4; i++) // insert only the first four !
- if (renumber.find(simplexVertices[i])==renumber.end())
- {
- renumber[simplexVertices[i]] = number_of_real_vertices++;
- factory.insertVertex(nodes[simplexVertices[i]]);
- }
- element_count++;
- break;
-
- default:
- fgets(buf,512,file); // skip rest of line if no tetrahedron
- }
- }
- if (verbose) std::cout << "number of real vertices = " << number_of_real_vertices << std::endl;
- if (verbose) std::cout << "number of boundary elements = " << boundary_element_count << std::endl;
- if (verbose) std::cout << "number of elements = " << element_count << std::endl;
- fscanf(file,"%s\n",buf);
- if (strcmp(buf,"$EndElements")!=0)
- DUNE_THROW(Dune::IOError, "expected $EndElements");
- boundary_id_to_physical_entity.resize(boundary_element_count);
- element_index_to_physical_entity.resize(element_count);
-
- //==============================================
- // Pass 2: Insert boundary segments and elements
- //==============================================
-
- // go to beginning of file
- rewind(file);
- boundary_element_count = 0;
- element_count = 0;
-
- // process header
- fscanf(file,"%s\n",buf);
- if (strcmp(buf,"$MeshFormat")!=0)
- DUNE_THROW(Dune::IOError, "expected $MeshFormat in first line");
- fscanf(file,"%d %d %d\n",&version_number,&file_type,&data_size);
- if (version_number!=2)
- DUNE_THROW(Dune::IOError, "can only read version_number==2");
- fscanf(file,"%s\n",buf);
- if (strcmp(buf,"$EndMeshFormat")!=0)
- DUNE_THROW(Dune::IOError, "expected $EndMeshFormat");
-
- // node section
- fscanf(file,"%s\n",buf);
- if (strcmp(buf,"$Nodes")!=0)
- DUNE_THROW(Dune::IOError, "expected $Nodes");
- fscanf(file,"%d\n",&number_of_nodes);
- for (int i=1; i<=number_of_nodes; i++)
- {
- int id;
- double x,y,z;
- fscanf(file,"%d %lg %lg %lg\n",&id,&x,&y,&z);
- if (id!=i)
- DUNE_THROW(Dune::IOError, "expected id " << i);
- }
- fscanf(file,"%s\n",buf);
- if (strcmp(buf,"$EndNodes")!=0)
- DUNE_THROW(Dune::IOError, "expected $EndNodes");
-
- // element section
- fscanf(file,"%s\n",buf);
- if (strcmp(buf,"$Elements")!=0)
- DUNE_THROW(Dune::IOError, "expected $Elements");
- fscanf(file,"%d\n",&number_of_elements);
-
- for (int i=1; i<=number_of_elements; i++)
- {
- int id, elm_type, number_of_tags;
- fscanf(file,"%d %d %d",&id,&elm_type,&number_of_tags);
- int physical_entity, elementary_entity, mesh_partition;
- for (int k=1; k<=number_of_tags; k++)
- {
- int blub;
- fscanf(file,"%d",&blub);
- if (k==1) physical_entity = blub;
- if (k==2) elementary_entity = blub;
- if (k==3) mesh_partition = blub;
- }
- std::vector<int> simplexVertices(10);
- std::vector<unsigned int> vertices(10);
- switch (elm_type)
- {
- case 2: // 3-node triangle
-
- // Read the vertices but don't do anything with them. Linear boundary segments
- // are the default anyways.
- simplexVertices.resize(3);
- fscanf(file,"%d %d %d\n",&(simplexVertices[0]),&(simplexVertices[1]),&(simplexVertices[2]));
- vertices.resize(3);
- for (int i=0; i<3; i++)
- vertices[i] = renumber[simplexVertices[i]]; // renumber vertices
-
- factory.insertBoundarySegment(vertices,new GmeshReaderLinearBoundarySegment<3>(nodes[simplexVertices[0]],nodes[simplexVertices[1]],nodes[simplexVertices[2]]));
- boundary_id_to_physical_entity[boundary_element_count] = physical_entity;
- boundary_element_count++;
- break;
-
- case 9: // 6-node triangle
- simplexVertices.resize(6);
- fscanf(file,"%d %d %d %d %d %d\n",&(simplexVertices[0]),&(simplexVertices[1]),&(simplexVertices[2]),
- &(simplexVertices[3]),&(simplexVertices[4]),&(simplexVertices[5]));
- vertices.resize(3);
- for (int i=0; i<3; i++)
- vertices[i] = renumber[simplexVertices[i]]; // renumber vertices first three vertices
-
- factory.insertBoundarySegment(vertices,new GmeshReaderQuadraticBoundarySegment<3>(nodes[simplexVertices[0]],nodes[simplexVertices[1]],nodes[simplexVertices[2]],
- nodes[simplexVertices[3]],nodes[simplexVertices[4]],nodes[simplexVertices[5]]));
- boundary_id_to_physical_entity[boundary_element_count] = physical_entity;
- boundary_element_count++;
- break;
-
- case 4: // 4-node tetrahedron
- simplexVertices.resize(4);
- fscanf(file,"%d %d %d %d\n",&(simplexVertices[0]),&(simplexVertices[1]),&(simplexVertices[2]),&(simplexVertices[3]));
- vertices.resize(4);
- for (int i=0; i<4; i++)
- vertices[i] = renumber[simplexVertices[i]]; // renumber vertices
- factory.insertElement(Dune::GeometryType(Dune::GeometryType::simplex,dim),vertices);
- element_index_to_physical_entity[element_count] = physical_entity;
- element_count++;
- break;
-
- case 11: // 10-node tetrahedron
- simplexVertices.resize(10);
- fscanf(file,"%d %d %d %d %d %d %d %d %d %d\n",&(simplexVertices[0]),&(simplexVertices[1]),&(simplexVertices[2]),
- &(simplexVertices[3]),&(simplexVertices[4]),&(simplexVertices[5]),
- &(simplexVertices[6]),&(simplexVertices[7]),&(simplexVertices[8]),&(simplexVertices[9]));
- vertices.resize(4);
- for (int i=0; i<4; i++)
- vertices[i] = renumber[simplexVertices[i]]; // renumber vertices
- factory.insertElement(Dune::GeometryType(Dune::GeometryType::simplex,dim),vertices);
- element_index_to_physical_entity[element_count] = physical_entity;
- element_count++;
- break;
-
- default:
- fgets(buf,512,file); // skip rest of line if no tetrahedron
- }
- }
- fscanf(file,"%s\n",buf);
- if (strcmp(buf,"$EndElements")!=0)
- DUNE_THROW(Dune::IOError, "expected $EndElements");
-
- fclose(file);
-
- return factory.createGrid();
- }
- };
-
- /**
- \ingroup Gmesh
-
- \brief Read Gmesh mesh file
-
- Read a .msh file generated using Gmesh and construct a grid using the grid factory interface.
- */
- template<typename GridType>
- class GmeshReader
- {
- public:
- /** \todo doc me */
- static GridType* read (const std::string& fileName, bool verbose = true)
- {
- // make a grid factory
- Dune::GridFactory<GridType> factory;
-
- std::vector<int> boundary_id_to_physical_entity;
- std::vector<int> element_index_to_physical_entity;
-
- return GmeshReaderImp<GridType,GridType::dimension>::read(factory,fileName,verbose,
- boundary_id_to_physical_entity,
- element_index_to_physical_entity);
- }
-
- /** \todo doc me */
- template<typename T>
- static GridType* read (const std::string& fileName, std::vector<T>& boundary_id_to_physical_entity,
- std::vector<T>& element_index_to_physical_entity, bool verbose = true)
- {
- // make a grid factory
- Dune::GridFactory<GridType> factory;
-
- return GmeshReaderImp<GridType,GridType::dimension>::read(factory,fileName,verbose,
- boundary_id_to_physical_entity,
- element_index_to_physical_entity);
- }
-
- /** \todo doc me */
- static GridType* read (GridType& grid, const std::string& fileName,
- bool verbose = true)
- {
- // make a grid factory
- Dune::GridFactory<GridType> factory(&grid);
-
- // store dummy
- std::vector<int> boundary_id_to_physical_entity;
- std::vector<int> element_index_to_physical_entity;
-
- return GmeshReaderImp<GridType,GridType::dimension>::read(factory,fileName,verbose,
- boundary_id_to_physical_entity,
- element_index_to_physical_entity);
- }
-
- /** \todo doc me */
- template<typename T>
- static GridType* read (GridType& grid, const std::string& fileName,
- std::vector<T>& boundary_id_to_physical_entity,
- std::vector<T>& element_index_to_physical_entity,
- bool verbose = true)
- {
- // make a grid factory
- Dune::GridFactory<GridType> factory(&grid);
-
- return GmeshReaderImp<GridType,GridType::dimension>::read(factory,fileName,verbose,
- boundary_id_to_physical_entity,
- element_index_to_physical_entity);
- }
- };
-
- /** \} */
-
-} // namespace Dune
-
-#endif
_______________________________________________
Dune-Commit mailing list
Dune-Commit at dune-project.org
http://lists.dune-project.org/mailman/listinfo/dune-commit
--
************************************************************************
* Oliver Sander ** email: sander at mi.fu-berlin.de *
* Freie Universität Berlin ** phone: + 49 (30) 838 75348 *
* Institut für Mathematik ** URL : page.mi.fu-berlin.de/~sander *
* Arnimallee 6 ** -------------------------------------*
* 14195 Berlin, Germany ** Member of MATHEON (www.matheon.de) *
************************************************************************
More information about the Dune
mailing list