## The Basics

The wave equation is a classic equation in partial differential equations (PDEs). It has a number of uses in physics and is often used to explain energy transfer through space. It is an example of a time-dependent equation, which simply means that the “solution”, or state of the system, changes over time. The PDE explains how the solution takes place. In the wave equation that we’ll simulate, we will use the notion of the Laplacian, which is a linear operator. The practical implication of the “linearity” of this problem is that there are good analytical insights into the problem and how it behaves, fast ways to solve it, and many different types of ways to go about solving it.

In this article, I’ll assume the reader has a decent understanding of calculus and linear algebra, but I’ll also try to explain things in as intuitive a way as I can..

The classical wave equation has the form

$$\frac{\partial^2}{\partial t^2}u = c \nabla^2 u$$

where $u$ is the “solution” which possibly changes value as we move around in space, and as we vary time. This means that $u$ is more properly written as a function of space and time, like $u(x,t)$, but for brevity, we’ll drop the extra notation. In order to simulate waves, we will need a way to describe the Laplacian $\nabla^2$, and prescribe a way to advance our current solution at a point in time to the next point in time, so we can simulate and animate this process.

## Waves on a Mesh

We can find descriptions of wave simulations all over the web, so why is this article any different? Most of the articles on the web assume that space is “flat”. In 2D, the term “flat” just refers to simulating the wave as something that moves on a flat plane, rather than on rolling hills for example. What this usually means programming wise, is that there is a regular (equally-spaced) grid of points and we calculate how the value of $u$ changes over time at each point in the grid in order to generate waves. While this is perfectly valid, it doesn’t allow us to do things like simulate waves on the surface of a sphere, which is obviously curved everywhere! We’ll see that there is actually a very nice method of simulating waves on a triangular mesh (surface) based on the idea of “flux”.

### The Laplacian, Divergence Theorem on a Triangular Mesh, and Flux

The Laplacian is a fundamental linear operator in PDE Theory and physics. It has deep connections to the notion of curvature. It actually has a very simple meaning though: It is simply the sum of all the second (non-mixed) partial derivatives. A slightly more rich definition, which probably sounds more confusing but is actually the one we need in this article, is that it is the divergence of the gradient. One simple (“layperson”) way to think about the Laplacian is the following: If the value is constant around a point, then the Laplacian there will be zero. This is like a flat plane, no curvature means zero Laplacian. If the value is change, but the rate at which it’s changing is a constant, then this will also be zero. This is like a hill with no valleys or crests. However, if there is a crest/valley, then there will be a non-zero Laplacian value there. It’s also worth noting that this value is signed depending on if the “bump” is locally pushing outwards or inwards relative to neighboring values.

Anyway… We mentioned that the Laplacian has a definition as the divergence of the gradient:

$$ \nabla^2 u \equiv \nabla \cdot (\nabla u)$$

which means that our wave equation must also be equivalent to

$$ \frac{\partial^2}{\partial t^2} u = \nabla \cdot (\nabla u) $$

The divergence theorem tells us that if we take a region of space and integrate the divergence of something in that region, it will be equal to the net flow across the boundary of that region:

$$ \int_\Omega \nabla \cdot \mathbf{g} = \int_{\partial\Omega} \mathbf{n} \cdot \mathbf{g} $$

Another way to say this is that if we imagine the net effect of adding up all the “compression” and “expansion” occurring within a space, that will be equal to adding up how much stuff crosses the boundary of that space. Practically though, this will lend us good insight into a way to solve the problem.

In our problem, we will assume that we have a value of $u$ known on each triangle in a manifold mesh. If there’s a question about what manifold means here, you can either google, ask in a comment, or basically imagine almost any model you would make by subdividing a sphere in Blender/Maya/etc and pulling things around. One fun thing about this method is that the mesh should be topologically manifold, but the immersion doesn’t have to be. What that means is we can do strange things like make surfaces that are manifold but don’t have any way to exist in reality, like a Klein Bottle, and actually simulate waves on them.

Since we have a value of $u$ on each triangle, we will want to know how to change that value over time, but in order to do that, we will want to know how to compute the Laplacian at each triangle. We will utilize our last result from the Divergence Theorem, which told us that the integral of the Laplacian was equal to adding up the gradient across the boundary. Integrating the entire wave equation over a single face $\Omega_i$

$$ \frac{\partial^2}{\partial t^2} \int_{\Omega_i} u_i = \frac{\partial^2}{\partial t^2} |\Omega_i| u_i = \int_{\partial\Omega_i} \mathbf{n} \cdot \nabla \mathbf{u} $$

$$\frac{\partial^2}{\partial t^2} u_i = |\Omega_i|^{-1} \int_{\partial\Omega_i} \mathbf{n} \cdot \nabla \mathbf{u}$$

To illustrate the “normals” for a the boundary of a face:

So now we need to calculate the gradient of the solution $u$ at the “boundary” of the triangle. The boundary of the triangle is just formed by the edges!. One simple estimator for the gradient across an edge is to calculate the difference between the value at the triangle across the edge $u_j$ and our value $u_i$ and to divide by the shortest distance between them.

Calculating this distance seems simple enough, right? Just get the centers of the triangles $C_i$, $C_j$ and compute $|C_i – C_j|$, right? No. The “straight line” distance is not correct, since our notion of distance now depends on the path lying on the surface. That means that if we’re given two triangles that share an edge, we need a path that goes from one triangle center to the other and that path must stay on the triangles. It’s “obvious” that this path must cross the edge at some point along the path. On each triangle, a portion of the path must go from the center to this special point on the edge. Since the triangles are themselves flat though, we know that this portion of the path will be a straight line. This results in the following: The minimum distance between the two centers is made by drawing a straight line from one center to a special point on the edge, and then continuing onward from that point on the edge to the center of the second triangle!

So this is actually pretty simple, but we now need to know what we should choose for center points, and also what the point on the edge, $e^*$ is that makes the total path length $|C_i-e^*|+|C_j-e^*|$ as small as possible. The proof is simple, but if our triangles have a property known as being “well-centered”, then they will have a special center called the circumcenter which lies inside the triangle. If all of the triangles have this property, then we can use them as the center points $C_i$, $C_j$, and the edge point to choose for the middle of the path is simply the middle of the edge. This is pretty simple to show with basic trig/geometry. Later perhaps. Side note: Another way to derive this would be to fix one of the triangles, compute a rotation that rotates the neighbouring triangle into the same plane, compute the straight line between them, and then find where this intersects the edge joining the two.

This is in fact all we need to compute the integral along the boundary of a triangle! For each triangle, we will compute the Laplacian there as the sum of gradients across the edges of the triangle. A natural result of this is that if triangles along the boundary of the surface (if there is a boundary at all) will have less than three edges contributing to the value of the Laplacian.

So, procedurally, we will need to precompute a few things in order to calculate the Laplacian: the circumcenters, the area of triangles, edge midpoints, etc. which can all be calculated pretty easily. From then on, the Laplacian at a triangle is computed as the sum of gradients across the edges:

For each Triangle Ti: L[i] = 0 # Laplacian at Ti For each Triangle Tj neighboring Ti: Let eij be the edge shared by Ti and Tj, having vertices a and b Let Ci be the circumcenter of Ti, Cj the circumcenter of Tj Let Ai be the area of Ti # midpoint of eij E = .5 * (eij.a + eij.b) # Compute the length of the path joining the centers D = |Ci - E| + |Cj + E| # Compute the flux across this edge F = (uj - ui) * |eij| / D L[i] += F / Ai

While this will work, it is a little slow… Neighboring triangles will compute the same value twice, just with opposite signs, so in practice they can be computed simultaneously, i.e. one per pair of neighbouring faces. Also, I found it much nicer to encode this into a matrix: Since $u_j$ and $u_i$ are the only unknowns at runtime and they only appear linearly in the expression for the Laplacian (this is expected since the Laplacian is linear!), we can write the Laplacian as a matrix $\mathbf{L}$, so the laplacian at each triangle is computed by just multiplying $\mathbf{Lu}$. The source code will make this more obvious! I have created a very small but working sparse matrix class for doing these calculations.

By the way, if face areas are neglected while you’re forming the Laplacian, the result may still “work”, but you’ll notice funky behaviour when a wave front crosses a triangle of a different size from it’s surrounding neighbours!

## Obligatory Demo

A demo of this idea is viewable at GitHub! I highly recommend you run this code in Chrome for all of the typical reasons.

If you have any questions or comments, please do write below!

Cool. For more updated reasons it also works very well in FF 17.0.1. at least.

Cool! It would be interesting to see some other equations on the mesh. Maybe you could try seeing if some of Tim Hutton’s work on reaction diffusion equations works well on surfaces:

https://www.youtube.com/watch?v=QJB1Jsk_oTE&list=UUVtErl9kecaaeeicumm7jXg&index=1

Another interesting thing could be to try computing the spectrum of the Laplacian operator, since it has many useful applications. Of course this isn’t anything new, but to my knowledge I’ve never seen a good WebGL demo that does it.

The reaction-diffusion video strongly reminds me of smooth life gliders. I’ll have to look into their work to see exactly what they’re doing, though I have done advection-reaction-diffusion work before.

I am currently writing a couple of articles on Octrees and simulating fluids with Vortons.

It would be neat to visualize the eigenspectrum of the Laplacian on an arbitrary mesh