dune-grid
2.2.0
|
00001 #ifndef DUNE_STARCD_READER_HH 00002 #define DUNE_STARCD_READER_HH 00003 00004 #include <dune/common/exceptions.hh> 00005 #include <dune/geometry/type.hh> 00006 #include <dune/grid/common/gridfactory.hh> 00007 #include <iostream> 00008 #include <fstream> 00009 00010 namespace Dune { 00011 00045 template <class GridType> 00046 class StarCDReader { 00047 00048 public: 00049 00055 static GridType* read(const std::string& fileName, bool verbose = true) 00056 { 00057 // extract the grid dimension 00058 const int dim = GridType::dimension; 00059 00060 // currently only dim = 3 is implemented 00061 if (dim != 3) 00062 DUNE_THROW(Dune::NotImplemented, 00063 "Reading Star-CD format is not implemented for dimension " << dim); 00064 00065 // set up the grid factory 00066 GridFactory<GridType> factory; 00067 00068 // set the name of the vertex file 00069 std::string vertexFileName = fileName + ".vrt"; 00070 00071 // set the vertex input stream 00072 std::ifstream vertexFile(vertexFileName.c_str()); 00073 if (!vertexFile) 00074 DUNE_THROW(Dune::IOError, "Could not open " << vertexFileName); 00075 00076 // read the vertices 00077 int dummyIdx; 00078 int numberOfVertices = 0; 00079 while (vertexFile >> dummyIdx) { 00080 numberOfVertices++; 00081 00082 Dune::FieldVector<double,dim> position; 00083 00084 for (int k = 0; k < dim; k++) 00085 vertexFile >> position[k]; 00086 00087 factory.insertVertex(position); 00088 } 00089 if (verbose) 00090 std::cout << numberOfVertices << " vertices read." << std::endl; 00091 00092 // set the name of the element file 00093 std::string elementFileName = fileName + ".cel"; 00094 00095 // set the element input stream 00096 std::ifstream elementFile(elementFileName.c_str()); 00097 if (!elementFile) 00098 DUNE_THROW(Dune::IOError, "Could not open " << elementFileName); 00099 00100 // read the elements 00101 int numberOfElements = 0; 00102 int numberOfSimplices = 0; 00103 int numberOfPyramids = 0; 00104 int numberOfPrisms = 0; 00105 int numberOfCubes = 0;; 00106 int maxNumberOfVertices = (int)pow(2, dim); 00107 int isVolume = 1; 00108 while (elementFile >> dummyIdx) { 00109 std::vector<unsigned int> vertices(maxNumberOfVertices); 00110 for (int k = 0; k < maxNumberOfVertices; k++) 00111 elementFile >> vertices[k]; 00112 00113 int boundaryId; 00114 elementFile >> boundaryId; 00115 00116 int volumeOrSurface[2]; 00117 elementFile >> volumeOrSurface[0] >> volumeOrSurface[1]; 00118 00119 if (volumeOrSurface[0] == isVolume) { 00120 numberOfElements++; 00121 00122 if (vertices[2] == vertices[3]) { // simplex or prism 00123 if (vertices[4] == vertices[5]) { // simplex 00124 numberOfSimplices++; 00125 std::vector<unsigned int> simplexVertices(4); 00126 for (int k = 0; k < 3; k++) 00127 simplexVertices[k] = vertices[k] - 1; 00128 simplexVertices[3] = vertices[4] - 1; 00129 factory.insertElement(Dune::GeometryType(Dune::GeometryType::simplex,dim), simplexVertices); 00130 } 00131 else { // prism 00132 numberOfPrisms++; 00133 std::vector<unsigned int> prismVertices(6); 00134 for (int k = 0; k < 3; k++) 00135 prismVertices[k] = vertices[k] - 1; 00136 for (int k = 3; k < 6; k++) 00137 prismVertices[k] = vertices[k+1] - 1; 00138 factory.insertElement(Dune::GeometryType(Dune::GeometryType::prism,dim), prismVertices); 00139 } 00140 } 00141 else { // cube or pyramid 00142 if (vertices[4] == vertices[5]) { // pyramid 00143 numberOfPyramids++; 00144 std::vector<unsigned int> pyramidVertices(5); 00145 for (int k = 0; k < 5; k++) 00146 pyramidVertices[k] = vertices[k] - 1; 00147 factory.insertElement(Dune::GeometryType(Dune::GeometryType::pyramid,dim), pyramidVertices); 00148 } 00149 else { // cube 00150 numberOfCubes++; 00151 std::vector<unsigned int> cubeVertices(8); 00152 for (int k = 0; k < 8; k++) 00153 cubeVertices[k] = vertices[k] - 1; 00154 std::swap(cubeVertices[2], cubeVertices[3]); 00155 std::swap(cubeVertices[6], cubeVertices[7]); 00156 factory.insertElement(Dune::GeometryType(Dune::GeometryType::cube,dim), cubeVertices); 00157 } 00158 } 00159 } 00160 } 00161 if (verbose) 00162 std::cout << numberOfElements << " elements read: " 00163 << numberOfSimplices << " simplices, " << numberOfPyramids << " pyramids, " 00164 << numberOfPrisms << " prisms, " << numberOfCubes << " cubes." << std::endl; 00165 00166 // finish off the construction of the grid object 00167 if (verbose) 00168 std::cout << "Starting createGrid() ... " << std::flush; 00169 00170 return factory.createGrid(); 00171 00172 } 00173 00174 }; 00175 00176 } 00177 00178 #endif