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.