From f18d3401bbd6a746cca6e9723bf29451ce6ca45c Mon Sep 17 00:00:00 2001
From: Kilian Weishaupt <kilian.weishaupt@iws.uni-stuttgart.de>
Date: Wed, 27 Oct 2021 13:15:24 +0200
Subject: [PATCH] [staggered][gridGeometry][scvf] Restructure insertion of
scvfs, add new boundary flag
* write scvfs directly into vector using global index and fixed local indices
* add flag in scvf for processor boundary
---
.../staggered/fvelementgeometry.hh | 115 ++++---
.../facecentered/staggered/fvgridgeometry.hh | 312 +++++++-----------
.../staggered/subcontrolvolumeface.hh | 24 +-
3 files changed, 210 insertions(+), 241 deletions(-)
diff --git a/dumux/discretization/facecentered/staggered/fvelementgeometry.hh b/dumux/discretization/facecentered/staggered/fvelementgeometry.hh
index 598e234290..58d008df61 100644
--- a/dumux/discretization/facecentered/staggered/fvelementgeometry.hh
+++ b/dumux/discretization/facecentered/staggered/fvelementgeometry.hh
@@ -24,16 +24,16 @@
#ifndef DUMUX_DISCRETIZATION_FACECENTERED_STAGGERED_FV_ELEMENT_GEOMETRY_HH
#define DUMUX_DISCRETIZATION_FACECENTERED_STAGGERED_FV_ELEMENT_GEOMETRY_HH
-#include <bitset>
-#include <type_traits>
#include <utility>
+#include <dune/common/rangeutilities.hh>
#include <dune/common/reservedvector.hh>
#include <dumux/common/indextraits.hh>
#include <dune/common/iteratorrange.hh>
#include <dumux/discretization/scvandscvfiterators.hh>
#include <dumux/discretization/facecentered/staggered/normalaxis.hh>
+#include <dumux/discretization/facecentered/staggered/consistentlyorientedgrid.hh>
namespace Dumux {
@@ -144,9 +144,7 @@ public:
//! number of sub control volumes in this fv element geometry
std::size_t numScvf() const
- {
- return scvfIndices_().size();
- }
+ { return scvfIndices_().size(); }
//! Returns whether one of the geometry's scvfs lies on a boundary
bool hasBoundaryScvf() const
@@ -244,7 +242,7 @@ public:
using GridGeometry = GG;
//! the maximum number of scvs per element
- static constexpr std::size_t maxNumElementScvs = 2*GridView::dimension;
+ static constexpr std::size_t maxNumElementScvs = numScvsPerElement;
FaceCenteredStaggeredFVElementGeometry(const GridGeometry& gridGeometry)
: gridGeometry_(&gridGeometry)
@@ -307,8 +305,8 @@ public:
const SubControlVolumeFace& lateralOrthogonalScvf(const SubControlVolumeFace& scvf) const
{
assert(scvf.isLateral());
- const auto otherGlobalIdx = gridGeometry().lateralOrthogonalScvf(scvf);
- return scvfs_[findLocalIndex_(otherGlobalIdx, scvfIndices_())];
+ const auto otherLocalIdx = GridGeometry::GeometryHelper::lateralOrthogonalScvfLocalIndex(scvf.localIndex());
+ return scvfs_[otherLocalIdx];
}
//! Return the frontal sub control volume face on a the boundary for a given sub control volume
@@ -360,12 +358,18 @@ public:
*/
FaceCenteredStaggeredFVElementGeometry bindElement(const Element& element) &&
{
- this->bindElement_(element);
+ typename GG::LocalIntersectionMapper localIsMapper;
+ localIsMapper.update(gridGeometry().gridView(), element);
+ this->bindElement_(element, localIsMapper);
return std::move(*this);
}
void bindElement(const Element& element) &
- { this->bindElement_(element); }
+ {
+ typename GG::LocalIntersectionMapper localIsMapper;
+ localIsMapper.update(gridGeometry().gridView(), element);
+ this->bindElement_(element, localIsMapper);
+ }
GridIndexType elementIndex() const
{ return eIdx_; }
@@ -386,15 +390,18 @@ private:
//! called by the local jacobian to prepare element assembly
void bind_(const Element& element)
{
- bindElement(element);
+ typename GG::LocalIntersectionMapper localIsMapper;
+ localIsMapper.update(gridGeometry().gridView(), element);
+
+ bindElement_(element, localIsMapper);
neighborScvIndices_.clear();
neighborScvs_.clear();
- LocalIndexType localScvIdx = 0;
for (const auto& intersection : intersections(gridGeometry().gridView(), element))
{
if (intersection.neighbor())
{
+ const auto localScvIdx = localIsMapper.realToRefIdx(intersection.indexInInside());
const auto localOppositeScvIdx = geometryHelper_.localOppositeIdx(localScvIdx);
const auto& neighborElement = intersection.outside();
const auto neighborElementIdx = gridGeometry().elementMapper().index(neighborElement);
@@ -404,12 +411,16 @@ private:
std::array<GridIndexType, numScvsPerElement> globalScvIndicesOfNeighborElement;
std::iota(globalScvIndicesOfNeighborElement.begin(), globalScvIndicesOfNeighborElement.end(), neighborElementIdx*numScvsPerElement);
- LocalIndexType localNeighborScvIdx = 0;
+ typename GG::LocalIntersectionMapper localNeighborIsMapper;
+ localNeighborIsMapper.update(gridGeometry().gridView(), neighborElement);
+
for (const auto& neighborIntersection : intersections(gridGeometry().gridView(), neighborElement))
{
+ const auto localNeighborScvIdx = localNeighborIsMapper.realToRefIdx(neighborIntersection.indexInInside());
if (localNeighborScvIdx != localScvIdx && localNeighborScvIdx != localOppositeScvIdx)
{
- const auto dofIndex = gridGeometry().intersectionMapper().globalIntersectionIndex(neighborElement, localNeighborScvIdx);
+
+ const auto dofIndex = gridGeometry().intersectionMapper().globalIntersectionIndex(neighborElement, neighborIntersection.indexInInside());
neighborScvs_.push_back(SubControlVolume(
neighborElementGeometry,
neighborIntersection.geometry(),
@@ -423,31 +434,28 @@ private:
neighborScvIndices_.push_back(globalScvIndicesOfNeighborElement[localNeighborScvIdx]);
}
-
- ++localNeighborScvIdx;
}
}
-
- ++localScvIdx;
}
}
//! Bind only element-local
- void bindElement_(const Element& element)
+ void bindElement_(const Element& element, const typename GG::LocalIntersectionMapper& localIsMapper)
{
elementPtr_ = &element;
eIdx_ = gridGeometry().elementMapper().index(element);
std::iota(scvIndicesOfElement_.begin(), scvIndicesOfElement_.end(), eIdx_*numScvsPerElement);
scvfs_.clear();
+ scvfs_.resize(minNumScvfsPerElement);
- LocalIndexType localScvIdx = 0;
- LocalIndexType localScvfIdx = 0;
const auto& elementGeometry = element.geometry();
for (const auto& intersection : intersections(gridGeometry().gridView(), element))
{
+ const auto localScvIdx = localIsMapper.realToRefIdx(intersection.indexInInside());
+ auto localScvfIdx = localScvIdx*(1 + numLateralScvfsPerScv);
const auto& intersectionGeometry = intersection.geometry();
- const auto dofIndex = gridGeometry().intersectionMapper().globalIntersectionIndex(element, localScvIdx);
+ const auto dofIndex = gridGeometry().intersectionMapper().globalIntersectionIndex(element, intersection.indexInInside());
scvs_[localScvIdx] = SubControlVolume(
elementGeometry,
@@ -462,7 +470,7 @@ private:
// the frontal sub control volume face at the element center
const auto localOppositeScvIdx = geometryHelper_.localOppositeIdx(localScvIdx);
- scvfs_.push_back(SubControlVolumeFace(
+ scvfs_[localScvfIdx] = SubControlVolumeFace(
elementGeometry,
intersectionGeometry,
std::array{scvIndicesOfElement_[localScvIdx], scvIndicesOfElement_[localOppositeScvIdx]},
@@ -470,21 +478,19 @@ private:
scvfIndices_()[localScvfIdx],
intersection.centerUnitOuterNormal(),
SubControlVolumeFace::FaceType::frontal,
- false
- ));
+ SubControlVolumeFace::BoundaryType::interior
+ );
++localScvfIdx;
// the lateral sub control volume faces
- const auto lateralFaceIndices = geometryHelper_.localLaterFaceIndices(localScvIdx);
- for (const auto lateralFacetIndex : lateralFaceIndices)
+ for (const auto lateralFacetIndex : Dune::transformedRangeView(geometryHelper_.localLaterFaceIndices(localScvIdx),
+ [&](auto&& idx) { return localIsMapper.refToRealIdx(idx) ;})
+ )
{
const auto& lateralIntersection = geometryHelper_.intersection(lateralFacetIndex, element);
const auto& lateralFacet = geometryHelper_.facet(lateralFacetIndex, element);
const auto& lateralFacetGeometry = lateralFacet.geometry();
- if (lateralIntersection.neighbor())
- geometryHelper_.update(element, lateralIntersection.outside());
-
// helper lambda to get the lateral scvf's global inside and outside scv indices
const auto globalScvIndicesForLateralFace = [&]
{
@@ -493,12 +499,10 @@ private:
if (lateralIntersection.neighbor())
{
const auto parallelElemIdx = gridGeometry().elementMapper().index(lateralIntersection.outside());
- const auto localOutsideScvIdx = geometryHelper_.localFaceIndexInOtherElement(localScvIdx);
-
// todo: could be done easier?
std::array<GridIndexType, numScvsPerElement> globalScvIndicesOfNeighborElement;
std::iota(globalScvIndicesOfNeighborElement.begin(), globalScvIndicesOfNeighborElement.end(), parallelElemIdx*numScvsPerElement);
- return globalScvIndicesOfNeighborElement[localOutsideScvIdx];
+ return globalScvIndicesOfNeighborElement[localScvIdx];
}
else
return gridGeometry().outsideVolVarIndex(scvfIndices_()[localScvfIdx]);
@@ -507,7 +511,17 @@ private:
return std::array{scvIndicesOfElement_[localScvIdx], globalOutsideScvIdx};
}();
- scvfs_.push_back(SubControlVolumeFace(
+ const auto boundaryType = [&]
+ {
+ if (onProcessorBoundary_(lateralIntersection))
+ return SubControlVolumeFace::BoundaryType::processorBoundary;
+ else if (onDomainBoundary_(lateralIntersection))
+ return SubControlVolumeFace::BoundaryType::physicalBoundary;
+ else
+ return SubControlVolumeFace::BoundaryType::interior;
+ }();
+
+ scvfs_[localScvfIdx] = SubControlVolumeFace(
elementGeometry,
intersectionGeometry,
lateralFacetGeometry,
@@ -516,23 +530,20 @@ private:
scvfIndices_()[localScvfIdx],
lateralIntersection.centerUnitOuterNormal(),
SubControlVolumeFace::FaceType::lateral,
- onDomainBoundary_(lateralIntersection)
- ));
+ boundaryType
+ );
++localScvfIdx;
}
-
-
-
- ++localScvIdx;
}
// do a second loop over all intersections to add frontal boundary faces
- localScvIdx = 0;
+ auto localScvfIdx = minNumScvfsPerElement;
for (const auto& intersection : intersections(gridGeometry().gridView(), element))
{
// the frontal sub control volume face at a domain boundary (coincides with element face)
if (onDomainBoundary_(intersection))
{
+ const auto localScvIdx = localIsMapper.realToRefIdx(intersection.indexInInside());
// the frontal sub control volume face at the boundary
scvfs_.push_back(SubControlVolumeFace(
element.geometry(),
@@ -542,12 +553,23 @@ private:
scvfIndices_()[localScvfIdx],
intersection.centerUnitOuterNormal(),
SubControlVolumeFace::FaceType::frontal,
- true
+ SubControlVolumeFace::BoundaryType::physicalBoundary
));
- ++localScvfIdx;
+ ++localScvfIdx;
}
+ }
- ++localScvIdx;
+ if constexpr (!ConsistentlyOrientedGrid<typename GridView::Grid>{})
+ {
+ static const bool makeConsistentlyOriented = getParam<bool>("Grid.MakeConsistentlyOriented", true);
+ if (makeConsistentlyOriented && scvfs_.size() > minNumScvfsPerElement)
+ {
+ // make sure frontal boundary scvfs are sorted correctly
+ std::sort(scvfs_.begin() + minNumScvfsPerElement, scvfs_.end(),
+ [](const auto& scvfLeft, const auto& scvfRight)
+ { return scvfLeft.insideScvIdx() < scvfRight.insideScvIdx(); }
+ );
+ }
}
}
@@ -568,6 +590,11 @@ private:
return !intersection.neighbor() && intersection.boundary();
}
+ bool onProcessorBoundary_(const typename GridView::Intersection& intersection) const
+ {
+ return !intersection.neighbor() && !intersection.boundary();
+ }
+
Dune::ReservedVector<SubControlVolumeFace, maxNumScvfsPerElement> scvfs_;
Dune::ReservedVector<SubControlVolume, numLateralScvfsPerElement> neighborScvs_;
diff --git a/dumux/discretization/facecentered/staggered/fvgridgeometry.hh b/dumux/discretization/facecentered/staggered/fvgridgeometry.hh
index 6a40b4e055..7158921f2f 100644
--- a/dumux/discretization/facecentered/staggered/fvgridgeometry.hh
+++ b/dumux/discretization/facecentered/staggered/fvgridgeometry.hh
@@ -26,6 +26,8 @@
#include <memory>
+#include <dune/common/rangeutilities.hh>
+
#include <dumux/common/defaultmappertraits.hh>
#include <dumux/common/indextraits.hh>
#include <dumux/common/intersectionmapper.hh>
@@ -42,6 +44,7 @@
#include <dumux/discretization/facecentered/staggered/geometryhelper.hh>
#include <dumux/discretization/facecentered/staggered/connectivitymap.hh>
#include <dumux/discretization/facecentered/staggered/normalaxis.hh>
+#include <dumux/discretization/facecentered/staggered/localintersectionindexmapper.hh>
namespace Dumux {
@@ -57,6 +60,7 @@ struct FaceCenteredStaggeredDefaultGridGeometryTraits : public DefaultMapperTrai
using SubControlVolume = FaceCenteredStaggeredSubControlVolume<GridView>;
using SubControlVolumeFace = FaceCenteredStaggeredSubControlVolumeFace<GridView>;
using IntersectionMapper = ConformingGridIntersectionMapper<GridView>;
+ using LocalIntersectionMapper = FaceCenteredStaggeredLocalIntersectionIndexMapper<GridView>;
using GeometryHelper = FaceCenteredStaggeredGeometryHelper<GridView, typename GridView::Grid>;
template<class GridGeometry>
@@ -139,6 +143,8 @@ public:
using GridView = GV;
//! export the geometry helper type
using GeometryHelper = typename Traits::GeometryHelper;
+ //! export the local intersection mapper
+ using LocalIntersectionMapper = typename Traits::LocalIntersectionMapper;
//! export static information
using StaticInformation = typename Traits::StaticInfo;
//! export the type of extrusion
@@ -203,9 +209,9 @@ public:
//! all scvs of the element-local fvGeometry.
auto scvs(const LocalView& fvGeometry) const
{
- const auto begin = scvs_.cbegin() + numScvsPerElement*fvGeometry.elementIndex();
+ auto begin = scvs_.cbegin() + numScvsPerElement*fvGeometry.elementIndex();
const auto end = begin + numScvsPerElement;
- return Dune::IteratorRange<decltype(begin)>(begin, end);
+ return Dune::IteratorRange<std::decay_t<decltype(begin)>>(begin, end);
}
//! Get a sub control volume face with a global scvf index
@@ -243,10 +249,6 @@ public:
const std::unordered_map<GridIndexType, GridIndexType>& periodicVertexMap() const
{ return periodicFaceMap_; }
- //! get the global index of the orthogonal face sharing a common entity
- GridIndexType lateralOrthogonalScvf(const SubControlVolumeFace& scvf) const
- { return lateralOrthogonalScvf_[scvf.index()]; }
-
private:
void update_()
@@ -255,15 +257,12 @@ private:
scvs_.clear();
scvfs_.clear();
scvfIndicesOfElement_.clear();
- lateralOrthogonalScvf_.clear();
intersectionMapper_.update(this->gridView());
// determine size of containers
const auto numElements = this->gridView().size(0);
scvfIndicesOfElement_.resize(numElements);
hasBoundaryScvf_.resize(numElements, false);
- // on frontal + maybe boundary per face and 2*(dim-1) laterals per frontal
- lateralOrthogonalScvf_.resize(numElements*(2*dim*(2 + 2*(dim-1))));
outSideBoundaryVolVarIdx_ = 0;
numBoundaryScv_ = 0;
@@ -271,103 +270,65 @@ private:
GeometryHelper geometryHelper(this->gridView());
- // get the global scv indices first
- GridIndexType scvfIdx = 0;
+ // get the global scvf indices first
+ GridIndexType numScvfs = 0;
for (const auto& element : elements(this->gridView()))
{
assert(numScvsPerElement == element.subEntities(1));
- std::vector<GridIndexType> scvfsIndexSet;
- scvfsIndexSet.reserve(1/*frontal in element*/ + 1 /*frontal on boundary*/ + numLateralScvfsPerScv);
-
- // keep track of frontal boundary scvfs
- std::size_t numFrontalBoundaryScvfs = 0;
-
- // a temporary map to store pairs of common entities
- std::unordered_map<GridIndexType, Dune::ReservedVector<GridIndexType, 2>> commonEntityIdxToScvfsMap;
-
- SmallLocalIndexType localScvIdx = 0;
for (const auto& intersection : intersections(this->gridView(), element))
{
- // the frontal sub control volume face at the element center
- scvfsIndexSet.push_back(scvfIdx++);
-
- if constexpr(dim > 1)
- {
- // the lateral sub control volume faces
- for (const auto lateralFacetIndex : geometryHelper.localLaterFaceIndices(localScvIdx))
- {
- const auto& lateralIntersection = geometryHelper.intersection(lateralFacetIndex, element);
- if (onProcessorBoundary_(lateralIntersection))
- continue;
+ // frontal scvf in element center
+ ++numScvfs;
- scvfsIndexSet.push_back(scvfIdx);
-
- const auto commonEntityIdx = geometryHelper.globalCommonEntityIndex(
- element, localScvIdx, lateralFacetIndex
- );
- commonEntityIdxToScvfsMap[commonEntityIdx].push_back(scvfIdx);
- ++scvfIdx;
- }
- }
+ // lateral scvfs
+ numScvfs += numLateralScvfsPerScv;
// handle physical domain boundary
if (onDomainBoundary_(intersection))
{
++numBoundaryScv_; // frontal face
numBoundaryScv_ += numLateralScvfsPerScv; // boundary scvs for lateral faces
- ++numFrontalBoundaryScvfs;
- }
- ++localScvIdx;
- }
-
- // add global indices of frontal boundary scvfs last
- for (std::size_t i = 0; i < numFrontalBoundaryScvfs; ++i)
- scvfsIndexSet.push_back(scvfIdx++);
-
- // Save the scv indices belonging to this element to build up fv element geometries fast
- const auto eIdx = this->elementMapper().index(element);
- scvfIndicesOfElement_[eIdx] = std::move(scvfsIndexSet);
-
- // create the bi-directional map
- for (const auto& [_, scvfs] : commonEntityIdxToScvfsMap)
- {
- // TODO: this maybe be less than 2 if there is a processor boundary
- // are we sure that the lateral orthogonal scvf is then never needed?
- if (scvfs.size() == 2)
- {
- lateralOrthogonalScvf_[scvfs[0]] = scvfs[1];
- lateralOrthogonalScvf_[scvfs[1]] = scvfs[0];
+ // frontal scvf at boundary
+ ++numScvfs;
}
}
}
- // reserve memory
- const auto numInteriorScvs = numElements*numScvsPerElement;
- scvs_.resize(numInteriorScvs);
- scvfs_.reserve((numScvsPerElement/*one interior frontal scvf per scv*/
- +numLateralScvfsPerElement)*numElements + numBoundaryScv_);
+ // allocate memory
+ const auto numScvs = numElements*numScvsPerElement;
+ scvs_.resize(numScvs);
+ scvfs_.reserve(numScvfs);
// Build the scvs and scv faces
+ std::size_t globalScvfIdx = 0;
for (const auto& element : elements(this->gridView()))
{
const auto eIdx = this->elementMapper().index(element);
- const auto& globalScvfIndices = scvfIndicesOfElement_[eIdx];
+ auto& globalScvfIndices = scvfIndicesOfElement_[eIdx];
+ globalScvfIndices.resize(minNumScvfsPerElement);
+ globalScvfIndices.reserve(maxNumScvfsPerElement);
auto getGlobalScvIdx = [&](const auto elementIdx, const auto localScvIdx)
{ return numScvsPerElement*elementIdx + localScvIdx; };
- SmallLocalIndexType localScvfIdx = 0;
+ LocalIntersectionMapper localIsMapper;
+ localIsMapper.update(this->gridView(), element);
+
for (const auto& intersection : intersections(this->gridView(), element))
{
- const auto localScvIdx = intersection.indexInInside();
+ const auto& intersectionUnitOuterNormal = intersection.centerUnitOuterNormal();
+ const auto localScvIdx = localIsMapper.realToRefIdx(intersection.indexInInside());
+ auto localScvfIdx = localScvIdx*(1 + numLateralScvfsPerScv);
+
const auto globalScvIdx = getGlobalScvIdx(eIdx, localScvIdx);
- const auto dofIndex = intersectionMapper().globalIntersectionIndex(element, localScvIdx);
+ const auto dofIndex = intersectionMapper().globalIntersectionIndex(element, intersection.indexInInside());
const auto localOppositeScvIdx = geometryHelper.localOppositeIdx(localScvIdx);
const auto& intersectionGeometry = intersection.geometry();
const auto& elementGeometry = element.geometry();
+ assert(localIsMapper.refToRealIdx(localScvIdx) == intersection.indexInInside());
// handle periodic boundaries
if (onPeriodicBoundary_(intersection))
@@ -384,7 +345,7 @@ private:
if (periodicFaceFound)
continue;
- if (Dune::FloatCmp::eq(intersection.centerUnitOuterNormal()*otherIntersection.centerUnitOuterNormal(), -1.0, 1e-7))
+ if (Dune::FloatCmp::eq(intersectionUnitOuterNormal*otherIntersection.centerUnitOuterNormal(), -1.0, 1e-7))
{
const auto periodicDofIdx = intersectionMapper().globalIntersectionIndex(otherElement, otherIntersectionLocalIdx);
periodicFaceMap_[dofIndex] = periodicDofIdx;
@@ -396,39 +357,37 @@ private:
}
// the sub control volume
- scvs_[globalScvIdx] = SubControlVolume(elementGeometry,
- intersectionGeometry,
- globalScvIdx,
- localScvIdx,
- dofIndex,
- Dumux::normalAxis(intersection.centerUnitOuterNormal()),
- this->elementMapper().index(element),
- onDomainBoundary_(intersection));
+ scvs_[globalScvIdx] = SubControlVolume(
+ elementGeometry,
+ intersectionGeometry,
+ globalScvIdx,
+ localScvIdx,
+ dofIndex,
+ Dumux::normalAxis(intersectionUnitOuterNormal),
+ this->elementMapper().index(element),
+ onDomainBoundary_(intersection)
+ );
// the frontal sub control volume face at the element center
scvfs_.emplace_back(elementGeometry,
- intersectionGeometry,
- std::array{globalScvIdx, getGlobalScvIdx(eIdx, localOppositeScvIdx)},
- localScvfIdx,
- globalScvfIndices[localScvfIdx],
- intersection.centerUnitOuterNormal(),
- SubControlVolumeFace::FaceType::frontal,
- false);
+ intersectionGeometry,
+ std::array{globalScvIdx, getGlobalScvIdx(eIdx, localOppositeScvIdx)},
+ localScvfIdx,
+ globalScvfIdx,
+ intersectionUnitOuterNormal,
+ SubControlVolumeFace::FaceType::frontal,
+ SubControlVolumeFace::BoundaryType::interior
+ );
+
+ globalScvfIndices[localScvfIdx] = globalScvfIdx++;
++localScvfIdx;
// the lateral sub control volume faces
- const auto lateralFaceIndices = geometryHelper.localLaterFaceIndices(localScvIdx);
- for (const auto lateralFacetIndex : lateralFaceIndices)
+ for (const auto lateralFacetIndex : Dune::transformedRangeView(geometryHelper.localLaterFaceIndices(localScvIdx),
+ [&](auto&& idx) { return localIsMapper.refToRealIdx(idx) ;})
+ )
{
const auto& lateralIntersection = geometryHelper.intersection(lateralFacetIndex, element);
- if (onProcessorBoundary_(lateralIntersection))
- continue;
-
- const auto& lateralFacet = geometryHelper.facet(lateralFacetIndex, element);
- const auto& lateralFacetGeometry = lateralFacet.geometry();
-
- if (lateralIntersection.neighbor()) // TODO: periodic?
- geometryHelper.update(element, lateralIntersection.outside());
// helper lambda to get the lateral scvf's global inside and outside scv indices
const auto globalScvIndicesForLateralFace = [&]
@@ -438,25 +397,40 @@ private:
if (lateralIntersection.neighbor())
{
const auto parallelElemIdx = this->elementMapper().index(lateralIntersection.outside());
- const auto localOutsideScvIdx = geometryHelper.localFaceIndexInOtherElement(localScvIdx);
- return getGlobalScvIdx(parallelElemIdx, localOutsideScvIdx);
+ return getGlobalScvIdx(parallelElemIdx, localScvIdx);
}
+ else if (onDomainBoundary_(lateralIntersection))
+ return numScvs + outSideBoundaryVolVarIdx_++;
else
- return numInteriorScvs + outSideBoundaryVolVarIdx_++;
+ return globalScvIdx; // fallback for parallel, won't be used anyway
}();
return std::array{globalScvIdx, globalOutsideScvIdx};
}();
- scvfs_.emplace_back(elementGeometry,
- intersectionGeometry,
- lateralFacetGeometry,
- globalScvIndicesForLateralFace, // TODO higher order
- localScvfIdx,
- globalScvfIndices[localScvfIdx],
- lateralIntersection.centerUnitOuterNormal(),
- SubControlVolumeFace::FaceType::lateral,
- onDomainBoundary_(lateralIntersection));
+ const auto boundaryType = [&]
+ {
+ if (onProcessorBoundary_(lateralIntersection))
+ return SubControlVolumeFace::BoundaryType::processorBoundary;
+ else if (onDomainBoundary_(lateralIntersection))
+ return SubControlVolumeFace::BoundaryType::physicalBoundary;
+ else
+ return SubControlVolumeFace::BoundaryType::interior;
+ }();
+
+ scvfs_.emplace_back(
+ elementGeometry,
+ intersectionGeometry,
+ geometryHelper.facet(lateralFacetIndex, element).geometry(),
+ globalScvIndicesForLateralFace, // TODO higher order
+ localScvfIdx,
+ globalScvfIdx,
+ lateralIntersection.centerUnitOuterNormal(),
+ SubControlVolumeFace::FaceType::lateral,
+ boundaryType
+ );
+
+ globalScvfIndices[localScvfIdx] = globalScvfIdx++;
++localScvfIdx;
if (onDomainBoundary_(lateralIntersection))
@@ -469,24 +443,30 @@ private:
} // end first loop over intersections
// do a second loop over all intersections to add frontal boundary faces
+ int localScvfIdx = minNumScvfsPerElement;
for (const auto& intersection : intersections(this->gridView(), element))
{
// the frontal sub control volume face at a domain boundary (coincides with element face)
if (onDomainBoundary_(intersection))
{
- const auto localScvIdx = intersection.indexInInside();
+ const auto localScvIdx = localIsMapper.realToRefIdx(intersection.indexInInside());
const auto globalScvIdx = getGlobalScvIdx(eIdx, localScvIdx);
++numBoundaryScvf_;
// the frontal sub control volume face at the boundary
- scvfs_.emplace_back(element.geometry(),
- intersection.geometry(),
- std::array{globalScvIdx, globalScvIdx}, // TODO outside boundary, periodic, parallel?
- localScvfIdx,
- globalScvfIndices[localScvfIdx],
- intersection.centerUnitOuterNormal(),
- SubControlVolumeFace::FaceType::frontal,
- true);
+ scvfs_.emplace_back(
+ element.geometry(),
+ intersection.geometry(),
+ std::array{globalScvIdx, globalScvIdx}, // TODO outside boundary, periodic, parallel?
+ localScvfIdx,
+ globalScvfIdx,
+ intersection.centerUnitOuterNormal(),
+ SubControlVolumeFace::FaceType::frontal,
+ SubControlVolumeFace::BoundaryType::physicalBoundary
+ );
+
+ globalScvfIndices.push_back(globalScvfIdx);
+ ++globalScvfIdx;
++localScvfIdx;
hasBoundaryScvf_[eIdx] = true;
}
@@ -494,9 +474,6 @@ private:
}
connectivityMap_.update(*this);
-
- // allow to free the non-needed memory
- lateralOrthogonalScvf_.resize(scvfs_.size());
}
bool onDomainBoundary_(const typename GridView::Intersection& intersection) const
@@ -525,7 +502,6 @@ private:
GridIndexType outSideBoundaryVolVarIdx_;
std::vector<bool> hasBoundaryScvf_;
- std::vector<GridIndexType> lateralOrthogonalScvf_;
std::vector<std::vector<GridIndexType>> scvfIndicesOfElement_;
// a map for periodic boundary vertices
@@ -576,6 +552,8 @@ public:
using GridView = GV;
//! export the geometry helper type
using GeometryHelper = typename Traits::GeometryHelper;
+ //! export the local intersection mapper
+ using LocalIntersectionMapper = typename Traits::LocalIntersectionMapper;
//! export static information
using StaticInformation = typename Traits::StaticInfo;
//! export the type of extrusion
@@ -667,10 +645,6 @@ public:
const std::unordered_map<GridIndexType, GridIndexType>& periodicVertexMap() const
{ return periodicFaceMap_; }
- //! get the global index of the orthogonal face sharing a common entity
- GridIndexType lateralOrthogonalScvf(const SubControlVolumeFace& scvf) const
- { return lateralOrthogonalScvf_[scvf.index()]; }
-
private:
void update_()
@@ -678,28 +652,25 @@ private:
intersectionMapper_.update(this->gridView());
// clear local data
- numScvs_ = 0;
numScvf_ = 0;
numBoundaryScv_ = 0;
numBoundaryScvf_ = 0;
hasBoundaryScvf_.clear();
scvfIndicesOfElement_.clear();
- lateralOrthogonalScvf_.clear();
outsideVolVarIndices_.clear();
// determine size of containers
const auto numElements = this->gridView().size(0);
scvfIndicesOfElement_.resize(numElements);
hasBoundaryScvf_.resize(numElements, false);
- // on frontal + maybe boundary per face and 2*(dim-1) laterals per frontal
- lateralOrthogonalScvf_.resize(numElements*(2*dim*(2 + 2*(dim-1))));
+ numScvs_ = numElements*numScvsPerElement;
GeometryHelper geometryHelper(this->gridView());
// get the global scv indices first
GridIndexType scvfIdx = 0;
- GridIndexType neighborVolVarIdx = numElements*numScvsPerElement;
+ GridIndexType neighborVolVarIdx = numScvs_;
for (const auto& element : elements(this->gridView()))
{
@@ -707,43 +678,43 @@ private:
assert(numScvsPerElement == element.subEntities(1));
// the element-wise index sets for finite volume geometry
- std::vector<GridIndexType> scvfsIndexSet;
- scvfsIndexSet.reserve(1/*frontal in element*/ + 1 /*frontal on boundary*/ + numLateralScvfsPerScv);
+ auto& globalScvfIndices = scvfIndicesOfElement_[eIdx];
+ globalScvfIndices.reserve(maxNumScvfsPerElement);
+ globalScvfIndices.resize(minNumScvfsPerElement);
// keep track of frontal boundary scvfs
std::size_t numFrontalBoundaryScvfs = 0;
- // a temporary map to store pairs of common entities
- std::unordered_map<GridIndexType, Dune::ReservedVector<GridIndexType, 2>> commonEntityIdxToScvfsMap;
+ using LocalIntersectionIndexMapper = FaceCenteredStaggeredLocalIntersectionIndexMapper<GridView>;
+ LocalIntersectionIndexMapper localIsMapper;
+ localIsMapper.update(this->gridView(), element);
- SmallLocalIndexType localScvIdx = 0;
for (const auto& intersection : intersections(this->gridView(), element))
{
- ++numScvs_;
+ const auto localScvIdx = localIsMapper.realToRefIdx(intersection.indexInInside());
+ auto localScvfIdx = localScvIdx*(1 + numLateralScvfsPerScv);
+ assert(localIsMapper.refToRealIdx(localScvIdx) == intersection.indexInInside());
// the frontal sub control volume face at the element center
- scvfsIndexSet.push_back(scvfIdx++);
+ globalScvfIndices[localScvfIdx] = scvfIdx++;
+ ++localScvfIdx;
if constexpr(dim > 1)
{
// the lateral sub control volume faces
- for (const auto lateralFacetIndex : geometryHelper.localLaterFaceIndices(localScvIdx))
+ for (const auto lateralFacetIndex : Dune::transformedRangeView(geometryHelper.localLaterFaceIndices(localScvIdx),
+ [&](auto idx) { return localIsMapper.refToRealIdx(idx) ;})
+ )
{
- const auto& lateralIntersection = geometryHelper.intersection(lateralFacetIndex, element);
- if (onProcessorBoundary_(lateralIntersection))
- continue;
-
- if (onDomainBoundary_(lateralIntersection))
+ if (onDomainBoundary_(geometryHelper.intersection(lateralFacetIndex, element)))
+ {
outsideVolVarIndices_[scvfIdx] = neighborVolVarIdx++;
+ ++numBoundaryScvf_;
+ hasBoundaryScvf_[eIdx] = true;
+ }
- scvfsIndexSet.push_back(scvfIdx);
-
- const auto commonEntityIdx = geometryHelper.globalCommonEntityIndex(
- element, localScvIdx, lateralFacetIndex
- );
- commonEntityIdxToScvfsMap[commonEntityIdx].push_back(scvfIdx);
-
- ++scvfIdx;
+ globalScvfIndices[localScvfIdx] = scvfIdx++;
+ ++localScvfIdx;
}
}
@@ -783,43 +754,11 @@ private:
++otherIntersectionLocalIdx;
}
}
-
- // the lateral sub control volume faces
- const auto lateralFaceIndices = geometryHelper.localLaterFaceIndices(localScvIdx);
- for (const auto lateralFacetIndex : lateralFaceIndices)
- {
- const auto& lateralIntersection = geometryHelper.intersection(lateralFacetIndex, element);
- if (onProcessorBoundary_(lateralIntersection))
- continue;
-
- if (onDomainBoundary_(lateralIntersection))
- {
- ++numBoundaryScvf_;
- hasBoundaryScvf_[eIdx] = true;
- }
- } // end loop over lateral facets
-
- ++localScvIdx;
}
// add global indices of frontal boundary scvfs last
for (std::size_t i = 0; i < numFrontalBoundaryScvfs; ++i)
- scvfsIndexSet.push_back(scvfIdx++);
-
- // Save the scv indices belonging to this element to build up fv element geometries fast
- scvfIndicesOfElement_[eIdx] = std::move(scvfsIndexSet);
-
- // create the bi-directional map
- for (const auto& [_, scvfs] : commonEntityIdxToScvfsMap)
- {
- // TODO: this maybe be less than 2 if there is a processor boundary
- // are we sure that the lateral orthogonal scvf is then never needed?
- if (scvfs.size() == 2)
- {
- lateralOrthogonalScvf_[scvfs[0]] = scvfs[1];
- lateralOrthogonalScvf_[scvfs[1]] = scvfs[0];
- }
- }
+ globalScvfIndices.push_back(scvfIdx++);
}
// set number of subcontrolvolume faces
@@ -854,7 +793,6 @@ private:
std::size_t numBoundaryScvf_;
std::vector<bool> hasBoundaryScvf_;
- std::vector<GridIndexType> lateralOrthogonalScvf_;
std::vector<std::vector<GridIndexType>> scvfIndicesOfElement_;
// a map for periodic boundary vertices
diff --git a/dumux/discretization/facecentered/staggered/subcontrolvolumeface.hh b/dumux/discretization/facecentered/staggered/subcontrolvolumeface.hh
index d7e6c6e4ba..e5029bfbeb 100644
--- a/dumux/discretization/facecentered/staggered/subcontrolvolumeface.hh
+++ b/dumux/discretization/facecentered/staggered/subcontrolvolumeface.hh
@@ -77,7 +77,8 @@ public:
using Traits = T;
using GlobalPosition = typename T::GlobalPosition;
- enum class FaceType {frontal, lateral};
+ enum class FaceType : SmallLocalIndexType {frontal, lateral};
+ enum class BoundaryType : SmallLocalIndexType {interior, physicalBoundary, processorBoundary};
FaceCenteredStaggeredSubControlVolumeFace() = default;
@@ -89,7 +90,7 @@ public:
const GridIndexType globalScvfIdx,
const GlobalPosition& unitOuterNormal,
const FaceType faceType,
- const bool boundary)
+ const BoundaryType boundaryType)
: globalScvIndices_(globalScvIndices)
, localScvfIdx_(localScvfIdx)
, globalScvfIdx_(globalScvfIdx)
@@ -97,17 +98,17 @@ public:
, normalAxis_(Dumux::normalAxis(unitOuterNormal))
, outerNormalSign_(sign(unitOuterNormal[normalAxis_]))
, faceType_(faceType)
- , boundary_(boundary)
+ , boundaryType_(boundaryType)
{
assert(faceType == FaceType::frontal);
- center_ = boundary ? intersectionGeometry.center() : elementGeometry.center();
+ center_ = boundary() ? intersectionGeometry.center() : elementGeometry.center();
ipGlobal_ = center_;
- if (!boundary)
+ if (!boundary())
outerNormalSign_ *= -1.0;
// the corners (coincide with intersection corners for boundary scvfs)
- const auto frontalOffSet = boundary ? GlobalPosition(0.0)
+ const auto frontalOffSet = boundary() ? GlobalPosition(0.0)
: intersectionGeometry.center() - elementGeometry.center();
for (int i = 0; i < corners_.size(); ++i)
@@ -124,7 +125,7 @@ public:
const GridIndexType globalScvfIdx,
const GlobalPosition& unitOuterNormal,
const FaceType faceType,
- const bool boundary)
+ const BoundaryType boundaryType)
: globalScvIndices_(globalScvIndices)
, localScvfIdx_(localScvfIdx)
, globalScvfIdx_(globalScvfIdx)
@@ -132,7 +133,7 @@ public:
, normalAxis_(Dumux::normalAxis(unitOuterNormal))
, outerNormalSign_(sign(unitOuterNormal[normalAxis_]))
, faceType_(faceType)
- , boundary_(boundary)
+ , boundaryType_(boundaryType)
{
assert(faceType == FaceType::lateral);
const auto shift = intersectionGeometry.center() - elementGeometry.center();
@@ -190,7 +191,10 @@ public:
{ return faceType_; }
bool boundary() const
- { return boundary_; }
+ { return boundaryType_ == BoundaryType::physicalBoundary; }
+
+ bool processorBoundary() const
+ { return boundaryType_ == BoundaryType::processorBoundary; }
bool isFrontal() const
{ return faceType_ == FaceType::frontal; }
@@ -232,7 +236,7 @@ private:
SmallLocalIndexType normalAxis_;
std::int_least8_t outerNormalSign_;
FaceType faceType_;
- bool boundary_;
+ BoundaryType boundaryType_;
};
} // end namespace Dumux
--
GitLab