Commit 056bc202 authored by Dennis Gläser's avatar Dennis Gläser
Browse files

[mpfa] delete (currently) not working schemes other than the o scheme

parent 2823993a
// -*- 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 the finite volume geometry vector for box models
* This builds up the sub control volumes and sub control volume faces
* for each element.
*/
#ifndef DUMUX_DISCRETIZATION_CC_MPFA_O_HYBRIDFPS_FV_GRID_GEOMETRY_HH
#define DUMUX_DISCRETIZATION_CC_MPFA_O_HYBRIDFPS_FV_GRID_GEOMETRY_HH
#include <dune/geometry/multilineargeometry.hh>
#include <dune/geometry/referenceelements.hh>
#include <dumux/common/elementmap.hh>
#include <dumux/discretization/cellcentered/mpfa/mpfageometryhelper.hh>
namespace Dumux
{
/*!
* \ingroup ImplicitModel
* \brief Base class for the finite volume geometry vector for mpfa models
* This builds up the sub control volumes and sub control volume faces
* for each element.
*/
template<class TypeTag, bool EnableFVElementGeometryCache>
class CCMpfaOHybridFpsFVGridGeometry
{};
// specialization in case the FVElementGeometries are stored
template<class TypeTag>
class CCMpfaOHybridFpsFVGridGeometry<TypeTag, true>
{
//! The local class needs access to the scv, scvfs and the fv element geometry
//! as they are globally cached
friend typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
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
CCMpfaOHybridFpsFVGridGeometry(const GridView gridView)
: gridView_(gridView), elementMap_(gridView), globalInteractionVolumeSeeds_(gridView) {}
//! The total number of sub control volumes
std::size_t numScv() const
{ return scvs_.size(); }
//! The total number of sun control volume faces
std::size_t numScvf() const
{ return scvfs_.size(); }
//! The total number of boundary sub control volume faces
std::size_t numBoundaryScvf() const
{ return numBoundaryScvf_; }
//! The total number of degrees of freedom
std::size_t numDofs() const
{ return this->gridView().size(0); }
// 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); }
//! 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 LocalIndexType eqIdx) const
{ return globalInteractionVolumeSeeds_.boundarySeed(scvf, eqIdx); }
//! 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 fpsVertices_[scvf.vertexIndex()]; }
//! update all fvElementGeometries (do this again after grid adaption)
void update(const Problem& problem)
{
problemPtr_ = &problem;
scvfs_.clear();
// reserve memory
IndexType numScvs = numDofs();
scvs_.resize(numScvs);
scvfs_.reserve(getNumScvf_());
scvfIndicesOfScv_.resize(numScvs);
fpsVertices_.resize(gridView_.size(dim), false);
boundaryVertices_.resize(gridView_.size(dim), false);
elementMap_.resize(numScvs);
// the quadrature point to be used on the scvf
const Scalar q = GET_PARAM_FROM_GROUP(TypeTag, Scalar, Mpfa, Q);
// Build the SCVs and SCV faces
IndexType scvfIdx = 0;
numBoundaryScvf_ = 0;
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());
// create sub control volume for this element
scvs_[eIdx] = SubControlVolume(std::move(elemGeometry), eIdx);
// The geometry helper class
MpfaGeometryHelper geomHelper(elemGeometry);
// The local scvf index set
std::vector<IndexType> scvfIndexSet;
scvfIndexSet.reserve(geomHelper.getNumLocalScvfs());
// determine if element has to be handled by the fps scheme
bool isFps = problem.spatialParams().isFpsElement(element);
// construct the sub control volume faces
for (const auto& intersection : intersections(gridView_, element))
{
// get the intersection geometry and some of its bools
auto isGeometry = intersection.geometry();
bool boundary = intersection.boundary();
bool neighbor = intersection.neighbor();
if (neighbor &&
element.partitionType() == Dune::GhostEntity &&
intersection.outside().partitionType() == Dune::GhostEntity)
continue;
// determine the outside volvar idx
IndexType nIdx;
if (neighbor)
nIdx = problem.elementMapper().index(intersection.outside());
else if (boundary)
nIdx = numScvs + numBoundaryScvf_++;
else
continue;
// make the scv faces of the intersection
for (unsigned int faceScvfIdx = 0; faceScvfIdx < isGeometry.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);
if (boundary || isFps || (!boundary && problem.spatialParams().isFpsElement(intersection.outside())))
fpsVertices_[vIdxGlobal] = true;
if (boundary)
boundaryVertices_[vIdxGlobal] = true;
// make the scv face
scvfIndexSet.push_back(scvfIdx);
scvfs_.emplace_back(SubControlVolumeFace(geomHelper,
geomHelper.getScvfCorners(isGeometry, faceScvfIdx),
intersection.centerUnitOuterNormal(),
vIdxGlobal,
scvfIdx,
std::array<IndexType, 2>({{eIdx, nIdx}}),
q,
boundary
));
// increment scvf counter
scvfIdx++;
}
}
// Save the scvf indices belonging to this scv to build up fv element geometries fast
scvfIndicesOfScv_[eIdx] = scvfIndexSet;
}
// Initialize the interaction volume seeds, this will also initialize the vector of boundary vertices
globalInteractionVolumeSeeds_.update(problem, fpsVertices_, boundaryVertices_);
}
/*!
* \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 CCMpfaOHybridFpsFVGridGeometry& global)
{ return FVElementGeometry(global); }
//! Get a sub control volume with a global scv index
const SubControlVolume& scv(IndexType scvIdx) const
{ return scvs_[scvIdx]; }
//! Get a sub control volume face with a global scvf index
const SubControlVolumeFace& scvf(IndexType scvfIdx) const
{ return scvfs_[scvfIdx]; }
//! Get the sub control volume face indices of an scv by global index
const std::vector<IndexType>& scvfIndicesOfScv(IndexType scvIdx) const
{ return scvfIndicesOfScv_[scvIdx]; }
template<int d = dim>
typename std::enable_if<d == 2, IndexType>::type getNumScvf_() const
{
Dune::GeometryType triangle, quadrilateral;
triangle.makeTriangle();
quadrilateral.makeQuadrilateral();
return gridView_.size(triangle)*6 + gridView_.size(quadrilateral)*8;
}
template<int d = dim>
typename std::enable_if<d == 3, IndexType>::type getNumScvf_() const
{
Dune::GeometryType simplex, pyramid, prism, cube;
simplex.makeTetrahedron();
pyramid.makePyramid();
prism.makePrism();
cube.makeHexahedron();
return gridView_.size(simplex)*12 + gridView_.size(pyramid)*16 + gridView_.size(prism)*18 + gridView_.size(cube)*24;
}
private:
const Problem& problem_() const
{ return *problemPtr_; }
const Problem* problemPtr_;
GridView gridView_;
Dumux::ElementMap<GridView> elementMap_;
std::vector<SubControlVolume> scvs_;
std::vector<SubControlVolumeFace> scvfs_;
std::vector<std::vector<IndexType>> scvfIndicesOfScv_;
std::vector<bool> fpsVertices_;
std::vector<bool> boundaryVertices_;
IndexType numBoundaryScvf_;
GlobalInteractionVolumeSeeds globalInteractionVolumeSeeds_;
};
} // end namespace
#endif
// -*- 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 the global interaction volumes of the mpfa-o method.
*/
#ifndef DUMUX_DISCRETIZATION_MPFA_O_HYBRIDFPS_GLOBALINTERACTIONVOLUMESEEDS_HH
#define DUMUX_DISCRETIZATION_MPFA_O_HYBRIDFPS_GLOBALINTERACTIONVOLUMESEEDS_HH
#include <dumux/discretization/cellcentered/mpfa/methods.hh>
#include <dumux/discretization/cellcentered/mpfa/facetypes.hh>
namespace Dumux
{
/*!
* \ingroup MpfaO
* \brief Base class for the creation and storage of the interaction volume seeds for mpfa methods.
*/
template<class TypeTag>
class CCMpfaOHybridFpsGlobalInteractionVolumeSeeds
{
using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
using Helper = typename GET_PROP_TYPE(TypeTag, MpfaHelper);
using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
using InteractionVolumeOMethod = typename GET_PROP_TYPE(TypeTag, InteractionVolume);
using OSeed = typename InteractionVolumeOMethod::Seed;
using InteractionVolumeFpsMethod = typename GET_PROP_TYPE(TypeTag, BoundaryInteractionVolume);
using FpsSeed = typename InteractionVolumeFpsMethod::Seed;
using Element = typename GridView::template Codim<0>::Entity;
using IndexType = typename GridView::IndexSet::IndexType;
using LocalIndexType = typename InteractionVolumeOMethod::LocalIndexType;
static const int dim = GridView::dimension;
static const int numEq = GET_PROP_VALUE(TypeTag, NumEq);
public:
CCMpfaOHybridFpsGlobalInteractionVolumeSeeds(const GridView gridView) : gridView_(gridView) {}
// initializes the interaction volumes or the seeds
void update(const Problem& problem, const std::vector<bool>& fpsVertices, const std::vector<bool>& boundaryVertices)
{
problemPtr_ = &problem;
fpsSeeds_.clear();
oSeeds_.clear();
// -1 indicates that the scvf has not been handled yet
auto numScvf = problem_().model().fvGridGeometry().numScvf();
scvfIndexMap_.resize(numScvf, -1);
// detect and handle the boundary first
initializeFpsSeeds_(fpsVertices, boundaryVertices);
initializeOSeeds_();
}
const OSeed& seed(const SubControlVolumeFace& scvf) const
{ return oSeeds_[scvfIndexMap_[scvf.index()]]; }
const FpsSeed& boundarySeed(const SubControlVolumeFace& scvf, const LocalIndexType eqIdx) const
{ return fpsSeeds_[scvfIndexMap_[scvf.index()]]; }
private:
void initializeFpsSeeds_(const std::vector<bool>& fpsVertices, const std::vector<bool>& boundaryVertices)
{
fpsSeeds_.reserve(fpsVertices.size());
IndexType fpsSeedIndex = 0;
for (const auto& element : elements(gridView_))
{
auto fvGeometry = localView(problem_().model().fvGridGeometry());
fvGeometry.bind(element);
for (const auto& scvf : scvfs(fvGeometry))
{
// skip the rest if we already handled this face or if face doesn't belong to an fps element
if (scvfIndexMap_[scvf.index()] != -1 || !fpsVertices[scvf.vertexIndex()])
continue;
// the fps interaction volume seed
auto fpsSeed = Helper::makeBoundaryInteractionVolumeSeed(problem_(), element, fvGeometry, scvf, boundaryVertices[scvf.vertexIndex()]);
// update the index map entries for the global scv faces in the interaction volume
for (const auto& localScvfSeed : fpsSeed.scvfSeeds())
for (const auto scvfIdxGlobal : localScvfSeed.globalScvfIndices())
scvfIndexMap_[scvfIdxGlobal] = fpsSeedIndex;
// store interaction volume and increment counter
fpsSeeds_.emplace_back(std::move(fpsSeed));
fpsSeedIndex++;
}
}
// shrink fps seed vector to actual size
fpsSeeds_.shrink_to_fit();
}
void initializeOSeeds_()
{
oSeeds_.reserve(problem_().gridView().size(dim));
IndexType oSeedIndex = 0;
for (const auto& element : elements(gridView_))
{
auto fvGeometry = localView(problem_().model().fvGridGeometry());
fvGeometry.bind(element);
for (const auto& scvf : scvfs(fvGeometry))
{
if (scvfIndexMap_[scvf.index()] != -1)
continue;
// the inner interaction volume seed
auto seed = 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] = oSeedIndex;
// store interaction volume and increment counter
oSeeds_.emplace_back(std::move(seed));
oSeedIndex++;
}
}
// shrink seed vector to actual size
oSeeds_.shrink_to_fit();
}
const Problem& problem_() const
{ return *problemPtr_; }
const Problem* problemPtr_;
GridView gridView_;
std::vector<IndexType> scvfIndexMap_;
std::vector<FpsSeed> fpsSeeds_;
std::vector<OSeed> oSeeds_;
};
} // end namespace
#endif
// -*- 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 Helper class to get the required information on an interaction volume.
*/
#ifndef DUMUX_DISCRETIZATION_CC_MPFAO_HYBRID_FPS_HELPER_HH
#define DUMUX_DISCRETIZATION_CC_MPFAO_HYBRID_FPS_HELPER_HH
#include <dumux/discretization/cellcentered/mpfa/omethod/helper.hh>
namespace Dumux
{
/*!
* \ingroup Mpfa
* \brief Helper class to get the required information on an interaction volume.
* Specialization for the Mpfa-O method in two dimensions.
*/
template<class TypeTag, int dim>
class MpfaHelperBase<TypeTag, MpfaMethods::oMethodFpsHybrid, dim> : public MpfaHelperBase<TypeTag, MpfaMethods::oMethod, dim>
{
using Base = MpfaHelperBase<TypeTag, MpfaMethods::oMethod, dim>;
using Implementation = typename GET_PROP_TYPE(TypeTag, MpfaHelper);
using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
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 GlobalIndexType = typename Base::GlobalIndexType;
using LocalIndexType = typename Base::LocalIndexType;
using GlobalIndexSet = typename Base::GlobalIndexSet;
using LocalIndexSet = typename Base::LocalIndexSet;
using InteractionVolumeSeed = typename InteractionVolume::Seed;
using ScvSeed = typename InteractionVolumeSeed::LocalScvSeed;
using ScvfSeed = typename InteractionVolumeSeed::LocalScvfSeed;
using Element = typename GridView::template Codim<0>::Entity;
public:
static InteractionVolumeSeed makeInnerInteractionVolumeSeed(const Problem& problem,
const Element& element,
const FVElementGeometry& fvGeometry,
const SubControlVolumeFace& scvf)
{
std::vector<ScvSeed> scvSeeds;
std::vector<ScvfSeed> scvfSeeds;
fillEntitySeeds_(scvSeeds, scvfSeeds, problem, element, fvGeometry, scvf, /*dummy*/0);
return InteractionVolumeSeed(std::move(scvSeeds), std::move(scvfSeeds), false);
}
static InteractionVolumeSeed makeBoundaryInteractionVolumeSeed(const Problem& problem,
const Element& element,
const FVElementGeometry& fvGeometry,
const SubControlVolumeFace& scvf,
const bool isBoundary)
{
std::vector<ScvSeed> scvSeeds;
std::vector<ScvfSeed> scvfSeeds;
fillEntitySeeds_(scvSeeds, scvfSeeds, problem, element, fvGeometry, scvf, isBoundary);
return InteractionVolumeSeed(std::move(scvSeeds), std::move(scvfSeeds), isBoundary);
}
private:
static void fillEntitySeeds_(std::vector<ScvSeed>& scvSeeds,
std::vector<ScvfSeed>& scvfSeeds,
const Problem& problem,
const Element& element,
const FVElementGeometry& fvGeometry,
const SubControlVolumeFace& scvf,
const bool isBoundary)
{
// Get the two scv faces in the first scv
auto scvfVector = Implementation::getScvFacesAtVertex(scvf.vertexIndex(), element, fvGeometry);
// The global index of the first scv of the interaction region
auto scvIdx0 = scvf.insideScvIdx();
// rotate counter clockwise and create the entities
Implementation::performRotation_(problem, scvfVector, scvSeeds, scvfSeeds, scvIdx0, /*dummy*/0);
if (isBoundary)
{
// the local scvf index of the second local face of the first local scv
LocalIndexType storeIdx = scvfSeeds.size();
// clockwise rotation until hitting the boundary again
Implementation::performRotation_(problem, scvfVector, scvSeeds, scvfSeeds, scvIdx0, /*dummy*/0, /*clockwise*/true);
// Finish by creating the first scv
scvSeeds.emplace(scvSeeds.begin(), ScvSeed(GlobalIndexSet({scvfVector[0]->index(), scvfVector[1]->index()}),
LocalIndexSet({0, storeIdx}),
scvIdx0));
}
else
// Finish by creating the first scv
scvSeeds.emplace(scvSeeds.begin(), ScvSeed(GlobalIndexSet({scvfVector[0]->index(), scvfVector[1]->index()}),
LocalIndexSet({0, static_cast<LocalIndexType>(scvfSeeds.size()-1)}),
scvIdx0));
}
};
} // end namespace
#endif
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*****************************************************************************