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