First stab at importing tracing routines

include/geom/intersection_tools.h

@@ -119,6 +119,30 @@ bool within_edge_on_side(const Elem & elem,
                          const bool linearize = false,
                          const Real tol = TOLERANCE);
 
+/**
+ * \returns The index of the vertex that the given point is at
+ * on the given element, invalid otherwise
+ * @param elem The element
+ * @param p The point
+ * @param tol The relative tolerance used in point comparison
+*/
+unsigned short at_vertex(const Elem & elem,
+                         const Point & p,
+                         const Real tol = TOLERANCE);
+
+/**
+ * \returns The index of the vertex that the given point is at
+ * on the given element on the given side, invalid otherwise
+ * @param elem The element
+ * @param p The point
+ * @param s The side
+ * @param tol The relative tolerance used in point comparison
+*/
+unsigned short at_vertex_on_side(const Elem & elem,
+                                 const Point & p,
+                                 const unsigned short s,
+                                 const Real tol = TOLERANCE);
+
 namespace detail
 {
 /**

src/geom/intersection_tools.C

@@ -150,6 +150,32 @@ bool within_edge(const Elem & elem,
   return false;
 }
 
+unsigned short at_vertex(const Elem & elem,
+                         const Point & p,
+                         const Real tol)
+{
+  for (const auto v : elem.vertex_index_range())
+    if (elem.point(v).relative_fuzzy_equals(p, tol))
+      return v;
+  return Elem::invalid_vertex;
+}
+
+unsigned short at_vertex_on_side(const Elem & elem,
+                                 const Point & p,
+                                 const unsigned short s,
+                                 const Real tol)
+{
+  libmesh_assert_less(s, elem.n_sides());
+  const unsigned int * side_node_map = elem.nodes_on_side_ptr(s);
+  for (const auto i : make_range(elem.n_vertices_on_side(s)))
+  {
+    const auto v = side_node_map[i];
+    if (elem.point(v).relative_fuzzy_equals(p, tol))
+      return v;
+  }
+  return Elem::invalid_vertex;
+}
+
 namespace detail
 {
 bool _within_edge(const Elem & elem,

tests/geom/intersection_tools_test.C

@@ -229,10 +229,56 @@ public:
       }
   }
 
+  void test_at_vertex()
+  {
+    LOG_UNIT_TEST;
+
+    // check locations at every node
+    for (const auto & mesh : _meshes)
+      for (const auto & elem : mesh->active_local_element_ptr_range())
+        for (const auto n : elem->node_index_range())
+        {
+          const auto & p = elem->point(n);
+          const auto is_vertex = elem->is_vertex(n);
+          const auto at_vertex = IntersectionTools::at_vertex(*elem, p);
+
+          if (is_vertex)
+            CPPUNIT_ASSERT_EQUAL(at_vertex, (unsigned short)n);
+          else
+            CPPUNIT_ASSERT_EQUAL(at_vertex, Elem::invalid_vertex);
+
+          for (const auto s : elem->side_index_range())
+          {
+            const auto at_vertex_on_side = IntersectionTools::at_vertex_on_side(*elem, p, s);
+            if (is_vertex && elem->is_node_on_side(n, s))
+              CPPUNIT_ASSERT_EQUAL(at_vertex_on_side, at_vertex);
+            else
+              CPPUNIT_ASSERT_EQUAL(at_vertex_on_side, Elem::invalid_vertex);
+          }
+        }
+
+    // check edge centroids
+    for (const auto & mesh : _meshes)
+      if (mesh->mesh_dimension() > 1)
+        for (const auto & elem : mesh->active_local_element_ptr_range())
+          for (const auto e : elem->edge_index_range())
+          {
+            const auto edge = elem->build_edge_ptr(e);
+            const auto at_vertex = IntersectionTools::at_vertex(*elem, edge->vertex_average());
+            CPPUNIT_ASSERT_EQUAL(at_vertex, Elem::invalid_vertex);
+          }
+
+    // check elem centroids
+    for (const auto & mesh : _meshes)
+      for (const auto & elem : mesh->active_local_element_ptr_range())
+        CPPUNIT_ASSERT_EQUAL(IntersectionTools::at_vertex(*elem, elem->centroid()), Elem::invalid_vertex);
+  }
+
 };
 
-#define MESHEDINTERSECTIONTOOLSTEST               \
-  CPPUNIT_TEST( test_within_edge );
+#define MESHEDINTERSECTIONTOOLSTEST  \
+  CPPUNIT_TEST( test_within_edge );  \
+  CPPUNIT_TEST( test_at_vertex );
 
 #define INSTANTIATE_MESHEDINTERSECTIONTOOLSTEST(elemtype)                                       \
   class MeshedIntersectionToolsTest_##elemtype : public MeshedIntersectionToolsTest<elemtype> { \