diff --git a/dumux/discretization/cellcentered/mpfa/dualgridindexset.hh b/dumux/discretization/cellcentered/mpfa/dualgridindexset.hh
index a93c7ce63970a589ca92cdf883dd9547556e1ee2..2288f3b49d00aa2d83089d6d30224ae23928162d 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 44b8f899fe6fdb14221087ea24a8bd92bf36fb6b..66928c912115b17854afe58988f2e4c00fd587f2 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 248c0129d9739c625ab49c728ab6f79c0414d3da..c355667c6ce89dae36982ac820648cc26cccf83d 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 ce805e818a37a5fa7f8c5395f284b6fd4d15a26f..aeb2b83479025d7084328b63d81079ad77a1b6d8 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 45bd3cc55a86d4b51e1c0102de357429722ef881..dfc14bfc61d970a2db2f3fbfbdb471e1ee27d455 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 c9521b8bbe4481579a1b7d376edfbda546523c61..1c6cd9db19025d1affe3ecc7fac2e4450fd80e60 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 b77ff57df69394e92020284290216af7897d874f..4ff6ebcdd9f677ee407e0d3ede15d603a76ee542 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 520823347a851b3c135c18374cd7145aa9da41f8..08a559b12ade704f4c1d5c9cb77aed080d413c71 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 c9262c4b141561ec9b65bdec6de69c0019461f51..1524324e2874c8101019fca8b60ceb4d615192f4 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 58c837af98aa0982817d2e04d1f706a75baa93c9..755dbeef0dd2fa0c3aa276a31b85803fbe924ab7 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_;
};