Coverage for src/beamme/space_time/beam_to_space_time.py: 95%
147 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"""Convert a beam to a space time surface mesh."""
24from typing import Callable as _Callable
25from typing import Dict as _Dict
26from typing import List as _List
27from typing import Tuple as _Tuple
28from typing import Type as _Type
29from typing import cast as _cast
31import numpy as _np
32import pyvista as _pv
33import vtk as _vtk
35from beamme.core.conf import bme as _bme
36from beamme.core.coupling import Coupling as _Coupling
37from beamme.core.element_volume import VolumeElement as _VolumeElement
38from beamme.core.geometry_set import GeometryName as _GeometryName
39from beamme.core.geometry_set import GeometrySet as _GeometrySet
40from beamme.core.geometry_set import GeometrySetNodes as _GeometrySetNodes
41from beamme.core.mesh import Mesh as _Mesh
42from beamme.core.mesh_utils import (
43 get_coupled_nodes_to_master_map as _get_coupled_nodes_to_master_map,
44)
45from beamme.core.node import NodeCosserat as _NodeCosserat
46from beamme.utils.nodes import get_nodal_coordinates as _get_nodal_coordinates
49class NodeCosseratSpaceTime(_NodeCosserat):
50 """A Cosserat node in space-time.
52 We add the 4th dimension time as a class variable.
53 """
55 node_type = _bme.node_type.space_time_cosserat
57 def __init__(self, coordinates, rotation, time, **kwargs):
58 super().__init__(coordinates, rotation, **kwargs)
59 self.time = time
62class SpaceTimeElement(_VolumeElement):
63 """A general beam space-time surface element."""
65 def __init__(self, nodes, **kwargs):
66 super().__init__(nodes=nodes, **kwargs)
69class SpaceTimeElementQuad4(SpaceTimeElement):
70 """A space-time element with 4 nodes."""
72 vtk_cell_type = _pv.CellType.QUAD
73 vtk_cell_type_legacy = _vtk.vtkQuad
74 vtk_topology = list(range(4))
77class SpaceTimeElementQuad9(SpaceTimeElement):
78 """A space-time element with 9 nodes."""
80 vtk_cell_type = _pv.CellType.BIQUADRATIC_QUAD
81 vtk_cell_type_legacy = _vtk.vtkQuadraticQuad
82 vtk_topology = list(range(9))
85def beam_to_space_time(
86 mesh_space_or_generator: _Mesh | _Callable[[float], _Mesh],
87 time_duration: float,
88 number_of_elements_in_time: int,
89 *,
90 time_start: float = 0.0,
91) -> _Tuple[_Mesh, _GeometryName]:
92 """Convert a beam mesh to a surface space-time mesh.
94 Args:
95 mesh_space_or_generator:
96 Either a fixed spatial Mesh object or a function that returns the
97 spatial mesh for a given time. If this is a generator, the topology
98 of the mesh at the initial time is chosen for all times, only the
99 positions and rotations are updated.
100 time_duration:
101 Total time increment to be solved with the space-time mesh
102 number_of_elements_in_time:
103 Number of elements in time direction
104 time_start:
105 Starting time for the space-time mesh. Can be used to create time slaps.
106 Returns:
107 Tuple (space_time_mesh, return_set)
108 - space_time_mesh:
109 The space time mesh. Be aware that translating / rotating this mesh
110 might lead to unexpected results.
111 - return_set:
112 The nodes sets to be returned for the space time mesh:
113 "start", "end", "left", "right", "surface"
114 """
116 # Get the "reference" spatial mesh
117 if callable(mesh_space_or_generator):
118 mesh_space_reference = mesh_space_or_generator(time_start)
119 else:
120 mesh_space_reference = mesh_space_or_generator
122 # Perform some sanity checks
123 element_types = {type(element) for element in mesh_space_reference.elements}
124 if not len(element_types) == 1:
125 raise ValueError(
126 f"Expected all elements to be of the same type, got {element_types}"
127 )
128 element_type = element_types.pop()
130 # Calculate global mesh properties
131 number_of_nodes_in_space = len(mesh_space_reference.nodes)
132 number_of_elements_in_space = len(mesh_space_reference.elements)
133 space_time_element_type: _Type[SpaceTimeElementQuad4] | _Type[SpaceTimeElementQuad9]
135 if len(element_type.nodes_create) == 2:
136 number_of_copies_in_time = number_of_elements_in_time + 1
137 time_increment_between_nodes = time_duration / number_of_elements_in_time
138 space_time_element_type = SpaceTimeElementQuad4
139 elif len(element_type.nodes_create) == 3:
140 number_of_copies_in_time = 2 * number_of_elements_in_time + 1
141 time_increment_between_nodes = time_duration / (2 * number_of_elements_in_time)
142 space_time_element_type = SpaceTimeElementQuad9
143 else:
144 raise TypeError(f"Got unexpected element type {element_type}")
146 # Number the nodes in the original mesh
147 for i_node, node in enumerate(mesh_space_reference.nodes):
148 node.i_global = i_node
150 # Get the nodes for the final space-time mesh
151 left_nodes = []
152 right_nodes = []
153 space_time_nodes = []
154 for i_mesh_space in range(number_of_copies_in_time):
155 time = time_increment_between_nodes * i_mesh_space + time_start
157 if callable(mesh_space_or_generator):
158 mesh_space_current_time = mesh_space_or_generator(time)
159 if (not len(mesh_space_current_time.nodes) == number_of_nodes_in_space) or (
160 not len(mesh_space_current_time.elements) == number_of_elements_in_space
161 ):
162 raise ValueError(
163 "The number of nodes and elements does not match for the generated "
164 "space time meshes."
165 )
166 else:
167 mesh_space_current_time = mesh_space_reference
169 space_time_nodes_to_add = [
170 NodeCosseratSpaceTime(
171 node.coordinates, node.rotation, time, arc_length=node.arc_length
172 )
173 for node in mesh_space_current_time.nodes
174 ]
175 space_time_nodes.extend(space_time_nodes_to_add)
177 if i_mesh_space == 0:
178 start_nodes = space_time_nodes_to_add
179 elif i_mesh_space == number_of_copies_in_time - 1:
180 end_nodes = space_time_nodes_to_add
182 left_nodes.append(space_time_nodes_to_add[0])
183 right_nodes.append(space_time_nodes_to_add[-1])
185 # Create the space time elements
186 space_time_elements = []
187 for i_element_time in range(number_of_elements_in_time):
188 for element in mesh_space_reference.elements:
189 element_node_ids = [node.i_global for node in element.nodes]
190 if space_time_element_type == SpaceTimeElementQuad4:
191 # Create the indices for the linear element
192 first_time_row_start_index = i_element_time * number_of_nodes_in_space
193 second_time_row_start_index = (
194 1 + i_element_time
195 ) * number_of_nodes_in_space
196 element_node_indices = [
197 first_time_row_start_index + element_node_ids[0],
198 first_time_row_start_index + element_node_ids[1],
199 second_time_row_start_index + element_node_ids[1],
200 second_time_row_start_index + element_node_ids[0],
201 ]
202 elif space_time_element_type == SpaceTimeElementQuad9:
203 # Create the indices for the quadratic element
204 first_time_row_start_index = (
205 2 * i_element_time * number_of_nodes_in_space
206 )
207 second_time_row_start_index = (
208 2 * i_element_time + 1
209 ) * number_of_nodes_in_space
210 third_time_row_start_index = (
211 2 * i_element_time + 2
212 ) * number_of_nodes_in_space
213 element_node_indices = [
214 first_time_row_start_index + element_node_ids[0],
215 first_time_row_start_index + element_node_ids[2],
216 third_time_row_start_index + element_node_ids[2],
217 third_time_row_start_index + element_node_ids[0],
218 first_time_row_start_index + element_node_ids[1],
219 second_time_row_start_index + element_node_ids[2],
220 third_time_row_start_index + element_node_ids[1],
221 second_time_row_start_index + element_node_ids[0],
222 second_time_row_start_index + element_node_ids[1],
223 ]
224 else:
225 raise TypeError(
226 f"Got unexpected space time element type {space_time_element_type}"
227 )
229 # Add the element to the mesh
230 space_time_elements.append(
231 space_time_element_type(
232 [space_time_nodes[i_node] for i_node in element_node_indices]
233 )
234 )
236 # Add joints to the space time mesh
237 space_time_couplings = []
238 for coupling in mesh_space_reference.boundary_conditions[
239 _bme.bc.point_coupling, _bme.geo.point
240 ]:
241 coupling_node_ids = [
242 node.i_global for node in coupling.geometry_set.get_points()
243 ]
244 for i_mesh_space in range(number_of_copies_in_time):
245 space_time_couplings.append(
246 _Coupling(
247 [
248 space_time_nodes[
249 node_id + i_mesh_space * number_of_nodes_in_space
250 ]
251 for node_id in coupling_node_ids
252 ],
253 coupling.bc_type,
254 coupling.data,
255 )
256 )
258 # Create the new mesh and add all the mesh items
259 space_time_mesh = _Mesh()
260 space_time_mesh.add(space_time_nodes)
261 space_time_mesh.add(space_time_elements)
262 space_time_mesh.add(space_time_couplings)
264 # Create the element sets
265 return_set = _GeometryName()
266 return_set["start"] = _GeometrySet(start_nodes)
267 return_set["end"] = _GeometrySet(end_nodes)
268 return_set["left"] = _GeometrySet(left_nodes)
269 return_set["right"] = _GeometrySet(right_nodes)
270 return_set["surface"] = _GeometrySetNodes(_bme.geo.surface, space_time_mesh.nodes)
272 return space_time_mesh, return_set
275def mesh_to_data_arrays(mesh: _Mesh):
276 """Get the relevant data arrays from the space time mesh."""
278 element_types = list(set([type(element) for element in mesh.elements]))
279 if len(element_types) > 1:
280 raise ValueError("Got more than a single element type, this is not supported")
281 elif not (
282 element_types[0] == SpaceTimeElementQuad4
283 or element_types[0] == SpaceTimeElementQuad9
284 ):
285 raise TypeError(
286 f"Expected either SpaceTimeElementQuad4 or SpaceTimeElementQuad9, got {element_types[0]}"
287 )
289 _, unique_nodes = _get_coupled_nodes_to_master_map(mesh, assign_i_global=True)
291 n_nodes = len(unique_nodes)
292 n_elements = len(mesh.elements)
293 n_nodes_per_element = len(mesh.elements[0].nodes)
295 coordinates = _get_nodal_coordinates(unique_nodes)
296 time = _np.zeros(n_nodes)
297 connectivity = _np.zeros((n_elements, n_nodes_per_element), dtype=int)
298 element_rotation_vectors = _np.zeros((n_elements, n_nodes_per_element, 3))
300 unique_nodes_casted_space_time = _cast(_List[NodeCosseratSpaceTime], unique_nodes)
301 for i_node, node in enumerate(unique_nodes_casted_space_time):
302 time[i_node] = node.time
304 for i_element, element in enumerate(mesh.elements):
305 for i_node, node in enumerate(element.nodes):
306 connectivity[i_element, i_node] = node.i_global
307 element_rotation_vectors[i_element, i_node, :] = (
308 node.rotation.get_rotation_vector()
309 )
311 geometry_sets = mesh.get_unique_geometry_sets()
312 node_sets: _Dict[str, _Dict] = {}
313 for value in geometry_sets.values():
314 for geometry_set in value:
315 node_ids = sorted(
316 list(set(node.i_global for node in geometry_set.get_all_nodes()))
317 )
318 node_set_data = {"node_ids": node_ids}
319 if geometry_set.name is not None:
320 node_set_data["name"] = geometry_set.name
321 node_sets[str(len(node_sets) + 1)] = node_set_data
323 return_dict = {
324 "coordinates": coordinates,
325 "time": time,
326 "connectivity": connectivity,
327 "element_rotation_vectors": element_rotation_vectors,
328 "node_sets": node_sets,
329 }
331 nodes_have_arc_length = {node.arc_length is not None for node in mesh.nodes}
332 if len(nodes_have_arc_length) > 1:
333 raise ValueError(
334 "Some nodes have an arc length, some don't. This is not supported."
335 )
336 if nodes_have_arc_length.pop():
337 # The arc length is added as an "element" property, since the same
338 # node can have a different arc length depending on the element
339 # (similar to the rotation vectors).
340 arc_length = _np.zeros((n_elements, n_nodes_per_element))
341 for i_element, element in enumerate(mesh.elements):
342 for i_node, node in enumerate(element.nodes):
343 connectivity[i_element, i_node] = node.i_global
344 arc_length[i_element, i_node] = node.arc_length
345 return_dict["arc_length"] = arc_length
347 return return_dict