Twisting torque is not constrained by FixedJoint and HingeJoint

[mstoelzle]

Describe the bug

The current implementation of the FixedJoint based a position spring does not constrain any twisting / torsional torques applied to a straight rod. This happens, because the current approach (essentially) computes an error between the tangential vector of rod 1 and tangential error of rod 2. As rotating around the z-axis / tangential vector does not change this vector, the twisting / torsion is not constrained by the FixedJoint.

In the example below, we have the classic example of two rods connected by a FixedJoint. They are both initially in straight configuration and we apply a UniformTorque to the second rod. For a twisting torque, we notice the faulty behaviour, while for a bending torque everything works as expected.

To Reproduce

__doc__ = """Fixed joint example, for detailed explanation refer to Zhang et. al. Nature Comm.  methods section."""

import matplotlib.pyplot as plt
import numpy as np
import sys

# FIXME without appending sys.path make it more generic
sys.path.append("../../")
from elastica import *
from examples.JointCases.external_force_class_for_joint_test import (
    EndpointForcesSinusoidal,
)
from examples.JointCases.joint_cases_callback import JointCasesCallback
from examples.JointCases.joint_cases_postprocessing import (
    plot_position,
    plot_video,
    plot_video_xy,
    plot_video_xz,
)

class FixedJointSimulator(
    BaseSystemCollection, Constraints, Connections, Forcing, CallBacks
):
    pass

fixed_joint_sim = FixedJointSimulator()

# setting up test params
n_elem = 10
direction = np.array([0.0, 0.0, 1.0])
normal = np.array([0.0, 1.0, 0.0])
roll_direction = np.cross(direction, normal)
base_length = 0.2
base_radius = 0.007
base_area = np.pi * base_radius ** 2
density = 1750
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
start_cylinder = start_rod_2 + direction * base_length

# Create rod 1
rod1 = CosseratRod.straight_rod(
    n_elem,
    start_rod_1,
    direction,
    normal,
    base_length,
    base_radius,
    density,
    nu,
    E,
    shear_modulus=shear_modulus,
)
fixed_joint_sim.append(rod1)
# Create rod 2
rod2 = CosseratRod.straight_rod(
    n_elem,
    start_rod_2,
    direction,
    normal,
    base_length,
    base_radius,
    density,
    nu,
    E,
    shear_modulus=shear_modulus,
)
fixed_joint_sim.append(rod2)

# Apply boundary conditions to rod1.
fixed_joint_sim.constrain(rod1).using(
    OneEndFixedBC, constrained_position_idx=(0,), constrained_director_idx=(0,)
)

# Connect rod 1 and rod 2
fixed_joint_sim.connect(
    first_rod=rod1, second_rod=rod2, first_connect_idx=-1, second_connect_idx=0
).using(FixedJoint, k=1e5, nu=0, kt=5e3)

# Add forces to rod2
fixed_joint_sim.add_forcing_to(rod2).using(
    UniformTorques,
    torque=5e-3,
    direction=np.array([0., 0., 1.])
)

pp_list_rod1 = defaultdict(list)
pp_list_rod2 = defaultdict(list)

fixed_joint_sim.collect_diagnostics(rod1).using(
    JointCasesCallback, step_skip=1000, callback_params=pp_list_rod1
)
fixed_joint_sim.collect_diagnostics(rod2).using(
    JointCasesCallback, step_skip=1000, callback_params=pp_list_rod2
)

fixed_joint_sim.finalize()
timestepper = PositionVerlet()
# timestepper = PEFRL()

final_time = 10
dl = base_length / n_elem
dt = 1e-5
total_steps = int(final_time / dt)
print("Total steps", total_steps)
integrate(timestepper, fixed_joint_sim, final_time, total_steps)

def plot_orientation_vs_time(title, time, directors):
    from scipy.spatial.transform import Rotation as R
    quat = []
    for t in range(len(time)):
        quat_t = R.from_matrix(directors[t].T).as_quat()
        quat.append(quat_t)
    quat = np.array(quat)

    plt.figure(num=title)
    plt.plot(time, quat[:, 0], label="x")
    plt.plot(time, quat[:, 1], label="y")
    plt.plot(time, quat[:, 2], label="z")
    plt.plot(time, quat[:, 3], label="w")
    plt.title(title)
    plt.legend()
    plt.xlabel("Time [s]")
    plt.ylabel("Quaternion")
    plt.show()

plot_orientation_vs_time("Orientation of last node of rod 1", pp_list_rod1["time"],
                         np.array(pp_list_rod1["director"])[..., -1])
plot_orientation_vs_time("Orientation of last node of rod 2", pp_list_rod2["time"],
                         np.array(pp_list_rod2["director"])[..., -1])

Environment

• Python=3.10.5, pip=22.1.2
• Numpy/numba version: 22.1.2
• PyElastica version: 0.2.3
• OS, Device: macOS / macBook Pro 15'' 2019

Expected behavior

We expect the FixedJoint to constrain all 6 DoF between the two rods: e.g. no translations or rotations. This should also include twisting / torsional movements in straight rod configuration.

Screenshots

Orientation of last node of rod 1 for torsion torque applied to rod 2:

Here, we would expect the rod 1 to rotate a bit around its z-axis to counteract the torsion torque applied at it's tip at the FixedJoint. However, it doesn't move at all.

Orientation of last node of rod 2 for torsion torque applied to rod 2:

Here, we would expect the rod 2 to rotate a bit around its z-axis to counteract the torsion torque. However, it is clearly not constrained at its base and just freely rotates.

Orientation of last node of rod 1 for bending torque applied to rod 2:

Here, we would expect the rod 1 to rotate a bit around its y-axis to counteract the torsion torque applied at it's tip at the FixedJoint. This is exactly what happens and after a few seconds the elastic systems reaches its steady state.

Orientation of last node of rod 12for bending torque applied to rod 2:

Here, we would expect the rod 2 to rotate a bit around its y-axis to counteract the torsion torque applied. This is exactly what happens and after a few seconds the elastic systems reaches its steady state.

[mstoelzle]

@xzhan139 I suggest moving to torsional springs for both the hinge joints and fixed joints. Could you please give me some hints about the previous implementation you have tried out?

[mstoelzle]

Do you guys want me to land the PR to address this issue together with #122 or in separate PR (would need to do this first then) ?

[xzhan139]

Hi, I think these two questions of yours are interconnected - yes, we opted for moment arm type of connection, not the torsional spring connection, and that's exactly why the twisting torques are not constrained in a fixed joint. This was not a critical issue in our example cases, but it would definitely be a problem if you would need the fixed joint to resist twisting torques, like the example you have here. So essentially, all you need is to switch to torsional spring connection. I am not sure if we have the torsional spring connection implemented in Python (I work primarily with C++), if not, this maybe something you can contribute (you can coordinate with Arman on this). But anyways, I can briefly describe how you can do this, similar with the current scheme - Let's start with constraining one direction of a joint, say, we want to constrain the angle between d_1 of both rod to be 90 degrees. What you need to do is to compute the current angle between the two directors, say 'angle_current', and compute the error 'angle_current - 90'. Then apply torque (scaled with error) with opposite directions to the nearest element at the connection on both rods. Be careful: - the coefficient (scaling factor on error) can be sensitive, especially when you have a complex interconnected structure; and make sure torques have the correct directions, so that the effect of both torques is to reduce the error, otherwise it will explode. Once you have one direction fixed, fixing another one you will have a hinge joint, and fixing all three directions you will have a complete fixed joint, resisting to torques in all directions.

[armantekinalp]

Hi @mstoelzle

You are right current implementation of fixed joint only constrains the bending and it allows twist. Currently PyElastica does not have a joint class uses torsional springs, so I think it will be a good addition.

Do you guys want me to land the PR to address this issue together with #122 or in separate PR (would need to do this first then) ?

I think let's make two different PRs so that each PR has minimal review.

[skim0119]

I'm adding this issue to docstring for FixedJoints. (aea196b)