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