Switch to fourcipp output

CONTRIBUTORS

@@ -6,6 +6,7 @@ Ivo Steinbrecher (@isteinbrecher)
 # Contributors (in alphabetical order)
 Janina Datz (@datzj)
 Martin Frank (@knarfnitram)
+Moritz Frey (@m-frey)
 Nora Hagmeyer (@NoraHagmeyer)
 Bishr Maradni (@BishrMaradni)
 Gabriela Loera (@eulovi)

pyproject.toml

@@ -13,7 +13,8 @@ license = "MIT"
 dependencies = [
   "execnet==1.9.0",
   "netCDF4",
-  "numpy"
+  "numpy",
+  "fourcipp @ git+https://github.com/4C-multiphysics/fourcipp.git@main"
 ]
 version = "0.0.1"
 
@@ -25,7 +26,8 @@ Issues = "https://github.com/imcs-compsim/cubitpy/issues/"
 dev = [
   "pre-commit",
   "pytest",
-  "pytest-cov"
+  "pytest-cov",
+  "deepdiff",
 ]
 
 [tool.pytest.ini_options]

src/cubitpy/cubitpy_to_dat.py → src/cubitpy/cubit_to_fourc_input.py

@@ -30,7 +30,7 @@
 from cubitpy.conf import cupy
 
 
-def add_node_sets(dat_lines, cubit, exo):
+def add_node_sets(cubit, exo, input_file):
     """Add the node sets contained in the cubit session/exo file to the
     dat_lines."""
 
@@ -65,51 +65,39 @@ def add_node_sets(dat_lines, cubit, exo):
         boundary_condition_map[bc_key].append(
             [len(node_sets[geometry_type]), bc_description, names[i_set]]
         )
+        if bc_section not in input_file.inlined.keys():
+            input_file[bc_section] = []
+        bc_description["E"] = len(node_sets[geometry_type])
 
-    # Write the boundary conditions
-    for (bc_section, geo), item in boundary_condition_map.items():
-        dat_lines.append("-" * 40 + bc_section)
-        for set_id, bc_description, name in item:
-            if not name == "":
-                dat_lines.append(f"// {name}")
-            dat_lines.append(f"E {set_id} {bc_description}")
+        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",
-        ],
+        [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:
-            dat_lines.append(section_name)
+            input_file[section_name] = []
             for i_set, node_set in enumerate(node_sets[geo]):
                 node_set.sort()
                 for i_node in node_set:
-                    dat_lines.append(f"NODE {i_node:6d} {set_label} {i_set+1}")
+                    input_file[section_name].append(
+                        {
+                            "type": "NODE",
+                            "node_id": i_node,
+                            "d_type": set_label,
+                            "d_id": i_set + 1,
+                        }
+                    )
 
 
-def get_element_connectivity_string(connectivity):
-    """Return the connectivity string for an element.
+def get_element_connectivity_list(connectivity):
+    """Return the connectivity list for an element.
 
     For hex27 we need a different ordering than the one we get from
-    cubit
+    cubit.
     """
 
     if len(connectivity) == 27:
@@ -143,34 +131,29 @@ def get_element_connectivity_string(connectivity):
             22,
             20,
         ]
-        return " ".join([f"{item:d}" for item in connectivity[ordering]])
+        return [connectivity[i] for i in ordering]
     else:
         # all other elements
-        return " ".join([f"{item:d}" for item in connectivity])
+        return connectivity.tolist()
 
 
-def cubit_to_dat(cubit):
-    """Convert a CubitPy session to a dat file that can be read with 4C."""
+def get_input_file_with_mesh(cubit):
+    """Return a copy of cubit.fourc_input with mesh data (nodes and elements)
+    added."""
 
     # Create exodus file
     os.makedirs(cupy.temp_dir, exist_ok=True)
     exo_path = os.path.join(cupy.temp_dir, "cubitpy.exo")
     cubit.export_exo(exo_path)
     exo = netCDF4.Dataset(exo_path)
 
-    dat_lines = []
-
-    # Add the header
-    for line in cubit.head.split("\n"):
-        dat_lines.append(line.strip())
-
+    # create a deep copy of the input_file
+    input_file = cubit.fourc_input.copy()
     # Add the node sets
-    add_node_sets(dat_lines, cubit, exo)
+    add_node_sets(cubit, exo, input_file)
 
     # Add the nodal data
-    dat_lines.append(
-        "-------------------------------------------------------NODE COORDS"
-    )
+    input_file["NODE COORDS"] = []
     if "coordz" in exo.variables:
         coordinates = np.array(
             [exo.variables["coord" + dim][:] for dim in ["x", "y", "z"]],
@@ -180,28 +163,36 @@ def cubit_to_dat(cubit):
         temp.append([0 for i in range(len(temp[0]))])
         coordinates = np.array(temp).transpose()
     for i, coordinate in enumerate(coordinates):
-        dat_lines.append(
-            f"NODE {i+1:9d} COORD {coordinate[0]: .16e} {coordinate[1]: .16e} {coordinate[2]: .16e}"
+        input_file["NODE COORDS"].append(
+            {
+                "COORD": [coordinate[0], coordinate[1], coordinate[2]],
+                "data": {"type": "NODE"},
+                "id": i + 1,
+            }
         )
 
     # Add the element connectivity
-    current_section = None
     connectivity_keys = [key for key in exo.variables.keys() if "connect" in key]
     connectivity_keys.sort()
     i_element = 0
     for i_block, key in enumerate(connectivity_keys):
-        ele_type, block_string = cubit.blocks[i_block]
-        block_section = ele_type.get_four_c_section()
-        if not block_section == current_section:
-            current_section = block_section
-            dat_lines.append(
-                f"------------------------------------------------{current_section} ELEMENTS"
-            )
+        ele_type, block_dict = cubit.blocks[i_block]
+        block_section = f"{ele_type.get_four_c_section()} ELEMENTS"
+        if block_section not in input_file.sections.keys():
+            input_file[block_section] = []
         for connectivity in exo.variables[key][:]:
-            connectivity_string = get_element_connectivity_string(connectivity)
-            dat_lines.append(
-                f"{i_element+1:9d} {ele_type.get_four_c_name()} {ele_type.get_four_c_type()} {connectivity_string} {block_string}"
+            input_file[block_section].append(
+                {
+                    "id": i_element + 1,
+                    "cell": {
+                        "connectivity": get_element_connectivity_list(connectivity),
+                        "type": ele_type.get_four_c_type(),
+                    },
+                    "data": {
+                        "type": ele_type.get_four_c_name(),
+                        **block_dict,
+                    },
+                }
             )
             i_element += 1
-
-    return dat_lines
+    return input_file

src/cubitpy/cubitpy.py

@@ -26,10 +26,12 @@
 import time
 import warnings
 
+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_wrapper.cubit_wrapper_host import CubitConnect
-from cubitpy.cubitpy_to_dat import cubit_to_dat
 
 
 class CubitPy(object):
@@ -62,13 +64,13 @@ def __init__(self, *, cubit_exe=None, **kwargs):
         # Set lists and counters for blocks and sets
         self._default_cubit_variables()
 
-        # Content of head file
-        self.head = ""
+        self.fourc_input = FourCInput()
 
     def _default_cubit_variables(self):
         """Set the default values for the lists and counters used in cubit."""
         self.blocks = []
         self.node_sets = []
+        self.fourc_input = FourCInput()
 
     def __getattr__(self, key, *args, **kwargs):
         """All calls to methods and attributes that are not in this object get
@@ -121,7 +123,7 @@ def _name_created_set(self, set_type, set_id, name, item):
             self.cubit.cmd('{} {} name "{}"'.format(set_type, set_id, rename_name))
 
     def add_element_type(
-        self, item, el_type, *, name=None, material="MAT 1", bc_description=None
+        self, item, el_type, *, name=None, material=None, bc_description=None
     ):
         """Add a block to cubit that contains the geometry in item. Also set
         the element type of block.
@@ -134,14 +136,20 @@ def add_element_type(
             Cubit element type.
         name: str
             Name of the block.
-        material: str
+        material: dict
             Material string of the block, will be the first part of the BC
             description.
-        bc_description: str
+        bc_description: dict
             Will be written after the material string. If this is not set, the
             default values for the given element type will be used.
         """
 
+        # default values
+        if material is None:
+            material = {"MAT": 1}
+        if bc_description is None:
+            bc_description = {}
+
         # Check that all blocks in cubit are created with this function.
         n_blocks = len(self.blocks)
         if not len(self.cubit.get_block_id_list()) == n_blocks:
@@ -180,11 +188,11 @@ def add_element_type(
         self._name_created_set("block", n_blocks + 1, name, item)
 
         # If the user does not give a bc_description, load the default one.
-        if bc_description is None:
+        if not bc_description:
             bc_description = el_type.get_default_four_c_description()
 
         # Add data that will be written to bc file.
-        self.blocks.append([el_type, " ".join([material, bc_description])])
+        self.blocks.append([el_type, material | bc_description])
 
     def reset_blocks(self):
         """This method deletes all blocks in Cubit and resets the counter in
@@ -203,7 +211,7 @@ def add_node_set(
         *,
         name=None,
         bc_type=None,
-        bc_description="NUMDOF 3 ONOFF 0 0 0 VAL 0 0 0 FUNCT 0 0 0",
+        bc_description=None,
         bc_section=None,
         geometry_type=None,
     ):
@@ -221,7 +229,7 @@ def add_node_set(
         bc_section: str
             Name of the section in the input file. Mutually exclusive with
             bc_type.
-        bc_description: str
+        bc_description: dict
             Definition of the boundary condition.
         geometry_type: cupy.geometry
             Directly set the geometry type, instead of obtaining it from the
@@ -267,6 +275,8 @@ def add_node_set(
             )
         if bc_section is None:
             bc_section = bc_type.get_dat_bc_section_header(geometry_type)
+        if bc_description is None:
+            bc_description = {}
         self.node_sets.append([bc_section, bc_description, geometry_type])
 
     def get_ids(self, geometry_type):
@@ -319,28 +329,22 @@ def export_exo(self, path):
         """Export the mesh."""
         self.cubit.cmd('export mesh "{}" dimension 3 overwrite'.format(path))
 
-    def create_dat(self, dat_path):
-        """Create the dat file an copy it to dat_path.
+    def write_input_file(self, yaml_path):
+        """Create the yaml file an save it in yaml_path.
 
         Args
         ----
-        dat_path: str
+        yaml_path: str
             Path where the input file file will be saved
         """
 
-        # Check if output path exists.
-        if os.path.isabs(dat_path):
-            dat_dir = os.path.dirname(dat_path)
-            if not os.path.exists(dat_dir):
-                raise ValueError("Path {} does not exist!".format(dat_dir))
-
-        with open(dat_path, "w") as the_file:
-            for line in self.get_dat_lines():
-                the_file.write(line + "\n")
+        # Check if output path exists
+        dat_dir = os.path.dirname(os.path.abspath(yaml_path))
+        if not os.path.exists(dat_dir):
+            raise ValueError("Path {} does not exist!".format(dat_dir))
 
-    def get_dat_lines(self):
-        """Return a list with all lines in this input file."""
-        return cubit_to_dat(self)
+        input_file = get_input_file_with_mesh(self)
+        input_file.dump(yaml_path)
 
     def group(self, **kwargs):
         """Reference a group in cubit.

src/cubitpy/cubitpy_types.py

@@ -174,10 +174,9 @@ def get_four_c_name(self):
         """Get the name of this element in 4C."""
 
         # Get the element type parameters.
-        if self == self.hex8sh:
-            return "SOLIDSH8"
-        elif (
-            self == self.hex8
+        if (
+            self == self.hex8sh
+            or self == self.hex8
             or self == self.hex20
             or self == self.hex27
             or self == self.tet10
@@ -264,15 +263,16 @@ def get_default_four_c_description(self):
             or self == self.tet10
             or self == self.wedge6
         ):
-            return "KINEM nonlinear"
+            return {"KINEM": "nonlinear"}
         elif self == self.hex8sh:
-            return "KINEM nonlinear EAS none ANS none THICKDIR auto"
+            return {"KINEM": "nonlinear", "TECH": "shell_eas_ans"}
+
         elif self == self.hex8_fluid or self == self.tet4_fluid:
-            return "NA ALE"
+            return {"NA": "ALE"}
         elif self == self.hex8_thermo or self == self.tet4_thermo:
-            return ""
+            return {}
         elif self == self.hex8_scatra or self == self.tet4_scatra:
-            return ""
+            return {}
         else:
             raise ValueError("Got wrong element type {}!".format(self))