diff --git a/dumux/discretization/staggered/fvelementgeometry.hh b/dumux/discretization/staggered/fvelementgeometry.hh
index 2a2213acac579eb4a7c2cb1695aa76a4d8f22437..1b98bcc29d61abbeae9a2e11d1dc1e66219c2c44 100644
--- a/dumux/discretization/staggered/fvelementgeometry.hh
+++ b/dumux/discretization/staggered/fvelementgeometry.hh
@@ -28,18 +28,30 @@
namespace Dumux {
+/*!
+ * \ingroup StaggeredDiscretization
+ * \brief Stencil-local finite volume geometry (scvs and scvfs) for staggered models
+ * This builds up the sub control volumes and sub control volume faces
+ * for each element in the local scope we are restricting to, e.g. stencil or element.
+ * \tparam GG the finite volume grid geometry type
+ * \tparam enableFVGridGeometryCache if the grid geometry is cached or not
+ * \note This class is specialized for versions with and without caching the fv geometries on the grid view
+ */
+template<class GG, bool enableFVGridGeometryCache>
+class StaggeredFVElementGeometry;
+
/*!
* \ingroup StaggeredDiscretization
* \brief Base class for the finite volume geometry vector for staggered models
* This locally builds up the sub control volumes and sub control volume faces
* for each element.
+ * Specialization for grid caching enabled
* \tparam GG the finite volume grid geometry type
- * \tparam enableFVGridGeometryCache if the grid geometry is cached or not
*/
-template<class GG, bool enableFVGridGeometryCache>
-class StaggeredFVElementGeometry : public CCTpfaFVElementGeometry<GG, enableFVGridGeometryCache>
+template<class GG>
+class StaggeredFVElementGeometry<GG, true> : public CCTpfaFVElementGeometry<GG, true>
{
- using ParentType = CCTpfaFVElementGeometry<GG, enableFVGridGeometryCache>;
+ using ParentType = CCTpfaFVElementGeometry<GG, true>;
using IndexType = typename GG::GridView::IndexSet::IndexType;
public:
//! export type of subcontrol volume face
@@ -60,6 +72,254 @@ public:
}
};
+/*!
+ * \ingroup StaggeredDiscretization
+ * \brief Base class for the finite volume geometry vector for staggered models
+ * This locally builds up the sub control volumes and sub control volume faces
+ * for each element.
+ * Specialization for grid caching enabled
+ * \tparam GG the finite volume grid geometry type
+ */
+template<class GG>
+class StaggeredFVElementGeometry<GG, false>
+{
+ using ThisType = StaggeredFVElementGeometry<GG, false>;
+ using GridView = typename GG::GridView;
+ using IndexType = typename GridView::IndexSet::IndexType;
+ using Element = typename GridView::template Codim<0>::Entity;
+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 FVGridGeometry = GG;
+
+ //! Constructor getting a auxiliary cell center of face specific FvGridGeometry type.
+ //! Needed for the multi-domain framework.
+ template<class CellCenterOrFaceFVGridGeometry>
+ StaggeredFVElementGeometry(const CellCenterOrFaceFVGridGeometry& fvGridGeometry)
+ : fvGridGeometryPtr_(&fvGridGeometry.actualfvGridGeometry()) {}
+
+ //! Constructor
+ StaggeredFVElementGeometry(const FVGridGeometry& fvGridGeometry)
+ : fvGridGeometryPtr_(&fvGridGeometry) {}
+
+ //! Get a sub control volume face with an element index and a local scvf index
+ const SubControlVolumeFace& scvf(IndexType eIdx, IndexType localScvfIdx) const
+ {
+ return scvf(this->fvGridGeometry().localToGlobalScvfIndex(eIdx, localScvfIdx));
+ }
+
+ //! Get an elment sub control volume with a global scv index
+ //! We separate element and neighbor scvs to speed up mapping
+ const SubControlVolume& scv(IndexType scvIdx) const
+ {
+ if (scvIdx == scvIndices_[0])
+ return scvs_[0];
+ else
+ return neighborScvs_[findLocalIndex_(scvIdx, neighborScvIndices_)];
+ }
+
+ //! Get an element sub control volume face with a global scvf index
+ //! We separate element and neighbor scvfs to speed up mapping
+ const SubControlVolumeFace& scvf(IndexType scvfIdx) const
+ {
+ auto it = std::find(scvfIndices_.begin(), scvfIndices_.end(), scvfIdx);
+ if (it != scvfIndices_.end())
+ return scvfs_[std::distance(scvfIndices_.begin(), it)];
+ else
+ return neighborScvfs_[findLocalIndex_(scvfIdx, neighborScvfIndices_)];
+ }
+
+ //! iterator range for sub control volumes. Iterates over
+ //! all scvs of the bound element (not including neighbor scvs)
+ //! 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 Dune::IteratorRange<typename std::array<SubControlVolume, 1>::const_iterator>
+ scvs(const ThisType& g)
+ {
+ using IteratorType = typename std::array<SubControlVolume, 1>::const_iterator;
+ return Dune::IteratorRange<IteratorType>(g.scvs_.begin(), g.scvs_.end());
+ }
+
+ //! iterator range for sub control volumes faces. Iterates over
+ //! all scvfs of the bound element (not including neighbor scvfs)
+ //! 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 Dune::IteratorRange<typename std::vector<SubControlVolumeFace>::const_iterator>
+ scvfs(const ThisType& g)
+ {
+ using IteratorType = typename std::vector<SubControlVolumeFace>::const_iterator;
+ return Dune::IteratorRange<IteratorType>(g.scvfs_.begin(), g.scvfs_.end());
+ }
+
+ //! number of sub control volumes in this fv element geometry
+ std::size_t numScv() const
+ { return scvs_.size(); }
+
+ //! number of sub control volumes in this fv element geometry
+ std::size_t numScvf() const
+ { return scvfs_.size(); }
+
+ //! Binding of an element preparing the geometries of the whole stencil
+ //! called by the local jacobian to prepare element assembly
+ void bind(const Element& element)
+ {
+ bindElement(element);
+
+ neighborScvs_.reserve(element.subEntities(1));
+ neighborScvfIndices_.reserve(element.subEntities(1));
+ neighborScvfs_.reserve(element.subEntities(1));
+
+ std::vector<IndexType> handledNeighbors;
+ handledNeighbors.reserve(element.subEntities(1));
+ for (const auto& intersection : intersections(fvGridGeometry().gridView(), element))
+ {
+ if (intersection.neighbor())
+ {
+ const auto outside = intersection.outside();
+ const auto outsideIdx = fvGridGeometry().elementMapper().index(outside);
+
+ // make outside geometries only if not done yet (could happen on non-conforming grids)
+ if ( std::find(handledNeighbors.begin(), handledNeighbors.end(), outsideIdx) == handledNeighbors.end() )
+ {
+ makeNeighborGeometries_(outside, outsideIdx);
+ handledNeighbors.push_back(outsideIdx);
+ }
+ }
+ }
+ }
+
+ //! Binding of an element preparing the geometries only inside the element
+ void bindElement(const Element& element)
+ {
+ clear_();
+ elementPtr_ = &element;
+ scvfs_.reserve(element.subEntities(1));
+ scvfIndices_.reserve(element.subEntities(1));
+ makeElementGeometries_(element);
+ }
+
+ //! The global finite volume geometry we are a restriction of
+ const FVGridGeometry& fvGridGeometry() const
+ { return *fvGridGeometryPtr_; }
+
+private:
+
+ //! create scvs and scvfs of the bound element
+ void makeElementGeometries_(const Element& element)
+ {
+ const auto eIdx = fvGridGeometry().elementMapper().index(element);
+ scvs_[0] = SubControlVolume(element.geometry(), eIdx);
+ scvIndices_[0] = eIdx;
+
+ const auto& scvFaceIndices = fvGridGeometry().scvfIndicesOfScv(eIdx);
+ const auto& neighborVolVarIndices = fvGridGeometry().neighborVolVarIndices(eIdx);
+
+ typename FVGridGeometry::GeometryHelper geometryHelper(element, fvGridGeometry().gridView());
+
+ int scvfCounter = 0;
+ for (const auto& intersection : intersections(fvGridGeometry().gridView(), element))
+ {
+ const auto& scvfNeighborVolVarIndex = neighborVolVarIndices[scvfCounter];
+
+ if (intersection.neighbor() || intersection.boundary())
+ {
+ geometryHelper.updateLocalFace(fvGridGeometry().intersectionMapper(), intersection);
+ std::vector<IndexType> scvIndices{eIdx, scvfNeighborVolVarIndex};
+ scvfs_.emplace_back(intersection,
+ intersection.geometry(),
+ scvFaceIndices[scvfCounter],
+ scvIndices,
+ geometryHelper);
+ scvfIndices_.emplace_back(scvFaceIndices[scvfCounter]);
+ scvfCounter++;
+ }
+ }
+ }
+
+ //! create the necessary scvs and scvfs of the neighbor elements to the bound elements
+ void makeNeighborGeometries_(const Element& element, const IndexType eIdx)
+ {
+ // using ScvfGridIndexStorage = typename SubControlVolumeFace::Traits::GridIndexStorage;
+
+ // create the neighbor scv
+ neighborScvs_.emplace_back(element.geometry(), eIdx);
+ neighborScvIndices_.push_back(eIdx);
+
+ const auto& scvFaceIndices = fvGridGeometry().scvfIndicesOfScv(eIdx);
+ const auto& neighborVolVarIndices = fvGridGeometry().neighborVolVarIndices(eIdx);
+
+ typename FVGridGeometry::GeometryHelper geometryHelper(element, fvGridGeometry().gridView());
+
+ int scvfCounter = 0;
+ for (const auto& intersection : intersections(fvGridGeometry().gridView(), element))
+ {
+ const auto& scvfNeighborVolVarIndex = neighborVolVarIndices[scvfCounter];
+ geometryHelper.updateLocalFace(fvGridGeometry().intersectionMapper(), intersection);
+
+ if (intersection.neighbor())
+ {
+ // only create subcontrol faces where the outside element is the bound element
+ if (intersection.outside() == *elementPtr_)
+ {
+ std::vector<IndexType> scvIndices{eIdx, scvfNeighborVolVarIndex};
+ neighborScvfs_.emplace_back(intersection,
+ intersection.geometry(),
+ scvFaceIndices[scvfCounter],
+ scvIndices,
+ geometryHelper);
+
+ neighborScvfIndices_.push_back(scvFaceIndices[scvfCounter]);
+ }
+ scvfCounter++;
+ }
+ else if (intersection.boundary())
+ scvfCounter++;
+ }
+ }
+
+ const IndexType findLocalIndex_(const IndexType idx,
+ const std::vector<IndexType>& indices) const
+ {
+ auto it = std::find(indices.begin(), indices.end(), idx);
+ assert(it != indices.end() && "Could not find the scv/scvf! Make sure to properly bind this class!");
+ return std::distance(indices.begin(), it);
+ }
+
+ //! Clear all local data
+ void clear_()
+ {
+ scvfIndices_.clear();
+ scvfs_.clear();
+
+ neighborScvIndices_.clear();
+ neighborScvfIndices_.clear();
+ neighborScvs_.clear();
+ neighborScvfs_.clear();
+ }
+
+ const Element* elementPtr_; //!< the element to which this fvgeometry is bound
+ const FVGridGeometry* fvGridGeometryPtr_; //!< the grid fvgeometry
+
+ // local storage after binding an element
+ std::array<IndexType, 1> scvIndices_;
+ std::array<SubControlVolume, 1> scvs_;
+
+ std::vector<IndexType> scvfIndices_;
+ std::vector<SubControlVolumeFace> scvfs_;
+
+ std::vector<IndexType> neighborScvIndices_;
+ std::vector<SubControlVolume> neighborScvs_;
+
+ std::vector<IndexType> neighborScvfIndices_;
+ std::vector<SubControlVolumeFace> neighborScvfs_;
+};
+
+
} // end namespace
#endif
diff --git a/dumux/discretization/staggered/fvgridgeometry.hh b/dumux/discretization/staggered/fvgridgeometry.hh
index 3fde8bfb994c58a68b6478a67f7be8d56ab1a9e9..dde1e9930ab4873bd4690364e6ca4692a0a5868f 100644
--- a/dumux/discretization/staggered/fvgridgeometry.hh
+++ b/dumux/discretization/staggered/fvgridgeometry.hh
@@ -223,7 +223,7 @@ public:
{
// Check if the overlap size is what we expect
if (!CheckOverlapSize<DiscretizationMethod::staggered>::isValid(gridView))
- DUNE_THROW(Dune::InvalidStateException, "The satggered discretization method needs at least an overlap of 1 for parallel computations. "
+ DUNE_THROW(Dune::InvalidStateException, "The staggered discretization method needs at least an overlap of 1 for parallel computations. "
<< " Set the parameter \"Grid.Overlap\" in the input file.");
}
@@ -314,8 +314,7 @@ public:
intersection.geometry(),
scvfIdx,
std::vector<IndexType>({eIdx, nIdx}),
- geometryHelper
- );
+ geometryHelper);
localToGlobalScvfIndices_[eIdx][localFaceIndex] = scvfIdx;
scvfsIndexSet.push_back(scvfIdx++);
}
@@ -326,8 +325,7 @@ public:
intersection.geometry(),
scvfIdx,
std::vector<IndexType>({eIdx, this->gridView().size(0) + numBoundaryScvf_++}),
- geometryHelper
- );
+ geometryHelper);
localToGlobalScvfIndices_[eIdx][localFaceIndex] = scvfIdx;
scvfsIndexSet.push_back(scvfIdx++);
}
@@ -421,8 +419,208 @@ private:
*/
template<class GV, class Traits>
class StaggeredFVGridGeometry<GV, false, Traits>
+: public BaseFVGridGeometry<StaggeredFVGridGeometry<GV, false, Traits>, GV, Traits>
{
- // TODO: implement without caching
+ using ThisType = StaggeredFVGridGeometry<GV, false, Traits>;
+ using ParentType = BaseFVGridGeometry<ThisType, GV, Traits>;
+ using IndexType = typename GV::IndexSet::IndexType;
+ using Element = typename GV::template Codim<0>::Entity;
+
+ using IntersectionMapper = typename Traits::IntersectionMapper;
+ using ConnectivityMap = typename Traits::template ConnectivityMap<ThisType>;
+
+public:
+ //! export discretization method
+ static constexpr DiscretizationMethod discMethod = DiscretizationMethod::staggered;
+
+ using GeometryHelper = typename Traits::GeometryHelper;
+
+ //! export the type of the fv element geometry (the local view type)
+ using LocalView = typename Traits::template LocalView<ThisType, false>;
+ //! 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 grid view type
+ using GridView = GV;
+ //! export the dof type indices
+ using DofTypeIndices = typename Traits::DofTypeIndices;
+
+ //! return a integral constant for cell center dofs
+ static constexpr auto cellCenterIdx()
+ { return typename DofTypeIndices::CellCenterIdx{}; }
+
+ //! return a integral constant for face dofs
+ static constexpr auto faceIdx()
+ { return typename DofTypeIndices::FaceIdx{}; }
+
+ using CellCenterFVGridGeometryType = CellCenterFVGridGeometry<ThisType>;
+ using FaceFVGridGeometryType = FaceFVGridGeometry<ThisType>;
+
+ using FVGridGeometryTuple = std::tuple< CellCenterFVGridGeometry<ThisType>, FaceFVGridGeometry<ThisType> >;
+
+ //! Constructor
+ StaggeredFVGridGeometry(const GridView& gridView)
+ : ParentType(gridView)
+ , intersectionMapper_(gridView)
+ {
+ // Check if the overlap size is what we expect
+ if (!CheckOverlapSize<DiscretizationMethod::staggered>::isValid(gridView))
+ DUNE_THROW(Dune::InvalidStateException, "The staggered discretization method needs at least an overlap of 1 for parallel computations. "
+ << " Set the parameter \"Grid.Overlap\" in the input file.");
+ }
+
+ //! update all fvElementGeometries (do this again after grid adaption)
+ void update()
+ {
+ // clear containers (necessary after grid refinement)
+ scvfIndicesOfScv_.clear();
+ intersectionMapper_.update();
+ neighborVolVarIndices_.clear();
+
+ numScvs_ = numCellCenterDofs();
+ numScvf_ = 0;
+ numBoundaryScvf_ = 0;
+ scvfIndicesOfScv_.resize(numScvs_);
+ localToGlobalScvfIndices_.resize(numScvs_);
+ neighborVolVarIndices_.resize(numScvs_);
+
+ // Build the scvs and scv faces
+ for (const auto& element : elements(this->gridView()))
+ {
+ auto eIdx = this->elementMapper().index(element);
+
+ // the element-wise index sets for finite volume geometry
+ auto numLocalFaces = intersectionMapper_.numFaces(element);
+ std::vector<IndexType> scvfsIndexSet;
+ scvfsIndexSet.reserve(numLocalFaces);
+ localToGlobalScvfIndices_[eIdx].resize(numLocalFaces);
+
+ std::vector<IndexType> neighborVolVarIndexSet;
+ neighborVolVarIndexSet.reserve(numLocalFaces);
+
+ for (const auto& intersection : intersections(this->gridView(), element))
+ {
+ const auto localFaceIndex = intersection.indexInInside();
+ localToGlobalScvfIndices_[eIdx][localFaceIndex] = numScvf_;
+ scvfsIndexSet.push_back(numScvf_++);
+
+ if (intersection.neighbor())
+ {
+ const auto nIdx = this->elementMapper().index(intersection.outside());
+ neighborVolVarIndexSet.emplace_back(nIdx);
+ }
+ else
+ neighborVolVarIndexSet.emplace_back(numScvs_ + numBoundaryScvf_++);
+ }
+
+ // Save the scvf indices belonging to this scv to build up fv element geometries fast
+ scvfIndicesOfScv_[eIdx] = scvfsIndexSet;
+ neighborVolVarIndices_[eIdx] = neighborVolVarIndexSet;
+ }
+
+ // build the connectivity map for an effecient assembly
+ connectivityMap_.update(*this);
+ }
+
+ //! The total number of sub control volumes
+ std::size_t numScv() const
+ {
+ return numScvs_;
+ }
+
+ //! The total number of sub control volume faces
+ std::size_t numScvf() const
+ {
+ return numScvf_;
+ }
+
+ //! The total number of boundary sub control volume faces
+ std::size_t numBoundaryScvf() const
+ {
+ return numBoundaryScvf_;
+ }
+
+ //! The total number of intersections
+ std::size_t numIntersections() const
+ {
+ return intersectionMapper_.numIntersections();
+ }
+
+ //! the total number of dofs
+ std::size_t numDofs() const
+ { return numCellCenterDofs() + numFaceDofs(); }
+
+ std::size_t numCellCenterDofs() const
+ { return this->gridView().size(0); }
+
+ std::size_t numFaceDofs() const
+ { return this->gridView().size(1); }
+
+ const std::vector<IndexType>& scvfIndicesOfScv(IndexType scvIdx) const
+ { return scvfIndicesOfScv_[scvIdx]; }
+
+ IndexType localToGlobalScvfIndex(IndexType eIdx, IndexType localScvfIdx) const
+ {
+ return localToGlobalScvfIndices_[eIdx][localScvfIdx];
+ }
+
+ /*!
+ * \brief Returns the connectivity map of which dofs have derivatives with respect
+ * to a given dof.
+ */
+ const ConnectivityMap &connectivityMap() const
+ { return connectivityMap_; }
+
+ //! Returns a pointer the cell center specific auxiliary class. Required for the multi-domain FVAssembler's ctor.
+ std::unique_ptr<CellCenterFVGridGeometry<ThisType>> cellCenterFVGridGeometryPtr() const
+ {
+ return std::make_unique<CellCenterFVGridGeometry<ThisType>>(this);
+ }
+
+ //! Returns a pointer the face specific auxiliary class. Required for the multi-domain FVAssembler's ctor.
+ std::unique_ptr<FaceFVGridGeometry<ThisType>> faceFVGridGeometryPtr() const
+ {
+ return std::make_unique<FaceFVGridGeometry<ThisType>>(this);
+ }
+
+ //! Return a copy of the cell center specific auxiliary class.
+ CellCenterFVGridGeometry<ThisType> cellCenterFVGridGeometry() const
+ {
+ return CellCenterFVGridGeometry<ThisType>(this);
+ }
+
+ //! Return a copy of the face specific auxiliary class.
+ FaceFVGridGeometry<ThisType> faceFVGridGeometry() const
+ {
+ return FaceFVGridGeometry<ThisType>(this);
+ }
+
+ //! Return a reference to the intersection mapper
+ const IntersectionMapper& intersectionMapper() const
+ {
+ return intersectionMapper_;
+ }
+
+ //! Return the neighbor volVar indices for all scvfs in the scv with index scvIdx
+ const std::vector<IndexType>& neighborVolVarIndices(IndexType scvIdx) const
+ { return neighborVolVarIndices_[scvIdx]; }
+
+private:
+
+ //! Information on the global number of geometries
+ IndexType numScvs_;
+ IndexType numScvf_;
+ IndexType numBoundaryScvf_;
+ std::vector<std::vector<IndexType>> localToGlobalScvfIndices_;
+ std::vector<std::vector<IndexType>> neighborVolVarIndices_;
+
+ // mappers
+ ConnectivityMap connectivityMap_;
+ IntersectionMapper intersectionMapper_;
+
+ //! vectors that store the global data
+ std::vector<std::vector<IndexType>> scvfIndicesOfScv_;
};
} // end namespace