Here're some notes to avoid (common? obvious?) pitfalls when implementing conformalized mean curvature flow ("Can Mean-Curvature Flow Be Made Non-Singular?" [Kazhdan et al. 2012]).

Naively building a smoothing operation out of the identity matrix (data term) and a discrete pointwise Laplacian (smoothness term) would lead to iterations looking something like this:

```
L = cotmatrix(V,F);
I = speye(size(L));
U = V;
while true
M = massmatrix(U,F,'barycentric');
U = (I-delta*M\L)\U;
tsurf(F,U)
drawnow;
end
```

This is vary badly scaled as we need to first invert the mass matrix, introducing numerical problems if our mesh is irregular, and then solve the system.

Algebraically we can fix these issues by multiplying both sides of the system by our mass matrix:

```
L = cotmatrix(V,F);
U = V;
while true
M = massmatrix(U,F,'barycentric');
U = (M-delta*L)\(M*U);
tsurf(F,U)
drawnow;
end
```

This solution might converge for very small delta values, but eventually the surface positions U will converge on a single point and our floating point representation will run out of precision. In our context, the problem will be that entries in our diagonal mass matrix will degenerate because triangles in our mesh are degenerating.

So for numerical reasons we must normalize the positions at every step. There are various ways we could do this, here I show normalizing the surface to have unit surface area:

```
L = cotmatrix(V,F);
U = V;
while true
M = massmatrix(U,F,'barycentric');
U = (M-delta*L)\(M*U);
U = U/sqrt(sum(doublearea(U,F))*0.5);
tsurf(F,U)
drawnow;
end
```

This is better, but still might only converge for small delta values and even then eventually diverge. The problem this time is a bit more subtle. Our mesh is unit area, but "floating" in space. Importantly, it's floating away from the origin. Thus to store the coordinates, we're *wasting* precision on each vertex by also storing this global translation. We eventually run into the same problem as before. Thus we should *recenter* the mesh around the origin at each step. Again, many ways to do this, but I show centering the centroid:

```
L = cotmatrix(V,F);
U = V;
while true
M = massmatrix(U,F,'barycentric');
U = (M-delta*L)\(M*U);
U = U/sqrt(sum(doublearea(U,F))*0.5);
U = bsxfun(@minus,U,centroid(U,F));
tsurf(F,U)
drawnow;
end
```

Depending on the quality of the original mesh we should now be able to take quite huge time steps and still see convergence (e.g. `delta = 1e12`

!).

**Important:** Normalizing the surface area and recentering at the origin are necessary for numerics, not just visualization.

**Note:** Convergence means that our mesh flows to the sphere. However, the mesh may still *slide around* on the sphere. The problem being that any MÃ¶bius transformation is *allowed* by the flow. Eventually vertices may group up at a pole on the sphere and our mass matrix will degenerate.