Workflow

When using PyElastica, users will setup a simulation in which they define a system of rods, define initial and boundary conditions on the rods, run the simulation, and then post-process the results. Here, we outline a typical template of using PyElastica.

Important

A note on notation: Like other FEA packages such as Abaqus, PyElastica does not enforce units. This means that you are required to make sure that all units for your input variables are consistent. When in doubt, SI units are always safe, however, if you have a very small length scale (\(\sim\) nm), then you may need to rescale your units to avoid needing prohibitively small time steps and/or roundoff errors.

1. Setup Simulation

from elastica.wrappers import (
    BaseSystemCollection,
    Connections,
    Constraints,
    Forcing,
    CallBacks
)

class SystemSimulator(
    BaseSystemCollection,
    Constraints, # Enabled to use boundary conditions 'OneEndFixedBC'
    Forcing,     # Enabled to use forcing 'GravityForces'
    Connections, # Enabled to use FixedJoint
    CallBacks    # Enabled to use callback
): 
    pass 

This simply combines all the wrappers previously imported together. If a wrapper is not needed for the simulation, it does not need to be added here.

Available components are:

Component

Note

BaseSystemCollection

Required for all simulator.

Constraints

Forcing

Connections

CallBacks

Note

We adopted a composition and mixin design paradigm in building elastica. The detail of the implementation is not important in using the package, but we left some references to read here.

2. Create Rods

Each rod has a number of physical parameters that need to be defined. These values then need to be assigned to the rod to create the object, and the rod needs to be added to the simulator.
from elastica.rod.cosserat_rod import CosseratRod

# Create rod
direction = np.array([0.0, 0.0, 1.0])
normal = np.array([0.0, 1.0, 0.0])
rod1 = CosseratRod.straight_rod(
    n_elements=50,                                # number of elements
    start=np.array([0.0, 0.0, 0.0]),              # Starting position of first node in rod
    direction=direction,                          # Direction the rod extends
    normal=normal,                                # normal vector of rod
    base_length=0.5,                              # original length of rod (m)
    base_radius=10e-2,                            # original radius of rod (m)
    density=1e3,                                  # density of rod (kg/m^3)
    nu=1e-3,                                      # Energy dissipation of rod
    youngs_modulus=1e7,                           # Elastic Modulus (Pa)
    poisson_ratio=0.5,                            # Poisson Ratio
)
rod2 = CosseratRod.straight_rod(
    n_elements=50,                                # number of elements
    start=np.array([0.0, 0.0, 0.5]),              # Starting position of first node in rod
    direction=direction,                          # Direction the rod extends
    normal=normal,                                # normal vector of rod
    base_length=0.5,                              # original length of rod (m)
    base_radius=10e-2,                            # original radius of rod (m)
    density=1e3,                                  # density of rod (kg/m^3)
    nu=1e-3,                                      # Energy dissipation of rod
    youngs_modulus=1e7,                           # Elastic Modulus (Pa)
    poisson_ratio=0.5,                            # Poisson Ratio
)

# Add rod to SystemSimulator
SystemSimulator.append(rod1)
SystemSimulator.append(rod2)

This can be repeated to create multiple rods. Supported geometries are listed in API documentation.

Note

The number of element (n_elements) and base_length determines the spatial discretization dx. More detail discussion is included here.

4. Define Boundary Conditions, Forcings, and Connections

Now that we have added all our rods to SystemSimulator, we need to apply the relevant boundary conditions. See this page for in-depth explanations and documentation.

As a simple example, to fix one end of a rod, we use the OneEndFixedBC boundary condition (which we imported in step 1 and apply it to the rod. Here we will be fixing the \(0^{\text{th}}\) node as well as the \(0^{\text{th}}\) element.

from elastica.boundary_conditions import OneEndFixedBC

SystemSimulator.constrain(rod1).using(
    OneEndFixedBC,                  # Displacement BC being applied
    constrained_position_idx=(0,),  # Node number to apply BC
    constrained_director_idx=(0,)   # Element number to apply BC
)

We have now fixed one end of the rod while leaving the other end free. We can also apply forces to free end using the EndpointForces. We can also add more complex forcings, such as friction, gravity, or torque throughout the rod. See this page for in-depth explanations and documentation.

from elastica.external_forces import EndpointForces

#Define 1x3 array of the applied forces
origin_force = np.array([0.0, 0.0, 0.0])
end_force = np.array([-15.0, 0.0, 0.0]) 
SystemSimulator.add_forcing_to(rod1).using(
    EndpointForces,                 # Traction BC being applied
    origin_force,                   # Force vector applied at first node
    end_force,                      # Force vector applied at last node
    ramp_up_time=final_time / 2.0   # Ramp up time 
)

One last condition we can define is the connections between rods. See this page for in-depth explanations and documentation.

from elastica.connections import FixedJoint

# Connect rod 1 and rod 2. '_connect_idx' specifies the node number that 
# the connection should be applied to. You are specifying the index of a 
# list so you can use -1 to access the last node. 
SystemSimulator.connect(
    first_rod  = rod1, 
    second_rod = rod2, 
    first_connect_idx  = -1, # Connect to the last node of the first rod. 
    second_connect_idx =  0  # Connect to first node of the second rod. 
    ).using(
        FixedJoint,  # Type of connection between rods
        k  = 1e5,    # Spring constant of force holding rods together (F = k*x)
        nu = 0,      # Energy dissipation of joint
        kt = 5e3     # Rotational stiffness of rod to avoid rods twisting
        )

5. Add Callback Functions (optional)

If you want to know what happens to the rod during the course of the simulation, you must collect data during the simulation. Here, we demonstrate how the callback function can be defined to export the data you need. There is a base class CallBackBaseClass that can help with this.

Note

PyElastica does not automatically saves the simulation result. If you do not define a callback function, you will only have the final state of the system at the end of the simulation.

from elastica.callback_functions import CallBackBaseClass

# MyCallBack class is derived from the base call back class.   
class MyCallBack(CallBackBaseClass):
    def __init__(self, step_skip: int, callback_params):
        CallBackBaseClass.__init__(self)
        self.every = step_skip
        self.callback_params = callback_params
    
    # This function is called every time step
    def make_callback(self, system, time, current_step: int):         
        if current_step % self.every == 0:
            # Save time, step number, position, orientation and velocity
            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["directors"].append(system.director_collection.copy())
            self.callback_params["velocity" ].append(system.velocity_collection.copy())
            return

# Create dictionary to hold data from callback function
callback_data_rod1, callback_data_rod2 = defaultdict(list), defaultdict(list)

# Add MyCallBack to SystemSimulator for each rod telling it how often to save data (step_skip)
SystemSimulator.collect_diagnostics(rod1).using(
    MyCallBack, step_skip=1000, callback_params=callback_data_rod1)
SystemSimulator.collect_diagnostics(rod2).using(
    MyCallBack, step_skip=1000, callback_params=callback_data_rod2)

You can define different callback functions for different rods and also have different data outputted at different time step intervals depending on your needs. See this page for more in-depth documentation.

6. Finalize Simulator

Now that we have finished defining our rods, the different boundary conditions and connections between them, and how often we want to save data, we have finished setting up the simulation. We now need to finalize the simulator by calling

SystemSimulator.finalize()

This goes through and collects all the rods and applied conditions, preparing the system for the simulation.

7. Set Timestepper

With our system now ready to be run, we need to define which time stepping algorithm to use. Currently, we suggest using the position Verlet algorithm. We also need to define how much time we want to simulate as well as either the time step (dt) or the number of total time steps we want to take. Once we have defined these things, we can run the simulation by calling integrate(), which will start the simulation.

We are still actively testing different integration and time-stepping techniques, PositionVerlet is the best default at this moment.

from elastica.timestepper.symplectic_steppers import PositionVerlet
from elastica.timestepper import integrate

timestepper = PositionVerlet()
final_time = 10   # seconds
dt = 1e-5         # seconds
total_steps = int(final_time / dt) 
integrate(timestepper, SystemSimulator, final_time, total_steps)

More documentation on timestepper and integrator is included here

8. Post Process

Once the simulation ends, it is time to analyze the data. If you defined a callback function, the data you outputted in available there (i.e. callback_data_rod1), otherwise you can access the final configuration of your system through your rod objects. For example, if you want the final position of one of your rods, you can get it from rod1.position_collection[:].