From aad87590e25c0411fd4098c4c4cb401bbcac2d79 Mon Sep 17 00:00:00 2001
From: DennisGlaeser <dennis.glaeser@iws.uni-stuttgart.de>
Date: Fri, 10 Feb 2017 18:04:23 +0100
Subject: [PATCH] [mixeddimension] add facet coupling framework
---
dumux/mixeddimension/CMakeLists.txt | 1 +
dumux/mixeddimension/facet/CMakeLists.txt | 6 +
.../facet/gmshdualfacetgridcreator.hh | 463 ++++++++++++++++++
.../mixeddimension/facet/mpfa/CMakeLists.txt | 10 +
.../facet/mpfa/couplingmanager.hh | 304 ++++++++++++
.../facet/mpfa/couplingmapper.hh | 257 ++++++++++
dumux/mixeddimension/facet/mpfa/darcyslaw.hh | 289 +++++++++++
dumux/mixeddimension/facet/mpfa/fickslaw.hh | 334 +++++++++++++
.../facet/mpfa/interactionvolume.hh | 149 ++++++
.../facet/mpfa/interiorboundarydata.hh | 201 ++++++++
dumux/mixeddimension/facet/mpfa/properties.hh | 72 +++
dumux/mixeddimension/facet/start.hh | 203 ++++++++
12 files changed, 2289 insertions(+)
create mode 100644 dumux/mixeddimension/facet/CMakeLists.txt
create mode 100644 dumux/mixeddimension/facet/gmshdualfacetgridcreator.hh
create mode 100644 dumux/mixeddimension/facet/mpfa/CMakeLists.txt
create mode 100644 dumux/mixeddimension/facet/mpfa/couplingmanager.hh
create mode 100644 dumux/mixeddimension/facet/mpfa/couplingmapper.hh
create mode 100644 dumux/mixeddimension/facet/mpfa/darcyslaw.hh
create mode 100644 dumux/mixeddimension/facet/mpfa/fickslaw.hh
create mode 100644 dumux/mixeddimension/facet/mpfa/interactionvolume.hh
create mode 100644 dumux/mixeddimension/facet/mpfa/interiorboundarydata.hh
create mode 100644 dumux/mixeddimension/facet/mpfa/properties.hh
create mode 100644 dumux/mixeddimension/facet/start.hh
diff --git a/dumux/mixeddimension/CMakeLists.txt b/dumux/mixeddimension/CMakeLists.txt
index a592084f9a..0eca574b98 100644
--- a/dumux/mixeddimension/CMakeLists.txt
+++ b/dumux/mixeddimension/CMakeLists.txt
@@ -1,4 +1,5 @@
add_subdirectory("embedded")
+add_subdirectory("facet")
add_subdirectory("glue")
install(FILES
diff --git a/dumux/mixeddimension/facet/CMakeLists.txt b/dumux/mixeddimension/facet/CMakeLists.txt
new file mode 100644
index 0000000000..a272c20bcc
--- /dev/null
+++ b/dumux/mixeddimension/facet/CMakeLists.txt
@@ -0,0 +1,6 @@
+add_subdirectory("mpfa")
+
+install(FILES
+gmshdualfacetgridcreator.hh
+start.hh
+DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dumux/mixeddimension/facet)
\ No newline at end of file
diff --git a/dumux/mixeddimension/facet/gmshdualfacetgridcreator.hh b/dumux/mixeddimension/facet/gmshdualfacetgridcreator.hh
new file mode 100644
index 0000000000..28d1433781
--- /dev/null
+++ b/dumux/mixeddimension/facet/gmshdualfacetgridcreator.hh
@@ -0,0 +1,463 @@
+/*****************************************************************************
+ * 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
+ *
+ * \brief Reads gmsh files where fractures are defined as lines in surfaces
+ * or surfaces in volumes and constructs a grid for the bulk part and a
+ * lower-dimensional grid for the lines/surfaces.
+ */
+
+#ifndef DUMUX_GMSHDUALFACETGRIDCREATOR_HH
+#define DUMUX_GMSHDUALFACETGRIDCREATOR_HH
+
+#include <algorithm>
+#include <iostream>
+#include <fstream>
+#include <sstream>
+#include <iomanip>
+
+#include <dune/common/version.hh>
+#include <dune/grid/common/mcmgmapper.hh>
+#include <dune/grid/common/gridfactory.hh>
+
+#include <dumux/common/math.hh>
+#include <dumux/common/propertysystem.hh>
+#include <dumux/common/parameters.hh>
+#include <dumux/common/valgrind.hh>
+
+namespace Dumux
+{
+
+namespace Properties
+{
+NEW_PROP_TAG(Scalar);
+}
+
+/*!
+ * \brief Reads gmsh files where fractures are defined as lines in surfaces
+ * or surfaces in volumes and constructs a grid for the bulk part and a
+ * lower-dimensional grid for the lines/surfaces.
+ */
+template <class TypeTag>
+class GmshDualFacetGridCreator
+{
+ using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+
+ // obtain the type tags of the sub problems
+ using BulkProblemTypeTag = typename GET_PROP_TYPE(TypeTag, BulkProblemTypeTag);
+ using LowDimProblemTypeTag = typename GET_PROP_TYPE(TypeTag, LowDimProblemTypeTag);
+
+ // obtain the grid types of the sub problems
+ using BulkGrid = typename GET_PROP_TYPE(BulkProblemTypeTag, Grid);
+ using LowDimGrid = typename GET_PROP_TYPE(LowDimProblemTypeTag, Grid);
+ using LowDimElement = typename LowDimGrid::template Codim<0>::Entity;
+
+ // The Bulk Element Mapper
+ using BulkElementMapper = typename GET_PROP_TYPE(BulkProblemTypeTag, ElementMapper);
+
+ static constexpr int bulkDim = BulkGrid::dimension;
+ static constexpr int lowDimDim = LowDimGrid::dimension;
+ static constexpr int dimWorld = BulkGrid::dimensionworld;
+ static constexpr int lowDimDimWorld = LowDimGrid::dimensionworld;
+ static_assert(dimWorld == lowDimDimWorld, "dimWorld cannot be different for the two sub-domains!");
+
+ static constexpr int invalidIndex = -5;
+
+ using GlobalPosition = Dune::FieldVector<double, dimWorld>;
+ using BulkIndexType = typename BulkGrid::LeafGridView::IndexSet::IndexType;
+ using LowDimIndexType = typename LowDimGrid::LeafGridView::IndexSet::IndexType;
+
+public:
+ /*!
+ * \brief Create the Grid.
+ * Specialization for bulk dim = 2 and lowDim dim = 1
+ */
+ template<class T = TypeTag>
+ static typename std::enable_if<bulkDim == 2 && lowDimDim == 1>::type
+ makeGrid()
+ {
+ const std::string fileName = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, std::string, Grid, File);
+
+ // set up the grid factories
+ auto& bulkFactory = bulkFactory_();
+ auto& lowDimFactory = lowDimFactory_();
+
+ // open the mesh file
+ std::cout << "Opening " << fileName << std::endl;
+ std::ifstream mshFile(fileName.c_str());
+ if (mshFile.fail())
+ DUNE_THROW(Dune::InvalidStateException, "Could not open the given .msh file. Make sure it exists");
+ if (getFileExtension_(fileName) != "msh")
+ DUNE_THROW(Dune::InvalidStateException, "Please provide a .msh file for this grid creator");
+
+ // read file until we get to the list of nodes
+ std::string line;
+ std::getline(mshFile, line);
+ while (line.compare(0, 6, "$Nodes") != 0)
+ std::getline(mshFile, line);
+
+ // get total number of nodes
+ std::getline(mshFile, line);
+ const auto numVertices = stringToNumber_(line);
+
+ // read vertices line by line
+ std::getline(mshFile, line);
+ std::vector<GlobalPosition> bulkVertices(numVertices);
+ unsigned int vertexCounter = 0;
+ while (line.compare(0, 9, "$EndNodes") != 0)
+ {
+ // read the coordinates of the vertex
+ std::string buf;
+ std::stringstream stream(line);
+
+ // drop first entry in line
+ stream >> buf;
+ std::vector<std::string> coords;
+ while (stream >> buf)
+ coords.push_back(buf);
+
+ // "create" a vertex
+ GlobalPosition v;
+ for (int i = 0; i < dimWorld; ++i)
+ v[i] = stringToNumber_(coords[i]);
+
+ // insert vertex into container and the bulk grid factory
+ bulkVertices[vertexCounter] = v;
+ bulkFactory.insertVertex(v);
+ std::getline(mshFile, line);
+ vertexCounter++;
+ }
+
+ // we should always find as many vertices as the mesh file states
+ assert(vertexCounter == numVertices);
+
+ // read file until we get to the list of elements
+ while(line.compare(0, 9, "$Elements") != 0)
+ std::getline(mshFile, line);
+
+ // get total number of elements
+ std::getline(mshFile, line);
+ const auto numElements = stringToNumber_(line);
+
+ // while inserting the low dim vertices, keep track of indices that have
+ // been inserted already and element vertex indices in "original" indexing
+ std::vector<GlobalPosition> lowDimVertices;
+ std::vector<LowDimIndexType> lowDimVertexIndices;
+ std::vector<std::vector<BulkIndexType>> lowDimElemVertices;
+
+ lowDimVertices.reserve(numVertices);
+ lowDimVertexIndices.reserve(numVertices);
+ lowDimElemVertices.reserve(numElements);
+
+ // the coupling data container lowDim -> bulk
+ auto& lowDimCouplingMap = lowDimCouplingElements_();
+
+ // read in the elements line by line
+ std::getline(mshFile, line);
+ BulkIndexType bulkElementCounter = 0;
+ LowDimIndexType lowDimElementCounter = 0;
+ while (line.compare(0, 12, "$EndElements") != 0)
+ {
+ // read the data of this element
+ std::string buf;
+ std::stringstream stream(line);
+ std::vector<unsigned int> elemData;
+ while (stream >> buf)
+ elemData.push_back(stringToNumber_(buf));
+
+ // get number of gmsh tags for this element
+ const auto noOfTags = elemData[2];
+
+ // get element type, make container of indexvertices and insert element
+ const auto elemType = elemData[1];
+
+ switch (elemType)
+ {
+ case 1: // 2-node line
+ {
+ // the geometry type for lines
+ static const Dune::GeometryType gt = [] () { Dune::GeometryType t; t.makeLine(); return t; } ();
+
+ // get the connected vertex indices
+ std::vector<LowDimIndexType> elemVertexIndices(2);
+
+ const auto vIdx1 = elemData[3 + noOfTags];
+ const auto vIdx2 = elemData[3 + noOfTags + 1];
+
+ elemVertexIndices[0] = findIndexInVector_(lowDimVertexIndices, vIdx1);
+ elemVertexIndices[1] = findIndexInVector_(lowDimVertexIndices, vIdx2);
+
+ if (elemVertexIndices[0] == invalidIndex)
+ {
+ const auto v1 = bulkVertices[vIdx1-1];
+
+ lowDimVertices.push_back(v1);
+ lowDimVertexIndices.push_back(vIdx1);
+ lowDimFactory.insertVertex(v1);
+
+ elemVertexIndices[0] = lowDimVertices.size()-1;
+ }
+ if (elemVertexIndices[1] == invalidIndex)
+ {
+ const auto v2 = bulkVertices[vIdx2-1];
+
+ lowDimVertices.push_back(v2);
+ lowDimVertexIndices.push_back(vIdx2);
+ lowDimFactory.insertVertex(v2);
+
+ elemVertexIndices[1] = lowDimVertices.size()-1;
+ }
+
+ // store the low dim elements' vertices
+ lowDimElemVertices.push_back(std::vector<BulkIndexType>({vIdx1, vIdx2}));
+
+ // insert element into the factory
+ lowDimFactory.insertElement(gt, elemVertexIndices);
+
+ // increase low dim element counter
+ lowDimElementCounter++;
+ break;
+ }
+
+ case 2: // 3-node triangle
+ {
+ // the geometry type for triangles
+ static const Dune::GeometryType gt = Dune::GeometryType(Dune::GeometryType::simplex,2);
+
+ // we know that the low dim elements have been read already
+ if (lowDimCouplingMap.empty())
+ lowDimCouplingMap.resize(lowDimElementCounter);
+
+ // get the vertex indices of this bulk element
+ std::vector<BulkIndexType> elemVertexIndices(3);
+ elemVertexIndices[0] = elemData[3 + noOfTags];
+ elemVertexIndices[1] = elemData[3 + noOfTags + 1];
+ elemVertexIndices[2] = elemData[3 + noOfTags + 2];
+
+ // get the insertion indices of the low dim elements connected to this element
+ const auto lowDimElemIndices = obtainLowDimConnections_(lowDimElemVertices, elemVertexIndices);
+
+ // if we found some, insert the connections in the map
+ for (auto lowDimElemIdx : lowDimElemIndices)
+ lowDimCouplingMap[lowDimElemIdx].push_back(bulkElementCounter);
+
+ // subtract 1 from the element indices (gmsh starts at index 1)
+ elemVertexIndices[0] -= 1;
+ elemVertexIndices[1] -= 1;
+ elemVertexIndices[2] -= 1;
+
+ // insert bulk element into the factory
+ bulkFactory.insertElement(gt, elemVertexIndices);
+
+ // increase bulk element counter
+ bulkElementCounter++;
+ break;
+ }
+
+ default:
+ DUNE_THROW(Dune::NotImplemented, "GmshDualFacetGridCreator for given element type");
+ }
+
+ // get next line
+ std::getline(mshFile, line);
+ }
+
+ std::cout << bulkElementCounter << " bulk elements and "
+ << lowDimElementCounter << " low dim elements have been created." << std::endl;
+
+
+ bulkGridPtr_() = std::shared_ptr<BulkGrid>(bulkFactory.createGrid());
+ lowDimGridPtr_() = std::shared_ptr<LowDimGrid>(lowDimFactory.createGrid());
+
+ // set up the map from insertion to actual global indices for bulk elements
+ BulkElementMapper elementMapper(bulkGrid().leafGridView());
+ auto& map = bulkInsertionToGlobalIdxMap_();
+ map.resize(bulkElementCounter);
+ for (const auto& element : elements(bulkGrid().leafGridView()))
+ map[bulkFactory.insertionIndex(element)] = elementMapper.index(element);
+ }
+
+ /*!
+ * \brief Returns a reference to the bulk grid.
+ */
+ static BulkGrid& bulkGrid()
+ { return *bulkGridPtr_(); }
+
+ /*!
+ * \brief Returns a reference to the lower dimensional facet grid.
+ */
+ static LowDimGrid& lowDimGrid()
+ { return *lowDimGridPtr_(); }
+
+ /*!
+ * \brief Distributes the grid on all processes of a parallel
+ * computation.
+ */
+ static void loadBalance()
+ {
+ bulkGrid().loadBalance();
+ lowDimGrid().loadBalance();
+ }
+
+ static std::vector<BulkIndexType> getCoupledBulkElementIndices(const LowDimElement& lowDimElement)
+ {
+ const auto& insertionIndices = lowDimCouplingElements_()[lowDimFactory_().insertionIndex(lowDimElement)];
+ const auto& bulkInsertionToGlobalMap = bulkInsertionToGlobalIdxMap_();
+
+ const auto n = insertionIndices.size();
+ std::vector<BulkIndexType> globalIndices(n);
+ for (unsigned int i = 0; i < n; ++i)
+ globalIndices[i] = bulkInsertionToGlobalMap[insertionIndices[i]];
+
+ return globalIndices;
+ }
+
+private:
+
+ /*!
+ * \brief Returns a reference to the vector containing lists of bulk element indices
+ * that are connected to the low dim elements. Note that the access occurs
+ * with the insertion index of a low dim element.
+ */
+ static std::vector<std::vector<BulkIndexType>>& lowDimCouplingElements_()
+ {
+ static std::vector<std::vector<BulkIndexType>> couplingMap_;
+ return couplingMap_;
+ }
+
+ /*!
+ * \brief Returns a reference to the map mapping bulk element insertion indices
+ * to actual global indices
+ */
+ static std::vector<BulkIndexType>& bulkInsertionToGlobalIdxMap_()
+ {
+ static std::vector<BulkIndexType> insertionIdxMap_;
+ return insertionIdxMap_;
+ }
+
+ /*!
+ * \brief Returns a reference to the bulk grid pointer (std::shared_ptr<BulkGrid>)
+ */
+ static std::shared_ptr<BulkGrid>& bulkGridPtr_()
+ {
+ static std::shared_ptr<BulkGrid> bulkGridPtr_;
+ return bulkGridPtr_;
+ }
+
+ /*!
+ * \brief Returns a reference to the lowDim grid pointer (std::shared_ptr<LowDimGrid>)
+ */
+ static std::shared_ptr<LowDimGrid>& lowDimGridPtr_()
+ {
+ static std::shared_ptr<LowDimGrid> lowDimGridPtr_;
+ return lowDimGridPtr_;
+ }
+
+ static Dune::GridFactory<BulkGrid>& bulkFactory_()
+ {
+ static Dune::GridFactory<BulkGrid> bulkFactory_;
+ return bulkFactory_;
+ }
+
+ static Dune::GridFactory<LowDimGrid>& lowDimFactory_()
+ {
+ static Dune::GridFactory<LowDimGrid> lowDimFactory_;
+ return lowDimFactory_;
+ }
+
+ /*!
+ * \brief Returns a list of low dim elements that are connected to this bulk element
+ */
+ static std::vector<LowDimIndexType> obtainLowDimConnections_(const std::vector<std::vector<LowDimIndexType>>& lowDimElemData,
+ const std::vector<BulkIndexType>& bulkElemVertices)
+ {
+ std::vector<LowDimIndexType> lowDimElemIndices = std::vector<LowDimIndexType>();
+
+ // check for each low dim element if it is contained
+ unsigned int lowDimElementIdx = 0;
+ for (const auto& lowDimElem : lowDimElemData)
+ {
+ const bool isContained = [&lowDimElem, &bulkElemVertices] ()
+ {
+ for (auto idx : lowDimElem)
+ if (findIndexInVector_(bulkElemVertices, idx) == invalidIndex)
+ return false;
+ return true;
+ } ();
+
+ // insert in list if low dim element is in bulk element
+ if (isContained)
+ lowDimElemIndices.push_back(lowDimElementIdx);
+
+ // keep track of the low dim element insertion index
+ lowDimElementIdx++;
+ }
+
+ return lowDimElemIndices;
+ }
+
+ /*!
+ * \brief Returns the local index of a given value in a given vector
+ */
+ template<typename IdxType1, typename IdxType2>
+ static int findIndexInVector_(const std::vector<IdxType1>& vector, const IdxType2 globalIdx)
+ {
+ auto it = std::find(vector.begin(), vector.end(), globalIdx);
+
+ // we use the -5 here to detect whether or not an index exists already
+ if (it != vector.end())
+ return std::distance(vector.begin(), it);
+ else
+ return invalidIndex;
+ }
+
+ /*!
+ * \brief Returns the filename extension of a given filename
+ */
+ static std::string getFileExtension_(const std::string& fileName)
+ {
+ std::size_t i = fileName.rfind('.', fileName.length());
+ if (i != std::string::npos)
+ return(fileName.substr(i+1, fileName.length() - i));
+ else
+ DUNE_THROW(Dune::IOError, "Please provide and extension for your grid file ('"<< fileName << "')!");
+ }
+
+ /*!
+ * \brief Turns a number in form of a string into a Scalar
+ */
+ static Scalar stringToNumber_(const std::string& numberAsString)
+ {
+ Scalar value;
+
+ std::stringstream stream(numberAsString);
+ stream >> value;
+
+ if (stream.fail())
+ {
+ std::runtime_error e(numberAsString);
+ throw e;
+ }
+
+ return value;
+ }
+};
+
+} // end namespace
+
+#endif // DUMUX_GMSHDUALFACETGRIDCREATOR_HH
diff --git a/dumux/mixeddimension/facet/mpfa/CMakeLists.txt b/dumux/mixeddimension/facet/mpfa/CMakeLists.txt
new file mode 100644
index 0000000000..22c0ab9868
--- /dev/null
+++ b/dumux/mixeddimension/facet/mpfa/CMakeLists.txt
@@ -0,0 +1,10 @@
+install(FILES
+couplingmanager.hh
+couplingmapper.hh
+darcyslaw.hh
+fickslaw.hh
+fourierslaw.hh
+interactionvolume.hh
+interiorboundarydata.hh
+properties.hh
+DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dumux/mixeddimension/facet)
\ No newline at end of file
diff --git a/dumux/mixeddimension/facet/mpfa/couplingmanager.hh b/dumux/mixeddimension/facet/mpfa/couplingmanager.hh
new file mode 100644
index 0000000000..6b6d65b0b0
--- /dev/null
+++ b/dumux/mixeddimension/facet/mpfa/couplingmanager.hh
@@ -0,0 +1,304 @@
+/*****************************************************************************
+ * 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 FacetCoupling
+ * \brief Manages the coupling between bulk elements and lower dimensional elements on the element facets
+ */
+
+#ifndef DUMUX_MIXEDDIMENSION_CCMPFA_FACET_COUPLINGMANAGER_HH
+#define DUMUX_MIXEDDIMENSION_CCMPFA_FACET_COUPLINGMANAGER_HH
+
+#include <dumux/implicit/properties.hh>
+#include <dumux/mixeddimension/facet/mpfa/couplingmapper.hh>
+
+namespace Dumux
+{
+
+/*!
+ * \ingroup FacetCoupling
+ * \brief Manages the interaction between the bulk and lower dimensional (facet) elements
+ */
+template<class TypeTag>
+class CCMpfaFacetCouplingManager
+{
+ using Implementation = typename GET_PROP_TYPE(TypeTag, CouplingManager);
+ using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+ using GlobalProblem = typename GET_PROP_TYPE(TypeTag, Problem);
+ using GridCreator = typename GET_PROP_TYPE(TypeTag, GridCreator);
+
+ using BulkProblemTypeTag = typename GET_PROP_TYPE(TypeTag, BulkProblemTypeTag);
+ using LowDimProblemTypeTag = typename GET_PROP_TYPE(TypeTag, LowDimProblemTypeTag);
+
+ using BulkGridView = typename GET_PROP_TYPE(BulkProblemTypeTag, GridView);
+ using BulkElement = typename BulkGridView::template Codim<0>::Entity;
+ using BulkIntersection = typename BulkGridView::Intersection;
+ using BulkIndexType = typename BulkGridView::IndexSet::IndexType;
+ using BulkProblem = typename GET_PROP_TYPE(BulkProblemTypeTag, Problem);
+ using BulkMpfaHelper = typename GET_PROP_TYPE(BulkProblemTypeTag, MpfaHelper);
+ using BulkLocalResidual = typename GET_PROP_TYPE(BulkProblemTypeTag, LocalResidual);
+ using BulkPrimaryVariables = typename GET_PROP_TYPE(BulkProblemTypeTag, PrimaryVariables);
+ using BulkFVElementGeometry = typename GET_PROP_TYPE(BulkProblemTypeTag, FVElementGeometry);
+ using BulkSubControlVolumeFace = typename GET_PROP_TYPE(BulkProblemTypeTag, SubControlVolumeFace);
+ using BulkElementBoundaryTypes = typename GET_PROP_TYPE(BulkProblemTypeTag, ElementBoundaryTypes);
+ using BulkElementVolumeVariables = typename GET_PROP_TYPE(BulkProblemTypeTag, ElementVolumeVariables);
+ using BulkElementFluxVariablesCache = typename GET_PROP_TYPE(BulkProblemTypeTag, ElementFluxVariablesCache);
+
+ using LowDimGridView = typename GET_PROP_TYPE(LowDimProblemTypeTag, GridView);
+ using LowDimElement = typename LowDimGridView::template Codim<0>::Entity;
+ using LowDimIndexType = typename LowDimGridView::IndexSet::IndexType;
+ using LowDimProblem = typename GET_PROP_TYPE(LowDimProblemTypeTag, Problem);
+ using LowDimLocalResidual = typename GET_PROP_TYPE(LowDimProblemTypeTag, LocalResidual);
+ using LowDimPrimaryVariables = typename GET_PROP_TYPE(LowDimProblemTypeTag, PrimaryVariables);
+ using LowDimSubControlVolume = typename GET_PROP_TYPE(LowDimProblemTypeTag, SubControlVolume);
+ using LowDimFVElementGeometry = typename GET_PROP_TYPE(LowDimProblemTypeTag, FVElementGeometry);
+ using LowDimElementBoundaryTypes = typename GET_PROP_TYPE(LowDimProblemTypeTag, ElementBoundaryTypes);
+ using LowDimVolumeVariables = typename GET_PROP_TYPE(LowDimProblemTypeTag, VolumeVariables);
+ using LowDimElementVolumeVariables = typename GET_PROP_TYPE(LowDimProblemTypeTag, ElementVolumeVariables);
+ using LowDimElementFluxVariablesCache = typename GET_PROP_TYPE(LowDimProblemTypeTag, ElementFluxVariablesCache);
+
+ static constexpr int dimWorld = BulkGridView::dimensionworld;
+ static constexpr int lowDimDimWorld = LowDimGridView::dimensionworld;
+ static_assert(dimWorld == lowDimDimWorld, "dimWorld cannot be different for the two sub-domains!");
+
+ using GlobalPosition = Dune::FieldVector<Scalar, dimWorld>;
+
+public:
+
+ /*!
+ * \brief Constructor
+ */
+ CCMpfaFacetCouplingManager(BulkProblem& bulkProblem, LowDimProblem& lowDimProblem)
+ : bulkProblemPtr_(&bulkProblem),
+ lowDimProblemPtr_(&lowDimProblem),
+ couplingMapper_(bulkProblem, lowDimProblem)
+ {
+ // initialize the local residuals
+ bulkLocalResidual_.init(bulkProblem);
+ lowDimLocalResidual_.init(lowDimProblem);
+ }
+
+ void preInit() {}
+
+ void postInit()
+ {
+ // abfrage ob no flow tip??
+ couplingMapper_.init();
+ }
+
+ //! evaluates if an intersection is on an interior boundary
+ bool isInteriorBoundary(const BulkElement& element, const BulkIntersection& is) const
+ {
+ // TODO: Facet elements that are on the bulk domain boundary
+ if (is.boundary())
+ return false;
+
+ const auto insideIdx = bulkProblem().elementMapper().index(element);
+ const auto outsideIdx = bulkProblem().elementMapper().index(is.outside());
+ return isInteriorBoundary_(insideIdx, outsideIdx);
+ }
+
+ //! evaluates if an scvf is on an interior boundary
+ bool isInteriorBoundary(const BulkElement& element, const BulkSubControlVolumeFace& scvf) const
+ {
+ // TODO: Facet elements that are on the bulk domain boundary
+ if (scvf.boundary())
+ return false;
+
+ if (couplingMapper_.isInitialized())
+ {
+ const auto& couplingData = couplingMapper_.getBulkCouplingData(element);
+
+ // if the element is coupled, check if the scvf is so as well
+ if (!couplingData.isCoupled)
+ return false;
+ return couplingData.getScvfCouplingData(scvf).first;
+ }
+
+ // If the mapper has not been initialized yet, use the private function
+ return isInteriorBoundary_(scvf.insideScvIdx(), scvf.outsideScvIdx());
+ }
+
+ const std::vector<LowDimIndexType>& couplingStencil(const BulkElement& element)
+ { return couplingMapper_.getBulkCouplingData(element).couplingStencil; }
+
+ const std::vector<BulkIndexType>& couplingStencil(const LowDimElement& element)
+ { return couplingMapper_.getLowDimCouplingData(element).couplingStencil; }
+
+ //! evaluate coupling residual for the derivative bulk DOF with respect to low dim DOF
+ //! we only need to evaluate the part of the residual that will be influence by the low dim DOF
+ BulkPrimaryVariables evalCouplingResidual(const BulkElement& element,
+ const BulkFVElementGeometry& fvGeometry,
+ const BulkElementVolumeVariables& curElemVolVars,
+ const BulkElementBoundaryTypes& elemBcTypes,
+ const BulkElementFluxVariablesCache& elemFluxVarsCache,
+ const LowDimElement& lowDimElement)
+ {
+ const auto& couplingData = couplingMapper_.getBulkCouplingData(element);
+
+ if (!couplingData.isCoupled)
+ return BulkPrimaryVariables(0.0);
+
+ // calculate the local residual of the bulk element
+ auto prevElemVolVars = localView(bulkProblem().model().prevGlobalVolVars());
+ prevElemVolVars.bindElement(element, fvGeometry, bulkProblem().model().prevSol());
+ bulkLocalResidual_.eval(element, fvGeometry, prevElemVolVars, curElemVolVars, elemBcTypes, elemFluxVarsCache);
+ return bulkLocalResidual_.residual(0);
+ }
+
+ //! evaluate coupling residual for the derivative low dim DOF with respect to bulk DOF
+ //! we only need to evaluate the part of the residual that will be influence by the bulk DOF
+ LowDimPrimaryVariables evalCouplingResidual(const LowDimElement& element,
+ const LowDimFVElementGeometry& fvGeometry,
+ const LowDimElementVolumeVariables& curElemVolVars,
+ const LowDimElementBoundaryTypes& elemBcTypes,
+ const LowDimElementFluxVariablesCache& elemFluxVarsCache,
+ const BulkElement& bulkElement)
+ {
+ const auto& couplingData = couplingMapper_.getLowDimCouplingData(element);
+
+ if (!couplingData.isCoupled)
+ return LowDimPrimaryVariables(0.0);
+
+ // calculate the source term of the low dim element
+ const auto eIdx = lowDimProblem().elementMapper().index(element);
+ const auto& scv = fvGeometry.scv(eIdx);
+
+ auto source = lowDimLocalResidual_.computeSource(element, fvGeometry, curElemVolVars, scv);
+ source *= -scv.volume()*curElemVolVars[scv].extrusionFactor();
+ return source;
+ }
+
+ //! evaluates the sources in the low dim domain coming from the bulk domain
+ //! i.e. the fluxes from the bulk domain into the given low dim element. We return bulk priVars here,
+ //! the low dim problem itself needs to transform that into suitable information
+ BulkPrimaryVariables evalSourcesFromBulk(const LowDimElement& element,
+ const LowDimFVElementGeometry& fvGeometry,
+ const LowDimElementVolumeVariables& curElemVolVars,
+ const LowDimSubControlVolume& scv) const
+ {
+ const auto& couplingData = couplingMapper_.getLowDimCouplingData(element);
+
+ if (!couplingData.isCoupled)
+ return BulkPrimaryVariables(0.0);
+
+ // evaluate the fluxes over each scvf in each bulk neighbor. Prepare the local views
+ // of the first element, the necessary quantities for all the scvfs will be in there,
+ // as the scvfs in the other elements are the "flipped" scvfs of the ones in the first element
+ const auto& firstPair = couplingData.elementScvfList[0];
+ const auto firstBulkElement = bulkProblem().model().globalFvGeometry().element(firstPair.first);
+
+ auto bulkFvGeometry = localView(bulkProblem().model().globalFvGeometry());
+ bulkFvGeometry.bind(firstBulkElement);
+
+ auto bulkElemVolVars = localView(bulkProblem().model().curGlobalVolVars());
+ bulkElemVolVars.bind(firstBulkElement, bulkFvGeometry, bulkProblem().model().curSol());
+
+ auto bulkElemFluxVarsCache = localView(bulkProblem().model().globalFluxVarsCache());
+ bulkElemFluxVarsCache.bind(firstBulkElement, bulkFvGeometry, bulkElemVolVars);
+
+ BulkPrimaryVariables flux(0.0);
+ for (const auto& pair : couplingData.elementScvfList)
+ {
+ const auto bulkElement = bulkProblem().model().globalFvGeometry().element(pair.first);
+ for (auto scvfIdx : pair.second)
+ flux += bulkLocalResidual_.computeFlux(bulkElement,
+ bulkFvGeometry,
+ bulkElemVolVars,
+ bulkFvGeometry.scvf(scvfIdx),
+ bulkElemFluxVarsCache);
+ }
+
+ // The source has to be formulated per volume
+ flux /= scv.volume()*curElemVolVars[scv].extrusionFactor();
+ return flux;
+ }
+
+ //! Compute lowDim volume variables for given bulk element and scvf
+ LowDimVolumeVariables lowDimVolVars(const BulkElement& element,
+ const BulkFVElementGeometry& fvGeometry,
+ const BulkSubControlVolumeFace& scvf) const
+ {
+ const auto& scvfCouplingData = couplingMapper_.getBulkCouplingData(element).getScvfCouplingData(scvf);
+
+ assert(scvfCouplingData.first && "no coupled facet element found for given scvf!");
+ const auto lowDimElement = lowDimProblem().model().globalFvGeometry().element(scvfCouplingData.second);
+ auto lowDimFvGeometry = localView(lowDimProblem().model().globalFvGeometry());
+ lowDimFvGeometry.bindElement(lowDimElement);
+
+ const auto& lowDimCurSol = lowDimProblem().model().curSol();
+ LowDimVolumeVariables lowDimVolVars;
+ lowDimVolVars.update(lowDimProblem().model().elementSolution(lowDimElement, lowDimCurSol),
+ lowDimProblem(),
+ lowDimElement,
+ lowDimFvGeometry.scv(scvfCouplingData.second));
+ return lowDimVolVars;
+ }
+
+ const BulkProblem& bulkProblem() const
+ { return *bulkProblemPtr_; }
+
+ const LowDimProblem& lowDimProblem() const
+ { return *lowDimProblemPtr_; }
+
+ const CCMpfaFacetCouplingMapper<TypeTag>& couplingMapper() const
+ { return couplingMapper_; }
+
+private:
+
+ //! evaluates if there is an interior boundary between two bulk elements
+ bool isInteriorBoundary_(BulkIndexType eIdxInside, BulkIndexType eIdxOutside) const
+ {
+ for (const auto& lowDimElement : elements(lowDimProblem().gridView()))
+ {
+ // The indices of bulk elements in which the actual low dim element is a facet
+ const auto& bulkElementIndices = GridCreator::getCoupledBulkElementIndices(lowDimElement);
+
+ // if inside & outside element are contained in them, this is a "coupling" intersection
+ if (BulkMpfaHelper::contains(bulkElementIndices, eIdxInside) && BulkMpfaHelper::contains(bulkElementIndices, eIdxOutside))
+ return true;
+ }
+
+ return false;
+ }
+
+ BulkProblem& bulkProblem()
+ { return *bulkProblemPtr_; }
+
+ LowDimProblem& lowDimProblem()
+ { return *lowDimProblemPtr_; }
+
+ //! Returns the implementation of the problem (i.e. static polymorphism)
+ Implementation &asImp_()
+ { return *static_cast<Implementation *>(this); }
+
+ //! \copydoc asImp_()
+ const Implementation &asImp_() const
+ { return *static_cast<const Implementation *>(this); }
+
+ BulkProblem *bulkProblemPtr_;
+ LowDimProblem *lowDimProblemPtr_;
+
+ mutable BulkLocalResidual bulkLocalResidual_;
+ mutable LowDimLocalResidual lowDimLocalResidual_;
+
+ CCMpfaFacetCouplingMapper<TypeTag> couplingMapper_;
+};
+
+} // end namespace
+
+#endif // DUMUX_FACETCOUPLINGMANAGER_HH
diff --git a/dumux/mixeddimension/facet/mpfa/couplingmapper.hh b/dumux/mixeddimension/facet/mpfa/couplingmapper.hh
new file mode 100644
index 0000000000..83168983a7
--- /dev/null
+++ b/dumux/mixeddimension/facet/mpfa/couplingmapper.hh
@@ -0,0 +1,257 @@
+/*****************************************************************************
+ * 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 Facet
+ * \brief Creates mappings of the entities from the bulk and the facet subdomain
+ */
+
+#ifndef DUMUX_MIXEDDIMENSION_CCMPFA_FACET_COUPLINGMAPPER_HH
+#define DUMUX_MIXEDDIMENSION_CCMPFA_FACET_COUPLINGMAPPER_HH
+
+#include <dune/common/version.hh>
+#include <dune/grid/common/mcmgmapper.hh>
+#include <dune/grid/common/gridfactory.hh>
+
+#include <dumux/common/math.hh>
+#include <dumux/common/propertysystem.hh>
+#include <dumux/common/parameters.hh>
+#include <dumux/common/valgrind.hh>
+
+namespace Dumux
+{
+
+namespace Properties
+{
+NEW_PROP_TAG(Scalar);
+NEW_PROP_TAG(BulkProblemTypeTag);
+NEW_PROP_TAG(LowDimProblemTypeTag);
+}
+
+/*!
+ * \ingroup FacetCoupling
+ * \brief Creates mappings of the entities from the bulk and the facet subdomain
+ */
+template<class TypeTag>
+class CCMpfaFacetCouplingMapper
+{
+ using GridCreator = typename GET_PROP_TYPE(TypeTag, GridCreator);
+
+ using BulkProblemTypeTag = typename GET_PROP_TYPE(TypeTag, BulkProblemTypeTag);
+ using LowDimProblemTypeTag = typename GET_PROP_TYPE(TypeTag, LowDimProblemTypeTag);
+
+ using BulkProblem = typename GET_PROP_TYPE(BulkProblemTypeTag, Problem);
+ using LowDimProblem = typename GET_PROP_TYPE(LowDimProblemTypeTag, Problem);
+
+ using BulkGridView = typename GET_PROP_TYPE(BulkProblemTypeTag, GridView);
+ using LowDimGridView = typename GET_PROP_TYPE(LowDimProblemTypeTag, GridView);
+
+ using BulkIndexType = typename BulkGridView::IndexSet::IndexType;
+ using LowDimIndexType = typename LowDimGridView::IndexSet::IndexType;
+
+ using BulkElement = typename BulkGridView::template Codim<0>::Entity;
+ using LowDimElement = typename LowDimGridView::template Codim<0>::Entity;
+
+ using BulkSubControlVolumeFace = typename GET_PROP_TYPE(BulkProblemTypeTag, SubControlVolumeFace);
+ using BulkMpfaHelper = typename GET_PROP_TYPE(BulkProblemTypeTag, MpfaHelper);
+
+ static constexpr int lowDimDim = LowDimGridView::dimension;
+
+public:
+
+ struct BulkDataInLowDim
+ {
+ // indicates whether or not this element is coupled
+ bool isCoupled;
+
+ // The list of bulk element indices connected to a low dim element
+ std::vector<BulkIndexType> couplingStencil;
+
+ // We store for each directly connected element a list of
+ // global scvf indices over which the two domains are coupled
+ std::vector<std::pair<BulkIndexType, std::vector<BulkIndexType>>> elementScvfList;
+
+ // The constructor
+ BulkDataInLowDim() : isCoupled(false) {}
+ };
+
+ struct LowDimDataInBulk
+ {
+ // indicates whether or not this element is coupled
+ bool isCoupled;
+
+ // The list of low dim elements with influence on this bulk element
+ std::vector<LowDimIndexType> couplingStencil;
+
+ // to each scvf of the bulk element we store
+ // the low dim element index that it is connected to
+ std::vector<std::pair<BulkIndexType, LowDimIndexType>> scvfToLowDimData;
+
+ // The constructor
+ LowDimDataInBulk() : isCoupled(false) {}
+
+ // For a given scvf of the element this method returns a pair of a bool and an index. The boolean
+ // indicates if this scvf does couple to a low dim element, the index is the low dim element index
+ std::pair<bool, LowDimIndexType> getScvfCouplingData(const BulkSubControlVolumeFace& scvf) const
+ { return getScvfCouplingData(scvf.index()); }
+
+ std::pair<bool, LowDimIndexType> getScvfCouplingData(BulkIndexType scvfIdx) const
+ {
+ for (const auto& pair : scvfToLowDimData)
+ if (pair.first == scvfIdx)
+ return std::make_pair(true, pair.second);
+ return std::make_pair(false, 0);
+ }
+ };
+
+ //! The constructor
+ CCMpfaFacetCouplingMapper(BulkProblem &bulkProblem, LowDimProblem &lowDimProblem)
+ : isInitialized_(false),
+ bulkProblemPtr_(&bulkProblem),
+ lowDimProblemPtr_(&lowDimProblem)
+ {}
+
+ /*!
+ * \brief initializes the maps
+ *
+ * \param bulkProblem The bulk sub problem
+ * \param lowDimProblem The lower-dimensional sub problem
+ */
+ void init()
+ {
+ const auto& bulkGridView = bulkProblem_().gridView();
+ const auto& lowDimGridView = lowDimProblem_().gridView();
+
+ lowDimCouplingData_.resize(lowDimGridView.size(0));
+ bulkCouplingData_.resize(bulkGridView.size(0));
+
+ for (const auto& lowDimElement : elements(lowDimGridView))
+ {
+ const auto lowDimElementGlobalIdx = lowDimProblem_().elementMapper().index(lowDimElement);
+
+ // The indices of bulk elements in which the actual low dim element is a facet
+ const auto& bulkElementIndices = GridCreator::getCoupledBulkElementIndices(lowDimElement);
+ if (bulkElementIndices.size() > 1)
+ {
+ lowDimCouplingData_[lowDimElementGlobalIdx].isCoupled = true;
+
+ for (auto bulkElementIdx : bulkElementIndices)
+ {
+ bulkCouplingData_[bulkElementIdx].isCoupled = true;
+
+ const auto bulkElement = bulkProblem_().model().globalFvGeometry().element(bulkElementIdx);
+ auto bulkFvGeometry = localView(bulkProblem_().model().globalFvGeometry());
+ bulkFvGeometry.bindElement(bulkElement);
+
+ // find the scvfs in the bulk element that "touch" the facet element
+ std::vector<BulkIndexType> scvfIndices;
+ for (const auto& scvf : scvfs(bulkFvGeometry))
+ {
+ if (scvf.boundary())
+ continue;
+
+ const auto anyOtherIdx = bulkElementIndices[0] == bulkElementIdx ?
+ bulkElementIndices[1] :
+ bulkElementIndices[0];
+
+ // check if any of the other bulk indices is in the outside indices of the scvf
+ if (BulkMpfaHelper::contains(scvf.outsideScvIndices(), anyOtherIdx))
+ {
+ scvfIndices.push_back(scvf.index());
+
+ // add data to the bulk data map
+ bulkCouplingData_[bulkElementIdx].scvfToLowDimData.emplace_back(std::make_pair(scvf.index(), lowDimElementGlobalIdx));
+
+ // also, insert scvf stencil to the coupling stencil of the low dim element
+ const auto scvfStencil = [&] ()
+ {
+ if (bulkProblem_().model().globalFvGeometry().isInBoundaryInteractionVolume(scvf))
+ return bulkProblem_().model().globalFvGeometry().boundaryInteractionVolumeSeed(scvf).globalScvIndices();
+ else
+ return bulkProblem_().model().globalFvGeometry().interactionVolumeSeed(scvf).globalScvIndices();
+ } ();
+
+ auto& lowDimCouplingStencil = lowDimCouplingData_[lowDimElementGlobalIdx].couplingStencil;
+ lowDimCouplingStencil.insert(lowDimCouplingStencil.end(), scvfStencil.begin(), scvfStencil.end());
+
+ // also, all the bulk elements in the scvf stencil will be coupled to the actual low dim element
+ for (auto bulkIdx : scvfStencil)
+ {
+ bulkCouplingData_[bulkIdx].isCoupled = true;
+ bulkCouplingData_[bulkIdx].couplingStencil.push_back(lowDimElementGlobalIdx);
+ }
+ }
+ }
+
+ // insert data into the low dim map
+ lowDimCouplingData_[lowDimElementGlobalIdx].elementScvfList.emplace_back(std::make_pair(bulkElementIdx, std::move(scvfIndices)));
+ }
+ }
+ else if (bulkElementIndices.size() == 1)
+ DUNE_THROW(Dune::NotImplemented, "Coupled facet elements on the bulk boundary");
+
+ // make the coupling stencil of the low dim element unique
+ auto& stencil = lowDimCouplingData_[lowDimElementGlobalIdx].couplingStencil;
+ std::sort(stencil.begin(), stencil.end());
+ stencil.erase(std::unique(stencil.begin(), stencil.end()), stencil.end());
+ }
+
+ // make the coupling stencils in the bulk map unique
+ for (auto& data : bulkCouplingData_)
+ {
+ auto& stencil = data.couplingStencil;
+ std::sort(stencil.begin(), stencil.end());
+ stencil.erase(std::unique(stencil.begin(), stencil.end()), stencil.end());
+ }
+
+ // state the initialization status
+ isInitialized_ = true;
+ }
+
+ const LowDimDataInBulk& getBulkCouplingData(const BulkElement& bulkElement) const
+ { return getBulkCouplingData(bulkProblem_().elementMapper().index(bulkElement)); }
+
+ const LowDimDataInBulk& getBulkCouplingData(BulkIndexType bulkElementIdx) const
+ { return bulkCouplingData_[bulkElementIdx]; }
+
+ const BulkDataInLowDim& getLowDimCouplingData(const LowDimElement& lowDimElement) const
+ { return getLowDimCouplingData(lowDimProblem_().elementMapper().index(lowDimElement)); }
+
+ const BulkDataInLowDim& getLowDimCouplingData(LowDimIndexType lowDimElementIdx) const
+ { return lowDimCouplingData_[lowDimElementIdx]; }
+
+ bool isInitialized() const
+ { return isInitialized_; }
+
+private:
+ const BulkProblem& bulkProblem_() const
+ { return *bulkProblemPtr_; }
+
+ const LowDimProblem& lowDimProblem_() const
+ { return *lowDimProblemPtr_; }
+
+ bool isInitialized_;
+ const BulkProblem* bulkProblemPtr_;
+ const LowDimProblem* lowDimProblemPtr_;
+
+ std::vector<LowDimDataInBulk> bulkCouplingData_;
+ std::vector<BulkDataInLowDim> lowDimCouplingData_;
+};
+
+} // end namespace
+
+#endif // DUMUX_FACETCOUPLINGMAPPER_HH
diff --git a/dumux/mixeddimension/facet/mpfa/darcyslaw.hh b/dumux/mixeddimension/facet/mpfa/darcyslaw.hh
new file mode 100644
index 0000000000..d1357c57de
--- /dev/null
+++ b/dumux/mixeddimension/facet/mpfa/darcyslaw.hh
@@ -0,0 +1,289 @@
+// -*- 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
+ * \brief This file contains the data which is required to calculate volume and mass fluxes of
+ * fluid phases over a face of a finite volume by means of the Darcy approximation in the
+ * presence of lower dimensional (coupled) elements living on the this domain's element facets.
+ */
+#ifndef DUMUX_DISCRETIZATION_CC_MPFA_FACET_DARCYS_LAW_HH
+#define DUMUX_DISCRETIZATION_CC_MPFA_FACET_DARCYS_LAW_HH
+
+#include <dumux/discretization/cellcentered/mpfa/darcyslaw.hh>
+#include <dumux/discretization/cellcentered/mpfa/facetypes.hh>
+
+namespace Dumux
+{
+
+/*!
+ * \ingroup DarcysLaw
+ * \brief Specialization of Darcy's Law for the CCMpfa method with lower dimensional
+ * elements living on the bulk elements' facets.
+ */
+template <class TypeTag>
+class CCMpfaFacetCouplingDarcysLaw : public DarcysLawImplementation<TypeTag, DiscretizationMethods::CCMpfa>
+{
+ using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+ using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
+ using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+ using Element = typename GridView::template Codim<0>::Entity;
+ using MpfaHelper = typename GET_PROP_TYPE(TypeTag, MpfaHelper);
+ using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
+ using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
+ using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
+ using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
+ using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);
+ using ElementFluxVariablesCache = typename GET_PROP_TYPE(TypeTag, ElementFluxVariablesCache);
+ using FluxVariablesCache = typename GET_PROP_TYPE(TypeTag, FluxVariablesCache);
+ using ElementSolutionVector = typename GET_PROP_TYPE(TypeTag, ElementSolutionVector);
+
+ // Always use the dynamic type for vectors (compatibility with the boundary)
+ using BoundaryInteractionVolume = typename GET_PROP_TYPE(TypeTag, BoundaryInteractionVolume);
+ using CoefficientVector = typename BoundaryInteractionVolume::Vector;
+
+ static constexpr int numPhases = GET_PROP_VALUE(TypeTag, NumPhases);
+ static constexpr bool facetCoupling = GET_PROP_VALUE(TypeTag, MpfaFacetCoupling);
+ static constexpr bool useTpfaBoundary = GET_PROP_VALUE(TypeTag, UseTpfaBoundary);
+ static constexpr bool enableInteriorBoundaries = GET_PROP_VALUE(TypeTag, EnableInteriorBoundaries);
+
+ //! The cache used in conjunction with the mpfa Darcy's Law
+ class MpfaFacetCouplingDarcysLawCache
+ {
+ // We always use the dynamic types here to be compatible on the boundary
+ using Stencil = typename BoundaryInteractionVolume::GlobalIndexSet;
+ using PositionVector = typename BoundaryInteractionVolume::PositionVector;
+
+ public:
+ //! update cached objects
+ template<class InteractionVolume>
+ void updateAdvection(const InteractionVolume& iv, const SubControlVolumeFace &scvf)
+ {
+ const auto& localFaceData = iv.getLocalFaceData(scvf);
+ // update the quantities that are equal for all phases
+ advectionVolVarsStencil_ = iv.volVarsStencil();
+ advectionVolVarsPositions_ = iv.volVarsPositions();
+ advectionTij_ = iv.getTransmissibilities(localFaceData);
+
+ // we will need the neumann flux transformation on interior Neumann boundaries
+ advectionCij_ = iv.getNeumannFluxTransformationCoefficients(localFaceData);
+
+ // The neumann fluxes always have to be set per phase
+ for (unsigned int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
+ phaseNeumannFluxes_[phaseIdx] = iv.getNeumannFlux(localFaceData, phaseIdx);
+ }
+
+ //! Returns the volume variables indices necessary for flux computation
+ //! This includes all participating boundary volume variables. Since we
+ //! do not allow mixed BC for the mpfa this is the same for all phases.
+ const Stencil& advectionVolVarsStencil() const
+ { return advectionVolVarsStencil_; }
+
+ //! Returns the position on which the volume variables live. This is
+ //! necessary as we need to evaluate gravity also for the boundary volvars
+ const PositionVector& advectionVolVarsPositions() const
+ { return advectionVolVarsPositions_; }
+
+ //! Returns the transmissibilities associated with the volume variables
+ //! All phases flow through the same rock, thus, tij are equal for all phases
+ const CoefficientVector& advectionTij() const
+ { return advectionTij_; }
+
+ //! Returns the vector of coefficients with which the vector of neumann boundary conditions
+ //! has to be multiplied in order to transform them on the scvf this cache belongs to
+ const CoefficientVector& advectionCij() const
+ { return advectionCij_; }
+
+ //! If the useTpfaBoundary property is set to false, the boundary conditions
+ //! are put into the local systems leading to possible contributions on all faces
+ Scalar advectionNeumannFlux(unsigned int phaseIdx) const
+ { return phaseNeumannFluxes_[phaseIdx]; }
+
+ private:
+ // Quantities associated with advection
+ Stencil advectionVolVarsStencil_;
+ PositionVector advectionVolVarsPositions_;
+ CoefficientVector advectionTij_;
+ CoefficientVector advectionCij_;
+ std::array<Scalar, numPhases> phaseNeumannFluxes_;
+ };
+
+public:
+ // state the discretization method this implementation belongs to
+ static const DiscretizationMethods myDiscretizationMethod = DiscretizationMethods::CCMpfa;
+
+ // state the new type for the corresponding cache
+ using Cache = MpfaFacetCouplingDarcysLawCache;
+
+ static Scalar flux(const Problem& problem,
+ const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars,
+ const SubControlVolumeFace& scvf,
+ const unsigned int phaseIdx,
+ const ElementFluxVariablesCache& elemFluxVarsCache)
+ {
+ static const bool gravity = GET_PARAM_FROM_GROUP(TypeTag, bool, Problem, EnableGravity);
+
+ const auto& fluxVarsCache = elemFluxVarsCache[scvf];
+ const auto& volVarsStencil = fluxVarsCache.advectionVolVarsStencil();
+ const auto& volVarsPositions = fluxVarsCache.advectionVolVarsPositions();
+ const auto& tij = fluxVarsCache.advectionTij();
+
+ const bool isInteriorBoundary = enableInteriorBoundaries && fluxVarsCache.isInteriorBoundary();
+
+ // Calculate the interface density for gravity evaluation
+ const auto rho = interpolateDensity(fvGeometry, elemVolVars, scvf, fluxVarsCache, phaseIdx, isInteriorBoundary);
+
+ // calculate Tij*pj
+ Scalar flux(0.0);
+ unsigned int localIdx = 0;
+ for (const auto volVarIdx : volVarsStencil)
+ {
+ const auto& volVars = elemVolVars[volVarIdx];
+ Scalar h = volVars.pressure(phaseIdx);
+
+ // if gravity is enabled, add gravitational acceleration
+ if (gravity)
+ {
+ // gravitational acceleration in the center of the actual element
+ const auto x = volVarsPositions[localIdx];
+ const auto g = problem.gravityAtPos(x);
+
+ h -= rho*(g*x);
+ }
+
+ flux += tij[localIdx++]*h;
+ }
+
+ // if no interior boundaries are present, return the flux
+ if (!enableInteriorBoundaries)
+ return useTpfaBoundary ? flux : flux + fluxVarsCache.advectionNeumannFlux(phaseIdx);
+
+ // Handle interior boundaries
+ flux += computeInteriorBoundaryContribution(problem, fvGeometry, elemVolVars, fluxVarsCache, phaseIdx, rho);
+
+ // return overall resulting flux
+ return useTpfaBoundary ? flux : flux + fluxVarsCache.advectionNeumannFlux(phaseIdx);
+ }
+
+ static Scalar interpolateDensity(const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars,
+ const SubControlVolumeFace& scvf,
+ const FluxVariablesCache& fluxVarsCache,
+ const unsigned int phaseIdx,
+ const bool isInteriorBoundary)
+ {
+ static const bool gravity = GET_PARAM_FROM_GROUP(TypeTag, bool, Problem, EnableGravity);
+
+ if (!gravity)
+ return Scalar(0.0);
+ else
+ {
+ // Return facet density on interior boundaries
+ if (isInteriorBoundary)
+ return fluxVarsCache.interiorBoundaryDataSelf().facetVolVars(fvGeometry).density(phaseIdx);
+
+ // use arithmetic mean of the densities around the scvf
+ if (!scvf.boundary())
+ {
+ Scalar rho = elemVolVars[scvf.insideScvIdx()].density(phaseIdx);
+ for (auto outsideIdx : scvf.outsideScvIndices())
+ rho += elemVolVars[outsideIdx].density(phaseIdx);
+ return rho/(scvf.outsideScvIndices().size()+1);
+ }
+ else
+ return elemVolVars[scvf.outsideScvIdx()].density(phaseIdx);
+ }
+ }
+
+ static Scalar computeInteriorBoundaryContribution(const Problem& problem,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars,
+ const FluxVariablesCache& fluxVarsCache,
+ unsigned int phaseIdx, Scalar rho)
+ {
+ static const bool gravity = GET_PARAM_FROM_GROUP(TypeTag, bool, Problem, EnableGravity);
+
+ // obtain the transmissibilites associated with all pressures
+ const auto& tij = fluxVarsCache.advectionTij();
+
+ // the interior dirichlet boundaries local indices start after
+ // the cell and the domain Dirichlet boundary pressures
+ const auto startIdx = fluxVarsCache.advectionVolVarsStencil().size();
+
+ // The vector of interior neumann fluxes
+ const auto& cij = fluxVarsCache.advectionCij();
+ CoefficientVector facetCouplingFluxes(cij.size(), 0.0);
+
+ // add interior Dirichlet boundary contributions
+ Scalar flux = 0.0;
+ for (auto&& data : fluxVarsCache.interiorBoundaryData())
+ {
+ // Add additional Dirichlet fluxes for interior Dirichlet faces
+ if (data.faceType() == MpfaFaceTypes::interiorDirichlet)
+ {
+ Scalar h = data.facetVolVars(fvGeometry).pressure(phaseIdx);
+
+ if (gravity)
+ {
+ const auto x = fvGeometry.scvf(data.scvfIndex()).ipGlobal();
+ const auto g = problem.gravityAtPos(x);
+
+ h -= rho*(g*x);
+ }
+
+ flux += tij[startIdx + data.localIndexInInteractionVolume()]*h;
+ }
+
+ // add neumann contributions
+ if (data.faceType() == MpfaFaceTypes::interiorNeumann)
+ {
+ // get the scvf corresponding to actual interior neumann face
+ const auto& curScvf = fvGeometry.scvf(data.scvfIndex());
+
+ // get the volvars of the actual interior neumann face
+ const auto facetVolVars = data.facetVolVars(fvGeometry, curScvf);
+
+ // calculate "leakage factor"
+ const auto n = curScvf.unitOuterNormal();
+ const auto v = [&] ()
+ {
+ auto res = n;
+ res *= -0.5*facetVolVars.extrusionFactor();
+ res += curScvf.ipGlobal();
+ res -= curScvf.facetCorner();
+ res /= res.two_norm2();
+ return res;
+ } ();
+
+ // add value to vector of interior neumann fluxes
+ facetCouplingFluxes[data.localIndexInInteractionVolume()] += facetVolVars.pressure(phaseIdx)*
+ curScvf.area()*
+ elemVolVars[curScvf.insideScvIdx()].extrusionFactor()*
+ MpfaHelper::nT_M_v(n, facetVolVars.permeability(), v);
+ }
+ }
+
+ return flux + cij*facetCouplingFluxes;
+ }
+};
+
+} // end namespace
+
+#endif
diff --git a/dumux/mixeddimension/facet/mpfa/fickslaw.hh b/dumux/mixeddimension/facet/mpfa/fickslaw.hh
new file mode 100644
index 0000000000..d97477569e
--- /dev/null
+++ b/dumux/mixeddimension/facet/mpfa/fickslaw.hh
@@ -0,0 +1,334 @@
+// -*- 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
+ * \brief This file contains the data which is required to calculate
+ * molar and mass fluxes of a component in a fluid phase over a face of a finite volume by means
+ * of Fick's Law for cell-centered MPFA models in the presence of lower dimensional (coupled) elements
+ * living on the this domain's element facets.
+ */
+#ifndef DUMUX_DISCRETIZATION_CC_MPFA_FACET_FICKS_LAW_HH
+#define DUMUX_DISCRETIZATION_CC_MPFA_FACET_FICKS_LAW_HH
+
+#include <dumux/discretization/cellcentered/mpfa/fickslaw.hh>
+#include <dumux/discretization/cellcentered/mpfa/facetypes.hh>
+
+namespace Dumux
+{
+/*!
+ * \ingroup CCMpfaFicksLaw
+ * \brief Specialization of Fick's Law for the CCMpfa method with lower dimensional
+ * elements living on the bulk elements' facets.
+ */
+template <class TypeTag>
+class CCMpfaFacetCouplingFicksLaw : public FicksLawImplementation<TypeTag, DiscretizationMethods::CCMpfa>
+{
+ using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+ using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
+ using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+ using MpfaHelper = typename GET_PROP_TYPE(TypeTag, MpfaHelper);
+ using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
+ using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
+ using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
+ using EffDiffModel = typename GET_PROP_TYPE(TypeTag, EffectiveDiffusivityModel);
+ using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
+ using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);
+ using ElementFluxVariablesCache = typename GET_PROP_TYPE(TypeTag, ElementFluxVariablesCache);
+ using FluxVariablesCache = typename GET_PROP_TYPE(TypeTag, FluxVariablesCache);
+
+ // Always use the dynamic type for vectors (compatibility with the boundary)
+ using BoundaryInteractionVolume = typename GET_PROP_TYPE(TypeTag, BoundaryInteractionVolume);
+ using CoefficientVector = typename BoundaryInteractionVolume::Vector;
+
+ using Element = typename GridView::template Codim<0>::Entity;
+ using IndexType = typename GridView::IndexSet::IndexType;
+
+ static constexpr int numPhases = GET_PROP_VALUE(TypeTag, NumPhases);
+ static constexpr bool useTpfaBoundary = GET_PROP_VALUE(TypeTag, UseTpfaBoundary);
+ static constexpr bool facetCoupling = GET_PROP_VALUE(TypeTag, MpfaFacetCoupling);
+ static constexpr bool enableInteriorBoundaries = GET_PROP_VALUE(TypeTag, EnableInteriorBoundaries);
+
+ //! The cache used in conjunction with the mpfa Fick's Law
+ class MpfaFicksLawCache
+ {
+ static const int numComponents = GET_PROP_VALUE(TypeTag, NumComponents);
+
+ // We always use the dynamic types here to be compatible on the boundary
+ using Stencil = typename BoundaryInteractionVolume::GlobalIndexSet;
+ using PositionVector = typename BoundaryInteractionVolume::PositionVector;
+
+ public:
+ //! The constructor. Initializes the Neumann flux to zero
+ MpfaFicksLawCache() { componentNeumannFluxes_.fill(0.0); }
+
+ // update cached objects for the diffusive fluxes
+ template<typename InteractionVolume>
+ void updateDiffusion(const InteractionVolume& iv, const SubControlVolumeFace &scvf,
+ unsigned int phaseIdx, unsigned int compIdx)
+ {
+ const auto& localFaceData = iv.getLocalFaceData(scvf);
+ diffusionTij_[phaseIdx][compIdx] = iv.getTransmissibilities(localFaceData);
+ // get the stencil only for the first call
+ if (phaseIdx == 0 && compIdx == 0)
+ diffusionVolVarsStencil_ = iv.volVarsStencil();
+
+ // we will need the neumann flux transformation on interior Neumann boundaries
+ diffusionCij_[phaseIdx][compIdx] = iv.getNeumannFluxTransformationCoefficients(localFaceData);
+
+ //! For compositional models, we associate neumann fluxes with the phases (main components)
+ //! This is done in the AdvectionCache. However, in single-phasic models we solve the phase AND
+ //! the component mass balance equations. Thus, in this case we have diffusive neumann contributions.
+ //! we assume compIdx = eqIdx
+ if (numPhases == 1 && phaseIdx != compIdx)
+ componentNeumannFluxes_[compIdx] = iv.getNeumannFlux(localFaceData, compIdx);
+
+ }
+
+ //! Returns the volume variables indices necessary for diffusive flux
+ //! computation. This includes all participating boundary volume variables
+ //! and it can be different for the phases & components.
+ const Stencil& diffusionVolVarsStencil(unsigned int phaseIdx, unsigned int compIdx) const
+ { return diffusionVolVarsStencil_; }
+
+ //! Returns the transmissibilities associated with the volume variables
+ //! This can be different for the phases & components.
+ const CoefficientVector& diffusionTij(unsigned int phaseIdx, unsigned int compIdx) const
+ { return diffusionTij_[phaseIdx][compIdx]; }
+
+ //! Returns the vector of coefficients with which the vector of neumann boundary conditions
+ //! has to be multiplied in order to transform them on the scvf this cache belongs to
+ const CoefficientVector& diffusionCij(unsigned int phaseIdx, unsigned int compIdx) const
+ { return diffusionCij_[phaseIdx][compIdx]; }
+
+ //! If the useTpfaBoundary property is set to false, the boundary conditions
+ //! are put into the local systems leading to possible contributions on all faces
+ Scalar componentNeumannFlux(unsigned int compIdx) const
+ {
+ assert(numPhases == 1);
+ return componentNeumannFluxes_[compIdx];
+ }
+
+ private:
+ // Quantities associated with molecular diffusion
+ Stencil diffusionVolVarsStencil_;
+ std::array< std::array<CoefficientVector, numComponents>, numPhases> diffusionTij_;
+ std::array< std::array<CoefficientVector, numComponents>, numPhases> diffusionCij_;
+
+ // diffusive neumann flux for single-phasic models
+ std::array<Scalar, numComponents> componentNeumannFluxes_;
+ };
+
+public:
+ // state the discretization method this implementation belongs to
+ static const DiscretizationMethods myDiscretizationMethod = DiscretizationMethods::CCMpfa;
+
+ // state the new type for the corresponding cache
+ using Cache = MpfaFicksLawCache;
+
+ static Scalar flux(const Problem& problem,
+ const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars,
+ const SubControlVolumeFace& scvf,
+ int phaseIdx, int compIdx,
+ const ElementFluxVariablesCache& elemFluxVarsCache,
+ bool useMoles = true)
+ {
+ const auto& fluxVarsCache = elemFluxVarsCache[scvf];
+ const auto& volVarsStencil = fluxVarsCache.diffusionVolVarsStencil(phaseIdx, compIdx);
+ const auto& tij = fluxVarsCache.diffusionTij(phaseIdx, compIdx);
+
+ const bool isInteriorBoundary = enableInteriorBoundaries && fluxVarsCache.isInteriorBoundary();
+
+ // get the scaling factor for the effective diffusive fluxes
+ const auto effFactor = computeEffectivityFactor(fvGeometry, elemVolVars, scvf, fluxVarsCache, phaseIdx, isInteriorBoundary);
+
+ // if factor is zero, the flux will end up zero anyway
+ if (effFactor == 0.0)
+ return 0.0;
+
+ // lambda functions depending on if we use mole or mass fractions
+ auto getX = [useMoles, phaseIdx, compIdx] (const auto& volVars)
+ { return useMoles ? volVars.moleFraction(phaseIdx, compIdx) : volVars.massFraction(phaseIdx, compIdx); };
+
+ auto getRho = [useMoles, phaseIdx] (const auto& volVars)
+ { return useMoles ? volVars.molarDensity(phaseIdx) : volVars.density(phaseIdx); };
+
+ // calculate the density at the interface
+ const auto rho = interpolateDensity(fvGeometry, elemVolVars, scvf, fluxVarsCache, getRho, isInteriorBoundary);
+
+ // calculate Tij*xj
+ Scalar flux(0.0);
+ unsigned int localIdx = 0;
+ for (const auto volVarIdx : volVarsStencil)
+ flux += tij[localIdx++]*getX(elemVolVars[volVarIdx]);
+
+ // if no interior boundaries are present, return effective mass flux
+ if (!enableInteriorBoundaries)
+ return useTpfaBoundary ? flux*rho*effFactor : flux*rho*effFactor + fluxVarsCache.componentNeumannFlux(compIdx);
+
+ // Handle interior boundaries
+ flux += computeInteriorBoundaryContribution(fvGeometry, elemVolVars, fluxVarsCache, getX, phaseIdx, compIdx);
+
+ // return overall resulting flux
+ return useTpfaBoundary ? flux*rho*effFactor : flux*rho*effFactor + fluxVarsCache.componentNeumannFlux(compIdx);
+ }
+
+ template<typename GetRhoFunction>
+ static Scalar interpolateDensity(const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars,
+ const SubControlVolumeFace& scvf,
+ const FluxVariablesCache& fluxVarsCache,
+ const GetRhoFunction& getRho,
+ const bool isInteriorBoundary)
+ {
+
+ // maybe use the density of the interior BC on the facet
+ if (isInteriorBoundary)
+ return getRho(fluxVarsCache.interiorBoundaryDataSelf().facetVolVars(fvGeometry));
+
+ // use arithmetic mean of the densities around the scvf
+ if (!scvf.boundary())
+ {
+ Scalar rho = getRho(elemVolVars[scvf.insideScvIdx()]);
+ for (auto outsideIdx : scvf.outsideScvIndices())
+ rho += getRho(elemVolVars[outsideIdx]);
+ return rho/(scvf.outsideScvIndices().size()+1);
+ }
+ else
+ return getRho(elemVolVars[scvf.outsideScvIdx()]);
+ }
+
+ //! Here we want to calculate the factors with which the diffusion coefficient has to be
+ //! scaled to get the effective diffusivity. For this we use the effective diffusivity with
+ //! a diffusion coefficient of 1.0 as input. Then we scale the transmissibilites during flux
+ //! calculation (above) with the harmonic average of the two factors
+ static Scalar computeEffectivityFactor(const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars,
+ const SubControlVolumeFace& scvf,
+ const FluxVariablesCache& fluxVarsCache,
+ const unsigned int phaseIdx,
+ const bool isInteriorBoundary)
+ {
+ if (isInteriorBoundary)
+ {
+ // use harmonic mean between the interior and the facet volvars
+ const auto& data = fluxVarsCache.interiorBoundaryDataSelf();
+
+ const auto& insideVolVars = elemVolVars[scvf.insideScvIdx()];
+ const auto factor = EffDiffModel::effectiveDiffusivity(insideVolVars.porosity(),
+ insideVolVars.saturation(phaseIdx),
+ /*Diffusion coefficient*/ 1.0);
+
+ const auto facetVolVars = data.facetVolVars(fvGeometry);
+ const auto outsideFactor = EffDiffModel::effectiveDiffusivity(facetVolVars.porosity(),
+ facetVolVars.saturation(phaseIdx),
+ /*Diffusion coefficient*/ 1.0);
+
+ return harmonicMean(factor, outsideFactor);
+ }
+
+ // use the harmonic mean between inside and outside
+ const auto& insideVolVars = elemVolVars[scvf.insideScvIdx()];
+ const auto factor = EffDiffModel::effectiveDiffusivity(insideVolVars.porosity(),
+ insideVolVars.saturation(phaseIdx),
+ /*Diffusion coefficient*/ 1.0);
+
+ if (!scvf.boundary())
+ {
+ // interpret outside factor as arithmetic mean
+ Scalar outsideFactor = 0.0;
+ for (auto outsideIdx : scvf.outsideScvIndices())
+ {
+ const auto& outsideVolVars = elemVolVars[outsideIdx];
+ outsideFactor += EffDiffModel::effectiveDiffusivity(outsideVolVars.porosity(),
+ outsideVolVars.saturation(phaseIdx),
+ /*Diffusion coefficient*/ 1.0);
+ }
+ outsideFactor /= scvf.outsideScvIndices().size();
+
+ // use the harmonic mean of the two
+ return harmonicMean(factor, outsideFactor);
+ }
+
+ return factor;
+ }
+
+ template<typename GetXFunction>
+ static Scalar computeInteriorBoundaryContribution(const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars,
+ const FluxVariablesCache& fluxVarsCache,
+ const GetXFunction& getX,
+ unsigned int phaseIdx, unsigned int compIdx)
+ {
+ // obtain the transmissibilites associated with all pressures
+ const auto& tij = fluxVarsCache.diffusionTij(phaseIdx, compIdx);
+
+ // the interior dirichlet boundaries local indices start after
+ // the cell and the domain Dirichlet boundary pressures
+ const auto startIdx = fluxVarsCache.diffusionVolVarsStencil(phaseIdx, compIdx).size();
+
+ // The vector of interior neumann fluxes
+ const auto& cij = fluxVarsCache.advectionCij();
+ CoefficientVector facetCouplingFluxes(cij.size(), 0.0);
+
+ // add interior Dirichlet boundary contributions
+ Scalar flux = 0.0;
+ for (auto&& data : fluxVarsCache.interiorBoundaryData())
+ {
+ // Add additional Dirichlet fluxes for interior Dirichlet faces
+ if (data.faceType() == MpfaFaceTypes::interiorDirichlet)
+ flux += tij[startIdx + data.localIndexInInteractionVolume()]*getX(data.facetVolVars(fvGeometry));
+
+ // add neumann contributions
+ if (data.faceType() == MpfaFaceTypes::interiorNeumann)
+ {
+ // get the scvf corresponding to actual interior neumann face
+ const auto& curScvf = fvGeometry.scvf(data.scvfIndex());
+
+ // get the volvars of the actual interior neumann face
+ const auto facetVolVars = data.facetVolVars(fvGeometry, curScvf);
+
+ // calculate "leakage factor"
+ const auto n = curScvf.unitOuterNormal();
+ const auto v = [&] ()
+ {
+ auto res = n;
+ res *= -0.5*facetVolVars.extrusionFactor();
+ res += curScvf.ipGlobal();
+ res -= curScvf.facetCorner();
+ res /= res.two_norm2();
+ return res;
+ } ();
+
+ // add value to vector of interior neumann fluxes
+ facetCouplingFluxes[data.localIndexInInteractionVolume()] += getX(facetVolVars)*
+ curScvf.area()*
+ elemVolVars[curScvf.insideScvIdx()].extrusionFactor()*
+ MpfaHelper::nT_M_v(n, facetVolVars.diffusionCoefficient(phaseIdx, compIdx), v);
+ }
+ }
+
+ return flux + cij*facetCouplingFluxes;
+ }
+};
+
+} // end namespace
+
+#endif
diff --git a/dumux/mixeddimension/facet/mpfa/interactionvolume.hh b/dumux/mixeddimension/facet/mpfa/interactionvolume.hh
new file mode 100644
index 0000000000..92410e3bf1
--- /dev/null
+++ b/dumux/mixeddimension/facet/mpfa/interactionvolume.hh
@@ -0,0 +1,149 @@
+// -*- 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
+ * \brief Base classes for interaction volumes of mpfa models with active coupling over the element facets.
+ */
+#ifndef DUMUX_MIXEDDIMENSION_FACET_INTERACTIONVOLUME_HH
+#define DUMUX_MIXEDDIMENSION_FACET_INTERACTIONVOLUME_HH
+
+#include <dumux/discretization/cellcentered/mpfa/interactionvolume.hh>
+#include <dumux/discretization/cellcentered/mpfa/facetypes.hh>
+#include <dumux/discretization/cellcentered/mpfa/methods.hh>
+
+namespace Dumux
+{
+// forward declaration of the implementation
+template<class TypeTag, MpfaMethods Method>
+class CCMpfaFacetCouplingInteractionVolumeImplementation;
+
+/*!
+ * \ingroup MixedDimension
+ * \brief Base class for the interaction volumes of the mpfa method with active coupling over the element facets.
+ */
+template<class TypeTag>
+using CCMpfaFacetCouplingInteractionVolume = CCMpfaFacetCouplingInteractionVolumeImplementation<TypeTag, GET_PROP_VALUE(TypeTag, MpfaMethod)>;
+
+// Per default, we inherit from the standard interaction volumes
+template<class TypeTag, MpfaMethods Method>
+class CCMpfaFacetCouplingInteractionVolumeImplementation : public CCMpfaInteractionVolumeImplementation<TypeTag, Method> {};
+
+// the o-method interaction volume is substituted by the one including data on the facet element's
+// tensorial quantities into the local system to be solved. This has to be used as boundary interaction volume
+template<class TypeTag>
+class CCMpfaFacetCouplingInteractionVolumeImplementation<TypeTag, MpfaMethods::oMethod>
+ : public CCMpfaInteractionVolumeImplementation<TypeTag, MpfaMethods::oMethod>
+{
+ using ParentType = CCMpfaInteractionVolumeImplementation<TypeTag, MpfaMethods::oMethod>;
+
+ using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+ using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+ using Element = typename GridView::template Codim<0>::Entity;
+ using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
+ using MpfaHelper = typename GET_PROP_TYPE(TypeTag, MpfaHelper);
+ using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
+ using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
+ using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
+ using InteriorBoundaryData = typename GET_PROP_TYPE(TypeTag, InteriorBoundaryData);
+ using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);
+
+ using LocalScvfType = typename ParentType::Traits::LocalScvfType;
+
+ static constexpr bool useTpfaBoundary = GET_PROP_VALUE(TypeTag, UseTpfaBoundary);
+
+public:
+ using typename ParentType::LocalIndexType;
+ using typename ParentType::Seed;
+
+ CCMpfaFacetCouplingInteractionVolumeImplementation(const Seed& seed,
+ const Problem& problem,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars)
+ : ParentType(seed, problem, fvGeometry, elemVolVars)
+ {}
+
+public:
+ // We include data on the tensorial quantities of the facet elements here
+ template<typename GetTensorFunction>
+ Scalar interiorNeumannTerm(const GetTensorFunction& getTensor,
+ const Element& element,
+ const LocalScvfType& localScvf,
+ const InteriorBoundaryData& data) const
+ {
+ // obtain the complete data on the facet element
+ const auto completeFacetData = data.completeCoupledFacetData(this->fvGeometry_());
+
+ // calculate "leakage factor"
+ const auto n = localScvf.unitOuterNormal();
+ const auto v = [&] ()
+ {
+ auto res = n;
+ res *= -0.5*completeFacetData.volVars().extrusionFactor();
+ res += localScvf.ip();
+ res -= localScvf.globalScvf().facetCorner();
+ res /= res.two_norm2();
+ return res;
+ } ();
+
+ // substract (n*T*v)*Area from diagonal matrix entry
+ const auto facetTensor = getTensor(completeFacetData.problem(),
+ completeFacetData.element(),
+ completeFacetData.volVars(),
+ completeFacetData.fvGeometry(),
+ completeFacetData.scv());
+
+ return localScvf.area()*
+ this->elemVolVars_()[localScvf.insideGlobalScvIndex()].extrusionFactor()*
+ MpfaHelper::nT_M_v(n, facetTensor, v);
+ }
+
+ void assembleNeumannFluxVector_()
+ {
+ // initialize the neumann fluxes vector to zero
+ this->neumannFluxes_.resize(this->fluxFaceIndexSet_.size(), PrimaryVariables(0.0));
+
+ if (!this->onDomainOrInteriorBoundary() || useTpfaBoundary)
+ return;
+
+ LocalIndexType fluxFaceIdx = 0;
+ for (auto localFluxFaceIdx : this->fluxFaceIndexSet_)
+ {
+ const auto& localScvf = this->localScvf_(localFluxFaceIdx);
+ const auto faceType = localScvf.faceType();
+
+ if (faceType == MpfaFaceTypes::neumann)
+ {
+ const auto& element = this->localElement_(localScvf.insideLocalScvIndex());
+ const auto& globalScvf = this->fvGeometry_().scvf(localScvf.insideGlobalScvfIndex());
+ auto neumannFlux = this->problem_().neumann(element, this->fvGeometry_(), this->elemVolVars_(), globalScvf);
+ neumannFlux *= globalScvf.area();
+ neumannFlux *= this->elemVolVars_()[globalScvf.insideScvIdx()].extrusionFactor();
+
+ // The flux is assumed to be prescribed in the form of -D*gradU
+ this->neumannFluxes_[fluxFaceIdx] = neumannFlux;
+ }
+
+ fluxFaceIdx++;
+ }
+ }
+};
+
+} // end namespace
+
+#endif
diff --git a/dumux/mixeddimension/facet/mpfa/interiorboundarydata.hh b/dumux/mixeddimension/facet/mpfa/interiorboundarydata.hh
new file mode 100644
index 0000000000..6a0fe64c23
--- /dev/null
+++ b/dumux/mixeddimension/facet/mpfa/interiorboundarydata.hh
@@ -0,0 +1,201 @@
+// -*- 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
+ * \brief A class to store info on interior boundaries
+ */
+#ifndef DUMUX_DISCRETIZATION_CC_MPFA_FACET_INTERIORBOUNDARYDATA_HH
+#define DUMUX_DISCRETIZATION_CC_MPFA_FACET_INTERIORBOUNDARYDATA_HH
+
+#include <dumux/discretization/cellcentered/mpfa/facetypes.hh>
+#include <dumux/mixeddimension/properties.hh>
+
+namespace Dumux
+{
+
+template<class TypeTag>
+class CCMpfaFacetCouplingInteriorBoundaryData
+{
+ // types associated with the low dim domain
+ using GlobalProblemTypeTag = typename GET_PROP_TYPE(TypeTag, GlobalProblemTypeTag);
+ using LowDimProblemTypeTag = typename GET_PROP_TYPE(GlobalProblemTypeTag, LowDimProblemTypeTag);
+
+ using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
+ using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+ using Element = typename GridView::template Codim<0>::Entity;
+ using SpatialParams = typename GET_PROP_TYPE(TypeTag, SpatialParams);
+ using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
+ using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
+ using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
+ using ElementSolutionVector = typename GET_PROP_TYPE(TypeTag, ElementSolutionVector);
+ using BoundaryInteractionVolume = typename GET_PROP_TYPE(TypeTag, BoundaryInteractionVolume);
+
+ using IndexType = typename GridView::IndexSet::IndexType;
+ using LocalIndexType = typename BoundaryInteractionVolume::LocalIndexType;
+
+ using LowDimProblem = typename GET_PROP_TYPE(LowDimProblemTypeTag, Problem);
+ using LowDimSpatialParams = typename GET_PROP_TYPE(LowDimProblemTypeTag, SpatialParams);
+ using LowDimVolumeVariables = typename GET_PROP_TYPE(LowDimProblemTypeTag, VolumeVariables);
+ using LowDimGridView = typename GET_PROP_TYPE(LowDimProblemTypeTag, GridView);
+ using LowDimElement = typename LowDimGridView::template Codim<0>::Entity;
+ using LowDimFVElementGeometry = typename GET_PROP_TYPE(LowDimProblemTypeTag, FVElementGeometry);
+ using LowDimSubControlVolume = typename GET_PROP_TYPE(LowDimProblemTypeTag, SubControlVolume);
+
+ //! Dummy type for the CompleteCoupledFacetData struct.
+ //! Implementations need to have at least the provided interfaces.
+ //! Note that the return types are also "wrong" (Here just to satisfy the compiler)
+ struct CompleteCoupledFacetData
+ {
+ const LowDimProblem& problem() const
+ { return *lowDimProblemPtr_; }
+
+ const LowDimSpatialParams& spatialParams() const
+ { return *lowDimSpatialParamsPtr_; }
+
+ const LowDimVolumeVariables& volVars() const
+ { return lowDimVolVars_; }
+
+ const LowDimElement& element() const
+ { return lowDimElement_;}
+
+ const LowDimFVElementGeometry& fvGeometry() const
+ { return lowDimFvGeometry_; }
+
+ const LowDimSubControlVolume& scv() const
+ { return lowDimScv_; }
+
+ CompleteCoupledFacetData(const LowDimProblem& problem,
+ const IndexType elementIndex,
+ LowDimElement&& element,
+ LowDimFVElementGeometry&& fvGeometry,
+ LowDimVolumeVariables&& volVars)
+ : lowDimProblemPtr_(&problem),
+ lowDimSpatialParamsPtr_(&problem.spatialParams()),
+ lowDimElement_(std::move(element)),
+ lowDimFvGeometry_(std::move(fvGeometry)),
+ lowDimVolVars_(std::move(volVars)),
+ lowDimScv_(lowDimFvGeometry_.scv(elementIndex))
+ {}
+
+ private:
+ const LowDimProblem* lowDimProblemPtr_;
+ const LowDimSpatialParams* lowDimSpatialParamsPtr_;
+
+ LowDimElement lowDimElement_;
+ LowDimFVElementGeometry lowDimFvGeometry_;
+ LowDimVolumeVariables lowDimVolVars_;
+ const LowDimSubControlVolume& lowDimScv_;
+ };
+
+public:
+ //! the constructor
+ CCMpfaFacetCouplingInteriorBoundaryData(const Problem& problem,
+ IndexType elementIndex,
+ IndexType scvfIndex,
+ LocalIndexType localIndex,
+ MpfaFaceTypes faceType)
+ : problemPtr_(&problem),
+ elementIndex_(elementIndex),
+ scvfIndex_(scvfIndex),
+ localIndex_(localIndex),
+ faceType_(faceType)
+ {}
+
+ //! returns the global index of the element/scv connected to the interior boundary
+ IndexType elementIndex() const
+ { return elementIndex_; }
+
+ //! returns the global index of the scvf connected to the interior boundary
+ IndexType scvfIndex() const
+ { return scvfIndex_; }
+
+ //! returns the local index i of the scvf within the interaction volume.
+ //! This is either:
+ //! - the i-th flux face index (interior neumann boundaries)
+ //! - the i-th interior dirichlet face (interior dirichlet boundaries)
+ LocalIndexType localIndexInInteractionVolume() const
+ { return localIndex_; }
+
+ //! returns the face type of this scvf
+ MpfaFaceTypes faceType() const
+ { return faceType_; }
+
+ //! returns the volume variables for interior dirichlet boundaries
+ LowDimVolumeVariables facetVolVars(const FVElementGeometry& fvGeometry) const
+ {
+ return problem_().couplingManager().lowDimVolVars(problem_().model().globalFvGeometry().element(elementIndex()),
+ fvGeometry,
+ fvGeometry.scvf(scvfIndex()));
+ }
+
+ //! returns the volume variables for interior dirichlet boundaries
+ LowDimVolumeVariables facetVolVars(const FVElementGeometry& fvGeometry, const SubControlVolumeFace& scvf) const
+ {
+ assert(scvf.index() == scvfIndex() && "calling facet volume variables for an scvf other than the bound one");
+ return problem_().couplingManager().lowDimVolVars(problem_().model().globalFvGeometry().element(elementIndex()),
+ fvGeometry,
+ scvf);
+ }
+
+ //! The following interface is here for compatibility reasonsto be overloaded for problems using facet coupling.
+ //! prepares all the necessary variables of the other domain.
+ //! Note that also an implementation of the CompleteFacetData structure has to be provided.
+ CompleteCoupledFacetData completeCoupledFacetData(const FVElementGeometry& fvGeometry) const
+ {
+ const auto& couplingMapper = problem_().couplingManager().couplingMapper();
+
+ // get coupling data for this scvf
+ const auto element = problem_().model().globalFvGeometry().element(elementIndex());
+ const auto& scvfCouplingData = couplingMapper.getBulkCouplingData(element).getScvfCouplingData(fvGeometry.scvf(scvfIndex()));
+
+ // obtain data necessary to fully instantiate the complete coupled facet data
+ assert(!scvfCouplingData.first && "no coupled facet element found for given scvf!");
+ const auto& lowDimProblem = problem_().couplingManager().lowDimProblem();
+
+ auto lowDimElement = lowDimProblem.model().globalFvGeometry().element(scvfCouplingData.second);
+ auto lowDimFvGeometry = localView(lowDimProblem.model().globalFvGeometry());
+ lowDimFvGeometry.bindElement(lowDimElement);
+
+ LowDimVolumeVariables lowDimVolVars;
+ lowDimVolVars.update(lowDimProblem.model().elementSolution(lowDimElement, lowDimProblem.model().curSol()),
+ lowDimProblem,
+ lowDimElement,
+ lowDimFvGeometry.scv(scvfCouplingData.second));
+
+ return CompleteCoupledFacetData(lowDimProblem,
+ scvfCouplingData.second,
+ std::move(lowDimElement),
+ std::move(lowDimFvGeometry),
+ std::move(lowDimVolVars));
+ }
+
+private:
+ const Problem& problem_() const
+ { return *problemPtr_; }
+
+ const Problem* problemPtr_;
+
+ IndexType elementIndex_;
+ IndexType scvfIndex_;
+ LocalIndexType localIndex_;
+ MpfaFaceTypes faceType_;
+};
+} // end namespace
+
+#endif
diff --git a/dumux/mixeddimension/facet/mpfa/properties.hh b/dumux/mixeddimension/facet/mpfa/properties.hh
new file mode 100644
index 0000000000..32e11462ac
--- /dev/null
+++ b/dumux/mixeddimension/facet/mpfa/properties.hh
@@ -0,0 +1,72 @@
+// -*- 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 MixedDimension
+ * \brief Base properties for the bulk problems in mixed dimensional models
+ * with a lower dimensional model living on the element facets.
+ */
+
+#ifndef DUMUX_FACET_MIXEDDIMENSION_PROPERTIES_HH
+#define DUMUX_FACET_MIXEDDIMENSION_PROPERTIES_HH
+
+#include <dumux/implicit/cellcentered/mpfa/properties.hh>
+
+#include <dumux/mixeddimension/subproblemproperties.hh>
+#include <dumux/mixeddimension/facet/mpfa/interactionvolume.hh>
+#include <dumux/mixeddimension/facet/mpfa/interiorboundarydata.hh>
+#include <dumux/mixeddimension/facet/mpfa/darcyslaw.hh>
+#include <dumux/mixeddimension/facet/mpfa/fickslaw.hh>
+
+namespace Dumux
+{
+
+namespace Properties
+{
+NEW_TYPE_TAG(FacetCouplingBulkMpfaModel, INHERITS_FROM(CCMpfaModel));
+
+//! The interaction volume class
+SET_TYPE_PROP(FacetCouplingBulkMpfaModel, InteractionVolume, CCMpfaFacetCouplingInteractionVolume<TypeTag>);
+
+//! The boundary interaction volume class (for methods other than the omethod)
+SET_TYPE_PROP(FacetCouplingBulkMpfaModel, BoundaryInteractionVolume, CCMpfaFacetCouplingInteractionVolumeImplementation<TypeTag, MpfaMethods::oMethod>);
+
+//! The interior boundary data class
+SET_TYPE_PROP(FacetCouplingBulkMpfaModel, InteriorBoundaryData, CCMpfaFacetCouplingInteriorBoundaryData<TypeTag>);
+
+//! Ẃe always enable interior boundaries
+SET_BOOL_PROP(FacetCouplingBulkMpfaModel, EnableInteriorBoundaries, true);
+
+//! Facet coupling is always true here
+SET_BOOL_PROP(FacetCouplingBulkMpfaModel, MpfaFacetCoupling, true);
+
+//! Darcy's Law
+SET_TYPE_PROP(FacetCouplingBulkMpfaModel, AdvectionType, CCMpfaFacetCouplingDarcysLaw<TypeTag>);
+
+//! Ficks's Law
+SET_TYPE_PROP(FacetCouplingBulkMpfaModel, MolecularDiffusionType, CCMpfaFacetCouplingFicksLaw<TypeTag>);
+
+//! Darcy's Law
+// SET_TYPE_PROP(FacetCouplingBulkMpfaModel, AdvectionType, CCMpfaFacetCouplingFouriersLaw<TypeTag>);
+
+}//end namespace Properties
+
+}//end namespace Dumux
+
+#endif
diff --git a/dumux/mixeddimension/facet/start.hh b/dumux/mixeddimension/facet/start.hh
new file mode 100644
index 0000000000..bd59fe1673
--- /dev/null
+++ b/dumux/mixeddimension/facet/start.hh
@@ -0,0 +1,203 @@
+// -*- 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
+ * \brief Provides a few default main functions for convenience.
+ */
+#ifndef DUMUX_MIXEDDIMENSION_FACET_START_HH
+#define DUMUX_MIXEDDIMENSION_FACET_START_HH
+
+#include <ctime>
+#include <iostream>
+
+#include <dune/common/parallel/mpihelper.hh>
+#include <dune/grid/io/file/dgfparser/dgfexception.hh>
+
+#include <dumux/common/propertysystem.hh>
+#include <dumux/common/parameters.hh>
+#include <dumux/common/valgrind.hh>
+#include <dumux/common/dumuxmessage.hh>
+#include <dumux/common/defaultusagemessage.hh>
+#include <dumux/common/parameterparser.hh>
+
+namespace Dumux
+{
+// forward declaration of property tags
+namespace Properties
+{
+NEW_PROP_TAG(Scalar);
+NEW_PROP_TAG(GridCreator);
+NEW_PROP_TAG(Problem);
+NEW_PROP_TAG(TimeManager);
+}
+
+/*!
+ * \ingroup Start
+ *
+ * \brief Provides a main function which reads in parameters from the
+ * command line and a parameter file.
+ *
+ * \tparam TypeTag The type tag of the problem which needs to be solved
+ *
+ * \param argc The 'argc' argument of the main function: count of arguments (1 if there are no arguments)
+ * \param argv The 'argv' argument of the main function: array of pointers to the argument strings
+ * \param usage Callback function for printing the usage message
+ */
+template <class TypeTag>
+int start_(int argc,
+ char **argv,
+ void (*usage)(const char *, const std::string &))
+{
+ // some aliases for better readability
+ using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+ using GridCreator = typename GET_PROP_TYPE(TypeTag, GridCreator);
+ using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
+ using TimeManager = typename GET_PROP_TYPE(TypeTag, TimeManager);
+ using ParameterTree = typename GET_PROP(TypeTag, ParameterTree);
+
+ // initialize MPI, finalize is done automatically on exit
+ const auto& mpiHelper = Dune::MPIHelper::instance(argc, argv);
+
+ // print dumux start message
+ if (mpiHelper.rank() == 0)
+ DumuxMessage::print(/*firstCall=*/true);
+
+ ////////////////////////////////////////////////////////////
+ // parse the command line arguments and input file
+ ////////////////////////////////////////////////////////////
+
+ // if the user just wanted to see the help / usage message show usage and stop program
+ if(!ParameterParser::parseCommandLineArguments(argc, argv, ParameterTree::tree(), usage))
+ {
+ usage(argv[0], defaultUsageMessage(argv[0]));
+ return 0;
+ }
+ // parse the input file into the parameter tree
+ // check first if the user provided an input file through the command line, if not use the default
+ const auto parameterFileName = ParameterTree::tree().hasKey("ParameterFile") ? GET_RUNTIME_PARAM(TypeTag, std::string, ParameterFile) : "";
+ ParameterParser::parseInputFile(argc, argv, ParameterTree::tree(), parameterFileName, usage);
+
+ ////////////////////////////////////////////////////////////
+ // check for some user debugging parameters
+ ////////////////////////////////////////////////////////////
+
+ bool printProps = false; // per default don't print all properties
+ if (ParameterTree::tree().hasKey("PrintProperties") || ParameterTree::tree().hasKey("TimeManager.PrintProperties"))
+ printProps = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, bool, TimeManager, PrintProperties);
+
+ if (printProps && mpiHelper.rank() == 0)
+ Properties::print<TypeTag>();
+
+ bool printParams = true; // per default print all properties
+ if (ParameterTree::tree().hasKey("PrintParameters") || ParameterTree::tree().hasKey("TimeManager.PrintParameters"))
+ printParams = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, bool, TimeManager, PrintParameters);
+
+ //////////////////////////////////////////////////////////////////////
+ // try to create a grid (from the given grid file or the input file)
+ /////////////////////////////////////////////////////////////////////
+
+ try { GridCreator::makeGrid(); }
+ catch (...) {
+ std::string usageMessage = "\n\t -> Creation of the grid failed! <- \n\n";
+ usageMessage += defaultUsageMessage(argv[0]);
+ usage(argv[0], usageMessage);
+ throw;
+ }
+ GridCreator::loadBalance();
+
+ //////////////////////////////////////////////////////////////////////
+ // run the simulation
+ /////////////////////////////////////////////////////////////////////
+
+ // read the initial time step and the end time (mandatory parameters)
+ auto tEnd = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, TimeManager, TEnd);
+ auto dt = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, TimeManager, DtInitial);
+
+ // check if we are about to restart a previously interrupted simulation
+ bool restart = false;
+ Scalar restartTime = 0;
+ if (ParameterTree::tree().hasKey("Restart") || ParameterTree::tree().hasKey("TimeManager.Restart"))
+ {
+ restart = true;
+ restartTime = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, TimeManager, Restart);
+ }
+
+ // instantiate and run the problem
+ TimeManager timeManager;
+ Problem problem(timeManager, GridCreator::bulkGrid().leafGridView(), GridCreator::lowDimGrid().leafGridView());
+ timeManager.init(problem, restartTime, dt, tEnd, restart);
+ timeManager.run();
+
+ // print dumux end message and maybe the parameters for debugging
+ if (mpiHelper.rank() == 0)
+ {
+ DumuxMessage::print(/*firstCall=*/false);
+
+ if (printParams)
+ Parameters::print<TypeTag>();
+ }
+
+ return 0;
+}
+
+/*!
+ * \ingroup Start
+ *
+ * \brief Provides a main function with error handling
+ *
+ * \tparam TypeTag The type tag of the problem which needs to be solved
+ *
+ * \param argc The number of command line arguments of the program
+ * \param argv The contents of the command line arguments of the program
+ * \param usage Callback function for printing the usage message
+ */
+template <class TypeTag>
+int start(int argc,
+ char **argv,
+ void (*usage)(const char *, const std::string &))
+{
+ try {
+ return start_<TypeTag>(argc, argv, usage);
+ }
+ catch (ParameterException &e) {
+ Parameters::print<TypeTag>();
+ std::cerr << std::endl << e << ". Abort!" << std::endl;
+ return 1;
+ }
+ catch (Dune::DGFException & e) {
+ std::cerr << "DGF exception thrown (" << e <<
+ "). Most likely, the DGF file name is wrong "
+ "or the DGF file is corrupted, "
+ "e.g. missing hash at end of file or wrong number (dimensions) of entries."
+ << std::endl;
+ return 2;
+ }
+ catch (Dune::Exception &e) {
+ std::cerr << "Dune reported error: " << e << std::endl;
+ return 3;
+ }
+ catch (...) {
+ std::cerr << "Unknown exception thrown!\n";
+ return 4;
+ }
+}
+
+} // end namespace Dumux
+
+#endif
--
GitLab