Coverage for src/meshpy/core/element_beam.py: 91%
120 statements
« prev ^ index » next coverage.py v7.8.0, created at 2025-04-28 04:21 +0000
1# The MIT License (MIT)
2#
3# Copyright (c) 2018-2025 MeshPy 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 in MeshPy."""
24from typing import Any as _Any
25from typing import Callable as _Callable
26from typing import List as _List
27from typing import Optional as _Optional
29import numpy as _np
30import vtk as _vtk
32from meshpy.core.conf import mpy as _mpy
33from meshpy.core.element import Element as _Element
34from meshpy.core.node import NodeCosserat as _NodeCosserat
35from meshpy.core.rotation import Rotation as _Rotation
36from meshpy.core.vtk_writer import add_point_data_node_sets as _add_point_data_node_sets
39class Beam(_Element):
40 """A base class for a beam element."""
42 # An array that defines the parameter positions of the element nodes,
43 # in ascending order.
44 nodes_create: _Any = []
46 # A list of valid material types for this element.
47 valid_material: _Any = []
49 # Coupling strings.
50 coupling_fix_string: _Optional[str] = None
51 coupling_joint_string: _Optional[str] = None
53 def __init__(self, material=None, nodes=None):
54 super().__init__(nodes=nodes, material=material)
56 def create_beam(
57 self,
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 ) -> _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 beam_function: function(xi)
74 Returns the position, rotation and (optionally) arc length of the
75 beam along the local coordinate xi. If no arc lengths is provided,
76 the returned value should simply be None.
77 start_node: Node
78 If this argument is given, this is the node of the beam at xi=-1.
79 end_node: Node
80 If this argument is given, this is the node of the beam at xi=1.
81 relative_twist: Rotation
82 Apply this relative rotation to all created nodes. This can be used to
83 "twist" the created beam to match the rotation of given nodes.
84 set_nodal_arc_length:
85 Flag if the arc length in the created nodes should be set.
86 nodal_arc_length_offset:
87 Offset of the stored nodal arc length w.r.t. to the one generated
88 by the function.
89 """
91 if len(self.nodes) > 0:
92 raise ValueError("The beam should not have any local nodes yet!")
94 def check_node(node, pos, rot, arc_length, name):
95 """Check if the given node matches with the position and rotation
96 and optionally also the arc length."""
98 if _np.linalg.norm(pos - node.coordinates) > _mpy.eps_pos:
99 raise ValueError(
100 f"{name} position does not match with function! Got {pos} from function but "
101 + f"given node value is {node.coordinates}"
102 )
103 if not node.rotation == rot:
104 raise ValueError(f"{name} rotation does not match with function!")
106 if arc_length is not None:
107 if _np.abs(node.arc_length - arc_length) > _mpy.eps_pos:
108 raise ValueError(
109 f"Arc lengths don't match, got {node.arc_length} and {arc_length}"
110 )
112 # Flags if nodes are given
113 has_start_node = start_node is not None
114 has_end_node = end_node is not None
116 # Loop over local nodes.
117 arc_length = None
118 for i, xi in enumerate(self.nodes_create):
119 # Get the position and rotation at xi
120 pos, rot, arc_length_from_function = beam_function(xi)
121 if relative_twist is not None:
122 rot = rot * relative_twist
123 if set_nodal_arc_length:
124 arc_length = arc_length_from_function + nodal_arc_length_offset
126 # Check if the position and rotation match existing nodes
127 if i == 0 and has_start_node:
128 check_node(start_node, pos, rot, arc_length, "start_node")
129 self.nodes = [start_node]
130 elif (i == len(self.nodes_create) - 1) and has_end_node:
131 check_node(end_node, pos, rot, arc_length, "end_node")
133 # Create the node
134 if (i > 0 or not has_start_node) and (
135 i < len(self.nodes_create) - 1 or not has_end_node
136 ):
137 is_middle_node = 0 < i < len(self.nodes_create) - 1
138 self.nodes.append(
139 _NodeCosserat(
140 pos, rot, is_middle_node=is_middle_node, arc_length=arc_length
141 )
142 )
144 # Get a list with the created nodes.
145 if has_start_node:
146 created_nodes = self.nodes[1:]
147 else:
148 created_nodes = self.nodes.copy()
150 # Add the end node to the beam.
151 if has_end_node:
152 self.nodes.append(end_node)
154 # Return the created nodes.
155 return created_nodes
157 @classmethod
158 def get_coupling_dict(cls, coupling_dof_type):
159 """Return the dict to couple this beam to another beam."""
161 match coupling_dof_type:
162 case _mpy.coupling_dof.joint:
163 if cls.coupling_joint_dict is None:
164 raise ValueError(f"Joint coupling is not implemented for {cls}")
165 return cls.coupling_joint_dict
166 case _mpy.coupling_dof.fix:
167 if cls.coupling_fix_dict is None:
168 raise ValueError("Fix coupling is not implemented for {cls}")
169 return cls.coupling_fix_dict
170 case _:
171 raise ValueError(
172 f'Coupling_dof_type "{coupling_dof_type}" is not implemented!'
173 )
175 def flip(self):
176 """Reverse the nodes of this element.
178 This is usually used when reflected.
179 """
180 self.nodes = [self.nodes[-1 - i] for i in range(len(self.nodes))]
182 def _check_material(self):
183 """Check if the linked material is valid for this type of beam
184 element."""
185 for material_type in self.valid_material:
186 if isinstance(self.material, material_type):
187 break
188 else:
189 raise TypeError(
190 f"Beam of type {type(self)} can not have a material of type {type(self.material)}!"
191 )
193 def get_vtk(
194 self,
195 vtk_writer_beam,
196 vtk_writer_solid,
197 *,
198 beam_centerline_visualization_segments=1,
199 **kwargs,
200 ):
201 """Add the representation of this element to the VTK writer as a poly
202 line.
204 Args
205 ----
206 vtk_writer_beam:
207 VTK writer for the beams.
208 vtk_writer_solid:
209 VTK writer for solid elements, not used in this method.
210 beam_centerline_visualization_segments: int
211 Number of segments to be used for visualization of beam centerline between successive
212 nodes. Default is 1, which means a straight line is drawn between the beam nodes. For
213 Values greater than 1, a Hermite interpolation of the centerline is assumed for
214 visualization purposes.
215 """
217 n_nodes = len(self.nodes)
218 n_segments = n_nodes - 1
219 n_additional_points_per_segment = beam_centerline_visualization_segments - 1
220 # Number of points (in addition to the nodes) to be used for output
221 n_additional_points = n_segments * n_additional_points_per_segment
222 n_points = n_nodes + n_additional_points
224 # Dictionary with cell data.
225 cell_data = self.vtk_cell_data.copy()
226 cell_data["cross_section_radius"] = self.material.radius
228 # Dictionary with point data.
229 point_data = {}
230 point_data["node_value"] = _np.zeros(n_points)
231 point_data["base_vector_1"] = _np.zeros((n_points, 3))
232 point_data["base_vector_2"] = _np.zeros((n_points, 3))
233 point_data["base_vector_3"] = _np.zeros((n_points, 3))
235 coordinates = _np.zeros((n_points, 3))
236 nodal_rotation_matrices = [
237 node.rotation.get_rotation_matrix() for node in self.nodes
238 ]
240 for i_node, (node, rotation_matrix) in enumerate(
241 zip(self.nodes, nodal_rotation_matrices)
242 ):
243 coordinates[i_node, :] = node.coordinates
244 if node.is_middle_node:
245 point_data["node_value"][i_node] = 0.5
246 else:
247 point_data["node_value"][i_node] = 1.0
249 point_data["base_vector_1"][i_node] = rotation_matrix[:, 0]
250 point_data["base_vector_2"][i_node] = rotation_matrix[:, 1]
251 point_data["base_vector_3"][i_node] = rotation_matrix[:, 2]
253 # We can output the element partner index, which is a debug quantity to help find elements
254 # with matching middle nodes. This is usually an indicator for an issue with the mesh.
255 # TODO: Check if we really need this partner index output any more
256 element_partner_indices = list(
257 [
258 node.element_partner_index
259 for node in self.nodes
260 if node.element_partner_index is not None
261 ]
262 )
263 if len(element_partner_indices) > 1:
264 raise ValueError(
265 "More than one element partner indices are currently not supported in the output"
266 "functionality"
267 )
268 elif len(element_partner_indices) == 1:
269 cell_data["partner_index"] = element_partner_indices[0] + 1
271 # Check if we have everything we need to write output or if we need to calculate additional
272 # points for a smooth beam visualization.
273 if beam_centerline_visualization_segments == 1:
274 point_connectivity = _np.arange(n_nodes)
275 else:
276 # We need the centerline shape function matrices, so calculate them once and use for
277 # all segments that we need. Drop the first and last value, since they represent the
278 # nodes which we have already added above.
279 xi = _np.linspace(-1, 1, beam_centerline_visualization_segments + 1)[1:-1]
280 hermite_shape_functions_pos = _np.array(
281 [
282 0.25 * (2.0 + xi) * (1.0 - xi) ** 2,
283 0.25 * (2.0 - xi) * (1.0 + xi) ** 2,
284 ]
285 ).transpose()
286 hermite_shape_functions_tan = _np.array(
287 [
288 0.125 * (1.0 + xi) * (1.0 - xi) ** 2,
289 -0.125 * (1.0 - xi) * (1.0 + xi) ** 2,
290 ]
291 ).transpose()
293 point_connectivity = _np.zeros(n_points, dtype=int)
295 for i_segment in range(n_segments):
296 positions = _np.array(
297 [
298 self.nodes[i_node].coordinates
299 for i_node in [i_segment, i_segment + 1]
300 ]
301 )
302 tangents = _np.array(
303 [
304 nodal_rotation_matrices[i_node][:, 0]
305 for i_node in [i_segment, i_segment + 1]
306 ]
307 )
308 length_factor = _np.linalg.norm(positions[1] - positions[0])
309 interpolated_coordinates = _np.dot(
310 hermite_shape_functions_pos, positions
311 ) + length_factor * _np.dot(hermite_shape_functions_tan, tangents)
313 index_first_point = (
314 n_nodes + i_segment * n_additional_points_per_segment
315 )
316 index_last_point = (
317 n_nodes + (i_segment + 1) * n_additional_points_per_segment
318 )
320 coordinates[index_first_point:index_last_point] = (
321 interpolated_coordinates
322 )
323 point_connectivity[
324 i_segment * beam_centerline_visualization_segments
325 ] = i_segment
326 point_connectivity[
327 (i_segment + 1) * beam_centerline_visualization_segments
328 ] = i_segment + 1
329 point_connectivity[
330 i_segment * beam_centerline_visualization_segments + 1 : (
331 i_segment + 1
332 )
333 * beam_centerline_visualization_segments
334 ] = _np.arange(index_first_point, index_last_point)
336 # Get the point data sets and add everything to the output file.
337 _add_point_data_node_sets(
338 point_data, self.nodes, extra_points=n_additional_points
339 )
340 indices = vtk_writer_beam.add_points(coordinates, point_data=point_data)
341 vtk_writer_beam.add_cell(
342 _vtk.vtkPolyLine, indices[point_connectivity], cell_data=cell_data
343 )