GazzolaLab / PyElastica

Python implementation of Elastica, an open-source software for the simulation of assemblies of slender, one-dimensional structures using Cosserat Rod theory.
https://www.cosseratrods.org
MIT License
223 stars 109 forks source link

Replicating side-winding results from your Nat-comm snake paper #306

Closed AkashVardhan7 closed 10 months ago

AkashVardhan7 commented 12 months ago

Hello again., I have been trying to replicate some of the side-winding results from your paper, and added a separate class to compute the torques from the lifting wave in the external_forces.py file like the MuscleTorques class, and then in the main simulation file, I add the forcing from the lifting wave using the liftingTorques class.. like it's done for the lateral wave with the MuscleTorques class..however it seems I will have to re-build/compile elastica again as I get the following error message AttributeError: module 'elastica' has no attribute 'LiftingTorques'. I am pasting the codes for these files...kindly take a peek and let me know if I am on the right track before I re-compile elastica to reflect the modification with the LiftingTorques.

here's the modified external_forces.py file with the modified Lifting wave doc = """ Numba implementation module for boundary condition implementations that apply external forces to the system."""

import numpy as np from elastica._linalg import _batch_matvec from elastica.typing import SystemType, RodType from elastica.utils import _bspline

from numba import njit from elastica._linalg import _batch_product_i_k_to_ik

class NoForces: """ This is the base class for external forcing boundary conditions applied to rod-like objects.

Notes
-----
Every new external forcing class must be derived
from NoForces class.

"""

def __init__(self):
    """
    NoForces class does not need any input parameters.
    """
    pass

def apply_forces(self, system: SystemType, time: np.float64 = 0.0):
    """Apply forces to a rod-like object.

    In NoForces class, this routine simply passes.

    Parameters
    ----------
    system : SystemType
        Rod or rigid-body object
    time : float
        The time of simulation.

    Returns
    -------

    """
    pass

def apply_torques(self, system: SystemType, time: np.float64 = 0.0):
    """Apply torques to a rod-like object.

    In NoForces class, this routine simply passes.

    Parameters
    ----------
    system : SystemType
        Rod or rigid-body object
    time : float
        The time of simulation.

    Returns
    -------

    """
    pass

class GravityForces(NoForces): """ This class applies a constant gravitational force to the entire rod.

    Attributes
    ----------
    acc_gravity: numpy.ndarray
        1D (dim) array containing data with 'float' type. Gravitational acceleration vector.

"""

def __init__(self, acc_gravity=np.array([0.0, -9.80665, 0.0])):
    """

    Parameters
    ----------
    acc_gravity: numpy.ndarray
        1D (dim) array containing data with 'float' type. Gravitational acceleration vector.

    """
    super(GravityForces, self).__init__()
    self.acc_gravity = acc_gravity

def apply_forces(self, system: SystemType, time=0.0):
    self.compute_gravity_forces(
        self.acc_gravity, system.mass, system.external_forces
    )

@staticmethod
@njit(cache=True)
def compute_gravity_forces(acc_gravity, mass, external_forces):
    """
    This function add gravitational forces on the nodes. We are
    using njit decorated function to increase the speed.

    Parameters
    ----------
    acc_gravity: numpy.ndarray
        1D (dim) array containing data with 'float' type. Gravitational acceleration vector.
    mass: numpy.ndarray
        1D (blocksize) array containing data with 'float' type. Mass on the nodes.
    external_forces: numpy.ndarray
        2D (dim, blocksize) array containing data with 'float' type. External force vector.

    Returns
    -------

    """
    inplace_addition(external_forces, _batch_product_i_k_to_ik(acc_gravity, mass))

class EndpointForces(NoForces): """ This class applies constant forces on the endpoint nodes.

    Attributes
    ----------
    start_force: numpy.ndarray
        1D (dim) array containing data with 'float' type. Force applied to first node of the system.
    end_force: numpy.ndarray
        1D (dim) array containing data with 'float' type. Force applied to last node of the system.
    ramp_up_time: float
        Applied forces are ramped up until ramp up time.

"""

def __init__(self, start_force, end_force, ramp_up_time):
    """

    Parameters
    ----------
    start_force: numpy.ndarray
        1D (dim) array containing data with 'float' type.
        Force applied to first node of the system.
    end_force: numpy.ndarray
        1D (dim) array containing data with 'float' type.
        Force applied to last node of the system.
    ramp_up_time: float
        Applied forces are ramped up until ramp up time.

    """
    super(EndpointForces, self).__init__()
    self.start_force = start_force
    self.end_force = end_force
    assert ramp_up_time > 0.0
    self.ramp_up_time = ramp_up_time

def apply_forces(self, system: SystemType, time=0.0):
    self.compute_end_point_forces(
        system.external_forces,
        self.start_force,
        self.end_force,
        time,
        self.ramp_up_time,
    )

@staticmethod
@njit(cache=True)
def compute_end_point_forces(
    external_forces, start_force, end_force, time, ramp_up_time
):
    """
    Compute end point forces that are applied on the rod using numba njit decorator.

    Parameters
    ----------
    external_forces: numpy.ndarray
        2D (dim, blocksize) array containing data with 'float' type. External force vector.
    start_force: numpy.ndarray
        1D (dim) array containing data with 'float' type.
    end_force: numpy.ndarray
        1D (dim) array containing data with 'float' type.
        Force applied to last node of the system.
    time: float
    ramp_up_time: float
        Applied forces are ramped up until ramp up time.

    Returns
    -------

    """
    factor = min(1.0, time / ramp_up_time)
    external_forces[..., 0] += start_force * factor
    external_forces[..., -1] += end_force * factor

class UniformTorques(NoForces): """ This class applies a uniform torque to the entire rod.

    Attributes
    ----------
    torque: numpy.ndarray
        2D (dim, 1) array containing data with 'float' type. Total torque applied to a rod-like object.

"""

def __init__(self, torque, direction=np.array([0.0, 0.0, 0.0])):
    """

    Parameters
    ----------
    torque: float
        Torque magnitude applied to a rod-like object.
    direction: numpy.ndarray
        1D (dim) array containing data with 'float' type.
        Direction in which torque applied.
    """
    super(UniformTorques, self).__init__()
    self.torque = torque * direction

def apply_torques(self, system: SystemType, time: np.float64 = 0.0):
    n_elems = system.n_elems
    torque_on_one_element = (
        _batch_product_i_k_to_ik(self.torque, np.ones((n_elems))) / n_elems
    )
    system.external_torques += _batch_matvec(
        system.director_collection, torque_on_one_element
    )

class UniformForces(NoForces): """ This class applies a uniform force to the entire rod.

    Attributes
    ----------
    force:  numpy.ndarray
        2D (dim, 1) array containing data with 'float' type. Total force applied to a rod-like object.
"""

def __init__(self, force, direction=np.array([0.0, 0.0, 0.0])):
    """

    Parameters
    ----------
    force: float
        Force magnitude applied to a rod-like object.
    direction: numpy.ndarray
        1D (dim) array containing data with 'float' type.
        Direction in which force applied.
    """
    super(UniformForces, self).__init__()
    self.force = (force * direction).reshape(3, 1)

def apply_forces(self, rod: RodType, time: np.float64 = 0.0):
    force_on_one_element = self.force / rod.n_elems

    rod.external_forces += force_on_one_element

    # Because mass of first and last node is half
    rod.external_forces[..., 0] -= 0.5 * force_on_one_element[:, 0]
    rod.external_forces[..., -1] -= 0.5 * force_on_one_element[:, 0]

class MuscleTorques(NoForces): """ This class applies muscle torques along the body. The applied muscle torques are treated as applied external forces. This class can apply muscle torques as a traveling wave with a beta spline or only as a traveling wave. For implementation details refer to Gazzola et. al. RSoS. (2018).

    Attributes
    ----------
    direction: numpy.ndarray
        2D (dim, 1) array containing data with 'float' type. Muscle torque direction.
    angular_frequency: float
        Angular frequency of traveling wave.
    wave_number: float
        Wave number of traveling wave.
    phase_shift: float
        Phase shift of traveling wave.
    ramp_up_time: float
        Applied muscle torques are ramped up until ramp up time.
    my_spline: numpy.ndarray
        1D (blocksize) array containing data with 'float' type. Generated spline.

"""

def __init__(
    self,
    base_length,
    b_coeff,
    period,
    wave_number,
    phase_shift,
    direction,
    rest_lengths,
    ramp_up_time,
    with_spline=False,
):
    """

    Parameters
    ----------
    base_length: float
        Rest length of the rod-like object.
    b_coeff: nump.ndarray
        1D array containing data with 'float' type.
        Beta coefficients for beta-spline.
    period: float
        Period of traveling wave.
    wave_number: float
        Wave number of traveling wave.
    phase_shift: float
        Phase shift of traveling wave.
    direction: numpy.ndarray
       1D (dim) array containing data with 'float' type. Muscle torque direction.
    ramp_up_time: float
        Applied muscle torques are ramped up until ramp up time.
    with_spline: boolean
        Option to use beta-spline.

    """
    super(MuscleTorques, self).__init__()

    self.direction = direction  # Direction torque applied
    self.angular_frequency = 2.0 * np.pi / period
    self.wave_number = wave_number
    self.phase_shift = phase_shift

    assert ramp_up_time > 0.0
    self.ramp_up_time = ramp_up_time

    # s is the position of nodes on the rod, we go from node=1 to node=nelem-1, because there is no
    # torques applied by first and last node on elements. Reason is that we cannot apply torque in an
    # infinitesimal segment at the beginning and end of rod, because there is no additional element
    # (at element=-1 or element=n_elem+1) to provide internal torques to cancel out an external
    # torque. This coupled with the requirement that the sum of all muscle torques has
    # to be zero results in this condition.
    self.s = np.cumsum(rest_lengths)
    self.s /= self.s[-1]

    if with_spline:
        assert b_coeff.size != 0, "Beta spline coefficient array (t_coeff) is empty"
        my_spline, ctr_pts, ctr_coeffs = _bspline(b_coeff)
        self.my_spline = my_spline(self.s)

    else:

        def constant_function(input):
            """
            Return array of ones same as the size of the input array. This
            function is called when Beta spline function is not used.

            Parameters
            ----------
            input

            Returns
            -------

            """
            return np.ones(input.shape)

        self.my_spline = constant_function(self.s)

def apply_torques(self, rod: RodType, time: np.float64 = 0.0):
    self.compute_muscle_torques(
        time,
        self.my_spline,
        self.s,
        self.angular_frequency,
        self.wave_number,
        self.phase_shift,
        self.ramp_up_time,
        self.direction,
        rod.director_collection,
        rod.external_torques,
    )

@staticmethod
@njit(cache=True)
def compute_muscle_torques(
    time,
    my_spline,
    s,
    angular_frequency,
    wave_number,
    phase_shift,
    ramp_up_time,
    direction,
    director_collection,
    external_torques,
):
    # Ramp up the muscle torque
    factor = min(1.0, time / ramp_up_time)
    # From the node 1 to node nelem-1
    # Magnitude of the torque. Am = beta(s) * sin(2pi*t/T + 2pi*s/lambda + phi)
    # There is an inconsistency with paper and Elastica cpp implementation. In paper sign in
    # front of wave number is positive, in Elastica cpp it is negative.
    torque_mag = (
        factor
        * my_spline
        * np.sin(angular_frequency * time - wave_number * s + phase_shift)
    )
    # Head and tail of the snake is opposite compared to elastica cpp. We need to iterate torque_mag
    # from last to first element.
    torque = _batch_product_i_k_to_ik(direction, torque_mag[::-1])
    inplace_addition(
        external_torques[..., 1:],
        _batch_matvec(director_collection, torque)[..., 1:],
    )
    inplace_substraction(
        external_torques[..., :-1],
        _batch_matvec(director_collection[..., :-1], torque[..., 1:]),
    )

class LiftingTorques(NoForces): """ This class applies Lifting torques along the body. The applied muscle torques are treated as applied external forces. This class can apply muscle torques as a traveling wave with a beta spline or only as a traveling wave. For implementation details refer to Gazzola et. al. RSoS. (2018).

    Attributes
    ----------
    direction: numpy.ndarray
        2D (dim, 1) array containing data with 'float' type. Muscle torque direction.
    angular_frequency: float
        Angular frequency of traveling wave.
    wave_number: float
        Wave number of traveling wave.
    phase_shift: float
        Phase shift of traveling wave.
    ramp_up_time: float
        Applied muscle torques are ramped up until ramp up time.
    my_spline: numpy.ndarray
        1D (blocksize) array containing data with 'float' type. Generated spline.

"""

def __init__(
        self,
        base_length,
        b_coeff,
        period,
        wave_number,
        Amplitude,
        lamBDA,
        phase_shift,
        direction,
        rest_lengths,
        ramp_up_time,
        with_spline=False,
):
    """

    Parameters
    ----------
    base_length: float
        Rest length of the rod-like object.
    b_coeff: nump.ndarray
        1D array containing data with 'float' type.
        Beta coefficients for beta-spline.
    period: float
        Period of traveling wave.
    wave_number: float
        Wave number of traveling wave.
    Amplitude : flot
    lamBDA : float
    phase_shift: float
        Phase shift of traveling wave.
    direction: numpy.ndarray
       1D (dim) array containing data with 'float' type. Muscle torque direction.
    ramp_up_time: float
        Applied muscle torques are ramped up until ramp up time.
    with_spline: boolean
        Option to use beta-spline.

    """
    super(LiftingTorques, self).__init__()

    self.direction = direction  # Direction torque applied
    self.angular_frequency = 2.0 * np.pi / period
    self.wave_number = wave_number
    self.phase_shift = phase_shift

    assert ramp_up_time > 0.0
    self.ramp_up_time = ramp_up_time

    # s is the position of nodes on the rod, we go from node=1 to node=nelem-1, because there is no
    # torques applied by first and last node on elements. Reason is that we cannot apply torque in an
    # infinitesimal segment at the beginning and end of rod, because there is no additional element
    # (at element=-1 or element=n_elem+1) to provide internal torques to cancel out an external
    # torque. This coupled with the requirement that the sum of all muscle torques has
    # to be zero results in this condition.
    self.s = np.cumsum(rest_lengths)
    self.s /= self.s[-1]

    if with_spline:
        assert b_coeff.size != 0, "Beta spline coefficient array (t_coeff) is empty"
        my_spline, ctr_pts, ctr_coeffs = _bspline(b_coeff)
        self.my_spline = my_spline(self.s)

    else:

        def constant_function(input):
            """
            Return array of ones same as the size of the input array. This
            function is called when Beta spline function is not used.

            Parameters
            ----------
            input

            Returns
            -------

            """
            return np.ones(input.shape)

        self.my_spline = constant_function(self.s)

def apply__lifting_torques(self, rod: RodType, time: np.float64 = 0.0):
    self.compute_lifting_torques(
        time,
        self.my_spline,
        self.s,
        self.angular_frequency,
        self.wave_number,
        self.Amplitude,
        self.lamBDA,
        self.phase_shift,
        self.ramp_up_time,
        self.direction,
        rod.director_collection,
        rod.external_torques,
    )

def compute_lifting_torques(
        time,
        my_spline,
        s,
        angular_frequency,
        wave_number,  # k pf muscle torque
        Amplitude,
        lamBDA,
        phase_shift,
        ramp_up_time,
        direction,
        director_collection,
        external_torques,
):
    # Ramp up the muscle torque
    factor = min(1.0, time / ramp_up_time)
    # From the node 1 to node nelem-1
    # Magnitude of the torque. Am = beta(s) * sin(2pi*t/T + 2pi*s/lambda + phi)
    # There is an inconsistency with paper and Elastica cpp implementation. In paper sign in
    # front of wave number is positive, in Elastica cpp it is negative.
    # torque_mag = (
    #    factor
    #    * my_spline
    #    * np.sin(angular_frequency * time - wave_number * s + phase_shift)
    # )
    kt = wave_number * lamBDA
    torque_mag = max(0, (
            factor
            * my_spline
            * Amplitude * np.cos(angular_frequency * time + kt * s + phase_shift) + 1
    ))
    # Head and tail of the snake is opposite compared to elastica cpp. We need to iterate torque_mag
    # from last to first element.
    torque = _batch_product_i_k_to_ik(direction, torque_mag[::-1])
    inplace_addition(
        external_torques[..., 1:],
        _batch_matvec(director_collection, torque)[..., 1:],
    )
    inplace_substraction(
        external_torques[..., :-1],
        _batch_matvec(director_collection[..., :-1], torque[..., 1:]),
    )

@njit(cache=True) def inplace_addition(external_force_or_torque, force_or_torque): """ This function does inplace addition. First argument external_force_or_torque is the system.external_forces or system.external_torques. Second argument force or torque vector to be added.

Parameters
----------
external_force_or_torque: numpy.ndarray
    2D (dim, blocksize) array containing data with 'float' type.
force_or_torque: numpy.ndarray
    2D (dim, blocksize) array containing data with 'float' type.

Returns
-------

"""
blocksize = force_or_torque.shape[1]
for i in range(3):
    for k in range(blocksize):
        external_force_or_torque[i, k] += force_or_torque[i, k]

@njit(cache=True) def inplace_substraction(external_force_or_torque, force_or_torque): """ This function does inplace substraction. First argument external_force_or_torque is the system.external_forces or system.external_torques. Second argument force or torque vector to be substracted. Parameters

external_force_or_torque: numpy.ndarray
    2D (dim, blocksize) array containing data with 'float' type.
force_or_torque: numpy.ndarray
    2D (dim, blocksize) array containing data with 'float' type.

Returns
-------

"""
blocksize = force_or_torque.shape[1]
for i in range(3):
    for k in range(blocksize):
        external_force_or_torque[i, k] -= force_or_torque[i, k]

class EndpointForcesSinusoidal(NoForces): """ This class applies sinusoidally varying forces to the ends of a rod. Forces are applied in a plane, which is defined by the tangent_direction and normal_direction.

    Attributes
    ----------
    start_force_mag: float
        Magnitude of the force that is applied to the start of the rod (node 0).
    end_force_mag: float
        Magnitude of the force that is applied to the end of the rod (node -1).
    ramp_up_time: float
        Applied forces are applied in the normal direction until time reaches ramp_up_time.
    normal_direction: np.ndarray
        An array (3,) contains type float.
        This is the normal direction of the rod.
    roll_direction: np.ndarray
        An array (3,) contains type float.
        This is the direction perpendicular to rod tangent, and rod normal.

    Notes
    -----
    In order to see example how to use this class, see joint examples.

"""

def __init__(
    self,
    start_force_mag,
    end_force_mag,
    ramp_up_time=0.0,
    tangent_direction=np.array([0, 0, 1]),
    normal_direction=np.array([0, 1, 0]),
):
    """

    Parameters
    ----------
    start_force_mag: float
        Magnitude of the force that is applied to the start of the system (node 0).
    end_force_mag: float
        Magnitude of the force that is applied to the end of the system (node -1).
    ramp_up_time: float
        Applied muscle torques are ramped up until ramp up time.
    tangent_direction: np.ndarray
        An array (3,) contains type float.
        This is the tangent direction of the system, or normal of the plane that forces applied.
    normal_direction: np.ndarray
        An array (3,) contains type float.
        This is the normal direction of the system.
    """
    super(EndpointForcesSinusoidal, self).__init__()
    # Start force
    self.start_force_mag = start_force_mag
    self.end_force_mag = end_force_mag

    # Applied force directions
    self.normal_direction = normal_direction
    self.roll_direction = np.cross(normal_direction, tangent_direction)

    assert ramp_up_time >= 0.0
    self.ramp_up_time = ramp_up_time

def apply_forces(self, system: SystemType, time=0.0):

    if time < self.ramp_up_time:
        # When time smaller than ramp up time apply the force in normal direction
        # First pull the rod upward or downward direction some time.
        start_force = -2.0 * self.start_force_mag * self.normal_direction
        end_force = -2.0 * self.end_force_mag * self.normal_direction

        system.external_forces[..., 0] += start_force
        system.external_forces[..., -1] += end_force

    else:
        # When time is greater than ramp up time, forces are applied in normal
        # and roll direction or forces are in a plane perpendicular to the
        # direction.

        # First force applied to start of the rod
        roll_forces_start = (
            self.start_force_mag
            * np.cos(0.5 * np.pi * (time - self.ramp_up_time))
            * self.roll_direction
        )
        normal_forces_start = (
            self.start_force_mag
            * np.sin(0.5 * np.pi * (time - self.ramp_up_time))
            * self.normal_direction
        )
        start_force = roll_forces_start + normal_forces_start
        # Now force applied to end of the rod
        roll_forces_end = (
            self.end_force_mag
            * np.cos(0.5 * np.pi * (time - self.ramp_up_time))
            * self.roll_direction
        )
        normal_forces_end = (
            self.end_force_mag
            * np.sin(0.5 * np.pi * (time - self.ramp_up_time))
            * self.normal_direction
        )
        end_force = roll_forces_end + normal_forces_end
        # Update external forces
        system.external_forces[..., 0] += start_force
        system.external_forces[..., -1] += end_force

here's how I call it in my simulation file

import os import numpy as np import elastica as ea

Apply Muscle force related parameters

Helper Function to check if a point is inside a patch

def is_point_inside_patch(point, patch_vertices): """ Check if a 3D point is inside a patch defined by four vertices.

Args:
    point (numpy.ndarray): 3D coordinates of the point.
    patch_vertices (list of numpy.ndarray): List of four 3D coordinates defining the patch.

Returns:
    bool: True if the point is inside the patch, False otherwise.
"""
# Calculate the normal vectors of the four triangles formed by the point and patch vertices
normals = []
for i in range(4):
    vertex1 = patch_vertices[i]
    vertex2 = patch_vertices[(i + 1) % 4]
    normal = np.cross(vertex2 - vertex1, point - vertex1)
    normals.append(normal)

# Check if the point is on the same side of all the patch's triangles (using dot product)
for i in range(4):
    if np.dot(normals[i], normals[(i + 1) % 4]) < 0:
        return False

return True

""" from snake_diffraction_postprocessing_3D import ( plot_snake_velocity, plot_video_3d, compute_projected_velocity, plot_curvature, ) """

class SnakeDiffractionSimulator( ea.BaseSystemCollection, ea.Constraints, ea.Forcing, ea.Damping, ea.CallBacks ): pass

def run_snake( b_coeff, PLOT_FIGURE=False, SAVE_FIGURE=False, SAVE_VIDEO=False, SAVE_RESULTS=False ):

Initialize the simulation class

snake_diffraction_sim = SnakeDiffractionSimulator()

# Simulation parameters
period = 15
final_time = (2.0 + 0.01) * period

# setting up test params
n_elem = 50
start = np.zeros((3,))
direction = np.array([0.0, 0.0, 1.0])
#lateral_direction = np.array([1.0,0.0,0.0])
normal = np.array([0.0, 1.0, 0.0])
base_length = 0.35
base_radius = base_length * 0.011
density = 1000
E = 1e6
poisson_ratio = 0.5
shear_modulus = E / (poisson_ratio + 1.0)

snake_body = ea.CosseratRod.straight_rod(
    n_elem,
    start,
    direction,
    normal,
    base_length,
    base_radius,
    density,
    youngs_modulus=E,
    shear_modulus=shear_modulus,
)

snake_diffraction_sim.append(snake_body)

# Add gravitational forces
gravitational_acc = -9.80665
snake_diffraction_sim.add_forcing_to(snake_body).using(
    ea.GravityForces, acc_gravity=np.array([0.0, gravitational_acc, 0.0])
)

# Add muscle torques only Lateral Wave right now
wave_length = b_coeff[-1]
snake_diffraction_sim.add_forcing_to(snake_body).using(
    ea.MuscleTorques,
    base_length=base_length,
    b_coeff=b_coeff[:-1],
    period=period,
    wave_number=2.0 * np.pi / (wave_length),
    phase_shift=0.0,
    rest_lengths=snake_body.rest_lengths,
    ramp_up_time=period,
    direction=normal,
    with_spline=True,
)

# Add Lifting Wave 
snake_diffraction_sim.add_forcing_to(snake_body).using(
  ea.LiftingTorques,
  base_length=base_length,
  b_coeff=b_coeff[:-1],
  period=period,
  wave_number=2.0 * np.pi / (wave_length),
  Amplitude=0.7,
  lamBDA=1,
  phase_shift=0.25,
  rest_lengths=snake_body.rest_lengths,
  ramp_up_time=period,
  direction=direction,
  with_spline=True,
)

# Add friction forces
# Uniform friction with ground
origin_plane = np.array([0.0, -base_radius, 0.0])
normal_plane = normal
slip_velocity_tol_Regular = 1e-6
slip_velocity_tol_HighFriction = 1e-4
p_factor = 30
froude = 0.1
mu = base_length / (period * period * np.abs(gravitational_acc) * froude)
kinetic_mu_array_Regular = np.array(
    [mu, 2 * mu, 10.0 * mu]
)  # [forward, backward, sideways]

static_mu_array_Regular = np.zeros(kinetic_mu_array_Regular.shape)

kinetic_mu_array_HighFriction = np.array(
    [mu * p_factor, 2 * mu * p_factor, 10.0 * mu * p_factor]
)

static_mu_array_HighFriction = np.zeros(kinetic_mu_array_HighFriction.shape)

patch1_A_coord = np.array([0.15, 0.0, 0.6])
patch1_B_coord = np.array([0.15, 0.0, 0.9])
patch1_C_coord = np.array([-0.15, 0.0, 0.9])
patch1_D_coord = np.array([-0.15, 0.0, 0.6])

patch_vertices = [
    patch1_A_coord,
    patch1_B_coord,
    patch1_C_coord,
    patch1_D_coord,
]

#print("Patch Vertices Are.")
#print(patch_vertices)

snake_positions = snake_body.position_collection
#print("snake position is...")
#print(snake_positions)
#print("coordinates of 23 elements")
#print(snake_positions[:, 23])

"""
# Add frictional forces to the snake based on the region of space
for i in range(n_elem):
    x_coord = snake_positions[0, i]  # x-coordinate of the element
    y_coord = snake_positions[1, i]  # y-coordinate of the element
    z_coord = snake_positions[2, i]  # z-coordinate of the element

    element_position = np.array([x_coord, y_coord, z_coord])  # Create a 3D position vector

    if is_point_inside_patch(element_position, patch_vertices):
        # Apply high friction forces
        mu = kinetic_mu_array_HighFriction
        static_mu = static_mu_array_HighFriction
        slip_velocity_tol = slip_velocity_tol_HighFriction
    else:
        # Apply regular friction forces
        mu = kinetic_mu_array_Regular
        static_mu = static_mu_array_Regular
        slip_velocity_tol = slip_velocity_tol_Regular

        snake_diffraction_sim.add_forcing_to(snake_body).using(
            ea.AnisotropicFrictionalPlane,
            k=1.0,
            nu=1e-6,
            plane_origin=origin_plane,
            plane_normal=normal_plane,
            slip_velocity_tol=slip_velocity_tol,
            static_mu_array=static_mu,
            kinetic_mu_array=mu,
        )
        """

# Add uniform frictional force to snake
snake_diffraction_sim.add_forcing_to(snake_body).using(
    ea.AnisotropicFrictionalPlane,
    k=1.0,
    nu=1e-6,
    plane_origin=origin_plane,
    plane_normal=normal_plane,
    slip_velocity_tol=slip_velocity_tol_Regular,
    static_mu_array=static_mu_array_Regular,
    kinetic_mu_array=kinetic_mu_array_Regular,
)

# add damping
damping_constant = 2e-3
time_step = 1e-4
snake_diffraction_sim.dampen(snake_body).using(
    ea.AnalyticalLinearDamper,
    damping_constant=damping_constant,
    time_step=time_step,
)

total_steps = int(final_time / time_step)
rendering_fps = 60
step_skip = int(1.0 / (rendering_fps * time_step))

# Add call backs
class SnakeDiffractionCallBack(ea.CallBackBaseClass):
    """
    Call back function for continuum snake
    """

    def __init__(self, step_skip: int, callback_params: dict):
        ea.CallBackBaseClass.__init__(self)
        self.every = step_skip
        self.callback_params = callback_params

    def make_callback(self, system, time, current_step: int):
        if current_step % self.every == 0:
            self.callback_params["time"].append(time)
            self.callback_params["step"].append(current_step)
            self.callback_params["position"].append(
                system.position_collection.copy()
            )
            self.callback_params["velocity"].append(
                system.velocity_collection.copy()
            )
            self.callback_params["avg_velocity"].append(
                system.compute_velocity_center_of_mass()
            )

            self.callback_params["center_of_mass"].append(
                system.compute_position_center_of_mass()
            )
            self.callback_params["curvature"].append(system.kappa.copy())

            return

pp_list = ea.defaultdict(list)
snake_diffraction_sim.collect_diagnostics(snake_body).using(
    SnakeDiffractionCallBack, step_skip=step_skip, callback_params=pp_list
)

snake_diffraction_sim.finalize()

timestepper = ea.PositionVerlet()
ea.integrate(timestepper, snake_diffraction_sim, final_time, total_steps)

if PLOT_FIGURE:
    filename_plot = "snake_diffraction_velocity_test6.png"
    plot_snake_velocity(pp_list, period, filename_plot, SAVE_FIGURE)
    plot_curvature(pp_list, snake_body.rest_lengths, period, SAVE_FIGURE)

    if SAVE_VIDEO:
        filename_video = 'snake_diffraction_test6.gif'

        # filename_video = "continuum_snake.mp4"
        plot_video(
            pp_list,
            video_name=filename_video,

            xlim=(0, 4),
            ylim=(-1, 1),
        )

        #filename_video = 'snake_diffraction_3d.gif'
        #plot_video_3d(
        #    [pp_list],
        #    video_name=filename_video,

        #    fps=rendering_fps,
        #    xlim=(-3, 3),
        #    ylim=(-3, 3),
         #   zlim=(-3, 3)
        #)

if SAVE_RESULTS:
    import pickle

    filename = "snake_diffraction_test6.dat"
    file = open(filename, "wb")
    pickle.dump(pp_list, file)
    file.close()

# Compute the average forward velocity. These will be used for optimization.
[_, _, avg_forward, avg_lateral] = compute_projected_velocity(pp_list, period)

return avg_forward, avg_lateral, pp_list

if name == "main":

# Options
PLOT_FIGURE = True
SAVE_FIGURE = True
SAVE_VIDEO = True
SAVE_RESULTS = False
CMA_OPTION = False

if CMA_OPTION:
    import cma

    SAVE_OPTIMIZED_COEFFICIENTS = False

    def optimize_snake(spline_coefficient):
        [avg_forward, _, _] = run_snake(
            spline_coefficient,
            PLOT_FIGURE=False,
            SAVE_FIGURE=False,
            SAVE_VIDEO=False,
            SAVE_RESULTS=False,
        )
        return -avg_forward

    # Optimize snake for forward velocity. In cma.fmin first input is function
    # to be optimized, second input is initial guess for coefficients you are optimizing
    # for and third input is standard deviation you initially set.
    optimized_spline_coefficients = cma.fmin(optimize_snake, 7 * [0], 0.5)

    # Save the optimized coefficients to a file
    filename_data = "optimized_coefficients.txt"
    if SAVE_OPTIMIZED_COEFFICIENTS:
        assert filename_data != "", "provide a file name for coefficients"
        np.savetxt(filename_data, optimized_spline_coefficients, delimiter=",")

else:
    # Add muscle forces on the rod
    if os.path.exists("optimized_coefficients.txt"):
        t_coeff_optimized = np.genfromtxt(
            "optimized_coefficients.txt", delimiter=","
        )
    else:
        wave_length = 1.0
        t_coeff_optimized = np.array(
            [3.4e-3, 3.3e-3, 4.2e-3, 2.6e-3, 3.6e-3, 3.5e-3]
        )
        t_coeff_optimized = np.hstack((t_coeff_optimized, wave_length))

    # run the simulation
    [avg_forward, avg_lateral, pp_list] = run_snake(
        t_coeff_optimized, PLOT_FIGURE, SAVE_FIGURE, SAVE_VIDEO, SAVE_RESULTS
    )

    print("average forward velocity:", avg_forward)
    print("average forward lateral:", avg_lateral)
Ali-7800 commented 12 months ago

Hi @AkashVardhan7 ,

In elastica we have an __init__.py file which is used for importing different functions/classes, if you take a look at it you can see that we are importing all the classes in external_forces.py like this:

from elastica.external_forces import ( NoForces, EndpointForces, GravityForces, UniformForces, UniformTorques, MuscleTorques, EndpointForcesSinusoidal, )

if you want use the class you just created the way you are doing it now, you have to add it there too. Alternatively, you can put your class in the simulation file with all the necessary imports and use it without the ea.. One last thing, in order for this class to work it must also have the method apply_torques so it can be called in the backend, right now you have it called apply__lifting_torques which would not work.

Let me know if you need more help!

AkashVardhan7 commented 12 months ago

I modified the init file to have the LiftingTorques in from elastica.external_forces import(),. and modified the def apply_torques(self, rod: RodType, time: np.float64 = 0.0): as you suggested in the external_forces.py... the error I get now has to di with the forcing.py... Here are the modified code snippets and the error message...

from elastica.external_forces import ( NoForces, EndpointForces, GravityForces, UniformForces, UniformTorques, MuscleTorques, LiftingTorques, EndpointForcesSinusoidal, )

def apply_torques(self, rod: RodType, time: np.float64 = 0.0): self.compute_lifting_torques( time, self.my_spline, self.s, self.angular_frequency, self.wave_number, self.Amplitude, self.lamBDA, self.phase_shift, self.ramp_up_time, self.direction, rod.director_collection, rod.external_torques, ) def apply_torques(self, rod: RodType, time: np.float64 = 0.0): self.compute_lifting_torques( time, self.my_spline, self.s, self.angular_frequency, self.wave_number, self.Amplitude, self.lamBDA, self.phase_shift, self.ramp_up_time, self.direction, rod.director_collection, rod.external_torques, )

/Users/avardhan7/Documents/OPEN_AI/pyElastica/snake_diffraction/snake_diffraction.py Traceback (most recent call last): File "C:\Users\avardhan7.conda\envs\pyELastica\lib\site-packages\elastica\modules\forcing.py", line 171, in call return self._forcing_cls(*self._args, **self._kwargs) TypeError: NoForces.init() got an unexpected keyword argument 'base_length'

During handling of the above exception, another exception occurred:

Traceback (most recent call last): File "c:\Users\avardhan7\Documents\OPEN_AI\pyElastica\snake_diffraction\snake_diffraction.py", line 379, in [avg_forward, avg_lateral, pp_list] = run_snake( File "c:\Users\avardhan7\Documents\OPEN_AI\pyElastica\snake_diffraction\snake_diffraction.py", line 280, in run_snake
snake_diffraction_sim.finalize() File "C:\Users\avardhan7.conda\envs\pyELastica\lib\site-packages\elastica\modules\base_system.py", line 163, in finalize
finalize() File "C:\Users\avardhan7.conda\envs\pyELastica\lib\site-packages\elastica\modules\forcing.py", line 59, in _finalize_forcing self._ext_forces_torques[:] = [ File "C:\Users\avardhan7.conda\envs\pyELastica\lib\site-packages\elastica\modules\forcing.py", line 60, in
(ext_force_torque.id(), ext_force_torque()) File "C:\Users\avardhan7.conda\envs\pyELastica\lib\site-packages\elastica\modules\forcing.py", line 173, in call
raise TypeError( TypeError: Unable to construct forcing class.\nDid you provide all necessary force properties? PS C:\Users\avardhan7\Documents\OPEN_AI\pyElastica\snake_diffraction>

Ali-7800 commented 12 months ago

It looks like the error is that you are passing base_length to a forcing class that does not take base_length, did you make any other changes to your code other than I told you about?

AkashVardhan7 commented 11 months ago

I didn't. The only changes I made are the following.

  1. Put the class LiftingTorques inside external_forces.py {when I tried running it from the main simulation after appropriately pasting all the necessary imports and also the class NoForces.py... I got the following message... Assertion Error: <class 'main.LiftingTorques'> is not a valid forcing. Did you forget to derive from NoForces?....so went back to running with having the class in external_forces.py}.

I am pasting this class again., which I modified from the class MuscleTorques...base_length is passed as an input to definit and then in the main simulation it is called by the LiftingTorques class by passing the same base_length value that is passed to the class MuscleTorques. class LiftingTorques(NoForces): """ This class applies Lifting torques along the body. The applied muscle torques are treated as applied external forces. This class can apply muscle torques as a traveling wave with a beta spline or only as a traveling wave. For implementation details refer to Gazzola et. al. RSoS. (2018).

Attributes
----------
direction: numpy.ndarray
    2D (dim, 1) array containing data with 'float' type. Muscle torque direction.
angular_frequency: float
    Angular frequency of traveling wave.
wave_number: float
    Wave number of traveling wave.
phase_shift: float
    Phase shift of traveling wave.
ramp_up_time: float
    Applied muscle torques are ramped up until ramp up time.
my_spline: numpy.ndarray
    1D (blocksize) array containing data with 'float' type. Generated spline.

"""

def init( self, base_length, b_coeff, period, wave_number, Amplitude, lamBDA, phase_shift, direction, rest_lengths, ramp_up_time, with_spline=False, ): """

Parameters
----------
base_length: float
    Rest length of the rod-like object.
b_coeff: nump.ndarray
    1D array containing data with 'float' type.
    Beta coefficients for beta-spline.
period: float
    Period of traveling wave.
wave_number: float
    Wave number of traveling wave.
Amplitude : flot
lamBDA : float
phase_shift: float
    Phase shift of traveling wave.
direction: numpy.ndarray
   1D (dim) array containing data with 'float' type. Muscle torque direction.
ramp_up_time: float
    Applied muscle torques are ramped up until ramp up time.
with_spline: boolean
    Option to use beta-spline.

"""
super(LiftingTorques, self).__init__()

self.direction = direction  # Direction torque applied
self.angular_frequency = 2.0 * np.pi / period
self.wave_number = wave_number
self.phase_shift = phase_shift

assert ramp_up_time > 0.0
self.ramp_up_time = ramp_up_time

# s is the position of nodes on the rod, we go from node=1 to node=nelem-1, because there is no
# torques applied by first and last node on elements. Reason is that we cannot apply torque in an
# infinitesimal segment at the beginning and end of rod, because there is no additional element
# (at element=-1 or element=n_elem+1) to provide internal torques to cancel out an external
# torque. This coupled with the requirement that the sum of all muscle torques has
# to be zero results in this condition.
self.s = np.cumsum(rest_lengths)
self.s /= self.s[-1]

if with_spline:
    assert b_coeff.size != 0, "Beta spline coefficient array (t_coeff) is empty"
    my_spline, ctr_pts, ctr_coeffs = _bspline(b_coeff)
    self.my_spline = my_spline(self.s)

else:

    def constant_function(input):
        """
        Return array of ones same as the size of the input array. This
        function is called when Beta spline function is not used.

        Parameters
        ----------
        input

        Returns
        -------

        """
        return np.ones(input.shape)

    self.my_spline = constant_function(self.s)

def apply_torques(self, rod: RodType, time: np.float64 = 0.0): self.compute_lifting_torques( time, self.my_spline, self.s, self.angular_frequency, self.wave_number, self.Amplitude, self.lamBDA, self.phase_shift, self.ramp_up_time, self.direction, rod.director_collection, rod.external_torques, )

def compute_lifting_torques( time, my_spline, s, angular_frequency, wave_number, # k pf muscle torque Amplitude, lamBDA, phase_shift, ramp_up_time, direction, director_collection, external_torques, ):

Ramp up the muscle torque

factor = min(1.0, time / ramp_up_time)
# From the node 1 to node nelem-1
# Magnitude of the torque. Am = beta(s) * sin(2pi*t/T + 2pi*s/lambda + phi)
# There is an inconsistency with paper and Elastica cpp implementation. In paper sign in
# front of wave number is positive, in Elastica cpp it is negative.
# torque_mag = (
#    factor
#    * my_spline
#    * np.sin(angular_frequency * time - wave_number * s + phase_shift)
# )
kt = wave_number * lamBDA
torque_mag = max(0, (
        factor
        * my_spline
        * Amplitude * np.cos(angular_frequency * time + kt * s + phase_shift) + 1
))
# Head and tail of the snake is opposite compared to elastica cpp. We need to iterate torque_mag
# from last to first element.
torque = _batch_product_i_k_to_ik(direction, torque_mag[::-1])
inplace_addition(
    external_torques[..., 1:],
    _batch_matvec(director_collection, torque)[..., 1:],
)
inplace_substraction(
    external_torques[..., :-1],
    _batch_matvec(director_collection[..., :-1], torque[..., 1:]),
)

Here's how I call the class in the main simulation...
 snake_diffraction_sim.add_forcing_to(snake_body).using(
  ea.LiftingTorques,
  base_length=base_length,
  b_coeff=b_coeff[:-1],
  period=period,
  wave_number=2.0 * np.pi / (wave_length),
  Amplitude=0.7,
  lamBDA=1,
  phase_shift=0.25,
  rest_lengths=snake_body.rest_lengths,
  ramp_up_time=period,
  direction=direction,
  with_spline=True,
)

and here' the addition of the class to the __init__.py file..
from elastica.external_forces import (
NoForces,
EndpointForces,
GravityForces,
UniformForces,
UniformTorques,
MuscleTorques,
LiftingTorques,
EndpointForcesSinusoidal,

)

Besides, these I made no changes to any of the files..

Ali-7800 commented 11 months ago

I think the problem is with the indentation of what is inside of your class, make sure that you indent everything in the class correctly. Let me know how it goes.

AkashVardhan7 commented 11 months ago

There was an indentation error, which I addressed, but now I am getting the following message,

Traceback (most recent call last): File "c:\Users\avardhan7\Documents\OPEN_AI\pyElastica\snake_diffraction\snake_diffraction.py", line 379, in [avg_forward, avg_lateral, pp_list] = run_snake( File "c:\Users\avardhan7\Documents\OPEN_AI\pyElastica\snake_diffraction\snake_diffraction.py", line 283, in run_snake ea.integrate(timestepper, snake_diffraction_sim, final_time, total_steps) File "C:\Users\avardhan7.conda\envs\pyELastica\lib\site-packages\elastica\timestepper__init__.py", line 106, in integrate time = do_step(StatefulStepper, stages_and_updates, System, time, dt) File "C:\Users\avardhan7.conda\envs\pyELastica\lib\site-packages\elastica\timestepper\symplectic_steppers.py", line 90, in do_step SystemCollection.synchronize(time) File "C:\Users\avardhan7.conda\envs\pyELastica\lib\site-packages\elastica\modules\base_system.py", line 175, in synchronize feature(time) File "C:\Users\avardhan7.conda\envs\pyELastica\lib\site-packages\elastica\modules\forcing.py", line 86, in _call_ext_forces_torques ext_force_torque.apply_torques(self._systems[sys_id], time, *args, **kwargs) File "C:\Users\avardhan7.conda\envs\pyELastica\lib\site-packages\elastica\external_forces.py", line 483, in apply_torques self.compute_lifting_torques( TypeError: LiftingTorques.compute_lifting_torques() takes 12 positional arguments but 13 were given PS C:\Users\avardhan7\Documents\OPEN_AI\pyElastica\snake_diffraction>

This is my current class for applying LiftingTorques

class LiftingTorques(NoForces): """ This class applies Lifting torques along the body. The applied muscle torques are treated as applied external forces. This class can apply muscle torques as a traveling wave with a beta spline or only as a traveling wave. For implementation details refer to Gazzola et. al. RSoS. (2018). """

def __init__(
    self,
    base_length,
    b_coeff,
    period,
    wave_number,
    Amplitude,
    lamBDA,
    phase_shift,
    direction,
    rest_lengths,
    ramp_up_time,
    with_spline=False
):
    super(LiftingTorques, self).__init__()

    self.direction = direction  # Direction torque applied
    self.Amplitude = Amplitude
    self.lamBDA = lamBDA
    self.angular_frequency = 2.0 * np.pi / period
    self.wave_number = wave_number
    self.phase_shift = phase_shift

    assert ramp_up_time > 0.0
    self.ramp_up_time = ramp_up_time

    # s is the position of nodes on the rod, we go from node=1 to node=nelem-1, because there is no
    # torques applied by first and last node on elements. Reason is that we cannot apply torque in an
    # infinitesimal segment at the beginning and end of rod, because there is no additional element
    # (at element=-1 or element=n_elem+1) to provide internal torques to cancel out an external
    # torque. This coupled with the requirement that the sum of all muscle torques has
    # to be zero results in this condition.
    self.s = np.cumsum(rest_lengths)
    self.s /= self.s[-1]

    if with_spline:
        assert b_coeff.size != 0, "Beta spline coefficient array (t_coeff) is empty"
        my_spline, ctr_pts, ctr_coeffs = _bspline(b_coeff)
        self.my_spline = my_spline(self.s)
    else:
        def constant_function(input):
            """
            Return array of ones same as the size of the input array. This
            function is called when Beta spline function is not used.

            Parameters
            ----------
            input

            Returns
            -------

            """
            return np.ones(input.shape)

        self.my_spline = constant_function(self.s)

def apply_torques(self, rod: RodType, time: np.float64 = 0.0):
    self.compute_lifting_torques(
        time,
        self.my_spline,
        self.s,
        self.angular_frequency,
        self.wave_number,
        self.Amplitude,
        self.lamBDA,
        self.phase_shift,
        self.ramp_up_time,
        self.direction,
        rod.director_collection,
        rod.external_torques,
    )

def compute_lifting_torques(
    time,
    my_spline,
    s,
    angular_frequency,
    wave_number,  # k pf muscle torque
    Amplitude,
    lamBDA,
    phase_shift,
    ramp_up_time,
    direction,
    director_collection,
    external_torques,
):
    # Ramp up the muscle torque
    factor = min(1.0, time / ramp_up_time)

    # From the node 1 to node nelem-1
    # Magnitude of the torque. Am = beta(s) * sin(2pi*t/T + 2pi*s/lambda + phi)
    # There is an inconsistency with paper and Elastica cpp implementation. In paper sign in
    # front of wave number is positive, in Elastica cpp it is negative.
    # torque_mag = (
    #    factor
    #    * my_spline
    #    * np.sin(angular_frequency * time - wave_number * s + phase_shift)
    # )
    kt = wave_number * lamBDA
    torque_mag = max(0, (
        factor
        * my_spline
        * Amplitude * np.cos(angular_frequency * time + kt * s + phase_shift) + 1
    ))

    # Head and tail of the snake is opposite compared to elastica cpp. We need to iterate torque_mag
    # from last to first element.
    torque = _batch_product_i_k_to_ik(direction, torque_mag[::-1])
    inplace_addition(
        external_torques[..., 1:],
        _batch_matvec(director_collection, torque)[..., 1:],
    )
    inplace_substraction(
        external_torques[..., :-1],
        _batch_matvec(director_collection[..., :-1], torque[..., 1:]),
    )
AkashVardhan7 commented 11 months ago

I have retained everything which was present in the class MuscleTorques with the addition of variables lamBDA and Amplitude... which represent those two parameters from the paper

AkashVardhan7 commented 11 months ago

Okay, so I got the simulation to run with the liftingtorques class after fiddling around with the indentation some more and using np.maximum to apply the liftingTorques rather than just max., but there is no snake in the final simulation and all I see is a blank.. when I comment out the (liftingTorques) part in the simulation and run it only with the muscle torques, I see the snake moving along a straight line as before.. clearly I am doing something wrong...

AkashVardhan7 commented 11 months ago

I put my code in a repo and shared it with you, it'd be helpful if you could try running it on your end to see what I am saying.

Ali-7800 commented 11 months ago

@AkashVardhan7 I think the problem is in torque_mag = np.maximum(0,(factor* my_spline* Amplitude* np.cos(angular_frequency * time + kt * s + phase_shift)+ 1 ),), I think it should be torque_mag = np.maximum(0,(factor* my_spline* (Amplitude* np.cos(angular_frequency * time + kt * s + phase_shift)+ 1) ),) instead.

AkashVardhan7 commented 11 months ago

image Ali I changed it to what you suggested, yet the problem persists. :(

armantekinalp commented 11 months ago

@xzhan139 can you also take a look at this issue.

AkashVardhan7 commented 11 months ago

I shared the repo with @xzhan139

AkashVardhan7 commented 11 months ago

Any further pointers on what I should be doing to make this work? I'd want to capture all the behaviors you reported in the manuscript before I can reliably start using this sim for my project...

xzhan139 commented 11 months ago

Hi can you share your repo again? I didn't get any invitation yet

AkashVardhan7 commented 11 months ago

https://github.com/AkashVardhan7/snake_diffraction_Elastica.git

Here's the link to the repo (I made it public so you should be able to access the code... I made modifications to externalforces.py and init.py as you can follow from the thread.. here's a snapshot of the invitation I sent you image

xzhan139 commented 11 months ago

Thanks for sharing! Sorry I missed your invitation on Saturday. A couple of things -

  1. In your lifting force Screenshot 2023-10-23 at 2 02 31 PM

    If I didn't misunderstand, "direction = direction" is twist direction, not lifting direction. Adding twist in general is very tricky as you have to be very careful with the torque amplitude and timestep, otherwise your simulation diverge very quickly. I guess that is why you don't see the snake, it exploded. Lifting is instead in binormal direction.

  2. When you are adding new force, it's generally recommended to start with a very small force value, just to test the algorithm implementation itself. Big force itself can make simulation unstable, so it's harder to figure out where the problem is.
  3. An alternative way of doing this - lifting wave and lateral wave are essential the same thing, just that the normal direction of the wave is different. You can eventually merge the two forcing class into one. But of course, this is after resolving the current issue you have. Please let me know if these help.
AkashVardhan7 commented 11 months ago

Thanks for your reply, I changed the direction of actuation to the binormal_direction = np.array([1.0,0.0,0.0]) and lowered the amplitude of the lifting wave to 0.07, I think there is still a force unbalance as you can see from the slipping snake in the video snake_diffraction_test7

AkashVardhan7 commented 11 months ago

image These are the lifting torque parameters with which I am running my simulation... lateral_direction is the same binormal direction vector pointing along x axis with unit vector [1.0,0.0,0.0]. The video shows that the snake behaves chaotically after slipping, which is indicating that there is an unbalance in the net forces somewhere., while I am running it with a time step of 1e-4.

https://github.com/GazzolaLab/PyElastica/assets/48886152/24f33303-cd8b-468c-8e0e-9c61ec7cc7e8

xzhan139 commented 11 months ago

Hi, this looks to me like a numerical unstability due to your parameter settings. As I can see from your code there are few parameters that are quite different from what we usually run with, for example you have a period = 15, which makes the snake run very slow. We will push a case in PyElastica soon, which will include functions for lifting wave etc. The parameters of the snake will be slightly different, but you can use it as a reference for your implementation. Thank you.

AkashVardhan7 commented 11 months ago

Thanks please let me know when you have done so, and I can modify the code accordingly..When I changed the period to 2 second which is what is used in the paper SI and the time step to 8e-6 again the value used from the period, I still find that snake_diffraction_test9 the numerical issue persists,

xzhan139 commented 11 months ago

couple of other things to try -

  1. Chang the number of rod elements from 50 to 25. This will make each element heavier, will in general enhance numerical stability.
  2. Increasing your damping constant to 1e-1.
  3. Initialize your snake in the air, remove gravity, remove lateral wave, and plot the snake sideview video. You should see the snake floating and doing lifting wave. If this is not stable, then it is something in your lifting wave forcing itself.
AkashVardhan7 commented 11 months ago

The damping and number of rods helped and made the simulation more stable, I am observing sidewinding like behavior but the trajectory still seems a little funky in terms of the dynamics. I will try turning off the lateral wave and plotting just the side view with snake initialized in air.

https://github.com/GazzolaLab/PyElastica/assets/48886152/7529ec8a-2e46-42b9-b1d1-806a4f823b2d

AkashVardhan7 commented 11 months ago

I think it's the phase offset and the wavelength of the lifting wave which is giving the snake this waveform, I am using the same wavelength as the lateral wave and a phase shift of 0.25, which is chosen from the region in the parameter space supposed to generate straight side winding for lambda = 1.

AkashVardhan7 commented 11 months ago

Here's the lifting wave without the lateral wave and gravity. snake_diffraction_test11

AkashVardhan7 commented 10 months ago

This is the best side-winding I got snake_diffraction_test13 When are you going to push your version of the code with sidewinding on github?

Ali-7800 commented 10 months ago

Hi @AkashVardhan7,

Sorry for the delay I had to resolve some issues, I just opened the pull request with the new continuum snake with lifting wave example, it is in review now.

Ali-7800 commented 10 months ago

Hi @AkashVardhan7,

The continuum snake with lifting wave case has been add, you can check it out in the update-0.3.2 branch

armantekinalp commented 10 months ago

Since we have merged PR, I am closing this issue. If you have questions, please reopen again.