From a4c6655fa780620ba1ebd57296f617a8fe20c98b Mon Sep 17 00:00:00 2001
From: Joachim Schoeberl <joachim.schoeberl@tuwien.ac.at>
Date: Fri, 29 Nov 2024 13:00:32 +0100
Subject: [PATCH] tolerance in CrossPointBarycentric

---
 libsrc/gprim/geom2d.cpp    | 5 ++++-
 libsrc/gprim/geom2d.hpp    | 2 +-
 libsrc/stlgeom/stltool.cpp | 2 +-
 3 files changed, 6 insertions(+), 3 deletions(-)

diff --git a/libsrc/gprim/geom2d.cpp b/libsrc/gprim/geom2d.cpp
index f50131aa..a6744979 100644
--- a/libsrc/gprim/geom2d.cpp
+++ b/libsrc/gprim/geom2d.cpp
@@ -128,7 +128,7 @@ Point2d CrossPoint (const Line2d & l1, const Line2d & l2)
 
 
 int CrossPointBarycentric (const Line2d & l1, const Line2d & l2,
-			   double & lam1, double & lam2)
+			   double & lam1, double & lam2, double eps)
 {
   // p = l1.1 + lam1 (l1.2-l1.1) = l2.1 + lam2 (l2.2-l2.1)
   double a11 = l1.p2.X() - l1.p1.X();
@@ -140,8 +140,11 @@ int CrossPointBarycentric (const Line2d & l1, const Line2d & l2,
   double b2 = l2.p1.Y() - l1.p1.Y();
   
   double det = a11*a22 - a12 * a21;
+  /*
   if (det == 0)
     return 1;
+  */
+  if (fabs (det) < eps * (fabs(a11*a22)+fabs(a12*a21))) return 1;
     
   lam1 = (a22 * b1 - a12 * b2) / det;
   lam2 = (a11 * b2 - a21 * b1) / det;
diff --git a/libsrc/gprim/geom2d.hpp b/libsrc/gprim/geom2d.hpp
index f28f968a..2c79ca4c 100644
--- a/libsrc/gprim/geom2d.hpp
+++ b/libsrc/gprim/geom2d.hpp
@@ -365,7 +365,7 @@ namespace netgen
     friend DLL_HEADER Point2d CrossPoint (const Line2d & l1, const Line2d & l2);
     /// returns 1 iff parallel
     friend int CrossPointBarycentric (const Line2d & l1, const Line2d & l2,
-				      double & lam1, double & lam2);
+				      double & lam1, double & lam2, double eps);
     
     ///
     friend int Parallel (const Line2d & l1, const Line2d & l2, double peps);
diff --git a/libsrc/stlgeom/stltool.cpp b/libsrc/stlgeom/stltool.cpp
index 67fec1bd..98dcfdc0 100644
--- a/libsrc/stlgeom/stltool.cpp
+++ b/libsrc/stlgeom/stltool.cpp
@@ -1369,7 +1369,7 @@ bool STLBoundary :: TestSegChartNV(const Point3d & p1, const Point3d& p2,
       Line2d l2 (sp1, sp2);
       double lam1, lam2;
 
-      int err = CrossPointBarycentric (l1, l2, lam1, lam2);
+      int err = CrossPointBarycentric (l1, l2, lam1, lam2, eps);
       bool in1 = (lam1 > eps) && (lam1 < 1-eps);
       bool on1 = (lam1 > -eps) && (lam1 < 1 + eps);
       bool in2 = (lam2 > eps) && (lam2 < 1-eps);