diff --git a/dumux/discretization/cellcentered/mpfa/darcyslaw.hh b/dumux/discretization/cellcentered/mpfa/darcyslaw.hh
index 4f5f9d79b84b728f9aa1d5f44b2a0fe496e8513d..50335bfd9628b24839945db19a9def02c72c3326 100644
--- a/dumux/discretization/cellcentered/mpfa/darcyslaw.hh
+++ b/dumux/discretization/cellcentered/mpfa/darcyslaw.hh
@@ -107,16 +107,14 @@ class DarcysLawImplementation<TypeTag, DiscretizationMethods::CCMpfa>
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 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 SecondaryIvCellVector = typename SecondaryInteractionVolume::Traits::MatVecTraits::CellVector;
+ using SecondaryIvTij = typename SecondaryInteractionVolume::Traits::MatVecTraits::TMatrix::row_type;
public:
// export the filler type
@@ -236,11 +234,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]; }
@@ -262,8 +260,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 84aca7ed72008982cba4e459729476ca3919caa1..b9b95cbc5ffd6da967f48d40601b9b77f41de4d6 100644
--- a/dumux/discretization/cellcentered/mpfa/dualgridindexset.hh
+++ b/dumux/discretization/cellcentered/mpfa/dualgridindexset.hh
@@ -38,7 +38,7 @@ namespace Dumux
* \brief Nodal index set for mpfa schemes, constructed
* around grid vertices.
*
- * \tparam GI The type used for indices on the grid
+ * \tparam GV The grid view type
* \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.
@@ -46,12 +46,16 @@ namespace Dumux
* This is only to be specified for network grids and defaults to 1
* for normal grids.
*/
-template< class GI, class LI, int dim, int maxE, int maxB = 2 >
+template< class GV, class LI, int dim, int maxE, int maxB = 2 >
class CCMpfaDualGridNodalIndexSet
{
+ using GI = typename GV::IndexSet::IndexType;
using DimIndexVector = Dune::ReservedVector<GI, dim>;
public:
+ //! Export the grid view type
+ using GridView = GV;
+
//! Export the used index types
using GridIndexType = GI;
using LocalIndexType = LI;
diff --git a/dumux/discretization/cellcentered/mpfa/fickslaw.hh b/dumux/discretization/cellcentered/mpfa/fickslaw.hh
index ae3979d58deb66b2702ce536db2a6ae3cda9d59a..3f68da0dde36f2484cc5ce41f81a90690edccb13 100644
--- a/dumux/discretization/cellcentered/mpfa/fickslaw.hh
+++ b/dumux/discretization/cellcentered/mpfa/fickslaw.hh
@@ -105,16 +105,14 @@ class FicksLawImplementation<TypeTag, DiscretizationMethods::CCMpfa>
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 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 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 +196,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 +215,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 8b0863abf7f9375255d8a4ea22212f02d5a3b871..2f4735da25a26c47143f5b33f3a7fe1679bc0e03 100644
--- a/dumux/discretization/cellcentered/mpfa/fluxvariablescachefiller.hh
+++ b/dumux/discretization/cellcentered/mpfa/fluxvariablescachefiller.hh
@@ -431,8 +431,8 @@ private:
using LambdaFactory = TensorLambdaFactory<TypeTag, DiscretizationMethods::CCMpfa>;
// get instance of the interaction volume-local assembler
- using IVTraits = typename InteractionVolume::Traits;
- using IvLocalAssembler = InteractionVolumeAssembler< IVTraits, InteractionVolume::MpfaMethod >;
+ 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
@@ -516,8 +516,8 @@ private:
using LambdaFactory = TensorLambdaFactory<TypeTag, DiscretizationMethods::CCMpfa>;
// get instance of the interaction volume-local assembler
- using IVTraits = typename InteractionVolume::Traits;
- using IvLocalAssembler = InteractionVolumeAssembler< IVTraits, InteractionVolume::MpfaMethod >;
+ 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
@@ -555,8 +555,8 @@ private:
using LambdaFactory = TensorLambdaFactory<TypeTag, DiscretizationMethods::CCMpfa>;
// get instance of the interaction volume-local assembler
- using IVTraits = typename InteractionVolume::Traits;
- using IvLocalAssembler = InteractionVolumeAssembler< IVTraits, InteractionVolume::MpfaMethod >;
+ static constexpr MpfaMethods M = InteractionVolume::MpfaMethod;
+ using IvLocalAssembler = InteractionVolumeAssembler< Problem, FVElementGeometry, ElementVolumeVariables, M >;
IvLocalAssembler localAssembler(problem(), fvGeometry(), elemVolVars());
if (forceUpdateAll || soldependentAdvection)
diff --git a/dumux/discretization/cellcentered/mpfa/fourierslaw.hh b/dumux/discretization/cellcentered/mpfa/fourierslaw.hh
index 56e3af6da2d0625c0aca32d27c8afbec6f26cb41..121de6a44c68bc791da278a8c49a626b231305a4 100644
--- a/dumux/discretization/cellcentered/mpfa/fourierslaw.hh
+++ b/dumux/discretization/cellcentered/mpfa/fourierslaw.hh
@@ -103,16 +103,14 @@ class FouriersLawImplementation<TypeTag, DiscretizationMethods::CCMpfa>
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 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 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 +186,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 +208,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/fvgridgeometry.hh b/dumux/discretization/cellcentered/mpfa/fvgridgeometry.hh
index a23d1674f7608d8b33e8531a429982d74023ac3a..ef2380333056368c0cbe57402df0934cf4abd0a4 100644
--- a/dumux/discretization/cellcentered/mpfa/fvgridgeometry.hh
+++ b/dumux/discretization/cellcentered/mpfa/fvgridgeometry.hh
@@ -346,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]; }
@@ -689,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]; }
diff --git a/dumux/discretization/cellcentered/mpfa/gridinteractionvolumeindexsets.hh b/dumux/discretization/cellcentered/mpfa/gridinteractionvolumeindexsets.hh
index 358e0b0fe1a2c7ec93329861d41835245293a866..9dea0f375c17ea722c5abdd68b0419e1f67e6950 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"
diff --git a/dumux/discretization/cellcentered/mpfa/helper.hh b/dumux/discretization/cellcentered/mpfa/helper.hh
index 9c7f47943e1a8a3bf2af30531446c9434156a1db..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;
@@ -457,9 +442,8 @@ public:
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);
}
}
@@ -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 PrimaryInteractionVolume = typename FVGridGeometry::GridIVIndexSets::PrimaryInteractionVolume;
+ using SecondaryInteractionVolume = typename FVGridGeometry::GridIVIndexSets::SecondaryInteractionVolume;
- 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 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 a4600b12b58ae2db93633b5c1936d066910cdbab..d55a189c3816ff7a69bfa918ea232372ce2ab572 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,11 @@ 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,10 +49,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 matrix/vector type traits to be used by the iv
+ * using MatVecTraits = ...;
* //! export the data handle type for this iv
* using DataHandle = ...;
* \endcode
@@ -78,25 +64,26 @@ 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 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)
@@ -124,16 +111,13 @@ public:
const typename Traits::IndexSet::GridStencilType& 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 >
diff --git a/dumux/discretization/cellcentered/mpfa/interactionvolumedatahandle.hh b/dumux/discretization/cellcentered/mpfa/interactionvolumedatahandle.hh
index 025154b7f9a042f990e029b7d2a61ddc6c731482..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,20 +201,17 @@ public:
{
for (unsigned int cIdx = 0; cIdx < numComponents; ++cIdx)
{
- using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
- if (cIdx == FluidSystem::getMainComponent(pIdx))
- 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;
@@ -230,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_]; }
@@ -243,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
@@ -270,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;
@@ -286,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 4629683fd96bf1c6e99e391cc92df4dd543bed6a..534291070f5c649bdd99b5dbeff4e25d4579c7e6 100644
--- a/dumux/discretization/cellcentered/mpfa/omethod/interactionvolume.hh
+++ b/dumux/discretization/cellcentered/mpfa/omethod/interactionvolume.hh
@@ -24,12 +24,14 @@
#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>
@@ -46,31 +48,44 @@ 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
+ * \tparam PhysicsTraits Traits class encapsulating info on the physical processes
+ */
+template< class NodalIndexSet, class Scalar, class PhysicsTraits >
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::GridView::dimension;
+ static constexpr int dimWorld = NodalIndexSet::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);
+ using GridView = typename NodalIndexSet::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 = Scalar;
//! export the type for the interaction volume index set
using IndexSet = CCMpfaOInteractionVolumeIndexSet< NodalIndexSet >;
//! export the type of interaction-volume local scvs
@@ -79,12 +94,10 @@ public:
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 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, PhysicsTraits >;
};
/*!
@@ -94,38 +107,38 @@ 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::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 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 +147,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 +163,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 +180,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 +208,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 +224,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 +244,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 +276,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,26 +291,24 @@ 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 constexpr 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
@@ -332,12 +344,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/localassembler.hh b/dumux/discretization/cellcentered/mpfa/omethod/localassembler.hh
index fd6700c7273d300f4433326e262320f653e77614..09f5eecde6e604cb663f3ac722ca05b481a0aade 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::ScalarType;
+ 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::ScalarType;
+ 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/properties.hh b/dumux/discretization/cellcentered/mpfa/properties.hh
index 06f3d8704a822edb5cd882f1d935ea00745607e5..220f4dbe53e6b4f04d9c475b85cce5b194650bab 100644
--- a/dumux/discretization/cellcentered/mpfa/properties.hh
+++ b/dumux/discretization/cellcentered/mpfa/properties.hh
@@ -81,10 +81,7 @@ public:
SET_PROP(CCMpfaModel, DualGridNodalIndexSet)
{
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;
+ using GV = typename GET_PROP_TYPE(TypeTag, GridView);
// per default, use uint8_t as iv-local index type
using LI = std::uint8_t;
@@ -95,32 +92,48 @@ private:
// 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 dim = GV::dimension;
static constexpr int maxE = dim == 3 ? 45 : 15;
public:
- using type = CCMpfaDualGridNodalIndexSet<GI, LI, dim, maxE, maxB>;
+ using type = CCMpfaDualGridNodalIndexSet<GV, LI, dim, maxE, maxB>;
};
//! The mpfa helper class
-SET_TYPE_PROP(CCMpfaModel, MpfaHelper, CCMpfaHelper<TypeTag>);
+SET_PROP(CCMpfaModel, MpfaHelper)
+{
+private:
+ using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
+public:
+ using type = CCMpfaHelper< FVGridGeometry >;
+};
//! 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 mpfa-o interaction volume as secondary type
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);
+
+ // set up the physics Traits
+ 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);
+ };
+ // use the default traits
+ using Traits = CCMpfaODefaultInteractionVolumeTraits< NodalIndexSet, Scalar, PhysicsTraits >;
public:
using type = CCMpfaOInteractionVolume< Traits >;
};
@@ -148,10 +161,10 @@ private:
PrimaryIV,
SecondaryIV >;
- template< class FVGridGeometry, bool enableCache >
- using LocalView = CCMpfaFVElementGeometry<FVGridGeometry, enableCache>;
+ template< class FVGridGeometry, bool enableGeomCache >
+ using LocalView = CCMpfaFVElementGeometry<FVGridGeometry, enableGeomCache>;
- //! Per default, we use the o-method and thus the simple assembly map
+ //! Use the correct connectivity map depending on mpfa scheme (obtain from primary iv)
template< class FVGridGeometry >
using ConnectivityMap = CCMpfaConnectivityMap<FVGridGeometry, FVGridGeometry::GridIVIndexSets::PrimaryInteractionVolume::MpfaMethod>;
};