Dissipation constant fix (#81)

elastica/memory_block/memory_block_rod.py

@@ -129,17 +129,26 @@ def allocate_block_variables_in_nodes(self, systems: Sequence):
 
         # Things in nodes that are scalars
         #             0 ("mass", float64[:]),
-        map_scalar_dofs_in_rod_nodes = {"mass": 0}
+        map_scalar_dofs_in_rod_nodes = {
+            "mass": 0,
+            "dissipation_constant_for_forces": 1,
+        }
         self.scalar_dofs_in_rod_nodes = np.zeros(
             (len(map_scalar_dofs_in_rod_nodes), self.n_nodes)
         )
-        self.mass = self.scalar_dofs_in_rod_nodes[0]
-        for system_idx, system in enumerate(systems):
-            start_idx = self.start_idx_in_rod_nodes[system_idx]
-            end_idx = self.end_idx_in_rod_nodes[system_idx]
-            self.mass[start_idx:end_idx] = system.mass.copy()
-            # create a view back to the rod after copying variable into the block structure
-            system.mass = np.ndarray.view(self.mass[start_idx:end_idx])
+        for k, v in map_scalar_dofs_in_rod_nodes.items():
+            self.__dict__[k] = np.lib.stride_tricks.as_strided(
+                self.scalar_dofs_in_rod_nodes[v], (self.n_nodes,)
+            )
+
+        for k, v in map_scalar_dofs_in_rod_nodes.items():
+            for system_idx, system in enumerate(systems):
+                start_idx = self.start_idx_in_rod_nodes[system_idx]
+                end_idx = self.end_idx_in_rod_nodes[system_idx]
+                self.__dict__[k][..., start_idx:end_idx] = system.__dict__[k].copy()
+                system.__dict__[k] = np.ndarray.view(
+                    self.__dict__[k][..., start_idx:end_idx]
+                )
 
         # Things in nodes that are vectors
         #             0 ("position_collection", float64[:, :]),
@@ -206,8 +215,7 @@ def allocate_block_variables_in_elements(self, systems: Sequence):
             "rest_lengths": 4,
             "dilatation": 5,
             "dilatation_rate": 6,
-            "dissipation_constant_for_forces": 7,
-            "dissipation_constant_for_torques": 8,
+            "dissipation_constant_for_torques": 7,
         }
         self.scalar_dofs_in_rod_elems = np.zeros(
             (len(map_scalar_dofs_in_rod_elems), self.n_elems)

elastica/rod/cosserat_rod.py

@@ -1,5 +1,6 @@
 __doc__ = """ Rod classes and implementation details """
 __all__ = ["CosseratRod"]
+
 import typing
 
 import numpy as np
@@ -21,7 +22,6 @@
     _difference,
     _average,
 )
-from elastica.interaction import node_to_element_pos_or_vel
 
 position_difference_kernel = _difference
 position_average = _average
@@ -158,7 +158,6 @@ def __init__(
         damping_forces,
         damping_torques,
     ):
-
         self.n_elems = n_elements
         self.position_collection = position
         self.velocity_collection = velocity
@@ -728,31 +727,20 @@ def _compute_damping_forces(
     damping_forces,
     velocity_collection,
     dissipation_constant_for_forces,
-    lengths,
-    ghost_elems_idx,
 ):
     """
-    Update <damping force> given <velocity and length>
+    Update <damping force> given <velocity>
     """
 
-    # Internal damping foces.
-    elemental_velocities = node_to_element_pos_or_vel(velocity_collection)
-
-    blocksize = elemental_velocities.shape[1]
-    elemental_damping_forces = np.zeros((3, blocksize))
+    # Internal damping forces.
+    blocksize = velocity_collection.shape[1]
 
     for i in range(3):
         for k in range(blocksize):
-            elemental_damping_forces[i, k] = (
-                dissipation_constant_for_forces[k]
-                * elemental_velocities[i, k]
-                * lengths[k]
+            damping_forces[i, k] = (
+                dissipation_constant_for_forces[k] * velocity_collection[i, k]
             )
 
-    damping_forces[:] = quadrature_kernel_for_block_structure(
-        elemental_damping_forces, ghost_elems_idx
-    )
-
 
 @numba.njit(cache=True)
 def _compute_internal_forces(
@@ -818,8 +806,6 @@ def _compute_internal_forces(
         damping_forces,
         velocity_collection,
         dissipation_constant_for_forces,
-        lengths,
-        ghost_elems_idx,
     )
 
     internal_forces[:] = (
@@ -830,18 +816,16 @@ def _compute_internal_forces(
 
 @numba.njit(cache=True)
 def _compute_damping_torques(
-    damping_torques, omega_collection, dissipation_constant_for_torques, lengths
+    damping_torques, omega_collection, dissipation_constant_for_torques
 ):
     """
-    Update <damping torque> given <angular velocity and length>.
+    Update <damping torque> given <angular velocity>.
     """
-    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]
-                * omega_collection[i, k]
-                * lengths[k]
+                dissipation_constant_for_torques[k] * omega_collection[i, k]
             )
 
 
@@ -925,7 +909,7 @@ def _compute_internal_torques(
     unsteady_dilatation = J_omega_upon_e * dilatation_rate / dilatation
 
     _compute_damping_torques(
-        damping_torques, omega_collection, dissipation_constant_for_torques, lengths
+        damping_torques, omega_collection, dissipation_constant_for_torques
     )
 
     blocksize = internal_torques.shape[1]

elastica/rod/factory_function.py

@@ -192,21 +192,27 @@ def allocate(
     mass[1:] += 0.5 * density * volume
 
     # Set dissipation constant or nu array
-    dissipation_constant_for_forces = np.zeros((n_elements))
+    dissipation_constant_for_forces = np.zeros((n_elements + 1))
     # Check if the user input nu is valid
     nu_temp = np.array(nu)
     _assert_dim(nu_temp, 2, "dissipation constant (nu) for forces)")
-    dissipation_constant_for_forces[:] = nu
+    dissipation_constant_for_forces[:] = nu * mass
     # Check if the elements of dissipation constant greater than tolerance
     assert np.all(
         dissipation_constant_for_forces >= 0.0
     ), " Dissipation constant(nu) has to be equal or greater than 0."
 
     # 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.copy()
+        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)"
     )

examples/AxialStretchingCase/axial_stretching.py

@@ -55,9 +55,9 @@ class StretchingBeamSimulator(BaseSystemCollection, Constraints, Forcing, CallBa
 normal = np.array([0.0, 1.0, 0.0])
 base_length = 1.0
 base_radius = 0.025
-base_area = np.pi * base_radius**2
+base_area = np.pi * base_radius ** 2
 density = 1000
-nu = 2.0
+nu = 1.0
 youngs_modulus = 1e4
 # For shear modulus of 1e4, nu is 99!
 poisson_ratio = 0.5

examples/Binder/1_Timoshenko_Beam.ipynb

@@ -81,7 +81,7 @@
     "n_elem = 100\n",
     "\n",
     "density = 1000\n",
-    "nu = 0.1\n",
+    "nu = 1e-4\n",
     "E = 1e6\n",
     "# For shear modulus of 1e4, nu is 99!\n",
     "poisson_ratio = 99\n",

examples/Binder/2_Slithering_Snake.ipynb

@@ -69,7 +69,7 @@
     "base_radius = base_length * 0.011\n",
     "base_area = np.pi * base_radius ** 2\n",
     "density = 1000\n",
-    "nu = 1e-4\n",
+    "nu = 2e-3\n",
     "E = 1e6\n",
     "poisson_ratio = 0.5\n",
     "shear_modulus = E / (poisson_ratio + 1.0)\n",

examples/ContinuumFlagellaCase/continuum_flagella.py

@@ -33,7 +33,7 @@ def run_flagella(
     base_length = 1.0
     base_radius = 0.025
     density = 1000
-    nu = 5.0
+    nu = 2.5
     E = 1e7
     poisson_ratio = 0.5
     shear_modulus = E / (poisson_ratio + 1.0)

examples/ContinuumSnakeCase/continuum_snake.py

@@ -41,7 +41,7 @@ def run_snake(
     base_length = 0.35
     base_radius = base_length * 0.011
     density = 1000
-    nu = 1e-4
+    nu = 2e-3
     E = 1e6
     poisson_ratio = 0.5
     shear_modulus = E / (poisson_ratio + 1.0)

examples/ExperimentalCases/ParallelConnectionExample/parallel_connection_example.py

@@ -36,7 +36,7 @@ class ParallelConnection(
 base_radius = 0.007
 base_area = np.pi * base_radius ** 2
 density = 1750
-nu = 1e-2
+nu = 4e-2
 E = 3e4
 poisson_ratio = 0.5
 shear_modulus = E / (poisson_ratio + 1.0)

examples/HelicalBucklingCase/convergence_helicalbuckling.py

@@ -39,7 +39,7 @@ def simulate_helicalbucklin_beam_with(
     base_radius = 0.35
     base_area = np.pi * base_radius ** 2
     density = 1.0 / (base_area)
-    nu = 0.01
+    nu = 0.01 / density / base_area
     E = 1e6
     slack = 3
     number_of_rotations = 27
@@ -60,12 +60,11 @@ def simulate_helicalbucklin_beam_with(
         density,
         nu,
         E,
-        poisson_ratio,
     )
     # TODO: CosseratRod has to be able to take shear matrix as input, we should change it as done below
 
-    shearable_rod.shear_matrix = shear_matrix
-    shearable_rod.bend_matrix = bend_matrix
+    shearable_rod.shear_matrix[:] = shear_matrix
+    shearable_rod.bend_matrix[:] = bend_matrix
 
     helicalbuckling_sim.append(shearable_rod)
     helicalbuckling_sim.constrain(shearable_rod).using(

examples/HelicalBucklingCase/helicalbuckling.py

@@ -32,7 +32,7 @@ class HelicalBucklingSimulator(BaseSystemCollection, Constraints, Forcing):
 base_radius = 0.35
 base_area = np.pi * base_radius ** 2
 density = 1.0 / (base_area)
-nu = 0.01
+nu = 0.01 / density / base_area
 E = 1e6
 slack = 3
 number_of_rotations = 27

examples/JointCases/fixed_joint.py

@@ -34,7 +34,7 @@ class FixedJointSimulator(
 base_radius = 0.007
 base_area = np.pi * base_radius ** 2
 density = 1750
-nu = 1e-1
+nu = 0.4
 E = 3e7
 poisson_ratio = 0.5
 shear_modulus = E / (poisson_ratio + 1.0)
@@ -137,7 +137,7 @@ def make_callback(self, system, time, current_step: int):
 
 PLOT_FIGURE = True
 SAVE_FIGURE = False
-PLOT_VIDEO = False
+PLOT_VIDEO = True
 
 # plotting results
 if PLOT_FIGURE:

examples/JointCases/hinge_joint.py

@@ -34,7 +34,7 @@ class HingeJointSimulator(
 base_radius = 0.007
 base_area = np.pi * base_radius ** 2
 density = 1750
-nu = 0.001
+nu = 4e-3
 E = 3e7
 poisson_ratio = 0.5
 shear_modulus = E / (poisson_ratio + 1.0)
@@ -140,7 +140,7 @@ def make_callback(self, system, time, current_step: int):
 
 PLOT_FIGURE = True
 SAVE_FIGURE = False
-PLOT_VIDEO = False
+PLOT_VIDEO = True
 
 # plotting results
 if PLOT_FIGURE:

examples/JointCases/spherical_joint.py

@@ -35,7 +35,7 @@ class SphericalJointSimulator(
 base_radius = 0.007
 base_area = np.pi * base_radius ** 2
 density = 1750
-nu = 1e-3
+nu = 4e-3
 E = 3e7
 poisson_ratio = 0.5
 shear_modulus = E / (poisson_ratio + 1.0)

examples/MuscularFlagella/muscular_flagella.py

@@ -40,6 +40,7 @@ class MuscularFlagellaSimulator(
 base_length_body = 1.927  # mm
 E = 3.86e6  # MPa
 poisson_ratio = 0.5
+shear_modulus = E / (poisson_ratio + 1.0)
 
 base_radius_head = 0.02  # mm
 base_radius_tail = 0.007  # mm
@@ -67,7 +68,7 @@ class MuscularFlagellaSimulator(
     density_body,
     nu_body,
     E,
-    poisson_ratio,
+    shear_modulus=shear_modulus,
 )
 
 # In order to match bending stiffness of the tail as given in below reference, recompute and
@@ -131,7 +132,8 @@ class MuscularFlagellaSimulator(
 base_radius_muscle = 0.01  # mm
 base_length_muscle = 0.10756
 E_muscle = 0.3e5  # MPa
-nu_muscle = 1e-6
+shear_modulus_muscle = E_muscle / (poisson_ratio + 1.0)
+nu_muscle = 1e-6 / density_muscle / (np.pi * base_radius_muscle ** 2)
 
 # Start position of the muscle is the 4th element position of body. Lets use the exact location, because this will
 # simplify the connection implementation.
@@ -152,7 +154,7 @@ class MuscularFlagellaSimulator(
     density_muscle,
     nu_muscle,
     E_muscle,
-    poisson_ratio,
+    shear_modulus=shear_modulus_muscle,
 )
 
 muscular_flagella_sim.append(flagella_muscle)

examples/MuscularSnake/muscular_snake.py

@@ -55,7 +55,7 @@ class MuscularSnakeSimulator(
     base_length_body,
     base_radius_body,
     density_body,
-    nu,
+    nu / density_body / (np.pi * base_radius_body ** 2),
     youngs_modulus=E,
     shear_modulus=shear_modulus,
 )
@@ -101,6 +101,7 @@ class MuscularSnakeSimulator(
 muscle_radius = np.zeros((n_elem_muscle_group_one_to_three))
 muscle_radius[:] = 0.003  # First set tendon radius for whole rod.
 muscle_radius[4 * 3 : 9 * 3] = 0.006  # Change the radius of muscle elements
+nu_muscle /= density_muscle * np.pi * 0.003 ** 2
 
 for i in range(int(n_muscle_fibers / 2)):
 

examples/RodContactCase/RodRodContact/rod_rod_contact_inclined_validation.py

@@ -1,4 +1,5 @@
 import sys
+import numpy as np
 
 sys.path.append("../../../")
 from elastica import *
@@ -29,7 +30,7 @@ class InclinedRodRodContact(
 base_radius = 0.01
 base_area = np.pi * base_radius ** 2
 density = 1750
-nu = 0.001
+nu = 0.002
 E = 3e5
 poisson_ratio = 0.5
 shear_modulus = E / (poisson_ratio + 1.0)

examples/RodContactCase/RodRodContact/rod_rod_contact_parallel_validation.py

@@ -1,4 +1,5 @@
 import sys
+import numpy as np
 
 sys.path.append("../../../")
 from elastica import *
@@ -29,7 +30,7 @@ class ParallelRodRodContact(
 base_radius = 0.01
 base_area = np.pi * base_radius ** 2
 density = 1750
-nu = 0.001
+nu = 0.002
 E = 3e5
 poisson_ratio = 0.5
 shear_modulus = E / (poisson_ratio + 1.0)