TetMesh Geodesics¶
[1]:
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.
[2]:
T = TetMesh.read_vtk("../data/cubeTetra.vtk")
# T.is_oriented()
T.orient_()
--> VTK format ...
--> DONE ( V: 9261 , T: 48000 )
Flipped 24000 tetrahedra
[2]:
np.int64(24000)
Laplace¶
Next we solve the Laplace eigenvalue problem to get 10 eigenvalues and -vectors/functions.
[3]:
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.
[4]:
# 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
[5]:
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.
[6]:
from lapy.diffgeo import compute_divergence, compute_gradient
grad = compute_gradient(T, evec[:, evnum])
divx = -compute_divergence(T, grad)
vfunc = Bi * divx
[7]:
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.