diff --git a/elastica/rod/cosserat_rod.py b/elastica/rod/cosserat_rod.py index eb0a0cf1d..690e3e955 100644 --- a/elastica/rod/cosserat_rod.py +++ b/elastica/rod/cosserat_rod.py @@ -817,16 +817,16 @@ def _compute_internal_forces( @numba.njit(cache=True) def _compute_damping_torques( - damping_torques, J_omega_upon_e, dissipation_constant_for_torques + damping_torques, omega_collection, dissipation_constant_for_torques ): """ - Update given . + Update given . """ - blocksize = damping_torques.shape[1] + blocksize = omega_collection.shape[1] for i in range(3): for k in range(blocksize): damping_torques[i, k] = ( - dissipation_constant_for_torques[k] * J_omega_upon_e[i, k] + dissipation_constant_for_torques[k] * omega_collection[i, k] ) @@ -910,7 +910,7 @@ def _compute_internal_torques( unsteady_dilatation = J_omega_upon_e * dilatation_rate / dilatation _compute_damping_torques( - damping_torques, J_omega_upon_e, dissipation_constant_for_torques + damping_torques, omega_collection, dissipation_constant_for_torques ) blocksize = internal_torques.shape[1] diff --git a/elastica/rod/factory_function.py b/elastica/rod/factory_function.py index 44541400f..5d74985d6 100644 --- a/elastica/rod/factory_function.py +++ b/elastica/rod/factory_function.py @@ -204,12 +204,15 @@ def allocate( # Custom nu for torques dissipation_constant_for_torques = np.zeros((n_elements)) + elemental_mass = (mass[1:] + mass[:-1]) / 2.0 if nu_for_torques is None: dissipation_constant_for_torques[:] = (dissipation_constant_for_forces / mass)[ :n_elements - ] + ] * elemental_mass else: - dissipation_constant_for_torques[:] = np.asarray(nu_for_torques) + dissipation_constant_for_torques[:] = ( + np.asarray(nu_for_torques) * elemental_mass + ) _assert_dim( dissipation_constant_for_torques, 2, "dissipation constant (nu) for torque)" ) diff --git a/tests/test_math/test_governing_equations.py b/tests/test_math/test_governing_equations.py index 8d5b7e985..71d849296 100644 --- a/tests/test_math/test_governing_equations.py +++ b/tests/test_math/test_governing_equations.py @@ -409,15 +409,17 @@ def test_compute_damping_forces_torques(self, n_elem, nu): # Construct velocity and omega test_rod.velocity_collection[:] = 1.0 test_rod.omega_collection[:] = 1.0 - J_omega_upon_e = test_rod.omega_collection.copy() # Compute damping forces and torques damping_forces = ( np.repeat(np.array([1.0, 1.0, 1.0])[:, np.newaxis], n_elem + 1, axis=1) * nu * test_rod.mass ) + elemental_mass = (test_rod.mass[1:] + test_rod.mass[:-1]) / 2.0 damping_torques = ( - np.repeat(np.array([1.0, 1.0, 1.0])[:, np.newaxis], n_elem, axis=1) * nu + np.repeat(np.array([1.0, 1.0, 1.0])[:, np.newaxis], n_elem, axis=1) + * nu + * elemental_mass ) # Compute geometry from state _compute_geometry_from_state( @@ -435,7 +437,7 @@ def test_compute_damping_forces_torques(self, n_elem, nu): ) _compute_damping_torques( test_rod.damping_torques, - J_omega_upon_e, + test_rod.omega_collection, test_rod.dissipation_constant_for_torques, ) diff --git a/tests/test_rod_initialisation.py b/tests/test_rod_initialisation.py index d7ef12aaf..9b0d84065 100644 --- a/tests/test_rod_initialisation.py +++ b/tests/test_rod_initialisation.py @@ -1191,7 +1191,8 @@ def test_constant_nu_for_torques(n_elems): shear_modulus=shear_modulus, nu_for_torques=nu_for_torques, ) - correct_nu = nu_for_torques + elemental_mass = (mockrod.mass[1:] + mockrod.mass[:-1]) / 2.0 + correct_nu = nu_for_torques * elemental_mass test_nu = mockrod.dissipation_constant_for_torques assert_allclose(correct_nu, test_nu, atol=Tolerance.atol()) @@ -1236,7 +1237,8 @@ def test_varying_nu_for_torques(n_elems): shear_modulus=shear_modulus, nu_for_torques=nu_for_torques, ) - correct_nu = nu_for_torques + elemental_mass = (mockrod.mass[1:] + mockrod.mass[:-1]) / 2.0 + correct_nu = nu_for_torques * elemental_mass test_nu = mockrod.dissipation_constant_for_torques assert_allclose(correct_nu, test_nu, atol=Tolerance.atol()) @@ -1319,7 +1321,8 @@ def test_constant_nu_for_torques_if_not_input(n_elems): youngs_modulus, shear_modulus=shear_modulus, ) - correct_nu = nu + elemental_mass = (mockrod.mass[1:] + mockrod.mass[:-1]) / 2.0 + correct_nu = nu * elemental_mass test_nu = mockrod.dissipation_constant_for_torques assert_allclose(correct_nu, test_nu, atol=Tolerance.atol()) @@ -1593,7 +1596,12 @@ def test_straight_rod(n_elems): nu * mockrod.mass, atol=Tolerance.atol(), ) - assert_allclose(mockrod.dissipation_constant_for_torques, nu, atol=Tolerance.atol()) + elemental_mass = (mockrod.mass[1:] + mockrod.mass[:-1]) / 2.0 + assert_allclose( + mockrod.dissipation_constant_for_torques, + nu * elemental_mass, + atol=Tolerance.atol(), + ) assert_allclose( rest_voronoi_lengths, mockrod.rest_voronoi_lengths, atol=Tolerance.atol() )