Coverage for src/meshpy/space_time/beam_to_space_time.py: 94%
142 statements
« prev ^ index » next coverage.py v7.8.0, created at 2025-04-28 04:21 +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 # For some space-time formulations it is required that the
167 # pre-processing provides the arc-length along a beam filament.
168 # Since the basic space meshes here contain nodes that don't have the
169 # arc_length attribute by default, we have to do the following check
170 # and branching.
171 # TODO: Think if it makes sense to somehow add this as a member to the
172 # default Node object.
173 nodes_have_arc_length_attribute = {
174 hasattr(node, "arc_length") for node in mesh_space_current_time.nodes
175 }
176 if not len(nodes_have_arc_length_attribute) == 1:
177 raise ValueError(
178 "There are some nodes in the mesh with the arc_length attribute and "
179 "some without it. This is not supported."
180 )
181 if nodes_have_arc_length_attribute.pop():
182 space_time_nodes_to_add = [
183 NodeCosseratSpaceTime(
184 node.coordinates, node.rotation, time, arc_length=node.arc_length
185 )
186 for node in mesh_space_current_time.nodes
187 ]
188 else:
189 space_time_nodes_to_add = [
190 NodeCosseratSpaceTime(node.coordinates, node.rotation, time)
191 for node in mesh_space_current_time.nodes
192 ]
193 space_time_nodes.extend(space_time_nodes_to_add)
195 if i_mesh_space == 0:
196 start_nodes = space_time_nodes_to_add
197 elif i_mesh_space == number_of_copies_in_time - 1:
198 end_nodes = space_time_nodes_to_add
200 left_nodes.append(space_time_nodes_to_add[0])
201 right_nodes.append(space_time_nodes_to_add[-1])
203 # Create the space time elements
204 space_time_elements = []
205 for i_element_time in range(number_of_elements_in_time):
206 for element in mesh_space_reference.elements:
207 element_node_ids = [node.i_global for node in element.nodes]
208 if space_time_element_type == SpaceTimeElementQuad4:
209 # Create the indices for the linear element
210 first_time_row_start_index = i_element_time * number_of_nodes_in_space
211 second_time_row_start_index = (
212 1 + i_element_time
213 ) * number_of_nodes_in_space
214 element_node_indices = [
215 first_time_row_start_index + element_node_ids[0],
216 first_time_row_start_index + element_node_ids[1],
217 second_time_row_start_index + element_node_ids[1],
218 second_time_row_start_index + element_node_ids[0],
219 ]
220 elif space_time_element_type == SpaceTimeElementQuad9:
221 # Create the indices for the quadratic element
222 first_time_row_start_index = (
223 2 * i_element_time * number_of_nodes_in_space
224 )
225 second_time_row_start_index = (
226 2 * i_element_time + 1
227 ) * number_of_nodes_in_space
228 third_time_row_start_index = (
229 2 * i_element_time + 2
230 ) * number_of_nodes_in_space
231 element_node_indices = [
232 first_time_row_start_index + element_node_ids[0],
233 first_time_row_start_index + element_node_ids[2],
234 third_time_row_start_index + element_node_ids[2],
235 third_time_row_start_index + element_node_ids[0],
236 first_time_row_start_index + element_node_ids[1],
237 second_time_row_start_index + element_node_ids[2],
238 third_time_row_start_index + element_node_ids[1],
239 second_time_row_start_index + element_node_ids[0],
240 second_time_row_start_index + element_node_ids[1],
241 ]
242 else:
243 raise TypeError(
244 f"Got unexpected space time element type {space_time_element_type}"
245 )
247 # Add the element to the mesh
248 space_time_elements.append(
249 space_time_element_type(
250 [space_time_nodes[i_node] for i_node in element_node_indices]
251 )
252 )
254 # Add joints to the space time mesh
255 space_time_couplings = []
256 for coupling in mesh_space_reference.boundary_conditions[
257 _mpy.bc.point_coupling, _mpy.geo.point
258 ]:
259 for i_mesh_space in range(number_of_copies_in_time):
260 coupling_node_ids = [
261 node.i_global for node in coupling.geometry_set.get_points()
262 ]
263 space_time_couplings.append(
264 _Coupling(
265 [
266 space_time_nodes[node_id + number_of_nodes_in_space]
267 for node_id in coupling_node_ids
268 ],
269 coupling.bc_type,
270 coupling.coupling_dof_type,
271 )
272 )
274 # Create the new mesh and add all the mesh items
275 space_time_mesh = _Mesh()
276 space_time_mesh.add(space_time_nodes)
277 space_time_mesh.add(space_time_elements)
278 space_time_mesh.add(space_time_couplings)
280 # Create the element sets
281 return_set = _GeometryName()
282 return_set["start"] = _GeometrySet(start_nodes)
283 return_set["end"] = _GeometrySet(end_nodes)
284 return_set["left"] = _GeometrySet(left_nodes)
285 return_set["right"] = _GeometrySet(right_nodes)
286 return_set["surface"] = _GeometrySetNodes(_mpy.geo.surface, space_time_mesh.nodes)
288 return space_time_mesh, return_set
291def mesh_to_data_arrays(mesh: _Mesh):
292 """Get the relevant data arrays from the space time mesh."""
294 element_types = list(set([type(element) for element in mesh.elements]))
295 if len(element_types) > 1:
296 raise ValueError("Got more than a single element type, this is not supported")
297 elif not (
298 element_types[0] == SpaceTimeElementQuad4
299 or element_types[0] == SpaceTimeElementQuad9
300 ):
301 raise TypeError(
302 f"Expected either SpaceTimeElementQuad4 or SpaceTimeElementQuad9, got {element_types[0]}"
303 )
305 _, raw_unique_nodes = _get_coupled_nodes_to_master_map(mesh, assign_i_global=True)
306 unique_nodes = _cast(_List[NodeCosseratSpaceTime], raw_unique_nodes)
308 n_nodes = len(unique_nodes)
309 n_elements = len(mesh.elements)
310 n_nodes_per_element = len(mesh.elements[0].nodes)
312 coordinates = _get_nodal_coordinates(unique_nodes)
313 time = _np.zeros(n_nodes)
314 connectivity = _np.zeros((n_elements, n_nodes_per_element), dtype=int)
315 element_rotation_vectors = _np.zeros((n_elements, n_nodes_per_element, 3))
317 for i_node, node in enumerate(unique_nodes):
318 time[i_node] = node.time
320 for i_element, element in enumerate(mesh.elements):
321 for i_node, node in enumerate(element.nodes):
322 connectivity[i_element, i_node] = node.i_global
323 element_rotation_vectors[i_element, i_node, :] = (
324 node.rotation.get_rotation_vector()
325 )
327 geometry_sets = mesh.get_unique_geometry_sets()
328 node_sets: _Dict[str, _List] = {}
329 for value in geometry_sets.values():
330 for geometry_set in value:
331 node_sets[str(len(node_sets) + 1)] = _np.array(
332 [node.i_global for node in geometry_set.get_all_nodes()]
333 )
335 return_dict = {
336 "coordinates": coordinates,
337 "time": time,
338 "connectivity": connectivity,
339 "element_rotation_vectors": element_rotation_vectors,
340 "node_sets": node_sets,
341 }
343 # We assume that either all nodes have an arc_length or no one.
344 # This is checked in the beam_to_space_time function.
345 if mesh.nodes[0].arc_length is not None:
346 # The arc length is added as an "element" property, since the same
347 # node can have a different arc length depending on the element
348 # (similar to the rotation vectors).
349 arc_length = _np.zeros((n_elements, n_nodes_per_element))
350 for i_element, element in enumerate(mesh.elements):
351 for i_node, node in enumerate(element.nodes):
352 connectivity[i_element, i_node] = node.i_global
353 arc_length[i_element, i_node] = node.arc_length
354 return_dict["arc_length"] = arc_length
356 return return_dict