WIP: avoid division by zero if 3 points are aligned

In case of polyhedron similar to a box with 2 coplanar faces (i.e. a box's face split by an edge on its middle)
TODO: add a test

diff --git a/src/INTERP_KERNEL/InterpolationUtils.hxx b/src/INTERP_KERNEL/InterpolationUtils.hxx
index 901b27a23d61bea88c7dec6a18b98e681590ca02..ce8ae50556626ff830853fbfaaa5c63cbbb01ad0 100644 (file)
--- a/src/INTERP_KERNEL/InterpolationUtils.hxx
+++ b/src/INTERP_KERNEL/InterpolationUtils.hxx
@@ -1306,7 +1306,10 @@ namespace INTERP_KERNEL
     double XA_cross_XB[3] = {XA[1]*XB[2]-XA[2]*XB[1], XA[2]*XB[0]-XA[0]*XB[2], XA[0]*XB[1]-XA[1]*XB[0]};
     // norm is equal to double the area of the triangle
     double norm = std::sqrt(XA_cross_XB[0]*XA_cross_XB[0]+XA_cross_XB[1]*XA_cross_XB[1]+XA_cross_XB[2]*XA_cross_XB[2]);
-
+    // if 3 points A B X are aligned, norm is 0
+    // return 0 to avoid division by 0
+    if (norm<std::numeric_limits<double>::epsilon())
+      return 0;
     return ( XA_cross_XB[0]*XC[0]+ XA_cross_XB[1]*XC[1] + XA_cross_XB[2]*XC[2] ) / norm;
   }