In this project you will implement a conjugate gradient solver and use it to solve a mean curvature flow problem. This will allow you to drastically smooth an input mesh in only a few steps. First, please, take a look at Lecture 6 notes, and Shewchuk's paper here.

The simplest way to start this homework for you is you implement mean curvature flow using explicit time stepping. This entails finding a mean curvature normal K at each vertex, and positioning each vertex such that:

x^{t+1} = x^{t} - dt * K^{t}

The mean curvature K is defined by the cotangent formula, as defined in equation 14 of this paper. The dt must be VERY small to prevent your mesh from exploding. The {t} values are the current values and the {t+1} values are the values after smoothing.

The real assignment for this homework, then, is to do implicit time stepping, using the backward Euler method (watch carefully where the time indices go):

x^{t+1} = x^{t} - dt * K^{t+1}

In order to solve this, you will need to solve a matrix equation. Ultimately this matrix equation can be solved by using the conjugate gradient technique, as described in Shewchuck's paper. In order to carry out the solution process, matrix multiplies are used. You can either create a sparse matrix class to do this, or "multiply" by your mesh. That is, your mesh can be treated as a matrix, if you think about it hard enough.

The equation that will be solved is presented in the implicit fairing paper:

(I - dt * K)X^{t+1} = X^{t}

You want to reformat this equation such that your matrix will feature coefficients, which are then multiplied by the positions X^{t+1} such that we get an equation that makes sense. You will ahve to work it out yourself, but feel free to alk to your TAs and friends about it to make sure you understand before you start galloping along. Here a re a few hints for you, though.

In the ith row, we want to find coefficients to fit the equation x_i^{t+1} - dt * K_i^t = x_i^{t}. K_i^t is defined in terms of cotangent coefficients. Thus, the matrix will have -dt/4A times the cotangent terms that correspond to x_j in the jth column. Again, if this isn't clear for you, make sure to work it out by hand.

One problem with this system is that it's not symmetric, since there's the 1/4A in the equations. Thus, if we multiply each row (both the row of the matrix and the right hand side vector) by 4A, we will be able to fix it so it's symmetric again, and can be solved by Shewchuck's method as is.

It is important to fix your boundaries, or find some smarter way to handle them, because otherwise your mesh will eat itself...