Commit 1fc1eb5e authored by Dennis Gläser's avatar Dennis Gläser
Browse files

[mpfa] use reserved vectors where possible

parent 2f8664f2
......@@ -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:
......
......@@ -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 >
......
......@@ -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_;
......
......@@ -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;
......
......@@ -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:
/*!
......
......@@ -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(); }
};
......
......@@ -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." );
......
......@@ -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
......
......@@ -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
......
......@@ -52,7 +52,8 @@ template< class TypeTag, int localSize >
struct CCMpfaODefaultStaticInteractionVolumeTraits
{
private: