A Quick and Dirty Introduction to Exterior Calculus — Part IV: Differential Operators

Originally we set out to develop exterior calculus. The objects we’ve looked at so far — \(k\)-forms, the wedge product \(\wedge\) and the Hodge star \(\star\) — actually describe a more general structure called an exterior algebra. To turn our algebra into a calculus, we also need to know how quantities change, as well as how to measure quantities. In other words, we need some tools for differentiation and integration. Let’s start with differentiation.

In our discussion of surfaces we briefly looked at the differential \(df\) of a surface \(f: M \rightarrow \mathbb{R}^3\), which tells us something about the way tangent vectors get “stretched out” as we move from the domain \(M\) to a curved surface sitting in \(\mathbb{R}^3\). More generally \(d\) is called the exterior derivative and is responsible for building up many of the differential operators in exterior calculus. The basic idea is that \(d\) tells us how quickly a \(k\)-form changes along every possible direction. But how exactly is it defined? So far we’ve seen only a high-level geometric description.

Div, Grad, and Curl

Before jumping into the exterior derivative, it’s worth reviewing what the basic vector derivatives \(\mathrm{div}\), \(\mathrm{grad}\), and \(\mathrm{curl}\) do, and more importantly, what they look like. The key player here is the operator \(\nabla\) (pronounced “nabla”) which can be expressed in coordinates as the vector of all partial derivatives:

\[ \nabla := \left( \frac{\partial}{\partial x^1}, \ldots, \frac{\partial}{\partial x^n} \right). \]

For instance, applying \(\nabla\) to a scalar function \(\phi: \mathbb{R}^n \rightarrow \mathbb{R}\) yields the gradient

\[ \nabla\phi = \left( \frac{\partial f}{\partial x^1}, \ldots, \frac{\partial f}{\partial x^n} \right), \]

which can be visualized as the direction of steepest ascent on some terrain:

We can also apply \(\nabla\) to a vector field \(X\) in two different ways. The dot product gives us the divergence

\[ \nabla \cdot X = \frac{\partial X^1}{\partial x^1} + \cdots + \frac{\partial X^n}{\partial x^n} \]

which measures how quickly the vector field is “spreading out”, and on \(\mathbb{R}^3\) the cross product gives us the curl

\[ \nabla \times X = \left( \frac{\partial X_3}{\partial x^2} - \frac{\partial X_2}{\partial x^3}, \frac{\partial X_1}{\partial x^3} - \frac{\partial X_3}{\partial x^1}, \frac{\partial X_2}{\partial x^1} - \frac{\partial X_1}{\partial x^2} \right), \]

which indicates how much a vector field is “spinning around.” For instance, here’s a pair of vector fields with a lot of divergence and a lot of curl, respectively:

(Note that in this case one field is just a 90-degree rotation of the other!) On a typical day it’s a lot more useful to think of \(\mathrm{div}\), \(\mathrm{grad}\) and \(\mathrm{curl}\) in terms of these kinds of pictures rather than the ugly expressions above.

Think Differential

Not surprisingly, we can express similar notions using exterior calculus. However, these notions will be a bit easier to generalize (for instance, what does “curl” mean for a vector field in \(\mathbb{R}^4\), where no cross product is defined?). Let’s first take a look at the exterior derivative of 0-forms (i.e., functions), which is often just called the differential. To keep things simple, we’ll start with real-valued functions \(\phi: \mathbb{R}^n \rightarrow \mathbb{R}\). In coordinates, the differential is defined as

\[ d\phi := \frac{\partial \phi}{\partial x^1} dx^1 + \cdots + \frac{\partial \phi}{\partial x^n} dx^n. \]

It’s important to note that the terms \(\frac{\partial \phi}{\partial x^i}\) actually correspond to partial derivatives of our function \(\phi\), whereas the terms \(dx^i\) simply denote an orthonormal basis for \(\mathbb{R}^n\). In other words, you can think of \(d\phi\) as just a list of all the partial derivatives of \(\phi\). Of course, this object looks a lot like the gradient \(\nabla \phi\) we saw just a moment ago. And indeed the two are closely related, except for the fact that \(\nabla \phi\) is a vector field and \(d\phi\) is a 1-form. More precisely,

\[ \nabla \phi = (d\phi)^\sharp. \]

Directional Derivatives

Another way to investigate the behavior of the exterior derivative is to see what happens when we stick a vector \(u\) into the 1-form \(df\). In coordinates we get something that looks like a dot product between the gradient of \(f\) and the vector \(u\):

\[ df(u) = \frac{\partial f}{\partial x^1} u^1 + \cdots + \frac{\partial f}{\partial x^n} u^n. \]

For instance, in \(\mathbb{R}^2\) we could stick in the unit vector \(u = (1,0)\) to get the partial derivative \(\frac{\partial f}{\partial x^1}\) along the first coordinate axis:

(Compare this picture to the picture of the gradient we saw above.) In general, \(df(u)\) represents the directional derivative of \(f\) along the direction \(u\). In other words, it tells us how quickly \(f\) changes if we take a short walk in the direction \(u\). Returning again to vector calculus notation, we have

\[ df(u) = u \cdot \nabla f. \]

Properties of the Exterior Derivative

How do derivatives of arbitrary \(k\)-forms behave? For one thing, we expect \(d\) to be linear — after all, a derivative is just the limit of a difference, and differences are certainly linear! What about the derivative of a wedge of two forms? Harkening back to good old-fashioned calculus, here’s a picture that explains the typical product rule \(\frac{\partial}{\partial x}(f(x)g(x)) = f’(x)g(x) + f(x)g’(x) \):

The dark region represents the value of \(fg\) at \(x\); the light blue region represents the change in this value as we move \(x\) some small distance \(h\). As \(h\) gets smaller and smaller, the contribution of the upper-right quadrant becomes negligible and we can write the derivative as the change in \(f\) times \(g\) plus the change in \(g\) times \(f\). (Can you make this argument more rigorous?) Since a \(k\)-form also measures a (signed) volume, this intuition also carries over to the exterior derivative of a wedge product. In particular, if \(\alpha\) is a \(k\)-form then \(d\) obeys the rule

\[ d(\alpha \wedge \beta) = d\alpha \wedge \beta + (-1)^{k}\alpha \wedge d\beta. \]

which says that the rate of change of the overall volume can be expressed in terms of changes in the constituent volumes, exactly as in the picture above.

Exterior Derivative of 1-Forms

To be a little more concrete, let’s see what happens when we differentiate a 1-form on \(\mathbb{R}^3\). Working things out in coordinates turns out to be a total mess, but in the end you may be pleasantly surprised with the simplicity of the outcome! (Later on we’ll see that these ideas can also be expressed quite nicely without coordinates using Stokes’ theorem, which paves the way to differentiation in the discrete setting.) Applying the linearity of \(d\), we have

&=& d(\alpha_1 dx^1 + \alpha_2 dx^2 + \alpha_3 dx^3) \\
&=& d(\alpha_1 dx^1) + d(\alpha_2 dx^2) + d(\alpha_3 dx^3).

Each term \(\alpha_j dx^j\) can really be thought of a wedge product \(\alpha_j \wedge dx^j\) between a 0-form \(\alpha_j\) and the corresponding basis 1-form \(dx^j\). Applying the exterior derivative to one of these terms we get

\[ d(\alpha_j \wedge dx^j) = (d\alpha_j) \wedge dx^j + \alpha_j \wedge \underbrace{(ddx^j)}_{=0} = \frac{\partial \alpha_j}{\partial x^i} dx^i \wedge dx^j. \]

To keep things short we used the Einstein summation convention here, but let’s really write out all the terms:

d\alpha &=& \frac{\partial \alpha_1}{\partial x^1} dx^1 \wedge dx^1 &+& \frac{\partial \alpha_1}{\partial x^2} dx^2 \wedge dx^1 &+& \frac{\partial \alpha_1}{\partial x^3} dx^3 \wedge dx^1 \\
&& \frac{\partial \alpha_2}{\partial x^1} dx^1 \wedge dx^2 &+& \frac{\partial \alpha_2}{\partial x^2} dx^2 \wedge dx^2 &+& \frac{\partial \alpha_2}{\partial x^3} dx^3 \wedge dx^2 \\
&& \frac{\partial \alpha_3}{\partial x^1} dx^1 \wedge dx^3 &+& \frac{\partial \alpha_3}{\partial x^2} dx^2 \wedge dx^3 &+& \frac{\partial \alpha_3}{\partial x^3} dx^3 \wedge dx^3. \\

Using the fact that \(\alpha \wedge \beta = -\beta \wedge \alpha\), we get a much simpler expression

d\alpha &=& ( \frac{\partial \alpha_3}{\partial x^2} - \frac{\partial \alpha_2}{\partial x^3} ) dx^2 \wedge dx^3 \\
&& ( \frac{\partial \alpha_1}{\partial x^3} - \frac{\partial \alpha_3}{\partial x^1} ) dx^3 \wedge dx^1 \\
&& ( \frac{\partial \alpha_2}{\partial x^1} - \frac{\partial \alpha_1}{\partial x^2} ) dx^1 \wedge dx^2. \\

Does this expression look familiar? If you look again at our review of vector derivatives, you’ll recognize that \(d\alpha\) basically looks like the curl of \(\alpha^\sharp\), except that it’s expressed as a 2-form instead of a vector field. Also remember (from our discussion of Hodge duality) that a 2-form and a 1-form are not so different here — geometrically they both specify some direction in \(\mathbb{R}^3\). Therefore, we can express the curl of any vector field \(X\) as

\[ \nabla \times X = \left( \star d X^\flat \right)^\sharp. \]

It’s worth stepping through the sequence of operations here to check that everything makes sense: \(\flat\) converts the vector field \(X\) into a 1-form \(X^\flat\); \(d\) computes something that looks like the curl, but expressed as a 2-form \(dX^\flat\); \(\star\) turns this 2-form into a 1-form \(\star d X^\flat\); and finally \(\sharp\) converts this 1-form back into the vector field \(\left( \star d X^\flat \right)^\sharp\). The take-home message here, though, is that the exterior derivative of a 1-form looks like the curl of a vector field.

So far we know how to express the gradient and the curl using \(d\). What about our other favorite vector derivative, the divergence? Instead of grinding through another tedious derivation, let’s make a simple geometric observation: in \(\mathbb{R}^2\) at least, we can determine the divergence of a vector field by rotating it by 90 degrees and computing its curl (consider the example we saw earlier). Moreover, in \(\mathbb{R}^2\) the Hodge star \(\star\) represents a rotation by 90 degrees, since it identifies any line with the direction orthogonal to that line:

Therefore, we might suspect that divergence can be computed by first applying the Hodge star, then applying the exterior derivative:

\[ \nabla \cdot X = \star d \star X^\flat. \]

The leftmost Hodge star accounts for the fact that \(d \star X^\flat\) is an \(n\)-form instead of a 0-form — in vector calculus divergence is viewed as a scalar quantity. Does this definition really work? Let’s give it a try in coordinates on \(\mathbb{R}^3\). First, we have

\begin{array}{rcl} \star X^\flat &=& \star( X_1 dx^1 + X_2 dx^2 + X_3 dx^3 ) \\
&=& X_1 dx^2 \wedge dx^3 + X_2 dx^3 \wedge dx^1 + X_3 dx^1 \wedge dx^2.

Differentiating we get

\begin{array}{rcl} d \star X^\flat &=& \frac{\partial X_1}{\partial x^1} dx^1 \wedge dx^2 \wedge dx^3 + \\
&& \frac{\partial X_2}{\partial x^2} dx^2 \wedge dx^3 \wedge dx^1 + \\
&& \frac{\partial X_3}{\partial x^3} dx^3 \wedge dx^1 \wedge dx^2, \\

but of course we can rearrange these wedge products to simply

\[ d \star X^\flat = \left( \frac{\partial X_1}{\partial x^1} + \frac{\partial X_2}{\partial x^2} + \frac{\partial X_3}{\partial x^3} \right) dx^1 \wedge dx^2 \wedge dx^3. \]

A final application of the Hodge star gives us the divergence

\[ \star d \star X^\flat = \frac{\partial X_1}{\partial x^1} + \frac{\partial X_2}{\partial x^2} + \frac{\partial X_3}{\partial x^3} \]

as desired.

In summary, for any scalar field \(\phi\) and vector field \(X\) we have

\nabla \phi &=& (d\phi)^\sharp \\
\nabla \times X &=& \left( \star d X^\flat \right)^\sharp \\
\nabla \cdot X &=& \star d \star X^\flat

One cute thing to notice here is that (in \(\mathbb{R}^3\)) \(\mathrm{grad}\), \(\mathrm{curl}\), and \(\mathrm{div}\) are more or less just \(d\) applied to a \(0-\), \(1-\) and \(2-\) form, respectively.

The Laplacian

Another key differential operator from vector calculus is the scalar Laplacian which (confusingly!) is often denoted by \(\Delta\) or \(\nabla^2\), and is defined as

\[ \Delta := \nabla \cdot \nabla, \]

i.e., the divergence of the gradient. Although the Laplacian may seem like just yet another in a long list of derivatives, it deserves your utmost respect: the Laplacian is central to the most fundamental of physical laws (any diffusion process and all forms of wave propagation, including the Schrödinger equation); its eigenvalues capture almost everything there is to know about a given piece of geometry (can you hear the shape of a drum?). Heavy tomes and entire lives have been devoted to the Laplacian, and in the discrete setting we’ll see that this one simple operator can be applied to a diverse array of tasks (surface parameterization, surface smoothing, vector field design and decomposition, distance computation, fluid simulation… you name it, we got it!).

Fortunately, now that we know how to write \(\mathrm{div}\), \(\mathrm{grad}\) and \(\mathrm{curl}\) using exterior calculus, expressing the scalar Laplacian is straightforward: \(\Delta = \star d \star d\). More generally, the \(k\)-form Laplacian is given by

\[ \Delta := \star d \star d + d \star d \star. \]

The name “Laplace-Beltrami” is used merely to indicate that the domain may have some amount of curvature (encapsulated by the Hodge star). Some people like to define the operator \(\delta := \star d \star\), called the codifferential, and write the Laplacian as \(\Delta = \delta d + d \delta\).

One question you might ask is: why is the Laplacian for 0-forms different from the general \(k\)-form Laplacian? Actually, it’s not — consider what happens when we apply the term \(d \star d \star\) to a 0-form \(\phi\): \(\star \phi\) is an \(n\)-form, and so \(d \star \phi\) must be an \((n+1)\)-form. But there are no \((n+1)\)-forms on an \(n\)-dimensional space! So this term is often omitted when writing the scalar Laplacian.

October 31, 2012 | Posted in: Notes | Comments Closed

A Quick and Dirty Introduction to Exterior Calculus — Part III: Hodge Duality

Previously we saw that a \(k\)-form measures the (signed) projected volume of a \(k\)-dimensional parallelpiped. For instance, a 2-form measures the area of a parallelogram projected onto some plane, as depicted above. But here’s a nice observation: a plane in \(\mathbb{R}^3\) can be described either by a pair of basis directions \((\alpha,\beta)\), or by a normal direction \(\gamma\). So rather than measuring projected area, we could instead measure how well the normal of a parallelogram \((u,v)\) lines up with the normal of our plane. In other words, we could look for a 1-form \(\gamma\) such that

\[ \gamma(u \times v) = \alpha \wedge \beta(u,v). \]

This observation captures the idea behind Hodge duality: a \(k\)-dimensional volume in an \(n\)-dimensional space can be specified either by \(k\) directions or by a complementary set of \((n-k)\) directions. There should therefore be some kind of natural correspondence between \(k\)-forms and \((n-k)\)-forms.

The Hodge Star

Let’s investigate this idea further by constructing an explicit basis for the space of 0-forms, 1-forms, 2-forms, etc. — to keep things manageable we’ll work with \(\mathbb{R}^3\) and its standard coordinate system \((x^1, x^2, x^3)\). 0-forms are easy: any 0-form can be thought of as some function times the constant 0-from, which we’ll denote “\(1\).” We’ve already seen the 1-form basis \(dx^1, dx^2, dx^3\), which looks like the standard orthonormal basis of a vector space:

What about 2-forms? Well, consider that any 2-form can be expressed as the wedge of two 1-forms:

\[ \alpha \wedge \beta = (\alpha_i dx^i) \wedge (\beta_j dx^j) = \alpha_i \beta_j dx^i \wedge dx^j. \]

In other words, any 2-form looks like some linear combination of the basis 2-forms \(dx^i \wedge dx^j\). How many of these bases are there? Initially it looks like there are a bunch of possibilities:

dx^1 \wedge dx^1 & dx^1 \wedge dx^2 & dx^1 \wedge dx^3 \\
dx^2 \wedge dx^1 & dx^2 \wedge dx^2 & dx^2 \wedge dx^3 \\
dx^3 \wedge dx^1 & dx^3 \wedge dx^2 & dx^3 \wedge dx^3 \\

But of course, not all of these guys are distinct: remember that the wedge product is antisymmetric (\(\alpha \wedge \beta = -\beta \wedge \alpha\)), which has the important consequence \(\alpha \wedge \alpha = 0\). So really our table looks more like this:

0 & dx^1 \wedge dx^2 & -dx^3 \wedge dx^1 \\
-dx^1 \wedge dx^2 & 0 & dx^2 \wedge dx^3 \\
dx^3 \wedge dx^1 & -dx^2 \wedge dx^3 & 0 \\

and we’re left with only three distinct bases: \(dx^2 \wedge dx^3\), \(dx^3 \wedge dx^1\), and \(dx^1 \wedge dx^2\). Geometrically all we’ve said is that there are three linearly-independent “planes” in \(\mathbb{R}^3\):

How about 3-form bases? We certainly have at least one:

\[ dx^1 \wedge dx^2 \wedge dx^3. \]

Are there any others? Again the antisymmetry of \(\wedge\) comes into play: many potential bases are just permutations of this first one:

\[ dx^2 \wedge dx^3 \wedge dx^1 = -dx^2 \wedge dx^1 \wedge dx^3 = dx^1 \wedge dx^2 \wedge dx^3, \]

and the rest vanish due to the appearance of repeated 1-forms:

\[ dx^2 \wedge dx^1 \wedge dx^2 = -dx^2 \wedge dx^2 \wedge dx^1 = 0 \wedge dx^1 = 0. \]

In general there is only one basis \(n\)-form \(dx^1 \wedge \cdots \wedge dx^n\), which measures the usual Euclidean volume of a parallelpiped:

Finally, what about 4-forms on \(\mathbb{R}^3\)? At this point it’s probably pretty easy to see that there are none, since we’d need to pick four distinct 1-form bases from a collection of only three. Geometrically: there are no four-dimensional volumes contained in \(\mathbb{R}^3\)! (Or volumes of any greater dimension, for that matter.) The complete list of \(k\)-form bases on \(\mathbb{R}^3\) is then

  • 0-form bases: 1
  • 1-form bases: \(dx^1\), \(dx^2\), \(dx^3\)
  • 2-form bases: \(dx^2 \wedge dx^3\), \(dx^3 \wedge dx^1\), \(dx^1 \wedge dx^2\)
  • 3-form bases: \(dx^1 \wedge dx^2 \wedge dx^3\),

which means the number of bases is 1, 3, 3, 1. In fact you may see a more general pattern here: the number of \(k\)-form bases on an \(n\)-dimensional space is given by the binomial coefficient

\[ \left( \begin{array}{c} n \\ k \end{array} \right) = \frac{n!}{k!(n-k)!} \]

(i.e., “\(n\) choose \(k\)”), since we want to pick \(k\) distinct 1-form bases and don’t care about the order. An important identity here is

\[ \left( \begin{array}{c} n \\ k \end{array} \right) = \left( \begin{array}{c} n \\ n-k \end{array} \right), \]

which, as anticipated, means that we have a one-to-one relationship between \(k\)-forms and \((n-k)\)-forms. In particular, we can identify any \(k\)-form with its complement. For example, on \(\mathbb{R}^3\) we have

\star\ 1 &=& dx^1 \wedge dx^2 \wedge dx^3 \\
\star\ dx^1 &=& dx^2 \wedge dx^3 \\
\star\ dx^2 &=& dx^3 \wedge dx^1 \\
\star\ dx^3 &=& dx^1 \wedge dx^2 \\
\star\ dx^1 \wedge dx^2 &=& dx^3 \\
\star\ dx^2 \wedge dx^3 &=& dx^1 \\
\star\ dx^3 \wedge dx^1 &=& dx^2 \\
\star\ dx^1 \wedge dx^2 \wedge dx^2 &=& 1

The map \(\star\) (pronounced “star”) is called the Hodge star and captures this idea that planes can be identified with their normals and so forth. More generally, on any flat space we have

\[ \star\ dx^{i_1} \wedge dx^{i_2} \wedge \cdots \wedge dx^{i_k} = dx^{i_{k+1}} \wedge dx^{i_{k+2}} \wedge \cdots \wedge dx^{i_n}, \]

where \((i_1, i_2, \ldots, i_n)\) is any even permutation of \((1, 2, \ldots, n)\).

The Volume Form

So far we’ve been talking about measuring volumes in flat spaces like \(\mathbb{R}^n\). But how do we take measurements in a curved space? Let’s think about our usual example of a surface \(f: \mathbb{R}^2 \supset M \rightarrow \mathbb{R}^3\). If we specify a region on our surface via a pair of unit orthogonal vectors \(u, v \in \mathbb{R}^2\), it’s clear that we don’t want the area \(dx^1 \wedge dx^2(u,v)=1\) since that just gives us the area in the plane. Instead, we want to know what a unit area looks like after it’s been “stretched-out” by the map \(f\). In particular, we said that the length of a vector \(df(u)\) can be expressed in terms of the metric \(g\):

\[ |df(u)| = \sqrt{df(u) \cdot df(u)} = \sqrt{g(u,u)}. \]

So the area we’re really interested in is the product of the lengths \(|df(u)||df(v)| = \sqrt{g(u,u)g(v,v)}\). When \(u\) and \(v\) are orthonormal the quantity \(\det(g) := g(u,u)g(v,v)-2g(u,v)\) is called the determinant of the metric, and can be used to define a 2-form \(\sqrt{\det(g)} dx^1 \wedge dx^2\) that measures any area on our surface. More generally, the \(n\)-form

\[ \omega := \sqrt{\det(g)} dx^1 \wedge \cdots \wedge dx^n \]

is called the volume form, and will play a key role when we talk about integration.

On curved spaces, we’d also like the Hodge star to capture the fact that volumes have been stretched out. For instance, it makes a certain amount of sense to identify the constant function \(1\) with the volume form \(\omega\), since \(\omega\) really represents unit volume on the curved space:

\[ \star 1 = \omega \]

The Inner Product on \(k\)-Forms

More generally we’ll ask that any \(n\)-form constructed from a pair of \(k\)-forms \(\alpha\) and \(\beta\) satisfies

\[ \alpha \wedge \star \beta = \langle\langle \alpha, \beta \rangle\rangle \omega, \]

where \(\langle\langle \alpha, \beta \rangle\rangle = \sum_i \alpha_i \beta_i\) is the inner product on \(k\)-forms. In fact, some authors use this relationship as the definition of the wedge product — in other words, they’ll start with something like, “the wedge product is the unique binary operation on \(k\)-forms such that \(\alpha \wedge \star \beta = \langle\langle \alpha, \beta \rangle\rangle \omega\),” and from there derive all the properties we’ve established above. This treatment is a bit abstract, and makes it far too easy to forget that the wedge product has an extraordinarily concrete geometric meaning. (It’s certainly not the way Hermann Grassmann thought about it when he invented exterior algebra!). In practice, however, this identity is quite useful. For instance, if \(u\) and \(v\) are vectors in \(\mathbb{R}^3\), then we can write

\[ u \cdot v = \star\left(u^\flat \wedge \star v^\flat\right), \]

i.e., on a flat space we can express the usual Euclidean inner product via the wedge product. Is it clear geometrically that this identity is true? Think about what it says: the Hodge star turns \(v\) into a plane with \(v\) as a normal. We then build a volume by extruding this plane along the direction \(u\). If \(u\) and \(v\) are nearly parallel the volume will be fairly large; if they’re nearly orthogonal the volume will be quite shallow. (But to be sure we really got it right, you should try verifying this identity in coordinates!) Similarly, we can express the Euclidean cross product as just

\[ u \times v = \star(u^\flat \wedge v^\flat)^\sharp, \]

i.e., we can create a plane with normal \(u \times v\) by wedging together the two basis vectors \(u\) and \(v\). (Again, working this fact out in coordinates may help soothe your paranoia.)

| Posted in: Notes | Comments Closed

Homework 2: Normals of Discrete Surfaces

[Homework 2 has been released (below) and is due next Thursday by 5pm. We suspect there will be a lot of questions about C++, libDDG, etc., so please do not hesitate to ask questions either in the comments (preferred) or via email. Finally, some of the exercises require knowledge of exterior calculus that we haven't quite seen yet in class -- we will cover this material next time, but for those who want to get an early start I will post some additional notes that you can read through at your leisure. Good luck!]

For a smooth surface in \(\mathbb{R}^3\), the normal direction is easy to define: it is the unique direction orthogonal to all tangent vectors — in other words, it’s the direction sticking “straight out” of the surface. For discrete surfaces the story is not so simple. If a mesh has planar faces (all vertices lie in a common plane) then of course the normal is well-defined: it is simply the normal of the plane. But if the polygon is nonplanar, or if we ask for the normal at a vertex, then it is not as clear how the normal should be defined.

In practice there are a number of different possibilities, which arise from different ways of looking at the smooth geometry. But before jumping in, let’s establish a few basic geometric facts.


Area of a Polygon

Here’s a simple question: how do you compute the area of a polygon in the plane? Suppose our polygon has vertices \(p_1, p_2, \ldots, p_n\). One way to compute the area is to stick another point \(q\) in the middle and sum up the areas of triangles \(q,p_i,p_{i+1}\) as done on the left:

A cute fact about life is that if we place \(q\) anywhere and sum up the signed triangle areas, we still recover the polygon area! (Signed area just means negative if our vertices are oriented clockwise; positive if they’re counter-clockwise.) You can get an idea of why this happens just by looking at the picture: positive triangles that cover “too much” area get accounted for by negative triangles.

The proof is an application of Stokes’ theorem — consider a different expression for the area \(A\) of a planar polygon \(P\):

\[ A = \int_P dx \wedge dy. \]

Noting that \(dx \wedge dy = d(xdy) = -d(ydx)\), we can also express the area as

\[ A = \frac{1}{2} \int_P d(xdy) - d(ydx) = \frac{1}{2} \int_{\partial P} xdy - ydx, \]

where we’ve applied Stokes’ theorem in the final step to convert our integral over the entire surface into an integral over just the boundary. Now suppose that our polygon vertices have coordinates \(p_i = (x_i,y_i)\). From here we can explicitly work out the boundary integral by summing up the integrals over each edge \(e_{ij}\):

\[ \int_{\partial P} xdy - ydx = \sum \int_{e_{ij}} xdy - ydx. \]

Since the coordinate functions \(x\) and \(y\) are linear along each edge (and their differentials \(dx\) and \(dy\) are therefore constant), we can write these integrals as

\sum \int_{e_{ij}} xdy - ydx
&=& \sum \frac{x_i+x_j}{2}(y_j-y_i) - \frac{y_i+y_j}{2}(x_j-x_i) \\
&=& \frac{1}{2} \sum (p_i + p_j) \times (p_j - p_i) \\
&=& \frac{1}{2} \sum p_i \times p_j - p_i \times p_i - p_j \times p_j - p_j \times p_i \\
&=& \sum p_i \times p_j.

In short, we’ve shown that the area of a polygon can be written as simply

\[ A = \frac{1}{2} \sum_i p_i \times p_j. \]

Exercise 2.1 Complete the proof by showing that for any point \(q\) the signed areas of triangles \((q,p_i,p_{i+1})\) sum to precisely the expression above.


Vector Area

A more general version of the situation we just looked at with polygon areas is the vector area of a surface patch \(f: M \rightarrow \mathbb{R}^3\), which is defined as the integral of the surface normal over the entire domain:

\[ N_\mathcal{V} := \int_M N dA. \]

A very nice property of the vector area is that it depends only on the shape of the boundary \(\partial M\) (as you will demonstrate in the next exercise). As a result, two surfaces that look very different (such as the ones above) can still have the same vector area — the physical intuition here is that the vector area measures the total flux through the boundary curve.

For a flat region the normal is constant over the surface and we get just the usual area times the unit normal vector. Things get more interesting when the surface is not flat — for instance, the vector area of a circular band is zero since opposing normals cancel each-other out:

Exercise 2.2 Using Stokes’ theorem, show that the vector area can be written as

\[ N_\mathcal{V} = \frac{1}{2} \int_{\partial M} f \wedge df, \]

where the product of two vectors in \(\mathbb{R}^3\) is given by the usual cross product \(\times\).


Gradient of Triangle Area

Here’s another fairly basic question: consider a triangle sitting in \(\mathbb{R}^3\), and imagine that we’re allowed to pull on one of its vertices \(p\). What’s the quickest way to increase its area \(A\)? In other words, what’s the gradient of \(A\) with respect to \(p\)?


Exercise 2.3 Show that the area gradient is given by

\[ \nabla_p A_\sigma = \frac{1}{2} \mathbf{u}^{\perp} \]

where \(\mathbf{u^\perp}\) is the edge vector across from \(p\) rotated by an angle \(\pi/2\) in the plane of the triangle (such that it points toward \(p\)). You should require only a few very simple geometric arguments — there’s no need to write things out in coordinates, etc.

Vertex Normals from Area Gradient

Ok, with these facts out of the way let’s take a look at some different ways to define vertex normals. There are essentially only two definitions that arise naturally from the smooth picture: the area gradient and the volume gradient; we’ll start with the former.

The area gradient asks, “which direction should we `push’ the surface in order to increase its total area \(A\) as quickly as possible?” Sliding all points tangentially along the surface clearly doesn’t change anything: we just end up with the same surface. In fact, the only thing we can do to increase surface area is move the surface in the normal direction. The idea, then, is to define the vertex normal as the gradient of area with respect to a given vertex.

Since we already know how to express the area gradient for a single triangle \(\sigma\), we can easily express the area gradient for the entire surface:

\[ \nabla_p A = \sum_\sigma \nabla A_\sigma. \]

Of course, a given vertex \(p\) influences only the areas of the triangles touching \(p\). So we can just sum up the area gradients over this small collection of triangles.


Exercise 2.4 Show that the gradient of surface area with respect to vertex \(p_i\) can be expressed as

\[ \nabla_p A = \frac{1}{2}\sum_j (\cot\alpha_j + \cot\beta_j)( p_j - p_i ) \]

where \(p_j\) is the coordinate of the \(j\)th neighbor of \(p_i\) and \(\alpha_j\) and \(\beta_j\) are the angles across from edge \((p_i,p_j)\).


Mean Curvature Vector

The expression for the area gradient derived in the last exercise shows up all over discrete differential geometry, and is often referred to as the cotan formula. Interestingly enough this same expression appears when taking a completely different approach to defining vertex normals, by way of the mean curvature vector \(HN\). In particular, for a smooth surface \(f: M \rightarrow \mathbb{R}^3\) we have

\[ \Delta f = 2HN \]

where \(H\) is the mean curvature, \(N\) is the unit surface normal (which we’d like to compute), and \(\Delta\) is the Laplace-Beltrami operator (see below). Therefore, another way to define vertex normals for a discrete surface is to simply apply a discrete Laplace operator to the vertex positions and normalize the resulting vector.

The question now becomes, “how do you discretize the Laplacian?” We’ll take a closer look at this question in the future, but the remarkable fact is that the most straightforward discretization of \(\Delta\) leads us right back to the cotan formula! In other words, the vertex normals we get from the mean curvature vector are precisely the same as the ones we get from the area gradient.

This whole story also helps us get better intuition for the Laplace-Beltrami operator \(\Delta\) itself. Unfortunately, there’s no really nice way to write \(\Delta\) — the standard coordinate formula you’ll find in a textbook on differential geometry is \(\Delta \phi = \frac{1}{\sqrt{|g|}} \frac{\partial}{\partial x^i}(\sqrt{|g|} g^{ij} \frac{\partial}{\partial x^j}\phi)\), where \(g\) is the metric. However, this obfuscated expression provides little intuition about what \(\Delta\) really does, and is damn-near useless when it comes to discretization since for a triangle mesh we never have a coordinate representation of \(g\)! Later, using exterior calculus, we’ll see that the (0-form) Laplacian can be expressed as \(\Delta = \star d \star d\), which leads to a fairly straightforward discretization. But for now, we’ll make use of a nice tool we learned about earlier: conformal parameterization. Remember that if \(f\) is a conformal map, then lengths on \(M\) and lengths on \(f(M)\) are related by a positive scaling \(e^{u}\). In other words, \(|df(X)| = e^{u}|X|\) for some real-valued function \(u\) on \(M\). Moreover, a conformal parameterization always exists — in other words, we don’t have to make any special assumptions about our geometry in order to use conformal coordinates in proofs or other calculations. The reason conformal coordinates are useful when talking about Laplace-Beltrami is that we can write \(\Delta\) as simply a rescaling of the standard Laplacian in the plane, i.e., as the sum of second partial derivatives divided by the metric scaling factor \(e^{2u}\):

\[ \Delta \phi = \frac{d(d\phi(X))(X) + d(d\phi(Y))(Y)}{e^{2u}}, \]

where \(X\) and \(Y\) are any pair of unit, orthogonal directions.

What’s the geometric meaning here? Remember that for a good old-fashioned function \(\phi: \mathbb{R} \rightarrow \mathbb{R}\) in 1D, second derivatives basically tell us about the curvature of a function, e.g., is it concave or convex?

Well, since \(\Delta\) is a sum of second derivatives, it’s no surprise that it tells us something about the mean curvature!


Exercise 2.5 Show that the relationship \(\Delta f = 2HN\) holds.


Vertex Normals from Volume Gradient

An alternative way to come up with normals is to look at the volume gradient. Suppose that our surface encloses some region of space with total volume \(\mathcal{V}\). As before, we know that sliding the surface along itself tangentially doesn’t really change anything: we end up with the same surface, which encloses the same region of space. Therefore, the quickest way to increase \(\mathcal{V}\) is to again move the surface in the normal direction. A somewhat surprising fact is that, in the discrete case, the volume gradient actually yields a different definition for vertex normals than the one we got from the area gradient. To express this gradient, we’ll use three-dimensional versions of of our “basic facts” from above.

First, much like we broke the area of a polygon into triangles, we’re going to decompose the volume enclosed by our surface into a collection of tetrahedra. Each tetrahedron includes exactly one face of our discrete surface, along with a new point \(q\). For instance, here’s what the volume might look like in the vicinity of a vertex \(p\):

Just as in the polygon case the location of \(q\) makes no difference, as long as we work with the signed volume of the tetrahedra. (Can you prove it?)

Next, what’s the volume gradient for a single tetrahedron? One way to write the volume of a tet is as

\[ \mathcal{V} = \frac{1}{3} A h, \]

where \(A\) is the area of the base triangle and \(h\) is the height. Then using the same kind of geometric reasoning as in the triangle case, we know that

\[ \nabla_p \mathcal{V} = \frac{1}{3} A N, \]

where \(N\) is the unit normal to the base.

To express the gradient of the enclosed volume with respect to a given vertex \(p\), we simply sum up the gradients for the tetrahedra containing \(p\):

\[ \nabla_p \mathcal{V} = \sum_i \mathcal{V}_i = \frac{1}{3} \sum_i A_i N_i. \]

At first glance this sum does not lead to a nice expression for \(\Delta_p \mathcal{V}\) — for instance, it uses the normals \(N_i\) of faces that have little to do with our surface geometry. However, remember that we can place \(q\) anywhere we please and still get the same expression for volume. In particular, if we put \(q\) directly on top of \(p\), then the \(N_i\) and \(A_i\) coincide with the normals and areas (respectively) of the faces containing \(p\) from our original surface:

Exercise 2.6 Show that the volume gradient points in the same direction as the vector area \(N_\mathcal{V}\) (i.e., show that they are the same up to a constant factor).


Other Definitions

So far we’ve only looked at definitions for vertex normals that arise from some smooth definition. This way of thinking captures the essential spirit of discrete differential geometry: relationships from the smooth setting should persist unperturbed in the discrete setting (e.g., \(\Delta f = 2HN\) should be true independent of whether \(\Delta\), \(H\), and \(N\) are smooth objects or discrete ones). Nonetheless, there are a number of common definitions for vertex normals that do not have a known origin in the smooth world. (Perhaps you can find one?)


Uniform Weighting

Perhaps the simplest way to get vertex normals is to just add up the neighboring face normals:

\[ N_U := \sum_i N_i \]

The main drawback to this approach is that two different tessellations of the same geometry can produce very different vertex normals, as illustrated above.

Tip-Angle Weights

A simple way to reduce dependence on the tessellation is to weigh face normals by their corresponding tip angles \(\theta\), i.e., the interior angles incident on the vertex of interest:

\[ N_\theta := \sum_i \theta_i N_i \]


Sphere-Inscribed Polytope

Here’s another interesting approach to vertex normals: consider the sphere \(S^2\) consisting of all points unit distance from the origin in \(\mathbb{R}^3\). A nice fact about the sphere is that the unit normal \(N\) at a point \(x \in S^2\) is simply the point itself! I.e., \(N(x) = x\). So if we start out with a polytope whose vertices all sit on the sphere, one reasonable way to define vertex normals is to simply use the vertex positions.

In fact, it’s not too hard to show that the direction of the normal at a vertex \(p_i\) can be expressed purely in terms of the edge vectors \(e_j = p_j – p_i\), where \(p_j\) are the immediate neighbors of \(p_i\). In particular, we have

\[ N_S = \frac{1}{c} \sum_{j=0}^{n-1} \frac{e_j \times e_{j+1}}{|e_j|^2 |e_{j+1}|^2} \]

where the constant \(c \in \mathbb{R}\) can be ignored since we’re only interested in the direction of the normal. (For a detailed derivation of this expression, see Max, “Weights for Computing Vertex Normals from Facet Normals.”) Of course, since this expression depends only on the edge vectors it can be evaluated on any mesh (not just those inscribed in a sphere).


Coding Assignment

Implement the following methods in libDDG:

  • Vertex::normalEquallyWeighted()
    Purpose: returns unit vertex normal using uniform weights \(N_U\)
  • Vertex::normalAreaWeighted()
    Purpose: returns unit vertex normal using face area weights \(N_\mathcal{V}\)
  • Vertex::normalAngleWeighted()
    Purpose: returns unit vertex normal using tip angle weights \(N_\theta\)
  • Vertex::normalMeanCurvature()
    Purpose: returns unit mean curvature normal \(\Delta f\)
  • Vertex::normalSphereInscribed()
    Purpose: returns unit sphere-inscribed normal \(N_S\)

(The definitions for these methods can be found in libddg/src/Vertex.cpp.)

Once you’ve successfully implemented these methods, test them out on the provided meshes. For convenience you can flip through the different methods using the keys [1]-[5]. You can also see what the mesh looks like by viewing it in “wireframe” mode. Do you notice that some definitions work better than others? When? Why? The only thing you need to submit for the coding assignment (via email) is your modified version of Vertex.cpp. (Of course, if you modify other files you should submit these too! For instance, you might find it convenient to add a method HalfEdge::cotan() that computes the cotangent of the angle across from a given half edge.)


Example meshes can be found here: ddg_hw2_meshes.zip

Code for Next Assignment Posted

Code for the next assignment can be downloaded here:


We’ll provide more details about the actual assignment later this week (including the due date), but for now just try to see if you can get the code to compile. Instructions on getting started can be found in the user guide — this precise sequence of steps may not work for you, but it should at least get you started. If you’re having trouble (or finding wild success!), feel free to gripe about it in the comments (below).

If you’ve never seen C++ before, don’t freak out — the path is pretty clearly laid out for you in the assignment.

Finally, we will be arranging a session (probably some evening later this week — with food!) to help setup the code and get everything compiling/linking. Please leave a comment below indicating which evenings/times work for you.