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

[mpfa][fvGeometry] implement classes for caching disabled

parent 9f7a78cd
......@@ -30,6 +30,8 @@
#include <dumux/discretization/scvandscvfiterators.hh>
#include <dumux/implicit/cellcentered/mpfa/properties.hh>
#include "mpfageometryhelper.hh"
namespace Dumux
{
/*!
......@@ -143,6 +145,7 @@ class CCMpfaFVElementGeometry<TypeTag, false>
{
using ThisType = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
using IndexType = typename GridView::IndexSet::IndexType;
using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
......@@ -153,6 +156,12 @@ class CCMpfaFVElementGeometry<TypeTag, false>
using ScvIterator = Dumux::ScvIterator<SubControlVolume, std::vector<IndexType>, ThisType>;
using ScvfIterator = Dumux::ScvfIterator<SubControlVolumeFace, std::vector<IndexType>, ThisType>;
static const int dim = GridView::dimension;
using CoordScalar = typename GridView::ctype;
using ReferenceElements = typename Dune::ReferenceElements<CoordScalar, dim>;
using MpfaGeometryHelper = Dumux::MpfaGeometryHelper<GridView, dim>;
public:
//! Constructor
CCMpfaFVElementGeometry(const GlobalFVGeometry& globalFvGeometry)
......@@ -216,13 +225,47 @@ public:
void bind(const Element& element)
{
bindElement(element);
for (const auto& intersection : intersections(globalFvGeometry().gridView(), element))
// reserve memory
auto numNeighbors = globalFvGeometry().problem_().model().stencils(element).neighborStencil().size();
auto numNeighborScvf = numNeighbors*dim*(2*dim-2);
neighborScvs_.reserve(numNeighbors);
neighborScvIndices_.reserve(numNeighbors);
neighborScvfIndices_.reserve(numNeighborScvf);
neighborScvfs_.reserve(numNeighborScvf);
// build scvfs connected to the interaction regions around the element vertices
std::vector<IndexType> finishedVertices;
finishedVertices.reserve(element.geometry().corners());
const auto eIdx = globalFvGeometry().problem_().elementMapper().index(element);
for (auto&& scvf : scvfs_)
{
neighborScvs_.reserve(element.subEntities(1));
neighborScvfIndices_.reserve(element.subEntities(1));
if (intersection.neighbor())
makeNeighborGeometries(intersection.outside());
if (existsLocalIndex(scvf.vertexIndex(), finishedVertices).second)
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 (eIdxGlobal != eIdx)
{
auto element = globalFvGeometry().element(eIdxGlobal);
makeNeighborGeometries(element, eIdxGlobal, scvfIndices);
}
}
finishedVertices.push_back(scvf.vertexIndex());
}
// free unused memory
neighborScvs_.shrink_to_fit();
neighborScvIndices_.shrink_to_fit();
neighborScvfIndices_.shrink_to_fit();
neighborScvfs_.shrink_to_fit();
}
//! Binding of an element preparing the geometries only inside the element
......@@ -242,62 +285,125 @@ private:
//! create scvs and scvfs of the bound element
void makeElementGeometries(const Element& element)
{
auto eIdx = globalFvGeometry().problem_().elementMapper().index(element);
// the problem
const auto& problem = globalFvGeometry().problem_();
auto eIdx = problem.elementMapper().index(element);
scvs_.emplace_back(element.geometry(), eIdx);
scvIndices_.emplace_back(eIdx);
const auto& scvFaceIndices = globalFvGeometry().scvfIndicesOfScv(eIdx);
const auto& neighborVolVarIndices = globalFvGeometry().neighborVolVarIndices(eIdx);
// the element geometry and reference element
auto g = element.geometry();
const auto& referenceElement = ReferenceElements::general(g.type());
// The geometry helper class
MpfaGeometryHelper geomHelper(g);
// the quadrature point to be used on the scvf
const Scalar q = GET_PARAM_FROM_GROUP(TypeTag, Scalar, Mpfa, Q);
// reserve memory for the scv faces
auto numLocalScvf = scvFaceIndices.size();
scvfIndices_.reserve(numLocalScvf);
scvfs_.reserve(numLocalScvf);
int scvfCounter = 0;
for (const auto& intersection : intersections(globalFvGeometry().gridView(), element))
for (const auto& is : intersections(globalFvGeometry().gridView(), element))
{
if (intersection.neighbor() || intersection.boundary())
auto isGeometry = is.geometry();
// make the scv faces of the intersection
for (unsigned int faceScvfIdx = 0; faceScvfIdx < isGeometry.corners(); ++faceScvfIdx)
{
scvfs_.emplace_back(intersection,
intersection.geometry(),
// get the global vertex index the scv face is connected to (mpfa-o method does not work for hanging nodes!)
const auto vIdxLocal = referenceElement.subEntity(is.indexInInside(), 1, faceScvfIdx, dim);
const auto vIdxGlobal = problem.vertexMapper().subIndex(element, vIdxLocal, dim);
// do not build scvfs connected to a processor boundary
if (globalFvGeometry().isGhostVertex(vIdxGlobal))
continue;
scvfs_.emplace_back(geomHelper,
geomHelper.getScvfCorners(isGeometry, faceScvfIdx),
is.centerUnitOuterNormal(),
vIdxGlobal,
scvFaceIndices[scvfCounter],
std::vector<IndexType>({eIdx, neighborVolVarIndices[scvfCounter]})
std::array<IndexType, 2>({{eIdx, neighborVolVarIndices[scvfCounter]}}),
q,
is.boundary()
);
scvfIndices_.emplace_back(scvFaceIndices[scvfCounter]);
scvfCounter++;
}
}
// free unused memory
scvfs_.shrink_to_fit();
scvfIndices_.shrink_to_fit();
}
//! create the necessary scvs and scvfs of the neighbor elements to the bound elements
void makeNeighborGeometries(const Element& element)
template<class IndexVector>
void makeNeighborGeometries(const Element& element, IndexType eIdxGlobal, const IndexVector& ivScvfs)
{
// create the neighbor scv
auto eIdx = globalFvGeometry().problem_().elementMapper().index(element);
neighborScvs_.emplace_back(element.geometry(), eIdx);
neighborScvIndices_.push_back(eIdx);
// create the neighbor scv if it doesn't exist yet
auto pair = existsLocalIndex(eIdxGlobal, neighborScvIndices_);
if (!pair.second)
{
neighborScvs_.emplace_back(element.geometry(), eIdxGlobal);
neighborScvIndices_.push_back(eIdxGlobal);
}
const auto& scvFaceIndices = globalFvGeometry().scvfIndicesOfScv(eIdx);
const auto& neighborVolVarIndices = globalFvGeometry().neighborVolVarIndices(eIdx);
const auto& scvFaceIndices = globalFvGeometry().scvfIndicesOfScv(eIdxGlobal);
const auto& neighborVolVarIndices = globalFvGeometry().neighborVolVarIndices(eIdxGlobal);
// the element geometry and reference element
auto g = element.geometry();
const auto& referenceElement = ReferenceElements::general(g.type());
// The geometry helper class
MpfaGeometryHelper geomHelper(g);
// the quadrature point to be used on the scvf
const Scalar q = GET_PARAM_FROM_GROUP(TypeTag, Scalar, Mpfa, Q);
int scvfCounter = 0;
for (const auto& intersection : intersections(globalFvGeometry().gridView(), element))
for (const auto& is : intersections(globalFvGeometry().gridView(), element))
{
if (intersection.neighbor())
auto isGeometry = is.geometry();
// make the scv faces of the intersection
for (unsigned int faceScvfIdx = 0; faceScvfIdx < isGeometry.corners(); ++faceScvfIdx)
{
// only create subcontrol faces where the outside element is the bound element
if (intersection.outside() == *elementPtr_)
// get the global vertex index the scv face is connected to (mpfa-o method does not work for hanging nodes!)
const auto vIdxLocal = referenceElement.subEntity(is.indexInInside(), 1, faceScvfIdx, dim);
const auto vIdxGlobal = globalFvGeometry().problem_().vertexMapper().subIndex(element, vIdxLocal, dim);
// do not build scvfs connected to a processor boundary
if (globalFvGeometry().isGhostVertex(vIdxGlobal))
continue;
// if the actual global scvf index is in the container, make scvf
if (existsLocalIndex(scvFaceIndices[scvfCounter], ivScvfs).second)
{
neighborScvfs_.emplace_back(intersection,
intersection.geometry(),
neighborScvfs_.emplace_back(geomHelper,
geomHelper.getScvfCorners(isGeometry, faceScvfIdx),
is.centerUnitOuterNormal(),
vIdxGlobal,
scvFaceIndices[scvfCounter],
std::vector<IndexType>({eIdx, neighborVolVarIndices[scvfCounter]})
std::array<IndexType, 2>({{eIdxGlobal, neighborVolVarIndices[scvfCounter]}}),
q,
is.boundary()
);
neighborScvfIndices_.push_back(scvFaceIndices[scvfCounter]);
neighborScvfIndices_.emplace_back(scvFaceIndices[scvfCounter]);
}
// always increase local counter
scvfCounter++;
}
else if (intersection.boundary())
{
scvfCounter++;
}
}
......@@ -312,6 +418,17 @@ private:
}
const std::pair<IndexType, bool> existsLocalIndex(const IndexType idx,
const std::vector<IndexType>& indices) const
{
auto it = std::find(indices.begin(), indices.end(), idx);
if (it != indices.end())
return std::make_pair<IndexType, bool>(std::distance(indices.begin(), it), true);
else
return std::make_pair<IndexType, bool>(0, false);
}
void clear()
{
scvIndices_.clear();
......
......@@ -335,6 +335,222 @@ private:
std::vector<std::vector<IndexType>> scvfIndicesOfScv_;
std::vector<bool> boundaryVertices_;
IndexType numBoundaryScvf_;
// the global interaction volume seeds
GlobalInteractionVolumeSeeds globalInteractionVolumeSeeds_;
};
// specialization in case the FVElementGeometries are not stored
template<class TypeTag>
class CCMpfaGlobalFVGeometryBase<TypeTag, false>
{
//! The local fvGeometry needs access to the problem
friend typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
using Implementation = typename GET_PROP_TYPE(TypeTag, GlobalFVGeometry);
using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
using GlobalInteractionVolumeSeeds = typename GET_PROP_TYPE(TypeTag, GlobalInteractionVolumeSeeds);
using BoundaryInteractionVolume = typename GET_PROP_TYPE(TypeTag, BoundaryInteractionVolume);
using InteractionVolume = typename GET_PROP_TYPE(TypeTag, InteractionVolume);
using InteractionVolumeSeed = typename InteractionVolume::Seed;
using BoundaryInteractionVolumeSeed = typename BoundaryInteractionVolume::Seed;
using Element = typename GridView::template Codim<0>::Entity;
using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
using CoordScalar = typename GridView::ctype;
using IndexType = typename GridView::IndexSet::IndexType;
using LocalIndexType = typename InteractionVolume::LocalIndexType;
static const int dim = GridView::dimension;
static const int dimWorld = GridView::dimensionworld;
using DimVector = Dune::FieldVector<Scalar, dimWorld>;
using ReferenceElements = typename Dune::ReferenceElements<CoordScalar, dim>;
using MpfaGeometryHelper = Dumux::MpfaGeometryHelper<GridView, dim>;
public:
//! Constructor
CCMpfaGlobalFVGeometryBase(const GridView gridView)
: gridView_(gridView), elementMap_(gridView), globalInteractionVolumeSeeds_(gridView_) {}
//! The total number of sub control volumes
std::size_t numScv() const
{ return numScvs_; }
//! The total number of sub control volume faces
std::size_t numScvf() const
{ return numScvf_; }
//! The total number of boundary sub control volume faces
std::size_t numBoundaryScvf() const
{ return numBoundaryScvf_; }
// Get an element from a sub control volume contained in it
Element element(const SubControlVolume& scv) const
{ return elementMap_.element(scv.elementIndex()); }
// Get an element from a global element index
Element element(IndexType eIdx) const
{ return elementMap_.element(eIdx); }
//! Return the gridView this global object lives on
const GridView& gridView() const
{ return gridView_; }
//! Get the inner interaction volume seed corresponding to an scvf
const InteractionVolumeSeed& interactionVolumeSeed(const SubControlVolumeFace& scvf) const
{ return globalInteractionVolumeSeeds_.seed(scvf); }
//! Get the boundary interaction volume seed corresponding to an scvf
const BoundaryInteractionVolumeSeed& boundaryInteractionVolumeSeed(const SubControlVolumeFace& scvf) const
{ return globalInteractionVolumeSeeds_.boundarySeed(scvf); }
//! Returns whether or not a scvf touches the boundary (has to be called before getting an interaction volume)
bool scvfTouchesBoundary(const SubControlVolumeFace& scvf) const
{ return boundaryVertices_[scvf.vertexIndex()]; }
//! Returns whether or not a vertex is on a processor boundary
bool isGhostVertex(const IndexType vIdxGlobal) const
{ return ghostVertices_[vIdxGlobal]; }
//! update all fvElementGeometries (do this again after grid adaption)
void update(const Problem& problem)
{
problemPtr_ = &problem;
// clear containers (necessary after grid refinement)
scvfIndicesOfScv_.clear();
neighborVolVarIndices_.clear();
boundaryVertices_.clear();
elementMap_.clear();
// reserve memory or resize the containers
numScvs_ = gridView_.size(0);
numScvf_ = 0;
numBoundaryScvf_ = 0;
elementMap_.resize(numScvs_);
scvfIndicesOfScv_.resize(numScvs_);
neighborVolVarIndices_.resize(numScvs_);
boundaryVertices_.resize(gridView_.size(dim), false);
// find vertices on processor boundaries
ghostVertices_ = CCMpfaGlobalFVGeometryHelper<TypeTag>::findGhostVertices(problem, gridView_);
// Store necessary info on SCVs and SCV faces
for (const auto& element : elements(gridView_))
{
auto eIdx = problem.elementMapper().index(element);
// fill the element map with seeds
elementMap_[eIdx] = element.seed();
// the element geometry and reference element
auto elemGeometry = element.geometry();
const auto& referenceElement = ReferenceElements::general(elemGeometry.type());
// The geometry helper class
MpfaGeometryHelper geomHelper(elemGeometry);
// the element-wise index sets for finite volume geometry
auto numLocalFaces = geomHelper.getNumLocalScvfs();
std::vector<IndexType> scvfsIndexSet(numLocalFaces);
std::vector<IndexType> neighborVolVarIndexSet(numLocalFaces);
IndexType localFaceIdx = 0;
// construct the sub control volume faces
for (const auto& intersection : intersections(gridView_, element))
{
// get some of the intersection bools
bool boundary = intersection.boundary();
bool neighbor = intersection.neighbor();
// determine the outside volvar idx
IndexType nIdx;
if (neighbor)
nIdx = problem.elementMapper().index(intersection.outside());
else if (boundary)
nIdx = numScvs_ + numBoundaryScvf_++;
// make the scv faces of the intersection
for (unsigned int faceScvfIdx = 0; faceScvfIdx < intersection.geometry().corners(); ++faceScvfIdx)
{
// get the global vertex index the scv face is connected to (mpfa-o method does not work for hanging nodes!)
const auto vIdxLocal = referenceElement.subEntity(intersection.indexInInside(), 1, faceScvfIdx, dim);
const auto vIdxGlobal = problem.vertexMapper().subIndex(element, vIdxLocal, dim);
// do not build scvfs connected to a processor boundary
if (ghostVertices_[vIdxGlobal])
continue;
// store info on which vertices are on the domain boundary
if (boundary)
boundaryVertices_[vIdxGlobal] = true;
// store information on the scv face
scvfsIndexSet[localFaceIdx] = numScvf_++;
neighborVolVarIndexSet[localFaceIdx++] = nIdx;
}
}
// store the sets of indices in the data container
scvfIndicesOfScv_[eIdx] = scvfsIndexSet;
neighborVolVarIndices_[eIdx] = neighborVolVarIndexSet;
}
// the number of actual boundary scvf is two times the number of boundary intersections
numBoundaryScvf_ *= 2;
// Initialize the interaction volume seeds
globalInteractionVolumeSeeds_.update(problem, boundaryVertices_);
}
const std::vector<IndexType>& scvfIndicesOfScv(IndexType scvIdx) const
{ return scvfIndicesOfScv_[scvIdx]; }
const std::vector<IndexType>& neighborVolVarIndices(IndexType scvIdx) const
{ return neighborVolVarIndices_[scvIdx]; }
/*!
* \brief Return a local restriction of this global object
* The local object is only functional after calling its bind/bindElement method
* This is a free function that will be found by means of ADL
*/
friend inline FVElementGeometry localView(const Implementation& global)
{ return FVElementGeometry(global); }
private:
//! Get the inner interaction volume seed corresponding to an scvf
const InteractionVolumeSeed& interactionVolumeSeed(const IndexType vIdxGlobal) const
{ return globalInteractionVolumeSeeds_.seed(vIdxGlobal); }
//! Get the boundary interaction volume seed corresponding to an scvf
const BoundaryInteractionVolumeSeed& boundaryInteractionVolumeSeed(const IndexType vIdxGlobal) const
{ return globalInteractionVolumeSeeds_.boundarySeed(vIdxGlobal); }
const Problem& problem_() const
{ return *problemPtr_; }
const Problem* problemPtr_;
GridView gridView_;
// Information on the global number of geometries
IndexType numScvs_;
IndexType numScvf_;
IndexType numBoundaryScvf_;
// vectors that store the global data
Dumux::ElementMap<GridView> elementMap_;
std::vector<std::vector<IndexType>> scvfIndicesOfScv_;
std::vector<std::vector<IndexType>> neighborVolVarIndices_;
std::vector<bool> boundaryVertices_;
std::vector<bool> ghostVertices_;
// the global interaction volume seeds
GlobalInteractionVolumeSeeds globalInteractionVolumeSeeds_;
};
......
......@@ -24,6 +24,7 @@
#define DUMUX_DISCRETIZATION_MPFA_GLOBALINTERACTIONVOLUMESEEDS_BASE_HH
#include <dumux/implicit/cellcentered/mpfa/properties.hh>
#include "globalfvgeometry.hh"
namespace Dumux
{
......@@ -34,6 +35,11 @@ namespace Dumux
template<class TypeTag>
class CCMpfaGlobalInteractionVolumeSeedsBase
{
//! the globalFvGeometry and the FVElementGeometry needs access to the private seed return functions
//! in case caching is disabled
friend CCMpfaGlobalFVGeometryBase<TypeTag, GET_PROP_VALUE(TypeTag, EnableGlobalFVGeometryCache)>;
friend typename GET_PROP_TYPE(TypeTag, GlobalFVGeometry);
using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
using Helper = typename GET_PROP_TYPE(TypeTag, MpfaHelper);
......@@ -48,7 +54,6 @@ class CCMpfaGlobalInteractionVolumeSeedsBase
using LocalIndexType = typename InteractionVolume::LocalIndexType;
static const int dim = GridView::dimension;
static const int numEq = GET_PROP_VALUE(TypeTag, NumEq);
public:
CCMpfaGlobalInteractionVolumeSeedsBase(const GridView gridView) : gridView_(gridView) {}
......@@ -75,6 +80,7 @@ public:
{ return boundarySeeds_[scvfIndexMap_[scvf.index()]]; }
private:
void initializeSeeds_(std::vector<bool>& boundaryVertices)
{
const auto numScvf = problem_().model().globalFvGeometry().numScvf();
......@@ -89,7 +95,7 @@ private:
for (const auto& element : elements(gridView_))
{
auto fvGeometry = localView(problem_().model().globalFvGeometry());
fvGeometry.bind(element);
fvGeometry.bindElement(element);
for (const auto& scvf : scvfs(fvGeometry))
{
// skip the rest if we already handled this face
......@@ -98,30 +104,26 @@ private:
if (boundaryVertices[scvf.vertexIndex()])
{
// the boundary interaction volume seed
auto seed = Helper::makeBoundaryInteractionVolumeSeed(problem_(), element, fvGeometry, scvf);
// make the boundary interaction volume seed
boundarySeeds_.emplace_back(Helper::makeBoundaryInteractionVolumeSeed(problem_(), element, fvGeometry, scvf));
// update the index map entries for the global scv faces in the interaction volume
for (const auto& localScvfSeed : seed.scvfSeeds())
for (const auto scvfIdxGlobal : localScvfSeed.globalScvfIndices())
scvfIndexMap_[scvfIdxGlobal] = boundarySeedIndex;
for (auto scvfIdxGlobal : boundarySeeds_.back().globalScvfIndices())
scvfIndexMap_[scvfIdxGlobal] = boundarySeedIndex;
// store interaction volume and increment counter
boundarySeeds_.emplace_back(std::move(seed));
// increment counter
boundarySeedIndex++;
}
else
{
// the inner interaction volume seed
auto seed = Helper::makeInnerInteractionVolumeSeed(problem_(), element, fvGeometry, scvf);
// make the inner interaction volume seed
seeds_.emplace_back(Helper::makeInnerInteractionVolumeSeed(problem_(), element, fvGeometry, scvf));
// update the index map entries for the global scv faces in the interaction volume
for (const auto& localScvf : seed.scvfSeeds())
for (const auto scvfIdxGlobal : localScvf.globalScvfIndices())
scvfIndexMap_[scvfIdxGlobal] = seedIndex;
for (auto scvfIdxGlobal : seeds_.back().globalScvfIndices())
scvfIndexMap_[scvfIdxGlobal] = seedIndex;
// store interaction volume and increment counter
seeds_.emplace_back(std::move(seed));
// increment counter
seedIndex++;
}
}
......
......@@ -36,6 +36,9 @@ namespace Dumux
template<class ScvSeedType, class ScvfSeedType>
class CCMpfaInteractionVolumeSeed
{
using GlobalScvIdxType = typename ScvSeedType::GlobalIndexType;
using GlobalScvfIdxType = typename ScvfSeedType::GlobalIndexType;
public:
using LocalScvSeed = ScvSeedType;
using LocalScvfSeed = ScvfSeedType;
......@@ -57,6 +60,30 @@ public:
const std::vector<LocalScvfSeed>& scvfSeeds() const
{ return scvfSeeds_; }
std::vector<GlobalScvIdxType> globalScvIndices() const
{
std::vector<GlobalScvIdxType> globalIndices;
globalIndices.reserve(scvSeeds().size());
for (auto&& localScvSeed : scvSeeds())