3.1-git
DUNE for Multi-{Phase, Component, Scale, Physics, ...} flow and transport in porous media
gmshreader.hh
Go to the documentation of this file.
1/*****************************************************************************
2 * See the file COPYING for full copying permissions. *
3 * *
4 * This program is free software: you can redistribute it and/or modify *
5 * it under the terms of the GNU General Public License as published by *
6 * the Free Software Foundation, either version 3 of the License, or *
7 * (at your option) any later version. *
8 * *
9 * This program is distributed in the hope that it will be useful, *
10 * but WITHOUT ANY WARRANTY; without even the implied warranty of *
11 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
12 * GNU General Public License for more details. *
13 * *
14 * You should have received a copy of the GNU General Public License *
15 * along with this program. If not, see <http://www.gnu.org/licenses/>. *
16 *****************************************************************************/
22#ifndef DUMUX_FACETCOUPLING_GMSH_READER_HH
23#define DUMUX_FACETCOUPLING_GMSH_READER_HH
24
25#include <algorithm>
26#include <cassert>
27#include <fstream>
28#include <iostream>
29#include <sstream>
30#include <typeinfo>
31#include <unordered_map>
32
33#include <dune/common/timer.hh>
34#include <dune/common/version.hh>
35#include <dune/common/fvector.hh>
36#include <dune/geometry/type.hh>
37
39
40namespace Dumux {
41
59template <class BulkGrid, int numGrids>
61{
62 // extract some necessary info from bulk grid
63 static constexpr int bulkDim = BulkGrid::dimension;
64 static constexpr int bulkDimWorld = BulkGrid::dimensionworld;
65 using ctype = typename BulkGrid::ctype;
67 using GlobalPosition = Dune::FieldVector<ctype, bulkDimWorld>;
68
69 // determine minimum dimension for which a grid is created
70 static constexpr int minGridDim = BulkGrid::dimension - numGrids + 1;
71 static_assert(minGridDim >= 1, "Grids with dim < 1 cannot be read!");
72
73 // structure to store data on an element
74 using VertexIndexSet = std::vector<GridIndexType>;
75 struct ElementData
76 {
77 Dune::GeometryType gt;
78 VertexIndexSet cornerIndices;
79 };
80
81public:
84 void read(const std::string& fileName, bool verbose = false)
85 {
86 read(fileName, 0, verbose);
87 }
88
90 void read(const std::string& fileName, std::size_t boundarySegThresh, bool verbose = false)
91 {
92 Dune::Timer watch;
93 if (verbose) std::cout << "Opening " << fileName << std::endl;
94 std::ifstream gridFile(fileName);
95 if (gridFile.fail())
96 DUNE_THROW(Dune::InvalidStateException, "Could not open the given .msh file. Make sure it exists");
97
98 // currently we only support version 2 file format
99 std::string line;
100 std::getline(gridFile, line);
101 if (line.find("$MeshFormat") == std::string::npos)
102 DUNE_THROW(Dune::InvalidStateException, "Expected $MeshFormat in the first line of the grid file!");
103 std::getline(gridFile, line);
104 if (line.find_first_of("2") != 0)
105 DUNE_THROW(Dune::InvalidStateException, "Currently only Gmsh mesh file format version 2 is supported!");
106
107 // read file until we get to the list of nodes
108 while (line.find("$Nodes") == std::string::npos)
109 std::getline(gridFile, line);
110
111 // read all vertices
112 std::getline(gridFile, line);
113 const auto numVertices = convertString<std::size_t>(line);
114 gridVertices_.resize(numVertices);
115
116 std::getline(gridFile, line);
117 std::size_t vertexCount = 0;
118 while (line.find("$EndNodes") == std::string::npos)
119 {
120 // drop first entry in line (vertex index) and read coordinates
121 std::istringstream stream(line);
122 std::string buf; stream >> buf;
123 GlobalPosition v;
124 for (auto& coord : v)
125 {
126 stream >> coord;
127 if (stream.fail()) DUNE_THROW(Dune::IOError, "Could not read vertex coordinate");
128 }
129
130 // insert vertex into container and move to next line
131 gridVertices_[vertexCount++] = v;
132 std::getline(gridFile, line);
133 }
134
135 // we should always find as many vertices as the mesh file states
136 if (vertexCount != numVertices)
137 DUNE_THROW(Dune::InvalidStateException, "Couldn't find as many vertices as stated in the .msh file");
138
139 // read file until we get to the list of elements
140 while(line.find("$Elements") == std::string::npos)
141 std::getline(gridFile, line);
142
143 // read elements
144 std::getline(gridFile, line);
145 const auto numElements = convertString<std::size_t>(line);
146
147 // keep track of number of elements
148 std::array<std::size_t, numGrids> elementCount;
149 std::fill(elementCount.begin(), elementCount.end(), 0);
150
151 // Construct maps that map bulk grid vertex
152 // indices to lowDim vertex indices. -1 indicates non-initialized status
153 std::size_t elemCount = 0;
154 std::array<std::size_t, numGrids> gridVertexCount;
155 std::array<std::vector<GridIndexType>, numGrids> gridVertexMap;
156 std::array<std::vector<bool>, numGrids> idxIsAssigned;
157 std::fill(gridVertexCount.begin(), gridVertexCount.end(), 0);
158 std::fill(gridVertexMap.begin(), gridVertexMap.end(), std::vector<GridIndexType>(vertexCount));
159 std::fill(idxIsAssigned.begin(), idxIsAssigned.end(), std::vector<bool>(vertexCount, false));
160 std::getline(gridFile, line);
161 while (line.find("$EndElements") == std::string::npos)
162 {
163 // pass all indices into vector
164 std::istringstream stream(line);
165 std::string buf;
166 std::vector<std::size_t> lineData;
167 while (stream >> buf) lineData.push_back(convertString<std::size_t>(buf));
168 assert(lineData.size() >= 4 && "Grid format erroneous or unsupported");
169
170 // obtain geometry type
171 const auto gt = obtainGeometryType( lineData[1] );
172 const std::size_t physicalIndex = lineData[3];
173 const auto geoDim = gt.dim();
174 const bool isBoundarySeg = geoDim != bulkDim && physicalIndex < boundarySegThresh;
175 if (geoDim >= minGridDim-1)
176 {
177 // insert boundary segment
178 if ((isBoundarySeg || geoDim == minGridDim-1))
179 {
180 const unsigned int nextLevelGridIdx = bulkDim-geoDim-1;
181
182 VertexIndexSet corners;
183 auto it = lineData.begin()+2+lineData[2]+1;
184 for (; it != lineData.end(); ++it)
185 {
186 *it -= 1; // gmsh indices start from 1
187
188 // insert map if vertex is not inserted yet
189 if (!idxIsAssigned[nextLevelGridIdx][*it])
190 {
191 gridVertexMap[nextLevelGridIdx][*it] = gridVertexCount[nextLevelGridIdx]++;
192 idxIsAssigned[nextLevelGridIdx][*it] = true;
193 vertexIndices_[nextLevelGridIdx].push_back(*it);
194 }
195
196 corners.push_back(gridVertexMap[nextLevelGridIdx][*it]);
197 }
198
199 // marker = physical entity index
200 boundaryMarkerMaps_[nextLevelGridIdx].push_back(physicalIndex);
201 boundarySegments_[nextLevelGridIdx].push_back(corners);
202 }
203
204 // insert element
205 else
206 {
207 const unsigned int gridIdx = bulkDim-geoDim;
208
209 VertexIndexSet corners;
210 auto it = lineData.begin()+2+lineData[2]+1;
211 for (; it != lineData.end(); ++it)
212 {
213 *it -= 1; // gmsh indices start from 1
214
215 // insert map if vertex is not inserted yet
216 if (!idxIsAssigned[gridIdx][*it])
217 {
218 gridVertexMap[gridIdx][*it] = gridVertexCount[gridIdx]++;
219 idxIsAssigned[gridIdx][*it] = true;
220 vertexIndices_[gridIdx].push_back(*it);
221 }
222
223 corners.push_back(gridVertexMap[gridIdx][*it]);
224 }
225
226 // add data to embedments/embeddings
227 if (geoDim > minGridDim)
228 {
229 const auto gridElemCount = elementData_[gridIdx].size();
230 const auto& embeddedVIndices = vertexIndices_[gridIdx+1];
231 const auto& embeddedIndicesAssigned = idxIsAssigned[gridIdx+1];
232
233 VertexIndexSet cornerIndicesGlobal(corners.size());
234 for (unsigned int i = 0; i < corners.size(); ++i)
235 cornerIndicesGlobal[i] = vertexIndices_[gridIdx][corners[i]];
236 addEmbeddings(cornerIndicesGlobal, gridIdx, gridElemCount, embeddedVIndices, embeddedIndicesAssigned);
237 }
238
239 // ensure dune-specific corner ordering
240 reorder(gt, corners);
241
242 // insert element data to grid's container
243 elementMarkerMaps_[gridIdx].push_back(physicalIndex);
244 elementData_[gridIdx].emplace_back(ElementData({gt, corners}));
245 }
246 }
247
248 // get next line
249 std::getline(gridFile, line);
250 elemCount++;
251 }
252
253 // make sure we read all elements
254 if (elemCount != numElements)
255 DUNE_THROW(Dune::InvalidStateException, "Didn't read as many elements as stated in the .msh file");
256
257 if (verbose)
258 {
259 std::cout << "Finished reading gmsh file" << std::endl;
260 for (std::size_t id = 0; id < numGrids; ++id)
261 {
262 std::cout << elementData_[id].size() << " "
263 << bulkDim-id << "-dimensional elements comprising of "
264 << gridVertexCount[id] << " vertices";
265 if (id < numGrids-1) std::cout << "," << std::endl;
266 }
267 std::cout << " have been read in " << watch.elapsed() << " seconds." << std::endl;
268 }
269 }
270
272 const std::vector<GlobalPosition>& gridVertices() const
273 { return gridVertices_; }
274
276 VertexIndexSet& vertexIndices(std::size_t id)
277 {
278 assert(id < numGrids && "Index exceeds number of grids provided");
279 return vertexIndices_[id];
280 }
281
283 const std::vector<ElementData>& elementData(std::size_t id) const
284 {
285 assert(id < numGrids && "Index exceeds number of grids provided");
286 return elementData_[id];
287 }
288
290 const std::vector<VertexIndexSet>& boundarySegmentData(std::size_t id) const
291 {
292 assert(id < numGrids && "Index exceeds number of grids provided");
293 return boundarySegments_[id];
294 }
295
297 std::vector<int>& elementMarkerMap(std::size_t id)
298 {
299 assert(id < numGrids && "Index exceeds number of grids provided");
300 return elementMarkerMaps_[id];
301 }
302
304 std::vector<int>& boundaryMarkerMap(std::size_t id)
305 {
306 assert(id < numGrids && "Index exceeds number of grids provided");
307 return boundaryMarkerMaps_[id];
308 }
309
311 std::unordered_map< GridIndexType, std::vector<GridIndexType> >& embeddedEntityMap(std::size_t id)
312 {
313 assert(id < numGrids && "Index exceeds number of grids provided");
314 return embeddedEntityMaps_[id];
315 }
316
318 std::unordered_map< GridIndexType, std::vector<GridIndexType> >& adjoinedEntityMap(std::size_t id)
319 {
320 assert(id < numGrids && "Index exceeds number of grids provided");
321 return adjoinedEntityMaps_[id];
322 }
323
324private:
326 template<class T>
327 T convertString(const std::string& string) const
328 {
329 T value;
330 std::istringstream stream(string);
331 stream >> value;
332 if (stream.fail())
333 DUNE_THROW(Dune::InvalidStateException, "Couldn't convert string: " << string << "to type: " << typeid(T).name());
334 return value;
335 }
336
338 Dune::GeometryType obtainGeometryType(std::size_t gmshElemType) const
339 {
340 // TODO: Version check with Dune 2.5!
341 switch (gmshElemType)
342 {
343 case 15: return Dune::GeometryTypes::vertex; // points
344 case 1: return Dune::GeometryTypes::line; // lines
345 case 2: return Dune::GeometryTypes::triangle; // triangle
346 case 3: return Dune::GeometryTypes::quadrilateral; // quadrilateral
347 case 4: return Dune::GeometryTypes::tetrahedron; // tetrahedron
348 case 5: return Dune::GeometryTypes::hexahedron; // hexahedron
349 default:
350 DUNE_THROW(Dune::NotImplemented, "FacetCoupling gmsh reader for gmsh element type " << gmshElemType);
351 }
352 }
353
355 void reorder(const Dune::GeometryType gt, VertexIndexSet& cornerIndices) const
356 {
357 // triangles, lines & tetrahedra need no reordering
358 if (gt == Dune::GeometryTypes::hexahedron)
359 {
360 using std::swap;
361 assert(cornerIndices.size() == 8);
362 swap(cornerIndices[2], cornerIndices[3]);
363 swap(cornerIndices[6], cornerIndices[7]);
364 }
365 else if (gt == Dune::GeometryTypes::quadrilateral)
366 {
367 using std::swap;
368 assert(cornerIndices.size() == 4);
369 swap(cornerIndices[2], cornerIndices[3]);
370 }
371 }
372
374 void addEmbeddings(const VertexIndexSet& globalCornerIndices,
375 unsigned int gridIdx,
376 std::size_t curElemIdx,
377 const std::vector<GridIndexType>& embeddedVIndices,
378 const std::vector<bool>& embededIndexIsAssigned)
379 {
380 const unsigned int embeddedGridIdx = gridIdx+1;
381
382 // check for embedments only if any of the corners exist in the embedded grid
383 for (auto globalCornerIdx : globalCornerIndices)
384 {
385 if (embededIndexIsAssigned[globalCornerIdx])
386 {
387 for (std::size_t i = 0; i < elementData_[embeddedGridIdx].size(); ++i)
388 {
389 const auto& e = elementData_[embeddedGridIdx][i];
390
391 auto vertIsContained = [&embeddedVIndices, &globalCornerIndices] (auto eCornerIdx)
392 {
393 return std::find(globalCornerIndices.begin(),
394 globalCornerIndices.end(),
395 embeddedVIndices[eCornerIdx]) != globalCornerIndices.end();
396 };
397
398 // if all corners are contained within this element, it is embedded
399 if ( std::all_of(e.cornerIndices.begin(), e.cornerIndices.end(), vertIsContained) )
400 {
401 embeddedEntityMaps_[gridIdx][curElemIdx].push_back(i);
402 adjoinedEntityMaps_[embeddedGridIdx][i].push_back(curElemIdx);
403 }
404 }
405
406 return;
407 }
408 }
409 }
410
412 std::vector<GlobalPosition> gridVertices_;
413 std::array<VertexIndexSet, numGrids> vertexIndices_;
414 std::array<std::vector<ElementData>, numGrids> elementData_;
415 std::array<std::vector<VertexIndexSet>, numGrids> boundarySegments_;
416
418 std::array< std::unordered_map< GridIndexType, std::vector<GridIndexType> >, numGrids > embeddedEntityMaps_;
419 std::array< std::unordered_map< GridIndexType, std::vector<GridIndexType> >, numGrids > adjoinedEntityMaps_;
420
422 std::array< std::vector<int>, numGrids > elementMarkerMaps_;
423 std::array< std::vector<int>, numGrids > boundaryMarkerMaps_;
424};
425
426} // end namespace Dumux
427
428#endif // DUMUX_FACETCOUPLING_GMSH_READER_HH
Defines the index types used for grid and local indices.
make the local view function available whenever we use the grid geometry
Definition: adapt.hh:29
typename GridView::IndexSet::IndexType GridIndex
Definition: indextraits.hh:39
Reads gmsh files where (n-1)-dimensional grids are defined on the faces or edges of n-dimensional gri...
Definition: gmshreader.hh:61
const std::vector< VertexIndexSet > & boundarySegmentData(std::size_t id) const
Returns the vector of read elements for a grid.
Definition: gmshreader.hh:290
std::unordered_map< GridIndexType, std::vector< GridIndexType > > & embeddedEntityMap(std::size_t id)
Returns the maps of the embedded entities.
Definition: gmshreader.hh:311
void read(const std::string &fileName, bool verbose=false)
Definition: gmshreader.hh:84
std::vector< int > & elementMarkerMap(std::size_t id)
Returns the maps of element markers.
Definition: gmshreader.hh:297
void read(const std::string &fileName, std::size_t boundarySegThresh, bool verbose=false)
Reads the data from a given mesh file.
Definition: gmshreader.hh:90
std::unordered_map< GridIndexType, std::vector< GridIndexType > > & adjoinedEntityMap(std::size_t id)
Returns the maps of the embedments.
Definition: gmshreader.hh:318
const std::vector< ElementData > & elementData(std::size_t id) const
Returns the vector of read elements for a grid.
Definition: gmshreader.hh:283
const std::vector< GlobalPosition > & gridVertices() const
Returns the vector with all grid vertices (entire hierarchy)
Definition: gmshreader.hh:272
std::vector< int > & boundaryMarkerMap(std::size_t id)
Returns the maps of domain markers.
Definition: gmshreader.hh:304
VertexIndexSet & vertexIndices(std::size_t id)
Returns a grid's vertex indices.
Definition: gmshreader.hh:276