# Basic convergence theory for the network element method

Julien Coatléven

► To cite this version:

Julien Coatléven. Basic convergence theory for the network element method. ESAIM: Mathematical Modelling and Numerical Analysis, 2021, 10.1051/m2an/2021062. hal-03832835

HAL Id: hal-03832835

https://hal.science/hal-03832835

Submitted on 28 Oct 2022

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
# BASIC CONVERGENCE THEORY FOR THE NETWORK ELEMENT METHOD*

JULIEN COATLÉVEN †‡

**Abstract.** A recent paper introduced the network element method (NEM) where the usual mesh was replaced by a discretization network. Using the associated network geometric coefficients and following the virtual element framework, a consistent and stable numerical scheme was proposed. The aim of the present paper is to derive a convergence theory for the NEM under mild assumptions on the exact problem. We also derive basic error estimates, which are sub-optimal in the sense that we have to assume more regularity than usual.

## INTRODUCTION

In recent years, the mimetic technology [10] has proved to be a very efficient tool to derive numerical schemes to handle probably all classical partial differential equations, with general coefficients, even on very distorted or exotic meshes. Many methods, and in particular the virtual element method (VEM, [8]) and the Hybrid-High-Order schemes (HHO, [24]) were developed following its principles, allowing to handle complex problems such as linear [23, 9] and non linear [3, 11] elasticity, parabolic problems [28], multiphase flow problems [5], Stokes problem [12], etc. Based on the success of those polygonal methods, in a recent paper [7] was explored the idea that we probably need less than a mesh to derive an efficient variational numerical method. This naturally led to the notion of discretization networks in [7], which is a common object in meshless methods (see [19, 27, 26]). The network element method was then derived by reproducing the VEM principles directly on the discretization network rather than on a mesh. Numerical examples illustrated the performance of the method, and the excepted convergence rates were observed in practice.

The present paper is an attempt to propose a basis of a convergence theory for the network element method, using again the elementary Poisson problem as a model problem. Notice that the consistency of the method is mainly inherited from the properties of the approximate geometry, as is usual for meshless methods based on discretization networks (see [20, 14, 29, 19, 27, 26]), while its stability comes from its VEM-like (and also discontinuous Galerkin like) formulation. Both were already studied in [7]. Thus, the major difficulty of the convergence analysis consists in going from the purely discrete world of degrees of freedom where the network element method is formulated, to the continuous world of Sobolev spaces. The key ingredient will be a family of functions forming a partition of unity and whose integral will replace the discrete weights of the method. In this way, they play the role usually assumed by a quadrature rule and will form what will naturally be called a quadrature family. This is a major difference with partition of unity finite element methods [22], where a partition of unity is used to decompose the functions and not the integral operator (see [6, 2]). Once this quadrature family is defined, thanks to the variational nature of the network element method the spirit of the convergence theory will be highly reminiscent of finite volume theory, discontinuous Galerkin or the unified framework of gradient discretization (see [16, 17, 13, 15]). In fact, once quadrature families are properly defined and their existence established, one could consider our convergence theory for minimal regularity solutions as establishing that the Gradient discretization framework properties are satisfied by the reconstruction operators

<footnote>*To cite this article, please refer to J. Coatléven, Basic convergence theory for the network element method, ESAIM: M2AN, Vol. 55 (2021), pp. 2503-2533</footnote>

<footnote>†IFP Énergies nouvelles, 1 et 4 avenue de Bois-Préau, 92852 Rueil-Malmaison, France</footnote>

<footnote>‡julien.coatleven@ifpen.fr</footnote>
associated to the quadrature family, with however many additional technicalities coming from the lack of a
mesh, this time very reminiscent of meshless techniques and in particular partition of unity based methods (see
[21]).

As explained in the original paper presenting the network element method [7], when compared to mesh based
methods the performance difference between them and the network element method comes from the cost com-
parison of mesh generation on one side and network and geometry generation on the other side. Comparing
network and geometry generation, it is clear that most of the cost lies in the geometry generation step, however
in best cases this can be done by simply solving *d* or *d* + 1 linear systems (see [7]). Establishing a convergence
theory allowing to understand which quality parameter is crucial to maintain convergence rates is consequently
very helpful for the long term goal of designing fast and robust geometry generation algorithms, which would
truly make the network element method more than a mathematical curiosity.

The paper will be organized as follows: in the first part of the paper (sections 1 and 2), we recall the definitions of a discretization network and the associated geometric weights, as well as the network element method itself. Then, in the second part of the paper (sections 3 to 5) we establish convergence results. Section 3 is devoted to the reconstruction of functions from network element degrees of freedom. In particular, a crucial existence result on quadrature families is established there, which constitutes the backbone of our convergence theory. Section 4 is devoted to convergence to minimal regularity solutions, to emphasize the robustness of the approach. Finally, section 5 deals with error estimates. Notice that once the core theorem of section 3 is established, the last two sections follow the general spirit of finite volume (or discontinuous Galerkin) theory, with additional technicalities specific to the network element method. For numerical experiments illustrating the behavior of the method in practice, we refer the reader to [7].

# 1 DISCRETIZATION NETWORKS AND NETWORK GEOMETRIES

## 1.1 Discretization networks

Let $\Omega$ be an open bounded connected subset of $\mathbb{R}^d$, $d \in \mathbb{N}\{0\}$, assumed to be at least Lipschitz. For any $\mathbf{x} \in \mathbb{R}^d$
and any $r > 0$, we denote $B(\mathbf{x}, r)$ the ball of radius $r$ centered at $\mathbf{x}$ for the usual euclidean norm $|\mathbf{x}|^2 = \sum_{i=1}^{d} x_i^2$.
Following [7, 19, 27, 26], a discretization network $\mathcal{N}$ of $\Omega$ is defined from two sets of points $\mathcal{P}_{\mathcal{T}}$ and $\mathcal{P}_{\mathcal{F}}$, by
setting $\mathcal{N} = \{\mathcal{T}, \mathcal{F}\}$, where :

* The set of cells $\mathcal{T}$ is a set of pairs $K = \{\mathbf{x}_K, r_K\}$, with $\mathbf{x}_K \in \mathcal{P}_{\mathcal{T}}$ strictly inside $\Omega$ and $r_K$ a strictly positive real number, for any $K \in \mathcal{T}$. We denote $h_K = 2r_K$.

* The set of interfaces, denoted $\mathcal{F}$, is a set of pairs $\sigma = \{\mathbf{x}_{\sigma}, \mathcal{T}_{\sigma}\}$, with $\mathbf{x}_{\sigma} \in \mathcal{P}_{\mathcal{F}}$ and $\mathcal{T}_{\sigma}$ a subset of $\mathcal{T}$. It is subdivided into two subsets, the set of boundary interfaces $\mathcal{F}_{\text{ext}}$ and the set of interior interfaces $\mathcal{F}_{\text{int}}$. The set of boundary interfaces $\mathcal{F}_{\text{ext}}$ is such that for all $\sigma \in \mathcal{F}_{\text{ext}}$, $\mathbf{x}_{\sigma}$ is a point in $\cap_{K \in \mathcal{T}_{\sigma}} B(\mathbf{x}_K, r_K) \cap \partial\Omega$. The set of interior interfaces $\mathcal{F}_{\text{int}}$ is such that for all $\sigma \in \mathcal{F}_{\text{int}}$, $\mathbf{x}_{\sigma}$ is a point in $\cap_{K \in \mathcal{T}_{\sigma}} B(\mathbf{x}_K, r_K) \cap \tilde{\Omega}$.

* For all $(K_1, K_2) \in \mathcal{T}^2$ such that $K_1 \neq K_2$, $\mathbf{x}_{K_1} \neq \mathbf{x}_{K_2}$. For all $(\sigma_1, \sigma_2) \in \mathcal{F}^2$ such that $\sigma_1 \neq \sigma_2$, $\mathbf{x}_{\sigma_1} \neq \mathbf{x}_{\sigma_2}$.

* $\Omega \subset \bigcup_{K \in \mathcal{T}} B(\mathbf{x}_K, r_K)$. For any $K \in \mathcal{T}$ such that $\partial\Omega \cap \overline{B(\mathbf{x}_K, r_K)} \neq \emptyset$, then $\mathcal{F}_K \cap \mathcal{F}_{\text{ext}} \neq \emptyset$. For any $(K, L) \in \mathcal{T}^2$ such that $B(\mathbf{x}_K, r_K) \cap B(\mathbf{x}_L, r_L) \neq \emptyset$, then there exists $\sigma \in \mathcal{F}$ such that $(K, L) \subset \mathcal{T}_{\sigma}$.

For any $K \in \mathcal{T}$, we also denote $\mathcal{F}_K = \{\sigma \in \mathcal{F} \mid K \in \mathcal{T}_{\sigma}\}$ (the interfaces of $K$), which implies that for any $\sigma \in \mathcal{F}$, $\mathcal{T}_{\sigma}$ denotes the cells connected to the interface $\sigma$ and satisfies $\mathcal{T}_{\sigma} = \{K \in \mathcal{T} \mid \sigma \in \mathcal{F}_K\}$. We denote $h = \max_{K \in \mathcal{T}} h_K$ and $\mathbb{P}_k(\mathbb{R}^d)$ the set of polynomials of order $k$. A network is said to be admissible if for any cell $K \in \mathcal{T}$, the set $(\mathbf{x}_{\sigma})_{\sigma \in \mathcal{F}_K}$ is unisolvent for first order polynomials (see [7] for details and figure 1 for an example of a network for a curved domain). We also recall the well known result that if $\Omega$ is Lipschitz then it satisfies
<figure><img src="image_2.png" /><figcaption>Figure 1: Example of network associated to a sectorial domain (orange triangles are interfaces, blue circles are cells, lines represent the connectivity)</figcaption></figure>

the cone condition for some angle $\tau$ and radius $r$ (see [1, 18]), i.e. for any $\boldsymbol{x} \in \Omega$, there exists $\boldsymbol{\xi} \in \mathbb{R}^d$ with $|\boldsymbol{\xi}| = 1$ such that $C(x, \boldsymbol{\xi}, \tau, r) \subset \Omega$ where $C(x, \boldsymbol{\xi}, \tau, r)$ denotes the cone:

$$C(x, \boldsymbol{\xi}, \tau, r) = B(\boldsymbol{x}, r) \cap \{\boldsymbol{y} \in \mathbb{R}^d \mid (\boldsymbol{y} - \boldsymbol{x})^T \boldsymbol{\xi} > |\boldsymbol{y} - \boldsymbol{x}| \cos \tau\}. \quad (1)$$

Still using the fact that $\Omega$ is assumed Lipschitz, using Stein's extension theorem [25] we also know that there exists an operator $E$ such that for any $k \ge 0$, there exists $C_{E,k} > 0$ such that for any $v \in H^k(\Omega)$, $Ev \in H^k(\mathbb{R}^d)$, $Ev = v$ in $\Omega$ and

$$|Ev|_{H^k(\mathbb{R}^d)} \le C_{E,k}|v|_{H^k(\Omega)},$$

and if $v \in H_0^1(\Omega)$, then $Ev = 0$ in $\mathbb{R}^d \setminus \Omega$. Finally, for any subset $\mathcal{O}$ of $\mathbb{R}^d$, we denote $\chi_{\mathcal{O}}$ the characteristic function of $\mathcal{O}$, i.e. $\chi(\boldsymbol{x}) = 1$ if $\boldsymbol{x} \in \mathcal{O}$ and $\chi(\boldsymbol{x}) = 0$ otherwise.

## 1.2 Network geometry

Following [7], A network geometry is defined as a set of coefficients:

$$\mathcal{G} = \left( (m_K)_{K \in \mathcal{T}}, (\eta_{K,\sigma})_{K \in \mathcal{T}, \sigma \in \mathcal{F}_K}, (\varepsilon_K^{0,i})_{K \in \mathcal{T}, 1 \le i \le d}, (\varepsilon_K^{1,ij})_{K \in \mathcal{T}, 1 \le i, j \le d}, (\varepsilon_\sigma^i)_{\sigma \in \mathcal{F}_{int}, 1 \le i \le d} \right).$$

The discrete measures $(m_K)_{K \in \mathcal{T}}$ are said to be admissible if and only if they satisfy

$$m_K > 0 \quad \text{for all } K \in \mathcal{T}, \qquad (2)$$

and

$$\sum_{K \in \mathcal{T}} m_K = |\Omega|, \qquad (3)$$

while the approximate consistency properties are given by

$$\sum_{\sigma \in \mathcal{F}_K} \eta_{K, \sigma}^i = m_K \varepsilon_K^{0, i} \quad \forall K \in \mathcal{T}, \forall 1 \le i \le d, \qquad (4)$$

and

$$\sum_{\sigma \in \mathcal{F}_K} \eta_{K, \sigma}^i (x_{\sigma}^j - x_K^j) = m_K (\delta_{ij} + \varepsilon_K^{1,ij}) \quad \forall K \in \mathcal{T}, \forall 1 \le i, j \le d, \qquad (5)$$
and the approximate compatibility (or conservation) properties by

$$ \sum_{K \in \mathcal{T}_{\sigma}} \eta_{K\sigma}^{i} = \varepsilon_{\sigma}^{i} \quad \forall \sigma \in \mathcal{F}_{\text{int}}, \forall 1 \le i \le d. \qquad (6) $$

A network geometry is said to be consistent if and only if it satisfies (4)-(5), and said to be conservative if and only if it satisfies (6). To measure the geometric approximation error, we introduce the constants $\theta_A > 0$ and $p \ge 1$, both independent on $h$ and such that:

$$ |\varepsilon_K^{0,i}| \le \theta_A h_K^p \quad \forall K \in \mathcal{T}, \forall 1 \le i \le d, \qquad (7) $$

and

$$ |\varepsilon_K^{1,ij}| \le \theta_A h_K^p \quad \forall K \in \mathcal{T}, \forall 1 \le i, j \le d, \qquad (8) $$

and

$$ |\varepsilon_{\sigma}^{i}| \le \theta_A \min_{K \in \mathcal{T}_{\sigma}} m_K h_K^p \quad \forall \sigma \in \mathcal{F}_{\text{int}}, \forall 1 \le i \le d. \qquad (9) $$

We denote $B_K = B(x_K, r_K)$, and for any $x \in \mathbb{R}^d$, we denote

$$ \mathcal{T}_{x}^{N} = \{K \in \mathcal{T} \mid x \in B_{K}\} \quad \text{and} \quad \eta_{N} = \sup_{x \in \mathbb{R}^{d}} \mathrm{card}(\mathcal{T}_{x}^{N}). $$

We say that a network geometry is admissible if and only if it is consistent and conservative and the family of measures is admissible. As soon as $N$ is an admissible network, existence of an admissible network geometry was established in [7].

**Remark 1.1.** Here, we have chosen to slightly simplify condition (3) regarding the original and more general notion of [7]:

$$ \sum_{K \in \mathcal{T}} m_K = (1 + \varepsilon_{\Omega}) |\Omega|, \qquad (10) $$

which allowed an additional error $\varepsilon_{\Omega}$ on the sum of the discrete measures. However, once one has computed measures $(\tilde{m}_K)_{K \in \mathcal{T}}$ satisfying the above approximate relation (10), it is always feasible to define:

$$ m_K = \frac{\tilde{m}_K |\Omega|}{\sum_{L \in \mathcal{T}} \tilde{m}_L} = \frac{\tilde{m}_K}{(1 + \varepsilon_{\Omega})}. $$

Indeed, we then have for the $\eta_{K,\sigma}$'s corresponding to those $(\tilde{m}_K)_{K \in \mathcal{T}}$:

$$ \sum_{\sigma \in \mathcal{F}_K} \eta_{K,\sigma}^{i}(x_{\sigma}^{j} - x_{K}^{j}) = \tilde{m}_K(\delta_{ij} + \varepsilon_K^{1,ij}) = m_K(1 + \varepsilon_{\Omega})(\delta_{ij} + \varepsilon_K^{1,ij}) = m_K(\delta_{ij} + \varepsilon_{\Omega}\delta_{ij} + (1 + \varepsilon_{\Omega})\varepsilon_K^{1,ij}). $$

The last term $\hat{\varepsilon}_K^{1,ij} = \varepsilon_{\Omega}\delta_{ij} + (1 + \varepsilon_{\Omega})\varepsilon_K^{1,ij}$ is bounded by $2\theta_A h_K^p + \theta_A^2 h_K^{2p} = \hat{\theta}_A h_K^p$, with $\hat{\theta}_A = 2\theta_A + \theta_A^2 h_K^p$. As the same holds for relations (4), we see that up to a modification of the value of $\theta_A$ using $\hat{\theta}_A$, we can always assume that $\sum_{K \in \mathcal{T}} m_K = |\Omega|$. In other words, if we are able to derive a convergence theory assuming the exact relation (3), then this convergence theory will also cover the more general notion of [7], which is the reason why we only consider the simplest version here.

# 2 THE NETWORK ELEMENT METHOD

## 2.1 Model problem

As in [7], to ease the understanding we consider the simplest possible model problem, i.e. the Poisson equation $-\Delta u = f$ on $\Omega$ with $f \in L^2(\Omega)$. We complement it with homogeneous Dirichlet boundary conditions $u = 0$ on
$\partial\Omega = \bar{\Omega}\setminus\Omega$, the boundary of the domain $\Omega$ assumed to be at least Lipschitz continuous. The associated weak solution is the unique $u \in H_0^1(\Omega)$ such that:

$$ \int_{\Omega} \nabla u \cdot \nabla v = \int_{\Omega} fv \quad \forall v \text{ in } H_0^1(\Omega) \Leftrightarrow a(u, v) = l(v) \quad \forall v \text{ in } H_0^1(\Omega). \quad (11) $$

## 2.2 Degrees of freedom and discrete variational formulation

The space of degrees of freedom is given by:

$$ X_N = \{(u_\sigma)_{\sigma \in \mathcal{F}} \mid u_\sigma \in \mathbb{R} \ \forall \ \sigma \in \mathcal{F}\} \quad \text{and} \quad X_{N,0} = \{U \in X_N \mid u_\sigma = 0 \text{ for all } \sigma \in \mathcal{F}_{\text{ext}}\}. $$

The local set of degrees of freedom associated to a cell is denoted

$$ X_{N,K} = \{(u_\sigma)_{\sigma \in \mathcal{F}_K} \mid u_\sigma \in \mathbb{R} \ \forall \ \sigma \in \mathcal{F}_K\}. $$

We denote $U = (u_\sigma)_{\sigma \in \mathcal{F}}$, and for any $U \in X_N$, $U_K = (u_\sigma)_{\sigma \in \mathcal{F}_K}$. To any cell $K \in \mathcal{T}$ is associated a point $\bar{x}_K$ such that:

$$ \bar{x}_K = \sum_{\sigma \in \mathcal{F}_K} \gamma_{K,\sigma} \bar{x}_\sigma \quad \text{where} \quad \sum_{\sigma \in \mathcal{F}_K} \gamma_{K,\sigma} = 1 $$

where the $(\gamma_{K,\sigma})_{\sigma \in \mathcal{F}_K}$ forms a barycentric interpolation for $\bar{x}_K$ from the interface points $(\bar{x}_\sigma)_{\sigma \in \mathcal{F}_K}$. Then we denote:

$$ \mathcal{M}_K(U_K) = \sum_{\sigma \in \mathcal{F}_K} \gamma_{K,\sigma} u_\sigma. $$

To any cell $K \in \mathcal{T}$, is associated the local reconstruction operator $\Pi_K$ defined by:

$$ \left| \begin{array}{lclcl} \Pi_K : & X_{N,K} & \to & \mathbb{P}_1(\mathbb{R}^d) & \\ & U_K & \to & \Pi_K(U_K) = \mathcal{M}_K(U_K) + \nabla_K(U_K) \cdot (\boldsymbol{x} - \bar{\boldsymbol{x}}_K), & \end{array} \right. \quad (12) $$

where

$$ \left| \begin{array}{lclcl} \nabla_K : & X_{N,K} & \to & \mathbb{P}_0(\mathbb{R}^d)^d & \\ & U_K & \to & \nabla_K(U_K) = \frac{1}{m_K} \sum_{\sigma \in \mathcal{F}_K} u_\sigma \eta_{K,\sigma}, & \end{array} \right. \quad (13) $$

with of course $\nabla\Pi_K(U_K) = \nabla_K(U_K)$, thus we will use either one notation or the other in the following. Finally, with a slight abuse of notation we extend the definition of $\Pi_K$, $\mathcal{M}_K$ and $\nabla_K$ to all $X_N$ by setting $\Pi_K(U) = \Pi_K(U_K)$, $\mathcal{M}_K(U) = \mathcal{M}_K(U_K)$ and $\nabla_K(U) = \nabla_K(U_K)$. For any $\varphi \in C^0(\mathbb{R}^d)$ (and more generally for any function for which it makes sense), we denote $\mathcal{D}_K(\varphi) = (\varphi(\bar{x}_\sigma))_{\sigma \in \mathcal{F}_K}$ the local set of degrees of freedom associated with $\varphi$, while $\mathcal{D}(\varphi) = (\varphi(\bar{x}_\sigma))_{\sigma \in \mathcal{F}}$ denotes the complete set of degrees of freedom associated with $\varphi$.

Let us recall some key ideas underlying the virtual element method and by this way also recall the ideas underlying the derivation of the network element method in [7]. Assume that we are given a true mesh of $\Omega$ whose set of cells is denoted $\mathcal{T}$ to make the analogy with networks more obvious. Using the virtual element projector $\Pi_K^{VEM}$ onto first order polynomial functions (and denoting $\pi_k$ the $L^2$ projection on polynomials of order $k$) the simplest first order virtual element method consist in solving

$$ \sum_{K \in \mathcal{T}} \int_K \nabla \Pi_K^{VEM}(u) \cdot \nabla \Pi_K^{VEM}(v) + \sum_{K \in \mathcal{T}} s_K^{VEM}(u - \Pi_K^{VEM}(u), v - \Pi_K^{VEM}(v)) = \sum_{K \in \mathcal{T}} \int_K f \pi_0(v), $$

where the first term handles the polynomial part of the unknown VEM function to ensure consistency and $s_K^{VEM}(u - \Pi_K^{VEM}(u), v - \Pi_K^{VEM}(v))$ is a stabilization bilinear form which only needs to scale with $h_K$ in the
same way than the term it replaces to preserve consistency. The network element mimics the principles of the virtual element method, but using a discretization network rather than a mesh. The discrete bilinear form is constructed by analogy on the discretization network and its associated network geometry using the discrete gradient $\nabla_K$ to handle the polynomial part (replacing the $\nabla_K^{VEM}(u)$ of the VEM) and then complemented with a stabilization term which has the same form as $s_K^{VEM}$ and has the correct scaling to maintain consistency. Consequently, the discrete counterpart $a_h : X_N \times X_N \to \mathbb{R}$ of the bilinear form $a(\cdot, \cdot)$ is defined by setting

$$a_h(U, V) = \sum_{K \in \mathcal{T}} a_h^K(U_K, V_K),$$

where $a_h^K : X_{N,K} \times X_{N,K} \to \mathbb{R}$ is given by

$$a_h^K(U_K, V_K) = m_K \nabla \Pi_K(U_K) \cdot \nabla \Pi_K(V_K) + s^K(U_K - \mathcal{D}_K(\Pi_K(U_K)), V_K - \mathcal{D}_K(\Pi_K(V_K))), \quad (14)$$

with $s^K$ a positive symmetric bilinear form on $X_{N,K} \times X_{N,K}$, such that

$$s^K(U_K, V_K) = m_K h_K^{-2} \sum_{\sigma \in \mathcal{F}_K} \sum_{\sigma' \in \mathcal{F}_K} S_{K, \sigma, \sigma'} u_\sigma v_{\sigma'}, \quad (15)$$

where $S_K = (S_{K,\sigma,\sigma'})_{\sigma,\sigma' \in \mathcal{F}_K}$ can be any symmetric positive definite matrix independent on the geometry $\mathcal{G}$ associated to the network, for which we denote

$$S_* = \inf_{K \in \mathcal{T}} \inf_{\boldsymbol{\xi} \in \mathbb{R}^{\mathrm{card}(\mathcal{F}_K)}, ||\boldsymbol{\xi}||=1} \boldsymbol{\xi}^T S_K \boldsymbol{\xi} \quad \text{and} \quad S^* = \sup_{K \in \mathcal{T}} \sup_{\boldsymbol{\xi} \in \mathbb{R}^{\mathrm{card}(\mathcal{F}_K)}, ||\boldsymbol{\xi}||=1} \boldsymbol{\xi}^T S_K \boldsymbol{\xi}.$$

For the right-hand side, assume that $f_K$ is an approximation of $f$ at $\bar{x}_K$ (for instance, one can use $f(\bar{x}_K)$ if $f$ is regular enough for this quantity to make sense, or $\frac{1}{|B_K|} \int_{B_K} Ef$), then we define a linear form $l_h : X_N \to \mathbb{R}$ by setting:

$$l_h(V) = \sum_{K \in \mathcal{T}} m_K f_K \mathcal{M}_K(V_K).$$

Then, the discretization by the network element method consists in finding a solution $U \in X_{N,0}$ of

$$a_h(U, V) = l_h(V) \quad \text{for all } V \in X_{N,0}. \quad (16)$$

## 2.3 Basic properties of the network element method

The spaces $X_{N,K}$ of degrees of freedom are endowed with the bilinear forms:

$$(U, V)_{0,K} = m_K \mathcal{M}_K(U_K) \mathcal{M}_K(V_K) + \sum_{\sigma \in \mathcal{F}_K} m_K (u_\sigma - \mathcal{M}_K(U_K))(v_\sigma - \mathcal{M}_K(V_K)),$$

and

$$(U, V)_{1,K} = \sum_{\sigma \in \mathcal{F}_K} m_K h_K^{-2} (u_\sigma - \mathcal{M}_K(U_K))(v_\sigma - \mathcal{M}_K(V_K)),$$

and the associated norm $||U||_{0,K}^2 = (U, U)_{0,K}$ and semi-norm $|U|_{1,K}^2 = (U, U)_{1,K}$, while we denote $||U||_{X,K}^2 = ||U||_{0,K}^2 + |U|_{1,K}^2$. We endow the space of degrees of freedom $X_N$ with the bilinear forms

$$(U, V)_0 = \sum_{K \in \mathcal{T}} (U, V)_{0,K} \quad \text{and} \quad (U, V)_1 = \sum_{K \in \mathcal{T}} (U, V)_{1,K},$$

and the associated norm $||U||_0^2 = (U, U)_0$ and semi-norm $|U|_1^2 = (U, U)_1$. Then we define:

$$(U, V)_X = (U, V)_0 + (U, V)_1 \quad \text{and} \quad ||U||_X^2 = (U, U)_X,$$
which are obviously a scalar product and its associated norm on $X_N$, making $X_N$ a Hilbert space. We recall now the measures of quality of the discretization network and its associated geometry:

$$ \theta_{\Pi} = \sup_{K \in \mathcal{T}} \sup_{\sigma \in \mathcal{F}_K} h_K \left| \frac{\eta_{K,\sigma}}{m_K} \right| \quad \text{and} \quad \theta_{\mathcal{M}} = \sup_{K \in \mathcal{T}} \sup_{\sigma \in \mathcal{F}_K} |\gamma_{K,\sigma}| \quad \text{and} \quad \theta_{\mathcal{F}} = \max_{K \in \mathcal{T}} \mathrm{card}(\mathcal{F}_K), $$

and

$$ \theta_{\mathcal{T}} = \sup_{K \in \mathcal{T}} \max \left( \frac{|B_K \cap \Omega|}{m_K}, \frac{m_K}{|B_K \cap \Omega|} \right), $$

and we denote $S_1^d = |B(0,1)|$ the volume of the unit ball in dimension $d$. Using the quality measures and the above norms, it was established in [7] that there exists $C_\nabla > 0$ depending only on the quality parameters and independent on $h$ such for any $U \in X_N$:

$$ m_K |\nabla_K(U)|^2 \le C_\nabla^2 (|U|_{1,K}^2 + m_K^2) |\mathcal{M}_K(U_K)|^2. \quad (17) $$

Moreover, there exists $C_a > 0$ depending only on $S^*$ and the quality parameters and independent on $h$ such that for any $(U, V) \in X_{N,K}^2$:

$$ a_h^K(U, V) \le C_a ||U||_{X,K} ||V||_{X,K}, \quad (18) $$

while for any $(U, V) \in X_N^2$:

$$ a_h(U, V) \le C_a ||U||_X ||V||_X. \quad (19) $$

There also exists $\alpha_a > 0$ depending only on $S_*$ and the quality parameters and independent on $h$ such that

$$ a_h^K(U, U) \ge \alpha_a |U|_{1,K}^2 \quad \text{for any } U \in X_{N,K} \quad \text{and} \quad a_h(U, U) \ge \alpha_a |U|_1^2 \quad \text{for any } U \in X_N. \quad (20) $$

Finally, for any $V \in X_N$

$$ l_h(V) \le C_f |V|_0 \quad \text{where} \quad C_f = \left( \sum_{K \in \mathcal{T}} m_K |f_K|^2 \right)^{\frac{1}{2}}. \quad (21) $$

Assume that $\Omega$ satisfies the cone condition with angle $\tau$ and radius $r$, and denote $\delta > 0$ the smallest real number such that for any $K \in \mathcal{T}$:

$$ \delta^{-1} r_K \le \min(r, r_K) \le \delta r_K. \quad (22) $$

Then, there exists $C_{P,X_N} > 0$ depending on $\tau$, $\delta$, $\eta_N$, $\theta_{\mathcal{T}}$ and $\Omega$ such that the following discrete Poincaré's inequality also holds, making existence, uniqueness and stability of the discrete solution an obvious consequence of Lax-Milgram's lemma:

$$ ||U||_0^2 \le C_{P,X_N} |U|_1^2. \quad (23) $$

### 3 RECONSTRUCTION OPERATORS FOR NETWORK ELEMENT SOLUTIONS

Our main objective in the present paper is to characterize the approximation properties of the network element method. The main difficulty comes from the fact that we have only worked in the degree of freedom (dof) space $X_N$, and in particular, we have not defined any function on $\Omega$ as it would be the case in classical variational methods, mesh-based or not. The first task is consequently to bridge this gap between the discrete and continuous worlds.
## 3.1 Quadrature families

To this end, we consider any family of functions $(\psi_K)_{K \in \mathcal{T}}$ such that for any $K \in \mathcal{T}$, $\psi_K \in L^\infty(\mathbb{R}^d)$ and:

$$ \int_{\Omega} \psi_K = m_K, \quad \sum_{K \in \mathcal{T}} \psi_K = 1 \quad \text{for a.e } \boldsymbol{x} \in \Omega, \quad \text{supp } \psi_K \subset B(\boldsymbol{x}_K, \rho_K) \supset B_K, \qquad (24) $$

and we denote $\mathcal{B}_K = B(\boldsymbol{x}_K, \rho_K)$. Such a family $(\psi_K)_{K \in \mathcal{T}}$ is called a quadrature family, while the set of all quadrature families is denoted $\Psi(\mathcal{N}, \mathcal{G})$. For any $x \in \mathbb{R}^d$, we denote

$$ \mathcal{T}_{\boldsymbol{x}}^{\mathcal{B}} = \{K \in \mathcal{T} \mid \boldsymbol{x} \in \mathcal{B}_K\} \quad \text{and} \quad \eta_{\psi} = \sup_{\boldsymbol{x} \in \mathbb{R}^d} \mathrm{card}(\mathcal{T}_{\boldsymbol{x}}^{\mathcal{B}}), \qquad (25) $$

$$ \kappa_{\psi} = \max \left( \max_{K \in \mathcal{T}} \frac{\rho_K}{r_K}, \max_{K \in \mathcal{T}} \frac{r_K}{\rho_K} \right) \quad \text{and} \quad M_{\psi} = \max_{K \in \mathcal{T}} \|\psi_K\|_{L^{\infty}(\Omega)}. \qquad (26) $$

We call $(\eta_{\psi}, \kappa_{\psi}, M_{\psi})$ the parameters of a quadrature family. From these definitions we immediately deduce that $\kappa_{\psi} > 0$, $M_{\psi} > 0$ and that

$$ \kappa_{\psi}^{-1} r_K \le \rho_K \le \kappa_{\psi} r_K \quad \text{and} \quad \|\psi_K\|_{L^{\infty}(\Omega)} \le M_{\psi}. \qquad (27) $$

Remark that the hypothesis $B_K \subset \mathcal{B}_K$ also immediately implies that $\kappa_{\psi} \ge 1$. If the domain $\Omega$ satisfies the cone condition with angle $\tau$ and radius $r$, and if $\delta > 0$ is defined as in (22), then noticing that $|C(0, \boldsymbol{\xi}, \tau, 1)|$ is in fact independent on $\boldsymbol{\xi}$ and denoting $|C(0, \boldsymbol{\xi}, \tau, 1)| = |C(0, \tau, 1)|$ this common value, we have

$$ |B_K \cap \Omega| \ge |C(\boldsymbol{x}_K, \boldsymbol{\xi}, \tau, \min(r, r_K))| = |C(0, \boldsymbol{\xi}, \tau, 1)| \min(r, r_K)^d \ge |C(0, \tau, 1)| \min(r, r_K)^d $$

leading to the useful inequality:

$$ \frac{1}{m_K} \int_{\Omega} |\psi_K| \le M_{\psi} \theta_{\tau} \frac{S_1^d \kappa_{\psi}^d \delta^d}{|C(0, \tau, 1)|}. \qquad (28) $$

We can now construct functions on $\mathbb{R}^d$ and in particular on $\Omega$ using the $\psi_K$'s by setting:

$$ \Pi_{\mathcal{T}}(U) = \sum_{K \in \mathcal{T}} \psi_K \mathcal{M}_K(U) \quad \text{and} \quad \nabla_{\mathcal{T}}(U) = \sum_{K \in \mathcal{T}} \psi_K \nabla_K(U) \quad \text{and} \quad \Pi_{\mathcal{N}}(U) = \sum_{K \in \mathcal{T}} \psi_K \Pi_K(U_K). $$

Those reconstructions will be the key to establish convergence results. Before turning to it, let us define

$$ \|U\|_{\mathcal{T}}^2 = \sum_{K \in \mathcal{T}} m_K \mathcal{M}_K(U)^2. $$

Obviously, for any $U \in X_{\mathcal{N}}$, $\|U\|_{\mathcal{T}} \le \|U\|_{X}$. Moreover:

**Lemma 3.1.** Let $(\mathcal{N}, \mathcal{G})$ be an admissible discretization network and an associated admissible network geometry. Assume that $\Omega$ satisfies the cone condition with angle $\tau$ and radius $r$. Then there exists $C > 0$ depending only on $\Omega$ and the quality parameters of the geometry and of the quadrature family such that for any $U \in X_{\mathcal{N}}$:

$$ \|\Pi_{\mathcal{T}}(U)\|_{L^2(\Omega)} \le C \|U\|_{\mathcal{T}} \quad \text{and} \quad \|\nabla_{\mathcal{T}}(U)\|_{L^2(\Omega)} \le C \|U\|_{X}. $$

*Proof.* By definition, we have:

$$ \|\Pi_{\mathcal{T}}(U)\|_{L^2(\Omega)}^2 = \int_{\Omega} \left( \sum_{K \in \mathcal{T}} \psi_K \mathcal{M}_K(U) \right)^2 \le \int_{\Omega} \left( \sum_{K \in \mathcal{T}} |\psi_K| \mathcal{M}_K(U)^2 \right) \left( \sum_{K \in \mathcal{T}} |\psi_K| \right). $$
Then notice that (24) implies that only $\eta_\psi$ terms are non zero in the second sum and thus using (28):

$$ ||\Pi_{\mathcal{T}}(\mathbf{U})||_{L^2(\Omega)}^2 \le \eta_{\psi} M_{\psi} \sum_{K \in \mathcal{T}} \left( \int_{\Omega} |\psi_K| \right) \mathcal{M}_K(\mathbf{U})^2 \le \eta_{\psi} M_{\psi}^2 \theta_{\mathcal{T}} \frac{S_1^d \kappa_{\psi}^d \delta^d}{|C(0, \tau, 1)|} ||\mathbf{U}||_{\mathcal{T}}^2. $$

Proceeding exactly in the same way, we obtain:

$$ ||\nabla_{\mathcal{T}}(\mathbf{U})||_{L^2(\Omega)}^2 \le \eta_{\psi} M_{\psi} \sum_{K \in \mathcal{T}} \left( \int_{\Omega} |\psi_K| \right) |\nabla_K(\mathbf{U})|^2, $$

and the second result follows using (17). □

## 3.2 Existence of quadrature families

To establish convergence and error estimates, we will not only need the existence of a quadrature family, but also some control over its parameters, independently of the mesh size. However as we have not used any specific partition in practice in the construction of the numerical scheme, this should be considered as a technical requirement to construct a convergence theory. Moreover, such a theory will in fact depend on optimal bounds for those constants. Deriving such optimal bounds over the entire set $\Psi(\mathcal{N}, \mathcal{G})$ using only properties of the point cloud is in fact a very difficult problem: one way to do it would consist in first defining what would be a relevant measure of optimality (clearly a compromise must be found between $(\kappa_\psi, \eta_\psi)$ and $M_\psi$) and then constructing either a minimizing sequence of quadrature families or a quadrature family reaching this optimal compromise. This is the reason why we introduce a more specific type of quadrature families, with the main advantage that it will provide a practical mean to compute an upper bound on those constants, and the obvious drawback that this bound could remain very pessimistic.

Consider any family $(\psi_\sigma)_{\sigma \in \mathcal{F}}$ such that $\text{supp} \psi_\sigma \subset B(\mathbf{x}_\sigma, r_\sigma)$ and

$$ \bar{\Omega} \subset \bigcup_{\sigma \in \mathcal{F}} \overline{B(\mathbf{x}_\sigma, r_\sigma)} \quad \text{and} \quad \sum_{\sigma \in \mathcal{F}} \psi_\sigma(\mathbf{x}) = 1 \quad \forall \mathbf{x} \in \Omega \quad \text{and} \quad 0 \le \psi_\sigma \le 1 \quad \forall \sigma \in \mathcal{F}. \quad (29) $$

Denoting $m_\sigma = \int_{\Omega} \psi_\sigma$, we consider a family of weights $\omega = (\omega_{K\sigma})_{K \in \mathcal{T}, \sigma \in \mathcal{F}_K}$ for which the family of functions $(\psi_K)_{K \in \mathcal{T}}$ defined by:

$$ \psi_K = \sum_{\sigma \in \mathcal{F}_K} \omega_{K,\sigma} \psi_\sigma, \quad (30) $$

is a quadrature family. To this end, simply injecting formula (30) in conditions (24) and (27) immediately leads to:

$$ \sum_{K \in \mathcal{T}} \sum_{\sigma \in \mathcal{F}_K} \omega_{K,\sigma} \psi_\sigma = 1 \quad \text{and} \quad \sum_{\sigma \in \mathcal{F}_K} \omega_{K,\sigma} m_\sigma = m_K \quad \text{where} \quad m_\sigma = \int_{\Omega} \psi_\sigma. $$

Rearranging the sums in the first condition, we get:

$$ \sum_{\sigma \in \mathcal{F}} \left( \sum_{K \in \mathcal{T}_\sigma} \omega_{K,\sigma} \right) \psi_\sigma = 1, $$

and thus as $(\psi_\sigma)_{\sigma \in \mathcal{F}}$ is a partition of unity it is sufficient to find a solution to the following linear system:

$$ A_{\psi} \boldsymbol{\omega} = \mathbf{1} \iff \left| \begin{array}{l} \displaystyle \sum_{K \in \mathcal{T}_{\sigma}} \omega_{K, \sigma} = 1 \quad \forall \sigma \in \mathcal{F} \\ \\ \displaystyle \sum_{\sigma \in \mathcal{F}_K} \omega_{K, \sigma} \frac{m_{\sigma}}{m_K} = 1 \quad \forall K \in \mathcal{T}, \end{array} \right. \quad (31) $$
to get a quadrature family. We say that such a quadrature family is an interface based quadrature family. The existence of families belonging to this subclass of quadrature family is the object of the following theorem, which plays a crucial role in our convergence theory. In fact, once the existence of this interface based quadrature family established, the spirit of the convergence theory will then be reminiscent of finite volume theory, with the major additional difficulty that we have to handle a discretization network rather than a mesh. The following proposition and the associated corollary can thus undoubtedly be considered as the main results of the present paper.

**Proposition 3.2.** If $\mathcal{N}$ is an admissible network and $\mathcal{G}$ an associated admissible geometry, then there exists a solution to (31) for any family $(m_\sigma)_{\sigma \in \mathcal{F}}$ such that

$$ \sum_{K \in \mathcal{T}} m_K = \sum_{\sigma \in \mathcal{F}} m_\sigma = |\Omega|. $$

*Proof.* Let $\mathbf{y} = ((y_\sigma)_{\sigma \in \mathcal{F}}, (y_K)_{K \in \mathcal{T}}) \in \mathbb{R}^{\mathrm{card}(\mathcal{F})+\mathrm{card}(\mathcal{T})}$ be such that $\mathbb{A}_\psi^T \mathbf{y} = 0$. We have:

$$ \omega^T \mathbb{A}_\psi^T \mathbf{y} = (\mathbb{A}_\psi \omega)^T \mathbf{y} = \sum_{\sigma \in \mathcal{F}} \sum_{K \in \mathcal{T}_\sigma} \omega_{K,\sigma} y_\sigma + \sum_{K \in \mathcal{T}} \sum_{\sigma \in \mathcal{F}_K} \omega_{K,\sigma} \frac{m_\sigma}{m_K} y_K = \sum_{K \in \mathcal{T}} \sum_{\sigma \in \mathcal{F}_K} \omega_{K,\sigma} \left( y_\sigma + \frac{m_\sigma}{m_K} y_K \right), $$

and thus $\mathbb{A}_\psi^T \mathbf{y} = 0$ is equivalent to:

$$ y_\sigma + \frac{m_\sigma}{m_K} y_K = 0 \quad \forall K \in \mathcal{T}, \forall \sigma \in \mathcal{F}_K. $$

Consequently, we have:

$$ \frac{y_\sigma}{m_\sigma} = - \frac{y_K}{m_K} \quad \forall K \in \mathcal{T}, \forall \sigma \in \mathcal{F}_K, $$

For any $(K, L) \in \mathcal{T}$, by definition of the discretization network there exists cells $(K_m)_{0 \le m \le d(K,L)}$ and interfaces $(\sigma_m)_{0 \le m \le d(K,L)-1}$ such that $K_0 = K$, $K_{d(K,L)} = L$ and $\{K_m, K_{m+1}\} \subset \mathcal{T}_{\sigma_m}$ for all $0 \le m \le d(K,L)-1$. Consequently:

$$ \frac{y_{\sigma_m}}{m_{\sigma_m}} = - \frac{y_{K_m}}{m_{K_m}} \quad \text{and} \quad \frac{y_{\sigma_m}}{m_{\sigma_m}} = - \frac{y_{K_{m+1}}}{m_{K_{m+1}}}. $$

Then there exists a constant $\alpha$ such that $y_K = \alpha m_K$ and $y_\sigma = -\alpha m_\sigma$ for all $K \in \mathcal{T}$ and all $\sigma \in \mathcal{F}$. We have, denoting $\mathbf{1}$ the vector with all components equal to 1:

$$ \mathbf{1}^T y = \sum_{\sigma \in \mathcal{F}} y_\sigma + \sum_{K \in \mathcal{T}} y_K = \alpha \left( \sum_{K \in \mathcal{T}_i} m_K - \sum_{\sigma \in \mathcal{F}_i} m_\sigma \right) = 0 \quad \text{as} \quad \sum_{K \in \mathcal{T}} m_K = \sum_{\sigma \in \mathcal{F}} m_\sigma = |\Omega|, $$

which from Fredholm alternative establishes the existence of a family $(\omega_{K\sigma})_{K \in \mathcal{T}, \sigma \in \mathcal{F}_K}$ satisfying (31). $\square$

As we will see in the following corollary (otherwise we refer the reader to [21]), a partition of unity satisfying (29) always exists. Thus as an immediate consequence of the above proposition we know that an interface based quadrature family also exists and thus $\Psi(\mathcal{N}, \mathcal{G})$ is non empty.

**Corollary 3.3.** Assume that $\Omega$ satisfies the cone condition with angle $\tau$ and radius $r$. Let $\mathcal{N}$ be an admissible network and $\mathcal{G}$ an associated admissible geometry. Assume that there exists $0 < \alpha < 1$ such that

$$ \bar{\Omega} \subset \bigcup_{\sigma \in \mathcal{F}} \overline{B(\mathbf{x}_\sigma, \alpha r_\sigma)} \quad \text{where} \quad r_\sigma = \max_{K \in \mathcal{T}_\sigma} r_K. \qquad (32) $$
Then there exists an interface based quadrature family for which $\kappa_{\psi}$ and $M_{\psi}$ are bounded by constants depending only on $\Omega$, $\theta_{\mathcal{N}}$, $\theta_{\mathcal{F}}$, $\alpha$ and $\theta_{\psi}$, where

$$ \theta_{\psi} = \inf_{(\omega_{K_{\sigma}})_{K \in \mathcal{T}}, \sigma \in \mathcal{F}_K \in \mathcal{A}_{\psi}^{-1}(1)} \sup_{K \in \mathcal{T}} \sup_{\sigma \in \mathcal{F}_K} |\omega_{K, \sigma}| \quad \text{and} \quad \theta_{\mathcal{N}} = \max \left( \sup_{K \in \mathcal{T}} \sup_{\sigma \in \mathcal{F}_K} \frac{r_{\sigma}}{r_K}, \left( \inf_{K \in \mathcal{T}} \inf_{\sigma \in \mathcal{F}_K} \frac{r_{\sigma}}{r_K} \right)^{-1} \right). $$

*Proof.* Consider a function $\zeta \in C_c^{\infty}(\mathbb{R})$ taking positive values between 0 and $c$, with $\zeta(0) = c$, $\zeta(-1) = 0$, $\zeta(1) = 0$, $\zeta(\alpha) > 0$ and compactly supported in $]-1, 1[$. One can use for instance the function:

$$ \zeta(z) = \begin{cases} \frac{1}{ce^{(z^2 - 1)}} & \text{for } |z| < 1 \\ 0 & \text{for } |z| \ge 1, \end{cases} $$

where $c$ is such that $\int_{\mathbb{R}} \zeta = 1$. Another possible choice is given by:

$$ \zeta(z) = \begin{cases} \frac{\zeta_*(z+1)}{\zeta_*(z+1) + \zeta_*(-z)} & \text{for } z \le 0 \\ \frac{\zeta_*(1-z)}{\zeta_*(1-z) + \zeta_*(z)} & \text{for } z \ge 0 \end{cases} \quad \text{where} \quad \zeta_*(z) = \begin{cases} e^{-1/z} & \text{for } z > 0 \\ 0 & \text{otherwise.} \end{cases} $$

Then the family $(\psi_{\sigma})_{\sigma \in \mathcal{F}}$ defined by setting for all $\sigma \in \mathcal{F}$:

$$ \psi_{\sigma}(\mathbf{x}) = \frac{\zeta_{\sigma}(\mathbf{x})}{\sum_{\sigma' \in \mathcal{F}} \zeta_{\sigma'}(\mathbf{x})} \quad \text{where} \quad \zeta_{\sigma}(\mathbf{x}) = \zeta\left(\frac{|\mathbf{x} - \mathbf{x}_{\sigma}|}{r_{\sigma}}\right), \qquad (33) $$

is a partition of unity of $\Omega$, as we get by construction using hypothesis (32):

$$ \bar{\Omega} \subset \bigcup_{\sigma \in \mathcal{F}} \overline{B(\mathbf{x}_{\sigma}, r_{\sigma})} \quad \text{and} \quad \sum_{\sigma \in \mathcal{F}} \psi_{\sigma}(\mathbf{x}) = 1 \quad \forall \mathbf{x} \in \Omega \quad \text{and} \quad 0 \le \psi_{\sigma} \le 1 \quad \forall \sigma \in \mathcal{F}. $$

First, we establish that this family also satisfies $||\psi_{\sigma}||_{L^{\infty}(\Omega)} \le C_{\psi}$ for some constant $C_{\psi}$ independent on $r_{\sigma}$. To this end, first remark that $0 \le \zeta_{\sigma} \le c$: and that for any $\mathbf{x} \in \Omega$, using hypothesis (32) there exists $\sigma \in \mathcal{F}$ such that $\mathbf{x} \in B(\mathbf{x}_{\sigma}, \alpha r_{\sigma})$, thus we have:

$$ \sum_{\sigma' \in \mathcal{F}} \zeta_{\sigma'}(\mathbf{x}) \ge \zeta_{\sigma}(\mathbf{x}) \ge \zeta(\alpha) > 0. $$

Immediately, this gives:

$$ ||\psi_{\sigma}||_{L^{\infty}(\Omega)} \le \frac{||\zeta||_{L^{\infty}(\mathbb{R})}}{\zeta(\alpha)} = C_{\psi} $$

It is immediate to see that the associated family $(\psi_{K})_{K \in \mathcal{T}}$ defined by (30) using an optimal solution to (31) satisfies $r_K \le \rho_K \le 2\theta_{\mathcal{N}}r_K$, as for all $\sigma \in \mathcal{F}_K$, $r_{\sigma} \le \theta_{\mathcal{N}}r_K$, $r_{\sigma} \ge r_K$, $\rho_K \le 2 \max_{\sigma \in \mathcal{F}_K} r_{\sigma}$, and that:

$$ ||\psi_K||_{L^{\infty}(\Omega)} \le \theta_{\mathcal{F}}\theta_{\psi}C_{\psi}. $$

□

**Remark 3.4.** Given a family $(\psi_{\sigma})_{\sigma \in \mathcal{F}}$, one can compute in practice the quantity $\theta_{\psi}$, or at the very least an upper bound, simply solving (31) and looking eventually for its optimal solution. Remark that contrary to the practice of the network element method itself, establishing this theoretical bound would require some numerical
integration techniques, most probably a quadrature rule. In some sense, we see that the network element method evacuates numerical integration from practice, and confines it to the computation of some theoretical stability bounds. To see that $\theta_{\psi}$ can be controlled in a huge number of cases, consider the case where we construct the discretization network from a mesh, using the cell barycenters and diameters to define the point cloud and the connectivity, and choosing the interfaces to be the vertices of the mesh. Then, the $(m_K)_{K \in \mathcal{T}}$ can be taken equal to the cell measure. If the mesh admits a simplicial submesh based on its vertices, which is a common requirement, then using those simplices it is easy to construct a partition of unity satisfying (29), and the $(\omega_{K,\sigma})_{K \in \mathcal{T}, \sigma \in \mathcal{F}_K}$ can be defined by $\omega_{K,\sigma} = \frac{1}{m_{\sigma}} \int_K \psi_{\sigma}$. Thus we get $\theta_{\psi} \le 1$.

# 4 CONVERGENCE TO MINIMAL REGULARITY SOLUTIONS

We first recall the following strong consistency result established in [7]. Let $(\mathcal{N}, \mathcal{G})$ be an admissible discretization network and an associated admissible network geometry. For any $\varphi \in C_c^1(\mathbb{R}^d)$, there exists $C_{\varphi} > 0$ depending only on $\varphi$ and the quality parameters and independent on $h$ such that for any $K \in \mathcal{T}$ and any $\mathbf{x} \in B(\mathbf{x}_K, \xi_K)$ where $\xi_K \le \kappa_{\xi} r_K$ with $\kappa_{\xi} \ge 1$:

$$|\varphi(\mathbf{x}) - \mathcal{M}_K(\mathbf{D}_K(\varphi))| \le C_{\varphi} \kappa_{\xi} h_K, \quad (34)$$

while for any $\varphi \in C_c^2(\mathbb{R}^d)$, there exists another $C_{\varphi} > 0$ depending only on $\varphi$ and the quality parameters and independent on $h$ such that for any $K \in \mathcal{T}$ and any $\mathbf{x} \in B(\mathbf{x}_K, \xi_K)$:

$$|\nabla \varphi(\mathbf{x}) - \nabla_K(\mathbf{D}_K(\varphi))| \le C_{\varphi} \kappa_{\xi} (h_K + h_K^p) \quad \text{and} \quad |\varphi(\mathbf{x}) - \Pi_K(\mathbf{D}_K(\varphi))| \le C_{\varphi} \kappa_{\xi}^2 (h_K^2 + h_K^{p+1}). \quad (35)$$

For any $\Phi \in C_c^2(\mathbb{R}^d)^d$, there exists $C_{\Phi} > 0$ depending only on $\Phi$ and the quality parameters and independent on $h$ such that, for any $K \in \mathcal{N}$ and any $\mathbf{x} \in B(\mathbf{x}_K, \xi_K)$:

$$|\operatorname{div} \Phi(\mathbf{x}) - \mathcal{D}\mathcal{I}\mathcal{V}_K(\mathbf{D}_K(\Phi))| \le C_{\Phi} \kappa_{\xi} (h_K + h_K^p) \quad \text{with} \quad \mathcal{D}\mathcal{I}\mathcal{V}_K(\mathbf{D}_K(\Phi)) = \frac{1}{m_K} \sum_{\sigma \in \mathcal{F}_K} \eta_{K,\sigma}^T \Phi(\mathbf{x}_{\sigma}). \quad (36)$$

To establish convergence for solutions with minimal $H^1$ regularity, we will first need to establish consistency results for the reconstruction operators $\Pi_{\mathcal{T}}$, $\nabla_{\mathcal{T}}$ and $\Pi_{\mathcal{N}}$, as well as the stabilization bilinear form $s_h$. Next, following the usual finite volume, discontinuous Galerkin our Gradient discretization procedure (see [16], [13] or [15]), we will establish a refined weak consistency result for the discrete gradient $\nabla_{\mathcal{T}}$ applied to sequences bounded in the $||\cdot||_X$ norm, finally allowing to establish convergence of the network element method. Again, no originality is claimed regarding the general guidelines of the proof, which are completely classical, however every of these classical steps will require a careful treatment because of the quadrature family which takes care of the lack of a mesh.

## 4.1 Global consistency results for smooth functions

Building on the local consistency results, we can derive consistency estimates for reconstructed functions:

**Lemma 4.1 (Global approximation property).** Let $(\mathcal{N}, \mathcal{G})$ be an admissible discretization network and an associated admissible network geometry. For any $\varphi \in C_c^1(\mathbb{R}^d)$, there exists $C_{\varphi}$ depending only on $\varphi$, $\theta_{\mathcal{F}}$, $\theta_{\mathcal{I}}$, $\theta_{\mathcal{M}}$, $\theta_{\mathcal{A}}$, $M_{\psi}$, $\kappa_{\psi}$ and $\eta_{\psi}$ such that:

$$||\Pi_{\mathcal{T}}(\mathbf{D}(\varphi)) - \varphi||_{L^2(\Omega)} \le C_{\varphi}h,$$

while for any $\varphi \in C_c^2(\mathbb{R}^d)$, there exists $C_{\varphi}$ depending only on $\varphi$, $\theta_{\mathcal{F}}$, $\theta_{\mathcal{I}}$, $\theta_{\mathcal{M}}$, $\theta_{\mathcal{A}}$, $M_{\psi}$, $\kappa_{\psi}$ and $\eta_{\psi}$ such that:

$$||\nabla_{\mathcal{T}}(\mathbf{D}(\varphi)) - \nabla_{\mathcal{T}}(\varphi)||_{L^2(\Omega)^d} \le C_{\varphi}(h + h^p) \quad \text{and} \quad ||\Pi_{\mathcal{N}}(\mathbf{D}(\varphi)) - \varphi||_{L^2(\Omega)} \le C_{\varphi}(h^2 + h^{p+1}).$$
Proof. For any $\varphi \in C_c^1(\mathbb{R}^d)$, we have as $\sum_{L \in \mathcal{T}} \psi_L = 1$:

$$ \varphi - \Pi_{\mathcal{T}}(\mathbf{D}(\varphi)) = \varphi - \sum_{L \in \mathcal{T}} \psi_L \mathcal{M}_L(\mathbf{D}_L(\varphi)) = \varphi - \mathcal{M}_K(\mathbf{D}_K(\varphi)) + \sum_{L \in \mathcal{T}} \psi_L(\mathcal{M}_K(\mathbf{D}_K(\varphi)) - \mathcal{M}_L(\mathbf{D}_L(\varphi))), $$

Then, as the family $(\mathcal{B}_K)_{K \in \mathcal{T}}$ is an open cover of $\Omega$, we get that:

$$ \begin{aligned} \|\Pi_{\mathcal{T}}(\mathbf{D}(\varphi)) - \varphi\|_{L^2(\Omega)}^2 &\le \sum_{K \in \mathcal{T}} \|\Pi_{\mathcal{T}}(\mathbf{D}(\varphi)) - \varphi\|_{L^2(\mathcal{B}_K \cap \Omega)}^2 \\ &\le 2 \sum_{K \in \mathcal{T}} \|\varphi - \mathcal{M}_K(\mathbf{D}_K(\varphi))\|_{L^2(\mathcal{B}_K \cap \Omega)}^2 + 2 \sum_{K \in \mathcal{T}} \left\| \sum_{L \in \mathcal{T}} \psi_L(\mathcal{M}_K(\mathbf{D}_K(\varphi)) - \mathcal{M}_L(\mathbf{D}_L(\varphi))) \right\|_{L^2(\mathcal{B}_K \cap \Omega)}^2. \end{aligned} \quad (37) $$

For the first term, we obviously have using (34):

$$ \sum_{K \in \mathcal{T}} \|\varphi - \mathcal{M}_K(\mathbf{D}_K(\varphi))\|_{L^2(\mathcal{B}_K \cap \Omega)}^2 \le \sum_{K \in \mathcal{T}, \mathcal{B}_K \cap \text{supp } \varphi \ne \emptyset} \int_{\mathcal{B}_K \cap \Omega} C_{\varphi}^2 \kappa_{\psi}^2 h_K^2 \le \eta_{\psi} \left( |\text{supp } \varphi| + \frac{S_1^d}{2^d} \kappa_{\psi}^d h^d \right) C_{\varphi}^2 \kappa_{\psi}^2 h^2. $$

Notice then that, using hypotheses (24) we get:

$$ \left| \sum_{L \in \mathcal{T}} \psi_L(\mathcal{M}_K(\mathbf{D}_K(\varphi)) - \mathcal{M}_L(\mathbf{D}_L(\varphi))) \right| \le M_{\psi} \sum_{L \in \mathcal{T}} |\mathcal{M}_K(\mathbf{D}_K(\varphi)) - \mathcal{M}_L(\mathbf{D}_L(\varphi))| \chi_{\mathcal{B}_L}. $$

Thus we obtain using Cauchy-Schwarz inequality:

$$ \begin{aligned} \left\| \sum_{L \in \mathcal{T}} \psi_L(\mathcal{M}_K(\mathbf{D}_K(\varphi)) - \mathcal{M}_L(\mathbf{D}_L(\varphi))) \right\|_{L^2(\mathcal{B}_K \cap \Omega)}^2 &\le \eta_{\psi} M_{\psi}^2 \int_{\mathcal{B}_K \cap \Omega} \sum_{L \in \mathcal{T}} |\mathcal{M}_K(\mathbf{D}_K(\varphi)) - \mathcal{M}_L(\mathbf{D}_L(\varphi)))|^2 \chi_{\mathcal{B}_L} \\ &\le \eta_{\psi} M_{\psi}^2 \int_{\Omega} \sum_{L \in \mathcal{T}} |\mathcal{M}_K(\mathbf{D}_K(\varphi)) - \mathcal{M}_L(\mathbf{D}_L(\varphi)))|^2 \chi_{\mathcal{B}_L} \chi_{\mathcal{B}_K}. \end{aligned} $$

Then the second term of (37) is bounded by

$$ \begin{aligned} & \sum_{K \in \mathcal{T}} \left\| \sum_{L \in \mathcal{T}} \psi_L(\mathcal{M}_K(\mathbf{D}_K(\varphi)) - \mathcal{M}_L(\mathbf{D}_L(\varphi))) \right\|_{L^2(\mathcal{B}_K \cap \Omega)}^2 \\ &\le \eta_{\psi} M_{\psi}^2 \int_{\Omega} \sum_{K \in \mathcal{T}} \sum_{L \in \mathcal{T}} |\mathcal{M}_K(\mathbf{D}_K(\varphi)) - \mathcal{M}_L(\mathbf{D}_L(\varphi)))|^2 \chi_{\mathcal{B}_L} \chi_{\mathcal{B}_K}. \end{aligned} $$

Noticing that $\mathcal{M}_K(\mathbf{D}_K(\varphi)) = 0$ if $\mathcal{B}_K \cap \text{supp } \varphi = \emptyset$, that

$$ |\mathcal{M}_K(\mathbf{D}_K(\varphi)) - \mathcal{M}_L(\mathbf{D}_L(\varphi))| \le |\mathcal{M}_K(\mathbf{D}_K(\varphi)) - \varphi| + |\varphi - \mathcal{M}_L(\mathbf{D}_L(\varphi))|, $$

and that (34) gives for any $\boldsymbol{x} \in \mathcal{B}_K \cap \mathcal{B}_L$:

$$ |\mathcal{M}_K(\mathbf{D}_K(\varphi)) - \mathcal{M}_L(\mathbf{D}_L(\varphi))| \le C_{\varphi} \kappa_{\psi}(h_K + h_L), $$

we finally obtain that:

$$ \sum_{K \in \mathcal{T}} \left\| \sum_{L \in \mathcal{T}} \psi_L(\mathcal{M}_K(\mathbf{D}_K(\varphi)) - \mathcal{M}_L(\mathbf{D}_L(\varphi))) \right\|_{L^2(\mathcal{B}_K \cap \Omega)}^2 $$
$$
\begin{align*}
&\le \eta_{\psi} M_{\psi}^2 \kappa_{\psi}^2 C_{\varphi}^2 \int_{\Omega} \sum_{K \in \mathcal{T} \cap \text{supp } \varphi \ne \emptyset} \sum_{L \in \mathcal{T} \cap \text{supp } \varphi \ne \emptyset} (h_K + h_L)^2 \chi_{\mathcal{B}_L} \chi_{\mathcal{B}_K} \\
&\le 4 \eta_{\psi}^2 M_{\psi}^2 \kappa_{\psi}^2 C_{\varphi}^2 h^2 \int_{\Omega} \sum_{K \in \mathcal{T}, \mathcal{B}_K \cap \text{supp } \varphi \ne \emptyset} \chi_{\mathcal{B}_K} \le 4 \eta_{\psi}^2 M_{\psi}^2 \kappa_{\psi}^2 C_{\varphi}^2 h^2 \left( |\text{supp } \varphi| + \frac{S_1^d}{2^d} \kappa_{\psi}^d h^d \right),
\end{align*}
$$

which concludes the proof of the first estimate. The same proof leads to the second and third estimates, replacing (34) by respectively the first and second estimates of (35). □

**Lemma 4.2** (Stabilization consistency). Let $(\mathcal{N}, \mathcal{G})$ be an admissible discretization network and an associated admissible network geometry. For any $\varphi \in C_c^2(\mathbb{R}^d)$, there exists $C_\varphi$ depending only on $\varphi, \Omega, \theta_F, \theta_\Pi, \theta_M, \theta_A, \theta_T$ and $\eta_N$ such that:

$$
s_h(\mathbf{D}_h(\varphi), \mathbf{D}_h(\varphi)) = \sum_{K \in \mathcal{T}} s^K(\mathbf{D}_K(\varphi) - \Pi_K(\mathbf{D}_K(\varphi)), \mathbf{D}_K(\varphi) - \Pi_K(\mathbf{D}_K(\varphi))) \le S^* C_\varphi (h + h^p)^2.
$$

*Proof.* Using (35), there exists a constant $C_\varphi$ such that for any $\sigma \in \mathcal{F}_K$:

$$
|\varphi(\mathbf{x}_{\sigma}) - \Pi_{K}(\mathbf{D}_{K}(\varphi))(\mathbf{x}_{\sigma})| \leq C_{\varphi}(h_{K}^{2} + h_{K}^{p+1}).
$$

Thus, we get:

$$
\begin{align*}
s_h(\mathbf{D}_h(\varphi), \mathbf{D}_h(\varphi)) &\le \theta_F S^* C_\varphi^2 \sum_{K \in \mathcal{T}, B_K \cap \text{supp } \varphi \ne \emptyset} m_K h_K^{-2} (h_K^2 + h_K^{(p+1)})^2 \\
&\le \theta_T \theta_F S^* C_\varphi^2 \sum_{K \in \mathcal{T}, B_K \cap \text{supp } \varphi \ne \emptyset} (h_K^2 + h_K^{(p+1)})^2 \int_{B_K \cap \Omega} 1 \le \eta_N \theta_T \theta_F S^* C_\varphi^2 \left( |\text{supp } \varphi| + \frac{S_1^d}{2^d} h^d \right) (h+h^p)^2.
\end{align*}
$$

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

## 4.2 Weak consistency of the reconstructed gradient $\nabla_T$

In the remaining of this section, we will consider a family $(\mathcal{N}_h, \mathcal{G}_h)_{h \in \mathcal{H}}$ of admissible discretization networks and associated admissible network geometries indexed by $h \in \mathcal{H}$, where $\mathcal{H}$ is a bounded at most countable subset of $\mathbb{R}^+$ with $0 \in \bar{\mathcal{H}}$. We will consider the case where there exists constants $\theta > 0$, $\eta > 0$, $M > 0$ and $\kappa > 0$ independent on $h$ and quadrature families $(\psi_h)_{h \in \mathcal{H}}$ associated with each $(\mathcal{N}_h, \mathcal{G}_h)$, such that $\forall h \in \mathcal{H}$, $\max(\theta_{\mathcal{F}_h}, \theta_{\Pi_h}, \theta_{\mathcal{A}_h}, \theta_{\mathcal{M}_h}) \le \theta$, $\max(\eta_{\mathcal{N}_h}, \eta_{\psi_h}) \le \eta$, $\kappa_\psi \le \kappa$ and $M_{\psi_h} \le M$. We call such a family an admissible discretization family. We have the following weak convergence property for discrete solutions on admissible discretization families:

**Lemma 4.3** (Weak consistency). Let $(\mathcal{N}_h, \mathcal{G}_h)_{h \in \mathcal{H}}$ be an admissible discretization family. Let $(U_h)_{h \in \mathcal{H}}$ be a family such that

• $U_h \in X_{\mathcal{N}_h, 0}$ for any $h \in \mathcal{H}$,

• There exists $C > 0$ independent on $h$ such that $||U_h||_{X_h} \le C$ for all $h \in \mathcal{H}$,

• There exists $u \in L^2(\Omega)$ such that $\Pi_{\mathcal{T}_h}(U_h) \to u$ weakly in $L^2(\Omega)$ when $h \to 0$.

Then we have, for any $\Phi \in C_c^\infty(\Omega)^d$:

$$
\int_{\Omega} u \operatorname{div} \Phi + \int_{\Omega} \nabla_{\mathcal{T}_h}(U_h) \cdot \Phi \to 0 \quad \text{when } h \to 0, \qquad (38)
$$

and also $u \in H_0^1(\Omega)$ and $\nabla_{\mathcal{T}_h}(U_h) \to \nabla u$ weakly in $L^2(\Omega)^d$ when $h \to 0$.
*Proof.* Let us start by extending $\Pi_{\mathcal{T}_h}(U_h)$ and $\nabla_{\mathcal{T}_h}(U_h)$ by 0 outside of $\Omega$, and denote those extensions respectively $\tilde{\Pi}_{\mathcal{T}_h}(U_h)$ and $\tilde{\nabla}_{\mathcal{T}_h}(U_h)$. By virtue of lemma 3.1, we know that up to a subsequence there exists $G \in L^2(\mathbb{R}^d)^d$ such that $\tilde{\nabla}_{\mathcal{T}_h}(U_h) \to G$ weakly in $L^2(\mathbb{R}^d)^d$ when $h \to 0$ and that $\tilde{\Pi}_{\mathcal{T}_h}(U_h) \to \tilde{u}$ weakly in $L^2(\mathbb{R}^d)$ when $h \to 0$, where $\tilde{u}$ denotes the extension by zero of $u$ outside $\Omega$. Then recalling that:

$$ \mathcal{D}\mathcal{I}\mathcal{V}_K(\mathcal{D}_K(\Phi)) = \frac{1}{m_K} \sum_{\sigma \in \mathcal{F}_K} \eta_{K,\sigma}^T \Phi(\mathbf{x}_{\sigma}), $$

notice that, for any $\Phi \in C_c^\infty(\mathbb{R}^d)^d$:

$$ \int_{\mathbb{R}^d} \tilde{\Pi}_{\mathcal{T}_h}(U_h) \operatorname{div} \Phi = \int_{\Omega} \left( \sum_{K \in \mathcal{T}_h} \psi_K \mathcal{M}_K(U_h) \right) \operatorname{div} \Phi = \sum_{K \in \mathcal{T}_h} \int_{\mathcal{B}_K \cap \Omega} \psi_K \mathcal{M}_K(U_h) \operatorname{div} \Phi, $$

which leads to:

$$ \int_{\mathbb{R}^d} \tilde{\Pi}_{\mathcal{T}_h}(U_h) \operatorname{div} \Phi = \sum_{K \in \mathcal{T}_h} \int_{\mathcal{B}_K \cap \Omega} \psi_K \mathcal{M}_K(U_h) \mathcal{D}\mathcal{I}\mathcal{V}_K(\mathcal{D}_K(\Phi)) \\ + \sum_{K \in \mathcal{T}_h} \int_{\mathcal{B}_K \cap \Omega} \psi_K \mathcal{M}_K(U_h) (\operatorname{div} \Phi - \mathcal{D}\mathcal{I}\mathcal{V}_K(\mathcal{D}_K(\Phi))). $$

Let us denote the previous identity $T_1+T_2$, with obvious notations. Focusing on $T_1$, we have using the definition of $\mathcal{D}\mathcal{I}\mathcal{V}_K$:

$$ T_1 = \sum_{K \in \mathcal{T}_h} \sum_{\sigma \in \mathcal{F}_K} \mathcal{M}_K(U_h) \eta_{K,\sigma} \cdot \Phi(\mathbf{x}_{\sigma}). $$

Then, using the approximate geometrical conservation property i.e. $\sum_{K \in \mathcal{T}_h} \eta_{K,\sigma} = \varepsilon_\sigma$ for $\sigma \in \mathcal{F}_{h,int}$ and the fact that $u_\sigma = 0$ for $\sigma \in \mathcal{F}_{h ext}$ as $U_h \in X_{N_h,0}$, this rewrites:

$$ \begin{aligned} T_1 &= \sum_{K \in \mathcal{T}_h} \sum_{\sigma \in \mathcal{F}_K} (\mathcal{M}_K(U_h) - u_\sigma) \eta_{K,\sigma} \cdot \Phi(\mathbf{x}_\sigma) + \sum_{K \in \mathcal{F}_{h, int}} \varepsilon_\sigma \cdot \Phi(\mathbf{x}_\sigma) u_\sigma \\ &= \sum_{K \in \mathcal{T}_h} \sum_{\sigma \in \mathcal{F}_K} (\mathcal{M}_K(U_h) - u_\sigma) \eta_{K,\sigma} \cdot \Phi(\mathbf{x}_K) + \sum_{K \in \mathcal{T}_h} \sum_{\sigma \in \mathcal{F}_K} (\mathcal{M}_K(U_h) - u_\sigma) \eta_{K,\sigma} \cdot (\Phi(\mathbf{x}_\sigma) - \Phi(\mathbf{x}_K)) \\ &\quad + \sum_{K \in \mathcal{F}_{h, int}} \varepsilon_\sigma \cdot \Phi(\mathbf{x}_\sigma) u_\sigma = T_{1,1} + T_{1,2} + T_{1,3}, \end{aligned} $$

with obvious notations. Then, recalling that:

$$ \nabla_K (U_K) = \frac{1}{m_K} \sum_{\sigma \in \mathcal{F}_K} u_\sigma \eta_{K,\sigma}, $$

the term $T_{1,1}$ rewrites:

$$ \begin{aligned} T_{1,1} &= - \sum_{K \in \mathcal{T}_h} m_K \nabla_K (U_h) \cdot \Phi(\mathbf{x}_K) = - \sum_{K \in \mathcal{T}_h} \int_{\mathcal{B}_K \cap \Omega} \psi_K \nabla_K (U_h) \cdot \Phi(\mathbf{x}_K) \\ &= - \int_{\mathbb{R}^d} \tilde{\nabla}_{\mathcal{T}_h}(U_h) \cdot \Phi(\mathbf{x}) - \sum_{K \in \mathcal{T}_h} \int_{\mathcal{B}_K \cap \Omega} \psi_K \nabla_K (U_h) \cdot (\Phi(\mathbf{x}_K) - \Phi(\mathbf{x})), \end{aligned} $$

which rewrites $T_{1,1} = T_{1,1,1}+T_{1,1,2}$ with obvious notations. For the second term, we have using Cauchy-Schwarz inequality:

$$ |T_{1,1,2}| = \left| \sum_{K \in \mathcal{T}_h} \int_{\mathcal{B}_K \cap \Omega} \psi_K \nabla_K (U_h) \cdot (\Phi(\mathbf{x}_K) - \Phi(\mathbf{x})) \right| $$
$$
\le h ||\nabla \Phi||_{L^{\infty}(\mathbb{R}^d)^{d \times d}} \left( \sum_{K \in \mathcal{T}_h} \int_{\mathcal{B}_K \cap \Omega} |\psi_K||\nabla_K(\mathbf{U}_h)|^2 \right)^{1/2} \left( \sum_{K \in \mathcal{T}_h, \mathcal{B}_K \cap \text{supp } \Phi \ne \emptyset} \int_{\mathcal{B}_K \cap \Omega} |\psi_K| \right)^{1/2} \\
\le C_{\nabla,h} h ||\mathbf{U}_h||_{X_h} ||\nabla \Phi||_{L^{\infty}(\mathbb{R}^d)^{d \times d}} \theta_{\mathcal{T}_h}^{1/2} M^{1/2} \left( \frac{S_1^d \kappa_{\psi}^d \delta^d}{|C(0, \tau, 1)|} \right)^{1/2} \left( \sum_{K \in \mathcal{T}_h, \mathcal{B}_K \cap \text{supp } \Phi \ne \emptyset} \int_{\mathcal{B}_K \cap \Omega} |\psi_K| \right)^{1/2},
$$

using (17) and (28), and thus as $C_{\nabla,h}$ is bounded by some $C_{\nabla} > 0$ independent on $h$:

$$
|T_{1,1,2}| \le \eta^{1/2} \theta^{1/2} M C_{\nabla} \left( \frac{S_1^d \kappa_{\psi}^d \delta^d}{|C(0, \tau, 1)|} \right)^{1/2} \left( |\text{supp } \Phi| + \frac{S_1^d}{2^d} \kappa_{\psi}^d h^d \right)^{1/2} h ||\mathbf{U}_h||_{X_h} ||\nabla \Phi||_{L^{\infty}(\mathbb{R}^d)^{d \times d}}.
$$

Now it just remains to bound $T_2$, $T_{1,2}$, $T_{1,3}$ and $T_{1,4}$. For $T_2$, using (36) and Cauchy-Schwarz inequality we immediately get:

$$
|T_2| \le C_{\Phi} \kappa_{\psi} (h + h^p) \left( \sum_{K \in \mathcal{T}_h} \int_{\mathcal{B}_K \cap \Omega} |\psi_K| |\mathcal{M}_K(\mathbf{U}_h)|^2 \right)^{1/2} \left( \sum_{K \in \mathcal{T}_h, \mathcal{B}_K \cap \text{supp } \Phi \ne \emptyset} \int_{\mathcal{B}_K \cap \Omega} |\psi_K| \right)^{1/2} \\
\le \eta^{1/2} C_{\Phi} \kappa_{\psi} M \left( \frac{S_1^d \kappa_{\psi}^d \delta^d}{|C(0, \tau, 1)|} \right)^{1/2} \left( |\text{supp } \Phi| + \frac{S_1^d}{2^d} \kappa_{\psi}^d h^d \right)^{1/2} (h + h^p) ||\mathbf{U}_h||_{X_h}.
$$

Next, for $T_{1,2}$ we get using Cauchy-Schwarz inequality once again:

$$
|T_{1,2}| \le h ||\nabla \Phi||_{L^{\infty}(\mathbb{R}^d)^{d \times d}} \left( \sum_{K \in \mathcal{T}_h} \sum_{\sigma \in \mathcal{F}_K} \int_{\mathcal{B}_K \cap \Omega} h_K^{-2} |\psi_K| |\mathcal{M}_K(\mathbf{U}_h) - u_\sigma|^2 \left| \frac{h_K \eta_{K,\sigma}}{m_K} \right|^2 \right)^{1/2} \\
\times \left( \sum_{K \in \mathcal{T}_h, \mathcal{B}_K \cap \text{supp } \Phi \neq \emptyset} \int_{\mathcal{B}_K \cap \Omega} |\psi_K| \right)^{1/2} \\
\le \eta^{1/2} \theta^{3/2} M \left( \frac{S_1^d \kappa_{\psi}^d \delta^d}{|C(0, \tau, 1)|} \right)^{1/2} \left( |\text{supp } \Phi| + \frac{S_1^d}{2^d} \kappa_{\psi}^d h^d \right)^{1/2} h ||\mathbf{U}_h||_{X_h} ||\nabla \Phi||_{L^{\infty}(\mathbb{R}^d)^{d \times d}}.
$$

Finally, for $T_{1,3}$, we have using (9) and Cauchy-Schwarz inequality:

$$
|T_{1,3}| \le \theta_A ||\Phi||_{L^\infty(\mathbb{R}^d)^d} \sum_{K \in \mathcal{T}_h} \sum_{\Phi \ne \emptyset} m_K h_K^p |u_\sigma| \\
\le \theta_A ||\Phi||_{L^\infty(\mathbb{R}^d)^d} h^p \left( \sum_{K \in \mathcal{T}_h} \sum_{\Phi \ne \emptyset} m_K \right)^{1/2} \left( \sum_{K \in \mathcal{T}_h} \sum_{\sigma \in \mathcal{F}_K} m_K |u_\sigma|^2 \right)^{1/2} \\
\le \theta_A ||\Phi||_{L^\infty(\mathbb{R}^d)^d} \theta_{\mathcal{F}}^{1/2} M_{\psi}^{1/2} \eta_{\psi}^{1/2} h^p \left( |\text{supp } \Phi| + \frac{S_1^d}{2^d} \kappa_{\psi}^d h^d \right)^{1/2} \left( \sum_{K \in \mathcal{T}_h} \sum_{\sigma \in \mathcal{F}_K} m_K |u_\sigma|^2 \right)^{1/2}.
$$

Then notice that:

$$
\sum_{K \in \mathcal{T}_h} \sum_{\sigma \in \mathcal{F}_K} m_K |u_\sigma|^2 \le 2 \sum_{K \in \mathcal{T}_h} \sum_{\sigma \in \mathcal{F}_K} m_K |u_\sigma - \mathcal{M}_K(\mathbf{U}_K)|^2 + 2 \sum_{K \in \mathcal{T}_h} \sum_{\sigma \in \mathcal{F}_K} m_K |\mathcal{M}_K(\mathbf{U}_K)|^2
$$
which immediately leads to:

$$|T_{1,3}| \le 2\theta_A ||\Phi||_{L^\infty(\mathbb{R}^d)} \theta_{\mathcal{F}}^{1/2} M_{\psi}^{1/2} \eta_{\psi}^{1/2} h^p \left( |\text{supp } \Phi| + \frac{S_1^d}{2^d} \kappa_{\psi}^d h^d \right)^{1/2} (||U_h||_{0,h}^2 + h^2 |U_h|_{1,h}^2)^{1/2}.$$

Consequently, using our hypothesis, for any $\Phi \in C_c^\infty(\mathbb{R}^d)^d$ there exists $C_\Phi(M, \theta, \eta, \kappa) > 0$ depending on $\theta, \eta, \kappa$ and $\Phi$ and the bound on the family $(U_h)_{h \in \mathcal{H}}$ such that, as soon as $h \le 1$:

$$\left| \int_{\mathbb{R}^d} \tilde{\Pi}_{\tau_h}(U_h) \operatorname{div} \Phi + \int_{\Omega} \tilde{\nabla}_{\tau_h}(U_h) \cdot \Phi \right| \le C_{\Phi}(M, \theta, \eta, \kappa)h,$$

which directly implies (38) and leads to $\nabla_{\tau_h}(U_h) \to \nabla \tilde{u}$ weakly in $L^2(\mathbb{R})^d$ when $h \to 0$. Thus $G = \nabla \tilde{u}$, which implies $\tilde{u} \in H^1(\mathbb{R}^d)$ and thus $u \in H_0^1(\Omega)$ and concludes the proof. $\square$

## 4.3 Convergence result for minimal regularity solutions

To conclude, we will need the following useful link between the discrete norm $|| \cdot ||_{X_h}$ and classical norms for regular functions:

**Lemma 4.4.** Let $(N_h, G_h)_{h \in \mathcal{H}}$ be an admissible discretization family. For any $\varphi \in C_c^\infty(\mathbb{R}^d)$, there exists $C > 0$ depending only on $\varphi, \theta$ and $\eta$ such that

$$||\mathcal{D}_h(\varphi)||_{X_h} \le C \left( |\text{supp } \varphi| + \frac{S_1^d}{2^d} h^d \right)^{1/2} \left( ||\varphi||_{L^\infty(\mathbb{R}^d)} + \sup_{|\alpha|=1} ||\partial^\alpha \varphi||_{L^\infty(\mathbb{R}^d)} \right).$$

*Proof.* First notice that by definition, we have:

$$||\mathcal{D}_h(\varphi)||_{X_h}^2 = \sum_{K \in \mathcal{T}_h, \text{supp } \varphi \cap B_K \neq \emptyset} m_K |\mathcal{M}_K(\mathcal{D}_K(\varphi))|^2 + \sum_{K \in \mathcal{T}_h, \text{supp } \varphi \cap B_K \neq \emptyset} \sum_{\sigma \in \mathcal{F}_K} m_K h_K^{-2} |\mathcal{M}_K(\mathcal{D}_K(\varphi)) - \varphi(\mathbf{x}_\sigma)|^2$$

Immediately, using Taylor's expansion, we know that:

$$\varphi(\mathbf{x}_{\sigma'}) = \varphi(\mathbf{x}_{\sigma}) + \sum_{|\alpha|=1} (\mathbf{x}_{\sigma'} - \mathbf{x}_{\sigma})^{\alpha} \int_{0}^{1} \partial^{\alpha} \varphi(\mathbf{x} + t(\mathbf{x}_{\sigma'} - \mathbf{x}_{\sigma})),$$

with:

$$\left| \sum_{|\alpha|=1} (\mathbf{x}_{\sigma'} - \mathbf{x}_{\sigma})^{\alpha} \int_{0}^{1} \partial^{\alpha} \varphi(\mathbf{x} + t(\mathbf{x}_{\sigma'} - \mathbf{x}_{\sigma})) \right| \le dh_K \sup_{|\alpha|=1} ||\partial^{\alpha} \varphi||_{L^{\infty}(\mathbb{R}^d)}.$$

Then, as $\sum_{\sigma \in \mathcal{F}_K} \gamma_{K,\sigma} = 1$, and as $\mathcal{M}_K(\mathcal{D}_K(\varphi)) = \sum_{\sigma' \in \mathcal{F}_K} \gamma_{K,\sigma'} \varphi(\mathbf{x}_{\sigma'})$ we immediately get that:

$$\begin{aligned} ||\mathcal{D}_h(\varphi)||_{X_h}^2 &\le \theta_{\mathcal{M}_h}^2 \theta_{\mathcal{F}_h}^2 \left( \sum_{K \in \mathcal{T}_h, \text{supp } \varphi \cap B_K \neq \emptyset} m_K \right) \left( ||\varphi||_{L^\infty(\mathbb{R}^d)}^2 + d^2 \theta_{\mathcal{F}_h} \sup_{|\alpha|=1} ||\partial^\alpha \varphi||_{L^\infty(\mathbb{R}^d)}^2 \right) \\ &\le \eta_{\mathcal{N}_h} \theta_{\mathcal{M}_h}^2 \theta_{\mathcal{F}_h}^2 \theta_{\tau_h} \left( |\text{supp } \varphi| + \frac{S_1^d}{2^d} h^d \right) \left( ||\varphi||_{L^\infty(\mathbb{R}^d)}^2 + d^2 \theta_{\mathcal{F}_h} \sup_{|\alpha|=1} ||\partial^\alpha \varphi||_{L^\infty(\mathbb{R}^d)}^2 \right), \end{aligned}$$

which concludes the proof. $\square$
Gathering lemma 4.1, 4.2, 4.3 and 4.4 and the stability (19) and coercivity (20) results of [7] up to some reformulation of our approximation extensively using the quadrature family, one could now recast our results inside the Gradient discretization framework (see [15]) and thus automatically obtain the convergence to minimal regularity solutions. However, this would require introducing many new concepts and notations, which is the reason why (and also for the sake of completeness) we provide a basic finite volume like proof of the convergence of the method:

**Proposition 4.5** (Convergence). Let $(\mathcal{N}_h, \mathcal{G}_h)_{h \in \mathcal{H}}$ be an admissible discretization family, and let $(\mathbf{U}_h)_{h \in \mathcal{H}}$ be the solution of the associated problem (16) for each $h \in \mathcal{H}$. Assume that:

$$ \sum_{K \in \mathcal{T}} \int_{\mathcal{B}_K \cap \Omega} |f_K - f|^2 \cdot \to 0 $$

Then $\Pi_{\mathcal{T}_h}(\mathbf{U}_h)$ strongly converges in $L^2(\Omega)$ to the solution $u$ of (11) when $h \to 0$. Moreover $\nabla_{\mathcal{T}_h}(\mathbf{U}_h)$ strongly converges in $L^2(\Omega)^d$ to $\nabla u$ when $h \to 0$

*Proof.* In the following, $C > 0$ denotes a constant independent on $h$ whose value can change from one line to another. From the hypothesis, we know that there exists some $h_0$ small enough such that:

$$ \sum_{K \in \mathcal{T}} \int_{\mathcal{B}_K \cap \Omega} |f_K - f|^2 \le ||f||_{L^2(\Omega)} $$

holds for any $h \le h_0$. Then, using (19)-(20)-(23), it is clear that

$$ ||\mathbf{U}_h||_{X_h} \le C \left( \sum_{K \in \mathcal{T}_h} m_K |f_K|^2 \right)^{\frac{1}{2}} \le C ||f||_{L^2(\Omega)} $$

Using lemma 3.1, we get $||\Pi_{\mathcal{T}_h}(\mathbf{U}_h)||_{L^2(\Omega)} \le C||f||_{L^2(\Omega)}$ and $||\nabla_{\mathcal{T}_h}(\mathbf{U}_h)||_{L^2(\Omega)^d} \le C||f||_{L^2(\Omega)}$. Consequently, up to subsequence, there exists $u \in L^2(\Omega)$ such that the hypothesis of lemma 4.3 are satisfied. Thus $u \in H_0^1(\Omega)$ and $\nabla_{\mathcal{T}_h}(\mathbf{U}_h) \to \nabla u$ weakly in $L^2(\Omega)^d$. Next, for any $\varphi \in C_c^\infty(\Omega)$, we have:

$$ a_h(\mathbf{U}_h, \mathbf{D}_h(\varphi)) = l_h(\mathbf{D}_h(\varphi)). $$

Remark that:

$$ a_h(\mathbf{U}_h, \mathbf{D}_h(\varphi)) = \int_{\Omega} \nabla_{\mathcal{T}_h}(\mathbf{U}_h) \cdot \nabla \varphi + s^h(\mathbf{U}_h, \mathbf{D}_h(\varphi)) + \sum_{K \in \mathcal{T}_h} \int_{\mathcal{B}_K \cap \Omega} \psi_K \nabla \Pi_K(\mathbf{U}_h) \cdot (\nabla \Pi_K(\mathbf{D}_K(\varphi)) - \nabla \varphi). $$

Using Cauchy-Schwarz inequality and lemma 4.2 we get:

$$ |s^h(\mathbf{U}_h, \mathbf{D}_h(\varphi))| \le |s^h(\mathbf{U}_h, \mathbf{U}_h)|^{1/2} |s^h(\mathbf{D}_h(\varphi), \mathbf{D}_h(\varphi))|^{1/2} \le C(h^2 + h^{p+1}) ||\mathbf{U}_h||_{X_h} \le C(h^2 + h^{p+1}) ||f||_{L^2(\Omega)}. $$

Estimate (35) gives:

$$ \left| \sum_{K \in \mathcal{T}_h} \int_{\mathcal{B}_K \cap \Omega} \psi_K \nabla \Pi_K(\mathbf{U}_h) \cdot (\nabla \Pi_K(\mathbf{D}_K(\varphi)) - \nabla \varphi) \right| \le Ch \left( |\text{supp } \varphi| + \frac{S_1^d}{2^d} \kappa_{\psi}^d h^d \right) ||\mathbf{U}_h||_{X_h}. $$

Thus using lemma 4.3 we get that:

$$ a_h(\mathbf{U}_h, \mathbf{D}_h(\varphi)) \to \int_{\Omega} \nabla u \cdot \nabla \varphi \quad \text{when } h \to 0. $$
Finally, for the right-hand side we have that:

$$l_h(\mathbf{D}_h(\varphi)) = \int_{\Omega} f \left( \sum_{K \in \mathcal{T}_h} \psi_K \mathcal{M}_K(\mathbf{D}_K(\varphi)) \right) + \sum_{K \in \mathcal{T}_h} \int_{\mathcal{B}_K \cap \Omega} \psi_K (f - f_K) \mathcal{M}_K(\mathbf{D}_K(\varphi)).$$

Using Cauchy-Schwarz inequality, we can bound the second term by:

$$\left| \sum_{K \in \mathcal{T}_h} \int_{\mathcal{B}_K} \psi_K (f - f_K) \mathcal{M}_K(\mathbf{D}_K(\varphi)) \right| \le C ||\varphi||_{L^\infty(\Omega)} \left( |\text{supp } \varphi| + \frac{S_1^d}{2^d} \kappa_\psi^d h^d \right)^{1/2} \left( \sum_{K \in \mathcal{T}_h} \int_{\mathcal{B}_K \cap \Omega} |f_K - f|^2 \right)^{1/2},$$

which immediately leads to:

$$l_h(\mathbf{D}_h(\varphi)) \to \int_{\Omega} f\varphi \quad \text{when } h \to 0.$$

Using the density of $C_c^\infty(\Omega)$ in $H_0^1(\Omega)$, this concludes the proof of the fact that for all $v \in H_0^1(\Omega)$:

$$\int_{\Omega} \nabla u \cdot \nabla v = \int_{\Omega} fv.$$

Next, let $\varphi \in C_c^\infty(\Omega)$. Using the triangular inequality, we have:

$$\int_{\Omega} |\nabla \tau_h (U_h) - \nabla u|^2 \le \int_{\Omega} |\nabla \tau_h (U_h) - \nabla \tau_h (\mathbf{D}_h(\varphi))|^2 + \int_{\Omega} |\nabla \tau_h (\mathbf{D}_h(\varphi)) - \nabla \varphi|^2 + \int_{\Omega} |\nabla \varphi - \nabla u|^2.$$

From lemma 4.1, we know that:

$$\int_{\Omega} |\nabla \tau_h (\mathbf{D}_h(\varphi)) - \nabla \varphi|^2 \to 0 \quad \text{when} \quad h \to 0.$$

Next, we have using lemma 3.1 and (20)

$$\int_{\Omega} |\nabla \tau_h (U_h - \mathbf{D}_h(\varphi))|^2 \le C ||U_h - \mathbf{D}_h(\varphi)||_{X_h}^2 \le C a_h (U_h - \mathbf{D}_h(\varphi), U_h - \mathbf{D}_h(\varphi)).$$

Immediately, we see that:

$$a_h(U_h - \mathbf{D}_h(\varphi), U_h - \mathbf{D}_h(\varphi)) = a_h(U_h, U_h) - 2a_h(U_h, \mathbf{D}_h(\varphi)) + a_h(\mathbf{D}_h(\varphi), \mathbf{D}_h(\varphi)).$$

We have already seen that:

$$a_h(U_h, \mathbf{D}_h(\varphi)) \to \int_{\Omega} \nabla u \cdot \nabla \varphi \quad \text{when } h \to 0.$$

From lemma 4.4, we know that the family $(||\mathbf{D}_h\varphi||_{X_h})_{h \in \mathcal{H}}$ is also bounded. Thus, we can apply the above reasoning to $\varphi$ and we get:

$$a_h(\mathbf{D}_h(\varphi), \mathbf{D}_h(\varphi)) \to \int_{\Omega} \nabla \varphi \cdot \nabla \varphi \quad \text{when } h \to 0.$$

Finally, we have as $U_h$ is solution of (16) that $a_h(U_h, U_h) = l_h(U_h)$. However

$$l_h(U_h) = \int_{\Omega} f \Pi_{\tau_h}(U_h) + \sum_{K \in \mathcal{T}_h} \int_{\mathcal{B}_K \cap \Omega} \psi_K (f - f_K) \mathcal{M}_K(U_K),$$
and as Cauchy-Schwarz inequality gives:

$$ \left| \sum_{K \in \mathcal{T}_h} \int_{B_K} \psi_K(f - f_K) \mathcal{M}_K(U_K) \right| \le C ||U_h||_0 \left( \sum_{K \in \mathcal{T}_h} \int_{B_K \cap \Omega} |f_K - f|^2 \right)^{1/2}, $$

and $\Pi_{\mathcal{T}_h}(U_h) \to u$ weakly in $L^2(\Omega)$, we deduce that:

$$ l_h(U_h) \to \int_{\Omega} fu = \int_{\Omega} |\nabla u|^2 \quad \text{when} \quad h \to 0. $$

Gathering all the previous results, we get that

$$ \limsup_{h \to 0} \int_{\Omega} |\nabla_{\mathcal{T}_h}(U_h - \mathcal{D}_h(\varphi))|^2 \le C ||u - \varphi||_{H^1(\Omega)}^2. $$

By density of $C^\infty(\Omega)$ in $H_0^1(\Omega)$, for any $\varepsilon > 0$ we can choose $\varphi$ such that

$$ ||\varphi - u||_{H^1(\Omega)}^2 \le \min \left( \frac{\varepsilon}{3}, \frac{\varepsilon}{3C} \right), $$

while for this fixed $\varphi$, we can choose $h_\varepsilon > 0$ such that for any $h \in \mathcal{H}$, $h \le h_\varepsilon$

$$ \int_{\Omega} |\nabla_{\mathcal{T}_h}(U_h - \mathcal{D}_h(\varphi))|^2 \le \frac{\varepsilon}{3} \quad \text{and} \quad \int_{\Omega} |\nabla_{\mathcal{T}_h}(\mathcal{D}_h(\varphi)) - \nabla \varphi|^2 \le \frac{\varepsilon}{3}. $$

Thus for any $\varepsilon > 0$, there exists $h_\varepsilon > 0$ such that for any $h \in \mathcal{H}$, $h \le h_\varepsilon$ we have:

$$ \int_{\Omega} |\nabla_{\mathcal{T}_h}(U_h) - \nabla u|^2 \le \varepsilon, $$

which implies that $\nabla_{\mathcal{T}_h}(U_h) \to \nabla u$ strongly in $L^2(\Omega)$ when $h \to 0$. Proceeding in exactly the same way, we show that $\Pi_{\mathcal{T}_h}(U_h) \to u$ strongly in $L^2(\Omega)$, which concludes the proof. $\square$

# 5 ERROR ESTIMATES FOR REGULAR PROBLEMS

The aim of this section is to provide explicit convergence rates when the solution $u$ of (11) is regular enough. To establish error estimates for solutions with Sobolev regularity, we will first need to refine the local consistency results for the NEM operators $\mathcal{M}_K$, $\nabla_K$ and $\Pi_K$. Then, building on those local results we will establish global consistency results similar to lemma (4.1) and (4.2), for functions with Sobolev regularity only. Finally, using those consistency results we will be able to establish our error estimates.

We recall the following useful result on Riesz potentials (see [4]): let $B$ be a ball of $\mathbb{R}^d$ of radius $\rho$, $f \in L^p(B)$, $p \ge 1$ and $m \ge 1$. Let $g$ be defined by:

$$ g(x) = \int_B |x - z|^{m-d} |f(z)| dz. $$

Then, there exists $C_{m,d} > 0$ depending only on $m$ and $d$ such that

$$ ||g||_{L^p(B)} \le C_{m,d} \rho^m ||f||_{L^p(B)}. \quad (39) $$

Another useful remark is the following: as $\Omega$ is Lipschitz, using Stein's extension theorem for any $k \ge 0$ and any $v \in H^k(\Omega)$, we have:

$$ \sum_{K \in \mathcal{T}} |v|_{H^k(B_K)}^2 \le \sum_{K \in \mathcal{T}} |Ev|_{H^k(B_K)}^2 = \sum_{K \in \mathcal{T}} \int_{\mathbb{R}^d} |Ev|^2 \chi_{B_K} = \int_{\mathbb{R}^d} |Ev|^2 \left( \sum_{K \in \mathcal{T}} \chi_K \right) \le \eta_\psi |Ev|_{H^k(\mathbb{R}^d)}^2, $$
and thus

$$ \sum_{K \in \mathcal{T}} |v|_{H^k(\mathcal{B}_K)}^2 \le \eta_\psi C_{E,k}^2 |v|_{H^k(\Omega)}^2, \quad (40) $$

where we recall that

$$ \mathcal{T}_{\mathbf{x}}^{\mathcal{B}} = \{K \in \mathcal{T} \mid \mathbf{x} \in \mathcal{B}_K\} \quad \text{and} \quad \eta_{\psi} = \sup_{\mathbf{x} \in \mathbb{R}^d} \mathrm{card}(\mathcal{T}_{\mathbf{x}}^{\mathcal{B}}). $$

## 5.1 Network element interpolation

For continuous functions, we have already defined degrees of freedom through the operator $\mathcal{D}: C^0(\bar{\Omega}) \to X_N$. To handle the case of functions that only belong to a Sobolev space, we define another operator, clearly inspired by the usual Clément finite element interpolant. To any $\sigma \in \mathcal{F}$, we associate a radius $r_\sigma > 0$ such that $B_\sigma \subset \mathcal{B}_K$, where we denote $B_\sigma = B(\mathbf{x}_\sigma, r_\sigma)$, as well as:

$$ \theta_{\mathcal{I}} = \max \left( \sup_{K \in \mathcal{T}} \sup_{\sigma \in \mathcal{F}_K} \frac{r_\sigma}{r_K}, \left( \inf_{K \in \mathcal{T}} \inf_{\sigma \in \mathcal{F}_K} \frac{r_\sigma}{r_K} \right)^{-1} \right). $$

We define the operator $\mathcal{I}: H^1(\Omega) \to X_N$ by setting $\mathcal{I}(v) = (\mathcal{I}_\sigma(v))_{\sigma \in \mathcal{F}}$ where:

$$ \mathcal{I}_{\sigma}(v) = \frac{1}{|B_{\sigma}|} \int_{B_{\sigma}} Ev \quad \text{for any } \sigma \in \mathcal{F}, \quad (41) $$

and we of course denote $\mathcal{I}_K(v) = (\mathcal{I}_\sigma(v))_{\sigma \in \mathcal{F}_K}$. We also introduce the operator $\mathcal{I}: H_0^1(\Omega) \to X_N$

$$ \mathcal{I}_{\sigma}^{0}(v) = \begin{cases} \frac{1}{|B_{\sigma}|} \int_{B_{\sigma}} Ev & \text{for any } \sigma \in \mathcal{F}_{\text{int}} \\ 0 & \text{for any } \sigma \in \mathcal{F}_{\text{ext}}. \end{cases} \quad (42) $$

## 5.2 Local consistency for the network element interpolation

**Proposition 5.1** (Local approximation results for network element interpolation). Assume that $\Omega$ is Lipschitz and satisfies the cone condition with angle $\tau$ and radius $r$. Let $\mathcal{N}$ be an admissible network and $\mathcal{G}$ an associated admissible geometry. Then, we have:

$$ \text{For any } v \in H^1(\Omega) \qquad \|v - \mathcal{M}_K(\mathcal{I}_K(v))\|_{L^2(\mathcal{B}_K \cap \Omega)} \le C h_K |Ev|_{H^1(\mathcal{B}_K)}. \quad (43) $$

$$ \text{For any } v \in H^2(\Omega): \qquad \|\nabla v - \nabla_K(\mathcal{I}_K(v))\|_{L^2(\mathcal{B}_K \cap \Omega)^d} \le C(h_K + h_K^p) \|Ev\|_{H^2(\mathcal{B}_K)}, \quad (44) $$

and

$$ \|v - \Pi_K(\mathcal{I}_K(v))\|_{L^2(\mathcal{B}_K \cap \Omega)} \le C(h_K^2 + h_K^{p+1}) \|Ev\|_{H^2(\mathcal{B}_K)}. \quad (45) $$

For any $\Phi \in H^2(\Omega)^d$:

$$ \|\mathrm{div}(\Phi) - \mathcal{D}\mathcal{I}\mathcal{V}_K(\mathcal{I}_K(\Phi))\|_{L^2(\mathcal{B}_K \cap \Omega)} \le C(h_K + h_K^p) \|\Phi\|_{H^2(\mathcal{B}_K)}, \quad (46) $$

where the constants $C > 0$ in the above result can vary from line to line but only depend on the quality parameters $\theta_A$, $\theta_T$, $\theta_F$, $\theta_M$, $\theta_\Pi$, $\eta_N$, $\eta_\psi$, $M_\psi$, $\theta_\mathcal{I}$, and not on $h$. The same results hold replacing $\mathcal{I}$ by $\mathcal{I}^0$ under the additional hypothesis that the functions $\varphi$ (resp. $\Phi$) belong to $H_0^1(\Omega)$ (resp. $H_0^1(\Omega)^d$).
Proof. First remark that by density of $C^\infty(\bar{\Omega})$ in $H^1(\Omega)$ and $H^2(\Omega)$, it suffices to establish the results for $\varphi \in C^\infty(\bar{\Omega})$. For any $\varphi \in C^\infty(\bar{\Omega})$, let us denote $\tilde{\varphi} = E\varphi$ to avoid repeating the notation E everywhere. Using Taylor's expansion formula, we have for any $\sigma \in \mathcal{F}_K$ and any $(\mathbf{x}, \mathbf{y}) \in \mathcal{B}_K^2$:

$$ \tilde{\varphi}(\mathbf{x}) = \tilde{\varphi}(\mathbf{y}) + \sum_{|\alpha|=1} (\mathbf{x} - \mathbf{y})^\alpha \int_0^1 \partial^\alpha \tilde{\varphi}(\mathbf{x} + t(\mathbf{y} - \mathbf{x}))dt, $$

which immediately gives, as $\sum_{\sigma \in \mathcal{F}_K} \gamma_{K,\sigma} = 1$:

$$ \tilde{\varphi}(\mathbf{x}) = \sum_{\sigma \in \mathcal{F}_K} \frac{\gamma_{K,\sigma}}{|B_\sigma|} \int_{B_\sigma} \tilde{\varphi}(\mathbf{x})dy = \mathcal{M}_K(\mathcal{I}_K(\varphi)) + \sum_{\sigma \in \mathcal{F}_K} \frac{\gamma_{K,\sigma}}{|B_\sigma|} \int_{B_\sigma} \sum_{|\alpha|=1} (\mathbf{x} - \mathbf{y})^\alpha \left( \int_0^1 \partial^\alpha \tilde{\varphi}(\mathbf{x} + t(\mathbf{y} - \mathbf{x}))dt \right) dy. $$

Next, denoting $R_{\mathcal{M}}\varphi(\mathbf{x}) = \tilde{\varphi}(\mathbf{x}) - \mathcal{M}_K(\mathcal{I}_K(\varphi))$, following [4] we define the change of variable $(\mathbf{y}, t) \to (\mathbf{z}, t)$ for which $\mathbf{z} = \mathbf{x} + t(\mathbf{y} - \mathbf{x})$ and $dtdy = t^{-d}dtdz$. The domain of integration for $(\mathbf{y}, t)$ is $B_\sigma \times ]0, 1[$. Its image by the above change of variable is:

$$ D_{\mathbf{x},\sigma} = \left\{ (\mathbf{z}, t) \mid t \in ]0, 1[, \left| \frac{1}{t}(\mathbf{z} - \mathbf{x}) + \mathbf{x} - \mathbf{x}_\sigma \right| < r_\sigma \right\}. $$

Notice that $(\mathbf{x} - \mathbf{y})^\alpha = t^{-m}(\mathbf{x} - \mathbf{z})^\alpha$ if $|\alpha| = m$ and that if $(\mathbf{z}, t) \in D_{\mathbf{x},\sigma}$ then:

$$ \left| \frac{1}{t} |\mathbf{z} - \mathbf{x}| - |\mathbf{x}_\sigma - \mathbf{x}| \right| \le \left| \frac{1}{t} (\mathbf{z} - \mathbf{x}) + \mathbf{x} - \mathbf{x}_\sigma \right| < r_\sigma \quad \text{and} \quad s(\mathbf{x}, \mathbf{z}) = \frac{|\mathbf{z} - \mathbf{x}|}{|\mathbf{x} - \mathbf{x}_\sigma| + r_\sigma} < t. \quad (47) $$

Using the above change of variable, we obtain:

$$ R_{\mathcal{M}}\varphi(\mathbf{x}) = \sum_{\sigma \in \mathcal{F}_K} \frac{\gamma_{K,\sigma}}{|B_\sigma|} \sum_{|\alpha|=1} \int_{D_{\mathbf{x},\sigma}} (\mathbf{x} - \mathbf{z})^\alpha \partial^\alpha \tilde{\varphi}(\mathbf{z}) t^{-d-1} dtdz. $$

The projection of $D_{\mathbf{x},\sigma}$ on the $\mathbf{z}$-space being the convex hull of $\{\mathbf{x}\} \cup B_\sigma$, denoted $C_{\mathbf{x},\sigma}$, applying Fubini-Tonelli's theorem we get:

$$ R_{\mathcal{M}}\varphi(\mathbf{x}) = \sum_{\sigma \in \mathcal{F}_K} \frac{\gamma_{K,\sigma}}{|B_\sigma|} \sum_{|\alpha|=1} \int_{C_{\mathbf{x},\sigma}} (\mathbf{x} - \mathbf{z})^\alpha \partial^\alpha \tilde{\varphi}(\mathbf{z}) \int_0^1 \chi_{D_{\mathbf{x},\sigma}}(\mathbf{z}, t) t^{-d-1} dtdz. $$

Using (47) we get:

$$ \left| \int_0^1 \chi_{D_{\mathbf{x},\sigma}}(\mathbf{z}, t) t^{-d-1} dtdz \right| \le \int_{s(\mathbf{x},\mathbf{z})}^1 t^{-d-1} dt \le \frac{1}{d} (s(\mathbf{x},\mathbf{z})^{-d} - 1) \le \frac{1}{d} (|\mathbf{x} - \mathbf{x}_\sigma| + r_\sigma)^d |\mathbf{x} - \mathbf{z}|^{-d}. $$

Thus, injecting this in the above expression for $R_{\mathcal{M}}\varphi(\mathbf{x})$, we get as $|\alpha| = 1$ and $|\mathbf{x} - \mathbf{x}_\sigma| \le \rho_K$:

$$ |R_{\mathcal{M}}\varphi(\mathbf{x})| \le \sum_{\sigma \in \mathcal{F}_K} \sum_{|\alpha|=1} \frac{\theta_{\mathcal{M}}}{dS_1^d r_\sigma^d} (\rho_K + r_\sigma)^d \int_{C_{\mathbf{x},\sigma}} |\partial^\alpha \tilde{\varphi}(\mathbf{z})| |\mathbf{x} - \mathbf{z}|^{1-d} dz. $$

Recall that $\theta_{\mathcal{I}}^{-1}r_K \le r_\sigma \le \theta_{\mathcal{I}}r_K$ and $\rho_K \le \kappa_\psi r_K$ and thus finally as $C_{\mathbf{x},\sigma} \subset \mathcal{B}_K$ by construction:

$$ |R_{\mathcal{M}}\varphi(\mathbf{x})| \le \frac{(d+1)\theta_{\mathcal{I}}^d \theta_{\mathcal{M}} \theta_{\mathcal{F}}}{S_1^d} (\kappa_\psi + \theta_{\mathcal{I}})^d \int_{\mathcal{B}_K} |\partial^\alpha \tilde{\varphi}(\mathbf{z})| |\mathbf{x} - \mathbf{z}|^{1-d} dz. $$
Consequently, applying (39) we obtain the first result (43):

$$ \int_{\mathcal{B}_K \cap \Omega} |\varphi - \mathcal{M}_K(\mathcal{I}_K(\varphi))|^2 \le C_{1,d}^2 \frac{((d+1)\theta_{\mathcal{I}}^d \theta_{\mathcal{M}} \theta_{\mathcal{F}\mathcal{K}\psi})^2}{4S_1^{d2}} (\kappa_{\psi} + \theta_{\mathcal{I}})^{2d} h_K^2 \|\nabla \tilde{\varphi}\|_{L^2(\mathcal{B}_K)}^2 $$

To establish (44) we proceed the same way however the geometric approximation errors require a specific treatment. Using again Taylor's expansion, we see that for any $\sigma \in \mathcal{F}_K$ and $(\mathbf{x}, \mathbf{y}) \in \mathcal{B}_K^2$:

$$ \tilde{\varphi}(\mathbf{x}) = \tilde{\varphi}(\mathbf{y}) + \nabla\tilde{\varphi}(\mathbf{y}) \cdot (\mathbf{x} - \mathbf{y}) + 2 \sum_{|\alpha|=2} \frac{(\mathbf{x} - \mathbf{y})^\alpha}{\alpha!} \int_0^1 t \partial^\alpha \tilde{\varphi}(\mathbf{x} + t(\mathbf{y} - \mathbf{x})) dt. $$

Then, using (4):

$$ \begin{aligned} \varepsilon_K^0 \tilde{\varphi}(\mathbf{x}) = & \frac{1}{m_K} \sum_{\sigma \in \mathcal{F}_K} \tilde{\varphi}(\mathbf{x}) \eta_{K,\sigma} = \frac{1}{m_K} \sum_{\sigma \in \mathcal{F}_K} \frac{1}{|B_\sigma|} \int_{B_\sigma} \tilde{\varphi}(\mathbf{y}) \eta_{K,\sigma} dy \\ & + \frac{1}{m_K} \sum_{\sigma \in \mathcal{F}_K} \frac{1}{|B_\sigma|} \int_{B_\sigma} \nabla \tilde{\varphi}(\mathbf{y}) \cdot (\mathbf{x} - \mathbf{y}) \eta_{K,\sigma} dy \\ & + \frac{2}{m_K} \sum_{\sigma \in \mathcal{F}_K} \frac{1}{|B_\sigma|} \int_{B_\sigma} \sum_{|\alpha|=2} \frac{(\mathbf{x} - \mathbf{y})^\alpha}{\alpha!} \int_0^1 t \partial^\alpha \tilde{\varphi}(\mathbf{x} + t(\mathbf{y} - \mathbf{x})) \eta_{K,\sigma} dt dy = \varepsilon_1 + \varepsilon_2 + \varepsilon_3, \end{aligned} $$

each $\mathcal{E}_i$ corresponding to one line in the above expression. By definition, the first term of the above expression is exactly:

$$ \mathcal{E}_1 = \frac{1}{m_K} \sum_{\sigma \in \mathcal{F}_K} \frac{1}{|B_\sigma|} \int_{B_\sigma} \tilde{\varphi}(\mathbf{y}) \eta_{K,\sigma} dy = \frac{1}{m_K} \sum_{\sigma \in \mathcal{F}_K} \mathcal{I}_\sigma(\tilde{\varphi}) \eta_{K,\sigma} dy = \nabla \Pi_K(\mathcal{I}_K(\varphi)). $$

From Taylor's expansion, we get:

$$ \nabla \tilde{\varphi}(\mathbf{y}) = \nabla \tilde{\varphi}(\mathbf{x}) + \sum_{|\alpha|=1} (\mathbf{y} - \mathbf{x})^\alpha \int_0^1 \partial^\alpha \nabla \tilde{\varphi}(\mathbf{x} + t(\mathbf{y} - \mathbf{x})) dt. $$

Using this expansion and the fact that $\frac{1}{|B_\sigma|} \int_{B_\sigma} \mathbf{y} dy = \mathbf{x}_\sigma$, the second term of the above expression rewrites:

$$ \begin{aligned} \mathcal{E}_2 = & \frac{1}{m_K} \sum_{\sigma \in \mathcal{F}_K} \frac{1}{|B_\sigma|} \int_{B_\sigma} \nabla \tilde{\varphi}(\mathbf{y}) \cdot (\mathbf{x} - \mathbf{y}) \eta_{K,\sigma} dy = \frac{1}{m_K} \sum_{\sigma \in \mathcal{F}_K} \nabla \tilde{\varphi}(\mathbf{x}) \cdot (\mathbf{x} - \mathbf{x}_\sigma) \eta_{K,\sigma} \\ & + \frac{1}{m_K} \sum_{\sigma \in \mathcal{F}_K} \frac{1}{|B_\sigma|} \int_{B_\sigma} \sum_{|\alpha|=1} (\mathbf{y} - \mathbf{x})^\alpha \int_0^1 \partial^\alpha \nabla \tilde{\varphi}(\mathbf{x} + t(\mathbf{y} - \mathbf{x})) \cdot (\mathbf{x} - \mathbf{y}) \eta_{K,\sigma} dt dy, \end{aligned} $$

and thus, using the first order approximate consistency properties (4) and (5):

$$ \begin{aligned} \mathcal{E}_2 = & -\nabla \tilde{\varphi}(\mathbf{x}) - \nabla \tilde{\varphi}(\mathbf{x}) \cdot (\mathbf{x}_K - \mathbf{x}) \varepsilon_K^0 - \sum_{i=1}^d \sum_{j=1}^d \varepsilon_K^{1,ij} \partial_{x_j} \tilde{\varphi}(\mathbf{x}) e_i \\ & + \frac{1}{m_K} \sum_{\sigma \in \mathcal{F}_K} \frac{1}{|B_\sigma|} \int_{B_\sigma} \sum_{|\alpha|=1} (\mathbf{y} - \mathbf{x})^\alpha \int_0^1 \partial^\alpha \nabla \tilde{\varphi}(\mathbf{x} + t(\mathbf{y} - \mathbf{x})) \cdot (\mathbf{x} - \mathbf{y}) \eta_{K,\sigma} dt dy. \end{aligned} $$

Gathering the previous results, we get that:

$$ R_{\Pi}(\varphi)(\mathbf{x}) = \nabla \tilde{\varphi}(\mathbf{x}) - \nabla_K(\mathcal{I}_K(\varphi)) = -(\tilde{\varphi}(\mathbf{x}) + \nabla \tilde{\varphi}(\mathbf{x}) \cdot (\mathbf{x}_K - \mathbf{x})) \varepsilon_K^0 - \sum_{i=1}^d \sum_{j=1}^d \varepsilon_K^{1,ij} \partial_{x_j} \tilde{\varphi}(\mathbf{x}) e_i $$
$$
+ \frac{1}{m_K} \sum_{\sigma \in \mathcal{F}_K} \frac{1}{|B_\sigma|} \int_{B_\sigma} \sum_{|\alpha|=1} (\mathbf{y} - \mathbf{x})^\alpha \int_0^1 \partial^\alpha \nabla \tilde{\varphi}(\mathbf{x} + t(\mathbf{y} - \mathbf{x})) \cdot (\mathbf{x} - \mathbf{y}) \eta_{K,\sigma} dt dy \\
+ \frac{2}{m_K} \sum_{\sigma \in \mathcal{F}_K} \frac{1}{|B_\sigma|} \int_{B_\sigma} \sum_{|\alpha|=2} \frac{(\mathbf{x} - \mathbf{y})^\alpha}{\alpha!} \int_0^1 t \partial^\alpha \tilde{\varphi}(\mathbf{x} + t(\mathbf{y} - \mathbf{x})) \eta_{K,\sigma} dt dy \\
= R_{\Pi,1}(\varphi)(\mathbf{x}) + R_{\Pi,2}(\varphi)(\mathbf{x}) + R_{\Pi,3}(\varphi)(\mathbf{x}),
$$

each $R_{\Pi,i}(\varphi)(\mathbf{x})$ corresponding to one line in the previous expression. Immediately, we see that:

$$
\begin{align*}
\| R_{\Pi,1}(\varphi)(\mathbf{x}) \|_{L^2(B_K \cap \Omega)}^2 &= \int_{B_K \cap \Omega} \left| \left( \tilde{\varphi}(\mathbf{x}) + \nabla \tilde{\varphi}(\mathbf{x}) \cdot (\mathbf{x}_K - \mathbf{x}) \right) \varepsilon_K^0 + \sum_{i=1}^d \sum_{j=1}^d \varepsilon_K^{1,ij} \partial_{x_j} \tilde{\varphi}(\mathbf{x}) e_i \right|^2 dx \\
&\le 2 \theta_A^2 h^{2p} \left( d \| \tilde{\varphi} \|_{L^2(B_K \cap \Omega)}^2 + (hd^{1/2} + d^2)^2 \| \nabla \tilde{\varphi} \|_{L^2(B_K \cap \Omega)}^2 \right).
\end{align*}
$$

Using the same change of variable as above and (39) to estimate $R_{\Pi,2}(\varphi)$ and $R_{\Pi,3}(\varphi)$, this finally leads to the second estimate (44). The proof of estimate (46) follows the same lines, while estimate (45) can be established by refining our expansion of $R_M(\varphi)$ through Taylor's expansion and proceeding as above to estimate the residual terms. As for the results involving $\mathcal{I}^0$, it suffices to notice that $0 = \varphi(\mathbf{x}_\sigma)$ for $\sigma \in \mathcal{F}_{\text{ext}}$ and $\varphi \in C_c^\infty(\Omega)$ and then proceed as above. $\square$

5.3 Global consistency for network element interpolation

**Proposition 5.2** (Global approximation results for network element interpolation). *Assume that Ω is Lipschitz and satisfies the cone condition with angle τ and radius r. Let N be an admissible network and G an associated admissible geometry. Then, we have:*

$$
\text{For any } v \in H^1(\Omega) \qquad \|v - \Pi_{\mathcal{T}}(\mathcal{I}(v))\|_{L^2(\Omega)} \le C_{E,1} Ch \|v\|_{H^1(\Omega)}, \qquad (48)
$$

$$
\text{and} \qquad \|\mathcal{I}(v)\|_{X} \le C_{E,1}C \|v\|_{H^1(\Omega)}. \qquad (49)
$$

$$
\text{For any } v \in H^2(\Omega) \qquad \|\nabla v - \nabla_{\mathcal{T}}(\mathcal{I}(v))\|_{L^2(\Omega)} \le C_{E,2} C(h+h^p) \|v\|_{H^2(\Omega)}, \qquad (50)
$$

$$
\text{and} \qquad \|v - \Pi_{\mathcal{N}}(\mathcal{I}(v))\|_{L^2(\Omega)} \le C_{E,2} C(h^2 + h^{p+1}) \|v\|_{H^2(\Omega)}. \qquad (51)
$$

$$
\text{For any } \Phi \in H^2(\Omega)^d: \qquad \| \mathrm{div}(\Phi) - \mathcal{D} \mathcal{I} \mathcal{V}_{\mathcal{T}}(\mathcal{I}(\Phi)) \|_{L^2(\Omega)} \le C_{E,2} C(h+h^p) \|\Phi\|_{H^2(\Omega)^d}, \qquad (52)
$$

where

$$
\mathcal{D}\mathcal{I}\mathcal{V}_{\mathcal{T}}(\mathcal{I}(\Phi)) = \sum_{K \in \mathcal{T}} \psi_K \mathcal{D}\mathcal{I}\mathcal{V}_K(\mathcal{I}_K(\Phi)),
$$

and where the constants $C > 0$ in the above result can vary from line to line but only depend on the quality parameters $\theta_A$, $\theta_T$, $\theta_F$, $\theta_M$, $\theta_\Pi$, $\eta_\mathcal{N}$, $\eta_\psi$, $M_\psi$, $\theta_\mathcal{I}$, and not on $h$. The same results hold replacing $\mathcal{I}$ by $\mathcal{I}^0$ under the additional hypothesis that the functions $\varphi$ (resp. $\Phi$) belong to $H_0^1(\Omega)$ (resp. $H_0^1(\Omega)^d$).

*Proof.* Using the notations of proposition 5.1 and proceeding again by density, let $\varphi \in C^\infty(\bar{\Omega})$. For (48), we proceed along the lines of lemma 4.1:

$$
\| \Pi_{\mathcal{T}}(\mathcal{I}(\varphi)) - \varphi \|_{L^2(\Omega)}^2 \le \sum_{K \in \mathcal{T}} \| \Pi_{\mathcal{T}}(\mathcal{I}(\varphi)) - \varphi \|_{L^2(B_K \cap \Omega)}^2 \le 2 \sum_{K \in \mathcal{T}} \| \varphi - M_K(\mathcal{I}_K(\varphi)) \|_{L^2(B_K \cap \Omega)}^2
$$
$$+2 \sum_{K \in \mathcal{T}} \left\| \sum_{L \in \mathcal{T}} \psi_L(\mathcal{M}_K(\mathcal{I}_K(\varphi)) - \mathcal{M}_L(\mathcal{I}_L(\varphi))) \right\|_{L^2(\mathcal{B}_K \cap \Omega)}^2 = 2 \sum_{K \in \mathcal{T}} (I_{K,1} + I_{K,2}),$$

with obvious notations. For the first term, using estimate (43) we obviously have for some $C > 0$ independent on $h$:

$$\sum_{K \in \mathcal{T}} I_{K,1} = \sum_{K \in \mathcal{T}} \|\varphi - \mathcal{M}_K(\mathcal{I}_K(\varphi))\|_{L^2(\mathcal{B}_K \cap \Omega)}^2 \le \sum_{K \in \mathcal{T}} C^2 h_K^2 \|\nabla \tilde{\varphi}\|_{L^2(\mathcal{B}_K)}^2 \le \eta_\psi C_{E,1}^2 C^2 h^2 \|\nabla \varphi\|_{L^2(\Omega)}^2.$$

Proceeding as in the proof of lemma 4.1 for the second term immediately leads to:

$$\left\| \sum_{L \in \mathcal{T}} \psi_L(\mathcal{M}_K(\mathcal{I}_K(\varphi)) - \mathcal{M}_L(\mathcal{I}_L(\varphi))) \right\|_{L^2(\mathcal{B}_K \cap \Omega)}^2 \le \eta_\psi M_\psi^2 \int_{\Omega} \sum_{L \in \mathcal{T}} |\mathcal{M}_K(\mathcal{I}_K(\varphi)) - \mathcal{M}_L(\mathcal{I}_L(\varphi))|^2 \chi_{\mathcal{B}_L} \chi_{\mathcal{B}_K}.$$

Then, notice that $|\mathcal{M}_K(\mathcal{I}_K(\varphi)) - \mathcal{M}_L(\mathcal{I}_L(\varphi))| \le |\mathcal{M}_K(\mathcal{I}_K(\varphi)) - \varphi| + |\varphi - \mathcal{M}_L(\mathcal{I}_L(\varphi))|$ which leads to, for some $C > 0$ independent on $h$ coming from estimate (43):

$$\begin{aligned} I_{K,2} &\le 2\eta_\psi M_\psi^2 \int_{\mathcal{B}_K \cap \Omega} |\mathcal{M}_K(\mathcal{I}_K(\varphi)) - \varphi|^2 \left(\sum_{L \in \mathcal{T}} \chi_L\right) + 2\eta_\psi M_\psi^2 \sum_{L \in \mathcal{T}} \sum_{\mathcal{B}_K \cap \mathcal{B}_L \ne \emptyset} \int_{\mathcal{B}_L \cap \Omega} |\varphi - \mathcal{M}_L(\mathcal{I}_L(\varphi))|^2 \chi_{\mathcal{B}_K} \\ &\le 2\eta_\psi^2 M_\psi^2 C h_K^2 \|\nabla \tilde{\varphi}\|_{L^2(\mathcal{B}_K)}^2 + 2\eta_\psi M_\psi^2 C \sum_{L \in \mathcal{T}} h_L^2 \int_{\mathcal{B}_L \cap \Omega} |\nabla \tilde{\varphi}|^2 \chi_{\mathcal{B}_K}. \end{aligned}$$

Summing over $K \in \mathcal{T}$, we get using Fubini-Tonelli's theorem:

$$\sum_{K \in \mathcal{T}} I_{K,2} \le 2\eta_\psi^3 M_\psi^2 C_{E,1}^2 C^2 h^2 \|\nabla \varphi\|_{L^2(\Omega)}^2 + 2\eta_\psi M_\psi^2 C \sum_{L \in \mathcal{T}} h_L^2 \int_{\mathcal{B}_L \cap \Omega} |\nabla \tilde{\varphi}|^2 \left(\int_{\Omega} \sum_{K \in \mathcal{T}} \chi_K\right),$$

and thus

$$\sum_{K \in \mathcal{T}} I_{K,2} \le 4\eta_\psi^3 M_\psi^2 C_{E,1}^2 C^2 h^2 \|\nabla \varphi\|_{L^2(\Omega)}^2.$$

Estimates (50) and (51) can be obtained proceeding the same way, following the lines of the proof of lemma 4.1 and using the local estimates established in the first part of the present proof. Estimate (49) can be obtained proceeding as in lemma 4.4 and using the above Taylor's expansions. $\square$

**Lemma 5.3** (Stabilization consistency for the network element interpolant). Let $(\mathcal{N}, \mathcal{G})$ be an admissible discretization network and an associated admissible network geometry. For any $v \in H^2(\Omega)$:

$$s_h(\mathcal{I}(v), \mathcal{I}(v)) \le S^* C_{E,2} C(h + h^p)^2 \|v\|_{H^2(\Omega)}^2.$$

where the constant $C > 0$ only depend on the quality parameters $\theta_A, \theta_T, \theta_F, \theta_M, \theta_\Pi, \eta_\psi, \eta_\psi, M_\psi, \theta_\Xi$, and not on $h$. The same result holds replacing $\mathcal{I}$ by $\mathcal{I}^0$ under the additional hypothesis that the functions $v$ belongs to $H_0^1(\Omega)$.

*Proof.* Again, by density of $C^\infty(\bar{\Omega})$ in $H^2(\Omega)$, it suffices to establish the result for $\varphi \in C^\infty(\bar{\Omega})$. We have:

$$s_h(\mathcal{I}(\varphi), \mathcal{I}(\varphi)) \le \theta_T S^* \sum_{K \in \mathcal{T}} \int_{B_K \cap \Omega} h_K^{-2} \sum_{\sigma \in \mathcal{F}_K} |\mathcal{I}_\sigma(\varphi) - \Pi_K(\mathcal{I}_K(\varphi))|(x_\sigma)|^2.$$

As $\Pi_K(\mathcal{I}_K(\varphi))$ is a first order polynomial, we have for any $\sigma \in \mathcal{F}_K$:

$$\frac{1}{|B_\sigma|} \int_{B_\sigma} \Pi_K(\mathcal{I}_K(\varphi))(x) = \Pi_K(\mathcal{I}_K(\varphi)(x_\sigma)).$$
Then using Cauchy-Schwarz inequality, for any $\sigma \in \mathcal{F}_K$:

$$|\mathcal{I}_{\sigma}(\varphi) - \Pi_{K}(\mathcal{I}_{K}(\varphi))(\mathbf{x}_{\sigma})|^2 \leq \frac{1}{|B_{\sigma}|} \int_{B_{\sigma}} |\tilde{\varphi} - \Pi_{K}(\mathcal{I}_{K}(\varphi)))|^2.$$

From the proof of proposition 5.1, using the notations defined there we know that Taylor's expansion gives for any $\mathbf{x} \in \mathcal{B}_K$:

$$\begin{aligned} \tilde{\varphi}(\mathbf{x}) - \Pi_K(\mathcal{I}_K(\varphi))(\mathbf{x}) ={}& R_{\Pi}(\varphi)(\mathbf{x}) \cdot (\mathbf{x} - \bar{\mathbf{x}}_K) \\ & + 2 \sum_{\sigma \in \mathcal{F}_K} \frac{\gamma_{K,\sigma}}{|B_{\sigma}|} \int_{B_{\sigma}} \sum_{|\alpha|=2} \frac{(\mathbf{x} - \mathbf{y})^{\alpha}}{\alpha!} \int_{0}^{1} t \partial^{\alpha} \tilde{\varphi}(\mathbf{x} + t(\mathbf{y} - \mathbf{x})) dtdy \\ & + \sum_{\sigma \in \mathcal{F}_K} \frac{\gamma_{K,\sigma}}{|B_{\sigma}|} \int_{B_{\sigma}} \sum_{|\alpha|=1} (\mathbf{y} - \mathbf{x})^{\alpha} \int_{0}^{1} \partial^{\alpha} \nabla \tilde{\varphi}(\mathbf{x} + t(\mathbf{y} - \mathbf{x})) \cdot (\mathbf{x} - \mathbf{y}) dtdy, \end{aligned}$$

and the results follows using the same techniques as in the proof of lemma (45). $\square$

## 5.4 Error estimates

We are now in position to establish error estimates for regular problems:

**Proposition 5.4.** Let $(\mathcal{N}, \mathcal{G})$ be an admissible discretization network and an associated admissible network geometry, and let $\mathbf{U}$ be the solution of the associated problem (16). Assume that there exists $C_f > 0$ such that:

$$\left( \sum_{K \in \mathcal{T}} \int_{B_K \cap \Omega} |f_K - f|^2 \right)^{1/2} \leq C_f h$$

and assume that the solution $u$ of (11) satisfies $u \in H^3(\Omega)$. Then, if $(\psi_K)_{K \in \mathcal{T}}$ is a quadrature family, we have the following error estimates:

$$\|\mathbf{U} - \mathcal{I}^0(u)\|_X \leq C(h + h^p), \qquad (53)$$

and

$$\|u - \Pi_{\mathcal{T}}(\mathbf{U})\|_{L^2(\Omega)} \leq C(h + h^p) \quad \text{and} \quad \|\nabla u - \nabla_{\mathcal{T}}(\mathbf{U})\|_{L^2(\Omega)} \leq C(h + h^p), \qquad (54)$$

and for any $K \in \mathcal{T}$

$$\|u - \mathcal{M}_K(\mathbf{U})\|_{L^2(B_K \cap \Omega)} \leq C(h + h^p) \quad \text{and} \quad \|\nabla u - \nabla_K(\mathbf{U})\|_{L^2(B_K \cap \Omega)^d} \leq C(h + h^p), \qquad (55)$$

where the constant $C > 0$ depends on $u$, $d$, $\theta_{\mathcal{F}}$, $\theta_{\mathcal{T}}$, $\theta_{\Pi}$, $\theta_{\mathcal{A}}$, $\theta_{\mathcal{M}}$, $M_{\psi}$, $\kappa_{\psi}$, $\eta_{\psi}$, $\theta_{\mathcal{I}}$, and $\Omega$ but not on $h$.

*Proof.* In the following, $C > 0$ denotes a constant that can depend on $u$, $d$, $\theta_{\mathcal{F}}$, $\theta_{\mathcal{T}}$, $\theta_{\Pi}$, $\theta_{\mathcal{A}}$, $\theta_{\mathcal{M}}$, $M_{\psi}$, $\kappa_{\psi}$, $\eta_{\psi}$, $\theta_{\mathcal{I}}$ and $\Omega$ whose value can change from one line to another. First, as we have already noticed in the proof of proposition 4.5, by (19)-(20) and (23), it is clear that:

$$\|\mathbf{U}_h\|_X \leq C \left( \sum_{K \in \mathcal{T}} m_K |f_K|^2 \right)^{\frac{1}{2}} \leq C(1+h) \|f\|_{L^2(\Omega)}.$$

As $\mathcal{I}^0(u) \in X_{\mathcal{N},0}$ by construction, we can use it in our discrete variational problem. Consequently for any $\mathbf{V} \in X_{\mathcal{N},0}$:

$$a_h(\mathcal{I}^0(u), \mathbf{V}) = \int_{\Omega} \nabla u \cdot \nabla_{\mathcal{T}}(\mathbf{V}) + \sum_{K \in \mathcal{T}} \int_{B_K \cap \Omega} \psi_K(\nabla \Pi_K(\mathcal{I}^0(u)) - \nabla u) \cdot \nabla \Pi_K(\mathbf{V}) + s^h(\mathcal{I}^0(u), \mathbf{V}).$$
Using Cauchy-Schwarz inequality and (19)-(20) and lemma 5.3, we get:

$$|s^h(\mathcal{I}^0(u), \mathbf{V})| \le |s^h(\mathcal{I}^0(u), \mathcal{I}^0(u))|^{1/2} |s^h(\mathbf{V}, \mathbf{V})|^{1/2} \le C(h + h^p) ||\mathbf{V}||_X |u|_{H^2(\Omega)}$$

while the same Cauchy-Schwarz inequality and (44) leads to:

$$\left| \sum_{K \in \mathcal{T}} \int_{\mathcal{B}_K \cap \Omega} \psi_K (\nabla \Pi_K(\mathcal{I}_K(u)) - \nabla u) \cdot \nabla \Pi_K(\mathbf{V}) \right| \le C(h + h^p) |u|_{H^2(\Omega)} ||\mathbf{V}||_X.$$

Then, using the fact that $U$ is solution of the discrete problem (16), we have that:

$$a_h(\mathbf{U}, \mathbf{V}) = l_h(\mathbf{V}) = \int_{\Omega} f \Pi_{\mathcal{T}}(\mathbf{V}) + \sum_{K \in \mathcal{T}} \int_{\mathcal{B}_K \cap \Omega} \psi_K (f_K - f) \mathcal{M}_K(\mathbf{V}).$$

Using Cauchy-Schwarz inequality, we can bound the second term by:

$$\left| \sum_{K \in \mathcal{T}} \int_{\mathcal{B}_K} \psi_K (f - f_K) \mathcal{M}_K(\mathbf{V}) \right| \le C ||\mathbf{V}||_0 \left( \sum_{K \in \mathcal{T}} \int_{\mathcal{B}_K \cap \Omega} |f_K - f|^2 \right)^{1/2}.$$

Combining the above results, we see that there exists $C > 0$ such that for all $V \in X_{\mathcal{N},0}$, we have:

$$|a_h(\mathbf{U} - \mathcal{I}^0(u), \mathbf{V})| \le \left| \int_{\Omega} \nabla u \cdot \nabla_{\mathcal{T}}(\mathbf{V}) - \int_{\Omega} f \Pi_{\mathcal{T}}(\mathbf{V}) \right| + C(h + h^p) ||\mathbf{V}||_X.$$

Consequently, it just remains to estimate:

$$\mathcal{R}(\mathbf{V}) = \int_{\Omega} \nabla u \cdot \nabla_{\mathcal{T}}(\mathbf{V}) - \int_{\Omega} f \Pi_{\mathcal{T}}(\mathbf{V}) = \int_{\Omega} \nabla u \cdot \nabla_{\mathcal{T}}(\mathbf{V}) + \int_{\Omega} \Delta u \Pi_{\mathcal{T}}(\mathbf{V}).$$

As $(\psi_K)_{K \in \mathcal{T}}$ is a quadrature family, we get:

$$\int_{\Omega} \Delta u \Pi_{\mathcal{T}}(\mathbf{V}) = \sum_{K \in \mathcal{T}} \int_{\Omega} \psi_K \mathrm{div}(\nabla u) \mathcal{M}_K(\mathbf{V}) \\ = \sum_{K \in \mathcal{T}} \sum_{\sigma \in \mathcal{F}_K} \eta_{K, \sigma} \cdot \mathcal{I}_{\sigma}(\nabla u) \mathcal{M}_K(\mathbf{V}) + \sum_{K \in \mathcal{T}} \int_{\Omega} \psi_K (\mathrm{div}(\nabla u) - \mathcal{D} \mathcal{I} \mathcal{V}_K(\mathcal{I}_K(\nabla u))) \mathcal{M}_K(\mathbf{V}).$$

As $v_\sigma = 0$ for any $\sigma \in \mathcal{F}_{\text{ext}}$ and as the geometry is approximately conservative, we get:

$$\int_{\Omega} \Delta u \Pi_{\mathcal{T}}(\mathbf{V}) = \sum_{K \in \mathcal{T}} \sum_{\sigma \in \mathcal{F}_K} \eta_{K, \sigma} \cdot \mathcal{I}_{\sigma}(\nabla u) (\mathcal{M}_K(\mathbf{V}) - v_{\sigma}) + \sum_{\sigma \in \mathcal{F}_{\text{int}}} \varepsilon_{\sigma} \cdot \mathcal{I}_{\sigma}(\nabla u) v_{\sigma} \\ + \sum_{K \in \mathcal{T}} \int_{\Omega} \psi_K (\mathrm{div}(\nabla u) - \mathcal{D} \mathcal{I} \mathcal{V}_K(\mathcal{I}_K(\nabla u))) \mathcal{M}_K(\mathbf{V}).$$

We denote:

$$G_K = \frac{1}{|B_K \cap \Omega|} \int_{B_K \cap \Omega} \nabla u.$$

then we get:

$$\int_{\Omega} \Delta u \Pi_{\mathcal{T}}(\mathbf{V}) = \sum_{K \in \mathcal{T}} \sum_{\sigma \in \mathcal{F}_K} \eta_{K, \sigma} \cdot (\mathcal{I}_{\sigma}(\nabla u) - G_K)(\mathcal{M}_K(\mathbf{V}) - v_{\sigma})$$
$$
\begin{align*}
& - \sum_{K \in \mathcal{T}} m_K G_K \cdot \nabla \Pi_K(\mathbf{V}) + \sum_{\sigma \in \mathcal{F}_{\text{int}}} \varepsilon_\sigma \cdot \mathcal{I}_\sigma(\nabla u) v_\sigma + \sum_{K \in \mathcal{T}} m_K \mathcal{M}_K(\mathbf{V}) G_K \cdot \varepsilon_K^0 \\
& \qquad + \sum_{K \in \mathcal{T}} \int_{\Omega} \psi_K(\text{div}(\nabla u) - \mathcal{D} \mathcal{I} \mathcal{V}_K(\mathbf{I}_K(\nabla u))) \mathcal{M}_K(\mathbf{V}) \\
&= \sum_{K \in \mathcal{T}} \sum_{\sigma \in \mathcal{F}_K} \eta_{K,\sigma} \cdot (\mathcal{I}_\sigma(\nabla u) - G_K)(\mathcal{M}_K(\mathbf{V}) - v_\sigma) + \int_{\Omega} \sum_{K \in \mathcal{T}} \psi_K(\nabla u - G_K) \cdot \nabla \Pi_K(\mathbf{V}) \\
& \qquad - \int_{\Omega} \nabla u \cdot \nabla_{\mathcal{T}}(\mathbf{V}) + \sum_{K \in \mathcal{T}} \int_{\Omega} \psi_K(\text{div}(\nabla u) - \mathcal{D} \mathcal{I} \mathcal{V}_K(\mathbf{I}_K(\nabla u))) \mathcal{M}_K(\mathbf{V}) \\
& \qquad + \sum_{\sigma \in \mathcal{F}_{\text{int}}} \varepsilon_\sigma \cdot \mathcal{I}_\sigma(\nabla u) v_\sigma + \sum_{K \in \mathcal{T}} m_K \mathcal{M}_K(\mathbf{V}) G_K \cdot \varepsilon_K^0,
\end{align*}
$$

and thus:

$$
\begin{align*}
\mathcal{R}(\mathbf{V}) = \sum_{K \in \mathcal{T}} \sum_{\sigma \in \mathcal{F}_K} \eta_{K,\sigma} \cdot (\mathcal{I}_\sigma(\nabla u) - G_K)(\mathcal{M}_K(\mathbf{V}) - v_\sigma) &+ \int_\Omega \sum_{K \in \mathcal{T}} \psi_K(\nabla u - G_K) \cdot \nabla \Pi_K(\mathbf{V}) \\
&+ \sum_{K \in \mathcal{T}} \int_\Omega \psi_K(\text{div}(\nabla u) - \mathcal{D} \mathcal{I} \mathcal{V}_K(\mathbf{I}_K(\nabla u))) \mathcal{M}_K(\mathbf{V}) + \sum_{\sigma \in \mathcal{F}_{\text{int}}} \varepsilon_\sigma \cdot \mathcal{I}_\sigma(\nabla u) v_\sigma + \sum_{K \in \mathcal{T}} m_K \mathcal{M}_K(\mathbf{V}) G_K \cdot \varepsilon_K^0.
\end{align*}
$$

We rewrite this last identity $\mathcal{R}(\mathbf{V}) = R_1 + R_2 + R_3 + R_4 + R_5$ with obvious notations. Proceeding as in the proof of proposition 5.2, it is clear as $\nabla u \in H^1(\Omega)$ that there exists $C > 0$ (applying Stein's extension theorem to $\nabla u$) such that:

$$
|\mathcal{I}_{\sigma}(\nabla u) - G_K| \le \frac{Ch}{|B_K \cap \Omega|} ||\nabla u||_{H^1(B_K \cap \Omega)} \quad \text{and} \quad \left|\int_{B_K \cap \Omega} |\nabla u - G_K|^2\right|^{1/2} \le Ch ||\nabla u||_{H^1(B_K \cap \Omega)}.
$$

Consequently, using Cauchy-Schwarz inequality we get that:

$$
\begin{align*}
|R_1| &\le \sum_{K \in \mathcal{T}} \sum_{\sigma \in \mathcal{F}_K} \frac{h_K |\eta_{K,\sigma}|}{m_K} m_K h_K^{-1} |\mathcal{I}_\sigma(\nabla u) - G_K| |\mathcal{M}_K(\mathbf{V}) - v_\sigma| \\
&\le \theta_\Pi \left( \sum_{K \in \mathcal{T}} m_K |\mathcal{I}_\sigma(\nabla u) - G_K|^2 \right)^{1/2} \left( \sum_{K \in \mathcal{T}} m_K h_K^{-2} |\mathcal{M}_K(\mathbf{V}) - v_\sigma|^2 \right)^{1/2} \\
&\le Ch \eta^{1/2} \theta_\Pi^{1/2} \theta_\Pi ||\nabla u||_{H^1(\Omega)} ||\mathbf{V}||_X,
\end{align*}
$$

and

$$
|R_2| \le \left( \sum_{K \in \mathcal{T}} \int_{B_K \cap \Omega} |\psi_K||\nabla u - G_K|^2 \right)^{1/2} \left( \sum_{K \in \mathcal{T}} \left( \frac{1}{m_K} \int_{\Omega} |\psi_K| \right) m_K ||\nabla \Pi_K(\mathbf{V})||^2 \right)^{1/2} \\
\le \left( \frac{\eta_N M_\psi^2 \theta_\Pi S_1^d \kappa_\psi^d \delta^d}{|C(0, \tau, 1)|} \right)^{1/2} Ch ||\nabla u||_{H^1(\Omega)} ||\mathbf{V}||_X,
$$

as well as:

$$
|R_3| \le C(h + h^p) |\nabla u|_{H^2(\Omega)^d} ||\mathbf{V}||_X.
$$

Next, using the convention $\varepsilon_\sigma = 0$ for $\sigma \in \mathcal{F}_{\text{ext}}$, we have:

$$
|R_4| = \left| \sum_{\sigma \in \mathcal{F}} \sum_{K \in \mathcal{T}_\sigma} \frac{\varepsilon_\sigma \cdot \mathcal{I}_\sigma(\nabla u) v_\sigma}{\mathrm{card}(\mathcal{T}_\sigma)} \right|
$$
$$ \leq \theta_A h^p \left( \sum_{K \in \mathcal{T}} \sum_{\sigma \in \mathcal{F}_K} \frac{m_K}{\mathrm{card}(\mathcal{T}_\sigma)} |\mathcal{I}_\sigma(\nabla u)| |v_\sigma - \mathcal{M}_K(\mathbf{V})| + \sum_{K \in \mathcal{T}} m_K |\mathcal{M}_K(\mathbf{V})| \left( \sum_{\sigma \in \mathcal{F}_K} \frac{1}{\mathrm{card}(\mathcal{T}_\sigma)} |\mathcal{I}_\sigma(\nabla u)| \right) \right), $$

and using Cauchy-Schwarz inequality, this leads to:

$$ |R_4| \leq \theta_A h^p (h |\mathbf{V}|_X + |\mathbf{V}|_0) \left( \sum_{K \in \mathcal{T}} \sum_{\sigma \in \mathcal{F}_K} m_K |\mathcal{I}_\sigma(\nabla u)|^2 \right)^{1/2}. $$

Then, notice that by definition of $\mathcal{B}_K$, we have:

$$ |\mathcal{I}_\sigma(\nabla u)|^2 \leq \frac{1}{|B_\sigma|} \int_{B_\sigma} |E(\nabla u)|^2 \leq \frac{|\mathcal{B}_K|}{|B_\sigma| |\mathcal{B}_K|} \int_{B_K} |E(\nabla u)|^2 \leq \frac{\theta_\mathcal{I}}{|\mathcal{B}_K|} \int_{B_K} |E(\nabla u)|^2. $$

As:

$$ \frac{m_K}{|\mathcal{B}_K|} = \frac{m_K}{|B_K \cap \Omega|} \frac{|B_K \cap \Omega|}{|\mathcal{B}_K|} \leq \theta_\mathcal{T} \frac{|B_K|}{|\mathcal{B}_K|} \leq \theta_\mathcal{T} \kappa_\psi^d, $$

and consequently $|R_4| \leq \eta_\psi^{1/2} \theta_\mathcal{T}^{1/2} \theta_\mathcal{A} \kappa_\psi^{d/2} h^p (h |\mathbf{V}|_X + |\mathbf{V}|_0) C_{E,0} ||\nabla u||_{L^2(\Omega)}$. Finally, we have using Cauchy-Schwarz inequality that:

$$ |R_5| \leq \theta_A h^p |\mathbf{V}|_0 \left( \sum_{K \in \mathcal{T}} m_K |G_K|^2 \right)^{1/2} \quad \text{and} \quad |G_K|^2 \leq \frac{1}{|B_K \cap \Omega|} \int_{B_K \cap \Omega} |E(\nabla u)|^2, $$

and thus:

$$ |R_5| \leq \theta_A h^p |\mathbf{V}|_0 \left( \sum_{K \in \mathcal{T}} m_K |G_K|^2 \right)^{1/2} \leq \theta_A \theta_\mathcal{T}^{1/2} h^p |\mathbf{V}|_0 \left( \sum_{K \in \mathcal{T}} \int_{B_K \cap \Omega} |E(\nabla u)|^2 \right)^{1/2}, $$

and consequently $|R_5| \leq \eta_N \theta_\mathcal{A} \theta_\mathcal{T}^{1/2} h^p |\mathbf{V}|_0 C_{E,0} ||\nabla u||_{L^2(\Omega)}$. Thus, there exists $C > 0$ such that $|\mathcal{R}(\mathbf{V})| \leq C(h + h^p)||\mathbf{V}||_X$. Using (19)-(20)-(23), we get:

$$ ||\mathbf{U} - \mathcal{I}^0(\mathbf{u})||_X^2 \leq C a_h (\mathbf{U} - \mathcal{I}^0(\mathbf{u}), \mathbf{U} - \mathcal{I}^0(\mathbf{u})), $$

and thus taking $\mathbf{V} = \mathbf{U} - \mathcal{I}^0(\mathbf{u})$, we obtain $||\mathbf{U} - \mathcal{I}^0(\mathbf{u})||_X \leq C(h + h^p)$. Finally remark that using the triangular inequality, we have:

$$ ||u - \Pi_\mathcal{T}(\mathbf{U})||_{L^2(\Omega)}^2 \leq ||u - \Pi_\mathcal{T}(\mathcal{I}^0(\mathbf{u}))||_{L^2(\Omega)}^2 + ||\Pi_\mathcal{T}(\mathcal{I}^0(\mathbf{u})) - \Pi_\mathcal{T}(\mathbf{U})||_{L^2(\Omega)}^2, $$

and

$$ ||\nabla u - \nabla_\mathcal{T}(\mathbf{U})||_{L^2(\Omega)}^2 \leq ||\nabla u - \nabla_\mathcal{T}(\mathcal{I}^0(\mathbf{u}))||_{L^2(\Omega)}^2 + ||\nabla_\mathcal{T}(\mathcal{I}^0(\mathbf{u})) - \nabla_\mathcal{T}(\mathbf{U})||_{L^2(\Omega)}^2, $$

and that lemma 3.1, and (19)-(20)-(23) gives us:

$$ ||\Pi_\mathcal{T}(\mathbf{U} - \mathcal{I}^0(\mathbf{u}))||_{L^2(\Omega)} \leq C ||\mathbf{U} - \mathcal{I}^0(\mathbf{u})||_X \quad \text{and} \quad ||\nabla_\mathcal{T}(\mathbf{U} - \mathcal{I}^0(\mathbf{u}))||_{L^2(\Omega)} \leq C ||\mathbf{U} - \mathcal{I}^0(\mathbf{u})||_X $$

Combining the above results with the interpolation results (48) and (50) consequently gives the desired estimates. To obtain the local estimates, remark that:

$$ ||u - \mathcal{M}_K(\mathbf{U})||_{L^2(B_K \cap \Omega)}^2 \leq ||u - \mathcal{M}_K(\mathcal{I}^0(\mathbf{u}))||_{L^2(B_K \cap \Omega)}^2 + ||\mathcal{M}_K(\mathcal{I}^0(\mathbf{u})) - \mathcal{M}_K(\mathbf{U})||_{L^2(B_K \cap \Omega)}^2, $$

and that:

$$ ||\mathcal{M}_K(\mathcal{I}^0(\mathbf{u})) - \mathcal{M}_K(\mathbf{U})||_{L^2(B_K \cap \Omega)}^2 = \int_{B_K \cap \Omega} |\mathcal{M}_K(\mathcal{I}^0(\mathbf{u})) - \mathcal{M}_K(\mathbf{U})|^2 = \frac{|B_K \cap \Omega|}{m_K} m_K |\mathcal{M}_K(\mathcal{I}^0(\mathbf{u})) - \mathcal{M}_K(\mathbf{U})|^2 $$
$$ \leq \theta_T |\mathcal{I}^0(u) - \mathcal{U}|_0^2 \leq \theta_T ||\mathcal{I}^0(u) - \mathcal{U}||_X^2 \leq C(h + h^p)^2, $$

and the result immediately follows from (43). In the same way, we have:

$$ \begin{aligned} & ||\nabla u - \nabla_K(\mathcal{U})||_{L^2(B_K \cap \Omega)^d}^2 \leq ||\nabla u - \nabla_K(\mathcal{I}^0(u))||_{L^2(B_K \cap \Omega)^d}^2 + ||\nabla_K(\mathcal{I}^0(u)) - \nabla_K(\mathcal{U})||_{L^2(B_K \cap \Omega)^d}^2 \\ & \leq ||\nabla u - \nabla_K(\mathcal{I}^0(u))||_{L^2(B_K \cap \Omega)^d}^2 + C \nabla \theta_T ||\mathcal{I}^0(u) - \mathcal{U}||_X^2, \end{aligned} $$

and the result follows from (44). □

A direct use of the estimates of proposition 5.4 would only provide a rate $h+h^p$ for the $L^2$ convergence of $\Pi_N(\mathcal{U})$ towards $u$. However, the estimates of propositions 5.1 and 5.2 suggest that one could achieve $\min(h^2, h^{p+1})$, and thus $L^2$ superconvergence if $p \geq 1$. This is moreover what is observed in practice (see [7]). However, the usual duality argument that is expected to lead to such a result is difficult to apply in our context. For this reason, we do not wish to elaborate any further on optimal $L^2$ convergence rates here. The above result is also sub-optimal in the sense that we require $u \in H^3(\Omega)$ instead of the usual $H^2(\Omega)$. This is due to the fact that in the above proof we use the strong form of the Poisson problem and the local consistency of the discrete divergence operator $DIV_K$ applied to $\nabla u$, which is probably sub-optimal. We nevertheless hope that the available results emphasize enough the link between quality parameters, geometrical approximation order and convergence rates.

# 6 CONCLUSION AND PERSPECTIVES

On the simplest possible model problem, we established convergence results and error estimates for the network element method. The error estimates are slightly sub-optimal as they require a solution belonging to $H^3$. The natural extension to heterogeneous and anisotropic diffusion tensors and reaction coefficients of the method and the associated convergence results will be the subject of a future paper. The results presented here could be improved in two ways: first by establishing a more explicit bound on the parameters of the quadrature family, probably through an estimation of $\theta_\psi$ using netowrk quality parameters, and secondly by establishing error estimates with $H^2$ regularity instead of $H^3$.

# REFERENCES

[1] R. Adams and J. Fournier. Sobolev Spaces, 2nd Edition. Elsevier, 2003.

[2] I. Babuška, U. Banerjee, and J.E. Osborn. Survey of meshless and generalized finite element methods: A unified approach. Acta Numerica, Vol.12, pp. 1-125, 2003.

[3] M. Botti, D.A. Di Pietro, and P. Sochala. A hybrid high-order method for nonlinear elasticity. SIAM J. Numer. Anal., 2017, Vol. 55(6), pp. 2687-2717, 2018.

[4] S. Brenner and R. Scott. The mathematical theory of finite element methods. Springer, 3rd ed., 2008.

[5] F. Chave, D.A. Di Pietro, and L. Formaggia. A hybrid high-order method for darcy flows in fractured porous media. SIAM J. Sci. Comput., Vol. 40(2), pp. 1063-1094, 2018.

[6] J.-S. Chen, M. Hillman, and S.-W Chi. Meshfree methods: Progress made after 20 years. Journal of Engineering Mechanics, Vol. 143(4), 2017.

[7] J. Coatléven. Principles of a network element method. Journal of Computational Physics, Vol. 433, Art. 110197, 2021.

[8] L. Beirao da Veiga, F. Brezzi, A. Cangiani, G. Manzini, L.D. Marini, and A. Russo. Basic principles of virtual element methods. Mathematical Models and Methods in Applied Sciences, Vol. 23(1), pp. 199-214, 2013.

[9] L. Beirao da Veiga, F. Brezzi, and L. D. Marini. Virtual elements for linear elasticity problems. SIAM J. Numer. Anal., Vol. 51(2), pp.794-812, 2013.

[10] L. Beirao da Veiga, K. Lipnikov, and G. Manzini. The Mimetic Finite Difference Method for Elliptic Problems. Springer, 2014.

[11] L. Beirao da Veiga, C. Lovadina, and D. Mora. A virtual element method for elastic and inelastic problems on polytope meshes. Computer Methods in Applied Mechanics and Engineering, Vol. 295, pp. 327-346, 2015.
[12] L. Beirao da Veiga, C. Lovadina, and G. Vacca. Divergence free virtual elements for the stokes problem on polygonal meshes. *ESAIM Math. Model. Numer. Anal., Vol. 51*(2), pp. 509-535, 2017.

[13] D. A. Di Pietro and A. Ern. *Mathematical aspects of discontinuous Galerkin methods*. Springer, 2012.

[14] O. Diyankov. Uncertain grid method for numerical solution of pdes. Technical report, NeurOK Software, 2008.

[15] J. Droniou, R. Eymard, T. Gallouët, C. Guichard, and R. Herbin. *The Gradient Discretisation Method*. Springer, 2018.

[16] R. Eymard, T. Gallouët, and R. Herbin. *Finite volume methods, In P.G. Ciarlet and J.-L. Lions editors, Techniques of scientific computing, Part III, Handbook of Numerical Analysis*, pp. 713-1020. North-Holland, Amsterdam, 2000.

[17] R. Eymard, T. Gallouët, and R. Herbin. Discretisation of heterogeneous and anisotropic diffusion problems on general nonconforming meshes sushi : a scheme using stabilisation and hybrid interfaces. *IMA J. Num. Anal.*, 30(4), pp 1009-1043, 2010.

[18] P. Grisvard. *Elliptic Problems in Nonsmooth Domains*. Pitman Publishing Inc, Massachusetts, 1985.

[19] A. Katz and A. Jameson. Edge-based meshless methods for compressible viscous flow with applications to overset grids. *Proceedings of the 38th Fluid Dynamics Conference and Exhibit, American Institute of Aeronautics and Astronautics*, 2008.

[20] A. Katz and A. Jameson. A meshless volume scheme. *Proceedings of 19th AIAA Computational Fluid Dynamics, Fluid Dynamics and Co-located Conferences*, pp. 2009-3534, American Institute of Aeronautics and Astronautics, 2009.

[21] J.M. Melenk. *On Approximation in Meshless Methods*, pages 65–141. Springer Berlin Heidelberg, Berlin, Heidelberg, 2005.

[22] J.M. Melenk and I. Babuška. The partition of unity finite element method: Basic theory and applications. *Computer Methods in Applied Mechanics and Engineering*, Vol. 139(1-4), pp. 289-314 Cambridge University Press, 1996.

[23] D.A. Di Pietro and A. Ern. A hybrid high-order locking-free method for linear elasticity on general meshes. *Comput. Meth. Appl. Mech. Engrg.*, Vol. 283, pp. 1-21, 2015.

[24] D.A. Di Pietro and A. Ern. Hybrid high-order methods for variable-diffusion problems on general meshes. *C. R. Acad. Sci. Paris, Ser. I*, Vol. 353, pp. 31-34, 2015.

[25] E.M. Stein. *Singular integrals and differentiability properties of functions*. Princeton University Press, 1970.

[26] N. Trask, P. Bochev, and M. Perego. A conservative, consistent, and scalable mesh-free mimetic method. *Journal of Computational Physics*, vol. 409 pp. 109-187, 2020.

[27] N. Trask, M. Perego, and P. Bochev. A high-order staggered meshless method for elliptic problems. *SIAM Journal on Scientific Computing*, Vol. 39(2), pp. 479-502, 2017.

[28] G. Vacca and L. Beirao da Veiga. Virtual element methods for parabolic problems on polygonal meshes. *Numer Methods Partial Differential Eq*, Vol.31(6), pp. 2110-2134, 2015.

[29] E. Kwan yu Chiu, Q. Wang, R. Hu, and A. Jameson. A conservative mesh-free scheme and generalized framework for conservation laws. *SIAM J. Sci. Comput.*, Vol. 34(6), pp.2896-2916, 2012.