From 0b7b974cc7894b22299e2796adef591125160259 Mon Sep 17 00:00:00 2001
From: DennisGlaeser <dennis.glaeser@iws.uni-stuttgart.de>
Date: Tue, 16 Jan 2018 15:21:06 +0100
Subject: [PATCH] [cc][scvf] allow for customizable type for storing neighbor
indices
This commit furthermore introduces default traits classes to be used for
the scvs and scvfs. Before, these were locally defined in the body of
the SET_PROP construction in the properties.
---
.../cellcentered/mpfa/dualgridindexset.hh | 3 +-
.../cellcentered/mpfa/fvelementgeometry.hh | 2 +-
.../cellcentered/mpfa/fvgridgeometry.hh | 50 ++++++-----------
.../cellcentered/mpfa/properties.hh | 55 ++-----------------
.../cellcentered/mpfa/subcontrolvolumeface.hh | 52 +++++++++++++++++-
.../cellcentered/subcontrolvolume.hh | 17 ++++++
.../cellcentered/tpfa/fvelementgeometry.hh | 29 +++++++---
.../cellcentered/tpfa/fvgridgeometry.hh | 44 ++++++++-------
.../cellcentered/tpfa/properties.hh | 44 ++-------------
.../cellcentered/tpfa/subcontrolvolumeface.hh | 44 ++++++++++++++-
10 files changed, 187 insertions(+), 153 deletions(-)
diff --git a/dumux/discretization/cellcentered/mpfa/dualgridindexset.hh b/dumux/discretization/cellcentered/mpfa/dualgridindexset.hh
index a93c7ce639..2288f3b49d 100644
--- a/dumux/discretization/cellcentered/mpfa/dualgridindexset.hh
+++ b/dumux/discretization/cellcentered/mpfa/dualgridindexset.hh
@@ -87,10 +87,11 @@ public:
}
//! Inserts scvf data
+ template<typename OutsideGridIndexStorage>
void insert(const bool boundary,
const GridIndexType scvfIdx,
const GridIndexType insideScvIdx,
- const std::vector<GridIndexType>& outsideScvIndices)
+ const OutsideGridIndexStorage& outsideScvIndices)
{
// this should always be called only once per scvf
assert(std::find(scvfIndices_.begin(), scvfIndices_.end(), scvfIdx ) == scvfIndices_.end()
diff --git a/dumux/discretization/cellcentered/mpfa/fvelementgeometry.hh b/dumux/discretization/cellcentered/mpfa/fvelementgeometry.hh
index 44b8f899fe..66928c9121 100644
--- a/dumux/discretization/cellcentered/mpfa/fvelementgeometry.hh
+++ b/dumux/discretization/cellcentered/mpfa/fvelementgeometry.hh
@@ -378,7 +378,7 @@ private:
if (dim < dimWorld)
{
const auto indexInInside = is.indexInInside();
- if(finishedFacets[indexInInside])
+ if (finishedFacets[indexInInside])
continue;
else
finishedFacets[indexInInside] = true;
diff --git a/dumux/discretization/cellcentered/mpfa/fvgridgeometry.hh b/dumux/discretization/cellcentered/mpfa/fvgridgeometry.hh
index 248c0129d9..c355667c6c 100644
--- a/dumux/discretization/cellcentered/mpfa/fvgridgeometry.hh
+++ b/dumux/discretization/cellcentered/mpfa/fvgridgeometry.hh
@@ -69,6 +69,7 @@ class CCMpfaFVGridGeometry<TypeTag, true> : public BaseFVGridGeometry<TypeTag>
using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
+ using ScvfOutsideGridIndexStorage = typename SubControlVolumeFace::Traits::OutsideGridIndexStorage;
using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
@@ -185,9 +186,9 @@ public:
std::vector<GridIndexType> scvfIndexSet;
scvfIndexSet.reserve(MpfaHelper::getNumLocalScvfs(elementGeometry.type()));
- // for network grids there might be multiple intersections with the same geometryInInside
+ // for network grids there might be multiple intersection with the same geometryInInside
// we indentify those by the indexInInside for now (assumes conforming grids at branching facets)
- std::vector<std::vector<GridIndexType>> outsideIndices;
+ std::vector<ScvfOutsideGridIndexStorage> outsideIndices;
if (dim < dimWorld)
{
outsideIndices.resize(element.subEntities(1));
@@ -195,8 +196,7 @@ public:
{
if (intersection.neighbor())
{
- const auto& outside = intersection.outside();
- auto nIdx = this->elementMapper().index(outside);
+ const auto nIdx = this->elementMapper().index( intersection.outside() );
outsideIndices[intersection.indexInInside()].push_back(nIdx);
}
}
@@ -257,18 +257,11 @@ public:
const auto& outsideScvIndices = [&] ()
{
if (!boundary)
- {
return dim == dimWorld ?
- std::vector<GridIndexType>({this->elementMapper().index(is.outside())}) :
+ ScvfOutsideGridIndexStorage({this->elementMapper().index(is.outside())}) :
outsideIndices[indexInInside];
- }
else
- {
- // compute boundary scv idx and increment counter
- const GridIndexType bIdx = numScvs + numBoundaryScvf_;
- numBoundaryScvf_++;
- return std::vector<GridIndexType>(1, bIdx);
- }
+ return ScvfOutsideGridIndexStorage({GridIndexType(numScvs) + numBoundaryScvf_++});
} ();
scvfIndexSet.push_back(scvfIdx);
@@ -382,7 +375,7 @@ private:
// containers storing the global data
std::vector<std::vector<GridIndexType>> scvfIndicesOfScv_;
std::vector<bool> secondaryInteractionVolumeVertices_;
- std::size_t numBoundaryScvf_;
+ GridIndexType numBoundaryScvf_;
// needed for embedded surface and network grids (dim < dimWorld)
std::vector<std::vector<GridIndexType>> flipScvfIndices_;
@@ -414,6 +407,7 @@ class CCMpfaFVGridGeometry<TypeTag, false> : public BaseFVGridGeometry<TypeTag>
using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
+ using ScvfOutsideGridIndexStorage = typename SubControlVolumeFace::Traits::OutsideGridIndexStorage;
using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
@@ -529,13 +523,13 @@ public:
// the element-wise index sets for finite volume geometry
const auto numLocalFaces = MpfaHelper::getNumLocalScvfs(elementGeometry.type());
std::vector<GridIndexType> scvfsIndexSet;
- std::vector< std::vector<GridIndexType> > neighborVolVarIndexSet;
+ std::vector<ScvfOutsideGridIndexStorage> neighborVolVarIndexSet;
scvfsIndexSet.reserve(numLocalFaces);
neighborVolVarIndexSet.reserve(numLocalFaces);
// for network grids there might be multiple intersections with the same geometryInInside
// we indentify those by the indexInInside for now (assumes conforming grids at branching facets)
- std::vector<std::vector<GridIndexType>> outsideIndices;
+ std::vector<ScvfOutsideGridIndexStorage> outsideIndices;
if (dim < dimWorld)
{
outsideIndices.resize(element.subEntities(1));
@@ -543,8 +537,7 @@ public:
{
if (intersection.neighbor())
{
- const auto& outside = intersection.outside();
- auto nIdx = this->elementMapper().index(outside);
+ auto nIdx = this->elementMapper().index( intersection.outside() );
outsideIndices[intersection.indexInInside()].push_back(nIdx);
}
}
@@ -595,18 +588,11 @@ public:
const auto& outsideScvIndices = [&] ()
{
if (!boundary)
- {
return dim == dimWorld ?
- std::vector<GridIndexType>({this->elementMapper().index(is.outside())}) :
+ ScvfOutsideGridIndexStorage({this->elementMapper().index(is.outside())}) :
outsideIndices[indexInInside];
- }
else
- {
- // compute boundary scv idx and increment counter
- const GridIndexType bIdx = numScvs_ + numBoundaryScvf_;
- numBoundaryScvf_++;
- return std::vector<GridIndexType>(1, bIdx);
- }
+ return ScvfOutsideGridIndexStorage({GridIndexType(numScvs_) + numBoundaryScvf_++});
} ();
// insert the scvf data into the dual grid index set
@@ -646,7 +632,7 @@ public:
{ return scvfIndicesOfScv_[scvIdx]; }
//! Returns the neighboring vol var indices for each scvf contained in an scv.
- const std::vector< std::vector<GridIndexType> >& neighborVolVarIndices(GridIndexType scvIdx) const
+ const std::vector<ScvfOutsideGridIndexStorage>& neighborVolVarIndices(GridIndexType scvIdx) const
{ return neighborVolVarIndices_[scvIdx]; }
//! Returns the connectivity map of which dofs
@@ -662,12 +648,12 @@ private:
// containers storing the global data
std::vector<std::vector<GridIndexType>> scvfIndicesOfScv_;
- std::vector< std::vector< std::vector<GridIndexType> > > neighborVolVarIndices_;
+ std::vector<std::vector<ScvfOutsideGridIndexStorage>> neighborVolVarIndices_;
std::vector<bool> secondaryInteractionVolumeVertices_;
std::vector<bool> isGhostVertex_;
- std::size_t numScvs_;
- std::size_t numScvf_;
- std::size_t numBoundaryScvf_;
+ GridIndexType numScvs_;
+ GridIndexType numScvf_;
+ GridIndexType numBoundaryScvf_;
// The grid interaction volume index set
GridIVIndexSets ivIndexSets_;
diff --git a/dumux/discretization/cellcentered/mpfa/properties.hh b/dumux/discretization/cellcentered/mpfa/properties.hh
index ce805e818a..aeb2b83479 100644
--- a/dumux/discretization/cellcentered/mpfa/properties.hh
+++ b/dumux/discretization/cellcentered/mpfa/properties.hh
@@ -25,9 +25,6 @@
#ifndef DUMUX_CC_MPFA_PROPERTIES_HH
#define DUMUX_CC_MPFA_PROPERTIES_HH
-#include <dune/common/fvector.hh>
-#include <dune/geometry/multilineargeometry.hh>
-
#include <dumux/common/properties.hh>
#include <dumux/assembly/cclocalresidual.hh>
@@ -165,60 +162,20 @@ SET_TYPE_PROP(CCMpfaModel, GridVolumeVariables, CCGridVolumeVariables<TypeTag, G
SET_PROP(CCMpfaModel, SubControlVolume)
{
private:
- using Grid = typename GET_PROP_TYPE(TypeTag, Grid);
- struct ScvGeometryTraits
- {
- using Geometry = typename Grid::template Codim<0>::Geometry;
- using GridIndexType = typename Grid::LeafGridView::IndexSet::IndexType;
- using LocalIndexType = unsigned int;
- using Scalar = typename Geometry::ctype;
- using GlobalPosition = Dune::FieldVector<Scalar, Geometry::coorddimension>;
- };
+ using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+ using Traits = CCDefaultScvGeometryTraits<GridView>;
public:
- using type = CCSubControlVolume<ScvGeometryTraits>;
+ using type = CCSubControlVolume<Traits>;
};
//! The sub-control volume face class
SET_PROP(CCMpfaModel, SubControlVolumeFace)
{
private:
- using Grid = typename GET_PROP_TYPE(TypeTag, Grid);
- static const int dim = Grid::dimension;
- static const int dimWorld = Grid::dimensionworld;
-
- // we use geometry traits that use static corner vectors to and a fixed geometry type
- template <class ct>
- struct ScvfMLGTraits : public Dune::MultiLinearGeometryTraits<ct>
- {
- // we use static vectors to store the corners as we know
- // the number of corners in advance (2^(dim-1) corners (1<<(dim-1))
- template< int mydim, int cdim >
- struct CornerStorage
- {
- using Type = std::array< Dune::FieldVector< ct, cdim >, (1<<(dim-1)) >;
- };
-
- // we know all scvfs will have the same geometry type
- template< int dim >
- struct hasSingleGeometryType
- {
- static const bool v = true;
- static const unsigned int topologyId = Dune::Impl::CubeTopology< dim >::type::id;
- };
- };
-
- struct ScvfGeometryTraits
- {
- using GridIndexType = typename Grid::LeafGridView::IndexSet::IndexType;
- using LocalIndexType = unsigned int;
- using Scalar = typename Grid::ctype;
- using Geometry = Dune::MultiLinearGeometry<Scalar, dim-1, dimWorld, ScvfMLGTraits<Scalar> >;
- using CornerStorage = typename ScvfMLGTraits<Scalar>::template CornerStorage<dim-1, dimWorld>::Type;
- using GlobalPosition = typename CornerStorage::value_type;
- };
-
+ using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+ using Traits = CCMpfaDefaultScvfGeometryTraits<GridView>;
public:
- using type = Dumux::CCMpfaSubControlVolumeFace< ScvfGeometryTraits>;
+ using type = Dumux::CCMpfaSubControlVolumeFace<Traits>;
};
//! Set the solution vector type for an element
diff --git a/dumux/discretization/cellcentered/mpfa/subcontrolvolumeface.hh b/dumux/discretization/cellcentered/mpfa/subcontrolvolumeface.hh
index 45bd3cc55a..dfc14bfc61 100644
--- a/dumux/discretization/cellcentered/mpfa/subcontrolvolumeface.hh
+++ b/dumux/discretization/cellcentered/mpfa/subcontrolvolumeface.hh
@@ -25,12 +25,57 @@
#define DUMUX_DISCRETIZATION_CC_MPFA_SUBCONTROLVOLUMEFACE_HH
#include <vector>
+#include <array>
#include <dune/common/version.hh>
+#include <dune/common/reservedvector.hh>
+#include <dune/common/fvector.hh>
#include <dune/geometry/type.hh>
+#include <dune/geometry/multilineargeometry.hh>
namespace Dumux
{
+template<class GridView>
+struct CCMpfaDefaultScvfGeometryTraits
+{
+ using Grid = typename GridView::Grid;
+
+ static const int dim = Grid::dimension;
+ static const int dimWorld = Grid::dimensionworld;
+
+ using Scalar = typename Grid::ctype;
+ using GridIndexType = typename Grid::LeafGridView::IndexSet::IndexType;
+ using LocalIndexType = unsigned int;
+ using OutsideGridIndexStorage = typename std::conditional_t< (dim<dimWorld),
+ std::vector<GridIndexType>,
+ Dune::ReservedVector<GridIndexType, 1> >;
+
+ // we use geometry traits that use static corner vectors to and a fixed geometry type
+ template <class ct>
+ struct ScvfMLGTraits : public Dune::MultiLinearGeometryTraits<ct>
+ {
+ // we use static vectors to store the corners as we know
+ // the number of corners in advance (2^(dim-1) corners (1<<(dim-1))
+ template< int mydim, int cdim >
+ struct CornerStorage
+ {
+ using Type = std::array< Dune::FieldVector< ct, cdim >, (1<<(dim-1)) >;
+ };
+
+ // we know all scvfs will have the same geometry type
+ template< int dim >
+ struct hasSingleGeometryType
+ {
+ static const bool v = true;
+ static const unsigned int topologyId = Dune::Impl::CubeTopology< dim >::type::id;
+ };
+ };
+
+ using Geometry = Dune::MultiLinearGeometry<Scalar, dim-1, dimWorld, ScvfMLGTraits<Scalar> >;
+ using CornerStorage = typename ScvfMLGTraits<Scalar>::template CornerStorage<dim-1, dimWorld>::Type;
+ using GlobalPosition = typename CornerStorage::value_type;
+};
+
/*!
* \ingroup CCMpfaDiscretization
* \brief Class for a sub control volume face in mpfa methods, i.e a part of the boundary
@@ -45,6 +90,7 @@ class CCMpfaSubControlVolumeFace
using Scalar = typename ScvfGeometryTraits::Scalar;
using GlobalPosition = typename ScvfGeometryTraits::GlobalPosition;
using CornerStorage = typename ScvfGeometryTraits::CornerStorage;
+ using OutsideGridIndexStorage = typename ScvfGeometryTraits::OutsideGridIndexStorage;
using Geometry = typename ScvfGeometryTraits::Geometry;
public:
@@ -73,7 +119,7 @@ public:
unsigned int vIdxLocal,
GridIndexType scvfIndex,
GridIndexType insideScvIdx,
- const std::vector<GridIndexType>& outsideScvIndices,
+ const OutsideGridIndexStorage& outsideScvIndices,
Scalar q,
bool boundary)
: boundary_(boundary)
@@ -122,7 +168,7 @@ public:
GridIndexType outsideScvIdx(int i = 0) const { return outsideScvIndices_[i]; }
//! returns the outside scv indices (can be more than one index for dim < dimWorld)
- const std::vector<GridIndexType>& outsideScvIndices() const { return outsideScvIndices_; }
+ const OutsideGridIndexStorage& outsideScvIndices() const { return outsideScvIndices_; }
//! Returns the number of corners
std::size_t corners() const { return corners_.size(); }
@@ -161,7 +207,7 @@ private:
GridIndexType vertexIndex_;
GridIndexType scvfIndex_;
GridIndexType insideScvIdx_;
- std::vector<GridIndexType> outsideScvIndices_;
+ OutsideGridIndexStorage outsideScvIndices_;
unsigned int vIdxInElement_;
CornerStorage corners_;
diff --git a/dumux/discretization/cellcentered/subcontrolvolume.hh b/dumux/discretization/cellcentered/subcontrolvolume.hh
index c9521b8bbe..1c6cd9db19 100644
--- a/dumux/discretization/cellcentered/subcontrolvolume.hh
+++ b/dumux/discretization/cellcentered/subcontrolvolume.hh
@@ -24,11 +24,28 @@
#ifndef DUMUX_DISCRETIZATION_CC_SUBCONTROLVOLUME_HH
#define DUMUX_DISCRETIZATION_CC_SUBCONTROLVOLUME_HH
+#include <dune/common/fvector.hh>
#include <dumux/discretization/subcontrolvolumebase.hh>
#include <dumux/common/optional.hh>
namespace Dumux
{
+
+/*!
+ * \ingroup CCDiscretization
+ * \brief Default traits class to be used for the sub-control volumes
+ * for the cell-centered finite volume scheme using TPFA
+ */
+template<class GridView>
+struct CCDefaultScvGeometryTraits
+{
+ using Geometry = typename GridView::template Codim<0>::Geometry;
+ using GridIndexType = typename GridView::IndexSet::IndexType;
+ using LocalIndexType = unsigned int;
+ using Scalar = typename GridView::ctype;
+ using GlobalPosition = Dune::FieldVector<Scalar, GridView::dimensionworld>;
+};
+
/*!
* \ingroup CCDiscretization
* \brief Sub control volumes for cell-centered discretization schemes
diff --git a/dumux/discretization/cellcentered/tpfa/fvelementgeometry.hh b/dumux/discretization/cellcentered/tpfa/fvelementgeometry.hh
index b77ff57df6..4ff6ebcdd9 100644
--- a/dumux/discretization/cellcentered/tpfa/fvelementgeometry.hh
+++ b/dumux/discretization/cellcentered/tpfa/fvelementgeometry.hh
@@ -26,6 +26,8 @@
#ifndef DUMUX_DISCRETIZATION_CCTPFA_FV_ELEMENT_GEOMETRY_HH
#define DUMUX_DISCRETIZATION_CCTPFA_FV_ELEMENT_GEOMETRY_HH
+#include <algorithm>
+
#include <dune/common/exceptions.hh>
#include <dune/common/iteratorrange.hh>
#include <dumux/common/properties.hh>
@@ -395,6 +397,8 @@ private:
//! create scvs and scvfs of the bound element
void makeElementGeometries(const Element& element)
{
+ using ScvfGridIndexStorage = typename SubControlVolumeFace::Traits::GridIndexStorage;
+
const auto eIdx = fvGridGeometry().elementMapper().index(element);
scvs_[0] = SubControlVolume(element.geometry(), eIdx);
scvIndices_[0] = eIdx;
@@ -416,10 +420,14 @@ private:
if (handledScvf[intersection.indexInInside()])
continue;
+ const auto& scvfNeighborVolVarIndices = neighborVolVarIndices[scvfCounter];
+
if (intersection.neighbor() || intersection.boundary())
{
- std::vector<IndexType> scvIndices({eIdx});
- scvIndices.insert(scvIndices.end(), neighborVolVarIndices[scvfCounter].begin(), neighborVolVarIndices[scvfCounter].end());
+ ScvfGridIndexStorage scvIndices;
+ scvIndices.resize(scvfNeighborVolVarIndices.size() + 1);
+ scvIndices[0] = eIdx;
+ std::copy(scvfNeighborVolVarIndices.begin(), scvfNeighborVolVarIndices.end(), scvIndices.begin()+1);
scvfs_.emplace_back(intersection,
intersection.geometry(),
scvFaceIndices[scvfCounter],
@@ -438,6 +446,8 @@ private:
//! 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);
@@ -459,6 +469,8 @@ private:
if (handledScvf[intersection.indexInInside()])
continue;
+ const auto& scvfNeighborVolVarIndices = neighborVolVarIndices[scvfCounter];
+
if (intersection.neighbor())
{
// only create subcontrol faces where the outside element is the bound element
@@ -466,8 +478,7 @@ private:
{
if (intersection.outside() == *elementPtr_)
{
- std::vector<IndexType> scvIndices({eIdx});
- scvIndices.insert(scvIndices.end(), neighborVolVarIndices[scvfCounter].begin(), neighborVolVarIndices[scvfCounter].end());
+ ScvfGridIndexStorage scvIndices({eIdx, scvfNeighborVolVarIndices[0]});
neighborScvfs_.emplace_back(intersection,
intersection.geometry(),
scvFaceIndices[scvfCounter],
@@ -483,12 +494,14 @@ private:
// (will be optimized away for dim == dimWorld)
else
{
- for (unsigned outsideScvIdx = 0; outsideScvIdx < neighborVolVarIndices[scvfCounter].size(); ++outsideScvIdx)
+ for (unsigned outsideScvIdx = 0; outsideScvIdx < scvfNeighborVolVarIndices.size(); ++outsideScvIdx)
{
- if (neighborVolVarIndices[scvfCounter][outsideScvIdx] == fvGridGeometry().elementMapper().index(*elementPtr_))
+ if (scvfNeighborVolVarIndices[outsideScvIdx] == fvGridGeometry().elementMapper().index(*elementPtr_))
{
- std::vector<IndexType> scvIndices({eIdx});
- scvIndices.insert(scvIndices.end(), neighborVolVarIndices[scvfCounter].begin(), neighborVolVarIndices[scvfCounter].end());
+ ScvfGridIndexStorage scvIndices;
+ scvIndices.resize(scvfNeighborVolVarIndices.size() + 1);
+ scvIndices[0] = eIdx;
+ std::copy(scvfNeighborVolVarIndices.begin(), scvfNeighborVolVarIndices.end(), scvIndices.begin()+1);
neighborScvfs_.emplace_back(intersection,
intersection.geometry(),
scvFaceIndices[scvfCounter],
diff --git a/dumux/discretization/cellcentered/tpfa/fvgridgeometry.hh b/dumux/discretization/cellcentered/tpfa/fvgridgeometry.hh
index 520823347a..08a559b12a 100644
--- a/dumux/discretization/cellcentered/tpfa/fvgridgeometry.hh
+++ b/dumux/discretization/cellcentered/tpfa/fvgridgeometry.hh
@@ -26,6 +26,8 @@
#ifndef DUMUX_DISCRETIZATION_CCTPFA_FV_GRID_GEOMETRY_HH
#define DUMUX_DISCRETIZATION_CCTPFA_FV_GRID_GEOMETRY_HH
+#include <algorithm>
+
#include <dune/common/version.hh>
#include <dumux/discretization/methods.hh>
@@ -60,6 +62,7 @@ class CCTpfaFVGridGeometry<TypeTag, true> : public BaseFVGridGeometry<TypeTag>
using IndexType = typename GridView::IndexSet::IndexType;
using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
+ using ScvfGridIndexStorage = typename SubControlVolumeFace::Traits::GridIndexStorage;
using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
using ElementMapper = typename GET_PROP_TYPE(TypeTag, ElementMapper);
using Element = typename GridView::template Codim<0>::Entity;
@@ -154,19 +157,21 @@ public:
// for network grids there might be multiple intersection with the same geometryInInside
// we indentify those by the indexInInside for now (assumes conforming grids at branching facets)
- std::vector<std::vector<IndexType>> outsideIndices;
+ std::vector<ScvfGridIndexStorage> outsideIndices;
if (dim < dimWorld)
{
+ //! first, push inside index in all neighbor sets
outsideIndices.resize(element.subEntities(1));
+ std::for_each(outsideIndices.begin(), outsideIndices.end(), [eIdx] (auto& nIndices) { nIndices.push_back(eIdx); });
+
+ // second, insert neighbors
for (const auto& intersection : intersections(this->gridView(), element))
{
if (intersection.neighbor())
{
- const auto& outside = intersection.outside();
- const auto nIdx = this->elementMapper().index(outside);
+ const auto nIdx = this->elementMapper().index( intersection.outside() );
outsideIndices[intersection.indexInInside()].push_back(nIdx);
}
-
}
}
@@ -181,7 +186,7 @@ public:
scvfs_.emplace_back(intersection,
intersection.geometry(),
scvfIdx,
- std::vector<IndexType>({eIdx, nIdx}),
+ ScvfGridIndexStorage({eIdx, nIdx}),
false);
scvfsIndexSet.push_back(scvfIdx++);
}
@@ -195,12 +200,10 @@ public:
continue;
else
{
- std::vector<IndexType> scvIndices({eIdx});
- scvIndices.insert(scvIndices.end(), outsideIndices[indexInInside].begin(), outsideIndices[indexInInside].end());
scvfs_.emplace_back(intersection,
intersection.geometry(),
scvfIdx,
- scvIndices,
+ outsideIndices[indexInInside],
false);
scvfsIndexSet.push_back(scvfIdx++);
outsideIndices[indexInInside].clear();
@@ -213,7 +216,7 @@ public:
scvfs_.emplace_back(intersection,
intersection.geometry(),
scvfIdx,
- std::vector<IndexType>({eIdx, this->gridView().size(0) + numBoundaryScvf_++}),
+ ScvfGridIndexStorage({eIdx, this->gridView().size(0) + numBoundaryScvf_++}),
true);
scvfsIndexSet.push_back(scvfIdx++);
}
@@ -320,6 +323,7 @@ class CCTpfaFVGridGeometry<TypeTag, false> : public BaseFVGridGeometry<TypeTag>
using IndexType = typename GridView::IndexSet::IndexType;
using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
+ using ScvfGridIndexStorage = typename SubControlVolumeFace::Traits::GridIndexStorage;
using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
using ElementMapper = typename GET_PROP_TYPE(TypeTag, ElementMapper);
using Element = typename GridView::template Codim<0>::Entity;
@@ -331,6 +335,10 @@ class CCTpfaFVGridGeometry<TypeTag, false> : public BaseFVGridGeometry<TypeTag>
using CoordScalar = typename GridView::ctype;
using GlobalPosition = Dune::FieldVector<CoordScalar, dimWorld>;
+ using NeighborVolVarIndices = typename std::conditional_t< (dim<dimWorld),
+ ScvfGridIndexStorage,
+ Dune::ReservedVector<IndexType, 1> >;
+
public:
//! export discretization method
static constexpr DiscretizationMethods discretizationMethod = DiscretizationMethods::CCTpfa;
@@ -399,13 +407,13 @@ public:
// the element-wise index sets for finite volume geometry
auto numLocalFaces = element.subEntities(1);
std::vector<IndexType> scvfsIndexSet;
- std::vector<std::vector<IndexType>> neighborVolVarIndexSet;
+ std::vector<NeighborVolVarIndices> neighborVolVarIndexSet;
scvfsIndexSet.reserve(numLocalFaces);
neighborVolVarIndexSet.reserve(numLocalFaces);
// for network grids there might be multiple intersection with the same geometryInInside
// we indentify those by the indexInInside for now (assumes conforming grids at branching facets)
- std::vector<std::vector<IndexType>> outsideIndices;
+ std::vector<NeighborVolVarIndices> outsideIndices;
if (dim < dimWorld)
{
outsideIndices.resize(numLocalFaces);
@@ -413,11 +421,9 @@ public:
{
if (intersection.neighbor())
{
- const auto& outside = intersection.outside();
- const auto nIdx = this->elementMapper().index(outside);
+ const auto nIdx = this->elementMapper().index(intersection.outside());
outsideIndices[intersection.indexInInside()].push_back(nIdx);
}
-
}
}
@@ -430,7 +436,7 @@ public:
{
scvfsIndexSet.push_back(numScvf_++);
const auto nIdx = this->elementMapper().index(intersection.outside());
- neighborVolVarIndexSet.push_back({nIdx});
+ neighborVolVarIndexSet.emplace_back(NeighborVolVarIndices({nIdx}));
}
// this is for network grids
// (will be optimized away of dim == dimWorld)
@@ -443,7 +449,7 @@ public:
else
{
scvfsIndexSet.push_back(numScvf_++);
- neighborVolVarIndexSet.push_back(outsideIndices[indexInInside]);
+ neighborVolVarIndexSet.emplace_back(std::move(outsideIndices[indexInInside]));
outsideIndices[indexInInside].clear();
}
}
@@ -452,7 +458,7 @@ public:
else if (intersection.boundary())
{
scvfsIndexSet.push_back(numScvf_++);
- neighborVolVarIndexSet.push_back({numScvs_ + numBoundaryScvf_++});
+ neighborVolVarIndexSet.emplace_back(NeighborVolVarIndices({numScvs_ + numBoundaryScvf_++}));
}
}
@@ -469,7 +475,7 @@ public:
{ return scvfIndicesOfScv_[scvIdx]; }
//! Return the neighbor volVar indices for all scvfs in the scv with index scvIdx
- const std::vector<std::vector<IndexType>>& neighborVolVarIndices(IndexType scvIdx) const
+ const std::vector<NeighborVolVarIndices>& neighborVolVarIndices(IndexType scvIdx) const
{ return neighborVolVarIndices_[scvIdx]; }
/*!
@@ -491,7 +497,7 @@ private:
//! vectors that store the global data
std::vector<std::vector<IndexType>> scvfIndicesOfScv_;
- std::vector<std::vector<std::vector<IndexType>>> neighborVolVarIndices_;
+ std::vector<std::vector<NeighborVolVarIndices>> neighborVolVarIndices_;
};
} // end namespace Dumux
diff --git a/dumux/discretization/cellcentered/tpfa/properties.hh b/dumux/discretization/cellcentered/tpfa/properties.hh
index c9262c4b14..1524324e28 100644
--- a/dumux/discretization/cellcentered/tpfa/properties.hh
+++ b/dumux/discretization/cellcentered/tpfa/properties.hh
@@ -99,52 +99,20 @@ public:
SET_PROP(CCTpfaModel, SubControlVolume)
{
private:
- using Grid = typename GET_PROP_TYPE(TypeTag, Grid);
- struct ScvGeometryTraits
- {
- using Geometry = typename Grid::template Codim<0>::Geometry;
- using GridIndexType = typename Grid::LeafGridView::IndexSet::IndexType;
- using LocalIndexType = unsigned int;
- using Scalar = typename Grid::ctype;
- using GlobalPosition = Dune::FieldVector<Scalar, Grid::dimensionworld>;
- };
+ using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+ using Traits = CCDefaultScvGeometryTraits<GridView>;
public:
- using type = CCSubControlVolume<ScvGeometryTraits>;
+ using type = CCSubControlVolume<Traits>;
};
//! The sub control volume face
SET_PROP(CCTpfaModel, SubControlVolumeFace)
{
private:
- using Grid = typename GET_PROP_TYPE(TypeTag, Grid);
- static constexpr int dim = Grid::dimension;
- static constexpr int dimWorld = Grid::dimensionworld;
-
- // we use geometry traits that use static corner vectors to and a fixed geometry type
- template <class ct>
- struct ScvfMLGTraits : public Dune::MultiLinearGeometryTraits<ct>
- {
- // we use static vectors to store the corners as we know
- // the number of corners in advance (2^(dim-1) corners (1<<(dim-1))
- template< int mydim, int cdim >
- struct CornerStorage
- {
- using Type = Dune::ReservedVector< Dune::FieldVector< ct, cdim >, (1<<(dim-1)) >;
- };
- };
-
- struct ScvfGeometryTraits
- {
- using GridIndexType = typename Grid::LeafGridView::IndexSet::IndexType;
- using LocalIndexType = unsigned int;
- using Scalar = typename Grid::ctype;
- using Geometry = Dune::MultiLinearGeometry<Scalar, dim-1, dimWorld, ScvfMLGTraits<Scalar> >;
- using CornerStorage = typename ScvfMLGTraits<Scalar>::template CornerStorage<dim-1, dimWorld>::Type;
- using GlobalPosition = typename CornerStorage::value_type;
- using BoundaryFlag = Dumux::BoundaryFlag<Grid>;
- };
+ using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+ using Traits = CCTpfaDefaultScvfGeometryTraits<GridView>;
public:
- using type = Dumux::CCTpfaSubControlVolumeFace<ScvfGeometryTraits>;
+ using type = Dumux::CCTpfaSubControlVolumeFace<Traits>;
};
//! Set the solution vector type for an element
diff --git a/dumux/discretization/cellcentered/tpfa/subcontrolvolumeface.hh b/dumux/discretization/cellcentered/tpfa/subcontrolvolumeface.hh
index 58c837af98..755dbeef0d 100644
--- a/dumux/discretization/cellcentered/tpfa/subcontrolvolumeface.hh
+++ b/dumux/discretization/cellcentered/tpfa/subcontrolvolumeface.hh
@@ -31,6 +31,45 @@
namespace Dumux
{
+/*!
+ * \ingroup CCDiscretization
+ * \brief Default traits class to be used for the sub-control volume faces
+ * for the cell-centered finite volume scheme using TPFA
+ */
+template<class GridView>
+struct CCTpfaDefaultScvfGeometryTraits
+{
+ using Grid = typename GridView::Grid;
+
+ static constexpr int dim = Grid::dimension;
+ static constexpr int dimWorld = Grid::dimensionworld;
+
+ using Scalar = typename Grid::ctype;
+ using GridIndexType = typename Grid::LeafGridView::IndexSet::IndexType;
+ using LocalIndexType = unsigned int;
+ using GridIndexStorage = typename std::conditional_t< (dim<dimWorld),
+ std::vector<GridIndexType>,
+ Dune::ReservedVector<GridIndexType, 2> >;
+
+ // we use geometry traits that use static corner vectors to and a fixed geometry type
+ template <class ct>
+ struct ScvfMLGTraits : public Dune::MultiLinearGeometryTraits<ct>
+ {
+ // we use static vectors to store the corners as we know
+ // the number of corners in advance (2^(dim-1) corners (1<<(dim-1))
+ template< int mydim, int cdim >
+ struct CornerStorage
+ {
+ using Type = Dune::ReservedVector< Dune::FieldVector< ct, cdim >, (1<<(dim-1)) >;
+ };
+ };
+
+ using Geometry = Dune::MultiLinearGeometry<Scalar, dim-1, dimWorld, ScvfMLGTraits<Scalar> >;
+ using CornerStorage = typename ScvfMLGTraits<Scalar>::template CornerStorage<dim-1, dimWorld>::Type;
+ using GlobalPosition = typename CornerStorage::value_type;
+ using BoundaryFlag = Dumux::BoundaryFlag<Grid>;
+};
+
/*!
* \ingroup CCTpfaDiscretization
* \brief The sub control volume face
@@ -43,6 +82,7 @@ class CCTpfaSubControlVolumeFace : public SubControlVolumeFaceBase<CCTpfaSubCont
using Scalar = typename ScvfGeometryTraits::Scalar;
using GlobalPosition = typename ScvfGeometryTraits::GlobalPosition;
using CornerStorage = typename ScvfGeometryTraits::CornerStorage;
+ using GridIndexStorage = typename ScvfGeometryTraits::GridIndexStorage;
using Geometry = typename ScvfGeometryTraits::Geometry;
using BoundaryFlag = typename ScvfGeometryTraits::BoundaryFlag;
@@ -66,7 +106,7 @@ public:
CCTpfaSubControlVolumeFace(const Intersection& is,
const typename Intersection::Geometry& isGeometry,
GridIndexType scvfIndex,
- const std::vector<GridIndexType>& scvIndices,
+ const GridIndexStorage& scvIndices,
bool isBoundary)
: ParentType()
, geomType_(isGeometry.type())
@@ -165,7 +205,7 @@ private:
GlobalPosition center_;
GlobalPosition unitOuterNormal_;
GridIndexType scvfIndex_;
- std::vector<GridIndexType> scvIndices_;
+ GridIndexStorage scvIndices_;
bool boundary_;
BoundaryFlag boundaryFlag_;
};
--
GitLab