| 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 EmbeddedCoupling | ||
| 10 | * \brief Coupling manager for 1d-3d mixed-dimension method with distributed sources | ||
| 11 | */ | ||
| 12 | #ifndef DUMUX_MULTIDOMAIN_EMBEDDED_COUPLINGMANAGER_1D3D_KERNEL_ROOTSOIL_HH | ||
| 13 | #define DUMUX_MULTIDOMAIN_EMBEDDED_COUPLINGMANAGER_1D3D_KERNEL_ROOTSOIL_HH | ||
| 14 | |||
| 15 | #include <random> | ||
| 16 | #include <vector> | ||
| 17 | #include <array> | ||
| 18 | #include <memory> | ||
| 19 | #include <type_traits> | ||
| 20 | #include <unordered_map> | ||
| 21 | |||
| 22 | #include <dune/common/timer.hh> | ||
| 23 | #include <dune/common/fvector.hh> | ||
| 24 | #include <dune/common/exceptions.hh> | ||
| 25 | #include <dune/geometry/quadraturerules.hh> | ||
| 26 | #include <dune/common/std/type_traits.hh> | ||
| 27 | |||
| 28 | #include <dumux/common/properties.hh> | ||
| 29 | #include <dumux/common/math.hh> | ||
| 30 | #include <dumux/common/indextraits.hh> | ||
| 31 | #include <dumux/geometry/distance.hh> | ||
| 32 | |||
| 33 | #include <dumux/multidomain/embedded/pointsourcedata.hh> | ||
| 34 | #include <dumux/multidomain/embedded/integrationpointsource.hh> | ||
| 35 | #include <dumux/multidomain/embedded/couplingmanagerbase.hh> | ||
| 36 | #include <dumux/multidomain/embedded/circlepoints.hh> | ||
| 37 | #include <dumux/multidomain/embedded/extendedsourcestencil.hh> | ||
| 38 | #include <dumux/multidomain/embedded/cylinderintegration.hh> | ||
| 39 | |||
| 40 | #include <dumux/multidomain/embedded/couplingmanager1d3d.hh> | ||
| 41 | |||
| 42 | namespace Dumux { | ||
| 43 | |||
| 44 | /*! | ||
| 45 | * \ingroup EmbeddedCoupling | ||
| 46 | * \brief Coupling manager for 1d-3d mixed-dimension method with distributed sources | ||
| 47 | * \tparam MDTraits the multidomain traits | ||
| 48 | * \tparam CouplingReconstruction the interface reconstruction operator | ||
| 49 | * | ||
| 50 | * Implements the coupling method described in Koch et al (2022), https://doi.org/10.1016/j.jcp.2021.110823 | ||
| 51 | */ | ||
| 52 | template<class MDTraits, class CouplingReconstruction> | ||
| 53 | class CouplingManagerRootSoilKernel | ||
| 54 | : public EmbeddedCouplingManagerBase<MDTraits, CouplingManagerRootSoilKernel<MDTraits, CouplingReconstruction>> | ||
| 55 | { | ||
| 56 | using ThisType = CouplingManagerRootSoilKernel<MDTraits, CouplingReconstruction>; | ||
| 57 | using ParentType = EmbeddedCouplingManagerBase<MDTraits, ThisType>; | ||
| 58 | using Scalar = typename MDTraits::Scalar; | ||
| 59 | using SolutionVector = typename MDTraits::SolutionVector; | ||
| 60 | using PointSourceData = typename ParentType::PointSourceTraits::PointSourceData; | ||
| 61 | |||
| 62 | static constexpr auto bulkIdx_ = typename MDTraits::template SubDomain<0>::Index(); | ||
| 63 | static constexpr auto lowDimIdx_ = typename MDTraits::template SubDomain<1>::Index(); | ||
| 64 | |||
| 65 | // the sub domain type aliases | ||
| 66 | template<std::size_t id> using SubDomainTypeTag = typename MDTraits::template SubDomain<id>::TypeTag; | ||
| 67 | template<std::size_t id> using Problem = GetPropType<SubDomainTypeTag<id>, Properties::Problem>; | ||
| 68 | template<std::size_t id> using GridGeometry = GetPropType<SubDomainTypeTag<id>, Properties::GridGeometry>; | ||
| 69 | template<std::size_t id> using GridView = typename GridGeometry<id>::GridView; | ||
| 70 | template<std::size_t id> using Element = typename GridView<id>::template Codim<0>::Entity; | ||
| 71 | template<std::size_t id> using GridIndex = typename IndexTraits<GridView<id>>::GridIndex; | ||
| 72 | |||
| 73 | using GlobalPosition = typename Element<bulkIdx_>::Geometry::GlobalCoordinate; | ||
| 74 | |||
| 75 | template<std::size_t id> | ||
| 76 | static constexpr bool isBox() | ||
| 77 | { return GridGeometry<id>::discMethod == DiscretizationMethods::box; } | ||
| 78 | |||
| 79 | static_assert(!isBox<bulkIdx_>() && !isBox<lowDimIdx_>(), "The kernel coupling method is only implemented for the tpfa method"); | ||
| 80 | static_assert(Dune::Capabilities::isCartesian<typename GridView<bulkIdx_>::Grid>::v, "The kernel coupling method is only implemented for structured grids"); | ||
| 81 | |||
| 82 | static constexpr int bulkDim = GridView<bulkIdx_>::dimension; | ||
| 83 | static constexpr int lowDimDim = GridView<lowDimIdx_>::dimension; | ||
| 84 | static constexpr int dimWorld = GridView<bulkIdx_>::dimensionworld; | ||
| 85 | |||
| 86 | // detect if a class (the spatial params) has a kernelWidthFactor() function | ||
| 87 | template <typename T, typename ...Ts> | ||
| 88 | using VariableKernelWidthDetector = decltype(std::declval<T>().kernelWidthFactor(std::declval<Ts>()...)); | ||
| 89 | |||
| 90 | template<class T, typename ...Args> | ||
| 91 | static constexpr bool hasKernelWidthFactor() | ||
| 92 | { return Dune::Std::is_detected<VariableKernelWidthDetector, T, Args...>::value; } | ||
| 93 | |||
| 94 | template <typename T, typename ...Ts> | ||
| 95 | using VariableMinKernelRadDetector = decltype(std::declval<T>().minKernelRadius(std::declval<Ts>()...)); | ||
| 96 | |||
| 97 | template<class T, typename ...Args> | ||
| 98 | static constexpr bool hasMinKernelRadius() | ||
| 99 | { return Dune::Std::is_detected<VariableMinKernelRadDetector, T, Args...>::value; } | ||
| 100 | |||
| 101 | using Reconstruction = CouplingReconstruction; | ||
| 102 | |||
| 103 | public: | ||
| 104 | using ParentType::ParentType; | ||
| 105 | |||
| 106 | static constexpr auto bulkIdx = typename MDTraits::template SubDomain<0>::Index(); | ||
| 107 | static constexpr auto lowDimIdx = typename MDTraits::template SubDomain<1>::Index(); | ||
| 108 | static constexpr auto couplingMode = Embedded1d3dCouplingMode::Kernel{}; | ||
| 109 | |||
| 110 | 1 | void init(std::shared_ptr<Problem<bulkIdx>> bulkProblem, | |
| 111 | std::shared_ptr<Problem<lowDimIdx>> lowDimProblem, | ||
| 112 | std::shared_ptr<Reconstruction> reconstruction, | ||
| 113 | const SolutionVector& curSol) | ||
| 114 | { | ||
| 115 |
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); |
| 116 | 1 | reconstruction_ = reconstruction; | |
| 117 | |||
| 118 | 1 | computeLowDimVolumeFractions(); | |
| 119 | |||
| 120 | 3 | const auto refinement = getParamFromGroup<int>(bulkProblem->paramGroup(), "Grid.Refinement", 0); | |
| 121 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
|
1 | if (refinement > 0) |
| 122 | ✗ | DUNE_THROW(Dune::NotImplemented, "The current intersection detection may likely fail for refined grids."); | |
| 123 | 1 | } | |
| 124 | |||
| 125 | /*! | ||
| 126 | * \brief extend the jacobian pattern of the diagonal block of domain i | ||
| 127 | * by those entries that are not already in the uncoupled pattern | ||
| 128 | * \note Such additional dependencies can arise from the coupling, e.g. if a coupling source | ||
| 129 | * term depends on a non-local average of a quantity of the same domain | ||
| 130 | */ | ||
| 131 | template<std::size_t id, class JacobianPattern> | ||
| 132 | ✗ | void extendJacobianPattern(Dune::index_constant<id> domainI, JacobianPattern& pattern) const | |
| 133 | { | ||
| 134 | ✗ | extendedSourceStencil_.extendJacobianPattern(*this, domainI, pattern); | |
| 135 | ✗ | } | |
| 136 | |||
| 137 | /*! | ||
| 138 | * \brief evaluate additional derivatives of the element residual of a domain with respect | ||
| 139 | * to dofs in the same domain that are not in the regular stencil (per default this is not the case) | ||
| 140 | * \note Such additional dependencies can arise from the coupling, e.g. if a coupling source | ||
| 141 | * term depends on a non-local average of a quantity of the same domain | ||
| 142 | * \note This is the same for box and cc | ||
| 143 | */ | ||
| 144 | template<std::size_t i, class LocalAssemblerI, class JacobianMatrixDiagBlock, class GridVariables> | ||
| 145 | ✗ | void evalAdditionalDomainDerivatives(Dune::index_constant<i> domainI, | |
| 146 | const LocalAssemblerI& localAssemblerI, | ||
| 147 | const typename LocalAssemblerI::LocalResidual::ElementResidualVector&, | ||
| 148 | JacobianMatrixDiagBlock& A, | ||
| 149 | GridVariables& gridVariables) | ||
| 150 | { | ||
| 151 | 5283450 | extendedSourceStencil_.evalAdditionalDomainDerivatives(*this, domainI, localAssemblerI, A, gridVariables); | |
| 152 | ✗ | } | |
| 153 | |||
| 154 | /* \brief Compute integration point point sources and associated data | ||
| 155 | * | ||
| 156 | * This method uses grid glue to intersect the given grids. Over each intersection | ||
| 157 | * we later need to integrate a source term. This method places point sources | ||
| 158 | * at each quadrature point and provides the point source with the necessary | ||
| 159 | * information to compute integrals (quadrature weight and integration element) | ||
| 160 | * \param order The order of the quadrature rule for integration of sources over an intersection | ||
| 161 | * \param verbose If the point source computation is verbose | ||
| 162 | */ | ||
| 163 | 1 | void computePointSourceData(std::size_t order = 1, bool verbose = false) | |
| 164 | { | ||
| 165 | // initialize the maps | ||
| 166 | // do some logging and profiling | ||
| 167 | 1 | Dune::Timer watch; | |
| 168 | 2 | std::cout << "[coupling] Initializing the integration point source data structures..." << std::endl; | |
| 169 | |||
| 170 | // prepare the internal data structures | ||
| 171 | 1 | prepareDataStructures_(); | |
| 172 | 2 | std::cout << "[coupling] Resized data structures." << std::endl; | |
| 173 | |||
| 174 |
2/4✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 1 times.
|
1 | const auto& bulkGridGeometry = this->problem(bulkIdx).gridGeometry(); |
| 175 |
2/4✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
|
1 | const auto& lowDimGridGeometry = this->problem(lowDimIdx).gridGeometry(); |
| 176 | |||
| 177 | // generate a bunch of random vectors and values for | ||
| 178 | // Monte-carlo integration on the cylinder defined by line and radius | ||
| 179 |
3/6✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
✓ Branch 3 taken 1 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 1 times.
✗ Branch 7 not taken.
|
1 | static const double characteristicRelativeLength = getParam<double>("MixedDimension.KernelIntegrationCRL", 0.1); |
| 180 | 1 | EmbeddedCoupling::CylinderIntegration<Scalar> cylIntegration(characteristicRelativeLength, 1); | |
| 181 | |||
| 182 |
3/6✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
✓ Branch 3 taken 1 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 1 times.
✗ Branch 7 not taken.
|
1 | static const bool writeIntegrationPointsToFile = getParam<bool>("MixedDimension.WriteIntegrationPointsToFile", false); |
| 183 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
|
1 | if (writeIntegrationPointsToFile) |
| 184 | { | ||
| 185 | ✗ | std::ofstream ipPointFile("kernel_points.log", std::ios::trunc); | |
| 186 | ✗ | ipPointFile << "x,y,z\n"; | |
| 187 | ✗ | std::cout << "[coupling] Initialized kernel_points.log." << std::endl; | |
| 188 | } | ||
| 189 | |||
| 190 |
4/4✓ Branch 0 taken 1329 times.
✓ Branch 1 taken 1 times.
✓ Branch 2 taken 1329 times.
✓ Branch 3 taken 1 times.
|
1332 | for (const auto& is : intersections(this->glue())) |
| 191 | { | ||
| 192 | // all inside elements are identical... | ||
| 193 | 1329 | const auto& inside = is.targetEntity(0); | |
| 194 | // get the intersection geometry for integrating over it | ||
| 195 |
1/2✓ Branch 1 taken 1329 times.
✗ Branch 2 not taken.
|
1329 | const auto intersectionGeometry = is.geometry(); |
| 196 | // get the Gaussian quadrature rule for the local intersection | ||
| 197 | 2658 | const auto lowDimElementIdx = lowDimGridGeometry.elementMapper().index(inside); | |
| 198 | |||
| 199 | // for each intersection integrate kernel and add: | ||
| 200 | // * 1d: a new point source | ||
| 201 | // * 3d: a new kernel volume source | ||
| 202 |
3/6✗ Branch 0 not taken.
✓ Branch 1 taken 1329 times.
✓ Branch 3 taken 1329 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 1329 times.
✗ Branch 7 not taken.
|
1329 | const auto radius = this->problem(lowDimIdx).spatialParams().radius(lowDimElementIdx); |
| 203 |
3/6✓ Branch 1 taken 1329 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1329 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 1329 times.
✗ Branch 8 not taken.
|
3987 | const auto kernelWidthFactor = kernelWidthFactor_(this->problem(lowDimIdx).spatialParams(), lowDimElementIdx); |
| 204 |
3/6✗ Branch 0 not taken.
✓ Branch 1 taken 1329 times.
✓ Branch 3 taken 1329 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 1329 times.
✗ Branch 7 not taken.
|
1329 | const auto minKernelRadius = minKernelRadius_(this->problem(lowDimIdx).spatialParams(), lowDimElementIdx); |
| 205 | using std::max; | ||
| 206 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 1329 times.
|
1329 | const auto kernelWidth = max(kernelWidthFactor*radius, minKernelRadius); |
| 207 | 1329 | const auto a = intersectionGeometry.corner(0); | |
| 208 | 1329 | const auto b = intersectionGeometry.corner(1); | |
| 209 |
1/2✓ Branch 1 taken 1329 times.
✗ Branch 2 not taken.
|
1329 | cylIntegration.setGeometry(a, b, kernelWidth); |
| 210 | |||
| 211 | // we can have multiple 3d elements per 1d intersection, for evaluation taking one of them should be fine | ||
| 212 |
4/4✓ Branch 0 taken 1329 times.
✓ Branch 1 taken 1329 times.
✓ Branch 2 taken 1329 times.
✓ Branch 3 taken 1329 times.
|
3987 | for (int outsideIdx = 0; outsideIdx < is.numDomainNeighbors(); ++outsideIdx) |
| 213 | { | ||
| 214 | // compute source id | ||
| 215 | // for each id we have | ||
| 216 | // (1) a point source for the 1d domain | ||
| 217 | // (2) a bulk source weight for the each element in the 3d domain (the fraction of the total source/sink for the 1d element with the point source) | ||
| 218 | // TODO: i.e. one integration of the kernel should be enough (for each of the weights it's weight/embeddings) | ||
| 219 | // (3) the flux scaling factor for each outside element, i.e. each id | ||
| 220 | 1329 | const auto id = this->idCounter_++; | |
| 221 | |||
| 222 | // compute the weights for the bulk volume sources | ||
| 223 |
1/2✓ Branch 1 taken 1329 times.
✗ Branch 2 not taken.
|
1329 | const auto& outside = is.domainEntity(outsideIdx); |
| 224 | 2658 | const auto bulkElementIdx = bulkGridGeometry.elementMapper().index(outside); | |
| 225 |
2/4✓ Branch 1 taken 1329 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1329 times.
✗ Branch 5 not taken.
|
2658 | const auto surfaceFactor = computeBulkSource(intersectionGeometry, radius, kernelWidth, id, lowDimElementIdx, bulkElementIdx, cylIntegration, is.numDomainNeighbors()); |
| 226 | |||
| 227 | // place a point source at the intersection center | ||
| 228 | 1329 | const auto center = intersectionGeometry.center(); | |
| 229 |
3/6✓ Branch 1 taken 1329 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1329 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 1329 times.
✗ Branch 8 not taken.
|
3987 | this->pointSources(lowDimIdx).emplace_back(center, id, surfaceFactor, intersectionGeometry.volume(), lowDimElementIdx); |
| 230 |
4/8✓ Branch 1 taken 1329 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1329 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 1329 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 1329 times.
✗ Branch 11 not taken.
|
5316 | this->pointSources(lowDimIdx).back().setEmbeddings(is.numDomainNeighbors()); |
| 231 | |||
| 232 | // pre compute additional data used for the evaluation of | ||
| 233 | // the actual solution dependent source term | ||
| 234 |
1/2✓ Branch 1 taken 1329 times.
✗ Branch 2 not taken.
|
2658 | PointSourceData psData; |
| 235 | |||
| 236 | // lowdim interpolation (evaluate at center) | ||
| 237 | if constexpr (isBox<lowDimIdx>()) | ||
| 238 | { | ||
| 239 | using ShapeValues = std::vector<Dune::FieldVector<Scalar, 1> >; | ||
| 240 | const auto lowDimGeometry = this->problem(lowDimIdx).gridGeometry().element(lowDimElementIdx).geometry(); | ||
| 241 | ShapeValues shapeValues; | ||
| 242 | this->getShapeValues(lowDimIdx, this->problem(lowDimIdx).gridGeometry(), lowDimGeometry, center, shapeValues); | ||
| 243 | psData.addLowDimInterpolation(shapeValues, this->vertexIndices(lowDimIdx, lowDimElementIdx), lowDimElementIdx); | ||
| 244 | } | ||
| 245 | else | ||
| 246 | { | ||
| 247 |
1/2✓ Branch 1 taken 1329 times.
✗ Branch 2 not taken.
|
1329 | psData.addLowDimInterpolation(lowDimElementIdx); |
| 248 | } | ||
| 249 | |||
| 250 | // bulk interpolation (evaluate at center) | ||
| 251 | if constexpr (isBox<bulkIdx>()) | ||
| 252 | { | ||
| 253 | using ShapeValues = std::vector<Dune::FieldVector<Scalar, 1> >; | ||
| 254 | const auto bulkGeometry = this->problem(bulkIdx).gridGeometry().element(bulkElementIdx).geometry(); | ||
| 255 | ShapeValues shapeValues; | ||
| 256 | this->getShapeValues(bulkIdx, this->problem(bulkIdx).gridGeometry(), bulkGeometry, center, shapeValues); | ||
| 257 | psData.addBulkInterpolation(shapeValues, this->vertexIndices(bulkIdx, bulkElementIdx), bulkElementIdx); | ||
| 258 | } | ||
| 259 | else | ||
| 260 | { | ||
| 261 |
1/2✓ Branch 1 taken 1329 times.
✗ Branch 2 not taken.
|
1329 | psData.addBulkInterpolation(bulkElementIdx); |
| 262 | } | ||
| 263 | |||
| 264 | // publish point source data in the global vector | ||
| 265 |
2/4✓ Branch 1 taken 1329 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1329 times.
✗ Branch 5 not taken.
|
2658 | this->pointSourceData().emplace_back(std::move(psData)); |
| 266 |
1/2✓ Branch 1 taken 1329 times.
✗ Branch 2 not taken.
|
1329 | const auto factor = getParam<Scalar>("Problem.Factor", 1.0); |
| 267 |
1/2✓ Branch 2 taken 1329 times.
✗ Branch 3 not taken.
|
3987 | const auto avgMinDist = factor*distancePointGeometry_(0.5*(a+b), outside.geometry()); |
| 268 | 2658 | this->averageDistanceToBulkCell().push_back(avgMinDist); | |
| 269 | |||
| 270 | // export the lowdim coupling stencil | ||
| 271 | // we insert all vertices / elements and make it unique later | ||
| 272 | if constexpr (isBox<bulkIdx>()) | ||
| 273 | { | ||
| 274 | const auto& vertices = this->vertexIndices(bulkIdx, bulkElementIdx); | ||
| 275 | this->couplingStencils(lowDimIdx)[lowDimElementIdx].insert(this->couplingStencils(lowDimIdx)[lowDimElementIdx].end(), | ||
| 276 | vertices.begin(), vertices.end()); | ||
| 277 | } | ||
| 278 | else | ||
| 279 | { | ||
| 280 |
3/6✓ Branch 1 taken 1329 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1329 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 1329 times.
✗ Branch 8 not taken.
|
2658 | this->couplingStencils(lowDimIdx)[lowDimElementIdx].push_back(bulkElementIdx); |
| 281 | } | ||
| 282 | } | ||
| 283 | } | ||
| 284 | |||
| 285 | // make the stencils unique | ||
| 286 | 1 | makeUniqueStencil_(); | |
| 287 | |||
| 288 |
3/6✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 1 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 1 times.
|
3 | if (!this->pointSources(bulkIdx).empty()) |
| 289 | ✗ | DUNE_THROW(Dune::InvalidStateException, "Kernel method shouldn't have point sources in the bulk domain but only volume sources!"); | |
| 290 | |||
| 291 |
2/4✓ Branch 3 taken 1 times.
✗ Branch 4 not taken.
✓ Branch 7 taken 1 times.
✗ Branch 8 not taken.
|
3 | std::cout << "[coupling] Finished preparing manager in " << watch.elapsed() << " seconds." << std::endl; |
| 292 | 1 | } | |
| 293 | |||
| 294 | //! Compute the low dim volume fraction in the bulk domain cells | ||
| 295 | 1 | void computeLowDimVolumeFractions() | |
| 296 | { | ||
| 297 | // resize the storage vector | ||
| 298 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
|
1 | lowDimVolumeInBulkElement_.resize(this->gridView(bulkIdx).size(0)); |
| 299 | // get references to the grid geometries | ||
| 300 |
2/4✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 1 times.
|
1 | const auto& lowDimGridGeometry = this->problem(lowDimIdx).gridGeometry(); |
| 301 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
|
1 | const auto& bulkGridGeometry = this->problem(bulkIdx).gridGeometry(); |
| 302 | |||
| 303 | // compute the low dim volume fractions | ||
| 304 |
4/4✓ Branch 0 taken 1329 times.
✓ Branch 1 taken 1 times.
✓ Branch 2 taken 1329 times.
✓ Branch 3 taken 1 times.
|
1332 | for (const auto& is : intersections(this->glue())) |
| 305 | { | ||
| 306 | // all inside elements are identical... | ||
| 307 | 1329 | const auto& inside = is.targetEntity(0); | |
| 308 | 1329 | const auto intersectionGeometry = is.geometry(); | |
| 309 | 2658 | const auto lowDimElementIdx = lowDimGridGeometry.elementMapper().index(inside); | |
| 310 | |||
| 311 | // compute the volume the low-dim domain occupies in the bulk domain if it were full-dimensional | ||
| 312 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 1329 times.
|
1329 | const auto radius = this->problem(lowDimIdx).spatialParams().radius(lowDimElementIdx); |
| 313 |
4/4✓ Branch 0 taken 1329 times.
✓ Branch 1 taken 1329 times.
✓ Branch 2 taken 1329 times.
✓ Branch 3 taken 1329 times.
|
3987 | for (int outsideIdx = 0; outsideIdx < is.numDomainNeighbors(); ++outsideIdx) |
| 314 | { | ||
| 315 | 1329 | const auto& outside = is.domainEntity(outsideIdx); | |
| 316 | 2658 | const auto bulkElementIdx = bulkGridGeometry.elementMapper().index(outside); | |
| 317 | 3987 | lowDimVolumeInBulkElement_[bulkElementIdx] += intersectionGeometry.volume()*M_PI*radius*radius; | |
| 318 | } | ||
| 319 | } | ||
| 320 | 1 | } | |
| 321 | |||
| 322 | //! return the average distance to the coupled bulk cell center | ||
| 323 | Scalar averageDistance(std::size_t id) const | ||
| 324 | 309506562 | { return averageDistanceToBulkCell_[id]; } | |
| 325 | |||
| 326 | //! Return the local radius at an integration point with identifier id | ||
| 327 | 154753281 | Scalar radius(std::size_t id) const | |
| 328 | { | ||
| 329 |
2/4✗ Branch 0 not taken.
✓ Branch 1 taken 154753281 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 154753281 times.
|
309506562 | const auto& data = this->pointSourceData()[id]; |
| 330 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 154753281 times.
|
154753281 | return this->problem(lowDimIdx).spatialParams().radius(data.lowDimElementIdx()); |
| 331 | } | ||
| 332 | |||
| 333 | //! Return the local radial permeability at an integration point with identifier id | ||
| 334 | 150023814 | Scalar Kr(std::size_t id) const | |
| 335 | { | ||
| 336 |
2/4✗ Branch 0 not taken.
✓ Branch 1 taken 150023814 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 150023814 times.
|
300047628 | const auto& data = this->pointSourceData()[id]; |
| 337 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 150023814 times.
|
150023814 | return this->problem(lowDimIdx).spatialParams().Kr(data.lowDimElementIdx()); |
| 338 | } | ||
| 339 | |||
| 340 | //! The volume the lower dimensional domain occupies in the bulk domain element | ||
| 341 | // For one-dimensional low dim domain we assume radial tubes | ||
| 342 | Scalar lowDimVolume(const Element<bulkIdx>& element) const | ||
| 343 | { | ||
| 344 | const auto eIdx = this->problem(bulkIdx).gridGeometry().elementMapper().index(element); | ||
| 345 | return lowDimVolumeInBulkElement_[eIdx]; | ||
| 346 | } | ||
| 347 | |||
| 348 | //! The volume fraction the lower dimensional domain occupies in the bulk domain element | ||
| 349 | // For one-dimensional low dim domain we assume radial tubes | ||
| 350 | Scalar lowDimVolumeFraction(const Element<bulkIdx>& element) const | ||
| 351 | { | ||
| 352 | const auto totalVolume = element.geometry().volume(); | ||
| 353 | return lowDimVolume(element) / totalVolume; | ||
| 354 | } | ||
| 355 | |||
| 356 | //! return all source ids for a bulk elements | ||
| 357 | const std::vector<std::size_t>& bulkSourceIds(GridIndex<bulkIdx> eIdx, int scvIdx = 0) const | ||
| 358 | 64209819 | { return bulkSourceIds_[eIdx][scvIdx]; } | |
| 359 | |||
| 360 | //! return all source ids for a bulk elements | ||
| 361 | const std::vector<Scalar>& bulkSourceWeights(GridIndex<bulkIdx> eIdx, int scvIdx = 0) const | ||
| 362 | 64209819 | { return bulkSourceWeights_[eIdx][scvIdx]; } | |
| 363 | |||
| 364 | const Reconstruction& reconstruction() const | ||
| 365 | 309506562 | { return *reconstruction_; } | |
| 366 | |||
| 367 | private: | ||
| 368 | |||
| 369 | std::vector<Scalar> averageDistanceToBulkCell_; | ||
| 370 | |||
| 371 | //! Return reference to average distances to bulk cell | ||
| 372 | std::vector<Scalar>& averageDistanceToBulkCell() | ||
| 373 |
1/2✓ Branch 1 taken 1329 times.
✗ Branch 2 not taken.
|
1329 | { return averageDistanceToBulkCell_; } |
| 374 | |||
| 375 | //! compute the kernel distributed sources and add stencils | ||
| 376 | template<class Line, class CylIntegration> | ||
| 377 | 1329 | Scalar computeBulkSource(const Line& line, const Scalar radius, const Scalar kernelWidth, | |
| 378 | std::size_t id, GridIndex<lowDimIdx> lowDimElementIdx, GridIndex<bulkIdx> coupledBulkElementIdx, | ||
| 379 | const CylIntegration& cylIntegration, int embeddings) | ||
| 380 | { | ||
| 381 | // Monte-carlo integration on the cylinder defined by line and radius | ||
| 382 |
6/10✓ Branch 0 taken 1 times.
✓ Branch 1 taken 1328 times.
✓ Branch 3 taken 1 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✓ Branch 6 taken 1 times.
✓ Branch 8 taken 1 times.
✗ Branch 9 not taken.
✓ Branch 11 taken 1 times.
✗ Branch 12 not taken.
|
1329 | static const auto bulkParamGroup = this->problem(bulkIdx).paramGroup(); |
| 383 |
4/6✓ Branch 0 taken 1 times.
✓ Branch 1 taken 1328 times.
✓ Branch 3 taken 1 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 1 times.
✗ Branch 7 not taken.
|
1329 | static const auto min = getParamFromGroup<GlobalPosition>(bulkParamGroup, "Grid.LowerLeft"); |
| 384 |
4/6✓ Branch 0 taken 1 times.
✓ Branch 1 taken 1328 times.
✓ Branch 3 taken 1 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 1 times.
✗ Branch 7 not taken.
|
1329 | static const auto max = getParamFromGroup<GlobalPosition>(bulkParamGroup, "Grid.UpperRight"); |
| 385 |
4/6✓ Branch 0 taken 1 times.
✓ Branch 1 taken 1328 times.
✓ Branch 3 taken 1 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 1 times.
✗ Branch 7 not taken.
|
1329 | static const auto cells = getParamFromGroup<std::array<int, bulkDim>>(bulkParamGroup, "Grid.Cells"); |
| 386 | 1329 | const auto cylSamples = cylIntegration.size(); | |
| 387 | 1329 | const auto& a = line.corner(0); | |
| 388 | 1329 | const auto& b = line.corner(1); | |
| 389 | |||
| 390 | // optionally write debugging / visual output of the integration points | ||
| 391 |
4/6✓ Branch 0 taken 1 times.
✓ Branch 1 taken 1328 times.
✓ Branch 3 taken 1 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 1 times.
✗ Branch 7 not taken.
|
1329 | static const auto writeIntegrationPointsToFile = getParam<bool>("MixedDimension.WriteIntegrationPointsToFile", false); |
| 392 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 1329 times.
|
1329 | if (writeIntegrationPointsToFile) |
| 393 | { | ||
| 394 | ✗ | std::ofstream ipPointFile("kernel_points.log", std::ios::app); | |
| 395 | ✗ | for (int i = 0; i < cylSamples; ++i) | |
| 396 | { | ||
| 397 | ✗ | const auto& point = cylIntegration.integrationPoint(i); | |
| 398 | ✗ | if (const bool hasIntersection = intersectsPointBoundingBox(point, min, max); hasIntersection) | |
| 399 | ✗ | ipPointFile << point[0] << "," << point[1] << "," << point[2] << '\n'; | |
| 400 | } | ||
| 401 | } | ||
| 402 | |||
| 403 | Scalar integral = 0.0; | ||
| 404 |
2/2✓ Branch 0 taken 410661 times.
✓ Branch 1 taken 1329 times.
|
411990 | for (int i = 0; i < cylSamples; ++i) |
| 405 | { | ||
| 406 | 410661 | const auto& point = cylIntegration.integrationPoint(i); | |
| 407 | // TODO: below only works for Cartesian grids with ijk numbering (e.g. level 0 YaspGrid (already fails when refined)) | ||
| 408 | // more general is the bounding box tree solution which always works, however it's much slower | ||
| 409 | //const auto bulkIndices = intersectingEntities(point, this->problem(bulkIdx).gridGeometry().boundingBoxTree(), true); | ||
| 410 |
2/2✓ Branch 1 taken 410249 times.
✓ Branch 2 taken 412 times.
|
410661 | if (const bool hasIntersection = intersectsPointBoundingBox(point, min, max); hasIntersection) |
| 411 | { | ||
| 412 | 410249 | const auto bulkElementIdx = intersectingEntityCartesianGrid(point, min, max, cells); | |
| 413 |
2/4✗ Branch 0 not taken.
✓ Branch 1 taken 410249 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 410249 times.
|
410249 | assert(bulkElementIdx < this->problem(bulkIdx).gridGeometry().gridView().size(0)); |
| 414 |
2/2✓ Branch 1 taken 409957 times.
✓ Branch 2 taken 292 times.
|
410249 | const auto localWeight = evalConstKernel_(a, b, point, radius, kernelWidth)*cylIntegration.integrationElement(i)/Scalar(embeddings); |
| 415 | 410249 | integral += localWeight; | |
| 416 |
14/14✓ Branch 0 taken 409957 times.
✓ Branch 1 taken 292 times.
✓ Branch 2 taken 409957 times.
✓ Branch 3 taken 292 times.
✓ Branch 4 taken 409957 times.
✓ Branch 5 taken 292 times.
✓ Branch 6 taken 409957 times.
✓ Branch 7 taken 292 times.
✓ Branch 8 taken 405451 times.
✓ Branch 9 taken 4506 times.
✓ Branch 10 taken 405451 times.
✓ Branch 11 taken 4506 times.
✓ Branch 12 taken 405451 times.
✓ Branch 13 taken 4506 times.
|
1640996 | if (!bulkSourceIds_[bulkElementIdx][0].empty() && id == bulkSourceIds_[bulkElementIdx][0].back()) |
| 417 | { | ||
| 418 | 1621804 | bulkSourceWeights_[bulkElementIdx][0].back() += localWeight; | |
| 419 | } | ||
| 420 | else | ||
| 421 | { | ||
| 422 | 9596 | bulkSourceIds_[bulkElementIdx][0].emplace_back(id); | |
| 423 | 14394 | bulkSourceWeights_[bulkElementIdx][0].emplace_back(localWeight); | |
| 424 | 4798 | addBulkSourceStencils_(bulkElementIdx, lowDimElementIdx, coupledBulkElementIdx); | |
| 425 | } | ||
| 426 | } | ||
| 427 | } | ||
| 428 | |||
| 429 | // the surface factor (which fraction of the source is inside the domain and needs to be considered) | ||
| 430 | 1329 | const auto length = (a-b).two_norm()/Scalar(embeddings); | |
| 431 | 1329 | return integral/length; | |
| 432 | } | ||
| 433 | |||
| 434 | 1 | void prepareDataStructures_() | |
| 435 | { | ||
| 436 | // clear all internal members like pointsource vectors and stencils | ||
| 437 | // initializes the point source id counter | ||
| 438 | 1 | this->clear(); | |
| 439 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
|
1 | bulkSourceIds_.clear(); |
| 440 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
|
1 | bulkSourceWeights_.clear(); |
| 441 | 2 | extendedSourceStencil_.stencil().clear(); | |
| 442 | |||
| 443 | // precompute the vertex indices for efficiency for the box method | ||
| 444 | 1 | this->precomputeVertexIndices(bulkIdx); | |
| 445 | 1 | this->precomputeVertexIndices(lowDimIdx); | |
| 446 | |||
| 447 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
|
1 | bulkSourceIds_.resize(this->gridView(bulkIdx).size(0)); |
| 448 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
|
1 | bulkSourceWeights_.resize(this->gridView(bulkIdx).size(0)); |
| 449 | |||
| 450 | // intersect the bounding box trees | ||
| 451 | 1 | this->glueGrids(); | |
| 452 | |||
| 453 | // reserve memory for data | ||
| 454 | 4 | this->pointSourceData().reserve(this->glue().size()); | |
| 455 | 3 | this->averageDistanceToBulkCell().reserve(this->glue().size()); | |
| 456 | |||
| 457 | // reserve memory for stencils | ||
| 458 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
|
1 | const auto numBulkElements = this->gridView(bulkIdx).size(0); |
| 459 |
2/2✓ Branch 0 taken 7680 times.
✓ Branch 1 taken 1 times.
|
7681 | for (GridIndex<bulkIdx> bulkElementIdx = 0; bulkElementIdx < numBulkElements; ++bulkElementIdx) |
| 460 | { | ||
| 461 | 15360 | this->couplingStencils(bulkIdx)[bulkElementIdx].reserve(10); | |
| 462 | 15360 | extendedSourceStencil_.stencil()[bulkElementIdx].reserve(10); | |
| 463 | 23040 | bulkSourceIds_[bulkElementIdx][0].reserve(10); | |
| 464 | 23040 | bulkSourceWeights_[bulkElementIdx][0].reserve(10); | |
| 465 | } | ||
| 466 | 1 | } | |
| 467 | |||
| 468 | //! Make the stencils unique | ||
| 469 | 1 | void makeUniqueStencil_() | |
| 470 | { | ||
| 471 | // make extra stencils unique | ||
| 472 |
2/2✓ Branch 0 taken 7680 times.
✓ Branch 1 taken 1 times.
|
7684 | for (auto&& stencil : extendedSourceStencil_.stencil()) |
| 473 | { | ||
| 474 | 23040 | std::sort(stencil.second.begin(), stencil.second.end()); | |
| 475 | 38400 | stencil.second.erase(std::unique(stencil.second.begin(), stencil.second.end()), stencil.second.end()); | |
| 476 | |||
| 477 | // remove the vertices element (box) | ||
| 478 | if constexpr (isBox<bulkIdx>()) | ||
| 479 | { | ||
| 480 | const auto& indices = this->vertexIndices(bulkIdx, stencil.first); | ||
| 481 | stencil.second.erase(std::remove_if(stencil.second.begin(), stencil.second.end(), | ||
| 482 | [&](auto i){ return std::find(indices.begin(), indices.end(), i) != indices.end(); }), | ||
| 483 | stencil.second.end()); | ||
| 484 | } | ||
| 485 | // remove the own element (cell-centered) | ||
| 486 | else | ||
| 487 | { | ||
| 488 | 38400 | stencil.second.erase(std::remove_if(stencil.second.begin(), stencil.second.end(), | |
| 489 |
9/10✓ Branch 0 taken 8 times.
✓ Branch 1 taken 95 times.
✓ Branch 2 taken 12 times.
✓ Branch 3 taken 83 times.
✓ Branch 4 taken 41 times.
✓ Branch 5 taken 112 times.
✓ Branch 6 taken 18 times.
✓ Branch 7 taken 136 times.
✓ Branch 8 taken 210 times.
✗ Branch 9 not taken.
|
715 | [&](auto i){ return i == stencil.first; }), |
| 490 | stencil.second.end()); | ||
| 491 | } | ||
| 492 | } | ||
| 493 | |||
| 494 | // make stencils unique | ||
| 495 | using namespace Dune::Hybrid; | ||
| 496 | 9 | forEach(integralRange(Dune::index_constant<2>{}), [&](const auto domainIdx) | |
| 497 | { | ||
| 498 |
4/4✓ Branch 0 taken 1170 times.
✓ Branch 1 taken 1 times.
✓ Branch 3 taken 7680 times.
✓ Branch 4 taken 1 times.
|
8858 | for (auto&& stencil : this->couplingStencils(domainIdx)) |
| 499 | { | ||
| 500 | 26550 | std::sort(stencil.second.begin(), stencil.second.end()); | |
| 501 | 44250 | stencil.second.erase(std::unique(stencil.second.begin(), stencil.second.end()), stencil.second.end()); | |
| 502 | } | ||
| 503 | }); | ||
| 504 | 1 | } | |
| 505 | |||
| 506 | //! add additional stencil entries for the bulk element | ||
| 507 | 4798 | void addBulkSourceStencils_(GridIndex<bulkIdx> bulkElementIdx, GridIndex<lowDimIdx> coupledLowDimElementIdx, GridIndex<bulkIdx> coupledBulkElementIdx) | |
| 508 | { | ||
| 509 | // add the lowdim element to the coupling stencil of this bulk element | ||
| 510 | if constexpr (isBox<lowDimIdx>()) | ||
| 511 | { | ||
| 512 | const auto& vertices = this->vertexIndices(lowDimIdx, coupledLowDimElementIdx); | ||
| 513 | this->couplingStencils(bulkIdx)[bulkElementIdx].insert(this->couplingStencils(bulkIdx)[bulkElementIdx].end(), | ||
| 514 | vertices.begin(), vertices.end()); | ||
| 515 | |||
| 516 | } | ||
| 517 | else | ||
| 518 | { | ||
| 519 | 9596 | auto& s = this->couplingStencils(bulkIdx)[bulkElementIdx]; | |
| 520 | 4798 | s.push_back(coupledLowDimElementIdx); | |
| 521 | } | ||
| 522 | |||
| 523 | // the extended source stencil, every 3d element with a source is coupled to | ||
| 524 | // the element/dofs where the 3d quantities are measured | ||
| 525 | if constexpr (isBox<bulkIdx>()) | ||
| 526 | { | ||
| 527 | const auto& vertices = this->vertexIndices(bulkIdx, coupledBulkElementIdx); | ||
| 528 | extendedSourceStencil_.stencil()[bulkElementIdx].insert(extendedSourceStencil_.stencil()[bulkElementIdx].end(), | ||
| 529 | vertices.begin(), vertices.end()); | ||
| 530 | } | ||
| 531 | else | ||
| 532 | { | ||
| 533 | 9596 | auto& s = extendedSourceStencil_.stencil()[bulkElementIdx]; | |
| 534 | 4798 | s.push_back(coupledBulkElementIdx); | |
| 535 | } | ||
| 536 | 4798 | } | |
| 537 | |||
| 538 | //! a cylindrical kernel around the segment a->b | ||
| 539 | 410249 | Scalar evalConstKernel_(const GlobalPosition& a, | |
| 540 | const GlobalPosition& b, | ||
| 541 | const GlobalPosition& point, | ||
| 542 | const Scalar R, | ||
| 543 | const Scalar rho) const noexcept | ||
| 544 | { | ||
| 545 | // projection of point onto line a + t*(b-a) | ||
| 546 | 410249 | const auto ab = b - a; | |
| 547 | 410249 | const auto t = (point - a)*ab/ab.two_norm2(); | |
| 548 | |||
| 549 | // return 0 if we are outside cylinder | ||
| 550 |
2/4✓ Branch 0 taken 410249 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 410249 times.
✗ Branch 3 not taken.
|
410249 | if (t < 0.0 || t > 1.0) |
| 551 | return 0.0; | ||
| 552 | |||
| 553 | // compute distance | ||
| 554 | 410249 | auto proj = a; proj.axpy(t, ab); | |
| 555 | 410249 | const auto radiusSquared = (proj - point).two_norm2(); | |
| 556 | |||
| 557 |
1/2✓ Branch 0 taken 410249 times.
✗ Branch 1 not taken.
|
410249 | if (radiusSquared > rho*rho) |
| 558 | return 0.0; | ||
| 559 | |||
| 560 | 410249 | return 1.0/(M_PI*rho*rho); | |
| 561 | } | ||
| 562 | |||
| 563 | /*! | ||
| 564 | * \brief Get the kernel width factor from the spatial params (if possible) | ||
| 565 | */ | ||
| 566 | template<class SpatialParams> | ||
| 567 | auto kernelWidthFactor_(const SpatialParams& spatialParams, unsigned int eIdx) | ||
| 568 | -> std::enable_if_t<hasKernelWidthFactor<SpatialParams, unsigned int>(), Scalar> | ||
| 569 | { return spatialParams.kernelWidthFactor(eIdx); } | ||
| 570 | |||
| 571 | /*! | ||
| 572 | * \brief Get the kernel width factor (constant) from the input file (if possible) | ||
| 573 | */ | ||
| 574 | template<class SpatialParams> | ||
| 575 | ✗ | auto kernelWidthFactor_(const SpatialParams& spatialParams, unsigned int eIdx) | |
| 576 | -> std::enable_if_t<!hasKernelWidthFactor<SpatialParams, unsigned int>(), Scalar> | ||
| 577 | { | ||
| 578 | ✗ | static const Scalar kernelWidthFactor = getParam<Scalar>("MixedDimension.KernelWidthFactor"); | |
| 579 | ✗ | return kernelWidthFactor; | |
| 580 | } | ||
| 581 | |||
| 582 | /*! | ||
| 583 | * \brief Get the minimal kernel radius from the spatial params (if possible) | ||
| 584 | */ | ||
| 585 | template<class SpatialParams> | ||
| 586 | auto minKernelRadius_(const SpatialParams& spatialParams, unsigned int eIdx) | ||
| 587 | -> std::enable_if_t<hasMinKernelRadius<SpatialParams, unsigned int>(), Scalar> | ||
| 588 | { return spatialParams.minKernelRadius(eIdx); } | ||
| 589 | |||
| 590 | /*! | ||
| 591 | * \brief Get the kernel width factor (constant) from the input file (if possible) | ||
| 592 | */ | ||
| 593 | template<class SpatialParams> | ||
| 594 | ✗ | auto minKernelRadius_(const SpatialParams& spatialParams, unsigned int eIdx) | |
| 595 | -> std::enable_if_t<!hasMinKernelRadius<SpatialParams, unsigned int>(), Scalar> | ||
| 596 | { | ||
| 597 | ✗ | static const Scalar minKernelRadius = getParam<Scalar>("MixedDimension.MinKernelRadius", 0.0); | |
| 598 | ✗ | return minKernelRadius; | |
| 599 | } | ||
| 600 | |||
| 601 | template<class Geometry> | ||
| 602 | 1329 | typename Geometry::ctype distancePointGeometry_(const typename Geometry::GlobalCoordinate& p, | |
| 603 | const Geometry& geometry, | ||
| 604 | std::size_t integrationOrder = 2) | ||
| 605 | { | ||
| 606 | 1329 | typename Geometry::ctype avgDist = 0.0; | |
| 607 | 6645 | Dune::FieldVector<double, 2> center = { geometry.center()[0], geometry.center()[1] }; | |
| 608 | 3987 | Dune::FieldVector<double, 2> p2 = { p[0], p[1] }; | |
| 609 | 1329 | const auto factor = getParam<double>("Problem.DistancePointGeometryFactor", 1.0); | |
| 610 | 2658 | avgDist = (center-p2).two_norm()*factor; | |
| 611 | 1329 | return avgDist; | |
| 612 | } | ||
| 613 | |||
| 614 | //! the extended source stencil object for kernel coupling | ||
| 615 | EmbeddedCoupling::ExtendedSourceStencil<ThisType> extendedSourceStencil_; | ||
| 616 | //! vector for the volume fraction of the lowdim domain in the bulk domain cells | ||
| 617 | std::vector<Scalar> lowDimVolumeInBulkElement_; | ||
| 618 | //! kernel sources to integrate for each bulk element | ||
| 619 | std::vector<std::array<std::vector<std::size_t>, isBox<bulkIdx>() ? 1<<bulkDim : 1>> bulkSourceIds_; | ||
| 620 | //! the integral of the kernel for each point source / integration point, i.e. weight for the source | ||
| 621 | std::vector<std::array<std::vector<Scalar>, isBox<bulkIdx>() ? 1<<bulkDim : 1>> bulkSourceWeights_; | ||
| 622 | |||
| 623 | //! source reconstruction scheme | ||
| 624 | std::shared_ptr<Reconstruction> reconstruction_; | ||
| 625 | }; | ||
| 626 | |||
| 627 | } // end namespace Dumux | ||
| 628 | |||
| 629 | #endif | ||
| 630 |