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