Changes to C0 polyhedron tetrahedralization loops

include/geom/cell_c0polyhedron.h

@@ -28,6 +28,8 @@
 namespace libMesh
 {
 
+class Node;
+
 /**
  * The \p C0Polyhedron is an element in 3D with an arbitrary (but fixed)
  * number of polygonal first-order (C0Polygon) sides.
@@ -55,9 +57,15 @@ class C0Polyhedron : public Polyhedron
    * We'll set up a simple default tetrahedralization here, but if users
    * want to adjust node positions later they'll want to retriangulate()
    * after.
+   * @param mid_elem_node a reference to the unique pointer set to the mid-element node,
+   *        only if the default (optimal) tetrahedralization algorithm fails, and the backup
+   *        trivial algorithm is used to tetrahedralize the element.
+   *        It's the caller's job to add this node to the mesh after!
+   * @param p pointer to the parent element of this one (useful for h-refinement)
    */
   explicit
   C0Polyhedron (const std::vector<std::shared_ptr<Polygon>> & sides,
+                std::unique_ptr<Node> & mid_elem_node,
                 Elem * p=nullptr);
 
   C0Polyhedron (C0Polyhedron &&) = delete;
@@ -87,7 +95,7 @@ class C0Polyhedron : public Polyhedron
    * \returns the number of vertices.  For the simple C0Polyhedron
    * every node is a vertex.
    */
-  virtual unsigned int n_vertices() const override final { return this->_nodelinks_data.size(); }
+  virtual unsigned int n_vertices() const override final { return this->_nodelinks_data.size() - _has_mid_elem_node; }
 
   /**
    * \returns the number of tetrahedra to break this into for
@@ -129,6 +137,8 @@ class C0Polyhedron : public Polyhedron
   virtual bool is_node_on_edge(const unsigned int n,
                                const unsigned int e) const override;
 
+  virtual std::vector<unsigned int> edges_adjacent_to_node(const unsigned int n) const override;
+
   /**
    * \returns \p true if the element map is definitely affine within
    * numerical tolerances.  We don't even share a master element from
@@ -170,6 +180,18 @@ class C0Polyhedron : public Polyhedron
 
   Point side_vertex_average_normal(const unsigned int s) const override final;
 
+  /**
+   * \returns the local side-indices of subelement sides of the
+   * polyhedron
+   * Each subelement here is a tetrahedron
+   * @param i the index of the sub-element
+   * The return array is indexed by sides of the subelement, and contains the
+   * index of the side of the polyhedron if the sub-element side is on that side,
+   * or libMesh::invalid_uint if the subelement side is internal to the polyhedron
+   * NOTE: this routine involves a loop over the sides of the polyhedron
+   */
+  virtual std::array<int, 4> subelement_sides_to_poly_sides(unsigned int i) const;
+
   /**
    * Create a triangulation (tetrahedralization) based on the current
    * sides' triangulations.
@@ -193,6 +215,9 @@ class C0Polyhedron : public Polyhedron
                                  const unsigned int /*k*/) const override
   { libmesh_not_implemented(); return 0; }
 
+  /// Whether we have a mid element node
+  bool _has_mid_elem_node;
+
 #endif // LIBMESH_ENABLE_AMR
 
 };

include/systems/generic_projector.h

@@ -1611,7 +1611,8 @@ void GenericProjector<FFunctor, GFunctor, FValue, ProjectionAction>::SortAndCopy
 
       // If we have non-vertex nodes, the first is an edge node, but
       // if we're in 2D we'll call that a side node
-      const bool has_edge_nodes = (n_nodes > n_vertices && dim > 2);
+      const bool has_poly_midnode = (elem->type() == C0POLYHEDRON) && (n_nodes == n_vertices + 1);
+      const bool has_edge_nodes = (n_nodes > n_vertices + has_poly_midnode && dim > 2);
 
       // If we have even more nodes, the next is a side node.
       const bool has_side_nodes =
@@ -1733,7 +1734,18 @@ void GenericProjector<FFunctor, GFunctor, FValue, ProjectionAction>::SortAndCopy
             }
         };
 
-      for (unsigned int v=0; v != n_vertices; ++v)
+      // Treat polyhedron midnode as a vertex
+      // NOTE: if we start having edge or side nodes on polyhedra, we need to use that +1 offset
+      // in the edge and side nodes code as well!
+#ifndef NDEBUG
+      for (auto v_num : this->projector.variables)
+        {
+          const auto & family = this->projector.system.variable(v_num).type().family;
+          // Add to the list once known to be correctly functioning with the midnode
+          libmesh_assert(!has_poly_midnode || (family == LAGRANGE || family == MONOMIAL || family == XYZ));
+        }
+#endif
+      for (const auto v : make_range(n_vertices + has_poly_midnode))
         {
           const Node * node = elem->node_ptr(v);
 

src/geom/cell_c0polyhedron.C

@@ -23,6 +23,7 @@
 #include "libmesh/mesh_tools.h"
 #include "libmesh/replicated_mesh.h"
 #include "libmesh/tensor_value.h"
+#include "libmesh/node.h"
 
 // C++ headers
 #include <numeric> // std::iota
@@ -34,10 +35,14 @@ namespace libMesh
 // C0Polyhedron class member functions
 
 C0Polyhedron::C0Polyhedron
-  (const std::vector<std::shared_ptr<Polygon>> & sides, Elem * p) :
-  Polyhedron(sides, p)
+  (const std::vector<std::shared_ptr<Polygon>> & sides, std::unique_ptr<Node> & mid_elem_node, Elem * p) :
+  Polyhedron(sides, p), _has_mid_elem_node(false)
 {
   this->retriangulate();
+
+  // Grab the last node
+  if (_has_mid_elem_node)
+    mid_elem_node.reset(this->_nodelinks_data.back());
 }
 
 
@@ -50,8 +55,14 @@ std::unique_ptr<Elem> C0Polyhedron::disconnected_clone() const
 {
   auto sides = this->side_clones();
 
-  std::unique_ptr<Elem> returnval =
-    std::make_unique<C0Polyhedron>(sides);
+  std::unique_ptr<Node> mid_elem_node;
+  std::unique_ptr<Elem> returnval = std::make_unique<C0Polyhedron>(sides, mid_elem_node);
+  // Swap out the new mid elem node with the previous one
+  if (mid_elem_node)
+  {
+    libmesh_assert(returnval->n_nodes() == this->n_nodes());
+    returnval->set_node(returnval->n_nodes() - 1, this->_nodelinks_data.back());
+  }
 
   returnval->set_id() = this->id();
 #ifdef LIBMESH_ENABLE_UNIQUE_ID
@@ -71,11 +82,14 @@ std::unique_ptr<Elem> C0Polyhedron::disconnected_clone() const
 
 
 
-bool C0Polyhedron::is_vertex(const unsigned int libmesh_dbg_var(i)) const
+bool C0Polyhedron::is_vertex(const unsigned int i) const
 {
-  libmesh_assert (i < this->n_nodes());
+  libmesh_assert_less (i, this->n_nodes());
 
-  return true;
+  if (i < this->n_vertices())
+    return true;
+  else
+    return false;
 }
 
 
@@ -155,6 +169,16 @@ bool C0Polyhedron::is_node_on_edge(const unsigned int n,
 
 
 
+std::vector<unsigned int>
+C0Polyhedron::edges_adjacent_to_node(const unsigned int n) const
+{
+  if (n == n_nodes() - 1)
+    return {};
+  return Polyhedron::edges_adjacent_to_node(n);
+}
+
+
+
 void C0Polyhedron::connectivity(const unsigned int /*sf*/,
                                 const IOPackage /*iop*/,
                                 std::vector<dof_id_type> & /*conn*/) const
@@ -270,6 +294,40 @@ C0Polyhedron::side_vertex_average_normal(const unsigned int s) const
 }
 
 
+std::array<int, 4>
+C0Polyhedron::subelement_sides_to_poly_sides(unsigned int tri_i) const
+{
+  std::array<int, 4> sides = {invalid_int, invalid_int, invalid_int, invalid_int};
+  libmesh_assert(tri_i < this->n_subelements());
+  const auto & tet_nodes = _triangulation[tri_i];
+  for (const auto poly_side : make_range(n_sides()))
+  {
+    const auto sd_nodes = nodes_on_side(poly_side);
+    bool zero_in = std::find(sd_nodes.begin(), sd_nodes.end(), tet_nodes[0]) != sd_nodes.end();
+    bool one_in = std::find(sd_nodes.begin(), sd_nodes.end(), tet_nodes[1]) != sd_nodes.end();
+    bool two_in = std::find(sd_nodes.begin(), sd_nodes.end(), tet_nodes[2]) != sd_nodes.end();
+    bool three_in = std::find(sd_nodes.begin(), sd_nodes.end(), tet_nodes[3]) != sd_nodes.end();
+
+    // Take advantage of the fact that the poly side can only contain one tet side
+    // Note: we could optimize by replacing the booleans with the searches above
+    if (zero_in)
+    {
+      if (one_in)
+      {
+        if (two_in)
+          sides[0] = poly_side;
+        else if (three_in)
+          sides[1] = poly_side;
+      }
+      else if (two_in && three_in)
+        sides[3] = poly_side;
+    }
+    else if (three_in && two_in && one_in)
+      sides[2] = poly_side;
+  }
+  return sides;
+}
+
 
 void C0Polyhedron::retriangulate()
 {
@@ -351,6 +409,13 @@ void C0Polyhedron::retriangulate()
   verify_surface();
 #endif
 
+  // First heuristic: try with no interior point
+  // This might not succeed, not every surface triangulation gives a tetrahedralization
+  // with no additional interior point
+  // But if it succeeds, it uses less tetrahedra to cut the polyhedron
+#ifdef LIBMESH_ENABLE_EXCEPTIONS
+  try
+  {
   // We'll have to edit this as we change the surface elements, but we
   // have a method to initialize it easily.
   std::vector<std::vector<dof_id_type>> nodes_to_elem_vec_map;
@@ -791,9 +856,54 @@ void C0Polyhedron::retriangulate()
       // eliminate again.
       nodes_by_geometry.erase(geometry_it);
     }
+    // At this point our surface should just have two triangles left.
+    libmesh_assert_equal_to(surface.n_elem(), 2);
+    _has_mid_elem_node = false;
+  }
+  // Failed without an interior point.
+  // Use a single vertex-average interior point and tetrahedralize around it
+  catch ( libMesh::LogicError & )
+#endif
+  {
+    // Clear the triangulation we started building
+    this->_triangulation.clear();
+
+    // Get the vertex-average, no need to triangulate for this
+    const auto v_avg = this->vertex_average();
+    if (!_has_mid_elem_node)
+    {
+      // Create the mid element node. Add it to nodelinks
+      // We'll attach a unique ptr to it later
+      Node * mid_elem_node = new Node(v_avg);
+      this->_nodelinks_data.push_back(mid_elem_node);
+      _has_mid_elem_node = true;
+    }
+    else
+    {
+      // We are triangulating for a second time, we have already added this mid-element node
+      libmesh_assert(n_vertices() == n_nodes() - 1);
+    }
 
-  // At this point our surface should just have two triangles left.
-  libmesh_assert_equal_to(surface.n_elem(), 2);
+    // Build the tetrahedralization with each of the triangles on each side
+    for (unsigned int s : make_range(this->n_sides()))
+    {
+      const auto & [side, inward_normal, node_map] = this->_sidelinks_data[s];
+
+      for (auto t : make_range(side->n_subtriangles()))
+      {
+        // Get all the nodes
+        const auto & n1 = node_map[side->subtriangle(t)[0]];
+        const auto & n2 = node_map[side->subtriangle(t)[1]];
+        const auto & n3 = node_map[side->subtriangle(t)[2]];
+
+        // The mid node is the last node in the _nodes array
+        if (!inward_normal)
+          this->add_tet((int)n1, (int)n2, (int)n3, this->n_nodes() - 1);
+        else
+          this->add_tet((int)n1, (int)n3, (int)n2, this->n_nodes() - 1);
+      }
+    }
+  }
 }
 
 
@@ -808,13 +918,15 @@ void C0Polyhedron::add_tet(int n1,
   libmesh_assert_less(n2, nn);
   libmesh_assert_less(n3, nn);
   libmesh_assert_less(n4, nn);
+#endif
 
   const Point v12 = this->point(n2) - this->point(n1);
   const Point v13 = this->point(n3) - this->point(n1);
   const Point v14 = this->point(n4) - this->point(n1);
   const Real six_vol = triple_product(v12, v13, v14);
-  libmesh_assert_greater(six_vol, Real(0));
-#endif
+  // We need to error on this in optimized modes to fall back onto the
+  // tetrahedralization with a mid node
+  libmesh_error_msg_if(six_vol <= 0, "Creating flat tet");
 
   this->_triangulation.push_back({n1, n2, n3, n4});
 }

src/geom/cell_polyhedron.C

@@ -93,6 +93,10 @@ Polyhedron::Polyhedron (const std::vector<std::shared_ptr<Polygon>> & sides,
           }
       }
 
+    // Plan room for an extra node so we don't invalidate this->_nodes when we add to
+    // nodelinks_data in the derived class
+    _nodelinks_data.reserve(_nodelinks_data.size() + 1);
+
     // Do the manual initialization that Elem::Elem and Cell::Cell
     // couldn't, now that we've resized both our vectors.  No need to
     // manually set nullptr, though, since std::vector does that.
@@ -684,6 +688,7 @@ std::array<Point, 4> Polyhedron::master_subelement (unsigned int i) const
 }
 
 
+
 std::tuple<unsigned int, Real, Real, Real>
 Polyhedron::subelement_coordinates (const Point & p, Real tol) const
 {

tests/fe/fe_lagrange_test.C

@@ -43,4 +43,7 @@ INSTANTIATE_FETEST(SECOND, LAGRANGE, PYRAMID14);
 INSTANTIATE_FETEST(THIRD, LAGRANGE, PYRAMID18);
 
 INSTANTIATE_FETEST(FIRST, LAGRANGE, C0POLYHEDRON);
+// We don't match the linear solution at this time
+// struct MidNode {};
+// INSTANTIATE_FETEST_CASE(FIRST, LAGRANGE, C0POLYHEDRON, MidNode);
 #endif

tests/fe/fe_monomial_test.C

@@ -66,4 +66,6 @@ INSTANTIATE_FETEST(FIFTH, MONOMIAL, PRISM18);
 INSTANTIATE_FETEST(FIFTH, MONOMIAL, PRISM21);
 
 INSTANTIATE_FETEST(FIRST, MONOMIAL, C0POLYHEDRON);
+struct MidNode {};
+INSTANTIATE_FETEST_CASE(FIRST, MONOMIAL, C0POLYHEDRON, MidNode);
 #endif