Add yaml with exo export

src/cubitpy/cubit_to_fourc_input.py

@@ -26,13 +26,14 @@
 
 import netCDF4
 import numpy as np
+from fourcipp.legacy_io.inline_dat import to_dat_string
 
 from cubitpy.conf import cupy
 
 
-def add_node_sets(cubit, exo, input_file):
+def add_node_sets(cubit, exo, input_file, write_topology_information=True):
     """Add the node sets contained in the cubit session/exo file to the
-    dat_lines."""
+    yaml file."""
 
     # If there are no node sets we can return immediately
     if len(cubit.node_sets) == 0:
@@ -69,28 +70,105 @@ def add_node_sets(cubit, exo, input_file):
             input_file[bc_section] = []
         bc_description["E"] = len(node_sets[geometry_type])
 
+        if not write_topology_information:
+            # when working with external .exo meshes, we do not write the
+            # topology information for the node sets explicitly, since 4C will
+            # deduce them based on the node set ids, when reading the .exo file
+            bc_description["ENTITY_TYPE"] = "node_set_id"
+
         input_file[bc_section].append(bc_description)
 
-    name_geometry_tuple = [
-        [cupy.geometry.vertex, "DNODE-NODE TOPOLOGY", "DNODE"],
-        [cupy.geometry.curve, "DLINE-NODE TOPOLOGY", "DLINE"],
-        [cupy.geometry.surface, "DSURF-NODE TOPOLOGY", "DSURFACE"],
-        [cupy.geometry.volume, "DVOL-NODE TOPOLOGY", "DVOL"],
-    ]
-    for geo, section_name, set_label in name_geometry_tuple:
-        if len(node_sets[geo]) > 0:
-            input_file[section_name] = []
-            for i_set, node_set in enumerate(node_sets[geo]):
-                node_set.sort()
-                for i_node in node_set:
-                    input_file[section_name].append(
-                        {
-                            "type": "NODE",
-                            "node_id": i_node,
-                            "d_type": set_label,
-                            "d_id": i_set + 1,
-                        }
-                    )
+    if write_topology_information:
+        # this is the default case: when the mesh is supposed to be contained
+        # in the .yaml file, we have to write the topology information of the
+        # node sets
+        name_geometry_tuple = [
+            [cupy.geometry.vertex, "DNODE-NODE TOPOLOGY", "DNODE"],
+            [cupy.geometry.curve, "DLINE-NODE TOPOLOGY", "DLINE"],
+            [cupy.geometry.surface, "DSURF-NODE TOPOLOGY", "DSURFACE"],
+            [cupy.geometry.volume, "DVOL-NODE TOPOLOGY", "DVOL"],
+        ]
+        for geo, section_name, set_label in name_geometry_tuple:
+            if len(node_sets[geo]) > 0:
+                input_file[section_name] = []
+                for i_set, node_set in enumerate(node_sets[geo]):
+                    node_set.sort()
+                    for i_node in node_set:
+                        input_file[section_name].append(
+                            {
+                                "type": "NODE",
+                                "node_id": i_node,
+                                "d_type": set_label,
+                                "d_id": i_set + 1,
+                            }
+                        )
+
+
+def add_exodus_geometry_section(cubit, rel_exo_file_path):
+    """Add the problem specific geometry section to the input file required to directly read the mesh from an exodus file.
+
+    This section contains information about all element blocks as well as the
+    path to the exo file that contains the mesh.
+
+    Args
+    ----
+    cubit: CubitPy
+        The python object for managing the current Cubit session.
+    rel_exo_file_path: str
+        The relative path (as seen from the yaml input file) to the exodus
+        file that contains the mesh.
+    """
+
+    # map the problem type to the keyword for the problem specific geometry
+    # section
+    problem_to_geometry_dict = {
+        "Structure": "STRUCTURE GEOMETRY",
+        "Fluid": "FLUID GEOMETRY",
+    }
+
+    # retrieve the problem type
+    problem_type = cubit.fourc_input["PROBLEM TYPE"]["PROBLEMTYPE"]
+    # check if the problem type is supported
+    problem_key = problem_to_geometry_dict.get(problem_type)
+    if problem_key is None:
+        raise ValueError(
+            f"I do not know how to generate a 'GEOMETRY' section for problem type '{problem_type}'. "
+            "Currently supported types are: "
+            f"{', '.join(problem_to_geometry_dict.keys())}."
+        )
+
+    # create the dictionary for the geometry section
+    geometry_section_dict = {}
+
+    # add the relative path to the exodus file
+    geometry_section_dict["FILE"] = rel_exo_file_path
+    # always show a detailed summary
+    geometry_section_dict["SHOW_INFO"] = "detailed_summary"
+
+    # add the element blocks
+    geometry_section_dict["ELEMENT_BLOCKS"] = []
+    # iterate over all blocks
+    block_ids = cubit.cubit.get_block_id_list()
+    for idx, bid in enumerate(block_ids):
+        # retrieve the data associated with the block
+        data = cubit.blocks[idx]
+        # retrieve the fourc name for the element
+        four_c_element_name = data[0].get_four_c_name()
+        # convert the material data from dict to string because 4C currently does not support a dict here
+        element_data_string = " ".join(
+            f"{key} {value}" for key, value in data[1].items()
+        )
+        # add block id, fourc element name and element data string to the element block dictionary
+        element_block_dict = {
+            "ID": bid,
+            "ELEMENT_NAME": four_c_element_name,
+            "ELEMENT_DATA": element_data_string,
+        }
+        # append the dictionary with the element block information to the element block list
+        geometry_section_dict["ELEMENT_BLOCKS"].append(element_block_dict)
+
+    # at long last, append the geometry section to the input file
+    cubit.fourc_input[problem_key] = geometry_section_dict
 
 
 def get_element_connectivity_list(connectivity):

src/cubitpy/cubitpy.py

@@ -26,11 +26,16 @@
 import time
 import warnings
 
+import netCDF4
 from fourcipp.fourc_input import FourCInput
 
 from cubitpy.conf import cupy
 from cubitpy.cubit_group import CubitGroup
-from cubitpy.cubit_to_fourc_input import get_input_file_with_mesh
+from cubitpy.cubit_to_fourc_input import (
+    add_node_sets,
+    add_problem_geometry_section,
+    get_input_file_with_mesh,
+)
 from cubitpy.cubit_wrapper.cubit_wrapper_host import CubitConnect
 
 
@@ -329,6 +334,38 @@ def export_exo(self, path):
         """Export the mesh."""
         self.cubit.cmd('export mesh "{}" dimension 3 overwrite'.format(path))
 
+    def dump_with_exo_mesh(self, path_stem):
+        """Create the 4C yaml input file and export the mesh in exo format.
+
+        The path_stem refers to the path where the input file and the mesh
+        should be saved, plus the prefix for the file names, but no endings.
+        The file endings will be added automatically, i.e., ``.4C.yaml`` for
+        the yaml input file and ``.exo`` for the mesh file.
+
+        Args
+        ----
+        path_stem: str
+            Path stem where the input file and the mesh will be saved.
+        """
+
+        # Check if output path exists
+        dat_dir = os.path.dirname(os.path.abspath(path_stem))
+        if not os.path.exists(dat_dir):
+            raise ValueError("Path {} does not exist!".format(dat_dir))
+
+        # Create the exodus file
+        exo_path = f"{path_stem}.exo"
+        self.export_exo(exo_path)
+        # parse the exodus file
+        exo = netCDF4.Dataset(exo_path)
+        # Add the node sets
+        add_node_sets(self, exo, self.fourc_input, write_topology_information=False)
+        # Add the problem geometry section
+        rel_exo_path = os.path.relpath(exo_path, start=dat_dir)
+        add_problem_geometry_section(self, rel_exo_path)
+        # Export the input file (now with nodesets) in YAML format
+        self.fourc_input.dump(f"{path_stem}.4C.yaml")
+
     def dump(self, yaml_path):
         """Create the yaml file and save it in under provided yaml_path.
 

tests/input-files-ref/test_yaml_with_exo_export.4C.yaml

@@ -0,0 +1,54 @@
+PROBLEM SIZE:
+  DIM: 3
+PROBLEM TYPE:
+  PROBLEMTYPE: "Structure"
+DESIGN SURF MORTAR CONTACT CONDITIONS 3D:
+  - InterfaceID: 1
+    Side: "Slave"
+    E: 1
+    ENTITY_TYPE: "node_set_id"
+  - InterfaceID: 1
+    Side: "Master"
+    E: 2
+    ENTITY_TYPE: "node_set_id"
+DESIGN SURF DIRICH CONDITIONS:
+  - NUMDOF: 3
+    ONOFF:
+      - 1
+      - 1
+      - 1
+    VAL:
+      - 0
+      - 0
+      - 0
+    FUNCT:
+      - null
+      - null
+      - null
+    E: 3
+    ENTITY_TYPE: "node_set_id"
+  - NUMDOF: 3
+    ONOFF:
+      - 1
+      - 0
+      - 0
+    VAL:
+      - -1.0
+      - 0.0
+      - 0.0
+    FUNCT:
+      - 1
+      - null
+      - null
+    E: 4
+    ENTITY_TYPE: "node_set_id"
+STRUCTURE GEOMETRY:
+  FILE: "test_yaml_with_exo_export.exo"
+  SHOW_INFO: "detailed_summary"
+  ELEMENT_BLOCKS:
+    - ID: 1
+      ELEMENT_NAME: "SOLID"
+      ELEMENT_DATA: "MAT 1 KINEM nonlinear"
+    - ID: 2
+      ELEMENT_NAME: "SOLID"
+      ELEMENT_DATA: "MAT 2 KINEM nonlinear"

tests/test_cubitpy.py

@@ -64,7 +64,13 @@ def check_tmp_dir():
 
 
 def compare_yaml(
-    cubit, *, base_name=None, additional_identifier=None, rtol=1.0e-12, atol=1.0e-12
+    cubit,
+    *,
+    base_name=None,
+    additional_identifier=None,
+    rtol=1.0e-12,
+    atol=1.0e-12,
+    dump_yaml=True,
 ):
     """Write and compare the YAML file from a Cubit object with the reference
     YAML file.
@@ -82,6 +88,12 @@ def compare_yaml(
         Relative tolerance for numerical differences.
     atol: float
         Absolute tolerance for numerical differences.
+    dump_yaml: bool
+        If True, this function will dump the YAML input to the (temporary)
+        output directory. This is the default behavior. If False, no YAML file
+        will be dumped from this function and it is assumed that a YAML file
+        with the corresponding name has already been generated prior to
+        invoking this function.
     """
     # Determine test name
     if base_name is None:
@@ -101,7 +113,8 @@ def compare_yaml(
     # File paths
     ref_file = os.path.join(testing_input, compare_name + ".4C.yaml")
     out_file = os.path.join(testing_temp, compare_name + ".4C.yaml")
-    cubit.dump(out_file)
+    if dump_yaml:
+        cubit.dump(out_file)
 
     ref_input_file = FourCInput.from_4C_yaml(ref_file)
     out_input_file = FourCInput.from_4C_yaml(out_file)
@@ -1759,3 +1772,102 @@ def test_extrude_artery_of_aneurysm():
     assert ref_volume == pytest.approx(
         cubit.get_meshed_volume_or_area("volume", [volume.id()]), 1e-5
     )
+
+
+def test_yaml_with_exo_export():
+    """Test if exporting a yaml file with an exodus mesh works."""
+    # Set up Cubit.
+    cubit = CubitPy()
+
+    # Initialize geometry
+    cubit.cmd("brick x 1 y 1 z 1")
+    cubit.cmd("brick x 5e-1 y 5e-1 z 5e-1")
+    cubit.cmd("move Volume 2 x 75e-2 y 0 z 0")
+    cubit.cmd("volume 1 size {1e-1}")
+    cubit.cmd("volume 2 size {1e-1}")
+
+    # mesh the two geometries
+    cubit.cmd("mesh volume 1")
+    cubit.cmd("mesh volume 2")
+
+    # Assign nodesets, required for boundary conditions
+    cubit.add_node_set(
+        cubit.group(add_value="add surface 6"),
+        name="slave",
+        bc_type=cupy.bc_type.solid_to_solid_contact,
+        bc_description={
+            "InterfaceID": 1,
+            "Side": "Slave",
+        },
+    )
+    cubit.add_node_set(
+        cubit.group(add_value="add surface 10"),
+        name="master",
+        bc_type=cupy.bc_type.solid_to_solid_contact,
+        bc_description={
+            "InterfaceID": 1,
+            "Side": "Master",
+        },
+    )
+    cubit.add_node_set(
+        cubit.group(add_value="add surface 4"),
+        name="wall",
+        bc_type=cupy.bc_type.dirichlet,
+        bc_description={
+            "NUMDOF": 3,
+            "ONOFF": [1, 1, 1],
+            "VAL": [0, 0, 0],
+            "FUNCT": [None, None, None],
+        },
+    )
+    cubit.add_node_set(
+        cubit.group(add_value="add surface 12"),
+        name="pushing",
+        bc_type=cupy.bc_type.dirichlet,
+        bc_description={
+            "NUMDOF": 3,
+            "ONOFF": [1, 0, 0],
+            "VAL": [-1.0, 0.0, 0.0],
+            "FUNCT": [1, None, None],
+        },
+    )
+
+    # Add the element types
+    cubit.add_element_type(
+        cubit.group(add_value="add volume 1"),
+        el_type=cupy.element_type.hex8,
+        material={
+            "MAT": 1,
+        },
+        bc_description={
+            "KINEM": "nonlinear",
+        },
+    )
+    cubit.add_element_type(
+        cubit.group(add_value="add volume 2"),
+        el_type=cupy.element_type.hex8,
+        material={
+            "MAT": 2,
+        },
+        bc_description={
+            "KINEM": "nonlinear",
+        },
+    )
+
+    cubit.fourc_input.combine_sections(
+        {
+            "PROBLEM SIZE": {"DIM": 3},
+            "PROBLEM TYPE": {"PROBLEMTYPE": "Structure"},
+        }
+    )
+
+    # export the yaml input with exo mesh into a temporary directory
+    check_tmp_dir()
+    out_file_stem = os.path.join(testing_temp, "test_yaml_with_exo_export")
+    cubit.dump_with_exo_mesh(out_file_stem)
+
+    # make sure the directory also contains the exo mesh
+    assert os.path.exists(f"{out_file_stem}.exo")
+
+    # Compare the input file created for 4C.
+    compare_yaml(cubit, dump_yaml=False)