Add functionality to convert exodus files to dictionaries

src/cubitpy/cubit_to_fourc_input.py

@@ -33,46 +33,12 @@
 from fourcipp.fourc_input import FourCInput
 
 from cubitpy.conf import GeometryType, cupy
+from cubitpy.exodus_utility import get_exo_info
 
 if TYPE_CHECKING:
     from cubitpy.cubitpy import CubitPy
 
 
-def get_exo_info(exo, entry_type) -> tuple[dict, dict]:
-    """Build mappings between Exodus IDs and Cubit IDs for blocks or
-    nodesets."""
-
-    if entry_type == "block":
-        exo_identifier = "eb"
-    elif entry_type == "nodeset":
-        exo_identifier = "ns"
-    else:
-        raise ValueError(f"Invalid entry type: {entry_type}")
-
-    if exo_identifier + "_names" not in exo.variables.keys():
-        return {}, {}
-
-    # List of explicitly given names
-    names = []
-    for line in exo.variables[exo_identifier + "_names"]:
-        name: str | None = str(netCDF4.chartostring(line))
-        if name == "":
-            name = None
-        names.append(name)
-
-    # Get information of all entries of the given type
-    cubit_id_to_info = {}
-    exo_id_to_info = {}
-    for exo_id, cubit_id in enumerate(
-        exo.variables[exo_identifier + "_prop1"][:].tolist()
-    ):
-        info = {"cubit_id": cubit_id, "exo_id": exo_id, "name": names[exo_id]}
-        cubit_id_to_info[cubit_id] = info.copy()
-        exo_id_to_info[exo_id] = info.copy()
-
-    return cubit_id_to_info, exo_id_to_info
-
-
 def add_node_sets_external_geometry(cubit: CubitPy, input_file: FourCInput) -> None:
     """Add a reference to the node sets contained in the cubit session/exo file
     to the yaml file."""

src/cubitpy/exodus_utility.py

@@ -0,0 +1,97 @@
+# The MIT License (MIT)
+#
+# Copyright (c) 2018-2026 CubitPy Authors
+#
+# Permission is hereby granted, free of charge, to any person obtaining a copy
+# of this software and associated documentation files (the "Software"), to deal
+# in the Software without restriction, including without limitation the rights
+# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+# copies of the Software, and to permit persons to whom the Software is
+# furnished to do so, subject to the following conditions:
+#
+# The above copyright notice and this permission notice shall be included in
+# all copies or substantial portions of the Software.
+#
+# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
+# THE SOFTWARE.
+"""Utility functions interacting with the exodus file format."""
+
+from pathlib import Path
+
+import netCDF4
+import numpy as np
+
+
+def get_exo_info(exo, entry_type) -> tuple[dict, dict]:
+    """Build mappings between Exodus IDs and Cubit IDs for blocks or
+    nodesets."""
+
+    if entry_type == "block":
+        exo_identifier = "eb"
+    elif entry_type == "nodeset":
+        exo_identifier = "ns"
+    else:
+        raise ValueError(f"Invalid entry type: {entry_type}")
+
+    if exo_identifier + "_names" not in exo.variables.keys():
+        return {}, {}
+
+    # List of explicitly given names
+    names = []
+    for line in exo.variables[exo_identifier + "_names"]:
+        name: str | None = str(netCDF4.chartostring(line))
+        if name == "":
+            name = None
+        names.append(name)
+
+    # Get information of all entries of the given type
+    cubit_id_to_info = {}
+    exo_id_to_info = {}
+    for exo_id, cubit_id in enumerate(
+        exo.variables[exo_identifier + "_prop1"][:].tolist()
+    ):
+        info = {"cubit_id": cubit_id, "exo_id": exo_id, "name": names[exo_id]}
+        cubit_id_to_info[cubit_id] = info.copy()
+        exo_id_to_info[exo_id] = info.copy()
+
+    return cubit_id_to_info, exo_id_to_info
+
+
+def convert_exodus_to_dict(exo_path: Path) -> dict:
+    """Load an exodus file from disk and convert it to a dictionary containing
+    all the relevant information.
+
+    This function can be used in testing to compare exodus files with
+    each other.
+    """
+    exo_data = {}
+    with netCDF4.Dataset(exo_path) as exo:
+        exo_data["coordinates"] = (
+            np.array(
+                [exo.variables["coord" + dim][:] for dim in ["x", "y", "z"]],
+            )
+            .transpose()
+            .tolist()
+        )
+
+        _, exo_block_id_to_info = get_exo_info(exo, "block")
+        exo_data["exo_block_id_to_info"] = exo_block_id_to_info
+        for exo_id in exo_block_id_to_info.keys():
+            connectivity_name = f"connect{exo_id + 1}"
+            connectivity = np.array(exo.variables[connectivity_name][:]).tolist()
+            exo_data[connectivity_name] = connectivity
+
+        _, exo_node_set_id_to_info = get_exo_info(exo, "nodeset")
+        exo_data["exo_node_set_id_to_info"] = exo_node_set_id_to_info
+        for exo_id in exo_node_set_id_to_info.keys():
+            node_set_name = f"node_ns{exo_id + 1}"
+            unique_node_set_ids = np.unique(np.array(exo.variables[node_set_name][:]))
+            # The correct data type here would be `set`. However, this can not be
+            # serialized by FourCIPP, so we return an ordered list of the unique IDs here.
+            exo_data[node_set_name] = unique_node_set_ids.tolist()
+    return exo_data

tests/test_cubitpy.py

@@ -42,6 +42,7 @@
     string_to_node_set_info,
 )
 from cubitpy.cubitpy import CubitPy
+from cubitpy.exodus_utility import convert_exodus_to_dict
 from cubitpy.geometry_creation_functions import (
     create_brick_by_corner_points,
     create_parametric_surface,
@@ -126,17 +127,23 @@ def compare_yaml(
     ref_file = Path(testing_input) / (compare_name + ".4C.yaml")
     out_file = Path(testing_temp) / (compare_name + ".4C.yaml")
 
-    if mesh_in_exo:
-        # dump the input script with the mesh in exodus format
-        cubit.dump(out_file, mesh_in_exo=True)
-        # make sure the directory also contains the exo mesh
-        out_file_stem = out_file.with_name(out_file.name.removesuffix(".4C.yaml"))
-        assert out_file_stem.with_suffix(".exo").exists()
-    else:
-        cubit.dump(out_file)
+    cubit.dump(out_file, mesh_in_exo=mesh_in_exo)
+
+    def get_input_file_with_data(input_file_path):
+        """Load the input file into a FourCIPP object and if the mesh is given
+        in exodus format, add this as a data structure for comparison."""
+        input_file = FourCInput.from_4C_yaml(input_file_path)
+        if mesh_in_exo:
+            exo_path = input_file_path.parent / (compare_name + ".exo")
+            exo_data = convert_exodus_to_dict(exo_path)
+            section_name = "STRUCTURE GEOMETRY"
+            structure_section = input_file.pop(section_name)
+            structure_section["FILE"] = exo_data
+            input_file.combine_sections({section_name: structure_section})
+        return input_file
 
-    ref_input_file = FourCInput.from_4C_yaml(ref_file)
-    out_input_file = FourCInput.from_4C_yaml(out_file)
+    ref_input_file = get_input_file_with_data(ref_file)
+    out_input_file = get_input_file_with_data(out_file)
 
     try:
         files_are_equal = ref_input_file.compare(
@@ -2503,3 +2510,54 @@ def test_cubit_warnings_and_errors():
         match="ERROR: All dimensions must be nonzero and positive. Entered values are:",
     ):
         cubit.cmd("brick x -10")
+
+
+def test_exodus_to_dict():
+    """Create a simple model and check the convert exodus to dict
+    functionality."""
+
+    # Create the brick with a single solid element
+    cubit = CubitPy()
+    create_brick(
+        cubit,
+        1,
+        2,
+        3,
+        mesh_interval=[1, 1, 1],
+        element_type=cupy.element_type.hex8,
+        name="cube",
+    )
+
+    # Add node sets to the top surface and to the whole volume of the brick, to make
+    # the exodus file more general.
+    cubit.add_node_set(cubit.group(add_value="add surface 2"), name="top")
+    cubit.add_node_set(cubit.group(add_value="add volume 1"), name="all")
+
+    input_file_path = Path(testing_temp) / "exo_to_dict.4C.yaml"
+    exo_path = Path(testing_temp) / "exo_to_dict.exo"
+    cubit.dump(input_file_path, mesh_in_exo=True)
+
+    exo_dict = convert_exodus_to_dict(exo_path)
+
+    reference_exo_dict = {
+        "coordinates": [
+            [-0.5, -1.0, 1.5],
+            [-0.5, -1.0, -1.5],
+            [-0.5, 1.0, -1.5],
+            [-0.5, 1.0, 1.5],
+            [0.5, -1.0, 1.5],
+            [0.5, -1.0, -1.5],
+            [0.5, 1.0, -1.5],
+            [0.5, 1.0, 1.5],
+        ],
+        "exo_block_id_to_info": {0: {"cubit_id": 1, "exo_id": 0, "name": "cube"}},
+        "connect1": [[1, 2, 3, 4, 5, 6, 7, 8]],
+        "exo_node_set_id_to_info": {
+            0: {"cubit_id": 1, "exo_id": 0, "name": "cpy_1_s_top"},
+            1: {"cubit_id": 2, "exo_id": 1, "name": "cpy_2_v_all"},
+        },
+        "node_ns1": [2, 3, 6, 7],
+        "node_ns2": [1, 2, 3, 4, 5, 6, 7, 8],
+    }
+
+    compare_nested_dicts_or_lists(exo_dict, reference_exo_dict)