Merge branch 'fix/point-in-geometry-2d-simplex' into 'master'

Fix/point in geometry 2d simplex

Closes #479

See merge request !1141

dumux/common/geometry/intersectspointsimplex.hh

@@ -82,18 +82,20 @@ bool intersectsPointSimplex(const Dune::FieldVector<ctype, dimworld>& point,
 {
     // adapted from the algorithm from from "Real-Time Collision Detection" by Christer Ericson,
     // published by Morgan Kaufmann Publishers, (c) 2005 Elsevier Inc. (Chapter 5.4.2)
-    static constexpr ctype eps_ = 1.0e-7;
+    constexpr ctype eps_ = 1.0e-7;
 
     // compute the normal of the triangle
     const auto v1 = p0 - p2;
     auto n = crossProduct(v1, p1 - p0);
-    n /= n.two_norm();
+    const ctype nnorm = n.two_norm();
+    const ctype eps4 = eps_*nnorm*nnorm; // compute an epsilon for later
+    n /= nnorm; // normalize
 
     // first check if we are in the plane of the triangle
     // if not we can return early
     using std::abs;
     auto x = p0 - point;
-    x /= x.two_norm();
+    x /= x.two_norm(); // normalize
 
     if (abs(x*n) > eps_)
         return false;
@@ -107,14 +109,16 @@ bool intersectsPointSimplex(const Dune::FieldVector<ctype, dimworld>& point,
     const auto u = crossProduct(b, c);
     const auto v = crossProduct(c, a);
 
-    // they have to point in the same direction
-    if (u*v < 0.0) return false;
+    // they have to point in the same direction or be orthogonal
+    if (u*v < 0.0 - eps4)
+        return false;
 
     // compute the normal vector for triangle P->C->A
     const auto w = crossProduct(a, b);
 
-    // it also has to point in the same direction
-    if (u*w < 0.0) return false;
+    // it also has to point in the same direction or be orthogonal
+    if (u*w < 0.0 - eps4)
+        return false;
 
     // now the point must be in the triangle (or on the faces)
     return true;

dumux/common/geometry/makegeometry.hh

@@ -24,6 +24,7 @@
 
 #include <vector>
 #include <array>
+#include <limits>
 #include <dune/common/fvector.hh>
 #include <dune/common/fmatrix.hh>
 #include <dune/geometry/multilineargeometry.hh>
@@ -34,17 +35,39 @@ namespace Dumux {
 
 //! Checks if four points lie within the same plane.
 template<class CoordScalar>
-bool pointsAreCoplanar(const std::vector<Dune::FieldVector<CoordScalar, 3>>& points, CoordScalar eps = 1e-20)
+bool pointsAreCoplanar(const std::vector<Dune::FieldVector<CoordScalar, 3>>& points, const CoordScalar scale)
 {
     assert(points.size() == 4);
     // (see "Real-Time Collision Detection" by Christer Ericson)
     Dune::FieldMatrix<CoordScalar, 4, 4> M;
     for(int i = 0; i < 3; ++i )
         M[i] = {points[0][i], points[1][i], points[2][i], points[3][i]};
-    M[3] = {1.0, 1.0, 1.0, 1.0};
+    M[3] = {1.0*scale, 1.0*scale, 1.0*scale, 1.0*scale};
 
     using std::abs;
-    return abs(M.determinant()) < eps;
+    return abs(M.determinant()) < 1.5e-7*scale*scale*scale*scale;
+}
+
+//! Checks if four points lie within the same plane.
+template<class CoordScalar>
+bool pointsAreCoplanar(const std::vector<Dune::FieldVector<CoordScalar, 3>>& points)
+{
+    Dune::FieldVector<CoordScalar, 3> bBoxMin(std::numeric_limits<CoordScalar>::max());
+    Dune::FieldVector<CoordScalar, 3> bBoxMax(std::numeric_limits<CoordScalar>::lowest());
+    for (const auto& p : points)
+    {
+        for (int i=0; i<3; i++)
+        {
+            using std::min;
+            using std::max;
+            bBoxMin[i] = min(bBoxMin[i], p[i]);
+            bBoxMax[i] = max(bBoxMax[i], p[i]);
+        }
+    }
+
+    const auto size = (bBoxMax - bBoxMin).two_norm();
+
+    return pointsAreCoplanar(points, size);
 }
 
 /*!
@@ -127,8 +150,24 @@ auto makeDuneQuadrilaterial(const std::vector<Dune::FieldVector<CoordScalar, 3>>
     if(!enableSanityCheck)
         return GeometryType(Dune::GeometryTypes::quadrilateral, points);
 
+    // compute size
+    Dune::FieldVector<CoordScalar, 3> bBoxMin(std::numeric_limits<CoordScalar>::max());
+    Dune::FieldVector<CoordScalar, 3> bBoxMax(std::numeric_limits<CoordScalar>::lowest());
+    for (const auto& p : points)
+    {
+        for (int i=0; i<3; i++)
+        {
+            using std::min;
+            using std::max;
+            bBoxMin[i] = min(bBoxMin[i], p[i]);
+            bBoxMax[i] = max(bBoxMax[i], p[i]);
+        }
+    }
+
+    const auto size = (bBoxMax - bBoxMin).two_norm();
+
     // otherwise, perform a number of checks and corrections
-    if(!pointsAreCoplanar(points))
+    if(!pointsAreCoplanar(points, size))
         DUNE_THROW(Dune::InvalidStateException, "Points do not lie within a plane");
 
     // make sure points conform with dune ordering
@@ -140,8 +179,8 @@ auto makeDuneQuadrilaterial(const std::vector<Dune::FieldVector<CoordScalar, 3>>
 
     const auto quadrilateral = GeometryType(Dune::GeometryTypes::quadrilateral, corners);
 
-    const auto eps = 1e-20;
-    if(quadrilateral.volume() < eps)
+    const auto eps = 1e-7;
+    if(quadrilateral.volume() < eps*size*size)
         DUNE_THROW(Dune::InvalidStateException, "Something went wrong, geometry has area of zero");
 
     return quadrilateral;

test/common/geometry/test_makegeometry.cc

@@ -79,7 +79,7 @@ void permutatePointsAndTest(const std::vector<GlobalPosition>& cornerPoints,
                         std::cout << "point " << p << " lies within the quadrilateral" << std::endl;
                 }
                 else
-                    DUNE_THROW(Dune::InvalidStateException, "Check for point inside geometry failed. Point " << p << " does not lie within the geometry!");
+                    DUNE_THROW(Dune::InvalidStateException, "False negative: Check for point " << p << " which is inside the geometry failed");
             }
 
             for(const auto& p : pointsOutsideGeometry)
@@ -90,7 +90,7 @@ void permutatePointsAndTest(const std::vector<GlobalPosition>& cornerPoints,
                         std::cout << "point " << p << " lies outside of the quadrilateral" << std::endl;
                 }
                 else
-                    DUNE_THROW(Dune::InvalidStateException, "Check for point outside geometry failed. Point " << p << " does lie within the geometry!");
+                    DUNE_THROW(Dune::InvalidStateException, "False positive: Check for point " << p << " which is outside the geometry failed");
             }
 
     } while(std::next_permutation(s.begin(), s.end()));
@@ -100,20 +100,21 @@ template<class GlobalPosition>
 void checkAxisAlignedGeometry(std::vector<GlobalPosition>& cornerPoints,
                               std::vector<GlobalPosition>& pointsWithinGeometry,
                               std::vector<GlobalPosition>& pointsOutsideGeometry,
-                              const int normalDirection)
+                              const int normalDirection,
+                              const typename GlobalPosition::value_type scale)
 {
     static const char dim[]= "xyz";
     std::cout << "testing for quadrilateral with normal in " << dim[normalDirection] << " direction" << std::endl;
 
     // check if points are coplanar
     if(!Dumux::pointsAreCoplanar(cornerPoints))
-        DUNE_THROW(Dune::InvalidStateException, "False positive, points are actually coplanar!");
+        DUNE_THROW(Dune::InvalidStateException, "False negative, points are actually coplanar!");
 
-    cornerPoints[0][normalDirection] += 1e-9;
+    cornerPoints[0][normalDirection] += 1e-5*scale; // we make them non-coplanar
     if(Dumux::pointsAreCoplanar(cornerPoints))
-        DUNE_THROW(Dune::InvalidStateException, "Points are not coplanar!");
+        DUNE_THROW(Dune::InvalidStateException, "False positive, points are actually not coplanar!");
 
-    cornerPoints[0][normalDirection] -= 1e-9;
+    cornerPoints[0][normalDirection] -= 1e-5*scale;
 
     permutatePointsAndTest(cornerPoints, pointsWithinGeometry, pointsOutsideGeometry);
 
@@ -122,19 +123,22 @@ void checkAxisAlignedGeometry(std::vector<GlobalPosition>& cornerPoints,
     std::vector<GlobalPosition> pointsWithinGeometryForRandomTest = { {0.6, 0.6, 0.6},
                                                                      {0.4, 0.4, 0.4} };
     for(auto& p : pointsWithinGeometryForRandomTest)
+    {
+        p *= scale;
         p[normalDirection] = 0.0;
+    }
 
     auto pointsOutsideGeometryForRandomTest = pointsOutsideGeometry;
-    pointsOutsideGeometryForRandomTest[0] = {1.5, 1.5, 1.5};
+    pointsOutsideGeometryForRandomTest[0] = {1.5*scale, 1.5*scale, 1.5*scale};
     pointsOutsideGeometryForRandomTest[0][normalDirection] = 0.0;
 
-    for(int i = 0; i< 10; i++)
-    {
-        // uniform random number generator
-        std::random_device rd;
-        std::mt19937 generator(rd());
-        std::uniform_real_distribution<> uniformdist(-0.3, 0.3);
+    // uniform random number generator
+    std::random_device rd;
+    std::mt19937 generator(rd());
+    std::uniform_real_distribution<> uniformdist(-0.3*scale, 0.3*scale);
 
+    for(int i = 0; i < 10; i++)
+    {
         for(auto&& p : cornerPoints)
         {
             for(int x = 0; x < p.size(); ++x)
@@ -148,8 +152,6 @@ void checkAxisAlignedGeometry(std::vector<GlobalPosition>& cornerPoints,
 
         cornerPoints = origCornerPoints;
     }
-
-
 }
 
 int main(int argc, char** argv) try
@@ -157,50 +159,59 @@ int main(int argc, char** argv) try
     using namespace Dumux;
     using GlobalPosition = Dune::FieldVector<double, 3>;
 
-    GlobalPosition p0 = {0,0,0};
-    GlobalPosition p1 = {1,0,0};
-    GlobalPosition p2 = {0,1,0};
-    GlobalPosition p3 = {1,1,0};
+    std::array<double, 3> scaling{{1e-12, 1.0, 1e12}};
 
-    std::vector<GlobalPosition> cornerPoints = {p0, p1, p2, p3};
+    for (const double scale : scaling)
+    {
+        const double size = 1.0*scale;
+        const double half = 0.5*scale;
+        const double small = 1e-3*scale;
 
-    std::vector<GlobalPosition> pointsWithinGeometry = { GlobalPosition{0.5, 0.5, 0.0},
-                                                         GlobalPosition{cornerPoints[0][0] + 1e-3, cornerPoints[0][1] + 1e-3, 0.0},
-                                                         GlobalPosition{cornerPoints[3][0] - 1e-3, cornerPoints[3][1] - 1e-3, 0.0} };
+        GlobalPosition p0 = {0, 0, 0};
+        GlobalPosition p1 = {size, 0, 0};
+        GlobalPosition p2 = {0, size, 0};
+        GlobalPosition p3 = {size, size, 0};
 
-    std::vector<GlobalPosition> pointsOutsideGeometry = { GlobalPosition{cornerPoints[0][0] - 1e-3, cornerPoints[0][1] - 1e-3, 0.0},
-                                                          GlobalPosition{0.5, 0.5, 1e-3} };
+        std::vector<GlobalPosition> cornerPoints = {p0, p1, p2, p3};
 
-    // do the checks for a quadrilateral parallel to the x and y axis
-    checkAxisAlignedGeometry(cornerPoints, pointsWithinGeometry, pointsOutsideGeometry, 2);
+        std::vector<GlobalPosition> pointsWithinGeometry = { GlobalPosition{half, half, 0.0},
+                                                             GlobalPosition{cornerPoints[0][0] + small, cornerPoints[0][1] + small, 0.0},
+                                                             GlobalPosition{cornerPoints[3][0] - small, cornerPoints[3][1] - small, 0.0} };
 
-    // rotate the quadrilateral to make it parallel to the other axes and test again
-    for(int i = 1; i >=0; --i)
-    {
-        for(auto& p : cornerPoints)
-            std::rotate(p.begin(), p.begin() + 1, p.end());
-        for(auto& p : pointsWithinGeometry)
-            std::rotate(p.begin(), p.begin() + 1, p.end());
-        for(auto& p : pointsOutsideGeometry)
-            std::rotate(p.begin(), p.begin() + 1, p.end());
-
-        checkAxisAlignedGeometry(cornerPoints, pointsWithinGeometry, pointsOutsideGeometry, i);
-    }
+        std::vector<GlobalPosition> pointsOutsideGeometry = { GlobalPosition{cornerPoints[0][0] - small, cornerPoints[0][1] - small, 0.0},
+                                                              GlobalPosition{half, half, small} };
 
-    std::cout << "testing for non axis-aligned quadrilateral" << std::endl;
+        // do the checks for a quadrilateral parallel to the x and y axis
+        checkAxisAlignedGeometry(cornerPoints, pointsWithinGeometry, pointsOutsideGeometry, 2, size);
 
-    cornerPoints[0][0] += 0.5;
-    cornerPoints[1][0] -= 0.5;
-    cornerPoints[2][0] += 0.5;
-    cornerPoints[3][0] -= 0.5;
+        // rotate the quadrilateral to make it parallel to the other axes and test again
+        for(int i = 1; i >=0; --i)
+        {
+            for(auto& p : cornerPoints)
+                std::rotate(p.begin(), p.begin() + 1, p.end());
+            for(auto& p : pointsWithinGeometry)
+                std::rotate(p.begin(), p.begin() + 1, p.end());
+            for(auto& p : pointsOutsideGeometry)
+                std::rotate(p.begin(), p.begin() + 1, p.end());
+
+            checkAxisAlignedGeometry(cornerPoints, pointsWithinGeometry, pointsOutsideGeometry, i, size);
+        }
 
-    GlobalPosition pointToCheck5 = {0.0, 0.5, 0.5};
+        std::cout << "testing for non axis-aligned quadrilateral" << std::endl;
 
-    pointsWithinGeometry = {pointToCheck5};
+        cornerPoints[0][0] += half;
+        cornerPoints[1][0] -= half;
+        cornerPoints[2][0] += half;
+        cornerPoints[3][0] -= half;
 
-    pointsOutsideGeometry[1] = pointsOutsideGeometry[0];
+        GlobalPosition pointToCheck5 = {0.0, half, half};
 
-    permutatePointsAndTest(cornerPoints, pointsWithinGeometry, pointsOutsideGeometry);
+        pointsWithinGeometry = {pointToCheck5};
+
+        pointsOutsideGeometry[1] = pointsOutsideGeometry[0];
+
+        permutatePointsAndTest(cornerPoints, pointsWithinGeometry, pointsOutsideGeometry);
+    }
 
     return 0;
 }