diff --git a/dumux/multidomain/CMakeLists.txt b/dumux/multidomain/CMakeLists.txt
index 366ff8b3000d71bbcdeee3175133d135028d4458..a9a4acd7ab6e3f5cda61873f120a89b11c80b3da 100644
--- a/dumux/multidomain/CMakeLists.txt
+++ b/dumux/multidomain/CMakeLists.txt
@@ -1,3 +1,4 @@
+add_subdirectory(boundary)
add_subdirectory(embedded)
install(FILES
diff --git a/dumux/multidomain/boundary/CMakeLists.txt b/dumux/multidomain/boundary/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..51f43731d8cd80d2701a5a520170f01de1c9d8b8
--- /dev/null
+++ b/dumux/multidomain/boundary/CMakeLists.txt
@@ -0,0 +1 @@
+add_subdirectory(darcydarcy)
diff --git a/dumux/multidomain/boundary/darcydarcy/CMakeLists.txt b/dumux/multidomain/boundary/darcydarcy/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..38ab77687c1c8c644aab5d337dab618e09f8b4ff
--- /dev/null
+++ b/dumux/multidomain/boundary/darcydarcy/CMakeLists.txt
@@ -0,0 +1,4 @@
+install(FILES
+couplingmanager.hh
+couplingmapper.hh
+DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dumux/multidomain/boundary/darcydarcy)
diff --git a/dumux/multidomain/boundary/darcydarcy/couplingmanager.hh b/dumux/multidomain/boundary/darcydarcy/couplingmanager.hh
new file mode 100644
index 0000000000000000000000000000000000000000..7e8f39cf43ed4c06b320d2a5fa479ae719ea2b9a
--- /dev/null
+++ b/dumux/multidomain/boundary/darcydarcy/couplingmanager.hh
@@ -0,0 +1,253 @@
+// -*- 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 BoundaryCoupling
+ * \brief Coupling manager for equal-dimension boundary coupling
+ */
+
+#ifndef DUMUX_MULTIDOMAIN_DARCYDARCY_BOUNDARY_COUPLINGMANAGER_HH
+#define DUMUX_MULTIDOMAIN_DARCYDARCY_BOUNDARY_COUPLINGMANAGER_HH
+
+#include <iostream>
+#include <vector>
+
+#include <dune/common/indices.hh>
+#include <dune/common/exceptions.hh>
+
+#include <dumux/common/properties.hh>
+#include <dumux/common/math.hh>
+#include <dumux/discretization/elementsolution.hh>
+#include <dumux/discretization/methods.hh>
+#include <dumux/discretization/cellcentered/tpfa/computetransmissibility.hh>
+#include <dumux/multidomain/couplingmanager.hh>
+#include <dumux/multidomain/boundary/darcydarcy/couplingmapper.hh>
+
+namespace Dumux {
+
+/*!
+ * \ingroup MultiDomain
+ * \ingroup BoundaryCoupling
+ * \brief Coupling manager for equal-dimension boundary coupling of darcy models
+ */
+template<class MDTraits>
+class DarcyDarcyBoundaryCouplingManager
+: public CouplingManager<MDTraits>
+{
+ using ParentType = CouplingManager<MDTraits>;
+
+ using Scalar = typename MDTraits::Scalar;
+ using SolutionVector = typename MDTraits::SolutionVector;
+
+ template<std::size_t i> using SubDomainTypeTag = typename MDTraits::template SubDomainTypeTag<i>;
+ template<std::size_t i> using Problem = typename GET_PROP_TYPE(SubDomainTypeTag<i>, Problem);
+ template<std::size_t i> using PrimaryVariables = typename GET_PROP_TYPE(SubDomainTypeTag<i>, PrimaryVariables);
+ template<std::size_t i> using NumEqVector = typename GET_PROP_TYPE(SubDomainTypeTag<i>, NumEqVector);
+ template<std::size_t i> using ElementVolumeVariables = typename GET_PROP_TYPE(SubDomainTypeTag<i>, GridVolumeVariables)::LocalView;
+ template<std::size_t i> using VolumeVariables = typename GET_PROP_TYPE(SubDomainTypeTag<i>, GridVolumeVariables)::VolumeVariables;
+ template<std::size_t i> using FVGridGeometry = typename MDTraits::template SubDomainFVGridGeometry<i>;
+ template<std::size_t i> using FVElementGeometry = typename FVGridGeometry<i>::LocalView;
+ template<std::size_t i> using SubControlVolumeFace = typename FVGridGeometry<i>::SubControlVolumeFace;
+ template<std::size_t i> using SubControlVolume = typename FVGridGeometry<i>::SubControlVolume;
+ template<std::size_t i> using GridView = typename FVGridGeometry<i>::GridView;
+ template<std::size_t i> using Element = typename GridView<i>::template Codim<0>::Entity;
+
+ template<std::size_t i>
+ static constexpr auto domainIdx()
+ { return typename MDTraits::template DomainIdx<i>{}; }
+
+ template<std::size_t i>
+ static constexpr bool isCCTpfa()
+ { return FVGridGeometry<i>::discMethod == DiscretizationMethod::cctpfa; }
+
+ using CouplingStencil = std::vector<std::size_t>;
+public:
+ //! export traits
+ using MultiDomainTraits = MDTraits;
+ //! export stencil types
+ using CouplingStencils = std::unordered_map<std::size_t, CouplingStencil>;
+
+ /*!
+ * \brief Methods to be accessed by main
+ */
+ // \{
+
+ void init(std::shared_ptr<Problem<domainIdx<0>()>> problem0,
+ std::shared_ptr<Problem<domainIdx<1>()>> problem1,
+ const SolutionVector& curSol)
+ {
+ static_assert(isCCTpfa<0>() && isCCTpfa<1>(), "Only cctpfa implemented!");
+
+ this->updateSolution(curSol);
+ problemTuple_ = std::make_tuple(problem0, problem1);
+ couplingMapper_.update(*this);
+ }
+
+ // \}
+
+ /*!
+ * \brief Methods to be accessed by the assembly
+ */
+ // \{
+
+ /*!
+ * \brief returns an iteratable container of all indices of degrees of freedom of domain j
+ * that couple with / influence the element residual of the given element of domain i
+ *
+ * \param domainI the domain index of domain i
+ * \param elementI the coupled element of domain Ã
+ * \param domainJ the domain index of domain j
+ */
+ template<std::size_t i, std::size_t j>
+ const CouplingStencil& couplingStencil(Dune::index_constant<i> domainI,
+ const Element<i>& element,
+ Dune::index_constant<j> domainJ) const
+ {
+ static_assert(i != j, "A domain cannot be coupled to itself!");
+ const auto eIdx = problem(domainI).fvGridGeometry().elementMapper().index(element);
+ return couplingMapper_.couplingStencil(domainI, eIdx, domainJ);
+ }
+
+ // \}
+
+ /*!
+ * \brief If the boundary entity is on a coupling boundary
+ * \param domainI the domain index for which to compute the flux
+ * \param scvf the sub control volume face
+ */
+ template<std::size_t i>
+ bool isCoupled(Dune::index_constant<i> domainI,
+ const SubControlVolumeFace<i>& scvf) const
+ {
+ return couplingMapper_.isCoupled(domainI, scvf);
+ }
+
+ /*!
+ * \brief Evaluate the boundary conditions for a coupled Neumann boundary segment.
+ *
+ * This is the method for the case where the Neumann condition is
+ * potentially solution dependent
+ *
+ * \param domainI the domain index for which to compute the flux
+ * \param element The finite element
+ * \param fvGeometry The finite-volume geometry
+ * \param elemVolVars All volume variables for the element
+ * \param scvf The sub control volume face
+ * \param phaseIdx the phase for which to compute the flux
+ *
+ * Negative values mean influx.
+ * E.g. for the mass balance that would the mass flux in \f$ [ kg / (m^2 \cdot s)] \f$.
+ */
+ template<std::size_t i>
+ Scalar advectiveFluxCoupling(Dune::index_constant<i> domainI,
+ const Element<i>& element,
+ const FVElementGeometry<i>& fvGeometry,
+ const ElementVolumeVariables<i>& elemVolVars,
+ const SubControlVolumeFace<i>& scvf,
+ int phaseIdx = 0) const
+ {
+ Scalar flux = 0.0;
+ static const bool enableGravity = getParamFromGroup<bool>(problem(domainI).paramGroup(), "Problem.EnableGravity");
+ constexpr auto otherDomainIdx = domainIdx<1-i>();
+
+ const auto& outsideElement = problem(otherDomainIdx).fvGridGeometry().element(couplingMapper_.outsideElementIndex(domainI, scvf));
+ auto fvGeometryOutside = localView(problem(otherDomainIdx).fvGridGeometry());
+ fvGeometryOutside.bindElement(outsideElement);
+
+ const auto& flipScvf = fvGeometryOutside.scvf(couplingMapper_.flipScvfIndex(domainI, scvf));
+ const auto& outsideScv = fvGeometryOutside.scv(flipScvf.insideScvIdx());
+ const auto outsideVolVars = volVars(otherDomainIdx, outsideElement, outsideScv);
+
+ const auto& insideScv = fvGeometry.scv(scvf.insideScvIdx());
+ const auto& insideVolVars = elemVolVars[insideScv];
+
+ const auto ti = computeTpfaTransmissibility(scvf, insideScv, insideVolVars.permeability(), insideVolVars.extrusionFactor());
+ const auto tj = computeTpfaTransmissibility(flipScvf, outsideScv, outsideVolVars.permeability(), outsideVolVars.extrusionFactor());
+ Scalar tij = 0.0;
+ if (ti*tj > 0.0)
+ tij = scvf.area()*(ti * tj)/(ti + tj);
+
+ if (enableGravity)
+ {
+ // do averaging for the density over all neighboring elements
+ const auto rho = (insideVolVars.density(phaseIdx) + outsideVolVars.density(phaseIdx))*0.5;
+ const auto& g = problem(domainI).gravityAtPos(scvf.ipGlobal());
+
+ //! compute alpha := n^T*K*g
+ const auto alpha_inside = vtmv(scvf.unitOuterNormal(), insideVolVars.permeability(), g)*insideVolVars.extrusionFactor();
+ const auto alpha_outside = vtmv(scvf.unitOuterNormal(), outsideVolVars.permeability(), g)*outsideVolVars.extrusionFactor();
+
+ // Obtain inside and outside pressures
+ const auto pInside = insideVolVars.pressure(0);
+ const auto pOutside = outsideVolVars.pressure(0);
+
+ flux = tij*(pInside - pOutside) + rho*scvf.area()*alpha_inside - rho*tij/tj*(alpha_inside - alpha_outside);
+
+ }
+ else
+ {
+ // Obtain inside and outside pressures
+ const auto pInside = insideVolVars.pressure(phaseIdx);
+ const auto pOutside = outsideVolVars.pressure(phaseIdx);
+
+ // return flux
+ flux = tij*(pInside - pOutside);
+ }
+
+ // upwind scheme
+ static const Scalar upwindWeight = getParam<Scalar>("Implicit.UpwindWeight");
+ auto upwindTerm = [phaseIdx](const auto& volVars){ return volVars.density(phaseIdx)*volVars.mobility(phaseIdx); };
+ if (std::signbit(flux)) // if sign of flux is negative
+ flux *= (upwindWeight*upwindTerm(outsideVolVars)
+ + (1.0 - upwindWeight)*upwindTerm(insideVolVars));
+ else
+ flux *= (upwindWeight*upwindTerm(insideVolVars)
+ + (1.0 - upwindWeight)*upwindTerm(outsideVolVars));
+
+ // make this a area-specific flux
+ flux /= scvf.area()*insideVolVars.extrusionFactor();
+ return flux;
+ }
+
+ //! Return a reference to the sub problem with domain index i
+ template<std::size_t i>
+ const Problem<i>& problem(Dune::index_constant<i> domainIdx) const
+ { return *std::get<i>(problemTuple_); }
+
+ //! Return the volume variables of domain i for a given element and scv
+ template<std::size_t i>
+ VolumeVariables<i> volVars(Dune::index_constant<i> domainI,
+ const Element<i>& element,
+ const SubControlVolume<i>& scv) const
+ {
+ VolumeVariables<i> volVars;
+ const auto elemSol = elementSolution(element, this->curSol()[domainI], problem(domainI).fvGridGeometry());
+ volVars.update(elemSol, problem(domainI), element, scv);
+ return volVars;
+ }
+
+private:
+ typename MDTraits::ProblemTuple problemTuple_;
+ DarcyDarcyBoundaryCouplingMapper<MDTraits> couplingMapper_;
+};
+
+} // end namespace Dumux
+
+#endif
diff --git a/dumux/multidomain/boundary/darcydarcy/couplingmapper.hh b/dumux/multidomain/boundary/darcydarcy/couplingmapper.hh
new file mode 100644
index 0000000000000000000000000000000000000000..079fbc02f2e5f7dcac0653e64716a6e986e08271
--- /dev/null
+++ b/dumux/multidomain/boundary/darcydarcy/couplingmapper.hh
@@ -0,0 +1,240 @@
+// -*- 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 BoundaryCoupling
+ * \copydoc Dumux::BoundaryCouplingMapper
+ */
+
+#ifndef DUMUX_MULTIDOMAIN_DARCYDARCY_BOUNDARY_COUPLINGMAPPER_HH
+#define DUMUX_MULTIDOMAIN_DARCYDARCY_BOUNDARY_COUPLINGMAPPER_HH
+
+#include <iostream>
+#include <unordered_map>
+#include <tuple>
+#include <vector>
+
+#include <dune/common/timer.hh>
+#include <dune/common/exceptions.hh>
+#include <dune/common/indices.hh>
+
+#include <dumux/common/geometry/intersectingentities.hh>
+#include <dumux/discretization/methods.hh>
+
+namespace Dumux {
+
+/*!
+ * \ingroup MultiDomain
+ * \ingroup BoundaryCoupling
+ * \brief the default mapper for conforming equal dimension boundary coupling between two domains (box or cc)
+ * \todo how to extend to arbitrary number of domains?
+ */
+template<class MDTraits>
+class DarcyDarcyBoundaryCouplingMapper
+{
+ using Scalar = typename MDTraits::Scalar;
+
+ template<std::size_t i> using FVGridGeometry = typename MDTraits::template SubDomainFVGridGeometry<i>;
+ template<std::size_t i> using SubControlVolumeFace = typename FVGridGeometry<i>::SubControlVolumeFace;
+ template<std::size_t i> using GridView = typename FVGridGeometry<i>::GridView;
+ template<std::size_t i> using Element = typename GridView<i>::template Codim<0>::Entity;
+
+ template<std::size_t i>
+ static constexpr auto domainIdx()
+ { return typename MDTraits::template DomainIdx<i>{}; }
+
+ template<std::size_t i>
+ static constexpr bool isCCTpfa()
+ { return FVGridGeometry<i>::discMethod == DiscretizationMethod::cctpfa; }
+
+ struct ScvfInfo
+ {
+ std::size_t eIdxOutside;
+ std::size_t flipScvfIdx;
+ };
+
+ using FlipScvfMapType = std::unordered_map<std::size_t, ScvfInfo>;
+ using MapType = std::unordered_map<std::size_t, std::vector<std::size_t>>;
+
+ static constexpr std::size_t numSD = MDTraits::numSubDomains;
+
+public:
+ /*!
+ * \brief Main update routine
+ */
+ template<class CouplingManager>
+ void update(const CouplingManager& couplingManager)
+ {
+ // TODO: Box and multiple domains
+ static_assert(numSD == 2, "More than two subdomains not implemented!");
+ static_assert(isCCTpfa<0>() && isCCTpfa<1>(), "Only cctpfa implemented!");
+
+ Dune::Timer watch;
+ std::cout << "Initializing the coupling map..." << std::endl;
+
+ for (std::size_t domIdx = 0; domIdx < numSD; ++domIdx)
+ {
+ stencils_[domIdx].clear();
+ scvfInfo_[domIdx].clear();
+ }
+
+ const auto& problem0 = couplingManager.problem(domainIdx<0>());
+ const auto& problem1 = couplingManager.problem(domainIdx<1>());
+ const auto& gg0 = problem0.fvGridGeometry();
+ const auto& gg1 = problem1.fvGridGeometry();
+
+ isCoupledScvf_[0].resize(gg0.numScvf(), false);
+ isCoupledScvf_[1].resize(gg1.numScvf(), false);
+
+ for (const auto& element0 : elements(gg0.gridView()))
+ {
+ auto fvGeometry0 = localView(gg0);
+ fvGeometry0.bindElement(element0);
+
+ for (const auto& scvf0 : scvfs(fvGeometry0))
+ {
+ // skip all non-boundaries
+ if (!scvf0.boundary())
+ continue;
+
+ // get elements intersecting with the scvf center
+ // for robustness add epsilon in unit outer normal direction
+ const auto eps = (scvf0.ipGlobal() - element0.geometry().corner(0)).two_norm()*1e-8;
+ auto globalPos = scvf0.ipGlobal(); globalPos.axpy(eps, scvf0.unitOuterNormal());
+ const auto indices = intersectingEntities(globalPos, gg1.boundingBoxTree());
+
+ // skip if no intersection was found
+ if (indices.empty())
+ continue;
+
+ // sanity check
+ if (indices.size() > 1)
+ DUNE_THROW(Dune::InvalidStateException, "Are you sure your grids is conforming at the boundary?");
+
+ // add the pair to the multimap
+ const auto eIdx0 = gg0.elementMapper().index(element0);
+ const auto eIdx1 = indices[0];
+ stencils_[0][eIdx0].push_back(eIdx1);
+
+ // mark the scvf and find and mark the flip scvf
+ isCoupledScvf_[0][scvf0.index()] = true;
+ const auto& element1 = gg1.element(eIdx1);
+ auto fvGeometry1 = localView(gg1);
+ fvGeometry1.bindElement(element1);
+
+ using std::abs;
+ for (const auto& scvf1 : scvfs(fvGeometry1))
+ if (scvf1.boundary())
+ if (abs(scvf1.unitOuterNormal()*scvf0.unitOuterNormal() + 1) < eps)
+ {
+ isCoupledScvf_[1][scvf1.index()] = true;
+ scvfInfo_[0][scvf0.index()] = ScvfInfo{eIdx1, scvf1.index()};
+ scvfInfo_[1][scvf1.index()] = ScvfInfo{eIdx0, scvf0.index()};
+ }
+ }
+ }
+
+ // create the inverse map for efficient access
+ for (const auto& entry : stencils_[0])
+ for (const auto idx : entry.second)
+ stencils_[1][idx].push_back(entry.first);
+
+ std::cout << "took " << watch.elapsed() << " seconds." << std::endl;
+ }
+
+ /*!
+ * \brief returns an iteratable container of all indices of degrees of freedom of domain j
+ * that couple with / influence the element residual of the given element of domain i
+ *
+ * \param domainI the domain index of domain i
+ * \param elementI the coupled element of domain Ã
+ * \param domainJ the domain index of domain j
+ *
+ * \note The element residual definition depends on the discretization scheme of domain i
+ * box: a container of the residuals of all sub control volumes
+ * cc : the residual of the (sub) control volume
+ * fem: the residual of the element
+ * \note This function has to be implemented by all coupling managers for all combinations of i and j
+ */
+ template<std::size_t i, std::size_t j>
+ const std::vector<std::size_t>& couplingStencil(Dune::index_constant<i> domainI,
+ const std::size_t eIdxI,
+ Dune::index_constant<j> domainJ) const
+ {
+ static_assert(i != j, "A domain cannot be coupled to itself!");
+ if (isCoupledElement(domainI, eIdxI))
+ return stencils_[i].at(eIdxI);
+ else
+ return emptyStencil_;
+ }
+
+ /*!
+ * \brief Return if an element residual with index eIdx of domain i is coupled to domain j
+ */
+ template<std::size_t i>
+ bool isCoupledElement(Dune::index_constant<i>, std::size_t eIdx) const
+ { return static_cast<bool>(stencils_[i].count(eIdx)); }
+
+ /*!
+ * \brief If the boundary entity is on a coupling boundary
+ * \param domainI the domain index for which to compute the flux
+ * \param scvf the sub control volume face
+ */
+ template<std::size_t i>
+ bool isCoupled(Dune::index_constant<i> domainI,
+ const SubControlVolumeFace<i>& scvf) const
+ {
+ return isCoupledScvf_[i].at(scvf.index());
+ }
+
+ /*!
+ * \brief Return the scvf index of the flipped scvf in the other domain
+ * \param domainI the domain index for which to compute the flux
+ * \param scvf the sub control volume face
+ */
+ template<std::size_t i>
+ std::size_t flipScvfIndex(Dune::index_constant<i> domainI,
+ const SubControlVolumeFace<i>& scvf) const
+ {
+ return scvfInfo_[i].at(scvf.index()).flipScvfIdx;
+ }
+
+ /*!
+ * \brief Return the outside element index (the element index of the other domain)
+ * \param domainI the domain index for which to compute the flux
+ * \param scvf the sub control volume face
+ */
+ template<std::size_t i>
+ std::size_t outsideElementIndex(Dune::index_constant<i> domainI,
+ const SubControlVolumeFace<i>& scvf) const
+ {
+ return scvfInfo_[i].at(scvf.index()).eIdxOutside;
+ }
+
+private:
+ std::array<MapType, numSD> stencils_;
+ std::vector<std::size_t> emptyStencil_;
+ std::array<std::vector<bool>, numSD> isCoupledScvf_;
+ std::array<FlipScvfMapType, numSD> scvfInfo_;
+};
+
+} // end namespace Dumux
+
+#endif