Skip to content
Snippets Groups Projects
Commit 1dbd052b authored by Dennis Gläser's avatar Dennis Gläser
Browse files

[box][geomhelper] use correct normals and positive volumes

parent 38418d52
No related branches found
No related tags found
Loading
This commit is part of merge request !639. Comments created here will be created in the context of that merge request.
...@@ -268,7 +268,7 @@ public: ...@@ -268,7 +268,7 @@ public:
GlobalPosition normal = Dumux::crossProduct(v3, t); GlobalPosition normal = Dumux::crossProduct(v3, t);
normal /= normal.two_norm(); 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 v = elementGeometry_.corner(scvIndices[1]) - elementGeometry_.corner(scvIndices[0]);
const auto s = v*normal; const auto s = v*normal;
if (std::signbit(s)) if (std::signbit(s))
...@@ -283,9 +283,17 @@ public: ...@@ -283,9 +283,17 @@ public:
normal(const ScvfCornerStorage& scvfCorners, normal(const ScvfCornerStorage& scvfCorners,
const std::vector<unsigned int>& scvIndices) const const std::vector<unsigned int>& scvIndices) const
{ {
//! obtain normal vector by 90° counter-clockwise rotation of t
const auto t = scvfCorners[1] - scvfCorners[0]; const auto t = scvfCorners[1] - scvfCorners[0];
GlobalPosition normal({-t[1], t[0]}); GlobalPosition normal({-t[1], t[0]});
normal /= normal.two_norm(); 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; return normal;
} }
...@@ -302,7 +310,10 @@ public: ...@@ -302,7 +310,10 @@ public:
typename std::enable_if<w == 2, Scalar>::type typename std::enable_if<w == 2, Scalar>::type
scvVolume(const ScvCornerStorage& p) const 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 //! get scvf area
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment