Merge pull request #735 from sebproell/exodus-numbering-options

Use the same internal numbering for direct Exodus input

Author: sebproell

src/core/io/src/4C_io_exodus.cpp

@@ -107,7 +107,8 @@ namespace
 /*----------------------------------------------------------------------*
  |  ctor (public)                                              maf 12/07|
  *----------------------------------------------------------------------*/
-Core::IO::Exodus::Mesh::Mesh(const std::string exofilename)
+Core::IO::Exodus::Mesh::Mesh(std::filesystem::path exodus_file, MeshParameters mesh_data)
+    : mesh_data_(mesh_data)
 {
   int error;
   int CPU_word_size, IO_word_size;
@@ -118,8 +119,8 @@ Core::IO::Exodus::Mesh::Mesh(const std::string exofilename)
 
   // open EXODUS II file
   int exo_handle =
-      ex_open(exofilename.c_str(), EX_READ, &CPU_word_size, &IO_word_size, &exoversion);
-  if (exo_handle <= 0) FOUR_C_THROW("Error while opening EXODUS II file {}", exofilename);
+      ex_open(exodus_file.c_str(), EX_READ, &CPU_word_size, &IO_word_size, &exoversion);
+  if (exo_handle <= 0) FOUR_C_THROW("Error while opening EXODUS II file {}", exodus_file.string());
 
   // read database parameters
   int num_elem_blk, num_node_sets, num_side_sets, num_nodes;
@@ -138,10 +139,12 @@ Core::IO::Exodus::Mesh::Mesh(const std::string exofilename)
     error = ex_get_coord(exo_handle, x.data(), y.data(), z.data());
     if (error != 0) FOUR_C_THROW("exo error returned");
 
+    FOUR_C_ASSERT_ALWAYS(
+        mesh_data_.node_start_id >= 0, "Node start id must be greater than or equal to 0");
+
     for (int i = 0; i < num_nodes; ++i)
     {
-      // Node numbering starts at one in 4C
-      nodes_[i + 1] = {x[i], y[i], z[i]};
+      nodes_[i + mesh_data_.node_start_id] = {x[i], y[i], z[i]};
     }
   }
 
@@ -175,13 +178,16 @@ Core::IO::Exodus::Mesh::Mesh(const std::string exofilename)
       error = ex_get_conn(exo_handle, EX_ELEM_BLOCK, ebids[i], allconn.data(), nullptr, nullptr);
       if (error != 0) FOUR_C_THROW("exo error returned");
       std::map<int, std::vector<int>> eleconn;
+
+      // Compare the desired start ID to Exodus' one-based indexing to get the offset.
+      const int node_offset = mesh_data_.node_start_id - 1;
       for (int j = 0; j < num_el_in_blk; ++j)
       {
         std::vector<int>& actconn = eleconn[j];
         actconn.reserve(num_nod_per_elem);
         for (int k = 0; k < num_nod_per_elem; ++k)
         {
-          actconn.push_back(allconn[k + j * num_nod_per_elem]);
+          actconn.push_back(allconn[k + j * num_nod_per_elem] + node_offset);
         }
         reorder_nodes_exodus_to_four_c(actconn, cell_type_from_exodus_string(ele_type));
       }
@@ -217,10 +223,12 @@ Core::IO::Exodus::Mesh::Mesh(const std::string exofilename)
       else if (error < 0)
         FOUR_C_THROW("error reading node set");
       std::set<int> nodes_in_set;
-      for (int j = 0; j < num_nodes_in_set; ++j) nodes_in_set.insert(node_set_node_list[j]);
+      // Compare the desired start ID to Exodus' one-based indexing to get the offset.
+      const int node_offset = mesh_data_.node_start_id - 1;
+      for (int j = 0; j < num_nodes_in_set; ++j)
+        nodes_in_set.insert(node_set_node_list[j] + node_offset);
       NodeSet actNodeSet(nodes_in_set, nodesetname);
 
-      // Add this NodeSet into Mesh map (here prelim due to pro names)
       node_sets_.insert(std::pair<int, NodeSet>(npropID[i], actNodeSet));
     }
   }  // end of nodeset section
@@ -362,48 +370,6 @@ void Core::IO::Exodus::Mesh::set_nsd(const int nsd)
 
   four_c_dim_ = nsd;
 }
-/*----------------------------------------------------------------------*/
-/*----------------------------------------------------------------------*/
-
-
-/*----------------------------------------------------------------------*/
-/*----------------------------------------------------------------------*/
-std::vector<double> Core::IO::Exodus::Mesh::normal(
-    const int head1, const int origin, const int head2) const
-{
-  std::vector<double> normal(3);
-  std::vector<double> h1 = get_node(head1);
-  std::vector<double> h2 = get_node(head2);
-  std::vector<double> o = get_node(origin);
-
-  normal[0] = ((h1[1] - o[1]) * (h2[2] - o[2]) - (h1[2] - o[2]) * (h2[1] - o[1]));
-  normal[1] = -((h1[0] - o[0]) * (h2[2] - o[2]) - (h1[2] - o[2]) * (h2[0] - o[0]));
-  normal[2] = ((h1[0] - o[0]) * (h2[1] - o[1]) - (h1[1] - o[1]) * (h2[0] - o[0]));
-
-  double length = sqrt(normal[0] * normal[0] + normal[1] * normal[1] + normal[2] * normal[2]);
-  normal[0] = normal[0] / length;
-  normal[1] = normal[1] / length;
-  normal[2] = normal[2] / length;
-
-  return normal;
-}
-
-/*----------------------------------------------------------------------*/
-/*----------------------------------------------------------------------*/
-std::vector<double> Core::IO::Exodus::Mesh::node_vec(const int tail, const int head) const
-{
-  std::vector<double> nv(3);
-  std::vector<double> t = get_node(tail);
-  std::vector<double> h = get_node(head);
-  nv[0] = h[0] - t[0];
-  nv[1] = h[1] - t[1];
-  nv[2] = h[2] - t[2];
-  double length = sqrt(nv[0] * nv[0] + nv[1] * nv[1] + nv[2] * nv[2]);
-  nv[0] = nv[0] / length;
-  nv[1] = nv[1] / length;
-  nv[2] = nv[2] / length;
-  return nv;
-}
 
 
 /*----------------------------------------------------------------------*/

src/core/io/src/4C_io_exodus.hpp

@@ -12,6 +12,7 @@
 
 #include "4C_fem_general_element.hpp"
 
+#include <filesystem>
 #include <set>
 #include <string>
 #include <vector>
@@ -25,14 +26,31 @@ namespace Core::IO::Exodus
   class NodeSet;
   class SideSet;
 
+  /**
+   * Additional parameters that are used in the constructor.
+   */
+  struct MeshParameters
+  {
+    /**
+     * The ID of the first node in the mesh. This defaults to 1, since this is the default in
+     * the Exodus II mesh format.
+     */
+    int node_start_id{1};
+  };
+
   /**
    * A class that stores the mesh information read from an Exodus file.
    */
   class Mesh
   {
    public:
-    //! constructor
-    Mesh(std::string exofilename);
+    /**
+     * @brief Read a mesh from an Exodus file.
+     *
+     * Read the data from the Exodus file and store it in the class. The optional @p
+     * mesh_data can be used to set options documented in the MeshData struct.
+     */
+    Mesh(std::filesystem::path exodus_file, MeshParameters mesh_data = {});
 
     //! Print mesh info
     void print(std::ostream& os, bool verbose = false) const;
@@ -76,12 +94,6 @@ namespace Core::IO::Exodus
     //! Get one SideSet
     SideSet get_side_set(const int id) const;
 
-    //! Get edge Normal at node
-    std::vector<double> normal(const int head1, const int origin, const int head2) const;
-
-    //! Get normalized Vector between 2 nodes
-    std::vector<double> node_vec(const int tail, const int head) const;
-
     //! Get Node map
     const std::map<int, std::vector<double>>& get_nodes() const { return nodes_; }
 
@@ -95,6 +107,8 @@ namespace Core::IO::Exodus
     void set_nsd(const int nsd);
 
    private:
+    MeshParameters mesh_data_;
+
     std::map<int, std::vector<double>> nodes_;
 
     std::map<int, ElementBlock> element_blocks_;

src/core/io/src/4C_io_meshreader.cpp

@@ -822,8 +822,11 @@ namespace
 
       const auto& geometry_data = data.group(exodus_reader.section_name);
       const auto& exodus_file = geometry_data.get<std::filesystem::path>("FILE");
-      exodus_reader.mesh_on_rank_zero =
-          std::make_unique<Core::IO::Exodus::Mesh>(exodus_file.string());
+      exodus_reader.mesh_on_rank_zero = std::make_unique<Core::IO::Exodus::Mesh>(
+          exodus_file.string(), Core::IO::Exodus::MeshParameters{
+                                    // We internally depend on node numbers starting at 0.
+                                    .node_start_id = 0,
+                                });
       const auto& mesh = *exodus_reader.mesh_on_rank_zero;
 
       // Initial implementation:

src/core/io/tests/4C_io_exodus_test.cpp

@@ -29,4 +29,13 @@ namespace
 
     mesh.print(std::cout, true);
   }
+
+  TEST(Exodus, NodeOffset)
+  {
+    Core::IO::Exodus::Mesh mesh(TESTING::get_support_file_path("test_files/exodus/cube.exo"),
+        Core::IO::Exodus::MeshParameters{.node_start_id = 100});
+    EXPECT_EQ(mesh.get_node(100), (std::vector<double>{-5.0, 0.0, 0.0}));
+    EXPECT_EQ(mesh.get_element_block(1).get_ele_nodes(1),
+        (std::vector<int>{108, 100, 103, 109, 110, 104, 107, 111}));
+  }
 }  // namespace

tests/input_files/contact3D_exodus_in.4C.yaml

@@ -82,13 +82,13 @@ DESIGN SURF DIRICH CONDITIONS:
 RESULT DESCRIPTION:
   - STRUCTURE:
       DIS: structure
-      NODE: 1270
+      NODE: 1269
       QUANTITY: dispx
       VALUE: -8.31974725251016095e-03
       TOLERANCE: 1e-8
   - STRUCTURE:
       DIS: structure
-      NODE: 1372
+      NODE: 1371
       QUANTITY: dispx
       VALUE: -8.36015173005985428e-03
       TOLERANCE: 1e-8

tests/input_files/fluid_exodus_in.4C.yaml

@@ -65,73 +65,73 @@ DESIGN POINT DIRICH CONDITIONS:
 RESULT DESCRIPTION:
   - FLUID:
       DIS: "fluid"
-      NODE: 2912
+      NODE: 2911
       QUANTITY: "velx"
       VALUE: 0
       TOLERANCE: 1e-08
   - FLUID:
       DIS: "fluid"
-      NODE: 2912
+      NODE: 2911
       QUANTITY: "vely"
       VALUE: 0
       TOLERANCE: 1e-08
   - FLUID:
       DIS: "fluid"
-      NODE: 2912
+      NODE: 2911
       QUANTITY: "pressure"
       VALUE: 2.316272235592883
       TOLERANCE: 1e-10
   - FLUID:
       DIS: "fluid"
-      NODE: 4034
+      NODE: 4033
       QUANTITY: "velx"
       VALUE: 0
       TOLERANCE: 1e-08
   - FLUID:
       DIS: "fluid"
-      NODE: 4034
+      NODE: 4033
       QUANTITY: "vely"
       VALUE: 0
       TOLERANCE: 1e-08
   - FLUID:
       DIS: "fluid"
-      NODE: 4034
+      NODE: 4033
       QUANTITY: "pressure"
       VALUE: 1.5110014506671248
       TOLERANCE: 1e-08
   - FLUID:
       DIS: "fluid"
-      NODE: 6667
+      NODE: 6666
       QUANTITY: "velx"
       VALUE: 0.8140348125434049
       TOLERANCE: 1e-08
   - FLUID:
       DIS: "fluid"
-      NODE: 6667
+      NODE: 6666
       QUANTITY: "vely"
       VALUE: -0.022743581282554975
       TOLERANCE: 1e-10
   - FLUID:
       DIS: "fluid"
-      NODE: 6667
+      NODE: 6666
       QUANTITY: "pressure"
       VALUE: 2.0581385928541125
       TOLERANCE: 1e-10
   - FLUID:
       DIS: "fluid"
-      NODE: 6998
+      NODE: 6997
       QUANTITY: "velx"
       VALUE: 0.9708173796517694
       TOLERANCE: 1e-08
   - FLUID:
       DIS: "fluid"
-      NODE: 6998
+      NODE: 6997
       QUANTITY: "vely"
       VALUE: -7.977045433472414e-16
       TOLERANCE: 1e-08
   - FLUID:
       DIS: "fluid"
-      NODE: 6998
+      NODE: 6997
       QUANTITY: "pressure"
       VALUE: 1.8450141345624358
       TOLERANCE: 1e-10