[multidomain][embedded] Add coupling managers for 1d3d and 2d3d embedded coupling

dumux/multidomain/CMakeLists.txt

+add_subdirectory(embedded)
+
 install(FILES
 couplingjacobianpattern.hh
 couplingmanager.hh

dumux/multidomain/embedded/CMakeLists.txt

+install(FILES
+circlepoints.hh
+couplingmanager1d3d.hh
+couplingmanager2d3d.hh
+couplingmanagerbase.hh
+extendedsourcestencil.hh
+integrationpointsource.hh
+mixeddimensionglue.hh
+pointsourcedata.hh
+DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dumux/multidomain/embedded)

dumux/multidomain/embedded/circlepoints.hh

+// -*- 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 2 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 Helper function to compute points on a circle
+ */
+
+#ifndef DUMUX_MULTIDOMAIN_EMBEDDED_CIRCLEPOINTS_HH
+#define DUMUX_MULTIDOMAIN_EMBEDDED_CIRCLEPOINTS_HH
+
+#include <vector>
+#include <cmath>
+
+namespace Dumux {
+namespace EmbeddedCoupling {
+
+/*!
+ * \ingroup EmbeddedCoupling
+ * \brief returns a vector of points on a circle
+ * \param center the circle center
+ * \param normal the normal to the circle plane
+ * \param radius the circle radius
+ * \param numPoints the number of points
+ */
+template<class GlobalPosition>
+std::vector<GlobalPosition> circlePoints(const GlobalPosition& center,
+                                         const GlobalPosition& normal,
+                                         const typename GlobalPosition::value_type radius,
+                                         const std::size_t numPoints = 20)
+{
+    using std::abs; using std::sin; using std::cos;
+    using ctype = typename GlobalPosition::value_type;
+
+    constexpr ctype eps = 1.5e-7;
+    static_assert(GlobalPosition::dimension == 3, "Only implemented for world dimension 3");
+
+    std::vector<GlobalPosition> points(numPoints);
+
+    // make sure n is a unit vector
+    auto n = normal;
+    n /= n.two_norm();
+
+    // caculate a vector u perpendicular to n
+    GlobalPosition u;
+    if (abs(n[0]) < eps && abs(n[1]) < eps)
+        if (abs(n[2]) < eps)
+            DUNE_THROW(Dune::MathError, "The normal vector has to be non-zero!");
+        else
+            u = {0, 1, 0};
+    else
+        u = {-n[1], n[0], 0};
+
+    u *= radius/u.two_norm();
+
+    // the circle parameterization is p(t) = r*cos(t)*u + r*sin(t)*(n x u) + c
+    auto tangent = crossProduct(u, n);
+    tangent *= radius/tangent.two_norm();
+
+    // the parameter with an offset
+    ctype t = 0 + 0.1;
+    // insert the vertices
+    for (std::size_t i = 0; i < numPoints; ++i)
+    {
+        points[i] = GlobalPosition({u[0]*cos(t) + tangent[0]*sin(t) + center[0],
+                                    u[1]*cos(t) + tangent[1]*sin(t) + center[1],
+                                    u[2]*cos(t) + tangent[2]*sin(t) + center[2]});
+
+        t += 2*M_PI/numPoints;
+
+        // periodic t
+        if(t > 2*M_PI)
+            t -= 2*M_PI;
+    }
+
+    return points;
+}
+
+} // end namespace EmbeddedCoupling
+} // end namespace Dumux
+
+#endif

dumux/multidomain/embedded/couplingmanager2d3d.hh

+// -*- 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 2 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 MultiDomain
+ * \ingroup EmbeddedCoupling
+ * \brief Coupling manager for embedded fractures
+ */
+
+#ifndef DUMUX_MULTIDOMAIN_EMBEDDED_COUPLINGMANAGER_2D3D_HH
+#define DUMUX_MULTIDOMAIN_EMBEDDED_COUPLINGMANAGER_2D3D_HH
+
+#include <dumux/multidomain/embedded/couplingmanagerbase.hh>
+
+namespace Dumux {
+
+/*!
+ * \ingroup MultiDomain
+ * \ingroup EmbeddedCoupling
+ * \brief Coupling manager for embedded fractures
+ * \note we just use the default coupling manager
+ */
+template<class MDTraits>
+class EmbeddedCouplingManager2d3d
+: public EmbeddedCouplingManagerBase<MDTraits,
+                                     EmbeddedCouplingManager2d3d<MDTraits>>
+{
+    using ThisType = EmbeddedCouplingManager2d3d<MDTraits>;
+    using ParentType = EmbeddedCouplingManagerBase<MDTraits, ThisType>;
+public:
+    using ParentType::ParentType;
+};
+
+} // end namespace Dumux
+
+#endif

dumux/multidomain/embedded/extendedsourcestencil.hh

+// -*- 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 2 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 Extended source stencil helper class for coupling managers
+ */
+
+#ifndef DUMUX_MULTIDOMAIN_EMBEDDED_EXTENDEDSOURCESTENCIL_HH
+#define DUMUX_MULTIDOMAIN_EMBEDDED_EXTENDEDSOURCESTENCIL_HH
+
+#include <vector>
+#include <dumux/common/properties.hh>
+
+namespace Dumux {
+namespace EmbeddedCoupling {
+
+/*!
+ * \ingroup EmbeddedCoupling
+ * \brief A class managing an extended source stencil
+ * \tparam CouplingManager the coupling manager type
+ */
+template<class CouplingManager>
+class ExtendedSourceStencil
+{
+    using MDTraits = typename CouplingManager::MultiDomainTraits;
+    using Scalar = typename MDTraits::Scalar;
+
+    template<std::size_t id> using SubDomainTypeTag = typename MDTraits::template SubDomainTypeTag<id>;
+    template<std::size_t id> using FVGridGeometry = typename GET_PROP_TYPE(SubDomainTypeTag<id>, FVGridGeometry);
+    template<std::size_t id> using GridView = typename FVGridGeometry<id>::GridView;
+    template<std::size_t id> using Element = typename GridView<id>::template Codim<0>::Entity;
+
+    static constexpr auto bulkIdx = typename MDTraits::template DomainIdx<0>();
+    static constexpr auto lowDimIdx = typename MDTraits::template DomainIdx<1>();
+
+    template<std::size_t id>
+    static constexpr bool isBox()
+    { return FVGridGeometry<id>::discMethod == DiscretizationMethod::box; }
+public:
+    /*!
+     * \brief extend the jacobian pattern of the diagonal block of domain i
+     *        by those entries that are not already in the uncoupled pattern
+     * \note Such additional dependencies can arise from the coupling, e.g. if a coupling source
+     *       term depends on a non-local average of a quantity of the same domain
+     */
+    template<std::size_t id, class JacobianPattern>
+    void extendJacobianPattern(const CouplingManager& couplingManager, Dune::index_constant<id> domainI, JacobianPattern& pattern) const
+    {
+        // add additional dof dependencies
+        for (const auto& element : elements(couplingManager.gridView(domainI)))
+        {
+            const auto& dofs = extendedSourceStencil_(couplingManager, domainI, element);
+
+            if (isBox<domainI>())
+            {
+                for (int i = 0; i < element.subEntities(GridView<domainI>::dimension); ++i)
+                    for (const auto globalJ : dofs)
+                        pattern.add(couplingManager.problem(domainI).fvGridGeometry().vertexMapper().subIndex(element, i, GridView<domainI>::dimension), globalJ);
+            }
+            else
+            {
+                const auto globalI = couplingManager.problem(domainI).fvGridGeometry().elementMapper().index(element);
+                for (const auto globalJ : dofs)
+                    pattern.add(globalI, globalJ);
+            }
+        }
+    }
+
+    /*!
+     * \brief evaluate additional derivatives of the element residual of a domain with respect
+     *        to dofs in the same domain that are not in the regular stencil (per default this is not the case)
+     * \note Such additional dependencies can arise from the coupling, e.g. if a coupling source
+     *       term depends on a non-local average of a quantity of the same domain
+     * \note This is the same for box and cc
+     */
+    template<std::size_t i, class LocalAssemblerI, class SolutionVector, class JacobianMatrixDiagBlock, class GridVariables>
+    void evalAdditionalDomainDerivatives(CouplingManager& couplingManager,
+                                         Dune::index_constant<i> domainI,
+                                         const LocalAssemblerI& localAssemblerI,
+                                         const SolutionVector& curSol,
+                                         JacobianMatrixDiagBlock& A,
+                                         GridVariables& gridVariables) const
+    {
+        constexpr auto numEq = std::decay_t<decltype(curSol[domainI][0])>::dimension;
+        const auto& elementI = localAssemblerI.element();
+
+        // only do something if we have an extended stencil
+        if (extendedSourceStencil_(couplingManager, domainI, elementI).empty())
+            return;
+
+        // compute the undeflected residual (source only!)
+        const auto origResidual = localAssemblerI.evalLocalSourceResidual(elementI);
+
+        // compute derivate for all additional dofs in the circle stencil
+        for (const auto dofIndex : extendedSourceStencil_(couplingManager, domainI, elementI))
+        {
+            auto partialDerivs = origResidual;
+            const auto origPriVars = curSol[domainI][dofIndex];
+
+            // calculate derivatives w.r.t to the privars at the dof at hand
+            for (int pvIdx = 0; pvIdx < numEq; pvIdx++)
+            {
+                // reset partial derivatives
+                partialDerivs = 0.0;
+
+                auto evalResiduals = [&](Scalar priVar)
+                {
+                    // update the coupling context (solution vector and recompute element residual)
+                    auto priVars = origPriVars;
+                    priVars[pvIdx] = priVar;
+                    couplingManager.updateCouplingContext(domainI, localAssemblerI, domainI, dofIndex, priVars, pvIdx);
+                    return localAssemblerI.evalLocalSourceResidual(elementI);
+                };
+
+                // derive the residuals numerically
+                static const int numDiffMethod = getParam<int>("Assembly.NumericDifferenceMethod");
+                NumericDifferentiation::partialDerivative(evalResiduals, curSol[domainI][dofIndex][pvIdx],
+                                                          partialDerivs, origResidual, numDiffMethod);
+
+                // update the global stiffness matrix with the current partial derivatives
+                for (auto&& scvJ : scvs(localAssemblerI.fvGeometry()))
+                {
+                    for (int eqIdx = 0; eqIdx < numEq; eqIdx++)
+                    {
+                        // A[i][col][eqIdx][pvIdx] is the rate of change of
+                        // the residual of equation 'eqIdx' at dof 'i'
+                        // depending on the primary variable 'pvIdx' at dof
+                        // 'col'.
+                        A[scvJ.dofIndex()][dofIndex][eqIdx][pvIdx] += partialDerivs[scvJ.indexInElement()][eqIdx];
+                    }
+                }
+
+                // restore the original coupling context
+                couplingManager.updateCouplingContext(domainI, localAssemblerI, domainI, dofIndex, origPriVars, pvIdx);
+            }
+        }
+    }
+
+    //! return a reference to the stencil
+    typename CouplingManager::CouplingStencils& stencil()
+    { return sourceStencils_; }
+
+private:
+    //! the extended source stencil for the bulk domain due to the source average
+    const std::vector<std::size_t>& extendedSourceStencil_(const CouplingManager& couplingManager, Dune::index_constant<0> bulkDomain, const Element<0>& bulkElement) const
+    {
+        const auto bulkElementIdx = couplingManager.problem(bulkIdx).fvGridGeometry().elementMapper().index(bulkElement);
+        if (sourceStencils_.count(bulkElementIdx))
+            return sourceStencils_.at(bulkElementIdx);
+        else
+            return couplingManager.emptyStencil();
+    }
+
+    //! the extended source stencil for the low dim domain is empty
+    const std::vector<std::size_t>& extendedSourceStencil_(const CouplingManager& couplingManager, Dune::index_constant<1> bulkDomain, const Element<1>& lowDimElement) const
+    { return couplingManager.emptyStencil(); }
+
+    //! the additional stencil for the kernel evaluations / circle averages
+    typename CouplingManager::CouplingStencils sourceStencils_;
+};
+
+} // end namespace EmbeddedCoupling
+} // end namespace Dumux
+
+#endif

dumux/multidomain/embedded/integrationpointsource.hh

+// -*- 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 2 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 An integration point source class,
+ *        i.e. sources located at a single point in space
+ *        associated with a quadrature point
+ */
+
+#ifndef DUMUX_INTEGRATION_POINTSOURCE_HH
+#define DUMUX_INTEGRATION_POINTSOURCE_HH
+
+#include <type_traits>
+#include <dune/common/reservedvector.hh>
+#include <dumux/common/pointsource.hh>
+#include <dumux/discretization/methods.hh>
+
+namespace Dumux {
+
+/*!
+ * \ingroup EmbeddedCoupling
+ * \brief An integration point source class with an identifier to attach data
+ *        and a quadrature weight and integration element
+ */
+template<class GlobalPosition, class SourceValues, typename IdType = std::size_t>
+class IntegrationPointSource : public IdPointSource<GlobalPosition, SourceValues, IdType>
+{
+    using ParentType = IdPointSource<GlobalPosition, SourceValues, IdType>;
+    using Scalar = std::decay_t<decltype(std::declval<SourceValues>()[0])>;
+
+public:
+    //! Constructor for integration point sources
+    IntegrationPointSource(GlobalPosition pos, SourceValues values, IdType id,
+                           Scalar qpweight, Scalar integrationElement,
+                           const std::vector<std::size_t>& elementIndices)
+      : ParentType(pos, values, id),
+        qpweight_(qpweight), integrationElement_(integrationElement),
+        elementIndices_(elementIndices) {}
+
+    //! Constructor for integration point sources, when there is no
+    // value known at the time of initialization
+    IntegrationPointSource(GlobalPosition pos, IdType id,
+                           Scalar qpweight, Scalar integrationElement,
+                           const std::vector<std::size_t>& elementIndices)
+      : ParentType(pos, id),
+        qpweight_(qpweight), integrationElement_(integrationElement),
+        elementIndices_(elementIndices) {}
+
+
+    Scalar quadratureWeight() const
+    {
+        return qpweight_;
+    }
+
+    Scalar integrationElement() const
+    {
+        return integrationElement_;
+    }
+
+    const std::vector<std::size_t>& elementIndices() const
+    {
+        return elementIndices_;
+    }
+
+    //! Convenience = operator overload modifying only the values
+    IntegrationPointSource& operator= (const SourceValues& values)
+    {
+        ParentType::operator=(values);
+        return *this;
+    }
+
+    //! Convenience = operator overload modifying only the values
+    IntegrationPointSource& operator= (Scalar s)
+    {
+        ParentType::operator=(s);
+        return *this;
+    }
+
+private:
+    Scalar qpweight_;
+    Scalar integrationElement_;
+    std::vector<std::size_t> elementIndices_;
+};
+
+/*!
+ * \ingroup EmbeddedCoupling
+ * \brief A helper class calculating a DOF-index to point source map
+ */
+class IntegrationPointSourceHelper
+{
+
+public:
+    //! calculate a DOF index to point source map from given vector of point sources
+    template<class FVGridGeometry, class PointSource, class PointSourceMap>
+    static void computePointSourceMap(const FVGridGeometry& fvGridGeometry,
+                                      std::vector<PointSource>& sources,
+                                      PointSourceMap& pointSourceMap)
+    {
+        for (auto&& source : sources)
+        {
+            // compute in which elements the point source falls
+            const auto& entities = source.elementIndices();
+            // split the source values equally among all concerned entities
+            source.setEmbeddings(source.embeddings()*entities.size());
+            // loop over all intersected elements
+            for (unsigned int eIdx : entities)
+            {
+                if (FVGridGeometry::discMethod == DiscretizationMethod::box)
+                {
+                    // check in which subcontrolvolume(s) we are
+                    const auto element = fvGridGeometry.boundingBoxTree().entitySet().entity(eIdx);
+                    auto fvGeometry = localView(fvGridGeometry);
+                    fvGeometry.bindElement(element);
+                    const auto globalPos = source.position();
+                    // loop over all sub control volumes and check if the point source is inside
+                    constexpr int dim = FVGridGeometry::GridView::dimension;
+                    Dune::ReservedVector<std::size_t, 1<<dim> scvIndices;
+                    for (auto&& scv : scvs(fvGeometry))
+                        if (intersectsPointGeometry(globalPos, scv.geometry()))
+                            scvIndices.push_back(scv.indexInElement());
+
+                    // for all scvs that where tested positiv add the point sources
+                    // to the element/scv to point source map
+                    for (auto scvIdx : scvIndices)
+                    {
+                        const auto key = std::make_pair(eIdx, scvIdx);
+                        if (pointSourceMap.count(key))
+                            pointSourceMap.at(key).push_back(source);
+                        else
+                            pointSourceMap.insert({key, {source}});
+                        // split equally on the number of matched scvs
+                        auto& s = pointSourceMap.at(key).back();
+                        s.setEmbeddings(scvIndices.size()*s.embeddings());
+                    }
+                }
+                else
+                {
+                    // add the pointsource to the DOF map
+                    const auto key = std::make_pair(eIdx, /*scvIdx=*/ 0);
+                    if (pointSourceMap.count(key))
+                        pointSourceMap.at(key).push_back(source);
+                    else
+                        pointSourceMap.insert({key, {source}});
+                }
+            }
+        }
+    }
+};
+
+} // end namespace Dumux
+
+#endif

dumux/multidomain/embedded/mixeddimensionglue.hh

+// -*- 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 2 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 MultiDomain
+ * \ingroup MixedDimension
+ * \brief A class glueing two grids of different dimension geometrically
+ *        Intersections are computed using axis-aligned bounding box trees
+ */
+
+#ifndef DUMUX_MULTIDOMAIN_MIXEDDIMENSION_GLUE_HH
+#define DUMUX_MULTIDOMAIN_MIXEDDIMENSION_GLUE_HH
+
+#include <iostream>
+#include <fstream>
+#include <string>
+#include <utility>
+
+#include <dune/common/timer.hh>
+#include <dune/common/iteratorrange.hh>
+#include <dune/geometry/affinegeometry.hh>
+#include <dune/geometry/type.hh>
+
+#include <dumux/common/geometry/geometricentityset.hh>
+#include <dumux/common/geometry/boundingboxtree.hh>
+#include <dumux/common/geometry/intersectingentities.hh>
+
+namespace Dumux {
+
+// forward declaration
+template<class BulkGridView, class LowDimGridView, class BulkMapper, class LowDimMapper>
+class MixedDimensionGlue;
+
+/*!
+ * \ingroup MixedDimension
+ * \brief Range generator to iterate with range-based for loops over all intersections
+ *        as follows: for (const auto& is : intersections(glue)) { ... }
+ */
+template<class BulkGridView, class LowDimGridView, class BulkMapper, class LowDimMapper>
+Dune::IteratorRange<typename MixedDimensionGlue<BulkGridView, LowDimGridView, BulkMapper, LowDimMapper>::Intersections::const_iterator>
+intersections(const MixedDimensionGlue<BulkGridView, LowDimGridView, BulkMapper, LowDimMapper>& glue)
+{ return {glue.ibegin(), glue.iend()}; }
+
+namespace Glue {
+
+/*!
+ * \ingroup MixedDimension
+ * \brief An intersection object representing an intersection
+ *        between two grid entites of different grids
+ */
+template<class BulkGridView, class LowDimGridView, class BulkMapper, class LowDimMapper>