Paolo Piersanti
email: ppiersan-c@my.cityu.edu.hk

Xiaoqin Shen
email: xqshen@xaut.edu.cn

Numerical methods for static shallow shells lying over an obstacle

Keywords: Shallow shell, Enriching operator, Non-conforming finite element method, Obstacle problems, Elliptic variational inequalities

In this paper a finite element analysis to approximate the solution of an obstacle problem for a static shallow shell confined in a half space is presented. To begin with, we establish, by relying on the properties of enriching operators, an estimate for the approximate bilinear form associated with the problem under consideration. Then, we conduct an error analysis and we prove the convergence of the proposed numerical scheme.

Introduction

The study of unilateral contact problems in elasticity arises in many applicative fields such as structural mechanics and civil engineering. Obstacle problems have lately been studied in, for instance, [START_REF] Ciarlet | An obstacle problem for elliptic membrane shells[END_REF][START_REF] Ciarlet | Un problème de confinement pour une coque membranaire linéairement élastique de type elliptique[END_REF][START_REF] Ciarlet | Obstacle problems for Koiter's shells[END_REF][START_REF] Ciarlet | A confinement problem for a linearly elastic Koiter's shell[END_REF][START_REF] Léger | A linearly elastic shell over an obstacle: The flexural case[END_REF][START_REF] Rodríguez-Arós | Mathematical justification of the obstacle problem for elastic elliptic membrane shells[END_REF].

The numerical analysis of obstacle problems has been arising the interest of many scientists since the late nineties. In this direction, a very direct and mathematically elegant approach is the one making use of enriching operators, the properties of which were studied by S. C. Brenner and her collaborators in the seminal papers [START_REF] Brenner | Two-level additive Schwarz preconditioners for nonconforming finite element methods[END_REF][START_REF] Brenner | A two-level additive Schwarz preconditioner for nonconforming plate elements[END_REF][START_REF] Brenner | Convergence of nonconforming multigrid methods without full elliptic regularity[END_REF][START_REF] Brenner | C 0 interior penalty method for fourth order elliptic boundary value problems on polygonal domains[END_REF]. These general theoretical results were then used to study finite element methods for obstacle problems, which can be found in [START_REF] Brenner | Finite element methods for the displacement obstacle problem of clamped plates[END_REF] and [START_REF] Brenner | A Morley finite element method for the displacement obstacle problem of clamped Kirchhoff plates[END_REF]. Nonconforming finite element methods for obstacle problem were also studied in [START_REF] Carstensen | Nonconforming fem for the obstacle problem[END_REF].

In this paper, we study the displacement of a static shallow shell lying over a planar obstacle from the numerical point of view, using a suitable finite element method. Shallow shells theory is extensively described in the books [START_REF] Ciarlet | Mathematical Elasticity[END_REF] and [START_REF] Ventsel | Thin Plates and Shells: Theory, Analysis, and Applications[END_REF]. According to this theory, the problem under examination is modelled in terms of a fourth-order differential operator (cf., e.g., [START_REF] Ciarlet | Mathematical Elasticity[END_REF]). The theory of finite element methods for fourth-order problems governed by variational inequalities has been investigated in, for instance, [START_REF] Brezzi | Error estimates for the finite element solution of variational inequalities[END_REF][START_REF] Falk | Error estimates for the approximation of a class of variational inequalities[END_REF][START_REF] Glowinski | Numerical Analysis of Variational Inequalities[END_REF][START_REF] Hlaváček | Solution of Variational Inequalities in Mechanics[END_REF].

Our mathematical model of an obstacle problem for a linearly elastic shallow shell in the static case is inspired by that of Léger and Miara (cf. [START_REF] Léger | Mathematical justification of the obstacle problem in the case of a shallow shell[END_REF] and [START_REF] Léger | Erratum to: Mathematical justification of the obstacle problem in the case of a shallow shell[END_REF]). To our knowledge, there is no reference on the study of numerical analysis of obstacle problems for linearly elastic shallow shells.

In this paper, we extend the method proposed in [START_REF] Brenner | Finite element methods for the displacement obstacle problem of clamped plates[END_REF] and [START_REF] Brenner | A Morley finite element method for the displacement obstacle problem of clamped Kirchhoff plates[END_REF] to derive error estimates for the solution to the obstacle problem under for linearly elastic shallow shells under consideration. The fact that the unknown is a vector field is the main difficulty to cope with in order for proving that the residual of the difference between the exact solution and the approximate solution approaches zero as the mesh size approaches zero.

In order to derive the sought convergence, it was necessary to improve and generalise a number of preparatory results for enriching operators (cf., e.g., [START_REF] Brenner | Two-level additive Schwarz preconditioners for nonconforming finite element methods[END_REF][START_REF] Brenner | A two-level additive Schwarz preconditioner for nonconforming plate elements[END_REF][START_REF] Brenner | Convergence of nonconforming multigrid methods without full elliptic regularity[END_REF][START_REF] Brenner | C 0 interior penalty method for fourth order elliptic boundary value problems on polygonal domains[END_REF]), and another number of preparatory results related to the convergence analysis of the scheme (cf., e.g., [START_REF] Brenner | Finite element methods for the displacement obstacle problem of clamped plates[END_REF]).

The paper is divided into four sections. In Section 2, we present some background and notation. In Section 3 we establish some properties of the enriching operator associated with the variational formulation of the problem under consideration and an estimate for Morley's triangle, used to approximate the transverse component of the displacement. In Section 4, following [START_REF] Brenner | Finite element methods for the displacement obstacle problem of clamped plates[END_REF] and [START_REF] Brenner | A Morley finite element method for the displacement obstacle problem of clamped Kirchhoff plates[END_REF], we introduce an intermediary problem and we prove some technical preparatory lemmas. Finally, in Section 5, the error estimate is derived as a result of an application of the previous results.

Background and notation

For an overview about the classical notions of differential geometry used in this paper see, e.g., [START_REF] Ciarlet | Mathematical Elasticity[END_REF] or [START_REF] Ciarlet | An Introduction to Differential Geometry with Applications to Elasticity[END_REF] while, for an overview about the classical notions of functional analysis used in this paper see, e.g., [START_REF] Ciarlet | Linear and Nonlinear Functional Analysis with Applications[END_REF]. Latin indices, except h, take their values in the set {1, 2, 3} while Greek indices, except ν and ε, take their values in the set {1, 2}. The notation δ αβ designates the Kronecker symbol. Given an open subset Ω of R n , where n ≥ 1, we denote the usual Lebesgue and Sobolev spaces by L 2 (Ω), H 1 (Ω), H 1 0 (Ω), H 2 (Ω), or H 2 0 (Ω); the notation D(Ω) designates the space of all functions that are infinitely differentiable over Ω and have compact supports in Ω; the notation P k (Ω) designates the space of all polynomials of degree ≤ k defined over Ω and the notation P k designates the space of all polynomials of degree ≤ k defined over R n . The Euclidean norm of any point x ∈ Ω is denoted by |x|. In what follows, the compact notation • m,p,Ω , where m ≥ 1 is an integer and p ≥ 1, designates the norm of the Sobolev space W m,p (Ω). The special notation

• m,Ω , where m ≥ 1 is an integer, denotes the norm of the space H m (Ω). If m = 0, then

• 0,Ω := • L 2 (Ω)
and, more generically,

• 0,p,Ω := • L p (Ω) for all p ≥ 1.

The special notation | • | m,Ω , where m ≥ 1 is an integer, denotes the standard semi-norm of the space H m (Ω).

Let ω ⊂ R 2 be a convex polygonal domain, namely a non-empty bounded open connected subset of R 2 with Lipschitz continuous boundary γ := ∂ω and such that ω is all on the same side of γ. Let y = (y α ) denote a generic point in ω and let ∂ α := ∂/∂y α and ∂ αβ := ∂ 2 /∂y α ∂y β .

Referring to [START_REF] Ciarlet | Justification of the two-dimensional equations of a linearly elastic shallow shell[END_REF] (see also Section 3.1 of [START_REF] Ciarlet | Mathematical Elasticity[END_REF] and see also [START_REF] Ventsel | Thin Plates and Shells: Theory, Analysis, and Applications[END_REF]), we recall the rigorous definition of a linearly elastic shallow shell (from now on shallow shell ). We assume that for each ε > 0, we are given a function θ ε ∈ C 3 (ω). We can then define the middle surface of the corresponding shallow shell having thinness equal to 2ε as follows:

ωε := {(y, θ ε (y)); y ∈ ω}.

A rigorous criterion for defining a shallow shell is provided by the existence of a function θ ∈ C 3 (ω), independent of ε, such that θ ε (y) = εθ(y) for all y ∈ ω.

(

This means that, up to an additive constant, the mapping θ ε : ω → R, measuring the deviation of the middle surface of the reference configuration of the shell from a plane, should be of the same order as the thinness of the shell. The shallow shells here considered are made of a homogeneous and isotropic material, they are clamped on their lateral boundary and they are subjected to both applied body forces and applied surface forces. The elastic behaviour of the shallow shell is then described by means of its two Lamé constants λ ≥ 0 and µ > 0 (cf., e.g., [START_REF] Ciarlet | Mathematical Elasticity. Vol. I: Three-Dimensional Elasticity[END_REF]).

In what follows ν denotes the outer unit normal vector field to the boundary γ and ∂ ν denotes the outer unit normal derivative operator along γ.

The function space over which the problem is posed is the following:

V (ω) := {η = (η i ) ∈ H 1 (ω) × H 1 (ω) × H 2 (ω); η i = ∂ ν η 3 = 0 on γ}.
We equip the space V (ω) with the norm • V (ω) defined as follows:

ξ V (ω) := ξ 1 1,ω + ξ 2 1,ω + ξ 3 2,ω for all ξ ∈ V (ω). The corresponding semi-norm | • | V (ω) is defined by |η| V (ω) := |η 1 | 1,ω + |η 2 | 1,ω + |η 3 | 2,ω for all η ∈ V (ω).
The obstacle problem studied in [START_REF] Léger | Mathematical justification of the obstacle problem in the case of a shallow shell[END_REF] and [START_REF] Léger | Erratum to: Mathematical justification of the obstacle problem in the case of a shallow shell[END_REF] is modelled by a set of variational equations and a set of variational inequalities and, besides, its solution is a Kirchhoff-Love field (see, for instance, Section 3.4 of [START_REF] Ciarlet | Mathematical Elasticity[END_REF]). As a result, we can "separate" the transverse component of the displacement vector field from the tangential components of the displacement vector field. We thus define the space associated with the tangential components by

V H (ω) := {η H = (η α ) ∈ H 1 (ω) × H 1 (ω); η α = 0 on γ},
and the space associated with the transverse component by

V 3 (ω) := H 2 0 (ω).
Observe that

V (ω) = V H (ω) × V 3 (ω).
The "physical" obstacle is represented by the plane x 3 = 0 and, in what follows, we assume that θ > 0 in ω. This implies θ ε > 0 in ω, i.e., the middle surface of the considered shallow shell is assumed to be above the obstacle and not in contact with the obstacle.

In what follows we state the scaled two-dimensional limit problem, which slightly differs from the one obtained in [START_REF] Léger | Mathematical justification of the obstacle problem in the case of a shallow shell[END_REF] and [START_REF] Léger | Erratum to: Mathematical justification of the obstacle problem in the case of a shallow shell[END_REF] as a result of a rigorous asymptotic analysis, where only the transverse component of the displacement is subjected to the geometrical constraint associated with the obstacle. More specifically, the transverse component of the displacement field belongs to the following non-empty, closed, and convex set (see [START_REF] Léger | Mathematical justification of the obstacle problem in the case of a shallow shell[END_REF]):

K 3 (ω) := {η 3 ∈ V 3 (ω); θ + η 3 ≥ 0 almost everywhere in ω}. (2)
By virtue of the Rellich-Kondrachov theorem, the compact embedding H 2 (ω) → → C 0 (ω) holds (the symbol " → →" denotes a compact embedding and the space C 0 (ω) is equipped with the sup-norm). Hence, by virtue of the fact that θ ∈ C 3 (ω), the set K 3 (ω) defined in (2) also takes the following form:

K 3 (ω) = {η 3 ∈ V 3 (ω); θ + η 3 ≥ 0 in ω}. (3)
Let Ω := ω × ]-1, 1[ , and let x = (x i ) denote a generic point in the set Ω. With each point x = (x i ) ∈ Ω, we associate the point x ε = (x ε i ) defined by

x ε α := x α = y α and x ε 3 := εx 3 , so that ∂ ε α = ∂ α and ∂ ε 3 = 1 ε ∂ 3 .
We assume that the shallow shell under consideration is subjected to applied body forces whose density per unit volume is defined by means of its covariant components

f ε i ∈ L 2 (ω × (-ε, ε))
and applied surface forces whose density per unit area is defined by means of its covariant components g +,ε i ∈ L 2 (ω × {ε}). Applied surface forces associated with the lower face of the reference configuration of the shallow shell are not to be considered since the obstacle is assumed to be rigid. We also assume that there exist functions f i ∈ L 2 (Ω) and g + i ∈ L 2 (ω × {1}) independent of ε such that the following assumptions on the data hold:

f ε α (x ε ) = ε 2 f α (x) at each x = (x i ) ∈ Ω, f ε 3 (x ε ) = ε 3 f 3 (x) at each x = (x i ) ∈ Ω, g +,ε α (x ε ) = ε 3 g + α (x) at each x = (x i ) ∈ ω × {1}, g +,ε 3 (x ε ) = ε 4 g + 3 (x) at each x = (x i ) ∈ ω × {1}.
We are now ready to state the scaled limit problem P(ω), which slightly differs from the one found in [START_REF] Léger | Mathematical justification of the obstacle problem in the case of a shallow shell[END_REF] and [START_REF] Léger | Erratum to: Mathematical justification of the obstacle problem in the case of a shallow shell[END_REF].

Problem P(ω) To find ζ = (ζ H , ζ 3 ) ∈ V H (ω) × K 3 (ω) satisfying the following variational inequalities - ω m αβ (ζ 3 )∂ αβ (η 3 -ζ 3 ) dy + ω n θ αβ (ζ)∂ α θ∂ β (η 3 -ζ 3 ) dy ≥ ω p 3 (η 3 -ζ 3 ) dy - ω s α ∂ α (η 3 -ζ 3 ) dy for all η 3 ∈ K 3 (ω),
and the following variational equations

ω n θ αβ (ζ)∂ β η α dy = ω p α η α dy for all η H = (η α ) ∈ V H (ω) where                            λ ≥ 0, µ > 0 are the Lamé constants, m αβ (ζ 3 ) := - 4λµ 3(λ + 2µ) ∆ζ 3 δ αβ - 4 3 µ∂ αβ ζ 3 , e θ αβ (ζ) := 1 2 (∂ α ζ β + ∂ β ζ α ) + 1 2 (∂ α θ∂ β ζ 3 + ∂ β θ∂ α ζ 3 ), n θ αβ (ζ) := 4λµ λ + 2µ e θ σσ (ζ)δ αβ + 4µe θ αβ (ζ), p i := 1 -1 f i dx 3 + g + i , s α := 1 -1 x 3 f α dx 3 + g + α .
Likewise, since θ ε ∈ C 3 (ω), we define the non-empty closed convex set K ε 3 (ω) by

K ε 3 (ω) := {η 3 ∈ V 3 (ω); θ ε + η 3 ≥ 0 in ω}. (4)
The next step consists in de-scaling problem P(ω). More specifically, the solution (ζ H , ζ 3 ) is de-scaled as follows (cf. [START_REF] Ciarlet | Mathematical Elasticity[END_REF])

ζ ε H = ε 2 ζ H in ω, ζ ε 3 = εζ 3 in ω.
Thanks to (1), if

ζ 3 ∈ K 3 (ω), then ζ ε 3 ∈ K ε 3 (ω)
. The de-scaled problem P ε (ω) can be thus stated and constitutes the point of departure of our numerical analysis.

Problem P ε (ω) To find ζ ε = (ζ ε H , ζ ε 3 ) ∈ V H (ω) × K ε 3 (ω) satisfying the following variational inequalities: - ω m ε αβ (ζ ε 3 )∂ αβ (η 3 -ζ ε 3 ) dy + ω n θ,ε αβ (ζ ε )(∂ α θ ε )∂ β (η 3 -ζ ε 3 ) dy ≥ ω p ε 3 (η 3 -ζ ε 3 ) dy - ω s ε α ∂ α (η 3 -ζ ε 3 ) dy for all η 3 ∈ K ε 3 (ω), (5)
and the following variational equations:

ω n θ,ε αβ (ζ ε )∂ β η α dy = ω p ε α η α dy for all η H = (η α ) ∈ V H (ω), (6)
where

                           λ ≥ 0, µ > 0 are the Lamé constants, m ε αβ (ζ ε 3 ) := -ε 3 4λµ 3(λ + 2µ) ∆ζ ε 3 δ αβ + 4 3 µ∂ αβ ζ ε 3 , e θ,ε αβ (ζ ε ) := 1 2 (∂ α ζ ε β + ∂ β ζ ε α ) + 1 2 (∂ α θ ε ∂ β ζ ε 3 + ∂ β θ ε ∂ α ζ ε 3 ), n θ,ε αβ (ζ ε ) := ε 4λµ λ + 2µ e θ,ε σσ (ζ ε )δ αβ + 4µe θ,ε αβ (ζ ε ) , p ε i := ε -ε f ε i dx ε 3 + g +,ε i , s ε α := ε -ε x ε 3 f ε α dx ε 3 + εg +,ε α . (7)
Clearly, ( 5) and ( 6) can can be combined into a single system of variational equations, whose left hand side is associated with the symmetric bilinear form b(•, •) given by (cf. Sections 3.5, 3.6 and 3.7 of [START_REF] Ciarlet | Mathematical Elasticity[END_REF])

b(ζ ε , η) = - ω m ε αβ (ζ ε 3 )∂ αβ η 3 dy + ω n θ,ε αβ (ζ ε )(∂ α θ ε )∂ β η 3 dy + ω n θ,ε αβ (ζ ε )∂ β η α dy. (8)
A straightforward computation shows that b(η, η) :=

ω 4λµ λ + 2µ ε 3 3 (∆η 3 ) 2 + ε(e θ,ε σσ (η)) 2 dy + 4µ    ε 3 3 α,β ∂ αβ η 3 2 0,ω + ε α,β e θ,ε αβ (η) 2 0,ω    , for all η ∈ V (ω).
Likewise, we associate the sum of the right hand sides of ( 5) and ( 6) with a linear and continuous form defined as follows:

(η) := ω p ε i η i dy - ω s ε α ∂ α η 3 dy for all η ∈ V (ω). (9)
The energy functional associated with the variational formulation in Problem (P ε (ω)) takes the following form:

J ε (η) = 1 2 b(η, η) -(η), for all η ∈ V H (ω) × K ε 3 (ω).
As a result, Problem (P ε (ω)) is equivalent to finding

ζ ε = (ζ ε H , ζ ε 3 ) ∈ V H (ω) × K ε 3 (ω) such that J ε (ζ ε ) = min{J ε (η); η ∈ V H (ω) × K ε 3 (ω)}.
The bilinear form b(•, •) is continuous, i.e., there exists a constant

M > 0 such that b(ξ, η) ≤ M ξ V (ω) η V (ω) for all ξ, η ∈ V (ω).
By Theorem 3.6-1 of [START_REF] Ciarlet | Mathematical Elasticity[END_REF], such a bilinear form b(•, •) is V (ω)-elliptic, i.e., there exists a constant α > 0 such that

b(η, η) ≥ α η 2 V (ω)
for all η ∈ V (ω).

As a result, problem P ε (ω) admits a unique solution

ζ ε = (ζ ε H , ζ ε 3 ) which belongs to V H (ω) × K ε 3 (ω) and satisfying b(ζ ε , η -ζ ε ) ≥ (η -ζ ε ) for all η = (η H , η 3 ) ∈ V H (ω) × K ε 3 (ω). ( 10
)

A finite element method for the obstacle problem

In this section we present a suitable finite element method to approximate the solution to problem P ε (ω). Following [START_REF] Ciarlet | The Finite Element Method for Elliptic Problems[END_REF] and [START_REF] Brenner | The mathematical theory of finite element methods[END_REF] (see also [START_REF] Chapelle | The finite element analysis of shells -fundamentals[END_REF], [START_REF] Chen | Current trends in scientific computing: ICM[END_REF], [START_REF] Ganesan | Finite Elements: Theory and Algorithms[END_REF] and [START_REF] Li | Finite element method and its applications[END_REF]), we recall some basic terminology and definitions. In what follows the letter h denotes a quantity approaching zero. For brevity, the same notation C (with or without subscripts) designates a positive constant independent of h, which can take different values at different places. We denote by (T h ) h>0 a family of triangulations of the polygonal domain ω made of triangles and we let T denote any element of such a family. Let us first recall, following [START_REF] Ciarlet | The Finite Element Method for Elliptic Problems[END_REF] and [START_REF] Brenner | The mathematical theory of finite element methods[END_REF], the rigorous definition of finite element in R n , where n ≥ 1 is an integer. A finite element in R n is a triple (T, P, N ) where:

(i) T is a closed subset of R n with non-empty interior and Lipschitz-continuous boundary, (ii) P is a finite dimensional space of real-valued functions defined over T , (iii) N is is a finite set of linearly independent linear forms N i , 1 ≤ i ≤ dim P , defined over the space P .

By definition, it is assumed that the set N is P -unisolvent in the following sense: given any real scalars α i , 1 ≤ i ≤ dim P , there exists a unique function g ∈ P which satisfies

N i (g) = α i , 1 ≤ i ≤ dim P.
It is henceforth assumed that the degrees of freedom, N i , lie in the dual space of a function space larger than P like, for instance, a Sobolev space (see [START_REF] Brenner | The mathematical theory of finite element methods[END_REF]). For brevity we shall conform our terminology to the one of [START_REF] Ciarlet | The Finite Element Method for Elliptic Problems[END_REF], calling the sole set T a finite element. Define the diameter of any finite element T as follows:

h T = diam T := max x,y∈T |x -y|.
Let us also define

ρ T := sup{diam B; B is a ball contained in T }.
A triangulation T h is said to be regular (cf., e.g., [START_REF] Ciarlet | The Finite Element Method for Elliptic Problems[END_REF]) if: (i) There exists a constant σ > 0, independent of h, such that

for all T ∈ T h , h T ρ T ≤ σ.
(ii) The quantity h := max{h T > 0; T ∈ T h } approaches zero.

A triangulation T h is said to satisfy an inverse assumption (cf., e.g., [START_REF] Ciarlet | The Finite Element Method for Elliptic Problems[END_REF]) if there exists a constant κ > 0 such that

for all T ∈ T h , h h T ≤ κ.
There is of course an ambiguity in the meaning of h, which was first regarded as a parameter associated with the considered family of triangulations, and which next denotes a geometrical entity. Nevertheless, in this paper, we have conformed to this standard notation (see [START_REF] Ciarlet | The Finite Element Method for Elliptic Problems[END_REF]). In the rest of this section, the parameter h is assumed to be fixed and we also assume that the triangulation T h under consideration is regular and satisfies the aforementioned inverse assumption. Let V h be the set of all of the nodal points of T h , let p denote any point of V h and let E h be the set of open edges of T h , in the sense that any edge e ∈ E h is isomorphic to the open interval (0, 1).

The forthcoming finite element analysis will be carried out using triangles of type (1) (see Figure 2.2.1 of [START_REF] Ciarlet | The Finite Element Method for Elliptic Problems[END_REF]) to approximate the tangential components of the displacement vector field and Morley's triangles (see [START_REF] Morley | The triangular equilibrium element in the solution of plate bending problems[END_REF] and also [START_REF] Ciarlet | The Finite Element Method for Elliptic Problems[END_REF]) to approximate the transverse component of the displacement vector field. In this case, the set V h consists of all the vertices and all the midpoints of the triangulation

T h . Let V 1,h and V 2,h be two finite dimensional spaces such that V 1,h × V 2,h ⊂ V H (ω)
and let (see ,e.g., [START_REF] Ciarlet | The Finite Element Method for Elliptic Problems[END_REF] and [START_REF] Brenner | A Morley finite element method for the displacement obstacle problem of clamped Kirchhoff plates[END_REF])

V 3,h := {η ∈ L 2 (ω); η T ∈ P 2 (T ), η is continuous at the vertices,
∂ ν η continuous at the midpoint of the edges} be the finite dimensional space associated with Morley's triangle. Define

V h := V 1,h × V 2,h × V 3,h .
We henceforth denote by η T the restriction of any function η ∈ L 2 (ω) to the finite element T . We denote by Ṽ3,h the subspace of V 3,h for which η(a k ) = 0, for all the vertices a k ∈ γ and

∂ ν η(b k ) = 0, for all the edges midpoints b k such that b k ∈ γ. Define the space Ṽh := V 1,h × V 2,h × Ṽ3,h .
Since Ṽ3,h is not contained in C 0 (ω) (see, e.g., [START_REF] Lascaux | Some nonconforming finite elements for the plate bending probelm[END_REF] and [START_REF] Ganesan | Finite Elements: Theory and Algorithms[END_REF]), we have

Ṽ3,h ⊂ V 3 (ω).
Define the space

H 2 (ω, T h ) := {η ∈ L 2 (ω); η T ∈ H 2 (T ) for all T ∈ T h } and introduce the semi-norm η ∈ H 2 (ω, T h ) → η h := T ∈T h |η| 2 2,T 1/2
, which becomes a norm over the space Ṽ3,h (cf. [START_REF] Brenner | A Morley finite element method for the displacement obstacle problem of clamped Kirchhoff plates[END_REF]). As a result, the mapping

η h ∈ Ṽh → η h := η 1,h 1,ω + η 2,h 1,ω + η 3,h h ,
is a norm over the space Ṽh .

Define the space

V 3 (ω) + Ṽ3,h := {ξ 3 = η 3 + η 3,h ; η 3 ∈ V 3 (ω) and η 3,h ∈ Ṽ3,h }.
Following [START_REF] Brenner | Finite element methods for the displacement obstacle problem of clamped plates[END_REF] and [START_REF] Ciarlet | The Finite Element Method for Elliptic Problems[END_REF], we define the approximate bilinear form b h (•, •), associated with the bilinear form b defined in [START_REF] Brezzi | Error estimates for the finite element solution of variational inequalities[END_REF], as follows:

b h : V H (ω) × (V 3 (ω) + Ṽ3,h ) × V H (ω) × (V 3 (ω) + Ṽ3,h ) → R is such that b h | V (ω)×V (ω) = b, i.e., b h (ξ, η) = b(ξ, η) for all ξ, η ∈ V (ω) and such that b h (ξ h , η h ) = - T ∈T h T m ε αβ (ξ 3,h )∂ αβ η 3,h dy + T ∈T h T n θ,ε αβ (ξ h )(∂ α θ ε )∂ β η 3,h dy + T ∈T h T n θ,ε αβ (ξ h )∂ β η α,h dy, for all ξ h , η h ∈ Ṽh .
Therefore, the bilinear form b h (•, •) is continuous over Ṽh , i.e., there exists M > 0, independent of h, such that

b h (ξ, η) ≤ M ξ η for all ξ, η ∈ Ṽh .
Besides, in view of Theorem 3.4-1 of [START_REF] Ciarlet | Mathematical Elasticity[END_REF] and the theory presented in [START_REF] Ciarlet | The Finite Element Method for Elliptic Problems[END_REF], the bilinear form b h (•, •) is Ṽh -elliptic, namely, there exists α > 0, independent of h, such that b h (η, η) ≥ α η 2 for all η ∈ Ṽh .

Let us now define the

V h interpolation operator Π h : C 0 (ω)×C 0 (ω)×H 2 (ω) → V h as follows Π h ξ := (Π 1,h ξ 1 , Π 2,h ξ 2 , Π 3,h ξ 3 ) for all ξ ∈ C 0 (ω) × C 0 (ω) × H 2 (ω),
where Π i,h is the standard V i,h interpolation operator (cf., e.g., [START_REF] Ciarlet | The Finite Element Method for Elliptic Problems[END_REF] and [START_REF] Brenner | The mathematical theory of finite element methods[END_REF]). It thus results that the interpolation operator Π h satisfies the following properties (Π j,h ξ j )(p) = ξ j (p) for all integers 1 ≤ j ≤ 3 and all vertices p ∈ V h , where ν e is outer unit normal vector to the edge e. Define the space

H(ω) := H 2 (ω) × H 2 (ω) × H 3 (ω)
and equip it with the norm

ξ ω := ξ 1 2,ω + ξ 2 2,ω + ξ 3 3,ω for all ξ ∈ H(ω).
An application of Theorem 3.2.1 of [START_REF] Ciarlet | The Finite Element Method for Elliptic Problems[END_REF] (see also Theorem 4.4.20 of [START_REF] Brenner | The mathematical theory of finite element methods[END_REF]) yields

ξ -Π h ξ ≤ Ch|ξ| ω , (11)
for all ξ ∈ H(ω) ∩ V (ω).

In order to provide the required estimates for the convergence of the numerical scheme, we make use of enriching operators. Enriching operators were first introduced in [START_REF] Brenner | Two-level additive Schwarz preconditioners for nonconforming finite element methods[END_REF] (see also [START_REF] Brenner | A two-level additive Schwarz preconditioner for nonconforming plate elements[END_REF], [START_REF] Brenner | Convergence of nonconforming multigrid methods without full elliptic regularity[END_REF] and [START_REF] Brenner | C 0 interior penalty method for fourth order elliptic boundary value problems on polygonal domains[END_REF]) and they play a key role in the study of obstacle problems for clamped plates (see [START_REF] Brenner | Finite element methods for the displacement obstacle problem of clamped plates[END_REF] and [START_REF] Brenner | A Morley finite element method for the displacement obstacle problem of clamped Kirchhoff plates[END_REF]). Following Example 2.2 of [START_REF] Brenner | Finite element methods for the displacement obstacle problem of clamped plates[END_REF], we recall that any enriching operator associated with conforming finite elements coincides with the canonical injection. We are going to connect Morley's triangle to the Hsieh-Clough-Tocher macro-element (from now onwards HCT macro-element), that we sketch below for sake of clarity, via an ad hoc enriching operator (for a complete overview on the properties of these finite elements and the meaning of the graphical symbols used for representing the various degrees of freedom, see Figures 6.1.3 and 6.2.3 of [START_REF] Ciarlet | The Finite Element Method for Elliptic Problems[END_REF]). The relation between the elements in Fig. 1 and Fig. 2 is due to the disposition of the vertices at which the pointwise evaluation of the shape functions occurs. Let us denote W 3,h the finite element space associated with the HCT macroelement and let us observe that, by the unisolvence of the HCT macro-element (cf. Theorem 6.1.2 of [START_REF] Ciarlet | The Finite Element Method for Elliptic Problems[END_REF]), the elements of W 3,h are completely determined by their values at the vertices, the values of their first derivatives at the vertices and the values of their normal derivatives at the midpoints of the sides of the triangular element. The reason why we have to make use of a nonconforming finite element to carry out the numerical analysis of the considered obstacle problem is due to the fact that it is not natural to assume the transverse component of the solution, i.e., the one which is affected by the geometrical constraint associated with the obstacle, to be more regular than H 3 (ω) (see, for instance, [START_REF] Frehse | Zum differenzierbarkeitsproblem bei variationsungleichungen höherer ordnung[END_REF] and [START_REF] Frehse | On the regularity of the solution of the biharmonic variational inequality[END_REF], [START_REF] Caffarelli | The obstacle problem for the biharmonic operator[END_REF], [START_REF] Caffarelli | The two-obstacle problem for the biharmonic operator[END_REF], [START_REF] Friedman | Variational principles and free-boundary problems[END_REF] and [START_REF] Brenner | Finite element methods for the displacement obstacle problem of clamped plates[END_REF]).

Let us thus define the enriching operator

E h : V 3,h → W 3,h by (cf. formula (3.2) of [5]) [N (E h η)] = 1 |T p | T ∈T p (N η T ), (12)
where p ∈ V h is any nodal point of the triangulation T h , N is any degree of freedom of the HCT macro-element associated with the nodal point p and T p is the set of triangles in T h sharing the nodal point p. Next, following [START_REF] Brenner | Convergence of nonconforming multigrid methods without full elliptic regularity[END_REF], we organise the proof of the already well-known properties of enriching operators in a series of lemmas (Lemmas 1-4).

Let us recall the definition of jump of the normal derivative across the edge e. Let η ∈ H 2 (ω) and let e ∈ E h such that e ⊂ ω. The jump of the normal derivative of η across the edge e is defined as follows

∂ ν η := ∂η + ∂ν e e - ∂η - ∂ν e e , (13)
where T + and T -are elements of T h that share the edge e, η ± is the restriction of η to T ± and ν e points from T + to T -(see Fig. 3 below). If e ⊂ γ, then we define the jump in this fashion:

∂ ν η := - ∂η ∂ν e e . ( 14
)
The proof of the next lemma relies on standard inverse estimates (cf., e.g., [START_REF] Ciarlet | The Finite Element Method for Elliptic Problems[END_REF]) and inverse trace inequalities (cf. formula (10.3.9) of [START_REF] Brenner | The mathematical theory of finite element methods[END_REF]).

Lemma 1 There exists a positive constant C such that

η -E h η 0,ω ≤ Ch 2 η h for all η ∈ V 3,h , (15)
|E h η| 2,ω ≤ C η h , for all η ∈ V 3,h . (16)
Proof Let us fix an arbitrary T ∈ T h and let N denote the set of degrees of freedom of the HCT macro-element T (cf. Fig. 2). For any η ∈ V 3,h , we have that (η -E h η)| K i is an element of P 3 (K i ), for all 1 ≤ i ≤ 3, where K i is a subtriangle of the HCT macro-element T (cf. Fig. 2). Using the same argument as in Theorem 3.1.5 of [START_REF] Ciarlet | The Finite Element Method for Elliptic Problems[END_REF], we infer the existence of a positive constant C for which

ξ 2 0,T ≤ Ch 2(1+|N |) N ∈N |N (ξ)| 2 , (17)
for all ξ ∈ C 1 (T ) such that ξ| K i ∈ P 3 (K i ), for all 1 ≤ i ≤ 3, where K i is a subtriangle of the HCT macro-element T (cf. Fig. 2) and |N | denotes the order of differentiation of the corresponding degree of freedom. In view of Remark 3.1.3 of [START_REF] Ciarlet | The Finite Element Method for Elliptic Problems[END_REF], it results that the diameter of T is of order O(h) and, therefore, |e| = O(h) as well. By [START_REF] Chapelle | The finite element analysis of shells -fundamentals[END_REF] and the continuity of η at the vertices of the triangulation, it results

N (η) = N (E h η) if |N | = 0. Therefore, letting ξ = η -E h η in (17) yields η -E h η 2 0,T ≤ Ch 4 N ∈N |N |=1 |N (η -E h η)| 2 . ( 18
)
Let us observe that if N is associated with the degree of freedom corresponding to the outer unit normal derivative at the midpoint of a side of the boundary γ then N (η -E h η) = 0.

Let N denote the degree of freedom corresponding to the evaluation of the outer unit normal derivative at the midpoint m e of an edge e ⊂ ω and let ν e denote one of the outer unit normal vectors to the edge e. By virtue of ( 12), a standard inverse estimate (Theorem 3.2.6 of [START_REF] Ciarlet | The Finite Element Method for Elliptic Problems[END_REF] with q = ∞, m = l = 2 and r = 2) and an inverse trace inequality, we obtain

|N (η -E h η)| 2 = |∂ ν e (η -E h η) (m e )| 2 = [∂ ν e η(m e ) -∂ ν e E h η(m e )] 2 = 1 2 ∂ ν e η + (m e ) - 1 2 ∂ ν e η -(m e ) 2 ≤ |e| 4 ∂ ν e η 0,e ≤ C T ∈T p |η| 2 2,T . (19)
Let N be the degree of freedom associated with the evaluation of any first order derivative at any vertex p ∈ V h . An arithmetic-geometric mean inequality yields

|N (η -E h η)| 2 ≤ |∇(η -E h η)(p)| 2 ≤ C T ,T ∈T p T and T share an edge |∇η T (p) -∇η T (p)| 2 . ( 20
)
An application of the mean value theorem (cf., e.g., Theorem 7.2-1 of [START_REF] Ciarlet | Linear and Nonlinear Functional Analysis with Applications[END_REF]) like on page 915 of [START_REF] Brenner | Two-level additive Schwarz preconditioners for nonconforming finite element methods[END_REF], standard inverse estimates (Theorem 3.2.6 of [START_REF] Ciarlet | The Finite Element Method for Elliptic Problems[END_REF] with r = m = l = 2 and q = ∞), an inverse trace inequality and the regularity of the triangulation, gives

|∇η T (p) -∇η T (p)| 2 ≤ C|∂ ν e η T (p) -∂ ν e η T (p)| 2 + C|∂ τ e η T (p) -∂ τ e η T (p)| 2 ≤ C|e| -1 ∂ ν η 2 0,e + C|e| 2 (|η T | 2 2,∞,T + |η T | 2 2,∞,T ) ≤ C T ∈T e |η| 2 2,T , (21)
where T e is the set of triangles sharing the edge e. Combining ( 18)-( 21), we obtain

η -E h η 2 0,T ≤ Ch 4 T ∈T T |η| 2 2,T . (22)
Estimate [START_REF] Ciarlet | Mathematical Elasticity. Vol. I: Three-Dimensional Elasticity[END_REF] follows by summing up ( 22) over all the triangles of T h . Estimate [START_REF] Ciarlet | Mathematical Elasticity[END_REF] follows by standard inverse estimates (Theorem 3.2.6 of [START_REF] Ciarlet | The Finite Element Method for Elliptic Problems[END_REF] with m = 2, l = 0, p = r = 2) and ( 22): Indeed,

|E h η| 2 2,ω ≤ C T ∈T h |η -E h η| 2 2,T + T ∈T h |η| 2 2,T ≤ T ∈T h h -4 η -E h η 2 0,T + |η| 2 2,T ≤ C η 2 h .
This completes the proof.

By virtue of an interpolation estimate (see, for instance, Theorem 3.1.5 of [START_REF] Ciarlet | The Finite Element Method for Elliptic Problems[END_REF] with m = q = p = k = 2), which holds by the fact that Morley's triangle is almost affine (cf., e.g., [START_REF] Lascaux | Some nonconforming finite elements for the plate bending probelm[END_REF]), and the standard trace theorem for Sobolev spaces defined over domains, we deduce that, for all T ∈ T h ,

∇(η -Π 3,h η) 0,∂T ≤ Ch|η| 3,T , (23)
where η ∈ H 3 (T ) ∩ V 3 (ω). By [START_REF] Chapelle | The finite element analysis of shells -fundamentals[END_REF], we deduce that E h η = η at the internal nodes of the triangulation (see also formula (6.11) of [START_REF] Brenner | Convergence of nonconforming multigrid methods without full elliptic regularity[END_REF]). As a result,

Π 3,h E h η = η for all η ∈ V 3,h .
The next preliminary result is inspired by Lemma 2 of [START_REF] Brenner | C 0 interior penalty method for fourth order elliptic boundary value problems on polygonal domains[END_REF], which is itself based on the unisolvence of the HCT macro-element (see, for instance, Theorem 6.1.2 of [START_REF] Ciarlet | The Finite Element Method for Elliptic Problems[END_REF]) and Bramble-Hilbert lemma (cf., e.g., Theorem 4.1.3 of [START_REF] Ciarlet | The Finite Element Method for Elliptic Problems[END_REF]). For convenience, we provide a complete proof.

Lemma 2 There exists a positive constant C solely depending on the regularity of the triangulation T h such that

2 m=0 h 2m |η -E h Π 3,h η| 2 m,T ≤ Ch 6 |η| 2 3,S T , (24)
for all T ∈ T h and all η ∈ H 3 (S T ) ∩ H 2 (ω), where S T is the polygon formed by all the triangles of T h sharing a vertex with T (cf. Fig. 4

below).

Proof Let T ∈ T h be an arbitrary element. Then we observe that (see Lemma 2 of [START_REF] Brenner | C 0 interior penalty method for fourth order elliptic boundary value problems on polygonal domains[END_REF]) the expression (E h Π 3,h η)| T is completely determined by η| S T and that the mapping

η| S T → (η -E h Π 3,h η)| T , is bounded from H 3 (S T ) into H 2 (T ).
Fig. 4: The polygon S T made of all the triangles of T T . Figure 9 of [START_REF] Brenner | A two-level additive Schwarz preconditioner for nonconforming plate elements[END_REF] Moreover, [START_REF] Chapelle | The finite element analysis of shells -fundamentals[END_REF] and the unisolvence of the HCT macro-element give q -E h Π 3,h q = 0 for all q ∈ P 2 (T ).

Thanks to [START_REF] Falk | Error estimates for the approximation of a class of variational inequalities[END_REF], we can apply Bramble-Hilbert lemma and infer the validity of [START_REF] Ciarlet | Obstacle problems for Koiter's shells[END_REF]. The proof is thus complete.

We now modify the definition of the enriching operator E h in order to incorporate the boundary conditions. As a result, we obtain the corresponding enriching operator Ẽh : Ṽ3,h → W3,h , where Ṽ3,h denotes the subspace of V 3,h whose degrees of freedom vanish along γ and W3,h := W 3,h ∩ H 2 0 (ω), i.e., the subspace of W 3,h whose degrees of freedom vanish along γ (cf. Example 6.1 of [START_REF] Brenner | Convergence of nonconforming multigrid methods without full elliptic regularity[END_REF] and [START_REF] Brenner | A Morley finite element method for the displacement obstacle problem of clamped Kirchhoff plates[END_REF]). Observe that E h η ∈ H 2 (ω) by the properties of the HCT macro-element. We now define Ẽh as follows, in such a way that Ẽh η ∈ H 2 0 (ω) for all η ∈ Ṽ3,h :

(i) The degrees of freedom of E h η and Ẽh η coincide in ω, (ii) The degrees of freedom of Ẽh η vanish on γ.

In the next lemma we prove the first properties of the modified enriching operator Ẽh . The proof resorts to standard inverse estimates (cf., e.g., [START_REF] Ciarlet | The Finite Element Method for Elliptic Problems[END_REF]), an inverse trace inequality (cf. formula (10.3.9) of [START_REF] Brenner | The mathematical theory of finite element methods[END_REF]) and Lemma 1.

Lemma 3 For each η ∈ Ṽ3,h there exists a positive constant C such that η -Ẽh η 0,ω ≤ Ch 2 η h , ( 26
) Ẽh η 2,ω ≤ C η h . ( 27
)
Proof Let us fix an arbitrary element η ∈ Ṽ3,h and let N be any degree of freedom associated with the HCT macro-element. Let us observe that if N is not related to any nodal point of γ, then, by property (i) in the definition of Ẽh , we obtain

N (E h η -Ẽh η) = 0.
On the other hand, for all T ∈ T h , if |N | = 0 and N is related to a nodal point of γ, then, by property (ii) in the definition of Ẽh and the fact that η ∈ Ṽ3,h , it follows that

N (E h η -Ẽh η) = N (E h η) = 0.
Let m e be the midpoint of an edge e ⊂ γ. By ( 14) and standard inverse estimates (Theorem 3.2.6 of [START_REF] Ciarlet | The Finite Element Method for Elliptic Problems[END_REF] with m = l = 0, q = ∞ and r = 2) we obtain

|∂ ν e (E h η -Ẽh η)(m e )| 2 = |∂ ν e η(m e )| 2 ≤ C |e| ∂ ν η 2 0,e ≤ C T ∈T e |η| 2 2,T , (28)
where the latter inequality holds true by virtue of an inverse trace inequality and the regularity of the triangulation T h .

Let p be a vertex on γ. Then p is the endpoint of an edge e * ⊂ γ. Let T * be an element of T h such that e * ⊂ T * . By ( 20), ( 21), an arithmetic-geometric mean inequality, standard inverse estimates (Theorem 3.2.6 of [START_REF] Ciarlet | The Finite Element Method for Elliptic Problems[END_REF] with m = l = r = 2 and q = ∞) and the regularity of the triangulation, we get

|∇(E h η -Ẽh η)(p)| 2 = |∇(E h η)(p)| 2 ≤ C(|∇(E h η -η T * )(p)| 2 + |∇η T * (p)| 2 ) ≤ C T ∈T T |η| 2 2,T , (29)
where E p is the set of edges of E h sharing p as a common vertex.

In conclusion, for all T ∈ T h , a combination of ( 22), ( 28) and ( 29) yields

η -Ẽh η 2 0,T ≤ Ch 4 T ∈T T |η| 2 2,T , (30)
which implies the inequality (26) after a summation over all the elements of T h . Using Poincaré-Friedrichs inequality, an arithmetic-geometric mean inequality, (26) and standard inverse estimates (Theorem 3.2.6 of [START_REF] Ciarlet | The Finite Element Method for Elliptic Problems[END_REF] with m = 2, l = 0, q = r = 2), we infer the validity of [START_REF] Frehse | On the regularity of the solution of the biharmonic variational inequality[END_REF]. Indeed,

Ẽh η 2 2,ω ≤ C| Ẽh η| 2 2,ω ≤ T ∈T h |η -Ẽh η| 2 2,T + |η| 2 2,T ≤ T ∈T h h -4 η -Ẽh η 2 0,T + |η| 2 2,T ≤ C η 2 h ,
which completes the proof.

By property (i) in the definition of Ẽh , we can easily observe that the following holds (cf. [START_REF] Brenner | C 0 interior penalty method for fourth order elliptic boundary value problems on polygonal domains[END_REF]):

Π 3,h Ẽh η = η for all η ∈ Ṽ3,h .

The next step consists in incorporating the boundary conditions into the above estimates.

Lemma 4 There exists a positive constant C solely depending on the regularity of the triangulation T h such that

2 m=0 h 2m |η -Ẽh Π 3,h η| 2 m,T ≤ Ch 6 |η| 2 3,S T , (31)
for all T ∈ T h and all η ∈ H 3 (S T ) ∩ V 3 (ω).

Proof Let us fix an arbitrary element η ∈ Ṽ3,h and let N be any degree of freedom associated with an internal vertex of the triangulation made of HCT macroelement. Like in Lemma 3, we have

N (E h η -Ẽh η) = 0.
Let e ∈ E h be an edge contained in γ and let m e be the midpoint of e, at which the normal derivative of η is evaluated. By virtue of ( 14), ( 23), ( 28) and a standard inverse estimate (Theorem 3.1.5 of [START_REF] Ciarlet | The Finite Element Method for Elliptic Problems[END_REF] with m = q = k = p = 2), we obtain

|∂ ν e (E h Π 3,h η -Ẽh Π 3,h η)(m e )| 2 = |∂ ν e (E h Π 3,h η)(m e )| 2 = |∂ ν e (Π 3,h η)(m e )| 2 ≤ C|e| -1 ∂ ν Π 3,h η 2 0,e ≤ Ch 2 T ∈T e |η| 2 3,T , (32)
for all η ∈ H 3 (S T ) ∩ V 3 (ω).

Similarly, for any vertex p ∈ γ, we have, by [START_REF] Ganesan | Finite Elements: Theory and Algorithms[END_REF], standard interpolation estimates (Theorem 3.1.5 of [START_REF] Ciarlet | The Finite Element Method for Elliptic Problems[END_REF]) and an inverse trace inequality (formula (10.3.9) of [START_REF] Brenner | The mathematical theory of finite element methods[END_REF])

|∇(E h Π 3,h η -Ẽh Π 3,h η)(p)| 2 ≤ T ∈T p |∂ τ e (E h Π 3,h η T )(p)| 2 + T ∈T p |∂ ν e (E h Π 3,h η T )(p)| 2 ≤ C T ∈T p |∂ τ e (Π 3,h η T -η T )(p)| 2 + C e∈E p |e| -1 ∂ ν Π 3,h η 2 0,e ≤ C T ∈T p |Π 3,h η T -η T | 2 2,T + C e∈E p |e| -1 ∂ ν (Π 3,h η -η) 2 0,e ≤ Ch 2 T ∈T p |η| 2 3,T . (33)
Summing over all the triangles of T h , we obtain the following estimate

2 m=0 h 2m |η -E h Π 3,h η| 2 m,T ≤ Ch 4 N ∈N |N |=1 |N (η -E h Π 3,h η)| 2 + 2 m=1 h 2m |η -E h Π 3,h η| 2 m,T ≤ Ch 6 |η| 2 3,S T , (34)
where the first inequality holds by [START_REF] Ciarlet | An Introduction to Differential Geometry with Applications to Elasticity[END_REF] and the latter inequality holds by [START_REF] Lascaux | Some nonconforming finite elements for the plate bending probelm[END_REF] and [START_REF] Léger | Mathematical justification of the obstacle problem in the case of a shallow shell[END_REF]. This completes the proof.

As a consequence of Lemmas 1-4 and standard inverse estimates (Theorem 3.2.6 of [START_REF] Ciarlet | The Finite Element Method for Elliptic Problems[END_REF] with m = 1, l = 0, q = r = 2), we obtain the following estimates for the enriching operator Ẽh (cf. Corollary 1 of [START_REF] Brenner | C 0 interior penalty method for fourth order elliptic boundary value problems on polygonal domains[END_REF]):

η -Ẽh η 0,ω + h T ∈T h |η -Ẽh η| 2 1,T 1/2 + h 2 | Ẽh η| 2,ω ≤ Ch 2 η h for all η ∈ Ṽ3,h , ( 35
) 2 m=0 h m |η -Ẽh Π 3,h η| m,ω ≤ Ch 3 |η| 3,ω for all η ∈ H 3 (ω) ∩ V 3 (ω). ( 36
)
Define the space

Wh := V 1,h × V 2,h × W3,h ,
and let us define the enriching operator Ẽh : Ṽh → Wh as follows:

Ẽh ξ := (ξ 1 , ξ 2 , Ẽh ξ 3 ) for all ξ ∈ Ṽh .

A direct application of ( 35) and ( 36) to (37) yields

η -Ẽh η 0,ω + h T ∈T h |η -Ẽh η| 2 1,T 1/2 + h 2 | Ẽh η| 2,ω ≤ Ch 2 η for all η ∈ Ṽh . (38)
Next, we prove a crucial estimate for b h (•, •) in the case where the transverse component of the displacement is approximated via Morley's triangles. The assumption that the solution ζ ε to Problem P ε (ω) is "more regular" is of paramount importance.

By virtue of the results proved in [START_REF] Caffarelli | The obstacle problem for the biharmonic operator[END_REF], [START_REF] Caffarelli | The two-obstacle problem for the biharmonic operator[END_REF], [START_REF] Frehse | On the regularity of the solution of the biharmonic variational inequality[END_REF] and [START_REF] Friedman | Variational principles and free-boundary problems[END_REF] and in order to make our analysis more general, we will derive the sought error estimate under the constraint that the transverse component of the solution of Problem P ε (ω) cannot be more regular than H 3 (ω).

In order to derive error estimate, we will have to assume that the solution of Problem P ε (ω) is more regular (cf., e.g., [START_REF] Ciarlet | The Finite Element Method for Elliptic Problems[END_REF]); in particular, we will assume that

ζ ε ∈ H(ω) ∩ V (ω).
The augmented regularity result for the tangential components is studied, for instance, in Section 8.7 of Chapter 2 of [START_REF] Lions | Quelques méthodes de résolution des problèmes aux limites non linéaires[END_REF], while the augmented regularity result for the transverse component is given for solutions of some fourth-order variational inequalities on pages 323-327 of [START_REF] Glowinski | Numerical Analysis of Variational Inequalities[END_REF], and is also recalled in [START_REF] Wang | C¡sup¿0¡/sup¿ discontinuous galerkin methods for a plate frictional contact problem[END_REF].

To prove the next result we follow Appendix B of [START_REF] Brenner | Finite element methods for the displacement obstacle problem of clamped plates[END_REF] and Lemma 4.2 of [START_REF] Brenner | A Morley finite element method for the displacement obstacle problem of clamped Kirchhoff plates[END_REF]. As a consequence of the trace properties (cf., e.g., Theorem 6.6-5 of [START_REF] Ciarlet | Linear and Nonlinear Functional Analysis with Applications[END_REF]) we can take into account the average along any edge e ∈ E h of a function f ∈ H 1 (ω) and denote it by f , viz.,

f := 1 |e| e f ds ∈ R.
Lemma 5 There exists a positive constant C such that the following estimate holds

|b h (ζ ε , η -Ẽh η)| ≤ Ch ζ ε ω η for all η ∈ Ṽh , ( 39
)
where

ζ ε ∈ H(ω) ∩ V (ω) is the solution to problem P ε (ω).
Proof Observe that if η 3 ∈ Ṽ3,h then ∇η 3 is continuous at the midpoints of the edges e ∈ E h and vanishes at the midpoints of the edges along γ. Indeed, after fixing an edge e ∈ E h , consider the restrictions η 3 | T + and η 3 | T -to the edge e, where T ± are, again, the elements of T h that share the edge e. Then ( Combining the definition and the properties of Ẽh , Green's formula, the midpoint rule, Cauchy-Schwarz inequality, inverse trace inequalities (formulas (10.3.8) and (10.3.9) of [START_REF] Brenner | The mathematical theory of finite element methods[END_REF]), Poincaré-Friedrichs inequality (Theorems 6.5-2 and 6.8-1 of [START_REF] Ciarlet | Linear and Nonlinear Functional Analysis with Applications[END_REF]) and [START_REF] Morley | The triangular equilibrium element in the solution of plate bending problems[END_REF], we obtain

η 3 | T + -η 3 | T -) ∈ P 2 (R 2 ) and (η 3 | T + -η 3 | T -)
b h (ζ ε , η -Ẽh η) = - T ∈T h T m ε αβ (ζ ε 3 )∂ αβ (η 3 -Ẽh η 3 ) dy + T ∈T h T (∂ α θ ε )n θ,ε αβ (ζ ε )∂ β (η 3 -Ẽh η 3 ) dy = T ∈T h T ∂ α (m ε αβ (ζ ε 3 ))∂ β (η 3 -Ẽh η 3 ) dy - T ∈T h ∂T m ε αβ (ζ ε 3 )∂ β (η 3 -Ẽh η 3 )ν α ds + T ∈T h T (∂ α θ ε )n θ,ε αβ (ζ ε )∂ β (η 3 -Ẽh η 3 ) dy = T ∈T h T ∂ α (m ε αβ (ζ ε 3 )) + (∂ α θ ε )n θ,ε αβ (ζ ε ) ∂ β (η 3 -Ẽh η 3 ) dy - e∈E h e m ε αβ (ζ ε 3 ) -m ε αβ (ζ ε 3 ) ∂ β (η 3 -Ẽh η 3 )ν α ds ≤ T ∈T h ∂ α m ε αβ (ζ ε 3 ) + (∂ α θ ε )n θ,ε αβ (ζ ε ) 0,T ∂ β (η 3 -Ẽh η 3 ) 0,T + e∈E h |e| -1 m ε αβ (ζ ε 3 ) -m ε αβ (ζ ε 3 ) 0,e 1/2 × e∈E h |e| ∂ β (η 3 -Ẽh η 3 )ν α 0,e 1/2 ≤ C ζ ε ω T ∈T h |η 3 -Ẽh η 3 | 2 1,T 1/2 ≤ Ch ζ ε ω η ,
where, in analogy with ( 13)and ( 14), we have

∂ β (η 3 -Ẽh η 3 )ν α := (∂ β (η 3 -Ẽh η 3 )ν α,+ ) + (∂ β (η 3 -Ẽh η 3 )ν α,-).
The proof is thus complete.

Let us observe that, by virtue of the definition of Ẽh (cf. ( 37)), the variational equations ( 6) do not give any contribution in the previous proof .

Having extended the properties of the enriching operator to our problem, we now prove, following [START_REF] Brenner | Finite element methods for the displacement obstacle problem of clamped plates[END_REF], a series of preparatory lemmas. Let us recall that, for all h > 0, the symbol V h designates the set of all of the nodal points of the triangulation T h . Let us define the functional J h and the set K ε 3,h as follows:

J h (η h ) := 1 2 b h (η h , η h ) -(η h ) for all η h ∈ V h , (40)
K ε 3,h := {η 3,h ∈ Ṽ3,h ; θ ε (p) + η 3,h (p) ≥ 0 for all p ∈ V h }, (41)
and let us then state the approximate problem P ε h corresponding to problem P ε (ω).

Problem P ε h To find ζ ε h ∈ Ṽh such that the transverse component ζ ε 3,h belongs to K ε 3,h and such that J h (ζ ε h ) = inf η h ∈ Ṽh η 3,h ∈K ε 3,h J h (η h ). ( 42
)
Since the bilinear form b h (•, •) is symmetric, continuous over the space Ṽh and it is Ṽh -elliptic, we infer that problem P ε h has a unique solution ζ ε h , which satisfies the variational inequalities

b h (ζ ε h , η h -ζ h ) ≥ (η h -ζ ε h ), (43)
for all η h ∈ Ṽh such that η 3,h ∈ K ε 3,h .

Lemma 6 Let ζ ε and ζ ε h respectively denote the solutions to problem P ε (ω) and problem P ε h . There exist two constants C 1 > 0 and C 2 > 0 such that

ζ ε -ζ ε h 2 ≤ C 1 ζ ε -Π h ζ ε 2 + C 2 [b h (ζ ε , Π h ζ ε -ζ ε h ) -(Π h ζ ε -ζ ε h )] . (44)
Proof Observe that Π h ζ ε belongs to the space Ṽh . By the continuity and the Ṽh -ellipticity of b h (•, •) and Young's inequality (see [START_REF] Young | On Classes of Summable Functions and their Fourier Series[END_REF]) we get

α Π h ζ ε -ζ ε h 2 ≤ b h (Π h ζ ε -ζ ε h , Π h ζ ε -ζ ε h ) ≤ M Π h ζ ε -ζ ε Π h ζ ε -ζ ε h + b h (ζ ε , Π h ζ ε -ζ ε h ) -(Π h ζ ε -ζ ε h ) ≤ M 2 M α Π h ζ ε -ζ ε 2 + α M Π h ζ ε -ζ ε h 2 + b h (ζ ε , Π h ζ ε -ζ ε h ) -(Π h ζ ε -ζ ε h ).
Letting C 1 := M 2 /α 2 and C 2 := α -1 , we obtain inequality (44).

An intermediary problem

In what follows, we shall estimate the term

[b h (ζ ε , Π h ζ ε -ζ ε h ) -(Π h ζ ε -ζ ε h )
] in order to apply the interpolation estimate [START_REF] Carstensen | Nonconforming fem for the obstacle problem[END_REF]. To this aim, we introduce an intermediary problem, since it is not easy to directly connect

K ε 3 (ω) to K ε 3,h . Define the set Kε 3,h (ω) := {η 3 ∈ H 2 0 (ω); θ ε (p) + η 3 (p) ≥ 0 for all p ∈ V h } (45)
and define

V h (ω) := H 1 0 (ω) × H 1 0 (ω) × Kε 3,h (ω) as to define the functional J : V h (ω) → R by J(η) := 1 2 b(η, η) -(η). (46)
Let us state the intermediary problem P ε h (ω) establishing the connection between problem P ε (ω) and problem

P ε h . Problem P ε h (ω) To find ζε h ∈ V h (ω) such that the transverse component ζε 3,h
belongs to Kε 3,h (ω) and such that

J( ζε h ) = inf η∈V h (ω) η 3 ∈ Kε 3,h (ω) J(η). (47)
Using properties (i) and (ii) of the enriching operator Ẽh , we immediately obtain that Ẽh η 3,h ∈ Kε 3,h (ω) for all η 3,h ∈ K ε 3,h . Using the symmetry and the continuity and the V (ω)-ellipticity of the bilinear form b(•, •), we infer that problem P ε h (ω) admits one and only one solution ζε h satisfying the following variational inequalities:

b( ζε h , η -ζε h ) ≥ (η -ζε h ) for all η ∈ V h (ω). ( 48
)
The aim of the next lemma, whose formulation is inspired by Lemma 3.1 of [START_REF] Brenner | Finite element methods for the displacement obstacle problem of clamped plates[END_REF], is to prove that the uniform boundedness of the family ( ζε h ) h>0 , where ζε h denotes the solution to problem P ε h (ω). The argument resorts on Cauchy-Schwarz inequality, Poincaré-Friedrichs inequality (Theorems 6.5-2 and 6.8-1 of [START_REF] Ciarlet | Linear and Nonlinear Functional Analysis with Applications[END_REF]) and Young's inequality (see [START_REF] Young | On Classes of Summable Functions and their Fourier Series[END_REF]).

Lemma 7 There exists a constant C > 0 such that ζε h V (ω) ≤ C for all h > 0.

(49)

Proof Fix h > 0. Since K ε 3 (ω) ⊂ Kε 3,h (ω), we infer that J( ζε h ) ≤ J(ζ ε ), where ζε h and ζ ε are respectively the solutions to problem P ε h (ω) and problem P ε (ω). Using Cauchy-Schwarz inequality, Poincaré-Friedrichs inequality and Young's inequality, we obtain

α 2 ζε h 2 V (ω) ≤ 1 2 b( ζε h , ζε h ) = J( ζε h ) + ( ζε h ) ≤ J( ζε h ) + C | ζε h | V (ω) ≤ J( ζε h ) + 1 α C 2 + α 4 | ζε h | 2 V (ω) ,
which in turn implies that ζε

h 2 V (ω) ≤ 4 α J( ζε h ) + C 2 α ,
from which the estimate (49) immediately follows.

The purpose of the following lemmas, whose formulations are respectively inspired by those of Lemmas 3.2-3.4 of [START_REF] Brenner | Finite element methods for the displacement obstacle problem of clamped plates[END_REF], is to estimate the distance between ζε h and

ζ ε ∈ V H (ω) × K ε 3 (ω).
In what follows, the symbol denotes weak convergences as h → 0. Strong convergences in the space C 0 (ω) are meant with respect to the sup-norm.

Lemma 8

The following convergences take place

ζε h ζ ε in V (ω), (50) ζε 3,h → ζ ε 3 in C 0 (ω). ( 51
)
Proof The uniform boundedness of ( ζε h ) h>0 proved in Lemma 7, yields the existence of an element ζ * ∈ V (ω) such that, up to passing to a subsequence, still denoted ( ζε

h ) h>0 ζε h ζ * in V (ω).
The functional J is clearly sequentially weakly lower semi-continuous. Hence,

J(ζ * ) ≤ lim inf h→0 J( ζε h ) ≤ J(ζ ε ),
where the latter inequality is derived in Lemma 7. By the Rellich-Kondrachov theorem, we infer that ζ * 3 ∈ C 0 (ω) and that ζε

3,h → ζ * 3 in C 0 (ω).
It remains to prove ζ * = ζ ε . To this end, by the uniqueness of the solution to problem P ε (ω), it suffices to show that

ζ * 3 ∈ K ε 3 (ω). Since ζε 3,h ∈ Kε 3,h (ω), then θ ε (p) + ζε 3,h (p) ≥ 0 for all p ∈ V h .
Besides, the following density with respect to the Euclidean norm

h>0 V h = ω,
yields, in conjunction with the previous inequality, that

θ(q) + ζ * 3 (q) = lim k→∞ q k ∈V h k θ(q k ) + ζ * 3 (q k ) = lim k→∞ q k ∈V h k lim h→0 θ(q k ) + ζ3,h (q k ) ≥ 0,
for all q ∈ ω. We have thus shown that ζ * 3 ∈ K ε 3 (ω). The convergences (50) and (51) immediately follow.

Let us denote C the contact zone for problem P ε (ω), i.e.,

C := {y ∈ ω; θ ε (y) + ζ ε 3 (y) = 0}.
The set C is compact in ω. Since the transverse component of ζ ε , solution to problem P ε (ω), belongs to the space H 2 0 (ω) and since θ ε > 0 in ω, it follows that C ∩ γ = ∅. For any ρ > 0, define the set The set C ρ is compact and such that, for sufficiently small ρ, C ρ ∩ γ = ∅. Moreover, we can choose ρ sufficiently small so that C 2ρ ∩ γ = ∅.

Lemma 9 There exist positive numbers h 0 and β 1 such that

θ ε (y) + ζε 3,h (y) ≥ β 1 if y ∈ ω and dist(y, C ) ≥ ρ,
for all h ≤ h 0 .

Proof Since (θ ε + ζ ε 3 ) > 0 outside the contact zone, then it is a fortiori > 0 in the compact set {y ∈ ω; dist(y, C ) ≥ ρ}. By virtue of (51) we immediately infer that

(θ ε + ζε 3,h ) → (θ ε + ζ ε 3 ) in C 0 (ω).
As a result, there exists h 0 > 0 such that

θ ε + ζε 3,h > 0 in {y ∈ ω; dist(y, C ) ≥ ρ},
and there thus exists β 1 > 0 such that

θ ε + ζε 3,h ≥ β 1 in {y ∈ ω; dist(y, C ) ≥ ρ}, for all h ≤ h 0 .
Following the ideas of [START_REF] Brenner | Finite element methods for the displacement obstacle problem of clamped plates[END_REF], we introduce the nodal interpolation operator for the conforming P 1 finite element associated with the triangulation T h and we denote it by I h . By definition of I h , it follows that ζε 3,h and I h ζε 3,h agree at the vertices of the conforming P 1 finite elements. By the linearity of I h we get

I h θ ε + I h ζε 3,h ≥ 0 in ω, (52)
since, again by the properties of I h , the functions I h θ ε and I h ζε 3,h are affine over ω. By standard interpolation estimates (Theorem 3.1.5 of [START_REF] Ciarlet | The Finite Element Method for Elliptic Problems[END_REF] with m = 0, p = 2, q = ∞, and k = 1), we infer

η -I h η 0,∞,ω ≤ Ch|η| 2,ω for all η ∈ H 2 (ω). ( 53
)
An application of ( 53) and (49) yields

ζε 3,h -I h ζε 3,h 0,∞,ω ≤ Ch. (54)
Since θ ε ∈ C 3 (ω), we deduce, by Taylor's theorem with integral remainder that there exists a positive constant C such that sup

y∈ω |θ ε -I h θ ε | ≤ Ch 2 , (55)
and such an estimate a fortiori holds for the norm

• L ∞ (ω) . Define δ h := ( ζε 3,h -I h ζε 3,h ) + (θ ε -I h θ ε ) 0,∞,ω .
In view of (54) and (55), it is straightforward to verify that there exists a positive constant C such that δ h ≤ Ch.

The proof of the next result is obtained by Lemmas 7-9.

Lemma 10 There exists a positive constant C such that

|ζ ε -ζε h | V (ω) ≤ Ch. (57)
Proof Let h 0 and β 1 be as in Lemma 9 and let us assume, without loss of generality, that h ≤ h 0 . By the property (52) of the nodal interpolation operator I h we have

θ ε + ζε 3,h ≥ θ ε + ζε 3,h -I h θ ε -I h ζε 3,h ≥ -δ h in C ρ .
By virtue of (56) it is also licit to assume δ h < β. Let f be an affine function defined over ω as follows

f = 1 in C ρ , (58)
f = 0 in ω \ C 2ρ . (59)
Let denote a mollifier whose support is a subset of C 2ρ . Define the function ϕ ∈ D(ω) by ϕ := * f.

It follows that 0 ≤ ϕ ≤ 1 in ω, (60)
ϕ = 1 in C ρ , (61)
ϕ = 0 in ω \ C 2ρ . (62)
We claim that the function ζε 3,h := ζε 3,h + δ h ϕ belongs to the set K ε 3 (ω). It is straightforward to verify that ζε 3,h = ∂ ν ζε 3,h = 0 on γ. It thus remains to show that (θ ε + ζε 3,h ) ≥ 0 in ω. To this aim, we will distinguish three cases:

Case 1: x ∈ ω \ C 2ρ . In this case, by virtue of (62), we get ζε 3,h = ζε 3,h and the conclusion immediately follows by Lemma 9.

Case 2: x ∈ C 2ρ \ C ρ . In this case, by virtue of (60) and Lemma 9, we get

θ ε + ζε 3,h = θ ε + ζε 3,h + δ h ϕ ≥ β 1 > 0.
Case 3: x ∈ C ρ . In this case, by virtue of (61), we get

θ ε + ζε 3,h = θ ε + ζε 3,h + δ h ≥ 0.
In conclusion, we have shown that ζε 3,h belongs to the set K ε 3 (ω). An application of (49) gives

J ( ζε 1,h , ζε 2,h , ζε 3,h ) = 1 2 b ζε h + (0, 0, δ h ϕ), ζε h + (0, 0, δ h ϕ) - ζε h + (0, 0, δ h ϕ) = 1 2 b( ζε h , ζε h ) -( ζε h ) + b ζε h , (0, 0, δ h ϕ) + 1 2 b ((0, 0, δ h ϕ), (0, 0, δ h ϕ)) -((0, 0, δ h ϕ)) = J( ζε h ) + b ζε h , (0, 0, δ h ϕ) -((0, 0, δ h ϕ)) + 1 2 b ((0, 0, δ h ϕ), (0, 0, δ h ϕ)) ≤ J( ζε h ) + Cδ h .
By the V (ω)-ellipticity of b(•, •) (cf. Theorem 3.6-1 of [START_REF] Ciarlet | Mathematical Elasticity[END_REF]), the intermediary inequality (48) and the fact that ζε 3,h is in K ε 3 (ω) (see Lemma 10), we obtain

α 2 ζ ε -ζε h V (ω) ≤ 1 2 b(ζ ε -ζε h , ζ ε -ζε h ) = 1 2 b(ζ ε , ζ ε ) -b( ζε h , ζ ε -ζε h ) - 1 2 b( ζε h , ζε h ) ≤ 1 2 b(ζ ε , ζ ε ) -(ζ ε -ζε h ) - 1 2 b( ζε h , ζε h ) = J(ζ ε ) -J( ζε h ) ≤ J ( ζε 1,h , ζε 2,h , ζε 3,h ) -J( ζε h ) ≤ Cδ h .
The conclusion immediately follows by (56).

Convergence analysis

The next lemma, inspired by Lemma 4.2 of [START_REF] Brenner | Finite element methods for the displacement obstacle problem of clamped plates[END_REF], provides an estimate for the term

[b h (ζ ε , Π h ζ ε -ζ ε h ) -(Π h ζ ε -ζ ε h )],
which thus allows us to complete the error analysis. The proof relies on Lemma 5, Lemma 10 and standard interpolation estimates (see, e.g., [START_REF] Ciarlet | The Finite Element Method for Elliptic Problems[END_REF]).

Lemma 11 There exists a positive constant C such that

b h (ζ ε , Π h ζ ε -ζ ε h ) -(Π h ζ ε -ζ ε h ) ≤ C √ h √ h + Π h ζ ε -ζ ε h . (63)
Proof Observe that we can write

b h (ζ ε , Π h ζ ε -ζ ε h ) = b h (ζ ε , Ẽh Π h ζ ε -Ẽh ζ ε h ) + b h (ζ ε , Π h ζ ε -ζ ε h -Ẽh (Π h ζ ε -ζ ε h )) = b(ζ ε , Ẽh (Π h ζ ε -ζ ε h )) + b h (ζ ε , Π h ζ ε -ζ ε h -Ẽh (Π h ζ ε -ζ ε h )) ≤ b(ζ ε , Ẽh (Π h ζ ε -ζ ε h )) + Ch ζ ε ω Π h ζ ε -ζ ε h ,
where the latter inequality holds by [START_REF] Rodríguez-Arós | Mathematical justification of the obstacle problem for elastic elliptic membrane shells[END_REF]. We have thus shown that there exists a constant C > 0 such that

b h (ζ ε , Π h ζ ε -ζ ε h ) ≤ b(ζ ε , Ẽh (Π h ζ ε -ζ ε h )) + Ch ζ ε ω Π h ζ ε -ζ ε h . (64)
Let us now estimate the term b(ζ ε , Ẽh (Π h ζ εζ ε h )). We first notice that we can write it in the more suitable equivalent form

b(ζ ε , Ẽh (Π h ζ ε -ζ ε h )) = b( ζε h , Ẽh (Π h ζ ε -ζ ε h )) + b(ζ ε -ζε h , Ẽh (Π h ζ ε -ζ ε h )). ( 65
)
Using ( 38), (57), the continuity of b(•, •) and Poincaré-Friedrichs inequality (Theorems 6.5-2 and 6.8-1 of [START_REF] Ciarlet | Linear and Nonlinear Functional Analysis with Applications[END_REF]), we can estimate the second term in the right hand side of (65) as follows

b(ζ ε -ζε h , Ẽh (Π h ζ ε -ζ ε h )) ≤ C|ζ ε -ζε h | V (ω) | Ẽh (Π h ζ ε -ζ ε h )| V (ω) ≤ C √ h Π h ζ ε -ζ ε h .
As a result, we obtain

b(ζ ε -ζε h , Ẽh (Π h ζ ε -ζ ε h )) ≤ C √ h Π h ζ ε -ζ ε h . (66)
Regarding the first term in the right hand side of (65), we observe that (48) yields

b( ζε h , Ẽh (Π h ζ ε -ζ ε h )) = b( ζε h , ζε h -Ẽh ζ ε h ) + b( ζε h , Ẽh Π h ζ ε -ζ ε h ) ≤ ( ζε h -Ẽh ζ ε h ) + b( ζε h , Ẽh Π h ζ ε -ζε h ). (67)
We note that

b( ζε h , Ẽh Π h ζ ε -ζε h ) = b( ζε h -ζ ε , Ẽh Π h ζ ε -ζε h ) + b(ζ ε , Ẽh Π h ζ ε -ζ ε ) + b(ζ ε , ζ ε -ζε h ). ( 68
)
We estimate the sum of the first two terms of the right hand side of (68) as follows |b( ζε

h -ζ ε , Ẽh Π h ζ ε -ζε h ) + b(ζ ε , Ẽh Π h ζ ε -ζ ε )| ≤ C|ζ ε -ζε h | V (ω) |( ζε h -ζ ε ) + (ζ ε -Ẽh Π h ζ ε )| V (ω) + C ζ ε ω |ζ ε -Ẽh Π h ζ ε | V (ω) ≤ C |ζ ε -ζε h | 2 V (ω) + |ζ ε -ζε h | V (ω) |ζ ε -Ẽh Π h ζ ε | V (ω) + ζ ε ω |ζ ε -Ẽh Π h ζ ε | V (ω) ≤ Ch, (69)
where the latter inequality is obtained by virtue of (36) (for the transverse component only), Poincaré-Friedrichs inequality (Theorems 6.5-2 and 6.8-1 of [START_REF] Ciarlet | Linear and Nonlinear Functional Analysis with Applications[END_REF]) and standard interpolation estimates (Theorem 3.1.5 of [START_REF] Ciarlet | The Finite Element Method for Elliptic Problems[END_REF] with m = k = 1 and p = q = 2). Let us assume, without loss of generality, that h is sufficiently small so that the definition of ζε 3,h is justified (see Lemma 10).

An application of ( 10), [START_REF] Li | Finite element method and its applications[END_REF] (for the transverse component only) and standard interpolation estimates (Theorem 3.1.5 of [START_REF] Ciarlet | The Finite Element Method for Elliptic Problems[END_REF] with m = 0, k = 1 and p = q = 2) yields In conclusion, we have shown that there exists a constant C > 0 such that

b(ζ ε , ζ ε -ζε h ) ≤ ( Ẽh Π h ζ ε -ζε h ) + Ch. ( 70
)
An application of ( 64)-( 70), Hölder's inequality and (38) yields

b h (ζ ε , Π h ζ ε -ζ ε h ) -(Π h ζ ε -ζ ε h ) ≤ b ζ ε , Ẽh (Π h ζ ε -ζ ε h ) + Ch ζ ε ω Π h ζ ε -ζ ε h -(Π h ζ ε -ζ ε h ) = b ζε h , Ẽh (Π h ζ ε -ζ ε h ) + b ζ ε -ζε h , Ẽh (Π h ζ ε -ζ ε h ) + Ch Π h ζ ε -ζ ε h -(Π h ζ ε -ζ ε h ) ≤ b ζε h , Ẽh (Π h ζ ε -ζ ε h ) + C √ h Π h ζ ε -ζ ε h -(Π h ζ ε -ζ ε h ) ≤ ( ζε h -Ẽh ζ ε h ) + b ζε h , Ẽh Π h ζ ε -ζε h + C √ h Π h ζ ε -ζ ε h -(Π h ζ ε -ζ ε h ) = b ζε h -ζ ε , Ẽh Π h ζ ε -ζε h + b ζ ε , Ẽh Π h ζ ε -ζ ε + b(ζ ε , ζ ε -ζε h ) + ( ζε h -Ẽh ζ ε h ) + C √ h Π h ζ ε -ζ ε h -(Π h ζ ε -ζ ε h ) ≤ Ch + C √ h Π h ζ ε -ζ ε h + b(ζ ε , ζ ε -ζε h ) + ( ζε h -Ẽh ζ ε h ) -(Π h ζ ε -ζ ε h ) ≤ Ch + C √ h Π h ζ ε -ζ ε h -(Π h ζ ε -ζ ε h ) + ( Ẽh (Π h ζ ε -ζ ε h )) ≤ Ch + C √ h Π h ζ ε -ζ ε h .
To sum up, we have shown that there exists C > 0 such that

b h (ζ ε , Π h ζ ε -ζ ε h ) -(Π h ζ ε -ζ ε h ) ≤ C √ h( √ h + Π h ζ ε -ζ ε h ),
which completes the proof.

We are now in a position to recover the error estimate in terms of the norm • , whose definition is recalled here below: η h := η 1,h 1,ω + η 2,h 1,ω + η 3,h h for all η h ∈ Ṽh .

The proof of the error estimate, which constitutes the main result of this paper, resorts to Lemma 6, Lemma 11 and Young's inequality (cf. [START_REF] Young | On Classes of Summable Functions and their Fourier Series[END_REF]).

Theorem 1 There exists a positive constant C such that

ζ ε -ζ ε h ≤ C √ h. (71)
Proof An application of Lemma 6, Lemma 11, (11) and Young's inequality yields

ζ ε -ζ ε h 2 ≤ C 1 Π h ζ ε -ζ ε 2 + C 2 [b(ζ ε , Π h ζ ε -ζ ε h ) -(Π h ζ ε -ζ ε h )] ≤ C 1 Π h ζ ε -ζ ε 2 + C √ h( √ h + Π h ζ ε -ζ ε h ) ≤ Ch + C √ h( √ h + Π h ζ ε -ζ ε + ζ ε -ζ ε h ) ≤ C h + √ h ζ ε -ζ ε h ≤ C h + Ch 2 + 1 2C ζ ε -ζ ε h 2 ≤ Ch + 1 2 ζ ε -ζ ε h 2 .
In conclusion, we obtain

ζ ε -ζ ε h 2 ≤ Ch,
and (71) is thus proved.

e∂

  ν e (Π 3,h ξ 3 ) ds = e ∂ ν e ξ 3 ds for all e ∈ E h ,

Fig. 1 :

 1 Fig. 1: Morley's triangle.Figure 6.2.3 of [14] Fig. 2: HCT macro-element.Figure 6.1.3 of [14]

  Fig. 1: Morley's triangle.Figure 6.2.3 of [14] Fig. 2: HCT macro-element.Figure 6.1.3 of [14]

  Fig. 1: Morley's triangle.Figure 6.2.3 of [14] Fig. 2: HCT macro-element.Figure 6.1.3 of [14]

Fig. 3 :

 3 Fig. 3: Configuration associated with the jump of the normal derivative across the edge e.

  vanishes at the endpoints of the edge e. As a result, by the mean value theorem and the properties of quadratic polynomials, we deduce that the tangential derivative along the edge e ∂ τ e (η 3 | T + -η 3 | T -) vanishes at the midpoint of e. By virtue of the decomposition of the gradient in terms of tangential and normal derivatives we get the continuity of the gradient at the midpoint of any edge e ∈ E h . The other property follows from the boundary conditions.

C

  ρ := {y ∈ ω; dist(y, C ) ≤ ρ}, where dist(y, C ) denotes the distance of any point y ∈ ω from the set C , i.e., dist(y, C ) := min x∈C |y -x|.

  b(ζ ε , ζ εζε h ) = b ζ ε , ζ ε -( ζε 1,h , ζε 2,h , ζε 3,h ) + b (ζ ε , (0, 0, δ h ϕ)) ≤ ζ ε -( ζε 1,h , ζε 2,h , ζε 3,h ) + δ h b (ζ ε , (0, 0, ϕ)) = ( Ẽh Π h ζ εζε h ) + (ζ ε -Ẽh Π h ζ ε ) -δ h [ ((0, 0, ϕ)) -b (ζ ε , (0, 0, ϕ))] ≤ ( Ẽh Π h ζ εζε h ) + Ch ζ ε ω -δ h [ ((0, 0, ϕ)) -b (ζ ε , (0, 0, ϕ))] ≤ ( Ẽh Π h ζ εζε h ) + Ch.

Acknowledgements

This paper was started while the first author was visiting the Xi'an University of Technology, as a result of an invitation from Professor Xiaoqin Shen (second author). The first author warmly thanks Professor Xiaoqin Shen and the Department of Applied Mathematics of the Xi'an University of Technology for their hospitality.

The authors are greatly indebted to Professor Philippe G. Ciarlet for his encouragement and guidance.

This paper is supported by the ERC advanced grant 668998 (OCLOC) under the EU's H2020 research program, the National Natural Science Foundation of China (NSFC. 11971379, 11571275, 11572244), the Natural Science Foundation of Shaanxi Province ( 2018JM1014).


Céline Fiot

Anne Laurent

Maguelonne Teisseire

Approximate Sequential Pattern for Incomplete Sequence Database Mining

Industrial databases often contains a large amount of unfilled information. When these data are mined for frequent sequences, incomplete data are, most of the time, deleted, which leads to an important loss of information. Extracted knowledge then becomes less representative of the database. Two techniques can then be investigated: either using only the available information or estimating the missing values. In this paper we propose an estimation-based approach that represents inclusion of an item within a record by a fuzzy set. Then the membership degree giving the item inclusion is used to compute the frequency of a sequence. Experiments run on various synthetic datasets show the feasibility and validity of our proposal as well in terms of quality as in terms of the robustness to the rate of missing values.

I. INTRODUCTION

For the last decade, data mining application fields have widened. Especially it is the case of sequential pattern mining approaches [START_REF] Agrawal | Mining sequential patterns[END_REF]. This data mining technique aims at discovering knowlegde from temporal databases. First designed to analyse customer behavior, they are now used in various industrial, medical, or biological fields. The searched datasets are sets of time-stamped records. Each of them is constituted of a set of values.

However these kind of data often contain imperfectionssuch as noise or incomplete data due to breakdowns or errors, for instance. Moreover, existing sequential patterns methods only allow the analysis of complete data, without considering incomplete records, which may represent an important loss of information. Besides, replacement of missing values or their estimation are often either too simplistic to perform unbiased results, or too time-consuming to be applied on large datasets. Thus, it seems necessary to implement a method for mining sequential patterns within incomplete databases. To solve this problem two tracks can be investigated. The concepts linked to sequential patterns can be adapted either to ignoremissing informationonly using the available one or to consider likely value of incomplete data.

In this paper we propose an estimation-based approach to discover frequent sequences within databases containing random missing values. This approach is instigated by an association rule mining method within incomplete databases [START_REF] Nayak | Approximate association rule mining[END_REF], [START_REF] Ng | Quantitative association rules over incomplete data[END_REF] and by a technique often used in machine learning or statistics: assessing missing value with respect to the values in the database. This principle consists in using available information, with a certain likelihood. Thus the whole dataset Céline Fiot , Anne Laurent and Maguelonne Teisseire are with the LIRMM, University of Montpellier II -CNRS, 34392 Montpellier, France (email: {fiot, laurent, teisseire}@lirmm.fr).

is used to discover the patterns, while taking into account uncertainty due to incomplete data instead of ignoring it, deleting it or simply replacing it by a single value. We here redefine this principle for sequential pattern mining, thus proposing to discover approximate maximal frequent sequences within incomplete time-stamped databases. For this reason, we adapted the concepts linked to sequential pattern discovery. Then we implemented the ApSPoID algorithm using a sequential pattern mining algorithm, PSP [START_REF] Masseglia | The PSP Approach for Mining Sequential Patterns[END_REF]. This algorithm was tested on synthetic datasets to show the validity of our approach.

The remainder of the paper is organized as follows. Section II introduces the methods for association rule mining in the presence of missing values, and the concepts linked to sequential pattern discovery. Section III details our approach to mine for approximate sequential patterns within incomplete data, and we describe our algorithm in Section IV. Section V is then dedicated to experiments that show the validity of our approach. Finally, we conclude in section VI with a discussion on further work opened by this proposal.

II. FROM SEQUENTIAL PATTERNS TO INCOMPLETENESS

Sequential patterns are often introduced as an extension of association rules, initially proposed in [START_REF] Agrawal | Mining association rules between sets of items in large databases[END_REF]. They highlight correlations between database records as well as their temporal relationship. Even so these algorithms do not mine incomplete records contained in the database. These missing values must then be removed either by deletion or replacement. Quality of results then depends on this preprocess. Moreover this step is often time-consuming. In order to reduce the preprocessing due to missing values and to improve the sequential pattern quality, we propose a method for sequential pattern mining within incomplete databases. This method is based on association rules approaches. In this section we first define the concepts linked to sequential pattern mining, then we detail our motivations before introducing techniques that allow association rule mining handling missing values.

A. Sequential Patterns

Sequential patterns are based on the idea of maximal frequent sequences. Let R be a set of objects records where each record R consists of three information elements: an object-id, a record timestamp and a set of attributes/items in the record. Let I = {i 1 , i 2 , ..., i m } be a set of items or attributes. An itemset is a non-empty set of attributes i k , denoted by (i 1 i 2 . . . i k ). It is a non-ordered representation. A sequence s is a non-empty ordered list of itemsets s p , denoted by < s 1 s 2 ...s p >.

A n-sequence is a sequence of n items (or of size n).

Example 1: Let us consider an example of market basket analysis. The object is a customer, and records are the transactions made by this customer. Timestamps are the date of transactions. If a customer purchases products e, a, k, u, and f according to the sequence s =< (e) (a k) (u) (f ) >, then all items of the sequence were bought separately, except products a and k which were purchased at the same time. In this example, s is a 5-sequence.

One sequence S =< s 1 s 2 ...s p > is a subsequence of another one

S ′ =< s ′ 1 s ′ 2 ...s ′ m > if there are integers l 1 < l 2 < ... < l p such that s 1 ⊆ s ′ l1 , s 2 ⊆ s ′ l2 , ..., s p ⊆ s ′ lp . Example 2: The sequence s ′ = <(a) (f)> is a subsequence of s because (a) ⊆ (a k) and (f) ⊆ (f). However, <(a) (k)> is not a subsequence of s.
All records from the same object o are grouped together and sorted in increasing order of their timestamp. They are called a data sequence. An object supports a sequence s if it is included within the data sequence of this object (s is a subsequence of the data sequence). The frequency of a sequence (f req(s)) is defined as the percentage of objects supporting s in the whole set of objects O. In order to decide whether a sequence is frequent or not, a minimum frequency value (minF req) is specified by the user and the sequence is said to be frequent if the condition f req(s) ≥ minF req holds. A sequence that may be frequent is a candidate sequence. Given a database of object records, the problem of sequential pattern mining is to find all maximal sequences of which the frequency is greater than a specified threshold (minF req) [START_REF] Agrawal | Mining sequential patterns[END_REF]. Each of these sequences represents a sequential pattern, also called a maximal frequent sequence.

Several extensions were proposed to consider incremental mining for sequential patterns [START_REF] Masseglia | Incremental mining of sequential patterns in large databases[END_REF], to handle numerical and quantitative values [START_REF] Hong | Mining Fuzzy Sequential Patterns from Multiple-Items Transactions[END_REF], [START_REF] Chen | Discovery of Fuzzy Sequential Patterns for Fuzzy Partitions in Quantitative Attributes[END_REF], [START_REF] Fiot | From crispness to fuzziness: Three algorithms for soft sequential pattern mining[END_REF] or to generalize sequential patterns with respect to various temporal parameters (timeinterval between events of a sequence, grouping several records into a single itemset...) [START_REF] Srikant | Mining sequential patterns: Generalizations and performance improvements[END_REF], [START_REF] Masseglia | Pre-Processing Time Constraints for Efficiently Mining Generalized Sequential Patterns[END_REF], [START_REF] Fiot | Des motifs séquentiels généralisés aux contraintes de temps étendues[END_REF]. However, no technique was proposed to deal with missing values while sequential pattern mining. For this reason, in the following sections, we propose an approach that can mine maximal frequent sequences from an incomplete sequence database.

B. Motivations

Let us consider the database given by Table I. The goal is to extract the sequential patterns with a minimum frequency equal to 50%. Now, let us consider the incomplete database given by Table II. Some information in the data sequences has not been filled in and these values are missing. Let us consider that these values are identified as missing and unfilled. In order to mine for patterns, previous approaches require the suppression of missing values.

Obj

. Sequence O1 (a b) (? ? c) (? b c) O2 (a) (? c) (b d) O3 (a b) (? c) (? ? c)
After preprocessing, the database within which sequential patterns are mined is given Table III.

Obj.

Sequence

O1 (a b) (c) (b c) O2 (a) (c) (b d) O3 (a b) (c) (c)
Resulting patterns for minF req = 50% are:

<(a b)(c)(c)> and <(a)(c)(b)>.
Note that only a short part of the database is used to discover information and only part of frequent patterns and frequent items of the complete dataset are found. In particular, item d is not frequent in the incomplete database. Therefore using the whole dataset without deleting part of it appears essential.

C. Association Rules and Missing Values

Some works were proposed to discover association rules within incomplete databases. The algorithm RAR proposed by [START_REF] Ragel | Treatment of missing values for association rules[END_REF] allow the user to consider incomplete data while association rule mining within incomplete relational databases, thanks to partial and temporary omission of such incomplete records. The whole database is not used to discover each rule but to generate the whole set of rules. This technique adapted to sequence mining in [START_REF] Fiot | Spoid: Extraction de motifs séquentiels pour les bases de données incomplètes[END_REF] gives interesting results. However, this approach requires the specification of a minimum significancy threshold and experiments show that if this threshold is not the optimal value, the method does not discover all the frequent patterns or too many parasitic patterns.

Therefore we develop a new approach, instigated by the proposals of [START_REF] Nayak | Approximate association rule mining[END_REF], [START_REF] Ng | Quantitative association rules over incomplete data[END_REF]. These methods implement an assessment systems based on a probabilistic distribution. With this approach one missing value can represent various values while mining for association rules. In order to consider those assessments, the notion of support (percentage of records in database that include the rule items) and thus confidence (probability for a record to contain the right part of the rule knowing it contains the left part) were reformulated.

From this principle we redefine in the next section the definitions linked to sequential pattern mining previously detailed.

III. APSPOID: DEALING WITH INCOMPLETENESS

In this section, we present the definitions we used for our algorithm to mine incomplete sequence databases for sequential patterns.

A. Approximate Sequential Patterns for Incomplete Databases

As deleting incomplete records leads to an important loss of information, we adapted an association rule mining approach, robust to missing values. In this section, we describe our method ApSPoID, based on the principles of the ∼AR algorithm [START_REF] Nayak | Approximate association rule mining[END_REF].

The main idea of our approach, as the one of the ∼AR method lies on the estimation of the value for unfilled items. While ∼AR computes the frequency of an incomplete itemset taking into account the probability of each likely value of a missing value, we propose here to use this estimation for the frequency computation of sequences.

Thus, while computing the item frequency of a sequence, for each missing value, several likely value may be checked. Each one of these possibilities is associated to a degree giving the likely appearance level. The support definition given in Section II is adapted in order to include this inclusion degree of the item.

Then inclusion of an item in an itemset is weighted by the certainty degree of its presence:

• If the item is in the itemset, the appearance level is equal to 1; • If the item is not in the itemset and this itemset does not contain missing value, the appearance level is nul; • Else, the appearance level of item i in the record r is in ]0, 1[. It will be denoted by µ r (i). Definition 1: The frequency of an itemset I is the average appearance level of the itemset for all the objects in the dataset. For each of them, the best occurrence is kept (i.e. the one being the most certain). The frequency is then computed by the formula (1) :

F req(I) = o∈O ⊥ r∈Ro ⊤ i∈I µ r (i) |O| (1)
Definition 2: Le frequency of a sequence S is the average appearance rate of the sequence for all the objects in the database. It is defined as the ratio of the sum of the best appearance level of the sequence for all the data sequences over the total number of data sequences. It is given by the formula (2):

F req(S) = o∈O F(o, S) |O| (2)
where L(S, o) is the appearance level of S in the data sequence of the object o. Definition 3: The appearance level of a sequence S for an object o corresponds to the best inclusion rate of the sequence S in the data sequence of o. This appearance level is computed by the ordered agregation of the itemset inclusion rate in the order of appearance in sequence S. It is given by the formula (3):

L(S, o) = ⊥ r∈Ro OrdAg <s1...si...s k > (⊤ j∈si µ r (j)) (3)
where k is the number of itemset in S.

B. Choosing Operators

Several operators could be used to agregate the item inclusion rates of an itemset or the itemset inclusion rate for a sequence. The method proposed by [START_REF] Nayak | Approximate association rule mining[END_REF] uses the average agregation for item inclusion rates. However we consider that this type of agregation can hide a low appearance level of several items in a record if the other ones are certain. Therefore we propose to use a t-norm operator for this computation. More especially we will implement our approach with the min operator for its idempotency property [START_REF] Fiot | From crispness to fuzziness: Three algorithms for soft sequential pattern mining[END_REF].

Regarding the ordered agregation, we propose to use the average agregation, each itemset having the same "interest" in the sequence. Uncertainty is indeed already taken into account by the use of the t-norm. It is then not necessary to consider it a second time at the sequence level [START_REF] Fiot | From crispness to fuzziness: Three algorithms for soft sequential pattern mining[END_REF].

C. Computing the Appearance Level

First step of our method consits in determining each possible value for missing values and the certainty degree of each one.

An aproximation of incomplete record composition should be determined. To do so, the dataset is scanned to compute appearance rate of each item. As the database is incomplete, the appearance rate is computed as the frequency with respect to the number of records using the disabling technique of RAR [START_REF] Ragel | Treatment of missing values for association rules[END_REF]. This computation results in the frequency distribution that will give the appearance level of an item in a record. This value distribution is here computed with respect to the appearance within a record, but it could be computed with respect to the appearance within data sequences or only considering complete records.

Thus, parsing the database for a candidate sequence, when an itemset is found complete in a record, its appearance level is 1. When it is partially found in an incomplete record, the frequency distribution is used to determine the appearance level of the missing part of the itemset.

We present some experiment results that show the consequences of choosing one distribution or another on discovered patterns.

IV. IMPLEMENTATION

In this section, we first apply our defintions and principles on a brief example. Then, we present how our algorithm runs.

A. Example

Consider the incomplete database given by Table II. minSup is equal to 50%.

First, the appearance rate distribution is computed. We use the formula given by [START_REF] Ragel | Treatment of missing values for association rules[END_REF] :

F req(i) = |{r ∈ R/i ∈ r}| |R| -|{r ∈ R/r.

isIncomplete(i)}|

For i = a : 3/(9-5) = 0.75. Table IV gives the frequency distribution for each value that could replace a missing one. Thus for the second record of the first object, the appearance level of c is 1, but a, b, d and e, it is given by the frequency distribution.

Then sequential patterns mining starts with the frequency computation for each item. Item a certainly appears in each data sequence then its frequency is 3/3=100%, the same for items b and c. The item d, could be supported by objects 1 and 3, iots frequency is then given by the following calculation: F req(< d >) = ⊥(0.25,0.25)+⊥(0.25,1)+⊥(0.25,0.25)

3 = 0.25+1+0.25 3 = 0.5
Now, consider sequence S = < (ab) (abcd) >. This candidate sequence cannot be supported by none of the data sequences because none of them includes an itemset composed of 4 items. Consider sequence S ′ = < (ab) (bc) >, it is supported by object 1, cannot be supported by object 2but could be supported by the third object. For the frequency of S ′ , we get:

F req(S ′ ) = 1 + 0 + (1 + ⊤(0.8, 1))/2 3 = 1 + 0.9 3 = 63%
Sequence S ′ is then frequent. Applying our method, discovered sequential patterns for minSup=50% are :

< (ab)(bc)(abc) >, < (ab)(ac)(abc) >, < (a)(bc)(bd) >, < (a)(ac)(bd) >, < (b)(c)(bd) > and < (a)(d)(b) >.
Even if these patterns are not exactly the one obtained on the complete database, they are closer to the one we should get than the one discovered using the preprocessed dataset. Experiments detailed in section V show that there exists a value of the minimum representativity for which the algorithm SPoID extracts the whole set of sequential patterns of the complete database from an incomplete one.

B. Algorithm

The algorithm SPoID runs similarly to the generate-prune sequential pattern mining algorithms. It consists in generating all the candidate k-sequences from the frequent (k-1)sequences. Then the database is scanned to count the number of data sequences that support each candidate sequence. The main difference stands in the counting of incomplete data sequences. This counting step uses the same counting principles used by the algorithm TotallyFuzzy for fuzzy sequential pattern mining [START_REF] Fiot | From crispness to fuzziness: Three algorithms for soft sequential pattern mining[END_REF]. It is described by Alg. 1. If the first itemset is found complete, the appearance level for the first itemset is 1, then the following itemset is searched. If the first itemset is partially found in an incomplete record, then appearance level for this itemset I is min i∈I µ r (i), then the following part of the sequence is searched. On the other hand, if the record only partially supports the itemset and if the number of missing value des not allow to complete it, the first itemset is not found. Then it is searched in the remainder of the data sequence. If the sequence is found, its appearance level is computed by an average agregation. If it is found several times, the best appearance level is kept. Once the whole database is parsed, the sum of appearance level is divided by the number of objects in database. Then all the candidate sequences having a frequency lower than minF req are pruned.

ApSPoID -

The temporal complexity of this algorithm is the same as the one of the algorithm TotallyFuzzy presented by [START_REF] Fiot | From crispness to fuzziness: Three algorithms for soft sequential pattern mining[END_REF]. We use the same kind of optimizations to reduce the number of database scans. The path structure described for TotallyFuzzy keeps in memory the most certain occurrence of a sequence and each possible start while scanning the database.

V. EXPERIMENTS

These experiments were carried out on a PC -Linux 2.6.7 OS, CPU 2.8 GHz with 512 MB of memory. The algorithm was implemented in Java on the PSP principle. In particular, the Prefix-Tree structure was used to store the candidate and frequent sequences. We used synthetic datasets randomly generated by a normal distribution. Then some items were randomly replaced by missing values. Sequential patterns were extracted from the complete database and from the preprocessed incomplete ones (i.e. incomplete databases in which incomplete records have been deleted). Then those patterns are compared to the one discovered by our algorithm SPoID. Results here detailed were obtained from several synthetic datasets containing around 2000 sequences of 20 transactions in average. Each transaction contains around 10 items chosen among 100.

Our analysis is based on the several counting:

• the total number of sequential patterns discovered by SPoID, • the number of sequential patterns discovered by SPoID, that are discovered in the complete database, • the number of wrong sequential patterns discovered by SPoID (that groups together the patterns that are not discovered in the complete database and the one not found by SPoID but should be). Table V sums up these notations.

A. ApSPoID vs. SPoID

Following figures give a comparison between patterns discovered by ApSPoID and those discovered by PSP after pre-processing (deletion of incomplete data) and also those discovered by SPoID ( [START_REF] Fiot | Spoid: Extraction de motifs séquentiels pour les bases de données incomplètes[END_REF], an existing approach for mining incomplete datasets. First, Figure 1 shows the evolution of the ratio β/τ . Ap-SPoID and SPoID have similar behaviors regarding the number of right discovered patterns. More especially, ApSPoID gives better results fronm 40% to 70% of incompleteness. On this interval, the results of SPoID are less interesting than those obtained on the range 10 to 20%, where SPoID discovers all the patterns of te complete database. It can also be noted that wathever the percentage of incompleteness, the pre-processing method is not powerful. Then figures 2(a) and 2(b) show a comparison of parasitic and missing patterns. Figure 2(a) presents the proportion of right patterns with respect to discovered patterns , according to the incompleteness percentage. This ratio is considered as good when it is close to 1 (each discovered patterns is also discovered in the complete database). However, these results should be compared to Figure 1. All the extacted patterns may indeed correspond to the complete database but an important amout of them may be missing. In our experiments, it is the case of PSP with pre-processing.

Note that ApSPoID results are better than those of SPoID from 40% of incompleteness. At this rate ApSPoID indeed discovers 90% of the patterns found in the complete database and those patterns represent around 60% of discovered patterns. Furthermore, when the database contains 50% of missing values, ApSPoID discovers 80% of the squential patterns of the complete database and those patterns represent 70% of discovered patterns while SPoID only find 70% of the complete database patetrns, only representing 40% of discovered patterns.

These observations are confirmed by Figure 2(b) that shows the evolution of β/delta. It shows that form 35% of incompleteness, the proportion of right patterns according to the different ones is clearly more important for ApSPoID than for SPoID. Moreover, the number of right pattern is still higher than the part of different patterns (this ratio should remain greater than 1), on the contrary to SPoID for which results are good until 30% of incomlpeteness and then collapse.

B. Choosing a distribution

We analyse the behavior of ApSPoID according to the frequency distribution for incomplete records. We show that it has an impact on obtained results.

We here present a comparison of ApSPoID performances used with two different distributions:

• One distribution is computed using the RAR algorithm (as in example, Section IV). It is used by ApSPoID -Dist 1 on Fig. 4. • A second distribution is found by considering the whole set of records to compute the frequency. It is used by ApSPoID -Dist 2 on Fig. 4. While the first distribution is rather optimistic (each missing value is ignored for computing the distribution), the second one is more pessimistic (each missing value is considered as different from the item of which the frequency is computed). Choosing one or the other has an impact on the final results. Fig. 3 shows that the second distribution does not allow ApSPoID to find the patterns of the complete database from 40% of incompleteness, whereas ApSPoID with the optimistic distribution finds an important amount of patterns until 60% of incompleteness. On the other hand, Fig. 4(a) shows that ApSPoID with the second distribution discovers less parasitic patterns.

From Fig. 4(b) we can conclude that the pessimistic distribution gives better results than the optimistic one until 40% of incompleteness. On the range 20 to 40% of incomplete data, results obtained by ApSPoID with the second distribution are similar to those of SPoID. Thus, choosing a distirbution more or less optimistic determines which quantity of sequential patterns found in the complete database will be retrieved by ApSPoID, but also the amount of parasitic patterns. A trade-off should then be made between the percentage of right patterns found and the percentage of parasitic patterns.

These experiments show that in general , ApSPoID is more reliable than SPoID from 30% of missing values and that these results are interesting until 60% of incomplete data. On the other hand patterns discovered by ApSPoID -Dist 1 include 80 to 90% of patterns contained in the complete database for an incompleteness rate between 10 to 40%. These results are almost as good as those of SPoID. However, using ApSPoID does not require to specify other parameter than minF req. Its use is then easier than the use of SPoID. Furthermore, results obtained by ApSPoID could be improved by refining the frequency distribution at each iteration, considering more and more precisely the whole records. Lastly, regarding the comparison of runtime, ApSPoID behaves similarly to SPoID.

VI. FURTHER WORK

Sequential pattern discovery is an interesting data mining method while learning knowledge from large time-stamped databases, such as industrial process databases. But this kind of databases often contains missing values. However there is only one mining techniques for discovering frequent sequences from incomplete databases. Therefore in this paper we proposed to adapt the original definitions related to sequential pattern mining in order to handle random incompleteness in data sequences. These new definitions enable the user to manage missing values directly during the mining task, then avoiding an heavy preprocessing step. Our method and algorithm ApSPoID has been implemented and tested on synthetic datasets. We have thus shown the robustness of our approach until an incompleteness rate around 50% , while the existing approaches give erroneous results since 10% of missing values. We now work on extending this approach in

  b) (b c d) (b c e) O2 (a) (b c) (b d) O3 (a b) (b c) (b c d) The sequential patterns obtained are: <(a b)(b c d)>, <(a b)(b c)(b c)> and <(a)(b c)(b d)>.

Fig. 2 .

 2 Fig. 2. (b): β/θ rate according to rate of database incompleteness;(b): β/δ rate according to rate of database incompleteness.

  Ratio #trouves/#dans base complete en fonction de % incompletude ApSPoID -Dist 1 ApSPoID -dist 2

Fig. 3 .

 3 Fig. 3. β/τ rate according to rate of database incompleteness.

Fig. 4 .

 4 Fig. 4. (a): β/θ rate according to rate of database incompleteness;(b): β/δ rate according to rate of database incompleteness.

TABLE III AFTER

 III MISSING VALUE DELETION.

  Input: O, a sequence database ; minFreq, user-specified minimum frequency Ouput: SP List, list of frequent sequences For each candidate sequence, data sequences are parsed for finding the itemsets of the candidate sequence. The subfunction findSequence (Alg. 2) is a recursive function that parses a data sequence t for finding a candidate sequence s. When the searched itemset is found then findSequence searches for the next itemset of s in the data sequence t. Else it keeps on looking for the same itemset in the remainder part of s.

	[Computation of the frequency distribution]
	float[] f reqDist ← calcDist() ;
	[Search for frequent items]
	C ← {i ∈ I} ; k = 1 ;
	F ← getFrepnRep(C,minF req) ;
	SP List.add(F ) ;
	[Search for frequent sequences, size ≥ 2]
	While (C = ∅) do k++ ; C ← generate(F ,k) ;
	For each candidate sequence s ∈ C do
	For each objet o ∈ O do
	float[] degTab; ag ← 0; id ← 1;
	For (j=1 to |Ro| do
	T rans ← {r End For
	tmpF req ← tmpF req + ag ;
	End For
	F req(s) ← tmpF req(s)/|(O)| ;
	If ((F req(s) < minF req)) Then
	prune(s) ;
	End If
	End For
	SP List.add(F ) ;
	End While return SP List ;
	Alg. 1: ApSPoID -Main algorithm.
	findSequence -Input: s, a candidate sequence ; t, a data sequence
		d, table of appearance level for the itemsets of s ;
		n, number of the currently searched itemset
	Ouput: d
	sT ail ← s\s.first() ; tT ail ← t\t.first() ;
	If ((sT ail = ∅)&(tT ail = ∅)) Then If ( min i∈s.first µ s.first (i)) Then
	d[n] ← min i∈s.first	µ s.first (i) ;
	findSequence(sT ail, tT ail, d, n + 1);
	Else
	findSequence(s, tT ail, d, n);
	End If
	Else
	return d ;
	End If

j , ...r |Ro| }; findSequence(s,T rans,degT ab) ; ag ← max(ag, i=id X i=1 degT ab[i]/id) ; Alg. 2: findSequence -data sequence parsing algorithm.

TABLE V NOTATIONS

 V FOR THE DIFFERENT KINDS OF PATTERNS.

	β # sequential patterns discovered by SPoID, also con-
		tained in the complete dataset
	δ	# different sequential patterns
	θ sequential patterns discovered by SPoID in the
		incomplete database
	τ	# sequential patterns discovered in the complete dataset
		Ratio #trouves/#dans base complete en fonction de % incompletude
		#trouves/#dans base complete	0.2 0.4 0.6 0.8 1		SPoID -optimal PSP + pretraitement ApSPoID
			0
			10	20	30	40	50	60	70
					% valeurs manquantes
	Fig. 1. β/τ rate according to rate of database incompleteness.

order to regard other kind of missing values (not randomly spread, for instance), after having detected the different kind of incomplete information. It seems also necessary to detect the attributes that should not be considered as incomplete but not existing. Lastly noise is also a common imperfection in real-life databases. It could then be interesting to consider it in a further version of our algorithm.