| 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 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 | 717136 | 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 | 1064322 | inline static unsigned int directionIndex(Vector&& vector) | |
| 122 | { | ||
| 123 | 1064322 | const auto eps = 1e-8; | |
| 124 |
1/32✗ 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 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
✗ Branch 15 not taken.
✗ Branch 16 not taken.
✗ Branch 17 not taken.
✗ Branch 18 not taken.
✗ Branch 19 not taken.
✗ Branch 20 not taken.
✗ Branch 21 not taken.
✗ Branch 22 not taken.
✗ Branch 23 not taken.
✗ Branch 24 not taken.
✗ Branch 25 not taken.
✗ Branch 26 not taken.
✗ Branch 27 not taken.
✗ Branch 28 not taken.
✗ Branch 29 not taken.
✗ Branch 34 not taken.
✓ Branch 35 taken 1064322 times.
|
4257288 | 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 | 89517 | FreeFlowStaggeredGeometryHelper(const Element& element, const GridView& gridView) | |
| 157 | : element_(element) | ||
| 158 | 254676 | , gridView_(gridView) | |
| 159 | 89517 | { } | |
| 160 | |||
| 161 | //! update the local face | ||
| 162 | template<class IntersectionMapper> | ||
| 163 | ✗ | void updateLocalFace(const IntersectionMapper&, const Intersection& intersection) | |
| 164 | { | ||
| 165 | ✗ | intersection_ = intersection; | |
| 166 | ✗ | fillAxisData_(); | |
| 167 | ✗ | fillPairData_(); | |
| 168 | ✗ | } | |
| 169 | |||
| 170 | /*! | ||
| 171 | * \brief Returns the local index of the face (i.e. the intersection) | ||
| 172 | */ | ||
| 173 | SmallLocalIndexType localFaceIndex() const | ||
| 174 | { | ||
| 175 |
4/4✓ Branch 2 taken 355660 times.
✓ Branch 3 taken 2658 times.
✓ Branch 4 taken 355660 times.
✓ Branch 5 taken 2658 times.
|
1433272 | return intersection_.indexInInside(); |
| 176 | } | ||
| 177 | |||
| 178 | /*! | ||
| 179 | * \brief Returns a copy of the axis data | ||
| 180 | */ | ||
| 181 | ✗ | AxisData axisData() const | |
| 182 | { | ||
| 183 | 18750 | return axisData_; | |
| 184 | } | ||
| 185 | |||
| 186 | /*! | ||
| 187 | * \brief Returns a copy of the pair data | ||
| 188 | */ | ||
| 189 | std::array<PairData, numPairs> pairData() const | ||
| 190 | { | ||
| 191 | 358318 | return pairData_; | |
| 192 | } | ||
| 193 | |||
| 194 | /*! | ||
| 195 | * \brief Returns the direction index of the primary facet (0 = x, 1 = y, 2 = z) | ||
| 196 | */ | ||
| 197 | 358318 | unsigned int directionIndex() const | |
| 198 | { | ||
| 199 | 688636 | 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 | ✗ | unsigned int directionIndex(const Intersection& intersection) const | |
| 206 | { | ||
| 207 | ✗ | 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 | 358318 | void fillAxisData_() | |
| 217 | { | ||
| 218 | if constexpr (useHigherOrder) | ||
| 219 | 18750 | fillOuterAxisData_(); | |
| 220 | |||
| 221 |
2/2✓ Branch 0 taken 179159 times.
✓ Branch 1 taken 179159 times.
|
358318 | const auto inIdx = intersection_.indexInInside(); |
| 222 | 716636 | const auto oppIdx = localOppositeIdx_(inIdx); | |
| 223 | |||
| 224 | // Set the self Dof | ||
| 225 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 358318 times.
|
358318 | 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 358318 times.
|
358318 | axisData_.oppositeDof = gridView_.indexSet().subIndex(intersection_.inside(), oppIdx, codimIntersection); |
| 229 | |||
| 230 | // Set the Self to Opposite Distance | ||
| 231 | 716636 | const auto self = getFacet_(inIdx, element_); | |
| 232 | 716636 | const auto opposite = getFacet_(oppIdx, element_); | |
| 233 | 1432522 | axisData_.selfToOppositeDistance = (self.geometry().center() - opposite.geometry().center()).two_norm(); | |
| 234 | |||
| 235 | 358318 | } | |
| 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 | 18800 | addForwardNeighborAxisData_(intersection_.outside(), 0, thisFaceLocalIdx, faceCenterPos); | |
| 252 | |||
| 253 | // Set the backward dofs | ||
| 254 | 37500 | const auto oppositeFaceLocalIdx = localOppositeIdx_(thisFaceLocalIdx); | |
| 255 | 18750 | const auto& oppositeFaceCenterPos = getFacet_(oppositeFaceLocalIdx, element_).geometry().center(); | |
| 256 |
10/13✓ Branch 1 taken 150 times.
✓ Branch 2 taken 21000 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 18000 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 45600 times.
✓ Branch 8 taken 400 times.
✓ Branch 9 taken 45600 times.
✓ Branch 10 taken 400 times.
✓ Branch 11 taken 18000 times.
✓ Branch 12 taken 27600 times.
✓ Branch 14 taken 28000 times.
✗ Branch 15 not taken.
|
102750 | for (const auto& intersection : intersections(gridView_, element_)) |
| 257 | { | ||
| 258 |
6/6✓ Branch 0 taken 18750 times.
✓ Branch 1 taken 29850 times.
✓ Branch 2 taken 18750 times.
✓ Branch 3 taken 29850 times.
✓ Branch 4 taken 18675 times.
✓ Branch 5 taken 75 times.
|
97200 | 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.
|
18800 | 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 | 18200 | const auto& forwardFacePos = getFacet_(localIntersectionIndex, neighbor).geometry().center(); | |
| 278 | 35800 | 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 | 18200 | const auto& backwardFacePos = getFacet_(localIntersectionIndex, neighbor).geometry().center(); | |
| 307 | 35800 | 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 | 358318 | void fillPairData_() | |
| 329 | { | ||
| 330 | // reset the pair data structs | ||
| 331 | 358318 | pairData_ = {}; | |
| 332 | |||
| 333 | // set basic global positions | ||
| 334 | 358318 | const auto& selfElementCenter = element_.geometry().center(); | |
| 335 | 358318 | const auto& selfFacetCenter = intersection_.geometry().center(); | |
| 336 | |||
| 337 | // get the inner lateral Dof Index | ||
| 338 | 358318 | SmallLocalIndexType numPairInnerLateralIdx = 0; | |
| 339 |
10/13✓ Branch 1 taken 330318 times.
✓ Branch 2 taken 1350772 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 9 taken 112000 times.
✓ Branch 10 taken 28000 times.
✓ Branch 11 taken 84000 times.
✓ Branch 12 taken 28000 times.
✓ Branch 14 taken 112000 times.
✗ Branch 15 not taken.
|
1849090 | for (const auto& innerElementIntersection : intersections(gridView_, element_)) |
| 340 | { | ||
| 341 |
4/4✓ Branch 0 taken 1076454 times.
✓ Branch 1 taken 358318 times.
✓ Branch 2 taken 1076454 times.
✓ Branch 3 taken 358318 times.
|
2869544 | if (facetIsNormal_(innerElementIntersection.indexInInside(), intersection_.indexInInside())) |
| 342 | { | ||
| 343 |
1/2✓ Branch 1 taken 56000 times.
✗ Branch 2 not taken.
|
718136 | const auto innerElementIntersectionIdx = innerElementIntersection.indexInInside(); |
| 344 |
1/2✓ Branch 1 taken 56000 times.
✗ Branch 2 not taken.
|
718136 | setLateralPairFirstInfo_(innerElementIntersectionIdx, element_, numPairInnerLateralIdx); |
| 345 | |||
| 346 | 1436272 | const auto distance = innerElementIntersection.geometry().center() - selfElementCenter; | |
| 347 | 1436272 | pairData_[numPairInnerLateralIdx].lateralStaggeredFaceCenter = selfFacetCenter+distance; | |
| 348 | |||
| 349 | 718136 | numPairInnerLateralIdx++; | |
| 350 | } | ||
| 351 | } | ||
| 352 | |||
| 353 | // get the outer lateral Dof Index | ||
| 354 | 358318 | SmallLocalIndexType numPairOuterLateralIdx = 0; | |
| 355 |
2/2✓ Branch 0 taken 355660 times.
✓ Branch 1 taken 2658 times.
|
358318 | if (intersection_.neighbor()) |
| 356 | { | ||
| 357 | // the direct neighbor element and the respective intersection index | ||
| 358 | 352402 | const auto& directNeighborElement = intersection_.outside(); | |
| 359 | |||
| 360 |
9/13✓ Branch 1 taken 325002 times.
✓ Branch 2 taken 1328608 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 9 taken 109600 times.
✓ Branch 10 taken 27400 times.
✓ Branch 12 taken 109600 times.
✗ Branch 13 not taken.
✓ Branch 15 taken 109600 times.
✗ Branch 16 not taken.
|
2143012 | for (const auto& directNeighborElementIntersection : intersections(gridView_, directNeighborElement)) |
| 361 | { | ||
| 362 | // skip the directly neighboring face itself and its opposing one | ||
| 363 |
5/5✓ Branch 0 taken 976206 times.
✓ Branch 1 taken 434602 times.
✓ Branch 2 taken 976206 times.
✓ Branch 3 taken 407202 times.
✓ Branch 4 taken 27400 times.
|
2712016 | if (facetIsNormal_(directNeighborElementIntersection.indexInInside(), intersection_.indexInOutside())) |
| 364 | { | ||
| 365 |
1/2✓ Branch 1 taken 54800 times.
✗ Branch 2 not taken.
|
706004 | const auto directNeighborElemIsIdx = directNeighborElementIntersection.indexInInside(); |
| 366 |
1/2✓ Branch 1 taken 54800 times.
✗ Branch 2 not taken.
|
706004 | setLateralPairSecondInfo_(directNeighborElemIsIdx, directNeighborElement, numPairOuterLateralIdx); |
| 367 | 706004 | numPairOuterLateralIdx++; | |
| 368 | } | ||
| 369 | } | ||
| 370 | } | ||
| 371 | else // intersection is on boundary | ||
| 372 | { | ||
| 373 |
10/13✓ Branch 1 taken 5316 times.
✓ Branch 2 taken 22164 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 9 taken 2400 times.
✓ Branch 10 taken 600 times.
✓ Branch 11 taken 1800 times.
✓ Branch 12 taken 600 times.
✓ Branch 14 taken 2400 times.
✗ Branch 15 not taken.
|
36396 | for (const auto& intersection : intersections(gridView_, element_)) |
| 374 | { | ||
| 375 |
4/4✓ Branch 0 taken 18048 times.
✓ Branch 1 taken 5916 times.
✓ Branch 2 taken 18048 times.
✓ Branch 3 taken 5916 times.
|
47928 | if (facetIsNormal_(intersection.indexInInside(), intersection_.indexInInside())) |
| 376 | { | ||
| 377 |
2/4✗ Branch 0 not taken.
✓ Branch 1 taken 12132 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 12132 times.
|
24264 | assert(!pairData_[numPairOuterLateralIdx].hasOuterLateral); |
| 378 | |||
| 379 | 12132 | const auto normalDistanceoffset = selfFacetCenter - selfElementCenter; | |
| 380 | 23664 | pairData_[numPairOuterLateralIdx].lateralDistance = normalDistanceoffset.two_norm(); | |
| 381 | 12132 | numPairOuterLateralIdx++; | |
| 382 | } | ||
| 383 | } | ||
| 384 | } | ||
| 385 | |||
| 386 | // fill the pair data with the parallel dofs and distances | ||
| 387 | 358318 | const auto parallelLocalIdx = intersection_.indexInInside(); | |
| 388 | 358318 | SmallLocalIndexType numPairParallelIdx = 0; | |
| 389 |
10/13✓ Branch 1 taken 330318 times.
✓ Branch 2 taken 1350772 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 9 taken 112000 times.
✓ Branch 10 taken 28000 times.
✓ Branch 11 taken 84000 times.
✓ Branch 12 taken 28000 times.
✓ Branch 14 taken 112000 times.
✗ Branch 15 not taken.
|
1849090 | for (const auto& intersection : intersections(gridView_, element_)) |
| 390 | { | ||
| 391 |
2/2✓ Branch 0 taken 1076454 times.
✓ Branch 1 taken 358318 times.
|
1434772 | if (facetIsNormal_(intersection.indexInInside(), parallelLocalIdx)) |
| 392 | { | ||
| 393 |
2/2✓ Branch 0 taken 712670 times.
✓ Branch 1 taken 5466 times.
|
718136 | if (intersection.neighbor()) |
| 394 | { | ||
| 395 | // recursively insert parallel neighbor faces into pair data | ||
| 396 | 706004 | const auto parallelAxisIdx = directionIndex(intersection); | |
| 397 |
2/2✓ Branch 0 taken 700786 times.
✓ Branch 1 taken 5218 times.
|
706004 | 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 700786 times.
✓ Branch 1 taken 5218 times.
|
706004 | if (intersection_.neighbor()) |
| 413 |
13/17✓ Branch 1 taken 53674 times.
✓ Branch 2 taken 640768 times.
✓ Branch 3 taken 2566912 times.
✓ Branch 4 taken 53674 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 53674 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 53674 times.
✗ Branch 11 not taken.
✓ Branch 12 taken 214696 times.
✓ Branch 13 taken 53674 times.
✓ Branch 14 taken 214696 times.
✓ Branch 15 taken 53674 times.
✓ Branch 16 taken 53674 times.
✓ Branch 17 taken 161022 times.
✓ Branch 19 taken 214696 times.
✗ Branch 20 not taken.
|
4864934 | for (const auto& outerIntersection : intersections(gridView_, intersection_.outside())) |
| 414 |
6/6✓ Branch 0 taken 694442 times.
✓ Branch 1 taken 2087166 times.
✓ Branch 2 taken 694442 times.
✓ Branch 3 taken 2087166 times.
✓ Branch 4 taken 694442 times.
✓ Branch 5 taken 2087166 times.
|
8344824 | if (intersection.indexInInside() == outerIntersection.indexInInside()) |
| 415 |
1/2✓ Branch 0 taken 694442 times.
✗ Branch 1 not taken.
|
694442 | if (!outerIntersection.neighbor()) |
| 416 | 100 | 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 |
2/2✓ Branch 0 taken 700786 times.
✓ Branch 1 taken 5218 times.
|
706004 | if (!intersection_.neighbor()) |
| 431 |
13/17✓ Branch 1 taken 1126 times.
✓ Branch 2 taken 10436 times.
✓ Branch 3 taken 42704 times.
✓ Branch 4 taken 1126 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 1126 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 1126 times.
✗ Branch 11 not taken.
✓ Branch 12 taken 4504 times.
✓ Branch 13 taken 1126 times.
✓ Branch 14 taken 4504 times.
✓ Branch 15 taken 1126 times.
✓ Branch 16 taken 1126 times.
✓ Branch 17 taken 3378 times.
✓ Branch 19 taken 4504 times.
✗ Branch 20 not taken.
|
81894 | for (const auto& outerIntersection : intersections(gridView_, intersection.outside())) |
| 432 |
6/6✓ Branch 0 taken 11562 times.
✓ Branch 1 taken 35646 times.
✓ Branch 2 taken 11562 times.
✓ Branch 3 taken 35646 times.
✓ Branch 4 taken 11562 times.
✓ Branch 5 taken 35646 times.
|
141624 | if (intersection_.indexInInside() == outerIntersection.indexInInside()) |
| 433 |
2/2✓ Branch 0 taken 6344 times.
✓ Branch 1 taken 5218 times.
|
11562 | if (outerIntersection.neighbor()) |
| 434 | 100 | 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.
|
1357208 | addParallelNeighborPairData_(intersection.outside(), 0, localLateralIntersectionIndex, parallelLocalIdx, parallelAxisIdx, numPairParallelIdx); |
| 438 | } | ||
| 439 | |||
| 440 | 718136 | ++numPairParallelIdx; | |
| 441 | } | ||
| 442 | } | ||
| 443 | 358318 | } | |
| 444 | |||
| 445 | // zero distance means direct neighbor to element_ | ||
| 446 | 742204 | 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 |
2/4✗ Branch 0 not taken.
✓ Branch 1 taken 668404 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 668404 times.
|
1484408 | pairData_[dataIdx].hasParallelNeighbor.set(distance, true); |
| 455 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 742204 times.
|
742204 | pairData_[dataIdx].parallelDofs[distance] = gridView_.indexSet().subIndex(neighbor, parallelLocalIdx, codimIntersection); |
| 456 |
4/4✓ Branch 1 taken 37600 times.
✓ Branch 2 taken 36200 times.
✓ Branch 3 taken 37600 times.
✓ Branch 4 taken 36200 times.
|
742204 | pairData_[dataIdx].parallelCellWidths[distance] = setParallelPairCellWidths_(neighbor, parallelAxisIdx); |
| 457 | |||
| 458 | // recursion end if we added all parallel data for the given stencil size (numParallelFaces) | ||
| 459 |
2/2✓ Branch 0 taken 37600 times.
✓ Branch 1 taken 36200 times.
|
1410608 | 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/13✗ Branch 1 not taken.
✓ Branch 2 taken 43600 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 35200 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 88000 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 88000 times.
✗ Branch 10 not taken.
✓ Branch 11 taken 35200 times.
✓ Branch 12 taken 52800 times.
✓ Branch 14 taken 52800 times.
✗ Branch 15 not taken.
|
202000 | for (const auto& lateralIntersection : intersections(gridView_, neighbor)) |
| 465 | { | ||
| 466 | // we only expect to find one intersection matching this condition | ||
| 467 |
4/4✓ Branch 0 taken 37600 times.
✓ Branch 1 taken 58800 times.
✓ Branch 2 taken 37600 times.
✓ Branch 3 taken 58800 times.
|
192800 | if (lateralIntersection.indexInInside() == localLateralIntersectionIndex) |
| 468 | { | ||
| 469 | // if there is no neighbor recursion ends here | ||
| 470 |
2/2✓ Branch 0 taken 37300 times.
✓ Branch 1 taken 300 times.
|
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.
|
38000 | 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 | ✗ | int localOppositeIdx_(const int idx) const | |
| 484 | { | ||
| 485 |
4/6✗ Branch 0 not taken.
✗ Branch 1 not taken.
✓ Branch 2 taken 179159 times.
✓ Branch 3 taken 179159 times.
✓ Branch 4 taken 9375 times.
✓ Branch 5 taken 9375 times.
|
3606430 | 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 | ✗ | bool facetIsNormal_(const int selfIdx, const int otherIdx) const | |
| 495 | { | ||
| 496 |
24/30✗ 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 taken 1076454 times.
✓ Branch 7 taken 358318 times.
✓ Branch 8 taken 538227 times.
✓ Branch 9 taken 538227 times.
✓ Branch 10 taken 718136 times.
✓ Branch 11 taken 358318 times.
✓ Branch 12 taken 1058406 times.
✓ Branch 13 taken 352402 times.
✓ Branch 14 taken 529203 times.
✓ Branch 15 taken 529203 times.
✓ Branch 16 taken 706004 times.
✓ Branch 17 taken 352402 times.
✓ Branch 18 taken 18048 times.
✓ Branch 19 taken 5916 times.
✓ Branch 20 taken 9024 times.
✓ Branch 21 taken 9024 times.
✓ Branch 22 taken 12132 times.
✓ Branch 23 taken 5916 times.
✓ Branch 24 taken 1076454 times.
✓ Branch 25 taken 358318 times.
✓ Branch 26 taken 538227 times.
✓ Branch 27 taken 538227 times.
✓ Branch 28 taken 718136 times.
✓ Branch 29 taken 358318 times.
|
7533678 | return !(selfIdx == otherIdx || localOppositeIdx_(selfIdx) == otherIdx); |
| 497 | } | ||
| 498 | |||
| 499 | ✗ | auto getFacet_(const int localFacetIdx, const Element& element) const | |
| 500 | { | ||
| 501 | 432218 | return element.template subEntity <1> (localFacetIdx); | |
| 502 | } | ||
| 503 | |||
| 504 | //! Sets the information about the lateral faces (within the element) | ||
| 505 | 718136 | 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 718136 times.
|
718136 | pairData_[numPairsIdx].lateralPair.first = gridView_.indexSet().subIndex(element, isIdx, codimIntersection); |
| 509 | |||
| 510 | // store the local lateral facet index | ||
| 511 | 1436272 | pairData_[numPairsIdx].localLateralFaceIdx = isIdx; | |
| 512 | 718136 | } | |
| 513 | |||
| 514 | //! Sets the information about the lateral faces (in the neighbor element) | ||
| 515 | 706004 | 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 706004 times.
|
706004 | pairData_[numPairsIdx].hasOuterLateral = true; |
| 519 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 706004 times.
|
706004 | pairData_[numPairsIdx].lateralPair.second = gridView_.indexSet().subIndex(element, isIdx, codimIntersection); |
| 520 | |||
| 521 | // set basic global positions | ||
| 522 | 706004 | const auto& selfFacetCenter = intersection_.geometry().center(); | |
| 523 | 706004 | const auto& selfElementCenter = element_.geometry().center(); | |
| 524 | |||
| 525 | 706004 | const auto& neighborElement = intersection_.outside(); | |
| 526 | 706004 | const auto& neighborElementCenter = neighborElement.geometry().center(); | |
| 527 | 1357208 | const auto& neighborFacetCenter = getFacet_(intersection_.indexInOutside(), neighborElement).geometry().center(); | |
| 528 | |||
| 529 | 1409608 | const Scalar insideLateralDistance = (selfFacetCenter - selfElementCenter).two_norm(); | |
| 530 | 1409608 | const Scalar outsideLateralDistance = (neighborFacetCenter - neighborElementCenter).two_norm(); | |
| 531 | |||
| 532 | 1412008 | pairData_[numPairsIdx].lateralDistance = insideLateralDistance + outsideLateralDistance; | |
| 533 | 706004 | } | |
| 534 | |||
| 535 | //! Sets the information about the parallel distances | ||
| 536 | 742204 | Scalar setParallelPairCellWidths_(const Element& element, const int parallelAxisIdx) const | |
| 537 | { | ||
| 538 |
2/6✓ Branch 1 taken 742204 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 742204 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
|
2226612 | std::vector<GlobalPosition> faces; |
| 539 |
1/2✓ Branch 1 taken 742204 times.
✗ Branch 2 not taken.
|
742204 | faces.reserve(numFacets); |
| 540 |
8/10✓ Branch 1 taken 742204 times.
✓ Branch 2 taken 2620416 times.
✓ Branch 4 taken 89200 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 89200 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 356800 times.
✓ Branch 10 taken 89200 times.
✓ Branch 11 taken 356800 times.
✓ Branch 12 taken 89200 times.
|
4254620 | for (const auto& intersection : intersections(gridView_, element)) |
| 541 |
2/5✓ Branch 2 taken 2977216 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 356800 times.
✗ Branch 6 not taken.
|
5954432 | faces.push_back(intersection.geometry().center()); |
| 542 | |||
| 543 |
3/4✓ Branch 0 taken 370000 times.
✓ Branch 1 taken 370804 times.
✓ Branch 2 taken 1400 times.
✗ Branch 3 not taken.
|
742204 | switch (parallelAxisIdx) |
| 544 | { | ||
| 545 | 370000 | case 0: | |
| 546 | 1478600 | return (faces[1] - faces[0]).two_norm(); | |
| 547 | 370804 | case 1: | |
| 548 | 1852620 | return (faces[3] - faces[2]).two_norm(); | |
| 549 | 1400 | case 2: | |
| 550 | { | ||
| 551 | ✗ | assert(dim == 3); | |
| 552 | 5600 | return (faces[5] - faces[4]).two_norm(); | |
| 553 | } | ||
| 554 | ✗ | default: | |
| 555 | ✗ | DUNE_THROW(Dune::InvalidStateException, "Something went terribly wrong"); | |
| 556 | } | ||
| 557 | } | ||
| 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 |