From 805340720770401ecb3fc1afef1628c2a33feec5 Mon Sep 17 00:00:00 2001 From: Kilian Weishaupt Date: Mon, 13 Jul 2020 09:58:03 +0200 Subject: [PATCH 01/25] [porenetworkmodel] Add files for 1p, 1pnc and 2p PNM --- dumux/CMakeLists.txt | 1 + dumux/discretization/CMakeLists.txt | 1 + .../discretization/porenetwork/CMakeLists.txt | 6 + .../porenetwork/fvelementgeometry.hh | 326 ++++++++ .../porenetwork/gridgeometry.hh | 760 ++++++++++++++++++ .../porenetwork/subcontrolvolume.hh | 143 ++++ .../porenetwork/subcontrolvolumeface.hh | 148 ++++ dumux/flux/CMakeLists.txt | 1 + dumux/flux/porenetwork/CMakeLists.txt | 5 + dumux/flux/porenetwork/advection.hh | 145 ++++ dumux/flux/porenetwork/fickslaw.hh | 114 +++ dumux/flux/porenetwork/fourierslaw.hh | 78 ++ dumux/io/CMakeLists.txt | 1 + dumux/io/grid/CMakeLists.txt | 2 + dumux/io/grid/porenetwork/CMakeLists.txt | 7 + dumux/io/grid/porenetwork/dgfwriter.hh | 95 +++ dumux/io/grid/porenetwork/griddata.hh | 452 +++++++++++ dumux/io/grid/porenetwork/gridmanager.hh | 361 +++++++++ .../porenetwork/parametersforgeneratedgrid.hh | 621 ++++++++++++++ .../structuredlatticegridcreator.hh | 607 ++++++++++++++ dumux/io/plotpnmmateriallaw.hh | 329 ++++++++ .../fluidmatrixinteractions/CMakeLists.txt | 1 + .../porenetwork/CMakeLists.txt | 10 + .../porenetwork/emptycache.hh | 36 + .../porenetwork/porenetworklocalrules.hh | 188 +++++ .../porenetworklocalrulesparams.hh | 79 ++ .../regularizedporenetworklocalrules.hh | 240 ++++++ .../regularizedporenetworklocalrulesparams.hh | 104 +++ .../thresholdcapillarypressures.hh | 78 ++ .../porenetwork/transmissibility1p.hh | 238 ++++++ .../porenetwork/transmissibility2p.hh | 212 +++++ dumux/material/spatialparams/CMakeLists.txt | 2 + .../spatialparams/porenetwork/CMakeLists.txt | 5 + .../porenetwork/porenetwork1p.hh | 93 +++ .../porenetwork/porenetwork2p.hh | 179 +++++ .../porenetwork/porenetworkbase.hh | 290 +++++++ dumux/porenetworkflow/1p/CMakeLists.txt | 6 + .../porenetworkflow/1p/fluxvariablescache.hh | 106 +++ dumux/porenetworkflow/1p/iofields.hh | 61 ++ dumux/porenetworkflow/1p/model.hh | 197 +++++ dumux/porenetworkflow/1p/volumevariables.hh | 82 ++ dumux/porenetworkflow/1pnc/CMakeLists.txt | 5 + dumux/porenetworkflow/1pnc/iofields.hh | 63 ++ dumux/porenetworkflow/1pnc/model.hh | 250 ++++++ dumux/porenetworkflow/1pnc/volumevariables.hh | 86 ++ dumux/porenetworkflow/2p/CMakeLists.txt | 14 + .../2p/elementfluxvariablescache.hh | 125 +++ .../porenetworkflow/2p/fluxvariablescache.hh | 198 +++++ .../2p/gridfluxvariablescache.hh | 186 +++++ dumux/porenetworkflow/2p/invasionstate.hh | 321 ++++++++ dumux/porenetworkflow/2p/iofields.hh | 86 ++ dumux/porenetworkflow/2p/model.hh | 267 ++++++ .../2p/multidomainnewtonsolver.hh | 329 ++++++++ .../2p/newtonconsistencychecks.hh | 166 ++++ dumux/porenetworkflow/2p/newtonsolver.hh | 107 +++ .../porenetworkflow/2p/static/CMakeLists.txt | 3 + .../2p/static/staticdrainge.hh | 159 ++++ dumux/porenetworkflow/2p/volumevariables.hh | 101 +++ dumux/porenetworkflow/CMakeLists.txt | 8 + dumux/porenetworkflow/common/CMakeLists.txt | 11 + dumux/porenetworkflow/common/boundaryflux.hh | 209 +++++ dumux/porenetworkflow/common/geometry.hh | 66 ++ dumux/porenetworkflow/common/iofields.hh | 60 ++ dumux/porenetworkflow/common/labels.hh | 49 ++ .../common/pnmvtkoutputmodule.hh | 181 +++++ .../porenetworkflow/common/poreproperties.hh | 96 +++ .../common/throatproperties.hh | 257 ++++++ dumux/porenetworkflow/common/utilities.hh | 129 +++ .../porenetworkflow/common/velocityoutput.hh | 143 ++++ dumux/porenetworkflow/properties.hh | 109 +++ 70 files changed, 10194 insertions(+) create mode 100644 dumux/discretization/porenetwork/CMakeLists.txt create mode 100644 dumux/discretization/porenetwork/fvelementgeometry.hh create mode 100644 dumux/discretization/porenetwork/gridgeometry.hh create mode 100644 dumux/discretization/porenetwork/subcontrolvolume.hh create mode 100644 dumux/discretization/porenetwork/subcontrolvolumeface.hh create mode 100644 dumux/flux/porenetwork/CMakeLists.txt create mode 100644 dumux/flux/porenetwork/advection.hh create mode 100644 dumux/flux/porenetwork/fickslaw.hh create mode 100644 dumux/flux/porenetwork/fourierslaw.hh create mode 100644 dumux/io/grid/porenetwork/CMakeLists.txt create mode 100644 dumux/io/grid/porenetwork/dgfwriter.hh create mode 100644 dumux/io/grid/porenetwork/griddata.hh create mode 100644 dumux/io/grid/porenetwork/gridmanager.hh create mode 100644 dumux/io/grid/porenetwork/parametersforgeneratedgrid.hh create mode 100644 dumux/io/grid/porenetwork/structuredlatticegridcreator.hh create mode 100644 dumux/io/plotpnmmateriallaw.hh create mode 100644 dumux/material/fluidmatrixinteractions/porenetwork/CMakeLists.txt create mode 100644 dumux/material/fluidmatrixinteractions/porenetwork/emptycache.hh create mode 100644 dumux/material/fluidmatrixinteractions/porenetwork/porenetworklocalrules.hh create mode 100644 dumux/material/fluidmatrixinteractions/porenetwork/porenetworklocalrulesparams.hh create mode 100644 dumux/material/fluidmatrixinteractions/porenetwork/regularizedporenetworklocalrules.hh create mode 100644 dumux/material/fluidmatrixinteractions/porenetwork/regularizedporenetworklocalrulesparams.hh create mode 100644 dumux/material/fluidmatrixinteractions/porenetwork/thresholdcapillarypressures.hh create mode 100644 dumux/material/fluidmatrixinteractions/porenetwork/transmissibility1p.hh create mode 100644 dumux/material/fluidmatrixinteractions/porenetwork/transmissibility2p.hh create mode 100644 dumux/material/spatialparams/porenetwork/CMakeLists.txt create mode 100644 dumux/material/spatialparams/porenetwork/porenetwork1p.hh create mode 100644 dumux/material/spatialparams/porenetwork/porenetwork2p.hh create mode 100644 dumux/material/spatialparams/porenetwork/porenetworkbase.hh create mode 100644 dumux/porenetworkflow/1p/CMakeLists.txt create mode 100644 dumux/porenetworkflow/1p/fluxvariablescache.hh create mode 100644 dumux/porenetworkflow/1p/iofields.hh create mode 100644 dumux/porenetworkflow/1p/model.hh create mode 100644 dumux/porenetworkflow/1p/volumevariables.hh create mode 100644 dumux/porenetworkflow/1pnc/CMakeLists.txt create mode 100644 dumux/porenetworkflow/1pnc/iofields.hh create mode 100644 dumux/porenetworkflow/1pnc/model.hh create mode 100644 dumux/porenetworkflow/1pnc/volumevariables.hh create mode 100644 dumux/porenetworkflow/2p/CMakeLists.txt create mode 100644 dumux/porenetworkflow/2p/elementfluxvariablescache.hh create mode 100644 dumux/porenetworkflow/2p/fluxvariablescache.hh create mode 100644 dumux/porenetworkflow/2p/gridfluxvariablescache.hh create mode 100644 dumux/porenetworkflow/2p/invasionstate.hh create mode 100644 dumux/porenetworkflow/2p/iofields.hh create mode 100644 dumux/porenetworkflow/2p/model.hh create mode 100644 dumux/porenetworkflow/2p/multidomainnewtonsolver.hh create mode 100644 dumux/porenetworkflow/2p/newtonconsistencychecks.hh create mode 100644 dumux/porenetworkflow/2p/newtonsolver.hh create mode 100644 dumux/porenetworkflow/2p/static/CMakeLists.txt create mode 100644 dumux/porenetworkflow/2p/static/staticdrainge.hh create mode 100644 dumux/porenetworkflow/2p/volumevariables.hh create mode 100644 dumux/porenetworkflow/CMakeLists.txt create mode 100644 dumux/porenetworkflow/common/CMakeLists.txt create mode 100644 dumux/porenetworkflow/common/boundaryflux.hh create mode 100644 dumux/porenetworkflow/common/geometry.hh create mode 100644 dumux/porenetworkflow/common/iofields.hh create mode 100644 dumux/porenetworkflow/common/labels.hh create mode 100644 dumux/porenetworkflow/common/pnmvtkoutputmodule.hh create mode 100644 dumux/porenetworkflow/common/poreproperties.hh create mode 100644 dumux/porenetworkflow/common/throatproperties.hh create mode 100644 dumux/porenetworkflow/common/utilities.hh create mode 100644 dumux/porenetworkflow/common/velocityoutput.hh create mode 100644 dumux/porenetworkflow/properties.hh diff --git a/dumux/CMakeLists.txt b/dumux/CMakeLists.txt index 819d019ad5..1aae30ec2a 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 a43e90fb37..6fbdf151ab 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 0000000000..27a356b4c2 --- /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 0000000000..8f131b7f01 --- /dev/null +++ b/dumux/discretization/porenetwork/fvelementgeometry.hh @@ -0,0 +1,326 @@ +// -*- 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 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_, + elementGeometry.corner(scvLocalIdx), + 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 0000000000..5b57388e47 --- /dev/null +++ b/dumux/discretization/porenetwork/gridgeometry.hh @@ -0,0 +1,760 @@ +// -*- 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 + +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"); + static const auto throatLabelIdx = gridData.parameterIndex("ThroatLabel"); + throatRadius_[eIdx] = params[throatRadiusIdx]; + throatLength_[eIdx] = params[throatLengthIdx]; + throatLabel_[eIdx] = params[throatLabelIdx]; + + 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