<jupyter_start><jupyter_text>Geodesics in Heat implementationThis notebook implements Geodesics in Heat [[Crane et al. 2014]](https://arxiv.org/pdf/1204.6216.pdf) for triangle and tet meshes.Compare to the C++ implementation in `experiments/geodesic_heat/main.cc`.<jupyter_code>import sys
sys.path.append('..')
import mesh, differential_operators, sparse_matrices, numpy as np
from tri_mesh_viewer import TriMeshViewer as Viewer
volMesh = mesh.Mesh('../../examples/meshes/3D_microstructure.msh', degree=1)
# Choose whether to work with the tet mesh or its boundary triangle mesh.
m = volMesh
#m = volMesh.boundaryMesh()
# Choose a timestep proportional to h^2 where h is the average edge length.
# (As discussed in section 3.2.4 of the paper)
c = 4 / np.sqrt(3)
t = c *  m.volume / m.numElements()
# Choose source vertex/vertices for computing distances
sourceVertices = [0]<jupyter_output><empty_output><jupyter_text>We have not yet bound the sparse matrix manipulation and solver functionality of MeshFEM, so we use scipy for now:<jupyter_code>import scipy, scipy.sparse, scipy.sparse.linalg<jupyter_output><empty_output><jupyter_text>Backwards Euler time stepping for heat equation $\frac{\mathrm{d}}{\mathrm{d}t} = \bigtriangleup u, \, \, u|_\gamma = 1 \, \forall t$:\begin{align}                         \frac{u_t - u_0}{t} &= \bigtriangleup u_t \\ \Longrightarrow \quad M \frac{u_t - u_0}{t} &= -L u_t    \quad \text{(positive FEM Laplacian discretizes $-\bigtriangleup$)} \\ \Longrightarrow \quad \underbrace{(M + t L)}_A u_t &= M u_0 \end{align}where $\gamma$ is the domain from which we wish to compute distances (here given by `sourceVertices`)<jupyter_code>L = differential_operators.laplacian(m).compressedColumn()
M = differential_operators.mass(m, lumped=False).compressedColumn()
A = L + t * M

mask = np.ones(m.numVertices(), dtype=bool)
mask[sourceVertices] = False

A_ff = A[:,  mask][mask, :]
A_fc = A[:, ~mask][mask, :]

# Solve (M + t L) u = 0 with the constraint u[sourceVertices] = 1
u = np.ones(m.numVertices())
u[mask] = scipy.sparse.linalg.spsolve(A_ff, -A_fc @ np.ones(len(sourceVertices)))
# Compute the heat gradients
g = differential_operators.gradient(m, u)
# Normalize the gradients to get an approximate gradient of the distance field
X = -g / np.linalg.norm(g, axis=1)[:, np.newaxis]<jupyter_output><empty_output><jupyter_text>Fit a scalar field's gradients to these normalized gradients $X$ by solving a Poisson equation:\begin{align}- \bigtriangleup \phi = -\nabla \cdot X \quad &\text{in } \Omega \\\frac{\mathrm{d} \phi}{\mathrm{d} {\bf n}} = {\bf n} \cdot X \quad &\text{on } \partial \Omega \\\phi = 0 \quad &\text{on } \gamma\end{align}<jupyter_code>divX = differential_operators.divergence(m, X)
L_ff = L[:, mask][mask, :]
heatDist = np.zeros(m.numVertices())
heatDist[mask] = scipy.sparse.linalg.spsolve(L_ff, divX[mask])<jupyter_output><empty_output><jupyter_text>Visualize the approximate distance field.<jupyter_code>view = Viewer(m, scalarField=heatDist)
view.show()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1390075583305948, 'fraction_numerical': 0.009530069010844562, 'mean_word_length': 4.980353634577603, 'pattern_counts': {'":': 0, '<': 21, '<?xml version=': 0, '>': 21, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '5', 'dir': 'jupyter-structured-clean-dedup', 'id': '149886', 'license': 'MIT', 'n_tokens_mistral': 1098, 'n_tokens_neox': 1046, 'n_words': 287, 'path': 'python/examples/GeodesicsInHeat.ipynb', 'repo_name': 'MeshFEM/MeshFEM', 'size': '3043'}