Coverage for src/beamme/four_c/input_file_dump_functions.py: 95%
160 statements
« prev ^ index » next coverage.py v7.14.1, created at 2026-06-08 11:03 +0000
1# The MIT License (MIT)
2#
3# Copyright (c) 2018-2026 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 file defines functions to dump mesh items for 4C."""
24from collections import defaultdict as _defaultdict
25from typing import Any as _Any
26from typing import KeysView as _KeysView
27from typing import List as _List
29import numpy as _np
30from fourcipp.fourc_input import FourCInput as _FourCInput
32from beamme.core.boundary_condition import BoundaryCondition as _BoundaryCondition
33from beamme.core.conf import Geometry as _Geometry
34from beamme.core.conf import bme as _bme
35from beamme.core.coupling import Coupling as _Coupling
36from beamme.core.function import Function as _Function
37from beamme.core.geometry_set import GeometrySetBase as _GeometrySetBase
38from beamme.core.material import Material as _Material
39from beamme.core.mesh import Mesh as _Mesh
40from beamme.core.mesh_representation import MeshRepresentation as _MeshRepresentation
41from beamme.core.mesh_representation import (
42 merge_mesh_representations as _merge_mesh_representations,
43)
44from beamme.core.mesh_representation import (
45 string_to_geometry_set_info as _string_to_geometry_set_info,
46)
47from beamme.core.nurbs_patch import NURBSPatch as _NURBSPatch
48from beamme.core.rotation import Rotation as _Rotation
49from beamme.four_c.element_data import FourCElementData as _FourCElementData
50from beamme.four_c.four_c_types import (
51 BeamKirchhoffParametrizationType as _BeamKirchhoffParametrizationType,
52)
53from beamme.four_c.four_c_types import BeamType as _BeamType
54from beamme.four_c.input_file_mappings import (
55 INPUT_FILE_MAPPINGS as _INPUT_FILE_MAPPINGS,
56)
57from beamme.four_c.material import (
58 get_material_to_i_global_mapping as _get_material_to_i_global_mapping,
59)
62def dump_function(function: _Function, i_global: int) -> dict[str, _Any]:
63 """Return the representation of a function in the 4C input file."""
64 return {f"FUNCT{i_global + 1}": function.data}
67def dump_coupling(coupling):
68 """Return the input file representation of the coupling condition."""
70 if isinstance(coupling.data, dict):
71 data = coupling.data
72 else:
73 # In this case we have to check which beams are connected to the node.
74 # TODO: Coupling also makes sense for different beam types, this can
75 # be implemented at some point.
76 nodes = coupling.geometry_set.get_points()
77 connected_elements = [
78 element for node in nodes for element in node.element_link
79 ]
80 element_types = {type(element) for element in connected_elements}
81 if len(element_types) > 1:
82 raise TypeError(
83 f"Expected a single connected type of beam elements, got {element_types}"
84 )
85 element_type = element_types.pop()
86 if element_type.beam_type is _BeamType.kirchhoff:
87 unique_parametrization_flags = {
88 type(element).kirchhoff_parametrization
89 for element in connected_elements
90 }
91 if (
92 len(unique_parametrization_flags) > 1
93 or not unique_parametrization_flags.pop()
94 == _BeamKirchhoffParametrizationType.rot
95 ):
96 raise TypeError(
97 "Couplings for Kirchhoff beams and tangent "
98 "based parametrization not yet implemented."
99 )
101 data = element_type.get_coupling_dict(coupling.data)
103 return {"E": coupling.geometry_set, **data}
106def dump_nurbs_patch_knotvectors(input_file, nurbs_patch: _NURBSPatch) -> None:
107 """Set the knot vectors of the NURBS patch in the input file."""
109 patch_data: dict[str, _Any] = {
110 "KNOT_VECTORS": [],
111 }
113 for dir_manifold in range(nurbs_patch.get_nurbs_dimension()):
114 knotvector = nurbs_patch.knot_vectors[dir_manifold]
115 num_knots = len(knotvector)
117 # Check the type of knot vector, in case that the multiplicity of the first and last
118 # knot vectors is not p + 1, then it is a closed (periodic) knot vector, otherwise it
119 # is an open (interpolated) knot vector.
120 knotvector_type = "Interpolated"
122 for i in range(nurbs_patch.polynomial_orders[dir_manifold] - 1):
123 if (abs(knotvector[i] - knotvector[i + 1]) > _bme.eps_knot_vector) or (
124 abs(knotvector[num_knots - 2 - i] - knotvector[num_knots - 1 - i])
125 > _bme.eps_knot_vector
126 ):
127 knotvector_type = "Periodic"
128 break
130 patch_data["KNOT_VECTORS"].append(
131 {
132 "DEGREE": nurbs_patch.polynomial_orders[dir_manifold],
133 "TYPE": knotvector_type,
134 "KNOTS": [
135 knot_vector_val
136 for knot_vector_val in nurbs_patch.knot_vectors[dir_manifold]
137 ],
138 }
139 )
141 if "STRUCTURE KNOTVECTORS" in input_file:
142 # Get all existing patches in the input file - they will be added to the
143 # input file again at the end of this function. By doing it this way, the
144 # FourCIPP type converter will be applied to the current patch.
145 # This also means that we apply the type converter again already existing
146 # patches. But, with the usual number of patches and data size, this
147 # should not lead to a measurable performance impact.
148 patches = input_file.pop("STRUCTURE KNOTVECTORS")["PATCHES"]
149 else:
150 patches = []
152 patch_data["ID"] = nurbs_patch
153 patches.append(patch_data)
154 input_file.add({"STRUCTURE KNOTVECTORS": {"PATCHES": patches}})
157def dump_mesh_to_input_file(input_file, mesh: _Mesh) -> None:
158 """Add a mesh to the input file.
160 Internally, we store the geometry information from the mesh in the mesh
161 representation of the input file. All other information, e.g., element
162 types, materials, boundary conditions and functions will be dumped to the
163 4C input file via FourCIPP.
165 Args:
166 input_file: The input file where we want to add the mesh information.
167 mesh: The mesh to be added to the input file.
168 """
170 # Compute starting index for element types
171 start_index_element_types = len(input_file.element_type_id_to_data)
173 # Compute starting index for NURBS patches
174 nurbs_patches = input_file.sections.get("STRUCTURE KNOTVECTORS", {}).get(
175 "PATCHES", []
176 )
177 start_index_nurbs_patches = len(nurbs_patches)
179 # Compute starting index for geometry sets
180 start_index_geometry_set = max(
181 (
182 entry["d_id"] # We don't need a + 1 here, as this index is index-1 based.
183 for section_name in _INPUT_FILE_MAPPINGS[
184 "geometry_sets_geometry_to_condition_name"
185 ].values()
186 for entry in input_file.sections.get(section_name, [])
187 ),
188 default=0,
189 )
190 for name in input_file.mesh_representation.point_data.keys():
191 info = _string_to_geometry_set_info(name)
192 if info is not None:
193 start_index_geometry_set = max(start_index_geometry_set, info.i_global + 1)
195 # Compute starting index for functions
196 start_index_functions = max(
197 (
198 int(section.split("FUNCT")[-1])
199 for section in input_file.sections
200 if section.startswith("FUNCT")
201 ),
202 default=0,
203 )
205 # Compute starting index for materials
206 start_index_materials = max(
207 (material["MAT"] for material in input_file.sections.get("MATERIALS", [])),
208 default=0,
209 )
211 # Get material to global index mapping for the materials in the mesh.
212 material_to_i_global = _get_material_to_i_global_mapping(mesh.materials)
214 # Get the mesh representation for the mesh.
215 (
216 mesh_representation,
217 mesh_element_type_id_to_data,
218 geometry_sets_to_i_global,
219 nurbs_patch_to_i_global,
220 ) = mesh.get_mesh_representation(material_to_i_global)
221 mesh_representation.offset_indices(
222 element_type_id_offset=start_index_element_types,
223 material_offset=start_index_materials,
224 geometry_set_offset=start_index_geometry_set,
225 )
227 # Add the new element types to the mapping in the input file.
228 for element_type_id, element_type_data in mesh_element_type_id_to_data.items():
229 input_file.element_type_id_to_data[
230 element_type_id + start_index_element_types
231 ] = element_type_data
232 input_file.mesh_representation = _merge_mesh_representations(
233 input_file.mesh_representation, mesh_representation
234 )
236 # Dump functions to the input file.
237 function_to_i_global: dict[_Function, int] = {}
238 for i_global, function in enumerate(mesh.functions, start=start_index_functions):
239 function_to_i_global[function] = i_global
240 input_file.add(dump_function(function, i_global))
241 if len(mesh.functions) != len(function_to_i_global):
242 raise ValueError("Functions are not unique!")
244 # Adapt the FourCIPP type converters such that the correct indices will be dumped.
245 input_file.fourc_input.type_converter.register_type(
246 _GeometrySetBase,
247 lambda _, obj: geometry_sets_to_i_global[obj] + 1 + start_index_geometry_set,
248 )
249 input_file.fourc_input.type_converter.register_type(
250 _Function, lambda _, obj: function_to_i_global[obj] + 1
251 )
252 input_file.fourc_input.type_converter.register_type(
253 _Material,
254 lambda _, obj: material_to_i_global[obj] + 1 + start_index_materials,
255 )
256 input_file.fourc_input.type_converter.register_type(
257 _NURBSPatch,
258 lambda _, obj: nurbs_patch_to_i_global[obj] + 1 + start_index_nurbs_patches,
259 )
261 def _dump(section_name: str, items: _List | _KeysView) -> None:
262 """Dump list of items to a section in the input file.
264 This function ensures that the dumped items will be appended to the
265 existing section in the input file, and that the full data is parsed
266 through the FourCIPP type converter again.
268 Args:
269 section_name: Name of the section
270 items: List of items to be dumped
271 """
272 if not items:
273 return
274 dumped: list[_Any] = []
275 for item in items:
276 dump_item = None
277 if hasattr(item, "dump_to_list"):
278 dump_item = item.dump_to_list()
279 elif isinstance(item, _BoundaryCondition):
280 dump_item = {"E": item.geometry_set, **item.data}
281 elif isinstance(item, _Coupling):
282 dump_item = dump_coupling(item)
283 else:
284 raise TypeError(f"Could not dump {item}")
285 dumped.append(dump_item)
287 # Go through FourCIPP to convert to native types
288 full_item_list = input_file.pop(section_name, [])
289 full_item_list.extend(dumped)
290 input_file.add({section_name: full_item_list})
292 # Dump materials
293 _dump("MATERIALS", material_to_i_global.keys())
295 # Dump couplings
296 # If there are couplings in the mesh, set the link between the nodes
297 # and elements, so the couplings can decide which DOFs they couple,
298 # depending on the type of the connected beam element.
299 if any(
300 mesh.boundary_conditions.get((key, _bme.geo.point), [])
301 for key in (_bme.bc.point_coupling, _bme.bc.point_coupling_penalty)
302 ):
303 is_linked_nodes = True
304 mesh.set_node_links()
305 else:
306 is_linked_nodes = False
308 # Dump boundary conditions
309 for (bc_key, geometry_key), bc_list in mesh.boundary_conditions.items():
310 if bc_list:
311 section = (
312 bc_key
313 if isinstance(bc_key, str)
314 else _INPUT_FILE_MAPPINGS["boundary_conditions"][(bc_key, geometry_key)]
315 )
316 _dump(section, bc_list)
318 # If we have previously set the node links, we unlink them here.
319 if is_linked_nodes:
320 mesh.unlink_nodes()
322 # Dump NURBS patch information.
323 for element in mesh.elements:
324 if isinstance(element, _NURBSPatch):
325 dump_nurbs_patch_knotvectors(input_file, element)
328def dump_mesh_representation_to_input_file_legacy(
329 fourc_input: _FourCInput,
330 mesh_representation: _MeshRepresentation,
331 element_type_id_to_data: dict[int, _FourCElementData],
332) -> None:
333 """Dump the information contained in the mesh representation to the 4C
334 input file via FourCIPP, in legacy format.
336 Args:
337 fourc_input: 4C input file via FourCIPP where the mesh information data will be dumped to.
338 mesh_representation: The mesh representation that is added to the input file.
339 element_type_id_to_data: The mapping between element type ID and the element type data.
340 """
342 # Compute the starting indices for the nodes and elements entities.
343 start_index_nodes = len(fourc_input.sections.get("NODE COORDS", []))
344 start_index_elements = sum(
345 len(fourc_input.sections.get(section, []))
346 for section in ("FLUID ELEMENTS", "STRUCTURE ELEMENTS")
347 )
349 def _dump(section_name: str, dictionary_list: _List):
350 """Append the given list of dictionaries to the section in the input
351 file."""
352 full_item_list = fourc_input.pop(section_name, [])
353 full_item_list.extend(dictionary_list)
354 fourc_input.combine_sections({section_name: full_item_list})
356 # Dump node information to the input file.
357 node_data_list = []
358 for i_node, (point, point_type, cp_weight) in enumerate(
359 zip(
360 mesh_representation.points,
361 mesh_representation.data_iterator("point_data", "point_type"),
362 mesh_representation.data_iterator("point_data", "control_point_weight"),
363 ),
364 start=start_index_nodes,
365 ):
366 node_type = _bme.node_type(point_type)
367 node_id = i_node + 1
368 if node_type == _bme.node_type.node or node_type == _bme.node_type.cosserat:
369 node_data_list.append(
370 {
371 "id": node_id,
372 "COORD": point,
373 "data": {"type": "NODE"},
374 }
375 )
376 elif node_type == _bme.node_type.control_point:
377 node_data_list.append(
378 {
379 "id": node_id,
380 "COORD": point,
381 "data": {
382 "type": "CP",
383 "weight": cp_weight,
384 },
385 }
386 )
387 else:
388 raise ValueError(f"Unknown node type {node_type} for node {i_node + 1}")
389 _dump("NODE COORDS", node_data_list)
391 # Dump element information to the input file.
392 element_list = []
393 for i_element, (
394 connectivity,
395 element_type_id,
396 element_material_id,
397 ) in enumerate(
398 zip(
399 mesh_representation.connectivity_iterator(),
400 mesh_representation.data_iterator("cell_data", "element_type_id"),
401 mesh_representation.data_iterator("cell_data", "material_id"),
402 )
403 ):
404 element_type_data = element_type_id_to_data[element_type_id]
406 node_ordering = _INPUT_FILE_MAPPINGS[
407 "four_c_cell_to_connectivity_mapping_from_vtk"
408 ].get(element_type_data.four_c_cell, None)
409 if node_ordering is not None:
410 connectivity = connectivity[node_ordering]
412 additional_element_data = None
413 if _INPUT_FILE_MAPPINGS["four_c_type_to_requires_triads"].get(
414 element_type_data.four_c_type, False
415 ):
416 # The numpy quaternion package can return rotation vectors outside
417 # the range of -pi to pi, which can cause issues in testing comparisons.
418 # To avoid this, we convert the rotations to BeamMe internal rotations
419 # and extract the rotation vectors again, ensuring they are in the
420 # correct range.
421 # This is super slow, but for now keep it, as a mesh based output is to
422 # be preferred when performance is of importance.
423 if "rotation_vector" not in mesh_representation.point_data:
424 raise KeyError(
425 f"Rotation vectors are required for type {element_type_data.four_c_type}, but "
426 "no rotation vectors found in the mesh representation!"
427 )
428 rotations = [
429 _Rotation.from_rotation_vector(rotation_vector)
430 for rotation_vector in mesh_representation.point_data[
431 "rotation_vector"
432 ][connectivity]
433 ]
434 additional_element_data = {
435 "TRIADS": _np.array(
436 [rotation.get_rotation_vector() for rotation in rotations]
437 ).ravel()
438 }
440 element_list.append(
441 element_type_data.get_legacy_dict(
442 element_id=start_index_elements + i_element,
443 connectivity=start_index_nodes + connectivity,
444 element_material_id=element_material_id,
445 additional_element_data=additional_element_data,
446 )
447 )
448 _dump("STRUCTURE ELEMENTS", element_list)
450 # Dump geometry sets to the input file.
451 # We first create a mapping from the geometry type to a dictionary which maps
452 # the global geometry set ID to the name of the corresponding data array in the
453 # mesh representation. This is required for the sorting later on.
454 geometry_type_to_geometry_sets: dict[_Geometry, dict[int, str]] = _defaultdict(dict)
455 for name in mesh_representation.point_data.keys():
456 info = _string_to_geometry_set_info(name)
457 if info is not None:
458 geometry_type_to_geometry_sets[info.geometry_type][info.i_global] = name
460 for geometry_type, id_to_name_map in geometry_type_to_geometry_sets.items():
461 geometry_set_list = []
462 # We sort the keys here to ensure that the geometry sets are dumped in the
463 # correct order.
464 sorted_ids = sorted(id_to_name_map.keys())
465 for i_global in sorted_ids:
466 name = id_to_name_map[i_global]
467 node_indices = _np.nonzero(mesh_representation.point_data[name])[0]
468 geometry_set_list.extend(
469 [
470 {
471 "type": "NODE",
472 "node_id": start_index_nodes + node_index + 1,
473 "d_type": _INPUT_FILE_MAPPINGS[
474 "geometry_sets_geometry_to_entry_name"
475 ][geometry_type],
476 "d_id": i_global + 1,
477 }
478 for node_index in node_indices
479 ]
480 )
481 _dump(
482 _INPUT_FILE_MAPPINGS["geometry_sets_geometry_to_condition_name"][
483 geometry_type
484 ],
485 geometry_set_list,
486 )