Compute Bi-harmonic weights over a triangular mesh
C++ code to solve the Biharmonic equation over a triangle mesh - 05/2019 - #Geometry
In this article I provide [ 
C++ code ] to produce a bi-harmonic weight map over a triangle mesh by solving the Bi-harmonic equation.
This is an extension of my more in depth article about harmonic weights, which I recommend to read first, unless you are already familiar with the topic.
The continuous form of the Bi-harmonic equation present itself as follows:
$$ \begin{array}{lll} \Delta \Delta \mathbb f & = & 0 \annotation{\text{Laplacian operator applied twice}} \\ \Delta^2 \mathbb f & = & 0 \annotation{\text{Is the bi-laplacian}} \\ \end{array} $$While its discrete form can be deduced in a straightforward manner, since we've seen in my article on harmonic weights that the Laplacian *operator* is described by \(\mathbf L \ \mathbf M^{-1}\):
$$ \begin{array}{lll} \mathbf M^{-1} \mathbf L \quad \mathbf M^{-1} \ \mathbf L \mathbf x & = & [ 0 ] \\ \\ \phantom{\mathbf M^{-1}} \mathbf L \quad \mathbf M^{-1} \mathbf L \ \mathbf x & = & [ 0 ] \annotation{\text{ multiply both sides by } \mathbf M}\\ \end{array} $$While simple this method may lead to numerical instabilities, this is especially true for higher order harmonices (tri-harmonic etc.) Alec Jacobson thesis may be a good starting point to solve these issues if they arise in your application
Also big thanks to Alec who took the time to write a quick note on the topic which I'll leave here:
There are two issues that you're running into. One is mathematical and the other
is numerical.
Here's a version of the narrative that skips a lot of the Finite Element Method
behind it, but will ultimately lead you to a good solution.
The best way to understand polyharmonic weights is to consider the energies
they minimize. In the smooth setting the correspondence goes like this:
∫ ‖∇u‖² dA <--> ∆u = 0
∫ (∆u)² dA <--> ∆²u = 0
∫ ‖∇∆u‖² dA <--> ∆³u = 0
...
In the discrete setting, our primary building blocks will be two matrices: L and
M. These matrices discretize the following energies:
∫ uv dA ~~~> u' M v
∫ ∇u⋅∇v dA ~~~> u' L v
L is the "usual" symmetric cotangent matrix. And M is the mass matrix (often a
diagonal matrix of vertex areas). Often we think of L as the Laplacian, but
notice that here it's acting as integrating the product of the gradients of the
two functions u and y. We recognize it as the Laplacian because Green's identity
tells us this equality:
-∫ u∆v dA + boundary terms = ∫ ∇u⋅∇v dA ~~~> u' L v
So now we can think of L as computing *integrated* Laplacian of y times u.
Skipping the FEM behind it, we can "unintegrate" this by applying M⁻¹ so that:
∆y ~~~> M^{-1}L y
Notice the integral is gone and here we're dealing with pointwise quantities.
Now we have enough to go back to the original energies and discretize each of
them:
∫ ‖∇u‖² dA ~~~> u'Lu
∫ (∆u)² dA ~~~> (M^{-1}L u)' M (M^{-1}L) u = u'LM^{-1}L u
∫ ‖∇∆u‖² dA ~~~> (M^{-1}L u)' L (M^{-1}L) u = u' L M^{-1}L M^{-1}L u
Minimize these is as easy as stripping off the u' and setting the equation to 0
on the right hand side.
However... multipying all these L's and M^{-1}s will start to create a really
poorly conditioned matrix. I've also observed that for triharmonic and higher
the results of using a direct solver like Cholesky can be poor. So a better
solution is to solve the equivalent but better conditioned problems. For example
for the biharmonic equation we should minimize:
min u'LM^{-1}L u
u
we can introduce a auxiliary vector of variables v and solve the constrained
problem:
min y M y subject to Lu = My
u,y
Using the Lagrange multiplier method, introduce another vector of variables λ
and solve the unconstrained saddle problem:
saddle y M y + λ(Lu - My)
u,y,λ
taking derivatives to u,y, and λ and setting them to zero reveals the solution:
/ 0 0 L \ / u \ / 0 \
| 0 M -M | | y | = | 0 |
\ L -M 0 / \ λ / \ 0 /
Notice that the middle equation just says that y=λ, so this can be factored out
leaving:
/ 0 L \ / u \ / 0 \
\ L -M / \ y / = \ 0 /
Now, you could further factor out y = M^{-1}L u, and this would lead to the
LM^{-1}L system above. But the point is to keep the system matrix bigger (and
sparser) and better conditioned.
Solving this "unflattened" system will (in my experience) produce better
results, purely due to better numerics. Algebraically the solutions are the
same.
You can similarly work this out for triharmonic and so on.
This is skipping the FEM theory which you could find in my thesis among other
places. Hopefully this helps.
-Alec
No comments