Vitoriano Ruas
email: vitoriano.ruas@upmc.fr

Accuracy enhancement for finite-element simulations in curved domains; application to fluid flow

Keywords: Curved domain, Dirichlet, finite elements, mixed, N -simplex, optimal order, Stokes, straight-edged AMS Subject Classification: 65N30, 76M10

Among a few known techniques the isoparametric version of the finite element method for meshes consisting of curved triangles or tetrahedra is the one most widely employed to solve PDEs with essential conditions prescribed on curved boundaries. It allows to recover optimal approximation properties that hold for elements of order greater than one in the energy norm for polytopic domains. However, besides a geometric complexity, this technique requires the manipulation of rational functions and the use of numerical integration. We consider a simple alternative to deal with Dirichlet boundary conditions that bypasses these drawbacks, without eroding qualitative approximation properties. In the present work we first recall the main principle this technique is based upon, by taking as a model the solution of the Poisson equation with quadratic Lagrange finite elements. Then we show that it extends very naturally to viscous incompressible flow problems. Although the technique applies to any higher order velocity-pressure pairing, as an illustration a thorough study thereof is conducted in the framework of the Stokes system solved by the classical Taylor-Hood method.

Introduction

Consider the finite-element solution of second order elliptic equations posed in curved domains with Dirichlet boundary conditions. It is well known that a considerable order lowering may occur if prescribed boundary values are shifted to nodes that are not mesh vertexes of an approximating polygon or polyhedron formed by the union of straight-edged N -simplexes of a fitted mesh. Over four decades ago some techniques were designed in order to remedy such a loss of accuracy, and possibly attain the same theoretical optimal orders as in the case of a polytopic domain, assuming that the solution is sufficiently smooth. Two examples of such attempts in the framework of two-dimensional problems are the interpolated boundary condition method by Nitsche and Scott (cf. [START_REF] Nitsche | The Mathematical Foundations of the Finite Element Method with Applications to Partial Differential Equations[END_REF] and [START_REF] Scott | Finite Element Techniques for Curved Boundaries[END_REF]), and the method introduced by Zlámal in [START_REF]Curved Elements in the Finite Element Method[END_REF] and extended by Žénišek [START_REF] Žénišek | Curved triangular finite C m -elements[END_REF]. Among such techniques the finite element method's isoparametric version is by far the one most widely in use since the sixties (cf. [START_REF] Zienkiewicz | The Finite Element Method in Engineering Science[END_REF]) in order to recover the lost optimality. One of the main reasons why it became so popular is the fact that the isoparametric technique applies to both two-and three-dimensional problems. We recall that this version of the finite element method is based on elements with curved boundary portions, aimed at better approximating a curved boundary than straight edges or plane faces. In this case the aforementioned shift of prescribed boundary values is avoided, since all nodes to which such values apply remain on the true boundary. The price to pay however is the manipulation of rational functions as both shape and trial functions defined upon the curved elements, and the resulting compulsory use of numerical integration. While on the one hand this is far from being an obstacle in most current applications such as linear problems with constant coefficients, numerical integration can be a delicate issue in more complex situations. The technique exploited in this work allows overcoming all such issues, since it is based only on polynomial algebra upon an ordinary (i.e. a straight-edged) N -simplex. Moreover, in contrast to the simple polygonal approach no erosion of the theoretical order of a given interpolation inherent to the method occurs, especially for methods which are not of the lowest possible order. In short, our technique is aimed at ensuring a theoretical order greater or equal to two in the natural norm, without the use of curved elements and interpolating functions other than piecewise polynomials. Actually the conception of the finite-element technique for solving boundary value problems with a smooth curvilinear boundary considered in this work is close to the interpolated (Dirichlet) boundary condition method studied in [START_REF] Brenner | The Mathematical Theory of Finite Element Methods[END_REF]. Though intuitive and known since the seventies, the latter technique has been of limited use so far. Among the reasons for this lies its difficult implementation, the lack of an extension to three-dimensional problems and restrictions on the choice of boundary nodal points to reach optimal convergence rates. In contrast our method is simple to implement in both two-and three-dimensional geometries. Moreover it is particularly handy, whenever a finite element method has normal component or normal derivative degrees of freedom as illustrated in [START_REF] Ruas | A simple alternative for accurate finite-element modeling in curved domains[END_REF]. Indeed when a method incorporates this type of degree of freedom the definition of isoparametric finite-element analogs is not always simple or clear (see e.g. [START_REF] Bertrand | First-order system least-squares on curved boundaries: higherorder Raviart-Thomas elements[END_REF]). It is important to point out that efficient finite-element techniques are known since long, to optimally handle boundary conditions other than Dirichlet's, such as Neumann or Robin boundary conditions prescribed on curved boundaries. In this respect the author refers for instance to the works by Barrett and Elliott [START_REF] Barrett | A Finite-element Method for Solving Elliptic Equations with Neumann Data on a Curved Boundary Using Unfitted Meshes[END_REF] and [START_REF] Barrett | Finite-Element Approximation of Elliptic Equations with a Neumann or Robin Condition on a Curved Boundary[END_REF], besides [START_REF] Strang | An Analysis of the Finite Element Method[END_REF] where a clear explanation on the issues brought about by Neumann conditions prescribed on curved boundaries is given.

The technique applied in this paper was introduced in [START_REF] Ruas | Optimal simplex finite-element approximations of arbitrary order in curved domains circumventing the isoparametric technique[END_REF], in connection with triangular Lagrange finite elements of any order k greater than one to solve the Poisson equation posed in a smooth curvilinear domain. In the subsequent work [START_REF] Ruas | Methods of arbitrary optimal order with tetrahedral finite-element meshes forming polyhedral approximations of curved domains[END_REF] the author addressed the case of tetrahedral Lagrange elements of arbitrary order for second order elliptic PDEs in the same class of domains. A synthesis of both papers is given in [START_REF] Ruas | Variational formulations yielding high-order finite-element solutions in smooth domains without curved elements[END_REF]. In [START_REF] Ruas | A Hermite method for Maxwell's equations[END_REF] the author and co-worker used the same approach to the solution of Maxwell's equations with a Hermite finite element method. In this work we push further such studies in accordance to the following outline. In Section 2 we recall the main results on the new technique provided in [START_REF] Ruas | Optimal simplex finite-element approximations of arbitrary order in curved domains circumventing the isoparametric technique[END_REF] and [START_REF] Ruas | Methods of arbitrary optimal order with tetrahedral finite-element meshes forming polyhedral approximations of curved domains[END_REF], restricted to the case k = 2. Numerical examples are given in Section 3 in connection with the two-dimensional Poisson equation as well. A rigorous study of this technique is carried out in Section 4 in the framework of the finite-element solution of the equations governing incompressible viscous flows. Although such a study applies to different mixed elements or formulations of these equations such as GLS and SUPG, we confine ourselves to the case of the popular Taylor-Hood element (cf. [START_REF] Hood | Navier-Stokes equation using mixed interpolation[END_REF]) as applied to the Stokes system in standard Galerkin formulation. In Section 5 we supply numerical examples illustrating the theoretical results of Section 4, and we conclude in Section 6 with a few comments.

Technique's short description

Referring to [START_REF] Ruas | Optimal simplex finite-element approximations of arbitrary order in curved domains circumventing the isoparametric technique[END_REF] and [START_REF] Ruas | Methods of arbitrary optimal order with tetrahedral finite-element meshes forming polyhedral approximations of curved domains[END_REF] for further details, here we describe our technique to solve boundary value problems with Dirichlet conditions prescribed on smooth curved boundaries, by solving a simple model problem as follows.

Let Ω be a bounded N -dimensional domain for N = 2 or N = 3, and Γ be its boundary assumed to be sufficiently smooth (Γ must be at least of the C 1 -class). Given a function f ∈ H 1 (Ω) we wish to find a function u ∈ H 3 (Ω) that solves -∆u = f in Ω with u = g on Γ assuming that g ∈ H 5/2 (Γ).

Now let P = {T h } h be a uniformly regular family of finite element meshes consisting of straight-edged triangle or tetrahedra, according to the space dimension, satisfying the usual compatibility conditions (see e.g. [START_REF] Ciarlet | The Finite Element Method for Elliptic Problems[END_REF]). Every element of T h is to be viewed as a closed set. Moreover each one of these meshes is assumed to fit Ω in such a way that all the vertexes of the polygon or the polyhedron ∪ T ∈T h T lie on Γ. We denote the interior of this union set by Ω h . Let Γ h be the boundary of Ω h , h T be the diameter of T ∈ T h and h := max T ∈T h h T . We make the very reasonable assumption that every mesh triangle has no more than one edge in Γ h and no mesh tetrahedron has either no more than one face or no more than one edge contained in Γ h . The subset of T h consisting of elements having at least one edge on Γ h is denoted by S h . Now let V h be the finite-element space consisting of continuous functions that vanish on Γ h , whose restriction to each T ∈ T h belongs to P 2 , where P k is the space of polynomials of degree less than or equal to k.

f ∈ L 2 (Ω h ) being an extension of f to Ω h \ Ω, we further set for u, v ∈ H 1 (Ω h ): āh (u, v) := Ω h grad u • grad v dx and L h (v) := Ω h f v dx. (1)
To make ideas clear, and without loss of essential aspects, let us consider the case where g ≡ 0. If we search for ūh such that, ūh ∈ V h and āh (ū

h , v) = L h (v) ∀v ∈ V h , (2)
it is well-known that the energy norm of u -ūh in Ω h , that is, u -ūh ē,h , where

• ē,h := Ω h |grad(•)| 2 dx 1/2 , (3)
will be only an O(h 1.5 ).

In order to recover the optimal O(h 2 )-convergence rate in the energy norm we propose the following. Let W h be a space consisting of functions defined in Ω h whose restriction to every element T ∈ T h is a polynomial of degree less than or equal to two, which are continuous at the vertexes of T and at the mid-points of the edges of T not contained in Γ h , ∀T ∈ T h . Besides this, a function w ∈ W h is required to vanish at all the mesh vertexes lying on Γ h , and at certain points of Γ belonging to the set ∆ S attached to an elements S ∈ S h containing the underlying portion of Γ, constructed as described below ∀S ∈ S h . We consider beforehand that the expression of w ∈ W h in every element S ∈ S h is extended to ∆ S \ Ω h . Referring to Figure 1, in the two-dimensional case ∆ S is the closed set delimited by Γ and the edge d of S contained in Γ h . For every S ∈ S h , the extension of w to ∆ S \ Ω h is required to vanish at a point P ∈ Γ located between two neighboring vertexes of S. To make implementation more straightforward P can be chosen to be the nearest intersection with Γ of the perpendicular to the edge d passing through its mid-point M . Henceforth ∆ S will be referred to as a skin.

Referring to Figures 2 and3 Now instead of solving (2) we search for u h such that,

u h ∈ W h and a h (u h , v) = L h (v) ∀v ∈ V h . ( 4
)
where

a h (w, v) := Ω h grad h w • grad v dx ∀w ∈ W h and ∀v ∈ V h , (5)
grad h : W h + H 1 (Ω h ) -→ [L 2 (Ω h )]
N being the broken gradient operator defined by,

[grad h w] |T = grad[w |T ] ∀T ∈ T h .
Of course if w is continuous the forms a h and āh coincide. This always happens in the two-dimensional case, since W h ⊂ C 0 (Ω h ). However this inclusion is false in three-dimension space because functions w ∈ W h are not necessarily continuous on the interfaces of two tetrahedra in S h (cf. [START_REF] Ruas | Methods of arbitrary optimal order with tetrahedral finite-element meshes forming polyhedral approximations of curved domains[END_REF]). For this reason in the sequel we will work with the broken energy norm • e,h given by,

• e,h :=

Ω h |grad h (•)| 2 dx 1/2 , (6)
and also with the reduced broken energy norm • ẽ,h given by,

• ẽ,h := Ωh |grad h (•)| 2 dx 1/2
where Ωh := Ω h ∩ Ω.

According to [START_REF] Ruas | Optimal simplex finite-element approximations of arbitrary order in curved domains circumventing the isoparametric technique[END_REF] and [START_REF] Ruas | Methods of arbitrary optimal order with tetrahedral finite-element meshes forming polyhedral approximations of curved domains[END_REF],

u h -u e,h (resp. u h -u ẽ,h ) is an O(h 2 ) if Ω is convex (resp. non-convex).
In the remainder of this section we give a detailed proof of these results. One of the keys to the problem is the following proposition:

Proposition 2.1 Provided h is sufficiently small, for every w ∈ W h there exists v ∈ V h such that w -v e,h ≤ C 1 h w e,h , a h (w, v) ≥ w 2 e,h /2 (8)
where C 1 is a strictly positive constant independent of both w and h.

Proof : Given w ∈ W h let v ∈ V h coincide with w at all mesh nodes, except those located on Γ which are not mesh vertexes. The latter are precisely the nodes P whose construction is illustrated in Figures 1 and3, at which w necessarily vanishes. Clearly enough (v -w) |T ≡ 0 for every T having no edge contained in Γ h . Therefore all we need are estimates of S |grad(v -w) |S | 2 dx for every S ∈ S h , i.e. having at least one edge d ⊂ Γ h . For the sake of brevity we only consider the case of an S ∈ S h having exactly one such an edge d; indeed the case of a tetrahedron having a face contained in Γ h can be handled as a trivial extension of the argument that follows.

Let then ϕ S ∈ P 2 be the Lagrange basis function that vanishes at all the five nodes of S for N = 2 or nine nodes for N = 3, which are not the mid-point of d. Denoting such a mid-point by M S we clearly have

[v -w] |S = w(M S )ϕ S . (9)
First we refer to Figures 1 and2 where the set ∆ S is illustrated, for S ∈ S h . Recalling that h S is the diameter of S, by construction there exists a mesh-independent constant C Γ such that

|w(M S )| ≤ C Γ h 2 S max x∈∆ S |grad w(x)|. ( 10
)
Since w is a polynomial in S there must exist another mesh-independent constant C δ such that (cf [START_REF] Strang | An Analysis of the Finite Element Method[END_REF]):

max x∈∆ S |grad w(x)| ≤ C δ max x∈S |grad w(x)|. (11)
Furthermore by a classical inverse inequality (cf. [START_REF] Ciarlet | The Finite Element Method for Elliptic Problems[END_REF]) there exists another mesh-independent constant C ι such that,

max x∈S |grad w(x)| ≤ C ι h -N/2 S S |grad w| 2 dx| 1/2 . ( 12
)
It follows from ( 9), ( 10), ( 11) and ( 12) that

S |grad(v -w)| 2 dx ≤ C 2 ι C 2 δ C 2 Γ h 4-N S S |grad w| 2 dx S |grad ϕ S | 2 dx. (13)
Combining [START_REF] Douglas | An Absolutely Stabilized Finite Element Method for the Stokes Problem[END_REF] with the obvious estimate

S |grad ϕ S | 2 dx ≤ C 2 ϕ h N -2 S ( 14
)
for a suitable mesh-independent constant C ϕ , we easily obtain

v -w e,h ≤ C 1 h w e,h with C 1 = C ϕ C ι C δ C Γ . (15)
Finally, noting that a h (w, v) = w 2 e,h +a h (w, v -w) ≥ w e,h ( w e,h -v -w e,h ), the lower bound,

a h (w, v) ≥ w 2 e,h /2 (16)
trivially follows from ( 15) for h = h 0 , where h 0 is the largest mesh size in the family

P such that 1 -C 1 h 0 is bounded below by 1/2.
As an immediate consequence of ( 8) a uniform inf-sup Babuška-Brezzi condition in connection with problem (4)-( 5) is satisfied. More precisely we have: Corollary 2.2 Provided h is sufficiently small, it holds for some α ≥ 1/3:

∀w ∈ W h = 0, sup v∈V h \{0} a h (w, v) w e,h v e,h ≥ α. (17)
Since obviously dim(V h ) = dim(W h ), and both a h and L h are uniformly bounded independently of h, the simple fact that (17) holds implies that ( 4) is uniquely solvable (cf. [START_REF] Cuminato | Unification of distance inequalities for linear variational problems[END_REF]).

Next we prove error estimates for problem [START_REF] Bertrand | First-order system least-squares on curved boundaries: higherorder Raviart-Thomas elements[END_REF]. In this aim we denote by | • | m,D the usual semi-norm of Sobolev space H m (D) in a bounded domain D ∈ N , for m ∈ N (cf. [START_REF] Adams | Sobolev spaces[END_REF]). First we have: Theorem 2.3 Assume that Ω is convex. Then provided h is sufficiently small, for a certain mesh-independent constant C 2 it holds:

u h -u e,h ≤ C 2 h 2 |u| 3,Ω . (18)
where • e,h is the discrete H 1 -semi-norm defined by [START_REF] Boffi | Three-dimensional finite element methods for the Stokes problem[END_REF].

Proof : From ( 17) we infer that

u h -w e,h ≤ 3 sup v∈V h \{0} a h (u h -w, v) v e,h ∀w ∈ W h . (19)
Let I h (u) ∈ W h be the standard interpolate of u at the mesh nodes associated with W h . Taking in (17) w = u h -I h (u), we add and subtract u in the first argument of a h . Thus by straightforward calculations,

u h -I h (u) e,h ≤ 3 u -I h (u) e,h + sup v∈V h \{0} a h (u h -u, v) v e,h . (20)
Noting that a h (u h , v) = L h (v) we come up with:

u h -I h (u) e,h ≤ 3 u -I h (u) e,h + sup v∈V h \{0} |a h (u, v) -L h (v)| v e,h . (21)
Since

Ω h ⊂ Ω if Ω is convex, we observe that a h (u, v) = Γ h v ∂u ∂n h dΓ h - Ω h
v∆u dx, where ∂u ∂n h is the outer normal derivative of u on Γ h . Noting that -∆u = f in Ω h and v ≡ 0 on Γ h , we trivially obtain,

u h -u e,h ≤ 4 u -I h (u) e,h . (22)
Then ( 18) is a consequence of standard estimates of the interpolation error in Sobolev norms (cf. [START_REF] Brenner | The Mathematical Theory of Finite Element Methods[END_REF]).

Remark 1 It is interesting to note that although our method is non-conforming for N = 3, in any case V h is a subspace of H 1 (Ω h ). Therefore the variational residual a h (u, v) -L h (v) vanishes if Ω is convex, in contrast to usual non-conforming methods.

Now we address the case of a non-convex Ω. Let us consider a smooth domain Ω close to Ω which strictly contains Ω ∪ Ω h for all h sufficiently small. More precisely, denoting by Γ the boundary of Ω we assume that meas( Γ) -meas(Γ) ≤ ε for ε sufficiently small. For the sake of simplicity henceforth we consider that f was extended by zero to the whole Ω \ Ω and still denote this extension by f . In doing so the following error estimate can be proved: Theorem 2.4 Assume that there exists a function ũ defined in Ω, satisfying:

• ũ|Ω = u; • ũ = 0 a.e. on Γ; • ũ ∈ H 3 ( Ω).
Then as long as h is sufficiently small it holds:

u h -u ẽ,h ≤ C2 [h 2 |ũ| 3, Ω + h 5/2 ∆ũ 0, Ω], ( 23
)
where C2 is a mesh-independent constant and • ẽ,h denotes the standard H 1 -semi-norm defined in [START_REF] Brenner | The Mathematical Theory of Finite Element Methods[END_REF].

Proof : According to (17) ∀w ∈ W h we have,

u h -w e,h ≤ 1 α sup v∈V h \{0} |a h (ũ, v) -F h (v)| + |a h (ũ -w, v)| v e,h . (24)
Since ũ ∈ H 3 ( Ω) we can apply First Green's identity to a h (ũ, v) thereby getting rid of integrals on portions of Γ.

Next we note that ∆ũ + f = 0 in every T ∈ T h \ S h ; this is also true of elements T not belonging to the subset R h of S h consisting of elements R such that R \ Ω is not restricted to a set of vertexes of Ω h . Finally we recall that ∆ũ + f vanishes identically in the set R ∩ Ω ∀R ∈ R h . Denoting by ∆R the interior of the set R \ Ω for R ∈ R h we can write:

|a h (ũ, v) -F h (v)| = R∈R h ∆R ∆ũ v dx ≤ R∈R h ∆ũ 0, ∆R v 0, ∆R . (25)
Now taking into account that v ≡ 0 on Γ h and recalling the constant C Γ defined in Proposition 8, it holds: grad v 0,R for a mesh-independent constant C I . Noticing that the measure of ∆R is bounded by a constant depending only on Ω times h N +1 R , after straightforward calculations we obtain for a certain mesh-independent constant C Z :

|v(x)| ≤ C Γ h 2 R |grad v| 0,∞,
∆ũ 0, ∆R v 0, ∆R ≤ C Z h 5/2 R ∆ũ 0, ∆R grad v 0,R ∀R ∈ R h . ( 26
)
Now plugging ( 26) into [START_REF] Scott | Finite Element Techniques for Curved Boundaries[END_REF] and applying the Cauchy-Schwarz inequality, we easily come up with,

|a h (ũ, v) -F h (v)| ≤ C Z h 5/2 ∆ũ 0, Ω v e,h . (27)
Finally plugging ( 27) into [START_REF] Ruas | Optimal Calculation of Solid-Body Deformations with Prescribed Degrees of Freedom over Smooth Boundaries[END_REF] and taking w = I h (ũ), we immediately establish the validity of error estimate [START_REF] Ruas | Variational formulations yielding high-order finite-element solutions in smooth domains without curved elements[END_REF].

Remark 2 There are many ways to ensure the existence of ũ satisfying the assumptions of Theorem 2.4, as long as Γ is as smooth as required. For instance we refer to [START_REF] Stein | Immersed boundary smooth extension: A high-order method for solving PDE on arbitrary smooth domains using Fourier spectral methods[END_REF] for an interesting construction of ũ.

Some numerics for the Poisson equation

Let us illustrate the performance of the new technique to handle Dirichlet conditions on curved boundaries. With this aim we first solve problems (2) and (4) in case Ω is the unit disk centered at the origin, and a uniformly regular family of meshes consisting of 8n 2 triangles for n = 2 m , with m = 1, 2, . . . is constructed in the way described in [START_REF] Ruas | Methods of arbitrary optimal order with tetrahedral finite-element meshes forming polyhedral approximations of curved domains[END_REF]. In these experiments we take f (x, y) = 9r where r = (x 2 + y 2 ) 1/2 , and hence the exact solution is given by u(x, y) = 1 -r 3 . Owing to symmetry only the quarter disk corresponding to x > 0 and y > 0 is taken into account in the computations, and therefore meshes containing 2n 2 elements are employed. For a fairer comparison we also supply results obtained for the same problem solved by the classical isoparametric technique. We denote the solution obtained with this method by ũh . Taking m = 2, 3, 4, 5, 6 and observing that h = 1/n, in Table 1 the quantities u h -u e,h , ũh -u e,h and ūh -u e,h for the resulting decreasing values of h are displayed.

u h -u e,h -→ 0.1329×10 -1 0.3343×10 -2 0.8381×10 -3 0.2097×10 -3 0.5245×10 -4 ūh -u e,h -→ 0.5434×10 -1 0.1969×10 -1 0.7042×10 -2 0.2503×10 -2 0.8870×10 -3 ũh -u e,h -→ 0.1559×10 -1 0.3837×10 -2 0.9477×10 -3 0.2353×10 -3 0.5861×10 -4
Table 1: Energy errors for a test-problem in a disk solved by methods (4), (2) and isoparametric FEs in the energy norm for the approach advocated in this paper, while the polygonal approach (2) yields only O(h 1.5 ) approximations in the same norm, as predicted in classical books (cf. [START_REF] Ciarlet | The Finite Element Method for Elliptic Problems[END_REF]). Of course the expected second order convergence of the isoparametric solution is also observed. However the new method is a little more accurate. Now in order to further illustrate the accuracy of method (4) in case Ω is not convex, we compare it again with method (2) by solving a problem whose exact solution is not axisymmetric. More specifically, here the domain described in polar coordinates (r, θ) is given by Ω := {(r, θ)| r ≤ [4 + cos(4θ)]/5}. Taking g ≡ 0 and f = 16r 2 -5.8r -11.2x 2 y 2 /r 3 the exact solution is the function u = r 3 -r 4 -1.6x 2 y 2 /r. Notice that f ∈ H 1 (Ω) and u ∈ H 3 (Ω).

Here again symmetry allows working with the computational domain corresponding to x > 0 and y > 0. The meshes also consist of 2n 2 elements, generated like in [START_REF] Ruas | A simple alternative for accurate finite-element modeling in curved domains[END_REF], by subdividing the radial coordinate r into n equal parts and the azimuthal coordinate θ ∈ (0, π/2) into 2n equal parts.

Observing that h = 1/n, we show in Table 2 the quantity u h -u e,h for n = 2 m , with m = 2, 3, 4, 5, 6. Moreover, in order to give a better idea of how effective our method is, we also supply the errors u h -u 0,h and |u h -u| ∞,h , where • 0,h and | • | ∞,h stand for the standard norm of L 2 (Ω h ) and the maximum absolute error at the nodal points, respectively. in Theorem 2.4. Even better news come from the observed convergence rates of three in the norm of L 2 (Ω h ) and of a little less than three in the

u h -u e,h -→ 0.1566×10 -1 0.4235×10 -2 0.1089×10 -2 0.2751×10 -3 0.6910×10 -4 u h -u 0,h -→ 0.4875×10 -3 0.5859×10 -4 0.7315×10 -5 0.9195×10 -6 0.1156×10 -6 |u h -u| ∞,h -→ 0.8110×10 -3 0.1790×10 -3 0.2745×10 -4 0.3701×10 -5 0.4915×10 -6
L ∞ -semi-norm | • | ∞,h .
4 Application to the Taylor-Hood element

The classical Taylor-Hood element was introduced in [START_REF] Hood | Navier-Stokes equation using mixed interpolation[END_REF] for the solution of the incompressible Navier-Stokes equations. It consists of continuous piecewise polynomial representations of both velocity and pressure in triangles or tetrahedra, of degree two for the former variable and of degree one for the latter. Second order convergence results for this method were established by Verfürth [START_REF] Verfürth | Error estimates for a mixed finite element approximation of the Stokes equations[END_REF] in the case of a polygonal domain and by Boffi [START_REF] Boffi | Stability of higher order triangular Hood-Taylor methods for the stationary Stokes equations[END_REF] in the case of polyhedrons. In this section we apply the method described in Section 2 in order to extend such results to the case of smooth curvilinear domains. This study will be restricted to the linearized form of the stationary incompressible Navier-Stokes equations, which governs incompressible viscous flows at a very low Reynolds number. More specifically our theory applies to the following Stokes system in a bounded domain Ω of N at least of the C 1 -class, for N = 2 or N = 3: Given a field f ∈ [H 1 (Ω)] N , and a velocity profile g defined on Γ assumed to belong to [H 5/2 (Γ)] N and to satisfy the conservation property Γ g • n ds = 0, where n is the unit outer normal vector on Γ, we wish to determine a velocity field u ∈ [H 1 (Ω)] N and a hydrostatic pressure p ∈ L 2 (Ω)/ , where A/B denotes the quotient between two vector spaces A and B, such that:

   -∆u + grad p = f div u = 0 in Ω u = g on Γ. (28)
A suitable regularity assumption on Ω, besides those applying to f and g, legitimately allows assuming in turn that u ∈ [H 3 (Ω)] N and p ∈ H 2 (Ω). Although all the results to be derived hereafter apply to the inhomogeneous case, in order to avoid non essential difficulties, we further restrict the analysis conducted in this section to the case where g ≡ -→ 0 .

Our working spaces here will be the pair (V h , W h ) of vector field spaces defined by

V h := [V h ] N and W h := [W h ] N , together with the function space Q h := Qh ∩ L 2 0 (Ω h ), with Qh := {q | q ∈ C 0 (Ω h ), q |T ∈ P 1 , ∀T ∈ T h }, where L 2 0 (Ω h ) = {q | q ∈ L 2 (Ω h ), Ω h q dx = 0}.
For the sake of simplicity, henceforth we denote by |w| 1,h the semi-norm grad h w 0,h of a field w ∈ W h + [H 1 (Ω h )] N , where • 0,h stands for standard norm of L 2 (Ω h ). We make the same assumptions as in Section 2 on a given family P of meshes T h of Ω into N -simplexes. In doing so we consider the extension by zero f of f to Ω h \ Ω, if applicable, and define the broken divergence operator div h :

W h + [H 1 (Ω h )] N -→ L 2 (Ω h ) by [div h w] |T = div w |T ∀T ∈ T h . We further set for w ∈ W h + [H 1 (Ω h )] N , v ∈ [H 1 (Ω h )] N and q ∈ L 2 (Ω h ):        c h (w, v) := Ω h grad h w : grad v dx b h (v, q) := -Ω h q div v dx d h (w, q) := -Ω h q div h w dx L h (v) := Ω h f • v dx. (29)
Now we pose the corresponding finite-element counterpart of (28) as:

   Find u h ∈ W h and p h ∈ Q h such that: c h (u h , v) + b h (v, p h ) = L h (v) ∀v ∈ V h , d h (u h , q) = 0 ∀q ∈ Q h . (30)
According to the classical theory of linear variational problems (see e.g. [START_REF] Cuminato | Unification of distance inequalities for linear variational problems[END_REF]) problem ( 30) is well-posed thanks to the validity of the underlying Babuška-Brezzi condition, or yet the inf-sup condition [START_REF] Zienkiewicz | The Finite Element Method in Engineering Science[END_REF], that is, Proposition 4.1 Provided h is sufficiently small, there exists a strictly positive constant A independent of h such that

inf (w,p)∈W h ×Q h \{( - → 0 ,0)} sup (v,q)∈V h ×Q h \{( - → 0 ,0)} c h (w, v) + b h (v, p) + d h (w, q) [|w| 2 1,h + p 2 0,h ] 1/2 [|v| 2 1,h + q 2 0,h ] 1/2 ≥ A. (31)
Proof: Let the pair (w, p) = ( -→ 0 , 0) be given in W h × Q h . First we observe that, since Taylor-Hood elements are uniformly stable, the following condition holds:

sup v∈V h \{ - → 0 } b h (v, p) |v| 1,h ≥ β p 0,h . (32)
for a constant β > 0 independent of both p and the mesh. Actually [START_REF]Curved Elements in the Finite Element Method[END_REF] is the consequence of well-known arguments (cf. [START_REF] Brezzi | Mixed and Hybrid Finite Element Methods[END_REF]), according to which there exist two mesh-independent constants C 3 and C 4 also independent of p, such that one can find

v 0 ∈ V h satisfying b h (v 0 , p) ≥ C 3 p 2 0,h |v 0 | 1,h ≤ C 4 p 0,h . (33)
Noticing that c h is nothing but a h applied to vector fields instead of functions, let v 1 ∈ V h satisfy the obvious vector analog of (8) for our given w ∈ W h . For a certain parameter η > 0 we define v := ηv 0 + v 1 and take q ≡ -p. From the obvious vector analog of [START_REF] Hazewinkel | Diffeomorphism[END_REF] we easily obtain

|v 1 | 1,h ≤ (1 + C 1 h)|w| 1,h ,
which together with (33) immediately yields:

   c h (w, v)+b h (v, p)+d h (w, q) ≥ |w| 2 1,h 2 -ηC 4 |w| 1,h p 0,h +b h (v 1 , p)-d h (w, p)+ηC 3 p 2 0,h and |v| 1,h + q 0,h ≤ (1 + C 1 h)|w| 1,h + (ηC 4 + 1) p 0,h . (34)
Next we note that

b h (v 1 , p) -d h (w, p) = Ω h p div h (w -v 1 )dx ≤ √ N |w -v 1 | 1,h p 0,h .
Thus using Young's inequality and recalling that C 1 h ≤ 1/2, from (34) we obtain,

     c h (w, v) + b h (v, p) + d h (w, q) ≥ |w| 2 1,h 4 - √ N |w -v 1 | 1,h p 0,h +(ηC 3 -η 2 C 2 4 ) p 2 0,h and |v| 1,h + q 0,h ≤ 3|w| 1,h 2 + (ηC 4 + 1) p 0,h . (35)
Moreover, plugging the natural vector version of ( 15) into the first inequality of (35), we derive

     c h (w, v) + b h (v, p) + d h (w, q) ≥ |w| 2 1,h 4 - √ N C 1 h|w| 1,h p 0,h +(ηC 3 -η 2 C 2 4 ) p 2 0,h and |v| 1,h + q 0,h ≤ 3|w| 1,h 2 + (ηC 4 + 1) p 0,h . (36)
Taking η = C 3 /(2C

c h (w, v) + b h (v, p) + d h (w, q) ≥ C 5 (|w| 2 1,h + p 2 0,h ) and [|v| 2 1,h + q 2 0,h ] 1/2 ≤ C 6 [|w| 2 1,h + p 2 0,h ] 1/2 . ( 37
)
with

C 6 = {9/4 + [C 3 /(2C 4 ) + 1] 2 } 1/2 .
In view of both inequalities in (37), as long as

h ≤ min[2C 5 / √ N , 1/2]/C 1 , ( 31
) holds with A = C 5 /C 6 .
Now we endeavor to derive error estimates for problem [START_REF] Žénišek | Curved triangular finite C m -elements[END_REF]. Essentially this task is not more complicated than the one carried out in Theorems 2.3 and 2.4. Indeed [START_REF] Žénišek | Curved triangular finite C m -elements[END_REF] can be rewritten as follows:

Find U h ∈ W h such that A h (U h , V ) = L h (V ) ∀V ∈ V, (38)
where

• U h = (u h , p h ); • V = (v, q); • W h := W h × Q h ; • V h := V h × Q h ; • A h ((w, p), (v, q)) := c h (w, v) + b h (v, p) + d h (w, q); • L h (V ) := L h (v).

Now we denote by • X,h the norm over {W

h + [H 1 0 (Ω h )] N } × L 2 (Ω h ), given by V X,h := [|v| 2 1,h + q 2 0,h ] 1/2 . ( 39
)
Then letting A play the same role as the constant 1/3 in ( 17), analogously to [START_REF] Ruas | A Hermite method for Maxwell's equations[END_REF] we obtain:

U h -W X,h ≤ 1 A sup V ∈V h \{O} A h (U h -W, V ) V X,h ∀W ∈ W h . ( 40
)
Finally noticing that here also the variational residual A h ((u, p), (v, q)) -L h ((v, q)) vanishes for every (v, q) ∈ V h if Ω is convex, using standard estimates for the interpolation error in Sobolev spaces, akin to Theorem 2.3, (40) leads to: Theorem 4.2 Provided h is small enough and Ω is convex, for a certain mesh-independent constant C it holds:

[|u -u h | 2 1,h + p -p h 2 0,h ] 1/2 ≤ Ch 2 [|u| 3,Ω + |p| 2,Ω ]. ( 41
)
The case where Ω is not convex can be treated quite in the same manner as in Section 2. The key to the problem is the existence of suitable extensions ũ of u and p of p to the domain Ω \ Ω, where Ω is defined in Section 2. More precisely, we extend f by zero to Ω \ Ω and still denote by f such an extension. However, naturally enough, more technicalities come into play here.

To begin with we need the following preliminary result:

Lemma 4.3 Let φ be a function in H 1 ( Ω) that vanishes on Γ. There exists a mesh-independent constant C X such that φ 0, ∆R ≤ C X h 2 R |φ| 1, ∆R ∀R ∈ R h . (42)
Proof : We refer to [START_REF] Hazewinkel | Diffeomorphism[END_REF] for the terminology and some properties of diffeomorphisms used in this proof.

Let us cover the whole non-convex region of Γ by a set of M overlapping local maps, say, ω i , i = 1, . . . , M . Owing to our regularity assumptions, there exists a C 1 -diffeomorphism F i that transforms ω i into a set ωi such that Γi := F i (Γ ∩ ω i ) is a line segment for N = 2 or a plane bounded set for N = 3. Without loss of generality we assume that the measure of Γi is not zero, and moreover that we can assign each R ∈ R h to a certain local map ω i , in such a way that R ⊂ ω i . We generically denote by x = ( t, n) the local Cartesian coordinate system of N with coordinates t = t for N = 2 or t = ( t1 , t2 ) for N = 3 along or upon Γi , and by n the coordinate along the axis orthogonal to Γi oriented from this manifold outwards the image of ω i ∩ Ω under F i . Let R ∈ R h and ∆R be the transformation of ∆R under F i for the appropriate i. We denote by φ the transformation of φ under F i defined in ωi . Since φ = 0 on Γi we can write

φ( t, n) = ν=n ν=0 ∂ φ ∂ n ( t, ν) dν.
Hence, we obtain successively,

∆R | φ| 2 dx ≤ ∆R ν=n ν=0 ∂ φ ∂ n ( t, ν) dν 2 d tdn, ∆R | φ| 2 dx ≤ ∆R l( t)    ν=l( t) ν=0 ∂ φ ∂ n ( t, ν) 2 dν    d tdn,
where l( t) is the width of ∆R measured in the direction normal to Γi from point ( t, 0). Then denoting by l the maximum of l( t) over ( t, 0) ∈ Γi ∩ ∆R , we trivially obtain,

∆R | φ| 2 dx ≤ l2 ∆R ∂ φ ∂ n ( t, n) 2 dnd t,
and further,

∆R | φ| 2 dx ≤ Ĉi h 2 R ∆R | grad φ| 2 dx
where grad(•) represents the gradient operator of a function defined in ωi , and the constant Ĉi depends only on Ω and ωi . Next we make straightforward changes of variables in the above integrals, thereby transforming them into integrals in ∆R , and observe that grad φ = F -1 i grad φ where F i is the Jacobian matrix of F i . From a basic property of diffeomorphisms the spectral norm of F i can be uniformly bounded above by a constant independent of the mesh, as much as the Jacobian of both F i and F -1 i . Finally taking the extrema over i of those constants and of Ĉi in the required senses, the result follows.

Now we have

Theorem 4.4 Assume that there exists ũ and p satisfying the following conditions:

• ũ|Ω = u and p|Ω = p • ũ = - → 0 a.e. on Γ; • ũ ∈ [H 3 ( Ω)] N and p ∈ H 2 ( Ω).
Then, as long as h is small enough, for a certain mesh-independent constant C it holds:

[|u -u h | 2 1, Ωh + p -p h 2 Ωh ] 1/2 ≤ C{h 2 [|ũ| 3, Ω + |p| 2, Ω] + h 5/2 [|ũ| 2, Ω + |p| 1, Ω]}, (43)
where Ωh := Ω h ∩ Ω.

Proof : The proof of this theorem is based on the same arguments as the proof of Theorem 2.4. Therefore we skip some details. First we set Ũ := (ũ, p). For every W = (w, r) ∈ W h we have:

U h -W X,h ≤ 1 A sup V =(v,q)∈V h =O |A h ( Ũ , V ) -L h (V )| + |A h (W, V ) -A h ( Ũ , V )| V X,h . (44)
The second term in the numerator of (44) can be handled in a standard manner by means of classical interpolation theory. This yields for a mesh-independent constant C J :

inf W ∈W h |A h (W, V ) -A h ( Ũ , V )| V X,h ≤ C J h 2 [|ũ| 3, Ω + |p| 2, Ω]. (45)
Next, thanks to the fact that ũ ∈ [H 3 ( Ω)] N and p ∈ H 2 ( Ω) we can write,

|A h ( Ũ , V ) -L h (V )| ≤ R∈R h ∆R [|∆ũ| + |grad p|]|v| dx + ∆R |q| |div ũ| dx . (46)
Similarly to Theorem 2.4 the summation of the first integral on the right hand side of (46) can be bounded above as follows, for a suitable mesh-independent constant C Y :

R∈R h ∆R [|∆ũ| + |grad p|]|v| dx ≤ C Y h 5/2 [ ∆ũ 0, Ω +|p| 1, Ω]|v| 1,h . (47)
On the other hand we have, ∆R |q| |div ũ| dx ≤ q 0, ∆R div ũ 0, ∆R .

Now, since div ũ vanishes on Γ, using Lemma 4.3, it holds for a certain mesh-independent constant C X :

div ũ 0, ∆R ≤ C X h 2 R |div ũ| 1, ∆R (49)
Moreover using the fact that meas( ∆R ) 48) and (49) we derive,

≤ C Q h N +1 R for a mesh-independent constant C Q , together with the inverse inequality q 0,∞,R ≤ C I h -N/2 R q 0,R , from (
∆R |q| |div ũ| dx ≤ C I C Q C X h 5/2 R |div ũ| 1, ∆R q 0,R . (50)
This trivially yields

R∈R h ∆R |q| |div ũ| dx ≤ √ N C I C Q C X h 5/2 q 0,h |ũ| 2, Ω. (51)
Finally combining (44), (45), ( 46), ( 47) and (51) we come up with (43).

As pointed out in Remark 2, the construction of a pair (ũ, p) satisfying the assumptions of Theorem 4.4 can be performed in different manners. In this respect we refer for instance to [START_REF] Stein | Immersed boundary smooth extension: A high-order method for solving PDE on arbitrary smooth domains using Fourier spectral methods[END_REF].

Numerical validation for confined rotating flows

One of the most remarkable applications of the method studied in the previous section is the simulation of confined rotating flows. Indeed in this case a viscous fluid adhere to the curved wall of the flow region, and thus handling the underlying Dirichlet boundary condition with a method of order higher than one requires the use of an accurate technique. In this section we present results obtained with ours, for two test-problems governed by the Stokes system.

In the tables of this section the acronym OCR stands for observed convergence rate.

Test-problem with a manufactured solution

First we apply the Taylor-Hood method combined with our technique to solve [START_REF] Strang | An Analysis of the Finite Element Method[END_REF] with a manufactured solution corresponding to the following data: Ω is the unit disk (centered at the origin), f = (8, 8)(x -y), and g ≡ -→ 0 . Prescribing p( √ 2/2, √ 2/2) = 0, the exact solution has polynomial expressions, namely u = (y, -x)(1-x 2 -y 2 ) and p = x 2 -y 2 . We use meshes constructed like in the first test-problem of Section 3, but here the computational domain is the whole disk. More specifically now we compute with (2n × 2n)-meshes containing 8n 2 triangles, each mesh being symmetric with respect to the axes x = 0 and y = 0, for n = 2 m with m = 2, 3, 4, 5. We recall that h = 1/n.

In order to discard any particularity inherent to the problem being solved, we compared the numerical solution with the one obtained by the simple polygonal approach. We display in Table 3 the velocity and pressure errors in the norms | • | 1,h and • 0,h for both approaches. The notations ūh and ph are employed to represent the velocity and pressure obtained by the simple polygonal approach. These results completely validate the analysis carried out in the previous sections for the case of a convex domain (cf. Theorem 4.2). It is no surprise that the polygonal approach does erode the order of the velocity approximation. It is interesting to note however that, at least in this test-case, this simple approach does not affect the pressure approximation.

2n

-→ 8 16 32 64 OCR

|u h -u| 1,h -→ 0.3585×10 -1 0.8833×10 -2 0.2157×10 -2 0.5299×10 -3 O(h 2 ) |ū h -u| 1,h -→ 0.7959×10 -1 0.2746×10 -1 0.9588×10 -2 0.3370×10 -2 O(h 1.5 ) p h -p 0,h -→ 0.4266×10 -1 0.1047×10 -1 0.2567×10 -2 0.6332×10 -3 O(h 2 ) ph -p 0,h -→ 0.4264×10 -1 0.1047×10 -1 0.2567×10 -2 0.6332×10 -3 O(h 2 )
Table 3: Errors for a test-(flow) problem in a disk solved by method [START_REF] Žénišek | Curved triangular finite C m -elements[END_REF] and the polygonal approach

Pseudo circular Couette flow

In order to check our method's performance in the case of a non-convex flow domain we used it to solve the problem described as follows.

Circular Couette flow of an incompressible viscous fluid with density ρ in a region comprised between two concentric cylinders, where the inner one of radius r i rotates at an angular velocity ω and the outer one with radius r e is kept fixed, is governed by the stationary Navier-Stokes equations with a zero body-force right hand side. As long as the Reynolds number is sufficiently low, the flow is laminar and the solution to the problem is given by u = (sinθ, -cosθ)u θ (r) where u θ (r) = ωr 2 i (r 2 e -r 2 )/[r(re 2 -ri 2 )] and p(r) = ρω 2 r 4 i /(r 2 e -r 2 i ) 2 [r 2 /2r 4 e /(2r 2 ) -2r 2 e log(r)] + c, c being a constant. If we enforce zero pressure on the outer wall, then c takes the value 2r 2 e log(r e )ρω 2 r 4 i /(r 2 e -r 2 i ) 2 . Although there is no particular difficulty to solve the Navier-Stokes equations with our method, in order to focus on our essentially validating goal, we apply it to a modified problem, in which the exact inertia term ρ[grad u]u with a minus sign is input as right hand side datum f . Of course the pair (u, p) is still the solution to the resulting Stokes system [START_REF] Strang | An Analysis of the Finite Element Method[END_REF] in the annulus Ω with inner radius r i and outer radius r e . The datum g in turn equals -→ 0 for r = r e , while its value for r = r i conforms to the given azimuthal velocity r i ω and a zero radial velocity.

Taking r e = 1, r i = 0.5, ω = 1 and ρ = 1, we proceeded to the numerical solution of thus defined (pseudo) circular Couette flow problem with the Taylor-Hood method combined with our technique to approximate the boundary conditions. In order to avoid non physical boundary conditions, computations were carried out for the whole annulus. With this aim we used again (2n × 2n) symmetric meshes, for n = 2 m , with m = 3, 4, 5, 6, constructed in the way described in the previous subsection, except for the fact that now the elements inside the disk with radius r i were disregarded. This yields meshes consisting of 6n 2 triangles, with h = 1/n. We display in Table 4 the velocity errors measured in the norms | • | 1,h and • 0,h , together with the pressure errors measured in the • 0,h -norm. It is interesting to note that the latter are decreasing at a rate faster than the O(h 2 ) observed in the test-problem of the previous sub-section. This seems to be due to the fact that in circular Couette flow the inertia term with a minus sign is nothing but the pressure gradient. On the other hand the velocity errors in the H 1 -semi-norm are in perfect agreement with the theoretical predictions. The velocity errors in the L 2 -norm in turn seem to decrease like an O(h 3 ), which is optimal.

2n

-→ 16 32 64 128 OCR |u h -u| 1,h -→ 0.1592×10 +0 0.4261×10 -1 0.1090×10 -1 0.2741×10 -2 O(h 2 ) u h -u 0,h -→ 0.3833×10 -2 0.5339×10 -3 0.6923×10 -4 0.8744×10 -5 O(h 3 ) p h -p 0,h -→ 0.1209×10 +0 0.2095×10 -1 0.3952×10 -2 0.6638×10 -3 O(h ≈2.5 )
Table 4: Errors for the pseudo circular Couette flow problem solved by method [START_REF] Žénišek | Curved triangular finite C m -elements[END_REF] equation as a basis for the solution of incompressible viscous flow problems by the popular Taylor-Hood element. However this second order mixed finite element was chosen here only for illustrative purposes. As a matter of fact our technique basically applies to most known reliable mixed methods of order greater than one, to solve this kind of problems, as long as velocity degrees of freedom must be prescribed at boundary points different from vertexes.

Let us be more specific about some of these methods.

1. If we stick to second-order methods based on the standard Galerkin formulation (such as Taylor-Hood elements), the convergence results that apply to the Crouzeix-Raviart method on triangles [START_REF] Crouzeix | Conforming and nonconforming finite element methods for solving the stationary Stokes equations[END_REF] for the polygonal case extend to method's obvious modification using our technique. Notice however that the application of this technique to the Crouzeix-Raviart method's extension to tetrahedra considered in [START_REF] Ruas | Finite Element Solution of 3D Viscous Flow Problems Using Non-standard Degrees of Freedom Japan[END_REF] must be the object of a specific study. This is because it employs certain mean values along element edges as velocity degrees of freedom, instead of nodal values.

2. Methods using a piecewise quadratic representation of the velocity combined with the Petrov-Galerkin formulation due to Franca & Hughes [START_REF] Franca | Stabilized Finite Element Methods[END_REF] or the one of Douglas & Wang [START_REF] Douglas | An Absolutely Stabilized Finite Element Method for the Stokes Problem[END_REF] can be combined with our technique quite in the same manner as Taylor-Hood elements. The final (second-order) qualitative results remain unchanged.

3. Any third-order method in the natural norms using a cubic velocity representation can also be optimally handled in association with our technique. This is true of Taylor-Hood element's cubic extension using the standard Galerkin formulation considered by Boffi [START_REF] Boffi | Three-dimensional finite element methods for the Stokes problem[END_REF], and also of the method in the Petrov-Galerkin formulation mentioned in the previous item. In both cases the analysis is based on the arguments developed in [START_REF] Ruas | Optimal simplex finite-element approximations of arbitrary order in curved domains circumventing the isoparametric technique[END_REF] and [START_REF] Ruas | Methods of arbitrary optimal order with tetrahedral finite-element meshes forming polyhedral approximations of curved domains[END_REF] for cubic Lagrange finite elements.

4. Methods of order k ≥ 4 in the natural norm, though of limited interest, can also be combined with our technique. More particularly this is the case of the generalized Taylor-Hood pairing consisting of the continuous P K -P k-1 velocity-pressure representation considered in [START_REF] Boffi | Three-dimensional finite element methods for the Stokes problem[END_REF], or yet the continuous P k -discontinuous P k-1 velocity-pressure method studied by Scott & Vogelius [START_REF] Scott | Norm estimates for a maximal right inverse of the divergence operator in spaces of piecewise polynomials[END_REF]. However here optimal convergence results hold under the condition that a numerical quadrature formula with a compatible order and without integration points in the interior of boundary edges or faces be employed to compute the right hand side term. We refer to [START_REF] Ruas | Optimal simplex finite-element approximations of arbitrary order in curved domains circumventing the isoparametric technique[END_REF] and [START_REF] Ruas | Methods of arbitrary optimal order with tetrahedral finite-element meshes forming polyhedral approximations of curved domains[END_REF] for more details about such a restriction, which also applies to isoparametric elements (cf. [START_REF] Ciarlet | The combined effect of curved boundaries and numerical integration in isoparametric finite element methods[END_REF]).

In conclusion the author emphasizes that the scope in Computational Engineering of the approach adopted in this work to handle Dirichlet conditions prescribed on curved boundaries is much wider than the one of classical techniques such as isoparametric finite elements. This was shown in [START_REF] Ruas | A Hermite method for Maxwell's equations[END_REF] and [START_REF] Ruas | A simple alternative for accurate finite-element modeling in curved domains[END_REF] in the framework of Maxwell's equations of Electromagnetism and deformations of elastic membranes in mixed formulation, respectively. Moreover, even in cases where the use of classical techniques is consolidated, our approach is at least as reliable and competitive in terms of accuracy. We would also like to point out that, as far as we can see, our technique has only two drawbacks: first of all it is necessary to solve a non symmetric problem, even when the original problem is symmetric. Moreover for each boundary element a small matrix has to be inverted in order to determine the local basis functions. However none of both issues are a real problem nowadays, taking into account the state-of-the-art of Computational Linear Algebra. A final remark on the choice of nodal points on Γ different from vertexes is in order. As one can easily infer from the analysis carried out throughout the paper, the construction of these nodes advocated in Section 2 is not compulsory at all. Actually, referring to Figures 1 and3, any other choice in ∆ S ∩ Γ for S ∈ S h will do. However intuitively we can say that these nodes should not be too close to the boundary vertexes of element S, since this may lead to a worse conditioning of the resulting linear system. Remark 3 Besides those considered in [START_REF] Ruas | A simple alternative for accurate finite-element modeling in curved domains[END_REF], applications to Solid Mechanics of the technique studied in this paper can be found in [START_REF] Ruas | Optimal Calculation of Solid-Body Deformations with Prescribed Degrees of Freedom over Smooth Boundaries[END_REF].

  (right), in the three-dimensional case, for every boundary edge d of an element S ∈ S h we first denote by Π the plane bisecting the dihedral formed by the faces F and F of S and another tetrahedron S ∈ S h respectively, whose intersection is d. We generically denote by δ d the closed subset of Π comprised between d and Γ, referred to as a plane skin hereafter, as depicted in Figure2for a tetrahedron S having a face F contained in Γ h . Referring to Figure3(left), in the case of such a tetrahedron we denote by ∆ S the closed subset of Ω ∪ Ω h delimited by the three plane skins δ d , Γ and F . For the other type of tetrahedrons S ∈ S h there will be only one such a plane skin δ d and we set ∆ S = δ d . Then for every S ∈ S h the extension of w to ∆ S \ Ω h is required to vanish at a point P ∈ Γ belonging to δ d located between the end-points of d, for every edge d of S contained in Γ h . Akin to the two-dimensional case P can be conveniently chosen to be the nearest intersection with Γ of the perpendicular to d in δ d passing through its mid-point M .

Figure 1 :

 1 Figure 1: Set ∆ S for triangles S in S h with their nodes P ∈ Γ associated with M ∈ Γ h

Figure 2 :Figure 3 :

 23 Figure 2: Plane skin δ d for two tetrahedra S, S having faces F, F ⊂ Γ h with F ∩ F = edge d

  ∆R , ∀x ∈ ∆R , where • 0,∞,D denotes the standard norm of L ∞ (D), D being a bounded open set of N . Next, using the same arguments as in the proof of Proposition 8, we derive the estimate |grad v| 0,∞, ∆R ≤ C I h -N/2 R

Table 1

 1

	confirms second order convergence

  Table 2 validates method's second order convergence in energy norm established

	n	-→	4	8	16	32	64

Table 2 :

 2 Errors in different senses for a test-problem in a non-convex domain solved by method[START_REF] Bertrand | First-order system least-squares on curved boundaries: higherorder Raviart-Thomas elements[END_REF]

Extensions to other mixed elements and final commentsIn the four previous sections we focused on the application of the technique introduced in[START_REF] Ruas | Optimal simplex finite-element approximations of arbitrary order in curved domains circumventing the isoparametric technique[END_REF] and[START_REF] Ruas | Methods of arbitrary optimal order with tetrahedral finite-element meshes forming polyhedral approximations of curved domains[END_REF] to solve boundary value problems in smooth curved domains with Dirichlet boundary conditions, in the particular case of quadratic Lagrange interpolation in N -simplexes. More specifically we considered the solution of the Poisson

Acknowledgment: The author gratefully acknowledges the financial support provided by CNPq through grant 307996/2008-5.