Merge pull request #39 from isteinbrecher/improve-rotations

Improve rotations

meshpy/mesh_creation_functions/__init__.py

@@ -36,7 +36,8 @@
 # Basic geometry functions
 from .beam_basic_geometry import (
     create_beam_mesh_line,
-    create_beam_mesh_arc_segment,
+    create_beam_mesh_arc_segment_via_axis,
+    create_beam_mesh_arc_segment_via_rotation,
     create_beam_mesh_arc_segment_2d,
     create_beam_mesh_line_at_node,
     create_beam_mesh_arc_at_node,
@@ -74,7 +75,8 @@
 __all__ = [
     # Base geometry.
     "create_beam_mesh_line",
-    "create_beam_mesh_arc_segment",
+    "create_beam_mesh_arc_segment_via_axis",
+    "create_beam_mesh_arc_segment_via_rotation",
     "create_beam_mesh_arc_segment_2d",
     "create_beam_mesh_line_at_node",
     "create_beam_mesh_arc_at_node",

meshpy/mesh_creation_functions/beam_basic_geometry.py

@@ -123,12 +123,18 @@ def beam_function(xi):
     )
 
 
-def create_beam_mesh_arc_segment(
+def create_beam_mesh_arc_segment_via_rotation(
     mesh, beam_object, material, center, axis_rotation, radius, angle, **kwargs
 ):
     """
     Generate a circular segment of beam elements.
 
+    The circular segment is defined via a rotation, specifying the "initial"
+    triad of the beam at the beginning of the arc.
+
+    This function exists for compatibility reasons with older MeshPy implementations.
+    The user is encouraged to use the newer implementation create_beam_mesh_arc_segment_via_axis
+
     Args
     ----
     mesh: Mesh
@@ -167,20 +173,99 @@ def create_beam_mesh_arc_segment(
         with all nodes of the line.
     """
 
+    # Convert the input to the one for create_beam_mesh_arc_segment_via_axis
+    axis = axis_rotation * [0, 0, 1]
+    start_point = center + radius * (axis_rotation * np.array([0, -1, 0]))
+    return create_beam_mesh_arc_segment_via_axis(
+        mesh, beam_object, material, axis, center, start_point, angle, **kwargs
+    )
+
+
+def create_beam_mesh_arc_segment_via_axis(
+    mesh,
+    beam_object,
+    material,
+    axis,
+    axis_point,
+    start_point,
+    angle,
+    *,
+    start_node=None,
+    **kwargs
+):
+    """
+    Generate a circular segment of beam elements.
+
+    The arc is defined via a rotation axis, a point on the rotation axis a starting
+    point, as well as the angle of the arc segment.
+
+    Args
+    ----
+    mesh: Mesh
+        Mesh that the arc 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 arc.
+    axis_point: np.array, list
+        Point lying on the rotation axis. Does not have to be the center of the arc.
+    start_point: np.array, list
+        Start point of the arc.
+    angle: float
+        The central angle of this segment in radians.
+
+    **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. This requires the option interval_length
+        to be set. 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.
+    """
+
     # The angle can not be negative with the current implementation.
     if angle <= 0.0:
-        raise ValueError("The angle for a beam segment has to be a positive number!")
+        raise ValueError(
+            "The angle for a beam arc segment has to be a positive number!"
+        )
+
+    # Shortest distance from the given point to the axis of rotation gives
+    # the "center" of the arc
+    axis = np.array(axis) / np.linalg.norm(axis)
+    diff = start_point - axis_point
+    distance = diff - np.dot(np.dot(diff, axis), axis)
+    radius = np.linalg.norm(distance)
+    center = start_point - distance
+
+    # Get the rotation at the start
+    if start_node is None:
+        tangent = np.cross(axis, distance)
+        tangent /= np.linalg.norm(tangent)
+        start_rotation = Rotation.from_rotation_matrix(
+            np.transpose(np.array([tangent, -distance / radius, axis]))
+        )
+    else:
+        start_rotation = get_single_node(start_node).rotation
 
     def get_beam_geometry(alpha, beta):
-        """
-        Return a function for the position and rotation along the beam
-        axis.
-        """
+        """Return a function for the position and rotation along the beam axis"""
 
         def beam_function(xi):
             phi = 0.5 * (xi + 1) * (beta - alpha) + alpha
-            rot = axis_rotation * Rotation([0, 0, 1], phi)
-            pos = center + radius * (rot * [0, -1, 0])
+            arc_rotation = Rotation(axis, phi)
+            rot = arc_rotation * start_rotation
+            pos = center + arc_rotation * distance
             return (pos, rot)
 
         return beam_function
@@ -193,6 +278,7 @@ def beam_function(xi):
         function_generator=get_beam_geometry,
         interval=[0.0, angle],
         interval_length=angle * radius,
+        start_node=start_node,
         **kwargs,
     )
 
@@ -246,25 +332,22 @@ def create_beam_mesh_arc_segment_2d(
 
     # Check if the beam is in clockwise or counter clockwise direction.
     angle = phi_end - phi_start
-    clockwise = not np.sign(angle) == 1
-
-    # Create rotation for the general arc segment function.
-    axis_rotation = Rotation([0, 0, 1], np.pi * 0.5 + phi_start)
-    if clockwise:
-        # If the beam is in counter clockwise direction, we have to rotate the
-        # rotation axis around its first basis vector - this will result in the
-        # arc facing in the other direction. Additionally we have to rotate
-        # around the z-axis for the angles to fit.
-        t1 = [-np.sin(phi_start), np.cos(phi_start), 0.0]
-        axis_rotation = Rotation([0, 0, 1], np.pi) * Rotation(t1, np.pi) * axis_rotation
-
-    return create_beam_mesh_arc_segment(
+    axis = np.array([0, 0, 1])
+    start_point = center + radius * (Rotation(axis, phi_start) * [1, 0, 0])
+
+    counter_clockwise = np.sign(angle) == 1
+    if not counter_clockwise:
+        # If the beam is not in counter clockwise direction, we have to flip
+        # the rotation axis.
+        axis = -1.0 * axis
+
+    return create_beam_mesh_arc_segment_via_axis(
         mesh,
         beam_object,
         material,
+        axis,
         center,
-        axis_rotation,
-        radius,
+        start_point,
         np.abs(angle),
         **kwargs,
     )
@@ -388,18 +471,13 @@ def create_beam_mesh_arc_at_node(
     center_direction *= 1.0 / np.linalg.norm(center_direction)
     center = start_node.coordinates - center_direction * radius
 
-    # Create rotation for the general arc segment function
-    axis_rotation = Rotation.from_rotation_matrix(
-        np.transpose(np.array([tangent, -center_direction, arc_axis_normal]))
-    )
-
-    return create_beam_mesh_arc_segment(
+    return create_beam_mesh_arc_segment_via_axis(
         mesh,
         beam_object,
         material,
+        arc_axis_normal,
         center,
-        axis_rotation,
-        radius,
+        start_node.coordinates,
         np.abs(angle),
         start_node=start_node,
         **kwargs,

meshpy/mesh_creation_functions/beam_stent.py

@@ -42,7 +42,10 @@
 from ..container import GeometryName
 from ..geometry_set import GeometrySet
 from ..utility import get_min_max_nodes
-from .beam_basic_geometry import create_beam_mesh_arc_segment, create_beam_mesh_line
+from .beam_basic_geometry import (
+    create_beam_mesh_arc_segment_via_rotation,
+    create_beam_mesh_line,
+)
 
 
 def create_stent_cell(
@@ -103,7 +106,7 @@ def add_line(pointa, pointb, n_el_line):
 
     def add_segment(center, axis_rotation, radius, angle, n_el_segment):
         """Shortcut to add arc segment."""
-        return create_beam_mesh_arc_segment(
+        return create_beam_mesh_arc_segment_via_rotation(
             mesh,
             beam_object,
             material,

meshpy/rotation.py

@@ -66,8 +66,10 @@ def __init__(self, *args):
             # Identity element.
             self.q[0] = 1
         elif len(args) == 1 and len(args[0]) == 4:
-            # Set from quaternion.
-            self.q[:] = np.array(args[0])
+            # Set directly from quaternion
+            # To avoid error accumulation, normalize the quaternion here
+            q = np.array(args[0])
+            self.q[:] = q / np.linalg.norm(q)
         elif len(args) == 2:
             # Set from vector and rotation angle.
             vector = args[0]
@@ -87,40 +89,29 @@ def __init__(self, *args):
     def from_rotation_matrix(cls, R):
         """
         Create the object from a rotation matrix.
-        The code is base on:
-            https://en.wikipedia.org/wiki/Rotation_matrix#Quaternion
+        The code is base on Spurriers algorithm:
+            R. A. Spurrier (1978): “Comment on “Singularity-free extraction of a quaternion from a direction-cosine matrix”
         """
 
         q = np.zeros(4)
         trace = np.trace(R)
-        if trace > 0:
-            r = np.sqrt(1 + trace)
-            s = 0.5 / r
-            q[0] = 0.5 * r
-            q[1] = (R[2, 1] - R[1, 2]) * s
-            q[2] = (R[0, 2] - R[2, 0]) * s
-            q[3] = (R[1, 0] - R[0, 1]) * s
-        elif R[0, 0] > R[1, 1] and R[0, 0] > R[2, 2]:
-            r = np.sqrt(1 + R[0, 0] - R[1, 1] - R[2, 2])
-            s = 0.5 / r
-            q[0] = (R[2, 1] - R[1, 2]) * s
-            q[1] = 0.5 * r
-            q[2] = (R[0, 1] + R[1, 0]) * s
-            q[3] = (R[0, 2] + R[2, 0]) * s
-        elif R[1, 1] > R[2, 2]:
-            r = np.sqrt(1 - R[0, 0] + R[1, 1] - R[2, 2])
-            s = 0.5 / r
-            q[0] = (R[0, 2] - R[2, 0]) * s
-            q[1] = (R[0, 1] + R[1, 0]) * s
-            q[2] = 0.5 * r
-            q[3] = (R[1, 2] + R[2, 1]) * s
+        values = [R[i, i] for i in range(3)]
+        values.append(trace)
+        arg_max = np.argmax(values)
+        if arg_max == 3:
+            q[0] = np.sqrt(trace + 1) * 0.5
+            q[1] = (R[2, 1] - R[1, 2]) / (4 * q[0])
+            q[2] = (R[0, 2] - R[2, 0]) / (4 * q[0])
+            q[3] = (R[1, 0] - R[0, 1]) / (4 * q[0])
         else:
-            r = np.sqrt(1 - R[0, 0] - R[1, 1] + R[2, 2])
-            s = 0.5 / r
-            q[0] = (R[1, 0] - R[0, 1]) * s
-            q[1] = (R[0, 2] + R[2, 0]) * s
-            q[2] = (R[1, 2] + R[2, 1]) * s
-            q[3] = 0.5 * r
+            i_index = arg_max
+            j_index = (i_index + 1) % 3
+            k_index = (i_index + 2) % 3
+            q_i = np.sqrt(R[i_index, i_index] * 0.5 + (1 - trace) * 0.25)
+            q[0] = (R[k_index, j_index] - R[j_index, k_index]) / (4 * q_i)
+            q[i_index + 1] = q_i
+            q[j_index + 1] = (R[j_index, i_index] + R[i_index, j_index]) / (4 * q_i)
+            q[k_index + 1] = (R[k_index, i_index] + R[i_index, k_index]) / (4 * q_i)
 
         return cls(q)
 

tests/performance_testing.py

@@ -178,26 +178,11 @@ class TestPerformance(object):
 
     # Set expected test times.
     expected_times = {}
-    expected_times["adonis"] = {
-        "cubitpy_create_solid": 3.2,
-        "meshpy_load_solid": 1.1,
-        "meshpy_load_solid_full": 3.5,
-        "meshpy_create_beams": 7.2,
-        "meshpy_rotate": 0.6,
-        "meshpy_translate": 0.6,
-        "meshpy_reflect": 0.7,
-        "meshpy_wrap_around_cylinder": 3.0,
-        "meshpy_wrap_around_cylinder_without_check": 0.9,
-        "meshpy_find_close_nodes": 2.0,
-        "meshpy_write_dat": 12.5,
-        "meshpy_write_vtk": 19,
-        "geometric_search_find_nodes_brute_force": 0.05,
-    }
     expected_times["ares.bauv.unibw-muenchen.de"] = {
         "cubitpy_create_solid": 4.0,
         "meshpy_load_solid": 0.9,
         "meshpy_load_solid_full": 3.0,
-        "meshpy_create_beams": 7.2,
+        "meshpy_create_beams": 9.0,
         "meshpy_rotate": 0.6,
         "meshpy_translate": 0.5,
         "meshpy_reflect": 0.7,
@@ -208,21 +193,6 @@ class TestPerformance(object):
         "meshpy_write_vtk": 24.0,
         "geometric_search_find_nodes_brute_force": 0.05,
     }
-    expected_times["sisyphos.bauv.unibw-muenchen.de"] = {
-        "cubitpy_create_solid": 3.0,
-        "meshpy_load_solid": 0.9,
-        "meshpy_load_solid_full": 2.8,
-        "meshpy_create_beams": 6.0,
-        "meshpy_rotate": 0.6,
-        "meshpy_translate": 0.5,
-        "meshpy_reflect": 0.7,
-        "meshpy_wrap_around_cylinder": 2.0,
-        "meshpy_wrap_around_cylinder_without_check": 0.7,
-        "meshpy_find_close_nodes": 1.6,
-        "meshpy_write_dat": 10.5,
-        "meshpy_write_vtk": 17.0,
-        "geometric_search_find_nodes_brute_force": 0.05,
-    }
 
     def __init__(self):
         """

tests/testing_abaqus.py

@@ -139,6 +139,8 @@ def test_abaqus_frame(self):
                 AbaqusBeamNormalDefinition.smallest_rotation_of_triad_at_first_node
             ),
             extension="inp",
+            split_string=",",
+            atol=1e-15,
         )