From 8feea16b6d666095902e1c8da15f49cfdda3fb09 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Dennis=20Gl=C3=A4ser?= <dennis.glaeser@iws.uni-stuttgart.de>
Date: Wed, 21 Nov 2018 11:56:45 +0100
Subject: [PATCH] [mpfa][iv] put omegas in data handle
---
.../mpfa/interactionvolumedatahandle.hh | 6 ++++
.../cellcentered/mpfa/localassemblerbase.hh | 2 +-
.../mpfa/omethod/interactionvolume.hh | 31 +++++--------------
.../mpfa/omethod/localassembler.hh | 17 ++++++----
.../mpfa/omethod/staticinteractionvolume.hh | 25 ++++-----------
5 files changed, 32 insertions(+), 49 deletions(-)
diff --git a/dumux/discretization/cellcentered/mpfa/interactionvolumedatahandle.hh b/dumux/discretization/cellcentered/mpfa/interactionvolumedatahandle.hh
index baee94234a..ebcebfd5a5 100644
--- a/dumux/discretization/cellcentered/mpfa/interactionvolumedatahandle.hh
+++ b/dumux/discretization/cellcentered/mpfa/interactionvolumedatahandle.hh
@@ -64,6 +64,7 @@ class MatrixDataHandleBase
using TMatrix = typename MVT::TMatrix;
using CellVector = typename MVT::CellVector;
using OutsideTij = std::vector< std::vector<CellVector> >;
+ using OmegaStorage = typename MVT::OmegaStorage;
public:
//! Access functions to context-dependent data
@@ -82,6 +83,9 @@ public:
const OutsideTij& tijOutside() const { return tijOutside_[contextIdx1_][contextIdx2_]; }
OutsideTij& tijOutside() { return tijOutside_[contextIdx1_][contextIdx2_]; }
+ const OmegaStorage& omegas() const { return wijk_[contextIdx1_][contextIdx2_]; }
+ OmegaStorage& omegas() { return wijk_[contextIdx1_][contextIdx2_]; }
+
//! functionality to set the context indices
void setContextIndex1(unsigned int idx) const { assert(idx < size1); contextIdx1_ = idx; }
void setContextIndex2(unsigned int idx) const { assert(idx < size2); contextIdx2_ = idx; }
@@ -91,6 +95,8 @@ protected:
mutable unsigned int contextIdx1_{0};
mutable unsigned int contextIdx2_{0};
+ std::array< std::array<OmegaStorage, size2>, size1 > wijk_; //!< The omega factors that form the entries of the matrices
+
std::array< std::array<TMatrix, size2>, size1 > T_; //!< The transmissibility matrix
std::array< std::array<AMatrix, size2>, size1 > A_; //!< Inverse of the iv-local system matrix
std::array< std::array<CMatrix, size2>, size1 > AB_; //!< A_ left multiplied to B
diff --git a/dumux/discretization/cellcentered/mpfa/localassemblerbase.hh b/dumux/discretization/cellcentered/mpfa/localassemblerbase.hh
index 20cbb2a1cc..56c350025d 100644
--- a/dumux/discretization/cellcentered/mpfa/localassemblerbase.hh
+++ b/dumux/discretization/cellcentered/mpfa/localassemblerbase.hh
@@ -266,7 +266,7 @@ class InteractionVolumeAssemblerBase
const auto localScvfIdx = localFaceData.ivLocalScvfIndex();
const auto idxInOutside = localFaceData.scvfLocalOutsideScvfIndex();
const auto& posLocalScv = iv.localScv(localScvIdx);
- const auto& wijk = iv.omegas()[localScvfIdx][idxInOutside+1];
+ const auto& wijk = handle.omegas()[localScvfIdx][idxInOutside+1];
// add contributions from all local directions
for (LocalIndexType localDir = 0; localDir < dim; localDir++)
diff --git a/dumux/discretization/cellcentered/mpfa/omethod/interactionvolume.hh b/dumux/discretization/cellcentered/mpfa/omethod/interactionvolume.hh
index 8b557cb505..ed91b8d447 100644
--- a/dumux/discretization/cellcentered/mpfa/omethod/interactionvolume.hh
+++ b/dumux/discretization/cellcentered/mpfa/omethod/interactionvolume.hh
@@ -67,9 +67,16 @@ private:
static constexpr int dim = NodalIndexSet::Traits::GridView::dimension;
static constexpr int dimWorld = NodalIndexSet::Traits::GridView::dimensionworld;
+ using DimVector = Dune::FieldVector<Scalar, dim>;
+ using FaceOmegas = typename std::conditional< (dim<dimWorld),
+ std::vector<DimVector>,
+ Dune::ReservedVector<DimVector, 2> >::type;
+
//! Matrix/Vector traits to be used by the data handle
struct MVTraits
{
+ using OmegaStorage = std::vector< FaceOmegas >;
+
using AMatrix = Dune::DynamicMatrix< Scalar >;
using BMatrix = Dune::DynamicMatrix< Scalar >;
using CMatrix = Dune::DynamicMatrix< Scalar >;
@@ -116,16 +123,6 @@ class CCMpfaOInteractionVolume
using LocalIndexType = typename IndexSet::LocalIndexType;
using Stencil = typename IndexSet::NodalGridStencilType;
- static constexpr int dim = GridView::dimension;
- static constexpr int dimWorld = GridView::dimensionworld;
-
- //! export scalar type from T matrix and define omegas
- using Scalar = typename Traits::MatVecTraits::TMatrix::value_type;
- using DimVector = Dune::FieldVector<Scalar, dim>;
- using FaceOmegas = typename std::conditional< (dim<dimWorld),
- std::vector<DimVector>,
- Dune::ReservedVector<DimVector, 2> >::type;
-
using LocalScvType = typename Traits::LocalScvType;
using LocalScvfType = typename Traits::LocalScvfType;
using LocalFaceData = typename Traits::LocalFaceData;
@@ -144,9 +141,6 @@ public:
GridIndexType volVarIndex() const { return volVarIndex_; }
};
- //! 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;
@@ -193,7 +187,6 @@ public:
numKnowns_ = numLocalScvs;
// set up quantitites related to sub-control volume faces
- wijk_.resize(numFaces_);
for (LocalIndexType faceIdxLocal = 0; faceIdxLocal < numFaces_; ++faceIdxLocal)
{
const auto& scvf = fvGeometry.scvf(indexSet.gridScvfIndex(faceIdxLocal));
@@ -206,7 +199,6 @@ public:
// 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())
@@ -230,7 +222,7 @@ public:
{
// loop over scvfs in outside scv until we find the one coinciding with current scvf
const auto outsideLocalScvIdx = neighborScvIndicesLocal[i];
- for (int coord = 0; coord < dim; ++coord)
+ for (int coord = 0; coord < GridView::dimension; ++coord)
{
if (indexSet.localScvfIndex(outsideLocalScvIdx, coord) == faceIdxLocal)
{
@@ -291,10 +283,6 @@ public:
const std::vector<DirichletData>& dirichletData() const
{ return dirichletData_; }
- //! returns container storing the transmissibilities for each face & coordinate
- const TransmissibilityStorage& omegas() const { return wijk_; }
- TransmissibilityStorage& omegas() { return wijk_; }
-
//! returns the number of interaction volumes living around a vertex
template< class NI >
static constexpr std::size_t numIVAtVertex(const NI& nodalIndexSet) { return 1; }
@@ -332,9 +320,6 @@ private:
std::vector<LocalFaceData> localFaceData_;
std::vector<DirichletData> dirichletData_;
- // The omega factors are stored during assembly of local system
- TransmissibilityStorage wijk_;
-
// sizes involved in the local system equations
std::size_t numFaces_;
std::size_t numUnknowns_;
diff --git a/dumux/discretization/cellcentered/mpfa/omethod/localassembler.hh b/dumux/discretization/cellcentered/mpfa/omethod/localassembler.hh
index dafa103d8c..a403964578 100644
--- a/dumux/discretization/cellcentered/mpfa/omethod/localassembler.hh
+++ b/dumux/discretization/cellcentered/mpfa/omethod/localassembler.hh
@@ -70,7 +70,7 @@ public:
template< class DataHandle, class IV, class TensorFunc >
void assembleMatrices(DataHandle& handle, IV& iv, const TensorFunc& getT)
{
- assembleLocalMatrices_(handle.A(), handle.AB(), handle.CA(), handle.T(), iv, getT);
+ assembleLocalMatrices_(handle.A(), handle.AB(), handle.CA(), handle.T(), handle.omegas(), iv, getT);
// maybe solve the local system
if (iv.numUnknowns() > 0)
@@ -113,7 +113,7 @@ public:
const auto scvfIdx = localFaceData.ivLocalScvfIndex();
const auto idxInOut = localFaceData.scvfLocalOutsideScvfIndex();
- const auto& wijk = iv.omegas()[scvfIdx][idxInOut+1];
+ const auto& wijk = handle.omegas()[scvfIdx][idxInOut+1];
auto& tij = tijOut[scvfIdx][idxInOut];
tij = 0.0;
@@ -194,12 +194,16 @@ private:
typename IV::Traits::MatVecTraits::BMatrix& B,
typename IV::Traits::MatVecTraits::CMatrix& C,
typename IV::Traits::MatVecTraits::DMatrix& D,
+ typename IV::Traits::MatVecTraits::OmegaStorage& wijk,
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;
+ // resize omegas
+ this->resizeVector_(wijk, iv.numFaces());
+
// if only Dirichlet faces are present in the iv,
// the matrices A, B & C are undefined and D = T
if (iv.numUnknowns() == 0)
@@ -222,15 +226,16 @@ private:
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());
+ this->resizeVector_(wijk[faceIdx], /*no outside scvs present*/1);
+ wijk[faceIdx][0] = computeMpfaTransmissibility(posLocalScv, curGlobalScvf, tensor, posVolVars.extrusionFactor());
const auto posScvLocalDofIdx = posLocalScv.localDofIndex();
for (LocalIndexType localDir = 0; localDir < dim; localDir++)
{
const auto& otherLocalScvf = iv.localScvf( posLocalScv.localScvfIndex(localDir) );
const auto otherLocalDofIdx = otherLocalScvf.localDofIndex();
- D[faceIdx][otherLocalDofIdx] -= wijk[localDir];
- D[faceIdx][posScvLocalDofIdx] += wijk[localDir];
+ D[faceIdx][otherLocalDofIdx] -= wijk[faceIdx][0][localDir];
+ D[faceIdx][posScvLocalDofIdx] += wijk[faceIdx][0][localDir];
}
}
}
@@ -242,7 +247,6 @@ private:
this->resizeMatrix_(C, iv.numFaces(), iv.numUnknowns()); C = 0.0;
this->resizeMatrix_(D, iv.numFaces(), iv.numKnowns()); D = 0.0;
- auto& wijk = iv.omegas();
for (LocalIndexType faceIdx = 0; faceIdx < iv.numFaces(); ++faceIdx)
{
const auto& curLocalScvf = iv.localScvf(faceIdx);
@@ -259,6 +263,7 @@ private:
const auto tensor = getT(this->problem(), posElement, posVolVars, this->fvGeometry(), posGlobalScv);
// the omega factors of the "positive" sub volume
+ this->resizeVector_(wijk[faceIdx], neighborScvIndices.size());
wijk[faceIdx][0] = computeMpfaTransmissibility(posLocalScv, curGlobalScvf, tensor, posVolVars.extrusionFactor());
// go over the coordinate directions in the positive sub volume
diff --git a/dumux/discretization/cellcentered/mpfa/omethod/staticinteractionvolume.hh b/dumux/discretization/cellcentered/mpfa/omethod/staticinteractionvolume.hh
index 82cbe5c1a7..25dc7c492f 100644
--- a/dumux/discretization/cellcentered/mpfa/omethod/staticinteractionvolume.hh
+++ b/dumux/discretization/cellcentered/mpfa/omethod/staticinteractionvolume.hh
@@ -68,9 +68,14 @@ private:
static constexpr int dim = NI::Traits::GridView::dimension;
static constexpr int dimWorld = NI::Traits::GridView::dimensionworld;
+ using DimVector = Dune::FieldVector<S, dim>;
+ using FaceOmegas = Dune::ReservedVector<DimVector, 2>;
+
//! Matrix/Vector traits to be used by the data handle
struct MVTraits
{
+ using OmegaStorage = std::array< FaceOmegas, F >;
+
using AMatrix = Dune::FieldMatrix< S, F, F >;
using BMatrix = Dune::FieldMatrix< S, F, C >;
using CMatrix = Dune::FieldMatrix< S, F, F >;
@@ -125,14 +130,6 @@ class CCMpfaOStaticInteractionVolume
using LocalIndexType = typename IndexSet::LocalIndexType;
using Stencil = typename IndexSet::NodalGridStencilType;
- static constexpr int dim = GridView::dimension;
- static constexpr int dimWorld = GridView::dimensionworld;
-
- //! export scalar type from T matrix and define omegas
- using Scalar = typename Traits::MatVecTraits::TMatrix::value_type;
- using DimVector = Dune::FieldVector<Scalar, dim>;
- using FaceOmegas = typename Dune::ReservedVector<DimVector, 2>;
-
using LocalScvType = typename Traits::LocalScvType;
using LocalScvfType = typename Traits::LocalScvfType;
using LocalFaceData = typename Traits::LocalFaceData;
@@ -144,9 +141,6 @@ public:
//! export the standard o-methods dirichlet data
using DirichletData = typename CCMpfaOInteractionVolume< Traits >::DirichletData;
- //! export the type used for transmissibility storage
- using TransmissibilityStorage = std::array< FaceOmegas, numScvf >;
-
//! publicly state the mpfa-scheme this interaction volume is associated with
static constexpr MpfaMethods MpfaMethod = MpfaMethods::oMethod;
@@ -190,7 +184,7 @@ public:
// add local face data objects for the outside face
const auto outsideLocalScvIdx = neighborScvIndicesLocal[1];
- for (int coord = 0; coord < dim; ++coord)
+ for (int coord = 0; coord < GridView::dimension; ++coord)
{
if (indexSet.localScvfIndex(outsideLocalScvIdx, coord) == faceIdxLocal)
{
@@ -250,10 +244,6 @@ public:
const std::array<DirichletData, 0>& dirichletData() const
{ return dirichletData_; }
- //! returns container storing the transmissibilities for each face & coordinate
- const TransmissibilityStorage& omegas() const { return wijk_; }
- TransmissibilityStorage& omegas() { return wijk_; }
-
//! returns the number of interaction volumes living around a vertex
template< class NI >
static constexpr std::size_t numIVAtVertex(const NI& nodalIndexSet)
@@ -291,9 +281,6 @@ private:
std::array<LocalScvfType, numScvf> scvfs_;
std::array<LocalFaceData, numScvf*2> localFaceData_;
- // The omega factors are stored during assembly of local system
- TransmissibilityStorage wijk_;
-
// Dummy dirichlet data container (compatibility with dynamic o-iv)
std::array<DirichletData, 0> dirichletData_;
};
--
GitLab