Add IntersectionTools and within_segment

Author: loganharbour

include/geom/intersection_tools.h

@@ -0,0 +1,75 @@
+// The libMesh Finite Element Library.
+// Copyright (C) 2002-2022 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
+
+// This library is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Lesser General Public
+// License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+
+// This library is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+// Lesser General Public License for more details.
+
+// You should have received a copy of the GNU Lesser General Public
+// License along with this library; if not, write to the Free Software
+// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
+
+
+
+#ifndef LIBMESH_INTERSECTIONTOOLS_H
+#define LIBMESH_INTERSECTIONTOOLS_H
+
+#include "libmesh/libmesh_common.h"
+
+namespace libMesh
+{
+
+class Point;
+
+namespace IntersectionTools
+{
+
+/**
+ * Enum that represents the result of the \p within_segment() methods.
+*/
+enum WithinSegmentResult : int {
+                         NOT_WITHIN = 0,
+                         AT_BEGINNING = 1,
+                         AT_END = 2,
+                         BETWEEN = 3 };
+
+/**
+ * Checks whether or not a point is within a line segment
+ * @param s1 The first point on the segment
+ * @param s2 The second point on the segment
+ * @param length The segment length (for optimization if it's already computed)
+ * @param p The point
+ * @param tol The tolerance to use
+ * @return Enum denoting whether or not the point is not within, at s1,
+ * at s2, or between [s1, s2]
+ */
+WithinSegmentResult within_segment(const Point & s1,
+                                   const Point & s2,
+                                   const Real length,
+                                   const Point & p,
+                                   const Real tol = TOLERANCE);
+
+/**
+ * Checks whether or not a point is within a line segment
+ * @param s1 The first point on the segment
+ * @param s2 The second point on the segment
+ * @param p The point
+ * @param tol The tolerance to use
+ * @return Enum denoting whether or not the point is not within, at s1,
+ * at s2, or between [s1, s2]
+ */
+WithinSegmentResult within_segment(const Point & s1,
+                                   const Point & s2,
+                                   const Point & p,
+                                   const Real tol = TOLERANCE);
+
+} // namespace IntersectionTools
+} // namespace libMesh
+
+#endif // LIBMESH_INTERSECTIONTOOLS_H

src/geom/intersection_tools.C

@@ -0,0 +1,67 @@
+// The libMesh Finite Element Library.
+// Copyright (C) 2002-2022 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
+
+// This library is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Lesser General Public
+// License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+
+// This library is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+// Lesser General Public License for more details.
+
+// You should have received a copy of the GNU Lesser General Public
+// License along with this library; if not, write to the Free Software
+// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
+
+
+
+#include "libmesh/intersection_tools.h"
+
+#include "libmesh/point.h"
+
+namespace libMesh::IntersectionTools
+{
+
+WithinSegmentResult within_segment(const Point & s1,
+                                   const Point & s2,
+                                   const Real length,
+                                   const Point & p,
+                                   const Real tol)
+{
+  libmesh_assert(!s1.absolute_fuzzy_equals(s2, tol));
+  libmesh_assert_less(std::abs((s1 - s2).norm() - length), tol);
+
+  const auto diff1 = p - s1;
+  const auto diff2 = p - s2;
+  const auto tol_scaled = tol * length;
+
+  if (diff1 * diff2 > tol_scaled)
+    return NOT_WITHIN;
+
+  const auto diff1_norm = diff1.norm();
+  if (diff1_norm < tol_scaled)
+    return AT_BEGINNING;
+
+  const auto diff2_norm = diff2.norm();
+  if (diff2_norm < tol_scaled)
+    return AT_END;
+
+  // whether or not p is _between_ [s1, s2]
+  if (std::abs(diff1_norm + diff2_norm - length) < tol_scaled)
+    return BETWEEN;
+  return NOT_WITHIN;
+}
+
+WithinSegmentResult within_segment(const Point & s1,
+                                   const Point & s2,
+                                   const Point & p,
+                                   const Real tol)
+{
+  return within_segment(s1, s2, (s1 - s2).norm(), p, tol);
+}
+
+} // namespace libMesh::IntersectionTools
+
+

tests/geom/intersection_tools_test.C

@@ -0,0 +1,52 @@
+#include "libmesh_cppunit.h"
+
+#include <libmesh/intersection_tools.h>
+#include <libmesh/point.h>
+#include <libmesh/int_range.h>
+
+using namespace libMesh;
+
+class IntersectionToolsTest : public CppUnit::TestCase
+{
+
+public:
+  LIBMESH_CPPUNIT_TEST_SUITE( IntersectionToolsTest );
+  CPPUNIT_TEST( within_segment );
+  CPPUNIT_TEST_SUITE_END();
+
+public:
+
+  void within_segment()
+  {
+    LOG_UNIT_TEST;
+
+    const Point s1(1.0, 2.0, 3.0);
+    const Point s2(2.0, 3.0, 4.0);
+    const auto length_vec = s2 - s1;
+    const auto length = length_vec.norm();
+    const auto s1_to_s2 = length_vec / length;
+
+    int segments = 10;
+    Real dx = (Real)1 / segments * length;
+    for (const auto i : make_range(-1, segments + 1))
+    {
+      const auto p = s1 + Real(i) * dx * s1_to_s2;
+      IntersectionTools::WithinSegmentResult within_result = IntersectionTools::WithinSegmentResult::NOT_WITHIN;
+      if (i == 0)
+        within_result = IntersectionTools::WithinSegmentResult::AT_BEGINNING;
+      else if (i > 0 && i < segments)
+        within_result = IntersectionTools::WithinSegmentResult::BETWEEN;
+      else if (i == segments)
+        within_result = IntersectionTools::WithinSegmentResult::AT_END;
+
+      CPPUNIT_ASSERT_EQUAL(IntersectionTools::within_segment(s1, s2, length, p), within_result);
+      CPPUNIT_ASSERT_EQUAL(IntersectionTools::within_segment(s1, s2, p), within_result);
+    }
+
+    CPPUNIT_ASSERT_EQUAL(IntersectionTools::within_segment(s1, s2, Point(9.9, 5, 3)),
+                         IntersectionTools::WithinSegmentResult::NOT_WITHIN);
+  }
+
+};
+
+CPPUNIT_TEST_SUITE_REGISTRATION( IntersectionToolsTest );