Coverage for src/meshpy/space_time/beam_to_space_time.py: 95%
140 statements
« prev ^ index » next coverage.py v7.9.0, created at 2025-06-13 04:26 +0000
1# The MIT License (MIT)
2#
3# Copyright (c) 2018-2025 MeshPy 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 MeshPy 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 Union as _Union
30from typing import cast as _cast
32import numpy as _np
33import vtk as _vtk
35from meshpy.core.conf import mpy as _mpy
36from meshpy.core.coupling import Coupling as _Coupling
37from meshpy.core.element_volume import VolumeElement as _VolumeElement
38from meshpy.core.geometry_set import GeometryName as _GeometryName
39from meshpy.core.geometry_set import GeometrySet as _GeometrySet
40from meshpy.core.geometry_set import GeometrySetNodes as _GeometrySetNodes
41from meshpy.core.mesh import Mesh as _Mesh
42from meshpy.core.mesh_utils import (
43 get_coupled_nodes_to_master_map as _get_coupled_nodes_to_master_map,
44)
45from meshpy.core.node import NodeCosserat as _NodeCosserat
46from meshpy.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 def __init__(self, coordinates, rotation, time, **kwargs):
56 super().__init__(coordinates, rotation, **kwargs)
57 self.time = time
60class SpaceTimeElement(_VolumeElement):
61 """A general beam space-time surface element."""
63 def __init__(self, nodes, **kwargs):
64 super().__init__(nodes=nodes, **kwargs)
67class SpaceTimeElementQuad4(SpaceTimeElement):
68 """A space-time element with 4 nodes."""
70 vtk_cell_type = _vtk.vtkQuad
71 vtk_topology = list(range(4))
74class SpaceTimeElementQuad9(SpaceTimeElement):
75 """A space-time element with 9 nodes."""
77 vtk_cell_type = _vtk.vtkQuadraticQuad
78 vtk_topology = list(range(9))
81def beam_to_space_time(
82 mesh_space_or_generator: _Union[_Mesh, _Callable[[float], _Mesh]],
83 time_duration: float,
84 number_of_elements_in_time: int,
85 *,
86 time_start: float = 0.0,
87) -> _Tuple[_Mesh, _GeometryName]:
88 """Convert a MeshPy beam mesh to a surface space-time mesh.
90 Args:
91 mesh_space_or_generator:
92 Either a fixed spatial Mesh object or a function that returns the
93 spatial mesh for a given time. If this is a generator, the topology
94 of the mesh at the initial time is chosen for all times, only the
95 positions and rotations are updated.
96 time_duration:
97 Total time increment to be solved with the space-time mesh
98 number_of_elements_in_time:
99 Number of elements in time direction
100 time_start:
101 Starting time for the space-time mesh. Can be used to create time slaps.
102 Returns:
103 Tuple (space_time_mesh, return_set)
104 - space_time_mesh:
105 The space time mesh. Be aware that translating / rotating this mesh
106 might lead to unexpected results.
107 - return_set:
108 The nodes sets to be returned for the space time mesh:
109 "start", "end", "left", "right", "surface"
110 """
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: _Union[
130 _Type[SpaceTimeElementQuad4], _Type[SpaceTimeElementQuad9]
131 ]
132 if len(element_type.nodes_create) == 2:
133 number_of_copies_in_time = number_of_elements_in_time + 1
134 time_increment_between_nodes = time_duration / number_of_elements_in_time
135 space_time_element_type = SpaceTimeElementQuad4
136 elif len(element_type.nodes_create) == 3:
137 number_of_copies_in_time = 2 * number_of_elements_in_time + 1
138 time_increment_between_nodes = time_duration / (2 * number_of_elements_in_time)
139 space_time_element_type = SpaceTimeElementQuad9
140 else:
141 raise TypeError(f"Got unexpected element type {element_type}")
143 # Number the nodes in the original mesh
144 for i_node, node in enumerate(mesh_space_reference.nodes):
145 node.i_global = i_node
147 # Get the nodes for the final space-time mesh
148 left_nodes = []
149 right_nodes = []
150 space_time_nodes = []
151 for i_mesh_space in range(number_of_copies_in_time):
152 time = time_increment_between_nodes * i_mesh_space + time_start
154 if callable(mesh_space_or_generator):
155 mesh_space_current_time = mesh_space_or_generator(time)
156 if (not len(mesh_space_current_time.nodes) == number_of_nodes_in_space) or (
157 not len(mesh_space_current_time.elements) == number_of_elements_in_space
158 ):
159 raise ValueError(
160 "The number of nodes and elements does not match for the generated "
161 "space time meshes."
162 )
163 else:
164 mesh_space_current_time = mesh_space_reference
166 space_time_nodes_to_add = [
167 NodeCosseratSpaceTime(
168 node.coordinates, node.rotation, time, arc_length=node.arc_length
169 )
170 for node in mesh_space_current_time.nodes
171 ]
172 space_time_nodes.extend(space_time_nodes_to_add)
174 if i_mesh_space == 0:
175 start_nodes = space_time_nodes_to_add
176 elif i_mesh_space == number_of_copies_in_time - 1:
177 end_nodes = space_time_nodes_to_add
179 left_nodes.append(space_time_nodes_to_add[0])
180 right_nodes.append(space_time_nodes_to_add[-1])
182 # Create the space time elements
183 space_time_elements = []
184 for i_element_time in range(number_of_elements_in_time):
185 for element in mesh_space_reference.elements:
186 element_node_ids = [node.i_global for node in element.nodes]
187 if space_time_element_type == SpaceTimeElementQuad4:
188 # Create the indices for the linear element
189 first_time_row_start_index = i_element_time * number_of_nodes_in_space
190 second_time_row_start_index = (
191 1 + i_element_time
192 ) * number_of_nodes_in_space
193 element_node_indices = [
194 first_time_row_start_index + element_node_ids[0],
195 first_time_row_start_index + element_node_ids[1],
196 second_time_row_start_index + element_node_ids[1],
197 second_time_row_start_index + element_node_ids[0],
198 ]
199 elif space_time_element_type == SpaceTimeElementQuad9:
200 # Create the indices for the quadratic element
201 first_time_row_start_index = (
202 2 * i_element_time * number_of_nodes_in_space
203 )
204 second_time_row_start_index = (
205 2 * i_element_time + 1
206 ) * number_of_nodes_in_space
207 third_time_row_start_index = (
208 2 * i_element_time + 2
209 ) * number_of_nodes_in_space
210 element_node_indices = [
211 first_time_row_start_index + element_node_ids[0],
212 first_time_row_start_index + element_node_ids[2],
213 third_time_row_start_index + element_node_ids[2],
214 third_time_row_start_index + element_node_ids[0],
215 first_time_row_start_index + element_node_ids[1],
216 second_time_row_start_index + element_node_ids[2],
217 third_time_row_start_index + element_node_ids[1],
218 second_time_row_start_index + element_node_ids[0],
219 second_time_row_start_index + element_node_ids[1],
220 ]
221 else:
222 raise TypeError(
223 f"Got unexpected space time element type {space_time_element_type}"
224 )
226 # Add the element to the mesh
227 space_time_elements.append(
228 space_time_element_type(
229 [space_time_nodes[i_node] for i_node in element_node_indices]
230 )
231 )
233 # Add joints to the space time mesh
234 space_time_couplings = []
235 for coupling in mesh_space_reference.boundary_conditions[
236 _mpy.bc.point_coupling, _mpy.geo.point
237 ]:
238 for i_mesh_space in range(number_of_copies_in_time):
239 coupling_node_ids = [
240 node.i_global for node in coupling.geometry_set.get_points()
241 ]
242 space_time_couplings.append(
243 _Coupling(
244 [
245 space_time_nodes[node_id + number_of_nodes_in_space]
246 for node_id in coupling_node_ids
247 ],
248 coupling.bc_type,
249 coupling.data,
250 )
251 )
253 # Create the new mesh and add all the mesh items
254 space_time_mesh = _Mesh()
255 space_time_mesh.add(space_time_nodes)
256 space_time_mesh.add(space_time_elements)
257 space_time_mesh.add(space_time_couplings)
259 # Create the element sets
260 return_set = _GeometryName()
261 return_set["start"] = _GeometrySet(start_nodes)
262 return_set["end"] = _GeometrySet(end_nodes)
263 return_set["left"] = _GeometrySet(left_nodes)
264 return_set["right"] = _GeometrySet(right_nodes)
265 return_set["surface"] = _GeometrySetNodes(_mpy.geo.surface, space_time_mesh.nodes)
267 return space_time_mesh, return_set
270def mesh_to_data_arrays(mesh: _Mesh):
271 """Get the relevant data arrays from the space time mesh."""
273 element_types = list(set([type(element) for element in mesh.elements]))
274 if len(element_types) > 1:
275 raise ValueError("Got more than a single element type, this is not supported")
276 elif not (
277 element_types[0] == SpaceTimeElementQuad4
278 or element_types[0] == SpaceTimeElementQuad9
279 ):
280 raise TypeError(
281 f"Expected either SpaceTimeElementQuad4 or SpaceTimeElementQuad9, got {element_types[0]}"
282 )
284 _, raw_unique_nodes = _get_coupled_nodes_to_master_map(mesh, assign_i_global=True)
285 unique_nodes = _cast(_List[NodeCosseratSpaceTime], raw_unique_nodes)
287 n_nodes = len(unique_nodes)
288 n_elements = len(mesh.elements)
289 n_nodes_per_element = len(mesh.elements[0].nodes)
291 coordinates = _get_nodal_coordinates(unique_nodes)
292 time = _np.zeros(n_nodes)
293 connectivity = _np.zeros((n_elements, n_nodes_per_element), dtype=int)
294 element_rotation_vectors = _np.zeros((n_elements, n_nodes_per_element, 3))
296 for i_node, node in enumerate(unique_nodes):
297 time[i_node] = node.time
299 for i_element, element in enumerate(mesh.elements):
300 for i_node, node in enumerate(element.nodes):
301 connectivity[i_element, i_node] = node.i_global
302 element_rotation_vectors[i_element, i_node, :] = (
303 node.rotation.get_rotation_vector()
304 )
306 geometry_sets = mesh.get_unique_geometry_sets()
307 node_sets: _Dict[str, _List] = {}
308 for value in geometry_sets.values():
309 for geometry_set in value:
310 node_sets[str(len(node_sets) + 1)] = _np.array(
311 [node.i_global for node in geometry_set.get_all_nodes()]
312 )
314 return_dict = {
315 "coordinates": coordinates,
316 "time": time,
317 "connectivity": connectivity,
318 "element_rotation_vectors": element_rotation_vectors,
319 "node_sets": node_sets,
320 }
322 nodes_have_arc_length = {node.arc_length is not None for node in mesh.nodes}
323 if len(nodes_have_arc_length) > 1:
324 raise ValueError(
325 "Some nodes have an arc length, some don't. This is not supported."
326 )
327 if nodes_have_arc_length.pop():
328 # The arc length is added as an "element" property, since the same
329 # node can have a different arc length depending on the element
330 # (similar to the rotation vectors).
331 arc_length = _np.zeros((n_elements, n_nodes_per_element))
332 for i_element, element in enumerate(mesh.elements):
333 for i_node, node in enumerate(element.nodes):
334 connectivity[i_element, i_node] = node.i_global
335 arc_length[i_element, i_node] = node.arc_length
336 return_dict["arc_length"] = arc_length
338 return return_dict