TriaMesh Geodesics

from lapy import TriaMesh

Load a triangle mesh of a flat square (OFF file).

T = TriaMesh.read_off("../data/square-mesh.off")
type(T).__name__

'TriaMesh'

# import plotting functions
# import plotly
# plotly.offline.init_notebook_mode(connected=True)
import plotly.io as pio

from lapy.plot import plot_tria_mesh

pio.renderers.default = "sphinx_gallery"

We now plot the triangle mesh with a function overlay of the triangle quality. Note, that this is a function attached to the triangles, not the vertices, so it is piecewise flat.

q = T.tria_qualities()
plot_tria_mesh(T, plot_edges=True, tfunc=q)

Laplace

Next, we check a few properties of eigenfunctions. For this we get the first few and also the stiffness matrix (A) and the lumped mass matrix (B, reduced to a diagonal), of which we can easily compute the inverse.

# compute first eigenfunction
from lapy import Solver

fem = Solver(T, lump=True)
eval, evec = fem.eigs()
vfunc = evec[:, 1]

# also get A,B (lumped), and inverse of B (easy as it is diagonal due to lumping)
A, B = fem.stiffness, fem.mass
Bi = B.copy()
Bi.data **= -1

TriaMesh with regular Laplace-Beltrami
Solver: spsolve (LU decomposition) ...

The mass matrix B represents the inner product so that the integral of the product of two functions x and y over the whole mesh is x B y’. The lumped mass matrix that we use here is a simplified version, as all off-diagonal elements are added to the diagonal. The entries on the diagonal represent the vertex areas and their sum is the total area of the mesh. For our unit square it will be 1.

B.sum()

np.float32(1.0)

Let’s see what happens when we apply the Laplace operator to an eigenfunction. Eigenfunctions are solutions to Delta f = - lambda f, so we should obtain a scaled version of f.

plot_tria_mesh(T, Bi * (A * vfunc), plot_edges=True)

This is the same as the corresponding eigenvalue times the eigenfunction.

plot_tria_mesh(T, eval[1] * vfunc, plot_edges=True)

Laplace is also defined as the -div(grad(f)). So first applying the gradient and then the divergence to an eigenfunction and then multiplying with inv(B) should yield the same result as above again. Note, that multiplying with inv(B) is necessary to get back from the integrated divergence to the original function.

from lapy.diffgeo import compute_divergence, compute_gradient

grad = compute_gradient(T, vfunc)
divx = -compute_divergence(T, grad)
plot_tria_mesh(T, Bi * divx, plot_edges=True)

Geodesics

Now we will replicate the idea of geodesics in heat, where first a heat diffusion is solved and massaged in the right way to yield an approximation of geodesics on the mesh. This also works on curved meshes, but for simplicity we keep using the square here. So let’s start with computing the heat diffusion from the boundary (with default time factor m=1).

from lapy import heat

bvert = T.boundary_loops()
u = heat.diffusion(T, bvert, m=1)

TriaMesh with regular Laplace-Beltrami
Matrix Format now:  csc
Solver: spsolve (LU decomposition) ...

We show some of the level sets. Note, that they are not evenly spaced and get steeper closer to the boundary.

plot_tria_mesh(T, u, plot_edges=True, plot_levels=True)

Next step is to compute the gradient (vector field) of the heat diffusion function and normalize all vectors to unit length.

import numpy as np

# compute gradient of heat diffusion
tfunc = compute_gradient(T, u)

# normalize gradient
X = -tfunc / np.sqrt((tfunc**2).sum(1))[:, np.newaxis]
X = np.nan_to_num(X)

Then we get the integrated divergence of the normalized gradients.

divx = compute_divergence(T, X)

Finally, to obtain the distance function, we need to solve a Poisson equation. The solution can be shifted arbitrary, so we need to subtract the minimum, which should be along the boundary of the mesh.

# compute distance
from scipy.sparse.linalg import splu

useCholmod = True
try:
    from sksparse.cholmod import cholesky
except ImportError:
    useCholmod = False

fem = Solver(T, lump=True)
A, B = fem.stiffness, fem.mass

H = -A
b0 = divx

# solve H x = b0
# we don't need the B matrix here, as divx is the integrated divergence
print("Matrix Format now: " + H.getformat())
if useCholmod:
    print("Solver: cholesky decomp - performance optimal ...")
    chol = cholesky(H)
    x = chol(b0)
else:
    print("Solver: spsolve (LU decomp) - performance not optimal ...")
    lu = splu(H)
    x = lu.solve(b0)

# remove shift
x = x - min(x)

TriaMesh with regular Laplace-Beltrami
Matrix Format now: csc
Solver: spsolve (LU decomp) - performance not optimal ...

In short, the idea is to first get a function (heat diffusion) that flows from a point or from a set of points (boundary) through the mesh, then normalize its gradient, compute the divergence and finally step backward through the Laplace operator to find a function that has this normalized gradient. If we look at it, we actually notice that level sets are equally spaced now.

plot_tria_mesh(T, x, plot_edges=True, plot_levels=True)

Nice, but only an approximation.

# max distance (smoothed)
(max(x), np.sqrt(2) / 2)

(np.float32(0.6049783), np.float64(0.7071067811865476))

Instead of computing the solution ourselves, we could simply employ the standard Poisson solver on inv(B) times divx.

vf = fem.poisson(-Bi * divx)
vf = vf - min(vf)
plot_tria_mesh(T, vf, plot_edges=True, plot_levels=True)

Matrix Format now: csc
Solver: spsolve (LU decomposition) ...

This should give us the same result as what we had earlier.

max(abs(vf - x))

np.float32(2.9802322e-07)

Or we can just call compute_geodesic_f which does all the work for us.

from lapy.diffgeo import compute_geodesic_f

gf = compute_geodesic_f(T, u)
plot_tria_mesh(T, gf, plot_edges=True, plot_levels=True)

TriaMesh with regular Laplace-Beltrami
Matrix Format now: csc
Solver: spsolve (LU decomposition) ...

And verify it is still the same result.

# again should be the same
max(abs(gf - x))

np.float32(0.0)

Similar to normalizing the gradients of a function, we can see if we can rotate it.

# testing if we can rotate the function
from lapy.diffgeo import compute_rotated_f

gf = compute_rotated_f(T, vf)
plot_tria_mesh(T, gf, plot_edges=True, plot_levels=True)

TriaMesh with regular Laplace-Beltrami
Matrix Format now: csc
Solver: spsolve (LU decomposition) ...

Mean Curvature Mapping

To demonstrate the mean curvature mapping to a sphere, we need to have a closed mesh. It should not have too many flat regions (so not the cube) as there is no cuvature.

# Load your mesh here and uncomment. The mesh should have not too many flat regions (not a cube)
# from lapy.diffgeo import tria_mean_curvature_flow
# from lapy.plot import plot_tria_mesh
# T = TriaIO.import_off("../data/???")
# T2 = tria_mean_curvature_flow(T)
# plot_tria_mesh(T2,plot_edges=True,plot_levels=True)