Coverage for src/meshpy/mesh_creation_functions/beam_basic_geometry.py: 91%
122 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 has functions to create basic geometry items with meshpy."""
24import warnings as _warnings
26import numpy as _np
28from meshpy.core.conf import mpy as _mpy
29from meshpy.core.mesh import Mesh as _Mesh
30from meshpy.core.rotation import Rotation as _Rotation
31from meshpy.mesh_creation_functions.beam_generic import (
32 create_beam_mesh_function as _create_beam_mesh_function,
33)
34from meshpy.utils.nodes import get_single_node as _get_single_node
37def create_beam_mesh_line(mesh, beam_class, material, start_point, end_point, **kwargs):
38 """Generate a straight line of beam elements.
40 Args
41 ----
42 mesh: Mesh
43 Mesh that the line will be added to.
44 beam_class: Beam
45 Class of beam that will be used for this line.
46 material: Material
47 Material for this line.
48 start_point, end_point: _np.array, list
49 3D-coordinates for the start and end point of the line.
51 **kwargs (for all of them look into create_beam_mesh_function)
52 ----
53 n_el: int
54 Number of equally spaced beam elements along the line. Defaults to 1.
55 Mutually exclusive with l_el.
56 l_el: float
57 Desired length of beam elements. Mutually exclusive with n_el.
58 Be aware, that this length might not be achieved, if the elements are
59 warped after they are created.
60 start_node: Node, GeometrySet
61 Node to use as the first node for this line. Use this if the line
62 is connected to other lines (angles have to be the same, otherwise
63 connections should be used). If a geometry set is given, it can
64 contain one, and one node only.
65 add_sets: bool
66 If this is true the sets are added to the mesh and then displayed
67 in eventual VTK output, even if they are not used for a boundary
68 condition or coupling.
70 Return
71 ----
72 return_set: GeometryName
73 Set with the 'start' and 'end' node of the line. Also a 'line' set
74 with all nodes of the line.
75 """
77 # Get geometrical values for this line.
78 start_point = _np.asarray(start_point)
79 end_point = _np.asarray(end_point)
80 direction = end_point - start_point
81 line_length = _np.linalg.norm(direction)
82 t1 = direction / line_length
84 # Check if the z or y axis are larger projected onto the direction.
85 # The tolerance is used here to ensure that round-off changes in the last digits of
86 # the floating point values don't switch the case. This increases the robustness in
87 # testing.
88 if abs(_np.dot(t1, [0, 0, 1])) < abs(_np.dot(t1, [0, 1, 0])) - _mpy.eps_quaternion:
89 t2 = [0, 0, 1]
90 else:
91 t2 = [0, 1, 0]
92 rotation = _Rotation.from_basis(t1, t2)
94 def get_beam_geometry(parameter_a, parameter_b):
95 """Return a function for the position along the beams axis."""
97 def beam_function(xi):
98 """Return a point on the beams axis for a given parameter
99 coordinate xi."""
100 point_a = start_point + parameter_a * direction
101 point_b = start_point + parameter_b * direction
102 pos = 0.5 * (1 - xi) * point_a + 0.5 * (1 + xi) * point_b
103 arc_length = (
104 0.5 * (1 - xi) * parameter_a + 0.5 * (1 + xi) * parameter_b
105 ) * line_length
106 return (pos, rotation, arc_length)
108 return beam_function
110 # Create the beam in the mesh
111 return _create_beam_mesh_function(
112 mesh,
113 beam_class=beam_class,
114 material=material,
115 function_generator=get_beam_geometry,
116 interval=[0.0, 1.0],
117 interval_length=line_length,
118 **kwargs,
119 )
122def create_beam_mesh_arc_segment_via_rotation(
123 mesh, beam_class, material, center, axis_rotation, radius, angle, **kwargs
124):
125 """Generate a circular segment of beam elements.
127 The circular segment is defined via a rotation, specifying the "initial"
128 triad of the beam at the beginning of the arc.
130 This function exists for compatibility reasons with older MeshPy implementations.
131 The user is encouraged to use the newer implementation create_beam_mesh_arc_segment_via_axis
133 Args
134 ----
135 mesh: Mesh
136 Mesh that the arc segment will be added to.
137 beam_class: Beam
138 Class of beam that will be used for this line.
139 material: Material
140 Material for this segment.
141 center: _np.array, list
142 Center of the arc.
143 axis_rotation: Rotation
144 This rotation defines the spatial orientation of the arc.
145 The 3rd base vector of this rotation is the rotation axis of the arc
146 segment. The segment starts in the direction of the 1st basis vector
147 and the starting point is along the 2nd basis vector.
148 radius: float
149 The radius of the segment.
150 angle: float
151 The central angle of this segment in radians.
153 **kwargs (for all of them look into create_beam_mesh_function)
154 ----
155 n_el: int
156 Number of equally spaced beam elements along the line. Defaults to 1.
157 Mutually exclusive with l_el.
158 l_el: float
159 Desired length of beam elements. Mutually exclusive with n_el.
160 Be aware, that this length might not be achieved, if the elements are
161 warped after they are created.
163 Return
164 ----
165 return_set: GeometryName
166 Set with the 'start' and 'end' node of the line. Also a 'line' set
167 with all nodes of the line.
168 """
170 # Convert the input to the one for create_beam_mesh_arc_segment_via_axis
171 axis = axis_rotation * [0, 0, 1]
172 start_point = center + radius * (axis_rotation * [0, -1, 0])
173 return create_beam_mesh_arc_segment_via_axis(
174 mesh, beam_class, material, axis, center, start_point, angle, **kwargs
175 )
178def create_beam_mesh_arc_segment_via_axis(
179 mesh,
180 beam_class,
181 material,
182 axis,
183 axis_point,
184 start_point,
185 angle,
186 **kwargs,
187):
188 """Generate a circular segment of beam elements.
190 The arc is defined via a rotation axis, a point on the rotation axis a starting
191 point, as well as the angle of the arc segment.
193 Args
194 ----
195 mesh: Mesh
196 Mesh that the arc segment will be added to.
197 beam_class: Beam
198 Class of beam that will be used for this line.
199 material: Material
200 Material for this segment.
201 axis: _np.array, list
202 Rotation axis of the arc.
203 axis_point: _np.array, list
204 Point lying on the rotation axis. Does not have to be the center of the arc.
205 start_point: _np.array, list
206 Start point of the arc.
207 angle: float
208 The central angle of this segment in radians.
210 **kwargs (for all of them look into create_beam_mesh_function)
211 ----
212 n_el: int
213 Number of equally spaced beam elements along the line. Defaults to 1.
214 Mutually exclusive with l_el.
215 l_el: float
216 Desired length of beam elements. Mutually exclusive with n_el.
217 Be aware, that this length might not be achieved, if the elements are
218 warped after they are created.
220 Return
221 ----
222 return_set: GeometryName
223 Set with the 'start' and 'end' node of the line. Also a 'line' set
224 with all nodes of the line.
225 """
227 # The angle can not be negative with the current implementation.
228 if angle <= 0.0:
229 raise ValueError(
230 "The angle for a beam arc segment has to be a positive number!"
231 )
233 # Shortest distance from the given point to the axis of rotation gives
234 # the "center" of the arc
235 axis = _np.asarray(axis)
236 axis_point = _np.asarray(axis_point)
237 start_point = _np.asarray(start_point)
239 axis = axis / _np.linalg.norm(axis)
240 diff = start_point - axis_point
241 distance = diff - _np.dot(_np.dot(diff, axis), axis)
242 radius = _np.linalg.norm(distance)
243 center = start_point - distance
245 # Get the rotation at the start
246 # No need to check the start node here, as eventual rotation offsets in
247 # tangential direction will be covered by the create beam functionality.
248 tangent = _np.cross(axis, distance)
249 tangent /= _np.linalg.norm(tangent)
250 start_rotation = _Rotation.from_rotation_matrix(
251 _np.transpose(_np.array([tangent, -distance / radius, axis]))
252 )
254 def get_beam_geometry(alpha, beta):
255 """Return a function for the position and rotation along the beam
256 axis."""
258 def beam_function(xi):
259 """Return a point and the triad on the beams axis for a given
260 parameter coordinate xi."""
261 phi = 0.5 * (xi + 1) * (beta - alpha) + alpha
262 arc_rotation = _Rotation(axis, phi)
263 rot = arc_rotation * start_rotation
264 pos = center + arc_rotation * distance
265 return (pos, rot, phi * radius)
267 return beam_function
269 # Create the beam in the mesh
270 return _create_beam_mesh_function(
271 mesh,
272 beam_class=beam_class,
273 material=material,
274 function_generator=get_beam_geometry,
275 interval=[0.0, angle],
276 interval_length=angle * radius,
277 **kwargs,
278 )
281def create_beam_mesh_arc_segment_2d(
282 mesh, beam_class, material, center, radius, phi_start, phi_end, **kwargs
283):
284 """Generate a circular segment of beam elements in the x-y plane.
286 Args
287 ----
288 mesh: Mesh
289 Mesh that the arc segment will be added to.
290 beam_class: Beam
291 Class of beam that will be used for this line.
292 material: Material
293 Material for this segment.
294 center: _np.array, list
295 Center of the arc. If the z component is not 0, an error will be
296 thrown.
297 radius: float
298 The radius of the segment.
299 phi_start, phi_end: float
300 The start and end angles of the arc w.r.t the x-axis. If the start
301 angle is larger than the end angle the beam faces in counter-clockwise
302 direction, and if the start angle is smaller than the end angle, the
303 beam faces in clockwise direction.
305 **kwargs (for all of them look into create_beam_mesh_function)
306 ----
307 n_el: int
308 Number of equally spaced beam elements along the line. Defaults to 1.
309 Mutually exclusive with l_el.
310 l_el: float
311 Desired length of beam elements. Mutually exclusive with n_el.
312 Be aware, that this length might not be achieved, if the elements are
313 warped after they are created.
315 Return
316 ----
317 return_set: GeometryName
318 Set with the 'start' and 'end' node of the line. Also a 'line' set
319 with all nodes of the line.
320 """
322 # The center point has to be on the x-y plane.
323 if _np.abs(center[2]) > _mpy.eps_pos:
324 raise ValueError("The z-value of center has to be 0!")
326 # Check if the beam is in clockwise or counter clockwise direction.
327 angle = phi_end - phi_start
328 axis = _np.array([0, 0, 1])
329 start_point = center + radius * (_Rotation(axis, phi_start) * [1, 0, 0])
331 counter_clockwise = _np.sign(angle) == 1
332 if not counter_clockwise:
333 # If the beam is not in counter clockwise direction, we have to flip
334 # the rotation axis.
335 axis = -1.0 * axis
337 return create_beam_mesh_arc_segment_via_axis(
338 mesh,
339 beam_class,
340 material,
341 axis,
342 center,
343 start_point,
344 _np.abs(angle),
345 **kwargs,
346 )
349def create_beam_mesh_line_at_node(
350 mesh, beam_class, material, start_node, length, **kwargs
351):
352 """Generate a straight line at a given node. The tangent will be the same
353 as at that node.
355 Args
356 ----
357 mesh: Mesh
358 Mesh that the arc segment will be added to.
359 beam_class: Beam
360 Class of beam that will be used for this line.
361 material: Material
362 Material for this segment.
363 start_node: _np.array, list
364 Point where the arc will continue.
365 length: float
366 Length of the line.
368 **kwargs (for all of them look into create_beam_mesh_function)
369 ----
370 n_el: int
371 Number of equally spaced beam elements along the line. Defaults to 1.
372 Mutually exclusive with l_el.
373 l_el: float
374 Desired length of beam elements. Mutually exclusive with n_el.
375 Be aware, that this length might not be achieved, if the elements are
376 warped after they are created.
378 Return
379 ----
380 return_set: GeometryName
381 Set with the 'start' and 'end' node of the line. Also a 'line' set
382 with all nodes of the line.
383 """
385 if length < 0:
386 raise ValueError("Length has to be positive!")
388 # Create the line starting from the given node
389 start_node = _get_single_node(start_node)
390 tangent = start_node.rotation * [1, 0, 0]
391 start_position = start_node.coordinates
392 end_position = start_position + tangent * length
394 return create_beam_mesh_line(
395 mesh,
396 beam_class,
397 material,
398 start_position,
399 end_position,
400 start_node=start_node,
401 **kwargs,
402 )
405def create_beam_mesh_arc_at_node(
406 mesh, beam_class, material, start_node, arc_axis_normal, radius, angle, **kwargs
407):
408 """Generate a circular segment starting at a given node. The arc will be
409 tangent to the given node.
411 Args
412 ----
413 mesh: Mesh
414 Mesh that the arc segment will be added to.
415 beam_class: Beam
416 Class of beam that will be used for this line.
417 material: Material
418 Material for this segment.
419 start_node: _np.array, list
420 Point where the arc will continue.
421 arc_axis_normal: 3d-vector
422 Rotation axis for the created arc.
423 radius: float
424 The radius of the arc segment.
425 angle: float
426 Angle of the arc. If the angle is negative, the arc will point in the
427 opposite direction, i.e., as if the arc_axis_normal would change sign.
429 **kwargs (for all of them look into create_beam_mesh_function)
430 ----
431 n_el: int
432 Number of equally spaced beam elements along the line. Defaults to 1.
433 Mutually exclusive with l_el.
434 l_el: float
435 Desired length of beam elements. Mutually exclusive with n_el.
436 Be aware, that this length might not be achieved, if the elements are
437 warped after they are created.
439 Return
440 ----
441 return_set: GeometryName
442 Set with the 'start' and 'end' node of the line. Also a 'line' set
443 with all nodes of the line.
444 """
446 # If the angle is negative, the normal is switched
447 arc_axis_normal = _np.asarray(arc_axis_normal)
448 if angle < 0:
449 arc_axis_normal = -1.0 * arc_axis_normal
451 # The normal has to be perpendicular to the start point tangent
452 start_node = _get_single_node(start_node)
453 tangent = start_node.rotation * [1, 0, 0]
454 if _np.abs(_np.dot(tangent, arc_axis_normal)) > _mpy.eps_pos:
455 raise ValueError(
456 "The normal has to be perpendicular to the tangent in the start node!"
457 )
459 # Get the center of the arc
460 center_direction = _np.cross(tangent, arc_axis_normal)
461 center_direction *= 1.0 / _np.linalg.norm(center_direction)
462 center = start_node.coordinates - center_direction * radius
464 return create_beam_mesh_arc_segment_via_axis(
465 mesh,
466 beam_class,
467 material,
468 arc_axis_normal,
469 center,
470 start_node.coordinates,
471 _np.abs(angle),
472 start_node=start_node,
473 **kwargs,
474 )
477def create_beam_mesh_helix(
478 mesh,
479 beam_class,
480 material,
481 axis_vector,
482 axis_point,
483 start_point,
484 *,
485 helix_angle=None,
486 height_helix=None,
487 turns=None,
488 warning_straight_line=True,
489 **kwargs,
490):
491 """Generate a helical segment starting at a given start point around a
492 predefined axis (defined by axis_vector and axis_point). The helical
493 segment is defined by a start_point and exactly two of the basic helical
494 quantities [helix_angle, height_helix, turns].
496 Args
497 ----
498 mesh: Mesh
499 Mesh that the helical segment will be added to.
500 beam_class: Beam
501 Class of beam that will be used for this line.
502 material: Material
503 Material for this segment.
504 axis_vector: _np.array, list
505 Vector for the orientation of the helical center axis.
506 axis_point: _np.array, list
507 Point lying on the helical center axis. Does not need to align with
508 bottom plane of helix.
509 start_point: _np.array, list
510 Start point of the helix. Defines the radius.
511 helix_angle: float
512 Angle of the helix (synonyms in literature: twist angle or pitch
513 angle).
514 height_helix: float
515 Height of helix.
516 turns: float
517 Number of turns.
518 warning_straight_line: bool
519 Warn if radius of helix is zero or helix angle is 90 degrees and
520 simple line is returned.
522 **kwargs (for all of them look into create_beam_mesh_function)
523 ----
524 n_el: int
525 Number of equally spaced beam elements along the line. Defaults to 1.
526 Mutually exclusive with l_el.
527 l_el: float
528 Desired length of beam elements. Mutually exclusive with n_el.
529 Be aware, that this length might not be achieved, if the elements are
530 warped after they are created.
532 Return
533 ----
534 return_set: GeometryName
535 Set with the 'start' and 'end' node of the line. Also a 'line' set
536 with all nodes of the line.
537 """
539 if [helix_angle, height_helix, turns].count(None) != 1:
540 raise ValueError(
541 "Exactly two arguments of [helix_angle, height_helix, turns]"
542 " must be provided!"
543 )
545 if helix_angle is not None and _np.isclose(_np.sin(helix_angle), 0.0):
546 raise ValueError(
547 "Helix angle of helix is 0 degrees! "
548 + "Change angle for feasible helix geometry!"
549 )
551 if height_helix is not None and _np.isclose(height_helix, 0.0):
552 raise ValueError(
553 "Height of helix is 0! Change height for feasible helix geometry!"
554 )
556 # determine radius of helix
557 axis_vector = _np.asarray(axis_vector)
558 axis_point = _np.asarray(axis_point)
559 start_point = _np.asarray(start_point)
561 axis_vector = axis_vector / _np.linalg.norm(axis_vector)
562 origin = axis_point + _np.dot(
563 _np.dot(start_point - axis_point, axis_vector), axis_vector
564 )
565 start_point_origin_vec = start_point - origin
566 radius = _np.linalg.norm(start_point_origin_vec)
568 # create temporary mesh to not alter original mesh
569 mesh_temp = _Mesh()
571 # return line if radius of helix is 0, helix angle is pi/2 or turns is 0
572 if (
573 _np.isclose(radius, 0)
574 or (helix_angle is not None and _np.isclose(_np.cos(helix_angle), 0.0))
575 or (turns is not None and _np.isclose(turns, 0.0))
576 ):
577 if height_helix is None:
578 raise ValueError(
579 "Radius of helix is 0, helix angle is 90 degrees or turns is 0! "
580 + "Fallback to simple line geometry but height cannot be "
581 + "determined based on helix angle and turns! Either switch one "
582 + "helix parameter to height of helix or change radius!"
583 )
585 if warning_straight_line:
586 _warnings.warn(
587 "Radius of helix is 0, helix angle is 90 degrees or turns is 0! "
588 + "Simple line geometry is returned!"
589 )
591 if helix_angle is not None and height_helix is not None:
592 end_point = start_point + height_helix * axis_vector * _np.sign(
593 _np.sin(helix_angle)
594 )
595 elif height_helix is not None and turns is not None:
596 end_point = start_point + height_helix * axis_vector
598 line_sets = create_beam_mesh_line(
599 mesh_temp,
600 beam_class,
601 material,
602 start_point=start_point,
603 end_point=end_point,
604 **kwargs,
605 )
607 # add line to mesh
608 mesh.add_mesh(mesh_temp)
610 return line_sets
612 # generate simple helix
613 if helix_angle and height_helix:
614 end_point = _np.array(
615 [
616 radius,
617 _np.sign(_np.sin(helix_angle)) * height_helix / _np.tan(helix_angle),
618 _np.sign(_np.sin(helix_angle)) * height_helix,
619 ]
620 )
621 elif helix_angle and turns:
622 end_point = _np.array(
623 [
624 radius,
625 _np.sign(_np.cos(helix_angle)) * 2 * _np.pi * radius * turns,
626 _np.sign(_np.cos(helix_angle))
627 * 2
628 * _np.pi
629 * radius
630 * _np.abs(turns)
631 * _np.tan(helix_angle),
632 ]
633 )
634 elif height_helix and turns:
635 end_point = _np.array(
636 [
637 radius,
638 2 * _np.pi * radius * turns,
639 height_helix,
640 ]
641 )
643 helix_sets = create_beam_mesh_line(
644 mesh_temp,
645 beam_class,
646 material,
647 start_point=[radius, 0, 0],
648 end_point=end_point,
649 **kwargs,
650 )
652 mesh_temp.wrap_around_cylinder()
654 # rotate and translate simple helix to align with necessary axis and starting point
655 mesh_temp.rotate(
656 _Rotation.from_basis(start_point_origin_vec, axis_vector)
657 * _Rotation([1, 0, 0], -_np.pi * 0.5)
658 )
659 mesh_temp.translate(-mesh_temp.nodes[0].coordinates + start_point)
661 # add helix to mesh
662 mesh.add_mesh(mesh_temp)
664 return helix_sets