TetMesh Geodesics

import numpy as np

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

from lapy import TetMesh
from lapy.plot import plot_tet_mesh

pio.renderers.default = "sphinx_gallery"

First we need a TetMesh, so lets open a cube with 48K tetrahedra and make sure it is oriented consistently.

T = TetMesh.read_vtk("../data/cubeTetra.vtk")
# T.is_oriented()
T.orient_()

--> VTK format         ...
 --> DONE ( V: 9261 , T: 48000 )

Flipped 24000 tetrahedra

np.int64(24000)

Laplace

Next we solve the Laplace eigenvalue problem to get 10 eigenvalues and -vectors/functions.

from lapy import Solver

fem = Solver(T, lump=True)

evals, evec = fem.eigs(10)

TetMesh with regular Laplace
Solver: spsolve (LU decomposition) ...

To better see the first non-constant function also in the interior we slice the cube at x<0.5.

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

evnum = 1
cutting = ["x<0.5"]
# also here we comment all plots to reduce file size
# uncomment and take a look
plot_tet_mesh(
    T,
    vfunc=evals[evnum] * evec[:, evnum],
    plot_edges=None,
    plot_levels=False,
    cutting=cutting,
    edge_color="rgb(50,50,50)",
    html_output=False,
    flatshading=True,
)

Mesh is oriented, nothing to do
Found 3040 triangles on boundary.

Similar to the triangle case, computing the - divergence of the gradient of an eigenfunctions (and multiplying with inv(B)) yields a scaled version of that function.

from lapy.diffgeo import compute_divergence, compute_gradient

grad = compute_gradient(T, evec[:, evnum])
divx = -compute_divergence(T, grad)
vfunc = Bi * divx

cutting = ["x<0.5"]
plot_tet_mesh(
    T,
    vfunc=vfunc,
    plot_edges=None,
    plot_levels=False,
    cutting=cutting,
    edge_color="rgb(50,50,50)",
    html_output=False,
    flatshading=True,
)

Mesh is oriented, nothing to do
Found 3040 triangles on boundary.

In fact, it is scaled by the eigenvalue.

np.max(np.abs(vfunc - (evals[evnum] * evec[:, evnum])))

np.float32(0.006122589)

Geodesics

Now we run a heat diffusion, applying initial heat to the boundary of the cube.

from lapy import heat

tria = T.boundary_tria()
bvert = np.unique(tria.t)

u = heat.diffusion(T, bvert, m=1)
cutting = ["x<0.5"]
plot_tet_mesh(
    T,
    vfunc=u,
    plot_edges=None,
    plot_levels=True,
    cutting=cutting,
    edge_color="rgb(50,50,50)",
    html_output=False,
    flatshading=True,
)

Found 4800 triangles on boundary.
TetMesh with regular Laplace
Matrix Format now:  csc
Solver: spsolve (LU decomposition) ...
Mesh is oriented, nothing to do
Found 3040 triangles on boundary.

You can see that we get level sets that are not evenly spaced and dense along the boundary. Next we compute the gradient of this heat diffusion, normalize it, and compute the divergence of this normalized gradient.

# get gradients
tfunc = compute_gradient(T, u)
# flip and normalize
X = -tfunc / np.sqrt((tfunc**2).sum(1))[:, np.newaxis]
X = np.nan_to_num(X)
# compute divergence
divx = compute_divergence(T, X)

Finally, we need to solve a Poisson equation to obtain a function that has these normalized gradients (and remove the remaining shift).

# compute distance
from scipy.sparse.linalg import splu

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

A, B = fem.stiffness, fem.mass  # computed above when creating Solver

H = A
b0 = -divx

# solve H x = b0
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)

x = x - np.min(x)

Matrix Format now: csc
Solver: spsolve (LU decomp) - performance not optimal ...

cutting = ["x<0.5"]
plot_tet_mesh(
    T,
    vfunc=x,
    plot_edges=None,
    plot_levels=True,
    cutting=cutting,
    edge_color="rgb(50,50,50)",
    html_output=False,
    flatshading=True,
)
max(x), 0.5 * np.sqrt(3.0)

Mesh is oriented, nothing to do
Found 3040 triangles on boundary.

(np.float32(0.69931495), np.float64(0.8660254037844386))

This results in equally spaced level sets. Instead of solving this manually, we can get the same by simply computing the heat diffusion and the distance function directly.

from lapy import heat
from lapy.diffgeo import compute_geodesic_f

tria = T.boundary_tria()
bvert = np.unique(tria.t)

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

gu = compute_geodesic_f(T, u)

cutting = ["x<0.5"]
plot_tet_mesh(
    T,
    vfunc=gu,
    plot_edges=None,
    plot_levels=True,
    cutting=cutting,
    edge_color="rgb(50,50,50)",
    html_output=False,
    flatshading=True,
)

Found 4800 triangles on boundary.
TetMesh with regular Laplace
Matrix Format now:  csc
Solver: spsolve (LU decomposition) ...
TetMesh with regular Laplace
Matrix Format now: csc
Solver: spsolve (LU decomposition) ...
Mesh is oriented, nothing to do
Found 3040 triangles on boundary.

Finally, we want to explore the gradient and divergence functions a little more. Here we construct the gradient of a function that computes the squared distance to each vertex (x^2+y2+z^2). As the color of each tetrahedon we set the z component of the gradient which should be 2z (or you could try any other value, such as the gradient length).

# test function is squared distance to each vertex
v1func = T.v[:, 0] * T.v[:, 0] + T.v[:, 1] * T.v[:, 1] + T.v[:, 2] * T.v[:, 2]

grad = compute_gradient(T, v1func)
# glength = np.sqrt(np.sum(grad * grad, axis=1))
# fcols=glength
fcols = grad[:, 2]
# cutting = ['x<0.5']
cutting = None
plot_tet_mesh(
    T,
    vfunc=None,
    tfunc=fcols,
    plot_edges=None,
    plot_levels=False,
    cutting=cutting,
    edge_color="rgb(50,50,50)",
    html_output=False,
)

Mesh is oriented, nothing to do
Found 4800 triangles on boundary.

Now let’s look at the divergence. While the gradient is constant for each tetrahedron, the divergence is a scalar function again, summing up the partial derivatives of the gradient components. In our case it should be 2+2+2=6.

divx = compute_divergence(T, grad)
divx2 = Bi * divx
cutting = ["z<0.5"]
plot_tet_mesh(
    T,
    vfunc=divx2,
    plot_edges=True,
    plot_levels=False,
    cutting=cutting,
    edge_color="rgb(50,50,50)",
    html_output=False,
    flatshading=True,
    caxis=[0, 8],
)

Mesh is oriented, nothing to do
Found 3040 triangles on boundary.

divx2[5000:5010]

array([5.9999948, 6.0000215, 6.0000215, 5.999988 , 6.000053 , 5.999975 ,
       5.9999676, 6.000024 , 6.000013 , 6.000008 ], dtype=float32)