Coverage for src/meshpy/mesh_creation_functions/beam_generic.py: 94%
123 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"""Generic function used to create all beams within meshpy."""
24from typing import Callable as _Callable
25from typing import Dict as _Dict
26from typing import List as _List
27from typing import Optional as _Optional
28from typing import Tuple as _Tuple
29from typing import Type as _Type
30from typing import Union as _Union
32import numpy as _np
34from meshpy.core.conf import mpy as _mpy
35from meshpy.core.element_beam import Beam as _Beam
36from meshpy.core.geometry_set import GeometryName as _GeometryName
37from meshpy.core.geometry_set import GeometrySet as _GeometrySet
38from meshpy.core.material import MaterialBeam as _MaterialBeam
39from meshpy.core.mesh import Mesh as _Mesh
40from meshpy.core.node import NodeCosserat as _NodeCosserat
41from meshpy.utils.nodes import get_single_node as _get_single_node
44def create_beam_mesh_function(
45 mesh: _Mesh,
46 *,
47 beam_class: _Type[_Beam],
48 material: _MaterialBeam,
49 function_generator: _Callable,
50 interval: _Tuple[float, float],
51 n_el: _Optional[int] = None,
52 l_el: _Optional[float] = None,
53 interval_length: _Optional[float] = None,
54 set_nodal_arc_length: bool = False,
55 nodal_arc_length_offset: _Optional[float] = None,
56 node_positions_of_elements: _Optional[_List[float]] = None,
57 start_node: _Optional[_Union[_NodeCosserat, _GeometrySet]] = None,
58 end_node: _Optional[_Union[_NodeCosserat, _GeometrySet]] = None,
59 close_beam: bool = False,
60 vtk_cell_data: _Optional[_Dict[str, _Tuple]] = None,
61) -> _GeometryName:
62 """Generic beam creation function.
64 Remark for given start and/or end nodes:
65 If the rotation does not match, but the tangent vector is the same,
66 the created beams triads are rotated so the physical problem stays
67 the same (for axi-symmetric beam cross-sections) but the nodes can
68 be reused.
70 Args:
71 mesh:
72 Mesh that the created beam(s) should be added to.
73 beam_class:
74 Class of beam that will be used for this line.
75 material:
76 Material for this line.
77 function_generator:
78 The function_generator has to take two variables, point_a and
79 point_b (both within the interval) and return a function(xi) that
80 calculates the position and rotation along the beam, where
81 point_a -> xi = -1 and point_b -> xi = 1.
83 Usually, the Jacobian of the returned position field should be a unit
84 vector. Otherwise, the nodes may be spaced in an undesired way. All
85 standard mesh creation functions in MeshPy fulfill this property.
86 interval:
87 Start and end values for interval that will be used to create the
88 beam.
89 n_el:
90 Number of equally spaced beam elements along the line. Defaults to 1.
91 Mutually exclusive with l_el
92 l_el:
93 Desired length of beam elements. This requires the option interval_length
94 to be set. Mutually exclusive with n_el. Be aware, that this length
95 might not be achieved, if the elements are warped after they are
96 created.
97 interval_length:
98 Total length of the interval. Is required when the option l_el is given.
99 set_nodal_arc_length:
100 Flag if the arc length along the beam filament is set in the created
101 nodes. It is ensured that the arc length is consistent with possible
102 given start/end nodes.
103 nodal_arc_length_offset:
104 Offset of the stored nodal arc length w.r.t. to the one generated by
105 the function. Defaults to 0, the arc length set in the start node, or
106 the arc length in the end node minus total length (such that the arc
107 length at the end node matches).
108 node_positions_of_elements:
109 A list of normalized positions (within [0,1] and in ascending order)
110 that define the boundaries of beam elements along the created curve.
111 The given values will be mapped to the actual `interval` given as an
112 argument to this function. These values specify where elements start
113 and end, additional internal nodes (such as midpoints in higher-order
114 elements) may be placed automatically.
115 start_node:
116 Node to use as the first node for this line. Use this if the line
117 is connected to other lines (angles have to be the same, otherwise
118 connections should be used). If a geometry set is given, it can
119 contain one, and one node only.
120 end_node:
121 If this is a Node or GeometrySet, the last node of the created beam
122 is set to that node.
123 If it is True the created beam is closed within itself.
124 close_beam:
125 If it is True the created beam is closed within itself (mutually
126 exclusive with end_node).
127 vtk_cell_data:
128 With this argument, a vtk cell data can be set for the elements
129 created within this function. This can be used to check which
130 elements are created by which function.
132 Returns:
133 Geometry sets with the 'start' and 'end' node of the curve. Also a 'line' set
134 with all nodes of the curve.
135 """
137 # Check for mutually exclusive parameters
138 n_given_arguments = sum(
139 1
140 for argument in [n_el, l_el, node_positions_of_elements]
141 if argument is not None
142 )
143 if n_given_arguments == 0:
144 # No arguments were given, use a single element per default
145 n_el = 1
146 elif n_given_arguments > 1:
147 raise ValueError(
148 'The arguments "n_el", "l_el" and "node_positions_of_elements" are mutually exclusive'
149 )
151 if close_beam and end_node is not None:
152 raise ValueError(
153 'The arguments "close_beam" and "end_node" are mutually exclusive'
154 )
156 if set_nodal_arc_length:
157 if close_beam:
158 raise ValueError(
159 "The flags 'set_nodal_arc_length' and 'close_beam' are mutually exclusive."
160 )
161 elif nodal_arc_length_offset is not None:
162 raise ValueError(
163 'Providing the argument "nodal_arc_length_offset" without setting '
164 '"set_nodal_arc_length" to True does not make sense.'
165 )
167 # Cases where we have equally spaced elements
168 if n_el is not None or l_el is not None:
169 if l_el is not None:
170 # Calculate the number of elements in case a desired element length is provided
171 if interval_length is None:
172 raise ValueError(
173 'The parameter "l_el" requires "interval_length" to be set.'
174 )
175 n_el = max([1, round(interval_length / l_el)])
176 elif n_el is None:
177 raise ValueError("n_el should not be None at this point")
179 node_positions_of_elements = [i_node / n_el for i_node in range(n_el + 1)]
180 # A list for the element node positions was provided
181 elif node_positions_of_elements is not None:
182 # Check that the given positions are in ascending order and start with 0 and end with 1
183 for index, value, name in zip([0, -1], [0, 1], ["First", "Last"]):
184 if not _np.isclose(
185 value,
186 node_positions_of_elements[index],
187 atol=1e-12,
188 rtol=0.0,
189 ):
190 raise ValueError(
191 f"{name} entry of node_positions_of_elements must be {value}, got {node_positions_of_elements[index]}"
192 )
193 if not all(
194 x < y
195 for x, y in zip(node_positions_of_elements, node_positions_of_elements[1:])
196 ):
197 raise ValueError(
198 f"The given node_positions_of_elements must be in ascending order. Got {node_positions_of_elements}"
199 )
201 # Get the scale the node positions to the interval coordinates
202 interval_node_positions_of_elements = interval[0] + (
203 interval[1] - interval[0]
204 ) * _np.asarray(node_positions_of_elements)
206 # We need to make sure we have the number of elements for the case a given end node
207 n_el = len(interval_node_positions_of_elements) - 1
209 # Make sure the material is in the mesh.
210 mesh.add_material(material)
212 # List with nodes and elements that will be added in the creation of
213 # this beam.
214 elements = []
215 nodes = []
217 def check_given_node(node):
218 """Check that the given node is already in the mesh."""
219 if node not in mesh.nodes:
220 raise ValueError("The given node is not in the current mesh")
222 def get_relative_twist(rotation_node, rotation_function, name):
223 """Check if the rotation at a node and the one returned by the function
224 match.
226 If not, check if the first basis vector of the triads is the
227 same. If that is the case, a simple relative twist can be
228 applied to ensure that the triad field is continuous. This
229 relative twist can lead to issues if the beam cross-section is
230 not double symmetric.
231 """
233 if rotation_node == rotation_function:
234 return None
235 elif not _mpy.allow_beam_rotation:
236 # The settings do not allow for a rotation of the beam
237 raise ValueError(
238 f"Given rotation of the {name} node does not match with given rotation function!"
239 )
240 else:
241 # Evaluate the relative rotation
242 # First check if the first basis vector is the same
243 relative_basis_1 = rotation_node.inv() * rotation_function * [1, 0, 0]
244 if _np.linalg.norm(relative_basis_1 - [1, 0, 0]) < _mpy.eps_quaternion:
245 # Calculate the relative rotation
246 return rotation_function.inv() * rotation_node
247 else:
248 raise ValueError(
249 f"The tangent of the {name} node does not match with the given function!"
250 )
252 # Position and rotation at the start and end of the interval
253 function_over_whole_interval = function_generator(*interval)
254 relative_twist_start = None
255 relative_twist_end = None
257 # If a start node is given, set this as the first node for this beam.
258 if start_node is not None:
259 start_node = _get_single_node(start_node)
260 nodes = [start_node]
261 check_given_node(start_node)
262 _, start_rotation, _ = function_over_whole_interval(-1.0)
263 relative_twist_start = get_relative_twist(
264 start_node.rotation, start_rotation, "start"
265 )
267 # If an end node is given, check what behavior is wanted.
268 if end_node is not None:
269 end_node = _get_single_node(end_node)
270 check_given_node(end_node)
271 _, end_rotation, _ = function_over_whole_interval(1.0)
272 relative_twist_end = get_relative_twist(end_node.rotation, end_rotation, "end")
274 # Get the start value for the arc length functionality
275 if set_nodal_arc_length:
276 if nodal_arc_length_offset is not None:
277 # Let's use the given value, if it does not match with possible given
278 # start or end nodes, the check in the create beam function will detect
279 # that.
280 pass
281 elif start_node is not None and start_node.arc_length is not None:
282 nodal_arc_length_offset = start_node.arc_length
283 elif end_node is not None and end_node.arc_length is not None:
284 nodal_arc_length_offset = end_node.arc_length - interval_length
285 else:
286 # Default value
287 nodal_arc_length_offset = 0.0
289 # Check if a relative twist has to be applied
290 if relative_twist_start is not None and relative_twist_end is not None:
291 if relative_twist_start == relative_twist_end:
292 relative_twist = relative_twist_start
293 else:
294 raise ValueError(
295 "The relative twist required for the start and end node do not match"
296 )
297 elif relative_twist_start is not None:
298 relative_twist = relative_twist_start
299 elif relative_twist_end is not None:
300 relative_twist = relative_twist_end
301 else:
302 relative_twist = None
304 # Create the beams.
305 for i_el in range(n_el):
306 # If the beam is closed with itself, set the end node to be the
307 # first node of the beam. This is done when the second element is
308 # created, as the first node already exists here.
309 if i_el == 1 and close_beam:
310 end_node = nodes[0]
312 # Get the function to create this beam element.
313 function = function_generator(
314 interval_node_positions_of_elements[i_el],
315 interval_node_positions_of_elements[i_el + 1],
316 )
318 # Set the start node for the created beam.
319 if start_node is not None or i_el > 0:
320 first_node = nodes[-1]
321 else:
322 first_node = None
324 # If an end node is given, set this one for the last element.
325 if end_node is not None and i_el == n_el - 1:
326 last_node = end_node
327 else:
328 last_node = None
330 element = beam_class(material=material)
331 elements.append(element)
332 nodes.extend(
333 element.create_beam(
334 function,
335 start_node=first_node,
336 end_node=last_node,
337 relative_twist=relative_twist,
338 set_nodal_arc_length=set_nodal_arc_length,
339 nodal_arc_length_offset=nodal_arc_length_offset,
340 )
341 )
343 # Set vtk cell data on created elements.
344 if vtk_cell_data is not None:
345 for data_name, data_value in vtk_cell_data.items():
346 for element in elements:
347 if data_name in element.vtk_cell_data.keys():
348 raise KeyError(
349 'The cell data "{}" already exists!'.format(data_name)
350 )
351 element.vtk_cell_data[data_name] = data_value
353 # Add items to the mesh
354 mesh.elements.extend(elements)
355 if start_node is None:
356 mesh.nodes.extend(nodes)
357 else:
358 mesh.nodes.extend(nodes[1:])
360 # Set the last node of the beam.
361 if end_node is None:
362 end_node = nodes[-1]
364 # Set the nodes that are at the beginning and end of line (for search
365 # of overlapping points)
366 nodes[0].is_end_node = True
367 end_node.is_end_node = True
369 # Create geometry sets that will be returned.
370 return_set = _GeometryName()
371 return_set["start"] = _GeometrySet(nodes[0])
372 return_set["end"] = _GeometrySet(end_node)
373 return_set["line"] = _GeometrySet(elements)
375 return return_set