Merge pull request #300 from isteinbrecher/arc_length_alternative

Add general functionality to compute arc lengths for created beam filaments (alternative to #298)

Author: isteinbrecher

src/meshpy/core/element_beam.py

@@ -22,6 +22,8 @@
 """This file defines the base beam element in MeshPy."""
 
 from typing import Any as _Any
+from typing import Callable as _Callable
+from typing import List as _List
 from typing import Optional as _Optional
 
 import numpy as _np
@@ -30,6 +32,7 @@
 from meshpy.core.conf import mpy as _mpy
 from meshpy.core.element import Element as _Element
 from meshpy.core.node import NodeCosserat as _NodeCosserat
+from meshpy.core.rotation import Rotation as _Rotation
 from meshpy.core.vtk_writer import add_point_data_node_sets as _add_point_data_node_sets
 
 
@@ -51,34 +54,46 @@ def __init__(self, material=None, nodes=None):
         super().__init__(nodes=nodes, material=material)
 
     def create_beam(
-        self, beam_function, *, start_node=None, end_node=None, relative_twist=None
-    ):
+        self,
+        beam_function: _Callable,
+        *,
+        start_node: _Optional[_NodeCosserat] = None,
+        end_node: _Optional[_NodeCosserat] = None,
+        relative_twist: _Optional[_Rotation] = None,
+        set_nodal_arc_length: bool = False,
+        nodal_arc_length_offset: _Optional[float] = None,
+    ) -> _List[_NodeCosserat]:
         """Create the nodes for this beam element. The function returns a list
         with the created nodes.
 
         In the case of start_node and end_node, it is checked, that the
         function and the node have the same coordinates and rotations.
 
-        Args
-        ----
-        beam_function: function(xi)
-            Returns the position and rotation of the beam along the local
-            coordinate xi.
-        start_node: Node
-            If this argument is given, this is the node of the beam at xi=-1.
-        end_node: Node
-            If this argument is given, this is the node of the beam at xi=1.
-        relative_twist: Rotation
-            Apply this relative rotation to all created nodes. This can be used to
-            "twist" the created beam to match the rotation of given nodes.
+        Args:
+            beam_function: function(xi)
+                Returns the position, rotation and (optionally) arc length of the
+                beam along the local coordinate xi. If no arc lengths is provided,
+                the returned value should simply be None.
+            start_node: Node
+                If this argument is given, this is the node of the beam at xi=-1.
+            end_node: Node
+                If this argument is given, this is the node of the beam at xi=1.
+            relative_twist: Rotation
+                Apply this relative rotation to all created nodes. This can be used to
+                "twist" the created beam to match the rotation of given nodes.
+            set_nodal_arc_length:
+                Flag if the arc length in the created nodes should be set.
+            nodal_arc_length_offset:
+                Offset of the stored nodal arc length w.r.t. to the one generated
+                by the function.
         """
 
         if len(self.nodes) > 0:
             raise ValueError("The beam should not have any local nodes yet!")
 
-        def check_node(node, pos, rot, name):
-            """Check if the given node matches with the position and
-            rotation."""
+        def check_node(node, pos, rot, arc_length, name):
+            """Check if the given node matches with the position and rotation
+            and optionally also the arc length."""
 
             if _np.linalg.norm(pos - node.coordinates) > _mpy.eps_pos:
                 raise ValueError(
@@ -88,31 +103,42 @@ def check_node(node, pos, rot, name):
             if not node.rotation == rot:
                 raise ValueError(f"{name} rotation does not match with function!")
 
+            if arc_length is not None:
+                if _np.abs(node.arc_length - arc_length) > _mpy.eps_pos:
+                    raise ValueError(
+                        f"Arc lengths don't match, got {node.arc_length} and {arc_length}"
+                    )
+
         # Flags if nodes are given
         has_start_node = start_node is not None
         has_end_node = end_node is not None
 
         # Loop over local nodes.
+        arc_length = None
         for i, xi in enumerate(self.nodes_create):
             # Get the position and rotation at xi
-            pos, rot = beam_function(xi)
+            pos, rot, arc_length_from_function = beam_function(xi)
             if relative_twist is not None:
                 rot = rot * relative_twist
+            if set_nodal_arc_length:
+                arc_length = arc_length_from_function + nodal_arc_length_offset
 
             # Check if the position and rotation match existing nodes
             if i == 0 and has_start_node:
-                check_node(start_node, pos, rot, "start_node")
+                check_node(start_node, pos, rot, arc_length, "start_node")
                 self.nodes = [start_node]
             elif (i == len(self.nodes_create) - 1) and has_end_node:
-                check_node(end_node, pos, rot, "end_node")
+                check_node(end_node, pos, rot, arc_length, "end_node")
 
             # Create the node
             if (i > 0 or not has_start_node) and (
                 i < len(self.nodes_create) - 1 or not has_end_node
             ):
                 is_middle_node = 0 < i < len(self.nodes_create) - 1
                 self.nodes.append(
-                    _NodeCosserat(pos, rot, is_middle_node=is_middle_node)
+                    _NodeCosserat(
+                        pos, rot, is_middle_node=is_middle_node, arc_length=arc_length
+                    )
                 )
 
         # Get a list with the created nodes.

src/meshpy/core/node.py

@@ -124,12 +124,15 @@ class NodeCosserat(Node):
     """This object represents a Cosserat node in the mesh, i.e., it contains
     three positions and three rotations."""
 
-    def __init__(self, coordinates, rotation, **kwargs):
+    def __init__(self, coordinates, rotation, *, arc_length=None, **kwargs):
         super().__init__(coordinates, **kwargs)
 
         # Rotation of this node.
         self.rotation = rotation.copy()
 
+        # Arc length along the filament that this beam is a part of
+        self.arc_length = arc_length
+
     def rotate(self, rotation, *, origin=None, only_rotate_triads=False):
         """Rotate this node.
 

src/meshpy/mesh_creation_functions/beam_basic_geometry.py

@@ -99,7 +99,11 @@ def beam_function(xi):
             coordinate xi."""
             point_a = start_point + parameter_a * direction
             point_b = start_point + parameter_b * direction
-            return (0.5 * (1 - xi) * point_a + 0.5 * (1 + xi) * point_b, rotation)
+            pos = 0.5 * (1 - xi) * point_a + 0.5 * (1 + xi) * point_b
+            arc_length = (
+                0.5 * (1 - xi) * parameter_a + 0.5 * (1 + xi) * parameter_b
+            ) * line_length
+            return (pos, rotation, arc_length)
 
         return beam_function
 
@@ -258,7 +262,7 @@ def beam_function(xi):
             arc_rotation = _Rotation(axis, phi)
             rot = arc_rotation * start_rotation
             pos = center + arc_rotation * distance
-            return (pos, rot)
+            return (pos, rot, phi * radius)
 
         return beam_function
 

src/meshpy/mesh_creation_functions/beam_curve.py

@@ -246,7 +246,7 @@ def get_beam_position_and_rotation_at_xi(xi):
                 self.S_start_newton = S
 
                 # Return the needed values for beam creation.
-                return (pos, rot)
+                return (pos, rot, S)
 
             return get_beam_position_and_rotation_at_xi
 

src/meshpy/mesh_creation_functions/beam_generic.py

@@ -51,10 +51,12 @@ def create_beam_mesh_function(
     n_el: _Optional[int] = None,
     l_el: _Optional[float] = None,
     interval_length: _Optional[float] = None,
+    set_nodal_arc_length: bool = False,
+    nodal_arc_length_offset: _Optional[float] = None,
     node_positions_of_elements: _Optional[_List[float]] = None,
     start_node: _Optional[_Union[_NodeCosserat, _GeometrySet]] = None,
     end_node: _Optional[_Union[_NodeCosserat, _GeometrySet]] = None,
-    close_beam: _Optional[bool] = False,
+    close_beam: bool = False,
     vtk_cell_data: _Optional[_Dict[str, _Tuple]] = None,
 ) -> _GeometryName:
     """Generic beam creation function.
@@ -77,6 +79,10 @@ def create_beam_mesh_function(
             point_b (both within the interval) and return a function(xi) that
             calculates the position and rotation along the beam, where
             point_a -> xi = -1 and point_b -> xi = 1.
+
+            Usually, the Jacobian of the returned position field should be a unit
+            vector. Otherwise, the nodes may be spaced in an undesired way. All
+            standard mesh creation functions in MeshPy fulfill this property.
         interval:
             Start and end values for interval that will be used to create the
             beam.
@@ -90,6 +96,15 @@ def create_beam_mesh_function(
             created.
         interval_length:
             Total length of the interval. Is required when the option l_el is given.
+        set_nodal_arc_length:
+            Flag if the arc length along the beam filament is set in the created
+            nodes. It is ensured that the arc length is consistent with possible
+            given start/end nodes.
+        nodal_arc_length_offset:
+            Offset of the stored nodal arc length w.r.t. to the one generated by
+            the function. Defaults to 0, the arc length set in the start node, or
+            the arc length in the end node minus total length (such that the arc
+            length at the end node matches).
         node_positions_of_elements:
             A list of normalized positions (within [0,1] and in ascending order)
             that define the boundaries of beam elements along the created curve.
@@ -138,6 +153,17 @@ def create_beam_mesh_function(
             'The arguments "close_beam" and "end_node" are mutually exclusive'
         )
 
+    if set_nodal_arc_length:
+        if close_beam:
+            raise ValueError(
+                "The flags 'set_nodal_arc_length' and 'close_beam' are mutually exclusive."
+            )
+    elif nodal_arc_length_offset is not None:
+        raise ValueError(
+            'Providing the argument "nodal_arc_length_offset" without setting '
+            '"set_nodal_arc_length" to True does not make sense.'
+        )
+
     # Cases where we have equally spaced elements
     if n_el is not None or l_el is not None:
         if l_el is not None:
@@ -150,13 +176,10 @@ def create_beam_mesh_function(
         elif n_el is None:
             raise ValueError("n_el should not be None at this point")
 
-        interval_node_positions_of_elements = [
-            interval[0] + i_node * (interval[1] - interval[0]) / n_el
-            for i_node in range(n_el + 1)
-        ]
+        node_positions_of_elements = [i_node / n_el for i_node in range(n_el + 1)]
     # A list for the element node positions was provided
     elif node_positions_of_elements is not None:
-        # Check that the given positions are in ascending order and start with 1 and end with 0
+        # Check that the given positions are in ascending order and start with 0 and end with 1
         for index, value, name in zip([0, -1], [0, 1], ["First", "Last"]):
             if not _np.isclose(
                 value,
@@ -174,9 +197,11 @@ def create_beam_mesh_function(
             raise ValueError(
                 f"The given node_positions_of_elements must be in ascending order. Got {node_positions_of_elements}"
             )
-        interval_node_positions_of_elements = interval[0] + (
-            interval[1] - interval[0]
-        ) * _np.asarray(node_positions_of_elements)
+
+    # Get the scale the node positions to the interval coordinates
+    interval_node_positions_of_elements = interval[0] + (
+        interval[1] - interval[0]
+    ) * _np.asarray(node_positions_of_elements)
 
     # We need to make sure we have the number of elements for the case a given end node
     n_el = len(interval_node_positions_of_elements) - 1
@@ -234,7 +259,7 @@ def get_relative_twist(rotation_node, rotation_function, name):
         start_node = _get_single_node(start_node)
         nodes = [start_node]
         check_given_node(start_node)
-        _, start_rotation = function_over_whole_interval(-1.0)
+        _, start_rotation, _ = function_over_whole_interval(-1.0)
         relative_twist_start = get_relative_twist(
             start_node.rotation, start_rotation, "start"
         )
@@ -243,9 +268,24 @@ def get_relative_twist(rotation_node, rotation_function, name):
     if end_node is not None:
         end_node = _get_single_node(end_node)
         check_given_node(end_node)
-        _, end_rotation = function_over_whole_interval(1.0)
+        _, end_rotation, _ = function_over_whole_interval(1.0)
         relative_twist_end = get_relative_twist(end_node.rotation, end_rotation, "end")
 
+    # Get the start value for the arc length functionality
+    if set_nodal_arc_length:
+        if nodal_arc_length_offset is not None:
+            # Let's use the given value, if it does not match with possible given
+            # start or end nodes, the check in the create beam function will detect
+            # that.
+            pass
+        elif start_node is not None and start_node.arc_length is not None:
+            nodal_arc_length_offset = start_node.arc_length
+        elif end_node is not None and end_node.arc_length is not None:
+            nodal_arc_length_offset = end_node.arc_length - interval_length
+        else:
+            # Default value
+            nodal_arc_length_offset = 0.0
+
     # Check if a relative twist has to be applied
     if relative_twist_start is not None and relative_twist_end is not None:
         if relative_twist_start == relative_twist_end:
@@ -262,7 +302,7 @@ def get_relative_twist(rotation_node, rotation_function, name):
         relative_twist = None
 
     # Create the beams.
-    for i_el in range(len(interval_node_positions_of_elements) - 1):
+    for i_el in range(n_el):
         # If the beam is closed with itself, set the end node to be the
         # first node of the beam. This is done when the second element is
         # created, as the first node already exists here.
@@ -295,6 +335,8 @@ def get_relative_twist(rotation_node, rotation_function, name):
                 start_node=first_node,
                 end_node=last_node,
                 relative_twist=relative_twist,
+                set_nodal_arc_length=set_nodal_arc_length,
+                nodal_arc_length_offset=nodal_arc_length_offset,
             )
         )
 

src/meshpy/space_time/beam_to_space_time.py

@@ -54,10 +54,9 @@ class NodeCosseratSpaceTime(_NodeCosserat):
     We add the 4th dimension time as a class variable.
     """
 
-    def __init__(self, coordinates, rotation, time, *, arc_length=None, **kwargs):
+    def __init__(self, coordinates, rotation, time, **kwargs):
         super().__init__(coordinates, rotation, **kwargs)
         self.time = time
-        self.arc_length = arc_length
 
     def _get_dat(self):
         """Return the line that represents this node in the input file.