Particle Spring Simulation Tutorial

A tutorial for the course Digital World 2015 demonstrating a bare simple particle spring simulation for Processing/Python. To run the code below, download Processing and install the Python editing mode. The tutorial contains a series of suggestions/challenges for improving this basic template.

#--
#-- Particle Spring Dynamics Simulation
#--

#-- Window Dimensions
#--
dx = 800
dy = 600

#-- Definition of Complex Types ---------------------------------------------
#--
#-- As not sure if you are familiar with such concepts as a class or a tuple, 
#-- I use below vanilla lists to express semantically complex data forms such 
#-- as a vector, a particle and a spring.
#--
#-- The convention used is: Any complex type of information will be stored in 
#-- a list each element thereof representing a property in orderly manner.
#-- 
#-- The name of the property will be globally defined as a constant containing 
#-- the index of the property value in the list.
#--
#-- Example: A point/vector is a complex entity which contains a number of 
#-- coordinates. For 2D vectors we need two entries the first one being X 
#-- and second Y coordinate. Below is the definition of X = 0 and Y = 1
#-- Therefore I can get/set the X coordinate of point P using P[X]
#--

#-- Vector Definition -------------------------------------------------------
#--
#-- A vector represents a position, displacement or force in two dimensions.
#-- Property definitions (X,Y coordinates) and basic algebra for vectors are
#-- expressed below as global functions.
#--
#-- CHALLENGE: Modify the definion to support 3D points/vectors (easy) 
#--
X = 0  
Y = 1  
    
#-- Create a New Vector
#--
def vnew( x, y ):
    return [x, y]

#-- Add Two Vectors
#--
def vadd( u, v ):
    return vnew( u[X] + v[X], u[Y] + v[Y] )

#-- Subtract Two Vectors
#--
def vsub( u, v ):
    return vnew( u[X] - v[X], u[Y] - v[Y] )

#-- Vector Scale
#-- 
def vmul( u, s ):
    return vnew( u[X] * s, u[Y] * s )

#-- Vector Dot Product
#--
def vdot( u, v ):
    return u[X] * v[X] + u[Y] * v[Y]

#-- Vector Length
#--
def vlen( u ):
    return sqrt( vdot( u, u ) )

#-- Particle Definition -----------------------------------------------------
#--
#-- A particle represent a dimensionless mass with internal dynamics state of 
#-- position and velocity. The force acting upon the particle can be also 
#-- though of as its acceleration. In other words the particle captures a set
#-- of information about a vector valued function and its first and second
#-- derivatives. 
#--
MASS     = 0 
POSITION = 1
VELOCITY = 2
FORCE    = 3

#-- Create Particle
#--
def pnew( m, p, v, f ):
    return [m, p, v, f]

#-- Spring Definition -------------------------------------------------------
#-- A spring represents an idealized elastic link between two particles that
#-- obays Hooke's Law of springyness: F = k * d, where F is the force magnitude
#-- k is the spring constant, and d is the displacement.
#--
#-- Alternatively you can think of this as: F = k * ( RestLength - CurrentLength )
#-- where RestLength is the length at which the spring excerts zero force to both
#-- particles (aka ends of the spring locations) and CurrentLength is the length
#-- which is either bigger (spring is stretched) or smaller (spring is compressed)
#-- than the rest length (if they are equal then F = 0).
#--
#-- In terms of properties we store the source and target locations (the end-points)
#-- of the spring such that we can connect them with other springs or fix them etc.
#-- We obviously need the spring constant and rest length. 
#--
#-- NOTE: The source and target properties are expressed as indices (integers)
#--       instead of particle objects (list of properties). The conventions is 
#--       made such that we can store multiple particles in one list (each having
#--       a known index) and define springs that associate/point-into the particle 
#--       collection by index (instead of duplicating particles per spring). 
#--
#-- CHALLENGE #1: Can we change the representation from the high-school notion
#--               of the spring, to the proper engineering idea of elasticity? (easy)
#--               HINT: google stress/strain and replace k with the elastic modulus e
#--
#-- CHALLENGE #2: The spring, apart from its direct reaction to loading expressed
#--               by Hooke's Law, it also has a force component, not used above, which
#--               captures its internal damping. Otherwise a spring would oscillate
#--               indefinately if think about it. Can express this also (medium)
#--               HINT: google linear spring damping
#--
SOURCE    = 0  #-- Source Particle INDEX
TARGET    = 1  #-- Target Particle INDEX
STIFFNESS = 2  #-- Spring Stiffness Constant
LENGTH    = 3  #-- Rest Length

def snew( s, t, k, l ):
    return [s, t, k, l]

#-- Particle and Spring Lists -----------------------------------------------
#--
particles = []
springs   = []

#-- Create Initial Model ----------------------------------------------------
#--
#-- This function is required by processing in order to prepare the program.
#-- It will be the first thing to be executed before any drawing takes place.
#-- In here we will construct a chain of particles and connect them with 
#-- springs. They will be stored in the global variables "particles" and
#-- "springs" for later use in the simulation/animation cycle aka draw( )
#--
#-- NOTE: There is a conceptual leap here... it is only a convention or some
#--       sort of approximation involved in the assumption that a sequentially
#--       connected set of particles represents a "physical" chain. In the
#--       fashion we may "define" a rectangular grid of connected particles
#--       as a "fabric" sheet. More on that later, meanwhile think of other 
#--       "things" that can be represented by such configurations of particles
#--       as springs.   
#--
def setup( ):    
    
    #-- Window Size
    #--
    size( dx, dy )
    frameRate( 24 )
    
    #-- Import Globals
    #--
    global particles
    global springs

    #-- Number of Particles
    #--
    count = 12
    
    #-- Construct Particles ---------------------------------
    #-- Create a sequence of particles along the x-direction
    #-- NOTE: We need to provide mass, possition and velocity
    #--       The last two items are important as they define
    #--       the dynamics state of the particle. You can set
    #--       velocity to non-zero if you like.
    #-- 
    for index in range( 0, count ):
        mass = 20.0
        
        x = index * 50 + 100
        y = 100
        
        position = vnew( x, y )
        velocity = vnew( 0, 0 )
        force    = vnew( 0, 0 )
        
        particles.append( pnew( mass, position, velocity, force ) )

    #-- Construct Springs -----------------------------------
    #-- Connect previously constructed particles in order
    #--
    for index in range( 1, count ):
        source_index = index - 1
        target_index = index - 0
        
        stiffness   = 15.0
        rest_length = 50.0
        
        springs.append( snew( source_index, target_index, stiffness, rest_length ) )
    
    #-- By own convention, particles with zero mass are assumed fixed in place
    #-- There are other ways to define such a constraints but for simplicity 
    #-- we define here as such. The first and last particles of the chain are fixed.
    #--
    particles[        0][MASS] = 0.0
    particles[count - 1][MASS] = 0.0
    
    
#-- Draw Window Contents -----------------------------------------------------
#-- This is also a standard function required by processing. It is executed 
#-- every some odd milliseconds and it is supposed to draw the screen (aka FPS)
#-- 
#-- In here we perform two tasks: 1. We run the simulation of particles and
#-- springs. 2. We draw them in the window as we are supposed.
#--
#-- Logic of Simulation: the steps taken for simulation are...
#-- 1. Force Calculation/Accumulation     F = m * a  -> 
#--                                       a = F / m
#--
#-- 2. Velocity update from acceleration  a = dv / dt -> 
#--                                       dv = a * dt ->  
#--                                       v_new - v_old = a * 1.0 ->  (!!! dv = v_new - v_old )
#--                                       v_new = v_old + a
#-- 
#-- 3. Position update from velocity      v  = dp / dt ->
#--                                       dp =  v * dt ->
#--                                       p_new - p_old = v * 1.0 ->  (!!! dp = p_new - p_old )
#--                                       p_new = p_old + v
#--
#-- 4. Enforce Constraints (such that particle dont fall off the screen)
#--
#-- NOTE: The above logic is called forward explicit Newton integration.
#--       The places with (!!!) are important because there are assumptions
#--       of approximating differential equations with their first linear
#--       term only (which will cause noticable inaccuracy, unless dt is small)
#--
#--       Newton integration is simple and fast but tends to spiral out of 
#--       control due to numerical errors (we are just following the tangent
#--       line forward instead of the curved true trajectory, picture it).
#--
#--       The second assumption is that the time step dt = 1.0 which is a
#--       simplification that couples the animation FPS with simulation time.
#--       You may slow down the simulation by changing the FPS or actually
#--       multiplying velocity and acceleration by a small scalar eg. 0.0001
#--       Scaling the time step (dt) is advisable as it will improve accuracy. 
#--                                                
#-- CHALLENGE: Add a proper time step component. HINT: basic vector math (easy)
#--
#-- CHALLENGE: There other more sophisticated integration methods which are 
#--            also explicit but more accurate eg. Mid-Point (easy), Verlet (easy)
#--            RK4 (tricky), as well as implicit methods using linear algebra (hard)
#--            
#-- CHALLENGE: So we know from theory that the shape of handing chain is
#--            expressed-by / has-a closed form equation known as the catenoid.
#--            Draw it over the simulated shape (same parameters) and compare the 
#--            results. What is the significance of the error? 
#-- 
#-- CHALLENGE: So we have been using abstract units for space (pixels), mass(units)
#--            time (fps). What if you apply real units eg. meters, kilograms, seconds
#--            such that you can compare the real shape of a physical chain with the
#--            simulation. How does the real shape (take photo and overlay), theory
#--            and simulated form compare?
#-- 
def draw( ):
    
    #-- Import Globals
    #--
    global particles
    global springs




    #-- Reset Force Accumulator ---------------------------------
    #--
    for particle in particles:
        particle[FORCE] = vnew( 0.0, 0.0 )
       
       
       
       
       
    #-- Accumulate Gravity Forces -------------------------------
    #-- Assuming earth gravity constant but note y is positive
    #-- because y-screen direction is pointing downwards already.
    #--
    GRAVITY = 9.81

    for particle in particles:
        gravity = vnew( 0, GRAVITY * particle[MASS] )        
        particle[FORCE] = vadd( particle[FORCE], gravity )
          
          
          
          
          
          
    #-- Accumulate Damping Forces -------------------------------
    #-- Fictional linear drag opposite to particle's velocity
    #-- for stability. Air drag is quadratic to velocity and
    #-- proportional to fluid density as well as surface area.
    #--
    #-- CHALLENGE: Change this unrealistic drag force to a proper
    #--            air drag force (easy). HINT: google drag force
    #--            when is doubt about areas & constants use 1.0 
    #--
    DAMPING = 2.00
    for particle in particles:
        damping = vmul( particle[VELOCITY], -DAMPING )        
        particle[FORCE] = vadd( particle[FORCE], damping )







    #-- Accumulate Spring Forces -------------------------------------------
    #--
    for spring in springs:
        
        #-- Extract Particle Indices
        #--
        source_index = spring[SOURCE]
        target_index = spring[TARGET]
        
        #-- Extract Particles (themselves)
        #--
        source_particle = particles[source_index]
        target_particle = particles[target_index]
        
        #-- Extract Positions
        #--
        source_position = source_particle[POSITION]
        target_position = target_particle[POSITION]
        
        #-- Hook's Law (over-simplified)
        #-- The proper formulation includes a damping
        #-- component proportional to velocity along the
        #-- spring direction. NOTE: I am dividing by the
        #-- distance between the spring's end-points aka
        #-- the current length, because the force magnitude
        #-- is assuming a normalized directional vector        
        #-- 
        stiffness   = spring[STIFFNESS]
        rest_length = spring[LENGTH]
        direction   = vsub( target_position, source_position )
        distance    = vlen( direction ) 
        magnitude   = stiffness * ( distance - rest_length ) / distance
        
        source_particle[FORCE] = vadd( source_particle[FORCE], vmul( direction,  magnitude ) )
        target_particle[FORCE] = vadd( target_particle[FORCE], vmul( direction, -magnitude ) )







    #-- Newton Integration -----------------------------------------
    #-- Discrete integration by assuming time step of one unit.
    #-- For better accuracy and stability either scale acceleration
    #-- and velocity vectors by a time step eg. dt = 0.001 or use
    #-- better integration approximation.
    #-- 
    for particle in particles:
        
        #-- Zero mass particles are ignored
        #--
        mass = particle[MASS]
        
        if( mass > 0.0 ):
            #-- F = ma <=> a = F / m
            #--
            acceleration = vmul( particle[FORCE], 1.0 / mass )
            
            #-- acceleration = dv/dt ~> v_new = v_old + acceleration * dt
            #--
            particle[VELOCITY] = vadd( particle[VELOCITY], acceleration )
            
            #-- velocity = dp/dt ~> p_new = p_old + velocity * dt
            #--
            particle[POSITION] = vadd( particle[POSITION], particle[VELOCITY] );









    #-- Constraint Enforcement -----------------------------------
    #--
    #-- Force particles to remain within window boundaries by
    #-- reversing their position to previous and flipping velocity
    #--
    #-- This is kind of a perfect elestic collision responce since
    #-- velocity is preserved 100%
    #--
    #-- The reason for reversing position to previous is because
    #-- the particle may have escaped significantly outside. So
    #-- better safe than sorry.
    #--
    #-- Constraints are better expressed as "force constraints"     
    #-- rather than "penalties". If you think about it, the code
    #-- below is not 100% realistic because every time we move
    #-- particles back and flip their velocity direction we are
    #-- fiddling with the energy/momentum of the system in ways 
    #-- that may violate constitutional laws of coservation.
    #--
    #-- CHALLENGE: Read about the notion of "virtual work" and
    #-- implement a constraint such as particle-on-circle (hard)
    #--
    for particle in particles:
        
        #-- Particle Properties
        #--
        position = particle[POSITION]
        velocity = particle[VELOCITY]
        radius   = particle[MASS] * 0.5  #-- Drawing artifact, see below
        
        #-- Bound Checking
        #-- 
        if( ( position[X] < radius ) or ( position[X] > ( dx - radius ) ) ):
            position[X] = position[X] - velocity[X]  
            velocity[X] =-velocity[X]
            
        if( ( position[Y] < radius ) or ( position[Y] > ( dy - radius ) ) ):
            position[Y] = position[Y] - velocity[Y] 
            velocity[Y] =-velocity[Y]
        
        
        
        
        
        
        
        
        
    #-- Prepare Frame -------------------------------------------
    #-- Clears the background, otherwise it will look psychedelic
    #-- try commenting out the background line below and enjoy
    #--
    fill( 0 )
    background( 255 )        
        
        
        
        
    #-- Draw Springs --------------------------------------------
    #--
    for spring in springs:
        
        source_index = spring[SOURCE]
        target_index = spring[TARGET]
        
        source_particle = particles[source_index]
        target_particle = particles[target_index]
        
        source_position = source_particle[POSITION]
        target_position = target_particle[POSITION]
                
        stroke( 0, 0, 255 )
        line( source_position[X], source_position[Y],
              target_position[X], target_position[Y] );
        
        
        
        
    #-- Draw Particles ---------------------------------------------
    #--
    noStroke( )
    index = 0
    for particle in particles:
        
        position = particle[POSITION]
        mass     = particle[MASS]

        #-- Fixed particles (with zero mass) have green color 
        #-- RGB=(0,255,0), unconstrained particles have red 
        #-- color (0,255,0) and radius = half their mass
        #--        
        if( mass > 0.0 ):
            radius = 0.5 * mass
            fill( 255, 0, 0 )
        else:
            fill( 0, 255, 0 )
            radius = 10
        
        ellipse( position[X], position[Y], radius, radius )
        
        #-- Display Particle Index
        #-- 
        fill( 0 )
        text( str( index ), position[X] + 10, position[Y] - 10 )
        index = index + 1