Merge pull request idaholab#28701 from DanielYankura/find_non_matchin…

…g_edges

framework/include/meshgenerators/MeshDiagnosticsGenerator.h

@@ -49,6 +49,8 @@ class MeshDiagnosticsGenerator : public MeshGenerator
   void checkNonConformalMeshFromAdaptivity(const std::unique_ptr<MeshBase> & mesh) const;
   /// Routine to check whether the Jacobians (elem and side) are not negative
   void checkLocalJacobians(const std::unique_ptr<MeshBase> & mesh) const;
+  //// Routine to check for non matching edges
+  void checkNonMatchingEdges(const std::unique_ptr<MeshBase> & mesh) const;
 
   /**
    * Utility routine to output the final diagnostics level in the desired mode
@@ -82,6 +84,10 @@ class MeshDiagnosticsGenerator : public MeshGenerator
   const MooseEnum _check_non_conformal_mesh;
   /// tolerance for detecting when meshes are not conformal
   const Real _non_conformality_tol;
+  //// whether to check for intersecting edges
+  const MooseEnum _check_non_matching_edges;
+  //// tolerance for detecting when edges intersect
+  const Real _non_matching_edge_tol;
   /// whether to check for the adaptivity of non-conformal meshes
   const MooseEnum _check_adaptivity_non_conformality;
   /// whether to check for negative jacobians in the domain

framework/include/utils/MeshBaseDiagnosticsUtils.h

@@ -24,4 +24,9 @@ void checkNonConformalMesh(const std::unique_ptr<libMesh::MeshBase> & mesh,
                            const unsigned int num_outputs,
                            const Real conformality_tol,
                            unsigned int & num_nonconformal_nodes);
+
+bool checkFirstOrderEdgeOverlap(const Elem & edge1,
+                                const Elem & edge2,
+                                Point & intersection_point,
+                                const Real intersection_tol);
 }

framework/src/meshgenerators/MeshDiagnosticsGenerator.C

@@ -20,6 +20,8 @@
 #include "libmesh/cell_tet4.h"
 #include "libmesh/face_quad4.h"
 #include "libmesh/cell_hex8.h"
+#include "libmesh/string_to_enum.h"
+#include "libmesh/enum_point_locator_type.h"
 
 registerMooseObject("MooseApp", MeshDiagnosticsGenerator);
 
@@ -64,6 +66,10 @@ MeshDiagnosticsGenerator::validParams()
   params.addParam<MooseEnum>("examine_non_conformality",
                              chk_option,
                              "whether to examine the conformality of elements in the mesh");
+  params.addParam<MooseEnum>("examine_non_matching_edges",
+                             chk_option,
+                             "Whether to check if there are any intersecting edges");
+  params.addParam<Real>("intersection_tol", TOLERANCE, "tolerence for intersecting edges");
   params.addParam<Real>("nonconformal_tol", TOLERANCE, "tolerance for element non-conformality");
   params.addParam<MooseEnum>(
       "search_for_adaptivity_nonconformality",
@@ -93,6 +99,8 @@ MeshDiagnosticsGenerator::MeshDiagnosticsGenerator(const InputParameters & param
     _check_non_planar_sides(getParam<MooseEnum>("examine_nonplanar_sides")),
     _check_non_conformal_mesh(getParam<MooseEnum>("examine_non_conformality")),
     _non_conformality_tol(getParam<Real>("nonconformal_tol")),
+    _check_non_matching_edges(getParam<MooseEnum>("examine_non_matching_edges")),
+    _non_matching_edge_tol(getParam<Real>("intersection_tol")),
     _check_adaptivity_non_conformality(
         getParam<MooseEnum>("search_for_adaptivity_nonconformality")),
     _check_local_jacobian(getParam<MooseEnum>("check_local_jacobian")),
@@ -111,7 +119,8 @@ MeshDiagnosticsGenerator::MeshDiagnosticsGenerator(const InputParameters & param
       _check_watertight_nodesets == "NO_CHECK" && _check_element_volumes == "NO_CHECK" &&
       _check_element_types == "NO_CHECK" && _check_element_overlap == "NO_CHECK" &&
       _check_non_planar_sides == "NO_CHECK" && _check_non_conformal_mesh == "NO_CHECK" &&
-      _check_adaptivity_non_conformality == "NO_CHECK" && _check_local_jacobian == "NO_CHECK")
+      _check_adaptivity_non_conformality == "NO_CHECK" && _check_local_jacobian == "NO_CHECK" &&
+      _check_non_matching_edges == "NO_CHECK")
     mooseError("You need to turn on at least one diagnostic. Did you misspell a parameter?");
 }
 
@@ -158,6 +167,9 @@ MeshDiagnosticsGenerator::generate()
   if (_check_local_jacobian != "NO_CHECK")
     checkLocalJacobians(mesh);
 
+  if (_check_non_matching_edges != "NO_CHECK")
+    checkNonMatchingEdges(mesh);
+
   return dynamic_pointer_cast<MeshBase>(mesh);
 }
 
@@ -1388,6 +1400,129 @@ MeshDiagnosticsGenerator::checkLocalJacobians(const std::unique_ptr<MeshBase> &
                  num_negative_side_qp_jacobians);
 }
 
+void
+MeshDiagnosticsGenerator::checkNonMatchingEdges(const std::unique_ptr<MeshBase> & mesh) const
+{
+  /*Algorithm Overview
+    1)Prechecks
+      a)This algorithm only works for 3D so check for that first
+    2)Loop
+      a)Loop through every element
+      b)For each element get the edges associated with it
+      c)For each edge check overlap with any edges of nearby elements
+      d)Have check to make sure the same pair of edges are not being tested twice for overlap
+    3)Overlap check
+      a)Shortest line that connects both lines is perpendicular to both lines
+      b)A good overview of the math for finding intersecting lines can be found
+    here->paulbourke.net/geometry/pointlineplane/
+  */
+  if (mesh->mesh_dimension() != 3)
+    mooseError("The edge intersection algorithm only works with 3D meshes");
+  if (!mesh->is_serial())
+    mooseError("Only serialized/replicated meshes are supported");
+  unsigned int num_intersecting_edges = 0;
+
+  // Create map of element to bounding box to avoing reinitializing the same bounding box multiple
+  // times
+  std::unordered_map<Elem *, BoundingBox> bounding_box_map;
+  for (const auto elem : mesh->active_element_ptr_range())
+  {
+    const auto boundingBox = elem->loose_bounding_box();
+    bounding_box_map.insert({elem, boundingBox});
+  }
+
+  std::unique_ptr<PointLocatorBase> point_locator = mesh->sub_point_locator();
+  std::set<std::array<dof_id_type, 4>> overlapping_edges_nodes;
+  for (const auto elem : mesh->active_element_ptr_range())
+  {
+    // loop through elem's nodes and find nearby elements with it
+    std::set<const Elem *> candidate_elems;
+    std::set<const Elem *> nearby_elems;
+    for (unsigned int i = 0; i < elem->n_nodes(); i++)
+    {
+      (*point_locator)(elem->point(i), candidate_elems);
+      nearby_elems.insert(candidate_elems.begin(), candidate_elems.end());
+    }
+    std::vector<std::unique_ptr<const Elem>> elem_edges(elem->n_edges());
+    for (auto i : elem->edge_index_range())
+      elem_edges[i] = elem->build_edge_ptr(i);
+    for (const auto other_elem : nearby_elems)
+    {
+      // If they're the same element then there's no need to check for overlap
+      if (elem->id() >= other_elem->id())
+        continue;
+
+      std::vector<std::unique_ptr<const Elem>> other_edges(other_elem->n_edges());
+      for (auto j : other_elem->edge_index_range())
+        other_edges[j] = other_elem->build_edge_ptr(j);
+      for (auto & edge : elem_edges)
+      {
+        for (auto & other_edge : other_edges)
+        {
+          // Get nodes from edges
+          const Node * n1 = edge->get_nodes()[0];
+          const Node * n2 = edge->get_nodes()[1];
+          const Node * n3 = other_edge->get_nodes()[0];
+          const Node * n4 = other_edge->get_nodes()[1];
+
+          // Create array<dof_id_type, 4> to check against set
+          std::array<dof_id_type, 4> node_id_array = {n1->id(), n2->id(), n3->id(), n4->id()};
+          std::sort(node_id_array.begin(), node_id_array.end());
+
+          // Check if the edges have already been added to our check_edges list
+          if (overlapping_edges_nodes.count(node_id_array))
+          {
+            continue;
+          }
+
+          // Check element/edge type
+          if (edge->type() != EDGE2)
+          {
+            std::string element_message = "Edge of type " + Utility::enum_to_string(edge->type()) +
+                                          " was found in cell " + std::to_string(elem->id()) +
+                                          " which is of type " +
+                                          Utility::enum_to_string(elem->type()) + '\n' +
+                                          "The edge intersection check only works for EDGE2 "
+                                          "elements.\nThis message will not be output again";
+            mooseDoOnce(_console << element_message << std::endl);
+            continue;
+          }
+          if (other_edge->type() != EDGE2)
+            continue;
+
+          // Now compare edge with other_edge
+          Point intersection_coords;
+          bool overlap = MeshBaseDiagnosticsUtils::checkFirstOrderEdgeOverlap(
+              *edge, *other_edge, intersection_coords, _non_matching_edge_tol);
+          if (overlap)
+          {
+            // Add the nodes that make up the 2 edges to the vector overlapping_edges_nodes
+            overlapping_edges_nodes.insert(node_id_array);
+            num_intersecting_edges += 2;
+            if (num_intersecting_edges < _num_outputs)
+            {
+              // Print error message
+              std::string elem_id = std::to_string(elem->id());
+              std::string other_elem_id = std::to_string(other_elem->id());
+              std::string x_coord = std::to_string(intersection_coords(0));
+              std::string y_coord = std::to_string(intersection_coords(1));
+              std::string z_coord = std::to_string(intersection_coords(2));
+              std::string message = "Intersecting edges found between elements " + elem_id +
+                                    " and " + other_elem_id + " near point (" + x_coord + ", " +
+                                    y_coord + ", " + z_coord + ")";
+              _console << message << std::endl;
+            }
+          }
+        }
+      }
+    }
+  }
+  diagnosticsLog("Number of intersecting element edges: " +
+                     Moose::stringify(num_intersecting_edges),
+                 _check_non_matching_edges,
+                 num_intersecting_edges);
+}
+
 void
 MeshDiagnosticsGenerator::diagnosticsLog(std::string msg,
                                          const MooseEnum & log_level,

framework/src/utils/MeshBaseDiagnosticsUtils.C

@@ -67,4 +67,66 @@ checkNonConformalMesh(const std::unique_ptr<MeshBase> & mesh,
   }
   pl->unset_close_to_point_tol();
 }
+
+bool
+checkFirstOrderEdgeOverlap(const Elem & edge1,
+                           const Elem & edge2,
+                           Point & intersection_point,
+                           const Real intersection_tol)
+{
+  // check that the two elements are of type EDGE2
+  mooseAssert(edge1.type() == EDGE2, "Elements must be of type EDGE2");
+  mooseAssert(edge2.type() == EDGE2, "Elements must be of type EDGE2");
+  // Get nodes from the two edges
+  const Point & p1 = edge1.point(0);
+  const Point & p2 = edge1.point(1);
+  const Point & p3 = edge2.point(0);
+  const Point & p4 = edge2.point(1);
+
+  // Check that the two edges are not sharing a node
+  if (p1 == p3 || p1 == p4 || p2 == p3 || p2 == p4)
+    return false;
+
+  /*
+    There's a chance that they overlap. Find shortest line that connects two edges and if its length
+    is close enough to 0 return true The shortest line between the two edges will be perpendicular
+    to both.
+  */
+  const auto d1343 = (p1 - p3) * (p4 - p3);
+  const auto d4321 = (p4 - p3) * (p2 - p1);
+  const auto d1321 = (p1 - p3) * (p2 - p1);
+  const auto d4343 = (p4 - p3) * (p4 - p3);
+  const auto d2121 = (p2 - p1) * (p2 - p1);
+
+  const auto denominator = d2121 * d4343 - d4321 * d4321;
+  const auto numerator = d1343 * d4321 - d1321 * d4343;
+
+  if (std::fabs(denominator) < intersection_tol)
+    // This indicates that the two lines are parallel so they don't intersect
+    return false;
+
+  const auto mua = numerator / denominator;
+  const auto mub = (d1343 + (mua * d4321)) / d4343;
+
+  // Use these values to solve for the two points that define the shortest line segment
+  const auto pa = p1 + mua * (p2 - p1);
+  const auto pb = p3 + mub * (p4 - p3);
+
+  // This method assume the two lines are infinite. This check to make sure na and nb are part of
+  // their respective line segments
+  if (mua < 0 || mua > 1)
+    return false;
+  if (mub < 0 || mub > 1)
+    return false;
+
+  // Calculate distance between these two nodes
+  const auto distance = (pa - pb).norm();
+  if (distance < intersection_tol)
+  {
+    intersection_point = pa;
+    return true;
+  }
+  else
+    return false;
+}
 }

test/tests/meshgenerators/mesh_diagnostics_generator/all_at_once.i

@@ -16,6 +16,7 @@
     examine_sidesets_orientation = WARNING
     check_for_watertight_sidesets = WARNING
     check_for_watertight_nodesets = WARNING
+    examine_non_matching_edges = WARNING
     search_for_adaptivity_nonconformality = WARNING
     check_local_jacobian = WARNING
   []

test/tests/meshgenerators/mesh_diagnostics_generator/intersecting_edge_test.i

@@ -0,0 +1,67 @@
+[Mesh]
+  [gmg1]
+    type = GeneratedMeshGenerator
+    dim = 3
+    nx = 2
+    ny = 2
+    nz = 2
+    xmin = -0.5
+    xmax = 0.5
+    ymin = -0.5
+    ymax = 0.5
+    zmin = -1.5
+    zmax = -0.5
+  []
+  [gmg2]
+    type = GeneratedMeshGenerator
+    dim = 3
+    nx = 2
+    ny = 2
+    nz = 2
+    xmin = -0.5
+    xmax = 0.5
+    ymin = -0.5
+    ymax = 0.5
+    zmin = -0.5
+    zmax = 0.5
+  []
+  [block_1]
+    type = ParsedSubdomainMeshGenerator
+    input = gmg1
+    combinatorial_geometry = 'z < -0.5'
+    block_id = 1
+  []
+  [block_2]
+    type = ParsedSubdomainMeshGenerator
+    input = gmg2
+    combinatorial_geometry = 'z > 0.5'
+    block_id = 2
+  []
+  [convert1]
+    type = ElementsToTetrahedronsConverter
+    input = 'block_1'
+  []
+  [convert2]
+    type = ElementsToTetrahedronsConverter
+    input = 'block_2'
+  []
+  [rotate]
+    type = TransformGenerator
+    input = convert2
+    transform = 'rotate'
+    vector_value = '180 0 0'
+  []
+  [stitch]
+    type = StitchedMeshGenerator
+    inputs = 'convert1 rotate'
+    stitch_boundaries_pairs = 'front back'
+  []
+  [diag]
+    type = MeshDiagnosticsGenerator
+    input = stitch
+    examine_non_matching_edges = INFO
+  []
+[]
+[Outputs]
+  exodus = true
+[]

test/tests/meshgenerators/mesh_diagnostics_generator/tests

@@ -13,6 +13,14 @@
       expect_out = 'Number of non-conformal nodes: 3'
       detail = 'element overlapping,'
     []
+    [intersecting_edges]
+      type = RunApp
+      input = 'intersecting_edge_test.i'
+      cli_args = '--mesh-only'
+      mesh_mode = replicated
+      expect_out = 'Number of intersecting element edges: 4'
+      detail = 'edges intersecting'
+    []
     [inconsistent_sidesets]
       type = RunApp
       input = 'consistent_domains.i'
@@ -269,7 +277,8 @@
                   Mesh/diag/examine_nonplanar_sides=NO_CHECK
                   Mesh/diag/examine_sidesets_orientation=NO_CHECK
                   Mesh/diag/check_for_watertight_sidesets=NO_CHECK
-                  Mesh/diag/check_for_watertight_nodesets=NO_CHECK
+                  Mesh/diag/check_for_watertight_nodesets=NO_CHECK   
+                  Mesh/diag/examine_non_matching_edges=NO_CHECK
                   Mesh/diag/search_for_adaptivity_nonconformality=NO_CHECK
                   Mesh/diag/check_local_jacobian=NO_CHECK --mesh-only"
       detail = 'a diagnostics object is created but no diagnostics are requested.'