Merge pull request #237 from isteinbrecher/improve-mesh-creation

Improve array conversions

meshpy/cosserat_curve/cosserat_curve.py

@@ -266,7 +266,7 @@ def get_centerline_positions_and_rotations(
         """
 
         # Get the points that are within the arc length of the given curve.
-        points_on_arc_length = np.array(points_on_arc_length)
+        points_on_arc_length = np.asarray(points_on_arc_length)
         points_in_bounds = np.logical_and(
             points_on_arc_length > self.point_arc_length[0],
             points_on_arc_length < self.point_arc_length[-1],

meshpy/mesh.py

@@ -331,7 +331,7 @@ def reflect(self, normal_vector, origin=None, flip_beams=False):
         """
 
         # Normalize the normal vector.
-        normal_vector = np.array(normal_vector / np.linalg.norm(normal_vector))
+        normal_vector = np.asarray(normal_vector) / np.linalg.norm(normal_vector)
 
         # Get array with all quaternions and positions for the nodes.
         pos = get_nodal_coordinates(self.nodes)
@@ -892,7 +892,7 @@ def display_pyvista(
                 directors = [
                     finite_element_nodes.glyph(
                         geom=arrow,
-                        orient=f"base_vector_{i+1}",
+                        orient=f"base_vector_{i + 1}",
                         scale="cross_section_radius",
                         factor=director_radius_scaling_factor,
                     )

meshpy/mesh_creation_functions/beam_basic_geometry.py

@@ -81,7 +81,9 @@ def create_beam_mesh_line(
     """
 
     # Get geometrical values for this line.
-    direction = np.array(end_point) - np.array(start_point)
+    start_point = np.asarray(start_point)
+    end_point = np.asarray(end_point)
+    direction = end_point - start_point
     line_length = np.linalg.norm(direction)
     t1 = direction / line_length
 
@@ -169,7 +171,7 @@ def create_beam_mesh_arc_segment_via_rotation(
 
     # 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]))
+    start_point = center + radius * (axis_rotation * [0, -1, 0])
     return create_beam_mesh_arc_segment_via_axis(
         mesh, beam_object, material, axis, center, start_point, angle, **kwargs
     )
@@ -234,7 +236,11 @@ def create_beam_mesh_arc_segment_via_axis(
 
     # 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)
+    axis = np.asarray(axis)
+    axis_point = np.asarray(axis_point)
+    start_point = np.asarray(start_point)
+
+    axis = axis / np.linalg.norm(axis)
     diff = start_point - axis_point
     distance = diff - np.dot(np.dot(diff, axis), axis)
     radius = np.linalg.norm(distance)
@@ -444,7 +450,7 @@ def create_beam_mesh_arc_at_node(
     """
 
     # If the angle is negative, the normal is switched
-    arc_axis_normal = np.array(arc_axis_normal)
+    arc_axis_normal = np.asarray(arc_axis_normal)
     if angle < 0:
         arc_axis_normal = -1.0 * arc_axis_normal
 
@@ -554,9 +560,13 @@ def create_beam_mesh_helix(
         )
 
     # determine radius of helix
-    axis_vector = np.array(axis_vector) / np.linalg.norm(np.array(axis_vector))
+    axis_vector = np.asarray(axis_vector)
+    axis_point = np.asarray(axis_point)
+    start_point = np.asarray(start_point)
+
+    axis_vector = axis_vector / np.linalg.norm(axis_vector)
     origin = axis_point + np.dot(
-        np.dot((np.array(start_point) - np.array(axis_point)), axis_vector), axis_vector
+        np.dot(start_point - axis_point, axis_vector), axis_vector
     )
     start_point_origin_vec = start_point - origin
     radius = np.linalg.norm(start_point_origin_vec)

meshpy/mesh_creation_functions/beam_fibers_in_rectangle.py

@@ -71,8 +71,8 @@ def _intersect_line_with_rectangle(
     """
 
     # Convert the input values to np.arrays.
-    start_line = np.array(start_line)
-    direction_line = np.array(direction_line)
+    start_line = np.asarray(start_line)
+    direction_line = np.asarray(direction_line)
 
     # Set definition for the boundary lines of the rectangle. The director is
     # chosen in a way, that the values [0, 1] for the line parameters alpha are

meshpy/rotation.py

@@ -71,7 +71,7 @@ def __init__(self, *args):
             self.q[0] = 1
         elif len(args) == 2:
             # Set from rotation axis and rotation angle.
-            axis = np.array(args[0])
+            axis = np.asarray(args[0])
             phi = args[1]
             norm = np.linalg.norm(axis)
             if norm < mpy.eps_quaternion:
@@ -95,9 +95,10 @@ def from_quaternion(cls, q, *, normalized=False):
             error accumulation.
         """
         rotation = object.__new__(cls)
-        rotation.q = np.array(q, dtype=float)
-        if not normalized:
-            rotation.q /= np.linalg.norm(rotation.q)
+        if normalized:
+            rotation.q = np.array(q)
+        else:
+            rotation.q = np.asarray(q) / np.linalg.norm(q)
         if (not rotation.q.ndim == 1) or (not len(rotation.q) == 4):
             raise ValueError("Got quaternion array with unexpected dimensions")
         return rotation
@@ -154,7 +155,7 @@ def from_rotation_vector(cls, rotation_vector):
         """Create the object from a rotation vector."""
 
         q = np.zeros(4)
-        rotation_vector = np.array(rotation_vector)
+        rotation_vector = np.asarray(rotation_vector)
         phi = np.linalg.norm(rotation_vector)
         q[0] = np.cos(0.5 * phi)
         if phi < mpy.eps_quaternion:
@@ -199,7 +200,7 @@ def get_rotation_matrix(self):
         return R
 
     def get_quaternion(self):
-        """Return the quaternion for this rotation, as tuple."""
+        """Return the quaternion for this rotation, as numpy array (copy)."""
 
         return np.array(self.q)
 
@@ -313,7 +314,7 @@ def __mul__(self, other):
             return Rotation.from_quaternion(added_rotation)
         elif isinstance(other, (list, np.ndarray)) and len(other) == 3:
             # Apply rotation to vector.
-            return np.dot(self.get_rotation_matrix(), np.array(other))
+            return np.dot(self.get_rotation_matrix(), np.asarray(other))
         raise NotImplementedError("Error, not implemented, does not make sense anyway!")
 
     def __eq__(self, other):
@@ -490,7 +491,7 @@ def smallest_rotation(q: Rotation, t):
 
     R_old = q.get_rotation_matrix()
     g1_old = R_old[:, 0]
-    g1 = np.array(t) / np.linalg.norm(t)
+    g1 = np.asarray(t) / np.linalg.norm(t)
 
     # Quaternion components of relative rotation
     q_rel = np.zeros(4)

tests/reference-files/test_mesh_creation_functions_arc_segment_via_axis.dat

@@ -0,0 +1,33 @@
+// -----------------------------------------------------------------------------
+// This input file was created with MeshPy.
+// Copyright (c) 2018-2025
+//     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 -1.0 POISSONRATIO 0.0 DENS 0.0 CROSSAREA 3.141592653589793 SHEARCORR 1 MOMINPOL 1.5707963267948966 MOMIN2 0.7853981633974483 MOMIN3 0.7853981633974483
+-------------------------------------------------------------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 0 0 0
+NODE 2 COORD 0.331792265387 0.0277135368741 0
+NODE 3 COORD 0.654389393592 0.110086107371 0
+NODE 4 COORD 0.958851077208 0.244834876219 0
+NODE 5 COORD 1.23673960614 0.428225478446 0
+NODE 6 COORD 1.48035370639 0.655175511834 0
+NODE 7 COORD 1.68294196962 0.919395388264 0
+--------------------------------------------------------------STRUCTURE ELEMENTS
+1 BEAM3R HERM2LINE3 1 3 2 MAT 1 TRIADS 0 0 0 0 0 0.333333333333 0 0 0.166666666667
+2 BEAM3R HERM2LINE3 3 5 4 MAT 1 TRIADS 0 0 0.333333333333 0 0 0.666666666667 0 0 0.5
+3 BEAM3R HERM2LINE3 5 7 6 MAT 1 TRIADS 0 0 0.666666666667 0 0 1 0 0 0.833333333333

tests/test_mesh_creation_functions.py

@@ -28,8 +28,6 @@
 # SOFTWARE.
 """This script is used to test the mesh creation functions."""
 
-import os
-
 import autograd.numpy as npAD
 import numpy as np
 import pytest
@@ -49,6 +47,7 @@
 from meshpy.mesh_creation_functions import (
     create_beam_mesh_arc_at_node,
     create_beam_mesh_arc_segment_2d,
+    create_beam_mesh_arc_segment_via_axis,
     create_beam_mesh_arc_segment_via_rotation,
     create_beam_mesh_curve,
     create_beam_mesh_from_nurbs,
@@ -136,10 +135,37 @@ def create_testing_nurbs_curve():
     )
 
 
-def test_mesh_creation_functions_arc_segment(
+def test_mesh_creation_functions_arc_segment_via_axis(
+    assert_results_equal, get_corresponding_reference_file_path
+):
+    """Create a circular segment via the axis method and compare it with the
+    reference file."""
+
+    # Create mesh
+    input_file = InputFile()
+    mat = MaterialReissner()
+    radius = 2.0
+    beam_set = create_beam_mesh_arc_segment_via_axis(
+        input_file,
+        Beam3rHerm2Line3,
+        mat,
+        [0, 0, 1],
+        [0, radius, 0],
+        [0, 0, 0],
+        1.0,
+        n_el=3,
+    )
+    input_file.add(beam_set)
+
+    # Check the output.
+    assert_results_equal(get_corresponding_reference_file_path(), input_file)
+
+
+def test_mesh_creation_functions_arc_segment_via_rotation(
     assert_results_equal, get_corresponding_reference_file_path
 ):
-    """Create a circular segment and compare it with the reference file."""
+    """Create a circular segment via the rotation method and compare it with
+    the reference file."""
 
     # Create input file.
     input_file = InputFile()