Particle Simulation

The objective of this tutorial is to develop a simple physics engine which can simulate dimensionless bodies. The physics knowledge required for this is actually high-school level and the aim is to practice on python programming. The computational concepts emphasized here pertain symbolic modeling of entities, such as particles and forces, which are already familiar. Advanced topics regarding simulation will be hinted about but not covered extensively.

A fantastic reference for an in depth introduction to the subject of particle and solid body dynamics, as well as advanced topics such as constraints and collision detection, can be found in the notes by David Baraff and Andrew Witkin from the 1999 Siggraph [📚]. Additional references are located at the end of the post. Here we will take a step-by-step introduction of dynamics for design-minded audiences.

Time Keeping & Static Variables

We will start with creating a clock component which will become the entry point of the simulation. The setup requires a python block with a (dt : float) and (reset : bool) input parameters. In addition, we need a timer block linked to the python block as seen below.

Block Diagram

The simulation needs to preserve state across executions of the script ie. the value of the time elapsed variable must be retained. We need a construct to express the idea of “if a variable is not defined then create it, otherwise update it”. But if we just assign the time variable in every clock execution, we would be stuck at a fixed value. To achieve the desired behavior, we need to store the variable’s contents in global scope such that in every execution we will increment its value instead of resetting it. This concept is called declaring a static variable. In python this concept is not supported directly, but we can emulate it as follows:

Python Definition

if( ( "time" in vars( ) ) and ( not reset ) ):
    time += dt
    print( time )
else:
    time = 0.0
    print( time )

Once a global variable is defined, by assignment in python, ie. outside a function at the top-most scope, it is registered within the built-in variables table maintained by the language’s run-time. This table can be obtained by using the vars( ) function. We can then check if a variable is already declared by using the “variable name as string, contained in variables table” construct. If the time variable is defined __and__ we are not in reset mode, then we can increment its value, otherwise we set it to zero.

Particle Class

The simulation operates on physical bodies with no dimensions, they are called particles. They can be thought of as points with mass and speed. A particle, at bare minimum, needs to know its position and velocity. Those two quantities define the kinematic state of the object. We can, and will do later, include additional attributes for each particle but for now those are sufficient.

Python Definition

class Particle:
    def __init__( self, position, velocity, mass ):
        self.Position = position
        self.Velocity = velocity
        self.SumForce = Vector3d.Zero
        self.Mass     = mass

The sum of forces property is used as temporary variable across simulation executions. It represents the cumulative force vector applied at an instance of time on the particle.

Maths and Physics Bits

The physics used, are based on Netwon’s laws of motion and calculus, in a simplified form of difference functions and summations often taught before proper calculus is introduced. What we want to compute is the positions at each time step ie. the geometry of the particle system. Positions change over time because forces cause the particles to move. Assuming no initial velocity, no forces means no motion. What we know at time zero, are the initial positions and the forces applied to particles, because have drawn and applied loads on them. Generally, we express an initial setup via computation or by converting an annotated 3d geometric model into a particle system. What we need is set of rules to transform this information into updated positions, so we can animate the system. We thus freeze time and observe the following:

Forces

We start from Newton’s law (F = m a) noting that (F) stands for the sum of forces applied to a body, (m) is its mass and (a) is its acceleration. So in each time instance, we aggregate forces to obtain the acceleration of particles by simply dividing by their mass (a = F / m).

Velocities

We know that acceleration is the change of velocity over time (a = dv / dt) or in discrete terms (a = Δv / Δt), where delta implies a sizable amount of time instead of an infinitesimal. Expanding we get (Δv = vnew – vold) and with a bit of rearrangement (vnew = vold + a Δt). We thus now have an expression for updating a particle’s velocity (v += a * dt). In fact it is quite remarkable that all calculus became mere subtractions and additions.

Positions

Similarly, velocity is the change of position over time (v = dp / dt) and in the same spirit (v = Δp / Δt). We finally arrive at (pnew = pold + v Δt). Which in computational terms becomes (p += v * dt). We have thus threaded a path between original positions and forces to position updates. We first compute the sum of forces, then we use the acceleration to update the particle’s velocity and we use the velocity to update the particle’s position.

The Simulation Class

The simulation class is a wrapper of all information we need to retain between script executions and the coordinator of the dynamic system. We will begin with a minimal definition and expand it as we introduce new concepts.

Python Definition

class Simulation:
    def __init__( self ):
        self.Particles   = []
        self.Forces      = []

    def Update( self, dt ):        
        for particle in self.Particles:
            particle.SumForce = Vector3d.Zero

        for force in self.Forces:
            force.Apply( )

        for particle in self.Particles:
            if( particle.Mass == 0 ): continue

            acceleration = particle.SumForce * ( 1.0 / particle.Mass )
            particle.Velocity += acceleration * dt
            particle.Position += particle.Velocity * dt

The simulation contains a list of particles and a list of forces (more on forces coming up). The key function of the simulation is to update particles’ states in discrete time steps. The process requires: (a) clearing all sums of forces for all particles, (b) accumulating all forces for all particles and (c) updating particle’s states. Clear and accumulate (a, b) is what we call force analysis and state update (c) is known as integration. The code above is literally a transcription of the math and physics into python.

Stationary Particles

In the code above, we use a trap (if particle.Mass == 0: continue) for particles with zero mass to avoid division by zero when computing the acceleration. This could have been enforced at a particle’s construction by throwing an error. However, we use this approach because it conveniently allows us to define particles fixed in space ie. if a particle’s update is skipped, it will not move from its original position.

Theory Bits

You may notice that the sequence of updating the particle’s state is a bit odd: We are changing the position using the velocity of the time step ahead. This is logically bogus, but there is tons of theory behind supporting this approach. What we are doing in the update loop is known as integration, here expressed as summation since we are working with discrete steps. The particular integrator scheme used here is call the semi-implicit Euler integrator or symplectic Euler integrator [📚] exactly because it mixes past and future.

The issue here is accuracy of simulation. Because we are using discrete time steps, we are notionally following a tangent line of a linearlized version of a system of differential equations. This mean we are not faithful against, most often, highly non-linear reality. The size of the time step controls the accuracy, ie. how far off reality we are. This can lead not only to wrong results and also what is known as instability ie. the system explodes [📚].

To keep the discussion short, symplectic integrators [📚] tend to be stable and fast, hence often used in video games and even scientific computing. If better accuracy and stability is required there are more involved schemes, using higher order of approximations, such as the Verlet [📚], Leapfrog [📚], Mid-point [📚], Forest-Ruth [📚] and Runge-Kutta [📚] methods. For additional information, look at the references section at the end of the post.

Desires and Forces

There a many types of forces which we can classify in terms of how many particles they involve: Unary forces are applied to particles individually without accounting others eg. gravity, point loads, viscous drag. Binary forces involve two particles which interact with one another as a pair eg. a linear spring can be considered as a force instead of a physical entity. N-ary forces involve interactions between multiple particles eg. planetary and electromagnetic forces.

We can classify forces in terms of their relationship to reality as natural and synthetic. It may come as a surprise but as long as we find a way to express how one or more particles interact with one another, we can formulate this in terms of position, velocity, mass, and any other attributes required, we can incorporate this into the simulation. Examples of natural forces include gravity and the spring forces while synthetic forces include swarm intelligence and computer graphics flames and explosions.

The key point is: forces express desires which all want to be achieved simultaneously. This may or may not be possible in every instance. It depends on the degrees of freedom of the system. But no matter the case, physics through the notion of energy, allows us to find an equilibrium, which often means a good compromise among conflicting causes. This is actually a remarkable way of thinking about computational problem solving.

Force of Gravity

Gravity is perhaps the simplest type of force we can implement. It follows Newton’s law of intelligent falling (F = mg) where (g) is the acceleration of gravity. The constructor requires a list of particles and the gravity constant. It constructs a force vector scaled against the world Z direction. The Apply( ) method accumulates gravity for each of the particles associated.

Python Definition

class Gravity:
    def __init__( self, particles, acceleration = 9.8 ):
        self.Particles    = particles
        self.Acceleration = Vector3d( 0.0, 0.0, -acceleration )

    def Apply( self ):
        for particle in self.Particles:
            particle.SumForce += self.Acceleration * particle.Mass

Viscous Drag Force

Damping represents a force exerted by a medium against a particle’s motion eg. moving through air is easier compared to water. The natural phenomenon taking place in reality involves quite some complicated physics including notions of a medium’s density, the shape of the object and its velocity.

Python Definition

class Damping:
    def __init__( self, particles, scaling = 1.0 ):
        self.Particles = particles
        self.Scaling   = scaling

    def Apply( self ):
        for particle in self.Particles:
            particle.SumForce += particle.Velocity * -self.Scaling

Here we take a massive shortcut by expressing this notion as a force linearly proportional against a particle’s velocity and packing all material and geometry information into a scaling factor. The objective is to have a force that dissipates energy so eventually the simulation can come down to a halt. In fact, adding drag also helps with issues pertaining numerical drift ie. accumulated errors due to inexact nature of floating point real numbers.

Penalty Constraints

Before putting everything together we need to introduce a form of physical constraints. The goal is to express various boundary conditions, such as a ground plane particles can bounce against, or in general rules that particles must obey unconditionally. This is the exact opposite of forces which describe desired behaviors.

There many ways to approach the subject matter, some more faithful to reality, some less. Here we will take the later, also known as penalty or projection constraints, because they are easy to implement and if used with care they can work perfectly fine.

The issue is that we will violate the laws of physics by heavy handedly putting particles back where they belong, ignoring the implied energy input that this would normally take. Just as an indication, to properly implement physically faithful constraints, we need to impose energy balance, ie. retain zero virtual work done at all times, which implies heavy use of calculus as well as tons of alchemy and black magic. Let’s keep it simple.

Ground Plane

The ground plane represents a horizontal surface at world Z = 0, where particles can bounce off. The constraint requires a list of particles and an energy loss factor. In the application method, we evaluate if a particle has crossed the zero plane and if so, we inverts it Z position coordinate, invert its velocity vector Z component and scale it down by the amount of loss.

Python Definition

class GroundPlane:
    def __init__( self, particles, loss = 1.0 ):
        self.Particles = particles
        self.Loss = loss
        
    def Apply( self ):
        for particle in self.Particles:
            if( particle.Position.Z < 0 ):
                particle.Position.Z *= -1
                particle.Velocity.Z *= -1
                particle.Velocity *= self.Loss

The way to consider this is as follows: We notice, after integration, that we are in the false position (Pf) ie. we have penetrated the ground. This implies that we missed a collision incident that took place at point (Pc). We need to restitute the particle’s position finding (Pr). This can be solved geometrically by reflecting the false position about the world XY plane.

You can check that flipping the sign of the particle’s Z coordinate does the trick (for this simple case of world XY plane). We also apply this to the velocity’s Z component. Finally, if we want to simulate loss of energy due to collision, we can take the loss factor below unity. Notice also, that we can also apply negative loss values, which would make the surface accelerate bouncing particles (think of mario cart boost strips).

Bouncing Particles Setup

Python Definition

from Rhino.Geometry import *
#------------------------------------------------------------------------------- Particle
class Particle:
    def __init__( self, position, velocity, mass ):
        self.Position = position
        self.Velocity = velocity
        self.SumForce = Vector3d.Zero
        self.Mass     = mass
        
    def Display( self ):
        return [self.Position]
#------------------------------------------------------------------------------- Gravity 
class Gravity:
    def __init__( self, particles, acceleration = 9.8 ):
        self.Particles    = particles
        self.Acceleration = Vector3d( 0, 0, -acceleration )
    def Apply( self ):
        for particle in self.Particles:
            particle.SumForce += particle.Mass * self.Acceleration
#------------------------------------------------------------------------------- Damping
class Damping:
    def __init__( self, particles, scaling = 1.0 ):
        self.Particles = particles
        self.Scaling   = scaling
    def Apply( self ):
        for particle in self.Particles:
            particle.SumForce += particle.Velocity * ( -self.Scaling )

#------------------------------------------------------------------------------- Ground
class Ground:
    def __init__( self, particles, loss = 1.0 ):
        self.Particles = particles
        self.Loss = loss
    def Apply( self ):
        for particle in self.Particles:
            if( particle.Position.Z < 0 ):
                particle.Position.Z *= -1
                particle.Velocity.Z *= -1
                particle.Velocity *= self.Loss
#------------------------------------------------------------------------------- Simulation
class Simulation:
    def __init__( self ):
        self.Particles   = []
        self.Forces      = []        
        self.Constraints = []
        
    def Update( self, dt ):        
        for particle in self.Particles:       #-- Zero All Sums of Forces
            particle.SumForce = Vector3d.Zero
            
        for force in self.Forces:             #-- Accumulate Forces
            force.Apply( )
            
        for particle in self.Particles:       #-- Symplectic Euler Integration
            if( particle.Mass == 0 ): continue
            acceleration = particle.SumForce * ( 1.0 / particle.Mass )
            particle.Velocity += acceleration * dt
            particle.Position += particle.Velocity * dt
            
        for constraint in self.Constraints:   #-- Apply Penalty Constraints
            constraint.Apply( )
            
    def KineticEnergy( self ):
        energy = 0.0
        for particle in self.Particles:
            energy += 0.5 * particle.Mass * particle.Velocity * particle.Velocity
        return energy
        
    def PotentialEnergy( self ):
        energy = 0.0
        for particle in self.Particles:
            energy += 9.8 * particle.Mass * particle.Position.Z
        return energy
        
    def Display( self ):
        #-- Geometry
        #--
        geometry = []        
        for particle in self.Particles:
            geometry += particle.Display( )
        
        #-- Messages
        #--
        ke = self.KineticEnergy( )
        pe = self.PotentialEnergy( )
        print( "Kinetic   {0}".format( ke      ) )
        print( "Potential {0}".format( pe      ) )
        print( "Total     {0}".format( ke + pe ) )
        
        return geometry
        
    def BouncingParticles( self ):
        #-- A number of particles along X-Axis with increasing mass
        #--
        for index in range( 10 ): 
            particle = Particle( 
                Point3d( index, 0, 100 ), 
                Vector3d.Zero, index + 1 )
            self.Particles.append( particle )
        
        #-- Setup forces
        #--
        gravity = Gravity( self.Particles )
        self.Forces.append( gravity )
        
        drag = Damping( self.Particles, 0.1 )
        self.Forces.append( drag )
        
        #-- Ground constraint
        #--
        ground = Ground( self.Particles, 0.5 )
        self.Constraints.append( ground )
#------------------------------------------------------------------------------- Entry 
if "simulation" in vars( ) and not reset:
    for iterations in range( 10 ): #-- Non-Visual Updates
        simulation.Update( dt )
    geometry = simulation.Display( )
else:                              #-- Reset & Construction
    simulation = Simulation( )
    simulation.BouncingParticles( )
    geometry = simulation.Display( )

The code above contains all of the pieces together including the necessary modifications to incorporate constraints and drawing/reporting.

Modeling Insights

There are some coding patterns used here that may be useful to consider when developing software that becomes longer than a couple of hundred lines.

  • It is generally a good idea to avoid polluting the global scope. Here, at the script’s entry point, we just minimally deal with one variable, namely simulation, and the script’s input / output parameters. The general rule is keep it as simple and small as possible.
  • It is also good practice to logically compose code into stand-alone units. Here we use classes to package functionality. Also note that no method is beyond a few lines of code. If you find yourself writing lengthy functions, perhaps you should decompose them into semantically reasonable chunks. This also helps with the need for excessive commenting. If functions are named clearly, then code logic reads naturally.
  • Finally, we define conventions regarding how classes behave. For example all forces and constraints implement an Apply( ) method with no parameters and all setup is done at construction. This way we can add all sorts of different type of forces into a common pool and invoke them serially without having to change call style per type of force. This is called polymorphic behavior coupled with dynamic dispatch. In python we approach this by following our own conventions religiously, in more advanced programming languages we define so called interfaces which enforce this type of consistency.

Elasticity Dynamics

Next we will expand the physics engine by introducing springs. Springs capture the behavior of elastic materials, which can take small deformations under load but return to their original shape when the load is released. Those will allow us to create deformable bodies including rope-like curves and cloth-like surfaces.

Python Definition

class Spring:
    def __init__( self, source, target, stiffness = 1.0, restlength = 1.0 ):
        self.Source     = source
        self.Target     = target
        self.Stiffness  = stiffness
        self.RestLength = restlength

    def Apply( self ):
        direction = self.Target.Position - self.Source.Position
        distance = direction.Length        
        magnitude = self.Stiffness * ( self.RestLength - distance )
        force = direction * ( magnitude / distance )
        self.Source.SumForce -= force
        self.Target.SumForce += force

    def Display( self ):
        return [Line( self.Source.Position, self.Target. Position )]

The linear spring force is a binary force between two particles. It follows the Hooke’s law (F = – k ΔL) where (k) is a scaling factor expressing the stiffness of the spring, due to its material and geometry, and (ΔL) its elongation or contraction length. The negative sign denotes that the spring resists deformation. This can be also expressed as (ΔL = Rest Length – Current Length). If the distance between particles equals its rest length then no force is exerted.

The constructor requires the start and end particles, along with the spring’s stiffness and rest length. The force method applies the law by first establishing a vector between the source and target particles, measures their distance and computes the force magnitude. The force vector is the inter-particle direction vector normalized, times the force magnitude. Particle’s sum of forces accumulates this force but their signs are opposite!

Single Spring Setup

Python Definition

from Rhino.Geometry import *
#------------------------------------------------------------------------------- Particle
class Particle:
    def __init__( self, position, velocity, mass ):
        self.Position = position
        self.Velocity = velocity
        self.SumForce = Vector3d.Zero
        self.Mass     = mass
        
    def Display( self ):
        return [self.Position]
#------------------------------------------------------------------------------- Gravity
class Gravity:
    def __init__( self, particles, acceleration = 9.8 ):
        self.Particles    = particles
        self.Acceleration = Vector3d( 0, 0, -acceleration )
        
    def Apply( self ):
        for particle in self.Particles:
            particle.SumForce += particle.Mass * self.Acceleration
    
    def Display( self ):
        return []
#------------------------------------------------------------------------------- Damping 
class Damping:
    def __init__( self, particles, scaling = 1.0 ):
        self.Particles = particles
        self.Scaling   = scaling

    def Apply( self ):
        for particle in self.Particles:
            particle.SumForce += particle.Velocity * ( -self.Scaling )

    def Display( self ):
        return []
#------------------------------------------------------------------------------- Spring
class Spring:
    def __init__( self, source, target, stiffness = 1.0, restlength = 1.0 ):
        self.Source     = source
        self.Target     = target
        self.Stiffness  = stiffness
        self.RestLength = restlength

    def Apply( self ):
        direction = self.Target.Position - self.Source.Position
        distance = direction.Length        
        magnitude = self.Stiffness * ( self.RestLength - distance )
        force = direction * ( magnitude / distance )
        self.Source.SumForce -= force
        self.Target.SumForce += force

    def Display( self ):
        return [Line( self.Source.Position, self.Target. Position )]
#------------------------------------------------------------------------------- Ground
class Ground:
    def __init__( self, particles, loss = 1.0 ):
        self.Particles = particles
        self.Loss = loss
        
    def Apply( self ):
        for particle in self.Particles:
            if( particle.Position.Z < 0 ):
                particle.Position.Z *= -1
                particle.Velocity.Z *= -1
                particle.Velocity *= self.Loss
                
    def Display( self ):
        return []
#------------------------------------------------------------------------------- Simulation
class Simulation:
    def __init__( self ):
        self.Particles   = []
        self.Forces      = []        
        self.Constraints = []
    
    def Analysis( self ):
        for particle in self.Particles:
            particle.SumForce = Vector3d.Zero
            
        for force in self.Forces:
            force.Apply( )
            
    def Integration( self ):
        for particle in self.Particles:
            if( particle.Mass == 0 ): continue
            acceleration = particle.SumForce * ( 1.0 / particle.Mass )
            particle.Velocity += acceleration * dt
            particle.Position += particle.Velocity * dt
            
    def Penalties( self ):
        for constraint in self.Constraints:
            constraint.Apply( )
    
    def Update( self, dt ):        
        self.Analysis( )
        self.Integration( )
        self.Penalties( )
            
    def KineticEnergy( self ):
        energy = 0.0
        for particle in self.Particles:
            energy += 0.5 * particle.Mass * particle.Velocity * particle.Velocity
        return energy
        
    def PotentialEnergy( self ):
        energy = 0.0
        for particle in self.Particles:
            energy += 9.8 * particle.Mass * particle.Position.Z
        return energy
        
    def Display( self ):
        #-- Geometry
        #--
        geometry = []        
        for particle in self.Particles:
            geometry += particle.Display( )
        
        for force in self.Forces:
            geometry += force.Display( )
        
        for constraint in self.Constraints:
            geometry += constraint.Display( )

        #-- Messages
        #--
        ke = self.KineticEnergy( )
        pe = self.PotentialEnergy( )
        print( "Kinetic   {0}".format( ke      ) )
        print( "Potential {0}".format( pe      ) )
        print( "Total     {0}".format( ke + pe ) )
        
        return geometry
        
    def BouncingParticles( self ):
        #-- A number of particles along X-Axis with increasing mass
        #--
        for index in range( 10 ): 
            particle = Particle( 
                Point3d( index, 0, 100 ), 
                Vector3d.Zero, index + 1 )
            self.Particles.append( particle )
        
        #-- Setup forces
        #--
        gravity = Gravity( self.Particles )
        self.Forces.append( gravity )
        
        drag = Damping( self.Particles, 0.1 )
        self.Forces.append( drag )
        
        #-- Ground constraint
        #--
        ground = Ground( self.Particles, 0.5 )
        self.Constraints.append( ground )
        
        return self
        
    def SingleSpring( self ):
        source = Particle( Point3d(  0, 0, 0 ), Vector3d.Zero, 0.0 )
        self.Particles.append( source )
        
        target = Particle( Point3d( 10, 0, 0 ), Vector3d.Zero, 1.0 )
        self.Particles.append( target )
        
        spring = Spring( source, target, 1.0, 5.0 )
        self.Forces.append( spring )
        
        drag = Damping( self.Particles, 0.5 )
        self.Forces.append( drag )    
        
        return self
#------------------------------------------------------------------------------- Entry 
if "simulation" in vars( ) and not reset:
    for iterations in range( 10 ): #-- Non-Visual Updates
        simulation.Update( dt )
    geometry = simulation.Display( )
else:                              #-- Reset & Construction
    simulation = Simulation( ).SingleSpring( )
    geometry = simulation.Display( )

The code is updated with a single spring setup method. All setup methods need to be integrated in the same fashion as part of the simulation class. There are quite a few modifications made along the way in the spirit of improving code quality ie. readability.

  • Motivated by the need to visualize springs, all classes acquired a Display( ) method which is supposed to return a list of geometric entities. Those will be compiled into a master geometry list by the simulation class and then sent out for visualization. For entities we don’t care visualizing we can return an empty list.
  • As exercise perhaps you can try to visualize various entities, for instance the force or velocity vectors of each particle and constraints such the ground plane (hint: use the bounding box of the particles).
  • The simulation steps are broken into semantically readable method calls, namely force analysis and integration, and comments were removed when the code logic is clear.
  • The setup functions now are expected to return (self) such that the entry point can be made a tiny bit neater by call chaining.

Hanging Chain Setup

We create a sequence of particles and link them together with springs. The first and last particles are assigned zero mass to become stationary. All springs have the same rest length and stiffness. The curve approximated by this sequence of particles is a catenary, ie. the curve approximated by a hanging rope or chain. Note how the start and end springs are the most elongated because they need to carry the weight of the rest.

Python Definition

    def HangingChain( self ):
        count = 10
        for index in range( count ):
            mass = 0.0 if( index == 0 or index == count - 1 ) else 1.0
            particle = Particle( Point3d( index, 0, 0 ), Vector3d.Zero, mass )
            self.Particles.append( particle )
        
        for index in range( 1, count ):
            source = self.Particles[index - 1]
            target = self.Particles[index - 0]
            
            spring = Spring( source, target, 15.0, 1.0 )
            self.Forces.append( spring )
        
        gravity = Gravity( self.Particles )
        self.Forces.append( gravity )
        
        drag = Damping( self.Particles, 0.5 )
        self.Forces.append( drag )    
        
        return self

Vault Setup

This setup applies the hanging chain logic in two dimensions. The result is a vault or rubber-like surface. Note that gravity has been inverted and the corners are made stationary.

Python Definition

    def VaultPatch( self ):
        nx = 12
        ny = 8
         
        for iy in range( ny ):
            for ix in range( nx ):
               mass = 0.0 if( ( ix % ( nx - 1 ) == 0 ) and 
                              ( iy % ( ny - 1 ) == 0 )) else 1.0
               particle = Particle( Point3d( ix, iy, 0 ), Vector3d.Zero, mass )
               self.Particles.append( particle )
        
        r = 1.0
        k = 200.0
        for iy in range( 1, ny ):
            for ix in range( 1, nx ):
                pa = self.Particles[( ix - 0 ) + ( iy - 0 ) * nx]
                pb = self.Particles[( ix - 1 ) + ( iy - 0 ) * nx]
                pc = self.Particles[( ix - 0 ) + ( iy - 1 ) * nx]
                pd = self.Particles[( ix - 1 ) + ( iy - 1 ) * nx]
                
                spring = Spring( pa, pc, k, r )
                self.Forces.append( spring )

                spring = Spring( pa, pb, k, r )
                self.Forces.append( spring )
                
                if( ix == 1 ):
                    spring = Spring( pd, pb, k, r )
                    self.Forces.append( spring )
                
                if( iy == 1 ):
                    spring = Spring( pd, pc, k, r )
                    self.Forces.append( spring )
               
        gravity = Gravity( self.Particles, -1 )
        self.Forces.append( gravity )
        
        drag = Damping( self.Particles, 0.5 )
        self.Forces.append( drag )    

        return self

Geometric Model to Simulation Setup

The objective of this section is to demonstrate how to import a geometric model and setup a dead-load simulation scenario. The design is comprised of 132 line segments forming a cantilevered truss-like structure. We will assume that red elements are stationary ie. a cantilever configuration. To recreate the geometry below, copy the model definition text and paste it in Rhino’s command line area. The macro instructions will draw and color-code the generated geometry.

Model Definition

-line 0.0,0.0,10.0 0.0,10.0,10.0 -sellast enter -properties o c o 0,0,0 enter enter
-line 10.0,0.0,10.0 10.0,10.0,10.0 -sellast enter -properties o c o 0,0,0 enter enter
-line 10.0,0.0,0.0 10.0,10.0,0.9 -sellast enter -properties o c o 0,0,0 enter enter
-line 10.0,0.0,0.0 10.0,10.0,10.0 -sellast enter -properties o c o 0,255,255 enter enter
-line 0.0,0.0,10.0 10.0,10.0,10.0 -sellast enter -properties o c o 0,255,0 enter enter
-line 0.0,0.0,0.0 10.0,10.0,0.9 -sellast enter -properties o c o 0,255,0 enter enter
-line 0.0,0.0,0.0 0.0,10.0,0.9 -sellast enter -properties o c o 0,0,0 enter enter
-line 0.0,10.0,0.9 10.0,10.0,0.9 -sellast enter -properties o c o 255,0,255 enter enter
-line 10.0,10.0,0.9 10.0,10.0,10.0 -sellast enter -properties o c o 255,0,255 enter enter
-line 10.0,10.0,10.0 0.0,10.0,10.0 -sellast enter -properties o c o 255,0,255 enter enter
-line 0.0,10.0,10.0 0.0,10.0,0.9 -sellast enter -properties o c o 255,0,255 enter enter
-line 0.0,0.0,0.0 10.0,0.0,0.0 -sellast enter -properties o c o 255,0,0 enter enter
-line 10.0,0.0,0.0 10.0,0.0,10.0 -sellast enter -properties o c o 255,0,0 enter enter
-line 10.0,0.0,10.0 0.0,0.0,10.0 -sellast enter -properties o c o 255,0,0 enter enter
-line 0.0,0.0,10.0 0.0,0.0,0.0 -sellast enter -properties o c o 255,0,0 enter enter
-line 0.0,10.0,10.0 0.0,20.0,10.0 -sellast enter -properties o c o 0,0,0 enter enter
-line 10.0,10.0,10.0 10.0,20.0,10.0 -sellast enter -properties o c o 0,0,0 enter enter
-line 10.0,10.0,0.9 10.0,20.0,1.719 -sellast enter -properties o c o 0,0,0 enter enter
-line 0.0,20.0,10.0 10.0,20.0,1.719 -sellast enter -properties o c o 255,255,255 enter enter
-line 0.0,10.0,0.9 0.0,20.0,1.719 -sellast enter -properties o c o 0,0,0 enter enter
-line 0.0,20.0,1.719 10.0,20.0,1.719 -sellast enter -properties o c o 255,0,255 enter enter
-line 10.0,20.0,1.719 10.0,20.0,10.0 -sellast enter -properties o c o 255,0,255 enter enter
-line 10.0,20.0,10.0 0.0,20.0,10.0 -sellast enter -properties o c o 255,0,255 enter enter
-line 0.0,20.0,10.0 0.0,20.0,1.719 -sellast enter -properties o c o 255,0,255 enter enter
-line 0.0,20.0,10.0 0.0,30.0,10.0 -sellast enter -properties o c o 0,0,0 enter enter
-line 10.0,20.0,10.0 10.0,30.0,10.0 -sellast enter -properties o c o 0,0,0 enter enter
-line 10.0,20.0,1.719 10.0,30.0,2.46429 -sellast enter -properties o c o 0,0,0 enter enter
-line 10.0,20.0,1.719 10.0,30.0,10.0 -sellast enter -properties o c o 0,255,255 enter enter
-line 0.0,20.0,10.0 10.0,30.0,10.0 -sellast enter -properties o c o 0,255,0 enter enter
-line 0.0,20.0,1.719 10.0,30.0,2.46429 -sellast enter -properties o c o 0,255,0 enter enter
-line 0.0,20.0,1.719 0.0,30.0,2.46429 -sellast enter -properties o c o 0,0,0 enter enter
-line 0.0,30.0,2.46429 10.0,30.0,2.46429 -sellast enter -properties o c o 255,0,255 enter enter
-line 10.0,30.0,2.46429 10.0,30.0,10.0 -sellast enter -properties o c o 255,0,255 enter enter
-line 10.0,30.0,10.0 0.0,30.0,10.0 -sellast enter -properties o c o 255,0,255 enter enter
-line 0.0,30.0,10.0 0.0,30.0,2.46429 -sellast enter -properties o c o 255,0,255 enter enter
-line 0.0,30.0,10.0 0.0,40.0,10.0 -sellast enter -properties o c o 0,0,0 enter enter
-line 10.0,30.0,10.0 10.0,40.0,10.0 -sellast enter -properties o c o 0,0,0 enter enter
-line 10.0,30.0,2.46429 10.0,40.0,3.1425039 -sellast enter -properties o c o 0,0,0 enter enter
-line 0.0,40.0,10.0 10.0,40.0,3.1425039 -sellast enter -properties o c o 255,255,255 enter enter
-line 0.0,30.0,2.46429 0.0,40.0,3.1425039 -sellast enter -properties o c o 0,0,0 enter enter
-line 0.0,40.0,3.1425039 10.0,40.0,3.1425039 -sellast enter -properties o c o 255,0,255 enter enter
-line 10.0,40.0,3.1425039 10.0,40.0,10.0 -sellast enter -properties o c o 255,0,255 enter enter
-line 10.0,40.0,10.0 0.0,40.0,10.0 -sellast enter -properties o c o 255,0,255 enter enter
-line 0.0,40.0,10.0 0.0,40.0,3.1425039 -sellast enter -properties o c o 255,0,255 enter enter
-line 0.0,40.0,10.0 0.0,50.0,10.0 -sellast enter -properties o c o 0,0,0 enter enter
-line 10.0,40.0,10.0 10.0,50.0,10.0 -sellast enter -properties o c o 0,0,0 enter enter
-line 10.0,40.0,3.1425039 10.0,50.0,3.759678549 -sellast enter -properties o c o 0,0,0 enter enter
-line 10.0,40.0,3.1425039 10.0,50.0,10.0 -sellast enter -properties o c o 0,255,255 enter enter
-line 0.0,40.0,10.0 10.0,50.0,10.0 -sellast enter -properties o c o 0,255,0 enter enter
-line 0.0,40.0,3.1425039 10.0,50.0,3.759678549 -sellast enter -properties o c o 0,255,0 enter enter
-line 0.0,40.0,3.1425039 0.0,50.0,3.759678549 -sellast enter -properties o c o 0,0,0 enter enter
-line 0.0,50.0,3.759678549 10.0,50.0,3.759678549 -sellast enter -properties o c o 255,0,255 enter enter
-line 10.0,50.0,3.759678549 10.0,50.0,10.0 -sellast enter -properties o c o 255,0,255 enter enter
-line 10.0,50.0,10.0 0.0,50.0,10.0 -sellast enter -properties o c o 255,0,255 enter enter
-line 0.0,50.0,10.0 0.0,50.0,3.759678549 -sellast enter -properties o c o 255,0,255 enter enter
-line 0.0,50.0,10.0 0.0,60.0,10.0 -sellast enter -properties o c o 0,0,0 enter enter
-line 10.0,50.0,10.0 10.0,60.0,10.0 -sellast enter -properties o c o 0,0,0 enter enter
-line 10.0,50.0,3.759678549 10.0,60.0,4.32130747959 -sellast enter -properties o c o 0,0,0 enter enter
-line 0.0,60.0,10.0 10.0,60.0,4.32130747959 -sellast enter -properties o c o 255,255,255 enter enter
-line 0.0,50.0,3.759678549 0.0,60.0,4.32130747959 -sellast enter -properties o c o 0,0,0 enter enter
-line 0.0,60.0,4.32130747959 10.0,60.0,4.32130747959 -sellast enter -properties o c o 255,0,255 enter enter
-line 10.0,60.0,4.32130747959 10.0,60.0,10.0 -sellast enter -properties o c o 255,0,255 enter enter
-line 10.0,60.0,10.0 0.0,60.0,10.0 -sellast enter -properties o c o 255,0,255 enter enter
-line 0.0,60.0,10.0 0.0,60.0,4.32130747959 -sellast enter -properties o c o 255,0,255 enter enter
-line 0.0,60.0,10.0 0.0,70.0,10.0 -sellast enter -properties o c o 0,0,0 enter enter
-line 10.0,60.0,10.0 10.0,70.0,10.0 -sellast enter -properties o c o 0,0,0 enter enter
-line 10.0,60.0,4.32130747959 10.0,70.0,4.83238980643 -sellast enter -properties o c o 0,0,0 enter enter
-line 10.0,60.0,4.32130747959 10.0,70.0,10.0 -sellast enter -properties o c o 0,255,255 enter enter
-line 0.0,60.0,10.0 10.0,70.0,10.0 -sellast enter -properties o c o 0,255,0 enter enter
-line 0.0,60.0,4.32130747959 10.0,70.0,4.83238980643 -sellast enter -properties o c o 0,255,0 enter enter
-line 0.0,60.0,4.32130747959 0.0,70.0,4.83238980643 -sellast enter -properties o c o 0,0,0 enter enter
-line 0.0,70.0,4.83238980643 10.0,70.0,4.83238980643 -sellast enter -properties o c o 255,0,255 enter enter
-line 10.0,70.0,4.83238980643 10.0,70.0,10.0 -sellast enter -properties o c o 255,0,255 enter enter
-line 10.0,70.0,10.0 0.0,70.0,10.0 -sellast enter -properties o c o 255,0,255 enter enter
-line 0.0,70.0,10.0 0.0,70.0,4.83238980643 -sellast enter -properties o c o 255,0,255 enter enter
-line 0.0,70.0,10.0 0.0,80.0,10.0 -sellast enter -properties o c o 0,0,0 enter enter
-line 10.0,70.0,10.0 10.0,80.0,10.0 -sellast enter -properties o c o 0,0,0 enter enter
-line 10.0,70.0,4.83238980643 10.0,80.0,5.29747472385 -sellast enter -properties o c o 0,0,0 enter enter
-line 0.0,80.0,10.0 10.0,80.0,5.29747472385 -sellast enter -properties o c o 255,255,255 enter enter
-line 0.0,70.0,4.83238980643 0.0,80.0,5.29747472385 -sellast enter -properties o c o 0,0,0 enter enter
-line 0.0,80.0,5.29747472385 10.0,80.0,5.29747472385 -sellast enter -properties o c o 255,0,255 enter enter
-line 10.0,80.0,5.29747472385 10.0,80.0,10.0 -sellast enter -properties o c o 255,0,255 enter enter
-line 10.0,80.0,10.0 0.0,80.0,10.0 -sellast enter -properties o c o 255,0,255 enter enter
-line 0.0,80.0,10.0 0.0,80.0,5.29747472385 -sellast enter -properties o c o 255,0,255 enter enter
-line 10.0,10.0,0.9 0.0,20.0,1.719 -sellast enter -properties o c o 0,255,0 enter enter
-line 10.0,10.0,10.0 0.0,20.0,10.0 -sellast enter -properties o c o 0,255,0 enter enter
-line 10.0,30.0,2.46429 0.0,40.0,3.1425039 -sellast enter -properties o c o 0,255,0 enter enter
-line 10.0,50.0,3.759678549 0.0,60.0,4.32130747959 -sellast enter -properties o c o 0,255,0 enter enter
-line 10.0,70.0,4.83238980643 0.0,80.0,5.29747472385 -sellast enter -properties o c o 0,255,0 enter enter
-line 10.0,30.0,10.0 0.0,40.0,10.0 -sellast enter -properties o c o 0,255,0 enter enter
-line 10.0,50.0,10.0 0.0,60.0,10.0 -sellast enter -properties o c o 0,255,0 enter enter
-line 10.0,70.0,10.0 0.0,80.0,10.0 -sellast enter -properties o c o 0,255,0 enter enter
-line 10.0,10.0,10.0 10.0,20.0,1.719 -sellast enter -properties o c o 0,255,255 enter enter
-line 10.0,30.0,10.0 10.0,40.0,3.1425039 -sellast enter -properties o c o 0,255,255 enter enter
-line 10.0,50.0,10.0 10.0,60.0,4.32130747959 -sellast enter -properties o c o 0,255,255 enter enter
-line 10.0,70.0,10.0 10.0,80.0,5.29747472385 -sellast enter -properties o c o 0,255,255 enter enter
-line 0.0,0.0,0.0 0.0,10.0,10.0 -sellast enter -properties o c o 0,255,255 enter enter
-line 0.0,20.0,1.719 0.0,30.0,10.0 -sellast enter -properties o c o 0,255,255 enter enter
-line 0.0,40.0,3.1425039 0.0,50.0,10.0 -sellast enter -properties o c o 0,255,255 enter enter
-line 0.0,60.0,4.32130747959 0.0,70.0,10.0 -sellast enter -properties o c o 0,255,255 enter enter
-line 0.0,10.0,10.0 0.0,20.0,1.719 -sellast enter -properties o c o 0,255,255 enter enter
-line 0.0,30.0,10.0 0.0,40.0,3.1425039 -sellast enter -properties o c o 0,255,255 enter enter
-line 0.0,50.0,10.0 0.0,60.0,4.32130747959 -sellast enter -properties o c o 0,255,255 enter enter
-line 0.0,70.0,10.0 0.0,80.0,5.29747472385 -sellast enter -properties o c o 0,255,255 enter enter
-line 10.0,10.0,10.0 0.0,10.0,0.9 -sellast enter -properties o c o 255,255,255 enter enter
-line 10.0,30.0,10.0 0.0,30.0,2.46429 -sellast enter -properties o c o 255,255,255 enter enter
-line 10.0,50.0,10.0 0.0,50.0,3.759678549 -sellast enter -properties o c o 255,255,255 enter enter
-line 10.0,70.0,10.0 0.0,70.0,4.83238980643 -sellast enter -properties o c o 255,255,255 enter enter
-line 0.0,10.0,0.9 10.0,0.0,0.0 -sellast enter -properties o c o 0,255,0 enter enter
-line 0.0,10.0,0.9 10.0,20.0,1.719 -sellast enter -properties o c o 0,255,0 enter enter
-line 0.0,30.0,2.46429 10.0,20.0,1.719 -sellast enter -properties o c o 0,255,0 enter enter
-line 0.0,30.0,2.46429 10.0,40.0,3.1425039 -sellast enter -properties o c o 0,255,0 enter enter
-line 0.0,50.0,3.759678549 10.0,40.0,3.1425039 -sellast enter -properties o c o 0,255,0 enter enter
-line 0.0,50.0,3.759678549 10.0,60.0,4.32130747959 -sellast enter -properties o c o 0,255,0 enter enter
-line 10.0,60.0,4.32130747959 0.0,70.0,4.83238980643 -sellast enter -properties o c o 0,255,0 enter enter
-line 0.0,70.0,4.83238980643 10.0,80.0,5.29747472385 -sellast enter -properties o c o 0,255,0 enter enter
-line 10.0,0.0,10.0 0.0,10.0,10.0 -sellast enter -properties o c o 0,255,0 enter enter
-line 0.0,10.0,10.0 10.0,20.0,10.0 -sellast enter -properties o c o 0,255,0 enter enter
-line 10.0,20.0,10.0 0.0,30.0,10.0 -sellast enter -properties o c o 0,255,0 enter enter
-line 0.0,30.0,10.0 10.0,40.0,10.0 -sellast enter -properties o c o 0,255,0 enter enter
-line 10.0,40.0,10.0 0.0,50.0,10.0 -sellast enter -properties o c o 0,255,0 enter enter
-line 0.0,50.0,10.0 10.0,60.0,10.0 -sellast enter -properties o c o 0,255,0 enter enter
-line 10.0,60.0,10.0 0.0,70.0,10.0 -sellast enter -properties o c o 0,255,0 enter enter
-line 0.0,70.0,10.0 10.0,80.0,10.0 -sellast enter -properties o c o 0,255,0 enter enter
-line 0.0,10.0,10.0 10.0,10.0,0.9 -sellast enter -properties o c o 255,255,255 enter enter
-line 10.0,20.0,10.0 0.0,20.0,1.719 -sellast enter -properties o c o 255,255,255 enter enter
-line 10.0,40.0,10.0 0.0,40.0,3.1425039 -sellast enter -properties o c o 255,255,255 enter enter
-line 0.0,50.0,10.0 10.0,50.0,3.759678549 -sellast enter -properties o c o 255,255,255 enter enter
-line 10.0,60.0,10.0 0.0,60.0,4.32130747959 -sellast enter -properties o c o 255,255,255 enter enter
-line 0.0,70.0,10.0 10.0,70.0,4.83238980643 -sellast enter -properties o c o 255,255,255 enter enter
-line 10.0,80.0,10.0 0.0,80.0,5.29747472385 -sellast enter -properties o c o 255,255,255 enter enter
-line 10.0,30.0,2.46429 0.0,30.0,10.0 -sellast enter -properties o c o 255,255,255 enter enter

We need to update the graph definition by incorporating a new (list of GUID) input parameter. The modeled lines are collected via a geometry, not a curve, input block. The reason for importing geometric entities from the 3d model via Global Unique Identifiers is because apart from their geometry we also need their display attributes ie. color, to discern stationary from dynamic elements. So this setup also presents direct geometry processing using document object attributes.

Python Definition

    def InitFromEdges( self, edges ):
        #-- Reference Active Document
        #--
        from System.Drawing import Color
        from Rhino import RhinoDoc
        document = RhinoDoc.ActiveDoc
        
        #-- Find Particle by Position
        #--
        def IndexOfParticle( point, epsilon = 0.001 ):
            count = len( self.Particles )
            for index in range( count ):
                particle = self.Particles[index]
                distance = point.DistanceTo( particle.Position )
                if( distance <= epsilon ): return index
            return -1
        
        def RetrieveOrAppendUnique( point ):
            index = IndexOfParticle( point )
            if( index < 0 ):
                particle = Particle( point, Vector3d.Zero, 1.0 )
                index = len( self.Particles )
                self.Particles.append( particle )
            return self.Particles[index]
        
        for edge in edges:
            #-- Rhino Object from GUID
            #--
            object = document.Objects.Find( edge )
            
            #-- Convert Points to Unique Particles
            #--
            p = object.Geometry.PointAtStart
            q = object.Geometry.PointAtEnd

            P = RetrieveOrAppendUnique( p )
            Q = RetrieveOrAppendUnique( q )
            
            #-- Red Edges -> Fixed Particles
            #--
            attributes = object.Attributes
            color = attributes.ObjectColor
            if( color == Color.Red ):
                P.Mass = Q.Mass = 0.0
                
            #-- Spring already at Rest
            #--
            spring = Spring( P, Q, 1000.0, p.DistanceTo( q ) )
            self.Forces.append( spring )
            
        gravity = Gravity( self.Particles, 10 )
        self.Forces.append( gravity )
        
        drag = Damping( self.Particles, 0.5 )
        self.Forces.append( drag )   
        
        return self

The code above requires a bit of explanation: We aim to convert incoming curves into springs and their start/end points into particles. However, if we just directly convert and collect particles from each line separately, they will not properly form common joints, where multiple lines meet, because we are collecting lots of duplicates. Instead, we need to create a pool of common unique particle nodes first and then associate springs with those.

This operation is known as topology extraction or node/vertex welding. The strategy is as follows: for each curve we will extracts its points and we will try to look up in the existing particles’ list, if there is any other particle already within very close distance away ie. within tolerance. If there is one, we will use that and if not we will create and insert a new in the particles list. This approach ensures that we create particles uniquely. In this fashion we can incrementally, convert all points into unique particle nodes and all lines into springs.

The rest of the code just retrieves the modeling object from the currently opened document, by its unique id, extracts its geometry and display attributes, and fixes by zeroing the mass of red colored particles.

Something to experiment with this code is: (a) the ratios between masses, stiffness, gravity and time step and (b) the truss’ internal bracing pattern. In the former case, it is surprising easy to have the simulation literally blow up due to instability. This is unfortunately common in scenarios where we need high stiffness to restrain deflection and a key reason for finite elements structural analysis using fully implicit integration scheme ie. stiffness matrix inversion. Playing with the bracing pattern is also interesting as it shows how easy is to acquire an off-axial bias, ie. twist, if the design is asymmetrically braced.

N-body Simulations

In this section we will develop a simple simulation strategy for planetary motion. The key concepts conveyed are: (a) creating forces between all possible particle pairs, (b) creating anti-explosion mechanisms, (c) creating geometry from particles’ motion.

Attraction Force

For each unique pair of particles, we will compute a force based on planetary motion physics (F = k m M / d^2) where (k) is an arbitrary gravitational constant or scaling factor, (m) and (M) are the masses of the two particles involved, and (d) is the distance between them. Of course the force is proportional to the inverse square of the distance, and here we use the third power because we also need to normalize the vector.

Python Definition

class Attraction:
    def __init__( self, particles, scale = 1.0 ):
        self.Particles = particles
        self.Scale     = scale

    def Apply( self ):
        for i in range( len( self.Particles ) - 1 ):
            p = self.Particles[i]

            for j in range( i + 1, len( self.Particles ) ):
                q = self.Particles[j]

                force = q.Position - p.Position
                force *= self.Scale * p.Mass * q.Mass / ( force.Length ** 3 )

                p.SumForce += force
                q.SumForce -= force

    def Display( self ):
        return []

The implementation looks not unlike a simple sorting algorithm, where we compare one element of a list with all others. Here we want to compute all pair particles once. Note that our force method became quadratic in terms of time complexity (n^2), ie. we have a pair of nested loops, unlike other forces seen earlier, which only require one pass through the list of affected particles (n). As a results, it is a rather expensive type of force which will slow down the simulation substantially as the number of particles used increases.

Python Definition

    def SolarSystem( self ):
        sun = Particle( Point3d.Origin, Vector3d.Zero, 10000.0 )
        self.Particles.append( sun )
        
        earth = Particle( Point3d( 100, 0, 0 ), Vector3d( 0, 10, 0 ), 10 )
        self.Particles.append( earth )

        attraction = Attraction( self.Particles )
        self.Forces.append( attraction )
        
        return self

The setup above shows how we can simulate planetary orbital motion between the sun and earth using only an initial tangential velocity and appropriate distance and masses. In fact it is quite fascinating to play with those parameters and observe that if the velocity is too low, earth falls into the sun, while if it is too high, it escapes sun’s gravity. Adding the moon is quite a challenge which may worth giving it a shot ie. to appreciate the universe.

Velocity Clamping

Apart from being slow, this simulation type also tends to explode. The reason is that attraction causes particles to come together, whereby their distance gets closer and closer to zero. Now, because forces are proportional to the inverse of the square distance, at near zero distance, forces become nearly infinite. Increasing the drag force will not solve this issue and it will make the simulation sluggish. We need to pull the breaks using a constraint that limits particles’ velocity to a maximum value.

class SpeedTrap:
    def __init__( self, particles, max = 100.0 ):
        self.Particles   = particles
        self.MaxVelocity = max
        
    def Apply( self ):
        for particle in self.Particles:
            velocity = particle.Velocity.Length
            if( velocity > self.MaxVelocity ):
                particle.Velocity *= self.MaxVelocity / velocity
                
    def Display( self ):
        return []

The ideal position to apply the max velocity logic is in fact not in the form of a constraint applied after integration, but actually right before integration, where we can also limit the amount of maximum force before we update the positions. After integration it does not make sense to limit forces anymore.

If the notion of a velocity limit is not intuitive enough, we can change it to distance limit quite easily, because from velocity we can compute the displacement over a time step (this is left as an exercise). Another hackish way to overcome such explosion issue to modify the force by adding a safety term ie. from (F = k m M / d^2) use (F = k m M / ( d^2 + e )), where (e) is a small number that prevents the denominator from ever becoming zero.

World Boundaries

Another way to deal with explosions is to create a world such that particles cannot escape its limits. The following constraint establishes a bounding box of certain size. Particles found outside of the box are clamped, ie. projected, onto it and their velocities are flipped. Optionally, we can enforce a energy reduction penalty for particles that tried to escape the world’s boundaries.

Python Definition

class BlockWorld:
    def __init__( self, particles, size = 1.0, decay = 1.0 ):
        self.Particles = particles
        self.Size      = size
        self.Decay     = decay
    
    def Clamp( self, value ):
        if( value < 0.0 ): 
            return 0.0
        if( value > self.Size ):
            return self.Size
        return value
    
    def Apply( self ):
        for particle in self.Particles:
            bounced = False
            
            for index in range( 3 ):
                coord = particle.Position[index]
                clamp = self.Clamp( coord )
                
                if( coord != clamp ):
                    particle.Position[index]  = clamp
                    particle.Velocity[index] *= -1
                    bounced = True
                    
            if( bounced != 0 ):
                particle.Velocity *= self.Decay
                
    def Display( self ):
        return [BoundingBox( 0, 0, 0, self.Size, self.Size, self.Size )]    

The setup code for this boxed planets world is below. Note how the particles go ballistic once attracted onto one another and reach a near zero distance.

Python Definition

    def Multibody( self ):
        from random import *
        scale = 100.0
        count = 32
        for index in range( count ):
            particle = Particle( 
                Point3d( scale * random( ), 
                         scale * random( ), 
                         scale * random( ) ),
                Vector3d.Zero, 1 + random( ) )
            self.Particles.append( particle )
        
        attraction = Attraction( self.Particles )
        self.Forces.append( attraction )
        
        drag = Damping( self.Particles, 0.1 )
        self.Forces.append( drag )           
        
        world = BlockWorld( self.Particles, scale, 1.0 )
        self.Constraints.append( world )           

        return self

Looped World

Another type of world constraint is what is known as a PacMan or toroidal form. The key difference between BlockWorld and TorusWorld is that the later has no boundaries. Particles live in a 3D looped environment, where exiting from one side, they enter on the opposite ie. like PacMan but in three directions. The implementation is quite simple as we use modular arithmetic to wrap around in all three cardinal directions.

Python Definition

class TorusWorld:
    def __init__( self, particles, size = 1.0 ):
        self.Particles = particles
        self.Size      = size
        
    def Apply( self ):
        for particle in self.Particles:
            particle.Position.X %= self.Size
            particle.Position.Y %= self.Size
            particle.Position.Z %= self.Size

    def Display( self ):
        return [BoundingBox( 0, 0, 0, self.Size, self.Size, self.Size )]

The result is quite visually interesting but in terms of simulation it tends to mess a lot of things up because many forces depend on Euclidean distance which this constraints grossly violates.

Spherical World

Yet another type of world we can create is a spherical surface where particles cannot escape. This can be generalized to any kind of surface, including planes or even Nurbs surfaces. In each time step we just need to project particles onto the target surface. Note that projecting points on a sphere or a plane is much more computationally efficient than using the ClosestPoint( ) method of a complex surface.

Python Definition

class SphereWorld:
    def __init__( self, particles, radius = 1.0 ):
        self.Particles = particles
        self.Radius    = radius
        
    def Apply( self ):
        for particle in self.Particles:
            distance = particle.Position.DistanceTo( Point3d.Origin )
            particle.Position *= self.Radius / distance

    def Display( self ):
        return [Sphere( Point3d.Origin, self.Radius )]

The setup for the this example is below. Note that the initial particle positions are randomly selected. They are not even on the sphere but the spherical world constraint will take care of this upon first evaluation.

Python Definition

    def Spherical( self ):
        from random import *
        scale = 100.0
        count = 64
        for index in range( count ):
            particle = Particle( 
                Point3d( scale * ( 1 - 2 * random( ) ), 
                         scale * ( 1 - 2 * random( ) ), 
                         scale * ( 1 - 2 * random( ) ) ),
                Vector3d( random( ) * 0, 
                          random( ) * 10, 
                          random( ) * 0 ), 
                10 )
            self.Particles.append( particle )
        
        attraction = Attraction( self.Particles )
        self.Forces.append( attraction )
        
        world = SphereWorld( self.Particles, scale )
        self.Constraints.append( world )           

        return self

Particle Tracers

What we wish to accomplish here, is to plot the sequence of positions of all particles overtime. This will produce traces of their trajectories. To make this happen we need to record points so we can eventually visualize them. Because each trajectory is associated with each particle, it makes sense to add a new property in the particle’s definition. We will call this the position history of each particle. We initialize the position history list, with the particle’s initial position.

Python Definition

class Particle:
    def __init__( self, position, mass = 1.0 ):
        self.Position = position
        self.Velocity = Vector3d.Zero
        self.SumForce = Vector3d.Zero
        self.Mass     = mass
        self.History  = [position]

In addition, we need somehow to trigger the recording of particle positions. This may be right after a particle has been updated. Conveniently we can craft this path tracing recorder as a constraint! The implementation below is quite elaborate so here are some key features: (a) We will only record a fixed number of points to avoid overloading memory, (b) We will record intermittently to avoid capturing too many points too close together, (c) We will not record points too close to one another to avoid creating degenerate polylines, and (d) We will reset the path if the particle makes a huge jump such as crossing a toroid world boundaries.

Python Definition

class PathTrace:
    def __init__( self, particles, count = 64, cycle = 10, min = 0.01, max = 10.0 ):
        self.Clock     = 0
        self.Cycle     = cycle
        self.Particles = particles
        self.Minimum   = min
        self.Maximum   = max
        self.Count     = count
        
    def Apply( self ):
        self.Clock += 1
        if( self.Clock % self.Cycle != 0 ):
            return
            
        for particle in self.Particles:
            distance = particle.Position.DistanceTo( particle.History[-1] )
            if( distance < self.Minimum ): #-- Too Close -> Ignore
                continue
                
            if( distance > self.Maximum ): #-- Too Far   -> Reset
                particle.History = [particle.Position]
                continue
                
            particle.History.append( particle.Position )
            if( len( particle.History ) > self.Count ):
                particle.History.pop( 0 )
                    
    def Display( self ):
        geometry = []
        for particle in self.Particles:
            if( len( particle.History ) > 2 ):
                geometry.append( PolylineCurve( particle.History ) )
        return geometry

The formation of couplings, clusters, bouncing around the world’s boundaries, the attraction between remote particles are now very vivid. The geometry created by particle trajectories is quite remarkable, if you consider that drawing such curves by standard modeling tools is almost impossible, or rather highly impractical. Physics makes it elegant and effortless in a way. As an exercise you may try change the particle’s visualization from a point to an arrow or paper plane (using its velocity) or use all particles for generate an evolving mesh or surface.

Swarm Intelligence

In the section we will look at swarm intelligence which emulates flocks of birds, schools of fish, etc. The key reference for this section is the paper “Steering Behaviors For Autonomous Characters” from GDC 1999 by Craig W. Reynolds [📚]. The key insight is that what we observe as complex and/or intelligent crowd behaviors can the result of actually very simple rules. There are three concepts we will implement to achieve this nature-inspired system: (a) Cohesion: particles want to be close to one another, (b) Separation: particles want to be sufficiently spaced apart, and (c) Alignment: particles want to flow along with one another.

All of those behaviors require a notion of proximity: Rules only apply for neighborhoods of particles. This allows for the formation of clusters. To implement this we will define the notion of nearest neighbors such that only particles up to some maximum distance can affect one another. Because of this requirement the computation becomes significantly more involved, even though the code is not too complicated.

Cohesion Force

It is a force that attracts particles not unlike the planetary force seen earlier except particles are not directly attracted to one another but by the center of their cluster of nearest neighbors. To find a cluster’s center we need to do the following: for each particle we will examine all other particles and if they are within close range we will accumulate their position such that eventually we can compute their geometric average ie. sum of all positions divided by their number. Then we will create a vector from the particle to the geometric center and scale it to create the required attraction force. Here I do not scale the force by the inverse of the square distance but just a constant scaling factor ie. anything goes.

Python Definition

class Cohesion:
    def __init__( self, particles, scale = 1.0, range = 1.0 ):
        self.Particles = particles
        self.Scale     = scale
        self.Range     = range
        
    def Apply( self ):
        for p in self.Particles:
            point = Point3d.Origin
            count = 0
            
            for q in self.Particles:
                if( p == q ): 
                    continue
                    
                vector = p.Position - q.Position
                length = vector.Length
                if( length > self.Range ): 
                    continue
                    
                point += q.Position
                count += 1
                
            if( count == 0 ):
                continue
            point *= 1.0 / count

            force = point - p.Position
            force.Unitize( )

            p.SumForce += force * self.Scale
            
    def Display( self ):
        return []

Separation Force

Separation force is not very different than cohesion. Here we look at again each particle and for all its nearest neighbors we accumulate vectors away from them. Eventually, we unitize the aggregate force vector, which is like averaging seen before, and scale appropriately. Again, the magnitude of the force is linear without weight scaling but feel free to change eg. particles with bigger mass can be scarier.

Python Definition

class Separation:
    def __init__( self, particles, scale = 1.0, range = 1.0 ):
        self.Particles = particles
        self.Scale     = scale
        self.Range     = range

    def Apply( self ):
        for p in self.Particles:
            force = Vector3d.Zero
            count = 0

            for q in self.Particles:
                if( p == q ): 
                    continue

                vector = p.Position - q.Position
                length = vector.Length
                if( length > self.Range ): 
                    continue

                force += vector
                count += 1

            if( count == 0 ):
                continue
            force.Unitize( )

            p.SumForce += force * self.Scale

    def Display( self ):
        return []

Alignment Force

To encourage particles to flow along similar paths we will use their velocity instead of position. In exactly the same fashion as with separation, we will average all nearest neighbors’ velocities and use this direction as force motivating a particle to also move along the aggregate flow direction.

class Alignment:
    def __init__( self, particles, scale = 1.0, range = 1.0 ):
        self.Particles = particles
        self.Scale     = scale
        self.Range     = range

    def Apply( self ):
        for p in self.Particles:
            angle = Vector3d.Zero
            count = 0

            for q in self.Particles:
                if( p == q ): 
                    continue

                vector = p.Position - q.Position
                length = vector.Length
                if( length > self.Range ): 
                    continue

                angle += q.Velocity
                count += 1

            if( count == 0 ):
                continue

            angle.Unitize( )

            p.SumForce += angle * self.Scale

    def Display( self ):
        return []

The setup code for this experiment is below. It is possible to replace the definition of the cohesion, separation and alignment forces with global variables and change the ratios of those behaviors during simulation. This makes it easier to visualize the effect on the swarm.

    def Swarms( self ):
        from random import *
        scale = 100
        count = 32
        for index in range( count ):
            particle = Particle( 
                Point3d( scale * random( ), 
                         scale * random( ), 
                         scale * random( ) ),
                Vector3d( random( ), 
                          random( ), 
                          random( ) ), 10.0 )
            self.Particles.append( particle )
        
        self.Forces.append( Cohesion   ( self.Particles, 0.25, 30 ) )
        self.Forces.append( Separation ( self.Particles, 0.30, 30 ) )
        self.Forces.append( Alignment  ( self.Particles, 0.10, 30 ) )
        self.Forces.append( Damping    ( self.Particles, 0.10     ) )

        self.Constraints.append( BlockWorld( self.Particles, scale ) )
        self.Constraints.append( PathTrace( self.Particles, 12 ) )
        
        return self 

Conclusions

This kind of concludes this lengthy post. There are tons of more interesting experiments to attempt using dynamics. Those will be progressively implemented and added in the near future. Examples include using the physics engine to: (a) solve geometric/mathematical problems without writing equations or modeling, (b) geometry processing such as smoothing, (c) create fire, snow, rain, water and fluids, (d) rationalizing building geometry, (e) creating origami structures, (f) creating agent based simulations such as people flow, (h) form-finding with membranes and thin films and many many more.

The goals was to demonstrate that with a couple of hundred lines of code we can create a system that blends knowledge from math, physics and computation and can be used for all sorts of fun and exciting applications.

References

  1. Hecker, C. (1996) Physics the Next Frontier [📚]
  2. Witkin, A. and Baraff, D. (1998) Physically Based Modeling [📚]
  3. Cline, M. (1999) Rigid Body Simulation with Contact and Constraints [📚]
  4. Omelyan, I. (1999) A new approach to the numerical integration of rotational motion for interacting rigid bodies [📚]
  5. Chin, S. and Kidwell, D. (2000) Higher Order Force Gradient Symplectic Algorithms [📚]
  6. Smith, R. (2002) How to make new joints in ODE [📚]
  7. Hairer E., Lubich, C., Wanner, G. (2003) Geometric Numerical Intergration Illustrated by the Stormer-Verlet Method [📚]
  8. Jakobsen, T. (2003) Advanced Character Physics [📚]
  9. Erleben, K. (2004) Stable, Robust, and Versatile Multibody Dynamics Animation. [📚]
  10. Catto, E. (2005) Iterative Dynamics with Temporal Coherence [📚]
  11. Muller, M., Heidelberger, B., Teschner, M., Gross, M. (2005) Meshless Deformations Based on Shape Matching [📚]
  12. McLachlan, R. and Quispel, R. (2006) Geometric Integrators for ODEs. [📚]
  13. Muller, M., Heidelberger, B., Hennix, M., Ratcliff, J. (2006) Position Based Dynamics [📚]
  14. Van Verth, J. (2011) Numerical Integration. [📚]
  15. Catto, E. (2007) Modeling and Solving Constraints. [📚]
  16. Daniel Chappuis (2013) Constraints Derivation for Rigid Body Simulation in 3D. [📚]
  17. Tonge R. (2013) Iterative Rigid Body Solvers [📚]
  18. Young, P. (2013) The leapfrog method and other “symplectic” algorithms for integrating Newton’s laws of motion [📚]
  19. Erwin Coumans (2014) Exploring MLCP solvers and Featherstone. [📚]
  20. Catto, E. (2014) Numerical Methods. [📚]
  21. Tomcin, R., Sibbing, D., Kobbelt, L. (2014) Efficient Enforcement of Hand Articulation Constraints [📚]
  22. Tamis, M. (2015) Comparison between Projected Gauss-Seidel and Sequential Impulse Solvers for Real-Time Physics Simulations [📚]