diff --git a/CHANGELOG.md b/CHANGELOG.md
index 1aa0908577a6a12ec3e243112f902d29079607c9..fdca7fada9833c2568db60c8543d7c120601ce32 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -40,6 +40,8 @@ for (const auto& element : elements(gridGeometry.gridView()))
fvGeometry.bind(element);
```
+- __Embedded coupling__: Add a coupling manager for the 1D-3D projection based scheme with resolved interface introduced in Koch 2021 (https://arxiv.org/abs/2106.06358)
+
### Immediate interface changes not allowing/requiring a deprecation period:
- __Virtual interface of GridDataTransfer__: The `GridDataTransfer` abstract base class now required the Grid type as a template argument. Furthermore, the `store` and `reconstruct` interface functions do now expect the grid as a function argument. This allows to correctly update grid geometries and corresponding mapper (see "Construction and update of GridGeometries changed" above in the changelog)
- `PengRobinsonMixture::computeMolarVolumes` has been removed without deprecation. It was used nowhere and did not translate.
diff --git a/dumux/multidomain/embedded/couplingmanager1d3d.hh b/dumux/multidomain/embedded/couplingmanager1d3d.hh
index 8845f9c82a102bd457629ba0285af936738c6743..a7e1bdfd91d14228d67b77f34b5bfd0833408b5c 100644
--- a/dumux/multidomain/embedded/couplingmanager1d3d.hh
+++ b/dumux/multidomain/embedded/couplingmanager1d3d.hh
@@ -41,5 +41,6 @@ class Embedded1d3dCouplingManager;
#include <dumux/multidomain/embedded/couplingmanager1d3d_average.hh>
#include <dumux/multidomain/embedded/couplingmanager1d3d_surface.hh>
#include <dumux/multidomain/embedded/couplingmanager1d3d_kernel.hh>
+#include <dumux/multidomain/embedded/couplingmanager1d3d_projection.hh>
#endif
diff --git a/dumux/multidomain/embedded/couplingmanager1d3d_projection.hh b/dumux/multidomain/embedded/couplingmanager1d3d_projection.hh
new file mode 100644
index 0000000000000000000000000000000000000000..ca29f716c47c5a4c98a3c2876d1030b112d2d46a
--- /dev/null
+++ b/dumux/multidomain/embedded/couplingmanager1d3d_projection.hh
@@ -0,0 +1,675 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ * See the file COPYING for full copying permissions. *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 3 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+/*!
+ * \file
+ * \author Timo Koch
+ * \ingroup EmbeddedCoupling
+ * \brief Coupling manager for low-dimensional domains embedded in the bulk
+ * domain with spatially resolved interface.
+ */
+
+#ifndef DUMUX_MULTIDOMAIN_EMBEDDED_COUPLINGMANAGER_1D3D_PROJECTION_HH
+#define DUMUX_MULTIDOMAIN_EMBEDDED_COUPLINGMANAGER_1D3D_PROJECTION_HH
+
+#include <vector>
+#include <algorithm>
+#include <string>
+
+#include <dune/common/timer.hh>
+#include <dune/common/exceptions.hh>
+#include <dune/common/fvector.hh>
+#include <dune/geometry/quadraturerules.hh>
+
+#include <dumux/common/tag.hh>
+#include <dumux/common/properties.hh>
+#include <dumux/common/math.hh>
+#include <dumux/common/indextraits.hh>
+
+#include <dumux/geometry/refinementquadraturerule.hh>
+#include <dumux/geometry/triangulation.hh>
+#include <dumux/geometry/grahamconvexhull.hh>
+
+#include <dumux/multidomain/embedded/couplingmanagerbase.hh>
+#include <dumux/multidomain/embedded/localrefinementquadrature.hh>
+
+namespace Dumux {
+
+// some implementation details of the projection coupling manager
+namespace Detail {
+
+/*!
+ * \ingroup EmbeddedCoupling
+ * \brief Segment representation of a 1d network grid
+ *
+ * We use separate segment representation which is fast to iterate over and which
+ * can be used to find the unique mapping from points on coupled bulk surface facets
+ * to the network centerline (i.e. a 1D dof on the centerline).
+ */
+template <class GridGeometry>
+class SegmentNetwork
+{
+ using GridView = typename GridGeometry::GridView;
+ using GlobalPosition = typename GridView::template Codim<0>::Geometry::GlobalCoordinate;
+ using GridIndex = typename IndexTraits<GridView>::GridIndex;
+
+ struct Segment
+ {
+ GlobalPosition a, b;
+ double r;
+ GridIndex eIdx;
+ };
+
+public:
+ template <class RadiusFunction>
+ SegmentNetwork(const GridGeometry& gg, const RadiusFunction& radiusFunction)
+ {
+ const auto& gridView = gg.gridView();
+ segments_.resize(gridView.size(0));
+ for (const auto &element : elements(gridView))
+ {
+ const auto geometry = element.geometry();
+ const auto eIdx = gg.elementMapper().index(element);
+ const auto radius = radiusFunction(eIdx);
+
+ segments_[eIdx] = Segment{ geometry.corner(0), geometry.corner(1), radius, eIdx };
+ }
+
+ std::cout << "-- Extracted " << segments_.size() << " segments.\n";
+ }
+
+ //! the segment radius
+ auto radius(std::size_t eIdx) const
+ { return segments_[eIdx].r; }
+
+ //! the segment with index eIdx
+ const auto& segment(std::size_t eIdx) const
+ { return segments_[eIdx]; }
+
+ //! the segments
+ const auto& segments() const
+ { return segments_; }
+
+ //! find the closest segment (the segment from which the signed distance to its surface is the closest to zero)
+ //! \note minDist is a positive distance
+ auto findClosestSegmentSurface(const GlobalPosition& globalPos) const
+ {
+ GridIndex eIdx = 0;
+ double minSignedDist = std::numeric_limits<double>::max();
+
+ for (const auto &segment : segments_)
+ {
+ const auto signedDist = computeSignedDistance_(globalPos, segment);
+ if (signedDist < minSignedDist)
+ {
+ minSignedDist = signedDist;
+ eIdx = segment.eIdx;
+ }
+ }
+
+ using std::abs;
+ return std::make_tuple(eIdx, abs(minSignedDist));
+ }
+
+ //! same as overload with taking a position but only look at the segments specified by the index vector
+ template <class IndexRange>
+ auto findClosestSegmentSurface(const GlobalPosition& globalPos,
+ const IndexRange& segIndices) const
+ {
+ GridIndex eIdx = 0;
+ double minSignedDist = std::numeric_limits<double>::max();
+
+ for (const auto index : segIndices)
+ {
+ const auto &segment = segments_[index];
+ const auto signedDist = computeSignedDistance_(globalPos, segment);
+ if (signedDist < minSignedDist)
+ {
+ minSignedDist = signedDist;
+ eIdx = segment.eIdx;
+ }
+ }
+
+ using std::abs;
+ return std::make_tuple(eIdx, abs(minSignedDist));
+ }
+
+ // Compute the projection point of p onto the segment connecting a->b
+ GlobalPosition projectionPointOnSegment(const GlobalPosition& p, std::size_t idx) const
+ {
+ const auto& segment = segments_[idx];
+ return projectionPointOnSegment_(p, segment.a, segment.b);
+ }
+
+private:
+ // Compute the projection of p onto the segment
+ // returns the closest end point if the orthogonal projection onto the line implied by the segment
+ // is outside the segment
+ GlobalPosition projectionPointOnSegment_(const GlobalPosition& p, const GlobalPosition& a, const GlobalPosition& b) const
+ {
+ const auto v = b - a;
+ const auto w = p - a;
+
+ const auto proj1 = v*w;
+ if (proj1 <= 0.0)
+ return a;
+
+ const auto proj2 = v.two_norm2();
+ if (proj2 <= proj1)
+ return b;
+
+ const auto t = proj1 / proj2;
+ auto x = a;
+ x.axpy(t, v);
+ return x;
+ }
+
+ // Compute the signed (!) distance to a cylinder segment surface and the projection point onto the centerline
+ auto computeSignedDistance_(const GlobalPosition& p, const Segment& segment) const
+ {
+ const auto projPoint = projectionPointOnSegment_(p, segment.a, segment.b);
+ return (p-projPoint).two_norm()-segment.r;
+ }
+
+ std::vector<Segment> segments_;
+};
+
+/*!
+ * \ingroup EmbeddedCoupling
+ * \brief Get the closest segment for a given surface point
+ * \note The point does not necessarily be exactly on the virtual surface
+ * implied by the networks radius function
+ */
+template<class Network>
+class NetworkIndicatorFunction
+{
+public:
+ NetworkIndicatorFunction(const Network& n)
+ : network_(n)
+ {}
+
+ template<class Pos>
+ auto operator()(const Pos& pos) const
+ {
+ const auto& [closestSegmentIdx, _] = network_.findClosestSegmentSurface(pos);
+ return closestSegmentIdx;
+ }
+
+ template <class Pos, class HintContainer>
+ auto operator()(const Pos& pos, const HintContainer& hints) const
+ {
+ const auto& [closestSegmentIdx, _] = network_.findClosestSegmentSurface(pos, hints);
+ return closestSegmentIdx;
+ }
+
+private:
+ const Network& network_;
+};
+
+/*!
+ * \ingroup EmbeddedCoupling
+ * \brief Simple legacy VTK writer for outputting debug data on the coupling interface
+ */
+class DebugIntersectionVTKOutput
+{
+public:
+ void push(const std::array<Dune::FieldVector<double, 3>, 3>& corners, double data)
+ {
+ vtkCells_.emplace_back(std::array<std::size_t, 3>{
+ { vtkPoints_.size(), vtkPoints_.size()+1, vtkPoints_.size()+2 }
+ });
+ for (const auto& c : corners)
+ vtkPoints_.push_back(c);
+ vtkCellData_.push_back(data);
+ }
+
+ void write(const std::string& name)
+ {
+ std::ofstream intersectionFile(name);
+ intersectionFile << "# vtk DataFile Version 2.0\n";
+ intersectionFile << "Really cool intersection data\n";
+ intersectionFile << "ASCII\n";
+ intersectionFile << "DATASET UNSTRUCTURED_GRID\n";
+ intersectionFile << "POINTS " << vtkPoints_.size() << " float\n";
+ for (const auto& p : vtkPoints_)
+ intersectionFile << p[0] << " " << p[1] << " " << p[2] << "\n";
+ intersectionFile << "\nCELLS " << vtkCells_.size() << " " << vtkCells_.size()*4 << "\n";
+ for (const auto& c : vtkCells_)
+ intersectionFile << "3 " << c[0] << " " << c[1] << " " << c[2] << "\n";
+ intersectionFile << "\nCELL_TYPES " << vtkCells_.size() << "\n";
+ for (int i = 0; i < vtkCells_.size(); ++i)
+ intersectionFile << "5\n";
+ intersectionFile << "\nCELL_DATA " << vtkCells_.size() << "\nSCALARS index float\nLOOKUP_TABLE default\n";
+ for (const auto& c : vtkCellData_)
+ intersectionFile << c << "\n";
+ }
+
+private:
+ std::vector<Dune::FieldVector<double, 3>> vtkPoints_;
+ std::vector<std::array<std::size_t, 3>> vtkCells_;
+ std::vector<double> vtkCellData_;
+};
+
+} // end namespace Detail
+
+namespace Embedded1d3dCouplingMode {
+struct Projection : public Utility::Tag<Projection> {
+ static std::string name() { return "projection"; }
+};
+
+inline constexpr Projection projection{};
+} // end namespace Embedded1d3dCouplingMode
+
+// forward declaration
+template<class MDTraits, class CouplingMode>
+class Embedded1d3dCouplingManager;
+
+/*!
+ * \ingroup EmbeddedCoupling
+ * \brief Manages the coupling between bulk elements and lower dimensional elements
+ *
+ * The method is described in detail in the article
+ * "Projection-based resolved interface mixed-dimension method for embedded tubular network systems"
+ * by Timo Koch (2021) available at https://arxiv.org/abs/2106.06358
+ *
+ * The network is represented as a line segment network. The bulk domain explicitly resolves the wall
+ * of the tube network (e.g. blood vessel wall, outer root wall), so this generally requires an unstructured grid.
+ * The coupling term coupling the PDEs on the network and in the bulk domain is intergrated over the
+ * coupling surface and each integration point couples to quantities evaluated at the closest point on the graph.
+ * There is a unique mapping from every point on the virtual tube surface to the tube centerline, therefore
+ * this coupling is determined by mapping to the closest point on the virtual surface and then evaluating
+ * the 1D quantity on the mapped centerline point. (The inverse mapping may be non-unique.)
+ *
+ * The original paper also includes an algorithm for computing exact intersection for the surface quadrature
+ * in case of simple straight vessels. However, as a comparison in the paper shows that the suggested
+ * approximate quadrature is exact enough (configured by parameter MixedDimension.SimplexIntegrationRefine
+ * specifying the number of (adaptive local) virtual refinements of a surface facet for integration), only
+ * the approximate general purpose algorithm is included in this implementation. Only one quadrature point
+ * will be added for each surface element and coupled 1D dof including all contribution evaluated by
+ * local refinement. That means the virtual refinement doesn't increase the number of integration points and
+ * additionally is also optimized by using local adaptive refinement only in the places where the index mapping
+ * to the 1D degree of freedom is still multivalent.
+ *
+ * Coupling sources are implement in terms of point sources at quadrature points to make it easy to reuse code
+ * written for other 1D-3D coupling manager modes.
+ *
+ * This algorithm can be configured by several parameters
+ * (although this is usually not necessary as there are sensible defaults):
+ *
+ * - MixedDimension.SimplexIntegrationRefine
+ * number of virtual refinement steps to determine coupled surface area
+ * - MixedDimension.EnableIntersectionOutput
+ * set to true to enable debug VTK output for intersections
+ * - MixedDimension.EstimateNumberOfPointSources
+ * provide an estimate for the expected number of coupling points for memory allocation
+ * - MixedDimension.ProjectionCoupledRadiusFactor
+ * threshold distance in which to search for coupled elements (specified as multiple of radius, default 0.1)
+ * - MixedDimension.ProjectionCoupledAngleFactor
+ * angle threshold in which to search for coupled elements (angle in radians from surface normal vector, default 0.3)
+ *
+ */
+template<class MDTraits>
+class Embedded1d3dCouplingManager<MDTraits, Embedded1d3dCouplingMode::Projection>
+: public EmbeddedCouplingManagerBase<MDTraits, Embedded1d3dCouplingManager<MDTraits, Embedded1d3dCouplingMode::Projection>>
+{
+ using ThisType = Embedded1d3dCouplingManager<MDTraits, Embedded1d3dCouplingMode::Projection>;
+ using ParentType = EmbeddedCouplingManagerBase<MDTraits, ThisType>;
+
+ using Scalar = typename MDTraits::Scalar;
+ using SolutionVector = typename MDTraits::SolutionVector;
+ using PointSourceData = typename ParentType::PointSourceTraits::PointSourceData;
+
+ static constexpr auto bulkIdx = typename MDTraits::template SubDomain<0>::Index();
+ static constexpr auto lowDimIdx = typename MDTraits::template SubDomain<1>::Index();
+
+ // the sub domain type aliases
+ template<std::size_t id> using SubDomainTypeTag = typename MDTraits::template SubDomain<id>::TypeTag;
+ template<std::size_t id> using Problem = GetPropType<SubDomainTypeTag<id>, Properties::Problem>;
+ template<std::size_t id> using GridGeometry = GetPropType<SubDomainTypeTag<id>, Properties::GridGeometry>;
+ template<std::size_t id> using GridView = typename GridGeometry<id>::GridView;
+ template<std::size_t id> using Element = typename GridView<id>::template Codim<0>::Entity;
+ template<std::size_t id> using GridIndex = typename IndexTraits<GridView<id>>::GridIndex;
+
+ static constexpr int lowDimDim = GridView<lowDimIdx>::dimension;
+ static constexpr int bulkDim = GridView<bulkIdx>::dimension;
+ static constexpr int dimWorld = GridView<bulkIdx>::dimensionworld;
+
+ template<std::size_t id>
+ static constexpr bool isBox()
+ { return GridGeometry<id>::discMethod == DiscretizationMethod::box; }
+
+ using GlobalPosition = typename Element<bulkIdx>::Geometry::GlobalCoordinate;
+
+ using SegmentNetwork = Detail::SegmentNetwork<GridGeometry<lowDimIdx>>;
+
+public:
+ static constexpr Embedded1d3dCouplingMode::Projection couplingMode{};
+
+ using ParentType::ParentType;
+
+ void init(std::shared_ptr<Problem<bulkIdx>> bulkProblem,
+ std::shared_ptr<Problem<lowDimIdx>> lowDimProblem,
+ const SolutionVector& curSol)
+ {
+ ParentType::init(bulkProblem, lowDimProblem, curSol);
+ }
+
+ /* \brief Compute integration point point sources and associated data
+ * This method computes a source term on the explicitly given surface of the network grid
+ * by integrating over the elements of the surface, projecting the integration points onto
+ * the centerlines to evaluate 1d quantities and set point sources (minimum distance projection)
+ *
+ * \param order Unused in this algorithm! (passed from default 1D3D coupling manager)
+ * \param verbose If the source computation is verbose
+ */
+ void computePointSourceData(std::size_t order = 1, bool verbose = false)
+ {
+ // initialize the maps
+ // do some logging and profiling
+ Dune::Timer watch;
+ std::cout << "Initializing the coupling manager (projection)" << std::endl;
+
+ // debug VTK output
+ const bool enableIntersectionOutput = getParam<bool>("MixedDimension.EnableIntersectionOutput", false);
+ std::unique_ptr<Detail::DebugIntersectionVTKOutput> vtkCoupledFaces =
+ enableIntersectionOutput ? std::make_unique<Detail::DebugIntersectionVTKOutput>() : nullptr;
+ std::unique_ptr<Detail::DebugIntersectionVTKOutput> vtkIntersections =
+ enableIntersectionOutput ? std::make_unique<Detail::DebugIntersectionVTKOutput>() : nullptr;
+
+ // temporarily represent the network domain as a list of segments
+ // a more efficient data structure
+ const auto& lowDimProblem = this->problem(lowDimIdx);
+ const auto& lowDimFvGridGeometry = lowDimProblem.gridGeometry();
+ const auto& lowDimGridView = lowDimFvGridGeometry.gridView();
+ localAreaFactor_.resize(lowDimGridView.size(0), 0.0);
+ const auto radiusFunc = [&](const GridIndex<lowDimIdx> eIdx) {
+ return lowDimProblem.spatialParams().radius(eIdx);
+ };
+ const auto network = SegmentNetwork{
+ lowDimFvGridGeometry, radiusFunc
+ };
+
+ // precompute the vertex indices for efficiency (box method)
+ this->precomputeVertexIndices(bulkIdx);
+ this->precomputeVertexIndices(lowDimIdx);
+
+ const auto& bulkProblem = this->problem(bulkIdx);
+ const auto& bulkFvGridGeometry = bulkProblem.gridGeometry();
+ const auto& bulkGridView = bulkFvGridGeometry.gridView();
+ bulkElementMarker_.assign(bulkGridView.size(0), 0);
+ bulkVertexMarker_.assign(bulkGridView.size(bulkDim), 0);
+ const GlobalPosition threshold(1e-2*(bulkFvGridGeometry.bBoxMax()-bulkFvGridGeometry.bBoxMin()).two_norm());
+ const auto bBoxMinSmall = bulkFvGridGeometry.bBoxMin() + threshold;
+ const auto bBoxMaxSmall = bulkFvGridGeometry.bBoxMax() - threshold;
+ auto insideBBox = [bBoxMin=bBoxMinSmall, bBoxMax=bBoxMaxSmall](const GlobalPosition& point) -> bool
+ {
+ static constexpr Scalar eps_ = 1.0e-7;
+ const Scalar eps0 = eps_*(bBoxMax[0] - bBoxMin[0]);
+ const Scalar eps1 = eps_*(bBoxMax[1] - bBoxMin[1]);
+ const Scalar eps2 = eps_*(bBoxMax[2] - bBoxMin[2]);
+ return (bBoxMin[0] - eps0 <= point[0] && point[0] <= bBoxMax[0] + eps0 &&
+ bBoxMin[1] - eps1 <= point[1] && point[1] <= bBoxMax[1] + eps1 &&
+ bBoxMin[2] - eps2 <= point[2] && point[2] <= bBoxMax[2] + eps2);
+ };
+
+ // setup simplex "quadrature" rule
+ static const auto quadSimplexRefineMaxLevel
+ = getParam<std::size_t>("MixedDimension.SimplexIntegrationRefine", 4);
+ const auto indicator = Detail::NetworkIndicatorFunction { network };
+
+ // estimate number of point sources for memory allocation
+ static const auto estimateNumberOfPointSources
+ = getParam<std::size_t>("MixedDimension.EstimateNumberOfPointSources", bulkFvGridGeometry.gridView().size(0));
+
+ this->pointSources(bulkIdx).reserve(estimateNumberOfPointSources);
+ this->pointSources(lowDimIdx).reserve(estimateNumberOfPointSources);
+ this->pointSourceData().reserve(estimateNumberOfPointSources);
+
+ // we go over the bulk surface intersection, check if we are coupled,
+ // and if this is the case, we add integration point sources using
+ // the local simplex refinement quadrature procedure (see paper Koch 2021)
+ for (const auto& element : elements(bulkGridView))
+ {
+ const auto bulkElementIdx = bulkFvGridGeometry.elementMapper().index(element);
+ const auto bulkGeometry = element.geometry();
+ const auto refElement = Dune::referenceElement(element);
+ for (const auto& intersection : intersections(bulkGridView, element))
+ {
+ // skip inner intersection
+ if (!intersection.boundary())
+ continue;
+
+ // check if we are close enough to the coupling interface
+ const auto [isCoupled, closestSegmentIdx, minDist] = isCoupled_(network, intersection, insideBBox);
+ if (!isCoupled)
+ continue;
+
+ // we found an intersection on the coupling interface: mark as coupled element
+ bulkElementMarker_[bulkElementIdx] = 1;
+
+ const auto interfaceGeometry = intersection.geometry();
+ for (int i = 0; i < interfaceGeometry.corners(); ++i)
+ {
+ const auto localVIdx = refElement.subEntity(intersection.indexInInside(), 1, i, bulkDim);
+ const auto vIdx = bulkFvGridGeometry.vertexMapper().subIndex(element, localVIdx, bulkDim);
+ bulkVertexMarker_[vIdx] = 1;
+ }
+
+ // debug graphical intersection output
+ if (enableIntersectionOutput)
+ vtkCoupledFaces->push({
+ { interfaceGeometry.corner(0), interfaceGeometry.corner(1), interfaceGeometry.corner(2)}
+ }, closestSegmentIdx);
+
+ const auto quadSimplex = LocalRefinementSimplexQuadrature{ interfaceGeometry, quadSimplexRefineMaxLevel, indicator };
+
+ // iterate over all quadrature points (children simplices)
+ for (const auto& qp : quadSimplex)
+ {
+ const auto surfacePoint = interfaceGeometry.global(qp.position());
+ const auto closestSegmentIdx = qp.indicator();
+ const auto projPoint = network.projectionPointOnSegment(surfacePoint, closestSegmentIdx);
+
+ addPointSourceAtIP_(bulkGeometry, interfaceGeometry, qp, bulkElementIdx, closestSegmentIdx, surfacePoint, projPoint);
+ addStencilEntries_(bulkElementIdx, closestSegmentIdx);
+ }
+ }
+ }
+
+ // make stencils unique
+ using namespace Dune::Hybrid;
+ forEach(integralRange(Dune::index_constant<2>{}), [&](const auto domainIdx)
+ {
+ for (auto&& stencil : this->couplingStencils(domainIdx))
+ {
+ std::sort(stencil.second.begin(), stencil.second.end());
+ stencil.second.erase(
+ std::unique(stencil.second.begin(), stencil.second.end()),
+ stencil.second.end()
+ );
+ }
+ });
+
+ // debug VTK output
+ if (enableIntersectionOutput)
+ vtkCoupledFaces->write("coupledfaces.vtk");
+
+ // compare theoretical cylinder surface to actual discrete surface
+ // of the coupled bulk interface facets
+ for (const auto& element : elements(lowDimGridView))
+ {
+ const auto length = element.geometry().volume();
+ const auto eIdx = lowDimFvGridGeometry.elementMapper().index(element);
+ const auto radius = this->problem(lowDimIdx).spatialParams().radius(eIdx);
+ const auto cylinderSurface = 2.0*M_PI*radius*length;
+ localAreaFactor_[eIdx] /= cylinderSurface;
+ localAreaFactor_[eIdx] -= 1.0;
+ }
+
+ this->pointSources(bulkIdx).shrink_to_fit();
+ this->pointSources(lowDimIdx).shrink_to_fit();
+ this->pointSourceData().shrink_to_fit();
+
+ std::cout << "-- Coupling at " << this->pointSourceData().size() << " integration points." << std::endl;
+ std::cout << "-- Finished initialization after " << watch.elapsed() << " seconds." << std::endl;
+ }
+
+ /*!
+ * \brief Methods to be accessed by the subproblems
+ */
+ // \{
+
+ //! Return a reference to the bulk problem
+ Scalar radius(std::size_t id) const
+ {
+ const auto& data = this->pointSourceData()[id];
+ return this->problem(lowDimIdx).spatialParams().radius(data.lowDimElementIdx());
+ }
+
+ // \}
+
+ //! Marker that is non-zero for bulk elements that are coupled
+ const std::vector<int>& bulkElementMarker() const
+ { return bulkElementMarker_; }
+
+ //! Marker that is non-zero for bulk vertices that belong to coupled surface facets
+ const std::vector<int>& bulkVertexMarker() const
+ { return bulkVertexMarker_; }
+
+ //! For each 1D element, the ratio of associated discrete bulk surface area
+ //! to the theoretical bulk surface area resulting from assuming a cylinder-shaped segment
+ //! (excluding cylinder caps).
+ const std::vector<Scalar>& localAreaFactor() const
+ { return localAreaFactor_; }
+
+private:
+ std::vector<int> bulkElementMarker_, bulkVertexMarker_;
+ std::vector<double> localAreaFactor_;
+
+ // add a point source containing all data needed to evaluate
+ // the coupling source term from both domains
+ template<class BulkGeometry, class InterfaceGeometry, class QP>
+ void addPointSourceAtIP_(const BulkGeometry& bulkGeometry,
+ const InterfaceGeometry& ifGeometry,
+ const QP& qp,
+ const GridIndex<bulkIdx> bulkElementIdx,
+ const GridIndex<lowDimIdx> lowDimElementIdx,
+ const GlobalPosition& surfacePoint,
+ const GlobalPosition& projPoint)
+ {
+ // add a new point source id (for the associated data)
+ const auto id = this->idCounter_++;
+
+ // add a bulk point source
+ const auto qpweight = qp.weight();
+ const auto integrationElement = ifGeometry.integrationElement(qp.position());
+ this->pointSources(bulkIdx).emplace_back(surfacePoint, id, qpweight, integrationElement, bulkElementIdx);
+
+ // find the corresponding projection point and 1d element
+ this->pointSources(lowDimIdx).emplace_back(projPoint, id, qpweight, integrationElement, lowDimElementIdx);
+ localAreaFactor_[lowDimElementIdx] += qpweight*integrationElement;
+
+ const auto& bulkFvGridGeometry = this->problem(bulkIdx).gridGeometry();
+ const auto& lowDimFvGridGeometry = this->problem(lowDimIdx).gridGeometry();
+
+ // pre compute additional data used for the evaluation of
+ // the actual solution dependent source term
+ PointSourceData psData;
+
+ if constexpr (isBox<lowDimIdx>())
+ {
+ std::vector<Dune::FieldVector<Scalar, 1> > shapeValues;
+ const auto lowDimGeometry = this->problem(lowDimIdx).gridGeometry().element(lowDimElementIdx).geometry();
+ this->getShapeValues(lowDimIdx, lowDimFvGridGeometry, lowDimGeometry, projPoint, shapeValues);
+ psData.addLowDimInterpolation(shapeValues, this->vertexIndices(lowDimIdx, lowDimElementIdx), lowDimElementIdx);
+ }
+ else
+ {
+ psData.addLowDimInterpolation(lowDimElementIdx);
+ }
+
+ if constexpr (isBox<bulkIdx>())
+ {
+ std::vector<Dune::FieldVector<Scalar, 1> > shapeValues;
+ this->getShapeValues(bulkIdx, bulkFvGridGeometry, bulkGeometry, surfacePoint, shapeValues);
+ psData.addBulkInterpolation(shapeValues, this->vertexIndices(bulkIdx, bulkElementIdx), bulkElementIdx);
+ }
+ else
+ {
+ psData.addBulkInterpolation(bulkElementIdx);
+ }
+
+ // publish point source data in the global vector
+ this->pointSourceData().emplace_back(std::move(psData));
+ }
+
+ void addStencilEntries_(const GridIndex<bulkIdx> bulkElementIdx, const GridIndex<lowDimIdx> lowDimElementIdx)
+ {
+ // export the bulk coupling stencil
+ if constexpr (isBox<lowDimIdx>())
+ {
+ const auto& vertices = this->vertexIndices(lowDimIdx, lowDimElementIdx);
+ this->couplingStencils(bulkIdx)[bulkElementIdx].insert(
+ this->couplingStencils(bulkIdx)[bulkElementIdx].end(),
+ vertices.begin(), vertices.end()
+ );
+
+ }
+ else
+ {
+ this->couplingStencils(bulkIdx)[bulkElementIdx].push_back(lowDimElementIdx);
+ }
+
+ // export the lowdim coupling stencil
+ if constexpr (isBox<bulkIdx>())
+ {
+ const auto& vertices = this->vertexIndices(bulkIdx, bulkElementIdx);
+ this->couplingStencils(lowDimIdx)[lowDimElementIdx].insert(
+ this->couplingStencils(lowDimIdx)[lowDimElementIdx].end(),
+ vertices.begin(), vertices.end()
+ );
+ }
+ else
+ {
+ this->couplingStencils(lowDimIdx)[lowDimElementIdx].push_back(bulkElementIdx);
+ }
+ }
+
+ template<class BoundaryIntersection, class BBoxCheck>
+ auto isCoupled_(const SegmentNetwork& network, const BoundaryIntersection& intersection, const BBoxCheck& insideBBox) const
+ {
+ const auto center = intersection.geometry().center();
+ const auto [eIdx, minDist] = network.findClosestSegmentSurface(center);
+
+ // all inner elements are coupled
+ if (insideBBox(center))
+ return std::make_tuple(true, eIdx, minDist);
+
+ // if the distance is under a certain threshold we are coupled (and the angle is not very big)
+ static const auto rFactor = getParam<Scalar>("MixedDimension.ProjectionCoupledRadiusFactor", 0.1);
+ static const auto spFactor = getParam<Scalar>("MixedDimension.ProjectionCoupledAngleFactor", 0.3);
+ const auto projPoint = network.projectionPointOnSegment(center, eIdx);
+ const auto distVec = projPoint - center;
+ const bool isCoupled = minDist < rFactor * network.radius(eIdx) && distVec * intersection.centerUnitOuterNormal() > distVec.two_norm() * spFactor;
+ return std::make_tuple(isCoupled, eIdx, minDist);
+ }
+};
+
+} // end namespace Dumux
+
+#endif
diff --git a/dumux/multidomain/embedded/localrefinementquadrature.hh b/dumux/multidomain/embedded/localrefinementquadrature.hh
new file mode 100644
index 0000000000000000000000000000000000000000..0b4ea574f674ba6c761544fc008fe3d8faa1fb42
--- /dev/null
+++ b/dumux/multidomain/embedded/localrefinementquadrature.hh
@@ -0,0 +1,184 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ * See the file COPYING for full copying permissions. *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 3 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+/*!
+ * \file
+ * \ingroup EmbeddedCoupling
+ * \brief A quadrature rule for surface coupling in the 1D-3D setting
+ */
+
+#ifndef DUMUX_MULTIDOMAIN_EMBEDDED_LOCAL_REFINEMENT_QUADRATURE_HH
+#define DUMUX_MULTIDOMAIN_EMBEDDED_LOCAL_REFINEMENT_QUADRATURE_HH
+
+#include <array>
+#include <vector>
+#include <utility>
+#include <algorithm>
+
+namespace Dumux {
+
+template<class Geometry, class IndicatorFunction>
+class LocalRefinementSimplexQuadrature
+{
+ using IndicatorTriple = std::array<std::size_t, 3>;
+ using GlobalPosition = typename Geometry::GlobalCoordinate;
+ using LocalPosition = typename Geometry::LocalCoordinate;
+ using Corners = std::array<GlobalPosition, 3>;
+
+ class QuadPoint
+ {
+ LocalPosition localPos_;
+ double weight_;
+ std::size_t indicator_;
+
+ public:
+ QuadPoint(LocalPosition&& localPos, double weight, std::size_t i)
+ : localPos_(weight*std::move(localPos))
+ , weight_(weight)
+ , indicator_(i)
+ {}
+
+ void add(LocalPosition&& localPos, double weight)
+ {
+ localPos_ += weight*localPos;
+ weight_ += weight;
+ }
+
+ void finalize()
+ {
+ localPos_ /= weight_;
+ }
+
+ const LocalPosition& position() const { return localPos_; }
+ double weight() const { return weight_; }
+ std::size_t indicator() const { return indicator_; }
+ };
+
+ class Triangle
+ {
+ Corners corners_;
+ public:
+ Triangle(const Corners& c) : corners_(c) {}
+ Triangle(Corners&& c) : corners_(std::move(c)) {}
+ const GlobalPosition& operator[](std::size_t i) const { return corners_[i]; }
+
+ GlobalPosition center() const
+ {
+ GlobalPosition center(0.0);
+ for (const auto& c : corners_)
+ center += c;
+ center /= corners_.size();
+ return center;
+ }
+
+ auto volume() const
+ {
+ const auto ab = corners_[1] - corners_[0];
+ const auto ac = corners_[2] - corners_[0];
+ return crossProduct(ab, ac).two_norm()*0.5;
+ }
+ };
+
+public:
+ LocalRefinementSimplexQuadrature(const Geometry& geo, std::size_t maxLevel, const IndicatorFunction& i)
+ : ind_(i)
+ , maxLevel_(maxLevel)
+ , geometry_(geo)
+ {
+ // evaluate indicator
+ const auto tri = Triangle{{ geo.corner(0), geo.corner(1), geo.corner(2) }};
+ const auto triple = IndicatorTriple{{ ind_(tri[0]), ind_(tri[1]), ind_(tri[2]) }};
+ volume_ = tri.volume();
+
+ // add new sub-qps recursively by adaptive refinement
+ addSimplex_(tri, 0, triple, triple);
+
+ for (auto& qp : qps_)
+ qp.finalize();
+ }
+
+ auto begin() const { return qps_.begin(); }
+ auto end() const { return qps_.end(); }
+ auto size() const { return qps_.size(); }
+
+private:
+ void addSimplex_(const Triangle& tri, std::size_t level, const IndicatorTriple& triple, const IndicatorTriple& childTriple)
+ {
+ // if all indicator are the same just add the triangle
+ if (std::all_of(childTriple.begin(), childTriple.end(), [a0=childTriple[0]](auto a){ return (a == a0); }))
+ {
+ // the volumes need to sum up to 0.5 (triangle reference element volume)
+ auto it = std::find_if(qps_.begin(), qps_.end(), [&](const auto& qp){ return (qp.indicator() == childTriple[0]); });
+ if (it != qps_.end())
+ it->add(geometry_.local(tri.center()), 0.5*tri.volume()/volume_);
+ else
+ qps_.emplace_back(geometry_.local(tri.center()), 0.5*tri.volume()/volume_, childTriple[0]);
+ }
+
+ // if they are different but the maximum level is reached, split volume in three
+ else if (level == maxLevel_)
+ {
+ for (int i = 0; i < 3; ++i)
+ {
+ // the volumes need to sum up to 0.5 (triangle reference element volume)
+ auto it = std::find_if(qps_.begin(), qps_.end(), [&](const auto& qp){ return (qp.indicator() == childTriple[i]); });
+ if (it != qps_.end())
+ it->add(geometry_.local(tri[i]), 0.5*tri.volume()/volume_/3.0);
+ else
+ qps_.emplace_back(geometry_.local(tri[i]), 0.5*tri.volume()/volume_/3.0, childTriple[i]);
+ }
+ }
+
+ // otherwise refine and go recursively over all children
+ else
+ {
+ const auto children = refine(tri);
+ for (const auto& c : children)
+ {
+ // make sure to recompute indicator every time since new indices might come up
+ const auto grandChildTriple = IndicatorTriple{{ ind_(c[0]), ind_(c[1]), ind_(c[2]) }};
+ addSimplex_(c, level+1, triple, grandChildTriple);
+ }
+ }
+ }
+
+ std::array<Triangle, 4> refine(const Triangle& tri) const
+ {
+ const auto p01 = 0.5*(tri[0] + tri[1]);
+ const auto p12 = 0.5*(tri[1] + tri[2]);
+ const auto p20 = 0.5*(tri[2] + tri[0]);
+
+ return {{
+ {{ tri[0], p01, p20 }},
+ {{ tri[1], p01, p12 }},
+ {{ tri[2], p12, p20 }},
+ {{ p12, p01, p20 }}
+ }};
+ }
+
+ const IndicatorFunction &ind_;
+ std::size_t maxLevel_;
+ const Geometry& geometry_;
+ double volume_;
+
+ std::vector<QuadPoint> qps_;
+};
+
+} // end namespace Dumux
+
+#endif