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

[tpfa][stencils] introduce non-symmetric stencil class

The stencils vector class is now set separately for the specific cell-centered
methods. Some methods might lead to non-symmetric global matrices, thus, the stencils
have to be set up differently.
parent 3ed76500
......@@ -23,17 +23,16 @@
#ifndef DUMUX_DISCRETIZATION_CC_STENCILS_HH
#define DUMUX_DISCRETIZATION_CC_STENCILS_HH
#include <set>
#include <dumux/implicit/cellcentered/properties.hh>
namespace Dumux
{
/*!
* \brief Element-related stencils
* \brief Element-related stencils for symmetric cc methods.
*/
template<class TypeTag>
class CCElementStencils
class CCSymmetricElementStencils
{
using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
......@@ -44,7 +43,7 @@ class CCElementStencils
// TODO a separate stencil class all stencils can derive from?
using Stencil = std::vector<IndexType>;
public:
//! Update the stencil. We expect a bound fvGeometry
//! Update the stencil.
void update(const Problem& problem,
const Element& element)
{
......@@ -93,14 +92,14 @@ private:
/*!
* \ingroup CCModel
* \brief The global stencil container class
* \brief The global stencil container class for symmetric cc methods
*/
template<class TypeTag>
class CCStencilsVector
class CCSymmetricStencilsVector
{
using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
using StencilType = CCElementStencils<TypeTag>;
using StencilType = CCSymmetricElementStencils<TypeTag>;
public:
void update(Problem& problem)
......@@ -123,7 +122,135 @@ public:
}
private:
std::vector<CCElementStencils<TypeTag>> elementStencils_;
std::vector<StencilType> elementStencils_;
const Problem* problemPtr_;
};
//! Forward declaration of the global stencil class for nonsymmetric cc methods
template<class TypeTag> class CCNonSymmetricStencilsVector;
/*!
* \brief Element-related stencils for non-symmetric cc methods
*/
template<class TypeTag>
class CCNonSymmetricElementStencils
{
using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
using FluxVariables = typename GET_PROP_TYPE(TypeTag, FluxVariables);
using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
using IndexType = typename GridView::IndexSet::IndexType;
using Element = typename GridView::template Codim<0>::Entity;
// TODO a separate stencil class all stencils can derive from?
using Stencil = std::vector<IndexType>;
public:
//! Update the stencil. We expect a bound fvGeometry
void update(const Problem& problem,
const Element& element,
CCNonSymmetricStencilsVector<TypeTag>& stencilsVector)
{
// the element index of our own here
auto eIdx = problem.elementMapper().index(element);
elementStencil_.clear();
// restrict the FvGeometry locally and bind to the element
auto fvGeometry = localView(problem.model().globalFvGeometry());
fvGeometry.bindElement(element);
// loop over sub control faces
for (auto&& scvf : scvfs(fvGeometry))
{
FluxVariables fluxVars;
const auto& stencil = fluxVars.computeStencil(problem, element, fvGeometry, scvf);
// insert stencil into the element stencil
elementStencil_.insert(elementStencil_.end(), stencil.begin(), stencil.end());
// insert our index in the neighbor stencils of the elements in the flux stencil
for (auto globalI : stencil)
{
if (globalI != eIdx)
stencilsVector[globalI].addNeighbor(eIdx);
}
}
// make values in element stencil unique
std::sort(elementStencil_.begin(), elementStencil_.end());
elementStencil_.erase(std::unique(elementStencil_.begin(), elementStencil_.end()), elementStencil_.end());
}
//! The full element stencil (all element this element is interacting with)
const Stencil& elementStencil() const
{
return elementStencil_;
}
//! The full element stencil without this element
const Stencil& neighborStencil() const
{
return neighborStencil_;
}
void addNeighbor(const IndexType nIdx)
{
neighborStencil_.push_back(nIdx);
}
void makeNeighborStencilUnique()
{
// make values in neighbor stencil unique
std::sort(neighborStencil_.begin(), neighborStencil_.end());
neighborStencil_.erase(std::unique(neighborStencil_.begin(), neighborStencil_.end()), neighborStencil_.end());
}
private:
Stencil elementStencil_;
Stencil neighborStencil_;
};
/*!
* \ingroup CCModel
* \brief The global stencil container class for non-symmetric cc methods
*/
template<class TypeTag>
class CCNonSymmetricStencilsVector
{
using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
using IndexType = typename GridView::IndexSet::IndexType;
using StencilType = CCNonSymmetricElementStencils<TypeTag>;
// The stencil type of an element has to access bracket operator
friend StencilType;
public:
void update(Problem& problem)
{
problemPtr_ = &problem;
elementStencils_.resize(problem.gridView().size(0));
for (const auto& element : elements(problem.gridView()))
{
auto eIdx = problem.elementMapper().index(element);
elementStencils_[eIdx].update(problem, element, *this);
}
for (auto&& stencil : elementStencils_)
stencil.makeNeighborStencilUnique();
}
//! overload for elements
template <class Entity>
typename std::enable_if<Entity::codimension == 0, const StencilType&>::type
get(const Entity& entity) const
{
return elementStencils_[problemPtr_->elementMapper().index(entity)];
}
private:
StencilType& operator[] (const IndexType globalIdx)
{ return elementStencils_[globalIdx]; }
std::vector<StencilType> elementStencils_;
const Problem* problemPtr_;
};
......
......@@ -30,7 +30,6 @@
#include <dumux/implicit/propertydefaults.hh>
#include <dumux/discretization/cellcentered/globalvolumevariables.hh>
#include <dumux/discretization/cellcentered/subcontrolvolume.hh>
#include <dumux/discretization/cellcentered/stencils.hh>
#include "elementboundarytypes.hh"
#include "localresidual.hh"
......@@ -44,7 +43,6 @@ namespace Dumux
// forward declaration
template<class TypeTag> class CCElementBoundaryTypes;
template<class TypeTag> class CCLocalResidual;
template<class TypeTag> class CCStencilsVector;
namespace Properties
{
......@@ -55,13 +53,10 @@ SET_TYPE_PROP(CCModel, ElementBoundaryTypes, CCElementBoundaryTypes<TypeTag>);
SET_TYPE_PROP(CCModel, DofMapper, typename GET_PROP_TYPE(TypeTag, ElementMapper));
//! The local jacobian operator
SET_TYPE_PROP(CCModel, LocalJacobian, Dumux::CCLocalJacobian<TypeTag>);
SET_TYPE_PROP(CCModel, LocalJacobian, CCLocalJacobian<TypeTag>);
//! Assembler for the global jacobian matrix
SET_TYPE_PROP(CCModel, JacobianAssembler, Dumux::CCAssembler<TypeTag>);
//! The stencil container
SET_TYPE_PROP(CCModel, StencilsVector, CCStencilsVector<TypeTag>);
SET_TYPE_PROP(CCModel, JacobianAssembler, CCAssembler<TypeTag>);
//! The sub control volume
SET_PROP(CCModel, SubControlVolume)
......@@ -71,11 +66,11 @@ private:
using ScvGeometry = typename Grid::template Codim<0>::Geometry;
using IndexType = typename Grid::LeafGridView::IndexSet::IndexType;
public:
typedef Dumux::CCSubControlVolume<ScvGeometry, IndexType> type;
using type = CCSubControlVolume<ScvGeometry, IndexType>;
};
//! The global current volume variables vector class
SET_TYPE_PROP(CCModel, GlobalVolumeVariables, Dumux::CCGlobalVolumeVariables<TypeTag, GET_PROP_VALUE(TypeTag, EnableGlobalVolumeVariablesCache)>);
SET_TYPE_PROP(CCModel, GlobalVolumeVariables, CCGlobalVolumeVariables<TypeTag, GET_PROP_VALUE(TypeTag, EnableGlobalVolumeVariablesCache)>);
//! Set the BaseLocalResidual to CCLocalResidual
SET_TYPE_PROP(CCModel, BaseLocalResidual, CCLocalResidual<TypeTag>);
......
......@@ -35,6 +35,7 @@
#include <dumux/discretization/cellcentered/tpfa/elementvolumevariables.hh>
#include <dumux/discretization/cellcentered/tpfa/elementfluxvariablescache.hh>
#include <dumux/discretization/cellcentered/tpfa/subcontrolvolumeface.hh>
#include <dumux/discretization/cellcentered/stencils.hh>
#include <dumux/implicit/cellcentered/properties.hh>
#include <dumux/discretization/methods.hh>
......@@ -60,10 +61,13 @@ SET_TYPE_PROP(CCTpfaModel, GlobalFluxVariablesCache, Dumux::CCTpfaGlobalFluxVari
SET_TYPE_PROP(CCTpfaModel, FVElementGeometry, CCTpfaFVElementGeometry<TypeTag, GET_PROP_VALUE(TypeTag, EnableGlobalFVGeometryCache)>);
//! The global previous volume variables vector class
SET_TYPE_PROP(CCModel, ElementVolumeVariables, Dumux::CCTpfaElementVolumeVariables<TypeTag, GET_PROP_VALUE(TypeTag, EnableGlobalVolumeVariablesCache)>);
SET_TYPE_PROP(CCTpfaModel, ElementVolumeVariables, Dumux::CCTpfaElementVolumeVariables<TypeTag, GET_PROP_VALUE(TypeTag, EnableGlobalVolumeVariablesCache)>);
//! The local flux variables cache vector class
SET_TYPE_PROP(CCModel, ElementFluxVariablesCache, Dumux::CCTpfaElementFluxVariablesCache<TypeTag, GET_PROP_VALUE(TypeTag, EnableGlobalFluxVariablesCache)>);
SET_TYPE_PROP(CCTpfaModel, ElementFluxVariablesCache, Dumux::CCTpfaElementFluxVariablesCache<TypeTag, GET_PROP_VALUE(TypeTag, EnableGlobalFluxVariablesCache)>);
//! The stencil container, tpfa leads to a symmetric global matrix
SET_TYPE_PROP(CCTpfaModel, StencilsVector, Dumux::CCSymmetricStencilsVector<TypeTag>);
SET_PROP(CCTpfaModel, SubControlVolumeFace)
{
......
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