From ae8ec14ed58855c81d005433d3966846576b2715 Mon Sep 17 00:00:00 2001
From: Ned Coltman <edward.coltman@iws.uni-stuttgart.de>
Date: Wed, 23 Jan 2019 17:14:20 +0100
Subject: [PATCH] [axisdata][pairdata] change axis and pair data structs to
store dofs and distances with arrays, not vectors
---
.../staggered/freeflow/connectivitymap.hh | 14 ++-
.../freeflow/staggeredgeometryhelper.hh | 93 +++++++------------
.../freeflow/subcontrolvolumeface.hh | 8 +-
.../staggered/fvgridgeometry.hh | 14 ++-
4 files changed, 64 insertions(+), 65 deletions(-)
diff --git a/dumux/discretization/staggered/freeflow/connectivitymap.hh b/dumux/discretization/staggered/freeflow/connectivitymap.hh
index fc0b40e355..e43321f1c0 100644
--- a/dumux/discretization/staggered/freeflow/connectivitymap.hh
+++ b/dumux/discretization/staggered/freeflow/connectivitymap.hh
@@ -92,6 +92,13 @@ public:
}
}
+ //! \brief Set the order needed by the scheme
+ void setStencilOrder(const int stencilOrder)
+ {
+ stencilOrder_ = stencilOrder;
+ }
+
+
//! Returns the stencil of a cell center dof w.r.t. other cell center dofs
const std::vector<GridIndexType>& operator() (CellCenterIdxType, CellCenterIdxType, const GridIndexType globalI) const
{
@@ -179,7 +186,7 @@ private:
{
if(stencil.empty())
{
- for(int i = 0; i < scvf.axisData().inAxisBackwardDofs.size(); i++)
+ for(int i = 0; i < stencilOrder_ - 1; i++)
{
if(scvf.hasBackwardNeighbor(i))
{
@@ -190,7 +197,7 @@ private:
stencil.push_back(scvf.axisData().selfDof);
stencil.push_back(scvf.axisData().oppositeDof);
- for(int i = 0; i < scvf.axisData().inAxisForwardDofs.size(); i++)
+ for(int i = 0; i < stencilOrder_ - 1; i++)
{
if(scvf.hasForwardNeighbor(i))
{
@@ -207,7 +214,7 @@ private:
stencil.push_back(data.normalPair.second);
// add parallel dofs
- for (int i = 0; i < data.parallelDofs.size(); i++)
+ for (int i = 0; i < stencilOrder_ /*data.parallelDofs.size()*/; i++)
{
if(!(data.parallelDofs[i] < 0))
{
@@ -221,6 +228,7 @@ private:
CellCenterToFaceMap cellCenterToFaceMap_;
FaceToCellCenterMap faceToCellCenterMap_;
FaceToFaceMap faceToFaceMap_;
+ int stencilOrder_;
};
} // end namespace Dumux
diff --git a/dumux/discretization/staggered/freeflow/staggeredgeometryhelper.hh b/dumux/discretization/staggered/freeflow/staggeredgeometryhelper.hh
index 7d87aa8825..6990c0daed 100644
--- a/dumux/discretization/staggered/freeflow/staggeredgeometryhelper.hh
+++ b/dumux/discretization/staggered/freeflow/staggeredgeometryhelper.hh
@@ -30,6 +30,7 @@
#include <dumux/common/math.hh>
#include <type_traits>
#include <algorithm>
+#include <array>
namespace Dumux
{
@@ -38,11 +39,11 @@ namespace Dumux
* \ingroup StaggeredDiscretization
* \brief Parallel Data stored per sub face
*/
-template<class Scalar, class GlobalPosition>
+template<class Scalar, class GlobalPosition, int geometryOrder>
struct PairData
{
- std::vector<int> parallelDofs;
- std::vector<Scalar> parallelDistances;
+ std::array<int, geometryOrder> parallelDofs;
+ std::array<Scalar, geometryOrder+1> parallelDistances;
std::pair<signed int,signed int> normalPair;
int localNormalFaceIdx;
Scalar normalDistance;
@@ -53,16 +54,16 @@ struct PairData
* \ingroup StaggeredDiscretization
* \brief In Axis Data stored per sub face
*/
-template<class Scalar>
+template<class Scalar, int geometryOrder>
struct AxisData
{
int selfDof;
int oppositeDof;
- std::vector<int> inAxisForwardDofs;
- std::vector<int> inAxisBackwardDofs;
+ std::array<int, geometryOrder-1> inAxisForwardDofs;
+ std::array<int, geometryOrder-1> inAxisBackwardDofs;
Scalar selfToOppositeDistance;
- std::vector<Scalar> inAxisForwardDistances;
- std::vector<Scalar> inAxisBackwardDistances;
+ std::array<Scalar, geometryOrder-1> inAxisForwardDistances;
+ std::array<Scalar, geometryOrder-1> inAxisBackwardDistances;
};
/*!
@@ -105,7 +106,6 @@ class FreeFlowStaggeredGeometryHelper
static constexpr int numFacetSubEntities = (dim == 2) ? 2 : 4;
static constexpr int numfacets = dimWorld * 2;
static constexpr int numPairs = 2 * (dimWorld - 1);
- static int order_;
public:
@@ -162,17 +162,7 @@ public:
return Dumux::directionIndex(std::move(intersection.centerUnitOuterNormal()));
}
- //! \brief Returns the order needed by the scheme
- static int order()
- {
- return order_;
- }
-
- //! \brief Set the order needed by the scheme
- static void setOrder(const int order)
- {
- order_ = order;
- }
+ static constexpr int geometryOrder = 2;
private:
/*!
@@ -180,24 +170,14 @@ private:
*/
void fillAxisData_()
{
- unsigned int numForwardBackwardAxisDofs = order_ - 1;
+ unsigned int numForwardBackwardAxisDofs = axisData_.inAxisForwardDofs.size();
const auto inIdx = intersection_.indexInInside();
const auto oppIdx = localOppositeIdx_(inIdx);
- // Clear the containers before filling them
- this->axisData_.inAxisForwardDofs.clear();
- this->axisData_.inAxisBackwardDofs.clear();
- this->axisData_.inAxisForwardDistances.clear();
- this->axisData_.inAxisBackwardDistances.clear();
-
- // initialize values that could remain unitialized if the intersection lies on a boundary
- for(int i = 0; i < numForwardBackwardAxisDofs; i++)
- {
- this->axisData_.inAxisForwardDofs.push_back(-1);
- this->axisData_.inAxisBackwardDofs.push_back(-1);
- this->axisData_.inAxisForwardDistances.push_back(0);
- this->axisData_.inAxisBackwardDistances.push_back(0);
- }
+ this->axisData_.inAxisForwardDofs.fill(-1);
+ this->axisData_.inAxisBackwardDofs.fill(-1);
+ this->axisData_.inAxisForwardDistances.fill(0);
+ this->axisData_.inAxisBackwardDistances.fill(0);
// Set the self Dof
this->axisData_.selfDof = gridView_.indexSet().subIndex(this->intersection_.inside(), inIdx, codimIntersection);
@@ -215,9 +195,8 @@ private:
auto selfFace = getFacet_(inIdx, element_);
if(intersection_.neighbor())
{
- if (order_ > 1)
- inAxisForwardElementStack.push(intersection_.outside());
- bool keepStackingForward = (inAxisForwardElementStack.size() < order_-1);
+ inAxisForwardElementStack.push(intersection_.outside());
+ bool keepStackingForward = (inAxisForwardElementStack.size() < numForwardBackwardAxisDofs);
while(keepStackingForward)
{
auto e = inAxisForwardElementStack.top();
@@ -228,7 +207,7 @@ private:
if( intersection.neighbor())
{
inAxisForwardElementStack.push(intersection.outside());
- keepStackingForward = (inAxisForwardElementStack.size() < order_-1);
+ keepStackingForward = (inAxisForwardElementStack.size() < numForwardBackwardAxisDofs);
}
else
{
@@ -238,6 +217,7 @@ private:
}
}
}
+
std::vector<GlobalPosition> forwardFaceCoordinates(inAxisForwardElementStack.size(), selfFace.geometry().center());
while(!inAxisForwardElementStack.empty())
{
@@ -247,6 +227,7 @@ private:
forwardFaceCoordinates[forwardIdx] = selfFace.geometry().center();
inAxisForwardElementStack.pop();
}
+
for(int i = 0; i< forwardFaceCoordinates.size(); i++)
{
if(i == 0)
@@ -268,9 +249,8 @@ private:
{
if(intersection.neighbor())
{
- if (order_ > 1)
- inAxisBackwardElementStack.push(intersection.outside());
- bool keepStackingBackward = (inAxisBackwardElementStack.size() < order_-1);
+ inAxisBackwardElementStack.push(intersection.outside());
+ bool keepStackingBackward = (inAxisBackwardElementStack.size() < numForwardBackwardAxisDofs);
while(keepStackingBackward)
{
auto e = inAxisBackwardElementStack.top();
@@ -282,7 +262,7 @@ private:
if( intersectionOut.neighbor())
{
inAxisBackwardElementStack.push(intersectionOut.outside());
- keepStackingBackward = (inAxisBackwardElementStack.size() < order_-1);
+ keepStackingBackward = (inAxisBackwardElementStack.size() < numForwardBackwardAxisDofs);
}
else
{
@@ -294,6 +274,7 @@ private:
}
}
}
+
std::vector<GlobalPosition> backwardFaceCoordinates(inAxisBackwardElementStack.size(), oppFace.geometry().center());
while(!inAxisBackwardElementStack.empty())
{
@@ -303,6 +284,7 @@ private:
backwardFaceCoordinates[backwardIdx] = oppFace.geometry().center();
inAxisBackwardElementStack.pop();
}
+
for(int i = 0; i< backwardFaceCoordinates.size(); i++)
{
if(i == 0)
@@ -324,20 +306,18 @@ private:
// initialize values that could remain unitialized if the intersection lies on a boundary
for(auto& data : pairData_)
{
- int numParallelDofs = order_;
-
// parallel Dofs
- data.parallelDofs.clear();
- data.parallelDofs.resize(numParallelDofs, -1);
+ data.parallelDofs.fill(-1);
// parallel Distances
- data.parallelDistances.clear();
- data.parallelDistances.resize(numParallelDofs + 1, 0.0);
+ data.parallelDistances.fill(0.0);
// outer normals
data.normalPair.second = -1;
}
+ unsigned int numParallelFaces = pairData_[0].parallelDistances.size();
+
// set basic global positions
const auto& elementCenter = this->element_.geometry().center();
const auto& FacetCenter = intersection_.geometry().center();
@@ -388,7 +368,6 @@ private:
}
}
-
// get the parallel Dofs
const int parallelLocalIdx = intersection_.indexInInside();
int numPairParallelIdx = 0;
@@ -402,7 +381,8 @@ private:
auto parallelAxisIdx = directionIndex(intersection);
auto localNormalIntersectionIndex = intersection.indexInInside();
parallelElementStack.push(element_);
- bool keepStacking = (parallelElementStack.size() < order_+1);
+
+ bool keepStacking = (parallelElementStack.size() < numParallelFaces);
while(keepStacking)
{
auto e = parallelElementStack.top();
@@ -413,7 +393,7 @@ private:
if( normalIntersection.neighbor() )
{
parallelElementStack.push(normalIntersection.outside());
- keepStacking = (parallelElementStack.size() < order_+1);
+ keepStacking = (parallelElementStack.size() < numParallelFaces);
}
else
{
@@ -422,6 +402,7 @@ private:
}
}
}
+
while(!parallelElementStack.empty())
{
if(parallelElementStack.size() > 1)
@@ -431,6 +412,7 @@ private:
this->pairData_[numPairParallelIdx].parallelDistances[parallelElementStack.size()-1] = setParallelPairDistances_(parallelElementStack.top(), parallelAxisIdx);
parallelElementStack.pop();
}
+
}
else
{
@@ -544,14 +526,11 @@ private:
const Element element_; //!< The respective element
const typename Element::Geometry elementGeometry_; //!< Reference to the element geometry
const GridView gridView_; //!< The grid view
- AxisData<Scalar> axisData_; //!< Data related to forward and backward faces
- std::array<PairData<Scalar, GlobalPosition>, numPairs> pairData_; //!< Collection of pair information related to normal and parallel faces
+ AxisData<Scalar, geometryOrder> axisData_; //!< Data related to forward and backward faces
+ std::array<PairData<Scalar, GlobalPosition, geometryOrder>, numPairs> pairData_; //!< Collection of pair information related to normal and parallel faces
std::vector<GlobalPosition> innerNormalFacePos_;
};
-template<class GridView>
-int FreeFlowStaggeredGeometryHelper<GridView>::order_ = 1;
-
} // end namespace Dumux
#endif
diff --git a/dumux/discretization/staggered/freeflow/subcontrolvolumeface.hh b/dumux/discretization/staggered/freeflow/subcontrolvolumeface.hh
index 0fc4fc8e9e..0fa6efc607 100644
--- a/dumux/discretization/staggered/freeflow/subcontrolvolumeface.hh
+++ b/dumux/discretization/staggered/freeflow/subcontrolvolumeface.hh
@@ -51,12 +51,14 @@ class FreeFlowStaggeredSubControlVolumeFace
using Geometry = typename T::Geometry;
using GridIndexType = typename IndexTraits<GV>::GridIndex;
using LocalIndexType = typename IndexTraits<GV>::LocalIndex;
+ using GeometryHelper = FreeFlowStaggeredGeometryHelper<GV>;
using Scalar = typename T::Scalar;
static const int dim = Geometry::mydimension;
static const int dimworld = Geometry::coorddimension;
static constexpr int numPairs = 2 * (dimworld - 1);
+ static constexpr int geometryOrder = GeometryHelper::geometryOrder;
public:
using GlobalPosition = typename T::GlobalPosition;
@@ -204,7 +206,7 @@ public:
}
//! Returns the data for one sub face
- const PairData<Scalar, GlobalPosition>& pairData(const int idx) const
+ const PairData<Scalar, GlobalPosition, geometryOrder>& pairData(const int idx) const
{
return pairData_[idx];
}
@@ -323,8 +325,8 @@ private:
int dofIdx_;
Scalar selfToOppositeDistance_;
- AxisData<Scalar> axisData_;
- std::array<PairData<Scalar, GlobalPosition>, numPairs> pairData_;
+ AxisData<Scalar, geometryOrder> axisData_;
+ std::array<PairData<Scalar, GlobalPosition, geometryOrder>, numPairs> pairData_;
int localFaceIdx_;
unsigned int dirIdx_;
diff --git a/dumux/discretization/staggered/fvgridgeometry.hh b/dumux/discretization/staggered/fvgridgeometry.hh
index 7ca4c9fac7..ed47bfce25 100644
--- a/dumux/discretization/staggered/fvgridgeometry.hh
+++ b/dumux/discretization/staggered/fvgridgeometry.hh
@@ -228,7 +228,9 @@ public:
DUNE_THROW(Dune::InvalidStateException, "The staggered discretization method needs at least an overlap of 1 for parallel computations. "
<< " Set the parameter \"Grid.Overlap\" in the input file.");
if (hasParamInGroup(paramGroup, "Discretization.TvdApproach"))
- GeometryHelper::setOrder(2);
+ stencilOrder_ = 2;
+ else
+ stencilOrder_ = 1;
}
//! The total number of sub control volumes
@@ -343,6 +345,7 @@ public:
}
// build the connectivity map for an effecient assembly
+ connectivityMap_.setStencilOrder(stencilOrder_);
connectivityMap_.update(*this);
}
@@ -414,6 +417,7 @@ private:
// mappers
ConnectivityMap connectivityMap_;
IntersectionMapper intersectionMapper_;
+ int stencilOrder_;
std::vector<SubControlVolume> scvs_;
std::vector<SubControlVolumeFace> scvfs_;
@@ -473,7 +477,7 @@ public:
using FVGridGeometryTuple = std::tuple< CellCenterFVGridGeometry<ThisType>, FaceFVGridGeometry<ThisType> >;
//! Constructor
- StaggeredFVGridGeometry(const GridView& gridView)
+ StaggeredFVGridGeometry(const GridView& gridView, const std::string& paramGroup = "")
: ParentType(gridView)
, intersectionMapper_(gridView)
{
@@ -481,6 +485,10 @@ public:
if (!CheckOverlapSize<DiscretizationMethod::staggered>::isValid(gridView))
DUNE_THROW(Dune::InvalidStateException, "The staggered discretization method needs at least an overlap of 1 for parallel computations. "
<< " Set the parameter \"Grid.Overlap\" in the input file.");
+ if (hasParamInGroup(paramGroup, "Discretization.TvdApproach"))
+ stencilOrder_ = 2;
+ else
+ stencilOrder_ = 1;
}
//! update all fvElementGeometries (do this again after grid adaption)
@@ -533,6 +541,7 @@ public:
}
// build the connectivity map for an effecient assembly
+ connectivityMap_.setStencilOrder(stencilOrder_);
connectivityMap_.update(*this);
}
@@ -631,6 +640,7 @@ private:
// mappers
ConnectivityMap connectivityMap_;
IntersectionMapper intersectionMapper_;
+ int stencilOrder_;
//! vectors that store the global data
std::vector<std::vector<GridIndexType>> scvfIndicesOfScv_;
--
GitLab