# Compact cell-based Compatible Discrete Operator diffusion scheme on Cartesian AMR mesh

Alban Vergnaud, A Lemoine, J Breil

► To cite this version:

Alban Vergnaud, A Lemoine, J Breil. Compact cell-based Compatible Discrete Operator diffusion scheme on Cartesian AMR mesh. Journal of Computational Physics, 2025, 533, 10.1016/j.jcp.2025.113982. hal-05056506

HAL Id: hal-05056506

https://normandie-univ.hal.science/hal-05056506v1

Submitted on 5 May 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

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

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

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

# Compact cell-based Compatible Discrete Operator diffusion scheme on Cartesian AMR mesh

A. Vergnaud<sup>a,b</sup>, A. Lemoine<sup>b</sup>, J. Breil<sup>a,b</sup>

## Abstract

In this paper we investigate cell-based Compatible Discrete Operator (CDO) scheme for elliptic problems. These mimetic schemes rely on mixed potential degrees of freedom at the cells and flux degrees of freedom at the faces. In this paper, we propose a compact rewriting of this scheme by eliminating a large part of the degrees of freedom associated with the faces, which greatly reduces the size of the linear systems to be inverted and therefore the computational cost. For this, we take advantage of Cartesian AMR (Adaptive Mesh Refinement) meshes, for which the discrete CDO Hodge operator can be made diagonal almost everywhere. We also propose a formulation of general Robin-type boundary conditions for these cell-based CDO schemes, also valid in the presence of Immersed Boundaries with a cut-cell approach. The proposed scheme and its performances is validated through various test cases.

**Keywords:** Compatible Discrete Operator, Mimetic scheme, Cell-based diffusion scheme, Cartesian AMR mesh, Immersed Boundaries

## 1. Introduction

Currently, numerical simulations of phenomena governed by the equations of continuum mechanics are increasingly complex. From an application and an industrial point of view, it is often necessary to use complex geometries to depict a correct physical behavior. Thereby, it is mandatory to generate meshes that faithfully represent these geometries. The meshes thus generated can present characteristics (such as small cells, bad aspect ratio, hanging nodes and/or coplanar faces) making the simulations less robust, or even unstable, for classical methods such as Finite Volume or Finite Element.

To overcome these constraints related to the generation of meshes, recent developments have been made on stable and robust numerical methods acting on arbitrary polyhedral meshes. These methods are named polyhedral discretizations methods. Examples of such methods include Mimetic Finite Difference (MFD) methods [14, 22, 35], Virtual Element method (VEM) [2, 3], Hybrid-Higher Order (HHO) method [23, 24], polyhedral finite element method [37, 42], discontinuous Galerkin (dG) method [15, 20], or Compatible Discrete Operator (CDO) methods [9, 10]. In the case of solving incompressible Navier-Stokes equations, the popular Marker-And-Cell (MAC) scheme by Harlow and Welch [32] (initially designed for Cartesian meshes with staggered grids for pressure and velocity fields) has also been extended to polyhedral meshes [6, 17, 31].

Among these polyhedral discretizations methods, many are those whose development has been driven by the concept of mimetic schemes (also called compatible schemes). This denomination refers to the fact that

<footnote><sup>a</sup>CEA CESTA, 15 Avenue des Sablières, CS 60001, 33116 Le Barp Cedex, France</footnote>

<footnote><sup>b</sup>Bordeaux INP, University of Bordeaux, CNRS, Art et Métiers Institute of Technology, INRAE, I2M Bordeaux, F-33400 Talence, France</footnote>
these schemes aim at preserving some properties of the continuous model problem at the discrete level, such as conservation laws or theorems of vector calculus. Such mimetic properties ensure intrinsically that the inferred discrete differential operators are robust and stable, even for complex polyhedral meshes (robustness and stability of the full scheme hinges on additional requirements such that reconstruction operators or time discretization properties). This type of schemes was initiated by the seminal work of Tonti [43] and first developments by Bossavit, in the field of electromagnetism [11, 12]. A key point of mimetic schemes is that the localisation of the degrees of freedom are based on the physical nature of the fields they represent: potentials are evaluated at points, gradients along lines, fluxes across surfaces and sources in volumes. These definitions allow the use of a powerful mathematical framework based on algebraic topology and differential geometry, useful to analyse and build mimetic schemes. We refer to the articles [1, 8, 9, 18, 19, 25, 26, 34, 38] for an non-exhaustive theoretical analysis and examples of mimetic schemes. Note that several classical numerical methods (Finite Element, Finite Volume or Finite Difference) can be made mimetic [38].

Recently, Bonelle [9, 10] proposed to unify several of these mimetic schemes under a single formalism and analysis framework, called Compatible Discrete Operator (CDO). The CDO framework relies on three main features: the introduction of a dual mesh, the construction of discrete differential operators and the use of discrete Hodge operators. The dual mesh is constructed so that its elements realize a one-to-one pairing with the elements of the primal mesh. Several methods can be used to build it and it plays a crucial part in other aspects of CDO schemes, such as the construction of discrete Hodge operators [9]. Through integration over mesh entities to define primal and dual degrees of freedom, it is possible to build discrete differential operators that can be exact (i.e. with no discretization error) and that preserve lots of properties of the continuous model, thanks to theorems of vector calculus (e.g. divergence theorem) [4, 9]. Finally, discrete Hodge operators are built to complete the scheme and to link degrees of freedom located on primal and dual meshes. Building discrete Hodge operators is the key point of CDO schemes as it is the only source of discretization error. There are several way to build discrete Hodge operators and to locate the degrees of freedom (on the primal or on the dual mesh), leading to a large family of CDO schemes with interesting properties. This paper does not aim to detail all these CDO schemes and the reader can refer to [9, 10] for an overview and a detailed analysis of all the schemes that can be generated. It should also be noted that recent extensions of CDO schemes to the Navier-Stokes equations have been proposed [6, 16, 39], some with connections with the MAC schemes mentioned above [6].

Among the many possible CDO schemes, schemes with degrees of freedom located at cells (also called cell-based schemes) are interesting for their similarities with MAC methods, for which the pressure degrees of freedom are evaluated at the cells of the mesh. In fact, Beltman [4] even shows an equivalence between this cell-based schemes and the MAC schemes for Cartesian meshes. This is why we decided to focus on these cell-based schemes, with a future perspective of applying them to solving the Navier-Stokes equations with projection methods. However, cell-based CDO schemes usually lead to saddle-point linear systems (due to the presence of both cells and faces degrees of freedom) which are quite expensive to solve [9], making them difficult to use in practice. One way to explicitly eliminate degrees of freedom associated to faces would be to have a diagonal discrete Hodge operator. Actually, a sufficient condition for having a local diagonal discrete Hodge operator is that primal and dual meshes are locally orthogonal and that the properties of the simulated medium are isotropic [27].

Alongside these cell-based schemes, we also want to use Cartesian AMR (Adaptive Mesh Refinement) meshes for their ability to finely and easily mesh areas of interest. This ability has also aroused the development of numerous numerical schemes specifically adapted to this type of AMR mesh. For the case of solving elliptic problems (which interest us in this paper) on AMR mesh, we can for example refer to [13, 40, 41]. These AMR meshes are relatively simple to generate and the resulting data structure is quite simple to handle. In this work, we will use the AMReX library [45] to generate these Cartesian AMR meshes and the associated data structures. It can be noticed that Cartesian AMR meshes are a particular case of polyhedral meshes (despite the presence of hanging nodes at the interfaces between two refinement levels), and so are totally usable for CDO schemes. Another advantage of this kind of meshes is related to the CDO cell-based schemes mentioned above. In fact, most of the cells are locally Cartesian (and so primal and dual meshes are locally orthogonal) which means that discrete Hodge operator is diagonal almost everywhere. It is then theoretically possible to explicitly eliminate most of the faces degrees of freedom, and so to reduce drastically
the size of the linear system. Following an idea outlined by Beltman [4, 5], the main contribution of this
paper is to propose a new and more compact formulation of cell-based CDO scheme by taking advantage of
AMR Cartesian meshes. The main difference with Beltman’s work is that we apply this idea to solve elliptic
problems using implicit linear solvers.

This paper is structured as follows. We first detail the CDO framework (discrete settings and operators) and its application to a cell-based elliptic problem. The elliptic problem is retained here for its use in both thermal and pressure equation in projection methods. General Robin boundary conditions are also presented. Then a compact rewriting of this scheme for Cartesian AMR meshes is proposed to reduce the computational cost, while maintaining the good performances of CDO schemes. Finally, this approach is validated on numerous test cases, with and without immersed boundary.

## 2. CDO framework

Before presenting the new compact formulation of the cell-based CDO scheme applied to Cartesian AMR meshes, we will first introduce the main concepts and formalism associated with CDO schemes in general: the discrete settings (meshes and degrees of freedom definitions), the discrete differential operators and the Hodge Operators. This section is not intended to be exhaustive, particularly regarding the mathematical analysis of the CDO schemes. Only the most important results are given here and the interested reader can refer to [4, 9] for more details.

### 2.1. Discrete Settings

#### 2.1.1. Primal mesh

The primal mesh is the mesh given by the mesh generator and which is used directly by practitioners. We give here a brief definition of the mesh used for CDO schemes, as well as the different associated notations.

The primal mesh $M = \{V, E, F, C\}$ is composed of a set of vertices $V$, a set of edges $E$, a set of faces $F$ and a set of cells $C$ (see Figure 1). $M$ is a discrete representation of the physical domain $\Omega$. Lowercase $x$ denotes a generic element of $X$ and $\#X$ the cardinality of the set $X$. For example, $f \in F$ refers to a given face of the mesh and $\#V$ gives the number of vertices of the mesh. We also define the subsets $A_x$ as the set of elements $a \in A$ linked to the entity $x$. For example, $E_c$ indicates the set of edges of the cell $c$ and $F_e$ is the set of faces having $e$ as an edge.

The boundary of the domain $\partial\Omega$ is discretised by the sets of vertices $V^B$, edges $E^B$ and faces $F^B$. In a complementary way, the sets of interior vertices, edges and faces are defined respectively by: $V^I = V \setminus V^B$, $E^I = E \setminus E^B$ and $F^I = F \setminus F^B$.

Generically, $|x|$ refers to the measure of an entity $x$. Thus the measure of a vertex $v$ is $|v| = 1$ by convention, $|e|$ is the length of the edge $e$, $|f|$ is the surface of the face $f$ and $|c|$ is the volume of cell $c$.

An important concept associated with meshes used in CDO formalism is the orientation of geometric entities. To each edge $e \in E$, a unit tangent vector $\tau_e$ is associated and to each face $f \in F$, a unit normal vector $\nu_f$ is associated (see Figure 1). The orientations $\tau_e$ and $\nu_f$ are chosen arbitrarily and fixed once and for all (except for a face at the boundary $f \in F^B$, for which an outward normal $\nu_f = \nu_{\partial\Omega}$ is chosen, by convention).

#### 2.1.2. Dual mesh

CDO schemes rely on the use of a dual mesh, built from the primal mesh defined previously. It has only an interest for the development and the analysis of the method but is totally implicit for a pure practitioner of CDO schemes. As for the primal mesh, we give here a brief definition of this dual mesh, as well as the associated notations. The dual mesh is defined by $\tilde{M} = \{\tilde{V}, \tilde{E}, \tilde{F}, \tilde{C}\}$, where $\tilde{V}$ collects the dual vertices $\tilde{v}$, $\tilde{E}$ the dual edges $\tilde{e}$, $\tilde{F}$ the dual faces $\tilde{f}$ and $\tilde{C}$ the dual cells $\tilde{c}$. All notations that use the tilde symbol refer to dual mesh.

There are several ways to construct the dual mesh. However, for the CDO formalism two constraints must be respected. First, the entities of the dual and primal meshes must be in bijection. More precisely,
<figure><img src="image_4.png" /><figcaption>Figure 1: Example of primal and dual meshes entities for 3D Cartesian AMR mesh.</figcaption></figure>

each vertex $v$ of the primal mesh is associated to a unique cell $\tilde{c}(v)$ of the dual mesh. Similarly, each edge $e$ of the primal mesh is associated to a unique dual face $\tilde{f}(e)$, each face $f$ of the primal mesh to a unique dual edge $\tilde{e}(f)$ and each cell $c$ of the primal mesh to a unique dual vertex $\tilde{v}(c)$. Typically, the dual mesh entities are constructed by connecting the centers of the primal mesh entities as depicted in Figure 1. An immediate consequence of this dual mesh definition is the equality of the cardinalities of the different sets: $\#\tilde{V} = \#C$, $\#\tilde{E} = \#F$, $\#\tilde{F} = \#E$ and $\#\tilde{C} = \#V$. Second, the orientations of the dual edges $\tilde{e}(f)$ and dual faces $\tilde{f}(e)$ are transferred from those of the primal mesh entities. In other words, the orientation of $\tau_{\tilde{e}(f)}$ is determined such that $\tau_{\tilde{e}(f)} \cdot \nu_f > 0$ and the orientation of $\nu_{\tilde{f}(e)}$ is determined such that $\nu_{\tilde{f}(e)} \cdot \tau_e > 0$. This orientation transfert can be seen in Figure 1.

The CDO formalism does not impose a unique way of constructing the dual mesh. Several choices are a priori possible (such as for example a construction based on a barycentric subdivision or Delaunay-Voronoï dual meshes), leading to different CDO schemes (a practical example of dual mesh building is given in section 2.1.4). For a polyhedral mesh, note that even if the edges and faces of the primal mesh are straight and plane, edges and faces of the dual mesh are not necessarily straight and plane (but piecewise straight and plane ; see for example the dual edge $\tilde{e}(f)$ depicted in Figure 1 which is piecewise straight).

The dual mesh construction described above does not provide closed geometric entities at the domain boundary. Indeed, a dual cell $\tilde{c}(v)$ such that $v \in V^B$, a dual face $\tilde{f}(e)$ such that $e \in E^B$ and a dual edge $\tilde{e}(f)$ such that $f \in F^B$ are not closed [9]. This phenomenon is illustrated in Figure 2. However, it is necessary to have closed geometric entities in order to impose the boundary conditions (see following sections). The dual mesh must therefore be extended at the boundary in order to have closed geometric entities. For this we will use the extension proposed by Beltman [4, 6], which consists of adding faces, edges and vertices to the dual mesh so that all the entities of the dual mesh thus extended are actually closed.
$$ \tilde{F} = \{ \tilde{f}_v = \tilde{c}(v) \cap \partial\Omega \mid v \in V^B \}, \quad \tilde{E} = \{ \tilde{e}_e = \tilde{f}(e) \cap \partial\Omega \mid e \in E^B \} \quad \text{and} \quad \tilde{V} = \{ \tilde{v}_f = \tilde{e}(f) \cap \partial\Omega \mid f \in F^B \} $$
denote the set of faces, edges and vertices that have been added to the dual mesh through this extension. See Figure 2 for an illustration of the different elements added to the dual mesh through this extension, as well as [4, 6] for a more detailed presentation of this point.

### 2.1.3. Regularity conditions for polyhedral meshes

It is important to question the regularity of the meshes, in the sense that CDO schemes based on these meshes had to have good properties. In particular, good convergence properties are wanted. For a given polyhedral mesh sequence $(M_n)_{n \ge 1}$ such that the size of $M_n$ tends toward 0 has $n \to \infty$, Bonelle [9] established the conditions that this sequence must respect to guarantee convergence properties of the CDO schemes. This condition can be stated as: for each mesh $M_n$ in the sequence, there exists a simplicial mesh
<figure><img src="image_5.png" /><figcaption>Figure 2: Extension of the dual mesh for a cubic cell at the boundary (hatched face belongs to $\partial\Omega$). Left: for $f \in F^B$, dual edge $\tilde{e}(f)$ is closed by adding the vertex $\bar{v}_f = \tilde{e}(f) \cap \partial\Omega$; Center: for $e \in E^B$, dual face $\tilde{f}(e)$ is closed by adding the edge $\bar{e}_e = \tilde{f}(e) \cap \partial\Omega$; Right: for $v \in V^B$, dual cell $\tilde{c}(v)$ is closed by adding the face $\tilde{f}_v = \tilde{c}(v) \cap \partial\Omega$.</figcaption></figure>

(i.e. a tetrahedral mesh in 3D and a triangular mesh in 2D) $\mathcal{G}(M_n, \tilde{M}_n)$ such that all cells, all faces and all edges of primal and dual meshes $\{M_n, \tilde{M}_n\}$ are a union of elements of $\mathcal{G}(M_n, \tilde{M}_n)$, and all aspect ratios of the cells of $\mathcal{G}(M_n, \tilde{M}_n)$ are bounded by a quantity independent of $n$.

This condition is verified for almost all meshes classically used in numerical methods. An immediate consequence is that the faces and edges of the primal and dual meshes must be piecewise planar and straight (although in practice the primal mesh has planar faces and straight edges). Also notice that this regularity condition is compatible with the presence of hanging nodes. Thus, the Cartesian AMR meshes that we will use subsequently respect this regularity condition.

### 2.1.4. Barycentric dual mesh

One of the most direct way to construct a dual mesh from any type of polyhedral mesh (including meshes containing hanging-nodes) is to use a barycentric subdivision of the primal mesh. For this, centroids of the geometric entities of the primal mesh are introduced. For an edge $e \in E$, a face $f \in F$ and a cell $c \in C$, these centroids are defined respectively by $\underline{x}_e = \frac{1}{|e|} \int_e \underline{x}$, $\underline{x}_f = \frac{1}{|f|} \int_f \underline{x}$ and $\underline{x}_c = \frac{1}{|c|} \int_c \underline{x}$. These different centroids are then connected in order to construct the associated dual mesh.

Finally, it should be noted that the primal mesh must verify a second regularity condition so that the construction of the barycentric dual mesh is possible. Indeed, each face $f \in F$ and each cell $c \in C$ must be star-shaped with respect to their centroid $\underline{x}_f$ and $\underline{x}_c$ respectively [9]. As previously, this property is verified for almost every meshes used in numerical methods.

Apart from its relative simplicity, this construction of the dual mesh is retained for the following because it notably allows a simple and direct construction of the discrete Hodge operators on each cell $c \in C$ (see section 2.3.1).

### 2.1.5. De Rham map and degrees of freedom

As mentioned in the introduction, a key point of mimetic methods in general and of CDO schemes in particular is the link between the physical nature of the different fields that are discretized and the geometric object on which the associated degrees of freedom are attached. A potential is traditionally evaluated at a point, a circulation along a line, a flux across a surface, and a source term inside a volume. So, degrees of freedom related to a potential are attached to vertices, those of a circulation to edges, those of a flux to faces, and those of a source term to cells.

Degrees of freedom in the CDO framework are defined through De Rham maps, and are simply the integration of the physical field over the geometrical entities to which they are attached. Thus, for fields of
potential p, circulation $\underline{u}$, flux $\underline{\phi}$ and source terms s, associated degrees of freedom on the primal mesh are defined as:

$$ \mathbf{p}_v = R_V(p)|_{v} = p(v), \qquad \forall v \in V, \qquad (1) $$

$$ \mathbf{u}_e = R_{\mathcal{E}}(\underline{u})|_{e} = \int_{e} \underline{u} \cdot \underline{\tau}_e, \qquad \forall e \in E, \qquad (2) $$

$$ \phi_f = R_F(\underline{\phi})|_f = \int_f \underline{\phi} \cdot \underline{\nu}_f, \qquad \forall f \in F, \qquad (3) $$

$$ s_c = R_C(s)|_c = \int_c s, \qquad \forall c \in C, \qquad (4) $$

with $V = \mathbb{R}^{\#V}$, $\mathcal{E} = \mathbb{R}^{\#E}$, $\mathcal{F} = \mathbb{R}^{\#F}$ and $C = \mathbb{R}^{\#C}$ the vector spaces containing the degrees of freedom (respectively at vertices, edges, faces and cells). So, vectors $\mathbf{p} \in V$, $\mathbf{u} \in \mathcal{E}$, $\phi \in \mathcal{F}$ and $s \in C$ contain respectively degrees of freedom at every vertices, edges, faces and cells. Note that de Rham maps act on smooth enough functions to be correctly defined [9, 44].

Of course, the definitions of De Rham maps for the primal mesh can be easily extended to the dual mesh. Thus, the components of the vectors containing the different degrees of freedom on the dual mesh are given by:

$$ \tilde{\mathbf{p}}_{\tilde{v}(c)} = R_{\tilde{V}}(p)|_{\tilde{v}(c)} = p(\tilde{v}(c)), \qquad \forall c \in C, \qquad (5) $$

$$ \tilde{\mathbf{u}}_{\tilde{e}(f)} = R_{\tilde{\mathcal{E}}}(\underline{u})|_{\tilde{e}(f)} = \int_{\tilde{e}(f)} \underline{u} \cdot \underline{\tau}_{\tilde{e}(f)}, \qquad \forall f \in F, \qquad (6) $$

$$ \tilde{\phi}_{\tilde{f}(e)} = R_{\tilde{\mathcal{F}}}(\underline{\phi})|_{\tilde{f}(e)} = \int_{\tilde{f}(e)} \underline{\phi} \cdot \underline{\nu}_{\tilde{f}(e)}, \qquad \forall e \in E, \qquad (7) $$

$$ \tilde{s}_{\tilde{c}(v)} = R_{\tilde{C}}(s)|_{\tilde{c}(v)} = \int_{\tilde{c}(v)} s, \qquad \forall v \in V. \qquad (8) $$

## 2.2. Discrete differential operators

The definitions given above for the degrees of freedom provides a natural way of defining the discrete differential operators in CDO framework, based on the theorems of vector calculus. We only present the gradient and divergence operators in the following, because they are the ones that are useful for solving an elliptic problem. The definition of the curl operator CURL, although based on the same construction principles, is not presented here as it not useful in the following of this work (see [9] for more details).

### 2.2.1. Discrete gradient

The discrete gradient operator is defined using the fundamental theorem of calculus. Indeed, the integration of a gradient along an edge $e$ (oriented by $\underline{\tau}_e$ from vertex $a$ to vertex $b$) is written as:

$$ \int_{e} \underline{\mathrm{grad}} p \cdot \underline{\tau}_{e} = p(b) - p(a) \Leftrightarrow \underline{\mathrm{grad}} u|_{e} = \mathbf{p}_{b} - \mathbf{p}_{a}, \qquad (9) $$

where gradient degrees of freedom $\underline{\mathrm{grad}} u$ attached to the edges (2) and potential degrees of freedom $\mathbf{p}$ attached to the vertices (1) appear naturally. Thus, the discrete gradient operator GRAD for the primal mesh is defined as follows:

$$ \mathrm{GRAD}: V \to \mathcal{E}, \qquad \mathrm{GRAD}(\mathbf{p})|_{e} = \sum_{v \in V_e} \iota_{v,e} \mathbf{p}_v, \qquad \forall \mathbf{p} \in V, \forall e \in E, \qquad (10) $$
with $\iota_{v,e} = +1$ if $\tau_e$ points toward $v$ and $\iota_{v,e} = -1$ otherwise. The GRAD operator can also be represented by a matrix $\mathbb{G}$ of size $\#E \times \#V$ and containing only coefficients 0, 1 or $-1$. We give an example of the construction of the matrix $\mathbb{G}$ in Figure 3.

<figure><img src="image_6.png" /><figcaption>Figure 3: Example of the discrete gradient matrix $\mathbb{G}$ for the sets of vertices $\{v_1, ..., v_8\}$ and edges $\{e_1, ..., e_{10}\}$ (arrows indicate edges arbitrary orientations).</figcaption></figure>

Similarly, the discrete gradient operator for the dual mesh is defined as follows:

$$ \widetilde{\text{GRAD}} : \tilde{\mathcal{V}} \to \tilde{\mathcal{E}}, \qquad \widetilde{\text{GRAD}}(\tilde{\mathbf{p}})|_{\tilde{\mathcal{E}}} = \sum_{\tilde{v} \in \tilde{\mathcal{V}}_{\tilde{\mathcal{E}}}} \iota_{\tilde{v}, \tilde{\mathcal{E}}} \tilde{\mathbf{p}}_{\tilde{v}}, \quad \forall \tilde{\mathbf{p}} \in \tilde{\mathcal{V}}, \forall \tilde{\mathcal{E}} \in \tilde{\mathcal{E}}, \quad (11) $$

with $\iota_{\tilde{v},\tilde{\mathcal{E}}} = +1$ if $\tau_{\tilde{\mathcal{E}}}$ points toward $\tilde{v}$ and $\iota_{\tilde{v},\tilde{\mathcal{E}}} = -1$ otherwise. $\widetilde{\text{GRAD}}$ is also represented by a matrix $\tilde{\mathbb{G}}$ of size $\#\tilde{\mathcal{E}} \times \#\tilde{\mathcal{V}} = \#F \times \#C$ and containing only coefficients 0, 1 or $-1$.

### 2.2.2. Discrete divergence

Similarly, the discrete divergence operator is defined using the divergence theorem. Indeed, the integration of a divergence over a cell $c$ is written as:

$$ \int_c \mathrm{div}(\underline{\phi}) = \int_{\partial c} \underline{\phi} \cdot \underline{\nu}_{\partial c} = \sum_{f \in F_c} \pm \int_f \underline{\phi} \cdot \underline{\nu}_f \quad \Leftrightarrow \quad \mathrm{div}\underline{\phi}|_c = \sum_{f \in F_c} \pm \phi_f, \quad (12) $$

where divergence degrees of freedom **div** attached to the cells (4) and fluxes degrees of freedom $\phi$ attached to the faces (3) appear naturally. The sign $\pm$ depends on the orientation previously given to each face $f$. Thus, the discrete divergence operator DIV for the primal mesh is defined as follows:

$$ \mathrm{DIV} : \mathcal{F} \to \mathcal{C}, \qquad \mathrm{DIV}(\phi)|_c = \sum_{f \in F_c} \iota_{f,c} \phi_f, \quad \forall \phi \in \mathcal{F}, \forall c \in \mathcal{C}, \quad (13) $$

with $\iota_{f,c} = +1$ if $\nu_f$ points outside $c$ and $\iota_{f,c} = -1$ otherwise. The DIV operator can also be represented by a matrix $\mathbb{D}$ of size $\#C \times \#F$ and containing only coefficients 0, 1 or $-1$. We give an example of the construction of the matrix $\mathbb{D}$ in Figure 4.

Similarly, the discrete divergence operator for the dual mesh is defined as follows:

$$ \widetilde{\mathrm{DIV}} : \widetilde{\mathcal{F}} \to \widetilde{\mathcal{C}}, \qquad \widetilde{\mathrm{DIV}}(\widetilde{\underline{\phi}})|_{\widetilde{\mathcal{C}}} = \sum_{\widetilde{f} \in \widetilde{F}_{\widetilde{\mathcal{C}}}} \iota_{\widetilde{f},\widetilde{\mathcal{C}}} \widetilde{\underline{\phi}}_{\widetilde{f}}, \quad \forall \widetilde{\underline{\phi}} \in \widetilde{\mathcal{F}}, \forall \widetilde{\mathcal{C}} \in \widetilde{\mathcal{C}}, \quad (14) $$

with $\iota_{\widetilde{f},\widetilde{\mathcal{C}}} = +1$ if $\nu_{\widetilde{f}}$ points outside $\widetilde{c}$ and $\iota_{\widetilde{f},\widetilde{\mathcal{C}}} = -1$ otherwise. $\widetilde{\mathrm{DIV}}$ is also represented by a matrix $\widetilde{\mathbb{D}}$ of size $\#\widetilde{\mathcal{C}} \times \#\widetilde{\mathcal{F}} = \#V \times \#E$ and containing only coefficients 0, 1 or $-1$.
<figure><img src="image_7.png" /><figcaption>Figure 4: 2D example of the discrete divergence matrix $\mathbb{D}$ for the sets of faces {$f_1, ..., f_{10}$} (arrows indicate faces arbitrary orientations) and cells {$c_1, c_2, c_3$}.</figcaption></figure>

### 2.2.3. Some properties of discrete differential operators

The definitions of discrete differential operators given above imply some interesting properties. First, the differential operators (10) and (13) are exact, in the sense that they do not involve any approximations or truncations. Thus, no numerical error is introduced by the use of these differential operators, and this is the first mimetic property of the CDO schemes. In CDO framework, for stationary and elliptic problem, the only source of numerical error will lie in the approximation of the discrete Hodge operators.

Secondly, as the operators can be represented by incidence matrices (i.e. matrices containing only 0, 1 and −1 coefficients), these discrete differential operators are therefore freed from any notion of metric. Only the relative connectivity of each entity (and therefore the general topology of the mesh) matters to define these operators. The notion of metric is only present in the Hodge operators presented in the following.

It should also be noted that the definitions given above remain perfectly valid for edges and faces of any shape (including non-straight edges or non-planar faces), because these definitions involve integrated quantities and that the theorems of vector calculus remain valid whatever the geometries considered.

Then, the differential operators verify some adjunction properties which can be summarized using matrix representation as follows: $\tilde{\mathbb{G}} = -\mathbb{D}^T$ and $\tilde{\mathbb{D}} = -\mathbb{G}^T$. These properties are due to the orientation transfer between the primal and the dual meshes, and which results in:

$$ \iota_{v,e} = -\iota_{\tilde{f}(e), \tilde{c}(v)}, \quad \forall e \in E, \forall v \in V_e, \quad (15) $$

$$ \iota_{f,c} = -\iota_{\tilde{v}(c), \tilde{c}(f)}, \quad \forall c \in C, \forall f \in F_c. \quad (16) $$

Finally, note that differential operators and De Rham maps commute in the following way: $R_{\varepsilon}(\text{grad}) = \text{GRAD}(R_{\nu})$, $R_{C}(\text{div}) = \text{DIV}(R_{\mathcal{F}})$, $R_{\tilde{\varepsilon}}(\text{grad}) = \tilde{\text{GRAD}}(R_{\tilde{\nu}})$ and $R_{\tilde{C}}(\text{div}) = \tilde{\text{DIV}}(R_{\tilde{\mathcal{F}}})$ (note that the last two relations regarding dual operators hold only if homogeneous boundary conditions are used [9]).

In addition to these properties, we can also mention the fact that the classical properties of composition between operators still hold exactly, despite the discretization. That is, $\text{CURL} \cdot \text{GRAD} = 0_{\mathcal{F}}$ and $\text{DIV} \cdot \text{CURL} = 0_{C}$ (there are similar properties for the operators defined on the dual mesh). This is another property that gives CDO schemes their mimetic character.

### 2.3. Discrete Hodge operators

An important point of the CDO framework is the concept of discrete Hodge operator [33]. These operators link the degrees of freedom which are based on entities in duality. In other words, given the degrees of freedom of a field on the primal mesh, the associated discrete Hodge operator gives the degrees of freedom of this field on the dual mesh. In a way, a discrete Hodge operator reconstructs a field on the dual mesh from its degrees of freedom on the primal mesh.
At continuous level, the fields $\underline{x}$ and $\underline{y}$, respectively evaluated on the primal and dual meshes, can be linked by a material parameter $\underline{\alpha}$ (assumed to be piecewise constant on each primal cell $c$):

$$ \underline{y} = \underline{\alpha} \underline{x}. \qquad (17) $$

This kind of relation, called closure relation in CDO formalism, typically represents a constitutive equation (for example Fourier's law for a thermal problem). The discrete Hodge operator associated with the relation (17) is then defined by:

$$ H_{\alpha}^{\chi \tilde{\mathcal{Y}}} : \mathcal{X} \to \tilde{\mathcal{Y}}, \qquad \tilde{\mathcal{Y}} = H_{\alpha}^{\chi \tilde{\mathcal{Y}}}(\underline{x}), \qquad \tilde{\mathcal{Y}} \in \tilde{\mathcal{Y}}, \underline{x} \in \mathcal{X}, \qquad (18) $$

where $\mathcal{X}$ and $\tilde{\mathcal{Y}}$ denote spaces in duality. $\underline{x}$ and $\tilde{\underline{x}}$ are the vectors containing the degrees of freedom associated with the fields $\underline{x}$ and $\tilde{\underline{y}}$.

In practice, discrete Hodge operator $H_{\alpha}^{\chi \tilde{\mathcal{Y}}}$ is obtained by a cellwise assembly of local Hodge operators $H_{\alpha}^{\chi_c \tilde{\mathcal{Y}}_c} : \mathcal{X}_c \to \tilde{\mathcal{Y}}_c$, defined at each cell $c$:

$$ H_{\alpha}^{\chi \tilde{\mathcal{Y}}} = \sum_{c \in \mathcal{C}} T_{X,c}^{*} \cdot H_{\alpha}^{\chi_c \tilde{\mathcal{Y}}_c} \cdot T_{X,c}, \qquad (19) $$

where $\mathcal{X}_c = \mathbb{R}^{\#X_c}$ and $\tilde{\mathcal{Y}}_c = \mathbb{R}^{\#\tilde{Y}_c}$ are the local restrictions to the cell $c$ of the vector spaces containing the degrees of freedom. More precisely, the local subsets $X_c$ and $\tilde{Y}_c$ are defined as:

$$ X_c = \{x \in X | x \subseteq \partial c\}, \qquad \tilde{Y}_c = \{\tilde{y}_c(x) = \tilde{y}(x) \cap c, x \in X_c\}. \qquad (20) $$

The operator $T_{X,c} : \mathbb{R}^{\#X} \to \mathbb{R}^{\#X_c}$ is the map from global to local degrees of freedom, and $T_{X,c}^*$ the adjoint operator. Thus $H_{\alpha}^{\chi \tilde{\mathcal{Y}}}$ is represented by a (large) sparse square matrix $\mathbb{H}_{\alpha}^{\chi \tilde{\mathcal{Y}}}$ of size $\#X = \#\tilde{\mathcal{Y}}$ and $H_{\alpha}^{\chi_c \tilde{\mathcal{Y}}_c}$ is represented by a (small) dense square matrix $\mathbb{H}_{\alpha}^{\chi_c \tilde{\mathcal{Y}}_c}$ of size $\#X_c = \#\tilde{\mathcal{Y}}_c$.

The design of the local Hodge operators $H_{\alpha}^{\chi_c \tilde{\mathcal{Y}}_c}$ must satisfy some properties to build a proper numerical scheme. The retained properties are [9]: a local symmetry; a local stability (i.e. the eigenvalues of the operator must be bounded and strictly positive); a local consistency (i.e. the operator must reproduce constant fields exactly). These properties are then propagated to the global discrete Hodge operator (19) by the cellwise assembly. These requirements subsequently ensure good properties of stability, convergence, consistency and robustness for the numerical schemes which will be constructed from these operators.

Unlike the differential operators defined above, the definition of the discrete Hodge operator includes a discretization error and takes into account the metric of the different meshes. It is also interesting to note that this discretization error of the numerical scheme occurs at the same level as the physical error made by modeling a constitutive equation.

Note also that only the Hodge operators acting from the primal mesh to the dual mesh are introduced here, and not from the dual mesh to the primal mesh. Indeed, such operators are more complex to design and to use in practice [9], although some reflections to circumvent these problems have been studied by Beltman [4].

### 2.3.1. Use of reconstruction operators

Several methods are used to construct local Hodge operators at each cell of the primal mesh and that fulfill the properties stated previously. In this section we present a construction of local Hodge operators based on the use of reconstruction operators. This construction has the advantage of being relatively easy to use and implement while retaining a certain generality.

The starting point is to introduce a local reconstruction operator $\underline{L}_{\mathcal{X}_c}$ on cell $c$:
$$ \underline{L}_{\mathcal{X}_c}(\mathbf{a})(\underline{x}) = \sum_{x \in \mathcal{X}_c} \mathbf{a}_x \underline{l}_{x,c}(\underline{x}), \quad \forall \mathbf{a} \in \mathcal{X}_c, \forall \underline{x} \in c, \qquad (21) $$

that gives the approximated reconstruction of a field $\underline{a}$ from its local degrees of freedom $\mathbf{a}_x$ at entities $x \in \mathcal{X}_c$ and at position $\underline{x}$. In this sense, reconstruction operator (21) acts as an approximate inverse of the de Rham map. Functions $\underline{l}_{x,c}$ are basis reconstruction functions (see below for an explicit expression of these functions). The local Hodge operator is then defined component by component as follows [9]:

$$ H_{\alpha}^{\mathcal{X}_c \bar{\mathcal{Y}}_c}|_{x', \bar{y}(x)} = \int_{c} \underline{l}_{x,c}(\underline{x}) \cdot \underline{\alpha} \underline{l}_{x',c}(\underline{x}), \quad \forall x, x' \in \mathcal{X}_c, \qquad (22) $$

where the integration is over the variable $\underline{x}$ and the functions $\underline{l}_{x,c}$ are reconstruction functions. Writing the discrete Hodge operator in this form also highlights the link between CDO formalism and Finite Element formalism [12], for which this Hodge operator can be seen as a discretisation of the bilinear form obtained by writing an elliptic problem in weak formulation. Reconstruction functions $\underline{l}_{x,c}$ can then be seen as the FE basis functions. In fact, (22) is the approximation of a bilinear form $f_{\alpha,c}$ that is usually used to solve elliptic problems in weak form (also see section 2.4):

$$ f_{\alpha,c}(\underline{a}, \underline{a}^*) = \int_{c} \underline{a} \cdot \underline{\alpha} \underline{a}^* \approx \int_{c} \underline{L}_{\mathcal{X}_c}(\underline{a}) \cdot \underline{\alpha} \underline{L}_{\mathcal{X}_c}(\underline{a}^*), \quad \forall \underline{a}, \underline{a}^* \in \mathcal{X}_c. \qquad (23) $$

(22) is then equivalent to the expression (23), but expressed component by component.

Reconstruction operators and associated reconstruction functions $\underline{l}_{x,c}$ are then designed to satisfy the properties of symmetry, stability and consistency stated above. Once again, several forms of reconstruction functions can be considered. In this work, we use the formulation proposed by Bonelle [9], based on a piecewise constant reconstruction on the barycentric subdivision of the polyhedral cell (see section 2.1.4). This formulation itself encompasses the reconstructions already proposed by Codecasa et al. [21] and Eymard et al. [27], respectively named DGA and SUSHI.

For the sake of concision we do not give here the expression of these reconstruction functions for all spaces in duality. We simply present the case that will interest us for the numerical schemes presented in the rest of this work: the operator $H_{\alpha}^{\mathcal{F}_c \bar{\mathcal{E}}_c}$ and therefore the functions $\underline{l}_{f,c}$. To define the flux reconstruction functions $\underline{l}_{f,c}$, a partition of the cell is considered: $c = \cup_{f \in F_c} \mathfrak{p}_{f,c}$ (see Figure 5 for an illustration of this partition).

<figure><img src="image_8.png" /><figcaption>Figure 5: Partition of cell c into subvolumes $\mathfrak{p}_{f,c}$ based on each faces f and on cell center.</figcaption></figure>

The local Hodge is then split into each subvolumes contribution:

$$ H_{\alpha}^{\mathcal{F}_c \bar{\mathcal{E}}_c}|_{f', \tilde{e}(f)} = \sum_{f'' \in F_c} \int_{\mathfrak{p}_{f''}, c} \underline{l}_{f,c}(\underline{x}) |_{\mathfrak{p}_{f''}, c} \cdot \underline{\alpha} \underline{l}_{f', c}(\underline{x}) |_{\mathfrak{p}_{f''}, c}, \quad \forall f, f' \in F_c, \qquad (24) $$
where the restriction of the functions on each $\mathfrak{p}_{f'',c}$ is given by:

$$ \underline{l}_{f,c}(\underline{x})|_{\mathfrak{p}_{f'',c}} = \frac{\tilde{e}_c(f)}{|c|} + \beta \left( \frac{\delta_{f,f''}}{|\mathfrak{p}_{f,c}|} \mathbb{I} - \frac{\tilde{e}_c(f'') \otimes f''}{|\mathfrak{p}_{f'',c}| |c|} \right) \tilde{e}_c(f) \quad (25) $$

$$ = \frac{\beta}{|\mathfrak{p}_{f,c}|} \tilde{e}_c(f) \delta_{f,f''} + \left( \mathbb{I} - \beta \frac{\tilde{e}_c(f'') \otimes f''}{|\mathfrak{p}_{f'',c}|} \right) \frac{\tilde{e}_c(f)}{|c|}, \quad \forall f'' \in F_c. \quad (26) $$

$\beta$ is a free parameter linked to the stabilization of the reconstruction operator, $\tilde{e}_c(f) = \int_{\tilde{e}_c(f)} \tilde{\tau}_{\tilde{e}_c(f)} = |\tilde{e}_c(f)| \tilde{\tau}_{\tilde{e}_c(f)}$ and $\underline{f} = \int_f \underline{\nu}_f = |f| \underline{\nu}_f$. The case $\beta = \frac{1}{d}$ corresponds to a DGA operator [21] while the case $\frac{1}{\sqrt{d}}$ corresponds to a SUSHI operator [27], where $d$ is the spatial dimension of the problem. The above integrals can be easily computed because the reconstruction functions are piecewise constant on each subvolume $\mathfrak{p}_{f,c}$. We will also refer to [9] for proofs that these functions indeed satisfy the properties required for discrete Hodge operators.

## 2.4. Application to cell-based elliptic problems

The formalism described above is now applied to solve an elliptic problem of unknown potential $p$:

$$ - \mathrm{div} \left( \underline{\kappa} \underline{\mathrm{grad}} \, p \right) = s, \quad (27) $$

on a domain $\Omega$ (boundary conditions are treated in a next section), with a source term $s$. The conductivity $\underline{\kappa}$ can be tensor-valued and is symmetric positive definite. This problem can be decomposed by introducing the gradient $\underline{g}$ and the flux $\underline{\phi}$ as:

$$ \underline{g} = \underline{\mathrm{grad}} \, p, \quad (28) $$

$$ \underline{\phi} = -\underline{\kappa} \underline{g} \Leftrightarrow \underline{g} = -\underline{\kappa}^{-1} \underline{\phi}, \quad (29) $$

$$ \mathrm{div} \, \underline{\phi} = s. \quad (30) $$

In the CDO formalism (28) and (30) are called topological laws, while (29) is referred as a closure relation. Topological laws are discretized using differential operators and do not introduce any discretization errors. The closure relations connect two quantities by means of a physical law (Fourier's law for our case if the potential is the temperature). These closure relations are discretized using discrete Hodge operators, which are a source of discretization error.

As pointed before, the potential solution $p$ to the elliptic problem (27) must be evaluated at points. However, we do not specify a priori whether it must be evaluated at the vertices of the primal mesh or at the vertices of the dual mesh. These two possible choices induce two families of CDO schemes [10]: vertex-based schemes (when the potential is evaluated at the vertices of the primal mesh) and cell-based schemes (when the potential is evaluated at the vertices of the dual mesh, and therefore at the cells of the primal mesh, by bijection between the two meshes). For the rest of our work we will only focus on cell-based schemes. Indeed, in order to apply the CDO framework to the resolution of Navier-Stokes equations using projection methods [30] on future works, we made this choice because there are links between CDO cell-based schemes and the MAC method (for which the pressure unknowns are located at the cells of the mesh), as already shown by Beltman [4].

To obtain the cell-based formulation of the elliptic problem, equations (28-30) are discretized. Using De Rham applications (1-8), the definition of the differential operators (13) (11), the commutative properties stated in section 2.2.3 and the discrete Hodge operator (18), the discretized problem is then algebraically written:
$$ \tilde{g} = \tilde{G}\tilde{p}, \qquad (31) $$

$$ \tilde{g} = -\mathbb{H}_{\kappa-1}^{\tilde{\mathcal{F}}\tilde{\mathcal{E}}} \phi, \qquad (32) $$

$$ \mathbb{D}\phi = s, \qquad (33) $$

where the integrated source term $s \in \mathcal{C}$ can be exact or computed by any quadrature rule (see section 4 for the consequences of this choice). This problem is a saddle-point problem: find $(\tilde{p}, \phi) \in \tilde{\mathcal{V}} \times \mathcal{F}$ such that,

$$ \left( -\mathbb{H}_{\kappa-1}^{\tilde{\mathcal{F}}\tilde{\mathcal{E}}} \quad -\tilde{G} \atop \mathbb{D} \quad 0 \right) \left( \begin{matrix} \phi \\ \tilde{p} \end{matrix} \right) = \left( \begin{matrix} 0 \\ s \end{matrix} \right). \qquad (34) $$

The minus sign in the first line is only to have a symmetric matrix, thanks to the adjunction property $\mathbb{D}^T = -\tilde{G}$. By construction of CDO operators, this system is sparse and have relatively small stencils, which gives it good properties for a possible inversion using an iterative linear solver. Indeed, for a cell c, the stencil is simply composed of the faces of the cell $f \in F_c$. For a face f, the stencil is composed of its neighboring cells $c \in C_f$ as well as the faces of these neighboring cells (and therefore in particular the face f itself).

Note that the discretization of the closure law (29) by the discrete Hodge operator can also be seen as an approximation of the bilinear form associated with the weak formulation of the problem as:

$$ -\int_c \mathrm{div} \left( \underline{\kappa} \underline{\mathrm{grad}} \, p \right) p^* = \int_c \left( \underline{\kappa} \underline{\mathrm{grad}} \, p \right) \cdot \underline{\mathrm{grad}} \, p^* = \int_c \underline{\phi} \underline{\cdot} \underline{\kappa}^{-1} \underline{\phi}^* = f_{\kappa^{-1},c} \left( \underline{\phi}, \underline{\phi}^* \right), \qquad (35) $$

where $p^*$ is a suitable test function and $\underline{\phi}^*$ its associated flux (see section 2.3.1).

### 2.4.1. Boundary conditions

In order to have a well-posed problem, boundary conditions must be added to the elliptic problem (27). We consider here the general case where these boundary conditions, imposed on the boundary $\partial\Omega$ of the computational domain, are Robin boundary conditions (for the sake of simplification, only isotropic conductivity is considered here):

$$ ap + b \underline{\mathrm{grad}} \, p \cdot \underline{\nu}_{\partial\Omega} = ap - \frac{b}{\kappa} \underline{\phi} \cdot \underline{\nu}_{\partial\Omega} = G, \qquad \text{on } \partial\Omega, \qquad (36) $$

where $a, b$ and $G$ are known functions defined on the boundary $\partial\Omega$ of the fluid domain $\Omega$. The particular cases ($a = 1, b = 0$) and ($a = 0, b = 1$) correspond to Dirichlet and Neumann boundary conditions respectively. $\underline{\nu}_{\partial\Omega}$ denotes the outward normal to the fluid domain. We present in the following a way to impose these Robin boundary conditions (36) in the system (34).

For a face $f \in F^B$ at the boundary $\partial\Omega$ of the fluid domain, flux degree of freedom $\phi_f$ (3) appears naturally by a simple integration of the expression (36):

$$ \int_f \left( ap - \frac{b}{\kappa} \underline{\phi} \cdot \underline{\nu}_f \right) = \int_f G \qquad \Rightarrow \qquad a_f |f| \bar{p}_{\bar{v}_f} - \frac{b_f}{\kappa_f} \phi_f \approx |f| G_f, \qquad (37) $$

$$ \bar{p}_{\bar{v}_f} \approx \frac{b_f}{a_f |f| \kappa_f} \phi_f + \frac{1}{a_f} G_f, \qquad (38) $$

where $\kappa_f, a_f, b_f$ and $G_f$ are respectively the values of the conductivity and the functions $a, b$ and $G$ evaluated at face $f$. As a first approximation, integrals were approximated using the values at the center of the face $f$ (the determination of boundary conditions at a higher order is left for future work). Here
$\nu_{\partial\Omega} = \nu_f$ by convention. $\tilde{p}$ designates the extension of the potential degrees of freedom at the vertices of the extension $\tilde{V}$ (cf. section 2.1.2), and therefore $\tilde{p}_{\tilde{v}_f}$ designates the potential at the vertex $\tilde{v}_f = \tilde{e}(f) \cap \partial\Omega$ (or equivalently at center of face $f$).

However, the definition of the gradient on the dual mesh (11) only involves the vertices $\tilde{v}$, and does not take into account the vertices $\tilde{v} \in \tilde{V}$ of the extension defined in section 2.1.2. For a face $f \in F^B$ at the boundary, the dual gradient can be extended by:

$$ \widetilde{\text{GRAD}}(\tilde{\boldsymbol{p}})|_{\tilde{e}(f)} = \iota_{\tilde{v}, \tilde{e}(f)} \tilde{p}_{\tilde{v}} + \tilde{p}_{\tilde{v}_f} \qquad \forall f \in F^B \quad (39) $$

$$ = \iota_{\tilde{v}, \tilde{e}(f)} \tilde{p}_{\tilde{v}} + \frac{b_f}{a_f |f| \kappa_f} \phi_f + \frac{1}{a_f} G_f \qquad \forall f \in F^B, \quad (40) $$

allowing to take into account the closure of the edge $\tilde{e}(f)$ by the vertex $\tilde{v}_f$. The sign + in front of the potential $\tilde{p}_{\tilde{v}_f}$ comes directly from the convention of an outward normal at the boundary of the fluid domain.

Thus the original saddle-point system (34) is modified in order to impose the boundary conditions:

$$ \left( -\mathbb{H}_{\kappa^{-1}}^{\tilde{\mathcal{E}}} - \mathbb{B} - \tilde{\mathbb{G}} \begin{pmatrix} \phi \\ \tilde{\boldsymbol{p}} \end{pmatrix} \right) = \begin{pmatrix} \boldsymbol{B} \\ \boldsymbol{s} \end{pmatrix}, \quad (41) $$

where the matrix $\mathbb{B}$ is a diagonal matrix with diagonal coefficients given by:

$$ \mathbb{B}_{f, \tilde{e}(f)} = \begin{cases} \frac{b_f}{a_f |f| \kappa_f} & \text{if } f \in F^B \\ 0 & \text{otherwise,} \end{cases} \quad (42) $$

and $\boldsymbol{B}$ is a vector whose coefficients are given by:

$$ \boldsymbol{B}_f = \begin{cases} \frac{1}{a_f} G_f & \text{if } f \in F^B \\ 0 & \text{otherwise.} \end{cases} \quad (43) $$

Note that for a face at the boundary, it is sufficient to multiply the associated line in the system by $a_f$ to avoid any divisions by 0. However, we keep this writing in order to maintain a symmetrical system.

### 3. Compact cell-based rewriting for Cartesian AMR meshes

In practice the saddle-point linear system (34) is too expensive to solve due to its size $\#F + \#C$, while only the potential degrees of freedom (in number $\#C$) are of interest. This system solves extra flux degrees of freedom while greatly increasing the size of the matrix to be inverted. A first solution was proposed by Bonelle [9] to reduce the size of the linear system. It consists of solving the saddle-point system using an augmented Lagrangian-Uzawa [29] algorithm, but still turns out to be relatively expensive because of the large number of iterations necessary to obtain a converged solution.

Another solution exists to reduce the size of the linear system while keeping unknowns at the cells. This is the Hybrid Cell-based CDO scheme, obtained by adding potential unknowns to the faces as Lagrange multiplier, which can be used to reduce the size of the system to $\#F$ thanks to static condensation [9]. The potential degrees of freedom at the cells of the primal mesh are determined in a post-processing step. The linear system to be inverted is of size $\#F$, which remains much more expensive to invert compared to a system ideally containing only degrees of freedom in the cells (and therefore of size $\#C$).

We propose another method to reduce the size of the linear system for cell-based schemes. As it will be seen in the following sections, this method only becomes attractive when a large part of the mesh is superadmissible (i.e. when primal and dual meshes are locally orthogonal), which is the case for AMR
Cartesian meshes that are use in this work. It is based on a rewriting of the system (34), already suggested by Beltman [4, 5]. But we propose to extend it to the resolution of elliptic problems and to cases where the linear system is inverted (in fact, Beltman does not invert the linear system and only use it in explicit time schemes for applications to the Navier-Stokes equations). Thus the rewriting proposed here no longer applies for any mesh, but nevertheless makes it possible to theoretically save computation time for specific problems, by the significant reduction in the size of the linear system.

### 3.1. Diagonal Hodge operator in the Cartesian domain

To easily invert the saddle-point system (34) it is sufficient to be able to invert the discrete Hodge operator $\mathbb{H}_{\kappa}^{\tilde{\mathcal{F}}\tilde{\mathcal{E}}}$, in order to explicitly determine the flux degrees of freedom $\phi$ and thus substitute and eliminate them. So, a sufficient condition to be able to directly invert the discrete Hodge operator would be that it is diagonal. But, for a SUSHI operator (i.e. for $\frac{1}{\sqrt{d}}$), [27] showed that a sufficient condition for having a local diagonal Hodge operator is to have an isotropic material property and a superadmissible mesh (this is the case for a locally Cartesian mesh, but also for Centroidal Voronoi Tessellation). We therefore choose to use SUSHI operators and an isotropic material property in the following.

In particular, note that for a uniform Cartesian cell (i.e. a $d$-cube without hanging-nodes) the local (and therefore diagonal) Hodge operator is simply given by:

$$ \mathbb{H}_{\alpha}^{\mathcal{F}_c \tilde{\mathcal{E}}_c} = \frac{\alpha_c}{2} \Delta x_c^{2-d} \mathbb{I}, \quad (44) $$

where $\Delta x_c$ is the size of the Cartesian cell $c$, $\alpha_c$ the material property in the cell $c$ and $\mathbb{I}$ the identity operator of size $\#F_c = 2d$. Thus, for a Cartesian AMR mesh, a simple expression is available for the discrete Hodge operator, for the majority of cells.

### 3.2. Degrees of freedom subsets in AMR Cartesian meshes

In order to distinguish areas where the local Hodge operator is diagonal from the rest of the AMR mesh, we define new subsets for the geometric entities, and therefore also for the associated degrees of freedom. The Cartesian and non-Cartesian subsets of the mesh are defined in the following way (see Figure 6 for an illustration of these definitions):

* The subset $C^C$ of Cartesian cells: the set of all $d$-cube cells that have not hanging-nodes;

* The subset $C^{NC}$ of non-Cartesian cells: the set of all cells that are not Cartesian (i.e. cut-cells and cells with hanging-nodes);

* The subset $F^{NC}$ of non-Cartesian faces: the set of faces of non-Cartesian cells (i.e. faces of which at least one neighboring cells is not Cartesian);

* The subset $F^C$ of Cartesian faces: the set of all faces that are not non-Cartesian faces (i.e. faces of which all neighboring cells are Cartesian).

These subsets satisfy: $C^C \cup C^{NC} = C$, $C^C \cap C^{NC} = \emptyset$, $F^C \cup F^{NC} = F$ and $F^C \cap F^{NC} = \emptyset$. Note that non-Cartesian cells $C^{NC}$ applied both for cells with hanging-nodes (like at the interface between two refinement level in Cartesian AMR meshes) and for cut-cells (for example in an Immersed Boundary approach). We also define $\phi^C$, the restriction of $\phi$ to Cartesian faces, and $\phi^{NC}$, the restriction of $\phi$ to non-Cartesian faces.
<figure><img src="image_9.png" /><figcaption>Figure 6: Definition of Cartesian and non-Cartesian cells and faces for an AMR Cartesian mesh with cut-cells.</figcaption></figure>

### 3.3. Compact rewriting

The partition of the mesh into Cartesian and non-Cartesian areas enables to decompose the saddle-point system (34) in the following way:

$$ \left( \begin{matrix} -\mathbb{H}_{\kappa^{-1}}^{\mathcal{F}\tilde{\mathcal{E}},C} & 0 & -\tilde{\mathbb{G}}^C \\ 0 & -\mathbb{H}_{\kappa^{-1}}^{\mathcal{F}\tilde{\mathcal{E}},NC} & -\tilde{\mathbb{G}}^{NC} \\ \mathbb{D}^C & \mathbb{D}^{NC} & 0 \end{matrix} \right) \left( \begin{matrix} \phi^C \\ \phi^{NC} \\ \tilde{\mathbf{p}} \end{matrix} \right) = \left( \begin{matrix} 0 \\ 0 \\ \mathbf{s} \end{matrix} \right), \quad (45) $$

where $\mathbb{D} = (\mathbb{D}^C \mathbb{D}^{NC})$ gives the restriction $\mathbb{D}^C$ of the divergence operator to Cartesian faces (operator of size $\#C \times \#F^C$) and the restriction $\mathbb{D}^{NC}$ of the divergence operator to non-Cartesian faces (operator of size $\#C \times \#F^{NC}$). The restrictions of the gradient operator are defined in an analogous way and adjunction properties still hold: $-\tilde{\mathbb{G}}^C = \mathbb{D}^{C^T}$ and $-\tilde{\mathbb{G}}^{NC} = \mathbb{D}^{NC^T}$. $\mathbb{H}_{\kappa^{-1}}^{\mathcal{F}\tilde{\mathcal{E}},C}$ and $\mathbb{H}_{\kappa^{-1}}^{\mathcal{F}\tilde{\mathcal{E}},NC}$ are respectively the restrictions of the global Hodge operator to Cartesian and non-Cartesian faces. Thanks to the fact that the local Hodge operator of a Cartesian cell (44) is diagonal, the assembly of the global Hodge operator (19) guarantees that $\mathbb{H}_{\kappa^{-1}}^{\mathcal{F}\tilde{\mathcal{E}},C}$ is indeed diagonal. In particular, the diagonal coefficient $h_f$ of matrix $\mathbb{H}_{\kappa^{-1}}^{\mathcal{F}\tilde{\mathcal{E}},C}$ associated with a Cartesian face $f \in F^C$ is written using (19) and (44) as:

$$ h_f = \sum_{c \in C_f} \frac{\Delta x_c^{2-d}}{2\kappa_c} = \frac{\Delta x_f^{2-d}}{2} \sum_{c \in C_f} \frac{1}{\kappa_c}, \quad (46) $$

where an harmonic mean of the conductivity appears naturally.

The restriction of the discrete Hodge operator to Cartesian faces being diagonal, one can easily invert it and therefore explicitly determine fluxes degrees of freedom at Cartesian faces: $\phi^C = -(\mathbb{H}_{\kappa^{-1}}^{\mathcal{F}\tilde{\mathcal{E}},C})^{-1} \tilde{\mathbb{G}}^C \tilde{\mathbf{p}}$. This then makes it possible to obtain a new formulation of the linear system:

$$ \left( \begin{matrix} -\mathbb{D}^C \left( \mathbb{H}_{\kappa^{-1}}^{\mathcal{F}\tilde{\mathcal{E}},C} \right)^{-1} \tilde{\mathbb{G}}^C & \mathbb{D}^{NC} \\ -\tilde{\mathbb{G}}^{NC} & -\mathbb{H}_{\kappa^{-1}}^{\mathcal{F}\tilde{\mathcal{E}},NC} \end{matrix} \right) \left( \begin{matrix} \tilde{\mathbf{p}} \\ \phi^{NC} \end{matrix} \right) = \left( \begin{matrix} \mathbf{s} \\ 0 \end{matrix} \right). \quad (47) $$
This system, always symmetric by construction, is therefore of size #C + #F$^{NC}$. We have therefore reduced the size of the linear system to be solved. When using a Cartesian AMR mesh, this reduction is even more significant if the characteristic size of the mesh tends towards 0 (for a fixed number of refinement levels). In this case, #C ≫ #F$^{NC}$ and the size of the system tends towards #C, which is a huge gain compared to the point-saddle system (34) or compared to the Hybrid Cell-based scheme.

It should also be noted that the system (47) remains essentially sparse, although the shape of the stencils has been modified. Indeed, for a Cartesian cell $c \in C^C$, the stencil is composed of the non-Cartesian faces $f \in F_c \cap F^{NC}$ of the cell c and the neighboring Cartesian cells of c (including c itself). For a non-Cartesian cell $c \in C^{NC}$, the stencil is composed of the faces $f \in F_c$. These faces are necessarily non-Cartesian, by definition of the sets introduced in section 3.2. Note that the non-Cartesian cell is not part of the stencil and that therefore the associated diagonal coefficient is zero. Finally, for a non-Cartesian face $f \in F^{NC}$, the stencil is composed of its neighboring cells $c \in C_f$ and the non-Cartesian faces of these neighboring cells (and therefore in particular the face f itself). So relatively small stencils are maintained, keeping good properties for a possible inversion using an iterative linear solver.

For the sake of concision, we do not present here the imposition of boundary conditions in the rewritten system (47). To do this, one must apply the decomposition and rewriting presented in this section directly to the saddle-point system containing the boundary conditions (41), which is quite straightforward given what has been presented here.

## 4. Numerical results

This section is devoted to applying the proposed scheme (47) to practical test cases, on AMR Cartesian meshes. To handle such meshes and the associated data structure, we use the AMReX library [45].

In the following sections, errors over a potential field $p$ are given by using a $L^2$ norm defined as:

$$Er_{L^2}(\tilde{p}) = \left( \frac{\sum_{c \in C} (\tilde{p}_{\tilde{v}(c)} - p|_{\tilde{v}(c)})^2 |c|}{\sum_{c \in C} (p|_{\tilde{v}(c)})^2 |c|} \right)^{1/2} \quad (48)$$

On few cases, error over a theoretical flux field $\phi$ is also studied and defined similarly as :

$$Er_{L^2}(\phi) = \left( \frac{\sum_{f \in F} (\phi_f - R_F(\phi)|_f)^2}{\sum_{f \in F} (R_F(\phi)|_f)^2} \right)^{1/2} \quad (49)$$

CDO linear systems are solved using an iterative solver. We chose to use a biconjugate gradient stabilized method (BiCGSTAB), together with an incomplete LU factorization preconditioner (more precisely an ILU(0) preconditioner). In practice, we used the Hypre [7, 28] library to perform this resolution and preconditioning. The tolerance of the solver is set to $\epsilon = 10^{-10}$. The influence of the solver and the preconditioner on the performance of this new method is not studied in this paper and left for future studies.

To estimate the performances of the new scheme (47), we use the estimation of the computational cost defined by Bonelle [9]:

$$\chi = NNZ \cdot n_{ite}, \quad (50)$$

with NNZ the number of non-zeros coefficients in the system matrix and $n_{ite}$ the number of iterations made by the iterative solver to reach a residual below the required tolerance. This definition is relevant to estimate the computational cost because the most costly operation of the iterative solver is the matrix-vector product, which cost $O(NNZ)$ at each iteration.
An important property of elliptic problems (27) is the maximum principle and we want to check if our discrete model (47) fulfills this principle. This principle states that, if $s \le 0$, then the potential solution $p$ has its maximum on the boundary (and if $s \ge 0$, then the potential solution $p$ has its minimum on the boundary) [36]. For all cases presented in this paper, we observed that this principle is respected in its discretized version. This result is only a numerical observation and does not constitute a formal proof.

## 4.1. Preliminary results

A key property that we want to verify first of all is the behavior of the proposed scheme for constant and/or linear potential solutions. For this, different mesh configurations (different number of refinement levels and/or presence of immersed boundaries) and boundary conditions (Dirichlet and/or Neumann and/or Robin) are tested. For all the tested configurations, we observed an error given by machine precision. This result follows from the linear consistency of the discrete Hodge operator [10] and has already been observed numerically in [9] for other CDO schemes.

Another result that needs to be checked is the comparison between the new compact cell-based scheme (47) and the classical cell-based scheme (34), both using SUSHI reconstruction operator. Although the linear systems associated with these two formulations are a priori different, identical results are observed, up to machine precision and linear solver tolerance $\epsilon$. This result is consistent with the fact that the compact scheme is a simple rewriting of the classical one.

## 4.2. Periodic test case

First, a periodic test case is considered. The chosen analytical potential solution is given over a unit square (or cube) computational domain by:

$$p(x) = \begin{cases} 1 + \cos(2\pi x) \cos(2\pi y) & , \text{ in } 2\text{D} \\ 1 + \cos(2\pi x) \cos(2\pi y) \cos(2\pi z) & , \text{ in } 3\text{D} \end{cases} \qquad (51)$$

The theoretical source term is computed according to (27) and no specific boundary conditions are imposed since this case is periodic. The solution is then defined up to a constant, and so we imposed the value at an arbitrary cell (here the bottom left cell) to have a well-posed problem. For this case, we investigate two ways of computing the numerical source term $s$. The integration (4) can be computed exactly or approximated by the value at the center of the cell, weighted by its volume. These two ways of computing the source term are studied below, for this specific test case. For all refinement boxes, their widths are half the width of the direct coarser box. An example of the computed solution is given in Figure 7.

In Figures 8 and 9, we plot error on potential field for 2D and 3D cases, for different number of refinement levels and for the two ways of computing the source term mentioned above. For all cases, a 2^nd^ order convergence is observed, even with an important number of refinement levels. These results also highlight the influence of the computation of the source term. In particular, more accurate results are obtained with an exact evaluation of the source term, which is expected. **One can note that increasing the number of refinement levels leads to a nearly constant error, for a given cell size at coarsest level. We assume that this is because most of the measured error is due to the coarsest level (which does not change when the number of refinement levels increases).** Notice that for a fully 2D Cartesian mesh, the solution obtained with an exact source term reached machine precision.

We now investigate the proportion of additional non-Cartesian face degrees of freedom that remain in the new linear system (47). In fact, from a computational point of view it is important that this proportion remain as low as possible to not increase the cost of the inversion of (47). Results are given in Figure 10. For all case, additional face degrees of freedom remained relatively low compared to the total number of degrees of freedom. Moreover, it is clear that this proportion tends towards zero as the number of cells in primal mesh increases. Thus, the size of the linear system to invert remains reasonable, ensuring reasonable computational cost compared to methods where only cell degrees of freedom are used.
<figure><img src="image_10.png" /><figcaption>Figure 7: Periodic test case: Computed solution with 3 refinement levels, in 2D (left) and 3D (right).</figcaption></figure>

### 4.3. Localize test case

The second test case consists in a peaked potential field given by :

$$p(x) = 1 + \exp\left(-\frac{\|x - x_c\|^2}{R^2}\right), \qquad (52)$$

with $x_c$ the position of the center of the peak and $R$ its characteristic size. The analytical source term is determined according to (27) and its discretized version $s$ is here simply evaluated as the value at the center of the cell, weighted by its volume. This solution decreases rapidly when going away from the peak center. Therefore, approximate uniform Dirichlet boundary conditions are sufficient for this case as machine precision is reached quickly. Here we chose $R = 0.025$ and a peak centered in a square domain of size $L = 1$.

Moreover, a lot of precision is needed around the potential peak while the rest of the domain can remain coarsely meshed. So this kind of test case is particularly interesting to highlight the interest of such an AMR capable scheme. In this case, the width of the finest refinement level is fixed to $5R$ and the coarser refinement level are determined so that a fixed number of cell lies between two refinement levels (here 8 cells are chosen). This ensures that maximum precision is given around the peak, without making unnecessary computations far from the area of interest. Therefore, no rate of convergence can be computed here because the position of the refinement interfaces changes for each mesh a priori. Here, we are only interested in relative error measurements. A rigorous convergence study can be performed for this test case as for the previous one. For the sake of concision, the results of this convergence study are not shown here but it should be noted that results similar to the previous case are obtained, with a second convergence rate observed. A visualization of the computed potential field is given in Figure 11.

In Figure 12, error on potential fields are plotted with respect to the estimated computational cost $\chi$ (50). These results show that it is always beneficial to increase the number of refinement level, with a view to optimizing cost/precision. So, even if the number of refinement level increases, computational costs do not explode and good performances are maintained. This shows the capability of the new formulation to deal with Cartesian AMR meshes.

### 4.4. Robin boundary conditions test case

We now investigate the performances of the proposed scheme with more complex Robin boundary conditions, as depicted in Figure 13. This case is more physical and consists in a square domain of length $L$,
<figure><img src="image_11.png" /><figcaption>Figure 8: Periodic test case: Error on 2D potential field for different number of refinement level and for mean and exact source term.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Figure 9: Periodic test case: Error on 3D potential field for different number of refinement level and for mean and exact source term.</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>Figure 10: Periodic test case: proportion of additional non-Cartesian face degrees of freedom (#$F^{NC}$) to total degrees of freedom (#$C$ + #$F^{NC}$), for different number of refinement level and in 2D and 3D.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Figure 11: Localize test case: Computed solution with 2 refinement levels, in 2D (left) and 3D (right).</figcaption></figure>

heated to temperature $T_0$ at $y = 0$ and $y = L$ and cooled by convection at $x = 0$ and $x = L$. No source term is used here. Potential $p$ corresponds here to a temperature field $T$. This problem is governed by the Biot number $Bi = \frac{hL}{\kappa}$, where $h$ is the heat transfer coefficient. $T_\infty$ is the temperature outside the computational domain.

The continuous solution to this case can be written as a function series:
<figure><img src="image_15.png" /><figcaption>Figure 12: Localize test case: Error on 2D and 3D potential fields for different numbers of refinement level, with respect to an estimation of the computational cost $\chi$ (50).</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>Figure 13: 2D test case with Robin boundary conditions.</figcaption></figure>

$$T(x, y) = T_0 + \sum_{k \in \mathbb{N}^{odd}} A_k \cosh \left( \lambda_k \left( x - \frac{L}{2} \right) \right) \cos \left( \lambda_k \left( y - \frac{L}{2} \right) \right), \quad (53)$$

$$\lambda_k = \frac{k\pi}{L}, \quad (54)$$

$$A_k = \frac{4Bi(T_\infty - T_0)(-1)^{\frac{k-1}{2}}}{k\pi(k\pi\sinh(\frac{k\pi}{2}) + B_{21}\cosh(\frac{k\pi}{2}))}. \quad (55)$$
This test case can be further restricted to $(x, y) \in [0, \frac{L}{2}]^2$, with zero Neumann boundary conditions for $x = \frac{L}{2}$ and $y = \frac{L}{2}$, thanks to the symmetries of the problem. L is set to $L = 1$ here, temperatures are chosen to $T_0 = 1.5$ and $T_\infty = 0$, and several Biot numbers are investigated. For practical reasons, only 200 terms are kept in the sum (53) to achieve sufficient convergence. Refinement areas are placed at the corners of the computational domain because this is where gradients are stronger (we also investigate refinement areas at the center of the computational domain, but we obtained similar convergence rate although increasing the error of the solution). Figure 14 shows the computed solution for two different Biot numbers, highlighting that the solution is stiffer with an important Biot number.

<figure><img src="image_17.png" /><figcaption>Figure 14: Robin boundary conditions test case: Computed solution for two different Biot numbers, $Bi = 2$ (left) and $Bi = 200$ (right) (computational domain restricted to $(x, y) \in [0, \frac{L}{2}]^2$ thanks to the symmetries of the problem).</figcaption></figure>

In Table 1 we give errors computed over the solution and the associated convergence rates, for various numbers of refinement levels and Biot numbers. A theoretical second order convergence rate is reached for the smoothest case (here for $Bi = 2$). For greater Biot numbers (and so for less smooth solution) the convergence rate is between one and two. This is easily explained by the fact that the analytical solution (53) then tends towards a discontinuous solution (at the corners of the domain), thereby limiting the order of the scheme. The number of refinement levels seems to have a lower influence on the convergence rate of the method. However, these results show that for a given mesh size (in other words for a given $\Delta x_0$) results are more accurate when the number of refinement levels increase. This shows again the interest of Cartesian AMR meshes.

### 4.5. Immersed Boundary test case

Finally, the proposed CDO scheme is tested with Immersed Boundary, with a cut-cell approach as described in section 3.2. For the sake of simplification, cases are limited here to 2D cases with a circle as an immersed boundary. Polygonal cut-cells are defined using the intersection points between the faces and the analytical shape of the circle. The geometric characteristics, the associated dual mesh and the local Hodge operator are then easily determined using the constructions presented above.

Here the length of the square domain and the radius of the immersed circle are respectively set to $L = 1$ and $R = 0.1$. The circle is centered in the domain and refinement levels are placed in the same way as for the periodic case.

A parabolic solution $p(x, y) = 1 - (x - \frac{1}{2})^2 - (y - \frac{1}{2})^2$ and a sinusoidal solution $p(x, y) = 1 + \sin(2\pi x)\sin(2\pi y)$ are studied. For the parabolic case, uniform Dirichlet boundary conditions are set on the circle and uniform
<table><caption>Table 1: Robin boundary conditions test case: Convergence results for various Biot numbers and numbers of refinement levels.</caption><thead><tr><th rowspan="2"># level</th><th rowspan="2">Δx<sub>0</sub></th><th rowspan="2">#C</th><th rowspan="2">#F<sup>NC</sup></th><th colspan="2">Bi=2</th><th colspan="2">Bi=20</th><th colspan="2">Bi=200</th></tr><tr><th>Error</th><th>Rate</th><th>Error</th><th>Rate</th><th>Error</th><th>Rate</th></tr></thead><tbody><tr><td rowspan="4">0</td><td>6.25e-02</td><td>256</td><td>0</td><td>5.61e-04</td><td>-</td><td>3.06e-03</td><td>-</td><td>4.19e-03</td><td>-</td></tr><tr><td>3.13e-02</td><td>1024</td><td>0</td><td>1.58e-04</td><td><strong>1.83</strong></td><td>1.09e-03</td><td><strong>1.50</strong></td><td>2.19e-03</td><td><strong>0.94</strong></td></tr><tr><td>1.56e-02</td><td>4096</td><td>0</td><td>4.37e-05</td><td><strong>1.86</strong></td><td>3.49e-04</td><td><strong>1.64</strong></td><td>1.06e-03</td><td><strong>1.05</strong></td></tr><tr><td>7.81e-03</td><td>16384</td><td>0</td><td>1.19e-05</td><td><strong>1.87</strong></td><td>1.06e-04</td><td><strong>1.72</strong></td><td>4.51e-04</td><td><strong>1.23</strong></td></tr><tr><td rowspan="4">1</td><td>6.25e-02</td><td>448</td><td>66</td><td>2.47e-04</td><td>-</td><td>1.11e-03</td><td>-</td><td>2.21e-03</td><td>-</td></tr><tr><td>3.13e-02</td><td>1792</td><td>130</td><td>6.41e-05</td><td><strong>1.95</strong></td><td>3.54e-04</td><td><strong>1.65</strong></td><td>1.06e-03</td><td><strong>1.06</strong></td></tr><tr><td>1.56e-02</td><td>7168</td><td>258</td><td>1.66e-05</td><td><strong>1.95</strong></td><td>1.07e-04</td><td><strong>1.73</strong></td><td>4.51e-04</td><td><strong>1.23</strong></td></tr><tr><td>7.81e-03</td><td>28672</td><td>514</td><td>4.48e-06</td><td><strong>1.89</strong></td><td>3.44e-05</td><td><strong>1.63</strong></td><td>2.18e-04</td><td><strong>1.05</strong></td></tr><tr><td rowspan="4">2</td><td>6.25e-02</td><td>640</td><td>132</td><td>2.04e-04</td><td>-</td><td>4.28e-04</td><td>-</td><td>1.10e-03</td><td>-</td></tr><tr><td>3.13e-02</td><td>2560</td><td>260</td><td>5.05e-05</td><td><strong>2.01</strong></td><td>1.21e-04</td><td><strong>1.82</strong></td><td>4.56e-04</td><td><strong>1.27</strong></td></tr><tr><td>1.56e-02</td><td>10240</td><td>516</td><td>1.26e-05</td><td><strong>2.00</strong></td><td>3.72e-05</td><td><strong>1.71</strong></td><td>2.18e-04</td><td><strong>1.06</strong></td></tr><tr><td>7.81e-03</td><td>40960</td><td>1028</td><td>3.21e-06</td><td><strong>1.97</strong></td><td>1.59e-05</td><td><strong>1.23</strong></td><td>1.34e-04</td><td><strong>0.70</strong></td></tr><tr><td rowspan="4">3</td><td>6.25e-02</td><td>832</td><td>198</td><td>2.00e-04</td><td>-</td><td>2.71e-04</td><td>-</td><td>5.33e-04</td><td>-</td></tr><tr><td>3.13e-02</td><td>3328</td><td>390</td><td>4.93e-05</td><td><strong>2.02</strong></td><td>6.95e-05</td><td><strong>1.97</strong></td><td>2.28e-04</td><td><strong>1.22</strong></td></tr><tr><td>1.56e-02</td><td>13312</td><td>774</td><td>1.23e-05</td><td><strong>2.01</strong></td><td>2.13e-05</td><td><strong>1.71</strong></td><td>1.35e-04</td><td><strong>0.76</strong></td></tr><tr><td>7.81e-03</td><td>53248</td><td>1542</td><td>3.06e-06</td><td><strong>2.00</strong></td><td>1.25e-05</td><td><strong>0.77</strong></td><td>1.15e-04</td><td><strong>0.24</strong></td></tr></tbody></table>

Neumann boundary conditions are set on the square. For the sinusoidal case, Dirichlet boundary conditions corresponding to the analytical solution are imposed on all the boundaries. An example of computed solutions is given in Figure 15. Convergence results are given in Table 2, for various numbers of refinement levels. In all cases the expected second order over potential is obtained, thus confirming the good performance of the CDO schemes. Actually, even in the presence of cut-cells sometimes of very small size, the theoretical order is reached. Regarding convergence rate over flux field, only a first order is reached, corroborating results obtained in [10] for cell-based CDO scheme on this type of AMR meshes.

<figure><img src="image_18.png" /><figcaption>Figure 15: Immersed Boundary test case: Computed solution for parabolic case (Left) and sinusoidal case (Right).</figcaption></figure>

For all these cases, good computational performances (stability, robustness, size of the linear system, number of iterations of the linear solver, etc.) are maintained, despite the approximation on the boundary
<table><caption>Table 2: Immersed Boundary test case: Convergence results over potential and flux, for two solutions and for various numbers of refinement levels.</caption><thead><tr><th rowspan="2"># level</th><th rowspan="2">Δx<sub>0</sub></th><th rowspan="2">#C</th><th rowspan="2">#F<sup>NC</sup></th><th colspan="4">Parabolic</th><th colspan="4">Sinusoidal</th></tr><tr><th>Er<sub>L<sup>2</sup></sub>(p̃)</th><th>Rate</th><th>Er<sub>L<sup>2</sup></sub>(φ)</th><th>Rate</th><th>Er<sub>L<sup>2</sup></sub>(p̃)</th><th>Rate</th><th>Er<sub>L<sup>2</sup></sub>(φ)</th><th>Rate</th></tr></thead><tbody><tr><td rowspan="4">0</td><td>6.25e-02</td><td>252</td><td>40</td><td>9.61e-04</td><td>-</td><td>6.12e-03</td><td>-</td><td>5.55e-03</td><td>-</td><td>5.67e-02</td><td>-</td></tr><tr><td>3.13e-02</td><td>1000</td><td>88</td><td>2.75e-04</td><td><b>1.81</b></td><td>1.49e-02</td><td><b>-1.27</b></td><td>1.40e-03</td><td><b>1.99</b></td><td>3.38e-02</td><td><b>0.74</b></td></tr><tr><td>1.56e-02</td><td>3992</td><td>160</td><td>7.47e-05</td><td><b>1.88</b></td><td>5.28e-03</td><td><b>1.50</b></td><td>3.56e-04</td><td><b>1.97</b></td><td>1.01e-02</td><td><b>1.75</b></td></tr><tr><td>7.81e-03</td><td>15924</td><td>304</td><td>1.83e-05</td><td><b>2.03</b></td><td>2.75e-03</td><td><b>0.94</b></td><td>8.84e-05</td><td><b>2.01</b></td><td>5.88e-03</td><td><b>0.78</b></td></tr><tr><td rowspan="4">1</td><td>6.25e-02</td><td>424</td><td>220</td><td>8.44e-04</td><td>-</td><td>1.49e-02</td><td>-</td><td>5.56e-03</td><td>-</td><td>3.38e-02</td><td>-</td></tr><tr><td>3.13e-02</td><td>1688</td><td>420</td><td>2.08e-04</td><td><b>2.02</b></td><td>5.28e-03</td><td><b>1.50</b></td><td>1.34e-03</td><td><b>2.06</b></td><td>1.01e-02</td><td><b>1.74</b></td></tr><tr><td>1.56e-02</td><td>6708</td><td>820</td><td>5.03e-05</td><td><b>2.04</b></td><td>2.75e-03</td><td><b>0.94</b></td><td>3.25e-04</td><td><b>2.04</b></td><td>5.88e-03</td><td><b>0.78</b></td></tr><tr><td>7.81e-03</td><td>26712</td><td>1644</td><td>1.24e-05</td><td><b>2.02</b></td><td>1.52e-03</td><td><b>0.86</b></td><td>8.03e-05</td><td><b>2.02</b></td><td>2.97e-03</td><td><b>0.98</b></td></tr><tr><td rowspan="4">2</td><td>6.25e-02</td><td>536</td><td>424</td><td>8.03e-04</td><td>-</td><td>5.33e-03</td><td>-</td><td>5.53e-03</td><td>-</td><td>1.04e-02</td><td>-</td></tr><tr><td>3.13e-02</td><td>2100</td><td>824</td><td>1.90e-04</td><td><b>2.08</b></td><td>2.76e-03</td><td><b>0.95</b></td><td>1.32e-03</td><td><b>2.06</b></td><td>5.94e-03</td><td><b>0.80</b></td></tr><tr><td>1.56e-02</td><td>8280</td><td>1648</td><td>4.62e-05</td><td><b>2.04</b></td><td>1.52e-03</td><td><b>0.86</b></td><td>3.23e-04</td><td><b>2.04</b></td><td>3.00e-03</td><td><b>0.99</b></td></tr><tr><td>7.81e-03</td><td>32920</td><td>3296</td><td>1.13e-05</td><td><b>2.03</b></td><td>8.39e-04</td><td><b>0.86</b></td><td>7.96e-05</td><td><b>2.02</b></td><td>1.76e-03</td><td><b>0.76</b></td></tr></tbody></table>

conditions (approximation of the geometry and approximation of the integrals on the boundary) and despite the presence of sometimes very small cut-cells.

## 5. Conclusions

In this paper, a rewriting of a cell-based CDO scheme for an elliptic problem has been proposed. It is based on the fact that the local Hodge operator is diagonal for Cartesian cells and for SUSHI Hodge operators. This allows to eliminate the flux degrees of freedom at Cartesian faces. The size of the linear system to be solved is therefore greatly reduced (at least if the majority of the cells are Cartesian), which reduce the theoretical computational cost. The proposed rewriting is therefore particularly suited to Cartesian AMR type meshes, for which only the degrees of freedom at faces close to the refinement interfaces and close to the immersed boundaries are preserved. The imposition of more general Robin boundary conditions has also been proposed for these cell-based CDO schemes, using the closure of the dual mesh.

This new approach was then validated on test cases of increasing complexity, both in 2D and 3D, for different boundary conditions, as well as in the presence of immersed boundaries. The expected second order is obtained for all cases, as well as an error up to machine precision for constant and linear solutions. It should be noted however that for solutions presenting strong gradients (or discontinuous solution) second order is no longer preserved, which is an expected behavior. In addition, the maximum principle at discrete level was also satisfied numerically for all cases. Finally, we observe the reduction in the size of the linear systems to be solved thanks to the use of Cartesian AMR meshes, thus validating the choices made in the proposed approach.

For future works, immersed boundaries in 3D and for more complex geometries will be considered. The influence of the choice of preconditioners and solvers on the performance of this new scheme is also left for future studies. Finally, the addition of CDO advection operators [6] will allow to consider an extension of this scheme to the Navier-Stokes equations.

## References

[1] Douglas N Arnold, Richard S Falk, and Ragnar Winther. Finite element exterior calculus, homological techniques, and applications. *Acta numerica*, 15:1–155, 2006.

[2] L Beirão da Veiga, Franco Brezzi, Andrea Cangiani, Gianmarco Manzini, L Donatella Marini, and Alessandro Russo. Basic principles of virtual element methods. *Mathematical Models and Methods in Applied Sciences*, 23(01):199–214, 2013.

[3] L Beirão da Veiga, Franco Brezzi, Luisa Donatella Marini, and Alessandro Russo. The hitchhiker's guide to the virtual element method. *Mathematical models and methods in applied sciences*, 24(08):1541–1573, 2014.
[4] René Beltman. *Mimetic discretizations of the incompressible Navier-Stokes equations for polyhedral meshes*. PhD thesis, Technische Universität Eindhoven, 2020.

[5] René Beltman, Martijn Anthonissen, and Barry Koren. Conservative mimetic cut-cell method for incompressible navier-stokes equations. In *European Conference on Numerical Mathematics and Advanced Applications*, pages 1035–1043. Springer, 2019.

[6] René Beltman, MJH Anthonissen, and Barry Koren. Conservative polytopal mimetic discretization of the incompressible navier-stokes equations. *Journal of Computational and Applied Mathematics*, 340:443–473, 2018.

[7] Michael W Berry, Kyle A Gallivan, Efstratios Gallopoulos, Ananth Grama, Bernard Philippe, Yousef Saad, and Faisal Saied. *High-performance scientific computing: algorithms and applications*. Springer Science & Business Media, 2012.

[8] Pavel B Bochev and James M Hyman. Principles of mimetic discretizations of differential operators. In *Compatible spatial discretizations*, pages 89–119. Springer, 2006.

[9] Jérôme Bonelle. *Opérateurs discrets compatibles pour la discrétisation sur maillages polyédriques des équations elliptiques et de Stokes*. PhD thesis, Paris Est, 2014.

[10] Jérôme Bonelle and Alexandre Ern. Analysis of compatible discrete operator schemes for elliptic problems on polyhedral meshes. *ESAIM: Mathematical Modelling and Numerical Analysis*, 48(2):553–581, 2014.

[11] Alain Bossavit. On the geometry of electromagnetism. *J. Japan Soc. Appl. Electromagn. & Mech*, 6:17–28, 1998.

[12] Alain Bossavit. Computational electromagnetism and geometry. *J. Japan Soc. Appl. Electromagn. & Mech*, 7-8:150–9, 1999.

[13] Jérôme Breil, Pascal Jacq, and Pierre-Henri Maire. A finite volume scheme for solving anisotropic diffusion on ale-amr grids. *International Journal for Numerical Methods in Fluids*, 84(12):737–756, 2017.

[14] Franco Brezzi, Konstantin Lipnikov, and Mikhail Shashkov. Convergence of the mimetic finite difference method for diffusion problems on polyhedral meshes. *SIAM Journal on Numerical Analysis*, 43(5):1872–1896, 2005.

[15] Andrea Cangiani, Zhaonan Dong, and Emmanuil Georgoulis. hp-version discontinuous galerkin methods on essentially arbitrarily-shaped elements. *Mathematics of Computation*, 91(333):1–35, 2022.

[16] Pierre Cantin and Alexandre Ern. Vertex-based compatible discrete operator schemes on polyhedral meshes for advection-diffusion equations. *Computational Methods in Applied Mathematics*, 16(2):187–212, 2016.

[17] Eric Chénier, Robert Eymard, Thierry Gallouët, and Raphaël Herbin. An extension of the mac scheme to locally refined meshes: convergence analysis for the full tensor time-dependent navier–stokes equations. *Calcolo*, 52:69–107, 2015.

[18] Snorre H Christiansen. A construction of spaces of compatible differential forms on cellular complexes. *Mathematical Models and Methods in Applied Sciences*, 18(05):739–757, 2008.

[19] Snorre H Christiansen, Hans Z Munthe-Kaas, and Brynjulf Owren. Topics in structure-preserving discretization. *Acta Numerica*, 20:1–119, 2011.

[20] Bernardo Cockburn, George E Karniadakis, and Chi-Wang Shu. *Discontinuous Galerkin methods: theory, computation and applications*, volume 11. Springer Science & Business Media, 2012.

[21] Lorenzo Codecasa, Ruben Specogna, and Francesco Trevisan. Base functions and discrete constitutive relations for staggered polyhedral grids. *Computer Methods in Applied Mechanics and Engineering*, 198(9-12):1117–1123, 2009.

[22] Lourenço Beirão da Veiga, Konstantin Lipnikov, and Gianmarco Manzini. *The mimetic finite difference method for elliptic problems*, volume 11. Springer, 2014.

[23] Daniele A Di Pietro and Alexandre Ern. A hybrid high-order locking-free method for linear elasticity on general meshes. *Computer Methods in Applied Mechanics and Engineering*, 283:1–21, 2015.

[24] Daniele A Di Pietro, Alexandre Ern, and Simon Lemaire. An arbitrary-order and compact-stencil discretization of diffusion on general meshes based on local reconstruction operators. *Computational Methods in Applied Mathematics*, 14(4):461–472, 2014.

[25] Jozef Dodziuk. Finite-difference approach to the hodge theory of harmonic forms. *American Journal of Mathematics*, 98(1):79–104, 1976.

[26] Jérôme Droniou, Robert Eymard, Thierry Gallouët, and Raphaël Herbin. A unified approach to mimetic finite difference, hybrid finite volume and mixed finite volume methods. *Mathematical Models and Methods in Applied Sciences*, 20(02):265–295, 2010.

[27] Robert Eymard, Thierry Gallouët, and Raphaël Herbin. Discretization of heterogeneous and anisotropic diffusion problems on general nonconforming meshes sushi: a scheme using stabilization and hybrid interfaces. *IMA Journal of Numerical Analysis*, 30(4):1009–1043, 2010.

[28] Robert D Falgout, Jim E Jones, and Ulrike Meier Yang. The design and implementation of hypre, a library of parallel high performance preconditioners. In *Numerical solution of partial differential equations on parallel computers*, pages 267–294. Springer, 2006.

[29] Michel Fortin and Roland Glowinski. *Augmented Lagrangian methods: applications to the numerical solution of boundary-value problems*. Elsevier, 2000.

[30] Jean-Luc Guermond, Peter Minev, and Jie Shen. An overview of projection methods for incompressible flows. *Computer methods in applied mechanics and engineering*, 195(44-47):6011–6045, 2006.

[31] CA Hall, JC Cavendish, and WH Frey. The dual variable method for solving fluid flow difference equations on delaunay triangulations. *Computers & fluids*, 20(2):145–164, 1991.

[32] Francis H Harlow and J Eddie Welch. Numerical calculation of time-dependent viscous incompressible flow of fluid with free surface. *The physics of fluids*, 8(12):2182–2189, 1965.

[33] Ralf Hiptmair. Discrete hodge operators: An algebraic perspective. *Progress In Electromagnetics Research*, 32(247269):27, 2001.

[34] James M Hyman and James C Scovel. Deriving mimetic difference approximations to differential operators using algebraic
topology. *Los Alamos National Laboratory*, 1988.

[35] Konstantin Lipnikov, Gianmarco Manzini, and Mikhail Shashkov. Mimetic finite difference method. *Journal of Computational Physics*, 257:1163–1227, 2014.

[36] Richard Liska and Mikhail Shashkov. Enforcing the discrete maximum principle for linear finite element solutions of second-order elliptic problems. *Commun. Comput. Phys.*, 3(4):852–877, 2008.

[37] Gianmarco Manzini, Alessandro Russo, and N Sukumar. New perspectives on polygonal and polyhedral finite element methods. *Mathematical Models and Methods in Applied Sciences*, 24(08):1665–1699, 2014.

[38] Claudio Mattiussi. The finite volume, finite element, and finite difference methods as numerical methods for physical field problems. In *Advances in Imaging and electron physics*, volume 113, pages 1–146. Elsevier, 2000.

[39] Riccardo Milani, Jérôme Bonelle, and Alexandre Ern. Artificial compressibility methods for the incompressible navier–stokes equations using lowest-order face-based schemes on polytopal meshes. *Computational Methods in Applied Mathematics*, 22(1):133–154, 2022.

[40] Stéphane Popinet. Gerris: a tree-based adaptive solver for the incompressible euler equations in complex geometries. *Journal of computational physics*, 190(2):572–600, 2003.

[41] Alice Raeli, Michel Bergmann, and Angelo Iollo. A finite-difference method for the variable coefficient poisson equation on hierarchical cartesian meshes. *Journal of computational Physics*, 355:59–77, 2018.

[42] Natarajan Sukumar and A21015991073 Tabarraei. Conforming polygonal finite elements. *International Journal for Numerical Methods in Engineering*, 61(12):2045–2066, 2004.

[43] Enzo Tonti et al. *On the formal structure of physical theories*. Istituto di matematica del Politecnico di Milano, 1975.

[44] Sabine Zaglmayr. High order finite element methods for electromagnetic field computation. 2006.

[45] Weiqun Zhang, Andrew Myers, Kevin Gott, Ann Almgren, and John Bell. Amrex: Block-structured adaptive mesh refinement for multiphysics applications. *The International Journal of High Performance Computing Applications*, 35(6):508–526, 2021.