Coverage for src/beamme/core/element_beam.py: 92%
127 statements
« prev ^ index » next coverage.py v7.10.7, created at 2025-09-29 11:30 +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 cell_data["cross_section_radius"] = self.material.radius
234 # Dictionary with point data.
235 point_data = {}
236 point_data["node_value"] = _np.zeros(n_points)
237 point_data["base_vector_1"] = _np.zeros((n_points, 3))
238 point_data["base_vector_2"] = _np.zeros((n_points, 3))
239 point_data["base_vector_3"] = _np.zeros((n_points, 3))
241 coordinates = _np.zeros((n_points, 3))
242 nodal_rotation_matrices = [
243 node.rotation.get_rotation_matrix() for node in self.nodes
244 ]
246 for i_node, (node, rotation_matrix) in enumerate(
247 zip(self.nodes, nodal_rotation_matrices)
248 ):
249 coordinates[i_node, :] = node.coordinates
250 if node.is_middle_node:
251 point_data["node_value"][i_node] = 0.5
252 else:
253 point_data["node_value"][i_node] = 1.0
255 point_data["base_vector_1"][i_node] = rotation_matrix[:, 0]
256 point_data["base_vector_2"][i_node] = rotation_matrix[:, 1]
257 point_data["base_vector_3"][i_node] = rotation_matrix[:, 2]
259 # We can output the element partner index, which is a debug quantity to help find elements
260 # with matching middle nodes. This is usually an indicator for an issue with the mesh.
261 # TODO: Check if we really need this partner index output any more
262 element_partner_indices = list(
263 [
264 node.element_partner_index
265 for node in self.nodes
266 if node.element_partner_index is not None
267 ]
268 )
269 if len(element_partner_indices) > 1:
270 raise ValueError(
271 "More than one element partner indices are currently not supported in the output"
272 "functionality"
273 )
274 elif len(element_partner_indices) == 1:
275 cell_data["partner_index"] = element_partner_indices[0] + 1
277 # Check if we have everything we need to write output or if we need to calculate additional
278 # points for a smooth beam visualization.
279 if beam_centerline_visualization_segments == 1:
280 point_connectivity = _np.arange(n_nodes)
281 else:
282 # We need the centerline shape function matrices, so calculate them once and use for
283 # all segments that we need. Drop the first and last value, since they represent the
284 # nodes which we have already added above.
285 xi = _np.linspace(-1, 1, beam_centerline_visualization_segments + 1)[1:-1]
286 hermite_shape_functions_pos = _np.array(
287 [
288 0.25 * (2.0 + xi) * (1.0 - xi) ** 2,
289 0.25 * (2.0 - xi) * (1.0 + xi) ** 2,
290 ]
291 ).transpose()
292 hermite_shape_functions_tan = _np.array(
293 [
294 0.125 * (1.0 + xi) * (1.0 - xi) ** 2,
295 -0.125 * (1.0 - xi) * (1.0 + xi) ** 2,
296 ]
297 ).transpose()
299 point_connectivity = _np.zeros(n_points, dtype=int)
301 for i_segment in range(n_segments):
302 positions = _np.array(
303 [
304 self.nodes[i_node].coordinates
305 for i_node in [i_segment, i_segment + 1]
306 ]
307 )
308 tangents = _np.array(
309 [
310 nodal_rotation_matrices[i_node][:, 0]
311 for i_node in [i_segment, i_segment + 1]
312 ]
313 )
314 length_factor = _np.linalg.norm(positions[1] - positions[0])
315 interpolated_coordinates = _np.dot(
316 hermite_shape_functions_pos, positions
317 ) + length_factor * _np.dot(hermite_shape_functions_tan, tangents)
319 index_first_point = (
320 n_nodes + i_segment * n_additional_points_per_segment
321 )
322 index_last_point = (
323 n_nodes + (i_segment + 1) * n_additional_points_per_segment
324 )
326 coordinates[index_first_point:index_last_point] = (
327 interpolated_coordinates
328 )
329 point_connectivity[
330 i_segment * beam_centerline_visualization_segments
331 ] = i_segment
332 point_connectivity[
333 (i_segment + 1) * beam_centerline_visualization_segments
334 ] = i_segment + 1
335 point_connectivity[
336 i_segment * beam_centerline_visualization_segments + 1 : (
337 i_segment + 1
338 )
339 * beam_centerline_visualization_segments
340 ] = _np.arange(index_first_point, index_last_point)
342 # Get the point data sets and add everything to the output file.
343 _add_point_data_node_sets(
344 point_data, self.nodes, extra_points=n_additional_points
345 )
346 indices = vtk_writer_beam.add_points(coordinates, point_data=point_data)
347 vtk_writer_beam.add_cell(
348 _vtk.vtkPolyLine, indices[point_connectivity], cell_data=cell_data
349 )
352def generate_beam_class(n_nodes: int):
353 """Return a class representing a general beam with n_nodes in BeamMe.
355 Args:
356 n_nodes: Number of equally spaced nodes along the beam centerline.
358 Returns:
359 A beam object that has n_nodes along the centerline.
360 """
362 # Define the class variable responsible for creating the nodes.
363 nodes_create = _np.linspace(-1, 1, num=n_nodes)
365 # Create the beam class which inherits from the base beam class.
366 return type(f"Beam{n_nodes}", (Beam,), {"nodes_create": nodes_create})
369Beam2 = generate_beam_class(2)
370Beam3 = generate_beam_class(3)
371Beam4 = generate_beam_class(4)
372Beam5 = generate_beam_class(5)