Skip to content
GitLab
Commit f18d3401 authored by Kilian Weishaupt's avatar Kilian Weishaupt
Browse files

[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
parent 430fe0e6
Expand all Hide whitespace changes
Inline Side-by-side
dumux/discretization/facecentered/staggered/fvelementgeometry.hh
View file @ f18d3401
......@@ -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_;
......
dumux/discretization/facecentered/staggered/subcontrolvolumeface.hh
View file @ f18d3401
......@@ -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
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment