dune-grid  2.2.0
gmshreader.hh
Go to the documentation of this file.
00001 // -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
00002 // vi: set et ts=4 sw=4 sts=4:
00003 
00004 #ifndef DUNE_GMSHREADER_HH
00005 #define DUNE_GMSHREADER_HH
00006 
00007 #include <cstdarg>
00008 #include <cstdio>
00009 #include <cstring>
00010 #include <fstream>
00011 #include <iostream>
00012 #include <map>
00013 #include <string>
00014 #include <vector>
00015 
00016 #include <stdio.h>
00017 
00018 #include <dune/common/exceptions.hh>
00019 #include <dune/common/fvector.hh>
00020 
00021 #include <dune/geometry/type.hh>
00022 
00023 #include <dune/grid/common/boundarysegment.hh>
00024 #include <dune/grid/common/gridfactory.hh>
00025 
00026 namespace Dune
00027 {
00028 
00034 
00035     struct GmshReaderOptions
00036     {
00037         enum GeometryOrder {
00039             firstOrder,
00041             secondOrder 
00042         };
00043     };
00044 
00045     namespace {
00046     
00047         // arbitrary dimension, implementation is in specialization
00048         template< int dimension, int dimWorld = dimension >
00049         class GmshReaderQuadraticBoundarySegment
00050         {
00051         };
00052     
00053         // quadratic boundary segments in 1d
00054         /*
00055           Note the points
00056 
00057           (0)   (alpha)   (1)
00058 
00059           are mapped to the points in global coordinates
00060 
00061           p0 p2 p1
00062 
00063           alpha is determined automatically from the given points.
00064         */
00065         template< int dimWorld >
00066         struct GmshReaderQuadraticBoundarySegment< 2, dimWorld >
00067             : public Dune::BoundarySegment< 2, dimWorld >
00068         {
00069             typedef Dune::FieldVector< double, dimWorld > GlobalVector;
00070 
00071             GmshReaderQuadraticBoundarySegment ( const GlobalVector &p0_, const GlobalVector &p1_, const GlobalVector &p2_)
00072                 : p0(p0_), p1(p1_), p2(p2_)
00073             {
00074                 GlobalVector d1 = p1;
00075                 d1 -= p0;
00076                 GlobalVector d2 = p2;
00077                 d2 -= p1;
00078 
00079                 alpha=d1.two_norm()/(d1.two_norm()+d2.two_norm());
00080                 if (alpha<1E-6 || alpha>1-1E-6)
00081                     DUNE_THROW(Dune::IOError, "ration in quadratic boundary segment bad");
00082             }
00083 
00084             virtual GlobalVector operator() ( const Dune::FieldVector<double,1> &local ) const
00085             {
00086                 GlobalVector y;
00087                 y = 0.0;
00088                 y.axpy((local[0]-alpha)*(local[0]-1.0)/alpha,p0);
00089                 y.axpy(local[0]*(local[0]-1.0)/(alpha*(alpha-1.0)),p1);
00090                 y.axpy(local[0]*(local[0]-alpha)/(1.0-alpha),p2);
00091                 return y;
00092             }
00093 
00094         private:
00095             GlobalVector p0,p1,p2;
00096             double alpha;
00097         };
00098 
00099 
00100         // quadratic boundary segments in 2d
00101         /* numbering of points corresponding to gmsh:
00102 
00103            2
00104 
00105            5  4
00106 
00107            0  3  1
00108 
00109            Note: The vertices 3, 4, 5 are not necessarily at the edge midpoints but can
00110            be placed with parameters alpha, beta , gamma at the following positions
00111            in local coordinates:
00112 
00113 
00114            2 = (0,1)
00115 
00116            5 = (0,beta)   4 = (1-gamma/sqrt(2),gamma/sqrt(2))
00117  
00118            0 = (0,0)      3 = (alpha,0)                        1 = (1,0)
00119 
00120            The parameters alpha, beta, gamma are determined from the given vertices in 
00121            global coordinates.
00122         */
00123         template<>
00124         class GmshReaderQuadraticBoundarySegment< 3, 3 >
00125             : public Dune::BoundarySegment< 3 >
00126         {
00127         public:
00128             GmshReaderQuadraticBoundarySegment (Dune::FieldVector<double,3> p0_, Dune::FieldVector<double,3> p1_, 
00129                 Dune::FieldVector<double,3> p2_, Dune::FieldVector<double,3> p3_,
00130                 Dune::FieldVector<double,3> p4_, Dune::FieldVector<double,3> p5_)
00131                 : p0(p0_), p1(p1_), p2(p2_), p3(p3_), p4(p4_), p5(p5_)
00132             {
00133                 sqrt2 = sqrt(2.0);
00134                 Dune::FieldVector<double,3> d1,d2;
00135 
00136                 d1 = p3; d1 -= p0;
00137                 d2 = p1; d2 -= p3;
00138                 alpha=d1.two_norm()/(d1.two_norm()+d2.two_norm());
00139                 if (alpha<1E-6 || alpha>1-1E-6)
00140                     DUNE_THROW(Dune::IOError, "alpha in quadratic boundary segment bad");
00141 
00142                 d1 = p5; d1 -= p0;
00143                 d2 = p2; d2 -= p5;
00144                 beta=d1.two_norm()/(d1.two_norm()+d2.two_norm());
00145                 if (beta<1E-6 || beta>1-1E-6)
00146                     DUNE_THROW(Dune::IOError, "beta in quadratic boundary segment bad");
00147 
00148                 d1 = p4; d1 -= p1;
00149                 d2 = p2; d2 -= p4;
00150                 gamma=sqrt2*(d1.two_norm()/(d1.two_norm()+d2.two_norm()));
00151                 if (gamma<1E-6 || gamma>1-1E-6)
00152                     DUNE_THROW(Dune::IOError, "gamma in quadratic boundary segment bad");
00153             }
00154 
00155             virtual Dune::FieldVector<double,3> operator() (const Dune::FieldVector<double,2>& local) const
00156             {
00157                 Dune::FieldVector<double,3> y;
00158                 y = 0.0;
00159                 y.axpy(phi0(local),p0);
00160                 y.axpy(phi1(local),p1);
00161                 y.axpy(phi2(local),p2);
00162                 y.axpy(phi3(local),p3);
00163                 y.axpy(phi4(local),p4);
00164                 y.axpy(phi5(local),p5);
00165                 return y;
00166             }
00167 
00168         private:
00169             // The six Lagrange basis function on the reference element
00170             // for the points given above
00171 
00172             double phi0 (const Dune::FieldVector<double,2>& local) const
00173             {
00174                 return (alpha*beta-beta*local[0]-alpha*local[1])*(1-local[0]-local[1])/(alpha*beta);
00175             }
00176             double phi3 (const Dune::FieldVector<double,2>& local) const
00177             {
00178                 return local[0]*(1-local[0]-local[1])/(alpha*(1-alpha));
00179             }
00180             double phi1 (const Dune::FieldVector<double,2>& local) const
00181             {
00182                 return local[0]*(gamma*local[0]-(sqrt2-gamma-sqrt2*alpha)*local[1]-alpha*gamma)/(gamma*(1-alpha));
00183             }
00184             double phi5 (const Dune::FieldVector<double,2>& local) const
00185             {
00186                 return local[1]*(1-local[0]-local[1])/(beta*(1-beta));
00187             }
00188             double phi4 (const Dune::FieldVector<double,2>& local) const
00189             {
00190                 return local[0]*local[1]/((1-gamma/sqrt2)*gamma/sqrt2);
00191             }
00192             double phi2 (const Dune::FieldVector<double,2>& local) const
00193             {
00194                 return local[1]*(beta*(1-gamma/sqrt2)-local[0]*(beta-gamma/sqrt2)-local[1]*(1-gamma/sqrt2))/((1-gamma/sqrt2)*(beta-1));
00195             }
00196 
00197             Dune::FieldVector<double,3> p0,p1,p2,p3,p4,p5;
00198             double alpha,beta,gamma,sqrt2;
00199         };
00200 
00201     } // end empty namespace
00202 
00204     template<typename GridType>
00205     class GmshReaderParser
00206     {
00207     protected:
00208         // private data
00209         Dune::GridFactory<GridType>& factory;
00210         bool verbose;
00211         bool insert_boundary_segments;
00212         unsigned int number_of_real_vertices;
00213         int boundary_element_count;
00214         int element_count;
00215         // read buffer
00216         char buf[512];
00217         std::string fileName;
00218         // exported data
00219         std::vector<int> boundary_id_to_physical_entity;
00220         std::vector<int> element_index_to_physical_entity;
00221 
00222         // static data
00223         static const int dim = GridType::dimension;
00224         static const int dimWorld = GridType::dimensionworld;
00225         dune_static_assert( (dimWorld <= 3), "GmshReader requires dimWorld <= 3." );
00226 
00227         // typedefs
00228         typedef FieldVector< double, dimWorld > GlobalVector;
00229 
00230         // don't use something like
00231         //   readfile(file, 1, "%s\n", buf);
00232         // to skip the rest of of the line -- that will only skip the next
00233         // whitespace-seperated word!  Use skipline() instead.
00234         void readfile(FILE * file, int cnt, const char * format,
00235                       void* t1, void* t2 = 0, void* t3 = 0, void* t4 = 0,
00236                       void* t5 = 0, void* t6 = 0, void* t7 = 0, void* t8 = 0,
00237                       void* t9 = 0, void* t10 = 0)
00238         {
00239             off_t pos = ftello(file);
00240             int c = fscanf(file, format, t1, t2, t3, t4, t5, t6, t7, t8, t9, t10);
00241             if (c != cnt)
00242                 DUNE_THROW(Dune::IOError, "Error parsing " << fileName << " "
00243                            "file pos " << pos
00244                     << ": Expected '" << format << "', only read " << c << " entries instead of " << cnt << ".");
00245         }
00246 
00247         // skip over the rest of the line, including the terminating newline
00248         void skipline(FILE * file)
00249         {
00250             int c;
00251             do {
00252                 c = std::fgetc(file);
00253             } while(c != '\n' && c != EOF);
00254         }
00255 
00256     public:
00257 
00258         GmshReaderParser(Dune::GridFactory<GridType>& _factory, bool v, bool i) :
00259             factory(_factory), verbose(v), insert_boundary_segments(i) {}
00260 
00261         std::vector<int> & boundaryIdMap()
00262         {
00263             return boundary_id_to_physical_entity;
00264         }
00265 
00266         std::vector<int> & elementIndexMap()
00267         {
00268             return element_index_to_physical_entity;
00269         }
00270       
00271         void read (const std::string& f)
00272         {
00273             if (verbose) std::cout << "Reading " << dim << "d Gmsh grid..." << std::endl;
00274 
00275             // open file name, we use C I/O
00276             fileName = f;
00277             FILE* file = fopen(fileName.c_str(),"r");
00278             if (file==0)
00279                 DUNE_THROW(Dune::IOError, "Could not open " << fileName);
00280 
00281             //=========================================
00282             // Header: Read vertices into vector
00283             //         Check vertices that are needed
00284             //=========================================
00285 
00286             number_of_real_vertices = 0;
00287             boundary_element_count = 0;
00288             element_count = 0;
00289       
00290             // process header
00291             double version_number;
00292             int file_type, data_size;
00293       
00294             readfile(file,1,"%s\n",buf);
00295             if (strcmp(buf,"$MeshFormat")!=0)
00296                 DUNE_THROW(Dune::IOError, "expected $MeshFormat in first line");
00297             readfile(file,3,"%lg %d %d\n",&version_number,&file_type,&data_size);
00298             if( (version_number < 2.0) || (version_number > 2.2) )
00299                 DUNE_THROW(Dune::IOError, "can only read Gmsh version 2 files");
00300             if (verbose) std::cout << "version " << version_number << " Gmsh file detected" << std::endl; 
00301             readfile(file,1,"%s\n",buf);
00302             if (strcmp(buf,"$EndMeshFormat")!=0)
00303                 DUNE_THROW(Dune::IOError, "expected $EndMeshFormat");
00304       
00305             // node section
00306             int number_of_nodes;
00307       
00308             readfile(file,1,"%s\n",buf);
00309             if (strcmp(buf,"$Nodes")!=0)
00310                 DUNE_THROW(Dune::IOError, "expected $Nodes");
00311             readfile(file,1,"%d\n",&number_of_nodes);
00312             if (verbose) std::cout << "file contains " << number_of_nodes << " nodes" << std::endl;
00313 
00314             // read nodes
00315             std::vector< GlobalVector > nodes( number_of_nodes+1 ); // store positions
00316             {
00317                 int id;
00318                 double x[ 3 ];
00319                 for( int i = 1; i <= number_of_nodes; ++i )
00320                 {
00321                     readfile(file,4, "%d %lg %lg %lg\n", &id, &x[ 0 ], &x[ 1 ], &x[ 2 ] );
00322                     if( id != i )
00323                         DUNE_THROW( Dune::IOError, "Expected id " << i << "(got id " << id << "." );
00324           
00325                     // just store node position
00326                     for( int j = 0; j < dimWorld; ++j )
00327                         nodes[ i ][ j ] = x[ j ];
00328                 }
00329                 readfile(file,1,"%s\n",buf);
00330                 if (strcmp(buf,"$EndNodes")!=0)
00331                     DUNE_THROW(Dune::IOError, "expected $EndNodes");
00332             }
00333   
00334             // element section
00335             readfile(file,1,"%s\n",buf);
00336             if (strcmp(buf,"$Elements")!=0)
00337                 DUNE_THROW(Dune::IOError, "expected $Elements");
00338             int number_of_elements;
00339             readfile(file,1,"%d\n",&number_of_elements);
00340             if (verbose) std::cout << "file contains " << number_of_elements << " elements" << std::endl;
00341 
00342             //=========================================
00343             // Pass 1: Renumber needed vertices
00344             //=========================================
00345 
00346             long section_element_offset = ftell(file);
00347             std::map<int,unsigned int> renumber;
00348             for (int i=1; i<=number_of_elements; i++)
00349             {
00350                 int id, elm_type, number_of_tags;
00351                 readfile(file,3,"%d %d %d ",&id,&elm_type,&number_of_tags);
00352                 for (int k=1; k<=number_of_tags; k++)
00353                 {
00354                     int blub;
00355                     readfile(file,1,"%d ",&blub);
00356                     // k == 1: physical entity (not used here)
00357                     // k == 2: elementary entity (not used here either)
00358                     // if version_number < 2.2:
00359                     //   k == 3: mesh partition 0
00360                     // else
00361                     //   k == 3: number of mesh partitions
00362                     //   k => 4: mesh partition k-4
00363                 }
00364                 pass1HandleElement(file, elm_type, renumber, nodes);
00365             }
00366             if (verbose) std::cout << "number of real vertices = " << number_of_real_vertices << std::endl;
00367             if (verbose) std::cout << "number of boundary elements = " << boundary_element_count << std::endl;
00368             if (verbose) std::cout << "number of elements = " << element_count << std::endl;
00369             readfile(file,1,"%s\n",buf);
00370             if (strcmp(buf,"$EndElements")!=0)
00371                 DUNE_THROW(Dune::IOError, "expected $EndElements");
00372             boundary_id_to_physical_entity.resize(boundary_element_count);
00373             element_index_to_physical_entity.resize(element_count);
00374 
00375             //==============================================
00376             // Pass 2: Insert boundary segments and elements
00377             //==============================================
00378 
00379             fseek(file, section_element_offset, SEEK_SET);
00380             boundary_element_count = 0;
00381             element_count = 0;
00382             for (int i=1; i<=number_of_elements; i++)
00383             {
00384                 int id, elm_type, number_of_tags;
00385                 readfile(file,3,"%d %d %d ",&id,&elm_type,&number_of_tags);
00386                 int physical_entity = -1;
00387                 std::vector<int> mesh_partitions;
00388                 if ( version_number < 2.2 )
00389                 {
00390                     mesh_partitions.resize(1);
00391                 }
00392                 for (int k=1; k<=number_of_tags; k++)
00393                 {
00394                     int blub;
00395                     readfile(file,1,"%d ",&blub);
00396                     if (k==1) physical_entity = blub;
00397                     // k == 2: elementary entity (not used here)
00398                     if ( version_number < 2.2 )
00399                     {
00400                         if (k==3) mesh_partitions[0] = blub;
00401                     }
00402                     else
00403                     {
00404                         if (k > 3)
00405                             mesh_partitions[k-4] = blub;
00406                         else
00407                             mesh_partitions.resize(blub);
00408                     }
00409                 }
00410                 pass2HandleElement(file, elm_type, renumber, nodes, physical_entity);
00411             }
00412             readfile(file,1,"%s\n",buf);
00413             if (strcmp(buf,"$EndElements")!=0)
00414                 DUNE_THROW(Dune::IOError, "expected $EndElements");
00415 
00416             fclose(file);
00417         }
00418         
00419                 // dimension dependent routines
00420         void pass1HandleElement(FILE* file, const int elm_type,
00421             std::map<int,unsigned int> & renumber,
00422             const std::vector< GlobalVector > & nodes)
00423         {
00424             // some data about gmsh elements
00425             const int nDofs[12]      = {-1, 2, 3, 4, 4, 8, 6, 5, 3, 6, -1, 10};
00426             const int nVertices[12]  = {-1, 2, 3, 4, 4, 8, 6, 5, 2, 3, -1, 4};
00427             const int elementDim[12] = {-1, 1, 2, 2, 3, 3, 3, 3, 1, 2, -1, 3};
00428             
00429             // test whether we support the element type
00430             if ( not (elm_type >= 0 && elm_type < 12   // index in suitable range?
00431                       && (elementDim[elm_type] == dim || elementDim[elm_type] == (dim-1) ) ) )   // real element or boundary element?
00432             {
00433                 skipline(file); // skip rest of line if element is unknown
00434                 return;
00435             }
00436             
00437             // The format string for parsing is n times '%d' in a row
00438             std::string formatString = "%d";
00439             for (int i=1; i<nDofs[elm_type]; i++)
00440                 formatString += " %d";
00441             formatString += "\n";
00442             
00443             // '10' is the largest number of dofs we may encounter in a .msh file
00444             std::vector<int> elementDofs(10);
00445 
00446             readfile(file,nDofs[elm_type], formatString.c_str(),
00447                      &(elementDofs[0]),&(elementDofs[1]),&(elementDofs[2]),
00448                      &(elementDofs[3]),&(elementDofs[4]),&(elementDofs[5]),
00449                      &(elementDofs[6]),&(elementDofs[7]),&(elementDofs[8]),
00450                      &(elementDofs[9]));
00451 
00452             // insert each vertex if it hasn't been inserted already
00453             for (int i=0; i<nVertices[elm_type]; i++)
00454                 if (renumber.find(elementDofs[i])==renumber.end()) 
00455                 {
00456                     renumber[elementDofs[i]] = number_of_real_vertices++;
00457                     factory.insertVertex(nodes[elementDofs[i]]);                  
00458                 }
00459                 
00460             // count elements and boundary elements
00461             if (elementDim[elm_type] == dim)
00462                 element_count++;
00463             else
00464                 boundary_element_count++;
00465             
00466         }
00467 
00468 
00469 
00470         // generic-case: This is not supposed to be used at runtime.
00471         template <class E, class V, class V2>
00472         void boundarysegment_insert(
00473                                     const V& nodes,
00474                                     const E& elementDofs,
00475                                     const V2& vertices
00476                                     ) 
00477         {
00478             DUNE_THROW(Dune::IOError, "tried to create a 3D boundary segment in a non-3D Grid");
00479         }
00480         
00481         // 3d-case:
00482         template <class E, class V>
00483         void boundarysegment_insert(
00484                                     const std::vector<FieldVector<double, 3> >& nodes,
00485                                     const E& elementDofs,
00486                                     const V& vertices
00487                                     ) 
00488         {
00489             array<FieldVector<double,dimWorld>, 6> v;
00490             for (int i=0; i<6; i++) 
00491                 for (int j=0; j<dimWorld; j++)
00492                     v[i][j] = nodes[elementDofs[i]][j];
00493             
00494             BoundarySegment<dim,dimWorld>* newBoundarySegment 
00495                 = (BoundarySegment<dim,dimWorld>*) new GmshReaderQuadraticBoundarySegment< 3, 3 >( v[0], v[1], v[2],
00496                                                                                                   v[3], v[4], v[5] );
00497             
00498             factory.insertBoundarySegment( vertices,
00499                                            shared_ptr<BoundarySegment<dim,dimWorld> >(newBoundarySegment) );
00500         }
00501         
00502         
00503 
00504 
00505 
00506         virtual void pass2HandleElement(FILE* file, const int elm_type,
00507             std::map<int,unsigned int> & renumber,
00508             const std::vector< GlobalVector > & nodes,
00509             const int physical_entity)
00510         {
00511             // some data about gmsh elements
00512             const int nDofs[12]      = {-1, 2, 3, 4, 4, 8, 6, 5, 3, 6, -1, 10};
00513             const int nVertices[12]  = {-1, 2, 3, 4, 4, 8, 6, 5, 2, 3, -1, 4};
00514             const int elementDim[12] = {-1, 1, 2, 2, 3, 3, 3, 3, 1, 2, -1, 3};
00515             
00516             // test whether we support the element type
00517             if ( not (elm_type >= 0 && elm_type < 12   // index in suitable range?
00518                       && (elementDim[elm_type] == dim || elementDim[elm_type] == (dim-1) ) ) )   // real element or boundary element?
00519             {
00520                 skipline(file); // skip rest of line if element is unknown
00521                 return;
00522             }
00523             
00524             // The format string for parsing is n times '%d' in a row
00525             std::string formatString = "%d";
00526             for (int i=1; i<nDofs[elm_type]; i++)
00527                 formatString += " %d";
00528             formatString += "\n";
00529             
00530             // '10' is the largest number of dofs we may encounter in a .msh file
00531             std::vector<int> elementDofs(10);
00532 
00533             readfile(file,nDofs[elm_type], formatString.c_str(),
00534                      &(elementDofs[0]),&(elementDofs[1]),&(elementDofs[2]),
00535                      &(elementDofs[3]),&(elementDofs[4]),&(elementDofs[5]),
00536                      &(elementDofs[6]),&(elementDofs[7]),&(elementDofs[8]),
00537                      &(elementDofs[9]));
00538 
00539             // correct differences between gmsh and Dune in the local vertex numbering
00540             switch (elm_type)
00541             {
00542                 case 3: // 4-node quadrilateral
00543                     std::swap(elementDofs[2],elementDofs[3]);
00544                     break;
00545                 case 5: // 8-node hexahedron
00546                     std::swap(elementDofs[2],elementDofs[3]);
00547                     std::swap(elementDofs[6],elementDofs[7]);
00548                     break;
00549                 case 7: // 5-node pyramid
00550                     std::swap(elementDofs[2],elementDofs[3]);
00551                     break;
00552             }
00553 
00554             // renumber corners to account for the explicitly given vertex
00555             // numbering in the file
00556             std::vector<unsigned int> vertices(nVertices[elm_type]);
00557 
00558             for (int i=0; i<nVertices[elm_type]; i++)
00559                 vertices[i] = renumber[elementDofs[i]]; 
00560 
00561             // If it is an element, insert it as such
00562             if (elementDim[elm_type] == dim) {
00563                 
00564                 switch (elm_type)
00565                 {
00566                     case 1: // 2-node line
00567                         factory.insertElement(Dune::GeometryType(Dune::GeometryType::simplex,dim),vertices);
00568                         break;
00569                     case 2: // 3-node triangle
00570                         factory.insertElement(Dune::GeometryType(Dune::GeometryType::simplex,dim),vertices);
00571                         break;
00572                     case 3: // 4-node quadrilateral
00573                         factory.insertElement(Dune::GeometryType(Dune::GeometryType::cube,dim),vertices);
00574                         break;
00575                     case 4: // 4-node tetrahedron
00576                         factory.insertElement(Dune::GeometryType(Dune::GeometryType::simplex,dim),vertices);
00577                         break;
00578                     case 5: // 8-node hexahedron
00579                         factory.insertElement(Dune::GeometryType(Dune::GeometryType::cube,dim),vertices);
00580                         break;
00581                     case 6: // 6-node prism
00582                         factory.insertElement(Dune::GeometryType(Dune::GeometryType::prism,dim),vertices);
00583                         break;
00584                     case 7: // 5-node pyramid
00585                         factory.insertElement(Dune::GeometryType(Dune::GeometryType::pyramid,dim),vertices);
00586                         break;
00587                     case 9: // 6-node triangle
00588                         factory.insertElement(Dune::GeometryType(Dune::GeometryType::simplex,dim),vertices);
00589                         break;
00590                     case 11: // 10-node tetrahedron
00591                         factory.insertElement(Dune::GeometryType(Dune::GeometryType::simplex,dim),vertices);
00592                         break;
00593                 }
00594 
00595             } else {
00596                 // it must be a boundary segment then
00597                 if (insert_boundary_segments) {
00598                     
00599                     switch (elm_type)
00600                     {
00601                         case 1: // 2-node line
00602                             factory.insertBoundarySegment(vertices);
00603                             break;
00604 
00605                         case 2: // 3-node triangle
00606                             factory.insertBoundarySegment(vertices);
00607                             break;
00608 
00609                         case 8: { // 3-node line
00610                             array<FieldVector<double,dimWorld>, 3> v;
00611                             for (int i=0; i<dimWorld; i++) {
00612                                 v[0][i] = nodes[elementDofs[0]][i];
00613                                 v[1][i] = nodes[elementDofs[2]][i];  // yes, the renumbering is intended!
00614                                 v[2][i] = nodes[elementDofs[1]][i];
00615                             }
00616                             BoundarySegment<dim,dimWorld>* newBoundarySegment 
00617                                 = (BoundarySegment<dim,dimWorld>*)new GmshReaderQuadraticBoundarySegment< 2, dimWorld >(v[0], v[1], v[2]);
00618                             factory.insertBoundarySegment(vertices, 
00619                                                           shared_ptr<BoundarySegment<dim,dimWorld> >(newBoundarySegment));
00620                             break;
00621                         }
00622                         case 9: { // 6-node triangle
00623                             boundarysegment_insert(nodes, elementDofs, vertices);
00624                             break;
00625                         }
00626                         
00627                     }
00628                     
00629                 }
00630             }
00631             
00632             // count elements and boundary elements
00633             if (elementDim[elm_type] == dim) {
00634                 element_index_to_physical_entity[element_count] = physical_entity;
00635                 element_count++;
00636             } else {
00637                 boundary_id_to_physical_entity[boundary_element_count] = physical_entity;
00638                 boundary_element_count++;
00639             }
00640 
00641         }
00642 
00643     };
00644 
00661     template<typename GridType>
00662     class GmshReader
00663     {
00664     public:
00665         typedef GridType Grid;
00666         
00668         static Grid* read (const std::string& fileName, bool verbose = true, bool insert_boundary_segments=true)
00669         { 
00670             // make a grid factory
00671             Dune::GridFactory<Grid> factory;
00672       
00673             // create parse object
00674             GmshReaderParser<Grid> parser(factory,verbose,insert_boundary_segments);
00675             parser.read(fileName);
00676 
00677             return factory.createGrid();
00678         }
00679 
00681         static Grid* read (const std::string& fileName,
00682             std::vector<int>& boundary_id_to_physical_entity,
00683             std::vector<int>& element_index_to_physical_entity,
00684             bool verbose = true, bool insert_boundary_segments=true)
00685         { 
00686             // make a grid factory
00687             Dune::GridFactory<Grid> factory;
00688 
00689             // create parse object
00690             GmshReaderParser<Grid> parser(factory,verbose,insert_boundary_segments);
00691             parser.read(fileName);
00692 
00693             boundary_id_to_physical_entity.swap(parser.boundaryIdMap());
00694             element_index_to_physical_entity.swap(parser.elementIndexMap());
00695 
00696             return factory.createGrid(); 
00697         }
00698 
00700         static void read (Dune::GridFactory<Grid>& factory, const std::string& fileName,
00701             bool verbose = true, bool insert_boundary_segments=true)
00702         { 
00703             // create parse object
00704             GmshReaderParser<Grid> parser(factory,verbose,insert_boundary_segments);
00705             parser.read(fileName);
00706         }
00707 
00709         static void read (Dune::GridFactory<Grid>& factory, 
00710             const std::string& fileName, 
00711             std::vector<int>& boundary_id_to_physical_entity,
00712             std::vector<int>& element_index_to_physical_entity,
00713             bool verbose = true, bool insert_boundary_segments=true)
00714         { 
00715             // create parse object
00716             GmshReaderParser<Grid> parser(factory,verbose,insert_boundary_segments);
00717             parser.read(fileName);
00718 
00719             boundary_id_to_physical_entity.swap(parser.boundaryIdMap());
00720             element_index_to_physical_entity.swap(parser.elementIndexMap());
00721         }
00722     };
00723 
00726 } // namespace Dune
00727 
00728 #endif