Dear PyElastica Team,

I am currently trying to add an Interaction plane (or any arbitrary obstacle) to the slithering snake example, so that the snake "bounces of" the plane and continues to move into a different direction. However, I wasn't able to find out the correct way to do it. I tried adding the following code into the example

side_plane = np.array([0.0, 0.0, 3.0]) side_normal = np.array([0.0, 0.0, 1.0]) surface_tol=1e-4

plane_sim = InteractionPlane(k=1.0,nu=1e-6,plane_origin=side_plane,plane_normal=normal_plane)

I do not obtain any error messages. Unfortunately, the code does not seem to have any effect on the snake. I am thankful for any tips. I am very new to the PyElastica (and Github) implementation, however I think this is a great project!

Features and Changes

This issue fixes #122 by implementing the following changes and features:

1. Enables the connection of Cosserat rods and rigid bodies in any arbitrary combination for the following joint types: FreeJoint, HingeJoint, and FixedJoint
2. The FreeJoint constrains the translational movement of both systems at the specified connection point. As the point of joint connection doesn't necessarily coincide with the center of mass of the rigid body or the node of the rod, the user can use the point_system_one and point_system_two parameters to specify how the connection point relates to the node / CoM in the local frame of the system.
3. For the HingeJoint, the link_direction of the second system is now extracted from the director_collection instead of subtracting node positions
4. The FixedJoint class now also accepts the point_system_one and point_system_two parameters and forwards them to the parent class FreeJoint

Examples of Usage

How to connect two Cosserat rods together using a spherical joint with a gap of 0.01 m in between.

simulator.connect(rod_one, rod_two).using(
    FreeJoint,
    k=1e4,
    nu=1,
    point_system_one=np.array([0.0, 0.0, 0.005]),
    point_system_two=np.array([0.0, 0.0, -0.005]),
)

How to connect the distal end of a CosseratRod with the base of a cylinder using a spherical joint.

simulator.connect(rod, cylinder).using(
    FreeJoint,
    k=1e4,
    nu=1,
    point_system_two=np.array([0.0, 0.0, -cylinder.length / 2.]),
)

How to connect a cosserat rod with the base of a cylinder using a fixed joint, where the cylinder is rotated by 45 degrees around the y-axis.

from scipy.spatial.transform import Rotation
simulator.connect(rod, cylinder).using(
      FixedJoint,
      k=1e5,
      nu=1e0,
      kt=1e3,
      nut=1e-3,
      point_system_two=np.array([0, 0, -cylinder.length / 2]),
      rest_rotation_matrix=Rotation.from_euler('y', np.pi / 4, degrees=False).as_matrix(),
)

Results from JointCases

Spherical joint:

System consisting of two rods and one cylinder connected with spherical joints. The sinusoidal force is applied to the second rod. Example can be run using:

python examples/JointCases/spherical_joint.py

https://user-images.githubusercontent.com/2360366/180027045-d9753a1a-56ec-4459-a719-a13c7398620e.mp4

Fixed joint:

System consisting of two rods and one cylinder connected with fixed joints. The sinusoidal force is applied to the second rod. Example can be run using:

python examples/JointCases/fixed_joint.py

https://user-images.githubusercontent.com/2360366/180449881-c806dd79-324b-40d9-9f57-b6257c7f8978.mp4

TODOs

• [x] Write tests for connecting rods with cylinders, particularly for non-zero point_system_one and point_system_two parameters
• [x] Add cylinder to spherical_joint.py example as third element.
• [x] Add cylinder to fixed_joint.py example as third element

Motivation

Currently, there exists only the option to either have a) no boundary conditions or b) fully fixed boundary conditions. The only functionality in between is that for example translations are allowed, but rotations are constrained and vice-versa. However in practice, there are many situations where we would like to allow certain translations (for example along a plane) or rotations (for example allow yawing of the end of the rod), but constrain all other DoF. The user of the library would have to implement this functionality all themselves inheriting from the FreeBC class. In many cases, such an implementation is not straight-forward. Thus, the library should provide the user the option to configure the BC themselves in a simple fashion. Although it obviously will never cover all the edge-cases a specific project might require, it will give the user an idea / hint / guidance on how to implement a similar solution.

Provided functionality

This PR add the following functionality through the ConfigurableFixedConstraint class:

• Using boolean arrays translational_constraint_selector and rotational_constraint_selector, the user can specify which Degrees of Freedom to allow / constrain at the specified nodes.
• These boolean arrays have to be specified in the inertial frame.
• The ConfigurableFixedConstraint will now only re-set the positions, velocities, directors and omegas to the saved, fixed values for the constrained dimensions. The other DoF are allowed to freely move / rotate.
• We can simplify the codebase by now having FixedConstraint be a child class of ConfigurableFixedConstraint. FixedConstraint will still constrain all dimension, e.g. translational_constraint_selector=np.array([True, True, True]) and rotational_constraint_selector=np.array([True, True, True])

Examples for usage

Below examples motivate the usage of this new class.

# How to fix all translational and rotational DoF except allowing twisting around z-axis in inertial frame:

simulator.constrain(rod).using(
    ConfigurableFixedConstraint,
    constrained_position_idx=(0,),
    constrained_director_idx=(0,),
    translational_constraint_selector=np.array([True, True, True]),
    rotational_constraint_selector=np.array([True, True, False]),
)

# How to allow the end of the rod to move in the x-y plane and allow all rotational DoF:

simulator.constrain(rod).using(
    ConfigurableFixedConstraint,
    constrained_position_idx=(-1,),
    translational_constraint_selector=np.array([True, True, False]),
)

I also added an example in examples/BoundaryConditionCases/configurable_fixed_constraint.py, which applies torsion to end tip of a rod. With a fully FixedConstraint the resulting orientation of the first node of the rod would look like this:

Fully Fixed BC:

When we now use rotational_constraint_selector=np.array([True, True, False]) (e.g. allow for yawing / rotation around the z-axis), the orientation of the first node of the rod looks like this:

Free Yawing:

TODO's

• [x] Tests for ConfigurableFixedConstraint class
• [x] Numba-compatible implementation of nb_constraint_rotational_values
• [x] Rotate omega_collection into the inertial frame before applying zero-speed constraints

Request for ideas / feedback

I would like to ask the maintainers specifically for feedback on the current prototype of nb_constraint_rotational_values. Obviously, this is currently not numba-compatible. But more importantly, I would like to ask for your thoughts if constraining the Euler angles is a suitable way to go and if you have other ideas to allow the rotation around certain axis and constrain the rotation around other axis?

Hello PyElastica,

I am a student and I am trying to apply the library 'PyElastica' in my project for a rod with a pre-curvature. But I could only create a straight rod. Could you please advise me if there is a way to make a pre-curved rod?

Thank you in advance for your help and congratulations for this beautiful work !

Émerson

Dear PyElastica,

thank you in advance. While using powerful PyElastica to model a soft robot, which is consisted of 3 parallel rods connected from each other. Luckily, I found a similar one from your paper: Topology, Geometry, and Mechanics of Strongly Stretched and Twisted Filaments: Solenoids, Plectonemes, and Artificial Muscle Fibers, but i can not find any reference code. So could you share the modeling code like boundary condition or how the bilateral rods connected to the centeral one. Thank you for your time and help.

BR Daphne

Hello! Is there any way to apply torque( at the end of the rod) in the same direction as rod is fixed. For example I have a rod in y direction and I would like to apply a torque in this direction. I have created a calss called "OnedirectionTorque"

class OnedirectionTorque(NoForces): def init(self,torque_end,direction): super(EndPointTorques, self).init() self.torque_end= (torque_end * direction).reshape(3, 1) # defined in global frame, shape(3,1) def apply_torques(self, system, time: np.float = 0.0): system.external_torques[..., -1] += system.director_collection[..., -1]@self.torque_end

I also have add it this piece of lines in examples/TimoshenkoBeamCase/timoshenko.py direction = np.array([0.0, 1.0, 0.0]) toque_end = -10 timoshenko_sim.add_forcing_to(shearable_rod).using( UniformTorques,toque_end,direction)

Is it implemented correctly

Thanks in advance

Hi PyElasitica teams,

I want to simulate the process of cable operation. I set one end of the cable fixed, applied an upward force in the middle of the cable, and applied gravity to the whole cable.

The unfixed end of the real cable will sag naturally under the influence of gravity, as shown in above figure, but in my simulation results, the unfixed end is still lifting until the whole cable is straightened, as shown in below figure.

Here is my code:

Hi, I'm developing a research project and found this library. However, I'm doing mainly cpp and need a lib to be able to integrate it into the current system.

After checking the code, it seems the main lib is written in python? Or do you have another version that is written in C/C++ elsewhere?

Thank you

Fixes #112 by providing an alternative add-on damping module (mixin class) for internal dissipation of Cosserat rods, with the following advantages:

1. Damping can be added as a Mixin class when it is needed for a rod, thus making damping modular.
2. The analytical version of the original internal damper ExponentialDamper is unconditionally stable, and removes any timestep restrictions coming in from damping.
3. With the incoming modularity, other damping models can be conveniently implemented.

Additionally, a DeprecationWarning is added when the damping constant is provided in the rod constructor (soon the option will be removed in the following release-0.4).

The examples will be refactored (with help from @armantekinalp), along with confirmation, that the dynamics of the individual cases stay the same:

• [x] axial stretching
• [x] binder stuff?
• [x] butterfly
• [x] continuum flagella
• [x] continuum snake
• [x] Experimental cases - parallel connection
• [x] flexible swinging pendulum
• [x] friction validation cases
• [x] helical buckling
• [x] joint cases
• [x] muscular flagella
• [x] muscular snake
• [x] restart example
• [x] rod-rod contact cases
• [x] rod self contact cases
• [x] timoshenko

@armantekinalp after finishing the confirmations, we should also check if we can raise timesteps in these cases.

Fixes #80, by

1. Adding momentum conserving node-to-element interpolation for velocity.
2. Renaming older node-to-element interpolation as to be used for position.
3. Specifically naming the node-to-element interpolation used for mass/forces (edge special treatment).

@armantekinalp I ran the snake and flagella, and they produce results as expected.

Fixes #131 for the FixedJoint class:

• Now any arbitrary torques applied to rods or rigid bodies connected by a FixedJoint are compensated (e.g. more general implementation)
• Rotation frame deviation is computed using rotation vectors
• Rotational spring-damper system is used to restrict the rotational mode
• Spring: spring constant kt scales the rotational frame deviation between systems one and two
• Damper: damping constant nut scales the difference between angular velocities of system one and two
• Introduce functionality of static rotation offset between systems one and two: records the initial relative rotation between the two systems and enforces this static rotation throughout the entire simulation. If use_static_rotation==False, the static rotation is set to an identity matrix and restoring torques are applied to fully align the local frames of both systems at the join.
• Example examples/JointCases/fixed_joint_torsion.py illustrates how this implementation generates torsional / twisting torques to prevent torsional deviations
• Use common diagnostic callback for all joint cases examples The attached plots demonstrates the results of the examples/JointCases/fixed_joint_torsion.py example:

Hello,

I am looking at the Timoschenko Beam example and I noticed that when one applies an end force to the beam the x-value of the solution doesn't change compared to the starting point. That would imply that the Beam gets longer over time. Is this intended or is there a workaround to this problem?

Thanks and kind regards, Matei

Hi,

I'm currently looking into different tools to do static and dynamic analysis of cosserat rods. Hence my question: Is it possible to do static analysis with PyElastica without time integration? And if not, is it possible without too much effort to include this feature with additional code by myself without manipulating PyElastica itself?

Thank you for your feedback. Bests

Hi,

I am trying to implement the helical motion benchmark from this paper:

https://www.researchgate.net/publication/329252597_A_Geometrically_Exact_Model_for_Soft_Continuum_Robots_The_Finite_Element_Deformation_Space_Formulation

To do this I need to define a follower force and follower torque at the tip of the rod. The force needs to remain perpendicular to the rod under deformation. This is my current attempt at defining a follower force.

class endpointTorque(NoForces):
    def __init__(
        self,
        torque=np.array([0.0, 0.0, 0.0]),
        rampUpTime = 0.0
    ):  
        self.torque=torque
        self.rampUpTime=rampUpTime

    def apply_torques(self, system, time: np.float64 = 0.0):
        if self.rampUpTime != 0:
            factor = min(1, time / self.rampUpTime)
            system.external_torques[..., -1] += system.director_collection[...,-1] @ (self.torque*factor)
        else:
            system.external_torques[..., -1] += system.director_collection[...,-1] @ self.torque



class endpointForce(NoForces):
    def __init__(
        self,
        force=np.array([0.0, 0.0, 0.0]),
        rampUpTime = 0.0
    ):  
        self.force=force
        self.rampUpTime = rampUpTime

    def apply_forces(self, system, time: np.float64 = 0.0):
        if self.rampUpTime != 0:
            factor = min(1, time / self.rampUpTime)
            system.external_forces[..., -1] += system.director_collection[...,-1] @ (self.force*factor)
        else:
            system.external_forces[..., -1] += system.director_collection[...,-1] @ self.force

Unfortunately, this doesn't result in the desired helical twisting of the rod. Instead, the rod is stretched in one direction.

Instead, I would expect something like this:

Do you know what I am doing wrong? Thank you very much.

Currently there is an issue with Windows CI builds, because of that we are removing Windows build in #206 . However, we should include it in future. @bhosale2 could you please attach the relevant links.

Currently, adding damping to the rod via AnalyticalLinearDamper results in the following warning message:

Per the warning to make a soft change, we should convert the warning message to an error message if the user provides internal damping other than None when initializing a Cosserat rod, and remove related functionality from the cosserat rod equations.

Thoughts @armantekinalp and @skim0119?

Currently, in example cases, the modules from elastica are imported directly as from elastica import *. This can become an issue when using elastica across packages (where elastica is a plugin), where modules from other packages can collide in namings with the modules from elastica. To avoid this, it may be helpful to use functionalities in elastica as:

import elastica as ea; ea.Damping(....) to make it explicit where the functionality comes from.

@skim0119 and @armantekinalp thoughts?