diff --git a/dumux/common/properties.hh b/dumux/common/properties.hh
index a0ccd5ffdf1d92ddba88b6463f44cc068ed68ed0..3ac71613f4bcfec3db3bc0f8000e45c5ddbb5482 100644
--- a/dumux/common/properties.hh
+++ b/dumux/common/properties.hh
@@ -91,12 +91,10 @@ NEW_PROP_TAG(ElementFluxVariablesCache); //!< A local vector of flux v
NEW_PROP_TAG(GridFluxVariablesCache); //!< The global vector of flux variable containers
NEW_PROP_TAG(EnableGridFluxVariablesCache); //!< specifies if data on flux vars should be saved (faster, but more memory consuming)
NEW_PROP_TAG(GridVariables); //!< The grid variables object managing variable data on the grid (volvars/fluxvars cache)
-NEW_PROP_TAG(MaxNumNeighborsPerScvf); //!< The maximum number of neighboring elements allowed per scvf (for static memory allocation)
/////////////////////////////////////////////////////////////////
// Additional properties used by the cell-centered mpfa schemes:
/////////////////////////////////////////////////////////////////
-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
diff --git a/dumux/discretization/cellcentered/mpfa/darcyslaw.hh b/dumux/discretization/cellcentered/mpfa/darcyslaw.hh
index 60e525c13d9666bbc8ed05fc5c1e43ace2a653ac..86e2e5502fbe5f558d3c765216e2c3ed6f28b015 100644
--- a/dumux/discretization/cellcentered/mpfa/darcyslaw.hh
+++ b/dumux/discretization/cellcentered/mpfa/darcyslaw.hh
@@ -32,8 +32,6 @@
#include <dumux/common/parameters.hh>
#include <dumux/discretization/methods.hh>
-#include <dumux/discretization/methods.hh>
-
namespace Dumux
{
//! forward declaration of the method-specific implementation
@@ -52,7 +50,8 @@ class DarcysLawImplementation<TypeTag, DiscretizationMethods::CCMpfa>
using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
using Element = typename GridView::template Codim<0>::Entity;
- using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry)::LocalView;
+ using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
+ using FVElementGeometry = typename FVGridGeometry::LocalView;
using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);
using ElementFluxVariablesCache = typename GET_PROP_TYPE(TypeTag, ElementFluxVariablesCache);
@@ -95,30 +94,22 @@ class DarcysLawImplementation<TypeTag, DiscretizationMethods::CCMpfa>
static constexpr int numPhases = GET_PROP_VALUE(TypeTag, NumPhases);
using DualGridNodalIndexSet = typename GET_PROP_TYPE(TypeTag, DualGridNodalIndexSet);
- using Stencil = typename DualGridNodalIndexSet::GridStencilType;
+ using Stencil = typename DualGridNodalIndexSet::NodalGridStencilType;
- using MpfaHelper = typename GET_PROP_TYPE(TypeTag, MpfaHelper);
+ using MpfaHelper = typename FVGridGeometry::MpfaHelper;
static constexpr bool considerSecondaryIVs = MpfaHelper::considerSecondaryIVs();
- // In the current implementation of the flux variables cache we cannot make a
- // disctinction between dynamic (e.g. mpfa-o unstructured) and static (e.g.mpfa-l)
- // matrix and vector types, as currently the cache class can only be templated
- // by a type tag (and there can only be one). Currently, pointers to both the
- // primary and secondary iv data is stored. Before accessing it has to be checked
- // whether or not the scvf is embedded in a secondary interaction volume.
using PrimaryInteractionVolume = typename GET_PROP_TYPE(TypeTag, PrimaryInteractionVolume);
using PrimaryIvLocalFaceData = typename PrimaryInteractionVolume::Traits::LocalFaceData;
- using PrimaryIvDataHandle = typename PrimaryInteractionVolume::Traits::DataHandle;
- using PrimaryIvVector = typename PrimaryInteractionVolume::Traits::Vector;
- using PrimaryIvMatrix = typename PrimaryInteractionVolume::Traits::Matrix;
- using PrimaryIvTij = typename PrimaryIvMatrix::row_type;
+ using PrimaryIvDataHandle = typename ElementFluxVariablesCache::PrimaryIvDataHandle;
+ using PrimaryIvCellVector = typename PrimaryInteractionVolume::Traits::MatVecTraits::CellVector;
+ using PrimaryIvTij = typename PrimaryInteractionVolume::Traits::MatVecTraits::TMatrix::row_type;
using SecondaryInteractionVolume = typename GET_PROP_TYPE(TypeTag, SecondaryInteractionVolume);
using SecondaryIvLocalFaceData = typename SecondaryInteractionVolume::Traits::LocalFaceData;
- using SecondaryIvDataHandle = typename SecondaryInteractionVolume::Traits::DataHandle;
- using SecondaryIvVector = typename SecondaryInteractionVolume::Traits::Vector;
- using SecondaryIvMatrix = typename SecondaryInteractionVolume::Traits::Matrix;
- using SecondaryIvTij = typename SecondaryIvMatrix::row_type;
+ using SecondaryIvDataHandle = typename ElementFluxVariablesCache::SecondaryIvDataHandle;
+ using SecondaryIvCellVector = typename SecondaryInteractionVolume::Traits::MatVecTraits::CellVector;
+ using SecondaryIvTij = typename SecondaryInteractionVolume::Traits::MatVecTraits::TMatrix::row_type;
public:
// export the filler type
@@ -238,11 +229,11 @@ class DarcysLawImplementation<TypeTag, DiscretizationMethods::CCMpfa>
//! Coefficients for the cell (& Dirichlet) unknowns in flux expressions (secondary type)
const SecondaryIvTij& advectionTijSecondaryIv() const { return *secondaryTij_; }
- //! The cell (& Dirichlet) pressures within this interaction volume (primary type)
- const PrimaryIvVector& pressuresPrimaryIv(unsigned int phaseIdx) const { return *primaryPj_[phaseIdx]; }
+ //! The cell (& maybe Dirichlet) pressures within this interaction volume (primary type)
+ const PrimaryIvCellVector& pressuresPrimaryIv(unsigned int phaseIdx) const { return *primaryPj_[phaseIdx]; }
- //! The cell (& Dirichlet) pressures within this interaction volume (secondary type)
- const SecondaryIvVector& pressuresSecondaryIv(unsigned int phaseIdx) const { return *secondaryPj_[phaseIdx]; }
+ //! The cell (& maybe Dirichlet) pressures within this interaction volume (secondary type)
+ const SecondaryIvCellVector& pressuresSecondaryIv(unsigned int phaseIdx) const { return *secondaryPj_[phaseIdx]; }
//! The gravitational acceleration for a phase on this scvf
Scalar gravity(unsigned int phaseIdx) const { return g_[phaseIdx]; }
@@ -264,8 +255,8 @@ class DarcysLawImplementation<TypeTag, DiscretizationMethods::CCMpfa>
const SecondaryIvTij* secondaryTij_;
//! The interaction-volume wide phase pressures pj
- std::array<const PrimaryIvVector*, numPhases> primaryPj_;
- std::array<const SecondaryIvVector*, numPhases> secondaryPj_;
+ std::array<const PrimaryIvCellVector*, numPhases> primaryPj_;
+ std::array<const SecondaryIvCellVector*, numPhases> secondaryPj_;
//! Gravitational flux contribution on this face
std::array< Scalar, numPhases > g_;
diff --git a/dumux/discretization/cellcentered/mpfa/dualgridindexset.hh b/dumux/discretization/cellcentered/mpfa/dualgridindexset.hh
index 2288f3b49d00aa2d83089d6d30224ae23928162d..59680aded654242d794f817ba0afad5c9138db51 100644
--- a/dumux/discretization/cellcentered/mpfa/dualgridindexset.hh
+++ b/dumux/discretization/cellcentered/mpfa/dualgridindexset.hh
@@ -38,40 +38,32 @@ namespace Dumux
* \brief Nodal index set for mpfa schemes, constructed
* around grid vertices.
*
- * \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.
+ * \tparam T The traits class to be used
*/
-template< class GI, class LI, int dim, int maxE, int maxB = 2 >
+template< class T >
class CCMpfaDualGridNodalIndexSet
{
- using DimIndexVector = Dune::ReservedVector<GI, dim>;
+ using LI = typename T::LocalIndexType;
+ using GI = typename T::GridIndexType;
+
+ using DimLocalIndexVector = Dune::ReservedVector<LI, T::GridView::dimension>;
+ using ScvfIndicesInScvStorage = typename T::template NodalScvDataStorage< DimLocalIndexVector >;
public:
- //! Export the used index types
- using GridIndexType = GI;
- using LocalIndexType = LI;
+ //! Export the traits type
+ using Traits = T;
- //! 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 index types used
+ using LocalIndexType = LI;
+ using GridIndexType = GI;
//! 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>;
+ using NodalGridStencilType = typename T::template NodalScvDataStorage< GI >;
+ using NodalLocalStencilType = typename T::template NodalScvDataStorage< LI >;
+ using NodalGridScvfStencilType = typename T::template NodalScvfDataStorage< GI >;
//! 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>;
+ using ScvfNeighborLocalIndexSet = typename T::template ScvfNeighborDataStorage< LI >;
//! Constructor
CCMpfaDualGridNodalIndexSet() : numBoundaryScvfs_(0) {}
@@ -80,49 +72,41 @@ public:
template<typename SubControlVolumeFace>
void insert(const SubControlVolumeFace& scvf)
{
- insert(scvf.boundary(),
- scvf.index(),
+ insert(scvf.index(),
scvf.insideScvIdx(),
- scvf.outsideScvIndices());
+ scvf.boundary());
}
//! Inserts scvf data
- template<typename OutsideGridIndexStorage>
- void insert(const bool boundary,
- const GridIndexType scvfIdx,
+ void insert(const GridIndexType scvfIdx,
const GridIndexType insideScvIdx,
- const OutsideGridIndexStorage& outsideScvIndices)
+ const bool boundary)
{
// this should always be called only once per scvf
- assert(std::find(scvfIndices_.begin(), scvfIndices_.end(), scvfIdx ) == scvfIndices_.end()
- && "scvf has already been inserted!");
+ assert(std::count(scvfIndices_.begin(), scvfIndices_.end(), scvfIdx ) == 0 && "scvf already inserted!");
// the local index of the scvf data about to be inserted
const LocalIndexType curScvfLocalIdx = scvfIndices_.size();
- // add grid scvf data
- ScvfNeighborIndexSet scvIndices;
- scvIndices.push_back(insideScvIdx);
- for (const auto& outsideIdx : outsideScvIndices)
- scvIndices.push_back(outsideIdx);
-
// if scvf is on boundary, increase counter
- if (boundary)
- numBoundaryScvfs_++;
+ if (boundary) numBoundaryScvfs_++;
// insert data on the new scv
scvfIndices_.push_back(scvfIdx);
scvfIsOnBoundary_.push_back(boundary);
- scvfNeighborScvIndices_.push_back(scvIndices);
- // if entry for the inside scv exists, append scvf local index, create otherwise
+ // if entry for the inside scv exists append data, create otherwise
auto it = std::find( scvIndices_.begin(), scvIndices_.end(), insideScvIdx );
if (it != scvIndices_.end())
- localScvfIndicesInScv_[ std::distance(scvIndices_.begin(), it) ].push_back(curScvfLocalIdx);
+ {
+ const auto scvIdxLocal = std::distance(scvIndices_.begin(), it);
+ scvfInsideScvIndices_.push_back(scvIdxLocal);
+ localScvfIndicesInScv_[scvIdxLocal].push_back(curScvfLocalIdx);
+ }
else
{
- localScvfIndicesInScv_.push_back({});
- localScvfIndicesInScv_.back().push_back(curScvfLocalIdx);
+ scvfInsideScvIndices_.push_back(scvIndices_.size());
+ localScvfIndicesInScv_.push_back({curScvfLocalIdx});
scvIndices_.push_back(insideScvIdx);
}
}
@@ -136,49 +120,64 @@ public:
//! returns the number of boundary scvfs around the node
std::size_t numBoundaryScvfs() const { return numBoundaryScvfs_; }
- //! returns the global scv indices connected to this dual grid node
- const GridStencilType& globalScvIndices() const { return scvIndices_; }
+ //! returns the grid scv indices connected to this dual grid node
+ const NodalGridStencilType& globalScvIndices() const { return scvIndices_; }
- //! returns the global scvf indices connected to this dual grid node
- const GridScvfStencilType& globalScvfIndices() const { return scvfIndices_; }
+ //! returns the grid scvf indices connected to this dual grid node
+ const NodalGridScvfStencilType& 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]; }
+ bool scvfIsOnBoundary(unsigned int i) const
+ {
+ assert(i < numScvfs());
+ return scvfIsOnBoundary_[i];
+ }
- //! returns the global scv idx of the i-th scv
- GridIndexType scvIdxGlobal(unsigned int i) const { return scvIndices_[i]; }
+ //! returns the grid scv idx of the i-th scv
+ GridIndexType scvIdxGlobal(unsigned int i) const
+ {
+ assert(i < numScvs());
+ return scvIndices_[i];
+ }
//! returns the index of the i-th scvf
- GridIndexType scvfIdxGlobal(unsigned int i) const { return scvfIndices_[i]; }
+ GridIndexType scvfIdxGlobal(unsigned int i) const
+ {
+ assert(i < numScvfs());
+ return scvfIndices_[i];
+ }
- //! returns the global index of the j-th scvf embedded in the i-th scv
+ //! returns the grid index of the j-th scvf embedded in the i-th scv
GridIndexType scvfIdxGlobal(unsigned int i, unsigned int j) const
{
- assert(j < dim); // only dim faces can be embedded in an scv
+ assert(i < numScvs());
+ assert(j < localScvfIndicesInScv_[i].size());
return scvfIndices_[ localScvfIndicesInScv_[i][j] ];
}
//! returns the node-local index of the j-th scvf embedded in the i-th scv
LocalIndexType scvfIdxLocal(unsigned int i, unsigned int j) const
{
- assert(j < dim); // only dim faces can be embedded in an scv
+ assert(i < numScvs());
+ assert(j < localScvfIndicesInScv_[i].size());
return localScvfIndicesInScv_[i][j];
}
- //! returns the indices of the neighboring scvs of the i-th scvf
- const ScvfNeighborIndexSet& neighboringScvIndices(unsigned int i) const
- { return scvfNeighborScvIndices_[i]; }
+ //! returns the node-local index of the inside scv of the i-th scvf
+ LocalIndexType insideScvIdxLocal(unsigned int i) const
+ {
+ assert(i < numScvfs());
+ return scvfInsideScvIndices_[i];
+ }
private:
- 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
+ NodalGridStencilType scvIndices_; //!< The indices of the scvs around a dual grid node
+ ScvfIndicesInScvStorage localScvfIndicesInScv_; //!< Maps to each scv a list of scvf indices embedded in it
- //! The neighboring scvs for the scvfs (order: 0 - inside, 1..n - outside)
- Dune::ReservedVector<ScvfNeighborIndexSet, maxNumScvfsAtNode> scvfNeighborScvIndices_;
+ std::size_t numBoundaryScvfs_; //!< stores how many boundary scvfs are embedded in this dual grid node
+ NodalGridScvfStencilType scvfIndices_; //!< the indices of the scvfs around a dual grid node
+ typename T::template NodalScvfDataStorage< bool > scvfIsOnBoundary_; //!< Maps to each scvf a boolean to indicate if it is on the boundary
+ typename T::template NodalScvfDataStorage< LI > scvfInsideScvIndices_; //!< The inside local scv index for each scvf
};
/*!
@@ -199,7 +198,7 @@ public:
//! Constructor taking a grid view
template< class GridView >
- CCMpfaDualGridIndexSet(const GridView& gridView) : nodalIndexSets_(gridView.size(NodalIndexSet::dimension)) {}
+ CCMpfaDualGridIndexSet(const GridView& gridView) : nodalIndexSets_(gridView.size(GridView::dimension)) {}
//! Access with an scvf
template< class SubControlVolumeFace >
diff --git a/dumux/discretization/cellcentered/mpfa/elementfluxvariablescache.hh b/dumux/discretization/cellcentered/mpfa/elementfluxvariablescache.hh
index f7709fe7878e8e776b4cf20f3649e6b26d786864..f1b766a4f995a507c42a19ebc27b1b7fdb5327e6 100644
--- a/dumux/discretization/cellcentered/mpfa/elementfluxvariablescache.hh
+++ b/dumux/discretization/cellcentered/mpfa/elementfluxvariablescache.hh
@@ -29,6 +29,8 @@
#include <dumux/common/properties.hh>
#include <dumux/common/parameters.hh>
+#include <dumux/discretization/cellcentered/mpfa/interactionvolumedatahandle.hh>
+
#include "fluxvariablescachefiller.hh"
#include "methods.hh"
@@ -63,8 +65,11 @@ class CCMpfaElementFluxVariablesCache<TypeTag, true>
using FluxVariablesCache = typename GET_PROP_TYPE(TypeTag, FluxVariablesCache);
using GridFluxVariablesCache = typename GET_PROP_TYPE(TypeTag, GridFluxVariablesCache);
-
public:
+ //! export the data handle types used by the grid-wide cache
+ using PrimaryIvDataHandle = typename GridFluxVariablesCache::PrimaryIvDataHandle;
+ using SecondaryIvDataHandle = typename GridFluxVariablesCache::SecondaryIvDataHandle;
+
//! The constructor
CCMpfaElementFluxVariablesCache(const GridFluxVariablesCache& global)
: gridFluxVarsCachePtr_(&global) {}
@@ -121,11 +126,26 @@ class CCMpfaElementFluxVariablesCache<TypeTag, false>
using FluxVariablesCacheFiller = CCMpfaFluxVariablesCacheFiller<TypeTag>;
using PrimaryInteractionVolume = typename GET_PROP_TYPE(TypeTag, PrimaryInteractionVolume);
- using PrimaryIvDataHandle = typename PrimaryInteractionVolume::Traits::DataHandle;
+ using PrimaryMatVecTraits = typename PrimaryInteractionVolume::Traits::MatVecTraits;
using SecondaryInteractionVolume = typename GET_PROP_TYPE(TypeTag, SecondaryInteractionVolume);
- using SecondaryIvDataHandle = typename SecondaryInteractionVolume::Traits::DataHandle;
+ using SecondaryMatVecTraits = typename SecondaryInteractionVolume::Traits::MatVecTraits;
+
+ //! physics traits class to define the data handles
+ struct PhysicsTraits
+ {
+ static constexpr bool enableAdvection = GET_PROP_VALUE(TypeTag, EnableAdvection);
+ static constexpr bool enableMolecularDiffusion = GET_PROP_VALUE(TypeTag, EnableMolecularDiffusion);
+ static constexpr bool enableHeatConduction = GET_PROP_VALUE(TypeTag, EnableEnergyBalance);
+
+ static constexpr int numPhases = GET_PROP_VALUE(TypeTag, NumPhases);
+ static constexpr int numComponents = GET_PROP_VALUE(TypeTag, NumComponents);
+ };
public:
+ //! export the data handle types used
+ using PrimaryIvDataHandle = InteractionVolumeDataHandle< PrimaryMatVecTraits, PhysicsTraits >;
+ using SecondaryIvDataHandle = InteractionVolumeDataHandle< SecondaryMatVecTraits, PhysicsTraits >;
+
CCMpfaElementFluxVariablesCache(const GridFluxVariablesCache& global)
: gridFluxVarsCachePtr_(&global) {}
diff --git a/dumux/discretization/cellcentered/mpfa/fickslaw.hh b/dumux/discretization/cellcentered/mpfa/fickslaw.hh
index ae3979d58deb66b2702ce536db2a6ae3cda9d59a..de03c053c1c5fa5bd584bc56afe2e7c71a5f74ce 100644
--- a/dumux/discretization/cellcentered/mpfa/fickslaw.hh
+++ b/dumux/discretization/cellcentered/mpfa/fickslaw.hh
@@ -46,7 +46,8 @@ class FicksLawImplementation<TypeTag, DiscretizationMethods::CCMpfa>
using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
using Element = typename GridView::template Codim<0>::Entity;
- using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry)::LocalView;
+ using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
+ using FVElementGeometry = typename FVGridGeometry::LocalView;
using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);
@@ -91,30 +92,22 @@ class FicksLawImplementation<TypeTag, DiscretizationMethods::CCMpfa>
class MpfaFicksLawCache
{
using DualGridNodalIndexSet = typename GET_PROP_TYPE(TypeTag, DualGridNodalIndexSet);
- using Stencil = typename DualGridNodalIndexSet::GridStencilType;
+ using Stencil = typename DualGridNodalIndexSet::NodalGridStencilType;
- using MpfaHelper = typename GET_PROP_TYPE(TypeTag, MpfaHelper);
+ using MpfaHelper = typename FVGridGeometry::MpfaHelper;
static constexpr bool considerSecondaryIVs = MpfaHelper::considerSecondaryIVs();
- // In the current implementation of the flux variables cache we cannot make a
- // disctinction between dynamic (e.g. mpfa-o unstructured) and static (e.g.mpfa-l)
- // matrix and vector types, as currently the cache class can only be templated
- // by a type tag (and there can only be one). Currently, pointers to both the
- // primary and secondary iv data is stored. Before accessing it has to be checked
- // whether or not the scvf is embedded in a secondary interaction volume.
using PrimaryInteractionVolume = typename GET_PROP_TYPE(TypeTag, PrimaryInteractionVolume);
using PrimaryIvLocalFaceData = typename PrimaryInteractionVolume::Traits::LocalFaceData;
- using PrimaryIvDataHandle = typename PrimaryInteractionVolume::Traits::DataHandle;
- using PrimaryIvVector = typename PrimaryInteractionVolume::Traits::Vector;
- using PrimaryIvMatrix = typename PrimaryInteractionVolume::Traits::Matrix;
- using PrimaryIvTij = typename PrimaryIvMatrix::row_type;
+ using PrimaryIvDataHandle = typename ElementFluxVariablesCache::PrimaryIvDataHandle;
+ using PrimaryIvCellVector = typename PrimaryInteractionVolume::Traits::MatVecTraits::CellVector;
+ using PrimaryIvTij = typename PrimaryInteractionVolume::Traits::MatVecTraits::TMatrix::row_type;
using SecondaryInteractionVolume = typename GET_PROP_TYPE(TypeTag, SecondaryInteractionVolume);
using SecondaryIvLocalFaceData = typename SecondaryInteractionVolume::Traits::LocalFaceData;
- using SecondaryIvDataHandle = typename SecondaryInteractionVolume::Traits::DataHandle;
- using SecondaryIvVector = typename SecondaryInteractionVolume::Traits::Vector;
- using SecondaryIvMatrix = typename SecondaryInteractionVolume::Traits::Matrix;
- using SecondaryIvTij = typename SecondaryIvMatrix::row_type;
+ using SecondaryIvDataHandle = typename ElementFluxVariablesCache::SecondaryIvDataHandle;
+ using SecondaryIvCellVector = typename SecondaryInteractionVolume::Traits::MatVecTraits::CellVector;
+ using SecondaryIvTij = typename SecondaryInteractionVolume::Traits::MatVecTraits::TMatrix::row_type;
static constexpr int dim = GridView::dimension;
static constexpr int dimWorld = GridView::dimensionworld;
@@ -198,12 +191,12 @@ class FicksLawImplementation<TypeTag, DiscretizationMethods::CCMpfa>
const SecondaryIvTij& diffusionTijSecondaryIv(unsigned int phaseIdx, unsigned int compIdx) const
{ return *secondaryTij_[phaseIdx][compIdx]; }
- //! The cell (& Dirichlet) mole fractions within this interaction volume (primary type)
- const PrimaryIvVector& moleFractionsPrimaryIv(unsigned int phaseIdx, unsigned int compIdx) const
+ //! The cell (& maybe Dirichlet) mole fractions within this interaction volume (primary type)
+ const PrimaryIvCellVector& moleFractionsPrimaryIv(unsigned int phaseIdx, unsigned int compIdx) const
{ return *primaryXj_[phaseIdx][compIdx]; }
- //! The cell (& Dirichlet) mole fractions within this interaction volume (secondary type)
- const SecondaryIvVector& moleFractionsSecondaryIv(unsigned int phaseIdx, unsigned int compIdx) const
+ //! The cell (& maybe Dirichlet) mole fractions within this interaction volume (secondary type)
+ const SecondaryIvCellVector& moleFractionsSecondaryIv(unsigned int phaseIdx, unsigned int compIdx) const
{ return *secondaryXj_[phaseIdx][compIdx]; }
//! The stencils corresponding to the transmissibilities
@@ -217,12 +210,12 @@ class FicksLawImplementation<TypeTag, DiscretizationMethods::CCMpfa>
std::array< std::array<const Stencil*, numComponents>, numPhases > stencil_;
//! The transmissibilities such that f = Tij*xj
- std::array< std::array<const PrimaryIvVector*, numComponents>, numPhases > primaryTij_;
- std::array< std::array<const SecondaryIvVector*, numComponents>, numPhases > secondaryTij_;
+ std::array< std::array<const PrimaryIvCellVector*, numComponents>, numPhases > primaryTij_;
+ std::array< std::array<const SecondaryIvCellVector*, numComponents>, numPhases > secondaryTij_;
//! The interaction-volume wide mole fractions xj
- std::array< std::array<const PrimaryIvVector*, numComponents>, numPhases > primaryXj_;
- std::array< std::array<const SecondaryIvVector*, numComponents>, numPhases > secondaryXj_;
+ std::array< std::array<const PrimaryIvCellVector*, numComponents>, numPhases > primaryXj_;
+ std::array< std::array<const SecondaryIvCellVector*, numComponents>, numPhases > secondaryXj_;
};
public:
diff --git a/dumux/discretization/cellcentered/mpfa/fluxvariablescachefiller.hh b/dumux/discretization/cellcentered/mpfa/fluxvariablescachefiller.hh
index 5a90cc053a209e2b0dd5f93f714733c626088d04..b2196d97c716c9657d12b30f28060eafc4d5ec85 100644
--- a/dumux/discretization/cellcentered/mpfa/fluxvariablescachefiller.hh
+++ b/dumux/discretization/cellcentered/mpfa/fluxvariablescachefiller.hh
@@ -43,20 +43,22 @@ template<class TypeTag>
class CCMpfaFluxVariablesCacheFiller
{
using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
- using MpfaHelper = typename GET_PROP_TYPE(TypeTag, MpfaHelper);
using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
using Element = typename GridView::template Codim<0>::Entity;
- using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry)::LocalView;
+ using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
+ using FVElementGeometry = typename FVGridGeometry::LocalView;
+ using MpfaHelper = typename FVGridGeometry::MpfaHelper;
using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);
+ using ElementFluxVariablesCache = typename GET_PROP_TYPE(TypeTag, ElementFluxVariablesCache);
using FluxVariablesCache = typename GET_PROP_TYPE(TypeTag, FluxVariablesCache);
using PrimaryInteractionVolume = typename GET_PROP_TYPE(TypeTag, PrimaryInteractionVolume);
- using PrimaryDataHandle = typename PrimaryInteractionVolume::Traits::DataHandle;
+ using PrimaryDataHandle = typename ElementFluxVariablesCache::PrimaryIvDataHandle;
using PrimaryLocalFaceData = typename PrimaryInteractionVolume::Traits::LocalFaceData;
using SecondaryInteractionVolume = typename GET_PROP_TYPE(TypeTag, SecondaryInteractionVolume);
- using SecondaryDataHandle = typename SecondaryInteractionVolume::Traits::DataHandle;
+ using SecondaryDataHandle = typename ElementFluxVariablesCache::SecondaryIvDataHandle;
using SecondaryLocalFaceData = typename SecondaryInteractionVolume::Traits::LocalFaceData;
static constexpr int dim = GridView::dimension;
@@ -252,84 +254,8 @@ private:
using AdvectionType = typename GET_PROP_TYPE(TypeTag, AdvectionType);
using AdvectionFiller = typename AdvectionType::Cache::Filler;
- static constexpr auto AdvectionMethod = AdvectionType::myDiscretizationMethod;
- using LambdaFactory = TensorLambdaFactory<TypeTag, AdvectionMethod>;
-
- // skip the following if advection doesn't use mpfa
- if (AdvectionMethod == DiscretizationMethods::CCMpfa)
- {
- // get instance of the interaction volume-local assembler
- using IVTraits = typename InteractionVolume::Traits;
- using IvLocalAssembler = InteractionVolumeAssembler< IVTraits, InteractionVolume::MpfaMethod >;
- IvLocalAssembler localAssembler(problem(), fvGeometry(), elemVolVars());
-
- // Use different assembly if gravity is enabled
- static const bool enableGravity = getParamFromGroup<bool>(GET_PROP_VALUE(TypeTag, ModelParameterGroup), "Problem.EnableGravity");
-
- // Assemble T only if permeability is sol-dependent or if update is forced
- if (forceUpdateAll || soldependentAdvection)
- {
- // distinguish between normal/surface grids (optimized away by compiler)
- if (dim < dimWorld)
- {
- if (enableGravity)
- localAssembler.assembleWithGravity( handle.advectionTout(),
- handle.advectionT(),
- handle.gravityOutside(),
- handle.gravity(),
- handle.advectionCA(),
- handle.advectionA(),
- iv,
- LambdaFactory::getAdvectionLambda() );
-
- else
- localAssembler.assemble( handle.advectionTout(),
- handle.advectionT(),
- iv,
- LambdaFactory::getAdvectionLambda() );
- }
-
- // normal grids
- else
- {
- if (enableGravity)
- localAssembler.assembleWithGravity( handle.advectionT(),
- handle.gravity(),
- handle.advectionCA(),
- iv,
- LambdaFactory::getAdvectionLambda() );
- else
- localAssembler.assemble( handle.advectionT(),
- iv,
- LambdaFactory::getAdvectionLambda() );
- }
- }
-
- // (maybe) only reassemble gravity vector
- else if (enableGravity)
- {
- if (dim == dimWorld)
- localAssembler.assembleGravity( handle.gravity(),
- iv,
- handle.advectionCA(),
- LambdaFactory::getAdvectionLambda() );
- else
- localAssembler.assembleGravity( handle.gravity(),
- handle.gravityOutside(),
- iv,
- handle.advectionCA(),
- handle.advectionA(),
- LambdaFactory::getAdvectionLambda() );
- }
-
- // assemble pressure vectors
- for (unsigned int pIdx = 0; pIdx < GET_PROP_VALUE(TypeTag, NumPhases); ++pIdx)
- {
- const auto& evv = &elemVolVars();
- auto getPressure = [&evv, pIdx] (auto volVarIdx) { return (evv->operator[](volVarIdx)).pressure(pIdx); };
- localAssembler.assemble(handle.pressures(pIdx), iv, getPressure);
- }
- }
+ // fill data in the handle
+ fillAdvectionHandle(iv, handle, forceUpdateAll);
// fill advection caches
for (unsigned int i = 0; i < iv.localFaceData().size(); ++i)
@@ -384,51 +310,19 @@ private:
using DiffusionType = typename GET_PROP_TYPE(TypeTag, MolecularDiffusionType);
using DiffusionFiller = typename DiffusionType::Cache::Filler;
- static constexpr auto DiffusionMethod = DiffusionType::myDiscretizationMethod;
- using LambdaFactory = TensorLambdaFactory<TypeTag, DiffusionMethod>;
-
static constexpr int numPhases = GET_PROP_VALUE(TypeTag, NumPhases);
static constexpr int numComponents = GET_PROP_VALUE(TypeTag, NumComponents);
- // get instance of the interaction volume-local assembler
- using IVTraits = typename InteractionVolume::Traits;
- using IvLocalAssembler = InteractionVolumeAssembler< IVTraits, InteractionVolume::MpfaMethod >;
- IvLocalAssembler localAssembler(problem(), fvGeometry(), elemVolVars());
-
for (unsigned int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
{
for (unsigned int compIdx = 0; compIdx < numComponents; ++compIdx)
{
- if (phaseIdx == compIdx)
+ using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
+ if (compIdx == FluidSystem::getMainComponent(phaseIdx))
continue;
- // solve the local system subject to the diffusion tensor (if uses mpfa)
- if (DiffusionMethod == DiscretizationMethods::CCMpfa)
- {
- // update the context in handle
- handle.setPhaseIndex(phaseIdx);
- handle.setComponentIndex(compIdx);
-
- // assemble T
- if (forceUpdateAll || soldependentDiffusion)
- {
- if (dim < dimWorld)
- localAssembler.assemble( handle.diffusionTout(),
- handle.diffusionT(),
- iv,
- LambdaFactory::getDiffusionLambda(phaseIdx, compIdx) );
- else
- localAssembler. assemble( handle.diffusionT(),
- iv,
- LambdaFactory::getDiffusionLambda(phaseIdx, compIdx) );
- }
-
- // assemble vector of mole fractions
- const auto& evv = &elemVolVars();
- auto getMoleFraction = [&evv, phaseIdx, compIdx] (auto volVarIdx)
- { return (evv->operator[](volVarIdx)).moleFraction(phaseIdx, compIdx); };
- localAssembler.assemble(handle.moleFractions(phaseIdx, compIdx), iv, getMoleFraction);
- }
+ // fill data in the handle
+ fillDiffusionHandle(iv, handle, forceUpdateAll, phaseIdx, compIdx);
// fill diffusion caches
for (unsigned int i = 0; i < iv.localFaceData().size(); ++i)
@@ -487,35 +381,8 @@ private:
using HeatConductionType = typename GET_PROP_TYPE(TypeTag, HeatConductionType);
using HeatConductionFiller = typename HeatConductionType::Cache::Filler;
- static constexpr auto HeatConductionMethod = HeatConductionType::myDiscretizationMethod;
- using LambdaFactory = TensorLambdaFactory<TypeTag, HeatConductionMethod>;
-
- // maybe solve the local system subject to fourier coefficient
- if (HeatConductionMethod == DiscretizationMethods::CCMpfa)
- {
- // get instance of the interaction volume-local assembler
- using IVTraits = typename InteractionVolume::Traits;
- using IvLocalAssembler = InteractionVolumeAssembler< IVTraits, InteractionVolume::MpfaMethod >;
- IvLocalAssembler localAssembler(problem(), fvGeometry(), elemVolVars());
-
- if (forceUpdateAll || soldependentAdvection)
- {
- if (dim < dimWorld)
- localAssembler.assemble( handle.heatConductionTout(),
- handle.heatConductionT(),
- iv,
- LambdaFactory::getHeatConductionLambda() );
- else
- localAssembler.assemble( handle.heatConductionT(),
- iv,
- LambdaFactory::getHeatConductionLambda() );
- }
-
- // assemble vector of temperatures
- const auto& evv = &elemVolVars();
- auto getMoleFraction = [&evv] (auto volVarIdx) { return (evv->operator[](volVarIdx)).temperature(); };
- localAssembler.assemble(handle.temperatures(), iv, getMoleFraction);
- }
+ // prepare data in handle
+ fillHeatConductionHandle(iv, handle, forceUpdateAll);
// fill heat conduction caches
for (unsigned int i = 0; i < iv.localFaceData().size(); ++i)
@@ -556,6 +423,184 @@ private:
bool forceUpdateAll = false)
{}
+ //! prepares the quantities necessary for advective fluxes in the handle
+ template< class InteractionVolume,
+ class DataHandle,
+ class AdvectionType = typename GET_PROP_TYPE(TypeTag, AdvectionType),
+ typename std::enable_if_t<AdvectionType::myDiscretizationMethod == DiscretizationMethods::CCMpfa, int> = 0 >
+ void fillAdvectionHandle(InteractionVolume& iv, DataHandle& handle, bool forceUpdateAll)
+ {
+ using LambdaFactory = TensorLambdaFactory<TypeTag, DiscretizationMethods::CCMpfa>;
+
+ // get instance of the interaction volume-local assembler
+ static constexpr MpfaMethods M = InteractionVolume::MpfaMethod;
+ using IvLocalAssembler = InteractionVolumeAssembler< Problem, FVElementGeometry, ElementVolumeVariables, M >;
+ IvLocalAssembler localAssembler(problem(), fvGeometry(), elemVolVars());
+
+ // Use different assembly if gravity is enabled
+ static const bool enableGravity = getParamFromGroup<bool>(GET_PROP_VALUE(TypeTag, ModelParameterGroup), "Problem.EnableGravity");
+
+ // Assemble T only if permeability is sol-dependent or if update is forced
+ if (forceUpdateAll || soldependentAdvection)
+ {
+ // distinguish between normal/surface grids (optimized away by compiler)
+ if (dim < dimWorld)
+ {
+ if (enableGravity)
+ localAssembler.assembleWithGravity( handle.advectionTout(),
+ handle.advectionT(),
+ handle.gravityOutside(),
+ handle.gravity(),
+ handle.advectionCA(),
+ handle.advectionA(),
+ iv,
+ LambdaFactory::getAdvectionLambda() );
+
+ else
+ localAssembler.assemble( handle.advectionTout(),
+ handle.advectionT(),
+ iv,
+ LambdaFactory::getAdvectionLambda() );
+ }
+
+ // normal grids
+ else
+ {
+ if (enableGravity)
+ localAssembler.assembleWithGravity( handle.advectionT(),
+ handle.gravity(),
+ handle.advectionCA(),
+ iv,
+ LambdaFactory::getAdvectionLambda() );
+ else
+ localAssembler.assemble( handle.advectionT(),
+ iv,
+ LambdaFactory::getAdvectionLambda() );
+ }
+ }
+
+ // (maybe) only reassemble gravity vector
+ else if (enableGravity)
+ {
+ if (dim == dimWorld)
+ localAssembler.assembleGravity( handle.gravity(),
+ iv,
+ handle.advectionCA(),
+ LambdaFactory::getAdvectionLambda() );
+ else
+ localAssembler.assembleGravity( handle.gravity(),
+ handle.gravityOutside(),
+ iv,
+ handle.advectionCA(),
+ handle.advectionA(),
+ LambdaFactory::getAdvectionLambda() );
+ }
+
+ // assemble pressure vectors
+ for (unsigned int pIdx = 0; pIdx < GET_PROP_VALUE(TypeTag, NumPhases); ++pIdx)
+ {
+ const auto& evv = &elemVolVars();
+ auto getPressure = [&evv, pIdx] (auto volVarIdx) { return (evv->operator[](volVarIdx)).pressure(pIdx); };
+ localAssembler.assemble(handle.pressures(pIdx), iv, getPressure);
+ }
+ }
+
+ //! prepares the quantities necessary for diffusive fluxes in the handle
+ template< class InteractionVolume,
+ class DataHandle,
+ class DiffusionType = typename GET_PROP_TYPE(TypeTag, MolecularDiffusionType),
+ typename std::enable_if_t<DiffusionType::myDiscretizationMethod == DiscretizationMethods::CCMpfa, int> = 0 >
+ void fillDiffusionHandle(InteractionVolume& iv,
+ DataHandle& handle,
+ bool forceUpdateAll,
+ int phaseIdx, int compIdx)
+ {
+ using LambdaFactory = TensorLambdaFactory<TypeTag, DiscretizationMethods::CCMpfa>;
+
+ // get instance of the interaction volume-local assembler
+ static constexpr MpfaMethods M = InteractionVolume::MpfaMethod;
+ using IvLocalAssembler = InteractionVolumeAssembler< Problem, FVElementGeometry, ElementVolumeVariables, M >;
+ IvLocalAssembler localAssembler(problem(), fvGeometry(), elemVolVars());
+
+ // solve the local system subject to the tensor and update the handle
+ handle.setPhaseIndex(phaseIdx);
+ handle.setComponentIndex(compIdx);
+
+ // assemble T
+ if (forceUpdateAll || soldependentDiffusion)
+ {
+ if (dim < dimWorld)
+ localAssembler.assemble( handle.diffusionTout(),
+ handle.diffusionT(),
+ iv,
+ LambdaFactory::getDiffusionLambda(phaseIdx, compIdx) );
+ else
+ localAssembler. assemble( handle.diffusionT(),
+ iv,
+ LambdaFactory::getDiffusionLambda(phaseIdx, compIdx) );
+ }
+
+ // assemble vector of mole fractions
+ const auto& evv = &elemVolVars();
+ auto getMoleFraction = [&evv, phaseIdx, compIdx] (auto volVarIdx)
+ { return (evv->operator[](volVarIdx)).moleFraction(phaseIdx, compIdx); };
+ localAssembler.assemble(handle.moleFractions(phaseIdx, compIdx), iv, getMoleFraction);
+ }
+
+ //! prepares the quantities necessary for conductive fluxes in the handle
+ template< class InteractionVolume,
+ class DataHandle,
+ class HeatConductionType = typename GET_PROP_TYPE(TypeTag, HeatConductionType),
+ typename std::enable_if_t<HeatConductionType::myDiscretizationMethod == DiscretizationMethods::CCMpfa, int> = 0 >
+ void fillHeatConductionHandle(InteractionVolume& iv, DataHandle& handle, bool forceUpdateAll)
+ {
+ using LambdaFactory = TensorLambdaFactory<TypeTag, DiscretizationMethods::CCMpfa>;
+
+ // get instance of the interaction volume-local assembler
+ static constexpr MpfaMethods M = InteractionVolume::MpfaMethod;
+ using IvLocalAssembler = InteractionVolumeAssembler< Problem, FVElementGeometry, ElementVolumeVariables, M >;
+ IvLocalAssembler localAssembler(problem(), fvGeometry(), elemVolVars());
+
+ if (forceUpdateAll || soldependentAdvection)
+ {
+ if (dim < dimWorld)
+ localAssembler.assemble( handle.heatConductionTout(),
+ handle.heatConductionT(),
+ iv,
+ LambdaFactory::getHeatConductionLambda() );
+ else
+ localAssembler.assemble( handle.heatConductionT(),
+ iv,
+ LambdaFactory::getHeatConductionLambda() );
+ }
+
+ // assemble vector of temperatures
+ const auto& evv = &elemVolVars();
+ auto getMoleFraction = [&evv] (auto volVarIdx) { return (evv->operator[](volVarIdx)).temperature(); };
+ localAssembler.assemble(handle.temperatures(), iv, getMoleFraction);
+ }
+
+ //! fill handle only when advection uses mpfa
+ template< class InteractionVolume,
+ class DataHandle,
+ class AdvectionType = typename GET_PROP_TYPE(TypeTag, AdvectionType),
+ typename std::enable_if_t<AdvectionType::myDiscretizationMethod != DiscretizationMethods::CCMpfa, int> = 0 >
+ void fillAdvectionHandle(InteractionVolume& iv, DataHandle& handle, bool forceUpdateAll) {}
+
+ //! fill handle only when diffusion uses mpfa
+ template< class InteractionVolume,
+ class DataHandle,
+ class DiffusionType = typename GET_PROP_TYPE(TypeTag, MolecularDiffusionType),
+ typename std::enable_if_t<DiffusionType::myDiscretizationMethod != DiscretizationMethods::CCMpfa, int> = 0 >
+ void fillDiffusionHandle(InteractionVolume& iv, DataHandle& handle, bool forceUpdateAll, int phaseIdx, int compIdx) {}
+
+ //! fill handle only when heat conduction uses mpfa
+ template< class InteractionVolume,
+ class DataHandle,
+ class HeatConductionType = typename GET_PROP_TYPE(TypeTag, HeatConductionType),
+ typename std::enable_if_t<HeatConductionType::myDiscretizationMethod != DiscretizationMethods::CCMpfa, int> = 0 >
+ void fillHeatConductionHandle(InteractionVolume& iv, DataHandle& handle, bool forceUpdateAll) {}
+
const Problem* problemPtr_;
const Element* elementPtr_;
const FVElementGeometry* fvGeometryPtr_;
diff --git a/dumux/discretization/cellcentered/mpfa/fourierslaw.hh b/dumux/discretization/cellcentered/mpfa/fourierslaw.hh
index 56e3af6da2d0625c0aca32d27c8afbec6f26cb41..aa153053a6447c327663794f6f982a052b5fdf90 100644
--- a/dumux/discretization/cellcentered/mpfa/fourierslaw.hh
+++ b/dumux/discretization/cellcentered/mpfa/fourierslaw.hh
@@ -49,7 +49,9 @@ class FouriersLawImplementation<TypeTag, DiscretizationMethods::CCMpfa>
using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
using Element = typename GridView::template Codim<0>::Entity;
- using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry)::LocalView;
+
+ using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
+ using FVElementGeometry = typename FVGridGeometry::LocalView;
using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);
using ElementFluxVarsCache = typename GET_PROP_TYPE(TypeTag, ElementFluxVariablesCache);
@@ -89,30 +91,22 @@ class FouriersLawImplementation<TypeTag, DiscretizationMethods::CCMpfa>
class MpfaFouriersLawCache
{
using DualGridNodalIndexSet = typename GET_PROP_TYPE(TypeTag, DualGridNodalIndexSet);
- using Stencil = typename DualGridNodalIndexSet::GridStencilType;
+ using Stencil = typename DualGridNodalIndexSet::NodalGridStencilType;
- using MpfaHelper = typename GET_PROP_TYPE(TypeTag, MpfaHelper);
+ using MpfaHelper = typename FVGridGeometry::MpfaHelper;
static constexpr bool considerSecondaryIVs = MpfaHelper::considerSecondaryIVs();
- // In the current implementation of the flux variables cache we cannot make a
- // disctinction between dynamic (e.g. mpfa-o unstructured) and static (e.g.mpfa-l)
- // matrix and vector types, as currently the cache class can only be templated
- // by a type tag (and there can only be one). Currently, pointers to both the
- // primary and secondary iv data is stored. Before accessing it has to be checked
- // whether or not the scvf is embedded in a secondary interaction volume.
using PrimaryInteractionVolume = typename GET_PROP_TYPE(TypeTag, PrimaryInteractionVolume);
using PrimaryIvLocalFaceData = typename PrimaryInteractionVolume::Traits::LocalFaceData;
- using PrimaryIvDataHandle = typename PrimaryInteractionVolume::Traits::DataHandle;
- using PrimaryIvVector = typename PrimaryInteractionVolume::Traits::Vector;
- using PrimaryIvMatrix = typename PrimaryInteractionVolume::Traits::Matrix;
- using PrimaryIvTij = typename PrimaryIvMatrix::row_type;
+ using PrimaryIvDataHandle = typename ElementFluxVarsCache::PrimaryIvDataHandle;
+ using PrimaryIvCellVector = typename PrimaryInteractionVolume::Traits::MatVecTraits::CellVector;
+ using PrimaryIvTij = typename PrimaryInteractionVolume::Traits::MatVecTraits::TMatrix::row_type;
using SecondaryInteractionVolume = typename GET_PROP_TYPE(TypeTag, SecondaryInteractionVolume);
using SecondaryIvLocalFaceData = typename SecondaryInteractionVolume::Traits::LocalFaceData;
- using SecondaryIvDataHandle = typename SecondaryInteractionVolume::Traits::DataHandle;
- using SecondaryIvVector = typename SecondaryInteractionVolume::Traits::Vector;
- using SecondaryIvMatrix = typename SecondaryInteractionVolume::Traits::Matrix;
- using SecondaryIvTij = typename SecondaryIvMatrix::row_type;
+ using SecondaryIvDataHandle = typename ElementFluxVarsCache::SecondaryIvDataHandle;
+ using SecondaryIvCellVector = typename SecondaryInteractionVolume::Traits::MatVecTraits::CellVector;
+ using SecondaryIvTij = typename SecondaryInteractionVolume::Traits::MatVecTraits::TMatrix::row_type;
static constexpr int dim = GridView::dimension;
static constexpr int dimWorld = GridView::dimensionworld;
@@ -188,10 +182,10 @@ class FouriersLawImplementation<TypeTag, DiscretizationMethods::CCMpfa>
const Stencil& heatConductionStencil() const { return *stencil_; }
//! The cell (& Dirichlet) temperatures within this interaction volume (primary type)
- const PrimaryIvVector& temperaturesPrimaryIv() const { return *primaryTj_; }
+ const PrimaryIvCellVector& temperaturesPrimaryIv() const { return *primaryTj_; }
//! The cell (& Dirichlet) temperatures within this interaction volume (secondary type)
- const SecondaryIvVector& temperaturesSecondaryIv() const { return *secondaryTj_; }
+ const SecondaryIvCellVector& temperaturesSecondaryIv() const { return *secondaryTj_; }
//! In the interaction volume-local system of eq we have one unknown per face.
//! On scvfs on this face, but in "outside" (neighbor) elements of it, we have
@@ -210,8 +204,8 @@ class FouriersLawImplementation<TypeTag, DiscretizationMethods::CCMpfa>
const SecondaryIvTij* secondaryTij_;
//! The interaction-volume wide temperature Tj
- const PrimaryIvVector* primaryTj_;
- const SecondaryIvVector* secondaryTj_;
+ const PrimaryIvCellVector* primaryTj_;
+ const SecondaryIvCellVector* secondaryTj_;
};
public:
diff --git a/dumux/discretization/cellcentered/mpfa/fvelementgeometry.hh b/dumux/discretization/cellcentered/mpfa/fvelementgeometry.hh
index 557b70dad7da825c71ab7501130d3f6c60dec7b1..1d0a9f19c829b3a1243924fbaec1046e23561144 100644
--- a/dumux/discretization/cellcentered/mpfa/fvelementgeometry.hh
+++ b/dumux/discretization/cellcentered/mpfa/fvelementgeometry.hh
@@ -222,22 +222,7 @@ public:
//! Note that e.g. the normals might be different in the case of surface grids
const SubControlVolumeFace& flipScvf(GridIndexType scvfIdx, unsigned int outsideScvIdx = 0) const
{
- auto it = std::find(scvfIndices_.begin(), scvfIndices_.end(), scvfIdx);
- if (it != scvfIndices_.end())
- {
- const auto localScvfIdx = std::distance(scvfIndices_.begin(), it);
- return neighborScvfs_[flipScvfIndices_[localScvfIdx][outsideScvIdx]];
- }
- else
- {
- const auto neighborScvfIdxLocal = findLocalIndex(scvfIdx, neighborScvfIndices_);
- const auto& scvf = neighborScvfs_[neighborScvfIdxLocal];
-
- if (scvf.outsideScvIdx(outsideScvIdx) == scvIndices_[0])
- return scvfs_[neighborFlipScvfIndices_[neighborScvfIdxLocal][outsideScvIdx]];
- else
- return neighborScvfs_[neighborFlipScvfIndices_[neighborScvfIdxLocal][outsideScvIdx]];
- }
+ return scvf( fvGridGeometry().flipScvfIdx(scvfIdx, outsideScvIdx) );
}
//! iterator range for sub control volumes. Iterates over
@@ -299,10 +284,6 @@ public:
dataJ.scvfsJ,
dataJ.additionalScvfs);
- // set up the scvf flip indices for network grids
- if (dim < dimWorld)
- makeFlipIndexSet();
-
// //! TODO Check if user added additional DOF dependencies, i.e. the residual of DOF globalI depends
// //! on additional DOFs not included in the discretization schemes' occupation pattern
// const auto& additionalDofDependencies = problem.getAdditionalDofDependencies(globalI);
@@ -524,98 +505,6 @@ private:
}
}
- //! find the "flip" indices for all scvfs
- void makeFlipIndexSet()
- {
- const auto numInsideScvfs = scvfIndices_.size();
- const auto numNeighborScvfs = neighborScvfIndices_.size();
-
- // first, handle the interior scvfs (flip scvf will be in outside scvfs)
- flipScvfIndices_.resize(numInsideScvfs);
- for (unsigned int insideFace = 0; insideFace < numInsideScvfs; ++insideFace)
- {
- const auto& scvf = scvfs_[insideFace];
- if (scvf.boundary())
- continue;
-
- const auto numOutsideScvs = scvf.numOutsideScvs();
- const auto vIdxGlobal = scvf.vertexIndex();
- const auto insideScvIdx = scvf.insideScvIdx();
-
- flipScvfIndices_[insideFace].resize(numOutsideScvs);
- for (unsigned int nIdx = 0; nIdx < numOutsideScvs; ++nIdx)
- {
- const auto outsideScvIdx = scvf.outsideScvIdx(nIdx);
- for (unsigned int outsideFace = 0; outsideFace < numNeighborScvfs; ++outsideFace)
- {
- const auto& outsideScvf = neighborScvfs_[outsideFace];
- // check if outside face is the flip face
- if (outsideScvf.insideScvIdx() == outsideScvIdx &&
- outsideScvf.vertexIndex() == vIdxGlobal &&
- MpfaHelper::vectorContainsValue(outsideScvf.outsideScvIndices(), insideScvIdx))
- {
- flipScvfIndices_[insideFace][nIdx] = outsideFace;
- // there is always only one flip face in an outside element
- break;
- }
- }
- }
- }
-
- // Now we look for the flip indices of the outside faces
- neighborFlipScvfIndices_.resize(numNeighborScvfs);
- for (unsigned int outsideFace = 0; outsideFace < numNeighborScvfs; ++outsideFace)
- {
- const auto& scvf = neighborScvfs_[outsideFace];
- if (scvf.boundary())
- continue;
-
- const auto numOutsideScvs = scvf.numOutsideScvs();
- const auto vIdxGlobal = scvf.vertexIndex();
- const auto insideScvIdx = scvf.insideScvIdx();
-
- neighborFlipScvfIndices_[outsideFace].resize(numOutsideScvs);
- for (unsigned int nIdx = 0; nIdx < numOutsideScvs; ++nIdx)
- {
- const auto neighborScvIdx = scvf.outsideScvIdx(nIdx);
-
- // if the neighbor scv idx is the index of the bound element,
- // then the flip scvf will be within the inside scvfs
- if (neighborScvIdx == scvIndices_[0])
- {
- for (unsigned int insideFace = 0; insideFace < numInsideScvfs; ++insideFace)
- {
- const auto& insideScvf = scvfs_[insideFace];
- // check if the face is the flip face
- if (insideScvf.vertexIndex() == vIdxGlobal &&
- MpfaHelper::vectorContainsValue(insideScvf.outsideScvIndices(), insideScvIdx))
- {
- neighborFlipScvfIndices_[outsideFace][nIdx] = insideFace;
- // there is always only one flip face in an outside element
- break;
- }
- }
- }
- else
- {
- for (unsigned int outsideFace2 = 0; outsideFace2 < numNeighborScvfs; ++outsideFace2)
- {
- const auto& outsideScvf = neighborScvfs_[outsideFace2];
- // check if outside face is the flip face
- if (outsideScvf.insideScvIdx() == neighborScvIdx &&
- outsideScvf.vertexIndex() == vIdxGlobal &&
- MpfaHelper::vectorContainsValue(outsideScvf.outsideScvIndices(), insideScvIdx))
- {
- neighborFlipScvfIndices_[outsideFace][nIdx] = outsideFace2;
- // there is always only one flip face in an outside element
- break;
- }
- }
- }
- }
- }
- }
-
//! map a global index to the local storage index
const unsigned int findLocalIndex(const GridIndexType idx,
const std::vector<GridIndexType>& indices) const
@@ -635,9 +524,6 @@ private:
neighborScvfIndices_.clear();
neighborScvs_.clear();
neighborScvfs_.clear();
-
- flipScvfIndices_.clear();
- neighborFlipScvfIndices_.clear();
}
const FVGridGeometry* fvGridGeometryPtr_;
@@ -652,10 +538,6 @@ private:
std::vector<GridIndexType> neighborScvfIndices_;
std::vector<SubControlVolume> neighborScvs_;
std::vector<SubControlVolumeFace> neighborScvfs_;
-
- // flip scvf index sets
- std::vector< std::vector<GridIndexType> > flipScvfIndices_;
- std::vector< std::vector<GridIndexType> > neighborFlipScvfIndices_;
};
} // end namespace
diff --git a/dumux/discretization/cellcentered/mpfa/fvgridgeometry.hh b/dumux/discretization/cellcentered/mpfa/fvgridgeometry.hh
index 36a380d9002300402425a7e8887cde12ef95804e..6a6d7f40b5f2750067f629ac747e126c49debbfa 100644
--- a/dumux/discretization/cellcentered/mpfa/fvgridgeometry.hh
+++ b/dumux/discretization/cellcentered/mpfa/fvgridgeometry.hh
@@ -71,9 +71,13 @@ class CCMpfaFVGridGeometry<GV, Traits, true>
using CoordScalar = typename GV::ctype;
using ReferenceElements = typename Dune::ReferenceElements<CoordScalar, dim>;
+ using ScvfOutsideGridIndexStorage = typename Traits::SubControlVolumeFace::Traits::OutsideGridIndexStorage;
+
public:
+ //! export the flip scvf index set type
+ using FlipScvfIndexSet = std::vector<ScvfOutsideGridIndexStorage>;
//! export the mpfa helper type
- using MpfaHelper = typename Traits::MpfaHelper;
+ using MpfaHelper = typename Traits::template MpfaHelper<ThisType>;
//! export the grid interaction volume index set type
using GridIVIndexSets = typename Traits::template GridIvIndexSets<ThisType>;
//! export the type to be used for indicators where to use the secondary ivs
@@ -192,7 +196,6 @@ 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)
- using ScvfOutsideGridIndexStorage = typename SubControlVolumeFace::Traits::OutsideGridIndexStorage;
std::vector<ScvfOutsideGridIndexStorage> outsideIndices;
if (dim < dimWorld)
{
@@ -301,34 +304,30 @@ public:
}
// Make the flip index set for network and surface grids
- if (dim < dimWorld)
+ flipScvfIndices_.resize(scvfs_.size());
+ for (const auto& scvf : scvfs_)
{
- flipScvfIndices_.resize(scvfs_.size());
- for (const auto& scvf : scvfs_)
- {
- if (scvf.boundary())
- continue;
+ if (scvf.boundary())
+ continue;
- const auto numOutsideScvs = scvf.numOutsideScvs();
- const auto vIdxGlobal = scvf.vertexIndex();
- const auto insideScvIdx = scvf.insideScvIdx();
+ const auto numOutsideScvs = scvf.numOutsideScvs();
+ const auto vIdxGlobal = scvf.vertexIndex();
+ const auto insideScvIdx = scvf.insideScvIdx();
- flipScvfIndices_[scvf.index()].resize(numOutsideScvs);
- for (std::size_t i = 0; i < numOutsideScvs; ++i)
+ flipScvfIndices_[scvf.index()].resize(numOutsideScvs);
+ for (std::size_t i = 0; i < numOutsideScvs; ++i)
+ {
+ const auto outsideScvIdx = scvf.outsideScvIdx(i);
+ for (auto outsideScvfIndex : scvfIndicesOfScv_[outsideScvIdx])
{
- const auto outsideScvIdx = scvf.outsideScvIdx(i);
- for (auto outsideScvfIndex : scvfIndicesOfScv_[outsideScvIdx])
+ const auto& outsideScvf = this->scvf(outsideScvfIndex);
+ if (outsideScvf.vertexIndex() == vIdxGlobal &&
+ MpfaHelper::vectorContainsValue(outsideScvf.outsideScvIndices(), insideScvIdx))
{
- const auto& outsideScvf = this->scvf(outsideScvfIndex);
- if (outsideScvf.vertexIndex() == vIdxGlobal &&
- MpfaHelper::vectorContainsValue(outsideScvf.outsideScvIndices(), insideScvIdx))
- {
- flipScvfIndices_[scvf.index()][i] = outsideScvfIndex;
- // there is always only one flip face in an outside element
- break;
- }
+ flipScvfIndices_[scvf.index()][i] = outsideScvfIndex;
+ // there is always only one flip face in an outside element
+ break;
}
-
}
}
}
@@ -347,6 +346,9 @@ public:
std::cout << "Initializing of the connectivity map took " << timer.elapsed() << " seconds." << std::endl;
}
+ //! Returns instance of the mpfa helper type
+ MpfaHelper mpfaHelper() const { return MpfaHelper(); }
+
//! Get a sub control volume with a global scv index
const SubControlVolume& scv(GridIndexType scvIdx) const { return scvs_[scvIdx]; }
@@ -360,15 +362,17 @@ public:
//! Returns the grid interaction volume index set class.
const GridIVIndexSets& gridInteractionVolumeIndexSets() const { return ivIndexSets_; }
+ //! Get the sub control volume face indices of an scv by global index
+ const std::vector<GridIndexType>& scvfIndicesOfScv(GridIndexType scvIdx) const { return scvfIndicesOfScv_[scvIdx]; }
+
+ //! Returns the flip scvf index set
+ const FlipScvfIndexSet& flipScvfIndexSet() const { return flipScvfIndices_; }
+
//! Get the scvf on the same face but from the other side
//! Note that e.g. the normals might be different in the case of surface grids
const SubControlVolumeFace& flipScvf(GridIndexType scvfIdx, unsigned int outsideScvfIdx = 0) const
{ return scvfs_[flipScvfIndices_[scvfIdx][outsideScvfIdx]]; }
- //! Get the sub control volume face indices of an scv by global index
- const std::vector<GridIndexType>& scvfIndicesOfScv(GridIndexType scvIdx) const
- { return scvfIndicesOfScv_[scvIdx]; }
-
private:
// connectivity map for efficient assembly
ConnectivityMap connectivityMap_;
@@ -383,7 +387,7 @@ private:
GridIndexType numBoundaryScvf_;
// needed for embedded surface and network grids (dim < dimWorld)
- std::vector<std::vector<GridIndexType>> flipScvfIndices_;
+ FlipScvfIndexSet flipScvfIndices_;
// The grid interaction volume index set
GridIVIndexSets ivIndexSets_;
@@ -419,8 +423,10 @@ class CCMpfaFVGridGeometry<GV, Traits, false>
using ScvfOutsideGridIndexStorage = typename Traits::SubControlVolumeFace::Traits::OutsideGridIndexStorage;
public:
+ //! export the flip scvf index set type
+ using FlipScvfIndexSet = std::vector<ScvfOutsideGridIndexStorage>;
//! export the mpfa helper type
- using MpfaHelper = typename Traits::MpfaHelper;
+ using MpfaHelper = typename Traits::template MpfaHelper<ThisType>;
//! export the grid interaction volume index set type
using GridIVIndexSets = typename Traits::template GridIvIndexSets<ThisType>;
//! export the type to be used for indicators where to use the secondary ivs
@@ -523,6 +529,13 @@ public:
// instantiate the dual grid index set (to be used for construction of interaction volumes)
typename GridIVIndexSets::DualGridIndexSet dualIdSet(this->gridView());
+ // keep track of boundary scvfs and scvf vertex indices in order to set up flip scvf index set
+ const auto maxNumScvfs = numScvs_*LocalView::maxNumElementScvfs;
+ std::vector<bool> scvfIsOnBoundary;
+ std::vector<GridIndexType> scvfVertexIndex;
+ scvfIsOnBoundary.reserve(maxNumScvfs);
+ scvfVertexIndex.reserve(maxNumScvfs);
+
// Build the SCVs and SCV faces
numScvf_ = 0;
numBoundaryScvf_ = 0;
@@ -609,10 +622,12 @@ public:
} ();
// insert the scvf data into the dual grid index set
- dualIdSet[vIdxGlobal].insert(boundary, numScvf_, eIdx, outsideScvIndices);
+ dualIdSet[vIdxGlobal].insert(numScvf_, eIdx, boundary);
// store information on the scv face
scvfsIndexSet.push_back(numScvf_++);
+ scvfIsOnBoundary.push_back(boundary);
+ scvfVertexIndex.push_back(vIdxGlobal);
neighborVolVarIndexSet.emplace_back(std::move(outsideScvIndices));
}
@@ -626,6 +641,43 @@ public:
neighborVolVarIndices_[eIdx] = neighborVolVarIndexSet;
}
+ // Make the flip scvf index set
+ flipScvfIndices_.resize(numScvf_);
+ for (std::size_t scvIdx = 0; scvIdx < numScvs_; ++scvIdx)
+ {
+ const auto& scvfIndices = scvfIndicesOfScv_[scvIdx];
+ for (unsigned int i = 0; i < scvfIndices.size(); ++i)
+ {
+ // boundary scvf have no flip scvfs
+ if (scvfIsOnBoundary[ scvfIndices[i] ])
+ continue;
+
+ const auto scvfIdx = scvfIndices[i];
+ const auto vIdxGlobal = scvfVertexIndex[scvfIdx];
+ const auto numOutsideScvs = neighborVolVarIndices_[scvIdx][i].size();
+
+ flipScvfIndices_[scvfIdx].resize(numOutsideScvs);
+ for (unsigned int j = 0; j < numOutsideScvs; ++j)
+ {
+ const auto outsideScvIdx = neighborVolVarIndices_[scvIdx][i][j];
+ const auto& outsideScvfIndices = scvfIndicesOfScv_[outsideScvIdx];
+ for (unsigned int k = 0; k < outsideScvfIndices.size(); ++k)
+ {
+ const auto outsideScvfIndex = outsideScvfIndices[k];
+ const auto outsideScvfVertexIndex = scvfVertexIndex[outsideScvfIndex];
+ const auto& outsideScvfNeighborIndices = neighborVolVarIndices_[outsideScvIdx][k];
+ if (outsideScvfVertexIndex == vIdxGlobal &&
+ MpfaHelper::vectorContainsValue(outsideScvfNeighborIndices, scvIdx))
+ {
+ flipScvfIndices_[scvfIdx][j] = outsideScvfIndex;
+ // there is always only one flip face in an outside element
+ break;
+ }
+ }
+ }
+ }
+ }
+
// building the geometries has finished
std::cout << "Initializing of the grid finite volume geometry took " << timer.elapsed() << " seconds." << std::endl;
@@ -640,6 +692,9 @@ public:
std::cout << "Initializing of the connectivity map took " << timer.elapsed() << " seconds." << std::endl;
}
+ //! Returns instance of the mpfa helper type
+ MpfaHelper mpfaHelper() const { return MpfaHelper(); }
+
//! Returns the sub control volume face indices of an scv by global index.
const std::vector<GridIndexType>& scvfIndicesOfScv(GridIndexType scvIdx) const
{ return scvfIndicesOfScv_[scvIdx]; }
@@ -648,6 +703,14 @@ public:
const std::vector<ScvfOutsideGridIndexStorage>& neighborVolVarIndices(GridIndexType scvIdx) const
{ return neighborVolVarIndices_[scvIdx]; }
+ //! Get the index scvf on the same face but from the other side
+ //! Note that e.g. the normals might be different in the case of surface grids
+ const GridIndexType flipScvfIdx(GridIndexType scvfIdx, unsigned int outsideScvfIdx = 0) const
+ { return flipScvfIndices_[scvfIdx][outsideScvfIdx]; }
+
+ //! Returns the flip scvf index set
+ const FlipScvfIndexSet& flipScvfIndexSet() const { return flipScvfIndices_; }
+
//! Returns the connectivity map of which dofs
//! have derivatives with respect to a given dof.
const ConnectivityMap& connectivityMap() const { return connectivityMap_; }
@@ -668,6 +731,9 @@ private:
GridIndexType numScvf_;
GridIndexType numBoundaryScvf_;
+ // needed for embedded surface and network grids (dim < dimWorld)
+ FlipScvfIndexSet flipScvfIndices_;
+
// The grid interaction volume index set
GridIVIndexSets ivIndexSets_;
diff --git a/dumux/discretization/cellcentered/mpfa/gridfluxvariablescache.hh b/dumux/discretization/cellcentered/mpfa/gridfluxvariablescache.hh
index 87392c2b4f71f4a911e799b90394a49f35bde083..a43048d8e8ad1db4874a1b7d0890c26d273ee8dc 100644
--- a/dumux/discretization/cellcentered/mpfa/gridfluxvariablescache.hh
+++ b/dumux/discretization/cellcentered/mpfa/gridfluxvariablescache.hh
@@ -25,6 +25,7 @@
#define DUMUX_DISCRETIZATION_CCMPFA_GRID_FLUXVARSCACHE_HH
#include <dumux/common/properties.hh>
+#include <dumux/discretization/cellcentered/mpfa/interactionvolumedatahandle.hh>
#include <dumux/discretization/cellcentered/mpfa/fluxvariablescachefiller.hh>
//! make the local view function available whenever we use this class
@@ -61,13 +62,28 @@ class CCMpfaGridFluxVariablesCache<TypeTag, true>
using FluxVariablesCache = typename GET_PROP_TYPE(TypeTag, FluxVariablesCache);
using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
+ using FluxVariablesCacheFiller = CCMpfaFluxVariablesCacheFiller<TypeTag>;
using PrimaryInteractionVolume = typename GET_PROP_TYPE(TypeTag, PrimaryInteractionVolume);
- using PrimaryIvDataHandle = typename PrimaryInteractionVolume::Traits::DataHandle;
+ using PrimaryMatVecTraits = typename PrimaryInteractionVolume::Traits::MatVecTraits;
using SecondaryInteractionVolume = typename GET_PROP_TYPE(TypeTag, SecondaryInteractionVolume);
- using SecondaryIvDataHandle = typename SecondaryInteractionVolume::Traits::DataHandle;
- using FluxVariablesCacheFiller = CCMpfaFluxVariablesCacheFiller<TypeTag>;
+ using SecondaryMatVecTraits = typename SecondaryInteractionVolume::Traits::MatVecTraits;
+
+ //! physics traits class to define the data handles
+ struct PhysicsTraits
+ {
+ static constexpr bool enableAdvection = GET_PROP_VALUE(TypeTag, EnableAdvection);
+ static constexpr bool enableMolecularDiffusion = GET_PROP_VALUE(TypeTag, EnableMolecularDiffusion);
+ static constexpr bool enableHeatConduction = GET_PROP_VALUE(TypeTag, EnableEnergyBalance);
+
+ static constexpr int numPhases = GET_PROP_VALUE(TypeTag, NumPhases);
+ static constexpr int numComponents = GET_PROP_VALUE(TypeTag, NumComponents);
+ };
public:
+ // export the data handle types used
+ using PrimaryIvDataHandle = InteractionVolumeDataHandle< PrimaryMatVecTraits, PhysicsTraits >;
+ using SecondaryIvDataHandle = InteractionVolumeDataHandle< SecondaryMatVecTraits, PhysicsTraits >;
+
//! export the type of the local view
using LocalView = typename GET_PROP_TYPE(TypeTag, ElementFluxVariablesCache);
@@ -254,6 +270,10 @@ public:
//! export the type of the local view
using LocalView = typename GET_PROP_TYPE(TypeTag, ElementFluxVariablesCache);
+ //! export the data handle types used by the local view
+ using PrimaryIvDataHandle = typename LocalView::PrimaryIvDataHandle;
+ using SecondaryIvDataHandle = typename LocalView::SecondaryIvDataHandle;
+
//! The constructor
CCMpfaGridFluxVariablesCache(const Problem& problem) : problemPtr_(&problem) {}
diff --git a/dumux/discretization/cellcentered/mpfa/gridinteractionvolumeindexsets.hh b/dumux/discretization/cellcentered/mpfa/gridinteractionvolumeindexsets.hh
index 7a6875ca5c1720aaaec215710a02812c877a4982..3f7020bc6a4dae2334d2497283ebf1b960b3c081 100644
--- a/dumux/discretization/cellcentered/mpfa/gridinteractionvolumeindexsets.hh
+++ b/dumux/discretization/cellcentered/mpfa/gridinteractionvolumeindexsets.hh
@@ -25,7 +25,6 @@
#define DUMUX_DISCRETIZATION_MPFA_O_GRIDINTERACTIONVOLUME_INDEXSETS_HH
#include <memory>
-#include <dumux/common/properties.hh>
#include "dualgridindexset.hh"
@@ -78,9 +77,9 @@ public:
{
const auto vIdxGlobal = fvGridGeometry.vertexMapper().index(vertex);
if (!fvGridGeometry.vertexUsesSecondaryInteractionVolume(vIdxGlobal))
- numPrimaryIV_ += PrimaryInteractionVolume::numInteractionVolumesAtVertex((*dualGridIndexSet_)[vIdxGlobal]);
+ numPrimaryIV_ += PrimaryInteractionVolume::numIVAtVertex((*dualGridIndexSet_)[vIdxGlobal]);
else
- numSecondaryIV_ += SecondaryInteractionVolume::numInteractionVolumesAtVertex((*dualGridIndexSet_)[vIdxGlobal]);
+ numSecondaryIV_ += SecondaryInteractionVolume::numIVAtVertex((*dualGridIndexSet_)[vIdxGlobal]);
}
// reserve memory
@@ -93,9 +92,15 @@ public:
{
const auto vIdxGlobal = fvGridGeometry.vertexMapper().index(vertex);
if (!fvGridGeometry.vertexUsesSecondaryInteractionVolume(vIdxGlobal))
- PrimaryInteractionVolume::addInteractionVolumeIndexSets(primaryIVIndexSets_, scvfIndexMap_, (*dualGridIndexSet_)[vIdxGlobal]);
+ PrimaryInteractionVolume::addIVIndexSets(primaryIVIndexSets_,
+ scvfIndexMap_,
+ (*dualGridIndexSet_)[vIdxGlobal],
+ fvGridGeometry.flipScvfIndexSet());
else
- SecondaryInteractionVolume::addInteractionVolumeIndexSets(secondaryIVIndexSets_, scvfIndexMap_, (*dualGridIndexSet_)[vIdxGlobal]);
+ SecondaryInteractionVolume::addIVIndexSets(secondaryIVIndexSets_,
+ scvfIndexMap_,
+ (*dualGridIndexSet_)[vIdxGlobal],
+ fvGridGeometry.flipScvfIndexSet());
}
}
diff --git a/dumux/discretization/cellcentered/mpfa/helper.hh b/dumux/discretization/cellcentered/mpfa/helper.hh
index 71aa0263bae23999892d91bbc0cfdf94446df8f3..dac6e546245690ca97c9da50c658cfe2050b0581 100644
--- a/dumux/discretization/cellcentered/mpfa/helper.hh
+++ b/dumux/discretization/cellcentered/mpfa/helper.hh
@@ -32,7 +32,6 @@
#include <dune/geometry/type.hh>
#include <dune/geometry/referenceelements.hh>
-#include <dumux/common/properties.hh>
#include <dumux/common/math.hh>
namespace Dumux {
@@ -41,34 +40,22 @@ namespace Dumux {
* \ingroup CCMpfaDiscretization
* \brief Dimension-specific helper class to get data required for mpfa scheme.
*/
-template<class TypeTag, int dim, int dimWorld>
+template<class FVGridGeometry, int dim, int dimWorld>
class MpfaDimensionHelper;
/*!
* \ingroup CCMpfaDiscretization
- * \brief Dimension-specific helper class to get data required for mpfa scheme.
- * Specialization for dim == 2 & dimWorld == 2
+ * \brief Dimension-specific mpfa helper class for dim == 2 & dimWorld == 2
*/
-template<class TypeTag>
-class MpfaDimensionHelper<TypeTag, /*dim*/2, /*dimWorld*/2>
+template<class FVGridGeometry>
+class MpfaDimensionHelper<FVGridGeometry, /*dim*/2, /*dimWorld*/2>
{
- using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
- using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
- using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry)::LocalView;
- using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
+ using GridView = typename FVGridGeometry::GridView;
+ using CoordScalar = typename GridView::ctype;
+ using GlobalPosition = Dune::FieldVector<CoordScalar, GridView::dimensionworld>;
+
+ using SubControlVolumeFace = typename FVGridGeometry::SubControlVolumeFace;
using ScvfCornerVector = typename SubControlVolumeFace::Traits::CornerStorage;
- using InteractionVolume = typename GET_PROP_TYPE(TypeTag, PrimaryInteractionVolume);
- using LocalScvType = typename InteractionVolume::Traits::LocalScvType;
- using ScvBasis = typename LocalScvType::LocalBasis;
-
- // We know that dim = 2 & dimworld = 2. However, the specialization for
- // dim = 2 & dimWorld = 3 reuses some methods of this class, thus, we
- // obtain the world dimension from the grid view here to get the
- // GlobalPosition vector right. Be picky about the dimension though.
- static_assert(GridView::dimension == 2, "The chosen mpfa helper expects a grid view with dim = 2");
- static const int dim = 2;
- static const int dimWorld = GridView::dimensionworld;
- using GlobalPosition = Dune::FieldVector<Scalar, dimWorld>;
// Container to store the positions of intersections required for scvf
// corner computation. In 2d, these are the center plus the two corners
@@ -80,9 +67,10 @@ public:
*
* \param scvBasis The basis of an scv
*/
+ template< class ScvBasis >
static ScvBasis calculateInnerNormals(const ScvBasis& scvBasis)
{
- static const Dune::FieldMatrix<Scalar, dim, dim> R = {{0.0, 1.0}, {-1.0, 0.0}};
+ static const Dune::FieldMatrix<CoordScalar, 2, 2> R = {{0.0, 1.0}, {-1.0, 0.0}};
ScvBasis innerNormals;
R.mv(scvBasis[1], innerNormals[0]);
@@ -103,10 +91,11 @@ public:
*
* \param scvBasis The basis of an scv
*/
- static typename LocalScvType::ctype calculateDetX(const ScvBasis& scvBasis)
+ template< class ScvBasis >
+ static CoordScalar calculateDetX(const ScvBasis& scvBasis)
{
using std::abs;
- return abs(crossProduct<Scalar>(scvBasis[0], scvBasis[1]));
+ return abs(crossProduct<CoordScalar>(scvBasis[0], scvBasis[1]));
}
/*!
@@ -139,8 +128,9 @@ public:
*
* \param scvBasis The basis of an scv
*/
+ template< class ScvBasis >
static bool isRightHandSystem(const ScvBasis& scvBasis)
- { return !std::signbit(crossProduct<Scalar>(scvBasis[0], scvBasis[1])); }
+ { return !std::signbit(crossProduct<CoordScalar>(scvBasis[0], scvBasis[1])); }
/*!
* \brief Returns a vector containing the positions on a given intersection
@@ -164,7 +154,8 @@ public:
p[0] = 0.0;
for (unsigned int c = 0; c < numCorners; ++c)
{
- p[c+1] = eg.global(refElement.position(refElement.subEntity(indexInElement, 1, c, dim), dim));
+ // codim = 1, dim = 2
+ p[c+1] = eg.global(refElement.position(refElement.subEntity(indexInElement, 1, c, 2), 2));
p[0] += p[c+1];
}
p[0] /= numCorners;
@@ -199,7 +190,7 @@ public:
*
* \param scvfCorners Container with the corners of the scvf
*/
- static Scalar getScvfArea(const ScvfCornerVector& scvfCorners)
+ static CoordScalar getScvfArea(const ScvfCornerVector& scvfCorners)
{ return (scvfCorners[1]-scvfCorners[0]).two_norm(); }
/*!
@@ -229,19 +220,15 @@ public:
/*!
* \ingroup CCMpfaDiscretization
- * \brief Dimension-specific helper class to get data required for mpfa scheme.
- * Specialization for dim == 2 & dimWorld == 2. Reuses some functionality
- * of the specialization for dim = dimWorld = 2
+ * \brief Dimension-specific mpfa helper class for dim == 2 & dimWorld == 2.
+ * Reuses some functionality of the specialization for dim = dimWorld = 2
*/
-template<class TypeTag>
-class MpfaDimensionHelper<TypeTag, /*dim*/2, /*dimWorld*/3>
- : public MpfaDimensionHelper<TypeTag, 2, 2>
+template<class FVGridGeometry>
+class MpfaDimensionHelper<FVGridGeometry, /*dim*/2, /*dimWorld*/3>
+ : public MpfaDimensionHelper<FVGridGeometry, 2, 2>
{
- using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
- using InteractionVolume = typename GET_PROP_TYPE(TypeTag, PrimaryInteractionVolume);
- using LocalScvType = typename InteractionVolume::Traits::LocalScvType;
- using ScvBasis = typename LocalScvType::LocalBasis;
-
+ using GridView = typename FVGridGeometry::GridView;
+ using CoordScalar = typename GridView::ctype;
public:
/*!
@@ -249,19 +236,20 @@ public:
*
* \param scvBasis The basis of an scv
*/
+ template< class ScvBasis >
static ScvBasis calculateInnerNormals(const ScvBasis& scvBasis)
{
// compute vector normal to the basis plane
const auto normal = [&] () {
- auto n = crossProduct<Scalar>(scvBasis[0], scvBasis[1]);
+ auto n = crossProduct<CoordScalar>(scvBasis[0], scvBasis[1]);
n /= n.two_norm();
return n;
} ();
// compute inner normals using the normal vector
ScvBasis innerNormals;
- innerNormals[0] = crossProduct<Scalar>(scvBasis[1], normal);
- innerNormals[1] = crossProduct<Scalar>(normal, scvBasis[0]);
+ innerNormals[0] = crossProduct<CoordScalar>(scvBasis[1], normal);
+ innerNormals[1] = crossProduct<CoordScalar>(normal, scvBasis[0]);
return innerNormals;
}
@@ -274,10 +262,11 @@ public:
*
* \param scvBasis The basis of an scv
*/
- static typename LocalScvType::ctype calculateDetX(const ScvBasis& scvBasis)
+ template< class ScvBasis >
+ static CoordScalar calculateDetX(const ScvBasis& scvBasis)
{
using std::abs;
- return abs(crossProduct<Scalar>(scvBasis[0], scvBasis[1]).two_norm());
+ return abs(crossProduct<CoordScalar>(scvBasis[0], scvBasis[1]).two_norm());
}
/*!
@@ -287,32 +276,24 @@ public:
*
* \param scvBasis The basis of an scv
*/
- static bool isRightHandSystem(const ScvBasis& scvBasis)
- { return true; }
+ template< class ScvBasis >
+ static constexpr bool isRightHandSystem(const ScvBasis& scvBasis) { return true; }
};
/*!
* \ingroup CCMpfaDiscretization
- * \brief Dimension-specific helper class to get data required for mpfa scheme.
- * Specialization for dim == 3 & dimWorld == 3.
+ * \brief Dimension-specific mpfa helper class for dim == 3 & dimWorld == 3.
+ *
*/
-template<class TypeTag>
-class MpfaDimensionHelper<TypeTag, /*dim*/3, /*dimWorld*/3>
+template<class FVGridGeometry>
+class MpfaDimensionHelper<FVGridGeometry, /*dim*/3, /*dimWorld*/3>
{
- using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
- using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
- using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry)::LocalView;
- using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
+ using SubControlVolumeFace = typename FVGridGeometry::SubControlVolumeFace;
using ScvfCornerVector = typename SubControlVolumeFace::Traits::CornerStorage;
- using InteractionVolume = typename GET_PROP_TYPE(TypeTag, PrimaryInteractionVolume);
- using LocalScvType = typename InteractionVolume::Traits::LocalScvType;
- using ScvBasis = typename LocalScvType::LocalBasis;
// Be picky about the dimensions
- static_assert(GridView::dimension == 3 && GridView::dimensionworld == 3,
- "The chosen mpfa helper expects a grid view with dim = 3 & dimWorld = 3");
- static const int dim = 3;
- static const int dimWorld = 3;
- using GlobalPosition = Dune::FieldVector<Scalar, dimWorld>;
+ using GridView = typename FVGridGeometry::GridView;
+ using CoordScalar = typename GridView::ctype;
+ using GlobalPosition = Dune::FieldVector<CoordScalar, /*dimworld*/3>;
// container to store the positions of intersections required for
// scvf corner computation. Maximum number of points needed is 9
@@ -326,13 +307,14 @@ public:
*
* \param scvBasis The basis of an scv
*/
+ template< class ScvBasis >
static ScvBasis calculateInnerNormals(const ScvBasis& scvBasis)
{
ScvBasis innerNormals;
- innerNormals[0] = crossProduct<Scalar>(scvBasis[1], scvBasis[2]);
- innerNormals[1] = crossProduct<Scalar>(scvBasis[2], scvBasis[0]);
- innerNormals[2] = crossProduct<Scalar>(scvBasis[0], scvBasis[1]);
+ innerNormals[0] = crossProduct<CoordScalar>(scvBasis[1], scvBasis[2]);
+ innerNormals[1] = crossProduct<CoordScalar>(scvBasis[2], scvBasis[0]);
+ innerNormals[2] = crossProduct<CoordScalar>(scvBasis[0], scvBasis[1]);
if (!isRightHandSystem(scvBasis))
{
@@ -350,10 +332,11 @@ public:
*
* \param scvBasis The basis of an scv
*/
- static typename LocalScvType::ctype calculateDetX(const ScvBasis& scvBasis)
+ template< class ScvBasis >
+ static CoordScalar calculateDetX(const ScvBasis& scvBasis)
{
using std::abs;
- return abs(tripleProduct<Scalar>(scvBasis[0], scvBasis[1], scvBasis[2]));
+ return abs(tripleProduct<CoordScalar>(scvBasis[0], scvBasis[1], scvBasis[2]));
}
/*!
@@ -396,8 +379,9 @@ public:
*
* \param scvBasis The basis of an scv
*/
+ template< class ScvBasis >
static bool isRightHandSystem(const ScvBasis& scvBasis)
- { return !std::signbit(tripleProduct<Scalar>(scvBasis[0], scvBasis[1], scvBasis[2])); }
+ { return !std::signbit(tripleProduct<CoordScalar>(scvBasis[0], scvBasis[1], scvBasis[2])); }
/*!
* \brief Returns a vector containing the positions on a given intersection
@@ -424,7 +408,8 @@ public:
p[0] = 0.0;
for (unsigned int c = 0; c < numCorners; ++c)
{
- p[c+1] = eg.global(refElement.position(refElement.subEntity(indexInElement, 1, c, dim), dim));
+ // codim = 1, dim = 3
+ p[c+1] = eg.global(refElement.position(refElement.subEntity(indexInElement, 1, c, 3), 3));
p[0] += p[c+1];
}
p[0] /= numCorners;
@@ -441,6 +426,7 @@ public:
p[numCorners+2] /= 2;
p[numCorners+3] = p[3] + p[2];
p[numCorners+3] /= 2;
+ return p;
}
case 4: // quadrilateral
{
@@ -453,14 +439,12 @@ public:
p[numCorners+3] /= 2;
p[numCorners+4] = p[4] + p[3];
p[numCorners+4] /= 2;
+ return p;
}
default:
- DUNE_THROW(Dune::NotImplemented, "Mpfa scvf corners for dim = " << dim
- << ", dimWorld = " << dimWorld
- << ", corners = " << numCorners);
+ DUNE_THROW(Dune::NotImplemented,
+ "Mpfa scvf corners for dim = 3, dimWorld = 3, corners = " << numCorners);
}
-
- return p;
}
/*!
@@ -515,9 +499,8 @@ public:
p[map[cornerIdx][3]]} );
}
default:
- DUNE_THROW(Dune::NotImplemented, "Mpfa scvf corners for dim = " << dim
- << ", dimWorld = " << dimWorld
- << ", corners = " << numIntersectionCorners);
+ DUNE_THROW(Dune::NotImplemented,
+ "Mpfa scvf corners for dim = 3, dimWorld = 3, corners = " << numIntersectionCorners);
}
}
@@ -526,7 +509,7 @@ public:
*
* \param scvfCorners Container with the corners of the scvf
*/
- static Scalar getScvfArea(const ScvfCornerVector& scvfCorners)
+ static CoordScalar getScvfArea(const ScvfCornerVector& scvfCorners)
{
// after Wolfram alpha quadrilateral area
return 0.5*Dumux::crossProduct(scvfCorners[3]-scvfCorners[0], scvfCorners[2]-scvfCorners[1]).two_norm();
@@ -566,26 +549,26 @@ public:
/*!
* \ingroup CCMpfaDiscretization
* \brief Helper class to get the required information on an interaction volume.
- * Specializations depending on the method and dimension are provided.
+ *
+ * \tparam FVGridGeometry The finite volume grid geometry
*/
-template<class TypeTag, int dim, int dimWorld>
-class CCMpfaHelperImplementation : public MpfaDimensionHelper<TypeTag, dim, dimWorld>
+template<class FVGridGeometry>
+class CCMpfaHelper : public MpfaDimensionHelper<FVGridGeometry,
+ FVGridGeometry::GridView::dimension,
+ FVGridGeometry::GridView::dimensionworld>
{
- using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
-
- using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
- using Element = typename GridView::template Codim<0>::Entity;
- using IndexType = typename GridView::IndexSet::IndexType;
-
- using PrimaryInteractionVolume = typename GET_PROP_TYPE(TypeTag, PrimaryInteractionVolume);
- using SecondaryInteractionVolume = typename GET_PROP_TYPE(TypeTag, SecondaryInteractionVolume);
+ using PrimaryInteractionVolume = typename FVGridGeometry::GridIVIndexSets::PrimaryInteractionVolume;
+ using SecondaryInteractionVolume = typename FVGridGeometry::GridIVIndexSets::SecondaryInteractionVolume;
- using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
using VertexMapper = typename FVGridGeometry::VertexMapper;
using FVElementGeometry = typename FVGridGeometry::LocalView;
using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
using ScvfCornerVector = typename SubControlVolumeFace::Traits::CornerStorage;
+ using GridView = typename FVGridGeometry::GridView;
+ static constexpr int dim = GridView::dimension;
+ static constexpr int dimWorld = GridView::dimensionworld;
+
using CoordScalar = typename GridView::ctype;
using GlobalPosition = Dune::FieldVector<CoordScalar, dimWorld>;
using ReferenceElements = typename Dune::ReferenceElements<CoordScalar, dim>;
@@ -597,14 +580,17 @@ public:
* \param scvfCorners Container with the corners of the scvf
* \param q Parameterization of the integration point on the scvf
*/
+ template< class Scalar >
static GlobalPosition getScvfIntegrationPoint(const ScvfCornerVector& scvfCorners, Scalar q)
{
- // scvfs in 3d are always quadrilaterals
// ordering -> first corner: facet center, last corner: vertex
if (q == 0.0)
return scvfCorners[0];
- const auto d = [&] () { auto tmp = scvfCorners.back() - scvfCorners.front(); tmp *= q; return tmp; } ();
- return scvfCorners[0] + d;
+
+ auto ip = scvfCorners.back() - scvfCorners.front();
+ ip *= q;
+ ip += scvfCorners[0];
+ return ip;
}
// returns a vector which maps true to each vertex on processor boundaries and false otherwise
@@ -617,7 +603,7 @@ public:
if (Dune::MPIHelper::getCollectiveCommunication().size() == 1)
return ghostVertices;
- // mpfa methods can not yet handle ghost cells
+ // mpfa methods cannot yet handle ghost cells
if (gridView.ghostSize(0) > 0)
DUNE_THROW(Dune::InvalidStateException, "Mpfa methods in parallel do not work with ghost cells. Use overlap cells instead.");
@@ -660,16 +646,6 @@ public:
{ return std::find(vector.begin(), vector.end(), value) != vector.end(); }
};
-/*!
- * \ingroup CCMpfaDiscretization
- * \brief Helper class to obtain data required for mpfa methods.
- * It inherits from dimension-, dimensionworld- and method-specific implementations.
- */
-template<class TypeTag>
-using CCMpfaHelper = CCMpfaHelperImplementation<TypeTag,
- GET_PROP_TYPE(TypeTag, GridView)::dimension,
- GET_PROP_TYPE(TypeTag, GridView)::dimensionworld>;
-
} // end namespace Dumux
#endif
diff --git a/dumux/discretization/cellcentered/mpfa/interactionvolumebase.hh b/dumux/discretization/cellcentered/mpfa/interactionvolumebase.hh
index d5e46c9b0ff4a71744be6ffe45a8dd158bcac447..95bee55df5e2d5f922aee9ddd15164bb7d4fd4c1 100644
--- a/dumux/discretization/cellcentered/mpfa/interactionvolumebase.hh
+++ b/dumux/discretization/cellcentered/mpfa/interactionvolumebase.hh
@@ -26,8 +26,6 @@
#include <dune/common/exceptions.hh>
-#include <dumux/common/properties.hh>
-
namespace Dumux
{
@@ -39,21 +37,9 @@ namespace Dumux
* The traits should contain the following type definitions:
*
* \code
- * //! export the problem type (needed for iv-local assembly)
- * using Problem = ...;
- * //! export the type of the local view on the finite volume grid geometry
- * using FVElementGeometry = ...;
- * //! export the type of the local view on the grid volume variables
- * using ElementVolumeVariables = ...;
* //! export the type of grid view
* using GridView = ...;
- * //! export the type used for scalar values
- * using ScalarType = ...;
- * //! export the type of the mpfa helper class
- * using MpfaHelper = ...;
* //! export the type used for local indices
- * using LocalIndexType = ...;
- * //! export the type for the interaction volume index set
* using IndexSet = ...;
* //! export the type of interaction-volume local scvs
* using LocalScvType = ...;
@@ -61,12 +47,8 @@ namespace Dumux
* using LocalScvfType = ...;
* //! export the type of used for the iv-local face data
* using LocalFaceData = ...;
- * //! export the type used for iv-local matrices
- * using Matrix = ...;
- * //! export the type used for iv-local vectors
- * using Vector = ...;
- * //! export the data handle type for this iv
- * using DataHandle = ...;
+ * //! export the matrix/vector type traits to be used by the iv
+ * using MatVecTraits = ...;
* \endcode
*/
@@ -78,25 +60,27 @@ namespace Dumux
* \tparam Impl The actual implementation of the interaction volume
* \tparam T The traits class to be used
*/
-template< class Impl, class T>
+template< class Impl, class T >
class CCMpfaInteractionVolumeBase
{
// Curiously recurring template pattern
Impl& asImp() { return static_cast<Impl&>(*this); }
const Impl& asImp() const { return static_cast<const Impl&>(*this); }
- using Problem = typename T::Problem;
- using FVElementGeometry = typename T::FVElementGeometry;
- using ElementVolumeVariables = typename T::ElementVolumeVariables;
-
using GridView = typename T::GridView;
using Element = typename GridView::template Codim<0>::Entity;
+ using NodalStencilType = typename T::IndexSet::NodalGridStencilType;
+ using LocalIndexType = typename T::IndexSet::LocalIndexType;
+ using LocalScvType = typename T::LocalScvType;
+ using LocalScvfType = typename T::LocalScvfType;
+
public:
//! state the traits type publicly
using Traits = T;
//! Prepares everything for the assembly
+ template< class Problem, class FVElementGeometry >
void setUpLocalScope(const typename Traits::IndexSet& indexSet,
const Problem& problem,
const FVElementGeometry& fvGeometry)
@@ -121,32 +105,32 @@ public:
{ DUNE_THROW(Dune::NotImplemented, "Interaction volume implementation does not provide a localFaceData() funtion"); }
//! returns the cell-stencil of this interaction volume
- const typename Traits::IndexSet::GridStencilType& stencil() const { return asImp().stencil(); }
+ const NodalStencilType& stencil() const { return asImp().stencil(); }
//! returns the local scvf entity corresponding to a given iv-local scvf idx
- const typename Traits::LocalScvfType& localScvf(typename Traits::LocalIndexType ivLocalScvfIdx) const
- { return asImp().localScvf(ivLocalScvfIdx); }
+ const LocalScvfType& localScvf(LocalIndexType ivLocalScvfIdx) const { return asImp().localScvf(ivLocalScvfIdx); }
//! returns the local scv entity corresponding to a given iv-local scv idx
- const typename Traits::LocalScvType& localScv(typename Traits::LocalIndexType ivLocalScvIdx) const
- { return asImp().localScv(ivLocalScvIdx); }
+ const LocalScvType& localScv(LocalIndexType ivLocalScvIdx) const { return asImp().localScv(ivLocalScvIdx); }
//! returns the element in which the scv with the given local idx is embedded in
- const Element& element(typename Traits::LocalIndexType ivLocalScvIdx) const
- { return asImp().element(); }
+ const Element& element(LocalIndexType ivLocalScvIdx) const { return asImp().element(); }
//! returns the number of interaction volumes living around a vertex
template< class NodalIndexSet >
- static std::size_t numInteractionVolumesAtVertex(const NodalIndexSet& nodalIndexSet)
- { return Impl::numInteractionVolumesAtVertex(nodalIndexSet); }
+ static std::size_t numIVAtVertex(const NodalIndexSet& nodalIndexSet) { return Impl::numIVAtVertex(nodalIndexSet); }
//! adds the iv index sets living around a vertex to a given container
//! and stores the the corresponding index in a map for each scvf
- template<class IvIndexSetContainer, class ScvfIndexMap, class NodalIndexSet>
- static void addInteractionVolumeIndexSets(IvIndexSetContainer& ivIndexSetContainer,
- ScvfIndexMap& scvfIndexMap,
- const NodalIndexSet& nodalIndexSet)
- { Impl::addInteractionVolumeIndexSets(ivIndexSetContainer, scvfIndexMap, nodalIndexSet); }
+ template< class IvIndexSetContainer,
+ class ScvfIndexMap,
+ class NodalIndexSet,
+ class FlipScvfIndexSet >
+ static void addIVIndexSets(IvIndexSetContainer& ivIndexSetContainer,
+ ScvfIndexMap& scvfIndexMap,
+ const NodalIndexSet& nodalIndexSet,
+ const FlipScvfIndexSet& flipScvfIndexSet)
+ { Impl::addIVIndexSets(ivIndexSetContainer, scvfIndexMap, nodalIndexSet, flipScvfIndexSet); }
};
} // end namespace Dumux
diff --git a/dumux/discretization/cellcentered/mpfa/interactionvolumedatahandle.hh b/dumux/discretization/cellcentered/mpfa/interactionvolumedatahandle.hh
index 7a41428485d51205c0db03158f21ed971346ac07..72e2ddc632184a61f1083d35e6d09501218890b4 100644
--- a/dumux/discretization/cellcentered/mpfa/interactionvolumedatahandle.hh
+++ b/dumux/discretization/cellcentered/mpfa/interactionvolumedatahandle.hh
@@ -26,7 +26,8 @@
#include <vector>
-#include <dumux/common/properties.hh>
+#include <dune/common/dynvector.hh>
+
#include <dumux/common/parameters.hh>
#include <dumux/common/matrixvectorhelper.hh>
@@ -42,26 +43,26 @@ public:
};
//! Data handle for quantities related to advection
-template<class TypeTag, class M, class V, class LI, bool EnableAdvection>
+template<class MatVecTraits, class PhysicsTraits, bool EnableAdvection>
class AdvectionDataHandle
{
// export matrix & vector types from interaction volume
- using Matrix = M;
- using Vector = V;
- using LocalIndexType = LI;
- using Scalar = typename Vector::value_type;
-
- using OutsideDataContainer = std::vector< std::vector<Vector> >;
+ using AMatrix = typename MatVecTraits::AMatrix;
+ using CMatrix = typename MatVecTraits::CMatrix;
+ using TMatrix = typename MatVecTraits::TMatrix;
+ using CellVector = typename MatVecTraits::CellVector;
+ using FaceVector = typename MatVecTraits::FaceVector;
+ using FaceScalar = typename FaceVector::value_type;
+ using OutsideGravityStorage = std::vector< Dune::DynamicVector<FaceScalar> >;
- static constexpr int dim = GET_PROP_TYPE(TypeTag, GridView)::dimension;
- static constexpr int dimWorld = GET_PROP_TYPE(TypeTag, GridView)::dimensionworld;
- static constexpr int numPhases = GET_PROP_VALUE(TypeTag, NumPhases);
+ static constexpr int numPhases = PhysicsTraits::numPhases;
public:
//! prepare data handle for subsequent fill (normal grids)
- template< class InteractionVolume, int d = dim, int dw = dimWorld, std::enable_if_t<(d==dw), int> = 0>
- void resize(const InteractionVolume& iv)
+ template< class IV,
+ std::enable_if_t<(IV::Traits::GridView::dimension==IV::Traits::GridView::dimensionworld), int> = 0>
+ void resize(const IV& iv)
{
// resize transmissibility matrix & pressure vectors
resizeMatrix(T_, iv.numFaces(), iv.numKnowns());
@@ -69,7 +70,7 @@ public:
resizeVector(p_[pIdx], iv.numKnowns());
// maybe resize gravity container
- static const bool enableGravity = getParamFromGroup<bool>(GET_PROP_VALUE(TypeTag, ModelParameterGroup), "Problem.EnableGravity");
+ static const bool enableGravity = getParam<bool>("Problem.EnableGravity");
if (enableGravity)
{
resizeMatrix(CA_, iv.numFaces(), iv.numUnknowns());
@@ -80,11 +81,13 @@ public:
//! prepare data handle for subsequent fill (surface grids)
- template< class InteractionVolume, int d = dim, int dw = dimWorld, std::enable_if_t<(d<dw), int> = 0>
- void resize(const InteractionVolume& iv)
+ template< class IV,
+ std::enable_if_t<(IV::Traits::GridView::dimension<IV::Traits::GridView::dimensionworld), int> = 0>
+ void resize(const IV& iv)
{
- static const bool enableGravity = getParamFromGroup<bool>(GET_PROP_VALUE(TypeTag, ModelParameterGroup), "Problem.EnableGravity");
+ using LocalIndexType = typename IV::Traits::IndexSet::LocalIndexType;
+ static const bool enableGravity = getParam<bool>("Problem.EnableGravity");
if (!enableGravity)
{
resizeMatrix(T_, iv.numFaces(), iv.numKnowns());
@@ -129,66 +132,61 @@ public:
}
//! Access to the iv-wide pressure of one phase
- const Vector& pressures(unsigned int pIdx) const { return p_[pIdx]; }
- Vector& pressures(unsigned int pIdx) { return p_[pIdx]; }
+ const CellVector& pressures(unsigned int pIdx) const { return p_[pIdx]; }
+ CellVector& pressures(unsigned int pIdx) { return p_[pIdx]; }
//! Access to the gravitational flux contributions for one phase
- const Vector& gravity(unsigned int pIdx) const { return g_[pIdx]; }
- Vector& gravity(unsigned int pIdx) { return g_[pIdx]; }
+ const FaceVector& gravity(unsigned int pIdx) const { return g_[pIdx]; }
+ FaceVector& gravity(unsigned int pIdx) { return g_[pIdx]; }
//! Access to the gravitational flux contributions for all phases
- const std::array< Vector, numPhases >& gravity() const { return g_; }
- std::array< Vector, numPhases >& gravity() { return g_; }
+ const std::array< FaceVector, numPhases >& gravity() const { return g_; }
+ std::array< FaceVector, numPhases >& gravity() { return g_; }
//! Projection matrix for gravitational acceleration
- const Matrix& advectionCA() const { return CA_; }
- Matrix& advectionCA() { return CA_; }
+ const CMatrix& advectionCA() const { return CA_; }
+ CMatrix& advectionCA() { return CA_; }
//! Additional projection matrix needed on surface grids
- const Matrix& advectionA() const { return A_; }
- Matrix& advectionA() { return A_; }
+ const AMatrix& advectionA() const { return A_; }
+ AMatrix& advectionA() { return A_; }
//! The transmissibilities associated with advective fluxes
- const Matrix& advectionT() const { return T_; }
- Matrix& advectionT() { return T_; }
+ const TMatrix& advectionT() const { return T_; }
+ TMatrix& advectionT() { return T_; }
//! The transmissibilities for "outside" faces (used on surface grids)
- const std::vector< std::vector<Vector> >& advectionTout() const { return outsideT_; }
- std::vector< std::vector<Vector> >& advectionTout() { return outsideT_; }
+ const std::vector< std::vector<CellVector> >& advectionTout() const { return outsideT_; }
+ std::vector< std::vector<CellVector> >& advectionTout() { return outsideT_; }
//! The gravitational acceleration for "outside" faces (used on surface grids)
- const std::array< std::vector<Vector>, numPhases >& gravityOutside() const { return outsideG_; }
- std::array< std::vector<Vector>, numPhases >& gravityOutside() { return outsideG_; }
+ const std::array< OutsideGravityStorage, numPhases >& gravityOutside() const { return outsideG_; }
+ std::array< OutsideGravityStorage, numPhases >& gravityOutside() { return outsideG_; }
//! The gravitational acceleration for one phase on "outside" faces (used on surface grids)
- const std::vector<Vector>& gravityOutside(unsigned int pIdx) const { return outsideG_[pIdx]; }
- std::vector<Vector>& gravityOutside(unsigned int pIdx) { return outsideG_[pIdx]; }
+ const OutsideGravityStorage& gravityOutside(unsigned int pIdx) const { return outsideG_[pIdx]; }
+ OutsideGravityStorage& gravityOutside(unsigned int pIdx) { return outsideG_[pIdx]; }
private:
- //! advection-related variables
- Matrix T_; //!< The transmissibilities such that f_i = T_ij*p_j
- Matrix CA_; //!< Matrix to project gravitational acceleration to all scvfs
- Matrix A_; //!< Matrix additionally needed for the projection on surface grids
- std::array< Vector, numPhases > p_; //!< The interaction volume-wide phase pressures
- std::array< Vector, numPhases > g_; //!< The gravitational acceleration at each scvf (only for enabled gravity)
- std::vector< std::vector<Vector> > outsideT_; //!< The transmissibilities for "outside" faces (only on surface grids)
- std::array< std::vector<Vector>, numPhases > outsideG_; //!< The gravitational acceleration on "outside" faces (only on surface grids)
+ TMatrix T_; //!< The transmissibilities such that f_i = T_ij*p_j
+ CMatrix CA_; //!< Matrix to project gravitational acceleration to all scvfs
+ AMatrix A_; //!< Matrix additionally needed for the projection on surface grids
+ std::array< CellVector, numPhases > p_; //!< The interaction volume-wide phase pressures
+ std::array< FaceVector, numPhases > g_; //!< The gravitational acceleration at each scvf (only for enabled gravity)
+ std::vector< std::vector<CellVector> > outsideT_; //!< The transmissibilities for "outside" faces (only on surface grids)
+ std::array< OutsideGravityStorage, numPhases > outsideG_; //!< The gravitational acceleration on "outside" faces (only on surface grids)
};
//! Data handle for quantities related to diffusion
-template<class TypeTag, class M, class V, class LI, bool EnableDiffusion>
+template<class MatVecTraits, class PhysicsTraits, bool EnableDiffusion>
class DiffusionDataHandle
{
- // export matrix & vector types from interaction volume
- using Matrix = M;
- using Vector = V;
- using OutsideTContainer = std::vector< std::vector<Vector> >;
-
- static constexpr int dim = GET_PROP_TYPE(TypeTag, GridView)::dimension;
- static constexpr int dimWorld = GET_PROP_TYPE(TypeTag, GridView)::dimensionworld;
+ using TMatrix = typename MatVecTraits::TMatrix;
+ using CellVector = typename MatVecTraits::CellVector;
+ using OutsideTContainer = std::vector< std::vector<CellVector> >;
- static constexpr int numPhases = GET_PROP_VALUE(TypeTag, NumPhases);
- static constexpr int numComponents = GET_PROP_VALUE(TypeTag, NumComponents);
+ static constexpr int numPhases = PhysicsTraits::numPhases;
+ static constexpr int numComponents = PhysicsTraits::numComponents;
public:
//! diffusion caches need to set phase and component index
@@ -203,19 +201,17 @@ public:
{
for (unsigned int cIdx = 0; cIdx < numComponents; ++cIdx)
{
- if (pIdx == cIdx)
- continue;
-
// resize transmissibility matrix & mole fraction vector
resizeMatrix(T_[pIdx][cIdx], iv.numFaces(), iv.numKnowns());
resizeVector(xj_[pIdx][cIdx], iv.numKnowns());
// resize outsideTij on surface grids
- if (dim < dimWorld)
+ using GridView = typename InteractionVolume::Traits::GridView;
+ if (int(GridView::dimension) < int(GridView::dimensionworld))
{
outsideT_[pIdx][cIdx].resize(iv.numFaces());
- using LocalIndexType = typename InteractionVolume::Traits::LocalIndexType;
+ using LocalIndexType = typename InteractionVolume::Traits::IndexSet::LocalIndexType;
for (LocalIndexType i = 0; i < iv.numFaces(); ++i)
{
const auto numNeighbors = iv.localScvf(i).neighboringLocalScvIndices().size() - 1;
@@ -229,12 +225,12 @@ public:
}
//! Access to the iv-wide mole fractions of a component in one phase
- const Vector& moleFractions(unsigned int pIdx, unsigned int compIdx) const { return xj_[pIdx][compIdx]; }
- Vector& moleFractions(unsigned int pIdx, unsigned int compIdx) { return xj_[pIdx][compIdx]; }
+ const CellVector& moleFractions(unsigned int pIdx, unsigned int compIdx) const { return xj_[pIdx][compIdx]; }
+ CellVector& moleFractions(unsigned int pIdx, unsigned int compIdx) { return xj_[pIdx][compIdx]; }
//! The transmissibilities associated with diffusive fluxes
- const Matrix& diffusionT() const { return T_[phaseIdx_][compIdx_]; }
- Matrix& diffusionT() { return T_[phaseIdx_][compIdx_]; }
+ const TMatrix& diffusionT() const { return T_[phaseIdx_][compIdx_]; }
+ TMatrix& diffusionT() { return T_[phaseIdx_][compIdx_]; }
//! The transmissibilities for "outside" faces (used on surface grids)
const OutsideTContainer& diffusionTout() const { return outsideT_[phaseIdx_][compIdx_]; }
@@ -242,22 +238,19 @@ public:
private:
//! diffusion-related variables
- unsigned int phaseIdx_; //!< The phase index set for the context
- unsigned int compIdx_; //!< The component index set for the context
- std::array< std::array<Matrix, numComponents>, numPhases > T_; //!< The transmissibilities such that f_i = T_ij*x_j
- std::array< std::array<Vector, numComponents>, numPhases > xj_; //!< The interaction volume-wide mole fractions
+ unsigned int phaseIdx_; //!< The phase index set for the context
+ unsigned int compIdx_; //!< The component index set for the context
+ std::array< std::array<TMatrix, numComponents>, numPhases > T_; //!< The transmissibilities such that f_i = T_ij*x_j
+ std::array< std::array<CellVector, numComponents>, numPhases > xj_; //!< The interaction volume-wide mole fractions
std::array< std::array<OutsideTContainer, numComponents>, numPhases> outsideT_;
};
//! Data handle for quantities related to heat conduction
-template<class TypeTag, class M, class V, class LI, bool EnableHeatConduction>
+template<class MatVecTraits, class PhysicsTraits, bool enableHeatConduction>
class HeatConductionDataHandle
{
- using Matrix = M;
- using Vector = V;
-
- static constexpr int dim = GET_PROP_TYPE(TypeTag, GridView)::dimension;
- static constexpr int dimWorld = GET_PROP_TYPE(TypeTag, GridView)::dimensionworld;
+ using TMatrix = typename MatVecTraits::TMatrix;
+ using CellVector = typename MatVecTraits::CellVector;
public:
//! prepare data handle for subsequent fill
@@ -269,11 +262,12 @@ public:
resizeVector(Tj_, iv.numKnowns());
//! resize outsideTij on surface grids
- if (dim < dimWorld)
+ using GridView = typename InteractionVolume::Traits::GridView;
+ if (int(GridView::dimension) < int(GridView::dimensionworld))
{
outsideT_.resize(iv.numFaces());
- using LocalIndexType = typename InteractionVolume::Traits::LocalIndexType;
+ using LocalIndexType = typename InteractionVolume::Traits::IndexSet::LocalIndexType;
for (LocalIndexType i = 0; i < iv.numFaces(); ++i)
{
const auto numNeighbors = iv.localScvf(i).neighboringLocalScvIndices().size() - 1;
@@ -285,49 +279,47 @@ public:
}
//! Access to the iv-wide temperatures
- const Vector& temperatures() const { return Tj_; }
- Vector& temperatures() { return Tj_; }
+ const CellVector& temperatures() const { return Tj_; }
+ CellVector& temperatures() { return Tj_; }
//! The transmissibilities associated with conductive fluxes
- const Matrix& heatConductionT() const { return T_; }
- Matrix& heatConductionT() { return T_; }
+ const TMatrix& heatConductionT() const { return T_; }
+ TMatrix& heatConductionT() { return T_; }
//! The transmissibilities for "outside" faces (used on surface grids)
- const std::vector< std::vector<Vector> >& heatConductionTout() const { return outsideT_; }
- std::vector< std::vector<Vector> >& heatConductionTout() { return outsideT_; }
+ const std::vector< std::vector<CellVector> >& heatConductionTout() const { return outsideT_; }
+ std::vector< std::vector<CellVector> >& heatConductionTout() { return outsideT_; }
private:
// heat conduction-related variables
- Matrix T_; //!< The transmissibilities such that f_i = T_ij*T_j
- Vector Tj_; //!< The interaction volume-wide temperatures
- std::vector< std::vector<Vector> > outsideT_; //!< The transmissibilities for "outside" faces (only necessary on surface grids)
+ TMatrix T_; //!< The transmissibilities such that f_i = T_ij*T_j
+ CellVector Tj_; //!< The interaction volume-wide temperatures
+ std::vector< std::vector<CellVector> > outsideT_; //!< The transmissibilities for "outside" faces (only necessary on surface grids)
};
//! Process-dependent data handles when related process is disabled
-template<class TypeTag, class M, class V, class LI>
-class AdvectionDataHandle<TypeTag, M, V, LI, false> : public EmptyDataHandle {};
-template<class TypeTag, class M, class V, class LI>
-class DiffusionDataHandle<TypeTag, M, V, LI, false> : public EmptyDataHandle {};
-template<class TypeTag, class M, class V, class LI>
-class HeatConductionDataHandle<TypeTag, M, V, LI, false> : public EmptyDataHandle {};
+template<class MatVecTraits, class PhysicsTraits>
+class AdvectionDataHandle<MatVecTraits, PhysicsTraits, false> : public EmptyDataHandle {};
+template<class MatVecTraits, class PhysicsTraits>
+class DiffusionDataHandle<MatVecTraits, PhysicsTraits, false> : public EmptyDataHandle {};
+template<class MatVecTraits, class PhysicsTraits>
+class HeatConductionDataHandle<MatVecTraits, PhysicsTraits, false> : public EmptyDataHandle {};
/*!
* \ingroup CCMpfaDiscretization
* \brief Class for the interaction volume data handle.
*
- * \tparam TypeTag The problem TypeTag
- * \tparam M The type used for iv-local matrices
- * \tparam V The type used for iv-local vectors
- * \tparam LI The type used for iv-local indexing
+ * \tparam MVT The matrix/vector traits collecting type information used by the iv
+ * \tparam PT The physics traits collecting data on the physical processes to be considered
*/
-template<class TypeTag, class M, class V, class LI>
-class InteractionVolumeDataHandle : public AdvectionDataHandle<TypeTag, M, V, LI, GET_PROP_VALUE(TypeTag, EnableAdvection)>,
- public DiffusionDataHandle<TypeTag, M, V, LI, GET_PROP_VALUE(TypeTag, EnableMolecularDiffusion)>,
- public HeatConductionDataHandle<TypeTag, M, V, LI, GET_PROP_VALUE(TypeTag, EnableEnergyBalance)>
+template<class MVT, class PT>
+class InteractionVolumeDataHandle : public AdvectionDataHandle<MVT, PT, PT::enableAdvection>,
+ public DiffusionDataHandle<MVT, PT, PT::enableMolecularDiffusion>,
+ public HeatConductionDataHandle<MVT, PT, PT::enableHeatConduction>
{
- using AdvectionHandle = AdvectionDataHandle<TypeTag, M, V, LI, GET_PROP_VALUE(TypeTag, EnableAdvection)>;
- using DiffusionHandle = DiffusionDataHandle<TypeTag, M, V, LI, GET_PROP_VALUE(TypeTag, EnableMolecularDiffusion)>;
- using HeatConductionHandle = HeatConductionDataHandle<TypeTag, M, V, LI, GET_PROP_VALUE(TypeTag, EnableEnergyBalance)>;
+ using AdvectionHandle = AdvectionDataHandle<MVT, PT, PT::enableAdvection>;
+ using DiffusionHandle = DiffusionDataHandle<MVT, PT, PT::enableMolecularDiffusion>;
+ using HeatConductionHandle = HeatConductionDataHandle<MVT, PT, PT::enableHeatConduction>;
public:
//! prepare data handles for subsequent fills
diff --git a/dumux/discretization/cellcentered/mpfa/localassembler.hh b/dumux/discretization/cellcentered/mpfa/localassembler.hh
index d5eb6bf8289f43f07a475a4cf63aaf24546e5b77..ad06867aee3fb0c3b0c2ca9b89eb32e9ad6e3e91 100644
--- a/dumux/discretization/cellcentered/mpfa/localassembler.hh
+++ b/dumux/discretization/cellcentered/mpfa/localassembler.hh
@@ -33,11 +33,11 @@
namespace Dumux
{
//! Forward declaration of the implementation
-template< class IVTraits, MpfaMethods M > class InteractionVolumeAssemblerImpl;
+template< class P, class EG, class EV, MpfaMethods M > class InteractionVolumeAssemblerImpl;
//! Alias to select the right implementation.
-template< class IVTraits, MpfaMethods M >
-using InteractionVolumeAssembler = InteractionVolumeAssemblerImpl< IVTraits, M >;
+template< class P, class EG, class EV, MpfaMethods M >
+using InteractionVolumeAssembler = InteractionVolumeAssemblerImpl< P, EG, EV, M >;
/*!
* \ingroup CCMpfaDiscretization
@@ -47,17 +47,16 @@ using InteractionVolumeAssembler = InteractionVolumeAssemblerImpl< IVTraits, M >
* for the available interaction volume implementations. these
* should derive from this base clases.
*
- * \tparam IVTraits The Traits class used by the interaction volume
+ * \tparam P The problem type
+ * \tparam EG The element finite volume geometry
+ * \tparam EV The element volume variables type
*/
-template< class IVTraits >
+template< class P, class EG, class EV >
class InteractionVolumeAssemblerBase
{
- using Matrix = typename IVTraits::Matrix;
- using Vector = typename IVTraits::Vector;
-
- using Problem = typename IVTraits::Problem;
- using FVElementGeometry = typename IVTraits::FVElementGeometry;
- using ElementVolumeVariables = typename IVTraits::ElementVolumeVariables;
+ using Problem = P;
+ using FVElementGeometry = EG;
+ using ElementVolumeVariables = EV;
public:
/*!
@@ -87,15 +86,15 @@ class InteractionVolumeAssemblerBase
* interaction volume-local transmissibility matrix.
*
* \tparam IV Interaction volume type implementation
- * \tparam GetTensor Lambda to obtain the tensor w.r.t.
- * which the local system is to be solved
+ * \tparam TensorFunc Lambda to obtain the tensor w.r.t.
+ * which the local system is to be solved
*
* \param T The transmissibility matrix to be assembled
* \param iv The interaction volume
- * \param getTensor Lambda to evaluate the scv-wise tensors
+ * \param getT Lambda to evaluate the scv-wise tensors
*/
- template< class IV, class GetTensor >
- void assemble(Matrix& T, IV& iv, const GetTensor& getTensor)
+ template< class IV, class TensorFunc >
+ void assemble(typename IV::Traits::MatVecTraits::MatVecTraits::TMatrix& T, IV& iv, const TensorFunc& getT)
{
DUNE_THROW(Dune::NotImplemented, "Implementation does not provide a assemble() function");
}
@@ -108,16 +107,16 @@ class InteractionVolumeAssemblerBase
*
* \tparam TOutside Container to store the "outside" transmissibilities
* \tparam IV Interaction volume type implementation
- * \tparam GetTensor Lambda to obtain the tensor w.r.t.
- * which the local system is to be solved
+ * \tparam TensorFunc Lambda to obtain the tensor w.r.t.
+ * which the local system is to be solved
*
* \param outsideTij tij on "outside" faces to be assembled
* \param T The transmissibility matrix tij to be assembled
* \param iv The mpfa-o interaction volume
- * \param getTensor Lambda to evaluate the scv-wise tensors
+ * \param getT Lambda to evaluate the scv-wise tensors
*/
- template< class TOutside, class IV, class GetTensor >
- void assemble(TOutside& outsideTij, Matrix& T, IV& iv, const GetTensor& getTensor)
+ template< class TOutside, class IV, class TensorFunc >
+ void assemble(TOutside& outsideTij, typename IV::Traits::MatVecTraits::TMatrix& T, IV& iv, const TensorFunc& getT)
{
DUNE_THROW(Dune::NotImplemented, "Implementation does not provide a assemble() function to be used on surface grids");
}
@@ -130,17 +129,21 @@ class InteractionVolumeAssemblerBase
* \tparam GC The type of container used to store the
* gravitational acceleration per scvf & phase
* \tparam IV The interaction volume type implementation
- * \tparam GetTensor Lambda to obtain the tensor w.r.t.
- * which the local system is to be solved
+ * \tparam TensorFunc Lambda to obtain the tensor w.r.t.
+ * which the local system is to be solved
*
* \param T The transmissibility matrix to be assembled
* \param g Container to assemble gravity per scvf & phase
* \param CA Matrix to store matrix product C*A^-1
* \param iv The interaction volume
- * \param getTensor Lambda to evaluate the scv-wise tensors
+ * \param getT Lambda to evaluate the scv-wise tensors
*/
- template< class GC, class IV, class GetTensor >
- void assembleWithGravity(Matrix& T, GC& g, Matrix& CA, IV& iv, const GetTensor& getTensor)
+ template< class GC, class IV, class TensorFunc >
+ void assembleWithGravity(typename IV::Traits::MatVecTraits::TMatrix& T,
+ GC& g,
+ typename IV::Traits::MatVecTraits::CMatrix& CA,
+ IV& iv,
+ const TensorFunc& getT)
{
DUNE_THROW(Dune::NotImplemented, "Implementation does not provide a assembleWithGravity() function");
}
@@ -158,8 +161,8 @@ class InteractionVolumeAssemblerBase
* \tparam GOut Type of container used to store gravity on "outside" faces
* \tparam TOutside Container to store the "outside" transmissibilities
* \tparam IV The interaction volume type implementation
- * \tparam GetTensor Lambda to obtain the tensor w.r.t.
- * which the local system is to be solved
+ * \tparam TensorFunc Lambda to obtain the tensor w.r.t.
+ * which the local system is to be solved
*
* \param outsideTij tij on "outside" faces to be assembled
* \param T The transmissibility matrix to be assembled
@@ -168,17 +171,17 @@ class InteractionVolumeAssemblerBase
* \param CA Matrix to store matrix product C*A^-1
* \param A Matrix to store the inverse A^-1
* \param iv The interaction volume
- * \param getTensor Lambda to evaluate the scv-wise tensors
+ * \param getT Lambda to evaluate the scv-wise tensors
*/
- template< class GC, class GOut, class TOutside, class IV, class GetTensor >
+ template< class GC, class GOut, class TOutside, class IV, class TensorFunc >
void assembleWithGravity(TOutside& outsideTij,
- Matrix& T,
+ typename IV::Traits::MatVecTraits::TMatrix& T,
GOut& outsideG,
GC& g,
- Matrix& CA,
- Matrix& A,
+ typename IV::Traits::MatVecTraits::CMatrix& CA,
+ typename IV::Traits::MatVecTraits::AMatrix& A,
IV& iv,
- const GetTensor& getTensor)
+ const TensorFunc& getT)
{
DUNE_THROW(Dune::NotImplemented, "Implementation does not provide a assembleWithGravity() function to be used on surface grids");
}
@@ -196,7 +199,7 @@ class InteractionVolumeAssemblerBase
* \param getU Lambda to obtain the desired cell value from grid indices
*/
template< class IV, class GetU >
- void assemble(Vector& u, const IV& iv, const GetU& getU)
+ void assemble(typename IV::Traits::MatVecTraits::CellVector& u, const IV& iv, const GetU& getU)
{
DUNE_THROW(Dune::NotImplemented, "Implementation does not provide a assemble() function for the cell unknowns");
}
@@ -211,17 +214,17 @@ class InteractionVolumeAssemblerBase
*
* \tparam GC The type of container used to store the
* gravitational acceleration per scvf & phase
- * \tparam GetTensor Lambda to obtain the tensor w.r.t.
- * which the local system is to be solved
+ * \tparam TensorFunc Lambda to obtain the tensor w.r.t.
+ * which the local system is to be solved
*
* \param g Container to assemble gravity per scvf & phase
* \param iv The interaction volume
* \param CA Projection matrix transforming the gravity terms in the local system of
* equations to the entire set of faces within the interaction volume
- * \param getTensor Lambda to evaluate scv-wise hydraulic conductivities
+ * \param getT Lambda to evaluate scv-wise hydraulic conductivities
*/
- template< class GC, class IV, class GetTensor >
- void assembleGravity(GC& g, const IV& iv, const Matrix& CA, const GetTensor& getTensor)
+ template< class GC, class IV, class TensorFunc >
+ void assembleGravity(GC& g, const IV& iv, const typename IV::Traits::MatVecTraits::CMatrix& CA, const TensorFunc& getT)
{
DUNE_THROW(Dune::NotImplemented, "Implementation does not provide a assembleGravity() function");
}
@@ -240,8 +243,8 @@ class InteractionVolumeAssemblerBase
* gravitational acceleration per scvf & phase
* \tparam GOut Type of container used to store gravity on "outside" faces
* \tparam IV The interaction volume type implementation
- * \tparam GetTensor Lambda to obtain the tensor w.r.t.
- * which the local system is to be solved
+ * \tparam TensorFunc Lambda to obtain the tensor w.r.t.
+ * which the local system is to be solved
*
* \param g Container to store gravity per scvf & phase
* \param outsideG Container to store gravity per "outside" scvf & phase
@@ -249,15 +252,15 @@ class InteractionVolumeAssemblerBase
* \param CA Projection matrix transforming the gravity terms in the local system of
* equations to the entire set of faces within the interaction volume
* \param A Matrix needed for the "reconstruction" of face unknowns as a function of gravity
- * \param getTensor Lambda to evaluate scv-wise hydraulic conductivities
+ * \param getT Lambda to evaluate scv-wise hydraulic conductivities
*/
- template< class GC, class GOut, class IV, class GetTensor >
+ template< class GC, class GOut, class IV, class TensorFunc >
void assembleGravity(GC& g,
GOut& outsideG,
const IV& iv,
- const Matrix& CA,
- const Matrix& A,
- const GetTensor& getTensor)
+ const typename IV::Traits::MatVecTraits::CMatrix& CA,
+ const typename IV::Traits::MatVecTraits::AMatrix& A,
+ const TensorFunc& getT)
{
DUNE_THROW(Dune::NotImplemented, "Implementation does not provide a assembleGravity() function to be used on surface grids");
}
diff --git a/dumux/discretization/cellcentered/mpfa/omethod/interactionvolume.hh b/dumux/discretization/cellcentered/mpfa/omethod/interactionvolume.hh
index be2e526ac16536fd6a71bc3d72877fc9dd503a4a..7549c2447956bc675a8198a0b0bab0d990dd705d 100644
--- a/dumux/discretization/cellcentered/mpfa/omethod/interactionvolume.hh
+++ b/dumux/discretization/cellcentered/mpfa/omethod/interactionvolume.hh
@@ -24,15 +24,16 @@
#ifndef DUMUX_DISCRETIZATION_CC_MPFA_O_INTERACTIONVOLUME_HH
#define DUMUX_DISCRETIZATION_CC_MPFA_O_INTERACTIONVOLUME_HH
+#include <type_traits>
+
#include <dune/common/dynmatrix.hh>
#include <dune/common/dynvector.hh>
#include <dune/common/fvector.hh>
+#include <dune/common/reservedvector.hh>
#include <dumux/common/math.hh>
-#include <dumux/common/properties.hh>
#include <dumux/common/matrixvectorhelper.hh>
-#include <dumux/discretization/cellcentered/mpfa/interactionvolumedatahandle.hh>
#include <dumux/discretization/cellcentered/mpfa/interactionvolumebase.hh>
#include <dumux/discretization/cellcentered/mpfa/dualgridindexset.hh>
#include <dumux/discretization/cellcentered/mpfa/localfacedata.hh>
@@ -46,45 +47,51 @@ namespace Dumux
//! Forward declaration of the o-method's interaction volume
template< class Traits > class CCMpfaOInteractionVolume;
-//! Specialization of the default interaction volume traits class for the mpfa-o method
-template< class TypeTag >
+/*!
+ * \ingroup CCMpfaDiscretization
+ * \brief The default interaction volume traits class for the mpfa-o method.
+ * This uses dynamic types types for matrices/vectors in order to work
+ * on general grids. For interaction volumes known at compile time use
+ * the static interaction volume implementation.
+ *
+ * \tparam NodalIndexSet The type used for the dual grid's nodal index sets
+ * \tparam Scalar The Type used for scalar values
+ */
+template< class NodalIndexSet, class Scalar >
struct CCMpfaODefaultInteractionVolumeTraits
{
private:
- 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;
+ using LocalIndexType = typename NodalIndexSet::LocalIndexType;
+
+ static constexpr int dim = NodalIndexSet::Traits::GridView::dimension;
+ static constexpr int dimWorld = NodalIndexSet::Traits::GridView::dimensionworld;
+
+ //! Matrix/Vector traits to be used by the data handle
+ struct MVTraits
+ {
+ using AMatrix = Dune::DynamicMatrix< Scalar >;
+ using BMatrix = Dune::DynamicMatrix< Scalar >;
+ using CMatrix = Dune::DynamicMatrix< Scalar >;
+ using DMatrix = Dune::DynamicMatrix< Scalar >;
+ using TMatrix = Dune::DynamicMatrix< Scalar >;
+ using CellVector = Dune::DynamicVector< Scalar >;
+ using FaceVector = Dune::DynamicVector< Scalar >;
+ };
public:
- //! export the problem type (needed for iv-local assembly)
- using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
- //! export the type of the local view on the finite volume grid geometry
- using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry)::LocalView;
- //! export the type of the local view on the grid volume variables
- using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);
//! export the type of grid view
- using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
- //! export the type used for scalar values
- using ScalarType = typename GET_PROP_TYPE(TypeTag, Scalar);
- //! 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 = typename NodalIndexSet::LocalIndexType;
+ using GridView = typename NodalIndexSet::Traits::GridView;
//! export the type for the interaction volume index set
using IndexSet = CCMpfaOInteractionVolumeIndexSet< NodalIndexSet >;
//! export the type of interaction-volume local scvs
- using LocalScvType = CCMpfaOInteractionVolumeLocalScv< IndexSet, ScalarType, dim, dimWorld >;
+ using LocalScvType = CCMpfaOInteractionVolumeLocalScv< IndexSet, Scalar, dim, dimWorld >;
//! export the type of interaction-volume local scvfs
using LocalScvfType = CCMpfaOInteractionVolumeLocalScvf< IndexSet >;
//! export the type of used for the iv-local face data
using LocalFaceData = InteractionVolumeLocalFaceData<GridIndexType, LocalIndexType>;
- //! export the type used for iv-local matrices
- using Matrix = Dune::DynamicMatrix< ScalarType >;
- //! export the type used for iv-local vectors
- using Vector = Dune::DynamicVector< ScalarType >;
- //! export the data handle type for this iv
- using DataHandle = InteractionVolumeDataHandle< TypeTag, Matrix, Vector, LocalIndexType >;
+ //! export the matrix/vector traits to be used by the iv
+ using MatVecTraits = MVTraits;
};
/*!
@@ -94,38 +101,41 @@ public:
* and can be used at boundaries and on unstructured grids.
*/
template< class Traits >
-class CCMpfaOInteractionVolume : public CCMpfaInteractionVolumeBase< CCMpfaOInteractionVolume<Traits>, Traits >
+class CCMpfaOInteractionVolume
+ : public CCMpfaInteractionVolumeBase< CCMpfaOInteractionVolume<Traits>, Traits >
{
- using ThisType = CCMpfaOInteractionVolume< Traits >;
- using ParentType = CCMpfaInteractionVolumeBase< CCMpfaOInteractionVolume<Traits>, Traits >;
-
- using Scalar = typename Traits::ScalarType;
- using Helper = typename Traits::MpfaHelper;
- using Problem = typename Traits::Problem;
- using FVElementGeometry = typename Traits::FVElementGeometry;
-
using GridView = typename Traits::GridView;
using Element = typename GridView::template Codim<0>::Entity;
+ using IndexSet = typename Traits::IndexSet;
+ using GridIndexType = typename IndexSet::GridIndexType;
+ using LocalIndexType = typename IndexSet::LocalIndexType;
+ using Stencil = typename IndexSet::NodalGridStencilType;
+
static constexpr int dim = GridView::dimension;
+ static constexpr int dimWorld = GridView::dimensionworld;
+
+ //! export scalar type from T matrix and define omegas
+ using Scalar = typename Traits::MatVecTraits::TMatrix::value_type;
using DimVector = Dune::FieldVector<Scalar, dim>;
+ using FaceOmegas = typename std::conditional< (dim<dimWorld),
+ std::vector<DimVector>,
+ Dune::ReservedVector<DimVector, 2> >::type;
+
+ using AMatrix = typename Traits::MatVecTraits::AMatrix;
+ using BMatrix = typename Traits::MatVecTraits::BMatrix;
+ using CMatrix = typename Traits::MatVecTraits::CMatrix;
- using Matrix = typename Traits::Matrix;
using LocalScvType = typename Traits::LocalScvType;
using LocalScvfType = typename Traits::LocalScvfType;
using LocalFaceData = typename Traits::LocalFaceData;
- using IndexSet = typename Traits::IndexSet;
- using GridIndexType = typename GridView::IndexSet::IndexType;
- using LocalIndexType = typename Traits::LocalIndexType;
- using Stencil = typename IndexSet::GridStencilType;
-
+public:
//! Data attached to scvf touching Dirichlet boundaries.
//! For the default o-scheme, we only store the corresponding vol vars index.
class DirichletData
{
GridIndexType volVarIndex_;
-
public:
//! Constructor
DirichletData(const GridIndexType index) : volVarIndex_(index) {}
@@ -134,11 +144,14 @@ class CCMpfaOInteractionVolume : public CCMpfaInteractionVolumeBase< CCMpfaOInte
GridIndexType volVarIndex() const { return volVarIndex_; }
};
-public:
+ //! export the type used for transmissibility storage
+ using TransmissibilityStorage = std::vector< FaceOmegas >;
+
//! publicly state the mpfa-scheme this interaction volume is associated with
static constexpr MpfaMethods MpfaMethod = MpfaMethods::oMethod;
//! Sets up the local scope for a given iv index set
+ template< class Problem, class FVElementGeometry >
void setUpLocalScope(const IndexSet& indexSet,
const Problem& problem,
const FVElementGeometry& fvGeometry)
@@ -147,7 +160,7 @@ public:
// index set of the dual grid's nodal index set
stencil_ = &indexSet.nodalIndexSet().globalScvIndices();
- // number of interaction-volume-local (= node-local for o-scheme) scvs/scvf
+ // number of interaction-volume-local scvs(=node-local for o-scheme) and scvfs
numFaces_ = indexSet.numFaces();
const auto numLocalScvs = indexSet.numScvs();
const auto numGlobalScvfs = indexSet.nodalIndexSet().numScvfs();
@@ -164,29 +177,26 @@ public:
dirichletData_.reserve(numFaces_);
localFaceData_.reserve(numGlobalScvfs);
- // set up quantities related to sub-control volumes
+ // set up stuff related to sub-control volumes
const auto& scvIndices = indexSet.globalScvIndices();
for (LocalIndexType scvIdxLocal = 0; scvIdxLocal < numLocalScvs; scvIdxLocal++)
{
- scvs_.emplace_back(Helper(),
+ elements_.emplace_back(fvGeometry.fvGridGeometry().element( scvIndices[scvIdxLocal] ));
+ scvs_.emplace_back(fvGeometry.fvGridGeometry().mpfaHelper(),
fvGeometry,
fvGeometry.scv( scvIndices[scvIdxLocal] ),
scvIdxLocal,
indexSet);
- elements_.emplace_back(fvGeometry.fvGridGeometry().element( scvIndices[scvIdxLocal] ));
}
// keep track of the number of unknowns etc
numUnknowns_ = 0;
numKnowns_ = numLocalScvs;
- // resize omega storage container
- wijk_.resize(numFaces_);
-
// set up quantitites related to sub-control volume faces
+ wijk_.resize(numFaces_);
for (LocalIndexType faceIdxLocal = 0; faceIdxLocal < numFaces_; ++faceIdxLocal)
{
- // get corresponding grid scvf
const auto& scvf = fvGeometry.scvf(indexSet.scvfIdxGlobal(faceIdxLocal));
// the neighboring scvs in local indices (order: 0 - inside scv, 1..n - outside scvs)
@@ -195,6 +205,10 @@ public:
// create local face data object for this face
localFaceData_.emplace_back(faceIdxLocal, neighborScvIndicesLocal[0], scvf.index());
+ // we will need as many omegas as scvs around the face
+ const auto numNeighborScvs = neighborScvIndicesLocal.size();
+ wijk_[faceIdxLocal].resize(numNeighborScvs);
+
// create iv-local scvf object
if (scvf.boundary())
{
@@ -207,18 +221,11 @@ public:
}
else
scvfs_.emplace_back(scvf, neighborScvIndicesLocal, numUnknowns_++, /*isDirichlet*/false);
-
- // on boundary faces we will only need one inside omega
- wijk_[faceIdxLocal].resize(1);
}
else
{
scvfs_.emplace_back(scvf, neighborScvIndicesLocal, numUnknowns_++, /*isDirichlet*/false);
- // we will need as many omegas as scvs around the face
- const auto numNeighborScvs = neighborScvIndicesLocal.size();
- wijk_[faceIdxLocal].resize(numNeighborScvs);
-
// add local face data objects for the outside faces
for (LocalIndexType i = 1; i < numNeighborScvs; ++i)
{
@@ -234,8 +241,12 @@ public:
outsideLocalScvIdx, // iv-local scv index
i-1, // scvf-local index in outside faces
flipScvf.index()); // global scvf index
+ break; // go to next outside face
}
}
+
+ // make sure we found it
+ assert(localFaceData_.back().ivLocalInsideScvIndex() == outsideLocalScvIdx);
}
}
}
@@ -262,13 +273,13 @@ public:
const Stencil& stencil() const { return *stencil_; }
//! returns the grid element corresponding to a given iv-local scv idx
- const Element& element(const LocalIndexType ivLocalScvIdx) const { return elements_[ivLocalScvIdx]; }
+ const Element& element(LocalIndexType ivLocalScvIdx) const { return elements_[ivLocalScvIdx]; }
//! returns the local scvf entity corresponding to a given iv-local scvf idx
- const LocalScvfType& localScvf(const LocalIndexType ivLocalScvfIdx) const { return scvfs_[ivLocalScvfIdx]; }
+ const LocalScvfType& localScvf(LocalIndexType ivLocalScvfIdx) const { return scvfs_[ivLocalScvfIdx]; }
//! returns the local scv entity corresponding to a given iv-local scv idx
- const LocalScvType& localScv(const LocalIndexType ivLocalScvIdx) const { return scvs_[ivLocalScvIdx]; }
+ const LocalScvType& localScv(LocalIndexType ivLocalScvIdx) const { return scvs_[ivLocalScvIdx]; }
//! returns a reference to the container with the local face data
const std::vector<LocalFaceData>& localFaceData() const { return localFaceData_; }
@@ -277,38 +288,41 @@ public:
const std::vector<DirichletData>& dirichletData() const { return dirichletData_; }
//! returns the matrix associated with face unknowns in local equation system
- const Matrix& A() const { return A_; }
- Matrix& A() { return A_; }
+ const AMatrix& A() const { return A_; }
+ AMatrix& A() { return A_; }
//! returns the matrix associated with cell unknowns in local equation system
- const Matrix& B() const { return B_; }
- Matrix& B() { return B_; }
+ const BMatrix& B() const { return B_; }
+ BMatrix& B() { return B_; }
//! returns the matrix associated with face unknowns in flux expressions
- const Matrix& C() const { return C_; }
- Matrix& C() { return C_; }
+ const CMatrix& C() const { return C_; }
+ CMatrix& C() { return C_; }
//! returns container storing the transmissibilities for each face & coordinate
- const std::vector< std::vector< DimVector > >& omegas() const { return wijk_; }
- std::vector< std::vector< DimVector > >& omegas() { return wijk_; }
+ const TransmissibilityStorage& omegas() const { return wijk_; }
+ TransmissibilityStorage& omegas() { return wijk_; }
//! returns the number of interaction volumes living around a vertex
- //! the mpfa-o scheme always constructs one iv per vertex
- template< class NodalIndexSet >
- static std::size_t numInteractionVolumesAtVertex(const NodalIndexSet& nodalIndexSet) { return 1; }
+ template< class NI >
+ static constexpr std::size_t numIVAtVertex(const NI& nodalIndexSet) { return 1; }
//! adds the iv index sets living around a vertex to a given container
//! and stores the the corresponding index in a map for each scvf
- template<class IvIndexSetContainer, class ScvfIndexMap, class NodalIndexSet>
- static void addInteractionVolumeIndexSets(IvIndexSetContainer& ivIndexSetContainer,
- ScvfIndexMap& scvfIndexMap,
- const NodalIndexSet& nodalIndexSet)
+ template< class IvIndexSetContainer,
+ class ScvfIndexMap,
+ class NodalIndexSet,
+ class FlipScvfIndexSet >
+ static void addIVIndexSets(IvIndexSetContainer& ivIndexSetContainer,
+ ScvfIndexMap& scvfIndexMap,
+ const NodalIndexSet& nodalIndexSet,
+ const FlipScvfIndexSet& flipScvfIndexSet)
{
// the global index of the iv index set that is about to be created
const auto curGlobalIndex = ivIndexSetContainer.size();
// make the one index set for this node
- ivIndexSetContainer.emplace_back(nodalIndexSet);
+ ivIndexSetContainer.emplace_back(nodalIndexSet, flipScvfIndexSet);
// store the index mapping
for (const auto scvfIdx : nodalIndexSet.globalScvfIndices())
@@ -327,12 +341,12 @@ private:
std::vector<DirichletData> dirichletData_;
// Matrices needed for computation of transmissibilities
- Matrix A_;
- Matrix B_;
- Matrix C_;
+ AMatrix A_;
+ BMatrix B_;
+ CMatrix C_;
// The omega factors are stored during assembly of local system
- std::vector< std::vector< DimVector > > wijk_;
+ TransmissibilityStorage wijk_;
// sizes involved in the local system equations
std::size_t numFaces_;
diff --git a/dumux/discretization/cellcentered/mpfa/omethod/interactionvolumeindexset.hh b/dumux/discretization/cellcentered/mpfa/omethod/interactionvolumeindexset.hh
index 1d63d804f259675a3619111b2056aa24dcd9430f..51894f6ce22f42b0c708117f76a06b45fe062243 100644
--- a/dumux/discretization/cellcentered/mpfa/omethod/interactionvolumeindexset.hh
+++ b/dumux/discretization/cellcentered/mpfa/omethod/interactionvolumeindexset.hh
@@ -47,22 +47,23 @@ public:
using LocalIndexType = typename DualGridNodalIndexSet::LocalIndexType;
using GridIndexType = typename DualGridNodalIndexSet::GridIndexType;
- // 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 stencil types used
+ using NodalGridStencilType = typename DualGridNodalIndexSet::NodalGridStencilType;
+ using NodalLocalStencilType = typename DualGridNodalIndexSet::NodalLocalStencilType;
+ using NodalGridScvfStencilType = typename DualGridNodalIndexSet::NodalGridScvfStencilType;
//! 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)
+ template< class FlipScvfIndexSet >
+ CCMpfaOInteractionVolumeIndexSet(const NodalIndexSet& nodalIndexSet, const FlipScvfIndexSet& flipScvfIndexSet)
+ : 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();
- // keeps track of which nodal scvfs have been handled already
+ // kee track of which nodal scvfs have been handled already
+ nodeToIvScvf_.resize(numNodalScvfs);
std::vector<bool> isHandled(numNodalScvfs, false);
// go over faces in nodal index set, check if iv-local face has been
@@ -77,6 +78,7 @@ public:
// for scvfs touching the boundary there are no "outside" scvfs
if (nodalIndexSet.scvfIsOnBoundary(i))
{
+ scvfNeighborScvLocalIndices_.push_back({nodalIndexSet.insideScvIdxLocal(i)});
nodeToIvScvf_[i] = ivToNodeScvf_.size();
isHandled[i] = true;
ivToNodeScvf_.push_back(i);
@@ -84,79 +86,46 @@ public:
continue;
}
- // We insert a new iv-local face and find all "outside" scvfs that map
- // to this face as well by comparing the set of neighboring scv indices.
- const auto scvIndices = [&nodalIndexSet, i] ()
- {
- auto tmp = nodalIndexSet.neighboringScvIndices(i);
- std::sort(tmp.begin(), tmp.end());
- return tmp;
- } ();
- const auto numNeighborsI = scvIndices.size();
-
- std::vector<LocalIndexType> outsideScvfs;
- for (LocalIndexType j = i+1; j < numNodalScvfs; ++j)
- {
- // a face that has been handled already cannot be an "outside" face here
- if (!isHandled[j])
- {
- const auto scvIndices2 = [&nodalIndexSet, j] ()
- {
- auto tmp = nodalIndexSet.neighboringScvIndices(j);
- std::sort(tmp.begin(), tmp.end());
- return tmp;
- } ();
-
- // if the sizes aren't equal, this can't be an "outside" face
- if (scvIndices2.size() != numNeighborsI)
- continue;
-
- // if the neighboring scv indices are the same, this is an "outside" face
- if (std::equal(scvIndices.begin(), scvIndices.end(), scvIndices2.begin()))
- outsideScvfs.push_back(j);
- }
- }
-
- // insert mappings
+ // insert a new iv-local face
const auto curIvLocalScvfIdx = ivToNodeScvf_.size();
nodeToIvScvf_[i] = curIvLocalScvfIdx;
isHandled[i] = true;
- for (const auto nodeLocalScvfIdx : outsideScvfs)
+
+ // construct local index sets
+ const auto& flipScvfIndices = flipScvfIndexSet[nodalIndexSet.scvfIdxGlobal(i)];
+ const auto numFlipIndices = flipScvfIndices.size();
+
+ ScvfNeighborLocalIndexSet neighborsLocal;
+ neighborsLocal.resize(numFlipIndices + 1);
+ neighborsLocal[0] = nodalIndexSet.insideScvIdxLocal(i);
+
+ // mappings for all flip scvf
+ for (unsigned int j = 0; j < numFlipIndices; ++j)
{
- nodeToIvScvf_[nodeLocalScvfIdx] = curIvLocalScvfIdx;
- isHandled[nodeLocalScvfIdx] = true;
+ const auto outsideScvfIdx = flipScvfIndices[j];
+ for (unsigned int nodeLocalIdx = 0; nodeLocalIdx < nodalIndexSet.numScvfs(); ++nodeLocalIdx)
+ {
+ if (nodalIndexSet.scvfIdxGlobal(nodeLocalIdx) == outsideScvfIdx)
+ {
+ neighborsLocal[j+1] = nodalIndexSet.insideScvIdxLocal(nodeLocalIdx);
+ nodeToIvScvf_[nodeLocalIdx] = curIvLocalScvfIdx;
+ isHandled[nodeLocalIdx] = true;
+ break; // go to next outside scvf
+ }
+ }
}
+
+ scvfNeighborScvLocalIndices_.push_back(neighborsLocal);
ivToNodeScvf_.push_back(i);
numFaces_++;
}
-
- // compute local neighboring scv indices for the iv-local scvfs
- scvfNeighborScvLocalIndices_.resize(numFaces_);
- for (unsigned int i = 0; i < numFaces_; ++i)
- {
- const auto& neighborsGlobal = nodalIndexSet_.neighboringScvIndices(ivToNodeScvf_[i]);
- const auto numNeighbors = nodalIndexSet_.scvfIsOnBoundary(ivToNodeScvf_[i]) ? 1 : neighborsGlobal.size();
-
- scvfNeighborScvLocalIndices_[i].resize(numNeighbors);
- for (unsigned int j = 0; j < numNeighbors; ++j)
- scvfNeighborScvLocalIndices_[i][j] = findLocalScvIdx_(neighborsGlobal[j]);
- }
}
//! returns the corresponding nodal index set
- const DualGridNodalIndexSet& nodalIndexSet() const { return nodalIndexSet_; }
-
- //! returns a global scvf idx for a given iv_local scvf index
- GridIndexType scvfIdxGlobal(LocalIndexType ivLocalScvfIdx) const
- { return nodalIndexSet_.scvfIdxGlobal( ivToNodeScvf_[ivLocalScvfIdx] ); }
+ const NodalIndexSet& nodalIndexSet() const { return nodalIndexSet_; }
- //! returns the iv-local scvf idx of the i-th scvfs embedded in a local scv
- LocalIndexType scvfIdxLocal(LocalIndexType scvIdxLocal, unsigned int i) const
- { return nodeToIvScvf_[ nodalIndexSet_.scvfIdxLocal(scvIdxLocal, i) ]; }
-
- //! returns the local indices of the neighboring scvs of an scvf
- const ScvfNeighborLocalIndexSet& neighboringLocalScvIndices(LocalIndexType ivLocalScvfIdx) const
- { return scvfNeighborScvLocalIndices_[ivLocalScvfIdx]; }
+ //! returns the global scv indices connected to this dual grid node
+ const NodalGridStencilType& globalScvIndices() const { return nodalIndexSet_.globalScvIndices(); }
//! returns the number of faces in the interaction volume
std::size_t numFaces() const { return numFaces_; }
@@ -164,30 +133,40 @@ public:
//! returns the number of scvs in the interaction volume
std::size_t numScvs() const { return nodalIndexSet_.numScvs(); }
- //! returns the global scv indices connected to this dual grid node
- const GridStencilType& globalScvIndices() const { return nodalIndexSet_.globalScvIndices(); }
+ //! returns a global scvf idx for a given iv-local scvf index
+ GridIndexType scvfIdxGlobal(LocalIndexType ivLocalScvfIdx) const
+ {
+ assert(ivLocalScvfIdx < numFaces());
+ return nodalIndexSet_.scvfIdxGlobal( ivToNodeScvf_[ivLocalScvfIdx] );
+ }
- //! returns the global scvf indices connected to this dual grid node
- const GridScvfStencilType& globalScvfIndices() const { return nodalIndexSet_.globalScvfIndices(); }
+ //! returns the iv-local scvf idx of the i-th scvf embedded in a local scv
+ LocalIndexType scvfIdxLocal(LocalIndexType scvIdxLocal, unsigned int i) const
+ {
+ assert(nodalIndexSet_.scvfIdxLocal(scvIdxLocal, i) < nodeToIvScvf_.size());
+ return nodeToIvScvf_[ nodalIndexSet_.scvfIdxLocal(scvIdxLocal, i) ];
+ }
-private:
- //! returns the local scv index to a given global scv index
- unsigned int findLocalScvIdx_(GridIndexType globalScvIdx) const
+ //! returns the local indices of the neighboring scvs of an scvf
+ const ScvfNeighborLocalIndexSet& neighboringLocalScvIndices(LocalIndexType ivLocalScvfIdx) const
{
- auto it = std::find( nodalIndexSet_.globalScvIndices().begin(), nodalIndexSet_.globalScvIndices().end(), globalScvIdx );
- assert(it != nodalIndexSet_.globalScvIndices().end() && "Global scv index not found in local container!");
- return std::distance(nodalIndexSet_.globalScvIndices().begin(), it);
+ assert(ivLocalScvfIdx < numFaces());
+ return scvfNeighborScvLocalIndices_[ivLocalScvfIdx];
}
- const DualGridNodalIndexSet& nodalIndexSet_;
+private:
+ using NI = NodalIndexSet;
std::size_t numFaces_;
- Dune::ReservedVector< LocalIndexType, NodalIndexSet::maxNumScvfsAtNode > ivToNodeScvf_;
- Dune::ReservedVector< LocalIndexType, NodalIndexSet::maxNumScvfsAtNode > nodeToIvScvf_;
-
+ const NI& nodalIndexSet_;
+ // Index maps from and to nodal index set. For the map to the
+ // nodal set we use the same storage type as we know the nodal
+ // has more faces, thus sufficient guaranteed here!
+ typename NI::Traits::template NodalScvfDataStorage< LocalIndexType > ivToNodeScvf_;
+ typename NI::Traits::template NodalScvfDataStorage< LocalIndexType > nodeToIvScvf_;
// maps to each scvf a list of neighbouring scv indices
// ordering: 0 - inside scv idx; 1..n - outside scv indices
- Dune::ReservedVector< ScvfNeighborLocalIndexSet, NodalIndexSet::maxNumScvfsAtNode > scvfNeighborScvLocalIndices_;
+ typename NI::Traits::template NodalScvfDataStorage< ScvfNeighborLocalIndexSet > scvfNeighborScvLocalIndices_;
};
} // end namespace Dumux
diff --git a/dumux/discretization/cellcentered/mpfa/omethod/localassembler.hh b/dumux/discretization/cellcentered/mpfa/omethod/localassembler.hh
index fd6700c7273d300f4433326e262320f653e77614..50116c73c6426a33c0309c797891c1dc7c801144 100644
--- a/dumux/discretization/cellcentered/mpfa/omethod/localassembler.hh
+++ b/dumux/discretization/cellcentered/mpfa/omethod/localassembler.hh
@@ -41,20 +41,15 @@ namespace Dumux
* \brief Specialization of the interaction volume-local
* assembler class for the schemes using an mpfa-o type assembly.
*
- * \tparam IVTraits The traits class used by the interaction volume implemetation
+ * \tparam P The problem type
+ * \tparam EG The element finite volume geometry
+ * \tparam EV The element volume variables type
*/
-template< class IVTraits >
-class InteractionVolumeAssemblerImpl< IVTraits, MpfaMethods::oMethod >
- : public InteractionVolumeAssemblerBase< IVTraits >
+template< class P, class EG, class EV >
+class InteractionVolumeAssemblerImpl< P, EG, EV, MpfaMethods::oMethod >
+ : public InteractionVolumeAssemblerBase< P, EG, EV >
{
- using ParentType = InteractionVolumeAssemblerBase< IVTraits >;
-
- using LocalIndexType = typename IVTraits::LocalIndexType;
- using Matrix = typename IVTraits::Matrix;
- using Vector = typename IVTraits::Vector;
-
- static constexpr int dim = IVTraits::LocalScvType::myDimension;
- static constexpr int dimWorld = IVTraits::LocalScvType::worldDimension;
+ using ParentType = InteractionVolumeAssemblerBase< P, EG, EV >;
public:
//! Use the constructor of the base class
@@ -65,18 +60,18 @@ public:
* interaction volume in an mpfa-o type way.
*
* \tparam IV The interaction volume type implementation
- * \tparam GetTensor Lambda to obtain the tensor w.r.t.
+ * \tparam TensorFunc Lambda to obtain the tensor w.r.t.
* which the local system is to be solved
*
* \param T The transmissibility matrix to be assembled
* \param iv The interaction volume
- * \param getTensor Lambda to evaluate the scv-wise tensors
+ * \param getT Lambda to evaluate the scv-wise tensors
*/
- template< class IV, class GetTensor >
- void assemble(Matrix& T, IV& iv, const GetTensor& getTensor)
+ template< class IV, class TensorFunc >
+ void assemble(typename IV::Traits::MatVecTraits::TMatrix& T, IV& iv, const TensorFunc& getT)
{
// assemble D into T directly
- assembleLocalMatrices_(iv.A(), iv.B(), iv.C(), T, iv, getTensor);
+ assembleLocalMatrices_(iv.A(), iv.B(), iv.C(), T, iv, getT);
// maybe solve the local system
if (iv.numUnknowns() > 0)
@@ -95,19 +90,19 @@ public:
*
* \tparam TOutside Container to store the "outside" transmissibilities
* \tparam IV The interaction volume type implementation
- * \tparam GetTensor Lambda to obtain the tensor w.r.t.
- * which the local system is to be solved
+ * \tparam TensorFunc Lambda to obtain the tensor w.r.t.
+ * which the local system is to be solved
*
* \param outsideTij tij on "outside" faces to be assembled
* \param T The transmissibility matrix tij to be assembled
* \param iv The interaction volume
- * \param getTensor Lambda to evaluate the scv-wise tensors
+ * \param getT Lambda to evaluate the scv-wise tensors
*/
- template< class TOutside, class IV, class GetTensor >
- void assemble(TOutside& outsideTij, Matrix& T, IV& iv, const GetTensor& getTensor)
+ template< class TOutside, class IV, class TensorFunc >
+ void assemble(TOutside& outsideTij, typename IV::Traits::MatVecTraits::TMatrix& T, IV& iv, const TensorFunc& getT)
{
// assemble D into T directly
- assembleLocalMatrices_(iv.A(), iv.B(), iv.C(), T, iv, getTensor);
+ assembleLocalMatrices_(iv.A(), iv.B(), iv.C(), T, iv, getT);
// maybe solve the local system
if (iv.numUnknowns() > 0)
@@ -135,7 +130,8 @@ public:
tij = 0.0;
// add contributions from all local directions
- for (LocalIndexType localDir = 0; localDir < dim; localDir++)
+ using LocalIndexType = typename IV::Traits::IndexSet::LocalIndexType;
+ for (LocalIndexType localDir = 0; localDir < IV::Traits::GridView::dimension; localDir++)
{
// the scvf corresponding to this local direction in the scv
const auto& curLocalScvf = iv.localScvf(posLocalScv.scvfIdxLocal(localDir));
@@ -166,20 +162,24 @@ public:
* \tparam GC The type of container used to store the
* gravitational acceleration per scvf & phase
* \tparam IV The interaction volume type implementation
- * \tparam GetTensor Lambda to obtain the tensor w.r.t.
- * which the local system is to be solved
+ * \tparam TensorFunc Lambda to obtain the tensor w.r.t.
+ * which the local system is to be solved
*
* \param T The transmissibility matrix to be assembled
* \param g Container to assemble gravity per scvf & phase
* \param CA Matrix to store matrix product C*A^-1
* \param iv The mpfa-o interaction volume
- * \param getTensor Lambda to evaluate the scv-wise tensors
+ * \param getT Lambda to evaluate the scv-wise tensors
*/
- template< class GC, class IV, class GetTensor >
- void assembleWithGravity(Matrix& T, GC& g, Matrix& CA, IV& iv, const GetTensor& getTensor)
+ template< class GC, class IV, class TensorFunc >
+ void assembleWithGravity(typename IV::Traits::MatVecTraits::TMatrix& T,
+ GC& g,
+ typename IV::Traits::MatVecTraits::CMatrix& CA,
+ IV& iv,
+ const TensorFunc& getT)
{
// assemble D into T & C into CA directly
- assembleLocalMatrices_(iv.A(), iv.B(), CA, T, iv, getTensor);
+ assembleLocalMatrices_(iv.A(), iv.B(), CA, T, iv, getT);
// maybe solve the local system
if (iv.numUnknowns() > 0)
@@ -191,7 +191,7 @@ public:
}
// assemble gravitational acceleration container (enforce usage of mpfa-o type version)
- assembleGravity(g, iv, CA, getTensor);
+ assembleGravity(g, iv, CA, getT);
}
/*!
@@ -206,8 +206,8 @@ public:
* \tparam GOut Type of container used to store gravity on "outside" faces
* \tparam TOutside Container to store the "outside" transmissibilities
* \tparam IV The interaction volume type implementation
- * \tparam GetTensor Lambda to obtain the tensor w.r.t.
- * which the local system is to be solved
+ * \tparam TensorFunc Lambda to obtain the tensor w.r.t.
+ * which the local system is to be solved
*
* \param outsideTij tij on "outside" faces to be assembled
* \param T The transmissibility matrix to be assembled
@@ -216,20 +216,20 @@ public:
* \param CA Matrix to store matrix product C*A^-1
* \param A Matrix to store the inverse A^-1
* \param iv The mpfa-o interaction volume
- * \param getTensor Lambda to evaluate the scv-wise tensors
+ * \param getT Lambda to evaluate the scv-wise tensors
*/
- template< class GC, class GOut, class TOutside, class IV, class GetTensor >
+ template< class GC, class GOut, class TOutside, class IV, class TensorFunc >
void assembleWithGravity(TOutside& outsideTij,
- Matrix& T,
+ typename IV::Traits::MatVecTraits::TMatrix& T,
GOut& outsideG,
GC& g,
- Matrix& CA,
- Matrix& A,
+ typename IV::Traits::MatVecTraits::CMatrix& CA,
+ typename IV::Traits::MatVecTraits::AMatrix& A,
IV& iv,
- const GetTensor& getTensor)
+ const TensorFunc& getT)
{
// assemble D into T directly
- assembleLocalMatrices_(iv.A(), iv.B(), iv.C(), T, iv, getTensor);
+ assembleLocalMatrices_(iv.A(), iv.B(), iv.C(), T, iv, getT);
// maybe solve the local system
if (iv.numUnknowns() > 0)
@@ -259,7 +259,8 @@ public:
tij = 0.0;
// add contributions from all local directions
- for (LocalIndexType localDir = 0; localDir < dim; localDir++)
+ using LocalIndexType = typename IV::Traits::IndexSet::LocalIndexType;
+ for (LocalIndexType localDir = 0; localDir < IV::Traits::GridView::dimension; localDir++)
{
// the scvf corresponding to this local direction in the scv
const auto& curLocalScvf = iv.localScvf(posLocalScv.scvfIdxLocal(localDir));
@@ -280,7 +281,7 @@ public:
}
}
- assembleGravity(g, outsideG, iv, CA, A, getTensor);
+ assembleGravity(g, outsideG, iv, CA, A, getT);
}
/*!
@@ -295,14 +296,15 @@ public:
* \param getU Lambda to obtain the desired cell/Dirichlet value from grid index
*/
template< class IV, class GetU >
- void assemble(Vector& u, const IV& iv, const GetU& getU)
+ void assemble(typename IV::Traits::MatVecTraits::CellVector& u, const IV& iv, const GetU& getU)
{
// put the cell pressures first
+ using LocalIndexType = typename IV::Traits::IndexSet::LocalIndexType;
for (LocalIndexType i = 0; i < iv.numScvs(); ++i)
u[i] = getU( iv.localScv(i).globalScvIndex() );
// Dirichlet BCs come afterwards
- unsigned int i = iv.numScvs();
+ LocalIndexType i = iv.numScvs();
for (const auto& data : iv.dirichletData())
u[i++] = getU( data.volVarIndex() );
}
@@ -317,37 +319,40 @@ public:
*
* \tparam GC The type of container used to store the
* gravitational acceleration per scvf & phase
- * \tparam GetTensor Lambda to obtain the tensor w.r.t.
- * which the local system is to be solved
+ * \tparam TensorFunc Lambda to obtain the tensor w.r.t.
+ * which the local system is to be solved
*
* \param g Container to assemble gravity per scvf & phase
* \param iv The mpfa-o interaction volume
* \param CA Projection matrix transforming the gravity terms in the local system of
* equations to the entire set of faces within the interaction volume
- * \param getTensor Lambda to evaluate scv-wise hydraulic conductivities
+ * \param getT Lambda to evaluate scv-wise hydraulic conductivities
*/
- template< class GC, class IV, class GetTensor >
- void assembleGravity(GC& g, const IV& iv, const Matrix& CA, const GetTensor& getTensor)
+ template< class GC, class IV, class TensorFunc >
+ void assembleGravity(GC& g,
+ const IV& iv,
+ const typename IV::Traits::MatVecTraits::CMatrix& CA,
+ const TensorFunc& getT)
{
//! We require the gravity container to be a two-dimensional vector/array type, structured as follows:
//! - first index adresses the respective phases
//! - second index adresses the face within the interaction volume
- // make sure CA matrix and the g vector to have the correct size already
- assert(std::all_of(g.begin(), g.end(), [&iv](const auto& v) { return v.size() == iv.numFaces(); })
- && "size of gravity vector does not match the number of iv-local faces!");
- assert(CA.rows() == iv.numFaces() && CA.cols() == iv.numUnknowns() && "Matrix CA does not have the correct size");
+ // make sure g vector and CA matrix have the correct sizes already
+ assert(std::all_of(g.begin(), g.end(), [&iv](const auto& v) { return v.size() == iv.numFaces(); }));
+ assert(CA.rows() == iv.numFaces() && CA.cols() == iv.numUnknowns());
//! For each face, we...
//! - arithmetically average the phase densities
//! - compute the term \f$ \alpha := A \rho \ \mathbf{n}^T \mathbf{K} \mathbf{g} \f$ in each neighboring cell
//! - compute \f$ \alpha^* = \alpha_{outside} - \alpha_{inside} \f$
- using Scalar = typename Vector::value_type;
+ using Scalar = typename IV::Traits::MatVecTraits::TMatrix::value_type;
+ using FaceVector = typename IV::Traits::MatVecTraits::FaceVector;
+ using LocalIndexType = typename IV::Traits::IndexSet::LocalIndexType;
// reset gravity containers to zero
const auto numPhases = g.size();
- std::vector< Vector > sum_alphas(numPhases);
-
+ std::vector< FaceVector > sum_alphas(numPhases);
for (unsigned int pIdx = 0; pIdx < numPhases; ++pIdx)
{
resizeVector(sum_alphas[pIdx], iv.numUnknowns());
@@ -368,11 +373,10 @@ public:
const auto& posGlobalScv = this->fvGeometry().scv(posLocalScv.globalScvIndex());
const auto& posVolVars = this->elemVolVars()[posGlobalScv];
const auto& posElement = iv.element(neighborScvIndices[0]);
- const auto tensor = getTensor(this->problem(), posElement, posVolVars, this->fvGeometry(), posGlobalScv);
+ const auto tensor = getT(this->problem(), posElement, posVolVars, this->fvGeometry(), posGlobalScv);
- // This function should never be called for surface grids,
- // for which there is the specialization of this function below
- assert(neighborScvIndices.size() <= 2 && "Scvf seems to have more than one outside scv!");
+ // On surface grids one should use the function specialization below
+ assert(neighborScvIndices.size() <= 2);
std::vector< Scalar > rho(numPhases);
const auto alpha_inside = posVolVars.extrusionFactor()*vtmv(curGlobalScvf.unitOuterNormal(), tensor, gravity);
@@ -388,7 +392,7 @@ public:
const auto& negGlobalScv = this->fvGeometry().scv(negLocalScv.globalScvIndex());
const auto& negVolVars = this->elemVolVars()[negGlobalScv];
const auto& negElement = iv.element( neighborScvIndices[1] );
- const auto negTensor = getTensor(this->problem(), negElement, negVolVars, this->fvGeometry(), negGlobalScv);
+ const auto negTensor = getT(this->problem(), negElement, negVolVars, this->fvGeometry(), negGlobalScv);
const auto sum_alpha = negVolVars.extrusionFactor()
* vtmv(curGlobalScvf.unitOuterNormal(), negTensor, gravity)
@@ -439,8 +443,8 @@ public:
* gravitational acceleration per scvf & phase
* \tparam GOut Type of container used to store gravity on "outside" faces
* \tparam IV The interaction volume type implementation
- * \tparam GetTensor Lambda to obtain the tensor w.r.t.
- * which the local system is to be solved
+ * \tparam TensorFunc Lambda to obtain the tensor w.r.t.
+ * which the local system is to be solved
*
* \param g Container to store gravity per scvf & phase
* \param outsideG Container to store gravity per "outside" scvf & phase
@@ -448,15 +452,15 @@ public:
* \param CA Projection matrix transforming the gravity terms in the local system of
* equations to the entire set of faces within the interaction volume
* \param A Matrix needed for the "reconstruction" of face unknowns as a function of gravity
- * \param getTensor Lambda to evaluate scv-wise hydraulic conductivities
+ * \param getT Lambda to evaluate scv-wise hydraulic conductivities
*/
- template< class GC, class GOut, class IV, class GetTensor >
+ template< class GC, class GOut, class IV, class TensorFunc >
void assembleGravity(GC& g,
GOut& outsideG,
const IV& iv,
- const Matrix& CA,
- const Matrix& A,
- const GetTensor& getTensor)
+ const typename IV::Traits::MatVecTraits::CMatrix& CA,
+ const typename IV::Traits::MatVecTraits::AMatrix& A,
+ const TensorFunc& getT)
{
//! We require the gravity container to be a two-dimensional vector/array type, structured as follows:
//! - first index adresses the respective phases
@@ -467,22 +471,21 @@ public:
//! - third index adresses the i-th "outside" face of the current face
// we require the CA matrix and the gravity containers to have the correct size already
- assert(CA.rows() == iv.numFaces() && CA.cols() == iv.numUnknowns() && "Matrix CA does not have the correct size");
- assert(std::all_of(g.begin(), g.end(), [&iv](const auto& v) { return v.size() == iv.numFaces(); })
- && "size of gravity container does not match the number of iv-local faces!");
- assert(std::all_of(outsideG.begin(), outsideG.end(), [&iv](const auto& v) { return v.size() == iv.numFaces(); })
- && "Outside gravity container does not match the number of iv-local faces!");
+ assert(CA.rows() == iv.numFaces() && CA.cols() == iv.numUnknowns());
+ assert(std::all_of(g.begin(), g.end(), [&iv](const auto& v) { return v.size() == iv.numFaces(); }));
+ assert(std::all_of(outsideG.begin(), outsideG.end(), [&iv](const auto& v) { return v.size() == iv.numFaces(); }));
//! For each face, we...
//! - arithmetically average the phase densities
//! - compute the term \f$ \alpha := \mathbf{A} \rho \ \mathbf{n}^T \mathbf{K} \mathbf{g} \f$ in each neighboring cell
//! - compute \f$ \alpha^* = \sum{\alpha_{outside, i}} - \alpha_{inside} \f$
- using Scalar = typename Vector::value_type;
+ using Scalar = typename IV::Traits::MatVecTraits::TMatrix::value_type;
+ using FaceVector = typename IV::Traits::MatVecTraits::FaceVector;
+ using LocalIndexType = typename IV::Traits::IndexSet::LocalIndexType;
// reset everything to zero
const auto numPhases = g.size();
- std::vector< Vector > sum_alphas(numPhases);
-
+ std::vector< FaceVector > sum_alphas(numPhases);
for (unsigned int pIdx = 0; pIdx < numPhases; ++pIdx)
{
resizeVector(sum_alphas[pIdx], iv.numUnknowns());
@@ -504,7 +507,7 @@ public:
const auto& posGlobalScv = this->fvGeometry().scv(posLocalScv.globalScvIndex());
const auto& posVolVars = this->elemVolVars()[posGlobalScv];
const auto& posElement = iv.element(neighborScvIndices[0]);
- const auto tensor = getTensor(this->problem(), posElement, posVolVars, this->fvGeometry(), posGlobalScv);
+ const auto tensor = getT(this->problem(), posElement, posVolVars, this->fvGeometry(), posGlobalScv);
const auto alpha_inside = posVolVars.extrusionFactor()*vtmv(curGlobalScvf.unitOuterNormal(), tensor, gravity);
const auto numOutsideFaces = curGlobalScvf.boundary() ? 0 : curGlobalScvf.numOutsideScvs();
@@ -527,11 +530,10 @@ public:
const auto& negGlobalScv = this->fvGeometry().scv(negLocalScv.globalScvIndex());
const auto& negVolVars = this->elemVolVars()[negGlobalScv];
const auto& negElement = iv.element( neighborScvIndices[idxInOutside] );
- const auto negTensor = getTensor(this->problem(), negElement, negVolVars, this->fvGeometry(), negGlobalScv);
+ const auto negTensor = getT(this->problem(), negElement, negVolVars, this->fvGeometry(), negGlobalScv);
const auto& flipScvf = this->fvGeometry().flipScvf(curGlobalScvf.index(), idxInOutside);
- alpha_outside[idxInOutside] = negVolVars.extrusionFactor()
- * vtmv(flipScvf.unitOuterNormal(), negTensor, gravity);
+ alpha_outside[idxInOutside] = negVolVars.extrusionFactor()*vtmv(flipScvf.unitOuterNormal(), negTensor, gravity);
for (unsigned int pIdx = 0; pIdx < numPhases; ++pIdx)
{
@@ -572,7 +574,7 @@ public:
{
CA.umv(sum_alphas[pIdx], g[pIdx]);
- Vector AG(iv.numUnknowns());
+ FaceVector AG(iv.numUnknowns());
A.mv(sum_alphas[pIdx], AG);
// compute gravitational accelerations
@@ -588,8 +590,11 @@ public:
const auto& posLocalScv = iv.localScv(localScvIdx);
const auto& wijk = iv.omegas()[localScvfIdx][idxInOutside+1];
+ // make sure the given outside gravity container has the right size
+ assert(outsideG[pIdx][localScvfIdx].size() == iv.localScvf(localScvfIdx).neighboringLocalScvIndices().size());
+
// add contributions from all local directions
- for (LocalIndexType localDir = 0; localDir < dim; localDir++)
+ for (LocalIndexType localDir = 0; localDir < IV::Traits::GridView::dimension; localDir++)
{
// the scvf corresponding to this local direction in the scv
const auto& curLocalScvf = iv.localScvf(posLocalScv.scvfIdxLocal(localDir));
@@ -616,22 +621,29 @@ private:
* \note The matrices are expected to have been resized beforehand.
*
* \tparam IV The interaction volume type implementation
- * \tparam GetTensor Lambda to obtain the tensor w.r.t.
- * which the local system is to be solved
+ * \tparam TensorFunc Lambda to obtain the tensor w.r.t.
+ * which the local system is to be solved
*
* \param A The A matrix of the iv-local equation system
* \param B The B matrix of the iv-local equation system
* \param C The C matrix of the iv-local flux expressions
* \param D The D matrix of the iv-local flux expressions
* \param iv The mpfa-o interaction volume
- * \param getTensor Lambda to evaluate the scv-wise tensors
+ * \param getT Lambda to evaluate the scv-wise tensors
*/
- template< class IV, class GetTensor >
- void assembleLocalMatrices_(Matrix& A, Matrix& B, Matrix& C, Matrix& D,
- IV& iv, const GetTensor& getTensor)
+ template< class IV, class TensorFunc >
+ void assembleLocalMatrices_(typename IV::Traits::MatVecTraits::AMatrix& A,
+ typename IV::Traits::MatVecTraits::BMatrix& B,
+ typename IV::Traits::MatVecTraits::CMatrix& C,
+ typename IV::Traits::MatVecTraits::DMatrix& D,
+ IV& iv, const TensorFunc& getT)
{
+ using LocalIndexType = typename IV::Traits::IndexSet::LocalIndexType;
+ static constexpr int dim = IV::Traits::GridView::dimension;
+ static constexpr int dimWorld = IV::Traits::GridView::dimensionworld;
+
// Matrix D is assumed to have the right size already
- assert(D.rows() == iv.numFaces() && D.cols() == iv.numKnowns() && "Matrix D does not have the correct size");
+ assert(D.rows() == iv.numFaces() && D.cols() == iv.numKnowns());
// if only Dirichlet faces are present in the iv,
// the matrices A, B & C are undefined and D = T
@@ -652,7 +664,7 @@ private:
const auto& posGlobalScv = this->fvGeometry().scv(posLocalScv.globalScvIndex());
const auto& posVolVars = this->elemVolVars()[posGlobalScv];
const auto& posElement = iv.element(neighborScvIndices[0]);
- const auto tensor = getTensor(this->problem(), posElement, posVolVars, this->fvGeometry(), posGlobalScv);
+ const auto tensor = getT(this->problem(), posElement, posVolVars, this->fvGeometry(), posGlobalScv);
// the omega factors of the "positive" sub volume
const auto wijk = computeMpfaTransmissibility(posLocalScv, curGlobalScvf, tensor, posVolVars.extrusionFactor());
@@ -670,9 +682,9 @@ private:
else
{
// we require the matrices A,B,C to have the correct size already
- assert(A.rows() == iv.numUnknowns() && A.cols() == iv.numUnknowns() && "Matrix A does not have the correct size");
- assert(B.rows() == iv.numUnknowns() && B.cols() == iv.numKnowns() && "Matrix B does not have the correct size");
- assert(C.rows() == iv.numFaces() && C.cols() == iv.numUnknowns() && "Matrix C does not have the correct size");
+ assert(A.rows() == iv.numUnknowns() && A.cols() == iv.numUnknowns());
+ assert(B.rows() == iv.numUnknowns() && B.cols() == iv.numKnowns());
+ assert(C.rows() == iv.numFaces() && C.cols() == iv.numUnknowns());
// reset matrices
A = 0.0;
@@ -694,7 +706,7 @@ private:
const auto& posGlobalScv = this->fvGeometry().scv(posLocalScv.globalScvIndex());
const auto& posVolVars = this->elemVolVars()[posGlobalScv];
const auto& posElement = iv.element(neighborScvIndices[0]);
- const auto tensor = getTensor(this->problem(), posElement, posVolVars, this->fvGeometry(), posGlobalScv);
+ const auto tensor = getT(this->problem(), posElement, posVolVars, this->fvGeometry(), posGlobalScv);
// the omega factors of the "positive" sub volume
wijk[faceIdx][0] = computeMpfaTransmissibility(posLocalScv, curGlobalScvf, tensor, posVolVars.extrusionFactor());
@@ -740,7 +752,7 @@ private:
const auto& negGlobalScv = this->fvGeometry().scv(negLocalScv.globalScvIndex());
const auto& negVolVars = this->elemVolVars()[negGlobalScv];
const auto& negElement = iv.element( neighborScvIndices[idxOnScvf] );
- const auto negTensor = getTensor(this->problem(), negElement, negVolVars, this->fvGeometry(), negGlobalScv);
+ const auto negTensor = getT(this->problem(), negElement, negVolVars, this->fvGeometry(), negGlobalScv);
// On surface grids, use outside face for "negative" transmissibility calculation
const auto& scvf = dim < dimWorld ? this->fvGeometry().flipScvf(curGlobalScvf.index(), idxInOutside)
diff --git a/dumux/discretization/cellcentered/mpfa/omethod/staticinteractionvolume.hh b/dumux/discretization/cellcentered/mpfa/omethod/staticinteractionvolume.hh
index 2304ac87383f7a48bcff8af92db2d87b91843bbf..7456191a3568c08f0be1eab130904c03ca59992c 100644
--- a/dumux/discretization/cellcentered/mpfa/omethod/staticinteractionvolume.hh
+++ b/dumux/discretization/cellcentered/mpfa/omethod/staticinteractionvolume.hh
@@ -30,15 +30,13 @@
#include <dune/common/fvector.hh>
#include <dumux/common/math.hh>
-#include <dumux/common/properties.hh>
#include <dumux/common/matrixvectorhelper.hh>
-#include <dumux/discretization/cellcentered/mpfa/interactionvolumedatahandle.hh>
#include <dumux/discretization/cellcentered/mpfa/interactionvolumebase.hh>
-#include <dumux/discretization/cellcentered/mpfa/omethod/interactionvolume.hh>
#include <dumux/discretization/cellcentered/mpfa/dualgridindexset.hh>
#include <dumux/discretization/cellcentered/mpfa/localfacedata.hh>
#include <dumux/discretization/cellcentered/mpfa/methods.hh>
+#include <dumux/discretization/cellcentered/mpfa/omethod/interactionvolume.hh>
#include "localsubcontrolentities.hh"
#include "interactionvolumeindexset.hh"
@@ -46,47 +44,59 @@
namespace Dumux
{
//! Forward declaration of the o-method's static interaction volume
-template< class Traits, int localSize > class CCMpfaOStaticInteractionVolume;
+template< class Traits > class CCMpfaOStaticInteractionVolume;
-//! Specialization of the default interaction volume traits class for the mpfa-o method
-template< class TypeTag, int localSize >
+/*!
+ * \ingroup CCMpfaDiscretization
+ * \brief The default interaction volume traits class for the mpfa-o method
+ * with known size of the interaction volumes at compile time. It uses
+ * statically sized containers for the iv-local data structures and static
+ * matrices and vectors.
+ *
+ * \tparam NI The type used for the dual grid's nodal index sets
+ * \tparam S The Type used for scalar values
+ * \tparam C The number of sub-control volumes (cells) in the ivs
+ * \tparam F The number of sub-control volume faces in the ivs
+ */
+template< class NI, class S, int C, int F >
struct CCMpfaODefaultStaticInteractionVolumeTraits
{
private:
- 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;
+ using GridIndexType = typename NI::GridIndexType;
+ using LocalIndexType = typename NI::LocalIndexType;
+
+ static constexpr int dim = NI::Traits::GridView::dimension;
+ static constexpr int dimWorld = NI::Traits::GridView::dimensionworld;
+
+ //! Matrix/Vector traits to be used by the data handle
+ struct MVTraits
+ {
+ using AMatrix = Dune::FieldMatrix< S, F, F >;
+ using BMatrix = Dune::FieldMatrix< S, F, C >;
+ using CMatrix = Dune::FieldMatrix< S, F, F >;
+ using DMatrix = Dune::FieldMatrix< S, F, C >;
+ using TMatrix = Dune::FieldMatrix< S, F, C >;
+ using CellVector = Dune::FieldVector< S, C >;
+ using FaceVector = Dune::FieldVector< S, F >;
+ };
public:
- //! export the problem type (needed for iv-local assembly)
- using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
- //! export the type of the local view on the finite volume grid geometry
- using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry)::LocalView;
- //! export the type of the local view on the grid volume variables
- using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);
//! export the type of grid view
- using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
- //! export the type used for scalar values
- using ScalarType = typename GET_PROP_TYPE(TypeTag, Scalar);
- //! 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 = typename NodalIndexSet::LocalIndexType;
+ using GridView = typename NI::Traits::GridView;
//! export the type for the interaction volume index set
- using IndexSet = CCMpfaOInteractionVolumeIndexSet< NodalIndexSet >;
+ using IndexSet = CCMpfaOInteractionVolumeIndexSet< NI >;
//! export the type of interaction-volume local scvs
- using LocalScvType = CCMpfaOInteractionVolumeLocalScv< IndexSet, ScalarType, dim, dimWorld >;
+ using LocalScvType = CCMpfaOInteractionVolumeLocalScv< IndexSet, S, dim, dimWorld >;
//! export the type of interaction-volume local scvfs
using LocalScvfType = CCMpfaOInteractionVolumeLocalScvf< IndexSet >;
//! export the type of used for the iv-local face data
using LocalFaceData = InteractionVolumeLocalFaceData<GridIndexType, LocalIndexType>;
- //! export the type used for iv-local matrices
- using Matrix = Dune::FieldMatrix< ScalarType, localSize, localSize >;
- //! export the type used for iv-local vectors
- using Vector = Dune::FieldVector< ScalarType, localSize >;
- //! export the data handle type for this iv
- using DataHandle = InteractionVolumeDataHandle< TypeTag, Matrix, Vector, LocalIndexType >;
+ //! export the matrix/vector traits to be used by the iv
+ using MatVecTraits = MVTraits;
+ //! export the number of scvs in the interaction volumes
+ static constexpr int numScvs = C;
+ //! export the number of scvfs in the interaction volumes
+ static constexpr int numScvfs = F;
};
/*!
@@ -98,182 +108,180 @@ public:
* vectors/matrices defined in the traits class in case static types are used.
*
* \tparam Traits The type traits class to be used
- * \tparam localSize The size of the local eq system
- * This also is the number of local scvs/scvfs
*/
-template< class Traits, int localSize >
-class CCMpfaOStaticInteractionVolume : public CCMpfaInteractionVolumeBase< CCMpfaOInteractionVolume<Traits>, Traits >
+template< class Traits >
+class CCMpfaOStaticInteractionVolume
+ : public CCMpfaInteractionVolumeBase< CCMpfaOStaticInteractionVolume<Traits>, Traits >
{
- using ParentType = CCMpfaInteractionVolumeBase< CCMpfaOInteractionVolume<Traits>, Traits >;
-
- using Scalar = typename Traits::ScalarType;
- using Helper = typename Traits::MpfaHelper;
- using Problem = typename Traits::Problem;
- using FVElementGeometry = typename Traits::FVElementGeometry;
-
using GridView = typename Traits::GridView;
using Element = typename GridView::template Codim<0>::Entity;
+ using IndexSet = typename Traits::IndexSet;
+ using GridIndexType = typename IndexSet::GridIndexType;
+ using LocalIndexType = typename IndexSet::LocalIndexType;
+ using Stencil = typename IndexSet::NodalGridStencilType;
+
static constexpr int dim = GridView::dimension;
+ static constexpr int dimWorld = GridView::dimensionworld;
+
+ //! export scalar type from T matrix and define omegas
+ using Scalar = typename Traits::MatVecTraits::TMatrix::value_type;
using DimVector = Dune::FieldVector<Scalar, dim>;
+ using FaceOmegas = typename Dune::ReservedVector<DimVector, 2>;
+
+ using AMatrix = typename Traits::MatVecTraits::AMatrix;
+ using BMatrix = typename Traits::MatVecTraits::BMatrix;
+ using CMatrix = typename Traits::MatVecTraits::CMatrix;
- using Matrix = typename Traits::Matrix;
using LocalScvType = typename Traits::LocalScvType;
using LocalScvfType = typename Traits::LocalScvfType;
using LocalFaceData = typename Traits::LocalFaceData;
- using IndexSet = typename Traits::IndexSet;
- using GridIndexType = typename GridView::IndexSet::IndexType;
- using LocalIndexType = typename Traits::LocalIndexType;
- using Stencil = typename IndexSet::GridStencilType;
-
- //! Data attached to scvf touching Dirichlet boundaries.
- //! For the default o-scheme, we only store the corresponding vol vars index.
- class DirichletData
- {
- GridIndexType volVarIndex_;
+ static constexpr int numScvf = Traits::numScvfs;
+ static constexpr int numScv = Traits::numScvs;
- public:
- //! Constructor
- DirichletData(const GridIndexType index) : volVarIndex_(index) {}
+public:
+ //! export the standard o-methods dirichlet data
+ using DirichletData = typename CCMpfaOInteractionVolume< Traits >::DirichletData;
- //! Return corresponding vol var index
- GridIndexType volVarIndex() const { return volVarIndex_; }
- };
+ //! export the type used for transmissibility storage
+ using TransmissibilityStorage = std::array< FaceOmegas, numScvf >;
-public:
//! publicly state the mpfa-scheme this interaction volume is associated with
static constexpr MpfaMethods MpfaMethod = MpfaMethods::oMethod;
//! Sets up the local scope for a given iv index set
+ template< class Problem, class FVElementGeometry >
void setUpLocalScope(const IndexSet& indexSet,
const Problem& problem,
const FVElementGeometry& fvGeometry)
{
// for the o-scheme, the stencil is equal to the scv
// index set of the dual grid's nodal index set
+ assert(indexSet.numScvs() == numScv);
stencil_ = &indexSet.nodalIndexSet().globalScvIndices();
- // set up local geometry containers
+ // set up stuff related to sub-control volumes
const auto& scvIndices = indexSet.globalScvIndices();
- for (LocalIndexType localIdx = 0; localIdx < localSize; localIdx++)
+ for (LocalIndexType scvIdxLocal = 0; scvIdxLocal < numScv; scvIdxLocal++)
{
- // scv-related quantities
- scvs_[localIdx] = LocalScvType(Helper(),
+ elements_[scvIdxLocal] = fvGeometry.fvGridGeometry().element( scvIndices[scvIdxLocal] );
+ scvs_[scvIdxLocal] = LocalScvType(fvGeometry.fvGridGeometry().mpfaHelper(),
fvGeometry,
- fvGeometry.scv( scvIndices[localIdx] ),
- localIdx,
+ fvGeometry.scv( scvIndices[scvIdxLocal] ),
+ scvIdxLocal,
indexSet);
- elements_[localIdx] = fvGeometry.fvGridGeometry().element( scvIndices[localIdx] );
+ }
- // scvf-related quantities
- const auto& scvf = fvGeometry.scvf(indexSet.scvfIdxGlobal(localIdx));
+ // set up quantitites related to sub-control volume faces
+ for (LocalIndexType faceIdxLocal = 0; faceIdxLocal < numScvf; ++faceIdxLocal)
+ {
+ const auto& scvf = fvGeometry.scvf(indexSet.scvfIdxGlobal(faceIdxLocal));
+
+ // the neighboring scvs in local indices (order: 0 - inside scv, 1..n - outside scvs)
+ const auto& neighborScvIndicesLocal = indexSet.neighboringLocalScvIndices(faceIdxLocal);
- // this interaction volume implementation does not work on boundaries
- assert(!scvf.boundary() && "static mpfa-o interaction volume cannot be used on boundaries");
+ // this does not work for network grids or on boundaries
+ assert(!scvf.boundary());
+ assert(neighborScvIndicesLocal.size() == 2);
- const auto& neighborScvIndicesLocal = indexSet.neighboringLocalScvIndices(localIdx);
- scvfs_[localIdx] = LocalScvfType(scvf, neighborScvIndicesLocal, localIdx, /*isDirichlet*/false);
- localFaceData_[localIdx*2] = LocalFaceData(localIdx, neighborScvIndicesLocal[0], scvf.index());
+ // create iv-local scvf objects
+ scvfs_[faceIdxLocal] = LocalScvfType(scvf, neighborScvIndicesLocal, faceIdxLocal, /*isDirichlet*/false);
+ localFaceData_[faceIdxLocal*2] = LocalFaceData(faceIdxLocal, neighborScvIndicesLocal[0], scvf.index());
+ // add local face data objects for the outside face
const auto outsideLocalScvIdx = neighborScvIndicesLocal[1];
- // loop over scvfs in outside scv until we find the one coinciding with current scvf
for (int coord = 0; coord < dim; ++coord)
{
- if (indexSet.scvfIdxLocal(outsideLocalScvIdx, coord) == localIdx)
+ if (indexSet.scvfIdxLocal(outsideLocalScvIdx, coord) == faceIdxLocal)
{
const auto globalScvfIdx = indexSet.nodalIndexSet().scvfIdxGlobal(outsideLocalScvIdx, coord);
const auto& flipScvf = fvGeometry.scvf(globalScvfIdx);
- localFaceData_[localIdx*2 + 1] = LocalFaceData(localIdx, // iv-local scvf idx
- outsideLocalScvIdx, // iv-local scv index
- 0, // scvf-local index in outside faces
- flipScvf.index()); // global scvf index
+ localFaceData_[faceIdxLocal*2+1] = LocalFaceData(faceIdxLocal, // iv-local scvf idx
+ outsideLocalScvIdx, // iv-local scv index
+ 0, // scvf-local index in outside faces
+ flipScvf.index()); // global scvf index
+ break; // go to next outside face
}
}
+
+ // make sure we found it
+ assert(localFaceData_[faceIdxLocal*2+1].ivLocalInsideScvIndex() == outsideLocalScvIdx);
}
// Maybe resize local matrices if dynamic types are used
- resizeMatrix(A_, localSize, localSize);
- resizeMatrix(B_, localSize, localSize);
- resizeMatrix(C_, localSize, localSize);
+ resizeMatrix(A_, numScvf, numScvf);
+ resizeMatrix(B_, numScvf, numScv);
+ resizeMatrix(C_, numScvf, numScv);
}
//! returns the number of primary scvfs of this interaction volume
- static constexpr std::size_t numFaces() { return localSize; }
+ static constexpr std::size_t numFaces() { return numScvf; }
//! returns the number of intermediate unknowns within this interaction volume
- static constexpr std::size_t numUnknowns() { return localSize; }
+ static constexpr std::size_t numUnknowns() { return numScvf; }
//! returns the number of (in this context) known solution values within this interaction volume
- static constexpr std::size_t numKnowns() { return localSize; }
+ static constexpr std::size_t numKnowns() { return numScv; }
//! returns the number of scvs embedded in this interaction volume
- static constexpr std::size_t numScvs() { return localSize; }
+ static constexpr std::size_t numScvs() { return numScv; }
//! returns the cell-stencil of this interaction volume
const Stencil& stencil() const { return *stencil_; }
//! returns the grid element corresponding to a given iv-local scv idx
- const Element& element(const LocalIndexType ivLocalScvIdx) const
- {
- assert(ivLocalScvIdx < localSize);
- return elements_[ivLocalScvIdx];
- }
+ const Element& element(LocalIndexType ivLocalScvIdx) const { return elements_[ivLocalScvIdx]; }
//! returns the local scvf entity corresponding to a given iv-local scvf idx
- const LocalScvfType& localScvf(const LocalIndexType ivLocalScvfIdx) const
- {
- assert(ivLocalScvfIdx < localSize);
- return scvfs_[ivLocalScvfIdx];
- }
+ const LocalScvfType& localScvf(LocalIndexType ivLocalScvfIdx) const { return scvfs_[ivLocalScvfIdx]; }
//! returns the local scv entity corresponding to a given iv-local scv idx
- const LocalScvType& localScv(const LocalIndexType ivLocalScvIdx) const
- {
- assert(ivLocalScvIdx < localSize);
- return scvs_[ivLocalScvIdx];
- }
+ const LocalScvType& localScv(LocalIndexType ivLocalScvIdx) const { return scvs_[ivLocalScvIdx]; }
//! returns a reference to the container with the local face data
- const std::array<LocalFaceData, 2*localSize>& localFaceData() const { return localFaceData_; }
+ const std::array<LocalFaceData, numScvf*2>& localFaceData() const { return localFaceData_; }
//! Returns a reference to the information container on Dirichlet BCs within this iv.
//! Here, we return an empty container as this implementation cannot be used on boundaries.
const std::array<DirichletData, 0>& dirichletData() const { return dirichletData_; }
//! returns the matrix associated with face unknowns in local equation system
- const Matrix& A() const { return A_; }
- Matrix& A() { return A_; }
+ const AMatrix& A() const { return A_; }
+ AMatrix& A() { return A_; }
//! returns the matrix associated with cell unknowns in local equation system
- const Matrix& B() const { return B_; }
- Matrix& B() { return B_; }
+ const BMatrix& B() const { return B_; }
+ BMatrix& B() { return B_; }
//! returns the matrix associated with face unknowns in flux expressions
- const Matrix& C() const { return C_; }
- Matrix& C() { return C_; }
+ const CMatrix& C() const { return C_; }
+ CMatrix& C() { return C_; }
//! returns container storing the transmissibilities for each face & coordinate
- const std::array< std::array<DimVector, 2>, localSize >& omegas() const { return wijk_; }
- std::array< std::array<DimVector, 2>, localSize >& omegas() { return wijk_; }
+ const TransmissibilityStorage& omegas() const { return wijk_; }
+ TransmissibilityStorage& omegas() { return wijk_; }
//! returns the number of interaction volumes living around a vertex
- //! the mpfa-o scheme always constructs one iv per vertex
- template< class NodalIndexSet >
- static std::size_t numInteractionVolumesAtVertex(const NodalIndexSet& nodalIndexSet) { return 1; }
+ template< class NI >
+ static constexpr std::size_t numIVAtVertex(const NI& nodalIndexSet) { return 1; }
//! adds the iv index sets living around a vertex to a given container
//! and stores the the corresponding index in a map for each scvf
- template<class IvIndexSetContainer, class ScvfIndexMap, class NodalIndexSet>
- static void addInteractionVolumeIndexSets(IvIndexSetContainer& ivIndexSetContainer,
- ScvfIndexMap& scvfIndexMap,
- const NodalIndexSet& nodalIndexSet)
+ template< class IvIndexSetContainer,
+ class ScvfIndexMap,
+ class NodalIndexSet,
+ class FlipScvfIndexSet >
+ static void addIVIndexSets(IvIndexSetContainer& ivIndexSetContainer,
+ ScvfIndexMap& scvfIndexMap,
+ const NodalIndexSet& nodalIndexSet,
+ const FlipScvfIndexSet& flipScvfIndexSet)
{
// the global index of the iv index set that is about to be created
const auto curGlobalIndex = ivIndexSetContainer.size();
// make the one index set for this node
- ivIndexSetContainer.emplace_back(nodalIndexSet);
+ ivIndexSetContainer.emplace_back(nodalIndexSet, flipScvfIndexSet);
// store the index mapping
for (const auto scvfIdx : nodalIndexSet.globalScvfIndices())
@@ -285,22 +293,21 @@ private:
const Stencil* stencil_;
// Variables defining the local scope
- std::array<Element, localSize> elements_;
- std::array<LocalScvType, localSize> scvs_;
- std::array<LocalScvfType, localSize> scvfs_;
- std::array<LocalFaceData, 2*localSize> localFaceData_;
-
- // Empty Dirichlet container to be compatible with dynamic assembly
- std::array<DirichletData, 0> dirichletData_;
+ std::array<Element, numScv> elements_;
+ std::array<LocalScvType, numScv> scvs_;
+ std::array<LocalScvfType, numScvf> scvfs_;
+ std::array<LocalFaceData, numScvf*2> localFaceData_;
// Matrices needed for computation of transmissibilities
- Matrix A_;
- Matrix B_;
- Matrix C_;
+ AMatrix A_;
+ BMatrix B_;
+ CMatrix C_;
// The omega factors are stored during assembly of local system
- // we assume one "outside" face per scvf (does not work on surface grids)
- std::array< std::array<DimVector, 2>, localSize > wijk_;
+ TransmissibilityStorage wijk_;
+
+ // Dummy dirichlet data container (compatibility with dynamic o-iv)
+ std::array<DirichletData, 0> dirichletData_;
};
} // end namespace
diff --git a/dumux/discretization/cellcentered/mpfa/properties.hh b/dumux/discretization/cellcentered/mpfa/properties.hh
index 06f3d8704a822edb5cd882f1d935ea00745607e5..3f9cba7bc470264286d22f48471dacbaa3a46c42 100644
--- a/dumux/discretization/cellcentered/mpfa/properties.hh
+++ b/dumux/discretization/cellcentered/mpfa/properties.hh
@@ -25,6 +25,8 @@
#ifndef DUMUX_CC_MPFA_PROPERTIES_HH
#define DUMUX_CC_MPFA_PROPERTIES_HH
+#include <dune/common/reservedvector.hh>
+
#include <dumux/common/properties.hh>
#include <dumux/common/defaultmappertraits.hh>
@@ -64,63 +66,49 @@ SET_PROP(CCMpfaModel, DiscretizationMethod)
static const DiscretizationMethods value = DiscretizationMethods::CCMpfa;
};
-//! Set the maximum admissible number of branches per scvf
-SET_PROP(CCMpfaModel, MaxNumNeighborsPerScvf)
-{
-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 constexpr std::size_t value = dim < dimWorld ? 9 : 2;
-};
-
//! Set the index set type used on the dual grid nodes
SET_PROP(CCMpfaModel, DualGridNodalIndexSet)
{
+ using GV = typename GET_PROP_TYPE(TypeTag, GridView);
+ static constexpr int dim = GV::dimension;
+ static constexpr int dimWorld = GV::dimensionworld;
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, MaxNumNeighborsPerScvf);
-
- // 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;
-
+ struct Traits
+ {
+ using GridView = GV;
+ using GridIndexType = typename GV::IndexSet::IndexType;
+ using LocalIndexType = std::uint8_t;
+
+ //! per default, we use dynamic data containers (iv size unknown)
+ template< class T > using NodalScvDataStorage = std::vector< T >;
+ template< class T > using NodalScvfDataStorage = std::vector< T >;
+
+ //! store data on neighbors of scvfs in static containers if possible
+ template< class T >
+ using ScvfNeighborDataStorage = typename std::conditional_t< (dim<dimWorld),
+ std::vector< T >,
+ Dune::ReservedVector< T, 2 > >;
+ };
public:
- using type = CCMpfaDualGridNodalIndexSet<GI, LI, dim, maxE, maxB>;
+ using type = CCMpfaDualGridNodalIndexSet< Traits >;
};
-//! The mpfa helper class
-SET_TYPE_PROP(CCMpfaModel, MpfaHelper, CCMpfaHelper<TypeTag>);
-
//! Per default, we use the dynamic mpfa-o interaction volume
SET_PROP(CCMpfaModel, PrimaryInteractionVolume)
{
-private:
- //! use the default traits
- using Traits = CCMpfaODefaultInteractionVolumeTraits< TypeTag >;
public:
- using type = CCMpfaOInteractionVolume< Traits >;
+ using type = typename GET_PROP_TYPE(TypeTag, SecondaryInteractionVolume);
};
-//! Per default, we use the dynamic mpfa-o interaction volume as secondary type
+//! Per default, we use the dynamic mpfa-o interaction volume on boundaries
SET_PROP(CCMpfaModel, SecondaryInteractionVolume)
{
private:
- //! use the default traits
- using Traits = CCMpfaODefaultInteractionVolumeTraits< TypeTag >;
+ using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+ using NodalIndexSet = typename GET_PROP_TYPE(TypeTag, DualGridNodalIndexSet);
+
+ // use the default traits
+ using Traits = CCMpfaODefaultInteractionVolumeTraits< NodalIndexSet, Scalar >;
public:
using type = CCMpfaOInteractionVolume< Traits >;
};
@@ -133,30 +121,27 @@ private:
using PrimaryIV = typename GET_PROP_TYPE(TypeTag, PrimaryInteractionVolume);
using SecondaryIV = typename GET_PROP_TYPE(TypeTag, SecondaryInteractionVolume);
- static constexpr bool enableCache = GET_PROP_VALUE(TypeTag, EnableFVGridGeometryCache);
-
struct Traits : public DefaultMapperTraits<GridView>
{
using SubControlVolume = CCSubControlVolume<GridView>;
using SubControlVolumeFace = CCMpfaSubControlVolumeFace<GridView>;
- using MpfaHelper = typename GET_PROP_TYPE(TypeTag, MpfaHelper);
using NodalIndexSet = typename GET_PROP_TYPE(TypeTag, DualGridNodalIndexSet);
- template< class FVGridGeometry >
- using GridIvIndexSets = CCMpfaGridInteractionVolumeIndexSets< FVGridGeometry,
- NodalIndexSet,
- PrimaryIV,
- SecondaryIV >;
+ //! type definitions depending on the FVGridGeometry itself
+ template< class FVGridGeom >
+ using MpfaHelper = CCMpfaHelper< FVGridGeom >;
+
+ template< class FVGridGeom >
+ using ConnectivityMap = CCMpfaConnectivityMap<FVGridGeom, FVGridGeom::GridIVIndexSets::PrimaryInteractionVolume::MpfaMethod>;
- template< class FVGridGeometry, bool enableCache >
- using LocalView = CCMpfaFVElementGeometry<FVGridGeometry, enableCache>;
+ template< class FVGridGeom >
+ using GridIvIndexSets = CCMpfaGridInteractionVolumeIndexSets< FVGridGeom, NodalIndexSet, PrimaryIV, SecondaryIV >;
- //! Per default, we use the o-method and thus the simple assembly map
- template< class FVGridGeometry >
- using ConnectivityMap = CCMpfaConnectivityMap<FVGridGeometry, FVGridGeometry::GridIVIndexSets::PrimaryInteractionVolume::MpfaMethod>;
+ template< class FVGridGeom, bool enableCache >
+ using LocalView = CCMpfaFVElementGeometry<FVGridGeom, enableCache>;
};
public:
- using type = CCMpfaFVGridGeometry<GridView, Traits, enableCache>;
+ using type = CCMpfaFVGridGeometry<GridView, Traits, GET_PROP_VALUE(TypeTag, EnableFVGridGeometryCache)>;
};
//! The global flux variables cache vector class
diff --git a/dumux/discretization/cellcentered/tpfa/properties.hh b/dumux/discretization/cellcentered/tpfa/properties.hh
index c0efc9cdb342547361e7d01e8e3a7b79048c2a28..8e9cb32413add7152cf826cd722ceb885fb6fa92 100644
--- a/dumux/discretization/cellcentered/tpfa/properties.hh
+++ b/dumux/discretization/cellcentered/tpfa/properties.hh
@@ -85,19 +85,6 @@ SET_TYPE_PROP(CCTpfaModel, ElementFluxVariablesCache, CCTpfaElementFluxVariables
//! The global current volume variables vector class
SET_TYPE_PROP(CCTpfaModel, GridVolumeVariables, CCGridVolumeVariables<TypeTag, GET_PROP_VALUE(TypeTag, EnableGridVolumeVariablesCache)>);
-//! Set the maximum admissible number of branches per scvf
-SET_PROP(CCTpfaModel, MaxNumNeighborsPerScvf)
-{
-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 constexpr std::size_t value = dim < dimWorld ? 9 : 2;
-};
-
//! Set the solution vector type for an element
SET_TYPE_PROP(CCTpfaModel, ElementSolutionVector, CCElementSolution<TypeTag>);
diff --git a/dumux/discretization/fluxstencil.hh b/dumux/discretization/fluxstencil.hh
index e34321cf03a10707a3a534099c732a2af21cf89b..b249dbe5a8afce8670470f82d633433a44024dc5 100644
--- a/dumux/discretization/fluxstencil.hh
+++ b/dumux/discretization/fluxstencil.hh
@@ -104,7 +104,7 @@ public:
using ScvfStencilIForJ = std::vector<IndexType>;
//! The flux stencil type
- using Stencil = typename NodalIndexSet::GridStencilType;
+ using Stencil = typename NodalIndexSet::NodalGridStencilType;
//! Returns the indices of the elements required for flux calculation on an scvf.
static const Stencil& stencil(const Element& element,
diff --git a/dumux/porousmediumflow/fluxvariablescache.hh b/dumux/porousmediumflow/fluxvariablescache.hh
index 57662df20fb9661fd1d5cc89f29f4df87e46378a..f775ba98a54191ed594749a48b5f2ac79da3b0a7 100644
--- a/dumux/porousmediumflow/fluxvariablescache.hh
+++ b/dumux/porousmediumflow/fluxvariablescache.hh
@@ -152,7 +152,7 @@ class PorousMediumFluxVariablesCacheImplementation<TypeTag, DiscretizationMethod
{
using GridIndexType = typename GET_PROP_TYPE(TypeTag, GridView)::IndexSet::IndexType;
- using MpfaHelper = typename GET_PROP_TYPE(TypeTag, MpfaHelper);
+ using MpfaHelper = typename GET_PROP_TYPE(TypeTag, FVGridGeometry)::MpfaHelper;
static constexpr bool considerSecondary = MpfaHelper::considerSecondaryIVs();
public:
//! Returns whether or not this cache has been updated