Mean Curvature Flow

Using the discrete Laplacian operator, we can push each vertex to the average position of its neighbors to smooth the mesh. This technique is known as the mean curvature flow. More formally, for each dimension $d$, we seek to find a solution $f^d: \mathcal{M} \to \mathbb{R}$ that satisfies the mean curvature flow, $\frac{\partial f^d}{\partial t} = \Delta f^d$. Then the final vertex positions $f$ is the concatenation of the solutions in each dimension, $f = [f^x, f^y, f^z]$.

Similar to the way we discretize the Laplacian operator, we can discretize this PDE using the hat functions $\psi_i$ as the basis. For example, to solve for $f$ at vertex $v_j$, we can integrate the PDE with the test function $\psi_j$ to obtain,

$\sum_i \frac{\partial f_i}{\partial t} \int \psi_i \cdot \psi_j \ dA = \sum_i f_i \int \Delta \psi_i \cdot \psi_j \ dA = - \sum_i f_i \int \nabla \psi_i \cdot \nabla \psi_j \ dA$

where in the last expression we use integration by parts and ignore the boundary terms.

We define a "mass" matrix $D$ whose entries are the integral: $D_{ij} = \int \psi_i \cdot \psi_j \ dA$. Intuitively, $D_{ij}$ decides how much weight we give for the edge $e_{ij}$ when solving for $f$.

As usual, we define the entries of Laplacian $L$ as the integral: $L_{ij} = -\int \nabla \psi_i \cdot \nabla \psi_j \ dA$.

Then we can simply the expression a bit:

$\sum_{i} \frac{\partial f_i}{\partial t} D_{ij} = \sum_{i} f_i L_{ij}$

We can discretize the time derivative using an unit time $\delta t$ to get $\frac{\partial f_i}{\partial t} = \frac{f^{t+1}_i - f^{t}_i}{ \delta t}$. Then we can convert the PDF into an implicit Euler integration scheme:

$\sum_i \left(D_{ij} - \delta t L_{ij}\right) f^{t+1}_i = \sum_i D_{ij} f^t$

Writing the above formula in matrix notation,

$\left(D - \delta t L\right) f^{t+1} = D f^t$

Note that this expression is a sparse linear system that can be efficiently solved using iterative numerical solvers such as conjugate gradient.