dune-grid  2.2.1
starcdreader.hh
Go to the documentation of this file.
1 #ifndef DUNE_STARCD_READER_HH
2 #define DUNE_STARCD_READER_HH
3 
4 #include <dune/common/exceptions.hh>
5 #include <dune/geometry/type.hh>
7 #include <iostream>
8 #include <fstream>
9 
10 namespace Dune {
11 
45  template <class GridType>
46  class StarCDReader {
47 
48  public:
49 
55  static GridType* read(const std::string& fileName, bool verbose = true)
56  {
57  // extract the grid dimension
58  const int dim = GridType::dimension;
59 
60  // currently only dim = 3 is implemented
61  if (dim != 3)
62  DUNE_THROW(Dune::NotImplemented,
63  "Reading Star-CD format is not implemented for dimension " << dim);
64 
65  // set up the grid factory
66  GridFactory<GridType> factory;
67 
68  // set the name of the vertex file
69  std::string vertexFileName = fileName + ".vrt";
70 
71  // set the vertex input stream
72  std::ifstream vertexFile(vertexFileName.c_str());
73  if (!vertexFile)
74  DUNE_THROW(Dune::IOError, "Could not open " << vertexFileName);
75 
76  // read the vertices
77  int dummyIdx;
78  int numberOfVertices = 0;
79  while (vertexFile >> dummyIdx) {
80  numberOfVertices++;
81 
82  Dune::FieldVector<double,dim> position;
83 
84  for (int k = 0; k < dim; k++)
85  vertexFile >> position[k];
86 
87  factory.insertVertex(position);
88  }
89  if (verbose)
90  std::cout << numberOfVertices << " vertices read." << std::endl;
91 
92  // set the name of the element file
93  std::string elementFileName = fileName + ".cel";
94 
95  // set the element input stream
96  std::ifstream elementFile(elementFileName.c_str());
97  if (!elementFile)
98  DUNE_THROW(Dune::IOError, "Could not open " << elementFileName);
99 
100  // read the elements
101  int numberOfElements = 0;
102  int numberOfSimplices = 0;
103  int numberOfPyramids = 0;
104  int numberOfPrisms = 0;
105  int numberOfCubes = 0;;
106  int maxNumberOfVertices = (int)pow(2, dim);
107  int isVolume = 1;
108  while (elementFile >> dummyIdx) {
109  std::vector<unsigned int> vertices(maxNumberOfVertices);
110  for (int k = 0; k < maxNumberOfVertices; k++)
111  elementFile >> vertices[k];
112 
113  int boundaryId;
114  elementFile >> boundaryId;
115 
116  int volumeOrSurface[2];
117  elementFile >> volumeOrSurface[0] >> volumeOrSurface[1];
118 
119  if (volumeOrSurface[0] == isVolume) {
120  numberOfElements++;
121 
122  if (vertices[2] == vertices[3]) { // simplex or prism
123  if (vertices[4] == vertices[5]) { // simplex
124  numberOfSimplices++;
125  std::vector<unsigned int> simplexVertices(4);
126  for (int k = 0; k < 3; k++)
127  simplexVertices[k] = vertices[k] - 1;
128  simplexVertices[3] = vertices[4] - 1;
129  factory.insertElement(Dune::GeometryType(Dune::GeometryType::simplex,dim), simplexVertices);
130  }
131  else { // prism
132  numberOfPrisms++;
133  std::vector<unsigned int> prismVertices(6);
134  for (int k = 0; k < 3; k++)
135  prismVertices[k] = vertices[k] - 1;
136  for (int k = 3; k < 6; k++)
137  prismVertices[k] = vertices[k+1] - 1;
138  factory.insertElement(Dune::GeometryType(Dune::GeometryType::prism,dim), prismVertices);
139  }
140  }
141  else { // cube or pyramid
142  if (vertices[4] == vertices[5]) { // pyramid
143  numberOfPyramids++;
144  std::vector<unsigned int> pyramidVertices(5);
145  for (int k = 0; k < 5; k++)
146  pyramidVertices[k] = vertices[k] - 1;
147  factory.insertElement(Dune::GeometryType(Dune::GeometryType::pyramid,dim), pyramidVertices);
148  }
149  else { // cube
150  numberOfCubes++;
151  std::vector<unsigned int> cubeVertices(8);
152  for (int k = 0; k < 8; k++)
153  cubeVertices[k] = vertices[k] - 1;
154  std::swap(cubeVertices[2], cubeVertices[3]);
155  std::swap(cubeVertices[6], cubeVertices[7]);
156  factory.insertElement(Dune::GeometryType(Dune::GeometryType::cube,dim), cubeVertices);
157  }
158  }
159  }
160  }
161  if (verbose)
162  std::cout << numberOfElements << " elements read: "
163  << numberOfSimplices << " simplices, " << numberOfPyramids << " pyramids, "
164  << numberOfPrisms << " prisms, " << numberOfCubes << " cubes." << std::endl;
165 
166  // finish off the construction of the grid object
167  if (verbose)
168  std::cout << "Starting createGrid() ... " << std::flush;
169 
170  return factory.createGrid();
171 
172  }
173 
174  };
175 
176 }
177 
178 #endif