# Network element methods for linear elasticity

Julien Coatléven

► To cite this version:

Julien Coatléven. Network element methods for linear elasticity. Comptes Rendus. Mécanique, 2023, 351 (S1), pp.331-356. 10.5802/crmeca.231. hal-04618647

HAL Id: hal-04618647

https://ifp.hal.science/hal-04618647

Submitted on 20 Jun 2024

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
# Network element methods for linear elasticity

Julien Coatléven

► To cite this version:

Julien Coatléven. Network element methods for linear elasticity. Comptes Rendus. Mécanique, 2023, 351 (S1), pp.331-356. 10.5802/crmeca.231. hal-04618647

HAL Id: hal-04618647

https://ifp.hal.science/hal-04618647v1

Submitted on 20 Jun 2024

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
# Comptes Rendus
## Mécanique

Julien Coatléven

**Network element methods for linear elasticity**

Volume 351, Special Issue S1 (2023), p. 331-356

Online since: 11 December 2023

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

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

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

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

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

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

This article is licensed under the

CREATIVE COMMONS ATTRIBUTION 4.0 INTERNATIONAL LICENSE.

http://creativecommons.org/licenses/by/4.0/
# Comptes Rendus
## Mécanique

Julien Coatléven

**Network element methods for linear elasticity**

Volume 351, Special Issue S1 (2023), p. 331-356

Online since: 11 December 2023

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

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

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

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

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

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

This article is licensed under the

CREATIVE COMMONS ATTRIBUTION 4.0 INTERNATIONAL LICENSE.

http://creativecommons.org/licenses/by/4.0/
differences (GFD, [6]), moving least-square based methods (see [7–9]) or methods based on radial basis functions (see [7, 10–12]). In contrast, methods based on meshfree basis functions such as the diffuse element method (see [13]), the element free Galerkin method (see [14, 15]), partition of unity based methods (see [16, 17]) theoretically inherit the built-in stability of weak form variational methods, however numerical integration of those basis functions can severely deteriorate their theoretical stability (see [7, 16]). The popular meshless local Petrov-Galerkin (MLPG, see [18, 19]) is an attempt to circumvent those numerical instabilities. For a general review on meshless methods and of recent progress made on them, we refer the reader to [20].

The recently introduced network element method (NEM) of [21] is a meshfree variational method that belongs to a third, less developped category of meshfree methods based soleley on a point cloud and a point connectivity, without any mesh at any stage of the method, not even for numerical integration. This point cloud and its connectivity of course form the network giving its name to the method. The literature on network based methods being mainly focused on finite volume like, strong form methods (see [22–27]), the NEM represents an attractive alternative as like methods based on meshfree basis functions, it inherits the stability of the variational approach. Moreover, thanks to its network based formulation, no numerical integration is required and thus stability is never a numerical issue (see [21]). Introduced on the Poisson problem, it has since been extended to heterogeneous and anisotropic diffusion reaction problems, either in non-conservative (see [28]) or conservative (see [29]) formulations. A convergence theory as well as error estimates were presented in [30]. The main practical challenge originally posed by the NEM was the computation of the so-called network geometry, i.e. of the family of weights that allows to derive polynomial reconstruction operators from the degrees of freedom. The original method for geometry generation of [21] was non-linear, however a linear bypass has been proposed in [31] along with fast network generation algorithms giving birth to a sufficiently fast network element workflow that competes with other meshfree methods.

The objective of the present article is the extension of the network element method to the more complex system of linear elasticity. As the network underlying the NEM generalizes in some sense the notion of mesh, if applied on a mesh it will also offer alternatives to existing first order methods for linear elasticity on mesh with general cell geometries, such as the first order virtual element method (see [32, 33]) or hybrid finite volume like methods (see [34–36]). Without any attempt at being exhaustive, let us nevertheless mention that linear elasticity problems have been of course investigated with both strong and weak form meshless methods (see [9, 37–39]), up to considering fretting [40] or fracture [41, 42] problems. For a quite recent review we refer the reader to [43]. All those methods are potentially subject to some numerical instabilities: the main advantage of the NEM for linear elasticity is, as in the case of diffusion problems, that its variational, numerical integration free nature will provide stability in the classical way: through a discrete counterpart of Korn's inequality.

The paper will be organized as follows: after recalling the notions of discretization network and associated network geometry on which the NEM is built, we propose a sufficient condition on networks that allow to establish a discrete counterpart of Korn's inequality. Building on it for the derivation of the stabilization term of the method, we then present a network element method for linear elasticity and establish its well-posedness. We then also present a flux based version that is conservative in terms of numerical tractions as well as mesh-based implementations of both the non-conservative and conservative NEMs for linear elasticity, that are interesting per se as they coincide with or are alternatives to existing mesh-based numerical methods. Finally, we present some numerical results on challenging problems of the literature, that illustrates the good behavior of the NEM for linear elasticity, even in the quasi-incompressible limit.
differences (GFD, [6]), moving least-square based methods (see [7–9]) or methods based on radial basis functions (see [7, 10–12]). In contrast, methods based on meshfree basis functions such as the diffuse element method (see [13]), the element free Galerkin method (see [14, 15]), partition of unity based methods (see [16, 17]) theoretically inherit the built-in stability of weak form variational methods, however numerical integration of those basis functions can severely deteriorate their theoretical stability (see [7, 16]). The popular meshless local Petrov–Galerkin (MLPG, see [18, 19]) is an attempt to circumvent those numerical instabilities. For a general review on meshless methods and of recent progress made on them, we refer the reader to [20].

The recently introduced network element method (NEM) of [21] is a meshfree variational method that belongs to a third, less developped category of meshfree methods based soleley on a point cloud and a point connectivity, without any mesh at any stage of the method, not even for numerical integration. This point cloud and its connectivity of course form the network giving its name to the method. The literature on network based methods being mainly focused on finite volume like, strong form methods (see [22–27]), the NEM represents an attractive alternative as like methods based on meshfree basis functions, it inherits the stability of the variational approach. Moreover, thanks to its network based formulation, no numerical integration is required and thus stability is never a numerical issue (see [21]). Introduced on the Poisson problem, it has since been extended to heterogeneous and anisotropic diffusion reaction problems, either in non-conservative (see [28]) or conservative (see [29]) formulations. A convergence theory as well as error estimates were presented in [30]. The main practical challenge originally posed by the NEM was the computation of the so-called network geometry, i.e. of the family of weights that allows to derive polynomial reconstruction operators from the degrees of freedom. The original method for geometry generation of [21] was non-linear, however a linear bypass has been proposed in [31] along with fast network generation algorithms giving birth to a sufficiently fast network element workflow that competes with other meshfree methods.

The objective of the present article is the extension of the network element method to the more complex system of linear elasticity. As the network underlying the NEM generalizes in some sense the notion of mesh, if applied on a mesh it will also offer alternatives to existing first order methods for linear elasticity on mesh with general cell geometries, such as the first order virtual element method (see [32, 33]) or hybrid finite volume like methods (see [34–36]). Without any attempt at being exhaustive, let us nevertheless mention that linear elasticity problems have been of course investigated with both strong and weak form meshless methods (see [9, 37–39]), up to considering fretting [40] or fracture [41, 42] problems. For a quite recent review we refer the reader to [43]. All those methods are potentially subject to some numerical instabilities: the main advantage of the NEM for linear elasticity is, as in the case of diffusion problems, that its variational, numerical integration free nature will provide stability in the classical way: through a discrete counterpart of Korn's inequality.

The paper will be organized as follows: after recalling the notions of discretization network and associated network geometry on which the NEM is built, we propose a sufficient condition on networks that allow to establish a discrete counterpart of Korn's inequality. Building on it for the derivation of the stabilization term of the method, we then present a network element method for linear elasticity and establish its well-posedness. We then also present a flux based version that is conservative in terms of numerical tractions as well as mesh-based implementations of both the non-conservative and conservative NEMs for linear elasticity, that are interesting per se as they coincide with or are alternatives to existing mesh-based numerical methods. Finally, we present some numerical results on challenging problems of the literature, that illustrates the good behavior of the NEM for linear elasticity, even in the quasi-incompressible limit.
## 2. Discretization networks and approximate geometries

In this section, we recall for the reader's convenience the notions of discretization network and associated approximate network geometry of [31] that where originally introduced in [21].

Let $\Omega$ be an open bounded connected subset of $\mathbb{R}^d$, $d \in \mathbb{N} \setminus \{0\}$, with boundary $\partial\Omega = \bar{\Omega} \setminus \Omega$ 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 open ball of radius $r$ centered at $\mathbf{x}$ for the usual euclidean norm $|\mathbf{x}|^2 = \sum_{i=1}^d x_i^2$. Following [21, 25-27], a discretization network $\mathcal{N}$ of $\Omega$ is defined from two sets of points $\mathcal{P}_{\mathcal{F}}$ and $\mathcal{P}_{\mathcal{F}'}$, by setting $\mathcal{N} = \{\mathcal{F}, \mathcal{F}'\}$:

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

(H2) 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{F}$. 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 \Omega$.

(H3) $\bar{\Omega} \subset \cup_{K \in \mathcal{F}} B(\mathbf{x}_K, r_K)$. For all $(K_1, K_2) \in \mathcal{F}^2$ such that $K_1 \neq K_2$, $\mathbf{x}_{K_1} \neq \mathbf{x}_{K_2}$.

(H4) To any $K \in \mathcal{F}$, we associate a Lipschitz open set $\Omega_K \subset B(\mathbf{x}_K, r_K)$ such that $\bar{\Omega} \subset \cup_{K \in \mathcal{F}} \bar{\Omega}_K$. For any $K \in \mathcal{F}$ such that $\partial\bar{\Omega} \cap \Omega_K \neq \emptyset$, there exists $\sigma \in \mathcal{F}_{\text{ext}}$ such that $K \in \mathcal{T}_\sigma$. For any $(K, L) \in \mathcal{F}^2$ such that $\bar{\Omega}_K \cap \bar{\Omega}_L \neq \emptyset$, there exists $\sigma \in \mathcal{F}$ such that $(K, L) \subset \mathcal{T}_\sigma$.

For any $K \in \mathcal{F}$, we denote $\mathcal{T}_K = \{\sigma \in \mathcal{F} \mid K \in \mathcal{T}_\sigma\}$ the set of interfaces of $K$, implying that for any $\sigma \in \mathcal{F}$ the set $\mathcal{T}_\sigma$ denotes the cells connected to the interface $\sigma$ and satisfies $\mathcal{T}_\sigma = \{K \in \mathcal{F} \mid \sigma \in \mathcal{T}_K\}$. We denote $B_K = B(\mathbf{x}_K, r_K)$, $h = \max_{K \in \mathcal{F}} 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{F}$ the set $(\mathbf{x}_\sigma)_{\sigma \in \mathcal{T}_K}$ is unisolvent for first order polynomials.

Important and obvious examples of discretization networks are mesh based discretization networks. Indeed, if $\mathcal{M}$ is a mesh, if we identify $\mathcal{F}$ with the set of cells of the of mesh $\mathcal{M}$ and $\mathcal{F}$ with the set of its faces (resp. its vertices), then choosing $\mathbf{x}_K$ to be the barycenter of each cell $K \in \mathcal{F}$ and $\mathbf{x}_\sigma$ to be the barycenter of each face $\sigma \in \mathcal{F}$ (resp. the position of each vertex), and finally choosing for each cell $K$ the open set $\Omega_K$ to coincide with cell $K$ we get an admissible discretization network.

A network geometry is defined as a set of coefficients:

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

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

$$ m_K > 0 \quad \text{for all } K \in \mathcal{F} \quad \text{and} \quad \sum_{K \in \mathcal{F}} m_K = |\Omega|, \qquad (1) $$

where $|\Omega|$ denotes the Lebesgue measure of $\Omega$, 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{F}, \forall 1 \le i \le d, \qquad (2) $$

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{F}, \forall 1 \le i, j \le d, \qquad (3) $$

where $\delta_{ij}$ is the Kronecker symbol, 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 (4) $$
## 2. Discretization networks and approximate geometries

In this section, we recall for the reader's convenience the notions of discretization network and associated approximate network geometry of [31] that where originally introduced in [21].

Let $\Omega$ be an open bounded connected subset of $\mathbb{R}^d$, $d \in \mathbb{N} \setminus \{0\}$, with boundary $\partial\Omega = \bar{\Omega} \setminus \Omega$ 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 open ball of radius $r$ centered at $\mathbf{x}$ for the usual euclidean norm $|\mathbf{x}|^2 = \sum_{i=1}^d x_i^2$. Following [21, 25-27], a discretization network $\mathcal{N}$ of $\Omega$ is defined from two sets of points $\mathcal{P}_{\mathcal{F}}$ and $\mathcal{P}_{\mathcal{F}'}$, by setting $\mathcal{N} = \{\mathcal{F}, \mathcal{F}'\}$:

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

(H2) 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{F}$. 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 \Omega$.

(H3) $\bar{\Omega} \subset \cup_{K \in \mathcal{F}} B(\mathbf{x}_K, r_K)$. For all $(K_1, K_2) \in \mathcal{F}^2$ such that $K_1 \neq K_2$, $\mathbf{x}_{K_1} \neq \mathbf{x}_{K_2}$.

(H4) To any $K \in \mathcal{F}$, we associate a Lipschitz open set $\Omega_K \subset B(\mathbf{x}_K, r_K)$ such that $\bar{\Omega} \subset \cup_{K \in \mathcal{F}} \bar{\Omega}_K$. For any $K \in \mathcal{F}$ such that $\partial\bar{\Omega} \cap \Omega_K \neq \emptyset$, there exists $\sigma \in \mathcal{F}_{\text{ext}}$ such that $K \in \mathcal{T}_\sigma$. For any $(K, L) \in \mathcal{F}^2$ such that $\bar{\Omega}_K \cap \bar{\Omega}_L \neq \emptyset$, there exists $\sigma \in \mathcal{F}$ such that $(K, L) \subset \mathcal{T}_\sigma$.

For any $K \in \mathcal{F}$, we denote $\mathcal{T}_K = \{\sigma \in \mathcal{F} \mid K \in \mathcal{T}_\sigma\}$ the set of interfaces of $K$, implying that for any $\sigma \in \mathcal{F}$ the set $\mathcal{T}_\sigma$ denotes the cells connected to the interface $\sigma$ and satisfies $\mathcal{T}_\sigma = \{K \in \mathcal{F} \mid \sigma \in \mathcal{T}_K\}$. We denote $B_K = B(\mathbf{x}_K, r_K)$, $h = \max_{K \in \mathcal{F}} 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{F}$ the set $(\mathbf{x}_\sigma)_{\sigma \in \mathcal{T}_K}$ is unisolvent for first order polynomials.

Important and obvious examples of discretization networks are mesh based discretization networks. Indeed, if $\mathcal{M}$ is a mesh, if we identify $\mathcal{F}$ with the set of cells of the of mesh $\mathcal{M}$ and $\mathcal{F}$ with the set of its faces (resp. its vertices), then choosing $\mathbf{x}_K$ to be the barycenter of each cell $K \in \mathcal{F}$ and $\mathbf{x}_\sigma$ to be the barycenter of each face $\sigma \in \mathcal{F}$ (resp. the position of each vertex), and finally choosing for each cell $K$ the open set $\Omega_K$ to coincide with cell $K$ we get an admissible discretization network.

A network geometry is defined as a set of coefficients:

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

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

$$ m_K > 0 \quad \text{for all } K \in \mathcal{F} \quad \text{and} \quad \sum_{K \in \mathcal{F}} m_K = |\Omega|, \qquad (1) $$

where $|\Omega|$ denotes the Lebesgue measure of $\Omega$, 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{F}, \forall 1 \le i \le d, \qquad (2) $$

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{F}, \forall 1 \le i, j \le d, \qquad (3) $$

where $\delta_{ij}$ is the Kronecker symbol, 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 (4) $$
The geometric approximation error is measured through the constants $\theta_{\mathcal{A}} > 0$ and $p \ge 1$ such that:

$$|\epsilon_K^{0,i}| \le \theta_{\mathcal{A}} h_K^p \quad \forall K \in \mathcal{T}, \forall 1 \le i \le d, \quad \text{and} \quad |\epsilon_K^{1,ij}| \le \theta_{\mathcal{A}} h_K^p \quad \forall K \in \mathcal{T}, \forall 1 \le i, j \le d, \quad (5)$$

and

$$|\epsilon_{\sigma}^{i}| \le \theta_{\mathcal{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. \quad (6)$$

We say that a network geometry is admissible if it satisfies (2)–(4) and the family of measures is admissible.

### 3. A basic network element method for linear elasticity

The aim of this section is to introduce the network element method discretization of a basic linear elasticity problem with constant Lamé coefficients. After introducing our model problem, we detail the set of degrees of freedom and the discrete reconstruction operators associated with the network element method. We then propose a discrete variational formulation, and next establish a discrete Korn inequality that naturally induces the stability of our variational formulation. We conclude this section by providing implementation details.

#### 3.1. Model problem

We consider the simplest possible linear elasticity model problem, i.e. the case of constant Lamé coefficients, set on $\Omega$:

$$-\mathrm{div} \sigma(\mathbf{u}) = \mathbf{f}, \quad (7)$$

with $\mathbf{f} \in L^2(\Omega)^d$ and

$$\sigma(\mathbf{u}) = 2\mu\epsilon(\mathbf{u}) + \lambda \operatorname{div} \mathbf{u} \mathbf{I}_d \quad \text{and} \quad \epsilon(\mathbf{u}) = \frac{1}{2}(\nabla \mathbf{u} + \nabla^T \mathbf{u}) \quad \text{and} \quad \operatorname{div} \sigma = \sum_{i=1}^{d} \sum_{j=1}^{d} \partial_{x_j} \sigma_{ij} e_i, \quad (8)$$

where $\mathbf{I}_d$ is the identity matrix in dimension $d$. For the sake of simplicity again, we complement it with homogeneous Dirichlet boundary conditions, i.e. $\mathbf{u} = \mathbf{0}$ on $\partial\Omega$. The weak solution associated to (7)–(8) is the unique $\mathbf{u} \in H_0^1(\Omega)^d$ such that

$$a(\mathbf{u}, \mathbf{v}) = 2\mu \int_{\Omega} \epsilon(\mathbf{u}) : \epsilon(\mathbf{v}) + \lambda \int_{\Omega} \operatorname{div} \mathbf{u} \operatorname{div} \mathbf{v} = \int_{\Omega} \mathbf{f} \cdot \mathbf{v} \quad \text{for all } \mathbf{v} \in H_0^1(\Omega)^d, \quad (9)$$

where we recall that $\mathbf{A}: \mathbf{B} = \operatorname{tr}(\mathbf{A}^T \mathbf{B}) = \sum_{i=1}^{d} \sum_{j=1}^{d} A_{ij} B_{ij}$.

#### 3.2. Degrees of freedom and network element operators

Let $\mathcal{N}$ be an admissible discretization network, endowed with an admissible network geometry $\mathcal{G}$. To $\mathcal{N}$ is associated the following space of degrees of freedom, which is simply a vector valued version of the original NEM dof space of [21]:

$$X_{\mathcal{N}}^{d} = \{(\mathbf{u}_{\sigma})_{\sigma \in \mathcal{F}} \mid \mathbf{u}_{\sigma} \in \mathbb{R}^{d} \ \forall \sigma \in \mathcal{F}\},$$

and denote $\mathbf{U} = (\mathbf{u}_{\sigma})_{\sigma \in \mathcal{F}}$. To account for homogeneous Dirichlet boundary condition, we also consider:

$$X_{\mathcal{N},0}^{d} = \{ \mathbf{U} \in X_{\mathcal{N}}^{d} \mid \mathbf{u}_{\sigma} = \mathbf{0} \text{ for all } \sigma \in \mathcal{F}_{\text{ext}} \},$$

and we define the local set of degrees of freedom associated to a cell by:

$$X_{\mathcal{N},K}^{d} = \{(\mathbf{u}_{\sigma})_{\sigma \in \mathcal{F}_{K}} \mid \mathbf{u}_{\sigma} \in \mathbb{R}^{d} \ \forall \sigma \in \mathcal{F}_{K}\}.$$
The geometric approximation error is measured through the constants $\theta_{\mathcal{A}} > 0$ and $p \ge 1$ such that:

$$|\epsilon_K^{0,i}| \le \theta_{\mathcal{A}} h_K^p \quad \forall K \in \mathcal{T}, \forall 1 \le i \le d, \quad \text{and} \quad |\epsilon_K^{1,ij}| \le \theta_{\mathcal{A}} h_K^p \quad \forall K \in \mathcal{T}, \forall 1 \le i, j \le d, \quad (5)$$

and

$$|\epsilon_{\sigma}^{i}| \le \theta_{\mathcal{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. \quad (6)$$

We say that a network geometry is admissible if it satisfies (2)–(4) and the family of measures is admissible.

### 3. A basic network element method for linear elasticity

The aim of this section is to introduce the network element method discretization of a basic linear elasticity problem with constant Lamé coefficients. After introducing our model problem, we detail the set of degrees of freedom and the discrete reconstruction operators associated with the network element method. We then propose a discrete variational formulation, and next establish a discrete Korn inequality that naturally induces the stability of our variational formulation. We conclude this section by providing implementation details.

#### 3.1. Model problem

We consider the simplest possible linear elasticity model problem, i.e. the case of constant Lamé coefficients, set on $\Omega$:

$$-\mathrm{div} \sigma(\mathbf{u}) = \mathbf{f}, \quad (7)$$

with $\mathbf{f} \in L^2(\Omega)^d$ and

$$\sigma(\mathbf{u}) = 2\mu\epsilon(\mathbf{u}) + \lambda \operatorname{div} \mathbf{u} \mathbf{I}_d \quad \text{and} \quad \epsilon(\mathbf{u}) = \frac{1}{2}(\nabla \mathbf{u} + \nabla^T \mathbf{u}) \quad \text{and} \quad \operatorname{div} \sigma = \sum_{i=1}^{d} \sum_{j=1}^{d} \partial_{x_j} \sigma_{ij} \mathbf{e}_i, \quad (8)$$

where $\mathbf{I}_d$ is the identity matrix in dimension $d$. For the sake of simplicity again, we complement it with homogeneous Dirichlet boundary conditions, i.e. $\mathbf{u} = \mathbf{0}$ on $\partial\Omega$. The weak solution associated to (7)–(8) is the unique $\mathbf{u} \in H_0^1(\Omega)^d$ such that

$$a(\mathbf{u}, \mathbf{v}) = 2\mu \int_{\Omega} \epsilon(\mathbf{u}) : \epsilon(\mathbf{v}) + \lambda \int_{\Omega} \operatorname{div} \mathbf{u} \operatorname{div} \mathbf{v} = \int_{\Omega} \mathbf{f} \cdot \mathbf{v} \quad \text{for all } \mathbf{v} \in H_0^1(\Omega)^d, \quad (9)$$

where we recall that $\mathbf{A}: \mathbf{B} = \operatorname{tr}(\mathbf{A}^T \mathbf{B}) = \sum_{i=1}^{d} \sum_{j=1}^{d} A_{ij} B_{ij}$.

#### 3.2. Degrees of freedom and network element operators

Let $\mathcal{N}$ be an admissible discretization network, endowed with an admissible network geometry $\mathcal{G}$. To $\mathcal{N}$ is associated the following space of degrees of freedom, which is simply a vector valued version of the original NEM dof space of [21]:

$$X_{\mathcal{N}}^{d} = \{(\mathbf{u}_{\sigma})_{\sigma \in \mathcal{F}} \mid \mathbf{u}_{\sigma} \in \mathbb{R}^{d} \ \forall \sigma \in \mathcal{F}\},$$

and denote $\mathbf{U} = (\mathbf{u}_{\sigma})_{\sigma \in \mathcal{F}}$. To account for homogeneous Dirichlet boundary condition, we also consider:

$$X_{\mathcal{N},0}^{d} = \{ \mathbf{U} \in X_{\mathcal{N}}^{d} \mid \mathbf{u}_{\sigma} = \mathbf{0} \text{ for all } \sigma \in \mathcal{F}_{\text{ext}} \},$$

and we define the local set of degrees of freedom associated to a cell by:

$$X_{\mathcal{N},K}^{d} = \{(\mathbf{u}_{\sigma})_{\sigma \in \mathcal{F}_{K}} \mid \mathbf{u}_{\sigma} \in \mathbb{R}^{d} \ \forall \sigma \in \mathcal{F}_{K}\}.$$
Of course, for any $\mathbf{U} \in X_{\mathcal{N}}^d$, we denote $\mathbf{U}_K = (\mathbf{u}_\sigma)_{\sigma \in \mathcal{F}_K}$. To any cell $K \in \mathcal{T}$, we associate the local $d$-dimensional reconstruction operator $\Pi_K^d$ defined by:

$$ \left| \begin{array}{l} \Pi_K^d : X_{\mathcal{N},K}^d \longrightarrow \mathbb{P}_1(\mathbb{R}^d)^d \\ \mathbf{U}_K \longrightarrow \Pi_K^d(\mathbf{U}_K) = \mathcal{M}_K^d(\mathbf{U}_K) + \nabla_K^d(\mathbf{U}_K)(\mathbf{x}-\mathbf{x}_K), \end{array} \right. \qquad (10) $$

where

$$ \left| \begin{array}{l} \nabla_K^d : X_{\mathcal{N},K}^d \longrightarrow \mathbb{P}_0(\mathbb{R}^d)^d \\ \mathbf{U}_K \longrightarrow \nabla_K^d(\mathbf{U}_K)_{ij} = \frac{1}{m_K} \sum_{\sigma \in \mathcal{F}_K} u_\sigma^i \eta_{K,\sigma}^j, \end{array} \right. \qquad (11) $$

and

$$ \mathcal{M}_K^d(\mathbf{U}_K) = \sum_{\sigma \in \mathcal{F}_K} \gamma_{K,\sigma} \mathbf{u}_\sigma, $$

with the $(\gamma_{K,\sigma})_{\sigma \in \mathcal{F}_K}$ form a barycentric interpolation of $\mathbf{x}_K$ from the $(\mathbf{x}_\sigma)_{\sigma \in \mathcal{F}_K}$:

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

We further define the discrete divergence operator by setting:

$$ \left| \begin{array}{l} \mathcal{D} \mathcal{F} \mathcal{V}_K : X_{\mathcal{N},K}^d \longrightarrow \mathbb{P}_0(\mathbb{R}^d) \\ \mathbf{U}_K \longrightarrow \mathcal{D} \mathcal{F} \mathcal{V}_K(\mathbf{U}_K) = \frac{1}{m_K} \sum_{i=1}^d \sum_{\sigma \in \mathcal{F}_K} u_\sigma^i \eta_{K,\sigma}^i. \end{array} \right. \qquad (12) $$

Notice that component-wise, the operators $\Pi_K^d$, $\nabla_K^d$ and $\mathcal{M}_K^d$ are identical to the usual scalar NEM operators $\Pi_K$, $\nabla_K$ and $\mathcal{M}_K$ of [21], i.e., denoting $\mathbf{U}_K^i = (u_\sigma^i)_{\sigma \in \mathcal{F}_K}$ and $\mathbf{U}^i = (u_\sigma^i)_{\sigma \in \mathcal{F}}$:

$$ \Pi_K^d(\mathbf{U}_K)^i = \Pi_K(\mathbf{U}_K^i) \quad \text{and} \quad \nabla_K^d(\mathbf{U}_K)^i = \nabla_K(\mathbf{U}_K^i) \quad \text{and} \quad \mathcal{M}_K^d(\mathbf{U}_K)^i = \mathcal{M}_K(\mathbf{U}_K^i). $$

Also notice that $\nabla \Pi_K^d(\mathbf{U}_K) = \nabla_K^d(\mathbf{U}_K)$ by construction, thus we will use both notations interchangeably. Finally, we denote:

$$ \boldsymbol{\epsilon}_K(\mathbf{U}_K) = \frac{1}{2} (\nabla_K^d(\mathbf{U}_K) + \nabla_K^d(\mathbf{U}_K)^T), \qquad (13) $$

which is of course intended to be an approximation of $\boldsymbol{\epsilon}(\mathbf{u})$.

### 3.3. Approximation by the network element method

The discretization by the network element method consists in finding a solution $\mathbf{U} \in X_{\mathcal{N},0}^d$ of

$$ a_h(\mathbf{U}, \mathbf{V}) = l_h(\mathbf{V}) \quad \text{for all } \mathbf{V} \in X_{\mathcal{N},0}^d, \qquad (14) $$

where for the right-hand side assuming that $f_K$ is an approximation of $f$ at $\mathbf{x}_K$ we define the linear form $l_h : X_{\mathcal{N}}^d \longrightarrow \mathbb{R}$ by setting:

$$ l_h(\mathbf{V}) = \sum_{K \in \mathcal{T}} m_K f_K \cdot \mathcal{M}_K^d(\mathbf{V}_K), \qquad (15) $$

while

$$ a_h(\mathbf{U}, \mathbf{V}) = \sum_{K \in \mathcal{T}} a_h^K(\mathbf{U}_K, \mathbf{V}_K) \quad \text{with} \quad a_h^K(\mathbf{U}_K, \mathbf{V}_K) = a_{\mu,h}^K(\mathbf{U}_K, \mathbf{V}_K) + a_{\lambda,h}^K(\mathbf{U}_K, \mathbf{V}_K), \qquad (16) $$

where the discrete bilinear form $a_{\mu,h}^K : X_{\mathcal{N},K}^d \times X_{\mathcal{N},K}^d \longrightarrow \mathbb{R}$ mimics the term $2\mu \int_{\Omega} \boldsymbol{\epsilon}(\mathbf{u}) : \boldsymbol{\epsilon}(\mathbf{v})$ and is given by

$$ a_{\mu,h}^K(\mathbf{U}_K, \mathbf{V}_K) = 2\mu m_K \boldsymbol{\epsilon}_K(\mathbf{U}_K) : \boldsymbol{\epsilon}_K(\mathbf{V}_K) + s^K(\mathbf{U}_K, \mathbf{V}_K), \qquad (17) $$
Of course, for any $\mathbf{U} \in X_{\mathcal{N}}^d$, we denote $\mathbf{U}_K = (\mathbf{u}_{\sigma})_{\sigma \in \mathcal{F}_K}$. To any cell $K \in \mathcal{T}$, we associate the local $d$-dimensional reconstruction operator $\Pi_K^d$ defined by:

$$ \left| \begin{array}{l} \Pi_K^d : X_{\mathcal{N},K}^d \longrightarrow \mathbb{P}_1(\mathbb{R}^d)^d \\ \mathbf{U}_K \longrightarrow \Pi_K^d(\mathbf{U}_K) = \mathcal{M}_K^d(\mathbf{U}_K) + \nabla_K^d(\mathbf{U}_K)(\mathbf{x}-\mathbf{x}_K), \end{array} \right. \qquad (10) $$

where

$$ \left| \begin{array}{l} \nabla_K^d : X_{\mathcal{N},K}^d \longrightarrow \mathbb{P}_0(\mathbb{R}^d)^d \\ \mathbf{U}_K \longrightarrow \nabla_K^d(\mathbf{U}_K)_{ij} = \frac{1}{m_K} \sum_{\sigma \in \mathcal{F}_K} u_{\sigma}^i \eta_{K,\sigma}^j, \end{array} \right. \qquad (11) $$

and

$$ \mathcal{M}_K^d(\mathbf{U}_K) = \sum_{\sigma \in \mathcal{F}_K} \gamma_{K,\sigma} \mathbf{u}_{\sigma}, $$

with the $(\gamma_{K,\sigma})_{\sigma \in \mathcal{F}_K}$ form a barycentric interpolation of $\mathbf{x}_K$ from the $(\mathbf{x}_{\sigma})_{\sigma \in \mathcal{F}_K}$:

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

We further define the discrete divergence operator by setting:

$$ \left| \begin{array}{l} \mathcal{D} \mathcal{F} \mathcal{V}_K : X_{\mathcal{N},K}^d \longrightarrow \mathbb{P}_0(\mathbb{R}^d) \\ \mathbf{U}_K \longrightarrow \mathcal{D} \mathcal{F} \mathcal{V}_K(\mathbf{U}_K) = \frac{1}{m_K} \sum_{i=1}^d \sum_{\sigma \in \mathcal{F}_K} u_{\sigma}^i \eta_{K,\sigma}^i. \end{array} \right. \qquad (12) $$

Notice that component-wise, the operators $\Pi_K^d$, $\nabla_K^d$ and $\mathcal{M}_K^d$ are identical to the usual scalar NEM operators $\Pi_K$, $\nabla_K$ and $\mathcal{M}_K$ of [21], i.e., denoting $\mathbf{U}_K^i = (u_{\sigma}^i)_{\sigma \in \mathcal{F}_K}$ and $\mathbf{U}^i = (u_{\sigma}^i)_{\sigma \in \mathcal{F}}$:

$$ \Pi_K^d(\mathbf{U}_K)^i = \Pi_K(\mathbf{U}_K^i) \quad \text{and} \quad \nabla_K^d(\mathbf{U}_K)^i = \nabla_K(\mathbf{U}_K^i) \quad \text{and} \quad \mathcal{M}_K^d(\mathbf{U}_K)^i = \mathcal{M}_K(\mathbf{U}_K^i). $$

Also notice that $\nabla \Pi_K^d(\mathbf{U}_K) = \nabla_K^d(\mathbf{U}_K)$ by construction, thus we will use both notations interchangeably. Finally, we denote:

$$ \boldsymbol{\epsilon}_K(\mathbf{U}_K) = \frac{1}{2} (\nabla_K^d(\mathbf{U}_K) + \nabla_K^d(\mathbf{U}_K)^T), \qquad (13) $$

which is of course intended to be an approximation of $\boldsymbol{\epsilon}(\mathbf{u})$.

### 3.3. Approximation by the network element method

The discretization by the network element method consists in finding a solution $\mathbf{U} \in X_{\mathcal{N},0}^d$ of

$$ a_h(\mathbf{U}, \mathbf{V}) = l_h(\mathbf{V}) \quad \text{for all } \mathbf{V} \in X_{\mathcal{N},0}^d, \qquad (14) $$

where for the right-hand side assuming that $f_K$ is an approximation of $f$ at $\mathbf{x}_K$ we define the linear form $l_h : X_{\mathcal{N}}^d \longrightarrow \mathbb{R}$ by setting:

$$ l_h(\mathbf{V}) = \sum_{K \in \mathcal{T}} m_K f_K \cdot \mathcal{M}_K^d(\mathbf{V}_K), \qquad (15) $$

while

$$ a_h(\mathbf{U}, \mathbf{V}) = \sum_{K \in \mathcal{T}} a_h^K(\mathbf{U}_K, \mathbf{V}_K) \quad \text{with} \quad a_h^K(\mathbf{U}_K, \mathbf{V}_K) = a_{\mu,h}^K(\mathbf{U}_K, \mathbf{V}_K) + a_{\lambda,h}^K(\mathbf{U}_K, \mathbf{V}_K), \qquad (16) $$

where the discrete bilinear form $a_{\mu,h}^K : X_{\mathcal{N},K}^d \times X_{\mathcal{N},K}^d \longrightarrow \mathbb{R}$ mimics the term $2\mu \int_{\Omega} \boldsymbol{\epsilon}(\mathbf{u}) : \boldsymbol{\epsilon}(\mathbf{v})$ and is given by

$$ a_{\mu,h}^K(\mathbf{U}_K, \mathbf{V}_K) = 2\mu m_K \boldsymbol{\epsilon}_K(\mathbf{U}_K) : \boldsymbol{\epsilon}_K(\mathbf{V}_K) + s^K(\mathbf{U}_K, \mathbf{V}_K), \qquad (17) $$
where $s^K$ is a positive symmetric bilinear form on $X_{\mathcal{N},K}^d \times X_{\mathcal{N},K'}^d$, such that

$$s^K(\mathbf{U}_K, \mathbf{V}_K) = 2\mu m_K h_K^{-2} \sum_{\sigma \in \mathcal{F}_K} \sum_{\sigma' \in \mathcal{F}_K} S_K^{\sigma, \sigma'} (\mathbf{u}_\sigma - \Pi_K^d(\mathbf{U}_K)(\mathbf{x}_\sigma)) \cdot (\mathbf{v}_{\sigma'} - \Pi_K^d(\mathbf{V}_K)(\mathbf{x}_{\sigma'})), \quad (18)$$

where $\mathcal{S}_K = (S_K^{\sigma, \sigma'})_{\sigma, \sigma' \in \mathcal{F}_K}$ can be any symmetric positive definite matrix independent of the geometry $\mathcal{G}$ associated to the network, such that there exists $S_* > 0$ and $S^* > 0$ for which, for any $K \in \mathcal{F}$ and any $(\xi_\sigma)_{\sigma \in \mathcal{F}_K} \in \mathbb{R}^{\text{card}(\mathcal{F}_K)}$:

$$S_* \sum_{\sigma \in \mathcal{F}_K} |\xi_\sigma|^2 \le \sum_{\sigma \in \mathcal{F}_K} \sum_{\sigma' \in \mathcal{F}_K} S_K^{\sigma, \sigma'} \xi_\sigma \xi_{\sigma'} \le S^* \sum_{\sigma \in \mathcal{F}_K} |\xi_\sigma|^2.$$

The second discrete bilinear form $a_{\lambda,h}^K : X_{\mathcal{N},K}^d \times X_{\mathcal{N},K'}^d \to \mathbb{R}$ this time mimics $\lambda \int_{\Omega} \text{div } \mathbf{u} \text{ div } \mathbf{v}$ and is defined by

$$a_{\lambda,h}^{K}(\mathbf{U}_{K}, \mathbf{V}_{K}) = \lambda m_{K} \mathcal{D} \mathcal{I} \mathcal{V}_{K}(\mathbf{U}_{K}) \mathcal{D} \mathcal{I} \mathcal{V}_{K}(\mathbf{V}_{K}). \quad (19)$$

### 3.4. A discrete Korn inequality

At the continuous level existence and uniqueness of a solution of the variational problem (9) is a consequence of the celebrated Korn inequality, valid for any $\mathbf{u} \in H^1(\Omega)^d$ with $\Omega$ at least Lipschitz (see for instance [44]):

$$\|\nabla \mathbf{u}\|_{L^2(\Omega)^{d \times d}}^2 \le C_{\mathcal{K}, \Omega} \left( \|\mathbf{u}\|_{L^2(\Omega)^d}^2 + \|\boldsymbol{\epsilon}(\mathbf{u})\|_{L^2(\Omega)^{d \times d}}^2 \right), \quad (20)$$

where $C_{\mathcal{K},\Omega}$ only depends on $d$ and $\Omega$. To derive a stable numerical method, we need a discrete counterpart of Korn's inequality relating $\boldsymbol{\epsilon}_K(\mathbf{U}_K)$ and $\nabla_K^d(\mathbf{U}_K)$. It is well known that this requires to correctly manage the set of rigid body motions:

$$RM(\mathbb{R}^d) = \{ \mathbf{v} \in \mathbb{P}^1(\mathbb{R}^d)^d \mid \nabla \mathbf{v} = -\nabla \mathbf{v}^T \}. \quad (21)$$

To this end, we need to add the following assumption on the network:

(H5) For any $(K, L) \in \mathcal{T}^2$, if $\mathcal{F}_K \cap \mathcal{F}_L \ne \emptyset$ then the set $(\mathbf{x}_\sigma)_{\sigma \in \mathcal{F}_K \cap \mathcal{F}_L}$ is unisolvent for $RM(\mathbb{R}^d)$. Moreover, the set $(\mathbf{x}_\sigma)_{\sigma \in \mathcal{F}_{\text{ext}}}$ is unisolvent for $RM(\mathbb{R}^d)$.

As we will ultimately deal with a NEM equivalent of piecewise $H^1$ functions, similarly to what is done for meshes (see [45]) our discrete norm must moreover take into account jumps on cell interfaces for a discrete Korn inequality to hold. This leads to the following result:

**Lemma 1 (Discrete Korn's inequality).** Assume that (H5) holds. Then, there exists $C_{\mathcal{K},\mathcal{N}} > 0$ independent of $h$ such that for any $\mathbf{U} \in X_{\mathcal{N}}^d$ we have:

$$\sum_{K \in \mathcal{T}} m_K |\nabla_K^d (\mathbf{U}_K)|^2 \le C_{\mathcal{K},\mathcal{N}} \sum_{K \in \mathcal{T}} m_K \left( |\boldsymbol{\epsilon}_K(\mathbf{U}_K)|^2 + h_K^{-2} \sum_{\sigma \in \mathcal{F}_K} |\mathbf{u}_\sigma - \Pi_K^d(\mathbf{U}_K)(\mathbf{x}_\sigma)|^2 + h_K^{-2} \sum_{\sigma \in \mathcal{F}_K \cap \mathcal{F}_{\text{ext}}} |\mathbf{u}_\sigma|^2 \right). \quad (22)$$

**Proof.** We proceed by contradiction. If the result does not hold true, then for any $n \in \mathbb{N}$ there exists $\mathbf{U}^n \in X_{\mathcal{N}}^d$ such that:

$$\sum_{K \in \mathcal{T}} m_K |\nabla_K^d (\mathbf{U}_K^n)|^2 > n \sum_{K \in \mathcal{T}} m_K \left( |\boldsymbol{\epsilon}_K (\mathbf{U}_K^n)|^2 + h_K^{-2} \sum_{\sigma \in \mathcal{F}_K} |\mathbf{u}_\sigma^n - \Pi_K^d (\mathbf{U}_K^n)(\mathbf{x}_\sigma)|^2 + h_K^{-2} \sum_{\sigma \in \mathcal{F}_K \cap \mathcal{F}_{\text{ext}}} |\mathbf{u}_\sigma^n|^2 \right).$$

We define $\tilde{\mathbf{U}}^n$ by setting

$$\tilde{u}_{\sigma}^{i,n} = \frac{u_{\sigma}^{i,n}}{\sum_{K \in \mathcal{T}} m_K |\nabla_K^d (\mathbf{U}_K^n)|^2} \quad \text{for all } \sigma \in \mathcal{F},$$
where $s^K$ is a positive symmetric bilinear form on $X_{\mathcal{N},K}^d \times X_{\mathcal{N},K}^d$, such that

$$s^K(\mathbf{U}_K, \mathbf{V}_K) = 2\mu m_K h_K^{-2} \sum_{\sigma \in \mathcal{F}_K} \sum_{\sigma' \in \mathcal{F}_K} S_K^{\sigma, \sigma'} (\mathbf{u}_\sigma - \Pi_K^d(\mathbf{U}_K)(\mathbf{x}_\sigma)) \cdot (\mathbf{v}_{\sigma'} - \Pi_K^d(\mathbf{V}_K)(\mathbf{x}_{\sigma'})), \quad (18)$$

where $\mathcal{S}_K = (S_K^{\sigma, \sigma'})_{\sigma, \sigma' \in \mathcal{F}_K}$ can be any symmetric positive definite matrix independent of the geometry $\mathcal{G}$ associated to the network, such that there exists $S_* > 0$ and $S^* > 0$ for which, for any $K \in \mathcal{F}$ and any $(\xi_\sigma)_{\sigma \in \mathcal{F}_K} \in \mathbb{R}^{\text{card}(\mathcal{F}_K)}$:

$$S_* \sum_{\sigma \in \mathcal{F}_K} |\xi_\sigma|^2 \le \sum_{\sigma \in \mathcal{F}_K} \sum_{\sigma' \in \mathcal{F}_K} S_K^{\sigma, \sigma'} \xi_\sigma \xi_{\sigma'} \le S^* \sum_{\sigma \in \mathcal{F}_K} |\xi_\sigma|^2.$$

The second discrete bilinear form $a_{\lambda,h}^K : X_{\mathcal{N},K}^d \times X_{\mathcal{N},K}^d \to \mathbb{R}$ this time mimics $\lambda \int_{\Omega} \text{div} \, \mathbf{u} \text{ div} \, \mathbf{v}$ and is defined by

$$a_{\lambda,h}^{K}(\mathbf{U}_{K}, \mathbf{V}_{K}) = \lambda m_{K} \mathcal{D} \mathcal{I} \mathcal{V}_{K}(\mathbf{U}_{K}) \mathcal{D} \mathcal{I} \mathcal{V}_{K}(\mathbf{V}_{K}). \quad (19)$$

### 3.4. A discrete Korn inequality

At the continuous level existence and uniqueness of a solution of the variational problem (9) is a consequence of the celebrated Korn inequality, valid for any $\mathbf{u} \in H^1(\Omega)^d$ with $\Omega$ at least Lipschitz (see for instance [44]):

$$\|\nabla \mathbf{u}\|_{L^2(\Omega)^{d \times d}}^2 \le C_{\mathcal{K},\Omega} \left( \|\mathbf{u}\|_{L^2(\Omega)^d}^2 + \|\boldsymbol{\epsilon}(\mathbf{u})\|_{L^2(\Omega)^{d \times d}}^2 \right), \quad (20)$$

where $C_{\mathcal{K},\Omega}$ only depends on $d$ and $\Omega$. To derive a stable numerical method, we need a discrete counterpart of Korn's inequality relating $\boldsymbol{\epsilon}_K(\mathbf{U}_K)$ and $\nabla_K^d(\mathbf{U}_K)$. It is well known that this requires to correctly manage the set of rigid body motions:

$$RM(\mathbb{R}^d) = \{ \mathbf{v} \in \mathbb{P}^1(\mathbb{R}^d)^d \mid \nabla \mathbf{v} = -\nabla \mathbf{v}^T \}. \quad (21)$$

To this end, we need to add the following assumption on the network:

(H5) For any $(K, L) \in \mathcal{T}^2$, if $\mathcal{F}_K \cap \mathcal{F}_L \ne \emptyset$ then the set $(\mathbf{x}_\sigma)_{\sigma \in \mathcal{F}_K \cap \mathcal{F}_L}$ is unisolvent for $RM(\mathbb{R}^d)$. Moreover, the set $(\mathbf{x}_\sigma)_{\sigma \in \mathcal{F}_{\text{ext}}}$ is unisolvent for $RM(\mathbb{R}^d)$.

As we will ultimately deal with a NEM equivalent of piecewise $H^1$ functions, similarly to what is done for meshes (see [45]) our discrete norm must moreover take into account jumps on cell interfaces for a discrete Korn inequality to hold. This leads to the following result:

**Lemma 1 (Discrete Korn's inequality).** Assume that (H5) holds. Then, there exists $C_{\mathcal{K},\mathcal{N}} > 0$ independent of $h$ such that for any $\mathbf{U} \in X_{\mathcal{N}}^d$ we have:

$$\sum_{K \in \mathcal{T}} m_K |\nabla_K^d (\mathbf{U}_K)|^2 \le C_{\mathcal{K},\mathcal{N}} \sum_{K \in \mathcal{T}} m_K \left( |\boldsymbol{\epsilon}_K(\mathbf{U}_K)|^2 + h_K^{-2} \sum_{\sigma \in \mathcal{F}_K} |\mathbf{u}_\sigma - \Pi_K^d(\mathbf{U}_K)(\mathbf{x}_\sigma)|^2 + h_K^{-2} \sum_{\sigma \in \mathcal{F}_K \cap \mathcal{F}_{\text{ext}}} |\mathbf{u}_\sigma|^2 \right). \quad (22)$$

**Proof.** We proceed by contradiction. If the result does not hold true, then for any $n \in \mathbb{N}$ there exists $\mathbf{U}^n \in X_{\mathcal{N}}^d$ such that:

$$\sum_{K \in \mathcal{T}} m_K |\nabla_K^d (\mathbf{U}_K^n)|^2 > n \sum_{K \in \mathcal{T}} m_K \left( |\boldsymbol{\epsilon}_K (\mathbf{U}_K^n)|^2 + h_K^{-2} \sum_{\sigma \in \mathcal{F}_K} |\mathbf{u}_\sigma^n - \Pi_K^d (\mathbf{U}_K^n)(\mathbf{x}_\sigma)|^2 + h_K^{-2} \sum_{\sigma \in \mathcal{F}_K \cap \mathcal{F}_{\text{ext}}} |\mathbf{u}_\sigma^n|^2 \right).$$

We define $\tilde{\mathbf{U}}^n$ by setting

$$\tilde{u}_{\sigma}^{i,n} = \frac{u_{\sigma}^{i,n}}{\sum_{K \in \mathcal{T}} m_K |\nabla_K^d (\mathbf{U}_K^n)|^2} \quad \text{for all } \sigma \in \mathcal{F},$$
thus leading to $\sum_{K \in \mathcal{T}} m_K |\nabla_K^d (\tilde{\mathbf{U}}_K^n)|^2 = 1$ and

$$ \sum_{K \in \mathcal{T}} m_K \left( |\boldsymbol{\epsilon}_K (\tilde{\mathbf{U}}_K^n)|^2 + h_K^{-2} \sum_{\sigma \in \mathcal{F}_K} |\tilde{\mathbf{u}}_\sigma^n - \Pi_K^d (\tilde{\mathbf{U}}_K^n) (\mathbf{x}_\sigma)|^2 + h_K^{-2} \sum_{\sigma \in \mathcal{F}_K \cap \mathcal{F}_{\text{ext}}} |\tilde{\mathbf{u}}_\sigma^n|^2 \right) < \frac{1}{n}. $$

The set $X_{\mathcal{N}}^d$ being finite dimensional, it is compact and thus we can extract a strongly convergent sub-sequence from $(\tilde{\mathbf{U}}_K^n)_{n \in \mathbb{N}}$, still denoted $(\tilde{\mathbf{U}}_K^n)$, with limit $\tilde{\mathbf{U}}$. Then, by construction:

$$ \sum_{K \in \mathcal{T}} m_K \left( |\boldsymbol{\epsilon}_K (\tilde{\mathbf{U}}_K)|^2 + h_K^{-2} \sum_{\sigma \in \mathcal{F}_K} |\tilde{\mathbf{u}}_\sigma - \Pi_K^d (\tilde{\mathbf{U}}_K) (\mathbf{x}_\sigma)|^2 + h_K^{-2} \sum_{\sigma \in \mathcal{F}_K \cap \mathcal{F}_{\text{ext}}} |\tilde{\mathbf{u}}_\sigma|^2 \right) = 0, $$

which implies that $\boldsymbol{\epsilon}_K(\tilde{\mathbf{U}}_K) = \mathbf{0}$ for all $K \in \mathcal{T}$ and thus $\Pi_K^d(\tilde{\mathbf{U}}_K)$ is a rigid body motion for all $K \in \mathcal{T}$. Moreover $\tilde{\mathbf{u}}_\sigma = \Pi_K^d(\tilde{\mathbf{U}}_K)(\mathbf{x}_\sigma)$ for all $K \in \mathcal{T}$ and all $\sigma \in \mathcal{F}_K$, which using (H4) and (H5) implies that all the $\Pi_K^d(\tilde{\mathbf{U}}_K)$ are equal to the same rigid body motion $\Pi_K^d(\tilde{\mathbf{U}}_K) = \mathbf{a} + R_v \cdot \mathbf{x}$ with $R_v$ a anti-symmetric matrix. Finally, as $\tilde{\mathbf{u}}_\sigma = \mathbf{0}$ for all $\sigma \in \mathcal{F}_{\text{ext}}$ and $\mathbf{a} + R_v \cdot \mathbf{x}_\sigma = \tilde{\mathbf{u}}_\sigma$, we get that $\Pi_K^d(\tilde{\mathbf{U}}_K) = \mathbf{0}$ for all $K \in \mathcal{T}$ and $\tilde{\mathbf{u}}_\sigma = \mathbf{0}$ for all $\sigma \in \mathcal{T}$. Thus, in particular $\nabla_K^d(\tilde{\mathbf{U}}_K) = \mathbf{0}$ for all $K \in \mathcal{T}$, which contradicts $\sum_{K \in \mathcal{T}} m_K |\nabla_K^d (\tilde{\mathbf{U}}_K)|^2 = 1$.

To conclude the proof, we use a scaling argument. For any $\beta > 0$, consider the set $\hat{\Omega}(\beta) = \beta^{-1}\Omega$ obtained through the variable change $\hat{\mathbf{x}} = \beta^{-1}\mathbf{x}$, and the associated network $\mathcal{N}(\beta)$ such that $\hat{\mathbf{x}}_K = \beta^{-1}\mathbf{x}_K$, $\hat{h}_K = \beta^{-1}h_K$ and $\hat{m}_K = \beta^{-d}m_K$ for all $K \in \mathcal{T}$, $\hat{\mathbf{x}}_\sigma = \beta^{-1}\mathbf{x}_\sigma$ for all $\sigma \in \mathcal{T}$, and finally $\hat{\eta}_{K,\sigma} = \beta^{1-d}\eta_{K,\sigma}$ for all $K \in \mathcal{T}$ and all $\sigma \in \mathcal{F}_K$. Notice that:

$$ \hat{\nabla}_K^d (\mathbf{U}_K)_{ij} = \frac{1}{\hat{m}_K} \sum_{\sigma \in \mathcal{F}_K} u_\sigma^i \hat{\eta}_{K,\sigma}^j = \frac{\beta^{1-d}}{\beta^{-d} m_K} \sum_{\sigma \in \mathcal{F}_K} u_\sigma^i \eta_{K,\sigma}^j = \beta \nabla_K^d (\mathbf{U}_K)_{ij}, $$

leading to $|\hat{\nabla}_K^d (\mathbf{U}_K)|^2 = \beta^{-2}|\nabla_K^d (\mathbf{U}_K)|^2$, $|\hat{\boldsymbol{\epsilon}}_K (\mathbf{U}_K)|^2 = \beta^{-2}|\boldsymbol{\epsilon}_K (\mathbf{U}_K)|^2$ and $\hat{\Pi}_K^d (\mathbf{U}_K)(\hat{\mathbf{x}}_\sigma) = \Pi_K^d (\mathbf{U}_K)(\mathbf{x}_\sigma)$. Now, on $\hat{\Omega}(\beta)$ and $\mathcal{N}(\beta)$, we deduce from the first part that:

$$ \sum_{K \in \mathcal{T}} \hat{m}_K |\hat{\nabla}_K^d (\mathbf{U}_K)|^2 \\ \le C_{\mathcal{K}, \mathcal{N}(\beta)} \sum_{K \in \mathcal{T}} \hat{m}_K \left( |\hat{\boldsymbol{\epsilon}}_K (\mathbf{U}_K)|^2 + \hat{h}_K^{-2} \sum_{\sigma \in \mathcal{F}_K} |\mathbf{u}_\sigma - \hat{\Pi}_K^d (\mathbf{U}_K) (\hat{\mathbf{x}}_\sigma)|^2 + \hat{h}_K^{-2} \sum_{\sigma \in \mathcal{F}_K \cap \mathcal{F}_{\text{ext}}} |\mathbf{u}_\sigma|^2 \right), $$

leading to, gathering the above scaling results:

$$ \sum_{K \in \mathcal{T}} \beta^{-d} m_K \beta^2 |\nabla_K^d (\mathbf{U}_K)|^2 \\ \le C_{\mathcal{K}, \mathcal{N}(\beta)} \sum_{K \in \mathcal{T}} \beta^{-d} m_K \left( \beta^2 |\boldsymbol{\epsilon}_K (\mathbf{U}_K)|^2 + \beta^2 \hat{h}_K^{-2} \sum_{\sigma \in \mathcal{F}_K} |\mathbf{u}_\sigma - \Pi_K^d (\mathbf{U}_K) (\mathbf{x}_\sigma)|^2 + \beta^2 \hat{h}_K^{-2} \sum_{\sigma \in \mathcal{F}_K \cap \mathcal{F}_{\text{ext}}} |\mathbf{u}_\sigma|^2 \right), $$

and thus $C_{\mathcal{K}, \mathcal{N}} \le C_{\mathcal{K}, \mathcal{N}(\beta)}$ for all $\beta$ (and in particular $\beta = h$), which concludes the proof. $\square$

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

For our model problem with Dirichlet boundary conditions, we obtain the following obvious corollary:

**Corollary 1 (Discrete Korn's first inequality).** Assume that (H5) holds. Then there exists $C_{\mathcal{K}, \mathcal{N}} > 0$ independent of $h$ such that for any $\mathbf{U} \in X_{\mathcal{N},0}^d$, we have:

$$ \sum_{K \in \mathcal{T}} m_K |\nabla_K^d (\mathbf{U}_K)|^2 \le C_{\mathcal{K}, \mathcal{N}} \sum_{K \in \mathcal{T}} m_K \left( |\boldsymbol{\epsilon}_K (\mathbf{U}_K)|^2 + h_K^{-2} \sum_{\sigma \in \mathcal{F}_K} |\mathbf{u}_\sigma - \Pi_K^d (\mathbf{U}_K) (\mathbf{x}_\sigma)|^2 \right). $$
thus leading to $\sum_{K \in \mathcal{T}} m_K |\nabla_K^d (\tilde{U}_K^n)|^2 = 1$ and

$$ \sum_{K \in \mathcal{T}} m_K \left( |\boldsymbol{\epsilon}_K (\tilde{\boldsymbol{U}}_K^n)|^2 + h_K^{-2} \sum_{\sigma \in \mathcal{F}_K} |\tilde{\boldsymbol{u}}_\sigma^n - \Pi_K^d (\tilde{\boldsymbol{U}}_K^n) (\boldsymbol{x}_\sigma)|^2 + h_K^{-2} \sum_{\sigma \in \mathcal{F}_K \cap \mathcal{F}_{\text{ext}}} |\tilde{\boldsymbol{u}}_\sigma^n|^2 \right) < \frac{1}{n}. $$

The set $X_{\mathcal{N}}^d$ being finite dimensional, it is compact and thus we can extract a strongly convergent sub-sequence from $(\tilde{\boldsymbol{U}}^n)_{n \in \mathbb{N}}$, still denoted $(\tilde{\boldsymbol{U}}^n)$, with limit $\tilde{\boldsymbol{U}}$. Then, by construction:

$$ \sum_{K \in \mathcal{T}} m_K \left( |\boldsymbol{\epsilon}_K (\tilde{\boldsymbol{U}}_K)|^2 + h_K^{-2} \sum_{\sigma \in \mathcal{F}_K} |\tilde{\boldsymbol{u}}_\sigma - \Pi_K^d (\tilde{\boldsymbol{U}}_K) (\boldsymbol{x}_\sigma)|^2 + h_K^{-2} \sum_{\sigma \in \mathcal{F}_K \cap \mathcal{F}_{\text{ext}}} |\tilde{\boldsymbol{u}}_\sigma|^2 \right) = 0, $$

which implies that $\boldsymbol{\epsilon}_K(\tilde{\boldsymbol{U}}_K) = \mathbf{0}$ for all $K \in \mathcal{T}$ and thus $\Pi_K^d(\tilde{\boldsymbol{U}}_K)$ is a rigid body motion for all $K \in \mathcal{T}$. Moreover $\tilde{\boldsymbol{u}}_\sigma = \Pi_K^d(\tilde{\boldsymbol{U}}_K)(\boldsymbol{x}_\sigma)$ for all $K \in \mathcal{T}$ and all $\sigma \in \mathcal{F}_K$, which using (H4) and (H5) implies that all the $\Pi_K^d(\tilde{\boldsymbol{U}}_K)$ are equal to the same rigid body motion $\Pi_K^d(\tilde{\boldsymbol{U}}_K) = \mathbf{a} + R_v \cdot \mathbf{x}$ with $R_v$ a anti-symmetric matrix. Finally, as $\tilde{\boldsymbol{u}}_\sigma = \mathbf{0}$ for all $\sigma \in \mathcal{F}_{\text{ext}}$ and $\mathbf{a}+R_v \cdot \mathbf{x}_\sigma = \tilde{\boldsymbol{u}}_\sigma$, we get that $\Pi_K^d(\tilde{\boldsymbol{U}}_K) = \mathbf{0}$ for all $K \in \mathcal{T}$ and $\tilde{\boldsymbol{u}}_\sigma = \mathbf{0}$ for all $\sigma \in \mathcal{T}$. Thus, in particular $\nabla_K^d(\tilde{\boldsymbol{U}}_K) = \mathbf{0}$ for all $K \in \mathcal{T}$, which contradicts $\sum_{K \in \mathcal{T}} m_K |\nabla_K^d (\tilde{\boldsymbol{U}}_K)|^2 = 1$.

To conclude the proof, we use a scaling argument. For any $\beta > 0$, consider the set $\hat{\Omega}(\beta) = \beta^{-1}\Omega$ obtained through the variable change $\hat{\mathbf{x}} = \beta^{-1}\mathbf{x}$, and the associated network $\mathcal{N}(\beta)$ such that $\hat{\mathbf{x}}_K = \beta^{-1}\mathbf{x}_K$, $\hat{h}_K = \beta^{-1}h_K$ and $\hat{m}_K = \beta^{-d}m_K$ for all $K \in \mathcal{T}$, $\hat{\mathbf{x}}_\sigma = \beta^{-1}\mathbf{x}_\sigma$ for all $\sigma \in \mathcal{T}$, and finally $\hat{\eta}_{K,\sigma} = \beta^{1-d}\eta_{K,\sigma}$ for all $K \in \mathcal{T}$ and all $\sigma \in \mathcal{F}_K$. Notice that:

$$ \hat{\nabla}_K^d (\mathbf{U}_K)_{ij} = \frac{1}{\hat{m}_K} \sum_{\sigma \in \mathcal{F}_K} u_\sigma^i \hat{\eta}_{K,\sigma}^j = \frac{\beta^{1-d}}{\beta^{-d} m_K} \sum_{\sigma \in \mathcal{F}_K} u_\sigma^i \eta_{K,\sigma}^j = \beta \nabla_K^d (\mathbf{U}_K)_{ij}, $$

leading to $|\hat{\nabla}_K^d (\mathbf{U}_K)|^2 = \beta^{-2}|\nabla_K^d (\mathbf{U}_K)|^2$, $|\hat{\boldsymbol{\epsilon}}_K (\mathbf{U}_K)|^2 = \beta^{-2}|\boldsymbol{\epsilon}_K (\mathbf{U}_K)|^2$ and $\hat{\Pi}_K^d (\mathbf{U}_K)(\hat{\mathbf{x}}_\sigma) = \Pi_K^d (\mathbf{U}_K)(\mathbf{x}_\sigma)$. Now, on $\hat{\Omega}(\beta)$ and $\mathcal{N}(\beta)$, we deduce from the first part that:

$$ \sum_{K \in \mathcal{T}} \hat{m}_K |\hat{\nabla}_K^d (\mathbf{U}_K)|^2 \\ \le C_{\mathcal{K}, \mathcal{N}(\beta)} \sum_{K \in \mathcal{T}} \hat{m}_K \left( |\hat{\boldsymbol{\epsilon}}_K (\mathbf{U}_K)|^2 + \hat{h}_K^{-2} \sum_{\sigma \in \mathcal{F}_K} |\mathbf{u}_\sigma - \hat{\Pi}_K^d (\mathbf{U}_K) (\hat{\mathbf{x}}_\sigma)|^2 + \hat{h}_K^{-2} \sum_{\sigma \in \mathcal{F}_K \cap \mathcal{F}_{\text{ext}}} |\mathbf{u}_\sigma|^2 \right), $$

leading to, gathering the above scaling results:

$$ \sum_{K \in \mathcal{T}} \beta^{-d} m_K \beta^2 |\nabla_K^d (\mathbf{U}_K)|^2 \\ \le C_{\mathcal{K}, \mathcal{N}(\beta)} \sum_{K \in \mathcal{T}} \beta^{-d} m_K \left( \beta^2 |\boldsymbol{\epsilon}_K (\mathbf{U}_K)|^2 + \beta^2 \hat{h}_K^{-2} \sum_{\sigma \in \mathcal{F}_K} |\mathbf{u}_\sigma - \Pi_K^d (\mathbf{U}_K) (\mathbf{x}_\sigma)|^2 + \beta^2 \hat{h}_K^{-2} \sum_{\sigma \in \mathcal{F}_K \cap \mathcal{F}_{\text{ext}}} |\mathbf{u}_\sigma|^2 \right), $$

and thus $C_{\mathcal{K}, \mathcal{N}} \le C_{\mathcal{K}, \mathcal{N}(\beta)}$ for all $\beta$ (and in particular $\beta = h$), which concludes the proof. $\square$

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

For our model problem with Dirichlet boundary conditions, we obtain the following obvious corollary:

**Corollary 1 (Discrete Korn's first inequality).** Assume that (H5) holds. Then there exists $C_{\mathcal{K}, \mathcal{N}} > 0$ independent of $h$ such that for any $\mathbf{U} \in X_{\mathcal{N},0}^d$, we have:

$$ \sum_{K \in \mathcal{T}} m_K |\nabla_K^d (\mathbf{U}_K)|^2 \le C_{\mathcal{K}, \mathcal{N}} \sum_{K \in \mathcal{T}} m_K \left( |\boldsymbol{\epsilon}_K (\mathbf{U}_K)|^2 + h_K^{-2} \sum_{\sigma \in \mathcal{F}_K} |\mathbf{u}_\sigma - \Pi_K^d (\mathbf{U}_K) (\mathbf{x}_\sigma)|^2 \right). $$
### 3.5. Stability properties of the network element method for linear elasticity

Let us recall the quality parameters of the discretization network and geometry (see [21, 30]):

$$ \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$. We endow the space of degrees of freedom $X_{\mathcal{N}}^d$ with the bilinear forms

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

and the associated norm $\|\mathbf{U}\|_0^2 = (\mathbf{U}, \mathbf{U})_0$ and semi-norm $|\mathbf{U}|_1^2 = (\mathbf{U}, \mathbf{U})_1$. where

$$ (\mathbf{U}, \mathbf{V})_{0,K} = m_K \mathcal{M}_K^d (\mathbf{U}_K) \cdot \mathcal{M}_K^d (\mathbf{V}_K) + \sum_{\sigma \in \mathcal{F}_K} m_K (\mathbf{u}_\sigma - \mathcal{M}_K^d (\mathbf{U}_K)) \cdot (\mathbf{v}_\sigma - \mathcal{M}_K^d (\mathbf{V}_K)), $$

and

$$ (\mathbf{U}, \mathbf{V})_{1,K} = \sum_{\sigma \in \mathcal{F}_K} m_K h_K^{-2} (\mathbf{u}_\sigma - \mathcal{M}_K^d (\mathbf{U}_K)) \cdot (\mathbf{v}_\sigma - \mathcal{M}_K^d (\mathbf{V}_K)). $$

Next, we define:

$$ (\mathbf{U}, \mathbf{V})_X = (\mathbf{U}, \mathbf{V})_0 + (\mathbf{U}, \mathbf{V})_1 \quad \text{and} \quad \|\mathbf{U}\|_X^2 = (\mathbf{U}, \mathbf{U})_X, $$

which are obviously a scalar product and its associated norm on $X_{\mathcal{N}}^d$, making $X_{\mathcal{N}}^d$ a Hilbert space. From [21, 30], we know that there exists some $C_{\nabla}(\theta_{\mathcal{F}}, \theta_{\Pi}, \theta_{\mathcal{A}})^2 > 0$ (simply denoted $C_{\nabla}$ in the remaining of the paper for simplicity) such that for any $\mathbf{U} \in X_{\mathcal{N}}$:

$$ \left( \sum_{K \in \mathcal{T}} m_K |\nabla_K^d (\mathbf{U}_K)|^2 \right)^{1/2} \le C_{\nabla} (\theta_{\mathcal{F}}, \theta_{\Pi}, \theta_{\mathcal{A}}) \|\mathbf{U}\|_X, \quad (23) $$

immediately implying that:

$$ \left( \sum_{K \in \mathcal{T}} m_K |\epsilon_K (\mathbf{U}_K)|^2 \right)^{1/2} \le C_{\nabla} (\theta_{\mathcal{F}}, \theta_{\Pi}, \theta_{\mathcal{A}}) \|\mathbf{U}\|_X. \quad (24) $$

The following result establishes the stability of the NEM for linear elasticity:

**Theorem 1.** Assume that (H5) holds. There exists $\gamma_a > 0$ depending only on $C_{\mathcal{K},\mathcal{N}}$, $S_*$, $\mu$, $\theta_{\mathcal{F}}$ such that for any $\mathbf{U} \in X_{\mathcal{N},0}^d$

$$ a_h(\mathbf{U}, \mathbf{U}) \ge \gamma_a \sum_{K \in \mathcal{T}} \sum_{\sigma \in \mathcal{F}_K} m_K h_K^{-2} |\mathbf{u}_\sigma - \mathcal{M}_K^d(\mathbf{U}_K)|^2. \quad (25) $$

Moreover, there exists $C_a > 0$ depending only on $S_*$, $\mu$, $\lambda$, $\theta_{\mathcal{F}}$, $\theta_{\Pi}$ and $\theta_{\mathcal{A}}$ such that for any $(\mathbf{U}, \mathbf{V}) \in X_{\mathcal{N}}^d \times X_{\mathcal{N}}^d$

$$ a_h(\mathbf{U}, \mathbf{V}) \le C_a \left( \sum_{K \in \mathcal{T}} \sum_{\sigma \in \mathcal{F}_K} m_K h_K^{-2} |\mathbf{u}_\sigma - \mathcal{M}_K^d(\mathbf{U}_K)|^2 \right)^{1/2} \left( \sum_{K \in \mathcal{T}} \sum_{\sigma \in \mathcal{F}_K} m_K h_K^{-2} |\mathbf{v}_\sigma - \mathcal{M}_K^d(\mathbf{V}_K)|^2 \right)^{1/2} . \quad (26) $$

**Proof.** By construction, we have:

$$ a_h^K(\mathbf{U}, \mathbf{U}) \ge a_{\mu,h}^K(\mathbf{U}, \mathbf{U}) = 2\mu m_K |\epsilon_K(\mathbf{U}_K)|^2 + 2\mu m_K h_K^{-2} S_* \sum_{\sigma \in \mathcal{F}_K} |\mathbf{u}_\sigma - \Pi_K^d(\mathbf{U}_K)(\mathbf{x}_\sigma)|^2. $$
### 3.5. Stability properties of the network element method for linear elasticity

Let us recall the quality parameters of the discretization network and geometry (see [21, 30]):

$$ \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{F}} = \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$. We endow the space of degrees of freedom $X_{\mathcal{N}}^d$ with the bilinear forms

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

and the associated norm $\|\mathbf{U}\|_0^2 = (\mathbf{U}, \mathbf{U})_0$ and semi-norm $|\mathbf{U}|_1^2 = (\mathbf{U}, \mathbf{U})_1$. where

$$ (\mathbf{U}, \mathbf{V})_{0,K} = m_K \mathcal{M}_K^d(\mathbf{U}_K) \cdot \mathcal{M}_K^d(\mathbf{V}_K) + \sum_{\sigma \in \mathcal{F}_K} m_K (\mathbf{u}_\sigma - \mathcal{M}_K^d(\mathbf{U}_K)) \cdot (\mathbf{v}_\sigma - \mathcal{M}_K^d(\mathbf{V}_K)), $$

and

$$ (\mathbf{U}, \mathbf{V})_{1,K} = \sum_{\sigma \in \mathcal{F}_K} m_K h_K^{-2} (\mathbf{u}_\sigma - \mathcal{M}_K^d(\mathbf{U}_K)) \cdot (\mathbf{v}_\sigma - \mathcal{M}_K^d(\mathbf{V}_K)). $$

Next, we define:

$$ (\mathbf{U}, \mathbf{V})_X = (\mathbf{U}, \mathbf{V})_0 + (\mathbf{U}, \mathbf{V})_1 \quad \text{and} \quad \|\mathbf{U}\|_X^2 = (\mathbf{U}, \mathbf{U})_X, $$

which are obviously a scalar product and its associated norm on $X_{\mathcal{N}}^d$, making $X_{\mathcal{N}}^d$ a Hilbert space. From [21, 30], we know that there exists some $C_{\nabla}(\theta_{\mathcal{F}}, \theta_{\Pi}, \theta_{\mathcal{A}})^2 > 0$ (simply denoted $C_{\nabla}$ in the remaining of the paper for simplicity) such that for any $\mathbf{U} \in X_{\mathcal{N}}$:

$$ \left( \sum_{K \in \mathcal{T}} m_K |\nabla_K^d (\mathbf{U}_K)|^2 \right)^{1/2} \le C_{\nabla} (\theta_{\mathcal{F}}, \theta_{\Pi}, \theta_{\mathcal{A}}) \|\mathbf{U}\|_X, \quad (23) $$

immediately implying that:

$$ \left( \sum_{K \in \mathcal{T}} m_K |\epsilon_K (\mathbf{U}_K)|^2 \right)^{1/2} \le C_{\nabla} (\theta_{\mathcal{F}}, \theta_{\Pi}, \theta_{\mathcal{A}}) \|\mathbf{U}\|_X. \quad (24) $$

The following result establishes the stability of the NEM for linear elasticity:

**Theorem 1.** Assume that (H5) holds. There exists $\gamma_a > 0$ depending only on $C_{\mathcal{K},\mathcal{N}}$, $S_*$, $\mu$, $\theta_{\mathcal{F}}$ such that for any $\mathbf{U} \in X_{\mathcal{N},0}^d$

$$ a_h(\mathbf{U}, \mathbf{U}) \ge \gamma_a \sum_{K \in \mathcal{T}} \sum_{\sigma \in \mathcal{F}_K} m_K h_K^{-2} |\mathbf{u}_\sigma - \mathcal{M}_K^d(\mathbf{U}_K)|^2. \quad (25) $$

Moreover, there exists $C_a > 0$ depending only on $S_*$, $\mu$, $\lambda$, $\theta_{\mathcal{F}}$, $\theta_{\Pi}$ and $\theta_{\mathcal{A}}$ such that for any $(\mathbf{U}, \mathbf{V}) \in X_{\mathcal{N}}^d \times X_{\mathcal{N}}^d$

$$ a_h(\mathbf{U}, \mathbf{V}) \le C_a \left( \sum_{K \in \mathcal{T}} \sum_{\sigma \in \mathcal{F}_K} m_K h_K^{-2} |\mathbf{u}_\sigma - \mathcal{M}_K^d(\mathbf{U}_K)|^2 \right)^{1/2} \left( \sum_{K \in \mathcal{T}} \sum_{\sigma \in \mathcal{F}_K} m_K h_K^{-2} |\mathbf{v}_\sigma - \mathcal{M}_K^d(\mathbf{V}_K)|^2 \right)^{1/2} . \quad (26) $$

**Proof.** By construction, we have:

$$ a_h^K(\mathbf{U}, \mathbf{U}) \ge a_{\mu,h}^K(\mathbf{U}, \mathbf{U}) = 2\mu m_K |\epsilon_K(\mathbf{U}_K)|^2 + 2\mu m_K h_K^{-2} S_* \sum_{\sigma \in \mathcal{F}_K} |\mathbf{u}_\sigma - \Pi_K^d(\mathbf{U}_K)(\mathbf{x}_\sigma)|^2. $$
Using the discrete Korn inequality of corollary 1, we see that:

$$ \sum_{K \in \mathcal{T}} \mu m_K |\boldsymbol{\epsilon}_K(\mathbf{U}_K)|^2 + \mu m_K h_K^{-2} S_* \sum_{K \in \mathcal{T}} \sum_{\sigma \in \mathcal{F}_K} |\boldsymbol{u}_\sigma - \Pi_K^d(\mathbf{U}_K)(\boldsymbol{x}_\sigma)|^2 \\ \geq \frac{\mu}{C_{\mathcal{K},\mathcal{N}}} \min(1, S_*) \sum_{K \in \mathcal{T}} m_K |\nabla_K^d(\mathbf{U})|^2, $$

leading to:

$$ a_h(\mathbf{U}, \mathbf{U}) = \sum_{K \in \mathcal{T}} a_h^K(\mathbf{U}, \mathbf{U}) \geq \sum_{K \in \mathcal{T}} \mu m_K |\boldsymbol{\epsilon}_K(\mathbf{U}_K)|^2 + \mu S_* \sum_{K \in \mathcal{T}} \sum_{\sigma \in \mathcal{F}_K} m_K h_K^{-2} |\boldsymbol{u}_\sigma - \Pi_K^d(\mathbf{U}_K)(\boldsymbol{x}_\sigma)|^2 \\ + \frac{\mu}{C_{\mathcal{K},\mathcal{N}}} \min(1, S_*) \sum_{K \in \mathcal{T}} m_K |\nabla_K^d(\mathbf{U})|^2. $$

Using the identity

$$ (a-b)^2 \geq \frac{\rho}{1+\rho} a^2 - \rho b^2 \quad \forall (a,b) \in \mathbb{R}^2, \forall \rho > -1, $$

we see that, for all $\rho > -1$

$$ |\boldsymbol{u}_\sigma - \Pi_K^d(\mathbf{U}_K)(\boldsymbol{x}_\sigma)|^2 \geq \frac{\rho}{1+\rho} \sum_{\sigma \in \mathcal{F}_K} |\boldsymbol{u}_\sigma - \mathcal{M}_K^d(\mathbf{U}_K)|^2 - \rho \theta_{\mathcal{F}} |\nabla_K^d(\mathbf{U}_K)|^2. $$

Gathering the previous results, we get:

$$ a_h(\mathbf{U}, \mathbf{U}) \geq \sum_{K \in \mathcal{T}} \mu m_K |\boldsymbol{\epsilon}_K(\mathbf{U}_K)|^2 + \frac{\mu \rho S_*}{1+\rho} \sum_{\sigma \in \mathcal{F}_K} m_K h_K^{-2} |\boldsymbol{u}_\sigma - \mathcal{M}_K^d(\mathbf{U}_K)|^2 \\ + \left( \frac{\min(1, S_*)}{C_{\mathcal{K},\mathcal{N}}} - \rho \theta_{\mathcal{F}} S_* \right) \mu m_K |\nabla_K^d(\mathbf{U})|^2, $$

and the first result follows with:

$$ \rho = \frac{\min(1, S_*)}{C_{\mathcal{K},\mathcal{N}} \theta_{\mathcal{F}} S_*} \quad \text{and} \quad \gamma_a = \frac{\mu S_* \min(1, S_*)}{C_{\mathcal{K},\mathcal{N}} \theta_{\mathcal{F}} S_* + \min(1, S_*)}. $$

The second result is an immediate consequence of estimates (23) and (24) and Cauchy-Schwarz inequality (see [21, 28, 30] for details in the scalar case). □

If $\Omega$ satisfies the cone condition then using (H4) and applying the results of [21, 30] componentwise and summing, we know that a discrete Poincaré inequality holds for the $\|\cdot\|_X$ and $\|\cdot\|_0$ norms. Existence, uniqueness and stability of the discrete solution is consequently an obvious consequence of Lax-Milgram's lemma.

## 3.6. Matrix formulation

By definition, we have:

$$ m_K \boldsymbol{\epsilon}_K(\mathbf{U}_K) : \boldsymbol{\epsilon}_K(\mathbf{V}_K) = \frac{1}{4m_K} \sum_{\sigma \in \mathcal{F}_K} \sum_{\sigma' \in \mathcal{F}_K} \sum_{i=1}^{d} \sum_{j=1}^{d} (u_{\sigma'}^i \eta_{K,\sigma'}^j + u_{\sigma'}^j \eta_{K,\sigma'}^i) (v_{\sigma}^i \eta_{K,\sigma}^j + v_{\sigma}^j \eta_{K,\sigma}^i) \\ = \frac{1}{2m_K} \sum_{\sigma \in \mathcal{F}_K} \sum_{\sigma' \in \mathcal{F}_K} \sum_{i=1}^{d} u_{\sigma'}^i v_{\sigma}^i \left( \sum_{j=1}^{d} \eta_{K,\sigma'}^j \eta_{K,\sigma}^j \right) \\ + \frac{1}{2m_K} \sum_{\sigma \in \mathcal{F}_K} \sum_{\sigma' \in \mathcal{F}_K} \sum_{i=1}^{d} \sum_{j=1}^{d} u_{\sigma'}^j v_{\sigma}^i \eta_{K,\sigma'}^i \eta_{K,\sigma}^j. $$

Thus, we get:

$$ a_{\mu,h}^{K}(\mathbf{U}_K, \mathbf{V}_K) = \sum_{\sigma \in \mathcal{F}_K} \sum_{i=1}^{d} \left( \sum_{\sigma' \in \mathcal{F}_K} \sum_{j=1}^{d} A_{\sigma,i,\sigma',j}^{K,\mu} u_{\sigma'}^{j} \right) v_{\sigma}^{i}, $$
Using the discrete Korn inequality of corollary 1, we see that:

$$ 
\begin{aligned}
& \sum_{K \in \mathcal{T}} \mu m_K |\epsilon_K(\mathbf{U}_K)|^2 + \mu m_K h_K^{-2} S_* \sum_{K \in \mathcal{T}} \sum_{\sigma \in \mathcal{F}_K} |\mathbf{u}_\sigma - \Pi_K^d(\mathbf{U}_K)(\mathbf{x}_\sigma)|^2 \\
& \geq \frac{\mu}{C_{\mathcal{K},\mathcal{N}}} \min(1, S_*) \sum_{K \in \mathcal{T}} m_K |\nabla_K^d(\mathbf{U})|^2,
\end{aligned}
$$

leading to:

$$ a_h(\mathbf{U}, \mathbf{U}) = \sum_{K \in \mathcal{T}} a_h^K(\mathbf{U}, \mathbf{U}) \geq \sum_{K \in \mathcal{T}} \mu m_K |\epsilon_K(\mathbf{U}_K)|^2 + \mu S_* \sum_{K \in \mathcal{T}} \sum_{\sigma \in \mathcal{F}_K} m_K h_K^{-2} |\mathbf{u}_\sigma - \Pi_K^d(\mathbf{U}_K)(\mathbf{x}_\sigma)|^2 \\ + \frac{\mu}{C_{\mathcal{K},\mathcal{N}}} \min(1, S_*) \sum_{K \in \mathcal{T}} m_K |\nabla_K^d(\mathbf{U})|^2. $$

Using the identity

$$ (a-b)^2 \geq \frac{\rho}{1+\rho} a^2 - \rho b^2 \quad \forall (a,b) \in \mathbb{R}^2, \forall \rho > -1, $$

we see that, for all $\rho > -1$

$$ |\mathbf{u}_\sigma - \Pi_K^d(\mathbf{U}_K)(\mathbf{x}_\sigma)|^2 \geq \frac{\rho}{1+\rho} \sum_{\sigma \in \mathcal{F}_K} |\mathbf{u}_\sigma - \mathcal{M}_K^d(\mathbf{U}_K)|^2 - \rho \theta_{\mathcal{F}} |\nabla_K^d(\mathbf{U}_K)|^2. $$

Gathering the previous results, we get:

$$ a_h(\mathbf{U}, \mathbf{U}) \geq \sum_{K \in \mathcal{T}} \mu m_K |\epsilon_K(\mathbf{U}_K)|^2 + \frac{\mu \rho S_*}{1+\rho} \sum_{\sigma \in \mathcal{F}_K} m_K h_K^{-2} |\mathbf{u}_\sigma - \mathcal{M}_K^d(\mathbf{U}_K)|^2 \\ + \left( \frac{\min(1, S_*)}{C_{\mathcal{K},\mathcal{N}}} - \rho \theta_{\mathcal{F}} S_* \right) \mu m_K |\nabla_K^d(\mathbf{U})|^2, $$

and the first result follows with:

$$ \rho = \frac{\min(1, S_*)}{C_{\mathcal{K},\mathcal{N}} \theta_{\mathcal{F}} S_*} \quad \text{and} \quad \gamma_a = \frac{\mu S_* \min(1, S_*)}{C_{\mathcal{K},\mathcal{N}} \theta_{\mathcal{F}} S_* + \min(1, S_*)}. $$

The second result is an immediate consequence of estimates (23) and (24) and Cauchy-Schwarz inequality (see [21, 28, 30] for details in the scalar case). □

If $\Omega$ satisfies the cone condition then using (H4) and applying the results of [21, 30] componentwise and summing, we know that a discrete Poincaré inequality holds for the $\|\cdot\|_X$ and $\|\cdot\|_0$ norms. Existence, uniqueness and stability of the discrete solution is consequently an obvious consequence of Lax-Milgram's lemma.

## 3.6. Matrix formulation

By definition, we have:

$$ 
\begin{aligned}
m_K \epsilon_K(\mathbf{U}_K) : \epsilon_K(\mathbf{V}_K) &= \frac{1}{4m_K} \sum_{\sigma \in \mathcal{F}_K} \sum_{\sigma' \in \mathcal{F}_K} \sum_{i=1}^d \sum_{j=1}^d (u_{\sigma'}^i \eta_{K,\sigma'}^j + u_{\sigma'}^j \eta_{K,\sigma'}^i) (v_{\sigma}^i \eta_{K,\sigma}^j + v_{\sigma}^j \eta_{K,\sigma}^i) \\
&= \frac{1}{2m_K} \sum_{\sigma \in \mathcal{F}_K} \sum_{\sigma' \in \mathcal{F}_K} \sum_{i=1}^d u_{\sigma'}^i v_{\sigma}^i \left( \sum_{j=1}^d \eta_{K,\sigma'}^j \eta_{K,\sigma}^j \right) \\
&\quad + \frac{1}{2m_K} \sum_{\sigma \in \mathcal{F}_K} \sum_{\sigma' \in \mathcal{F}_K} \sum_{i=1}^d \sum_{j=1}^d u_{\sigma'}^j v_{\sigma}^i \eta_{K,\sigma'}^i \eta_{K,\sigma}^j.
\end{aligned}
$$

Thus, we get:

$$ a_{\mu,h}^{K}(\mathbf{U}_K, \mathbf{V}_K) = \sum_{\sigma \in \mathcal{F}_K} \sum_{i=1}^{d} \left( \sum_{\sigma' \in \mathcal{F}_K} \sum_{j=1}^{d} A_{\sigma,i,\sigma',j}^{K,\mu} u_{\sigma'}^{j} \right) v_{\sigma}^{i}, $$
where

$$A_{\sigma,i,\sigma',j}^{K,\mu} = \frac{\mu}{m_K} \left( \eta_{K,\sigma'}^{i} \eta_{K,\sigma}^{j} + \delta_{ij} \left( \sum_{k=1}^{d} \eta_{K,\sigma}^{k} \eta_{K,\sigma'}^{k} \right) \right).$$

Meanwhile, we have:

$$m_K \mathcal{D} \mathcal{I} \mathcal{V}_K (\mathbf{U}_K) \mathcal{D} \mathcal{I} \mathcal{V}_K (\mathbf{V}_K) = \frac{1}{m_K} \sum_{\sigma \in \mathcal{F}_K} \sum_{\sigma' \in \mathcal{F}_K} \left( \sum_{j=1}^{d} u_{\sigma'}^{j} \eta_{K,\sigma'}^{j} \right) \left( \sum_{i=1}^{d} v_{\sigma}^{i} \eta_{K,\sigma}^{i} \right),$$

leading to:

$$a_{\lambda,h}^{K}(\mathbf{U}_{K}, \mathbf{V}_{K}) = \sum_{\sigma \in \mathcal{F}_{K}} \sum_{i=1}^{d} \left( \sum_{\sigma' \in \mathcal{F}_{K}} \sum_{j=1}^{d} A_{\sigma,i,\sigma',j}^{K,\lambda} u_{\sigma'}^{j} \right) v_{\sigma}^{i} \quad \text{where} \quad A_{\sigma,i,\sigma',j}^{K,\lambda} = \frac{\lambda}{m_{K}} \eta_{K,\sigma}^{i} \eta_{K,\sigma'}^{j}.$$

Finally, since:

$$\mathbf{u}_{\sigma} - \Pi_{K}^{d}(\mathbf{U}_{K})(\mathbf{x}_{\sigma}) = \sum_{\sigma'' \in \mathcal{F}_{K}} \left( \delta_{\sigma\sigma''} - \gamma_{K,\sigma''} - \frac{1}{m_{K}} \eta_{K,\sigma''} \cdot (\mathbf{x}_{\sigma} - \mathbf{x}_{K}) \right) u_{\sigma''} = \sum_{\sigma'' \in \mathcal{F}_{K}} y_{K,\sigma''}^{\sigma} u_{\sigma''},$$

and

$$\mathbf{v}_{\sigma'} - \Pi_{K}^{d}(\mathbf{V}_{K})(\mathbf{x}_{\sigma'}) = \sum_{\sigma''' \in \mathcal{F}_{K}} \left( \delta_{\sigma'\sigma'''} - \gamma_{K,\sigma'''} - \frac{1}{m_{K}} \eta_{K,\sigma''} \cdot (\mathbf{x}_{\sigma'} - \mathbf{x}_{K}) \right) v_{\sigma'''} = \sum_{\sigma''' \in \mathcal{F}_{K}} y_{K,\sigma''}^{\sigma'} v_{\sigma'''}.$$

where we have denoted

$$y_{K,\sigma''}^{\sigma} = \delta_{\sigma\sigma''} - \gamma_{K,\sigma''} - \frac{1}{m_K} \eta_{K,\sigma''} \cdot (\mathbf{x}_{\sigma} - \mathbf{x}_K) \quad \text{and} \quad \mathbf{y}_{K,\sigma''} = (y_{K,\sigma''}^{\sigma})_{\sigma \in \mathcal{F}_K},$$

we see that:

$$\sum_{\sigma \in \mathcal{F}_K} \sum_{\sigma' \in \mathcal{F}_K} S_K^{\sigma, \sigma'} (\mathbf{u}_\sigma - \Pi_K^d (\mathbf{U}_K) (\mathbf{x}_\sigma)) \cdot (\mathbf{v}_{\sigma'} - \Pi_K^d (\mathbf{V}_K) (\mathbf{x}_{\sigma'})) \\ = \sum_{i=1}^d \sum_{\sigma \in \mathcal{F}_K} \sum_{\sigma' \in \mathcal{F}_K} \sum_{\sigma'' \in \mathcal{F}_K} \sum_{\sigma''' \in \mathcal{F}_K} S_K^{\sigma, \sigma'} y_{K,\sigma''}^{\sigma} y_{K,\sigma'''}^{\sigma'} u_{\sigma''}^{i} v_{\sigma'''}^{i},$$

leading to (notice that we have inverted the role of the pairs of indices $(\sigma, \sigma')$ and $(\sigma'', \sigma''')$):

$$s^K(\mathbf{U}_K, \mathbf{V}_K) = \sum_{\sigma \in \mathcal{F}_K} \sum_{i=1}^{d} \left( \sum_{\sigma' \in \mathcal{F}_K} \sum_{j=1}^{d} S_{\sigma,i,\sigma',j}^{K,\mu} u_{\sigma'}^{j} \right) v_{\sigma}^{i},$$

where we have denoted:

$$S_{\sigma,i,\sigma',j}^{K,\mu} = 2\mu\delta_{ij}m_K h_K^{-2} \mathbf{y}_{K,\sigma}^{T} \mathcal{S}_K \mathbf{y}_{K,\sigma'} \quad \text{where} \quad \mathcal{S}_K = (S_K^{\sigma'',\sigma'''})^{\sigma'',\sigma'''} \in \mathcal{F}_K \quad \text{and} \quad \mathbf{y}_{K,\sigma} = (y_{K,\sigma}^{\sigma''})^{\sigma'',\sigma'''} \in \mathcal{F}_K.$$

Thus, introducing the matrix $\mathcal{A}_K = (A_{\sigma,i,\sigma',j}^K)_{1 \le i,j \le d, (\sigma,\sigma') \in \mathcal{F}_K \times \mathcal{F}_K}$ defined by

$$A_{\sigma,i,\sigma',j}^{K,\lambda} = A_{\sigma,i,\sigma',j}^{K,\mu} + A_{\sigma,i,\sigma',j}^{K,\lambda} + S_{\sigma,i,\sigma',j}^{K,\mu}$$

and rewriting the right hand side as:

$$l_h(\mathbf{V}) = \sum_{K \in \mathcal{F}} \sum_{i=1}^{d} \sum_{\sigma \in \mathcal{F}_K} \mathbb{L}_{\sigma,i}^{K} v_{\sigma,i} \quad \text{with} \quad \mathbb{L}_{\sigma,i}^{K} = m_K f_K^{i} \gamma_{K,\sigma},$$

we see that (14) is equivalent to solving:

$$\sum_{K \in \mathcal{F}} \sum_{i=1}^{d} \sum_{j=1}^{d} \sum_{\sigma \in \mathcal{F}_K} \sum_{\sigma' \in \mathcal{F}_K} \mathcal{A}_{\sigma,i,\sigma',j}^{K} u_{\sigma'}^{j} v_{\sigma}^{i} = \sum_{K \in \mathcal{F}} \sum_{i=1}^{d} \sum_{\sigma \in \mathcal{F}_K} \mathbb{L}_{\sigma,i}^{K} v_{\sigma}^{i}, \quad (27)$$

complemented by the boundary conditions $u_{\sigma}^{i} = 0$ and $v_{\sigma}^{i} = 0$ for any $\sigma \in \mathcal{F}_{\text{ext}}$ and any $1 \le i \le d$. Following [21] to get an even more familiar version of this system, let us introduce the square
where

$$A_{\sigma,i,\sigma',j}^{K,\mu} = \frac{\mu}{m_K} \left( \eta_{K,\sigma'}^{i} \eta_{K,\sigma}^{j} + \delta_{ij} \left( \sum_{k=1}^{d} \eta_{K,\sigma}^{k} \eta_{K,\sigma'}^{k} \right) \right).$$

Meanwhile, we have:

$$m_K \mathcal{D} \mathcal{I} \mathcal{V}_K (\mathbf{U}_K) \mathcal{D} \mathcal{I} \mathcal{V}_K (\mathbf{V}_K) = \frac{1}{m_K} \sum_{\sigma \in \mathcal{F}_K} \sum_{\sigma' \in \mathcal{F}_K} \left( \sum_{j=1}^{d} u_{\sigma'}^{j} \eta_{K,\sigma'}^{j} \right) \left( \sum_{i=1}^{d} v_{\sigma}^{i} \eta_{K,\sigma}^{i} \right),$$

leading to:

$$a_{\lambda,h}^{K}(\mathbf{U}_{K}, \mathbf{V}_{K}) = \sum_{\sigma \in \mathcal{F}_{K}} \sum_{i=1}^{d} \left( \sum_{\sigma' \in \mathcal{F}_{K}} \sum_{j=1}^{d} A_{\sigma,i,\sigma',j}^{K,\lambda} u_{\sigma'}^{j} \right) v_{\sigma}^{i} \quad \text{where} \quad A_{\sigma,i,\sigma',j}^{K,\lambda} = \frac{\lambda}{m_{K}} \eta_{K,\sigma}^{i} \eta_{K,\sigma'}^{j}.$$

Finally, since:

$$\mathbf{u}_{\sigma} - \Pi_{K}^{d}(\mathbf{U}_{K})(\mathbf{x}_{\sigma}) = \sum_{\sigma'' \in \mathcal{F}_{K}} \left( \delta_{\sigma\sigma''} - \gamma_{K,\sigma''} - \frac{1}{m_{K}} \eta_{K,\sigma''} \cdot (\mathbf{x}_{\sigma} - \mathbf{x}_{K}) \right) u_{\sigma''} = \sum_{\sigma'' \in \mathcal{F}_{K}} y_{K,\sigma''}^{\sigma} u_{\sigma''},$$

and

$$\mathbf{v}_{\sigma'} - \Pi_{K}^{d}(\mathbf{V}_{K})(\mathbf{x}_{\sigma'}) = \sum_{\sigma''' \in \mathcal{F}_{K}} \left( \delta_{\sigma'\sigma'''} - \gamma_{K,\sigma'''} - \frac{1}{m_{K}} \eta_{K,\sigma'''} \cdot (\mathbf{x}_{\sigma'} - \mathbf{x}_{K}) \right) v_{\sigma'''} = \sum_{\sigma''' \in \mathcal{F}_{K}} y_{K,\sigma'''}^{\sigma'} v_{\sigma'''},$$

where we have denoted

$$y_{K,\sigma''}^{\sigma} = \delta_{\sigma\sigma''} - \gamma_{K,\sigma''} - \frac{1}{m_K} \eta_{K,\sigma''} \cdot (\mathbf{x}_{\sigma} - \mathbf{x}_K) \quad \text{and} \quad \mathbf{y}_{K,\sigma''} = (y_{K,\sigma''}^{\sigma})_{\sigma \in \mathcal{F}_K},$$

we see that:

$$\sum_{\sigma \in \mathcal{F}_K} \sum_{\sigma' \in \mathcal{F}_K} S_K^{\sigma, \sigma'} (\mathbf{u}_\sigma - \Pi_K^d (\mathbf{U}_K) (\mathbf{x}_\sigma)) \cdot (\mathbf{v}_{\sigma'} - \Pi_K^d (\mathbf{V}_K) (\mathbf{x}_{\sigma'})) \\ = \sum_{i=1}^d \sum_{\sigma \in \mathcal{F}_K} \sum_{\sigma' \in \mathcal{F}_K} \sum_{\sigma'' \in \mathcal{F}_K} \sum_{\sigma''' \in \mathcal{F}_K} S_K^{\sigma, \sigma'} y_{K,\sigma''}^{\sigma} y_{K,\sigma'''}^{\sigma'} u_{\sigma''}^{i} v_{\sigma'''}^{i},$$

leading to (notice that we have inverted the role of the pairs of indices $(\sigma, \sigma')$ and $(\sigma'', \sigma''')$):

$$s^K(\mathbf{U}_K, \mathbf{V}_K) = \sum_{\sigma \in \mathcal{F}_K} \sum_{i=1}^{d} \left( \sum_{\sigma' \in \mathcal{F}_K} \sum_{j=1}^{d} S_{\sigma,i,\sigma',j}^{K,\mu} u_{\sigma'}^{j} \right) v_{\sigma}^{i},$$

where we have denoted:

$$S_{\sigma,i,\sigma',j}^{K,\mu} = 2\mu\delta_{ij}m_K h_K^{-2} \mathbf{y}_{K,\sigma}^{T} \mathcal{S}_K \mathbf{y}_{K,\sigma'} \text{ where } \mathcal{S}_K = (S_K^{\sigma'',\sigma'''})^{\sigma'',\sigma'''} \in \mathcal{F}_K \text{ and } \mathbf{y}_{K,\sigma} = (y_{K,\sigma}^{\sigma''})^{\sigma'',\sigma'''} \in \mathcal{F}_K.$$

Thus, introducing the matrix $\mathcal{A}_K = (A_{\sigma,i,\sigma',j}^K)_{1 \le i,j \le d, (\sigma,\sigma') \in \mathcal{F}_K \times \mathcal{F}_K}$ defined by

$$A_{\sigma,i,\sigma',j}^{K,\lambda} = A_{\sigma,i,\sigma',j}^{K,\mu} + A_{\sigma,i,\sigma',j}^{K,\mu} + S_{\sigma,i,\sigma',j}^{K,\mu}$$

and rewriting the right hand side as:

$$l_h(\mathbf{V}) = \sum_{K \in \mathcal{F}} \sum_{i=1}^{d} \sum_{\sigma \in \mathcal{F}_K} \mathbb{L}_{\sigma,i}^{K} v_{\sigma,i} \quad \text{with} \quad \mathbb{L}_{\sigma,i}^{K} = m_K f_K^{i} \gamma_{K,\sigma},$$

we see that (14) is equivalent to solving:

$$\sum_{K \in \mathcal{F}} \sum_{i=1}^{d} \sum_{j=1}^{d} \sum_{\sigma \in \mathcal{F}_K} \sum_{\sigma' \in \mathcal{F}_K} \mathcal{A}_{\sigma,i,\sigma',j}^{K} u_{\sigma'}^{j} v_{\sigma}^{i} = \sum_{K \in \mathcal{F}} \sum_{i=1}^{d} \sum_{\sigma \in \mathcal{F}_K} \mathbb{L}_{\sigma,i}^{K} v_{\sigma}^{i}, \quad (27)$$

complemented by the boundary conditions $u_{\sigma}^{i} = 0$ and $v_{\sigma}^{i} = 0$ for any $\sigma \in \mathcal{F}_{\text{ext}}$ and any $1 \le i \le d$. Following [21] to get an even more familiar version of this system, let us introduce the square
matrix $\mathcal{A}$ of size $d$ card($\mathcal{F}$) as well as $F \in X_N$ such that and for any $\sigma, \sigma' \in \mathcal{F}_{\text{int}} \times \mathcal{F}$ and any $1 \le i, j, \le d$:

$$ \mathcal{A}_{\sigma,i,\sigma',j} = \sum_{K \in \mathcal{F}_{\sigma} \cap \mathcal{F}_{\sigma'}} \mathcal{A}_{\sigma,i,\sigma',j}^{K} \quad \text{and} \quad F_{\sigma,i} = \sum_{K \in \mathcal{F}_{\sigma}} \mathbb{L}_{\sigma,i}^{K}, $$

and for any $\sigma, \sigma' \in \mathcal{F}_{\text{ext}} \times \mathcal{F}$ and any $1 \le i, j \le d$:

$$ \mathcal{A}_{\sigma,i,\sigma',j} = \delta_{\sigma',\sigma} \quad \text{and} \quad F_{\sigma,i} = 0. $$

Problem (14) is then equivalent to solving the linear system

$$ \mathcal{A} \mathbf{U} = \mathbf{F}. $$

Consequently both the assembly of the system or the treatment of the boundary degrees of freedom can be performed in a finite element fashion, looping over all cells $K \in \mathcal{T}$ and incrementally constructing the global system.

# 4. Flux based and mesh based versions

In this section, we explore variations around our network element method for linear elasticity. We propose a flux based version of our NEM and provide possible implementations for mesh based discretization networks. In particular, we explore mesh based networks with interfaces related to faces and show how to obtain a compact stencil method, to be compared with the method of [35].

## 4.1. A flux based version

Following [29], we can derive a flux based alternative formulation by enriching the set of degrees of freedom with cell unknowns:

$$ X_{\mathcal{N}}^{d} = \left\{ \left( (\mathbf{u}_{\mathcal{K}})_{K \in \mathcal{T}}, (u_{\sigma})_{\sigma \in \mathcal{F}} \right) \middle| \mathbf{u}_{\mathcal{K}} \in \mathbb{R}^{d} \quad \forall K \in \mathcal{T} \text{ and } \mathbf{u}_{\sigma} \in \mathbb{R}^{d} \quad \forall \sigma \in \mathcal{F} \right\}, $$

with of course

$$ X_{\mathcal{N},0}^{d} = \{ \mathbf{U} \in X_{\mathcal{N}}^{d} \mid \mathbf{u}_{\sigma} = 0 \text{ for all } \sigma \in \mathcal{F}_{\text{ext}} \}, $$

and

$$ X_{\mathcal{N},K}^{d} = \left\{ (\mathbf{u}_{\mathcal{K}}, (\mathbf{u}_{\sigma})_{\sigma \in \mathcal{F}_{\mathcal{K}}}) \mid \mathbf{u}_{\mathcal{K}} \in \mathbb{R}^{d} \text{ and } \mathbf{u}_{\sigma} \in \mathbb{R}^{d} \quad \forall \sigma \in \mathcal{F}_{\mathcal{K}} \right\}. $$

The modifications of the network element method consist in choosing $\mathcal{M}_{\mathcal{K}}^{d}(\mathbf{U}_{\mathcal{K}}) = \mathbf{u}_{\mathcal{K}}$ as well as:

$$ \nabla_{\mathcal{K}}^{d}(\mathbf{U}_{\mathcal{K}})_{ij} = \frac{1}{m_{\mathcal{K}}} \sum_{\sigma \in \mathcal{F}_{\mathcal{K}}} (u_{\sigma}^{i} - u_{\mathcal{K}}^{i}) \eta_{K,\sigma}^{j} \quad \text{and} \quad \mathcal{D} \mathcal{F} \mathcal{V}_{\mathcal{K}}(\mathbf{U}_{\mathcal{K}}) = \frac{1}{m_{\mathcal{K}}} \sum_{i=1}^{d} \sum_{\sigma \in \mathcal{F}_{\mathcal{K}}} (u_{\sigma}^{i} - u_{\mathcal{K}}^{i}) \eta_{K,\sigma}^{i}. $$

The remaining steps for defining the flux based formulation are then exactly the same as for (14), using (17)–(19) and (15). A careful look at the difference between the two formulations reveals that the flux formulation mainly consists in replacing any term $u_{\sigma}^{i}$ by $u_{\sigma}^{i} - u_{\mathcal{K}}^{i}$ in the discrete gradients $\nabla_{\mathcal{K}}^{d}$ and $\epsilon_{\mathcal{K}}$. A straightforward computation thus leads to:

$$ a_{h}^{K}(\mathbf{U}_{\mathcal{K}}, \mathbf{V}_{\mathcal{K}}) = \sum_{i=1}^{d} \sum_{\sigma \in \mathcal{F}_{\mathcal{K}}} \left( \sum_{j=1}^{d} \sum_{\sigma' \in \mathcal{F}_{\mathcal{K}}} A_{\sigma,i,\sigma',j}^{K} (u_{\sigma'},j - u_{\mathcal{K}}^{j}) \right) (v_{\sigma,i} - v_{\mathcal{K}}^{i}), $$

with the very same definition for $A_{\sigma,i,\sigma',j}^{K}$ as before except that we now have:

$$ y_{K,\sigma}^{\prime\prime} = \delta_{\sigma^{\prime\prime},\sigma} - \frac{1}{m_K} \eta_{K,\sigma} \cdot (\mathbf{x}_{\sigma^{\prime\prime}} - \mathbf{x}_K), $$
matrix $\mathcal{A}$ of size $d$ card($\mathcal{F}$) as well as $F \in X_N$ such that and for any $\sigma, \sigma' \in \mathcal{F}_{\text{int}} \times \mathcal{F}$ and any $1 \le i, j, \le d$:

$$ \mathcal{A}_{\sigma,i,\sigma',j} = \sum_{K \in \mathcal{F}_{\sigma} \cap \mathcal{F}_{\sigma'}} \mathcal{A}_{\sigma,i,\sigma',j}^{K} \quad \text{and} \quad F_{\sigma,i} = \sum_{K \in \mathcal{F}_{\sigma}} \mathbb{L}_{\sigma,i}^{K}, $$

and for any $\sigma, \sigma' \in \mathcal{F}_{\text{ext}} \times \mathcal{F}$ and any $1 \le i, j \le d$:

$$ \mathcal{A}_{\sigma,i,\sigma',j} = \delta_{\sigma',\sigma} \quad \text{and} \quad F_{\sigma,i} = 0. $$

Problem (14) is then equivalent to solving the linear system

$$ \mathcal{A}\mathbf{U} = \mathbf{F}. $$

Consequently both the assembly of the system or the treatment of the boundary degrees of freedom can be performed in a finite element fashion, looping over all cells $K \in \mathcal{T}$ and incrementally constructing the global system.

# 4. Flux based and mesh based versions

In this section, we explore variations around our network element method for linear elasticity. We propose a flux based version of our NEM and provide possible implementations for mesh based discretization networks. In particular, we explore mesh based networks with interfaces related to faces and show how to obtain a compact stencil method, to be compared with the method of [35].

## 4.1. A flux based version

Following [29], we can derive a flux based alternative formulation by enriching the set of degrees of freedom with cell unknowns:

$$ X_{\mathcal{N}}^{d} = \left\{ \left( (\mathbf{u}_{\mathcal{K}})_{K \in \mathcal{T}}, (u_{\sigma})_{\sigma \in \mathcal{F}} \right) \middle| \mathbf{u}_{\mathcal{K}} \in \mathbb{R}^{d} \quad \forall K \in \mathcal{T} \text{ and } \mathbf{u}_{\sigma} \in \mathbb{R}^{d} \quad \forall \sigma \in \mathcal{F} \right\}, $$

with of course

$$ X_{\mathcal{N},0}^{d} = \{ \mathbf{U} \in X_{\mathcal{N}}^{d} \mid \mathbf{u}_{\sigma} = 0 \text{ for all } \sigma \in \mathcal{F}_{\text{ext}} \}, $$

and

$$ X_{\mathcal{N},K}^{d} = \left\{ (\mathbf{u}_{\mathcal{K}}, (\mathbf{u}_{\sigma})_{\sigma \in \mathcal{F}_{\mathcal{K}}}) \mid \mathbf{u}_{\mathcal{K}} \in \mathbb{R}^{d} \text{ and } \mathbf{u}_{\sigma} \in \mathbb{R}^{d} \quad \forall \sigma \in \mathcal{F}_{\mathcal{K}} \right\}. $$

The modifications of the network element method consist in choosing $\mathcal{M}_{\mathcal{K}}^{d}(\mathbf{U}_{\mathcal{K}}) = \mathbf{u}_{\mathcal{K}}$ as well as:

$$ \nabla_{\mathcal{K}}^{d}(\mathbf{U}_{\mathcal{K}})_{ij} = \frac{1}{m_{\mathcal{K}}} \sum_{\sigma \in \mathcal{F}_{\mathcal{K}}} (u_{\sigma}^{i} - u_{\mathcal{K}}^{i}) \eta_{K,\sigma}^{j} \quad \text{and} \quad \mathcal{D} \mathcal{F} \mathcal{V}_{\mathcal{K}}(\mathbf{U}_{\mathcal{K}}) = \frac{1}{m_{\mathcal{K}}} \sum_{i=1}^{d} \sum_{\sigma \in \mathcal{F}_{\mathcal{K}}} (u_{\sigma}^{i} - u_{\mathcal{K}}^{i}) \eta_{K,\sigma}^{i}. $$

The remaining steps for defining the flux based formulation are then exactly the same as for (14), using (17)–(19) and (15). A careful look at the difference between the two formulations reveals that the flux formulation mainly consists in replacing any term $u_{\sigma}^{i}$ by $u_{\sigma}^{i} - u_{\mathcal{K}}^{i}$ in the discrete gradients $\nabla_{\mathcal{K}}^{d}$ and $\epsilon_{\mathcal{K}}$. A straightforward computation thus leads to:

$$ a_{h}^{K}(\mathbf{U}_{\mathcal{K}}, \mathbf{V}_{\mathcal{K}}) = \sum_{i=1}^{d} \sum_{\sigma \in \mathcal{F}_{\mathcal{K}}} \left( \sum_{j=1}^{d} \sum_{\sigma' \in \mathcal{F}_{\mathcal{K}}} A_{\sigma,i,\sigma',j}^{K} (u_{\sigma'}, j - u_{\mathcal{K}}^{j}) \right) (v_{\sigma,i} - v_{\mathcal{K}}^{i}), $$

with the very same definition for $A_{\sigma,i,\sigma',j}^{K}$ as before except that we now have:

$$ y_{K,\sigma}^{\prime\prime} = \delta_{\sigma^{\prime\prime}\sigma} - \frac{1}{m_K} \eta_{K,\sigma} \cdot (\mathbf{x}_{\sigma^{\prime\prime}} - \mathbf{x}_K), $$
as the coefficient $\gamma_{K,\sigma}$ disappears thanks to the new definition of $\mathcal{M}_K^d$. Naturally, we define flux functions by setting:

$$F_{K,\sigma,i}(\mathbf{U}_K) = \sum_{j=1}^{d} \sum_{\sigma' \in \mathcal{F}_K} A_{\sigma,i,\sigma',j}^{K}(u_K^{j} - u_{\sigma',j}), \quad (28)$$

and problem (14) can now be rewritten, taking one degree of freedom equal to one and all the others equal to zero:

$$\left| \begin{array}{ll} \sum_{\sigma \in \mathcal{F}_K} F_{K,\sigma,i}(\mathbf{U}_K) = f_K^i & \text{for all } K \in \mathcal{T} \text{ and all } 1 \le i \le d, \\ \sum_{K \in \mathcal{T}_{\sigma}} F_{K,\sigma,i}(\mathbf{U}_K) = 0 & \text{for all } \sigma \in \mathcal{F}_{\text{int}} \text{ and all } 1 \le i \le d, \\ u_{\sigma}^i = 0 & \text{for all } \sigma \in \mathcal{F}_{\text{ext}} \text{ and all } 1 \le i \le d. \end{array} \right. \quad (29)$$

## 4.2. Mesh based discretizations

As announced in the introduction, the proposed scheme is also interesting in the case of mesh based discretization, thanks to our direct analogy between meshes and networks. We consider two ways of identifying meshes with networks, using some subfaces of the faces or directly the vertices of the mesh.

### 4.2.1. Interfaces from faces

Given a mesh $\mathcal{M}$ of $\Omega$, as explained in Section 2 the most natural analogy between meshes and networks consists in identifying $\mathcal{T}$ with the set of cells of the mesh and $\mathcal{F}$ with the set of faces of the mesh, and moreover assuming that $\mathbf{x}_K$ is the barycenter of cell $K$ and $\mathbf{x}_\sigma$ is the barycenter of face $\sigma$. Doing so we identify the pair $K = \{\mathbf{x}_K, r_K\}$ of the network definition with the open set $K$ corresponding to a cell of the mesh $\mathcal{M}$, and the same holds for interfaces and faces. Then, an admissible discrete geometry is simply given by

$$m_K = |K| \quad \text{for any } K \in \mathcal{T} \quad \text{and} \quad \eta_{K,\sigma} = |\sigma|\mathbf{n}_{K,\sigma} \quad \text{for any } K \in \mathcal{T} \text{ and any } \sigma \in \mathcal{F}_K, \quad (30)$$

where $|K|$ and $|\sigma|$ denotes the Lebesgue measures of cell $K$ and face $\sigma$, and $\mathbf{n}_{K,\sigma}$ the unit normal to face $\sigma$, outgward for $K$. Using this analogy, the definition of the scheme is now the same for mesh based discretizations as for networks. Unfortunately, for many meshes hypothesis (H5) will not be satisfied and even if (H5) is only a sufficient condition it is well known that the discrete Korn inequality (22) will not hold true (see for instance [35]). A remedy for this is proposed in [35] and consists in adding for any interface $\sigma \in \mathcal{F}$ a jump term on the first order polynomial $(\Pi_K^d(\mathbf{U}_K))_{K \in \mathcal{T}}$, thus requiring integration of first order polynomials on $\sigma$. This cannot be straightforwardly extended in the network element context, which is the reason why here we propose another remedy that has the advantage of maintaining a compact stencil. Its main drawback being that the number of interfaces will increase.

Momentarily denoting with an index $\mathcal{M}$ the mesh sets for clarity, for any mesh face $\sigma_{\mathcal{M}} \in \mathcal{F}_{\mathcal{M},\text{int}}$, we introduce a set $\mathcal{S}\mathcal{F}_{\sigma_{\mathcal{M}}}$ of at least $d$ subfaces such that $\overline{\sigma_{\mathcal{M}}} = \cup_{f \in \mathcal{S}\mathcal{F}_{\sigma_{\mathcal{M}}}} \overline{f}$. In dimension 2 it simply consists in splitting each face from its barycenter into two sub-segments, while in dimension 3 if the face is star-shaped, which will generally be the case, one can also split it into the usual triangles by joining the face barycenter and its edges. Setting $\mathcal{S}\mathcal{F}_{\sigma_{\mathcal{M}}}$ = $\{\sigma_{\mathcal{M}}\}$ for any boundary face $\sigma_{\mathcal{M}} \in \mathcal{F}_{\mathcal{M},\text{ext}}$ and denoting $\mathcal{S}\mathcal{F}_{\mathcal{M}} = \cup_{\sigma_{\mathcal{M}} \in \mathcal{F}_{\mathcal{M}}} \mathcal{S}\mathcal{F}_{\sigma_{\mathcal{M}}}$, we now identify the set of interfaces $\mathcal{F}$ with the set of mesh subfaces $\mathcal{S}\mathcal{F}_{\mathcal{M}}$. Then, denoting $\mathcal{F}_K = \cup_{\sigma_{\mathcal{M}} \in \mathcal{F}_{\mathcal{M},K}} \mathcal{S}\mathcal{F}_{\sigma_{\mathcal{M}}}$, the set of subfaces of cell $K$, we can now define:

$$m_K = |K| \quad \text{for any } K \in \mathcal{T} \quad \text{and} \quad \eta_{K,\sigma} = |\sigma|\mathbf{n}_{K,\sigma} \quad \text{for any } K \in \mathcal{T} \text{ and any } \sigma \in \mathcal{F}_K, \quad (31)$$
as the coefficient $\gamma_{K,\sigma}$ disappears thanks to the new definition of $\mathcal{M}_K^d$. Naturally, we define flux functions by setting:

$$F_{K,\sigma,i}(\mathbf{U}_K) = \sum_{j=1}^{d} \sum_{\sigma' \in \mathcal{F}_K} A_{\sigma,i,\sigma',j}^{K}(u_K^{j} - u_{\sigma',j}), \quad (28)$$

and problem (14) can now be rewritten, taking one degree of freedom equal to one and all the others equal to zero:

$$\left| \begin{array}{ll} \sum_{\sigma \in \mathcal{F}_K} F_{K,\sigma,i}(\mathbf{U}_K) = f_K^i & \text{for all } K \in \mathcal{T} \text{ and all } 1 \le i \le d, \\ \sum_{K \in \mathcal{T}_{\sigma}} F_{K,\sigma,i}(\mathbf{U}_K) = 0 & \text{for all } \sigma \in \mathcal{F}_{\text{int}} \text{ and all } 1 \le i \le d, \\ u_{\sigma}^i = 0 & \text{for all } \sigma \in \mathcal{F}_{\text{ext}} \text{ and all } 1 \le i \le d. \end{array} \right. \quad (29)$$

## 4.2. Mesh based discretizations

As announced in the introduction, the proposed scheme is also interesting in the case of mesh based discretization, thanks to our direct analogy between meshes and networks. We consider two ways of identifying meshes with networks, using some subfaces of the faces or directly the vertices of the mesh.

### 4.2.1. Interfaces from faces

Given a mesh $\mathcal{M}$ of $\Omega$, as explained in Section 2 the most natural analogy between meshes and networks consists in identify $\mathcal{T}$ with the set of cells of the mesh and $\mathcal{F}$ with the set of faces of the mesh, and moreover assuming that $\mathbf{x}_K$ is the barycenter of cell $K$ and $\mathbf{x}_\sigma$ is the barycenter of face $\sigma$. Doing so we identify the pair $K = \{\mathbf{x}_K, r_K\}$ of the network definition with the open set $K$ corresponding to a cell of the mesh $\mathcal{M}$, and the same holds for interfaces and faces. Then, an admissible discrete geometry is simply given by

$$m_K = |K| \quad \text{for any } K \in \mathcal{T} \quad \text{and} \quad \eta_{K,\sigma} = |\sigma|\mathbf{n}_{K,\sigma} \quad \text{for any } K \in \mathcal{T} \text{ and any } \sigma \in \mathcal{F}_K, \quad (30)$$

where $|K|$ and $|\sigma|$ denotes the Lebesgue measures of cell $K$ and face $\sigma$, and $\mathbf{n}_{K,\sigma}$ the unit normal to face $\sigma$, outgward for $K$. Using this analogy, the definition of the scheme is now the same for mesh based discretizations as for networks. Unfortunately, for many meshes hypothesis (H5) will not be satisfied and even if (H5) is only a sufficient condition it is well known that the discrete Korn inequality (22) will not hold true (see for instance [35]). A remedy for this is proposed in [35] and consists in adding for any interface $\sigma \in \mathcal{F}$ a jump term on the first order polynomial $(\Pi_K^d(\mathbf{U}_K))_{K \in \mathcal{T}}$, thus requiring integration of first order polynomials on $\sigma$. This cannot be straightforwardly extended in the network element context, which is the reason why here we propose another remedy that has the advantage of maintaining a compact stencil. Its main drawback being that the number of interfaces will increase.

Momentarily denoting with an index $\mathcal{M}$ the mesh sets for clarity, for any mesh face $\sigma_{\mathcal{M}} \in \mathcal{F}_{\mathcal{M},\text{int}}$, we introduce a set $\mathcal{S}\mathcal{F}_{\sigma_{\mathcal{M}}}$ of at least $d$ subfaces such that $\overline{\sigma_{\mathcal{M}}} = \cup_{f \in \mathcal{S}\mathcal{F}_{\sigma_{\mathcal{M}}}} \overline{f}$. In dimension 2 it simply consists in splitting each face from its barycenter into two sub-segments, while in dimension 3 if the face is star-shaped, which will generally be the case, one can also split it into the usual triangles by joining the face barycenter and its edges. Setting $\mathcal{S}\mathcal{F}_{\sigma_{\mathcal{M}}}$ = $\{\sigma_{\mathcal{M}}\}$ for any boundary face $\sigma_{\mathcal{M}} \in \mathcal{F}_{\mathcal{M},\text{ext}}$ and denoting $\mathcal{S}\mathcal{F}_{\mathcal{M}} = \cup_{\sigma_{\mathcal{M}} \in \mathcal{F}_{\mathcal{M}}} \mathcal{S}\mathcal{F}_{\sigma_{\mathcal{M}}}$, we now identify the set of interfaces $\mathcal{F}$ with the set of mesh subfaces $\mathcal{S}\mathcal{F}_{\mathcal{M}}$. Then, denoting $\mathcal{F}_K = \cup_{\sigma_{\mathcal{M}} \in \mathcal{F}_{\mathcal{M},K}} \mathcal{S}\mathcal{F}_{\sigma_{\mathcal{M}}}$, the set of subfaces of cell $K$, we can now define:

$$m_K = |K| \quad \text{for any } K \in \mathcal{T} \quad \text{and} \quad \eta_{K,\sigma} = |\sigma|\mathbf{n}_{K,\sigma} \quad \text{for any } K \in \mathcal{T} \text{ and any } \sigma \in \mathcal{F}_K, \quad (31)$$
where $\sigma \in \mathcal{F}$ now coincides with a mesh subface $f \in \mathcal{F}\mathcal{F}_{\mathcal{M}}$, to obtain a mesh based discretization network that satisfies (H5). The discrete gradient and divergence operators then coincide with the classical hybrid gradient and divergence operators:

$$ \nabla_K^d (\mathbf{U}_K)_{ij} = \frac{1}{|K|} \sum_{\sigma \in \mathcal{F}_K} |\sigma| u_{\sigma}^i \mathbf{n}_{K,\sigma}^j \quad \text{and} \quad \mathcal{D}\mathcal{F}\mathcal{V}_K (\mathbf{U}_K) = \frac{1}{|K|} \sum_{i=1}^d \sum_{\sigma \in \mathcal{F}_K} |\sigma| u_{\sigma}^i \mathbf{n}_{K,\sigma}^i. $$

Compared to existing approaches for face based first order discretizations of linear elasticity problems on meshes, the main interest of our formulation lies in the stabilization term that only involves the interfaces of a cell, avoiding the need to enlarge the stencil of the method to recover stability contrary to the alternative first order approaches of the literature also based on the hybrid gradient of [35]. The drawback is that we have more unknowns in our final linear system than for the approach of [35]. Also notice that by adding buffer variables located at the quadrature points of the additional jump term of the formulation of [35], one could rewrite the scheme of [35] on a stencil as compact as ours, increasing the number of unknowns in a similar way than in our case.

Finally, it is obvious that one can also enrich the set of degrees of freedom with cell centered unknowns in the mesh based context, defining again $\mathcal{M}_K^d(\mathbf{U}_K) = \mathbf{u}_K$ and proceeding as in the previous section, to obtain a small stencil, unconditionally stable flux based formulation for mesh based discretizations (note that there is no real need to modify the discrete gradient and divergence in this case, as $\sum_{\sigma \in \mathcal{F}_K} |\sigma| \mathbf{n}_{K,\sigma} = 0$).

**Remark 1.** Notice that our results are not in contradiction with those of [35]: indeed, in [35] it was established that no discrete Korn's inequality for a discrete norm equivalent to $\|\cdot\|_X$ can exists if the interfaces are defined as mesh faces in general, as depending on the mesh configuration there might exists non zero piecewise rigid body motions with value zero for norm $\|\cdot\|_X$ and on the boundary. Here, as our interfaces based on subfaces satisfy (H5), this not possible: any piecewise rigid body motion that cancels the corresponding subface based norm $\|\cdot\|_X$ is in fact a single body motion over the whole domain with zero value on the boundary, and is thus zero.

### 4.2.2. Interfaces from vertices

As was also explained in Section 2, given a mesh $\mathcal{M}$ of $\Omega$ another very natural analogy between meshes and networks consists in again identifying $\mathcal{T}$ with the set of cells of the mesh but this time identifying $\mathcal{F}$ with the set of vertices of the mesh, and moreover assuming that $\mathbf{x}_K$ is the barycenter of cell $K$ and $\mathbf{x}_\sigma$ is position of vertex $\sigma$. In this case, the consistency conditions (2), (3) and conservation property (4) on the $\eta_{K,\sigma}$'s are exactly satisfied, i.e. without the error terms, by the coefficients $\mathbf{g}_{K,s}$ of the vertex averaged schemes (VAG) of [46]. Thus, setting:

$$ m_K = |K| \quad \text{for any } K \in \mathcal{T} $$

and

$$ \eta_{K,\sigma} = \mathbf{g}_{K,s} \quad \text{for any } K \in \mathcal{T} \text{ and } \sigma \in \mathcal{F}_K \text{ corresponding to vertex } s, \quad (32) $$

we obtain an admissible network geometry. Denoting $\mathcal{V}_K$ the set of vertices of cell $K$, the discrete gradient and divergence operators then coincide with the VAG gradient and divergence operators:

$$ \nabla_K^d (\mathbf{U}_K)_{ij} = \frac{1}{|K|} \sum_{s \in \mathcal{V}_K} u_s^i \mathbf{g}_{K,s}^j \quad \text{and} \quad \mathcal{D}\mathcal{F}\mathcal{V}_K (\mathbf{U}_K) = \frac{1}{|K|} \sum_{i=1}^d \sum_{s \in \mathcal{V}_K} u_s^i \mathbf{g}_{K,s}^i. $$

If the $\mathbf{g}_{K,s}$ are computed using virtual element techniques as in [47], then the resulting scheme coincides with the first order virtual element scheme for linear elasticity of [33]. Once again one can obtain a flux based formulation by adding cell centered unknowns, setting again $\mathcal{M}_K^d(\mathbf{U}_K) = \mathbf{u}_K$. In this case, if the $\mathbf{g}_{K,s}$ are computed using virtual element techniques as in [47] we
where $\sigma \in \mathcal{F}$ now coincides with a mesh subface $f \in \mathcal{S}\mathcal{F}\mathcal{M}$, to obtain a mesh based discretization network that satisfies (H5). The discrete gradient and divergence operators then coincide with the classical hybrid gradient and divergence operators:

$$ \nabla_K^d (\mathbf{U}_K)_{ij} = \frac{1}{|K|} \sum_{\sigma \in \mathcal{F}_K} |\sigma| u_{\sigma}^i \mathbf{n}_{K,\sigma}^j \quad \text{and} \quad \mathcal{D}\mathcal{F}\mathcal{V}_K (\mathbf{U}_K) = \frac{1}{|K|} \sum_{i=1}^d \sum_{\sigma \in \mathcal{F}_K} |\sigma| u_{\sigma}^i \mathbf{n}_{K,\sigma}^i. $$

Compared to existing approaches for face based first order discretizations of linear elasticity problems on meshes, the main interest of our formulation lies in the stabilization term that only involves the interfaces of a cell, avoiding the need to enlarge the stencil of the method to recover stability contrary to the alternative first order approaches of the literature also based on the hybrid gradient of [35]. The drawback is that we have more unknowns in our final linear system than for the approach of [35]. Also notice that by adding buffer variables located at the quadrature points of the additional jump term of the formulation of [35], one could rewrite the scheme of [35] on a stencil as compact as ours, increasing the number of unknowns in a similar way than in our case.

Finally, it is obvious that one can also enrich the set of degrees of freedom with cell centered unknowns in the mesh based context, defining again $\mathcal{M}_K^d(\mathbf{U}_K) = \mathbf{u}_K$ and proceeding as in the previous section, to obtain a small stencil, unconditionally stable flux based formulation for mesh based discretizations (note that there is no real need to modify the discrete gradient and divergence in this case, as $\sum_{\sigma \in \mathcal{F}_K} |\sigma| \mathbf{n}_{K,\sigma} = 0$).

**Remark 1.** Notice that our results are not in contradiction with those of [35]: indeed, in [35] it was established that no discrete Korn's inequality for a discrete norm equivalent to $\|\cdot\|_X$ can exists if the interfaces are defined as mesh faces in general, as depending on the mesh configuration there might exists non zero piecewise rigid body motions with value zero for norm $\|\cdot\|_X$ and on the boundary. Here, as our interfaces based on subfaces satisfy (H5), this not possible: any piecewise rigid body motion that cancels the corresponding subface based norm $\|\cdot\|_X$ is in fact a single body motion over the whole domain with zero value on the boundary, and is thus zero.

### 4.2.2. Interfaces from vertices

As was also explained in Section 2, given a mesh $\mathcal{M}$ of $\Omega$ another very natural analogy between meshes and networks consists in again identifying $\mathcal{T}$ with the set of cells of the mesh but this time identifying $\mathcal{F}$ with the set of vertices of the mesh, and moreover assuming that $\mathbf{x}_K$ is the barycenter of cell $K$ and $\mathbf{x}_\sigma$ is position of vertex $\sigma$. In this case, the consistency conditions (2), (3) and conservation property (4) on the $\eta_{K,\sigma}$'s are exactly satisfied, i.e. without the error terms, by the coefficients $\mathbf{g}_{K,s}$ of the vertex averaged schemes (VAG) of [46]. Thus, setting:

$$ m_K = |K| \quad \text{for any } K \in \mathcal{T} $$

and

$$ \eta_{K,\sigma} = \mathbf{g}_{K,s} \quad \text{for any } K \in \mathcal{T} \text{ and } \sigma \in \mathcal{F}_K \text{ corresponding to vertex } s, \quad (32) $$

we obtain an admissible network geometry. Denoting $\mathcal{V}_K$ the set of vertices of cell $K$, the discrete gradient and divergence operators then coincide with the VAG gradient and divergence operators:

$$ \nabla_K^d (\mathbf{U}_K)_{ij} = \frac{1}{|K|} \sum_{s \in \mathcal{V}_K} u_s^i \mathbf{g}_{K,s}^j \quad \text{and} \quad \mathcal{D}\mathcal{F}\mathcal{V}_K (\mathbf{U}_K) = \frac{1}{|K|} \sum_{i=1}^d \sum_{s \in \mathcal{V}_K} u_s^i \mathbf{g}_{K,s}^i. $$

If the $\mathbf{g}_{K,s}$ are computed using virtual element techniques as in [47], then the resulting scheme coincides with the first order virtual element scheme for linear elasticity of [33]. Once again one can obtain a flux based formulation by adding cell centered unknowns, setting again $\mathcal{M}_K^d(\mathbf{U}_K) = \mathbf{u}_K$. In this case, if the $\mathbf{g}_{K,s}$ are computed using virtual element techniques as in [47] we
obtain a flux based first order virtual element scheme for linear elasticity. More generally, for coefficients $g_{K,s}$ satisfying the requirements of VAG schemes, we obtain VAG linear elasticity schemes, flux based or not. Finally, notice that for vertices condition (H5) is automatically satisfied, providing another point of view on the reason why for face based schemes one needs to add extra stabilization terms while this is not necessary for vertex based schemes.

## 5. Error estimates

Many theoretical frameworks are available in the literature to establish error estimates for consistent approximations of symmetric and coercive problems. Without trying to be exhaustive, one could use for instance the classical two-point flux finite volume theory (see [48]), the theory of discontinuous Galerkin methods (see [49]), or the more recent gradient discretization framework (see [50]). For the diffusion problem which can be seen as the scalar version of our linear elasticity problem, details for the network element method can be found in [28, 30]. Apart from the technicalities arising from the fact that we do not necessarily have an underlying mesh, the main difficulty of the analysis lies in establishing a weak consistency estimate for the discrete gradient. However as we have already noticed, our vector operators coincide component-wise with their scalar counterparts:

$$ \Pi_K^d(\mathbf{U}_K)^i = \Pi_K(\mathbf{U}_K^i) \quad \text{and} \quad \nabla_K^d(\mathbf{U}_K)^i = \nabla_K(\mathbf{U}_K^i) \quad \text{and} \quad \mathcal{M}_K^d(\mathbf{U}_K)^i = \mathcal{M}_K(\mathbf{U}_K^i). $$

As the scalar versions have all the desired strong consistency properties it is clear that our vector ones will also satisfy the vector counterparts. Thus reproducing the proofs of [28, 30] it is equally clear that one can establish similar results for the linear elasticity problem. Unfortunately, even to simply state those results precisely, we need to recall some definitions linked to the network element framework. This is the reason why we start this section by recalling some notations and results from the network element framework, and then ends it by our general error estimates for which we establish the weak consistency estimate of the discrete gradient of vectors.

### 5.1. Quadrature families and network element interpolation

In [30], a family of functions $(\psi_K)_{K \in \mathcal{T}}$ is called a quadrature family if and only if 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 } \mathbf{x} \in \Omega, \quad \text{supp } \psi_K \subset \mathcal{B}_K = B(\mathbf{x}_K, \rho_K) \supset B_K. \quad (33) $$

For any $\mathbf{x} \in \mathbb{R}^d$, we denote

$$ \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} \text{card}(\mathcal{T}_{\mathbf{x}}^{\mathcal{B}}), \quad (34) $$

$$ \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)}. \quad (35) $$

True functions are reconstructed on $\mathbb{R}^d$ and $\Omega$ from the dofs using the $\psi_K$'s through:

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

Given a quadrature family $\psi_{\mathcal{T}} = (\psi_K)_{K \in \mathcal{T}}$, a family of functions $\zeta_{\mathcal{T}} = (\zeta_K)_{K \in \mathcal{T}}$ is called a family of weights adapted to $\psi_{\mathcal{T}}$ if and only if for any $K \in \mathcal{T}$, $\zeta_K \in L^\infty(\mathbb{R}^d)$ and:

$$ \int_{\Omega} \zeta_K = m_K, \quad \text{for a.e } \mathbf{x} \in \Omega, \quad \text{supp } \zeta_K \subset \mathcal{B}_K, \quad M_{\zeta} = \max_{K \in \mathcal{T}} \|\zeta_K\|_{L^{\infty}(\Omega)}. \quad (36) $$
obtain a flux based first order virtual element scheme for linear elasticity. More generally, for coefficients $g_{K,s}$ satisfying the requirements of VAG schemes, we obtain VAG linear elasticity schemes, flux based or not. Finally, notice that for vertices condition (H5) is automatically satisfied, providing another point of view on the reason why for face based schemes one needs to add extra stabilization terms while this is not necessary for vertex based schemes.

## 5. Error estimates

Many theoretical frameworks are available in the literature to establish error estimates for consistent approximations of symmetric and coercive problems. Without trying to be exhaustive, one could use for instance the classical two-point flux finite volume theory (see [48]), the theory of discontinuous Galerkin methods (see [49]), or the more recent gradient discretization framework (see [50]). For the diffusion problem which can be seen as the scalar version of our linear elasticity problem, details for the network element method can be found in [28, 30]. Apart from the technicalities arising from the fact that we do not necessarily have an underlying mesh, the main difficulty of the analysis lies in establishing a weak consistency estimate for the discrete gradient. However as we have already noticed, our vector operators coincide component-wise with their scalar counterparts:

$$ \Pi_K^d(\mathbf{U}_K)^i = \Pi_K(\mathbf{U}_K^i) \quad \text{and} \quad \nabla_K^d(\mathbf{U}_K)^i = \nabla_K(\mathbf{U}_K^i) \quad \text{and} \quad \mathcal{M}_K^d(\mathbf{U}_K)^i = \mathcal{M}_K(\mathbf{U}_K^i). $$

As the scalar versions have all the desired strong consistency properties it is clear that our vector ones will also satisfy the vector counterparts. Thus reproducing the proofs of [28, 30] it is equally clear that one can establish similar results for the linear elasticity problem. Unfortunately, even to simply state those results precisely, we need to recall some definitions linked to the network element framework. This is the reason why we start this section by recalling some notations and results from the network element framework, and then ends it by our general error estimates for which we establish the weak consistency estimate of the discrete gradient of vectors.

### 5.1. Quadrature families and network element interpolation

In [30], a family of functions $(\psi_K)_{K \in \mathcal{T}}$ is called a quadrature family if and only if 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 } \mathbf{x} \in \Omega, \quad \text{supp } \psi_K \subset \mathcal{B}_K = B(\mathbf{x}_K, \rho_K) \supset B_K. \quad (33) $$

For any $x \in \mathbb{R}^d$, we denote

$$ \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} \text{card}(\mathcal{T}_{\mathbf{x}}^{\mathcal{B}}), \quad (34) $$

$$ \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)}. \quad (35) $$

True functions are reconstructed on $\mathbb{R}^d$ and $\Omega$ from the dofs using the $\psi_K$'s through:

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

Given a quadrature family $\psi_{\mathcal{T}} = (\psi_K)_{K \in \mathcal{T}}$, a family of functions $\zeta_{\mathcal{T}} = (\zeta_K)_{K \in \mathcal{T}}$ is called a family of weights adapted to $\psi_{\mathcal{T}}$ if and only if for any $K \in \mathcal{T}$, $\zeta_K \in L^\infty(\mathbb{R}^d)$ and:

$$ \int_{\Omega} \zeta_K = m_K, \quad \text{for a.e } \mathbf{x} \in \Omega, \quad \text{supp } \zeta_K \subset \mathcal{B}_K, \quad M_{\zeta} = \max_{K \in \mathcal{T}} \|\zeta_K\|_{L^{\infty}(\Omega)}. \quad (36) $$
As $\Omega$ is assumed Lipschitz, using Stein's extension theorem [51] 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$. We define the operator $\mathcal{I} : H^1(\Omega) \to X_{\mathcal{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}, \qquad (37)$$

and to any $\sigma \in \mathcal{F}$ is associated a radius $r_\sigma > 0$ such that $B_\sigma \subset \mathcal{B}_K$, where we denote $B_\sigma = B(x_\sigma, r_\sigma)$, as well as:

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

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}^0 : H_0^1(\Omega) \to X_{\mathcal{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} \qquad (38)$$

For vectors $v \in H^1(\Omega)^d$ (resp. tensors in $\sigma \in H^1(\Omega)^{d \times d}$), we of course denote $\mathcal{I}(v) = (\mathcal{I}(v^i))_{1 \le i \le d}$ (resp. $\mathcal{I}(\sigma) = (\mathcal{I}(\sigma_{ij}))_{1 \le i,j \le d}$) i.e. the operator $\mathcal{I}$ applied component by component, the same being done for vectors $v \in H_0^1(\Omega)^d$ and operator $\mathcal{I}^0$.

## 5.2. Statement of the error estimates

We can now state our main error estimation result:

**Proposition 1.** Let $(\mathcal{N}, \mathcal{G})$ be a discretization network and associated admissible geometry for which (H5) holds, and let $\mathbf{U}$ be the solution of the associated problem (14). Then, if $\psi_{\mathcal{F}} = (\psi_K)_{K \in \mathcal{F}}$ is a quadrature family and $\zeta_{\mathcal{F}} = (\zeta_K)_{K \in \mathcal{F}}$ is a family of weights adapted to $\psi_{\mathcal{F}}$ we have the following error estimates:

$$\| \mathbf{U} - \mathcal{I}^0(u) \|_{X} \le C \mathcal{E}_h(\mathbf{u}, \mathbf{f}), \qquad (39)$$

and

$$\| \mathbf{u} - \Pi_{\mathcal{F}}^d(\mathbf{U}) \|_{L^2(\Omega)} \le C \mathcal{E}_h(\mathbf{u}, \mathbf{f}) \quad \text{and} \quad \| \nabla \mathbf{u} - \nabla_{\mathcal{F}}^d(\mathbf{U}) \|_{L^2(\Omega)} \le C \mathcal{E}_h(\mathbf{u}, \mathbf{f}), \qquad (40)$$

where the constant $C > 0$ depends on $u, \mu, \lambda, d, \theta_{\mathcal{F}}, \theta_{\mathcal{G}}, \theta_{\Pi}, \theta_{\mathcal{A}}, \theta_{\mathcal{M}}, M_{\zeta}, M_{\psi}, \kappa_{\psi}, \eta_{\psi}, \theta_{\mathcal{I}}, C_{\mathcal{K},\mathcal{N}}$ and $\Omega$ but not on $h$, and where the abstract error is given by

$$\begin{aligned} \mathcal{E}_h(\mathbf{u}, \mathbf{f}) = & \left( \sum_{K \in \mathcal{F}} \| \boldsymbol{\epsilon}_K(\mathcal{I}^0(\mathbf{u})) - \boldsymbol{\epsilon}(\mathbf{u}) \|_{L^2(\mathcal{B}_K \cap \Omega)}^{2 \times d} \right)^{1/2} \\ & + \left( \sum_{K \in \mathcal{F}} \| \mathcal{M}_K^d(\mathcal{I}^0(\mathbf{u})) - \mathbf{u} \|_{L^2(\mathcal{B}_K \cap \Omega)}^{2} \right)^{1/2} \\ & + \left( \sum_{K \in \mathcal{F}} h_K^{-2} \| \Pi_K^d(\mathcal{I}^0(\mathbf{u})) - \mathbf{u} \|_{L^2(\mathcal{B}_K \cap \Omega)}^{2} \right)^{1/2} \\ & + \left( \sum_{K \in \mathcal{F}} \| \mathcal{I}_K(\boldsymbol{\sigma}(\mathbf{u})) - \boldsymbol{\sigma}(\mathbf{u}) \|_{L^2(\mathcal{B}_K \cap \Omega)}^{2} \right)^{1/2} \\ & + \left( \sum_{K \in \mathcal{F}} \sum_{\sigma \in \mathcal{F}_K} \| \mathcal{I}_{\sigma}(\boldsymbol{\sigma}(\mathbf{u})) - \boldsymbol{\sigma}(\mathbf{u}) \|_{L^2(\mathcal{B}_K \cap \Omega)}^{2} \right)^{1/2} \\ & + h^p + \mathcal{E}_h^S(\zeta_{\mathcal{F}}, \mathbf{u}, \mathbf{f}), \end{aligned}$$
As $\Omega$ is assumed Lipschitz, using Stein's extension theorem [51] 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$. We define the operator $\mathcal{I} : H^1(\Omega) \to X_{\mathcal{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}, \qquad (37)$$

and to any $\sigma \in \mathcal{F}$ is associated a radius $r_\sigma > 0$ such that $B_\sigma \subset \mathcal{B}_K$, where we denote $B_\sigma = B(x_\sigma, r_\sigma)$, as well as:

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

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}^0 : H_0^1(\Omega) \to X_{\mathcal{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} \qquad (38)$$

For vectors $v \in H^1(\Omega)^d$ (resp. tensors in $\sigma \in H^1(\Omega)^{d \times d}$), we of course denote $\mathcal{I}(v) = (\mathcal{I}(v^i))_{1 \le i \le d}$ (resp. $\mathcal{I}(\sigma) = (\mathcal{I}(\sigma_{ij}))_{1 \le i,j \le d}$) i.e. the operator $\mathcal{I}$ applied component by component, the same being done for vectors $v \in H_0^1(\Omega)^d$ and operator $\mathcal{I}^0$.

## 5.2. Statement of the error estimates

We can now state our main error estimation result:

**Proposition 1.** Let $(\mathcal{N}, \mathcal{G})$ be a discretization network and associated admissible geometry for which (H5) holds, and let $\mathbf{U}$ be the solution of the associated problem (14). Then, if $\psi_{\mathcal{F}} = (\psi_K)_{K \in \mathcal{F}}$ is a quadrature family and $\zeta_{\mathcal{F}} = (\zeta_K)_{K \in \mathcal{F}}$ is a family of weights adapted to $\psi_{\mathcal{F}}$ we have the following error estimates:

$$\| \mathbf{U} - \mathcal{I}^0(u) \|_{X} \le C \mathcal{E}_h(\mathbf{u}, \mathbf{f}), \qquad (39)$$

and

$$\| \mathbf{u} - \Pi_{\mathcal{F}}^d(\mathbf{U}) \|_{L^2(\Omega)} \le C \mathcal{E}_h(\mathbf{u}, \mathbf{f}) \quad \text{and} \quad \| \nabla \mathbf{u} - \nabla_{\mathcal{F}}^d(\mathbf{U}) \|_{L^2(\Omega)} \le C \mathcal{E}_h(\mathbf{u}, \mathbf{f}), \qquad (40)$$

where the constant $C > 0$ depends on $u, \mu, \lambda, d, \theta_{\mathcal{F}}, \theta_{\mathcal{G}}, \theta_{\Pi}, \theta_{\mathcal{A}}, \theta_{\mathcal{M}}, M_{\zeta}, M_{\psi}, \kappa_{\psi}, \eta_{\psi}, \theta_{\mathcal{I}}, C_{\mathcal{K},\mathcal{N}}$ and $\Omega$ but not on $h$, and where the abstract error is given by

$$\begin{aligned} \mathcal{E}_h(\mathbf{u}, \mathbf{f}) = & \left( \sum_{K \in \mathcal{F}} \| \boldsymbol{\epsilon}_K(\mathcal{I}^0(\mathbf{u})) - \boldsymbol{\epsilon}(\mathbf{u}) \|_{L^2(\mathcal{B}_K \cap \Omega)}^{2 \times d} \right)^{1/2} \\ & + \left( \sum_{K \in \mathcal{F}} \| \mathcal{M}_K^d(\mathcal{I}^0(\mathbf{u})) - \mathbf{u} \|_{L^2(\mathcal{B}_K \cap \Omega)}^{2} \right)^{1/2} \\ & + \left( \sum_{K \in \mathcal{F}} h_K^{-2} \| \Pi_K^d(\mathcal{I}^0(\mathbf{u})) - \mathbf{u} \|_{L^2(\mathcal{B}_K \cap \Omega)}^{2} \right)^{1/2} \\ & + \left( \sum_{K \in \mathcal{F}} \| \mathcal{I}_K(\boldsymbol{\sigma}(\mathbf{u})) - \boldsymbol{\sigma}(\mathbf{u}) \|_{L^2(\mathcal{B}_K \cap \Omega)}^{2} \right)^{1/2} \\ & + \left( \sum_{K \in \mathcal{F}} \sum_{\sigma \in \mathcal{F}_K} \| \mathcal{I}_{\sigma}(\boldsymbol{\sigma}(\mathbf{u})) - \boldsymbol{\sigma}(\mathbf{u}) \|_{L^2(\mathcal{B}_K \cap \Omega)}^{2} \right)^{1/2} \\ & + h^p + \mathcal{E}_h^S(\zeta_{\mathcal{F}}, \mathbf{u}, \mathbf{f}), \end{aligned}$$
where

$$
\mathcal{E}_h^S(\zeta_{\mathcal{F}}, \mathbf{u}, \mathbf{f}) = \left( \sum_{K \in \mathcal{F}} \int_{\mathcal{B}_K \cap \Omega} |\zeta_K| |\mathbf{f}_K - \mathbf{f}|^2 \right)^{1/2} + \left( \sum_{K \in \mathcal{F}} \int_{\mathcal{B}_K \cap \Omega} |\zeta_K| |\mathcal{D} \mathcal{I} \mathcal{V}_K^d(\mathcal{I}(\boldsymbol{\sigma}(\mathbf{u}))) - \operatorname{div}(\boldsymbol{\sigma}(\mathbf{u}))|^2 \right)^{1/2},
$$

and

$$
\mathcal{I}_K(\boldsymbol{\sigma}(\mathbf{u})) = \frac{1}{|\mathcal{B}_K \cap \Omega|} \int_{\mathcal{B}_K \cap \Omega} \boldsymbol{\sigma}(\mathbf{u}) \quad \forall K \in \mathcal{F},
$$

and

$$
\mathcal{D} \mathcal{I} \mathcal{V}_K^d(\mathcal{I}(\boldsymbol{\sigma})) = \frac{1}{m_K} \sum_{i=1}^d \sum_{j=1}^d \sum_{\boldsymbol{\sigma} \in \mathcal{F}_K} \mathcal{I}_{\boldsymbol{\sigma}}(\boldsymbol{\sigma}_{ij}) \eta_{K,\sigma}^j e_i = \frac{1}{m_K} \sum_{\boldsymbol{\sigma} \in \mathcal{F}_K} \mathcal{I}_{\boldsymbol{\sigma}}(\boldsymbol{\sigma}) \eta_{K,\sigma}.
$$

Moreover, for any $K \in \mathcal{T}$

$$
\| \mathbf{u} - \mathcal{M}_K^d(\mathbf{U}) \|_{L^2(B_K \cap \Omega)^d} \le C \mathcal{E}_h(\mathbf{u}, \mathbf{f}) \quad \text{and} \quad \| \nabla \mathbf{u} - \nabla_K^d(\mathbf{U}) \|_{L^2(B_K \cap \Omega)^{d \times d}} \le C \mathcal{E}_h(\mathbf{u}, \mathbf{f}). \quad (41)
$$

As announced at the beginning of this section, the proof of proposition 1 is a straightforward
adaptation of its counterpart in the scalar case [28, 30]. Those proofs in turn can be mostly con-
sidered to coincide with well established frameworks, apart from the network element additional
difficulty of requiring quadrature families to bridge the gap between the dof space and the func-
tional spaces. However a key point is that in the network element framework, the second term of
$\mathcal{E}_h^S(\zeta_{\mathcal{F}}, \mathbf{u}, \mathbf{f})$, which measures in some sense the weighted error on the strong from of prob-
lem (7) cannot be eliminated by a careful choice of $\mathcal{I}$, contrary to what happens in general for
mesh-based discretizations for which the discrete divergence is equal to the projection of the ex-
act divergence on constant functions.The origin of this additional term is the following weak con-
sistency estimate for the discrete gradient, which is the only part of the proof of proposition 1 we
will detail here as it is the only true difference with the scalar case:

**Lemma 2 (Weak consistency estimate).** Let $(\mathcal{N}, \mathcal{G})$ be an admissible network, and associated admissible geometry. For any $\mathbf{V} \in X_{\mathcal{N},0}$ and any $\boldsymbol{\sigma} \in H(\operatorname{div}, \Omega)$, any quadrature family $\psi_{\mathcal{F}} = (\psi_K)_{K \in \mathcal{F}}$ and any family of weights $\zeta_{\mathcal{F}} = (\zeta_K)_{K \in \mathcal{F}}$ adapted to $\psi_{\mathcal{F}}$, we denote:

$$
\mathcal{R}(\psi_{\mathcal{F}}, \zeta_{\mathcal{F}}, \boldsymbol{\sigma}, \mathbf{V}) = \int_{\Omega} \boldsymbol{\sigma} : \nabla_{\mathcal{F}}^d(\mathbf{V}) + \sum_{K \in \mathcal{F}} \int_{\mathcal{B}_K \cap \Omega} \zeta_K \operatorname{div}(\boldsymbol{\sigma}) \cdot \mathcal{M}_K^d(\mathbf{V}). \quad (42)
$$

Then, there exists $C > 0$ depending on $\boldsymbol{\sigma}$, $d$, $\theta_{\mathcal{F}}$, $\theta_{\mathcal{F}}$, $\theta_{\Pi}$, $\theta_{\mathcal{A}}$, $\theta_{\mathcal{M}}$, $M_{\zeta}$, $M_{\psi}$, $\kappa_{\psi}$, $\eta_{\psi}$, $\theta_{\mathcal{I}}$ and $\Omega$ but
not on $h$ such that

$$
|\mathcal{R}(\psi_{\mathcal{F}}, \zeta_{\mathcal{F}}, \boldsymbol{\sigma}, \mathbf{V})| \le C \mathcal{E}_h^{\mathcal{R}}(\zeta_{\mathcal{F}}, \boldsymbol{\sigma}) \|\mathbf{V}\|_X, \quad (43)
$$

where

$$
\begin{align}
\mathcal{E}_h^{\mathcal{R}}(\zeta_{\mathcal{F}}, \boldsymbol{\sigma}) &= \left( \sum_{K \in \mathcal{F}} \| \mathcal{I}_K(\boldsymbol{\sigma}) - \boldsymbol{\sigma} \|_{L^2(\mathcal{B}_K \cap \Omega)^{d \times d}}^2 \right)^{1/2} \nonumber \\
&\quad + \left( \sum_{K \in \mathcal{F}} \sum_{\boldsymbol{\sigma} \in \mathcal{F}_K} \| \mathcal{I}_{\boldsymbol{\sigma}}(\boldsymbol{\sigma}) - \boldsymbol{\sigma} \|_{L^2(\mathcal{B}_K \cap \Omega)^{d \times d}}^2 \right)^{1/2} \nonumber \\
&\quad + h^p \nonumber \\
&\quad + \left( \sum_{K \in \mathcal{F}} \int_{\mathcal{B}_K \cap \Omega} |\zeta_K| |\mathcal{D} \mathcal{I} \mathcal{V}_K^{\mathcal{G}}(\mathcal{I}(\boldsymbol{\sigma})) - \operatorname{div}(\boldsymbol{\sigma})|^2 \right)^{1/2}. \tag{44}
\end{align}
$$

**Proof.** We start by writing

$$
\sum_{K \in \mathcal{F}} \int_{\mathcal{B}_K \cap \Omega} \zeta_K \operatorname{div}(\boldsymbol{\sigma}) \cdot \mathcal{M}_K^d(\mathbf{V}) \\
= \sum_{K \in \mathcal{F}} \int_{\mathcal{B}_K \cap \Omega} \zeta_K \left( \operatorname{div}(\boldsymbol{\sigma}) - \mathcal{D} \mathcal{I} \mathcal{V}_K^d(\boldsymbol{\mathcal{I}}(\boldsymbol{\sigma})) \right) \cdot \mathcal{M}_K^d(\mathbf{V}) + \sum_{K \in \mathcal{F}} \sum_{\boldsymbol{\sigma} \in \mathcal{F}_K} (\mathcal{I}_{\boldsymbol{\sigma}}(\boldsymbol{\sigma}) \eta_{K,\sigma}) \cdot \mathcal{M}_K^d(\mathbf{V}),
$$
where

$$
\mathcal{E}_h^S(\zeta_{\mathcal{F}}, \mathbf{u}, \mathbf{f}) = \left( \sum_{K \in \mathcal{F}} \int_{\mathcal{B}_K \cap \Omega} |\zeta_K| |\mathbf{f}_K - \mathbf{f}|^2 \right)^{1/2} + \left( \sum_{K \in \mathcal{F}} \int_{\mathcal{B}_K \cap \Omega} |\zeta_K| |\mathcal{D} \mathcal{I} \mathcal{V}_K^d(\mathcal{I}(\boldsymbol{\sigma}(\mathbf{u}))) - \operatorname{div}(\boldsymbol{\sigma}(\mathbf{u}))|^2 \right)^{1/2},
$$

and

$$
\mathcal{I}_K(\boldsymbol{\sigma}(\mathbf{u})) = \frac{1}{|\mathcal{B}_K \cap \Omega|} \int_{\mathcal{B}_K \cap \Omega} \boldsymbol{\sigma}(\mathbf{u}) \quad \forall K \in \mathcal{F},
$$

and

$$
\mathcal{D} \mathcal{I} \mathcal{V}_K^d(\mathcal{I}(\boldsymbol{\sigma})) = \frac{1}{m_K} \sum_{i=1}^d \sum_{j=1}^d \sum_{\boldsymbol{\sigma} \in \mathcal{F}_K} \mathcal{I}_{\boldsymbol{\sigma}}(\boldsymbol{\sigma}_{ij}) \eta_{K,\sigma}^j e_i = \frac{1}{m_K} \sum_{\boldsymbol{\sigma} \in \mathcal{F}_K} \mathcal{I}_{\boldsymbol{\sigma}}(\boldsymbol{\sigma}) \eta_{K,\sigma}.
$$

Moreover, for any $K \in \mathcal{F}$

$$
\| \mathbf{u} - \mathcal{M}_K^d(\mathbf{U}) \|_{L^2(B_K \cap \Omega)^d} \le C \mathcal{E}_h(\mathbf{u}, \mathbf{f}) \quad \text{and} \quad \| \nabla \mathbf{u} - \nabla_K^d(\mathbf{U}) \|_{L^2(B_K \cap \Omega)^{d \times d}} \le C \mathcal{E}_h(\mathbf{u}, \mathbf{f}). \quad (41)
$$

As announced at the beginning of this section, the proof of proposition 1 is a straightforward
adaptation of its counterpart in the scalar case [28, 30]. Those proofs in turn can be mostly con-
sidered to coincide with well established frameworks, apart from the network element additional
difficulty of requiring quadrature families to bridge the gap between the dof space and the func-
tional spaces. However a key point is that in the network element framework, the second term of
$\mathcal{E}_h^S(\zeta_{\mathcal{F}}, \mathbf{u}, \mathbf{f})$, which measures in some sense the weighted error on the strong from of prob-
lem (7) cannot be eliminated by a careful choice of $\mathcal{I}$, contrary to what happens in general for
mesh-based discretizations for which the discrete divergence is equal to the projection of the ex-
act divergence on constant functions.The origin of this additional term is the following weak con-
sistency estimate for the discrete gradient, which is the only part of the proof of proposition 1 we
will detail here as it is the only true difference with the scalar case:

**Lemma 2 (Weak consistency estimate).** Let $(\mathcal{N}, \mathcal{G})$ be an admissible network, and associated admissible geometry. For any $\mathbf{V} \in X_{\mathcal{N},0}$ and any $\boldsymbol{\sigma} \in H(\operatorname{div}, \Omega)$, any quadrature family $\psi_{\mathcal{F}} = (\psi_K)_{K \in \mathcal{F}}$ and any family of weights $\zeta_{\mathcal{F}} = (\zeta_K)_{K \in \mathcal{F}}$ adapted to $\psi_{\mathcal{F}}$, we denote:

$$
\mathcal{R}(\psi_{\mathcal{F}}, \zeta_{\mathcal{F}}, \boldsymbol{\sigma}, \mathbf{V}) = \int_{\Omega} \boldsymbol{\sigma} : \nabla_{\mathcal{F}}^d(\mathbf{V}) + \sum_{K \in \mathcal{F}} \int_{\mathcal{B}_K \cap \Omega} \zeta_K \operatorname{div}(\boldsymbol{\sigma}) \cdot \mathcal{M}_K^d(\mathbf{V}). \quad (42)
$$

Then, there exists $C > 0$ depending on $\boldsymbol{\sigma}$, $d$, $\theta_{\mathcal{F}}$, $\theta_{\mathcal{F}}$, $\theta_{\Pi}$, $\theta_{\mathcal{A}}$, $\theta_{\mathcal{M}}$, $M_{\zeta}$, $M_{\psi}$, $\kappa_{\psi}$, $\eta_{\psi}$, $\theta_{\mathcal{I}}$ and $\Omega$ but
not on $h$ such that

$$
|\mathcal{R}(\psi_{\mathcal{F}}, \zeta_{\mathcal{F}}, \boldsymbol{\sigma}, \mathbf{V})| \le C \mathcal{E}_h^{\mathcal{R}}(\zeta_{\mathcal{F}}, \boldsymbol{\sigma}) \|\mathbf{V}\|_X, \quad (43)
$$

where

$$
\begin{align}
\mathcal{E}_h^{\mathcal{R}}(\zeta_{\mathcal{F}}, \boldsymbol{\sigma}) &= \left( \sum_{K \in \mathcal{F}} \| \mathcal{I}_K(\boldsymbol{\sigma}) - \boldsymbol{\sigma} \|_{L^2(\mathcal{B}_K \cap \Omega)^{d \times d}}^2 \right)^{1/2} \nonumber \\
&\quad + \left( \sum_{K \in \mathcal{F}} \sum_{\boldsymbol{\sigma} \in \mathcal{F}_K} \| \mathcal{I}_{\boldsymbol{\sigma}}(\boldsymbol{\sigma}) - \boldsymbol{\sigma} \|_{L^2(\mathcal{B}_K \cap \Omega)^{d \times d}}^2 \right)^{1/2} \nonumber \\
&\quad + h^p \nonumber \\
&\quad + \left( \sum_{K \in \mathcal{F}} \int_{\mathcal{B}_K \cap \Omega} |\zeta_K| |\mathcal{D} \mathcal{I} \mathcal{V}_K^d(\mathcal{I}(\boldsymbol{\sigma})) - \operatorname{div}(\boldsymbol{\sigma})|^2 \right)^{1/2}. \tag{44}
\end{align}
$$

**Proof.** We start by writing

$$
\sum_{K \in \mathcal{F}} \int_{\mathcal{B}_K \cap \Omega} \zeta_K \operatorname{div}(\boldsymbol{\sigma}) \cdot \mathcal{M}_K^d(\mathbf{V}) \\
= \sum_{K \in \mathcal{F}} \int_{\mathcal{B}_K \cap \Omega} \zeta_K \left( \operatorname{div}(\boldsymbol{\sigma}) - \mathcal{D} \mathcal{I} \mathcal{V}_K^d(\boldsymbol{\mathcal{I}}(\boldsymbol{\sigma})) \right) \cdot \mathcal{M}_K^d(\mathbf{V}) + \sum_{K \in \mathcal{F}} \sum_{\boldsymbol{\sigma} \in \mathcal{F}_K} (\mathcal{I}_{\boldsymbol{\sigma}}(\boldsymbol{\sigma}) \eta_{K,\sigma}) \cdot \mathcal{M}_K^d(\mathbf{V}),
$$
since by definition:

$$ \int_{\mathcal{B}_K \cap \Omega} \zeta_K = \int_{\Omega} \zeta_K = m_K \quad \text{and} \quad \mathcal{D} \mathcal{I} \mathcal{V}_K^d(\mathcal{I}(\boldsymbol{\sigma})) = \frac{1}{m_K} \sum_{\sigma \in \mathcal{F}_K} \mathcal{I}_{\sigma}(\boldsymbol{\sigma}) \boldsymbol{\eta}_{K, \sigma}. $$

As $v_\sigma = 0$ for any $\sigma \in \mathcal{F}_{\text{ext}}$ and using the approximate conservation properties $\sum_{K \in \mathcal{F}_\sigma} \boldsymbol{\eta}_{K, \sigma} = \boldsymbol{\varepsilon}_\sigma$ of the geometry, we get:

$$ \sum_{K \in \mathcal{F}} \sum_{\sigma \in \mathcal{F}_K} (\mathcal{I}_{\sigma}(\boldsymbol{\sigma}) \boldsymbol{\eta}_{K, \sigma}) \cdot \mathcal{M}_K^d(\mathbf{V}) = \sum_{K \in \mathcal{F}} \sum_{\sigma \in \mathcal{F}_K} (\mathcal{I}_{\sigma}(\boldsymbol{\sigma}) \boldsymbol{\eta}_{K, \sigma}) \cdot (\mathcal{M}_K^d(\mathbf{V}) - \mathbf{v}_{\sigma}) + \sum_{\sigma \in \mathcal{F}_{\text{int}}} (\mathcal{I}_{\sigma}(\boldsymbol{\sigma}) \boldsymbol{\varepsilon}_{\sigma}) \cdot \mathbf{v}_{\sigma}, $$

then we get:

$$ \int_{\Omega} \mathrm{div}(\boldsymbol{\sigma}) \cdot \Pi_{\mathcal{F}}^{d}(\mathbf{V}) = \sum_{K \in \mathcal{F}} \int_{\mathcal{B}_K \cap \Omega} \zeta_K (\mathrm{div}(\boldsymbol{\sigma}) - \mathcal{D} \mathcal{I} \mathcal{V}_K(\mathcal{I}(\boldsymbol{\sigma}))) \cdot \mathcal{M}_K^d(\mathbf{V}) \\ + \sum_{K \in \mathcal{F}} \sum_{\sigma \in \mathcal{F}_K} ((\mathcal{I}_{\sigma}(\boldsymbol{\sigma}) - \mathcal{I}_K(\boldsymbol{\sigma})) \boldsymbol{\eta}_{K, \sigma}) \cdot (\mathcal{M}_K^d(\mathbf{V}) - \mathbf{v}_{\sigma}) \\ - \sum_{K \in \mathcal{F}} \sum_{\sigma \in \mathcal{F}_K} (\mathcal{I}_K(\boldsymbol{\sigma}) \boldsymbol{\eta}_{K, \sigma}) \cdot \mathbf{v}_{\sigma} + \sum_{\sigma \in \mathcal{F}_{\text{int}}} (\mathcal{I}_{\sigma}(\boldsymbol{\sigma}) \boldsymbol{\varepsilon}_{\sigma}) \cdot \mathbf{v}_{\sigma} \\ + \sum_{K \in \mathcal{F}} m_K (\mathcal{I}_K(\boldsymbol{\sigma}) \boldsymbol{\varepsilon}_K^0) \cdot \mathcal{M}_K^d(\mathbf{V}), $$

using this time the zero order approximate consistency property:

$$ \sum_{\sigma \in \mathcal{F}_K} \boldsymbol{\eta}_{K, \sigma} = \boldsymbol{\varepsilon}_K^0. $$

Notice that:

$$ - \sum_{K \in \mathcal{F}} \sum_{\sigma \in \mathcal{F}_K} (\mathcal{I}_K(\boldsymbol{\sigma}) \boldsymbol{\eta}_{K, \sigma}) \cdot \mathbf{v}_{\sigma} = - \sum_{K \in \mathcal{F}} \sum_{\sigma \in \mathcal{F}_K} \sum_{i=1}^{d} \sum_{j=1}^{d} \mathcal{I}_K(\sigma_{ij}) \boldsymbol{\eta}_{K, \sigma}^{j} v_{\sigma}^{i} = - \sum_{K \in \mathcal{F}} m_K \mathcal{I}_K(\boldsymbol{\sigma}) : \nabla \Pi_K^d(\mathbf{V}). $$

Consequently, we have using the definition of $\nabla_{\mathcal{F}}^d(\mathbf{V})$:

$$ \int_{\Omega} \mathrm{div}(\boldsymbol{\sigma}) \cdot \Pi_{\mathcal{F}}^{d}(\mathbf{V}) = \sum_{K \in \mathcal{F}} \int_{\mathcal{B}_K \cap \Omega} \zeta_K (\mathrm{div}(\boldsymbol{\sigma}) - \mathcal{D} \mathcal{I} \mathcal{V}_K(\mathcal{I}(\boldsymbol{\sigma}))) \cdot \mathcal{M}_K^d(\mathbf{V}) \\ + \sum_{K \in \mathcal{F}} \sum_{\sigma \in \mathcal{F}_K} ((\mathcal{I}_{\sigma}(\boldsymbol{\sigma}) - \mathcal{I}_K(\boldsymbol{\sigma})) \boldsymbol{\eta}_{K, \sigma}) \cdot (\mathcal{M}_K^d(\mathbf{V}) - \mathbf{v}_{\sigma}) \\ - \sum_{K \in \mathcal{F}} m_K \mathcal{I}_K(\boldsymbol{\sigma}) : \nabla \Pi_K^d(\mathbf{V}) + \sum_{\sigma \in \mathcal{F}_{\text{int}}} (\mathcal{I}_{\sigma}(\boldsymbol{\sigma}) \boldsymbol{\varepsilon}_{\sigma}) \cdot \mathbf{v}_{\sigma} \\ + \sum_{K \in \mathcal{F}} m_K (\mathcal{I}_K(\boldsymbol{\sigma}) \boldsymbol{\varepsilon}_K^0) \cdot \mathcal{M}_K^d(\mathbf{V}) \\ = \sum_{K \in \mathcal{F}} \int_{\mathcal{B}_K \cap \Omega} \zeta_K (\mathrm{div}(\boldsymbol{\sigma}) - \mathcal{D} \mathcal{I} \mathcal{V}_K(\mathcal{I}(\boldsymbol{\sigma}))) \cdot \mathcal{M}_K^d(\mathbf{V}) \\ + \sum_{K \in \mathcal{F}} \sum_{\sigma \in \mathcal{F}_K} ((\mathcal{I}_{\sigma}(\boldsymbol{\sigma}) - \mathcal{I}_K(\boldsymbol{\sigma})) \boldsymbol{\eta}_{K, \sigma}) \cdot (\mathcal{M}_K^d(\mathbf{V}) - \mathbf{v}_{\sigma}) \\ + \int_{\Omega} \sum_{K \in \mathcal{F}} \psi_K(\boldsymbol{\sigma} - \mathcal{I}_K(\boldsymbol{\sigma})) \cdot \nabla \Pi_K^d(\mathbf{V}) \\ - \int_{\Omega} \boldsymbol{\sigma} : \nabla_{\mathcal{F}}^d(\mathbf{V}) + \sum_{\sigma \in \mathcal{F}_{\text{int}}} (\mathcal{I}_{\sigma}(\boldsymbol{\sigma}) \boldsymbol{\varepsilon}_{\sigma}) \cdot \mathbf{v}_{\sigma} + \sum_{K \in \mathcal{F}} m_K (\mathcal{I}_K(\boldsymbol{\sigma}) \boldsymbol{\varepsilon}_K^0) \cdot \mathcal{M}_K^d(\mathbf{V}), $$
since by definition:

$$ \int_{\mathcal{B}_K \cap \Omega} \zeta_K = \int_{\Omega} \zeta_K = m_K \quad \text{and} \quad \mathcal{D} \mathcal{I} \mathcal{V}_K^d(\mathcal{I}(\boldsymbol{\sigma})) = \frac{1}{m_K} \sum_{\sigma \in \mathcal{F}_K} \mathcal{I}_{\sigma}(\boldsymbol{\sigma}) \boldsymbol{\eta}_{K, \sigma}. $$

As $v_{\sigma} = 0$ for any $\sigma \in \mathcal{F}_{\text{ext}}$ and using the approximate conservation properties $\sum_{K \in \mathcal{F}_{\sigma}} \boldsymbol{\eta}_{K, \sigma} = \boldsymbol{\varepsilon}_{\sigma}$ of the geometry, we get:

$$ \sum_{K \in \mathcal{F}} \sum_{\sigma \in \mathcal{F}_K} (\mathcal{I}_{\sigma}(\boldsymbol{\sigma}) \boldsymbol{\eta}_{K, \sigma}) \cdot \mathcal{M}_K^d(\mathbf{V}) = \sum_{K \in \mathcal{F}} \sum_{\sigma \in \mathcal{F}_K} (\mathcal{I}_{\sigma}(\boldsymbol{\sigma}) \boldsymbol{\eta}_{K, \sigma}) \cdot (\mathcal{M}_K^d(\mathbf{V}) - \mathbf{v}_{\sigma}) + \sum_{\sigma \in \mathcal{F}_{\text{int}}} (\mathcal{I}_{\sigma}(\boldsymbol{\sigma}) \boldsymbol{\varepsilon}_{\sigma}) \cdot \mathbf{v}_{\sigma}, $$

then we get:

$$ \int_{\Omega} \mathrm{div}(\boldsymbol{\sigma}) \cdot \Pi_{\mathcal{F}}^{d}(\mathbf{V}) = \sum_{K \in \mathcal{F}} \int_{\mathcal{B}_K \cap \Omega} \zeta_K (\mathrm{div}(\boldsymbol{\sigma}) - \mathcal{D} \mathcal{I} \mathcal{V}_K(\mathcal{I}(\boldsymbol{\sigma}))) \cdot \mathcal{M}_K^d(\mathbf{V}) \\ + \sum_{K \in \mathcal{F}} \sum_{\sigma \in \mathcal{F}_K} ((\mathcal{I}_{\sigma}(\boldsymbol{\sigma}) - \mathcal{I}_K(\boldsymbol{\sigma})) \boldsymbol{\eta}_{K, \sigma}) \cdot (\mathcal{M}_K^d(\mathbf{V}) - \mathbf{v}_{\sigma}) \\ - \sum_{K \in \mathcal{F}} \sum_{\sigma \in \mathcal{F}_K} (\mathcal{I}_K(\boldsymbol{\sigma}) \boldsymbol{\eta}_{K, \sigma}) \cdot \mathbf{v}_{\sigma} + \sum_{\sigma \in \mathcal{F}_{\text{int}}} (\mathcal{I}_{\sigma}(\boldsymbol{\sigma}) \boldsymbol{\varepsilon}_{\sigma}) \cdot \mathbf{v}_{\sigma} \\ + \sum_{K \in \mathcal{F}} m_K (\mathcal{I}_K(\boldsymbol{\sigma}) \boldsymbol{\varepsilon}_K^0) \cdot \mathcal{M}_K^d(\mathbf{V}), $$

using this time the zero order approximate consistency property:

$$ \sum_{\sigma \in \mathcal{F}_K} \boldsymbol{\eta}_{K, \sigma} = \boldsymbol{\varepsilon}_K^0. $$

Notice that:

$$ - \sum_{K \in \mathcal{F}} \sum_{\sigma \in \mathcal{F}_K} (\mathcal{I}_K(\boldsymbol{\sigma}) \boldsymbol{\eta}_{K, \sigma}) \cdot \mathbf{v}_{\sigma} = - \sum_{K \in \mathcal{F}} \sum_{\sigma \in \mathcal{F}_K} \sum_{i=1}^{d} \sum_{j=1}^{d} \mathcal{I}_K(\sigma_{ij}) \boldsymbol{\eta}_{K, \sigma}^{j} v_{\sigma}^{i} = - \sum_{K \in \mathcal{F}} m_K \mathcal{I}_K(\boldsymbol{\sigma}) : \nabla \Pi_K^d(\mathbf{V}). $$

Consequently, we have using the definition of $\nabla_{\mathcal{F}}^d(\mathbf{V})$:

$$ \int_{\Omega} \mathrm{div}(\boldsymbol{\sigma}) \cdot \Pi_{\mathcal{F}}^{d}(\mathbf{V}) = \sum_{K \in \mathcal{F}} \int_{\mathcal{B}_K \cap \Omega} \zeta_K (\mathrm{div}(\boldsymbol{\sigma}) - \mathcal{D} \mathcal{I} \mathcal{V}_K(\mathcal{I}(\boldsymbol{\sigma}))) \cdot \mathcal{M}_K^d(\mathbf{V}) \\ + \sum_{K \in \mathcal{F}} \sum_{\sigma \in \mathcal{F}_K} ((\mathcal{I}_{\sigma}(\boldsymbol{\sigma}) - \mathcal{I}_K(\boldsymbol{\sigma})) \boldsymbol{\eta}_{K, \sigma}) \cdot (\mathcal{M}_K^d(\mathbf{V}) - \mathbf{v}_{\sigma}) \\ - \sum_{K \in \mathcal{F}} m_K \mathcal{I}_K(\boldsymbol{\sigma}) : \nabla \Pi_K^d(\mathbf{V}) + \sum_{\sigma \in \mathcal{F}_{\text{int}}} (\mathcal{I}_{\sigma}(\boldsymbol{\sigma}) \boldsymbol{\varepsilon}_{\sigma}) \cdot \mathbf{v}_{\sigma} \\ + \sum_{K \in \mathcal{F}} m_K (\mathcal{I}_K(\boldsymbol{\sigma}) \boldsymbol{\varepsilon}_K^0) \cdot \mathcal{M}_K^d(\mathbf{V}) \\ = \sum_{K \in \mathcal{F}} \int_{\mathcal{B}_K \cap \Omega} \zeta_K (\mathrm{div}(\boldsymbol{\sigma}) - \mathcal{D} \mathcal{I} \mathcal{V}_K(\mathcal{I}(\boldsymbol{\sigma}))) \cdot \mathcal{M}_K^d(\mathbf{V}) \\ + \sum_{K \in \mathcal{F}} \sum_{\sigma \in \mathcal{F}_K} ((\mathcal{I}_{\sigma}(\boldsymbol{\sigma}) - \mathcal{I}_K(\boldsymbol{\sigma})) \boldsymbol{\eta}_{K, \sigma}) \cdot (\mathcal{M}_K^d(\mathbf{V}) - \mathbf{v}_{\sigma}) \\ + \int_{\Omega} \sum_{K \in \mathcal{F}} \psi_K(\boldsymbol{\sigma} - \mathcal{I}_K(\boldsymbol{\sigma})) \cdot \nabla \Pi_K^d(\mathbf{V}) \\ - \int_{\Omega} \boldsymbol{\sigma} : \nabla_{\mathcal{F}}^d(\mathbf{V}) + \sum_{\sigma \in \mathcal{F}_{\text{int}}} (\mathcal{I}_{\sigma}(\boldsymbol{\sigma}) \boldsymbol{\varepsilon}_{\sigma}) \cdot \mathbf{v}_{\sigma} + \sum_{K \in \mathcal{F}} m_K (\mathcal{I}_K(\boldsymbol{\sigma}) \boldsymbol{\varepsilon}_K^0) \cdot \mathcal{M}_K^d(\mathbf{V}), $$
and thus:

$$
\begin{align*}
\mathcal{R}(\boldsymbol{\psi}_{\mathcal{F}}, \zeta_{\mathcal{F}}, \boldsymbol{\sigma}, \mathbf{V}) &= \sum_{K \in \mathcal{F}} \int_{\mathcal{B}_K \cap \Omega} \zeta_K (\mathrm{div}(\boldsymbol{\sigma}) - \mathcal{D} \mathcal{I} \mathcal{V}_K(\boldsymbol{\mathcal{I}}(\boldsymbol{\sigma}))) \cdot \mathcal{M}_K^d(\mathbf{V}) \\
&\quad + \sum_{K \in \mathcal{F}} \sum_{\sigma \in \mathcal{F}_K} ((\mathcal{I}_{\sigma}(\boldsymbol{\sigma}) - \mathcal{I}_K(\boldsymbol{\sigma}))\boldsymbol{\eta}_{K,\sigma}) \cdot (\mathcal{M}_K^d(\mathbf{V}) - \mathbf{v}_{\sigma}) \\
&\quad + \int_{\Omega} \sum_{K \in \mathcal{F}} \psi_K(\boldsymbol{\sigma} - \mathcal{I}_K(\boldsymbol{\sigma})) : \nabla \Pi_K^d(\mathbf{V}) \\
&\quad + \sum_{\sigma \in \mathcal{F}_{\mathrm{int}}} (\mathcal{I}_{\sigma}(\boldsymbol{\sigma})\boldsymbol{\varepsilon}_{\sigma}) \cdot \mathbf{v}_{\sigma} + \sum_{K \in \mathcal{F}} m_K(\mathcal{I}_K(\boldsymbol{\sigma})\boldsymbol{\varepsilon}_K^0) \cdot \mathcal{M}_K^d(\mathbf{V}).
\end{align*}
$$

We rewrite this last identity $\mathcal{R}(\boldsymbol{\psi}_{\mathcal{F}}, \zeta_{\mathcal{F}}, \boldsymbol{\sigma}, \mathbf{V}) = R_1 + R_2 + R_3 + R_4 + R_5$ with obvious notations. Using Cauchy–Schwarz inequality we get that:

$$
|R_1| \le C \left( \sum_{K \in \mathcal{F}} \int_{\mathcal{B}_K \cap \Omega} |\zeta_K| |\mathcal{D} \mathcal{I} \mathcal{V}_K(\boldsymbol{\mathcal{I}}(\boldsymbol{\sigma})) - \mathrm{div}(\boldsymbol{\sigma})|^2 \right)^{1/2} \|\mathbf{V}\|_0.
$$

Still using Cauchy–Schwarz inequality, we get

$$
\begin{align*}
|R_2| &\le \sum_{K \in \mathcal{F}} \sum_{\sigma \in \mathcal{F}_K} \frac{h_K |\eta_{K,\sigma}|}{m_K} m_K h_K^{-1} |\mathcal{I}_\sigma(\boldsymbol{\sigma}) - \mathcal{I}_K(\boldsymbol{\sigma})| |\mathcal{M}_K^d(\mathbf{V}) - \mathbf{v}_\sigma| \\
&\le \theta_\Pi \left( \sum_{K \in \mathcal{F}} \sum_{\sigma \in \mathcal{F}_K} m_K |\mathcal{I}_\sigma(\boldsymbol{\sigma}) - \mathcal{I}_K(\boldsymbol{\sigma})|^2 \right)^{1/2} \left( \sum_{K \in \mathcal{F}} \sum_{\sigma \in \mathcal{F}_K} m_K h_K^{-2} |\mathcal{M}_K^d(\mathbf{V}) - \mathbf{v}_\sigma|^2 \right)^{1/2} \\
&\le \theta_\Pi M_\psi^{1/2} \left( \sum_{K \in \mathcal{F}} \sum_{\sigma \in \mathcal{F}_K} \|\mathcal{I}_\sigma(\boldsymbol{\sigma}) - \boldsymbol{\sigma}\|_{L^2(\mathcal{B}_K \cap \Omega)^{d \times d}}^2 \right)^{1/2} \|\mathbf{V}\|_X.
\end{align*}
$$

Now, since

$$
|R_3| \le \left( \sum_{K \in \mathcal{F}} \int_{\mathcal{B}_K \cap \Omega} |\psi_K| |\boldsymbol{\sigma} - \mathcal{I}_K(\boldsymbol{\sigma})|^2 \right)^{1/2} \left( \sum_{K \in \mathcal{F}} \left( \frac{1}{m_K} \int_{\Omega} |\psi_K| \right) m_K \|\nabla \Pi_K^d(\mathbf{V})\|^2 \right)^{1/2},
$$

as $\Omega$ is Lipschitz and satisfies the cone condition with angle $\tau$ and radius $r$, defining $\delta > 0$ as the smallest real number such that for any $K \in \mathcal{F}$:

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

we easily get:

$$
\frac{1}{m_K} \int_{\Omega} |\psi_K| \le M_\psi \theta_{\mathcal{F}} \frac{S_1^d \kappa_\psi^d \delta^d}{|C(0, \tau, 1)|} \quad \text{as} \quad \frac{|\mathcal{B}_K \cap \Omega|}{m_K} \le \theta_{\mathcal{F}} \frac{|\mathcal{B}_K|}{|B_K \cap \Omega|} \le \theta_{\mathcal{F}} \frac{S_1^d \kappa_\psi^d \delta^d}{|C(0, \tau, 1)|}, \quad (46)
$$

with $S_1^d = |B(0, 1)|$ the volume of the unit ball, from which we deduce:

$$
|R_3| \le C_{\mathcal{V}} \left( \frac{M_{\psi}^2 \theta_{\mathcal{F}} S_1^d \kappa_{\psi}^d \delta^d}{|C(0, \tau, 1)|} \right)^{1/2} \left( \sum_{K \in \mathcal{F}} \| \mathcal{I}_K(\boldsymbol{\sigma}) - \boldsymbol{\sigma} \|_{L^2(\mathcal{B}_K \cap \Omega)^{d \times d}}^2 \right)^{1/2} \| \mathbf{V} \|_X.
$$

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

$$
|R_4| = \left| \sum_{\sigma \in \mathcal{F}} \sum_{K \in \mathcal{F}_\sigma} \frac{(\mathcal{I}_\sigma(\boldsymbol{\sigma})\boldsymbol{\varepsilon}_\sigma) \cdot \mathbf{v}_\sigma}{\mathrm{card}(\mathcal{F}_\sigma)} \right| \\
\le \theta_{\mathcal{A}} h^p \left( \sum_{K \in \mathcal{F}} \sum_{\sigma \in \mathcal{F}_K} \frac{m_K}{\mathrm{card}(\mathcal{F}_\sigma)} |\mathcal{I}_\sigma(\boldsymbol{\sigma})| |\mathbf{v}_\sigma - \mathcal{M}_K^d(\mathbf{V})| + \sum_{K \in \mathcal{F}} m_K |\mathcal{M}_K^d(\mathbf{V})| \left( \sum_{\sigma \in \mathcal{F}_K} \frac{1}{\mathrm{card}(\mathcal{F}_\sigma)} |\mathcal{I}_\sigma(\boldsymbol{\sigma})| \right) \right),
$$

and using Cauchy–Schwarz inequality, this leads to:

$$
|R_4| \le \theta_{\mathcal{A}} h^p (h |\mathbf{V}|_X + |\mathbf{V}|_0) \left( \sum_{K \in \mathcal{F}} \sum_{\sigma \in \mathcal{F}_K} m_K |\mathcal{I}_\sigma(\boldsymbol{\sigma})|^2 \right)^{1/2} .
$$
and thus:

$$
\begin{align*}
\mathcal{R}(\boldsymbol{\psi}_{\mathcal{F}}, \zeta_{\mathcal{F}}, \boldsymbol{\sigma}, \mathbf{V}) ={}& \sum_{K \in \mathcal{F}} \int_{\mathcal{B}_K \cap \Omega} \zeta_K (\mathrm{div}(\boldsymbol{\sigma}) - \mathcal{D} \mathcal{I} \mathcal{V}_K(\boldsymbol{\mathcal{I}}(\boldsymbol{\sigma}))) \cdot \mathcal{M}_K^d(\mathbf{V}) \\
& + \sum_{K \in \mathcal{F}} \sum_{\sigma \in \mathcal{F}_K} ((\mathcal{I}_{\sigma}(\boldsymbol{\sigma}) - \mathcal{I}_K(\boldsymbol{\sigma})) \eta_{K,\sigma}) \cdot (\mathcal{M}_K^d(\mathbf{V}) - \mathbf{v}_{\sigma}) \\
& + \int_{\Omega} \sum_{K \in \mathcal{F}} \psi_K(\boldsymbol{\sigma} - \mathcal{I}_K(\boldsymbol{\sigma})) : \nabla \Pi_K^d(\mathbf{V}) \\
& + \sum_{\sigma \in \mathcal{F}_{\mathrm{int}}} (\mathcal{I}_{\sigma}(\boldsymbol{\sigma}) \boldsymbol{\varepsilon}_{\sigma}) \cdot \mathbf{v}_{\sigma} + \sum_{K \in \mathcal{F}} m_K (\mathcal{I}_K(\boldsymbol{\sigma}) \boldsymbol{\varepsilon}_K^0) \cdot \mathcal{M}_K^d(\mathbf{V}).
\end{align*}
$$

We rewrite this last identity $\mathcal{R}(\boldsymbol{\psi}_{\mathcal{F}}, \zeta_{\mathcal{F}}, \boldsymbol{\sigma}, \mathbf{V}) = R_1 + R_2 + R_3 + R_4 + R_5$ with obvious notations. Using Cauchy–Schwarz inequality we get that:

$$
|R_1| \le C \left( \sum_{K \in \mathcal{F}} \int_{\mathcal{B}_K \cap \Omega} |\zeta_K| |\mathcal{D} \mathcal{I} \mathcal{V}_K(\boldsymbol{\mathcal{I}}(\boldsymbol{\sigma})) - \mathrm{div}(\boldsymbol{\sigma})|^2 \right)^{1/2} \|\mathbf{V}\|_0.
$$

Still using Cauchy–Schwarz inequality, we get

$$
\begin{align*}
|R_2| &\le \sum_{K \in \mathcal{F}} \sum_{\sigma \in \mathcal{F}_K} \frac{h_K |\eta_{K,\sigma}|}{m_K} m_K h_K^{-1} |\mathcal{I}_\sigma(\boldsymbol{\sigma}) - \mathcal{I}_K(\boldsymbol{\sigma})| |\mathcal{M}_K^d(\mathbf{V}) - \mathbf{v}_\sigma| \\
&\le \theta_\Pi \left( \sum_{K \in \mathcal{F}} \sum_{\sigma \in \mathcal{F}_K} m_K |\mathcal{I}_\sigma(\boldsymbol{\sigma}) - \mathcal{I}_K(\boldsymbol{\sigma})|^2 \right)^{1/2} \left( \sum_{K \in \mathcal{F}} \sum_{\sigma \in \mathcal{F}_K} m_K h_K^{-2} |\mathcal{M}_K^d(\mathbf{V}) - \mathbf{v}_\sigma|^2 \right)^{1/2} \\
&\le \theta_\Pi M_\psi^{1/2} \left( \sum_{K \in \mathcal{F}} \sum_{\sigma \in \mathcal{F}_K} \|\mathcal{I}_\sigma(\boldsymbol{\sigma}) - \boldsymbol{\sigma}\|_{L^2(\mathcal{B}_K \cap \Omega)^{d \times d}}^2 \right)^{1/2} \|\mathbf{V}\|_X.
\end{align*}
$$

Now, since

$$
|R_3| \le \left( \sum_{K \in \mathcal{F}} \int_{\mathcal{B}_K \cap \Omega} |\psi_K| |\boldsymbol{\sigma} - \mathcal{I}_K(\boldsymbol{\sigma})|^2 \right)^{1/2} \left( \sum_{K \in \mathcal{F}} \left( \frac{1}{m_K} \int_{\Omega} |\psi_K| \right) m_K \|\nabla \Pi_K^d(\mathbf{V})\|^2 \right)^{1/2},
$$

as $\Omega$ is Lipschitz and satisfies the cone condition with angle $\tau$ and radius $r$, defining $\delta > 0$ as the smallest real number such that for any $K \in \mathcal{F}$:

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

we easily get:

$$
\frac{1}{m_K} \int_{\Omega} |\psi_K| \le M_\psi \theta_{\mathcal{F}} \frac{S_1^d \kappa_\psi^d \delta^d}{|C(0, \tau, 1)|} \quad \text{as} \quad \frac{|\mathcal{B}_K \cap \Omega|}{m_K} \le \theta_{\mathcal{F}} \frac{|\mathcal{B}_K|}{|B_K \cap \Omega|} \le \theta_{\mathcal{F}} \frac{S_1^d \kappa_\psi^d \delta^d}{|C(0, \tau, 1)|}, \quad (46)
$$

with $S_1^d = |B(0, 1)|$ the volume of the unit ball, from which we deduce:

$$
|R_3| \le C_V \left( \frac{M_\psi^2 \theta_{\mathcal{F}} S_1^d \kappa_\psi^d \delta^d}{|C(0, \tau, 1)|} \right)^{1/2} \left( \sum_{K \in \mathcal{F}} \| \mathcal{I}_K(\boldsymbol{\sigma}) - \boldsymbol{\sigma} \|_{L^2(\mathcal{B}_K \cap \Omega)^{d \times d}}^2 \right)^{1/2} \| \mathbf{V} \|_X.
$$

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

$$
|R_4| = \left| \sum_{\sigma \in \mathcal{F}} \sum_{K \in \mathcal{F}_\sigma} \frac{(\mathcal{I}_\sigma(\boldsymbol{\sigma}) \boldsymbol{\varepsilon}_\sigma) \cdot \mathbf{v}_\sigma}{\mathrm{card}(\mathcal{F}_\sigma)} \right| \\
\le \theta_{\mathcal{A}} h^p \left( \sum_{K \in \mathcal{F}} \sum_{\sigma \in \mathcal{F}_K} \frac{m_K}{\mathrm{card}(\mathcal{F}_\sigma)} |\mathcal{I}_\sigma(\boldsymbol{\sigma})| |\mathbf{v}_\sigma - \mathcal{M}_K^d(\mathbf{V})| + \sum_{K \in \mathcal{F}} m_K |\mathcal{M}_K^d(\mathbf{V})| \left( \sum_{\sigma \in \mathcal{F}_K} \frac{1}{\mathrm{card}(\mathcal{F}_\sigma)} |\mathcal{I}_\sigma(\boldsymbol{\sigma})| \right) \right),
$$

and using Cauchy–Schwarz inequality, this leads to:

$$
|R_4| \le \theta_{\mathcal{A}} h^p (h |\mathbf{V}|_X + |\mathbf{V}|_0) \left( \sum_{K \in \mathcal{F}} \sum_{\sigma \in \mathcal{F}_K} m_K |\mathcal{I}_\sigma(\boldsymbol{\sigma})|^2 \right)^{1/2} .
$$
Then, notice that by definition of $\mathcal{B}_K$, we have:

$$|\mathcal{I}_{\sigma}(\boldsymbol{\sigma})|^2 \le \frac{|\mathcal{B}_K|}{|B_{\sigma}||\mathcal{B}_K|} \int_{\mathcal{B}_K} |E(\boldsymbol{\sigma})|^2 \le \frac{\theta_{\mathcal{F}}^d}{|\mathcal{B}_K|} \int_{\mathcal{B}_K} |E(\boldsymbol{\sigma})|^2.$$

As:

$$\frac{m_K}{|\mathcal{B}_K|} = \frac{m_K}{|B_K \cap \Omega|} \frac{|B_K \cap \Omega|}{|\mathcal{B}_K|} \le \theta_{\mathcal{F}} \frac{|B_K|}{|\mathcal{B}_K|} \le \theta_{\mathcal{F}},$$

we consequently get:

$$|R_4| \le \eta_{\psi}^{1/2} \theta_{\mathcal{F}}^{1/2} \theta_{\mathcal{F}}^{1/2} \theta_{\mathcal{F}}^{d/2} \theta_{\mathcal{A}} h^p (h|\mathbf{V}|_X + |\mathbf{V}|_0) C_{E,0} \|\boldsymbol{\sigma}\|_{L^2(\Omega)^{d \times d}}.$$

Finally, we have using Cauchy-Schwarz inequality that:

$$|R_5| \le \theta_{\mathcal{A}} h^p |\mathbf{V}|_0 \left( \sum_{K \in \mathcal{F}} m_K |\mathcal{I}_K(\boldsymbol{\sigma})|^2 \right)^{1/2} \le \theta_{\mathcal{F}}^{1/2} \eta_{\psi}^{d/2} \theta_{\mathcal{A}} h^p |\mathbf{V}|_0 \|\boldsymbol{\sigma}\|_{L^2(\Omega)^{d \times d}},$$

as $(m_K/|\mathcal{B}_K \cap \Omega|) \le \theta_{\mathcal{F}}(|B_K \cap \Omega|/|\mathcal{B}_K \cap \Omega|) \le \theta_{\mathcal{F}}$. Gathering the estimates on the $R_i$'s, we do obtain that there exists $C > 0$ such that $|\mathcal{R}(\psi_{\mathcal{F}}, \zeta_{\mathcal{F}}, \boldsymbol{\sigma}, \mathbf{V})| \le C \mathcal{E}_h^{\mathcal{R}}(\zeta_{\mathcal{F}}, \boldsymbol{\sigma}) \|\mathbf{V}\|_X$. $\square$

Proceeding exactly as in the scalar case (see [30]), one could easily establish that for any $\mathbf{v} \in H^1(\Omega)^d$:

$$\|\mathbf{v} - \mathcal{M}_K^d (\mathcal{I}_K(\mathbf{v}))\|_{L^2(\mathcal{B}_K \cap \Omega)^d} \le C h_K |E\mathbf{v}|_{H^1(\mathcal{B}_K)^d}, \quad (47)$$

for any $\mathbf{v} \in H^2(\Omega)^d$:

$$\|\nabla \mathbf{v} - \nabla_K^d (\mathcal{I}_K(\mathbf{v}))\|_{L^2(\mathcal{B}_K \cap \Omega)^{d \times d}} \le C(h_K + h_K^p) \|E\mathbf{v}\|_{H^2(\mathcal{B}_K)^d}, \quad (48)$$

and

$$\|\mathbf{v} - \Pi_K^d (\mathcal{I}_K(\mathbf{v}))\|_{L^2(\mathcal{B}_K \cap \Omega)^d} \le C(h_K^2 + h_K^{p+1}) \|E\mathbf{v}\|_{H^2(\mathcal{B}_K)^d}, \quad (49)$$

for any $\mathbf{v} \in H^2(\Omega)^{d \times d}$

$$\|\mathrm{div} \, \mathbf{v} - \mathcal{D} \mathcal{I} \mathcal{V}_K^d (\mathcal{I}_K(\mathbf{v}))\|_{L^2(\mathcal{B}_K \cap \Omega)^d} \le C(h_K + h_K^p) \|E\mathbf{v}\|_{H^2(\mathcal{B}_K)^{d \times d}}, \quad (50)$$

for any $\mathbf{v} \in H^1(\Omega)^d$

$$\|\mathbf{v} - \Pi_{\mathcal{F}}^d (\mathcal{I}(\mathbf{v}))\|_{L^2(\Omega)^d} \le C_{E,1} C(h+h^p) \|\mathbf{v}\|_{H^1(\Omega)^d}, \quad (51)$$

for any $\mathbf{v} \in H^2(\Omega)^d$

$$\|\nabla \mathbf{v} - \nabla_{\mathcal{F}}^d (\mathcal{I}(\mathbf{v}))\|_{L^2(\Omega)^{d \times d}} \le C_{E,2} C(h^2 + h^p) \|\mathbf{v}\|_{H^2(\Omega)^d}, \quad (52)$$

and

$$\|\mathbf{v} - \Pi_{\mathcal{N}}^d (\mathcal{I}(\mathbf{v}))\|_{L^2(\Omega)^d} \le C_{E,2} C(h^2 + h^{p+1}) \|\mathbf{v}\|_{H^2(\Omega)^d}, \quad (53)$$

where the constants $C > 0$ in the above result can vary from line to line but only depend on the quality parameters $\theta_{\mathcal{A}}$, $\theta_{\mathcal{F}}$, $\theta_{\mathcal{A},M}$, $\theta_{\Pi}$, $\eta$, $\eta_{\psi}$, $M_{\psi}$, $C_{\psi}$, $\theta_{\mathcal{F}}$ and not on $h$. The same results hold when replacing $\mathcal{I}$ by $\mathcal{I}^0$ if $\mathbf{v} \in H_0^1(\Omega)^d$. This immediately implies that the interpolation terms of estimates (39) and (40) will have the expected behavior regarding the regularity of the exact solution (provided $p \ge 1$), leading to the same convergence order as mesh based methods. For the extra term $\mathcal{E}_h^{\mathcal{S}}(\zeta_{\mathcal{F}}, \mathbf{u}, \mathbf{f})$ specific to the NEM, the regularity assumption leading to an optimal error estimate is more subtle. The first thing to notice is that for $\mathbf{f}$ we need more regularity than $L^2(\Omega)^d$, which is in fact the same as for approximating the right-hand side in the case of mesh-based methods. Thus assuming that $\mathbf{f}$ has enough regularity (say $\mathbf{f} \in H^1$) we see that thanks to the weights $\zeta_{\mathcal{F}}$ we can cancel the extra error where $\mathbf{u}$ does not have the required regularity (say $\mathbf{u} \in H^3$), for instance near the boundaries or where $\mathbf{f}$ is not regular enough. However, if the global order of convergence is preserved by doing so should $\zeta_{\mathcal{F}}$ cancel in too many places this would increase the constant $M_{\zeta}$ and consequently the overall constant involved in estimates (39) and (40).

Finally, let us mention that if the constants involved in the above estimates as well as the extra term $\mathcal{E}_h^{\mathcal{S}}(\zeta_{\mathcal{F}}, \mathbf{u}, \mathbf{f})$ depend on the quadrature and weight families $\psi_{\mathcal{F}}$ and $\zeta_{\mathcal{F}}$, they hold for any
Then, notice that by definition of $\mathcal{B}_K$, we have:

$$|\mathcal{I}_{\sigma}(\boldsymbol{\sigma})|^2 \le \frac{|\mathcal{B}_K|}{|B_{\sigma}||\mathcal{B}_K|} \int_{\mathcal{B}_K} |E(\boldsymbol{\sigma})|^2 \le \frac{\theta_{\mathcal{I}}^d}{|\mathcal{B}_K|} \int_{\mathcal{B}_K} |E(\boldsymbol{\sigma})|^2.$$

As:

$$\frac{m_K}{|\mathcal{B}_K|} = \frac{m_K}{|B_K \cap \Omega|} \frac{|B_K \cap \Omega|}{|\mathcal{B}_K|} \le \theta_{\mathcal{F}} \frac{|B_K|}{|\mathcal{B}_K|} \le \theta_{\mathcal{F}},$$

we consequently get:

$$|R_4| \le \eta_{\psi}^{1/2} \theta_{\mathcal{F}}^{1/2} \theta_{\mathcal{F}}^{1/2} \theta_{\mathcal{F}}^{d/2} \theta_{\mathcal{A}} h^p (h|\mathbf{V}|_X + |\mathbf{V}|_0) C_{E,0} \|\boldsymbol{\sigma}\|_{L^2(\Omega)^{d \times d}}.$$

Finally, we have using Cauchy–Schwarz inequality that:

$$|R_5| \le \theta_{\mathcal{A}} h^p |\mathbf{V}|_0 \left( \sum_{K \in \mathcal{F}} m_K |\mathcal{I}_K(\boldsymbol{\sigma})|^2 \right)^{1/2} \le \theta_{\mathcal{F}}^{1/2} \eta_{\psi}^{d/2} \theta_{\mathcal{A}} h^p |\mathbf{V}|_0 \|\boldsymbol{\sigma}\|_{L^2(\Omega)^{d \times d}},$$

as $(m_K/|\mathcal{B}_K \cap \Omega|) \le \theta_{\mathcal{F}}(|B_K \cap \Omega|/|\mathcal{B}_K \cap \Omega|) \le \theta_{\mathcal{F}}$. Gathering the estimates on the $R_i$'s, we do obtain that there exists $C > 0$ such that $|\mathcal{R}(\psi_{\mathcal{F}}, \zeta_{\mathcal{F}}, \boldsymbol{\sigma}, \mathbf{V})| \le C \mathcal{E}_h^{\mathcal{R}}(\zeta_{\mathcal{F}}, \boldsymbol{\sigma}) \|\mathbf{V}\|_X$. $\square$

Proceeding exactly as in the scalar case (see [30]), one could easily establish that for any $\mathbf{v} \in H^1(\Omega)^d$:

$$\|\mathbf{v} - \mathcal{M}_K^d(\mathcal{I}_K(\mathbf{v}))\|_{L^2(\mathcal{B}_K \cap \Omega)^d} \le C h_K |E\mathbf{v}|_{H^1(\mathcal{B}_K)^d}, \quad (47)$$

for any $\mathbf{v} \in H^2(\Omega)^d$:

$$\|\nabla \mathbf{v} - \nabla_K^d (\mathcal{I}_K(\mathbf{v}))\|_{L^2(\mathcal{B}_K \cap \Omega)^{d \times d}} \le C(h_K + h_K^p) \|E\mathbf{v}\|_{H^2(\mathcal{B}_K)^d}, \quad (48)$$

and

$$\|\mathbf{v} - \Pi_K^d (\mathcal{I}_K(\mathbf{v}))\|_{L^2(\mathcal{B}_K \cap \Omega)^d} \le C(h_K^2 + h_K^{p+1}) \|E\mathbf{v}\|_{H^2(\mathcal{B}_K)^d}, \quad (49)$$

for any $\mathbf{v} \in H^2(\Omega)^{d \times d}$

$$\|\mathrm{div} \, \mathbf{v} - \mathcal{D} \mathcal{I} \mathcal{V}_K^d (\mathcal{I}_K(\mathbf{v}))\|_{L^2(\mathcal{B}_K \cap \Omega)^d} \le C(h_K + h_K^p) \|E\mathbf{v}\|_{H^2(\mathcal{B}_K)^{d \times d}}, \quad (50)$$

for any $\mathbf{v} \in H^1(\Omega)^d$

$$\|\mathbf{v} - \Pi_{\mathcal{F}}^d (\mathcal{I}(\mathbf{v}))\|_{L^2(\Omega)^d} \le C_{E,1} C(h+h^p) \|\mathbf{v}\|_{H^1(\Omega)^d}, \quad (51)$$

for any $\mathbf{v} \in H^2(\Omega)^d$

$$\|\nabla \mathbf{v} - \nabla_{\mathcal{F}}^d (\mathcal{I}(\mathbf{v}))\|_{L^2(\Omega)^{d \times d}} \le C_{E,2} C(h^2 + h^p) \|\mathbf{v}\|_{H^2(\Omega)^d}, \quad (52)$$

and

$$\|\mathbf{v} - \Pi_{\mathcal{N}}^d (\mathcal{I}(\mathbf{v}))\|_{L^2(\Omega)^d} \le C_{E,2} C(h^2 + h^{p+1}) \|\mathbf{v}\|_{H^2(\Omega)^d}, \quad (53)$$

where the constants $C > 0$ in the above result can vary from line to line but only depend on the quality parameters $\theta_{\mathcal{A}}$, $\theta_{\mathcal{F}}$, $\theta_{\mathcal{A},M}$, $\theta_{\Pi}$, $\eta$, $\eta_{\psi}$, $M_{\psi}$, $C_{\psi}$, $\theta_{\mathcal{F}}$ and not on $h$. The same results hold when replacing $\mathcal{I}$ by $\mathcal{I}^0$ if $\mathbf{v} \in H_0^1(\Omega)^d$. This immediately implies that the interpolation terms of estimates (39) and (40) will have the expected behavior regarding the regularity of the exact solution (provided $p \ge 1$), leading to the same convergence order as mesh based methods. For the extra term $\mathcal{E}_h^{\mathcal{S}}(\zeta_{\mathcal{F}}, \mathbf{u}, \mathbf{f})$ specific to the NEM, the regularity assumption leading to an optimal error estimate is more subtle. The first thing to notice is that for $\mathbf{f}$ we need more regularity than $L^2(\Omega)^d$, which is in fact the same as for approximating the right-hand side in the case of mesh-based methods. Thus assuming that $\mathbf{f}$ has enough regularity (say $\mathbf{f} \in H^1$) we see that thanks to the weights $\zeta_{\mathcal{F}}$ we can cancel the extra error where $\mathbf{u}$ does not have the required regularity (say $\mathbf{u} \in H^3$), for instance near the boundaries or where $\mathbf{f}$ is not regular enough. However, if the global order of convergence is preserved by doing so should $\zeta_{\mathcal{F}}$ cancel in too many places this would increase the constant $M_{\zeta}$ and consequently the overall constant involved in estimates (39) and (40).

Finally, let us mention that if the constants involved in the above estimates as well as the extra term $\mathcal{E}_h^{\mathcal{S}}(\zeta_{\mathcal{F}}, \mathbf{u}, \mathbf{f})$ depend on the quadrature and weight families $\psi_{\mathcal{F}}$ and $\zeta_{\mathcal{F}}$, they hold for any
pair $(\psi_{\mathcal{F}}, \zeta_{\mathcal{F}})$. Such pairs are mainly a theoretical ingredient useful to bridge the gap between $X_{\mathcal{N}}^d$ and $L^2(\Omega)^d$: in practice they are never computed and not even chosen. Thus, we could had a minimum over all the pairs $(\psi_{\mathcal{F}}, \zeta_{\mathcal{F}})$ in estimates (39) and (40). If necessary, the constants corresponding to $\psi_{\mathcal{F}}$ can be estimated using the derivation of $\psi_{\mathcal{F}}$ given in [30] that is used there to establish the existence of quadrature families (and which can be adapted to $\zeta_{\mathcal{F}}$). This derivation directly relates $\eta_{\psi}$, $M_{\psi}$ and $C_{\psi}$ to the other quality parameters.

# 6. Numerical exploration

## 6.1. Standard materials

From [36] on $\Omega = [0, 1]^2$ we consider the sinusoidal analytic solution:

$$u^i = \prod_{k=1}^{d} \sin(4\pi x^k) \quad \text{for all } i = 1 \dots d \quad \text{with } \lambda = \mu = 1,$$

with right-hand side:

$$f = -\operatorname{div} \sigma(\mathbf{u}) = \sum_{i=1}^{d} 64\pi^2 \prod_{k=1}^{d} \sin(4\pi x^k) e_i - \sum_{i=1}^{d} \sum_{j \neq i} 32\pi^2 \cos(4\pi x^j) \cos(4\pi x^i) e_j,$$

as well as the exponential analytic solution from [34]:

$$u^i = \exp \left( \cos \left( \sum_{j=1}^{d} C_{ij} x^j \right) \right) \quad \text{with } \lambda = \mu = 1 \quad \text{and} \quad C = \begin{pmatrix} 1 & 1 \\ 2 & -1 \end{pmatrix},$$

and second member:

$$\begin{aligned} f &= -\operatorname{div}\sigma(\mathbf{u}) \\ &= \sum_{j=1}^{d} \left( \sum_{i=1}^{d} 2C_{ij}C_{ii} \left( \cos\left(\sum_{m=1}^{d} C_{im}x^m\right) - \sin^2\left(\sum_{m=1}^{d} C_{im}x^m\right) \right) \exp\left(\cos\left(\sum_{m=1}^{d} C_{im}x^m\right)\right) \right) e_j \\ &\quad + \sum_{j=1}^{d} \left( \sum_{i=1}^{d} C_{ji}^2 \left( \cos\left(\sum_{m=1}^{d} C_{jm}x^m\right) - \sin^2\left(\sum_{m=1}^{d} C_{jm}x^m\right) \right) \exp\left(\cos\left(\sum_{m=1}^{d} C_{jm}x^m\right)\right) \right) e_j \end{aligned}$$

The discretization networks for domain $\Omega$ are generated using the method of [31] that relies on the node generator of [52], using the interface enrichment technique described in [31] to ensure that hypothesis (H5) is satisfied. The obtained cell nodes are displayed on the left side of Figure 1. On the right side of this figure, we also display a distorted version of the same network as described in [31], that we will use to assess the robustness of the NEM for linear elasticity with respect to network distortion. The network geometries are of course generated by solving the linear system introduced in [31], and we set $p = 1$.

Convergence curves for the sinusoidal and exponential test cases are displayed on Figures 2 and 3, while approximate convergence orders computed from those curves are gathered in Table 1, for both the basic NEM and its flux based counterpart (CNEM). For comparison purposes, on the same figures we display convergence curves and orders for the first order virtual element method (VEM) [33] and its flux based counterpart (CVEM) computed as explained in Section 4.2.2 using the approach of [47], both applied on a Delaunay mesh sequence of the unit square. Notice that on this simplicial mesh, the virtual element method degenerates into a $\mathbb{P}^1$ finite element scheme. We see on those results that the NEM and the flux based CNEM perform particularly well, leading to optimal convergence rates.
pair $(\psi_{\mathcal{F}}, \zeta_{\mathcal{F}})$. Such pairs are mainly a theoretical ingredient useful to bridge the gap between $X_{\mathcal{N}}^d$ and $L^2(\Omega)^d$: in practice they are never computed and not even chosen. Thus, we could had a minimum over all the pairs $(\psi_{\mathcal{F}}, \zeta_{\mathcal{F}})$ in estimates (39) and (40). If necessary, the constants corresponding to $\psi_{\mathcal{F}}$ can be estimated using the derivation of $\psi_{\mathcal{F}}$ given in [30] that is used there to establish the existence of quadrature families (and which can be adapted to $\zeta_{\mathcal{F}}$). This derivation directly relates $\eta_{\psi}$, $M_{\psi}$ and $C_{\psi}$ to the other quality parameters.

# 6. Numerical exploration

## 6.1. Standard materials

From [36] on $\Omega = [0, 1]^2$ we consider the sinusoidal analytic solution:

$$u^i = \prod_{k=1}^{d} \sin(4\pi x^k) \quad \text{for all } i = 1 \dots d \quad \text{with } \lambda = \mu = 1,$$

with right-hand side:

$$\mathbf{f} = -\operatorname{div} \sigma(\mathbf{u}) = \sum_{i=1}^{d} 64\pi^2 \prod_{k=1}^{d} \sin(4\pi x^k) e_i - \sum_{i=1}^{d} \sum_{j \neq i} 32\pi^2 \cos(4\pi x^j) \cos(4\pi x^i) e_j,$$

as well as the exponential analytic solution from [34]:

$$u^i = \exp \left( \cos \left( \sum_{j=1}^{d} C_{ij} x^j \right) \right) \quad \text{with } \lambda = \mu = 1 \quad \text{and} \quad C = \begin{pmatrix} 1 & 1 \\ 2 & -1 \end{pmatrix},$$

and second member:

$$\begin{aligned} \mathbf{f} = -\operatorname{div} \sigma(\mathbf{u}) &= \sum_{j=1}^{d} \left( \sum_{i=1}^{d} 2C_{ij}C_{ii} \left( \cos\left(\sum_{m=1}^{d} C_{im}x^m\right) - \sin^2\left(\sum_{m=1}^{d} C_{im}x^m\right) \right) \exp\left(\cos\left(\sum_{m=1}^{d} C_{im}x^m\right)\right) \right) e_j \\ &\quad + \sum_{j=1}^{d} \left( \sum_{i=1}^{d} C_{ji}^2 \left( \cos\left(\sum_{m=1}^{d} C_{jm}x^m\right) - \sin^2\left(\sum_{m=1}^{d} C_{jm}x^m\right) \right) \exp\left(\cos\left(\sum_{m=1}^{d} C_{jm}x^m\right)\right) \right) e_j \end{aligned}$$

The discretization networks for domain $\Omega$ are generated using the method of [31] that relies on the node generator of [52], using the interface enrichment technique described in [31] to ensure that hypothesis (H5) is satisfied. The obtained cell nodes are displayed on the left side of Figure 1. On the right side of this figure, we also display a distorted version of the same network as described in [31], that we will use to assess the robustness of the NEM for linear elasticity with respect to network distortion. The network geometries are of course generated by solving the linear system introduced in [31], and we set $p = 1$.

Convergence curves for the sinusoidal and exponential test cases are displayed on Figures 2 and 3, while approximate convergence orders computed from those curves are gathered in Table 1, for both the basic NEM and its flux based counterpart (CNEM). For comparison purposes, on the same figures we display convergence curves and orders for the first order virtual element method (VEM) [33] and its flux based counterpart (CVEM) computed as explained in Section 4.2.2 using the approach of [47], both applied on a Delaunay mesh sequence of the unit square. Notice that on this simplicial mesh, the virtual element method degenerates into a $\mathbb{P}^1$ finite element scheme. We see on those results that the NEM and the flux based CNEM perform particularly well, leading to optimal convergence rates.
<figure><img src="image_4.png" /><figcaption>**Figure 1.** Cell nodes in dimension 2: left reference cloud, right distorted cloud (reproduced from [31]).</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 2.** Convergence curves for the sinusoidal test case.</figcaption></figure>

<table><caption>**Table 1.** Approximate orders of convergence</caption><thead><tr><th></th><th>Sinusoidal</th><th>Exponential</th></tr></thead><tbody><tr><td>NEM</td><td>2.07</td><td>2.27</td></tr><tr><td>CNEM</td><td>2.31</td><td>2.30</td></tr><tr><td>NEM distorted</td><td>2.11</td><td>2.08</td></tr><tr><td>CNEM distorted</td><td>2.41</td><td>2.24</td></tr><tr><td>VEM</td><td>2.00</td><td>2.01</td></tr><tr><td>CVEM</td><td>2.08</td><td>2.09</td></tr></tbody></table>
<figure><img src="image_4.png" /><figcaption>**Figure 1.** Cell nodes in dimension 2: left reference cloud, right distorted cloud (reproduced from [31]).</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 2.** Convergence curves for the sinusoidal test case.</figcaption></figure>

<table><caption>**Table 1.** Approximate orders of convergence</caption><thead><tr><th></th><th>Sinusoidal</th><th>Exponential</th></tr></thead><tbody><tr><td>NEM</td><td>2.07</td><td>2.27</td></tr><tr><td>CNEM</td><td>2.31</td><td>2.30</td></tr><tr><td>NEM distorted</td><td>2.11</td><td>2.08</td></tr><tr><td>CNEM distorted</td><td>2.41</td><td>2.24</td></tr><tr><td>VEM</td><td>2.00</td><td>2.01</td></tr><tr><td>CVEM</td><td>2.08</td><td>2.09</td></tr></tbody></table>
<figure><img src="image_6.png" /><figcaption>Figure 3. Convergence curves for the exponential test case.</figcaption></figure>

## 6.2. Quasi-incompressible materials

The aim of this subsection is to assess the robustness of the NEM with respect to quasi-incompressibility. To this end, we consider the following test case from [36] on $\Omega = [0, 1]^2$ with $v \in \{0.25, 0.4995, 0.4999995\}$ and

$$ \mu = 1 \quad \text{and} \quad v = \frac{\lambda}{2(\lambda + \mu)}, $$

$$ u^x = \cos\left(\frac{2\pi}{\lambda}x\right)\sin(2\pi y) \quad \text{and} \quad u^y = \sin(2\pi x)\cos\left(\frac{2\pi}{\lambda}y\right), $$

with second member:

$$ f^x = 4\pi^2 \left( \mu + \frac{1}{\lambda} + \frac{2\mu}{\lambda^2} \right) \cos\left(\frac{2\pi}{\lambda}x\right) \sin(2\pi y) + 4\pi^2 \left( 1 + \frac{\mu}{\lambda} \right) \cos(2\pi x) \sin\left(\frac{2\pi}{\lambda}y\right), $$

$$ f^y = 4\pi^2 \left( \mu + \frac{1}{\lambda} + \frac{2\mu}{\lambda^2} \right) \sin(2\pi x) \cos\left(\frac{2\pi}{\lambda}y\right) + 4\pi^2 \left( 1 + \frac{\mu}{\lambda} \right) \sin\left(\frac{2\pi}{\lambda}x\right) \cos(2\pi y). $$

Similarly to the observations of [36], from Figure 4 we recover the fact that the $\mathbb{P}^1$ finite element (i.e. the first order VEM scheme) suffers when $\lambda$ grows. For $\lambda \approx 10 \times 10^6$ ($v = 0.4999995$) numerical locking is observed, true convergence being only recovered for very small $h$. The NEM and CNEM outperforms the VEM scheme in this situation on the regular point cloud, however if convergence is maintained on the distorted point cloud the initial error on coarse networks is much larger than even for the locking VEM scheme. In the moderate regime $\lambda \approx 10 \times 10^3$ ($v = 0.4995$), the NEM and CNEM clearly outperforms the VEM scheme, even on the distorted point cloud, while in the neutral regime $\lambda = 1$ ($v = 0.25$) we recover results similar to the previous analytic test cases. Notice that the network generation algorithm from [31] produces well connected networks, with each cell node connected to all its immediate neighbors. Thus we have in general $\text{card}(\mathcal{F}_K) > d + 1$ which probably explains why the NEM avoids numerical locking: there is enough degrees of freedom per cell to represent accurately divergence free displacements.
<figure><img src="image_6.png" /><figcaption>Figure 3. Convergence curves for the exponential test case.</figcaption></figure>

## 6.2. Quasi-incompressible materials

The aim of this subsection is to assess the robustness of the NEM with respect to quasi-incompressibility. To this end, we consider the following test case from [36] on $\Omega = [0, 1]^2$ with $v \in \{0.25, 0.4995, 0.4999995\}$ and

$$ \mu = 1 \quad \text{and} \quad v = \frac{\lambda}{2(\lambda + \mu)}, $$

$$ u^x = \cos\left(\frac{2\pi}{\lambda}x\right)\sin(2\pi y) \quad \text{and} \quad u^y = \sin(2\pi x)\cos\left(\frac{2\pi}{\lambda}y\right), $$

with second member:

$$ f^x = 4\pi^2 \left( \mu + \frac{1}{\lambda} + \frac{2\mu}{\lambda^2} \right) \cos\left(\frac{2\pi}{\lambda}x\right) \sin(2\pi y) + 4\pi^2 \left( 1 + \frac{\mu}{\lambda} \right) \cos(2\pi x) \sin\left(\frac{2\pi}{\lambda}y\right), $$

$$ f^y = 4\pi^2 \left( \mu + \frac{1}{\lambda} + \frac{2\mu}{\lambda^2} \right) \sin(2\pi x) \cos\left(\frac{2\pi}{\lambda}y\right) + 4\pi^2 \left( 1 + \frac{\mu}{\lambda} \right) \sin\left(\frac{2\pi}{\lambda}x\right) \cos(2\pi y). $$

Similarly to the observations of [36], from Figure 4 we recover the fact that the $\mathbb{P}^1$ finite element (i.e. the first order VEM scheme) suffers when $\lambda$ grows. For $\lambda \approx 10 \times 10^6$ ($v = 0.4999995$) numerical locking is observed, true convergence being only recovered for very small $h$. The NEM and CNEM outperforms the VEM scheme in this situation on the regular point cloud, however if convergence is maintained on the distorted point cloud the initial error on coarse networks is much larger than even for the locking VEM scheme. In the moderate regime $\lambda \approx 10 \times 10^3$ ($v = 0.4995$), the NEM and CNEM clearly outperforms the VEM scheme, even on the distorted point cloud, while in the neutral regime $\lambda = 1$ ($v = 0.25$) we recover results similar to the previous analytic test cases. Notice that the network generation algorithm from [31] produces well connected networks, with each cell node connected to all its immediate neighbors. Thus we have in general $\text{card}(\mathcal{F}_K) > d + 1$ which probably explains why the NEM avoids numerical locking: there is enough degrees of freedom per cell to represent accurately divergence free displacements.
<figure><img src="image_7.png" /><figcaption>**Figure 4.** Convergence curves for incompressible analytical test case.</figcaption></figure>

<table><caption>**Table 2.** Approximate orders of convergence on Delaunay meshes</caption><thead><tr><th></th><th>Sinusoidal</th><th>Exponential</th></tr></thead><tbody><tr><td>NEM</td><td>2.09</td><td>2.14</td></tr><tr><td>CNEM</td><td>2.09</td><td>2.09</td></tr><tr><td>VEM</td><td>2.00</td><td>2.01</td></tr><tr><td>CVEM</td><td>2.08</td><td>2.09</td></tr></tbody></table>

### 6.3. The NEM for linear elasticity as a subface-based scheme on meshes

To conclude this numerical section, we display on Figures 5 and 6 and on Table 2 a comparison of convergence curves between the VEM schemes and the small stencil, subface-based NEM scheme described in Section 4.2.1, for the sinusoidal and exponential analytic test cases.

We see that the proposed small-stencil mesh-based NEM schemes are perfectly stable and convergent, and thus are a competitive alternative to the hybrid scheme proposed in [35].
<figure><img src="image_7.png" /><figcaption>Figure 4. Convergence curves for incompressible analytical test case.</figcaption></figure>

<table><caption>Table 2. Approximate orders of convergence on Delaunay meshes</caption><thead><tr><th></th><th>Sinusoidal</th><th>Exponential</th></tr></thead><tbody><tr><th>NEM</th><td>2.09</td><td>2.14</td></tr><tr><th>CNEM</th><td>2.09</td><td>2.09</td></tr><tr><th>VEM</th><td>2.00</td><td>2.01</td></tr><tr><th>CVEM</th><td>2.08</td><td>2.09</td></tr></tbody></table>

### 6.3. The NEM for linear elasticity as a subface-based scheme on meshes

To conclude this numerical section, we display on Figures 5 and 6 and on Table 2 a comparison of convergence curves between the VEM schemes and the small stencil, subface-based NEM scheme described in Section 4.2.1, for the sinusoidal and exponential analytic test cases.

We see that the proposed small-stencil mesh-based NEM schemes are perfectly stable and convergent, and thus are a competitive alternative to the hybrid scheme proposed in [35].
<figure><img src="image_8.png" /><figcaption>**Figure 5.** Convergence curves for the sinusoidal test case.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 6.** Convergence curves for the exponential test case.</figcaption></figure>

## 7. Conclusion and perspectives

We have presented an extension to the system of linear elasticity of the meshless network element method. Its stability relies on a sufficient condition on the network's connectivity that is easily enforced using existing network generation algorithms. We have then presented a modified version that enjoys a finite volume like flux formulation. Used in a mesh-based context, the NEM schemes also offer small stencil alternatives to existing first order face based schemes, at the expense of using unknowns associated with at least *d* subfaces of each face instead of
<figure><img src="image_8.png" /><figcaption>**Figure 5.** Convergence curves for the sinusoidal test case.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 6.** Convergence curves for the exponential test case.</figcaption></figure>

## 7. Conclusion and perspectives

We have presented an extension to the system of linear elasticity of the meshless network element method. Its stability relies on a sufficient condition on the network's connectivity that is easily enforced using existing network generation algorithms. We have then presented a modified version that enjoys a finite volume like flux formulation. Used in a mesh-based context, the NEM schemes also offer small stencil alternatives to existing first order face based schemes, at the expense of using unknowns associated with at least *d* subfaces of each face instead of
one unknown per face. Numerical results illustrate the good behavior of the method, even on distorted networks or in the incompressible limit.

## Declaration of interests

The authors do not work for, advise, own shares in, or receive funds from any organization that could benefit from this article, and have declared no affiliations other than their research organizations.

## Dedication

The manuscript was written through contributions of all authors. All authors have given approval to the final version of the manuscript.

## References

[1] R. Glowinski, T.-W. Pan, J. Périaux, "A fictitious domain method for Dirichlet problem and applications", *Comput. Methods Appl. Mech. Eng.* **111** (1994), no. 3-4, p. 283-303.

[2] V. Girault, R. Glowinski, "Error analysis of a fictitious domain method applied to a Dirichlet problem", *Japan. J. Ind. Appl. Math.* **12** (1995), no. 3, p. 487-514.

[3] L. Lucy, "A numerical approach to the testing of the fission hypothesis", *Astron. J.* **82** (1977), p. 1013-1024.

[4] R. Gingold, J. Monaghan, "Smoothed particle hydrodynamics: theory and application to non-spherical stars", *Mon. Not. R. Astron. Soc.* **181** (1977), p. 375-389.

[5] W. Liu, S. Jun, Y. Zhang, "Reproducing kernel particle methods", *Int. J. Numer. Methods Fluids* **20** (1995), p. 1081-1106.

[6] P. Jensen, "Finite difference techniques for variable grids", *Comput. Struct.* **2** (1972), p. 17-29.

[7] J. Melenk, *On Approximation in Meshless Methods*, Springer, Berlin, Heidelberg, 2005, 65-141 pages.

[8] E. Oñate, S. Idelsohn, O. Zienkiewicz, R. Taylor, "A finite point method in computational mechanics, Applications to convective transport and fluid flow", *Int. J. Numer. Methods Eng.* **39** (1996), p. 3839-3866.

[9] R. Trobec, G. Kosec, *Parallel Scientific Computing: Theory, Algorithms, and Applications of Mesh Based and Meshless Methods*, Springer, Cham, 2015.

[10] R. Hardy, "Multiquadric equations of topography and other irregular surfaces", *J. Geophys. Res.* **76** (1971), p. 1905-1915.

[11] R. Hardy, "Theory and applications of the multiquadric-biharmonic method 20 years of discovery 1968-1988", *Comput. Math. Appl.* **19** (1990), p. 163-208.

[12] C. Franke, R. Schaback, "Solving partial differential equations by collocation using radial basis functions", *Appl. Math. Comput.* **93** (1998), p. 73-82.

[13] B. Nayroles, G. Touzot, P. Villon, "Generalizing the finite element method: Diffuse approximation and diffuse elements", *Comput. Mech.* **10** (1992), p. 307-318.

[14] T. Belytschko, Y. Lu, L. Gu, "Element-free Galerkin methods", *Int. J. Numer. Methods Eng.* **37** (1994), p. 229-256.

[15] Y. Lu, T. Belytschko, L. Gu, "A new implementation of the element free Galerkin method", *Comput. Methods Appl. Mech. Eng.* **113** (1994), p. 397-414.

[16] J. Melenk, I. Babuška, "The partition of unity finite element method: Basic theory and applications", *Comput. Methods Appl. Mech. Eng.* **139** (1996), no. 1-4, p. 289-314.

[17] J. Melenk, I. Babuška, "The partition of unity method", *Int. J. Numer. Methods Eng.* **40** (1997), p. 727-758.

[18] S. Atluri, T. Zhu, "A new meshless local Petrov-Galerkin (MLPG) approach in computational mechanics", *Comput. Mech.* **22** (1998), no. 2, p. 117-127.

[19] S. Atluri, S. Shen, "The meshless local Petrov-Galerkin (MLPG) method: a simple and less-costly alternative to the finite element and boundary element methods", *Comput. Model Eng. Sci.* **3** (2002), no. 1, p. 11-51.

[20] J.-S. Chen, M. Hillman, S.-W. Chi, "Meshfree methods: progress made after 20 years", *J. Eng. Mech.* **143** (2017), no. 4, article no. 04017001.

[21] J. Coatléven, "Principles of a network element method", *J. Comput. Phys.* **433** (2021), article no. 110197.

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

[23] O. Diyankov, "Uncertain grid method for numerical solution of PDEs", Tech. report, NeurOK Software, 2008.
one unknown per face. Numerical results illustrate the good behavior of the method, even on distorted networks or in the incompressible limit.

## Declaration of interests

The authors do not work for, advise, own shares in, or receive funds from any organization that could benefit from this article, and have declared no affiliations other than their research organizations.

## Dedication

The manuscript was written through contributions of all authors. All authors have given approval to the final version of the manuscript.

## References

[1] R. Glowinski, T.-W. Pan, J. Périaux, "A fictitious domain method for Dirichlet problem and applications", *Comput. Methods Appl. Mech. Eng.* **111** (1994), no. 3-4, p. 283-303.

[2] V. Girault, R. Glowinski, "Error analysis of a fictitious domain method applied to a Dirichlet problem", *Japan. J. Ind. Appl. Math.* **12** (1995), no. 3, p. 487-514.

[3] L. Lucy, "A numerical approach to the testing of the fission hypothesis", *Astron. J.* **82** (1977), p. 1013-1024.

[4] R. Gingold, J. Monaghan, "Smoothed particle hydrodynamics: theory and application to non-spherical stars", *Mon. Not. R. Astron. Soc.* **181** (1977), p. 375-389.

[5] W. Liu, S. Jun, Y. Zhang, "Reproducing kernel particle methods", *Int. J. Numer. Methods Fluids* **20** (1995), p. 1081-1106.

[6] P. Jensen, "Finite difference techniques for variable grids", *Comput. Struct.* **2** (1972), p. 17-29.

[7] J. Melenk, *On Approximation in Meshless Methods*, Springer, Berlin, Heidelberg, 2005, 65-141 pages.

[8] E. Oñate, S. Idelsohn, O. Zienkiewicz, R. Taylor, "A finite point method in computational mechanics, Applications to convective transport and fluid flow", *Int. J. Numer. Methods Eng.* **39** (1996), p. 3839-3866.

[9] R. Trobec, G. Kosec, *Parallel Scientific Computing: Theory, Algorithms, and Applications of Mesh Based and Meshless Methods*, Springer, Cham, 2015.

[10] R. Hardy, "Multiquadric equations of topography and other irregular surfaces", *J. Geophys. Res.* **76** (1971), p. 1905-1915.

[11] R. Hardy, "Theory and applications of the multiquadric-biharmonic method 20 years of discovery 1968-1988", *Comput. Math. Appl.* **19** (1990), p. 163-208.

[12] C. Franke, R. Schaback, "Solving partial differential equations by collocation using radial basis functions", *Appl. Math. Comput.* **93** (1998), p. 73-82.

[13] B. Nayroles, G. Touzot, P. Villon, "Generalizing the finite element method: Diffuse approximation and diffuse elements", *Comput. Mech.* **10** (1992), p. 307-318.

[14] T. Belytschko, Y. Lu, L. Gu, "Element-free Galerkin methods", *Int. J. Numer. Methods Eng.* **37** (1994), p. 229-256.

[15] Y. Lu, T. Belytschko, L. Gu, "A new implementation of the element free Galerkin method", *Comput. Methods Appl. Mech. Eng.* **113** (1994), p. 397-414.

[16] J. Melenk, I. Babuška, "The partition of unity finite element method: Basic theory and applications", *Comput. Methods Appl. Mech. Eng.* **139** (1996), no. 1-4, p. 289-314.

[17] J. Melenk, I. Babuška, "The partition of unity method", *Int. J. Numer. Methods Eng.* **40** (1997), p. 727-758.

[18] S. Atluri, T. Zhu, "A new meshless local Petrov-Galerkin (MLPG) approach in computational mechanics", *Comput. Mech.* **22** (1998), no. 2, p. 117-127.

[19] S. Atluri, S. Shen, "The meshless local Petrov-Galerkin (MLPG) method: a simple and less-costly alternative to the finite element and boundary element methods", *Comput. Model Eng. Sci.* **3** (2002), no. 1, p. 11-51.

[20] J.-S. Chen, M. Hillman, S.-W. Chi, "Meshfree methods: progress made after 20 years", *J. Eng. Mech.* **143** (2017), no. 4, article no. 04017001.

[21] J. Coatléven, "Principles of a network element method", *J. Comput. Phys.* **433** (2021), article no. 110197.

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

[23] O. Diyankov, "Uncertain grid method for numerical solution of PDEs", Tech. report, NeurOK Software, 2008.
[24] E. K. yu Chiu, Q. Wang, R. Hu, A. Jameson, "A conservative mesh-free scheme and generalized framework for conservation laws", *SIAM J. Sci. Comput.* **34** (2012), no. 6, p. 2896-2916.

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

[26] N. Trask, M. Perego, P. Bochev, "A high-order staggered meshless method for elliptic problems", *SIAM J. Sci. Comput.* **39** (2017), no. 2, p. 479-502.

[27] N. Trask, P. Bochev, M. Perego, "A conservative, consistent, and scalable mesh-free mimetic method", *J. Comput. Phys.* **409** (2020), article no. 109187.

[28] J. Coatléven, "A network element method for heterogeneous and anisotropic diffusion-reaction problems", *J. Comput. Phys.* **470** (2022), article no. 111597.

[29] J. Coatléven, "A conservative network element method for diffusion-reaction problems", *ESAIM Math. Model. Numer. Anal.* **57** (2023), no. 4, p. 2007-2040.

[30] J. Coatléven, "Basic convergence theory for the network element method", *ESAIM Math. Model. Numer. Anal.* **55** (2021), no. 5, p. 2503-2533.

[31] J. Coatléven, "On network and geometry generation for the network element method", 2022, preprint.

[32] L. B. da Veiga, F. Brezzi, A. Cangiani, G. Manzini, L. Marini, A. Russo, "Basic principles of virtual element methods", *Math. Models Methods Appl. Sci.* **23** (2013), no. 1, p. 199-214.

[33] L. B. da Veiga, F. Brezzi, L. D. Marini, "Virtual elements for linear elasticity problems", *SIAM J. Numer. Anal.* **51** (2013), no. 2, p. 794-812.

[34] D. A. Di Pietro, R. Eymard, S. Lemaire, R. Masson, "Hybrid finite volume discretization of linear elasticity models on general meshes", in *Finite Vol. for Complex App. VI Problems & Perspectives. Springer Proceedings in Mathematics*, vol. 4, Springer, Berlin, Heidelberg, 2011.

[35] M. Botti, D. Di Pietro, A. Guglielmana, "A low-order nonconforming method for linear elasticity on general meshes", *Comput. Methods Appl. Mech. Eng.* **354** (2019), p. 96-118.

[36] S. Lemaire, *Discrétisations non-conformes d'un modèle poromécanique sur maillages généraux*, Phd thesis, Univerisé Paris-Est, 2013.

[37] B. Mavrič, B. Sarler, "Local radial basis function collocation method for linear thermoelasticity in two dimensions", *Int J. Numer. Method Heat Fluid Flow* **25** (2015), p. 1488-1510.

[38] J. Slak, G. Kosec, "Refined meshless local strong form solution of Cauchy-Navier equation on an irregular domain", *Eng. Anal. Boundary Elem.* **100** (2019), p. 3-13.

[39] J. Slak, G. Kosec, "Adaptive radial basis function-generated finite differences method for contact problems", *Int. J. Numer. Method Eng.* **119** (2019), no. 7, p. 661-686.

[40] G. Kosec, J. Slak, M. Depolli, R. Trobec, K. Pereira, S. Tomar, T. Jacquemin, S. Bordas, M. A. Wahab, "Weak and strong form meshless methods for linear elastic problem under fretting contact conditions", *Tribol. Int.* **138** (2019), p. 392-402.

[41] T. Oliveira, W. Vélez, E. Santana, T. Araújo, F. Mendonça, A. Portela, "A local mesh free method for linear elasticity and fracture mechanics", *Eng. Anal. Bound. Elem.* **101** (2019), p. 221-242.

[42] L. Sang-Ho, Y. Young-Cheol, "Meshfree point collocation method for elasticity and crack problems", *Int. J. Numer. Methods Eng.* **61** (2004), p. 22-48.

[43] V. Nguyen, M. Duflot, "Meshless methods: a review and computer implementation aspects", *Math. Comput. Sim.* **79** (2008), no. 3, p. 763-813.

[44] W. McLean, *Strongly Elliptic Systems and Boundary Integral Equations*, Cambridge University Press, Cambridge, UK, 2000.

[45] S. C. Brenner, "Korn's inequalities for piecewise $H^1$ vector fields", *Math. Comput.* **73** (2004), no. 247, p. 1067-1087.

[46] R. Eymard, C. Guichard, R. Herbin, "Small-stencil 3D schemes for diffusive flows in porous media", *ESAIM Math. Model. Numer. Anal.* **46** (2011), no. 2, p. 265-290.

[47] J. Coatléven, "A virtual volume method for heterogeneous and anisotropic diffusion-reaction problems on general meshes", *ESAIM Math. Model. Numer. Anal.* **51** (2017), p. 797-824.

[48] R. Eymard, T. Gallouët, R. Herbin, "Finite volume methods", in *Techniques of Scientific Computing Part III* (P. G. Ciarlet, J.-L. Lions, eds.), Handbook of Numerical Analysis, North-Holland, Amsterdam, 2000, p. 713-1020.

[49] D. A. Di Pietro, A. Ern, *Mathematical Aspects of Discontinuous Galerkin Methods*, Springer, Berlin, Heidelberg, 2012.

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

[51] E. Stein, *Singular Integrals and Differentiability Properties of Functions*, Princeton University Press, Princeton, NJ, 1970.

[52] J. Slak, G. Kosec, "On generation of node distributions for meshless PDE discretizations", *SIAM J. Sci. Comput.* **41** (2019), no. 5, p. A3202-A3229.
[24] E. K. yu Chiu, Q. Wang, R. Hu, A. Jameson, "A conservative mesh-free scheme and generalized framework for conservation laws", *SIAM J. Sci. Comput.* **34** (2012), no. 6, p. 2896-2916.

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

[26] N. Trask, M. Perego, P. Bochev, "A high-order staggered meshless method for elliptic problems", *SIAM J. Sci. Comput.* **39** (2017), no. 2, p. 479-502.

[27] N. Trask, P. Bochev, M. Perego, "A conservative, consistent, and scalable mesh-free mimetic method", *J. Comput. Phys.* **409** (2020), article no. 109187.

[28] J. Coatléven, "A network element method for heterogeneous and anisotropic diffusion-reaction problems", *J. Comput. Phys.* **470** (2022), article no. 111597.

[29] J. Coatléven, "A conservative network element method for diffusion-reaction problems", *ESAIM Math. Model. Numer. Anal.* **57** (2023), no. 4, p. 2007-2040.

[30] J. Coatléven, "Basic convergence theory for the network element method", *ESAIM Math. Model. Numer. Anal.* **55** (2021), no. 5, p. 2503-2533.

[31] J. Coatléven, "On network and geometry generation for the network element method", 2022, preprint.

[32] L. B. da Veiga, F. Brezzi, A. Cangiani, G. Manzini, L. Marini, A. Russo, "Basic principles of virtual element methods", *Math. Models Methods Appl. Sci.* **23** (2013), no. 1, p. 199-214.

[33] L. B. da Veiga, F. Brezzi, L. D. Marini, "Virtual elements for linear elasticity problems", *SIAM J. Numer. Anal.* **51** (2013), no. 2, p. 794-812.

[34] D. A. Di Pietro, R. Eymard, S. Lemaire, R. Masson, "Hybrid finite volume discretization of linear elasticity models on general meshes", in *Finite Vol. for Complex App. VI Problems & Perspectives. Springer Proceedings in Mathematics*, vol. 4, Springer, Berlin, Heidelberg, 2011.

[35] M. Botti, D. Di Pietro, A. Guglielmana, "A low-order nonconforming method for linear elasticity on general meshes", *Comput. Methods Appl. Mech. Eng.* **354** (2019), p. 96-118.

[36] S. Lemaire, *Discrétisations non-conformes d'un modèle poromécanique sur maillages généraux*, Phd thesis, Univerisé Paris-Est, 2013.

[37] B. Mavrič, B. Sarler, "Local radial basis function collocation method for linear thermoelasticity in two dimensions", *Int J. Numer. Method Heat Fluid Flow* **25** (2015), p. 1488-1510.

[38] J. Slak, G. Kosec, "Refined meshless local strong form solution of Cauchy-Navier equation on an irregular domain", *Eng. Anal. Boundary Elem.* **100** (2019), p. 3-13.

[39] J. Slak, G. Kosec, "Adaptive radial basis function-generated finite differences method for contact problems", *Int. J. Numer. Method Eng.* **119** (2019), no. 7, p. 661-686.

[40] G. Kosec, J. Slak, M. Depolli, R. Trobec, K. Pereira, S. Tomar, T. Jacquemin, S. Bordas, M. A. Wahab, "Weak and strong form meshless methods for linear elastic problem under fretting contact conditions", *Tribol. Int.* **138** (2019), p. 392-402.

[41] T. Oliveira, W. Vélez, E. Santana, T. Araújo, F. Mendonça, A. Portela, "A local mesh free method for linear elasticity and fracture mechanics", *Eng. Anal. Bound. Elem.* **101** (2019), p. 221-242.

[42] L. Sang-Ho, Y. Young-Cheol, "Meshfree point collocation method for elasticity and crack problems", *Int. J. Numer. Methods Eng.* **61** (2004), p. 22-48.

[43] V. Nguyen, M. Duflot, "Meshless methods: a review and computer implementation aspects", *Math. Comput. Sim.* **79** (2008), no. 3, p. 763-813.

[44] W. McLean, *Strongly Elliptic Systems and Boundary Integral Equations*, Cambridge University Press, Cambridge, UK, 2000.

[45] S. C. Brenner, "Korn's inequalities for piecewise $H^1$ vector fields", *Math. Comput.* **73** (2004), no. 247, p. 1067-1087.

[46] R. Eymard, C. Guichard, R. Herbin, "Small-stencil 3D schemes for diffusive flows in porous media", *ESAIM Math. Model. Numer. Anal.* **46** (2011), no. 2, p. 265-290.

[47] J. Coatléven, "A virtual volume method for heterogeneous and anisotropic diffusion-reaction problems on general meshes", *ESAIM Math. Model. Numer. Anal.* **51** (2017), p. 797-824.

[48] R. Eymard, T. Gallouët, R. Herbin, "Finite volume methods", in *Techniques of Scientific Computing Part III* (P. G. Ciarlet, J.-L. Lions, eds.), Handbook of Numerical Analysis, North-Holland, Amsterdam, 2000, p. 713-1020.

[49] D. A. Di Pietro, A. Ern, *Mathematical Aspects of Discontinuous Galerkin Methods*, Springer, Berlin, Heidelberg, 2012.

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

[51] E. Stein, *Singular Integrals and Differentiability Properties of Functions*, Princeton University Press, Princeton, NJ, 1970.

[52] J. Slak, G. Kosec, "On generation of node distributions for meshless PDE discretizations", *SIAM J. Sci. Comput.* **41** (2019), no. 5, p. A3202-A3229.