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

Merge branch 'feature/gravity-mpfa' into 'next'

Feature/gravity mpfa

See merge request !710
parents 0b3198f9 26dfeee8
add_subdirectory("omethod")
install(FILES
computetransmissibility
connectivitymap.hh
darcyslaw.hh
dualgridindexset.hh
......@@ -18,6 +19,8 @@ helper.hh
interactionvolumebase.hh
interactionvolumedatahandle.hh
interactionvolume.hh
localassembler.hh
localfacedata.hh
methods.hh
properties.hh
subcontrolvolumeface.hh
......
// -*- 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
* \ingroup CCMpfaDiscretization
* \brief This file contains free functions to evaluate the transmissibilities
* associated with flux evaluations across sub-control volume faces
* in the context of cell-centered Mpfa schemes.
*/
#ifndef DUMUX_DISCRETIZATION_CC_MPFA_COMPUTE_TRANSMISSIBILITY_HH
#define DUMUX_DISCRETIZATION_CC_MPFA_COMPUTE_TRANSMISSIBILITY_HH
#include <dune/common/typetraits.hh>
#include <dune/common/fvector.hh>
#include <dumux/common/math.hh>
namespace Dumux
{
/*!
* \ingroup CCMpfaDiscretization
* \brief Free function to evaluate the Mpfa transmissibility associated
* with the flux (in the form of flux = T*gradU) across a
* sub-control volume face stemming from a given sub-control
* volume with corresponding tensor T.
*
* \param scv The iv-local sub-control volume
* \param scvf The grid sub-control volume face
* \param T The tensor living in the scv
* \param extrusionFactor The extrusion factor of the scv
*/
template<class Scv, class Scvf, class Tensor>
Dune::FieldVector< typename Tensor::field_type, Scv::myDimension >
computeMpfaTransmissibility(const Scv& scv,
const Scvf& scvf,
const Tensor& T,
typename Scv::ctype extrusionFactor)
{
Dune::FieldVector< typename Tensor::field_type, Scv::myDimension > wijk;
for (unsigned int dir = 0; dir < Scv::myDimension; ++dir)
wijk[dir] = vtmv(scvf.unitOuterNormal(), T, scv.nu(dir));
wijk *= scvf.area()*extrusionFactor;
wijk /= scv.detX();
return wijk;
}
/*!
* \ingroup CCMpfaDiscretization
* \brief Free function to evaluate the Mpfa transmissibility associated
* with the flux (in the form of flux = T*gradU) across a
* sub-control volume face stemming from a given sub-control
* volume with corresponding tensor T, where T is a scalar.
*
* \param scv The iv-local sub-control volume
* \param scvf The grid sub-control volume face
* \param t the scalar quantity living in the scv
* \param extrusionFactor The extrusion factor of the scv
*/
template< class Scv,
class Scvf,
class Tensor,
std::enable_if_t< Dune::IsNumber<Tensor>::value, int > = 1 >
Dune::FieldVector<Tensor, Scv::myDimension>
computeMpfaTransmissibility(const Scv& scv,
const Scvf& scvf,
Tensor t,
typename Scv::ctype extrusionFactor)
{
Dune::FieldVector<Tensor, Scv::myDimension> wijk;
for (unsigned int dir = 0; dir < Scv::myDimension; ++dir)
wijk[dir] = vtmv(scvf.unitOuterNormal(), t, scv.nu(dir));
wijk *= scvf.area()*extrusionFactor;
wijk /= scv.detX();
return wijk;
}
} // end namespace Dumux
#endif
......@@ -18,8 +18,11 @@
*****************************************************************************/
/*!
* \file
* \ingroup CCMpfaDiscretization
* \brief Stores the face indices corresponding to the neighbors of an element
* that contribute to the derivative calculation
* that contribute to the derivative calculation. Depending on if an
* mpfa scheme leads to a symmetric/unsymmetric sparsity pattern, the
* adequate implementation of the connectiviy map is chosen.
*/
#ifndef DUMUX_CC_MPFA_CONNECTIVITY_MAP_HH
#define DUMUX_CC_MPFA_CONNECTIVITY_MAP_HH
......@@ -35,15 +38,15 @@ namespace Dumux
template<class TypeTag, MpfaMethods method>
class CCMpfaConnectivityMapImplementation;
// //! The Assembly map for models using mpfa methods
//! The Assembly map for models using mpfa methods
template<class TypeTag>
using CCMpfaConnectivityMap = CCMpfaConnectivityMapImplementation<TypeTag, GET_PROP_VALUE(TypeTag, MpfaMethod)>;
//! The default is the general assembly map for mpfa schemes
//! The default is the general assembly map for mpfa schemes (non-symmetric schemes)
template<class TypeTag, MpfaMethods method>
class CCMpfaConnectivityMapImplementation : public CCMpfaGeneralConnectivityMap<TypeTag> {};
//! The o-method can use the simple assembly map
//! The o-method can use the simple (symmetric) assembly map
template<class TypeTag>
class CCMpfaConnectivityMapImplementation<TypeTag, MpfaMethods::oMethod> : public CCSimpleConnectivityMap<TypeTag> {};
}
......
......@@ -18,27 +18,31 @@
*****************************************************************************/
/*!
* \file
* \brief This file contains the class which is required to calculate
* volume and mass fluxes of fluid phases over a face of a finite volume by means
* of the Darcy approximation. This specializations is for cell-centered schemes
* using multi-point flux approximation.
* \ingroup CCMpfaDiscretization
* \brief Darcy's Law for cell-centered finite volume schemes
* with multi-point flux approximation.
*/
#ifndef DUMUX_DISCRETIZATION_CC_MPFA_DARCYS_LAW_HH
#define DUMUX_DISCRETIZATION_CC_MPFA_DARCYS_LAW_HH
#include <dune/common/densevector.hh>
#include <dune/common/densematrix.hh>
#include <dumux/common/properties.hh>
#include <dumux/common/parameters.hh>
#include <dumux/discretization/methods.hh>
#include <dumux/discretization/methods.hh>
namespace Dumux
{
// forward declarations
//! forward declaration of the method-specific implementation
template<class TypeTag, DiscretizationMethods discMethod>
class DarcysLawImplementation;
/*!
* \ingroup Mpfa
* \brief Specialization of Darcy's Law for the CCMpfa method.
* \ingroup CCMpfaDiscretization
* \brief Darcy's law for cell-centered finite volume schemes with multi-point flux approximation.
*/
template <class TypeTag>
class DarcysLawImplementation<TypeTag, DiscretizationMethods::CCMpfa>
......@@ -47,21 +51,13 @@ class DarcysLawImplementation<TypeTag, DiscretizationMethods::CCMpfa>
using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
using Element = typename GridView::template Codim<0>::Entity;
using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);
using ElementFluxVariablesCache = typename GET_PROP_TYPE(TypeTag, ElementFluxVariablesCache);
using FluxVariablesCache = typename GET_PROP_TYPE(TypeTag, FluxVariablesCache);
// Always use the dynamic type for vectors (compatibility with the boundary)
using PrimaryInteractionVolume = typename GET_PROP_TYPE(TypeTag, PrimaryInteractionVolume);
using CoefficientVector = typename PrimaryInteractionVolume::Traits::DynamicVector;
using DataHandle = typename PrimaryInteractionVolume::Traits::DataHandle;
static constexpr int numPhases = GET_PROP_VALUE(TypeTag, NumPhases);
static constexpr int dim = GridView::dimension;
static constexpr int dimWorld = GridView::dimensionworld;
//! Class that fills the cache corresponding to mpfa Darcy's Law
class MpfaDarcysLawCacheFiller
{
......@@ -77,14 +73,16 @@ class DarcysLawImplementation<TypeTag, DiscretizationMethods::CCMpfa>
const SubControlVolumeFace& scvf,
const FluxVariablesCacheFiller& fluxVarsCacheFiller)
{
// get interaction volume from the flux vars cache filler & upate the cache
// get interaction volume related data from the filler class & upate the cache
if (fvGeometry.fvGridGeometry().vertexUsesSecondaryInteractionVolume(scvf.vertexIndex()))
scvfFluxVarsCache.updateAdvection(fluxVarsCacheFiller.secondaryInteractionVolume(),
fluxVarsCacheFiller.dataHandle(),
fluxVarsCacheFiller.secondaryIvLocalFaceData(),
fluxVarsCacheFiller.secondaryIvDataHandle(),
scvf);
else
scvfFluxVarsCache.updateAdvection(fluxVarsCacheFiller.primaryInteractionVolume(),
fluxVarsCacheFiller.dataHandle(),
fluxVarsCacheFiller.primaryIvLocalFaceData(),
fluxVarsCacheFiller.primaryIvDataHandle(),
scvf);
}
};
......@@ -92,54 +90,137 @@ class DarcysLawImplementation<TypeTag, DiscretizationMethods::CCMpfa>
//! The cache used in conjunction with the mpfa Darcy's Law
class MpfaDarcysLawCache
{
// We always use the dynamic types here to be compatible on the boundary
using Stencil = typename PrimaryInteractionVolume::Traits::DynamicGlobalIndexContainer;
using DirichletDataContainer = typename PrimaryInteractionVolume::DirichletDataContainer;
static constexpr int dim = GridView::dimension;
static constexpr int dimWorld = GridView::dimensionworld;
static constexpr int numPhases = GET_PROP_VALUE(TypeTag, NumPhases);
using GridIndexType = typename GridView::IndexSet::IndexType;
// In the current implementation of the flux variables cache we cannot
// make a disctinction between dynamic (mpfa-o) and static (mpfa-l)
// matrix and vector types, as currently the cache class can only be templated
// by a type tag (and there can only be one). We use a dynamic vector here to
// make sure it works in case one of the two used interaction volume types uses
// dynamic types performance is thus lowered for schemes using static types.
// TODO: this has to be checked thoroughly as soon as a scheme using static types
// is implemented. One idea to overcome the performance drop could be only
// storing the iv-local index here and obtain tij always from the datahandle
// of the fluxVarsCacheContainer
using Vector = Dune::DynamicVector< Scalar >;
using Matrix = Dune::DynamicMatrix< Scalar >;
using Stencil = std::vector< GridIndexType >;
using PrimaryInteractionVolume = typename GET_PROP_TYPE(TypeTag, PrimaryInteractionVolume);
using PrimaryIvVector = typename PrimaryInteractionVolume::Traits::Vector;
using PrimaryIvMatrix = typename PrimaryInteractionVolume::Traits::Matrix;
using PrimaryStencil = typename PrimaryInteractionVolume::Traits::Stencil;
static_assert( std::is_convertible<PrimaryIvVector*, Vector*>::value,
"The vector type used in primary interaction volumes is not convertible to Dune::DynamicVector!" );
static_assert( std::is_convertible<PrimaryIvMatrix*, Matrix*>::value,
"The matrix type used in primary interaction volumes is not convertible to Dune::DynamicMatrix!" );
static_assert( std::is_convertible<PrimaryStencil*, Stencil*>::value,
"The stencil type used in primary interaction volumes is not convertible to std::vector<GridIndexType>!" );
using SecondaryInteractionVolume = typename GET_PROP_TYPE(TypeTag, SecondaryInteractionVolume);
using SecondaryIvVector = typename SecondaryInteractionVolume::Traits::Vector;
using SecondaryIvMatrix = typename SecondaryInteractionVolume::Traits::Matrix;
using SecondaryStencil = typename SecondaryInteractionVolume::Traits::Stencil;
static_assert( std::is_convertible<SecondaryIvVector*, Vector*>::value,
"The vector type used in secondary interaction volumes is not convertible to Dune::DynamicVector!" );
static_assert( std::is_convertible<SecondaryIvMatrix*, Matrix*>::value,
"The matrix type used in secondary interaction volumes is not convertible to Dune::DynamicMatrix!" );
static_assert( std::is_convertible<SecondaryStencil*, Stencil*>::value,
"The stencil type used in secondary interaction volumes is not convertible to std::vector<GridIndexType>!" );
public:
// export the filler type
using Filler = MpfaDarcysLawCacheFiller;
//! update cached objects
template<class InteractionVolume>
void updateAdvection(const InteractionVolume& iv, const DataHandle& dataHandle, const SubControlVolumeFace &scvf)
/*!
* \brief Update cached objects (transmissibilities and gravity)
*
* \tparam InteractionVolume The (mpfa scheme-specific) interaction volume
* \tparam LocalFaceData The object used to store iv-local info on an scvf
* \tparam DataHandle The object used to store transmissibility matrices etc.
*
* \param iv The interaction volume this scvf is embedded in
* \param localFaceData iv-local info on this scvf
* \param dataHandle Transmissibility matrix & gravity data of this iv
* \param scvf The sub-control volume face
*/
template< class InteractionVolume, class LocalFaceData, class DataHandle >
void updateAdvection(const InteractionVolume& iv,
const LocalFaceData& localFaceData,
const DataHandle& dataHandle,
const SubControlVolumeFace &scvf)
{
const auto& localFaceData = iv.getLocalFaceData(scvf);
stencil_ = &iv.stencil();
switchFluxSign_ = localFaceData.isOutside();
for (unsigned int pIdx = 0; pIdx < numPhases; ++pIdx)
pj_[pIdx] = &dataHandle.pressures(pIdx);
// update the quantities that are equal for all phases
advectionSwitchFluxSign_ = localFaceData.isOutside();
advectionVolVarsStencil_ = &dataHandle.volVarsStencil();
advectionDirichletData_ = &dataHandle.dirichletData();
static const bool enableGravity = getParamFromGroup<bool>(GET_PROP_VALUE(TypeTag, ModelParameterGroup), "Problem.EnableGravity");
const auto ivLocalIdx = localFaceData.ivLocalScvfIndex();
// the transmissibilities on surface grids have to be obtained from the outside
// standard grids
if (dim == dimWorld)
advectionTij_ = &dataHandle.T()[localFaceData.ivLocalScvfIndex()];
{
Tij_ = &dataHandle.advectionT()[ivLocalIdx];
if (enableGravity)
for (unsigned int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
g_[phaseIdx] = dataHandle.gravity(phaseIdx)[ivLocalIdx];
}
// surface grids
else
advectionTij_ = localFaceData.isOutside() ?
&dataHandle.outsideTij()[localFaceData.ivLocalOutsideScvfIndex()] :
&dataHandle.T()[localFaceData.ivLocalScvfIndex()];
{
if (!localFaceData.isOutside())
{
Tij_ = &dataHandle.advectionT()[ivLocalIdx];
if (enableGravity)
for (unsigned int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
g_[phaseIdx] = dataHandle.gravity(phaseIdx)[ivLocalIdx];
}
else
{
const auto idxInOutsideFaces = localFaceData.scvfLocalOutsideScvfIndex();
Tij_ = &dataHandle.advectionTout()[ivLocalIdx][idxInOutsideFaces];
if (enableGravity)
for (unsigned int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
g_[phaseIdx] = dataHandle.gravityOutside(phaseIdx)[ivLocalIdx][idxInOutsideFaces];
}
}
}
//! Returns the stencil for advective scvf flux computation
const Stencil& advectionVolVarsStencil() const { return *advectionVolVarsStencil_; }
//! The stencil corresponding to the transmissibilities
const Stencil& advectionStencil() const { return *stencil_; }
//! Coefficients for the cell (& Dirichlet) unknowns in flux expressions
const Vector& advectionTij() const { return *Tij_; }
//! Returns the transmissibilities associated with the volume variables
//! All phases flow through the same rock, thus, tij are equal for all phases
const CoefficientVector& advectionTij() const { return *advectionTij_; }
//! The cell (& Dirichlet) pressures within this interaction volume
const Vector& pressures(unsigned int phaseIdx) const { return *pj_[phaseIdx]; }
//! On faces that are "outside" w.r.t. a face in the interaction volume,
//! we have to take the negative value of the fluxes, i.e. multiply by -1.0
bool advectionSwitchFluxSign() const { return advectionSwitchFluxSign_; }
//! The gravitational acceleration for a phase on this scvf
Scalar gravity(unsigned int phaseIdx) const { return g_[phaseIdx]; }
//! Returns the data on dirichlet boundary conditions affecting
//! the flux computation on this face
const DirichletDataContainer& advectionDirichletData() const { return *advectionDirichletData_; }
//! In the interaction volume-local system of eq we have one unknown per face.
//! On scvfs on this face, but in "outside" (neighbor) elements of it, we have
//! to take the negative value of the fluxes due to the flipped normal vector.
//! This function returns whether or not this scvf is an "outside" face in the iv.
bool advectionSwitchFluxSign() const { return switchFluxSign_; }
private:
bool advectionSwitchFluxSign_;
const Stencil* advectionVolVarsStencil_;
const CoefficientVector* advectionTij_;
const DirichletDataContainer* advectionDirichletData_;
bool switchFluxSign_;
const Stencil* stencil_; //!< The stencil, i.e. the grid indices j
const Vector* Tij_; //!< The transmissibilities such that f = Tij*pj
std::array< Scalar, numPhases > g_; //!< Gravitational flux contribution on this face
std::array<const Vector*, numPhases> pj_; //!< The interaction-volume wide phase pressures pj
};
public:
......@@ -149,6 +230,7 @@ public:
// export the type for the corresponding cache
using Cache = MpfaDarcysLawCache;
//! Compute the advective flux across an scvf
static Scalar flux(const Problem& problem,
const Element& element,
const FVElementGeometry& fvGeometry,
......@@ -157,78 +239,23 @@ public:
const unsigned int phaseIdx,
const ElementFluxVariablesCache& elemFluxVarsCache)
{
static const bool gravity = getParamFromGroup<bool>(GET_PROP_VALUE(TypeTag, ModelParameterGroup), "Problem.EnableGravity");
// Calculate the interface density for gravity evaluation
const auto rho = interpolateDensity(elemVolVars, scvf, phaseIdx);
// prepare computations
unsigned int i = 0;
Scalar scvfFlux(0.0);
const auto& fluxVarsCache = elemFluxVarsCache[scvf];
const auto& tij = fluxVarsCache.advectionTij();
const auto& pj = fluxVarsCache.pressures(phaseIdx);
// add contributions from cell-centered unknowns
for (const auto volVarIdx : fluxVarsCache.advectionVolVarsStencil())
{
const auto& volVars = elemVolVars[volVarIdx];
Scalar h = volVars.pressure(phaseIdx);
// if gravity is enabled, add gravitational acceleration
if (gravity)
{
// gravitational acceleration in the center of the actual element
const auto x = fvGeometry.scv(volVarIdx).center();
const auto g = problem.gravityAtPos(x);
h -= rho*(g*x);
}
scvfFlux += tij[i++]*h;
}
// add contributions from possible dirichlet boundary conditions
for (const auto& d : fluxVarsCache.advectionDirichletData())
{
const auto& volVars = elemVolVars[d.volVarIndex()];
Scalar h = volVars.pressure(phaseIdx);
// maybe add gravitational acceleration
if (gravity)
{
const auto x = d.ipGlobal();
const auto g = problem.gravityAtPos(x);
h -= rho*(g*x);
}
// compute t_ij*p_j
Scalar scvfFlux = tij*pj;
scvfFlux += tij[i++]*h;
}
// maybe add gravitational acceleration
static const bool enableGravity = getParamFromGroup<bool>(GET_PROP_VALUE(TypeTag, ModelParameterGroup), "Problem.EnableGravity");
if (enableGravity)
scvfFlux += fluxVarsCache.gravity(phaseIdx);
// return the flux (maybe adjust the sign)
return fluxVarsCache.advectionSwitchFluxSign() ? -scvfFlux : scvfFlux;
}
// switch the sign if necessary
if (fluxVarsCache.advectionSwitchFluxSign())
scvfFlux *= -1.0;
private:
static Scalar interpolateDensity(const ElementVolumeVariables& elemVolVars,
const SubControlVolumeFace& scvf,
const unsigned int phaseIdx)
{
static const bool gravity = getParamFromGroup<bool>(GET_PROP_VALUE(TypeTag, ModelParameterGroup), "Problem.EnableGravity");
if (!gravity)
return Scalar(0.0);
else
{
// use arithmetic mean of the densities around the scvf
if (!scvf.boundary())
{
Scalar rho = elemVolVars[scvf.insideScvIdx()].density(phaseIdx);
for (const auto outsideIdx : scvf.outsideScvIndices())
rho += elemVolVars[outsideIdx].density(phaseIdx);
return rho/(scvf.outsideScvIndices().size()+1);
}
else
return elemVolVars[scvf.outsideScvIdx()].density(phaseIdx);
}
return scvfFlux;
}
};
......
......@@ -18,39 +18,40 @@
*****************************************************************************/
/*!
* \file
* \brief Base class for the index sets of the dual grid in mpfa schemes.
* \ingroup CCMpfaDiscretization
* \brief Class for the index sets of the dual grid in mpfa schemes.
*/
#ifndef DUMUX_DISCRETIZATION_MPFA_DUALGRID_INDEXSET_BASE_HH
#define DUMUX_DISCRETIZATION_MPFA_DUALGRID_INDEXSET_BASE_HH
#include <dumux/discretization/cellcentered/mpfa/fvgridgeometry.hh>
#ifndef DUMUX_DISCRETIZATION_MPFA_DUALGRID_INDEX_SET_HH
#define DUMUX_DISCRETIZATION_MPFA_DUALGRID_INDEX_SET_HH
namespace Dumux
{
/*!
* \ingroup Mpfa
* \brief Nodal index set for the dual grid of mpfa schemes.
* \ingroup CCMpfaDiscretization
* \brief Nodal index set for mpfa schemes, constructed
* around grid vertices.
*
* \tparam GI The type used for indices on the grid
* \tparam LI The type used for indexing in interaction volumes
* \tparam dim The dimension of the grid
*/
template<class TypeTag>
template< class GI, class LI, int dim >
class DualGridNodalIndexSet
{
using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
using IndexType = typename GridView::IndexSet::IndexType;
using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
using PrimaryInteractionVolume = typename GET_PROP_TYPE(TypeTag, PrimaryInteractionVolume);
using LocalIndexType = typename PrimaryInteractionVolume::Traits::LocalIndexType;
using LocalIndexContainer = typename PrimaryInteractionVolume::Traits::DynamicLocalIndexContainer;
using GlobalIndexContainer = typename PrimaryInteractionVolume::Traits::DynamicGlobalIndexContainer;
public:
using GridIndexType = GI;
using LocalIndexType = LI;
static const int dim = GridView::dimension;
// we use dynamic containers to store indices here
using GridIndexContainer = std::vector<GridIndexType>;
using LocalIndexContainer = std::vector<LocalIndexType>;
public:
//! constructor
//! Constructor
DualGridNodalIndexSet() : numBoundaryScvfs_(0) {}
// Inserts data for a given scvf. This should only called once per scvf!
//! Inserts data for a given scvf
template<typename SubControlVolumeFace>
void insert(const SubControlVolumeFace& scvf)
{
insert(scvf.boundary(),
......@@ -59,20 +60,21 @@ public:
scvf.outsideScvIndices());
}
// Inserts data for a given scvf. This should only called once per scvf!
//! Inserts scvf data
void insert(const bool boundary,
const IndexType scvfIdx,
const IndexType insideScvIdx,
const std::vector<IndexType>& outsideScvIndices)
const GridIndexType scvfIdx,
const GridIndexType insideScvIdx,
const std::vector<GridIndexType>& outsideScvIndices)
{
// this should always be called only once
assert( std::find(scvfIndices_.begin(), scvfIndices_.end(), scvfIdx ) == scvfIndices_.end() && "scvf has already been inserted!");
// this should always be called only once per scvf
assert(std::find(scvfIndices_.begin(), scvfIndices_.end(), scvfIdx ) == scvfIndices_.end()
&& "scvf has already been inserted!");
// the local index of the scvf data about to be inserted
const LocalIndexType curScvfLocalIdx = scvfIndices_.size();
// add global scvf data
GlobalIndexContainer scvIndices;
// add grid scvf data
GridIndexContainer scvIndices;
scvIndices.reserve(outsideScvIndices.size()+1);
scvIndices.push_back(insideScvIdx);