# A unified enrichment approach of the standard three-node triangular element

Francesco Dell'Accio, Filomena Di Tommaso, Allal Guessab, Federico Nudo

► To cite this version:

Francesco Dell'Accio, Filomena Di Tommaso, Allal Guessab, Federico Nudo. A unified enrichment approach of the standard three-node triangular element. Applied Numerical Mathematics: an IMACS journal, 2023, 187, pp.1-23. 10.1016/j.apnum.2023.02.001. hal-03991658

HAL Id: hal-03991658

https://hal.science/hal-03991658v1

Submitted on 6 Sep 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
# A unified enrichment approach of the standard three-node triangular element

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

Check for
updates

Francesco Dell'Accio<sup>a,b,*</sup>, Filomena Di Tommaso<sup>a,b</sup>, Allal Guessab<sup>c</sup>, Federico Nudo<sup>a,c</sup>

## ARTICLE INFO

**Article history:**

Received 15 November 2022

Received in revised form 26 January 2023

Accepted 3 February 2023

Available online 9 February 2023

**Keywords:**

Enriched finite element method

Error bound

Polynomial enrichment

Non-polynomial enrichment

Triangular linear element

## ABSTRACT

The aim of this paper is to unify the ideas and to extend to a more general setting the work done in (Dell'Accio et al., 2023 [14]) for a polynomial enrichment of the standard three-node triangular element (triangular linear element) using line integrals and quadratic polynomials. More precisely, we introduce a new class of nonconforming finite elements by enriching the class of linear polynomial functions with additional functions which are not necessarily polynomials. We provide a simple condition on the enrichment functions, which is both necessary and sufficient, that guarantees the existence of a family of such enriched elements. Several sets of admissible enriched functions that satisfy the admissibility condition are also provided, together with the explicit expression of the related approximation error. Our main result shows that the approximation error can be decomposed into two parts: the first one is related to the linear triangular element while the second one depends on the enrichment functions. This representation of the approximation error allows us to derive sharp error bounds in both $L^\infty$ and $L^1$ norms, with explicit constants, for continuously differentiable functions with Lipschitz continuous gradients. These bounds are proportional to the second and the fourth power of the circumcircle radius of the triangle, respectively. We also provide explicit expressions of these bounds in terms of the circumcircle diameter and the sum of squares of the triangle edge lengths.

© 2023 The Authors. Published by Elsevier B.V. on behalf of IMACS. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/).

## 1. Introduction

The finite element method (FEM) has achieved great success in many fields of science and technology. It is a very useful and powerful tool for calculating accurate approximations of solutions of partial differential equations [8]. One of the simplest finite elements is the traditional low-order finite element using three-node triangular element in two-dimensional Euclidean space [8]. It is defined locally as the triplet

$$ (T, P_1(T), \Sigma_T), $$

<footnote><sup>a</sup> Department of Mathematics and Computer Science, University of Calabria, Rende (CS), Italy</footnote>

<footnote><sup>b</sup> Istituto per le Applicazioni del Calcolo "Mauro Picone", Naples Branch, C.N.R. National Research Council of Italy, Napoli, Italy</footnote>

<footnote><sup>c</sup> Laboratoire de Mathématiques et de leurs Applications, UMR CNRS 4152, Université de Pau et des Pays de l'Adour (UPPA), 64000 Pau, France</footnote>

<footnote>* Corresponding author at: Department of Mathematics and Computer Science, University of Calabria, Rende (CS), Italy.</footnote>

<footnote>E-mail addresses: francesco.dellaccio@unical.it (F. Dell'Accio), filomena.ditommaso@unical.it (F. Di Tommaso), allal.guessab@univ-pau.fr (A. Guessab), federico.nudo@unical.it (F. Nudo).</footnote>
where $T$ is any non degenerate triangle of a triangulation, $P_1(T)$ is the space of bivariate linear polynomials and $\Sigma_T$ is the set of functional evaluations at the vertices of $T$, also called degrees of freedom of the finite element. The three-node triangular element is widely used in the applications, however it does not yield satisfactory results due to the low order of approximation of its trial functions. Basing on this simple element, alternative approaches have been proposed in an effort to obtain more effective finite element methods. One approach is to enrich the traditional triangular element $(T, P_1(T), \Sigma_T)$ with special enrichment functions and additional appropriate degrees of freedom, which are added to $P_1(T)$ and $\Sigma_T$, respectively, in order to obtain a more accurate element (for an overview of the relevant literature see, for instance, [2–4,21]). This is a natural way to improve the performance of finite element analysis and it is very important in the applications, see [10,18]. Motivated by the polynomial enrichment proposed in [14], which is applied to improve the triangular Shepard method (see [6,12,13,16,27] for other approaches) and by recent works on finite element approximations via the use of enrichment functions [2,4,21], in this paper we propose a general polynomial and non-polynomial approach for the enrichment of the standard triangular linear element. While in theory this approach has considerable generality, we focus on the two-dimensional case, from which both one-dimensional and three-dimensional cases can similarly be developed. More precisely, the goal of the current study is to unify the ideas and to extend to a more general setting the polynomial enrichment presented in [14], by using more general polynomial and non-polynomial enrichment functions.

The paper is organized as follows. Section 2 begins by presenting the traditional three-node triangular element and by recalling some classical results and notations of the finite element method, used throughout the paper. Then we propose a generalization, using both polynomial and non-polynomial functions, of the polynomial enrichment for the standard triangular linear element presented in [14]. Furthermore we provide explicit expressions of the basis functions in terms of barycentric coordinates and an explicit error representation. Thereafter, in Section 3 we present the general framework and extend the results of Section 2 to non-polynomial enrichments of the standard triangular linear element. Here an essential problem is to establish necessary and sufficient conditions, on the enrichment functions, to guarantee the existence of well-defined enriched finite elements. This will be done by deriving, in Theorem 3.1, a general characterization result, which significantly extends and improves the corresponding characterization Theorem of the polynomial enrichment case introduced in [14]. In Section 3.1 we derive an explicit error representation formula for the non-polynomial enriched approximant, in Proposition 3.4, and provide explicit bounds for the approximation error, both in $L^1$-norm and $L^\infty$-norm, in Corollary 3.12 and Corollary 3.19, respectively. We also list several relevant a priori error estimates. We discuss the quality of the presented bounds and obtain, in Corollary 3.19, sharp error bounds in $L^1$-norm, while, in Corollary 3.21, we extend this result for some special choices of admissible enrichment functions. In fact, we show that both in the $L^\infty$-norm and $L^1$-norm cases these bounds are respectively proportional to the square and the fourth power of the circumcircle of the triangle, respectively, independently of admissible enrichment functions. To our knowledge, the form of the bounds and the computational approach based on Leibniz and Weitzenböck inequalities, are new in the context of the error estimate.

## 2. Polynomial enrichment of the standard triangular linear element

### 2.1. The standard triangular linear element

Given a non-degenerate triangle $T \subset \mathbb{R}^2$ with vertices $v_i = (x_i, y_i)$, we denote by $\Gamma_i$ the edge of $T$ opposite to the vertex $v_i$, $i = 1, 2, 3$ and by $\{\lambda_i(x)\}_i$ the barycentric coordinates of $T$, defined as follows

$$ \lambda_1(\mathbf{x}) := \frac{\begin{vmatrix} 1 & 1 & 1 \\ x & x_2 & x_3 \\ y & y_2 & y_3 \end{vmatrix}}{\begin{vmatrix} 1 & 1 & 1 \\ x_1 & x_2 & x_3 \\ y_1 & y_2 & y_3 \end{vmatrix}}, \quad \lambda_2(\mathbf{x}) := \frac{\begin{vmatrix} 1 & 1 & 1 \\ x_1 & x & x_3 \\ y_1 & y & y_3 \end{vmatrix}}{\begin{vmatrix} 1 & 1 & 1 \\ x_1 & x_2 & x_3 \\ y_1 & y_2 & y_3 \end{vmatrix}}, \quad \lambda_3(\mathbf{x}) := \frac{\begin{vmatrix} 1 & 1 & 1 \\ x_1 & x_2 & x \\ y_1 & y_2 & y \end{vmatrix}}{\begin{vmatrix} 1 & 1 & 1 \\ x_1 & x_2 & x_3 \\ y_1 & y_2 & y_3 \end{vmatrix}} $$

where $\mathbf{x} = (x, y) \in \mathbb{R}^2$. We denote by $\langle \cdot, \cdot \rangle$ the scalar product of the Euclidean space $\mathbb{R}^2$, and by $\|\cdot\|$ the relative norm. Then

$$ |\Gamma_1| = \|\mathbf{v}_2 - \mathbf{v}_3\|, \quad |\Gamma_2| = \|\mathbf{v}_1 - \mathbf{v}_3\|, \quad |\Gamma_3| = \|\mathbf{v}_1 - \mathbf{v}_2\|, $$

are the triangle $T$ edge lengths. We also denote by $|T|$ the area of the triangle $T$. Let

$$ P_1(T) := \operatorname{span}\{\lambda_1, \lambda_2, \lambda_3\} $$

be the space of linear polynomials in $\mathbb{R}^2$. The standard triangular linear element $P_1$ is the triple

$$ P_1 := (T, P_1(T), \Sigma_T), $$

where $\Sigma_T$ is the set of point evaluation functionals at the vertices of $T$, that is

$$ \Sigma_T := \{L_1(f) := f(\mathbf{v}_1), L_2(f) := f(\mathbf{v}_2), L_3(f) := f(\mathbf{v}_3)\}. $$
## 2.2. The general problem for polynomial enrichment of the standard triangular linear element

In this Section we generalize the quadratic polynomial enrichment of the standard triangular linear element proposed in [14] and based on the set of linear functionals, called degrees of freedom

$$ \tilde{\Sigma}_T := \{L_i, I_i, i = 1, 2, 3\}, $$

with

$$ L_i(f) = f(v_i), \quad f \in \tilde{P}_2(T), \quad i = 1, 2, 3, \qquad (1) $$

$$ I_i(f) := \frac{1}{|\Gamma_i|} \int_{\Gamma_i} f d\sigma, \quad f \in \tilde{P}_2(T), \quad i = 1, 2, 3, \qquad (2) $$

to the case of polynomials of degree greater than 2. We start with the case of cubic polynomials by fixing three linear polynomials $l_1, l_2, l_3$, which satisfy

$$ l_1(v_{12}) = l_2(v_{23}) = l_3(v_{13}) = 1, \qquad (3) $$

where $v_{12}, v_{23}$ and $v_{13}$ are the midpoints of the sides of $T$, and precisely

$$ v_{12} := \frac{v_1 + v_2}{2} \in \Gamma_3, \quad v_{23} := \frac{v_2 + v_3}{2} \in \Gamma_1, \quad v_{13} := \frac{v_1 + v_3}{2} \in \Gamma_2. $$

We introduce the enriched space $\tilde{P}_2(T)$ as follows:

$$ \tilde{P}_2(T) := P_1(T) \oplus \operatorname{span}\{l_1\lambda_2\lambda_3, l_2\lambda_1\lambda_3, l_3\lambda_1\lambda_2\} \qquad (4) $$

and we denote by AF3 the enriched polynomial element

$$ \mathrm{AF3} := (T, \tilde{P}_2(T), \tilde{\Sigma}_T). \qquad (5) $$

The following technical lemma will be useful in establishing the proof of Theorem 2.3 about the unisolvence of the element AF3.

**Lemma 2.1.** Let $p \in \tilde{P}_2(T)$, then we have

$$ I_1(p) = \frac{1}{6}(p(v_2) + 4p(v_{23}) + p(v_3)), \qquad (6) $$

$$ I_2(p) = \frac{1}{6}(p(v_1) + 4p(v_{13}) + p(v_3)), \qquad (7) $$

$$ I_3(p) = \frac{1}{6}(p(v_1) + 4p(v_{12}) + p(v_2)). \qquad (8) $$

**Proof.** Let us prove the identity (6). We parametrize the edge $\Gamma_1$ by

$$ t \to tv_3 + (1-t)v_2, \quad t \in [0, 1]. $$

Since the restriction of $p$ to the edge $\Gamma_1$ is a cubic polynomial in $t$, the Simpson's rule provides exact results for the integral, and then by equation (2) we get

$$ I_1(p) = \frac{1}{|\Gamma_1|} \int_{\Gamma_1} p d\sigma = \int_{0}^{1} p(tv_3 + (1-t)v_2)dt = \frac{1}{6}(p(v_2) + 4p(v_{23}) + p(v_3)). $$

The same argument can be used to prove identities (7) and (8). $\square$

**Remark 2.2.** Note that, if $p \in \tilde{P}_2(T)$ and $L_j(p) = I_j(p) = 0, j = 1, 2, 3$, then, by Lemma 2.1, we get

$$ p(v_{23}) = p(v_{13}) = p(v_{12}) = 0, \qquad (9) $$

that is, $p$ vanishes at the midpoints of the sides of $T$.

**Theorem 2.3.** The triple AF3 = $(T, \tilde{P}_2(T), \tilde{\Sigma}_T)$ is a unisolvent element.
**Proof.** We have to show that $\tilde{\Sigma}_T$ is $\tilde{P}_2(T)$-unisolvent, i.e., if $p \in \tilde{P}_2(T)$ and

$$L_j(p) = 0, \quad j = 1, 2, 3, \tag{10}$$

$$I_j(p) = 0, \quad j = 1, 2, 3, \tag{11}$$

then $p \equiv 0$ [11, Ch. 2]. Let $p \in \tilde{P}_2(T)$ and assume that (10) and (11) hold. By definition, $p$ can be represented as

$$p = \alpha_1 \lambda_1 + \alpha_2 \lambda_2 + \alpha_3 \lambda_3 + \beta_1 l_1 \lambda_2 \lambda_3 + \beta_2 l_2 \lambda_1 \lambda_3 + \beta_3 l_3 \lambda_1 \lambda_2 \tag{12}$$

for some constants $\alpha_i, \beta_i \in \mathbb{R}$, $i = 1, 2, 3$. Since the barycentric coordinates satisfy Lagrange property, that is $L_j(\lambda_i) = \lambda_i(v_j) = \delta_{ij}$, where $\delta_{ij}$ is the Kronecker delta symbol, then (10) implies $\alpha_1 = \alpha_2 = \alpha_3 = 0$ and $p$ reduces to

$$p = \beta_1 l_1 \lambda_2 \lambda_3 + \beta_2 l_2 \lambda_1 \lambda_3 + \beta_3 l_3 \lambda_1 \lambda_2.$$

On the other hand, $I_1(p) = 0$ and $L_2(p) = L_3(p) = 0$ imply, by (6), $p(v_{23}) = 0$. Now, since $\lambda_2(v_{23}) = \lambda_3(v_{23}) = \frac{1}{2}$ and $\lambda_1(v_{23}) = 0$, we get

$$0 = \frac{4}{6} p(v_{23}) = \frac{\beta_1}{6} l_1(v_{23}) = \frac{\beta_1}{6},$$

where in the last equality, we used the conditions (3). Thus $\beta_1 = 0$. A similar argument can be used to show that $\beta_2 = \beta_3 = 0$, therefore $p \equiv 0$ and the proof of the theorem is completed. □

**Remark 2.4.** We notice that if we take $l_1 = l_2 = l_3 = 1$, then the enriched space $\tilde{P}_2(T)$ becomes the standard space of quadratic polynomials $P_2(T)$ and, in this sense, the AF3 element (5) generalizes the element $(T, P_2(T), \tilde{\Sigma}_T)$ introduced in [14].

**Remark 2.5.** The enriched function space $\tilde{P}_2(T)$ satisfies the following properties:

i) it is a subspace of $P_3(T)$, the space of cubic polynomials, and hence the restriction of any $p \in \tilde{P}_2(T)$ to each side of $T$ is a cubic polynomial in one variable;

ii) it contains the set of linear polynomials;

iii) the nonlinear terms in the expression (12) of $p \in \tilde{P}_2(T)$ vanish at the vertices of $T$.

Theorem 2.3 can be stated equivalently by saying that the functionals of $\tilde{\Sigma}_T$ are linearly independent in the dual space $\tilde{P}_2(T)^*$. Then, there exists a basis $\{\varphi_i, \phi_i, i = 1, 2, 3\}$ of $\tilde{P}_2(T)$ which satisfies

$$L_i(\varphi_j) = \delta_{ij}, \quad I_i(\varphi_j) = 0, \quad i, j = 1, 2, 3, \tag{13}$$

$$L_i(\phi_j) = 0, \quad I_i(\phi_j) = \delta_{ij}, \quad i, j = 1, 2, 3. \tag{14}$$

**Theorem 2.6.** The basis functions $\varphi_i, \phi_i, i = 1, 2, 3$, of $(T, \tilde{P}_2(T), \tilde{\Sigma}_T)$ satisfying (13) and (14) have the following expressions

$$\varphi_1 = \lambda_1(1 - 3l_3\lambda_2 - 3l_2\lambda_3), \quad \varphi_2 = \lambda_2(1 - 3l_3\lambda_1 - 3l_1\lambda_3), \quad \varphi_3 = \lambda_3(1 - 3l_2\lambda_1 - 3l_1\lambda_2), \tag{15}$$

$$\varphi_1 = 6l_1\lambda_2\lambda_3, \quad \varphi_2 = 6l_2\lambda_1\lambda_3, \quad \varphi_3 = 6l_3\lambda_1\lambda_2. \tag{16}$$

**Proof.** Let us prove the first of identities (15). The element $\varphi_1 \in \tilde{P}_2(T)$ can be represented as

$$\varphi_1 = \alpha_1 \lambda_1 + \alpha_2 \lambda_2 + \alpha_3 \lambda_3 + \beta_1 l_1 \lambda_2 \lambda_3 + \beta_2 l_2 \lambda_1 \lambda_3 + \beta_3 l_3 \lambda_1 \lambda_2, \tag{17}$$

for some constants $\alpha_i, \beta_i \in \mathbb{R}$, $i = 1, 2, 3$. Since $L_1(\varphi_1) = 1$, $L_2(\varphi_1) = L_3(\varphi_1) = 0$, the Lagrange property of the barycentric coordinates implies $\alpha_1 = 1$, $\alpha_2 = \alpha_3 = 0$ and therefore

$$\varphi_1 = \lambda_1 + \beta_1 l_1 \lambda_2 \lambda_3 + \beta_2 l_2 \lambda_1 \lambda_3 + \beta_3 l_3 \lambda_1 \lambda_2.$$

Moreover, since $\lambda_1 = 0$ on $\Gamma_1$, from $I_1(\varphi_1) = 0$, by applying the Simpson's rule, we get

$$0 = \frac{1}{|\Gamma_1|} \int_{\Gamma_1} \beta_1 l_1 \lambda_2 \lambda_3 d\sigma = \frac{4}{6} \beta_1 l_1 (v_{23}) \lambda_2 (v_{23}) \lambda_3 (v_{23}) = \frac{\beta_1}{6},$$

which readily gives $\beta_1 = 0$. Similarly, since $\lambda_2 = 0$ on $\Gamma_2$, from $I_2(\varphi_1) = 0$ we get
$$ 
\begin{aligned}
0 &= \frac{1}{|\Gamma_2|} \int_{\Gamma_2} (\lambda_1 + \beta_2 l_2 \lambda_1 \lambda_3) d\sigma \\
&= \frac{1}{6} \left( 1 + 4 \left( \frac{1}{2} + \beta_2 l_2 (\mathbf{v}_{13}) \lambda_1 (\mathbf{v}_{13}) \lambda_3 (\mathbf{v}_{13}) \right) \right) \\
&= \frac{1}{6} (3 + \beta_2),
\end{aligned}
$$

that implies $\beta_2 = -3$. Finally, since $\lambda_3 = 0$ on $\Gamma_3$, from $l_3(\phi_1) = 0$ we get $\beta_3 = -3$, and then the first of (15) is proved. We now show that the first of identities (16) holds. Since $l_i(\phi_1) = 0$, $i = 1, 2, 3$, the Lagrange property of the barycentric coordinates implies that $\phi_1 \in \text{span}\{l_1\lambda_2\lambda_3, l_2\lambda_1\lambda_3, l_3\lambda_1\lambda_2\}$ and then it can be written as

$$ \phi_1 = \beta_1 l_1 \lambda_2 \lambda_3 + \beta_2 l_2 \lambda_1 \lambda_3 + \beta_3 l_3 \lambda_1 \lambda_2, $$

for some constants $\beta_i$, $i = 1, 2, 3$. Moreover, from $l_1(\phi_1) = 1$ we get

$$ 1 = \frac{1}{|\Gamma_1|} \int_{\Gamma_1} \beta_1 l_1 \lambda_2 \lambda_3 d\sigma = \frac{4}{6} \beta_1 l_1 (\mathbf{v}_{23}) \lambda_2 (\mathbf{v}_{23}) \lambda_3 (\mathbf{v}_{23}) = \frac{\beta_1}{6}, $$

hence $\beta_1 = 6$. From $l_2(\phi_1) = l_3(\phi_1) = 0$, we also get $\beta_2 = \beta_3 = 0$ and then the first of (16) is proved. We can obtain the expressions for the other functions by using symmetry arguments. $\square$

**Remark 2.7.** It is easily seen that conditions (3), which are sufficient for the existence of the enrichment (5) of the element $(T, P_1(T), \Sigma_T)$, can be replaced by the more general ones

$$ l_1(\mathbf{v}_{12}) \neq 0, \quad l_2(\mathbf{v}_{23}) \neq 0, \quad l_3(\mathbf{v}_{13}) \neq 0. \qquad (18) $$

As it has become clear during the discussion, general conditions (18) are also necessary.

More generally, let us consider three interior points $\mathbf{x}_1, \mathbf{x}_2, \mathbf{x}_3$, of the sides $\Gamma_1, \Gamma_2, \Gamma_3$ of the triangle $T$, respectively, and three linear polynomials $l_1, l_2, l_3$. We call these points (and polynomials) admissible if they generate a unisolvent enriched element, possibly with different basis functions. In analogy with (18), we introduce the nonvanishing conditions

$$ l_1(\mathbf{x}_1) \neq 0, \quad l_2(\mathbf{x}_2) \neq 0, \quad l_3(\mathbf{x}_3) \neq 0. \qquad (19) $$

The question arises whether conditions (19) are sufficient (and necessary) to generate a unisolvent enriched element. As we will see below, the answer to previous question is positive, but the differentiability class of the enrichment functions will depend on the position of the point $\mathbf{x}_1, \mathbf{x}_2, \mathbf{x}_3$. To prove this result, we use appropriate 1-point Gauss quadrature rules, instead of Simpson's rule, used in the case of the midpoints of the sides of $T$. The idea is rather simple and dates back to paper [19]. By assumption, there exist real numbers $\alpha_i, \beta_i > 1$, $i = 1, 2, 3$, such that

$$ 
\begin{aligned}
\mathbf{x}_1 &= \frac{\alpha_1}{\alpha_1 + \beta_1} \mathbf{v}_2 + \frac{\beta_1}{\alpha_1 + \beta_1} \mathbf{v}_3, \\
\mathbf{x}_2 &= \frac{\alpha_2}{\alpha_2 + \beta_2} \mathbf{v}_1 + \frac{\beta_2}{\alpha_2 + \beta_2} \mathbf{v}_3, \\
\mathbf{x}_3 &= \frac{\alpha_3}{\alpha_3 + \beta_3} \mathbf{v}_1 + \frac{\beta_3}{\alpha_3 + \beta_3} \mathbf{v}_2.
\end{aligned}
\qquad (20) $$

As in Section 2, for the sake of simplicity, we can assume that

$$ l_1(\mathbf{x}_1) = l_2(\mathbf{x}_2) = l_3(\mathbf{x}_3) = 1. \qquad (21) $$

We introduce the following more general enriched space

$$ \tilde{S}_2(T) = P_1(T) \oplus \text{span} \left\{ l_1 \lambda_2^{\alpha_1-1} \lambda_3^{\beta_1-1}, l_2 \lambda_1^{\alpha_2-1} \lambda_3^{\beta_2-1}, l_3 \lambda_1^{\alpha_3-1} \lambda_2^{\beta_3-1} \right\}, \qquad (22) $$

which includes $\tilde{P}_2(T)$ as particular case. In order to test whether or not $\tilde{\Sigma}_T$ is $\tilde{S}_2(T)$-unisolvent and to compute the basis of the enriched space, we recall the classical Euler beta function

$$ B(\alpha, \beta) = \int_{0}^{1} t^{\alpha-1} (1-t)^{\beta-1} dt, \quad \alpha, \beta > 0, \qquad (23) $$
which satisfies the key property

$$B(\alpha + 1, \beta) = \frac{\alpha}{\alpha + \beta} B(\alpha, \beta), \quad \alpha, \beta > 0. \tag{24}$$

The beta function is connected to the gamma function through the equation [1]

$$B(\alpha, \beta) = \frac{\Gamma(\alpha)\Gamma(\beta)}{\Gamma(\alpha + \beta)}, \tag{25}$$

and hence for positive integers $\alpha$ and $\beta$, we have

$$B(\alpha, \beta) = \frac{(\alpha - 1)!(\beta - 1)!}{(\alpha + \beta - 1)!}. \tag{26}$$

The following Proposition is crucial to prove that $\tilde{\Sigma}_T$ is $\tilde{S}_2(T)$-unisolvent.

**Proposition 2.8.** Under nonvanishing conditions (21), the enrichment functions of $\tilde{S}_2(T)$ satisfy the following delta properties

$$\begin{align}
\frac{1}{|\Gamma_i|} \int_{\Gamma_i} l_1 \lambda_2^{\alpha_1-1} \lambda_3^{\beta_1-1} d\sigma &= \delta_{1i} B(\alpha_1, \beta_1), \nonumber \\
\frac{1}{|\Gamma_i|} \int_{\Gamma_i} l_2 \lambda_1^{\alpha_2-1} \lambda_3^{\beta_2-1} d\sigma &= \delta_{2i} B(\alpha_2, \beta_2), \tag{27} \\
\frac{1}{|\Gamma_i|} \int_{\Gamma_i} l_3 \lambda_1^{\alpha_3-1} \lambda_3^{\beta_3-1} d\sigma &= \delta_{3i} B(\alpha_3, \beta_3), \nonumber
\end{align}$$

for each $i = 1, 2, 3$.

**Proof.** Let us prove the first of identities (27). Since $\alpha_1 > 1$ and $\beta_1 > 1$ then $\lambda_2^{\alpha_1-1} \lambda_3^{\beta_1-1}$ vanishes on $\Gamma_2$ and $\Gamma_3$, therefore the first of identities (27) holds for $i = 2, 3$. For $i = 1$, using the fact that $\lambda_2$ and $\lambda_3$ are affine functions, we get

$$\begin{align*}
\frac{1}{|\Gamma_1|} \int_{\Gamma_1} l_1 \lambda_2^{\alpha_1-1} \lambda_3^{\beta_1-1} d\sigma &= \int_0^1 (l_1 \lambda_2^{\alpha_1-1} \lambda_3^{\beta_1-1}) (t \mathbf{v}_2 + (1-t) \mathbf{v}_3) dt \\
&= \int_0^1 l_1(t \mathbf{v}_2 + (1-t) \mathbf{v}_3) t^{\alpha_1-1} (1-t)^{\beta_1-1} dt.
\end{align*}$$

For the 1-point Gauss quadrature in the interval $[0, 1]$ associated with the weight function $w_{\alpha_1, \beta_1}(t) = t^{\alpha_1-1}(1-t)^{\beta_1-1}$, the node is located at the point $\frac{\alpha_1}{\alpha_1+\beta_1}$, while the corresponding weight is equal to $B(\alpha_1, \beta_1)$ [25, Sect. 3.1]. Indeed, it suffices to determine the orthogonal polynomial $q(t)$ of degree 1 relative to the weight function $w_{\alpha_1, \beta_1}$ on the interval $[0, 1]$ and using (23) and (24), we get

$$q(t) = t - \frac{\alpha_1}{\alpha_1 + \beta_1}. \tag{28}$$

Since we assumed by (21) that $l_1(\mathbf{x}_1) = 1$, the result then follows from the exactness of the 1-point Gauss-Jacobi quadrature for polynomials of degree 1. $\square$

The following theorem extends the Theorem 2.3 to the case of the general configuration of three points (20).

**Theorem 2.9.** Let the linear polynomials $l_1, l_2, l_3$ satisfy the nonvanishing conditions (21) in the points $\mathbf{x}_1, \mathbf{x}_2, \mathbf{x}_3$. Then the element $(T, \tilde{S}_2(T), \tilde{\Sigma}_T)$ is unisolvent.

**Proof.** As the dimension of $\tilde{S}_2(T)$ is equal to the cardinality of $\tilde{\Sigma}_T$, it suffices to show that $\psi \in \tilde{S}_2(T)$ is identically zero if all the degrees of freedom (1) and (2) vanish when applied to $\psi$. The proof follows the same argument of Theorem 2.3 by using the general identities given in Proposition 2.8. It is therefore omitted here. $\square$
**Remark 2.10.** The nonvanishing conditions (19) are also necessary. Indeed, let us assume that conditions (19) do not hold, and without loss of generality, we can assume that $l_1(\mathbf{x}_1) = 0$. The function $e_1 = l_1\lambda_2^{\alpha_1-1}\lambda_3^{\beta_1-1} \in \tilde{S}_2(T)$, satisfies $e_1(\mathbf{v}_i) = 0, i = 1, 2, 3$ and

$$ \frac{1}{|\Gamma_2|} \int_{\Gamma_2} e_1 d\sigma = \frac{1}{|\Gamma_3|} \int_{\Gamma_3} e_1 d\sigma = 0, $$

since $e_1$ vanishes both on $\Gamma_2$ and $\Gamma_3$. We also have

$$ \frac{1}{|\Gamma_1|} \int_{\Gamma_1} e_1 d\sigma = \int_{0}^{1} t^{\alpha_1-1}(1-t)^{\beta_1-1} l_1(t\mathbf{v}_2 + (1-t)\mathbf{v}_3) dt = B(\alpha_1, \beta_1)l_1(\mathbf{x}_1) = 0. $$

Then $L_i(e_1) = I_j(e_1) = 0, i, j = 1, 2, 3$, and therefore $\tilde{\Sigma}_T$ is not $\tilde{S}_2(T)$-unisolvent since $e_1 \neq 0$.

By assuming that the nonvanishing conditions (19) hold, we introduce the following notations

$$ e_1 = l_1 \lambda_2^{\alpha_1-1} \lambda_3^{\beta_1-1}, \quad e_2 = l_2 \lambda_1^{\alpha_2-1} \lambda_3^{\beta_2-1}, \quad e_3 = l_3 \lambda_1^{\alpha_3-1} \lambda_2^{\beta_3-1}, \quad \gamma_i = B(\alpha_i, \beta_i), \quad i = 1, 2, 3. \quad (29) $$

Using the delta properties of the enriched terms given in Proposition 2.8, we can provide simple but elegant expressions for the basis functions $\varphi_i, \phi_i, i = 1, 2, 3$ of $(T, \tilde{S}_2(T), \tilde{\Sigma}_T)$. The proof of the following theorem follows the same argument of Theorem 2.6 and it is omitted

**Theorem 2.11.** The basis functions of $(T, \tilde{S}_2(T), \tilde{\Sigma}_T)$ satisfying (13) and (14) have the following expressions

$$ \varphi_i = \lambda_i - \sum_{j=1}^{3} \frac{e_j}{2\gamma_j} (1 - \delta_{ij}), \quad i = 1, 2, 3, \quad (30) $$

$$ \phi_i = \frac{e_i}{\gamma_i}, \quad i = 1, 2, 3. \quad (31) $$

## 2.3. An explicit error representation

For the enriched space $\tilde{S}_2(T)$, we introduce the approximation operator

$$ \begin{aligned} \Pi : C(T) & \rightarrow \tilde{S}_2(T) \\ f & \mapsto \sum_{i=1}^{3} L_i(f) \varphi_i + \sum_{i=1}^{3} L_i(f) \phi_i, \end{aligned} \quad (32) $$

where $\varphi_i, \phi_i, i = 1, 2, 3$, are the basis functions defined in Theorem 2.11. We are interested in evaluating or estimating the approximation error

$$ E[f] := f - \Pi[f]. \quad (33) $$

The following result shows that the error (33) can be decomposed in two parts: the first one is related to the linear triangular element while the second one depends on the enriched functions $e_i, i = 1, 2, 3$. To short the notation, in this representation we make use of a generalization of the classical trapezoidal formula to the case of line integrals. More precisely, for each $j = 1, 2, 3$, we set

$$ \mathcal{L}_j := \frac{1}{2} \sum_{\substack{i=1 \\ i \neq j}}^{3} L_i \quad (34) $$

and

$$ \mathcal{E}_j^{\text{tra}} := \mathcal{L}_j - I_j, \quad (35) $$

where $L_i$ and $I_j, i, j = 1, 2, 3$ are defined in (1) and (2), respectively.
**Proposition 2.12.** Let us assume that the nonvanishing conditions (21) hold and let $e_i$, $i = 1, 2, 3$, be the functions defined in (29). Then, for all $f \in C(T)$, the approximation error at any point $\mathbf{x} \in T$ is given by

$$E[f](\mathbf{x}) = E^{\text{lin}}[f](\mathbf{x}) + E^{\text{enr}}[f](\mathbf{x}), \qquad (36)$$

where

$$E^{\text{lin}}[f](\mathbf{x}) = f(\mathbf{x}) - \sum_{i=1}^{3} L_i(f) \lambda_i(\mathbf{x}), \qquad (37)$$

and

$$E^{\text{enr}}[f](\mathbf{x}) = \sum_{j=1}^{3} \frac{e_j(\mathbf{x})}{\gamma_j} \varepsilon_j^{\text{tra}}(f). \qquad (38)$$

**Proof.** By (33) and (32)

$$E[f] = f - \sum_{i=1}^{3} L_i(f) \varphi_i - \sum_{i=1}^{3} I_i(f) \phi_i.$$

By (30) and by changing the order of the summation, we get

$$\begin{aligned} \sum_{i=1}^{3} L_i(f) \varphi_i &= \sum_{i=1}^{3} L_i(f) \left( \lambda_i - \sum_{j=1}^{3} \frac{e_j}{2\gamma_j} (1 - \delta_{ij}) \right), \\ &= \sum_{i=1}^{3} L_i(f) \lambda_i - \sum_{j=1}^{3} \frac{e_j}{2\gamma_j} \sum_{i=1}^{3} (1 - \delta_{ij}) L_i(f), \\ &= \sum_{i=1}^{3} L_i(f) \lambda_i - \sum_{j=1, i \neq j}^{3} \frac{e_j}{2\gamma_j} \sum_{i=1, i \neq j}^{3} L_i(f), \\ &= \sum_{i=1}^{3} L_i(f) \lambda_i - \sum_{j=1, i=1}^{3} \frac{e_j}{\gamma_j} \mathcal{L}_j(f). \end{aligned}$$

Therefore, for all $\mathbf{x} \in T$, we get

$$\begin{aligned} E[f](\mathbf{x}) &= f(\mathbf{x}) - \sum_{i=1}^{3} L_i(f) \varphi_i(\mathbf{x}) - \sum_{j=1}^{3} I_j(f) \phi_j(\mathbf{x}) \\ &= f(\mathbf{x}) - \sum_{i=1}^{3} L_i(f) \lambda_i(\mathbf{x}) + \sum_{j=1, i=1}^{3} \frac{e_j(\mathbf{x})}{\gamma_j} \mathcal{L}_j(f) - \sum_{j=1}^{3} I_j(f) \phi_j(\mathbf{x}) \\ &= E^{\text{lin}}[f](\mathbf{x}) + \sum_{j=1}^{3} \frac{e_j(\mathbf{x})}{\gamma_j} (\mathcal{L}_j(f) - I_j(f)), \end{aligned}$$

as desired. $\square$

# 3. Non-polynomial enrichment of the standard triangular linear element

## 3.1. Non-polynomial enrichment

In this Section we introduce a more general enrichment of the standard triangular linear element $P_1$ based on three linearly independent arbitrary enrichment functions $e_1, e_2, e_3$. As before, we assume that these functions satisfy the vanishing conditions on the vertices

$$e_1(\mathbf{v}_i) = e_2(\mathbf{v}_i) = e_3(\mathbf{v}_i) = 0, \quad i = 1, 2, 3. \qquad (39)$$

The new element, denoted by GF3, is
$$ (T, G_2(T), \tilde{\Sigma}_T), $$

where

$$ G_2(T) = P_1(T) \oplus \operatorname{span}\{e_1, e_2, e_3\}. \qquad (40) $$

The motivation for the introduction of this new enriched element lies in the possibility to capture, through the new enrichment functions, features of the function to be approximated that cannot be accurately captured by previously considered basis. It is worthwhile to note that, in the present situation, we cannot apply previously developed approaches, based on the use of Simpson's rule or 1-point Gauss quadrature rule. The following theorem gives necessary and sufficient conditions on the enrichment functions $e_1, e_2, e_3$, such that $\tilde{\Sigma}_T$ is $G_2(T)$-unisolvent.

**Theorem 3.1.** Let be

$$ N = \begin{bmatrix} I_1(e_1) & I_1(e_2) & I_1(e_3) \\ I_2(e_1) & I_2(e_2) & I_2(e_3) \\ I_3(e_1) & I_3(e_2) & I_3(e_3) \end{bmatrix}, \qquad (41) $$

then the triple $(T, G_2(T), \tilde{\Sigma}_T)$ is a unisolvent element if and only if

$$ \det(N) \neq 0. \qquad (42) $$

**Proof.** Let assume that $\det(N) \neq 0$. Let $f \in G_2(T)$ be a function satisfying

$$ L_i(f) = f(v_i) = 0, \quad i = 1, 2, 3, \qquad (43) $$

$$ I_i(f) = \frac{1}{|\Gamma_i|} \int_{\Gamma_i} f d\sigma = 0, \quad i = 1, 2, 3. \qquad (44) $$

By (40), $f$ can be decomposed into the sum of a linear polynomial $p \in P_1(T)$ and an enriched part, that is

$$ f = p + \alpha_1 e_1 + \alpha_2 e_2 + \alpha_3 e_3, \quad \alpha_i \in \mathbb{R}, \quad i = 1, 2, 3. $$

The vanishing conditions (39) for the enrichment functions imply that

$$ p(v_i) = 0, \quad i = 1, 2, 3, $$

and therefore, since $p$ is linear, $p \equiv 0$ and $f$ coincides with its enriched part

$$ f = \alpha_1 e_1 + \alpha_2 e_2 + \alpha_3 e_3. $$

Consequently, equations (43) and (44) imply that

$$ \begin{bmatrix} I_1(e_1) & I_1(e_2) & I_1(e_3) \\ I_2(e_1) & I_2(e_2) & I_2(e_3) \\ I_3(e_1) & I_3(e_2) & I_3(e_3) \end{bmatrix} \begin{bmatrix} \alpha_1 \\ \alpha_2 \\ \alpha_3 \end{bmatrix} = \begin{bmatrix} 0 \\ 0 \\ 0 \end{bmatrix}. \qquad (45) $$

As the determinant of $N$ is non zero, this system has the unique solution $\alpha_1 = \alpha_2 = \alpha_3 = 0$. Hence $f = 0$.

Now, let assume that $\det(N) = 0$. Therefore, there exist three real numbers $\gamma_1, \gamma_2, \gamma_3$, not all zero, for which the function

$$ e := \sum_{i=1}^{3} \gamma_i e_i $$

satisfies

$$ I_j(e) = 0, \quad j = 1, 2, 3. $$

Moreover the vanishing conditions (39) imply that

$$ L_j(e) = 0, \quad j = 1, 2, 3, $$

therefore, we can exhibit a linear combination of the basis functions of $G_2(T)$ with coefficient not all zero in which all degrees of freedom vanish. Then $\tilde{\Sigma}_T$ is not $G_2(T)$-unisolvent. $\square$

There exists a large class of admissible enrichment functions, as the following example shows.
**Example 3.2.** Let us consider the three functions

$$ 
\begin{aligned}
e_1 &= (1 - \lambda_1)^{\alpha_1 - 1} \lambda_2^{\beta_1 - 1} \lambda_3^{\gamma_1 - 1}, \\
e_2 &= (1 - \lambda_2)^{\alpha_2 - 1} \lambda_1^{\beta_2 - 1} \lambda_3^{\gamma_2 - 1}, \\
e_3 &= (1 - \lambda_3)^{\alpha_3 - 1} \lambda_1^{\beta_3 - 1} \lambda_2^{\gamma_3 - 1},
\end{aligned} 
\quad (46) $$

with $\alpha_i, \beta_i, \gamma_i > 1, i = 1, 2, 3$. It is easy to see that $e_1, e_2, e_3$ satisfy conditions (39). Moreover, by using (27) we get

$$ I_i(e_j) = \frac{1}{|\Gamma_i|} \int_{\Gamma_i} e_j d\sigma = \delta_{ij} B(\beta_i, \gamma_i), \quad i, j = 1, 2, 3 \quad (47) $$

and then the matrix $N$ defined in (42) is a diagonal matrix with determinant different from zero. This implies that the functions $e_i, i = 1, 2, 3$, are linearly independent and therefore, by Theorem 3.1 we can enrich $P_1$ to the unisolvent element $(T, G_2(T), \Sigma_T)$ by using (46) as enrichment functions.

As in Section 2, we can obtain simple explicit expressions for the basis $\{\varphi_i, \phi_i, i = 1, 2, 3\}$ of $G_2(T)$, associated to a unisolvent GF3 element, which satisfy (13), (14). To this aim we denote by

$$ N^{-1} = [c_1 \ c_2 \ c_3], \quad (48) $$

the inverse matrix of $N$, where $c_i \in \mathbb{R}^3, i = 1, 2, 3$, are column vectors.

**Theorem 3.3.** The basis functions $\varphi_i, \phi_i, i = 1, 2, 3$, of $(T, G_2(T), \Sigma_T)$ satisfying (13) and (14) have the following expressions

$$ \varphi_i = \lambda_i - \frac{1}{2} \sum_{\substack{j=1 \\ j \neq i}}^{3} \langle \mathbf{e}, \mathbf{c}_j \rangle, \quad i = 1, 2, 3, \quad (49) $$

$$ \phi_i = \langle \mathbf{e}, \mathbf{c}_i \rangle, \quad i = 1, 2, 3, \quad (50) $$

where

$$ \mathbf{e} = [e_1, e_2, e_3]^T. \quad (51) $$

**Proof.** We start by showing the validity of (49) in the case $i = 1$. We set

$$ I_j(\mathbf{e}) = [I_j(e_1), I_j(e_2), I_j(e_3)]^T, \quad j = 1, 2, 3 $$

hence, since $NN^{-1} = I$, we easily get

$$ \langle I_j(\mathbf{e}), \mathbf{c}_i \rangle = \delta_{ij}. \quad (52) $$

Moreover, by Lemma 2.1, we get

$$ I_j(\lambda_i) = \frac{1}{2}(1 - \delta_{ij}). \quad (53) $$

As an element of $G_2(T)$, $\varphi_1$ can be represented as

$$ \varphi_1 = p + \beta_1 e_1 + \beta_2 e_2 + \beta_3 e_3 = p + \langle \mathbf{e}, \beta \rangle, \quad (54) $$

where $p \in P_1(T)$ and $\beta = [\beta_1, \beta_2, \beta_3]^T \in \mathbb{R}^3$. By using (13) and the vanishing conditions (39) we have

$$ \lambda_1(\mathbf{v}_j) = \delta_{1j} = L_j(\varphi_1) = p(\mathbf{v}_j), \quad j = 1, 2, 3, $$

so that $p = \lambda_1$, since they are linear polynomials. Therefore (54) becomes

$$ \varphi_1 = \lambda_1 + \langle \mathbf{e}, \beta \rangle, \quad (55) $$

and by applying $I_j, j = 1, 2, 3$, to both members of (55), by (13) and (53), we get

$$ 0 = \frac{1}{2}(1 - \delta_{1j}) + \langle I_j(\mathbf{e}), \beta \rangle, \quad j = 1, 2, 3, $$

or, in matrix form,
$$ \begin{bmatrix} 0 \\ 0 \\ 0 \end{bmatrix} = \frac{1}{2} \begin{bmatrix} 0 \\ 1 \\ 1 \end{bmatrix} + N\beta. $$

Consequently

$$ \beta = -\frac{1}{2} N^{-1} \begin{bmatrix} 0 \\ 1 \\ 1 \end{bmatrix} = -\frac{1}{2} (c_2 + c_3), $$

which, substituted in (55), gives the required expression (49) for $\varphi_1$. Now we show the validity of (50) for $i = 1$. We proceed in analogy to previous case and then we set

$$ \phi_1 = q + \langle \mathbf{e}, \mathbf{\gamma} \rangle, $$

where $q \in P_1(T)$ and $[\gamma_1, \gamma_2, \gamma_3]^T \in \mathbb{R}^3$. Since $\phi_1(\mathbf{v}_j) = 0$ for $j = 1, 2, 3$, this function can be expressed as

$$ \phi_1 = \langle \mathbf{e}, \mathbf{\gamma} \rangle. \qquad (56) $$

By applying $I_i$, $i = 1, 2, 3$, to both members of (56), by (14), we get

$$ \delta_{i1} = \langle I_i(\mathbf{e}), \mathbf{\gamma} \rangle, \quad i = 1, 2, 3, $$

or, in matrix form,

$$ \begin{bmatrix} 1 \\ 0 \\ 0 \end{bmatrix} = N\mathbf{\gamma}. $$

Consequently

$$ \mathbf{\gamma} = N^{-1} \begin{bmatrix} 1 \\ 0 \\ 0 \end{bmatrix} = \mathbf{c}_1, $$

which, substituted in (56), gives the required expression (50) for $\phi_1$. The expression of the other functions can be obtained using symmetry arguments. $\square$

## 3.2. An explicit error representation

In analogy to the case of polynomial enrichment described in previous section, we are interested in evaluating or estimating the error

$$ E[f] = f - \Pi[f] \qquad (57) $$

of the approximation operator

$$ \Pi[f] = \sum_{i=1}^{3} L_i(f)\varphi_i + \sum_{i=1}^{3} I_i(f)\phi_i, \qquad (58) $$

where the basis functions $\varphi_i$, $\phi_i$, $i = 1, 2, 3$ are now given as in (49) and (50). As before, we start by proving a decomposition of the error $E[f]$ as a sum of the error of the (non-enriched) linear triangular element plus an additional term which depends both on the enrichment functions $e_i$, $i = 1, 2, 3$ and the error (35) of the generalization of the classical trapezoidal formula to the case of line integrals (34). This representation will play an important role in the derivation of explicit bounds for the error in $L^1$-norm.

**Proposition 3.4.** Let $e_i \in C(T)$, $i = 1, 2, 3$ be enrichment functions such that $\tilde{\Sigma}_T$ is $G_2(T)$-unisolvent. Then, for all $f \in C(T)$, the approximation error $E[f]$ at any point $\mathbf{x} \in T$ is given by

$$ E[f](\mathbf{x}) = E^{\text{lin}}[f](\mathbf{x}) + \sum_{j=1}^{3} \langle \mathbf{e}(\mathbf{x}), \mathbf{c}_j \rangle \varepsilon_j^{\text{tra}}(f), \qquad (59) $$

where $E^{\text{lin}}[f](\mathbf{x})$, $\varepsilon_j^{\text{tra}}(f)$, $\mathbf{e}(\mathbf{x})$ and $\mathbf{c}_j$, $j = 1, 2, 3$ are defined as in (37), (35), (51) and (48), respectively.
**Proof.** By (57) and (58)

$$E[f] = f - \sum_{i=1}^{3} L_i(f)\varphi_i - \sum_{i=1}^{3} I_i(f)\phi_i.$$

By (49) and by changing the order of the summation, we get

$$\begin{aligned}
\sum_{i=1}^{3} L_i(f)\varphi_i &= \sum_{i=1}^{3} L_i(f) \left( \lambda_i - \frac{1}{2} \sum_{\substack{j=1 \\ j \neq i}}^{3} \langle \mathbf{e}, \mathbf{c}_j \rangle \right) \\
&= \sum_{i=1}^{3} L_i(f)\lambda_i - \sum_{i=1}^{3} L_i(f) \frac{1}{2} \sum_{j=1}^{3} \langle \mathbf{e}, \mathbf{c}_j \rangle (1 - \delta_{ij}) \\
&= \sum_{i=1}^{3} L_i(f)\lambda_i - \sum_{j=1}^{3} \langle \mathbf{e}, \mathbf{c}_j \rangle \frac{1}{2} \sum_{i=1}^{3} L_i(f)(1 - \delta_{ij}) \\
&= \sum_{i=1}^{3} L_i(f)\lambda_i - \sum_{j=1}^{3} \langle \mathbf{e}, \mathbf{c}_j \rangle \frac{1}{2} \sum_{\substack{i=1 \\ i \neq j}}^{3} L_i(f).
\end{aligned}$$

Consequently, for each $\mathbf{x} \in T$, we get

$$\begin{aligned}
E[f](\mathbf{x}) &= f(\mathbf{x}) - \sum_{i=1}^{3} L_i(f)\varphi_i(\mathbf{x}) - \sum_{j=1}^{3} I_j(f)\phi_j(\mathbf{x}) \\
&= f(\mathbf{x}) - \left( \sum_{i=1}^{3} L_i(f)\lambda_i(\mathbf{x}) - \sum_{j=1}^{3} \langle \mathbf{e}, \mathbf{c}_j \rangle \frac{1}{2} \sum_{\substack{i=1 \\ i \neq j}}^{3} L_i(f) \right) - \sum_{j=1}^{3} I_j(f) \langle \mathbf{e}(\mathbf{x}), \mathbf{c}_j \rangle \\
&= E^{\text{lin}}[f](\mathbf{x}) + \sum_{j=1}^{3} \langle \mathbf{e}(\mathbf{x}), \mathbf{c}_j \rangle \left( \frac{1}{2} \sum_{\substack{i=1 \\ i \neq j}}^{3} L_i(f) - I_j(f) \right) \\
&= E^{\text{lin}}[f](\mathbf{x}) + \sum_{j=1}^{3} \langle \mathbf{e}(\mathbf{x}), \mathbf{c}_j \rangle (\mathcal{L}_j(f) - I_j(f)),
\end{aligned}
\quad (60)$$

as required. □

**Remark 3.5.** It may be interesting to compare the proposed approximation operator $\Pi$ and the standard three node triangular approximation operator

$$\Pi^{\text{lin}}[f](\mathbf{x}) = \sum_{i=1}^{3} L_i(f)\lambda_i(\mathbf{x}). \qquad (61)$$

Using the representation of the error (60), the operator $\Pi$ can be formulated in an subtractive form more convenient for practical computation

$$\Pi[f](\mathbf{x}) = \Pi^{\text{lin}}[f](\mathbf{x}) - E^{\text{enr}}[f](\mathbf{x}), \quad (62)$$

where

$$E^{\text{enr}}[f](\mathbf{x}) := \sum_{j=1}^{3} \langle \mathbf{e}(\mathbf{x}), \mathbf{c}_j \rangle (\mathcal{L}_j(f) - I_j(f)).$$

Indeed, as shown in equation (62), the operator $\Pi$ may be computed by simply subtracting the enriched approximation operator $E^{\text{enr}}$ from the operator $\Pi^{\text{lin}}$, so that the two contributions can be evaluated separately.
### 3.3. Error bounds

The decomposition (59) is the key result to get the estimate of the error $E[f]$ in the case of a particular class of functions with continuous gradient. As usually, we say that $f$ is continuously differentiable on $T$ if it is continuously differentiable on an open set containing $T$. Other useful terminology and notations are clarified in the following.

**Definition 3.6.** A differentiable function $f$ is said to have a Lipschitz continuous gradient on $T$, if there exists a constant $\rho > 0$ such that

$$ \|\nabla f(\mathbf{x}) - \nabla f(\mathbf{y})\| \le \rho \|\mathbf{x} - \mathbf{y}\|, \quad \forall \mathbf{x}, \mathbf{y} \in T. \qquad (63) $$

By $C^{1,1}(\Omega)$ we denote the subclass of all functions $f$ which are continuously differentiable with Lipschitz continuous gradient on $T$. We call the smallest possible $\rho$ such that (63) holds Lipschitz constant for $\nabla f$ and we denote it by $L(\nabla f)$. The following result (see [17, Theorem 2.3]) will be useful in the following.

**Theorem 3.7.** Let $A: C^1(T) \to C(T)$ be a linear operator. The following statements are equivalent:

(i) for every convex function $g \in C^1(T)$, we have

$$ g(\mathbf{x}) \le A[g](\mathbf{x}), \quad \mathbf{x} \in T; \qquad (64) $$

(ii) for every $f \in C^{1,1}(T)$ we have

$$ |f(\mathbf{x}) - A[f](\mathbf{x})| \le \frac{L(\nabla f)}{2} (A[\|\cdot\|^2](\mathbf{x}) - \|\mathbf{x}\|^2), \quad \mathbf{x} \in T. \qquad (65) $$

Equality is attained for all functions of the form

$$ f(\mathbf{x}) = a(\mathbf{x}) + c \|\mathbf{x}\|^2, $$

where $c \in \mathbb{R}$ and $a(\mathbf{x})$ is any affine function.

**Remark 3.8.** We notice that the results of Theorem 3.7 hold true, with the needed changes, in the case of standard simplex in $\mathbb{R}^d$, $d \in \mathbb{N}$ (see [17, Theorem 2.3]).

Since each $\mathbf{x} \in T$ can be expressed as $\mathbf{x} = \sum_{i=1}^{3} \lambda_i(\mathbf{x})\mathbf{v}_i$, then for every convex function $f$ on $T$ we have

$$ f(\mathbf{x}) \le \sum_{i=1}^{3} \lambda_i(\mathbf{x}) f(\mathbf{v}_i) := \Pi_1[f](\mathbf{x}), \quad \mathbf{x} \in T, \qquad (66) $$

that is, the linear interpolation operator $\Pi_1$ satisfies condition (64), i.e. it approximates every convex function from above. Then, the following result holds.

**Theorem 3.9.** For every $f \in C^{1,1}(T)$ we have

$$ |E^{\text{lin}}[f](\mathbf{x})| = |f - \Pi_1[f]| \le \frac{L(\nabla f)}{2} \sum_{i=1}^{3} \lambda_i(\mathbf{x}) \|\mathbf{x} - \mathbf{v}_i\|^2, \quad \mathbf{x} \in T. $$

Equality is attained for all functions $f$ of the form

$$ f(\mathbf{x}) := a(\mathbf{x}) + c \|\mathbf{x}\|^2, $$

where $c \in \mathbb{R}$ and $a(\mathbf{x})$ is any affine function.

**Proof.** From (66) we can apply Theorem 3.7 to the linear operator $\Pi_1$ so that (65) becomes

$$ |f(\mathbf{x}) - \Pi_1[f](\mathbf{x})| \le \frac{L(\nabla f)}{2} \left( \sum_{i=1}^{3} \lambda_i(\mathbf{x}) \|\mathbf{v}_i\|^2 - \|\mathbf{x}\|^2 \right), \quad \mathbf{x} \in T. $$

It remains to show that
$$ \sum_{i=1}^{3} \lambda_{i}(\mathbf{x}) \|\mathbf{v}_{i}\|^{2} - \|\mathbf{x}\|^{2} = \sum_{i=1}^{3} \lambda_{i}(\mathbf{x}) \|\mathbf{x} - \mathbf{v}_{i}\|^{2}. $$

Indeed, we have

$$ \|\mathbf{x} - \mathbf{v}_i\|^2 = \|\mathbf{x}\|^2 - 2 \langle \mathbf{x}, \mathbf{v}_i \rangle + \|\mathbf{v}_i\|^2, \quad i = 1, 2, 3. $$

By multiplying each of the above equalities by $\lambda_i(\mathbf{x})$ and summing over all $i = 1, 2, 3$, we immediately get

$$ \sum_{i=1}^{3} \lambda_{i}(\mathbf{x}) \|\mathbf{x} - \mathbf{v}_{i}\|^{2} = \sum_{i=1}^{3} \lambda_{i}(\mathbf{x}) \|\mathbf{x}\|^{2} - 2 \left\langle \mathbf{x}, \sum_{i=1}^{3} \lambda_{i}(\mathbf{x}) \mathbf{v}_{i} \right\rangle + \sum_{i=1}^{3} \lambda_{i}(\mathbf{x}) \|\mathbf{v}_{i}\|^{2}. $$

The desired result now follows from the partition of unity and the linear precision properties of the barycentric coordinates. The last statement of the theorem follows directly by Theorem 3.7. □

With reference to the formula (59), it remains to bound the error $\mathcal{E}_j^{\text{tra}}(f)$ defined in (35).

**Theorem 3.10.** For every $f \in C^{1,1}(T)$ we have

$$ \left| \mathcal{L}_j(f) - \frac{1}{|\Gamma_j|} \int_{\Gamma_j} f d\sigma \right| \le \frac{L(\nabla f)}{12} |\Gamma_j|^2, \quad j = 1, 2, 3. \qquad (67) $$

Equality in (67) is attained for all functions of the form

$$ f(\mathbf{x}) := a(\mathbf{x}) + c \|\mathbf{x}\|^2, \qquad (68) $$

where $c \in \mathbb{R}$ and $a(\mathbf{x})$ is any affine function.

**Proof.** We prove (67) in the particular case $j = 1$ since the remaining cases can be proved by analogy. Let us denote by $\tilde{f}$ the map

$$ \tilde{f}(t) = f((1-t)\mathbf{v}_2 + t\mathbf{v}_3), \quad t \in [0, 1] \qquad (69) $$

and by

$$ \tilde{\Pi}_1[\tilde{f}](t) = (1-t)\tilde{f}(0) + t\tilde{f}(1) \qquad (70) $$

its linear interpolant at the end points of the interval $[0, 1]$. Therefore, we have

$$ \begin{aligned} \mathcal{E}_1^{\text{tra}}[f] &= \frac{1}{2} (f(\mathbf{v}_2) + f(\mathbf{v}_3)) - \frac{1}{|\Gamma_1|} \int_{\Gamma_1} f d\sigma \\ &= \frac{1}{2} (\tilde{f}(0) + \tilde{f}(1)) - \int_0^1 \tilde{f}(t) dt \\ &= \int_0^1 (\tilde{\Pi}_1[\tilde{f}](t) - \tilde{f}(t)) dt. \end{aligned} \qquad (71) $$

Consequently

$$ |\mathcal{E}_1^{\text{tra}}[f]| \le \int_0^1 |\tilde{\Pi}_1[\tilde{f}](t) - \tilde{f}(t)| dt. \qquad (72) $$

Since $\tilde{\Pi}_1$ approximates from above any convex function on $[0, 1]$, from Remark 3.8 we have

$$ |\tilde{\Pi}_1[\tilde{f}](t) - \tilde{f}(t)| \le \frac{L(\tilde{f}')}{2} (\tilde{\Pi}_1[|\cdot|^2](t) - |t|^2) = \frac{L(\tilde{f}')}{2} (t - t^2), $$

and therefore, from (72) we get
$$|\varepsilon_1^{\text{tra}}[f]| \le \frac{L(\tilde{f}')}{12}. \qquad (73)$$

Moreover, for each $s, t \in [0, 1]$, we have

$$\begin{aligned} |\tilde{f}'(s) - \tilde{f}'(t)| &= |\langle \nabla f(s\mathbf{v}_2 + (1-s)\mathbf{v}_3) - \nabla f(t\mathbf{v}_2 + (1-t)\mathbf{v}_3), \mathbf{v}_2 - \mathbf{v}_3 \rangle| \\ &\le \|\nabla f(s\mathbf{v}_2 + (1-s)\mathbf{v}_3) - \nabla f(t\mathbf{v}_2 + (1-t)\mathbf{v}_3)\| \|\mathbf{v}_2 - \mathbf{v}_3\| \\ &\le L(\nabla f) \|\mathbf{v}_2 - \mathbf{v}_3\|^2 |s-t|, \end{aligned} \qquad (74)$$

and then, by definition of Lipschitz constant, we get

$$L(\tilde{f}') \le L(\nabla f) \|\mathbf{v}_2 - \mathbf{v}_3\|^2 = L(\nabla f)|\Gamma_1|^2 \qquad (75)$$

which concludes the proof of the inequality (67). Let assume that the function $f(\mathbf{x})$ has the form (68) with $c = 0$. Therefore the function $\tilde{f}$ defined in (69) is affine and then $\tilde{\Pi}_1[\tilde{f}] = \tilde{f}$ and $L(\tilde{f}') = 0$. Consequently, (67) holds with equality in this case. Now let assume that $f(\mathbf{x}) = \|\mathbf{x}\|^2$. In this case the function in (69)

$$\tilde{f}(t) = \|t\mathbf{v}_2 + (1-t)\mathbf{v}_3\|^2, \quad t \in [0, 1],$$

is a univariate quadratic polynomial satisfying

$$\tilde{f}''(t) = 2 \|\mathbf{v}_2 - \mathbf{v}_3\|^2. \qquad (76)$$

We expand $\tilde{f}(0)$ and $\tilde{f}(1)$ in Taylor series centered in $t$ and we get, from (76),

$$\tilde{f}(0) = \tilde{f}(t) - t\tilde{f}'(t) + t^2 \|\mathbf{v}_2 - \mathbf{v}_3\|^2, \qquad (77)$$

$$\tilde{f}(1) = \tilde{f}(t) + (1-t)\tilde{f}'(t) + (1-t)^2 \|\mathbf{v}_2 - \mathbf{v}_3\|^2. \qquad (78)$$

Multiplying (77) by $1-t$, (78) by $t$ and summing yields

$$(1-t)\tilde{f}(0) + t\tilde{f}(1) = \tilde{f}(t) + ((1-t)t^2 + t(1-t)^2) \|\mathbf{v}_2 - \mathbf{v}_3\|^2, \qquad (79)$$

or equivalently

$$\tilde{\Pi}_1[\tilde{f}](t) - \tilde{f}(t) = ((1-t)t^2 + t(1-t)^2) \|\mathbf{v}_2 - \mathbf{v}_3\|^2. \qquad (80)$$

Hence, by integrating (80), we get from (71)

$$\begin{aligned} \varepsilon_1^{\text{tra}}[\|\cdot\|^2] &= \|\mathbf{v}_2 - \mathbf{v}_3\|^2 \int_0^1 \left((1-t)t^2 + t(1-t)^2\right) dt, \\ &= \|\mathbf{v}_2 - \mathbf{v}_3\|^2 (B(3, 2) + B(2, 3)), \\ &= \frac{1}{6} \|\mathbf{v}_2 - \mathbf{v}_3\|^2, \end{aligned}$$

which is exactly the term on the right hand side in (67), since $L(\nabla (\|\cdot\|^2)) = 2$. Finally, if $f(\mathbf{x})$ has the general form (68), then the equality (67) easily follows. □

Combining Proposition 3.4, Theorem 3.9 and Theorem 3.10, we arrive at the main error estimate.

**Theorem 3.11.** For any $v \in C^{1,1}(T)$, the following explicit error estimate holds

$$|f(\mathbf{x}) - \Pi[f](\mathbf{x})| \le \frac{L(\nabla f)}{2} \left( \sum_{i=1}^{3} \lambda_i(\mathbf{x}) \|\mathbf{x} - \mathbf{v}_i\|^2 + \frac{1}{6} \sum_{i=1}^{3} |\Gamma_i|^2 |\langle c_i, e(\mathbf{x}) \rangle| \right), \quad \mathbf{x} \in T. \qquad (81)$$
### 3.4. The $L^\infty$ error estimate

The bound of the point-wise error given in Theorem 3.11 allows us to get a bound in $L^\infty$-norm

$$ \|f\|_{L^\infty(T)} = \max_{\mathbf{x} \in T} |f(\mathbf{x})| $$

of the error $E[f]$ defined in (57). As shown in the following result, this bound is proportional to the square of the radius of the circumcircle of the triangle $T$, by a constant factor which depends on $f$ and on the enrichment functions $e_i$, $i = 1, 2, 3$.

**Corollary 3.12.** For any $f \in C^{1,1}(T)$, the following explicit error estimate holds

$$ \|f - \Pi[f]\|_{L^\infty(T)} \le \frac{L(\nabla f)}{2} \left( 1 + \frac{3}{2} \max_{i=1,2,3} \|\langle \mathbf{c}_i, \mathbf{e} \rangle\|_{L^\infty(T)} \right) R^2, \quad (82) $$

where $\mathbf{c}_i$, $i = 1, 2, 3$ and $\mathbf{e}$ are defined in (48) and (51), respectively, and $R$ is the circumradius of $T$.

**Proof.** We bound the two terms inside the bracket appearing on the right hand side of equation (81) separately. With regard to the first term, we show that

$$ \sup_{\mathbf{x} \in T} \sum_{i=1}^{3} \lambda_i(\mathbf{x}) \|\mathbf{x} - \mathbf{v}_i\|^2 \le R^2. \quad (83) $$

Indeed, if $\mathbf{c}$ is the circumcenter of $T$, by using the linear precision and partition of unity properties of barycentric coordinates, we have

$$ \begin{aligned} \sum_{i=1}^{3} \lambda_i(\mathbf{x}) \|\mathbf{x} - \mathbf{v}_i\|^2 &= \sum_{i=1}^{3} \lambda_i(\mathbf{x}) \|\mathbf{x} - \mathbf{c} - (\mathbf{v}_i - \mathbf{c})\|^2 \\ &= \sum_{i=1}^{3} \lambda_i(\mathbf{x}) \left( \|\mathbf{x} - \mathbf{c}\|^2 - 2 \langle \mathbf{x} - \mathbf{c}, \mathbf{v}_i - \mathbf{c} \rangle + \|\mathbf{v}_i - \mathbf{c}\|^2 \right) \\ &= \|\mathbf{x} - \mathbf{c}\|^2 - 2 \left\langle \mathbf{x} - \mathbf{c}, \sum_{i=1}^{3} \lambda_i(\mathbf{x}) \mathbf{v}_i - \mathbf{c} \right\rangle + \sum_{i=1}^{3} \lambda_i(\mathbf{x}) \|\mathbf{v}_i - \mathbf{c}\|^2 \\ &= \|\mathbf{x} - \mathbf{c}\|^2 - 2 \langle \mathbf{x} - \mathbf{c}, \mathbf{x} - \mathbf{c} \rangle + R^2 \\ &= R^2 - \|\mathbf{x} - \mathbf{c}\|^2. \end{aligned} $$

Therefore, for each $\mathbf{x} \in T$

$$ \sum_{i=1}^{3} \lambda_i(\mathbf{x}) \|\mathbf{x} - \mathbf{v}_i\|^2 \le R^2 - \min_{\mathbf{x} \in T} \|\mathbf{x} - \mathbf{c}\|^2 \le R^2, \quad \mathbf{x} \in T, \quad (84) $$

and so (83) is valid. With regard to the second term inside the bracket, we use Leibniz's inequality to bound the sum of the squares of edge lengths of the triangle in terms of its circumradius [24], i.e.

$$ \sum_{i=1}^{3} |\Gamma_i|^2 \le 9R^2. \quad (85) $$

Combining the last estimate with (83) gives us the desired inequality (82). □

**Remark 3.13.** We notice that, for $\mathbf{x} = \mathbf{c}$ we get

$$ \sum_{i=1}^{3} \lambda_i(\mathbf{c}) \|\mathbf{c} - \mathbf{v}_i\|^2 = R^2, $$

then inequality (84) holds with equality if and only if the circumcenter $\mathbf{c}$ belongs to $T$. This property characterizes acute triangles, i.e. triangles having all angles less than $\pi/2$.

**Remark 3.14.** Let $X$ be a set of discrete points in a general position, e.g. a set of scattered points, and assume that we need to triangulate $X$. Inequality (82) suggests us the use of a triangulation $\mathcal{T}$ which minimizes the maximum of the squares of the circumradii over all triangles of $\mathcal{T}$. The Delaunay triangulation has this property.
**Remark 3.15.** For any triangle $T$, it is possible to obtain a more precise $L^\infty$-error bound in terms of the square of the circum-radius minus a nonnegative quantity. Indeed, denoting by $\mathbf{c}, \mathbf{b}$ and $R$ the circumcenter, the barycenter and the circumradius of $T$, respectively, then the following bound holds, see [22, Theorem 2.4.4]

$$ \frac{1}{9} \sum_{i=1}^{3} |\Gamma_i|^2 = R^2 - \|\mathbf{b} - \mathbf{c}\|^2 \le R^2 - \min_{\mathbf{x} \in T} \|\mathbf{x} - \mathbf{c}\|^2. $$

The last step follows because $\mathbf{b}$ belong to $T$. Therefore, in Corollary 3.12, $R^2$ can be replaced by the smaller value $R^2 - \min_{\mathbf{x} \in T} \|\mathbf{c} - \mathbf{x}\|^2$.

**Remark 3.16.** Let $\Gamma_3$ be the longest side of $T$ and let $\mathbf{v}_{12}$ be its midpoint. The line segment joining $\mathbf{c}$ with $\mathbf{v}_{12}$ is a perpendicular bisector of $\Gamma_3$. Denoting by $h = \sup_{\mathbf{v}, \mathbf{w} \in T} \|\mathbf{v} - \mathbf{w}\|$ the diameter of $T$, by the Pythagorean Theorem we get

$$ \left(\frac{h}{2}\right)^2 + \min_{\mathbf{x} \in T} \|\mathbf{x} - \mathbf{c}\|^2 = R^2, $$

and then

$$ R^2 - \min_{\mathbf{x} \in T} \|\mathbf{x} - \mathbf{c}\|^2 = \frac{1}{4}h^2. $$

From Remark 3.15 and Remark 3.16, the Corollary 3.12 becomes

**Corollary 3.17.** For any $f \in C^{1,1}(T)$, the following explicit error estimate holds

$$ \|f - \Pi[f]\|_{L^{\infty}(T)} \le \frac{L(\nabla f)}{8} \left( 1 + \frac{3}{2} \max_{i=1,2,3} \|\langle \mathbf{c}_i, \mathbf{e} \rangle\|_{L^{\infty}(T)} \right) h^2, \quad (86) $$

where $h$ is the diameter of $T$.

### 3.5. The $L^1$ error estimate

The bound of the point-wise error given in Theorem 3.11 allows us to get a bound in $L^1$-norm

$$ \|f\|_{L^1(T)} = \int_T |f(\mathbf{x})| d\mathbf{x} $$

of the error $E[f]$ defined in (57).

**Theorem 3.18.** For any $f \in C^{1,1}(T)$, the following explicit error estimate holds

$$ \|f - \Pi[f]\|_{L^1(T)} \le \frac{L(\nabla f)}{24} (1 + e_{max}) |T| \omega(T), \quad (87) $$

where $|T|$ is the area of the triangle $T$,

$$ e_{max} := \max_{i=1,2,3} \frac{2}{|T|} \|\langle \mathbf{c}_i, \mathbf{e} \rangle\|_{L^1(T)} \quad (88) $$

and

$$ \omega(T) := \sum_{i=1}^{3} |\Gamma_i|^2. $$

**Proof.** We start by proving the identity

$$ \int_T \sum_{i=1}^{3} \lambda_i(\mathbf{x}) \|\mathbf{x} - \mathbf{v}_i\|^2 d\mathbf{x} = \frac{|T|}{12} \sum_{i=1}^{3} |\Gamma_i|^2. \quad (89) $$

Since $\sum_{i=1}^3 \lambda_i(\mathbf{x}) \|\mathbf{x} - \mathbf{v}_i\|^2$ is a quadratic polynomial which vanishes at all vertices of $T$, by [20, Theorem 5.1], we get
$$ \int_T \sum_{i=1}^{3} \lambda_i(\mathbf{x}) \|\mathbf{x} - \mathbf{v}_i\|^2 d\mathbf{x} = \frac{|T|}{4} \sum_{i=1}^{3} \|\mathbf{b} - \mathbf{v}_i\|^2, \qquad (90) $$

where $\mathbf{b}$ is the barycenter of $T$. Moreover, since

$$ \|\mathbf{v}_i - \mathbf{v}_j\|^2 = \|\mathbf{v}_i\|^2 + \|\mathbf{v}_j\|^2 - 2\langle \mathbf{v}_i, \mathbf{v}_j \rangle, \quad i, j = 1, 2, 3, $$

$$ \|\mathbf{b} - \mathbf{v}_j\|^2 = \|\mathbf{b}\|^2 + \|\mathbf{v}_j\|^2 - 2\langle \mathbf{b}, \mathbf{v}_j \rangle, \quad j = 1, 2, 3, $$

$$ \|\mathbf{v}_i - \mathbf{v}_j\|^2 - \|\mathbf{b} - \mathbf{v}_j\|^2 = \|\mathbf{v}_i\|^2 - \|\mathbf{b}\|^2 + 2\langle \mathbf{b} - \mathbf{v}_i, \mathbf{v}_j \rangle \quad i, j = 1, 2, 3, $$

and then, by summing over all $j = 1, 2, 3$, both members of the last equality, we immediately get

$$ \begin{aligned} \sum_{j=1}^{3} \|\mathbf{v}_i - \mathbf{v}_j\|^2 - \sum_{j=1}^{3} \|\mathbf{b} - \mathbf{v}_j\|^2 &= 3 (\|\mathbf{v}_i\|^2 - \|\mathbf{b}\|^2) + 6 \langle \mathbf{b} - \mathbf{v}_i, \mathbf{b} \rangle \\ &= 3 (\|\mathbf{v}_i\|^2 - \|\mathbf{b}\|^2 + 2 \langle \mathbf{b} - \mathbf{v}_i, \mathbf{b} \rangle) \\ &= 3 \|\mathbf{b} - \mathbf{v}_i\|^2, \quad i = 1, 2, 3, \end{aligned} $$

or equivalently,

$$ \sum_{j=1}^{3} \|\mathbf{v}_i - \mathbf{v}_j\|^2 = \sum_{j=1}^{3} \|\mathbf{b} - \mathbf{v}_j\|^2 + 3 \|\mathbf{b} - \mathbf{v}_i\|^2, \quad i = 1, 2, 3. $$

Now we sum both members of the above equality over all $i = 1, 2, 3$ and we get

$$ \sum_{i=1}^{3} \sum_{j=1}^{3} \|\mathbf{v}_i - \mathbf{v}_j\|^2 = 6 \sum_{i=1}^{3} \|\mathbf{b} - \mathbf{v}_i\|^2, $$

or equivalently

$$ \sum_{i=1}^{3} \|\mathbf{b} - \mathbf{v}_i\|^2 = \frac{1}{6} \sum_{i=1}^{3} \sum_{j=1}^{3} \|\mathbf{v}_i - \mathbf{v}_j\|^2 = \frac{1}{3} \sum_{i=1}^{3} |\Gamma_i|^2. \qquad (91) $$

The identity (89) follows by substituting (91) in (90). Finally, by integrating (81) and by using (89), we get (87) and then the thesis. □

In the following, we further analyze the case of the enrichment functions are

$$ e_1 = \lambda_2^{\alpha_1-1} \lambda_3^{\beta_1-1}, \quad e_2 = \lambda_1^{\alpha_2-1} \lambda_3^{\beta_2-1}, \quad e_3 = \lambda_1^{\alpha_3-1} \lambda_2^{\beta_3-1}, \qquad (92) $$

with $\alpha_i, \beta_i > 1$. In particular, we determine the 'best' parameters $\alpha_i, \beta_i$, in order to minimize the approximation error (87).

**Theorem 3.19.** Let $e_i$, $i = 1, 2, 3$, be the enrichment functions defined in (92). Then, for any $f \in C^{1,1}(T)$, we get

$$ \|f - \Pi[f]\|_{L^1(T)} \le \frac{L(\nabla f)}{24} \left(1 + \frac{4}{\mu}\right) |T| \omega(T), \qquad (93) $$

where $\mu = \min_{i=1,2,3} (\alpha_i + \beta_i)$.

**Proof.** Let us compute $e_{max}$ defined in (88), relative to the enrichment functions (92). In this case, from (27), the matrix $N$ defined in (41) is a diagonal matrix and

$$ N^{-1} = [c_1, c_2, c_3] = \begin{bmatrix} \frac{1}{B(\alpha_1, \beta_1)} & 0 & 0 \\ 0 & \frac{1}{B(\alpha_2, \beta_2)} & 0 \\ 0 & 0 & \frac{1}{B(\alpha_3, \beta_3)} \end{bmatrix}. $$

Consequently we get
$$ 
\begin{aligned}
\frac{2}{|T|} \int_T |\langle \mathbf{c}_1, \mathbf{e}(\mathbf{x}) \rangle| d\mathbf{x} &= \frac{2}{|T| B(\alpha_1, \beta_1)} \int_T \lambda_2^{\alpha_1-1}(\mathbf{x}) \lambda_3^{\beta_1-1}(\mathbf{x}) d\mathbf{x}, \\
\frac{2}{|T|} \int_T |\langle \mathbf{c}_2, \mathbf{e}(\mathbf{x}) \rangle| d\mathbf{x} &= \frac{2}{|T| B(\alpha_2, \beta_2)} \int_T \lambda_1^{\alpha_2-1}(\mathbf{x}) \lambda_3^{\beta_2-1}(\mathbf{x}) d\mathbf{x}, \\
\frac{2}{|T|} \int_T |\langle \mathbf{c}_3, \mathbf{e}(\mathbf{x}) \rangle| d\mathbf{x} &= \frac{2}{|T| B(\alpha_3, \beta_3)} \int_T \lambda_1^{\alpha_3-1}(\mathbf{x}) \lambda_2^{\beta_3-1}(\mathbf{x}) d\mathbf{x}.
\end{aligned}
\quad (94) $$

In order to compute the integrals at the right side of (94) we denote by

$$ \hat{T} = \{ \hat{\mathbf{x}} = (\hat{x}, \hat{y}) \in \mathbb{R}^2, \hat{x} \ge 0, \hat{y} \ge 0, \hat{x} + \hat{y} \le 1 \}, \quad (95) $$

the triangle of vertices $\hat{\mathbf{v}}_1 = (1, 0)$, $\hat{\mathbf{v}}_2 = (0, 1)$ and $\hat{\mathbf{v}}_3 = (0, 0)$ and by $F: \hat{T} \to T$ the affine map defined as

$$ F(\hat{\mathbf{x}}) = \hat{\lambda}_1(\hat{\mathbf{x}}) \mathbf{v}_1 + \hat{\lambda}_2(\hat{\mathbf{x}}) \mathbf{v}_2 + \hat{\lambda}_3(\hat{\mathbf{x}}) \mathbf{v}_3, $$

where $\hat{\lambda}_1(\hat{\mathbf{x}})$, $\hat{\lambda}_2(\hat{\mathbf{x}})$, $\hat{\lambda}_3(\hat{\mathbf{x}})$ are the barycentric coordinates of $\hat{T}$. Since $F$ preserves the ratio of areas of two triangles, by definition of barycenter coordinates we get

$$ \lambda_i \circ F(\hat{\mathbf{x}}) = \hat{\lambda}_i(\hat{\mathbf{x}}). $$

Therefore, by using $F$ as change of variables, we get

$$ 
\begin{aligned}
\int_T \lambda_2^{\alpha_1-1}(\mathbf{x}) \lambda_3^{\beta_1-1}(\mathbf{x}) d\mathbf{x} &= 2|T| \int_{\hat{T}} \hat{\lambda}_2^{\alpha_1-1}(\hat{\mathbf{x}}) \hat{\lambda}_3^{\beta_1-1}(\hat{\mathbf{x}}) d\hat{\mathbf{x}}, \\
\int_T \lambda_1^{\alpha_2-1}(\mathbf{x}) \lambda_3^{\beta_2-1}(\mathbf{x}) d\mathbf{x} &= 2|T| \int_{\hat{T}} \hat{\lambda}_1^{\alpha_2-1}(\hat{\mathbf{x}}) \hat{\lambda}_3^{\beta_2-1}(\hat{\mathbf{x}}) d\hat{\mathbf{x}}, \\
\int_T \lambda_1^{\alpha_3-1}(\mathbf{x}) \lambda_2^{\beta_3-1}(\mathbf{x}) d\mathbf{x} &= 2|T| \int_{\hat{T}} \hat{\lambda}_1^{\alpha_3-1}(\hat{\mathbf{x}}) \hat{\lambda}_2^{\beta_3-1}(\hat{\mathbf{x}}) d\hat{\mathbf{x}},
\end{aligned}
\quad (96) $$

where $2|T|$ is the Jacobian determinant of $F$. The integrals at right member of equalities (96) can be computed by using the powerful formula presented in [5]

$$ \int_{\hat{T}} \hat{\lambda}_1^{\alpha}(\hat{\mathbf{x}}) \hat{\lambda}_2^{\beta}(\hat{\mathbf{x}}) \hat{\lambda}_3^{\gamma}(\hat{\mathbf{x}}) d\hat{\mathbf{x}} = \frac{\Gamma(\alpha + 1)\Gamma(\beta + 1)\Gamma(\gamma + 1)}{\Gamma(\alpha + \beta + \gamma + 3)}, \quad (97) $$

valid for each $\alpha, \beta, \gamma > -1$. By combining equalities (94), (96), (97) and (25), we get (93). □

Now we consider the more general enrichment functions, already introduced in Example 3.2,

$$ \tilde{e}_i := (1 - \lambda_i)^{\gamma_i} e_i, \quad i = 1, 2, 3, \quad \gamma_i \ge 0, \quad (98) $$

where $e_i$, $i = 1, 2, 3$, are defined in (92). In analogy to the Theorem 3.19, in Theorem 3.21 we determine, by using the enrichment functions (98), the 'best' parameters $\alpha_i$, $\beta_i$, $\gamma_i$ which minimize the approximation error (87). To this aim we consider the triangle $\hat{T}$ with vertices $\hat{\mathbf{v}}_1 = (1, 0)$, $\hat{\mathbf{v}}_2 = (0, 1)$ and $\hat{\mathbf{v}}_3 = (0, 0)$ and barycentric coordinates $\hat{\lambda}_1(\hat{\mathbf{x}})$, $\hat{\lambda}_2(\hat{\mathbf{x}})$, $\hat{\lambda}_3(\hat{\mathbf{x}})$.

**Lemma 3.20.** For each $i = 1, 2, 3$, the following equality holds

$$ \int_{\hat{T}} (1 - \hat{\lambda}_i(\hat{\mathbf{x}}))^{\gamma_i} \hat{\lambda}_1^{\alpha_1-1}(\hat{\mathbf{x}}) \hat{\lambda}_2^{\alpha_2-1}(\hat{\mathbf{x}}) \hat{\lambda}_3^{\alpha_3-1}(\hat{\mathbf{x}}) d\hat{\mathbf{x}} = \frac{\Gamma(\alpha_1)\Gamma(\alpha_2)\Gamma(\alpha_3)}{\Gamma(\sum_{j=1, j \neq i}^3 \alpha_j)} \mu_i, \quad (99) $$

where

$$ \mu_i := \frac{\Gamma(\gamma_i + \sum_{j=1, j \neq i}^3 \alpha_j)}{\Gamma(\gamma_i + \sum_{j=1}^3 \alpha_j)}, \quad \gamma_i \ge 0, \quad \alpha_i > 1, \quad i = 1, 2, 3. $$
**Proof.** It is well known that the barycentric coordinates of $\hat{T}$ in terms of Cartesian coordinates are

$$ \hat{\lambda}_1(\hat{\mathbf{x}}) = \hat{\mathbf{x}}, \quad \hat{\lambda}_2(\hat{\mathbf{x}}) = \hat{\mathbf{y}}, \quad \hat{\lambda}_3(\hat{\mathbf{x}}) = 1 - \hat{\mathbf{x}} - \hat{\mathbf{y}}. $$

We use them to prove the formula (99) for $i = 1$. We note that

$$ \int_{\hat{T}} (1 - \hat{x})^{\gamma_1 \hat{x}^{\alpha_1 - 1}} \hat{y}^{\alpha_2 - 1} (1 - \hat{x} - \hat{y})^{\alpha_3 - 1} d\hat{x} d\hat{y} = \int_{0}^{1} (1 - \hat{x})^{\gamma_1 \hat{x}^{\alpha_1 - 1}} I(\hat{x}) d\hat{x}, \quad (100) $$

where

$$ I(\hat{\mathbf{x}}) = \int_{0}^{1-\hat{x}} \hat{y}^{\alpha_2-1} (1-\hat{x}-\hat{y})^{\alpha_3-1} d\hat{y}. \quad (101) $$

In order to compute the integral (101) we consider the change of variable

$$ \hat{z} = \frac{\hat{y}}{1 - \hat{x}}, $$

and then

$$ I(\hat{\mathbf{x}}) = (1 - \hat{x})^{\alpha_2 + \alpha_3 - 1} \int_{0}^{1} \hat{z}^{\alpha_2 - 1} (1 - \hat{z})^{\alpha_3 - 1} d\hat{z} \quad (102) $$

$$ = (1 - \hat{x})^{\alpha_2 + \alpha_3 - 1} \frac{\Gamma(\alpha_2)\Gamma(\alpha_3)}{\Gamma(\alpha_2 + \alpha_3)}. \quad (103) $$

Finally, the right side of (100) becomes

$$ \int_{0}^{1} \hat{x}^{\alpha_1 - 1} (1 - \hat{x})^{\gamma_1 + \alpha_2 + \alpha_3 - 1} dx = \frac{\Gamma(\alpha_1)\Gamma(\gamma_1 + \alpha_2 + \alpha_3)}{\Gamma(\alpha_1 + \alpha_2 + \alpha_3 + \gamma_1)}. \quad (104) $$

The result follows by combining (100), (103) and (104). By the same arguments we can prove formula (99) for $i = 2$. In order to prove (99) for $i = 3$, we consider the coordinate transformations

$$ \pi : \hat{T} \to \hat{T}, \quad \pi(x, y) = (y, x), $$

$$ \theta : \hat{T} \to \hat{T}, \quad \theta(x, y) = (1 - x - y, y) =: (u, v), $$

and the change of variables $\pi \circ \theta$ whose Jacobian determinant is equal to 1. Then we get

$$ \int_{\hat{T}} (x + y)^{\gamma_1} x^{\alpha_1 - 1} y^{\alpha_2 - 1} (1 - x - y)^{\alpha_3 - 1} dxdy = \int_{\hat{T}} (1 - u)^{\gamma_1} G(u, v) dudv, $$

where $G(u, v) := u^{\alpha_3 - 1} v^{\alpha_2 - 1} (1 - u - v)^{\alpha_1 - 1}$. The result now follows from the case $i = 1$. $\square$

**Theorem 3.21.** Let $\tilde{e}_i$, $i = 1, 2, 3$, be the enrichment functions defined in (98). Then, for any $f \in C^{1,1}(T)$, we get

$$ \|f - \Pi[f]\|_{L^1(T)} \le \frac{L(\nabla f)}{24} \left(1 + \frac{4}{\mu}\right) |T| \omega(T), \quad (105) $$

where $\mu = \min_{i=1,2,3} (\alpha_i + \beta_i + \gamma_i)$.

**Proof.** Let us compute $e_{max}$ defined in (88), relative to the enrichment functions (98). By (47), the matrix $N$ defined in (41) is a diagonal matrix and

$$ N^{-1} = [c_1, c_2, c_3] = \begin{bmatrix} \frac{1}{B(\alpha_1, \beta_1)} & 0 & 0 \\ 0 & \frac{1}{B(\alpha_2, \beta_2)} & 0 \\ 0 & 0 & \frac{1}{B(\alpha_3, \beta_3)} \end{bmatrix}. $$

Consequently we get
$$ \frac{2}{|T|} \int_T |\langle \mathbf{c}_1, \mathbf{e}(\mathbf{x}) \rangle| d\mathbf{x} = \frac{2}{|T| B(\alpha_1, \beta_1)} \int_T (1 - \lambda_1(\mathbf{x}))^{\gamma_1} \lambda_2^{\alpha_1-1}(\mathbf{x}) \lambda_3^{\beta_1-1}(\mathbf{x}) d\mathbf{x}, $$

$$ \frac{2}{|T|} \int_T |\langle \mathbf{c}_2, \mathbf{e}(\mathbf{x}) \rangle| d\mathbf{x} = \frac{2}{|T| B(\alpha_2, \beta_2)} \int_T (1 - \lambda_2(\mathbf{x}))^{\gamma_2} \lambda_1^{\alpha_2-1}(\mathbf{x}) \lambda_3^{\beta_2-1}(\mathbf{x}) d\mathbf{x}, $$

$$ \frac{2}{|T|} \int_T |\langle \mathbf{c}_3, \mathbf{e}(\mathbf{x}) \rangle| d\mathbf{x} = \frac{2}{|T| B(\alpha_3, \beta_3)} \int_T (1 - \lambda_3(\mathbf{x}))^{\gamma_3} \lambda_1^{\alpha_3-1}(\mathbf{x}) \lambda_2^{\beta_3-1}(\mathbf{x}) d\mathbf{x}. $$

By using the same strategy of the proof of Theorem 3.19, we get

$$ \int_T (1 - \lambda_1(\mathbf{x}))^{\gamma_1} \lambda_2^{\alpha_1-1}(\mathbf{x}) \lambda_3^{\beta_1-1}(\mathbf{x}) d\mathbf{x} = 2|T| \int_{\hat{T}} (1 - \hat{\lambda}_1(\hat{\mathbf{x}}))^{\gamma_1} \hat{\lambda}_2^{\alpha_1-1}(\hat{\mathbf{x}}) \hat{\lambda}_3^{\beta_1-1}(\hat{\mathbf{x}}) d\hat{\mathbf{x}}, $$

$$ \int_T (1 - \lambda_2(\mathbf{x}))^{\gamma_2} \lambda_1^{\alpha_2-1}(\mathbf{x}) \lambda_3^{\beta_2-1}(\mathbf{x}) d\mathbf{x} = 2|T| \int_{\hat{T}} (1 - \hat{\lambda}_2(\hat{\mathbf{x}}))^{\gamma_2} \hat{\lambda}_1^{\alpha_2-1}(\hat{\mathbf{x}}) \hat{\lambda}_3^{\beta_2-1}(\hat{\mathbf{x}}) d\hat{\mathbf{x}}, $$

$$ \int_T (1 - \lambda_3(\mathbf{x}))^{\gamma_3} \lambda_1^{\alpha_3-1}(\mathbf{x}) \lambda_2^{\beta_3-1}(\mathbf{x}) d\mathbf{x} = 2|T| \int_{\hat{T}} (1 - \hat{\lambda}_3(\hat{\mathbf{x}}))^{\gamma_3} \hat{\lambda}_1^{\alpha_3-1}(\hat{\mathbf{x}}) \hat{\lambda}_2^{\beta_3-1}(\hat{\mathbf{x}}) d\hat{\mathbf{x}}. $$

The result follows from Lemma 3.20. □

The following result shows that the bound (105) is proportional to the fourth power of the circumradius of the triangle T.

**Corollary 3.22.** Let $e_i$, $i = 1, 2, 3$, be the enrichment functions defined in (98). Then, for any $f \in C^{1,1}(T)$, we get

$$ \|f - \Pi[f]\|_{L^1(T)} \le \frac{9\sqrt{3}L(\nabla f)}{32} \left(1 + \frac{4}{\mu}\right) R^4, $$

where R is the circumradius of T and $\mu = \min_{i=1,2,3} (\alpha_i + \beta_i + \gamma_i)$.

**Proof.** By Theorem 3.21 and Leibniz's inequality, it suffices to bound the area of T in terms of its circumradius. To this end, we use Weitzenböck's inequality [24], which states that for any triangle T of sides $\Gamma_1, \Gamma_2, \Gamma_3$, the following inequality holds

$$ 4\sqrt{3}|T| \le \sum_{i=1}^{3} |\Gamma_i|^2. $$

With the aid of this estimate and Leibniz's inequality, it is easily seen that

$$ |T| \le \frac{9}{4\sqrt{3}} R^2. \qquad (106) $$

The result follows by (85), (106) and (105). □

By Remark 3.15, it is possible to obtain a more precise $L^1$-error bound in terms of the diameter of T.

**Corollary 3.23.** Let $e_i$, $i = 1, 2, 3$, be the enrichment functions defined in (98). Then, for any $f \in C^{1,1}(T)$, we have

$$ \|f - \Pi[f]\|_{L^1(T)} \le \frac{9\sqrt{3}L(\nabla f)}{512} \left(1 + \frac{4}{\mu}\right) h^4, $$

where h is the diameter of T and $\mu = \min_{i=1,2,3} (\alpha_i + \beta_i + \gamma_i)$.

### 3.6. Practical consideration

Let $X_n = \{\mathbf{x}_i, i = 1, \dots, n\}$ be a set of $n$ scattered data in $\mathbb{R}^2$ and let $\mathcal{T}_m = \{T_\alpha, \alpha = 1, \dots, m\}$ be a triangulation of $X_n$. We denote by $C = \text{conv}(X_n)$ the convex hull of $X_n$. We define the global approximation operator $\Pi_{\mathcal{T}_m}$ by setting, for any $f \in C^{1,1}(C)$ and $\mathbf{x} \in C$,
$$ \Pi_{\mathcal{T}_m}[f](\mathbf{x}) = \Pi[f, T_\alpha](\mathbf{x}), \quad \text{if } \mathbf{x} \in T_\alpha, \alpha = 1, \dots, m, $$

where $\Pi[f, T_\alpha]$ is the approximation operator defined in (58), based on the triangular element $(T_\alpha, G_2(T_\alpha), \tilde{\Sigma}_{T_\alpha})$. In this case, Theorem 3.19 gives the following global $L^1$-error bound

$$ \begin{aligned} \|f - \Pi_{\mathcal{T}_m}\|_{L^1(C)} &= \int_C |f(\mathbf{x}) - \Pi_{\mathcal{T}_m}[f](\mathbf{x})| \, d\mathbf{x} \\ &= \sum_{\alpha=1}^{m} \int_{T_\alpha} |f(\mathbf{x}) - \Pi[f, T_\alpha](\mathbf{x})| \, d\mathbf{x} \\ &\le \frac{L(\nabla f)}{24} \left(1 + \frac{4}{\mu}\right) \sum_{\alpha=1}^{m} |T_\alpha| \omega(T_\alpha). \end{aligned} \quad (107) $$

For any triangulation $\mathcal{T}_m$ we denote by

$$ E_1(\mathcal{T}_m) := \sum_{\alpha=1}^{m} |T| \omega(T_\alpha). $$

The global error bound (107) is proportional to $E_1$, by a constant factor which is independent on the triangulation $\mathcal{T}_m$. In analogy, by using the results of Corollaries 3.12 and 3.22 we denote by

$$ E_2(\mathcal{T}_m) := \sum_{\alpha=1}^{m} R_{\alpha}^2, $$

$$ E_3(\mathcal{T}_m) := \sum_{\alpha=1}^{m} R_{\alpha}^4, $$

where $R_\alpha$ is the circumradius of the triangle $T_\alpha \in \mathcal{T}_m$. By using the optimality results of Delaunay triangulation, which can be found in [15] and [17], it is possible to prove the following result.

**Theorem 3.24.** Let $X_n$ be a set of n scattered data in $\mathbb{R}^2$ and let $\mathcal{T}_m$ be a triangulation of $X_n$. Then $E_i(\mathcal{T}_m)$, $i = 1, 2, 3$ achieve their minimum if and only if $\mathcal{T}_m$ is the Delaunay triangulation.

In this paper we have introduced a new class of nonconforming finite elements by enriching the standard triangular linear element with additional functions which are not necessarily polynomials. We provide a simple condition on the enrichment functions, which is both necessary and sufficient, that guarantees the existence of a family of such enriched elements. We have also derived an error bound in both $L^\infty$ and $L^1$ norms. Future researches can be carried out in order to extend the proposed here technique to triangular element with one curved side (see for example [7,9,23,26]).

## Acknowledgements

This research has been achieved as part of RITA "Research ITalian network on Approximation" and as part of the UMI group "Teoria dell'Approssimazione e Applicazioni". The research was supported by GNCS-INdAM 2022 projects and by the grant "Bando Professori visitatori 2022" which has allowed the visit of Prof. Allal Guessab to the Department of Mathematics and Computer Science of the University of Calabria in the spring 2022. The first, the second and the last author are members of the INdAM Research group GNCS.

## References

[1] M. Abramowitz, I.A. Stegun, Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables, vol. 55, US Government Printing Office, 1948.

[2] B. Achchab, K. Bouihat, A. Guessab, G. Schmeisser, A general approach to the construction of nonconforming finite elements on convex polytopes, Appl. Math. Comput. 268 (1) (2015) 916–923.

[3] M. Bachar, A. Guessab, A simple necessary and sufficient condition for the enrichment of the Crouzeix-Raviart element, Appl. Anal. Discrete Math. 10 (2) (2016) 378–393.

[4] M. Bachar, A. Guessab, Characterization of the existence of an enriched linear finite element approximation using biorthogonal systems, Results Math. 70 (3–4) (2016) 401–413.

[5] B.C. Carlson, Special Functions of Applied Mathematics, Academic Press, New York, NY (USA), 1977.

[6] T. Cătinaş, The bivariate Shepard operator of Bernoulli type, Calcolo 44 (4) (2007) 189–202.

[7] T. Cătinaş, Nielson interpolation operators on an arbitrary triangle with one curved side, BIT Numer. Math. 61 (3) (2021) 757–770.

[8] P.G. Ciarlet, The Finite Element Method for Elliptic Problems, SIAM, 2002.

[9] G. Coman, T. Cătinaş, Interpolation operators on a triangle with one curved side, BIT Numer. Math. 50 (2) (2010) 243–267.

[10] M. Crouzeix, P.-A. Raviart, Conforming and nonconforming finite element methods for solving the stationary Stokes equations I, Rev. Fr. Autom. Inform. Rech. Opér. Math. 7 (R3) (1973) 33–75.

[11] P.J. Davis, Interpolation and Approximation, Courier Corporation, 1975.
[12] F. Dell'Accio, F. Di Tommaso, Bivariate Shepard-Bernoulli operators, Math. Comput. Simul. 141 (2017) 65-82.

[13] F. Dell'Accio, F. Di Tommaso, O. Nouisser, B. Zerroudi, Increasing the approximation order of the triangular Shepard method, Appl. Numer. Math. 126 (2018) 78-91.

[14] F. Dell'Accio, F. Di Tommaso, A. Guessab, F. Nudo, On the improvement of the triangular Shepard method by non conformal polynomial elements, Appl. Numer. Math. 184 (2023) 446-460.

[15] N.P. Dolbilin, H. Edelsbrunner, O.R. Musin, On the optimality of functionals over triangulations of Delaunay sets, Russ. Math. Surv. 67 (4) (2012) 781-788.

[16] R. Franke, G. Nielson, Smooth interpolation of large sets of scattered data, Int. J. Numer. Methods Eng. 15 (11) (1980) 1691-1704.

[17] A. Guessab, Approximations of differentiable convex functions on arbitrary convex polytopes, Appl. Math. Comput. 240 (2014) 326-338.

[18] A. Guessab, Sharp Approximations Base on Delaunay Triangulations and Voronoi Diagrams, NSU Publishing and Printing Center, Novosibirsk, 2022.

[19] A. Guessab, Q. Rahman, Quadrature formulae and polynomial inequalities, J. Approx. Theory 90 (2) (1997) 255-282.

[20] A. Guessab, G. Schmeisser, Convexity results and sharp error estimates in approximate multivariate integration, Math. Comput. 73 (247) (2004) 1365-1384.

[21] A. Guessab, Y. Zaim, A Unified and General Framework for Enriching Finite Element Approximations, Springer, 2017.

[22] R.B. Manfrino, J.A.G. Ortega, R.V. Delgado, Inequalities: A Mathematical Olympiad Approach, Birkhäuser Verlag AG, 2009.

[23] R. McLeod, A. Mitchell, The construction of basis functions for curved elements in the finite element method, IMA J. Appl. Math. 10 (3) (1972) 382-393.

[24] A.S. Posamentier, I. Lehmann, The Secrets of Triangles: A Mathematical Journey, Prometheus Books, 2012.

[25] A. Ralston, P. Rabinowitz, A First Course in Numerical Analysis, Courier Corporation, 2001.

[26] M. Zlámal, Curved elements in the finite element method. I, SIAM J. Numer. Anal. 10 (1) (1973) 229-240.

[27] C. Zuppa, Error estimates for modified local Shepard's interpolation formula, Appl. Numer. Math. 49 (2) (2004) 245-259.