Repository: GazzolaLab/PyElastica

Issue:
Typo in `compute_shear_energy` method of Cosserat Rod.
Although the code to compute shear energy function is correct, naming of one of the variables is wrong. `shear_internal_torques` is not correct and misleading, we should call this variable as `shear_internal_forces` . 

https://github.com/GazzolaLab/PyElastica/blob/fee87b9da22e310ff925c16fdc839bf8405c51a4/elastica/rod/cosserat_rod.py#L516

Discussion:

Patch:
diff --git a/Makefile b/Makefile
index 4c032cb..a162076 100644
--- a/Makefile
+++ b/Makefile
@@ -1,6 +1,6 @@
 black:
 	@black --version
-	@black --required-version 21.12b0 elastica tests examples
+	@black --required-version 21.12b0 elastica tests
 
 
 black_check:
diff --git a/RELEASE.md b/RELEASE.md
index c224808..5e74ac3 100644
--- a/RELEASE.md
+++ b/RELEASE.md
@@ -1,15 +1,3 @@
-# Release Note (version 0.2.3)
-
-
-## Developer Note
-The major updates are knot theory module added to the Cosserat rod as *mixin*, and muscular snake example is added. 
-
-## Notable Changes
-- #70: Knot theory module to compute topological quantities.
-- #71: Reorganize rod constructor warning messages and collect messages in log.  
-- #72: Muscular snake example is added.
----
-
 # Release Note (version 0.2.2)
 
 ## Developer Note
diff --git a/elastica/_elastica_numba.py b/elastica/_elastica_numba.py
new file mode 100644
index 0000000..39b4a6d
--- /dev/null
+++ b/elastica/_elastica_numba.py
@@ -0,0 +1,8 @@
+class _elastica_numba:
+    """The purpose is to throw deprecation error to people previously
+    using _elastica_numba module. Please remove this in the future.
+    """
+
+    raise ImportError(
+        "The module _elastica_numba is moved to main implementation. Please import the modules without _elastica_numba."
+    )
diff --git a/elastica/_elastica_numba/_calculus.py b/elastica/_elastica_numba/_calculus.py
deleted file mode 100644
index 1edc7cf..0000000
--- a/elastica/_elastica_numba/_calculus.py
+++ /dev/null
@@ -1,24 +0,0 @@
-import warnings
-from elastica._calculus import (
-    _trapezoidal,
-    _two_point_difference,
-    _difference,
-    _average,
-    _clip_array,
-    _isnan_check,
-    _get_zero_array,
-    _trapezoidal_for_block_structure,
-    _two_point_difference_for_block_structure,
-)
-
-warnings.warn(
-    "The numba-implementation is included in the default elastica module. Please import without _elastica_numba.",
-    DeprecationWarning,
-)
-
-position_difference_kernel = _difference
-position_average = _average
-quadrature_kernel = _trapezoidal
-difference_kernel = _two_point_difference
-quadrature_kernel_for_block_structure = _trapezoidal_for_block_structure
-difference_kernel_for_block_structure = _two_point_difference_for_block_structure
diff --git a/elastica/_elastica_numba/_external_forces.py b/elastica/_elastica_numba/_external_forces.py
deleted file mode 100644
index 5d53a12..0000000
--- a/elastica/_elastica_numba/_external_forces.py
+++ /dev/null
@@ -1,14 +0,0 @@
-import warnings
-from elastica.external_forces import (
-    NoForces,
-    GravityForces,
-    EndpointForces,
-    UniformTorques,
-    UniformForces,
-    MuscleTorques,
-)
-
-warnings.warn(
-    "The numba-implementation is included in the default elastica module. Please import without _elastica_numba.",
-    DeprecationWarning,
-)
diff --git a/elastica/_elastica_numba/_interaction.py b/elastica/_elastica_numba/_interaction.py
deleted file mode 100644
index 8865650..0000000
--- a/elastica/_elastica_numba/_interaction.py
+++ /dev/null
@@ -1,17 +0,0 @@
-import warnings
-from elastica.interaction import (
-    AnisotropicFrictionalPlane,
-    SlenderBodyTheory,
-    # AnisotropicFrictionalPlaneRigidBody,
-    # InteractionPlane,
-    # InteractionPlaneRigidBody,
-    elements_to_nodes_inplace,
-    node_to_element_pos_or_vel,
-    nodes_to_elements,
-    sum_over_elements,
-)
-
-warnings.warn(
-    "The numba-implementation is included in the default elastica module. Please import without _elastica_numba.",
-    DeprecationWarning,
-)
diff --git a/elastica/_elastica_numba/_joint.py b/elastica/_elastica_numba/_joint.py
deleted file mode 100644
index 89a692e..0000000
--- a/elastica/_elastica_numba/_joint.py
+++ /dev/null
@@ -1,14 +0,0 @@
-import warnings
-
-from elastica.joint import (
-    FreeJoint,
-    HingeJoint,
-    FixedJoint,
-    ExternalContact,
-    SelfContact,
-)
-
-warnings.warn(
-    "The numba-implementation is included in the default elastica module. Please import without _elastica_numba.",
-    DeprecationWarning,
-)
diff --git a/elastica/_elastica_numba/_linalg.py b/elastica/_elastica_numba/_linalg.py
deleted file mode 100644
index 445bba7..0000000
--- a/elastica/_elastica_numba/_linalg.py
+++ /dev/null
@@ -1,20 +0,0 @@
-import warnings
-from elastica._linalg import (
-    levi_civita_tensor,
-    _batch_matvec,
-    _batch_matmul,
-    _batch_cross,
-    _batch_vec_oneD_vec_cross,
-    _batch_dot,
-    _batch_norm,
-    _batch_product_i_k_to_ik,
-    _batch_product_i_ik_to_k,
-    _batch_product_k_ik_to_ik,
-    _batch_vector_sum,
-    _batch_matrix_transpose,
-)
-
-warnings.warn(
-    "The numba-implementation is included in the default elastica module. Please import without _elastica_numba.",
-    DeprecationWarning,
-)
diff --git a/elastica/_elastica_numba/_rod/__init__.py b/elastica/_elastica_numba/_rod/__init__.py
deleted file mode 100644
index e69de29..0000000
diff --git a/elastica/_elastica_numba/_rod/_cosserat_rod.py b/elastica/_elastica_numba/_rod/_cosserat_rod.py
deleted file mode 100644
index 1aa9d75..0000000
--- a/elastica/_elastica_numba/_rod/_cosserat_rod.py
+++ /dev/null
@@ -1,7 +0,0 @@
-import warnings
-from elastica.rod.cosserat_rod import CosseratRod
-
-warnings.warn(
-    "The numba-implementation is included in the default elastica module. Please import without _elastica_numba.",
-    DeprecationWarning,
-)
diff --git a/elastica/_elastica_numba/_rod/_data_structures.py b/elastica/_elastica_numba/_rod/_data_structures.py
deleted file mode 100644
index 0f741c4..0000000
--- a/elastica/_elastica_numba/_rod/_data_structures.py
+++ /dev/null
@@ -1,16 +0,0 @@
-import warnings
-import numpy as np
-from elastica._rotations import _get_rotation_matrix, _rotate
-from elastica.rod.data_structures import (
-    _RodSymplecticStepperMixin,
-    _bootstrap_from_data,
-    _State,
-    _DerivativeState,
-    _KinematicState,
-    _DynamicState,
-)
-
-warnings.warn(
-    "The numba-implementation is included in the default elastica module. Please import without _elastica_numba.",
-    DeprecationWarning,
-)
diff --git a/elastica/_elastica_numba/_rotations.py b/elastica/_elastica_numba/_rotations.py
deleted file mode 100644
index f3a0d7b..0000000
--- a/elastica/_elastica_numba/_rotations.py
+++ /dev/null
@@ -1,7 +0,0 @@
-import warnings
-from elastica._rotations import _get_rotation_matrix, _inv_rotate, _rotate
-
-warnings.warn(
-    "The numba-implementation is included in the default elastica module. Please import without _elastica_numba.",
-    DeprecationWarning,
-)
diff --git a/elastica/_elastica_numba/_timestepper/__init__.py b/elastica/_elastica_numba/_timestepper/__init__.py
deleted file mode 100644
index 416d02a..0000000
--- a/elastica/_elastica_numba/_timestepper/__init__.py
+++ /dev/null
@@ -1,10 +0,0 @@
-__doc__ = """Timestepping utilities to be used with Rod and RigidBody classes of Elastica Numba implementation"""
-
-import warnings
-from elastica.timestepper import extend_stepper_interface, integrate
-
-
-warnings.warn(
-    "The numba-implementation is included in the default elastica module. Please import without _elastica_numba.",
-    DeprecationWarning,
-)
diff --git a/elastica/_elastica_numba/_timestepper/_explicit_steppers.py b/elastica/_elastica_numba/_timestepper/_explicit_steppers.py
deleted file mode 100644
index 560c2bb..0000000
--- a/elastica/_elastica_numba/_timestepper/_explicit_steppers.py
+++ /dev/null
@@ -1,14 +0,0 @@
-__doc__ = """Explicit timesteppers  and concepts of Elastica Numba implementation"""
-
-import warnings
-from elastica.timestepper.explicit_steppers import (
-    ExplicitStepperTag,
-    RungeKutta4,
-    EulerForward,
-)
-
-
-warnings.warn(
-    "The numba-implementation is included in the default elastica module. Please import without _elastica_numba.",
-    DeprecationWarning,
-)
diff --git a/elastica/_elastica_numba/_timestepper/_symplectic_steppers.py b/elastica/_elastica_numba/_timestepper/_symplectic_steppers.py
deleted file mode 100644
index d6c9c36..0000000
--- a/elastica/_elastica_numba/_timestepper/_symplectic_steppers.py
+++ /dev/null
@@ -1,14 +0,0 @@
-__doc__ = """Symplectic timesteppers and concepts of Elastica Numba implementation"""
-
-import warnings
-from elastica.timestepper.symplectic_steppers import (
-    SymplecticStepperTag,
-    PositionVerlet,
-    PEFRL,
-)
-
-
-warnings.warn(
-    "The numba-implementation is included in the default elastica module. Please import without _elastica_numba.",
-    DeprecationWarning,
-)
diff --git a/elastica/rod/cosserat_rod.py b/elastica/rod/cosserat_rod.py
index fa0eb12..d74f48e 100644
--- a/elastica/rod/cosserat_rod.py
+++ b/elastica/rod/cosserat_rod.py
@@ -514,11 +514,11 @@ class CosseratRod(RodBase, KnotTheory):
         """
 
         sigma_diff = self.sigma - self.rest_sigma
-        shear_internal_torques = _batch_matvec(self.shear_matrix, sigma_diff)
+        shear_internal_forces = _batch_matvec(self.shear_matrix, sigma_diff)
 
         return (
             0.5
-            * (_batch_dot(sigma_diff, shear_internal_torques) * self.rest_lengths).sum()
+            * (_batch_dot(sigma_diff, shear_internal_forces) * self.rest_lengths).sum()
         )
 
 
diff --git a/elastica/rod/factory_function.py b/elastica/rod/factory_function.py
index 267f7a0..7cf1378 100644
--- a/elastica/rod/factory_function.py
+++ b/elastica/rod/factory_function.py
@@ -145,7 +145,7 @@ def allocate(
 
     # Value taken based on best correlation for Poisson ratio = 0.5, from
     # "On Timoshenko's correction for shear in vibrating beams" by Kaneko, 1975
-    alpha_c = 0.964
+    alpha_c = 27.0 / 28.0
     shear_matrix = np.zeros(
         (MaxDimension.value(), MaxDimension.value(), n_elements), np.float64
     )
diff --git a/elastica/version.py b/elastica/version.py
index ae50f1f..867e4b1 100644
--- a/elastica/version.py
+++ b/elastica/version.py
@@ -1,1 +1,1 @@
-VERSION = "0.2.3"
+VERSION = "0.3.0"
diff --git a/examples/ExperimentalCases/ParallelConnectionExample/parallel_connection_example.py b/examples/ExperimentalCases/ParallelConnectionExample/parallel_connection_example.py
index 5ca8d3a..08db51c 100644
--- a/examples/ExperimentalCases/ParallelConnectionExample/parallel_connection_example.py
+++ b/examples/ExperimentalCases/ParallelConnectionExample/parallel_connection_example.py
@@ -6,10 +6,7 @@ import sys
 # FIXME without appending sys.path make it more generic
 sys.path.append("../../../../")
 from elastica import *
-from elastica.experimental.connection_contact_joint.parallel_connection import (
-    get_connection_vector_straight_straight_rod,
-    SurfaceJointSideBySide,
-)
+from elastica.experimental.connection_contact_joint.parallel_connection import get_connection_vector_straight_straight_rod, SurfaceJointSideBySide
 from elastica._calculus import difference_kernel
 from examples.JointCases.joint_cases_postprocessing import (
     plot_position,
@@ -85,11 +82,7 @@ parallel_connection_sim.constrain(rod_two).using(
 
 # Apply a contraction force on rod one.
 class ContractionForce(NoForces):
-    def __init__(
-        self,
-        ramp,
-        force_mag,
-    ):
+    def __init__(self, ramp, force_mag, ):
         self.ramp = ramp
         self.force_mag = force_mag
 
@@ -97,42 +90,26 @@ class ContractionForce(NoForces):
         # Ramp the force
         factor = min(1.0, time / self.ramp)
 
-        system.external_forces[:] -= factor * difference_kernel(
-            self.force_mag * system.tangents
-        )
-
+        system.external_forces[:] -= factor * difference_kernel(self.force_mag * system.tangents)
 
-parallel_connection_sim.add_forcing_to(rod_one).using(
-    ContractionForce, ramp=0.5, force_mag=1.0
-)
+parallel_connection_sim.add_forcing_to(rod_one).using(ContractionForce, ramp=0.5, force_mag=1. )
 
 # Connect rod 1 and rod 2
 (
     rod_one_direction_vec_in_material_frame,
     rod_two_direction_vec_in_material_frame,
     offset_btw_rods,
-) = get_connection_vector_straight_straight_rod(
-    rod_one, rod_two, (0, n_elem), (0, n_elem)
-)
+) = get_connection_vector_straight_straight_rod(rod_one, rod_two, (0, n_elem), (0, n_elem))
 
 for i in range(n_elem):
     parallel_connection_sim.connect(
         first_rod=rod_one, second_rod=rod_two, first_connect_idx=i, second_connect_idx=i
     ).using(
-        SurfaceJointSideBySide,
-        k=1e2,
-        nu=1e-5,
-        k_repulsive=1e3,
-        rod_one_direction_vec_in_material_frame=rod_one_direction_vec_in_material_frame[
-            :, i
-        ],
-        rod_two_direction_vec_in_material_frame=rod_two_direction_vec_in_material_frame[
-            :, i
-        ],
-        offset_btw_rods=offset_btw_rods[i],
+        SurfaceJointSideBySide, k=1E2, nu=1E-5, k_repulsive = 1E3, rod_one_direction_vec_in_material_frame=rod_one_direction_vec_in_material_frame[:,i],
+        rod_two_direction_vec_in_material_frame = rod_two_direction_vec_in_material_frame[:,i],
+        offset_btw_rods = offset_btw_rods[i]
     )  # k=kg/s2 nu=kg/s 1e-2
 
-
 class ParallelConnecitonCallback(CallBackBaseClass):
     """
     Call back function for parallel connection
diff --git a/examples/MuscularSnake/muscle_forces.py b/examples/MuscularSnake/muscle_forces.py
index 3a59c6d..b30c4b6 100644
--- a/examples/MuscularSnake/muscle_forces.py
+++ b/examples/MuscularSnake/muscle_forces.py
@@ -37,7 +37,7 @@ class MuscleForces(NoForces):
         side_of_body,
         muscle_start_end_index,
         step,
-        post_processing,
+            post_processing,
     ):
         """
 
@@ -76,7 +76,7 @@ class MuscleForces(NoForces):
 
         self.post_processing = post_processing
         self.step = step
-        self.counter = 0
+        self.counter=0
 
     def apply_forces(self, system, time: np.float = 0.0):
         forces = self._apply_forces(
@@ -90,16 +90,15 @@ class MuscleForces(NoForces):
             system.external_forces,
         )
 
-        if self.counter % self.step == 0:
+        if self.counter % self.step ==0:
             self.post_processing["time"].append(time)
             self.post_processing["step"].append(self.counter)
             self.post_processing["external_forces"].append(forces.copy())
-            self.post_processing["element_position"].append(
-                np.cumsum(system.lengths).copy()
-            )
+            self.post_processing["element_position"].append(np.cumsum(system.lengths).copy())
 
         self.counter += 1
 
+
     @staticmethod
     @njit(cache=True)
     def _apply_forces(
@@ -131,4 +130,4 @@ class MuscleForces(NoForces):
         ] += difference_kernel(muscle_forces)
 
         external_forces += forces
-        return forces
+        return forces
\ No newline at end of file
diff --git a/examples/MuscularSnake/muscular_snake.py b/examples/MuscularSnake/muscular_snake.py
index d7495b7..9a0bfcb 100644
--- a/examples/MuscularSnake/muscular_snake.py
+++ b/examples/MuscularSnake/muscular_snake.py
@@ -260,8 +260,8 @@ for idx, rod_two in enumerate(muscle_rod_list):
     ) = get_connection_vector_straight_straight_rod(
         rod_one,
         rod_two,
-        (muscle_start_connection_index[idx], muscle_end_connection_index[idx]),
-        (0, rod_two.n_elems),
+        (muscle_start_connection_index[idx],muscle_end_connection_index[idx]),
+        (0, rod_two.n_elems)
     )
     straight_straight_rod_connection_list.append(
         [
diff --git a/examples/MuscularSnake/post_processing.py b/examples/MuscularSnake/post_processing.py
index 44b36eb..8de47eb 100644
--- a/examples/MuscularSnake/post_processing.py
+++ b/examples/MuscularSnake/post_processing.py
@@ -12,16 +12,16 @@ from typing import Dict, Sequence
 
 
 def plot_video_with_surface(
-    rods_history: Sequence[Dict],
-    video_name="video.mp4",
-    fps=60,
-    step=1,
-    vis2D=True,
-    **kwargs,
+        rods_history: Sequence[Dict],
+        video_name="video.mp4",
+        fps=60,
+        step=1,
+        vis2D=True,
+        **kwargs,
 ):
     plt.rcParams.update({"font.size": 22})
 
-    folder_name = kwargs.get("folder_name", "")
+    folder_name = kwargs.get("folder_name","")
 
     # 2d case <always 2d case for now>
     import matplotlib.animation as animation
@@ -115,7 +115,7 @@ def plot_video_with_surface(
                 ax.add_artist(sphere_artists[sphere_idx])
 
         # ax.set_aspect("equal")
-        video_name_3D = folder_name + "3D_" + video_name
+        video_name_3D = folder_name+"3D_" + video_name
 
         with writer.saving(fig, video_name_3D, dpi):
             with plt.style.context("seaborn-whitegrid"):
@@ -127,7 +127,7 @@ def plot_video_with_surface(
                         )
                         if not inst_position.shape[1] == inst_radius.shape[0]:
                             inst_position = 0.5 * (
-                                inst_position[..., 1:] + inst_position[..., :-1]
+                                    inst_position[..., 1:] + inst_position[..., :-1]
                             )
 
                         rod_scatters[rod_idx]._offsets3d = (
@@ -206,7 +206,7 @@ def plot_video_with_surface(
                 ax.add_artist(sphere_artists[sphere_idx])
 
         ax.set_aspect("equal")
-        video_name_2D = folder_name + "2D_xy_" + video_name
+        video_name_2D = folder_name+"2D_xy_" + video_name
 
         with writer.saving(fig, video_name_2D, dpi):
             with plt.style.context("seaborn-whitegrid"):
@@ -218,7 +218,7 @@ def plot_video_with_surface(
                         )
                         if not inst_position.shape[1] == inst_radius.shape[0]:
                             inst_position = 0.5 * (
-                                inst_position[..., 1:] + inst_position[..., :-1]
+                                    inst_position[..., 1:] + inst_position[..., :-1]
                             )
 
                         rod_lines[rod_idx].set_xdata(inst_position[0])
@@ -297,7 +297,7 @@ def plot_video_with_surface(
                 ax.add_artist(sphere_artists[sphere_idx])
 
         ax.set_aspect("equal")
-        video_name_2D = folder_name + "2D_zy_" + video_name
+        video_name_2D = folder_name+"2D_zy_" + video_name
 
         with writer.saving(fig, video_name_2D, dpi):
             with plt.style.context("seaborn-whitegrid"):
@@ -309,7 +309,7 @@ def plot_video_with_surface(
                         )
                         if not inst_position.shape[1] == inst_radius.shape[0]:
                             inst_position = 0.5 * (
-                                inst_position[..., 1:] + inst_position[..., :-1]
+                                    inst_position[..., 1:] + inst_position[..., :-1]
                             )
 
                         rod_lines[rod_idx].set_xdata(inst_position[2])
@@ -390,7 +390,7 @@ def plot_video_with_surface(
                 ax.add_artist(sphere_artists[sphere_idx])
 
         ax.set_aspect("equal")
-        video_name_2D = folder_name + "2D_xz_" + video_name
+        video_name_2D = folder_name+"2D_xz_" + video_name
 
         with writer.saving(fig, video_name_2D, dpi):
             with plt.style.context("seaborn-whitegrid"):
@@ -402,7 +402,7 @@ def plot_video_with_surface(
                         )
                         if not inst_position.shape[1] == inst_radius.shape[0]:
                             inst_position = 0.5 * (
-                                inst_position[..., 1:] + inst_position[..., :-1]
+                                    inst_position[..., 1:] + inst_position[..., :-1]
                             )
 
                         rod_lines[rod_idx].set_xdata(inst_position[0])
@@ -437,11 +437,10 @@ def plot_video_with_surface(
         # See https://github.com/matplotlib/matplotlib/issues/8560/
         plt.close(plt.gcf())
 
-
 def plot_snake_velocity(
-    plot_params: dict,
-    period,
-    filename="slithering_snake_velocity.png",
+        plot_params: dict,
+        period,
+        filename="slithering_snake_velocity.png",
 ):
     time_per_period = np.array(plot_params["time"]) / period
     avg_velocity = np.array(plot_params["avg_velocity"])
@@ -470,7 +469,6 @@ def plot_snake_velocity(
     fig.legend(prop={"size": 20})
     fig.savefig(filename)
 
-
 def compute_projected_velocity(plot_params: dict, period):
 
     time_per_period = np.array(plot_params["time"]) / period
@@ -496,7 +494,7 @@ def compute_projected_velocity(plot_params: dict, period):
             center_of_mass[(i + 1) * period_step : (i + 2) * period_step]
             - center_of_mass[(i + 0) * period_step : (i + 1) * period_step],
             axis=0,
-        )
+            )
 
     # Average the rod directions over multiple periods and get the direction of the rod.
     direction_of_rod = np.mean(center_of_mass_averaged_over_one_period, axis=0)
@@ -528,3 +526,4 @@ def compute_projected_velocity(plot_params: dict, period):
         average_velocity_over_simulation[0],
         average_velocity_over_simulation[1],
     )
+
diff --git a/examples/RestartExample/restart_example.py b/examples/RestartExample/restart_example.py
index e97cbb3..816b028 100644
--- a/examples/RestartExample/restart_example.py
+++ b/examples/RestartExample/restart_example.py
@@ -3,11 +3,9 @@ This script is an example to how to use Pyelastica restart functionality.
 """
 
 import sys
-
 sys.path.append("../../")
 from elastica import *
 
-
 class RestartExampleSimulator(BaseSystemCollection, Constraints, Forcing, CallBacks):
     pass
 
@@ -64,7 +62,7 @@ restart_example_simulator.finalize()
 
 # After finalize you can load restart file. This step is important for current implementation of restart functions,
 # it is required to load restart files after the finalize step.
-LOAD_FROM_RESTART = False
+LOAD_FROM_RESTART = False 
 SAVE_DATA_RESTART = True
 restart_file_location = "data/"
 
@@ -78,16 +76,11 @@ dl = base_length / n_elem
 dt = 0.01 * dl
 total_steps = int(final_time / dt)
 
-time = integrate(
-    timestepper,
-    restart_example_simulator,
-    final_time,
-    total_steps,
-    restart_time=restart_time,
-)
+time = integrate(timestepper, restart_example_simulator, final_time, total_steps, restart_time=restart_time)
 
 # Save all the systems appended on the simulator class. Since in this example have only one system, under the
 # `restart_file_location` directory there is one file called system_0.npz . For each system appended on the simulator
 # separate system_#.npz file will be created.
 if SAVE_DATA_RESTART:
     save_state(restart_example_simulator, restart_file_location, time, True)
+
diff --git a/examples/RodContactCase/RodRodContact/rod_rod_contact_inclined_validation.py b/examples/RodContactCase/RodRodContact/rod_rod_contact_inclined_validation.py
index 8f71add..139bbe8 100644
--- a/examples/RodContactCase/RodRodContact/rod_rod_contact_inclined_validation.py
+++ b/examples/RodContactCase/RodRodContact/rod_rod_contact_inclined_validation.py
@@ -35,9 +35,7 @@ poisson_ratio = 0.5
 shear_modulus = E / (poisson_ratio + 1.0)
 
 # Rod orientations
-start = np.zeros(
-    3,
-)
+start = np.zeros(3,)
 inclination = np.deg2rad(30)
 direction = np.array([0.0, np.cos(inclination), np.sin(inclination)])
 normal = np.array([0.0, -np.sin(inclination), np.cos(inclination)])
@@ -58,7 +56,7 @@ rod_one = CosseratRod.straight_rod(
     nu,
     E,
     poisson_ratio,
-    shear_modulus=shear_modulus,
+    shear_modulus = shear_modulus,
 )
 
 rod_one.velocity_collection[:] += 0.05 * -normal.reshape(3, 1)
@@ -96,7 +94,9 @@ inclined_rod_rod_contact_sim.connect(rod_one, rod_two).using(
 
 # Add call backs
 class RodCallBack(CallBackBaseClass):
-    """ """
+    """
+
+    """
 
     def __init__(self, step_skip: int, callback_params: dict):
         CallBackBaseClass.__init__(self)
@@ -130,9 +130,7 @@ post_processing_dict_rod1 = defaultdict(
 )  # list which collected data will be append
 # set the diagnostics for rod and collect data
 inclined_rod_rod_contact_sim.collect_diagnostics(rod_one).using(
-    RodCallBack,
-    step_skip=step_skip,
-    callback_params=post_processing_dict_rod1,
+    RodCallBack, step_skip=step_skip, callback_params=post_processing_dict_rod1,
 )
 
 post_processing_dict_rod2 = defaultdict(
@@ -140,9 +138,7 @@ post_processing_dict_rod2 = defaultdict(
 )  # list which collected data will be append
 # set the diagnostics for rod and collect data
 inclined_rod_rod_contact_sim.collect_diagnostics(rod_two).using(
-    RodCallBack,
-    step_skip=step_skip,
-    callback_params=post_processing_dict_rod2,
+    RodCallBack, step_skip=step_skip, callback_params=post_processing_dict_rod2,
 )
 
 inclined_rod_rod_contact_sim.finalize()
diff --git a/examples/RodContactCase/RodRodContact/rod_rod_contact_parallel_validation.py b/examples/RodContactCase/RodRodContact/rod_rod_contact_parallel_validation.py
index 09b18ca..7ba4b44 100644
--- a/examples/RodContactCase/RodRodContact/rod_rod_contact_parallel_validation.py
+++ b/examples/RodContactCase/RodRodContact/rod_rod_contact_parallel_validation.py
@@ -35,9 +35,7 @@ poisson_ratio = 0.5
 shear_modulus = E / (poisson_ratio + 1.0)
 
 # Rod orientations
-start = np.zeros(
-    3,
-)
+start = np.zeros(3,)
 inclination = np.deg2rad(0)
 direction = np.array([0.0, np.cos(inclination), np.sin(inclination)])
 normal = np.array([0.0, -np.sin(inclination), np.cos(inclination)])
@@ -57,7 +55,7 @@ rod_one = CosseratRod.straight_rod(
     density,
     nu,
     E,
-    shear_modulus=shear_modulus,
+    shear_modulus = shear_modulus,
 )
 
 rod_one.velocity_collection[:] += 0.05 * -normal.reshape(3, 1)
@@ -91,7 +89,9 @@ parallel_rod_rod_contact_sim.connect(rod_one, rod_two).using(
 
 # Add call backs
 class RodCallBack(CallBackBaseClass):
-    """ """
+    """
+
+    """
 
     def __init__(self, step_skip: int, callback_params: dict):
         CallBackBaseClass.__init__(self)
@@ -125,9 +125,7 @@ post_processing_dict_rod1 = defaultdict(
 )  # list which collected data will be append
 # set the diagnostics for rod and collect data
 parallel_rod_rod_contact_sim.collect_diagnostics(rod_one).using(
-    RodCallBack,
-    step_skip=step_skip,
-    callback_params=post_processing_dict_rod1,
+    RodCallBack, step_skip=step_skip, callback_params=post_processing_dict_rod1,
 )
 
 post_processing_dict_rod2 = defaultdict(
@@ -135,9 +133,7 @@ post_processing_dict_rod2 = defaultdict(
 )  # list which collected data will be append
 # set the diagnostics for rod and collect data
 parallel_rod_rod_contact_sim.collect_diagnostics(rod_two).using(
-    RodCallBack,
-    step_skip=step_skip,
-    callback_params=post_processing_dict_rod2,
+    RodCallBack, step_skip=step_skip, callback_params=post_processing_dict_rod2,
 )
 
 parallel_rod_rod_contact_sim.finalize()


Test patch:
diff --git a/tests/test_rod_initialisation.py b/tests/test_rod_initialisation.py
index 097b2c0..0ae9d2f 100644
--- a/tests/test_rod_initialisation.py
+++ b/tests/test_rod_initialisation.py
@@ -152,7 +152,7 @@ class MockRodForTest:
             density,
             nu,
             youngs_modulus,
-            alpha_c=0.964,
+            alpha_c=(27.0 / 28.0),
             *args,
             **kwargs
         )
@@ -741,8 +741,8 @@ def test_shear_matrix_for_varying_shear_modulus(n_elems, shear_modulus):
     np.fill_diagonal(
         correct_shear_matrix[:],
         [
-            0.964 * shear_modulus * base_area,
-            0.964 * shear_modulus * base_area,
+            (27.0 / 28.0) * shear_modulus * base_area,
+            (27.0 / 28.0) * shear_modulus * base_area,
             youngs_modulus * base_area,
         ],
     )
@@ -850,10 +850,10 @@ def test_inertia_shear_bend_matrices_for_varying_radius():
 
     correct_shear_matrix = np.array(
         [
-            [23296.11783, 23296.11783, 31415.92654],
-            [93184.47129, 93184.47129, 125663.7061],
-            [209665.06048, 209665.06048, 282743.33882],
-            [372737.88191, 372737.88191, 502654.82],
+            [23303.02243, 23303.02243, 31415.92654],
+            [93212.08972, 93212.08972, 125663.7061],
+            [209727.20187, 209727.20187, 282743.33882],
+            [372848.35889, 372848.35889, 502654.82],
         ]
     )
 
@@ -1547,7 +1547,9 @@ def test_straight_rod(n_elems):
     inv_mass_second_moment_of_inertia = np.linalg.inv(mass_second_moment_of_inertia)
     # Shear/Stretch matrix
     shear_matrix = np.zeros((3, 3), np.float64)
-    np.fill_diagonal(shear_matrix, [0.964 * G * A0, 0.964 * G * A0, E * A0])
+    np.fill_diagonal(
+        shear_matrix, [(27.0 / 28.0) * G * A0, (27.0 / 28.0) * G * A0, E * A0]
+    )
     # Bend/Twist matrix
     bend_matrix = np.zeros((3, 3), np.float64)
     np.fill_diagonal(bend_matrix, [E * I0_1, E * I0_2, G * I0_3])