| 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 | * \author Timo Koch | ||
| 10 | * \ingroup EmbeddedCoupling | ||
| 11 | * \brief Coupling manager for low-dimensional domains embedded in the bulk | ||
| 12 | * domain with spatially resolved interface. | ||
| 13 | */ | ||
| 14 | |||
| 15 | #ifndef DUMUX_MULTIDOMAIN_EMBEDDED_COUPLINGMANAGER_1D3D_PROJECTION_HH | ||
| 16 | #define DUMUX_MULTIDOMAIN_EMBEDDED_COUPLINGMANAGER_1D3D_PROJECTION_HH | ||
| 17 | |||
| 18 | #include <vector> | ||
| 19 | #include <algorithm> | ||
| 20 | #include <string> | ||
| 21 | |||
| 22 | #include <dune/common/timer.hh> | ||
| 23 | #include <dune/common/exceptions.hh> | ||
| 24 | #include <dune/common/fvector.hh> | ||
| 25 | #include <dune/geometry/quadraturerules.hh> | ||
| 26 | |||
| 27 | #include <dumux/common/tag.hh> | ||
| 28 | #include <dumux/common/properties.hh> | ||
| 29 | #include <dumux/common/math.hh> | ||
| 30 | #include <dumux/common/indextraits.hh> | ||
| 31 | |||
| 32 | #include <dumux/geometry/refinementquadraturerule.hh> | ||
| 33 | #include <dumux/geometry/triangulation.hh> | ||
| 34 | #include <dumux/geometry/grahamconvexhull.hh> | ||
| 35 | |||
| 36 | #include <dumux/multidomain/embedded/couplingmanagerbase.hh> | ||
| 37 | #include <dumux/multidomain/embedded/localrefinementquadrature.hh> | ||
| 38 | |||
| 39 | namespace Dumux { | ||
| 40 | |||
| 41 | // some implementation details of the projection coupling manager | ||
| 42 | namespace Detail { | ||
| 43 | |||
| 44 | /*! | ||
| 45 | * \ingroup EmbeddedCoupling | ||
| 46 | * \brief Segment representation of a 1d network grid | ||
| 47 | * | ||
| 48 | * We use separate segment representation which is fast to iterate over and which | ||
| 49 | * can be used to find the unique mapping from points on coupled bulk surface facets | ||
| 50 | * to the network centerline (i.e. a 1D dof on the centerline). | ||
| 51 | */ | ||
| 52 | template <class GridGeometry> | ||
| 53 | 2 | class SegmentNetwork | |
| 54 | { | ||
| 55 | using GridView = typename GridGeometry::GridView; | ||
| 56 | using GlobalPosition = typename GridView::template Codim<0>::Geometry::GlobalCoordinate; | ||
| 57 | using GridIndex = typename IndexTraits<GridView>::GridIndex; | ||
| 58 | |||
| 59 | 132 | struct Segment | |
| 60 | { | ||
| 61 | GlobalPosition a, b; | ||
| 62 | double r; | ||
| 63 | GridIndex eIdx; | ||
| 64 | }; | ||
| 65 | |||
| 66 | public: | ||
| 67 | template <class RadiusFunction> | ||
| 68 | 1 | SegmentNetwork(const GridGeometry& gg, const RadiusFunction& radiusFunction) | |
| 69 | 1 | { | |
| 70 | 1 | const auto& gridView = gg.gridView(); | |
| 71 |
1/2✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
|
1 | segments_.resize(gridView.size(0)); |
| 72 |
5/6✓ Branch 1 taken 44 times.
✓ Branch 2 taken 1 times.
✓ Branch 3 taken 44 times.
✓ Branch 4 taken 1 times.
✓ Branch 6 taken 44 times.
✗ Branch 7 not taken.
|
89 | for (const auto &element : elements(gridView)) |
| 73 | { | ||
| 74 |
1/2✓ Branch 1 taken 44 times.
✗ Branch 2 not taken.
|
44 | const auto geometry = element.geometry(); |
| 75 | 88 | const auto eIdx = gg.elementMapper().index(element); | |
| 76 | 88 | const auto radius = radiusFunction(eIdx); | |
| 77 | |||
| 78 | 132 | segments_[eIdx] = Segment{ geometry.corner(0), geometry.corner(1), radius, eIdx }; | |
| 79 | } | ||
| 80 | |||
| 81 |
3/6✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 1 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 1 times.
✗ Branch 9 not taken.
|
3 | std::cout << "-- Extracted " << segments_.size() << " segments.\n"; |
| 82 | 1 | } | |
| 83 | |||
| 84 | //! the segment radius | ||
| 85 | auto radius(std::size_t eIdx) const | ||
| 86 | 1644 | { return segments_[eIdx].r; } | |
| 87 | |||
| 88 | //! the segment with index eIdx | ||
| 89 | const auto& segment(std::size_t eIdx) const | ||
| 90 | { return segments_[eIdx]; } | ||
| 91 | |||
| 92 | //! the segments | ||
| 93 | const auto& segments() const | ||
| 94 | { return segments_; } | ||
| 95 | |||
| 96 | //! find the closest segment (the segment from which the signed distance to its surface is the closest to zero) | ||
| 97 | //! \note minDist is a positive distance | ||
| 98 | 134892 | auto findClosestSegmentSurface(const GlobalPosition& globalPos) const | |
| 99 | { | ||
| 100 | 134892 | GridIndex eIdx = 0; | |
| 101 | 134892 | double minSignedDist = std::numeric_limits<double>::max(); | |
| 102 | |||
| 103 |
4/4✓ Branch 0 taken 5935248 times.
✓ Branch 1 taken 134892 times.
✓ Branch 2 taken 5935248 times.
✓ Branch 3 taken 134892 times.
|
6339924 | for (const auto &segment : segments_) |
| 104 | { | ||
| 105 | 5935248 | const auto signedDist = computeSignedDistance_(globalPos, segment); | |
| 106 |
2/2✓ Branch 0 taken 2081694 times.
✓ Branch 1 taken 3853554 times.
|
5935248 | if (signedDist < minSignedDist) |
| 107 | { | ||
| 108 | 2081694 | minSignedDist = signedDist; | |
| 109 | 2081694 | eIdx = segment.eIdx; | |
| 110 | } | ||
| 111 | } | ||
| 112 | |||
| 113 | using std::abs; | ||
| 114 | 404676 | return std::make_tuple(eIdx, abs(minSignedDist)); | |
| 115 | } | ||
| 116 | |||
| 117 | //! same as overload with taking a position but only look at the segments specified by the index vector | ||
| 118 | template <class IndexRange> | ||
| 119 | auto findClosestSegmentSurface(const GlobalPosition& globalPos, | ||
| 120 | const IndexRange& segIndices) const | ||
| 121 | { | ||
| 122 | GridIndex eIdx = 0; | ||
| 123 | double minSignedDist = std::numeric_limits<double>::max(); | ||
| 124 | |||
| 125 | for (const auto index : segIndices) | ||
| 126 | { | ||
| 127 | const auto &segment = segments_[index]; | ||
| 128 | const auto signedDist = computeSignedDistance_(globalPos, segment); | ||
| 129 | if (signedDist < minSignedDist) | ||
| 130 | { | ||
| 131 | minSignedDist = signedDist; | ||
| 132 | eIdx = segment.eIdx; | ||
| 133 | } | ||
| 134 | } | ||
| 135 | |||
| 136 | using std::abs; | ||
| 137 | return std::make_tuple(eIdx, abs(minSignedDist)); | ||
| 138 | } | ||
| 139 | |||
| 140 | // Compute the projection point of p onto the segment connecting a->b | ||
| 141 | GlobalPosition projectionPointOnSegment(const GlobalPosition& p, std::size_t idx) const | ||
| 142 | { | ||
| 143 | 6148 | const auto& segment = segments_[idx]; | |
| 144 | 3074 | return projectionPointOnSegment_(p, segment.a, segment.b); | |
| 145 | } | ||
| 146 | |||
| 147 | private: | ||
| 148 | // Compute the projection of p onto the segment | ||
| 149 | // returns the closest end point if the orthogonal projection onto the line implied by the segment | ||
| 150 | // is outside the segment | ||
| 151 | 5938322 | GlobalPosition projectionPointOnSegment_(const GlobalPosition& p, const GlobalPosition& a, const GlobalPosition& b) const | |
| 152 | { | ||
| 153 | 5938322 | const auto v = b - a; | |
| 154 | 5938322 | const auto w = p - a; | |
| 155 | |||
| 156 | 5938322 | const auto proj1 = v*w; | |
| 157 |
2/2✓ Branch 0 taken 2641205 times.
✓ Branch 1 taken 3297117 times.
|
5938322 | if (proj1 <= 0.0) |
| 158 | 3297117 | return a; | |
| 159 | |||
| 160 | 2641205 | const auto proj2 = v.two_norm2(); | |
| 161 |
2/2✓ Branch 0 taken 2368097 times.
✓ Branch 1 taken 273108 times.
|
2641205 | if (proj2 <= proj1) |
| 162 | 2368097 | return b; | |
| 163 | |||
| 164 | 273108 | const auto t = proj1 / proj2; | |
| 165 | 273108 | auto x = a; | |
| 166 | 273108 | x.axpy(t, v); | |
| 167 | 273108 | return x; | |
| 168 | } | ||
| 169 | |||
| 170 | // Compute the signed (!) distance to a cylinder segment surface and the projection point onto the centerline | ||
| 171 | 5935248 | auto computeSignedDistance_(const GlobalPosition& p, const Segment& segment) const | |
| 172 | { | ||
| 173 | 5935248 | const auto projPoint = projectionPointOnSegment_(p, segment.a, segment.b); | |
| 174 | 5935248 | return (p-projPoint).two_norm()-segment.r; | |
| 175 | } | ||
| 176 | |||
| 177 | std::vector<Segment> segments_; | ||
| 178 | }; | ||
| 179 | |||
| 180 | /*! | ||
| 181 | * \ingroup EmbeddedCoupling | ||
| 182 | * \brief Get the closest segment for a given surface point | ||
| 183 | * \note The point does not necessarily be exactly on the virtual surface | ||
| 184 | * implied by the networks radius function | ||
| 185 | */ | ||
| 186 | template<class Network> | ||
| 187 | class NetworkIndicatorFunction | ||
| 188 | { | ||
| 189 | public: | ||
| 190 | 1 | NetworkIndicatorFunction(const Network& n) | |
| 191 | 1 | : network_(n) | |
| 192 | {} | ||
| 193 | |||
| 194 | template<class Pos> | ||
| 195 | ✗ | auto operator()(const Pos& pos) const | |
| 196 | { | ||
| 197 | 44204 | const auto& [closestSegmentIdx, _] = network_.findClosestSegmentSurface(pos); | |
| 198 | 44204 | return closestSegmentIdx; | |
| 199 | } | ||
| 200 | |||
| 201 | template <class Pos, class HintContainer> | ||
| 202 | auto operator()(const Pos& pos, const HintContainer& hints) const | ||
| 203 | { | ||
| 204 | const auto& [closestSegmentIdx, _] = network_.findClosestSegmentSurface(pos, hints); | ||
| 205 | return closestSegmentIdx; | ||
| 206 | } | ||
| 207 | |||
| 208 | private: | ||
| 209 | const Network& network_; | ||
| 210 | }; | ||
| 211 | |||
| 212 | /*! | ||
| 213 | * \ingroup EmbeddedCoupling | ||
| 214 | * \brief Simple legacy VTK writer for outputting debug data on the coupling interface | ||
| 215 | */ | ||
| 216 | ✗ | class DebugIntersectionVTKOutput | |
| 217 | { | ||
| 218 | public: | ||
| 219 | ✗ | void push(const std::array<Dune::FieldVector<double, 3>, 3>& corners, double data) | |
| 220 | { | ||
| 221 | ✗ | vtkCells_.emplace_back(std::array<std::size_t, 3>{ | |
| 222 | ✗ | { vtkPoints_.size(), vtkPoints_.size()+1, vtkPoints_.size()+2 } | |
| 223 | ✗ | }); | |
| 224 | ✗ | for (const auto& c : corners) | |
| 225 | ✗ | vtkPoints_.push_back(c); | |
| 226 | ✗ | vtkCellData_.push_back(data); | |
| 227 | ✗ | } | |
| 228 | |||
| 229 | ✗ | void write(const std::string& name) | |
| 230 | { | ||
| 231 | ✗ | std::ofstream intersectionFile(name); | |
| 232 | ✗ | intersectionFile << "# vtk DataFile Version 2.0\n"; | |
| 233 | ✗ | intersectionFile << "Really cool intersection data\n"; | |
| 234 | ✗ | intersectionFile << "ASCII\n"; | |
| 235 | ✗ | intersectionFile << "DATASET UNSTRUCTURED_GRID\n"; | |
| 236 | ✗ | intersectionFile << "POINTS " << vtkPoints_.size() << " float\n"; | |
| 237 | ✗ | for (const auto& p : vtkPoints_) | |
| 238 | ✗ | intersectionFile << p[0] << " " << p[1] << " " << p[2] << "\n"; | |
| 239 | ✗ | intersectionFile << "\nCELLS " << vtkCells_.size() << " " << vtkCells_.size()*4 << "\n"; | |
| 240 | ✗ | for (const auto& c : vtkCells_) | |
| 241 | ✗ | intersectionFile << "3 " << c[0] << " " << c[1] << " " << c[2] << "\n"; | |
| 242 | ✗ | intersectionFile << "\nCELL_TYPES " << vtkCells_.size() << "\n"; | |
| 243 | ✗ | for (int i = 0; i < vtkCells_.size(); ++i) | |
| 244 | ✗ | intersectionFile << "5\n"; | |
| 245 | ✗ | intersectionFile << "\nCELL_DATA " << vtkCells_.size() << "\nSCALARS index float\nLOOKUP_TABLE default\n"; | |
| 246 | ✗ | for (const auto& c : vtkCellData_) | |
| 247 | ✗ | intersectionFile << c << "\n"; | |
| 248 | ✗ | } | |
| 249 | |||
| 250 | private: | ||
| 251 | std::vector<Dune::FieldVector<double, 3>> vtkPoints_; | ||
| 252 | std::vector<std::array<std::size_t, 3>> vtkCells_; | ||
| 253 | std::vector<double> vtkCellData_; | ||
| 254 | }; | ||
| 255 | |||
| 256 | } // end namespace Detail | ||
| 257 | |||
| 258 | namespace Embedded1d3dCouplingMode { | ||
| 259 | struct Projection : public Utility::Tag<Projection> { | ||
| 260 | static std::string name() { return "projection"; } | ||
| 261 | }; | ||
| 262 | |||
| 263 | inline constexpr Projection projection{}; | ||
| 264 | } // end namespace Embedded1d3dCouplingMode | ||
| 265 | |||
| 266 | // forward declaration | ||
| 267 | template<class MDTraits, class CouplingMode> | ||
| 268 | class Embedded1d3dCouplingManager; | ||
| 269 | |||
| 270 | /*! | ||
| 271 | * \ingroup EmbeddedCoupling | ||
| 272 | * \brief Manages the coupling between bulk elements and lower dimensional elements | ||
| 273 | * | ||
| 274 | * The method is described in detail in the article | ||
| 275 | * "Projection-based resolved interface mixed-dimension method for embedded tubular network systems" | ||
| 276 | * by Timo Koch (2021) available at https://arxiv.org/abs/2106.06358 | ||
| 277 | * | ||
| 278 | * The network is represented as a line segment network. The bulk domain explicitly resolves the wall | ||
| 279 | * of the tube network (e.g. blood vessel wall, outer root wall), so this generally requires an unstructured grid. | ||
| 280 | * The coupling term coupling the PDEs on the network and in the bulk domain is integrated over the | ||
| 281 | * coupling surface and each integration point couples to quantities evaluated at the closest point on the graph. | ||
| 282 | * There is a unique mapping from every point on the virtual tube surface to the tube centerline, therefore | ||
| 283 | * this coupling is determined by mapping to the closest point on the virtual surface and then evaluating | ||
| 284 | * the 1D quantity on the mapped centerline point. (The inverse mapping may be non-unique.) | ||
| 285 | * | ||
| 286 | * The original paper also includes an algorithm for computing exact intersection for the surface quadrature | ||
| 287 | * in case of simple straight vessels. However, as a comparison in the paper shows that the suggested | ||
| 288 | * approximate quadrature is exact enough (configured by parameter MixedDimension.Projection.SimplexIntegrationRefine | ||
| 289 | * specifying the number of (adaptive local) virtual refinements of a surface facet for integration), only | ||
| 290 | * the approximate general purpose algorithm is included in this implementation. Only one quadrature point | ||
| 291 | * will be added for each surface element and coupled 1D dof including all contribution evaluated by | ||
| 292 | * local refinement. That means the virtual refinement doesn't increase the number of integration points and | ||
| 293 | * additionally is also optimized by using local adaptive refinement only in the places where the index mapping | ||
| 294 | * to the 1D degree of freedom is still multivalent. | ||
| 295 | * | ||
| 296 | * Coupling sources are implement in terms of point sources at quadrature points to make it easy to reuse code | ||
| 297 | * written for other 1D-3D coupling manager modes. | ||
| 298 | * | ||
| 299 | * This algorithm can be configured by several parameters | ||
| 300 | * (although this is usually not necessary as there are sensible defaults): | ||
| 301 | * | ||
| 302 | * - MixedDimension.Projection.SimplexIntegrationRefine | ||
| 303 | * number of virtual refinement steps to determine coupled surface area | ||
| 304 | * - MixedDimension.Projection.EnableIntersectionOutput | ||
| 305 | * set to true to enable debug VTK output for intersections | ||
| 306 | * - MixedDimension.Projection.EstimateNumberOfPointSources | ||
| 307 | * provide an estimate for the expected number of coupling points for memory allocation | ||
| 308 | * - MixedDimension.Projection.CoupledRadiusFactor | ||
| 309 | * threshold distance in which to search for coupled elements (specified as multiple of radius, default 0.1) | ||
| 310 | * - MixedDimension.Projection.CoupledAngleFactor | ||
| 311 | * angle threshold in which to search for coupled elements (angle in radians from surface normal vector, default 0.3) | ||
| 312 | * - MixedDimension.Projection.ConsiderFacesWithinBoundingBoxCoupled | ||
| 313 | * determines if all 3D boundary facets within the mesh bounding box should be considered as coupling faces | ||
| 314 | * - MixedDimension.Projection.CoupledBoundingBoxShrinkingFactor | ||
| 315 | * if ConsiderFacesWithinBoundingBoxCoupled=true shrink the bounding box in all directions by this factor | ||
| 316 | * | ||
| 317 | */ | ||
| 318 | template<class MDTraits> | ||
| 319 | class Embedded1d3dCouplingManager<MDTraits, Embedded1d3dCouplingMode::Projection> | ||
| 320 | : public EmbeddedCouplingManagerBase<MDTraits, Embedded1d3dCouplingManager<MDTraits, Embedded1d3dCouplingMode::Projection>> | ||
| 321 | { | ||
| 322 | using ThisType = Embedded1d3dCouplingManager<MDTraits, Embedded1d3dCouplingMode::Projection>; | ||
| 323 | using ParentType = EmbeddedCouplingManagerBase<MDTraits, ThisType>; | ||
| 324 | |||
| 325 | using Scalar = typename MDTraits::Scalar; | ||
| 326 | using SolutionVector = typename MDTraits::SolutionVector; | ||
| 327 | using PointSourceData = typename ParentType::PointSourceTraits::PointSourceData; | ||
| 328 | |||
| 329 | static constexpr auto bulkIdx = typename MDTraits::template SubDomain<0>::Index(); | ||
| 330 | static constexpr auto lowDimIdx = typename MDTraits::template SubDomain<1>::Index(); | ||
| 331 | |||
| 332 | // the sub domain type aliases | ||
| 333 | template<std::size_t id> using SubDomainTypeTag = typename MDTraits::template SubDomain<id>::TypeTag; | ||
| 334 | template<std::size_t id> using Problem = GetPropType<SubDomainTypeTag<id>, Properties::Problem>; | ||
| 335 | template<std::size_t id> using GridGeometry = GetPropType<SubDomainTypeTag<id>, Properties::GridGeometry>; | ||
| 336 | template<std::size_t id> using GridView = typename GridGeometry<id>::GridView; | ||
| 337 | template<std::size_t id> using Element = typename GridView<id>::template Codim<0>::Entity; | ||
| 338 | template<std::size_t id> using GridIndex = typename IndexTraits<GridView<id>>::GridIndex; | ||
| 339 | |||
| 340 | static constexpr int lowDimDim = GridView<lowDimIdx>::dimension; | ||
| 341 | static constexpr int bulkDim = GridView<bulkIdx>::dimension; | ||
| 342 | static constexpr int dimWorld = GridView<bulkIdx>::dimensionworld; | ||
| 343 | |||
| 344 | template<std::size_t id> | ||
| 345 | static constexpr bool isBox() | ||
| 346 | { return GridGeometry<id>::discMethod == DiscretizationMethods::box; } | ||
| 347 | |||
| 348 | using GlobalPosition = typename Element<bulkIdx>::Geometry::GlobalCoordinate; | ||
| 349 | |||
| 350 | using SegmentNetwork = Detail::SegmentNetwork<GridGeometry<lowDimIdx>>; | ||
| 351 | |||
| 352 | public: | ||
| 353 | static constexpr Embedded1d3dCouplingMode::Projection couplingMode{}; | ||
| 354 | |||
| 355 |
5/10✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 1 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 1 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 1 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 1 times.
✗ Branch 10 not taken.
|
1 | using ParentType::ParentType; |
| 356 | |||
| 357 | 1 | void init(std::shared_ptr<Problem<bulkIdx>> bulkProblem, | |
| 358 | std::shared_ptr<Problem<lowDimIdx>> lowDimProblem, | ||
| 359 | const SolutionVector& curSol) | ||
| 360 | { | ||
| 361 |
3/10✓ Branch 3 taken 1 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 1 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 1 times.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
|
2 | ParentType::init(bulkProblem, lowDimProblem, curSol); |
| 362 | 1 | } | |
| 363 | |||
| 364 | /* \brief Compute integration point point sources and associated data | ||
| 365 | * This method computes a source term on the explicitly given surface of the network grid | ||
| 366 | * by integrating over the elements of the surface, projecting the integration points onto | ||
| 367 | * the centerlines to evaluate 1d quantities and set point sources (minimum distance projection) | ||
| 368 | * | ||
| 369 | * \param order Unused in this algorithm! (passed from default 1D3D coupling manager) | ||
| 370 | * \param verbose If the source computation is verbose | ||
| 371 | */ | ||
| 372 | 1 | void computePointSourceData(std::size_t order = 1, bool verbose = false) | |
| 373 | { | ||
| 374 | // initialize the maps | ||
| 375 | // do some logging and profiling | ||
| 376 | 1 | Dune::Timer watch; | |
| 377 | 2 | std::cout << "Initializing the coupling manager (projection)" << std::endl; | |
| 378 | |||
| 379 | // debug VTK output | ||
| 380 | 1 | const bool enableIntersectionOutput = getParam<bool>("MixedDimension.Projection.EnableIntersectionOutput", false); | |
| 381 |
1/4✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
|
1 | std::unique_ptr<Detail::DebugIntersectionVTKOutput> vtkCoupledFaces = |
| 382 | enableIntersectionOutput ? std::make_unique<Detail::DebugIntersectionVTKOutput>() : nullptr; | ||
| 383 |
2/6✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✓ Branch 6 taken 1 times.
|
2 | std::unique_ptr<Detail::DebugIntersectionVTKOutput> vtkIntersections = |
| 384 | enableIntersectionOutput ? std::make_unique<Detail::DebugIntersectionVTKOutput>() : nullptr; | ||
| 385 | |||
| 386 | // temporarily represent the network domain as a list of segments | ||
| 387 | // a more efficient data structure | ||
| 388 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
|
1 | const auto& lowDimProblem = this->problem(lowDimIdx); |
| 389 | 1 | const auto& lowDimFvGridGeometry = lowDimProblem.gridGeometry(); | |
| 390 | 1 | const auto& lowDimGridView = lowDimFvGridGeometry.gridView(); | |
| 391 |
1/4✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
|
1 | localAreaFactor_.resize(lowDimGridView.size(0), 0.0); |
| 392 | 1 | const auto radiusFunc = [&](const GridIndex<lowDimIdx> eIdx) { | |
| 393 | 132 | return lowDimProblem.spatialParams().radius(eIdx); | |
| 394 | }; | ||
| 395 |
2/4✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 1 times.
|
2 | const auto network = SegmentNetwork{ |
| 396 | lowDimFvGridGeometry, radiusFunc | ||
| 397 | }; | ||
| 398 | |||
| 399 | // precompute the vertex indices for efficiency (box method) | ||
| 400 |
1/2✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
|
1 | this->precomputeVertexIndices(bulkIdx); |
| 401 | 1 | this->precomputeVertexIndices(lowDimIdx); | |
| 402 | |||
| 403 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
|
1 | const auto& bulkProblem = this->problem(bulkIdx); |
| 404 |
1/2✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
|
1 | const auto& bulkFvGridGeometry = bulkProblem.gridGeometry(); |
| 405 |
1/2✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
|
1 | const auto& bulkGridView = bulkFvGridGeometry.gridView(); |
| 406 |
2/4✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
|
1 | bulkElementMarker_.assign(bulkGridView.size(0), 0); |
| 407 |
2/6✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
|
1 | bulkVertexMarker_.assign(bulkGridView.size(bulkDim), 0); |
| 408 | |||
| 409 | // construct a bounding box that may be used to determine coupled faces | ||
| 410 | // the shrinking factor shrinks this bounding box | ||
| 411 |
2/4✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
✓ Branch 3 taken 1 times.
✗ Branch 4 not taken.
|
2 | static const auto bBoxShrinking |
| 412 |
1/2✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
|
2 | = getParam<Scalar>("MixedDimension.Projection.CoupledBoundingBoxShrinkingFactor", 1e-2); |
| 413 | 2 | const GlobalPosition threshold( | |
| 414 | 4 | bBoxShrinking*( bulkFvGridGeometry.bBoxMax()-bulkFvGridGeometry.bBoxMin() ).two_norm() | |
| 415 | ); | ||
| 416 | 2 | const auto bBoxMinSmall = bulkFvGridGeometry.bBoxMin() + threshold; | |
| 417 | 2 | const auto bBoxMaxSmall = bulkFvGridGeometry.bBoxMax() - threshold; | |
| 418 | 2281 | auto insideBBox = [bBoxMin=bBoxMinSmall, bBoxMax=bBoxMaxSmall](const GlobalPosition& point) -> bool | |
| 419 | { | ||
| 420 | static constexpr Scalar eps_ = 1.0e-7; | ||
| 421 |
4/4✓ Branch 0 taken 2147 times.
✓ Branch 1 taken 133 times.
✓ Branch 2 taken 2147 times.
✓ Branch 3 taken 133 times.
|
4560 | const Scalar eps0 = eps_*(bBoxMax[0] - bBoxMin[0]); |
| 422 |
4/4✓ Branch 0 taken 2147 times.
✓ Branch 1 taken 133 times.
✓ Branch 2 taken 2147 times.
✓ Branch 3 taken 133 times.
|
4560 | const Scalar eps1 = eps_*(bBoxMax[1] - bBoxMin[1]); |
| 423 |
4/4✓ Branch 0 taken 2147 times.
✓ Branch 1 taken 133 times.
✓ Branch 2 taken 2147 times.
✓ Branch 3 taken 133 times.
|
4560 | const Scalar eps2 = eps_*(bBoxMax[2] - bBoxMin[2]); |
| 424 |
10/10✓ Branch 0 taken 2147 times.
✓ Branch 1 taken 133 times.
✓ Branch 2 taken 2147 times.
✓ Branch 3 taken 133 times.
✓ Branch 4 taken 2035 times.
✓ Branch 5 taken 112 times.
✓ Branch 6 taken 2035 times.
✓ Branch 7 taken 112 times.
✓ Branch 8 taken 2035 times.
✓ Branch 9 taken 112 times.
|
4560 | return (bBoxMin[0] - eps0 <= point[0] && point[0] <= bBoxMax[0] + eps0 && |
| 425 |
12/12✓ Branch 0 taken 1914 times.
✓ Branch 1 taken 121 times.
✓ Branch 2 taken 1914 times.
✓ Branch 3 taken 121 times.
✓ Branch 4 taken 1914 times.
✓ Branch 5 taken 121 times.
✓ Branch 6 taken 1776 times.
✓ Branch 7 taken 138 times.
✓ Branch 8 taken 1776 times.
✓ Branch 9 taken 138 times.
✓ Branch 10 taken 1776 times.
✓ Branch 11 taken 138 times.
|
6105 | bBoxMin[1] - eps1 <= point[1] && point[1] <= bBoxMax[1] + eps1 && |
| 426 |
14/14✓ Branch 0 taken 2147 times.
✓ Branch 1 taken 133 times.
✓ Branch 2 taken 1744 times.
✓ Branch 3 taken 32 times.
✓ Branch 4 taken 1744 times.
✓ Branch 5 taken 32 times.
✓ Branch 6 taken 1744 times.
✓ Branch 7 taken 32 times.
✓ Branch 8 taken 286 times.
✓ Branch 9 taken 1458 times.
✓ Branch 10 taken 286 times.
✓ Branch 11 taken 1458 times.
✓ Branch 12 taken 286 times.
✓ Branch 13 taken 1458 times.
|
7608 | bBoxMin[2] - eps2 <= point[2] && point[2] <= bBoxMax[2] + eps2); |
| 427 | }; | ||
| 428 | |||
| 429 | // setup simplex "quadrature" rule | ||
| 430 |
2/4✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
✓ Branch 3 taken 1 times.
✗ Branch 4 not taken.
|
2 | static const auto quadSimplexRefineMaxLevel |
| 431 |
1/2✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
|
2 | = getParam<std::size_t>("MixedDimension.Projection.SimplexIntegrationRefine", 4); |
| 432 | 2 | const auto indicator = Detail::NetworkIndicatorFunction { network }; | |
| 433 | |||
| 434 | // estimate number of point sources for memory allocation | ||
| 435 |
2/4✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
✓ Branch 3 taken 1 times.
✗ Branch 4 not taken.
|
2 | static const auto estimateNumberOfPointSources |
| 436 |
3/6✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 1 times.
✗ Branch 8 not taken.
|
3 | = getParam<std::size_t>("MixedDimension.Projection.EstimateNumberOfPointSources", bulkFvGridGeometry.gridView().size(0)); |
| 437 | |||
| 438 |
2/4✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
|
2 | this->pointSources(bulkIdx).reserve(estimateNumberOfPointSources); |
| 439 |
2/4✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
|
2 | this->pointSources(lowDimIdx).reserve(estimateNumberOfPointSources); |
| 440 |
2/4✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
|
2 | this->pointSourceData().reserve(estimateNumberOfPointSources); |
| 441 | |||
| 442 | // we go over the bulk surface intersection, check if we are coupled, | ||
| 443 | // and if this is the case, we add integration point sources using | ||
| 444 | // the local simplex refinement quadrature procedure (see paper Koch 2021) | ||
| 445 |
8/12✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 1 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 8224 times.
✓ Branch 10 taken 1 times.
✓ Branch 11 taken 8224 times.
✓ Branch 12 taken 1 times.
✓ Branch 14 taken 8224 times.
✗ Branch 15 not taken.
|
16449 | for (const auto& element : elements(bulkGridView)) |
| 446 | { | ||
| 447 |
2/4✓ Branch 1 taken 8224 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 8224 times.
✗ Branch 5 not taken.
|
16448 | const auto bulkElementIdx = bulkFvGridGeometry.elementMapper().index(element); |
| 448 |
2/4✓ Branch 1 taken 8224 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 8224 times.
✗ Branch 5 not taken.
|
16448 | const auto bulkGeometry = element.geometry(); |
| 449 |
1/2✓ Branch 1 taken 8224 times.
✗ Branch 2 not taken.
|
8224 | const auto refElement = Dune::referenceElement(element); |
| 450 |
3/8✓ Branch 1 taken 8224 times.
✗ Branch 2 not taken.
✓ Branch 5 taken 41120 times.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 11 taken 32896 times.
✗ Branch 12 not taken.
|
90464 | for (const auto& intersection : intersections(bulkGridView, element)) |
| 451 | { | ||
| 452 | // skip inner intersection | ||
| 453 |
2/2✓ Branch 1 taken 2280 times.
✓ Branch 2 taken 30616 times.
|
32896 | if (!intersection.boundary()) |
| 454 | 31384 | continue; | |
| 455 | |||
| 456 | // check if we are close enough to the coupling interface | ||
| 457 |
7/8✓ Branch 1 taken 2280 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 1512 times.
✓ Branch 4 taken 768 times.
✓ Branch 5 taken 1512 times.
✓ Branch 6 taken 768 times.
✓ Branch 7 taken 1512 times.
✓ Branch 8 taken 768 times.
|
2280 | const auto [isCoupled, closestSegmentIdx, minDist] = isCoupled_(network, intersection, insideBBox); |
| 458 |
2/2✓ Branch 0 taken 1512 times.
✓ Branch 1 taken 768 times.
|
2280 | if (!isCoupled) |
| 459 | continue; | ||
| 460 | |||
| 461 | // we found an intersection on the coupling interface: mark as coupled element | ||
| 462 |
1/2✓ Branch 1 taken 1512 times.
✗ Branch 2 not taken.
|
1512 | bulkElementMarker_[bulkElementIdx] = 1; |
| 463 | |||
| 464 |
1/2✓ Branch 1 taken 1512 times.
✗ Branch 2 not taken.
|
3024 | const auto interfaceGeometry = intersection.geometry(); |
| 465 |
4/4✓ Branch 0 taken 4536 times.
✓ Branch 1 taken 1512 times.
✓ Branch 2 taken 4536 times.
✓ Branch 3 taken 1512 times.
|
12096 | for (int i = 0; i < interfaceGeometry.corners(); ++i) |
| 466 | { | ||
| 467 |
1/2✓ Branch 1 taken 4536 times.
✗ Branch 2 not taken.
|
4536 | const auto localVIdx = refElement.subEntity(intersection.indexInInside(), 1, i, bulkDim); |
| 468 |
2/4✓ Branch 1 taken 4536 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 4536 times.
✗ Branch 5 not taken.
|
9072 | const auto vIdx = bulkFvGridGeometry.vertexMapper().subIndex(element, localVIdx, bulkDim); |
| 469 | 9072 | bulkVertexMarker_[vIdx] = 1; | |
| 470 | } | ||
| 471 | |||
| 472 | // debug graphical intersection output | ||
| 473 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 1512 times.
|
1512 | if (enableIntersectionOutput) |
| 474 | ✗ | vtkCoupledFaces->push({ | |
| 475 | { interfaceGeometry.corner(0), interfaceGeometry.corner(1), interfaceGeometry.corner(2)} | ||
| 476 | }, closestSegmentIdx); | ||
| 477 | |||
| 478 |
2/6✓ Branch 1 taken 1512 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 1512 times.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
|
4536 | const auto quadSimplex = LocalRefinementSimplexQuadrature{ interfaceGeometry, quadSimplexRefineMaxLevel, indicator }; |
| 479 | |||
| 480 | // iterate over all quadrature points (children simplices) | ||
| 481 |
4/4✓ Branch 0 taken 2252 times.
✓ Branch 1 taken 1512 times.
✓ Branch 2 taken 2252 times.
✓ Branch 3 taken 1512 times.
|
6788 | for (const auto& qp : quadSimplex) |
| 482 | { | ||
| 483 | 4504 | const auto surfacePoint = interfaceGeometry.global(qp.position()); | |
| 484 | 2252 | const auto closestSegmentIdx = qp.indicator(); | |
| 485 | 2252 | const auto projPoint = network.projectionPointOnSegment(surfacePoint, closestSegmentIdx); | |
| 486 | |||
| 487 |
1/2✓ Branch 1 taken 2252 times.
✗ Branch 2 not taken.
|
2252 | addPointSourceAtIP_(bulkGeometry, interfaceGeometry, qp, bulkElementIdx, closestSegmentIdx, surfacePoint, projPoint); |
| 488 |
1/2✓ Branch 1 taken 2252 times.
✗ Branch 2 not taken.
|
2252 | addStencilEntries_(bulkElementIdx, closestSegmentIdx); |
| 489 | } | ||
| 490 | } | ||
| 491 | } | ||
| 492 | |||
| 493 | // make stencils unique | ||
| 494 | using namespace Dune::Hybrid; | ||
| 495 | 9 | forEach(integralRange(Dune::index_constant<2>{}), [&](const auto domainIdx) | |
| 496 | { | ||
| 497 |
4/4✓ Branch 0 taken 44 times.
✓ Branch 1 taken 1 times.
✓ Branch 3 taken 1456 times.
✓ Branch 4 taken 1 times.
|
1508 | for (auto&& stencil : this->couplingStencils(domainIdx)) |
| 498 | { | ||
| 499 | 4500 | std::sort(stencil.second.begin(), stencil.second.end()); | |
| 500 | 3000 | stencil.second.erase( | |
| 501 | 4500 | std::unique(stencil.second.begin(), stencil.second.end()), | |
| 502 | stencil.second.end() | ||
| 503 | ); | ||
| 504 | } | ||
| 505 | }); | ||
| 506 | |||
| 507 | // debug VTK output | ||
| 508 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
|
1 | if (enableIntersectionOutput) |
| 509 | ✗ | vtkCoupledFaces->write("coupledfaces.vtk"); | |
| 510 | |||
| 511 | // compare theoretical cylinder surface to actual discrete surface | ||
| 512 | // of the coupled bulk interface facets | ||
| 513 |
5/6✓ Branch 1 taken 44 times.
✓ Branch 2 taken 1 times.
✓ Branch 3 taken 44 times.
✓ Branch 4 taken 1 times.
✓ Branch 6 taken 44 times.
✗ Branch 7 not taken.
|
89 | for (const auto& element : elements(lowDimGridView)) |
| 514 | { | ||
| 515 |
1/2✓ Branch 1 taken 44 times.
✗ Branch 2 not taken.
|
44 | const auto length = element.geometry().volume(); |
| 516 | 88 | const auto eIdx = lowDimFvGridGeometry.elementMapper().index(element); | |
| 517 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 44 times.
|
44 | const auto radius = this->problem(lowDimIdx).spatialParams().radius(eIdx); |
| 518 | 44 | const auto cylinderSurface = 2.0*M_PI*radius*length; | |
| 519 | 44 | localAreaFactor_[eIdx] /= cylinderSurface; | |
| 520 | 44 | localAreaFactor_[eIdx] -= 1.0; | |
| 521 | } | ||
| 522 | |||
| 523 |
2/4✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
|
2 | this->pointSources(bulkIdx).shrink_to_fit(); |
| 524 |
2/4✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
|
2 | this->pointSources(lowDimIdx).shrink_to_fit(); |
| 525 |
1/2✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
|
1 | this->pointSourceData().shrink_to_fit(); |
| 526 | |||
| 527 |
4/8✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 1 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 1 times.
✗ Branch 9 not taken.
✓ Branch 12 taken 1 times.
✗ Branch 13 not taken.
|
5 | std::cout << "-- Coupling at " << this->pointSourceData().size() << " integration points." << std::endl; |
| 528 |
3/6✓ Branch 3 taken 1 times.
✗ Branch 4 not taken.
✓ Branch 7 taken 1 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 1 times.
✗ Branch 10 not taken.
|
3 | std::cout << "-- Finished initialization after " << watch.elapsed() << " seconds." << std::endl; |
| 529 | 1 | } | |
| 530 | |||
| 531 | /*! | ||
| 532 | * \brief Methods to be accessed by the subproblems | ||
| 533 | */ | ||
| 534 | // \{ | ||
| 535 | |||
| 536 | //! Return a reference to the bulk problem | ||
| 537 | Scalar radius(std::size_t id) const | ||
| 538 | { | ||
| 539 | const auto& data = this->pointSourceData()[id]; | ||
| 540 | return this->problem(lowDimIdx).spatialParams().radius(data.lowDimElementIdx()); | ||
| 541 | } | ||
| 542 | |||
| 543 | // \} | ||
| 544 | |||
| 545 | //! Marker that is non-zero for bulk elements that are coupled | ||
| 546 | const std::vector<int>& bulkElementMarker() const | ||
| 547 |
1/2✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
|
1 | { return bulkElementMarker_; } |
| 548 | |||
| 549 | //! Marker that is non-zero for bulk vertices that belong to coupled surface facets | ||
| 550 | const std::vector<int>& bulkVertexMarker() const | ||
| 551 | { return bulkVertexMarker_; } | ||
| 552 | |||
| 553 | //! For each 1D element, the ratio of associated discrete bulk surface area | ||
| 554 | //! to the theoretical bulk surface area resulting from assuming a cylinder-shaped segment | ||
| 555 | //! (excluding cylinder caps). | ||
| 556 | const std::vector<Scalar>& localAreaFactor() const | ||
| 557 | { return localAreaFactor_; } | ||
| 558 | |||
| 559 | private: | ||
| 560 | std::vector<int> bulkElementMarker_, bulkVertexMarker_; | ||
| 561 | std::vector<double> localAreaFactor_; | ||
| 562 | |||
| 563 | // add a point source containing all data needed to evaluate | ||
| 564 | // the coupling source term from both domains | ||
| 565 | template<class BulkGeometry, class InterfaceGeometry, class QP> | ||
| 566 | 2252 | void addPointSourceAtIP_(const BulkGeometry& bulkGeometry, | |
| 567 | const InterfaceGeometry& ifGeometry, | ||
| 568 | const QP& qp, | ||
| 569 | const GridIndex<bulkIdx> bulkElementIdx, | ||
| 570 | const GridIndex<lowDimIdx> lowDimElementIdx, | ||
| 571 | const GlobalPosition& surfacePoint, | ||
| 572 | const GlobalPosition& projPoint) | ||
| 573 | { | ||
| 574 | // add a new point source id (for the associated data) | ||
| 575 | 2252 | const auto id = this->idCounter_++; | |
| 576 | |||
| 577 | // add a bulk point source | ||
| 578 | 2252 | const auto qpweight = qp.weight(); | |
| 579 | 4504 | const auto integrationElement = ifGeometry.integrationElement(qp.position()); | |
| 580 | 4504 | this->pointSources(bulkIdx).emplace_back(surfacePoint, id, qpweight, integrationElement, bulkElementIdx); | |
| 581 | |||
| 582 | // find the corresponding projection point and 1d element | ||
| 583 | 4504 | this->pointSources(lowDimIdx).emplace_back(projPoint, id, qpweight, integrationElement, lowDimElementIdx); | |
| 584 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 2252 times.
|
2252 | localAreaFactor_[lowDimElementIdx] += qpweight*integrationElement; |
| 585 | |||
| 586 |
2/4✗ Branch 0 not taken.
✓ Branch 1 taken 2252 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 2252 times.
|
2252 | const auto& bulkFvGridGeometry = this->problem(bulkIdx).gridGeometry(); |
| 587 |
2/4✗ Branch 0 not taken.
✓ Branch 1 taken 2252 times.
✓ Branch 3 taken 2252 times.
✗ Branch 4 not taken.
|
2252 | const auto& lowDimFvGridGeometry = this->problem(lowDimIdx).gridGeometry(); |
| 588 | |||
| 589 | // pre compute additional data used for the evaluation of | ||
| 590 | // the actual solution dependent source term | ||
| 591 |
1/2✓ Branch 1 taken 2252 times.
✗ Branch 2 not taken.
|
4504 | PointSourceData psData; |
| 592 | |||
| 593 | if constexpr (isBox<lowDimIdx>()) | ||
| 594 | { | ||
| 595 | std::vector<Dune::FieldVector<Scalar, 1> > shapeValues; | ||
| 596 | const auto lowDimGeometry = this->problem(lowDimIdx).gridGeometry().element(lowDimElementIdx).geometry(); | ||
| 597 | this->getShapeValues(lowDimIdx, lowDimFvGridGeometry, lowDimGeometry, projPoint, shapeValues); | ||
| 598 | psData.addLowDimInterpolation(shapeValues, this->vertexIndices(lowDimIdx, lowDimElementIdx), lowDimElementIdx); | ||
| 599 | } | ||
| 600 | else | ||
| 601 | { | ||
| 602 |
1/2✓ Branch 1 taken 2252 times.
✗ Branch 2 not taken.
|
2252 | psData.addLowDimInterpolation(lowDimElementIdx); |
| 603 | } | ||
| 604 | |||
| 605 | if constexpr (isBox<bulkIdx>()) | ||
| 606 | { | ||
| 607 |
1/4✓ Branch 1 taken 2252 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
|
2252 | std::vector<Dune::FieldVector<Scalar, 1> > shapeValues; |
| 608 |
1/2✓ Branch 1 taken 2252 times.
✗ Branch 2 not taken.
|
2252 | this->getShapeValues(bulkIdx, bulkFvGridGeometry, bulkGeometry, surfacePoint, shapeValues); |
| 609 |
3/6✓ Branch 1 taken 2252 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2252 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 2252 times.
✗ Branch 7 not taken.
|
6756 | psData.addBulkInterpolation(shapeValues, this->vertexIndices(bulkIdx, bulkElementIdx), bulkElementIdx); |
| 610 | } | ||
| 611 | else | ||
| 612 | { | ||
| 613 | psData.addBulkInterpolation(bulkElementIdx); | ||
| 614 | } | ||
| 615 | |||
| 616 | // publish point source data in the global vector | ||
| 617 |
2/4✓ Branch 1 taken 2252 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2252 times.
✗ Branch 5 not taken.
|
4504 | this->pointSourceData().emplace_back(std::move(psData)); |
| 618 | 2252 | } | |
| 619 | |||
| 620 | 2252 | void addStencilEntries_(const GridIndex<bulkIdx> bulkElementIdx, const GridIndex<lowDimIdx> lowDimElementIdx) | |
| 621 | { | ||
| 622 | // export the bulk coupling stencil | ||
| 623 | if constexpr (isBox<lowDimIdx>()) | ||
| 624 | { | ||
| 625 | const auto& vertices = this->vertexIndices(lowDimIdx, lowDimElementIdx); | ||
| 626 | this->couplingStencils(bulkIdx)[bulkElementIdx].insert( | ||
| 627 | this->couplingStencils(bulkIdx)[bulkElementIdx].end(), | ||
| 628 | vertices.begin(), vertices.end() | ||
| 629 | ); | ||
| 630 | |||
| 631 | } | ||
| 632 | else | ||
| 633 | { | ||
| 634 | 4504 | this->couplingStencils(bulkIdx)[bulkElementIdx].push_back(lowDimElementIdx); | |
| 635 | } | ||
| 636 | |||
| 637 | // export the lowdim coupling stencil | ||
| 638 | if constexpr (isBox<bulkIdx>()) | ||
| 639 | { | ||
| 640 | 2252 | const auto& vertices = this->vertexIndices(bulkIdx, bulkElementIdx); | |
| 641 | 4504 | this->couplingStencils(lowDimIdx)[lowDimElementIdx].insert( | |
| 642 | 4504 | this->couplingStencils(lowDimIdx)[lowDimElementIdx].end(), | |
| 643 | vertices.begin(), vertices.end() | ||
| 644 | ); | ||
| 645 | } | ||
| 646 | else | ||
| 647 | { | ||
| 648 | this->couplingStencils(lowDimIdx)[lowDimElementIdx].push_back(bulkElementIdx); | ||
| 649 | } | ||
| 650 | 2252 | } | |
| 651 | |||
| 652 | template<class BoundaryIntersection, class BBoxCheck> | ||
| 653 | 2280 | auto isCoupled_(const SegmentNetwork& network, const BoundaryIntersection& intersection, const BBoxCheck& insideBBox) const | |
| 654 | { | ||
| 655 |
1/2✓ Branch 2 taken 2280 times.
✗ Branch 3 not taken.
|
2280 | const auto center = intersection.geometry().center(); |
| 656 |
4/4✓ Branch 1 taken 1 times.
✓ Branch 2 taken 2279 times.
✓ Branch 3 taken 1 times.
✓ Branch 4 taken 2279 times.
|
2280 | const auto [eIdx, minDist] = network.findClosestSegmentSurface(center); |
| 657 | |||
| 658 | // all elements in the bounding box are coupled if this is enabled | ||
| 659 | // this is a simplified check that can be enabled for box-shaped domains | ||
| 660 |
4/6✓ Branch 0 taken 1 times.
✓ Branch 1 taken 2279 times.
✓ Branch 3 taken 1 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 1 times.
✗ Branch 7 not taken.
|
2280 | static const bool considerBBox = getParam<bool>("MixedDimension.Projection.ConsiderFacesWithinBoundingBoxCoupled", false); |
| 661 |
3/4✓ Branch 0 taken 2280 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 1458 times.
✓ Branch 3 taken 822 times.
|
2280 | if (considerBBox && insideBBox(center)) |
| 662 | 1458 | return std::make_tuple(true, eIdx, minDist); | |
| 663 | |||
| 664 | // general coupling face detection algorithm: | ||
| 665 | // if the distance is under a certain threshold we are coupled (and the angle is not very big) | ||
| 666 |
4/6✓ Branch 0 taken 1 times.
✓ Branch 1 taken 821 times.
✓ Branch 3 taken 1 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 1 times.
✗ Branch 7 not taken.
|
822 | static const auto rFactor = getParam<Scalar>("MixedDimension.Projection.CoupledRadiusFactor", 0.1); |
| 667 |
4/6✓ Branch 0 taken 1 times.
✓ Branch 1 taken 821 times.
✓ Branch 3 taken 1 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 1 times.
✗ Branch 7 not taken.
|
822 | static const auto spFactor = getParam<Scalar>("MixedDimension.Projection.CoupledAngleFactor", 0.3); |
| 668 | 822 | const auto projPoint = network.projectionPointOnSegment(center, eIdx); | |
| 669 | 822 | const auto distVec = projPoint - center; | |
| 670 |
5/6✓ Branch 0 taken 54 times.
✓ Branch 1 taken 768 times.
✓ Branch 2 taken 54 times.
✓ Branch 3 taken 768 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 54 times.
|
1698 | const bool isCoupled = minDist < rFactor * network.radius(eIdx) && distVec * intersection.centerUnitOuterNormal() > distVec.two_norm() * spFactor; |
| 671 | 822 | return std::make_tuple(isCoupled, eIdx, minDist); | |
| 672 | } | ||
| 673 | }; | ||
| 674 | |||
| 675 | //! we support multithreaded assembly | ||
| 676 | template<class MDTraits> | ||
| 677 | struct CouplingManagerSupportsMultithreadedAssembly<Embedded1d3dCouplingManager<MDTraits, Embedded1d3dCouplingMode::Projection>> | ||
| 678 | : public std::true_type {}; | ||
| 679 | |||
| 680 | } // end namespace Dumux | ||
| 681 | |||
| 682 | #endif | ||
| 683 |