diff --git a/dumux/CMakeLists.txt b/dumux/CMakeLists.txt
index 819d019ad5f37d819bd128854c2892abba2684c9..1aae30ec2a946ece00c478cd893d618c9d53c002 100644
--- a/dumux/CMakeLists.txt
+++ b/dumux/CMakeLists.txt
@@ -12,4 +12,5 @@ add_subdirectory(material)
add_subdirectory(multidomain)
add_subdirectory(nonlinear)
add_subdirectory(parallel)
+add_subdirectory(porenetworkflow)
add_subdirectory(porousmediumflow)
diff --git a/dumux/discretization/CMakeLists.txt b/dumux/discretization/CMakeLists.txt
index a43e90fb378a589d1ad6eb942091d1764ad5af9a..6fbdf151abcf8e514b6221fe8a4701f389595d22 100644
--- a/dumux/discretization/CMakeLists.txt
+++ b/dumux/discretization/CMakeLists.txt
@@ -1,6 +1,7 @@
add_subdirectory(box)
add_subdirectory(cellcentered)
add_subdirectory(fem)
+add_subdirectory(porenetwork)
add_subdirectory(projection)
add_subdirectory(staggered)
diff --git a/dumux/discretization/porenetwork/CMakeLists.txt b/dumux/discretization/porenetwork/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..27a356b4c2a9c6a29d6c5b3eb4e32e65f1390d2f
--- /dev/null
+++ b/dumux/discretization/porenetwork/CMakeLists.txt
@@ -0,0 +1,6 @@
+install(FILES
+fvelementgeometry.hh
+gridgeometry.hh
+subcontrolvolume.hh
+subcontrolvolumeface.hh
+DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dumux/discretization/porenetwork)
diff --git a/dumux/discretization/porenetwork/fvelementgeometry.hh b/dumux/discretization/porenetwork/fvelementgeometry.hh
new file mode 100644
index 0000000000000000000000000000000000000000..735fe27e5c5cd7ee951f2ae41ea60b25a65e6968
--- /dev/null
+++ b/dumux/discretization/porenetwork/fvelementgeometry.hh
@@ -0,0 +1,330 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ * See the file COPYING for full copying permissions. *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 3 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see . *
+ *****************************************************************************/
+/*!
+ * \file
+ * \ingroup PoreNetworkDiscretization
+ * \brief Base class for the local geometry for porenetworks
+ */
+#ifndef DUMUX_DISCRETIZATION_PNM_FV_ELEMENT_GEOMETRY_HH
+#define DUMUX_DISCRETIZATION_PNM_FV_ELEMENT_GEOMETRY_HH
+
+#include
+
+#include
+#include
+
+namespace Dumux {
+
+/*!
+ * \ingroup PoreNetworkDiscretization
+ * \brief Base class for the local geometry for porenetworks
+ * \tparam GG the finite volume grid geometry type
+ * \tparam enableFVGridGeometryCache if the grid geometry is cached or not
+ */
+template
+class PNMFVElementGeometry;
+
+//! specialization in case the FVElementGeometries are stored
+template
+class PNMFVElementGeometry
+{
+ using GridView = typename GG::GridView;
+ static constexpr int dim = GridView::dimension;
+ static constexpr int dimWorld = GridView::dimensionworld;
+ using GridIndexType = typename IndexTraits::GridIndex;
+ using LocalIndexType = typename IndexTraits::LocalIndex;
+ using Element = typename GridView::template Codim<0>::Entity;
+ using CoordScalar = typename GridView::ctype;
+ using FeLocalBasis = typename GG::FeCache::FiniteElementType::Traits::LocalBasisType;
+public:
+ //! export type of subcontrol volume
+ using SubControlVolume = typename GG::SubControlVolume;
+ //! export type of subcontrol volume face
+ using SubControlVolumeFace = typename GG::SubControlVolumeFace;
+ //! export type of finite volume grid geometry
+ using GridGeometry = GG;
+ using FVGridGeometry [[deprecated("Use GridGeometry")]] = GG;
+ //! the maximum number of scvs per element
+ static constexpr std::size_t maxNumElementScvs = 2;
+
+ //! Constructor
+ PNMFVElementGeometry(const GridGeometry& gridGeometry)
+ : gridGeometryPtr_(&gridGeometry) {}
+
+ //! Get a sub control volume with a local scv index
+ const SubControlVolume& scv(LocalIndexType scvIdx) const
+ {
+ return gridGeometry().scvs(eIdx_)[scvIdx];
+ }
+
+ //! Get a sub control volume face with a local scvf index
+ const SubControlVolumeFace& scvf(LocalIndexType scvfIdx) const
+ {
+ return gridGeometry().scvfs(eIdx_)[scvfIdx];
+ }
+
+ //! iterator range for sub control volumes. Iterates over
+ //! all scvs of the bound element.
+ //! This is a free function found by means of ADL
+ //! To iterate over all sub control volumes of this FVElementGeometry use
+ //! for (auto&& scv : scvs(fvGeometry))
+ friend inline auto scvs(const PNMFVElementGeometry& fvGeometry)
+ {
+ const auto& g = fvGeometry.gridGeometry();
+ using Iter = typename std::decay_t::const_iterator;
+ return Dune::IteratorRange(g.scvs(fvGeometry.eIdx_).begin(), g.scvs(fvGeometry.eIdx_).end());
+ }
+
+ //! iterator range for sub control volumes faces. Iterates over
+ //! all scvfs of the bound element.
+ //! This is a free function found by means of ADL
+ //! To iterate over all sub control volume faces of this FVElementGeometry use
+ //! for (auto&& scvf : scvfs(fvGeometry))
+ friend inline auto scvfs(const PNMFVElementGeometry& fvGeometry)
+ {
+ const auto& g = fvGeometry.gridGeometry();
+ using Iter = typename std::decay_t::const_iterator;
+ return Dune::IteratorRange(g.scvfs(fvGeometry.eIdx_).begin(), g.scvfs(fvGeometry.eIdx_).end());
+ }
+
+ //! Get a local finite element basis
+ const FeLocalBasis& feLocalBasis() const
+ {
+ return gridGeometry().feCache().get(elementPtr_->geometry().type()).localBasis();
+ }
+
+ //! The total number of sub control volumes
+ std::size_t numScv() const
+ {
+ return gridGeometry().scvs(eIdx_).size();
+ }
+
+ //! The total number of sub control volume faces
+ std::size_t numScvf() const
+ {
+ return gridGeometry().scvfs(eIdx_).size();
+ }
+
+ //! this function is for compatibility reasons with cc methods
+ //! The box stencil is always element-local so bind and bindElement
+ //! are identical.
+ void bind(const Element& element)
+ {
+ this->bindElement(element);
+ }
+
+ //! Binding of an element, has to be called before using the fvgeometries
+ //! Prepares all the volume variables within the element
+ //! For compatibility reasons with the FVGeometry cache being disabled
+ void bindElement(const Element& element)
+ {
+ elementPtr_ = &element;
+ eIdx_ = gridGeometry().elementMapper().index(element);
+ }
+
+ //! The global finite volume geometry we are a restriction of
+ [[deprecated("Use gridGeometry")]]
+ const GridGeometry& fvGridGeometry() const
+ { return *gridGeometryPtr_; }
+
+ //! The global finite volume geometry we are a restriction of
+ const GridGeometry& gridGeometry() const
+ { return *gridGeometryPtr_; }
+
+ //! Returns whether one of the geometry's scvfs lies on a boundary
+ bool hasBoundaryScvf() const
+ { return gridGeometry().hasBoundaryScvf(eIdx_); }
+
+private:
+ const Element* elementPtr_;
+ const GridGeometry* gridGeometryPtr_;
+
+ GridIndexType eIdx_;
+};
+
+//! specialization in case the FVElementGeometries are not stored
+template
+class PNMFVElementGeometry
+{
+ using GridView = typename GG::GridView;
+ static constexpr int dim = GridView::dimension;
+ static constexpr int dimWorld = GridView::dimensionworld;
+ using GridIndexType = typename IndexTraits::GridIndex;
+ using LocalIndexType = typename IndexTraits::LocalIndex;
+ using Element = typename GridView::template Codim<0>::Entity;
+ using CoordScalar = typename GridView::ctype;
+ using FeLocalBasis = typename GG::FeCache::FiniteElementType::Traits::LocalBasisType;
+
+public:
+ //! export type of subcontrol volume
+ using SubControlVolume = typename GG::SubControlVolume;
+ //! export type of subcontrol volume face
+ using SubControlVolumeFace = typename GG::SubControlVolumeFace;
+ //! export type of finite volume grid geometry
+ using GridGeometry = GG;
+ using FVGridGeometry [[deprecated("Use GridGeometry")]] = GG;
+ //! the maximum number of scvs per element
+ static constexpr std::size_t maxNumElementScvs = 2;
+
+ //! Constructor
+ PNMFVElementGeometry(const GridGeometry& gridGeometry)
+ : gridGeometryPtr_(&gridGeometry) {}
+
+ //! Get a sub control volume with a local scv index
+ const SubControlVolume& scv(LocalIndexType scvIdx) const
+ {
+ return scvs_[scvIdx];
+ }
+
+ //! Get a sub control volume face with a local scvf index
+ const SubControlVolumeFace& scvf(LocalIndexType scvfIdx) const
+ {
+ return scvfs_[0];
+ }
+
+ //! iterator range for sub control volumes. Iterates over
+ //! all scvs of the bound element.
+ //! This is a free function found by means of ADL
+ //! To iterate over all sub control volumes of this FVElementGeometry use
+ //! for (auto&& scv : scvs(fvGeometry))
+ friend inline auto scvs(const PNMFVElementGeometry& fvGeometry)
+ {
+ using Iter = typename std::decay_t::const_iterator;
+ return Dune::IteratorRange(fvGeometry.scvs_.begin(), fvGeometry.scvs_.end());
+ }
+
+ //! iterator range for sub control volumes faces. Iterates over
+ //! all scvfs of the bound element.
+ //! This is a free function found by means of ADL
+ //! To iterate over all sub control volume faces of this FVElementGeometry use
+ //! for (auto&& scvf : scvfs(fvGeometry))
+ friend inline auto scvfs(const PNMFVElementGeometry& fvGeometry)
+ {
+ using Iter = typename std::decay_t::const_iterator;
+ return Dune::IteratorRange(fvGeometry.scvfs_.begin(), fvGeometry.scvfs_.end());
+ }
+
+ //! Get a local finite element basis
+ const FeLocalBasis& feLocalBasis() const
+ {
+ return gridGeometry().feCache().get(elementPtr_->geometry().type()).localBasis();
+ }
+
+ //! The total number of sub control volumes
+ std::size_t numScv() const
+ {
+ return scvs_.size();
+ }
+
+ //! The total number of sub control volume faces
+ std::size_t numScvf() const
+ {
+ return scvfs_.size();
+ }
+
+ //! this function is for compatibility reasons with cc methods
+ //! The box stencil is always element-local so bind and bindElement
+ //! are identical.
+ void bind(const Element& element)
+ {
+ this->bindElement(element);
+ }
+
+ //! Binding of an element, has to be called before using the fvgeometries
+ //! Prepares all the volume variables within the element
+ //! For compatibility reasons with the FVGeometry cache being disabled
+ void bindElement(const Element& element)
+ {
+ elementPtr_ = &element;
+ eIdx_ = gridGeometry().elementMapper().index(element);
+ makeElementGeometries(element);
+ }
+
+ //! The global finite volume geometry we are a restriction of
+ [[deprecated("Use gridGeometry")]]
+ const GridGeometry& fvGridGeometry() const
+ { return *gridGeometryPtr_; }
+
+ //! The global finite volume geometry we are a restriction of
+ const GridGeometry& gridGeometry() const
+ { return *gridGeometryPtr_; }
+
+ //! Returns whether one of the geometry's scvfs lies on a boundary
+ bool hasBoundaryScvf() const
+ { return hasBoundaryScvf_; }
+
+private:
+
+ void makeElementGeometries(const Element& element)
+ {
+ hasBoundaryScvf_ = false;
+
+ // get the element geometry
+ auto elementGeometry = element.geometry();
+
+ // construct the sub control volumes
+ for (LocalIndexType scvLocalIdx = 0; scvLocalIdx < elementGeometry.corners(); ++scvLocalIdx)
+ {
+ // get asssociated dof index
+ const auto dofIdxGlobal = gridGeometry().vertexMapper().subIndex(element, scvLocalIdx, dim);
+
+ // get the corners
+ auto corners = std::array{elementGeometry.corner(scvLocalIdx), elementGeometry.center()};
+
+ // get the fractional volume asssociated with the scv
+ const auto volume = gridGeometry().poreVolume(dofIdxGlobal) / gridGeometry().coordinationNumber(dofIdxGlobal);
+
+ // add scv to the local container
+ scvs_[scvLocalIdx] = SubControlVolume(dofIdxGlobal,
+ scvLocalIdx,
+ eIdx_,
+ std::move(corners),
+ volume);
+
+ if (gridGeometry().poreLabel(dofIdxGlobal) > 0)
+ hasBoundaryScvf_ = true;
+ }
+
+ // construct the inner sub control volume face
+ auto unitOuterNormal = elementGeometry.corner(1)-elementGeometry.corner(0);
+ unitOuterNormal /= unitOuterNormal.two_norm();
+ LocalIndexType scvfLocalIdx = 0;
+ scvfs_[0] = SubControlVolumeFace(elementGeometry.center(),
+ std::move(unitOuterNormal),
+ scvfLocalIdx++,
+ std::vector({0, 1}));
+ }
+
+ //! The bound element
+ const Element* elementPtr_;
+ GridIndexType eIdx_;
+
+ //! The global geometry this is a restriction of
+ const GridGeometry* gridGeometryPtr_;
+
+ //! vectors to store the geometries locally after binding an element
+ std::array scvs_;
+ std::array scvfs_;
+
+ bool hasBoundaryScvf_ = false;
+};
+
+} // end namespace Dumux
+
+#endif
diff --git a/dumux/discretization/porenetwork/gridgeometry.hh b/dumux/discretization/porenetwork/gridgeometry.hh
new file mode 100644
index 0000000000000000000000000000000000000000..ff17acacc6ad2f9f1d74a96a0ee3ae9228969f63
--- /dev/null
+++ b/dumux/discretization/porenetwork/gridgeometry.hh
@@ -0,0 +1,802 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ * See the file COPYING for full copying permissions. *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 3 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see . *
+ *****************************************************************************/
+/*!
+ * \file
+ * \ingroup PoreNetworkDiscretization
+ * \brief Base class for the finite volume geometry for porenetwork models
+ */
+#ifndef DUMUX_DISCRETIZATION_PNM_GRID_GEOMETRY_HH
+#define DUMUX_DISCRETIZATION_PNM_GRID_GEOMETRY_HH
+
+#include
+#include
+#include
+
+#include
+#include
+
+#include
+#include
+#include
+#include
+
+#include
+#include
+#include
+#include
+#include
+#include
+#include
+
+namespace Dumux {
+
+/*!
+ * \ingroup PoreNetworkDiscretization
+ * \brief Base class for geometry data extraction from the grid data format
+ */
+template
+class DefaultPNMData
+{
+ using GridIndex = typename IndexTraits::GridIndex;
+ using SmallLocalIndex = typename IndexTraits::SmallLocalIndex;
+ using Label = std::int_least8_t;
+ using Vertex = typename GridView::template Codim::Entity;
+ using Element = typename GridView::template Codim<0>::Entity;
+
+ static const int dim = GridView::dimension;
+
+public:
+
+ template
+ void update(const GridView& gridView, const GridData& gridData)
+ {
+ coordinationNumber_ = gridData.getCoordinationNumbers();
+
+ const auto numThroats = gridView.size(0);
+ throatRadius_.resize(numThroats);
+ throatLength_.resize(numThroats);
+ throatLabel_.resize(numThroats);
+ throatCrossSectionalArea_.resize(numThroats);
+ throatShapeFactor_.resize(numThroats);
+
+ useSameGeometryForAllPores_ = true;
+ useSameShapeForAllThroats_ = true;
+
+ // first check if the same geometry shall be used for all entities ...
+ if (hasParamInGroup(gridData.paramGroup(), "Grid.ThroatCrossSectionShape"))
+ {
+ const auto throatGeometryInput = getParamFromGroup(gridData.paramGroup(), "Grid.ThroatCrossSectionShape");
+ const auto throatGeometry = Throat::shapeFromString(throatGeometryInput);
+ throatGeometry_.resize(1);
+ throatGeometry_[0] = throatGeometry;
+
+ std::cout << "Using '" << throatGeometryInput << "' as cross-sectional shape for all throats." << std::endl;
+ }
+ else // .. otherwise, get the corresponding parameter index from the grid data and resize the respective vector
+ {
+ std::cout << "Reading shape factors for throats from grid data." << std::endl;
+ useSameShapeForAllThroats_ = false;
+ throatGeometry_.resize(numThroats);
+ }
+
+ // get the vertex parameters
+ const auto numPores = gridView.size(dim);
+ poreRadius_.resize(numPores);
+ poreLabel_.resize(numPores);
+ poreVolume_.resize(numPores);
+
+ // first check if the same geometry shall be used for all entities ...
+ if (hasParamInGroup(gridData.paramGroup(), "Grid.PoreGeometry"))
+ {
+ const auto poreGeometryInput = getParamFromGroup(gridData.paramGroup(), "Grid.PoreGeometry");
+ poreGeometry_.resize(1);
+ poreGeometry_[0] = Pore::shapeFromString(poreGeometryInput);;
+
+ std::cout << "Using '" << poreGeometryInput << "' as geometry for all pores." << std::endl;
+ }
+ else // .. otherwise, get the corresponding parameter index from the grid data and resize the respective vector
+ {
+ std::cout << "Reading pore shapes from grid data." << std::endl;
+ useSameGeometryForAllPores_ = false;
+ poreGeometry_.resize(numPores);
+ }
+
+
+ for (const auto& vertex : vertices(gridView))
+ {
+ static const auto poreRadiusIdx = gridData.parameterIndex("PoreRadius");
+ static const auto poreLabelIdx = gridData.parameterIndex("PoreLabel");
+ const auto vIdx = gridView.indexSet().index(vertex);
+ const auto& params = gridData.parameters(vertex);
+ poreRadius_[vIdx] = params[poreRadiusIdx];
+ assert(poreRadius_[vIdx] > 0.0);
+ poreLabel_[vIdx] = params[poreLabelIdx];
+
+ if (!useSameGeometryForAllPores())
+ poreGeometry_[vIdx] = getPoreGeometry_(gridData, vertex);
+
+ poreVolume_[vIdx] = getPoreVolume_(gridData, vertex, vIdx);
+ }
+
+ for (const auto& element : elements(gridView))
+ {
+ const int eIdx = gridView.indexSet().index(element);
+ const auto& params = gridData.parameters(element);
+ static const auto throatRadiusIdx = gridData.parameterIndex("ThroatRadius");
+ static const auto throatLengthIdx = gridData.parameterIndex("ThroatLength");
+ throatRadius_[eIdx] = params[throatRadiusIdx];
+ throatLength_[eIdx] = params[throatLengthIdx];
+
+ // use a default value if no throat label is given by the grid
+ static const bool gridHasThroatLabel = gridData.gridHasElementParameter("ThroatLabel");
+ if (gridHasThroatLabel)
+ {
+ static const auto throatLabelIdx = gridData.parameterIndex("ThroatLabel");
+ throatLabel_[eIdx] = params[throatLabelIdx];
+ }
+ else
+ {
+ const auto vIdx0 = gridView.indexSet().subIndex(element, 0, dim);
+ const auto vIdx1 = gridView.indexSet().subIndex(element, 1, dim);
+
+ const auto poreLabel0 = poreLabel(vIdx0);
+ const auto poreLabel1 = poreLabel(vIdx1);
+
+ if (poreLabel0 >= 0 && poreLabel1 >= 0)
+ {
+ std::cout << "\n Warning: Throat "
+ << eIdx << " connects two boundary pores with different pore labels. Using the greater pore label as throat label.\n"
+ << "Set the throat labels explicitly in your grid file, if needed." << std::endl;
+ }
+
+ using std::max;
+ throatLabel_[eIdx] = max(poreLabel0, poreLabel1);
+ }
+
+ if (!useSameShapeForAllThroats())
+ {
+ static const auto throatShapeFactorIdx = gridData.parameterIndex("ThroatShapeFactor");
+ static const auto throatAreaIdx = gridData.parameterIndex("ThroatCrossSectionalArea");
+ throatShapeFactor_[eIdx] = params[throatShapeFactorIdx];
+ throatGeometry_[eIdx] = Throat::shape(throatShapeFactor_[eIdx]);
+ throatCrossSectionalArea_[eIdx] = params[throatAreaIdx];
+ }
+ else
+ {
+ throatCrossSectionalArea_[eIdx] = getThroatCrossSectionalArea_(gridData, element, eIdx);
+ throatShapeFactor_[eIdx] = getThroatShapeFactor_(gridData, element, eIdx);
+ }
+
+ assert(throatRadius_[eIdx] > 0.0);
+ assert(throatLength_[eIdx] > 0.0);
+ assert(throatCrossSectionalArea_[eIdx] > 0.0);
+
+ static const bool addThroatVolumeToPoreVolume = getParamFromGroup(gridData.paramGroup(), "Grid.AddThroatVolumeToPoreVolume", false);
+ if (addThroatVolumeToPoreVolume)
+ {
+ for (int vIdxLocal = 0; vIdxLocal < 2; ++vIdxLocal)
+ {
+ const auto vIdx = gridView.indexSet().subIndex(element, vIdxLocal, dim);
+ poreVolume_[vIdx] += 0.5 * throatCrossSectionalArea_[eIdx] * throatLength_[eIdx];
+ }
+ }
+ }
+
+ maybeResizeContainers_();
+ }
+
+ //! Returns the pore label (e.g. used for setting BCs)
+ Label poreLabel(const GridIndex dofIdxGlobal) const
+ { return poreLabel_[dofIdxGlobal]; }
+
+ //! Returns the vector of pore labels
+ const std::vector& poreLabel() const
+ { return poreLabel_; }
+
+ //! Returns the radius of the pore
+ Scalar poreRadius(const GridIndex dofIdxGlobal) const
+ { return poreRadius_[dofIdxGlobal]; }
+
+ //! Returns the vector of pore radii
+ const std::vector& poreRadius() const
+ { return poreRadius_; }
+
+ //! Returns the volume of the pore
+ Scalar poreVolume(const GridIndex dofIdxGlobal) const
+ { return poreVolume_[dofIdxGlobal]; }
+
+ //! Returns the vector of pore volumes
+ const std::vector& poreVolume() const
+ { return poreVolume_; }
+
+ //! Returns the radius of the throat
+ Scalar throatRadius(const GridIndex eIdx) const
+ { return throatRadius_[eIdx]; }
+
+ //! Returns the vector of throat radii
+ const std::vector& throatRadius() const
+ { return throatRadius_; }
+
+ //! Returns the length of the throat
+ Scalar throatLength(const GridIndex eIdx) const
+ { return throatLength_[eIdx]; }
+
+ //! Returns the vector of throat lengths
+ const std::vector& throatLength() const
+ { return throatLength_; }
+
+ //! Returns an index indicating if a throat is touching the domain boundary
+ Label throatLabel(const GridIndex eIdx) const
+ { return throatLabel_[eIdx]; }
+
+ //! Returns the vector of throat labels
+ const std::vector& throatLabel() const
+ { return throatLabel_; }
+
+ //! Returns the number of throats connected to a pore (coordination number)
+ SmallLocalIndex coordinationNumber(const GridIndex dofIdxGlobal) const
+ { return coordinationNumber_[dofIdxGlobal]; }
+
+ //! Returns the vector of coordination numbers
+ const std::vector& coordinationNumber() const
+ { return coordinationNumber_; }
+
+ //! the geometry of the pore
+ Pore::Shape poreGeometry(const GridIndex vIdx) const
+ { return useSameGeometryForAllPores() ? poreGeometry_[0] : poreGeometry_[vIdx]; }
+
+ //! Returns the vector of pore geometries
+ const std::vector& poreGeometry() const
+ {
+ if (useSameGeometryForAllPores())
+ {
+ // if a vector of pore geometries is requested (e.g., for vtk output),
+ // resize the container and fill it with the same value everywhere
+ const auto poreGeo = poreGeometry_[0];
+ poreGeometry_.resize(poreRadius_.size(), poreGeo);
+ }
+
+ return poreGeometry_;
+ }
+
+ //! Returns the throat's cross-sectional shape
+ Throat::Shape throatCrossSectionShape(const GridIndex eIdx) const
+ { return useSameShapeForAllThroats() ? throatGeometry_[0] : throatGeometry_[eIdx]; }
+
+ //! Returns the vector of cross-sectional shapes
+ const std::vector& throatCrossSectionShape() const
+ {
+ if (useSameShapeForAllThroats())
+ {
+ // if a vector of throat cross section shapes is requested (e.g., for vtk output),
+ // resize the container and fill it with the same value everywhere
+ const auto throatShape = throatGeometry_[0];
+ throatGeometry_.resize(throatRadius_.size(), throatShape);
+ }
+
+ return throatGeometry_;
+ }
+
+ //! Returns the throat's cross-sectional area
+ Scalar throatCrossSectionalArea(const GridIndex eIdx) const
+ { return throatCrossSectionalArea_[eIdx]; }
+
+ //! Returns the vector of throat cross-sectional areas
+ const std::vector& throatCrossSectionalArea() const
+ { return throatCrossSectionalArea_; }
+
+ //! Returns the throat's shape factor
+ Scalar throatShapeFactor(const GridIndex eIdx) const
+ { return useSameShapeForAllThroats() ? throatShapeFactor_[0] : throatShapeFactor_[eIdx]; }
+
+ //! Returns the vector of throat shape factors
+ const std::vector& throatShapeFactor() const
+ {
+ if (useSameShapeForAllThroats())
+ {
+ // if a vector of throat shape factors is requested (e.g., for vtk output),
+ // resize the container and fill it with the same value everywhere
+ const auto shapeFactor = throatShapeFactor_[0];
+ throatShapeFactor_.resize(throatRadius_.size(), shapeFactor);
+ }
+
+ return throatShapeFactor_;
+ }
+
+ //! Returns whether all pores feature the same shape
+ bool useSameGeometryForAllPores() const
+ { return useSameGeometryForAllPores_; }
+
+ //! Returns whether all throats feature the same cross-sectional shape
+ bool useSameShapeForAllThroats() const
+ { return useSameShapeForAllThroats_; }
+
+private:
+
+ //! determine the pore geometry provided as scalar value by the grid file
+ template
+ Pore::Shape getPoreGeometry_(const GridData& gridData, const Vertex& vertex) const
+ {
+ static const auto poreGeometryIdx = gridData.parameterIndex("PoreGeometry");
+ using T = std::underlying_type_t;
+ const auto poreGeometryValue = static_cast(gridData.parameters(vertex)[poreGeometryIdx]);
+ return static_cast(poreGeometryValue);
+ }
+
+ //! automatically determine the pore volume if not provided by the grid file
+ template
+ Scalar getPoreVolume_(const GridData& gridData, const Vertex& vertex, const std::size_t vIdx) const
+ {
+ static const bool gridHasPoreVolume = gridData.gridHasVertexParameter("PoreVolume");
+
+ if (gridHasPoreVolume)
+ {
+ static const auto poreVolumeIdx = gridData.parameterIndex("PoreVolume");
+ return gridData.parameters(vertex)[poreVolumeIdx];
+ }
+ else
+ {
+ if (poreGeometry(vIdx) == Pore::Shape::cylinder)
+ {
+ static const Scalar fixedHeight = getParamFromGroup(gridData.paramGroup(), "Grid.PoreHeight", -1.0);
+ const Scalar h = fixedHeight > 0.0 ? fixedHeight : gridData.getParameter(vertex, "PoreHeight");
+ return Pore::volume(Pore::Shape::cylinder, poreRadius(vIdx), h);
+ }
+ else
+ return Pore::volume(poreGeometry(vIdx), poreRadius(vIdx));
+ }
+ }
+
+ //! automatically determine throat cross-sectional area if not provided by the grid file
+ template
+ Scalar getThroatCrossSectionalArea_(const GridData& gridData, const Element& element, const std::size_t eIdx) const
+ {
+ static const bool gridHasThroatCrossSectionalArea = gridData.gridHasVertexParameter("ThroatCrossSectionalArea");
+ if (gridHasThroatCrossSectionalArea)
+ {
+ static const auto throatAreaIdx = gridData.parameterIndex("ThroatCrossSectionalArea");
+ return gridData.parameters(element)[throatAreaIdx];
+ }
+ else
+ {
+ if (const auto shape = throatCrossSectionShape(eIdx); shape == Throat::Shape::rectangle)
+ {
+ static const auto throatHeight = getParamFromGroup(gridData.paramGroup(), "Grid.ThroatHeight");
+ return Throat::totalCrossSectionalAreaForRectangle(throatRadius_[eIdx], throatHeight);
+ }
+ else
+ return Throat::totalCrossSectionalArea(shape, throatRadius_[eIdx]);
+ }
+ }
+
+ //! automatically determine throat shape factor if not provided by the grid file
+ template
+ Scalar getThroatShapeFactor_(const GridData& gridData, const Element& element, const std::size_t eIdx) const
+ {
+ static const bool gridHasThroatShapeFactor = gridData.gridHasVertexParameter("ThroatShapeFactor");
+ if (gridHasThroatShapeFactor)
+ {
+ static const auto throatShapeFactorIdx = gridData.parameterIndex("ThroatShapeFactor");
+ return gridData.parameters(element)[throatShapeFactorIdx];
+ }
+ else
+ {
+ if (const auto shape = throatCrossSectionShape(eIdx); shape == Throat::Shape::rectangle)
+ {
+ static const auto throatHeight = getParamFromGroup(gridData.paramGroup(), "Grid.ThroatHeight");
+ return Throat::shapeFactorRectangle(throatRadius_[eIdx], throatHeight);
+ }
+ else if (shape == Throat::Shape::polygon || shape == Throat::Shape::scaleneTriangle)
+ {
+ static const auto shapeFactor = getParamFromGroup(gridData.paramGroup(), "Grid.ThroatShapeFactor");
+ return shapeFactor;
+ }
+ else
+ return Throat::shapeFactor(shape, throatRadius_[eIdx]);
+ }
+ }
+
+ void maybeResizeContainers_()
+ {
+ // check if all throat might have the same shape in order to save some memory
+ if (!useSameShapeForAllThroats() &&
+ std::adjacent_find(throatGeometry_.begin(), throatGeometry_.end(), std::not_equal_to() ) == throatGeometry_.end())
+ {
+ std::cout << "All throats feature the same shape, resizing containers" << std::endl;
+ useSameShapeForAllThroats_ = true;
+ const Scalar shapeFactor = throatShapeFactor_[0];
+ const auto throatGeometry = throatGeometry_[0];
+ throatShapeFactor_.resize(1);
+ throatGeometry_.resize(1);
+ throatShapeFactor_[0] = shapeFactor;
+ throatGeometry_[0] = throatGeometry;
+ }
+
+ // check if all throat might have the same shape in order to save some memory
+ if (!useSameGeometryForAllPores() &&
+ std::adjacent_find(poreGeometry_.begin(), poreGeometry_.end(), std::not_equal_to() ) == poreGeometry_.end())
+ {
+ std::cout << "All pores feature the same shape, resizing containers" << std::endl;
+ useSameGeometryForAllPores_ = true;
+ const auto poreGeometry = poreGeometry_[0];
+ poreGeometry_.resize(1);
+ poreGeometry_[0] = poreGeometry;
+ }
+ }
+
+ mutable std::vector poreGeometry_;
+ std::vector poreRadius_;
+ std::vector poreVolume_;
+ std::vector poreLabel_; // 0:no, 1:general, 2:coupling1, 3:coupling2, 4:inlet, 5:outlet
+ std::vector coordinationNumber_;
+ mutable std::vector throatGeometry_;
+ mutable std::vector throatShapeFactor_;
+ std::vector throatRadius_;
+ std::vector throatLength_;
+ std::vector throatLabel_; // 0:no, 1:general, 2:coupling1, 3:coupling2, 4:inlet, 5:outlet
+ std::vector throatCrossSectionalArea_;
+ bool useSameGeometryForAllPores_;
+ bool useSameShapeForAllThroats_;
+};
+
+/*!
+ * \ingroup PoreNetworkDiscretization
+ * \brief The default traits
+ * \tparam the grid view type
+ */
+template>
+struct PNMDefaultGridGeometryTraits
+: public MapperTraits
+{
+ using SubControlVolume = PNMSubControlVolume;
+ using SubControlVolumeFace = PNMSubControlVolumeFace;
+
+ template
+ using LocalView = PNMFVElementGeometry;
+
+ using PNMData = DefaultPNMData;
+};
+
+/*!
+ * \ingroup PoreNetworkDiscretization
+ * \brief Base class for the finite volume geometry for porenetwork models
+ * \note This class is specialized for versions with and without caching the fv geometries on the grid view
+ */
+template >
+class PNMGridGeometry;
+
+/*!
+ * \ingroup PoreNetworkDiscretization
+ * \brief Base class for the finite volume geometry for porenetwork models
+ * \note For caching enabled we store the fv geometries for the whole grid view which is memory intensive but faster
+ */
+template
+class PNMGridGeometry
+: public BaseGridGeometry
+, public Traits::PNMData
+{
+ using ThisType = PNMGridGeometry;
+ using ParentType = BaseGridGeometry;
+ using GridIndexType = typename IndexTraits::GridIndex;
+ using LocalIndexType = typename IndexTraits::LocalIndex;
+ using PNMData = typename Traits::PNMData;
+
+ using Element = typename GV::template Codim<0>::Entity;
+ using CoordScalar = typename GV::ctype;
+ static const int dim = GV::dimension;
+ static const int dimWorld = GV::dimensionworld;
+
+public:
+ //! export discretization method
+ static constexpr DiscretizationMethod discMethod = DiscretizationMethod::box;
+
+ //! export the type of the fv element geometry (the local view type)
+ using LocalView = typename Traits::template LocalView;
+ //! export the type of sub control volume
+ using SubControlVolume = typename Traits::SubControlVolume;
+ //! export the type of sub control volume
+ using SubControlVolumeFace = typename Traits::SubControlVolumeFace;
+ //! export the type of extrusion
+ using Extrusion = Extrusion_t;
+ //! export dof mapper type
+ using DofMapper = typename Traits::VertexMapper;
+ //! export the finite element cache type
+ using FeCache = Dune::PQkLocalFiniteElementCache;
+ //! export the grid view type
+ using GridView = GV;
+
+ //! Constructor
+ PNMGridGeometry(const GridView gridView)
+ : ParentType(gridView)
+ {
+ static_assert(GridView::dimension == 1, "Porenetwork model only allow GridView::dimension == 1!");
+ }
+
+ //! the vertex mapper is the dofMapper
+ //! this is convenience to have better chance to have the same main files for box/tpfa/mpfa...
+ const DofMapper& dofMapper() const
+ { return this->vertexMapper(); }
+
+ //! The total number of sub control volumes
+ std::size_t numScv() const
+ { return numScv_; }
+
+ //! The total number of sun control volume faces
+ std::size_t numScvf() const
+ { return numScvf_; }
+
+ //! The total number of boundary sub control volume faces
+ //! For compatibility reasons with cc methods
+ std::size_t numBoundaryScvf() const
+ { return 0; }
+
+ //! The total number of degrees of freedom
+ std::size_t numDofs() const
+ { return this->vertexMapper().size(); }
+
+ //! update all fvElementGeometries (do this again after grid adaption)
+ template
+ void update(const GridData& gridData)
+ {
+ ParentType::update();
+ PNMData::update(this->gridView(), gridData);
+
+ scvs_.clear();
+ scvfs_.clear();
+
+ auto numElements = this->gridView().size(0);
+ scvs_.resize(numElements);
+ scvfs_.resize(numElements);
+ hasBoundaryScvf_.resize(numElements, false);
+
+ boundaryDofIndices_.assign(numDofs(), false);
+
+ numScvf_ = numElements;
+ numScv_ = 2*numElements;
+
+ // Build the SCV and SCV faces
+ for (const auto& element : elements(this->gridView()))
+ {
+ // get the element geometry
+ auto eIdx = this->elementMapper().index(element);
+ auto elementGeometry = element.geometry();
+
+ // construct the sub control volumes
+ for (LocalIndexType scvLocalIdx = 0; scvLocalIdx < elementGeometry.corners(); ++scvLocalIdx)
+ {
+ const auto dofIdxGlobal = this->vertexMapper().subIndex(element, scvLocalIdx, dim);
+
+ // get the corners
+ auto corners = std::array{elementGeometry.corner(scvLocalIdx), elementGeometry.center()};
+
+ // get the fractional volume asssociated with the scv
+ const auto volume = this->poreVolume(dofIdxGlobal) / this->coordinationNumber(dofIdxGlobal);
+
+ scvs_[eIdx][scvLocalIdx] = SubControlVolume(dofIdxGlobal,
+ scvLocalIdx,
+ eIdx,
+ std::move(corners),
+ volume);
+
+ if (this->poreLabel(dofIdxGlobal) > 0)
+ {
+ if (boundaryDofIndices_[dofIdxGlobal])
+ continue;
+
+ boundaryDofIndices_[dofIdxGlobal] = true;
+ hasBoundaryScvf_[eIdx] = true;
+ }
+ }
+
+ // construct the inner sub control volume face
+ auto unitOuterNormal = elementGeometry.corner(1)-elementGeometry.corner(0);
+ unitOuterNormal /= unitOuterNormal.two_norm();
+ LocalIndexType scvfLocalIdx = 0;
+ scvfs_[eIdx][0] = SubControlVolumeFace(elementGeometry.center(),
+ std::move(unitOuterNormal),
+ scvfLocalIdx++,
+ std::vector({0, 1}));
+ }
+ }
+
+ //! The finite element cache for creating local FE bases
+ const FeCache& feCache() const
+ { return feCache_; }
+
+ //! Get the local scvs for an element
+ const std::array& scvs(GridIndexType eIdx) const
+ { return scvs_[eIdx]; }
+
+ //! Get the local scvfs for an element
+ const std::array& scvfs(GridIndexType eIdx) const
+ { return scvfs_[eIdx]; }
+
+ //! If a vertex / d.o.f. is on the boundary
+ bool dofOnBoundary(GridIndexType dofIdx) const
+ { return boundaryDofIndices_[dofIdx]; }
+
+ //! If a vertex / d.o.f. is on a periodic boundary (not implemented)
+ bool dofOnPeriodicBoundary(GridIndexType dofIdx) const
+ { return false; }
+
+ //! The index of the vertex / d.o.f. on the other side of the periodic boundary
+ GridIndexType periodicallyMappedDof(GridIndexType dofIdx) const
+ { DUNE_THROW(Dune::NotImplemented, "Periodic boundaries"); }
+
+ //! Returns the map between dofs across periodic boundaries
+ std::unordered_map periodicVertexMap() const
+ { return std::unordered_map{}; }
+
+ //! Returns whether one of the geometry's scvfs lies on a boundary
+ bool hasBoundaryScvf(GridIndexType eIdx) const
+ { return hasBoundaryScvf_[eIdx]; }
+
+private:
+ const FeCache feCache_;
+
+ std::vector> scvs_;
+ std::vector> scvfs_;
+
+ std::size_t numScv_;
+ std::size_t numScvf_;
+
+ // vertices on the boudary
+ std::vector boundaryDofIndices_;
+ std::vector hasBoundaryScvf_;
+};
+
+/*!
+ * \ingroup PoreNetworkDiscretization
+ * \brief Base class for the finite volume geometry for porenetwork models
+ * \note For caching disabled we store only some essential index maps to build up local systems on-demand in
+ * the corresponding FVElementGeometry
+ */
+template
+class PNMGridGeometry
+: public BaseGridGeometry
+, public Traits::PNMData
+{
+ using ThisType = PNMGridGeometry;
+ using ParentType = BaseGridGeometry;
+ using GridIndexType = typename IndexTraits::GridIndex;
+ using LocalIndexType = typename IndexTraits::LocalIndex;
+ using PNMData = typename Traits::PNMData;
+
+ static const int dim = GV::dimension;
+ static const int dimWorld = GV::dimensionworld;
+
+ using Element = typename GV::template Codim<0>::Entity;
+ using CoordScalar = typename GV::ctype;
+
+public:
+ //! export discretization method
+ static constexpr DiscretizationMethod discMethod = DiscretizationMethod::box;
+
+ //! export the type of the fv element geometry (the local view type)
+ using LocalView = typename Traits::template LocalView;
+ //! export the type of sub control volume
+ using SubControlVolume = typename Traits::SubControlVolume;
+ //! export the type of sub control volume
+ using SubControlVolumeFace = typename Traits::SubControlVolumeFace;
+ //! export the type of extrusion
+ using Extrusion = Extrusion_t;
+ //! export dof mapper type
+ using DofMapper = typename Traits::VertexMapper;
+ //! export the finite element cache type
+ using FeCache = Dune::PQkLocalFiniteElementCache;
+ //! export the grid view type
+ using GridView = GV;
+
+ //! Constructor
+ PNMGridGeometry(const GridView gridView)
+ : ParentType(gridView)
+ {
+ static_assert(GridView::dimension == 1, "Porenetwork model only allow GridView::dimension == 1!");
+ }
+
+ //! the vertex mapper is the dofMapper
+ //! this is convenience to have better chance to have the same main files for box/tpfa/mpfa...
+ const DofMapper& dofMapper() const
+ { return this->vertexMapper(); }
+
+ //! The total number of sub control volumes
+ std::size_t numScv() const
+ { return numScv_; }
+
+ //! The total number of sun control volume faces
+ std::size_t numScvf() const
+ { return numScvf_; }
+
+ //! The total number of boundary sub control volume faces
+ //! For compatibility reasons with cc methods
+ std::size_t numBoundaryScvf() const
+ { return 0; }
+
+ //! The total number of degrees of freedom
+ std::size_t numDofs() const
+ { return this->vertexMapper().size(); }
+
+ //! update all fvElementGeometries (do this again after grid adaption)
+ template
+ void update(const GridData& gridData)
+ {
+ ParentType::update();
+ PNMData::update(this->gridView(), gridData);
+
+ boundaryDofIndices_.assign(numDofs(), false);
+
+ // save global data on the grid's scvs and scvfs
+ numScvf_ = this->gridView().size(0);
+ numScv_ = 2*numScvf_;
+
+ for (const auto& element : elements(this->gridView()))
+ {
+ // treat boundaries
+ for (LocalIndexType vIdxLocal = 0; vIdxLocal < 2; ++vIdxLocal)
+ {
+ const auto vIdxGlobal = this->vertexMapper().subIndex(element, vIdxLocal, dim);
+ if (this->poreLabel(vIdxGlobal) > 0)
+ {
+ if (boundaryDofIndices_[vIdxGlobal])
+ continue;
+
+ boundaryDofIndices_[vIdxGlobal] = true;
+ }
+ }
+ }
+ }
+
+ //! The finite element cache for creating local FE bases
+ const FeCache& feCache() const
+ { return feCache_; }
+
+ //! If a vertex / d.o.f. is on the boundary
+ bool dofOnBoundary(GridIndexType dofIdx) const
+ { return boundaryDofIndices_[dofIdx]; }
+
+ //! If a vertex / d.o.f. is on a periodic boundary (not implemented)
+ bool dofOnPeriodicBoundary(GridIndexType dofIdx) const
+ { return false; }
+
+ //! The index of the vertex / d.o.f. on the other side of the periodic boundary
+ GridIndexType periodicallyMappedDof(GridIndexType dofIdx) const
+ { DUNE_THROW(Dune::NotImplemented, "Periodic boundaries"); }
+
+ //! Returns the map between dofs across periodic boundaries
+ std::unordered_map periodicVertexMap() const
+ { return std::unordered_map{}; }
+
+private:
+
+ const FeCache feCache_;
+
+ // Information on the global number of geometries
+ std::size_t numScv_;
+ std::size_t numScvf_;
+
+ // vertices on the boudary
+ std::vector boundaryDofIndices_;
+};
+
+} // end namespace Dumux
+
+#endif
diff --git a/dumux/discretization/porenetwork/subcontrolvolume.hh b/dumux/discretization/porenetwork/subcontrolvolume.hh
new file mode 100644
index 0000000000000000000000000000000000000000..be7bf270c5f7e3115733cc0e1eb4c81e9e46f8bc
--- /dev/null
+++ b/dumux/discretization/porenetwork/subcontrolvolume.hh
@@ -0,0 +1,150 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ * See the file COPYING for full copying permissions. *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 3 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see . *
+ *****************************************************************************/
+/*!
+ * \file
+ * \ingroup PoreNetworkDiscretization
+ * \brief the sub control volume for porenetworks
+ */
+#ifndef DUMUX_DISCRETIZATION_PNM_SUBCONTROLVOLUME_HH
+#define DUMUX_DISCRETIZATION_PNM_SUBCONTROLVOLUME_HH
+
+#include
+#include
+#include
+#include
+
+namespace Dumux {
+
+/*!
+ * \ingroup PoreNetworkDiscretization
+ * \brief Default traits class
+ * \tparam GV the type of the grid view
+ */
+template
+struct PNMDefaultScvGeometryTraits
+{
+ using Grid = typename GridView::Grid;
+
+ static const int dim = Grid::dimension;
+ static const int dimWorld = Grid::dimensionworld;
+
+ using GridIndexType = typename IndexTraits::GridIndex;
+ using LocalIndexType = typename IndexTraits::LocalIndex;
+ using Scalar = typename Grid::ctype;
+ using GlobalPosition = Dune::FieldVector;
+ using CornerStorage = std::array;
+ using Geometry = Dune::AffineGeometry;
+};
+
+/*!
+ * \ingroup PoreNetworkDiscretization
+ * \brief the sub control volume for porenetworks
+ * \tparam GV the type of the grid view
+ * \tparam T the scv geometry traits
+ */
+template >
+class PNMSubControlVolume
+: public SubControlVolumeBase, T>
+{
+ using ThisType = PNMSubControlVolume;
+ using ParentType = SubControlVolumeBase;
+ using GridIndexType = typename T::GridIndexType;
+ using LocalIndexType = typename T::LocalIndexType;
+ using Scalar = typename T::Scalar;
+ using CornerStorage = typename T::CornerStorage;
+ using Geometry = typename T::Geometry;
+
+public:
+ //! export the type used for global coordinates
+ using GlobalPosition = typename T::GlobalPosition;
+ //! state the traits public and thus export all types
+ using Traits = T;
+
+ //! The default constructor
+ PNMSubControlVolume() = default;
+
+ // the contructor in the box case
+ template
+ PNMSubControlVolume(GridIndexType dofIndex,
+ LocalIndexType scvIdx,
+ GridIndexType elementIndex,
+ Corners&& corners,
+ const Scalar volume)
+ : center_((corners[0]+corners[1])/2.0),
+ corners_(std::forward(corners)),
+ volume_(volume),
+ elementIndex_(elementIndex),
+ localDofIdx_(scvIdx),
+ dofIndex_(dofIndex)
+ {}
+
+ //! The center of the sub control volume (return pore center)
+ const GlobalPosition& center() const
+ { return center_; }
+
+ //! The volume of the sub control volume (part of a pore)
+ Scalar volume() const
+ { return volume_; }
+
+ //! The element-local index of the dof this scv is embedded in
+ LocalIndexType localDofIndex() const
+ { return localDofIdx_; }
+
+ //! The element-local index of this scv.
+ LocalIndexType indexInElement() const
+ { return localDofIdx_; }
+
+ //! The index of the dof this scv is embedded in
+ GridIndexType dofIndex() const
+ { return dofIndex_; }
+
+ // The position of the dof this scv is embedded in
+ const GlobalPosition& dofPosition() const
+ { return corners_[0]; }
+
+ //! The global index of the element this scv is embedded in
+ GridIndexType elementIndex() const
+ { return elementIndex_; }
+
+ //! Return the corner for the given local index
+ const GlobalPosition& corner(LocalIndexType localIdx) const
+ {
+ assert(localIdx < corners_.size() && "provided index exceeds the number of corners");
+ return corners_[localIdx];
+ }
+
+ //! The geometry of the sub control volume e.g. for integration
+ Geometry geometry() const
+ {
+ return Geometry(Dune::GeometryTypes::simplex(1), corners_);
+ }
+
+private:
+ GlobalPosition center_;
+ CornerStorage corners_;
+ Scalar volume_;
+ GridIndexType elementIndex_;
+ LocalIndexType localDofIdx_;
+ GridIndexType dofIndex_;
+};
+
+} // end namespace Dumux
+
+#endif
diff --git a/dumux/discretization/porenetwork/subcontrolvolumeface.hh b/dumux/discretization/porenetwork/subcontrolvolumeface.hh
new file mode 100644
index 0000000000000000000000000000000000000000..72daef3a0d46926cc41f2545f25bd34b7751ff20
--- /dev/null
+++ b/dumux/discretization/porenetwork/subcontrolvolumeface.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 3 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see . *
+ *****************************************************************************/
+/*!
+ * \file
+ * \ingroup PoreNetworkDiscretization
+ * \brief Base class for a sub control volume face
+ */
+#ifndef DUMUX_DISCRETIZATION_PNM_SUBCONTROLVOLUMEFACE_HH
+#define DUMUX_DISCRETIZATION_PNM_SUBCONTROLVOLUMEFACE_HH
+
+#include
+#include
+
+namespace Dumux {
+
+/*!
+ * \ingroup PoreNetworkDiscretization
+ * \brief Default traits class
+ * \tparam GV the type of the grid view
+ */
+template
+struct PNMDefaultScvfGeometryTraits
+{
+ using Grid = typename GridView::Grid;
+ static constexpr int dim = Grid::dimension;
+ static constexpr int dimWorld = Grid::dimensionworld;
+
+ using GridIndexType = typename IndexTraits::GridIndex;
+ using LocalIndexType = typename IndexTraits::LocalIndex;
+
+ using Scalar = typename Grid::ctype;
+ using GlobalPosition = Dune::FieldVector;
+ using CornerStorage = std::array, 1>;
+};
+
+/*!
+ * \ingroup PoreNetworkDiscretization
+ * \brief Class for a sub control volume face for porenetworks
+ * \tparam GV the type of the grid view
+ * \tparam T the scvf geometry traits
+ */
+template >
+class PNMSubControlVolumeFace
+: public SubControlVolumeFaceBase, T>
+{
+ using ThisType = PNMSubControlVolumeFace;
+ using ParentType = SubControlVolumeFaceBase;
+ using GridIndexType = typename T::GridIndexType;
+ using LocalIndexType = typename T::LocalIndexType;
+ using Scalar = typename T::Scalar;
+
+public:
+ //! export the type used for global coordinates
+ using GlobalPosition = typename T::GlobalPosition;
+ //! state the traits public and thus export all types
+ using Traits = T;
+
+ //! The default constructor
+ PNMSubControlVolumeFace() = default;
+
+ //! Constructor for inner scvfs
+ PNMSubControlVolumeFace(const GlobalPosition& center,
+ const GlobalPosition& unitOuterNormal,
+ GridIndexType scvfIndex,
+ std::vector&& scvIndices)
+ : center_(center),
+ unitOuterNormal_(unitOuterNormal),
+ area_(1.0), // TODO throat cross-section area
+ scvfIndex_(scvfIndex),
+ scvIndices_(std::move(scvIndices))
+ {}
+
+ //! The center of the sub control volume face
+ const GlobalPosition& center() const
+ {
+ return center_;
+ }
+
+ //! The integration point for flux evaluations in global coordinates
+ const GlobalPosition& ipGlobal() const
+ {
+ return center_;
+ }
+
+ //! The area of the sub control volume face
+ Scalar area() const
+ {
+ return area_;
+ }
+
+ //! We assume to always have a pore body and not a pore throat at the boundary.
+ bool boundary() const
+ {
+ return false;
+ }
+
+ const GlobalPosition& unitOuterNormal() const
+ {
+ return unitOuterNormal_;
+ }
+
+ //! index of the inside sub control volume for spatial param evaluation
+ LocalIndexType insideScvIdx() const
+ {
+ return scvIndices_[0];
+ }
+
+ //! index of the outside sub control volume for spatial param evaluation
+ // This results in undefined behaviour if boundary is true
+ LocalIndexType outsideScvIdx() const
+ {
+ assert(!boundary());
+ return scvIndices_[1];
+ }
+
+ //! The local index of this sub control volume face
+ GridIndexType index() const
+ {
+ return scvfIndex_;
+ }
+
+private:
+ GlobalPosition center_;
+ GlobalPosition unitOuterNormal_;
+ Scalar area_;
+ GridIndexType scvfIndex_;
+ std::vector scvIndices_;
+};
+
+} // end namespace Dumux
+
+#endif
diff --git a/dumux/flux/CMakeLists.txt b/dumux/flux/CMakeLists.txt
index 6ec79ed53dd6e07e0553b3f76a1110ff88534cf6..76b85c4062d2e079da6bdad8aa441d944dba9d4b 100644
--- a/dumux/flux/CMakeLists.txt
+++ b/dumux/flux/CMakeLists.txt
@@ -1,6 +1,7 @@
add_subdirectory(box)
add_subdirectory(ccmpfa)
add_subdirectory(cctpfa)
+add_subdirectory(porenetwork)
add_subdirectory(shallowwater)
add_subdirectory(staggered)
diff --git a/dumux/flux/porenetwork/CMakeLists.txt b/dumux/flux/porenetwork/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..bb91055eeca8076bf7ca01c912e9d91f2af28b2d
--- /dev/null
+++ b/dumux/flux/porenetwork/CMakeLists.txt
@@ -0,0 +1,5 @@
+install(FILES
+advection.hh
+fickslaw.hh
+fourierslaw.hh
+DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dumux/flux/porenetwork)
diff --git a/dumux/flux/porenetwork/advection.hh b/dumux/flux/porenetwork/advection.hh
new file mode 100644
index 0000000000000000000000000000000000000000..154a5f23ce9a9b987a1e2b66669d7bfd111ce912
--- /dev/null
+++ b/dumux/flux/porenetwork/advection.hh
@@ -0,0 +1,154 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ * See the file COPYING for full copying permissions. *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 3 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see . *
+ *****************************************************************************/
+/*!
+ * \file
+ * \brief This file contains the data which is required to calculate
+ * the fluxes of the pore network model over a face of a finite volume.
+ */
+#ifndef DUMUX_FLUX_PNM_ADVECTION_HH
+#define DUMUX_FLUX_PNM_ADVECTION_HH
+
+namespace Dumux
+{
+
+namespace Detail {
+
+template
+struct Transmissibility : public TransmissibilityLawTypes...
+{};
+
+}
+
+/*!
+ * \file
+ * \ingroup PoreNetworkFlux
+ * \brief Hagen–Poiseuille-type flux law to describe the advective flux for pore-network models.
+ */
+template
+class PoreNetworkCreepingFlow
+{
+
+public:
+ //! Export the Scalar type
+ using Scalar = ScalarT;
+
+ //! Export the transmissibility law
+ using Transmissibility = Detail::Transmissibility;
+
+ /*!
+ * \brief Returns the advective flux of a fluid phase
+ * across the given sub-control volume face (corresponding to a pore throat).
+ * \note The flux is given in N*m, and can be converted
+ * into a volume flux (m^3/s) or mass flux (kg/s) by applying an upwind scheme
+ * for the mobility (1/viscosity) or the product of density and mobility, respectively.
+ */
+ template
+ static Scalar flux(const Problem& problem,
+ const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars,
+ const SubControlVolumeFace& scvf,
+ const int phaseIdx,
+ const ElemFluxVarsCache& elemFluxVarsCache)
+ {
+ const auto& fluxVarsCache = elemFluxVarsCache[scvf];
+
+ // Get the inside and outside volume variables
+ const auto& insideScv = fvGeometry.scv(scvf.insideScvIdx());
+ const auto& outsideScv = fvGeometry.scv(scvf.outsideScvIdx());
+ const auto& insideVolVars = elemVolVars[insideScv];
+ const auto& outsideVolVars = elemVolVars[outsideScv];
+
+ // calculate the pressure difference
+ const Scalar deltaP = insideVolVars.pressure(phaseIdx) - outsideVolVars.pressure(phaseIdx);
+ const Scalar transmissibility = fluxVarsCache.transmissibility(phaseIdx);
+
+ Scalar volumeFlow = transmissibility*deltaP;
+
+ // add gravity term
+ static const bool enableGravity = getParamFromGroup(problem.paramGroup(), "Problem.EnableGravity");
+ if (enableGravity)
+ {
+ const Scalar rho = 0.5*insideVolVars.density(phaseIdx) + 0.5*outsideVolVars.density(phaseIdx);
+ const Scalar g = problem.spatialParams().gravity(scvf.center()) * scvf.unitOuterNormal();
+
+ // The transmissibility is with respect to the effective throat length (potentially dropping the pore body radii).
+ // For gravity, we need to consider the total throat length (i.e., the cell-center to cell-center distance).
+ // This might cause some inconsistencies TODO: is there a better way?
+ volumeFlow += transmissibility * element.geometry().volume() * rho * g;
+ }
+
+ return volumeFlow;
+ }
+
+ /*!
+ * \brief Returns the throat conductivity
+ */
+ template
+ static Scalar calculateTransmissibility(const Problem& problem,
+ const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const typename FVElementGeometry::SubControlVolumeFace& scvf,
+ const ElementVolumeVariables& elemVolVars,
+ const FluxVariablesCache& fluxVarsCache,
+ const int phaseIdx)
+ {
+ if constexpr (ElementVolumeVariables::VolumeVariables::numFluidPhases() == 1)
+ return Transmissibility::singlePhaseTransmissibility(problem, element, fvGeometry, scvf, elemVolVars, fluxVarsCache, phaseIdx);
+ else
+ {
+ static_assert(ElementVolumeVariables::VolumeVariables::numFluidPhases() == 2);
+
+ const auto& spatialParams = problem.spatialParams();
+ using FluidSystem = typename ElementVolumeVariables::VolumeVariables::FluidSystem;
+ const int wPhaseIdx = spatialParams.template wettingPhase(element, elemVolVars);
+ const bool invaded = fluxVarsCache.invaded();
+
+ if (phaseIdx == wPhaseIdx)
+ {
+ return invaded ? Transmissibility::wettingLayerTransmissibility(element, fvGeometry, scvf, fluxVarsCache)
+ : Transmissibility::singlePhaseTransmissibility(problem, element, fvGeometry, scvf, elemVolVars, fluxVarsCache, phaseIdx);
+ }
+ else // non-wetting phase
+ {
+ return invaded ? Transmissibility::nonWettingPhaseTransmissibility(element, fvGeometry, scvf, fluxVarsCache)
+ : 0.0;
+ }
+ }
+ }
+
+ template
+ static std::array calculateTransmissibilities(const Problem& problem,
+ const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars,
+ const typename FVElementGeometry::SubControlVolumeFace& scvf,
+ const FluxVariablesCache& fluxVarsCache)
+ {
+ static_assert(ElementVolumeVariables::VolumeVariables::numFluidPhases() == 1);
+ const Scalar t = calculateTransmissibility(problem, element, fvGeometry, scvf, elemVolVars, fluxVarsCache, 0);
+ return std::array{t, -t};
+ }
+};
+
+
+} // end namespace
+
+#endif
diff --git a/dumux/flux/porenetwork/fickslaw.hh b/dumux/flux/porenetwork/fickslaw.hh
new file mode 100644
index 0000000000000000000000000000000000000000..038da6ef28b61570c8bb1b81e877cc3af54e35d2
--- /dev/null
+++ b/dumux/flux/porenetwork/fickslaw.hh
@@ -0,0 +1,114 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ * See the file COPYING for full copying permissions. *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 3 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see . *
+ *****************************************************************************/
+/*!
+ * \file
+ * \brief This file contains the data which is required to calculate
+ * diffusive mass fluxes due to molecular diffusion with Fick's law.
+ */
+#ifndef DUMUX_FLUX_PNM_FICKS_LAW_HH
+#define DUMUX_FLUX_PNM_FICKS_LAW_HH
+
+#include
+#include
+#include
+#include
+
+namespace Dumux
+{
+
+/*!
+ * \ingroup PoreNetworkFlux
+ * \brief Specialization of Fick's Law for the pore-network model.
+ */
+template
+class PNMFicksLaw
+{
+public:
+ using DiffusionCoefficientsContainer = FickianDiffusionCoefficients;
+
+ //return the reference system
+ static constexpr ReferenceSystemFormulation referenceSystemFormulation()
+ { return referenceSystem; }
+
+ template
+ static auto flux(const Problem& problem,
+ const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars,
+ const typename FVElementGeometry::SubControlVolumeFace& scvf,
+ const int phaseIdx,
+ const ElementFluxVariablesCache& elemFluxVarsCache)
+ {
+ Dune::FieldVector componentFlux(0.0);
+
+ // get inside and outside diffusion tensors and calculate the harmonic mean
+ const auto& insideVolVars = elemVolVars[scvf.insideScvIdx()];
+ const auto& outsideVolVars = elemVolVars[scvf.outsideScvIdx()];
+
+ const auto& fluxVarsCache = elemFluxVarsCache[scvf];
+
+ const Scalar density = 0.5 * (massOrMolarDensity(insideVolVars, referenceSystem, phaseIdx) + massOrMolarDensity(outsideVolVars, referenceSystem, phaseIdx));
+ const Scalar throatLength = fluxVarsCache.throatLength();
+ const Scalar phaseCrossSectionalArea = fluxVarsCache.throatCrossSectionalArea(phaseIdx);
+
+ for (int compIdx = 0; compIdx < numComponents; compIdx++)
+ {
+ if(compIdx == phaseIdx)
+ continue;
+
+ auto insideD = getDiffusionCoefficient_(phaseIdx, compIdx, insideVolVars);
+ auto outsideD = getDiffusionCoefficient_(phaseIdx, compIdx, outsideVolVars);
+
+ // scale by extrusion factor
+ insideD *= insideVolVars.extrusionFactor();
+ outsideD *= outsideVolVars.extrusionFactor();
+
+ // the resulting averaged diffusion coefficient
+ const auto D = harmonicMean(insideD, outsideD);
+
+ const Scalar insideMoleFraction = massOrMoleFraction(insideVolVars, referenceSystem, phaseIdx, compIdx);
+ const Scalar outsideMoleFraction = massOrMoleFraction(outsideVolVars, referenceSystem, phaseIdx, compIdx);
+
+ componentFlux[compIdx] = density * (insideMoleFraction - outsideMoleFraction) / throatLength * D * phaseCrossSectionalArea;
+ componentFlux[phaseIdx] -= componentFlux[compIdx];
+ }
+ return componentFlux;
+ }
+private:
+
+ template
+ static auto getDiffusionCoefficient_(const int phaseIdx, const int compIdx,
+ const VolumeVariables& volVars)
+ {
+ using FluidSystem = typename VolumeVariables::FluidSystem;
+
+ if constexpr (!FluidSystem::isTracerFluidSystem())
+ {
+ const auto mainCompIdx = FluidSystem::getMainComponent(phaseIdx);
+ return volVars.diffusionCoefficient(phaseIdx, mainCompIdx, compIdx);
+ }
+ else
+ return volVars.diffusionCoefficient(0, 0, compIdx);
+ }
+};
+} // end namespace
+
+#endif
diff --git a/dumux/flux/porenetwork/fourierslaw.hh b/dumux/flux/porenetwork/fourierslaw.hh
new file mode 100644
index 0000000000000000000000000000000000000000..38f1b438cfd69e28a194d27058576baacd8d75b6
--- /dev/null
+++ b/dumux/flux/porenetwork/fourierslaw.hh
@@ -0,0 +1,78 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ * See the file COPYING for full copying permissions. *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 3 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see . *
+ *****************************************************************************/
+/*!
+ * \file
+ * \brief This file contains the data which is required to calculate
+ * diffusive heat fluxes with Fourier's law.
+ */
+#ifndef DUMUX_FLUX_PNM_FOURIERS_LAW_HH
+#define DUMUX_FLUX_PNM_FOURIERS_LAW_HH
+
+#include
+
+namespace Dumux
+{
+
+ /*!
+ * \ingroup PoreNetworkFlux
+ * \brief Specialization of Fouriers's Law for the pore-network model.
+ */
+struct PNMFouriersLaw
+{
+
+ template
+ static auto flux(const Problem& problem,
+ const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars,
+ const typename FVElementGeometry::SubControlVolumeFace& scvf,
+ const ElementFluxVariablesCache& elemFluxVarsCache)
+ {
+ const auto& insideScv = fvGeometry.scv(scvf.insideScvIdx());
+ const auto& outsideScv = fvGeometry.scv(scvf.outsideScvIdx());
+ const auto& insideVolVars = elemVolVars[insideScv];
+ const auto& outsideVolVars = elemVolVars[outsideScv];
+ const auto& fluxVarsCache = elemFluxVarsCache[scvf];
+
+ static constexpr auto numPhases = ElementVolumeVariables::VolumeVariables::numFluidPhases();
+ using Scalar = typename ElementVolumeVariables::VolumeVariables::PrimaryVariables::value_type;
+ Scalar heatflux = 0;
+
+ for (int phaseIdx = 0; phaseIdx < numPhases; phaseIdx++)
+ {
+ auto insideThermalConducitivity = insideVolVars.fluidThermalConductivity(phaseIdx);
+ auto outsideThermalConducitivity = outsideVolVars.fluidThermalConductivity(phaseIdx);
+
+ auto thermalConductivity = Dumux::harmonicMean(insideThermalConducitivity, outsideThermalConducitivity);
+ auto area = fluxVarsCache.throatCrossSectionalArea(phaseIdx);
+
+ // calculate the temperature gradient
+ const Scalar deltaT = insideVolVars.temperature() - outsideVolVars.temperature();
+ const Scalar gradT = deltaT/fluxVarsCache.throatLength();
+
+ heatflux += thermalConductivity*gradT*area;
+ }
+ return heatflux;
+ }
+
+};
+} // end namespace
+
+#endif
diff --git a/dumux/io/CMakeLists.txt b/dumux/io/CMakeLists.txt
index a4080825f4862cfbe236a471bd0545bed008c305..43c86dc3849fabc92d791a00b6692d9edcff5551 100644
--- a/dumux/io/CMakeLists.txt
+++ b/dumux/io/CMakeLists.txt
@@ -12,6 +12,7 @@ name.hh
ploteffectivediffusivitymodel.hh
plotmateriallaw.hh
plotmateriallaw3p.hh
+plotpnmmateriallaw.hh
plotthermalconductivitymodel.hh
pointcloudvtkwriter.hh
rasterimagereader.hh
diff --git a/dumux/io/grid/CMakeLists.txt b/dumux/io/grid/CMakeLists.txt
index a5c54d8162c9c2e55c41047a7f3022980f739aed..7e21464ae142dda4fd79cab7b87ac6e7069484c9 100644
--- a/dumux/io/grid/CMakeLists.txt
+++ b/dumux/io/grid/CMakeLists.txt
@@ -1,3 +1,5 @@
+add_subdirectory(porenetwork)
+
install(FILES
cakegridmanager.hh
cpgridmanager.hh
diff --git a/dumux/io/grid/porenetwork/CMakeLists.txt b/dumux/io/grid/porenetwork/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..d5e9ae72a09fb52b35e1907e4a593f122ec1dac9
--- /dev/null
+++ b/dumux/io/grid/porenetwork/CMakeLists.txt
@@ -0,0 +1,7 @@
+install(FILES
+dgfwriter.hh
+griddata.hh
+gridmanager.hh
+parametersforgeneratedgrid.hh
+structuredlatticegridcreator.hh
+DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dumux/io/grid/porenetwork)
diff --git a/dumux/io/grid/porenetwork/dgfwriter.hh b/dumux/io/grid/porenetwork/dgfwriter.hh
new file mode 100644
index 0000000000000000000000000000000000000000..0475bd5d4dc4ba291a78dae17d53bab72dcf47e3
--- /dev/null
+++ b/dumux/io/grid/porenetwork/dgfwriter.hh
@@ -0,0 +1,95 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ * See the file COPYING for full copying permissions. *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 3 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see . *
+ *****************************************************************************/
+/*!
+ * \file
+ * \brief Write pore-network grids with attached data to dgf file
+ */
+#ifndef DUMUX_PORE_NETWORK_DGF_WRITER_HH
+#define DUMUX_PORE_NETWORK_DGF_WRITER_HH
+
+#include
+#include
+
+namespace Dumux {
+
+ /*!
+ * \ingroup PoreNetworkModel
+ * \brief Write pore-network grids with attached data to dgf file
+ */
+template
+inline static void writeDgf(const std::string& fileName, const GridView& gridView, const GridData& gridData)
+{
+ const auto someElement = *(elements(gridView).begin());
+ const auto someVertex = *(vertices(gridView).begin());
+ const auto numVertexParams = gridData.parameters(someVertex).size();
+ const auto numElementParams = gridData.parameters(someElement).size();
+
+ std::ofstream dgfFile;
+ dgfFile.open(fileName);
+ dgfFile << "DGF\nVertex % Coordinates, volumes and boundary flags of the pore bodies\nparameters " << numVertexParams << "\n";
+ dgfFile << "% Vertex parameters: ";
+ for (const auto& p : gridData.vertexParameterNames())
+ dgfFile << p << " ";
+ dgfFile << "\n% Element parameters: ";
+ for (const auto& p : gridData.elementParameterNames())
+ dgfFile << p << " ";
+ dgfFile << std::endl;
+
+ for (const auto& vertex : vertices(gridView))
+ {
+ dgfFile << vertex.geometry().center() << " ";
+
+ const auto& params = gridData.parameters(vertex);
+ for (int i = 0; i < params.size(); ++i)
+ {
+ dgfFile << params[i];
+
+ if (i < params.size() - 1)
+ dgfFile << " ";
+ }
+
+ dgfFile << std::endl;
+ }
+
+ dgfFile << "#\nSIMPLEX % Connections of the pore bodies (pore throats)\nparameters " << numElementParams << "\n";
+
+ for (const auto& element : elements(gridView))
+ {
+ dgfFile << gridView.indexSet().subIndex(element, 0, 1) << " ";
+ dgfFile << gridView.indexSet().subIndex(element, 1, 1) << " ";
+
+ const auto& params = gridData.parameters(element);
+ for (int i = 0; i < params.size(); ++i)
+ {
+ dgfFile << params[i];
+
+ if (i < params.size() - 1)
+ dgfFile << " ";
+ }
+
+ dgfFile << std::endl;
+ }
+
+ dgfFile << "#";
+ dgfFile.close();
+}
+
+}
+
+#endif
diff --git a/dumux/io/grid/porenetwork/griddata.hh b/dumux/io/grid/porenetwork/griddata.hh
new file mode 100644
index 0000000000000000000000000000000000000000..4dc28fdec4aef1f794584b17b0e36af6b8b64cf4
--- /dev/null
+++ b/dumux/io/grid/porenetwork/griddata.hh
@@ -0,0 +1,452 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ * See the file COPYING for full copying permissions. *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 3 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see . *
+ *****************************************************************************/
+/*!
+ * \file
+ * \ingroup InputOutput
+ * \brief Class for grid data attached to dgf or gmsh grid files
+ */
+#ifndef DUMUX_IO_PORENETWORKGRID_DATA_HH
+#define DUMUX_IO_PORENETWORKGRID_DATA_HH
+
+#include
+#include
+#include
+#include
+#include
+#include
+
+#include
+#include
+#include
+#include
+#include
+
+// FoamGrid specific includes
+#if HAVE_DUNE_FOAMGRID
+#include
+#include
+#endif
+
+#include
+#include
+#include
+
+#include "parametersforgeneratedgrid.hh"
+
+namespace Dumux {
+
+/*!
+ * \ingroup InputOutput
+ * \brief Class for grid data attached to dgf or gmsh grid files
+ */
+template
+class PoreNetworkGridData
+{
+ static constexpr int dim = Grid::dimension;
+ static constexpr int dimWorld = Grid::dimensionworld;
+ using Intersection = typename Grid::LeafIntersection;
+ using Element = typename Grid::template Codim<0>::Entity;
+ using Vertex = typename Grid::template Codim::Entity;
+ using GridView = typename Grid::LeafGridView;
+ using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
+ using SmallLocalIndex = typename IndexTraits::SmallLocalIndex;
+ using StringVector = std::vector;
+
+ using Scalar = double;
+ using PersistentParameterContainer = Dune::PersistentContainer>;
+ using ParametersForGeneratedGrid = Dumux::ParametersForGeneratedGrid;
+
+public:
+
+ //! constructor for dgf grid data
+ PoreNetworkGridData(Dune::GridPtr grid, const std::string& paramGroup)
+ : dgfGrid_(grid)
+ , isDgfData_(true)
+ , paramGroup_(paramGroup)
+ , numSubregions_(0)
+ {
+ setParameterIndices_();
+ }
+
+ //! constructor for non-dgf grid data
+ PoreNetworkGridData(std::shared_ptr grid, const std::string& paramGroup)
+ : factoryGrid_(grid)
+ , isDgfData_(false)
+ , paramGroup_(paramGroup)
+ {
+ numSubregions_ = getParamFromGroup(paramGroup_, "Grid.NumSubregions", 0);
+ setParameterIndices_();
+ parametersForGeneratedGrid_ = std::make_unique(gridView_(), paramGroup);
+ }
+
+ /*!
+ * \brief Call the parameters function of the DGF grid pointer if available for vertex data
+ * \note You can only pass vertices that exist on level 0!
+ */
+ const std::vector& parameters(const Vertex& vertex) const
+ {
+ if (isDgfData_ && !useCopiedDgfData_)
+ return dgfGrid_.parameters(vertex);
+ else
+ {
+ assert(!(*vertexParameters_)[vertex].empty() && "No parameters available.");
+ return (*vertexParameters_)[vertex];
+ }
+ }
+
+ /*!
+ * \brief Call the parameters function of the DGF grid pointer if available for element data
+ */
+ const std::vector& parameters(const Element& element) const
+ {
+ if (isDgfData_ && !useCopiedDgfData_)
+ {
+ if (element.hasFather())
+ {
+ auto level0Element = element;
+ while(level0Element.hasFather())
+ level0Element = level0Element.father();
+
+ return dgfGrid_.parameters(level0Element);
+ }
+ else
+ {
+ return dgfGrid_.parameters(element);
+ }
+ }
+ else
+ {
+ assert(!(*elementParameters_)[element].empty() && "No parameters available.");
+ return (*elementParameters_)[element];
+ }
+ }
+
+ /*!
+ * \brief Call the parameters function of the DGF grid pointer if available
+ */
+ template
+ const Dune::DGFBoundaryParameter::type& parameters(const Dune::Intersection& intersection) const
+ {
+ if (isDgfData_)
+ return dgfGrid_.parameters(intersection);
+ else
+ DUNE_THROW(Dune::InvalidStateException, "The parameters method is only available if the grid was constructed with a DGF file.");
+ }
+
+ /*!
+ * \brief Returns the coordination numbers for all pore bodies.
+ */
+ std::vector getCoordinationNumbers() const
+ {
+ std::vector coordinationNumbers(gridView_().size(dim), 0);
+
+ for (const auto &element : elements(gridView_()))
+ {
+ for (SmallLocalIndex vIdxLocal = 0; vIdxLocal < 2; ++vIdxLocal)
+ {
+ const auto vIdxGlobal = gridView_().indexSet().subIndex(element, vIdxLocal, dim);
+ coordinationNumbers[vIdxGlobal] += 1;
+ }
+ }
+
+ if (std::any_of(coordinationNumbers.begin(), coordinationNumbers.end(), [](auto i){ return i == 0; }))
+ DUNE_THROW(Dune::InvalidStateException, "One of the pores is not connected to another pore. SanitizeGrid will not help in this case. Check your grid file");
+
+ return coordinationNumbers;
+ }
+
+ /*!
+ * \brief Assign parameters for generically created grids
+ */
+ void assignParameters()
+ {
+ if (isDgfData_)
+ DUNE_THROW(Dune::InvalidStateException, "Assigning parameter not possible for dgf gids");
+
+ const auto numVertexParams = vertexParameterNames_.size();
+ const auto numElementParams = elementParameterNames_.size();
+ vertexParameters_ = makeParamContainer_(*factoryGrid_, numVertexParams, 1/*codim*/);
+ elementParameters_ = makeParamContainer_(*factoryGrid_, numElementParams, 0/*codim*/);
+
+ auto setParamHelper = [&](const auto& entity, const std::string& param, const Scalar value)
+ {
+ setParameter_(entity, param, value);
+ };
+
+ auto getParamHelper = [&](const auto& entity, const std::string& param)
+ {
+ return getParameter(entity, param);
+ };
+
+ parametersForGeneratedGrid_->assignParameters(setParamHelper, getParamHelper, numSubregions_);
+ }
+
+ void resizeParameterContainers()
+ {
+ // resize the parameters
+ vertexParameters_->resize();
+ elementParameters_->resize();
+ vertexParameters_->shrinkToFit();
+ elementParameters_->shrinkToFit();
+ }
+
+ void copyDgfData()
+ {
+ if (!isDgfData_)
+ DUNE_THROW(Dune::InvalidStateException, "copying dgf data only works when a dgf grid is actually used");
+
+ useCopiedDgfData_ = true;
+ const auto someVertex = *(vertices(gridView_()).begin());
+ const auto someElement = *(elements(gridView_()).begin());
+ const auto numVertexParams = dgfGrid_.parameters(someVertex).size();
+ const auto numElementParams = dgfGrid_.parameters(someElement).size();
+ vertexParameters_ = makeParamContainer_(*dgfGrid_, numVertexParams, 1);
+ elementParameters_ = makeParamContainer_(*dgfGrid_, numElementParams, 0);
+
+ for (const auto& element : elements(gridView_()))
+ {
+ for (int i = 0; i < numElementParams; ++i)
+ (*elementParameters_)[element][i] = dgfGrid_.parameters(element)[i];
+ }
+
+ for (const auto& vertex : vertices(gridView_()))
+ {
+ for (int i = 0; i < numVertexParams; ++i)
+ (*vertexParameters_)[vertex][i] = dgfGrid_.parameters(vertex)[i];
+ }
+ }
+
+ /*!
+ * \brief Return the index for a given parameter name
+ */
+ int parameterIndex(const std::string& paramName) const
+ {
+ // make sure the string is present in the map, throw a Dumux exception otherwise (and not a std one)
+ // the [] operator can't be used here due to const correctness
+ if (parameterIndex_.count(paramName))
+ return parameterIndex_.at(paramName);
+ else
+ {
+ std::string msg;
+ if (paramName.find("Throat") != std::string::npos)
+ msg = "Make sure to include it in the vector of parameter names ElementParameters = " + paramName + " ... ...";
+ else if (paramName.find("Pore") != std::string::npos)
+ msg = "Make sure to include it in the vector of parameter names VertexParameters = " + paramName + " ... ...";
+
+ DUNE_THROW(Dumux::ParameterException, paramName << " not set in the input file. \n" << msg);
+ }
+ }
+
+ /*!
+ * \brief Return the parameter group
+ */
+ const std::string& paramGroup() const
+ { return paramGroup_; }
+
+ /*!
+ * \brief Return if a given element parameter is provided by the grid
+ */
+ bool gridHasElementParameter(const std::string& param) const
+ {
+ return std::any_of(elementParameterNames_.begin(), elementParameterNames_.end(), [¶m]( const auto& i ){ return (i == param); });
+ }
+
+ /*!
+ * \brief Return if a given vertex parameter is provided by the grid
+ */
+ bool gridHasVertexParameter(const std::string& param) const
+ {
+ return std::any_of(vertexParameterNames_.begin(), vertexParameterNames_.end(), [¶m]( const auto& i ){ return (i == param); });
+ }
+
+ /*!
+ * \brief Returns the value of an element parameter
+ */
+ Scalar getParameter(const Element& element, const std::string& param) const
+ { return (*elementParameters_)[element][parameterIndex(param)]; }
+
+ /*!
+ * \brief Returns the value of an vertex parameter
+ */
+ Scalar getParameter(const Vertex& vertex, const std::string& param) const
+ { return (*vertexParameters_)[vertex][parameterIndex(param)]; }
+
+ /*!
+ * \brief Returns the pore label at a given position for a generic grid.
+ * This is needed by the grid creator in case not all parameters are initialized yet.
+ */
+ int poreLabelAtPosForGenericGrid(const GlobalPosition& pos) const
+ { return parametersForGeneratedGrid_->boundaryFaceMarkerAtPos(pos); }
+
+ /*!
+ * \brief Returns the names of the vertex parameters
+ */
+ const std::vector& vertexParameterNames() const
+ { return vertexParameterNames_; }
+
+ /*!
+ * \brief Returns the names of the element parameters
+ */
+ const std::vector& elementParameterNames() const
+ { return elementParameterNames_; }
+
+private:
+
+ void setParameter_(const Element& element, const std::string& param, const Scalar value)
+ { (*elementParameters_)[element][parameterIndex(param)] = value; }
+
+ void setParameter_(const Vertex& vertex, const std::string& param, const Scalar value)
+ { (*vertexParameters_)[vertex][parameterIndex(param)] = value; }
+
+ void setParameterIndices_()
+ {
+ if (!isDgfData_)
+ {
+ vertexParameterNames_ = StringVector{"PoreRadius", "PoreVolume", "PoreLabel"};
+ elementParameterNames_ = StringVector{"ThroatRadius", "ThroatLength", "ThroatLabel"};
+ if (numSubregions_ > 0)
+ {
+ vertexParameterNames_.push_back("PoreRegionId");
+ elementParameterNames_.push_back("ThroatRegionId");
+ }
+ }
+ else // DGF grid data
+ {
+ // treat vertex parameter names
+ if (auto inputFileVertexParameterNames = getParamFromGroup(paramGroup_, "Grid.VertexParameters", StringVector{}); !inputFileVertexParameterNames.empty())
+ {
+ // TODO Remove at some point
+ std::cout << "\n***\nWARNING: Setting Grid.VertexParameters in the input file is deprecated. Set '% Vertex parameters: Param1 Param2 ...' in the dgf file instead\n***\n" << std::endl;
+ vertexParameterNames_ = std::move(inputFileVertexParameterNames);
+ }
+ else
+ {
+ if (auto dgfFileVertexParameterNames = dgfFileParameterNames_("Vertex"); dgfFileVertexParameterNames.empty())
+ DUNE_THROW(Dune::InvalidStateException, "No vertex parameter names specified in dgf file. Set '% Vertex parameters: Param1 Param2 ...'");
+ else
+ vertexParameterNames_ = std::move(dgfFileVertexParameterNames);
+ }
+
+ // treat element parameter names
+ if (auto inputFileElementParameterNames = getParamFromGroup(paramGroup_, "Grid.ElementParameters", StringVector{}); !inputFileElementParameterNames.empty())
+ {
+ // TODO Remove at some point
+ std::cout << "\n***\nWARNING: Setting Grid.ElementParameters in the input file is deprecated. Set '% Element parameters: Param1 Param2 ...' in the dgf file instead\n***\n" << std::endl;
+ elementParameterNames_ = std::move(inputFileElementParameterNames);
+ }
+ else
+ {
+ if (auto dgfFileElementParameterNames = dgfFileParameterNames_("Element"); dgfFileElementParameterNames.empty())
+ DUNE_THROW(Dune::InvalidStateException, "No element parameter names specified in dgf file. Set '% Element parameters: Param1 Param2 ...'");
+ else
+ elementParameterNames_ = std::move(dgfFileElementParameterNames);
+ }
+
+ // make sure that the number of specified parameters matches with the dgf file
+ if (const auto& someElement = *(elements(gridView_()).begin()); elementParameterNames_.size() != dgfGrid_.nofParameters(someElement))
+ DUNE_THROW(Dune::InvalidStateException, "Number of user-specified element parameters (" << elementParameterNames_.size()
+ << ") does not match number of element paramters in dgf file (" << dgfGrid_.nofParameters(someElement) << ")");
+
+ if (const auto& someVertex = *(vertices(gridView_()).begin()); vertexParameterNames_.size() != dgfGrid_.nofParameters(someVertex))
+ DUNE_THROW(Dune::InvalidStateException, "Number of user-specified vertex parameters (" << vertexParameterNames_.size()
+ << ") does not match number of vertex paramters in dgf file (" << dgfGrid_.nofParameters(someVertex) << ")");
+ }
+
+ for (int i = 0; i < vertexParameterNames_.size(); ++i)
+ {
+ std::cout << vertexParameterNames_[i] << " is vertex parameter " << i << std::endl;
+ parameterIndex_[vertexParameterNames_[i]] = i;
+ }
+
+ for (int i = 0; i < elementParameterNames_.size(); ++i)
+ {
+ std::cout << elementParameterNames_[i] << " is element parameter " << i << std::endl;
+ parameterIndex_[elementParameterNames_[i]] = i;
+ }
+ }
+
+ /*!
+ * \brief Initializes and returns a container for vertex (codim dim) or element (codim 0) data
+ *
+ * \param grid The grid
+ * \param numParams The number of paramters
+ * \param codim The codimension
+ */
+ auto makeParamContainer_(const Grid& grid, int numParams, int codim) const
+ {
+ auto parameters = std::make_unique(grid, codim);
+ (*parameters).resize();
+ for (auto&& v : (*parameters))
+ v.resize(numParams);
+ return parameters;
+ }
+
+ StringVector dgfFileParameterNames_(const std::string& entity) const
+ {
+ std::ifstream gridFile(getParamFromGroup(paramGroup_, "Grid.File"));
+ std::string line;
+ while (getline(gridFile, line))
+ {
+ if (line.find(entity + " parameters:", 0) != std::string::npos)
+ {
+ std::string args = line.substr(line.find(":")+1, std::string::npos);
+ StringVector paramNames;
+ std::istringstream iss(args);
+ std::string item;
+ while (std::getline(iss, item, ' '))
+ if (!item.empty())
+ *std::back_inserter(paramNames)++ = item;
+
+ return paramNames;
+ }
+ }
+ return StringVector();
+ }
+
+ /*!
+ * \brief Return the gridView this grid geometry object lives on
+ */
+ const GridView gridView_() const
+ { return isDgfData_ ? dgfGrid_->leafGridView() : factoryGrid_->leafGridView(); }
+
+ // dgf grid data
+ Dune::GridPtr dgfGrid_;
+
+ std::shared_ptr factoryGrid_;
+
+ bool isDgfData_ = false;
+ bool useCopiedDgfData_ = false;
+ std::string paramGroup_;
+
+ std::unique_ptr parametersForGeneratedGrid_;
+
+ std::vector vertexParameterNames_;
+ std::vector elementParameterNames_;
+
+ std::unique_ptr vertexParameters_;
+ std::unique_ptr elementParameters_;
+
+ std::size_t numSubregions_;
+
+ std::unordered_map parameterIndex_;
+};
+
+} // namespace Dumux
+
+#endif
diff --git a/dumux/io/grid/porenetwork/gridmanager.hh b/dumux/io/grid/porenetwork/gridmanager.hh
new file mode 100644
index 0000000000000000000000000000000000000000..880095b6ddba005cf5fb7e097acadfdc6f55b74a
--- /dev/null
+++ b/dumux/io/grid/porenetwork/gridmanager.hh
@@ -0,0 +1,361 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ * See the file COPYING for full copying permissions. *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 3 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see . *
+ *****************************************************************************/
+/*!
+ * \file
+ * \brief Provides a grid manager for pore-network models
+ */
+#ifndef DUMUX_PORE_NETWORK_GRID_MANAGER_HH
+#define DUMUX_PORE_NETWORK_GRID_MANAGER_HH
+
+#include
+
+#include
+#include
+#include
+
+// FoamGrid specific includes
+#if HAVE_DUNE_FOAMGRID
+#include
+#include
+#else
+static_assert(false, "dune-foamgrid required!");
+#endif
+
+#include
+
+#include "griddata.hh"
+#include "structuredlatticegridcreator.hh"
+
+namespace Dumux {
+
+/*!
+ * \ingroup PoreNetworkModel
+ * \brief A grid manager for pore-network models
+ */
+template>>
+class PoreNetworkGridManager
+{
+ using GridType = Dune::FoamGrid<1, dimWorld>;
+ using Element = typename GridType::template Codim<0>::Entity;
+ using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
+
+public:
+ using Grid = GridType;
+ using GridData = GData;
+
+ void init(const std::string& paramGroup = "")
+ {
+ Dune::Timer timer;
+ paramGroup_ = paramGroup;
+
+ // First try to create it from a DGF file in GridParameterGroup.File
+ if (hasParamInGroup(paramGroup, "Grid.File"))
+ {
+ makeGridFromDgfFile(getParamFromGroup(paramGroup, "Grid.File"));
+ loadBalance();
+ }
+ else // no grid file found
+ {
+ try
+ {
+ createOneDGrid_();
+ }
+ catch(Dune::RangeError& e)
+ {
+ makeGridFromStructuredLattice();
+ }
+
+ loadBalance();
+ }
+
+ timer.stop();
+ const auto gridType = enableDgfGridPointer_ ? "grid from dgf file" : "generic grid from structured lattice";
+ std::cout << "Creating " << gridType << " with " << grid().leafGridView().size(0) << " elements and "
+ << grid().leafGridView().size(1) << " vertices took " << timer.elapsed() << " seconds." << std::endl;
+ }
+
+ /*!
+ * \brief Make a grid from a DGF file.
+ */
+ void makeGridFromDgfFile(const std::string& fileName)
+ {
+ // We found a file in the input file...does it have a supported extension?
+ const std::string extension = getFileExtension(fileName);
+ if (extension != "dgf")
+ DUNE_THROW(Dune::IOError, "Grid type " << Dune::className() << " only supports DGF (*.dgf) but the specified filename has extension: *."<< extension);
+
+ enableDgfGridPointer_ = true;
+ dgfGridPtr() = Dune::GridPtr(fileName.c_str(), Dune::MPIHelper::getCommunicator());
+ gridData_ = std::make_shared(dgfGridPtr_, paramGroup_);
+
+ if (getParamFromGroup(paramGroup_, "Grid.Sanitize", false))
+ sanitizeGrid();
+ }
+
+ /*!
+ * \brief Make a grid based on a structured lattice which allows to randomly delete elements based on Raoof et al. 2009
+ */
+ void makeGridFromStructuredLattice()
+ {
+ StructuredLatticeGridCreator creator;
+ creator.init(paramGroup_);
+ gridPtr() = creator.gridPtr();
+
+ gridData_ = std::make_shared(gridPtr_, paramGroup_);
+
+ if (getParamFromGroup(paramGroup_, "Grid.Sanitize", true))
+ sanitizeGrid();
+ else
+ std::cout << "\nWARNING: Set Grid.Sanitize = true in order to remove insular patches of elements not connected to the boundary." << std::endl;
+
+ gridData_->assignParameters();
+ }
+
+ /*!
+ * \brief Returns the filename extension of a given filename
+ */
+ std::string getFileExtension(const std::string& fileName) const
+ {
+ 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 << "')!");
+ }
+ return "";
+ }
+
+ /*!
+ * \brief Returns a reference to the grid.
+ */
+ Grid& grid()
+ { return enableDgfGridPointer_ ? *dgfGridPtr() : *gridPtr(); }
+
+ /*!
+ * \brief Returns the grid data.
+ */
+ std::shared_ptr getGridData() const
+ {
+ if (!gridData_)
+ DUNE_THROW(Dune::IOError, "No grid data available");
+
+ return gridData_;
+ }
+
+ /*!
+ * \brief Call loadBalance() function of the grid.
+ */
+ void loadBalance()
+ {
+ if (Dune::MPIHelper::getCollectiveCommunication().size() > 1)
+ {
+ // if we may have dgf parameters use load balancing of the dgf pointer
+ if (enableDgfGridPointer_)
+ {
+ dgfGridPtr().loadBalance();
+ // update the grid data
+ gridData_ = std::make_shared(dgfGridPtr(), paramGroup_);
+ }
+
+ else
+ gridPtr()->loadBalance();
+ }
+ }
+
+ /*!
+ * \brief post-processing to remove insular groups of elements that are not connected to a Dirichlet boundary
+ */
+ void sanitizeGrid()
+ {
+ // evaluate the coordination numbers to check if all pores are connected to at least one throat
+ gridData_->getCoordinationNumbers();
+
+ // for dgf grid, copy the data to peristent containers
+ if (enableDgfGridPointer_)
+ gridData_->copyDgfData();
+
+ // pruning -- remove elements not connected to a Dirichlet boundary (marker == 1)
+ const auto pruningSeedIndices = getParamFromGroup>(paramGroup_, "Grid.PruningSeedIndices", std::vector{1});
+ const auto gridView = grid().leafGridView();
+ std::vector elementIsConnected(gridView.size(0), false);
+
+ auto boundaryIdx = [&](const auto& vertex)
+ {
+ if (enableDgfGridPointer_)
+ return static_cast(dgfGridPtr_.parameters(vertex)[gridData_->parameterIndex("PoreLabel")]);
+ else
+ return static_cast(gridData_->poreLabelAtPosForGenericGrid(vertex.geometry().center()));
+ };
+
+ for (const auto& element : elements(gridView))
+ {
+ const auto eIdx = gridView.indexSet().index(element);
+ if (elementIsConnected[eIdx])
+ continue;
+
+ // try to find a seed from which to start the search process
+ bool isSeed = false;
+ bool hasNeighbor = false;
+ for (const auto& intersection : intersections(gridView, element))
+ {
+ auto vertex = element.template subEntity<1>(intersection.indexInInside());
+ // seed found
+ if (std::any_of(pruningSeedIndices.begin(), pruningSeedIndices.end(),
+ [&]( const int i ){ return i == boundaryIdx(vertex); }))
+ {
+ isSeed = true;
+ // break;
+ }
+
+ if (intersection.neighbor())
+ hasNeighbor = true;
+ }
+
+ if (!hasNeighbor)
+ continue;
+
+ if (isSeed)
+ {
+ elementIsConnected[eIdx] = true;
+
+ // use iteration instead of recursion here because the recursion can get too deep
+ std::stack elementStack;
+ elementStack.push(element);
+ while (!elementStack.empty())
+ {
+ auto e = elementStack.top();
+ elementStack.pop();
+ for (const auto& intersection : intersections(gridView, e))
+ {
+ if (!intersection.boundary())
+ {
+ auto outside = intersection.outside();
+ auto nIdx = gridView.indexSet().index(outside);
+ if (!elementIsConnected[nIdx])
+ {
+ elementIsConnected[nIdx] = true;
+ elementStack.push(outside);
+ }
+ }
+ }
+ }
+ }
+ }
+
+ if (std::none_of(elementIsConnected.begin(), elementIsConnected.end(), [](const bool i){return i;}))
+ DUNE_THROW(Dune::InvalidStateException, "sanitize() deleted all elements! Check your boundary face indices. "
+ << "If the problem persisits, add at least one of your boundary face indices to PruningSeedIndices");
+
+ // remove the elements in the grid
+ std::size_t numDeletedElements = 0;
+ grid().preGrow();
+ for (const auto& element : elements(gridView))
+ {
+ const auto eIdx = gridView.indexSet().index(element);
+ if (!elementIsConnected[eIdx])
+ {
+ grid().removeElement(element);
+ ++numDeletedElements;
+ }
+ }
+ // triggers the grid growth process
+ grid().grow();
+ grid().postGrow();
+
+ // resize the parameters for dgf grids
+ if (enableDgfGridPointer_)
+ gridData_->resizeParameterContainers();
+
+ if (numDeletedElements > 0)
+ std::cout << "\nDeleted " << numDeletedElements << " elements not connected to a boundary." << std::endl;
+ }
+
+protected:
+
+ /*!
+ * \brief Returns a reference to the grid pointer (std::shared_ptr)
+ */
+ std::shared_ptr& gridPtr()
+ {
+ if(!enableDgfGridPointer_)
+ return gridPtr_;
+ else
+ DUNE_THROW(Dune::InvalidStateException, "You are using DGF. To get the grid pointer use method dgfGridPtr()!");
+ }
+
+ /*!
+ * \brief Returns a reference to the DGF grid pointer (Dune::GridPtr).
+ */
+ Dune::GridPtr& dgfGridPtr()
+ {
+ if(enableDgfGridPointer_)
+ return dgfGridPtr_;
+ else
+ DUNE_THROW(Dune::InvalidStateException, "The DGF grid pointer is only available if the grid was constructed with a DGF file!");
+ }
+
+ /*!
+ * \brief A state variable if the DGF Dune::GridPtr has been enabled.
+ * It is always enabled if a DGF grid file was used to create the grid.
+ */
+ bool enableDgfGridPointer_ = false;
+
+ std::shared_ptr gridPtr_;
+ Dune::GridPtr dgfGridPtr_;
+
+ std::shared_ptr gridData_;
+
+ std::string paramGroup_;
+
+private:
+ void createOneDGrid_()
+ {
+ const auto lowerLeft = getParamFromGroup