# How to get a normalized gradient with FreeFem++?

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);
}

• Is it possible to compute $\nabla\cdot X$ by hand then evaluate it for a given $u$ in this software and solve the Poisson equation for $\phi$? Commented Apr 2, 2023 at 21:39
• No, because that term cannot be computed analytically only approximated numerically. Commented Apr 2, 2023 at 21:53
• I mean that you can compute $\nabla \cdot (\nabla u/\lvert\nabla u\rvert)$ with pencil and paper to get an expression involving the first and second derivatives of $u$, then this expression can be evaluated numerically using the FEM software and you can solve the Poisson equation numerically for $\phi$ Commented Apr 3, 2023 at 0:13

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.