| 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 CCMpfaDiscretization | ||
| 10 | * \brief Helper class to get data required for mpfa scheme. | ||
| 11 | */ | ||
| 12 | #ifndef DUMUX_DISCRETIZATION_CC_MPFA_HELPER_HH | ||
| 13 | #define DUMUX_DISCRETIZATION_CC_MPFA_HELPER_HH | ||
| 14 | |||
| 15 | #include <dune/common/exceptions.hh> | ||
| 16 | #include <dune/common/fvector.hh> | ||
| 17 | #include <dune/common/fmatrix.hh> | ||
| 18 | #include <dune/common/parallel/mpihelper.hh> | ||
| 19 | #include <dune/geometry/type.hh> | ||
| 20 | |||
| 21 | #include <dumux/common/math.hh> | ||
| 22 | |||
| 23 | namespace Dumux { | ||
| 24 | |||
| 25 | /*! | ||
| 26 | * \ingroup CCMpfaDiscretization | ||
| 27 | * \brief Dimension-specific helper class to get data required for mpfa scheme. | ||
| 28 | */ | ||
| 29 | template<class GridGeometry, int dim, int dimWorld> | ||
| 30 | class MpfaDimensionHelper; | ||
| 31 | |||
| 32 | /*! | ||
| 33 | * \ingroup CCMpfaDiscretization | ||
| 34 | * \brief Dimension-specific mpfa helper class for dim == 2 & dimWorld == 2 | ||
| 35 | */ | ||
| 36 | template<class GridGeometry> | ||
| 37 | class MpfaDimensionHelper<GridGeometry, /*dim*/2, /*dimWorld*/2> | ||
| 38 | { | ||
| 39 | using GridView = typename GridGeometry::GridView; | ||
| 40 | using CoordScalar = typename GridView::ctype; | ||
| 41 | using Element = typename GridView::template Codim<0>::Entity; | ||
| 42 | |||
| 43 | using GlobalPosition = typename Element::Geometry::GlobalCoordinate; | ||
| 44 | |||
| 45 | using SubControlVolumeFace = typename GridGeometry::SubControlVolumeFace; | ||
| 46 | using ScvfCornerVector = typename SubControlVolumeFace::Traits::CornerStorage; | ||
| 47 | |||
| 48 | // Container to store the positions of intersections required for scvf | ||
| 49 | // corner computation. In 2d, these are the center plus the two corners | ||
| 50 | using ScvfPositionsOnIntersection = std::array<GlobalPosition, 3>; | ||
| 51 | |||
| 52 | public: | ||
| 53 | /*! | ||
| 54 | * \brief Calculates the inner normal vectors to a given scv basis. | ||
| 55 | * \param scvBasis The basis of an scv | ||
| 56 | */ | ||
| 57 | template< class ScvBasis > | ||
| 58 | 34663321 | static ScvBasis calculateInnerNormals(const ScvBasis& scvBasis) | |
| 59 | { | ||
| 60 |
4/4✓ Branch 0 taken 38 times.
✓ Branch 1 taken 34663283 times.
✓ Branch 3 taken 37 times.
✓ Branch 4 taken 1 times.
|
34663321 | static const Dune::FieldMatrix<CoordScalar, 2, 2> R = {{0.0, 1.0}, {-1.0, 0.0}}; |
| 61 | |||
| 62 | 34663321 | ScvBasis innerNormals; | |
| 63 | 103989963 | R.mv(scvBasis[1], innerNormals[0]); | |
| 64 | 34663321 | R.mv(scvBasis[0], innerNormals[1]); | |
| 65 | |||
| 66 | // adjust sign depending on basis being a RHS | ||
| 67 | 69326642 | if (!isRightHandSystem(scvBasis)) | |
| 68 | innerNormals[0] *= -1.0; | ||
| 69 | else | ||
| 70 | innerNormals[1] *= -1.0; | ||
| 71 | |||
| 72 | 34663321 | return innerNormals; | |
| 73 | } | ||
| 74 | |||
| 75 | /*! | ||
| 76 | * \brief Calculates the determinant of an scv basis. | ||
| 77 | * This is equal to the cross product for dim = dimWorld = 2 | ||
| 78 | * \param scvBasis The basis of an scv | ||
| 79 | */ | ||
| 80 | template< class ScvBasis > | ||
| 81 | static CoordScalar calculateDetX(const ScvBasis& scvBasis) | ||
| 82 | { | ||
| 83 | using std::abs; | ||
| 84 | 173316605 | return abs(crossProduct<CoordScalar>(scvBasis[0], scvBasis[1])); | |
| 85 | } | ||
| 86 | |||
| 87 | /*! | ||
| 88 | * \brief Returns the global number of scvfs in the grid. Assumes the grid to be made up of only | ||
| 89 | * basic geometry types. Overlad this function if you want to use different geometry types. | ||
| 90 | * \param gridView The grid view to be checked | ||
| 91 | */ | ||
| 92 | 17 | static std::size_t getGlobalNumScvf(const GridView& gridView) | |
| 93 | { | ||
| 94 |
1/2✗ Branch 3 not taken.
✓ Branch 4 taken 17 times.
|
17 | assert(gridView.size(Dune::GeometryTypes::triangle) |
| 95 | + gridView.size(Dune::GeometryTypes::quadrilateral) == gridView.size(0)); | ||
| 96 | |||
| 97 | return gridView.size(Dune::GeometryTypes::triangle) | ||
| 98 | 17 | * getNumLocalScvfs(Dune::GeometryTypes::triangle) | |
| 99 | + gridView.size(Dune::GeometryTypes::quadrilateral) | ||
| 100 | 17 | * getNumLocalScvfs(Dune::GeometryTypes::quadrilateral); | |
| 101 | } | ||
| 102 | |||
| 103 | /*! | ||
| 104 | * \brief Checks whether or not a given scv basis forms a right hand system. | ||
| 105 | * \param scvBasis The basis of an scv | ||
| 106 | */ | ||
| 107 | template< class ScvBasis > | ||
| 108 | static bool isRightHandSystem(const ScvBasis& scvBasis) | ||
| 109 |
20/20✓ Branch 0 taken 13628205 times.
✓ Branch 1 taken 21035116 times.
✓ Branch 2 taken 13628205 times.
✓ Branch 3 taken 21035116 times.
✓ Branch 4 taken 13628205 times.
✓ Branch 5 taken 21035116 times.
✓ Branch 6 taken 13628205 times.
✓ Branch 7 taken 21035116 times.
✓ Branch 8 taken 13628205 times.
✓ Branch 9 taken 21035116 times.
✓ Branch 10 taken 200 times.
✓ Branch 11 taken 200 times.
✓ Branch 12 taken 200 times.
✓ Branch 13 taken 200 times.
✓ Branch 14 taken 200 times.
✓ Branch 15 taken 200 times.
✓ Branch 16 taken 200 times.
✓ Branch 17 taken 200 times.
✓ Branch 18 taken 200 times.
✓ Branch 19 taken 200 times.
|
173318605 | { return !std::signbit(crossProduct<CoordScalar>(scvBasis[0], scvBasis[1])); } |
| 110 | |||
| 111 | /*! | ||
| 112 | * \brief Returns a vector containing the positions on a given intersection | ||
| 113 | * that are relevant for scvf corner computation. | ||
| 114 | * Ordering -> 1: facet center, 2: the two facet corners | ||
| 115 | * | ||
| 116 | * \param eg Geometry of the element the facet is embedded in | ||
| 117 | * \param refElement Reference element of the element the facet is embedded in | ||
| 118 | * \param indexInElement The local index of the facet in the element | ||
| 119 | * \param numCorners The number of corners on the facet | ||
| 120 | */ | ||
| 121 | template<class ElementGeometry, class ReferenceElement> | ||
| 122 | 55148863 | static ScvfPositionsOnIntersection computeScvfCornersOnIntersection(const ElementGeometry& eg, | |
| 123 | const ReferenceElement& refElement, | ||
| 124 | unsigned int indexInElement, | ||
| 125 | unsigned int numCorners) | ||
| 126 | { | ||
| 127 | 55148863 | ScvfPositionsOnIntersection p; | |
| 128 | |||
| 129 | // compute facet center and corners | ||
| 130 | 110297726 | p[0] = 0.0; | |
| 131 |
2/2✓ Branch 0 taken 110297726 times.
✓ Branch 1 taken 55148863 times.
|
165446589 | for (unsigned int c = 0; c < numCorners; ++c) |
| 132 | { | ||
| 133 | // codim = 1, dim = 2 | ||
| 134 | 119089698 | p[c+1] = eg.global(refElement.position(refElement.subEntity(indexInElement, 1, c, 2), 2)); | |
| 135 | 330893178 | p[0] += p[c+1]; | |
| 136 | } | ||
| 137 | 110297726 | p[0] /= numCorners; | |
| 138 | |||
| 139 | 55148863 | return p; | |
| 140 | } | ||
| 141 | |||
| 142 | /*! | ||
| 143 | * \brief Returns the corners of the sub control volume face constructed | ||
| 144 | * in a corner (vertex) of an intersection. | ||
| 145 | * | ||
| 146 | * \param p Container with all scvf corners of the intersection | ||
| 147 | * \param numIntersectionCorners Number of corners of the intersection (required in 3d) | ||
| 148 | * \param cornerIdx Local vertex index on the intersection | ||
| 149 | */ | ||
| 150 | ✗ | static ScvfCornerVector getScvfCorners(const ScvfPositionsOnIntersection& p, | |
| 151 | unsigned int numIntersectionCorners, | ||
| 152 | unsigned int cornerIdx) | ||
| 153 | { | ||
| 154 | // make sure the given input is admissible | ||
| 155 | ✗ | assert(cornerIdx < 2 && "provided index exceeds the number of corners of facets in 2d"); | |
| 156 | |||
| 157 | // create & return the scvf corner vector | ||
| 158 | ✗ | return cornerIdx == 0 ? ScvfCornerVector({p[0], p[1]}) | |
| 159 | ✗ | : ScvfCornerVector({p[0], p[2]}); | |
| 160 | } | ||
| 161 | |||
| 162 | /*! | ||
| 163 | * \brief Calculates the area of an scvf. | ||
| 164 | * \param scvfCorners Container with the corners of the scvf | ||
| 165 | */ | ||
| 166 | 51738064 | static CoordScalar computeScvfArea(const ScvfCornerVector& scvfCorners) | |
| 167 | 203923908 | { return (scvfCorners[1]-scvfCorners[0]).two_norm(); } | |
| 168 | |||
| 169 | /*! | ||
| 170 | * \brief Calculates the number of scvfs in a given element geometry type. | ||
| 171 | * \param gt The element geometry type | ||
| 172 | */ | ||
| 173 | 79782 | static std::size_t getNumLocalScvfs(const Dune::GeometryType& gt) | |
| 174 | { | ||
| 175 |
1/2✓ Branch 0 taken 79782 times.
✗ Branch 1 not taken.
|
79782 | if (gt == Dune::GeometryTypes::triangle) |
| 176 | return 6; | ||
| 177 |
1/2✓ Branch 0 taken 69915 times.
✗ Branch 1 not taken.
|
69915 | else if (gt == Dune::GeometryTypes::quadrilateral) |
| 178 | return 8; | ||
| 179 | else | ||
| 180 | ✗ | DUNE_THROW(Dune::NotImplemented, "Mpfa for 2d geometry type " << gt); | |
| 181 | } | ||
| 182 | }; | ||
| 183 | |||
| 184 | /*! | ||
| 185 | * \ingroup CCMpfaDiscretization | ||
| 186 | * \brief Dimension-specific mpfa helper class for dim == 2 & dimWorld == 3. | ||
| 187 | * Reuses some functionality of the specialization for dim = dimWorld = 2 | ||
| 188 | */ | ||
| 189 | template<class GridGeometry> | ||
| 190 | class MpfaDimensionHelper<GridGeometry, /*dim*/2, /*dimWorld*/3> | ||
| 191 | : public MpfaDimensionHelper<GridGeometry, 2, 2> | ||
| 192 | { | ||
| 193 | using GridView = typename GridGeometry::GridView; | ||
| 194 | using CoordScalar = typename GridView::ctype; | ||
| 195 | public: | ||
| 196 | |||
| 197 | /*! | ||
| 198 | * \brief Calculates the inner normal vectors to a given scv basis. | ||
| 199 | * \param scvBasis The basis of an scv | ||
| 200 | */ | ||
| 201 | template< class ScvBasis > | ||
| 202 | 2707270 | static ScvBasis calculateInnerNormals(const ScvBasis& scvBasis) | |
| 203 | { | ||
| 204 | // compute vector normal to the basis plane | ||
| 205 | 5414540 | const auto normal = [&scvBasis] () | |
| 206 | { | ||
| 207 | 21658160 | auto n = crossProduct<CoordScalar>(scvBasis[0], scvBasis[1]); | |
| 208 | 10829080 | n /= n.two_norm(); | |
| 209 | 5414540 | return n; | |
| 210 | 13536350 | } (); | |
| 211 | |||
| 212 | // compute inner normals using the normal vector | ||
| 213 | 2707270 | ScvBasis innerNormals; | |
| 214 | 5414540 | innerNormals[0] = crossProduct<CoordScalar>(scvBasis[1], normal); | |
| 215 | 5414540 | innerNormals[1] = crossProduct<CoordScalar>(normal, scvBasis[0]); | |
| 216 | |||
| 217 | 2707270 | return innerNormals; | |
| 218 | } | ||
| 219 | |||
| 220 | /*! | ||
| 221 | * \brief Calculates the determinant of an scv basis. | ||
| 222 | * For dim = 2 < dimWorld = 3 this is actually not the determinant of the | ||
| 223 | * basis but it is simply the area of the parallelofram spanned by the | ||
| 224 | * basis vectors. | ||
| 225 | * | ||
| 226 | * \param scvBasis The basis of an scv | ||
| 227 | */ | ||
| 228 | template< class ScvBasis > | ||
| 229 | 2707270 | static CoordScalar calculateDetX(const ScvBasis& scvBasis) | |
| 230 | { | ||
| 231 | using std::abs; | ||
| 232 | 8121810 | return abs(crossProduct<CoordScalar>(scvBasis[0], scvBasis[1]).two_norm()); | |
| 233 | } | ||
| 234 | |||
| 235 | /*! | ||
| 236 | * \brief Checks whether or not a given scv basis forms a right hand system. | ||
| 237 | * Note that for dim = 2 < dimWorld = 3 the bases forming a right hand system | ||
| 238 | * are not unique. | ||
| 239 | * | ||
| 240 | * \param scvBasis The basis of an scv | ||
| 241 | */ | ||
| 242 | template< class ScvBasis > | ||
| 243 | static constexpr bool isRightHandSystem(const ScvBasis& scvBasis) { return true; } | ||
| 244 | }; | ||
| 245 | /*! | ||
| 246 | * \ingroup CCMpfaDiscretization | ||
| 247 | * \brief Dimension-specific mpfa helper class for dim == 3 & dimWorld == 3. | ||
| 248 | * | ||
| 249 | */ | ||
| 250 | template<class GridGeometry> | ||
| 251 | class MpfaDimensionHelper<GridGeometry, /*dim*/3, /*dimWorld*/3> | ||
| 252 | { | ||
| 253 | using SubControlVolumeFace = typename GridGeometry::SubControlVolumeFace; | ||
| 254 | using ScvfCornerVector = typename SubControlVolumeFace::Traits::CornerStorage; | ||
| 255 | |||
| 256 | // Be picky about the dimensions | ||
| 257 | using GridView = typename GridGeometry::GridView; | ||
| 258 | using Element = typename GridView::template Codim<0>::Entity; | ||
| 259 | using CoordScalar = typename GridView::ctype; | ||
| 260 | using GlobalPosition = typename Element::Geometry::GlobalCoordinate; | ||
| 261 | |||
| 262 | // container to store the positions of intersections required for | ||
| 263 | // scvf corner computation. Maximum number of points needed is 9 | ||
| 264 | // for the supported geometry types (quadrilateral facet) | ||
| 265 | using ScvfPositionsOnIntersection = std::array<GlobalPosition, 9>; | ||
| 266 | |||
| 267 | public: | ||
| 268 | |||
| 269 | /*! | ||
| 270 | * \brief Calculates the inner normal vectors to a given scv basis. | ||
| 271 | * | ||
| 272 | * \param scvBasis The basis of an scv | ||
| 273 | */ | ||
| 274 | template< class ScvBasis > | ||
| 275 | static ScvBasis calculateInnerNormals(const ScvBasis& scvBasis) | ||
| 276 | { | ||
| 277 | ScvBasis innerNormals; | ||
| 278 | |||
| 279 | innerNormals[0] = crossProduct<CoordScalar>(scvBasis[1], scvBasis[2]); | ||
| 280 | innerNormals[1] = crossProduct<CoordScalar>(scvBasis[2], scvBasis[0]); | ||
| 281 | innerNormals[2] = crossProduct<CoordScalar>(scvBasis[0], scvBasis[1]); | ||
| 282 | |||
| 283 | if (!isRightHandSystem(scvBasis)) | ||
| 284 | std::for_each(innerNormals.begin(), innerNormals.end(), [] (auto& n) { n *= -1.0; }); | ||
| 285 | |||
| 286 | return innerNormals; | ||
| 287 | } | ||
| 288 | |||
| 289 | /*! | ||
| 290 | * \brief Calculates the determinant of an scv basis. | ||
| 291 | * This is equal to the cross product for dim = dimWorld = 2 | ||
| 292 | * | ||
| 293 | * \param scvBasis The basis of an scv | ||
| 294 | */ | ||
| 295 | template< class ScvBasis > | ||
| 296 | static CoordScalar calculateDetX(const ScvBasis& scvBasis) | ||
| 297 | { | ||
| 298 | using std::abs; | ||
| 299 | return abs(tripleProduct<CoordScalar>(scvBasis[0], scvBasis[1], scvBasis[2])); | ||
| 300 | } | ||
| 301 | |||
| 302 | /*! | ||
| 303 | * \brief Returns the global number of scvfs in the grid. Assumes the grid to be made up of only | ||
| 304 | * basic geometry types. Overload this function if you want to use different geometry types. | ||
| 305 | * | ||
| 306 | * \param gridView The grid view to be checked | ||
| 307 | */ | ||
| 308 | 2 | static std::size_t getGlobalNumScvf(const GridView& gridView) | |
| 309 | { | ||
| 310 |
1/2✗ Branch 5 not taken.
✓ Branch 6 taken 2 times.
|
2 | assert(gridView.size(Dune::GeometryTypes::tetrahedron) |
| 311 | + gridView.size(Dune::GeometryTypes::pyramid) | ||
| 312 | + gridView.size(Dune::GeometryTypes::prism) | ||
| 313 | + gridView.size(Dune::GeometryTypes::hexahedron) == gridView.size(0)); | ||
| 314 | |||
| 315 | return gridView.size(Dune::GeometryTypes::tetrahedron) | ||
| 316 | 2 | * getNumLocalScvfs(Dune::GeometryTypes::tetrahedron) | |
| 317 | 2 | + gridView.size(Dune::GeometryTypes::pyramid) | |
| 318 | 2 | * getNumLocalScvfs(Dune::GeometryTypes::pyramid) | |
| 319 | + gridView.size(Dune::GeometryTypes::prism) | ||
| 320 | 2 | * getNumLocalScvfs(Dune::GeometryTypes::prism) | |
| 321 | 2 | + gridView.size(Dune::GeometryTypes::hexahedron) | |
| 322 | 2 | * getNumLocalScvfs(Dune::GeometryTypes::hexahedron); | |
| 323 | } | ||
| 324 | |||
| 325 | /*! | ||
| 326 | * \brief Checks whether or not a given scv basis forms a right hand system. | ||
| 327 | * \param scvBasis The basis of an scv | ||
| 328 | */ | ||
| 329 | template< class ScvBasis > | ||
| 330 | static bool isRightHandSystem(const ScvBasis& scvBasis) | ||
| 331 | { return !std::signbit(tripleProduct<CoordScalar>(scvBasis[0], scvBasis[1], scvBasis[2])); } | ||
| 332 | |||
| 333 | /*! | ||
| 334 | * \brief Returns a vector containing the positions on a given intersection | ||
| 335 | * that are relevant for scvf corner computation. | ||
| 336 | * Ordering -> 1: facet center, 2: facet corners, 3: edge centers | ||
| 337 | * | ||
| 338 | * \param eg Geometry of the element the facet is embedded in | ||
| 339 | * \param refElement Reference element of the element the facet is embedded in | ||
| 340 | * \param indexInElement The local index of the facet in the element | ||
| 341 | * \param numCorners The number of corners on the facet | ||
| 342 | */ | ||
| 343 | template<class ElementGeometry, class ReferenceElement> | ||
| 344 | 3928 | static ScvfPositionsOnIntersection computeScvfCornersOnIntersection(const ElementGeometry& eg, | |
| 345 | const ReferenceElement& refElement, | ||
| 346 | unsigned int indexInElement, | ||
| 347 | unsigned int numCorners) | ||
| 348 | { | ||
| 349 | // The size of ScvfPositionsOnIntersection doesn't allow for faces with more than four corners! | ||
| 350 | 3928 | ScvfPositionsOnIntersection p; | |
| 351 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 3928 times.
|
3928 | if (numCorners > 4) |
| 352 | ✗ | DUNE_THROW(Dune::InvalidStateException, "Mpfa implementation cannot handle faces with more than 4 corners"); | |
| 353 | |||
| 354 | // compute facet center and corners | ||
| 355 | p[0] = 0.0; | ||
| 356 |
2/2✓ Branch 0 taken 11784 times.
✓ Branch 1 taken 3928 times.
|
15712 | for (unsigned int c = 0; c < numCorners; ++c) |
| 357 | { | ||
| 358 | // codim = 1, dim = 3 | ||
| 359 | 11784 | p[c+1] = eg.global(refElement.position(refElement.subEntity(indexInElement, 1, c, 3), 3)); | |
| 360 | 35352 | p[0] += p[c+1]; | |
| 361 | } | ||
| 362 | 7856 | p[0] /= numCorners; | |
| 363 | |||
| 364 | // proceed according to number of corners | ||
| 365 |
1/3✓ Branch 0 taken 3928 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
|
3928 | switch (numCorners) |
| 366 | { | ||
| 367 | 3928 | case 3: // triangle | |
| 368 | { | ||
| 369 | // add edge midpoints, there are 3 for triangles | ||
| 370 | 15712 | p[numCorners+1] = p[2] + p[1]; | |
| 371 | 7856 | p[numCorners+1] /= 2; | |
| 372 | 15712 | p[numCorners+2] = p[3] + p[1]; | |
| 373 | 7856 | p[numCorners+2] /= 2; | |
| 374 | 15712 | p[numCorners+3] = p[3] + p[2]; | |
| 375 | 7856 | p[numCorners+3] /= 2; | |
| 376 | return p; | ||
| 377 | } | ||
| 378 | ✗ | case 4: // quadrilateral | |
| 379 | { | ||
| 380 | // add edge midpoints, there are 4 for quadrilaterals | ||
| 381 | ✗ | p[numCorners+1] = p[3] + p[1]; | |
| 382 | ✗ | p[numCorners+1] /= 2; | |
| 383 | ✗ | p[numCorners+2] = p[4] + p[2]; | |
| 384 | ✗ | p[numCorners+2] /= 2; | |
| 385 | ✗ | p[numCorners+3] = p[2] + p[1]; | |
| 386 | ✗ | p[numCorners+3] /= 2; | |
| 387 | ✗ | p[numCorners+4] = p[4] + p[3]; | |
| 388 | ✗ | p[numCorners+4] /= 2; | |
| 389 | return p; | ||
| 390 | } | ||
| 391 | ✗ | default: | |
| 392 | ✗ | DUNE_THROW(Dune::NotImplemented, | |
| 393 | "Mpfa scvf corners for dim = 3, dimWorld = 3, corners = " << numCorners); | ||
| 394 | } | ||
| 395 | } | ||
| 396 | |||
| 397 | /*! | ||
| 398 | * \brief Returns the corners of the sub control volume face constructed | ||
| 399 | * in a corner (vertex) of an intersection. | ||
| 400 | * | ||
| 401 | * \param p Container with all scvf corners of the intersection | ||
| 402 | * \param numIntersectionCorners Number of corners of the intersection | ||
| 403 | * \param cornerIdx Local vertex index on the intersection | ||
| 404 | */ | ||
| 405 | 11784 | static ScvfCornerVector getScvfCorners(const ScvfPositionsOnIntersection& p, | |
| 406 | unsigned int numIntersectionCorners, | ||
| 407 | unsigned int cornerIdx) | ||
| 408 | { | ||
| 409 | // proceed according to number of corners | ||
| 410 | // we assume the ordering according to the above method computeScvfCornersOnIntersection() | ||
| 411 |
1/3✓ Branch 0 taken 11784 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
|
11784 | switch (numIntersectionCorners) |
| 412 | { | ||
| 413 | 11784 | case 3: // triangle | |
| 414 | { | ||
| 415 | // Only build the maps the first time we encounter a triangle | ||
| 416 | static const std::uint8_t vo = 1; // vertex offset in point vector p | ||
| 417 | static const std::uint8_t eo = 4; // edge offset in point vector p | ||
| 418 | static const std::uint8_t map[3][4] = | ||
| 419 | { | ||
| 420 | {0, eo+1, eo+0, vo+0}, | ||
| 421 | {0, eo+0, eo+2, vo+1}, | ||
| 422 | {0, eo+2, eo+1, vo+2} | ||
| 423 | }; | ||
| 424 | |||
| 425 | 11784 | return ScvfCornerVector( {p[map[cornerIdx][0]], | |
| 426 | 23568 | p[map[cornerIdx][1]], | |
| 427 | 23568 | p[map[cornerIdx][2]], | |
| 428 | 47136 | p[map[cornerIdx][3]]} ); | |
| 429 | } | ||
| 430 | ✗ | case 4: // quadrilateral | |
| 431 | { | ||
| 432 | // Only build the maps the first time we encounter a quadrilateral | ||
| 433 | static const std::uint8_t vo = 1; // vertex offset in point vector p | ||
| 434 | static const std::uint8_t eo = 5; // face offset in point vector p | ||
| 435 | static const std::uint8_t map[4][4] = | ||
| 436 | { | ||
| 437 | {0, eo+0, eo+2, vo+0}, | ||
| 438 | {0, eo+2, eo+1, vo+1}, | ||
| 439 | {0, eo+3, eo+0, vo+2}, | ||
| 440 | {0, eo+1, eo+3, vo+3} | ||
| 441 | }; | ||
| 442 | |||
| 443 | ✗ | return ScvfCornerVector( {p[map[cornerIdx][0]], | |
| 444 | ✗ | p[map[cornerIdx][1]], | |
| 445 | ✗ | p[map[cornerIdx][2]], | |
| 446 | ✗ | p[map[cornerIdx][3]]} ); | |
| 447 | } | ||
| 448 | ✗ | default: | |
| 449 | ✗ | DUNE_THROW(Dune::NotImplemented, | |
| 450 | "Mpfa scvf corners for dim = 3, dimWorld = 3, corners = " << numIntersectionCorners); | ||
| 451 | } | ||
| 452 | } | ||
| 453 | |||
| 454 | /*! | ||
| 455 | * \brief Calculates the area of an scvf. | ||
| 456 | * \param scvfCorners Container with the corners of the scvf | ||
| 457 | */ | ||
| 458 | 11784 | static CoordScalar computeScvfArea(const ScvfCornerVector& scvfCorners) | |
| 459 | { | ||
| 460 | // after Wolfram alpha quadrilateral area | ||
| 461 | 82488 | return 0.5*Dumux::crossProduct(scvfCorners[3]-scvfCorners[0], scvfCorners[2]-scvfCorners[1]).two_norm(); | |
| 462 | } | ||
| 463 | |||
| 464 | /*! | ||
| 465 | * \brief Calculates the number of scvfs in a given element geometry type. | ||
| 466 | * | ||
| 467 | * \param gt The element geometry type | ||
| 468 | */ | ||
| 469 | 990 | static std::size_t getNumLocalScvfs(const Dune::GeometryType& gt) | |
| 470 | { | ||
| 471 |
1/2✓ Branch 0 taken 990 times.
✗ Branch 1 not taken.
|
990 | if (gt == Dune::GeometryTypes::tetrahedron) |
| 472 | return 12; | ||
| 473 |
1/2✓ Branch 0 taken 6 times.
✗ Branch 1 not taken.
|
6 | else if (gt == Dune::GeometryTypes::pyramid) |
| 474 | return 16; | ||
| 475 |
1/2✓ Branch 0 taken 4 times.
✗ Branch 1 not taken.
|
4 | else if (gt == Dune::GeometryTypes::prism) |
| 476 | return 18; | ||
| 477 |
1/2✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
|
2 | else if (gt == Dune::GeometryTypes::hexahedron) |
| 478 | return 24; | ||
| 479 | else | ||
| 480 | ✗ | DUNE_THROW(Dune::NotImplemented, "Mpfa for 3d geometry type " << gt); | |
| 481 | } | ||
| 482 | }; | ||
| 483 | |||
| 484 | /*! | ||
| 485 | * \ingroup CCMpfaDiscretization | ||
| 486 | * \brief Helper class to get the required information on an interaction volume. | ||
| 487 | * | ||
| 488 | * \tparam GridGeometry The finite volume grid geometry | ||
| 489 | */ | ||
| 490 | template<class GridGeometry> | ||
| 491 | class CCMpfaHelper : public MpfaDimensionHelper<GridGeometry, | ||
| 492 | GridGeometry::GridView::dimension, | ||
| 493 | GridGeometry::GridView::dimensionworld> | ||
| 494 | { | ||
| 495 | using PrimaryInteractionVolume = typename GridGeometry::GridIVIndexSets::PrimaryInteractionVolume; | ||
| 496 | using SecondaryInteractionVolume = typename GridGeometry::GridIVIndexSets::SecondaryInteractionVolume; | ||
| 497 | |||
| 498 | using VertexMapper = typename GridGeometry::VertexMapper; | ||
| 499 | using FVElementGeometry = typename GridGeometry::LocalView; | ||
| 500 | using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace; | ||
| 501 | using ScvfCornerVector = typename SubControlVolumeFace::Traits::CornerStorage; | ||
| 502 | |||
| 503 | using GridView = typename GridGeometry::GridView; | ||
| 504 | static constexpr int dim = GridView::dimension; | ||
| 505 | |||
| 506 | using Element = typename GridView::template Codim<0>::Entity; | ||
| 507 | using GlobalPosition = typename Element::Geometry::GlobalCoordinate; | ||
| 508 | using CoordScalar = typename GridView::ctype; | ||
| 509 | |||
| 510 | public: | ||
| 511 | /*! | ||
| 512 | * \brief Calculates the integration point on an scvf. | ||
| 513 | * | ||
| 514 | * \param scvfCorners Container with the corners of the scvf | ||
| 515 | * \param q Parameterization of the integration point on the scvf | ||
| 516 | */ | ||
| 517 | template< class Scalar > | ||
| 518 | 51749848 | static GlobalPosition getScvfIntegrationPoint(const ScvfCornerVector& scvfCorners, Scalar q) | |
| 519 | { | ||
| 520 | // ordering -> first corner: facet center, last corner: vertex | ||
| 521 |
1/2✓ Branch 0 taken 51749848 times.
✗ Branch 1 not taken.
|
51749848 | if (q == 0.0) |
| 522 | 103499696 | return scvfCorners[0]; | |
| 523 | |||
| 524 | ✗ | auto ip = scvfCorners.back() - scvfCorners.front(); | |
| 525 | ✗ | ip *= q; | |
| 526 | ✗ | ip += scvfCorners[0]; | |
| 527 | ✗ | return ip; | |
| 528 | } | ||
| 529 | |||
| 530 | /*! | ||
| 531 | * \brief Returns a vector which maps true to each vertex on processor boundaries and false otherwise | ||
| 532 | * \todo TODO: The name of this function is not so good, as these are not ghost vertices according | ||
| 533 | * to the Dune definition of ghost entities. Moreover, it should be tried to make MPFA work | ||
| 534 | * also with ghost entities. | ||
| 535 | */ | ||
| 536 | 50 | static std::vector<bool> findGhostVertices(const GridView& gridView, const VertexMapper& vertexMapper) | |
| 537 | { | ||
| 538 |
4/8✓ Branch 1 taken 32 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 50 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 32 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 32 times.
✗ Branch 11 not taken.
|
100 | std::vector<bool> ghostVertices(gridView.size(dim), false); |
| 539 | |||
| 540 | // if not run in parallel, skip the rest | ||
| 541 |
5/6✓ Branch 1 taken 50 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 2 times.
✓ Branch 4 taken 48 times.
✓ Branch 5 taken 2 times.
✓ Branch 6 taken 48 times.
|
50 | if (Dune::MPIHelper::getCommunication().size() == 1) |
| 542 | return ghostVertices; | ||
| 543 | |||
| 544 | // mpfa methods cannot yet handle ghost cells | ||
| 545 |
1/5✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
|
2 | if (gridView.ghostSize(0) > 0) |
| 546 | ✗ | DUNE_THROW(Dune::InvalidStateException, "Mpfa methods in parallel do not work with ghost cells. Use overlap cells instead."); | |
| 547 | |||
| 548 | // mpfa methods have to have overlapping cells | ||
| 549 |
2/5✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 2 times.
|
2 | if (gridView.overlapSize(0) == 0) |
| 550 | ✗ | DUNE_THROW(Dune::InvalidStateException, "Grid no overlapping cells. This is required by mpfa methods in parallel."); | |
| 551 | |||
| 552 |
2/4✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 418 times.
✗ Branch 4 not taken.
|
834 | for (const auto& element : elements(gridView)) |
| 553 | { | ||
| 554 |
2/2✓ Branch 1 taken 416 times.
✓ Branch 2 taken 1664 times.
|
2496 | for (const auto& is : intersections(gridView, element)) |
| 555 | { | ||
| 556 |
3/4✓ Branch 0 taken 1606 times.
✓ Branch 1 taken 58 times.
✓ Branch 2 taken 116 times.
✗ Branch 3 not taken.
|
1722 | if (!is.neighbor() && !is.boundary()) |
| 557 | { | ||
| 558 |
1/2✓ Branch 1 taken 32 times.
✗ Branch 2 not taken.
|
32 | const auto refElement = referenceElement(element); |
| 559 | |||
| 560 |
4/4✓ Branch 1 taken 64 times.
✓ Branch 2 taken 32 times.
✓ Branch 3 taken 64 times.
✓ Branch 4 taken 32 times.
|
96 | for (int isVertex = 0; isVertex < is.geometry().corners(); ++isVertex) |
| 561 | { | ||
| 562 | 128 | const auto vIdxLocal = refElement.subEntity(is.indexInInside(), 1, isVertex, dim); | |
| 563 |
1/2✓ Branch 1 taken 64 times.
✗ Branch 2 not taken.
|
64 | const auto vIdxGlobal = vertexMapper.subIndex(element, vIdxLocal, dim); |
| 564 | 192 | ghostVertices[vIdxGlobal] = true; | |
| 565 | } | ||
| 566 | } | ||
| 567 | } | ||
| 568 | } | ||
| 569 | |||
| 570 | 2 | return ghostVertices; | |
| 571 | } | ||
| 572 | |||
| 573 | //! Returns whether or not secondary interaction volumes have to be considered in the model. | ||
| 574 | //! This is always the case when the specified types for the interaction volumes differ. | ||
| 575 | static constexpr bool considerSecondaryIVs() | ||
| 576 | { return !std::is_same<PrimaryInteractionVolume, SecondaryInteractionVolume>::value; } | ||
| 577 | |||
| 578 | //! returns whether or not a value exists in a vector | ||
| 579 | template<typename V, typename T> | ||
| 580 | static bool vectorContainsValue(const V& vector, const T value) | ||
| 581 |
19/24✓ Branch 0 taken 257944 times.
✓ Branch 1 taken 506292 times.
✓ Branch 2 taken 257944 times.
✓ Branch 3 taken 506292 times.
✓ Branch 4 taken 312520 times.
✓ Branch 5 taken 615376 times.
✓ Branch 6 taken 312520 times.
✓ Branch 7 taken 615376 times.
✓ Branch 8 taken 257944 times.
✓ Branch 9 taken 506292 times.
✓ Branch 12 taken 3406072 times.
✓ Branch 13 taken 2061704 times.
✓ Branch 14 taken 58528302 times.
✓ Branch 15 taken 40137020 times.
✓ Branch 16 taken 55122230 times.
✓ Branch 17 taken 38075316 times.
✗ Branch 18 not taken.
✗ Branch 19 not taken.
✗ Branch 20 not taken.
✓ Branch 21 taken 3406072 times.
✗ Branch 22 not taken.
✓ Branch 23 taken 58528302 times.
✗ Branch 24 not taken.
✓ Branch 25 taken 55122230 times.
|
399137108 | { return std::find(vector.begin(), vector.end(), value) != vector.end(); } |
| 582 | }; | ||
| 583 | |||
| 584 | } // end namespace Dumux | ||
| 585 | |||
| 586 | #endif | ||
| 587 |