Coverage for src/beamme/cosserat_curve/warping_along_cosserat_curve.py: 97%
101 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"""This file contains functionality to warp an existing mesh along a 1D
23curve."""
25from typing import Tuple as _Tuple
27import numpy as _np
28import quaternion as _quaternion
29from numpy.typing import NDArray as _NDArray
31from beamme.core.boundary_condition import BoundaryCondition as _BoundaryCondition
32from beamme.core.conf import bme as _bme
33from beamme.core.geometry_set import GeometrySet as _GeometrySet
34from beamme.core.mesh import Mesh as _Mesh
35from beamme.core.node import Node as _Node
36from beamme.core.node import NodeCosserat as _NodeCosserat
37from beamme.core.rotation import Rotation as _Rotation
38from beamme.cosserat_curve.cosserat_curve import CosseratCurve as _CosseratCurve
39from beamme.four_c.function_utility import (
40 create_linear_interpolation_function as _create_linear_interpolation_function,
41)
42from beamme.geometric_search.find_close_points import (
43 find_close_points as _find_close_points,
44)
47def get_arc_length_and_cross_section_coordinates(
48 coordinates: _np.ndarray, origin: _np.ndarray, reference_rotation: _Rotation
49) -> _Tuple[float, _np.ndarray]:
50 """Return the arc length and the cross section coordinates for a coordinate
51 system defined by the reference rotation and the origin.
53 Args
54 ----
55 coordinates:
56 Point coordinates in R3
57 origin:
58 Origin of the coordinate system
59 reference_rotation:
60 Rotation of the coordinate system. The first basis vector is the arc
61 length direction.
62 """
64 transformed_coordinates = reference_rotation.inv() * (coordinates - origin)
65 centerline_position = transformed_coordinates[0]
66 cross_section_coordinates = [0.0, *transformed_coordinates[1:]]
67 return centerline_position, cross_section_coordinates
70def get_mesh_transformation(
71 curve: _CosseratCurve,
72 nodes: list[_Node],
73 *,
74 origin=[0.0, 0.0, 0.0],
75 reference_rotation=_Rotation(),
76 n_steps: int = 10,
77 initial_configuration: bool = True,
78 **kwargs,
79) -> _Tuple[_np.ndarray, _NDArray[_quaternion.quaternion]]:
80 """Generate a list of positions for each node that describe the
81 transformation of the nodes from the given configuration to the Cosserat
82 curve.
84 Args
85 ----
86 curve:
87 Curve to warp the mesh to
88 nodes:
89 Optional, if this is given only warp the given nodes. Per default all nodes in
90 the mesh are warped.
91 origin:
92 Origin of the coordinate system
93 reference_rotation:
94 Rotation of the coordinate system. The first basis vector is the arc
95 length direction.
96 n_steps:
97 Number of steps to get from the unwrapped configuration to the final configuration.
98 initial_configuration:
99 If the initial, unwrapped configuration (factor=0) should also be added to the
100 results.
101 kwargs:
102 Keyword arguments passed to CosseratCurve.get_centerline_positions_and_rotations
104 Return
105 ----
106 positions: list(_np.array(n_nodes x 3))
107 A list for each time step containing the position of all nodes for that time step
108 relative_rotations: list(list(Rotation))
109 A list for each time step containing the relative rotations for all nodes at that
110 time step
111 """
113 # Define the factors for which we will generate the positions and rotations
114 factors = _np.linspace(0.0, 1.0, n_steps + 1)
115 if initial_configuration:
116 n_output_steps = n_steps + 1
117 else:
118 n_output_steps = n_steps
119 factors = _np.delete(factors, 0)
121 # Create output arrays
122 n_nodes = len(nodes)
123 positions = _np.zeros((n_output_steps, n_nodes, 3))
124 relative_rotations = _np.zeros(
125 (n_output_steps, n_nodes), dtype=_quaternion.quaternion
126 )
128 # Get all arc lengths and cross section positions
129 arc_lengths = _np.zeros((n_nodes, 1))
130 cross_section_coordinates = [None] * n_nodes
131 for i_node, node in enumerate(nodes):
132 (
133 arc_lengths[i_node],
134 cross_section_coordinates[i_node],
135 ) = get_arc_length_and_cross_section_coordinates(
136 node.coordinates, origin, reference_rotation
137 )
139 # Get unique arc length points
140 has_partner, n_partner = _find_close_points(arc_lengths)
141 arc_lengths_unique = [None] * n_partner
142 has_partner_total = [-2] * len(arc_lengths)
143 for i in range(len(arc_lengths)):
144 partner_id = has_partner[i]
145 if partner_id == -1:
146 has_partner_total[i] = len(arc_lengths_unique)
147 arc_lengths_unique.append(arc_lengths[i][0])
148 else:
149 if arc_lengths_unique[partner_id] is None:
150 arc_lengths_unique[partner_id] = arc_lengths[i][0]
151 has_partner_total[i] = partner_id
153 n_total = len(arc_lengths_unique)
154 arc_lengths_unique = _np.array(arc_lengths_unique)
155 arc_lengths_sorted_index = _np.argsort(arc_lengths_unique)
156 arc_lengths_sorted = arc_lengths_unique[arc_lengths_sorted_index]
157 arc_lengths_sorted_index_inv = [-2 for i in range(n_total)]
158 for i in range(n_total):
159 arc_lengths_sorted_index_inv[arc_lengths_sorted_index[i]] = i
160 point_to_unique = []
161 for partner in has_partner_total:
162 point_to_unique.append(arc_lengths_sorted_index_inv[partner])
164 # Get all configurations for the unique points
165 positions_for_all_steps = []
166 quaternions_for_all_steps = []
168 for factor in factors:
169 sol_r, sol_q = curve.get_centerline_positions_and_rotations(
170 arc_lengths_sorted, factor=factor, **kwargs
171 )
172 positions_for_all_steps.append(sol_r)
173 quaternions_for_all_steps.append(sol_q)
175 # Get data required for the rigid body motion
176 curve_start_pos, curve_start_rot = curve.get_centerline_position_and_rotation(0.0)
177 rigid_body_translation = curve_start_pos - origin
178 rigid_body_rotation = curve_start_rot
180 # Loop over nodes and map them to the new configuration
181 for i_node, node in enumerate(nodes):
182 if not isinstance(node, _Node):
183 raise TypeError(
184 "All nodes in the mesh have to be derived from the base Node object"
185 )
187 node_unique_id = point_to_unique[i_node]
188 cross_section_position = cross_section_coordinates[i_node]
190 # Check that the arc length coordinates match
191 if (
192 _np.abs(arc_lengths[i_node] - arc_lengths_sorted[node_unique_id])
193 > _bme.eps_pos
194 ):
195 raise ValueError("Arc lengths do not match")
197 # Create the functions that describe the deformation
198 for i_step, factor in enumerate(factors):
199 centerline_pos = positions_for_all_steps[i_step][node_unique_id]
200 centerline_relative_pos = _quaternion.rotate_vectors(
201 curve_start_rot.conjugate(), centerline_pos - curve_start_pos
202 )
203 centerline_rotation = quaternions_for_all_steps[i_step][node_unique_id]
204 centerline_relative_rotation = (
205 curve_start_rot.conjugate() * centerline_rotation
206 )
208 rigid_body_rotation_for_factor = _quaternion.slerp_evaluate(
209 reference_rotation.get_numpy_quaternion(), rigid_body_rotation, factor
210 )
212 current_pos = (
213 _quaternion.rotate_vectors(
214 rigid_body_rotation_for_factor,
215 (
216 centerline_relative_pos
217 + _quaternion.rotate_vectors(
218 centerline_relative_rotation, cross_section_position
219 )
220 ),
221 )
222 + origin
223 + factor * rigid_body_translation
224 )
226 positions[i_step, i_node] = current_pos
227 relative_rotations[i_step, i_node] = (
228 rigid_body_rotation_for_factor
229 * centerline_relative_rotation
230 * reference_rotation.get_numpy_quaternion().conjugate()
231 )
233 return positions, relative_rotations
236def create_transform_boundary_conditions(
237 mesh: _Mesh,
238 curve: _CosseratCurve,
239 *,
240 nodes: list[_Node] | None = None,
241 t_end: float = 1.0,
242 n_steps: int = 10,
243 n_dof_per_node: int = 3,
244 **kwargs,
245) -> None:
246 """Create the Dirichlet boundary conditions that enforce the warping. The
247 warped object is assumed to align with the x-axis in the reference
248 configuration.
250 Args
251 ----
252 mesh:
253 Mesh to be warped
254 curve:
255 Curve to warp the mesh to
256 nodes:
257 Optional, if this is given only warp the given nodes. Per default all nodes in
258 the mesh are warped.
259 n_steps:
260 Number of steps to apply the warping condition
261 t_end:
262 End time for applying the warping boundary conditions
263 n_dof_per_node:
264 Number of DOF per node in 4C (is needed to correctly define the boundary conditions)
265 kwargs:
266 Keyword arguments passed to get_mesh_transformation
267 """
269 # If no nodes are given, use all nodes in the mesh
270 if nodes is None:
271 nodes = mesh.nodes
273 time_values = _np.linspace(0.0, t_end, n_steps + 1)
275 # Get positions and rotations for each step
276 positions, _ = get_mesh_transformation(curve, nodes, n_steps=n_steps, **kwargs)
278 # Loop over nodes and map them to the new configuration
279 for i_node, node in enumerate(nodes):
280 # Create the functions that describe the deformation
281 reference_position = node.coordinates
282 displacement_values = _np.array(
283 [
284 positions[i_step][i_node] - reference_position
285 for i_step in range(n_steps + 1)
286 ]
287 )
288 fun_pos = [
289 _create_linear_interpolation_function(
290 time_values, displacement_values[:, i_dir]
291 )
292 for i_dir in range(3)
293 ]
294 for fun in fun_pos:
295 mesh.add(fun)
296 n_additional_dof = n_dof_per_node - 3
297 mesh.add(
298 _BoundaryCondition(
299 _GeometrySet(node),
300 {
301 "NUMDOF": n_dof_per_node,
302 "ONOFF": [1] * 3 + [0] * n_additional_dof,
303 "VAL": [1.0] * 3 + [0.0] * n_additional_dof,
304 "FUNCT": fun_pos + [None] * n_additional_dof,
305 "TAG": "monitor_reaction",
306 },
307 bc_type=_bme.bc.dirichlet,
308 )
309 )
312def warp_mesh_along_curve(
313 mesh: _Mesh,
314 curve: _CosseratCurve,
315 *,
316 origin=[0.0, 0.0, 0.0],
317 reference_rotation=_Rotation(),
318) -> None:
319 """Warp an existing mesh along the given curve.
321 The reference coordinates for the transformation are defined by the
322 given origin and rotation, where the first basis vector of the triad
323 defines the centerline axis.
324 """
326 pos, rot = get_mesh_transformation(
327 curve,
328 mesh.nodes,
329 origin=origin,
330 reference_rotation=reference_rotation,
331 n_steps=1,
332 initial_configuration=False,
333 )
335 # Loop over nodes and map them to the new configuration
336 for i_node, node in enumerate(mesh.nodes):
337 if not isinstance(node, _Node):
338 raise TypeError(
339 "All nodes in the mesh have to be derived from the base Node object"
340 )
342 new_pos = pos[0, i_node]
343 node.coordinates = new_pos
344 if isinstance(node, _NodeCosserat):
345 node.rotation = _Rotation.from_quaternion(rot[0, i_node]) * node.rotation