From b8fcee85c3be1ef6ae8ff94e01e62e02fa2b2ab3 Mon Sep 17 00:00:00 2001
From: DennisGlaeser <dennis.glaeser@iws.uni-stuttgart.de>
Date: Thu, 1 Feb 2018 13:44:07 +0100
Subject: [PATCH] [mpfa-o][iv] make static iv more flexible
The traits class now has to contain both the number of scvs and scvfs
in the ivs. This allows the applicability in 3d where numScvs != numScvfs.
---
.../mpfa/omethod/staticinteractionvolume.hh | 262 +++++++++---------
1 file changed, 136 insertions(+), 126 deletions(-)
diff --git a/dumux/discretization/cellcentered/mpfa/omethod/staticinteractionvolume.hh b/dumux/discretization/cellcentered/mpfa/omethod/staticinteractionvolume.hh
index 2304ac8738..ae77cdb463 100644
--- a/dumux/discretization/cellcentered/mpfa/omethod/staticinteractionvolume.hh
+++ b/dumux/discretization/cellcentered/mpfa/omethod/staticinteractionvolume.hh
@@ -30,15 +30,14 @@
#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 +45,64 @@
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 PT Traits class encapsulating info on the physical processes
+ * \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, class PT, 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::GridView::dimension;
+ static constexpr int dimWorld = NI::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);
+ using GridView = typename NI::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 ScalarType = S;
//! 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 >;
//! 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 matrix/vector traits to be used by the iv
+ using MatVecTraits = MVTraits;
//! export the data handle type for this iv
- using DataHandle = InteractionVolumeDataHandle< TypeTag, Matrix, Vector, LocalIndexType >;
+ using DataHandle = InteractionVolumeDataHandle< MatVecTraits, PT >;
+ //! 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 +114,177 @@ 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::GridStencilType;
+
static constexpr int dim = GridView::dimension;
- using DimVector = Dune::FieldVector<Scalar, dim>;
+ static constexpr int dimWorld = GridView::dimensionworld;
+ using DimVector = Dune::FieldVector<typename Traits::ScalarType, 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;
+ static constexpr int numScvf = Traits::numScvfs;
+ static constexpr int numScv = Traits::numScvs;
- //! 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) {}
+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 numInteractionVolumesAtVertex(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>
+ template< class IvIndexSetContainer,
+ class ScvfIndexMap,
+ class NodalIndexSet,
+ class FlipScvfIndexSet >
static void addInteractionVolumeIndexSets(IvIndexSetContainer& ivIndexSetContainer,
ScvfIndexMap& scvfIndexMap,
- const NodalIndexSet& nodalIndexSet)
+ 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 +296,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
--
GitLab