Coverage for src/beamme/space_time/beam_to_space_time.py: 96%
148 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"""Convert a beam to a space time surface mesh."""
24from collections.abc import Callable as _Callable
25from typing import cast as _cast
27import numpy as _np
28import pyvista as _pv
30from beamme.core.conf import bme as _bme
31from beamme.core.coupling import Coupling as _Coupling
32from beamme.core.element_volume import VolumeElement as _VolumeElement
33from beamme.core.geometry_set import GeometryName as _GeometryName
34from beamme.core.geometry_set import GeometrySet as _GeometrySet
35from beamme.core.geometry_set import GeometrySetBase as _GeometrySetBase
36from beamme.core.geometry_set import GeometrySetNodes as _GeometrySetNodes
37from beamme.core.mesh import Mesh as _Mesh
38from beamme.core.mesh_representation import MeshRepresentation as _MeshRepresentation
39from beamme.core.mesh_utils import (
40 apply_nodal_coupling_to_mesh_representation as _apply_nodal_coupling_to_mesh_representation,
41)
42from beamme.core.node import Node as _Node
43from beamme.core.node import NodeCosserat as _NodeCosserat
46class NodeCosseratSpaceTime(_NodeCosserat):
47 """A Cosserat node in space-time.
49 We add the 4th dimension time as a class variable.
50 """
52 node_type = _bme.node_type.space_time_cosserat
54 def __init__(self, coordinates, rotation, time, **kwargs):
55 super().__init__(coordinates, rotation, **kwargs)
56 self.time = time
59class SpaceTimeElement(_VolumeElement):
60 """A general beam space-time surface element."""
62 def __init__(self, nodes, **kwargs):
63 super().__init__(nodes=nodes, **kwargs)
66class SpaceTimeElementQuad4(SpaceTimeElement):
67 """A space-time element with 4 nodes."""
69 element_type = _bme.element_type.space_time_beam
70 vtk_cell_type = _pv.CellType.QUAD
71 data = {}
74class SpaceTimeElementQuad9(SpaceTimeElement):
75 """A space-time element with 9 nodes."""
77 element_type = _bme.element_type.space_time_beam
78 vtk_cell_type = _pv.CellType.BIQUADRATIC_QUAD
79 data = {}
82def beam_to_space_time(
83 mesh_space_or_generator: _Mesh | _Callable[[float], _Mesh],
84 time_duration: float,
85 number_of_elements_in_time: int,
86 *,
87 time_start: float = 0.0,
88) -> tuple[_Mesh, _GeometryName]:
89 """Convert a beam mesh to a surface space-time mesh.
91 Args:
92 mesh_space_or_generator:
93 Either a fixed spatial Mesh object or a function that returns the
94 spatial mesh for a given time. If this is a generator, the topology
95 of the mesh at the initial time is chosen for all times, only the
96 positions and rotations are updated.
97 time_duration:
98 Total time increment to be solved with the space-time mesh
99 number_of_elements_in_time:
100 Number of elements in time direction
101 time_start:
102 Starting time for the space-time mesh. Can be used to create time slaps.
103 Returns:
104 Tuple (space_time_mesh, return_set)
105 - space_time_mesh:
106 The space time mesh. Be aware that translating / rotating this mesh
107 might lead to unexpected results.
108 - return_set:
109 The nodes sets to be returned for the space time mesh:
110 "start", "end", "surface"
111 """
112 # Get the "reference" spatial mesh
113 if callable(mesh_space_or_generator):
114 mesh_space_reference = mesh_space_or_generator(time_start)
115 else:
116 mesh_space_reference = mesh_space_or_generator
118 # Perform some sanity checks
119 element_types = {type(element) for element in mesh_space_reference.elements}
120 if not len(element_types) == 1:
121 raise ValueError(
122 f"Expected all elements to be of the same type, got {element_types}"
123 )
124 element_type = element_types.pop()
126 # Calculate global mesh properties
127 number_of_nodes_in_space = len(mesh_space_reference.nodes)
128 number_of_elements_in_space = len(mesh_space_reference.elements)
129 space_time_element_type: type[SpaceTimeElementQuad4] | type[SpaceTimeElementQuad9]
131 if len(element_type.nodes_create) == 2:
132 number_of_copies_in_time = number_of_elements_in_time + 1
133 time_increment_between_nodes = time_duration / number_of_elements_in_time
134 space_time_element_type = SpaceTimeElementQuad4
135 elif len(element_type.nodes_create) == 3:
136 number_of_copies_in_time = 2 * number_of_elements_in_time + 1
137 time_increment_between_nodes = time_duration / (2 * number_of_elements_in_time)
138 space_time_element_type = SpaceTimeElementQuad9
139 else:
140 raise TypeError(f"Got unexpected element type {element_type}")
142 # Number nodes and elements in the original mesh
143 for i_node, node in enumerate(mesh_space_reference.nodes):
144 node.i_global = i_node
145 for i_element, element in enumerate(mesh_space_reference.elements):
146 element.i_global = i_element
148 # Get the nodes for the final space-time mesh
149 space_time_nodes = []
150 start_nodes: list[_Node] = []
151 end_nodes: list[_Node] = []
152 for i_mesh_space in range(number_of_copies_in_time):
153 time = time_increment_between_nodes * i_mesh_space + time_start
155 if callable(mesh_space_or_generator):
156 mesh_space_current_time = mesh_space_or_generator(time)
157 if (not len(mesh_space_current_time.nodes) == number_of_nodes_in_space) or (
158 not len(mesh_space_current_time.elements) == number_of_elements_in_space
159 ):
160 raise ValueError(
161 "The number of nodes and elements does not match for the generated "
162 "space time meshes."
163 )
164 else:
165 mesh_space_current_time = mesh_space_reference
167 space_time_nodes_to_add = [
168 NodeCosseratSpaceTime(
169 node.coordinates, node.rotation, time, arc_length=node.arc_length
170 )
171 for node in mesh_space_current_time.nodes
172 ]
173 space_time_nodes.extend(space_time_nodes_to_add)
175 if i_mesh_space == 0:
176 start_nodes.extend(space_time_nodes_to_add)
177 elif i_mesh_space == number_of_copies_in_time - 1:
178 end_nodes.extend(space_time_nodes_to_add)
180 # Create the space time elements
181 space_time_elements = []
182 for i_element_time in range(number_of_elements_in_time):
183 for element in mesh_space_reference.elements:
184 element_node_ids = [node.i_global for node in element.nodes]
185 if space_time_element_type == SpaceTimeElementQuad4:
186 # Create the indices for the linear element
187 first_time_row_start_index = i_element_time * number_of_nodes_in_space
188 second_time_row_start_index = (
189 1 + i_element_time
190 ) * number_of_nodes_in_space
191 element_node_indices = [
192 first_time_row_start_index + element_node_ids[0],
193 first_time_row_start_index + element_node_ids[1],
194 second_time_row_start_index + element_node_ids[1],
195 second_time_row_start_index + element_node_ids[0],
196 ]
197 elif space_time_element_type == SpaceTimeElementQuad9:
198 # Create the indices for the quadratic element
199 first_time_row_start_index = (
200 2 * i_element_time * number_of_nodes_in_space
201 )
202 second_time_row_start_index = (
203 2 * i_element_time + 1
204 ) * number_of_nodes_in_space
205 third_time_row_start_index = (
206 2 * i_element_time + 2
207 ) * number_of_nodes_in_space
208 element_node_indices = [
209 first_time_row_start_index + element_node_ids[0],
210 first_time_row_start_index + element_node_ids[2],
211 third_time_row_start_index + element_node_ids[2],
212 third_time_row_start_index + element_node_ids[0],
213 first_time_row_start_index + element_node_ids[1],
214 second_time_row_start_index + element_node_ids[2],
215 third_time_row_start_index + element_node_ids[1],
216 second_time_row_start_index + element_node_ids[0],
217 second_time_row_start_index + element_node_ids[1],
218 ]
219 else:
220 raise TypeError(
221 f"Got unexpected space time element type {space_time_element_type}"
222 )
224 # Add the element to the mesh
225 space_time_elements.append(
226 space_time_element_type(
227 [space_time_nodes[i_node] for i_node in element_node_indices]
228 )
229 )
231 # Add joints to the space time mesh
232 space_time_couplings = []
233 coupling_geometry_sets = set()
234 for coupling in mesh_space_reference.boundary_conditions[
235 _bme.bc.point_coupling, _bme.geo.point
236 ]:
237 coupling_set = coupling.geometry_set
238 coupling_geometry_sets.add(coupling_set)
239 coupling_node_ids = [node.i_global for node in coupling_set.get_points()]
240 for i_mesh_space in range(number_of_copies_in_time):
241 space_time_couplings.append(
242 _Coupling(
243 [
244 space_time_nodes[
245 node_id + i_mesh_space * number_of_nodes_in_space
246 ]
247 for node_id in coupling_node_ids
248 ],
249 coupling.bc_type,
250 coupling.data,
251 )
252 )
254 # Convert geometry sets to the space time mesh
255 raise_geometry_type = {
256 _bme.geo.point: _bme.geo.line,
257 _bme.geo.line: _bme.geo.surface,
258 _bme.geo.surface: _bme.geo.volume,
259 }
260 all_sets_in_space = mesh_space_reference.get_unique_geometry_sets()
261 space_time_geometry_sets: list[_GeometrySetBase] = []
262 for geometry_type, geometry_sets in all_sets_in_space.items():
263 for geometry_set in geometry_sets:
264 if geometry_set in coupling_geometry_sets:
265 # The coupling geometry sets are already handled above, so we skip them here.
266 continue
268 if isinstance(geometry_set, _GeometrySet) and (
269 geometry_type == _bme.geo.line
270 or geometry_type == _bme.geo.surface
271 or geometry_type == _bme.geo.volume
272 ):
273 raised_geometry_set_elements = []
274 for element in geometry_set.get_geometry_objects():
275 for i_element_row_in_time in range(number_of_elements_in_time):
276 raised_geometry_set_elements.append(
277 space_time_elements[
278 _cast(int, element.i_global)
279 + i_element_row_in_time * number_of_elements_in_space
280 ]
281 )
282 space_time_geometry_sets.append(
283 _GeometrySet(
284 raised_geometry_set_elements,
285 name=geometry_set.name,
286 )
287 )
289 else:
290 geometry_set_nodes = geometry_set.get_all_nodes()
291 raised_geometry_set_nodes = []
292 for node in geometry_set_nodes:
293 for i_mesh_space in range(number_of_copies_in_time):
294 raised_geometry_set_nodes.append(
295 space_time_nodes[
296 node.i_global + i_mesh_space * number_of_nodes_in_space
297 ]
298 )
300 geometry_type_raised = raise_geometry_type[geometry_type]
301 space_time_geometry_sets.append(
302 _GeometrySetNodes(
303 geometry_type_raised,
304 raised_geometry_set_nodes,
305 name=geometry_set.name,
306 )
307 )
309 # Create the new mesh and add all the mesh items
310 space_time_mesh = _Mesh()
311 space_time_mesh.add(space_time_nodes)
312 space_time_mesh.add(space_time_elements)
313 space_time_mesh.add(space_time_couplings)
314 space_time_mesh.add(space_time_geometry_sets)
316 # Create the element sets
317 return_set = _GeometryName()
318 return_set["start"] = _GeometrySetNodes(_bme.geo.line, start_nodes)
319 return_set["end"] = _GeometrySetNodes(_bme.geo.line, end_nodes)
320 return_set["surface"] = _GeometrySetNodes(_bme.geo.surface, space_time_mesh.nodes)
322 return space_time_mesh, return_set
325def get_space_time_mesh_representation(mesh: _Mesh) -> _MeshRepresentation:
326 """Get the mesh representation for the space time mesh.
328 Compared to the standard mesh representation, coupled nodes are represented by the
329 same node. This requires some additional element data arrays which are added by
330 this function.
332 Args:
333 mesh: The space time mesh.
335 Returns:
336 The mesh representation for the space time mesh.
337 """
338 element_types = list(set([type(element) for element in mesh.elements]))
339 if len(element_types) > 1:
340 raise ValueError("Got more than a single element type, this is not supported")
341 elif not (
342 element_types[0] == SpaceTimeElementQuad4
343 or element_types[0] == SpaceTimeElementQuad9
344 ):
345 raise TypeError(
346 f"Expected either SpaceTimeElementQuad4 or SpaceTimeElementQuad9, got {element_types[0]}"
347 )
349 # Number of nodes per element
350 n_nodes_per_element = len(mesh.elements[0].nodes)
352 # Get the mesh representation
353 (mesh_representation, _, geometry_sets_to_i_global, _) = (
354 mesh.get_mesh_representation()
355 )
357 # Get the element rotation vectors and arc length values. This has to be done before
358 # the coupled nodes are removed.
359 point_rotation_vectors = mesh_representation.point_data["rotation_vector"]
360 element_rotation_vectors = _np.zeros(
361 (mesh_representation.n_cells, n_nodes_per_element * 3)
362 )
363 element_arc_lengths = None
364 if "arc_length" in mesh_representation.point_data:
365 point_arc_lengths = mesh_representation.point_data["arc_length"]
366 element_arc_lengths = _np.zeros(
367 (mesh_representation.n_cells, n_nodes_per_element)
368 )
370 for i_element, connectivity in enumerate(
371 mesh_representation.connectivity_iterator()
372 ):
373 for i_local, i_global in enumerate(connectivity):
374 element_rotation_vectors[i_element, i_local * 3 : (i_local + 1) * 3] = (
375 point_rotation_vectors[i_global]
376 )
377 if element_arc_lengths is not None:
378 element_arc_lengths[i_element, i_local] = point_arc_lengths[i_global]
380 # Add the element data arrays
381 mesh_representation.cell_data["rotation_vector"] = element_rotation_vectors
382 if element_arc_lengths is not None:
383 mesh_representation.cell_data["arc_length"] = element_arc_lengths
385 # Apply the coupling by explicitly replacing the coupled nodes.
386 _apply_nodal_coupling_to_mesh_representation(
387 mesh_representation,
388 geometry_sets_to_i_global,
389 mesh.boundary_conditions[_bme.bc.point_coupling, _bme.geo.point],
390 )
392 return mesh_representation