dune-grid
2.2.0
|
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