| 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-FileCopyrightText: 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 StaggeredDiscretization | ||
| 10 | * \copydoc Dumux::FreeFlowStaggeredGeometryHelper | ||
| 11 | */ | ||
| 12 | #ifndef DUMUX_DISCRETIZATION_STAGGERED_GEOMETRY_HELPER_HH | ||
| 13 | #define DUMUX_DISCRETIZATION_STAGGERED_GEOMETRY_HELPER_HH | ||
| 14 | |||
| 15 | #include <dune/geometry/multilineargeometry.hh> | ||
| 16 | |||
| 17 | #include <dumux/common/math.hh> | ||
| 18 | #include <dumux/common/indextraits.hh> | ||
| 19 | #include <type_traits> | ||
| 20 | #include <algorithm> | ||
| 21 | #include <array> | ||
| 22 | #include <bitset> | ||
| 23 | |||
| 24 | namespace Dumux { | ||
| 25 | namespace Detail { | ||
| 26 | |||
| 27 | /*! | ||
| 28 | * \ingroup StaggeredDiscretization | ||
| 29 | * \brief Parallel Data stored per sub face | ||
| 30 | * | ||
| 31 | * ------------ | ||
| 32 | * | | | ||
| 33 | * | | | ||
| 34 | * | | | ||
| 35 | * ----------------------- | ||
| 36 | * | yyyyyyyy s | | ||
| 37 | * | yyyyyyyy s | | ||
| 38 | * | yyyyyyyy s | | ||
| 39 | * ----------------------- | ||
| 40 | * In this corner geometry, hasParallelNeighbor will return true for subcontrolvolumeface s belonging to the | ||
| 41 | * element filled by 'y's, but hasParallelNeighbor will return false for the subcontrolvolumeface that has the | ||
| 42 | * same dofIndex. We name this situation hasHalfParallelNeighbor. | ||
| 43 | * | ||
| 44 | * ------------ | ||
| 45 | * | yyyyyyyy s | ||
| 46 | * | yyyyyyyy s | ||
| 47 | * | yyyyyyyy s | ||
| 48 | * ----------------------- | ||
| 49 | * | | | | ||
| 50 | * | | | | ||
| 51 | * | | | | ||
| 52 | * ----------------------- | ||
| 53 | * In this corner geometry, hasParallelNeighbor will return true for subcontrolvolumeface s belonging to the | ||
| 54 | * element filled by 'y's. However, as there also might be a boundary velocity value known at the corner, which | ||
| 55 | * can be used instead of the standard parallel velocity in some cases, we want to identify this situation. We | ||
| 56 | * name it cornerParallelNeighbor. | ||
| 57 | */ | ||
| 58 | template<class GridView, int upwindSchemeOrder> | ||
| 59 | 187254 | struct PairData | |
| 60 | { | ||
| 61 | using Scalar = typename GridView::ctype; | ||
| 62 | using GlobalPosition = typename GridView::template Codim<0>::Entity::Geometry::GlobalCoordinate; | ||
| 63 | using GridIndexType = typename IndexTraits<GridView>::GridIndex; | ||
| 64 | using SmallLocalIndexType = typename IndexTraits<GridView>::SmallLocalIndex; | ||
| 65 | |||
| 66 | std::bitset<upwindSchemeOrder> hasParallelNeighbor; | ||
| 67 | bool hasHalfParallelNeighbor = false; | ||
| 68 | bool hasCornerParallelNeighbor = false; | ||
| 69 | std::array<GridIndexType, upwindSchemeOrder> parallelDofs; | ||
| 70 | std::array<Scalar, upwindSchemeOrder> parallelCellWidths; | ||
| 71 | bool hasOuterLateral = false; | ||
| 72 | std::pair<GridIndexType, GridIndexType> lateralPair; | ||
| 73 | SmallLocalIndexType localLateralFaceIdx; | ||
| 74 | Scalar lateralDistance; | ||
| 75 | GlobalPosition lateralStaggeredFaceCenter; | ||
| 76 | }; | ||
| 77 | |||
| 78 | /*! | ||
| 79 | * \ingroup StaggeredDiscretization | ||
| 80 | * \brief In Axis Data stored per sub face | ||
| 81 | */ | ||
| 82 | template<class GridView, int upwindSchemeOrder> | ||
| 83 | 4625 | struct AxisData | |
| 84 | { | ||
| 85 | using Scalar = typename GridView::ctype; | ||
| 86 | using GridIndexType = typename IndexTraits<GridView>::GridIndex; | ||
| 87 | |||
| 88 | GridIndexType selfDof; | ||
| 89 | GridIndexType oppositeDof; | ||
| 90 | std::bitset<upwindSchemeOrder-1> hasForwardNeighbor; | ||
| 91 | std::bitset<upwindSchemeOrder-1> hasBackwardNeighbor; | ||
| 92 | std::array<GridIndexType, upwindSchemeOrder-1> inAxisForwardDofs; | ||
| 93 | std::array<GridIndexType, upwindSchemeOrder-1> inAxisBackwardDofs; | ||
| 94 | Scalar selfToOppositeDistance; | ||
| 95 | std::array<Scalar, upwindSchemeOrder-1> inAxisForwardDistances; | ||
| 96 | std::array<Scalar, upwindSchemeOrder-1> inAxisBackwardDistances; | ||
| 97 | }; | ||
| 98 | |||
| 99 | /*! | ||
| 100 | * \ingroup StaggeredDiscretization | ||
| 101 | * \brief In Axis Data stored per sub face for first-order scheme | ||
| 102 | */ | ||
| 103 | template<class GridView> | ||
| 104 | struct AxisData<GridView, 1> | ||
| 105 | { | ||
| 106 | using Scalar = typename GridView::ctype; | ||
| 107 | using GridIndexType = typename IndexTraits<GridView>::GridIndex; | ||
| 108 | |||
| 109 | GridIndexType selfDof; | ||
| 110 | GridIndexType oppositeDof; | ||
| 111 | Scalar selfToOppositeDistance; | ||
| 112 | }; | ||
| 113 | |||
| 114 | } // namespace Detail | ||
| 115 | |||
| 116 | /*! | ||
| 117 | * \ingroup StaggeredDiscretization | ||
| 118 | * \brief Returns the direction index of the facet (0 = x, 1 = y, 2 = z) | ||
| 119 | */ | ||
| 120 | template<class Vector> | ||
| 121 | 1088150 | inline static unsigned int directionIndex(Vector&& vector) | |
| 122 | { | ||
| 123 | 1088150 | const auto eps = 1e-8; | |
| 124 |
6/16✗ Branch 0 not taken.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✓ Branch 8 taken 1050 times.
✓ Branch 9 taken 2100 times.
✓ Branch 10 taken 543028 times.
✓ Branch 11 taken 544072 times.
✓ Branch 12 taken 544072 times.
✗ Branch 13 not taken.
✗ Branch 15 not taken.
✓ Branch 16 taken 1088150 times.
|
2722472 | return std::find_if(vector.begin(), vector.end(), [eps](const auto& x) { return std::abs(x) > eps; } ) - vector.begin(); |
| 125 | } | ||
| 126 | |||
| 127 | /*! | ||
| 128 | * \ingroup StaggeredDiscretization | ||
| 129 | * \brief Helper class constructing the dual grid finite volume geometries | ||
| 130 | * for the free flow staggered discretization method. | ||
| 131 | */ | ||
| 132 | template<class GridView, int upwindSchemeOrder> | ||
| 133 | class FreeFlowStaggeredGeometryHelper | ||
| 134 | { | ||
| 135 | using Scalar = typename GridView::ctype; | ||
| 136 | static constexpr auto dim = GridView::dimension; | ||
| 137 | |||
| 138 | using Element = typename GridView::template Codim<0>::Entity; | ||
| 139 | using Intersection = typename GridView::Intersection; | ||
| 140 | using GlobalPosition = typename Element::Geometry::GlobalCoordinate; | ||
| 141 | |||
| 142 | using GridIndexType = typename IndexTraits<GridView>::GridIndex; | ||
| 143 | using SmallLocalIndexType = typename IndexTraits<GridView>::SmallLocalIndex; | ||
| 144 | |||
| 145 | //TODO include assert that checks for quad geometry | ||
| 146 | static constexpr auto codimIntersection = 1; | ||
| 147 | static constexpr auto numFacets = dim * 2; | ||
| 148 | static constexpr auto numPairs = 2 * (dim - 1); | ||
| 149 | |||
| 150 | static constexpr bool useHigherOrder = upwindSchemeOrder > 1; | ||
| 151 | |||
| 152 | public: | ||
| 153 | using PairData = Detail::PairData<GridView, upwindSchemeOrder>; | ||
| 154 | using AxisData = Detail::AxisData<GridView, upwindSchemeOrder>; | ||
| 155 | |||
| 156 | 93502 | FreeFlowStaggeredGeometryHelper(const Element& element, const GridView& gridView) | |
| 157 | 93502 | : element_(element) | |
| 158 | 179879 | , gridView_(gridView) | |
| 159 | 93502 | { } | |
| 160 | |||
| 161 | //! update the local face | ||
| 162 | template<class IntersectionMapper> | ||
| 163 | 374258 | void updateLocalFace(const IntersectionMapper&, const Intersection& intersection) | |
| 164 | { | ||
| 165 |
1/2✓ Branch 1 taken 750 times.
✗ Branch 2 not taken.
|
374258 | intersection_ = intersection; |
| 166 |
1/2✓ Branch 1 taken 750 times.
✗ Branch 2 not taken.
|
374258 | fillAxisData_(); |
| 167 |
2/3✓ Branch 1 taken 311118 times.
✓ Branch 2 taken 19780 times.
✗ Branch 3 not taken.
|
374258 | fillPairData_(); |
| 168 | 358898 | } | |
| 169 | |||
| 170 | /*! | ||
| 171 | * \brief Returns the local index of the face (i.e. the intersection) | ||
| 172 | */ | ||
| 173 | 701856 | SmallLocalIndexType localFaceIndex() const | |
| 174 | { | ||
| 175 |
4/4✓ Branch 0 taken 375 times.
✓ Branch 1 taken 336280 times.
✓ Branch 2 taken 3138 times.
✓ Branch 3 taken 75 times.
|
701856 | return intersection_.indexInInside(); |
| 176 | } | ||
| 177 | |||
| 178 | /*! | ||
| 179 | * \brief Returns a copy of the axis data | ||
| 180 | */ | ||
| 181 | 362738 | AxisData axisData() const | |
| 182 | { | ||
| 183 |
2/2✓ Branch 0 taken 375 times.
✓ Branch 1 taken 375 times.
|
362738 | return axisData_; |
| 184 | } | ||
| 185 | |||
| 186 | /*! | ||
| 187 | * \brief Returns a copy of the pair data | ||
| 188 | */ | ||
| 189 | 362738 | std::array<PairData, numPairs> pairData() const | |
| 190 | { | ||
| 191 |
2/2✓ Branch 0 taken 375 times.
✓ Branch 1 taken 375 times.
|
362738 | return pairData_; |
| 192 | } | ||
| 193 | |||
| 194 | /*! | ||
| 195 | * \brief Returns the direction index of the primary facet (0 = x, 1 = y, 2 = z) | ||
| 196 | */ | ||
| 197 | 362738 | unsigned int directionIndex() const | |
| 198 | { | ||
| 199 |
2/2✓ Branch 0 taken 167369 times.
✓ Branch 1 taken 167369 times.
|
530107 | return Dumux::directionIndex(intersection_.centerUnitOuterNormal()); |
| 200 | } | ||
| 201 | |||
| 202 | /*! | ||
| 203 | * \brief Returns the direction index of the facet passed as an argument (0 = x, 1 = y, 2 = z) | ||
| 204 | */ | ||
| 205 |
2/2✓ Branch 0 taken 334146 times.
✓ Branch 1 taken 334066 times.
|
725412 | unsigned int directionIndex(const Intersection& intersection) const |
| 206 | { | ||
| 207 | 725412 | return Dumux::directionIndex(std::move(intersection.centerUnitOuterNormal())); | |
| 208 | } | ||
| 209 | |||
| 210 | private: | ||
| 211 | |||
| 212 | /*! | ||
| 213 | * \brief Fills all entries of the in axis data | ||
| 214 | * Calls a function to extend the axis data for higher order upwind methods if required. | ||
| 215 | */ | ||
| 216 | 374258 | void fillAxisData_() | |
| 217 | { | ||
| 218 | if constexpr (useHigherOrder) | ||
| 219 | 18750 | fillOuterAxisData_(); | |
| 220 | |||
| 221 |
2/2✓ Branch 0 taken 187129 times.
✓ Branch 1 taken 187129 times.
|
374258 | const auto inIdx = intersection_.indexInInside(); |
| 222 | 374258 | const auto oppIdx = localOppositeIdx_(inIdx); | |
| 223 | |||
| 224 | // Set the self Dof | ||
| 225 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 374258 times.
|
374258 | axisData_.selfDof = gridView_.indexSet().subIndex(intersection_.inside(), inIdx, codimIntersection); |
| 226 | |||
| 227 | // Set the opposite Dof | ||
| 228 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 28000 times.
|
374258 | axisData_.oppositeDof = gridView_.indexSet().subIndex(intersection_.inside(), oppIdx, codimIntersection); |
| 229 | |||
| 230 | // Set the Self to Opposite Distance | ||
| 231 | 374258 | const auto self = getFacet_(inIdx, element_); | |
| 232 | 374258 | const auto opposite = getFacet_(oppIdx, element_); | |
| 233 | 1497032 | axisData_.selfToOppositeDistance = (self.geometry().center() - opposite.geometry().center()).two_norm(); | |
| 234 | |||
| 235 | 374258 | } | |
| 236 | |||
| 237 | /*! | ||
| 238 | * \brief Fills the axis data with the outer dofs and distances. | ||
| 239 | * This function builds and extended stencil (see AxisData for upwindSchemeOrder > 1) | ||
| 240 | * and is therefore only called when higher order upwinding methods are prescribed. | ||
| 241 | */ | ||
| 242 | 18750 | void fillOuterAxisData_() | |
| 243 | { | ||
| 244 | // reset the axis data struct | ||
| 245 | 18750 | axisData_ = {}; | |
| 246 | |||
| 247 | // Set the forward dofs | ||
| 248 | 18750 | const auto thisFaceLocalIdx = intersection_.indexInInside(); | |
| 249 | 18750 | const auto& faceCenterPos = getFacet_(thisFaceLocalIdx, element_).geometry().center(); | |
| 250 |
2/2✓ Branch 0 taken 18675 times.
✓ Branch 1 taken 75 times.
|
18750 | if (intersection_.neighbor()) |
| 251 | 18200 | addForwardNeighborAxisData_(intersection_.outside(), 0, thisFaceLocalIdx, faceCenterPos); | |
| 252 | |||
| 253 | // Set the backward dofs | ||
| 254 | 18750 | const auto oppositeFaceLocalIdx = localOppositeIdx_(thisFaceLocalIdx); | |
| 255 | 18750 | const auto& oppositeFaceCenterPos = getFacet_(oppositeFaceLocalIdx, element_).geometry().center(); | |
| 256 |
8/10✓ Branch 2 taken 21000 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 20250 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 18000 times.
✓ Branch 8 taken 27600 times.
✓ Branch 9 taken 45600 times.
✓ Branch 10 taken 400 times.
✓ Branch 1 taken 150 times.
✓ Branch 4 taken 750 times.
|
121150 | for (const auto& intersection : intersections(gridView_, element_)) |
| 257 | { | ||
| 258 |
3/4✓ Branch 0 taken 18750 times.
✓ Branch 1 taken 29850 times.
✓ Branch 3 taken 28000 times.
✗ Branch 4 not taken.
|
49150 | if (intersection.indexInInside() == oppositeFaceLocalIdx && intersection.neighbor()) |
| 259 | { | ||
| 260 |
2/4✓ Branch 1 taken 17600 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 17600 times.
✗ Branch 5 not taken.
|
18200 | addBackwardNeighborAxisData_(intersection.outside(), 0, oppositeFaceLocalIdx, oppositeFaceCenterPos); |
| 261 | 18200 | break; | |
| 262 | } | ||
| 263 | } | ||
| 264 | 18750 | } | |
| 265 | |||
| 266 | /*! | ||
| 267 | * \brief Recursively add axis data in forward direction | ||
| 268 | */ | ||
| 269 | 18200 | void addForwardNeighborAxisData_(const Element& neighbor, | |
| 270 | const std::size_t distance, | ||
| 271 | const unsigned int localIntersectionIndex, | ||
| 272 | const GlobalPosition& faceCenterPos) | ||
| 273 | { | ||
| 274 | // fill the forward axis data for this element in the stencil | ||
| 275 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 18200 times.
|
18200 | axisData_.hasForwardNeighbor.set(distance, true); |
| 276 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 18200 times.
|
18200 | axisData_.inAxisForwardDofs[distance] = gridView_.indexSet().subIndex(neighbor, localIntersectionIndex, codimIntersection); |
| 277 | 36400 | const auto& forwardFacePos = getFacet_(localIntersectionIndex, neighbor).geometry().center(); | |
| 278 | 18200 | axisData_.inAxisForwardDistances[distance] = (forwardFacePos - faceCenterPos).two_norm(); | |
| 279 | |||
| 280 | // recursion end if we added all axis data for the given stencil size (numParallelFaces) | ||
| 281 | 18200 | const auto numForwardBackwardAxisDofs = axisData_.inAxisForwardDofs.size(); | |
| 282 | if (distance >= numForwardBackwardAxisDofs-1) | ||
| 283 | 600 | return; | |
| 284 | |||
| 285 | for (const auto& intersection : intersections(gridView_, neighbor)) | ||
| 286 | { | ||
| 287 | if (intersection.indexInInside() == localIntersectionIndex && intersection.neighbor()) | ||
| 288 | { | ||
| 289 | addForwardNeighborAxisData_(intersection.outside(), distance+1, localIntersectionIndex, forwardFacePos); | ||
| 290 | break; | ||
| 291 | } | ||
| 292 | } | ||
| 293 | } | ||
| 294 | |||
| 295 | /*! | ||
| 296 | * \brief Recursively add axis data in backward direction | ||
| 297 | */ | ||
| 298 | 18200 | void addBackwardNeighborAxisData_(const Element& neighbor, | |
| 299 | const std::size_t distance, | ||
| 300 | const unsigned int localIntersectionIndex, | ||
| 301 | const GlobalPosition& faceCenterPos) | ||
| 302 | { | ||
| 303 | // fill the forward axis data for this element in the stencil | ||
| 304 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 18200 times.
|
18200 | axisData_.hasBackwardNeighbor.set(distance, true); |
| 305 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 18200 times.
|
18200 | axisData_.inAxisBackwardDofs[distance] = gridView_.indexSet().subIndex(neighbor, localIntersectionIndex, codimIntersection); |
| 306 | 36400 | const auto& backwardFacePos = getFacet_(localIntersectionIndex, neighbor).geometry().center(); | |
| 307 | 18200 | axisData_.inAxisBackwardDistances[distance] = (backwardFacePos - faceCenterPos).two_norm(); | |
| 308 | |||
| 309 | // recursion end if we added all axis data for the given stencil size (numParallelFaces) | ||
| 310 | 18200 | const auto numForwardBackwardAxisDofs = axisData_.inAxisForwardDofs.size(); | |
| 311 | if (distance >= numForwardBackwardAxisDofs-1) | ||
| 312 | 600 | return; | |
| 313 | |||
| 314 | for (const auto& intersection : intersections(gridView_, neighbor)) | ||
| 315 | { | ||
| 316 | if (intersection.indexInInside() == localIntersectionIndex && intersection.neighbor()) | ||
| 317 | { | ||
| 318 | addBackwardNeighborAxisData_(intersection.outside(), distance+1, localIntersectionIndex, backwardFacePos); | ||
| 319 | break; | ||
| 320 | } | ||
| 321 | } | ||
| 322 | } | ||
| 323 | |||
| 324 | /*! | ||
| 325 | * \brief Fills the pair data with the lateral dofs and distances | ||
| 326 | * and calls a further function to collect the parallel dofs and distances | ||
| 327 | */ | ||
| 328 | 374258 | void fillPairData_() | |
| 329 | { | ||
| 330 | // reset the pair data structs | ||
| 331 | 374258 | pairData_ = {}; | |
| 332 | |||
| 333 | // set basic global positions | ||
| 334 | 374258 | const auto& selfElementCenter = element_.geometry().center(); | |
| 335 | 374258 | const auto& selfFacetCenter = intersection_.geometry().center(); | |
| 336 | |||
| 337 | // get the inner lateral Dof Index | ||
| 338 | 374258 | SmallLocalIndexType numPairInnerLateralIdx = 0; | |
| 339 |
5/7✓ Branch 2 taken 1414532 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 28000 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 112000 times.
✓ Branch 8 taken 28000 times.
✓ Branch 1 taken 346258 times.
|
1984790 | for (const auto& innerElementIntersection : intersections(gridView_, element_)) |
| 340 | { | ||
| 341 |
5/7✓ Branch 0 taken 84000 times.
✓ Branch 1 taken 1068274 times.
✓ Branch 3 taken 56000 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 112000 times.
✗ Branch 7 not taken.
✓ Branch 2 taken 346258 times.
|
3372822 | if (facetIsNormal_(innerElementIntersection.indexInInside(), intersection_.indexInInside())) |
| 342 | { | ||
| 343 | 750016 | const auto innerElementIntersectionIdx = innerElementIntersection.indexInInside(); | |
| 344 |
1/2✓ Branch 1 taken 56000 times.
✗ Branch 2 not taken.
|
750016 | setLateralPairFirstInfo_(innerElementIntersectionIdx, element_, numPairInnerLateralIdx); |
| 345 | |||
| 346 | 2250048 | const auto distance = innerElementIntersection.geometry().center() - selfElementCenter; | |
| 347 | 750016 | pairData_[numPairInnerLateralIdx].lateralStaggeredFaceCenter = selfFacetCenter+distance; | |
| 348 | |||
| 349 | 750016 | numPairInnerLateralIdx++; | |
| 350 | } | ||
| 351 | } | ||
| 352 | |||
| 353 | // get the outer lateral Dof Index | ||
| 354 | 374258 | SmallLocalIndexType numPairOuterLateralIdx = 0; | |
| 355 |
2/2✓ Branch 0 taken 368462 times.
✓ Branch 1 taken 5796 times.
|
374258 | if (intersection_.neighbor()) |
| 356 | { | ||
| 357 | // the direct neighbor element and the respective intersection index | ||
| 358 | 362106 | const auto& directNeighborElement = intersection_.outside(); | |
| 359 | |||
| 360 |
5/7✓ Branch 2 taken 1367424 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 27400 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 109600 times.
✓ Branch 8 taken 27400 times.
✓ Branch 1 taken 334706 times.
|
2256036 | for (const auto& directNeighborElementIntersection : intersections(gridView_, directNeighborElement)) |
| 361 | { | ||
| 362 | // skip the directly neighboring face itself and its opposing one | ||
| 363 |
6/8✓ Branch 1 taken 1114918 times.
✓ Branch 2 taken 334706 times.
✓ Branch 3 taken 82200 times.
✓ Branch 4 taken 27400 times.
✓ Branch 6 taken 54800 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 109600 times.
✗ Branch 10 not taken.
|
3262554 | if (facetIsNormal_(directNeighborElementIntersection.indexInInside(), intersection_.indexInOutside())) |
| 364 | { | ||
| 365 | 725412 | const auto directNeighborElemIsIdx = directNeighborElementIntersection.indexInInside(); | |
| 366 |
1/2✓ Branch 1 taken 54800 times.
✗ Branch 2 not taken.
|
725412 | setLateralPairSecondInfo_(directNeighborElemIsIdx, directNeighborElement, numPairOuterLateralIdx); |
| 367 | 725412 | numPairOuterLateralIdx++; | |
| 368 | } | ||
| 369 | } | ||
| 370 | } | ||
| 371 | else // intersection is on boundary | ||
| 372 | { | ||
| 373 |
5/7✓ Branch 2 taken 47108 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 600 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 2400 times.
✓ Branch 8 taken 600 times.
✓ Branch 1 taken 11552 times.
|
75012 | for (const auto& intersection : intersections(gridView_, element_)) |
| 374 | { | ||
| 375 |
3/4✓ Branch 0 taken 36756 times.
✓ Branch 1 taken 12152 times.
✓ Branch 3 taken 2400 times.
✗ Branch 4 not taken.
|
85664 | if (facetIsNormal_(intersection.indexInInside(), intersection_.indexInInside())) |
| 376 | { | ||
| 377 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 24604 times.
|
24604 | assert(!pairData_[numPairOuterLateralIdx].hasOuterLateral); |
| 378 | |||
| 379 | 24604 | const auto normalDistanceoffset = selfFacetCenter - selfElementCenter; | |
| 380 | 24604 | pairData_[numPairOuterLateralIdx].lateralDistance = normalDistanceoffset.two_norm(); | |
| 381 | 24604 | numPairOuterLateralIdx++; | |
| 382 | } | ||
| 383 | } | ||
| 384 | } | ||
| 385 | |||
| 386 | // fill the pair data with the parallel dofs and distances | ||
| 387 | 374258 | const auto parallelLocalIdx = intersection_.indexInInside(); | |
| 388 | 374258 | SmallLocalIndexType numPairParallelIdx = 0; | |
| 389 |
7/9✓ Branch 2 taken 1414532 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 28000 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 84000 times.
✓ Branch 8 taken 28000 times.
✓ Branch 9 taken 112000 times.
✓ Branch 10 taken 28000 times.
✓ Branch 1 taken 346258 times.
|
1984790 | for (const auto& intersection : intersections(gridView_, element_)) |
| 390 | { | ||
| 391 |
3/4✓ Branch 0 taken 738274 times.
✓ Branch 1 taken 11742 times.
✓ Branch 3 taken 112000 times.
✗ Branch 4 not taken.
|
2248548 | if (facetIsNormal_(intersection.indexInInside(), parallelLocalIdx)) |
| 392 | { | ||
| 393 |
2/2✓ Branch 0 taken 1200 times.
✓ Branch 1 taken 1200 times.
|
1472428 | if (intersection.neighbor()) |
| 394 | { | ||
| 395 | // recursively insert parallel neighbor faces into pair data | ||
| 396 |
3/3✓ Branch 0 taken 2160 times.
✓ Branch 1 taken 713210 times.
✓ Branch 2 taken 10042 times.
|
725412 | const auto parallelAxisIdx = directionIndex(intersection); |
| 397 | 725412 | const auto localLateralIntersectionIndex = intersection.indexInInside(); | |
| 398 | |||
| 399 | /* | ||
| 400 | * ------------ | ||
| 401 | * | | | ||
| 402 | * | | | ||
| 403 | * | | | ||
| 404 | * iiiiiiiiiii*bbbbbbbbbbb | ||
| 405 | * | o zzzzzzzz | | ||
| 406 | * | o zzzzzzzz | | ||
| 407 | * | o zzzzzzzz | | ||
| 408 | * ----------------------- | ||
| 409 | * | ||
| 410 | * i:intersection,o:intersection_, b: outerIntersection, z: intersection_.outside() | ||
| 411 | */ | ||
| 412 |
2/2✓ Branch 0 taken 715130 times.
✓ Branch 1 taken 10282 times.
|
725412 | if (intersection_.neighbor()) |
| 413 |
10/12✓ Branch 1 taken 703762 times.
✓ Branch 2 taken 2604192 times.
✓ Branch 4 taken 703762 times.
✓ Branch 5 taken 1954104 times.
✓ Branch 7 taken 53674 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 53674 times.
✗ Branch 11 not taken.
✓ Branch 12 taken 53674 times.
✓ Branch 13 taken 161022 times.
✓ Branch 14 taken 214696 times.
✓ Branch 15 taken 53674 times.
|
4387434 | for (const auto& outerIntersection : intersections(gridView_, intersection_.outside())) |
| 414 |
2/2✓ Branch 0 taken 703762 times.
✓ Branch 1 taken 2115126 times.
|
2818888 | if (intersection.indexInInside() == outerIntersection.indexInInside()) |
| 415 |
1/2✓ Branch 1 taken 214696 times.
✗ Branch 2 not taken.
|
2818888 | if (!outerIntersection.neighbor()) |
| 416 | 50 | pairData_[numPairParallelIdx].hasHalfParallelNeighbor = true; | |
| 417 | |||
| 418 | /* ------------ | ||
| 419 | * | o | ||
| 420 | * | o | ||
| 421 | * | o | ||
| 422 | * iiiiiiiiiii------------ | ||
| 423 | * | zzzzzzzz b | | ||
| 424 | * | zzzzzzzz b | | ||
| 425 | * | zzzzzzzz b | | ||
| 426 | * ----------------------- | ||
| 427 | * | ||
| 428 | * i:intersection,o:intersection_, b: outerIntersection, z: intersection.outside() | ||
| 429 | */ | ||
| 430 |
3/4✓ Branch 0 taken 715130 times.
✓ Branch 1 taken 10282 times.
✓ Branch 3 taken 1126 times.
✗ Branch 4 not taken.
|
736780 | if (!intersection_.neighbor()) |
| 431 |
6/8✓ Branch 1 taken 21650 times.
✓ Branch 2 taken 83056 times.
✓ Branch 4 taken 1126 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 1126 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 4504 times.
✓ Branch 10 taken 1126 times.
|
134238 | for (const auto& outerIntersection : intersections(gridView_, intersection.outside())) |
| 432 |
2/2✓ Branch 0 taken 21650 times.
✓ Branch 1 taken 65910 times.
|
87560 | if (intersection_.indexInInside() == outerIntersection.indexInInside()) |
| 433 |
1/2✓ Branch 1 taken 4504 times.
✗ Branch 2 not taken.
|
87560 | if (outerIntersection.neighbor()) |
| 434 | 50 | pairData_[numPairParallelIdx].hasCornerParallelNeighbor = true; | |
| 435 | |||
| 436 | |||
| 437 |
2/4✓ Branch 1 taken 54800 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 54800 times.
✗ Branch 5 not taken.
|
725412 | addParallelNeighborPairData_(intersection.outside(), 0, localLateralIntersectionIndex, parallelLocalIdx, parallelAxisIdx, numPairParallelIdx); |
| 438 | } | ||
| 439 | |||
| 440 | 750016 | ++numPairParallelIdx; | |
| 441 | } | ||
| 442 | } | ||
| 443 | 374258 | } | |
| 444 | |||
| 445 | // zero distance means direct neighbor to element_ | ||
| 446 | 761612 | void addParallelNeighborPairData_(const Element& neighbor, | |
| 447 | const std::size_t distance, | ||
| 448 | const unsigned int localLateralIntersectionIndex, | ||
| 449 | const unsigned int parallelLocalIdx, | ||
| 450 | const unsigned int parallelAxisIdx, | ||
| 451 | const SmallLocalIndexType dataIdx) | ||
| 452 | { | ||
| 453 | // fill the pair data for this element in the stencil | ||
| 454 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 687812 times.
|
761612 | pairData_[dataIdx].hasParallelNeighbor.set(distance, true); |
| 455 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 761612 times.
|
761612 | pairData_[dataIdx].parallelDofs[distance] = gridView_.indexSet().subIndex(neighbor, parallelLocalIdx, codimIntersection); |
| 456 |
2/2✓ Branch 1 taken 37600 times.
✓ Branch 2 taken 36200 times.
|
761612 | pairData_[dataIdx].parallelCellWidths[distance] = setParallelPairCellWidths_(neighbor, parallelAxisIdx); |
| 457 | |||
| 458 | // recursion end if we added all parallel data for the given stencil size (numParallelFaces) | ||
| 459 | 761612 | const auto numParallelFaces = pairData_[0].parallelCellWidths.size(); | |
| 460 |
2/2✓ Branch 0 taken 37600 times.
✓ Branch 1 taken 36200 times.
|
73800 | if (distance >= numParallelFaces-1) |
| 461 | return; | ||
| 462 | |||
| 463 | // continue with neighbor's neighbor element in the same direction, if we find one | ||
| 464 |
7/12✓ Branch 2 taken 43600 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 41200 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 35200 times.
✓ Branch 8 taken 52800 times.
✓ Branch 10 taken 52800 times.
✗ Branch 11 not taken.
✓ Branch 12 taken 88000 times.
✗ Branch 13 not taken.
✗ Branch 1 not taken.
✓ Branch 4 taken 2400 times.
|
239600 | for (const auto& lateralIntersection : intersections(gridView_, neighbor)) |
| 465 | { | ||
| 466 | // we only expect to find one intersection matching this condition | ||
| 467 |
2/2✓ Branch 0 taken 37600 times.
✓ Branch 1 taken 58800 times.
|
96400 | if (lateralIntersection.indexInInside() == localLateralIntersectionIndex) |
| 468 | { | ||
| 469 | // if there is no neighbor recursion ends here | ||
| 470 | 37600 | if (lateralIntersection.neighbor()) | |
| 471 |
2/4✓ Branch 1 taken 34400 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 34400 times.
✗ Branch 5 not taken.
|
36200 | addParallelNeighborPairData_(lateralIntersection.outside(), distance+1, |
| 472 | localLateralIntersectionIndex, parallelLocalIdx, parallelAxisIdx, dataIdx); | ||
| 473 | break; | ||
| 474 | } | ||
| 475 | } | ||
| 476 | } | ||
| 477 | |||
| 478 | /*! | ||
| 479 | * \brief Returns the local opposing intersection index | ||
| 480 | * | ||
| 481 | * \param idx The local index of the intersection itself | ||
| 482 | */ | ||
| 483 | 3765830 | int localOppositeIdx_(const int idx) const | |
| 484 | { | ||
| 485 |
12/12✓ Branch 0 taken 542637 times.
✓ Branch 1 taken 542637 times.
✓ Branch 2 taken 544754 times.
✓ Branch 3 taken 544714 times.
✓ Branch 4 taken 45323 times.
✓ Branch 5 taken 45283 times.
✓ Branch 6 taken 534237 times.
✓ Branch 7 taken 534237 times.
✓ Branch 8 taken 206629 times.
✓ Branch 9 taken 206629 times.
✓ Branch 10 taken 9375 times.
✓ Branch 11 taken 9375 times.
|
3765830 | return (idx % 2) ? (idx - 1) : (idx + 1); |
| 486 | } | ||
| 487 | |||
| 488 | /*! | ||
| 489 | * \brief Returns true if the intersection lies normal to another given intersection | ||
| 490 | * | ||
| 491 | * \param selfIdx The local index of the intersection itself | ||
| 492 | * \param otherIdx The local index of the other intersection | ||
| 493 | */ | ||
| 494 | 4495596 | bool facetIsNormal_(const int selfIdx, const int otherIdx) const | |
| 495 | { | ||
| 496 |
16/16✓ Branch 0 taken 778016 times.
✓ Branch 1 taken 374258 times.
✓ Branch 2 taken 1096274 times.
✓ Branch 3 taken 374258 times.
✓ Branch 4 taken 752812 times.
✓ Branch 5 taken 362106 times.
✓ Branch 6 taken 1060118 times.
✓ Branch 7 taken 362106 times.
✓ Branch 8 taken 25204 times.
✓ Branch 9 taken 12152 times.
✓ Branch 10 taken 36156 times.
✓ Branch 11 taken 12152 times.
✓ Branch 12 taken 778016 times.
✓ Branch 13 taken 374258 times.
✓ Branch 14 taken 56000 times.
✓ Branch 15 taken 28000 times.
|
6481886 | return !(selfIdx == otherIdx || localOppositeIdx_(selfIdx) == otherIdx); |
| 497 | } | ||
| 498 | |||
| 499 | 1118770 | auto getFacet_(const int localFacetIdx, const Element& element) const | |
| 500 | { | ||
| 501 | 1173570 | return element.template subEntity <1> (localFacetIdx); | |
| 502 | } | ||
| 503 | |||
| 504 | //! Sets the information about the lateral faces (within the element) | ||
| 505 | 750016 | void setLateralPairFirstInfo_(const int isIdx, const Element& element, const int numPairsIdx) | |
| 506 | { | ||
| 507 | // store the inner lateral dofIdx | ||
| 508 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 750016 times.
|
750016 | pairData_[numPairsIdx].lateralPair.first = gridView_.indexSet().subIndex(element, isIdx, codimIntersection); |
| 509 | |||
| 510 | // store the local lateral facet index | ||
| 511 | 750016 | pairData_[numPairsIdx].localLateralFaceIdx = isIdx; | |
| 512 | 750016 | } | |
| 513 | |||
| 514 | //! Sets the information about the lateral faces (in the neighbor element) | ||
| 515 | 725412 | void setLateralPairSecondInfo_(const int isIdx, const Element& element, const int numPairsIdx) | |
| 516 | { | ||
| 517 | // store the dofIdx | ||
| 518 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 725412 times.
|
725412 | pairData_[numPairsIdx].hasOuterLateral = true; |
| 519 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 725412 times.
|
725412 | pairData_[numPairsIdx].lateralPair.second = gridView_.indexSet().subIndex(element, isIdx, codimIntersection); |
| 520 | |||
| 521 | // set basic global positions | ||
| 522 | 725412 | const auto& selfFacetCenter = intersection_.geometry().center(); | |
| 523 | 725412 | const auto& selfElementCenter = element_.geometry().center(); | |
| 524 | |||
| 525 | 725412 | const auto& neighborElement = intersection_.outside(); | |
| 526 | 725412 | const auto& neighborElementCenter = neighborElement.geometry().center(); | |
| 527 | 1450824 | const auto& neighborFacetCenter = getFacet_(intersection_.indexInOutside(), neighborElement).geometry().center(); | |
| 528 | |||
| 529 | 1450824 | const Scalar insideLateralDistance = (selfFacetCenter - selfElementCenter).two_norm(); | |
| 530 | 725412 | const Scalar outsideLateralDistance = (neighborFacetCenter - neighborElementCenter).two_norm(); | |
| 531 | |||
| 532 | 725412 | pairData_[numPairsIdx].lateralDistance = insideLateralDistance + outsideLateralDistance; | |
| 533 | 725412 | } | |
| 534 | |||
| 535 | //! Sets the information about the parallel distances | ||
| 536 | 761612 | Scalar setParallelPairCellWidths_(const Element& element, const int parallelAxisIdx) const | |
| 537 | { | ||
| 538 | 761612 | std::vector<GlobalPosition> faces; | |
| 539 |
1/2✓ Branch 1 taken 761612 times.
✗ Branch 2 not taken.
|
761612 | faces.reserve(numFacets); |
| 540 |
6/8✓ Branch 1 taken 761612 times.
✓ Branch 2 taken 2698048 times.
✓ Branch 4 taken 89200 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 89200 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 356800 times.
✓ Branch 11 taken 89200 times.
|
4173260 | for (const auto& intersection : intersections(gridView_, element)) |
| 541 |
2/5✓ Branch 1 taken 3054848 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 356800 times.
✗ Branch 5 not taken.
✗ Branch 3 not taken.
|
6109696 | faces.push_back(intersection.geometry().center()); |
| 542 | |||
| 543 |
3/4✓ Branch 0 taken 379584 times.
✓ Branch 1 taken 380628 times.
✓ Branch 2 taken 1400 times.
✗ Branch 3 not taken.
|
761612 | switch (parallelAxisIdx) |
| 544 | { | ||
| 545 | 379584 | case 0: | |
| 546 | 1138752 | return (faces[1] - faces[0]).two_norm(); | |
| 547 | 380628 | case 1: | |
| 548 | 1141884 | return (faces[3] - faces[2]).two_norm(); | |
| 549 | 1400 | case 2: | |
| 550 | { | ||
| 551 | ✗ | assert(dim == 3); | |
| 552 | 4200 | return (faces[5] - faces[4]).two_norm(); | |
| 553 | } | ||
| 554 | ✗ | default: | |
| 555 |
1/22✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
✗ Branch 16 not taken.
✗ Branch 17 not taken.
✗ Branch 19 not taken.
✗ Branch 20 not taken.
✗ Branch 22 not taken.
✗ Branch 23 not taken.
✗ Branch 25 not taken.
✗ Branch 26 not taken.
✗ Branch 28 not taken.
✗ Branch 29 not taken.
✓ Branch 33 taken 757412 times.
✗ Branch 34 not taken.
|
757412 | DUNE_THROW(Dune::InvalidStateException, "Something went terribly wrong"); |
| 556 | } | ||
| 557 | 761612 | } | |
| 558 | |||
| 559 | Intersection intersection_; //!< The intersection of interest | ||
| 560 | const Element& element_; //!< The respective element | ||
| 561 | const GridView gridView_; //!< The grid view | ||
| 562 | AxisData axisData_; //!< Data related to forward and backward faces | ||
| 563 | std::array<PairData, numPairs> pairData_; //!< Collection of pair information related to lateral and parallel faces | ||
| 564 | }; | ||
| 565 | |||
| 566 | } // end namespace Dumux | ||
| 567 | |||
| 568 | #endif | ||
| 569 |