Merge pull request #60 from davidrudlstorfer/add_helical_base_geometry

Add helical beam basic geometry function

CONTRIBUTORS

@@ -8,3 +8,4 @@ Dao Viet Anh
 Nora Hagmeyer (@NoraHagmeyer)
 Matthias Mayr (@mayrmt)
 Gabriela Loera (@eulovi)
+David Rudlstorfer (@davidrudlstorfer)

meshpy/mesh_creation_functions/__init__.py

@@ -41,6 +41,7 @@
     create_beam_mesh_arc_segment_2d,
     create_beam_mesh_line_at_node,
     create_beam_mesh_arc_at_node,
+    create_beam_mesh_helix,
 )
 
 # Parametric curve.
@@ -82,6 +83,7 @@
     "create_beam_mesh_arc_segment_2d",
     "create_beam_mesh_line_at_node",
     "create_beam_mesh_arc_at_node",
+    "create_beam_mesh_helix",
     # Parametric curve.
     "create_beam_mesh_curve",
     # Honeycomb.

meshpy/mesh_creation_functions/beam_basic_geometry.py

@@ -34,6 +34,7 @@
 
 # Python packages.
 import numpy as np
+import warnings
 
 # Meshpy modules.
 from ..conf import mpy
@@ -91,7 +92,10 @@ def create_beam_mesh_line(
     t1 = direction / line_length
 
     # Check if the z or y axis are larger projected onto the direction.
-    if abs(np.dot(t1, [0, 0, 1])) < abs(np.dot(t1, [0, 1, 0])):
+    # The tolerance is used here to ensure that round-off changes in the last digits of
+    # the floating point values don't switch the case. This increases the robustness in
+    # testing.
+    if abs(np.dot(t1, [0, 0, 1])) < abs(np.dot(t1, [0, 1, 0])) - mpy.eps_quaternion:
         t2 = [0, 0, 1]
     else:
         t2 = [0, 1, 0]
@@ -476,3 +480,150 @@ def create_beam_mesh_arc_at_node(
         start_node=start_node,
         **kwargs,
     )
+
+
+def create_beam_mesh_helix(
+    mesh,
+    beam_object,
+    material,
+    axis,
+    axis_point,
+    start_point,
+    twist_angle,
+    *,
+    height_helix=None,
+    turns=None,
+    warning_straight_line=True,
+    **kwargs
+):
+    """
+    Generate a helical segment starting at a given start point around a
+    predefined axis.
+
+    Args
+    ----
+    mesh: Mesh
+        Mesh that the helical segment will be added to.
+    beam_object: Beam
+        Class of beam that will be used for this line.
+    material: Material
+        Material for this segment.
+    axis: np.array, list
+        Rotation axis of the helix.
+    axis_point: np.array, list
+        Point lying on the rotation axis. Does not need to align with
+        bottom plane of helix.
+    start_point: np.array, list
+        Start point of the helix. Defines the radius.
+    twist_angle: float
+        Twist angle / Pitch of helix between 0 and 2pi.
+    height_helix: float
+        Height of helix. Mutually exclusive with number of turns.
+    turns: float
+        Number of turns. Mutually exclusive with height of helix.
+    warning_straight_line: bool
+        Warn if radius of helix is zero or twist angle is 90 degrees and
+        simple line is returned.
+
+    **kwargs (for all of them look into create_beam_mesh_function)
+    ----
+    n_el: int
+        Number of equally spaced beam elements along the line. Defaults to 1.
+        Mutually exclusive with l_el.
+    l_el: float
+        Desired length of beam elements. Mutually exclusive with n_el.
+        Be aware, that this length might not be achieved, if the elements are
+        warped after they are created.
+
+    Return
+    ----
+    return_set: GeometryName
+        Set with the 'start' and 'end' node of the line. Also a 'line' set
+        with all nodes of the line.
+    """
+
+    if height_helix is None and turns is None:
+        raise ValueError("Either provide height_helix or turns!")
+    elif height_helix is not None and turns is not None:
+        raise ValueError("Only provide height_helix OR turns!")
+
+    if np.isclose(np.sin(twist_angle), 0.0):
+        raise ValueError(
+            "Twist angle of helix is 0 degrees! "
+            + "Change angle for feasible helix geometry!"
+        )
+
+    # determine radius of helix
+    axis = np.array(axis) / np.linalg.norm(np.array(axis))
+    origin = axis_point + np.dot(
+        np.dot((np.array(start_point) - np.array(axis_point)), axis), axis
+    )
+    start_point_origin_vec = start_point - origin
+    radius = np.linalg.norm(start_point_origin_vec)
+
+    # return line if radius of helix is 0 or twist angle is np.pi/2
+    if np.isclose(radius, 0) or np.isclose(np.cos(twist_angle), 0.0):
+        if turns:
+            raise ValueError(
+                "Radius of helix is 0 or twist angle is 90 degrees! "
+                + "Height of helix can not be determined through turns! "
+                + "Either switch to height of helix or change radius!"
+            )
+
+        if warning_straight_line:
+            warnings.warn(
+                "Radius of helix is 0 or twist angle is 90 degrees! "
+                + "Simple line geometry is returned!"
+            )
+
+        return create_beam_mesh_line(
+            mesh,
+            beam_object,
+            material,
+            start_point=start_point,
+            end_point=start_point + height_helix * axis * np.sign(np.sin(twist_angle)),
+            **kwargs,
+        )
+
+    # generate simple helix
+    if height_helix:
+        end_point = np.array(
+            [
+                radius,
+                np.sign(np.sin(twist_angle)) * height_helix / np.tan(twist_angle),
+                np.sign(np.sin(twist_angle)) * height_helix,
+            ]
+        )
+    elif turns:
+        end_point = np.array(
+            [
+                radius,
+                np.sign(np.cos(twist_angle)) * 2 * np.pi * radius * turns,
+                np.sign(np.cos(twist_angle))
+                * 2
+                * np.pi
+                * radius
+                * turns
+                * np.tan(twist_angle),
+            ]
+        )
+
+    helix = create_beam_mesh_line(
+        mesh,
+        beam_object,
+        material,
+        start_point=[radius, 0, 0],
+        end_point=end_point,
+        **kwargs,
+    )
+
+    mesh.wrap_around_cylinder()
+
+    # rotate and translate simple helix to align with neccessary axis and starting point
+    mesh.rotate(
+        Rotation.from_basis(start_point_origin_vec, axis)
+        * Rotation([1, 0, 0], -np.pi * 0.5)
+    )
+    mesh.translate(-mesh.nodes[0].coordinates + start_point)
+
+    return helix

tests/reference-files/test_mesh_creation_functions_helix_no_rotation_reference.dat

@@ -0,0 +1,35 @@
+// -----------------------------------------------------------------------------
+// This input file was created with MeshPy.
+// Copyright (c) 2018-2024
+//     Ivo Steinbrecher
+//     Institute for Mathematics and Computer-Based Simulation
+//     Universitaet der Bundeswehr Muenchen
+//     https://www.unibw.de/imcs-en
+// -----------------------------------------------------------------------------
+-----------------------------------------------------------------------MATERIALS
+MAT 1 MAT_BeamReissnerElastHyper YOUNG 100000.0 POISSONRATIO 0.0 DENS 0.0 CROSSAREA 0.7853981633974483 SHEARCORR 1.0 MOMINPOL 0.09817477042468103 MOMIN2 0.04908738521234052 MOMIN3 0.04908738521234052
+-------------------------------------------------------------DNODE-NODE TOPOLOGY
+NODE 7 DNODE 1
+NODE 1 DNODE 2
+-------------------------------------------------------------DLINE-NODE TOPOLOGY
+NODE 1 DLINE 1
+NODE 2 DLINE 1
+NODE 3 DLINE 1
+NODE 4 DLINE 1
+NODE 5 DLINE 1
+NODE 6 DLINE 1
+NODE 7 DLINE 1
+---------------------------------------------------------------------NODE COORDS
+NODE 1 COORD 2 0 0
+NODE 2 COORD 1.34482448817 1.48035370639 1.66666666667
+NODE 3 COORD -0.191447096029 1.9908159155 3.33333333333
+NODE 4 COORD -1.60228723109 1.19694428821 5
+NODE 5 COORD -1.96334800942 -0.381135925751 6.66666666667
+NODE 6 COORD -1.03807125077 -1.70950521448 8.33333333333
+NODE 7 COORD 0.567324370926 -1.91784854933 10
+--------------------------------------------------------------STRUCTURE ELEMENTS
+1 BEAM3R HERM2LINE3 1 3 2 MAT 1 TRIADS -0.613943125569 -1.48218982027 0.613943125569 0.825845001377 -1.75132059139 1.82233126689 0.0412003090647 -1.71867210284 1.24441771648
+2 BEAM3R HERM2LINE3 3 5 4 MAT 1 TRIADS 0.825845001377 -1.75132059139 1.82233126689 -1.76693362712 0.540443031687 -1.63156168219 -1.62734522345 1.40867035905 -2.14678720617
+3 BEAM3R HERM2LINE3 5 7 6 MAT 1 TRIADS -1.76693362712 0.540443031687 -1.63156168219 -1.36158461459 -0.809718268376 -0.390228435716 -1.6574070681 -0.199576520935 -1.03084186572
+------------------------------------------------------------------FLUID ELEMENTS
+-----------------------------------------------------------------------------END

tests/reference-files/test_mesh_creation_functions_helix_radius_zero_reference.dat

@@ -0,0 +1,40 @@
+// -----------------------------------------------------------------------------
+// This input file was created with MeshPy.
+// Copyright (c) 2018-2024
+//     Ivo Steinbrecher
+//     Institute for Mathematics and Computer-Based Simulation
+//     Universitaet der Bundeswehr Muenchen
+//     https://www.unibw.de/imcs-en
+// -----------------------------------------------------------------------------
+-----------------------------------------------------------------------MATERIALS
+MAT 1 MAT_BeamReissnerElastHyper YOUNG 100000.0 POISSONRATIO 0.0 DENS 0.0 CROSSAREA 0.7853981633974483 SHEARCORR 1.0 MOMINPOL 0.09817477042468103 MOMIN2 0.04908738521234052 MOMIN3 0.04908738521234052
+-------------------------------------------------------------DNODE-NODE TOPOLOGY
+NODE 9 DNODE 1
+NODE 1 DNODE 2
+-------------------------------------------------------------DLINE-NODE TOPOLOGY
+NODE 1 DLINE 1
+NODE 2 DLINE 1
+NODE 3 DLINE 1
+NODE 4 DLINE 1
+NODE 5 DLINE 1
+NODE 6 DLINE 1
+NODE 7 DLINE 1
+NODE 8 DLINE 1
+NODE 9 DLINE 1
+---------------------------------------------------------------------NODE COORDS
+NODE 1 COORD 1 1 1
+NODE 2 COORD 6.7735026919 6.7735026919 6.7735026919
+NODE 3 COORD 12.5470053838 12.5470053838 12.5470053838
+NODE 4 COORD 18.3205080757 18.3205080757 18.3205080757
+NODE 5 COORD 24.0940107676 24.0940107676 24.0940107676
+NODE 6 COORD 29.8675134595 29.8675134595 29.8675134595
+NODE 7 COORD 35.6410161514 35.6410161514 35.6410161514
+NODE 8 COORD 41.4145188433 41.4145188433 41.4145188433
+NODE 9 COORD 47.1880215352 47.1880215352 47.1880215352
+--------------------------------------------------------------STRUCTURE ELEMENTS
+1 BEAM3R HERM2LINE3 1 3 2 MAT 1 TRIADS -0.241535826833 -0.75993556603 0.58311906894 -0.241535826833 -0.75993556603 0.58311906894 -0.241535826833 -0.75993556603 0.58311906894
+2 BEAM3R HERM2LINE3 3 5 4 MAT 1 TRIADS -0.241535826833 -0.75993556603 0.58311906894 -0.241535826833 -0.75993556603 0.58311906894 -0.241535826833 -0.75993556603 0.58311906894
+3 BEAM3R HERM2LINE3 5 7 6 MAT 1 TRIADS -0.241535826833 -0.75993556603 0.58311906894 -0.241535826833 -0.75993556603 0.58311906894 -0.241535826833 -0.75993556603 0.58311906894
+4 BEAM3R HERM2LINE3 7 9 8 MAT 1 TRIADS -0.241535826833 -0.75993556603 0.58311906894 -0.241535826833 -0.75993556603 0.58311906894 -0.241535826833 -0.75993556603 0.58311906894
+------------------------------------------------------------------FLUID ELEMENTS
+-----------------------------------------------------------------------------END

tests/reference-files/test_mesh_creation_functions_helix_rotation_offset_reference.dat

@@ -0,0 +1,40 @@
+// -----------------------------------------------------------------------------
+// This input file was created with MeshPy.
+// Copyright (c) 2018-2024
+//     Ivo Steinbrecher
+//     Institute for Mathematics and Computer-Based Simulation
+//     Universitaet der Bundeswehr Muenchen
+//     https://www.unibw.de/imcs-en
+// -----------------------------------------------------------------------------
+-----------------------------------------------------------------------MATERIALS
+MAT 1 MAT_BeamReissnerElastHyper YOUNG 100000.0 POISSONRATIO 0.0 DENS 0.0 CROSSAREA 0.7853981633974483 SHEARCORR 1.0 MOMINPOL 0.09817477042468103 MOMIN2 0.04908738521234052 MOMIN3 0.04908738521234052
+-------------------------------------------------------------DNODE-NODE TOPOLOGY
+NODE 9 DNODE 1
+NODE 1 DNODE 2
+-------------------------------------------------------------DLINE-NODE TOPOLOGY
+NODE 1 DLINE 1
+NODE 2 DLINE 1
+NODE 3 DLINE 1
+NODE 4 DLINE 1
+NODE 5 DLINE 1
+NODE 6 DLINE 1
+NODE 7 DLINE 1
+NODE 8 DLINE 1
+NODE 9 DLINE 1
+---------------------------------------------------------------------NODE COORDS
+NODE 1 COORD 3 0 0
+NODE 2 COORD 2.9899942382 2.42676917578 -0.25169990452
+NODE 3 COORD 2.05197680159 4.33763484933 0.940515367999
+NODE 4 COORD 1.40034341356 4.86248321003 3.23236390479
+NODE 5 COORD 2.03994442247 4.14534044222 5.47496917316
+NODE 6 COORD 4.03084223457 3.23898968807 6.55548562467
+NODE 7 COORD 6.44436403582 3.33465636213 6.21136065882
+NODE 8 COORD 8.04260574482 4.89039614134 5.22244268006
+NODE 9 COORD 8.18419850451 7.29593855857 4.84037101261
+--------------------------------------------------------------STRUCTURE ELEMENTS
+1 BEAM3R HERM2LINE3 1 3 2 MAT 1 TRIADS 1.50130718108 1.35661672402 0.477171338627 -1.21863105979 -2.05469395851 -0.864081375082 1.75047263615 2.16242430293 0.846693286851
+2 BEAM3R HERM2LINE3 3 5 4 MAT 1 TRIADS -1.21863105979 -2.05469395851 -0.864081375082 0.0863366257001 -0.868861065578 -0.445997447356 -0.52331096555 -1.49406240657 -0.676920827978
+3 BEAM3R HERM2LINE3 5 7 6 MAT 1 TRIADS 0.0863366257001 -0.868861065578 -0.445997447356 1.08651489404 0.517671463843 0.11327033771 0.622783907829 -0.196088893105 -0.181481669134
+4 BEAM3R HERM2LINE3 7 9 8 MAT 1 TRIADS 1.08651489404 0.517671463843 0.11327033771 1.73107566054 2.07231804729 0.804238861004 1.46615833297 1.27206928616 0.439606672431
+------------------------------------------------------------------FLUID ELEMENTS
+-----------------------------------------------------------------------------END

tests/reference-files/test_mesh_creation_functions_helix_twist_angle_right_angle_reference.dat

@@ -0,0 +1,30 @@
+// -----------------------------------------------------------------------------
+// This input file was created with MeshPy.
+// Copyright (c) 2018-2024
+//     Ivo Steinbrecher
+//     Institute for Mathematics and Computer-Based Simulation
+//     Universitaet der Bundeswehr Muenchen
+//     https://www.unibw.de/imcs-en
+// -----------------------------------------------------------------------------
+-----------------------------------------------------------------------MATERIALS
+MAT 1 MAT_BeamReissnerElastHyper YOUNG 100000.0 POISSONRATIO 0.0 DENS 0.0 CROSSAREA 0.7853981633974483 SHEARCORR 1.0 MOMINPOL 0.09817477042468103 MOMIN2 0.04908738521234052 MOMIN3 0.04908738521234052
+-------------------------------------------------------------DNODE-NODE TOPOLOGY
+NODE 5 DNODE 1
+NODE 1 DNODE 2
+-------------------------------------------------------------DLINE-NODE TOPOLOGY
+NODE 1 DLINE 1
+NODE 2 DLINE 1
+NODE 3 DLINE 1
+NODE 4 DLINE 1
+NODE 5 DLINE 1
+---------------------------------------------------------------------NODE COORDS
+NODE 1 COORD 2 2 1
+NODE 2 COORD 3.44337567297 3.44337567297 2.44337567297
+NODE 3 COORD 4.88675134595 4.88675134595 3.88675134595
+NODE 4 COORD 6.33012701892 6.33012701892 5.33012701892
+NODE 5 COORD 7.7735026919 7.7735026919 6.7735026919
+--------------------------------------------------------------STRUCTURE ELEMENTS
+1 BEAM3R HERM2LINE3 1 3 2 MAT 1 TRIADS -0.241535826833 -0.75993556603 0.58311906894 -0.241535826833 -0.75993556603 0.58311906894 -0.241535826833 -0.75993556603 0.58311906894
+2 BEAM3R HERM2LINE3 3 5 4 MAT 1 TRIADS -0.241535826833 -0.75993556603 0.58311906894 -0.241535826833 -0.75993556603 0.58311906894 -0.241535826833 -0.75993556603 0.58311906894
+------------------------------------------------------------------FLUID ELEMENTS
+-----------------------------------------------------------------------------END