Coverage for src/meshpy/cosserat_curve/warping_along_cosserat

1# The MIT License (MIT)

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

104

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 """

113

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)

121

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 )

128

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 )

139

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

153

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])

164

165 # Get all configurations for the unique points

166 positions_for_all_steps = []

167 quaternions_for_all_steps = []

168

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)

175

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

180

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 )

187

188 node_unique_id = point_to_unique[i_node]

189 cross_section_position = cross_section_coordinates[i_node]

190

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")

197

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 )

208

209 rigid_body_rotation_for_factor = _quaternion.slerp_evaluate(

210 _quaternion.from_float_array(reference_rotation.q),

211 rigid_body_rotation,

212 factor,

213 )

214

215 current_pos = (

216 _quaternion.rotate_vectors(

217 rigid_body_rotation_for_factor,

218 (

219 centerline_relative_pos

220 + _quaternion.rotate_vectors(

221 centerline_relative_rotation, cross_section_position

222 )

223 ),

224 )

225 + origin

226 + factor * rigid_body_translation

227 )

228

229 positions[i_step, i_node] = current_pos

230 relative_rotations[i_step, i_node] = (

231 rigid_body_rotation_for_factor

232 * centerline_relative_rotation

233 * _quaternion.from_float_array(reference_rotation.q).conjugate()

234 )

235

236 return positions, relative_rotations

237

238

239def create_transform_boundary_conditions(

240 mesh: _Mesh,

241 curve: _CosseratCurve,

242 *,

243 nodes: _Optional[list[_Node]] = None,

244 t_end: float = 1.0,

245 n_steps: int = 10,

246 n_dof_per_node: int = 3,

247 **kwargs,

248) -> None:

249 """Create the Dirichlet boundary conditions that enforce the warping. The

250 warped object is assumed to align with the z-axis in the reference

251 configuration.

252

253 Args

254 ----

255 mesh:

256 Mesh to be warped

257 curve:

258 Curve to warp the mesh to

259 nodes:

260 Optional, if this is given only warp the given nodes. Per default all nodes in

261 the mesh are warped.

262 n_steps:

263 Number of steps to apply the warping condition

264 t_end:

265 End time for applying the warping boundary conditions

266 n_dof_per_node:

267 Number of DOF per node in 4C (is needed to correctly define the boundary conditions)

268 kwargs:

269 Keyword arguments passed to get_mesh_transformation

270 """

271

272 # If no nodes are given, use all nodes in the mesh

273 if nodes is None:

274 nodes = mesh.nodes

275

276 time_values = _np.linspace(0.0, t_end, n_steps + 1)

277

278 # Get positions and rotations for each step

279 positions, _ = get_mesh_transformation(curve, nodes, n_steps=n_steps, **kwargs)

280

281 # Loop over nodes and map them to the new configuration

282 for i_node, node in enumerate(nodes):

283 # Create the functions that describe the deformation

284 reference_position = node.coordinates

285 displacement_values = _np.array(

286 [

287 positions[i_step][i_node] - reference_position

288 for i_step in range(n_steps + 1)

289 ]

290 )

291 fun_pos = [

292 _create_linear_interpolation_function(

293 time_values, displacement_values[:, i_dir]

294 )

295 for i_dir in range(3)

296 ]

297 for fun in fun_pos:

298 mesh.add(fun)

299 n_additional_dof = n_dof_per_node - 3

300 mesh.add(

301 _BoundaryCondition(

302 _GeometrySet(node),

303 {

304 "NUMDOF": n_dof_per_node,

305 "ONOFF": [1] * 3 + [0] * n_additional_dof,

306 "VAL": [1.0] * 3 + [0.0] * n_additional_dof,

307 "FUNCT": fun_pos + [None] * n_additional_dof,

308 "TAG": "monitor_reaction",

309 },

310 bc_type=_mpy.bc.dirichlet,

311 )

312 )

313

314

315def warp_mesh_along_curve(

316 mesh: _Mesh,

317 curve: _CosseratCurve,

318 *,

319 origin=[0.0, 0.0, 0.0],

320 reference_rotation=_Rotation(),

321) -> None:

322 """Warp an existing mesh along the given curve.

323

324 The reference coordinates for the transformation are defined by the

325 given origin and rotation, where the first basis vector of the triad

326 defines the centerline axis.

327 """

328

329 pos, rot = get_mesh_transformation(

330 curve,

331 mesh.nodes,

332 origin=origin,

333 reference_rotation=reference_rotation,

334 n_steps=1,

335 initial_configuration=False,

336 )

337

338 # Loop over nodes and map them to the new configuration

339 for i_node, node in enumerate(mesh.nodes):

340 if not isinstance(node, _Node):

341 raise TypeError(

342 "All nodes in the mesh have to be derived from the base Node object"

343 )

344

345 new_pos = pos[0, i_node]

346 node.coordinates = new_pos

347 if isinstance(node, _NodeCosserat):

348 node.rotation = (

349 _Rotation.from_quaternion(_quaternion.as_float_array(rot[0, i_node]))

350 * node.rotation

351 )

Coverage for src/meshpy/cosserat_curve/warping_along_cosserat_curve.py: 97%

102 statements