TriaMesh GeodesicsΒΆ
[1]:
from lapy import TriaMesh
Load a triangle mesh of a flat square (OFF file).
[2]:
T = TriaMesh.read_off("../data/square-mesh.off")
type(T).__name__
[2]:
'TriaMesh'
[3]:
# 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.
[4]:
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.
[5]:
# 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.
[6]:
B.sum()
[6]:
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.
[7]:
plot_tria_mesh(T,Bi*(A*vfunc),plot_edges=True)
This is the same as the corresponding eigenvalue times the eigenfunction.
[8]:
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.
[9]:
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).
[10]:
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.
[11]:
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.
[12]:
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.
[13]:
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.
[14]:
# 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.
[15]:
plot_tria_mesh(T,x,plot_edges=True,plot_levels=True)