formatting

elastica/joint.py

@@ -1056,8 +1056,9 @@ class SelfContact(FreeJoint):
 
     """
 
-    def __init__(self, k, nu):
-        super().__init__(k, nu)
+    def __init__(self, k, nu, **kwargs):
+        # super().__init__(k, nu)
+        super().__init__(np.array(k), np.array(nu))
 
     def apply_forces(self, rod_one: RodType, index_one, rod_two: SystemType, index_two):
         # del index_one, index_two

elastica/rod/knot_theory.py

@@ -251,16 +251,29 @@ def _compute_twist(center_line, normal_collection):
         alpha /= 2 * np.pi
         gamma /= 2 * np.pi
         twist_temp = alpha + gamma
+
+        # Check if there is any Nan. Nan's appear when rod tangents are parallel to each other.
+        idx = np.where(np.isnan(twist_temp))[0]
+        for i in idx:
+            # compute angle between two vectors
+            dot_product = np.dot(
+                projection_of_normal_collection[:, i],
+                projection_of_normal_collection[:, i + 1],
+            )
+            angle = np.arccos(dot_product)
+            cross = np.cross(
+                projection_of_normal_collection[:, i],
+                projection_of_normal_collection[:, i + 1],
+            )
+            angle *= np.sign(np.dot(tangent[:, i], cross))
+            twist_temp[i] = angle / (2 * np.pi)
+
         # Make sure twist is between (-1/2 to 1/2) as defined in pg 313 Klenin & Langowski 2000
         idx = np.where(twist_temp > 0.5)[0]
         twist_temp[idx] -= 1
         idx = np.where(twist_temp < -0.5)[0]
         twist_temp[idx] += 1
 
-        # Check if there is any Nan. Nan's appear when rod tangents are parallel to each other.
-        idx = np.where(np.isnan(twist_temp))[0]
-        twist_temp[idx] = 0.0
-
         local_twist[k, :] = twist_temp
         total_twist[k] = np.sum(twist_temp)
 

examples/ArtificialMusclesCases/MuscleTuningCases/muscle_actuation_amount_tuner.py

@@ -3,53 +3,148 @@
 from elastica import *
 from examples.ArtificialMusclesCases import *
 from examples.ArtificialMusclesCases.SingleMuscleCases.single_muscle_simulator import (
-    muscle_contraction_simulation,
+    muscle_stretching_actuation_simulator,
 )
 
 
-tuning_sim_settings_dict = {
+# how the muscle is constrained
+constrain_settings_dict = {
+    "isometric_test": False,
+    "isobaric_test": True,
+    "constrain_start_time": 0.0,
+}
+
+# what external forcing it is experiencing
+forcing_settings_dict = {
+    "force_mag": 25.0,
+    "desired_muscle_strain": 1.0,
+    "start_force_direction": np.array([0.0, 0.0, -1.0]),
+    "end_force_direction": np.array([0.0, 0.0, 1.0]),
+    "ramp_up_time": 1.0,
+    "forcing_start_time": 0.0,
+}
+
+# how is it actuated
+actuation_settings_dict = {
+    "actuation": True,
+    "actuation_duration": 6.0,
+    "start_time": 0.0,
+}
+
+# self contact settings
+self_contact_settings_dict = {
+    "self_contact": True,
+    "contact_radius_ratio": 1 / 10,
+    "k_val": 0e0,
+    "nu": 0.0,
+    "k_repulsive_val": 4e0,  # 4e0
+    "separation_distance": 1e-3,
+}
+
+# fiber-fiber interaction settings
+# fiber_connection_settings_dict = {
+#     "connection_range": 1,
+#     "k_val": 0.0,
+#     "nu": 0.0,
+#     "k_repulsive_val": 0.0,
+#     "friction_coefficient": 0.0,
+#     "velocity_damping_coefficient": 1e5,
+# }
+
+fiber_connection_settings_dict = {
+    "connection_range": 2,
+    "k_val": 2.5e2 / 16,
+    "nu": 0.0,
+    "k_repulsive_val": 1e2,
+    "friction_coefficient": 5e-1,
+    "velocity_damping_coefficient": 1e5,
+}
+
+# convert settings dicts to classes for ease
+constrain_settings = Dict2Class(constrain_settings_dict)
+forcing_settings = Dict2Class(forcing_settings_dict)
+actuation_settings = Dict2Class(actuation_settings_dict)
+self_contact_settings = Dict2Class(self_contact_settings_dict)
+fiber_connection_settings = Dict2Class(fiber_connection_settings_dict)
+
+
+sim_settings_dict = {
     "sim_name": "Tuning",
     "final_time": 10,  # seconds
-    "untwisting_start_time": 0,  # seconds
-    "time_untwisting": 1,  # seconds
     "rendering_fps": 20,
-    "contraction": True,
+    "LaplaceFilter": False,
+    "LaplaceFilterOrder": 7,
     "plot_video": True,
     "save_data": False,
     "return_data": True,
     "povray_viz": False,
-    "isometric_test": True,
-    "isobaric_test": False,
-    "self_contact": False,
-    "muscle_strain": 0.0,
-    "force_mag": 0.0,
+    "save_folder": "",  # do not change
+    "additional_identifier": "",
 }
+sim_settings = Dict2Class(sim_settings_dict)
 
+strain_difference = 10
+tuning_muscle = Liuyang_monocoil(experimental_data=True)
+# tuning_muscle = Samuel_supercoil_stl(experimental_data=False)
 
-tuning_sim_settings = Dict2Class(tuning_sim_settings_dict)
 
-force_difference = 10
-tuning_muscle = Liuyang_monocoil()
-k = 0.01
+n_muscles = 1
+for i in range(len(tuning_muscle.geometry.n_ply_per_coil_level)):
+    n_muscles *= tuning_muscle.geometry.n_ply_per_coil_level[i]
+base_area = np.pi * (tuning_muscle.geometry.fiber_radius ** 2)
+forcing_settings.force_mag = tuning_muscle.active_force_experimental[0] / (
+    1e-3
+    * n_muscles
+    * tuning_muscle.properties.youngs_modulus
+    * tuning_muscle.sim_settings.E_scale
+    * base_area
+)
+# forcing_settings.force_mag = 0.0
 
+k = 1e-3
 
-while abs(force_difference) > 1e-6:
+current_path = os.getcwd()
+sim_settings.save_folder = os.path.join(
+    current_path, sim_settings.sim_name + "/" + tuning_muscle.name + "/data"
+)
+# tuning_muscle.sim_settings.actuation_kappa_change = 0.1
+while abs(strain_difference) > 1e-6:
     print(
         "Current Kappa Change:" + str(tuning_muscle.sim_settings.actuation_kappa_change)
     )
-    data = muscle_contraction_simulation(
-        input_muscle=tuning_muscle, sim_settings=tuning_sim_settings
+    data = muscle_stretching_actuation_simulator(
+        input_muscle=tuning_muscle,  # the type of muscle used
+        constrain_settings=constrain_settings,  # how the muscle is constrained
+        forcing_settings=forcing_settings,  # what external forcing it is experiencing
+        actuation_settings=actuation_settings,  # how is it actuated
+        self_contact_settings=self_contact_settings,  # self contact settings
+        fiber_connection_settings=fiber_connection_settings,  # fiber-fiber interaction settings
+        sim_settings=sim_settings,
     )
 
-    internal_force = np.zeros_like(np.array(data[0]["internal_force"]))
+    # internal_force = np.zeros_like(np.array(data[0]["internal_force"]))
+    # for muscle_rods in data:
+    #     internal_force += np.array(muscle_rods["internal_force"])
+
+    # force_difference = (
+    #     tuning_muscle.total_force_experimental[0]
+    #     - internal_force[-1, 2, 0] * tuning_muscle.sim_settings.E_scale
+    # )
+
+    centerline_position = np.zeros_like(np.array(data[0]["position"]))
     for muscle_rods in data:
-        internal_force += np.array(muscle_rods["internal_force"])
+        centerline_position += np.array(muscle_rods["position"]) / n_muscles
 
-    force_difference = (
-        tuning_muscle.total_force_experimental[0]
-        - internal_force[-1, 2, 0] * tuning_muscle.sim_settings.E_scale
+    n = centerline_position.shape[-1]
+    strain_difference = np.dot(
+        tuning_muscle.geometry.direction,
+        (
+            (centerline_position[-1, :, -1] - centerline_position[-1, :, 0])
+            - (centerline_position[0, :, -1] - centerline_position[0, :, 0])
+        )
+        / (centerline_position[0, :, -1] - centerline_position[0, :, 0]),
     )
-    print("Current Force Difference:" + str(force_difference))
-    tuning_muscle.sim_settings.actuation_kappa_change += k * force_difference
+    print("Current Strain Difference:" + str(strain_difference))
+    tuning_muscle.sim_settings.actuation_kappa_change += k * strain_difference
 
 print("Tuned Kappa Change: " + str(tuning_muscle.sim_settings.actuation_kappa_change))