# Computing discrete laplacian matrix for mesh fairing

I asked this question on the math stack exchange and got an answer, but I am just as utterly confused as before. My fundamental goal is to actually construct the matrix, that is, a series of steps I can tell an actual computer to build it up.

The gist of it is that there is a 2004 paper called "An Intuitive Framework for Real-Time Freeform Modeling". I am interested in equation (5) on the paper.

I am trying to build that matrix, but it is not obvious to me from the paper alone what the actual coefficients need to be. The answer that I got in the linked post seems correct but it spends most of the time talking about the derivation of the specific matrix, but right now I don't care about the FEM theory used to derive the matrix, all I am trying to do is implement the algorithm. I just want to know what coefficients go where.

I know that many of the rows of $$L$$ (matrix in the paper) will have the coefficients specified by the discrete Laplace-Beltrami operator.

My current understanding is as follows:

Assume that your vertices are sorted such that $$k = m + n$$ where $$k$$ is the total number of vertices, $$m$$ the number of free vertices and $$n$$ the number of fixed vertices. We will assume that $$i < m$$ means that the vertex $$v_i$$ is free, and fixed otherwise.

The way to construct the matrix $$L$$ for a single level (i.e. the left most image in figure 2) would be:

for vertex in free_vertices:
for neighbour in beighbours of vertex:
wij = cotan_weights(vertex, neighbour)
W[vertex, neighbour] = wij
W[vertex, vertex] += -wij
area = vertex_area(vertex)
M[vertex, vertex] = 1.0 / vertex_area

for vertex in fixed_vertices:
M[vertex, vertex] = 1.0
W[vertex, vertex] = 1.0


And then the full Laplacian is $$L = M \times W$$, and it has dimensions $$k\times k$$. Is this at all correct? If this is incorrect, what is the right way to construct the matrix?

## Additional info

I am trying to get this to work in the simplest possible example I could think of, this simple small cylinder:

In this case the top and bottom vertices are fixed the middle are allowed to move.

In this case I get the following matrix, if I try it with my formulation:

-4.2, 0.7, 0.0, 0.7, 0.0, 0.0, 0.0, 1.4, 0.0, 0.0, 0.0, 1.4,
0.7, -4.2, 0.7, 0.0, 0.0, 0.0, 1.4, 0.0, 0.0, 0.0, 1.4, 0.0,
0.0, 0.7, -4.2, 0.7, 0.0, 1.4, 0.0, 0.0, 0.0, 1.4, 0.0, 0.0,
0.7, 0.0, 0.7, -4.2, 1.4, 0.0, 0.0, 0.0, 1.4, 0.0, 0.0, 0.0,
0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0,
0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0,
0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0,
0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0,
0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0,
0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0,


Which is yielding incorrect results:


x in: [0.0, 0.0, 0.0, 0.0, 0.0, -1.0, -0.0, 1.0, 0.0, -1.0, -0.0, 1.0]
x out: [0.0, -21.0, 21.0, 0.0, 0.0, -1.0, -0.0, 1.0, 0.0, -1.0, -0.0, 1.0]
y in: [0.0, 0.0, 0.0, 0.0, -1.0, -0.0, 1.0, 0.0, -1.0, -0.0, 1.0, 0.0]
y out: [28.0, -44.2, 3.5, -0.7, -1.0, -0.0, 1.0, 0.0, -1.0, -0.0, 1.0, 0.0]
z in: [0.0, 0.0, 0.0, 0.0, 2.0, 2.0, 2.0, 2.0, 0.0, 0.0, 0.0, 0.0]
z out: [-28.0, 23.3, -24.5, 0.7, 2.0, 2.0, 2.0, 2.0, 0.0, 0.0, 0.0, 0.0]


I know that libigl implements this in example 401]4, which is exactly what i am trying to achieve. But unfortunately the implementation of the harmonic subroutine is very impregnable and I have not been able to reverse engineer what the matrix is actually supposed to look like.

• How did you decide that your solution is incorrect, because of the large values? Note that your problem is nonlinear. As far as the definition of harmonic goes, it does: $Q = -L$, and then has a loop that concatenates $Q = Q M^{-1}(-L)$. So you basically get $Q = -L(-M^{-1}L)^{k-1}$, $L$ being the cotangent stiffness matrix and $M$ being the mass matrix. Then instead of specifying the boundary conditions explicitly as you do, they use a Lagrange multipliers approach derived from the problem: $\min_u u^TQu, \,:\, Cu = Cb$, which yields $\max_{\lambda}\min_u u^TQu +\lambda (Cb-Cu)$. Commented Sep 21, 2023 at 11:34
• The matrix $A$ is a function of position, so you actually have a non-linear problem. To solve it you can instead try to find the steady state of the time-evolution equation $\partial_t p = \Delta_S p$ with Dirichlet BCs. Discretising time and space in a semi-implicit manner yields: $\frac{M^kp^{k+1}−M^kp^k}{\tau}=L^kp^{k+1}$, thus you can now consecutively solve the problems $(M^k−\tau L^k)p^{k+1}=M^kp^k$ (for the bc indices keep $p^{k+1} = p^k$). You can try a stopping criterion $\|p^{k+1}−p^k\|<\epsilon$. Note that the superscripts denote iterate indices and not powers. Commented Sep 21, 2023 at 13:18

## 1 Answer

I checked your solution, it's wrong. I ran the following MATLAB code:

A = zeros(12,12);
A(1,:) = [-4.2, 0.7, 0.0, 0.7, 0.0, 0.0, 0.0, 1.4, 0.0, 0.0, 0.0, 1.4];
A(2,:) = [0.7, -4.2, 0.7, 0.0, 0.0, 0.0, 1.4, 0.0, 0.0, 0.0, 1.4, 0.0];
A(3,:) = [0.0, 0.7, -4.2, 0.7, 0.0, 1.4, 0.0, 0.0, 0.0, 1.4, 0.0, 0.0];
A(4,:) = [0.7, 0.0, 0.7, -4.2, 1.4, 0.0, 0.0, 0.0, 1.4, 0.0, 0.0, 0.0];
for i=5:12
A(i,i) = 1;
end

x_in = [0.0, 0.0, 0.0, 0.0, 0.0, -1.0, -0.0, 1.0, 0.0, -1.0, -0.0, 1.0]';
y_in = [0.0, 0.0, 0.0, 0.0, -1.0, -0.0, 1.0, 0.0, -1.0, -0.0, 1.0, 0.0]';
z_in = [0.0, 0.0, 0.0, 0.0, 2.0, 2.0, 2.0, 2.0, 0.0, 0.0, 0.0, 0.0]';
x_out = (A\x_in);
y_out = (A\y_in);
z_out = (A\z_in);


The resulting solution is:

x_out = [0.666666666666667, 0,  -0.666666666666667, 0,  0,  -1, 0,  1,  0,  -1, 0,  1]
y_out = [0, 0.666666666666667,  0,  -0.666666666666667, -1, 0,  1,  0,  -1, 0,  1,  0]
z_out = [1.00000000000000,  1.00000000000000,   1.00000000000000,   1.00000000000000,   2,  2,  2,  2,  0,  0,  0,  0]


This clearly doesn't match what you have. It is very likely that your solver has a bug or doesn't work for the matrix above.