# Evaluating a distance function

R. Abgrall

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

##  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{split}\|\nabla u\|-1=0&\text{~{}}\mathbf{x}\in{\Omega}\\\
u(\mathbf{x})=0&\text{ if and only if }\mathbf{x}\in\Gamma_{D}\\\
u(\mathbf{x})\geq 0.\end{split}$ (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 [12, 1]. This can be improved by
using a fast marching technique in the spirit of [13], 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 [11] and the reference
therein. Other references, where elliptic problems are considered, can be
found in [5, 3]. They solve (1) by minimizing $\big{(}\|\nabla u\|-1)^{2}$
using a variational formulation with continuous or discontinuous finite
element that needs to be penalized to 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{split}\|\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{split}$ (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 [4]
for example. For the sake of completeness, we recall the notion of viscosity
solution for (2): Let $\varphi\in C^{1}(\overline{\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\geq
0,$ and if $\mathbf{x}_{0}\in\Gamma_{D}$
$\min\big{(}\|\nabla\varphi(\mathbf{x}_{0})\|-1,u(x_{0})\big{)}\geq 0$ while,
if $\mathbf{x}_{0}\in\Gamma_{S}$ $\|\nabla\varphi(\mathbf{x}_{0})\|-1\geq 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\leq 0,$ if
$\mathbf{x}_{0}\in\Gamma_{D}$
$\max\big{(}\|\nabla\varphi(\mathbf{x}_{0})\|-1,u(x_{0})\big{)}\leq 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{split}\|\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{split}$ (3)

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

$\begin{split}\|\nabla u\|^{2}-1=\nu\Delta u&\quad\mathbf{x}\in\Omega\\\
u=0&\quad\mathbf{x}\in\Gamma_{D}\\\
u=+\infty&\quad\mathbf{x}\in\Gamma_{S}\end{split}$ (4)

##  3 Rewriting the problem

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

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

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

$\dfrac{\partial u}{\partial t}+\frac{1}{2}\|\nabla u\|^{2}=\nu\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(\mathbf{x},t)=-2\nu\log\big{(}\varphi(\mathbf{x},t)\big{)}.$ (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(\mathbf{x}+h)=u(\mathbf{x})+\nabla u(\mathbf{x})\cdot h+o(h),$

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

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

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

$\nabla
u(\mathbf{x},t)=-2\nu\dfrac{\nabla\varphi(\mathbf{x},t)}{\varphi(\mathbf{x},t)},\qquad\dfrac{\partial
u}{\partial t}=-2\nu\dfrac{\dfrac{\partial\varphi(\mathbf{x},t)}{\partial
t}}{\varphi(\mathbf{x},t)}$

and

$D^{2}u(\mathbf{x},t)=-2\nu\dfrac{D^{2}\varphi(\mathbf{x},t)}{\varphi(\mathbf{x},t)}+2\nu\dfrac{\nabla\varphi(x,t)\otimes\nabla\varphi(x,t)}{\varphi(\mathbf{x},t)^{2}},$

so that

$\Delta
u=\text{trace}\big{(}D^{2}u(\mathbf{x},t)\big{)}=-2\nu\dfrac{\Delta\varphi(\mathbf{x},t)}{\varphi(\mathbf{x},t)}+2\nu\dfrac{\|\nabla\varphi(x,t)\|^{2}}{\varphi(\mathbf{x},t)^{2}}.$

Hence plugging this into the Burgers equation, we have

$\begin{split}\dfrac{\partial u}{\partial t}+\frac{1}{2}\|\nabla
u\|^{2}-\nu\Delta
u&=-2\nu\dfrac{\dfrac{\partial\varphi(\mathbf{x},t)}{\partial
t}}{\varphi(\mathbf{x},t)}+\frac{1}{2}\bigg{(}4\nu^{2}\dfrac{\|\nabla\varphi(\mathbf{x},t)\|^{2}}{\varphi(\mathbf{x},t)^{2}}\bigg{)}\\\
&\qquad\qquad-\nu\bigg{(}-2\nu\dfrac{\Delta\varphi(\mathbf{x},t)}{\varphi(\mathbf{x},t)}+2\nu\dfrac{\|\nabla\varphi(x,t)\|^{2}}{\varphi(\mathbf{x},t)^{2}}\bigg{)}\\\
&=-\frac{2\nu}{\varphi(\mathbf{x},t)}\bigg{(}\dfrac{\partial\varphi(\mathbf{x},t)}{\partial
t}-\nu\Delta\varphi(\mathbf{x},t)\bigg{)}\end{split}$

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

$\dfrac{\partial\varphi(\mathbf{x},t)}{\partial
t}-\nu\Delta\varphi(\mathbf{x},t)=0$

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

Unfortunately, the time dependant problem

$\dfrac{\partial u}{\partial t}+\|\nabla u\|^{2}-1=\nu\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=\nu\Delta u$

for which we again use the change of variable

$u(\mathbf{x})=\alpha\log\big{(}\varphi(\mathbf{x})\big{)}$

with $\alpha$ to be determined. We get

$\begin{split}\|\nabla u\|^{2}-1-\nu\Delta
u&=\alpha^{2}\dfrac{\|\nabla\varphi(\mathbf{x},t)\|^{2}}{\varphi(\mathbf{x},t)^{2}}-1-\nu\bigg{(}\alpha\dfrac{\Delta\varphi(\mathbf{x},t)}{\varphi(\mathbf{x},t)}-\alpha\dfrac{\|\nabla\varphi(x,t)\|^{2}}{\varphi(\mathbf{x},t)^{2}}\bigg{)}\\\
&=\frac{-1}{\varphi(\mathbf{x},t)}\bigg{(}\alpha\nu\Delta\varphi(\mathbf{x},t)+\varphi(\mathbf{x},t)\bigg{)}+\dfrac{\alpha^{2}+\nu\alpha}{\varphi(\mathbf{x},t)^{2}}\|\nabla\varphi(\mathbf{x},t)\|^{2}\end{split}$
(6)

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

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

with the boundary condition

$\varphi(\mathbf{x},t)=1,\qquad\mathbf{x}\in\Gamma_{D}.$ (8)

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

$\nu^{2}\Delta\varphi(\mathbf{x},t)\leq\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 Neuman
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 consists of triangles (or tetrahedrons)
that are generically denoted by $K$. The vertices of the triangulation are
denoted by $\mathbf{a}_{i}$, $i=1,\ldots,n_{s}$. The number of element is
$n_{e}$. For any vertex $\mathbf{a}_{i}$, $\mathcal{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}$,

$\nu^{2}\int_{\Omega}\nabla\varphi^{h}\cdot\nabla\psi^{h}\;d\mathbf{x}+\int_{\Omega}\varphi^{h}\psi^{h}\;d\mathbf{x}=0$
(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 [2].
There is defined a numerical Hamiltonian $\mathcal{H}$ which role is to
translate the viscosity inequality

$\|\nabla u\|-1-\nu^{2}\Delta u\geq 0$

on $\Gamma_{S}$ in the limit $\nu\rightarrow 0$. There are several possible
versions, but the best (because the gradient of the numerical solution is
controlled, see [2]) is to use a Godunov Hamiltonian which amounts to write
for any vertex $\mathbf{a}_{i}\in\Gamma_{S}$, that

$\max\limits_{j\in\mathcal{V}(i)}\big{(}\frac{u_{i}-u_{j}}{\|a_{i}a_{j}\|}-1\big{)}=0$

that is the distance function $d$ satisfies

$u_{i}=\min\limits_{j\in\mathcal{V}(i)}\big{(}u_{j}+\|a_{i}a_{j}\|\big{)}.$

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

$u_{i}=\exp\big{(}-\frac{\varphi_{i}}{\nu}\big{)},\quad\varphi_{i}=\max\limits_{j\in\mathcal{V}(i)}\big{(}-\nu\log
u_{j}+\|a_{i}a_{j}\|\big{)}.$ (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}=\dfrac{\mathbf{a}_{j}\mathbf{a}_{k}^{\bot}}{2|K|}$

where, for any vector $\mathbf{x}$, $\mathbf{x}^{\bot}$ is orthogonal to
$\mathbf{x}$ such that the basis $(\mathbf{x},\mathbf{x}^{\bot})$ 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+\nu^{2}R\varphi=0$ (10a) with the boundary condition
$\varphi_{i}=1,\text{ for any }\mathbf{a}_{i}\in\Gamma_{D}$ (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},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. 1.

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

$\widetilde{\varphi^{n+1}}=1\text{ on }\Gamma_{D}\text{ and
}\widetilde{\varphi^{n+1}}=\varphi^{n}\text{ on }\Gamma_{S}.$ (11)

  2. 2.

Then we set

$\begin{split}\varphi^{n+1}&=\widetilde{\varphi^{n+1}}\text{ on
}\Omega\backslash\Gamma_{S}\\\
\varphi_{i}^{n+1}&=\exp\big{(}-\frac{v_{i}}{\nu}\big{)},\quad
v_{i}=\max\limits_{j\in\mathcal{V}(ji}\big{(}-\nu\log\widetilde{\varphi^{n+1}}_{j}+\|a_{i}a_{j}\|\big{)}\text{
on }\Gamma_{S}.\end{split}$ (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_{j}^{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 111a quicker way to see this is to write
$R_{ii}=\frac{\|\mathbf{a}_{j}-\mathbf{a}_{k}\|^{2}}{2|K|}$.. 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}\big{(}\cot\alpha_{k}^{K^{+}}+\cot\alpha_{k}^{K^{-}}\big{)}}$

and it is known that $R_{ij}\leq 0$ for $i\neq j$ if and only if

$\alpha_{k}^{K^{+}}+\alpha_{k}^{K^{-}}\leq\pi,$ (13)

see figure 1.

Figure 1:  Notations.

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\limits_{\mathbf{x}\in\Gamma_{S}}\varphi^{n}(\mathbf{x})\leq\widetilde{\varphi}^{n+1}\leq
1},$

and from (12), we see that

$0\leq\min\limits_{\mathbf{x}\in\Gamma_{S}}\varphi^{n+1}(\mathbf{x})\leq{\varphi}^{n+1}\leq
1.$

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

######  Proposition 4.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 is 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 296vertices, 587 520 elements
generated by GMSH [7] 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}\neq\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\leq\|\mathbf{x}\|\leq 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 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}|$

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$

Figure 3:  Results for the two Naca problem. The mesh and a zoom of the mesh
is displayed.

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\leq\|\mathbf{x}\|\leq
2{,}x_{1}\geq 0\\}.$ (14a) We take $\Gamma_{D}=\\{(x_{1},0),1\leq x_{1}\leq
2\\}\text{ , and }\Gamma_{S}=\partial\Omega\backslash\Gamma_{S}.$ (14b)

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

Figure 4:  Solution for the geometry and the boundary conditions defined by
(14). The values range between $0$ and $3.907$, 30 isolines.

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

Figure 5:  Convergence to the steady solution for several meshes.

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$,

$\dfrac{\partial u}{\partial n}=\nabla
u\cdot{\mathbf{n}}=-\nu\dfrac{\nabla\varphi\cdot{\mathbf{n}}}{\varphi}.$

Thus, to set $\dfrac{\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)=\left\\{\begin{array}[]{ll}y&\text{ if }1\leq x\leq 2\text{ and
}(x,y)\in\Omega\\\ d_{1}+\theta&\text{ else,}\end{array}\right.$

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

$d_{1}=\|\mathbf{x}-\mathbf{p}\|\text{ with
}\mathbf{p}=\frac{\mathbf{x}}{\|\mathbf{x}\|^{2}}-\sqrt{1-\frac{1}{\|\mathbf{x}\|^{2}}}\frac{\mathbf{x}^{\bot}}{\|\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 $\nu=0.01$.

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

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 $\nu$ 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

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 [10, 8, 6]. This plus the handling of [11] as editor of JCP has
motivated the present paper.

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 $\nu$ 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 [4] 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.

## References

  * [1] R. Abgrall.  Numerical discretization of the first-order Hamilton-Jacobi equation on triangular meshes.  Commun. Pure Appl. Math., 49(12):1339–1373, 1996.
  * [2] R. Abgrall.  Construction of simple, stable, and convergent high order schemes for steady first order Hamilton-Jacobi equations.  SIAM J. Sci. Comput., 31(4):2419–2446, 2009.
  * [3] Thomas Adams, Stefano Giani, and William M. Coombs.  A high-order elliptic pde based level set reinitialisation method using a discontinuous galerkin discretisation.  Journal of Computational Physics, 379(373-391), 2019.
  * [4] Guy Barles.  Solutions de viscosité des équations de Hamilton-Jacobi, volume 17 of Mathématiques & Applications (Berlin) [Mathematics & Applications].  Springer-Verlag, Paris, 1994.
  * [5] C. Basting and D. Kuzmin.  A minimization-based finite element formulation for interface-preserving level set reinitialization.  Computing, 95(S1):13–25, 2013.
  * [6] Alexandre Caboussat, Roland Glowinski, and Dimitrios Gourzoulidis.  A least-squares/relaxation method for the numerical solution of the three-dimensional elliptic Monge-Ampère equation.  J. Sci. Comput., 77(1):53–78, 2018.
  * [7] C. Geuzaine and J.-F. Remacle.  Gmsh: a three-dimensional finite element mesh generator with built-in pre- and post-processing facilities.  International Journal for Numerical Methods in Engineering, 79(11):1309–1331, 2009.
  * [8] Roland Glowinski, Hao Liu, Shingyu Leung, and Jianliang Qian.  A finite element/operator-splitting method for the numerical solution of the two dimensional elliptic Monge-Ampère equation.  J. Sci. Comput., 79(1):1–47, 2019.
  * [9] Pascal Hénon, Pierre Ramet, and Jean Roman.  Pastix: A high-performance parallel direct solver for sparse symmetric definite systems.  Parallel Computing, 28(2):301–321, 2002.  https://hal.inria.fr/inria-00346017.
  * [10] Hao Liu, Roland Glowinski, Shingyu Leung, and Jianliang Qian.  A finite element/operator-splitting method for the numerical solution of the three dimensional Monge-Ampère equation.  J. Sci. Comput., 81(3):2271–2302, 2019.
  * [11] H. Nishikawa.  A hyperbolic poisson solver for wall distance computation on irregular triangular grids.  Journal of Computational Physics, 110599, 2021.
  * [12] Stanley Osher and Chi-Wang Shu.  High-order essentially nonoscillatory schemes for Hamilton-Jacobi equations.  SIAM J. Numer. Anal., 28(4):907–922, 1991.
  * [13] J. Sethian.  Fast marching methods.  SIAM Rev., 41(2):199–235, 1999.