Add Polygon/Polyhedron NSIDED/NFACED ExodusII output

include/mesh/exodusII_io_helper.h

@@ -225,6 +225,11 @@ class ExodusII_IO_Helper : public ParallelObject
    */
   void read_elem_in_block(int block);
 
+  /**
+   * Reads NSIDED face blocks used by NFACED element blocks.
+   */
+  void read_face_blocks();
+
   /**
    * Read in edge blocks, storing information in the BoundaryInfo object.
    */
@@ -703,6 +708,13 @@ class ExodusII_IO_Helper : public ParallelObject
   // each block must equal num_edge.
   int & num_edge_blk;
 
+  // Total number of faces
+  int & num_face;
+
+  // Total number of face blocks. The sum of the number of faces in
+  // each block must equal num_face.
+  int & num_face_blk;
+
   // Total number of node sets
   int & num_node_sets;
 
@@ -748,6 +760,15 @@ class ExodusII_IO_Helper : public ParallelObject
   // Vector of nodes in an element
   std::vector<int> connect;
 
+  // For NSIDED element blocks, the number of nodes in each element
+  std::vector<int> elem_node_counts;
+
+  // For NFACED element blocks, the number of faces in each element
+  std::vector<int> elem_face_counts;
+
+  // For NFACED element blocks, the nodes in each Exodus face
+  std::vector<std::vector<int>> c0polyhedron_face_connect;
+
   // Vector of the sideset IDs
   std::vector<int> ss_ids;
 

include/mesh/exodus_header_info.h

@@ -57,7 +57,8 @@ class ExodusHeaderInfo
     // Pack individual integers into vector
     std::vector<int> buffer =
       {num_dim, num_elem, num_elem_blk, num_node_sets,
-       num_side_sets, num_elem_sets, num_edge_blk, num_edge};
+       num_side_sets, num_elem_sets, num_edge_blk, num_edge,
+       num_face_blk, num_face};
 
     // broadcast integers
     comm.broadcast(buffer);
@@ -72,6 +73,8 @@ class ExodusHeaderInfo
     num_elem_sets = buffer[ctr++];
     num_edge_blk  = buffer[ctr++];
     num_edge      = buffer[ctr++];
+    num_face_blk  = buffer[ctr++];
+    num_face      = buffer[ctr++];
   }
 
   std::vector<char> title;
@@ -84,6 +87,8 @@ class ExodusHeaderInfo
   int num_elem_sets;
   int num_edge_blk;
   int num_edge;
+  int num_face_blk;
+  int num_face;
 };
 
 } // namespace libMesh

src/geom/cell_c0polyhedron.C

@@ -18,6 +18,7 @@
 // Local includes
 #include "libmesh/cell_c0polyhedron.h"
 
+#include "libmesh/enum_io_package.h"
 #include "libmesh/enum_order.h"
 #include "libmesh/face_polygon.h"
 #include "libmesh/mesh_tools.h"
@@ -188,11 +189,44 @@ C0Polyhedron::edges_adjacent_to_node(const unsigned int n) const
 
 
 
-void C0Polyhedron::connectivity(const unsigned int /*sf*/,
-                                const IOPackage /*iop*/,
-                                std::vector<dof_id_type> & /*conn*/) const
+void C0Polyhedron::connectivity(const unsigned int sf,
+                                const IOPackage iop,
+                                std::vector<dof_id_type> & conn) const
 {
-  libmesh_not_implemented();
+  libmesh_assert_less (sf, this->n_sub_elem());
+  libmesh_assert_not_equal_to (iop, INVALID_IO_PACKAGE);
+
+  const auto & subtet = this->_triangulation[sf];
+
+  switch (iop)
+    {
+    case TECPLOT:
+      {
+        conn.resize(8);
+        conn[0] = this->node_id(subtet[0])+1;
+        conn[1] = this->node_id(subtet[1])+1;
+        conn[2] = this->node_id(subtet[2])+1;
+        conn[3] = this->node_id(subtet[2])+1;
+        conn[4] = this->node_id(subtet[3])+1;
+        conn[5] = this->node_id(subtet[3])+1;
+        conn[6] = this->node_id(subtet[3])+1;
+        conn[7] = this->node_id(subtet[3])+1;
+        return;
+      }
+
+    case VTK:
+      {
+        conn.resize(4);
+        conn[0] = this->node_id(subtet[0]);
+        conn[1] = this->node_id(subtet[1]);
+        conn[2] = this->node_id(subtet[2]);
+        conn[3] = this->node_id(subtet[3]);
+        return;
+      }
+
+    default:
+      libmesh_error_msg("Unsupported IO package " << iop);
+    }
 }
 
 

src/mesh/exodusII_io.C

@@ -20,13 +20,15 @@
 #include "libmesh/exodusII_io.h"
 
 #include "libmesh/boundary_info.h"
+#include "libmesh/cell_c0polyhedron.h"
 #include "libmesh/dof_map.h"
 #include "libmesh/dyna_io.h"  // ElementDefinition for BEX
 #include "libmesh/enum_elem_type.h"
 #include "libmesh/elem.h"
 #include "libmesh/enum_to_string.h"
 #include "libmesh/equation_systems.h"
 #include "libmesh/exodusII_io_helper.h"
+#include "libmesh/face_c0polygon.h"
 #include "libmesh/fpe_disabler.h"
 #include "libmesh/int_range.h"
 #include "libmesh/libmesh_logging.h"
@@ -452,20 +454,89 @@ void ExodusII_IO::read (const std::string & fname)
       // Set any relevant node/edge maps for this element
       const std::string type_str (exio_helper->get_elem_type());
       const auto & conv = exio_helper->get_conversion(type_str);
+      const bool is_c0polygon = (conv.libmesh_elem_type() == C0POLYGON);
+      const bool is_c0polyhedron = (conv.libmesh_elem_type() == C0POLYHEDRON);
+      std::size_t c0polygon_connect_offset = 0;
+      std::size_t c0polyhedron_connect_offset = 0;
 
       // Loop over all the faces in this block
       int jmax = nelem_last_block+exio_helper->num_elem_this_blk;
       for (int j=nelem_last_block; j<jmax; j++)
         {
-          auto uelem = Elem::build(conv.libmesh_elem_type());
-
           const int elem_num = j - nelem_last_block;
+          std::unique_ptr<Elem> uelem;
+
+          if (is_c0polygon)
+            {
+              const int n_elem_nodes = exio_helper->elem_node_counts[elem_num];
+              libmesh_error_msg_if(n_elem_nodes < 3,
+                                   "Error: Exodus NSIDED block element "
+                                   << elem_num
+                                   << " has only " << n_elem_nodes << " nodes.");
+
+              uelem = std::make_unique<C0Polygon>(cast_int<unsigned int>(n_elem_nodes));
+            }
+          else if (is_c0polyhedron)
+            {
+              const int n_faces = exio_helper->elem_face_counts[elem_num];
+              libmesh_error_msg_if(n_faces < 4,
+                                   "Error: Exodus NFACED block element "
+                                   << elem_num
+                                   << " has only " << n_faces << " faces.");
+
+              std::vector<std::shared_ptr<Polygon>>
+                sides(cast_int<std::size_t>(n_faces));
+
+              for (auto s : index_range(sides))
+                {
+                  libmesh_assert_less(c0polyhedron_connect_offset,
+                                      exio_helper->connect.size());
+                  const int exodus_face_id =
+                    exio_helper->connect[c0polyhedron_connect_offset++];
+
+                  libmesh_error_msg_if
+                    (exodus_face_id <= 0 ||
+                     exodus_face_id >
+                     cast_int<int>(exio_helper->c0polyhedron_face_connect.size()),
+                     "Error: Exodus NFACED block element "
+                     << elem_num
+                     << " references face ID " << exodus_face_id
+                     << ", but the file contains "
+                     << exio_helper->c0polyhedron_face_connect.size()
+                     << " faces.");
+
+                  const auto & face_nodes =
+                    exio_helper->c0polyhedron_face_connect[exodus_face_id - 1];
+                  auto side = std::make_shared<C0Polygon>
+                    (cast_int<unsigned int>(face_nodes.size()));
+
+                  for (auto n : index_range(face_nodes))
+                    {
+                      const auto libmesh_node_id =
+                        exio_helper->get_libmesh_node_id(face_nodes[n]);
+                      side->set_node(n, mesh.node_ptr(libmesh_node_id));
+                    }
+
+                  sides[s] = std::move(side);
+                }
+
+              std::unique_ptr<Node> mid_elem_node;
+              uelem = std::make_unique<C0Polyhedron>(sides, mid_elem_node);
+              if (mid_elem_node)
+                {
+                  Node * added_node = mesh.add_node(std::move(mid_elem_node));
+                  if (added_node->id() >= n_nodes)
+                    n_nodes = added_node->id() + 1;
+                }
+            }
+          else
+            uelem = Elem::build(conv.libmesh_elem_type());
 
           // Make sure that Exodus's number of nodes per Elem matches
           // the number of Nodes for this type of Elem. We only check
           // this for the first Elem in each block, since these values
           // are the same for every Elem in the block.
-          if (!elem_num)
+          if (!is_c0polygon && !is_c0polyhedron && !elem_num)
             libmesh_error_msg_if(exio_helper->num_nodes_per_elem != static_cast<int>(uelem->n_nodes()),
                                  "Error: Exodus file says "
                                  << exio_helper->num_nodes_per_elem
@@ -549,12 +620,26 @@ void ExodusII_IO::read (const std::string & fname)
           // If we don't have any Bezier extraction operators, this
           // is easy: we've already built all our nodes and just need
           // to link to them.
-          if (exio_helper->bex_cv_conn.empty())
+          if (is_c0polyhedron)
+            {
+              // Node pointers were already assigned while constructing
+              // the polygonal sides above.
+            }
+          else if (exio_helper->bex_cv_conn.empty())
             {
-              for (int k=0; k<exio_helper->num_nodes_per_elem; k++)
+              const int n_nodes_this_elem =
+                is_c0polygon ?
+                exio_helper->elem_node_counts[elem_num] :
+                exio_helper->num_nodes_per_elem;
+
+              for (int k=0; k<n_nodes_this_elem; k++)
                 {
                   // Get index into this block's connectivity array
-                  int gi = elem_num * exio_helper->num_nodes_per_elem + conv.get_node_map(k);
+                  std::size_t gi =
+                    is_c0polygon ?
+                    c0polygon_connect_offset++ :
+                    elem_num * exio_helper->num_nodes_per_elem + conv.get_node_map(k);
+                  libmesh_assert_less(gi, exio_helper->connect.size());
 
                   // Get the 1-based Exodus node id from the "connect" array
                   auto exodus_node_id = exio_helper->connect[gi];
@@ -689,6 +774,11 @@ void ExodusII_IO::read (const std::string & fname)
             }
         }
 
+      libmesh_assert(!is_c0polygon ||
+                     c0polygon_connect_offset == exio_helper->connect.size());
+      libmesh_assert(!is_c0polyhedron ||
+                     c0polyhedron_connect_offset == exio_helper->connect.size());
+
       // running sum of # of elements per block,
       // (should equal total number of elements in the end)
       nelem_last_block += exio_helper->num_elem_this_blk;