mathr / blog / #

This toy physics demo featuring inflatable donuts is from April last year, but I didn't get around to blogging about it before. The inline code here is all from the physics tick function, and involves calculating all the forces on each of the vertices of the mesh approximating each torus.

The first step is to find the internal spring forces within each mesh in isolation. For each of the P meshes, we loop over each of its M*N vertices, considering a B*B neighbourhood of vertices on the surface. We accumulate the difference between the vertex and its neighbour vertex in dx, with the length of this vector in ldx. The difference between this length and the target length s->l for this link determines a force, which we accumulate in t->f.

/* intra-mesh spring forces */
#pragma omp parallel for
for (int p = 0; p < P; ++p) {
  for (int m = 0; m < M; ++m) {
    for (int n = 0; n < N; ++n) {
      for (int i = -B/2; i <= B/2; ++i) {
        for (int j = -B/2; j <= B/2; ++j) {
          int mm = (m + i + M) % M;
          int nn = (n + j + N) % N;
          if (i == 0 && j == 0) {
            nn = (n + N / 2) % N;
          }
          R dx[D];
          R ldx = 0;
          for (int d = 0; d < D; ++d) {
            dx[d] = s->x[p][mm][nn][d] - s->x[p][m][n][d];
            ldx += dx[d] * dx[d];
          }
          if (!(ldx > 0)) { continue; }
          ldx = sqrt(ldx);
          R k = stiffness * (s->l[p][m][n][i+B/2][j+B/2] - ldx) / ldx;
          for (int d = 0; d < D; ++d) {
            t->f[p][m][n][d] -= dx[d] * k;
          }
        }
      }
    }
  }
}

To simplify the calculation of the forces between different meshes, we squash the torus mesh down to a loop. We accumulate the center of mass of each segment in t->x and the center of mass of the whole torus in t->x0.

/* mesh center circle */
#pragma omp parallel for
for (int p = 0; p < P; ++p) {
  for (int m = 0; m < M; ++m) {
    for (int n = 0; n < N; ++n) {
      for (int d = 0; d < D; ++d) {
        t->x[p][m][d] += s->x[p][m][n][d] * s->m[p][m][n];
        t->x0[p][d] += s->x[p][m][n][d] * s->m[p][m][n];
      }
      t->m[p][m] += s->m[p][m][n];
      t->m0[p] += s->m[p][m][n];
    }
    for (int d = 0; d < D; ++d) {
      t->x[p][m][d] /= t->m[p][m];
      t->x0[p][d] /= t->m0[p];
    }
  }
}

We also use the previously computed center of mass of segments to model the air pressure inside the torus - this is a constant outwards force on each vertex from the center of its segment.

/* inflate */
#pragma omp parallel for
for (int p = 0; p < P; ++p) {
  for (int m = 0; m < M; ++m) {
    for (int n = 0; n < N; ++n) {
      R dx[D];
      R ldx = 0;
      for (int d = 0; d < D; ++d) {
        dx[d] = s->x[p][m][n][d] - t->x[p][m][d];
        ldx += dx[d] * dx[d];
      }
      ldx = sqrt(ldx);
      for (int d = 0; d < D; ++d) {
        t->f[p][m][n][d] += pressure * dx[d] / ldx;
      }
    }
  }
}

We now check for collisions. We only need to consider nearby meshes, which we test for assuming the meshes are bounded by spheres (s->R,s->r are the initial radii of each torus mesh). If this test passes, we loop over nearby segments. Assuming they are circular, the amount of overlap is scaled by a power law to turn it into a repulsive force between the two segments, stored in t->g.

/* inter mesh forces */
#pragma omp parallel for
for (int p = 0; p < P; ++p) {
  for (int q = p + 1; q < P; ++q) {
    R l = 0;
    for (int d = 0; d < D; ++d) {
      l += (t->x0[p][d] - t->x0[q][d]) * (t->x0[p][d] - t->x0[q][d]);
    }
    l = sqrt(l);
    if (l < s->R[p] + s->r[p] + s->R[q] + s->r[q] + 0.01) {
      for (int mp = 0; mp < M; ++mp) {
        for (int mq = 0; mq < M; ++mq) {
          R dx[D];
          R ldx = 0;
          for (int d = 0; d < D; ++d) {
            dx[d] = t->x[p][mp][d] - t->x[q][mq][d];
            ldx += dx[d] * dx[d];
          }
          if (!(ldx > 0)) { continue; }
          ldx = sqrt(ldx);
          R k = pow(fmax(0.01, ldx - (s->r[p] + s->r[q])), -repulsion) / ldx;
          for (int d = 0; d < D; ++d) {
            R f = dx[d] * k;
            t->g[p][mp][d] += f;
            t->g[q][mq][d] -= f;
          }
        }
      }
    }
  }
}

For smooth lighting we need to calculate the surface normals s->n at each vertex, which we assume point away from the center of the corresponding segment.

/* normals */
#pragma omp parallel for
for (int p = 0; p < P; ++p) {
  for (int m = 0; m < M; ++m) {
    for (int n = 0; n < N; ++n) {
      R l = 0;
      for (int d = 0; d < D; ++d) {
        s->n[p][m][n][d] = s->x[p][m][n][d] - t->x[p][m][d];
        l += s->n[p][m][n][d] * s->n[p][m][n][d];
      }
      l = sqrt(l);
      for (int d = 0; d < D; ++d) {
        s->n[p][m][n][d] /= l;
      }
    }
  }
}

Almost all the forces are now gathered, there just remains gravity and bouyancy which both act in the vertical direction - gravity is a constant, and bouyancy is proportional to depth below the surface of the liquid. For each vertex of each mesh, force f is turned to acceleration a by dividing by mass and accumulated to vertex velocity s->v. Friction with the environment (useful to ensure stability of the simulation) slows down each vertex, and finally the vertex position s->x is updated.

/* accumulate */
#pragma omp parallel for
for (int p = 0; p < P; ++p) {
  for (int m = 0; m < M; ++m) {
    for (int n = 0; n < N; ++n) {
      for (int d = 0; d < D; ++d) {
        R f = t->f[p][m][n][d] + t->g[p][m][d];
        R a = f / s->m[p][m][n];
        if (d == D - 1) {
          a -= gravity + bouyancy * fmin(0, s->x[p][m][n][d]);
        }
        s->v[p][m][n][d] += a * dt;
        s->v[p][m][n][d] *= resistance;
        s->x[p][m][n][d] += s->v[p][m][n][d] * dt;
      }
    }
  }
}

You can download the full rubber dinghy rapids source code. You can compile it like:

gcc -std=c99 -Wall -pedantic -Wextra -march=native -O3 -fopenmp -o dinghy dinghy.c -lglfw -lGLEW -lGL