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