Merge pull request #89 from isteinbrecher/improve_beam_nurbs_function

Improve beam nurbs function

meshpy/mesh_creation_functions/beam_nurbs.py

@@ -41,7 +41,73 @@
 from .beam_curve import create_beam_mesh_curve
 
 
-def create_beam_mesh_from_nurbs(mesh, beam_object, material, curve, tol=None, **kwargs):
+def get_nurbs_curve_function_and_jacobian_for_integration(
+    curve, curve_start, curve_end, tol=None
+):
+    """Return function objects for evaluating the curve and the derivative. These functions are
+    used in the curve integration. It can happen that the integration algorithm has to evaluate
+    the curve outside of the defined domain. This usually leads to errors in common NURBS
+    packages. Therefore, we check for this evaluation outside of the parameter domain here and
+    perform a linear extrapolation.
+
+    Args
+    ----
+    curve: geomdl object
+        Curve that is used to describe the beam centerline.
+    curve_start, curve_end: (float)
+        Start and end parameter coordinate for the NURBS curve.
+    tol: float
+        Tolerance for checking if point is close to the start or end of the
+        interval. If None is given, use the default tolerance from mpy.
+    """
+
+    if tol is None:
+        tol = mpy.eps_pos
+
+    def eval_r(t):
+        """Evaluate the position along the curve."""
+        return np.array(curve.derivatives(u=t, order=0)[0])
+
+    def eval_rp(t):
+        """Evaluate the derivative along the curve."""
+        return np.array(curve.derivatives(u=t, order=1)[1])
+
+    def function(t):
+        """Convert the curve to a function that can be used for beam
+        generation."""
+
+        if curve_start <= t <= curve_end:
+            return eval_r(t)
+        elif t < curve_start and np.abs(t - curve_start) < tol:
+            diff = t - curve_start
+            return eval_r(curve_start) + diff * eval_rp(curve_start)
+        elif t > curve_end and np.abs(t - curve_end) < tol:
+            diff = t - curve_end
+            return eval_r(curve_end) + diff * eval_rp(curve_end)
+        raise ValueError(
+            "Can not evaluate the curve function outside of the interval (plus tolerances).\n"
+            f"Abs diff start: {np.abs(curve_start - t)}\nAbs diff end: {np.abs(t - curve_end)}"
+        )
+
+    def jacobian(t):
+        """Convert the spline to a Jacobian function that can be used for curve
+        generation. There is no tolerance here, since the integration algorithms sometimes
+        evaluate the derivative far outside the interval."""
+
+        if curve_start <= t <= curve_end:
+            return eval_rp(t)
+        elif t < curve_start:
+            return eval_rp(curve_start)
+        elif curve_end < t:
+            return eval_rp(curve_end)
+        raise ValueError("Should not happen")
+
+    return function, jacobian
+
+
+def create_beam_mesh_from_nurbs(
+    mesh, beam_object, material, curve, *, tol=None, **kwargs
+):
     """
     Generate a beam from a NURBS curve.
 
@@ -57,7 +123,7 @@ def create_beam_mesh_from_nurbs(mesh, beam_object, material, curve, tol=None, **
         Curve that is used to describe the beam centerline.
     tol: float
         Tolerance for checking if point is close to the start or end of the
-        interval.
+        interval. If None is given, use the default tolerance from mpy.
 
     **kwargs (for all of them look into create_beam_mesh_function)
     ----
@@ -76,63 +142,19 @@ def create_beam_mesh_from_nurbs(mesh, beam_object, material, curve, tol=None, **
         with all nodes of the curve.
     """
 
-    # Get the tolerance.
-    if tol is None:
-        tol = mpy.eps_pos
-
-    # Get start and end values of the curve parameter space.
     curve_start = np.min(curve.knotvector)
     curve_end = np.max(curve.knotvector)
+    function, jacobian = get_nurbs_curve_function_and_jacobian_for_integration(
+        curve, curve_start, curve_end, tol=tol
+    )
 
-    def eval_r(t):
-        """Evaluate the position along the curve."""
-        return np.array(curve.derivatives(u=t, order=0)[0])
-
-    def eval_rp(t):
-        """Evaluate the derivative along the curve."""
-        return np.array(curve.derivatives(u=t, order=1)[1])
-
-    def function(t):
-        """Convert the curve to a function that can be used for beam
-        generation."""
-
-        # Due to numeric tolerances it is possible that the position has to be
-        # evaluated outside the interval.
-        if curve_start <= t and t <= curve_end:
-            return eval_r(t)
-        elif np.abs(curve_start - t) < tol:
-            return eval_r(curve_start)
-        elif np.abs(t - curve_end) < tol:
-            return eval_r(curve_end)
-        else:
-            raise ValueError(
-                (
-                    "Can not evaluate the curve function outside of the interval (plus tolerances)."
-                    "\nAbs diff start: {}\nAbs diff end: {}"
-                ).format(np.abs(curve_start - t), np.abs(t - curve_end))
-            )
-
-    def jacobian(t):
-        """Convert the spline to a Jacobian function that can be used for curve
-        generation."""
-
-        # In the curve integration it is possible that the Jacobian is
-        # evaluated outside the interval. In that case use the values at the
-        # limits of the interval.
-        if curve_start <= t and t <= curve_end:
-            return eval_rp(t)
-        elif t < curve_start:
-            return eval_rp(curve_start)
-        elif curve_end < t:
-            return eval_rp(curve_end)
-
-    # Create the beams.
+    # Create the beams
     return create_beam_mesh_curve(
         mesh,
         beam_object,
         material,
         function,
         [curve_start, curve_end],
         function_derivative=jacobian,
-        **kwargs
+        **kwargs,
     )

tests/testing_mesh_creation_functions.py

@@ -37,6 +37,8 @@
 import numpy as np
 import autograd.numpy as npAD
 import os
+from geomdl import NURBS
+from geomdl import utilities
 
 # Meshpy imports.
 from meshpy import (
@@ -68,6 +70,9 @@
     create_beam_mesh_from_nurbs,
     create_beam_mesh_helix,
 )
+from meshpy.mesh_creation_functions.beam_nurbs import (
+    get_nurbs_curve_function_and_jacobian_for_integration,
+)
 
 # Testing imports.
 from utils import testing_input, compare_test_result, compare_strings
@@ -127,6 +132,16 @@ def helix(t):
     return helix
 
 
+def create_testing_nurbs_curve():
+    """Create a NURBS curve used for testing"""
+
+    curve = NURBS.Curve()
+    curve.degree = 2
+    curve.ctrlpts = [[0, 0, 0], [1, 2, -1], [2, 0, 0]]
+    curve.knotvector = utilities.generate_knot_vector(curve.degree, len(curve.ctrlpts))
+    return curve
+
+
 class TestMeshCreationFunctions(unittest.TestCase):
     """
     Test the mesh creation functions.
@@ -357,18 +372,8 @@ def test_mesh_creation_functions_nurbs(self):
         Test the create_beam_mesh_from_nurbs function.
         """
 
-        # Setup the NURBS curve.
-        from geomdl import NURBS
-        from geomdl import utilities
-
-        curve = NURBS.Curve()
-        curve.degree = 2
-        curve.ctrlpts = [[0, 0, 0], [1, 2, -1], [2, 0, 0]]
-        curve.knotvector = utilities.generate_knot_vector(
-            curve.degree, len(curve.ctrlpts)
-        )
-
         # Create beam elements.
+        curve = create_testing_nurbs_curve()
         mat = MaterialReissner(radius=0.05)
         mesh = Mesh()
         create_beam_mesh_from_nurbs(mesh, Beam3rHerm2Line3, mat, curve, n_el=3)
@@ -378,6 +383,32 @@ def test_mesh_creation_functions_nurbs(self):
         input_file.add(mesh)
         compare_test_result(self, input_file.get_string(header=False))
 
+    def test_mesh_creation_functions_nurbs_unit(self):
+        """Unittest the function and jacobian creation in the create_beam_mesh_from_nurbs function"""
+
+        curve = create_testing_nurbs_curve()
+        curve_start = np.min(curve.knotvector)
+        curve_end = np.max(curve.knotvector)
+        r, dr = get_nurbs_curve_function_and_jacobian_for_integration(
+            curve, curve_start, curve_end, tol=10
+        )
+
+        t_values = [5.0 / 7.0, -0.3, 1.2]
+        results_r = [
+            [1.4285714285714286, 0.8163265306122449, -0.4081632653061225],
+            [-0.6, -1.2, 0.6],
+            [2.4, -0.8, 0.4],
+        ]
+        results_dr = [
+            [2.0, -1.7142857142857144, 0.8571428571428572],
+            [2.0, 4.0, -2.0],
+            [2.0, -4.0, 2.0],
+        ]
+
+        for t, result_r, result_dr in zip(t_values, results_r, results_dr):
+            self.assertTrue(np.allclose(r(t), result_r, atol=mpy.eps_pos))
+            self.assertTrue(np.allclose(dr(t), result_dr, atol=mpy.eps_pos))
+
     def test_mesh_creation_functions_node_continuation(self):
         """Test that the node continuation function work as expected."""