Fun Equations of the Torus (FEniCS part 2)

The story so far: we wished to create a horseshoe-like initial condition on the torus. We needed to (1) come up with the equation for the bounding curve, and (2) distinguish between the inside and outside. One tool is perfect for the job: level sets (another common name for them is isosurfaces, and some more dimension-specific terms like level contours or level surfaces). Surfaces (or curves) of the form f(x,y,z) = c (or f(x,y) = c). In order to do this, we first start with the classic parametrization of the torus as a surface of revolution:

x = (R+r cos(u))cos(θ)
y = (R+r cos(u))sin(θ)
z = r sin(u),

with 0 ≤ u ≤ 2π and 0 ≤ θ ≤ 2π. What this is saying is that we take a circle of radius r in a vertical plane (its angle coordinate is u) and bring it out so its center is at radius R from the origin and still lying in the xy-plane (that's why the R only affects the term with r cos(u)). Then we revolve it around the z-axis, which is where the cos(θ) and sin(θ) come in. We need to come up with an equation in terms of x, y, and z.

If we square these three coordinates and add, we get

x² + y² + z² = (R+r cos(u))² + r²sin(u)² = R²+2Rr cos(u) + r².

However, r cos(u) is simply the extra distance one has to go from the radius R to the coordinates (x,y,0), that is, r cos(u) =  ρ-R, where ρ² = x² + y² is the usual polar coordinate (click to enlarge):

This gives us

x² + y² + z² = (R+r cos(u))² + r²sin(u)² = R²+2R(ρ-R)+ r² = r²+ 2Rρ – R².

Moving everything on the other side, except 2Rρ, and dividing by 2R, we finally get

(x² + y² + z² + R² – r²)/(2R) = ρ,

or

(x² + y² + z² + R² – r²)²/(4R²)= x² + y², which finally gives us

f(x,y,z) = (x² + y² + z² + R² – r²)²/(4R²) – x² – y² = 0.

This is the Cartesian equation of the torus (really, a Cartesian equation, because there are several). Of course to be really proper, we should include dependence on R and r in there too (by including that in Mathematica, I got to see how the surfaces change as I vary each parameter). The torus plus its solid interior (a donut!) is given by replacing it with a ≤ rather than a = (note also we can determine which side is inside and which side is outside by computing the differential df, also variously known as the gradient vector–it points in the direction of increasing f). To get nearby shapes, we can consider other level sets f(x,y,z) = c for c near to 0. The picture introducing this post is precisely this, for c in the interval [-1,1] with increments of 1/4. Notice that the largest two of these has no hole (look carefully at the upper center part of the picture). So it's almost nested tori.

The interesting thing about using a level set description of a family of surfaces is that nothing special happens when some topological transition (shape with a hole becoming a shape without a hole) occurs. Well, almost nothing special—the transition point is precisely when c is not a regular value of f, that is, there are one or more points on the level set for which the gradient df vanishes at p. Such ps are called critical points (and is the starting point of Morse theory). Making use of level sets and tracking various geometric dynamics (surfaces, interfaces, etc) using these concepts in numerical methods is the level set method. The trickiness, however, is finding the right level set function to use, and it can be more of an art than a science (although I don't necessarily think the strict separation of art and science is that great of a thing. But that's a story for another day).

Back to our original problem, however. We wanted to see what the cross section looks like and how it changes as we slice a bagel at different angles. This is quite simple to do in terms of level sets (which is why we went through all the trouble to derive it). We can imagine the slicing plane to be z = mx + b. The slope m is the tangent of the angle α between the slicing plane and the horizontal xy-plane. And, of course, we wish it to lie on the surface of the torus, so that we have a system of equations

z = mx + b,
f(x,y,z) = (x² + y² + z² + R² – r²)²/(4R²) – x² – y² = 0.

The key is, then, to substitute one into the other, getting

g(x,y) =  (x² + y² + (mx + b)² + R² – r²)²/(4R²) – x² – y² = 0.

If we consider this as a function of x and y only, that is, a subset of the plane, it gives the xy-projection of the 3D slice. To get the interior of the cross section, we use a ≤ instead of a = as before. The values chosen here are r = 1, R = 1.61803398874 (the Golden Ratio), m = 1/2 and b = 1/3. And that gives us our inside vs. the outside of the domain we used in the heat equation. If we had wanted to make it more circular, we could eliminate the projection by rotating the entire picture by the angle α. This makes the equations much more complicated, however.