Merge branch 'fix/epsilon-2d-2d-intersections' into 'master'

Fix/epsilon 2d 2d intersections

See merge request !2918

dumux/geometry/geometryintersection.hh

@@ -877,14 +877,13 @@ public:
         const auto b2 = geo2.corner(2) - geo2.corner(0);
         const auto n2 = crossProduct(a2, b2);
 
-        using std::max;
-        using std::sqrt;
+        // compute an epsilon scaled with larger edge length
         const auto a1Norm2 = a1.two_norm2();
         const auto b1Norm2 = b1.two_norm2();
+
+        using std::max;
         const auto maxNorm2 = max(a1Norm2, b1Norm2);
         const auto eps2 = maxNorm2*eps_;
-        const auto eps = sqrt(maxNorm2)*eps_;
-        const auto eps3 = eps*eps2;
 
         // early exit if polygons are not parallel
         using std::abs;
@@ -896,6 +895,10 @@ public:
         if (d1.two_norm2() < eps2)
             d1 = geo2.corner(1) - geo1.corner(0);
 
+        using std::sqrt;
+        const auto maxNorm = sqrt(maxNorm2);
+        const auto eps3 = maxNorm*eps2;
+
         if (abs(d1*n2) > eps3)
             return false;
 
@@ -960,6 +963,7 @@ public:
             return false;
 
         // remove duplicates
+        const auto eps = maxNorm*eps_;
         std::sort(points.begin(), points.end(), [eps] (const auto& a, const auto& b) -> bool
         {
             using std::abs;

test/geometry/test_2d2d_intersection.cc

@@ -172,6 +172,28 @@ void testNonParallelPolygons(std::vector<bool>& returns)
     std::cout << std::endl;
 }
 
+// unit test discovered in a simulation, that originally failed on commit 51fc9c1e3
+void testSimulationCase(std::vector<bool>& returns)
+{
+    using Position = Dune::FieldVector<double, 3>;
+    const auto quad1 = makeQuadrilateral(
+        Position{{9.49264705847428835739e-01, 1.40073529420942810564e-01, 4.00000000000000077716e-01}},
+        Position{{9.44669117627101373458e-01, 1.42095588240539000280e-01, 4.00000000000000077716e-01}},
+        Position{{9.49632352923713174420e-01, 1.45036764710471721695e-01, 4.00000000000000022204e-01}},
+        Position{{9.46446078418066827354e-01, 1.44730392160359544462e-01, 4.00000000000000077716e-01}}
+    );
+    const auto quad2 = makeQuadrilateral(
+        Position{{9.39999999999999946709e-01, 1.39999999999999985567e-01, 4.00000000000000022204e-01}},
+        Position{{9.59999999999999964473e-01, 1.39999999999999985567e-01, 4.00000000000000022204e-01}},
+        Position{{9.39999999999999946709e-01, 1.60000000000000003331e-01, 4.00000000000000022204e-01}},
+        Position{{9.59999999999999964473e-01, 1.60000000000000003331e-01, 4.00000000000000022204e-01}}
+    );
+
+    returns.push_back(testPolygonIntersection<3>(quad1, quad2, true));
+
+    std::cout << std::endl;
+}
+
 #endif
 
 int main(int argc, char* argv[])
@@ -190,6 +212,9 @@ int main(int argc, char* argv[])
     std::cout << "Testing non-parallel polygons in 3d space" << std::endl;
     testNonParallelPolygons(returns);
 
+    std::cout << "Testing 3d unit test from actual simulation" << std::endl;
+    testSimulationCase(returns);
+
     // TODO: implement and test point and segment intersections
 
     // determine the exit code