Commit eb258652 authored by Dennis Gläser's avatar Dennis Gläser Committed by Timo Koch
Browse files

make subcontrolvolume entities base classes simple interfaces

The scv and scvf base classes now get the actual implementation
as a template parameter and serve as an interface for scv and scvf
implementations. Furthermore, the box scv and scvf does not store a geometry
anymore, but only the corner points. The geometry is built on the fly if requested.
The BoxGeometryHelper is adjusted to this new concept, as well as the global and
element-fvgeometry.

(cherry picked from commit 8b55b5d6db4414e57072ed05b8903a5ec49b0ee5)
parent 452c0e03
......@@ -104,6 +104,8 @@ private:
using ScvfGeometry = Dune::MultiLinearGeometry<Scalar, dim-1, dimWorld>;
using GlobalPosition = typename ScvGeometry::GlobalCoordinate;
using PointVector = std::vector<GlobalPosition>;
using Element = typename GridView::template Codim<0>::Entity;
using Intersection = typename GridView::Intersection;
......@@ -112,129 +114,129 @@ private:
public:
BoxGeometryHelper(const typename Element::Geometry& geometry)
: vertexOffset(1),
faceOffset(1 + geometry.corners())
{
// extract the corners of the sub control volumes
const auto& referenceElement = ReferenceElements::general(geometry.type());
// the element center
p.emplace_back(geometry.center());
// vertices
for (int i = 0; i < geometry.corners(); ++i)
v.emplace_back(geometry.corner(i));
p.emplace_back(geometry.corner(i));
// face midpoints
for (int i = 0; i < referenceElement.size(1); ++i)
f.emplace_back(geometry.global(referenceElement.position(i, 1)));
c = geometry.center();
p.emplace_back(geometry.global(referenceElement.position(i, 1)));
corners = geometry.corners();
}
//! Create the sub control volume geometries
std::vector<ScvGeometry> createScvGeometries()
//! Create a vector with the scv corners
PointVector getScvCorners(unsigned int localScvIdx)
{
std::vector<ScvGeometry> scvGeometries;
// sub control volume geometries in 2D are always quadrilaterals
Dune::GeometryType scvGeometryType;
scvGeometryType.makeQuadrilateral();
// proceed according to number of corners
// proceed according to number of corners of the element
switch (corners)
{
case 3: // triangle
{
scvGeometries.reserve(3);
// the sub control volumes
scvGeometries.emplace_back(scvGeometryType, std::vector<GlobalPosition>({v[0], f[0], f[1], c}));
scvGeometries.emplace_back(scvGeometryType, std::vector<GlobalPosition>({f[0], v[1], c, f[2]}));
scvGeometries.emplace_back(scvGeometryType, std::vector<GlobalPosition>({f[1], c, v[2], f[2]}));
break;
static const unsigned int triangleScvMap[3][4] =
{
{vertexOffset+0, faceOffset+0, faceOffset+1, 0},
{vertexOffset+1, faceOffset+2, faceOffset+0, 0},
{vertexOffset+2, faceOffset+1, faceOffset+2, 0}
};
return PointVector( {p[triangleScvMap[localScvIdx][0]],
p[triangleScvMap[localScvIdx][1]],
p[triangleScvMap[localScvIdx][2]],
p[triangleScvMap[localScvIdx][3]]} );
}
case 4: // quadrilateral
{
scvGeometries.reserve(4);
// the sub control volumes
scvGeometries.emplace_back(scvGeometryType, std::vector<GlobalPosition>({v[0], f[2], f[0], c}));
scvGeometries.emplace_back(scvGeometryType, std::vector<GlobalPosition>({f[2], v[1], c, f[1]}));
scvGeometries.emplace_back(scvGeometryType, std::vector<GlobalPosition>({f[0], c, v[2], f[3]}));
scvGeometries.emplace_back(scvGeometryType, std::vector<GlobalPosition>({c, f[1], f[3], v[3]}));
break;
static const unsigned int quadScvMap[4][4] =
{
{vertexOffset+0, faceOffset+2, faceOffset+0, 0},
{vertexOffset+1, faceOffset+1, faceOffset+2, 0},
{vertexOffset+2, faceOffset+0, faceOffset+3, 0},
{vertexOffset+3, faceOffset+3, faceOffset+1, 0}
};
return PointVector( {p[quadScvMap[localScvIdx][0]],
p[quadScvMap[localScvIdx][1]],
p[quadScvMap[localScvIdx][2]],
p[quadScvMap[localScvIdx][3]]} );
}
default:
DUNE_THROW(Dune::NotImplemented, "Box scv geometries for dim=" << dim
<< " dimWorld=" << dimWorld
<< " corners=" << corners);
}
return scvGeometries;
}
//! Create the sub control volume face geometries
std::vector<ScvfGeometry> createScvfGeometries()
//! Create a vector with the corners of sub control volume faces
PointVector getScvfCorners(unsigned int localScvfIdx)
{
std::vector<ScvfGeometry> scvfGeometries;
// sub control volume face geometries in 2D are always lines
Dune::GeometryType scvfGeometryType;
scvfGeometryType.makeLine();
// proceed according to number of corners
switch (corners)
{
case 3: // triangle
{
scvfGeometries.reserve(3);
// the sub control volume faces
scvfGeometries.emplace_back(scvfGeometryType, std::vector<GlobalPosition>({f[0], c}));
scvfGeometries.emplace_back(scvfGeometryType, std::vector<GlobalPosition>({f[1], c}));
scvfGeometries.emplace_back(scvfGeometryType, std::vector<GlobalPosition>({f[2], c}));
break;
static const unsigned int triangleScvfMap[3][2] =
{
{0, faceOffset+0},
{faceOffset+1, 0},
{0, faceOffset+2}
};
return PointVector( {p[triangleScvfMap[localScvfIdx][0]],
p[triangleScvfMap[localScvfIdx][1]]} );
}
case 4: // quadrilateral
{
scvfGeometries.reserve(4);
// the sub control volume faces
scvfGeometries.emplace_back(scvfGeometryType, std::vector<GlobalPosition>({f[0], c}));
scvfGeometries.emplace_back(scvfGeometryType, std::vector<GlobalPosition>({f[1], c}));
scvfGeometries.emplace_back(scvfGeometryType, std::vector<GlobalPosition>({f[2], c}));
scvfGeometries.emplace_back(scvfGeometryType, std::vector<GlobalPosition>({f[3], c}));
break;
static const unsigned int quadScvfMap[4][2] =
{
{faceOffset+0, 0},
{0, faceOffset+1},
{0, faceOffset+2},
{faceOffset+3, 0}
};
return PointVector( {p[quadScvfMap[localScvfIdx][0]],
p[quadScvfMap[localScvfIdx][1]]} );
}
default:
DUNE_THROW(Dune::NotImplemented, "Box scv geometries for dim=" << dim
DUNE_THROW(Dune::NotImplemented, "Box scvf geometries for dim=" << dim
<< " dimWorld=" << dimWorld
<< " corners=" << corners);
}
return scvfGeometries;
}
//! Create the sub control volume face geometries on the boundary
static std::vector<ScvfGeometry> createBoundaryScvfGeometries(const typename Intersection::Geometry& geometry)
static PointVector getBoundaryScvfCorners(const typename Intersection::Geometry& geometry,
unsigned int idxOnIntersection)
{
return {ScvfGeometry(Dune::GeometryType(1), std::vector<GlobalPosition>({geometry.corner(0), geometry.center()})),
ScvfGeometry(Dune::GeometryType(1), std::vector<GlobalPosition>({geometry.center(), geometry.corner(1)}))};
if (idxOnIntersection == 0)
return PointVector({geometry.corner(0), geometry.center()});
else if (idxOnIntersection == 1)
return PointVector({geometry.center(), geometry.corner(1)});
else
DUNE_THROW(Dune::InvalidStateException, "local index exceeds the number of corners of 2d intersections");
}
//! get scvf normal vector for dim == 2, dimworld == 3
template <int w = dimWorld>
static typename std::enable_if<w == 3, GlobalPosition>::type
normal(const typename Element::Geometry& geometry,
const ScvfGeometry& scvfGeometry)
const PointVector& scvfCorners)
{
const auto v1 = geometry.corner(1) - geometry.corner(0);
const auto v2 = geometry.corner(2) - geometry.corner(0);
const auto v3 = Dumux::crossProduct(v1, v2);
const auto t = scvfGeometry.corner(1) - scvfGeometry.corner(0);
const auto t = scvfCorners[1] - scvfCorners[0];
GlobalPosition normal = Dumux::crossProduct(v3, t);
normal /= normal.two_norm();
return normal;
......@@ -244,18 +246,44 @@ public:
template <int w = dimWorld>
static typename std::enable_if<w == 2, GlobalPosition>::type
normal(const typename Element::Geometry& geometry,
const ScvfGeometry& scvfGeometry)
const PointVector& scvfCorners)
{
const auto t = scvfGeometry.corner(1) - scvfGeometry.corner(0);
const auto t = scvfCorners[1] - scvfCorners[0];
GlobalPosition normal({-t[1], t[0]});
normal /= normal.two_norm();
return normal;
}
//! get scv volume for dim == 2, dimworld == 3
template <int w = dimWorld>
static typename std::enable_if<w == 3, Scalar>::type
volume(const PointVector& scvCorners)
{
const auto v1 = scvCorners[1] - scvCorners[0];
const auto v2 = scvCorners[2] - scvCorners[0];
return Dumux::crossProduct(v1, v2).two_norm();
}
//! get scv volume for dim == 2, dimworld == 2
template <int w = dimWorld>
static typename std::enable_if<w == 2, Scalar>::type
volume(const PointVector& scvCorners)
{
const auto v1 = scvCorners[1] - scvCorners[0];
const auto v2 = scvCorners[2] - scvCorners[0];
return Dumux::crossProduct(v1, v2);
}
//! get scvf area
static Scalar area(const PointVector& scvfCorners)
{
return (scvfCorners[1] - scvfCorners[0]).two_norm();
}
private:
std::vector<GlobalPosition> v; // vertices
std::vector<GlobalPosition> f; // face midpoints
GlobalPosition c; // element center
std::vector<GlobalPosition> p; // the points needed for construction of the geometries
std::size_t corners; // number of element corners
const unsigned int vertexOffset; // offset for vertex positions in point vector
const unsigned int faceOffset; // offset for face center positions in point vector
};
//! A class to create sub control volume and sub control volume face geometries per element
......
......@@ -283,48 +283,48 @@ private:
// get the sub control volume geometries of this element
GeometryHelper geometryHelper(elementGeometry);
auto scvGeometries = geometryHelper.createScvGeometries();
auto scvfGeometries = geometryHelper.createScvfGeometries();
// construct the sub control volumes
IndexType scvLocalIdx = 0;
for (auto&& scvGeometry : scvGeometries)
for (unsigned int scvLocalIdx = 0; scvLocalIdx < elementGeometry.corners(); ++scvLocalIdx)
{
// get asssociated dof index
auto dofIdxGlobal = globalFvGeometry().problem_().vertexMapper().subIndex(element, scvLocalIdx, dim);
// get the corners and the volume of the scv
auto scvCorners = geometryHelper.getScvCorners(scvLocalIdx);
auto volume = geometryHelper.volume(scvCorners);
// add scv to the local container
scvs_.emplace_back(std::move(scvGeometry),
scvs_.emplace_back(std::move(scvCorners),
volume,
scvLocalIdx,
eIdx,
dofIdxGlobal);
// increment local counter
scvLocalIdx++;
}
// construct the sub control volume faces
IndexType scvfLocalIdx = 0;
for (auto&& scvfGeometry : scvfGeometries)
unsigned int scvfLocalIdx = 0;
for (; scvfLocalIdx < element.subEntities(1); ++scvfLocalIdx)
{
// find the local scv indices this scvf is connsected to
std::vector<IndexType> localScvIndices({static_cast<IndexType>(referenceElement.subEntity(scvfLocalIdx, dim-1, 0, dim)),
static_cast<IndexType>(referenceElement.subEntity(scvfLocalIdx, dim-1, 1, dim))});
// compute the scvf normal unit outer normal
auto normal = geometryHelper.normal(elementGeometry, scvfGeometry);
// get the corner points, the area, and the unit normal vector of the scv face
auto scvfCorners = geometryHelper.getScvfCorners(scvfLocalIdx);
auto area = geometryHelper.area(scvfCorners);
auto normal = geometryHelper.normal(elementGeometry, scvfCorners);
const auto v = elementGeometry.corner(localScvIndices[1]) - elementGeometry.corner(localScvIndices[0]);
const auto s = v*normal;
if (std::signbit(s))
normal *= -1;
scvfs_.emplace_back(std::move(scvfGeometry),
scvfs_.emplace_back(std::move(scvfCorners),
normal,
area,
scvfLocalIdx,
localScvIndices,
false);
// increment local counter
scvfLocalIdx++;
}
// construct the sub control volume faces on the domain boundary
......@@ -333,18 +333,21 @@ private:
if (intersection.boundary())
{
auto isGeometry = intersection.geometry();
auto boundaryScvfGeometries = geometryHelper.createBoundaryScvfGeometries(isGeometry);
IndexType boundaryFaceLocalIdx = 0;
for (auto&& scvfGeometry : boundaryScvfGeometries)
for (unsigned int isScvfLocalIdx = 0; isScvfLocalIdx < isGeometry.corners(); ++isScvfLocalIdx)
{
// find the scv this scvf is connected to
std::vector<IndexType> localScvIndices =
{static_cast<IndexType>(referenceElement.subEntity(intersection.indexInInside(), 1,
boundaryFaceLocalIdx++, dim))};
isScvfLocalIdx, dim))};
// get the corner points and the area of the boundary scv face
auto scvfCorners = geometryHelper.getBoundaryScvfCorners(isGeometry, isScvfLocalIdx);
auto area = geometryHelper.area(scvfCorners);
scvfs_.emplace_back(std::move(scvfGeometry),
scvfs_.emplace_back(std::move(scvfCorners),
intersection.centerUnitOuterNormal(),
area,
scvfLocalIdx,
localScvIndices,
true);
......
......@@ -135,49 +135,50 @@ public:
auto elementGeometry = element.geometry();
const auto& referenceElement = ReferenceElements::general(elementGeometry.type());
// get the sub control volume geometries of this element
// instantiate the geometry helper
GeometryHelper geometryHelper(elementGeometry);
auto scvGeometries = geometryHelper.createScvGeometries();
auto scvfGeometries = geometryHelper.createScvfGeometries();
// construct the sub control volumes
IndexType scvLocalIdx = 0;
scvs_[eIdx].reserve(scvGeometries.size());
for (auto&& scvGeometry : scvGeometries)
scvs_[eIdx].reserve(elementGeometry.corners());
for (unsigned int scvLocalIdx = 0; scvLocalIdx < elementGeometry.corners(); ++scvLocalIdx)
{
auto dofIdxGlobal = problem.vertexMapper().subIndex(element, scvLocalIdx, dim);
scvs_[eIdx].emplace_back(std::move(scvGeometry),
// get the corners and the volume of the scv
auto scvCorners = geometryHelper.getScvCorners(scvLocalIdx);
auto volume = geometryHelper.volume(scvCorners);
scvs_[eIdx].emplace_back(std::move(scvCorners),
volume,
scvLocalIdx,
eIdx,
dofIdxGlobal);
// increment local counter
scvLocalIdx++;
}
// construct the sub control volume faces
IndexType scvfLocalIdx = 0;
scvfs_[eIdx].reserve(scvfGeometries.size());
for (auto&& scvfGeometry : scvfGeometries)
auto numInteriorScvfs = element.subEntities(1);
unsigned int scvfLocalIdx = 0;
scvfs_[eIdx].reserve(numInteriorScvfs);
for (; scvfLocalIdx < numInteriorScvfs; ++scvfLocalIdx)
{
// find the global and local scv indices this scvf is belonging to
std::vector<IndexType> localScvIndices({static_cast<IndexType>(referenceElement.subEntity(scvfLocalIdx, dim-1, 0, dim)),
static_cast<IndexType>(referenceElement.subEntity(scvfLocalIdx, dim-1, 1, dim))});
// compute the scvf normal unit outer normal
auto normal = geometryHelper.normal(elementGeometry, scvfGeometry);
// get the corner points, the area, and the unit normal vector of the scv face
auto scvfCorners = geometryHelper.getScvfCorners(scvfLocalIdx);
auto area = geometryHelper.area(scvfCorners);
auto normal = geometryHelper.normal(elementGeometry, scvfCorners);
const auto v = elementGeometry.corner(localScvIndices[1]) - elementGeometry.corner(localScvIndices[0]);
const auto s = v*normal;
if (std::signbit(s))
normal *= -1;
scvfs_[eIdx].emplace_back(std::move(scvfGeometry),
scvfs_[eIdx].emplace_back(std::move(scvfCorners),
normal,
area,
scvfLocalIdx,
localScvIndices,
false);
// increment local counter
scvfLocalIdx++;
}
// construct the sub control volume faces on the domain boundary
......@@ -189,18 +190,21 @@ public:
// count
numScvf_ += isGeometry.corners();
numBoundaryScvf_ += isGeometry.corners();
auto boundaryScvfGeometries = geometryHelper.createBoundaryScvfGeometries(isGeometry);
IndexType isScvfLocalIdx = 0;
for (auto&& scvfGeometry : boundaryScvfGeometries)
for (unsigned int isScvfLocalIdx = 0; isScvfLocalIdx < isGeometry.corners(); ++isScvfLocalIdx)
{
// find the scvs this scvf is belonging to
std::vector<IndexType> localScvIndices =
{static_cast<IndexType>(referenceElement.subEntity(intersection.indexInInside(), 1,
isScvfLocalIdx++, dim))};
isScvfLocalIdx, dim))};
// get the corner points and the area of the boundary scv face
auto scvfCorners = geometryHelper.getBoundaryScvfCorners(isGeometry, isScvfLocalIdx);
auto area = geometryHelper.area(scvfCorners);
scvfs_[eIdx].emplace_back(std::move(scvfGeometry),
scvfs_[eIdx].emplace_back(std::move(scvfCorners),
intersection.centerUnitOuterNormal(),
area,
scvfLocalIdx,
localScvIndices,
true);
......
......@@ -25,49 +25,99 @@
#include <dune/common/fvector.hh>
#include <dumux/discretization/subcontrolvolumebase.hh>
#include <dumux/common/math.hh>
namespace Dumux
{
template<class G, typename I>
class BoxSubControlVolume : public SubControlVolumeBase<G, I>
class BoxSubControlVolume : public SubControlVolumeBase<BoxSubControlVolume<G, I>, G, I>
{
public:
// exported types
using ParentType = SubControlVolumeBase<BoxSubControlVolume<G, I>, G, I>;
using Geometry = G;
using IndexType = I;
private:
using Scalar = typename Geometry::ctype;
enum { dim = Geometry::mydimension };
enum { dimworld = Geometry::coorddimension };
using GlobalPosition = Dune::FieldVector<Scalar, dimworld>;
public:
// the contructor in the box case
BoxSubControlVolume(Geometry&& geometry,
BoxSubControlVolume(std::vector<GlobalPosition>&& corners,
Scalar volume,
IndexType scvIdx,
IndexType elementIndex,
IndexType dofIndex)
: SubControlVolumeBase<G, I>(std::move(geometry), std::move(elementIndex)),
: ParentType(),
corners_(std::move(corners)),
center_(0.0),
volume_(volume),
elementIndex_(std::move(elementIndex)),
scvIdx_(std::move(scvIdx)),
dofIndex_(std::move(dofIndex)) {}
dofIndex_(std::move(dofIndex))
{
// compute center point
for (const auto& corner : corners_)
center_ += corner;
center_ /= corners_.size();
}
//! The center of the sub control volume
GlobalPosition center() const
{
return center_;
}
//! The local index of this scv
//! The volume of the sub control volume
Scalar volume() const
{
return volume_;
}
//! The geometry of the sub control volume
// e.g. for integration
Geometry geometry() const
{
return Geometry(Dune::GeometryType(Dune::GeometryType::cube, dim), corners_);
}
//! The global index of this scv
IndexType index() const
{
return scvIdx_;
}
//! The index of the dof this scv is embedded in
IndexType dofIndex() const
{
return dofIndex_;
}
// The position of the dof this scv is embedded in
GlobalPosition dofPosition() const
{
return this->geometry().corner(scvIdx_);
// The corner list is defined such that the first entry is the vertex itself
return corner(0);
}
//! The global index of the element this scv is embedded in
IndexType elementIndex() const
{
return elementIndex_;
}
//! Return the corner for the given local index
GlobalPosition corner(unsigned int localIdx) const
{
assert(localIdx < corners_.size() && "provided index exceeds the number of corners");
return corners_[localIdx];
}
private:
std::vector<GlobalPosition> corners_;
GlobalPosition center_;
Scalar volume_;
IndexType elementIndex_;
IndexType scvIdx_;
IndexType dofIndex_;
};
......
......@@ -35,9 +35,13 @@ namespace Dumux
* \brief Class for a sub control volume face in the box method, i.e a part of the boundary
* of a sub control volume we compute fluxes on. We simply use the base class here.
*/
template<class Geometry, typename IndexType>
class BoxSubControlVolumeFace : public SubControlVolumeFaceBase<Geometry, IndexType>
template<class G, typename I>
class BoxSubControlVolumeFace : public SubControlVolumeFaceBase<BoxSubControlVolumeFace<G, I>, G, I>
{
using ParentType = SubControlVolumeFaceBase<BoxSubControlVolumeFace<G, I>, G, I>;
using Geometry = G;
using IndexType = I;
using Scalar = typename Geometry::ctype;
static const int dim = Geometry::mydimension;
static const int dimworld = Geometry::coorddimension;
......@@ -46,12 +50,96 @@ class BoxSubControlVolumeFace : public SubControlVolumeFaceBase<Geometry, IndexT
using LocalPosition = Dune::FieldVector<Scalar, dim>;
public:
BoxSubControlVolumeFace(Geometry&& geometry,
BoxSubControlVolumeFace(std::vector<GlobalPosition>&& corners,
const GlobalPosition& unitOuterNormal,
Scalar area,
IndexType scvfIndex,
const std::vector<IndexType>& scvIndices,
bool boundary = false)
: SubControlVolumeFaceBase<Geometry, IndexType>(std::move(geometry), geometry.center(), unitOuterNormal, scvfIndex, scvIndices, boundary) {}
: ParentType(),
corners_(std::move(corners)),