Discrete Laplacian on a triangle mesh

Hi everybody,

I have a triangle mesh with N nodes and T triangles given by the following matrices

n(N,3): which defines the nodes coordinates

t(T,3): which is the connectivity matrix

If I note L as the Discrete Laplacian on my triangle mesh, and I code this

function1= n(:,1).^2;

function2= L*function1;

function2 should be an array of all twos, shouldn't it?

Unfortunately, although I have found some Discrete Laplacian matlab functions in Matlabcentral none of them retrieve a result where function2 should be an array of all twos.

I would appreciate any tip in this matter about how to compute the Discrete Laplacian

Thanks in advance.

Cheers.

Answer 1

[Alec Jacobson]

L discretizes the integrated discrete Laplacian in the sense that:

x' * L * x discretizes Dirichlet energy: ∫ ‖∇x‖² dx. Via Green's idenity this energy is related to the Laplacian as:

∫ ‖∇x‖² dx = -∫ x∆x dx + boundary terms

So, there are at least 3 reasons you won't get all twos:

1) If n,t is a curved surface (as in, all vertices are not on the same plane) then L also encodes the non-Euclidean metric,

2) if n,t has a boundary, then for vertices along the boundary you'll get the integrate Laplacian + boundary terms

3) L computes integrated values, to convert to intergral average, pointwise values you can "unintegrate" (i.e., divide by local area) by multiplying with the inverse of the mass matrix: M \ L * x

Answer 2

[Perry Mason]

Thanks for your answer. I really appreciate.

In order to understand what you said I produced the following uniform spherical (to evoid boundaries) mesh (of 0.9 m radius)

with N=3840 nodes given the vector n(N:,3) and 1280 triangular elements, where the connectivity matrix is t(T,3).

By using gptoolbox, I try to compute the Laplacian () over the spherical surface of the function with the following code (as you suggested in your previous answer)

L = -cotmatrix(n,t);

M = massmatrix(n,t);

f1=n(:,1).^2;

f2=M\L*f1;

I was expecting to be an array of all twos, but this is what I get

Once more I would appreciate any comment or correction about what I am doing.

Cheers

[Alec Jacobson]

See #1 in the answer above. If your surface is curved, then you don't have a flat metric and won't get just 2s.

Try this in the Euclidean domain via:

[V,F] = create_regular_grid(10,10);
L = cotmatrix(V,F);
M = massmatrix(V,F);
M\(L*f1)

everything is 2.0 except the boundary.