From c320536d2d28fc7c636b0e231bbb1af1c60ec695 Mon Sep 17 00:00:00 2001
From: BRUNO LATHUILIERE <bruno.lathuiliere@edf.fr>
Date: Wed, 21 Feb 2024 12:58:25 +0100
Subject: [PATCH] [TetraIntersect] Avoid multiple divisions in
 calculatePolygonBarycenter

+ more accurate (if no overflow) and faster
+ cosmetics floating point fixes
---
 src/INTERP_KERNEL/TransformedTriangle.cxx     | 21 +++++++++----------
 src/INTERP_KERNEL/TransformedTriangle.hxx     |  2 +-
 src/INTERP_KERNEL/TransformedTriangleMath.cxx | 10 ++++-----
 3 files changed, 15 insertions(+), 18 deletions(-)

diff --git a/src/INTERP_KERNEL/TransformedTriangle.cxx b/src/INTERP_KERNEL/TransformedTriangle.cxx
index 65ada4f2c..4b1892045 100644
--- a/src/INTERP_KERNEL/TransformedTriangle.cxx
+++ b/src/INTERP_KERNEL/TransformedTriangle.cxx
@@ -646,20 +646,19 @@ namespace INTERP_KERNEL
       const std::size_t m = polygon.size();
 
       for(int j = 0 ; j < 3 ; ++j)
-        {
-          barycenter[j] = 0.0;
-        }
+	barycenter[j] = 0.0;
+
+      for(std::size_t i = 0 ; i < m ; ++i)
+	{
+	  const double* pt = polygon[i];
+	  for(int j = 0 ; j < 3 ; ++j)
+	    barycenter[j] += pt[j];
+	}
 
       if(m != 0)
         {
-          for(std::size_t i = 0 ; i < m ; ++i)
-            {
-              const double* pt = polygon[i];
-              for(int j = 0 ; j < 3 ; ++j)
-                {
-                  barycenter[j] += pt[j] / double(m);
-                }
-            }
+	  for(int j = 0 ; j < 3 ; ++j)
+	    barycenter[j] = barycenter[j] / double(m);
         }
       LOG(3,"Barycenter is " << vToStr(barycenter));
     }
diff --git a/src/INTERP_KERNEL/TransformedTriangle.hxx b/src/INTERP_KERNEL/TransformedTriangle.hxx
index b2ce44315..14d396597 100644
--- a/src/INTERP_KERNEL/TransformedTriangle.hxx
+++ b/src/INTERP_KERNEL/TransformedTriangle.hxx
@@ -512,7 +512,7 @@ namespace INTERP_KERNEL
       if (std::fabs(t1+t2) < THRESHOLD_F*MULT_PREC_F)
         return false;
 
-    return (t1*t2 <= 0.0) && !epsilonEqual(t1,t2, (double)MULT_PREC_F);
+    return (t1*t2 <= 0.0) && !epsilonEqual(t1,t2, MULT_PREC_F);
   }
 
   inline bool TransformedTriangle::testFacetSurroundsSegment(const TriSegment seg, const TetraFacet facet) const
diff --git a/src/INTERP_KERNEL/TransformedTriangleMath.cxx b/src/INTERP_KERNEL/TransformedTriangleMath.cxx
index 30ecbc621..7b3eba4e8 100644
--- a/src/INTERP_KERNEL/TransformedTriangleMath.cxx
+++ b/src/INTERP_KERNEL/TransformedTriangleMath.cxx
@@ -331,7 +331,7 @@ namespace INTERP_KERNEL
           {
             // this value will not be used - we set it to whatever
             LOG(6, "Triple product not calculated for corner " << corner );
-            _tripleProducts[corner] = -3.14159265;
+            _tripleProducts[corner] = std::nan("triplep");
             _validTP[corner] = false;
           }
         anglesForRows.clear();
@@ -390,11 +390,9 @@ namespace INTERP_KERNEL
     
     //? is this more stable? -> no subtraction
     //    return asin( dotProd / ( lenNormal * lenEdgeVec ) ) + 3.141592625358979 / 2.0;
-    double tmp=dotProd / ( lenNormal * lenEdgeVec );
-    tmp=std::max(tmp,-1.);
-    tmp=std::min(tmp,1.);
-    return atan(1.0)*4.0 - acos(tmp);
-
+    const double tmp=dotProd / ( lenNormal * lenEdgeVec );
+    const double safe_tmp=std::max(std::min(tmp,1.),-1.);
+    return M_PI - std::acos(safe_tmp);
   }
 
   /**
-- 
2.39.2