beamme.four_c.input_file

  1# The MIT License (MIT)
  2#
  3# Copyright (c) 2018-2025 BeamMe Authors
  4#
  5# Permission is hereby granted, free of charge, to any person obtaining a copy
  6# of this software and associated documentation files (the "Software"), to deal
  7# in the Software without restriction, including without limitation the rights
  8# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
  9# copies of the Software, and to permit persons to whom the Software is
 10# furnished to do so, subject to the following conditions:
 11#
 12# The above copyright notice and this permission notice shall be included in
 13# all copies or substantial portions of the Software.
 14#
 15# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
 16# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 17# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
 18# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 19# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
 20# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
 21# THE SOFTWARE.
 22"""This module defines the classes that are used to create an input file for
 234C."""
 24
 25import os as _os
 26import sys as _sys
 27from datetime import datetime as _datetime
 28from pathlib import Path as _Path
 29from typing import Dict as _Dict
 30from typing import List as _List
 31
 32from fourcipp.fourc_input import FourCInput as _FourCInput
 33
 34from beamme.core.boundary_condition import BoundaryCondition as _BoundaryCondition
 35from beamme.core.conf import mpy as _mpy
 36from beamme.core.coupling import Coupling as _Coupling
 37from beamme.core.function import Function as _Function
 38from beamme.core.geometry_set import GeometrySet as _GeometrySet
 39from beamme.core.geometry_set import GeometrySetNodes as _GeometrySetNodes
 40from beamme.core.material import Material as _Material
 41from beamme.core.mesh import Mesh as _Mesh
 42from beamme.core.node import Node as _Node
 43from beamme.core.nurbs_patch import NURBSPatch as _NURBSPatch
 44from beamme.four_c.input_file_mappings import (
 45    INPUT_FILE_MAPPINGS as _INPUT_FILE_MAPPINGS,
 46)
 47from beamme.utils.environment import cubitpy_is_available as _cubitpy_is_available
 48from beamme.utils.environment import get_git_data as _get_git_data
 49
 50if _cubitpy_is_available():
 51    import cubitpy as _cubitpy
 52
 53
 54def get_geometry_set_indices_from_section(
 55    section_list: _List, *, append_node_ids: bool = True
 56) -> _Dict:
 57    """Return a dictionary with the geometry set ID as keys and the node IDs as
 58    values.
 59
 60    Args:
 61        section_list: A list with the legacy strings for the geometry pair
 62        append_node_ids: If the node IDs shall be appended, or only the
 63            dict with the keys should be returned.
 64    """
 65
 66    geometry_set_dict: _Dict[int, _List[int]] = {}
 67    for entry in section_list:
 68        id_geometry_set = entry["d_id"]
 69        index_node = entry["node_id"] - 1
 70        if id_geometry_set not in geometry_set_dict:
 71            geometry_set_dict[id_geometry_set] = []
 72        if append_node_ids:
 73            geometry_set_dict[id_geometry_set].append(index_node)
 74
 75    return geometry_set_dict
 76
 77
 78def _dump_coupling(coupling):
 79    """Return the input file representation of the coupling condition."""
 80
 81    # TODO: Move this to a better place / gather all dump functions for general
 82    # BeamMe items in a file or so.
 83
 84    if isinstance(coupling.data, dict):
 85        data = coupling.data
 86    else:
 87        # In this case we have to check which beams are connected to the node.
 88        # TODO: Coupling also makes sense for different beam types, this can
 89        # be implemented at some point.
 90        nodes = coupling.geometry_set.get_points()
 91        connected_elements = [
 92            element for node in nodes for element in node.element_link
 93        ]
 94        element_types = {type(element) for element in connected_elements}
 95        if len(element_types) > 1:
 96            raise TypeError(
 97                f"Expected a single connected type of beam elements, got {element_types}"
 98            )
 99        element_type = element_types.pop()
100        if element_type.beam_type is _mpy.beam.kirchhoff:
101            rotvec = {element.rotvec for element in connected_elements}
102            if len(rotvec) > 1 or not rotvec.pop():
103                raise TypeError(
104                    "Couplings for Kirchhoff beams and rotvec==False not yet implemented."
105                )
106
107        data = element_type.get_coupling_dict(coupling.data)
108
109    return {"E": coupling.geometry_set.i_global, **data}
110
111
112class InputFile(_FourCInput):
113    """An item that represents a complete 4C input file."""
114
115    def __init__(self, sections=None):
116        """Initialize the input file."""
117
118        super().__init__(sections=sections)
119
120        # Contents of NOX xml file.
121        self.nox_xml_contents = ""
122
123        # Register converters to directly convert non-primitive types
124        # to native Python types via the FourCIPP type converter.
125        self.type_converter.register_numpy_types()
126        self.type_converter.register_type(
127            (_Function, _Material, _Node), lambda converter, obj: obj.i_global
128        )
129
130    def add(self, object_to_add, **kwargs):
131        """Add a mesh or a dictionary to the input file.
132
133        Args:
134            object: The object to be added. This can be a mesh or a dictionary.
135            **kwargs: Additional arguments to be passed to the add method.
136        """
137
138        if isinstance(object_to_add, _Mesh):
139            self.add_mesh_to_input_file(mesh=object_to_add, **kwargs)
140
141        else:
142            super().combine_sections(object_to_add)
143
144    def dump(
145        self,
146        input_file_path: str | _Path,
147        *,
148        nox_xml_file: str | None = None,
149        add_header_default: bool = True,
150        add_header_information: bool = True,
151        add_footer_application_script: bool = True,
152        sort_sections=False,
153        validate=True,
154        validate_sections_only: bool = False,
155    ):
156        """Write the input file to disk.
157
158        Args:
159            input_file_path:
160                Path to the input file that should be created.
161            nox_xml_file:
162                If this is a string, the NOX xml file will be created with this
163                name. If this is None, the NOX xml file will be created with the
164                name of the input file with the extension ".nox.xml".
165            add_header_default:
166                Prepend the default header comment to the input file.
167            add_header_information:
168                If the information header should be exported to the input file
169                Contains creation date, git details of BeamMe, CubitPy and
170                original application which created the input file if available.
171            add_footer_application_script:
172                Append the application script which creates the input files as a
173                comment at the end of the input file.
174            sort_sections:
175                Sort sections alphabetically with FourCIPP.
176            validate:
177                Validate if the created input file is compatible with 4C with FourCIPP.
178            validate_sections_only:
179                Validate each section independently. Required sections are no longer
180                required, but the sections must be valid.
181        """
182
183        # Make sure the given input file is a Path instance.
184        input_file_path = _Path(input_file_path)
185
186        if self.nox_xml_contents:
187            if nox_xml_file is None:
188                nox_xml_file = input_file_path.name.split(".")[0] + ".nox.xml"
189
190            self["STRUCT NOX/Status Test"] = {"XML File": nox_xml_file}
191
192            # Write the xml file to the disc.
193            with open(input_file_path.parent / nox_xml_file, "w") as xml_file:
194                xml_file.write(self.nox_xml_contents)
195
196        # Add information header to the input file
197        if add_header_information:
198            self.add({"TITLE": self._get_header()})
199
200        super().dump(
201            input_file_path=input_file_path,
202            sort_sections=sort_sections,
203            validate=validate,
204            validate_sections_only=validate_sections_only,
205            convert_to_native_types=False,  # conversion already happens during add()
206        )
207
208        if add_header_default or add_footer_application_script:
209            with open(input_file_path, "r") as input_file:
210                lines = input_file.readlines()
211
212                if add_header_default:
213                    lines = [
214                        "# " + line + "\n" for line in _mpy.input_file_header
215                    ] + lines
216
217                if add_footer_application_script:
218                    application_path = _Path(_sys.argv[0]).resolve()
219                    lines += self._get_application_script(application_path)
220
221                with open(input_file_path, "w") as input_file:
222                    input_file.writelines(lines)
223
224    def add_mesh_to_input_file(self, mesh: _Mesh):
225        """Add a mesh to the input file.
226
227        Args:
228            mesh: The mesh to be added to the input file.
229        """
230
231        # Perform some checks on the mesh.
232        if _mpy.check_overlapping_elements:
233            mesh.check_overlapping_elements()
234
235        def _get_global_start_geometry_set(dictionary):
236            """Get the indices for the first "real" BeamMe geometry sets."""
237
238            start_indices_geometry_set = {}
239            for geometry_type, section_name in _INPUT_FILE_MAPPINGS[
240                "geometry_sets"
241            ].items():
242                max_geometry_set_id = 0
243                if section_name in dictionary:
244                    section_list = dictionary[section_name]
245                    if len(section_list) > 0:
246                        geometry_set_dict = get_geometry_set_indices_from_section(
247                            section_list, append_node_ids=False
248                        )
249                        max_geometry_set_id = max(geometry_set_dict.keys())
250                start_indices_geometry_set[geometry_type] = max_geometry_set_id
251            return start_indices_geometry_set
252
253        def _get_global_start_node():
254            """Get the index for the first "real" BeamMe node."""
255
256            return len(self.sections.get("NODE COORDS", []))
257
258        def _get_global_start_element():
259            """Get the index for the first "real" BeamMe element."""
260
261            return sum(
262                len(self.sections.get(section, []))
263                for section in ["FLUID ELEMENTS", "STRUCTURE ELEMENTS"]
264            )
265
266        def _get_global_start_material():
267            """Get the index for the first "real" BeamMe material.
268
269            We have to account for materials imported from yaml files
270            that have arbitrary numbering.
271            """
272
273            # Get the maximum material index in materials imported from a yaml file
274            max_material_id = 0
275            section_name = "MATERIALS"
276            if section_name in self.sections:
277                for material in self.sections[section_name]:
278                    max_material_id = max(max_material_id, material["MAT"])
279            return max_material_id
280
281        def _get_global_start_function():
282            """Get the index for the first "real" BeamMe function."""
283
284            max_function_id = 0
285            for section_name in self.sections.keys():
286                if section_name.startswith("FUNCT"):
287                    max_function_id = max(
288                        max_function_id, int(section_name.split("FUNCT")[-1])
289                    )
290            return max_function_id
291
292        def _set_i_global(data_list, *, start_index=0):
293            """Set i_global in every item of data_list."""
294
295            # A check is performed that every entry in data_list is unique.
296            if len(data_list) != len(set(data_list)):
297                raise ValueError("Elements in data_list are not unique!")
298
299            # Set the values for i_global.
300            for i, item in enumerate(data_list):
301                # TODO make i_global index-0 based
302                item.i_global = i + 1 + start_index
303
304        def _set_i_global_elements(element_list, *, start_index=0):
305            """Set i_global in every item of element_list."""
306
307            # A check is performed that every entry in element_list is unique.
308            if len(element_list) != len(set(element_list)):
309                raise ValueError("Elements in element_list are not unique!")
310
311            # Set the values for i_global.
312            i = start_index
313            i_nurbs_patch = 0
314            for item in element_list:
315                # As a NURBS patch can be defined with more elements, an offset is applied to the
316                # rest of the items
317                # TODO make i_global index-0 based
318                item.i_global = i + 1
319                if isinstance(item, _NURBSPatch):
320                    item.n_nurbs_patch = i_nurbs_patch + 1
321                    offset = item.get_number_elements()
322                    i += offset
323                    i_nurbs_patch += 1
324                else:
325                    i += 1
326
327        def _dump_mesh_items(section_name, data_list):
328            """Output a section name and apply either the default dump or the
329            specialized the dump_to_list for each list item."""
330
331            # Do not write section if no content is available
332            if len(data_list) == 0:
333                return
334
335            list = []
336
337            for item in data_list:
338                if (
339                    isinstance(item, _GeometrySet)
340                    or isinstance(item, _GeometrySetNodes)
341                    or isinstance(item, _NURBSPatch)
342                ):
343                    list.extend(item.dump_to_list())
344                elif hasattr(item, "dump_to_list"):
345                    list.append(item.dump_to_list())
346                elif isinstance(item, _BoundaryCondition):
347                    list.append(
348                        {
349                            "E": item.geometry_set.i_global,
350                            **item.data,
351                        }
352                    )
353
354                elif isinstance(item, _Coupling):
355                    list.append(_dump_coupling(item))
356                else:
357                    raise TypeError(f"Could not dump {item}")
358
359            # If section already exists, retrieve from input file and
360            # add newly. We always need to go through fourcipp to convert
361            # the data types correctly.
362            if section_name in self.sections:
363                existing_entries = self.pop(section_name)
364                existing_entries.extend(list)
365                list = existing_entries
366
367            self.add({section_name: list})
368
369        # Add sets from couplings and boundary conditions to a temp container.
370        mesh.unlink_nodes()
371        start_indices_geometry_set = _get_global_start_geometry_set(self.sections)
372        mesh_sets = mesh.get_unique_geometry_sets(
373            geometry_set_start_indices=start_indices_geometry_set
374        )
375
376        # Assign global indices to all entries.
377        start_index_nodes = _get_global_start_node()
378        _set_i_global(mesh.nodes, start_index=start_index_nodes)
379
380        start_index_elements = _get_global_start_element()
381        _set_i_global_elements(mesh.elements, start_index=start_index_elements)
382
383        start_index_materials = _get_global_start_material()
384        _set_i_global(mesh.materials, start_index=start_index_materials)
385
386        start_index_functions = _get_global_start_function()
387        _set_i_global(mesh.functions, start_index=start_index_functions)
388
389        # Add material data to the input file.
390        _dump_mesh_items("MATERIALS", mesh.materials)
391
392        # Add the functions.
393        for function in mesh.functions:
394            self.add({f"FUNCT{function.i_global}": function.data})
395
396        # If there are couplings in the mesh, set the link between the nodes
397        # and elements, so the couplings can decide which DOFs they couple,
398        # depending on the type of the connected beam element.
399        def get_number_of_coupling_conditions(key):
400            """Return the number of coupling conditions in the mesh."""
401            if (key, _mpy.geo.point) in mesh.boundary_conditions.keys():
402                return len(mesh.boundary_conditions[key, _mpy.geo.point])
403            else:
404                return 0
405
406        if (
407            get_number_of_coupling_conditions(_mpy.bc.point_coupling)
408            + get_number_of_coupling_conditions(_mpy.bc.point_coupling_penalty)
409            > 0
410        ):
411            mesh.set_node_links()
412
413        # Add the boundary conditions.
414        for (bc_key, geom_key), bc_list in mesh.boundary_conditions.items():
415            if len(bc_list) > 0:
416                section_name = (
417                    bc_key
418                    if isinstance(bc_key, str)
419                    else _INPUT_FILE_MAPPINGS["boundary_conditions"][(bc_key, geom_key)]
420                )
421                _dump_mesh_items(section_name, bc_list)
422
423        # Add additional element sections, e.g., for NURBS knot vectors.
424        for element in mesh.elements:
425            element.dump_element_specific_section(self)
426
427        # Add the geometry sets.
428        for geom_key, item in mesh_sets.items():
429            if len(item) > 0:
430                _dump_mesh_items(_INPUT_FILE_MAPPINGS["geometry_sets"][geom_key], item)
431
432        # Add the nodes and elements.
433        _dump_mesh_items("NODE COORDS", mesh.nodes)
434        _dump_mesh_items("STRUCTURE ELEMENTS", mesh.elements)
435        # TODO: reset all links and counters set in this method.
436
437    def _get_header(self) -> dict:
438        """Return the information header for the current BeamMe run.
439
440        Returns:
441            A dictionary with the header information.
442        """
443
444        header: dict = {"BeamMe": {}}
445
446        header["BeamMe"]["creation_date"] = _datetime.now().isoformat(
447            sep=" ", timespec="seconds"
448        )
449
450        # application which created the input file
451        application_path = _Path(_sys.argv[0]).resolve()
452        header["BeamMe"]["Application"] = {"path": str(application_path)}
453
454        application_git_sha, application_git_date = _get_git_data(
455            application_path.parent
456        )
457        if application_git_sha is not None and application_git_date is not None:
458            header["BeamMe"]["Application"].update(
459                {
460                    "git_sha": application_git_sha,
461                    "git_date": application_git_date,
462                }
463            )
464
465        # BeamMe information
466        beamme_git_sha, beamme_git_date = _get_git_data(
467            _Path(__file__).resolve().parent
468        )
469        if beamme_git_sha is not None and beamme_git_date is not None:
470            header["BeamMe"]["BeamMe"] = {
471                "git_SHA": beamme_git_sha,
472                "git_date": beamme_git_date,
473            }
474
475        # CubitPy information
476        if _cubitpy_is_available():
477            cubitpy_git_sha, cubitpy_git_date = _get_git_data(
478                _os.path.dirname(_cubitpy.__file__)
479            )
480
481            if cubitpy_git_sha is not None and cubitpy_git_date is not None:
482                header["BeamMe"]["CubitPy"] = {
483                    "git_SHA": cubitpy_git_sha,
484                    "git_date": cubitpy_git_date,
485                }
486
487        return header
488
489    def _get_application_script(self, application_path: _Path) -> list[str]:
490        """Get the script that created this input file.
491
492        Args:
493            application_path: Path to the script that created this input file.
494        Returns:
495            A list of strings with the script that created this input file.
496        """
497
498        application_script_lines = [
499            "# Application script which created this input file:\n"
500        ]
501
502        with open(application_path) as script_file:
503            application_script_lines.extend("# " + line for line in script_file)
504
505        return application_script_lines