Task 2: Unit test for specific functions part 3

elastica/joint.py

@@ -726,7 +726,7 @@ def _calculate_contact_forces_self_rod(
     edge_collection_rod_one = _batch_product_k_ik_to_ik(length_rod, tangent_rod)
 
     for i in range(n_points_rod):
-        skip = 1 + np.ceil(0.8 * np.pi * radius_rod[i] / length_rod[i])
+        skip = int(1 + np.ceil(0.8 * np.pi * radius_rod[i] / length_rod[i]))
         for j in range(i - skip, -1, -1):
             radii_sum = radius_rod[i] + radius_rod[j]
             length_sum = length_rod[i] + length_rod[j]

tests/test_contact_specific_functions.py

@@ -8,6 +8,8 @@
     _prune_using_aabbs_rod_rigid_body,
     _prune_using_aabbs_rod_rod,
     _calculate_contact_forces_rod_rigid_body,
+    _calculate_contact_forces_rod_rod,
+    _calculate_contact_forces_self_rod,
 )
 
 
@@ -224,3 +226,195 @@ def test_claculate_contact_forces_rod_rigid_body():
         rtol=tol,
         atol=tol,
     )
+
+
+def test_calculate_contact_forces_rod_rod():
+    "Function to test the calculate contact forces rod rod function"
+
+    "Testing function with analytically verified values"
+
+    tol = 1e-5
+
+    "initializing rod 1 parameters"
+    rod_one_position_collection = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9]])
+    rod_one_radius_collection = np.array([1.3, 0.5])
+    rod_one_length_collection = np.array([2.1, 0.9])
+    rod_one_tangent_collection = np.array([[0, -0.29], [-0.029, 1.12], [-1.04, -1.165]])
+
+    rod_one_internal_forces = np.array(
+        [
+            [0.59155762, -0.0059393, 0.39725841],
+            [-2.45916096, 0.81593963, 1.19209269],
+            [1.46996099, -0.52774132, -0.08287433],
+        ]
+    )
+    rod_one_external_forces = np.array(
+        [
+            [-1.55308218, 1.00714652, 2.04033497],
+            [-0.35890683, -0.22667242, -0.31385005],
+            [-0.66969949, -0.67150666, 0.17822428],
+        ]
+    )
+
+    rod_one_velocity_collection = np.array(
+        [
+            [-0.91917535, -1.04929632, -0.32208191],
+            [0.13459384, -1.77261144, -0.28656239],
+            [1.63574433, 0.31741648, 1.43917516],
+        ]
+    )
+
+    "initializing rod 2 parameters"
+    rod_two_position_collection = np.array(
+        [[1.0000001, 2.0000034, 3], [4, 5, 6], [7, 8, 9]]
+    )
+    rod_two_radius_collection = np.array([2.3, 0.4])
+    rod_two_length_collection = np.array([1.5, 1.0])
+    rod_two_tangent_collection = np.array([[0, -0.39], [-0.29, 1.00], [-1.0, -0.165]])
+
+    rod_two_internal_forces = np.array(
+        [
+            [-0.46552762, 1.16896583, 1.06695832],
+            [0.35104219, 0.5868144, -0.79542854],
+            [1.90896989, -1.7709093, 1.21209849],
+        ]
+    )
+    rod_two_external_forces = np.array(
+        [
+            [0.12615472, -0.12539237, -1.01333332],
+            [0.7193244, 1.12085914, 1.57535336],
+            [1.79006353, 0.20498294, 1.11384582],
+        ]
+    )
+
+    rod_two_velocity_collection = np.array(
+        [
+            [-1.09200111, -0.87381467, 0.44073417],
+            [0.93961598, 1.25513012, 0.09103194],
+            [1.02214026, -0.78790631, -0.74019659],
+        ]
+    )
+
+    "initializing constants"
+    k = 1.0
+    nu = 0.5
+
+    "Function call"
+    _calculate_contact_forces_rod_rod(
+        rod_one_position_collection[..., :-1],
+        rod_one_radius_collection,
+        rod_one_length_collection,
+        rod_one_tangent_collection,
+        rod_one_velocity_collection,
+        rod_one_internal_forces,
+        rod_one_external_forces,
+        rod_two_position_collection[..., :-1],
+        rod_two_radius_collection,
+        rod_two_length_collection,
+        rod_two_tangent_collection,
+        rod_two_velocity_collection,
+        rod_two_internal_forces,
+        rod_two_external_forces,
+        k,
+        nu,
+    )
+
+    "Test values"
+    assert_allclose(
+        rod_one_external_forces,
+        np.array(
+            [
+                [-2.752861, -2.349062, 1.083684],
+                [-0.358907, 0.627338, 0.54016],
+                [-0.669699, 0.387652, 1.237383],
+            ]
+        ),
+        rtol=tol,
+        atol=tol,
+    )
+
+    assert_allclose(
+        rod_two_external_forces,
+        np.array(
+            [
+                [0.23889, 2.687294, 1.573883],
+                [-0.565314, -0.37547, 2.648301],
+                [0.825641, -1.336387, 1.501321],
+            ]
+        ),
+        rtol=tol,
+        atol=tol,
+    )
+
+
+def test_calculate_contact_forces_self_rod():
+    "Function to test the calculate contact forces self rod function"
+
+    "Testing function with analytically verified values"
+
+    tol = 1e-5
+
+    "initializing rod parameters"
+    rod_position_collection = np.array(
+        [[1, 2, 3, 4, 5], [6, 7, 8, 9, 10], [11, 12, 13, 14, 15]]
+    )
+    rod_radius_collection = np.array(
+        [5.12952263, 0.19493561, 0.72904683, 0.21676892, 0.18705351]
+    )
+    rod_length_collection = np.array(
+        [0.89050271, 0.7076813, 0.21078658, 0.95372826, 0.86037329]
+    )
+
+    rod_tangent_collection = np.array(
+        [
+            [-0.92511524, 0.19778438, 0.42897242, 0.28430662, -0.57637184],
+            [1.73622348, 1.55757074, -0.4606567, -1.30228854, -0.34647765],
+            [0.61561226, -0.86529598, -0.9180072, 0.99279484, -2.09424394],
+        ]
+    )
+
+    rod_external_forces = np.array(
+        [
+            [-0.77306421, -0.25648047, -0.93419262, -0.77665042, -0.33345937],
+            [-1.04834225, -1.92250527, -1.46505411, -0.71192403, -0.99256648],
+            [0.33465609, 1.22871475, 0.06250578, -0.49531749, 0.58044695],
+        ]
+    )
+
+    rod_velocity_collection = np.array(
+        [
+            [0.70544082, 0.05240655, -1.93283144, -0.35381074, -0.1305802],
+            [-0.15193337, -0.16199143, 0.94085659, 0.53076711, 2.15766298],
+            [-0.60888955, 0.36362709, 1.31370542, -0.7457939, -0.78005834],
+        ]
+    )
+
+    "initializing constants"
+    k = 1.0
+    nu = 0.5
+
+    "Function call"
+    _calculate_contact_forces_self_rod(
+        rod_position_collection[..., :-1],
+        rod_radius_collection,
+        rod_length_collection,
+        rod_tangent_collection,
+        rod_velocity_collection,
+        rod_external_forces,
+        k,
+        nu,
+    )
+
+    "Test values"
+    assert_allclose(
+        rod_external_forces,
+        np.array(
+            [
+                [-0.976629, -0.663611, -0.934193, -0.471303, -0.028112],
+                [-1.125742, -2.077304, -1.465054, -0.595825, -0.876467],
+                [0.204132, 0.967667, 0.062506, -0.299531, 0.776233],
+            ]
+        ),
+        rtol=tol,
+        atol=tol,
+    )