## Geodesics in Heat implementation
This 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`.


```python
import sys
sys.path.append('..')
import mesh, differential_operators, sparse_matrices, numpy as np
from tri_mesh_viewer import TriMeshViewer as Viewer
```


```python
volMesh = mesh.Mesh('../../examples/meshes/3D_microstructure.msh', degree=1)
```


```python
# Choose whether to work with the tet mesh or its boundary triangle mesh.
m = volMesh
#m = volMesh.boundaryMesh()
```


```python
# 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]
```

We have not yet bound the sparse matrix manipulation and solver functionality of MeshFEM, so we use scipy for now:


```python
import scipy, scipy.sparse, scipy.sparse.linalg
```

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`)


```python
L = differential_operators.laplacian(m).compressedColumn()
M = differential_operators.mass(m, lumped=False).compressedColumn()
A = L + t * M

mask = np.ones(m.numVertices(), dtype=np.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)))
```


```python
# 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]
```

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}


```python
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]) 
```

Visualize the approximate distance field.


```python
view = Viewer(m, scalarField=heatDist)
view.show()
```

528f0604fe5ef985714f31edc4def5ca94d0edbf
ipynb
Jupyter Notebook
python/examples/GeodesicsInHeat.ipynb
pbedenbaugh/MeshFEM
742d609d4851582ffb9c5616774fc2ef489e2f88
MIT
python/examples/GeodesicsInHeat.ipynb
pbedenbaugh/MeshFEM
742d609d4851582ffb9c5616774fc2ef489e2f88
MIT
python/examples/GeodesicsInHeat.ipynb
pbedenbaugh/MeshFEM
742d609d4851582ffb9c5616774fc2ef489e2f88
MIT
Qwen/Qwen-72B
1. YES
2. YES