Timoschenko Beam example-length of Beam remains the same

[matei1996]

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

[armantekinalp]

Hello @matei1996

Yes this behavior is expected. Cosserat rods can also extend (stretch) in addition to bend.

[matei1996]

Ok, thanks. I am trying to simulate a wire using the Cosserat rods. However for this I would need that the rod doesn't stretch when I apply a force onto it. Do you know how I can achieve this? I am trying to reconstruct this image.

However I can only achieve this using the Cosserat Rod.

I would need the rod to have a more circular figure.

[armantekinalp]

For that you can modify stretch stiffness (EA) of the rod. After creating rod just do

rod.shear_matrix[2,2,:] *= 1000

This would increase the stretch stiffness of the rod by 1000 times. You can play around with this parameter.

[matei1996]

That doesn't seem to work I added my code below class TimoshenkoBeamSimulator(BaseSystemCollection, Constraints, Forcing): pass

timoshenko_sim = TimoshenkoBeamSimulator()

# setting up test params
n_elem = 100 #number of elements in which rod is divided

#Material properties
density = 1000
nu = 0.1
E = 1e6
# For shear modulus of 1e4, nu is 99!
poisson_ratio = 99
shear_modulus = E / (poisson_ratio + 1.0)

#Geometry of the rod; location/orientation of rod and length and radius
start = np.zeros((3,))
direction = np.array([0.0, 0.0, 1.0])
normal = np.array([0.0, 1.0, 0.0])
base_length = 3.0
base_radius = 0.25
base_area = np.pi * base_radius ** 2

shearable_rod = CosseratRod.straight_rod(
n_elem, 
start,
direction,
normal,
base_length,
base_radius,
density,
nu,
E,
shear_modulus=shear_modulus,
)
shearable_rod.shear_matrix[2,2,:] *= 10

timoshenko_sim.append(shearable_rod) # needs to be included otherwise it is not added to the simulation

#Boundary conditions

timoshenko_sim.constrain(shearable_rod).using(
    OneEndFixedRod, constrained_position_idx=(0,), constrained_director_idx=(0,)
)
print("One end of the rod is now fixed in place")

#Apply force to the end
origin_force = np.array([0.0, 0.0, 0.0])
end_force = np.array([-15, 0.0, 0.0])
ramp_up_time = 5.0

timoshenko_sim.add_forcing_to(shearable_rod).using(
    EndpointForces, origin_force, end_force, ramp_up_time=ramp_up_time
)
print("Forces added to the rod")

#Finalize system
timoshenko_sim.finalize()
print("System finalized")

#Define run time of the system
final_time = 20.0
dl = base_length / n_elem
dt = 0.01 * dl
total_steps = int(final_time / dt)
print("Total steps to take", total_steps)

timestepper = PositionVerlet()

#Calculate solution
integrate(timestepper, timoshenko_sim, final_time, total_steps) 

The image that I am getting is

The blue rod is the one generated with the code above. Ideally I would expect that the right end point is a little shifted to the left, the base length of my rod is 3.0, therefore the x-value of my right hand point after applying some force onto it, shouldn't also be 3.0 but smaller than that.

[matei1996]

Hello again,

I just went over the documentation again and noticed that the force isn't applied from above on the rod, but rather the rod is pulled into that direction. Now I understand why the rod doesn't change its x-value. However, I am interested into applying a force from above onto the beam. Is that possible in this model?

[armantekinalp]

Forces are always applied on the centerline of the rod (on the nodes). I think what you are trying to do apply a force always perpendicular to the tangent vector of the rod at any point. Is that correct? Then what you can do is using the local frame of the rod (directors) you can define a force vector which is always perpendicular to the rod tangent.

I suggest defining force vector in the direction of rod normal or binormal direction, which is d1 or d2 . You can slice the director_collection array to get the corresponding values for d1 and d2.

[matei1996]

Hello, sorry for the late reply, I was on holidays. Exactly I need a force vector that is perpendicular to the tangent vector of the last node (i.e. where the force is applied). However, I do not understand yet how to do it. After I generate the rod I can use the director_collection array to get the values of d1 and d2. And then I apply the force in the corresponding direction? But do I not have to do that at every time step? So I need a loop to do that?

[matei1996]

Hello, is it possible to give starting positions for all the nodes? So for example I apply a force onto my rod. Then I want to use "position_collection" to determine the position of all the nodes. And then start a new experiment with those values.

A more simple explanation: I want to apply a force on to a rod of -3 for example. And after it finished I want to another force of a different magnitude on the rod which is in the last position of the first experiment.

Or rather are there objects of the type "non_straight_rod"

[armantekinalp]

After you are done, you can save position_collection and director_collection of the rod into an npz file. In the next experiment you can load them and pass position and directors to the rod at the construction step. straight_rod method takes positions and directors as keyword arguments.

[matei1996]

So I am generating a rod via a function and then I save the results using savez_compressed

shearable_rod = GenTimoShear_base()
x=shearable_rod.position_collection
y=shearable_rod.director_collection

savez_compressed('data.npz',x,y)
data_path='data.npz'

Afterwards I load the data into my function that generates the second rod with all the other parameters

#Geometry of the rod; location/orientation of rod and length and radius
    start = np.zeros((3,))
    direction = np.array([0.0, 0.0, 1.0])
    load_data = load(data_path)
    position = load_data['arr_0']
    directors = load_data['arr_1']
    normal = np.array([0.0, 1.0, 0.0])
    base_length = 3.0
    base_radius = 0.25
    base_area = np.pi * base_radius ** 2

and then construct the rod

    shearable_rod = CosseratRod.straight_rod(
    n_elem, 
    start,
    direction,
    normal,
    base_length,
    base_radius,
    density,
    nu,
    E,
    shear_modulus=shear_modulus,
    position=position,
    directors=directors,
    )

However I get the following error message:

Do you have an idea what might cause this error? Are the values of 'start' and 'direction' supposed to take some specific values based on the ones I give in 'position' and 'directors'?

[armantekinalp]

I think your directors and tangents are not aligned because there is a shear on your rod. One other way is to set rod directors just after the construction. You can do shearable_rod.director_collection[:] = directors[:].

[matei1996]

Thanks that worked! Now I still have the problem, that I can't apply the force perpendicular to the tangent vector of the last node. I use:

new_force = unshearable_rod.director_collection[:,1,-1]

timoshenko_sim.add_forcing_to(unshearable_rod).using(
    EndpointForces, origin_force, -3*new_force, ramp_up_time=ramp_up_time
)

And tried different scalings, as well as using the d_1 vector (that would be unshearable_rod.director_collection[:,0,-1]) right?)

I also tried to use arbitrary parameters for the end force. However, the x-coordinate isn't changing no matter what I try. (see picture below). Is "new_force" the correct vector?

nging

[armantekinalp]

In PyElastica d1,d2,d3 are in row order instead of column. For example to access the last elements d1 (normal vector) shearable_rod.director_collection[0,:,-1] and for d3 (in case no shear same as tangent) shearable_rod.director_collection[2,:,-1]

[matei1996]

Ok, I applied it and there are still some things not working out. I started however to implement a function that computes the rod position for a short amount of time. Then saves the position and directions and afterwards starts a new experiment with those values and applying a force that is perpendicular to the tangent vector of the last node. Below is the code

def GenTimoShearDyn_new(x):

    # setting up test params
    n_elem = 100 #number of elements in which rod is divided

    #Material properties
    density = 1000
    nu = 0.1
    E = 1e6
    # For shear modulus of 1e4, nu is 99!
    poisson_ratio = 99
    shear_modulus = E / (poisson_ratio + 1.0)

    #Geometry of the rod; location/orientation of rod and length and radius
    start = np.zeros((3,))
    #start = np.array([0.0, 0.0, 0.0])
    direction = np.array([0.0, 0.0, 1.0])
    normal = np.array([0.0, 1.0, 0.0])
    base_length = 9.2
    base_radius = 0.25
    base_area = np.pi * base_radius ** 2

    #position = rod.position_collection
    #directors = rod.director_collection

    unshearable_start = np.array([0.0, -1.0, 0.0])
    #Define run time of system
    #final_time = 10.0
    dl = base_length / n_elem
    dt = 0.01 * dl
    final_time= dt
    total_steps = int(final_time / dt)
    print("Total steps to take", total_steps)

    time = 0.0

    class BeamSimulator(BaseSystemCollection, Constraints, Forcing):
        pass

    dynamic_update_sim = BeamSimulator()

    unshearable_rod_new = CosseratRod.straight_rod(
        n_elem,
        unshearable_start,
        direction,
        normal,
        base_length,
        base_radius,
        density,
        nu,
        E,
        poisson_ratio=-0.85,
        #position=position,
    )
    #unshearable_rod.director_collection[:] = directors[:]

    dynamic_update_sim.append(unshearable_rod_new)
    dynamic_update_sim.constrain(unshearable_rod_new).using(
        OneEndFixedRod, constrained_position_idx=(0,), constrained_director_idx=(0,)
    )
        #Apply force to the end
    origin_force = np.array([0.0, 0.0, 0.0])
    end_force = np.array([x[0], 0.0, 0.0]) #-10.0
    ramp_up_time = 0.1

    dynamic_update_sim.add_forcing_to(unshearable_rod_new).using(
    EndpointForces, origin_force, end_force, ramp_up_time=ramp_up_time
    )    

    dynamic_update_sim.finalize()

    #Define run time of system
    #evolve_for_time = 10.0
    evolve_for_time = dt
    update_interval = 1.0e-1

    ft=5
    final_time = 1
    dl = base_length / n_elem
    dt = 0.01 * dl
    total_steps = int(final_time / dt)
    t=0
    print("Total steps to take", total_steps)

    timestepper = PositionVerlet()

    #Calculate solution
    integrate(timestepper, dynamic_update_sim, final_time, total_steps)

    #Save computed values
    x=unshearable_rod_new.position_collection
    y=unshearable_rod_new.director_collection
    savez_compressed('data.npz',x,y)
    data_path='data.npz'

    def plot_timoshenko_dynamic_test(unshearable_rod):
        import matplotlib.pyplot as plt
        from IPython import display

        #ax.clear()
        fig = plt.figure(figsize=(5, 4), frameon=True, dpi=150)
        ax = fig.add_subplot(111)
        #ax.grid(b=True, which="major", color="grey", linestyle="-", linewidth=0.25)

        ax.plot(
            unshearable_rod.position_collection[2, :],
            unshearable_rod.position_collection[0, :],
            "k-",
            #label="unshearable rod",
        )

        #ax.legend(prop={"size": 12}, loc="lower left")
        #ax.set_ylabel("Y Position (m)", fontsize=12)
        #ax.set_xlabel("X Position (m)", fontsize=12)
        #ax.set_title("Simulation Time: %0.2f seconds" % time)
        ax.set_xlim([0.0, 9.2])
        ax.set_ylim([-0.01, 0.045])#0.045
        #ax.set_ylim([x[1],0.045])
        plt.axis('off')

    %matplotlib inline
    import matplotlib.pyplot as plt
    from IPython import display

    plot_timoshenko_dynamic_test(unshearable_rod_new)

    while t<ft:

        #Load old values
        load_data = load(data_path)
        position = load_data['arr_0']
        directors = load_data['arr_1']

        dynamic_update_sim = BeamSimulator()

        start = position[:,0]
        direction = directors[0,:,0]

        unshearable_rod_new = CosseratRod.straight_rod(
            n_elem,
            start,
            direction,
            normal,
            base_length,
            base_radius,
            density,
            nu,
            E,
            poisson_ratio=-0.85,
            position=position,
            #directors=directors,
        )
        unshearable_rod_new.director_collection[:] = directors[:]

        #Apply force to the end
        origin_force = np.array([0.0, 0.0, 0.0])
        #end_force = np.array([x[0], 0.0, 0.0]) #-10.0
        ramp_up_time = 0.1

        dynamic_update_sim.append(unshearable_rod_new)
        dynamic_update_sim.constrain(unshearable_rod_new).using(
            OneEndFixedRod, constrained_position_idx=(0,), constrained_director_idx=(0,)
        )

        #Apply force in perpendicular direction of the tangent of the last node
        force=unshearable_rod.director_collection[0,:,-1]

        dynamic_update_sim.add_forcing_to(unshearable_rod_new).using(
        EndpointForces, origin_force, force, ramp_up_time=ramp_up_time
        )

        dynamic_update_sim.finalize()
            #Define run time of system

        final_time = 1
        total_steps = int(final_time / dt)
        print("Total steps to take", total_steps)

        timestepper = PositionVerlet()

        #Calculate solution
        integrate(timestepper, dynamic_update_sim, final_time, total_steps)
        t=t+dt

        x=unshearable_rod_new.position_collection
        y=unshearable_rod_new.director_collection
        savez_compressed('data.npz',x,y)
        data_path='data.npz'

        plot_timoshenko_dynamic_test(unshearable_rod_new)

    #segFrames()
    return dynamic_update_sim

However, after the first iteration I get again the error that the rod normal and tangent are not perpendicular to each other. Can you help me find the issue?

[armantekinalp]

Hi @matei1996

Maybe use restart functionality in PyElastica? It might be a better option. Restart function will load the all rod parameters from saved file. Will this be useful?

Here is the example case

https://github.com/GazzolaLab/PyElastica/blob/master/examples/RestartExample/restart_example.py#L81

[matei1996]

Hi, I think the restart function is what I am looking for. However, there are still some uncertainties. Below my Code

    def plot_timoshenko_dynamic_test(shearable_rod):
        import matplotlib.pyplot as plt
        from IPython import display

        #ax.clear()
        fig = plt.figure(figsize=(5, 4), frameon=True, dpi=150)
        ax = fig.add_subplot(111)
        #ax.grid(b=True, which="major", color="grey", linestyle="-", linewidth=0.25)

        ax.plot(
            shearable_rod.position_collection[2, :],
            shearable_rod.position_collection[0, :],
            "k-",
            #label="unshearable rod",
        )

        #ax.legend(prop={"size": 12}, loc="lower left")
        #ax.set_ylabel("Y Position (m)", fontsize=12)
        #ax.set_xlabel("X Position (m)", fontsize=12)
        #ax.set_title("Simulation Time: %0.2f seconds" % time)
        ax.set_xlim([0.0, 9.2])
        ax.set_ylim([-0.01, 0.045])#0.045
        #ax.set_ylim([x[1],0.045])
        plt.axis('off')

    %matplotlib inline
    import matplotlib.pyplot as plt
    from IPython import display
    """

    def run_and_update_plot(simulator, dt, start_time, stop_time, ax):
        from elastica.timestepper import extend_stepper_interface
        from elastica.timestepper.symplectic_steppers import PositionVerlet

        timestepper = PositionVerlet()
        do_step, stages_and_updates = extend_stepper_interface(timestepper, simulator)

        n_steps = int((stop_time - start_time) / dt)
        time = start_time
        for i in range(n_steps):
            time = do_step(timestepper, stages_and_updates, simulator, time, dt)
        plot_timoshenko_dynamic(shearable_rod_new, unshearable_rod_new, end_force, time, ax)
        return time

    def plot_timoshenko_dynamic(shearable_rod, time, ax):
        import matplotlib.pyplot as plt
        from IPython import display

        ax.clear()
        #ax.grid(b=True, which="major", color="grey", linestyle="-", linewidth=0.25)

        ax.plot(
            shearable_rod.position_collection[2, :],
            shearable_rod.position_collection[0, :],
            "k-",
            #label="shearable rod",
        )

        #ax.legend(prop={"size": 12}, loc="lower left")
        #ax.set_ylabel("Y Position (m)", fontsize=12)
        #ax.set_xlabel("X Position (m)", fontsize=12)
        #ax.set_title("Simulation Time: %0.2f seconds" % time)
        ax.set_xlim([0.0, 9.2])
        ax.set_ylim([-0.01, 0.045])#0.045
        #ax.set_ylim([x[1],0.045])
        plt.axis('off')
    """

    class RestartExampleSimulator(BaseSystemCollection, Constraints, Forcing, CallBacks):
        pass

    restart_example_simulator = RestartExampleSimulator()

    # setting up test params
    final_time = 20
    n_elem = 50
    start = np.zeros((3,))
    direction = np.array([0.0, 0.0, 1.0])
    normal = np.array([0.0, 1.0, 0.0])
    base_length = 3.0
    base_radius = 0.25
    base_area = np.pi * base_radius ** 2
    density = 300
    E = 1e6
    # For shear modulus of 1e4, nu is 99!
    poisson_ratio = 99
    shear_modulus = E / (poisson_ratio + 1.0)

    # Create rod
    shearable_rod = CosseratRod.straight_rod(
        n_elem,
        start,
        direction,
        normal,
        base_length,
        base_radius,
        density,
        0.0,  # internal damping constant, deprecated in v0.3.0
        E,
        shear_modulus=shear_modulus,
    )

    # Append rod to simulator
    # In this example we have one Cosserat rod, but you can add multiple Cosserat rods or rigid bodies to the simulator.
    restart_example_simulator.append(shearable_rod)

    # Constrain one end of the rod.
    restart_example_simulator.constrain(shearable_rod).using(
        OneEndFixedBC, constrained_position_idx=(0,), constrained_director_idx=(0,)
    )

    # Add end point forces to rod
    #end_force = np.array([1.0, 0.0, 0.0])
    end_force = shearable_rod.director_collection[1,:,-1]
    restart_example_simulator.add_forcing_to(shearable_rod).using(
        EndpointForces, 0.0 * end_force, -5*end_force, ramp_up_time=final_time / 10.0
    )

    # add damping
    #damping_constant = 10.0
    dl = base_length / n_elem
    dt = 0.01 * dl
    """
    restart_example_simulator.dampen(shearable_rod).using(
        AnalyticalLinearDamper,
        damping_constant=damping_constant,
        time_step=dt,
    )
    """
    # Finalize simulation
    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
    SAVE_DATA_RESTART = True
    restart_file_location = "data/"

    if LOAD_FROM_RESTART:
        restart_time = load_state(restart_example_simulator, restart_file_location, True)
    else:
        restart_time = np.float64(0.0)

    timestepper = PositionVerlet()
    total_steps = int(final_time / dt)

    time = integrate(
        timestepper,
        restart_example_simulator,
        final_time,
        total_steps,
        restart_time=restart_time,
    )

    #plot results dynamically

    plot_timoshenko_dynamic_test(shearable_rod)

    # 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)
    i=0
    while i<10:

        # 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 = True
        SAVE_DATA_RESTART = True
        restart_file_location = "data/"

        if LOAD_FROM_RESTART:
            restart_time = load_state(restart_example_simulator, restart_file_location, True)
        else:
            restart_time = np.float64(0.0)

        timestepper = PositionVerlet()
        total_steps = int(final_time / dt)

        time = integrate(
            timestepper,
            restart_example_simulator,
            final_time,
            total_steps,
            restart_time=restart_time,
        )
        print(shearable_rod.position_collection[:,-1])
        plot_timoshenko_dynamic_test(shearable_rod)

        # 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)
        #return dynamic_update_sim
        i=i+1

First, I used as my force vector, which should be perpendicular to the tangent at the end node the vector, shearable_rod.director_collection[1,:,-1]. As I understand this changes the y and x coordinates, below I have some outputs of the end_position of my rods. shearable_rod.director_collection[0,:,-1] changes z and x coordinates (The coordinate axes in the np array correspond [y,z,x] right?). However, the visual effect on my x-coordinate is very minimal, whereas the effect on the y-coordinate is very big. Contrary to the that the effect on the absolute value is small for both entries.

and visually

Might there be something wrong with my plot function?

[bhosale2]

Hi @matei1996 after briefly going through the conversation above, it seems like your issue is centered around implementing a derived functionality for your case study, and not regarding any issue in pyelastica API or functionality. Also given that there has been some difficulty in progressing towards the solution of your problem with the above conversations, I would recommend we tackle this via another route, such as a conversation or a call with one of our junior developers, who is well versed with pyelastica. Thoughts @armantekinalp?

[armantekinalp]

@bhosale2 it is a good idea. @Ali-7800 will take over the issue and help @matei1996 .

Thanks

[bhosale2]

Great. @matei1996 please reach out to @Ali-7800, after which we will close this issue.

[matei1996]

Hi, sounds like a good idea. @Ali-7800 can you send me an email and then we could meet up online? My Email is matei.hanu(at)fu-berlin.de

[bhosale2]

Hi @matei1996, @Ali-7800 has sent you a mail, and both of you can continue this issue via private communication. As such, Im closing this issue.