diff --git a/dumux/common/properties.hh b/dumux/common/properties.hh
index 4899e913ebb381a2709c9706823b7a343c2ba63b..cb6e6de4c836c9237fa0ef1b280063f831d9e457 100644
--- a/dumux/common/properties.hh
+++ b/dumux/common/properties.hh
@@ -107,7 +107,7 @@ NEW_PROP_TAG(MpfaMethod); //!< Specifies the mpfa metho
NEW_PROP_TAG(MpfaHelper); //!< A Helper class depending on the mpfa method and grid dimension
NEW_PROP_TAG(PrimaryInteractionVolume); //!< The primary interaction volume type
NEW_PROP_TAG(SecondaryInteractionVolume); //!< The secondary interaction volume type used e.g. on the boundaries
-
+NEW_PROP_TAG(DualGridNodalIndexSet); //!< The type used for the nodal index sets of the dual grid
/////////////////////////////////////////////////////////////
// Properties used by models involving flow in porous media:
diff --git a/dumux/discretization/cellcentered/mpfa/dualgridindexset.hh b/dumux/discretization/cellcentered/mpfa/dualgridindexset.hh
index 52817810a3787d72f3259f47afd5935c80cae7d2..a93c7ce63970a589ca92cdf883dd9547556e1ee2 100644
--- a/dumux/discretization/cellcentered/mpfa/dualgridindexset.hh
+++ b/dumux/discretization/cellcentered/mpfa/dualgridindexset.hh
@@ -28,6 +28,8 @@
#include <vector>
#include <algorithm>
+#include <dune/common/reservedvector.hh>
+
namespace Dumux
{
@@ -39,20 +41,40 @@ namespace Dumux
* \tparam GI The type used for indices on the grid
* \tparam LI The type used for indexing in interaction volumes
* \tparam dim The dimension of the grid
+ * \tparam maxE The maximum admissible number of elements around vertices.
+ * \tparam maxB The maximum admissible number of branches on intersections.
+ * This is only to be specified for network grids and defaults to 1
+ * for normal grids.
*/
-template< class GI, class LI, int dim >
-class DualGridNodalIndexSet
+template< class GI, class LI, int dim, int maxE, int maxB = 2 >
+class CCMpfaDualGridNodalIndexSet
{
+ using DimIndexVector = Dune::ReservedVector<GI, dim>;
+
public:
+ //! Export the used index types
using GridIndexType = GI;
using LocalIndexType = LI;
- // we use dynamic containers to store indices here
- using GridIndexContainer = std::vector<GridIndexType>;
- using LocalIndexContainer = std::vector<LocalIndexType>;
+ //! Export the specified maximum admissible sizes
+ static constexpr int dimension = dim;
+ static constexpr int maxBranches = maxB;
+ static constexpr int maxNumElementsAtNode = maxE*(maxBranches-1);
+ static constexpr int maxNumScvfsAtNode = maxNumElementsAtNode*dim;
+
+ //! Export the stencil types used
+ using GridStencilType = Dune::ReservedVector<GridIndexType, maxNumElementsAtNode>;
+ using LocalStencilType = Dune::ReservedVector<LocalIndexType, maxNumElementsAtNode>;
+
+ //! Export the type used for storing the global scvf indices at this node
+ using GridScvfStencilType = Dune::ReservedVector<GridIndexType, maxNumScvfsAtNode>;
+
+ //! Data structure to store the neighboring scv indices of an scvf (grid/local indices)
+ using ScvfNeighborIndexSet = Dune::ReservedVector<GridIndexType, maxBranches>;
+ using ScvfNeighborLocalIndexSet = Dune::ReservedVector<LocalIndexType, maxBranches>;
//! Constructor
- DualGridNodalIndexSet() : numBoundaryScvfs_(0) {}
+ CCMpfaDualGridNodalIndexSet() : numBoundaryScvfs_(0) {}
//! Inserts data for a given scvf
template<typename SubControlVolumeFace>
@@ -78,10 +100,10 @@ public:
const LocalIndexType curScvfLocalIdx = scvfIndices_.size();
// add grid scvf data
- GridIndexContainer scvIndices;
- scvIndices.reserve(outsideScvIndices.size()+1);
+ ScvfNeighborIndexSet scvIndices;
scvIndices.push_back(insideScvIdx);
- scvIndices.insert(scvIndices.end(), outsideScvIndices.begin(), outsideScvIndices.end());
+ for (const auto& outsideIdx : outsideScvIndices)
+ scvIndices.push_back(outsideIdx);
// if scvf is on boundary, increase counter
if (boundary)
@@ -90,7 +112,7 @@ public:
// insert data on the new scv
scvfIndices_.push_back(scvfIdx);
scvfIsOnBoundary_.push_back(boundary);
- scvfNeighborScvIndices_.emplace_back( std::move(scvIndices) );
+ scvfNeighborScvIndices_.push_back(scvIndices);
// if entry for the inside scv exists, append scvf local index, create otherwise
auto it = std::find( scvIndices_.begin(), scvIndices_.end(), insideScvIdx );
@@ -98,10 +120,8 @@ public:
localScvfIndicesInScv_[ std::distance(scvIndices_.begin(), it) ].push_back(curScvfLocalIdx);
else
{
- LocalIndexContainer localScvfIndices;
- localScvfIndices.reserve(dim);
- localScvfIndices.push_back(curScvfLocalIdx);
- localScvfIndicesInScv_.emplace_back( std::move(localScvfIndices) );
+ localScvfIndicesInScv_.push_back({});
+ localScvfIndicesInScv_.back().push_back(curScvfLocalIdx);
scvIndices_.push_back(insideScvIdx);
}
}
@@ -116,10 +136,10 @@ public:
std::size_t numBoundaryScvfs() const { return numBoundaryScvfs_; }
//! returns the global scv indices connected to this dual grid node
- const GridIndexContainer& globalScvIndices() const { return scvIndices_; }
+ const GridStencilType& globalScvIndices() const { return scvIndices_; }
//! returns the global scvf indices connected to this dual grid node
- const GridIndexContainer& globalScvfIndices() const { return scvfIndices_; }
+ const GridScvfStencilType& globalScvfIndices() const { return scvfIndices_; }
//! returns whether or not the i-th scvf is on a domain boundary
bool scvfIsOnBoundary(unsigned int i) const { return scvfIsOnBoundary_[i]; }
@@ -145,39 +165,40 @@ public:
}
//! returns the indices of the neighboring scvs of the i-th scvf
- const GridIndexContainer& neighboringScvIndices(unsigned int i) const
+ const ScvfNeighborIndexSet& neighboringScvIndices(unsigned int i) const
{ return scvfNeighborScvIndices_[i]; }
private:
- GridIndexContainer scvIndices_; //!< The indices of the scvs around a dual grid node
- std::vector<LocalIndexContainer> localScvfIndicesInScv_; //!< Maps to each scv a list of scvf indices embedded in it
+ GridStencilType scvIndices_; //!< The indices of the scvs around a dual grid node
+ Dune::ReservedVector<DimIndexVector, maxNumElementsAtNode> localScvfIndicesInScv_; //!< Maps to each scv a list of scvf indices embedded in it
+
+ GridScvfStencilType scvfIndices_; //!< the indices of the scvfs around a dual grid node
+ std::size_t numBoundaryScvfs_; //!< stores how many boundary scvfs are embedded in this dual grid node
+ Dune::ReservedVector<bool, maxNumScvfsAtNode> scvfIsOnBoundary_; //!< Maps to each scvf a boolean to indicate if it is on the boundary
- GridIndexContainer scvfIndices_; //!< the indices of the scvfs around a dual grid node
- std::vector<bool> scvfIsOnBoundary_; //!< Maps to each scvf a boolean to indicate if it is on the boundary
- std::vector<GridIndexContainer> scvfNeighborScvIndices_; //!< The neighboring scvs for the scvfs (order: 0 - inside, 1..n - outside)
- std::size_t numBoundaryScvfs_; //!< stores how many boundary scvfs are embedded in this dual grid node
+ //! The neighboring scvs for the scvfs (order: 0 - inside, 1..n - outside)
+ Dune::ReservedVector<ScvfNeighborIndexSet, maxNumScvfsAtNode> scvfNeighborScvIndices_;
};
/*!
* \ingroup CCMpfaDiscretization
* \brief Class for the index sets of the dual grid in mpfa schemes.
*
- * \tparam GridIndexType The type used for indices on the grid
- * \tparam LocalIndexType The type used for indexing in interaction volumes
- * \tparam dim The dimension of the grid
+ * \tparam NI The type used for the nodal index sets.
*/
-template< class GridIndexType, class LocalIndexType, int dim >
+template< class NI >
class CCMpfaDualGridIndexSet
{
public:
- using NodalIndexSet = DualGridNodalIndexSet< GridIndexType, LocalIndexType, dim >;
+ using NodalIndexSet = NI;
+ using GridIndexType = typename NodalIndexSet::GridIndexType;
//! Default constructor should not be used
CCMpfaDualGridIndexSet() = delete;
//! Constructor taking a grid view
template< class GridView >
- CCMpfaDualGridIndexSet(const GridView& gridView) : nodalIndexSets_(gridView.size(dim)) {}
+ CCMpfaDualGridIndexSet(const GridView& gridView) : nodalIndexSets_(gridView.size(NodalIndexSet::dimension)) {}
//! Access with an scvf
template< class SubControlVolumeFace >
diff --git a/dumux/discretization/cellcentered/mpfa/fvelementgeometry.hh b/dumux/discretization/cellcentered/mpfa/fvelementgeometry.hh
index f3fbca56476ed92c90b1194dbb6e5c5f1e1e73ff..44b8f899fe6fdb14221087ea24a8bd92bf36fb6b 100644
--- a/dumux/discretization/cellcentered/mpfa/fvelementgeometry.hh
+++ b/dumux/discretization/cellcentered/mpfa/fvelementgeometry.hh
@@ -59,16 +59,22 @@ class CCMpfaFVElementGeometry<TypeTag, true>
using ThisType = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
using Element = typename GridView::template Codim<0>::Entity;
-
using GridIndexType = typename GridView::IndexSet::IndexType;
+
+ static constexpr int dim = GridView::dimension;
+
+public:
+ //! export type of subcontrol volume
using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
+ //! export type of subcontrol volume face
using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
+ //! export type of finite volume grid geometry
using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
+ //! the maximum number of scvs per element
+ static constexpr std::size_t maxNumElementScvs = 1;
+ //! the maximum number of scvfs per element (use cubes for maximum)
+ static constexpr std::size_t maxNumElementScvfs = dim == 3 ? 24 : 8;
- using ScvIterator = Dumux::ScvIterator<SubControlVolume, std::vector<GridIndexType>, ThisType>;
- using ScvfIterator = Dumux::ScvfIterator<SubControlVolumeFace, std::vector<GridIndexType>, ThisType>;
-
-public:
//! Constructor
CCMpfaFVElementGeometry(const FVGridGeometry& fvGridGeometry)
: fvGridGeometryPtr_(&fvGridGeometry) {}
@@ -97,9 +103,10 @@ public:
//! 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<ScvIterator>
+ friend inline Dune::IteratorRange< ScvIterator<SubControlVolume, std::array<GridIndexType, 1>, ThisType> >
scvs(const CCMpfaFVElementGeometry& fvGeometry)
{
+ using ScvIterator = Dumux::ScvIterator<SubControlVolume, std::array<GridIndexType, 1>, ThisType>;
return Dune::IteratorRange<ScvIterator>(ScvIterator(fvGeometry.scvIndices_.begin(), fvGeometry),
ScvIterator(fvGeometry.scvIndices_.end(), fvGeometry));
}
@@ -109,11 +116,12 @@ public:
//! 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<ScvfIterator>
+ friend inline Dune::IteratorRange< ScvfIterator<SubControlVolumeFace, std::vector<GridIndexType>, ThisType> >
scvfs(const CCMpfaFVElementGeometry& fvGeometry)
{
const auto& g = fvGeometry.fvGridGeometry();
const auto scvIdx = fvGeometry.scvIndices_[0];
+ using ScvfIterator = Dumux::ScvfIterator<SubControlVolumeFace, std::vector<GridIndexType>, ThisType>;
return Dune::IteratorRange<ScvfIterator>(ScvfIterator(g.scvfIndicesOfScv(scvIdx).begin(), fvGeometry),
ScvfIterator(g.scvfIndicesOfScv(scvIdx).end(), fvGeometry));
}
@@ -139,7 +147,7 @@ public:
//! Bind only element-local
void bindElement(const Element& element)
{
- scvIndices_ = std::vector<GridIndexType>({fvGridGeometry().elementMapper().index(element)});
+ scvIndices_[0] = fvGridGeometry().elementMapper().index(element);
}
//! The global finite volume geometry we are a restriction of
@@ -148,7 +156,7 @@ public:
private:
- std::vector<GridIndexType> scvIndices_;
+ std::array<GridIndexType, 1> scvIndices_;
const FVGridGeometry* fvGridGeometryPtr_;
};
@@ -168,12 +176,6 @@ class CCMpfaFVElementGeometry<TypeTag, false>
using GridIndexType = typename GridView::IndexSet::IndexType;
using MpfaHelper = typename GET_PROP_TYPE(TypeTag, MpfaHelper);
- using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
- using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
- using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
-
- using ScvIterator = Dumux::ScvIterator<SubControlVolume, std::vector<GridIndexType>, ThisType>;
- using ScvfIterator = Dumux::ScvfIterator<SubControlVolumeFace, std::vector<GridIndexType>, ThisType>;
static const int dim = GridView::dimension;
static const int dimWorld = GridView::dimensionworld;
@@ -182,6 +184,17 @@ class CCMpfaFVElementGeometry<TypeTag, false>
using GlobalPosition = Dune::FieldVector<Scalar, dimWorld>;
public:
+ //! export type of subcontrol volume
+ using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
+ //! export type of subcontrol volume face
+ using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
+ //! export type of finite volume grid geometry
+ using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
+ //! the maximum number of scvs per element
+ static constexpr std::size_t maxNumElementScvs = 1;
+ //! the maximum number of scvfs per element (use cubes for maximum)
+ static constexpr std::size_t maxNumElementScvfs = dim == 3 ? 24 : 8;
+
//! Constructor
CCMpfaFVElementGeometry(const FVGridGeometry& fvGridGeometry)
: fvGridGeometryPtr_(&fvGridGeometry) {}
@@ -234,11 +247,11 @@ public:
//! 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::vector<SubControlVolume>::const_iterator>
+ friend inline Dune::IteratorRange<typename std::array<SubControlVolume, 1>::const_iterator>
scvs(const ThisType& g)
{
- using Iter = typename std::vector<SubControlVolume>::const_iterator;
- return Dune::IteratorRange<Iter>(g.scvs_.begin(), g.scvs_.end());
+ 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
@@ -249,8 +262,8 @@ public:
friend inline Dune::IteratorRange<typename std::vector<SubControlVolumeFace>::const_iterator>
scvfs(const ThisType& g)
{
- using Iter = typename std::vector<SubControlVolumeFace>::const_iterator;
- return Dune::IteratorRange<Iter>(g.scvfs_.begin(), g.scvfs_.end());
+ 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
@@ -336,8 +349,8 @@ private:
{
// make the scv
const auto eIdx = fvGridGeometry().elementMapper().index(element);
- scvs_.emplace_back(element.geometry(), eIdx);
- scvIndices_.emplace_back(eIdx);
+ scvs_[0] = SubControlVolume(element.geometry(), eIdx);
+ scvIndices_[0] = eIdx;
// get data on the scv faces
const auto& scvFaceIndices = fvGridGeometry().scvfIndicesOfScv(eIdx);
@@ -617,9 +630,7 @@ private:
//! clear all containers
void clear()
{
- scvIndices_.clear();
scvfIndices_.clear();
- scvs_.clear();
scvfs_.clear();
neighborScvIndices_.clear();
@@ -634,9 +645,9 @@ private:
const FVGridGeometry* fvGridGeometryPtr_;
// local storage after binding an element
- std::vector<GridIndexType> scvIndices_;
+ std::array<GridIndexType, 1> scvIndices_;
std::vector<GridIndexType> scvfIndices_;
- std::vector<SubControlVolume> scvs_;
+ std::array<SubControlVolume, 1> scvs_;
std::vector<SubControlVolumeFace> scvfs_;
std::vector<GridIndexType> neighborScvIndices_;
diff --git a/dumux/discretization/cellcentered/mpfa/fvgridgeometry.hh b/dumux/discretization/cellcentered/mpfa/fvgridgeometry.hh
index 0b1e09f5d9722e5a3babc41ab78ad4a2b3adb46d..dfba3592e209700e45d4642e397f517b5433f180 100644
--- a/dumux/discretization/cellcentered/mpfa/fvgridgeometry.hh
+++ b/dumux/discretization/cellcentered/mpfa/fvgridgeometry.hh
@@ -64,6 +64,7 @@ class CCMpfaFVGridGeometry<TypeTag, true> : public BaseFVGridGeometry<TypeTag>
using MpfaHelper = typename GET_PROP_TYPE(TypeTag, MpfaHelper);
using PrimaryInteractionVolume = typename GET_PROP_TYPE(TypeTag, PrimaryInteractionVolume);
+ using DualGridNodalIndexSet = typename GET_PROP_TYPE(TypeTag, DualGridNodalIndexSet);
using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
@@ -164,7 +165,7 @@ public:
const auto isGhostVertex = MpfaHelper::findGhostVertices(this->gridView(), this->vertexMapper());
// instantiate the dual grid index set (to be used for construction of interaction volumes)
- CCMpfaDualGridIndexSet<GridIndexType, LocalIndexType, dim> dualIdSet(this->gridView());
+ CCMpfaDualGridIndexSet< DualGridNodalIndexSet > dualIdSet(this->gridView());
// Build the SCVs and SCV faces
GridIndexType scvfIdx = 0;
@@ -405,6 +406,7 @@ class CCMpfaFVGridGeometry<TypeTag, false> : public BaseFVGridGeometry<TypeTag>
using MpfaHelper = typename GET_PROP_TYPE(TypeTag, MpfaHelper);
using PrimaryInteractionVolume = typename GET_PROP_TYPE(TypeTag, PrimaryInteractionVolume);
+ using DualGridNodalIndexSet = typename GET_PROP_TYPE(TypeTag, DualGridNodalIndexSet);
using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
@@ -505,7 +507,7 @@ public:
isGhostVertex_ = MpfaHelper::findGhostVertices(this->gridView(), this->vertexMapper());
// instantiate the dual grid index set (to be used for construction of interaction volumes)
- CCMpfaDualGridIndexSet<GridIndexType, LocalIndexType, dim> dualIdSet(this->gridView());
+ CCMpfaDualGridIndexSet< DualGridNodalIndexSet > dualIdSet(this->gridView());
// Build the SCVs and SCV faces
numScvf_ = 0;
diff --git a/dumux/discretization/cellcentered/mpfa/gridinteractionvolumeindexsets.hh b/dumux/discretization/cellcentered/mpfa/gridinteractionvolumeindexsets.hh
index 5a7a7882f0ddfcfe07e9b4ed0005244d9e0ef336..7ace724febac7f73d204653a76a66ea72fc4ac84 100644
--- a/dumux/discretization/cellcentered/mpfa/gridinteractionvolumeindexsets.hh
+++ b/dumux/discretization/cellcentered/mpfa/gridinteractionvolumeindexsets.hh
@@ -43,6 +43,7 @@ class CCMpfaGridInteractionVolumeIndexSets
using GridIndexType = typename GridView::IndexSet::IndexType;
using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
+ using DualGridNodalIndexSet = typename GET_PROP_TYPE(TypeTag, DualGridNodalIndexSet);
using PrimaryIV = typename GET_PROP_TYPE(TypeTag, PrimaryInteractionVolume);
using PrimaryIVIndexSet = typename PrimaryIV::Traits::IndexSet;
@@ -51,7 +52,7 @@ class CCMpfaGridInteractionVolumeIndexSets
static constexpr int dim = GridView::dimension;
using LocalIndexType = typename PrimaryIV::Traits::LocalIndexType;
- using DualGridIndexSet = CCMpfaDualGridIndexSet< GridIndexType, LocalIndexType, dim>;
+ using DualGridIndexSet = CCMpfaDualGridIndexSet< DualGridNodalIndexSet >;
public:
/*!
diff --git a/dumux/discretization/cellcentered/mpfa/helper.hh b/dumux/discretization/cellcentered/mpfa/helper.hh
index f129138540b5be817f3789ac4880a4f75802f17c..5a32421a9ed42b371e90f08cf5fd5fae8155fb0c 100644
--- a/dumux/discretization/cellcentered/mpfa/helper.hh
+++ b/dumux/discretization/cellcentered/mpfa/helper.hh
@@ -654,8 +654,8 @@ public:
{ return !std::is_same<PrimaryInteractionVolume, SecondaryInteractionVolume>::value; }
//! returns whether or not a value exists in a vector
- template<typename V1, typename V2>
- static bool vectorContainsValue(const std::vector<V1>& vector, const V2 value)
+ template<typename V, typename T>
+ static bool vectorContainsValue(const V& vector, const T value)
{ return std::find(vector.begin(), vector.end(), value) != vector.end(); }
};
diff --git a/dumux/discretization/cellcentered/mpfa/omethod/interactionvolume.hh b/dumux/discretization/cellcentered/mpfa/omethod/interactionvolume.hh
index 298ec506bd6461955c8b68fe6ee4c25d3417d8ab..0a8f754d6fdbc29fbf4dac6e41cc65bcc1cedf27 100644
--- a/dumux/discretization/cellcentered/mpfa/omethod/interactionvolume.hh
+++ b/dumux/discretization/cellcentered/mpfa/omethod/interactionvolume.hh
@@ -51,7 +51,8 @@ template< class TypeTag >
struct CCMpfaODefaultInteractionVolumeTraits
{
private:
- using GridIndexType = typename GET_PROP_TYPE(TypeTag, GridView)::IndexSet::IndexType;
+ using NodalIndexSet = typename GET_PROP_TYPE(TypeTag, DualGridNodalIndexSet);
+ using GridIndexType = typename NodalIndexSet::GridIndexType;
static constexpr int dim = GET_PROP_TYPE(TypeTag, GridView)::dimension;
static constexpr int dimWorld = GET_PROP_TYPE(TypeTag, GridView)::dimensionworld;
@@ -69,9 +70,9 @@ public:
//! export the type of the mpfa helper class
using MpfaHelper = typename GET_PROP_TYPE(TypeTag, MpfaHelper);
//! export the type used for local indices
- using LocalIndexType = std::uint8_t;
+ using LocalIndexType = typename NodalIndexSet::LocalIndexType;
//! export the type for the interaction volume index set
- using IndexSet = CCMpfaOInteractionVolumeIndexSet< DualGridNodalIndexSet<GridIndexType, LocalIndexType, dim> >;
+ using IndexSet = CCMpfaOInteractionVolumeIndexSet< NodalIndexSet >;
//! export the type of interaction-volume local scvs
using LocalScvType = CCMpfaOInteractionVolumeLocalScv< IndexSet, ScalarType, dim, dimWorld >;
//! export the type of interaction-volume local scvfs
@@ -83,7 +84,7 @@ public:
//! export the type used for iv-local vectors
using Vector = Dune::DynamicVector< ScalarType >;
//! export the type used for the iv-stencils
- using Stencil = std::vector< GridIndexType >;
+ using Stencil = typename NodalIndexSet::GridStencilType;
//! export the data handle type for this iv
using DataHandle = InteractionVolumeDataHandle< TypeTag, Matrix, Vector, LocalIndexType >;
};
@@ -122,10 +123,10 @@ class CCMpfaOInteractionVolume : public CCMpfaInteractionVolumeBase< CCMpfaOInte
using Stencil = typename Traits::Stencil;
//! For the o method, the interaction volume stencil can be taken directly
- //! from the nodal index set, which always uses dynamic types to be compatible
- //! on the boundaries and unstructured grids. Thus, we have to make sure that
+ //! from the nodal index set, which always uses Dune::ReservedVector with a maximum
+ //! admissible number of elements per node. Thus, we have to make sure that
//! the type set for the stencils in the traits is castable.
- static_assert( std::is_convertible<Stencil*, typename IndexSet::GridIndexContainer*>::value,
+ static_assert( std::is_convertible<Stencil*, typename IndexSet::GridStencilType*>::value,
"The o-method uses the (dynamic) nodal index set's stencil as the interaction volume stencil. "
"Using a different type is not permissive here." );
diff --git a/dumux/discretization/cellcentered/mpfa/omethod/interactionvolumeindexset.hh b/dumux/discretization/cellcentered/mpfa/omethod/interactionvolumeindexset.hh
index 84a26f7a4e4a253f905b1ca6405f6b1e2a2bdd65..1d63d804f259675a3619111b2056aa24dcd9430f 100644
--- a/dumux/discretization/cellcentered/mpfa/omethod/interactionvolumeindexset.hh
+++ b/dumux/discretization/cellcentered/mpfa/omethod/interactionvolumeindexset.hh
@@ -24,6 +24,8 @@
#ifndef DUMUX_DISCRETIZATION_MPFA_O_INTERACTIONVOLUME_INDEXSET_HH
#define DUMUX_DISCRETIZATION_MPFA_O_INTERACTIONVOLUME_INDEXSET_HH
+#include <dune/common/reservedvector.hh>
+
#include <dumux/discretization/cellcentered/mpfa/dualgridindexset.hh>
namespace Dumux
@@ -38,46 +40,27 @@ template< class DualGridNodalIndexSet >
class CCMpfaOInteractionVolumeIndexSet
{
public:
+ //! Export the type used for the nodal grid index sets
+ using NodalIndexSet = DualGridNodalIndexSet;
+
+ //! Export the types used for local/grid indices
using LocalIndexType = typename DualGridNodalIndexSet::LocalIndexType;
using GridIndexType = typename DualGridNodalIndexSet::GridIndexType;
- using LocalIndexContainer = typename DualGridNodalIndexSet::LocalIndexContainer;
- using GridIndexContainer = typename DualGridNodalIndexSet::GridIndexContainer;
-
- /*!
- * \brief The constructor
- * \note The actual type used for the nodal index sets might be different, as maybe
- * a different type for the local indexes is used. We therefore template this
- * constructor. However, a static assertion enforces you to use the same LocalIndexType
- * in the traits for both the secondary and the primary interaction volume traits.
- *
- * \tparam NodalIndexSet Possibly differing type for the DualGridNodalIndexSet
- */
- template< class NodalIndexSet >
- CCMpfaOInteractionVolumeIndexSet(const NodalIndexSet& nodalIndexSet)
- : nodalIndexSet_( static_cast<const DualGridNodalIndexSet&>(nodalIndexSet) )
- {
- // make sure the index types used are the same in order to avoid any losses due to type conversion
- static_assert(std::is_same<GridIndexType, typename NodalIndexSet::GridIndexType>::value,
- "Provided nodal index set does not use the same type for grid indices as the given template argument");
- static_assert(std::is_same<LocalIndexType, typename NodalIndexSet::LocalIndexType>::value,
- "Provided nodal index set does not use the same type for local indices as the given template argument");
+ // Export the types used for local/grid stencils
+ using LocalStencilType = typename DualGridNodalIndexSet::LocalStencilType;
+ using GridStencilType = typename DualGridNodalIndexSet::GridStencilType;
+ using GridScvfStencilType = typename DualGridNodalIndexSet::GridScvfStencilType;
+ //! Export the type used for the neighbor scv index sets of the scvfs
+ using ScvfNeighborLocalIndexSet = typename DualGridNodalIndexSet::ScvfNeighborLocalIndexSet;
+
+ //! The constructor
+ CCMpfaOInteractionVolumeIndexSet(const NodalIndexSet& nodalIndexSet) : nodalIndexSet_(nodalIndexSet)
+ {
// determine the number of iv-local faces for memory reservation
// note that this might be a vast overestimation on surface grids!
const auto numNodalScvfs = nodalIndexSet.numScvfs();
- const auto numBoundaryScvfs = nodalIndexSet.numBoundaryScvfs();
- const std::size_t numFaceEstimate = numBoundaryScvfs + (numNodalScvfs-numBoundaryScvfs)/2;
-
- // make sure we found a reasonable number of faces
- assert((numNodalScvfs-numBoundaryScvfs)%2 == 0);
-
- // index transformation from interaction-volume-local to node-local
- ivToNodeScvf_.reserve(numFaceEstimate);
- nodeToIvScvf_.resize(numNodalScvfs);
-
- // the local neighboring scv indices of the faces
- scvfNeighborScvLocalIndices_.reserve(numFaceEstimate);
// keeps track of which nodal scvfs have been handled already
std::vector<bool> isHandled(numNodalScvfs, false);
@@ -172,7 +155,7 @@ public:
{ return nodeToIvScvf_[ nodalIndexSet_.scvfIdxLocal(scvIdxLocal, i) ]; }
//! returns the local indices of the neighboring scvs of an scvf
- const LocalIndexContainer& neighboringLocalScvIndices(LocalIndexType ivLocalScvfIdx) const
+ const ScvfNeighborLocalIndexSet& neighboringLocalScvIndices(LocalIndexType ivLocalScvfIdx) const
{ return scvfNeighborScvLocalIndices_[ivLocalScvfIdx]; }
//! returns the number of faces in the interaction volume
@@ -182,10 +165,10 @@ public:
std::size_t numScvs() const { return nodalIndexSet_.numScvs(); }
//! returns the global scv indices connected to this dual grid node
- const GridIndexContainer& globalScvIndices() const { return nodalIndexSet_.globalScvIndices(); }
+ const GridStencilType& globalScvIndices() const { return nodalIndexSet_.globalScvIndices(); }
//! returns the global scvf indices connected to this dual grid node
- const GridIndexContainer& globalScvfIndices() const { return nodalIndexSet_.globalScvfIndices(); }
+ const GridScvfStencilType& globalScvfIndices() const { return nodalIndexSet_.globalScvfIndices(); }
private:
//! returns the local scv index to a given global scv index
@@ -199,11 +182,12 @@ private:
const DualGridNodalIndexSet& nodalIndexSet_;
std::size_t numFaces_;
- std::vector<LocalIndexType> ivToNodeScvf_;
- std::vector<LocalIndexType> nodeToIvScvf_;
+ Dune::ReservedVector< LocalIndexType, NodalIndexSet::maxNumScvfsAtNode > ivToNodeScvf_;
+ Dune::ReservedVector< LocalIndexType, NodalIndexSet::maxNumScvfsAtNode > nodeToIvScvf_;
+
// maps to each scvf a list of neighbouring scv indices
// ordering: 0 - inside scv idx; 1..n - outside scv indices
- std::vector< LocalIndexContainer > scvfNeighborScvLocalIndices_;
+ Dune::ReservedVector< ScvfNeighborLocalIndexSet, NodalIndexSet::maxNumScvfsAtNode > scvfNeighborScvLocalIndices_;
};
} // end namespace Dumux
diff --git a/dumux/discretization/cellcentered/mpfa/omethod/localsubcontrolentities.hh b/dumux/discretization/cellcentered/mpfa/omethod/localsubcontrolentities.hh
index 7ddded18c3abba12cc578707e0358606da02597d..510ace43c52a0393a4acfa1003eb7d3482409b6a 100644
--- a/dumux/discretization/cellcentered/mpfa/omethod/localsubcontrolentities.hh
+++ b/dumux/discretization/cellcentered/mpfa/omethod/localsubcontrolentities.hh
@@ -134,7 +134,7 @@ private:
template< class IvIndexSet >
struct CCMpfaOInteractionVolumeLocalScvf
{
- using LocalIndexContainer = typename IvIndexSet::LocalIndexContainer;
+ using ScvfNeighborLocalIndexSet = typename IvIndexSet::ScvfNeighborLocalIndexSet;
public:
// export index types
@@ -154,7 +154,7 @@ public:
*/
template< class SubControlVolumeFace >
CCMpfaOInteractionVolumeLocalScvf(const SubControlVolumeFace& scvf,
- const LocalIndexContainer& localScvIndices,
+ const ScvfNeighborLocalIndexSet& localScvIndices,
const LocalIndexType localDofIdx,
const bool isDirichlet)
: isDirichlet_(isDirichlet)
@@ -171,7 +171,7 @@ public:
GridIndexType globalScvfIndex() const { return scvfIdxGlobal_; }
//! Returns the local indices of the scvs neighboring this scvf
- const LocalIndexContainer& neighboringLocalScvIndices() const { return *neighborScvIndicesLocal_; }
+ const ScvfNeighborLocalIndexSet& neighboringLocalScvIndices() const { return *neighborScvIndicesLocal_; }
//! states if this is scvf is on a Dirichlet boundary
bool isDirichlet() const { return isDirichlet_; }
@@ -180,7 +180,7 @@ private:
bool isDirichlet_;
GridIndexType scvfIdxGlobal_;
LocalIndexType localDofIndex_;
- const LocalIndexContainer* neighborScvIndicesLocal_;
+ const ScvfNeighborLocalIndexSet* neighborScvIndicesLocal_;
};
} // end namespace Dumux
diff --git a/dumux/discretization/cellcentered/mpfa/omethod/staticinteractionvolume.hh b/dumux/discretization/cellcentered/mpfa/omethod/staticinteractionvolume.hh
index d4401e178678120e2d3c807aa683281c052f6a35..c57a82f3ac745683675b6674c03e764af349c7a6 100644
--- a/dumux/discretization/cellcentered/mpfa/omethod/staticinteractionvolume.hh
+++ b/dumux/discretization/cellcentered/mpfa/omethod/staticinteractionvolume.hh
@@ -52,7 +52,8 @@ template< class TypeTag, int localSize >
struct CCMpfaODefaultStaticInteractionVolumeTraits
{
private:
- using GridIndexType = typename GET_PROP_TYPE(TypeTag, GridView)::IndexSet::IndexType;
+ using NodalIndexSet = typename GET_PROP_TYPE(TypeTag, DualGridNodalIndexSet);
+ using GridIndexType = typename NodalIndexSet::GridIndexType;
static constexpr int dim = GET_PROP_TYPE(TypeTag, GridView)::dimension;
static constexpr int dimWorld = GET_PROP_TYPE(TypeTag, GridView)::dimensionworld;
@@ -70,9 +71,9 @@ public:
//! export the type of the mpfa helper class
using MpfaHelper = typename GET_PROP_TYPE(TypeTag, MpfaHelper);
//! export the type used for local indices
- using LocalIndexType = std::uint8_t;
+ using LocalIndexType = typename NodalIndexSet::LocalIndexType;
//! export the type for the interaction volume index set
- using IndexSet = CCMpfaOInteractionVolumeIndexSet< DualGridNodalIndexSet<GridIndexType, LocalIndexType, dim> >;
+ using IndexSet = CCMpfaOInteractionVolumeIndexSet< NodalIndexSet >;
//! export the type of interaction-volume local scvs
using LocalScvType = CCMpfaOInteractionVolumeLocalScv< IndexSet, ScalarType, dim, dimWorld >;
//! export the type of interaction-volume local scvfs
@@ -84,7 +85,7 @@ public:
//! export the type used for iv-local vectors
using Vector = Dune::FieldVector< ScalarType, localSize >;
//! export the type used for the iv-stencils
- using Stencil = std::vector< GridIndexType >;
+ using Stencil = typename NodalIndexSet::GridStencilType;
//! export the data handle type for this iv
using DataHandle = InteractionVolumeDataHandle< TypeTag, Matrix, Vector, LocalIndexType >;
};
@@ -132,7 +133,7 @@ class CCMpfaOStaticInteractionVolume : public CCMpfaInteractionVolumeBase< CCMpf
//! from the nodal index set, which always uses dynamic types to be compatible
//! on the boundaries and unstructured grids. Thus, we have to make sure that
//! the type set for the stencils in the traits is castable.
- static_assert( std::is_convertible<Stencil*, typename IndexSet::GridIndexContainer*>::value,
+ static_assert( std::is_convertible<Stencil*, typename IndexSet::GridStencilType*>::value,
"The o-method uses the (dynamic) nodal index set's stencil as the interaction volume stencil. "
"Using a different type is not permissive here." );
diff --git a/dumux/discretization/cellcentered/mpfa/properties.hh b/dumux/discretization/cellcentered/mpfa/properties.hh
index 1a5e96d98a389f00c0976233008b749cc5f4e1d0..1b23351b354a8c10f336069bbce11669a95c409e 100644
--- a/dumux/discretization/cellcentered/mpfa/properties.hh
+++ b/dumux/discretization/cellcentered/mpfa/properties.hh
@@ -47,10 +47,10 @@
#include <dumux/discretization/cellcentered/mpfa/elementvolumevariables.hh>
#include <dumux/discretization/cellcentered/mpfa/elementfluxvariablescache.hh>
#include <dumux/discretization/cellcentered/mpfa/subcontrolvolumeface.hh>
+#include <dumux/discretization/cellcentered/mpfa/dualgridindexset.hh>
#include <dumux/discretization/cellcentered/mpfa/helper.hh>
#include <dumux/discretization/cellcentered/mpfa/omethod/interactionvolume.hh>
-#include <dumux/discretization/cellcentered/mpfa/omethod/staticinteractionvolume.hh>
namespace Dumux
{
@@ -71,6 +71,44 @@ SET_PROP(CCMpfaModel, MpfaMethod)
static const MpfaMethods value = GET_PROP_TYPE(TypeTag, PrimaryInteractionVolume)::MpfaMethod;
};
+//! Set the maximum admissible number of branches per scvf
+SET_PROP(CCMpfaModel, MaxNumBranchesPerScvf)
+{
+private:
+ using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+ static constexpr int dim = GridView::dimension;
+ static constexpr int dimWorld = GridView::dimensionworld;
+
+public:
+ // Per default, we allow for 8 neighbors on network/surface grids
+ static const std::size_t value = dim < dimWorld ? 9 : 2;
+};
+
+//! Set the index set type used on the dual grid nodes
+SET_PROP(CCMpfaModel, DualGridNodalIndexSet)
+{
+private:
+ using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+
+ // Use the grid view's index type for grid indices
+ using GI = typename GridView::IndexSet::IndexType;
+
+ // per default, use uint8_t as iv-local index type
+ using LI = std::uint8_t;
+
+ // the specified maximum admissible number of branches per scvf
+ static constexpr int maxB = GET_PROP_VALUE(TypeTag, MaxNumBranchesPerScvf);
+
+ // maximum admissible number of elements around a node
+ // if for a given grid this number is still not high enough,
+ // overwrite this property in your problem with a higher number
+ static constexpr int dim = GridView::dimension;
+ static constexpr int maxE = dim == 3 ? 45 : 15;
+
+public:
+ using type = CCMpfaDualGridNodalIndexSet<GI, LI, dim, maxE, maxB>;
+};
+
//! The mpfa helper class
SET_TYPE_PROP(CCMpfaModel, MpfaHelper, CCMpfaHelper<TypeTag>);
diff --git a/dumux/discretization/cellcentered/tpfa/properties.hh b/dumux/discretization/cellcentered/tpfa/properties.hh
index 4e77be6a2f5e28d9573d0f3d902b364b464acaab..008770527745e92c7791bd11ac82717d2c6ec539 100644
--- a/dumux/discretization/cellcentered/tpfa/properties.hh
+++ b/dumux/discretization/cellcentered/tpfa/properties.hh
@@ -91,7 +91,7 @@ private:
static constexpr int dimWorld = GridView::dimensionworld;
public:
- //! Per default, we allow for 8 neighbors on network/surface grids
+ // Per default, we allow for 8 neighbors on network/surface grids
static const std::size_t value = dim < dimWorld ? 9 : 2;
};
diff --git a/dumux/discretization/fluxstencil.hh b/dumux/discretization/fluxstencil.hh
index 6ee56e89d2050398bc4905ca2fd30ed37677e16e..6c7f329ddc19970d15fa0f37557ab3101011f239 100644
--- a/dumux/discretization/fluxstencil.hh
+++ b/dumux/discretization/fluxstencil.hh
@@ -100,13 +100,36 @@ class FluxStencilImplementation<TypeTag, DiscretizationMethods::CCMpfa>
using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
using Element = typename GridView::template Codim<0>::Entity;
using IndexType = typename GridView::IndexSet::IndexType;
+ static constexpr int dim = GridView::dimension;
+
+ // Use the stencil type of the primary interaction volume
+ using NodalIndexSet = typename GET_PROP_TYPE(TypeTag, DualGridNodalIndexSet);
public:
- using Stencil = std::vector<IndexType>;
+ //! The maximum number of elements in a flux stencil (equal to number of elements at node)
+ static constexpr int maxFluxStencilSize = NodalIndexSet::maxNumElementsAtNode;
- static Stencil stencil(const Element& element,
- const FVElementGeometry& fvGeometry,
- const SubControlVolumeFace& scvf)
+ //! States how many scvfs of an element J might have an element I in the flux stencil
+ //! We use cubes here for determining the maximum. Only basic geometries are not supported.
+ static constexpr int maxNumScvfJForI = dim == 3 ? 12 : 4;
+
+ //! The flux stencil type
+ using Stencil = typename NodalIndexSet::GridStencilType;
+
+ /*
+ * \brief Returns a set of grid element indices that participate in the
+ * flux calculations on a given scvf.
+ *
+ * \note The interaction volume index sets must use the same type for the
+ * stencils as the nodal index set. If not, the compiler will complain here.
+ *
+ * \param element The grid element
+ * \param fvGeometry The finite volume geometry of this element
+ * \param scvf The sub-control volume face embedded in this element
+ */
+ static const Stencil& stencil(const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const SubControlVolumeFace& scvf)
{
const auto& fvGridGeometry = fvGeometry.fvGridGeometry();