How to get a normalized gradient with FreeFem++?

Question

I am trying to use FreeFem++ to solve the heat geodesics algorithm.

The algorithm is:

• solve $\dot u = \Delta u$ at a specific time $t$.
• compute $X = \frac{\nabla u_t}{|\nabla u_t|}$
• solve $\Delta\phi = \nabla \cdot X$

I am having a hard time doing step 2. The documentation does talk about gradients, but it seems to use them in the context of solving linear systems?

This is my current snippet that tries to solve the heat part of the Varadahn equation.

// Parameters
func u0 = 0; // initial condition: u is zero everywhere
real T = 1.0; // final time
real dt = 0.01; // time step

// Define mesh boundary
border C(t=-pi, -pi + 0.3){x=cos(t) -1/(1 + 4 * t * t); y=sin(t);}
border B(t=-pi + 0.3, pi){x=cos(t) -1/(1 + 4 * t * t); y=sin(t);}

// The triangulated domain Th is on the left side of its boundary
mesh pmesh = buildmesh(C(500) + B(500));

// Fespace
fespace Vh(pmesh, P1);
Vh u, v, uold;
real S;

// Problem
problem heatFlow(u,v) =
    int2d(pmesh)( u*v/dt + ( dx(u)*dx(v) + dy(u)*dy(v) ) )
    - int2d(pmesh)(- uold*v/dt)
    + on(C, u=1) + on(B, u=0); // Dirichlet B.C.: on Curve C the field variable u is zero.

// apply the initial condition:
u = u0;

// Time iterations
ofstream ff("sol.dat");
for(real t = 0; t < T; t += dt)
{
    uold = u; //equivalent to u^{n-1} = u^n
    S = t;
    heatFlow;
    plot(u, fill=true);
}

Answer 1

I actually believe you can do this without needing to compute any complicated derivatives at all and only using basically what you have already. In FEM, you need to solve the Poisson equation $\Delta \phi = \nabla\cdot X$ via the weak form. However, writing this equation in the weak form yields $$ \begin{aligned} \Delta \phi &= \nabla\cdot\left(\frac{\nabla u}{|\nabla u|}\right) \\ \implies \int_\Omega\nabla\phi\cdot\nabla v \ dx &= \int_\Omega\frac{\nabla u\cdot\nabla v}{|\nabla u|} \ dx - \int_{\partial\Omega}\frac{\nabla u\cdot\hat{n}}{|\nabla u|}v \ ds \ \ \ \ \forall v\in V, \end{aligned} $$

where I have multiplied through by $-1$ to makes the bilinear forms positive and $\hat{n}$ is the outer normal to the domain. This could all be rewritten explicitly in terms of the $x$ and $y$-derivatives as $$ \int_{\Omega} \phi_xv_x + \phi_y v_y \ dxdy = \int_{\Omega}\frac{u_xv_x + u_y v_y}{\sqrt{u_x^2 + u_y^2}} \ dxdy - \int_{\partial\Omega} \frac{u_x\hat{n}_1 + u_y\hat{n}_2}{\sqrt{u_x^2 + u_y^2}}v \ ds. $$

The LHS is a bilinear form in $\phi$ and $v$ and the RHS is a linear form on $v$. So given a function $u$, one should be able to construct these integrals and solve the FEM system for $\phi$. You will need to figure out how the boundary integral works in your software, but everything else is pretty straightforward.

Additionally, when numerically dealing with terms like $\nabla u/|\nabla u|$, it is often necessary to regularize the problem by instead using $\nabla u /\sqrt{|\nabla u|^2+\varepsilon}$ for some small $\varepsilon$. This keeps your code from dividing by zero at the cost of introducing a tunable hyperparameter.