Merge pull request #3491 from jwpeterson/match_lower_order_edges

ExodusII_IO: Allow low-order edge blocks to refer to higher-order Elem edges

Author: jwpeterson

src/mesh/exodusII_io_helper.C

@@ -1444,6 +1444,30 @@ void ExodusII_IO_Helper::read_edge_blocks(MeshBase & mesh)
         // Creates vector if not already there
         auto & vec = edge_map[edge_key];
         vec.emplace_back(elem->id(), e);
+
+        // If edge_ptr is a higher-order Elem (EDGE3 or higher) then also add
+        // a map entry for the lower-order (EDGE2) element which has matching
+        // vertices. This allows us to match lower-order edge blocks to edges
+        // of higher-order 3D elems (e.g. HEX20, TET10) and simplifies the
+        // definition of edge blocks.
+        if (edge_ptr->default_order() != FIRST)
+          {
+            // Construct a temporary low-order edge so that we can compute its key()
+            auto low_order_edge =
+              Elem::build(Elem::first_order_equivalent_type(edge_ptr->type()));
+
+            // Assign node pointers to low-order edge
+            for (unsigned int v=0; v<edge_ptr->n_vertices(); ++v)
+              low_order_edge->set_node(v) = edge_ptr->node_ptr(v);
+
+            // Compute the key for the temporary low-order edge we just built
+            dof_id_type low_order_edge_key = low_order_edge->key();
+
+            // Add this key to the map associated with the same (elem,
+            // edge) pair as the higher-order edge
+            auto & low_order_vec = edge_map[low_order_edge_key];
+            low_order_vec.emplace_back(elem->id(), e);
+          }
       }
 
   // Get reference to the mesh's BoundaryInfo object, as we will be
@@ -1533,9 +1557,12 @@ void ExodusII_IO_Helper::read_edge_blocks(MeshBase & mesh)
                     build_edge_ptr(edge_ptr, elem_edge_pair.second);
 
                   // Determine whether this candidate edge is a "real" match,
-                  // i.e. also has the same orientation.
+                  // i.e. also has the same orientation. Note that here we
+                  // only check that the vertices match regardless of how many
+                  // nodes the edge has, which allows us to match a lower-order
+                  // edge to a higher-order Elem.
                   bool is_match = true;
-                  for (int n=0; n<num_nodes_per_edge; ++n)
+                  for (unsigned int n=0; n<edge_ptr->n_vertices(); ++n)
                     if (edge_ptr->node_id(n) != edge->node_id(n))
                       {
                         is_match = false;

tests/Makefile.am

@@ -144,6 +144,7 @@ data = meshes/1_quad.bxt.gz \
        meshes/bad_64bit_elem_integers.e \
        meshes/good_32bit_elem_integers.e \
        meshes/mesh_assign_test_mesh.xda \
+       meshes/mesh_with_low_order_edge_blocks.e \
        meshes/shark_tooth_tri6.xda.gz \
        meshes/tetgen_one_tet10.ele \
        meshes/tetgen_one_tet10.node

tests/Makefile.in

@@ -264,6 +264,7 @@ am__unit_tests_dbg_SOURCES_DIST = driver.C libmesh_cppunit.h \
 	meshes/bad_64bit_elem_integers.e \
 	meshes/good_32bit_elem_integers.e \
 	meshes/mesh_assign_test_mesh.xda \
+	meshes/mesh_with_low_order_edge_blocks.e \
 	meshes/shark_tooth_tri6.xda.gz meshes/tetgen_one_tet10.ele \
 	meshes/tetgen_one_tet10.node fparser/autodiff.C
 am__dirstamp = $(am__leading_dot)dirstamp
@@ -453,6 +454,7 @@ am__unit_tests_devel_SOURCES_DIST = driver.C libmesh_cppunit.h \
 	meshes/bad_64bit_elem_integers.e \
 	meshes/good_32bit_elem_integers.e \
 	meshes/mesh_assign_test_mesh.xda \
+	meshes/mesh_with_low_order_edge_blocks.e \
 	meshes/shark_tooth_tri6.xda.gz meshes/tetgen_one_tet10.ele \
 	meshes/tetgen_one_tet10.node fparser/autodiff.C
 @LIBMESH_ENABLE_FPARSER_TRUE@am__objects_4 = fparser/unit_tests_devel-autodiff.$(OBJEXT)
@@ -636,6 +638,7 @@ am__unit_tests_oprof_SOURCES_DIST = driver.C libmesh_cppunit.h \
 	meshes/bad_64bit_elem_integers.e \
 	meshes/good_32bit_elem_integers.e \
 	meshes/mesh_assign_test_mesh.xda \
+	meshes/mesh_with_low_order_edge_blocks.e \
 	meshes/shark_tooth_tri6.xda.gz meshes/tetgen_one_tet10.ele \
 	meshes/tetgen_one_tet10.node fparser/autodiff.C
 @LIBMESH_ENABLE_FPARSER_TRUE@am__objects_6 = fparser/unit_tests_oprof-autodiff.$(OBJEXT)
@@ -819,6 +822,7 @@ am__unit_tests_opt_SOURCES_DIST = driver.C libmesh_cppunit.h \
 	meshes/bad_64bit_elem_integers.e \
 	meshes/good_32bit_elem_integers.e \
 	meshes/mesh_assign_test_mesh.xda \
+	meshes/mesh_with_low_order_edge_blocks.e \
 	meshes/shark_tooth_tri6.xda.gz meshes/tetgen_one_tet10.ele \
 	meshes/tetgen_one_tet10.node fparser/autodiff.C
 @LIBMESH_ENABLE_FPARSER_TRUE@am__objects_8 = fparser/unit_tests_opt-autodiff.$(OBJEXT)
@@ -1002,6 +1006,7 @@ am__unit_tests_prof_SOURCES_DIST = driver.C libmesh_cppunit.h \
 	meshes/bad_64bit_elem_integers.e \
 	meshes/good_32bit_elem_integers.e \
 	meshes/mesh_assign_test_mesh.xda \
+	meshes/mesh_with_low_order_edge_blocks.e \
 	meshes/shark_tooth_tri6.xda.gz meshes/tetgen_one_tet10.ele \
 	meshes/tetgen_one_tet10.node fparser/autodiff.C
 @LIBMESH_ENABLE_FPARSER_TRUE@am__objects_10 = fparser/unit_tests_prof-autodiff.$(OBJEXT)
@@ -2168,6 +2173,7 @@ data = meshes/1_quad.bxt.gz \
        meshes/bad_64bit_elem_integers.e \
        meshes/good_32bit_elem_integers.e \
        meshes/mesh_assign_test_mesh.xda \
+       meshes/mesh_with_low_order_edge_blocks.e \
        meshes/shark_tooth_tri6.xda.gz \
        meshes/tetgen_one_tet10.ele \
        meshes/tetgen_one_tet10.node

tests/mesh/mesh_input.C

@@ -105,6 +105,7 @@ public:
   CPPUNIT_TEST( testExodusReadHeader );
 #if LIBMESH_DIM > 2
   CPPUNIT_TEST( testExodusIGASidesets );
+  CPPUNIT_TEST( testLowOrderEdgeBlocks );
 #endif
 #ifndef LIBMESH_USE_COMPLEX_NUMBERS
   CPPUNIT_TEST( testExodusCopyElementVectorDistributed );
@@ -230,6 +231,44 @@ public:
     }
   }
 
+  void testLowOrderEdgeBlocks ()
+  {
+    LOG_UNIT_TEST;
+
+    Mesh mesh(*TestCommWorld);
+    ExodusII_IO exii(mesh);
+
+    if (mesh.processor_id() == 0)
+      exii.read("meshes/mesh_with_low_order_edge_blocks.e");
+
+    MeshCommunication().broadcast(mesh);
+    mesh.prepare_for_use();
+
+    // Check that we see the boundary ids we expect
+    BoundaryInfo & bi = mesh.get_boundary_info();
+
+    // On a ReplicatedMesh, check that the number of edge boundary
+    // conditions is as expected. The real test is that we can read
+    // this file in at all. Prior to the changes in #3491, the Exodus
+    // reader threw an exception while trying to read this mesh.
+    if (mesh.is_serial())
+    {
+      // Mesh has 15 boundary ids total (including edge and side ids).
+      CPPUNIT_ASSERT_EQUAL(static_cast<std::size_t>(15), bi.n_boundary_ids());
+
+      // We can binary_search() the build_edge_list() which is sorted
+      // in lexicographical order before it's returned.
+      auto edge_list = bi.build_edge_list();
+
+      // Search for some tuples we expect to be present
+      CPPUNIT_ASSERT(std::binary_search(edge_list.begin(), edge_list.end(), std::make_tuple(4, 1, 8007)));
+      CPPUNIT_ASSERT(std::binary_search(edge_list.begin(), edge_list.end(), std::make_tuple(10, 6, 8001)));
+
+      // And make sure we don't have entries we shouldn't have
+      CPPUNIT_ASSERT(!std::binary_search(edge_list.begin(), edge_list.end(), std::make_tuple(1, 8, 8009)));
+      CPPUNIT_ASSERT(!std::binary_search(edge_list.begin(), edge_list.end(), std::make_tuple(2, 10, 9011)));
+    }
+  }
 
   void testExodusIGASidesets ()
   {