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