| 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 NavierStokesModel | ||
| 10 | * \copydoc Dumux::FluxOverAxisAlignedSurface | ||
| 11 | */ | ||
| 12 | #ifndef DUMUX_FREELOW_NAVIERSTOKES_FLUX_OVER_AXISALIGNED_SURFACE_HH | ||
| 13 | #define DUMUX_FREELOW_NAVIERSTOKES_FLUX_OVER_AXISALIGNED_SURFACE_HH | ||
| 14 | |||
| 15 | #include <algorithm> | ||
| 16 | #include <type_traits> | ||
| 17 | #include <vector> | ||
| 18 | |||
| 19 | #include <dune/common/exceptions.hh> | ||
| 20 | #include <dune/geometry/axisalignedcubegeometry.hh> | ||
| 21 | |||
| 22 | #include <dumux/common/parameters.hh> | ||
| 23 | #include <dumux/geometry/diameter.hh> | ||
| 24 | #include <dumux/geometry/distance.hh> | ||
| 25 | #include <dumux/geometry/intersectspointgeometry.hh> | ||
| 26 | #include <dumux/geometry/geometricentityset.hh> | ||
| 27 | #include <dumux/geometry/intersectingentities.hh> | ||
| 28 | |||
| 29 | namespace Dumux { | ||
| 30 | |||
| 31 | /*! | ||
| 32 | * \ingroup NavierStokesModel | ||
| 33 | * \brief Class used to calculate fluxes over axis-aligned surfaces. | ||
| 34 | */ | ||
| 35 | template<class GridVariables, class SolutionVector, class LocalResidual> | ||
| 36 |
1/2✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
|
9 | class FluxOverAxisAlignedSurface |
| 37 | { | ||
| 38 | using Scalar = typename GridVariables::Scalar; | ||
| 39 | using GridGeometry = typename GridVariables::GridGeometry; | ||
| 40 | using FVElementGeometry = typename GridGeometry::LocalView; | ||
| 41 | using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace; | ||
| 42 | using GridView = typename GridGeometry::GridView; | ||
| 43 | using VolumeVariables = typename GridVariables::VolumeVariables; | ||
| 44 | using Element = typename GridView::template Codim<0>::Entity; | ||
| 45 | using NumEqVector = typename LocalResidual::ElementResidualVector::value_type; | ||
| 46 | |||
| 47 | static constexpr auto dim = GridView::dimension; | ||
| 48 | static constexpr auto dimWorld = GridView::dimensionworld; | ||
| 49 | |||
| 50 | static_assert(dim > 1, "Only implemented for dim > 1"); | ||
| 51 | |||
| 52 | using GlobalPosition = typename Element::Geometry::GlobalCoordinate; | ||
| 53 | |||
| 54 | // in 2D, the surface is represented as a line | ||
| 55 | using SurfaceT = Dune::AxisAlignedCubeGeometry<Scalar, (dim == 2 ? 1 : 2), dimWorld>; | ||
| 56 | |||
| 57 | struct SurfaceData | ||
| 58 | { | ||
| 59 | SurfaceT surface; | ||
| 60 | std::size_t normalDirectionIndex; | ||
| 61 | NumEqVector flux; | ||
| 62 | }; | ||
| 63 | |||
| 64 | public: | ||
| 65 | |||
| 66 | using Surface = SurfaceT; | ||
| 67 | |||
| 68 | /*! | ||
| 69 | * \brief The constructor | ||
| 70 | */ | ||
| 71 | 9 | FluxOverAxisAlignedSurface(const GridVariables& gridVariables, | |
| 72 | const SolutionVector& sol, | ||
| 73 | const LocalResidual& localResidual, | ||
| 74 | bool nonIntersectingSurfaceIsError = false) | ||
| 75 | : gridVariables_(gridVariables) | ||
| 76 | , sol_(sol) | ||
| 77 | , localResidual_(localResidual) | ||
| 78 |
1/2✓ Branch 1 taken 9 times.
✗ Branch 2 not taken.
|
9 | , nonIntersectingSurfaceIsError_(nonIntersectingSurfaceIsError) |
| 79 | { | ||
| 80 |
1/2✓ Branch 1 taken 9 times.
✗ Branch 2 not taken.
|
9 | verbose_ = getParamFromGroup<bool>(problem_().paramGroup(), "FluxOverAxisAlignedSurface.Verbose", false); |
| 81 | 9 | } | |
| 82 | |||
| 83 | /*! | ||
| 84 | * \brief Add an axis-aligned surface with a given name | ||
| 85 | * | ||
| 86 | * \param name The name of the surface | ||
| 87 | * \param surface The surface to add | ||
| 88 | */ | ||
| 89 | template<class T> | ||
| 90 | void addAxisAlignedSurface(const std::string& name, T&& surface) | ||
| 91 | { | ||
| 92 | static_assert(std::is_same_v<std::decay_t<T>, Surface>); | ||
| 93 | surfaces_.emplace(std::make_pair( | ||
| 94 | name, std::make_pair(surface, NumEqVector(0.0)) | ||
| 95 | )); | ||
| 96 | } | ||
| 97 | |||
| 98 | /*! | ||
| 99 | * \brief Add an axis-aligned surface (segment in 2D) with a given name, specifying the surface's corner points. | ||
| 100 | * | ||
| 101 | * \param name The name of the surface | ||
| 102 | * \param lowerLeft Lower left corner of surface | ||
| 103 | * \param upperRight Upper right corner of surface | ||
| 104 | */ | ||
| 105 | 19 | void addAxisAlignedSurface(const std::string& name, | |
| 106 | const GlobalPosition& lowerLeft, | ||
| 107 | const GlobalPosition& upperRight) | ||
| 108 | { | ||
| 109 | using std::abs; | ||
| 110 | 19 | const GlobalPosition v = upperRight - lowerLeft; | |
| 111 |
8/28✗ 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 taken 9 times.
✗ Branch 17 not taken.
✓ Branch 18 taken 9 times.
✗ Branch 19 not taken.
✓ Branch 20 taken 3 times.
✓ Branch 21 taken 7 times.
✓ Branch 22 taken 3 times.
✓ Branch 23 taken 7 times.
✓ Branch 24 taken 7 times.
✗ Branch 25 not taken.
✓ Branch 26 taken 7 times.
✗ Branch 27 not taken.
|
109 | const auto it = std::find_if(v.begin(), v.end(), [](const auto& x){ return abs(x) < 1e-20; }); |
| 112 |
2/4✗ Branch 0 not taken.
✓ Branch 1 taken 19 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 19 times.
|
38 | if (it == v.end()) |
| 113 | ✗ | DUNE_THROW(Dune::InvalidStateException, "Surface is not axis-parallel!"); | |
| 114 | |||
| 115 |
2/4✗ Branch 0 not taken.
✓ Branch 1 taken 19 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 19 times.
|
38 | const std::size_t normalDirectionIndex = std::distance(v.begin(), it); |
| 116 | 38 | auto inSurfaceAxes = std::move(std::bitset<dimWorld>{}.set()); | |
| 117 | 19 | inSurfaceAxes.set(normalDirectionIndex, false); | |
| 118 | 19 | auto surface = Surface(lowerLeft, upperRight, inSurfaceAxes); | |
| 119 | |||
| 120 |
3/8✓ Branch 1 taken 19 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 19 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 19 times.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
|
19 | surfaces_.emplace(std::make_pair( |
| 121 | name, | ||
| 122 | SurfaceData{ | ||
| 123 |
1/2✓ Branch 1 taken 19 times.
✗ Branch 2 not taken.
|
19 | std::move(surface), normalDirectionIndex, NumEqVector(0.0) |
| 124 | } | ||
| 125 | )); | ||
| 126 | 19 | } | |
| 127 | |||
| 128 | /*! | ||
| 129 | * \brief Add an axis-aligned plane (line in 2D) with a given name, specifying the planes's center and normal. | ||
| 130 | * | ||
| 131 | * \param name The name of the plane | ||
| 132 | * \param center Center point of the plane | ||
| 133 | * \param normalDirectionIndex Index of the plane's normal axis (0=x, 1=y, 2=z) | ||
| 134 | */ | ||
| 135 | 6 | void addAxisAlignedPlane(const std::string& name, | |
| 136 | const GlobalPosition& center, | ||
| 137 | const std::size_t normalDirectionIndex) | ||
| 138 | { | ||
| 139 | 12 | GlobalPosition lowerLeft = gridVariables_.gridGeometry().bBoxMin(); | |
| 140 | 12 | GlobalPosition upperRight = gridVariables_.gridGeometry().bBoxMax(); | |
| 141 | |||
| 142 | 12 | lowerLeft[normalDirectionIndex] = center[normalDirectionIndex]; | |
| 143 | 12 | upperRight[normalDirectionIndex] = center[normalDirectionIndex]; | |
| 144 | |||
| 145 | 12 | auto inSurfaceAxes = std::move(std::bitset<dimWorld>{}.set()); | |
| 146 | 6 | inSurfaceAxes.set(normalDirectionIndex, false); | |
| 147 | 6 | auto surface = Surface(lowerLeft, upperRight, inSurfaceAxes); | |
| 148 | |||
| 149 |
3/8✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 6 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 6 times.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
|
6 | surfaces_.emplace(std::make_pair( |
| 150 | name, | ||
| 151 | SurfaceData{ | ||
| 152 |
1/2✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
|
6 | std::move(surface), normalDirectionIndex, NumEqVector(0.0) |
| 153 | } | ||
| 154 | )); | ||
| 155 | 6 | } | |
| 156 | |||
| 157 | /*! | ||
| 158 | * \brief Calculate the fluxes over all surfaces. | ||
| 159 | */ | ||
| 160 | void calculateAllFluxes() | ||
| 161 | { | ||
| 162 | 9 | auto fluxType = [this](const auto& element, | |
| 163 | const auto& fvGeometry, | ||
| 164 | const auto& elemVolVars, | ||
| 165 | const auto& scvf, | ||
| 166 | ✗ | const auto& elemFluxVarsCache) | |
| 167 | { | ||
| 168 | return localResidual_.evalFlux( | ||
| 169 | problem_(), element, fvGeometry, elemVolVars, elemFluxVarsCache, scvf | ||
| 170 | 2428 | ); | |
| 171 | }; | ||
| 172 | |||
| 173 |
2/4✓ Branch 1 taken 8 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
|
9 | calculateFluxes(fluxType); |
| 174 | } | ||
| 175 | |||
| 176 | /*! | ||
| 177 | * \brief Calculate the fluxes over all surfaces for a given flux type. | ||
| 178 | * | ||
| 179 | * \param fluxType The flux type. This can be a lambda of the following form: | ||
| 180 | * [](const auto& element, | ||
| 181 | const auto& fvGeometry, | ||
| 182 | const auto& elemVolVars, | ||
| 183 | const auto& scvf, | ||
| 184 | const auto& elemFluxVarsCache) | ||
| 185 | { return ... ; } | ||
| 186 | */ | ||
| 187 | template<class FluxType> | ||
| 188 | 9 | void calculateFluxes(const FluxType& fluxType) | |
| 189 | { | ||
| 190 | // make sure to reset all the values of the surfaces, in case this method has been called already before | ||
| 191 |
2/2✓ Branch 0 taken 25 times.
✓ Branch 1 taken 9 times.
|
52 | for (auto& surface : surfaces_) |
| 192 | 50 | surface.second.flux = 0.0; | |
| 193 | |||
| 194 | 9 | snapSurfaceToClosestFace_(); | |
| 195 | 9 | calculateFluxes_(fluxType); | |
| 196 | 9 | } | |
| 197 | |||
| 198 | /*! | ||
| 199 | * \brief Return the flux over given surface | ||
| 200 | * | ||
| 201 | * \param name The name of the surface | ||
| 202 | */ | ||
| 203 | const auto& flux(const std::string& name) const | ||
| 204 | { | ||
| 205 |
9/15✓ Branch 1 taken 8 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 8 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 7 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 1 times.
✓ Branch 10 taken 6 times.
✗ Branch 11 not taken.
✓ Branch 12 taken 3 times.
✓ Branch 13 taken 4 times.
✗ Branch 14 not taken.
✓ Branch 15 taken 3 times.
✓ Branch 16 taken 4 times.
✗ Branch 17 not taken.
|
19 | return surfaces_.at(name).flux; |
| 206 | } | ||
| 207 | |||
| 208 | /*! | ||
| 209 | * \brief Provides access to all surfaces. | ||
| 210 | */ | ||
| 211 | const std::map<std::string, SurfaceData>& surfaces() const | ||
| 212 | { return surfaces_; } | ||
| 213 | |||
| 214 | /*! | ||
| 215 | * \brief Prints all fluxes. | ||
| 216 | */ | ||
| 217 | 4 | void printAllFluxes() const | |
| 218 | { | ||
| 219 |
2/2✓ Branch 0 taken 12 times.
✓ Branch 1 taken 4 times.
|
24 | for (const auto& [name, data] : surfaces_) |
| 220 | 36 | std::cout << "Flux over surface " << name << ": " << data.flux << std::endl; | |
| 221 | 4 | } | |
| 222 | |||
| 223 | /*! | ||
| 224 | * \brief Set if non-intersecting surfaces are treated as error. | ||
| 225 | */ | ||
| 226 | void setNonIntersectingSurfaceIsError(bool isError = true) | ||
| 227 | { nonIntersectingSurfaceIsError_ = isError; } | ||
| 228 | |||
| 229 | private: | ||
| 230 | |||
| 231 | template<class FluxType> | ||
| 232 | 9 | void calculateFluxes_(const FluxType& fluxType) | |
| 233 | { | ||
| 234 |
1/2✓ Branch 1 taken 9 times.
✗ Branch 2 not taken.
|
9 | auto fvGeometry = localView(problem_().gridGeometry()); |
| 235 |
2/4✓ Branch 1 taken 9 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 9 times.
✗ Branch 5 not taken.
|
18 | auto elemVolVars = localView(gridVariables_.curGridVolVars()); |
| 236 |
2/4✓ Branch 1 taken 9 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 9 times.
✗ Branch 5 not taken.
|
18 | auto elemFluxVarsCache = localView(gridVariables_.gridFluxVarsCache()); |
| 237 | |||
| 238 |
4/8✓ Branch 1 taken 9 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 9 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 22409 times.
✗ Branch 7 not taken.
✓ Branch 11 taken 992 times.
✗ Branch 12 not taken.
|
44827 | for (const auto& element : elements(problem_().gridGeometry().gridView())) |
| 239 | { | ||
| 240 |
1/2✓ Branch 1 taken 992 times.
✗ Branch 2 not taken.
|
22400 | fvGeometry.bindElement(element); |
| 241 |
0/2✗ Branch 1 not taken.
✗ Branch 2 not taken.
|
22400 | elemVolVars.bindElement(element, fvGeometry, sol_); |
| 242 |
0/2✗ Branch 1 not taken.
✗ Branch 2 not taken.
|
22400 | elemFluxVarsCache.bindElement(element, fvGeometry, elemVolVars); |
| 243 | |||
| 244 |
4/4✓ Branch 0 taken 91464 times.
✓ Branch 1 taken 22400 times.
✓ Branch 2 taken 91464 times.
✓ Branch 3 taken 22400 times.
|
136264 | for (const auto& scvf : scvfs(fvGeometry)) |
| 245 | { | ||
| 246 | // iterate through all surfaces and check if the flux at the given position | ||
| 247 | // should be accounted for in the respective surface | ||
| 248 |
2/2✓ Branch 0 taken 192288 times.
✓ Branch 1 taken 91464 times.
|
466680 | for (auto& [name, surfaceData] : surfaces_) |
| 249 | { | ||
| 250 |
2/2✓ Branch 1 taken 1214 times.
✓ Branch 2 taken 191074 times.
|
192288 | if (considerScvf_(scvf, surfaceData)) |
| 251 | { | ||
| 252 | 2428 | const auto result = fluxType(element, fvGeometry, elemVolVars, scvf, elemFluxVarsCache); | |
| 253 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 1214 times.
|
1214 | surfaceData.flux += result; |
| 254 | |||
| 255 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 1214 times.
|
1214 | if (verbose_) |
| 256 | ✗ | std::cout << "At element " << problem_().gridGeometry().elementMapper().index(element) | |
| 257 | ✗ | << ": Flux at face " << scvf.ipGlobal() << ": " << result << " (" << name << ")" << std::endl; | |
| 258 | } | ||
| 259 | } | ||
| 260 | } | ||
| 261 | } | ||
| 262 | 9 | } | |
| 263 | |||
| 264 | //! Check whether a scvf should be considered for the flux calculation | ||
| 265 | ✗ | bool considerScvf_(const SubControlVolumeFace& scvf, const SurfaceData& SurfaceData) const | |
| 266 | { | ||
| 267 | // In order to avoid considering scvfs at the same element intersection (and hence, the corresponding flux) twice, | ||
| 268 | // only use those with a unit outer normal pointing towards positive coordinate direction, | ||
| 269 | // unless the scvf lies on a boundary (then there is no second scvf). | ||
| 270 | ✗ | if (scvf.boundary() || !std::signbit(scvf.unitOuterNormal()[SurfaceData.normalDirectionIndex])) | |
| 271 | ✗ | return intersectsPointGeometry(scvf.ipGlobal(), SurfaceData.surface); | |
| 272 | else | ||
| 273 | ✗ | return false; | |
| 274 | } | ||
| 275 | |||
| 276 | 9 | void snapSurfaceToClosestFace_() | |
| 277 | { | ||
| 278 | using GeometriesEntitySet = Dumux::GeometriesEntitySet<Surface>; | ||
| 279 | 18 | const auto gridView = problem_().gridGeometry().gridView(); | |
| 280 | |||
| 281 |
5/8✓ Branch 0 taken 25 times.
✓ Branch 1 taken 9 times.
✓ Branch 3 taken 25 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 25 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 25 times.
✗ Branch 10 not taken.
|
52 | for (auto& [name, surfaceData] : surfaces_) |
| 282 | { | ||
| 283 |
1/2✓ Branch 1 taken 25 times.
✗ Branch 2 not taken.
|
50 | GeometriesEntitySet entitySet({surfaceData.surface}); |
| 284 |
5/14✓ Branch 1 taken 25 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 25 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 25 times.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✓ Branch 10 taken 25 times.
✓ Branch 12 taken 25 times.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
✗ Branch 15 not taken.
✗ Branch 17 not taken.
✗ Branch 18 not taken.
|
75 | Dumux::BoundingBoxTree<GeometriesEntitySet> geometriesTree(std::make_shared<GeometriesEntitySet>(entitySet)); |
| 285 |
1/2✓ Branch 1 taken 25 times.
✗ Branch 2 not taken.
|
25 | const auto intersectingElements = intersectingEntities( |
| 286 |
1/2✓ Branch 1 taken 25 times.
✗ Branch 2 not taken.
|
25 | problem_().gridGeometry().boundingBoxTree(), geometriesTree |
| 287 | ); | ||
| 288 | |||
| 289 |
2/4✗ Branch 0 not taken.
✓ Branch 1 taken 25 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 25 times.
|
50 | if (intersectingElements.empty()) |
| 290 | { | ||
| 291 | ✗ | if (!nonIntersectingSurfaceIsError_) | |
| 292 | ✗ | continue; | |
| 293 | |||
| 294 | ✗ | std::cout << "surface boundaries: " << std::endl; | |
| 295 | ✗ | printSurfaceBoundaries_(surfaceData.surface); | |
| 296 | |||
| 297 | ✗ | DUNE_THROW(Dune::InvalidStateException, "surface " << name << " does not intersect with any element"); | |
| 298 | } | ||
| 299 | |||
| 300 |
2/7✓ Branch 1 taken 25 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 13 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
|
75 | std::vector<std::size_t> sortedResults; |
| 301 |
3/5✓ Branch 1 taken 13 times.
✓ Branch 2 taken 12 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 13 times.
✗ Branch 5 not taken.
|
25 | sortedResults.reserve(gridView.size(0)); |
| 302 | |||
| 303 |
4/4✓ Branch 0 taken 3976 times.
✓ Branch 1 taken 25 times.
✓ Branch 2 taken 3976 times.
✓ Branch 3 taken 25 times.
|
8027 | for (const auto& i : intersectingElements) |
| 304 |
1/4✓ Branch 1 taken 3976 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
|
3976 | sortedResults.push_back(i.first()); |
| 305 | |||
| 306 | 75 | std::sort(sortedResults.begin(), sortedResults.end()); | |
| 307 | 75 | sortedResults.erase(std::unique( | |
| 308 | sortedResults.begin(), sortedResults.end() | ||
| 309 | 50 | ), sortedResults.end()); | |
| 310 | |||
| 311 | // pick the first intersecting element and make sure the surface snaps to the closest face with the same (or opposite facing) normal vector | ||
| 312 | 25 | GlobalPosition normalVector(0.0); | |
| 313 |
1/2✓ Branch 1 taken 25 times.
✗ Branch 2 not taken.
|
25 | normalVector[surfaceData.normalDirectionIndex] = 1.0; |
| 314 | |||
| 315 |
1/2✓ Branch 1 taken 25 times.
✗ Branch 2 not taken.
|
25 | const auto& firstIntersectingElement = problem_().gridGeometry().element(sortedResults[0]); |
| 316 | 25 | Scalar distance = std::numeric_limits<Scalar>::max(); | |
| 317 | 25 | bool snappingOcurred = false; | |
| 318 | |||
| 319 | 25 | GlobalPosition surfaceLowerLeft = surfaceData.surface.corner(0); | |
| 320 | 25 | GlobalPosition surfaceUpperRight = surfaceData.surface.corner(3); | |
| 321 | |||
| 322 | 25 | bool surfaceAlreadyOnFaces = false; | |
| 323 |
11/14✓ Branch 1 taken 25 times.
✓ Branch 2 taken 52 times.
✓ Branch 4 taken 12 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 12 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 66 times.
✓ Branch 10 taken 12 times.
✓ Branch 11 taken 66 times.
✓ Branch 12 taken 12 times.
✓ Branch 13 taken 34 times.
✓ Branch 14 taken 32 times.
✓ Branch 16 taken 66 times.
✗ Branch 17 not taken.
|
167 | for (const auto& intersection : intersections(gridView, firstIntersectingElement)) |
| 324 | { | ||
| 325 |
2/2✓ Branch 0 taken 79 times.
✓ Branch 1 taken 39 times.
|
118 | if (surfaceAlreadyOnFaces) |
| 326 | continue; | ||
| 327 | |||
| 328 | using std::abs; | ||
| 329 |
7/7✓ Branch 0 taken 23 times.
✓ Branch 1 taken 42 times.
✓ Branch 2 taken 37 times.
✓ Branch 3 taken 42 times.
✓ Branch 4 taken 37 times.
✓ Branch 5 taken 42 times.
✓ Branch 6 taken 14 times.
|
361 | if (abs(1.0 - abs(normalVector * intersection.centerUnitOuterNormal())) < 1e-8) |
| 330 | { | ||
| 331 | |||
| 332 | ✗ | const auto getDistance = [](const auto& p, const auto& geo) | |
| 333 | { | ||
| 334 | if constexpr (dim == 2) | ||
| 335 | ✗ | return distancePointSegment(p, geo); | |
| 336 | else | ||
| 337 | ✗ | return distancePointPolygon(p, geo); | |
| 338 | }; | ||
| 339 | |||
| 340 | 43 | const auto& geo = intersection.geometry(); | |
| 341 |
2/2✓ Branch 2 taken 26 times.
✓ Branch 3 taken 17 times.
|
86 | if (const Scalar d = getDistance(geo.center(), surfaceData.surface); d < 1e-8 * diameter(geo)) |
| 342 | { | ||
| 343 | // no snapping required, face already lies on surface | ||
| 344 | surfaceAlreadyOnFaces = true; | ||
| 345 | snappingOcurred = false; | ||
| 346 | } | ||
| 347 |
2/2✓ Branch 0 taken 18 times.
✓ Branch 1 taken 8 times.
|
26 | else if (d < distance) |
| 348 | { | ||
| 349 | distance = d; | ||
| 350 | snappingOcurred = true; | ||
| 351 | |||
| 352 | // move the surface boundaries | ||
| 353 |
2/2✓ Branch 0 taken 48 times.
✓ Branch 1 taken 18 times.
|
66 | for (int i = 0; i < surfaceData.surface.corners(); ++i) |
| 354 | { | ||
| 355 | 48 | const auto& faceCenter = geo.center(); | |
| 356 | 96 | surfaceLowerLeft[surfaceData.normalDirectionIndex] = faceCenter[surfaceData.normalDirectionIndex]; | |
| 357 | 144 | surfaceUpperRight[surfaceData.normalDirectionIndex] = faceCenter[surfaceData.normalDirectionIndex]; | |
| 358 | } | ||
| 359 | } | ||
| 360 | } | ||
| 361 | } | ||
| 362 | |||
| 363 |
2/2✓ Branch 0 taken 8 times.
✓ Branch 1 taken 17 times.
|
25 | if (snappingOcurred) |
| 364 | { | ||
| 365 |
2/4✓ Branch 2 taken 8 times.
✗ Branch 3 not taken.
✓ Branch 6 taken 8 times.
✗ Branch 7 not taken.
|
24 | std::cout << "\n\nSurface '" << name << "' was automatically snapped to the closest faces" << std::endl; |
| 366 |
1/2✓ Branch 2 taken 8 times.
✗ Branch 3 not taken.
|
16 | std::cout << "Old surface boundaries: " << std::endl; |
| 367 |
1/2✓ Branch 1 taken 8 times.
✗ Branch 2 not taken.
|
8 | printSurfaceBoundaries_(surfaceData.surface); |
| 368 | |||
| 369 | // overwrite the old surface with the new boundaries | ||
| 370 |
2/4✓ Branch 1 taken 8 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 8 times.
✗ Branch 5 not taken.
|
16 | auto inSurfaceAxes = std::move(std::bitset<dimWorld>{}.set()); |
| 371 |
1/2✓ Branch 1 taken 8 times.
✗ Branch 2 not taken.
|
8 | inSurfaceAxes.set(surfaceData.normalDirectionIndex, false); |
| 372 | 8 | surfaceData.surface = Surface{surfaceLowerLeft, surfaceUpperRight, inSurfaceAxes}; | |
| 373 | |||
| 374 |
1/2✓ Branch 2 taken 8 times.
✗ Branch 3 not taken.
|
16 | std::cout << "New surface boundaries: " << std::endl; |
| 375 |
1/2✓ Branch 1 taken 8 times.
✗ Branch 2 not taken.
|
8 | printSurfaceBoundaries_(surfaceData.surface); |
| 376 |
1/2✓ Branch 1 taken 8 times.
✗ Branch 2 not taken.
|
8 | std::cout << std::endl; |
| 377 | } | ||
| 378 | } | ||
| 379 | 9 | } | |
| 380 | |||
| 381 | 16 | void printSurfaceBoundaries_(const Surface& surface) const | |
| 382 | { | ||
| 383 |
2/2✓ Branch 0 taken 44 times.
✓ Branch 1 taken 16 times.
|
60 | for (int i = 0; i < surface.corners(); ++i) |
| 384 | 88 | std::cout << surface.corner(i) << std::endl; | |
| 385 | 16 | } | |
| 386 | |||
| 387 |
9/18✓ Branch 1 taken 9 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 9 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 1214 times.
✗ Branch 8 not taken.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
✓ Branch 12 taken 13 times.
✗ Branch 13 not taken.
✓ Branch 14 taken 12 times.
✓ Branch 15 taken 13 times.
✗ Branch 16 not taken.
✓ Branch 17 taken 12 times.
✓ Branch 18 taken 5 times.
✗ Branch 19 not taken.
✓ Branch 20 taken 4 times.
✗ Branch 21 not taken.
|
1300 | const auto& problem_() const { return gridVariables_.curGridVolVars().problem(); } |
| 388 | |||
| 389 | std::map<std::string, SurfaceData> surfaces_; | ||
| 390 | const GridVariables& gridVariables_; | ||
| 391 | const SolutionVector& sol_; | ||
| 392 | const LocalResidual localResidual_; // store a copy of the local residual | ||
| 393 | bool verbose_; | ||
| 394 | bool nonIntersectingSurfaceIsError_; | ||
| 395 | }; | ||
| 396 | |||
| 397 | } // end namespace Dumux | ||
| 398 | |||
| 399 | #endif | ||
| 400 |