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