Coverage for src/beamme/mesh_creation_functions/beam_generic.py: 94%
146 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"""Generic function for beam creation."""
24from collections.abc import Callable as _Callable
25from typing import Any as _Any
27import numpy as _np
29from beamme.core.conf import bme as _bme
30from beamme.core.element_beam import Beam as _Beam
31from beamme.core.geometry_set import GeometryName as _GeometryName
32from beamme.core.geometry_set import GeometrySet as _GeometrySet
33from beamme.core.material import MaterialBeamBase as _MaterialBeamBase
34from beamme.core.mesh import Mesh as _Mesh
35from beamme.core.node import NodeCosserat as _NodeCosserat
36from beamme.core.rotation import Rotation as _Rotation
37from beamme.utils.nodes import get_single_node as _get_single_node
40def _get_interval_node_positions_of_elements(
41 interval: tuple[float, float],
42 n_el: int | None,
43 l_el: float | None,
44 node_positions_of_elements: list[float] | None,
45 interval_length: float | None,
46) -> _np.ndarray:
47 """Get the node positions within the interval [0,1].
49 Args:
50 interval:
51 Start and end values for interval that will be used to create the
52 beam filament.
53 n_el:
54 Number of equally spaced beam elements along the line. Defaults to 1.
55 Mutually exclusive with l_el
56 l_el:
57 Desired length of beam elements. This requires the option interval_length
58 to be set. Mutually exclusive with n_el. Be aware, that this length
59 might not be achieved, if the elements are warped after they are
60 created.
61 node_positions_of_elements:
62 A list of normalized positions (within [0,1] and in ascending order)
63 that define the boundaries of beam elements along the created curve.
64 The given values will be mapped to the actual `interval` given as an
65 argument to this function. These values specify where elements start
66 and end, additional internal nodes (such as midpoints in higher-order
67 elements) may be placed automatically.
68 interval_length:
69 Approximation of the total length of the interval. Is required when
70 the option `l_el` is given.
72 Returns:
73 Numpy array with the node positions within the interval.
74 """
75 # Check for mutually exclusive parameters
76 n_given_arguments = sum(
77 1
78 for argument in [n_el, l_el, node_positions_of_elements]
79 if argument is not None
80 )
81 if n_given_arguments == 0:
82 # No arguments were given, use a single element per default
83 n_el = 1
84 elif n_given_arguments > 1:
85 raise ValueError(
86 'The arguments "n_el", "l_el" and "node_positions_of_elements" are mutually exclusive'
87 )
89 # Cases where we have equally spaced elements
90 if n_el is not None or l_el is not None:
91 if l_el is not None:
92 # Calculate the number of elements in case a desired element length is provided
93 if interval_length is None:
94 raise ValueError(
95 'The parameter "l_el" requires "interval_length" to be set.'
96 )
97 n_el = max([1, round(interval_length / l_el)])
98 elif n_el is None:
99 raise ValueError("n_el should not be None at this point")
101 node_positions_of_elements = [i_node / n_el for i_node in range(n_el + 1)]
102 # A list for the element node positions was provided
103 elif node_positions_of_elements is not None:
104 # Check that the given positions are in ascending order and start with 0 and end with 1
105 for index, value, name in zip([0, -1], [0, 1], ["First", "Last"]):
106 if not _np.isclose(
107 value,
108 node_positions_of_elements[index],
109 atol=1e-12,
110 rtol=0.0,
111 ):
112 raise ValueError(
113 f"{name} entry of node_positions_of_elements must be {value}, got {node_positions_of_elements[index]}"
114 )
115 if not all(
116 x < y
117 for x, y in zip(node_positions_of_elements, node_positions_of_elements[1:])
118 ):
119 raise ValueError(
120 f"The given node_positions_of_elements must be in ascending order. Got {node_positions_of_elements}"
121 )
122 else:
123 raise ValueError(
124 'One of the parameters "n_el", "l_el" or "node_positions_of_elements" has to be provided.'
125 )
127 return interval[0] + (interval[1] - interval[0]) * _np.asarray(
128 node_positions_of_elements
129 )
132def _get_interval_nodal_positions(
133 interval_node_positions_of_elements: _np.ndarray, nodes_create: list[float]
134) -> tuple[_np.ndarray, _np.ndarray]:
135 """Return the nodal positions along the interval depending on the element
136 formulation.
138 Args:
139 interval_node_positions_of_elements:
140 Numpy array with the element boundary node positions within the interval of the filament
141 nodes_create:
142 List with the FE parameter coordinates (in the interval [-1,1]) of the
143 element nodes.
145 Returns:
146 evaluation_positions:
147 Numpy array with the interval positions of all nodes along the beam filament (ordered).
148 middle_node_flags:
149 Numpy array with flags that indicate if a node is an element internal node.
150 """
151 if (
152 _np.abs(nodes_create[0] + 1.0) > _bme.eps_parameter_space
153 or _np.abs(nodes_create[-1] - 1.0) > _bme.eps_parameter_space
154 ):
155 raise ValueError(
156 "The first and last entry of nodes_create must be -1 and 1, respectively."
157 )
159 middle_node_coordinates = nodes_create[1:-1]
160 n_middle_nodes = len(middle_node_coordinates)
161 n_el = len(interval_node_positions_of_elements) - 1
162 n_nodes = n_el * n_middle_nodes + (n_el + 1)
164 evaluation_positions = _np.zeros(n_nodes)
165 evaluation_positions[:: n_middle_nodes + 1] = interval_node_positions_of_elements
167 interval_start_positions = interval_node_positions_of_elements[:-1]
168 interval_end_positions = interval_node_positions_of_elements[1:]
169 interval_length = interval_end_positions - interval_start_positions
171 for i in range(n_middle_nodes):
172 nodes_create_position = 0.5 * (middle_node_coordinates[i] + 1.0)
173 evaluation_positions[i + 1 :: n_middle_nodes + 1] = (
174 interval_start_positions + nodes_create_position * interval_length
175 )
177 middle_node_flags = _np.ones(n_nodes, dtype=bool)
178 middle_node_flags[:: n_middle_nodes + 1] = False
180 return evaluation_positions, middle_node_flags
183def _evaluate_positions_and_rotations(
184 beam_function: _Callable[[float], tuple[_np.ndarray, _Rotation, float | None]],
185 evaluation_positions: _np.ndarray,
186) -> tuple[_np.ndarray, list[_Rotation], _np.ndarray]:
187 """Evaluate positions, rotations and arc lengths along the filament, also return a
188 flag indicating middle nodes.
190 Args:
191 beam_function:
192 The `beam_function` has to take one variable s (from `evaluation_positions`)
193 and return the position, rotation and arc-length along the beam.
194 evaluation_positions:
195 Numpy array with the node positions within the interval of the filament.
197 Returns:
198 coordinates:
199 Numpy array with the coordinates of all nodes along the beam.
200 rotations:
201 List with the rotations of all nodes along the beam.
202 arc_lengths:
203 Numpy array with the arc lengths of all nodes along the beam.
204 """
205 n_nodes = len(evaluation_positions)
206 coordinates = _np.zeros((n_nodes, 3))
207 rotations: list[_Rotation] = []
208 arc_lengths = _np.zeros(n_nodes)
210 for i_node, evaluation_position in enumerate(evaluation_positions):
211 position, rotation, arc_length = beam_function(evaluation_position)
212 coordinates[i_node, :] = position
213 rotations.append(rotation)
214 arc_lengths[i_node] = arc_length
216 return coordinates, rotations, arc_lengths
219def _check_given_node_and_return_relative_twist(
220 mesh: _Mesh,
221 node: _NodeCosserat,
222 position_from_function: _np.ndarray,
223 rotation_from_function: _Rotation,
224 name: str,
225) -> _Rotation | None:
226 """Perform some checks for given nodes and return relative twist if necessary.
228 If the rotations do not match, check if the first basis vector of the triads is the same. If that is the case, a simple relative twist can be applied to ensure that the triad field is continuous. This relative twist can lead to issues if the beam cross-section is not double symmetric.
230 Args:
231 mesh: Mesh in which to check if the given nodes already exist.
232 node: Given node that should be used at the start or end of the beam.
233 position_from_function: Position at the start or end of the beam as given
234 by the beam function.
235 rotation_from_function: Rotation at the start or end of the beam as given
236 by the beam function.
237 name: Name of the node ("start" or "end") for better error messages.
239 Returns:
240 relative_twist:
241 If the rotation of the given node does not match with the one from the
242 function, but the tangent is the same, the relative twist that has to
243 be applied to the rotation field is returned. If no relative twist is
244 necessary, None is returned.
245 """
246 if node not in mesh.nodes:
247 raise ValueError("The given node is not in the current mesh")
249 if _np.linalg.norm(position_from_function - node.coordinates) > _bme.eps_pos:
250 raise ValueError(
251 f"The position of the given {name} node does not match with the position from the function!"
252 )
254 if rotation_from_function == node.rotation:
255 return None
256 elif not _bme.allow_beam_rotation:
257 raise ValueError(
258 f"The rotation of the given {name} node does not match with the rotation from the function!"
259 )
260 else:
261 # Evaluate the relative rotation
262 # First check if the first basis vector is the same
263 relative_basis_1 = node.rotation.inv() * rotation_from_function * [1, 0, 0]
264 if _np.linalg.norm(relative_basis_1 - [1, 0, 0]) < _bme.eps_quaternion:
265 # Calculate the relative rotation
266 return rotation_from_function.inv() * node.rotation
267 else:
268 raise ValueError(
269 f"The tangent of the {name} node does not match with the given function!"
270 )
273def create_beam_mesh_generic(
274 mesh: _Mesh,
275 *,
276 beam_class: type[_Beam],
277 material: _MaterialBeamBase,
278 beam_function: _Any,
279 interval: tuple[float, float],
280 beam_function_evaluate_positions_and_rotations: bool = False,
281 n_el: int | None = None,
282 l_el: float | None = None,
283 node_positions_of_elements: list[float] | None = None,
284 interval_length: float | None = None,
285 set_nodal_arc_length: bool = False,
286 nodal_arc_length_offset: float | None = None,
287 start_node: _NodeCosserat | _GeometrySet | None = None,
288 end_node: _NodeCosserat | _GeometrySet | None = None,
289 close_beam: bool = False,
290) -> _GeometryName:
291 """Generic beam creation function.
293 Remark for given start and/or end nodes:
294 If the rotation does not match, but the tangent vector is the same,
295 the created beams triads are rotated so the physical problem stays
296 the same (for axi-symmetric beam cross-sections) but the nodes can
297 be reused.
299 Args:
300 mesh:
301 Mesh that the created beam(s) should be added to.
302 beam_class:
303 Class of beam that will be used for this line.
304 material:
305 Material for this line.
306 beam_function:
307 The beam_function has to return the position along the beam centerline
308 for any point in the given `interval`.
310 Usually, the Jacobian of the returned position field should be a unit
311 vector. Otherwise, the nodes may be spaced in an undesired way. All
312 standard mesh creation functions fulfill this property.
313 interval:
314 Start and end values for interval that will be used to create the
315 beam.
316 beam_function_evaluate_positions_and_rotations:
317 Flag to indicate if the beam_function already provides an efficient
318 evaluation of all positions and rotations (and arc lengths) at once.
319 If this is True, the beam_function has to provide a method
320 `evaluate_positions_and_rotations(evaluation_positions, middle_node_flags)`
321 that returns the positions, rotations and arc lengths for all given
322 evaluation positions at once. This can speed up the creation of beams
323 significantly, especially for complex beam functions.
324 n_el:
325 Number of equally spaced beam elements along the line. Defaults to 1.
326 Mutually exclusive with l_el
327 l_el:
328 Desired length of beam elements. This requires the option `interval_length`
329 to be set. Mutually exclusive with n_el. Be aware, that this length
330 might not be achieved, if the elements are warped after they are
331 created.
332 node_positions_of_elements:
333 A list of normalized positions (within [0,1] and in ascending order)
334 that define the boundaries of beam elements along the created curve.
335 The given values will be mapped to the actual `interval` given as an
336 argument to this function. These values specify where elements start
337 and end, additional internal nodes (such as midpoints in higher-order
338 elements) may be placed automatically.
339 interval_length:
340 Approximation of the total length of the interval. Is required when
341 the option `l_el` is given.
342 set_nodal_arc_length:
343 Flag if the arc length along the beam filament is set in the created
344 nodes. It is ensured that the arc length is consistent with possible
345 given start/end nodes.
346 nodal_arc_length_offset:
347 Offset of the stored nodal arc length w.r.t. to the one generated by
348 the function. Defaults to 0, the arc length set in the start node, or
349 the arc length in the end node minus total length (such that the arc
350 length at the end node matches).
351 start_node:
352 Node to use as the first node for this line. Use this if the line
353 is connected to other lines (angles have to be the same, otherwise
354 connections should be used). If a geometry set is given, it can
355 contain one, and one node only.
356 end_node:
357 If this is a Node or GeometrySet, the last node of the created beam
358 is set to that node.
359 If it is True the created beam is closed within itself.
360 close_beam:
361 If it is True the created beam is closed within itself (mutually
362 exclusive with end_node).
364 Returns:
365 Geometry sets with the 'start' and 'end' node of the curve. Also a 'line' set
366 with all nodes of the curve.
367 """
368 if close_beam and end_node is not None:
369 raise ValueError(
370 'The arguments "close_beam" and "end_node" are mutually exclusive'
371 )
373 if set_nodal_arc_length:
374 if close_beam:
375 raise ValueError(
376 "The flags 'set_nodal_arc_length' and 'close_beam' are mutually exclusive."
377 )
378 elif nodal_arc_length_offset is not None:
379 raise ValueError(
380 'Providing the argument "nodal_arc_length_offset" without setting '
381 '"set_nodal_arc_length" to True does not make sense.'
382 )
384 # Get element boundary node positions within the given interval
385 interval_node_positions_of_elements = _get_interval_node_positions_of_elements(
386 interval, n_el, l_el, node_positions_of_elements, interval_length
387 )
388 n_el = len(interval_node_positions_of_elements) - 1
390 # Get the nodal positions in the interval for all nodes (depending on the element formulation).
391 evaluation_positions, middle_node_flags = _get_interval_nodal_positions(
392 interval_node_positions_of_elements, beam_class.nodes_create
393 )
395 # Evaluate the centerline position and the rotation for all beam nodes
396 if not beam_function_evaluate_positions_and_rotations:
397 coordinates, rotations, arc_lengths = _evaluate_positions_and_rotations(
398 beam_function, evaluation_positions
399 )
400 else:
401 coordinates, rotations, arc_lengths = (
402 beam_function.evaluate_positions_and_rotations(
403 evaluation_positions, middle_node_flags
404 )
405 )
407 # Make sure the material is in the mesh.
408 mesh.add_material(material)
410 # Inspect given nodes and get relative twists if necessary
411 relative_twist_start = None
412 if start_node is not None:
413 start_node = _get_single_node(start_node)
414 relative_twist_start = _check_given_node_and_return_relative_twist(
415 mesh, start_node, coordinates[0], rotations[0], "start"
416 )
418 # If an end node is given, check what behavior is wanted.
419 relative_twist_end = None
420 if end_node is not None:
421 end_node = _get_single_node(end_node)
422 relative_twist_end = _check_given_node_and_return_relative_twist(
423 mesh, end_node, coordinates[-1], rotations[-1], "end"
424 )
426 # Check if a relative twist has to be applied
427 relative_twist_list = [
428 twist
429 for twist in [relative_twist_start, relative_twist_end]
430 if twist is not None
431 ]
432 if len(relative_twist_list) == 2:
433 if not relative_twist_list[0] == relative_twist_list[1]:
434 raise ValueError(
435 "The relative twist required for the start and end node do not match"
436 )
437 if len(relative_twist_list) > 0:
438 relative_twist = relative_twist_list[0]
439 for i_rot, rotation in enumerate(rotations):
440 rotations[i_rot] = rotation * relative_twist
442 # Get the start value for the arc length functionality
443 if set_nodal_arc_length:
444 if nodal_arc_length_offset is not None:
445 # Let's use the given value, the later check will detect if this
446 # does not match the given nodes.
447 pass
448 elif start_node is not None and start_node.arc_length is not None:
449 nodal_arc_length_offset = start_node.arc_length
450 elif end_node is not None and end_node.arc_length is not None:
451 nodal_arc_length_offset = end_node.arc_length - arc_lengths[-1]
452 else:
453 # Default value
454 nodal_arc_length_offset = 0.0
455 arc_lengths += nodal_arc_length_offset
457 if start_node is not None:
458 if _np.abs(start_node.arc_length - arc_lengths[0]) > _bme.eps_pos:
459 raise ValueError(
460 "The arc length at the start node does not match with "
461 "the calculated one!"
462 )
463 if end_node is not None:
464 if _np.abs(end_node.arc_length - arc_lengths[-1]) > _bme.eps_pos:
465 raise ValueError(
466 "The arc length at the end node does not match with "
467 "the calculated one!"
468 )
469 else:
470 arc_lengths = [None] * len(arc_lengths)
472 # Create the nodes and add the new ones to the mesh
473 nodes = [
474 _NodeCosserat(pos, rot, is_middle_node=middle_node_flag, arc_length=arc_length)
475 for pos, rot, arc_length, middle_node_flag in zip(
476 coordinates, rotations, arc_lengths, middle_node_flags
477 )
478 ]
479 if start_node is not None:
480 nodes[0] = start_node
481 if close_beam:
482 nodes[-1] = nodes[0]
483 elif end_node is not None:
484 nodes[-1] = end_node
485 start_slice = 1 if start_node is not None else None
486 end_slice = -1 if end_node is not None or close_beam else None
487 mesh.nodes.extend(nodes[start_slice:end_slice])
489 # Create the beam elements and assign the nodes
490 nodes_per_element = len(beam_class.nodes_create)
491 elements: list[_Beam] = []
492 for i_el in range(n_el):
493 beam = beam_class(
494 material=material,
495 nodes=nodes[
496 i_el * (nodes_per_element - 1) : (i_el + 1) * (nodes_per_element - 1)
497 + 1
498 ],
499 )
500 elements.append(beam)
502 # Add items to the mesh
503 mesh.elements.extend(elements)
505 # Set the nodes that are at the beginning and end of line (for search
506 # of overlapping points)
507 nodes[0].is_end_node = True
508 nodes[-1].is_end_node = True
510 # Create geometry sets that will be returned.
511 return_set = _GeometryName()
512 return_set["start"] = _GeometrySet(nodes[0])
513 return_set["end"] = _GeometrySet(nodes[-1])
514 return_set["line"] = _GeometrySet(elements)
516 return return_set