Coverage for src/beamme/four_c/input_file_dump_functions.py: 93%
246 statements
« prev ^ index » next coverage.py v7.17.2, created at 2026-11-05 08:51 +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 collections.abc import KeysView as _KeysView
26from typing import Any as _Any
28import numpy as _np
29import pyvista as _pv
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 GeometrySetInfo as _GeometrySetInfo
41from beamme.core.mesh_representation import MeshRepresentation as _MeshRepresentation
42from beamme.core.mesh_representation import (
43 merge_mesh_representations as _merge_mesh_representations,
44)
45from beamme.core.nurbs_patch import NURBSPatch as _NURBSPatch
46from beamme.core.rotation import Rotation as _Rotation
47from beamme.four_c.boundary_condition_data import (
48 FourCBoundaryConditionData as _FourCBoundaryConditionData,
49)
50from beamme.four_c.element_data import FourCElementData as _FourCElementData
51from beamme.four_c.four_c_types import (
52 BeamKirchhoffParametrizationType as _BeamKirchhoffParametrizationType,
53)
54from beamme.four_c.four_c_types import BeamType as _BeamType
55from beamme.four_c.input_file_mappings import (
56 INPUT_FILE_MAPPINGS as _INPUT_FILE_MAPPINGS,
57)
58from beamme.four_c.material import (
59 get_material_to_i_global_mapping as _get_material_to_i_global_mapping,
60)
61from beamme.utils.data_structures import (
62 create_inverse_mapping as _create_inverse_mapping,
63)
66def get_geometry_set_infos_with_nodal_flag_vector(
67 mesh_representation: _MeshRepresentation,
68) -> list[_GeometrySetInfo]:
69 """Get list of geometry set infos with nodal flag vectors.
71 Args:
72 mesh_representation: The mesh representation from which we want to get the
73 geometry set infos.
75 Returns:
76 A list of geometry set infos. If the geometry set is a node-set we return it as
77 is, if it is a cell-set we compute the nodes connected to the cells and add
78 also return this vector.
79 """
80 geometry_set_infos = mesh_representation.get_geometry_set_infos()
81 for info in geometry_set_infos:
82 if info.point_flag_vector is None:
83 unique_values = _np.unique(
84 _np.concatenate(
85 list(
86 mesh_representation.connectivity_iterator(
87 element_indices=_np.flatnonzero(info.cell_flag_vector)
88 )
89 )
90 )
91 )
92 point_flag_vector = _np.zeros(mesh_representation.n_points, dtype=int)
93 point_flag_vector[unique_values] = 1
94 info.point_flag_vector = point_flag_vector
96 return geometry_set_infos
99def dump_function(function: _Function, i_global: int) -> dict[str, _Any]:
100 """Return the representation of a function in the 4C input file."""
101 return {f"FUNCT{i_global + 1}": function.data}
104def dump_coupling(coupling):
105 """Return the input file representation of the coupling condition."""
106 if isinstance(coupling.data, dict):
107 data = coupling.data
108 else:
109 # In this case we have to check which beams are connected to the node.
110 # TODO: Coupling also makes sense for different beam types, this can
111 # be implemented at some point.
112 nodes = coupling.geometry_set.get_points()
113 connected_elements = [
114 element for node in nodes for element in node.element_link
115 ]
116 element_types = {type(element) for element in connected_elements}
117 if len(element_types) > 1:
118 raise TypeError(
119 f"Expected a single connected type of beam elements, got {element_types}"
120 )
121 element_type = element_types.pop()
122 if element_type.beam_type is _BeamType.kirchhoff:
123 unique_parametrization_flags = {
124 type(element).kirchhoff_parametrization
125 for element in connected_elements
126 }
127 if (
128 len(unique_parametrization_flags) > 1
129 or not unique_parametrization_flags.pop()
130 == _BeamKirchhoffParametrizationType.rot
131 ):
132 raise TypeError(
133 "Couplings for Kirchhoff beams and tangent "
134 "based parametrization not yet implemented."
135 )
137 data = element_type.get_coupling_dict(coupling.data)
139 return data
142def dump_nurbs_patch_knotvectors(input_file, nurbs_patch: _NURBSPatch) -> None:
143 """Set the knot vectors of the NURBS patch in the input file."""
144 patch_data: dict[str, _Any] = {
145 "KNOT_VECTORS": [],
146 }
148 for dir_manifold in range(nurbs_patch.get_nurbs_dimension()):
149 knotvector = nurbs_patch.knot_vectors[dir_manifold]
150 num_knots = len(knotvector)
152 # Check the type of knot vector, in case that the multiplicity of the first and last
153 # knot vectors is not p + 1, then it is a closed (periodic) knot vector, otherwise it
154 # is an open (interpolated) knot vector.
155 knotvector_type = "Interpolated"
157 for i in range(nurbs_patch.polynomial_orders[dir_manifold] - 1):
158 if (abs(knotvector[i] - knotvector[i + 1]) > _bme.eps_knot_vector) or (
159 abs(knotvector[num_knots - 2 - i] - knotvector[num_knots - 1 - i])
160 > _bme.eps_knot_vector
161 ):
162 knotvector_type = "Periodic"
163 break
165 patch_data["KNOT_VECTORS"].append(
166 {
167 "DEGREE": nurbs_patch.polynomial_orders[dir_manifold],
168 "TYPE": knotvector_type,
169 "KNOTS": [
170 knot_vector_val
171 for knot_vector_val in nurbs_patch.knot_vectors[dir_manifold]
172 ],
173 }
174 )
176 if "STRUCTURE KNOTVECTORS" in input_file:
177 # Get all existing patches in the input file - they will be added to the
178 # input file again at the end of this function. By doing it this way, the
179 # FourCIPP type converter will be applied to the current patch.
180 # This also means that we apply the type converter again already existing
181 # patches. But, with the usual number of patches and data size, this
182 # should not lead to a measurable performance impact.
183 patches = input_file.pop("STRUCTURE KNOTVECTORS")["PATCHES"]
184 else:
185 patches = []
187 patch_data["ID"] = nurbs_patch
188 patches.append(patch_data)
189 input_file.add({"STRUCTURE KNOTVECTORS": {"PATCHES": patches}})
192def dump_mesh_to_input_file(input_file, mesh: _Mesh) -> None:
193 """Add a mesh to the input file.
195 Internally, we store the geometry information from the mesh in the mesh
196 representation of the input file. All other information, e.g., element
197 types, materials, boundary conditions and functions will be dumped to the
198 4C input file via FourCIPP.
200 Args:
201 input_file: The input file where we want to add the mesh information.
202 mesh: The mesh to be added to the input file.
203 """
204 # Compute starting index for element types
205 start_index_element_types = len(input_file.element_type_id_to_data)
207 # Compute starting index for NURBS patches
208 nurbs_patches = input_file.sections.get("STRUCTURE KNOTVECTORS", {}).get(
209 "PATCHES", []
210 )
211 start_index_nurbs_patches = len(nurbs_patches)
213 # Compute starting index for geometry sets
214 start_index_geometry_set = max(
215 (
216 entry["d_id"] # We don't need a + 1 here, as this index is index-1 based.
217 for section_name in _INPUT_FILE_MAPPINGS[
218 "geometry_sets_geometry_to_condition_name"
219 ].values()
220 for entry in input_file.sections.get(section_name, [])
221 ),
222 default=0,
223 )
224 for info in input_file.mesh_representation.get_geometry_set_infos():
225 start_index_geometry_set = max(start_index_geometry_set, info.i_global + 1)
227 # Compute starting index for functions
228 start_index_functions = max(
229 (
230 int(section.split("FUNCT")[-1])
231 for section in input_file.sections
232 if section.startswith("FUNCT")
233 ),
234 default=0,
235 )
237 # Compute starting index for materials
238 start_index_materials = max(
239 (material["MAT"] for material in input_file.sections.get("MATERIALS", [])),
240 default=0,
241 )
243 # Get material to global index mapping for the materials in the mesh.
244 material_to_i_global = _get_material_to_i_global_mapping(mesh.materials)
246 # Get the mesh representation for the mesh.
247 (
248 mesh_representation,
249 mesh_element_type_id_to_data,
250 geometry_sets_to_i_global,
251 nurbs_patch_to_i_global,
252 ) = mesh.get_mesh_representation(material_to_i_global)
253 mesh_representation.offset_indices(
254 element_type_id_offset=start_index_element_types,
255 material_offset=start_index_materials,
256 geometry_set_offset=start_index_geometry_set,
257 )
259 # Add the new element types to the mapping in the input file.
260 for element_type_id, element_type_data in mesh_element_type_id_to_data.items():
261 input_file.element_type_id_to_data[
262 element_type_id + start_index_element_types
263 ] = element_type_data
264 input_file.mesh_representation = _merge_mesh_representations(
265 input_file.mesh_representation, mesh_representation
266 )
268 # Dump functions to the input file.
269 function_to_i_global: dict[_Function, int] = {}
270 for i_global, function in enumerate(mesh.functions, start=start_index_functions):
271 function_to_i_global[function] = i_global
272 input_file.add(dump_function(function, i_global))
273 if len(mesh.functions) != len(function_to_i_global):
274 raise ValueError("Functions are not unique!")
276 # Adapt the FourCIPP type converters such that the correct indices will be dumped.
277 input_file.fourc_input.type_converter.register_type(
278 _GeometrySetBase,
279 lambda _, obj: geometry_sets_to_i_global[obj] + 1 + start_index_geometry_set,
280 )
281 input_file.fourc_input.type_converter.register_type(
282 _Function, lambda _, obj: function_to_i_global[obj] + 1
283 )
284 input_file.fourc_input.type_converter.register_type(
285 _Material,
286 lambda _, obj: material_to_i_global[obj] + 1 + start_index_materials,
287 )
288 input_file.fourc_input.type_converter.register_type(
289 _NURBSPatch,
290 lambda _, obj: nurbs_patch_to_i_global[obj] + 1 + start_index_nurbs_patches,
291 )
293 def _dump(section_name: str, items: list | _KeysView) -> None:
294 """Dump list of items to a section in the input file.
296 This function ensures that the dumped items will be appended to the
297 existing section in the input file, and that the full data is parsed
298 through the FourCIPP type converter again.
300 Args:
301 section_name: Name of the section
302 items: List of items to be dumped
303 """
304 if not items:
305 return
306 dumped: list[_Any] = []
307 for item in items:
308 dump_item = None
309 if hasattr(item, "dump_to_list"):
310 dump_item = item.dump_to_list()
311 elif isinstance(item, _BoundaryCondition):
312 dump_item = {"E": item.geometry_set, **item.data}
313 elif isinstance(item, _Coupling):
314 dump_item = dump_coupling(item)
315 else:
316 raise TypeError(f"Could not dump {item}")
317 dumped.append(dump_item)
319 # Go through FourCIPP to convert to native types
320 full_item_list = input_file.pop(section_name, [])
321 full_item_list.extend(dumped)
322 input_file.add({section_name: full_item_list})
324 # Dump materials
325 _dump("MATERIALS", material_to_i_global.keys())
327 # Dump couplings
328 # If there are couplings in the mesh, set the link between the nodes
329 # and elements, so the couplings can decide which DOFs they couple,
330 # depending on the type of the connected beam element.
331 if any(
332 mesh.boundary_conditions.get((key, _bme.geo.point), [])
333 for key in (_bme.bc.point_coupling, _bme.bc.point_coupling_penalty)
334 ):
335 is_linked_nodes = True
336 mesh.set_node_links()
337 else:
338 is_linked_nodes = False
340 # Dump boundary conditions
341 for (bc_key, geometry_key), bc_list in mesh.boundary_conditions.items():
342 if bc_list:
343 if isinstance(bc_key, str):
344 section = bc_key
345 else:
346 section = _INPUT_FILE_MAPPINGS["boundary_conditions"][
347 (bc_key, geometry_key)
348 ]
350 for boundary_condition in bc_list:
351 if isinstance(boundary_condition, _BoundaryCondition):
352 bc_data = boundary_condition.data
353 elif isinstance(boundary_condition, _Coupling):
354 bc_data = dump_coupling(boundary_condition)
355 else:
356 raise TypeError(
357 f"Got unexpected type {type(boundary_condition)} for boundary condition"
358 )
359 input_file.boundary_conditions[section].append(
360 _FourCBoundaryConditionData(
361 geometry_set_id=input_file.fourc_input.type_converter(
362 boundary_condition.geometry_set
363 ),
364 data=input_file.fourc_input.type_converter(bc_data),
365 )
366 )
368 # If we have previously set the node links, we unlink them here.
369 if is_linked_nodes:
370 mesh.unlink_nodes()
372 # Dump NURBS patch information.
373 for element in mesh.elements:
374 if isinstance(element, _NURBSPatch):
375 dump_nurbs_patch_knotvectors(input_file, element)
378def dump_mesh_representation_to_input_file_yaml(
379 fourc_input: _FourCInput,
380 mesh_representation: _MeshRepresentation,
381 element_type_id_to_data: dict[int, _FourCElementData],
382 boundary_conditions: dict[str, list[_FourCBoundaryConditionData]],
383) -> None:
384 """Dump the information contained in the mesh representation to the 4C input file
385 via FourCIPP, in yaml format.
387 Args:
388 fourc_input: 4C input file via FourCIPP where the mesh information data will be dumped to.
389 mesh_representation: The mesh representation that is added to the input file.
390 element_type_id_to_data: The mapping between element type ID and the element type data.
391 boundary_conditions: The boundary conditions in the input file.
392 """
393 # Compute the starting indices for the nodes and elements entities.
394 start_index_nodes = len(fourc_input.sections.get("NODE COORDS", []))
395 start_index_elements = sum(
396 len(fourc_input.sections.get(section, []))
397 for section in ("FLUID ELEMENTS", "STRUCTURE ELEMENTS")
398 )
400 def _dump(section_name: str, dictionary_list: list):
401 """Append the given list of dictionaries to the section in the input file."""
402 if len(dictionary_list) == 0:
403 return
404 full_item_list = fourc_input.pop(section_name, [])
405 full_item_list.extend(dictionary_list)
406 fourc_input.combine_sections({section_name: full_item_list})
408 # Dump node information to the input file.
409 node_data_list = []
410 for i_node, (point, point_type, cp_weight) in enumerate(
411 zip(
412 mesh_representation.points,
413 mesh_representation.data_iterator("point_data", "point_type"),
414 mesh_representation.data_iterator("point_data", "control_point_weight"),
415 ),
416 start=start_index_nodes,
417 ):
418 node_type = _bme.node_type(point_type)
419 node_id = i_node + 1
420 if node_type == _bme.node_type.node or node_type == _bme.node_type.cosserat:
421 node_data_list.append(
422 {
423 "id": node_id,
424 "COORD": point,
425 "data": {"type": "NODE"},
426 }
427 )
428 elif node_type == _bme.node_type.control_point:
429 node_data_list.append(
430 {
431 "id": node_id,
432 "COORD": point,
433 "data": {
434 "type": "CP",
435 "weight": cp_weight,
436 },
437 }
438 )
439 else:
440 raise ValueError(f"Unknown node type {node_type} for node {i_node + 1}")
441 _dump("NODE COORDS", node_data_list)
443 # Dump element information to the input file.
444 element_list = []
445 for i_element, (
446 connectivity,
447 element_type_id,
448 element_material_id,
449 ) in enumerate(
450 zip(
451 mesh_representation.connectivity_iterator(),
452 mesh_representation.data_iterator("cell_data", "element_type_id"),
453 mesh_representation.data_iterator("cell_data", "material_id"),
454 )
455 ):
456 element_type_data = element_type_id_to_data[int(element_type_id)]
458 node_ordering = _INPUT_FILE_MAPPINGS[
459 "four_c_cell_to_connectivity_mapping_from_vtk"
460 ].get(element_type_data.four_c_cell, None)
461 if node_ordering is not None:
462 connectivity = connectivity[node_ordering]
464 additional_element_data = None
465 if _INPUT_FILE_MAPPINGS["four_c_type_to_requires_triads"].get(
466 element_type_data.four_c_type, False
467 ):
468 # The numpy quaternion package can return rotation vectors outside
469 # the range of -pi to pi, which can cause issues in testing comparisons.
470 # To avoid this, we convert the rotations to BeamMe internal rotations
471 # and extract the rotation vectors again, ensuring they are in the
472 # correct range.
473 # This is super slow, but for now keep it, as a mesh based output is to
474 # be preferred when performance is of importance.
475 if "rotation_vector" not in mesh_representation.point_data:
476 raise KeyError(
477 f"Rotation vectors are required for type {element_type_data.four_c_type}, but "
478 "no rotation vectors found in the mesh representation!"
479 )
480 rotations = [
481 _Rotation.from_rotation_vector(rotation_vector)
482 for rotation_vector in mesh_representation.point_data[
483 "rotation_vector"
484 ][connectivity]
485 ]
486 additional_element_data = {
487 "TRIADS": _np.array(
488 [rotation.get_rotation_vector() for rotation in rotations]
489 ).ravel()
490 }
492 element_list.append(
493 element_type_data.get_yaml_dict(
494 element_id=start_index_elements + i_element,
495 connectivity=start_index_nodes + connectivity,
496 element_material_id=element_material_id,
497 additional_element_data=additional_element_data,
498 )
499 )
500 _dump("STRUCTURE ELEMENTS", element_list)
502 # Dump geometry sets to the input file.
503 # We first create a mapping from the geometry type to a dictionary which maps the
504 # global geometry set ID to the name of the info object. This is required for the
505 # sorting later on.
506 geometry_set_infos = get_geometry_set_infos_with_nodal_flag_vector(
507 mesh_representation
508 )
509 geometry_type_to_geometry_infos: dict[_Geometry, dict[int, _GeometrySetInfo]] = (
510 _defaultdict(dict)
511 )
512 for geometry_set_info in geometry_set_infos:
513 geometry_type_to_geometry_infos[geometry_set_info.geometry_type][
514 geometry_set_info.i_global
515 ] = geometry_set_info
516 for geometry_type, id_to_info in geometry_type_to_geometry_infos.items():
517 geometry_set_list = []
518 # We sort the keys here to ensure that the geometry sets are dumped in the
519 # correct order.
520 sorted_ids = sorted(id_to_info.keys())
521 for i_global in sorted_ids:
522 info = id_to_info[i_global]
523 node_indices = _np.nonzero(info.point_flag_vector)[0]
524 geometry_set_list.extend(
525 [
526 {
527 "type": "NODE",
528 "node_id": start_index_nodes + node_index + 1,
529 "d_type": _INPUT_FILE_MAPPINGS[
530 "geometry_sets_geometry_to_entry_name"
531 ][geometry_type],
532 "d_id": i_global + 1,
533 }
534 for node_index in node_indices
535 ]
536 )
537 _dump(
538 _INPUT_FILE_MAPPINGS["geometry_sets_geometry_to_condition_name"][
539 geometry_type
540 ],
541 geometry_set_list,
542 )
544 # Add the boundary conditions to the input file.
545 for section, input_file_bc in boundary_conditions.items():
546 bc_list = fourc_input.pop(section, [])
547 for boundary_condition in input_file_bc:
548 bc_list.append(boundary_condition.dump_to_input_file_yaml())
549 fourc_input.combine_sections({section: bc_list})
552def dump_mesh_representation_to_input_file_vtu(
553 fourc_input: _FourCInput,
554 mesh_representation: _MeshRepresentation,
555 element_type_id_to_data: dict[int, _FourCElementData],
556 boundary_conditions: dict[str, list[_FourCBoundaryConditionData]],
557) -> _pv.UnstructuredGrid:
558 """Dump the input file to a vtu mesh data format.
560 This function does two things:
561 - Create the vtu file containing the mesh information from the mesh
562 representation.
563 - Dump the information required for a vtu mesh format to the FourCIPP
564 input file.
566 Args:
567 fourc_input: FourCIPP input file where the mesh information data will be
568 dumped to.
569 mesh_representation: The mesh representation that is added to the input file.
570 element_type_id_to_data: The mapping between element type ID and the element
571 type data.
572 boundary_conditions: The boundary conditions in the input file.
574 Returns:
575 The unstructured grid containing the vtu mesh.
576 """
577 # VTU output can not be combined with yaml output.
578 n_yaml_nodes = len(fourc_input.sections.get("NODE COORDS", []))
579 n_yaml_elements = sum(
580 len(fourc_input.sections.get(section, []))
581 for section in ("FLUID ELEMENTS", "STRUCTURE ELEMENTS")
582 )
583 if n_yaml_nodes > 0 or n_yaml_elements > 0:
584 raise ValueError(
585 "Mesh output in `vtu` format is not possible if there are yaml nodes or "
586 "elements in the input file."
587 )
589 # VTU output can not be combined with existing STRUCTURE GEOMETRY section.
590 if "STRUCTURE GEOMETRY" in fourc_input:
591 raise ValueError(
592 "Mesh output in `vtu` format is not possible if there is already a "
593 "STRUCTURE GEOMETRY section in the input file."
594 )
596 # Get a pyvista grid representing the mesh representation
597 grid = mesh_representation.get_pyvista_grid()
599 # Check which element types need triads.
600 element_type_ids_with_triads: set[int] = set()
601 triad_n_nodes: set[int] = set()
602 requires_triads = _INPUT_FILE_MAPPINGS["four_c_type_to_requires_triads"]
603 cell_to_type_and_n_nodes = _INPUT_FILE_MAPPINGS[
604 "four_c_cell_to_element_type_and_n_nodes"
605 ]
606 for (
607 i_element_type,
608 element_data,
609 ) in element_type_id_to_data.items():
610 if requires_triads.get(element_data.four_c_type, False):
611 n_nodes = cell_to_type_and_n_nodes[element_data.four_c_cell][1]
612 element_type_ids_with_triads.add(i_element_type)
613 triad_n_nodes.add(n_nodes)
615 # Get the indices of the elements where we need triad information.
616 triad_info_indices = _np.flatnonzero(
617 _np.isin(
618 mesh_representation.cell_data["element_type_id"],
619 list(element_type_ids_with_triads),
620 )
621 )
623 # Extract the triad information for each element where we need it.
624 rotation_vectors = mesh_representation.point_data["rotation_vector"]
625 element_rotation_vectors = {
626 n_nodes: _np.zeros((mesh_representation.n_cells, 3 * n_nodes))
627 for n_nodes in triad_n_nodes
628 }
629 beam_connectivity_mapping = _INPUT_FILE_MAPPINGS["beam_vtk_mapping_to_four_c"]
630 for i_element, connectivity in zip(
631 triad_info_indices,
632 mesh_representation.connectivity_iterator(element_indices=triad_info_indices),
633 ):
634 # The order of the element rotation vectors has to be consistent with the
635 # connectivity in 4C, not the one we write to vtu. Thus, we have to reorder
636 # the connectivity here.
637 n_nodes = len(connectivity)
638 element_rotation_vectors[n_nodes][i_element] = rotation_vectors[
639 connectivity[beam_connectivity_mapping[n_nodes]]
640 ].ravel()
641 for n_nodes, element_rotation_vector in element_rotation_vectors.items():
642 grid.cell_data[f"element_rotation_vector_{n_nodes}"] = element_rotation_vector
644 # Get the block IDs for each cell (elements with the same material and element data are in the same block).
645 element_type_and_material_to_block_id: dict[tuple[int, int], int] = {}
646 cell_block_id = _np.empty(mesh_representation.n_cells, dtype=int)
647 for i_cell, (element_type_id, material_id) in enumerate(
648 zip(
649 mesh_representation.data_iterator("cell_data", "element_type_id"),
650 mesh_representation.data_iterator("cell_data", "material_id"),
651 )
652 ):
653 key = (int(element_type_id), int(material_id))
654 if key not in element_type_and_material_to_block_id:
655 element_type_and_material_to_block_id[key] = len(
656 element_type_and_material_to_block_id
657 )
658 cell_block_id[i_cell] = element_type_and_material_to_block_id[key]
659 grid.cell_data["block_id"] = cell_block_id
661 # Add the element blocks to the input file.
662 cell_block_id_to_element_type_and_material = _create_inverse_mapping(
663 element_type_and_material_to_block_id
664 )
665 element_blocks = []
666 for i_block in sorted(cell_block_id_to_element_type_and_material.keys()):
667 element_type_id, material_id = cell_block_id_to_element_type_and_material[
668 i_block
669 ]
670 element_data = element_type_id_to_data[element_type_id]
671 additional_data = None
672 has_triads = requires_triads.get(element_data.four_c_type, False)
673 if has_triads:
674 n_nodes = cell_to_type_and_n_nodes[element_data.four_c_cell][1]
675 additional_data = {
676 "NODAL_ROTATION_VECTORS": (f"element_rotation_vector_{n_nodes}")
677 }
678 element_blocks.append(
679 element_data.get_block_dict(i_block, material_id, additional_data)
680 )
681 fourc_input.combine_sections(
682 {
683 "STRUCTURE GEOMETRY": {
684 "FILE": "TO_BE_DEFINED",
685 "SHOW_INFO": "detailed_summary",
686 "ELEMENT_BLOCKS": element_blocks,
687 }
688 }
689 )
691 # Add geometry sets from the mesh representation to the VTU grid.
692 geometry_set_infos = get_geometry_set_infos_with_nodal_flag_vector(
693 mesh_representation
694 )
695 geometry_sets_in_mr: dict[int, _GeometrySetInfo] = {}
696 for geometry_set_info in geometry_set_infos:
697 geometry_set_id = geometry_set_info.i_global
698 geometry_sets_in_mr[geometry_set_id] = geometry_set_info
699 if geometry_set_info.name is not None:
700 data_array_name = geometry_set_info.name
701 else:
702 data_array_name = f"point_set_{geometry_set_id}"
703 if data_array_name in grid.point_data:
704 raise ValueError(
705 f"Data array name {data_array_name} for geometry set {geometry_set_id} "
706 "already exists in the VTU grid. Please ensure that the geometry set names "
707 "in the mesh representation do not conflict with existing data array names."
708 )
709 grid.point_data[data_array_name] = geometry_set_info.point_flag_vector
711 # Add definition of the boundary conditions (geometry sets defined in the vtu
712 # mesh) to the yaml input file.
713 for section, input_file_bc in boundary_conditions.items():
714 bc_list = fourc_input.pop(section, [])
715 for boundary_condition in input_file_bc:
716 bc_list.append(
717 boundary_condition.dump_to_input_file_vtu(geometry_sets_in_mr)
718 )
719 fourc_input.combine_sections({section: bc_list})
721 # Return the vtu grid.
722 return grid