Coverage for src/meshpy/mesh_creation_functions/beam_curve.py: 96%
81 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"""This file has functions to create a beam from a parametric curve."""
24import autograd.numpy as _npAD
25import numpy as _np
26import scipy.integrate as _integrate
27import scipy.optimize as _optimize
28from autograd import jacobian as _jacobian
30from meshpy.core.conf import mpy as _mpy
31from meshpy.core.rotation import Rotation as _Rotation
32from meshpy.core.rotation import smallest_rotation as _smallest_rotation
33from meshpy.mesh_creation_functions.beam_generic import (
34 create_beam_mesh_function as _create_beam_mesh_function,
35)
38def create_beam_mesh_curve(
39 mesh,
40 beam_class,
41 material,
42 function,
43 interval,
44 *,
45 output_length=False,
46 function_derivative=None,
47 function_rotation=None,
48 **kwargs,
49):
50 """Generate a beam from a parametric curve. Integration along the beam is
51 performed with scipy, and if the gradient is not explicitly provided, it is
52 calculated with the numpy wrapper autograd.
54 Args
55 ----
56 mesh: Mesh
57 Mesh that the curve will be added to.
58 beam_class: Beam
59 Class of beam that will be used for this line.
60 material: Material
61 Material for this line.
62 function: function
63 3D-parametric curve that represents the beam axis. If only a 2D
64 point is returned, the triad creation is simplified. If
65 mathematical functions are used, they have to come from the wrapper
66 autograd.numpy.
67 interval: [start end]
68 Start and end values for the parameter of the curve.
69 output_length: bool
70 If this is true, the function returns a tuple containing the created
71 sets and the total arc length along the integrated function.
72 function_derivative: function -> R3
73 Explicitly provide the jacobian of the centerline position.
74 function_rotation: function -> Rotation
75 If this argument is given, the triads are computed with this
76 function, on the same interval as the position function. Must
77 return a Rotation object.
78 If no function_rotation is given, the rotation of the first node
79 is calculated automatically and all subsequent nodal rotations
80 are calculated based on a smallest rotation mapping onto the curve
81 tangent vector.
83 **kwargs (for all of them look into create_beam_mesh_function)
84 ----
85 n_el: int
86 Number of equally spaced beam elements along the line. Defaults to 1.
87 Mutually exclusive with l_el.
88 l_el: float
89 Desired length of beam elements. Mutually exclusive with n_el.
90 Be aware, that this length might not be achieved, if the elements are
91 warped after they are created.
93 Return
94 ----
95 return_set: GeometryName
96 Set with the 'start' and 'end' node of the curve. Also a 'line' set
97 with all nodes of the curve.
98 """
100 # Check size of position function
101 if len(function(interval[0])) == 2:
102 is_3d_curve = False
103 elif len(function(interval[0])) == 3:
104 is_3d_curve = True
105 else:
106 raise ValueError("Function must return either 2d or 3d curve!")
108 # Check rotation function.
109 if function_rotation is None:
110 is_rot_funct = False
111 else:
112 is_rot_funct = True
114 # Check that the position is an np.array
115 if not isinstance(function(interval[0]), _np.ndarray):
116 raise TypeError(
117 "Function must be of type np.ndarray, got {}!".format(
118 type(function(interval[0]))
119 )
120 )
122 # Get the derivative of the position function and the increment along
123 # the curve.
124 if function_derivative is None:
125 rp = _jacobian(function)
126 else:
127 rp = function_derivative
129 # Check which one of the boundaries is larger.
130 if interval[0] > interval[1]:
131 # In this case rp needs to be negated.
132 rp_positive = rp
134 def rp(t):
135 """Return the inverted tangent vector."""
136 return -(rp_positive(t))
138 def ds(t):
139 """Increment along the curve."""
140 return _npAD.linalg.norm(rp(t))
142 def S(t, start_t=None, start_S=None):
143 """Function that integrates the length until a parameter value.
145 A speedup can be achieved by giving start_t and start_S, the
146 parameter and Length at a known point.
147 """
148 if start_t is None and start_S is None:
149 st = interval[0]
150 sS = 0
151 elif start_t is not None and start_S is not None:
152 st = start_t
153 sS = start_S
154 else:
155 raise ValueError("Input parameters are wrong!")
156 return _integrate.quad(ds, st, t)[0] + sS
158 def get_t_along_curve(arc_length, t0, **kwargs):
159 """Calculate the parameter t where the length along the curve is
160 arc_length.
162 t0 is the start point for the Newton iteration.
163 """
164 t_root = _optimize.newton(lambda t: S(t, **kwargs) - arc_length, t0, fprime=ds)
165 return t_root
167 class BeamFunctions:
168 """This object manages the creation of functions which are used to
169 create the beam nodes.
171 By wrapping this in this object, we can store some data and
172 speed up the numerical integration along the curve.
173 """
175 def __init__(self):
176 """Initialize the object."""
177 self._reset_start_values()
178 self.last_triad = None
180 if is_3d_curve:
181 r_prime = rp(interval[0])
182 if abs(_np.dot(r_prime, [0, 0, 1])) < abs(_np.dot(r_prime, [0, 1, 0])):
183 t2_temp = [0, 0, 1]
184 else:
185 t2_temp = [0, 1, 0]
186 self.last_triad = _Rotation.from_basis(r_prime, t2_temp)
188 def _reset_start_values(self):
189 """Reset the stored start values for the next Newton iteration."""
190 self.t_start_newton = interval[0]
191 self.S_start_newton = 0.0
193 def __call__(self, length_a, length_b):
194 """This object is called with the interval limits.
196 This method returns a function that evaluates the position
197 and rotation within this interval.
198 """
200 # In case the interval is not continuous with the last one, we reset the start
201 # values for the Newton iteration here
202 if length_a < self.S_start_newton - _mpy.eps_pos:
203 self._reset_start_values()
205 # Length of the beam element in physical space.
206 L = length_b - length_a
208 def get_beam_position_and_rotation_at_xi(xi):
209 """Evaluate the beam position and rotation at xi.
211 xi is the beam element parameter coordinate, i.e., xi =
212 [-1, 1].
213 """
214 # Parameter for xi.
215 S = length_a + 0.5 * (xi + 1) * L
216 t = get_t_along_curve(
217 S,
218 self.t_start_newton,
219 start_t=self.t_start_newton,
220 start_S=self.S_start_newton,
221 )
223 # Position at xi.
224 if is_3d_curve:
225 pos = function(t)
226 else:
227 pos = _np.zeros(3)
228 pos[:2] = function(t)
230 # Rotation at xi.
231 if is_rot_funct:
232 rot = function_rotation(t)
233 else:
234 r_prime = rp(t)
235 if is_3d_curve:
236 # Create the next triad via the smallest rotation mapping based
237 # on the last triad.
238 rot = _smallest_rotation(self.last_triad, r_prime)
239 self.last_triad = rot.copy()
240 else:
241 # The rotation simplifies in the 2d case.
242 rot = _Rotation([0, 0, 1], _np.arctan2(r_prime[1], r_prime[0]))
244 # Set start values for the next iteration
245 self.t_start_newton = t
246 self.S_start_newton = S
248 # Return the needed values for beam creation.
249 return (pos, rot, S)
251 return get_beam_position_and_rotation_at_xi
253 # Now create the beam.
254 # Get the length of the whole segment.
255 length = S(interval[1])
257 # Create the beam in the mesh
258 created_sets = _create_beam_mesh_function(
259 mesh,
260 beam_class=beam_class,
261 material=material,
262 function_generator=BeamFunctions(),
263 interval=[0.0, length],
264 interval_length=length,
265 **kwargs,
266 )
268 if output_length:
269 return (created_sets, length)
270 else:
271 return created_sets