| Line | Branch | Exec | Source |
|---|---|---|---|
| 1 | // -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- | ||
| 2 | // vi: set et ts=4 sw=4 sts=4: | ||
| 3 | // | ||
| 4 | // SPDX-FileCopyrightInfo: Copyright © DuMux Project contributors, see AUTHORS.md in root folder | ||
| 5 | // SPDX-License-Identifier: GPL-3.0-or-later | ||
| 6 | // | ||
| 7 | /*! | ||
| 8 | * \file | ||
| 9 | * \ingroup FacetCoupling | ||
| 10 | * \copydoc Dumux::FacetCouplingGmshReader | ||
| 11 | */ | ||
| 12 | #ifndef DUMUX_FACETCOUPLING_GMSH_READER_HH | ||
| 13 | #define DUMUX_FACETCOUPLING_GMSH_READER_HH | ||
| 14 | |||
| 15 | #include <algorithm> | ||
| 16 | #include <cassert> | ||
| 17 | #include <fstream> | ||
| 18 | #include <iostream> | ||
| 19 | #include <sstream> | ||
| 20 | #include <typeinfo> | ||
| 21 | #include <unordered_map> | ||
| 22 | |||
| 23 | #include <dune/common/timer.hh> | ||
| 24 | #include <dune/common/fvector.hh> | ||
| 25 | #include <dune/geometry/type.hh> | ||
| 26 | |||
| 27 | #include <dumux/common/indextraits.hh> | ||
| 28 | |||
| 29 | namespace Dumux { | ||
| 30 | |||
| 31 | /*! | ||
| 32 | * \ingroup FacetCoupling | ||
| 33 | * \brief Reads gmsh files where (n-1)-dimensional grids are defined on the faces | ||
| 34 | * or edges of n-dimensional grids. | ||
| 35 | * | ||
| 36 | * \note Lower-dimensional entities appearing in the grid file are interpreted | ||
| 37 | * either as parts of a lower-dimensional grid living on the sub-entities of | ||
| 38 | * the grids with higher dimension, or as boundary segments. Per default, we | ||
| 39 | * consider all entities as part of the lower-dimensional grids. If you want | ||
| 40 | * to specify boundary segments as well, provide a threshold physical entity | ||
| 41 | * index. All entities with physical entity indices below this threshold will | ||
| 42 | * then be interpreted as boundary segments. Use respective physical entity | ||
| 43 | * indexing in your grid file in that case. | ||
| 44 | * | ||
| 45 | * \tparam BulkGrid The type of the highest-dimensional grid in the hierarchy | ||
| 46 | * \tparam numGrids The number of grids to be considered in the hierarchy | ||
| 47 | */ | ||
| 48 | template <class BulkGrid, int numGrids> | ||
| 49 | class FacetCouplingGmshReader | ||
| 50 | { | ||
| 51 | // extract some necessary info from bulk grid | ||
| 52 | static constexpr int bulkDim = BulkGrid::dimension; | ||
| 53 | static constexpr int bulkDimWorld = BulkGrid::dimensionworld; | ||
| 54 | using ctype = typename BulkGrid::ctype; | ||
| 55 | using GridIndexType = typename IndexTraits< typename BulkGrid::LeafGridView >::GridIndex; | ||
| 56 | using GlobalPosition = Dune::FieldVector<ctype, bulkDimWorld>; | ||
| 57 | |||
| 58 | // determine minimum dimension for which a grid is created | ||
| 59 | static constexpr int minGridDim = BulkGrid::dimension - numGrids + 1; | ||
| 60 | static_assert(minGridDim >= 1, "Grids with dim < 1 cannot be read!"); | ||
| 61 | |||
| 62 | // structure to store data on an element | ||
| 63 | using VertexIndexSet = std::vector<GridIndexType>; | ||
| 64 |
2/10✗ Branch 0 not taken.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 237742 times.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✓ Branch 10 taken 164037 times.
✗ Branch 11 not taken.
|
1369374 | struct ElementData |
| 65 | { | ||
| 66 | Dune::GeometryType gt; | ||
| 67 | VertexIndexSet cornerIndices; | ||
| 68 | }; | ||
| 69 | |||
| 70 | public: | ||
| 71 | //! Reads the data from a given mesh file | ||
| 72 | //! Use this routine if you don't specify boundary segments in the grid file | ||
| 73 | void read(const std::string& fileName, bool verbose = false) | ||
| 74 | { | ||
| 75 | read(fileName, 0, verbose); | ||
| 76 | } | ||
| 77 | |||
| 78 | //! Reads the data from a given mesh file | ||
| 79 | 74 | void read(const std::string& fileName, std::size_t boundarySegThresh, bool verbose = false) | |
| 80 | { | ||
| 81 | 74 | Dune::Timer watch; | |
| 82 |
2/2✓ Branch 0 taken 8 times.
✓ Branch 1 taken 66 times.
|
82 | if (verbose) std::cout << "Opening " << fileName << std::endl; |
| 83 | 74 | std::ifstream gridFile(fileName); | |
| 84 |
2/4✗ Branch 0 not taken.
✓ Branch 1 taken 74 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 74 times.
|
148 | if (gridFile.fail()) |
| 85 | ✗ | DUNE_THROW(Dune::InvalidStateException, "Could not open the given .msh file. Make sure it exists"); | |
| 86 | |||
| 87 | // currently we only support version 2 file format | ||
| 88 |
1/4✓ Branch 1 taken 74 times.
✗ Branch 2 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
|
148 | std::string line; |
| 89 |
1/2✓ Branch 1 taken 74 times.
✗ Branch 2 not taken.
|
74 | std::getline(gridFile, line); |
| 90 |
1/2✗ Branch 1 not taken.
✓ Branch 2 taken 74 times.
|
148 | if (line.find("$MeshFormat") == std::string::npos) |
| 91 | ✗ | DUNE_THROW(Dune::InvalidStateException, "Expected $MeshFormat in the first line of the grid file!"); | |
| 92 |
1/2✓ Branch 1 taken 74 times.
✗ Branch 2 not taken.
|
74 | std::getline(gridFile, line); |
| 93 |
1/2✗ Branch 1 not taken.
✓ Branch 2 taken 74 times.
|
148 | if (line.find_first_of("2") != 0) |
| 94 | ✗ | DUNE_THROW(Dune::InvalidStateException, "Currently only Gmsh mesh file format version 2 is supported!"); | |
| 95 | |||
| 96 | // read file until we get to the list of nodes | ||
| 97 |
2/2✓ Branch 1 taken 148 times.
✓ Branch 2 taken 74 times.
|
444 | while (line.find("$Nodes") == std::string::npos) |
| 98 |
1/2✓ Branch 1 taken 148 times.
✗ Branch 2 not taken.
|
148 | std::getline(gridFile, line); |
| 99 | |||
| 100 | // read all vertices | ||
| 101 |
1/2✓ Branch 1 taken 74 times.
✗ Branch 2 not taken.
|
74 | std::getline(gridFile, line); |
| 102 |
1/2✓ Branch 1 taken 74 times.
✗ Branch 2 not taken.
|
74 | const auto numVertices = convertString<std::size_t>(line); |
| 103 |
1/2✓ Branch 1 taken 74 times.
✗ Branch 2 not taken.
|
74 | gridVertices_.resize(numVertices); |
| 104 | |||
| 105 |
1/2✓ Branch 1 taken 74 times.
✗ Branch 2 not taken.
|
74 | std::getline(gridFile, line); |
| 106 | std::size_t vertexCount = 0; | ||
| 107 |
2/2✓ Branch 1 taken 147998 times.
✓ Branch 2 taken 74 times.
|
296144 | while (line.find("$EndNodes") == std::string::npos) |
| 108 | { | ||
| 109 | // drop first entry in line (vertex index) and read coordinates | ||
| 110 |
1/2✓ Branch 1 taken 147998 times.
✗ Branch 2 not taken.
|
147998 | std::istringstream stream(line); |
| 111 |
3/8✓ Branch 1 taken 147998 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 147998 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 147998 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
|
443994 | std::string buf; stream >> buf; |
| 112 | 147998 | GlobalPosition v; | |
| 113 |
3/4✓ Branch 0 taken 147998 times.
✓ Branch 1 taken 319143 times.
✓ Branch 3 taken 319143 times.
✗ Branch 4 not taken.
|
615139 | for (auto& coord : v) |
| 114 | { | ||
| 115 |
1/2✓ Branch 1 taken 319143 times.
✗ Branch 2 not taken.
|
319143 | stream >> coord; |
| 116 |
2/22✗ Branch 0 not taken.
✓ Branch 1 taken 319143 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 319143 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✗ Branch 17 not taken.
✗ Branch 18 not taken.
✗ Branch 21 not taken.
✗ Branch 22 not taken.
✗ Branch 24 not taken.
✗ Branch 25 not taken.
✗ Branch 27 not taken.
✗ Branch 28 not taken.
✗ Branch 29 not taken.
✗ Branch 30 not taken.
✗ Branch 32 not taken.
✗ Branch 33 not taken.
✗ Branch 35 not taken.
✗ Branch 36 not taken.
|
638286 | if (stream.fail()) DUNE_THROW(Dune::IOError, "Could not read vertex coordinate"); |
| 117 | } | ||
| 118 | |||
| 119 | // insert vertex into container and move to next line | ||
| 120 |
1/2✓ Branch 1 taken 147998 times.
✗ Branch 2 not taken.
|
147998 | gridVertices_[vertexCount++] = v; |
| 121 |
1/2✓ Branch 1 taken 147998 times.
✗ Branch 2 not taken.
|
147998 | std::getline(gridFile, line); |
| 122 | } | ||
| 123 | |||
| 124 | // we should always find as many vertices as the mesh file states | ||
| 125 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 74 times.
|
74 | if (vertexCount != numVertices) |
| 126 | ✗ | DUNE_THROW(Dune::InvalidStateException, "Couldn't find as many vertices as stated in the .msh file"); | |
| 127 | |||
| 128 | // read file until we get to the list of elements | ||
| 129 |
2/2✓ Branch 1 taken 74 times.
✓ Branch 2 taken 74 times.
|
296 | while(line.find("$Elements") == std::string::npos) |
| 130 |
1/2✓ Branch 1 taken 74 times.
✗ Branch 2 not taken.
|
74 | std::getline(gridFile, line); |
| 131 | |||
| 132 | // read elements | ||
| 133 |
1/2✓ Branch 1 taken 74 times.
✗ Branch 2 not taken.
|
74 | std::getline(gridFile, line); |
| 134 |
1/2✓ Branch 1 taken 74 times.
✗ Branch 2 not taken.
|
74 | const auto numElements = convertString<std::size_t>(line); |
| 135 | |||
| 136 | // keep track of number of elements | ||
| 137 | std::array<std::size_t, numGrids> elementCount; | ||
| 138 | std::fill(elementCount.begin(), elementCount.end(), 0); | ||
| 139 | |||
| 140 | // Construct maps that map bulk grid vertex | ||
| 141 | // indices to lowDim vertex indices. -1 indicates non-initialized status | ||
| 142 | std::size_t elemCount = 0; | ||
| 143 | std::array<std::size_t, numGrids> gridVertexCount; | ||
| 144 |
1/2✓ Branch 0 taken 74 times.
✗ Branch 1 not taken.
|
74 | std::array<std::vector<GridIndexType>, numGrids> gridVertexMap; |
| 145 | 74 | std::array<std::vector<bool>, numGrids> idxIsAssigned; | |
| 146 | 74 | std::fill(gridVertexCount.begin(), gridVertexCount.end(), 0); | |
| 147 |
4/10✓ Branch 1 taken 74 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 74 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 74 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 74 times.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
|
296 | std::fill(gridVertexMap.begin(), gridVertexMap.end(), std::vector<GridIndexType>(vertexCount)); |
| 148 |
3/8✓ Branch 1 taken 74 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 74 times.
✗ Branch 5 not taken.
✓ Branch 8 taken 74 times.
✗ Branch 9 not taken.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
|
222 | std::fill(idxIsAssigned.begin(), idxIsAssigned.end(), std::vector<bool>(vertexCount, false)); |
| 149 |
1/2✓ Branch 1 taken 74 times.
✗ Branch 2 not taken.
|
74 | std::getline(gridFile, line); |
| 150 |
2/2✓ Branch 1 taken 165959 times.
✓ Branch 2 taken 74 times.
|
332066 | while (line.find("$EndElements") == std::string::npos) |
| 151 | { | ||
| 152 | // pass all indices into vector | ||
| 153 |
1/2✓ Branch 1 taken 165959 times.
✗ Branch 2 not taken.
|
165959 | std::istringstream stream(line); |
| 154 |
0/2✗ Branch 1 not taken.
✗ Branch 2 not taken.
|
331918 | std::string buf; |
| 155 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 165959 times.
|
331918 | std::vector<std::size_t> lineData; |
| 156 |
7/12✓ Branch 1 taken 1633926 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 1467967 times.
✓ Branch 4 taken 165959 times.
✓ Branch 5 taken 1467967 times.
✓ Branch 6 taken 165959 times.
✓ Branch 8 taken 1467967 times.
✗ Branch 9 not taken.
✓ Branch 11 taken 1467967 times.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
|
1633926 | while (stream >> buf) lineData.push_back(convertString<std::size_t>(buf)); |
| 157 |
2/4✗ Branch 0 not taken.
✓ Branch 1 taken 165959 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 165959 times.
|
331918 | assert(lineData.size() >= 4 && "Grid format erroneous or unsupported"); |
| 158 | |||
| 159 | // obtain geometry type | ||
| 160 |
1/2✓ Branch 1 taken 165959 times.
✗ Branch 2 not taken.
|
165959 | const auto gt = obtainGeometryType( lineData[1] ); |
| 161 | 165959 | const std::size_t physicalIndex = lineData[3]; | |
| 162 |
2/2✓ Branch 0 taken 7013 times.
✓ Branch 1 taken 158946 times.
|
165959 | const auto geoDim = gt.dim(); |
| 163 | 165959 | const bool isBoundarySeg = geoDim != bulkDim && physicalIndex < boundarySegThresh; | |
| 164 |
1/2✓ Branch 0 taken 5093 times.
✗ Branch 1 not taken.
|
5093 | if (geoDim >= minGridDim-1) |
| 165 | { | ||
| 166 | // insert boundary segment | ||
| 167 |
2/2✓ Branch 0 taken 1922 times.
✓ Branch 1 taken 164037 times.
|
165959 | if ((isBoundarySeg || geoDim == minGridDim-1)) |
| 168 | { | ||
| 169 | 1922 | const unsigned int nextLevelGridIdx = bulkDim-geoDim-1; | |
| 170 | |||
| 171 |
1/2✓ Branch 0 taken 1922 times.
✗ Branch 1 not taken.
|
5766 | VertexIndexSet corners; |
| 172 | 9610 | auto it = lineData.begin()+2+lineData[2]+1; | |
| 173 |
6/6✓ Branch 0 taken 5764 times.
✓ Branch 1 taken 1922 times.
✓ Branch 2 taken 5764 times.
✓ Branch 3 taken 1922 times.
✓ Branch 4 taken 5764 times.
✓ Branch 5 taken 1922 times.
|
11530 | for (; it != lineData.end(); ++it) |
| 174 | { | ||
| 175 | 5764 | *it -= 1; // gmsh indices start from 1 | |
| 176 | |||
| 177 | // insert map if vertex is not inserted yet | ||
| 178 |
8/8✓ Branch 0 taken 979 times.
✓ Branch 1 taken 4785 times.
✓ Branch 2 taken 979 times.
✓ Branch 3 taken 4785 times.
✓ Branch 4 taken 979 times.
✓ Branch 5 taken 4785 times.
✓ Branch 6 taken 979 times.
✓ Branch 7 taken 4785 times.
|
23056 | if (!idxIsAssigned[nextLevelGridIdx][*it]) |
| 179 | { | ||
| 180 |
3/6✓ Branch 1 taken 979 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 979 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 979 times.
✗ Branch 8 not taken.
|
2937 | gridVertexMap[nextLevelGridIdx][*it] = gridVertexCount[nextLevelGridIdx]++; |
| 181 |
3/6✓ Branch 1 taken 979 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 979 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 979 times.
✗ Branch 8 not taken.
|
2937 | idxIsAssigned[nextLevelGridIdx][*it] = true; |
| 182 |
2/4✓ Branch 1 taken 979 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 979 times.
✗ Branch 5 not taken.
|
1958 | vertexIndices_[nextLevelGridIdx].push_back(*it); |
| 183 | } | ||
| 184 | |||
| 185 |
3/6✓ Branch 1 taken 5764 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 5764 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 5764 times.
✗ Branch 8 not taken.
|
17292 | corners.push_back(gridVertexMap[nextLevelGridIdx][*it]); |
| 186 | } | ||
| 187 | |||
| 188 | // marker = physical entity index | ||
| 189 |
2/6✓ Branch 1 taken 1922 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1922 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
|
3844 | boundaryMarkerMaps_[nextLevelGridIdx].push_back(physicalIndex); |
| 190 |
2/4✓ Branch 1 taken 1922 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1922 times.
✗ Branch 5 not taken.
|
3844 | boundarySegments_[nextLevelGridIdx].push_back(corners); |
| 191 | } | ||
| 192 | |||
| 193 | // insert element | ||
| 194 | else | ||
| 195 | { | ||
| 196 | 164037 | const unsigned int gridIdx = bulkDim-geoDim; | |
| 197 | |||
| 198 | 328074 | VertexIndexSet corners; | |
| 199 | 820185 | auto it = lineData.begin()+2+lineData[2]+1; | |
| 200 |
6/6✓ Branch 0 taken 632408 times.
✓ Branch 1 taken 164037 times.
✓ Branch 2 taken 632408 times.
✓ Branch 3 taken 164037 times.
✓ Branch 4 taken 632408 times.
✓ Branch 5 taken 164037 times.
|
1124519 | for (; it != lineData.end(); ++it) |
| 201 | { | ||
| 202 | 632408 | *it -= 1; // gmsh indices start from 1 | |
| 203 | |||
| 204 | // insert map if vertex is not inserted yet | ||
| 205 |
8/8✓ Branch 0 taken 149909 times.
✓ Branch 1 taken 482499 times.
✓ Branch 2 taken 149909 times.
✓ Branch 3 taken 482499 times.
✓ Branch 4 taken 149909 times.
✓ Branch 5 taken 482499 times.
✓ Branch 6 taken 149909 times.
✓ Branch 7 taken 482499 times.
|
2529632 | if (!idxIsAssigned[gridIdx][*it]) |
| 206 | { | ||
| 207 |
3/6✓ Branch 1 taken 149909 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 149909 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 149909 times.
✗ Branch 8 not taken.
|
449727 | gridVertexMap[gridIdx][*it] = gridVertexCount[gridIdx]++; |
| 208 |
3/6✓ Branch 1 taken 149909 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 149909 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 149909 times.
✗ Branch 8 not taken.
|
449727 | idxIsAssigned[gridIdx][*it] = true; |
| 209 |
2/4✓ Branch 1 taken 149909 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 149909 times.
✗ Branch 5 not taken.
|
299818 | vertexIndices_[gridIdx].push_back(*it); |
| 210 | } | ||
| 211 | |||
| 212 |
3/6✓ Branch 1 taken 632408 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 632408 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 632408 times.
✗ Branch 8 not taken.
|
1897224 | corners.push_back(gridVertexMap[gridIdx][*it]); |
| 213 | } | ||
| 214 | |||
| 215 | // add data to embedments/embeddings | ||
| 216 |
2/2✓ Branch 0 taken 161425 times.
✓ Branch 1 taken 2612 times.
|
164037 | if (geoDim > minGridDim) |
| 217 | { | ||
| 218 |
2/4✓ Branch 1 taken 161425 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 161425 times.
✗ Branch 5 not taken.
|
322850 | const auto gridElemCount = elementData_[gridIdx].size(); |
| 219 |
1/2✓ Branch 1 taken 161425 times.
✗ Branch 2 not taken.
|
161425 | const auto& embeddedVIndices = vertexIndices_[gridIdx+1]; |
| 220 |
1/2✓ Branch 1 taken 161425 times.
✗ Branch 2 not taken.
|
161425 | const auto& embeddedIndicesAssigned = idxIsAssigned[gridIdx+1]; |
| 221 | |||
| 222 |
4/10✓ Branch 1 taken 161425 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 161425 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 161425 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 161425 times.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
|
807125 | VertexIndexSet cornerIndicesGlobal(corners.size()); |
| 223 |
4/4✓ Branch 0 taken 626932 times.
✓ Branch 1 taken 161425 times.
✓ Branch 2 taken 626932 times.
✓ Branch 3 taken 161425 times.
|
1415289 | for (unsigned int i = 0; i < corners.size(); ++i) |
| 224 | 3134660 | cornerIndicesGlobal[i] = vertexIndices_[gridIdx][corners[i]]; | |
| 225 |
1/2✓ Branch 1 taken 161425 times.
✗ Branch 2 not taken.
|
161425 | addEmbeddings(cornerIndicesGlobal, gridIdx, gridElemCount, embeddedVIndices, embeddedIndicesAssigned); |
| 226 | } | ||
| 227 | |||
| 228 | // ensure dune-specific corner ordering | ||
| 229 | 164037 | reorder(gt, corners); | |
| 230 | |||
| 231 | // insert element data to grid's container | ||
| 232 |
2/4✓ Branch 1 taken 164037 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 164037 times.
✗ Branch 5 not taken.
|
328074 | elementMarkerMaps_[gridIdx].push_back(physicalIndex); |
| 233 |
5/14✓ Branch 1 taken 164037 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 164037 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 164037 times.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✓ Branch 10 taken 164037 times.
✓ Branch 11 taken 164037 times.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
✗ Branch 15 not taken.
✗ Branch 16 not taken.
|
328074 | elementData_[gridIdx].emplace_back(ElementData({gt, corners})); |
| 234 | } | ||
| 235 | } | ||
| 236 | |||
| 237 | // get next line | ||
| 238 |
1/2✓ Branch 1 taken 165959 times.
✗ Branch 2 not taken.
|
165959 | std::getline(gridFile, line); |
| 239 |
1/2✓ Branch 0 taken 165959 times.
✗ Branch 1 not taken.
|
165959 | elemCount++; |
| 240 | } | ||
| 241 | |||
| 242 | // make sure we read all elements | ||
| 243 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 74 times.
|
74 | if (elemCount != numElements) |
| 244 | ✗ | DUNE_THROW(Dune::InvalidStateException, "Didn't read as many elements as stated in the .msh file"); | |
| 245 | |||
| 246 |
2/2✓ Branch 0 taken 8 times.
✓ Branch 1 taken 66 times.
|
74 | if (verbose) |
| 247 | { | ||
| 248 |
1/2✓ Branch 2 taken 8 times.
✗ Branch 3 not taken.
|
16 | std::cout << "Finished reading gmsh file" << std::endl; |
| 249 |
2/2✓ Branch 0 taken 24 times.
✓ Branch 1 taken 8 times.
|
32 | for (std::size_t id = 0; id < numGrids; ++id) |
| 250 | { | ||
| 251 |
3/6✓ Branch 1 taken 24 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 24 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 24 times.
✗ Branch 8 not taken.
|
72 | std::cout << elementData_[id].size() << " " |
| 252 |
1/2✓ Branch 1 taken 24 times.
✗ Branch 2 not taken.
|
24 | << bulkDim-id << "-dimensional elements comprising of " |
| 253 |
3/6✓ Branch 1 taken 24 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 24 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 24 times.
✗ Branch 8 not taken.
|
48 | << gridVertexCount[id] << " vertices"; |
| 254 |
3/4✓ Branch 0 taken 16 times.
✓ Branch 1 taken 8 times.
✓ Branch 4 taken 16 times.
✗ Branch 5 not taken.
|
40 | if (id < numGrids-1) std::cout << "," << std::endl; |
| 255 | } | ||
| 256 |
2/4✓ Branch 3 taken 8 times.
✗ Branch 4 not taken.
✓ Branch 7 taken 8 times.
✗ Branch 8 not taken.
|
24 | std::cout << " have been read in " << watch.elapsed() << " seconds." << std::endl; |
| 257 | } | ||
| 258 | 74 | } | |
| 259 | |||
| 260 | //! Returns the vector with all grid vertices (entire hierarchy) | ||
| 261 | const std::vector<GlobalPosition>& gridVertices() const | ||
| 262 |
1/2✓ Branch 1 taken 64 times.
✗ Branch 2 not taken.
|
64 | { return gridVertices_; } |
| 263 | |||
| 264 | //! Returns a grid's vertex indices | ||
| 265 | VertexIndexSet& vertexIndices(std::size_t id) | ||
| 266 | { | ||
| 267 | assert(id < numGrids && "Index exceeds number of grids provided"); | ||
| 268 | 632 | return vertexIndices_[id]; | |
| 269 | } | ||
| 270 | |||
| 271 | //! Returns the vector of read elements for a grid | ||
| 272 | const std::vector<ElementData>& elementData(std::size_t id) const | ||
| 273 | { | ||
| 274 | assert(id < numGrids && "Index exceeds number of grids provided"); | ||
| 275 | 316 | return elementData_[id]; | |
| 276 | } | ||
| 277 | |||
| 278 | //! Returns the vector of read elements for a grid | ||
| 279 | const std::vector<VertexIndexSet>& boundarySegmentData(std::size_t id) const | ||
| 280 | { | ||
| 281 | assert(id < numGrids && "Index exceeds number of grids provided"); | ||
| 282 | 48 | return boundarySegments_[id]; | |
| 283 | } | ||
| 284 | |||
| 285 | //! Returns the maps of element markers | ||
| 286 | std::vector<int>& elementMarkerMap(std::size_t id) | ||
| 287 | { | ||
| 288 | assert(id < numGrids && "Index exceeds number of grids provided"); | ||
| 289 | 48 | return elementMarkerMaps_[id]; | |
| 290 | } | ||
| 291 | |||
| 292 | //! Returns the maps of domain markers | ||
| 293 | std::vector<int>& boundaryMarkerMap(std::size_t id) | ||
| 294 | { | ||
| 295 | assert(id < numGrids && "Index exceeds number of grids provided"); | ||
| 296 | 48 | return boundaryMarkerMaps_[id]; | |
| 297 | } | ||
| 298 | |||
| 299 | //! Returns the maps of the embedded entities | ||
| 300 | std::unordered_map< GridIndexType, std::vector<GridIndexType> >& embeddedEntityMap(std::size_t id) | ||
| 301 | { | ||
| 302 | assert(id < numGrids && "Index exceeds number of grids provided"); | ||
| 303 | 316 | return embeddedEntityMaps_[id]; | |
| 304 | } | ||
| 305 | |||
| 306 | //! Returns the maps of the embedments | ||
| 307 | std::unordered_map< GridIndexType, std::vector<GridIndexType> >& adjoinedEntityMap(std::size_t id) | ||
| 308 | { | ||
| 309 | assert(id < numGrids && "Index exceeds number of grids provided"); | ||
| 310 | 316 | return adjoinedEntityMaps_[id]; | |
| 311 | } | ||
| 312 | |||
| 313 | private: | ||
| 314 | //! converts a value contained in a string | ||
| 315 | template<class T> | ||
| 316 | 1468115 | T convertString(const std::string& string) const | |
| 317 | { | ||
| 318 | T value; | ||
| 319 | 1468115 | std::istringstream stream(string); | |
| 320 |
1/2✓ Branch 1 taken 1468115 times.
✗ Branch 2 not taken.
|
1468115 | stream >> value; |
| 321 |
2/4✗ Branch 0 not taken.
✓ Branch 1 taken 1468115 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 1468115 times.
|
2936230 | if (stream.fail()) |
| 322 | ✗ | DUNE_THROW(Dune::InvalidStateException, "Couldn't convert string: " << string << "to type: " << typeid(T).name()); | |
| 323 | 1468115 | return value; | |
| 324 | } | ||
| 325 | |||
| 326 | //! obtain Dune::GeometryType from a given gmsh element type | ||
| 327 | 165959 | Dune::GeometryType obtainGeometryType(std::size_t gmshElemType) const | |
| 328 | { | ||
| 329 |
4/7✗ Branch 0 not taken.
✓ Branch 1 taken 2362 times.
✓ Branch 2 taken 20940 times.
✓ Branch 3 taken 130036 times.
✓ Branch 4 taken 12621 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
|
165959 | switch (gmshElemType) |
| 330 | { | ||
| 331 | ✗ | case 15: return Dune::GeometryTypes::vertex; // points | |
| 332 | 2362 | case 1: return Dune::GeometryTypes::line; // lines | |
| 333 | 20940 | case 2: return Dune::GeometryTypes::triangle; // triangle | |
| 334 | 130036 | case 3: return Dune::GeometryTypes::quadrilateral; // quadrilateral | |
| 335 | 12621 | case 4: return Dune::GeometryTypes::tetrahedron; // tetrahedron | |
| 336 | ✗ | case 5: return Dune::GeometryTypes::hexahedron; // hexahedron | |
| 337 | ✗ | default: | |
| 338 | ✗ | DUNE_THROW(Dune::NotImplemented, "FacetCoupling gmsh reader for gmsh element type " << gmshElemType); | |
| 339 | } | ||
| 340 | } | ||
| 341 | |||
| 342 | //! reorders in a dune way a set of given element corners in gmsh ordering | ||
| 343 | 164037 | void reorder(const Dune::GeometryType gt, VertexIndexSet& cornerIndices) const | |
| 344 | { | ||
| 345 | // triangles, lines & tetrahedra need no reordering | ||
| 346 |
1/2✓ Branch 0 taken 164037 times.
✗ Branch 1 not taken.
|
164037 | if (gt == Dune::GeometryTypes::hexahedron) |
| 347 | { | ||
| 348 | using std::swap; | ||
| 349 | ✗ | assert(cornerIndices.size() == 8); | |
| 350 | ✗ | swap(cornerIndices[2], cornerIndices[3]); | |
| 351 | ✗ | swap(cornerIndices[6], cornerIndices[7]); | |
| 352 | } | ||
| 353 |
1/2✓ Branch 0 taken 164037 times.
✗ Branch 1 not taken.
|
164037 | else if (gt == Dune::GeometryTypes::quadrilateral) |
| 354 | { | ||
| 355 | using std::swap; | ||
| 356 |
2/4✗ Branch 0 not taken.
✓ Branch 1 taken 130036 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 130036 times.
|
260072 | assert(cornerIndices.size() == 4); |
| 357 | 390108 | swap(cornerIndices[2], cornerIndices[3]); | |
| 358 | } | ||
| 359 | 164037 | } | |
| 360 | |||
| 361 | //! adds embeddings/embedments to the map for a given element | ||
| 362 | 161425 | void addEmbeddings(const VertexIndexSet& globalCornerIndices, | |
| 363 | unsigned int gridIdx, | ||
| 364 | std::size_t curElemIdx, | ||
| 365 | const std::vector<GridIndexType>& embeddedVIndices, | ||
| 366 | const std::vector<bool>& embededIndexIsAssigned) | ||
| 367 | { | ||
| 368 | 161425 | const unsigned int embeddedGridIdx = gridIdx+1; | |
| 369 | |||
| 370 | // check for embedments only if any of the corners exist in the embedded grid | ||
| 371 |
4/4✓ Branch 0 taken 604009 times.
✓ Branch 1 taken 150510 times.
✓ Branch 2 taken 604009 times.
✓ Branch 3 taken 150510 times.
|
1077369 | for (auto globalCornerIdx : globalCornerIndices) |
| 372 | { | ||
| 373 |
4/4✓ Branch 0 taken 10915 times.
✓ Branch 1 taken 593094 times.
✓ Branch 2 taken 10915 times.
✓ Branch 3 taken 593094 times.
|
1208018 | if (embededIndexIsAssigned[globalCornerIdx]) |
| 374 | { | ||
| 375 |
6/6✓ Branch 0 taken 715524 times.
✓ Branch 1 taken 10915 times.
✓ Branch 2 taken 715524 times.
✓ Branch 3 taken 10915 times.
✓ Branch 4 taken 715524 times.
✓ Branch 5 taken 10915 times.
|
748269 | for (std::size_t i = 0; i < elementData_[embeddedGridIdx].size(); ++i) |
| 376 | { | ||
| 377 | 1431048 | const auto& e = elementData_[embeddedGridIdx][i]; | |
| 378 | |||
| 379 | 2204140 | auto vertIsContained = [&embeddedVIndices, &globalCornerIndices] (auto eCornerIdx) | |
| 380 | { | ||
| 381 | 1488616 | return std::find(globalCornerIndices.begin(), | |
| 382 | globalCornerIndices.end(), | ||
| 383 | 2232924 | embeddedVIndices[eCornerIdx]) != globalCornerIndices.end(); | |
| 384 | }; | ||
| 385 | |||
| 386 | // if all corners are contained within this element, it is embedded | ||
| 387 |
2/2✓ Branch 3 taken 6044 times.
✓ Branch 4 taken 709480 times.
|
2146572 | if ( std::all_of(e.cornerIndices.begin(), e.cornerIndices.end(), vertIsContained) ) |
| 388 | { | ||
| 389 |
3/6✓ Branch 1 taken 6044 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 6044 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 6044 times.
✗ Branch 8 not taken.
|
12088 | embeddedEntityMaps_[gridIdx][curElemIdx].push_back(i); |
| 390 |
3/6✓ Branch 1 taken 6044 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 6044 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 6044 times.
✗ Branch 8 not taken.
|
12088 | adjoinedEntityMaps_[embeddedGridIdx][i].push_back(curElemIdx); |
| 391 | } | ||
| 392 | } | ||
| 393 | |||
| 394 | return; | ||
| 395 | } | ||
| 396 | } | ||
| 397 | } | ||
| 398 | |||
| 399 | //! data on grid entities | ||
| 400 | std::vector<GlobalPosition> gridVertices_; | ||
| 401 | std::array<VertexIndexSet, numGrids> vertexIndices_; | ||
| 402 | std::array<std::vector<ElementData>, numGrids> elementData_; | ||
| 403 | std::array<std::vector<VertexIndexSet>, numGrids> boundarySegments_; | ||
| 404 | |||
| 405 | //! data on connectivity between the grids | ||
| 406 | std::array< std::unordered_map< GridIndexType, std::vector<GridIndexType> >, numGrids > embeddedEntityMaps_; | ||
| 407 | std::array< std::unordered_map< GridIndexType, std::vector<GridIndexType> >, numGrids > adjoinedEntityMaps_; | ||
| 408 | |||
| 409 | //! data on domain and boundary markers | ||
| 410 | std::array< std::vector<int>, numGrids > elementMarkerMaps_; | ||
| 411 | std::array< std::vector<int>, numGrids > boundaryMarkerMaps_; | ||
| 412 | }; | ||
| 413 | |||
| 414 | } // end namespace Dumux | ||
| 415 | |||
| 416 | #endif | ||
| 417 |