diff --git a/dumux/discretization/box/boxgeometryhelper.hh b/dumux/discretization/box/boxgeometryhelper.hh
index 61877fbe936ceeef5c1162740a140266ea286797..ce5e31cac9912b807f17af8cef7c055d57c5f7b2 100644
--- a/dumux/discretization/box/boxgeometryhelper.hh
+++ b/dumux/discretization/box/boxgeometryhelper.hh
@@ -268,7 +268,7 @@ public:
GlobalPosition normal = Dumux::crossProduct(v3, t);
normal /= normal.two_norm();
- // TODO can this be done easier?, e.g. always ensure the right direction?
+ //! ensure the right direction of the normal
const auto v = elementGeometry_.corner(scvIndices[1]) - elementGeometry_.corner(scvIndices[0]);
const auto s = v*normal;
if (std::signbit(s))
@@ -283,9 +283,17 @@ public:
normal(const ScvfCornerStorage& scvfCorners,
const std::vector<unsigned int>& scvIndices) const
{
+ //! obtain normal vector by 90° counter-clockwise rotation of t
const auto t = scvfCorners[1] - scvfCorners[0];
GlobalPosition normal({-t[1], t[0]});
normal /= normal.two_norm();
+
+ //! ensure the right direction of the normal
+ const auto v = elementGeometry_.corner(scvIndices[1]) - elementGeometry_.corner(scvIndices[0]);
+ const auto s = v*normal;
+ if (std::signbit(s))
+ normal *= -1;
+
return normal;
}
@@ -302,7 +310,10 @@ public:
typename std::enable_if<w == 2, Scalar>::type
scvVolume(const ScvCornerStorage& p) const
{
- return 0.5*Dumux::crossProduct(p[3]-p[0], p[2]-p[1]);
+ //! make sure we are using positive volumes
+ //! Cross product of diagonals might be negative, depending on element orientation
+ using std::abs;
+ return 0.5*abs(Dumux::crossProduct(p[3]-p[0], p[2]-p[1]));
}
//! get scvf area
diff --git a/dumux/discretization/box/darcyslaw.hh b/dumux/discretization/box/darcyslaw.hh
index 4d8eadc42fba11995ce80d3d96e8887b2e350c27..e39e95ad0efe4616af3a9b5e8e775fc1a6e30ea5 100644
--- a/dumux/discretization/box/darcyslaw.hh
+++ b/dumux/discretization/box/darcyslaw.hh
@@ -111,8 +111,7 @@ public:
gradP.axpy(-rho, problem.gravityAtPos(scvf.center()));
// apply the permeability and return the flux
- const auto KGradP = applyPermeability_(K, gradP);
- return -1.0*(KGradP*scvf.unitOuterNormal())*scvf.area();
+ return -1.0*vtmv(scvf.unitOuterNormal(), K, gradP)*scvf.area();
}
// compute transmissibilities ti for analytical jacobians
@@ -140,26 +139,10 @@ public:
std::vector<Scalar> ti(fvGeometry.numScv());
for (const auto& scv : scvs(fvGeometry))
ti[scv.indexInElement()] =
- -1.0*(applyPermeability_(K, fluxVarCache.gradN(scv.indexInElement()))
- *scvf.unitOuterNormal())*scvf.area();
+ -1.0*scvf.area()*vtmv(scvf.unitOuterNormal(), K, fluxVarCache.gradN(scv.indexInElement()));
return ti;
}
-
-private:
- inline static GlobalPosition applyPermeability_(const DimWorldMatrix& K, const GlobalPosition& gradI)
- {
- GlobalPosition result(0.0);
- K.mv(gradI, result);
- return result;
- }
-
- inline static GlobalPosition applyPermeability_(const Scalar k, const GlobalPosition& gradI)
- {
- GlobalPosition result(gradI);
- result *= k;
- return result;
- }
};
} // end namespace Dumux
diff --git a/dumux/discretization/box/fvelementgeometry.hh b/dumux/discretization/box/fvelementgeometry.hh
index e27f68ed74dd81578ecf240b2f0173b4c535351a..e67bc26c22c4285d1a63d87d2d6b46830e8ef545 100644
--- a/dumux/discretization/box/fvelementgeometry.hh
+++ b/dumux/discretization/box/fvelementgeometry.hh
@@ -170,6 +170,7 @@ class BoxFVElementGeometry<TypeTag, false>
using IndexType = typename GridView::IndexSet::IndexType;
using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
+ using LocalIndexType = typename SubControlVolumeFace::Traits::LocalIndexType;
using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
@@ -283,10 +284,10 @@ private:
// construct the sub control volumes
scvs_.resize(elementGeometry.corners());
- for (unsigned int scvLocalIdx = 0; scvLocalIdx < elementGeometry.corners(); ++scvLocalIdx)
+ for (LocalIndexType scvLocalIdx = 0; scvLocalIdx < elementGeometry.corners(); ++scvLocalIdx)
{
// get asssociated dof index
- auto dofIdxGlobal = fvGridGeometry().vertexMapper().subIndex(element, scvLocalIdx, dim);
+ const auto dofIdxGlobal = fvGridGeometry().vertexMapper().subIndex(element, scvLocalIdx, dim);
// add scv to the local container
scvs_[scvLocalIdx] = SubControlVolume(geometryHelper,
@@ -303,14 +304,14 @@ private:
for (; scvfLocalIdx < numInnerScvf; ++scvfLocalIdx)
{
// find the local scv indices this scvf is connected 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))});
+ std::vector<LocalIndexType> localScvIndices({static_cast<LocalIndexType>(referenceElement.subEntity(scvfLocalIdx, dim-1, 0, dim)),
+ static_cast<LocalIndexType>(referenceElement.subEntity(scvfLocalIdx, dim-1, 1, dim))});
scvfs_[scvfLocalIdx] = SubControlVolumeFace(geometryHelper,
element,
elementGeometry,
scvfLocalIdx,
- localScvIndices,
+ std::move(localScvIndices),
false);
}
@@ -324,15 +325,15 @@ private:
for (unsigned int isScvfLocalIdx = 0; isScvfLocalIdx < isGeometry.corners(); ++isScvfLocalIdx)
{
// find the scv this scvf is connected to
- IndexType insideScvIdx = static_cast<IndexType>(referenceElement.subEntity(intersection.indexInInside(), 1, isScvfLocalIdx, dim));
- std::vector<IndexType> localScvIndices = {insideScvIdx, insideScvIdx};
+ const LocalIndexType insideScvIdx = static_cast<LocalIndexType>(referenceElement.subEntity(intersection.indexInInside(), 1, isScvfLocalIdx, dim));
+ std::vector<LocalIndexType> localScvIndices = {insideScvIdx, insideScvIdx};
scvfs_.emplace_back(geometryHelper,
intersection,
isGeometry,
isScvfLocalIdx,
scvfLocalIdx,
- localScvIndices,
+ std::move(localScvIndices),
true);
// increment local counter
diff --git a/dumux/discretization/box/fvgridgeometry.hh b/dumux/discretization/box/fvgridgeometry.hh
index eddffd98d797187d29060406595b993f422d3d2f..76801c88b90b41f35d5ef5b9a48a8bd8d1afbb21 100644
--- a/dumux/discretization/box/fvgridgeometry.hh
+++ b/dumux/discretization/box/fvgridgeometry.hh
@@ -54,6 +54,7 @@ class BoxFVGridGeometry<TypeTag, true> : public BaseFVGridGeometry<TypeTag>
using IndexType = typename GridView::IndexSet::IndexType;
using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
+ using LocalIndexType = typename SubControlVolumeFace::Traits::LocalIndexType;
using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
using VertexMapper = typename GET_PROP_TYPE(TypeTag, VertexMapper);
using Element = typename GridView::template Codim<0>::Entity;
@@ -132,9 +133,9 @@ public:
// construct the sub control volumes
scvs_[eIdx].resize(elementGeometry.corners());
- for (unsigned int scvLocalIdx = 0; scvLocalIdx < elementGeometry.corners(); ++scvLocalIdx)
+ for (LocalIndexType scvLocalIdx = 0; scvLocalIdx < elementGeometry.corners(); ++scvLocalIdx)
{
- auto dofIdxGlobal = this->vertexMapper().subIndex(element, scvLocalIdx, dim);
+ const auto dofIdxGlobal = this->vertexMapper().subIndex(element, scvLocalIdx, dim);
scvs_[eIdx][scvLocalIdx] = SubControlVolume(geometryHelper,
scvLocalIdx,
@@ -143,19 +144,19 @@ public:
}
// construct the sub control volume faces
- unsigned int scvfLocalIdx = 0;
+ LocalIndexType scvfLocalIdx = 0;
scvfs_[eIdx].resize(element.subEntities(dim-1));
for (; scvfLocalIdx < element.subEntities(dim-1); ++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))});
+ std::vector<LocalIndexType> localScvIndices({static_cast<LocalIndexType>(referenceElement.subEntity(scvfLocalIdx, dim-1, 0, dim)),
+ static_cast<LocalIndexType>(referenceElement.subEntity(scvfLocalIdx, dim-1, 1, dim))});
scvfs_[eIdx][scvfLocalIdx] = SubControlVolumeFace(geometryHelper,
element,
elementGeometry,
scvfLocalIdx,
- localScvIndices,
+ std::move(localScvIndices),
false);
}
@@ -172,15 +173,15 @@ public:
for (unsigned int isScvfLocalIdx = 0; isScvfLocalIdx < isGeometry.corners(); ++isScvfLocalIdx)
{
// find the scvs this scvf is belonging to
- IndexType insideScvIdx = static_cast<IndexType>(referenceElement.subEntity(intersection.indexInInside(), 1, isScvfLocalIdx, dim));
- std::vector<IndexType> localScvIndices = {insideScvIdx, insideScvIdx};
+ const LocalIndexType insideScvIdx = static_cast<LocalIndexType>(referenceElement.subEntity(intersection.indexInInside(), 1, isScvfLocalIdx, dim));
+ std::vector<LocalIndexType> localScvIndices = {insideScvIdx, insideScvIdx};
scvfs_[eIdx].emplace_back(geometryHelper,
intersection,
isGeometry,
isScvfLocalIdx,
scvfLocalIdx,
- localScvIndices,
+ std::move(localScvIndices),
true);
// increment local counter
diff --git a/dumux/discretization/box/subcontrolvolumeface.hh b/dumux/discretization/box/subcontrolvolumeface.hh
index f69e48f62e2a6385273c61cfcc546631b8c6ff21..ab3eb636136417e482f4986f7e24ff7a8018cea0 100644
--- a/dumux/discretization/box/subcontrolvolumeface.hh
+++ b/dumux/discretization/box/subcontrolvolumeface.hh
@@ -58,14 +58,14 @@ public:
const Element& element,
const typename Element::Geometry& elemGeometry,
GridIndexType scvfIndex,
- const std::vector<LocalIndexType>& scvIndices,
+ std::vector<LocalIndexType>&& scvIndices,
bool boundary = false)
: corners_(geometryHelper.getScvfCorners(scvfIndex)),
center_(0.0),
unitOuterNormal_(geometryHelper.normal(corners_, scvIndices)),
area_(geometryHelper.scvfArea(corners_)),
scvfIndex_(scvfIndex),
- scvIndices_(scvIndices),
+ scvIndices_(std::move(scvIndices)),
boundary_(boundary)
{
for (const auto& corner : corners_)
@@ -80,14 +80,14 @@ public:
const typename Intersection::Geometry& isGeometry,
LocalIndexType indexInIntersection,
GridIndexType scvfIndex,
- const std::vector<LocalIndexType>& scvIndices,
+ std::vector<LocalIndexType>&& scvIndices,
bool boundary = false)
: corners_(geometryHelper.getBoundaryScvfCorners(isGeometry, indexInIntersection)),
center_(0.0),
unitOuterNormal_(intersection.centerUnitOuterNormal()),
area_(geometryHelper.scvfArea(corners_)),
scvfIndex_(scvfIndex),
- scvIndices_(scvIndices),
+ scvIndices_(std:: move(scvIndices)),
boundary_(boundary)
{
for (const auto& corner : corners_)