Coverage for src / beamme / core / element_beam.py: 91%
87 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"""This file defines the base beam element."""
24from typing import Any as _Any
26import numpy as _np
27import vtk as _vtk
29from beamme.core.conf import bme as _bme
30from beamme.core.element import Element as _Element
31from beamme.core.vtk_writer import add_point_data_node_sets as _add_point_data_node_sets
34class Beam(_Element):
35 """A base class for a beam element."""
37 # An array that defines the parameter positions of the element nodes,
38 # in ascending order.
39 nodes_create: _Any = []
41 # A list of valid material types for this element.
42 valid_materials: _Any = []
44 def __init__(self, material=None, nodes=None):
45 super().__init__(nodes=nodes, material=material)
47 @classmethod
48 def get_coupling_dict(cls, coupling_dof_type):
49 """Return the dict to couple this beam to another beam."""
51 match coupling_dof_type:
52 case _bme.coupling_dof.joint:
53 if cls.coupling_joint_dict is None:
54 raise ValueError(f"Joint coupling is not implemented for {cls}")
55 return cls.coupling_joint_dict
56 case _bme.coupling_dof.fix:
57 if cls.coupling_fix_dict is None:
58 raise ValueError("Fix coupling is not implemented for {cls}")
59 return cls.coupling_fix_dict
60 case _:
61 raise ValueError(
62 f'Coupling_dof_type "{coupling_dof_type}" is not implemented!'
63 )
65 def flip(self):
66 """Reverse the nodes of this element.
68 This is usually used when reflected.
69 """
70 self.nodes = [self.nodes[-1 - i] for i in range(len(self.nodes))]
72 def _check_material(self):
73 """Check if the linked material is valid for this type of beam
74 element."""
75 for material_type in type(self).valid_materials:
76 if isinstance(self.material, material_type):
77 break
78 else:
79 raise TypeError(
80 f"Beam of type {type(self)} can not have a material of type {type(self.material)}!"
81 )
83 def get_vtk(
84 self,
85 vtk_writer_beam,
86 vtk_writer_solid,
87 *,
88 beam_centerline_visualization_segments=1,
89 **kwargs,
90 ):
91 """Add the representation of this element to the VTK writer as a poly
92 line.
94 Args
95 ----
96 vtk_writer_beam:
97 VTK writer for the beams.
98 vtk_writer_solid:
99 VTK writer for solid elements, not used in this method.
100 beam_centerline_visualization_segments: int
101 Number of segments to be used for visualization of beam centerline between successive
102 nodes. Default is 1, which means a straight line is drawn between the beam nodes. For
103 Values greater than 1, a Hermite interpolation of the centerline is assumed for
104 visualization purposes.
105 """
107 n_nodes = len(self.nodes)
108 n_segments = n_nodes - 1
109 n_additional_points_per_segment = beam_centerline_visualization_segments - 1
110 # Number of points (in addition to the nodes) to be used for output
111 n_additional_points = n_segments * n_additional_points_per_segment
112 n_points = n_nodes + n_additional_points
114 # Dictionary with cell data.
115 cell_data = self.vtk_cell_data.copy()
116 if self.material.radius is not None:
117 cell_data["cross_section_radius"] = self.material.radius
119 # Dictionary with point data.
120 point_data = {}
121 point_data["node_value"] = _np.zeros(n_points)
122 point_data["base_vector_1"] = _np.zeros((n_points, 3))
123 point_data["base_vector_2"] = _np.zeros((n_points, 3))
124 point_data["base_vector_3"] = _np.zeros((n_points, 3))
126 coordinates = _np.zeros((n_points, 3))
127 nodal_rotation_matrices = [
128 node.rotation.get_rotation_matrix() for node in self.nodes
129 ]
131 for i_node, (node, rotation_matrix) in enumerate(
132 zip(self.nodes, nodal_rotation_matrices)
133 ):
134 coordinates[i_node, :] = node.coordinates
135 if node.is_middle_node:
136 point_data["node_value"][i_node] = 0.5
137 else:
138 point_data["node_value"][i_node] = 1.0
140 point_data["base_vector_1"][i_node] = rotation_matrix[:, 0]
141 point_data["base_vector_2"][i_node] = rotation_matrix[:, 1]
142 point_data["base_vector_3"][i_node] = rotation_matrix[:, 2]
144 # We can output the element partner index, which is a debug quantity to help find elements
145 # with matching middle nodes. This is usually an indicator for an issue with the mesh.
146 # TODO: Check if we really need this partner index output any more
147 element_partner_indices = list(
148 [
149 node.element_partner_index
150 for node in self.nodes
151 if node.element_partner_index is not None
152 ]
153 )
154 if len(element_partner_indices) > 1:
155 raise ValueError(
156 "More than one element partner indices are currently not supported in the output"
157 "functionality"
158 )
159 elif len(element_partner_indices) == 1:
160 cell_data["partner_index"] = element_partner_indices[0] + 1
162 # Check if we have everything we need to write output or if we need to calculate additional
163 # points for a smooth beam visualization.
164 if beam_centerline_visualization_segments == 1:
165 point_connectivity = _np.arange(n_nodes)
166 else:
167 # We need the centerline shape function matrices, so calculate them once and use for
168 # all segments that we need. Drop the first and last value, since they represent the
169 # nodes which we have already added above.
170 xi = _np.linspace(-1, 1, beam_centerline_visualization_segments + 1)[1:-1]
171 hermite_shape_functions_pos = _np.array(
172 [
173 0.25 * (2.0 + xi) * (1.0 - xi) ** 2,
174 0.25 * (2.0 - xi) * (1.0 + xi) ** 2,
175 ]
176 ).transpose()
177 hermite_shape_functions_tan = _np.array(
178 [
179 0.125 * (1.0 + xi) * (1.0 - xi) ** 2,
180 -0.125 * (1.0 - xi) * (1.0 + xi) ** 2,
181 ]
182 ).transpose()
184 point_connectivity = _np.zeros(n_points, dtype=int)
186 for i_segment in range(n_segments):
187 positions = _np.array(
188 [
189 self.nodes[i_node].coordinates
190 for i_node in [i_segment, i_segment + 1]
191 ]
192 )
193 tangents = _np.array(
194 [
195 nodal_rotation_matrices[i_node][:, 0]
196 for i_node in [i_segment, i_segment + 1]
197 ]
198 )
199 length_factor = _np.linalg.norm(positions[1] - positions[0])
200 interpolated_coordinates = _np.dot(
201 hermite_shape_functions_pos, positions
202 ) + length_factor * _np.dot(hermite_shape_functions_tan, tangents)
204 index_first_point = (
205 n_nodes + i_segment * n_additional_points_per_segment
206 )
207 index_last_point = (
208 n_nodes + (i_segment + 1) * n_additional_points_per_segment
209 )
211 coordinates[index_first_point:index_last_point] = (
212 interpolated_coordinates
213 )
214 point_connectivity[
215 i_segment * beam_centerline_visualization_segments
216 ] = i_segment
217 point_connectivity[
218 (i_segment + 1) * beam_centerline_visualization_segments
219 ] = i_segment + 1
220 point_connectivity[
221 i_segment * beam_centerline_visualization_segments + 1 : (
222 i_segment + 1
223 )
224 * beam_centerline_visualization_segments
225 ] = _np.arange(index_first_point, index_last_point)
227 # Get the point data sets and add everything to the output file.
228 _add_point_data_node_sets(
229 point_data, self.nodes, extra_points=n_additional_points
230 )
231 indices = vtk_writer_beam.add_points(coordinates, point_data=point_data)
232 vtk_writer_beam.add_cell(
233 _vtk.vtkPolyLine, indices[point_connectivity], cell_data=cell_data
234 )
237def generate_beam_class(n_nodes: int):
238 """Return a class representing a general beam with n_nodes in BeamMe.
240 Args:
241 n_nodes: Number of equally spaced nodes along the beam centerline.
243 Returns:
244 A beam object that has n_nodes along the centerline.
245 """
247 # Define the class variable responsible for creating the nodes.
248 nodes_create = _np.linspace(-1, 1, num=n_nodes)
250 # Create the beam class which inherits from the base beam class.
251 return type(f"Beam{n_nodes}", (Beam,), {"nodes_create": nodes_create})
254Beam2 = generate_beam_class(2)
255Beam3 = generate_beam_class(3)
256Beam4 = generate_beam_class(4)
257Beam5 = generate_beam_class(5)