| 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 FreeFlowPoreNetworkCoupling | ||
| 10 | * \copydoc Dumux::StaggeredFreeFlowPoreNetworkCouplingMapper | ||
| 11 | */ | ||
| 12 | |||
| 13 | #ifndef DUMUX_MULTIDOMAIN_FREEFLOW_POROUSMEDIUM_COUPLINGMAPPER_HH | ||
| 14 | #define DUMUX_MULTIDOMAIN_FREEFLOW_POROUSMEDIUM_COUPLINGMAPPER_HH | ||
| 15 | |||
| 16 | #include <algorithm> | ||
| 17 | #include <iostream> | ||
| 18 | #include <map> | ||
| 19 | #include <unordered_map> | ||
| 20 | #include <tuple> | ||
| 21 | #include <vector> | ||
| 22 | |||
| 23 | #include <dune/common/timer.hh> | ||
| 24 | #include <dune/common/exceptions.hh> | ||
| 25 | #include <dune/common/indices.hh> | ||
| 26 | #include <dune/common/reservedvector.hh> | ||
| 27 | #include <dune/geometry/axisalignedcubegeometry.hh> | ||
| 28 | |||
| 29 | #include <dumux/geometry/diameter.hh> | ||
| 30 | #include <dumux/geometry/intersectingentities.hh> | ||
| 31 | #include <dumux/discretization/method.hh> | ||
| 32 | #include <dumux/discretization/facecentered/staggered/normalaxis.hh> | ||
| 33 | |||
| 34 | namespace Dumux { | ||
| 35 | |||
| 36 | /*! | ||
| 37 | * \ingroup FreeFlowPoreNetworkCoupling | ||
| 38 | * \brief Coupling mapper for staggered free-flow and pore-network models. | ||
| 39 | */ | ||
| 40 | // template<class MDTraits, class CouplingManager> | ||
| 41 | class StaggeredFreeFlowPoreNetworkCouplingMapper | ||
| 42 | { | ||
| 43 | using MapType = std::unordered_map<std::size_t, std::vector<std::size_t>>; | ||
| 44 | |||
| 45 | public: | ||
| 46 | /*! | ||
| 47 | * \brief Main update routine | ||
| 48 | */ | ||
| 49 | template<class FreeFlowMomentumGridGeometry, class FreeFlowMassGridGeometry, class PoreNetworkGridGeometry> | ||
| 50 | 6 | void update(const FreeFlowMomentumGridGeometry& ffMomentumGridGeometry, | |
| 51 | const FreeFlowMassGridGeometry& ffMassGridGeometry, | ||
| 52 | const PoreNetworkGridGeometry& pnmGridGeometry) | ||
| 53 | { | ||
| 54 | 6 | clear_(); | |
| 55 | 6 | resize_(ffMomentumGridGeometry, pnmGridGeometry); | |
| 56 | 6 | Dune::Timer watch; | |
| 57 | 6 | std::cout << "Initializing the coupling map..." << std::endl; | |
| 58 | |||
| 59 | 6 | auto ffFvGeometry = localView(ffMomentumGridGeometry); | |
| 60 | 6 | auto pnmFvGeometry = localView(pnmGridGeometry); | |
| 61 | |||
| 62 | using GlobalPosition = typename FreeFlowMomentumGridGeometry::GridView::template Codim<0>::Entity::Geometry::GlobalCoordinate; | ||
| 63 | |||
| 64 |
2/2✓ Branch 2 taken 15 times.
✓ Branch 3 taken 6 times.
|
36 | for (const auto& pnmElement : elements(pnmGridGeometry.gridView())) |
| 65 | { | ||
| 66 | 15 | const auto pnmElementIdx = pnmGridGeometry.elementMapper().index(pnmElement); | |
| 67 | 15 | pnmFvGeometry.bindElement(pnmElement); | |
| 68 |
4/4✓ Branch 0 taken 7 times.
✓ Branch 1 taken 23 times.
✓ Branch 4 taken 30 times.
✓ Branch 5 taken 15 times.
|
53 | for (const auto& pnmScv : scvs(pnmFvGeometry)) |
| 69 | { | ||
| 70 | // skip the dof if it is not on the boundary | ||
| 71 |
2/2✓ Branch 0 taken 7 times.
✓ Branch 1 taken 23 times.
|
30 | if (!pnmGridGeometry.dofOnBoundary(pnmScv.dofIndex())) |
| 72 | 22 | continue; | |
| 73 | |||
| 74 | // get the intersection bulk element | ||
| 75 | 23 | const auto pnmPos = pnmScv.dofPosition(); | |
| 76 |
1/2✓ Branch 1 taken 23 times.
✗ Branch 2 not taken.
|
23 | const auto pnmDofIdx = pnmScv.dofIndex(); |
| 77 |
1/2✓ Branch 1 taken 23 times.
✗ Branch 2 not taken.
|
23 | const auto& otherPNMScv = pnmFvGeometry.scv(1 - pnmScv.indexInElement()); |
| 78 | 23 | const auto otherPNMScvDofIdx = otherPNMScv.dofIndex(); | |
| 79 | |||
| 80 | // check for intersections, skip if no intersection was found | ||
| 81 |
1/2✓ Branch 1 taken 23 times.
✗ Branch 2 not taken.
|
23 | const auto directlyCoupledFreeFlowElements = intersectingEntities(pnmPos, ffMomentumGridGeometry.boundingBoxTree()); |
| 82 |
2/2✓ Branch 0 taken 15 times.
✓ Branch 1 taken 8 times.
|
23 | if (directlyCoupledFreeFlowElements.empty()) |
| 83 | continue; | ||
| 84 | else | ||
| 85 |
1/2✓ Branch 1 taken 8 times.
✗ Branch 2 not taken.
|
8 | isCoupledPNMDof_[pnmDofIdx] = true; |
| 86 | |||
| 87 | // determine the normal direction of the local coupling interface heuristically: | ||
| 88 | // find all element intersections touching the pore and take the normal which | ||
| 89 | // occurs most frequently | ||
| 90 | 8 | const std::size_t couplingNormalDirectionIndex = [&] | |
| 91 | { | ||
| 92 | using Key = std::pair<std::size_t, bool>; | ||
| 93 | 8 | std::map<Key, std::size_t> result; | |
| 94 |
2/2✓ Branch 0 taken 8 times.
✓ Branch 1 taken 8 times.
|
16 | for (const auto eIdx : directlyCoupledFreeFlowElements) |
| 95 | { | ||
| 96 |
3/4✓ Branch 1 taken 8 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 8 times.
✓ Branch 5 taken 32 times.
|
56 | for (const auto& intersection : intersections(ffMomentumGridGeometry.gridView(), ffMomentumGridGeometry.element(eIdx))) |
| 97 | { | ||
| 98 |
2/2✓ Branch 1 taken 8 times.
✓ Branch 2 taken 24 times.
|
32 | if (intersectsPointGeometry(pnmPos, intersection.geometry())) |
| 99 | { | ||
| 100 | 8 | const auto& normal = intersection.centerUnitOuterNormal(); | |
| 101 | 8 | const auto normalAxis = Dumux::normalAxis(normal); | |
| 102 |
1/2✓ Branch 1 taken 8 times.
✗ Branch 2 not taken.
|
8 | const Key key = std::make_pair(normalAxis, std::signbit(normal[normalAxis])); |
| 103 |
1/2✓ Branch 1 taken 8 times.
✗ Branch 2 not taken.
|
8 | ++result[key]; |
| 104 | } | ||
| 105 | } | ||
| 106 | } | ||
| 107 | |||
| 108 | // TODO how to properly handle this corner (literally) case | ||
| 109 |
2/4✓ Branch 0 taken 8 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 8 times.
✗ Branch 3 not taken.
|
8 | if (directlyCoupledFreeFlowElements.size() == 1 && result.size() > 1) |
| 110 |
0/20✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
✗ Branch 14 not taken.
✗ Branch 15 not taken.
✗ Branch 17 not taken.
✗ Branch 18 not taken.
✗ Branch 20 not taken.
✗ Branch 21 not taken.
✗ Branch 23 not taken.
✗ Branch 24 not taken.
✗ Branch 26 not taken.
✗ Branch 27 not taken.
✗ Branch 29 not taken.
✗ Branch 30 not taken.
|
8 | DUNE_THROW(Dune::InvalidStateException, "Pore may not intersect with faces of different orientation when coupled to only one element"); |
| 111 | |||
| 112 |
0/2✗ Branch 0 not taken.
✗ Branch 1 not taken.
|
8 | return std::max_element(result.begin(), result.end(), [](const auto& x, const auto& y) { return x.second < y.second;})->first.first; |
| 113 |
1/2✓ Branch 1 taken 8 times.
✗ Branch 2 not taken.
|
16 | }(); |
| 114 | |||
| 115 | using Scalar = typename FreeFlowMomentumGridGeometry::GridView::ctype; | ||
| 116 | 8 | const Scalar couplingPoreRadius = pnmGridGeometry.poreInscribedRadius(pnmDofIdx); | |
| 117 | 16 | GlobalPosition lowerLeft = pnmPos - GlobalPosition(couplingPoreRadius); | |
| 118 | 8 | lowerLeft[couplingNormalDirectionIndex] = pnmPos[couplingNormalDirectionIndex]; | |
| 119 |
1/2✓ Branch 1 taken 8 times.
✗ Branch 2 not taken.
|
16 | GlobalPosition upperRight = pnmPos + GlobalPosition(couplingPoreRadius); |
| 120 |
1/2✓ Branch 1 taken 8 times.
✗ Branch 2 not taken.
|
8 | upperRight[couplingNormalDirectionIndex] = pnmPos[couplingNormalDirectionIndex]; |
| 121 | |||
| 122 |
1/2✓ Branch 1 taken 8 times.
✗ Branch 2 not taken.
|
8 | auto axes = std::move(std::bitset<FreeFlowMomentumGridGeometry::Grid::dimensionworld>{}.set()); |
| 123 |
1/2✓ Branch 1 taken 8 times.
✗ Branch 2 not taken.
|
8 | axes.set(couplingNormalDirectionIndex, false); |
| 124 | |||
| 125 | using PoreIntersectionGeometryType = Dune::AxisAlignedCubeGeometry<Scalar, | ||
| 126 | FreeFlowMomentumGridGeometry::GridView::dimension-1, | ||
| 127 | FreeFlowMomentumGridGeometry::GridView::dimensionworld>; | ||
| 128 | |||
| 129 | 8 | PoreIntersectionGeometryType poreIntersectionGeometry(lowerLeft, upperRight, axes); | |
| 130 |
1/2✓ Branch 1 taken 8 times.
✗ Branch 2 not taken.
|
8 | const auto allCoupledFreeFlowElements = intersectingEntities(std::move(poreIntersectionGeometry), ffMomentumGridGeometry.boundingBoxTree()); |
| 131 | |||
| 132 | |||
| 133 |
3/4✓ Branch 1 taken 28 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 28 times.
✓ Branch 4 taken 8 times.
|
36 | for (const auto& ffElementInfo : allCoupledFreeFlowElements) |
| 134 | { | ||
| 135 | 28 | const auto freeFlowElementIndex = ffElementInfo.second(); | |
| 136 |
2/4✓ Branch 1 taken 28 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 28 times.
✗ Branch 5 not taken.
|
28 | pnmElementToFreeFlowElementsMap_[pnmElementIdx].push_back(freeFlowElementIndex); |
| 137 |
1/2✓ Branch 1 taken 28 times.
✗ Branch 2 not taken.
|
28 | freeFlowElementToPNMElementMap_[freeFlowElementIndex] = pnmElementIdx; |
| 138 | |||
| 139 |
2/4✓ Branch 1 taken 28 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 28 times.
✗ Branch 5 not taken.
|
28 | pnmToFreeFlowMassStencils_[pnmElementIdx].push_back(freeFlowElementIndex); |
| 140 |
2/4✓ Branch 1 taken 28 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 28 times.
✗ Branch 5 not taken.
|
28 | freeFlowMassToPNMStencils_[freeFlowElementIndex].push_back(pnmDofIdx); |
| 141 | |||
| 142 |
1/2✓ Branch 1 taken 28 times.
✗ Branch 2 not taken.
|
56 | ffFvGeometry.bindElement(ffMomentumGridGeometry.element(freeFlowElementIndex)); |
| 143 |
1/2✓ Branch 1 taken 28 times.
✗ Branch 2 not taken.
|
28 | const auto coupledFreeFlowMomentumDofIndices = coupledFFMomentumDofs_(ffFvGeometry, ffMassGridGeometry, pnmPos, couplingPoreRadius, couplingNormalDirectionIndex); |
| 144 | |||
| 145 |
2/4✓ Branch 1 taken 28 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 28 times.
✗ Branch 5 not taken.
|
28 | pnmToFreeFlowMomentumStencils_[pnmElementIdx].push_back(coupledFreeFlowMomentumDofIndices.coupledFrontalDof); |
| 146 |
2/4✓ Branch 1 taken 28 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 28 times.
✗ Branch 5 not taken.
|
28 | freeFlowMomentumToPNMStencils_[coupledFreeFlowMomentumDofIndices.coupledFrontalDof].push_back(pnmDofIdx); |
| 147 |
2/4✓ Branch 1 taken 28 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 28 times.
✗ Branch 5 not taken.
|
28 | freeFlowMomentumToPNMStencils_[coupledFreeFlowMomentumDofIndices.coupledFrontalDof].push_back(otherPNMScvDofIdx); |
| 148 | |||
| 149 |
2/4✗ Branch 0 not taken.
✓ Branch 1 taken 28 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 28 times.
|
28 | isCoupledFreeFlowMomentumDof_[coupledFreeFlowMomentumDofIndices.coupledFrontalDof] = true; |
| 150 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 28 times.
|
28 | isCoupledFreeFlowMomentumDofOnInterface_[coupledFreeFlowMomentumDofIndices.coupledFrontalDof] = true; |
| 151 | |||
| 152 | // treat the coupled ff dofs not directly on interface | ||
| 153 |
2/2✓ Branch 0 taken 56 times.
✓ Branch 1 taken 28 times.
|
84 | for (const auto ffDofIdx : coupledFreeFlowMomentumDofIndices.coupledLateralDofs) |
| 154 | { | ||
| 155 |
2/4✓ Branch 1 taken 56 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 56 times.
✗ Branch 5 not taken.
|
56 | freeFlowMomentumToPNMStencils_[ffDofIdx].push_back(pnmDofIdx); |
| 156 |
2/4✓ Branch 1 taken 56 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 56 times.
✗ Branch 5 not taken.
|
56 | freeFlowMomentumToPNMStencils_[ffDofIdx].push_back(otherPNMScvDofIdx); |
| 157 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 56 times.
|
56 | isCoupledFreeFlowMomentumDof_[ffDofIdx] = true; |
| 158 | } | ||
| 159 | } | ||
| 160 | } | ||
| 161 | } | ||
| 162 | |||
| 163 | 6 | std::cout << "took " << watch.elapsed() << " seconds." << std::endl; | |
| 164 | 6 | } | |
| 165 | |||
| 166 | /*! | ||
| 167 | * \brief returns an iterable container of all indices of degrees of freedom of domain j | ||
| 168 | * that couple with / influence the element residual of the given element of domain i | ||
| 169 | * \param eIdxI the index of the coupled element of domain í | ||
| 170 | */ | ||
| 171 | 218 | const std::vector<std::size_t>& poreNetworkToFreeFlowMomentumCouplingStencil(const std::size_t eIdxI) const | |
| 172 | { | ||
| 173 | 436 | if (isCoupledPoreNetworkElement(eIdxI)) | |
| 174 | 190 | return pnmToFreeFlowMomentumStencils_.at(eIdxI); | |
| 175 | else | ||
| 176 | 28 | return emptyStencil_; | |
| 177 | } | ||
| 178 | |||
| 179 | /*! | ||
| 180 | * \brief returns an iterable container of all indices of degrees of freedom of domain j | ||
| 181 | * that couple with / influence the element residual of the given element of domain i | ||
| 182 | * \param eIdxI the index of the coupled element of domain í | ||
| 183 | */ | ||
| 184 | 221 | const std::vector<std::size_t>& poreNetworkToFreeFlowMassCouplingStencil(const std::size_t eIdxI) const | |
| 185 | { | ||
| 186 | 442 | if (isCoupledPoreNetworkElement(eIdxI)) | |
| 187 | 193 | return pnmToFreeFlowMassStencils_.at(eIdxI); | |
| 188 | else | ||
| 189 | 28 | return emptyStencil_; | |
| 190 | } | ||
| 191 | |||
| 192 | /*! | ||
| 193 | * \brief returns an iterable container of all indices of degrees of freedom of domain j | ||
| 194 | * that couple with / influence the element residual of the given element of domain i | ||
| 195 | * \param eIdxI the index of the coupled element of domain i | ||
| 196 | */ | ||
| 197 | 75842 | const std::vector<std::size_t>& freeFlowMassToPoreNetworkCouplingStencil(const std::size_t eIdxI) const | |
| 198 | { | ||
| 199 | 151684 | if (isCoupledFreeFlowElement(eIdxI)) | |
| 200 | 584 | return freeFlowMassToPNMStencils_.at(eIdxI); | |
| 201 | else | ||
| 202 | 75258 | return emptyStencil_; | |
| 203 | } | ||
| 204 | |||
| 205 | /*! | ||
| 206 | * \brief returns an iterable container of all indices of degrees of freedom of domain j | ||
| 207 | * that couple with / influence the element residual of the given element of domain i | ||
| 208 | * \param dofIndex the degree of freedom index | ||
| 209 | */ | ||
| 210 | 303368 | const std::vector<std::size_t>& freeFlowMomentumToPoreNetworkCouplingStencil(const std::size_t dofIndex) const | |
| 211 | { | ||
| 212 | 606736 | if (isCoupledFreeFlowMomentumDof(dofIndex)) | |
| 213 | 2104 | return freeFlowMomentumToPNMStencils_.at(dofIndex); | |
| 214 | else | ||
| 215 | 301264 | return emptyStencil_; | |
| 216 | } | ||
| 217 | |||
| 218 | /*! | ||
| 219 | * \brief Return if an element residual with index eIdx of domain i is coupled to domain j | ||
| 220 | */ | ||
| 221 | 4843239 | bool isCoupledFreeFlowElement(std::size_t eIdx) const | |
| 222 | { | ||
| 223 |
5/6✗ Branch 2 not taken.
✓ Branch 3 taken 72998 times.
✓ Branch 5 taken 21160 times.
✓ Branch 6 taken 2960 times.
✓ Branch 8 taken 72998 times.
✓ Branch 9 taken 565 times.
|
5090318 | return static_cast<bool>(freeFlowElementToPNMElementMap_.count(eIdx)); |
| 224 | } | ||
| 225 | |||
| 226 | /*! | ||
| 227 | * \brief Return if an element residual with index eIdx of domain i is coupled to domain j | ||
| 228 | */ | ||
| 229 | 595544 | bool isCoupledFreeFlowMomentumDof(std::size_t dofIdx) const | |
| 230 | { | ||
| 231 |
7/8✓ Branch 0 taken 138 times.
✓ Branch 1 taken 46 times.
✓ Branch 2 taken 336 times.
✓ Branch 3 taken 291656 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 303368 times.
✓ Branch 6 taken 2104 times.
✓ Branch 7 taken 301264 times.
|
595544 | return isCoupledFreeFlowMomentumDof_[dofIdx]; |
| 232 | } | ||
| 233 | |||
| 234 | /*! | ||
| 235 | * \brief Return if an element residual with index eIdx of domain i is coupled to domain j | ||
| 236 | */ | ||
| 237 | 3188 | bool isCoupledPoreNetworkElement(std::size_t eIdx) const | |
| 238 | { | ||
| 239 |
7/8✓ Branch 1 taken 21 times.
✓ Branch 2 taken 14 times.
✓ Branch 4 taken 52 times.
✓ Branch 5 taken 88 times.
✓ Branch 7 taken 1522 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 1039 times.
✓ Branch 11 taken 13 times.
|
6376 | return static_cast<bool>(pnmElementToFreeFlowElementsMap_.count(eIdx)); |
| 240 | } | ||
| 241 | |||
| 242 | /*! | ||
| 243 | * \brief Return if an element residual with index eIdx of domain i is coupled to domain j | ||
| 244 | */ | ||
| 245 | 11006 | bool isCoupledPoreNetworkDof(std::size_t dofIdx) const | |
| 246 | { | ||
| 247 |
8/8✓ Branch 0 taken 892 times.
✓ Branch 1 taken 892 times.
✓ Branch 2 taken 3611 times.
✓ Branch 3 taken 3611 times.
✓ Branch 4 taken 745 times.
✓ Branch 5 taken 881 times.
✓ Branch 6 taken 201 times.
✓ Branch 7 taken 173 times.
|
11006 | return isCoupledPNMDof_[dofIdx]; |
| 248 | } | ||
| 249 | |||
| 250 | 918946 | bool isCoupledFreeFlowMomentumScvf(std::size_t scvfIdx) const | |
| 251 | { | ||
| 252 | 1837892 | return isCoupledFrontalFreeFlowMomentumScvf_.count(scvfIdx); | |
| 253 | } | ||
| 254 | |||
| 255 | 800964 | bool isCoupledFreeFlowMomentumLateralScvf(std::size_t scvfIdx) const | |
| 256 | { | ||
| 257 | 1601928 | return isCoupledLateralFreeFlowMomentumScvf_.count(scvfIdx); | |
| 258 | } | ||
| 259 | |||
| 260 | 343482 | bool isCoupledFreeFlowMassScvf(std::size_t scvfIdx) const | |
| 261 | { | ||
| 262 |
5/6✓ Branch 1 taken 55 times.
✗ Branch 2 not taken.
✓ Branch 5 taken 89795 times.
✓ Branch 6 taken 39016 times.
✓ Branch 8 taken 118930 times.
✓ Branch 9 taken 77650 times.
|
686964 | return isCoupledFreeFlowMassScvf_.count(scvfIdx); |
| 263 | } | ||
| 264 | |||
| 265 | const auto& pnmElementToFreeFlowElementsMap() const | ||
| 266 | 31 | { return pnmElementToFreeFlowElementsMap_;} | |
| 267 | |||
| 268 | const auto& freeFlowElementToPNMElementMap() const | ||
| 269 |
1/2✓ Branch 1 taken 892 times.
✗ Branch 2 not taken.
|
1457 | { return freeFlowElementToPNMElementMap_; } |
| 270 | |||
| 271 | private: | ||
| 272 | |||
| 273 | 6 | void clear_() | |
| 274 | { | ||
| 275 | 6 | pnmElementToFreeFlowElementsMap_.clear(); | |
| 276 | 6 | freeFlowElementToPNMElementMap_.clear(); | |
| 277 | 6 | isCoupledPNMDof_.clear(); | |
| 278 | 6 | isCoupledFreeFlowMomentumDof_.clear(); | |
| 279 | 6 | isCoupledFreeFlowMomentumDofOnInterface_.clear(); | |
| 280 | 6 | pnmToFreeFlowMassStencils_.clear(); | |
| 281 | 6 | pnmToFreeFlowMomentumStencils_.clear(); | |
| 282 | 6 | freeFlowMassToPNMStencils_.clear(); | |
| 283 | 6 | freeFlowMomentumToPNMStencils_.clear(); | |
| 284 | 6 | } | |
| 285 | |||
| 286 | template<class FreeFlowMomentumGridGeometry, class PoreNetworkGridGeometry> | ||
| 287 | 6 | void resize_(const FreeFlowMomentumGridGeometry& ffMomentumGridGeometry, | |
| 288 | const PoreNetworkGridGeometry& pnmGridGeometry) | ||
| 289 | { | ||
| 290 | 6 | const auto numPNMDofs = pnmGridGeometry.numDofs(); | |
| 291 | 6 | const auto numFreeFlowMomentumDofs = ffMomentumGridGeometry.numDofs(); | |
| 292 | 6 | isCoupledPNMDof_.resize(numPNMDofs, false); | |
| 293 | 6 | isCoupledFreeFlowMomentumDof_.resize(numFreeFlowMomentumDofs, false); | |
| 294 | 6 | isCoupledFreeFlowMomentumDofOnInterface_.resize(numFreeFlowMomentumDofs, false); | |
| 295 | 6 | } | |
| 296 | |||
| 297 | //! get the free-flow momentum dofs that are coupled to the PNM dofs | ||
| 298 | template<class FVElementGeometry, class FreeFlowMassGridGeometry, class GlobalPosition, class Scalar> | ||
| 299 | 28 | auto coupledFFMomentumDofs_(const FVElementGeometry& fvGeometry, | |
| 300 | const FreeFlowMassGridGeometry& ffMassGridGeometry, | ||
| 301 | const GlobalPosition& pnmPos, | ||
| 302 | const Scalar couplingPoreRadius, | ||
| 303 | const int couplingInterfaceDirectionIdx) | ||
| 304 | { | ||
| 305 | |||
| 306 | 28 | struct Result | |
| 307 | { | ||
| 308 | Dune::ReservedVector<std::size_t, FVElementGeometry::maxNumElementScvs> coupledLateralDofs; | ||
| 309 | std::size_t coupledFrontalDof; | ||
| 310 | 28 | } result; | |
| 311 | |||
| 312 | |||
| 313 | using std::abs; | ||
| 314 |
2/2✓ Branch 2 taken 112 times.
✓ Branch 3 taken 28 times.
|
140 | for (const auto& scv : scvs(fvGeometry)) |
| 315 | { | ||
| 316 | 112 | const Scalar eps = diameter(fvGeometry.geometry(scv))*1e-6; // TODO | |
| 317 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 112 times.
|
112 | assert(eps < couplingPoreRadius); |
| 318 | |||
| 319 |
2/2✓ Branch 0 taken 56 times.
✓ Branch 1 taken 56 times.
|
112 | if (scv.dofAxis() == couplingInterfaceDirectionIdx) // the free flow dofs that lie within the coupling interface |
| 320 | { | ||
| 321 |
2/2✓ Branch 0 taken 28 times.
✓ Branch 1 taken 28 times.
|
56 | if (abs(scv.dofPosition()[couplingInterfaceDirectionIdx] - pnmPos[couplingInterfaceDirectionIdx]) < eps) |
| 322 | { | ||
| 323 | 28 | result.coupledFrontalDof = scv.dofIndex(); | |
| 324 | |||
| 325 | // treat scvfs | ||
| 326 |
4/4✓ Branch 1 taken 56 times.
✓ Branch 2 taken 56 times.
✓ Branch 4 taken 112 times.
✓ Branch 5 taken 28 times.
|
364 | for (const auto& scvf : scvfs(fvGeometry, scv)) |
| 327 | { | ||
| 328 | // add lateral faces "standing" on coupling interface | ||
| 329 |
3/4✓ Branch 0 taken 56 times.
✓ Branch 1 taken 56 times.
✓ Branch 2 taken 56 times.
✗ Branch 3 not taken.
|
112 | if (scvf.isLateral() && !scvf.boundary()) |
| 330 | 56 | isCoupledLateralFreeFlowMomentumScvf_[scvf.index()] = true; | |
| 331 |
3/4✓ Branch 0 taken 56 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 28 times.
✓ Branch 3 taken 28 times.
|
56 | else if (scvf.isFrontal() && scvf.boundary()) // add face lying on interface |
| 332 | { | ||
| 333 | 28 | isCoupledFrontalFreeFlowMomentumScvf_[scvf.index()] = true; | |
| 334 | |||
| 335 | 28 | const auto& element = ffMassGridGeometry.element(fvGeometry.elementIndex()); // this local variable is needed to prevent a memory error | |
| 336 | 56 | const auto ffMassFVGeometry = localView(ffMassGridGeometry).bindElement(element); | |
| 337 |
4/4✓ Branch 0 taken 28 times.
✓ Branch 1 taken 84 times.
✓ Branch 2 taken 112 times.
✓ Branch 3 taken 28 times.
|
140 | for (const auto& ffMassScvf : scvfs(ffMassFVGeometry)) |
| 338 | { | ||
| 339 |
2/2✓ Branch 0 taken 28 times.
✓ Branch 1 taken 84 times.
|
112 | if (abs(ffMassScvf.center()[couplingInterfaceDirectionIdx] - pnmPos[couplingInterfaceDirectionIdx]) < eps) |
| 340 | 28 | isCoupledFreeFlowMassScvf_[ffMassScvf.index()] = true; | |
| 341 | } | ||
| 342 | } | ||
| 343 | } | ||
| 344 | } | ||
| 345 | } | ||
| 346 | else // the free flow dofs perpendicular to the coupling interface | ||
| 347 | { | ||
| 348 | bool isCoupledDof = false; | ||
| 349 | |||
| 350 |
2/2✓ Branch 0 taken 112 times.
✓ Branch 1 taken 56 times.
|
168 | for (int dimIdx = 0; dimIdx < GlobalPosition::dimension; ++dimIdx) |
| 351 | { | ||
| 352 |
3/4✓ Branch 0 taken 56 times.
✓ Branch 1 taken 56 times.
✓ Branch 2 taken 56 times.
✗ Branch 3 not taken.
|
112 | if (dimIdx == couplingInterfaceDirectionIdx || scv.boundary()) |
| 353 | 56 | continue; | |
| 354 | |||
| 355 | 56 | isCoupledDof = abs(scv.dofPosition()[dimIdx] - pnmPos[dimIdx]) < couplingPoreRadius + eps; | |
| 356 | } | ||
| 357 | |||
| 358 |
1/2✓ Branch 0 taken 56 times.
✗ Branch 1 not taken.
|
56 | if (isCoupledDof) |
| 359 | { | ||
| 360 | 56 | result.coupledLateralDofs.push_back(scv.dofIndex()); | |
| 361 | |||
| 362 | // treat scvfs | ||
| 363 |
4/4✓ Branch 1 taken 112 times.
✓ Branch 2 taken 56 times.
✓ Branch 4 taken 168 times.
✓ Branch 5 taken 56 times.
|
504 | for (const auto& scvf : scvfs(fvGeometry, scv)) |
| 364 | { | ||
| 365 |
4/4✓ Branch 0 taken 112 times.
✓ Branch 1 taken 56 times.
✓ Branch 2 taken 56 times.
✓ Branch 3 taken 56 times.
|
168 | if (scvf.isLateral() && scvf.boundary()) |
| 366 | { | ||
| 367 | // add lateral scvfs lying on interface | ||
| 368 |
1/2✓ Branch 0 taken 56 times.
✗ Branch 1 not taken.
|
56 | if (abs(scvf.ipGlobal()[couplingInterfaceDirectionIdx] - pnmPos[couplingInterfaceDirectionIdx]) < eps) |
| 369 | 56 | isCoupledLateralFreeFlowMomentumScvf_[scvf.index()] = true; | |
| 370 | } | ||
| 371 | } | ||
| 372 | } | ||
| 373 | } | ||
| 374 | } | ||
| 375 | |||
| 376 | 28 | return result; | |
| 377 | } | ||
| 378 | |||
| 379 | std::vector<std::size_t> emptyStencil_; | ||
| 380 | |||
| 381 | std::unordered_map<std::size_t, std::vector<std::size_t>> pnmElementToFreeFlowElementsMap_; | ||
| 382 | std::unordered_map<std::size_t, std::size_t> freeFlowElementToPNMElementMap_; | ||
| 383 | |||
| 384 | std::vector<bool> isCoupledPNMDof_; | ||
| 385 | std::vector<bool> isCoupledFreeFlowMomentumDof_; | ||
| 386 | std::vector<bool> isCoupledFreeFlowMomentumDofOnInterface_; | ||
| 387 | |||
| 388 | std::unordered_map<std::size_t, bool> isCoupledLateralFreeFlowMomentumScvf_; | ||
| 389 | std::unordered_map<std::size_t, bool> isCoupledFrontalFreeFlowMomentumScvf_; | ||
| 390 | std::unordered_map<std::size_t, bool> isCoupledFreeFlowMassScvf_; | ||
| 391 | |||
| 392 | MapType pnmToFreeFlowMassStencils_; | ||
| 393 | MapType pnmToFreeFlowMomentumStencils_; | ||
| 394 | MapType freeFlowMassToPNMStencils_; | ||
| 395 | MapType freeFlowMomentumToPNMStencils_; | ||
| 396 | }; | ||
| 397 | |||
| 398 | } // end namespace Dumux | ||
| 399 | |||
| 400 | #endif | ||
| 401 |