From 93369e7446fbe56d937751250585a5afd86ed20e Mon Sep 17 00:00:00 2001 From: Seung Hyun Kim Date: Sun, 14 Jul 2024 17:56:36 -0500 Subject: [PATCH 1/4] hotfix: twist angle for parallel elements --- elastica/rod/knot_theory.py | 9 ++++++++- 1 file changed, 8 insertions(+), 1 deletion(-) diff --git a/elastica/rod/knot_theory.py b/elastica/rod/knot_theory.py index e9688423..e6fc27c8 100644 --- a/elastica/rod/knot_theory.py +++ b/elastica/rod/knot_theory.py @@ -259,7 +259,14 @@ def _compute_twist(center_line, normal_collection): # Check if there is any Nan. Nan's appear when rod tangents are parallel to each other. idx = np.where(np.isnan(twist_temp))[0] - twist_temp[idx] = 0.0 + for i in idx: + # compute angle between two vectors + angle = np.arccos( + (normal_collection[k, 0, i] * normal_collection[k, 0, i + 1]) + + (normal_collection[k, 1, i] * normal_collection[k, 1, i + 1]) + + (normal_collection[k, 2, i] * normal_collection[k, 2, i + 1]) + ) + twist_temp[i] = angle / (2 * np.pi) local_twist[k, :] = twist_temp total_twist[k] = np.sum(twist_temp) From 8de13976194a5a1a6a46f2fcd03459b067ac7f6f Mon Sep 17 00:00:00 2001 From: Seung Hyun Kim Date: Sun, 14 Jul 2024 20:12:21 -0500 Subject: [PATCH 2/4] Update elastica/rod/knot_theory.py Co-authored-by: Arman Tekinalp <53585636+armantekinalp@users.noreply.github.com> --- elastica/rod/knot_theory.py | 16 +++++++++++----- 1 file changed, 11 insertions(+), 5 deletions(-) diff --git a/elastica/rod/knot_theory.py b/elastica/rod/knot_theory.py index e6fc27c8..c9dea94d 100644 --- a/elastica/rod/knot_theory.py +++ b/elastica/rod/knot_theory.py @@ -257,17 +257,23 @@ def _compute_twist(center_line, normal_collection): idx = np.where(twist_temp < -0.5)[0] twist_temp[idx] += 1 + # Check if there is any Nan. Nan's appear when rod tangents are parallel to each other. + idx = np.where(np.isnan(twist_temp))[0] # Check if there is any Nan. Nan's appear when rod tangents are parallel to each other. idx = np.where(np.isnan(twist_temp))[0] for i in idx: # compute angle between two vectors - angle = np.arccos( - (normal_collection[k, 0, i] * normal_collection[k, 0, i + 1]) + - (normal_collection[k, 1, i] * normal_collection[k, 1, i + 1]) + - (normal_collection[k, 2, i] * normal_collection[k, 2, i + 1]) - ) + dot_product = np.dot(projection_of_normal_collection[:, i], projection_of_normal_collection[:, i+1]) + angle = np.arccos(dot_product) twist_temp[i] = angle / (2 * np.pi) + # Make sure twist is between (-1/2 to 1/2) as defined in pg 313 Klenin & Langowski 2000 + idx = np.where(twist_temp > 0.5)[0] + twist_temp[idx] -= 1 + idx = np.where(twist_temp < -0.5)[0] + twist_temp[idx] += 1 + + local_twist[k, :] = twist_temp total_twist[k] = np.sum(twist_temp) From b5cd5327d48e342a44fd6cb447d52df19b60a73d Mon Sep 17 00:00:00 2001 From: Seung Hyun Kim Date: Sun, 14 Jul 2024 20:52:55 -0500 Subject: [PATCH 3/4] add test-case for pure-twist --- tests/test_rod/test_knot_theory.py | 47 ++++++++++++++++++++++++++++++ 1 file changed, 47 insertions(+) diff --git a/tests/test_rod/test_knot_theory.py b/tests/test_rod/test_knot_theory.py index a5630d73..9b822bbe 100644 --- a/tests/test_rod/test_knot_theory.py +++ b/tests/test_rod/test_knot_theory.py @@ -265,3 +265,50 @@ def test_knot_theory_compute_additional_segment_none_case( assert_allclose(new_center_line.shape, [timesteps, 3, n_elem]) assert_allclose(beginning_direction.shape[0], timesteps) assert_allclose(end_direction.shape[0], timesteps) + + +@pytest.mark.parametrize("n_elem", [1, 2, 10, 50, 100]) +@pytest.mark.parametrize( + "diff_angle", + [ + 0.0, + np.pi * 0.1, + np.pi * 0.5, + np.pi * 0.9, + np.pi, + -np.pi * 0.1, + -np.pi * 0.5, + -np.pi * 0.9, + -np.pi, + ], +) +def test_knot_theory_compute_pure_twist(n_elem, diff_angle): + angle = np.arange(n_elem) * diff_angle + + # setting up test params + normal = np.array([0.0, 1.0, 0.0]) + base_length = 1.0 + + position_collection = np.zeros((3, n_elem + 1)) + position_collection[2, :] = np.linspace(0, base_length, n_elem + 1) + normal_collection = np.zeros((3, n_elem)) + + for i in range(n_elem): + alpha = angle[i] + rot_matrix = np.array( + [ + [np.cos(alpha), -np.sin(alpha), 0.0], + [np.sin(alpha), np.cos(alpha), 0.0], + [0.0, 0.0, 1.0], + ] + ) + normal_temp = rot_matrix @ normal + normal_temp /= np.linalg.norm(normal_temp) + + normal_collection[:, i] = normal_temp + + total_twist, local_twist = compute_twist( + position_collection[None, ...], normal_collection[None, ...] + ) + np.testing.assert_allclose(local_twist[0], np.diff(angle / (2 * np.pi))) + np.testing.assert_allclose(total_twist[0], angle[-1] / (2 * np.pi)) From 04119263ca61ddc01d1e878de7bc5df6359e7f44 Mon Sep 17 00:00:00 2001 From: Seung Hyun Kim Date: Sun, 14 Jul 2024 20:54:36 -0500 Subject: [PATCH 4/4] fix: cosine angle computed with correct sign --- elastica/rod/knot_theory.py | 18 +++++++++--------- tests/test_rod/test_knot_theory.py | 4 ++-- 2 files changed, 11 insertions(+), 11 deletions(-) diff --git a/elastica/rod/knot_theory.py b/elastica/rod/knot_theory.py index c9dea94d..1ac8f353 100644 --- a/elastica/rod/knot_theory.py +++ b/elastica/rod/knot_theory.py @@ -251,20 +251,21 @@ def _compute_twist(center_line, normal_collection): alpha /= 2 * np.pi gamma /= 2 * np.pi twist_temp = alpha + gamma - # Make sure twist is between (-1/2 to 1/2) as defined in pg 313 Klenin & Langowski 2000 - idx = np.where(twist_temp > 0.5)[0] - twist_temp[idx] -= 1 - idx = np.where(twist_temp < -0.5)[0] - twist_temp[idx] += 1 - # Check if there is any Nan. Nan's appear when rod tangents are parallel to each other. - idx = np.where(np.isnan(twist_temp))[0] # Check if there is any Nan. Nan's appear when rod tangents are parallel to each other. idx = np.where(np.isnan(twist_temp))[0] for i in idx: # compute angle between two vectors - dot_product = np.dot(projection_of_normal_collection[:, i], projection_of_normal_collection[:, i+1]) + dot_product = np.dot( + projection_of_normal_collection[:, i], + projection_of_normal_collection[:, i + 1], + ) angle = np.arccos(dot_product) + cross = np.cross( + projection_of_normal_collection[:, i], + projection_of_normal_collection[:, i + 1], + ) + angle *= np.sign(np.dot(tangent[:, i], cross)) twist_temp[i] = angle / (2 * np.pi) # Make sure twist is between (-1/2 to 1/2) as defined in pg 313 Klenin & Langowski 2000 @@ -273,7 +274,6 @@ def _compute_twist(center_line, normal_collection): idx = np.where(twist_temp < -0.5)[0] twist_temp[idx] += 1 - local_twist[k, :] = twist_temp total_twist[k] = np.sum(twist_temp) diff --git a/tests/test_rod/test_knot_theory.py b/tests/test_rod/test_knot_theory.py index 9b822bbe..6911fe0a 100644 --- a/tests/test_rod/test_knot_theory.py +++ b/tests/test_rod/test_knot_theory.py @@ -275,11 +275,11 @@ def test_knot_theory_compute_additional_segment_none_case( np.pi * 0.1, np.pi * 0.5, np.pi * 0.9, - np.pi, + np.pi * 0.999, -np.pi * 0.1, -np.pi * 0.5, -np.pi * 0.9, - -np.pi, + -np.pi * 0.999, ], ) def test_knot_theory_compute_pure_twist(n_elem, diff_angle):