diff --git a/dumux/discretization/cellcentered/mpfa/darcyslaw.hh b/dumux/discretization/cellcentered/mpfa/darcyslaw.hh
index 50bf47b25a56cefb6e836247d1d307ebec9fce20..568c88befa292f696f3eeaf49f5752fd60861638 100644
--- a/dumux/discretization/cellcentered/mpfa/darcyslaw.hh
+++ b/dumux/discretization/cellcentered/mpfa/darcyslaw.hh
@@ -33,6 +33,7 @@
#include <dumux/common/parameters.hh>
#include <dumux/implicit/properties.hh>
+#include <dumux/discretization/cellcentered/mpfa/methods.hh>
namespace Dumux
@@ -53,6 +54,7 @@ class DarcysLawImplementation<TypeTag, DiscretizationMethods::CCMpfa>
{
using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
+ using InteractionVolume = typename GET_PROP_TYPE(TypeTag, InteractionVolume);
using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
@@ -64,6 +66,8 @@ class DarcysLawImplementation<TypeTag, DiscretizationMethods::CCMpfa>
using IndexType = typename GridView::IndexSet::IndexType;
using Stencil = std::vector<IndexType>;
+ static const MpfaMethods method = GET_PROP_VALUE(TypeTag, MpfaMethod);
+
public:
static Scalar flux(const Problem& problem,
@@ -79,19 +83,10 @@ public:
const auto& volVarsPositions = fluxVarsCache.advectionVolVarsPositions(phaseIdx);
const auto& tij = fluxVarsCache.advectionTij(phaseIdx);
- // interface density as arithmetic mean of the two neighbors (when gravity is on)
+ // interface density as arithmetic mean of the neighbors (when gravity is on)
Scalar rho;
if (gravity)
- {
- if (!scvf.boundary())
- {
- rho = elemVolVars[scvf.outsideScvIdx()].density(phaseIdx);
- rho += elemVolVars[scvf.insideScvIdx()].density(phaseIdx);
- rho /= 2.0;
- }
- else
- rho = elemVolVars[scvf.outsideScvIdx()].density(phaseIdx);
- }
+ rho = interpolateDensity(elemVolVars, scvf, phaseIdx);
// calculate Tij*pj
Scalar flux(0.0);
@@ -110,7 +105,6 @@ public:
h -= rho*(g*x);
}
-
flux += tij[localIdx++]*h;
}
@@ -130,6 +124,27 @@ public:
else
return globalFvGeometry.interactionVolumeSeed(scvf).globalScvIndices();
}
+
+private:
+ static Scalar interpolateDensity(const ElementVolumeVariables& elemVolVars,
+ const SubControlVolumeFace& scvf,
+ const unsigned int phaseIdx)
+ {
+ // use arithmetic mean of the densities around the scvf
+ if (!scvf.boundary())
+ {
+ const auto& outsideScvIndices = scvf.outsideScvIndices();
+
+ Scalar rho = elemVolVars[scvf.insideScvIdx()].density(phaseIdx);
+ for (auto outsideIdx : outsideScvIndices)
+ rho += elemVolVars[outsideIdx].density(phaseIdx);
+ rho /= outsideScvIndices.size();
+
+ return rho;
+ }
+ else
+ return elemVolVars[scvf.outsideScvIdx()].density(phaseIdx);
+ }
};
} // end namespace
diff --git a/dumux/discretization/cellcentered/mpfa/elementfluxvariablescache.hh b/dumux/discretization/cellcentered/mpfa/elementfluxvariablescache.hh
index ed0c3e10e681e8b522b7106aa7ba548289cdd9a6..0d8352e5e0da8f5dc234182539d630d42c4e4947 100644
--- a/dumux/discretization/cellcentered/mpfa/elementfluxvariablescache.hh
+++ b/dumux/discretization/cellcentered/mpfa/elementfluxvariablescache.hh
@@ -44,6 +44,9 @@ class CCMpfaElementFluxVariablesCache;
template<class TypeTag>
class CCMpfaElementFluxVariablesCache<TypeTag, true>
{
+ // the local jacobian needs to be able to update the cache during assembly
+ friend typename GET_PROP_TYPE(TypeTag, LocalJacobian);
+
using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
using IndexType = typename GridView::IndexSet::IndexType;
@@ -68,11 +71,6 @@ public:
const FVElementGeometry& fvGeometry,
const ElementVolumeVariables& elemVolVars) {}
- // Specialization for the global caching being enabled - do nothing here
- void update(const Element& element,
- const FVElementGeometry& fvGeometry,
- const ElementVolumeVariables& elemVolVars) {}
-
// Specialization for the global caching being enabled - do nothing here
void bindScvf(const Element& element,
const FVElementGeometry& fvGeometry,
@@ -89,6 +87,11 @@ public:
private:
const GlobalFluxVariablesCache* globalFluxVarsCachePtr_;
+
+ // Specialization for the global caching being enabled - do nothing here
+ void update(const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars) {}
};
/*!
@@ -98,6 +101,9 @@ private:
template<class TypeTag>
class CCMpfaElementFluxVariablesCache<TypeTag, false>
{
+ // the local jacobian needs to be able to update the cache during assembly
+ friend typename GET_PROP_TYPE(TypeTag, LocalJacobian);
+
using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
using IndexType = typename GridView::IndexSet::IndexType;
@@ -119,6 +125,7 @@ public:
const FVElementGeometry& fvGeometry,
const ElementVolumeVariables& elemVolVars)
{
+ // TODO
DUNE_THROW(Dune::InvalidStateException, "Does local flux var cache binding make sense in general?");
}
@@ -129,52 +136,53 @@ public:
const ElementVolumeVariables& elemVolVars)
{
clear();
+
const auto& problem = globalFluxVarsCache().problem_();
const auto& globalFvGeometry = problem.model().globalFvGeometry();
+ const auto& neighborStencil = problem.model().stencils(element).neighborStencil();
+ const auto& assemblyMap = problem.model().localJacobian().assemblyMap();
+ const auto globalI = problem.elementMapper().index(element);
- // // first prepare all the global indices
+ // reserve initial guess of memory (won't be enough though - several scvfs per neighbor will be required)
+ globalScvfIndices_.reserve(fvGeometry.numScvf() + assemblyMap[globalI].size());
+
+ // first add all the indices inside the element
for (auto&& scvf : scvfs(fvGeometry))
- {
- const bool boundary = globalFvGeometry.scvfTouchesBoundary(scvf);
+ globalScvfIndices_.push_back(scvf.index());
+
+ // for the indices in the neighbors, use assembly map of the local jacobian
+ for (unsigned int j = 0; j < neighborStencil.size(); ++j)
+ for (auto fluxVarIdx : assemblyMap[globalI][j])
+ globalScvfIndices_.push_back(fluxVarIdx);
- 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());
globalScvfIndices_.erase(std::unique(globalScvfIndices_.begin(), globalScvfIndices_.end()), globalScvfIndices_.end());
- // prepare all the caches of the scvfs inside the corresponding interaction volume using helper class
+ // prepare all the caches of the scvfs inside the corresponding interaction volumes using helper class
fluxVarsCache_.resize(globalScvfIndices_.size());
for (auto&& scvf : scvfs(fvGeometry))
{
if (!(*this)[scvf].isUpdated())
FluxVariablesCacheFiller::fillFluxVarCache(problem, element, fvGeometry, elemVolVars, scvf, *this);
}
- }
- // This function updates the transmissibilities after the solution has been deflected during jacobian assembly
- void update(const Element& element,
- const FVElementGeometry& fvGeometry,
- const ElementVolumeVariables& elemVolVars)
- {
- for (auto&& scvf : scvfs(fvGeometry))
- (*this)[scvf].setOutdated();
-
- for (auto&& scvf : scvfs(fvGeometry))
+ // prepare the caches in the remaining neighbors
+ unsigned int j = 0;
+ for (auto globalJ : neighborStencil)
{
- if (!(*this)[scvf].isUpdated())
- FluxVariablesCacheFiller::updateFluxVarCache(globalFluxVarsCache().problem_(), element, fvGeometry, elemVolVars, scvf, *this);
+ for (auto fluxVarIdx : assemblyMap[globalI][j])
+ {
+ const auto& scvf = fvGeometry.scvf(fluxVarIdx);
+ if (!(*this)[scvf].isUpdated())
+ {
+ auto elementJ = globalFvGeometry.element(globalJ);
+ FluxVariablesCacheFiller::fillFluxVarCache(problem, elementJ, fvGeometry, elemVolVars, scvf, *this);
+ }
+ }
+
+ // increment counter
+ j++;
}
}
@@ -183,6 +191,7 @@ public:
const ElementVolumeVariables& elemVolVars,
const SubControlVolumeFace& scvf)
{
+ // TODO
DUNE_THROW(Dune::InvalidStateException, "Does binding of one scvf make sense in general?");
}
@@ -212,11 +221,25 @@ private:
globalScvfIndices_.clear();
}
+ // This function updates the transmissibilities after the solution has been deflected during jacobian assembly
+ void update(const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars)
+ {
+ for (auto&& scvf : scvfs(fvGeometry))
+ (*this)[scvf].setUpdateStatus(false);
+
+ for (auto&& scvf : scvfs(fvGeometry))
+ {
+ if (!(*this)[scvf].isUpdated())
+ FluxVariablesCacheFiller::updateFluxVarCache(globalFluxVarsCache().problem_(), element, fvGeometry, elemVolVars, scvf, *this);
+ }
+ }
+
// get index of scvf in the local container
int getLocalScvfIdx_(const int scvfIdx) const
{
auto it = std::lower_bound(globalScvfIndices_.begin(), globalScvfIndices_.end(), scvfIdx);
- assert(it != globalScvfIndices_.end() && "Could not find the flux vars cache for scvfIdx");
assert(globalScvfIndices_[std::distance(globalScvfIndices_.begin(), it)] == scvfIdx && "Could not find the flux vars cache for scvfIdx");
return std::distance(globalScvfIndices_.begin(), it);
}
diff --git a/dumux/discretization/cellcentered/mpfa/elementvolumevariables.hh b/dumux/discretization/cellcentered/mpfa/elementvolumevariables.hh
index 1df49787cfc5aa9561468ed171e8a7609a1ae7b5..8dc18e813da59e936bb80cf5f7c8e866d22df51d 100644
--- a/dumux/discretization/cellcentered/mpfa/elementvolumevariables.hh
+++ b/dumux/discretization/cellcentered/mpfa/elementvolumevariables.hh
@@ -91,6 +91,7 @@ class CCMpfaElementVolumeVariables<TypeTag, /*enableGlobalVolVarsCache*/false>
{
using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+ using MpfaHelper = typename GET_PROP_TYPE(TypeTag, MpfaHelper);
using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector);
using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
using GlobalVolumeVariables = typename GET_PROP_TYPE(TypeTag, GlobalVolumeVariables);
@@ -114,7 +115,7 @@ public:
const SolutionVector& sol)
{
const auto& problem = globalVolVars().problem_();
- auto eIdx = problem.elementMapper().index(element);
+ const auto& globalFvGeometry = problem.model().globalFvGeometry();
// stencil information
const auto& neighborStencil = problem.model().stencils(element).neighborStencil();
@@ -126,6 +127,7 @@ public:
int localIdx = 0;
// update the volume variables of the element at hand
+ auto eIdx = problem.elementMapper().index(element);
auto&& scvI = fvGeometry.scv(eIdx);
volumeVariables_[localIdx].update(sol[eIdx], problem, element, scvI);
volVarIndices_[localIdx] = scvI.index();
@@ -134,86 +136,92 @@ public:
// Update the volume variables of the neighboring elements
for (auto globalJ : neighborStencil)
{
- const auto& elementJ = fvGeometry.globalFvGeometry().element(globalJ);
+ const auto& elementJ = globalFvGeometry.element(globalJ);
auto&& scvJ = fvGeometry.scv(globalJ);
volumeVariables_[localIdx].update(sol[globalJ], problem, elementJ, scvJ);
volVarIndices_[localIdx] = scvJ.index();
++localIdx;
}
- std::size_t neighborScvfEstimate = neighborStencil.size()*dim*(2*dim-2);
- std::vector<IndexType> finishedBoundaries;
- finishedBoundaries.reserve(neighborScvfEstimate);
+ // if the element is connected to a boundary, additionally prepare BCs
+ bool boundary = element.hasBoundaryIntersections();
+ if (!boundary)
+ for (auto&& scvf : scvfs(fvGeometry))
+ if (globalFvGeometry.scvfTouchesBoundary(scvf))
+ {
+ boundary = true;
+ break;
+ }
- // if the element is connected to a boundary, additionally treat the BC
- if (element.hasBoundaryIntersections())
+ if (boundary)
{
- // reserve enough space for the boundary volume variables
- volumeVariables_.reserve(numDofs+neighborScvfEstimate);
- volVarIndices_.reserve(numDofs+neighborScvfEstimate);
-
- for (auto&& scvf : scvfs(fvGeometry))
+ // reserve memory (12 is the case of 3d hexahedron with one face on boundary)
+ std::size_t estimate = 12;
+ std::vector<IndexType> finishedBoundaries;
+ finishedBoundaries.reserve(estimate);
+ volumeVariables_.reserve(numDofs+estimate);
+ volVarIndices_.reserve(numDofs+estimate);
+
+ // treat the BCs inside the element
+ if (element.hasBoundaryIntersections())
{
- // if we are not on a boundary, skip to the next scvf
- if (!scvf.boundary())
- continue;
-
- // boundary volume variables
- VolumeVariables dirichletVolVars;
- const auto dirichletPriVars = problem.dirichlet(element, scvf);
- dirichletVolVars.update(dirichletPriVars, problem, element, scvI);
- volumeVariables_.emplace_back(std::move(dirichletVolVars));
-
- // boundary vol var index
- auto bVolVarIdx = scvf.outsideScvIdx();
- volVarIndices_.push_back(bVolVarIdx);
- finishedBoundaries.push_back(bVolVarIdx);
+ for (auto&& scvf : scvfs(fvGeometry))
+ {
+ // if we are not on a boundary, skip to the next scvf
+ if (!scvf.boundary())
+ continue;
+
+ // boundary volume variables
+ VolumeVariables dirichletVolVars;
+ const auto dirichletPriVars = problem.dirichlet(element, scvf);
+ dirichletVolVars.update(dirichletPriVars, problem, element, scvI);
+ volumeVariables_.emplace_back(std::move(dirichletVolVars));
+
+ // boundary vol var index
+ auto bVolVarIdx = scvf.outsideScvIdx();
+ volVarIndices_.push_back(bVolVarIdx);
+ finishedBoundaries.push_back(bVolVarIdx);
+ }
}
- }
- // Update boundary volume variables in the neighbors
- const auto& globalFvGeometry = problem.model().globalFvGeometry();
- for (auto&& scvf : scvfs(fvGeometry))
- {
- // reserve enough space for the boundary volume variables
- 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);
- if (!boundary)
- continue;
-
- // loop over all the scvfs in the interaction region
- const auto& ivSeed = globalFvGeometry.boundaryInteractionVolumeSeed(scvf);
- for (auto scvfIdx : ivSeed.globalScvfIndices())
+ // Update boundary volume variables in the neighbors
+ for (auto&& scvf : scvfs(fvGeometry))
{
- auto&& ivScvf = fvGeometry.scvf(scvfIdx);
-
- if (!ivScvf.boundary() || contains(ivScvf.outsideScvIdx(), finishedBoundaries))
+ // skip the rest if the scvf does not touch a boundary
+ if (!globalFvGeometry.scvfTouchesBoundary(scvf))
continue;
- // that means we are on a not yet handled boundary scvf
- auto insideScvIdx = ivScvf.insideScvIdx();
- auto&& ivScv = fvGeometry.scv(insideScvIdx);
- auto ivElement = globalFvGeometry.element(insideScvIdx);
-
- // boundary volume variables
- VolumeVariables dirichletVolVars;
- const auto dirichletPriVars = problem.dirichlet(ivElement, ivScvf);
- dirichletVolVars.update(dirichletPriVars, problem, ivElement, ivScv);
- volumeVariables_.emplace_back(std::move(dirichletVolVars));
-
- // boundary vol var index
- auto bVolVarIdx = ivScvf.outsideScvIdx();
- volVarIndices_.push_back(bVolVarIdx);
- finishedBoundaries.push_back(bVolVarIdx);
+ // 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);
+
+ if (!ivScvf.boundary() || MpfaHelper::contains(finishedBoundaries, ivScvf.outsideScvIdx()))
+ continue;
+
+ // that means we are on a not yet handled boundary scvf
+ auto insideScvIdx = ivScvf.insideScvIdx();
+ auto&& ivScv = fvGeometry.scv(insideScvIdx);
+ auto ivElement = globalFvGeometry.element(insideScvIdx);
+
+ // boundary volume variables
+ VolumeVariables dirichletVolVars;
+ const auto dirichletPriVars = problem.dirichlet(ivElement, ivScvf);
+ dirichletVolVars.update(dirichletPriVars, problem, ivElement, ivScv);
+ volumeVariables_.emplace_back(std::move(dirichletVolVars));
+
+ // boundary vol var index
+ auto bVolVarIdx = ivScvf.outsideScvIdx();
+ volVarIndices_.push_back(bVolVarIdx);
+ finishedBoundaries.push_back(bVolVarIdx);
+ }
}
- }
- // free unused memory
- volumeVariables_.shrink_to_fit();
- volVarIndices_.shrink_to_fit();
+ // free unused memory
+ volumeVariables_.shrink_to_fit();
+ volVarIndices_.shrink_to_fit();
+ }
}
// Binding of an element, prepares only the volume variables of the element
@@ -258,10 +266,6 @@ private:
return std::distance(volVarIndices_.begin(), it);
}
- const bool contains(const IndexType idx,
- const std::vector<IndexType>& indices) const
- { return std::find(indices.begin(), indices.end(), idx) != indices.end(); }
-
std::vector<IndexType> volVarIndices_;
std::vector<VolumeVariables> volumeVariables_;
};
diff --git a/dumux/discretization/cellcentered/mpfa/facetypes.hh b/dumux/discretization/cellcentered/mpfa/facetypes.hh
index f1768e84be166ee38cbb12f11d05524614f54cd6..8da25323436a11e6e15353e9f3853f38ef69c04e 100644
--- a/dumux/discretization/cellcentered/mpfa/facetypes.hh
+++ b/dumux/discretization/cellcentered/mpfa/facetypes.hh
@@ -29,9 +29,7 @@ namespace Dumux
{
interior,
neumann,
- dirichlet,
- interiorNeumann,
- interiorDirichlet
+ dirichlet
};
} // end namespace
diff --git a/dumux/discretization/cellcentered/mpfa/fluxvariablescachefiller.hh b/dumux/discretization/cellcentered/mpfa/fluxvariablescachefiller.hh
index 1db15ed006c0aaed16f68d33cdcd04bb3d9d3db9..bc27e2f94d96071e683f9d5a9d25d962f555cb1a 100644
--- a/dumux/discretization/cellcentered/mpfa/fluxvariablescachefiller.hh
+++ b/dumux/discretization/cellcentered/mpfa/fluxvariablescachefiller.hh
@@ -108,7 +108,7 @@ public:
{
if (!updateNeumannOnly)
cache.updateAdvection(problem, element, fvGeometry, elemVolVars, scvf, iv);
- cache.setUpdated();
+ cache.setUpdateStatus(true);
}
// update flux variable caches of the other scvfs of the interaction volume
@@ -124,7 +124,7 @@ public:
{
if (!updateNeumannOnly)
cacheJ.updateAdvection(problem, elementJ, fvGeometry, elemVolVars, scvfJ, iv);
- cacheJ.setUpdated();
+ cacheJ.setUpdateStatus(true);
}
}
}
@@ -140,7 +140,7 @@ public:
const auto scvfIdx = scvf.index();
auto& cache = fluxVarsCache[scvfIdx];
cache.updateAdvection(problem, element, fvGeometry, elemVolVars, scvf, iv);
- cache.setUpdated();
+ cache.setUpdateStatus(true);
// update flux variable caches of the other scvfs of the interaction volume
for (const auto scvfIdxJ : iv.globalScvfs())
@@ -151,7 +151,7 @@ public:
const auto elementJ = problem.model().globalFvGeometry().element(scvfJ.insideScvIdx());
auto& cacheJ = fluxVarsCache[scvfIdxJ];
cacheJ.updateAdvection(problem, elementJ, fvGeometry, elemVolVars, scvfJ, iv);
- cacheJ.setUpdated();
+ cacheJ.setUpdateStatus(true);
}
}
}
@@ -168,7 +168,7 @@ public:
{
if (!solDependentParams)
{
- fluxVarsCache[scvf.index()].setUpdated();
+ fluxVarsCache[scvf.index()].setUpdateStatus(true);
// if we do not use tpfa on the boundary, we have to update the neumann fluxes
if (!useTpfaBoundary)
diff --git a/dumux/discretization/cellcentered/mpfa/fvelementgeometry.hh b/dumux/discretization/cellcentered/mpfa/fvelementgeometry.hh
index 83982b7d3b5f5670c5d3bb9fd48b39160d1cd737..9dcd2ff583d4f0175ee101840c45b7feea8d6451 100644
--- a/dumux/discretization/cellcentered/mpfa/fvelementgeometry.hh
+++ b/dumux/discretization/cellcentered/mpfa/fvelementgeometry.hh
@@ -30,8 +30,6 @@
#include <dumux/discretization/scvandscvfiterators.hh>
#include <dumux/implicit/cellcentered/mpfa/properties.hh>
-#include "mpfageometryhelper.hh"
-
namespace Dumux
{
/*!
@@ -148,6 +146,7 @@ class CCMpfaFVElementGeometry<TypeTag, false>
using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
using IndexType = typename GridView::IndexSet::IndexType;
+ using MpfaHelper = typename GET_PROP_TYPE(TypeTag, MpfaHelper);
using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
using Element = typename GridView::template Codim<0>::Entity;
@@ -157,10 +156,10 @@ class CCMpfaFVElementGeometry<TypeTag, false>
using ScvfIterator = Dumux::ScvfIterator<SubControlVolumeFace, std::vector<IndexType>, ThisType>;
static const int dim = GridView::dimension;
+ static const int dimWorld = GridView::dimensionworld;
using CoordScalar = typename GridView::ctype;
using ReferenceElements = typename Dune::ReferenceElements<CoordScalar, dim>;
-
- using MpfaGeometryHelper = Dumux::MpfaGeometryHelper<GridView, dim>;
+ using GlobalPosition = Dune::FieldVector<Scalar, dimWorld>;
public:
//! Constructor
@@ -226,36 +225,34 @@ public:
{
bindElement(element);
+ // get some references for convenience
+ const auto& problem = globalFvGeometry().problem_();
+ const auto& assemblyMap = problem.model().localJacobian().assemblyMap();
+ const auto& neighborStencil = problem.model().stencils(element).neighborStencil();
+ const auto globalI = problem.elementMapper().index(element);
+
// reserve memory
- auto numNeighbors = globalFvGeometry().problem_().model().stencils(element).neighborStencil().size();
- auto numNeighborScvf = numNeighbors*dim*(2*dim-2);
+ const auto numNeighbors = neighborStencil.size();
+ const auto estimate = numNeighbors*dim*(4*dim-4); // estimate holds for quadrilaterals in 2D/3D, overestimates else
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_)
+ neighborScvfIndices_.reserve(estimate);
+ neighborScvfs_.reserve(estimate);
+
+ // build scvfs in neighbors
+ // use the assembly map to determine which faces are necessary
+ unsigned int j = 0;
+ for (auto globalJ : neighborStencil)
{
- // skip the rest if scvf has been handled already
- if (contains(scvf.vertexIndex(), finishedVertices))
- continue;
+ // get data on the neighbor
+ auto elementJ = globalFvGeometry().element(globalJ);
+ const auto& scvfIndices = assemblyMap[globalI][j];
- if (globalFvGeometry().scvfTouchesBoundary(scvf))
- {
- const auto& ivSeed = globalFvGeometry().boundaryInteractionVolumeSeed(scvf);
- handleInteractionRegion(ivSeed);
- }
- else
- {
- const auto& ivSeed = globalFvGeometry().interactionVolumeSeed(scvf);
- handleInteractionRegion(ivSeed);
- }
+ // make neighbor geometries
+ makeNeighborGeometries(elementJ, globalJ, scvfIndices);
- finishedVertices.push_back(scvf.vertexIndex());
+ // increment counter
+ j++;
}
// free unused memory
@@ -294,12 +291,8 @@ private:
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);
+ // Instantiate helper class to pass it to scvf constructors
+ MpfaHelper helper;
// the quadrature point to be used on the scvf
const Scalar q = GET_PARAM_FROM_GROUP(TypeTag, Scalar, Mpfa, Q);
@@ -309,114 +302,139 @@ private:
scvfIndices_.reserve(numLocalScvf);
scvfs_.reserve(numLocalScvf);
+ // for dim < dimWorld we'll have multiple intersections at one point
+ // we store here the centers of those we have handled already
+ std::vector<GlobalPosition> finishedCenters;
+
int scvfCounter = 0;
for (const auto& is : intersections(globalFvGeometry().gridView(), element))
{
- auto isGeometry = is.geometry();
+ // if we are dealing with a lower dimensional network
+ // only make a new scvf if we haven't handled it yet
+ if (dim < dimWorld)
+ {
+ auto isCenter = is.geometry().center();
+ if(MpfaHelper::contains(finishedCenters, isCenter))
+ continue;
+ else
+ finishedCenters.push_back(isCenter);
+ }
+
+ // get the intersection corners according to generic numbering
+ auto numCorners = is.geometry().corners();
+ std::vector<GlobalPosition> isCorners(numCorners);
+
+ // if outside level > inside level, use the outside element in the following
+ bool useNeighbor = is.neighbor() && is.outside().level() > element.level();
+ const auto& e = useNeighbor ? is.outside() : element;
+ const auto indexInElement = useNeighbor ? is.indexInOutside() : is.indexInInside();
+ const auto eg = e.geometry();
+ const auto& refElement = ReferenceElements::general(eg.type());
- // make the scv faces of the intersection
- for (unsigned int faceScvfIdx = 0; faceScvfIdx < isGeometry.corners(); ++faceScvfIdx)
+ for (unsigned int c = 0; c < numCorners; ++c)
+ isCorners[c] = eg.global(refElement.position(refElement.subEntity(indexInElement, 1, c, dim), dim));
+
+ // make the scv faces belonging to each corner of the intersection
+ for (unsigned int c = 0; c < numCorners; ++c)
{
- // 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);
+ // get the global vertex index the scv face is connected to
+ auto vIdxLocal = refElement.subEntity(indexInElement, 1, c, dim);
+ auto vIdxGlobal = problem.vertexMapper().subIndex(e, vIdxLocal, dim);
// do not build scvfs connected to a processor boundary
if (globalFvGeometry().isGhostVertex(vIdxGlobal))
continue;
- scvfs_.emplace_back(geomHelper,
- geomHelper.getScvfCorners(isGeometry, faceScvfIdx),
+ scvfs_.emplace_back(helper,
+ helper.getScvfCorners(isCorners, c),
is.centerUnitOuterNormal(),
vIdxGlobal,
vIdxLocal,
scvFaceIndices[scvfCounter],
- std::array<IndexType, 2>({{eIdx, neighborVolVarIndices[scvfCounter]}}),
+ eIdx,
+ neighborVolVarIndices[scvfCounter],
q,
is.boundary()
);
scvfIndices_.emplace_back(scvFaceIndices[scvfCounter]);
scvfCounter++;
-
}
}
-
- // free unused memory
- scvfs_.shrink_to_fit();
- 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
+ //! create the scv and necessary scvfs of a neighboring element
template<typename IndexVector>
- void makeNeighborGeometries(const Element& element, IndexType eIdxGlobal, const IndexVector& ivScvfs)
+ void makeNeighborGeometries(const Element& element, IndexType eIdxGlobal, const IndexVector& scvfIndices)
{
// create the neighbor scv if it doesn't exist yet
- if (!contains(eIdxGlobal, neighborScvIndices_))
- {
- neighborScvs_.emplace_back(element.geometry(), eIdxGlobal);
- neighborScvIndices_.push_back(eIdxGlobal);
- }
+ neighborScvs_.emplace_back(element.geometry(), eIdxGlobal);
+ neighborScvIndices_.push_back(eIdxGlobal);
+ // get data on the scv faces
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);
+ // Instantiate a helper class to pass it to scvf constructors
+ MpfaHelper helper;
// the quadrature point to be used on the scvf
const Scalar q = GET_PARAM_FROM_GROUP(TypeTag, Scalar, Mpfa, Q);
+ // for dim < dimWorld we'll have multiple intersections at one point
+ // we store here the centers of those we have handled already
+ std::vector<GlobalPosition> finishedCenters;
+
int scvfCounter = 0;
for (const auto& is : intersections(globalFvGeometry().gridView(), element))
{
- auto isGeometry = is.geometry();
+ // if we are dealing with a lower dimensional network
+ // only make a new scvf if we haven't handled it yet
+ if (dim < dimWorld)
+ {
+ auto isCenter = is.geometry().center();
+ if(MpfaHelper::contains(finishedCenters, isCenter))
+ continue;
+ else
+ finishedCenters.push_back(isCenter);
+ }
+
+ // get the intersection corners according to generic numbering
+ auto numCorners = is.geometry().corners();
+ std::vector<GlobalPosition> isCorners(numCorners);
+
+ // if outside level > inside level, use the outside element in the following
+ bool useNeighbor = is.neighbor() && is.outside().level() > element.level();
+ const auto& e = useNeighbor ? is.outside() : element;
+ const auto indexInElement = useNeighbor ? is.indexInOutside() : is.indexInInside();
+ const auto eg = e.geometry();
+ const auto& refElement = ReferenceElements::general(eg.type());
- // make the scv faces of the intersection
- for (unsigned int faceScvfIdx = 0; faceScvfIdx < isGeometry.corners(); ++faceScvfIdx)
+ for (unsigned int c = 0; c < numCorners; ++c)
+ isCorners[c] = eg.global(refElement.position(refElement.subEntity(indexInElement, 1, c, dim), dim));
+
+ // make the scv faces belonging to each corner of the intersection
+ for (unsigned int c = 0; c < numCorners; ++c)
{
- // 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);
+ // get the global vertex index the scv face is connected to
+ auto vIdxLocal = refElement.subEntity(indexInElement, 1, c, dim);
+ auto vIdxGlobal = globalFvGeometry().problem_().vertexMapper().subIndex(e, 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 (contains(scvFaceIndices[scvfCounter], ivScvfs))
+ // only build the scvf if it is in the list of necessary indices
+ if (MpfaHelper::contains(scvfIndices, scvFaceIndices[scvfCounter]))
{
- neighborScvfs_.emplace_back(geomHelper,
- geomHelper.getScvfCorners(isGeometry, faceScvfIdx),
+ neighborScvfs_.emplace_back(helper,
+ helper.getScvfCorners(isCorners, c),
is.centerUnitOuterNormal(),
vIdxGlobal,
vIdxLocal,
scvFaceIndices[scvfCounter],
- std::array<IndexType, 2>({{eIdxGlobal, neighborVolVarIndices[scvfCounter]}}),
+ eIdxGlobal,
+ neighborVolVarIndices[scvfCounter],
q,
is.boundary()
);
@@ -424,6 +442,7 @@ private:
neighborScvfIndices_.emplace_back(scvFaceIndices[scvfCounter]);
}
+ // increment counter either way
scvfCounter++;
}
}
@@ -437,11 +456,6 @@ private:
return std::distance(indices.begin(), it);
}
- //! Returns whether or not an index is stored in a given vector
- const bool contains(const IndexType idx,
- const std::vector<IndexType>& indices) const
- { return std::find(indices.begin(), indices.end(), idx) != indices.end(); }
-
void clear()
{
scvIndices_.clear();
diff --git a/dumux/discretization/cellcentered/mpfa/globalfvgeometry.hh b/dumux/discretization/cellcentered/mpfa/globalfvgeometry.hh
index e8c8a148fc47cca4ef2253e1dcee27a94044db54..e8e02feadc1cad1f7c31ed1c70cdf00908ab6949 100644
--- a/dumux/discretization/cellcentered/mpfa/globalfvgeometry.hh
+++ b/dumux/discretization/cellcentered/mpfa/globalfvgeometry.hh
@@ -25,14 +25,23 @@
#ifndef DUMUX_DISCRETIZATION_CC_MPFA_GLOBALFVGEOMETRY_HH
#define DUMUX_DISCRETIZATION_CC_MPFA_GLOBALFVGEOMETRY_HH
+#include <dumux/discretization/cellcentered/mpfa/globalfvgeometrybase.hh>
#include "methods.hh"
namespace Dumux
{
-//! forward declaration of the method-specific specialization
+//! Forward declaration of the method-specific specialization
+//! By default, the methods simply inherit from the base class
+//! Methods that use a different implementation have to provide this below
template<class TypeTag, MpfaMethods method, bool EnableFVElementGeometryCache>
-class CCMpfaGlobalFVGeometryImplementation
-{};
+class CCMpfaGlobalFVGeometryImplementation : public CCMpfaGlobalFVGeometryBase<TypeTag, EnableFVElementGeometryCache>
+{
+ using ParentType = CCMpfaGlobalFVGeometryBase<TypeTag, EnableFVElementGeometryCache>;
+ using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+
+public:
+ CCMpfaGlobalFVGeometryImplementation(const GridView& gridView) : ParentType(gridView) {}
+};
/*!
* \ingroup ImplicitModel
@@ -45,8 +54,6 @@ using CCMpfaGlobalFVGeometry = CCMpfaGlobalFVGeometryImplementation<TypeTag, GET
} // end namespace
-// the specializations of this class for the available methods have to be included here
-#include <dumux/discretization/cellcentered/mpfa/omethod/globalfvgeometry.hh>
-#include <dumux/discretization/cellcentered/mpfa/omethodfps/globalfvgeometry.hh>
+// further specializations of this class for the available methods have to be included here
#endif
diff --git a/dumux/discretization/cellcentered/mpfa/globalfvgeometrybase.hh b/dumux/discretization/cellcentered/mpfa/globalfvgeometrybase.hh
index 0f9e80e1ff92c7ea99324e4fa77ec0ef5049ed30..ef5420f94e8f38d68ef10e1a0dc75bbfc15749c4 100644
--- a/dumux/discretization/cellcentered/mpfa/globalfvgeometrybase.hh
+++ b/dumux/discretization/cellcentered/mpfa/globalfvgeometrybase.hh
@@ -18,7 +18,7 @@
*****************************************************************************/
/*!
* \file
- * \brief Base class for the finite volume geometry vector for box models
+ * \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.
*/
@@ -31,82 +31,8 @@
#include <dumux/implicit/cellcentered/mpfa/properties.hh>
#include <dumux/common/elementmap.hh>
-#include "mpfageometryhelper.hh"
-
namespace Dumux
{
-//! Helper class for the global fv geometries
-template<class TypeTag>
-class CCMpfaGlobalFVGeometryHelper
-{
- using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
- using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
- using CoordScalar = typename GridView::ctype;
-
- static const int dim = GridView::dimension;
- using ReferenceElements = typename Dune::ReferenceElements<CoordScalar, dim>;
-
-public:
- static std::vector<bool> findGhostVertices(const Problem& problem, const GridView& gridView)
- {
- std::vector<bool> ghostVertices(gridView.size(dim), false);
-
- // if not run in parallel, skip the rest
- if (Dune::MPIHelper::getCollectiveCommunication().size() == 1)
- return ghostVertices;
-
- // mpfa methods can not yet handle ghost cells
- if (gridView.ghostSize(0) > 0)
- DUNE_THROW(Dune::InvalidStateException, "Mpfa methods in parallel do not work with ghost cells. Use overlap cells instead.");
-
- // mpfa methods have to have overlapping cells
- if (gridView.overlapSize(0) == 0)
- DUNE_THROW(Dune::InvalidStateException, "Grid no overlaping cells. This is required by mpfa methods in parallel.");
-
- for (const auto& element : elements(gridView))
- {
- for (const auto& is : intersections(gridView, element))
- {
- if (!is.neighbor() && !is.boundary())
- {
- const auto& refElement = ReferenceElements::general(element.geometry().type());
- for (unsigned int isVertex = 0; isVertex < is.geometry().corners(); ++isVertex)
- {
- const auto vIdxLocal = refElement.subEntity(is.indexInInside(), 1, isVertex, dim);
- const auto vIdxGlobal = problem.vertexMapper().subIndex(element, vIdxLocal, dim);
-
- ghostVertices[vIdxGlobal] = true;
- }
- }
- }
- }
-
- return ghostVertices;
- }
-
- template<int d = dim>
- static typename std::enable_if<d == 2, std::size_t>::type getGlobalNumScvf(const GridView& gridView)
- {
- Dune::GeometryType triangle, quadrilateral;
- triangle.makeTriangle();
- quadrilateral.makeQuadrilateral();
-
- return gridView.size(triangle)*6 + gridView.size(quadrilateral)*8;
- }
-
- template<int d = dim>
- static typename std::enable_if<d == 3, std::size_t>::type getGlobalNumScvf(const GridView& gridView)
- {
- 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;
- }
-};
-
/*!
* \ingroup ImplicitModel
* \brief Base class for the finite volume geometry vector for mpfa models
@@ -121,13 +47,16 @@ class CCMpfaGlobalFVGeometryBase
template<class TypeTag>
class CCMpfaGlobalFVGeometryBase<TypeTag, true>
{
- //! The local class needs access to the scv, scvfs and the fv element geometry
- //! as they are globally cached
+ using Implementation = typename GET_PROP_TYPE(TypeTag, GlobalFVGeometry);
+
+ //! The actual implementation might overwrite the update() routine
+ friend Implementation;
+ //! The local class 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 MpfaHelper = typename GET_PROP_TYPE(TypeTag, MpfaHelper);
using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
@@ -146,9 +75,8 @@ class CCMpfaGlobalFVGeometryBase<TypeTag, true>
static const int dim = GridView::dimension;
static const int dimWorld = GridView::dimensionworld;
- using DimVector = Dune::FieldVector<Scalar, dimWorld>;
+ using GlobalPosition = Dune::FieldVector<Scalar, dimWorld>;
using ReferenceElements = typename Dune::ReferenceElements<CoordScalar, dim>;
- using MpfaGeometryHelper = Dumux::MpfaGeometryHelper<GridView, dim>;
public:
//! Constructor
@@ -201,7 +129,7 @@ public:
// reserve memory
std::size_t numScvs = gridView_.size(0);
- std::size_t numScvfs = CCMpfaGlobalFVGeometryHelper<TypeTag>::getGlobalNumScvf(gridView_);
+ std::size_t numScvfs = MpfaHelper::getGlobalNumScvf(gridView_);
scvs_.resize(numScvs);
scvfs_.reserve(numScvfs);
scvfIndicesOfScv_.resize(numScvs);
@@ -209,11 +137,14 @@ public:
elementMap_.resize(numScvs);
// find vertices on processor boundaries
- auto ghostVertices = CCMpfaGlobalFVGeometryHelper<TypeTag>::findGhostVertices(problem, gridView_);
+ auto isGhostVertex = MpfaHelper::findGhostVertices(problem, gridView_);
// the quadrature point to be used on the scvf
const Scalar q = GET_PARAM_FROM_GROUP(TypeTag, Scalar, Mpfa, Q);
+ // Instantiate a helper class to pass it to the scvf constructors
+ MpfaHelper helper;
+
// Build the SCVs and SCV faces
IndexType scvfIdx = 0;
numBoundaryScvf_ = 0;
@@ -221,48 +152,89 @@ public:
{
auto eIdx = problem.elementMapper().index(element);
+ // the element geometry
+ auto elementGeometry = element.geometry();
+
// 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);
-
- // create sub control volume for this element
- scvs_[eIdx] = SubControlVolume(std::move(elemGeometry), eIdx);
-
// The local scvf index set
std::vector<IndexType> scvfIndexSet;
- scvfIndexSet.reserve(geomHelper.getNumLocalScvfs());
+ scvfIndexSet.reserve(MpfaHelper::getNumLocalScvfs(elementGeometry.type()));
+
+ // for dim < dimWorld we'll have multiple intersections at one point
+ // we store here the centers of those we have handled already
+ std::vector<GlobalPosition> finishedCenters;
// construct the sub control volume faces
- for (const auto& intersection : intersections(gridView_, element))
+ for (const auto& is : intersections(gridView_, element))
{
- // get the intersection geometry and some of its bools
- auto isGeometry = intersection.geometry();
- bool boundary = intersection.boundary();
-
- // determine the outside volvar idx
- // we initialize nIdx so that the compiler does not complain in overlap elements in parallel
- // nIdx will not be used in this case anyway (see continue call below)
- IndexType nIdx = eIdx;
- if (intersection.neighbor())
- nIdx = problem.elementMapper().index(intersection.outside());
+ // if we are dealing with a lower dimensional network
+ // only make a new scvf if we haven't handled it yet
+ if (dim < dimWorld)
+ if(MpfaHelper::contains(finishedCenters, is.geometry().center()))
+ continue;
+
+ // get some of the intersection's bools
+ bool boundary = is.boundary();
+ bool neighbor = is.neighbor();
+
+ // determine the outside volvar indices
+ std::vector<IndexType> nIndices;
+ if (neighbor)
+ {
+ nIndices = std::vector<IndexType>( {problem.elementMapper().index(is.outside())} );
+
+ // check for further neighbors in case of lower dimensional networks
+ if (dim < dimWorld)
+ {
+ auto insideIndex = is.indexInInside();
+ auto isCenter = is.geometry().center();
+
+ // find further possible intersections at the same point (bifurcations etc)
+ for (const auto& nextIs : intersections(gridView_, element))
+ {
+ if (nextIs.indexInInside() == insideIndex && nextIs.geometry().center() == isCenter)
+ {
+ auto nIdx = problem.elementMapper().index(nextIs.outside());
+ if (nIdx != nIndices[0])
+ nIndices.push_back(nIdx);
+ }
+ }
+
+ // store this center as a handled one
+ finishedCenters.push_back(isCenter);
+ }
+ }
else if (boundary)
- nIdx = numScvs + numBoundaryScvf_++;
+ {
+ nIndices.resize(1);
+ nIndices[0] = numScvs + numBoundaryScvf_++;
+ }
- // make the scv faces of the intersection
- for (unsigned int faceScvfIdx = 0; faceScvfIdx < isGeometry.corners(); ++faceScvfIdx)
+ // get the intersection corners according to generic numbering
+ auto numCorners = is.geometry().corners();
+ std::vector<GlobalPosition> isCorners(numCorners);
+
+ // if outside level > inside level, use the outside element in the following
+ bool useNeighbor = neighbor && is.outside().level() > element.level();
+ const auto& e = useNeighbor ? is.outside() : element;
+ const auto indexInElement = useNeighbor ? is.indexInOutside() : is.indexInInside();
+ const auto eg = e.geometry();
+ const auto& refElement = ReferenceElements::general(eg.type());
+
+ for (unsigned int c = 0; c < numCorners; ++c)
+ isCorners[c] = eg.global(refElement.position(refElement.subEntity(indexInElement, 1, c, dim), dim));
+
+ // make the scv faces belonging to each corner of the intersection
+ for (unsigned int c = 0; c < numCorners; ++c)
{
- // 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);
+ // get the global vertex index the scv face is connected to
+ auto vIdxLocal = refElement.subEntity(indexInElement, 1, c, dim);
+ auto vIdxGlobal = problem.vertexMapper().subIndex(e, vIdxLocal, dim);
// do not build scvfs connected to a processor boundary
- if (ghostVertices[vIdxGlobal])
+ if (isGhostVertex[vIdxGlobal])
continue;
// store info on which vertices are on the domain boundary
@@ -271,13 +243,14 @@ public:
// make the scv face
scvfIndexSet.push_back(scvfIdx);
- scvfs_.emplace_back(geomHelper,
- geomHelper.getScvfCorners(isGeometry, faceScvfIdx),
- intersection.centerUnitOuterNormal(),
+ scvfs_.emplace_back(helper,
+ helper.getScvfCorners(isCorners, c),
+ is.centerUnitOuterNormal(),
vIdxGlobal,
vIdxLocal,
scvfIdx,
- std::array<IndexType, 2>({{eIdx, nIdx}}),
+ eIdx,
+ nIndices,
q,
boundary
);
@@ -287,6 +260,9 @@ public:
}
}
+ // create sub control volume for this element
+ scvs_[eIdx] = SubControlVolume(std::move(elementGeometry), eIdx);
+
// Save the scvf indices belonging to this scv to build up fv element geometries fast
scvfIndicesOfScv_[eIdx] = scvfIndexSet;
}
@@ -321,6 +297,10 @@ public:
const std::vector<IndexType>& scvfIndicesOfScv(IndexType scvIdx) const
{ return scvfIndicesOfScv_[scvIdx]; }
+ //! Return a const reference to the grid view
+ const GridView& gridView() const
+ { return gridView_; }
+
private:
const Problem& problem_() const
@@ -353,6 +333,7 @@ class CCMpfaGlobalFVGeometryBase<TypeTag, false>
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 MpfaHelper = typename GET_PROP_TYPE(TypeTag, MpfaHelper);
using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
@@ -371,9 +352,8 @@ class CCMpfaGlobalFVGeometryBase<TypeTag, false>
static const int dim = GridView::dimension;
static const int dimWorld = GridView::dimensionworld;
- using DimVector = Dune::FieldVector<Scalar, dimWorld>;
+ using GlobalPosition = Dune::FieldVector<Scalar, dimWorld>;
using ReferenceElements = typename Dune::ReferenceElements<CoordScalar, dim>;
- using MpfaGeometryHelper = Dumux::MpfaGeometryHelper<GridView, dim>;
public:
//! Constructor
@@ -441,7 +421,7 @@ public:
boundaryVertices_.resize(gridView_.size(dim), false);
// find vertices on processor boundaries
- ghostVertices_ = CCMpfaGlobalFVGeometryHelper<TypeTag>::findGhostVertices(problem, gridView_);
+ ghostVertices_ = MpfaHelper::findGhostVertices(problem, gridView_);
// Store necessary info on SCVs and SCV faces
for (const auto& element : elements(gridView_))
@@ -452,38 +432,72 @@ public:
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);
+ auto eg = element.geometry();
+ const auto& referenceElement = ReferenceElements::general(eg.type());
// 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;
+ auto numLocalFaces = MpfaHelper::getNumLocalScvfs(eg.type());
+ std::vector<IndexType> scvfsIndexSet;
+ std::vector< std::vector<IndexType> > neighborVolVarIndexSet;
+ scvfsIndexSet.reserve(numLocalFaces);
+ neighborVolVarIndexSet.reserve(numLocalFaces);
+
+ // for dim < dimWorld we'll have multiple intersections at one point
+ // we store here the centers of those we have handled already
+ std::vector<GlobalPosition> finishedCenters;
+
+ unsigned int localFaceIdx = 0;
// construct the sub control volume faces
- for (const auto& intersection : intersections(gridView_, element))
+ for (const auto& is : intersections(gridView_, element))
{
+ // if we are dealing with a lower dimensional network
+ // only make a new scvf if we haven't handled it yet
+ if (dim < dimWorld)
+ if(MpfaHelper::contains(finishedCenters, is.geometry().center()))
+ continue;
+
// get some of the intersection bools
- bool boundary = intersection.boundary();
- bool neighbor = intersection.neighbor();
+ bool boundary = is.boundary();
+ bool neighbor = is.neighbor();
- // determine the outside volvar idx
- // we initialize nIdx so that the compiler does not complain in overlap elements in parallel
- // nIdx will not be used in this case anyway (see continue call below)
- IndexType nIdx = eIdx;
+ // determine the outside volvar indices
+ std::vector<IndexType> nIndices;
if (neighbor)
- nIdx = problem.elementMapper().index(intersection.outside());
+ {
+ nIndices = std::vector<IndexType>( {problem.elementMapper().index(is.outside())} );
+
+ // check for further neighbors in case of lower dimensional networks
+ if (dim < dimWorld)
+ {
+ auto insideIndex = is.indexInInside();
+ auto isCenter = is.geometry().center();
+
+ // find further possible intersections at the same point (bifurcations etc)
+ for (const auto& nextIs : intersections(gridView_, element))
+ {
+ if (nextIs.indexInInside() == insideIndex && nextIs.geometry().center() == isCenter)
+ {
+ auto nIdx = problem.elementMapper().index(nextIs.outside());
+ if (nIdx != nIndices[0])
+ nIndices.push_back(nIdx);
+ }
+ }
+
+ // store this center as a handled one
+ finishedCenters.push_back(isCenter);
+ }
+ }
else if (boundary)
- nIdx = numScvs_ + numBoundaryScvf_++;
+ {
+ nIndices.resize(1);
+ nIndices[0] = numScvs_ + numBoundaryScvf_++;
+ }
// make the scv faces of the intersection
- for (unsigned int faceScvfIdx = 0; faceScvfIdx < intersection.geometry().corners(); ++faceScvfIdx)
+ for (unsigned int faceScvfIdx = 0; faceScvfIdx < is.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 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
@@ -495,8 +509,11 @@ public:
boundaryVertices_[vIdxGlobal] = true;
// store information on the scv face
- scvfsIndexSet[localFaceIdx] = numScvf_++;
- neighborVolVarIndexSet[localFaceIdx++] = nIdx;
+ scvfsIndexSet.push_back(numScvf_++);
+ neighborVolVarIndexSet.push_back(nIndices);
+
+ // increment counter
+ localFaceIdx++;
}
}
@@ -515,7 +532,7 @@ public:
const std::vector<IndexType>& scvfIndicesOfScv(IndexType scvIdx) const
{ return scvfIndicesOfScv_[scvIdx]; }
- const std::vector<IndexType>& neighborVolVarIndices(IndexType scvIdx) const
+ const std::vector< std::vector<IndexType> >& neighborVolVarIndices(IndexType scvIdx) const
{ return neighborVolVarIndices_[scvIdx]; }
/*!
@@ -543,7 +560,7 @@ private:
// vectors that store the global data
Dumux::ElementMap<GridView> elementMap_;
std::vector<std::vector<IndexType>> scvfIndicesOfScv_;
- std::vector<std::vector<IndexType>> neighborVolVarIndices_;
+ std::vector< std::vector< std::vector<IndexType> > > neighborVolVarIndices_;
std::vector<bool> boundaryVertices_;
std::vector<bool> ghostVertices_;
diff --git a/dumux/discretization/cellcentered/mpfa/globalinteractionvolumeseeds.hh b/dumux/discretization/cellcentered/mpfa/globalinteractionvolumeseeds.hh
index d1db5c93a2fc7482bf080f42d58682cd758d2bdd..87af94939467c06d328631046c2f2655949dabf4 100644
--- a/dumux/discretization/cellcentered/mpfa/globalinteractionvolumeseeds.hh
+++ b/dumux/discretization/cellcentered/mpfa/globalinteractionvolumeseeds.hh
@@ -24,14 +24,23 @@
#define DUMUX_DISCRETIZATION_MPFA_GLOBALINTERACTIONVOLUMESEEDS_HH
#include <dumux/implicit/cellcentered/mpfa/properties.hh>
+#include <dumux/discretization/cellcentered/mpfa/globalinteractionvolumeseedsbase.hh>
#include "methods.hh"
namespace Dumux
{
-//! forward declaration of the actual implementation
+//! forward declaration of the actual method-specific implementation
+//! By default we simply inherit from the base class
+//! Actual implementations for other methods have to be provided below
template<class TypeTag, MpfaMethods method>
-class CCMpfaGlobalInteractionVolumeSeedsImplementation
-{};
+class CCMpfaGlobalInteractionVolumeSeedsImplementation : public CCMpfaGlobalInteractionVolumeSeedsBase<TypeTag>
+{
+ using ParentType = CCMpfaGlobalInteractionVolumeSeedsBase<TypeTag>;
+ using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+
+public:
+ CCMpfaGlobalInteractionVolumeSeedsImplementation(const GridView gridView) : ParentType(gridView) {}
+};
/*!
* \ingroup Mpfa
@@ -42,8 +51,7 @@ using CCMpfaGlobalInteractionVolumeSeeds = CCMpfaGlobalInteractionVolumeSeedsImp
} // end namespace
-// the specializations of this class for the available methods have to be included here
-#include <dumux/discretization/cellcentered/mpfa/omethod/globalinteractionvolumeseeds.hh>
-#include <dumux/discretization/cellcentered/mpfa/omethodfps/globalinteractionvolumeseeds.hh>
+// the specializations of this class differing from the default have to be included here
+#include <dumux/discretization/cellcentered/mpfa/lmethod/globalinteractionvolumeseeds.hh>
#endif
diff --git a/dumux/discretization/cellcentered/mpfa/globalinteractionvolumeseedsbase.hh b/dumux/discretization/cellcentered/mpfa/globalinteractionvolumeseedsbase.hh
index b036f9717f49530f0680b5f1cfe5edb9a07a891d..73ec5c9a1275821ca75407f65e44ae70ca65e742 100644
--- a/dumux/discretization/cellcentered/mpfa/globalinteractionvolumeseedsbase.hh
+++ b/dumux/discretization/cellcentered/mpfa/globalinteractionvolumeseedsbase.hh
@@ -35,6 +35,10 @@ namespace Dumux
template<class TypeTag>
class CCMpfaGlobalInteractionVolumeSeedsBase
{
+ using Implementation = typename GET_PROP_TYPE(TypeTag, GlobalInteractionVolumeSeeds);
+ // the actual implementation needs to be friend
+ friend Implementation;
+
using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
using Helper = typename GET_PROP_TYPE(TypeTag, MpfaHelper);
@@ -58,12 +62,12 @@ public:
seeds_.clear();
boundarySeeds_.clear();
- // -1 indicates that the scvf has not been handled yet
auto numScvf = problem_().model().globalFvGeometry().numScvf();
- scvfIndexMap_.resize(numScvf, -1);
+ scvfIndexMap_.resize(numScvf);
+ std::vector<bool> isFaceHandled(numScvf, false);
- // detect and handle the boundary first
- initializeSeeds_(boundaryVertices);
+ // initialize the seeds according to the mpfa method
+ asImp_().initializeSeeds_(boundaryVertices, isFaceHandled);
}
const InteractionVolumeSeed& seed(const SubControlVolumeFace& scvf) const
@@ -74,7 +78,8 @@ public:
private:
- void initializeSeeds_(std::vector<bool>& boundaryVertices)
+ void initializeSeeds_(std::vector<bool>& boundaryVertices,
+ std::vector<bool>& isFaceHandled)
{
const auto numScvf = problem_().model().globalFvGeometry().numScvf();
const auto numBoundaryScvf = problem_().model().globalFvGeometry().numBoundaryScvf();
@@ -92,7 +97,7 @@ private:
for (const auto& scvf : scvfs(fvGeometry))
{
// skip the rest if we already handled this face
- if (scvfIndexMap_[scvf.index()] != -1)
+ if (isFaceHandled[scvf.index()])
continue;
if (boundaryVertices[scvf.vertexIndex()])
@@ -102,7 +107,10 @@ private:
// update the index map entries for the global scv faces in the interaction volume
for (auto scvfIdxGlobal : boundarySeeds_.back().globalScvfIndices())
+ {
scvfIndexMap_[scvfIdxGlobal] = boundarySeedIndex;
+ isFaceHandled[scvfIdxGlobal] = true;
+ }
// increment counter
boundarySeedIndex++;
@@ -114,7 +122,10 @@ private:
// update the index map entries for the global scv faces in the interaction volume
for (auto scvfIdxGlobal : seeds_.back().globalScvfIndices())
+ {
scvfIndexMap_[scvfIdxGlobal] = seedIndex;
+ isFaceHandled[scvfIdxGlobal] = true;
+ }
// increment counter
seedIndex++;
@@ -130,6 +141,12 @@ private:
const Problem& problem_() const
{ return *problemPtr_; }
+ const Implementation& asImp_() const
+ { return *static_cast<Implementation*>(this); }
+
+ Implementation& asImp_()
+ { return *static_cast<Implementation*>(this); }
+
const Problem* problemPtr_;
GridView gridView_;
std::vector<IndexType> scvfIndexMap_;
diff --git a/dumux/discretization/cellcentered/mpfa/helper.hh b/dumux/discretization/cellcentered/mpfa/helper.hh
index 164ab6d0e642105c311e3fec2ce976a926941e98..a2f04306595427d2c25b550bfac154d4a18974dc 100644
--- a/dumux/discretization/cellcentered/mpfa/helper.hh
+++ b/dumux/discretization/cellcentered/mpfa/helper.hh
@@ -29,19 +29,17 @@
namespace Dumux
{
-// Mpfa method-specific implementation of the helper class
-template<class TypeTag, MpfaMethods Method, int dim>
-class MpfaHelperBase {};
+// Mpfa method-specific implementation of the helper class (again dimension-dependent)
+template<class TypeTag, MpfaMethods Method, int dim, int dimWorld>
+class MpfaMethodHelper;
// dimension-specific implementation of the helper class (common for all methods)
-template<class TypeTag, int dim>
-class MpfaHelperImplementation {};
+template<class TypeTag, int dim, int dimWorld>
+class MpfaDimensionHelper;
-// Specialization for dim == 2
+// Specialization for dim == 2, dimWorld == 2
template<class TypeTag>
-class MpfaHelperImplementation<TypeTag, 2> : public MpfaHelperBase<TypeTag,
- GET_PROP_VALUE(TypeTag, MpfaMethod),
- 2>
+class MpfaDimensionHelper<TypeTag, /*dim*/2, /*dimWorld*/2>
{
using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
@@ -51,15 +49,19 @@ class MpfaHelperImplementation<TypeTag, 2> : public MpfaHelperBase<TypeTag,
using LocalIndexType = typename InteractionVolume::LocalIndexType;
+ // We know that dim = 2 and dimworld = 2
static const int dim = GridView::dimension;
static const int dimWorld = GridView::dimensionworld;
- using DimVector = Dune::FieldVector<Scalar, dimWorld>;
+ using GlobalPosition = Dune::FieldVector<Scalar, dimWorld>;
using ScvfVector = std::array<const SubControlVolumeFace*, dim>;
- using LocalBasis = std::array<DimVector, dim>;
+ using LocalBasis = std::array<GlobalPosition, dim>;
+
+ using PointVector = std::vector<GlobalPosition>;
+ using FaceReferenceElements = typename Dune::ReferenceElements<Scalar, dim-1>;
public:
- // Gets the common scv face in the outer element and the other scv face sharing the same vertex
+ // Gets the two scv faces in the outer element, that share the vertex
// orders them to form a right hand system, local indices can be deduced from on the rotational direction
static ScvfVector getCommonAndNextScvFace(const SubControlVolumeFace& outsideScvf,
const FVElementGeometry& fvGeometry,
@@ -67,19 +69,29 @@ public:
{
LocalIndexType commonFaceIdx = clockWise ? 0 : 1;
LocalIndexType nextFaceIdx = clockWise ? 1 : 0;
+ auto vIdxGlobal = outsideScvf.vertexIndex();
+ unsigned int count = 0;
ScvfVector scvfVector;
for (const auto& scvf : scvfs(fvGeometry))
{
- if (outsideScvf.vertexIndex() == scvf.vertexIndex())
+ if (scvf.vertexIndex() == vIdxGlobal)
{
if (scvf.outsideScvIdx() == outsideScvf.insideScvIdx())
+ {
scvfVector[commonFaceIdx] = &scvf;
+ count++;
+ }
else
+ {
scvfVector[nextFaceIdx] = &scvf;
+ count++;
+ }
}
}
+ // make sure we found two faces
+ assert(count == 2 && "did not find two scv faces sharing the vertex in the outside element");
return scvfVector;
}
@@ -98,32 +110,74 @@ public:
return innerNormals;
}
- // calculates the determinant of the local basis
+ // the determinant of the local basis is equal to the cross product for dim = dimWorld = 2
static Scalar calculateDetX(const LocalBasis& localBasis)
{
- static const Dune::FieldMatrix<Scalar, dim, dim> R = {{0.0, 1.0}, {-1.0, 0.0}};
// make sure the basis forms a right hand system
assert(isRightHandSystem(localBasis) > 0 && "Local basis does not form a right hand system");
+ return crossProduct<Scalar>(localBasis[0], localBasis[1]);
+ }
- DimVector tmp(0.0);
- R.mv(localBasis[1], tmp);
+ // returns the global number of scvfs in the grid
+ static std::size_t getGlobalNumScvf(const GridView& gridView)
+ {
+ Dune::GeometryType triangle, quadrilateral;
+ triangle.makeTriangle();
+ quadrilateral.makeQuadrilateral();
- return tmp*localBasis[0];
+ return gridView.size(triangle)*6 + gridView.size(quadrilateral)*8;
}
+ //! Check whether or not the local basis forms a right hand system
static bool isRightHandSystem(const LocalBasis& localBasis)
+ { return !std::signbit(crossProduct<Scalar>(localBasis[0], localBasis[1])); }
+
+ //! get sub control volume face corners on a given face geometry for the given local index
+ static PointVector getScvfCorners(const PointVector& isCorners, unsigned int cornerIdx)
{
- if (std::signbit(crossProduct<Scalar>(localBasis[0], localBasis[1])))
- return false;
- return true;
+ // compute intersection center (in 2d intersections have two corners)
+ GlobalPosition center = isCorners[0] + isCorners[1];
+ center /= 2;
+
+ if (cornerIdx == 0)
+ return PointVector({center, isCorners[0]});
+ else if (cornerIdx == 1)
+ return PointVector({center, isCorners[1]});
+ else
+ DUNE_THROW(Dune::InvalidStateException, "local index exceeds the number of corners of 2d intersections");
+ }
+
+ //! calculate integration point on an scvf
+ static GlobalPosition getScvfIntegrationPoint(const PointVector& scvfCorners, Scalar q)
+ {
+ // scvfs in 3d are always quadrilaterals
+ // ordering -> first corner: facet center, last corner: vertex
+ if (q == 0.0)
+ return scvfCorners[0];
+
+ auto d = scvfCorners[1] - scvfCorners[0];
+ d *= q;
+ return scvfCorners[0] + d;
+ }
+
+ //! calculate the area of a scvf
+ static Scalar getScvfArea(const PointVector& scvfCorners)
+ { return (scvfCorners[1]-scvfCorners[0]).two_norm(); }
+
+ static std::size_t getNumLocalScvfs(const Dune::GeometryType gt)
+ {
+ if (gt == Dune::GeometryType(Dune::GeometryType::simplex, 2))
+ return 6;
+ else if (gt == Dune::GeometryType(Dune::GeometryType::cube, 2))
+ return 8;
+ else
+ DUNE_THROW(Dune::InvalidStateException, "unknown 2d geometry type " << gt);
}
};
-// Specialization for dim == 3
+// Specialization for dim == 2, dimWorld == 3
template<class TypeTag>
-class MpfaHelperImplementation<TypeTag, 3> : public MpfaHelperBase<TypeTag,
- GET_PROP_VALUE(TypeTag, MpfaMethod),
- 3>
+class MpfaDimensionHelper<TypeTag, /*dim*/2, /*dimWorld*/3>
{
using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
@@ -133,12 +187,91 @@ class MpfaHelperImplementation<TypeTag, 3> : public MpfaHelperBase<TypeTag,
using LocalIndexType = typename InteractionVolume::LocalIndexType;
+ // We know that dim = 2 and dimworld = 3
static const int dim = GridView::dimension;
static const int dimWorld = GridView::dimensionworld;
- using DimVector = Dune::FieldVector<Scalar, dimWorld>;
+ using GlobalPosition = Dune::FieldVector<Scalar, dimWorld>;
using ScvfVector = std::array<const SubControlVolumeFace*, dim>;
- using LocalBasis = std::array<DimVector, dim>;
+ using LocalBasis = std::array<GlobalPosition, dim>;
+ using PointVector = std::vector<GlobalPosition>;
+
+public:
+ // calculates the inner normal vectors
+ static LocalBasis calculateInnerNormals(const LocalBasis& localBasis)
+ {
+ // make sure the basis forms a right hand system
+ assert(isRightHandSystem(localBasis) > 0 && "Local basis does not form a right hand system");
+
+ // compute vector normal to the basis plane
+ auto normal = crossProduct<Scalar>(localBasis[0], localBasis[1]);
+ normal /= normal.two_norm();
+
+ // compute inner normals using the normal vector
+ LocalBasis innerNormals;
+ innerNormals[0] = crossProduct<Scalar>(localBasis[1], normal);
+ innerNormals[1] = crossProduct<Scalar>(normal, localBasis[0]);
+
+ return innerNormals;
+ }
+
+ // This is actually not the determinant of the basis for dim = 2 < dimWorld = 3
+ // It simply is the area of the parallelogram spanned by the basis vectors
+ static Scalar calculateDetX(const LocalBasis& localBasis)
+ { return std::abs( crossProduct<Scalar>(localBasis[0], localBasis[1]).two_norm() ); }
+
+ // returns the global number of scvfs in the grid
+ static std::size_t getGlobalNumScvf(const GridView& gridView)
+ {
+ Dune::GeometryType triangle, quadrilateral;
+ triangle.makeTriangle();
+ quadrilateral.makeQuadrilateral();
+
+ return gridView.size(triangle)*6 + gridView.size(quadrilateral)*8;
+ }
+
+ //! For 2d in 3d there is no unique basis forming a right hand system
+ static bool isRightHandSystem(const LocalBasis& localBasis)
+ { return true; }
+
+ //! get sub control volume face corners on a given face geometry for the given local index
+ static PointVector getScvfCorners(const PointVector& isCorners, unsigned int cornerIdx)
+ { return MpfaDimensionHelper<TypeTag, dim, dim>::getScvfCorners(isCorners, cornerIdx); }
+
+ //! calculate integration point on an scvf
+ static GlobalPosition getScvfIntegrationPoint(const PointVector& scvfCorners, Scalar q)
+ { return MpfaDimensionHelper<TypeTag, dim, dim>::getScvfIntegrationPoint(scvfCorners, q); }
+
+ //! calculate the area of a scvf
+ static Scalar getScvfArea(const PointVector& scvfCorners)
+ { return MpfaDimensionHelper<TypeTag, dim, dim>::getScvfArea(scvfCorners); }
+
+ static std::size_t getNumLocalScvfs(const Dune::GeometryType gt)
+ { return MpfaDimensionHelper<TypeTag, dim, dim>::getNumLocalScvfs(gt); }
+};
+
+// Specialization for dim == 3 (dimWorld has to be 3)
+template<class TypeTag>
+class MpfaDimensionHelper<TypeTag, /*dim*/3, /*dimWorld*/3>
+{
+ using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+ 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 LocalIndexType = typename InteractionVolume::LocalIndexType;
+
+ // We know that dim = 3 and dimworld = 3
+ static const int dim = GridView::dimension;
+ static const int dimWorld = GridView::dimensionworld;
+
+ using GlobalPosition = Dune::FieldVector<Scalar, dimWorld>;
+ using ScvfVector = std::array<const SubControlVolumeFace*, dim>;
+ using LocalBasis = std::array<GlobalPosition, dim>;
+
+ using PointVector = std::vector<GlobalPosition>;
+ using FaceReferenceElements = typename Dune::ReferenceElements<Scalar, dim-1>;
public:
@@ -161,11 +294,134 @@ public:
return tripleProduct<Scalar>(localBasis[0], localBasis[1], localBasis[2]);
}
+ // returns the global number of scvfs in the grid
+ static std::size_t getGlobalNumScvf(const GridView& gridView)
+ {
+ 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;
+ }
+
+ //! Check whether or not the local basis forms a right hand system
static bool isRightHandSystem(const LocalBasis& localBasis)
+ { return !std::signbit(tripleProduct<Scalar>(localBasis[0], localBasis[1], localBasis[2])); }
+
+ //! get sub control volume face corners on a given face geometry for the given local index
+ static PointVector getScvfCorners(const PointVector& isCorners, unsigned int cornerIdx)
+ {
+ // maximum number of necessary points is 9 (for quadrilateral)
+ GlobalPosition p[9];
+ auto numCorners = isCorners.size();
+
+ // fill point vector with center and corners
+ p[0] = 0.0;
+ for (int i = 0; i < numCorners; ++i)
+ {
+ p[0] += isCorners[i];
+ p[i+1] = isCorners[i];
+ }
+ p[0] /= numCorners;
+
+ // proceed according to number of corners
+ switch (numCorners)
+ {
+ case 3: // triangle
+ {
+ // add edge midpoints, there are 3 for triangles
+ p[numCorners+1] = isCorners[1] + isCorners[0];
+ p[numCorners+1] /= 2;
+ p[numCorners+2] = isCorners[2] + isCorners[0];
+ p[numCorners+2] /= 2;
+ p[numCorners+3] = isCorners[2] + isCorners[1];
+ p[numCorners+3] /= 2;
+
+ //! Only build the maps the first time we encounter a triangle
+ static const std::uint8_t vo = 1; //! vertex offset in point vector p
+ static const std::uint8_t eo = 4; //! edge offset in point vector p
+ static const std::uint8_t map[3][4] =
+ {
+ {0, eo+1, eo+0, vo+0},
+ {0, eo+0, eo+2, vo+1},
+ {0, eo+2, eo+1, vo+2}
+ };
+
+ return PointVector( {p[map[cornerIdx][0]],
+ p[map[cornerIdx][1]],
+ p[map[cornerIdx][2]],
+ p[map[cornerIdx][3]]} );
+ }
+ case 4: // quadrilateral
+ {
+ // add edge midpoints, there are 4 for quadrilaterals
+ p[numCorners+1] = isCorners[2] + isCorners[0];
+ p[numCorners+1] /= 2;
+ p[numCorners+2] = isCorners[3] + isCorners[1];
+ p[numCorners+2] /= 2;
+ p[numCorners+3] = isCorners[1] + isCorners[0];
+ p[numCorners+3] /= 2;
+ p[numCorners+4] = isCorners[3] + isCorners[2];
+ p[numCorners+4] /= 2;
+
+ //! Only build the maps the first time we encounter a quadrilateral
+ static const std::uint8_t vo = 1; //! vertex offset in point vector p
+ static const std::uint8_t eo = 5; //! face offset in point vector p
+ static const std::uint8_t map[4][4] =
+ {
+ {0, eo+0, eo+2, vo+0},
+ {0, eo+2, eo+1, vo+1},
+ {0, eo+3, eo+0, vo+2},
+ {0, eo+1, eo+3, vo+3}
+ };
+
+ return PointVector( {p[map[cornerIdx][0]],
+ p[map[cornerIdx][1]],
+ p[map[cornerIdx][2]],
+ p[map[cornerIdx][3]]} );
+ }
+ default:
+ DUNE_THROW(Dune::NotImplemented, "Mpfa scvf corners for dim=" << dim
+ << " dimWorld=" << dimWorld
+ << " corners=" << numCorners);
+ }
+ }
+
+ //! calculate integration point on a scvf
+ static GlobalPosition getScvfIntegrationPoint(const PointVector& scvfCorners, Scalar q)
+ {
+ // scvfs in 3d are always quadrilaterals
+ // ordering -> first corner: facet center, last corner: vertex
+ if (q == 0.0)
+ return scvfCorners[0];
+
+ auto d = scvfCorners[3] - scvfCorners[0];
+ d *= q;
+ return scvfCorners[0] + d;
+ }
+
+ //! calculate scvf area
+ static Scalar getScvfArea(const PointVector& scvfCorners)
+ {
+ // after Wolfram alpha quadrilateral area
+ return 0.5*Dumux::crossProduct(scvfCorners[3]-scvfCorners[0], scvfCorners[2]-scvfCorners[1]).two_norm();
+ }
+
+ //! determine the number of scvfs of an element
+ static std::size_t getNumLocalScvfs(const Dune::GeometryType gt)
{
- if (std::signbit(tripleProduct<Scalar>(localBasis[0], localBasis[1], localBasis[2])))
- return false;
- return true;
+ if (gt == Dune::GeometryType(Dune::GeometryType::simplex, 3))
+ return 12;
+ else if (gt == Dune::GeometryType(Dune::GeometryType::pyramid, 3))
+ return 16;
+ else if (gt == Dune::GeometryType(Dune::GeometryType::prism, 3))
+ return 18;
+ else if (gt == Dune::GeometryType(Dune::GeometryType::cube, 3))
+ return 24;
+ else
+ DUNE_THROW(Dune::InvalidStateException, "unknown 3d geometry type " << gt);
}
};
@@ -174,11 +430,11 @@ public:
* \brief Helper class to get the required information on an interaction volume.
* Specializations depending on the method and dimension are provided.
*/
-template<class TypeTag>
-class MpfaHelper : public MpfaHelperImplementation<TypeTag, GET_PROP_TYPE(TypeTag, GridView)::dimension>
+template<class TypeTag, MpfaMethods Method, int dim, int dimWorld>
+class CCMpfaHelperImplementation : public MpfaDimensionHelper<TypeTag, dim, dimWorld>,
+ public MpfaMethodHelper<TypeTag, Method, dim, dimWorld>
{
- using ParentType = MpfaHelperImplementation<TypeTag, GET_PROP_TYPE(TypeTag, GridView)::dimension>;
-
+ using Implementation = typename GET_PROP_TYPE(TypeTag, MpfaHelper);
using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
@@ -190,15 +446,15 @@ class MpfaHelper : public MpfaHelperImplementation<TypeTag, GET_PROP_TYPE(TypeTa
using GlobalIndexType = 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 GlobalPosition = Dune::FieldVector<Scalar, dimWorld>;
using ScvfVector = std::array<const SubControlVolumeFace*, dim>;
- using LocalBasis = std::array<DimVector, dim>;
+ using LocalBasis = std::array<GlobalPosition, dim>;
+
+ using CoordScalar = typename GridView::ctype;
+ using ReferenceElements = typename Dune::ReferenceElements<CoordScalar, dim>;
public:
- // returns shared pointers to the two scv faces that share a vertex in the order of a right hand system
+ // returns shared pointers to the scv faces that share a vertex in the order of a right hand system
static ScvfVector getScvFacesAtVertex(const GlobalIndexType vIdxGlobal,
const Element& element,
const FVElementGeometry& fvGeometry)
@@ -225,7 +481,7 @@ public:
assert(count == dim);
// sort the scv faces to form a right hand system
- if (!ParentType::isRightHandSystem(basisVectors))
+ if (!Implementation::isRightHandSystem(basisVectors))
std::swap(scvfVector[0], scvfVector[1]);
return scvfVector;
@@ -235,8 +491,20 @@ public:
static LocalIndexType getCommonFaceLocalIndex(const SubControlVolumeFace& outsideScvf,
const ScvfVector& insideScvFaces)
{
+ auto outsideScvIdx = outsideScvf.insideScvIdx();
for (int i = 0; i < insideScvFaces.size(); i++)
- if (insideScvFaces[i]->outsideScvIdx() == outsideScvf.insideScvIdx())
+ if (contains(insideScvFaces[i]->outsideScvIndices(), outsideScvIdx))
+ return i;
+
+ DUNE_THROW(Dune::InvalidStateException, "Could not find corresponding scvf in the provided vector of scvfs.");
+ }
+
+ // Finds the local index in an ScvfVector that corresponds to a given scvf
+ static LocalIndexType getLocalFaceIndex(const SubControlVolumeFace& scvf,
+ const ScvfVector& scvFaces)
+ {
+ for (int i = 0; i < scvFaces.size(); i++)
+ if (scvFaces[i]->index() == scvf.index())
return i;
DUNE_THROW(Dune::InvalidStateException, "Could not find corresponding scvf in the provided vector of scvfs.");
@@ -264,10 +532,66 @@ public:
DUNE_THROW(Dune::InvalidStateException, "unknown boundary condition type");
}
+
+ // returns a vector, which maps a bool (true if ghost vertex) to each vertex index
+ static std::vector<bool> findGhostVertices(const Problem& problem, const GridView& gridView)
+ {
+ std::vector<bool> ghostVertices(gridView.size(dim), false);
+
+ // if not run in parallel, skip the rest
+ if (Dune::MPIHelper::getCollectiveCommunication().size() == 1)
+ return ghostVertices;
+
+ // mpfa methods can not yet handle ghost cells
+ if (gridView.ghostSize(0) > 0)
+ DUNE_THROW(Dune::InvalidStateException, "Mpfa methods in parallel do not work with ghost cells. Use overlap cells instead.");
+
+ // mpfa methods have to have overlapping cells
+ if (gridView.overlapSize(0) == 0)
+ DUNE_THROW(Dune::InvalidStateException, "Grid no overlaping cells. This is required by mpfa methods in parallel.");
+
+ for (const auto& element : elements(gridView))
+ {
+ for (const auto& is : intersections(gridView, element))
+ {
+ if (!is.neighbor() && !is.boundary())
+ {
+ const auto& refElement = ReferenceElements::general(element.geometry().type());
+ for (unsigned int isVertex = 0; isVertex < is.geometry().corners(); ++isVertex)
+ {
+ const auto vIdxLocal = refElement.subEntity(is.indexInInside(), 1, isVertex, dim);
+ const auto vIdxGlobal = problem.vertexMapper().subIndex(element, vIdxLocal, dim);
+
+ ghostVertices[vIdxGlobal] = true;
+ }
+ }
+ }
+ }
+
+ return ghostVertices;
+ }
+
+ //! returns whether or not a value exists in a vector
+ template<typename V1, typename V2>
+ static bool contains(const std::vector<V1>& vector, const V2 value)
+ { return std::find(vector.begin(), vector.end(), value) != vector.end(); }
};
+
+/*!
+ * \ingroup Mpfa
+ * \brief Helper class for the mpfa methods for the construction of the interaction regions etc.
+ * It inherits from dimension-, dimensionworld- and method-specific implementations.
+ */
+template<class TypeTag>
+using CCMpfaHelper = CCMpfaHelperImplementation<TypeTag,
+ GET_PROP_VALUE(TypeTag, MpfaMethod),
+ GET_PROP_TYPE(TypeTag, GridView)::dimension,
+ GET_PROP_TYPE(TypeTag, GridView)::dimensionworld>;
+
} // end namespace
// The implemented helper classes need to be included here
+#include <dumux/discretization/cellcentered/mpfa/lmethod/helper.hh>
#include <dumux/discretization/cellcentered/mpfa/omethod/helper.hh>
#include <dumux/discretization/cellcentered/mpfa/omethodfps/helper.hh>
diff --git a/dumux/discretization/cellcentered/mpfa/hybridfps/globalfvgeometry.hh b/dumux/discretization/cellcentered/mpfa/hybridfps/globalfvgeometry.hh
new file mode 100644
index 0000000000000000000000000000000000000000..ba4526ae8fb618bc0c8db48f60f5a0cab4ca99b7
--- /dev/null
+++ b/dumux/discretization/cellcentered/mpfa/hybridfps/globalfvgeometry.hh
@@ -0,0 +1,281 @@
+// -*- 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_GLOBALFVGEOMETRY_HH
+#define DUMUX_DISCRETIZATION_CC_MPFA_O_HYBRIDFPS_GLOBALFVGEOMETRY_HH
+
+#include <dune/geometry/multilineargeometry.hh>
+#include <dune/geometry/referenceelements.hh>
+
+#include <dumux/implicit/cellcentered/mpfa/properties.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 CCMpfaOHybridFpsGlobalFVGeometry
+{};
+
+// specialization in case the FVElementGeometries are stored
+template<class TypeTag>
+class CCMpfaOHybridFpsGlobalFVGeometry<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
+ CCMpfaOHybridFpsGlobalFVGeometry(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_; }
+
+ // 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 = gridView_.size(0);
+ 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 CCMpfaOHybridFpsGlobalFVGeometry& 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
diff --git a/dumux/discretization/cellcentered/mpfa/hybridfps/globalinteractionvolumeseeds.hh b/dumux/discretization/cellcentered/mpfa/hybridfps/globalinteractionvolumeseeds.hh
new file mode 100644
index 0000000000000000000000000000000000000000..f14ac446c7739e0168a90c273208ab8cdd0d9bee
--- /dev/null
+++ b/dumux/discretization/cellcentered/mpfa/hybridfps/globalinteractionvolumeseeds.hh
@@ -0,0 +1,158 @@
+// -*- 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/implicit/cellcentered/mpfa/properties.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().globalFvGeometry().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().globalFvGeometry());
+ 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().globalFvGeometry());
+ 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
diff --git a/dumux/discretization/cellcentered/mpfa/hybridfps/helper.hh b/dumux/discretization/cellcentered/mpfa/hybridfps/helper.hh
new file mode 100644
index 0000000000000000000000000000000000000000..695fc18d6a33da953caafbf4e67284e42d4a895c
--- /dev/null
+++ b/dumux/discretization/cellcentered/mpfa/hybridfps/helper.hh
@@ -0,0 +1,125 @@
+// -*- 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
diff --git a/dumux/discretization/cellcentered/mpfa/hybridfps/interactionvolume.hh b/dumux/discretization/cellcentered/mpfa/hybridfps/interactionvolume.hh
new file mode 100644
index 0000000000000000000000000000000000000000..6e0b3242680b82ee10ae501de918ba080159786a
--- /dev/null
+++ b/dumux/discretization/cellcentered/mpfa/hybridfps/interactionvolume.hh
@@ -0,0 +1,72 @@
+// -*- 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 classes for interaction volume of mpfa methods.
+ */
+#ifndef DUMUX_DISCRETIZATION_CC_MPFA_O_HYBRID_FPS_INTERACTIONVOLUME_HH
+#define DUMUX_DISCRETIZATION_CC_MPFA_O_HYBRID_FPS_INTERACTIONVOLUME_HH
+
+#include <dumux/common/math.hh>
+
+#include <dune/localfunctions/lagrange/pqkfactory.hh>
+#include <dumux/implicit/cellcentered/mpfa/properties.hh>
+
+#include <dumux/discretization/cellcentered/mpfa/omethodfps/interactionvolume.hh>
+
+namespace Dumux
+{
+//! Specialization of the interaction volume traits class
+template<class TypeTag>
+class CCMpfaOHybridFpsInteractionVolumeTraits : public CCMpfaOInteractionVolumeTraits<TypeTag>
+{
+public:
+ // we will use the o method's interaction volume where necessary
+ using BoundaryInteractionVolume = CCMpfaInteractionVolumeImplementation<TypeTag, MpfaMethods::oMethodFps>;
+};
+
+/*!
+ * \ingroup Mpfa
+ * \brief Base class for the interaction volumes of the mpfa-o method with full pressure support.
+ */
+template<class TypeTag>
+class CCMpfaInteractionVolumeImplementation<TypeTag, MpfaMethods::oMethodFpsHybrid> : public CCMpfaOInteractionVolume<TypeTag, CCMpfaOHybridFpsInteractionVolumeTraits<TypeTag>>
+{
+ using Traits = CCMpfaOHybridFpsInteractionVolumeTraits<TypeTag>;
+ using ParentType = CCMpfaOInteractionVolume<TypeTag, Traits>;
+
+ using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
+ using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
+ using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);
+
+ using IVSeed = typename Traits::Seed;
+
+public:
+
+ CCMpfaInteractionVolumeImplementation(const IVSeed& seed,
+ const Problem& problem,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars)
+ : ParentType(seed, problem, fvGeometry, elemVolVars)
+ {}
+};
+
+} // end namespace
+
+#endif
diff --git a/dumux/discretization/cellcentered/mpfa/interactionvolume.hh b/dumux/discretization/cellcentered/mpfa/interactionvolume.hh
index 2ef7555a9a0a2e083a4f412764816120200b46c6..8d8f233b95d10045126e4b55e0c4542659f85c19 100644
--- a/dumux/discretization/cellcentered/mpfa/interactionvolume.hh
+++ b/dumux/discretization/cellcentered/mpfa/interactionvolume.hh
@@ -42,6 +42,7 @@ using CCMpfaInteractionVolume = CCMpfaInteractionVolumeImplementation<TypeTag, G
} // end namespace
// the specializations of this class for the available methods have to be included here
+#include <dumux/discretization/cellcentered/mpfa/lmethod/interactionvolume.hh>
#include <dumux/discretization/cellcentered/mpfa/omethod/interactionvolume.hh>
#include <dumux/discretization/cellcentered/mpfa/omethodfps/interactionvolume.hh>
diff --git a/dumux/discretization/cellcentered/mpfa/interactionvolumebase.hh b/dumux/discretization/cellcentered/mpfa/interactionvolumebase.hh
index 4af5b4591cd409dc570cb4cd90fe137ab4c90f55..aa02689f9d739bc5efa6fb7cd23871dc35f34f63 100644
--- a/dumux/discretization/cellcentered/mpfa/interactionvolumebase.hh
+++ b/dumux/discretization/cellcentered/mpfa/interactionvolumebase.hh
@@ -27,6 +27,26 @@
namespace Dumux
{
+//! Base class for the interaction volume traits
+template<class TypeTag>
+class CCMpfaInteractionVolumeTraitsBase
+{
+ using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+ using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+
+ static const int dim = GridView::dimension;
+ static const int dimWorld = GridView::dimensionworld;
+ using GlobalPosition = Dune::FieldVector<Scalar, dimWorld>;
+
+public:
+ using LocalIndexType = std::uint8_t;
+ using LocalIndexSet = std::vector<LocalIndexType>;
+ using GlobalIndexType = typename GridView::IndexSet::IndexType;
+ using GlobalIndexSet = std::vector<GlobalIndexType>;
+
+ using Tensor = Dune::FieldMatrix<Scalar, dim, dim>;
+};
+
/*!
* \ingroup Mpfa
* \brief Base class for the interaction volumes of mpfa methods.
@@ -43,13 +63,26 @@ public:
using BoundaryInteractionVolume = typename Traits::BoundaryInteractionVolume;
using LocalIndexType = typename Traits::LocalIndexType;
using LocalIndexSet = typename Traits::LocalIndexSet;
- using LocalIndexPair = typename Traits::LocalIndexPair;
using GlobalIndexType = typename Traits::GlobalIndexType;
using GlobalIndexSet = typename Traits::GlobalIndexSet;
using Vector = typename Traits::Vector;
using PositionVector = typename Traits::PositionVector;
using Seed = typename Traits::Seed;
+ struct LocalFaceData
+ {
+ LocalIndexType localScvfIndex;
+ LocalIndexType localScvIndex;
+ bool isOutside;
+
+ LocalFaceData(const LocalIndexType faceIndex,
+ const LocalIndexType scvIndex,
+ bool isOut)
+ : localScvfIndex(faceIndex),
+ localScvIndex(scvIndex),
+ isOutside(isOut) {}
+ };
+
//! solves the local equation system for the computation of the transmissibilities
template<typename GetTensorFunction>
void solveLocalSystem(const GetTensorFunction& getTensor)
@@ -68,15 +101,15 @@ public:
{ DUNE_THROW(Dune::NotImplemented, "Actual interaction volume implementation does not provide a globalScvfs() method."); }
//! returns the local index of an scvf in the IV and a boolean whether or not it is on the negative side of the local scvf (flux has to be inverted)
- LocalIndexPair getLocalIndexPair(const SubControlVolumeFace& scvf) const
- { DUNE_THROW(Dune::NotImplemented, "Actual interaction volume implementation does not provide a getLocalIndexPair() method."); }
+ LocalFaceData getLocalFaceData(const SubControlVolumeFace& scvf) const
+ { DUNE_THROW(Dune::NotImplemented, "Actual interaction volume implementation does not provide a getLocalFaceData() method."); }
//! returns the transmissibilities corresponding to a local scvf
- Vector getTransmissibilities(const LocalIndexPair& localIndexPair) const
+ Vector getTransmissibilities(const LocalFaceData& localFaceData) const
{ DUNE_THROW(Dune::NotImplemented, "Actual interaction volume implementation does not provide a getTransmissibilities() method."); }
//! returns the neumann flux corresponding to a local scvf
- Scalar getNeumannFlux(LocalIndexPair& localIndexPair) const
+ Scalar getNeumannFlux(const LocalFaceData& localFaceData) const
{ DUNE_THROW(Dune::NotImplemented, "Actual interaction volume implementation does not provide a getNeumannFlux() method."); }
//! returns the local index in a vector for a given global index
diff --git a/dumux/discretization/cellcentered/mpfa/lmethod/globalinteractionvolumeseeds.hh b/dumux/discretization/cellcentered/mpfa/lmethod/globalinteractionvolumeseeds.hh
new file mode 100644
index 0000000000000000000000000000000000000000..86778212839be40182a94e40aa358a1e2a97f14f
--- /dev/null
+++ b/dumux/discretization/cellcentered/mpfa/lmethod/globalinteractionvolumeseeds.hh
@@ -0,0 +1,131 @@
+// -*- 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 volume seeds of mpfa methods.
+ */
+#ifndef DUMUX_DISCRETIZATION_MPFA_L_GLOBALINTERACTIONVOLUMESEEDS_HH
+#define DUMUX_DISCRETIZATION_MPFA_L_GLOBALINTERACTIONVOLUMESEEDS_HH
+
+#include <dumux/discretization/cellcentered/mpfa/globalinteractionvolumeseedsbase.hh>
+#include <dumux/discretization/cellcentered/mpfa/methods.hh>
+
+namespace Dumux
+{
+/*!
+ * \ingroup Mpfa
+ * \brief Specialization of the class for the mpfa-l method.
+ */
+template<class TypeTag>
+class CCMpfaGlobalInteractionVolumeSeedsImplementation<TypeTag, MpfaMethods::lMethod> : public CCMpfaGlobalInteractionVolumeSeedsBase<TypeTag>
+{
+ using ParentType = CCMpfaGlobalInteractionVolumeSeedsBase<TypeTag>;
+ // the parent needs to be friend to access the private initialization method
+ friend ParentType;
+
+ using InteractionVolume = typename GET_PROP_TYPE(TypeTag, InteractionVolume);
+ using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
+ using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
+ using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
+ using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+ using Helper = typename GET_PROP_TYPE(TypeTag, MpfaHelper);
+
+ using IndexType = typename InteractionVolume::GlobalIndexSet::value_type;
+ using Element = typename GridView::template Codim<0>::Entity;
+
+public:
+ CCMpfaGlobalInteractionVolumeSeedsImplementation(const GridView gridView) : ParentType(gridView) {}
+
+private:
+ void initializeSeeds_(std::vector<bool>& boundaryVertices, std::vector<bool>& isFaceHandled)
+ {
+ const auto numScvf = this->problem_().model().globalFvGeometry().numScvf();
+ const auto numBoundaryScvf = this->problem_().model().globalFvGeometry().numBoundaryScvf();
+
+ // reserve memory
+ this->seeds_.reserve(numScvf - numBoundaryScvf);
+ this->boundarySeeds_.reserve(numBoundaryScvf);
+
+ IndexType boundarySeedIndex = 0;
+ IndexType seedIndex = 0;
+ for (const auto& element : elements(this->gridView_))
+ {
+ auto fvGeometry = localView(this->problem_().model().globalFvGeometry());
+ fvGeometry.bindElement(element);
+ for (const auto& scvf : scvfs(fvGeometry))
+ {
+ // skip the rest if we already handled this face
+ if (isFaceHandled[scvf.index()])
+ continue;
+
+ if (boundaryVertices[scvf.vertexIndex()])
+ {
+ // make the boundary interaction volume seed
+ this->boundarySeeds_.emplace_back(Helper::makeBoundaryInteractionVolumeSeed(this->problem_(), element, fvGeometry, scvf));
+
+ // update the index map entries for the global scv faces in the interaction volume
+ for (auto scvfIdxGlobal : this->boundarySeeds_.back().globalScvfIndices())
+ {
+ this->scvfIndexMap_[scvfIdxGlobal] = boundarySeedIndex;
+ isFaceHandled[scvfIdxGlobal] = true;
+ }
+
+ // increment counter
+ boundarySeedIndex++;
+ }
+ else
+ {
+ // make the inner interaction volume seed only if we are on highest level of all connected elements
+ if (isLocalMaxLevel(element, scvf))
+ {
+ this->seeds_.emplace_back(Helper::makeInnerInteractionVolumeSeed(this->problem_(), element, fvGeometry, scvf));
+
+ // update the index map entries for the global scv faces in the interaction volume
+ for (auto scvfIdxGlobal : this->seeds_.back().globalScvfIndices())
+ {
+ this->scvfIndexMap_[scvfIdxGlobal] = seedIndex;
+ isFaceHandled[scvfIdxGlobal] = true;
+ }
+
+ // increment counter
+ seedIndex++;
+ }
+ }
+ }
+ }
+
+ // shrink vectors to actual size
+ this->seeds_.shrink_to_fit();
+ this->boundarySeeds_.shrink_to_fit();
+ }
+
+ bool isLocalMaxLevel(const Element& element, const SubControlVolumeFace& scvf) const
+ {
+ auto inLevel = element.level();
+ for (auto outsideIdx : scvf.outsideScvIndices())
+ if (this->problem_().model().globalFvGeometry().element(outsideIdx).level() > inLevel)
+ return false;
+ return true;
+ }
+};
+
+} // end namespace
+
+
+#endif
diff --git a/dumux/discretization/cellcentered/mpfa/lmethod/helper.hh b/dumux/discretization/cellcentered/mpfa/lmethod/helper.hh
new file mode 100644
index 0000000000000000000000000000000000000000..34e5699952acf4ce3afb117137ef10935e84b679
--- /dev/null
+++ b/dumux/discretization/cellcentered/mpfa/lmethod/helper.hh
@@ -0,0 +1,183 @@
+// -*- 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_MPFA_L_HELPER_HH
+#define DUMUX_DISCRETIZATION_CC_MPFA_L_HELPER_HH
+
+#include <dumux/common/math.hh>
+#include <dumux/discretization/cellcentered/mpfa/facetypes.hh>
+#include <dumux/discretization/cellcentered/mpfa/methods.hh>
+
+#include "localsubcontrolentityseeds.hh"
+
+namespace Dumux
+{
+/*!
+ * \ingroup Mpfa
+ * \brief Helper class to get the required information on an interaction volume.
+ * Specialization for the Mpfa-L method in two dimensions.
+ */
+template<class TypeTag>
+class MpfaMethodHelper<TypeTag, MpfaMethods::lMethod, /*dim*/2, /*dimWorld*/2>
+{
+ using Implementation = typename GET_PROP_TYPE(TypeTag, MpfaHelper);
+
+ static const int dim = 2;
+ static const int dimWorld = 2;
+
+ // The mpfa-o helper class used to construct the boundary interaction volume seeds
+ using oMethodHelper = CCMpfaHelperImplementation<TypeTag, MpfaMethods::oMethod, 2, 2>;
+
+ 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 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);
+
+ using Element = typename GridView::template Codim<0>::Entity;
+
+ using InteractionVolumeSeed = typename InteractionVolume::Seed;
+ using ScvSeed = typename InteractionVolumeSeed::LocalScvSeed;
+ using OuterScvSeed = typename InteractionVolumeSeed::LocalOuterScvSeed;
+ using BoundaryInteractionVolumeSeed = typename BoundaryInteractionVolume::Seed;
+
+ using GlobalIndexSet = typename InteractionVolume::GlobalIndexSet;
+ using GlobalIndexType = typename InteractionVolume::GlobalIndexType;
+ using LocalIndexSet = typename InteractionVolume::LocalIndexSet;
+ using LocalIndexType = typename InteractionVolume::LocalIndexType;
+ using Matrix = typename InteractionVolume::Matrix;
+public:
+ static InteractionVolumeSeed makeInnerInteractionVolumeSeed(const Problem& problem,
+ const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const SubControlVolumeFace& scvf)
+ {
+ std::vector<ScvSeed> scvSeeds;
+ std::vector<OuterScvSeed> outerScvSeeds;
+ std::vector<GlobalIndexType> globalScvfIndices(2);
+
+ // we'll have maximal 2 ScvSeeds and 2 OuterScvSeeds
+ scvSeeds.reserve(2);
+ outerScvSeeds.reserve(2);
+
+ fillEntitySeeds_(scvSeeds, outerScvSeeds, globalScvfIndices, problem, element, fvGeometry, scvf);
+
+ // free unused memory
+ scvSeeds.shrink_to_fit();
+ outerScvSeeds.shrink_to_fit();
+
+ // return interaction volume seed
+ return InteractionVolumeSeed(std::move(scvSeeds), std::move(outerScvSeeds), std::move(globalScvfIndices));
+ }
+
+ static BoundaryInteractionVolumeSeed makeBoundaryInteractionVolumeSeed(const Problem& problem,
+ const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const SubControlVolumeFace& scvf)
+ { return oMethodHelper::makeBoundaryInteractionVolumeSeed(problem, element, fvGeometry, scvf); }
+
+ template<class InteractionRegion>
+ static LocalIndexType selectionCriterion(const InteractionRegion& I1,
+ const InteractionRegion& I2,
+ const Matrix& M1,
+ const Matrix& M2)
+ {
+ Scalar eps = 1e-10;
+ Scalar t11, t12, t21, t22;
+ t11 = M1[I1.contiFaceLocalIdx][0];
+ t21 = M2[I2.contiFaceLocalIdx][0];
+ if (I1.contiFaceLocalIdx == 0)
+ {
+ t12 = M1[I1.contiFaceLocalIdx][1];
+ t22 = M2[I2.contiFaceLocalIdx][2];
+ }
+ else
+ {
+ t12 = M1[I1.contiFaceLocalIdx][2];
+ t22 = M2[I2.contiFaceLocalIdx][1];
+ }
+
+ Scalar s1 = std::abs(t11-t12);
+ Scalar s2 = std::abs(t22-t21);
+
+ if (s1 < s2 + eps*s1)
+ return 0;
+ else
+ return 1;
+ }
+
+private:
+ static void fillEntitySeeds_(std::vector<ScvSeed>& scvSeeds,
+ std::vector<OuterScvSeed>& outerScvSeeds,
+ std::vector<GlobalIndexType>& globalScvfIndices,
+ const Problem& problem,
+ const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const SubControlVolumeFace& scvf)
+ {
+ // make the first scv seed, we know this element will NOT be on the lowest local level
+ auto scvfVector = Implementation::getScvFacesAtVertex(scvf.vertexIndex(), element, fvGeometry);
+ auto localScvfIdx = Implementation::getLocalFaceIndex(scvf, scvfVector);
+ auto globalScvIdx = scvf.insideScvIdx();
+ scvSeeds.emplace_back( GlobalIndexSet({scvfVector[0]->index(), scvfVector[1]->index()}),
+ globalScvIdx,
+ localScvfIdx );
+
+ // get the surrounding elements and "outside" data
+ LocalIndexType otherScvfIdx = localScvfIdx == 1 ? 0 : 1;
+ auto e2 = problem.model().globalFvGeometry().element(scvf.outsideScvIdx());
+ auto e3 = problem.model().globalFvGeometry().element(scvfVector[otherScvfIdx]->outsideScvIdx());
+ auto e2Geometry = localView(problem.model().globalFvGeometry());
+ auto e3Geometry = localView(problem.model().globalFvGeometry());
+ e2Geometry.bindElement(e2);
+ e3Geometry.bindElement(e3);
+ auto e2Scvfs = Implementation::getCommonAndNextScvFace(scvf, e2Geometry, /*clockwise?*/localScvfIdx == 1);
+ auto e3Scvfs = Implementation::getCommonAndNextScvFace(*scvfVector[otherScvfIdx], e3Geometry, /*clockwise?*/localScvfIdx == 0);
+
+ // we now know the two faces for which flux calculation will happen using this iv seed
+ globalScvfIndices[0] = scvf.index();
+ globalScvfIndices[1] = e2Scvfs[otherScvfIdx]->index();
+
+ // scv seed for e2, we know the local common scvf index will be otherScvfIdx in 2d
+ scvSeeds.emplace_back( GlobalIndexSet({e2Scvfs[0]->index(), e2Scvfs[1]->index()}),
+ e2Scvfs[0]->insideScvIdx(),
+ otherScvfIdx );
+
+ // Outer seed for e3, we know the local common scvf index will be localScvfIdx in 2d
+ const auto& f3 = *e3Scvfs[localScvfIdx];
+ outerScvSeeds.emplace_back(f3.insideScvIdx(), f3.index());
+
+ // Outer seed for outside of e2, we know the local scvf index here will be localScvfIdx in 2d
+ auto e4 = problem.model().globalFvGeometry().element(e2Scvfs[localScvfIdx]->outsideScvIdx());
+ auto e4Geometry = localView(problem.model().globalFvGeometry());
+ e4Geometry.bindElement(e4);
+ auto e4Scvfs = Implementation::getCommonAndNextScvFace(*e2Scvfs[localScvfIdx], e4Geometry, /*clockwise?*/localScvfIdx == 1);
+ const auto& f4 = *e4Scvfs[otherScvfIdx];
+ outerScvSeeds.emplace_back(f4.insideScvIdx(), f4.index());
+ }
+};
+
+} // end namespace
+
+#endif
diff --git a/dumux/discretization/cellcentered/mpfa/lmethod/interactionregions.hh b/dumux/discretization/cellcentered/mpfa/lmethod/interactionregions.hh
new file mode 100644
index 0000000000000000000000000000000000000000..2a9a4287fcee3d73ca46574cabfda23419ffd1f4
--- /dev/null
+++ b/dumux/discretization/cellcentered/mpfa/lmethod/interactionregions.hh
@@ -0,0 +1,129 @@
+// -*- 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 classes for interaction volume of mpfa methods.
+ */
+#ifndef DUMUX_DISCRETIZATION_CC_MPFA_L_INTERACTIONREGIONS_HH
+#define DUMUX_DISCRETIZATION_CC_MPFA_L_INTERACTIONREGIONS_HH
+
+namespace Dumux
+{
+/*!
+ * \ingroup Mpfa
+ * \brief Base class for the interaction regions of the mpfa-l method.
+ */
+template<class TypeTag>
+struct InteractionRegion
+{
+private:
+ using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+ 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 InteractionVolume = typename GET_PROP_TYPE(TypeTag, InteractionVolume);
+
+ static const int dim = GridView::dimension;
+ static const int dimWorld = GridView::dimensionworld;
+
+ using Element = typename GridView::template Codim<0>::Entity;
+ using LocalIndexType = typename InteractionVolume::LocalIndexType;
+ using GlobalIndexType = typename InteractionVolume::GlobalIndexType;
+ using GlobalPosition = Dune::FieldVector<Scalar, dimWorld>;
+
+public:
+ LocalIndexType contiFaceLocalIdx;
+ std::vector<GlobalIndexType> scvIndices;
+ std::vector<GlobalPosition> scvCenters;
+ std::vector<GlobalPosition> normal;
+ std::vector<GlobalPosition> nu;
+ std::vector<Scalar> detX;
+ std::vector<Element> elements;
+
+ std::array<GlobalIndexType, 2> globalScvfs;
+
+ // Constructor for dim == 2
+ template<class ScvSeedType, class OuterScvSeedType>
+ InteractionRegion(const Problem& problem,
+ const FVElementGeometry& fvGeometry,
+ const ScvSeedType& scvSeed,
+ const OuterScvSeedType& outerSeed1,
+ const OuterScvSeedType& outerSeed2,
+ const Element& element1,
+ const Element& element2,
+ const Element& element3)
+ {
+ contiFaceLocalIdx = scvSeed.contiFaceLocalIdx();
+ globalScvfs[0] = scvSeed.globalScvfIndices()[contiFaceLocalIdx];
+ globalScvfs[1] = contiFaceLocalIdx == 0 ? outerSeed1.globalScvfIndex() : outerSeed2.globalScvfIndex();
+
+ elements.resize(3);
+ elements[0] = element1;
+ elements[1] = element2;
+ elements[2] = element3;
+
+ // The participating sub control entities
+ auto&& scv1 = fvGeometry.scv(scvSeed.globalIndex());
+ auto&& scv2 = fvGeometry.scv(outerSeed1.globalIndex());
+ auto&& scv3 = fvGeometry.scv(outerSeed2.globalIndex());
+ auto&& scvf1 = fvGeometry.scvf(scvSeed.globalScvfIndices()[0]);
+ auto&& scvf2 = fvGeometry.scvf(scvSeed.globalScvfIndices()[1]);
+
+ // The necessary coordinates and normals
+ GlobalPosition v = scvf1.vertexCorner();
+ GlobalPosition f1 = scvf1.facetCorner();
+ GlobalPosition f2 = scvf2.facetCorner();
+
+ scvCenters.resize(3);
+ scvCenters[0] = scv1.center();
+ scvCenters[1] = scv2.center();
+ scvCenters[2] = scv3.center();
+
+ normal.resize(2);
+ normal[0] = scvf1.unitOuterNormal();
+ normal[0] *= scvf1.area();
+ normal[1] = scvf2.unitOuterNormal();
+ normal[1] *= scvf2.area();
+
+ // indices of the participating scvs
+ scvIndices.resize(3);
+ scvIndices[0] = scv1.index();
+ scvIndices[1] = scv2.index();
+ scvIndices[2] = scv3.index();
+
+ // calculate nus and detXs
+ static const Dune::FieldMatrix<Scalar, dim, dim> R = {{0.0, 1.0}, {-1.0, 0.0}};
+ nu.resize(7);
+ R.mv(f2-scvCenters[0], nu[0]);
+ R.mv(scvCenters[0]-f1, nu[1]);
+ R.mv(f1-scvCenters[1], nu[2]);
+ R.mv(scvCenters[1]-v, nu[3]);
+ R.mv(v-scvCenters[2], nu[4]);
+ R.mv(scvCenters[2]-f2, nu[5]);
+ R.mv(v-scvCenters[0], nu[6]);
+
+ detX.resize(3);
+ detX[0] = (f1-scvCenters[0])*nu[0];
+ detX[1] = (v-scvCenters[1])*nu[2];
+ detX[2] = (f2-scvCenters[2])*nu[4];
+ }
+};
+}; //end namespace
+
+#endif
\ No newline at end of file
diff --git a/dumux/discretization/cellcentered/mpfa/lmethod/interactionvolume.hh b/dumux/discretization/cellcentered/mpfa/lmethod/interactionvolume.hh
new file mode 100644
index 0000000000000000000000000000000000000000..e597c04bb7ab6cd77b3de86247a655a1685bb57f
--- /dev/null
+++ b/dumux/discretization/cellcentered/mpfa/lmethod/interactionvolume.hh
@@ -0,0 +1,364 @@
+// -*- 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 classes for interaction volume of mpfa methods.
+ */
+#ifndef DUMUX_DISCRETIZATION_CC_MPFA_L_INTERACTIONVOLUME_HH
+#define DUMUX_DISCRETIZATION_CC_MPFA_L_INTERACTIONVOLUME_HH
+
+#include <dumux/common/math.hh>
+#include <dumux/implicit/cellcentered/mpfa/properties.hh>
+#include <dumux/discretization/cellcentered/mpfa/interactionvolumebase.hh>
+#include <dumux/discretization/cellcentered/mpfa/omethod/interactionvolume.hh>
+#include <dumux/discretization/cellcentered/mpfa/facetypes.hh>
+#include <dumux/discretization/cellcentered/mpfa/methods.hh>
+
+#include "interactionvolumeseed.hh"
+#include "interactionregions.hh"
+
+namespace Dumux
+{
+//! Specialization of the interaction volume traits class for the mpfa-o method
+template<class TypeTag>
+class CCMpfaLInteractionVolumeTraits : public CCMpfaInteractionVolumeTraitsBase<TypeTag>
+{
+ using BaseTraits = CCMpfaInteractionVolumeTraitsBase<TypeTag>;
+
+ using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+ using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+
+ static const int dim = GridView::dimension;
+ static const int dimWorld = GridView::dimensionworld;
+ using GlobalPosition = Dune::FieldVector<Scalar, dimWorld>;
+
+public:
+ using BoundaryInteractionVolume = CCMpfaInteractionVolumeImplementation<TypeTag, MpfaMethods::oMethod>;
+
+ // in the interaction region there will be dim faces and dim+1 cell pressures
+ using PositionVector = std::vector<GlobalPosition>;
+ using Matrix = Dune::FieldMatrix<Scalar, dim, dim+1>;
+ using Vector = Dune::FieldVector<Scalar, dim+1>;
+
+ using typename BaseTraits::LocalIndexSet;
+ using typename BaseTraits::GlobalIndexSet;
+ using Seed = CCMpfaLInteractionVolumeSeed<GlobalIndexSet, LocalIndexSet>;
+};
+
+/*!
+ * \ingroup Mpfa
+ * \brief Base class for the interaction volumes of the mpfa-l method.
+ */
+template<class TypeTag>
+class CCMpfaInteractionVolumeImplementation<TypeTag, MpfaMethods::lMethod> : public CCMpfaInteractionVolumeBase<TypeTag, CCMpfaLInteractionVolumeTraits<TypeTag>>
+{
+ // The interaction volume implementation has to be friend
+ friend typename GET_PROP_TYPE(TypeTag, InteractionVolume);
+ using Traits = CCMpfaLInteractionVolumeTraits<TypeTag>;
+ using ParentType = CCMpfaInteractionVolumeBase<TypeTag, Traits>;
+
+ 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 MpfaHelper = typename GET_PROP_TYPE(TypeTag, MpfaHelper);
+ using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
+ using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
+ using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);
+
+ static const int dim = GridView::dimension;
+ static const int dimWorld = GridView::dimensionworld;
+ using Element = typename GridView::template Codim<0>::Entity;
+ using GlobalPosition = Dune::FieldVector<Scalar, dimWorld>;
+ using InteractionRegion = InteractionRegion<TypeTag>;
+
+ using Vector = typename Traits::Vector;
+ using Tensor = typename Traits::Tensor;
+
+public:
+ using Matrix = typename Traits::Matrix;
+ using typename ParentType::LocalIndexType;
+ using typename ParentType::LocalIndexSet;
+ using typename ParentType::LocalFaceData;
+ using typename ParentType::GlobalIndexType;
+ using typename ParentType::GlobalIndexSet;
+ using typename ParentType::PositionVector;
+ using typename ParentType::Seed;
+
+ using ScvSeedType = typename Seed::LocalScvSeed;
+ using OuterScvSeedType = typename Seed::LocalOuterScvSeed;
+
+ CCMpfaInteractionVolumeImplementation(const Seed& seed,
+ const Problem& problem,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars)
+ : problemPtr_(&problem),
+ fvGeometryPtr_(&fvGeometry),
+ elemVolVarsPtr_(&elemVolVars),
+ regionUnique_(seed.isUnique()),
+ systemSolved_(false)
+ {
+ // set up the possible interaction regions
+ setupInteractionRegions_(seed, fvGeometry);
+ }
+
+ template<typename GetTensorFunction>
+ void solveLocalSystem(const GetTensorFunction& getTensor)
+ {
+ if (regionUnique_)
+ {
+ const auto& ir = interactionRegions_[0];
+ T_ = assembleMatrix_(getTensor, ir);
+ postSolve_(ir);
+ }
+ else
+ {
+ std::array<Matrix, 2> M;
+ M[0] = assembleMatrix_(getTensor, interactionRegions_[0]);
+ M[1] = assembleMatrix_(getTensor, interactionRegions_[1]);
+ const auto id = MpfaHelper::selectionCriterion(interactionRegions_[0], interactionRegions_[1], M[0], M[1]);
+
+ const auto& ir = interactionRegions_[id];
+ T_ = M[id];
+ postSolve_(ir);
+ }
+ }
+
+ //! returns the local index pair. This holds the scvf's local index and a boolean whether or not the flux has to be inverted
+ LocalFaceData getLocalFaceData(const SubControlVolumeFace& scvf) const
+ {
+ assert(systemSolved_ && "Scvf Indices not set yet. You have to call solveLocalSystem() beforehand.");
+ assert(scvf.index() == globalScvfIndices_[0] || scvf.index() == globalScvfIndices_[1] && "The provided scvf is not the flux face of the interaction volume.");
+
+ if (globalScvfIndices_[0] == scvf.index())
+ return LocalFaceData(contiFaceLocalIdx_, /*dummy*/0, false);
+ else
+ return LocalFaceData(contiFaceLocalIdx_, /*dummy*/0, true);
+ }
+
+ //! returns the transmissibilities corresponding to the bound scvf
+ Vector getTransmissibilities(const LocalFaceData& localFaceData) const
+ {
+ assert(systemSolved_ && "Transmissibilities not calculated yet. You have to call solveLocalSystem() beforehand.");
+
+ auto tij = T_[localFaceData.localScvfIndex];
+ if (localFaceData.isOutside)
+ tij *= -1.0;
+ return tij;
+ }
+
+ const GlobalIndexSet& volVarsStencil() const
+ {
+ assert(systemSolved_ && "volVarsStencil not set yet. You have to call solveLocalSystem() beforehand.");
+ return volVarsStencil_;
+ }
+
+ const PositionVector& volVarsPositions() const
+ {
+ assert(systemSolved_ && "globalScvfs not set yet. You have to call solveLocalSystem() beforehand.");
+ return volVarsPositions_;
+ }
+
+ const GlobalIndexSet& globalScvfs() const
+ {
+ assert(systemSolved_ && "globalScvfs not set yet. You have to call solveLocalSystem() beforehand.");
+ return globalScvfIndices_;
+ }
+
+ const Matrix& matrix() const
+ { return T_; }
+
+private:
+ //! Assembles and solves the local equation system
+ //! Specialization for dim = 2
+ template<typename GetTensorFunction, int d = dim>
+ typename std::enable_if<d == 2, Matrix>::type
+ assembleMatrix_(const GetTensorFunction& getTensor, const InteractionRegion& ir)
+ {
+ // the elements the scvs live in
+ const auto& e1 = ir.elements[0];
+ const auto& e2 = ir.elements[1];
+ const auto& e3 = ir.elements[2];
+
+ // the corresponding scvs
+ const auto scv1 = fvGeometry_().scv(ir.scvIndices[0]);
+ const auto scv2 = fvGeometry_().scv(ir.scvIndices[1]);
+ const auto scv3 = fvGeometry_().scv(ir.scvIndices[2]);
+
+ // Get diffusion tensors in the three scvs
+ const auto T1 = getTensor(e1, elemVolVars_()[scv1], scv1);
+ const auto T2 = getTensor(e2, elemVolVars_()[scv2], scv2);
+ const auto T3 = getTensor(e3, elemVolVars_()[scv3], scv3);
+
+ // required omega factors
+ Scalar w111 = calculateOmega_(ir.normal[0], ir.nu[0], ir.detX[0], T1);
+ Scalar w112 = calculateOmega_(ir.normal[0], ir.nu[1], ir.detX[0], T1);
+ Scalar w123 = calculateOmega_(ir.normal[0], ir.nu[2], ir.detX[1], T2);
+ Scalar w124 = calculateOmega_(ir.normal[0], ir.nu[3], ir.detX[1], T2);
+
+ Scalar w211 = calculateOmega_(ir.normal[1], ir.nu[0], ir.detX[0], T1);
+ Scalar w212 = calculateOmega_(ir.normal[1], ir.nu[1], ir.detX[0], T1);
+ Scalar w235 = calculateOmega_(ir.normal[1], ir.nu[4], ir.detX[2], T3);
+ Scalar w236 = calculateOmega_(ir.normal[1], ir.nu[5], ir.detX[2], T3);
+
+ Scalar xi711 = calculateXi_(ir.nu[6], ir.nu[0], ir.detX[0]);
+ Scalar xi712 = calculateXi_(ir.nu[6], ir.nu[1], ir.detX[0]);
+
+ Dune::FieldMatrix<Scalar, dim, dim> C, A;
+ Dune::FieldMatrix<Scalar, dim, dim+1> B;
+
+ C[0][0] = -w111;
+ C[0][1] = -w112;
+ C[1][0] = -w211;
+ C[1][1] = -w212;
+
+ A[0][0] = w111 - w124 - w123*xi711;
+ A[0][1] = w112 - w123*xi712;
+ A[1][0] = w211 - w236*xi711;
+ A[1][1] = w212 - w235 - w236*xi712;
+
+ B[0][0] = w111 + w112 + w123*(1 - xi711 -xi712);
+ B[0][1] = -w123 - w124;
+ B[0][2] = 0.0;
+ B[1][0] = w211 + w212 + w236*(1 - xi711 - xi712);
+ B[1][1] = 0.0;
+ B[1][2] = -w235 - w236;
+
+ // T = CA^-1B + D
+ A.invert();
+ auto T = (A.leftmultiply(C)).rightmultiplyany(B);
+ T[0][0] += w111 + w112;
+ T[1][0] += w211 + w212;
+
+ return T;
+ }
+
+ //! sets up the interaction regions for later transmissibility matrix calculation
+ //! Specialization for dim = 2
+ template<int d = dim>
+ typename std::enable_if<d == 2>::type
+ setupInteractionRegions_(const Seed& seed, const FVElementGeometry& fvGeometry)
+ {
+ const auto& globalFvGeometry = problem_().model().globalFvGeometry();
+
+ if (regionUnique_)
+ {
+ interactionRegions_.reserve(1);
+ const auto& scvSeed1 = seed.scvSeed(0);
+ const auto& scvSeed2 = seed.outerScvSeed(0);
+ const auto& scvSeed3 = seed.outerScvSeed(1);
+ auto e1 = globalFvGeometry.element(scvSeed1.globalIndex());
+ auto e2 = globalFvGeometry.element(scvSeed2.globalIndex());
+ auto e3 = globalFvGeometry.element(scvSeed3.globalIndex());
+ interactionRegions_.emplace_back(problem_(), fvGeometry, scvSeed1, scvSeed2, scvSeed3, e1, e2, e3);
+
+ }
+ else
+ {
+ interactionRegions_.reserve(2);
+ const auto& scvSeed1 = seed.scvSeed(0);
+ const auto& scvSeed2 = seed.scvSeed(1);
+ const auto& outerScvSeed1 = seed.outerScvSeed(0);
+ const auto& outerScvSeed2 = seed.outerScvSeed(1);
+ auto e1 = globalFvGeometry.element(scvSeed1.globalIndex());
+ auto e2 = globalFvGeometry.element(scvSeed2.globalIndex());
+ auto e3 = globalFvGeometry.element(outerScvSeed1.globalIndex());
+ auto e4 = globalFvGeometry.element(outerScvSeed2.globalIndex());
+
+ // scvSeed1 is the one the seed construction began at
+ if (scvSeed1.contiFaceLocalIdx() == 0)
+ {
+ interactionRegions_.emplace_back(problem_(), fvGeometry, scvSeed1, OuterScvSeedType(scvSeed2), outerScvSeed1, e1, e2, e3);
+ interactionRegions_.emplace_back(problem_(), fvGeometry, scvSeed2, outerScvSeed2, OuterScvSeedType(scvSeed1), e2, e4, e1);
+ }
+ else
+ {
+ interactionRegions_.emplace_back(problem_(), fvGeometry, scvSeed1, outerScvSeed1, OuterScvSeedType(scvSeed2), e1, e3, e2);
+ interactionRegions_.emplace_back(problem_(), fvGeometry, scvSeed2, OuterScvSeedType(scvSeed1), outerScvSeed2, e2, e1, e4);
+ }
+ }
+ }
+
+ void postSolve_(const InteractionRegion& ir)
+ {
+ globalScvfIndices_.resize(2);
+ globalScvfIndices_[0] = ir.globalScvfs[0];
+ globalScvfIndices_[1] = ir.globalScvfs[1];
+ volVarsStencil_ = ir.scvIndices;
+ volVarsPositions_ = ir.scvCenters;
+ contiFaceLocalIdx_ = ir.contiFaceLocalIdx;
+ systemSolved_ = true;
+ }
+
+ Scalar calculateOmega_(const GlobalPosition& normal,
+ const GlobalPosition& nu,
+ const Scalar detX,
+ const Tensor& T) const
+ {
+ GlobalPosition tmp;
+ T.mv(nu, tmp);
+ return (tmp*normal)/detX;
+ }
+
+ Scalar calculateOmega_(const GlobalPosition& normal,
+ const GlobalPosition& nu,
+ const Scalar detX,
+ const Scalar t) const
+ { return (normal*nu)*t/detX; }
+
+ //! Specialization for dim = 2
+ template<int d = dim>
+ typename std::enable_if<d == 2, Scalar>::type
+ calculateXi_(const GlobalPosition& nu1,
+ const GlobalPosition& nu2,
+ const Scalar detX) const
+ { return crossProduct<Scalar>(nu1, nu2)/detX; }
+
+ const Seed& seed_() const
+ { return seedPtr_; }
+
+ const Problem& problem_() const
+ { return *problemPtr_; }
+
+ const FVElementGeometry& fvGeometry_() const
+ { return *fvGeometryPtr_; }
+
+ const ElementVolumeVariables& elemVolVars_() const
+ { return *elemVolVarsPtr_; }
+
+ const Seed* seedPtr_;
+ const Problem* problemPtr_;
+ const FVElementGeometry* fvGeometryPtr_;
+ const ElementVolumeVariables* elemVolVarsPtr_;
+
+ bool regionUnique_;
+ bool systemSolved_;
+
+ LocalIndexType contiFaceLocalIdx_;
+ GlobalIndexSet globalScvfIndices_;
+ GlobalIndexSet volVarsStencil_;
+ PositionVector volVarsPositions_;
+
+ std::vector<InteractionRegion> interactionRegions_;
+
+ Matrix T_;
+};
+
+} // end namespace
+
+#endif
diff --git a/dumux/discretization/cellcentered/mpfa/lmethod/interactionvolumeseed.hh b/dumux/discretization/cellcentered/mpfa/lmethod/interactionvolumeseed.hh
new file mode 100644
index 0000000000000000000000000000000000000000..808abbece40e3feaf9b6c02aca04a6fbb26eb59f
--- /dev/null
+++ b/dumux/discretization/cellcentered/mpfa/lmethod/interactionvolumeseed.hh
@@ -0,0 +1,98 @@
+// -*- 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 interaction volume seeds of mpfa methods.
+ */
+#ifndef DUMUX_DISCRETIZATION_CC_MPFA_L_INTERACTIONVOLUMESEED_HH
+#define DUMUX_DISCRETIZATION_CC_MPFA_L_INTERACTIONVOLUMESEED_HH
+
+#include <dumux/implicit/cellcentered/mpfa/properties.hh>
+#include "localsubcontrolentityseeds.hh"
+
+namespace Dumux
+{
+
+/*!
+ * \ingroup Mpfa
+ * \brief Base class for the interaction volume seed of mpfa methods
+ */
+template<typename G, typename L>
+class CCMpfaLInteractionVolumeSeed
+{
+ using GlobalIndexSet = G;
+ using GlobalIndexType = typename G::value_type;
+ using LocalIndexSet = L;
+ using LocalIndexType = typename L::value_type;
+
+public:
+ using LocalScvSeed = CCMpfaLCentralLocalScvSeed<G, L>;
+ using LocalOuterScvSeed = CCMpfaLOuterLocalScvSeed<G>;
+
+ CCMpfaLInteractionVolumeSeed(std::vector<LocalScvSeed>&& scvSeeds,
+ std::vector<LocalOuterScvSeed>&& outerScvSeeds,
+ std::vector<GlobalIndexType>&& globalScvfIndices)
+ : scvSeeds_(std::move(scvSeeds)),
+ outerScvSeeds_(std::move(outerScvSeeds)),
+ globalScvfIndices_(std::move(globalScvfIndices)) {}
+
+ const std::vector<LocalScvSeed>& scvSeeds() const
+ { return scvSeeds_; }
+
+ const LocalScvSeed& scvSeed(const LocalIndexType idx) const
+ { return scvSeeds_[idx]; }
+
+ const std::vector<LocalOuterScvSeed>& outerScvSeeds() const
+ { return outerScvSeeds_; }
+
+ const LocalOuterScvSeed& outerScvSeed(const LocalIndexType idx) const
+ { return outerScvSeeds_[idx]; }
+
+ std::vector<GlobalIndexType> globalScvIndices() const
+ {
+ std::vector<GlobalIndexType> globalIndices;
+ globalIndices.reserve(scvSeeds().size() + outerScvSeeds().size());
+
+ for (auto&& localScvSeed : scvSeeds())
+ globalIndices.push_back(localScvSeed.globalIndex());
+
+ for (auto&& localScvSeed : outerScvSeeds())
+ globalIndices.push_back(localScvSeed.globalIndex());
+
+ // make the entries unique
+ std::sort(globalIndices.begin(), globalIndices.end());
+ globalIndices.erase(std::unique(globalIndices.begin(), globalIndices.end()), globalIndices.end());
+
+ return globalIndices;
+ }
+
+ const std::vector<GlobalIndexType>& globalScvfIndices() const
+ { return globalScvfIndices_; }
+
+ bool isUnique() const
+ { return scvSeeds_.size() == 1; }
+
+private:
+ std::vector<LocalScvSeed> scvSeeds_;
+ std::vector<LocalOuterScvSeed> outerScvSeeds_;
+ std::vector<GlobalIndexType> globalScvfIndices_;
+};
+} // end namespace
+
+#endif
diff --git a/dumux/discretization/cellcentered/mpfa/lmethod/localsubcontrolentityseeds.hh b/dumux/discretization/cellcentered/mpfa/lmethod/localsubcontrolentityseeds.hh
new file mode 100644
index 0000000000000000000000000000000000000000..739ed82e5fc17e2514a17334203d230ea91ecd78
--- /dev/null
+++ b/dumux/discretization/cellcentered/mpfa/lmethod/localsubcontrolentityseeds.hh
@@ -0,0 +1,98 @@
+// -*- 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 sub control entity seeds of the mpfa-o method.
+ */
+#ifndef DUMUX_DISCRETIZATION_CC_MPFA_L_LOCALSUBCONTROLENTITYSEEDS_HH
+#define DUMUX_DISCRETIZATION_CC_MPFA_L_LOCALSUBCONTROLENTITYSEEDS_HH
+
+#include <dumux/implicit/cellcentered/mpfa/properties.hh>
+#include <dumux/discretization/cellcentered/mpfa/facetypes.hh>
+
+namespace Dumux
+{
+//! The scv seed class for scvs which could be the central scv in the actual interaction region
+template<typename G, typename L>
+class CCMpfaLCentralLocalScvSeed
+{
+ using GlobalIndexSet = G;
+ using LocalIndexSet = L;
+public:
+ using GlobalIndexType = typename GlobalIndexSet::value_type;
+ using LocalIndexType = typename LocalIndexSet::value_type;
+
+ //! constructor fully defining the scv seed
+ CCMpfaLCentralLocalScvSeed(GlobalIndexSet&& globalScvfIndices,
+ GlobalIndexType globalScvIndex,
+ LocalIndexType contiFaceLocalIdx)
+ : contiFaceLocalIdx_(contiFaceLocalIdx),
+ globalScvIndex_(globalScvIndex),
+ globalScvfIndices_(std::move(globalScvfIndices)) {}
+
+ LocalIndexType contiFaceLocalIdx() const
+ { return contiFaceLocalIdx_; }
+
+ GlobalIndexType globalIndex() const
+ { return globalScvIndex_; }
+
+ const GlobalIndexSet& globalScvfIndices() const
+ { return globalScvfIndices_; }
+
+private:
+ LocalIndexType contiFaceLocalIdx_;
+ GlobalIndexType globalScvIndex_;
+ GlobalIndexSet globalScvfIndices_;
+};
+
+//! The scv seed class for scvs which could be the central scv in the actual interaction region
+template<typename G>
+class CCMpfaLOuterLocalScvSeed
+{
+ using GlobalIndexSet = G;
+public:
+ using GlobalIndexType = typename G::value_type;
+
+ //! Constructor fully initializing the members
+ CCMpfaLOuterLocalScvSeed(GlobalIndexType globalScvIndex,
+ GlobalIndexType globalScvfIndex)
+ : scvIndexGlobal_(globalScvIndex),
+ scvfIndexGlobal_(globalScvfIndex) {}
+
+ //! Construct from central scv seed
+ template<typename GI, typename LI>
+ CCMpfaLOuterLocalScvSeed(const CCMpfaLCentralLocalScvSeed<GI, LI>& scvSeed)
+ : scvIndexGlobal_(scvSeed.globalIndex()),
+ scvfIndexGlobal_(scvSeed.globalScvfIndices()[scvSeed.contiFaceLocalIdx()]) {}
+
+
+ const GlobalIndexType globalIndex() const
+ { return scvIndexGlobal_; }
+
+ const GlobalIndexType globalScvfIndex() const
+ { return scvfIndexGlobal_; }
+
+private:
+ GlobalIndexType scvIndexGlobal_;
+ GlobalIndexType scvfIndexGlobal_;
+};
+
+} // end namespace
+
+#endif
diff --git a/dumux/discretization/cellcentered/mpfa/lmethod/subcontrolvolumeface.hh b/dumux/discretization/cellcentered/mpfa/lmethod/subcontrolvolumeface.hh
new file mode 100644
index 0000000000000000000000000000000000000000..13e90f8edc57ffaaea75353a59fe54a100b7250d
--- /dev/null
+++ b/dumux/discretization/cellcentered/mpfa/lmethod/subcontrolvolumeface.hh
@@ -0,0 +1,77 @@
+// -*- 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 Class for an mpfa-o sub control volume face
+ */
+#ifndef DUMUX_DISCRETIZATION_CC_MPFA_L_SUBCONTROLVOLUMEFACE_HH
+#define DUMUX_DISCRETIZATION_CC_MPFA_L_SUBCONTROLVOLUMEFACE_HH
+
+#include <dumux/discretization/cellcentered/mpfa/methods.hh>
+#include <dumux/discretization/cellcentered/mpfa/subcontrolvolumefacebase.hh>
+
+namespace Dumux
+{
+
+/*!
+ * \ingroup Discretization
+ * \brief Class for a sub control volume face in the mpfa-l method. We simply inherit from the base class here.
+ */
+template<class G, typename I>
+class CCMpfaSubControlVolumeFace<MpfaMethods::lMethod, G, I> : public CCMpfaSubControlVolumeFaceBase<G, I>
+{
+ using ParentType = CCMpfaSubControlVolumeFaceBase<G, I>;
+ using Geometry = G;
+ using IndexType = I;
+
+ using Scalar = typename Geometry::ctype;
+ static const int dim = Geometry::mydimension;
+ static const int dimworld = Geometry::coorddimension;
+
+ using GlobalPosition = Dune::FieldVector<Scalar, dimworld>;
+
+public:
+ //! We do not use the localIndex variable.
+ //! It is here to satisfy the general mpfa scvf interface.
+ template<class MpfaHelper>
+ CCMpfaSubControlVolumeFace(const MpfaHelper& helper,
+ std::vector<GlobalPosition>&& corners,
+ GlobalPosition&& unitOuterNormal,
+ IndexType vertexIndex,
+ unsigned int localIndex,
+ IndexType scvfIndex,
+ IndexType insideScvIdx,
+ const std::vector<IndexType>& outsideScvIndices,
+ Scalar q,
+ bool boundary)
+ : ParentType(helper,
+ std::move(corners),
+ std::move(unitOuterNormal),
+ vertexIndex,
+ localIndex,
+ scvfIndex,
+ insideScvIdx,
+ outsideScvIndices,
+ 0.0, // q should be always zero for the mpfa-l method,
+ boundary) {}
+};
+
+} // end namespace
+
+#endif
diff --git a/dumux/discretization/cellcentered/mpfa/methods.hh b/dumux/discretization/cellcentered/mpfa/methods.hh
index 9309a2568b99b9ec74e1e7d6b9bf4c7e85216318..07d47b78ac6b716f63ff92e0601ee172008ec03f 100644
--- a/dumux/discretization/cellcentered/mpfa/methods.hh
+++ b/dumux/discretization/cellcentered/mpfa/methods.hh
@@ -27,6 +27,7 @@ namespace Dumux
{
enum class MpfaMethods : unsigned int
{
+ lMethod,
oMethod,
oMethodFps
};
diff --git a/dumux/discretization/cellcentered/mpfa/mpfageometryhelper.hh b/dumux/discretization/cellcentered/mpfa/mpfageometryhelper.hh
deleted file mode 100644
index c98451ae93e5e595a02e684cfcb38252b503e824..0000000000000000000000000000000000000000
--- a/dumux/discretization/cellcentered/mpfa/mpfageometryhelper.hh
+++ /dev/null
@@ -1,236 +0,0 @@
-// -*- 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_GEOMETRYHELPER_HH
-#define DUMUX_DISCRETIZATION_CC_MPFA_GEOMETRYHELPER_HH
-
-#include <dune/geometry/multilineargeometry.hh>
-#include <dune/geometry/referenceelements.hh>
-
-#include <dumux/common/math.hh>
-
-namespace Dumux
-{
-//! A class to create sub control volume face geometries per intersection
-template <class GridView, int dim>
-class MpfaGeometryHelper
-{};
-
-//! Specialization for dim == 2
-template <class GridView>
-class MpfaGeometryHelper<GridView, 2>
-{
-private:
- using Scalar = typename GridView::ctype;
- static const int dim = GridView::dimension;
- static const int dimWorld = GridView::dimensionworld;
-
- using GlobalPosition = Dune::FieldVector<Scalar, dim>;
- using Element = typename GridView::template Codim<0>::Entity;
- using Intersection = typename GridView::Intersection;
-
- using FaceReferenceElements = typename Dune::ReferenceElements<Scalar, dim-1>;
-
-public:
- using PointVector = std::vector<GlobalPosition>;
-
- MpfaGeometryHelper(const typename Element::Geometry& elemGeom) : gt(elemGeom.type()) {}
-
- //! get sub control volume face corners of an intersection for the given local index
- static PointVector getScvfCorners(const typename Intersection::Geometry& geometry,
- unsigned int indexOnIntersection)
- {
- if (indexOnIntersection == 0)
- return PointVector({geometry.center(), geometry.corner(0)});
- else if (indexOnIntersection == 1)
- return PointVector({geometry.center(), geometry.corner(1)});
- else
- DUNE_THROW(Dune::InvalidStateException, "local index exceeds the number of corners of 2d intersections");
- }
-
- static GlobalPosition getScvfIntegrationPoint(const PointVector& scvfCorners, Scalar q)
- {
- auto d = scvfCorners[1];
- auto ip = scvfCorners[0];
- d -= ip;
- d *= q;
- ip += d;
-
- return ip;
- }
-
- static Scalar getScvfArea(const PointVector& scvfCorners)
- {
- return (scvfCorners[1]-scvfCorners[0]).two_norm();
- }
-
- std::size_t getNumLocalScvfs()
- {
- Dune::GeometryType triangle, quadrilateral;
- triangle.makeTriangle();
- quadrilateral.makeQuadrilateral();
-
- if (gt == triangle)
- return 6;
- else if (gt == quadrilateral)
- return 8;
- else
- DUNE_THROW(Dune::InvalidStateException, "unknown 2d geometry type " << gt);
- }
-
-private:
- const Dune::GeometryType gt;
-};
-
-//! Specialization for dim == 3
-template <class GridView>
-class MpfaGeometryHelper<GridView, 3>
-{
-private:
- using Scalar = typename GridView::ctype;
- static const int dim = GridView::dimension;
- static const int dimWorld = GridView::dimensionworld;
-
- using GlobalPosition = Dune::FieldVector<Scalar, dim>;
- using Element = typename GridView::template Codim<0>::Entity;
- using Intersection = typename GridView::Intersection;
-
- using FaceReferenceElements = typename Dune::ReferenceElements<Scalar, dim-1>;
-
-public:
- using PointVector = std::vector<GlobalPosition>;
-
- MpfaGeometryHelper(const typename Element::Geometry& elemGeom) : gt(elemGeom.type()) {}
-
- //! get sub control volume face corners of an intersection for the given local index
- static PointVector getScvfCorners(const typename Intersection::Geometry& geometry,
- unsigned int indexOnIntersection)
- {
- // extract the corners of the sub control volumes
- const auto& referenceElement = FaceReferenceElements::general(geometry.type());
-
- // maximum number of necessary points is 9 (for quadrilateral)
- GlobalPosition p[9];
- auto corners = geometry.corners();
-
- // the intersection center
- p[0] = geometry.center();
-
- // vertices
- for (int i = 0; i < corners; ++i)
- p[i+1] = geometry.corner(i);
-
- // edge midpoints
- for (int i = 0; i < referenceElement.size(1); ++i)
- p[i+corners+1] = geometry.global(referenceElement.position(i, 1));
-
- // proceed according to number of corners
- switch (corners)
- {
- case 3: // triangle
- {
- //! Only build the maps the first time we encounter a triangle
- static const std::uint8_t vo = 1; //! vertex offset in point vector p
- static const std::uint8_t eo = 4; //! edge offset in point vector p
- static const std::uint8_t map[3][4] =
- {
- {0, eo+1, eo+0, vo+0},
- {0, eo+0, eo+2, vo+1},
- {0, eo+2, eo+1, vo+2}
- };
-
- return PointVector( {p[map[indexOnIntersection][0]],
- p[map[indexOnIntersection][1]],
- p[map[indexOnIntersection][2]],
- p[map[indexOnIntersection][3]]} );
- }
- case 4: // quadrilateral
- {
- //! Only build the maps the first time we encounter a quadrilateral
- static const std::uint8_t vo = 1; //! vertex offset in point vector p
- static const std::uint8_t eo = 5; //! face offset in point vector p
- static const std::uint8_t map[4][4] =
- {
- {0, eo+0, eo+2, vo+0},
- {0, eo+2, eo+1, vo+1},
- {0, eo+3, eo+0, vo+2},
- {0, eo+1, eo+3, vo+3}
- };
-
- return PointVector( {p[map[indexOnIntersection][0]],
- p[map[indexOnIntersection][1]],
- p[map[indexOnIntersection][2]],
- p[map[indexOnIntersection][3]]} );
- }
- default:
- DUNE_THROW(Dune::NotImplemented, "Mpfa scvf corners for dim=" << dim
- << " dimWorld=" << dimWorld
- << " corners=" << corners);
- }
- }
-
- static GlobalPosition getScvfIntegrationPoint(const PointVector& scvfCorners, Scalar q)
- {
- auto d = scvfCorners[3];
- auto ip = scvfCorners[0];
- d -= ip;
- d *= q;
- ip += d;
-
- return ip;
- }
-
- static Scalar getScvfArea(const PointVector& scvfCorners)
- {
- // after Wolfram alpha quadrilateral area
- return 0.5*Dumux::crossProduct(scvfCorners[3]-scvfCorners[0], scvfCorners[2]-scvfCorners[1]).two_norm();
- }
-
- std::size_t getNumLocalScvfs()
- {
- Dune::GeometryType tetrahedron, pyramid, prism, hexahedron;
- tetrahedron.makeTetrahedron();
- pyramid.makePyramid();
- prism.makePrism();
- hexahedron.makeHexahedron();
-
- if (gt == tetrahedron)
- return 12;
- else if (gt == pyramid)
- return 16;
- else if (gt == prism)
- return 18;
- else if (gt == hexahedron)
- return 24;
- else
- DUNE_THROW(Dune::InvalidStateException, "unknown 3d geometry type " << gt);
- }
-
-private:
- const Dune::GeometryType gt; // the geometry type of the element
-};
-
-} // end namespace
-
-#endif
\ No newline at end of file
diff --git a/dumux/discretization/cellcentered/mpfa/omethod/globalfvgeometry.hh b/dumux/discretization/cellcentered/mpfa/omethod/globalfvgeometry.hh
deleted file mode 100644
index 1e8288198553c7f29e4562c6a633f710bb90543e..0000000000000000000000000000000000000000
--- a/dumux/discretization/cellcentered/mpfa/omethod/globalfvgeometry.hh
+++ /dev/null
@@ -1,46 +0,0 @@
-// -*- 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_GLOBALFVGEOMETRY_HH
-#define DUMUX_DISCRETIZATION_CC_MPFA_O_GLOBALFVGEOMETRY_HH
-
-#include <dumux/discretization/cellcentered/mpfa/globalfvgeometrybase.hh>
-#include <dumux/discretization/cellcentered/mpfa/methods.hh>
-
-namespace Dumux
-{
-//! Specialization for the mpfa-o scheme.
-template<class TypeTag, bool EnableCache>
-class CCMpfaGlobalFVGeometryImplementation<TypeTag, MpfaMethods::oMethod, EnableCache> : public CCMpfaGlobalFVGeometryBase<TypeTag, EnableCache>
-{
- using ParentType = CCMpfaGlobalFVGeometryBase<TypeTag, EnableCache>;
- using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
-
-public:
- CCMpfaGlobalFVGeometryImplementation(const GridView gridView) : ParentType(gridView) {}
-};
-
-} // end namespace
-
-#endif
diff --git a/dumux/discretization/cellcentered/mpfa/omethod/helper.hh b/dumux/discretization/cellcentered/mpfa/omethod/helper.hh
index 2de4395fcaaaeec48dcc61b8ee90f1855315a0aa..d12415c2da7aad1940094008b73b67bbc86b0980 100644
--- a/dumux/discretization/cellcentered/mpfa/omethod/helper.hh
+++ b/dumux/discretization/cellcentered/mpfa/omethod/helper.hh
@@ -35,38 +35,38 @@ 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.
+ * Specialization for the Mpfa-O method in two dimensions embedded in a 2d world.
*/
template<class TypeTag>
-class MpfaHelperBase<TypeTag, MpfaMethods::oMethod, 2>
+class MpfaMethodHelper<TypeTag, MpfaMethods::oMethod, /*dim*/2, /*dimWorld*/2>
{
using Implementation = typename GET_PROP_TYPE(TypeTag, MpfaHelper);
+
+ static const int dim = 2;
+ static const int dimWorld = 2;
+
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 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;
- using Element = typename GridView::template Codim<0>::Entity;
- using GlobalIndexType = typename GridView::IndexSet::IndexType;
- using LocalIndexType = typename InteractionVolume::LocalIndexType;
+ // The o method is used on the boundaries also for other methods. These will use this helper
+ // class to set up the seeds on the boundary. Therefore we extract all seed information from
+ // the boundary interaction volume here to be compatible with other mpfa methods where seed types
+ // most probably differ. We use the fact here, that for the o-method, boundary and interior
+ // interaction volumes are identical, which is always given.
+ using InteractionVolume = typename GET_PROP_TYPE(TypeTag, BoundaryInteractionVolume);
using InteractionVolumeSeed = typename InteractionVolume::Seed;
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>;
+ using Element = typename GridView::template Codim<0>::Entity;
+ using GlobalIndexType = typename InteractionVolume::GlobalIndexType;
+ using LocalIndexType = typename InteractionVolume::LocalIndexType;
- using DimVector = Dune::FieldVector<Scalar, dimWorld>;
+ using GlobalIndexSet = typename InteractionVolume::GlobalIndexSet;
+ using LocalIndexSet = typename InteractionVolume::LocalIndexSet;
public:
static InteractionVolumeSeed makeInnerInteractionVolumeSeed(const Problem& problem,
const Element& element,
@@ -89,13 +89,13 @@ public:
return InteractionVolumeSeed(std::move(scvSeeds), std::move(scvfSeeds), false);
}
- static BoundaryInteractionVolumeSeed makeBoundaryInteractionVolumeSeed(const Problem& problem,
- const Element& element,
- const FVElementGeometry& fvGeometry,
- const SubControlVolumeFace& scvf)
+ static InteractionVolumeSeed makeBoundaryInteractionVolumeSeed(const Problem& problem,
+ const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const SubControlVolumeFace& scvf)
{
- std::vector<BoundaryScvSeed> scvSeeds;
- std::vector<BoundaryScvfSeed> scvfSeeds;
+ std::vector<ScvSeed> scvSeeds;
+ std::vector<ScvfSeed> scvfSeeds;
// reserve sufficient memory
scvSeeds.reserve(20);
@@ -107,13 +107,13 @@ public:
scvSeeds.shrink_to_fit();
scvfSeeds.shrink_to_fit();
- return BoundaryInteractionVolumeSeed(std::move(scvSeeds), std::move(scvfSeeds), true);
+ return InteractionVolumeSeed(std::move(scvSeeds), std::move(scvfSeeds), true);
}
private:
- template<typename ScvSeed, typename ScvfSeed>
- static void fillEntitySeeds_(std::vector<ScvSeed>& scvSeeds,
- std::vector<ScvfSeed>& scvfSeeds,
+ template<typename ScvSeedType, typename ScvfSeedType>
+ static void fillEntitySeeds_(std::vector<ScvSeedType>& scvSeeds,
+ std::vector<ScvfSeedType>& scvfSeeds,
const Problem& problem,
const Element& element,
const FVElementGeometry& fvGeometry,
@@ -140,23 +140,25 @@ private:
performRotation_(problem, scvfVector, scvSeeds, scvfSeeds, scvIdx0, /*clockwise*/true);
// Finish by creating the first scv
- scvSeeds.emplace(scvSeeds.begin(), ScvSeed(GlobalIndexSet({scvfVector[0]->index(), scvfVector[1]->index()}),
- LocalIndexSet({0, storeIdx}),
- scvIdx0));
+ scvSeeds.emplace(scvSeeds.begin(), 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));
+ scvSeeds.emplace(scvSeeds.begin(), GlobalIndexSet({scvfVector[0]->index(), scvfVector[1]->index()}),
+ LocalIndexSet({0, static_cast<LocalIndexType>(scvfSeeds.size()-1)}),
+ scvIdx0);
}
- // clockwise rotation and construction of local scv & scv face entities
- template<typename ScvSeed, typename ScvfSeed, typename ScvfPointerVector>
+ // in 2d we can make use of knowing the basis orientation (right hand system)
+ // that way we can find the interaction regions more rapidly by clockwise and
+ // counter clockwise rotation and construction of local scv & scv face entities
+ template<typename ScvSeedType, typename ScvfSeedType, typename ScvfPointerVector>
static void performRotation_(const Problem& problem,
const ScvfPointerVector& scvfVector,
- std::vector<ScvSeed>& scvSeeds,
- std::vector<ScvfSeed>& scvfSeeds,
+ std::vector<ScvSeedType>& scvSeeds,
+ std::vector<ScvfSeedType>& scvfSeeds,
const GlobalIndexType scvIdx0,
const bool clockWise = false)
{
@@ -176,23 +178,22 @@ private:
while (!finished)
{
// Get some indices beforehand
- GlobalIndexType globalScvfIdx = curScvf.index();
- GlobalIndexType insideGlobalScvIdx = curScvf.insideScvIdx();
GlobalIndexType outsideGlobalScvIdx = curScvf.outsideScvIdx();
LocalIndexType insideLocalScvIdx = firstIteration ? 0 : localScvIdx;
// the current element inside of the scv face
- auto insideElement = problem.model().globalFvGeometry().element(insideGlobalScvIdx);
+ auto insideElement = problem.model().globalFvGeometry().element(curScvf.insideScvIdx());
auto faceType = Implementation::getMpfaFaceType(problem, insideElement, curScvf);
// if the face touches the boundary, create a boundary scvf entity
if (curScvf.boundary())
{
assert(faceType == MpfaFaceTypes::neumann || faceType == MpfaFaceTypes::dirichlet);
- scvfSeeds.emplace_back(ScvfSeed(curScvf,
- LocalIndexSet({insideLocalScvIdx}),
- GlobalIndexSet({globalScvfIdx}),
- faceType));
+ scvfSeeds.emplace_back( curScvf,
+ insideLocalScvIdx,
+ LocalIndexSet(),
+ GlobalIndexSet(),
+ faceType );
// rotation loop is finished
finished = true; return;
}
@@ -201,10 +202,11 @@ private:
if (outsideGlobalScvIdx == scvIdx0)
{
// create scv face entity for the last face of the loop
- scvfSeeds.emplace_back(ScvfSeed(curScvf,
- LocalIndexSet({insideLocalScvIdx, 0}),
- GlobalIndexSet({globalScvfIdx, scvfVector[1]->index()}),
- faceType));
+ scvfSeeds.emplace_back( curScvf,
+ insideLocalScvIdx,
+ LocalIndexSet({0}),
+ GlobalIndexSet({scvfVector[1]->index()}),
+ faceType );
// rotation loop is finished
finished = true; return;
@@ -222,13 +224,12 @@ private:
auto& nextScvf = *outsideScvfVector[nextFaceCoordIdx];
// create local scv face entity of the current scvf
- GlobalIndexType commonGlobalScvfIdx = commonScvf.index();
LocalIndexType outsideLocalScvIdx = localScvIdx+1;
-
- scvfSeeds.emplace_back(ScvfSeed(curScvf,
- LocalIndexSet({insideLocalScvIdx, outsideLocalScvIdx}),
- GlobalIndexSet({globalScvfIdx, commonGlobalScvfIdx}),
- faceType));
+ scvfSeeds.emplace_back( curScvf,
+ insideLocalScvIdx,
+ LocalIndexSet({outsideLocalScvIdx}),
+ GlobalIndexSet({commonScvf.index()}),
+ faceType );
localScvfIdx++;
// create index set storing the two local scvf indices
@@ -238,7 +239,7 @@ private:
// create "outside" scv
GlobalIndexSet globalScvfIndices({outsideScvfVector[0]->index(), outsideScvfVector[1]->index()});
- scvSeeds.emplace_back(ScvSeed(std::move(globalScvfIndices), std::move(localScvfs), curScvf.outsideScvIdx()));
+ scvSeeds.emplace_back(std::move(globalScvfIndices), std::move(localScvfs), outsideGlobalScvIdx);
localScvIdx++;
// create the next scvf in the following iteration
@@ -251,38 +252,292 @@ private:
/*!
* \ingroup Mpfa
* \brief Helper class to get the required information on an interaction volume.
- * Specialization for the Mpfa-O method in three dimensions.
+ * Specialization for the Mpfa-O method in two dimensions embedded in a 3d world.
*/
template<class TypeTag>
-class MpfaHelperBase<TypeTag, MpfaMethods::oMethod, 3>
+class MpfaMethodHelper<TypeTag, MpfaMethods::oMethod, /*dim*/2, /*dimWorld*/3>
{
using Implementation = typename GET_PROP_TYPE(TypeTag, MpfaHelper);
+
+ static const int dim = 2;
+ static const int dimWorld = 3;
+
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 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;
+ // The o method is used on the boundaries also for other methods. These will use this helper
+ // class to set up the seeds on the boundary. Therefore we extract all seed information from
+ // the boundary interaction volume here to be compatible with other mpfa methods where seed types
+ // most probably differ. We use the fact here, that for the o-method, boundary and interior
+ // interaction volumes are identical, which is always given.
+ using InteractionVolume = typename GET_PROP_TYPE(TypeTag, BoundaryInteractionVolume);
+ using InteractionVolumeSeed = typename InteractionVolume::Seed;
+ using ScvSeed = typename InteractionVolumeSeed::LocalScvSeed;
+ using ScvfSeed = typename InteractionVolumeSeed::LocalScvfSeed;
+
using Element = typename GridView::template Codim<0>::Entity;
- using GlobalIndexType = typename GridView::IndexSet::IndexType;
+ using GlobalIndexType = typename InteractionVolume::GlobalIndexType;
using LocalIndexType = typename InteractionVolume::LocalIndexType;
+ using GlobalIndexSet = typename InteractionVolume::GlobalIndexSet;
+ using LocalIndexSet = typename InteractionVolume::LocalIndexSet;
+public:
+ static InteractionVolumeSeed makeInnerInteractionVolumeSeed(const Problem& problem,
+ const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const SubControlVolumeFace& scvf)
+ {
+ std::vector<ScvSeed> scvSeeds;
+ std::vector<ScvfSeed> scvfSeeds;
+
+ // reserve sufficient memory
+ scvSeeds.reserve(20);
+ scvfSeeds.reserve(20);
+
+ // The vertex index around which we construct the interaction volume
+ auto vIdxGlobal = scvf.vertexIndex();
+
+ // Get the two scv faces in the scv
+ auto scvfVector = Implementation::getScvFacesAtVertex(vIdxGlobal, element, fvGeometry);
+
+ // fill the entity seed data
+ fillEntitySeeds_(scvSeeds, scvfSeeds, problem, element, fvGeometry, scvfVector, vIdxGlobal);
+
+ // shrink containers to necessary size
+ scvSeeds.shrink_to_fit();
+ scvfSeeds.shrink_to_fit();
+
+ for (auto& scvfSeed : scvfSeeds)
+ scvfSeed.makeOutsideDataUnique();
+
+ return InteractionVolumeSeed(std::move(scvSeeds), std::move(scvfSeeds), false);
+ }
+
+ static InteractionVolumeSeed makeBoundaryInteractionVolumeSeed(const Problem& problem,
+ const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const SubControlVolumeFace& scvf)
+ {
+ std::vector<ScvSeed> scvSeeds;
+ std::vector<ScvfSeed> scvfSeeds;
+
+ // reserve sufficient memory
+ scvSeeds.reserve(20);
+ scvfSeeds.reserve(20);
+
+ // The vertex index around which we construct the interaction volume
+ auto vIdxGlobal = scvf.vertexIndex();
+
+ // Get the two scv faces in the scv
+ auto scvfVector = Implementation::getScvFacesAtVertex(vIdxGlobal, element, fvGeometry);
+
+ // fill the entity seed data
+ fillEntitySeeds_(scvSeeds, scvfSeeds, problem, element, fvGeometry, scvfVector, vIdxGlobal);
+
+ // shrink containers to necessary size
+ scvSeeds.shrink_to_fit();
+ scvfSeeds.shrink_to_fit();
+
+ for (auto& scvfSeed : scvfSeeds)
+ scvfSeed.makeOutsideDataUnique();
+
+ return InteractionVolumeSeed(std::move(scvSeeds), std::move(scvfSeeds), true);
+ }
+
+private:
+ template<class ScvSeedType, class ScvfSeedType, class ScvfVector>
+ static void fillEntitySeeds_(std::vector<ScvSeedType>& scvSeeds,
+ std::vector<ScvfSeedType>& scvfSeeds,
+ const Problem& problem,
+ const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const ScvfVector& scvfVector,
+ const GlobalIndexType vIdxGlobal)
+ {
+ // make the scv without knowing the local scvf indices yet
+ // take the inside scv index from the first scvf (is the same for all)
+ scvSeeds.emplace_back( GlobalIndexSet({scvfVector[0]->index(),
+ scvfVector[1]->index()}),
+ scvfVector[0]->insideScvIdx());
+
+ // make the scvf seeds for the two scvfs connected to the scv
+ auto& actualScvSeed = scvSeeds.back();
+ LocalIndexType actualLocalScvIdx = scvSeeds.size()-1;
+
+ for (int coordDir = 0; coordDir < dim; ++coordDir)
+ {
+ const auto& actualScvf = *scvfVector[coordDir];
+
+ // if scvf is on a boundary, we create the scvfSeed and make no neighbor
+ if (actualScvf.boundary())
+ {
+ // set the local scvfIndex of the face that is about to created
+ actualScvSeed.setLocalScvfIndex(coordDir, scvfSeeds.size());
+
+ // create the scvf seed
+ scvfSeeds.emplace_back( actualScvf,
+ actualLocalScvIdx,
+ LocalIndexSet(),
+ GlobalIndexSet(),
+ Implementation::getMpfaFaceType(problem, element, actualScvf) );
+ }
+ else
+ {
+ // we loop over all neighbors of this face
+ for (auto outsideGlobalScvIdx : actualScvf.outsideScvIndices())
+ {
+ // get outside element, fvgeometry etc.
+ auto outsideElement = problem.model().globalFvGeometry().element(outsideGlobalScvIdx);
+ auto outsideFvGeometry = localView(problem.model().globalFvGeometry());
+ outsideFvGeometry.bindElement(outsideElement);
+
+ // find scvf in outside corresponding to the actual scvf
+ auto outsideScvfVector = Implementation::getScvFacesAtVertex(vIdxGlobal, outsideElement, outsideFvGeometry);
+ auto commonFaceLocalIdx = Implementation::getCommonFaceLocalIndex(actualScvf, outsideScvfVector);
+ const auto& outsideScvf = *outsideScvfVector[commonFaceLocalIdx];
+ const auto outsideScvfIdx = outsideScvf.index();
+
+ // check if the outside scv already exists and get its local index
+ bool outsideScvExists = false;
+ LocalIndexType outsideLocalScvIdx = 0;
+ for (auto&& scvSeed : scvSeeds)
+ {
+ if (scvSeed.globalIndex() == outsideGlobalScvIdx)
+ {
+ outsideScvExists = true;
+ break;
+ }
+ // keep track of local index
+ outsideLocalScvIdx++;
+ }
+
+ // check if face already exists inside or in any neighbor
+ bool outsideScvfExists = false;
+ bool insideScvfExists = false;
+ LocalIndexType outsideLocalScvfIdx = 0;
+ LocalIndexType insideLocalScvfIdx = 0;
+ for (auto&& scvfSeed : scvfSeeds)
+ {
+ if (Implementation::contains(actualScvf.outsideScvIndices(), scvfSeed.insideGlobalScvIndex()) &&
+ Implementation::contains(scvfSeed.outsideGlobalScvIndices(), actualScvf.insideScvIdx()))
+ outsideScvfExists = true;
+ else
+ if (!outsideScvfExists)
+ outsideLocalScvfIdx++;
+
+ if (scvfSeed.insideGlobalScvIndex() == actualScvf.insideScvIdx() &&
+ std::equal(scvfSeed.outsideGlobalScvIndices().begin(), scvfSeed.outsideGlobalScvIndices().end(), actualScvf.outsideScvIndices().begin()))
+ insideScvfExists = true;
+ else
+ if (!insideScvfExists)
+ insideLocalScvfIdx++;
+ }
+
+ // we should never have two local scv faces existing inside and in an outside element
+ assert(!(outsideScvfExists && insideScvfExists) && "The scv face seems to exist twice!");
+
+ // outside scv has to be created
+ if (!outsideScvExists)
+ {
+ // The face does not yet exist
+ if (!insideScvfExists && !outsideScvfExists)
+ {
+ // set the local scvfIndex of the face that is about to created
+ actualScvSeed.setLocalScvfIndex(coordDir, scvfSeeds.size());
+
+ // create scvf seed
+ scvfSeeds.emplace_back(actualScvf,
+ actualLocalScvIdx,
+ Implementation::getMpfaFaceType(problem, element, actualScvf));
+
+ // pass the actual outside indices to the new scvf seed
+ scvfSeeds.back().addOutsideScvfIndex(outsideScvfIdx);
+ scvfSeeds.back().addOutsideLocalScvIndex(static_cast<LocalIndexType>(scvSeeds.size()));
+ }
+ else if (insideScvfExists && !outsideScvfExists)
+ {
+ // pass info on outside to the inside scvf seed
+ scvfSeeds[insideLocalScvfIdx].addOutsideScvfIndex(outsideScvfIdx);
+ scvfSeeds[insideLocalScvfIdx].addOutsideLocalScvIndex(static_cast<LocalIndexType>(scvSeeds.size()));
+ }
+ else if (!insideScvfExists && outsideScvfExists)
+ {
+ // pass info on outside to the inside scvf seed
+ scvfSeeds[outsideLocalScvfIdx].addOutsideScvfIndex(actualScvf.index());
+ scvfSeeds[outsideLocalScvfIdx].addOutsideLocalScvIndex(static_cast<LocalIndexType>(scvSeeds.size()));
+ }
+
+ // make outside scv by recursion
+ fillEntitySeeds_(scvSeeds, scvfSeeds, problem, outsideElement, outsideFvGeometry, outsideScvfVector, vIdxGlobal);
+ }
+ // outside scv exists, but no corresponding local scv face yet. Make it from inside.
+ else if (outsideScvExists && !outsideScvfExists && !insideScvfExists)
+ {
+ // set the local scvfIndex of the face that is about to created
+ actualScvSeed.setLocalScvfIndex(coordDir, scvfSeeds.size());
+
+ // create scvf seed
+ scvfSeeds.emplace_back(actualScvf,
+ actualLocalScvIdx,
+ Implementation::getMpfaFaceType(problem, element, actualScvf));
+
+ // pass the actual outside indices to the new scvf seed
+ scvfSeeds.back().addOutsideScvfIndex(outsideScvfIdx);
+ scvfSeeds.back().addOutsideLocalScvIndex(outsideLocalScvIdx);
+ }
+ else if (outsideScvExists && !insideScvfExists && outsideScvfExists)
+ {
+ // set the local scvfIndex of the found outside local scv face seed
+ actualScvSeed.setLocalScvfIndex(coordDir, outsideLocalScvfIdx);
+
+ // pass info on outside to the inside found local scvf seed
+ scvfSeeds[outsideLocalScvfIdx].addOutsideScvfIndex(outsideScvfIdx);
+ scvfSeeds[outsideLocalScvfIdx].addOutsideLocalScvIndex(actualLocalScvIdx);
+ }
+ }
+ }
+ }
+ }
+};
+
+/*!
+ * \ingroup Mpfa
+ * \brief Helper class to get the required information on an interaction volume.
+ * Specialization for the Mpfa-O method in three dimensions.
+ */
+template<class TypeTag>
+class MpfaMethodHelper<TypeTag, MpfaMethods::oMethod, /*dim*/3, /*dimWorld*/3>
+{
+ using Implementation = typename GET_PROP_TYPE(TypeTag, MpfaHelper);
+
+ static const int dim = 3;
+ static const int dimWorld = 3;
+
+ 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 FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
+ using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
+
+ // The o method is used on the boundaries also for other methods. These will use this helper
+ // class to set up the seeds on the boundary. Therefore we extract all seed information from
+ // the boundary interaction volume here to be compatible with other mpfa methods where seed types
+ // most probably differ. We use the fact here, that for the o-method, boundary and interior
+ // interaction volumes are identical, which is always given.
+ using InteractionVolume = typename GET_PROP_TYPE(TypeTag, BoundaryInteractionVolume);
using InteractionVolumeSeed = typename InteractionVolume::Seed;
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>;
+ using Element = typename GridView::template Codim<0>::Entity;
+ using GlobalIndexType = typename InteractionVolume::GlobalIndexType;
+ using LocalIndexType = typename InteractionVolume::LocalIndexType;
- using DimVector = Dune::FieldVector<Scalar, dimWorld>;
+ using GlobalIndexSet = typename InteractionVolume::GlobalIndexSet;
+ using LocalIndexSet = typename InteractionVolume::LocalIndexSet;
public:
static InteractionVolumeSeed makeInnerInteractionVolumeSeed(const Problem& problem,
@@ -300,8 +555,11 @@ public:
// The vertex index around which we construct the interaction volume
auto vIdxGlobal = scvf.vertexIndex();
+ // Get the three scv faces in the scv
+ auto scvfVector = Implementation::getScvFacesAtVertex(vIdxGlobal, element, fvGeometry);
+
// create the scv entity seeds
- fillEntitySeeds_(scvSeeds, scvfSeeds, problem, element, fvGeometry, scvf, vIdxGlobal);
+ fillEntitySeeds_(scvSeeds, scvfSeeds, problem, element, fvGeometry, scvfVector, vIdxGlobal);
// shrink containers to necessary size
scvSeeds.shrink_to_fit();
@@ -310,13 +568,13 @@ public:
return InteractionVolumeSeed(std::move(scvSeeds), std::move(scvfSeeds), false);
}
- static BoundaryInteractionVolumeSeed makeBoundaryInteractionVolumeSeed(const Problem& problem,
- const Element& element,
- const FVElementGeometry& fvGeometry,
- const SubControlVolumeFace& scvf)
+ static InteractionVolumeSeed makeBoundaryInteractionVolumeSeed(const Problem& problem,
+ const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const SubControlVolumeFace& scvf)
{
- std::vector<BoundaryScvSeed> scvSeeds;
- std::vector<BoundaryScvfSeed> scvfSeeds;
+ std::vector<ScvSeed> scvSeeds;
+ std::vector<ScvfSeed> scvfSeeds;
// reserve sufficient memory
scvSeeds.reserve(50);
@@ -325,37 +583,35 @@ public:
// The vertex index around which we construct the interaction volume
auto vIdxGlobal = scvf.vertexIndex();
+ // Get the three scv faces in the scv
+ auto scvfVector = Implementation::getScvFacesAtVertex(vIdxGlobal, element, fvGeometry);
+
// create the scv entity seeds
- fillEntitySeeds_(scvSeeds, scvfSeeds, problem, element, fvGeometry, scvf, vIdxGlobal);
+ fillEntitySeeds_(scvSeeds, scvfSeeds, problem, element, fvGeometry, scvfVector, vIdxGlobal);
// shrink containers to necessary size
scvSeeds.shrink_to_fit();
scvfSeeds.shrink_to_fit();
- return BoundaryInteractionVolumeSeed(std::move(scvSeeds), std::move(scvfSeeds), true);
+ return InteractionVolumeSeed(std::move(scvSeeds), std::move(scvfSeeds), true);
}
private:
- template<typename ScvSeed, typename ScvfSeed>
- static void fillEntitySeeds_(std::vector<ScvSeed>& scvSeeds,
- std::vector<ScvfSeed>& scvfSeeds,
+ template<class ScvSeedType, class ScvfSeedType, class ScvfVector>
+ static void fillEntitySeeds_(std::vector<ScvSeedType>& scvSeeds,
+ std::vector<ScvfSeedType>& scvfSeeds,
const Problem& problem,
const Element& element,
const FVElementGeometry& fvGeometry,
- const SubControlVolumeFace& scvf,
+ const ScvfVector& scvfVector,
const GlobalIndexType vIdxGlobal)
{
- // Get the three scv faces in the scv
- auto scvfVector = Implementation::getScvFacesAtVertex(vIdxGlobal, element, fvGeometry);
-
- // global scvIdx and global scvf indices
- auto globalScvIdx = scvf.insideScvIdx();
- GlobalIndexSet globalScvfIndices( {scvfVector[0]->index(),
- scvfVector[1]->index(),
- scvfVector[2]->index()} );
-
- // make the scv
- scvSeeds.emplace_back(ScvSeed(std::move(globalScvfIndices), globalScvIdx));
+ // make the scv without knowing the local scvf indices yet
+ // take the inside scv index from the first scvf (is the same for all)
+ scvSeeds.emplace_back( GlobalIndexSet({scvfVector[0]->index(),
+ scvfVector[1]->index(),
+ scvfVector[2]->index()}),
+ scvfVector[0]->insideScvIdx());
// make the scvf seeds for the three scvfs connected to the scv
auto& actualScvSeed = scvSeeds.back();
@@ -372,10 +628,11 @@ private:
actualScvSeed.setLocalScvfIndex(coordDir, scvfSeeds.size());
// create the scvf seed
- auto faceType = Implementation::getMpfaFaceType(problem, element, *scvfVector[coordDir]);
+ auto faceType = Implementation::getMpfaFaceType(problem, element, actualScvf);
scvfSeeds.emplace_back( actualScvf,
- LocalIndexSet({actualLocalScvIdx}),
- GlobalIndexSet({scvfVector[coordDir]->index()}),
+ actualLocalScvIdx,
+ LocalIndexSet(),
+ GlobalIndexSet(),
faceType );
}
else
@@ -410,17 +667,18 @@ private:
// find scvf in outside corresponding to the actual scvf
auto outsideScvfVector = Implementation::getScvFacesAtVertex(vIdxGlobal, outsideElement, outsideFvGeometry);
- auto commonFaceLocalIdx = Implementation::getCommonFaceLocalIndex(*scvfVector[coordDir], outsideScvfVector);
+ auto commonFaceLocalIdx = Implementation::getCommonFaceLocalIndex(actualScvf, outsideScvfVector);
const auto& outsideScvf = *outsideScvfVector[commonFaceLocalIdx];
// create scvf seed
- auto faceType = Implementation::getMpfaFaceType(problem, element, actualScvf);
- LocalIndexSet scvIndicesLocal({actualLocalScvIdx, static_cast<LocalIndexType>(scvSeeds.size())});
- GlobalIndexSet scvfIndicesGlobal({globalScvfIndex, outsideScvf.index()});
- scvfSeeds.emplace_back(actualScvf, std::move(scvIndicesLocal), std::move(scvfIndicesGlobal), faceType);
+ scvfSeeds.emplace_back(actualScvf,
+ actualLocalScvIdx,
+ LocalIndexSet( {static_cast<LocalIndexType>(scvSeeds.size())} ),
+ GlobalIndexSet( {outsideScvf.index()} ),
+ Implementation::getMpfaFaceType(problem, element, actualScvf));
// make outside scv by recursion
- fillEntitySeeds_(scvSeeds, scvfSeeds, problem, outsideElement, outsideFvGeometry, outsideScvf, vIdxGlobal);
+ fillEntitySeeds_(scvSeeds, scvfSeeds, problem, outsideElement, outsideFvGeometry, outsideScvfVector, vIdxGlobal);
}
// we have to find out if it is necessary to create a new scvf
else
@@ -456,14 +714,15 @@ private:
// find scvf in outside corresponding to the actual scvf
auto outsideScvfVector = Implementation::getScvFacesAtVertex(vIdxGlobal, outsideElement, outsideFvGeometry);
- auto commonFaceLocalIdx = Implementation::getCommonFaceLocalIndex(*scvfVector[coordDir], outsideScvfVector);
+ auto commonFaceLocalIdx = Implementation::getCommonFaceLocalIndex(actualScvf, outsideScvfVector);
const auto& outsideScvf = *outsideScvfVector[commonFaceLocalIdx];
- // some data on the face
- auto faceType = Implementation::getMpfaFaceType(problem, element, *scvfVector[coordDir]);
- LocalIndexSet scvIndicesLocal( {actualLocalScvIdx, outsideLocalScvIdx} );
- GlobalIndexSet scvfIndicesGlobal( {globalScvfIndex, outsideScvf.index()} );
- scvfSeeds.emplace_back(*scvfVector[coordDir], std::move(scvIndicesLocal), std::move(scvfIndicesGlobal), faceType);
+ // make scv face seed
+ scvfSeeds.emplace_back(actualScvf,
+ actualLocalScvIdx,
+ LocalIndexSet( {outsideLocalScvIdx} ),
+ GlobalIndexSet( {outsideScvf.index()} ),
+ Implementation::getMpfaFaceType(problem, element, actualScvf));
}
}
}
diff --git a/dumux/discretization/cellcentered/mpfa/omethod/interactionvolume.hh b/dumux/discretization/cellcentered/mpfa/omethod/interactionvolume.hh
index dc55c25bf95eb127806203f7fbe181b977980c76..245ba09e53332a44c0935394ce600eb49a43c020 100644
--- a/dumux/discretization/cellcentered/mpfa/omethod/interactionvolume.hh
+++ b/dumux/discretization/cellcentered/mpfa/omethod/interactionvolume.hh
@@ -36,35 +36,28 @@ namespace Dumux
{
//! Specialization of the interaction volume traits class for the mpfa-o method
template<class TypeTag>
-class CCMpfaOInteractionVolumeTraits
+class CCMpfaOInteractionVolumeTraits : public CCMpfaInteractionVolumeTraitsBase<TypeTag>
{
+ using BaseTraits = CCMpfaInteractionVolumeTraitsBase<TypeTag>;
using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
static const int dim = GridView::dimension;
static const int dimWorld = GridView::dimensionworld;
-public:
- using BoundaryInteractionVolume = CCMpfaInteractionVolumeImplementation<TypeTag, MpfaMethods::oMethod>;
-
using GlobalPosition = Dune::FieldVector<Scalar, dimWorld>;
- using DimVector = Dune::FieldVector<Scalar, dim>;
- using Tensor = Dune::FieldMatrix<Scalar, dim, dim>;
- using LocalIndexType = std::uint8_t;
- using LocalIndexSet = std::vector<LocalIndexType>;
- using LocalIndexPair = std::pair<LocalIndexType, bool>;
- using GlobalIndexType = typename GridView::IndexSet::IndexType;
- using GlobalIndexSet = std::vector<GlobalIndexType>;
+public:
+ using BoundaryInteractionVolume = CCMpfaInteractionVolumeImplementation<TypeTag, MpfaMethods::oMethod>;
+ using PositionVector = std::vector<GlobalPosition>;
using Matrix = Dune::DynamicMatrix<Scalar>;
using Vector = typename Matrix::row_type;
- using PositionVector = std::vector<GlobalPosition>;
-
using LocalScvType = CCMpfaOLocalScv<TypeTag>;
using LocalScvfType = CCMpfaOLocalScvf<TypeTag>;
-
- using Seed = CCMpfaOInteractionVolumeSeed<GlobalIndexType, LocalIndexType>;
+ using typename BaseTraits::LocalIndexSet;
+ using typename BaseTraits::GlobalIndexSet;
+ using Seed = CCMpfaOInteractionVolumeSeed<GlobalIndexSet, LocalIndexSet, dim, dimWorld>;
};
//! Forward declaration of the mpfa-o interaction volume
@@ -101,7 +94,8 @@ public:
template<class TypeTag, class Traits>
class CCMpfaOInteractionVolume : public CCMpfaInteractionVolumeBase<TypeTag, Traits>
{
- // The interaction volume implementation has to be friend
+ // The interaction volume implementation has to be friend,
+ // because some methods use the mpfa-o interaction volume as base
friend typename GET_PROP_TYPE(TypeTag, InteractionVolume);
using ParentType = CCMpfaInteractionVolumeBase<TypeTag, Traits>;
@@ -114,10 +108,11 @@ class CCMpfaOInteractionVolume : public CCMpfaInteractionVolumeBase<TypeTag, Tra
using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);
static const int dim = GridView::dimension;
+ static const int dimWorld = GridView::dimensionworld;
using Element = typename GridView::template Codim<0>::Entity;
+ using DimVector = Dune::FieldVector<Scalar, dim>;
+ using GlobalPosition = Dune::FieldVector<Scalar, dimWorld>;
- using DimVector = typename Traits::DimVector;
- using GlobalPosition = typename Traits::GlobalPosition;
using DynamicVector = typename Traits::Vector;
using DynamicMatrix = typename Traits::Matrix;
using Tensor = typename Traits::Tensor;
@@ -128,7 +123,7 @@ class CCMpfaOInteractionVolume : public CCMpfaInteractionVolumeBase<TypeTag, Tra
public:
using typename ParentType::LocalIndexType;
using typename ParentType::LocalIndexSet;
- using typename ParentType::LocalIndexPair;
+ using typename ParentType::LocalFaceData;
using typename ParentType::GlobalIndexSet;
using typename ParentType::PositionVector;
using typename ParentType::Seed;
@@ -238,7 +233,11 @@ public:
}
}
- LocalIndexPair getLocalIndexPair(const SubControlVolumeFace& scvf) const
+ //! gets local data on a global scvf within the interaction volume
+ //! specialization for dim = dimWorld
+ template<int d = dim, int dw = dimWorld>
+ typename std::enable_if< (d == dw), LocalFaceData >::type
+ getLocalFaceData(const SubControlVolumeFace& scvf) const
{
auto scvfGlobalIdx = scvf.index();
@@ -246,10 +245,10 @@ public:
for (const auto& localScvf : localScvfs_)
{
if (localScvf.insideGlobalScvfIndex() == scvfGlobalIdx)
- return std::make_pair(localIdx, false);
+ return LocalFaceData(localIdx, localScvf.insideLocalScvIndex(), false);
if (!localScvf.boundary() && localScvf.outsideGlobalScvfIndex() == scvfGlobalIdx)
- return std::make_pair(localIdx, true);
+ return LocalFaceData(localIdx, localScvf.outsideLocalScvIndex(), true);
localIdx++;
}
@@ -257,22 +256,115 @@ public:
DUNE_THROW(Dune::InvalidStateException, "Could not find the local scv face in the interaction volume for the given scvf with index: " << scvf.index());
}
- DynamicVector getTransmissibilities(const LocalIndexPair& localIndexPair) const
+ //! gets local data on a global scvf within the interaction volume
+ //! specialization for dim < dimWorld
+ template<int d = dim, int dw = dimWorld>
+ typename std::enable_if< (d < dw), LocalFaceData >::type
+ getLocalFaceData(const SubControlVolumeFace& scvf) const
{
- auto tij = T_[localIndexPair.first];
+ const auto scvfGlobalIdx = scvf.index();
- if (localIndexPair.second)
+ LocalIndexType localFaceIdx = 0;
+ for (const auto& localScvf : localScvfs_)
+ {
+ if (localScvf.insideGlobalScvfIndex() == scvfGlobalIdx)
+ return LocalFaceData(localFaceIdx, localScvf.insideLocalScvIndex(), false);
+
+ if (!localScvf.boundary() && MpfaHelper::contains(localScvf.outsideGlobalScvfIndices(), scvfGlobalIdx))
+ {
+ if (localScvf.outsideLocalScvIndices().size() == 1)
+ return LocalFaceData(localFaceIdx, localScvf.outsideLocalScvIndex(), true);
+
+ // Now we have to find wich local scv is the one the global scvf is embedded in
+ const auto scvGlobalIdx = scvf.insideScvIdx();
+ for (auto localScvIdx : localScvf.outsideLocalScvIndices())
+ if (localScv_(localScvIdx).globalIndex() == scvGlobalIdx)
+ return LocalFaceData(localFaceIdx, localScvIdx, true);
+
+ DUNE_THROW(Dune::InvalidStateException, "Could not find the local scv in the interaction volume for the given scvf with index: " << scvfGlobalIdx);
+ }
+
+ localFaceIdx++;
+ }
+
+ DUNE_THROW(Dune::InvalidStateException, "Could not find the local scv face in the interaction volume for the given scvf with index: " << scvfGlobalIdx);
+ }
+
+ //! gets the transmissibilities for a sub-control volume face within the interaction volume
+ //! specialization for dim == dimWorld
+ template<int d = dim, int dw = dimWorld>
+ typename std::enable_if< (d == dw), DynamicVector >::type
+ getTransmissibilities(const LocalFaceData& localFaceData) const
+ {
+ auto tij = T_[localFaceData.localScvfIndex];
+ if (localFaceData.isOutside)
tij *= -1.0;
return tij;
}
- Scalar getNeumannFlux(const LocalIndexPair& localIndexPair) const
+ //! gets the transmissibilities for a sub-control volume face within the interaction volume
+ //! specialization for dim < dimWorld
+ template<int d = dim, int dw = dimWorld>
+ typename std::enable_if< (d < dw), DynamicVector >::type
+ getTransmissibilities(const LocalFaceData& localFaceData) const
{
- if (!onBoundary() || fluxScvfIndexSet_().size() == 0 || GET_PROP_VALUE(TypeTag, UseTpfaBoundary))
+ // if we are not on a branching point, do the usual
+ if (localScvf_(localFaceData.localScvfIndex).outsideLocalScvIndices().size() == 1)
+ {
+ auto tij = T_[localFaceData.localScvfIndex];
+ if (localFaceData.isOutside)
+ tij *= -1.0;
+ return tij;
+ }
+
+ // This means we are on a branching point. If we come from the inside, simply return tij
+ if (!localFaceData.isOutside)
+ return T_[localFaceData.localScvfIndex];
+
+ // otherwise compute the outside transmissibilities
+ DynamicVector tij(volVarsStencil().size(), 0.0);
+
+ // get the local scv and iterate over local coordinates
+ const std::size_t numLocalScvs = localScvs_.size();
+ const auto& localScv = localScv_(localFaceData.localScvIndex);
+ const auto& localScvf = localScvf_(localFaceData.localScvfIndex);
+
+ auto idxInOutside = this->findLocalIndex(localScvf.outsideLocalScvIndices(), localFaceData.localScvIndex);
+ const auto& wijk = wijk_[localFaceData.localScvfIndex][idxInOutside+1];
+ for (int localDir = 0; localDir < dim; localDir++)
+ {
+ auto localScvfIdx = localScv.localScvfIndex(localDir);
+ const auto& localScvf = localScvf_(localScvfIdx);
+
+ // add entries from the face unknowns
+ if (localScvf.faceType() != MpfaFaceTypes::dirichlet)
+ {
+ auto fluxFaceIndex = this->findLocalIndex(fluxScvfIndexSet_(), localScvfIdx);
+ auto tmp = AinvB_[fluxFaceIndex];
+ tmp *= wijk[localDir];
+
+ tij += tmp;
+ }
+ else
+ {
+ auto idxInDiriFaces = this->findLocalIndex(dirichletScvfIndexSet_(), localScvfIdx);
+ tij[numLocalScvs + idxInDiriFaces] += wijk[localDir];
+ }
+
+ // add entry from the scv unknown
+ tij[localFaceData.localScvIndex] -= wijk[localDir];
+ }
+
+ return tij;
+ }
+
+ Scalar getNeumannFlux(const LocalFaceData& localFaceData) const
+ {
+ if (!onBoundary() || GET_PROP_VALUE(TypeTag, UseTpfaBoundary) || fluxScvfIndexSet_().size() == 0 )
return 0.0;
- auto flux = CAinv_[localIndexPair.first] * neumannFluxes_;
- if (localIndexPair.second)
+ auto flux = CAinv_[localFaceData.localScvfIndex] * neumannFluxes_;
+ if (localFaceData.isOutside)
return -1.0*flux;
return flux;
}
@@ -339,6 +431,9 @@ private:
{
const std::size_t numLocalScvs = localScvs_.size();
+ // reserve space for the omegas
+ wijk_.resize(localScvfs_.size());
+
// loop over the local faces
LocalIndexType rowIdx = 0;
for (const auto& localScvf : localScvfs_)
@@ -391,40 +486,49 @@ private:
B[idxInFluxFaces][posLocalScvIdx] += posWijk[localDir];
}
+ // store the omegas
+ wijk_[rowIdx].emplace_back(std::move(posWijk));
+
// If face is not on a boundary or an interior dirichlet face, add entries for the "negative" scv
if (faceType == MpfaFaceTypes::interior)
{
- const auto negLocalScvIdx = localScvf.outsideLocalScvIndex();
- const auto& negLocalScv = localScv_(negLocalScvIdx);
- const auto& negGlobalScv = fvGeometry_().scv(negLocalScv.globalIndex());
- auto negElement = localElement_(negLocalScvIdx);;
- auto negTensor = getTensor(negElement, elemVolVars_()[negGlobalScv], negGlobalScv);
-
- // the omega factors of the "negative" sub volume
- auto negWijk = calculateOmegas_(negLocalScv, localScvf, negTensor);
- negWijk *= problem_().boxExtrusionFactor(negElement, negGlobalScv);
-
- // Check local directions of negative sub volume
- for (int localDir = 0; localDir < dim; localDir++)
+ // loop over all the outside neighbors of this face and add entries
+ for (auto negLocalScvIdx : localScvf.outsideLocalScvIndices())
{
- auto curLocalScvfIdx = negLocalScv.localScvfIndex(localDir);
- const auto& curLocalScvf = localScvf_(curLocalScvfIdx);
+ const auto& negLocalScv = localScv_(negLocalScvIdx);
+ const auto& negGlobalScv = fvGeometry_().scv(negLocalScv.globalIndex());
+ auto negElement = localElement_(negLocalScvIdx);;
+ auto negTensor = getTensor(negElement, elemVolVars_()[negGlobalScv], negGlobalScv);
- if (curLocalScvf.faceType() != MpfaFaceTypes::dirichlet)
- {
- // we need the index of the current local scvf in the flux face indices
- auto curIdxInFluxFaces = this->findLocalIndex(fluxScvfIndexSet_(), curLocalScvfIdx);
- A[idxInFluxFaces][curIdxInFluxFaces] -= negWijk[localDir];
- }
- else
+ // the omega factors of the "negative" sub volume
+ auto negWijk = calculateOmegas_(negLocalScv, localScvf, negTensor);
+ negWijk *= problem_().boxExtrusionFactor(negElement, negGlobalScv);
+
+ // Check local directions of negative sub volume
+ for (int localDir = 0; localDir < dim; localDir++)
{
- // the current face is a Dirichlet face and creates entries in B
- auto curIdxInDiriFaces = this->findLocalIndex(dirichletScvfIndexSet_(), curLocalScvfIdx);
- B[idxInFluxFaces][numLocalScvs + curIdxInDiriFaces] += negWijk[localDir];
+ auto curLocalScvfIdx = negLocalScv.localScvfIndex(localDir);
+ const auto& curLocalScvf = localScvf_(curLocalScvfIdx);
+
+ if (curLocalScvf.faceType() != MpfaFaceTypes::dirichlet)
+ {
+ // we need the index of the current local scvf in the flux face indices
+ auto curIdxInFluxFaces = this->findLocalIndex(fluxScvfIndexSet_(), curLocalScvfIdx);
+ A[idxInFluxFaces][curIdxInFluxFaces] -= negWijk[localDir];
+ }
+ else
+ {
+ // the current face is a Dirichlet face and creates entries in B
+ auto curIdxInDiriFaces = this->findLocalIndex(dirichletScvfIndexSet_(), curLocalScvfIdx);
+ B[idxInFluxFaces][numLocalScvs + curIdxInDiriFaces] += negWijk[localDir];
+ }
+
+ // add entries to matrix B
+ B[idxInFluxFaces][negLocalScvIdx] -= negWijk[localDir];
}
- // add entries to matrix B
- B[idxInFluxFaces][negLocalScvIdx] -= negWijk[localDir];
+ // store the omegas
+ wijk_[rowIdx].emplace_back(std::move(negWijk));
}
}
// go to the next face
@@ -444,6 +548,9 @@ private:
AinvB_.resize(0, 0);
CAinv_.resize(0, 0);
+ // resize the omegas
+ wijk_.resize(numFaces);
+
// Loop over all the faces, in this case these are all dirichlet boundaries
LocalIndexType rowIdx = 0;
for (const auto& localScvf : localScvfs_)
@@ -468,6 +575,9 @@ private:
T_[rowIdx][posLocalScvIdx] -= posWijk[localDir];
}
+ // store the omegas
+ wijk_[rowIdx].emplace_back(std::move(posWijk));
+
// go to the next face
rowIdx++;
}
@@ -477,7 +587,7 @@ private:
const LocalScvfType& localScvf,
const Tensor& T) const
{
- GlobalPosition wijk;
+ DimVector wijk;
GlobalPosition tmp;
for (int dir = 0; dir < dim; ++dir)
{
@@ -495,7 +605,7 @@ private:
const LocalScvfType& localScvf,
const Scalar t) const
{
- GlobalPosition wijk;
+ DimVector wijk;
GlobalPosition tmp(localScvf.unitOuterNormal());
tmp *= t;
@@ -522,7 +632,6 @@ private:
const ElementVolumeVariables* elemVolVarsPtr_;
bool onBoundary_;
- LocalIndexType eqIdx_;
std::vector<Element> localElements_;
std::vector<LocalScvType> localScvs_;
@@ -535,6 +644,8 @@ private:
LocalIndexSet fluxFaceIndexSet_;
LocalIndexSet dirichletFaceIndexSet_;
+ std::vector< std::vector< DimVector > > wijk_;
+
DynamicMatrix T_;
DynamicMatrix AinvB_;
DynamicMatrix CAinv_;
diff --git a/dumux/discretization/cellcentered/mpfa/omethod/interactionvolumeseed.hh b/dumux/discretization/cellcentered/mpfa/omethod/interactionvolumeseed.hh
index 0387ae55bc66e1415824cae0c479fe867e46ed97..bdcdc831fc42678eb57b8b7eca62ac20b8689857 100644
--- a/dumux/discretization/cellcentered/mpfa/omethod/interactionvolumeseed.hh
+++ b/dumux/discretization/cellcentered/mpfa/omethod/interactionvolumeseed.hh
@@ -31,12 +31,13 @@ namespace Dumux
/*!
* \ingroup Mpfa
- * \brief Base class for the interaction volume seed of mpfa methods
+ * \brief Class for the interaction volume seed of the mpfa-o method
*/
-template<typename G, typename L>
+template<typename G, typename L, int dim, int dimWorld>
class CCMpfaOInteractionVolumeSeed
{
- using GlobalIndexType = G;
+ using GlobalIndexSet = G;
+ using GlobalIndexType = typename GlobalIndexSet::value_type;
public:
using LocalScvSeed = CCMpfaOLocalScvSeed<G, L>;
@@ -47,8 +48,7 @@ public:
bool onBoundary)
: onBoundary_(onBoundary),
scvSeeds_(std::move(scvSeeds)),
- scvfSeeds_(std::move(scvfSeeds))
- {}
+ scvfSeeds_(std::move(scvfSeeds)) {}
bool onBoundary() const
{ return onBoundary_; }
@@ -76,8 +76,11 @@ public:
globalIndices.reserve(scvfSeeds().size() * 2);
for (auto&& localScvfSeed : scvfSeeds())
- for (auto scvfIdxGlobal : localScvfSeed.globalScvfIndices())
+ {
+ globalIndices.push_back(localScvfSeed.insideGlobalScvfIndex());
+ for (auto scvfIdxGlobal : localScvfSeed.outsideGlobalScvfIndices())
globalIndices.push_back(scvfIdxGlobal);
+ }
globalIndices.shrink_to_fit();
return globalIndices;
diff --git a/dumux/discretization/cellcentered/mpfa/omethod/localsubcontrolentities.hh b/dumux/discretization/cellcentered/mpfa/omethod/localsubcontrolentities.hh
index 1e78961d0b3ddcd18445bcccef35ce507d42e19d..ce4046ae6e6afebaff9ebc1fefc1d98f6078f497 100644
--- a/dumux/discretization/cellcentered/mpfa/omethod/localsubcontrolentities.hh
+++ b/dumux/discretization/cellcentered/mpfa/omethod/localsubcontrolentities.hh
@@ -36,17 +36,18 @@ class CCMpfaOLocalScv
using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
using Helper = typename GET_PROP_TYPE(TypeTag, MpfaHelper);
using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
- using InteractionVolume = typename GET_PROP_TYPE(TypeTag, InteractionVolume);
using Element = typename GridView::template Codim<0>::Entity;
- using LocalScvSeed = typename InteractionVolume::Seed::LocalScvSeed;
- using GlobalIndexType = typename GridView::IndexSet::IndexType;
+ // we use the seed types of the boundary interaction volume to be compatible with other mpfa
+ // methods that use o-type interaction volumes on the boundary but differing ones inside the domain.
+ using InteractionVolume = typename GET_PROP_TYPE(TypeTag, BoundaryInteractionVolume);
+ using LocalScvSeed = typename InteractionVolume::Seed::LocalScvSeed;
+ using GlobalIndexType = typename InteractionVolume::GlobalIndexType;
using LocalIndexType = typename InteractionVolume::LocalIndexType;
static const int dim = GridView::dimension;
static const int dimWorld = GridView::dimensionworld;
-
using GlobalPosition = Dune::FieldVector<Scalar, dimWorld>;
using LocalBasis = std::array<GlobalPosition, dim>;
@@ -113,15 +114,19 @@ struct CCMpfaOLocalScvf
using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
using Helper = typename GET_PROP_TYPE(TypeTag, MpfaHelper);
- using InteractionVolume = typename GET_PROP_TYPE(TypeTag, InteractionVolume);
using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
using Element = typename GridView::template Codim<0>::Entity;
- using LocalScvfSeed = typename InteractionVolume::Seed::LocalScvfSeed;
- using GlobalIndexType = typename GridView::IndexSet::IndexType;
+ // we use the seed types of the boundary interaction volume to be compatible with other mpfa
+ // methods that use o-type interaction volumes on the boundary but differing ones inside the domain.
+ using InteractionVolume = typename GET_PROP_TYPE(TypeTag, BoundaryInteractionVolume);
+ using LocalScvfSeed = typename InteractionVolume::Seed::LocalScvfSeed;
+ using GlobalIndexType = typename InteractionVolume::GlobalIndexType;
+ using GlobalIndexSet = typename InteractionVolume::GlobalIndexSet;
using LocalIndexType = typename InteractionVolume::LocalIndexType;
+ using LocalIndexSet = typename InteractionVolume::LocalIndexSet;
static const int dim = GridView::dimension;
static const int dimWorld = GridView::dimensionworld;
@@ -131,36 +136,50 @@ public:
CCMpfaOLocalScvf(const LocalScvfSeed& scvfSeed,
const SubControlVolumeFace& scvf)
: seedPtr_(&scvfSeed),
- center_(scvf.center()),
ip_(scvf.ipGlobal()),
normal_(scvf.unitOuterNormal()),
area_(scvf.area())
{}
- GlobalIndexType insideGlobalScvfIndex() const
+ const GlobalIndexType insideGlobalScvfIndex() const
{ return scvfSeed_().insideGlobalScvfIndex(); }
- GlobalIndexType outsideGlobalScvfIndex() const
- { return scvfSeed_().outsideGlobalScvfIndex(); }
+ const GlobalIndexSet& outsideGlobalScvfIndices() const
+ { return scvfSeed_().outsideGlobalScvfIndices(); }
- LocalIndexType insideLocalScvIndex() const
+ const GlobalIndexType outsideGlobalScvfIndex() const
+ {
+ assert(scvfSeed_().outsideGlobalScvfIndices().size() == 1 && "outside scvf index not uniquely defined");
+ return scvfSeed_().outsideGlobalScvfIndices()[0];
+ }
+
+ const LocalIndexType insideLocalScvIndex() const
{ return scvfSeed_().insideLocalScvIndex(); }
- LocalIndexType outsideLocalScvIndex() const
- { return scvfSeed_().outsideLocalScvIndex(); }
+ const LocalIndexSet& outsideLocalScvIndices() const
+ { return scvfSeed_().outsideLocalScvIndices(); }
+
+ const LocalIndexType outsideLocalScvIndex() const
+ {
+ assert(scvfSeed_().outsideLocalScvIndices().size() == 1 && "outside local scv index not uniquely defined");
+ return scvfSeed_().outsideLocalScvIndices()[0];
+ }
- GlobalIndexType insideGlobalScvIndex() const
+ const GlobalIndexType insideGlobalScvIndex() const
{ return scvfSeed_().insideGlobalScvIndex(); }
- GlobalIndexType outsideGlobalScvIndex() const
- { return scvfSeed_().outsideGlobalScvIndex();}
+ const GlobalIndexSet& outsideGlobalScvIndices() const
+ { return scvfSeed_().outsideGlobalScvIndices(); }
+
+ const GlobalIndexType outsideGlobalScvIndex() const
+ {
+ assert(scvfSeed_().outsideGlobalScvIndices().size() == 1 && "outside scv index not uniquely defined");
+ return scvfSeed_().outsideGlobalScvIndices()[0];
+ }
MpfaFaceTypes faceType() const
{ return scvfSeed_().faceType(); }
- GlobalPosition center() const
- { return center_; }
-
GlobalPosition ip() const
{ return ip_; }
@@ -178,7 +197,6 @@ private:
{ return *seedPtr_; }
const LocalScvfSeed* seedPtr_;
- GlobalPosition center_;
GlobalPosition ip_;
GlobalPosition normal_;
Scalar area_;
diff --git a/dumux/discretization/cellcentered/mpfa/omethod/localsubcontrolentityseeds.hh b/dumux/discretization/cellcentered/mpfa/omethod/localsubcontrolentityseeds.hh
index add3fdca34f91efeb85d1eb7de424a32621f50c7..49fbc1ad2b5e99c4d87108e7f1a6ca87b71243bf 100644
--- a/dumux/discretization/cellcentered/mpfa/omethod/localsubcontrolentityseeds.hh
+++ b/dumux/discretization/cellcentered/mpfa/omethod/localsubcontrolentityseeds.hh
@@ -32,10 +32,10 @@ template<typename G, typename L>
class CCMpfaOLocalScvSeed
{
public:
- using GlobalIndexType = G;
- using LocalIndexType = L;
- using GlobalIndexSet = std::vector<GlobalIndexType>;
- using LocalIndexSet = std::vector<LocalIndexType>;
+ using GlobalIndexSet = G;
+ using LocalIndexSet = L;
+ using GlobalIndexType = typename GlobalIndexSet::value_type;
+ using LocalIndexType = typename LocalIndexSet::value_type;
//! constructor fully defining the scv seed
CCMpfaOLocalScvSeed(GlobalIndexSet&& globalScvfIndices,
@@ -43,8 +43,7 @@ public:
GlobalIndexType globalScvIndex)
: globalScvIndex_(globalScvIndex),
localScvfIndices_(std::move(localScvfIndices)),
- globalScvfIndices_(std::move(globalScvfIndices))
- {}
+ globalScvfIndices_(std::move(globalScvfIndices)) {}
//! Constructor when the local scvf indices are not known at this point
CCMpfaOLocalScvSeed(GlobalIndexSet&& globalScvfIndices,
@@ -75,60 +74,78 @@ private:
GlobalIndexSet globalScvfIndices_;
};
-
+/*!
+ * \ingroup Mpfa
+ * \brief Class for a sub-control volume face seed of the mpfa-o method.
+ *
+ * \param G the global index set type
+ * \param L the local index set type
+ */
template<typename G, typename L>
class CCMpfaOLocalScvfSeed
{
public:
- using GlobalIndexType = G;
- using LocalIndexType = L;
- using GlobalIndexSet = std::vector<GlobalIndexType>;
- using LocalIndexSet = std::vector<LocalIndexType>;
+ using GlobalIndexSet = G;
+ using LocalIndexSet = L;
+ using GlobalIndexType = typename GlobalIndexSet::value_type;
+ using LocalIndexType = typename LocalIndexSet::value_type;
+ //! constructor fully defining the scv face seed
+ template<class SubControlVolumeFace>
+ CCMpfaOLocalScvfSeed(const SubControlVolumeFace& scvf,
+ LocalIndexType insideLocalScvIndex,
+ LocalIndexSet&& outsideLocalScvIndices,
+ GlobalIndexSet&& outsideGlobalScvfIndices,
+ const MpfaFaceTypes faceType)
+ : boundary_(scvf.boundary()),
+ insideScvLocalIdx_(insideLocalScvIndex),
+ outsideLocalScvIndices_(std::move(outsideLocalScvIndices)),
+ insideScvGlobalIdx_(scvf.insideScvIdx()),
+ outsideGlobalScvIndices_(scvf.outsideScvIndices()),
+ insideScvfGlobalIdx_(scvf.index()),
+ outsideGlobalScvfIndices_(std::move(outsideGlobalScvfIndices)),
+ faceType_(faceType) {}
+
+ //! Constructor when the outside indices are not known at this point
template<class SubControlVolumeFace>
CCMpfaOLocalScvfSeed(const SubControlVolumeFace& scvf,
- LocalIndexSet&& scvIndicesLocal,
- GlobalIndexSet&& scvfIndicesGlobal,
+ LocalIndexType insideLocalScvIndex,
const MpfaFaceTypes faceType)
: boundary_(scvf.boundary()),
- scvIndicesGlobal_({scvf.insideScvIdx(), scvf.outsideScvIdx()}),
- scvfIndicesGlobal_(std::move(scvfIndicesGlobal)),
- scvIndicesLocal_(std::move(scvIndicesLocal)),
+ insideScvLocalIdx_(insideLocalScvIndex),
+ insideScvGlobalIdx_(scvf.insideScvIdx()),
+ outsideGlobalScvIndices_(scvf.outsideScvIndices()),
+ insideScvfGlobalIdx_(scvf.index()),
faceType_(faceType)
- {}
-
- const GlobalIndexSet& globalScvfIndices() const
- { return scvfIndicesGlobal_; }
-
- const LocalIndexSet& localScvIndices() const
- { return scvIndicesLocal_; }
+ {
+ auto size = outsideGlobalScvIndices_.size();
+ outsideLocalScvIndices_.reserve(size);
+ outsideGlobalScvfIndices_.reserve(size);
+ }
- const GlobalIndexSet& globalScvIndices() const
- { return scvIndicesGlobal_; }
+ const GlobalIndexType insideGlobalScvfIndex() const
+ { return insideScvfGlobalIdx_; }
- LocalIndexType insideLocalScvIndex() const
- { return scvIndicesLocal_[0]; }
+ const GlobalIndexSet& outsideGlobalScvfIndices() const
+ { return outsideGlobalScvfIndices_; }
- LocalIndexType outsideLocalScvIndex() const
+ const GlobalIndexType outsideGlobalScvfIndex() const
{
- assert(!boundary() && "There is no outside local scv for boundary scvfs");
- return scvIndicesLocal_[1];
+ assert(outsideGlobalScvfIndices_.size() == 1 && "outside global scvf index not uniquely defined");
+ return outsideGlobalScvfIndices_[0];
}
- GlobalIndexType insideGlobalScvIndex() const
- { return scvIndicesGlobal_[0]; }
+ const LocalIndexType insideLocalScvIndex() const
+ { return insideScvLocalIdx_; }
- GlobalIndexType outsideGlobalScvIndex() const
- { return scvIndicesGlobal_[1]; }
+ const LocalIndexSet& outsideLocalScvIndices() const
+ { return outsideLocalScvIndices_; }
- GlobalIndexType insideGlobalScvfIndex() const
- { return scvfIndicesGlobal_[0]; }
+ const GlobalIndexType insideGlobalScvIndex() const
+ { return insideScvGlobalIdx_; }
- GlobalIndexType outsideGlobalScvfIndex() const
- {
- assert(!boundary() && "There is no outside global scvf for boundary scvfs");
- return scvfIndicesGlobal_[1];
- }
+ const GlobalIndexSet& outsideGlobalScvIndices() const
+ { return outsideGlobalScvIndices_; }
MpfaFaceTypes faceType() const
{ return faceType_; }
@@ -136,11 +153,37 @@ public:
bool boundary() const
{ return boundary_; }
+ void addOutsideScvfIndex(const GlobalIndexType index)
+ {
+ outsideGlobalScvfIndices_.push_back(index);
+ }
+
+ void addOutsideLocalScvIndex(const LocalIndexType index)
+ {
+ outsideLocalScvIndices_.push_back(index);
+ }
+
+ void makeOutsideDataUnique()
+ {
+ std::sort(outsideGlobalScvfIndices_.begin(), outsideGlobalScvfIndices_.end());
+ outsideGlobalScvfIndices_.erase(std::unique(outsideGlobalScvfIndices_.begin(), outsideGlobalScvfIndices_.end()), outsideGlobalScvfIndices_.end());
+
+ std::sort(outsideLocalScvIndices_.begin(), outsideLocalScvIndices_.end());
+ outsideLocalScvIndices_.erase(std::unique(outsideLocalScvIndices_.begin(), outsideLocalScvIndices_.end()), outsideLocalScvIndices_.end());
+ }
+
private:
bool boundary_;
- GlobalIndexSet scvIndicesGlobal_;
- GlobalIndexSet scvfIndicesGlobal_;
- LocalIndexSet scvIndicesLocal_;
+
+ LocalIndexType insideScvLocalIdx_;
+ LocalIndexSet outsideLocalScvIndices_;
+
+ GlobalIndexType insideScvGlobalIdx_;
+ GlobalIndexSet outsideGlobalScvIndices_;
+
+ GlobalIndexType insideScvfGlobalIdx_;
+ GlobalIndexSet outsideGlobalScvfIndices_;
+
MpfaFaceTypes faceType_;
};
} // end namespace
diff --git a/dumux/discretization/cellcentered/mpfa/omethod/subcontrolvolumeface.hh b/dumux/discretization/cellcentered/mpfa/omethod/subcontrolvolumeface.hh
index 5bc232954fcf6f677d402b67d1782482d3ae8bbf..992f3364269731a39cea040e4fd27df62d4a7bc4 100644
--- a/dumux/discretization/cellcentered/mpfa/omethod/subcontrolvolumeface.hh
+++ b/dumux/discretization/cellcentered/mpfa/omethod/subcontrolvolumeface.hh
@@ -47,7 +47,7 @@ class CCMpfaSubControlVolumeFace<MpfaMethods::oMethod, G, I> : public CCMpfaSubC
using GlobalPosition = Dune::FieldVector<Scalar, dimworld>;
public:
- //! We do not use the vertexIndexInElement variable.
+ //! We do not use the localIndex variable.
//! It is here to satisfy the general mpfa scvf interface.
template<class MpfaGeometryHelper>
CCMpfaSubControlVolumeFace(const MpfaGeometryHelper& geomHelper,
@@ -56,19 +56,20 @@ public:
IndexType vertexIndex,
unsigned int localIndex,
IndexType scvfIndex,
- std::array<IndexType, 2>&& scvIndices,
+ IndexType insideScvIdx,
+ const std::vector<IndexType>& outsideScvIndices,
Scalar q,
bool boundary)
: ParentType(geomHelper,
- std::move(corners),
- std::move(unitOuterNormal),
+ std::forward<std::vector<GlobalPosition>>(corners),
+ std::forward<GlobalPosition>(unitOuterNormal),
vertexIndex,
localIndex,
scvfIndex,
- std::move(scvIndices),
+ insideScvIdx,
+ outsideScvIndices,
q,
- boundary)
- {}
+ boundary) {}
};
} // end namespace
diff --git a/dumux/discretization/cellcentered/mpfa/omethodfps/globalfvgeometry.hh b/dumux/discretization/cellcentered/mpfa/omethodfps/globalfvgeometry.hh
deleted file mode 100644
index d7a01c290242302191e5d28e49487ec0b9fd890a..0000000000000000000000000000000000000000
--- a/dumux/discretization/cellcentered/mpfa/omethodfps/globalfvgeometry.hh
+++ /dev/null
@@ -1,46 +0,0 @@
-// -*- 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_FPS_GLOBALFVGEOMETRY_HH
-#define DUMUX_DISCRETIZATION_CC_MPFA_O_FPS_GLOBALFVGEOMETRY_HH
-
-#include <dumux/discretization/cellcentered/mpfa/globalfvgeometrybase.hh>
-#include <dumux/discretization/cellcentered/mpfa/methods.hh>
-
-namespace Dumux
-{
-//! Specialization for the mpfa-o scheme with full pressure support.
-template<class TypeTag, bool EnableCache>
-class CCMpfaGlobalFVGeometryImplementation<TypeTag, MpfaMethods::oMethodFps, EnableCache> : public CCMpfaGlobalFVGeometryBase<TypeTag, EnableCache>
-{
- using ParentType = CCMpfaGlobalFVGeometryBase<TypeTag, EnableCache>;
- using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
-
-public:
- CCMpfaGlobalFVGeometryImplementation(const GridView gridView) : ParentType(gridView) {}
-};
-
-} // end namespace
-
-#endif
diff --git a/dumux/discretization/cellcentered/mpfa/omethodfps/helper.hh b/dumux/discretization/cellcentered/mpfa/omethodfps/helper.hh
index 3a00caf30fe06b5d025d05813921cc104b25370f..d0e6d3faa15dfac346e93162a9d7ee1b0aca231d 100644
--- a/dumux/discretization/cellcentered/mpfa/omethodfps/helper.hh
+++ b/dumux/discretization/cellcentered/mpfa/omethodfps/helper.hh
@@ -32,8 +32,8 @@ namespace Dumux
* \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::oMethodFps, dim> : public MpfaHelperBase<TypeTag, MpfaMethods::oMethod, dim>
+template<class TypeTag, int dim, int dimWorld>
+class MpfaMethodHelper<TypeTag, MpfaMethods::oMethodFps, dim, dimWorld> : public MpfaMethodHelper<TypeTag, MpfaMethods::oMethod, dim, dimWorld>
{};
} // end namespace
diff --git a/dumux/discretization/cellcentered/mpfa/omethodfps/interactionvolume.hh b/dumux/discretization/cellcentered/mpfa/omethodfps/interactionvolume.hh
index 62fc0bc6ca112937068976dd0ef059a9f18d7c79..d5483f18d34b53dfd6fc1fc5ab13de2b2eb7641c 100644
--- a/dumux/discretization/cellcentered/mpfa/omethodfps/interactionvolume.hh
+++ b/dumux/discretization/cellcentered/mpfa/omethodfps/interactionvolume.hh
@@ -67,6 +67,7 @@ class CCMpfaInteractionVolumeImplementation<TypeTag, MpfaMethods::oMethodFps> :
using LocalScvfType = typename Traits::LocalScvfType;
static const int dim = GridView::dimension;
+ static const int dimWorld = GridView::dimensionworld;
using CoordScalar = typename GridView::ctype;
using FeCache = Dune::PQkLocalFiniteElementCache<CoordScalar, Scalar, dim, 1>;
using FeLocalBasis = typename FeCache::FiniteElementType::Traits::LocalBasisType;
@@ -74,8 +75,8 @@ class CCMpfaInteractionVolumeImplementation<TypeTag, MpfaMethods::oMethodFps> :
using ShapeValue = typename Dune::FieldVector<Scalar, 1>;
using JacobianInverseTransposed = typename LocalScvType::Geometry::JacobianInverseTransposed;
- using GlobalPosition = typename Traits::GlobalPosition;
using LocalPosition = typename LocalScvType::Geometry::LocalCoordinate;
+ using GlobalPosition = Dune::FieldVector<Scalar, dimWorld>;
using DynamicVector = typename Traits::Vector;
using DynamicMatrix = typename Traits::Matrix;
using Tensor = typename Traits::Tensor;
@@ -101,7 +102,7 @@ class CCMpfaInteractionVolumeImplementation<TypeTag, MpfaMethods::oMethodFps> :
public:
using typename ParentType::LocalIndexType;
- using typename ParentType::LocalIndexPair;
+ using typename ParentType::LocalFaceData;
using typename ParentType::Seed;
CCMpfaInteractionVolumeImplementation(const Seed& seed,
@@ -234,34 +235,36 @@ private:
LocalMatrixContainer& mc)
{
// get diffusion tensor in "negative" sub volume
- const auto localScvIdx = localScvf.outsideLocalScvIndex();
- const auto& localScv = this->localScv_(localScvIdx);
- const auto& globalScv = this->fvGeometry_().scv(localScv.globalIndex());
- const auto& element = this->localElement_(localScvIdx);;
- auto D = makeTensor_(getTensor(element, this->elemVolVars_()[globalScv], globalScv));
+ for (auto localScvIdx : localScvf.outsideLocalScvIndices())
+ {
+ const auto& localScv = this->localScv_(localScvIdx);
+ const auto& globalScv = this->fvGeometry_().scv(localScv.globalIndex());
+ const auto& element = this->localElement_(localScvIdx);;
+ auto D = makeTensor_(getTensor(element, this->elemVolVars_()[globalScv], globalScv));
- // the local finite element bases of the scvs
- const auto& localBasis = feCache_.get(localScv.geometry().type()).localBasis();
+ // the local finite element bases of the scvs
+ const auto& localBasis = feCache_.get(localScv.geometry().type()).localBasis();
- // the normal and local integration point
- // On the ref element, normal vector points in the direction of local coordinate
- auto normalDir = localScv.getScvfIdxInScv(localScvfIdx);
- auto ipLocal = localScv.geometry().local(localScvf.ip());
+ // the normal and local integration point
+ // On the ref element, normal vector points in the direction of local coordinate
+ auto normalDir = localScv.getScvfIdxInScv(localScvfIdx);
+ auto ipLocal = localScv.geometry().local(localScvf.ip());
- // find normals and integration points in the two scvs for condition of zero divergence
- LocalIndexType divEqNormalDir = normalDir == 1 ? 0 : 1;
- LocalPosition divEqIpLocal(0.0);
- divEqIpLocal[divEqNormalDir] = divEqNormalDir == 1 ? c_ : 1.0 - (1.0-c_)*p_;
- divEqIpLocal[normalDir] = divEqNormalDir == 1 ? c_ + (1.0-c_)*p_ : c_;
+ // find normals and integration points in the two scvs for condition of zero divergence
+ LocalIndexType divEqNormalDir = normalDir == 1 ? 0 : 1;
+ LocalPosition divEqIpLocal(0.0);
+ divEqIpLocal[divEqNormalDir] = divEqNormalDir == 1 ? c_ : 1.0 - (1.0-c_)*p_;
+ divEqIpLocal[normalDir] = divEqNormalDir == 1 ? c_ + (1.0-c_)*p_ : c_;
- // does the face has an unknown associated with it?
- bool isFluxFace = localScvf.faceType() != MpfaFaceTypes::dirichlet;
+ // does the face has an unknown associated with it?
+ bool isFluxFace = localScvf.faceType() != MpfaFaceTypes::dirichlet;
- // assemble coefficients for the face fluxes
- addFaceFluxCoefficients_(localScv, localBasis, D, localScvfIdx, ipLocal, normalDir, mc, isFluxFace, true);
+ // assemble coefficients for the face fluxes
+ addFaceFluxCoefficients_(localScv, localBasis, D, localScvfIdx, ipLocal, normalDir, mc, isFluxFace, true);
- // assemble matrix entries for the condition of zero divergence
- addDivEquationCoefficients_(localScv, localBasis, D, divEqIpLocal, divEqNormalDir, mc);
+ // assemble matrix entries for the condition of zero divergence
+ addDivEquationCoefficients_(localScv, localBasis, D, divEqIpLocal, divEqNormalDir, mc);
+ }
}
void addFaceFluxCoefficients_(const LocalScvType& localScv,
diff --git a/dumux/discretization/cellcentered/mpfa/omethodfps/mpfafpsgeometryhelper.hh b/dumux/discretization/cellcentered/mpfa/omethodfps/mpfafpsgeometryhelper.hh
index b993c7e22806172ba4fb2fd64622ba6719b1ced8..ba427928048b4a7ad2a7d3bad68111971d9e9055 100644
--- a/dumux/discretization/cellcentered/mpfa/omethodfps/mpfafpsgeometryhelper.hh
+++ b/dumux/discretization/cellcentered/mpfa/omethodfps/mpfafpsgeometryhelper.hh
@@ -34,11 +34,11 @@
namespace Dumux
{
-//! Create sub control volumes and sub control volume face geometries
+//! A class to mainly extract the scv corners of the scvs inside an element
template<class GridView, int dim>
class MpfaFpsGeometryHelper;
-//! A class to create sub control volume and sub control volume face geometries per element
+//! Specialization for dim = 2
template <class GridView>
class MpfaFpsGeometryHelper<GridView, 2>
{
@@ -134,121 +134,6 @@ private:
std::size_t corners_; // number of element corners
};
-//! A class to create sub control volume and sub control volume face geometries per element
-template <class GridView>
-class MpfaFpsGeometryHelper<GridView, 3>
-{
- using Scalar = typename GridView::ctype;
- static const int dim = GridView::dimension;
- static const int dimWorld = GridView::dimensionworld;
-
- using ScvGeometry = Dune::CachedMultiLinearGeometry<Scalar, dim, dimWorld>;
- using ScvfGeometry = Dune::CachedMultiLinearGeometry<Scalar, dim-1, dimWorld>;
-
- using GlobalPosition = typename ScvGeometry::GlobalCoordinate;
- using PointVector = std::vector<GlobalPosition>;
-
- using Element = typename GridView::template Codim<0>::Entity;
- using Intersection = typename GridView::Intersection;
-
- using ReferenceElements = typename Dune::ReferenceElements<Scalar, dim>;
- using FaceReferenceElements = typename Dune::ReferenceElements<Scalar, dim-1>;
-
- //! the maximum number of helper points used to construct the geometries
- //! Using a statically sized point array is much faster than dynamic allocation
- static constexpr int maxPoints = 27;
-public:
- MpfaFpsGeometryHelper(const typename Element::Geometry& geometry)
- : elementGeometry_(geometry), corners_(geometry.corners())
- {
- // extract the corners of the sub control volumes
- const auto& referenceElement = ReferenceElements::general(geometry.type());
-
- // the element center
- p[0] = geometry.center();
-
- // vertices
- for (int i = 0; i < corners_; ++i)
- p[i+1] = geometry.corner(i);
-
- // edge midpoints
- for (int i = 0; i < referenceElement.size(dim-1); ++i)
- p[i+corners_+1] = geometry.global(referenceElement.position(i, dim-1));
-
- // face midpoints
- for (int i = 0; i < referenceElement.size(1); ++i)
- p[i+corners_+1+referenceElement.size(dim-1)] = geometry.global(referenceElement.position(i, 1));
- }
-
- //! Create a vector with the scv corners
- PointVector getScvCorners(unsigned int localScvIdx) const
- {
- // proceed according to number of corners of the element
- switch (corners_)
- {
- case 4: // tetrahedron
- {
- //! Only build the maps the first time we encounter a tetrahedron
- static const std::uint8_t vo = 1; //! vertex offset in point vector p
- static const std::uint8_t eo = 5; //! edge offset in point vector p
- static const std::uint8_t fo = 11; //! face offset in point vector p
- static const std::uint8_t map[4][8] =
- {
- {vo+0, eo+0, eo+1, fo+0, eo+3, fo+1, fo+2, 0},
- {vo+1, eo+2, eo+0, fo+0, eo+4, fo+3, fo+1, 0},
- {vo+2, eo+1, eo+2, fo+0, eo+5, fo+2, fo+3, 0},
- {vo+3, eo+3, eo+5, fo+2, eo+4, fo+1, fo+3, 0}
- };
-
- return PointVector( {p[map[localScvIdx][0]],
- p[map[localScvIdx][1]],
- p[map[localScvIdx][2]],
- p[map[localScvIdx][3]],
- p[map[localScvIdx][4]],
- p[map[localScvIdx][5]],
- p[map[localScvIdx][6]],
- p[map[localScvIdx][7]]} );
- }
- case 8: // hexahedron
- {
- //! Only build the maps the first time we encounter a quadrilateral
- static const std::uint8_t vo = 1; //! vertex offset in point vector p
- static const std::uint8_t eo = 9; //! edge offset in point vector p
- static const std::uint8_t fo = 21; //! face offset in point vector p
- static const std::uint8_t map[8][8] =
- {
- {vo+0, eo+6, eo+4, fo+4, eo+0, fo+2, fo+0, 0},
- {vo+1, eo+5, eo+6, fo+4, eo+1, fo+1, fo+2, 0},
- {vo+2, eo+4, eo+7, fo+4, eo+2, fo+0, fo+3, 0},
- {vo+3, eo+7, eo+5, fo+4, eo+3, fo+3, fo+1, 0},
- {vo+4, eo+8, eo+10, fo+5, eo+0, fo+0, fo+2, 0},
- {vo+5, eo+10, eo+9, fo+5, eo+1, fo+2, fo+1, 0},
- {vo+6, eo+11, eo+8, fo+5, eo+2, fo+3, fo+0, 0},
- {vo+7, eo+9, eo+11, fo+5, eo+3, fo+1, fo+3, 0},
- };
-
- return PointVector( {p[map[localScvIdx][0]],
- p[map[localScvIdx][1]],
- p[map[localScvIdx][2]],
- p[map[localScvIdx][3]],
- p[map[localScvIdx][4]],
- p[map[localScvIdx][5]],
- p[map[localScvIdx][6]],
- p[map[localScvIdx][7]]} );
- }
- default:
- DUNE_THROW(Dune::NotImplemented, "Mpfa-Fps scv corners for dim=" << dim
- << " dimWorld=" << dimWorld
- << " corners=" << corners_);
- }
- }
-
-private:
- const typename Element::Geometry& elementGeometry_; //! Reference to the element geometry
- GlobalPosition p[maxPoints]; // the points needed for construction of the scv/scvf geometries
- std::size_t corners_; // number of element corners
-};
-
} // end namespace Dumux
#endif
diff --git a/dumux/discretization/cellcentered/mpfa/omethodfps/subcontrolvolumeface.hh b/dumux/discretization/cellcentered/mpfa/omethodfps/subcontrolvolumeface.hh
index f98927a6d6f6864abdd0d1fdf445f50036549f28..cffb6a13cca3bd57db01622c84c4de7125b55a64 100644
--- a/dumux/discretization/cellcentered/mpfa/omethodfps/subcontrolvolumeface.hh
+++ b/dumux/discretization/cellcentered/mpfa/omethodfps/subcontrolvolumeface.hh
@@ -24,7 +24,7 @@
#define DUMUX_DISCRETIZATION_CC_MPFA_O_FPS_SUBCONTROLVOLUMEFACE_HH
#include <dumux/discretization/cellcentered/mpfa/methods.hh>
-#include <dumux/discretization/cellcentered/mpfa/omethod/subcontrolvolumeface.hh>
+#include <dumux/discretization/cellcentered/mpfa/subcontrolvolumefacebase.hh>
namespace Dumux
{
@@ -34,9 +34,9 @@ namespace Dumux
* \brief Class for a sub control volume face in the mpfao-fps method.
*/
template<class G, typename I>
-class CCMpfaSubControlVolumeFace<MpfaMethods::oMethodFps, G, I> : public CCMpfaSubControlVolumeFace<MpfaMethods::oMethod, G, I>
+class CCMpfaSubControlVolumeFace<MpfaMethods::oMethodFps, G, I> : public CCMpfaSubControlVolumeFaceBase<G, I>
{
- using ParentType = CCMpfaSubControlVolumeFace<MpfaMethods::oMethod, G, I>;
+ using ParentType = CCMpfaSubControlVolumeFaceBase<G, I>;
using Geometry = G;
using IndexType = I;
@@ -47,23 +47,25 @@ class CCMpfaSubControlVolumeFace<MpfaMethods::oMethodFps, G, I> : public CCMpfaS
using GlobalPosition = Dune::FieldVector<Scalar, dimworld>;
public:
- template<class MpfaGeometryHelper>
- CCMpfaSubControlVolumeFace(const MpfaGeometryHelper& geomHelper,
+ template<class MpfaHelper>
+ CCMpfaSubControlVolumeFace(const MpfaHelper& helper,
std::vector<GlobalPosition>&& corners,
GlobalPosition&& unitOuterNormal,
IndexType vIdxGlobal,
unsigned int vIdxLocal,
IndexType scvfIndex,
- std::array<IndexType, 2>&& scvIndices,
+ IndexType insideScvIdx,
+ const std::vector<IndexType>& outsideScvIndices,
Scalar q,
bool boundary)
- : ParentType(geomHelper,
+ : ParentType(helper,
std::move(corners),
std::move(unitOuterNormal),
vIdxGlobal,
vIdxLocal,
scvfIndex,
- std::move(scvIndices),
+ insideScvIdx,
+ outsideScvIndices,
q,
boundary),
vIdxInElement_(vIdxLocal)
diff --git a/dumux/discretization/cellcentered/mpfa/omethod/globalinteractionvolumeseeds.hh b/dumux/discretization/cellcentered/mpfa/stencils.hh
similarity index 63%
rename from dumux/discretization/cellcentered/mpfa/omethod/globalinteractionvolumeseeds.hh
rename to dumux/discretization/cellcentered/mpfa/stencils.hh
index 9ee05b9332fe9ab30eadca7884166befbb6c32db..734951864ea814c656c4bfad1d40b752c0e4e6a7 100644
--- a/dumux/discretization/cellcentered/mpfa/omethod/globalinteractionvolumeseeds.hh
+++ b/dumux/discretization/cellcentered/mpfa/stencils.hh
@@ -18,31 +18,24 @@
*****************************************************************************/
/*!
* \file
- * \brief Base class for the global interaction volume seeds of mpfa methods.
+ * \brief Implements the notion of stencils for cell-centered models
*/
-#ifndef DUMUX_DISCRETIZATION_MPFA_O_GLOBALINTERACTIONVOLUMESEEDS_HH
-#define DUMUX_DISCRETIZATION_MPFA_O_GLOBALINTERACTIONVOLUMESEEDS_HH
+#ifndef DUMUX_DISCRETIZATION_CC_MPFA_STENCILS_HH
+#define DUMUX_DISCRETIZATION_CC_MPFA_STENCILS_HH
-#include <dumux/discretization/cellcentered/mpfa/globalinteractionvolumeseedsbase.hh>
-#include <dumux/discretization/cellcentered/mpfa/methods.hh>
+#include <dumux/discretization/cellcentered/stencils.hh>
namespace Dumux
{
-/*!
- * \ingroup Mpfa
- * \brief Specialization of the class for the mpfa-o method.
- */
-template<class TypeTag>
-class CCMpfaGlobalInteractionVolumeSeedsImplementation<TypeTag, MpfaMethods::oMethod> : public CCMpfaGlobalInteractionVolumeSeedsBase<TypeTag>
-{
- using ParentType = CCMpfaGlobalInteractionVolumeSeedsBase<TypeTag>;
- using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+//! Method-specific implementation of the stencils vector
+//! By default, we inherit from the symmetric stencils vector class
+//! For non-symmetric mpfa methods a corresponding specialization has to be provided below
+template<class TypeTag, MpfaMethods method>
+class CCMpfaStencilsVectorImplementation : public CCSymmetricStencilsVector<TypeTag> {};
-public:
- CCMpfaGlobalInteractionVolumeSeedsImplementation(const GridView gridView) : ParentType(gridView) {}
-};
+template<class TypeTag>
+using CCMpfaStencilsVector = CCMpfaStencilsVectorImplementation<TypeTag, GET_PROP_VALUE(TypeTag, MpfaMethod)>;
} // end namespace
-
#endif
diff --git a/dumux/discretization/cellcentered/mpfa/subcontrolvolumeface.hh b/dumux/discretization/cellcentered/mpfa/subcontrolvolumeface.hh
index 947df913a9dc69cbbe830be1d28613d369f4d2da..5d1f96917818b2e6cf6d2182e143b350158e367c 100644
--- a/dumux/discretization/cellcentered/mpfa/subcontrolvolumeface.hh
+++ b/dumux/discretization/cellcentered/mpfa/subcontrolvolumeface.hh
@@ -53,11 +53,10 @@ class CCMpfaSubControlVolumeFace : public CCMpfaSubControlVolumeFaceImplementati
public:
/*!
- * \brief Constructor interface of an mpfa sub control volume face
+ * \brief Constructor
*
- * We define a general constructor for scv faces in mpfa methods.
- * This is necessary so that the class instantiating scv faces can work
- * method-independent.
+ * We do not use the localIndex here. Its meaning can vary depending on the
+ * implementation (i.e. mpfa method) and is handled by the implementation itself.
*
* \param geomHelper The mpfa geometry helper
* \param corners The corners of the scv face
@@ -65,26 +64,38 @@ public:
* \param vIdxGlobal The global vertex index the scvf is connected to
* \param localIndex Some element local index (the local vertex index in mpfao-fps)
* \param scvfIndex The global index of this scv face
- * \param scvIndices The inside and outside scv indices connected to this face
+ * \param insideScvIdx The inside scv index connected to this face
+ * \param outsideScvIndices The outside scv indices connected to this face
* \param q The parameterization of the quadrature point on the scvf for flux calculation
* \param boundary Boolean to specify whether or not the scvf is on a boundary
*/
- template<class MpfaGeometryHelper>
- CCMpfaSubControlVolumeFace(const MpfaGeometryHelper& geomHelper,
+ template<class MpfaHelper>
+ CCMpfaSubControlVolumeFace(const MpfaHelper& helper,
std::vector<GlobalPosition>&& corners,
GlobalPosition&& unitOuterNormal,
IndexType vIdxGlobal,
unsigned int localIndex,
IndexType scvfIndex,
- std::array<IndexType, 2>&& scvIndices,
+ IndexType insideScvIdx,
+ const std::vector<IndexType>& outsideScvIndices,
Scalar q,
bool boundary)
- : ParentType(geomHelper, corners, unitOuterNormal, vIdxGlobal, localIndex, scvfIndex, scvIndices, q, boundary)
+ : ParentType(helper,
+ std::forward<std::vector<GlobalPosition>>(corners),
+ std::forward<GlobalPosition>(unitOuterNormal),
+ vIdxGlobal,
+ localIndex,
+ scvfIndex,
+ insideScvIdx,
+ outsideScvIndices,
+ q,
+ boundary)
{}
};
} // end namespace
//! The available implementations should be included here
+#include <dumux/discretization/cellcentered/mpfa/lmethod/subcontrolvolumeface.hh>
#include <dumux/discretization/cellcentered/mpfa/omethod/subcontrolvolumeface.hh>
#include <dumux/discretization/cellcentered/mpfa/omethodfps/subcontrolvolumeface.hh>
diff --git a/dumux/discretization/cellcentered/mpfa/subcontrolvolumefacebase.hh b/dumux/discretization/cellcentered/mpfa/subcontrolvolumefacebase.hh
index 57dc6591b31ffb5867e5a16833bcb0101c6dcce0..ca4598fd8ee1c20fb2fa9fe9e93178611fbabd41 100644
--- a/dumux/discretization/cellcentered/mpfa/subcontrolvolumefacebase.hh
+++ b/dumux/discretization/cellcentered/mpfa/subcontrolvolumefacebase.hh
@@ -60,25 +60,28 @@ public:
* \param vIdxGlobal The global vertex index the scvf is connected to
* \param localIndex Some element local index (the local vertex index in mpfao-fps)
* \param scvfIndex The global index of this scv face
- * \param scvIndices The inside and outside scv indices connected to this face
+ * \param insideScvIdx The inside scv index connected to this face
+ * \param outsideScvIndices The outside scv indices connected to this face
* \param q The parameterization of the quadrature point on the scvf for flux calculation
* \param boundary Boolean to specify whether or not the scvf is on a boundary
*/
- template<class MpfaGeometryHelper>
- CCMpfaSubControlVolumeFaceBase(const MpfaGeometryHelper& geomHelper,
+ template<class MpfaHelper>
+ CCMpfaSubControlVolumeFaceBase(const MpfaHelper& helper,
std::vector<GlobalPosition>&& corners,
GlobalPosition&& unitOuterNormal,
IndexType vertexIndex,
unsigned int localIndex,
IndexType scvfIndex,
- std::array<IndexType, 2>&& scvIndices,
+ IndexType insideScvIdx,
+ const std::vector<IndexType>& outsideScvIndices,
Scalar q,
bool boundary)
: ParentType(),
boundary_(boundary),
vertexIndex_(vertexIndex),
scvfIndex_(scvfIndex),
- scvIndices_(std::move(scvIndices)),
+ insideScvIdx_(insideScvIdx),
+ outsideScvIndices_(outsideScvIndices),
corners_(std::move(corners)),
center_(0.0),
unitOuterNormal_(std::move(unitOuterNormal))
@@ -86,8 +89,8 @@ public:
for (const auto& corner : corners_)
center_ += corner;
center_ /= corners_.size();
- ipGlobal_ = geomHelper.getScvfIntegrationPoint(corners_, q);
- area_ = geomHelper.getScvfArea(corners_);
+ ipGlobal_ = helper.getScvfIntegrationPoint(corners_, q);
+ area_ = helper.getScvfArea(corners_);
}
//! The center of the sub control volume face
@@ -111,11 +114,20 @@ public:
//! index of the inside sub control volume for spatial param evaluation
IndexType insideScvIdx() const
- { return scvIndices_[0]; }
+ { return insideScvIdx_; }
//! index of the outside sub control volume for spatial param evaluation
+ //! returns in undefined behaviour is boundary is true
+ //! is not uniquely defined for grids where dim < dimWorld and shouldn't be used in that case.
IndexType outsideScvIdx() const
- { return scvIndices_[1]; }
+ {
+ assert(outsideScvIndices_.size() == 1 && "outside scv index not uniquely defined");
+ return outsideScvIndices_[0];
+ }
+
+ //! returns the outside scv indices (can be more than one index for dim < dimWorld)
+ std::vector<IndexType> outsideScvIndices() const
+ { return outsideScvIndices_; }
//! The global index of this sub control volume face
IndexType index() const
@@ -152,7 +164,8 @@ private:
bool boundary_;
IndexType vertexIndex_;
IndexType scvfIndex_;
- std::array<IndexType, 2> scvIndices_;
+ IndexType insideScvIdx_;
+ std::vector<IndexType> outsideScvIndices_;
std::vector<GlobalPosition> corners_;
GlobalPosition center_;
diff --git a/dumux/discretization/cellcentered/mpfa/omethodfps/globalinteractionvolumeseeds.hh b/dumux/implicit/cellcentered/mpfa/assemblymap.hh
similarity index 61%
rename from dumux/discretization/cellcentered/mpfa/omethodfps/globalinteractionvolumeseeds.hh
rename to dumux/implicit/cellcentered/mpfa/assemblymap.hh
index cd1ce2fb5f7b1fb97682fae7a1fafa5c3780dc3e..cf9b374fba30f5b29447b9a43a31bcd151e5011d 100644
--- a/dumux/discretization/cellcentered/mpfa/omethodfps/globalinteractionvolumeseeds.hh
+++ b/dumux/implicit/cellcentered/mpfa/assemblymap.hh
@@ -18,31 +18,27 @@
*****************************************************************************/
/*!
* \file
- * \brief Base class for the global interaction volume seeds of mpfa methods.
+ * \brief Stores the face indices corresponding to the neighbors of an element
+ * that contribute to the derivative calculation
*/
-#ifndef DUMUX_DISCRETIZATION_MPFA_O_FPS_GLOBALINTERACTIONVOLUMESEEDS_HH
-#define DUMUX_DISCRETIZATION_MPFA_O_FPS_GLOBALINTERACTIONVOLUMESEEDS_HH
+#ifndef DUMUX_CC_MPFA_ASSEMBLY_MAP_HH
+#define DUMUX_CC_MPFA_ASSEMBLY_MAP_HH
-#include <dumux/discretization/cellcentered/mpfa/globalinteractionvolumeseedsbase.hh>
-#include <dumux/discretization/cellcentered/mpfa/methods.hh>
+#include <dumux/implicit/cellcentered/assemblymap.hh>
namespace Dumux
{
-/*!
- * \ingroup Mpfa
- * \brief Specialization of the class for the mpfa-o method with full pressure support.
- */
-template<class TypeTag>
-class CCMpfaGlobalInteractionVolumeSeedsImplementation<TypeTag, MpfaMethods::oMethodFps> : public CCMpfaGlobalInteractionVolumeSeedsBase<TypeTag>
-{
- using ParentType = CCMpfaGlobalInteractionVolumeSeedsBase<TypeTag>;
- using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
-
-public:
- CCMpfaGlobalInteractionVolumeSeedsImplementation(const GridView gridView) : ParentType(gridView) {}
-};
+//! Forward declaration of method specific implementation of the assembly map
+template<class TypeTag, MpfaMethods method>
+class CCMpfaAssemblyMapImplementation;
-} // end namespace
+// //! The Assembly map for models using mpfa methods
+template<class TypeTag>
+using CCMpfaAssemblyMap = CCMpfaAssemblyMapImplementation<TypeTag, GET_PROP_VALUE(TypeTag, MpfaMethod)>;
+//! The default is simply the CCAssemblyMap
+template<class TypeTag, MpfaMethods method>
+class CCMpfaAssemblyMapImplementation : public CCAssemblyMap<TypeTag> {};
+}
#endif
diff --git a/dumux/implicit/cellcentered/mpfa/propertydefaults.hh b/dumux/implicit/cellcentered/mpfa/propertydefaults.hh
index fda0c3b10a71233bb9da282cad35d12cbac11748..9ff21ee167e7f90fa00b8933d457bc31eea491d6 100644
--- a/dumux/implicit/cellcentered/mpfa/propertydefaults.hh
+++ b/dumux/implicit/cellcentered/mpfa/propertydefaults.hh
@@ -39,6 +39,7 @@
#include <dumux/discretization/cellcentered/mpfa/elementfluxvariablescache.hh>
#include <dumux/discretization/cellcentered/mpfa/subcontrolvolumeface.hh>
#include <dumux/discretization/cellcentered/mpfa/helper.hh>
+#include <dumux/discretization/cellcentered/mpfa/stencils.hh>
#include <dumux/discretization/cellcentered/mpfa/interactionvolume.hh>
#include <dumux/discretization/cellcentered/mpfa/globalinteractionvolumeseeds.hh>
#include <dumux/implicit/cellcentered/mpfa/localresidual.hh>
@@ -46,9 +47,6 @@
namespace Dumux {
-// forward declarations
-template<class TypeTag> class CCElementBoundaryTypes;
-
namespace Properties {
//! Set the corresponding discretization method property
SET_PROP(CCMpfaModel, DiscretizationMethod)
@@ -56,7 +54,7 @@ SET_PROP(CCMpfaModel, DiscretizationMethod)
static const DiscretizationMethods value = DiscretizationMethods::CCMpfa;
};
-//! Set the BaseLocalResidual to CCLocalResidual
+//! Set the BaseLocalResidual to CCMpfaLocalResidual
SET_TYPE_PROP(CCMpfaModel, BaseLocalResidual, CCMpfaLocalResidual<TypeTag>);
//! By default we set the o-method as the Mpfa method of choice
@@ -66,7 +64,7 @@ SET_PROP(CCMpfaModel, MpfaMethod)
};
//! The mpfa helper class
-SET_TYPE_PROP(CCMpfaModel, MpfaHelper, MpfaHelper<TypeTag>);
+SET_TYPE_PROP(CCMpfaModel, MpfaHelper, CCMpfaHelper<TypeTag>);
//! The interaction volume class
SET_TYPE_PROP(CCMpfaModel, InteractionVolume, CCMpfaInteractionVolume<TypeTag>);
@@ -95,6 +93,9 @@ SET_TYPE_PROP(CCMpfaModel, ElementVolumeVariables, Dumux::CCMpfaElementVolumeVar
//! The local flux variables cache vector class
SET_TYPE_PROP(CCMpfaModel, ElementFluxVariablesCache, Dumux::CCMpfaElementFluxVariablesCache<TypeTag, GET_PROP_VALUE(TypeTag, EnableGlobalFluxVariablesCache)>);
+//! The stencil container, some mpfa methods might lead to non-symmetric global systems
+SET_TYPE_PROP(CCMpfaModel, StencilsVector, Dumux::CCMpfaStencilsVector<TypeTag>);
+
SET_PROP(CCMpfaModel, SubControlVolumeFace)
{
private:
diff --git a/dumux/porousmediumflow/implicit/fluxvariablescache.hh b/dumux/porousmediumflow/implicit/fluxvariablescache.hh
index 3e7de05f2ea6703e84a35d8fb17d97a6850854d5..6d6725745629b843350a39b5559d8b7b629b19c2 100644
--- a/dumux/porousmediumflow/implicit/fluxvariablescache.hh
+++ b/dumux/porousmediumflow/implicit/fluxvariablescache.hh
@@ -152,13 +152,13 @@ private:
// forward declaration of the base class of the mpfa flux variables cache
template<class TypeTag, bool EnableAdvection, bool EnableMolecularDiffusion, bool EnableEnergyBalance>
-class PorousMediumMpfaFluxVariablesCache
+class MpfaPorousMediumFluxVariablesCache
{};
// specialization for cell centered mpfa methods
template<class TypeTag>
class PorousMediumFluxVariablesCacheImplementation<TypeTag, DiscretizationMethods::CCMpfa>
- : public PorousMediumMpfaFluxVariablesCache<TypeTag,
+ : public MpfaPorousMediumFluxVariablesCache<TypeTag,
GET_PROP_VALUE(TypeTag, EnableAdvection),
GET_PROP_VALUE(TypeTag, EnableMolecularDiffusion),
GET_PROP_VALUE(TypeTag, EnableEnergyBalance)>
@@ -166,7 +166,7 @@ class PorousMediumFluxVariablesCacheImplementation<TypeTag, DiscretizationMethod
// specialization for the case of pure advection
template<class TypeTag>
-class PorousMediumMpfaFluxVariablesCache<TypeTag, true, false, false>
+class MpfaPorousMediumFluxVariablesCache<TypeTag, true, false, false>
{
using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
@@ -192,11 +192,11 @@ class PorousMediumMpfaFluxVariablesCache<TypeTag, true, false, false>
public:
// the constructor
- PorousMediumMpfaFluxVariablesCache() : isUpdated_(false)
+ MpfaPorousMediumFluxVariablesCache() : isUpdated_(false)
{
- // We have to initialize to zero (for inner interaction volumes)
- for (auto& nFlux : phaseNeumannFluxes_)
- nFlux = 0.0;
+ // We have to initialize the neumann fluxes to zero (for inner interaction volumes)
+ for (auto& flux : phaseNeumannFluxes_)
+ flux = 0.0;
}
// update cached objects
@@ -208,19 +208,10 @@ public:
const SubControlVolumeFace &scvf,
const InteractionVolume& interactionVolume)
{
- const auto& localIndexPair = interactionVolume.getLocalIndexPair(scvf);
- const auto& volVarsStencil = interactionVolume.volVarsStencil();
- const auto& volVarsPositions = interactionVolume.volVarsPositions();
-
- // the types coming from the inner interaction volumes might differ (thus, = assignment is not possible)
- const auto numVolVars = volVarsStencil.size();
- volVarsStencil_.clear();
- volVarsStencil_.reserve(numVolVars);
- volVarsPositions_.clear();
- volVarsPositions_.reserve(numVolVars);
- volVarsStencil_.insert(volVarsStencil_.begin(), volVarsStencil.begin(), volVarsStencil.end());
- volVarsPositions_.insert(volVarsPositions_.begin(), volVarsPositions.begin(), volVarsPositions.end());
- tij_ = interactionVolume.getTransmissibilities(localIndexPair);
+ const auto& localFaceData = interactionVolume.getLocalFaceData(scvf);
+ volVarsStencil_ = interactionVolume.volVarsStencil();
+ volVarsPositions_ = interactionVolume.volVarsPositions();
+ tij_ = interactionVolume.getTransmissibilities(localFaceData);
}
template<typename InteractionVolume>
@@ -232,8 +223,8 @@ public:
const InteractionVolume& interactionVolume,
const unsigned int phaseIdx)
{
- const auto& localIndexPair = interactionVolume.getLocalIndexPair(scvf);
- phaseNeumannFluxes_[phaseIdx] = interactionVolume.getNeumannFlux(localIndexPair);
+ const auto& localFaceData = interactionVolume.getLocalFaceData(scvf);
+ phaseNeumannFluxes_[phaseIdx] = interactionVolume.getNeumannFlux(localFaceData);
}
const Stencil& advectionVolVarsStencil(const unsigned int phaseIdx) const
@@ -251,14 +242,9 @@ public:
bool isUpdated() const
{ return isUpdated_; }
- void setUpdated()
+ void setUpdateStatus(const bool status)
{
- isUpdated_ = true;
- }
-
- void setOutdated()
- {
- isUpdated_ = false;
+ isUpdated_ = status;
}
private: