Commit 0fac260a authored by Dennis Gläser's avatar Dennis Gläser Committed by Timo Koch
Browse files

[mpfa][IAVolume] define interaction volume interface

It is introduced a base interaction volume which defines the interface.
Additionally, a traits class is introduced. This defines the types of
vectors and matrices etc. to be used. The base class provides public typedefs
to be extracted from outside.
parent 7b58e77a
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*****************************************************************************
* See the file COPYING for full copying permissions. *
* *
* This program is free software: you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published by *
* the Free Software Foundation, either version 2 of the License, or *
* (at your option) any later version. *
* *
* This program is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU General Public License for more details. *
* *
* You should have received a copy of the GNU General Public License *
* along with this program. If not, see <http://www.gnu.org/licenses/>. *
*****************************************************************************/
/*!
* \file
* \brief Base class for interaction volumes of mpfa methods. Defines the interface.
*/
#ifndef DUMUX_DISCRETIZATION_CC_MPFA_INTERACTIONVOLUMEBASE_HH
#define DUMUX_DISCRETIZATION_CC_MPFA_INTERACTIONVOLUMEBASE_HH
#include <dumux/discretization/cellcentered/mpfa/interactionvolumeseed.hh>
#include <dumux/discretization/cellcentered/mpfa/methods.hh>
namespace Dumux
{
//! Implementation of the traits class for interaction volumes
template<class TypeTag, MpfaMethods method, int dim>
class CCMpfaInteractionVolumeTraitsImplementation
{};
//! traits class for interaction volumes
template<class TypeTag>
using CCMpfaInteractionVolumeTraits = CCMpfaInteractionVolumeTraitsImplementation<TypeTag,
GET_PROP_VALUE(TypeTag, MpfaMethod),
GET_PROP_TYPE(TypeTag, GridView)::dimension>;
/*!
* \ingroup Mpfa
* \brief Base class for the interaction volumes of mpfa methods.
* It defines the interface. Actual implementations should derive from this class.
*/
template<class TypeTag>
class CCMpfaInteractionVolumeBase
{
using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
using Traits = CCMpfaInteractionVolumeTraits<TypeTag>;
public:
// some types to be exported
using LocalIndexType = typename Traits::LocalIndexType;
using LocalIndexSet = typename Traits::LocalIndexSet;
using LocalIndexPair = typename Traits::LocalIndexPair;
using GlobalIndexType = typename Traits::GlobalIndexType;
using GlobalIndexSet = typename Traits::GlobalIndexSet;
using Vector = typename Traits::Vector;
using PositionVector = typename Traits::PositionVector;
using Seed = typename Traits::Seed;
//! solves the local equation system for the computation of the transmissibilities
template<typename GetTensorFunction>
void solveLocalSystem(const GetTensorFunction& getTensor)
{ DUNE_THROW(Dune::NotImplemented, "Actual interaction volume implementation does not provide a solveLocalSystem() method."); }
//! returns the dof indices in the stencil of the interaction volume
const GlobalIndexSet& stencil() const
{ DUNE_THROW(Dune::NotImplemented, "Actual interaction volume implementation does not provide a stencil() method."); }
//! returns the indices of the volvars in the stencil of the interaction volume
const GlobalIndexSet& volVarsStencil() const
{ DUNE_THROW(Dune::NotImplemented, "Actual interaction volume implementation does not provide a volVarsStencil() method."); }
//! returns the positions corresponding to the volvars in the stencil of the interaction volume (cell centers or scvf.ipGlobal() on boundary)
const PositionVector& volVarsPositions() const
{ DUNE_THROW(Dune::NotImplemented, "Actual interaction volume implementation does not provide a volVarsPositions() method."); }
//! returns a list of global scvf indices that are connected to this interaction volume
GlobalIndexSet globalScvfs() const
{ DUNE_THROW(Dune::NotImplemented, "Actual interaction volume implementation does not provide a globalScvfs() method."); }
//! returns the local index of an scvf in the IV and a boolean whether or not it is on the negative side of the local scvf (flux has to be inverted)
LocalIndexPair getLocalIndexPair(const SubControlVolumeFace& scvf) const
{ DUNE_THROW(Dune::NotImplemented, "Actual interaction volume implementation does not provide a getLocalIndexPair() method."); }
//! returns the transmissibilities corresponding to a local scvf
Vector getTransmissibilities(const LocalIndexPair& localIndexPair) const
{ DUNE_THROW(Dune::NotImplemented, "Actual interaction volume implementation does not provide a getTransmissibilities() method."); }
//! returns the neumann flux corresponding to a local scvf
Scalar getNeumannFlux(LocalIndexPair& localIndexPair) const
{ DUNE_THROW(Dune::NotImplemented, "Actual interaction volume implementation does not provide a getNeumannFlux() method."); }
//! returns the local index in a vector for a given global index
template<typename IdxType1, typename IdxType2>
LocalIndexType findLocalIndex(const std::vector<IdxType1>& vector, const IdxType2 globalIdx) const
{
auto it = std::find(vector.begin(), vector.end(), globalIdx);
assert(it != vector.end() && "could not find local index in the vector for the given global index!");
return std::distance(vector.begin(), it);
}
};
} // end namespace
#endif
......@@ -25,6 +25,7 @@
#include <dumux/common/math.hh>
#include <dumux/implicit/cellcentered/mpfa/properties.hh>
#include <dumux/discretization/cellcentered/mpfa/interactionvolumebase.hh>
#include <dumux/discretization/cellcentered/mpfa/interactionvolumeseed.hh>
#include <dumux/discretization/cellcentered/mpfa/facetypes.hh>
#include <dumux/discretization/cellcentered/mpfa/methods.hh>
......@@ -34,12 +35,40 @@
namespace Dumux
{
//! Specialization of the interaction volume traits clas
template<class TypeTag, int dim>
class CCMpfaInteractionVolumeTraitsImplementation<TypeTag, MpfaMethods::oMethod, dim>
{
using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
static const int dimWorld = GridView::dimensionworld;
public:
using GlobalPosition = Dune::FieldVector<Scalar, dimWorld>;
using DimVector = Dune::FieldVector<Scalar, dim>;
using Tensor = Dune::FieldMatrix<Scalar, dim, dim>;
using LocalIndexType = std::uint8_t;
using LocalIndexSet = std::vector<LocalIndexType>;
using LocalIndexPair = std::pair<LocalIndexType, bool>;
using GlobalIndexType = typename GridView::IndexSet::IndexType;
using GlobalIndexSet = std::vector<GlobalIndexType>;
using Matrix = Dune::DynamicMatrix<Scalar>;
using Vector = typename Matrix::row_type;
using PositionVector = std::vector<GlobalPosition>;
using LocalScvSeedType = CCMpfaOLocalScvSeed<GlobalIndexType, LocalIndexType>;
using LocalScvfSeedType = CCMpfaOLocalScvfSeed<GlobalIndexType, LocalIndexType>;
using Seed = CCMpfaInteractionVolumeSeed<LocalScvSeedType, LocalScvfSeedType>;
};
/*!
* \ingroup Mpfa
* \brief Base class for the interaction volumes of the mpfa-o method
*/
template<class TypeTag>
class CCMpfaInteractionVolumeImplementation<TypeTag, MpfaMethods::oMethod>
class CCMpfaInteractionVolumeImplementation<TypeTag, MpfaMethods::oMethod> : public CCMpfaInteractionVolumeBase<TypeTag>
{
using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
......@@ -49,40 +78,37 @@ class CCMpfaInteractionVolumeImplementation<TypeTag, MpfaMethods::oMethod>
using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);
static const int dim = GridView::dimension;
using Element = typename GridView::template Codim<0>::Entity;
using GlobalIndexType = typename GridView::IndexSet::IndexType;
using LocalIdxType = std::uint8_t;
using LocalScvSeedType = CCMpfaOLocalScvSeed<GlobalIndexType, LocalIdxType>;
using LocalScvfSeedType = CCMpfaOLocalScvfSeed<GlobalIndexType, LocalIdxType>;
// choose names different from those in the base class
// base class typedefs are public to be used by other classes
using Traits = CCMpfaInteractionVolumeTraits<TypeTag>;
using GlobalIdxType = typename Traits::GlobalIndexType;
using LocalIdxType = typename Traits::LocalIndexType;
using GlobalIdxSet = typename Traits::GlobalIndexSet;
using LocalIdxSet = typename Traits::LocalIndexSet;
using LocalIdxPair = typename Traits::LocalIndexPair;
using IVSeed = typename Traits::Seed;
using GlobalPosition = typename Traits::GlobalPosition;
using Vector = typename Traits::Vector;
using Matrix = typename Traits::Matrix;
using Tensor = typename Traits::Tensor;
using LocalScvType = CCMpfaOLocalScv<TypeTag>;
using LocalScvfType = CCMpfaOLocalScvf<TypeTag>;
static const int dim = GridView::dimension;
using DimVector = Dune::FieldVector<Scalar, dim>;
using Tensor = Dune::FieldMatrix<Scalar, dim, dim>;
using Stencil = std::vector<GlobalIndexType>;
using LocalIndexSet = std::vector<LocalIdxType>;
using Matrix = Dune::DynamicMatrix<Scalar>;
using Vector = typename Matrix::row_type;
public:
// the interaction volume seed type to be exported
using LocalIndexType = LocalIdxType;
using Seed = CCMpfaInteractionVolumeSeed<LocalScvSeedType, LocalScvfSeedType>;
CCMpfaInteractionVolumeImplementation(const Seed& seed,
CCMpfaInteractionVolumeImplementation(const IVSeed& seed,
const Problem& problem,
const FVElementGeometry& fvGeometry,
const ElementVolumeVariables& elemVolVars)
: problemRef_(problem),
fvGeometryRef_(fvGeometry),
elemVolVarsRef_(elemVolVars),
onBoundary_(seed.onBoundary())
onBoundary_(seed.onBoundary()),
hasInteriorBoundary_(false)
{
// create local sub control entities from the seed
createLocalEntities_(seed);
......@@ -100,25 +126,28 @@ public:
volVarsPositions_.push_back(localScv.center());
}
// eventually add dirichlet vol var indices
// eventually add dirichlet vol var indices and set up local index sets of flux and dirichlet faces
LocalIdxType localScvfIdx = 0;
for (const auto& localScvf : localScvfs_)
{
if (localScvf.faceType() == MpfaFaceTypes::dirichlet || localScvf.faceType() == MpfaFaceTypes::interiorDirichlet)
auto faceType = localScvf.faceType();
// eventually add vol var index and corresponding position
if (faceType == MpfaFaceTypes::dirichlet || faceType == MpfaFaceTypes::interiorDirichlet)
{
volVarsStencil_.push_back(localScvf.outsideGlobalScvIndex());
volVarsPositions_.push_back(localScvf.ip());
}
}
// set up local index sets of flux and dirichlet faces
LocalIndexType localScvfIdx = 0;
for (const auto& localScvf : localScvfs_)
{
auto faceType = localScvf.faceType();
// fill local scvf type index sets accordingly
if (faceType != MpfaFaceTypes::dirichlet && faceType != MpfaFaceTypes::interiorDirichlet)
fluxFaceIndexSet_.push_back(localScvfIdx++);
else
dirichletFaceIndexSet_.push_back(localScvfIdx++);
// if face is on an interior boundary, save this information
if (!hasInteriorBoundary_ && (faceType == MpfaFaceTypes::interiorDirichlet || faceType == MpfaFaceTypes::interiorNeumann))
hasInteriorBoundary_ = true;
}
neumannFluxes_ = Vector(fluxFaceIndexSet_.size(), 0.0);
......@@ -158,7 +187,7 @@ public:
const unsigned int eqIdx)
{
LocalIndexType fluxFaceIdx = 0;
LocalIdxType fluxFaceIdx = 0;
for (const auto localFluxFaceIdx : fluxFaceIndexSet_)
{
const auto& localScvf = localScvf_(localFluxFaceIdx);
......@@ -192,35 +221,46 @@ public:
bool onBoundary() const
{ return onBoundary_; }
const Stencil& stencil() const
bool hasInteriorBoundary() const
{ return hasInteriorBoundary_; }
const GlobalIdxSet& stencil() const
{ return stencil_; }
const Stencil& volVarsStencil() const
const GlobalIdxSet& volVarsStencil() const
{ return volVarsStencil_; }
const std::vector<DimVector>& volVarsPositions() const
const std::vector<GlobalPosition>& volVarsPositions() const
{ return volVarsPositions_; }
Stencil globalScvfs() const
GlobalIdxSet globalScvfs() const
{
Stencil stencil;
stencil.reserve(localScvfs_.size());
GlobalIdxSet globalScvfs;
globalScvfs.reserve(localScvfs_.size());
for (const auto& localScvf : localScvfs_)
{
stencil.push_back(localScvf.insideGlobalScvfIndex());
globalScvfs.push_back(localScvf.insideGlobalScvfIndex());
if (!localScvf.boundary())
stencil.push_back(localScvf.outsideGlobalScvfIndex());
globalScvfs.push_back(localScvf.outsideGlobalScvfIndex());
}
// if interior boundaries are present we have to make the entries unique
if (hasInteriorBoundary())
{
// make values in elementstencil unique
std::sort(globalScvfs.begin(), globalScvfs.end());
globalScvfs.erase(std::unique(globalScvfs.begin(), globalScvfs.end()), globalScvfs.end());
}
return stencil;
return globalScvfs;
}
std::pair<LocalIndexType, bool> getLocalIndexPair(const SubControlVolumeFace& scvf) const
LocalIdxPair getLocalIndexPair(const SubControlVolumeFace& scvf) const
{
auto scvfGlobalIdx = scvf.index();
LocalIndexType localIdx = 0;
LocalIdxType localIdx = 0;
for (const auto& localScvf : localScvfs_)
{
if (localScvf.insideGlobalScvfIndex() == scvfGlobalIdx)
......@@ -235,7 +275,7 @@ public:
DUNE_THROW(Dune::InvalidStateException, "Could not find the local scv face in the interaction volume for the given scvf with index: " << scvf.index());
}
Vector getTransmissibilities(const std::pair<LocalIndexType, bool>& localIndexPair) const
Vector getTransmissibilities(const std::pair<LocalIdxType, bool>& localIndexPair) const
{
auto tij = T_[localIndexPair.first];
......@@ -244,7 +284,7 @@ public:
return tij;
}
Scalar getNeumannFlux(const std::pair<LocalIndexType, bool>& localIndexPair) const
Scalar getNeumannFlux(const std::pair<LocalIdxType, bool>& localIndexPair) const
{
if (fluxScvfIndexSet_().size() == 0)
return 0.0;
......@@ -257,22 +297,22 @@ public:
private:
const LocalScvfType& localScvf_(const LocalIndexType localScvfIdx) const
const LocalScvfType& localScvf_(const LocalIdxType localScvfIdx) const
{ return localScvfs_[localScvfIdx]; }
const LocalScvType& localScv_(const LocalIndexType localScvIdx) const
const LocalScvType& localScv_(const LocalIdxType localScvIdx) const
{ return localScvs_[localScvIdx]; }
const LocalIndexSet& fluxScvfIndexSet_() const
const LocalIdxSet& fluxScvfIndexSet_() const
{ return fluxFaceIndexSet_; }
const LocalIndexSet& dirichletScvfIndexSet_() const
const LocalIdxSet& dirichletScvfIndexSet_() const
{ return dirichletFaceIndexSet_; }
const Element& localElement_(const LocalIndexType localScvIdx) const
const Element& localElement_(const LocalIdxType localScvIdx) const
{ return localElements_[localScvIdx]; }
void createLocalEntities_(const Seed& seed)
void createLocalEntities_(const IVSeed& seed)
{
auto numScvs = seed.scvSeeds().size();
auto numScvfs = seed.scvfSeeds().size();
......@@ -304,12 +344,12 @@ private:
const std::size_t numLocalScvs = localScvs_.size();
// loop over the local faces
LocalIndexType rowIdx = 0;
LocalIdxType rowIdx = 0;
for (const auto& localScvf : localScvfs_)
{
auto faceType = localScvf.faceType();
bool hasUnknown = faceType != MpfaFaceTypes::dirichlet && faceType != MpfaFaceTypes::interiorDirichlet;
LocalIndexType idxInFluxFaces = hasUnknown ? findLocalIndex_(fluxScvfIndexSet_(), rowIdx) : -1;
LocalIdxType idxInFluxFaces = hasUnknown ? this->findLocalIndex(fluxScvfIndexSet_(), rowIdx) : -1;
// get diffusion tensor in "positive" sub volume
const auto posLocalScvIdx = localScvf.insideLocalScvIndex();
......@@ -334,7 +374,7 @@ private:
if (curFaceHasUnknown)
{
// we need the index of the current local scvf in the flux face indices
auto curIdxInFluxFaces = findLocalIndex_(fluxScvfIndexSet_(), curLocalScvfIdx);
auto curIdxInFluxFaces = this->findLocalIndex(fluxScvfIndexSet_(), curLocalScvfIdx);
C[rowIdx][curIdxInFluxFaces] += posWijk[localDir];
if (hasUnknown)
......@@ -353,7 +393,7 @@ private:
else
{
// the current face is a Dirichlet face and creates entries in D & eventually B
auto curIdxInDiriFaces = findLocalIndex_(dirichletScvfIndexSet_(), curLocalScvfIdx);
auto curIdxInDiriFaces = this->findLocalIndex(dirichletScvfIndexSet_(), curLocalScvfIdx);
D[rowIdx][numLocalScvs + curIdxInDiriFaces] += posWijk[localDir];
if (hasUnknown)
......@@ -396,7 +436,7 @@ private:
if (curFaceHasUnknown)
{
// we need the index of the current local scvf in the flux face indices
auto curIdxInFluxFaces = findLocalIndex_(fluxScvfIndexSet_(), curLocalScvfIdx);
auto curIdxInFluxFaces = this->findLocalIndex(fluxScvfIndexSet_(), curLocalScvfIdx);
if (faceType == MpfaFaceTypes::interiorNeumann)
A[idxInFluxFaces][curIdxInFluxFaces] -= (1.0-xi)*negWijk[localDir];
......@@ -406,7 +446,7 @@ private:
else
{
// the current face is a Dirichlet face and creates entries in B
auto curIdxInDiriFaces = findLocalIndex_(dirichletScvfIndexSet_(), curLocalScvfIdx);
auto curIdxInDiriFaces = this->findLocalIndex(dirichletScvfIndexSet_(), curLocalScvfIdx);
B[idxInFluxFaces][numLocalScvs + curIdxInDiriFaces] += negWijk[localDir];
}
......@@ -435,7 +475,7 @@ private:
CAinv_.resize(0, 0);
// Loop over all the faces, in this case these are all dirichlet boundaries
LocalIndexType rowIdx = 0;
LocalIdxType rowIdx = 0;
for (const auto& localScvf : localScvfs_)
{
// get diffusion tensor in "positive" sub volume
......@@ -452,7 +492,7 @@ private:
for (int localDir = 0; localDir < dim; localDir++)
{
auto curLocalScvfIdx = posLocalScv.localScvfIndex(localDir);
auto curIdxInDiriFaces = findLocalIndex_(dirichletScvfIndexSet_(), curLocalScvfIdx);
auto curIdxInDiriFaces = this->findLocalIndex(dirichletScvfIndexSet_(), curLocalScvfIdx);
T_[rowIdx][numLocalScvs + curIdxInDiriFaces] += posWijk[localDir];
T_[rowIdx][posLocalScvIdx] -= posWijk[localDir];
......@@ -463,12 +503,12 @@ private:
}
}
DimVector calculateOmegas_(const LocalScvType& localScv,
GlobalPosition calculateOmegas_(const LocalScvType& localScv,
const LocalScvfType& localScvf,
const Tensor& T) const
{
DimVector wijk;
DimVector tmp;
GlobalPosition wijk;
GlobalPosition tmp;
for (int dir = 0; dir < dim; ++dir)
{
T.mv(localScv.innerNormal(dir), tmp);
......@@ -481,12 +521,12 @@ private:
return wijk;
}
DimVector calculateOmegas_(const LocalScvType& localScv,
GlobalPosition calculateOmegas_(const LocalScvType& localScv,
const LocalScvfType& localScvf,
const Scalar t) const
{
DimVector wijk;
DimVector tmp(localScvf.unitOuterNormal());
GlobalPosition wijk;
GlobalPosition tmp(localScvf.unitOuterNormal());
tmp *= t;
for (int dir = 0; dir < dim; ++dir)
......@@ -498,14 +538,6 @@ private:
return wijk;
}
template<typename IdxType1, typename IdxType2>
LocalIndexType findLocalIndex_(const std::vector<IdxType1>& vector, const IdxType2 globalIdx) const
{
auto it = std::find(vector.begin(), vector.end(), globalIdx);
assert(it != vector.end() && "could not find local index in the vector for the given global index!");
return std::distance(vector.begin(), it);
}
const Problem& problem_() const
{ return problemRef_; }
......@@ -520,19 +552,20 @@ private:
const ElementVolumeVariables& elemVolVarsRef_;
bool onBoundary_;
LocalIndexType eqIdx_;
bool hasInteriorBoundary_;
LocalIdxType eqIdx_;
std::vector<Element> localElements_;
std::vector<LocalScvType> localScvs_;
std::vector<LocalScvfType> localScvfs_;
Stencil globalScvfIndices_;
Stencil stencil_;
Stencil volVarsStencil_;
std::vector<DimVector> volVarsPositions_;
GlobalIdxSet globalScvfIndices_;
GlobalIdxSet stencil_;
GlobalIdxSet volVarsStencil_;
std::vector<GlobalPosition> volVarsPositions_;
LocalIndexSet fluxFaceIndexSet_;
LocalIndexSet dirichletFaceIndexSet_;
LocalIdxSet fluxFaceIndexSet_;
LocalIdxSet dirichletFaceIndexSet_;
Matrix T_;
Matrix AinvB_;
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment