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;
    plot(u, fill=true);
  • $\begingroup$ 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$? $\endgroup$
    – whpowell96
    Commented Apr 2, 2023 at 21:39
  • $\begingroup$ No, because that term cannot be computed analytically only approximated numerically. $\endgroup$
    – Makogan
    Commented Apr 2, 2023 at 21:53
  • $\begingroup$ 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$ $\endgroup$
    – whpowell96
    Commented Apr 3, 2023 at 0:13

1 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.

  • 1
    $\begingroup$ You will want to integrate by parts on ever cell, rather than globally, to deal with the nonsmoothness of the ratio term. $\endgroup$ Commented Apr 3, 2023 at 16:22

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.