Add IntersectionTools::within_edge_on_side

libMesh · Nov 10, 2022 · 76ac32d · 76ac32d
1 parent 8f7d2f7
commit 76ac32d
Show file tree

Hide file tree

Showing 3 changed files with 244 additions and 0 deletions.
diff --git a/include/geom/intersection_tools.h b/include/geom/intersection_tools.h
@@ -25,6 +25,8 @@
 namespace libMesh
 {
 
+class Elem;
+class ElemExtrema;
 class Point;
 
 namespace IntersectionTools
@@ -82,6 +84,31 @@ bool collinear(const Point & p1,
                const Point & p3,
                const Real tol = TOLERANCE);
 
+/**
+ * \returns True if the given point is contained within an edge
+ * on the given side of an element
+ * @param elem The element
+ * @param p The point
+ * @param s The side
+ * @param extrema "Extrema" to be filled with the edge that the point
+ * is within, if any (must be initially invalid)
+ * @param linearize Whether or not to "linearize" the check, if this
+ * is set to false and edges are found to not be collinear, an error
+ * is thrown
+ *
+ * \p extrema will be set to an "at vertex" state if the point is
+ * both within the edge _and_ at a vertex.
+ *
+ * This method is only implemented for three-dimensional, finite
+ * elements.
+*/
+bool within_edge_on_side(const Elem & elem,
+                         const Point & p,
+                         const unsigned short s,
+                         ElemExtrema & extrema,
+                         const bool linearize = false,
+                         const Real tol = TOLERANCE);
+
 } // namespace IntersectionTools
 } // namespace libMesh
 

diff --git a/src/geom/intersection_tools.C b/src/geom/intersection_tools.C
@@ -21,6 +21,8 @@
 
 #include "libmesh/point.h"
 #include "libmesh/int_range.h"
+#include "libmesh/elem.h"
+#include "libmesh/elem_extrema.h"
 
 namespace libMesh::IntersectionTools
 {
@@ -76,6 +78,71 @@ bool collinear(const Point & p1,
   return true;
 }
 
+bool within_edge_on_side(const Elem & elem,
+                         const Point & p,
+                         const unsigned short s,
+                         ElemExtrema & extrema,
+                         const bool linearize,
+                         const Real tol)
+{
+  libmesh_assert_less(s, elem.n_sides());
+  libmesh_assert(extrema.is_invalid());
+  libmesh_assert_equal_to(elem.dim(), 3);
+
+  // For higher order than linear without linearization, make sure
+  // that our edges are collinear
+  if (elem.default_order() > 1 && !linearize)
+    for (const auto e : elem.edge_index_range())
+    {
+      // we should expect 3 edges for our higher order elems
+      libmesh_assert_equal_to(elem.n_nodes_on_edge(e), 3);
+
+      const unsigned int * edge_nodes_map = elem.nodes_on_edge_ptr(e);
+      if (!collinear(elem.point(edge_nodes_map[0]),
+                     elem.point(edge_nodes_map[1]),
+                     elem.point(edge_nodes_map[2])))
+        libmesh_error_msg("Failed to use Cell::without_edge_on_side without linearization "
+                          "because an edge was found that is not collinear.");
+    }
+
+  const auto vs = elem.n_vertices_on_side(s);
+  const unsigned int * side_nodes_map = elem.nodes_on_side_ptr(s);
+
+  // side_nodes_map will point to something like [v0, v1, v2, v3]
+  // With the loop below, we're going to check (in this order):
+  // [v3 -> v0], [v0 -> v1], [v1 -> v2], [v2 -> v3]
+  auto last_v = side_nodes_map[vs - 1];
+  for (const auto side_v : make_range(vs))
+  {
+    const auto other_v = side_nodes_map[side_v];
+    const auto within_result = within_segment(elem.point(last_v), elem.point(other_v), p, tol);
+    if (within_result == NOT_WITHIN)
+    {
+      last_v = side_nodes_map[side_v];
+      continue;
+    }
+
+    if (within_result == BETWEEN)
+    {
+      extrema.set_edge(last_v, other_v);
+      libmesh_assert(extrema.build_edge(elem)->close_to_point(p, tol));
+    }
+    else if (within_result == AT_BEGINNING)
+    {
+      extrema.set_vertex(last_v);
+      libmesh_assert(elem.point(last_v).absolute_fuzzy_equals(p, tol));
+    }
+    else
+    {
+      extrema.set_vertex(other_v);
+      libmesh_assert(elem.point(other_v).absolute_fuzzy_equals(p, tol));
+    }
+    return true;
+  }
+
+  return false;
+}
+
 } // namespace libMesh::IntersectionTools
 
 
diff --git a/tests/geom/intersection_tools_test.C b/tests/geom/intersection_tools_test.C
@@ -1,8 +1,15 @@
 #include "libmesh_cppunit.h"
 
+#include "test_comm.h"
+
 #include <libmesh/intersection_tools.h>
 #include <libmesh/point.h>
 #include <libmesh/int_range.h>
+#include <libmesh/enum_elem_type.h>
+#include <libmesh/mesh.h>
+#include <libmesh/mesh_generation.h>
+#include <libmesh/elem.h>
+#include <libmesh/elem_extrema.h>
 
 using namespace libMesh;
 
@@ -80,3 +87,146 @@ public:
 };
 
 CPPUNIT_TEST_SUITE_REGISTRATION( IntersectionToolsTest );
+
+template <ElemType elem_type>
+class MeshedIntersectionToolsTest : public CppUnit::TestCase {
+
+private:
+  std::unique_ptr<Mesh> _mesh;
+
+protected:
+  std::string libmesh_suite_name;
+
+public:
+  void setUp()
+  {
+    const Real minpos = 1.5, maxpos = 4.86;
+    const unsigned int N = 3;
+
+    _mesh = std::make_unique<Mesh>(*TestCommWorld);
+    auto test_elem = Elem::build(elem_type);
+
+    const unsigned int dim = test_elem->dim();
+    const unsigned int use_y = dim > 1;
+    const unsigned int use_z = dim > 2;
+
+    MeshTools::Generation::build_cube (*_mesh,
+                                       N, N*use_y, N*use_z,
+                                       minpos, maxpos,
+                                       minpos, use_y*maxpos,
+                                       minpos, use_z*maxpos,
+                                       elem_type);
+  }
+
+  void test_within_edge_on_side()
+  {
+    LOG_UNIT_TEST;
+
+    if (_mesh->mesh_dimension() != 3)
+      return;
+
+    // check locations at every node
+    for (const auto & elem : _mesh->active_local_element_ptr_range())
+      for (const auto s : elem->side_index_range())
+        for (const auto e : elem->edge_index_range())
+          for (const auto n : elem->nodes_on_edge(e))
+          {
+            ElemExtrema extrema;
+            const auto within = IntersectionTools::within_edge_on_side(*elem,
+                                                                        elem->point(n),
+                                                                        s,
+                                                                        extrema);
+
+            CPPUNIT_ASSERT_EQUAL(within, elem->is_node_on_side(n, s));
+            if (elem->is_node_on_side(n, s))
+            {
+              CPPUNIT_ASSERT_EQUAL(elem->is_vertex(n), extrema.at_vertex(n));
+              CPPUNIT_ASSERT_EQUAL(elem->is_vertex(n), !extrema.at_edge(*elem, e));
+            }
+          }
+
+    // cut edges into segments
+    for (const auto & elem : _mesh->active_local_element_ptr_range())
+      for (const auto e : elem->edge_index_range())
+        for (const auto s : elem->side_index_range())
+          if (elem->is_edge_on_side(e, s))
+          {
+            const auto nodes_on_edge = elem->nodes_on_edge(e);
+            const auto & p1 = elem->point(nodes_on_edge[0]);
+            const auto & p2 = elem->point(nodes_on_edge[1]);
+            const auto length_vec = p2 - p1;
+            const auto length = length_vec.norm();
+            const auto p1_to_p2 = length_vec / length;
+
+            int segments = 5;
+            Real dx = (Real)1 / segments * length;
+            for (const auto i : make_range(-1, segments + 1))
+            {
+              const auto p = p1 + Real(i) * dx * p1_to_p2;
+              ElemExtrema extrema;
+              const auto within = IntersectionTools::within_edge_on_side(*elem,
+                                                                         p,
+                                                                         s,
+                                                                         extrema);
+
+              CPPUNIT_ASSERT_EQUAL(within, i >= 0 && i <= segments);
+              CPPUNIT_ASSERT_EQUAL(extrema.at_vertex(nodes_on_edge[0]), i == 0);
+              CPPUNIT_ASSERT_EQUAL(extrema.at_vertex(nodes_on_edge[1]), i == segments);
+              CPPUNIT_ASSERT_EQUAL(extrema.at_edge(*elem, e), i > 0 && i < segments);
+            }
+          }
+
+    // check elem centroids
+    for (const auto & elem : _mesh->active_local_element_ptr_range())
+      for (const auto s : elem->side_index_range())
+      {
+        ElemExtrema extrema;
+        CPPUNIT_ASSERT(!IntersectionTools::within_edge_on_side(*elem,
+                                                               elem->vertex_average(),
+                                                               s,
+                                                               extrema));
+      }
+  }
+
+};
+
+#define MESHEDINTERSECTIONTOOLSTEST               \
+  CPPUNIT_TEST( test_within_edge_on_side );
+
+#define INSTANTIATE_MESHEDINTERSECTIONTOOLSTEST(elemtype)                                       \
+  class MeshedIntersectionToolsTest_##elemtype : public MeshedIntersectionToolsTest<elemtype> { \
+  public:                                                                                       \
+  MeshedIntersectionToolsTest_##elemtype() :                                                    \
+    MeshedIntersectionToolsTest<elemtype>() {                                                   \
+    if (unitlog->summarized_logs_enabled())                                                     \
+      this->libmesh_suite_name = "MeshedIntersectionToolsTest";                                 \
+    else                                                                                        \
+      this->libmesh_suite_name = "MeshedIntersectionToolsTest_" #elemtype;                      \
+  }                                                                                             \
+  CPPUNIT_TEST_SUITE( MeshedIntersectionToolsTest_##elemtype );                                 \
+  MESHEDINTERSECTIONTOOLSTEST;                                                                  \
+  CPPUNIT_TEST_SUITE_END();                                                                     \
+  };                                                                                            \
+                                                                                                \
+  CPPUNIT_TEST_SUITE_REGISTRATION( MeshedIntersectionToolsTest_##elemtype )
+
+
+#if LIBMESH_DIM > 2
+INSTANTIATE_MESHEDINTERSECTIONTOOLSTEST(TET4);
+INSTANTIATE_MESHEDINTERSECTIONTOOLSTEST(TET10);
+INSTANTIATE_MESHEDINTERSECTIONTOOLSTEST(TET14);
+
+INSTANTIATE_MESHEDINTERSECTIONTOOLSTEST(HEX8);
+INSTANTIATE_MESHEDINTERSECTIONTOOLSTEST(HEX20);
+INSTANTIATE_MESHEDINTERSECTIONTOOLSTEST(HEX27);
+
+INSTANTIATE_MESHEDINTERSECTIONTOOLSTEST(PRISM6);
+INSTANTIATE_MESHEDINTERSECTIONTOOLSTEST(PRISM15);
+INSTANTIATE_MESHEDINTERSECTIONTOOLSTEST(PRISM18);
+INSTANTIATE_MESHEDINTERSECTIONTOOLSTEST(PRISM20);
+INSTANTIATE_MESHEDINTERSECTIONTOOLSTEST(PRISM21);
+
+INSTANTIATE_MESHEDINTERSECTIONTOOLSTEST(PYRAMID5);
+INSTANTIATE_MESHEDINTERSECTIONTOOLSTEST(PYRAMID13);
+INSTANTIATE_MESHEDINTERSECTIONTOOLSTEST(PYRAMID14);
+#endif // LIBMESH_DIM > 2