[Dune] gmsh vs. gmesh
Peter Bastian
peter.bastian at iwr.uni-heidelberg.de
Fri Jun 5 15:18:46 CEST 2009
Oliver,
damn, you are right. I thought I had it wrong all the time. Then
obviously you changed it and I forgot that you changed it. I will change
it one more time ...
Peter
Am Freitag, den 05.06.2009, 13:32 +0200 schrieb Oliver Sander:
> 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
>
>
--
Peter Bastian
Interdisziplinäres Zentrum für Wissenschaftliches Rechnen
Im Neuenheimer Feld 368
D-69120 Heidelberg
Telefon +49 (6221) 54-8261/8249
Telefax +49 (6221) 54-8884
Email: peter.bastian at iwr.uni-heidelberg.de
More information about the Dune
mailing list