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

[mpfa] allow for different interaction volume seed types

Up to now the types for the inner and boundary interaction volume seeds
were always identical. For the mpfa-l and other methods, it would be more efficient
if these types could differ, as e.g. the sizes of the containers are known at compile time
for the l method. This commit introduces the capability of the types to differ.
parent d6ce2b51
......@@ -136,11 +136,19 @@ public:
for (auto&& scvf : scvfs(fvGeometry))
{
const bool boundary = globalFvGeometry.scvfTouchesBoundary(scvf);
const auto& ivSeed = boundary ? globalFvGeometry.boundaryInteractionVolumeSeed(scvf) : globalFvGeometry.interactionVolumeSeed(scvf);
// loop over all the scvfs in the interaction region
for (auto scvfIdx : ivSeed.globalScvfIndices())
globalScvfIndices_.push_back(scvfIdx);
if (boundary)
{
const auto& ivSeed = globalFvGeometry.boundaryInteractionVolumeSeed(scvf);
for (auto scvfIdx : ivSeed.globalScvfIndices())
globalScvfIndices_.push_back(scvfIdx);
}
else
{
const auto& ivSeed = globalFvGeometry.interactionVolumeSeed(scvf);
for (auto scvfIdx : ivSeed.globalScvfIndices())
globalScvfIndices_.push_back(scvfIdx);
}
}
// make global indices unique
std::sort(globalScvfIndices_.begin(), globalScvfIndices_.end());
......
......@@ -181,13 +181,13 @@ public:
volumeVariables_.reserve(numDofs+neighborScvfEstimate);
volVarIndices_.reserve(numDofs+neighborScvfEstimate);
// skip the rest if the scvf does not touch a boundary
const bool boundary = globalFvGeometry.scvfTouchesBoundary(scvf);
const auto& ivSeed = boundary ? globalFvGeometry.boundaryInteractionVolumeSeed(scvf) : globalFvGeometry.interactionVolumeSeed(scvf);
if (!ivSeed.onBoundary())
if (!boundary)
continue;
// loop over all the scvfs in the interaction region
const auto& ivSeed = globalFvGeometry.boundaryInteractionVolumeSeed(scvf);
for (auto scvfIdx : ivSeed.globalScvfIndices())
{
auto&& ivScvf = fvGeometry.scvf(scvfIdx);
......
......@@ -72,8 +72,9 @@ public:
FluxVarsCacheVector& fluxVarsCache,
const bool updateNeumannOnly = false)
{
// if we only want to update the neumann fluxes, skip the rest if we don't touch the boundary
const bool boundary = problem.model().globalFvGeometry().scvfTouchesBoundary(scvf);
// if we only want to update the neumann fluxes, skip the rest if we don't touch the boundary
if (updateNeumannOnly && !boundary)
return;
......@@ -84,10 +85,10 @@ public:
const SubControlVolume& scv)
{ return prob->spatialParams().intrinsicPermeability(scv, volVars); };
// get the interaction volume seed and update flux var caches
const auto& seed = boundary ? problem.model().globalFvGeometry().boundaryInteractionVolumeSeed(scvf) : problem.model().globalFvGeometry().interactionVolumeSeed(scvf);
// update flux var caches
if (boundary)
{
const auto& seed = problem.model().globalFvGeometry().boundaryInteractionVolumeSeed(scvf);
BoundaryInteractionVolume iv(seed, problem, fvGeometry, elemVolVars);
iv.solveLocalSystem(permFunction);
......@@ -131,6 +132,7 @@ public:
}
else
{
const auto& seed = problem.model().globalFvGeometry().interactionVolumeSeed(scvf);
InteractionVolume iv(seed, problem, fvGeometry, elemVolVars);
iv.solveLocalSystem(permFunction);
......
......@@ -240,22 +240,19 @@ public:
const auto eIdx = globalFvGeometry().problem_().elementMapper().index(element);
for (auto&& scvf : scvfs_)
{
// skip the rest if scvf has been handled already
if (contains(scvf.vertexIndex(), finishedVertices))
continue;
const bool boundary = globalFvGeometry().scvfTouchesBoundary(scvf);
const auto& ivSeed = boundary ? globalFvGeometry().boundaryInteractionVolumeSeed(scvf) : globalFvGeometry().interactionVolumeSeed(scvf);
auto scvfIndices = ivSeed.globalScvfIndices();
// make the scv for each element in the interaction region
// and the corresponding scv faces in that scv
for (auto eIdxGlobal : ivSeed.globalScvIndices())
if (globalFvGeometry().scvfTouchesBoundary(scvf))
{
if (eIdxGlobal != eIdx)
{
auto element = globalFvGeometry().element(eIdxGlobal);
makeNeighborGeometries(element, eIdxGlobal, scvfIndices);
}
const auto& ivSeed = globalFvGeometry().boundaryInteractionVolumeSeed(scvf);
handleInteractionRegion(ivSeed);
}
else
{
const auto& ivSeed = globalFvGeometry().interactionVolumeSeed(scvf);
handleInteractionRegion(ivSeed);
}
finishedVertices.push_back(scvf.vertexIndex());
......@@ -288,10 +285,12 @@ private:
// the problem
const auto& problem = globalFvGeometry().problem_();
// make the scv
auto eIdx = problem.elementMapper().index(element);
scvs_.emplace_back(element.geometry(), eIdx);
scvIndices_.emplace_back(eIdx);
// get data on the scv faces
const auto& scvFaceIndices = globalFvGeometry().scvfIndicesOfScv(eIdx);
const auto& neighborVolVarIndices = globalFvGeometry().neighborVolVarIndices(eIdx);
......@@ -347,8 +346,28 @@ private:
scvfIndices_.shrink_to_fit();
}
template<typename Seed>
void handleInteractionRegion(const Seed& ivSeed)
{
// the element index of the element which we are binding to
auto eIdx = globalFvGeometry().problem_().elementMapper().index(*elementPtr_);
// the scvf indices in the interaction region
auto scvfIndices = ivSeed.globalScvfIndices();
// make the scvs/scvfs for each element in the interaction region
for (auto eIdxGlobal : ivSeed.globalScvIndices())
{
if (eIdxGlobal != eIdx)
{
auto element = globalFvGeometry().element(eIdxGlobal);
makeNeighborGeometries(element, eIdxGlobal, scvfIndices);
}
}
}
//! create the necessary scvs and scvfs of the neighbor elements to the bound elements
template<class IndexVector>
template<typename IndexVector>
void makeNeighborGeometries(const Element& element, IndexType eIdxGlobal, const IndexVector& ivScvfs)
{
// create the neighbor scv if it doesn't exist yet
......
......@@ -23,7 +23,6 @@
#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
......
......@@ -47,6 +47,7 @@ class MpfaHelperBase<TypeTag, MpfaMethods::oMethod, 2>
using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
using InteractionVolume = typename GET_PROP_TYPE(TypeTag, InteractionVolume);
using BoundaryInteractionVolume = typename GET_PROP_TYPE(TypeTag, BoundaryInteractionVolume);
static const int dim = GridView::dimension;
static const int dimWorld = GridView::dimensionworld;
......@@ -58,6 +59,10 @@ class MpfaHelperBase<TypeTag, MpfaMethods::oMethod, 2>
using ScvSeed = typename InteractionVolumeSeed::LocalScvSeed;
using ScvfSeed = typename InteractionVolumeSeed::LocalScvfSeed;
using BoundaryInteractionVolumeSeed = typename BoundaryInteractionVolume::Seed;
using BoundaryScvSeed = typename BoundaryInteractionVolumeSeed::LocalScvSeed;
using BoundaryScvfSeed = typename BoundaryInteractionVolumeSeed::LocalScvfSeed;
using GlobalIndexSet = std::vector<GlobalIndexType>;
using LocalIndexSet = std::vector<LocalIndexType>;
......@@ -84,13 +89,13 @@ public:
return InteractionVolumeSeed(std::move(scvSeeds), std::move(scvfSeeds), false);
}
static InteractionVolumeSeed makeBoundaryInteractionVolumeSeed(const Problem& problem,
const Element& element,
const FVElementGeometry& fvGeometry,
const SubControlVolumeFace& scvf)
static BoundaryInteractionVolumeSeed makeBoundaryInteractionVolumeSeed(const Problem& problem,
const Element& element,
const FVElementGeometry& fvGeometry,
const SubControlVolumeFace& scvf)
{
std::vector<ScvSeed> scvSeeds;
std::vector<ScvfSeed> scvfSeeds;
std::vector<BoundaryScvSeed> scvSeeds;
std::vector<BoundaryScvfSeed> scvfSeeds;
// reserve sufficient memory
scvSeeds.reserve(20);
......@@ -102,10 +107,11 @@ public:
scvSeeds.shrink_to_fit();
scvfSeeds.shrink_to_fit();
return InteractionVolumeSeed(std::move(scvSeeds), std::move(scvfSeeds), true);
return BoundaryInteractionVolumeSeed(std::move(scvSeeds), std::move(scvfSeeds), true);
}
private:
template<typename ScvSeed, typename ScvfSeed>
static void fillEntitySeeds_(std::vector<ScvSeed>& scvSeeds,
std::vector<ScvfSeed>& scvfSeeds,
const Problem& problem,
......@@ -146,7 +152,7 @@ private:
}
// clockwise rotation and construction of local scv & scv face entities
template<class ScvfPointerVector>
template<typename ScvSeed, typename ScvfSeed, typename ScvfPointerVector>
static void performRotation_(const Problem& problem,
const ScvfPointerVector& scvfVector,
std::vector<ScvSeed>& scvSeeds,
......@@ -257,6 +263,7 @@ class MpfaHelperBase<TypeTag, MpfaMethods::oMethod, 3>
using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
using InteractionVolume = typename GET_PROP_TYPE(TypeTag, InteractionVolume);
using BoundaryInteractionVolume = typename GET_PROP_TYPE(TypeTag, BoundaryInteractionVolume);
static const int dim = GridView::dimension;
static const int dimWorld = GridView::dimensionworld;
......@@ -268,6 +275,10 @@ class MpfaHelperBase<TypeTag, MpfaMethods::oMethod, 3>
using ScvSeed = typename InteractionVolumeSeed::LocalScvSeed;
using ScvfSeed = typename InteractionVolumeSeed::LocalScvfSeed;
using BoundaryInteractionVolumeSeed = typename BoundaryInteractionVolume::Seed;
using BoundaryScvSeed = typename BoundaryInteractionVolumeSeed::LocalScvSeed;
using BoundaryScvfSeed = typename BoundaryInteractionVolumeSeed::LocalScvfSeed;
using GlobalIndexSet = std::vector<GlobalIndexType>;
using LocalIndexSet = std::vector<LocalIndexType>;
......@@ -299,13 +310,13 @@ public:
return InteractionVolumeSeed(std::move(scvSeeds), std::move(scvfSeeds), false);
}
static InteractionVolumeSeed makeBoundaryInteractionVolumeSeed(const Problem& problem,
const Element& element,
const FVElementGeometry& fvGeometry,
const SubControlVolumeFace& scvf)
static BoundaryInteractionVolumeSeed makeBoundaryInteractionVolumeSeed(const Problem& problem,
const Element& element,
const FVElementGeometry& fvGeometry,
const SubControlVolumeFace& scvf)
{
std::vector<ScvSeed> scvSeeds;
std::vector<ScvfSeed> scvfSeeds;
std::vector<BoundaryScvSeed> scvSeeds;
std::vector<BoundaryScvfSeed> scvfSeeds;
// reserve sufficient memory
scvSeeds.reserve(50);
......@@ -321,10 +332,11 @@ public:
scvSeeds.shrink_to_fit();
scvfSeeds.shrink_to_fit();
return InteractionVolumeSeed(std::move(scvSeeds), std::move(scvfSeeds), true);
return BoundaryInteractionVolumeSeed(std::move(scvSeeds), std::move(scvfSeeds), true);
}
private:
template<typename ScvSeed, typename ScvfSeed>
static void fillEntitySeeds_(std::vector<ScvSeed>& scvSeeds,
std::vector<ScvfSeed>& scvfSeeds,
const Problem& problem,
......
......@@ -26,16 +26,15 @@
#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>
#include "localsubcontrolentityseeds.hh"
#include "interactionvolumeseed.hh"
#include "localsubcontrolentities.hh"
namespace Dumux
{
//! Specialization of the interaction volume traits clas
//! Specialization of the interaction volume traits class for the mpfa-o method
template<class TypeTag>
class CCMpfaOInteractionVolumeTraits
{
......@@ -65,9 +64,7 @@ public:
using LocalScvType = CCMpfaOLocalScv<TypeTag>;
using LocalScvfType = CCMpfaOLocalScvf<TypeTag>;
using LocalScvSeedType = CCMpfaOLocalScvSeed<GlobalIndexType, LocalIndexType>;
using LocalScvfSeedType = CCMpfaOLocalScvfSeed<GlobalIndexType, LocalIndexType>;
using Seed = CCMpfaInteractionVolumeSeed<LocalScvSeedType, LocalScvfSeedType>;
using Seed = CCMpfaOInteractionVolumeSeed<GlobalIndexType, LocalIndexType>;
};
//! Forward declaration of the mpfa-o interaction volume
......
......@@ -20,11 +20,11 @@
* \file
* \brief Base class for interaction volume seeds of mpfa methods.
*/
#ifndef DUMUX_DISCRETIZATION_CC_MPFA_INTERACTIONVOLUMESEED_HH
#define DUMUX_DISCRETIZATION_CC_MPFA_INTERACTIONVOLUMESEED_HH
#ifndef DUMUX_DISCRETIZATION_CC_MPFA_O_INTERACTIONVOLUMESEED_HH
#define DUMUX_DISCRETIZATION_CC_MPFA_O_INTERACTIONVOLUMESEED_HH
#include <dumux/implicit/cellcentered/mpfa/properties.hh>
#include "facetypes.hh"
#include "localsubcontrolentityseeds.hh"
namespace Dumux
{
......@@ -33,19 +33,18 @@ namespace Dumux
* \ingroup Mpfa
* \brief Base class for the interaction volume seed of mpfa methods
*/
template<class ScvSeedType, class ScvfSeedType>
class CCMpfaInteractionVolumeSeed
template<typename G, typename L>
class CCMpfaOInteractionVolumeSeed
{
using GlobalScvIdxType = typename ScvSeedType::GlobalIndexType;
using GlobalScvfIdxType = typename ScvfSeedType::GlobalIndexType;
using GlobalIndexType = G;
public:
using LocalScvSeed = ScvSeedType;
using LocalScvfSeed = ScvfSeedType;
using LocalScvSeed = CCMpfaOLocalScvSeed<G, L>;
using LocalScvfSeed = CCMpfaOLocalScvfSeed<G, L>;
CCMpfaInteractionVolumeSeed(std::vector<LocalScvSeed>&& scvSeeds,
std::vector<LocalScvfSeed>&& scvfSeeds,
bool onBoundary)
CCMpfaOInteractionVolumeSeed(std::vector<LocalScvSeed>&& scvSeeds,
std::vector<LocalScvfSeed>&& scvfSeeds,
bool onBoundary)
: onBoundary_(onBoundary),
scvSeeds_(std::move(scvSeeds)),
scvfSeeds_(std::move(scvfSeeds))
......@@ -60,9 +59,9 @@ public:
const std::vector<LocalScvfSeed>& scvfSeeds() const
{ return scvfSeeds_; }
std::vector<GlobalScvIdxType> globalScvIndices() const
std::vector<GlobalIndexType> globalScvIndices() const
{
std::vector<GlobalScvIdxType> globalIndices;
std::vector<GlobalIndexType> globalIndices;
globalIndices.reserve(scvSeeds().size());
for (auto&& localScvSeed : scvSeeds())
......@@ -71,9 +70,9 @@ public:
return globalIndices;
}
std::vector<GlobalScvfIdxType> globalScvfIndices() const
std::vector<GlobalIndexType> globalScvfIndices() const
{
std::vector<GlobalScvfIdxType> globalIndices;
std::vector<GlobalIndexType> globalIndices;
globalIndices.reserve(scvfSeeds().size() * 2);
for (auto&& localScvfSeed : scvfSeeds())
......
......@@ -29,7 +29,6 @@
#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>
......
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