Coverage for src/beamme/mesh_creation_functions/beam_parametric_curve.py: 98%
160 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 has functions to create a beam from a parametric curve."""
24from collections.abc import Callable as _Callable
26import numpy as _np
27import scipy.integrate as _integrate
28from autograd import jacobian as _jacobian
29from scipy.integrate import quad as _quad
30from scipy.interpolate import interp1d as _interp1d
32from beamme.core.element_beam import Beam as _Beam
33from beamme.core.geometry_set import GeometryName as _GeometryName
34from beamme.core.material import MaterialBeamBase as _MaterialBeamBase
35from beamme.core.mesh import Mesh as _Mesh
36from beamme.core.rotation import Rotation as _Rotation
37from beamme.core.rotation import smallest_rotation as _smallest_rotation
38from beamme.mesh_creation_functions.beam_generic import (
39 create_beam_mesh_generic as _create_beam_mesh_generic,
40)
43class _ArcLengthEvaluation:
44 """Class to allow evaluation of the arc length S(t) and the inverse mapping t(S).
46 This class uses precomputed samples to interpolate between arc length values. This
47 is much more efficient than root finding algorithms and should provide a suitable
48 accuracy.
49 """
51 def __init__(
52 self,
53 function_derivative: _Callable,
54 interval: tuple[float, float],
55 n_precomputed_intervals: int = 100,
56 scipy_integrate: bool = True,
57 scipy_integrate_points: list[float] | None = None,
58 method: str = "arc-length",
59 ) -> None:
60 """Initialize the arc length evaluator and precomputed samples.
62 Args:
63 function_derivative:
64 Function that returns the tangent vector along the curve for a given
65 parameter value t.
66 interval:
67 Start and end values for the parameter coordinate of the curve,
68 must be in ascending order.
69 n_precomputed_intervals:
70 Number of intervals to use for the pre-computation. More intervals
71 lead to higher accuracy.
72 scipy_integrate:
73 If true, the one single arc length integration is performed with scipy's
74 quad function. This is an adaptive method and gives good estimates on
75 subdivisions of the total interval.
76 scipy_integrate_points:
77 If scipy_integrate is true, this can be used to provide additional points
78 for the adaptive integration. This can be useful if there are known points
79 along the curve where Jacobian has a kink.
80 method:
81 Method to use for evaluation of nodal positions along the curve:
82 - "arc-length": Uniform spacing along the arc-length of the curve. This means
83 that elements along a curve will have equal length in space.
84 - "parametric": Uniform spacing along the parameter coordinate of the curve.
85 This means that elements along a curve will have equal length in parameter
86 space, but not necessarily in physical space.
87 - "parametric_consistent_middle_nodes": Same as `parametric`, but element
88 middle nodes are adjusted to be consistent with the arc-length mapping.
89 This means that elements along a curve will have equal length in parameter
90 space, but the middle nodes are placed such that the elements themselves are
91 not distorted.
92 """
93 if not scipy_integrate and scipy_integrate_points is not None:
94 raise ValueError(
95 "scipy_integrate_points cannot be provided if scipy_integrate is False!"
96 )
98 self.function_derivative = function_derivative
99 self.interval = interval
100 self.n_precomputed_intervals = n_precomputed_intervals
101 self.scipy_integrate = scipy_integrate
102 self.scipy_integrate_points = scipy_integrate_points
103 self.method = method
105 self._compute_samples()
106 self._compute_interpolation_functions()
108 def _compute_samples(self) -> None:
109 """Compute the samples for the arc length mapping.
111 This function computes the arc length S(t) at a set of sample points along the
112 parameter coordinate t with the accumulative Simpson integration.
113 """
114 if self.scipy_integrate:
115 ds_dt = lambda t: _np.linalg.norm(self.function_derivative([t])[0])
116 integral = _quad(
117 ds_dt,
118 self.interval[0],
119 self.interval[1],
120 points=self.scipy_integrate_points,
121 full_output=True,
122 )
123 integral_data = integral[2]
124 n_segments = integral_data["last"]
126 # Sort the segments such that they are in ascending order along the
127 # integration integral.
128 interval_sort = _np.argsort(integral_data["alist"][:n_segments])
129 intervals_integral = integral_data["rlist"][interval_sort]
130 intervals_left_end = integral_data["alist"][interval_sort]
131 intervals_right_end = integral_data["blist"][interval_sort]
133 self.t_grid = _np.zeros(n_segments * self.n_precomputed_intervals + 1)
134 self.S_grid = _np.zeros(n_segments * self.n_precomputed_intervals + 1)
135 local_integral_start = 0.0
136 for i_segment in range(n_segments):
137 segment_a = intervals_left_end[i_segment]
138 segment_b = intervals_right_end[i_segment]
140 global_start_index = i_segment * self.n_precomputed_intervals
141 global_end_index = (i_segment + 1) * self.n_precomputed_intervals + 1
143 t_local = _np.linspace(
144 segment_a, segment_b, self.n_precomputed_intervals + 1
145 )
146 self.t_grid[global_start_index:global_end_index] = t_local
148 local_integral = intervals_integral[i_segment]
150 tangents = self.function_derivative(t_local)
151 ds_at_t_samples = _np.linalg.norm(tangents, axis=1)
152 S_local = _integrate.cumulative_simpson(
153 y=ds_at_t_samples, x=t_local, initial=0.0
154 )
155 S_local *= local_integral / S_local[-1]
156 S_local += local_integral_start
158 self.S_grid[global_start_index:global_end_index] = S_local
160 local_integral_start += local_integral
162 else:
163 # Uniform grid in t for sampling.
164 self.t_grid = _np.linspace(
165 self.interval[0], self.interval[1], self.n_precomputed_intervals + 1
166 )
168 # Evaluate ds at all grid points.
169 tangents = self.function_derivative(self.t_grid)
170 ds_at_t_samples = _np.linalg.norm(tangents, axis=1)
172 self.S_grid = _integrate.cumulative_simpson(
173 y=ds_at_t_samples, x=self.t_grid, initial=0.0
174 )
176 def _compute_interpolation_functions(self) -> None:
177 """Setup the interpolation functions for S(t) and t(S)."""
178 self.S_from_t = _interp1d(
179 self.t_grid,
180 self.S_grid,
181 kind="cubic",
182 fill_value="extrapolate",
183 assume_sorted=True,
184 )
185 self.t_from_S = _interp1d(
186 self.S_grid,
187 self.t_grid,
188 kind="cubic",
189 fill_value="extrapolate",
190 assume_sorted=True,
191 )
193 def approximate_total_arc_length(self) -> float:
194 """Approximate the total arc length along the curve.
196 This value is only needed to choose the number of elements along the curve.
197 """
198 return self.S_grid[-1]
200 def get_total_arc_length(self) -> float:
201 """Get the total arc length along the curve.
203 This function might return a different arc-length than
204 `approximate_total_arc_length`, if the integral is adaptively refined in
205 `evaluate_all`.
206 """
207 return self.S_grid[-1]
209 def evaluate_all(
210 self, evaluation_points: _np.ndarray, middle_node_flags: _np.ndarray
211 ) -> tuple[_np.ndarray, _np.ndarray]:
212 """Evaluate the parameter coordinates corresponding to each node.
214 Args:
215 evaluation_points:
216 Evaluation points in the interval [0, 1]. Depending on the integration method,
217 this can be in normalized parameter space or arc-length space.
218 middle_node_flags:
219 Boolean array indicating which of the evaluation points
220 are middle nodes (True) and which are nodal points (False).
222 Returns:
223 t_evaluate:
224 Parameter coordinates along the curve for each evaluation point.
225 S_evaluate:
226 Arc-length coordinates along the curve for each evaluation point.
227 """
228 # Todo: Check if it makes sense to adaptively refine the arc-length
229 # integration here.
231 if self.method == "arc-length":
232 S_evaluate = evaluation_points * self.get_total_arc_length()
233 t_evaluate = self.t_from_S(S_evaluate)
235 elif self.method == "parametric":
236 interval_length = self.interval[1] - self.interval[0]
237 t_evaluate = self.interval[0] + evaluation_points * interval_length
238 S_evaluate = self.S_from_t(t_evaluate)
240 elif self.method == "parametric_consistent_middle_nodes":
241 interval_length = self.interval[1] - self.interval[0]
242 is_nodal_point = ~middle_node_flags
243 nodal_evaluation_points = evaluation_points[is_nodal_point]
245 n_el = len(nodal_evaluation_points) - 1
246 middle_nodes = int(((len(evaluation_points) - 1) / n_el) - 1)
248 t_evaluate = _np.zeros_like(evaluation_points)
249 S_evaluate = _np.zeros_like(evaluation_points)
251 # Evaluate the nodal points based on direct parametric mapping.
252 t_evaluate[is_nodal_point] = (
253 self.interval[0] + nodal_evaluation_points * interval_length
254 )
255 S_evaluate[is_nodal_point] = self.S_from_t(t_evaluate[is_nodal_point])
257 # For the middle nodes, do an interpolation in arc-length space.
258 for i_interval in range(n_el):
259 eval_a = nodal_evaluation_points[i_interval]
260 eval_b = nodal_evaluation_points[i_interval + 1]
262 S_a = S_evaluate[i_interval * (middle_nodes + 1)]
263 S_b = S_evaluate[(i_interval + 1) * (middle_nodes + 1)]
265 for i_middle in range(middle_nodes):
266 index_node = i_interval * (middle_nodes + 1) + i_middle + 1
267 eval_middle_node = evaluation_points[index_node]
268 factor = (eval_middle_node - eval_a) / (eval_b - eval_a)
269 S = S_a + factor * (S_b - S_a)
270 t_evaluate[index_node] = self.t_from_S(S)
271 S_evaluate[index_node] = S
273 else:
274 raise ValueError(f"Unknown method {self.method} for arc-length evaluation!")
276 return (t_evaluate, S_evaluate)
279def create_beam_mesh_parametric_curve(
280 mesh: _Mesh,
281 beam_class: type[_Beam],
282 material: _MaterialBeamBase,
283 function: _Callable,
284 interval: tuple[float, float],
285 *,
286 output_length: bool | None = False,
287 function_derivative: _Callable | None = None,
288 function_rotation: _Callable | None = None,
289 vectorized: bool = False,
290 arc_length_integrator_kwargs: dict | None = None,
291 **kwargs,
292) -> _GeometryName | tuple[_GeometryName, float]:
293 """Generate a beam from a parametric curve.
295 Integration along the beam is performed with scipy, and if the gradient is
296 not explicitly provided, it is calculated with the numpy wrapper autograd.
298 Args
299 ----
300 mesh: Mesh
301 Mesh that the curve will be added to.
302 beam_class: Beam
303 Class of beam that will be used for this line.
304 material: Material
305 Material for this line.
306 function: function
307 3D-parametric curve that represents the beam axis. If only a 2D
308 point is returned, the triad creation is simplified. If
309 mathematical functions are used, they have to come from the wrapper
310 autograd.numpy.
311 interval: [start end]
312 Start and end values for the parameter of the curve, must be in ascending
313 order.
314 output_length: bool
315 If this is true, the function returns a tuple containing the created
316 sets and the total arc length along the integrated function.
317 function_derivative: function -> R3
318 Explicitly provide the jacobian of the centerline position.
319 function_rotation: function -> Rotation
320 If this argument is given, the triads are computed with this
321 function, on the same interval as the position function. Must
322 return a Rotation object.
323 If no function_rotation is given, the rotation of the first node
324 is calculated automatically and all subsequent nodal rotations
325 are calculated based on a smallest rotation mapping onto the curve
326 tangent vector.
327 vectorized:
328 If true, the function and function_derivative are assumed to be
329 vectorized, i.e., they can take arrays of parameter values and return
330 arrays of positions/tangents.
331 arc_length_integrator_kwargs:
332 Additional arguments for the arc-length integrator.
334 **kwargs (for all of them look into create_beam_mesh_function)
335 ----
336 n_el: int
337 Number of equally spaced beam elements along the line. Defaults to 1.
338 Mutually exclusive with l_el.
339 l_el: float
340 Desired length of beam elements. Mutually exclusive with n_el.
341 Be aware, that this length might not be achieved, if the elements are
342 warped after they are created.
344 Return
345 ----
346 return_set: GeometryName
347 Set with the 'start' and 'end' node of the curve. Also a 'line' set
348 with all nodes of the curve.
349 """
350 # Set default values for optional arguments.
351 if arc_length_integrator_kwargs is None:
352 arc_length_integrator_kwargs = {}
354 # To avoid issues with automatic differentiation, we need to ensure that the interval
355 # values are of type float.
356 interval_array = _np.asarray(interval, dtype=float)
358 # Validate interval shape, length and order.
359 if interval_array.ndim != 1:
360 raise ValueError(
361 f"Interval must be a 1D sequence of exactly two values, got array with shape {interval_array.shape}."
362 )
363 if interval_array.size != 2:
364 raise ValueError(
365 f"Interval must contain exactly two values, got {interval_array.size}."
366 )
367 if interval_array[0] >= interval_array[1]:
368 raise ValueError(f"Interval must be in ascending order, got {interval_array}.")
370 # Ensure that the function can be evaluated for multiple values of t at once.
371 if not vectorized:
372 original_function = function
374 def function(t_array):
375 """Perform a vectorized evaluation of the function."""
376 return _np.array([original_function(t) for t in t_array])
378 if function_derivative is not None:
379 original_function_derivative = function_derivative
380 else:
381 # If no function derivative is given, we can use autograd to compute it. We need to make sure that the
382 # autograd jacobian can also handle vectorized inputs.
383 original_function_derivative = _jacobian(original_function)
385 def function_derivative(t_array):
386 """Perform a vectorized evaluation of the function derivative."""
387 return _np.array([original_function_derivative(t) for t in t_array])
389 if function_derivative is None:
390 raise ValueError(
391 "Function derivative could not be determined! For vectorized inputs, "
392 "the function derivative must be explicitly provided."
393 )
395 # Check size and type of position function
396 test_evaluation_of_function = function([interval_array[0]])[0]
398 if len(test_evaluation_of_function) == 2:
399 is_3d_curve = False
400 elif len(test_evaluation_of_function) == 3:
401 is_3d_curve = True
402 else:
403 raise ValueError("Function must return either 2d or 3d curve!")
405 # Setup the arc length integration object.
406 arc_length_evaluator = _ArcLengthEvaluation(
407 function_derivative, interval_array, **arc_length_integrator_kwargs
408 )
410 class _BeamFunctionGenerator:
411 """This class manages the creation the actual beam nodes and rotations."""
413 def __init__(
414 self,
415 interval_array: _np.ndarray,
416 function: _Callable,
417 function_derivative: _Callable,
418 function_rotation: _Callable | None,
419 is_3d_curve: bool,
420 ):
421 """Initialize the object and define a starting triad."""
422 self.interval_array = interval_array
423 self.function = function
424 self.function_derivative = function_derivative
425 self.function_rotation = function_rotation
426 self.is_3d_curve = is_3d_curve
428 if self.is_3d_curve:
429 r_prime = self.function_derivative([self.interval_array[0]])[0]
430 if abs(_np.dot(r_prime, [0, 0, 1])) < abs(_np.dot(r_prime, [0, 1, 0])):
431 t2_temp = [0, 0, 1]
432 else:
433 t2_temp = [0, 1, 0]
434 self.last_created_triad = _Rotation.from_basis(r_prime, t2_temp)
436 def evaluate_positions_and_rotations(
437 self, evaluation_positions: _np.ndarray, middle_node_flags: _np.ndarray
438 ) -> tuple[_np.ndarray, list[_Rotation], _np.ndarray]:
439 """This function evaluates the positions and rotations for given node
440 positions within the interval [0,1].
442 Args:
443 evaluation_positions:
444 Node positions in the interval [0, 1]. Depending on the integration method,
445 this can be in normalized parameter space or arc-length space.
446 middle_node_flags:
447 Boolean array indicating which of the node positions
448 are middle nodes (True) and which are inter element nodes (False).
450 Returns:
451 coordinates:
452 Physical coordinates of all nodes points along the curve.
453 rotations:
454 Rotations at all nodes points along the curve.
455 """
456 # Get the nodal parameter coordinates and the nodal arc-lengths.
457 t_evaluate, S_evaluate = arc_length_evaluator.evaluate_all(
458 evaluation_positions, middle_node_flags
459 )
461 # Position at S.
462 coordinates = _np.zeros((len(evaluation_positions), 3))
463 function_evaluation = self.function(t_evaluate)
464 if self.is_3d_curve:
465 coordinates[:, :] = function_evaluation
466 else:
467 coordinates[:, :2] = function_evaluation
469 # Rotation at S.
470 rotations = []
471 if self.function_rotation is not None:
472 for t in t_evaluate:
473 rotations.append(self.function_rotation(t))
475 else:
476 tangents = self.function_derivative(t_evaluate)
477 for r_prime in tangents:
478 if self.is_3d_curve:
479 # Create the next triad via the smallest rotation mapping based
480 # on the last triad.
481 rot = _smallest_rotation(self.last_created_triad, r_prime)
482 self.last_created_triad = rot.copy()
483 else:
484 # The rotation simplifies in the 2d case.
485 rot = _Rotation([0, 0, 1], _np.arctan2(r_prime[1], r_prime[0]))
486 rotations.append(rot)
488 # Return the needed values for beam creation.
489 return (coordinates, rotations, S_evaluate)
491 # Create the beam in the mesh
492 created_sets = _create_beam_mesh_generic(
493 mesh,
494 beam_class=beam_class,
495 material=material,
496 beam_function_evaluate_positions_and_rotations=True,
497 beam_function=_BeamFunctionGenerator(
498 interval_array,
499 function,
500 function_derivative,
501 function_rotation,
502 is_3d_curve,
503 ),
504 interval=(0.0, 1.0),
505 interval_length=arc_length_evaluator.approximate_total_arc_length(),
506 **kwargs,
507 )
509 if output_length:
510 return (created_sets, arc_length_evaluator.get_total_arc_length())
511 else:
512 return created_sets