# Comptes Rendus
## Mécanique

Rémi Abgrall

**Evaluating a distance function**

Volume 351, Special Issue S1 (2023), p. 5-15

Online since: 27 February 2023

**Part of Special Issue:** The scientific legacy of Roland Glowinski

**Guest editors:** Grégoire Allaire (CMAP, Ecole Polytechnique, Institut Polytechnique de Paris, Palaiseau, France),

Jean-Michel Coron (Laboratoire Jacques-Louis Lions, Sorbonne Université)

and Vivette Girault (Laboratoire Jacques-Louis Lions, Sorbonne Université)

https://doi.org/10.5802/crmeca.155

<figure><img src="image_1.png" /></figure>

This article is licensed under the

CREATIVE COMMONS ATTRIBUTION 4.0 INTERNATIONAL LICENSE.

http://creativecommons.org/licenses/by/4.0/
<figure><img src="image_2.png" /></figure>

# The scientific legacy of Roland Glowinski / L'héritage scientifique de Roland Glowinski

## Evaluating a distance function

Rémi Abgrall*,a

a Institut für Mathematik, Universität Zürich, Winterthurerstrasse 190, CH-8057

Zürich, Switzerland

E-mail: Remi.abgrall@math.uzh.ch

**Abstract.** Computing the distance function to some surface or line is a problem that occurs very frequently. There are several ways of computing a relevant approximation of this function, using for example technique originating from the approximation of Hamilton Jacobi problems, or the fast sweeping method. Here we make a link with some elliptic problem and propose a very fast way to approximate the distance function.

**Keywords.** Hamilton jacobi equation, distance function, Hopf-Cole transformation, elliptic problem, linear solver.

Manuscript received 28 October 2022, revised 4 November 2022 and 15 November 2022, accepted 23 November 2022.

## 1. Introduction

In many cases, one has to evaluate the distance function to a surface $\Gamma_D$ which is part of the boundary of an open set $\Omega \in \mathbb{R}^d$. An example in fluid mechanics is that of turbulence modelling: in some models, one of the parameters in the evaluation of the turbulent viscosity is the distance to the airfoil. Other examples can be found in medical image processing, surface reconstruction, etc. To evaluate the distance function, there are many ways. One technique is to numerically evaluate the viscosity solution of the Eikonal equation,

$$ \begin{array}{ll} \|\nabla u\| - 1 = 0 & x \in \Omega \\ u(x) = 0 & \text{if and only if } x \in \Gamma_D \\ u(x) \ge 0. & \end{array} \qquad (1) $$

In this formulation, only one boundary condition is prescribed, the Dirichlet one on $\Gamma_D$, nothing is said for the other parts of $\partial\Omega$. Numerical techniques for this can be found in [1, 2]. This can be improved by using a fast marching technique in the spirit of [3], or technique coming from computer sciences. From time to time, one can see in the literature methods that compute an approximation of the distance function as the solution of some Laplace equation, an example can be found in [4] and the reference therein. Other references, where elliptic problems are considered, can be found in [5, 6]. They solve (1) by minimizing $(\|\nabla u\| - 1)^2$ using a variational formulation with continuous or discontinuous finite element that needs to be penalized to

<footnote>* Corresponding author.</footnote>
enforce the Dirichlet boundary conditions. It is interesting to notice that their formulation is relatively close to ours, using completely different paths.

There is no obvious links between (1), which is of hyperbolic nature, and an elliptic problem. The purpose of this small note is to provide a link, via the Hopf–Cole transform, and to provide a very fast algorithm (compared to explicit algorithms for computing the viscosity solution of (1) to evaluate the solution of (1)), at least if one accept a small viscosity term (which will nevertheless be present in any numerical method of PDE origin), and discretisation errors.

The format of this note is the following. First, we recall the Hopf–Cole transform and show how it can be applied to the steady Eikonal equation. This leads to an elliptic problem on a function constructed from the distance. We discuss the boundary conditions for this problem, and in particular for the part of $\partial\Omega$ which is not $\Gamma_D$. Then we provide a numerical method, and show its behaviour on unstructured triangular meshes.

## 2. The problem

We want to solve the following problem: Let $\Omega \subset \mathbb{R}^n$ be open, and we set $\partial\Omega = \Gamma_D \cup \Gamma_S$, $\Gamma_D \cap \Gamma_S$ of empty interior. We want to compute the distance function to $\Gamma_D$. We consider the problem of finding the viscosity solution of

$$ 
\begin{array}{l}
\|\nabla u\| - 1 = 0 \quad \mathbf{x} \in \Omega \\
u = 0 \quad \mathbf{x} \in \Gamma_D \\
u = +\infty \quad \mathbf{x} \in \Gamma_S
\end{array}
\qquad (2) $$

Of course $\Gamma_S$ can be empty, but $\Gamma_D$ is never empty by assumption. On $\Gamma_S$ we have set Soner type boundary condition, see [7] for example. For the sake of completeness, we recall the notion of viscosity solution for (2): Let $\varphi \in C^1(\bar{\Omega})$, and $\mathbf{x}_0$ a point where $u - \varphi$ reaches a local minimum: (2) means that if $\mathbf{x}_0 \in \Omega$, $\|\nabla\varphi(\mathbf{x}_0)\| - 1 \ge 0$, and if $\mathbf{x}_0 \in \Gamma_D \min(\|\nabla\varphi(\mathbf{x}_0)\| - 1, u(\mathbf{x}_0)) \ge 0$ while, if $\mathbf{x}_0 \in \Gamma_S$, $\|\nabla\varphi(\mathbf{x}_0)\| - 1 \ge 0$. If $\mathbf{x}_0$ is a minimum of $u - \varphi$, we get: if $\mathbf{x}_0 \in \Omega$, $\|\nabla\varphi(\mathbf{x}_0)\| - 1 \le 0$, if $\mathbf{x}_0 \in \Gamma_D \max(\|\nabla\varphi(\mathbf{x}_0)\| - 1, u(\mathbf{x}_0)) \le 0$ and there is no condition on $\Gamma_S$.

Here we propose a method where we solve a viscous regularisation of (2), or more precisely of

$$ 
\begin{array}{l}
\|\nabla u\|^2 - 1 = 0 \quad \mathbf{x} \in \Omega \\
u = 0 \quad \mathbf{x} \in \Gamma_D \\
u = +\infty \quad \mathbf{x} \in \Gamma_S
\end{array}
\qquad (3) $$

since the two problems have the same solutions, we consider, for $v > 0$, the problem (in the viscosity sense, see [7] for second order problems)

$$ 
\begin{array}{l}
\|\nabla u\|^2 - 1 = v\Delta u \quad \mathbf{x} \in \Omega \\
u = 0 \quad \mathbf{x} \in \Gamma_D \\
u = +\infty \quad \mathbf{x} \in \Gamma_S
\end{array}
\qquad (4) $$

## 3. Rewriting the problem

If, instead of looking at the steady problem (4), we consider the unsteady one,

$$ \frac{\partial u}{\partial t} + \|\nabla u\|^2 - 1 = v\Delta u $$

with the same initial and boundary conditions, this is “almost” the viscous Burgers equation,

$$ \frac{\partial u}{\partial t} + \frac{1}{2} \|\nabla u\|^2 = v\Delta u. $$
for which it is very well known that we can transform it into the heat equation by using the Hopf-Cole transform,

$$u(x, t) = -2v \log(\varphi(x, t)). \qquad (5)$$

The proof is classical (though done in one dimension in most textbooks), but we nevertheless repeat it.

Our notations will be: $\nabla u$ represents the first derivative of the function $u$: for any $h \in \mathbb{R}^d$

$$u(x+h) = u(x) + \nabla u(x) \cdot h + o(h),$$

and $D^2 u$ represents the second derivative (the Hessian) of $u$:

$$\nabla u(x+h) = \nabla u(x) + D^2 u(x) \cdot h + o(h).$$

With this in mind, we have, from (5), that

$$\nabla u(x, t) = -2v \frac{\nabla \varphi(x, t)}{\varphi(x, t)}, \qquad \frac{\partial u}{\partial t} = -2v \frac{\frac{\partial \varphi(x, t)}{\partial t}}{\varphi(x, t)}$$

and

$$D^2 u(x, t) = -2v \frac{D^2 \varphi(x, t)}{\varphi(x, t)} + 2v \frac{\nabla \varphi(x, t) \otimes \nabla \varphi(x, t)}{\varphi(x, t)^2},$$

so that

$$\Delta u = \operatorname{trace}(D^2 u(x, t)) = -2v \frac{\Delta \varphi(x, t)}{\varphi(x, t)} + 2v \frac{\|\nabla \varphi(x, t)\|^2}{\varphi(x, t)^2}.$$

Hence plugging this into the Burgers equation, we have

$$\begin{aligned} \frac{\partial u}{\partial t} + \frac{1}{2} \|\nabla u\|^2 - v \Delta u &= -2v \frac{\frac{\partial \varphi(x, t)}{\partial t}}{\varphi(x, t)} + \frac{1}{2} \left( 4v^2 \frac{\|\nabla \varphi(x, t)\|^2}{\varphi(x, t)^2} \right) \\ &\quad - v \left( -2v \frac{\Delta \varphi(x, t)}{\varphi(x, t)} + 2v \frac{\|\nabla \varphi(x, t)\|^2}{\varphi(x, t)^2} \right) \\ &= - \frac{2v}{\varphi(x, t)} \left( \frac{\partial \varphi(x, t)}{\partial t} - v \Delta \varphi(x, t) \right) \end{aligned}$$

so that in the end we see that $\varphi$ needs to satisfy

$$\frac{\partial \varphi(x, t)}{\partial t} - v \Delta \varphi(x, t) = 0$$

with $\varphi = 1$ on the Dirichlet boundary and $\varphi \ge 0$ on $\Omega$.

Unfortunately, the time dependent problem

$$\frac{\partial u}{\partial t} + \|\nabla u\|^2 - 1 = v \Delta u$$

does not go through as well, but this is not an issue because this is not the problem we want to solve. We want to solve

$$\|\nabla u\|^2 - 1 = v \Delta u$$

for which we again use the change of variable

$$u(x) = \alpha \log(\varphi(x))$$
with $\alpha$ to be determined. We get

$$ \begin{aligned} \|\nabla u\|^2 - 1 - v\Delta u &= \alpha^2 \frac{\|\nabla\varphi(\mathbf{x}, t)\|^2}{\varphi(\mathbf{x}, t)^2} - 1 - v \left( \alpha \frac{\Delta\varphi(\mathbf{x}, t)}{\varphi(\mathbf{x}, t)} - \alpha \frac{\|\nabla\varphi(\mathbf{x}, t)\|^2}{\varphi(\mathbf{x}, t)^2} \right) \\ &= \frac{-1}{\varphi(\mathbf{x}, t)} \left( \alpha v \Delta\varphi(\mathbf{x}, t) + \varphi(\mathbf{x}, t) \right) + \frac{\alpha^2 + v\alpha}{\varphi(\mathbf{x}, t)^2} \|\nabla\varphi(\mathbf{x}, t)\|^2 \end{aligned} \quad (6) $$

so we take $\alpha = -v$ and we need to solve in $\Omega$

$$ v^2 \Delta \varphi(\mathbf{x}, t) = \varphi(\mathbf{x}, t). \quad (7) $$

with the boundary condition

$$ \varphi(\mathbf{x}, t) = 1, \quad \mathbf{x} \in \Gamma_D. \quad (8) $$

On $\Gamma_S$, inspired by what is done for the inviscid problem, and using (6), we see that a condition is

$$ v^2 \Delta \varphi(\mathbf{x}, t) \le \varphi(\mathbf{x}, t) $$

on $\Gamma_S$. This looks a bit like an obstacle problem, but this is not exactly the same (because the "obstacle" is at the boundary). In the next section, inspired by what is done for the Eikonal equation, we will propose a discretisation of this kind of condition. In the numerical section, we will also compare this boundary condition with more natural ones, such as a Neumann condition on the distance function, on $\Gamma_S$.

# 4. Numerical discretisation

## 4.1. Formulation

We consider a triangulation of the polygonal domain $\Omega$ that respects $\Gamma_D$ and $\Gamma_S$. They consist of triangles (or tetrahedrons) that are generically denoted by $K$. The vertices of the triangulation are denoted by $\mathbf{a}_i$, $i = 1, \dots, n_s$. The number of element is $n_e$. For any vertex $\mathbf{a}_i$, $V(j)$ is the set of vertices connected to $\mathbf{a}_i$ by an edge of the triangulation. Often, we make the identification between a vertex $\mathbf{a}_i$ and its index $i$. For the sake of simplicity we only consider the two dimensional case, the three dimensional one can be done in a similar way. The approximation space is

$$ V^h = \{\psi \in H^1(\Omega), \forall K, \psi|_{K} \in \mathbb{P}^1(K)\} \cap \{\psi = 1 \text{ on } \Gamma_D\}. $$

The trial space is

$$ W^h = \{\psi \in H^1(\Omega), \forall K, \psi|_{K} \in \mathbb{P}^1(K)\} \cap \{\psi = 0 \text{ on } \partial\Omega\}. $$

We write the problem as: find $\varphi^h \in V^h$ such that for any $\psi^h \in W^h$,

$$ v^2 \int_{\Omega} \nabla \varphi^h \cdot \nabla \psi^h \, d\mathbf{x} + \int_{\Omega} \varphi^h \psi^h \, d\mathbf{x} = 0 \quad (9a) $$

coupled to boundary conditions on $\Gamma_S$. If $\Gamma_S = \emptyset$, there is nothing more to do.

In the case when $\Gamma_S \neq \emptyset$, we define the boundary conditions according to what is done for the Eikonal equations, see for example [8]. There is defined a numerical Hamiltonian $\mathcal{H}$ which role is to translate the viscosity inequality

$$ \|\nabla u\| - 1 - v^2 \Delta u \ge 0 $$

on $\Gamma_S$ in the limit $v \to 0$. There are several possible versions, but the best (because the gradient of the numerical solution is controlled, see [8]) is to use a Godunov Hamiltonian which amounts to write for any vertex $\mathbf{a}_i \in \Gamma_S$, that

$$ \max_{j \in V(i)} \left( \frac{u_i - u_j}{\|a_i a_j\|} - 1 \right) = 0 $$
that is the distance function $d$ satisfies

$$u_i = \min_{j \in \mathcal{V}(i)} (u_j + \|a_i a_j\|).$$

Keeping in mind that the solution of (9) is related to the solution of (4) by $d = -v \log \varphi$, we will consider the following implementation for the Soner boundary condition: for $\mathbf{a}_i \in \Gamma_S$,

$$u_i = \exp\left(-\frac{\varphi_i}{v}\right), \quad \varphi_i = \max_{j \in \mathcal{V}(i)} (-v \log u_j + \|a_i a_j\|). \qquad (9b)$$

## 4.2. Numerical procedure

We use the following notations. A triangle $K$ has 3 vertices, denoted by $a_i$, $a_j$, $a_k$. We assume that the elements are oriented positively. The gradient of the basis function, $\theta_i$ associated to the vertex $a_i$ is

$$\nabla \theta_{i|K} = \frac{\mathbf{a}_j \mathbf{a}_k^{\perp}}{2|K|}$$

where, for any vector $\mathbf{x}$, $\mathbf{x}^{\perp}$ is orthogonal to $\mathbf{x}$ such that the basis ($\mathbf{x}$, $\mathbf{x}^{\perp}$) is direct. As usual, the angle at $\mathbf{a}_i$ in $K$ is denoted by $\alpha_i^K$.

The variational formulation in $\Omega$ leads to

$$M\varphi + v^2 R\varphi = 0 \qquad (10a)$$

with the boundary condition

$$\varphi_i = 1, \text{ for any } \mathbf{a}_i \in \Gamma_D \qquad (10b)$$

and (9b) on $\Gamma_S$ when this set is not empty. Here, $M$ is the mass matrix and $R$ the rigidity matrix,

$$M_{ij} = \int_{\Omega} \theta_i \theta_j \, d\mathbf{x}, \quad R_{ij} = \int \nabla \theta_i \cdot \nabla \theta_j \, d\mathbf{x}.$$

If $\Gamma_S = \emptyset$, this can be solved by an iterative or a direct solver. Here we have chosen the direct solver PastiX [9]. If $\Gamma_S \neq \emptyset$ the problem becomes non linear. In that case we use an Uzawa-type procedure: we construct a sequence of functions by initialising with $\varphi^0 = 1$ on $\Omega$, and from $\varphi^n$ we construct $\varphi^{n+1}$ by setting:

(1) We compute $\tilde{\varphi}^{n+1}$ solution of (10a) with the Dirichlet boundary

$$\tilde{\varphi}^{n+1} = 1 \text{ on } \Gamma_D \text{ and } \tilde{\varphi}^{n+1} = \varphi^n \text{ on } \Gamma_S. \qquad (11)$$

(2) Then we set

$$\varphi^{n+1} = \widetilde{\varphi^{n+1}} \text{ on } \Omega \setminus \Gamma_S \\ \varphi_i^{n+1} = \exp\left(-\frac{v_i}{v}\right), \quad v_i = \max_{j \in \mathcal{V}(i)} \left(-v \log \widetilde{\varphi^{n+1}}_j + \|a_i a_j\|\right) \text{ on } \Gamma_S. \qquad (12)$$

It is well known that for $\mathbb{P}^1$ approximation on triangular elements, the contribution of $K$ for $R$ is

$$\frac{1}{2} \begin{pmatrix} \cot \alpha_j^K + \cot \alpha_k^K & -\cot \alpha_k^K & -\cot \alpha_i^K \\ -\cot \alpha_k^K & \cot \alpha_i^K + \cot \alpha_k^K & -\cot \alpha_i^K \\ -\cot \alpha_j^K & -\cot \alpha_i^K & \cot \alpha_i^K + \cot \alpha_j^K \end{pmatrix}$$

Since $\cot \alpha + \cot \beta = \frac{\sin(\alpha+\beta)}{\sin \alpha \sin \beta}$, and since $\alpha_i^K + \alpha_j^K + \alpha_k^K = \pi$, the diagonal terms are positive $^1$. The term $R_{ij}$ for two adjacent points is, since $\mathbf{a}_i$ and $\mathbf{a}_j$ defines the common edges between two triangles $K^+$ and $K^-$,

$$R_{ij} = \int_{K^+} \langle \nabla \theta_i, \nabla \theta_j \rangle + \int_{K^-} \langle \nabla \theta_i, \nabla \theta_j \rangle = -\frac{1}{2} (\cot \alpha_k^{K+} + \cot \alpha_k^{K-})$$

<footnote>$^1$A quicker way to see this is to write $R_{ii} = \frac{\|\mathbf{a}_j - \mathbf{a}_k\|^2}{2|K|}$.</footnote>
and it is known that $R_{ij} \le 0$ for $i \ne j$ if and only if

$$ \alpha_k^{K^+} + \alpha_k^{K^-} \le \pi, \qquad (13) $$

see Figure 1.

<figure><img src="image_3.png" /><figcaption>**Figure 1.** Notations.</figcaption></figure>

The mass matrix $M$ is a matrix with positive entries, and the contribution $M_K$ of element $K$ to it is

$$ |K| \begin{pmatrix} 1/6 & 1/12 & 1/12 \\ 1/12 & 1/6 & 1/6 \\ 1/12 & 1/12 & 1/6 \end{pmatrix} $$

If for any vertex, the condition (13) is met, the solution of the auxiliary problem satisfies

$$ \min_{\mathbf{x} \in \Gamma_S} \varphi^n(\mathbf{x}) \le \tilde{\varphi}^{n+1} \le 1, $$

and from (12), we see that

$$ 0 \le \min_{\mathbf{x} \in \Gamma_S} \varphi^{n+1}(\mathbf{x}) \le \varphi^{n+1} \le 1. $$

Similarly, we can show that the sequence is monotone increasing, and since $\varphi^n \ge 0$, it is convergent. The monotone increasing nature comes from $\varphi^1 \le 1 = \varphi^0$ and then we proceed by induction. We have (using the discrete maximum principle thanks to the condition (13)) that $\varphi^{n+1} \ge \varphi^n$ and then using (12), we see that $\varphi^{n+1} \le \tilde{\varphi}^{n+1}$ on $\Gamma_S$. We have thus shown:

**Proposition 1.** If the variational formulation of the Laplace operator satisfies a maximum principle, the sequence $(\varphi^n)_{n \in \mathbb{N}}$ converges. This is in particular true if the triangulation satisfies the angle condition (13).

## 5. Numerical examples

All the calculations have been done on an Imac with 3.5GHz Quad-Core Intel core i7 processors with the version 6.0.2 of PastiX [9]. To report the performance of the solver, for a mesh with 295 296 vertices, 587 520 elements generated by GMSH [10] using the frontal Delaunay option, the symbolic factorisation takes 0.47 s, the evaluation of the non zeros entries of the matrices and the second hand side takes 0.11 s, and the solution takes 9.16 seconds. The averaged maximal band-with of the matrix was 171 781, its maximal band-with is 292 530. The computations have been done sequential. We do only calculations when $\Gamma_S \ne \emptyset$ because they are a priori more complicated. The symbolic factorisation is done once for all (for a given mesh) if the Uzawa-type method is needed.
The first test is the evaluation of the distance function on the annulus $\{\mathbf{x}, 1 \le \|\mathbf{x}\| \le 2\}$. The viscosity is set to $\nu = 0.1$. The Dirichlet condition is set for the inner circle, and the Soner condition on the outer circle. The solution is displayed on Figure 2-a, while the error to the true solution is on Figure 2-b.

<figure><img src="image_4.png" /><figcaption>**Figure 2.** Results for the distance in an annulus. On (a) we have the solution, and on (b) we have $| - \nu \log u_i - d_i |$.</figcaption></figure>

The second case is that of the distance to a body made of two NACA airfoils. The Dirichlet condition is set on the airfoils, and the Soner one on the outer boundary. The results are displayed on Figure 3 as well as the mesh. Here again, $\nu = 0.1$.

We also have considered a case where $\Gamma_D$ and $\Gamma_S$ are not disjoint. The example under consideration is

$$ \Omega = \{ \mathbf{x} = (x_1, x_2) \in \mathbb{R}^2, 1 \le \|\mathbf{x}\| \le 2, x_1 \ge 0 \}. \quad (14a) $$

We take

$$ \Gamma_D = \{(x_1, 0), 1 \le x_1 \le 2\}, \quad \text{and} \quad \Gamma_S = \partial\Omega \setminus \Gamma_S. \quad (14b) $$

We take $\nu = 0.01$. We provide the result on a fine mesh (151 713 points and 301 568 elements) on Figure 4.

In Figure 5, we compare the iterative convergence of the algorithm for several meshes (2473, 9657 and 151 713 vertices).

We observe that we need 65 iterations for the coarse mesh, 130 for the medium one and 485 for the fine one. The ratio of mesh points w.r.t. the coarse one is 1:4:61, and for the iteration the ratios 1:2:7.5, so we see that the cost evolves like $h^{-1}$ (the mesh is very regular). For a classical explicit hyperbolic solver, the cost scales like the number of points, so $h^{-2}$. In this case, and others that we have computed (such as the airfoil case, the conclusion is similar, or even better.

To end this paragraph, let us comment a bit on the boundary condition of Soner type. Since we expect that approximation of the distance has a gradient of norm approximately equal to unity, one may wonder why a Neuman type boundary condition on the distance, say $\nabla u \cdot \mathbf{n} = 1$, would not fit. Indeed, this has been our starting point for imposing boundary conditions on $\Gamma_S$.

If we do that, we first have, since $u = -\nu \log \varphi$,

$$ \frac{\partial u}{\partial n} = \nabla u \cdot \mathbf{n} = -\nu \frac{\nabla \varphi \cdot \mathbf{n}}{\varphi}. $$
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Results for the two NACA problem. The mesh and a zoom of the mesh is displayed.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 4.** Solution for the geometry and the boundary conditions defined by (14). The values range between 0 and 3.907, 30 isolines.</figcaption></figure>
Thus, to set $\frac{\partial u}{\partial n} = 1$, we are let to Robin type conditions on $\varphi$. We have compared our formulation and this one on the fine mesh of the previous test case to the exact solution

$$u(x, y) = \begin{cases} y & \text{if } 1 \le x \le 2 \text{ and } (x, y) \in \Omega \\ d_1 + \theta & \text{else,} \end{cases}$$

where, if $\mathbf{x} = (x, y)$,

$$d_1 = \|\mathbf{x} - \mathbf{p}\| \quad \text{with} \quad \mathbf{p} = \frac{\mathbf{x}}{\|\mathbf{x}\|^2} - \sqrt{1 - \frac{1}{\|\mathbf{x}\|^2} \frac{\mathbf{x}^{\perp}}{\|\mathbf{x}\|}}$$

and $\theta$ is the arclength on the inner circle between $\mathbf{p}$ and (1,0). The results are displayed on Figure 6 with $v = 0.01$.

<figure><img src="image_7.png" /><figcaption>Figure 5. Convergence to the steady solution for several meshes.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Figure 6. (a): isolines for the exact solution (red) and the solution with the Hopf–Cole transform, (b): isolines for the exact solution(red) and the solution with Robin boundary conditions (black). The background colour is that of the exact solution, in both case we have 30 isolines between 0 and $u = 3.826$ (the maximum value of the distance on this domain). The solution with the Neumann condition get values larger than 3.9, this explains why the isolines stop.</figcaption></figure>
From the figure it is clear that the solution with the Neumann condition is completely off compared with our method, though this one is not perfect (because of $v$ not small enough). However, the solution with the Robin condition could be used as an initial guess, this has not been done in our code.

## 6. Conclusion

A final remark is that it is certainly possible to establish rigorous error bounds between the distance function and what is computed here using approximation results between the viscous regularisation of Hamilton Jacobi equations and standard $L^\infty$ error estimates on $\mathbb{P}^1$ approximation of elliptic equations. This has not been done here because we were motivated by designing a working algorithm.

This algorithm has its own drawbacks. The first one is that when $v$ becomes very small, the problem becomes stiffer and stiffer. When the domain is large, the actual value of the solution of the elliptic problem becomes extremely small. It is interesting to note links with large deviation problems (see the last chapter of [7] where exactly the same PDE is studied, for completely different reasons). However, if one comes back to the initial motivation of this work (finding the distance function for turbulence modelling), our experience is that the computation close to the Dirichlet boundary is very reliable, and that applying the Dirichlet condition on all boundaries is enough to get a good approximation.

When the mesh is too distorted so that a discrete maximum principle does not apply, the solution can be slightly above 1 (so that the 'distance' would be negative): in that case, the solution provided by this method can be used as a good initial condition to a 'traditional ' Hamilton Jacobi problem.

## Conflicts of interest

The authors declare no competing financial interest.

## Acknowledgments

This paper is dedicated to Roland Glowinski. He always have been very nice to the author. Roland has also worked a lot on Hamilton Jacobi equations, for example in [11–13]. This plus the handling of [4] as editor of JCP has motivated the present paper.

## References

[1] S. Osher, C.-W. Shu, "High-order essentially nonoscillatory schemes for Hamilton-Jacobi equations", *SIAM J. Numer. Anal.* **28** (1991), no. 4, p. 907-922.

[2] R. Abgrall, "Numerical discretization of the first-order Hamilton-Jacobi equation on triangular meshes", *Commun. Pure Appl. Math.* **49** (1996), no. 12, p. 1339-1373.

[3] J. A. Sethian, "Fast marching methods", *SIAM Rev.* **41** (1999), no. 2, p. 199-235.

[4] H. Nishikawa, "A hyperbolic Poisson solver for wall distance computation on irregular triangular grids", *J. Comput. Phys.* **445** (2021), article no. 110599.

[5] C. Basting, D. Kuzmin, "A minimization-based finite element formulation for interface-preserving level set reinitialization", *Computing* **95** (2013), no. S1, p. 13-25.

[6] T. Adams, S. Giani, W. M. Coombs, "A high-order elliptic PDE based level set reinitialisation method using a discontinuous Galerkin discretisation", *J. Comput. Phys.* **379** (2019), p. 373-391.

[7] G. Barles, *Solutions de viscosité des équations de Hamilton-Jacobi*, Mathématiques & Applications (Berlin), vol. 17, Springer, 1994.
[8] R. Abgrall, "Construction of simple, stable, and convergent high order schemes for steady first order Hamilton-Jacobi equations", *SIAM J. Sci. Comput.* **31** (2009), no. 4, p. 2419-2446.

[9] P. Hénon, P. Ramet, J. Roman, "PaStiX: A High-Performance Parallel Direct Solver for Sparse Symmetric Definite Systems", *Parallel Comput.* **28** (2002), no. 2, p. 301-321.

[10] C. Geuzaine, J.-F. Remacle, "Gmsh: a 3-D finite element mesh generator with built-in pre- and post-processing facilities", *Int. J. Numer. Methods Eng.* **79** (2009), no. 11, p. 1309-1331.

[11] H. Liu, R. Glowinski, S. Leung, J. Qian, "A finite element/operator-splitting method for the numerical solution of the three dimensional Monge-Ampère equation", *J. Sci. Comput.* **81** (2019), no. 3, p. 2271-2302.

[12] R. Glowinski, H. Liu, S. Leung, J. Qian, "A finite element/operator-splitting method for the numerical solution of the two dimensional elliptic Monge-Ampère equation", *J. Sci. Comput.* **79** (2019), no. 1, p. 1-47.

[13] A. Caboussat, R. Glowinski, D. Gourzoulidis, "A least-squares/relaxation method for the numerical solution of the three-dimensional elliptic Monge-Ampère equation", *J. Sci. Comput.* **77** (2018), no. 1, p. 53-78.