Make collinear more robust with a relative tolerance

libMesh · Nov 14, 2022 · be46863 · be46863
1 parent 65d923a
commit be46863
Show file tree

Hide file tree

Showing 2 changed files with 10 additions and 7 deletions.
diff --git a/include/geom/intersection_tools.h b/include/geom/intersection_tools.h
@@ -60,7 +60,8 @@ WithinSegmentResult within_segment(const Point & s1,
  * @param p1 The first point
  * @param p2 The second point
  * @param p3 The third point
- * @param tol The tolerance to use
+ * @param tol The tolerance to use (absolute tolerance of the cosine
+ * of the angle between [p1 -> p2] and [p1 -> p3])
  * @return Whether or not the given points are collinear
  */
 bool collinear(const Point & p1,

diff --git a/src/geom/intersection_tools.C b/src/geom/intersection_tools.C
@@ -63,12 +63,14 @@ bool collinear(const Point & p1,
                const Point & p3,
                const Real tol)
 {
-  // (p1 -> p2) X (p1 - > p3) should be == the zero vector
-  const auto p1p2_cross_p1p3 = (p2 - p1).cross(p3 - p1);
-  for (const auto i : make_range(LIBMESH_DIM))
-    if (std::abs(p1p2_cross_p1p3(i)) > tol)
-      return false;
-  return true;
+  const auto p2_p1 = p1 - p2;
+  const auto p2_p1_norm = p2_p1.norm();
+  const auto p3_p2 = p2 - p3;
+  const auto p3_p2_norm = p3_p2.norm();
+  const auto denom = p2_p1_norm * p3_p2_norm;
+  if (denom == 0) // same points
+    return true;
+  return (std::abs(p2_p1 * p3_p2) / denom) > ((Real)1 - tol);
 }
 
 bool within_edge_on_side(const Elem & elem,