Add heuristics to tetrahedralize hexahedrons with high valence polygons on two sides

that arise from extrusion
- dont create a 0 vol elements if we have 3 surrounding nodes and everybody is coplanar
- try to avoid this situation by improving as many potential 0-volume neighbor tets as possible
  when selecting which tet to construct

Author: GiudGiud

src/geom/cell_c0polyhedron.C

@@ -373,7 +373,7 @@ void C0Polyhedron::retriangulate()
   // sometimes with a similarly-sorted vector of surrounding elements
   // to go with it.
   auto surroundings_of =
-    [&nodes_to_elem_map, & surface]
+    [& nodes_to_elem_map, & surface]
     (const Node & node,
      std::vector<Elem *> * surrounding_elems)
     {
@@ -439,7 +439,11 @@ void C0Polyhedron::retriangulate()
             v02 = static_cast<Point>(*surrounding_nodes[n+1]) - node,
             v03 = static_cast<Point>(*surrounding_nodes[n+2]) - node;
 
-          total_solid_angle += solid_angle(v01, v02, v03);
+          // Solid angles are all negative because the surface triangulation triangles were pointing
+          // outwards and the 4th node is inward because the polyhedron is convex
+          // Some are zero when all the nodes selected are from the same (planar) polygon
+          // Make it positive to store smallest angle first
+          total_solid_angle += std::abs(solid_angle(v01, v02, v03));
         }
 
       return std::make_pair(n_surrounding, total_solid_angle);
@@ -456,8 +460,11 @@ void C0Polyhedron::retriangulate()
   // in parallel.
   typedef std::multimap<std::pair<int, Real>, Node*> node_map_type;
   node_map_type nodes_by_geometry;
+  // Keep track of index in map
   std::map<Node *, node_map_type::iterator> node_index;
 
+  // The map is first sorted by the valence
+  // then by the angle. Max angle, if they are negative (positive right now), first
   for (auto node : surface.node_ptr_range())
     node_index[node] =
       nodes_by_geometry.emplace(geometry_at(*node), node);
@@ -466,6 +473,9 @@ void C0Polyhedron::retriangulate()
   // each vertex, and an inner loop to remove multiple tetrahedra in
   // cases where the vertex has more than 3 neighboring triangles.
 
+  // One of the checks does not work if we skip adding a 0 volume tet
+  bool has_skipped_adding_tets_and_retriangulating = false;
+
   // We'll be done when there are only three "unremoved" nodes left,
   // so they don't actually enclose any volume.
   for (auto i : make_range(nodes_by_geometry.size()-3))
@@ -525,8 +535,9 @@ void C0Polyhedron::retriangulate()
       std::vector<Real> local_tet_quality(n_surrounding, 1);
 
       // From our center node with N surrounding nodes we can make N-2
-      // tetrahedra.  The first N-3 each replace two surface tets with
-      // two new surface tets.
+      // tetrahedra.
+      // For the first N-3 tets we build, we may have to replace two surface
+      // triangles with two new surface triangles.
       //
       // My first idea was to greedily pick nodes with the smallest
       // local (solid) angles to get the best quality.  This works in
@@ -539,15 +550,28 @@ void C0Polyhedron::retriangulate()
       // My third idea is to try to fix the lowest quality tets first,
       // by picking cases where they have higher quality neighbors,
       // and creating those neighbors so as to change them.
+      //
+      // A fourth idea was to make sure to fix as many neighbor tets as
+      // possible.
 
       auto find_new_tet_nodes =
-        [& local_tet_quality, & find_valid_nodes_around]
+        [& local_tet_quality, & find_valid_nodes_around, & surrounding_nodes]
         ()
       {
         unsigned int jbest = 0;
         auto [jminus, jplus] = find_valid_nodes_around(jbest);
         Real qneighbest = std::min(local_tet_quality[jminus],
                                    local_tet_quality[jplus]);
+        // Count number of bad neighbors
+        auto num_bad_neigh_best = (local_tet_quality[jminus] <= 0) + (local_tet_quality[jplus] <= 0);
+
+        // Count the number of valid surrounding nodes
+        unsigned int n_valid_surrounding = 0;
+        for (const auto surr_node : surrounding_nodes)
+          if (surr_node)
+            n_valid_surrounding++;
+
+        // Expected qualities post removing jbest from the surrounding nodes
         for (auto j : make_range(std::size_t(1),
                                  local_tet_quality.size()))
           {
@@ -565,77 +589,88 @@ void C0Polyhedron::retriangulate()
                 qneighj > 0)
               continue;
 
+            const auto num_bad_neigh = (local_tet_quality[jminus] <= 0) +
+                                       (local_tet_quality[jplus] <= 0);
             // We want to try for the best possible fix.
-            if ((local_tet_quality[j] - qneighj) >
-                (local_tet_quality[jbest] - qneighj))
+            // - The more bad (volume <= 0) neighbors we fix the better
+            // - If same number of bad neighbors, pick best quality one
+            if ((num_bad_neigh >= num_bad_neigh_best) &&
+                ((local_tet_quality[j] - qneighj) > (local_tet_quality[jbest] - qneighj)))
               {
                 jbest = j;
                 qneighbest = qneighj;
+                num_bad_neigh_best = num_bad_neigh;
               }
           }
 
+        // If there are only 3 neighors and the tet is bad, we will skip it downstream
         libmesh_error_msg_if
-          (local_tet_quality[jbest] <= 0,
+          (n_valid_surrounding > 3 && local_tet_quality[jbest] <= 0,
            "Cannot build non-singular non-inverted tet");
 
         std::tie(jminus, jplus) = find_valid_nodes_around(jbest);
 
         return std::make_tuple(jbest, jminus, jplus);
       };
 
-      if (n_surrounding > 3)
-        {
-          // We'll be searching local_tet_quality even after tet
-          // extraction disconnects us from some nodes; when we do we
-          // don't want to get one.
-          constexpr Real far_node = -1e6;
-
-          // Vectors from the center node to each of its surrounding
-          // nodes are helpful for calculating prospective tet
-          // quality.
-          std::vector<Point> v0s(n_surrounding);
-          for (auto j : make_range(n_surrounding))
-            v0s[j] = *(Point *)surrounding_nodes[j] - *node;
-
-          // Find the tet quality we'd potentially get from each
-          // possible choice of tet
-          auto local_tet_quality_of =
-            [& surrounding_nodes, & v0s, & find_valid_nodes_around]
-            (unsigned int j)
-          {
-            auto [jminus, jplus] = find_valid_nodes_around(j);
-
-            // Anything proportional to the ratio of volume to
-            // total-edge-length-cubed should peak for perfect tets
-            // but hit 0 for pancakes and slivers.
-
-            const Real total_len =
-              v0s[j].norm() + v0s[jminus].norm() + v0s[jplus].norm() +
-              (*(Point *)surrounding_nodes[jplus] -
-               *(Point *)surrounding_nodes[j]).norm() +
-              (*(Point *)surrounding_nodes[j] -
-               *(Point *)surrounding_nodes[jminus]).norm() +
-              (*(Point *)surrounding_nodes[jminus] -
-               *(Point *)surrounding_nodes[jplus]).norm();
-
-            // Orientation here is tricky.  Think of the triple
-            // product as (v1 cross v2) dot v3, with right hand rule.
-            const Real six_vol =
-              triple_product(v0s[jminus], v0s[jplus], v0s[j]);
+      // Vectors from the center node to each of its surrounding
+      // nodes are helpful for calculating prospective tet
+      // quality.
+      std::vector<Point> v0s(n_surrounding);
+      for (auto j : make_range(n_surrounding))
+        v0s[j] = *(Point *)surrounding_nodes[j] - *node;
+
+      // Find the tet quality we'd potentially get from each
+      // possible choice of tet
+      auto local_tet_quality_of =
+        [& surrounding_nodes, & v0s, & find_valid_nodes_around]
+        (unsigned int j)
+      {
+        auto [jminus, jplus] = find_valid_nodes_around(j);
+
+        // Anything proportional to the ratio of volume to
+        // total-edge-length-cubed should peak for perfect tets
+        // but hit 0 for pancakes and slivers.
+
+        const Real total_len =
+          v0s[j].norm() + v0s[jminus].norm() + v0s[jplus].norm() +
+          (*(Point *)surrounding_nodes[jplus] -
+            *(Point *)surrounding_nodes[j]).norm() +
+          (*(Point *)surrounding_nodes[j] -
+            *(Point *)surrounding_nodes[jminus]).norm() +
+          (*(Point *)surrounding_nodes[jminus] -
+            *(Point *)surrounding_nodes[jplus]).norm();
+
+        // Orientation here is tricky.  Think of the triple
+        // product as (v1 cross v2) dot v3, with right hand rule.
+        const Real six_vol =
+          triple_product(v0s[jminus], v0s[jplus], v0s[j]);
+
+        return six_vol / (total_len * total_len * total_len);
+      };
 
-            return six_vol / (total_len * total_len * total_len);
-          };
+      // Get the quality of the tets for every node around the node:
+      // we have chosen 1 node. If we chose another node for the tet, there
+      // are only two other nodes that are also neighbors, hence 1 tet per node
+      // Note: we always need to know the quality, because as we delete vertices
+      // we could end up with a a node with 3 surrounding nodes, but all coplanar!
+      for (auto j : make_range(n_surrounding))
+        local_tet_quality[j] = local_tet_quality_of(j);
 
-          for (auto j : make_range(n_surrounding))
-            local_tet_quality[j] = local_tet_quality_of(j);
+      // We'll be searching local_tet_quality even after tet
+      // extraction disconnects us from some nodes; when we do we
+      // don't want to re-use it so we mark the quality as -1e6
+      constexpr Real far_node = -1e6;
 
-          // If we have N surrounding nodes, we can make N tets and
+      if (n_surrounding > 3)
+        {
+          // If we have N surrounding nodes, we can make N-3 tets and
           // that'll bring us back to the 3-surrounding-node case to
           // finish.
           for (auto t : make_range(n_surrounding-3))
             {
               libmesh_ignore(t);
-
+              
               auto [jbest, jminus, jplus] = find_new_tet_nodes();
 
               // Turn these four nodes into a tet
@@ -697,6 +732,7 @@ void C0Polyhedron::retriangulate()
               // unchanged until we're out of this inner loop; let's
               // recalculate it here, and then we'll be done with it.
               Node * & nbestref = surrounding_nodes[jbest];
+              libmesh_assert(nbestref);
               nodes_by_geometry.erase(node_index[nbestref]);
               node_index[nbestref] =
                 nodes_by_geometry.emplace(geometry_at(*nbestref), nbestref);
@@ -724,6 +760,18 @@ void C0Polyhedron::retriangulate()
       // counterclockwise" we get from the lambda.
       auto [j2, j1, j3] = find_new_tet_nodes();
 
+      // If the final tet has zero volume, skip, we have covered all we need for this node
+      // This occurs when considering a polygonal prism.
+      // The opposite node on the other extruded face was deleted first
+      // and depending on the triangulation, this leaves the node only attached
+      // to 3 other coplanar nodes on one of the side
+      if (local_tet_quality[j2] == 0 || local_tet_quality[j2] == far_node)
+      {
+        nodes_by_geometry.erase(geometry_it);
+        has_skipped_adding_tets_and_retriangulating = true;
+        continue;
+      }
+
       // Turn these last four nodes into a tet
       Node * n1 = surrounding_nodes[j1],
            * n2 = surrounding_nodes[j2],
@@ -793,7 +841,8 @@ void C0Polyhedron::retriangulate()
     }
 
   // At this point our surface should just have two triangles left.
-  libmesh_assert_equal_to(surface.n_elem(), 2);
+  if (!has_skipped_adding_tets_and_retriangulating)
+    libmesh_assert_equal_to(surface.n_elem(), 2);
 }
 
 

src/geom/cell_polyhedron.C

@@ -115,6 +115,7 @@ Polyhedron::Polyhedron (const std::vector<std::shared_ptr<Polygon>> & sides,
       {
         const Polygon & side = *sides[s];
         const Point x_i = side.point(0);
+        // opposite of normal
         const Point n_i =
           (side.point(1) - side.point(0)).cross
           (side.point(0) - side.point(side.n_sides()-1)).unit();
@@ -123,7 +124,7 @@ Polyhedron::Polyhedron (const std::vector<std::shared_ptr<Polygon>> & sides,
         inward_normal = (n_i * (center - x_i) > TOLERANCE);
       }
 
-    // We're betting a lot on "our polyhedra are all convex", so let's
+    // We're betting a lot on "our polygon sides are all convex", so let's
     // check that if we have time.
 #ifdef DEBUG
     for (unsigned int s : index_range(sides))