# High order recovery of geometric interfaces from cell-average data

Albert Cohen,\* Olga Mula,† Agustin Somacal\*

January 8, 2024

## Abstract

We consider the problem of recovering characteristic functions $u := \chi_{\Omega}$ from cell-average data on a coarse grid, and where $\Omega$ is a compact set of $\mathbb{R}^d$. This task arises in very different contexts such as image processing, inverse problems, and the accurate treatment of interfaces in finite volume schemes. While linear recovery methods are known to perform poorly, nonlinear strategies based on local reconstructions of the jump interface $\Gamma := \partial\Omega$ by geometrically simpler interfaces may offer significant improvements. We study two main families of local reconstruction schemes, the first one based on nonlinear least-squares fitting, the second one based on the explicit computation of a polynomial-shaped curve fitting the data, which yields simpler numerical computations and high order geometric fitting. For each of them, we derive a general theoretical framework which allows us to control the recovery error by the error of best approximation up to a fixed multiplicative constant. Numerical tests in 2d illustrate the expected approximation order of these strategies. Several extensions are discussed, in particular the treatment of piecewise smooth interfaces with corners.

# 1 Introduction

## 1.1 Reconstruction from cell-averages

We consider the problem of reconstructing a function $u : D \to \mathbb{R}$ defined on a multivariate domain $D \subset \mathbb{R}^d$ from cell averages

$$a_T(u) := \frac{1}{|T|} \int_T u(x)dx, \quad T \in \mathcal{T}, \qquad (1)$$

over a partition $\mathcal{T}$ of $D$. This task occurs in various contexts, the most notable ones being:

1. **Image processing:** here $u$ is the light intensity of an image and $\mathcal{T}$ represents a grid of pixels in dimension $d = 2$ or voxels in dimension $d = 3$. Various processing tasks are facilitated by the reconstruction of the image at the continuous level, for example when applying operations that are not naturally compatible with the pixel grid such as rotations, or when changing the format of the pixel grid such as in super-resolution.

2. **Hyperbolic transport PDE's:** here $u$ is a solution to such an equation and $\mathcal{T}$ is a computational grid, typically in dimension $d = 1, 2$ or $3$. Finite volume schemes evolve the cell average data by computing at each time step the numerical fluxes at the interfaces between each cell. Several such schemes are based on an intermediate step that reconstructs simple approximations to $u$ on each cell and compute the numerical fluxes by applying the transport operator to these approximations.

3. **Inverse Problems:** numerous inversion tasks can be formulated as the recovery of a function from observational data, and this data could typically come in the form of local averages of the type (1).

In this paper, we consider functions defined on the unit cube

$$D := [0, 1]^d,$$

<footnote>\*Laboratoire Jacques-Louis Lions, Sorbonne Université, 4 place Jussieu, 75005 Paris, France (albert.cohen@sorbonne-universite.fr and agustin.somacal@sorbonne-universite.fr)</footnote>

<footnote>†TU Eindhoven, Department of Mathematics and Computer Science, 5612 AZ Eindhoven, Netherlands (o.mula@tue.nl)</footnote>
and partitions $\mathcal{T}_h$ of $D$ based on uniform cartesian meshes, that is, consisting of cells of the form

$$T = h(k + D), \quad k := (k_1, \dots, k_d) \in \{0, \dots, l-1\}^d,$$

where $h := \frac{1}{l} > 0$ is the side-length of each cell in $\mathcal{T}_h$, for some $l > 1$. The cardinality of the partition is therefore

$$n := \#(\mathcal{T}_h) = l^d = h^{-d}.$$

We are thus interested in reconstruction operators $R$ that return an approximation $\tilde{u} = R(a)$ to $u$ from the $n$-dimensional vector $a = a(u) = (a_T(u))_{T \in \mathcal{T}_h} \in \mathbb{R}^n$. The most trivial one is the piecewise constant function

$$\tilde{u} := \sum_{T \in \mathcal{T}_h} a_T(u) \chi_T, \qquad (2)$$

which is for example used in the Godunov finite volume scheme. Elementary arguments show that this reconstruction is first order accurate: if $u \in W^{1,p}(D)$, one has

$$\|u - \tilde{u}\|_{L^p(D)} \leq C h \|\nabla u\|_{L^p(D)} \sim n^{-\frac{1}{d}},$$

where $C$ is a fixed constant, and the exponent in this estimate cannot be improved for smoother functions.

A simple way to raise the order of accuracy is by reconstructing on each cell polynomials of higher degree using neighbouring cell averages. For example, in the univariate case $d = 1$ and for some fixed $m \geq 1$, we associate to each interval $T_k := [kh, (k+1)h]$ the centered stencil consisting of the cells $T_l$ for $l = k - m, \dots, k + m$. Then, there exists a unique polynomial $p_k \in \mathbb{P}_{2m}$ such that

$$a_{T_l}(p_k) = a_{T_l}(u), \quad l = k - m, \dots, k + m.$$

We then define a piecewise polynomial reconstruction by

$$\tilde{u} := \sum_{k=0, \dots, n-1} p_k \chi_{T_k}.$$

This strategy can be generalized to higher dimension $d > 1$ in a straightforward manner: for each cell $T$ we consider the stencil $S = S_T$ of $(2m + 1)^d$ cells centered around $T$ and define the piecewise polynomial reconstruction

$$\tilde{u} := \sum_{T \in \mathcal{T}_h} p_T \chi_T,$$

where $p_T$ is the unique polynomial of degree $2m$ in each variable such that

$$a_{\tilde{T}}(p_T) = a_{\tilde{T}}(u), \quad \tilde{T} \in S_T.$$

For example in the bivariate case $d = 2$, we may use $3 \times 3$ stencils to reconstruct bi-quadratic polynomials on each cell. Standard approximation theory arguments show that these local reconstruction operators now satisfy accuracy estimates of the form

$$\|u - \tilde{u}\|_{L^p(D)} \leq C h^r |u|_{W^{r,p}(D)} \sim n^{-\frac{r}{d}},$$

for $r \leq 2m + 1$.

**Remark 1.1.** Polynomials of odd degree can also be constructed by using non-centered stencils. Also note that non-centered stencils need to be used when approaching the boundary of $D$, but this does not affect the above estimate.

These classical methods are therefore efficient to reconstruct smooth functions with a rate of accuracy that optimally reflects their amount of smoothness. Unfortunately they are doomed to perform poorly in the case of functions $u$ that are piecewise smooth with jump discontinuities across hypersurfaces. For example, a piecewise constant reconstruction will have $O(1)$ error on each cell that is crossed by the interface. Since the amount of such cells is of order $h^{1-d}$, we cannot expect a reconstruction error better than

$$\|u - \tilde{u}\|_{L^p} \gtrsim (h^d h^{1-d})^{\frac{1}{p}} = h^{-\frac{1}{p}} = n^{-\frac{1}{dp}}. \qquad (3)$$
<figure><img src="image_1.png" /><figcaption>Figure 1: ENO-SR in one dimension: the jump point $X_0$ is identified by matching the average on the singular cell with the piecewise polynomial reconstruction.</figcaption></figure>

In particular, this reconstruction has first order accuracy $O(h)$ in the $L^1$ norm.

The use of higher order polynomial cannot improve this rate. In fact, a fundamental obstruction is the fact that all the above methods produce approximations $\tilde{u}$ that depend linearly on $a(u)$, and therefore belong to a linear space of dimension $n$. The bottleneck of such methods for a given class of functions $\mathcal{K}$ is therefore given by the so-called Kolmogorov $n$-width defined by

$$d_n(\mathcal{K})_{L^p} := \inf_{\dim(V_n)=n} \sup_{u \in \mathcal{K}} \min_{v \in V_n} \|u-v\|_{L^p}.$$

Then, it can be shown that for very simple classes $\mathcal{K}$ of discontinuous functions such as those of the form $u = \chi_H$, where $H$ is any half-space passing through $D$, the $n$-width in $L^p$ precisely behaves like $n^{-1/d_p}$, see [4]. In summary, any linear method cannot do much better than the low order piecewise constant method, even for interfaces that are infinitely smooth.

Improving the accuracy in the reconstruction of piecewise smooth functions from cell averages therefore motivates the development and study of nonlinear reconstruction strategies, which is at the heart of this work. We first recall the main existing approaches.

## 1.2 Reconstruction of discontinuous interfaces

One first approach aiming to tackle jump discontinuities while maintaining high order approximation in the smooth regions was proposed for the univariate case $d=1$ by Ami Harten in terms of ENO (Essentially Non Oscillatory) and ENO-SR (Subcell Resolution). The ENO strategy [6] is based on selecting for each cell $T$ a stencil $S_T$ that should not include the cell $T^*$ which contains the point of jump discontinuity. This is achieved by choosing among the stencils that contains $T$ the one where the cell average values have the least numerical variation. For $T \neq T^*$ such adaptively selected stencils will tend to avoid $T^*$.

As a consequence high order reconstruction can be preserved in all cells where $u$ is smooth, and in addition this yields for free a singularity detection mechanism which identifies the singular cell $T^*$ that is avoided from both side by the stencil selection. The ENO-SR strategy [5] then consists in reconstructing in this singular cell by extending the polynomials fitted on both sides until a point for which the resulting average match the observed average. The position of this point can therefore be identified by solving a simple algebraic equation. This strategy is very effective in the univariate case as illustrated in Figure 1.

While the ENO stencil selection can be generalized to higher dimension, the ENO-SR strategy does not have a straightforward multivariate version. This is due to the fact that, instead of a single point, the jump discontinuity to be identified is now a hypersurface (curve in 2d, surface in 3d, ...) that cannot be described by finitely many parameters. The approximate recovery of geometric interfaces from cell-average has been the object of continuous investigation, with the particular focus on functions of the form

$$u = \chi_{\Omega},$$

where $\Omega \subset D$ is a set with boundary $\Gamma = \partial\Omega$ having a certain smoothness. For such characteristic functions,
that could for example represent a two phase flow without any mixing or the evolution of a front, the whole difficulty is concentrated in the recovery of $\Gamma$ since the smooth parts are the trivial constant functions 0 or 1.

In this paper, we denote by

$$S_h := \{T \in \mathcal{T}_h : |T \cap \Omega| \neq 0 \text{ and } |T \cap \Omega^c| \neq 0\},$$

the set of cells that are crossed by the interface $\Gamma$ in a non trivial way. Such cells are termed as singular, the other as regular. Let us note that singular cells are characterized by the property

$$0 < a_T(u) < 1,$$

and can thus be identified from the cell-average data.

Practical computational strategies, often termed as volume of fluid methods, consist in using the cell-average data to reconstruct a local approximation of the interface that can be described by finitely many parameters, such as lines in 2d or planes in 3d. This idea was introduced in [12], and was significantly improved in [14] with the LVIRA algorithm that consists in reconstructing in each singular cell a linear interface whose parameters are found by least-square minimization of the difference between exact and reconstructed cell averages on centered 3 × 3 stencil, in the 2d case.

Noe that this continuous least-square minimization is performed over a nonlinear set. This induces a substantial computational burden that could be avoided through a variant, the ELVIRA algorithm, in which the line selection is made between 6 possible configurations explicitly computed from the cell averages. One main result is the fact that this reconstruction returns precisely the true interface if this one is indeed a straight line. We refer to [13] for a comparative survey on these reconstruction algorithms and to [7, 16, 15, 17] for improvements and applications in the domain of 2d and 3d fluid mechanics.

<figure><img src="image_2.png" /><figcaption>Figure 2: Local approximation of a smooth interface by a line interface.</figcaption></figure>

Intuitively the advantage of locally fitting a line or plane interface is that it has the ability to better approximate the interface $\Gamma$ if it is smooth, so that it is expected to improve the low order of accuracy (3) of linear methods. More precisely, if $\Gamma$ has $C^2$ regularity, on each cell of side-length $h$, it can be approximated with Hausdorff distance $\mathcal{O}(h^2)$ by a line in 2d, a plane in 3d, etc. Therefore, if the locally reconstructed linear interface is optimally fitted, the $\mathcal{O}(1)$ error is observed on a strip of volume $\mathcal{O}(h^{2+d-1})$, see Figure 2. Since the amount of singular cell is of order $\#(S_h) \sim h^{1-d}$, we may hope for a reconstruction error with improved order,

$$\|u - \tilde{u}\|_{L^p} \sim (h^{d+1}h^{1-d})^{\frac{1}{p}} = h^{-\frac{2}{p}} = n^{-\frac{2}{dp}}, \quad (4)$$

that is, we double the rate compared to (3). In particular, we obtain second order accuracy $\mathcal{O}(h^2)$ in the $L^1$ norm.

However, to our knowledge such estimates have never been rigourously proved for the aforementioned method. In addition the approximation power of linear interface is also limited and we cannot hope to improve the above rate for interfaces that are smoother than $C^2$. In this context, the objective of this paper is twofold:

1. introduce a theoretical framework for the rigorous convergence analysis of local interface reconstructions from cell averages,

2. within this framework develop reconstruction methods going beyond linear interfaces and provably achieving higher order of accuracy.
## 1.3 Outline

In this paper, we will essentially work in the bivariate case $d = 2$ which makes the exposition simpler while most of our discussion can be carried over to higher dimension.

The recovery methods that we study are local: on each cell $T$ identified as singular, the unknown function $u = \chi_\Omega$ is approximated by a simpler $\tilde{u} = \chi_{\Omega_T}$ picked from a family $V_n$ that can be described by $n$ parameters and that enjoy certain approximation properties for interfaces having prescribed smoothness. This approximation is computed from the cell averages of a rectangular stencil $S_T$ of $m \ge n$ cells centered around $T$. We begin §2 by giving examples of such families and discussing their approximation properties for prior classes $K_s$ of $\chi_\Omega$ associated to sets $\Omega$ with $C^s$ boundaries for $s > 1$. Our ultimate goal is to develop recovery schemes such that the error between $u$ and $\tilde{u}$ is near optimal in the sense of being bounded (up to a multiplicative constant) by the error of best approximation of $u$ by elements from $V_n$. This, in particular means that the recovery should be exact if the true function $u$ belongs to $V_n$.

We introduce in §3 a first class of recovery strategies which we call Optimization Based Edge Reconstruction Algorithms (OBERA). Similar to LVIRA it is based on least-square fitting of simpler interfaces such as lines, but also circles or polynomials that allow to raise the order of accuracy. We show that near optimality of this recovery is ensured by an inverse stability inequality which can in particular be established for line edges and $3 \times 3$ stencils. Unfortunately, this property seems more difficult to prove when raising the order of accuracy, which also leads to more difficult nonlinear optimization procedures.

As a more manageable alternative, we consider recovery methods that are based on the identification of a certain preferred orientation - vertical or horizontal - for describing the interface by a function in the vicinity of each such cell. This leads us to the second class of recovery methods, termed as Algorithms for Edge Reconstruction using Oriented Stencils (AEROS) which is discussed in §4. It avoids continuous optimization by finding an edge described by a polynomial function $y = p_T(x)$ or $x = p_T(y)$ after having used the previously discussed orientation mechanism. The polynomial $p_T$ is identified by simple linear equations. We show that this process satisfies an exactness and stability property that leads to the near optimal recovery bound. In addition, the rate of convergence can be raised to an arbitrarily high order by raising the degree of $p_T$, however at the price of using larger stencils. The analysis of the orientation selection mechanism, based on the Sobel filter, is postponed to the Appendix, where it is shown that it correctly classifies the cells when $h$ is sufficiently small.

All these methods are numerically tested and compared in §5. The differences in terms of convergence rates are confirmed, and we also compare the computational costs, which are by far less for AEROS than OBERA. We also discuss and test the specific recovery of corner singularities in the interface. Finally, in the case of linear transport, we illustrate the propagation of error when using finite volume schemes based on these local recovery strategies.

An open-source python framework<fnref n="1" /> is made available to show the methods presented here but specially to allow an easy way of creating, testing, and comparing new subcell resolution and interface reconstruction methods without the need to re-implement everything from scratch.

# 2 Numerical analysis of local recovery methods

## 2.1 Local approximation by nonlinear families

The methods that we study in this paper for the recovery of the unknown $u = \chi_\Omega$ are based on local approximations of $u$ on each cell $T \in \mathcal{T}_h$ that is identified as singular by simpler characteristic functions picked from an $n$ parameter family $V_n$.

Let us give three examples that will be used further. We stress that in all such examples $V_n$ is not an $n$-dimensional linear space, but instead should be thought as an $n$-dimensional nonlinear manifold.

**Example 1: linear interfaces.** These are functions of the form $v = \chi_H$ where $H$ is a half-plane with a line interface $L = \partial H$. Such functions are described by $n = 2$ parameters. One convenient description is by the pair $(r, \theta)$, where $r \ge 0$ is the offset distance between the center $z_T$ of the cell $T$ of interest and the linear interface and $\theta \in [0, 2\pi[$ is the angle between this line and the horizontal axis. In other words, in this case we

<footnote n="1">¹https://github.com/agussomacal/SubCellResolution</footnote>
use

$$V_2 := \{v_{r,\theta} := \chi_{\{(z-z_T, e_\theta) \le r\}}, \theta \in [0, 2\pi[; r > 0\},$$

where $e_\theta = (-\sin(\theta), \cos(\theta))$. Of course, the $d$-dimensional generalization by half-spaces is a $d$-parameter family where the unit normal vector $e_\theta$ lives on the $d-1$-dimensional unit sphere.

**Example 2: circular interfaces.** These are functions of the form $v = \chi_D$ or $v = \chi_{D^c}$ where $D$ is a disc with circular interface, and $D^c$ its complement. It is easily seen that the corresponding space $V_3$ is now a 3 parameter family, and its $d$-dimensional generalization by characteristic function of balls and their complements is a $d+1$ parameter family. The idea of using circles instead of lines is to increase the approximation capability. However, as we will see, the next family turns out to be more effective both from the point of view of theoretical analysis and computational simplicity.

**Example 3: oriented graphs.** These are functions of the form $v = \chi_P$ where $P$ is the subgraph or the epigraph of a function $p \in W_n$, either applied to the coordinate $x$ or $y$, where $W_n$ is a linear space of univariate functions. In other words, $P$ is given by one among the four equations

$$y \le p(x), \quad y \ge p(x), \quad x \le p(y), \quad x \ge p(y). \tag{5}$$

Of course, the corresponding space $V_n$ is an $n$-parameter family that is not a linear space, while $W_n$ is. Note that the linear interfaces of Example 1 are a particular case where $W_2$ is the set of affine functions. Raising the order of accuracy will be achieved by taking for $W_n$ the space of polynomials of degree $n-1$. The $d$-dimensional generalisation is obtained by taking for $W_n$ a linear space of functions of $d-1$ variables and considering $P$ to be defined by one of the $2d$ equations

$$x_i \le p(x_1, \dots, x_{i-1}, x_{i+1}, \dots, x_d), \quad x_i \ge p(x_1, \dots, x_{i-1}, x_{i+1}, \dots, x_d), \quad i = 1, \dots, d,$$

for some $p \in W_n$. In particular $W_n$ could be a space of multivariate polynomials.

**Example 4: piecewise linear interface.** These are functions of the form $u = \chi_{H_1 \cap H_2}$ where $H_1$ and $H_2$ are two half planes. Therefore the interface consists of two half lines that touch at a corner point $x_0$. The corresponding space $V_4$ is a 4 parameter family, for example by considering the coordinates $x_0 = (x_1, x_2)$ and the angles $\theta_1$ and $\theta_2$ of the normal vectors to the two lines. The goal of this family is to better approximate piecewise smooth interfaces that have corner singularities. Note that Example 1 may be viewed as a particular case where the two half-planes coincide and there is no corner point.

Given such a family $V_n$ and a set $S \subset D$, we denote by

$$e_n(u)_S = \min_{v \in V_n} \|u - v\|_{L^1(S)},$$

the error of best approximation on $S$ measured in the $L^1$-norm.

**Remark 2.1.** Throughout this paper, we shall systematically measure error in $L^1$ norm which is the most natural since in the case of $u = \chi_\Omega$ and $\tilde{u} = \chi_{\tilde{\Omega}}$, this error is simply the area of the symmetric difference between domains, that is,

$$\|u - v\|_{L^1(D)} = |\Omega \Delta \tilde{\Omega}| = |\Omega \cup \tilde{\Omega} - \Omega \cap \tilde{\Omega}|.$$

Note however that estimates in $L^p$ norms can be derived from $L^1$ estimates in a straightforward manner since

$$\|u - v\|_{L^p(D)} = |\Omega \Delta \tilde{\Omega}|^{1/p} = \|u - v\|_{L^1(D)}^{1/p}.$$

In our analysis of the reconstruction error on a singular cell $T$, we will need to estimate the local approximation error on a stencil $S = S_T$ that consists of finitely many cells surrounding $T$. We shall systematically consider rectangular stencils of symmetric shape centered around $T$, so in the $2d$ case they are of size

$$m = (2k + 1) \times (2l + 1),$$

for some fixed $k, l \ge 1$.
The order of magnitude $e_n(u)_S$ both depends on the type of family $V_n$ that one uses and on the smoothness property of the boundary $\Gamma = \partial\Omega$. We describe these properties by introducing prior classes of characteristic functions $\chi_\Omega$ of sets $\Omega \subset D$ with boundary of a prescribed Hölder smoothness. There exists several equivalent definitions of a $C^s$ domain. We follow the approach from [11], that expresses the fact that the boundary can locally be described by graphs of $C^s$ functions (see also Chapter 4 of [1]).

**Definition 2.2.** Let $s > 0$. A domain $\Omega \subset \mathbb{R}^2$ is of class $C^s$ if and only if there exists an $R > 0$, $P > 1$ and $M > 0$, such that for any point $z_0 \in \partial\Omega$, the following holds: there exists an orthonormal system $(e_1, e_2)$ and a function $\psi \in C^s([-R, R])$ with $\|\psi\|_{C^s} \le M$, taking its value in $[-PR, PR]$ and such that

$$z \in \Omega \iff z_2 \le \psi(z_1),$$

for any $z = z_0 + z_1 e_1 + z_2 e_2$ with $|z_1| \le R$ and $|z_2| \le PR$.

Here, we have used the usual definition

$$\|\psi\|_{C^s} = \sup_{0 \le k \le \lfloor s \rfloor} \|\psi^{(k)}\|_{L^\infty([-R, R])} + \sup_{s, t \in [-R, R]} |s-t|^{\lfloor s \rfloor - s} |\psi^{(\lfloor s \rfloor)}(s) - \psi^{(\lfloor s \rfloor)}(t)|,$$

for the Hölder norm. In the case of integer smoothness, we use the convention that $C^s$ denotes functions with Lipschitz derivatives up to order $s-1$, so that in particular the case $s = 1$ corresponds to domains with Lipschitz boundaries. This definition naturally extends to domains of $\mathbb{R}^d$ with $d > 2$ with $\psi$ now being a $C^s$ function of $d-1$ variables.

We can immediately derive a first local approximation error estimate for the space $V_2$ of linear interfaces from the above Example 1: let $u = \chi_\Omega$ with $\Omega$ a domain of class $C^s$. Then, if $S$ is a $2k+1 \times 2l+1$ stencil centered around a cell $T$ that is crossed by the interface, we apply the above Definition 2.2 taking $z_0 = z_T$ the center of $T$. We assume that the sidelength $h$ is small enough so that the stencil $S$ is contained in the rectangle $\{z = z_0 + z_1 e_1 + z_2 e_2 : , |z_1| \le R, |z_2| \le PR\}$, in which $\partial\Omega$ is described by the graph of the $C^s$ function $\psi$. For $z \in S$, we have in addition that $|z_1| \le C_0 h \le R$ where $C_0$ depends only on $l$ and $k$. Using a Taylor expansion and the smoothness of $\psi$ we find that there exists an affine function $a$ such that

$$\|\psi - a\|_{L^\infty([-C_0h, C_0h])} \le C_1 h^r, \quad r := \min\{s, 2\}, \qquad (6)$$

where $C_1$ depends on $C_0$ and the bound $M$ on the $C^s$ norm of $\psi$. For example, we can take for $a$ the Taylor polynomials of order 1 at $z_1 = 0$ when $s > 1$, which corresponds to match the tangent of the interface at the point $z_0 + \psi(0)e_2$, or of order 0 when $s \le 1$. Then, taking $v = \chi_H \in V_2$, where

$$H := \{z_2 \le a(z_1)\},$$

is the corresponding half-space, it follows that

$$\|u - v\|_{L^1(S)} = |S \cap (\Omega \Delta H)| \le Ch^{r+1},$$

where $C$ depends on $(M, l, k)$. In summary, for the local approximation error of $C^s$ domains by a linear interface we have

$$e_n(u)_S \le Ch^{r+1}, \quad r := \min\{s, 2\}. \qquad (7)$$

The same reasoning in $d$ dimensions delivers a local approximation estimate of order $h^{r+d-1}$.

One way to raise this order of accuracy for smoother domains is to use approximation by circular interfaces from Example 2 since this allows us to locally match the curvature in addition to the tangent. In turn we reach a similar estimate of order $h^{r+1}$ however with $r := \min\{s, 3\}$. One more systematic way of raising the order arbitrarily high is to use approximation by oriented subgraphs from Example 3. This approach is central to the AEROS strategies discussed in §5 and we thus discuss it below in more detail.

We begin with the observation that when $s > 1$, the unit tangent vector varies continuously on $\partial\Omega$ if $\Omega$ is of class $C^s$. It follows that locally around any point $z_0$, this vector remains away either from the horizontal vector (1,0) or from the vertical vector (1,0). This allows us to locally describe the boundary by graphs of functions of the standard cartesian coordinates, as expressed by the following alternate definition of $C^s$ domains.
**Definition 2.3.** Let $s > 1$. A domain $\Omega \subset \mathbb{R}^2$ is of class $C^s$ if and only if there exists an $R > 0$, $P > 1$ and $M > 0$, such that for any point $z_0 = (x_0, y_0) \in \partial\Omega$, the following holds: there exists a function $\psi \in C^s([-R, R])$ with $\|\psi\|_{C^s} \le M$, taking its value in $[-PR, PR]$ and such that the membership in $\Omega$ of a point $z = (x, y)$ is equivalent to one of the two equations

$$y \le \psi(x), \quad y \ge \psi(x), \tag{8}$$

when $|x - x_0| \le R$ and $|y - y_0| \le PR$, or one of the two equations

$$x \le \psi(y), \quad x \ge \psi(y), \tag{9}$$

when $|y - y_0| \le R$ and $|x - x_0| \le PR$.

The generalization of this alternate definition to higher dimension is straightforward by considering equations one of the form

$$x_i \le \psi(x_1, \dots, x_{i-1}, x_{i+1}, \dots, x_d) \quad \text{or} \quad x_i \ge \psi(x_1, \dots, x_{i-1}, x_{i+1}, \dots, x_d), \quad i = 1, \dots, d,$$

with $\psi$ a $C^s$ function of $d-1$ variables defined on $[-R, R]^{d-1}$.

**Remark 2.4.** We stress that this definition is only valid for $s > 1$ and not for less smooth domains such as Lipschitz domains. For example if $\Omega$ is a rectangle with side oriented along principal axes, then no such local parametrization can be derived if $z_0$ is a corner point.

Consider now the family $V_n$ of oriented subgraphs from Example 3, associated with the linear space $W_n = \mathbb{P}_{n-1}$ of univariate polynomials of degree $n-1$. Let $u = \chi_\Omega$ with $\Omega$ a domain of class $C^s$ for some $s > 1$. Then, if $S$ is a $(2k+1) \times (2l+1)$ stencil centered around a cell $T$ that is crossed by the interface, we apply the above Definition 2.3 taking $z_0 = z_T$ the center of $T$. Without loss of generality, assume for example that the description of $\Omega$ near $z_0$ is by the equation

$$y \le \psi(x),$$

for $z = (x, y)$ in the rectangle $\{|x - x_0| \le R, |y - y_0| \le PR\}$. We assume that the sidelength $h$ is small enough so that the stencil $S$ is contained in this rectangle. For $z = (x, y) \in S$ we thus have

$$|x - x_0| \le C_0 h \le R, \quad C_0 = k + \frac{1}{2}.$$

Using Taylor formula and the smoothness of $\psi$ we find that there exists a polynomial $p \in \mathbb{P}_{n-1}$ such that

$$\|\psi - p\|_{L^\infty([-C_0h, C_0h])} \le C_1 h^r, \quad r := \min\{s, n\}, \tag{10}$$

where $C_1$ depends on $C_0$ and the bound $M$ on the $C^s$ norm of $\psi$. For example, we can take for $p$ the Taylor polynomials of order $\tilde{n} = \min\{\lceil s \rceil - 1, n - 1\}$ at $x = x_0$. Therefore, taking $v$, where

$$v := \chi_{\{y \le p(x)\}} \in V_n,$$

the corresponding subgraph, it follows that

$$\|u - v\|_{L^1(S)} = 2C_0hC_1h^r = Ch^{r+1}.$$

where $C$ depends on $(M, k)$. We treat the other cases $y \ge \psi(x)$, $x \le \psi(y)$ and $x \ge \psi(y)$ in a similar manner. In summary, for the local approximation error of $C^s$ domains by polynomial oriented subgraphs, we have

$$e_n(u)_S \le Ch^{r+1}, \quad r := \min\{s, n\}. \tag{11}$$

The same reasoning in $d$ dimensions delivers a local approximation estimate of order $h^{r+d-1}$.
## 2.2 Near optimal recovery from cell averages

The recovery methods that we develop in §3 and §4 are based on recovering on each singular cell $T$ an element $\tilde{u}_T \in V_n$ where $V_n$ is a given nonlinear family, based on the data of the cell-averages

$$a_S(u) = (a_{\tilde{T}}(u))_{\tilde{T} \in S},$$

where $S = S_T$ is a rectangular stencil centered around $T$. It can therefore be summarized by a local nonlinear recovery operator

$$R_T : \mathbb{R}^m \to V_n$$

where $m = (2k+1) \times (2l+1)$ is the size of the stencil, such that

$$\tilde{u}_T = R_T(a_S(u)).$$

We are interested in deriving a favorable comparison between the local recovery error $\|u - \tilde{u}_T\|_{L^1(T)}$ and the error of best approximation by $V_n$ whose magnitude can be estimated depending on the amount of smoothness of the boundary, as previously discussed.

**Definition 2.5.** The local recovery procedure is said to be near-optimal over a class of function $\mathcal{K}$ if there exists a fixed constant $C$ so that one has

$$\|u - R_T(a_S(u))\|_{L^1(T)} \le C e_n(u)_S, \quad (12)$$

for all $u$ in this class. In particular $C$ should be independent of the considered singular cell $T$ and mesh size $h$.

**Remark 2.6.** In the above definition, the recovery error on $T$ is bounded by the approximation error on the larger stencil $S$. This is due to the fact that the recovery operator $R_T$ acts on the cell averages $a_{\tilde{T}}(u)$ for all cells $\tilde{T} \subset S$.

On regular cells $T \in \mathcal{T}_h \setminus S_h$, that is, such that $a_T(u) = 0$ or $a_T(u) = 1$, we simply define the reconstruction by the constant value

$$\tilde{u}_T = a_T(u),$$

which is then the exact value of $u$. The global reconstruction of $u$ from the cell-averages $(a_T(u))_{T \in \mathcal{T}_h}$ is given by the function

$$\tilde{u} = \sum_{T \in \mathcal{T}_h} \tilde{u}_T \chi_T.$$

The global $L^1$ error can be estimated by aggregating all local error estimates, which thus gives

$$\|u - \tilde{u}\|_{L^1(D)} = \sum_{T \in S_h} \|u - \tilde{u}_T\|_{L^1(T)} \le C \sum_{T \in S_h} e_n(u)_{S_T}. \quad (13)$$

where $C$ is the stability constant in (12) and where $S_T$ denotes the stencil centered at $T$ which is used in the reconstruction of $\tilde{u}_T$.

If $\Omega$ is a $C^s$ domain with $s \ge 1$, that is, at least a Lipschitz domain, one has the cardinality estimate

$$\#(S_h) \le Ch^{-1} \quad (14)$$

where $C$ depends on the length of $\Gamma$. We may thus derive a global error estimate by combining (13), (14) and the local approximation estimates (7) and (11): is $\Omega$ a $C^s$ domain with $s \ge 1$, we obtain

$$\|u - \tilde{u}\|_{L^1(D)} \le Ch^r, \quad (15)$$

with $r := \min\{s, 2\}$ when using local recovery by linear interfaces and $r := \min\{s, n\}$ when using local recovery by polynomial subgraphs of degree $n-1$. This estimate generalizes to the higher dimensional case, combining local approximation estimates with the cardinality estimate $\#(S_h) \le Ch^{1-d}$.

Our central objective is now to propose local recovery methods that provably satisfy the near optimal recovery bound (12). For this, we start by remarking that this bound implies the property

$$R_T(a_S(v)) = v, \quad \forall v \in V_n, \quad (16)$$
<figure><img src="image_3.png" /><figcaption>Figure 3: Cases of non-injectivity: two elements of $V_n$ having same averages on a 3 × 3 stencil.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Figure 4: Illustration of the restrictions ensuring injectivity</figcaption></figure>

that is, the recovery is exact for elements from $V_n$. This property itself implies that any element from $V_n$ should be exactly characterized by its cell-average on the stencil $S$. In other words, the averaging operator

$$v \mapsto a_S(v) = (a_{\tilde{T}}(v))_{\tilde{T} \in S},$$

should be injective from $V_n$ to $\mathbb{R}^m$. For this to hold, the classes $V_n$ from Examples 1, 2, 3, 4 need to be restricted.

In the case of linear interfaces, if $H$ and $\tilde{H}$ are two half-spaces that contain the stencil $S$, the functions $v = \chi_H$ and $\tilde{v} = \chi_{\tilde{H}}$ obviously have the same cell-average vector $a_S$ with component identically equal to 1. Asking that the linear interface passes through the stencil is thus necessary but not sufficient as illustrated on Figure (3a): two lines passing only through a corner cell may result in identical cell averages. The correct restriction on $V_n$ in this case is obtained by imposing that the linear interface passes through the central cell, which in the case of a 3 × 3 stencil suffices to ensure injectivity as recalled in §4 and illustrated on Figure 4a. An important observation is that this type of restriction does not affect the local error estimates (7) since we are precisely considering a stencil $S$ centered at a cell $T$ that contains the interface. Therefore, the tangent of the interface at the point $z_0 + \psi(0)e_2$ delivering this estimate satisfies this restriction.

In the case of circular interface, asking that the disk intercepts the central cell is not sufficient as illustrated on Figure (3b): two discs $D$ and $\tilde{D}$ of equal size and contained in the central cell will result in identical cell averages for $v = \chi_D$ and $\tilde{v} = \chi_{\tilde{D}}$. In this case, an additional restriction should be that the radius of the disc is sufficiently large compared to the size of the stencil, for example by imposing that the disc center is not contained in $S$, as illustrated on Figure 4b.

In the case of subgraphs of polynomial functions, again two subgraphs passing through several cells of the stencil might have the same cell averages. This typically occurs when the polynomials are too peaky, as illustrated on Figure (3c). In this case, the additional restriction should be that the range of $p$ remains inside the stencil. By this we mean, say for a subgraph of the type $y \le p(x)$ and a rectangular stencil $S = [a, b] \times [c, d]$, one has that $p([a, b]) \subset [c, d]$, as illustrated in Figure (4c). Similarly to linear interfaces, the local error estimate (11) is not affected by such a restriction, as we discuss in §5.

A similar type of observation shows that there is no hope to uniquely characterize a piecewise linear interface from the cell averages on a given stencil if it is too peaky, that is, the opening angle of the cone embraced by the two lines cannot be arbitrarily close to 0 or $2\pi$. In other words, corners cannot be arbitrarily acute or obtuse.
# 3 Reconstruction by optimization (OBERA)

## 3.1 Presentation of the method

Optimization-Based Edge Reconstruction Algorithms (OBERA) consists in recovering on each singular cell $T \in S_h$ a recovery $\tilde{u}_T \in V_n$ by a best fit of the available cell-average data on the stencil $S = S_T$. For this purpose, we solve a minimization problem of the form

$$ \tilde{u}_T = R_T(a_S(u)) \in \arg\min_{v \in V_n} \|a_S(v) - a_S(u)\|^2 $$

where $\| \cdot \|$ is a given norm on $\mathbb{R}^m$ where $m := \#(S)$ is the size of the stencil. In a practical implementation, one first simple choice is to use the Euclidean $\ell^2$ norm, that is, minimize the loss function

$$ \mathcal{L}(u, v) := \sum_{\tilde{T} \in S} |a_{\tilde{T}}(u) - a_{\tilde{T}}(v)|^2, $$

over all $v \in V_n$. The case of linear interface corresponds to the LVIRA method [14]. Note that $v \in V_n$ is defined through an appropriate parametrization as in Example 1, 2, 3 and 4,

$$ \mu \in \mathcal{M} \subset \mathbb{R}^n \mapsto v_{\mu} \in V_n, $$

where $\mathcal{M}$ is the restricted range of the parameter $\mu$ that defines the family $V_n$. Therefore the optimization is done in practice by searching for

$$ \mu^* \in \arg\min_{\mu \in \mathcal{M}} \mathcal{L}(u, v_{\mu}) $$

and taking $\tilde{u}_T = v_{\mu^*}$.

It is interesting to note that this recovery method is not consistent in the sense that it does not guarantee that

$$ a_T(\tilde{u}_T) = a_T(u), \qquad (17) $$

a property that is required, typically in finite volume methods since it reflects the conservation of mass. In order to restore this property, one possibility is to define the recovery by solving the constrained optimization problem

$$ \tilde{u}_T \in \arg\min_{v \in V_n} \{\mathcal{L}(u, v) : a_T(v) = a_T(u)\} \qquad (18) $$

In practice, this can be emulated by modifying the loss function into

$$ \mathcal{L}(u, v) := K |a_T(u) - a_T(v)|^2 + \sum_{\tilde{T} \in S, \tilde{T} \neq T} |a_{\tilde{T}}(u) - a_{\tilde{T}}(v)|^2, \qquad (19) $$

and taking $K \gg 1$ (in our numerical tests we took $K = 100$). As explained below, this constrained recovery satisfies similar error bounds as its unconstrained counterpart.

The practical difficulty of the OBERA lies in the quick and accurate computation of the cell averages $a_{S_T}(v)$ for any given cell $T$ and $v \in V_n$, that is, have a fast evaluation procedure for the parameter to average map

$$ \mu \in \mathcal{M} \subset \mathbb{R}^n \mapsto a_{S_T}(v_{\mu}) \in \mathbb{R}^m, $$

as it is needed to calculate $\mathcal{L}(u, v_{\mu})$ at each iteration of the optimization algorithm. While analytic expressions are easily available for linear interfaces, they become more difficult to derive for more general curves like polynomials or circle interfaces. In such cases, the exact computation of $a_S(v_{\mu})$ is possible if for each cell $\tilde{T} \in S_T$ one is able to identify the points where the curved interface crosses its boundary. We followed this approach in our numerical implementation. Another option relies on quadrature methods as used in [3], however at the expense of potentially many evaluations of $v$. Finally, another perspective is to use machine learning methods in order to derive a cheaply computable surrogate of the parameter to average map.

Even with such tools in hand, the computation of $a_S(v_{\mu})$ for the many parameter values $\mu$ that are explored through the optimization process results in a time overhead that one would like to avoid. For linear interfaces, this was achieved by the ELVIRA method, as it only requires 6 calculations of $a_S(v_{\mu})$ in order to decide which $\mu^*$ should be retained. This can also be avoided for more general interfaces having higher order geometric approximations by the AEROS approach that we present in §4.
## 3.2 Analysis of the recovery error

In order to prove that the recovery error is near optimal in the $L^1(S)$ norm, we follow a general strategy introduced in [4] which is based on comparing the continuous $L^1(S)$ norm of functions and the discrete $\ell^1(\mathbb{R}^m)$ norm of their cell-averages.

In the 2d case, one obviously has on the one hand the inequality

$$h^2 \|a_S(v)\|_{\ell^1} \le \|v\|_{L^1(S)}, \quad v \in L^1(S), \qquad (20)$$

which is obtained by summing up

$$h^2 |a_{\tilde{T}}(v)| = \left| \int_{\tilde{T}} v \right| \le \|v\|_{L^1(\tilde{T})},$$

over all $\tilde{T} \in S$. This property reflects the stability of the averaging operator, between the continuous and the conveniently normalized discrete norm. Note that in more general dimension $d$ the normalizing factor is $h^d$.

Conversely, we say that the family $V_n$ satisfies an inverse stability property if there exists a constant $C$ independent of $h$ such that

$$\|v - \tilde{v}\|_{L^1(S)} \le Ch^2 \|a_S(v) - a_S(\tilde{v})\|_{\ell^1}, \quad \forall v, \tilde{v} \in V_n. \qquad (21)$$

We stress that such a property cannot hold for general pairs of integrable functions, their membership in $V_n$ is critical. Note that this property is a more quantitative version of the injectivity of the map $v \mapsto a_S(v)$ from $V_n$ to $\mathbb{R}^m$. Its validity is thus conditioned to a proper restriction of the classes $V_n$ from the various Examples 1, 2, 3, and 4, as already explained in §2.2. In the particular case of linear edges the following result was proved in [4].

**Theorem 3.1.** Let $S$ be the 3 × 3 stencil centered at $T$ and let $V_2$ be the family of linear interfaces from Example 1, with the restriction that the linear interfaces passes through $T$. Then (21) holds and the best constant is $C = \frac{3}{2}$.

The above stability and inverse stability property allow us to assess the recovery error in the following way. We first write that for any $v \in V_n$

$$\|u - \tilde{u}_T\|_{L^1(T)} \le \|u - v\|_{L^1(T)} + \|v - \tilde{u}_T\|_{L^1(T)} \le \|u - v\|_{L^1(T)} + Ch^2 \|a_S(v) - a_S(\tilde{u}_T)\|_{\ell^1},$$

where we have used (21). We then have

$$\|a_S(v) - a_S(\tilde{u}_T)\|_{\ell^1} \le \sqrt{m} \|a_S(v) - a_S(\tilde{u}_T)\|_{\ell^2} \le 2\sqrt{m} \|a_S(v) - a_S(u)\|_{\ell^2},$$

where the first inequality is Cauchy-Schwartz and the second comes by triangle inequality and the $\ell^2$ minimization property of $\tilde{u}_T$. Finally, we have

$$\|a_S(v) - a_S(u)\|_{\ell^2} \le \|a_S(v) - a_S(u)\|_{\ell^1} \le h^{-1} \|u - v\|_{L^1(S)},$$

by using (20). Combining all these, and using that $v \in V_n$ is arbitrary, we obtain that

$$\|u - \tilde{u}_T\|_{L^1(T)} \le (1 + 2C\sqrt{m})e_n(u)_S,$$

which is summarized in the following.

**Theorem 3.2.** Under (20) and (21), the recovery by $\ell^2$ minimization is near optimal in $L^1$ norm with multiplicative constant $1 + 2C\sqrt{m}$. In the case of linear interfaces using 3 × 3 stencils (LVIRA), this constant is $1 + 2\frac{3}{2}\sqrt{9} = 10$.

As observed in §2.2, near optimal local recovery allows us to derive convergence rates for smooth domains in terms of the global error estimate (15). This gives the following.

**Corollary 3.3.** If $\Omega$ is a $C^s$ domain with $s \ge 1$, the LVIRA method which is OBERA recovery based on linear interfaces converges in $L^1$ norm at rate $O(h^r)$ with $r = \min\{s, 2\}$.
<figure><img src="image_5.png" /><figcaption>Figure 5: By shifting the linear interface L one achieves average consistency on T by a linear interface L^c having the same order of accuracy.</figcaption></figure>

**Remark 3.4.** Note that if we were using a more general $ℓ^p$ norm for the data fitting, we would obtain a similar result with constant $1 + C2m^{1-\frac{1}{p}}$. The fact that we do not need to restrict ourselves just to $p = 2$ had been mentioned in [13] concerning the LVIRA method.

As previously remarked, we can modify the OBERA approach in order to impose the consistency condition (17), by solving the constrained optimization problem (18), that is optimizing inside the subfamily

$$ \tilde{V}_n := \{v \in V_n : a_T(v) = a_T(u)\} \subset V_n $$

It is obvious that if the inverse stability property (21) is valid for $V_n$, it is also valid for the smaller set $\tilde{V}_n$. We thus reach a similar estimate

$$ \|u - \tilde{u}_T\|_{L^1(T)} \le (1 + 2C\sqrt{m})\tilde{e}_n(u)_S, $$

where $\tilde{e}_n(u)_S$ is the error of best approximation of $u$ in $L^1(S)$ norm from the element of $\tilde{V}_n$, that is best approximation from $V_n$ under the consistency constraint.

We thus need to understand if $\tilde{e}_n(u)_S$ satisfies similar size estimates as $e_n(u)_S$. While this cannot be ensured in full generality (for example the set $\tilde{V}_n$ could be empty), the following simple argument shows that this indeed holds in the particular case of linear interface: the estimate (6) shows that the function $\psi$ describing the interface in the stencil $S$ satisfies

$$ a^{-} \le \psi \le a^{+}, $$

where $a^{-} = a - C_1h^r$ and $a^{+} = a + C_1h^r$ are two affine functions that parametrize two linear interfaces $L^-$ and $L^+$ that circumscribe $\Gamma$ in $S$, as illustrated on Figure 5. For the corresponding halfplanes $H^-$ and $H^+$ thus satisfy

$$ a_T(\chi_{H^-}) \le a_T(u) \quad \text{and} \quad a_T(\chi_{H^+}) \ge a_T(u). $$

Therefore, by sliding continuously a linear interface between $L^-$ and $L^+$, which corresponds to the affine function $a + t$ when $t$ varies in $[-C_1h^r, C_1h^r]$, there exists an intermediate interface $L^c$ corresponding to a particular $t^c$ and halfplane $H^c$ for which one has

$$ a_T(\chi_{H^c}) = a_T(u). $$

Therefore $v = \chi_{H^c} \in \tilde{V}_n$, and since one also has

$$ \|\psi - a^c\|_{L^\infty([-C_0h, C_0h])} \le C_1h^r, \quad r := \min\{s, 2\}, \qquad (22) $$

we reach the same estimate for $\tilde{e}_n(u)_S$ as the one obtained for $e_n(u)_S$.

From the theoretical perspective, an open problem is to establish the inverse stability bound (21) for nonlinear families $V_n$ offering higher order approximation properties than the linear interface, for which the proof of Theorem 3.1 is already quite involved. This, together with the already mentioned computational complexity of the optimization process, leads us to give up on OBERA for higher order geometrical approximation in favor of the AEROS approach that we next discuss.
<figure><img src="image_6.png" /><figcaption>Figure 6: On the left $h = 1/10$ and on the right $h = 1/30$. The singular cells identified as horizontally and vertically oriented are pictured in blue and green. For $h = 1/10$, there exists a singular cell $T$ (indicated in red) for which the interface cannot be described by a graph $y = \psi(x)$ or $x = \psi(y)$ in any stencil $S_T$ of width 3 that contains $T$. This is no more the case for $h = 1/30$ when using an adaptive selection of the stencil.</figcaption></figure>

# 4 Reconstruction on oriented stencils (AEROS)

## 4.1 Presentation of the method

Algorithms for Edge Reconstruction using Oriented Stencils (AEROS) are based on recovering an element from the family $V_n$ presented in Example 3. We thus recover on each singular cell $T \in S_h$ a domain having one of the four subgraph or epigraph forms (5), that is, an interface having one of the two Cartesian forms $y = p(x)$ or $x = p(y)$, where $p \in W_n$ a linear space of dimension $n$.

More specifically we consider for some fixed $k \ge 1$ the space

$$W_{2k+1} = \mathbb{P}_{2k},$$

of polynomials of even degree $2k$. We then use stencils $S_T$ containing $T$ of the form $(2k+1) \times L$ or $L \times (2k+1)$, when reconstructing an interface of the form $y = p(x)$ or $x = p(y)$ respectively, for some $L > 0$. As we further explain, the value of $L$ and the exact positioning of $T$ inside $S_T$ may depend on the considered cell $T$.

As a first step we need to identify for each $T \in S_h$ the exact orientation of the subgraph or epigraph that we decide to use. The decision must be based on the available data of the cell averages. We have already noticed that if $\Omega$ is a $C^s$ domain for $s > 1$, then it can itself be locally described by (at least) one of the four forms with a function $\psi \in C^s$ describing the interface, as expressed by (8) and (9) in Definition 2.3. Our objective is that our choice of form in the recovery is consistent with the form of the exact interface over the stencil $S_T$ for each cell $T$.

We thus need to identify for each $T$ an orientation $y = \psi(x)$ or $x = \psi(y)$ for the exact interface over the stencil $S_T$. Let us immediately observe that this is only possible if $h$ is below a certain critical resolution $h^* = h^*(\Omega)$ that depends on the amount of variation of the tangent to the interface, as shown on Figure 6 for the case $k = 1$ (stencils of widths 3).

More precisely, when $\Omega$ is a $C^s$ domain with $s > 1$, the variation of the slope of the tangent to the interface between two points $z$ and $z'$ is controlled by a bound of the form $M|z - z'|^r$ where $r := \min\{1, s-1\}$ and $M$ the bound on the $C^s$ norm of the functions $\psi$ that describe the interface. Therefore a given orientation, say $y \le \psi(x)$, can be maintained on a stencil $S_T$ of width $2k+1$ in the $x$ direction provided that $h \le h^* \sim k^{-1}M^{-1/r}$.

For identifying the orientation, we introduce a selection mechanism based on a numerical gradient computed by the Sobel filter. We denote by $T_\varepsilon$ with $\varepsilon = (\varepsilon_x, \varepsilon_y) \in \{-1, 0, 1\}^2$ the cells in the $3 \times 3$ stencil centered around $T = T_{0,0}$ where $\varepsilon_x$ and $\varepsilon_y$ indicate the amount of displacement by $h$ from $T$ in the $x$ and $y$ direction, respectively. We then define the numerical gradient

$$G_T = (H_T, V_T),$$
with horizontal component

$$H_T := 2a_{T_{1,0}} + a_{T_{1,1}} + a_{T_{1,-1}} - (2a_{T_{-1,0}} + a_{T_{-1,1}} + a_{T_{-1,-1}})$$

and vertical component

$$V_T := 2a_{T_{0,1}} + a_{T_{1,1}} + a_{T_{-1,1}} - (2a_{T_{0,-1}} + a_{T_{1,-1}} + a_{T_{-1,-1}})$$

where $a_{T_ε} = a_{T_ε}(u)$ are the available cell-average.

The selection mechanism is based on comparing the absolute size of $H_T$ and $V_T$ and examining their sign. More precisely:

1. If $|V_T| \ge |H_T|$ and if $V_T \le 0$, we search for a subgraph of the form $y \le p(x)$.

2. If $|V_T| \ge |H_T|$ and if $V_T > 0$, we search for an epigraph of the form $y \ge p(x)$.

3. If $|V_T| < |H_T|$ and if $H_T \le 0$, we search for a subgraph of the form $x \le p(y)$.

4. If $|V_T| < |H_T|$ and if $H_T > 0$, we search for an epigraph of the form $x \ge p(y)$.

One important result is that this selection mechanism correctly detects the orientation of the exact interface for $h$ sufficiently small, as also illustrated on Figure 6.

**Theorem 4.1.** Let $\Omega$ be a $C^s$ domain for some $s > 1$, then there exists $h^* = h^*(\Omega)$ such that under the assumption $h < h^*$, the following holds for any $T \in S_h$: in each of the above cases (1, 2, 3, 4) of the selection mechanism, the exact domain $\Omega$ can be described by an equation of the same form with $p$ replaced by a function $\psi \in C^s$ over the stencil $S_T$ centered at $T$ and of size $(2k+1) \times (2l+1)$ in case (1, 2) or $(2l+1) \times (2k+1)$ in case (2, 3) with $l = k+2$. The graph of $\psi$ remains confined in $S_T$ in the following sense: denoting by $I \times J := \cup \{\tilde{T} : \tilde{T} \in S_T\}$ the total support of $S_T$, one has $\psi(I) \subset J$ in case (1, 2) and $\psi(J) \subset I$ in case (3, 4).

We postpone the proof of this result to the Appendix, and proceed with the description and error analysis of the recovery method. We place ourselves in case 1 without loss of generality since all other cases are dealt with similarly up to an obvious exchange of $x$ and $y$ or change of sign in one of these variables.

According to the above theorem, we are ensured that $\Omega$ is characterized by the equation $y \le \psi(x)$ when $(x, y) \in S_T$ where $S_T$ is a stencil of size $(2k+1) \times (2l+1)$ with $l = k+2$. The choice $l = k+2$ is conservative and our numerical experiments show that it can sometimes be reduced while maintaining the property that the graph of $\psi$ remains confined in $S_T$. In practice we use the following adaptive strategy to use a stencil of minimal vertical side.

By convention, we denote by $(i, j)$ the coordinates of a generic cell $\tilde{T}$ when the lower left corner of $\tilde{T}$ is $(ih, jh)$. Let $(i_T, j_T)$ be the coordinates of the singular cell $T$. We explore the neighboring cells by defining

$$l^{-} := \min\{l > 0 : a_{\tilde{T}}(u) = 1, i = i_{T} - k, \dots, i_{T} + k, j = j_{T} - l - 1\},$$

which is the smallest lower shift below which we find a row of non-singular cells. Likewise, we define

$$l^{+} := \min\{l > 0 : a_{\tilde{T}}(u) = 0, i = i_{T} - k, \dots, i_{T} + k, j = j_{T} + l + 1\}.$$

Then, we take for $S_T$ the stencil that consists of cells $\tilde{T}$ of coordinates $(i, j)$ for $i = i_T - k, \dots, i_T + k$ and $j = j_T - l^{-}, \dots, j_T + l^{+}$. This stencil has size $(2k+1) \times (1+l^{-}+l^{+})$ and is centered around the cell $T$ horizontally but not vertically. From the definition of $l^{-}$ and $l^{+}$ we have the guarantee that the graph of $\psi$ remains confined in $S_T$. One option to further adapt the stencil $S_T$ is to allow that it is also not centered horizontally but still contains $T$. This leads to $(2k+1) \times (1+l^{-}+l^{+})$ stencils corresponding to values $i = i_T - k^{-}, \dots, i_T + k^{+}$ and $j = j_T - l^{-}, \dots, j_T + l^{+}$, where $k^{-}, k^{+} \ge 0$ are such that $k^{-} + k^{+} = 2k$ and are selected so to minimize the vertical size $1+l^{-}+l^{+}$.

Once the stencil $S_T$ has been selected, the polynomial $p_T \in \mathbb{P}_{2k}$ is constructed as follows. For each $i = i_T + k^{-}, \dots, i_T + k^{+}$, we denote by $R_i$ the column that consists of the cells $\tilde{T} \in S_T$ with first coordinate equal to $i$ and define the corresponding column average

$$a_i(u) = \frac{1}{l^{-} + l^{+} + 1} \sum_{\tilde{T} \in R_i} a_{\tilde{T}}(u).$$
Since the graph of $\psi$ remains confined in $S_T$, it follows that $a_i(u)$ can be identified to the univariate cell average of $\psi$ on the interval $[ih, (i+1)h]$ after having subtracted the base elevation of the stencil, that is,

$$a_i(u) + b_T = \frac{1}{h} \int_{ih}^{(i+1)h} \psi(x)dx, \quad i = i_T + k^-, \dots, i_T + k^+, \quad b_T := (j_T - l^-)h$$

We then define $p_T \in \mathbb{P}_{2k}$ as the unique polynomial of degree at most $2k$ that agrees with the observed averages of $\psi$, that is, such that

$$\frac{1}{h} \int_{ih}^{(i+1)h} p_T(x)dx = a_i(u) + b_T \quad i = i_T + k^-, \dots, i_T + k^+.$$

The polynomial $p_T$ is sometimes called the interpolant of averages, and its existence and uniqueness is standard, similar to the more usual Lagrange interpolant of point values. In particular, it is easily checked that the $\mathbb{P}_{2k}$ interpolant of the averages of a function $v$ on $2k+1$ adjacent intervals is the derivative of the $\mathbb{P}_{2k+1}$ Lagrange interpolant at the $2k+2$ interval endpoints for the primitive of $v$.

Quite remarkably, although we are locally approximating $u$ by a nonlinear family, we observe that the recovery map

$$L_T : \psi \mapsto p_T,$$

is linear, and that the AEROS recovery approach amounts to solving a simple $(2k+1) \times (2k+1)$ linear system, resulting in substantial computational saving compared to the OBERA approach.

## 4.2 Analysis of the recovery error

In order to assess the recovery error, we first observe that the above described strategy has the property of exact recovery for polynomials

$$\psi \in \mathbb{P}_{2k} \implies p_T = \psi, \qquad (23)$$

due to the uniqueness of the interpolant of averages. In other words, AEROS recovers on $T$ the true interface if it is described by a polynomial of degree $2k$ on $S_T$.

Our next observation is that the linear application $L_T$ is stable in the max norm over the relevant interval $I_T = [(i_T - k^-)h, \dots, (i_T + k^+)h]$, with stability constant that does not depend on $h$. This can be proved by making the affine change of variable

$$x = \varphi(\hat{x}) = h(i_T + \hat{x}),$$

that maps the reference interval $\hat{I} := [-k^-, \dots, k^+]$ onto $I_T$. Then, it is readily checked that

$$L\psi \circ \varphi = \hat{L}(\psi \circ \varphi),$$

where $\hat{L}$ is the average interpolant for the intervals of size 1 contained in $\hat{I}$. Therefore, we may write

$$\|L\psi\|_{L^\infty(I_T)} = \|L\psi \circ \varphi\|_{L^\infty(\hat{I})} = \|\hat{L}(\psi \circ \varphi)\|_{L^\infty(\hat{I})} \le C \|\psi \circ \varphi\|_{L^\infty(\hat{I})} = \hat{C} \|\psi\|_{L^\infty(I_T)}$$

where the constant $\hat{C}$ is the norm of $\hat{L}$ acting on $L^\infty(\hat{I})$. This constant only depends on $k$.

We are now in position to obtain an error estimate by writing for all $p \in \mathbb{P}_{2k}$,

$$\|\psi - p_T\|_{L^\infty(I_T)} \le \|\psi - p\|_{L^\infty(I_T)} + \|p_T - p\|_{L^\infty(I_T)} = \|\psi - p\|_{L^\infty(I_T)} + \|L(\psi - p)\|_{L^\infty(I_T)} \le (1 + \hat{C}) \|\psi - p\|_{L^\infty(I_T)},$$

where we have combined exact recovery of polynomials and uniform stability. Since $p$ is arbitrary we have obtained the following result.

**Theorem 4.2.** The AEROS recovery of the interface based on polynomials of degree $2k$ satisfies for each singular cell the near optimality property

$$\|\psi - p_T\|_{L^\infty(I_T)} \le (1 + \hat{C}) \min_{p \in \mathbb{P}_{2k}} \|\psi - p\|_{L^\infty(I_T)}. \qquad (24)$$
<figure><img src="image_7.png" /><figcaption>Figure 7: Circular domain and its cell average data at different scales of refinement.</figcaption></figure>

This error bound of $\psi$ in $L^\infty$ readily induces an $L^1(T)$ error bound between the recovery

$$R_T(a_S(u)) = u_T := \chi_{\{y \le p_T(x)\}}$$

and $u$ by multiplying by the width of the cell. This gives

$$\|u - R_T(a_S(u))\|_{L^1(T)} \le (1 + \hat{C})h \min_{p \in \mathbb{P}_{2k}} \|\psi - p\|_{L^\infty(I_T)} \quad (25)$$

Note that the right-side is not $e_n(u)_S$ and we thus have not obtained the near-optimal recovery property in the form (12). Nevertheless we can derive from (25) the same convergence rates since these are based on the Taylor polynomial approximation error (10), which thus yields the following result.

**Theorem 4.3.** Let $\Omega$ be a $C^s$ domain for some $s \ge 1$. The AEROS recovery of the interface based on polynomial of degree $2k$ satisfies for each singular cell $T$ a local error bound of the form

$$\|u - R_T(a_S(u))\|_{L^1(T)} \le Ch^{r+1}, \quad r := \min\{s, 2k+1\}, \quad (26)$$

and the global error bound (15) of order $O(h^r)$ for the same value of $r$.

# 5 Numerical experiments

## 5.1 Recovery of smooth domains

In this section, we compare various recovery strategies in terms of:

1. visual aspects,

2. quantitative rate of convergence,

3. computational time.

In order to draw the second comparison, we first consider the simple case of a domain $\Omega$ with circular boundary (see Figure 7) for which the recovery error can be computed within machine precision.

The following eight reconstruction strategies are compared:

1. **Piecewise Constant**: as mentioned in the introduction, the simplest linear method that one can come up with is $\tilde{u} := \sum_{T \in \mathcal{T}_h} a_T(u)\chi_T$, that is, on each cell $T$ the value given by the observed cell average $a_T(u)$.

2. **OBERA Linear**: we apply the minimization strategy described in §3 for linear interfaces as in Example 1, with loss function $\mathcal{L}(u, v) = \|a_S(u) - a_S(v)\|$ using $\ell^1$ norm and a $3 \times 3$ stencil $S$.

3. **OBERA-W Linear**: we apply the same approach but enforcing area consistency on $T$ through the weighted loss function (19) with $K = 100$.
4. **ELVIRA:** following [13], the ELVIRA method consists on replacing, still for linear interfaces, the OBERA continuous optimization strategy by providing only 6 combinations of parameters $\mu$ from which to choose $\mu^*$, the minimizer of $\mathcal{L}$ with $\ell^2$ norm instead of $\ell^1$. The 6 alternatives are obtained by proposing, for the interface's slope, the 6 possible finite differences estimations of a 3 × 3 stencil.

5. **ELVIRA-WO:** we apply ELVIRA but choosing first an *orientation*, as in AEROS, which allows to reduce the choice to 3 alternatives. In addition, we work with the modified *weighted loss function*, with $K = 100$, to favour area consistency on $T$.

6. **OBERA Quadratic:** we apply the minimization strategy for quadratic interfaces after choosing an orientation to have a Cartesian parametrization of the interface, that is, $v_\mu = \chi_P$, as in Example 3, and $p \in \mathbb{P}_2$ the space of univariate polynomials of degree 2. Here, we also use a 3 × 3 fixed stencil and enforce area consistency on $T$ through the same modified loss function (19) with $K = 100$.

7. **AEROS Quadratic:** we apply the AEROS reconstruction strategy for quadratic interfaces with the adaptive method, described in §4.1, to build stencils of width 3 minimal height.

8. **AEROS Quartic:** we apply the AEROS strategy now with polynomials of degree 4, therefore with stencils of width 5.

Figure 8 and Figure 9. display a detail of the recovery with $h = 1/10$ and $h = 1/20$ respectively in order to compare the visual quality.

As to linear interfaces, we clearly notice the relevance of enforcing area consistency on the cell $T$ of interest by appropriately modifying the loss function $\mathcal{L}$. Although the four methods benefit from similar convergence rates of $O(h^2)$ as expected from (3.3) and (4.3), see Figure 10, the reconstruction error improvement is of an order of magnitude by this simple change. When area consistency is not imposed, the linear interfaces are pushed towards the interior of the circle as the curvature of the original domain points inwardly.

As to the AEROS strategies, we notice that for $h = 1/10$, they have difficulties to reconstruct the interface on cells where $\Gamma$ is rapidly passing from a situation where an horizontal orientation is preferred to another in which a vertical one is better. This effect if particulary evident for the case of AEROS Quartic as the method needs a stencil with 5 columns. At this scale we are above the critical scale (see Figure 6) in which for some cells there is no stencil of the needed width allowing the curve to be described as a graph. This problem disappears for $h = 1/20$, for which these methods have the best visual quality.

On Figure 10, we see that in terms of convergence we obtain for both AEROS the expected rates from (4.3): for quadratics we get $O(h^3)$ and we almost get $O(h^5)$ for polynomials of degree 4. As mentioned before and graphically shown in Figure 8, AEROS Quartic breaks down when the scale of the discretization is above the critical scale which, for this particular, example is around $h = 1/20$.

Regarding the computational time per cell taken by each algorithm we observe on Table 1 that OBERA strategies are two orders of magnitude slower than any AEROS approach. At the same time, although ELVIRA methods are faster than OBERA, they still are an order slower than AEROS due to the bottleneck of having to compute the stencil cell averages to compare the 6 or 3 alternatives under evaluation. This limit is justified by the fact that choosing an orientation, as in ELVIRA-WO, cuts by half the computing time of the overall algorithm, while the improvement in accuracy is achieved by modifying the loss function.

In summary, all three comparisons in terms of visual aspect, order of convergence and computational time are in favor of the AEROS strategy provided that $h$ is below the critical scale.

Finally, we show in Figure 11 how the best linear interface method, OBERA-W Linear, and the three higher order methods compare when used to reconstruct an arbitrary, still smooth, domain. We see that by passing from the linear interface method to AEROS Quadratic the reconstruction becomes smoother while still suffering from some imperfections, notably in regions where there is a stronger change in the orientation of the curve. This is slightly improved by the optimization done in OBERA Quadratic method. An even smoother result is obtained with AEROS Quartic at the expense of some small deviations again in regions where the orientation is rapidly changing.
<figure><img src="image_8.png" /><figcaption>Figure 8: Reconstruction of a portion of the circle by different methods for a scale of $h = 1/10$.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Figure 9: Reconstruction of a portion of the circle by different methods for a scale of $h = 1/20$.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>Figure 10: Convergence for different reconstruction models. The convergence rates (in parenthesis) are estimated using values $1/h > 30$.</figcaption></figure>

<table><caption>Table 1: Average time (in seconds) taken to find the parameters of the interface by the different tested models.</caption><tr><td>OBERA Linear</td><td><b>0.8</b></td></tr><tr><td>OBERA-W Linear</td><td><b>0.7</b></td></tr><tr><td>OBERA Quadratic</td><td><b>0.3</b></td></tr><tr><td>ELVIRA</td><td><b>0.04</b></td></tr><tr><td>ELVIRA-WO</td><td><b>0.02</b></td></tr><tr><td>AEROS Quadratic</td><td><b>0.003</b></td></tr><tr><td>AEROS Quartic</td><td><b>0.003</b></td></tr></table>
<figure><img src="image_11.png" /><figcaption>Figure 11: Reconstruction of a smooth domain by different methods for a scale of $h = 1/30$.</figcaption></figure>
## 5.2 The treatment of corner domains

The previously tested strategies achieve to reconstruct, at their corresponding orders, different smooth domains with Hölder smoothness $s > 1$. This excludes the case of domains with piecewise smooth boundaries for which the presence of corners call for a specific treatment. The simplest option that we study here is to use local recovery by the approximation family $V_4$ of piecewise linear interface from Example 4. We propose, in what follows, two methods to deal with vertices of any angle, bearing in mind the limitations expressed in Figure 4c.

**AEROS Vertex:** The first method applies the AEROS strategy for $v \in \hat{V}_4$ a restriction of $V_4$ where $\pi/2 < \theta_1 < 3\pi/2$ and $-\pi/2 < \theta_2 < \pi/2$, i.e. elements $v$ whose interface can be written as an oriented graph, which is a particular instance of Example 3 where instead of searching $p$ in a space of univariate polynomials, we pick it into $W_4$ the space of piecewise functions with one breakpoint contained in $T$ or its immediate neighbours. This excludes the possibility of reconstructing a rectangle whose sides are parallel to the mesh but it applies to the case of the same rectangle slightly rotated. Under this restrictions it is possible, although lengthy, to extract explicit equations that allow us to derive a finite set of admissible parameters $\mu = (x_1, x_2, \theta_1, \theta_2)$ from the observed cell averaged vector $a_S$. We compare each proposed local approximation $v \in \hat{V}_4$ as in ELVIRA or OBERA, that is, by means of their associated loss $\mathcal{L}(u, v)$ while retaining at the end the one achieving the minimal value between the many possibilities. This same model selection strategy can be used to aggregate other competing models, like quadratic interfaces. This has the effect of keeping higher order models when the interface is locally smooth, while taking corners into account as illustrated on Figure 12. By this approach we avoid defining a vertex detection mechanism at the expense of computational overhead as we now need to compute for each cell many losses, which was already the time bottleneck for ELVIRA method.

**Tangent Extension Method (TEM):** The above restriction that the interface with a vertex needs to be an oriented graph could be limiting in some applications but it can be removed at the expense of complexifying the reconstruction procedure. Our second proposed method deals with this aspect and consists in the following steps:

1. Associate to each singular cell $T \in S_h$ some reconstruction $v_T$ stemming from any of the local interface reconstruction methods discussed so far.

2. For each cell where the presence of a vertex is suspected (eventually for all $T \in S_h$) we search for two singular cells $T_1, T_2 \in S_h$ satisfying the following

* $T \neq T_1 \neq T_2 \neq T$

* $S_{T_1} \cap \{T\} = \emptyset$

* $S_{T_2} \cap \{T\} = \emptyset$

where $S_{v_{T_i}}$ denotes the stencil $S$ used by a given smooth interface reconstruction method (for example linear or quadratic) to produce local approximations $v_{T_i}$.

3. Take the parameterized interfaces $\Gamma_1$ and $\Gamma_2$, associated to $v_{T_1}$ and $v_{T_2}$ respectively, and do an order 1 Taylor expansion at an intermediate point between cells $(T, T_1)$ and $(T, T_2)$ respectively. This yields the parameters of the two half planes $H_1$ and $H_2$ of Example 4 needed to define $v \in V_4$.

4. Finally, we compare the new local approximation $v \in V_4$ with the existing one $v_T$, as explained above, retaining only the one whose associated loss, $\mathcal{L}(u, v)$ or $\mathcal{L}(u, v_T)$, is minimal.

This procedure will reconstruct exactly corners when the interface is a line along both directions, but it will not produce area-consistent reconstructions on cell $T$ otherwise. In this regard, AEROS Vertex, being based on AEROS strategy, will yield interfaces that are, though not cell-consistent as one could get with OBERA, at least column-consistent as long as we remain in the interpolation case. This is ensured if the stencil width equals the number of parameters of the approximating class which in the case of $\hat{V}_4$ is guaranteed by using 4-width stencils.

Figure 12 displays the successive improvements in the reconstruction when combining the different strategies described so far:
* On Figure 12a, we use the AEROS Quadratic method. We observe that it recovers the interface in a satisfactory manner only far enough from corners.

* On Figure 12b, we first find a curve for each singular cell using AEROS Quadratic and then TEM. In this case, some of the problems are solved, in particular corners with a 90° angle and parallel to the grid.

* On Figure 12c, we add to the previous method the first proposed approach, based on AEROS for vertices. We obtain almost perfect results except for some cells where the quadratic approximation of the interface given by AEROS Quadratic was not replaced by a better one.

* On Figure 12d, this last issue is addressed by aggregating, before applying any of the vertex mechanisms described before, an ELVIRA-WO strategy to offer an alternative when AEROS Quadratic is too much affected by the presence of a nearby corner.
<figure><img src="image_12.png" /><figcaption>Figure 12: Reconstruction of a domain with corners from cell-averages of scale $h = 1/30$.</figcaption></figure>
## 5.3 Finite volume evolution in time

Finally, we use the interface recovery methods presented so far as a constituent part of a finite volume solver with the objective of reducing numerical dissipation. We study the particular case of a simple linear transport PDE:

$$ \frac{\partial u}{\partial t} + b \cdot \nabla u = 0, $$

in the unit square domain $D = [0, 1]^2$ with periodic boundary conditions and initial condition being a piecewise constant function $u^0: D \to \{0, 1\}$ with $\Omega$ limited by a smooth interface as in Figure 11 or an interface having corners as in Figure 12.

As in a large class finite volume schemes, the reconstruction is used at each step to compute the flux that updates the averages at the next step. For simplicity of the presentation we have set a constant (both in space and time) velocity field $b = (h/4, 0)$ and worked with unit time steps $\Delta t = 1$ and coarse grids of size $h = 1/30$, so that the CFL condition is maintained. In this case, the numerical flux induced by a local reconstruction $u_{i,j}^k$ on a cell $T$ of coordinate $(i, j)$ at time step $k$ takes the form

$$ \mathcal{F}(u_{i,j}^k) := \frac{1}{|R_T|} \int_{R_T} u_{i,j}^k(x) dx $$

where $R_T = [(i+1)h - b, (i+1)h] \times [jh, (j+1)h]$. The finite volume approximation at the next time step $k+1$ is then given by the updated cell-average

$$ a_{i,j}^{k+1} = a_{i,j}^{k} + \mathcal{F}(\tilde{u}_{i-1,j}^{k}) - \mathcal{F}(\tilde{u}_{i,j}^{k}). $$

Figure 13a displays the evolution of the $L^1$ error between the exact solution and the reconstruction, for the time evolution of a smooth domain. In this case, all three methods ELVIRA-WO, AEROS Quadratic and ELVIRA-WO + AEROS Quadratic + TEM + AEROS Vertex behave similarly with an error that is maintained at the same level for all tested times showing that numerical dissipation has been avoided on the three cases. This contrasts with the piecewise constant reconstruction that corresponds to the standard upwind scheme.

Figure 13b shows the effect of the presence of corners in the interface of $\Omega$ in terms of a slow but accumulative deterioration in both methods ELVIRA-WO, AEROS Quadratic as they are not designed to treat vertices. In contrast, ELVIRA-WO + AEROS Quadratic + TEM + AEROS Vertex, as before, keeps its error on the same level for all the times tested.

# 6 Conclusion and perspectives

In this work, we have presented several interface recovery methods. For the two main classes OBERA and AEROS, we have provided general analysis strategies for establishing convergence rates that depend on the geometric smoothness of the interface. From a practical perspective, the methods can be combined with the aggregation strategy outlined in the previous section, and we have made them available in the open-source python package $^2$.

Several natural perspectives are foreseen (and we have explored some of them already in the open-source package):

1. Address the reconstruction of more general piecewise smooth functions with jump discontinuities across geometrically smooth or piecewise smooth interfaces. This requires a proper adaptation of the interface recovery strategies, combined with a high-order treatment of the smooth part of the function corresponding to the cell $T \notin S_h$. The latter can be done by using polynomial reconstructions on stencils not containing cells of $S_h$ following the standard ENO strategy.

2. Study the use of machine learning techniques trained on sufficiently rich sets of interfaces for performing certain tasks in an automated and hopefully more efficient manner. Such tasks include, for example, the fast reconstruction of the parameter $\mu$ from the cell averages, the identification of cells that may contain vertices, or the direct access to the numerical flux in the case of finite volume scheme, as proposed for

<footnote>$^2$https://github.com/agussomacal/SubCellResolution</footnote>
<figure><img src="image_13.png" /><figcaption>Figure 13: Time evolution of the finite volume scheme error for $h = 1/30$.</figcaption></figure>

example in [3] for vertices forming angles of $\phi = 90^{\circ}$. It should however be noted that, as opposed to the approaches that we developed in this paper, the machine learning-based approach does not offer rigorous convergence guarantees. Also, our attempts to beat the AEROS method with machine learning strategies were so far unsuccessful, both from the accuracy point of view, and the runtime point of view. We provide our implementation of these strategies in the Python package.

## References

[1] R. Adams and J. Fournier *Sobolev Spaces*, Academic Press, 2003.

[2] F. Arandiga, A. Cohen, R. Donat, N. Dyn, and B. Matei *Approximation of piecewise smooth functions and images by edge-adapted (ENO-EA) nonlinear multiresolution techniques*, Applied and Computational Harmonic Analysis, 24, 225–250, 2008.

[3] B. Després and H. Jourdren, *Machine Learning design of Volume of Fluid schemes for compressible flows*, Journal of Computational Physics 408, 2020.

[4] A. Cohen, M. Dolbeault, O. Mula, A. Somacal. *Nonlinear approximation spaces for inverse problems*, Analysis and Applications (AA) 21, 01, 217–253, 2023.

[5] A. Harten, *ENO schemes with subcell resolution*, Journal of Computational Physics 83, 148-184, 1989.

[6] A. Harten, B. Engquist, S. Osher, and C. Chakravarthy, *Uniformly high order accurate essentially non-oscillatory schemes*, III, Journal of Computational Physics 71, 231-303, 1987.

[7] C. W. Hirt and B. D. Nichols. *Volume of Fluid (VOF) Method for the Dynamics of Free Boundaries*, Journal Of Computational Physics 39, 201-225 (1981).

[8] S. Ioffe and C. Szegedy. *Batch Normalization: Accelerating Deep Network Training by Reducing Internal Covariate Shift*, arXiv, 2015, http://arxiv.org/abs/1502.03167.

[9] D. Kroon. *Numerical Optimization of Kernel-Based Image Derivatives*, Short Paper University Twente, 2009.

[10] X.-D. Liu, S. Osher and T. Chan. *Weighted Essentially Non-oscillatory Schemes*, Journal of Computational Physics 115, 200–212, 1994.

[11] J. Necas, *Les méthodes directes en théorie des équations elliptiques*, Masson, 1967.
[12] W. F. Noh and P. Woodward. SLIC (Simple Line Interface Calculation), A. I. van der Vooren and P. J. Zandbergen, Eds, Lecture Notes in Physics, Springer New York, 1976, p. 330.

[13] J.E. Pilliod and E. G. Puckett, Second-order accurate volume-of-fluid algorithms for tracking material interfaces, Journal of Computational Physics 199, 465–502, 2004.

[14] E. G. Puckett, A volume-of-fluid interface tracking algorithm with applications to computing shock wave refraction, proceedings of the fourth international symposium on Computational Fluid Dynamics, 933–938, 1991.

[15] W. J. Rider and D. B. Kothe. Reconstructing Volume Tracking, Journal of Computational Physics, Vol. 141, No. 2, 1998, pp. 112-152. doi:10.1006/jcph.1998.5906.

[16] M. Rudman. Volume-Tracking Methods for Interfacial Flow Calculations, International Journal for Numerical Methods in Fluids, Vol. 24, No. 7, 1997, pp. 671-691, doi:10.1002/(SICI)1097-0363(19970415)24:7<671::AID-FLD508>3.0.CO;2-9.

[17] R. Scardovelli and S. Zaleski. Direct Numerical Simulation of Free-Surface and Interfacial Flow, Annual Review of Fluid Mechanics Vol. 31, 1999, pp. 567-603. doi:10.1146/annurev.fluid.31.1.567.

[18] D. Gueyffier, J. Li, A. Nadim, R. Scardovelli, and S. Zaleski. Volume-of-Fluid Interface Tracking with Smoothed Surface Stress Methods for Three-Dimensional Flows., Journal of Computational Physics 152, no. 2 (July 1999): 423–56. doi:10.1006/jcph.1998.6168.

[19] Virtanen, et al. SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python, Nature Methods 17, 261–272, 2020.

# A The orientation test

In this appendix we give the proof of Theorem 4.1, which is based on:

1. First studying the case where the interface $\Gamma$ is a line over the $3 \times 3$ stencil where the numerical gradient $G_T = (H_T, V_T)$ is computed.

2. Second applying a perturbation argument in the case of a general $C^s$ interface which locally deviates from a line in a quantitatively controlled manner.

## A.1 The case of a linear interface

We assume here that, over the $3 \times 3$ stencil $S$ centered at $T$, the interface $\Gamma$ is a line crossing $T$. Therefore, the restriction of $u$ to $S$ is of the form

$$u|_{S} = v_{r,\theta} := \chi_{\{\langle z - z_{T}, e_{\theta}^{\pm} \rangle \le r\}}.$$

where $z_T$ is the center of $T$ and $e_\theta = (-\sin(\theta), \cos(\theta))$ with $\theta$ the angle between $\Gamma$ and the horizontal line, that is, $e_\theta$ is the unit normal vector to $\Gamma$ pointing to the outward direction where $u|_S = 0$. The following result shows that the orientation test discriminates exactly if the direction of $\Gamma$ is closer to horizontal or vertical. Its proof uses elementary geometrical arguments, which are only sketched using pictures in order to avoid cumbersome analytic developments.

**Theorem A.1.** If $G_T = (H_T, V_T)$ is the numerical gradient based on the Sobel filter for the above function $v_{r,\theta}$, then the following holds:

* $|V_T| > |H_T|$ if and only if $\theta \in [0, \pi/4[ \cup [3\pi/4, 5\pi/4[ \cup [7\pi/4, 2\pi[$

* $|H_T| > |V_T|$ if and only if $\theta \in ]\pi/4, 3\pi/4[ \cup [5\pi/4, 7\pi/4[$

* $|H_T| = |V_T|$ if and only if $\theta \in \{\pi/4, 3\pi/4, 5\pi/4, 7\pi/4\}$

In addition
• $V_T > 0$ if and only if $\theta \in ]\pi/2, 3\pi/2[$ and $V_T < 0$ if and only if $\theta \in [0, \pi/2[ \cup ]3\pi/2, 2\pi[$

• $H_T > 0$ if and only if $\theta \in ]0, \pi[$ and $H_T < 0$ if and only if $\theta \in ]\pi, 2\pi[$.

**Proof:** Without loss of generality, we only consider the case where $\theta \in [0, \pi/4]$ since all other cases $[k\pi/4, (k+1)\pi/4]$ for $k = 1, \dots, 7$ are treated in a similar way. In order to understand the effect of $\theta$ on the values of $H_T$ and $V_T$, we parametrize the function $v_{r,\theta}$ differently: we fix $\bar{z}_T$ to be the point crossed by $\Gamma$ on the descending diagonal of $T$ (which exists and is unique when $\theta \in [0, \pi/4]$) and study $H_T$ and $V_T$ for the function

$$v_{\theta} := \chi_{\{\langle z - \bar{z}_T, e_{\theta}^{\perp} \rangle \le 0\}},$$

as we let $\theta$ vary. By scale invariance, we may assume that we work with cells of side-length equal to 1 without affecting $H_T$ and $V_T$.

Figure 14 (left) pictures the value of $V_T$ as the difference between areas of the portions of cells from the upper and lower rows crossed by the half-plane below $\Gamma$ with weight 2 for central cells and 1 for left and right cells. This difference is strictly negative for all $\theta \in [0, \pi/4]$. Its value at $\theta = 0$ is equal $-4$. As $\theta$ grows towards $\pi/4$ it first stays equal to $-4$ until it starts strictly increasing for some value $\theta^* \in [0, \pi/4[$ that depends on the position of $\bar{z}_T$ on the diagonal. This monotonic growth can be checked by observing that for $0 \le \theta_1 < \theta_2 \le \pi/4$, one has $V_T(v_{\theta_2}) - V_T(v_{\theta_1}) = V_T(v_{\theta_2} - v_{\theta_1})$ and the function $v_{\theta_2} - v_{\theta_1}$ is supported in a symmetric cone $K_{\theta_1, \theta_2}$ centered at $\bar{z}_T$ and has value 1 on the right and $-1$ on the left. Thus $V_T(v_{\theta_2}) - V_T(v_{\theta_1})$ is the sum of the areas of the portions of cells from the upper and lower rows intersected by $K_{\theta_1, \theta_2}$ with weight 2 for central cells and 1 for left and right cells, which is strictly positive if $\theta_2 > \theta^*$.

Figure 14 (center) pictures the value of $H_T$ as the difference between areas of the portions of cells from the right and left columns crossed by the half-plane below $\Gamma$ with weight 2 for central cells and 1 for lower and upper cells. This difference is null when $\theta = 0$ and increases strictly as $\theta$ grows from 0 to $\pi/4$. Once again, the strictly monotonic growth is due to the fact that $V_T(v_{\theta_2}) - V_T(v_{\theta_1})$ is the sum of the areas of the portions of cells from the left and right columns crossed by $K_{\theta_1, \theta_2}$ with weight 2 for central cells and 1 for left and right cells, which is strictly positive whenever $0 \le \theta_1 < \theta_2 \le \pi/4$.

<figure><img src="image_14.png" /><figcaption>Figure 14: Dependence of $V_T$ (left), $H_T$ (center), $|V_T| - |H_T|$ (right) as $\theta$ varies in $[0, \pi/4]$.</figcaption></figure>

This demonstrates the second statement regarding signs of $V_T$ and $H_T$.

On the other hand, by symmetry, we note that $V_T(v_\theta) = -H_T(v_\theta)$ where $\tilde{\theta} = \pi/2 - \theta$ with same base point $z_T$. Therefore $|V_T| - |H_T| = -V_T - H_T = H_T(v_\theta - v_\theta)$ and $v_\theta - v_\theta$ is supported in a symmetric cone $K_{\theta, \tilde{\theta}}$ centered at $\bar{z}_T$ with value 1 on the right and $-1$ on the left. As pictured on Figure 14 (right), this quantity is the sum of the areas of the portions of cells from the right and left column crossed by $K_{\theta, \tilde{\theta}}$ with weight 2 for central cells and 1 for lower and upper cells. These quantities are null when $\theta = \pi/4$ since the cone is then restricted to a line, and increases strictly as $\theta$ decreases from $\theta/4$ to 0 since the cone is opening.

This demonstrates the first statement regarding comparison between $|V_T|$ and $|H_T|$.

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

We will use the following direct consequence of this result, which is obtained by compactness since $V_T$ and $H_T$ are continuous with respect to $\theta$ and $\bar{z}_T$: for any $0 < \delta < \pi/4$, there exists a $\gamma = \gamma(\delta) > 0$ such that

$$\theta \in [\pi/4 + \delta, 3\pi/4 - \delta] \cup [5\pi/4 + \delta, 7\pi/4 - \delta] \implies |H_T| \ge |V_T| + \gamma \quad (27)$$

and

$$\theta \in [3\pi/4 - \delta, 5\pi/4 + \delta] \implies V_T \ge \gamma, \quad (28)$$
## A.2 A perturbation analysis

We next turn to the proof of Theorem 4.1, focusing without loss of generality on case 1. Let us assume that the interface $\Gamma = \partial\Omega$ has $C^s$ smoothness for some $s > 1$, and let $T$ be a singular cell crossed by $\Gamma$ and $S$ the $3 \times 3$ stencil centered at $T$.

We now fix $\tilde{z}_T$ to be one point of $\Gamma \cap T$ and $e_\theta^\perp = (-\sin(\theta), \cos(\theta))$ be the outer unit normal to $\Omega$ at $\tilde{z}_T$. The function

$$v_\theta := \chi_{\{(z-\tilde{z}_T, e_\theta^\perp) \le 0\}},$$

is a perturbation of $u$ as pictured on Figure 15 where the line interface $L$ is the tangent to $\Omega$ at $\tilde{z}_T$. Since $\Omega$ has $C^s$ smoothness, the deviation between the curved interface $\Gamma$ and its tangent $L$ has area of order $O(h^{s+1})$ over $S$, and therefore

$$|a_{\tilde{T}}(u) - a_{\tilde{T}}(v_\theta)| \le Ch^{s-1}, \quad \tilde{T} \in S,$$

for some fixed constant $C$ that only depends on the $C^s$ norm of the graph that locally characterizes $\Gamma$. In turn, up to multiplying the constant $C$ by 8, one has

$$|H_T(u) - H_T(v_\theta)| \le Ch^{s-1} \quad \text{and} \quad |V_T(u) - V_T(v_\theta)| \le Ch^{s-1} \qquad (29)$$

<figure><img src="image_16.png" /><figcaption>Figure 15: Approximation of a smooth interface by its tangent</figcaption></figure>

We now take any $0 < \delta < \pi/4$ arbitrarily small and consider the quantity $\gamma = \gamma(\delta)$ such that (27) and (28) are valid. For $h \le h_0$ small enough, we are ensured that

$$Ch^{s-1} \le \frac{\gamma}{3},$$

where $C$ is the constant in (29). Therefore, if $\theta \notin [0, \pi/4+\delta] \cup [7\pi/4-\delta, 2\pi[$, or equivalently $\theta \notin [-\pi/4-\delta, \pi/4+\delta]$, we obtain either by (27) that

$$|H_T(u)| \ge |H_T(v_\theta)| - \frac{\gamma}{3} \ge |V_T(v_\theta)| + \frac{2\gamma}{3} \ge |V_T(u)| + \frac{\gamma}{3} > |V_T(u)|,$$

or by (28) that

$$V_T(u) \ge V_T(v_\theta) - \frac{\gamma}{3} \ge \frac{2\gamma}{3} > 0.$$

Therefore, we conclude for case 1 that if $|H_T(u)| \le |V_T(u)|$ and $V_T(u) \le 0$, the angle $\theta$ of the tangent line $L$ necessarily lies in $[-\pi/4 - \delta, \pi/4 + \delta]$. In other words, the points $z = (x, y)$ in the half plane $\{\langle z - \tilde{z}_T, e_\theta^\perp \rangle \le 0\}$ can be characterized by an equation of the form

$$y \le a(x),$$

where $a$ is affine and $|a'(x)| \le 1 + \varepsilon$ such that $1 + \varepsilon = \tan(\pi/4 + \delta)$, with $0 \le \varepsilon \le \delta$ and $\varepsilon \sim \delta$ when $\delta$ is small. On the other hand the interface $\Gamma$ is described on $T$ by an equation of the form

$$y \le \psi(x),$$

where due to the $C^s$ smoothness of $\Gamma$, one has an estimate of the form

$$|\psi'(x) - a'(x)| \le Ch^{s-1},$$
and the same will hold on a stencil $S_T$ of width $2k + 1$ up to enlarging the value of $C$, so that

$$|\psi'(x)| \le 1 + \varepsilon + Ch^{s-1}, \quad x \in I,$$

where $I$ is the horizontal support of $S_T$. Therefore, we can find $h^* = h^*(\Omega)$ such that if $h \le h^*$ we are thus ensured that

$$|\psi'(x)| \le \frac{k+2}{k+1}, \quad x \in I.$$

This implies that the graph of $\psi$ remains confined in $S_T$ if we choose it to be of size $(2k + 1) \times (2l + 1)$ with $l = k + 2$, which concludes the proof of the theorem.