Merge pull request #3084 from roystgnr/iga_multiblock_fix

Exodus IGA fix for multiblock files

Author: roystgnr

src/mesh/exodusII_io.C

@@ -46,6 +46,28 @@
 #include <map>
 #include <sstream>
 
+namespace
+{
+  using namespace libMesh;
+
+  const std::vector<Real> & bex_constraint_vec(std::size_t i,
+                                               const ExodusII_IO_Helper & helper)
+  {
+    std::size_t vec_offset = 0;
+    for (auto block_num : index_range(helper.bex_dense_constraint_vecs))
+    {
+      const auto & vecblock = helper.bex_dense_constraint_vecs[block_num];
+      libmesh_assert_greater_equal(i, vec_offset);
+      if (i - vec_offset < vecblock.size())
+        return vecblock[i - vec_offset];
+
+      vec_offset += vecblock.size();
+    }
+
+    libmesh_error_msg("Requested BEX coefficient vector " << i << " not found");
+  }
+}
+
 namespace libMesh
 {
 
@@ -456,9 +478,7 @@ void ExodusII_IO::read (const std::string & fname)
                   const auto & my_constraint_rows = exio_helper->bex_cv_conn[elem_num];
                   const unsigned long elem_coef_vec_index =
                     my_constraint_rows[spline_node_index] - 1; // Exodus isn't 0-based
-                  const auto & my_vec =
-                    libmesh_vector_at(exio_helper->bex_dense_constraint_vecs[i],
-                                      elem_coef_vec_index);
+                  const auto & my_vec = bex_constraint_vec(elem_coef_vec_index, *exio_helper);
                   for (auto elem_node_index :
                        make_range(elem->n_nodes()))
                     {
@@ -484,9 +504,11 @@ void ExodusII_IO::read (const std::string & fname)
                       const unsigned long elem_coef_vec_index =
                         my_constraint_rows[spline_node_index] - 1; // Exodus isn't 0-based
 
+                      auto & coef_vec =
+                        bex_constraint_vec(elem_coef_vec_index, *exio_helper);
+
                       const Real coef =
-                        libmesh_vector_at(exio_helper->bex_dense_constraint_vecs[i],
-                                          elem_coef_vec_index)[elem_node_index];
+                        libmesh_vector_at(coef_vec, elem_node_index);
 
                       const int gi = (elem_num)*exio_helper->num_nodes_per_elem +
                         spline_node_index;
@@ -522,9 +544,10 @@ void ExodusII_IO::read (const std::string & fname)
 
                           const Node & spline_node = mesh.node_ref(libmesh_node_id);
 
+                          auto & coef_vec =
+                            bex_constraint_vec(elem_coef_vec_index, *exio_helper);
                           const Real coef =
-                            libmesh_vector_at(exio_helper->bex_dense_constraint_vecs[i],
-                                              elem_coef_vec_index)[elem_node_index];
+                            libmesh_vector_at(coef_vec, elem_node_index);
 
                           // We don't need to store 0 entries;
                           // constraint_rows is a sparse structure.

src/mesh/exodusII_io_helper.C

@@ -909,6 +909,9 @@ void ExodusII_IO_Helper::read_bex_cv_blocks()
               exII::ex_get_attribute(ex_id, &attr);
 
               bex_num_dense_cv_blocks = value_count / 2;
+
+              libmesh_error_msg_if(bex_num_dense_cv_blocks > 1,
+                                   "Found more than 1 dense bex CV block; unsure how to handle that");
             }
         }
 

tests/mesh/mesh_input.C

@@ -67,6 +67,7 @@ public:
 #endif // !LIBMESH_USE_COMPLEX_NUMBERS
 
   CPPUNIT_TEST( testExodusFileMappingsPlateWithHole);
+  CPPUNIT_TEST( testExodusFileMappingsTwoBlocks);
   CPPUNIT_TEST( testExodusFileMappingsCyl3d);
 #endif // LIBMESH_HAVE_EXODUS_API
 
@@ -712,8 +713,13 @@ public:
   {
     auto locator = mesh.sub_point_locator();
 
-    const std::set<subdomain_id_type> manifold_subdomain { 0 };
-    const std::set<subdomain_id_type> nodeelem_subdomain { 1 };
+    // A multiblock Exodus file will have a subdomain for each block;
+    // we stick the nodeelem subdomain at the end.
+    const subdomain_id_type nodeelems = mesh.n_subdomains()-1;
+    const std::set<subdomain_id_type> nodeelem_subdomain { nodeelems };
+    std::set<subdomain_id_type> manifold_subdomain;
+    for (auto i : make_range(nodeelems))
+      manifold_subdomain.insert(i);
 
     for (auto & elem : mesh.element_ptr_range())
       {
@@ -999,6 +1005,14 @@ public:
                             1.791640772215248, 1.413679237529376});
   }
 
+  void testExodusFileMappingsTwoBlocks ()
+  {
+    testExodusFileMappings("meshes/two_quads_two_blocks.e",
+    // Regression values for sin_x_plus_cos_y
+                           {2.03496953073072, 1.97996853164955,
+                            1.18462134113435, 1.03085301158959});
+  }
+
   void testExodusFileMappingsCyl3d ()
   {
     testExodusFileMappings("meshes/PressurizedCyl3d_Patch1_8Elem.e",