diff --git a/dumux/discretization/cellcentered/mpfa/fvelementgeometry.hh b/dumux/discretization/cellcentered/mpfa/fvelementgeometry.hh
index e79a41db5eca80335665f9585cb48af9f002a811..3be824293ef4db33eb2fb33d195895e6209847e9 100644
--- a/dumux/discretization/cellcentered/mpfa/fvelementgeometry.hh
+++ b/dumux/discretization/cellcentered/mpfa/fvelementgeometry.hh
@@ -261,18 +261,19 @@ public:
// reserve memory
const auto numNeighbors = assemblyMapI.size();
- const auto estimate = neighborScvfEstimate_(element, numNeighbors);
+ const auto numNeighborScvfs = numNeighborScvfs_(assemblyMapI);
neighborScvs_.reserve(numNeighbors);
neighborScvIndices_.reserve(numNeighbors);
- neighborScvfIndices_.reserve(estimate);
- neighborScvfs_.reserve(estimate);
+ neighborScvfIndices_.reserve(numNeighborScvfs);
+ neighborScvfs_.reserve(numNeighborScvfs);
// make neighbor geometries
// use the assembly map to determine which faces are necessary
for (auto&& dataJ : assemblyMapI)
makeNeighborGeometries(globalFvGeometry().element(dataJ.globalJ),
dataJ.globalJ,
- dataJ.scvfsJ);
+ dataJ.scvfsJ,
+ dataJ.additionalScvfs);
// set up the scvf flip indices for network grids
if (dim < dimWorld)
@@ -309,31 +310,13 @@ private:
//! Estimates the number of neighboring scvfs that have to be prepared
//! The estimate holds for inner elements, overestimats on the boundary
- int neighborScvfEstimate_(const Element& element, int numNeighbors)
+ template<class DataJ>
+ unsigned int numNeighborScvfs_(const std::vector<DataJ>& dataJVector)
{
- auto type = element.geometry().type();
-
- if (type == Dune::GeometryType(Dune::GeometryType::simplex, dim))
- {
- if (dim == 2)
- return 6 + ( (numNeighbors-3 > 0) ? (numNeighbors-3)*2 : 0 );
- if (dim == 3)
- return 12 + ( (numNeighbors-4 > 0) ? (numNeighbors-4)*3 : 0 );
- }
- else if (type == Dune::GeometryType(Dune::GeometryType::cube, dim))
- {
- if (dim == 2)
- return 16;
- if (dim == 3)
- return 120;
- }
- else if (type == Dune::GeometryType(Dune::GeometryType::pyramid, dim))
- return 16 + ( (numNeighbors-4 > 0) ? (numNeighbors-4)*3 : 0 );
- else if (type == Dune::GeometryType(Dune::GeometryType::prism, dim))
- return 18 + ( (numNeighbors-4 > 0) ? (numNeighbors-4)*3 : 0 );
- else
- DUNE_THROW(Dune::NotImplemented, "Mpfa for given geometry type.");
-
+ unsigned int numNeighborScvfs = 0;
+ for (const auto& dataJ : dataJVector)
+ numNeighborScvfs += dataJ.scvfsJ.size() + dataJ.additionalScvfs.size();
+ return numNeighborScvfs;
}
//! create scvs and scvfs of the bound element
@@ -424,7 +407,10 @@ private:
//! create the scv and necessary scvfs of a neighboring element
template<typename IndexVector>
- void makeNeighborGeometries(const Element& element, IndexType eIdxGlobal, const IndexVector& scvfIndices)
+ void makeNeighborGeometries(const Element& element,
+ IndexType eIdxGlobal,
+ const IndexVector& scvfIndices,
+ const IndexVector& additionalScvfs)
{
// create the neighbor scv if it doesn't exist yet
neighborScvs_.emplace_back(element.geometry(), eIdxGlobal);
@@ -483,7 +469,8 @@ private:
continue;
// only build the scvf if it is in the list of necessary indices
- if (MpfaHelper::contains(scvfIndices, scvFaceIndices[scvfCounter]))
+ if (MpfaHelper::contains(scvfIndices, scvFaceIndices[scvfCounter])
+ || MpfaHelper::contains(additionalScvfs, scvFaceIndices[scvfCounter]))
{
neighborScvfs_.emplace_back(MpfaHelper(),
MpfaHelper::getScvfCorners(isCorners, c),
diff --git a/dumux/discretization/cellcentered/mpfa/interactionvolumebase.hh b/dumux/discretization/cellcentered/mpfa/interactionvolumebase.hh
index 32af446ad2a96ed62aeadc991d5cd6a2222e5ba8..144e7881fb81cad5e8b12ebd5c2195c0be1c2683 100644
--- a/dumux/discretization/cellcentered/mpfa/interactionvolumebase.hh
+++ b/dumux/discretization/cellcentered/mpfa/interactionvolumebase.hh
@@ -91,6 +91,8 @@ public:
isOutside(isOut) {}
};
+ using GlobalLocalFaceDataPair = std::pair<const SubControlVolumeFace*, LocalFaceData>;
+
//! solves the local equation system for the computation of the transmissibilities
template<typename GetTensorFunction>
void solveLocalSystem(const GetTensorFunction& getTensor)
@@ -104,10 +106,6 @@ public:
const PositionVector& volVarsPositions() const
{ DUNE_THROW(Dune::NotImplemented, "Actual interaction volume implementation does not provide a volVarsPositions() method."); }
- //! returns a list of global scvf indices that are connected to this interaction volume
- const GlobalIndexSet& globalScvfs() const
- { 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)
LocalFaceData getLocalFaceData(const SubControlVolumeFace& scvf) const
{ DUNE_THROW(Dune::NotImplemented, "Actual interaction volume implementation does not provide a getLocalFaceData() method."); }
@@ -128,6 +126,10 @@ public:
assert(it != vector.end() && "could not find local index in the vector for the given global index!");
return std::distance(vector.begin(), it);
}
+
+ //! returns GlobalLocalFaceDataPair objects for the scvfs involved in this interaction volume
+ const std::vector<GlobalLocalFaceDataPair>& globalLocalScvfPairedData() const
+ { DUNE_THROW(Dune::NotImplemented, "Actual interaction volume implementation does not provide a globalLocalScvfPairedData() method."); }
};
} // end namespace
diff --git a/dumux/discretization/cellcentered/mpfa/lmethod/globalinteractionvolumeseeds.hh b/dumux/discretization/cellcentered/mpfa/lmethod/globalinteractionvolumeseeds.hh
index 8a3da8af4f4bdc17d2e75722b525f3dc270592e9..d9fe0f67c7a7c67b03f904fe58100e5f5c151d66 100644
--- a/dumux/discretization/cellcentered/mpfa/lmethod/globalinteractionvolumeseeds.hh
+++ b/dumux/discretization/cellcentered/mpfa/lmethod/globalinteractionvolumeseeds.hh
@@ -141,12 +141,7 @@ public:
}
bool isLocalMaxLevel_(const Element& element, const SubControlVolumeFace& scvf) const
- {
- auto inLevel = element.level();
- if (this->problem().model().globalFvGeometry().element(scvf.outsideScvIdx()).level() > inLevel)
- return false;
- return true;
- }
+ { return this->problem().model().globalFvGeometry().element(scvf.outsideScvIdx()).level() <= element.level(); }
};
} // end namespace
diff --git a/dumux/discretization/cellcentered/mpfa/lmethod/helper.hh b/dumux/discretization/cellcentered/mpfa/lmethod/helper.hh
index f5d88a68972c485f9b8111d617c90e6b3f744612..1ca199a5d2ae2d837d6a01611d03dbae7f7fa480 100644
--- a/dumux/discretization/cellcentered/mpfa/lmethod/helper.hh
+++ b/dumux/discretization/cellcentered/mpfa/lmethod/helper.hh
@@ -98,7 +98,6 @@ public:
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];
@@ -116,7 +115,7 @@ public:
Scalar s1 = std::abs(t11-t12);
Scalar s2 = std::abs(t22-t21);
- if (s1 < s2 + eps*s1)
+ if (s1 < s2)
return 0;
else
return 1;
@@ -132,16 +131,16 @@ private:
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);
+ const auto scvfVector = Implementation::getScvFacesAtVertex(scvf.vertexIndex(), element, fvGeometry);
+ const auto localScvfIdx = Implementation::getLocalFaceIndex(scvf, scvfVector);
scvSeeds.emplace_back( GlobalIndexSet({scvfVector[0]->index(), scvfVector[1]->index()}),
scvf.insideScvIdx(),
localScvfIdx );
// get the surrounding elements and "outside" data
LocalIndexType otherScvfIdx = 1-localScvfIdx;
- auto e2 = problem.model().globalFvGeometry().element(scvf.outsideScvIdx());
- auto e3 = problem.model().globalFvGeometry().element(scvfVector[otherScvfIdx]->outsideScvIdx());
+ const auto e2 = problem.model().globalFvGeometry().element(scvf.outsideScvIdx());
+ const auto e3 = problem.model().globalFvGeometry().element(scvfVector[otherScvfIdx]->outsideScvIdx());
auto e2Geometry = localView(problem.model().globalFvGeometry());
auto e3Geometry = localView(problem.model().globalFvGeometry());
@@ -166,10 +165,10 @@ private:
e3Scvfs[localScvfIdx]->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());
+ const 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 e4Scvfs = Implementation::getCommonAndNextScvFace(*e2Scvfs[localScvfIdx], e4Geometry, /*clockwise?*/localScvfIdx == 1);
outerScvSeeds.emplace_back(e4Scvfs[otherScvfIdx]->insideScvIdx(),
e4Scvfs[otherScvfIdx]->index());
diff --git a/dumux/discretization/cellcentered/mpfa/lmethod/interactionregions.hh b/dumux/discretization/cellcentered/mpfa/lmethod/interactionregions.hh
index fc4c039ea1d955f95c7aa93de31bbc978e0217dd..17380bc1a95bd01334e766af94ce86e5109bbbe2 100644
--- a/dumux/discretization/cellcentered/mpfa/lmethod/interactionregions.hh
+++ b/dumux/discretization/cellcentered/mpfa/lmethod/interactionregions.hh
@@ -58,22 +58,22 @@ public:
std::vector<Scalar> detX;
std::vector<Element> elements;
- std::array<GlobalIndexType, 2> globalScvfs;
+ std::array<GlobalIndexType, 2> globalScvfIndices;
// Constructor signature for dim == 2
- template<class ScvSeedType, class OuterScvSeedType>
+ template<class ScvSeed, class OuterScvSeed>
InteractionRegion(const Problem& problem,
const FVElementGeometry& fvGeometry,
- const ScvSeedType& scvSeed,
- const OuterScvSeedType& outerSeed1,
- const OuterScvSeedType& outerSeed2,
+ const ScvSeed& scvSeed,
+ const OuterScvSeed& outerSeed1,
+ const OuterScvSeed& 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();
+ globalScvfIndices[0] = scvSeed.globalScvfIndices()[contiFaceLocalIdx];
+ globalScvfIndices[1] = contiFaceLocalIdx == 0 ? outerSeed1.globalScvfIndex() : outerSeed2.globalScvfIndex();
elements.resize(3);
elements[0] = element1;
@@ -81,11 +81,11 @@ public:
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]);
+ const auto& scv1 = fvGeometry.scv(scvSeed.globalIndex());
+ const auto& scv2 = fvGeometry.scv(outerSeed1.globalIndex());
+ const auto& scv3 = fvGeometry.scv(outerSeed2.globalIndex());
+ const auto& scvf1 = fvGeometry.scvf(scvSeed.globalScvfIndices()[0]);
+ const auto& scvf2 = fvGeometry.scvf(scvSeed.globalScvfIndices()[1]);
// The necessary coordinates and normals
GlobalPosition v = scvf1.vertexCorner();
diff --git a/dumux/discretization/cellcentered/mpfa/lmethod/interactionvolume.hh b/dumux/discretization/cellcentered/mpfa/lmethod/interactionvolume.hh
index 51e7929705249eb21b8656195fdd125552499394..6ab3030efd62216e261933b8b7cb2e4d56fe7e6c 100644
--- a/dumux/discretization/cellcentered/mpfa/lmethod/interactionvolume.hh
+++ b/dumux/discretization/cellcentered/mpfa/lmethod/interactionvolume.hh
@@ -80,6 +80,7 @@ class CCMpfaInteractionVolumeImplementation<TypeTag, MpfaMethods::lMethod> : pub
using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);
+ using InteriorBoundaryData = typename GET_PROP_TYPE(TypeTag, InteriorBoundaryData);
static const int dim = GridView::dimension;
static const int dimWorld = GridView::dimensionworld;
@@ -93,6 +94,7 @@ class CCMpfaInteractionVolumeImplementation<TypeTag, MpfaMethods::lMethod> : pub
public:
using Matrix = typename Traits::Matrix;
+ using typename ParentType::GlobalLocalFaceDataPair;
using typename ParentType::LocalIndexType;
using typename ParentType::LocalIndexSet;
using typename ParentType::LocalFaceData;
@@ -108,7 +110,8 @@ public:
const Problem& problem,
const FVElementGeometry& fvGeometry,
const ElementVolumeVariables& elemVolVars)
- : problemPtr_(&problem),
+ : seedPtr_(&seed),
+ problemPtr_(&problem),
fvGeometryPtr_(&fvGeometry),
elemVolVarsPtr_(&elemVolVars),
regionUnique_(seed.isUnique()),
@@ -125,7 +128,9 @@ public:
{
const auto& ir = interactionRegions_[0];
T_ = assembleMatrix_(getTensor, ir);
- postSolve_(ir);
+ interactionRegionId_ = 0;
+ updateGlobalLocalFaceData();
+ systemSolved_ = true;
}
else
{
@@ -133,22 +138,25 @@ public:
M[0] = assembleMatrix_(getTensor, interactionRegions_[0]);
M[1] = assembleMatrix_(getTensor, interactionRegions_[1]);
- auto id = MpfaHelper::selectionCriterion(interactionRegions_[0], interactionRegions_[1], M[0], M[1]);
- postSolve_(interactionRegions_[id]);
- T_ = M[id];
+ interactionRegionId_ = MpfaHelper::selectionCriterion(interactionRegions_[0], interactionRegions_[1], M[0], M[1]);
+ updateGlobalLocalFaceData();
+ systemSolved_ = true;
+ T_ = M[interactionRegionId_];
}
}
//! 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
+ const 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.");
+ assert((scvf.index() == globalLocalScvfPairedData_[0].first->index()
+ || scvf.index() == globalLocalScvfPairedData_[1].first->index())
+ && "The provided scvf is not the flux face of the interaction volume.");
- if (globalScvfIndices_[0] == scvf.index())
- return LocalFaceData(contiFaceLocalIdx_, /*dummy*/0, false);
+ if (globalLocalScvfPairedData_[0].first->index() == scvf.index())
+ return globalLocalScvfPairedData_[0].second;
else
- return LocalFaceData(contiFaceLocalIdx_, /*dummy*/0, true);
+ return globalLocalScvfPairedData_[1].second;
}
//! returns the transmissibilities corresponding to the bound scvf
@@ -164,21 +172,18 @@ public:
const GlobalIndexSet& volVarsStencil() const
{
- assert(systemSolved_ && "volVarsStencil not set yet. You have to call solveLocalSystem() beforehand.");
- return volVarsStencil_;
+ assert(systemSolved_ && "stencil not set yet. You have to call solveLocalSystem() beforehand.");
+ return interactionRegion().scvIndices;
}
const PositionVector& volVarsPositions() const
{
- assert(systemSolved_ && "globalScvfs not set yet. You have to call solveLocalSystem() beforehand.");
- return volVarsPositions_;
+ assert(systemSolved_ && "stencil not set yet. You have to call solveLocalSystem() beforehand.");
+ return interactionRegion().scvCenters;
}
- const GlobalIndexSet& globalScvfs() const
- {
- assert(systemSolved_ && "globalScvfs not set yet. You have to call solveLocalSystem() beforehand.");
- return globalScvfIndices_;
- }
+ const std::vector<GlobalLocalFaceDataPair>& globalLocalScvfPairedData() const
+ { return globalLocalScvfPairedData_; }
//! Boundaries will be treated by a different mpfa method (e.g. o method). Thus, on
//! faces in l-method interaction volumes there will never be a Neumann flux contribution.
@@ -192,6 +197,13 @@ public:
const Matrix& matrix() const
{ return T_; }
+ const InteractionRegion& interactionRegion() const
+ { return interactionRegions_[interactionRegionId_]; }
+
+ // The l-method can't treat interior boundaries
+ const std::vector<InteriorBoundaryData> interiorBoundaryData() const
+ { return std::vector<InteriorBoundaryData>(); }
+
private:
//! Assembles and solves the local equation system
//! Specialization for dim = 2
@@ -210,9 +222,9 @@ private:
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);
+ const auto T1 = getTensor(problem_(), e1, elemVolVars_()[scv1], fvGeometry_(), scv1);
+ const auto T2 = getTensor(problem_(), e2, elemVolVars_()[scv2], fvGeometry_(), scv2);
+ const auto T3 = getTensor(problem_(), e3, elemVolVars_()[scv3], fvGeometry_(), scv3);
// required omega factors
Scalar w111 = calculateOmega_(ir.normal[0], ir.nu[0], ir.detX[0], T1);
@@ -248,7 +260,7 @@ private:
B[1][1] = 0.0;
B[1][2] = -w235 - w236;
- // T = CA^-1B + D
+ // T = C*A^-1*B + D
A.invert();
auto T = (A.leftmultiply(C)).rightmultiplyany(B);
T[0][0] += w111 + w112;
@@ -265,53 +277,57 @@ private:
{
const auto& globalFvGeometry = problem_().model().globalFvGeometry();
- if (regionUnique_)
+ interactionRegions_.reserve(2);
+ auto&& scvSeed1 = seed.scvSeed(0);
+ auto&& scvSeed2 = seed.scvSeed(1);
+ auto&& outerScvSeed1 = seed.outerScvSeed(0);
+ 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_.reserve(1);
- auto&& scvSeed1 = seed.scvSeed(0);
- auto&& scvSeed2 = seed.outerScvSeed(0);
- 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);
-
+ 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_.reserve(2);
- auto&& scvSeed1 = seed.scvSeed(0);
- auto&& scvSeed2 = seed.scvSeed(1);
- auto&& outerScvSeed1 = seed.outerScvSeed(0);
- 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);
- }
+ interactionRegions_.emplace_back(problem_(), fvGeometry, scvSeed1, outerScvSeed1, OuterScvSeedType(scvSeed2), e1, e3, e2);
+ interactionRegions_.emplace_back(problem_(), fvGeometry, scvSeed2, OuterScvSeedType(scvSeed1), outerScvSeed2, e2, e1, e4);
}
+
+ // as an initial guess we set the first interaction region as the one used
+ // this is updated after solution of the local system
+ interactionRegionId_ = 0;
+ updateGlobalLocalFaceData();
}
- void postSolve_(const InteractionRegion& ir)
+ //! set the global&local face data according to the current choice for the interaction region
+ void updateGlobalLocalFaceData()
{
- globalScvfIndices_.resize(2);
- globalScvfIndices_[0] = ir.globalScvfs[0];
- globalScvfIndices_[1] = ir.globalScvfs[1];
- volVarsStencil_ = ir.scvIndices;
- volVarsPositions_ = ir.scvCenters;
- contiFaceLocalIdx_ = ir.contiFaceLocalIdx;
- systemSolved_ = true;
+ const auto numPairs = globalLocalScvfPairedData_.size();
+
+ // if no data has been stored yet, initialize
+ if (numPairs == 0)
+ {
+ globalLocalScvfPairedData_.clear();
+ globalLocalScvfPairedData_.reserve(2);
+ globalLocalScvfPairedData_.emplace_back(&fvGeometry_().scvf(seed_().globalScvfIndices()[0]),
+ LocalFaceData(interactionRegion().contiFaceLocalIdx, /*dummy*/0, false));
+ globalLocalScvfPairedData_.emplace_back(&fvGeometry_().scvf(seed_().globalScvfIndices()[1]),
+ LocalFaceData(interactionRegion().contiFaceLocalIdx, /*dummy*/0, true));
+ }
+ // if order was swapped, change bools indicating which face is the "outside" version of the local face
+ else if (globalLocalScvfPairedData_[1].first->index() == interactionRegion().globalScvfIndices[0])
+ {
+ globalLocalScvfPairedData_[0].second.localScvfIndex = interactionRegion().contiFaceLocalIdx;
+ globalLocalScvfPairedData_[1].second.localScvfIndex = interactionRegion().contiFaceLocalIdx;
+ globalLocalScvfPairedData_[0].second.isOutside = true;
+ globalLocalScvfPairedData_[1].second.isOutside = false;
+ }
}
Scalar calculateOmega_(const GlobalPosition& normal,
@@ -347,7 +363,7 @@ private:
}
const Seed& seed_() const
- { return seedPtr_; }
+ { return *seedPtr_; }
const Problem& problem_() const
{ return *problemPtr_; }
@@ -366,12 +382,9 @@ private:
bool regionUnique_;
bool systemSolved_;
- LocalIndexType contiFaceLocalIdx_;
- GlobalIndexSet globalScvfIndices_;
- GlobalIndexSet volVarsStencil_;
- PositionVector volVarsPositions_;
-
+ LocalIndexType interactionRegionId_;
std::vector<InteractionRegion> interactionRegions_;
+ std::vector<GlobalLocalFaceDataPair> globalLocalScvfPairedData_;
Matrix T_;
};
diff --git a/dumux/discretization/cellcentered/mpfa/lmethod/subcontrolvolumeface.hh b/dumux/discretization/cellcentered/mpfa/lmethod/subcontrolvolumeface.hh
index 13e90f8edc57ffaaea75353a59fe54a100b7250d..9468f8c2fb2b12d043f76ee6cbbf20de45f2b95e 100644
--- a/dumux/discretization/cellcentered/mpfa/lmethod/subcontrolvolumeface.hh
+++ b/dumux/discretization/cellcentered/mpfa/lmethod/subcontrolvolumeface.hh
@@ -61,8 +61,8 @@ public:
Scalar q,
bool boundary)
: ParentType(helper,
- std::move(corners),
- std::move(unitOuterNormal),
+ std::forward<std::vector<GlobalPosition>>(corners),
+ std::forward<GlobalPosition>(unitOuterNormal),
vertexIndex,
localIndex,
scvfIndex,
diff --git a/dumux/discretization/cellcentered/mpfa/omethod/interactionvolume.hh b/dumux/discretization/cellcentered/mpfa/omethod/interactionvolume.hh
index 444fc6a823bf5866652f14127c7957da6f390c29..73991700b7ae86df988dcd21ecf8c5476b41f66f 100644
--- a/dumux/discretization/cellcentered/mpfa/omethod/interactionvolume.hh
+++ b/dumux/discretization/cellcentered/mpfa/omethod/interactionvolume.hh
@@ -128,6 +128,7 @@ class CCMpfaOInteractionVolume : public CCMpfaInteractionVolumeBase<TypeTag, Tra
using LocalScvfType = typename Traits::LocalScvfType;
public:
+ using typename ParentType::GlobalLocalFaceDataPair;
using typename ParentType::LocalIndexType;
using typename ParentType::LocalIndexSet;
using typename ParentType::LocalFaceData;
@@ -135,9 +136,6 @@ public:
using typename ParentType::PositionVector;
using typename ParentType::Seed;
- // structure to store global and local face data
- using GlobalLocalFaceDataPair = std::pair<const SubControlVolumeFace*, LocalFaceData>;
-
CCMpfaOInteractionVolume(const Seed& seed,
const Problem& problem,
const FVElementGeometry& fvGeometry,
diff --git a/dumux/implicit/cellcentered/assemblymap.hh b/dumux/implicit/cellcentered/assemblymap.hh
index 58c40166350bd8dbb09184e0f71891be73101beb..eaf1ee205d8ffba74960161ea98ad8f6cdfdc363 100644
--- a/dumux/implicit/cellcentered/assemblymap.hh
+++ b/dumux/implicit/cellcentered/assemblymap.hh
@@ -22,8 +22,8 @@
* that contribute to the derivative calculation. This is used for
* finite-volume schemes with symmetric sparsity pattern in the global matrix.
*/
-#ifndef DUMUX_CC_SYMMETRIC_ASSEMBLY_MAP_HH
-#define DUMUX_CC_SYMMETRIC_ASSEMBLY_MAP_HH
+#ifndef DUMUX_CC_ASSEMBLY_MAP_HH
+#define DUMUX_CC_ASSEMBLY_MAP_HH
#include <dune/istl/bcrsmatrix.hh>
@@ -32,8 +32,18 @@
namespace Dumux
{
+/*!
+ * \ingroup CellCentered
+ * \brief A simple version of the assembly map for cellcentered schemes.
+ * This implementation works for schemes in which for a given cell I only
+ * those cells J have to be prepared in whose stencil the cell I appears.
+ * This means that for the flux calculations in the cells J (in order to compute
+ * the derivatives with respect to cell I), we do not need data on any additional cells J
+ * to compute these fluxes. The same holds for scvfs in the cells J, i.e. we need only those
+ * scvfs in the cells J in which the cell I is in the stencil.
+ */
template<class TypeTag>
-class CCSymmetricAssemblyMap
+class CCSimpleAssemblyMap
{
using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
@@ -45,6 +55,9 @@ class CCSymmetricAssemblyMap
{
IndexType globalJ;
std::vector<IndexType> scvfsJ;
+ // A list of additional scvfs is needed for compatibility
+ // reasons with more complex assembly maps (see mpfa)
+ std::vector<IndexType> additionalScvfs;
};
using Map = std::vector<std::vector<DataJ>>;
diff --git a/dumux/implicit/cellcentered/mpfa/assemblymap.hh b/dumux/implicit/cellcentered/mpfa/assemblymap.hh
index cf9b374fba30f5b29447b9a43a31bcd151e5011d..76b6bca57a3928ff6ec4e86903acff330dae3132 100644
--- a/dumux/implicit/cellcentered/mpfa/assemblymap.hh
+++ b/dumux/implicit/cellcentered/mpfa/assemblymap.hh
@@ -25,6 +25,7 @@
#define DUMUX_CC_MPFA_ASSEMBLY_MAP_HH
#include <dumux/implicit/cellcentered/assemblymap.hh>
+#include <dumux/implicit/cellcentered/mpfa/generalassemblymap.hh>
namespace Dumux
{
@@ -36,9 +37,13 @@ class CCMpfaAssemblyMapImplementation;
template<class TypeTag>
using CCMpfaAssemblyMap = CCMpfaAssemblyMapImplementation<TypeTag, GET_PROP_VALUE(TypeTag, MpfaMethod)>;
-//! The default is simply the CCAssemblyMap
+//! The default is the general assembly map for mpfa schemes
template<class TypeTag, MpfaMethods method>
-class CCMpfaAssemblyMapImplementation : public CCAssemblyMap<TypeTag> {};
+class CCMpfaAssemblyMapImplementation : public CCMpfaGeneralAssemblyMap<TypeTag> {};
+
+//! The o-method can use the simple assembly map
+template<class TypeTag>
+class CCMpfaAssemblyMapImplementation<TypeTag, MpfaMethods::oMethod> : public CCSimpleAssemblyMap<TypeTag> {};
}
#endif
diff --git a/dumux/implicit/cellcentered/mpfa/generalassemblymap.hh b/dumux/implicit/cellcentered/mpfa/generalassemblymap.hh
new file mode 100644
index 0000000000000000000000000000000000000000..c69911e08de49ca1dc55e01f785ea6bc8f3421e7
--- /dev/null
+++ b/dumux/implicit/cellcentered/mpfa/generalassemblymap.hh
@@ -0,0 +1,201 @@
+// -*- 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 Stores the face indices corresponding to the neighbors of an element
+ * that contribute to the derivative calculation
+ */
+#ifndef DUMUX_CC_MPFA_GENERAL_ASSEMBLY_MAP_HH
+#define DUMUX_CC_MPFA_GENERAL_ASSEMBLY_MAP_HH
+
+#include <dumux/implicit/cellcentered/assemblymap.hh>
+
+namespace Dumux
+{
+
+/*!
+ * \ingroup CellCentered
+ * \brief General version of the assembly map for cellcentered schemes. To each
+ * cell I we store a list of cells J that are needed to compute the fluxes
+ * in these cells J that depend on cell I. Furthermore, we store for each cell J
+ * a list of scvfs in which cell I is in the stencil, as well as additional scvfs
+ * that are also required to set up the transmissibilities.
+ */
+template<class TypeTag>
+class CCMpfaGeneralAssemblyMap
+{
+ using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
+ using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+ using Element = typename GridView::template Codim<0>::Entity;
+ using MpfaHelper = typename GET_PROP_TYPE(TypeTag, MpfaHelper);
+ using FluxVariables = typename GET_PROP_TYPE(TypeTag, FluxVariables);
+ using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
+
+ using IndexType = typename GridView::IndexSet::IndexType;
+
+ // To each cell "globalI" there will be a list of "globalJ", in which globalI is part
+ // of the stencil. We save the scvfs over which fluxes depend on globalI and a list of
+ // additional scvfs which are needed temporarily to set up the transmissibilities of the scvfsJ
+ struct DataJ
+ {
+ IndexType globalJ;
+ std::vector<IndexType> scvfsJ;
+ std::vector<IndexType> additionalScvfs;
+ };
+
+ using Map = std::vector<std::vector<DataJ>>;
+
+public:
+
+ /*!
+ * \brief Initialize the AssemblyMap object.
+ *
+ * \param problem The problem which we want to simulate.
+ */
+ void init(const Problem& problem)
+ {
+ map_.resize(problem.gridView().size(0));
+ for (const auto& element : elements(problem.gridView()))
+ {
+ // We are looking for the elements I, for which this element J is in the flux stencil
+ auto globalJ = problem.elementMapper().index(element);
+
+ auto fvGeometry = localView(problem.model().globalFvGeometry());
+ fvGeometry.bindElement(element);
+
+ // obtain the data of J in elements I
+ std::vector<std::pair<IndexType, std::vector<DataJ>>> dataJForI;
+
+ // loop over sub control faces
+ for (auto&& scvf : scvfs(fvGeometry))
+ {
+ FluxVariables fluxVars;
+ const auto& stencil = fluxVars.computeStencil(problem, element, fvGeometry, scvf);
+
+ // insert our index in the neighbor stencils of the elements in the flux stencil
+ for (auto globalI : stencil)
+ {
+ if (globalI == globalJ)
+ continue;
+
+ auto it = std::find_if(dataJForI.begin(),
+ dataJForI.end(),
+ [globalI](const auto& pair) { return pair.first == globalI; });
+
+ if (it != dataJForI.end())
+ {
+ // get the data J which corresponds to the actual global J
+ // This will be the first entry, as we do so in the else statement (see below)
+ auto& globalJDataJ = it->second[0];
+
+ // insert actual scvf index in the list of scvfs which couple I and J
+ globalJDataJ.scvfsJ.push_back(scvf.index());
+
+ // Also, all scvfs connected to a vertex together with the actual scvf
+ // land in the list of additional scvfs. Of that list we will delete those
+ // that are already in the list of scvfsJ later...
+ const auto scvfVectorAtVertex = MpfaHelper::getScvFacesAtVertex(scvf.vertexIndex(), element, fvGeometry);
+ std::vector<IndexType> scvfIndicesAtVertex(scvfVectorAtVertex.size());
+ for (unsigned int i = 0; i < scvfVectorAtVertex.size(); ++i)
+ scvfIndicesAtVertex[i] = scvfVectorAtVertex[i]->index();
+ globalJDataJ.additionalScvfs.insert(globalJDataJ.additionalScvfs.end(),
+ scvfIndicesAtVertex.begin(),
+ scvfIndicesAtVertex.end());
+
+ // all the other dofs in the stencil have to appear as "globalJ" to globalI as well
+ for (auto globalJ2 : stencil)
+ {
+ if (globalJ2 == globalJ || globalJ2 == globalI)
+ continue;
+
+ auto it2 = std::find_if(it->second.begin(),
+ it->second.end(),
+ [globalJ2](const auto& dataJ) { return dataJ.globalJ == globalJ2; });
+
+ // if entry for globalJ2 does not exist yet, add globalJ2 to the J-data of globalI
+ // with an empty set of scvfs over which I and J are coupled (i.e. they aren't coupled)
+ if (it2 == it->second.end())
+ it->second.push_back(DataJ({globalJ2, std::vector<IndexType>()}));
+ }
+ }
+ else
+ {
+ // No DataJ for globalI exists yet. Make it and insert data on the actual
+ // global J as first entry in the vector of DataJs belonging to globalI
+ dataJForI.emplace_back(std::make_pair(globalI,
+ std::vector<DataJ>({DataJ({globalJ, std::vector<IndexType>({scvf.index()})})})));
+
+ // Also, all scvfs connected to a vertex together with the actual scvf
+ // land in the list of additional scvfs. Of that list we will delete those
+ // that are already in the list of scvfsJ later...
+ const auto scvfVectorAtVertex = MpfaHelper::getScvFacesAtVertex(scvf.vertexIndex(), element, fvGeometry);
+ std::vector<IndexType> scvfIndicesAtVertex(scvfVectorAtVertex.size());
+ for (unsigned int i = 0; i < scvfVectorAtVertex.size(); ++i)
+ scvfIndicesAtVertex[i] = scvfVectorAtVertex[i]->index();
+ dataJForI.back().second[0].additionalScvfs.insert(dataJForI.back().second[0].additionalScvfs.end(),
+ scvfIndicesAtVertex.begin(),
+ scvfIndicesAtVertex.end());
+
+ // all the other dofs in the stencil will be "globalJ" to globalI as well
+ for (auto globalJ2 : stencil)
+ if (globalJ2 != globalJ && globalJ2 != globalI)
+ dataJForI.back().second.push_back(DataJ({globalJ2, std::vector<IndexType>()}));
+ }
+ }
+ }
+
+ // Insert the data into the global map
+ for (auto&& pair : dataJForI)
+ {
+ // obtain the corresponding entry in the map
+ auto& dataJVector = map_[pair.first];
+ for (auto&& dataJ : pair.second)
+ {
+ // delete those additionalScvfs indices that are already in the list of scvfs
+ dataJ.additionalScvfs.erase(std::remove_if(dataJ.additionalScvfs.begin(),
+ dataJ.additionalScvfs.end(),
+ [&dataJ] (const auto& idx) { return MpfaHelper::contains(dataJ.scvfsJ, idx); }),
+ dataJ.additionalScvfs.end());
+
+ // if entry for j exists in the map already add scvf and additional scvf indices, create otherwise
+ auto it = std::find_if(dataJVector.begin(),
+ dataJVector.end(),
+ [&dataJ](const auto& dataJofMap) { return dataJofMap.globalJ == dataJ.globalJ; });
+
+ if (it != dataJVector.end())
+ {
+ it->scvfsJ.insert(it->scvfsJ.end(), dataJ.scvfsJ.begin(), dataJ.scvfsJ.end());
+ it->additionalScvfs.insert(it->additionalScvfs.end(), dataJ.additionalScvfs.begin(), dataJ.additionalScvfs.end());
+ }
+ else
+ dataJVector.emplace_back(std::move(dataJ));
+ }
+ }
+ }
+ }
+
+ const std::vector<DataJ>& operator[] (const IndexType globalI) const
+ { return map_[globalI]; }
+
+private:
+ Map map_;
+};
+}
+
+#endif
diff --git a/dumux/implicit/cellcentered/mpfa/propertydefaults.hh b/dumux/implicit/cellcentered/mpfa/propertydefaults.hh
index 483f831e8367bbe97157c1a43e74944574674f9d..0f0a0dbc26d5b7728dc8aa0605630c2f8dfca79d 100644
--- a/dumux/implicit/cellcentered/mpfa/propertydefaults.hh
+++ b/dumux/implicit/cellcentered/mpfa/propertydefaults.hh
@@ -43,6 +43,7 @@
#include <dumux/discretization/cellcentered/mpfa/globalinteractionvolumeseeds.hh>
#include <dumux/discretization/cellcentered/mpfa/interiorboundarydata.hh>
#include <dumux/implicit/cellcentered/mpfa/localresidual.hh>
+#include <dumux/implicit/cellcentered/mpfa/assemblymap.hh>
#include <dumux/implicit/cellcentered/properties.hh>
namespace Dumux {
@@ -66,6 +67,9 @@ SET_PROP(CCMpfaModel, MpfaMethod)
//! The mpfa helper class
SET_TYPE_PROP(CCMpfaModel, MpfaHelper, CCMpfaHelper<TypeTag>);
+//! The assembly map for mpfa schemes
+SET_TYPE_PROP(CCMpfaModel, AssemblyMap, CCMpfaAssemblyMap<TypeTag>);
+
//! The interaction volume class
SET_TYPE_PROP(CCMpfaModel, InteractionVolume, CCMpfaInteractionVolume<TypeTag>);
diff --git a/dumux/implicit/cellcentered/propertydefaults.hh b/dumux/implicit/cellcentered/propertydefaults.hh
index 702259c31e8b32641980917887ad55f61c822ec8..99fbb8925c9fadc81dfe84ac82ad9ae3929d5e99 100644
--- a/dumux/implicit/cellcentered/propertydefaults.hh
+++ b/dumux/implicit/cellcentered/propertydefaults.hh
@@ -76,8 +76,8 @@ SET_TYPE_PROP(CCModel, GlobalVolumeVariables, CCGlobalVolumeVariables<TypeTag, G
//! Set the BaseLocalResidual to CCLocalResidual
SET_TYPE_PROP(CCModel, BaseLocalResidual, CCLocalResidual<TypeTag>);
-//! Set the AssemblyMap to the CCAssemblyMap (default: symmetric occupation pattern)
-SET_TYPE_PROP(CCModel, AssemblyMap, Dumux::CCSymmetricAssemblyMap<TypeTag>);
+//! Set the AssemblyMap to the CCAssemblyMap per default
+SET_TYPE_PROP(CCModel, AssemblyMap, Dumux::CCSimpleAssemblyMap<TypeTag>);
// By default, we disable interior boundaries
SET_BOOL_PROP(CCModel, EnableInteriorBoundaries, false);