Divergence of gradient of temperature in heat method for geodesic distance approximation

Question

The is a paper from some years ago describing how to approximate the geodesic distance on a surface.

The method solve the unsteady heat equation for a certain amount of $t$, then computes the normalized gradient of the temperature field and then performs a second solve of the heat equation, now using the divergence of the normalized gradient as heat input

$$ \Delta u = \frac{dT}{dt} $$ $$ X = \frac{\nabla u}{|\nabla u|} $$ $$ \Delta d = \nabla \cdot X $$

With this algorithm, d is an approximation of the real geodesic distance field.

I have been trying to implement this algorithm using FEM with linear triangular elements. However, I am facing an issue. Since my $u$ approximation is a piecewise linear field, the gradient inside a triangle is constant. This means that $\nabla \cdot X$ is 0, and when applying the FEM for the second solve, the integral of $q$

There is an implementation of this algorithm by the CGAL library (https://doc.cgal.org/latest/Heat_method_3/index.html), which uses piecewise linear functions to approximate the $u$ field, so I am totally confused, as I don't understand how they can compute the divergence of the gradient.

How should I approximate the gradient of the $u$ field so I can then compute the divergence correctly?

Answer 1

It's important to keep track of on which aspects of the mesh the various discretized quantities are defined over. The space discretized $u$ for example is defined on the vertices of the mesh, and from these vertex values a discretized normalized gradient $X := \frac{\nabla u}{|\nabla u|}$ is calculated on each triangular facet. This is what you are saying in that the gradient is constant within each triangle, which is true; because of this if we were to define the divergence $\nabla \cdot X$ also on each face we would always get zero. But notice that when solving the poisson equation in step 1 it was useful to have both sides be vertex valued (this is what is used in finite element solvers). We might therefore try to see if we can define a suitable vertex valued discretization of $\nabla \cdot X$. The best way to do this is to think of it as an average or integrated version of the divergence based on the divergence theorem; this is the heart of the finite volume method. So at each vertex $v$ we define the discretized $\nabla\cdot X$ there to be the average divergence of the composite cell made up of all triangles containing $v$ as a vertex. By the divergence theorem this is equivalent to a line integral along the edges of these trangles not containing $v$, but since the line integral around any individual triangle is zero due to there being a constant $X$ on each facet this is equivalent to evaluating the line integrals along all the edges of the adjacent triangles that do contain $v$. This is how it is calculated in the link you share. This kind of averaged divergence is quite common since it leads to a conservative scheme, one where there is cross talk among neighboring cells to ensure the quantities change in a cooperative way

Answer 2

According to the divergence theorem:

$$ \int_\Omega \psi \nabla \cdot \boldsymbol{\Gamma} d\Omega = \int_{\partial\Omega} \psi (\boldsymbol{\Gamma} \cdot \boldsymbol{n} )dS - \int_\Omega \boldsymbol{\Gamma} \cdot \nabla \psi d\Omega $$

So I can obtain the weak formulation for the second solve, and integrate the constant vector field $\boldsymbol{X}$ using the divergence theorem, which is just a matter of taking into account the integral on the triangle, but also, the integral along the boundaries of the triangles.

In practice, I have found that the boundary integral, at least for the problem I am solving, can be neglected.