diff --git a/dumux/discretization/cellcentered/mpfa/fvelementgeometry.hh b/dumux/discretization/cellcentered/mpfa/fvelementgeometry.hh
index 557b70dad7da825c71ab7501130d3f6c60dec7b1..1d0a9f19c829b3a1243924fbaec1046e23561144 100644
--- a/dumux/discretization/cellcentered/mpfa/fvelementgeometry.hh
+++ b/dumux/discretization/cellcentered/mpfa/fvelementgeometry.hh
@@ -222,22 +222,7 @@ public:
//! Note that e.g. the normals might be different in the case of surface grids
const SubControlVolumeFace& flipScvf(GridIndexType scvfIdx, unsigned int outsideScvIdx = 0) const
{
- auto it = std::find(scvfIndices_.begin(), scvfIndices_.end(), scvfIdx);
- if (it != scvfIndices_.end())
- {
- const auto localScvfIdx = std::distance(scvfIndices_.begin(), it);
- return neighborScvfs_[flipScvfIndices_[localScvfIdx][outsideScvIdx]];
- }
- else
- {
- const auto neighborScvfIdxLocal = findLocalIndex(scvfIdx, neighborScvfIndices_);
- const auto& scvf = neighborScvfs_[neighborScvfIdxLocal];
-
- if (scvf.outsideScvIdx(outsideScvIdx) == scvIndices_[0])
- return scvfs_[neighborFlipScvfIndices_[neighborScvfIdxLocal][outsideScvIdx]];
- else
- return neighborScvfs_[neighborFlipScvfIndices_[neighborScvfIdxLocal][outsideScvIdx]];
- }
+ return scvf( fvGridGeometry().flipScvfIdx(scvfIdx, outsideScvIdx) );
}
//! iterator range for sub control volumes. Iterates over
@@ -299,10 +284,6 @@ public:
dataJ.scvfsJ,
dataJ.additionalScvfs);
- // set up the scvf flip indices for network grids
- if (dim < dimWorld)
- makeFlipIndexSet();
-
// //! TODO Check if user added additional DOF dependencies, i.e. the residual of DOF globalI depends
// //! on additional DOFs not included in the discretization schemes' occupation pattern
// const auto& additionalDofDependencies = problem.getAdditionalDofDependencies(globalI);
@@ -524,98 +505,6 @@ private:
}
}
- //! find the "flip" indices for all scvfs
- void makeFlipIndexSet()
- {
- const auto numInsideScvfs = scvfIndices_.size();
- const auto numNeighborScvfs = neighborScvfIndices_.size();
-
- // first, handle the interior scvfs (flip scvf will be in outside scvfs)
- flipScvfIndices_.resize(numInsideScvfs);
- for (unsigned int insideFace = 0; insideFace < numInsideScvfs; ++insideFace)
- {
- const auto& scvf = scvfs_[insideFace];
- if (scvf.boundary())
- continue;
-
- const auto numOutsideScvs = scvf.numOutsideScvs();
- const auto vIdxGlobal = scvf.vertexIndex();
- const auto insideScvIdx = scvf.insideScvIdx();
-
- flipScvfIndices_[insideFace].resize(numOutsideScvs);
- for (unsigned int nIdx = 0; nIdx < numOutsideScvs; ++nIdx)
- {
- const auto outsideScvIdx = scvf.outsideScvIdx(nIdx);
- for (unsigned int outsideFace = 0; outsideFace < numNeighborScvfs; ++outsideFace)
- {
- const auto& outsideScvf = neighborScvfs_[outsideFace];
- // check if outside face is the flip face
- if (outsideScvf.insideScvIdx() == outsideScvIdx &&
- outsideScvf.vertexIndex() == vIdxGlobal &&
- MpfaHelper::vectorContainsValue(outsideScvf.outsideScvIndices(), insideScvIdx))
- {
- flipScvfIndices_[insideFace][nIdx] = outsideFace;
- // there is always only one flip face in an outside element
- break;
- }
- }
- }
- }
-
- // Now we look for the flip indices of the outside faces
- neighborFlipScvfIndices_.resize(numNeighborScvfs);
- for (unsigned int outsideFace = 0; outsideFace < numNeighborScvfs; ++outsideFace)
- {
- const auto& scvf = neighborScvfs_[outsideFace];
- if (scvf.boundary())
- continue;
-
- const auto numOutsideScvs = scvf.numOutsideScvs();
- const auto vIdxGlobal = scvf.vertexIndex();
- const auto insideScvIdx = scvf.insideScvIdx();
-
- neighborFlipScvfIndices_[outsideFace].resize(numOutsideScvs);
- for (unsigned int nIdx = 0; nIdx < numOutsideScvs; ++nIdx)
- {
- const auto neighborScvIdx = scvf.outsideScvIdx(nIdx);
-
- // if the neighbor scv idx is the index of the bound element,
- // then the flip scvf will be within the inside scvfs
- if (neighborScvIdx == scvIndices_[0])
- {
- for (unsigned int insideFace = 0; insideFace < numInsideScvfs; ++insideFace)
- {
- const auto& insideScvf = scvfs_[insideFace];
- // check if the face is the flip face
- if (insideScvf.vertexIndex() == vIdxGlobal &&
- MpfaHelper::vectorContainsValue(insideScvf.outsideScvIndices(), insideScvIdx))
- {
- neighborFlipScvfIndices_[outsideFace][nIdx] = insideFace;
- // there is always only one flip face in an outside element
- break;
- }
- }
- }
- else
- {
- for (unsigned int outsideFace2 = 0; outsideFace2 < numNeighborScvfs; ++outsideFace2)
- {
- const auto& outsideScvf = neighborScvfs_[outsideFace2];
- // check if outside face is the flip face
- if (outsideScvf.insideScvIdx() == neighborScvIdx &&
- outsideScvf.vertexIndex() == vIdxGlobal &&
- MpfaHelper::vectorContainsValue(outsideScvf.outsideScvIndices(), insideScvIdx))
- {
- neighborFlipScvfIndices_[outsideFace][nIdx] = outsideFace2;
- // there is always only one flip face in an outside element
- break;
- }
- }
- }
- }
- }
- }
-
//! map a global index to the local storage index
const unsigned int findLocalIndex(const GridIndexType idx,
const std::vector<GridIndexType>& indices) const
@@ -635,9 +524,6 @@ private:
neighborScvfIndices_.clear();
neighborScvs_.clear();
neighborScvfs_.clear();
-
- flipScvfIndices_.clear();
- neighborFlipScvfIndices_.clear();
}
const FVGridGeometry* fvGridGeometryPtr_;
@@ -652,10 +538,6 @@ private:
std::vector<GridIndexType> neighborScvfIndices_;
std::vector<SubControlVolume> neighborScvs_;
std::vector<SubControlVolumeFace> neighborScvfs_;
-
- // flip scvf index sets
- std::vector< std::vector<GridIndexType> > flipScvfIndices_;
- std::vector< std::vector<GridIndexType> > neighborFlipScvfIndices_;
};
} // end namespace
diff --git a/dumux/discretization/cellcentered/mpfa/fvgridgeometry.hh b/dumux/discretization/cellcentered/mpfa/fvgridgeometry.hh
index aee1736c8868e503e56aa8c6ae4188d41118c71c..e56f3b1dd7b7c797409f57fb93ecd1b58e289c73 100644
--- a/dumux/discretization/cellcentered/mpfa/fvgridgeometry.hh
+++ b/dumux/discretization/cellcentered/mpfa/fvgridgeometry.hh
@@ -329,7 +329,6 @@ public:
break;
}
}
-
}
}
}
@@ -524,6 +523,13 @@ public:
// instantiate the dual grid index set (to be used for construction of interaction volumes)
typename GridIVIndexSets::DualGridIndexSet dualIdSet(this->gridView());
+ // keep track of boundary scvfs and scvf vertex indices in order to set up flip scvf index set
+ const auto maxNumScvfs = numScvs_*LocalView::maxNumElementScvfs;
+ std::vector<bool> scvfIsOnBoundary;
+ std::vector<GridIndexType> scvfVertexIndex;
+ scvfIsOnBoundary.reserve(maxNumScvfs);
+ scvfVertexIndex.reserve(maxNumScvfs);
+
// Build the SCVs and SCV faces
numScvf_ = 0;
numBoundaryScvf_ = 0;
@@ -614,6 +620,8 @@ public:
// store information on the scv face
scvfsIndexSet.push_back(numScvf_++);
+ scvfIsOnBoundary.push_back(boundary);
+ scvfVertexIndex.push_back(vIdxGlobal);
neighborVolVarIndexSet.emplace_back(std::move(outsideScvIndices));
}
@@ -627,6 +635,46 @@ public:
neighborVolVarIndices_[eIdx] = neighborVolVarIndexSet;
}
+ // Make the flip index set for network and surface grids
+ if (dim < dimWorld)
+ {
+ flipScvfIndices_.resize(numScvf_);
+ for (std::size_t scvIdx = 0; scvIdx < numScvs_; ++scvIdx)
+ {
+ const auto& scvfIndices = scvfIndicesOfScv_[scvIdx];
+ for (unsigned int i = 0; i < scvfIndices.size(); ++i)
+ {
+ // boundary scvf have no flip scvfs
+ if (scvfIsOnBoundary[ scvfIndices[i] ])
+ continue;
+
+ const auto scvfIdx = scvfIndices[i];
+ const auto vIdxGlobal = scvfVertexIndex[scvfIdx];
+ const auto numOutsideScvs = neighborVolVarIndices_[scvIdx][i].size();
+
+ flipScvfIndices_[scvfIdx].resize(numOutsideScvs);
+ for (unsigned int j = 0; j < numOutsideScvs; ++j)
+ {
+ const auto outsideScvIdx = neighborVolVarIndices_[scvIdx][i][j];
+ const auto& outsideScvfIndices = scvfIndicesOfScv_[outsideScvIdx];
+ for (unsigned int k = 0; k < outsideScvfIndices.size(); ++k)
+ {
+ const auto outsideScvfIndex = outsideScvfIndices[k];
+ const auto outsideScvfVertexIndex = scvfVertexIndex[outsideScvfIndex];
+ const auto& outsideScvfNeighborIndices = neighborVolVarIndices_[outsideScvIdx][k];
+ if (outsideScvfVertexIndex == vIdxGlobal &&
+ MpfaHelper::vectorContainsValue(outsideScvfNeighborIndices, scvIdx))
+ {
+ flipScvfIndices_[scvfIdx][j] = outsideScvfIndex;
+ // there is always only one flip face in an outside element
+ break;
+ }
+ }
+ }
+ }
+ }
+ }
+
// building the geometries has finished
std::cout << "Initializing of the grid finite volume geometry took " << timer.elapsed() << " seconds." << std::endl;
@@ -649,6 +697,11 @@ public:
const std::vector<ScvfOutsideGridIndexStorage>& neighborVolVarIndices(GridIndexType scvIdx) const
{ return neighborVolVarIndices_[scvIdx]; }
+ //! Get the index scvf on the same face but from the other side
+ //! Note that e.g. the normals might be different in the case of surface grids
+ const GridIndexType flipScvfIdx(GridIndexType scvfIdx, unsigned int outsideScvfIdx = 0) const
+ { return flipScvfIndices_[scvfIdx][outsideScvfIdx]; }
+
//! Returns the connectivity map of which dofs
//! have derivatives with respect to a given dof.
const ConnectivityMap& connectivityMap() const { return connectivityMap_; }
@@ -669,6 +722,9 @@ private:
GridIndexType numScvf_;
GridIndexType numBoundaryScvf_;
+ // needed for embedded surface and network grids (dim < dimWorld)
+ std::vector<ScvfOutsideGridIndexStorage> flipScvfIndices_;
+
// The grid interaction volume index set
GridIVIndexSets ivIndexSets_;