From 52a0bccc64f9d035982c3e25165d607d577009b6 Mon Sep 17 00:00:00 2001 From: armantekinalp Date: Sun, 21 Nov 2021 17:00:15 -0600 Subject: [PATCH 1/8] update:np.float to np.float64 numpy depreciation warnining We changed the np.float in the code to np.float64 since numpy will stop using np.float in future relases. --- elastica/_elastica_numpy/_external_forces.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/elastica/_elastica_numpy/_external_forces.py b/elastica/_elastica_numpy/_external_forces.py index 18f118beb..b480eb07e 100644 --- a/elastica/_elastica_numpy/_external_forces.py +++ b/elastica/_elastica_numpy/_external_forces.py @@ -22,7 +22,7 @@ def __init__(self): """ pass - def apply_forces(self, system, time: np.float = 0.0): + def apply_forces(self, system, time: np.float64 = 0.0): """Apply forces to a rod-like object. In NoForces class, this routine simply passes. @@ -42,7 +42,7 @@ def apply_forces(self, system, time: np.float = 0.0): pass - def apply_torques(self, system, time: np.float = 0.0): + def apply_torques(self, system, time: np.float64 = 0.0): """Apply torques to a rod-like object. In NoForces class, this routine simply passes. @@ -156,7 +156,7 @@ def __init__(self, torque, direction=np.array([0.0, 0.0, 0.0])): super(UniformTorques, self).__init__() self.torque = (torque * direction).reshape(3, 1) - def apply_torques(self, system, time: np.float = 0.0): + def apply_torques(self, system, time: np.float64 = 0.0): torque_on_one_element = self.torque / system.n_elems system.external_torques += _batch_matvec( system.director_collection, torque_on_one_element @@ -187,7 +187,7 @@ def __init__(self, force, direction=np.array([0.0, 0.0, 0.0])): super(UniformForces, self).__init__() self.force = (force * direction).reshape(3, 1) - def apply_forces(self, system, time: np.float = 0.0): + def apply_forces(self, system, time: np.float64 = 0.0): force_on_one_element = self.force / system.n_elems system.external_forces += force_on_one_element @@ -299,7 +299,7 @@ def constant_function(input): self.my_spline = constant_function(self.s) - def apply_torques(self, system, time: np.float = 0.0): + def apply_torques(self, system, time: np.float64 = 0.0): # Ramp up the muscle torque factor = min(1.0, time / self.ramp_up_time) From 90d94a1d45c09916f9a9ac691d6d42fe392349f4 Mon Sep 17 00:00:00 2001 From: armantekinalp Date: Sun, 21 Nov 2021 17:40:47 -0600 Subject: [PATCH 2/8] rename a class to fix test function warnings This commit fixes some of the Elastica test function warnings. Previously we named a dummy rod as a TestRod, which was assumed by pytests as a test function and a warninig is raised. We rename the class as DummyTestRod. --- .../test_boundary_conditions_nb.py | 8 ++--- .../test_interaction_nb.py | 6 ++-- tests/test_elastica_numba/test_rod_nb.py | 2 +- .../test_boundary_conditions_np.py | 8 ++--- .../test_interaction_np.py | 4 +-- tests/test_elastica_numpy/test_rod_np.py | 2 +- tests/test_joint.py | 33 ++++++++++--------- tests/test_rod.py | 2 +- 8 files changed, 34 insertions(+), 31 deletions(-) diff --git a/tests/test_elastica_numba/test_boundary_conditions_nb.py b/tests/test_elastica_numba/test_boundary_conditions_nb.py index 0e2a9e78f..2bd2b5755 100644 --- a/tests/test_elastica_numba/test_boundary_conditions_nb.py +++ b/tests/test_elastica_numba/test_boundary_conditions_nb.py @@ -5,7 +5,7 @@ # System imports import numpy as np -from test_rod_nb import TestRod +from test_rod_nb import DummyTestRod from elastica._elastica_numba._boundary_conditions import ( FreeRod, OneEndFixedRod, @@ -18,7 +18,7 @@ # tests free rod boundary conditions def test_free_rod(): - test_rod = TestRod() + test_rod = DummyTestRod() free_rod = FreeRod() test_position_collection = np.random.rand(3, 20) test_rod.position_collection = ( @@ -55,7 +55,7 @@ def test_free_rod(): def test_one_end_fixed_rod(): - test_rod = TestRod() + test_rod = DummyTestRod() start_position_collection = np.random.rand(3) start_director_collection = np.random.rand(3, 3) fixed_rod = OneEndFixedRod(start_position_collection, start_director_collection) @@ -106,7 +106,7 @@ def test_helical_buckling_bc(): end_position_collection = np.array([100.0, 0.0, 0.0]) end_director_collection = np.identity(3, float) - test_rod = TestRod() + test_rod = DummyTestRod() test_position_collection = np.random.rand(3, 20) test_position_collection[..., 0] = start_position_collection diff --git a/tests/test_elastica_numba/test_interaction_nb.py b/tests/test_elastica_numba/test_interaction_nb.py index f3edcef9a..dced86239 100644 --- a/tests/test_elastica_numba/test_interaction_nb.py +++ b/tests/test_elastica_numba/test_interaction_nb.py @@ -12,11 +12,11 @@ SlenderBodyTheory, ) -# from test_rod import TestRod -from test_rod_nb import TestRod +# from test_rod import DummyTestRod +from test_rod_nb import DummyTestRod -class BaseRodClass(TestRod): +class BaseRodClass(DummyTestRod): def __init__(self, n_elem): """ This class initialize a straight rod, diff --git a/tests/test_elastica_numba/test_rod_nb.py b/tests/test_elastica_numba/test_rod_nb.py index 567ab3a90..147a58409 100644 --- a/tests/test_elastica_numba/test_rod_nb.py +++ b/tests/test_elastica_numba/test_rod_nb.py @@ -8,7 +8,7 @@ from elastica.utils import MaxDimension -class TestRod: +class DummyTestRod: def __init__(self): bs = 32 self.position_collection = np.random.randn(MaxDimension.value(), bs) diff --git a/tests/test_elastica_numpy/test_boundary_conditions_np.py b/tests/test_elastica_numpy/test_boundary_conditions_np.py index 1a243c13a..90a01e32c 100644 --- a/tests/test_elastica_numpy/test_boundary_conditions_np.py +++ b/tests/test_elastica_numpy/test_boundary_conditions_np.py @@ -5,7 +5,7 @@ # System imports import numpy as np -from test_rod_np import TestRod +from test_rod_np import DummyTestRod from elastica._elastica_numpy._boundary_conditions import ( FreeRod, OneEndFixedRod, @@ -18,7 +18,7 @@ # tests free rod boundary conditions def test_free_rod(): - test_rod = TestRod() + test_rod = DummyTestRod() free_rod = FreeRod() test_position_collection = np.random.rand(3, 20) test_rod.position_collection = ( @@ -55,7 +55,7 @@ def test_free_rod(): def test_one_end_fixed_rod(): - test_rod = TestRod() + test_rod = DummyTestRod() start_position_collection = np.random.rand(3) start_director_collection = np.random.rand(3, 3) fixed_rod = OneEndFixedRod(start_position_collection, start_director_collection) @@ -106,7 +106,7 @@ def test_helical_buckling_bc(): end_position_collection = np.array([100.0, 0.0, 0.0]) end_director_collection = np.identity(3, float) - test_rod = TestRod() + test_rod = DummyTestRod() test_position_collection = np.random.rand(3, 20) test_position_collection[..., 0] = start_position_collection diff --git a/tests/test_elastica_numpy/test_interaction_np.py b/tests/test_elastica_numpy/test_interaction_np.py index bb811736f..757d62872 100644 --- a/tests/test_elastica_numpy/test_interaction_np.py +++ b/tests/test_elastica_numpy/test_interaction_np.py @@ -13,10 +13,10 @@ # node_to_element_pos_or_vel, SlenderBodyTheory, ) -from test_rod_np import TestRod +from test_rod_np import DummyTestRod -class BaseRodClass(TestRod): +class BaseRodClass(DummyTestRod): def __init__(self, n_elem): """ This class initialize a straight rod, diff --git a/tests/test_elastica_numpy/test_rod_np.py b/tests/test_elastica_numpy/test_rod_np.py index 16004690d..58df41d4a 100644 --- a/tests/test_elastica_numpy/test_rod_np.py +++ b/tests/test_elastica_numpy/test_rod_np.py @@ -8,7 +8,7 @@ from elastica.utils import MaxDimension -class TestRod: +class DummyTestRod: def __init__(self): bs = 32 self.position_collection = np.random.randn(MaxDimension.value(), bs) diff --git a/tests/test_joint.py b/tests/test_joint.py index e50bb54be..66e37b122 100644 --- a/tests/test_joint.py +++ b/tests/test_joint.py @@ -23,9 +23,10 @@ def test_freejoint(): nu = 0.1 # Youngs Modulus [Pa] - E = 1e6 + youngs_modulus = 1e6 # poisson ratio poisson_ratio = 0.5 + shear_modulus = youngs_modulus / (poisson_ratio + 1.0) # Origin of the rod origin1 = np.array([0.0, 0.0, 0.0]) @@ -44,8 +45,8 @@ def test_freejoint(): base_radius, density, nu, - E, - poisson_ratio, + youngs_modulus, + shear_modulus=shear_modulus, ) rod2 = CosseratRod.straight_rod( n, @@ -56,8 +57,8 @@ def test_freejoint(): base_radius, density, nu, - E, - poisson_ratio, + youngs_modulus, + shear_modulus=shear_modulus, ) # Stiffness between points @@ -120,9 +121,10 @@ def test_hingejoint(): nu = 0.1 # Youngs Modulus [Pa] - E = 1e6 + youngs_modulus = 1e6 # poisson ratio poisson_ratio = 0.5 + shear_modulus = youngs_modulus / (poisson_ratio + 1.0) # Origin of the rod origin1 = np.array([0.0, 0.0, 0.0]) @@ -141,8 +143,8 @@ def test_hingejoint(): base_radius, density, nu, - E, - poisson_ratio, + youngs_modulus, + shear_modulus=shear_modulus, ) rod2 = CosseratRod.straight_rod( n, @@ -153,8 +155,8 @@ def test_hingejoint(): base_radius, density, nu, - E, - poisson_ratio, + youngs_modulus, + shear_modulus=shear_modulus, ) # Rod velocity @@ -235,9 +237,10 @@ def test_fixedjoint(): nu = 0.1 # Youngs Modulus [Pa] - E = 1e6 + youngs_modulus = 1e6 # poisson ratio poisson_ratio = 0.5 + shear_modulus = youngs_modulus / (poisson_ratio + 1.0) # Origin of the rod origin1 = np.array([0.0, 0.0, 0.0]) @@ -256,8 +259,8 @@ def test_fixedjoint(): base_radius, density, nu, - E, - poisson_ratio, + youngs_modulus, + shear_modulus=shear_modulus, ) rod2 = CosseratRod.straight_rod( n, @@ -268,8 +271,8 @@ def test_fixedjoint(): base_radius, density, nu, - E, - poisson_ratio, + youngs_modulus, + shear_modulus=shear_modulus, ) # Rod velocity diff --git a/tests/test_rod.py b/tests/test_rod.py index 2d12783f6..a8f065823 100644 --- a/tests/test_rod.py +++ b/tests/test_rod.py @@ -5,7 +5,7 @@ from elastica.utils import MaxDimension -class TestRod: +class DummyTestRod: def __init__(self): bs = 32 self.position_collection = np.random.randn(MaxDimension.value(), bs) From c44951fd860752c145f73ff836f8c184251980e6 Mon Sep 17 00:00:00 2001 From: armantekinalp Date: Sun, 21 Nov 2021 17:41:55 -0600 Subject: [PATCH 3/8] fix warnings raised due to np.float depreciation warning --- elastica/_elastica_numba/_external_forces.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/elastica/_elastica_numba/_external_forces.py b/elastica/_elastica_numba/_external_forces.py index 7a08d3e6e..9216711f1 100644 --- a/elastica/_elastica_numba/_external_forces.py +++ b/elastica/_elastica_numba/_external_forces.py @@ -26,7 +26,7 @@ def __init__(self): """ pass - def apply_forces(self, system, time: np.float = 0.0): + def apply_forces(self, system, time: np.float64 = 0.0): """Apply forces to a rod-like object. In NoForces class, this routine simply passes. @@ -46,7 +46,7 @@ def apply_forces(self, system, time: np.float = 0.0): pass - def apply_torques(self, system, time: np.float = 0.0): + def apply_torques(self, system, time: np.float64 = 0.0): """Apply torques to a rod-like object. In NoForces class, this routine simply passes. @@ -221,7 +221,7 @@ def __init__(self, torque, direction=np.array([0.0, 0.0, 0.0])): super(UniformTorques, self).__init__() self.torque = torque * direction - def apply_torques(self, system, time: np.float = 0.0): + def apply_torques(self, system, time: np.float64 = 0.0): n_elems = system.n_elems torque_on_one_element = ( _batch_product_i_k_to_ik(self.torque, np.ones((n_elems))) / n_elems @@ -255,7 +255,7 @@ def __init__(self, force, direction=np.array([0.0, 0.0, 0.0])): super(UniformForces, self).__init__() self.force = (force * direction).reshape(3, 1) - def apply_forces(self, system, time: np.float = 0.0): + def apply_forces(self, system, time: np.float64 = 0.0): force_on_one_element = self.force / system.n_elems system.external_forces += force_on_one_element @@ -368,7 +368,7 @@ def constant_function(input): self.my_spline = constant_function(self.s) - def apply_torques(self, system, time: np.float = 0.0): + def apply_torques(self, system, time: np.float64 = 0.0): self.compute_muscle_torques( time, self.my_spline, From 2d3653b97bc211eb2b39f4c9b8325839bcaff923 Mon Sep 17 00:00:00 2001 From: armantekinalp Date: Sun, 21 Nov 2021 17:44:41 -0600 Subject: [PATCH 4/8] enhancement: add user warnings if poisson ratio is used instead of shear modulus This commit changes the elastica rod initialization such that if user does not pass shear modulus as an argument during the initialization of rod then a warning message raised. We also changed the necessary test functions according to new warning and add more test functions. --- .../_elastica_numba/_rod/_cosserat_rod.py | 4 - .../_elastica_numpy/_rod/_cosserat_rod.py | 4 - elastica/rod/factory_function.py | 45 ++- .../test_governing_equations_nb.py | 3 +- .../test_governing_equations_np.py | 2 + tests/test_rod_initialisation.py | 310 ++++++++++++++++-- 6 files changed, 325 insertions(+), 43 deletions(-) diff --git a/elastica/_elastica_numba/_rod/_cosserat_rod.py b/elastica/_elastica_numba/_rod/_cosserat_rod.py index 5a2053645..5de1de3ee 100644 --- a/elastica/_elastica_numba/_rod/_cosserat_rod.py +++ b/elastica/_elastica_numba/_rod/_cosserat_rod.py @@ -133,8 +133,6 @@ def straight_rod( density, nu, youngs_modulus, - poisson_ratio, - alpha_c=4.0 / 3.0, *args, **kwargs ): @@ -182,8 +180,6 @@ def straight_rod( density, nu, youngs_modulus, - poisson_ratio, - alpha_c=4.0 / 3.0, *args, **kwargs ) diff --git a/elastica/_elastica_numpy/_rod/_cosserat_rod.py b/elastica/_elastica_numpy/_rod/_cosserat_rod.py index a73ab2f3c..017014f08 100644 --- a/elastica/_elastica_numpy/_rod/_cosserat_rod.py +++ b/elastica/_elastica_numpy/_rod/_cosserat_rod.py @@ -106,8 +106,6 @@ def straight_rod( density, nu, youngs_modulus, - poisson_ratio, - alpha_c=4.0 / 3.0, *args, **kwargs ): @@ -155,8 +153,6 @@ def straight_rod( density, nu, youngs_modulus, - poisson_ratio, - alpha_c=4.0 / 3.0, *args, **kwargs ) diff --git a/elastica/rod/factory_function.py b/elastica/rod/factory_function.py index c48f10ffd..b8d55f5cf 100644 --- a/elastica/rod/factory_function.py +++ b/elastica/rod/factory_function.py @@ -1,10 +1,9 @@ __doc__ = """ Factory function to allocate variables for Cosserat Rod""" __all__ = ["allocate"] +import warnings import numpy as np from numpy.testing import assert_allclose - from elastica.utils import MaxDimension, Tolerance - from elastica._linalg import _batch_cross, _batch_norm, _batch_dot @@ -18,8 +17,8 @@ def allocate( density, nu, youngs_modulus, - poisson_ratio, - alpha_c=4.0 / 3.0, + # poisson_ratio, + # alpha_c=4.0 / 3.0, *args, **kwargs ): @@ -216,7 +215,43 @@ def allocate( ) # Shear/Stretch matrix - shear_modulus = youngs_modulus / (poisson_ratio + 1.0) + if kwargs.__contains__("shear_modulus"): + shear_modulus = kwargs.get("shear_modulus") + if kwargs.__contains__("poisson_ratio"): + poisson_ratio = kwargs.get("poisson_ratio") + message = ( + " Poisson ratio ( " + + str(poisson_ratio) + + " ) given in kwargs is not used. \n" + + "Since shear modulus ( " + + str(shear_modulus) + + " ) is given in kwargs. \n" + ) + warnings.warn(message, category=UserWarning) + else: + if kwargs.__contains__("poisson_ratio"): + poisson_ratio = kwargs.get("poisson_ratio") + else: + message = "Shear modulus or poisson ratio cannot be found in kwargs. Poisson ratio is taken as 0.5" + warnings.warn(message, category=UserWarning) + poisson_ratio = 0.5 + + shear_modulus = youngs_modulus / (poisson_ratio + 1.0) + + message = ( + "Shear modulus cannot be found in kwargs. \n" + "Poisson ratio " + + str(poisson_ratio) + + " is used to compute shear modulus " + + str(shear_modulus) + + ", " + "using the equation: shear_modulus = youngs_modulus / (poisson_ratio + 1.0)" + ) + + warnings.warn(message, category=UserWarning) + + alpha_c = kwargs.get("alpha_c", 4.0 / 3.0) + shear_matrix = np.zeros( (MaxDimension.value(), MaxDimension.value(), n_elements), np.float64 ) diff --git a/tests/test_elastica_numba/test_governing_equations_nb.py b/tests/test_elastica_numba/test_governing_equations_nb.py index 50f8aada0..7957b00c6 100644 --- a/tests/test_elastica_numba/test_governing_equations_nb.py +++ b/tests/test_elastica_numba/test_governing_equations_nb.py @@ -37,6 +37,7 @@ def __init__(self, n_elem, nu=0.0): self.nu = nu self.E = 1 self.poisson_ratio = 0.5 + self.shear_modulus = self.E / (self.poisson_ratio + 1.0) def constructor(n_elem, nu=0.0): @@ -52,7 +53,7 @@ def constructor(n_elem, nu=0.0): cls.density, cls.nu, cls.E, - cls.poisson_ratio, + shear_modulus=cls.shear_modulus, ) return cls, rod diff --git a/tests/test_elastica_numpy/test_governing_equations_np.py b/tests/test_elastica_numpy/test_governing_equations_np.py index 2874c5b31..3e3bef24c 100644 --- a/tests/test_elastica_numpy/test_governing_equations_np.py +++ b/tests/test_elastica_numpy/test_governing_equations_np.py @@ -22,6 +22,7 @@ def __init__(self, n_elem, nu=0.0): self.nu = nu self.E = 1 self.poisson_ratio = 0.5 + self.shear_modulus = self.E / (self.poisson_ratio + 1.0) # @pytest.fixture # (scope="function") @@ -39,6 +40,7 @@ def constructor(n_elem, nu=0.0): cls.nu, cls.E, cls.poisson_ratio, + shear_modulus=cls.shear_modulus, ) return cls, rod diff --git a/tests/test_rod_initialisation.py b/tests/test_rod_initialisation.py index 7032f2a3d..c7c6901e8 100644 --- a/tests/test_rod_initialisation.py +++ b/tests/test_rod_initialisation.py @@ -94,8 +94,7 @@ def straight_rod( density, nu, youngs_modulus, - poisson_ratio, - alpha_c=4.0 / 3.0, + # poisson_ratio, *args, **kwargs ): @@ -143,7 +142,7 @@ def straight_rod( density, nu, youngs_modulus, - poisson_ratio, + # poisson_ratio, alpha_c=4.0 / 3.0, *args, **kwargs @@ -213,6 +212,7 @@ def test_input_and_output_position_array(n_elems): correct_position[0] = np.random.randn(n_elems + 1) correct_position[1] = np.random.randn(n_elems + 1) correct_position[..., 0] = start + shear_modulus = youngs_modulus / (poisson_ratio + 1.0) mockrod = MockRodForTest.straight_rod( n_elems, start, @@ -223,7 +223,7 @@ def test_input_and_output_position_array(n_elems): density, nu, youngs_modulus, - poisson_ratio, + shear_modulus=shear_modulus, position=correct_position, ) test_position = mockrod.position_collection @@ -253,6 +253,7 @@ def test_input_and_position_array_for_different_start(n_elems): nu = 0.1 youngs_modulus = 1e6 poisson_ratio = 0.3 + shear_modulus = youngs_modulus / (poisson_ratio + 1.0) # Check if the input position vector start position is different than the user defined start position correct_position = np.random.randn(3, n_elems + 1) @@ -266,7 +267,7 @@ def test_input_and_position_array_for_different_start(n_elems): density, nu, youngs_modulus, - poisson_ratio, + shear_modulus=shear_modulus, position=correct_position, ) test_position = mockrod.position_collection @@ -291,6 +292,7 @@ def test_compute_position_array_using_user_inputs(): nu = 0.1 youngs_modulus = 1e6 poisson_ratio = 0.3 + shear_modulus = youngs_modulus / (poisson_ratio + 1.0) # Check if without input position vector, output position vector is valid mockrod = MockRodForTest.straight_rod( n_elems, @@ -302,7 +304,7 @@ def test_compute_position_array_using_user_inputs(): density, nu, youngs_modulus, - poisson_ratio, + shear_modulus=shear_modulus, ) correct_position = np.zeros((3, n_elems + 1)) correct_position[0, :] = np.array([0.0, 0.25, 0.5, 0.75, 1.0]) @@ -328,6 +330,7 @@ def test_compute_directors_matrix_using_user_inputs(n_elems): nu = 0.1 youngs_modulus = 1e6 poisson_ratio = 0.3 + shear_modulus = youngs_modulus / (poisson_ratio + 1.0) # Check directors, if we dont input any directors, computed ones should be valid correct_directors = np.zeros((MaxDimension.value(), MaxDimension.value(), n_elems)) binormal = np.cross(direction, normal) @@ -349,7 +352,7 @@ def test_compute_directors_matrix_using_user_inputs(n_elems): density, nu, youngs_modulus, - poisson_ratio, + shear_modulus=shear_modulus, ) test_directors = mockrod.director_collection assert_allclose(correct_directors, test_directors, atol=Tolerance.atol()) @@ -377,6 +380,7 @@ def test_directors_using_input_position_array(n_elems): nu = 0.1 youngs_modulus = 1e6 poisson_ratio = 0.3 + shear_modulus = youngs_modulus / (poisson_ratio + 1.0) # Check directors, give position as input and let allocate function to compute directors. input_position = np.zeros((3, n_elems + 1)) input_position[0, :] = np.linspace(start[0], start[0] + base_length, n_elems + 1) @@ -401,7 +405,7 @@ def test_directors_using_input_position_array(n_elems): density, nu, youngs_modulus, - poisson_ratio, + shear_modulus=shear_modulus, position=input_position, ) test_directors = mockrod.director_collection @@ -431,6 +435,7 @@ def test_directors_using_input_directory_array(n_elems): nu = 0.1 youngs_modulus = 1e6 poisson_ratio = 0.3 + shear_modulus = youngs_modulus / (poisson_ratio + 1.0) # Check directors, give position as input and let allocate function to compute directors. input_position = np.zeros((3, n_elems + 1)) input_position[0, :] = np.linspace(start[0], start[0] + base_length, n_elems + 1) @@ -455,7 +460,7 @@ def test_directors_using_input_directory_array(n_elems): density, nu, youngs_modulus, - poisson_ratio, + shear_modulus=shear_modulus, position=input_position, directors=correct_directors, ) @@ -482,6 +487,7 @@ def test_director_if_d3_cross_d2_notequal_to_d1(): nu = 0.1 youngs_modulus = 1e6 poisson_ratio = 0.3 + shear_modulus = youngs_modulus / (poisson_ratio + 1.0) # Check directors, give directors as input and check their validity. # Let the assertion fail by setting d3=d2 for the input director input_directors = np.zeros((MaxDimension.value(), MaxDimension.value(), n_elems)) @@ -503,7 +509,7 @@ def test_director_if_d3_cross_d2_notequal_to_d1(): density, nu, youngs_modulus, - poisson_ratio, + shear_modulus=shear_modulus, directors=input_directors, ) @@ -528,6 +534,7 @@ def test_director_if_tangent_and_d3_are_not_same(): nu = 0.1 youngs_modulus = 1e6 poisson_ratio = 0.3 + shear_modulus = youngs_modulus / (poisson_ratio + 1.0) position = np.zeros((3, n_elems + 1)) end = start + direction * base_length @@ -556,7 +563,7 @@ def test_director_if_tangent_and_d3_are_not_same(): density, nu, youngs_modulus, - poisson_ratio, + shear_modulus=shear_modulus, position=position, directors=input_directors, ) @@ -584,6 +591,7 @@ def test_compute_radius_using_base_radius(n_elems): nu = 0.1 youngs_modulus = 1e6 poisson_ratio = 0.3 + shear_modulus = youngs_modulus / (poisson_ratio + 1.0) mockrod = MockRodForTest.straight_rod( n_elems, @@ -595,7 +603,7 @@ def test_compute_radius_using_base_radius(n_elems): density, nu, youngs_modulus, - poisson_ratio, + shear_modulus=shear_modulus, ) correct_radius = base_radius * np.ones((n_elems)) test_radius = mockrod.radius @@ -624,6 +632,7 @@ def test_radius_using_user_defined_radius(n_elems): nu = 0.1 youngs_modulus = 1e6 poisson_ratio = 0.3 + shear_modulus = youngs_modulus / (poisson_ratio + 1.0) mockrod = MockRodForTest.straight_rod( n_elems, @@ -635,7 +644,7 @@ def test_radius_using_user_defined_radius(n_elems): density, nu, youngs_modulus, - poisson_ratio, + shear_modulus=shear_modulus, ) correct_radius = base_radius test_radius = mockrod.radius @@ -665,6 +674,7 @@ def test_radius_not_correct_radius_shape(n_elems): nu = 0.1 youngs_modulus = 1e6 poisson_ratio = 0.3 + shear_modulus = youngs_modulus / (poisson_ratio + 1.0) MockRodForTest.straight_rod( n_elems, start, @@ -675,10 +685,238 @@ def test_radius_not_correct_radius_shape(n_elems): density, nu, youngs_modulus, - poisson_ratio, + shear_modulus=shear_modulus, ) +@pytest.mark.parametrize("n_elems", [5, 10, 50]) +@pytest.mark.parametrize("shear_modulus", [5e3, 10e3, 50e3]) +def test_shear_matrix_for_varying_shear_modulus(n_elems, shear_modulus): + """ + This test, is checking if for user defined shear modulus and validity of shear matrix. + + Returns + ------- + + """ + start = np.array([0.0, 0.0, 0.0]) + direction = np.array([1.0, 0.0, 0.0]) + normal = np.array([0.0, 0.0, 1.0]) + base_length = 1.0 + base_radius = 0.1 + density = 1000 + nu = 0.1 + youngs_modulus = 1e6 + base_area = np.pi * base_radius ** 2 + + mockrod = MockRodForTest.straight_rod( + n_elems, + start, + direction, + normal, + base_length, + base_radius, + density, + nu, + youngs_modulus, + shear_modulus=shear_modulus, + ) + + test_shear_matrix = mockrod.shear_matrix + + correct_shear_matrix = np.zeros((3, 3)) + np.fill_diagonal( + correct_shear_matrix[:], + [ + 4 / 3 * shear_modulus * base_area, + 4 / 3 * shear_modulus * base_area, + youngs_modulus * base_area, + ], + ) + + for k in range(n_elems): + assert_allclose( + correct_shear_matrix, + test_shear_matrix[..., k], + atol=Tolerance.atol(), + ) + + +@pytest.mark.parametrize("n_elems", [5, 10, 50]) +@pytest.mark.parametrize("shear_modulus", [5e3, 10e3, 50e3]) +def test_shear_matrix_for_varying_shear_modulus_warning_message_check_if_poisson_ratio_defined( + n_elems, shear_modulus +): + """ + This test, is checking if for user defined shear modulus and validity of shear matrix, + if the poisson ratio is defined. We expect if poisson ratio and shear modulus defined together then Elastica will + raise a UserWarning message and use the user defined shear modulus. + + Returns + ------- + + """ + start = np.array([0.0, 0.0, 0.0]) + direction = np.array([1.0, 0.0, 0.0]) + normal = np.array([0.0, 0.0, 1.0]) + base_length = 1.0 + base_radius = 0.1 + density = 1000 + nu = 0.1 + youngs_modulus = 1e6 + poisson_ratio = 0.3 + base_area = np.pi * base_radius ** 2 + + with pytest.warns(UserWarning): + mockrod = MockRodForTest.straight_rod( + n_elems, + start, + direction, + normal, + base_length, + base_radius, + density, + nu, + youngs_modulus, + shear_modulus=shear_modulus, + poisson_ratio=poisson_ratio, + ) + + test_shear_matrix = mockrod.shear_matrix + + correct_shear_matrix = np.zeros((3, 3)) + np.fill_diagonal( + correct_shear_matrix[:], + [ + 4 / 3 * shear_modulus * base_area, + 4 / 3 * shear_modulus * base_area, + youngs_modulus * base_area, + ], + ) + + for k in range(n_elems): + assert_allclose( + correct_shear_matrix, + test_shear_matrix[..., k], + atol=Tolerance.atol(), + ) + + +@pytest.mark.parametrize("n_elems", [5, 10, 50]) +@pytest.mark.parametrize("poisson_ratio", [0.5, 0.3, 1.0]) +def test_shear_matrix_for_varying_poisson_ratio_warning_message_check_if_no_shear_modulus_defined( + n_elems, poisson_ratio +): + """ + This test, is checking if for user defined poisson ratio and validity of shear matrix. + We expect if shear modulus is not defined then a UserWarninig is raised and shear matrix will be computed. + + Returns + ------- + + """ + start = np.array([0.0, 0.0, 0.0]) + direction = np.array([1.0, 0.0, 0.0]) + normal = np.array([0.0, 0.0, 1.0]) + base_length = 1.0 + base_radius = 0.1 + density = 1000 + nu = 0.1 + youngs_modulus = 1e6 + base_area = np.pi * base_radius ** 2 + + with pytest.warns(UserWarning): + mockrod = MockRodForTest.straight_rod( + n_elems, + start, + direction, + normal, + base_length, + base_radius, + density, + nu, + youngs_modulus, + poisson_ratio=poisson_ratio, + ) + + test_shear_matrix = mockrod.shear_matrix + + shear_modulus = youngs_modulus / (1 + poisson_ratio) + correct_shear_matrix = np.zeros((3, 3)) + np.fill_diagonal( + correct_shear_matrix[:], + [ + 4 / 3 * shear_modulus * base_area, + 4 / 3 * shear_modulus * base_area, + youngs_modulus * base_area, + ], + ) + + for k in range(n_elems): + assert_allclose( + correct_shear_matrix, + test_shear_matrix[..., k], + atol=Tolerance.atol(), + ) + + +@pytest.mark.parametrize("n_elems", [5, 10, 50]) +def test_shear_matrix_for_no_shear_modulus_or_poisson_ratio_defined_warning_message_check( + n_elems, +): + """ + This test, checks validity of shear matrix if there is no user defined shear modulus and poisson ratio. + Then Elastica uses poisson ratio of 0.5 and computes shear matrix. We expect a UserWarning to be raised. + + Returns + ------- + + """ + start = np.array([0.0, 0.0, 0.0]) + direction = np.array([1.0, 0.0, 0.0]) + normal = np.array([0.0, 0.0, 1.0]) + base_length = 1.0 + base_radius = 0.1 + density = 1000 + nu = 0.1 + youngs_modulus = 1e6 + poisson_ratio = 0.5 + base_area = np.pi * base_radius ** 2 + + with pytest.warns(UserWarning): + mockrod = MockRodForTest.straight_rod( + n_elems, + start, + direction, + normal, + base_length, + base_radius, + density, + nu, + youngs_modulus, + ) + + test_shear_matrix = mockrod.shear_matrix + + shear_modulus = youngs_modulus / (1 + poisson_ratio) + correct_shear_matrix = np.zeros((3, 3)) + np.fill_diagonal( + correct_shear_matrix[:], + [ + 4 / 3 * shear_modulus * base_area, + 4 / 3 * shear_modulus * base_area, + youngs_modulus * base_area, + ], + ) + + for k in range(n_elems): + assert_allclose( + correct_shear_matrix, + test_shear_matrix[..., k], + atol=Tolerance.atol(), + ) + + def test_inertia_shear_bend_matrices_for_varying_radius(): """ This test, is checking if for user defined varying radius, validity @@ -698,6 +936,7 @@ def test_inertia_shear_bend_matrices_for_varying_radius(): nu = 0.1 youngs_modulus = 1e6 poisson_ratio = 0.3 + shear_modulus = youngs_modulus / (poisson_ratio + 1.0) mockrod = MockRodForTest.straight_rod( n_elems, @@ -709,7 +948,7 @@ def test_inertia_shear_bend_matrices_for_varying_radius(): density, nu, youngs_modulus, - poisson_ratio, + shear_modulus=shear_modulus, ) correct_mass_second_moment_of_inertia = np.array( @@ -807,6 +1046,7 @@ def test_constant_density(n_elems): nu = 0.1 youngs_modulus = 1e6 poisson_ratio = 0.3 + shear_modulus = youngs_modulus / (poisson_ratio + 1.0) mockrod = MockRodForTest.straight_rod( n_elems, @@ -818,7 +1058,7 @@ def test_constant_density(n_elems): density, nu, youngs_modulus, - poisson_ratio, + shear_modulus=shear_modulus, ) correct_mass = density * np.pi * base_radius ** 2 * base_length / n_elems test_mass = mockrod.mass @@ -850,6 +1090,7 @@ def test_varying_density(n_elems): nu = 0.1 youngs_modulus = 1e6 poisson_ratio = 0.3 + shear_modulus = youngs_modulus / (poisson_ratio + 1.0) mockrod = MockRodForTest.straight_rod( n_elems, @@ -861,7 +1102,7 @@ def test_varying_density(n_elems): density, nu, youngs_modulus, - poisson_ratio, + shear_modulus=shear_modulus, ) volume = np.pi * base_radius ** 2 * base_length / n_elems correct_mass = np.zeros(n_elems + 1) @@ -895,6 +1136,7 @@ def test_density_invalid_shape(n_elems): nu = 0.1 youngs_modulus = 1e6 poisson_ratio = 0.3 + shear_modulus = youngs_modulus / (poisson_ratio + 1.0) MockRodForTest.straight_rod( n_elems, start, @@ -905,7 +1147,7 @@ def test_density_invalid_shape(n_elems): density, nu, youngs_modulus, - poisson_ratio, + shear_modulus=shear_modulus, ) @@ -932,6 +1174,7 @@ def test_constant_nu_for_forces(n_elems): nu = 0.1 youngs_modulus = 1e6 poisson_ratio = 0.3 + shear_modulus = youngs_modulus / (poisson_ratio + 1.0) mockrod = MockRodForTest.straight_rod( n_elems, @@ -943,7 +1186,7 @@ def test_constant_nu_for_forces(n_elems): density, nu, youngs_modulus, - poisson_ratio, + shear_modulus=shear_modulus, ) correct_nu = nu test_nu = mockrod.dissipation_constant_for_forces @@ -973,6 +1216,7 @@ def test_varying_nu_for_forces(n_elems): nu = np.linspace(0.1, 1.0, n_elems) youngs_modulus = 1e6 poisson_ratio = 0.3 + shear_modulus = youngs_modulus / (poisson_ratio + 1.0) mockrod = MockRodForTest.straight_rod( n_elems, @@ -984,7 +1228,7 @@ def test_varying_nu_for_forces(n_elems): density, nu, youngs_modulus, - poisson_ratio, + shear_modulus=shear_modulus, ) correct_nu = nu test_nu = mockrod.dissipation_constant_for_forces @@ -1014,6 +1258,7 @@ def test_nu_for_forces_invalid_shape(n_elems): nu = np.linspace(0.1, 1.0, n_elems).reshape(1, n_elems) youngs_modulus = 1e6 poisson_ratio = 0.3 + shear_modulus = youngs_modulus / (poisson_ratio + 1.0) MockRodForTest.straight_rod( n_elems, start, @@ -1024,7 +1269,7 @@ def test_nu_for_forces_invalid_shape(n_elems): density, nu, youngs_modulus, - poisson_ratio, + shear_modulus=shear_modulus, ) @@ -1052,6 +1297,7 @@ def test_constant_nu_for_torques(n_elems): nu_for_torques = 0.1 youngs_modulus = 1e6 poisson_ratio = 0.3 + shear_modulus = youngs_modulus / (poisson_ratio + 1.0) mockrod = MockRodForTest.straight_rod( n_elems, @@ -1063,7 +1309,7 @@ def test_constant_nu_for_torques(n_elems): density, nu_for_forces, youngs_modulus, - poisson_ratio, + shear_modulus=shear_modulus, nu_for_torques=nu_for_torques, ) correct_nu = nu_for_torques @@ -1096,6 +1342,7 @@ def test_varying_nu_for_torques(n_elems): nu_for_torques = np.linspace(0.1, 1.0, n_elems) youngs_modulus = 1e6 poisson_ratio = 0.3 + shear_modulus = youngs_modulus / (poisson_ratio + 1.0) mockrod = MockRodForTest.straight_rod( n_elems, @@ -1107,7 +1354,7 @@ def test_varying_nu_for_torques(n_elems): density, nu, youngs_modulus, - poisson_ratio, + shear_modulus=shear_modulus, nu_for_torques=nu_for_torques, ) correct_nu = nu_for_torques @@ -1139,6 +1386,7 @@ def test_nu_for_torques_invalid_shape(n_elems): nu_for_torques = np.linspace(0.1, 1.0, n_elems).reshape(1, n_elems) youngs_modulus = 1e6 poisson_ratio = 0.3 + shear_modulus = youngs_modulus / (poisson_ratio + 1.0) MockRodForTest.straight_rod( n_elems, start, @@ -1149,7 +1397,7 @@ def test_nu_for_torques_invalid_shape(n_elems): density, nu, youngs_modulus, - poisson_ratio, + shear_modulus=shear_modulus, nu_for_torques=nu_for_torques, ) @@ -1178,6 +1426,7 @@ def test_constant_nu_for_torques_if_not_input(n_elems): nu = 0.2 youngs_modulus = 1e6 poisson_ratio = 0.3 + shear_modulus = youngs_modulus / (poisson_ratio + 1.0) mockrod = MockRodForTest.straight_rod( n_elems, @@ -1189,7 +1438,7 @@ def test_constant_nu_for_torques_if_not_input(n_elems): density, nu, youngs_modulus, - poisson_ratio, + shear_modulus=shear_modulus, ) correct_nu = nu test_nu = mockrod.dissipation_constant_for_torques @@ -1218,6 +1467,7 @@ def test_rest_sigma_and_kappa_user_input(n_elems): nu = 0.1 youngs_modulus = 1e6 poisson_ratio = 0.3 + shear_modulus = youngs_modulus / (poisson_ratio + 1.0) input_rest_sigma = np.random.randn(3, n_elems) input_rest_kappa = np.random.randn(3, n_elems - 1) @@ -1232,7 +1482,7 @@ def test_rest_sigma_and_kappa_user_input(n_elems): density, nu, youngs_modulus, - poisson_ratio, + shear_modulus=shear_modulus, rest_sigma=input_rest_sigma, rest_kappa=input_rest_kappa, ) @@ -1269,6 +1519,7 @@ def test_rest_sigma_and_kappa_invalid_shape(n_elems): nu = 0.1 youngs_modulus = 1e6 poisson_ratio = 0.3 + shear_modulus = youngs_modulus / (poisson_ratio + 1.0) input_rest_sigma = np.random.randn(3, n_elems).reshape(n_elems, 3) input_rest_kappa = np.random.randn(3, n_elems - 1).reshape(n_elems - 1, 3) @@ -1283,7 +1534,7 @@ def test_rest_sigma_and_kappa_invalid_shape(n_elems): density, nu, youngs_modulus, - poisson_ratio, + shear_modulus=shear_modulus, rest_sigma=input_rest_sigma, rest_kappa=input_rest_kappa, ) @@ -1312,6 +1563,7 @@ def test_validity_of_allocated(n_elems): nu = 0.1 youngs_modulus = 1e6 poisson_ratio = 0.3 + shear_modulus = youngs_modulus / (poisson_ratio + 1.0) mockrod = MockRodForTest.straight_rod( n_elems, start, @@ -1322,7 +1574,7 @@ def test_validity_of_allocated(n_elems): density, nu, youngs_modulus, - poisson_ratio, + shear_modulus=shear_modulus, ) assert_allclose(n_elems, mockrod.n_elems, atol=Tolerance.atol()) @@ -1431,7 +1683,7 @@ def test_straight_rod(n_elems): density, nu, E, - poisson_ratio, + shear_modulus=G, ) # checking origin and length of rod assert_allclose(mockrod.position_collection[..., 0], start, atol=Tolerance.atol()) From 9718e6c2618629060d37c4e07eb1380ec85cd389 Mon Sep 17 00:00:00 2001 From: armantekinalp Date: Sun, 21 Nov 2021 18:51:29 -0600 Subject: [PATCH 5/8] update: change the initialization of example cases. With this commit during the initialization step of cosserat rods, we remove the poisson ratio and use shear modulus. --- examples/Binder/1_Timoshenko_Beam.ipynb | 8 ++++---- examples/Binder/2_Slithering_Snake.ipynb | 3 ++- examples/ContinuumFlagellaCase/continuum_flagella.py | 3 ++- examples/ContinuumSnakeCase/continuum_snake.py | 5 ++--- examples/FrictionValidationCases/axial_friction.py | 3 ++- .../rolling_friction_initial_velocity.py | 3 ++- .../rolling_friction_on_inclined_plane.py | 3 ++- .../FrictionValidationCases/rolling_friction_torque.py | 3 ++- examples/HelicalBucklingCase/helicalbuckling.py | 9 +++++---- examples/JointCases/fixed_joint.py | 5 +++-- examples/JointCases/hinge_joint.py | 5 +++-- examples/JointCases/spherical_joint.py | 5 +++-- examples/TimoshenkoBeamCase/timoshenko.py | 8 +++++--- 13 files changed, 37 insertions(+), 26 deletions(-) diff --git a/examples/Binder/1_Timoshenko_Beam.ipynb b/examples/Binder/1_Timoshenko_Beam.ipynb index db1553b62..346ca6b62 100644 --- a/examples/Binder/1_Timoshenko_Beam.ipynb +++ b/examples/Binder/1_Timoshenko_Beam.ipynb @@ -85,7 +85,7 @@ "E = 1e6\n", "# For shear modulus of 1e4, nu is 99!\n", "poisson_ratio = 99\n", - "\n", + "shear_modulus = E / (poisson_ratio + 1.0)\n", "\n", "start = np.zeros((3,))\n", "direction = np.array([0.0, 0.0, 1.0])\n", @@ -118,7 +118,7 @@ " density,\n", " nu,\n", " E,\n", - " poisson_ratio,\n", + " shear_modulus=shear_modulus,\n", ")\n", "\n", "timoshenko_sim.append(shearable_rod)" @@ -197,7 +197,7 @@ " nu,\n", " E,\n", " # Unshearable rod needs G -> inf, which is achievable with a poisson ratio of -1.0\n", - " poisson_ratio=-0.85,\n", + " shear_modulus = E / (-0.85 + 1.0)\n", ")\n", "\n", "timoshenko_sim.append(unshearable_rod)\n", @@ -576,4 +576,4 @@ }, "nbformat": 4, "nbformat_minor": 2 -} +} \ No newline at end of file diff --git a/examples/Binder/2_Slithering_Snake.ipynb b/examples/Binder/2_Slithering_Snake.ipynb index 511905bba..855433810 100644 --- a/examples/Binder/2_Slithering_Snake.ipynb +++ b/examples/Binder/2_Slithering_Snake.ipynb @@ -72,6 +72,7 @@ "nu = 1e-4\n", "E = 1e6\n", "poisson_ratio = 0.5\n", + "shear_modulus = E / (poisson_ratio + 1.0)\n", "\n", "# Create rod\n", "shearable_rod = CosseratRod.straight_rod(\n", @@ -84,7 +85,7 @@ " density,\n", " nu,\n", " E,\n", - " poisson_ratio,\n", + " shear_modulus = shear_modulus,\n", ")\n", "\n", "# Add rod to the snake system\n", diff --git a/examples/ContinuumFlagellaCase/continuum_flagella.py b/examples/ContinuumFlagellaCase/continuum_flagella.py index 803d3771e..bad9bc444 100644 --- a/examples/ContinuumFlagellaCase/continuum_flagella.py +++ b/examples/ContinuumFlagellaCase/continuum_flagella.py @@ -36,6 +36,7 @@ def run_flagella( nu = 5.0 E = 1e7 poisson_ratio = 0.5 + shear_modulus = E / (poisson_ratio + 1.0) shearable_rod = CosseratRod.straight_rod( n_elem, @@ -47,7 +48,7 @@ def run_flagella( density, nu, E, - poisson_ratio, + shear_modulus=shear_modulus, ) flagella_sim.append(shearable_rod) diff --git a/examples/ContinuumSnakeCase/continuum_snake.py b/examples/ContinuumSnakeCase/continuum_snake.py index 2e5c09ce5..a984af0e1 100644 --- a/examples/ContinuumSnakeCase/continuum_snake.py +++ b/examples/ContinuumSnakeCase/continuum_snake.py @@ -5,8 +5,6 @@ sys.path.append("../../") from elastica import * -# from collections import defaultdict -from examples.SnakeFrictionCase.external_force_class_for_snake import MuscleTorques from examples.ContinuumSnakeCase.continuum_snake_postprocessing import ( plot_snake_velocity, plot_video, @@ -44,6 +42,7 @@ def run_snake( nu = 1e-4 E = 1e6 poisson_ratio = 0.5 + shear_modulus = E / (poisson_ratio + 1.0) shearable_rod = CosseratRod.straight_rod( n_elem, @@ -55,7 +54,7 @@ def run_snake( density, nu, E, - poisson_ratio, + shear_modulus=shear_modulus, ) snake_sim.append(shearable_rod) diff --git a/examples/FrictionValidationCases/axial_friction.py b/examples/FrictionValidationCases/axial_friction.py index 5e931a035..dca214821 100644 --- a/examples/FrictionValidationCases/axial_friction.py +++ b/examples/FrictionValidationCases/axial_friction.py @@ -39,6 +39,7 @@ def simulate_axial_friction_with(force=0.0): E = 1e5 # For shear modulus of 2E/3 poisson_ratio = 0.5 + shear_modulus = E / (poisson_ratio + 1.0) # Set shear matrix shear_matrix = np.repeat(1e4 * np.identity((3))[:, :, np.newaxis], n_elem, axis=2) @@ -53,7 +54,7 @@ def simulate_axial_friction_with(force=0.0): density, nu, E, - poisson_ratio, + shear_modulus=shear_modulus, ) # TODO: CosseratRod has to be able to take shear matrix as input, we should change it as done below diff --git a/examples/FrictionValidationCases/rolling_friction_initial_velocity.py b/examples/FrictionValidationCases/rolling_friction_initial_velocity.py index 866b71737..8da68522b 100644 --- a/examples/FrictionValidationCases/rolling_friction_initial_velocity.py +++ b/examples/FrictionValidationCases/rolling_friction_initial_velocity.py @@ -42,6 +42,7 @@ def simulate_rolling_friction_initial_velocity_with(IFactor=0.0): E = 1e9 # For shear modulus of 2E/3 poisson_ratio = 0.5 + shear_modulus = E / (poisson_ratio + 1.0) # Set shear matrix shear_matrix = np.repeat(1e4 * np.identity((3))[:, :, np.newaxis], n_elem, axis=2) @@ -56,7 +57,7 @@ def simulate_rolling_friction_initial_velocity_with(IFactor=0.0): density, nu, E, - poisson_ratio, + shear_modulus=shear_modulus, ) # TODO: CosseratRod has to be able to take shear matrix as input, we should change it as done below diff --git a/examples/FrictionValidationCases/rolling_friction_on_inclined_plane.py b/examples/FrictionValidationCases/rolling_friction_on_inclined_plane.py index 0a8c191a7..21f4539b8 100644 --- a/examples/FrictionValidationCases/rolling_friction_on_inclined_plane.py +++ b/examples/FrictionValidationCases/rolling_friction_on_inclined_plane.py @@ -42,6 +42,7 @@ def simulate_rolling_friction_on_inclined_plane_with(alpha_s=0.0): E = 1e9 # For shear modulus of 2E/3 poisson_ratio = 0.5 + shear_modulus = E / (poisson_ratio + 1.0) # Set shear matrix shear_matrix = np.repeat(1e4 * np.identity((3))[:, :, np.newaxis], n_elem, axis=2) @@ -56,7 +57,7 @@ def simulate_rolling_friction_on_inclined_plane_with(alpha_s=0.0): density, nu, E, - poisson_ratio, + shear_modulus=shear_modulus, ) # TODO: CosseratRod has to be able to take shear matrix as input, we should change it as done below diff --git a/examples/FrictionValidationCases/rolling_friction_torque.py b/examples/FrictionValidationCases/rolling_friction_torque.py index 5db84e50e..8e0e96486 100644 --- a/examples/FrictionValidationCases/rolling_friction_torque.py +++ b/examples/FrictionValidationCases/rolling_friction_torque.py @@ -40,6 +40,7 @@ def simulate_rolling_friction_torque_with(C_s=0.0): E = 1e9 # For shear modulus of 2E/3 poisson_ratio = 0.5 + shear_modulus = E / (poisson_ratio + 1.0) # Set shear matrix shear_matrix = np.repeat(1e4 * np.identity((3))[:, :, np.newaxis], n_elem, axis=2) @@ -54,7 +55,7 @@ def simulate_rolling_friction_torque_with(C_s=0.0): density, nu, E, - poisson_ratio, + shear_modulus=shear_modulus, ) # TODO: CosseratRod has to be able to take shear matrix as input, we should change it as done below diff --git a/examples/HelicalBucklingCase/helicalbuckling.py b/examples/HelicalBucklingCase/helicalbuckling.py index 60e326010..89b70859a 100644 --- a/examples/HelicalBucklingCase/helicalbuckling.py +++ b/examples/HelicalBucklingCase/helicalbuckling.py @@ -36,9 +36,10 @@ class HelicalBucklingSimulator(BaseSystemCollection, Constraints, Forcing): E = 1e6 slack = 3 number_of_rotations = 27 -# For shear modulus of 1e4, nu is 99! -poisson_ratio = 99 -shear_matrix = np.repeat(1e5 * np.identity((3))[:, :, np.newaxis], n_elem, axis=2) +# For shear modulus of 1e5, nu is 99! +poisson_ratio = 9 +shear_modulus = E / (poisson_ratio + 1.0) +shear_matrix = np.repeat(shear_modulus * np.identity((3))[:, :, np.newaxis], n_elem, axis=2) temp_bend_matrix = np.zeros((3, 3)) np.fill_diagonal(temp_bend_matrix, [1.345, 1.345, 0.789]) bend_matrix = np.repeat(temp_bend_matrix[:, :, np.newaxis], n_elem - 1, axis=2) @@ -53,7 +54,7 @@ class HelicalBucklingSimulator(BaseSystemCollection, Constraints, Forcing): density, nu, E, - poisson_ratio, + shear_modulus=shear_modulus, ) # TODO: CosseratRod has to be able to take shear matrix as input, we should change it as done below diff --git a/examples/JointCases/fixed_joint.py b/examples/JointCases/fixed_joint.py index 977864bd9..2a2aedfc3 100644 --- a/examples/JointCases/fixed_joint.py +++ b/examples/JointCases/fixed_joint.py @@ -37,6 +37,7 @@ class FixedJointSimulator( nu = 1e-1 E = 3e7 poisson_ratio = 0.5 +shear_modulus = E / (poisson_ratio + 1.0) start_rod_1 = np.zeros((3,)) start_rod_2 = start_rod_1 + direction * base_length @@ -52,7 +53,7 @@ class FixedJointSimulator( density, nu, E, - poisson_ratio, + shear_modulus=shear_modulus, ) fixed_joint_sim.append(rod1) # Create rod 2 @@ -66,7 +67,7 @@ class FixedJointSimulator( density, nu, E, - poisson_ratio, + shear_modulus=shear_modulus, ) fixed_joint_sim.append(rod2) diff --git a/examples/JointCases/hinge_joint.py b/examples/JointCases/hinge_joint.py index eb53d8143..0871d108c 100644 --- a/examples/JointCases/hinge_joint.py +++ b/examples/JointCases/hinge_joint.py @@ -37,6 +37,7 @@ class HingeJointSimulator( nu = 0.001 E = 3e7 poisson_ratio = 0.5 +shear_modulus = E / (poisson_ratio + 1.0) start_rod_1 = np.zeros((3,)) start_rod_2 = start_rod_1 + direction * base_length @@ -52,7 +53,7 @@ class HingeJointSimulator( density, nu, E, - poisson_ratio, + shear_modulus=shear_modulus, ) hinge_joint_sim.append(rod1) # Create rod 2 @@ -66,7 +67,7 @@ class HingeJointSimulator( density, nu, E, - poisson_ratio, + shear_modulus=shear_modulus, ) hinge_joint_sim.append(rod2) diff --git a/examples/JointCases/spherical_joint.py b/examples/JointCases/spherical_joint.py index 5b8e1357d..89e8728cb 100644 --- a/examples/JointCases/spherical_joint.py +++ b/examples/JointCases/spherical_joint.py @@ -38,6 +38,7 @@ class SphericalJointSimulator( nu = 1e-3 E = 3e7 poisson_ratio = 0.5 +shear_modulus = E / (poisson_ratio + 1.0) start_rod_1 = np.zeros((3,)) start_rod_2 = start_rod_1 + direction * base_length @@ -53,7 +54,7 @@ class SphericalJointSimulator( density, nu, E, - poisson_ratio, + shear_modulus=shear_modulus, ) spherical_joint_sim.append(rod1) # Create rod 2 @@ -67,7 +68,7 @@ class SphericalJointSimulator( density, nu, E, - poisson_ratio, + shear_modulus=shear_modulus, ) spherical_joint_sim.append(rod2) diff --git a/examples/TimoshenkoBeamCase/timoshenko.py b/examples/TimoshenkoBeamCase/timoshenko.py index 7cf7d592f..df41ad4ca 100644 --- a/examples/TimoshenkoBeamCase/timoshenko.py +++ b/examples/TimoshenkoBeamCase/timoshenko.py @@ -21,7 +21,7 @@ class TimoshenkoBeamSimulator(BaseSystemCollection, Constraints, Forcing): PLOT_FIGURE = True SAVE_FIGURE = False SAVE_RESULTS = False -ADD_UNSHEARABLE_ROD = False +ADD_UNSHEARABLE_ROD = True # setting up test params n_elem = 100 @@ -36,6 +36,7 @@ class TimoshenkoBeamSimulator(BaseSystemCollection, Constraints, Forcing): E = 1e6 # For shear modulus of 1e4, nu is 99! poisson_ratio = 99 +shear_modulus = E / (poisson_ratio + 1.0) shearable_rod = CosseratRod.straight_rod( n_elem, @@ -47,7 +48,7 @@ class TimoshenkoBeamSimulator(BaseSystemCollection, Constraints, Forcing): density, nu, E, - poisson_ratio, + shear_modulus=shear_modulus, ) timoshenko_sim.append(shearable_rod) @@ -64,6 +65,7 @@ class TimoshenkoBeamSimulator(BaseSystemCollection, Constraints, Forcing): if ADD_UNSHEARABLE_ROD: # Start into the plane unshearable_start = np.array([0.0, -1.0, 0.0]) + shear_modulus = E / (-0.7 + 1.0) unshearable_rod = CosseratRod.straight_rod( n_elem, unshearable_start, @@ -75,7 +77,7 @@ class TimoshenkoBeamSimulator(BaseSystemCollection, Constraints, Forcing): nu, E, # Unshearable rod needs G -> inf, which is achievable with -ve poisson ratio - poisson_ratio=-0.7, + shear_modulus=shear_modulus, ) timoshenko_sim.append(unshearable_rod) From 5fd7c851208f17ed26c12495e14009711752d905 Mon Sep 17 00:00:00 2001 From: armantekinalp Date: Sun, 21 Nov 2021 19:46:11 -0600 Subject: [PATCH 6/8] black formatting --- examples/Binder/1_Timoshenko_Beam.ipynb | 2 +- examples/Binder/2_Slithering_Snake.ipynb | 2 +- examples/HelicalBucklingCase/helicalbuckling.py | 4 +++- 3 files changed, 5 insertions(+), 3 deletions(-) diff --git a/examples/Binder/1_Timoshenko_Beam.ipynb b/examples/Binder/1_Timoshenko_Beam.ipynb index 346ca6b62..0c0cf2562 100644 --- a/examples/Binder/1_Timoshenko_Beam.ipynb +++ b/examples/Binder/1_Timoshenko_Beam.ipynb @@ -197,7 +197,7 @@ " nu,\n", " E,\n", " # Unshearable rod needs G -> inf, which is achievable with a poisson ratio of -1.0\n", - " shear_modulus = E / (-0.85 + 1.0)\n", + " shear_modulus=E / (-0.85 + 1.0),\n", ")\n", "\n", "timoshenko_sim.append(unshearable_rod)\n", diff --git a/examples/Binder/2_Slithering_Snake.ipynb b/examples/Binder/2_Slithering_Snake.ipynb index 855433810..406ebc4af 100644 --- a/examples/Binder/2_Slithering_Snake.ipynb +++ b/examples/Binder/2_Slithering_Snake.ipynb @@ -85,7 +85,7 @@ " density,\n", " nu,\n", " E,\n", - " shear_modulus = shear_modulus,\n", + " shear_modulus=shear_modulus,\n", ")\n", "\n", "# Add rod to the snake system\n", diff --git a/examples/HelicalBucklingCase/helicalbuckling.py b/examples/HelicalBucklingCase/helicalbuckling.py index 89b70859a..e17d23109 100644 --- a/examples/HelicalBucklingCase/helicalbuckling.py +++ b/examples/HelicalBucklingCase/helicalbuckling.py @@ -39,7 +39,9 @@ class HelicalBucklingSimulator(BaseSystemCollection, Constraints, Forcing): # For shear modulus of 1e5, nu is 99! poisson_ratio = 9 shear_modulus = E / (poisson_ratio + 1.0) -shear_matrix = np.repeat(shear_modulus * np.identity((3))[:, :, np.newaxis], n_elem, axis=2) +shear_matrix = np.repeat( + shear_modulus * np.identity((3))[:, :, np.newaxis], n_elem, axis=2 +) temp_bend_matrix = np.zeros((3, 3)) np.fill_diagonal(temp_bend_matrix, [1.345, 1.345, 0.789]) bend_matrix = np.repeat(temp_bend_matrix[:, :, np.newaxis], n_elem - 1, axis=2) From cd81c904ea2abf57ab18a56d157d4f32a5bc24a5 Mon Sep 17 00:00:00 2001 From: armantekinalp Date: Sun, 21 Nov 2021 19:47:53 -0600 Subject: [PATCH 7/8] enhancement: warning functions for shear modulus and poisson ratio moved to other function. This commit adds a new function which contains the warning functions of the shear modulus and poisson ratio, since allocate function was very complex previously. --- elastica/rod/factory_function.py | 112 +++++++++++++++++++++---------- 1 file changed, 78 insertions(+), 34 deletions(-) diff --git a/elastica/rod/factory_function.py b/elastica/rod/factory_function.py index b8d55f5cf..29b13e33b 100644 --- a/elastica/rod/factory_function.py +++ b/elastica/rod/factory_function.py @@ -215,40 +215,7 @@ def allocate( ) # Shear/Stretch matrix - if kwargs.__contains__("shear_modulus"): - shear_modulus = kwargs.get("shear_modulus") - if kwargs.__contains__("poisson_ratio"): - poisson_ratio = kwargs.get("poisson_ratio") - message = ( - " Poisson ratio ( " - + str(poisson_ratio) - + " ) given in kwargs is not used. \n" - + "Since shear modulus ( " - + str(shear_modulus) - + " ) is given in kwargs. \n" - ) - warnings.warn(message, category=UserWarning) - else: - if kwargs.__contains__("poisson_ratio"): - poisson_ratio = kwargs.get("poisson_ratio") - else: - message = "Shear modulus or poisson ratio cannot be found in kwargs. Poisson ratio is taken as 0.5" - warnings.warn(message, category=UserWarning) - poisson_ratio = 0.5 - - shear_modulus = youngs_modulus / (poisson_ratio + 1.0) - - message = ( - "Shear modulus cannot be found in kwargs. \n" - "Poisson ratio " - + str(poisson_ratio) - + " is used to compute shear modulus " - + str(shear_modulus) - + ", " - "using the equation: shear_modulus = youngs_modulus / (poisson_ratio + 1.0)" - ) - - warnings.warn(message, category=UserWarning) + shear_modulus = get_shear_modulus(youngs_modulus, kwargs) alpha_c = kwargs.get("alpha_c", 4.0 / 3.0) @@ -423,3 +390,80 @@ def allocate( damping_forces, damping_torques, ) + + +def get_shear_modulus(youngs_modulus, kwargs): + """ + From the kwargs get shear modulus, or compute it. This function contains warnining messages. + + Parameters + ---------- + youngs_modulus : float + + kwargs + + Returns + ------- + + """ + # Shear/Stretch matrix + if kwargs.__contains__("shear_modulus"): + # User set shear modulus use that. + shear_modulus = kwargs.get("shear_modulus") + + if kwargs.__contains__("shear_modulus") and kwargs.__contains__("poisson_ratio"): + # User set shear modulus and also a poisson ratio. Do not use poisson ratio and raise warning. + shear_modulus = kwargs.get("shear_modulus") + poisson_ratio = kwargs.get("poisson_ratio") + message = ( + " Poisson ratio ( " + + str(poisson_ratio) + + " ) given in kwargs is not used. \n" + + "Since shear modulus ( " + + str(shear_modulus) + + " ) is given in kwargs. \n" + ) + warnings.warn(message, category=UserWarning) + + if kwargs.__contains__("poisson_ratio") and ( + not kwargs.__contains__("shear_modulus") + ): + # User set poisson ratio but not the shear modulus. Then use poisson ratio to compute shear modulus, and + # raise a warning. + poisson_ratio = kwargs.get("poisson_ratio") + shear_modulus = youngs_modulus / (poisson_ratio + 1.0) + + message = ( + "Shear modulus cannot be found in kwargs. \n" + "Poisson ratio " + + str(poisson_ratio) + + " is used to compute shear modulus " + + str(shear_modulus) + + ", " + "using the equation: shear_modulus = youngs_modulus / (poisson_ratio + 1.0)" + ) + warnings.warn(message, category=UserWarning) + + if not ( + kwargs.__contains__("poisson_ratio") or kwargs.__contains__("shear_modulus") + ): + # If user does not set poisson ratio or shear modulus, then take poisson ratio as 0.5 and raise warning. + message = "Shear modulus or poisson ratio cannot be found in kwargs. Poisson ratio is taken as 0.5" + warnings.warn(message, category=UserWarning) + poisson_ratio = 0.5 + + shear_modulus = youngs_modulus / (poisson_ratio + 1.0) + + message = ( + "Shear modulus cannot be found in kwargs. \n" + "Poisson ratio " + + str(poisson_ratio) + + " is used to compute shear modulus " + + str(shear_modulus) + + ", " + "using the equation: shear_modulus = youngs_modulus / (poisson_ratio + 1.0)" + ) + + warnings.warn(message, category=UserWarning) + + return shear_modulus From 0e25e24c69a3268058814582ffc8e544296b9880 Mon Sep 17 00:00:00 2001 From: armantekinalp Date: Sun, 28 Nov 2021 20:15:51 -0600 Subject: [PATCH 8/8] enhancement:factor_function get_shear_modulus warnings are updated This commit updates the warning messages for `get_shear_modulus` function. --- elastica/rod/factory_function.py | 24 ++++++++++++------------ 1 file changed, 12 insertions(+), 12 deletions(-) diff --git a/elastica/rod/factory_function.py b/elastica/rod/factory_function.py index 29b13e33b..01c0316e5 100644 --- a/elastica/rod/factory_function.py +++ b/elastica/rod/factory_function.py @@ -416,12 +416,13 @@ def get_shear_modulus(youngs_modulus, kwargs): shear_modulus = kwargs.get("shear_modulus") poisson_ratio = kwargs.get("poisson_ratio") message = ( - " Poisson ratio ( " + "Both a Poisson ratio and a shear modulus are provided. " + "The Poisson ratio is only used to compute a shear modulus " + "so the provided Poisson ratio of ( " + str(poisson_ratio) - + " ) given in kwargs is not used. \n" - + "Since shear modulus ( " + + " ) is being ignored in favor of the provided shear modulus ( " + str(shear_modulus) - + " ) is given in kwargs. \n" + + " ). \n" ) warnings.warn(message, category=UserWarning) @@ -434,13 +435,14 @@ def get_shear_modulus(youngs_modulus, kwargs): shear_modulus = youngs_modulus / (poisson_ratio + 1.0) message = ( - "Shear modulus cannot be found in kwargs. \n" - "Poisson ratio " + "The given Poisson ratio of " + str(poisson_ratio) - + " is used to compute shear modulus " + + " is used only to compute a shear modulus of " + str(shear_modulus) + ", " - "using the equation: shear_modulus = youngs_modulus / (poisson_ratio + 1.0)" + "using the equation: shear_modulus = youngs_modulus / (poisson_ratio + 1.0). " + "Use of a Poisson ratio will be depreciated in a future release. " + "It is encouraged that you discontinue using a Poisson ratio and instead directly provide the shear_modulus. \n" ) warnings.warn(message, category=UserWarning) @@ -448,20 +450,18 @@ def get_shear_modulus(youngs_modulus, kwargs): kwargs.__contains__("poisson_ratio") or kwargs.__contains__("shear_modulus") ): # If user does not set poisson ratio or shear modulus, then take poisson ratio as 0.5 and raise warning. - message = "Shear modulus or poisson ratio cannot be found in kwargs. Poisson ratio is taken as 0.5" - warnings.warn(message, category=UserWarning) poisson_ratio = 0.5 shear_modulus = youngs_modulus / (poisson_ratio + 1.0) message = ( - "Shear modulus cannot be found in kwargs. \n" + "Shear modulus cannot be found in kwargs. " "Poisson ratio " + str(poisson_ratio) + " is used to compute shear modulus " + str(shear_modulus) + ", " - "using the equation: shear_modulus = youngs_modulus / (poisson_ratio + 1.0)" + "using the equation: shear_modulus = youngs_modulus / (poisson_ratio + 1.0)." ) warnings.warn(message, category=UserWarning)