Quasi-Interpolation in a Space of C 2 Sextic Splines over Powell–Sabin Triangulations

: In this work, we study quasi-interpolation in a space of sextic splines deﬁned over Powell– Sabin triangulations. These spline functions are of class C 2 on the whole domain but fourth-order regularity is required at vertices and C 3 regularity is imposed across the edges of the reﬁned triangulation and also at the interior point chosen to deﬁne the reﬁnement. An algorithm is proposed to deﬁne the Powell–Sabin triangles with a small area and diameter needed to construct a normalized basis. Quasi-interpolation operators which reproduce sextic polynomials are constructed after deriving Marsden’s identity from a more explicit version of the control polynomials introduced some years ago in the literature. Finally, some tests show the good performance of these operators. By taking into account the multi-afﬁne property of the polar form, the claim followed.


Introduction
Spline functions over triangulations have been, and are, the object of intense research for their role in the Approximation Theory and for dealing with a wide variety of problems of practical interest, among which the approximation of scattered data and the numerical solution of partial differential equations occupy a prominent place.
It is well known that the requirement of a regularity C m of a spline on a given triangulation implies that the degree must be greater than or equal to 4m + 1 [1]. As in practice, it is essential to use splines of the lowest degree for a given class, different finite elements obtained by subdiving every triangle have been introduced and analyzed in the literature, among them the Clough-Tocher (CT-) and Powell-Sabin (PS-) refinements (see [2] and [3], respectively).
The PS-refinement was proposed in [3] for contour plotting. A first subdivision into six triangles is achieved by selecting an inner point in every triangle and connecting it with similar points in the adjacent triangles as well as with the three vertices. The inner point of a boundary triangle is joined to a point over a boundary edge when no adjacent triangle is available. From this PS6-split, a PS12-split is easily derived by joining in every triangle the three points lying on the edges of the triangle that the previous construction produces [4]. Furthermore, in this case, PS-splines of higher degree and smoothness have been constructed [5].
The application of splines in various fields requires efficient algorithms for constructing locally supported bases for the spline spaces. The B-spline representation of bivariate C 1 quadratic splines achieved by Dierckx [19] was essential in the development of spline spaces on PS partitions and applications. The method proposed by P. Dierckx is completely geometrical, it is reduced to finding a set of PS2 triangles that must contain a number of specified points. Linear and quadratic programming problems are the standard methods proposed by many authors in the literature [19,20]. The main idea of both methods is to minimize the area of a triangle without imposing any condition concerning the diameter of the sought triangles. Moreover, the quadratic problem only provides local maxima. In order to avoid these limitations, we present an algorithm that aims to produce PS6 triangles with a smaller area and diameter, and compare it with the one proposed in [21].
The study of spline function spaces on Powell-Sabin partitions obtained by a refinement into six sub-triangles, has attracted a great interest in the scientific community since its introduction. The cubic case was considered in [20,22]. Spaces of quintic splines were analyzed in [23] and more recently in [24], among others. In [25] and [26], normalized bases for PS-splines of degree 3r − 1 are defined and super-splines of arbitrary degree are given, respectively.
Quasi-interpolation over Powell-Sabin triangulations for specific spaces have been also studied in depth [11][12][13][14][15][16][17]24], as well as for a family of spaces [27]. The construction of such operators is based on establishing a Marsden's identity. It is a powerful tool that allows to write the monomials in terms of the corresponding B-spline-like functions (Bsplines for short). In this view, we establish a general Marsden's identity in the subspace of sextic splines from an easy approach based on a version of the control polynomials different from the one used in [25].
In this paper, we revise a subspace of C 2 sextic PS6 splines obtained by imposing additional smoothness requirements at the interior points of the triangulation chosen to construct the sub-triangulation and also across some edges of the refined triangulation. This subspace of splines was studied in [26], where it is shown that every spline is uniquely determined by its values at the vertices of the initial triangulation and the interior points and those of its partial derivatives up to the fourth order at the vertices.
The paper is organized as follows: In Section 2, we recall some general concepts of polynomials on triangles and recall the notion of control polynomials. In Section 3, we recall some results concerning the space of sextic PS6 splines as well as the Hermite interpolation problems needed to obtain the B-splines forming a normalized basis, and provide explicit expressions for the Bernstein-Bézier coefficients of the B-splines in order to provide a fully worked out construction. Furthermore, the Marsden's identity is stated. In Section 4, an algorithm for determining the set of PS triangles needed to define the B-splines is proposed, which aims to obtain PS splines with a small area and diameter. In Section 5, control polynomials are used to define quasi-interpolation operators having an optimal approximation order. Furthermore, some tests are carried out to illustrate the performance of such operators. A section of conclusions is also included.

Bernstein-Bézier, Polar Forms and Control Polynomials
Firstly, we introduce some notation, results relative to polynomials defined on triangles and the notion of polar forms.
Consider a non-degenerated triangle T := V 1 , V 2 , V 3 with vertices V i := (x i , y i ), i = 1, 2, 3. It is well-known that every point V := (x, y) ∈ R 2 can be uniquely expressed as where the barycentric coordinates τ := (τ 1 , τ 2 , τ 3 ) with respect to T are the unique solution of the system   x 1 x 2 x 3 y 1 y 2 y 3 1 1 1 Hereafter, we denote by P d the linear space of bivariate polynomials of degree less than or equal to d. Any bivariate polynomial p ∈ P d has a unique representation in barycentric coordinates 3 3 are the Bernstein-Bézier polynomials of degree d with respect to T. The coefficients b β are called the Bézier ordinates of the polynomial p with respect to the triangle T, and b(τ) is stated to be the Bernstein-Bézier form of p. It may be represented by associating each coefficient b β with the domain points ξ β determined by the barycentric coordinates with respect to triangle T (see Figure 1). The points ξ β , b β ∈ R 3 are the control points of the so called B-net for the surface of equation z = p(x, y). This surface is tangent at the vertices of T to the linear piecewise function defined by the B-net. The graph of the surface is contained in the convex hull of the control points and p can be easily bounded from them. The conversion of the Bézier form to a different triangleT can be neatly expressed in terms of the polar form [28,29]. Recall that the polar form B[p d ] of a polynomial p d ∈ P d is completely characterized by three properties: it is symmetric, multi-affine and diagonal, i.e., For further use, let us recall the following restricted version of Lemma 4.1 given in [11]. Lemma 1. Let d 1 and d 2 be two positive integers, with d 2 ≤ d 1 . Then, for any polynomial p ∈ P d 1 and any points V 1 , . . . , is a polynomial of degree ≤ d 2 . Moreover, for any points W 1 , . . . , W d 2 in R 2 it holds The behavior of the controlled spline function at any vertex can be detected from the behavior of control polynomials at the same vertex [20]. Now, by using the relationship between polynomials and their blossoms, we obtained a result that will allow to define the control polynomial that was the main tool for establishing Marsden's identity which is the key for building quasi-interpolation schemes based on a C 2 sextic PS-spline space. The notation ∂ a,b f (A) was used for the partial derivative ∂ a+b f ∂x a ∂y b (A) with a + b ≥ 0. The following result, that defines the control polynomial of degree d 2 at the vertex V 1 of a polynomial p of degree d 1 , is an alternative way to establish Marsden's identity. Proposition 1. Let d 1 and d 2 be two positive integers, with d 2 ≤ d 1 . Let p ∈ P d 1 and V 1 ∈ R 2 . For a given real number 0 < θ < 1, we define the polynomial q of degree d 2 by A general result is proved in Theorem 1 in [27], where more detailed information is given.
Hereafter, D r (V 1 ) denotes the disk of radius r around the vertex V 1 of a triangle T = V 1 , V 2 , V 3 . It is the subset of domain points ξ β defined as In the following section, we present a completely worked out construction of a normalized basis of the subspace of splines over a Powell-Sabin partition introduced in [26].

Explicit Construction of a B-Spline Basis for a Space of Powell-Sabin Super Splines
Let Ω be a polygonal domain in R 2 , and let ∆ : be a regular triangulation of Ω. We denoted by V i := (x i , y i ) T , i = 1, . . . , nv, the vertices of the given triangulation. Let ∆ PS be a PS-refinement of ∆, which divides each macro triangle T j into six micro-triangles (see Figure 2). This partition was defined algorithmically as follows:

1.
Choose an interior point Z j in each triangle T j . If two triangles T i and T j have a common edge, then the line joining Z i and Z j should intersect the common edge at some point R ij .

2.
Join each point Z j to the vertices of T j .

3.
For each edge of the triangle T j : (a) which is common to a triangle T i , join Z j to R ij ; (b) which belongs to the boundary ∂Ω, join Z j to an arbitrary point on that edge.
The space of sextic piecewise polynomials on ∆ PS with global C 2 continuity was defined as: This kind of spline space can be considered as the Hermitian finite element [30]. Then, we considered a particular subspace of S 2 6 (Ω, ∆ PS ) introduced in [26].
and E * be, respectively, the subsets of vertices in ∆, split points in ∆ PS , and edges in ∆ PS that connect a split point Z i to a point R ij . As given in [26], the space of PS splines is defined as Each C 2 (Ω) function s is of class C 4 at any vertex in V and of class C 3 at any split point in Z and across any edge in E * . In [26], by using minimal determining sets, it was proved that for given values f a,b i , i = 1, . . . , nv, and g k , k = 1, . . . , nt, there exists a unique spline s ∈ S 2,4,3 6 (Ω, ∆ PS ) such that Therefore, the dimension of the space S 2,4,3 6 (Ω, ∆ PS ) was equal to 15nv + nt. A procedure for the construction of a normalized basis for the space S 2,4,3 6 (Ω, ∆ PS ) was then based on the solution of the above Hermite interpolation problem for appropriate values f a,b i and g k (see [26]). Non-negative and locally supported basis functions B v i,j and B t k with respect to vertices and triangles, respectively, that form a partition of unity were defined, and any s ∈ S 2,4,3 6 (Ω, ∆ PS ) could be represented as In what follows, we gave a fully elaborate construction of such a normalized basis [25,26]. For every vertex V i , let M i := ∪ T∈∆,V i ∈ T T be the molecule of vertex V i , i.e., the union of all triangles in ∆ containing V i . For all vertices V , ∈ Λ i (e.g., Λ i is the set of indices for the vertices that form an edge in ∆ with V i ), lying on the boundary of M i , let The points V i and S i, , ∈ Λ i , were stated to be PS6 points associated with the vertex V i . 3 ) be a triangle containing the PS6 points of V i . It will be called the PS6 triangle. Denote by B 4 t i ,mnl , m + n + l = 4, the Bernstein polynomials of degree four with respect to t i , and define, for all 0 ≤ a + b ≤ 4, the values They are used to define the B-splines B v i,j and B t k as follows.

Vertex B-Spline
Every B-spline B v i,j , 1 ≤ j ≤ 15, relative to the vertex V i was defined as the unique solution of a particular Hermite interpolation with conditions given by (3). Firstly, all f a,b were equal to zero except for = i, and f a,b i = α a,b i,j . Moreover, if V i was a vertex of a triangle T k := V 1 , V 2 , V 3 , then g k was equal to a value β k i,j to be precise later and the remaining values of g were all equal to zero. The spline defined in this way was zero outside the molecule M i of vertex V i . Next, we computed the BB-coefficients of B v i,j relative to the triangles determining their support. For the sake of simplicity, we only computed the BB-coefficients of the B-spline B v 1,j relative to the vertex V 1 of triangle T k . The corresponding Bézier ordinates are schematically represented in Figure 3.  From the definition of B v 1,j , many BB-coefficients were equal to zero. Figure 4 shows the refinement of T k and we assumed that the points indicated in the figure had the following barycentric coordinates: Because of the C 4 smoothness of the spline at V 1 , the ordinates c 1 , c 2 , . . . , c 25 were uniquely determined by the values α a,b 1,j , 0 ≤ a + b ≤ 4. The ordinates c 26 , . . . , c 34 were obtained by the C 3 smoothness across the edge R 12 , Z .  , respectively. Before subdivision, their BB-coefficients were respectively, wherê Therefore, we obtained The values c 35 , . . . , c 43 were determined using a similar method. They were given by the following expressions: The remaining Bézier ordinates had to be chosen in such a way that the B-spline was C 3 continuous at split point Z. Therefore, let us first define the points and let p 3 ∈ P 3 be the polynomial of degree 3 defined over the triangle T = w 1 , w 2 , w 3 with ordinates We obtained c 44 = λ 2 12 c 13 + 2λ 12 λ 21ĉ20 , c 45 = 12λ 3 c 13 + 3λ 2 12 λ 21ĉ20 , c 46 = λ 2 12ĉ 20 , The choice β k 1,j = c 68 provided the values needed to completely define the B-spline B v 1,j . Figure 5 shows typical plots of the fifteen C 2 sextic B-splines associated with a vertex of the triangulation.

Triangle B-Spline
For the sake of simplicity, we denoted by b l the B-ordinates with respect to a triangle (see Figure 6). The B-spline B t k with respect to triangle T k was defined as the spline satisfying conditions (3) with all f a,b i equal to zero, g k = β k and the remaining values of g equal to zero. It vanished outside T k . In order to specify the value of β k , we looked at the Bernstein-Bézier representation of the B-spline B t k . We consider, again, the macro-triangle T k = V 1 , V 2 , V 3 , as above.
Let us define again a polynomial p 3 ∈ P 3 of degree 3 defined on the triangle T = w 1 , w 2 , w 3 , where w i were defined in (6), and having the following B-ordinates: Furthermore, as in the above subsection, we obtained From the construction, it was clear that all the Bézier ordinates were non-negative. Then, the B-spline B t k was non-negative. We could choose β k = 6z 1 z 2 z 3 . For each vertex V i and each triangle T k , we defined points Q i,β := X i,β , Y i,β with β := (β 1 , β 2 , β 3 ), |β| := β 1 + β 2 + β 3 = 4 and Q t k := X t k , Y t k in such a way that the reproduction of the monomials x and y held, i.e., then (9) and (10) hold.
Proof. For all (x, y) ∈ t i , we had Using (5) and (11), we obtained (9). Now, to prove (10), we needed to show that Recall that, in the construction of B-splines in the above section, the value of a PS6 spline at a split point Z was computed through a particular cubic polynomial evaluated at the split point. We considered again the macro-triangle T k = V 1 , V 2 , V 3 . The two cubic polynomials corresponding to the two PS6 splines in the equations (9) and (10) were denoted by p x,3 (τ) and p y,3 (τ). They were defined on the triangle with the vertices given in (6). The Bézier ordinates of p x,3 were given by the following expressions: By the definition of Q t k , it held Therefore, it was clear that p x,3 (τ) = τ 1 b x 300 + τ 2 b x 030 + τ 3 b x 003 , and (12) followed. Hence, (9) was proved. Equality (10) could be proved in a similar way.

Nearly Optimal PS6 Triangles
The construction of a normalized PS6 basis of S 2,4,3 6 (Ω, ∆ PS ) was reduced to finding a set of PS6 triangles that had to contain a number of specified points. The set of PS6 triangles was not uniquely defined for a given refinement [31]. One possibility for their construction was to calculate triangles of minimal area, the so-called optimal PS triangles introduced by P. Dierckx [19]. Computationally, this problem led to a quadratic programming problem. From a practical point of view, other choices may have been more appropriate. An alternative (and easier to implement) solution is given in [31], where the sides of the PS triangle are obtained by connecting neighboring PS-points in a suitable way. This technique was adopted and improved in [21]. A particular choice of the PS6 triangles could also simplify the treatment of boundary conditions. For quasi-interpolation (see [17]), the corners of each PS6 triangle were preferred to be chosen on edges of the triangulation.
We recalled the standard method proposed in the literature [19,20,25] to construct PS6 triangles and, then, we introduced a novel procedure.

Quadratic Programming Problem
Consider points Q i,j = (X i,j , Y i,j ), j = 1, 2, 3, yielding a PS6 triangle relative to the vertex V i = (x i , y i ) and triplets Γ i,j , Γ x i,j , Γ y i,j , j = 1, 2, 3, satisfying the following equality: The area of the PS6 triangle being then, maximizing the objective function Γ x i,1 Γ y i,2 − Γ y i,1 Γ x i,2 is one approach to obtain a triangle of smallest area. Additional constraints are needed to obtain a PS6 triangle containing all PS6 points with respect to V i . The classical construction due to Dierckx was then summarized in the next result.

Proposition 3.
The construction of an optimal PS6 triangle t i with respect to vertex V i is equivalent to the following quadratic programming problem: find a set of triplets The objective function of the optimization problem can be written as max 1 2 x T A x, where The eigenvalues of the matrix A were −1, −1, 1, 1, 0 and 0, so that A was indefinite. As pointed out in [19], "since the Hessian matrix of the objective function is not negative (semi-) definite, appropriate software can only find a local maximum". Therefore, we could not guarantee that the quadratic optimization problem has a unique solution, which led to a scenario of local solutions.
The technique for determining PS6 triangles is not unique. One option for their construction is to calculate a triangle with a minimal area. Although the quadratic program of P. Dierckx [19] produces excellent results, it can also produce PS6 triangle with quite large diameters. Therefore, in order to overcome the limitation of the above optimization problem, namely, the appearance of pre-degenerated triangles, i.e., triangles with a minimal area and long diameters which impact negatively the quality of the approximation, we proposed an algorithm yielding a PS6 triangle with a diameter as small as possible.

Algorithm for Determining a Triangle Containing a Set of Points
i=0 be the interiors of the seven regions obtained by extending the edges of T indefinitely (see Figure 8). Then, for each fixed 0 ≤ i ≤ 6, the barycentric coordinates of all the points in Ω i have constant signs. In particular, a point lies in the interior of T if and only if its barycentric coordinates are positive.
The algorithm proposed here to define a triangle containing the points A i , 1 ≤ i ≤ n, started from an initial triangle and built step by step triangles so that triangle T j := A j 1 , A j 2 , A j 3 , j ≥ 2, obtained at the j th step of the algorithm contained the points A 1 ,. . . , A j−1 . Denote by Ω j k , k = 0, 1, 2, 3, and Ω j k,k+1 , k = 1, 2, 3, the seven regions obtained by dividing the plane through T j (see Figure 8). More precisely, the procedure described in Algorithm 1 was carried out to determine a triangle from the previous one. Figure 13 shows the PS6 triangles produced by the algorithm when applied to the PS6 points close to those used in [21]. They had two or three edges in common with the convex hull of the PS6 points.
Next, we gave a result needed to determine triangles having a nearly minimal area.

Algorithm 1 DETERMINING THE TRIANGLE T j+1 FROM T j
Require: compute the barycentric coordinates of A j with respect to T j and select the region where A j is located.
if A j ∈ Ω j 0 then A j is in T j , perform T j+1 ← T j and move to the next point A j+1 3. Define T j+1 as the triangle of minimum area among T 1 j+1 and T 2 j+1 . The same process is used if A j belongs to Ω Lemma 2. Let a, A 1 , A 2 , A 3 and A 4 be five points in R 2 . If a ∈ T ijk := A i , A j , A k for i, j, k = 1, 2, 3, 4 and i = j = k, then, a is in the triangle obtained by applying the algorithm using T ijk and A l , l = i = j = k.
Proof. For the sake of simplicity, consider only one of the four different triangles which can be obtained from four points. Let T 134 := A 1 , A 3 , A 4 be a triangle containing a. By applying the algorithm proposed here to T 134 and A 2 , we can distinguish the following scenarios: • If A 2 ∈ T 134 , then, the resulting triangle will be T 134 itself. • If A 2 / ∈ T 134 , then the obtained triangle will contain T 134 .
In both cases the resulting triangle will contain T 134 , so will contain also a. The proof is complete.
From Lemma 2, at step j in the algorithm, we used the four triangles obtained by a permutation of the vertices of T j and A j , and we chose the triangle of the small diameter among the four ones. Figure 9 shows the PS6 triangles provided by the proposed algorithm for the considered triangulation. It can be noticed that the resulting triangles passed through at least three PS6 points. They had near minimal areas and smaller diameters.
As stated before, the quadratic optimization problem proposed by P. Dierckx [19] can produce PS6 triangles with quite large diameters, and the algorithm proposed here aimed to avoid this problem even though the resulting triangles had no minimal areas. Figure 10 shows the results provided by the Dierckx's method and the algorithm for minimizing the diameter when a near degenerate vertex was considered. Figure 11 shows the results obtained when the time of execution of both algorithms was examined. The time required by Dierckx' algorithm was more than 30 times longer than that required by the proposed algorithm.  Other algorithms for determing PS triangles have also been described in the literature. As stated before, in [21], after Proposition 1, the authors outline an algorithm that produces PS triangles sharing two or three edges with the convex hull of the PS points. Next, we compare it with our Algorithm. To conduct that, we considered PS points such as those in Figure 1 in [21]. They are represented in Figure 12. Algorithm 1 provided the PS6 triangles shown in Figure 13. Each of them was produced from a choice of an initial triangle. On the left side, we showed those obtained after three steps starting from the small dark triangle. We saw that these PS triangles shared two or three sides with the convex hull of the PS points. On the right side, we showed two other PS triangles produced by the algorithm after four steps. They also shared two or three sides with the convex hull. The results provided by the algorithm in [21] and Algorithm 1 were similar, although the later one did not need to compute the convex hull of the PS points.

Quasi-Interpolation Schemes with Optimal Approximation Order
In this section, we gave proof of the Marsden's identity for the space S 2,4,3 6 (Ω, ∆ PS ), expressing any spline s in this space as a linear combination of the normalized sextic Powell-Sabin B-splines defined above. The coefficients in that combination were given in terms of the polar forms of s. Therefore, Proposition 1 facilitates the establishment of Marsden's identity in comparison with other existing methods (e.g., matrix inverse [20]).
Here, we used the same notation as in Section 3. Let Q i,j , j = 1, 2, 3, be the vertices of a PS6 triangle t i with respect to V i . Definẽ We had the following result.
Corollary 1. For any p ∈ P 6 it holds where V k1 , V k2 and V k3 are the vertices of the macro triangle containing Z k .

Proof. Define
We proved that from which the equality s = p follows. It was clear that Then, we used the notion of control polynomial developed in Section 2. Let be the control polynomial of degree 4 of p at the vertex V i . We could writeq vi on the PS-triangle t i as According to Lemma 1, Using Proposition 1, we deduced that Then, it sufficed to prove that s(Z k ) = p(Z k ). Without the loss of generality, we proved the equality only for one triangle in ∆. Let T = V 1 , V 2 , V 3 be a triangle in ∆ with split point Z 1 . Then, Similarly, By taking into account the multi-affine property of the polar form, the claim followed.
Then, we stated the following result, whose proof followed the idea used in [12] in dealing with quadratic Powell-Sabin splines. (Ω, ∆ PS ), it holds where s i := s |t i stands for the restriction of s to the triangle t i in ∆ PS ands k is the restriction of s to a triangle t k = V k1 , V k2 , V k3 containing Z k .
Proof. Consider a spline s in S 2,4,3 6 (Ω, ∆ PS ). Let t i be a triangle in ∆ PS having V i as a vertex. Let s i be the restriction of s to t i , i.e., the sextic polynomial such that Let p i be the restriction of s on t i . From Corollary 1, it was clear that for all (x, y) ∈ Ω and r = 1, . . . , n v it held Then, Therefore, Define, Then, for all r = 1, . . . , n v , we obtained Similarly, we obtained Since every element in S 2,4,3 6 (Ω, ∆ PS ) was uniquely determined by its values and derivative values up to order four at the vertices of ∆, then the claim followed and the proof was completed.
Marsden's identity is a useful tool for constructing quasi-interpolants to enough regular functions (see [12] and references therein for details). We used it to define differential quasi-interpolants in S 2,4,3 6 (Ω, ∆ PS ). Only an outline of the construction was given here.
Lef f ∈ C 6 (Ω) and L j i := L j i,x , L j i,y , i = 1, . . . , nv, j = 1, . . . , 15, be some fixed points lying in the union of all triangles in ∆ having V i as vertex. Let us suppose that they formed an unisolvent scheme in P 6 R 2 , and let p j i be the Taylor polynomial of f of degree six at L j i , i.e., p Let p k be the Taylor polynomial of degree 6 at point L k in the support of B t k . Define Let Q f be the quasi-interpolant defined by (16) and (15). Then, the quasi-interpolation operator Q : C 6 (Ω) → S 2,4,3 6 (Ω, ∆ PS ) defined such that Q( f ) := Q f was exact on P 6 , i.e., Q(p) = p for all p ∈ P 6 .
Moreover, if each L j i belonged to a triangle τ j i in ∆ PS with V i as a vertex, then Q(s) = s for any spline s ∈ S 2,4,3 6 (Ω, ∆ PS ). Which meant that the quasi-interpolation operator Q was the projector.
The quasi-interpolation error was estimated as max ,k=1,..., 50 where x i and y j were equally spaced points in (0, 1). The numerical convergence order (NCO) was given by the rate where E(m) stands for the estimated error associated with ∆ m . The estimated errors and NCOs for the functions f 1 and f 2 are shown in Table 1. They confirmed the theoretical results. Figure 15 shows two of the meshes used to define quasi-interpolants for the test functions f 1 and f 2 . Figure 16 shows the plots of the splines Q f 1 and Q f 2 for the two above meshes.

Conclusions
In this paper, a fully carried out construction of a normalized basis of the space S 2,4,3 6 (Ω, ∆ PS ) introduced in [26] was given and an algorithm was proposed and compared with two others in the literature. Furthermore, an efficient manner to establish Marsden's identity was detailed from which quasi-interpolation operators with optimal approximation order were defined. Some tests showed the good performance of these operators.