From e89389cd51630856d31420f9e8bf2c7b90c464f8 Mon Sep 17 00:00:00 2001
From: "Dennis.Glaeser" <dennis.glaeser@iws.uni-stuttgart.de>
Date: Fri, 17 Apr 2020 14:45:06 +0200
Subject: [PATCH] [geoisection] improve 3d-3d segment algorithm
---
dumux/common/geometry/geometryintersection.hh | 23 +++++++++----------
1 file changed, 11 insertions(+), 12 deletions(-)
diff --git a/dumux/common/geometry/geometryintersection.hh b/dumux/common/geometry/geometryintersection.hh
index 479c00534d..8ae5add778 100644
--- a/dumux/common/geometry/geometryintersection.hh
+++ b/dumux/common/geometry/geometryintersection.hh
@@ -1344,28 +1344,28 @@ public:
const auto& a = geo1.corner(0);
const auto& b = geo1.corner(1);
+ const auto ab = b-a;
+
const auto& p = geo2.corner(0);
const auto& q = geo2.corner(1);
-
- const auto ab = b-a;
const auto pq = q-p;
- const auto abNorm = ab.two_norm();
- const auto pqNorm = pq.two_norm();
+
+ const auto abNorm2 = ab.two_norm2();
+ const auto pqNorm2 = pq.two_norm2();
using std::max;
- const auto maxNorm = max(abNorm, pqNorm);
- const auto eps2 = eps_*maxNorm*maxNorm;
+ const auto eps2 = eps_*max(abNorm2, pqNorm2);
// if the segment are not parallel there is no segment intersection
using std::abs;
- if (!(abs(abs(ab*pq) - abNorm*pqNorm) < eps2))
+ if (crossProduct(ab, pq).two_norm2() > eps2*eps2)
return false;
const auto ap = (p-a);
const auto aq = (q-a);
// points have to be colinear
- if (!(crossProduct(ap, aq).two_norm() < eps2))
+ if (crossProduct(ap, aq).two_norm2() > eps2*eps2)
return false;
// scaled local coordinates
@@ -1373,7 +1373,6 @@ public:
// and do it in the very end if we find an intersection
auto tp = ap*ab;
auto tq = aq*ab;
- const auto abnorm2 = abNorm*abNorm;
// make sure they are sorted
using std::swap;
@@ -1381,13 +1380,13 @@ public:
swap(tp, tq);
using std::min; using std::max;
- tp = min(abnorm2, max(0.0, tp));
- tq = max(0.0, min(abnorm2, tq));
+ tp = min(abNorm2, max(0.0, tp));
+ tq = max(0.0, min(abNorm2, tq));
if (abs(tp-tq) < eps2)
return false;
- intersection = Intersection({geo1.global(tp/abnorm2), geo1.global(tq/abnorm2)});
+ intersection = Intersection({geo1.global(tp/abNorm2), geo1.global(tq/abNorm2)});
return true;
}
};
--
GitLab