Cloth simulation

This uses a simple spring-mass system configured in a grid to simulate the internal forces on a piece of flexible material. Combined with gravity this will cause the material to fall and to slip off objects.

The cloth is a 32x32 grid of points, and each point is tied to the eight points around it by springs. Fifteen times a frame, every point is moved a little by gravity and a little by each of its springs. Five spheres bob up and down underneath, and any point that ends up inside a sphere or below the floor is pushed back out.

Setting up

The program is a plain Win32 window with an OpenGL context. The window and GL set-up are from a small Windows OpenGL helper library written by Mark Lord at Gremlin Interactive, and the cloth itself is all in test.cpp, in three functions: Initialise(), GameDraw() and GameMain().

It opens an 800x600 window called "OpenGL Test Program". There's a full-screen version in the same function, and you pick between them with a commenting trick:

//// windowed version
//*
	int width=800;
	int height=600;
	...
/*/
//// fullscreen version
	int width=GetSystemMetrics(SM_CXSCREEN);
	int height=GetSystemMetrics(SM_CYSCREEN);
	...
//*/

As it stands, //* is just a line comment, so the windowed block is compiled. Then /*/ opens a comment that runs down to //*/, which skips the full-screen block. Delete the first slash and it all flips: /* comments out the windowed block, the */ inside /*/ closes it, and the last line is a plain comment again.

The pixel format asks for a double buffer, 16-bit colour and a 16-bit depth buffer. When the window gets its size, WM_SIZE calls GL_Initialise(), which sets up the viewport and a 60 degree perspective, and turns on the depth test:

	// enable back face culling - cull counter clock wise
//	glEnable (GL_CULL_FACE);
	glCullFace (GL_BACK);
	glFrontFace (GL_CW);

	// set up the viewport dimensions
	glViewport (0, 0, width, height);

	// load a default projection matrix at 60 degree F.O.V.
	glMatrixMode (GL_PROJECTION);
	glLoadIdentity ();
	gluPerspective (60.0, (GLfloat)width/(GLfloat)height, 1.0, 30000.0);
	...
	glDisable (GL_BLEND);
  glEnable (GL_DEPTH_TEST);

Back-face culling is set up but its glEnable is commented out, which is just as well here: a cloth has two sides and you see both of them. There's no lighting, no texturing and no blending. Everything is drawn in flat colours, and the depth test sorts out what's in front.

Each frame, GameMain() clears the screen, moves the camera 5 units back, turns the scene by however far the mouse has been dragged, and calls GameDraw(), which does all the work:

void GameMain()
{
	glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
	glLoadIdentity();
	glTranslatef(0.f, 0.f, -5.f);

	glRotatef(-r2*0.2, 1.0, 0.0, 0.0);
	glRotatef(-r1*0.2, 0.0, 1.0, 0.0);

	GameDraw();
	SwapBuffers (gl_dc);
}

The cloth

The cloth starts as a flat, level sheet of 32x32 points, 0.1 apart, at a height of 1.5:

for(j=0; j<32; j++) for(i=0; i<32; i++) {
  x[j][i] = (i-16)*0.1;
  y[j][i] = 1.5;
  z[j][i] = (j-16)*0.1;
}

for(j=0; j<5; j++) for(i=0; i<5; i++) {
  dist[j][i] = sqrt((2-i)*(2-i) + (2-j)*(2-j));
}

The positions are kept in three arrays, x, y and z. They're declared 64x64, but only 32x32 of each is used. With i-16 the sheet runs from -1.6 to 1.5 across and back, so it's half a square off centre.

Every point is joined to the eight points around it. The four along the grid lines rest at 0.1 apart, and the four diagonals at 0.1·√2, about 0.141. The springs aren't stored anywhere. A spring is just "the point next to me", and its rest length is looked up from the dist table by where that neighbour is: dist[2+jj-j][2+ii-i].

The table is 5x5, so it has room for neighbours two squares away, but the loops only ever reach one square out and only the middle 3x3 gets used. Points on the edge of the cloth simply have fewer neighbours: five along a side, three in a corner.

Moving the points

There are no velocities. Each pass works out how far every point should move, from gravity and its springs, and moves it straight there. Here's the inner loop for one point:

mx = 0; my = -0.001; mz = 0;

for(jj=j-1; jj<=j+1; jj++) if(jj>=0 && jj<32) {
  for(ii=i-1; ii<=i+1; ii++) if(ii>=0 && ii<32) if(ii != i || jj != j) {
    vx = x[jj][ii] - x[j][i];
    vy = y[jj][ii] - y[j][i];
    vz = z[jj][ii] - z[j][i];

    d = vx*vx+vy*vy+vz*vz;
    dd = 0.03 * (sqrt(d)-0.1*dist[2+jj-j][2+ii-i])/(d);

    mx += vx*dd;
    my += vy*dd;
    mz += vz*dd;
  }
}

nx[j][i] = x[j][i] + mx; ny[j][i] = y[j][i] + my; nz[j][i] = z[j][i] + mz;

The move, m, starts as gravity: 0.001 straight down. Then for each neighbour, v is the vector to it and d is its length squared, so with len = √d the point moves along v by 0.03 × (len - rest) / len². Since v is len long, that's a distance of 0.03 × (len - rest) / len, towards the neighbour if the spring is stretched and away from it if it's squashed.

Let's say a neighbour along the grid is 0.12 away, so its spring is stretched by 0.02. Then 0.03 × 0.02 / 0.12 = 0.005, and the point moves 0.005 towards it, a quarter of the stretch. The neighbour does the same in its own update, so between them they close half the gap in one pass. Here are a few more lengths for an edge spring:

Hooke's law would move the point in proportion to the stretch. This divides by the length once more, which makes the spring lopsided. However far it's stretched, a spring can never pull a point more than 0.03 in a pass, but a squashed one pushes harder and harder as the two points get close, so neighbours don't collapse onto each other.

The new position goes into nx, ny and nz (new x, not a normal), not back into x, y and z. Every point in a pass is worked out from where its neighbours were at the end of the last pass, and only when all 1,024 are done are the new positions copied back:

for(k=0; k<15; k++) {

for(j=0; j<32; j++) {
  for(i=0; i<32; i++) {
    ...
  }
}

for(j=0; j<32; j++) for(i=0; i<32; i++) {
  x[j][i] = nx[j][i];
  y[j][i] = ny[j][i];
  z[j][i] = nz[j][i];
}

}

That way the order the loop visits the points doesn't matter. If each point were moved in place, the ones late in the loop would see their neighbours' new positions and the early ones wouldn't, and the result would depend on which corner the loop started from.

Doing fifteen passes a frame does two things. A pass only looks one point away, so a tug on one edge of the cloth travels one point per pass. Fifteen passes let it spread fifteen points, half way across the cloth, in a single frame, which makes the cloth much stiffer than one pass would. And gravity is applied fifteen times, so the cloth falls 0.015 a frame.

With no velocities, a falling sheet doesn't speed up. The flat starting sheet has every spring at its rest length, so it just moves down 0.015 every frame, and with nothing in the way it would take 300 frames to fall from 1.5 to the floor at -3.

Why it acts like cloth

Moving each point by the force on it, rather than speeding it up, is like everything moving through thick syrup. Nothing swings or bounces, and a point stops as soon as the forces on it cancel out. But where it stops is the right place, because the forces cancel in the same places they would for a real cloth hanging still.

Take a strip hanging straight down. The bottom point has one spring above it, and it stops when that spring pulls up 0.001 to cancel gravity: 0.03 × (len - 0.1) / len = 0.001, which is a length of about 0.1034. The point above that has to hold up its own 0.001 as well as the 0.001 pulling down from below, so its spring stretches to about 0.1071, and so on up the strip, just as the top of a real rope carries the weight of everything below it.

The springs give it the rest of its character:

• The four edge springs stop it stretching.
• The four diagonals stop it shearing, so the squares can't squash into thin diamonds.
• Nothing stops it bending. There are no springs to the points two away, so it folds along any line of points with no resistance at all, which is what makes it cloth rather than a sheet of metal.

Springs to the points two away would give it some stiffness against bending, and the 5x5 table has room for their rest lengths, but the loops stop at one.

There's also nothing to stop the cloth passing through itself. Once it has slid off the spheres it settles as a crumpled heap on the floor rather than lying flat.

Collisions

Once a point's new position is worked out, it's checked against the floor and then against each sphere:

//  Floor collision

if(ny[j][i] < -2.999) ny[j][i] = -2.999;

//  Sphere collision

for(kk = 0; kk<spheres; kk++) {
  d = (nx[j][i]-sx[kk])*(nx[j][i]-sx[kk]) + (ny[j][i]-sy[kk])*(ny[j][i]-sy[kk]) + (nz[j][i]-sz[kk])*(nz[j][i]-sz[kk]);
  if(d < sr[kk]*sr[kk]) {
    d = sr[kk]/sqrt(d);
    nx[j][i] = (nx[j][i]-sx[kk])*d + sx[kk];
    ny[j][i] = (ny[j][i]-sy[kk])*d + sy[kk];
    nz[j][i] = (nz[j][i]-sz[kk])*d + sz[kk];
  }
}

The floor is just a height. Anything below -2.999 is put back at -2.999, a thousandth above the floor squares drawn at -3, so the cloth is drawn on top of them rather than mixed in. As far as the physics goes the floor carries on for ever in every direction, though only a 4x4 board is drawn.

For a sphere, d is the squared distance from the centre. If that's less than the radius squared, the point is inside, and it's moved straight out along the line from the centre until it's exactly on the surface. sr[kk]/sqrt(d) is how much the offset from the centre needs scaling by: a point 0.4 from the centre of a sphere of radius 0.5 gets its offset scaled by 0.5 / 0.4 = 1.25, which puts it at 0.5.

That's all the collision there is, and it's enough to make the cloth slip. Let's follow a point resting on the side of a sphere:

Gravity moves it 0.001 straight down, a little way into the sphere. The push moves it back out along the radius. The two don't cancel: what's left over is a move along the surface, downhill. At an angle θ from the top of the sphere it's about 0.001 × sin θ a pass, nothing at the very top and the full 0.001 on the side.

There's no friction, as nothing ever holds a point in place, so the only thing keeping the cloth on a sphere is the rest of the cloth. Cloth over the top of a sphere stays there while the pull on each side balances, and slides off as soon as one side is heavier.

The spheres are drawn with a radius of 0.49 while the collision uses 0.5, so the cloth sits just outside the surface rather than cutting into it. And because the check uses where the spheres are this frame, a sphere coming up underneath the cloth pushes up every point it reaches and lifts the cloth with it.

The moving spheres

The five spheres all have a radius of 0.5 and only ever move up and down. Each frame they start again from their rest positions, the counter a goes up by one, and each one's height is set from a wave:

for(i=0; i<5; i++) {
  sx[i] = ssx[i];
  sy[i] = ssy[i];
  sz[i] = ssz[i];
  sr[i] = ssr[i];
}

a ++;

sy[0] = ssy[0] - 1.16*cos(-a*3.14159*0.0007722-0.5947);
sy[1] = ssy[1] - 1.67*cos(-a*3.14159*0.0006224+0.1344);
sy[2] = ssy[2] - 1.42*sin(-a*3.14159*0.0004629+0.6143);
sy[3] = ssy[3] + 1.55*cos(-a*3.14159*0.0003726-0.8941);
sy[4] = ssy[4] - 1.71*cos(-a*3.14159*0.0001722-0.7228);

Each wave has its own size, speed and starting point. The angle goes up by π times the rate each frame, so a full cycle takes 2 / rate frames:

None of the rates is a simple multiple of another, so the spheres keep drifting in and out of step and the scene doesn't repeat for a long time. Sphere 3 goes down to -3.25, three quarters of the way into the floor, and sphere 1 comes up to 1.27, which puts its top above the cloth's starting height.

Several of them are close enough sideways to overlap, and as there's no collision between the spheres, they simply pass through each other.

Drawing

Everything is drawn as quads, each one a four-vertex GL_TRIANGLE_STRIP with a single glColor3f, in a checkerboard of two colours picked by ((i^j)&1). With no lighting, a quad is one flat colour whichever way it faces, so it's the checks that show the shape: you can see the cloth fold because its squares get squashed and turned.

Spheres

A sphere is built from a 16x16 grid of points. i steps down from the top in sixteenths of π, and j goes round in eighths of π:

for(j=0; j<16; j++) for(i=0; i<16; i++) {
  x = sr*sin(i*3.14159/16);
  y = sr*cos(i*3.14159/16);
  z = 0;

  //  Rotate about y axis

  double nx = z*sin(j*3.14159/8) + x*cos(j*3.14159/8);
  double ny = y;
  double nz = z*cos(j*3.14159/8) - x*sin(j*3.14159/8);

  sphere_x[j][i] = sx+nx;
  sphere_y[j][i] = sy+ny;
  sphere_z[j][i] = sz+nz;
}

...

  ii = (i+1)&15;
  jj = (j+1)&15;

  ...

  glBegin(GL_TRIANGLE_STRIP);
    glColor3f(r, g, b);
    glVertex3f(sphere_x[ j][ i], sphere_y[ j][ i], sphere_z[ j][ i]);
    glVertex3f(sphere_x[ j][ii], sphere_y[ j][ii], sphere_z[ j][ii]);
    glVertex3f(sphere_x[jj][ i], sphere_y[jj][ i], sphere_z[jj][ i]);
    glVertex3f(sphere_x[jj][ii], sphere_y[jj][ii], sphere_z[jj][ii]);
  glEnd();

The next column, (j+1)&15, wraps from 15 back to 0 to close the circle, which is right. The next row, (i+1)&15, wraps the same way, which isn't. The rows only reach 15π/16, just short of the bottom pole, so the last band of quads joins the bottom ring back up to row 0, the top pole, straight through the middle of the sphere:

From the outside you never see it, as the cone is inside the sphere and the depth test hides it. The only place it shows is underneath, where the small gap left for the bottom cap is filled by the inside of the cone. Row 0 is a single point too, so the quads round the top are really triangles.

The spheres are red and white, and drawn at the radius less 0.01, as covered in Collisions.

The floor

The floor is an 8x8 board of squares 0.5 across at a height of -3, running from -2 to 2 each way, in light and dark grey (0.9 and 0.1). It never moves, but it's sent to OpenGL again every frame like everything else.

The cloth

The cloth is 31x31 quads between its 32x32 points, in yellow and blue:

for(j=0; j<31; j++) {
  for(i=0; i<31; i++) {

    if(((i^j)&1) == 0) {
      r = 1.0; g = 1.0; b = 0.0;
    } else {
      r = 0.0; g = 0.0; b = 1.0;
    }

    glBegin(GL_TRIANGLE_STRIP);
      glColor3f(r, g, b); glVertex3f(x[j  ][i  ], y[j  ][i  ], z[j  ][i  ]);
      glColor3f(r, g, b); glVertex3f(x[j  ][i+1], y[j  ][i+1], z[j  ][i+1]);
      glColor3f(r, g, b); glVertex3f(x[j+1][i  ], y[j+1][i  ], z[j+1][i  ]);
      glColor3f(r, g, b); glVertex3f(x[j+1][i+1], y[j+1][i+1], z[j+1][i+1]);
    glEnd();

  }
}

A strip of four vertices makes two triangles, the first three vertices and the last three, so when the four corners of a square aren't flat it bends along the diagonal from [j][i+1] to [j+1][i]. The cloth is drawn straight from the simulation's own arrays, with nothing smoothed or added.

Each frame

WinMain() calls Initialise() once and then just loops: handle any waiting Windows messages, then call GameMain(). There's no timer. Each time round, the spheres move one step and the cloth gets its fifteen passes, so it all runs as fast as the machine can draw it, and a faster machine runs it faster.

The physics and the drawing share one function, GameDraw(), and the window's WM_PAINT handler calls it too, so a repaint also moves everything on a step, without clearing the screen or placing the camera first.

Controls

Drag with the left mouse button to turn the view, and press Escape to quit. While the button is down, every mouse move adds the distance moved to r1 and r2:

case WM_MOUSEMOVE:
  POINT pt;
  GetCursorPos(&pt);

  if(lbutton) {
    r1 += (mx-pt.x);
    r2 += (my-pt.y);
  }

  mx = pt.x;
  my = pt.y;

  return 0;

GameMain() turns the scene by 0.2 degrees per pixel: r2 tilts it about the x axis and r1 spins it about the y axis. The last position is recorded on every move, button or not, so the view doesn't jump when you press it. The button press calls SetCapture(), so the moves keep coming even when the mouse leaves the window, and the position comes from GetCursorPos(), in screen pixels.

Escape arrives as a WM_CHAR and posts the window a close message, which destroys it and ends the loop.

The WebGL version

The demo at the top is a WebGL port. The physics in its demo.js is the C++ line for line, in double precision with Float64Array, so the cloth falls and slides just as the original does.

The drawing needs something different, as WebGL has no fixed-function pipeline. This program only ever uses glColor3f and the depth test, though, so a shader that passes the colour straight through is all it takes:

attribute vec3 aPos;
attribute vec3 aCol;
uniform mat4 uMVP;
varying vec3 vCol;
void main() {
  vCol = aCol;
  gl_Position = uMVP * vec4(aPos, 1.0);
}

precision mediump float;
varying vec3 vCol;
void main() {
  gl_FragColor = vec4(vCol, 1.0);
}

Each four-vertex strip becomes the same two triangles, so the cloth's squares bend along the same diagonal, and the spheres, floor and cloth all go to the GPU in one buffer with one draw call. The camera is the same gluPerspective, translate and two rotations, built as matrices in JavaScript.

A few things can't be the same:

• Escape can't close a browser tab, so it stops the program instead, and a click starts it again from Initialise(), like running the program again.
• The original asked for 16-bit colour and a 16-bit depth buffer. The browser gives 8 bits per colour channel and usually a finer depth buffer, so the port doesn't reproduce the coarser 16-bit output.
• The original has no clock and runs a step per frame as fast as it can. The port runs 60 steps a second whatever the screen's refresh rate, and catches up by up to 15 steps if a frame is late. At that rate the sphere cycles in the table take from 43 seconds for sphere 0 to over three minutes for sphere 4. As in the original, the first frame you see comes after two steps, because the window's first WM_PAINT runs an extra one.
• The canvas keeps the shape of the original window's 792x573 drawing area, and mouse moves are scaled so that a drag across the canvas turns the view as far as a drag across that window did.