Approach:
To build this constrained particle system we must represent the system of particles in terms of vectors, matrices and constraints. Once we have this representation we can find how external forces interact with the constrained particle system prior to using euler step integration to move the particles.
Representation:
The state of the particle system is contained in the vector q. q is a column vector containing the x and y components of each particle such that:
q = [x1,y1,x2,y2...xn,yn]^
(^ denotes transpose ) Constraints can be expressed in terms of function that are zero when the constraints are satisfied:
C(q) = 0
These constraints can be compiled into constraint vector C such that :
C(q) = [C1(q)..Cn(q)]^ = 0
Dynamics
From this representation we can understand how particles move under the influence of forces. The forces acting on the system are compiled into a single vector Q:
Q = [ Fx1,Fy1,Fx2,Fy2,...Fxn,Fyn]^
Where Fxi,Fyi are the x and y components of the force acting on the ith particle. Since C is a function of q it describes a mapping from q (position) space to C space. In C-space we consider only the subset of loci that C(q) =0 as the space of valid particle locations or constraint space. Ideal we would like the particle system to stay in constraint space, however we can allow some fluctuations in and out of constraint space as long as the particles system state tends to a point on constraint space. This motion in constraint space can thus be described as a higher dimensional harmonic oscillator:
(1) C'' + ksC + kdC' = 0
Where ks is the spring coefficient and kd is the damping term and (') denotes a time derivative. To solve the dynamics of particles we need to relate this oscillator to q-space. This is done using the chain rule:
(2) C' = Jq'
Where J = dC/dq is the Jacobean. C'' can then be represent as:
(3) C'' = J'q'+Jq'' = J'q' + JW(Qvalid)
In (3) we related particle acceleration q'' to newtonian forces with W, the inverse mass matrix- a diagonal matrix whose elements are 1/mass. This just follows from F = ma ==> a =F/m. Note that Qvalid is a force that acts tangentially to constraint space. In general however a force vector will act in arbitrary directions thus Qvalid can be expressed as:
(4) Qvalid = Qnet +Qhat
Where Qhat is a corrective force that cancels out any illegal components of Qnet. Qhat must act perpendicularly to constraint space as it must cancel any components of Qnet that drive q from constraint space. Thus Qhat must be in the space spanned by the rows of J since they are vectors locally perpendicular to constraint space. Thus:
(5) Qhat = J^(lamda)
Combining (1) ... (5) gives:
(6) (JWJ^) lamda = -J'q'=JWQ-ksC -kdC'.
With this equation we can solve for lamda and find Qhat with (5) to cancel any forces acting in directions violation acting out of constraint space.
Constraint Used:
In this simulation we used a distance constraint to keep two particles a fixed distance from each other and a nail constraint to keep a particle in a fixed location:
Distance constraint (^ denotes transpose) :
The necessary equations conserning the distance constrains are:
Cdist = ( xa-xb)^ (xa-xb).
dCdist/dxa = 2(xa-xb)
dCdist/dxb = -2(xa-xb)
Cdist'= 2(va-vb)^ (xa-xb)
dCdist'/dxa = 2(va-vb)
dCdist'/dxb = -2(va-vb)
Nail:
Here a particle has a nail at location p:
Cnail = ( xa-p)^ (xa-p).
dCnail/dxa = 2(xa-p)
Cnail'= 2(va)^ (xa-p)
dCnail'/dxa = 2(va)
External Forces:
The external force used in this simulation include that of gravity, damping,and spring forces. Gravity is a constant force applied in the y direction on each particle. Damping acts as -kdq' to continually slow the particles. We use damped hookian springs forces acting on points in this system. A spring attached to 2 particles a and b will act on a according to :
Fa = [-ks(|L| - r) + ks L'*L/|L||] L/|L|
Where L = pa - pb; and L' = va - vb. I.e the relative positions and locations of the particles.
Constructing the matrix.
From these constraints the C, C,' vectors , J and J' matrices can be computed:
C = [Ci];
C' = [C'i];
Jij = [ dCi/dxj]
J'ij = [dC'i/dxj]
From these we can find Qhat ( from (5),(6) and using a matrix solver to find lamda) and cancel the illegal compents of Q. Then the particle system can be stepped forward one using euler step integration.
Implementation Details
This system was imlemented using arrays and vector/matrix operations. At each time step C,C',J,J' were calculated and Qhat was found. q and qDot could then be updated using euler steps:
To make the nail constraint more stable a scaled identity matrix was added to JWJ^ prior to solving for lamda. This made the sytem solveable. This value added is called diagonal drag in the simulation and can be changed from the menu. Also to increase system stability Qhat was scaled down. This has the affect of making the constraints less rigid, and thus preventing it form exploding (see imlementaion section for details). This value can be changed from the elastisity menu in the program. (Low causes Qhat to be scaled by 1 and thus turns off this feature.)