Coverage for src / beamme / core / element_beam.py: 92%
128 statements
« prev ^ index » next coverage.py v7.12.0, created at 2025-11-21 12:57 +0000
1# The MIT License (MIT)
2#
3# Copyright (c) 2018-2025 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 defines the base beam element."""
24from typing import Any as _Any
25from typing import Callable as _Callable
26from typing import List as _List
27from typing import Optional as _Optional
28from typing import Type as _Type
30import numpy as _np
31import vtk as _vtk
33from beamme.core.conf import bme as _bme
34from beamme.core.element import Element as _Element
35from beamme.core.material import MaterialBeamBase as _MaterialBeamBase
36from beamme.core.node import NodeCosserat as _NodeCosserat
37from beamme.core.rotation import Rotation as _Rotation
38from beamme.core.vtk_writer import add_point_data_node_sets as _add_point_data_node_sets
41class Beam(_Element):
42 """A base class for a beam element."""
44 # An array that defines the parameter positions of the element nodes,
45 # in ascending order.
46 nodes_create: _Any = []
48 # A list of valid material types for this element.
49 valid_materials: _Any = []
51 def __init__(self, material=None, nodes=None):
52 super().__init__(nodes=nodes, material=material)
54 @classmethod
55 def create_beam(
56 cls: _Type["Beam"],
57 material: _MaterialBeamBase,
58 beam_function: _Callable,
59 *,
60 start_node: _Optional[_NodeCosserat] = None,
61 end_node: _Optional[_NodeCosserat] = None,
62 relative_twist: _Optional[_Rotation] = None,
63 set_nodal_arc_length: bool = False,
64 nodal_arc_length_offset: _Optional[float] = None,
65 ) -> tuple["Beam", _List[_NodeCosserat]]:
66 """Create the nodes for this beam element. The function returns a list
67 with the created nodes.
69 In the case of start_node and end_node, it is checked, that the
70 function and the node have the same coordinates and rotations.
72 Args:
73 material: The material to be used with this beam.
74 beam_function:
75 Returns the position, rotation and (optionally) arc length of the
76 beam along the local coordinate xi. If no arc lengths is provided,
77 the returned value should simply be None.
78 start_node:
79 If this argument is given, this is the node of the beam at xi=-1.
80 end_node:
81 If this argument is given, this is the node of the beam at xi=1.
82 relative_twist:
83 Apply this relative rotation to all created nodes. This can be used to
84 "twist" the created beam to match the rotation of given nodes.
85 set_nodal_arc_length:
86 Flag if the arc length in the created nodes should be set.
87 nodal_arc_length_offset:
88 Offset of the stored nodal arc length w.r.t. to the one generated
89 by the function.
91 Returns:
92 The created beam object and the newly created nodes. The created nodes
93 are returned, because we can then directly add them to a `Mesh` without
94 having to check if they are already in the mesh.
95 """
97 def check_node(node, pos, rot, arc_length, name):
98 """Check if the given node matches with the position and rotation
99 and optionally also the arc length."""
101 if _np.linalg.norm(pos - node.coordinates) > _bme.eps_pos:
102 raise ValueError(
103 f"{name} position does not match with function! Got {pos} from function but "
104 + f"given node value is {node.coordinates}"
105 )
106 if not node.rotation == rot:
107 raise ValueError(f"{name} rotation does not match with function!")
109 if arc_length is not None:
110 if _np.abs(node.arc_length - arc_length) > _bme.eps_pos:
111 raise ValueError(
112 f"Arc lengths don't match, got {node.arc_length} and {arc_length}"
113 )
115 # Flags if nodes are given
116 has_start_node = start_node is not None
117 has_end_node = end_node is not None
119 # Create the element
120 beam = cls(material=material)
122 # Loop over local nodes.
123 arc_length = None
124 for i, xi in enumerate(cls.nodes_create):
125 # Get the position and rotation at xi
126 pos, rot, arc_length_from_function = beam_function(xi)
127 if relative_twist is not None:
128 rot = rot * relative_twist
129 if set_nodal_arc_length:
130 arc_length = arc_length_from_function + nodal_arc_length_offset
132 # Check if the position and rotation match existing nodes
133 if i == 0 and has_start_node:
134 check_node(start_node, pos, rot, arc_length, "start_node")
135 beam.nodes = [start_node]
136 elif (i == len(cls.nodes_create) - 1) and has_end_node:
137 check_node(end_node, pos, rot, arc_length, "end_node")
139 # Create the node
140 if (i > 0 or not has_start_node) and (
141 i < len(cls.nodes_create) - 1 or not has_end_node
142 ):
143 is_middle_node = 0 < i < len(cls.nodes_create) - 1
144 beam.nodes.append(
145 _NodeCosserat(
146 pos, rot, is_middle_node=is_middle_node, arc_length=arc_length
147 )
148 )
150 # Get a list with the created nodes.
151 if has_start_node:
152 created_nodes = beam.nodes[1:]
153 else:
154 created_nodes = beam.nodes.copy()
156 # Add the end node to the beam.
157 if has_end_node:
158 beam.nodes.append(end_node)
160 # Return the created beam object and the created nodes.
161 return beam, created_nodes
163 @classmethod
164 def get_coupling_dict(cls, coupling_dof_type):
165 """Return the dict to couple this beam to another beam."""
167 match coupling_dof_type:
168 case _bme.coupling_dof.joint:
169 if cls.coupling_joint_dict is None:
170 raise ValueError(f"Joint coupling is not implemented for {cls}")
171 return cls.coupling_joint_dict
172 case _bme.coupling_dof.fix:
173 if cls.coupling_fix_dict is None:
174 raise ValueError("Fix coupling is not implemented for {cls}")
175 return cls.coupling_fix_dict
176 case _:
177 raise ValueError(
178 f'Coupling_dof_type "{coupling_dof_type}" is not implemented!'
179 )
181 def flip(self):
182 """Reverse the nodes of this element.
184 This is usually used when reflected.
185 """
186 self.nodes = [self.nodes[-1 - i] for i in range(len(self.nodes))]
188 def _check_material(self):
189 """Check if the linked material is valid for this type of beam
190 element."""
191 for material_type in type(self).valid_materials:
192 if isinstance(self.material, material_type):
193 break
194 else:
195 raise TypeError(
196 f"Beam of type {type(self)} can not have a material of type {type(self.material)}!"
197 )
199 def get_vtk(
200 self,
201 vtk_writer_beam,
202 vtk_writer_solid,
203 *,
204 beam_centerline_visualization_segments=1,
205 **kwargs,
206 ):
207 """Add the representation of this element to the VTK writer as a poly
208 line.
210 Args
211 ----
212 vtk_writer_beam:
213 VTK writer for the beams.
214 vtk_writer_solid:
215 VTK writer for solid elements, not used in this method.
216 beam_centerline_visualization_segments: int
217 Number of segments to be used for visualization of beam centerline between successive
218 nodes. Default is 1, which means a straight line is drawn between the beam nodes. For
219 Values greater than 1, a Hermite interpolation of the centerline is assumed for
220 visualization purposes.
221 """
223 n_nodes = len(self.nodes)
224 n_segments = n_nodes - 1
225 n_additional_points_per_segment = beam_centerline_visualization_segments - 1
226 # Number of points (in addition to the nodes) to be used for output
227 n_additional_points = n_segments * n_additional_points_per_segment
228 n_points = n_nodes + n_additional_points
230 # Dictionary with cell data.
231 cell_data = self.vtk_cell_data.copy()
232 if self.material.radius is not None:
233 cell_data["cross_section_radius"] = self.material.radius
235 # Dictionary with point data.
236 point_data = {}
237 point_data["node_value"] = _np.zeros(n_points)
238 point_data["base_vector_1"] = _np.zeros((n_points, 3))
239 point_data["base_vector_2"] = _np.zeros((n_points, 3))
240 point_data["base_vector_3"] = _np.zeros((n_points, 3))
242 coordinates = _np.zeros((n_points, 3))
243 nodal_rotation_matrices = [
244 node.rotation.get_rotation_matrix() for node in self.nodes
245 ]
247 for i_node, (node, rotation_matrix) in enumerate(
248 zip(self.nodes, nodal_rotation_matrices)
249 ):
250 coordinates[i_node, :] = node.coordinates
251 if node.is_middle_node:
252 point_data["node_value"][i_node] = 0.5
253 else:
254 point_data["node_value"][i_node] = 1.0
256 point_data["base_vector_1"][i_node] = rotation_matrix[:, 0]
257 point_data["base_vector_2"][i_node] = rotation_matrix[:, 1]
258 point_data["base_vector_3"][i_node] = rotation_matrix[:, 2]
260 # We can output the element partner index, which is a debug quantity to help find elements
261 # with matching middle nodes. This is usually an indicator for an issue with the mesh.
262 # TODO: Check if we really need this partner index output any more
263 element_partner_indices = list(
264 [
265 node.element_partner_index
266 for node in self.nodes
267 if node.element_partner_index is not None
268 ]
269 )
270 if len(element_partner_indices) > 1:
271 raise ValueError(
272 "More than one element partner indices are currently not supported in the output"
273 "functionality"
274 )
275 elif len(element_partner_indices) == 1:
276 cell_data["partner_index"] = element_partner_indices[0] + 1
278 # Check if we have everything we need to write output or if we need to calculate additional
279 # points for a smooth beam visualization.
280 if beam_centerline_visualization_segments == 1:
281 point_connectivity = _np.arange(n_nodes)
282 else:
283 # We need the centerline shape function matrices, so calculate them once and use for
284 # all segments that we need. Drop the first and last value, since they represent the
285 # nodes which we have already added above.
286 xi = _np.linspace(-1, 1, beam_centerline_visualization_segments + 1)[1:-1]
287 hermite_shape_functions_pos = _np.array(
288 [
289 0.25 * (2.0 + xi) * (1.0 - xi) ** 2,
290 0.25 * (2.0 - xi) * (1.0 + xi) ** 2,
291 ]
292 ).transpose()
293 hermite_shape_functions_tan = _np.array(
294 [
295 0.125 * (1.0 + xi) * (1.0 - xi) ** 2,
296 -0.125 * (1.0 - xi) * (1.0 + xi) ** 2,
297 ]
298 ).transpose()
300 point_connectivity = _np.zeros(n_points, dtype=int)
302 for i_segment in range(n_segments):
303 positions = _np.array(
304 [
305 self.nodes[i_node].coordinates
306 for i_node in [i_segment, i_segment + 1]
307 ]
308 )
309 tangents = _np.array(
310 [
311 nodal_rotation_matrices[i_node][:, 0]
312 for i_node in [i_segment, i_segment + 1]
313 ]
314 )
315 length_factor = _np.linalg.norm(positions[1] - positions[0])
316 interpolated_coordinates = _np.dot(
317 hermite_shape_functions_pos, positions
318 ) + length_factor * _np.dot(hermite_shape_functions_tan, tangents)
320 index_first_point = (
321 n_nodes + i_segment * n_additional_points_per_segment
322 )
323 index_last_point = (
324 n_nodes + (i_segment + 1) * n_additional_points_per_segment
325 )
327 coordinates[index_first_point:index_last_point] = (
328 interpolated_coordinates
329 )
330 point_connectivity[
331 i_segment * beam_centerline_visualization_segments
332 ] = i_segment
333 point_connectivity[
334 (i_segment + 1) * beam_centerline_visualization_segments
335 ] = i_segment + 1
336 point_connectivity[
337 i_segment * beam_centerline_visualization_segments + 1 : (
338 i_segment + 1
339 )
340 * beam_centerline_visualization_segments
341 ] = _np.arange(index_first_point, index_last_point)
343 # Get the point data sets and add everything to the output file.
344 _add_point_data_node_sets(
345 point_data, self.nodes, extra_points=n_additional_points
346 )
347 indices = vtk_writer_beam.add_points(coordinates, point_data=point_data)
348 vtk_writer_beam.add_cell(
349 _vtk.vtkPolyLine, indices[point_connectivity], cell_data=cell_data
350 )
353def generate_beam_class(n_nodes: int):
354 """Return a class representing a general beam with n_nodes in BeamMe.
356 Args:
357 n_nodes: Number of equally spaced nodes along the beam centerline.
359 Returns:
360 A beam object that has n_nodes along the centerline.
361 """
363 # Define the class variable responsible for creating the nodes.
364 nodes_create = _np.linspace(-1, 1, num=n_nodes)
366 # Create the beam class which inherits from the base beam class.
367 return type(f"Beam{n_nodes}", (Beam,), {"nodes_create": nodes_create})
370Beam2 = generate_beam_class(2)
371Beam3 = generate_beam_class(3)
372Beam4 = generate_beam_class(4)
373Beam5 = generate_beam_class(5)