From f03d50e7259bbcd97c76ef1f24a87a6d37002b28 Mon Sep 17 00:00:00 2001 From: DennisGlaeser Date: Mon, 28 Nov 2016 16:25:47 +0100 Subject: [PATCH 1/5] [mpfa][elemvolvars] only make boundary vol vars on dirichlet boundaries --- .../mpfa/elementvolumevariables.hh | 51 +++++++++++-------- 1 file changed, 30 insertions(+), 21 deletions(-) diff --git a/dumux/discretization/cellcentered/mpfa/elementvolumevariables.hh b/dumux/discretization/cellcentered/mpfa/elementvolumevariables.hh index 8dc18e813d..03e8db6eb0 100644 --- a/dumux/discretization/cellcentered/mpfa/elementvolumevariables.hh +++ b/dumux/discretization/cellcentered/mpfa/elementvolumevariables.hh @@ -24,6 +24,7 @@ #define DUMUX_DISCRETIZATION_CCMPFA_ELEMENT_VOLUMEVARIABLES_HH #include +#include "facetypes.hh" namespace Dumux { @@ -171,16 +172,20 @@ public: 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); + // only proceed if we are on a dirichlet boundary + if (MpfaHelper::getMpfaFaceType(problem, element, scvf) == MpfaFaceTypes::dirichlet) + { + // 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); + } } } @@ -202,19 +207,23 @@ public: // 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); + // only proceed if we are on a dirichlet boundary + if (MpfaHelper::getMpfaFaceType(problem, ivElement, ivScvf) == MpfaFaceTypes::dirichlet) + { + // boundary volume variables + VolumeVariables dirichletVolVars; + auto&& ivScv = fvGeometry.scv(insideScvIdx); + 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); + } } } -- GitLab From e3d70e357f6c531d3ce74d39dd24b1733303b299 Mon Sep 17 00:00:00 2001 From: DennisGlaeser Date: Mon, 28 Nov 2016 16:26:37 +0100 Subject: [PATCH 2/5] [mpfa][globalfvgeom] use indexInInside() to detect branching points we don't want to compare positions for the detection of branching points. Note that the new implementation fails for hanging nodes on a 2d grid embedded in a 3d domain. --- .../cellcentered/mpfa/globalfvgeometrybase.hh | 20 +++++++++---------- 1 file changed, 10 insertions(+), 10 deletions(-) diff --git a/dumux/discretization/cellcentered/mpfa/globalfvgeometrybase.hh b/dumux/discretization/cellcentered/mpfa/globalfvgeometrybase.hh index ef5420f94e..5d8c27ebb8 100644 --- a/dumux/discretization/cellcentered/mpfa/globalfvgeometrybase.hh +++ b/dumux/discretization/cellcentered/mpfa/globalfvgeometrybase.hh @@ -164,15 +164,16 @@ public: // for dim < dimWorld we'll have multiple intersections at one point // we store here the centers of those we have handled already - std::vector finishedCenters; + std::vector finishedIndices; // construct the sub control volume faces 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 + // Note that for hanging nodes on the lower dimension network, it doesn't work if (dim < dimWorld) - if(MpfaHelper::contains(finishedCenters, is.geometry().center())) + if(MpfaHelper::contains(finishedIndices, is.indexInInside())) continue; // get some of the intersection's bools @@ -189,12 +190,11 @@ public: 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) + if (nextIs.indexInInside() == insideIndex) { auto nIdx = problem.elementMapper().index(nextIs.outside()); if (nIdx != nIndices[0]) @@ -203,7 +203,7 @@ public: } // store this center as a handled one - finishedCenters.push_back(isCenter); + finishedIndices.push_back(insideIndex); } } else if (boundary) @@ -444,7 +444,7 @@ public: // for dim < dimWorld we'll have multiple intersections at one point // we store here the centers of those we have handled already - std::vector finishedCenters; + std::vector finishedIndices; unsigned int localFaceIdx = 0; // construct the sub control volume faces @@ -452,8 +452,9 @@ public: { // if we are dealing with a lower dimensional network // only make a new scvf if we haven't handled it yet + // Note that for hanging nodes on the lower dimension network, it doesn't work if (dim < dimWorld) - if(MpfaHelper::contains(finishedCenters, is.geometry().center())) + if(MpfaHelper::contains(finishedIndices, is.indexInInside())) continue; // get some of the intersection bools @@ -470,12 +471,11 @@ public: 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) + if (nextIs.indexInInside() == insideIndex) { auto nIdx = problem.elementMapper().index(nextIs.outside()); if (nIdx != nIndices[0]) @@ -484,7 +484,7 @@ public: } // store this center as a handled one - finishedCenters.push_back(isCenter); + finishedIndices.push_back(insideIndex); } } else if (boundary) -- GitLab From 6518c179d0ff90a0ea78530436508f4ff7a9447d Mon Sep 17 00:00:00 2001 From: DennisGlaeser Date: Mon, 28 Nov 2016 16:33:41 +0100 Subject: [PATCH 3/5] [mpfa][localseeds] merge functions to add outside data --- .../cellcentered/mpfa/omethod/helper.hh | 16 ++++++---------- .../mpfa/omethod/localsubcontrolentityseeds.hh | 11 ++++------- 2 files changed, 10 insertions(+), 17 deletions(-) diff --git a/dumux/discretization/cellcentered/mpfa/omethod/helper.hh b/dumux/discretization/cellcentered/mpfa/omethod/helper.hh index d12415c2da..91fda3c12a 100644 --- a/dumux/discretization/cellcentered/mpfa/omethod/helper.hh +++ b/dumux/discretization/cellcentered/mpfa/omethod/helper.hh @@ -429,6 +429,7 @@ private: outsideLocalScvfIdx++; if (scvfSeed.insideGlobalScvIndex() == actualScvf.insideScvIdx() && + scvfSeed.outsideGlobalScvIndices().size() == actualScvf.outsideScvIndices().size() && std::equal(scvfSeed.outsideGlobalScvIndices().begin(), scvfSeed.outsideGlobalScvIndices().end(), actualScvf.outsideScvIndices().begin())) insideScvfExists = true; else @@ -454,20 +455,17 @@ private: Implementation::getMpfaFaceType(problem, element, actualScvf)); // pass the actual outside indices to the new scvf seed - scvfSeeds.back().addOutsideScvfIndex(outsideScvfIdx); - scvfSeeds.back().addOutsideLocalScvIndex(static_cast(scvSeeds.size())); + scvfSeeds.back().addOutsideData(outsideScvfIdx, static_cast(scvSeeds.size())); } else if (insideScvfExists && !outsideScvfExists) { // pass info on outside to the inside scvf seed - scvfSeeds[insideLocalScvfIdx].addOutsideScvfIndex(outsideScvfIdx); - scvfSeeds[insideLocalScvfIdx].addOutsideLocalScvIndex(static_cast(scvSeeds.size())); + scvfSeeds[insideLocalScvfIdx].addOutsideData(outsideScvfIdx, static_cast(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(scvSeeds.size())); + scvfSeeds[outsideLocalScvfIdx].addOutsideData(actualScvf.index(), static_cast(scvSeeds.size())); } // make outside scv by recursion @@ -485,8 +483,7 @@ private: Implementation::getMpfaFaceType(problem, element, actualScvf)); // pass the actual outside indices to the new scvf seed - scvfSeeds.back().addOutsideScvfIndex(outsideScvfIdx); - scvfSeeds.back().addOutsideLocalScvIndex(outsideLocalScvIdx); + scvfSeeds.back().addOutsideData(outsideScvfIdx, outsideLocalScvIdx); } else if (outsideScvExists && !insideScvfExists && outsideScvfExists) { @@ -494,8 +491,7 @@ private: actualScvSeed.setLocalScvfIndex(coordDir, outsideLocalScvfIdx); // pass info on outside to the inside found local scvf seed - scvfSeeds[outsideLocalScvfIdx].addOutsideScvfIndex(outsideScvfIdx); - scvfSeeds[outsideLocalScvfIdx].addOutsideLocalScvIndex(actualLocalScvIdx); + scvfSeeds[outsideLocalScvfIdx].addOutsideData(actualScvf.index(), actualLocalScvIdx); } } } diff --git a/dumux/discretization/cellcentered/mpfa/omethod/localsubcontrolentityseeds.hh b/dumux/discretization/cellcentered/mpfa/omethod/localsubcontrolentityseeds.hh index 49fbc1ad2b..28289b0a7e 100644 --- a/dumux/discretization/cellcentered/mpfa/omethod/localsubcontrolentityseeds.hh +++ b/dumux/discretization/cellcentered/mpfa/omethod/localsubcontrolentityseeds.hh @@ -153,14 +153,11 @@ public: bool boundary() const { return boundary_; } - void addOutsideScvfIndex(const GlobalIndexType index) + void addOutsideData(const GlobalIndexType outsideGlobalScvfIndex, + const LocalIndexType outsideLocalScvIndex) { - outsideGlobalScvfIndices_.push_back(index); - } - - void addOutsideLocalScvIndex(const LocalIndexType index) - { - outsideLocalScvIndices_.push_back(index); + outsideLocalScvIndices_.push_back(outsideLocalScvIndex); + outsideGlobalScvfIndices_.push_back(outsideGlobalScvfIndex); } void makeOutsideDataUnique() -- GitLab From f0703fc5e9cb3a3748f2883cce29f3b60700fc07 Mon Sep 17 00:00:00 2001 From: DennisGlaeser Date: Mon, 28 Nov 2016 16:35:28 +0100 Subject: [PATCH 4/5] [mpfa][omethod] use correct outside normal vector in case of lower dimensional grids --- .../mpfa/omethod/interactionvolume.hh | 43 ++++++++++++++----- .../omethod/localsubcontrolentityseeds.hh | 22 +++++++--- 2 files changed, 49 insertions(+), 16 deletions(-) diff --git a/dumux/discretization/cellcentered/mpfa/omethod/interactionvolume.hh b/dumux/discretization/cellcentered/mpfa/omethod/interactionvolume.hh index 245ba09e53..7091b5213c 100644 --- a/dumux/discretization/cellcentered/mpfa/omethod/interactionvolume.hh +++ b/dumux/discretization/cellcentered/mpfa/omethod/interactionvolume.hh @@ -450,7 +450,7 @@ private: auto tensor = getTensor(element, elemVolVars_()[posGlobalScv], posGlobalScv); // the omega factors of the "positive" sub volume - auto posWijk = calculateOmegas_(posLocalScv, localScvf, tensor); + auto posWijk = calculateOmegas_(posLocalScv, localScvf.unitOuterNormal(), localScvf.area(), tensor); posWijk *= problem_().boxExtrusionFactor(element, posGlobalScv); // Check the local directions of the positive sub volume @@ -493,15 +493,30 @@ private: if (faceType == MpfaFaceTypes::interior) { // loop over all the outside neighbors of this face and add entries + unsigned int indexInOutsideData = 0; for (auto negLocalScvIdx : localScvf.outsideLocalScvIndices()) { const auto& negLocalScv = localScv_(negLocalScvIdx); const auto& negGlobalScv = fvGeometry_().scv(negLocalScv.globalIndex()); - auto negElement = localElement_(negLocalScvIdx);; + 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); + DimVector negWijk; + + // if dim < dimWorld, use outside normal vector + if (dim < dimWorld) + { + // outside scvf + auto&& outsideScvf = fvGeometry_().scvf(localScvf.outsideGlobalScvfIndices()[indexInOutsideData]); + auto negNormal = outsideScvf.unitOuterNormal(); + negNormal *= -1.0; + negWijk = calculateOmegas_(negLocalScv, negNormal, localScvf.area(), negTensor); + } + else + negWijk = calculateOmegas_(negLocalScv, localScvf.unitOuterNormal(), localScvf.area(), negTensor); + + // scale by extrusion factpr negWijk *= problem_().boxExtrusionFactor(negElement, negGlobalScv); // Check local directions of negative sub volume @@ -527,8 +542,12 @@ private: B[idxInFluxFaces][negLocalScvIdx] -= negWijk[localDir]; } - // store the omegas + // store the omegas (negative, because of normal vector sign switch) + negWijk *= -1.0; wijk_[rowIdx].emplace_back(std::move(negWijk)); + + // increment counter in outside data + indexInOutsideData++; } } // go to the next face @@ -563,7 +582,7 @@ private: auto tensor = getTensor(element, elemVolVars_()[posGlobalScv], posGlobalScv); // the omega factors of the "positive" sub volume - auto posWijk = calculateOmegas_(posLocalScv, localScvf, tensor); + auto posWijk = calculateOmegas_(posLocalScv, localScvf.unitOuterNormal(), localScvf.area(), tensor); posWijk *= problem_().boxExtrusionFactor(element, posGlobalScv); for (int localDir = 0; localDir < dim; localDir++) @@ -584,7 +603,8 @@ private: } DimVector calculateOmegas_(const LocalScvType& localScv, - const LocalScvfType& localScvf, + const GlobalPosition normal, + const Scalar area, const Tensor& T) const { DimVector wijk; @@ -592,9 +612,9 @@ private: for (int dir = 0; dir < dim; ++dir) { T.mv(localScv.innerNormal(dir), tmp); - wijk[dir] = tmp*localScvf.unitOuterNormal(); + wijk[dir] = tmp*normal; } - wijk *= localScvf.area(); + wijk *= area; wijk /= localScv.detX(); wijk *= -1.0; @@ -602,16 +622,17 @@ private: } DimVector calculateOmegas_(const LocalScvType& localScv, - const LocalScvfType& localScvf, + const GlobalPosition normal, + const Scalar area, const Scalar t) const { DimVector wijk; - GlobalPosition tmp(localScvf.unitOuterNormal()); + GlobalPosition tmp(normal); tmp *= t; for (int dir = 0; dir < dim; ++dir) wijk[dir] = tmp*localScv.innerNormal(dir); - wijk *= localScvf.area(); + wijk *= area; wijk /= localScv.detX(); wijk *= -1.0; diff --git a/dumux/discretization/cellcentered/mpfa/omethod/localsubcontrolentityseeds.hh b/dumux/discretization/cellcentered/mpfa/omethod/localsubcontrolentityseeds.hh index 28289b0a7e..9567cab120 100644 --- a/dumux/discretization/cellcentered/mpfa/omethod/localsubcontrolentityseeds.hh +++ b/dumux/discretization/cellcentered/mpfa/omethod/localsubcontrolentityseeds.hh @@ -160,13 +160,25 @@ public: outsideGlobalScvfIndices_.push_back(outsideGlobalScvfIndex); } + // for grids with dim < dimWorld, some outside indices might be doubled + // we want to make the outside indices unique, but, the i-th outside global scvf face + // should correspond to the j-th outside local scv.Therefore we apply the same operations on both containers 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()); + for (auto scvIt = outsideLocalScvIndices_.begin(); scvIt != outsideLocalScvIndices_.end(); ++scvIt) + { + auto scvfIt = outsideGlobalScvfIndices_.begin(); + for (auto scvIt2 = scvIt+1; scvIt2 != outsideLocalScvIndices_.end(); ++scvIt2) + { + if (*scvIt2 == *scvIt) + { + outsideLocalScvIndices_.erase(scvIt2); + outsideGlobalScvfIndices_.erase(scvfIt); + break; + } + ++scvfIt; + } + } } private: -- GitLab From 18f65640edae0a367670ee8526e564e152e1bb15 Mon Sep 17 00:00:00 2001 From: DennisGlaeser Date: Mon, 28 Nov 2016 17:35:04 +0100 Subject: [PATCH 5/5] [mpfa][omethod][helper] forward to 2d helper when scvf does not touch a branching point for dim < dimWorld (i.e. 2 < 3), the algorithm for the setup of the interactionvolume seed is more complex than the one for dim = dimWorld = 2. This commit introduces a forwarding to the more simple algorithm if the face does not touch a branching point. --- .../cellcentered/mpfa/globalfvgeometrybase.hh | 30 +++++++++++++++++-- .../cellcentered/mpfa/helper.hh | 5 ++++ .../cellcentered/mpfa/omethod/helper.hh | 8 +++++ 3 files changed, 41 insertions(+), 2 deletions(-) diff --git a/dumux/discretization/cellcentered/mpfa/globalfvgeometrybase.hh b/dumux/discretization/cellcentered/mpfa/globalfvgeometrybase.hh index 5d8c27ebb8..d92b5e3150 100644 --- a/dumux/discretization/cellcentered/mpfa/globalfvgeometrybase.hh +++ b/dumux/discretization/cellcentered/mpfa/globalfvgeometrybase.hh @@ -111,10 +111,14 @@ public: const BoundaryInteractionVolumeSeed& boundaryInteractionVolumeSeed(const SubControlVolumeFace& scvf) const { return globalInteractionVolumeSeeds_.boundarySeed(scvf); } - //! Returns whether or not a scvf touches the boundary (has to be called before getting an interaction volume) + //! Returns whether or not an scvf touches the boundary (has to be called before getting an interaction volume) bool scvfTouchesBoundary(const SubControlVolumeFace& scvf) const { return boundaryVertices_[scvf.vertexIndex()]; } + //! Returns whether or not an scvf touches a branching point (for dim < dimWorld) + bool scvfTouchesBranchingPoint(const SubControlVolumeFace& scvf) const + { return branchingVertices_[scvf.vertexIndex()]; } + //! update all fvElementGeometries (do this again after grid adaption) void update(const Problem& problem) { @@ -136,6 +140,10 @@ public: boundaryVertices_.resize(gridView_.size(dim), false); elementMap_.resize(numScvs); + // keep track of branching points + if (dim < dimWorld) + branchingVertices_.resize(gridView_.size(dim), false); + // find vertices on processor boundaries auto isGhostVertex = MpfaHelper::findGhostVertices(problem, gridView_); @@ -241,6 +249,10 @@ public: if (boundary) boundaryVertices_[vIdxGlobal] = true; + // is vertex on a branching point? + if (dim < dimWorld && nIndices.size() > 1) + branchingVertices_[vIdxGlobal] = true; + // make the scv face scvfIndexSet.push_back(scvfIdx); scvfs_.emplace_back(helper, @@ -317,6 +329,7 @@ private: // vectors that store the global data std::vector> scvfIndicesOfScv_; std::vector boundaryVertices_; + std::vector branchingVertices_; IndexType numBoundaryScvf_; // the global interaction volume seeds @@ -392,7 +405,7 @@ public: const BoundaryInteractionVolumeSeed& boundaryInteractionVolumeSeed(const SubControlVolumeFace& scvf) const { return globalInteractionVolumeSeeds_.boundarySeed(scvf); } - //! Returns whether or not a scvf touches the boundary (has to be called before getting an interaction volume) + //! Returns whether or not an scvf touches the boundary (has to be called before getting an interaction volume) bool scvfTouchesBoundary(const SubControlVolumeFace& scvf) const { return boundaryVertices_[scvf.vertexIndex()]; } @@ -400,6 +413,10 @@ public: bool isGhostVertex(const IndexType vIdxGlobal) const { return ghostVertices_[vIdxGlobal]; } + //! Returns whether or not an scvf touches a branching point (for dim < dimWorld) + bool scvfTouchesBranchingPoint(const SubControlVolumeFace& scvf) const + { return branchingVertices_[scvf.vertexIndex()]; } + //! update all fvElementGeometries (do this again after grid adaption) void update(const Problem& problem) { @@ -420,6 +437,10 @@ public: neighborVolVarIndices_.resize(numScvs_); boundaryVertices_.resize(gridView_.size(dim), false); + // keep track of branching points + if (dim < dimWorld) + branchingVertices_.resize(gridView_.size(dim), false); + // find vertices on processor boundaries ghostVertices_ = MpfaHelper::findGhostVertices(problem, gridView_); @@ -508,6 +529,10 @@ public: if (boundary) boundaryVertices_[vIdxGlobal] = true; + // is vertex on a branching point? + if (dim < dimWorld && nIndices.size() > 1) + branchingVertices_[vIdxGlobal] = true; + // store information on the scv face scvfsIndexSet.push_back(numScvf_++); neighborVolVarIndexSet.push_back(nIndices); @@ -563,6 +588,7 @@ private: std::vector< std::vector< std::vector > > neighborVolVarIndices_; std::vector boundaryVertices_; std::vector ghostVertices_; + std::vector branchingVertices_; // the global interaction volume seeds GlobalInteractionVolumeSeeds globalInteractionVolumeSeeds_; diff --git a/dumux/discretization/cellcentered/mpfa/helper.hh b/dumux/discretization/cellcentered/mpfa/helper.hh index a2f0430659..2b131c515d 100644 --- a/dumux/discretization/cellcentered/mpfa/helper.hh +++ b/dumux/discretization/cellcentered/mpfa/helper.hh @@ -248,6 +248,11 @@ public: static std::size_t getNumLocalScvfs(const Dune::GeometryType gt) { return MpfaDimensionHelper::getNumLocalScvfs(gt); } + + static ScvfVector getCommonAndNextScvFace(const SubControlVolumeFace& outsideScvf, + const FVElementGeometry& fvGeometry, + const bool clockWise) + { return MpfaDimensionHelper::getCommonAndNextScvFace(outsideScvf, fvGeometry, clockWise); } }; // Specialization for dim == 3 (dimWorld has to be 3) diff --git a/dumux/discretization/cellcentered/mpfa/omethod/helper.hh b/dumux/discretization/cellcentered/mpfa/omethod/helper.hh index 91fda3c12a..ae505e43ba 100644 --- a/dumux/discretization/cellcentered/mpfa/omethod/helper.hh +++ b/dumux/discretization/cellcentered/mpfa/omethod/helper.hh @@ -290,6 +290,10 @@ public: const FVElementGeometry& fvGeometry, const SubControlVolumeFace& scvf) { + // if the scvf does not touch a branching point, use simplified algorithm to create interaction volume seed + if (!problem.model().globalFvGeometry().scvfTouchesBranchingPoint(scvf)) + return MpfaMethodHelper::makeInnerInteractionVolumeSeed(problem, element, fvGeometry, scvf); + std::vector scvSeeds; std::vector scvfSeeds; @@ -321,6 +325,10 @@ public: const FVElementGeometry& fvGeometry, const SubControlVolumeFace& scvf) { + // if the scvf does not touch a branching point, use simplified algorithm to create interaction volume seed + if (!problem.model().globalFvGeometry().scvfTouchesBranchingPoint(scvf)) + return MpfaMethodHelper::makeBoundaryInteractionVolumeSeed(problem, element, fvGeometry, scvf); + std::vector scvSeeds; std::vector scvfSeeds; -- GitLab