Vitoriano Ruas
email: vitoriano.ruas@upmc.fr

Marco Antonio

Silva Ramos

Efficiency of nonparametric finite elements for optimal-order enforcement of Dirichlet conditions on curvilinear boundaries

Keywords: Curved Domain, Dirichlet, Finite Element, Iterative solution, Optimal Order, Petrov-Galerkin, Straight-edged simplex

 and the standard isoparametric one. The study of this symmetrization is completed by an optimal error estimation in the broken H 1 -norm for the nonparametric version of the nonconforming method, which had not been addressed in previous work.

Introduction

In the past decades Petrov-Galerkin formulations of boundary value problems showed to be a powerful tool to overcome difficulties brought about by the space discretization of certain types of partial differential equations. A significant illustration is provided by the families of methods proposed by Franca and Hughes and collaborators in the late eighties for the finite-element modeling of various problems in Continuum Mechanics, in particular as a popular alternative to Galerkin methods for viscous incompressible flow (see e.g. [START_REF] Franca | Stabilized Finite Element Methods[END_REF]). The outstanding contributions in the seventies of Babuška (see e.g. [START_REF] Babuška | The finite element method with Lagrange multipliers[END_REF]) and Brezzi [START_REF] Brezzi | On the existence, uniqueness and approximation of saddle-point problems arising from Lagrange multipliers[END_REF], among other authors, were decisive to provide a theoretical background that allowed to formally justify the reliability of Petrov-Galerkin formulations, namely, the so-called inf-sup condition.

In a series of papers published since 2017 (cf. Ruas [START_REF] Ruas | Optimal simplex finite-element approximations of arbitrary order in curved domains circumventing the isoparametric technique[END_REF][START_REF] Ruas | A simple alternative for accurate finite-element modeling in curved domains[END_REF] and Ruas and Silva Ramos [START_REF] Ruas | A Hermite Method for Maxwell's Equations[END_REF]) a nonparametric technique of the Petrov-Galerkin type was introduced, in order to enhance the accuracy of higher order finite element methods to solve boundary value problems with Dirichlet conditions, posed in smooth curved domains. In contrast to parametric elements, it employs straight-edged triangular or tetrahedral meshes fitting the domain. In order to attain best-possible orders greater than one, piecewise polynomial trial-functions are employed, which interpolate the Dirichlet conditions at points of the true boundary. In the two-dimensional case this kind of trial-functions is similar to the one also employed as test-functions by the method known as interpolated boundary conditions studied in [START_REF] Brenner | The Mathematical Theory of Finite Element Methods[END_REF]. However, in spite of being very intuitive and known since the seventies (cf. [START_REF] Nitsche | On Dirichlet problems using subspaces with nearly zero boundary conditions[END_REF] and [START_REF] Scott | Finite Element Techniques for Curved Boundaries[END_REF]), the lack of an extension to three-dimensional problems seems to have inhibited its use among practitioners. In contrast, the test-functions for our method are defined upon the degrees of freedom associated with the underlying finite element method for the mesh forming a polytope equal to the union of straight-edged simplexes. This polytope fits the curved domain in such a manner that all of its vertexes lie on the boundary of the latter. In doing so the integration domain is restricted to this polytope, thereby rendering method's implementation straightforward in both two-and three-dimensional geometries. Moreover only polynomial algebra is necessary, while best-order approximations can be obtained for non-restrictive choices of boundary nodal points.

Generally speaking, the Petrov-Galerkin methodology studied in this work is designed to enforce Dirichlet conditions in the form of prescribed boundary degrees of freedom of various types, in connection with methods of order greater than one in problem's natural norm, for a wide spectrum of boundary value problems. According to numerous numerical experiments reported in previous papers, including those cited above, it showed to be fully reliable in different contexts. It also appeared to be superior to well known techniques to tackle the same kind of problem, in case they exist. For instance in [START_REF] Ruas | Optimal Dirichlet-condition enforcement on curved boundaries for Lagrange and Hermite FEM with straight-edged simplexes[END_REF] and [START_REF] Ruas | Optimal-rate finite-element solution of Dirichlet problems in curved domains with straight-edged tetrahedra[END_REF] comparisons of this method with the isoparametric version of the finite element method for second order boundary value problems revealed that the former is more accurate than the latter. As a matter of fact, as far as the authors can see, the new method's only real demerit is the fact that non symmetric linear systems have to be solved, even when the problem at hand is self-adjoint. The primary aim of this paper is to show that, in such a case, an efficient symmetrization of the solution procedure can be achieved by means of a fast converging iterative method. So far the nonparametric approach considered in this work was only studied as applied to finite elements, which are conforming in the case of polytopic domains. However our technique to handle Dirichlet conditions prescribed on curved boundaries has a wide scope of applicability. This feature is exemplified here by applying such a symmetrization procedure to the solution of the three-dimensional Poisson equation by a nonconforming quadratic finite element with degrees of freedom other than nodal values. This method is based on the same type of piecewise quadratic interpolation as the one introduced in [START_REF] Ruas | Finite element solution of 3D viscous flow problems using non standard degrees of freedom[END_REF], in order to represent the velocity in the framework of the stable solution of incompressible viscous flow problems. Actually the corresponding velocity representation enriched by the quartic bubble-functions of the tetrahedra combined with a discontinuous piecewise linear pressure in each tetrahedron, is a sort of nonconforming three-dimensional analog of the popular conforming Crouzeix-Raviart mixed finite element [START_REF] Crouzeix | Conforming and nonconforming finite element methods for solving the stationary Stokes equations I[END_REF] for solving viscous flow problems in two-dimension space. After carrying out a numerical validation of the Petrov-Galerkin approach combined with the symmetrization procedure for a nonparametric version of this nonconforming quadratic method, its efficiency as compared to the isoparametric and nonparametric versions of the conforming quadratic element is examined. An error estimation in the broken H 1 -norm for the nonconforming method in the case of a curved domain completes these studies.

An outline of the paper is as follows. Section 2 is devoted to some preliminaries, in which we first recall the model Poisson equation in a smooth three-dimensional domain and present some pertaining notations; several definitions, notations and assumptions related to the finite element meshes are also introduced therein. In Section 3 we describe our technique to handle the Dirichlet boundary conditions for the model problem, in connection with the nonconforming quadratic finite element method; the underlying approximate problem is posed and corresponding stability and well-posedness results are given. In Section 4 we address the symmetrization solution procedure and validate the resulting numerical scheme. In Section 5 the performance of such a scheme is compared with the one of the asymmetric solution procedure, both extended to the standard conforming quadratic Lagrange element. Error estimates for the nonconforming method in the Petrov-Galerkin formulation to treat curved boundaries are given in Section 6. Finally in Section 7 we draw some conclusions from the whole work.

Preliminaries

In this section we specify the model problem considered in this work and supply some material to be used in the sequel.

The model problem and pertaining notations

Let us consider as a model the Poisson equation with Dirichlet boundary conditions in a threedimensional domain Ω with boundary Γ having suitable regularity properties, that is,

{ -∆u = f in Ω u = g on Γ, (1)
where f and g are given functions defined in Ω and on Γ.

For quadratic finite element methods our technique is most effective in case u ∈ H 3 (Ω). In order to make sure that u possesses such a regularity property we shall assume that f ∈ H 1 (Ω) and g ∈ H 5/2 (Γ) (cf. [START_REF] Adams | Sobolev Spaces[END_REF]). We observe that, owing to the Sobolev Embedding Theorem [START_REF] Adams | Sobolev Spaces[END_REF], g is necessarily continuous. We must further assume that Γ is at least of the C 1 -class. Actually, more than this, we make the assumption that the principal curvatures of Γ (cf. [START_REF] Cartan | Formes différentielles[END_REF]) are uniquely defined almost everywhere. Notice that in doing so we are not requiring that Γ be of the C 2 -class.

Throughout this article ∥ • ∥ 0 stands for the standard norm of L 2 (Ω). Furthermore ∥ • ∥ r,D and |•| r,D represent, respectively, the standard norm and semi-norm of Sobolev space H r (D) (cf. [START_REF] Adams | Sobolev Spaces[END_REF]), for r ∈ ℜ + with H 0 (D) = L 2 (D), D being any bounded subset of ℜ 3 . We also denote by

∥ • ∥ m,p,D the usual norm of W m,p (D) for m ∈ IN * and p ∈ [1, ∞] \ {2} with W 0,p (D) = L p (D). Whenever D is Ω the subscript , D is dropped.

Meshes and related notions

Let us be given a mesh T h consisting of straight-edged tetrahedra 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 this mesh is assumed to fit Ω in such a way that all the vertexes of the polyhedron ∪ T ∈T h T lie on Γ. We denote the interior of this union set by Ω h and define Ωh := Ω ∩ Ω h together with Ω ′ h := Ω ∪ Ω h . The boundaries of Ω h and Ωh are respectively denoted by Γ h and Γh and moreover Γ ′ h := Ωh ∩ Γ. T h is assumed to belong to a regular family of partitions in the sense of [START_REF] Ciarlet | The Finite Element Method for Elliptic Problems[END_REF], though not necessarily quasi-uniform. The boundary of every ∀T ∈ T h is represented by ∂T , while h T is the diameter of T and h := max T ∈T h h T . We make the non essential and yet reasonable assumption that any element in T h have at most either one edge or one face contained in Γ h . Let S h be the subset of T h consisting of tetrahedra having one face on Γ h and R h be the subset of T h \ S h consisting of tetrahedra having exactly one edge on Γ h . We further set O h := S h ∪ R h . Notice that, owing to our initial assumption, the interior of any tetrahedron in T h \ O h has an empty intersection with Γ h . For every T ∈ S h we denote by O T the vertex of T not belonging to Γ. Finally we introduce the notations ∥ • ∥ 0,h (resp. ∥ • ∥ 0,h ) for the standard norm of L 2 (Ω h ) (resp. L 2 ( Ωh )).

Remark 1 Even though for practical purposes this is by no means necessary, in all the constructions and analyzes given hereafter, we shall assume that the mesh is sufficiently fine. We refer to [START_REF] Ruas | Optimal-rate finite-element solution of Dirichlet problems in curved domains with straight-edged tetrahedra[END_REF] for a precise quantification of the assumed smallness of h.

We also need some definitions and auxiliary results regarding the set

(Ω \ Ω h ) ∪ (Ω h \ Ω).
With every edge e of the mesh contained in Γ h we associate a closed plane set δ e containing e, delimited by Γ and e itself. The plane of δ e can be arbitrarily chosen about e. However for better results it should be close to the bisector of the faces of the pair of elements in S h intersecting at e, which can eventually be a face shared by both. Such a choice will be assumed throughout this work. We also define δe := δ e ∩ Ω. In Figure 1 we illustrate one out of three plane sets δ e corresponding to the edges of the faces F T and F T ′ contained in Γ h of tetrahedra T and T ′ belonging to S h . More precisely δ e is depicted for the edge e common to F T and F T ′ .

Further, for every T ∈ S h , we define a closed set ∆ T delimited by Γ, ∂T and the plane sets δe associated with the edges of F T , as illustrated in Figure 1. In this manner we can assert that, if Ω is convex, Ω h is a proper subset of Ω and Ω is the union of the disjoint sets Ω h and ∪ T ∈S h ∆ T .

Otherwise Ω h \ Ω is a nonempty set containing subsets of T ∈ S h whose volume is an O(h 4 T ) and subsets of T ∈ R h whose volume is an O(h 5 T ), both types of subsets corresponding to non-convex portions of Γ. Whatever the case, the above configurations are of merely academic interest and carry no practical meaning, as much as the sets

T ∆ := T ∪∆ T ∀T ∈ S h or T ∆ := T ∪δ e ∀T ∈ R h , T := T ∩ Ω ∀T ∈ O h and ∆ ′ T := ∆ T \ Ω.
Referring to Figures 2 and3 for illustrations in particular cases, T h is supposed to fulfill the following reasonable conditions: Assumption + : h is small enough for the intersection P with Γ of the half line s with origin at O T passing through any point M ∈ F T to be uniquely defined ∀T ∈ S h .

Assumption ++ : h is small enough for the intersection Q ∈ δ e with Γ of the half line r perpendicular to e with origin at any point N ∈ e to be uniquely defined.

We recall a result formally established in [START_REF] Ruas | Optimal-rate finite-element solution of Dirichlet problems in curved domains with straight-edged tetrahedra[END_REF], according to which there exists a mesh-

independent constant C Γ such that length(M P ) ≤ C Γ h 2 T and length(N Q) ≤ C Γ h 2 T .
3 A nonconforming method with mean-value degrees of freedom

In this section we apply our technique to handle Dirichlet conditions on curved boundaries to a nonconforming method with degrees of freedom other than function nodal values. Incidentally we note that for many well known nonconforming finite element methods the construction of an isoparametric counterpart brings no improvement. This does not prevent suitable parametric elements from being successfully employed in this case. However to the best of author's knowledge studies in this direction are incipient. This fact motivates us to show here that our technique for handling curvilinear boundaries can be optimally extended in a straightforward manner to finite element methods, which are nonconforming even in the case of polytopes. We

T Є S h = A 1 A 2 A 3 A 4 A 1 e = A 1 A 2 A 3 A 6 A 2 A 4 A 5 T' Є S h = A 1 A 2 A 5 A 6 T'' Є R h = A 1 A 2 A 4 A 5 δ e Γ Γ h Δ T ' Δ T Figure 1: Sets ∆ T , ∆ T ′
, δ e for T, T ′ ∈ S h having a common edge e and a tetrahedron T ′′ in R h use such a nonconforming approach to solve the model problem ( 1), confining ourselves to the case of homogeneous boundary conditions for the sake of simplicity, though without any loss of essential aspects.

To begin with we recall the space V h of test-functions defined in Ω h , associated with the method under consideration. F and e being a face and an edge of a tetrahedron T ∈ T h respectively, we denote by M the centroid of F , by A and B the end-points of e and by N the mid-point of e. Now any function v ∈ V h restricted to every T is a polynomial of degree less than or equal to two, defined upon the following set of degrees degrees of freedom:

• The four values µ F (v) of v at the centroids M of F ; • The six mean values ν e (v) along e, where ν e (v) = 0.4v(N ) + 0.3[v(A) + v(B)].
∀v ∈ V h and ∀F and e, we require that both µ F (v) and ν e (v) coincide for all tetrahedra of the mesh sharing the face F or the edge e; moreover we require that both µ F (v) and ν e (v) vanish whenever F or e is contained in Γ h . Clearly enough these requirements are not sufficient to ensure the continuity in Ω h of a function in V h , and hence this space is not a subspace of H 1 0 (Ω h ). The set of local canonical quadratic basis functions in a tetrahedron T ∈ T h associated with the above degrees of freedom can be found in [START_REF] Ruas | Finite element solution of 3D viscous flow problems using non standard degrees of freedom[END_REF]. It is noteworthy that the gradients of all of them are an O(h -1 T ). This is a key property for the proof of Lemma 3.1 hereafter.

Similarly to the case of the standard Lagrangian piecewise quadratic elements, we define the trial-function space W h in the same way as V h , except for the fact that the degrees of freedom associated with faces F and edges e contained in Γ h are modified as follows: For a given function w ∈ W h , µ F (w) is replaced by µ ′ F (w) defined to be the value of w at the point P lying in the nearest intersection with Γ of the perpendicular to F passing through the centroid M of F as depicted in Figure 2; referring to Figure 3, ν e (w) is replaced by e , e ⊂ Γ h . If Q were replaced by the mid-point N of e, it is clear that the result would hold true, according to the properties of the interpolation under consideration (cf. [START_REF] Ruas | Finite element solution of 3D viscous flow problems using non standard degrees of freedom[END_REF]). The vector ⃗ a of coefficients a i for i = 1, 2, . . . , 10 of the underlying canonical basis functions φ j ∈ P 2 (T ) for 1 ≤ i ≤ 10 would be precisely b i for 1 ≤ i ≤ 10. Denoting the associated degrees of freedom (for a straight-edged tetrahedron) by π i , where π i is some µ F for 1 ≤ i ≤ 4 and π i is some ν e for 5 ≤ i ≤ 10, we assume that the corresponding canonical basis functions φ i are numbered accordingly. This means that the matrix K whose entries are

ν ′ e (w) := 0.4w(Q) + 0.3[w(A) + w(B)],
k ij := π i (φ j ) is the identity matrix. Let π ′ i = π i if 1 ≤ i ≤ 9
and π ′ 10 (w) be given by 0.4w(Q) + 0.3[w(A) + w(B)] for any w ∈ P 2 (T ), where A and B are the end-points of e ⊂ Γ h . We must establish that the 10 × 10 linear system of algebraic equations K ′ ⃗ a = ⃗ b is uniquely solvable, where K ′ is the matrix with entries k

′ ij := π ′ i (φ j ). Clearly we have K ′ = K + E K , where the entries of E K are e ij := π ′ i (φ j ) -π i (φ j ).
At this point we recall the [START_REF] Ruas | Optimal-rate finite-element solution of Dirichlet problems in curved domains with straight-edged tetrahedra[END_REF] we know that ∥ grad φ j ∥ 0,∞,T ∆ ≤ C ∞ ∥ grad φ j ∥ 0,∞,T for a suitable mesh-independent constant C ∞ . Moreover from standard arguments we know that the latter norm in turn is bounded above by another mesh-independent constant times h -1 T . In short we have

C Γ h 2 T . From Rolle's Theorem it follows that ∀ i, j, |e ij | ≤ 0.4C Γ h 2 T ∥ grad φ j ∥ 0,∞,T ∆ . From Lemma 3.2 of
|e ij | ≤ C E h T ∀ i, j,
where C E is independent of T . Hence the matrix K ′ equals the identity matrix plus an O(h T ) matrix E K . Therefore K ′ is an invertible matrix, as long as h is sufficiently small. The case of an element T ∈ S h can be dealt with as a mere variant of the above argument, and in this respect we also refer to Lemma 3.3 of [START_REF] Ruas | Optimal-rate finite-element solution of Dirichlet problems in curved domains with straight-edged tetrahedra[END_REF]. Lemma 3.1 allows us to assert that W h is indeed a nonempty function space, whose dimension equals the one of V h . Before pursuing we introduce the broken gradient operator grad h for any function w defined in Ω h which is continuously differentiable in every T ∈ T h , given by

[grad h w] |T ≡ grad w |T ∀T ∈ T h . Now if u is a function in H 2 (Ω) ∩ H 1 0 (Ω), we can define I h (u) ∈ W h to be the function given by µ F (I h (u)) = µ F (u)
and ν e (I h (u)) = ν e (u) for all the faces F and edges e of tetrahedra in T h not contained in Γ h . From standard interpolation results it is not difficult to establish that I h enjoys the following property: There exists a mesh-independent constant

C P such that ∀u ∈ H 3 (Ω) ∩ H 1 0 (Ω) it holds, ∥ grad h (u -I h (u)) ∥ 0,h ≤ C P h 2 |u| 3 . ( 2
)
Extending f by zero in Ω h \ Ω and still denoting the resulting function by f , the following problem is considered to approximate (1):

                     Find u h ∈ W h such that a h (u h , v) = L h (v) ∀v ∈ V h , where a h (w, v) := ∫ Ω h grad h w • grad h v, for w ∈ W h + H 1 (Ω h ), v ∈ V h and L h (v) := ∫ Ω h f v ∀v ∈ V h .
(

The matrix associated with (3) is a sparse band matrix whose sparsity structure is the same as for the standard Galerkin FEM, in which the spaces of trial functions and test functions coincide. However here such a matrix is non symmetric, since the basis functions of W h and V h are the same only for nodes not belonging to elements in O h . Hence the stability and well-posedness of problem (3) are not trivial issues, which we next address.

Proposition 3.2 If h is sufficiently small there exists a constant α > 0 independent of h such that,

∀w ∈ W h ̸ = 0, sup v∈V h \{0} a h (w, v) ∥ grad h w ∥ 0,h ∥ grad h v ∥ 0,h ≥ α. ( 4
)
Proof. Given w ∈ W h , let v be the unique function in V h such that all its degrees of freedom attached to a face or an edge of the mesh not contained in Γ h coincide with those of w. Notice that by construction µ F (v) = 0 and ν e (v) = 0 as long as F or e is contained in Γ h . For a given T ∈ O h we denote by m T the number of degrees of freedom {π T i } m T i=1 of V h attached to a face F or an edge e contained in Γ h . Clearly enough we have

a h (w, v) = ∑ T ∈T h ∫ T |grad w| 2 - ∑ T ∈O h ∫ T grad w • grad r T (w), (5)
where r T (w) = ∑ m T i=1 π T i (w)φ T i , φ T i being the canonical basis function of the space P 2 (T ) associated with the degree of freedom π T i . Now from standard results it holds ∥ grad φ

T i ∥ 0,T ≤ C φ h 1/2 T
where C φ is a mesh independent constant. Referring to Figures 2 and3, since

w(P ) = µ ′ F (w) = 0 (resp. 0.4w(Q) + 0.3[w(A) + w(B)] = ν ′ (w) = 0)
, where F (resp. e) generically represent a face (resp. an edge) of T contained in Γ h , in accordance with the definition of W h , a simple Taylor expansion about P (resp. Q) allows us to conclude that |w(M )| (resp. |w(N )|) are bounded above by l ∥ grad w ∥ 0,∞,T ∆ , where l = length(P M ) (resp. length(QN )), or yet that |w(M )| (resp.

|w(N )|) is bounded above by C Γ h 2 T ∥ grad w ∥ 0,∞,T ∆ . On the other hand from Lemma 2.2 of [18] it holds ∥ grad w ∥ 0,∞,T ∆ ≤ C J h -3/2 T
∥ grad w ∥ 0,T for a mesh-independent constant C J . Plugging all those estimates into (5), since m T ≤ 4, we obtain:

a h (w, v) ≥ ∫ Ω h |grad h w| 2 -4C φ C J C Γ h ∑ T ∈O h ∥ grad w ∥ 2 0,T . ( 6
)
Then it holds with

c := 4C φ C J C Γ , ( 7
)
a h (w, v) ≥ (1 -ch) ∥ grad h w ∥ 2 0,h . ( 8
)
Now using arguments in all similar to those employed above, we easily conclude that

∥ grad h v ∥ 0,h ≤∥ grad h w ∥ 0,h + ∥ grad h v -grad h w ∥ 0,h ≤ (1 + ch) ∥ grad h w ∥ 0,h . (9)
Combining ( 8) and ( 9), provided h ≤ (2c) -1 we establish (4) with α = 1/3. Proposition 3.3 Provided h is sufficiently small, problem (3) has a unique solution.

Proof. From well-known results (cf. [START_REF] Babuška | The finite element method with Lagrange multipliers[END_REF], [START_REF] Brezzi | On the existence, uniqueness and approximation of saddle-point problems arising from Lagrange multipliers[END_REF] and [START_REF] Cuminato | Unification of distance inequalities for linear variational problems[END_REF]) this is an immediate consequence of Proposition 3.2 and of the fact that V h and W h have the same dimension.

Symmetrization of the solution procedure

Since (3) is not a symmetric problem we can use the following iterative procedure to solve it as a sequence of symmetric problems. First of all let n h be the dimension of both V h and W h , that is the total number of degrees of freedom of both spaces not assigned to zero beforehand. Let also ∥ • ∥ 0,∞,h be the norm of either V h or W h defined to be the maximum absolute value of their n h degrees of freedom. For every v ∈ V h we denote by Π W (v) the function of W h whose degrees of freedom coincide with those of v. Similarly for every w ∈ W h we denote by Π V (w) the function of V h whose degrees of freedom coincide with those of w. Now we consider the following symmetric problem,

{ Find ū0 h ∈ V h such that a h (ū 0 h , v) = L h (v) ∀v ∈ V h , (10)
which is clearly uniquely solvable. Defining

u 0 h := Π W (ū 0 h ) ∈ W h , (11)
solve successively for n = 1, 2, . . . the problems,

   Find u n h ∈ W h := Π W (ū n h ) where ūn h ∈ V h is the unique solution of a h (ū n h , v) = a h (ū n-1 h , v) -a h (u n-1 h , v) + L h (v) ∀v ∈ V h , (12)
until ∥ u n h -u n-1 h ∥ 0,∞,h is less than a small tolerance ε. Since the matrix associated with ( 12) is a symmetric positive definite matrix for the standard Galerkin FEM, the stability and well-posedness of ( 12) is guaranteed. It is also a band matrix with the same sparsity structure within its band as the matrix associated with [START_REF] Brenner | The Mathematical Theory of Finite Element Methods[END_REF]. Let us study the convergence of the above iterative procedure. With this aim we first set ūh = Π V (u h ) and note that, 12) with ( 13) we have:

a h (ū h , v) = a h (ū h , v) -a h (u h , v) + L h (v) ∀v ∈ V h . ( 13
) ∀n ≥ 0, let w n h := u n h -u h ∈ W h and wn := ūn h -ūh ∈ V h . Combining (
a h ( wn h , v) = a h ( wn-1 h , v) -a h (w n-1 h , v) ∀v ∈ V h . ( 14
)
We next establish that, provided h is sufficiently small, ∥ grad h w n h ∥ 0,h tends to zero roughly as fast as b(h)(2ch) n as n goes to infinity, where c fulfills 2ch ≤ 1 and b(h) is an O(h).

Since wn-1

h only differs from w n-1 h in elements in O h we have, a h ( wn-1 h , v) -a h (w n-1 h , v) = ∑ T ∈O h ∫ T grad( wn-1 h -w n-1 h ) • grad v. ( 15
)
Using the same arguments leading to [START_REF] Ciarlet | The Finite Element Method for Elliptic Problems[END_REF] together with [START_REF] Ciarlet | The Mathematical Foundations of the Finite Element Method with Applications to Partial Differential Equations[END_REF], we obtain successively,

a h ( wn-1 h , v) -a h (w n-1 h , v) ≤ ch ∑ T ∈O h ∥ grad w n-1 h ∥ 0,T ∥ grad v ∥ 0,T , ( 16
)
a h ( wn-1 h , v) -a h (w n-1 h , v) ≤ ch ∥ grad h w n-1 h ∥ 0,h ∥ grad h v ∥ 0,h . ( 17
)
Taking v = wn h in both ( 14) and ( 17) we come up with,

∥ grad h wn h ∥ 0,h ≤ ch ∥ grad h w n-1 h ∥ 0,h (18)
Now noting that ∀m ≥ 0 ∥ grad h w m h ∥ 0,h ≤∥ grad h wm h ∥ 0,h + ∥ grad h (w m h -wm h ) ∥ 0,h , similarly to [START_REF] Ruas | Optimal Dirichlet-condition enforcement on curved boundaries for Lagrange and Hermite FEM with straight-edged simplexes[END_REF], as long as h is less than 1/c, we easily conclude that

∥ grad h w m h ∥ 0,h ≤ (1 -ch) -1 ∥ grad h wm h ∥ 0,h ∀m ≥ 0. (19)
Plugging [START_REF] Scott | Finite Element Techniques for Curved Boundaries[END_REF] with m = n -1 into (18) we establish that,

∥ grad h wn h ∥ 0,h ≤ ρ(h) ∥ grad h wn-1 h ∥ 0,h ∀n > 1 with ρ(h) = ch/(1 -ch). ( 20
)
Assuming that h < 1/(2c) the fraction ρ(h) will be less than one, and hence the quantity ∥ grad h wn h ∥ 0,h will decrease by a factor of ρ(h) at every iteration. Actually using again [START_REF] Scott | Finite Element Techniques for Curved Boundaries[END_REF], this time with m = n, and noting that 1 -ch ≥ 1/2 by assumption, we have,

∥ grad h w n h ∥ 0,h ≤ σ(h)(2ch) n with σ(h) :=∥ grad h (ū 0 h -ūh ) ∥ 0,h /2 ∀n > 1. (21)
Observing that ū0

h -ūh = (ū 0 h -u) + (u -u h ) + (u h -ūh )
and that the orders of magnitude of the norms ∥ grad h • ∥ 0,h of the terms in parentheses on the right hand side are respectively O(h 3/2 ), O(h 2 ) and O(h), we can assert that σ(h) is bounded above by a coefficient b(h), whose order of magnitude is at most an O(h). All this advocates in favor of a faster convergence of the iterations [START_REF] Nitsche | On Dirichlet problems using subspaces with nearly zero boundary conditions[END_REF], the smaller h.

Let us check the efficiency of the iterative symmetrization procedure ( 10)-( 11)-( 12) by solving a test-problem with successively refined meshes. The solution of the linear system resulting from ( 12) is computed by means of both Cholesky's method with BMS (band matrix storage) and the CG (conjugate gradient) method by storing only the non zero coefficients of the matrix, i.e., with VSMS (very sparse matrix storage). In the model problem Ω is the ellipsoid of equation x 2 /a 2 + y 2 /b 2 + z 2 < 1 in a cartesian coordinate system (x, y, z), whose origin is its center, with a = 0.6 and b = 0.8. We take an exact solution given by u(x, y, z) = (1 -

x 2 /a 2 -y 2 /b 2 - z 2 )(1 -x 2 /b 2 -y 2 /a 2 -z 2 )
, so that f = -∆u. The computations are carried out only for the octant corresponding to non negative values of the coordinates, with a family of quasi-uniform meshes consisting of 6p 3 tetrahedra for an integer p ≥ 1. For each value of p the mesh of the ellipsoid is the transformation of the uniform mesh of a unit cube with 6p 3 tetrahedra having edges parallel to the line x = y = z, by suitably mapping the set of vertexes of the latter given in cartesian coordinates into the one of the actual mesh expressed in spherical coordinates. In this manner we have h ≃ p -1 . In Tables 1 and2 we show the number of iterations m necessary to satisfy tolerances of ε = 10 -5 and ε = 10 -7 , for increasing values of p, using the Cholesky and the CG solver, respectively. The smaller value of ε in the latter case is due to the observation that this tolerance must be compatible with the necessarily small one in the convergence test of the CG method.

∥ u m h -u m-1 h ∥ 0,∞,h -→ 0.62431E-5 0.51154E-5 0.56915E-5 0.12739E-5 0.92735E-5 Table 1: Number of iterations m such that ∥ u m h -u m-1 h
∥ 0,∞,h < 10 -5 using Cholesky's method

According to Table 1 the number of iterations necessary for convergence of the symmetrization procedure decreases indeed with the mesh size. Table 2 in turn points in the opposite direction, but this effect can be credited to the combination of two iterative procedures. Nevertheless such a behavior is far from being a drawback, as seen in the next section.

∥ u m h -u m-1 h ∥ 0,∞,h -→ 0.14937E-7 0.42120E-7 0.24374E-7 0.93131E-7 0.65410E-7 Table 2: Number of iterations m such that ∥ u m h -u m-1 h
∥ 0,∞,h < 10 -7 using the CG method

Comparative study

In this section we further investigate the iterative solution procedure of symmetric problems posed in the nonparametric Petrov-Galerkin variational form of the type [START_REF] Nitsche | On Dirichlet problems using subspaces with nearly zero boundary conditions[END_REF]. More particularly we carry out a comparative study thereof with the direct solution of the underlying non symmetric linear system. In this framework two approaches are assessed: The direct solution of the non symmetric system performed by either Crout's method with partial pivoting and BMS or the GMRES method with VSMS; the iterative solution with a symmetric positive definite matrix performed by either Cholesky's method with BMS or the CG method with VSMS, resp. Additionally we extend such numerical comparisons to the classical conforming quadratic finite element in both the nonparametric Petrov-Galerkin form and the isoparametric version. A Lenovo T440s laptop was employed in all the computations reported below.

Iterative vs. direct solution of (3)

First of all we compare the performance of the nonconforming quadratic method studied in Sections 3 and 4 to approximate (1) using both solution strategies. With this aim we take the same test-problem as in the previous section. Depending on whether the iterative symmetrization procedure is employed or not, in this comparison we use both direct solvers with BMS, namely, Cholesky's method and Crout's method, and the iterative solvers GC and GMRES with VSMS. We supply in Table 3 the total processing (CPU) time in seconds for successively refined meshes, using direct solvers for both the iterative symmetrization procedure with a tolerance equal to 10 -5 and the direct solution. Similarly, we display in Table 4 the total CPU time in seconds for successively refined meshes for both the scheme [START_REF] Nitsche | On Dirichlet problems using subspaces with nearly zero boundary conditions[END_REF] and the direct solution, using now the iterative solvers and a tolerance equal to 10 -7 for all iterative procedures. From the above results we infer the overwhelming superiority of the iterative symmetrization strategy over the direct solution, since the ratio of CPU times for the former and the latter decreases as the mesh is refined. This effect is significantly magnified in case direct solvers are used, in which unreasonable processing times for barely intermediate meshes are either reported or omitted (cf. Table 3). On the other hand, the fact that an increasing number of iterations is necessary for convergence of the procedure [START_REF] Nitsche | On Dirichlet problems using subspaces with nearly zero boundary conditions[END_REF] as the mesh is refined (cf. Table 2), is probably the cause of a lesser discrepancy of CPU times in case iterative solvers are used, as one can see in Table 4. Notice however that in spite of these observations, the drastic reduction of matrix storage for iterative solution methods advocates in favor of them, as compared to direct ones.

Additional comparisons involving second order finite-element methods

Keeping the same test-problem as above, we pursue the performance evaluation of our iterative scheme as compared to the direct solution, taking also classical conforming quadratic finite elements. In order to have better insight on the merits of the nonparametric Petrov-Galerkin formulation, such comparisons are extended to isoparametric finite elements of the same order in standard Galerkin formulation.

Referring to [START_REF] Ruas | Optimal-rate finite-element solution of Dirichlet problems in curved domains with straight-edged tetrahedra[END_REF], let ǔh represent the approximate solution of (1) related to mesh T h obtained by the nonparametric Petrov-Galerkin approach, in connection with conforming Lagrange quadratic finite elements.

To begin with we illustrate the strength of the nonconforming approach, by displaying in Tables 5 through 8 data related to u h and ǔh , respectively, for different values of p, that is h. Besides the errors measured in three different manners, we supply a degree of freedom (DOF) count for both FEMs. The number of iterations necessary to satisfy the stop criterion for the iterative scheme of the type ( 12) is still denoted by m for the nonconforming method and by m for the conforming method. The corresponding total number of DOFs are denoted by M and M respectively. Similarly to the previous subsection in the stop criterion the tolerance ε applies to the maximum absolute value of the difference between DOFs in two successive iterations.

The results given in Tables 5 and6 were obtained with a Cholesky solver for ε = 10 -5 . As one infers from Tables 5 and6 the methods under experimentation are both of the third order in L 2 (Ω h ) and of the second order in the broken (semi)norm of H 1 (Ω h ) as expected or predicted either in [START_REF] Ruas | Optimal-rate finite-element solution of Dirichlet problems in curved domains with straight-edged tetrahedra[END_REF] or in Section 6 hereafter. Both methods are also fairly equivalent from the point of view of accuracy in these norms. On the other hand there is a clear advantage of the nonconforming method over the conforming method in terms of DOF ("pointwise") errors.

p -→ 2 4 8 16
∥ u -u h ∥ 0,h -→ 0.64013E-2 0.88793E-3 0.11467E-3 0.14585E-4 ∥ grad h (u -u h ) ∥ 0,h -→ 0.11549E+0 0.34444E-1 0.91569E-2 0.23477E-2 ∥ u -u h ∥ 0,∞,h -→ 0.27816E-1 0.39244E-2 0.62257E-3 0.85775E-4 m -→ 7 5 4 3 M -→ 218 1,
∥ u -ǔh ∥ 0,h -→ 0.70568E-2 0.95648E-3 0.12203E-3 0.15445E-4 ∥ grad h (u -ǔh ) ∥ 0,h -→ 0.11772E+0 0.353106E-1 0.94375E-2 0.24253E-2 ∥ u -ǔh ∥ 0,∞,h -→ 0.36064E-1 0.69394E-2 0.10616E-2 0.14339E-3 m -→ 6 5 4
The results in Tables 7 and8 were obtained by using a CG solver with VSMS. Here we took ε = 10 -7 , which is also the tolerance employed in the stop criterion for the CG method. 7 and8 confirm roughly the same orders of both FEM observed in Tables 5 and6, and the slightly better accuracy of the nonconforming method except for the L 2 -norm of the error for the finest mesh. Notice that the number of iterations necessary for convergence of the conforming method decreases smoothly as the mesh is refined, as expected, in contrast to the nonconforming method. This could explain the more significant deterioration of the accuracy in the L 2 -norm observed for the latter method, as compared to the former. Whatever the case, such an effect advocates in favor of direct solvers, since in this case there is no need to adjust a tolerance to optimally fit the one of the iterative symmetrization scheme itself. However it turns out that iterative solvers are in principle less time consuming for a given mesh, while requiring much less storage.

p -→ 3 6 12 24 ∥ u -u h ∥ 0,h -→ 0.20505E-2 0.26893E-3 0.34373E-4 0.56906E-5 ∥ grad h (u -u h ) ∥ 0,h -→ 0.58245E-1 0.15976E-1 0.41403E-2 0.10525E-2 ∥ u -u h ∥ 0,∞,h -→ 0.82123E-2 0.13700E-2 0.19521E-3 0.41366E-4 m -→
∥ u -ǔh ∥ 0,h -→ 0.22262E-2 0.28733E-3 0.36439E-4 0.48589E-5 ∥ grad h (u -ǔh ) ∥ 0,h -→ 0.59516E-1 0.16437E-1 0.42741E-2 0.10872E-2 ∥ u -ǔh ∥ 0,∞,h -→ 0.13889E-1 0.23689E-2 0.33214E-3 0.43462E-4 m -→ 7 6 6
It is also interesting to watch the behavior of both methods in terms of CPU time, when the direct and the iterative solving approaches are employed. Tables 9 and10 supply the CPU times for the conforming quadratic method with successively refined meshes, similarly to Tables 3 and4 respectively, for the nonconforming method.

It is noticeable here again the great superiority of the iterative approach from the point of view of processing time. A quick comparison of Tables 3 and4 with Tables 9 and 10 also indicates that the nonconforming method is much more time consuming than the conforming method. However this is no surprise since there are more than twice as many degrees of freedom for the latter with respect to the former for the same mesh. Next we compare the nonparametric Petrov-Galerkin formulation for the conforming quadratic element with the corresponding isoparametric formulation, whose optimal second order in the H 1 -norm was established in [START_REF] Ciarlet | The Mathematical Foundations of the Finite Element Method with Applications to Partial Differential Equations[END_REF]. We denote by ũh the approximate solution to (1) determined by the isoparametric technique for the same mesh as ǔh Naturally enough ũh is computed using the Cholesky's method with BMS and the CG method with VSMS. However for a more fair comparison with the nonparametric approach in terms CPU time, we also compute the isoparametric solution using Crout's method and the GMRES method, without taking into account symmetry. This is because in the case of non symmetric problems the use of both Cholesky's method and the CG method has to be discarded.

CPU times necessary to determine ũh with direct solvers are displayed in Table 11 11 and12 confirm the great superiority in terms of CPU time, of methods whose use is restricted to symmetric positive definite matrices, over methods applying to any regular matrix. In particular Cholesky's method is much better than Crout's method as shown in Table 11. Moreover, resorting to Tables 9 and11, it turns out that both approaches are fairly equivalent in terms of CPU, with a slight advantage of isoparametric elements over nonparametric elements. This contradicts observations in the opposite sense in the two-dimensional case (cf. [START_REF] Ruas | Optimal simplex finite-element approximations of arbitrary order in curved domains circumventing the isoparametric technique[END_REF]). On the other hand, if one compares the solutions using Cholesky's method, isoparametric elements perform a little better only for the coarser meshes, while the contrary occurs in an increasingly significant manner as the mesh is refined. Such a behavior is noteworthy taking into account that iterations are necessary for the nonparametric approach, in contrast to the isoparametric approach. This seems to advocate in favor of the former, and could be due to its better matrix conditioning. We push further our numerical study by comparing the solutions determined by the nonparametric and the isoparametric approaches in terms of accuracy. In Table 13 the errors for the isoparametric solution computed by Cholesky's method are given in three different measures.

p -→ 2 4 8 16
∥ u -ũh ∥ 0,h -→ 0.75220E-2 0.10564E-2 0.13173E-3 0.16185E-4 ∥ grad h (u -ũh ) ∥ 0,h -→ 0.13931E+0 0.39089E-1 0.10015E-1 0.25061E-2 ∥ u -ũh ∥ 0,∞,h -→ 0.40980E-1 0.79148E-2 0.12384E-2 0.16897E-3
Table 13: Errors for the conforming quadratic FEM in isoparametric formulation

Comparing the results displayed in Tables 6 and13, we figure out that the nonparametric approach is a little more accurate than the isoparametric approach in all respects. Taking into account the previous observations, together with the two-dimensional experiments reported in [START_REF] Ruas | Optimal simplex finite-element approximations of arbitrary order in curved domains circumventing the isoparametric technique[END_REF] we are inclined to conclude that the former is definitively superior to the latter.

To conclude we comment on the cost of storage in the experiments reported in this section. First we note that the DOFs were numbered in a standard sequential manner for uniform meshes of a cube. In doing so the number of unknowns (NU) for the nonconforming method and the conforming method are 19p 3 and 8p 3 , respectively. This also leads to band matrices for both methods, whose half band width (HBW) for large values of p is asymptotically equal to 19p 2 for the nonconforming method and to 8p 2 for the conforming method. It follows that the direct solvers handle arrays whose total number of entries (TNE) are roughly 19 2 p 5 and 8 2 p 5 , respectively. This explains the growing discrepancy in CPU time to run direct solvers for both methods with the same mesh, as p increases (cf. Tables 7 and8). On the other hand, in case iterative solvers are used, arrays with TNE asymptotically equal to 19λ N C p 3 and 8λ C p 3 are handled for the nonconforming method and the conforming method, where λ N C ≃ 11 and λ C ≃ 14 and the subscripts NC and C stand for nonconforming and conforming. This is the reason why the ratios between CPU times to run both methods with the same mesh using iterative solvers are smaller, as shown in Tables 9 and10. Just to give an overview of the matrix storage required to run both finite element methods, we supply in self-explanatory Tables 14 and15 the above key figures as p varies, for direct and iterative solvers, respectively.

                          |a h (u ′ , v) -L h (v)| = c h (u ′ , v) + d h (u ′ , v)
where

c h (u ′ , v) = ∑ T ∈T h ∫ ∂T v ∂u ′ ∂n T and d h (u ′ , v) = - ∑ T ∈Q h ∫ ∆ ′ T ∆u ′ v. ( 27
)
c h (u ′ , v) can be estimated by means of standard arguments for nonconforming finite elements. More specifically in the case under study (cf. [START_REF] Ruas | Finite element solution of 3D viscous flow problems using non standard degrees of freedom[END_REF]) an estimate of the same nature as (24) applies to c h , i.e.,

c h (u ′ , v) ≤ C R h 2 |u ′ | 3,Ω h ∥ grad h v ∥ 0,h . ( 28
)
As for bilinear form d h first we observe that,

d h (u ′ , v) ≤ ∑ T ∈Q h [volume(∆ ′ T )] 1/2 ∥ ∆u ′ ∥ 0,∆ ′ T ∥ v ∥ 0,∞,∆ ′ T . ( 29
)
Since µ F (v) = 0 for all faces F contained in Γ h , there exists a mesh-independent constant

C ′ Γ such that ∥ v ∥ 0,∞,∆ ′ T ≤∥ v ∥ 0,∞,T ≤ C ′ Γ h T ∥ grad v ∥ 0,∞,T . (30)
Using the well-known inverse inequality (see e.g. [START_REF] Verfürth | A Posteriori Error Estimation Techniques for Finite Element Methods[END_REF]),

∥ w ∥ 0,∞,T ≤ C I h -3/2 T ∥ w ∥ 0,T ∀w ∈ P 2 (T ), (31)
where C I is a mesh-independent constant, like in Theorem 5.8 of [START_REF] Ruas | Optimal-rate finite-element solution of Dirichlet problems in curved domains with straight-edged tetrahedra[END_REF], the following result derives from (30),

∥ v ∥ 0,∞,∆ ′ T ≤ C ′ Γ C I h -1/2 T ∥ grad v ∥ 0,T . (32) Noticing that volume(∆ ′ T ) is bounded by h 4 T multiplied by a constant C Ω depending only on Ω, for both T ∈ S h ∩ Q h and T ∈ R h ∩ Q h , from straightforward calculations it follows that, ∥ ∆u ′ ∥ 0,∆ ′ T ≤ [C Ω ] 1/4 h T [ ∫ ∆ ′ T (∆u ′ ) 4 ] 1/4 ∀T ∈ Q h . ( 33
)
Then combining (29), (30), ( 32) and (33), applying the Cauchy-Schwarz inequality to the summation over T , and setting

C S := [C Ω ] 3/4 C ′ Γ C I we come up with, d h (u ′ , v) ≤ C S h 2    ∑ T ∈Q h h T [ ∫ ∆ ′ T (∆u ′ ) 4 ] 1/2    1/2 ∥ grad h v ∥ 0,h . ( 34
)
Applying againn the Cauchy-Schwarz inequality to the summation on the right hand side of (34) we readily obtain,

d h (u ′ , v) ≤ C S h 2   ∑ T ∈Q h h 2 T   1/4   ∑ T ∈Q h ∫ ∆ ′ T (∆u ′ ) 4   1/4 ∥ grad h v ∥ 0,h . ( 35
)
Noticing that there exists a constant ĈΓ such that

  ∑ T ∈Q h h 2 T   1/2 ≤ ĈΓ independently of h, (36)
we come up with,

d h (u ′ , v) ≤ C S [ ĈΓ ] 1/2 h 2 ∥ ∆u ′ ∥ 0,4,Ω h ∥ grad h v ∥ 0,h . ( 37
) Since H 1 (Ω ′ ) is continuously embedded in L 4 (Ω ′ ) (cf. [1]
), from (37) we infer the existence of a mesh-independent constant C R such that

d h (u ′ , v) ≤ C R h 2 ∥ ∆u ′ ∥ 1,Ω ′ ∥ grad h v ∥ 0,h , (38)
Now we plug (28) and ( 38) into (27), and then the resulting inequality into (26). Finally using the trivial variant of (2) according to which

∥ grad h (u ′ -I h (u ′ )) ∥ 0,h ≤ C ′ P h 2 |u ′ | 3,Ω ′ (39)
for a suitable mesh-independent constant C ′ P together with the triangle inequality, the result follows.

Conclusions

The authors believe to have undoubtedly demonstrated that the nonparametric Petrov-Galerkin formulation studied in this work is a very efficient universal tool to solve boundary value problems posed in curved domains with Dirichlet boundary conditions. This assertion is supported by several evidences presented throughout the article. The conclusions of the experimentation carried out in this work can be summarized as follows:

1. First of all we emphasize that, although the nonparametric formulation leads to non symmetric linear systems, even when the problem at hand is self-adjoint, in practical terms this fact is not a real demerit. Indeed, we saw that an easy-to-implement iterative procedure can be used to solve the system, thereby generating a fast-converging sequence of solutions of symmetric systems with a fixed matrix (to be factorized once for all before it starts, in the case of a direct solver). It turns out that this solution procedure is much less time consuming than the direct solution. Moreover we observed that it can perform better with respect to methods whose system matrix is symmetric anyway, such as the isoparametric formulation of self-adjoint problems.

2. Error estimates for the nonparametric formulation can be proved using the well established theory of linear variational problems. We should emphasize that this is not at all restricted to the nonconforming method studied in Section 6. Indeed a similar analysis applies to many other classes of methods, such as Lagrange FEM of any order higher than one, as shown in [START_REF] Ruas | Optimal Dirichlet-condition enforcement on curved boundaries for Lagrange and Hermite FEM with straight-edged simplexes[END_REF] and [START_REF] Ruas | Optimal-rate finite-element solution of Dirichlet problems in curved domains with straight-edged tetrahedra[END_REF], or yet Hermite FEM for biharmonic equations (cf. [START_REF] Ruas | Optimal Dirichlet-condition enforcement on curved boundaries for Lagrange and Hermite FEM with straight-edged simplexes[END_REF]).

3. The use of nonparametric shape and test functions allows for flexible constructions, in the sense that they are well adapted to several types of degrees of freedom, in contrast to classical formulations. In this work this property was exemplified more particularly for mean-value degrees of freedom associated with a nonconforming quadratic tetrahedral element, which adds to many other cases already addressed in [START_REF] Ruas | A simple alternative for accurate finite-element modeling in curved domains[END_REF], [START_REF] Ruas | A Hermite Method for Maxwell's Equations[END_REF], [START_REF] Ruas | Optimal-rate finite-element solution of Dirichlet problems in curved domains with straight-edged tetrahedra[END_REF] and [START_REF] Ruas | Optimal Dirichlet-condition enforcement on curved boundaries for Lagrange and Hermite FEM with straight-edged simplexes[END_REF].

4. The nonparametric Petrov-Galerkin formulation appeared to be more accurate than classical techniques for the same purpose, such as the isoparametric version of the finite element method, in case the latter exists.

Finally we note that some observations listed above had already been reported in the validation sections of previous publications such as [START_REF] Ruas | Optimal simplex finite-element approximations of arbitrary order in curved domains circumventing the isoparametric technique[END_REF], [START_REF] Ruas | A simple alternative for accurate finite-element modeling in curved domains[END_REF], [START_REF] Ruas | A Hermite Method for Maxwell's Equations[END_REF], [START_REF] Ruas | Optimal Dirichlet-condition enforcement on curved boundaries for Lagrange and Hermite FEM with straight-edged simplexes[END_REF] and [START_REF] Ruas | Optimal-rate finite-element solution of Dirichlet problems in curved domains with straight-edged tetrahedra[END_REF]. However here the authors focused on a systematic efficiency study of the nonparametric formulation. Nevertheless they are aware of the fact that more experimentation with this new technique is necessary, in order to evaluate it in contexts other than those considered in this article. For this reason they intend to push further this kind of study in future work.

Remark 2 Besides direct methods known for roughly one hundred years or more, the numerical experimentation in this work was carried out by means of two iterative methods widely in use to solve linear systems, namely, the conjugate gradient method and the GMRES method. As a by-product of our studies, the globally great superiority of iterative methods over direct methods was highlighted once more. This is particularly due to the fact that, in principle, the former are significantly less time consuming than the latter, while enabling practitioners to work with much finer meshes. In the authors' view, both advantages largely make up for the eventual need to adjust numerical parameters or to call on side techniques for improving the convergence and/or the accuracy of iterative methods. Among them lies preconditioning, but we declined to use this technique here in order to avoid deviation from our main validation and comparison goals. This is also because preconditioning may fail, depending on the kind of technique and the FEM in use, or yet bring about little improvement of performance, owing to a substantial increase of computational effort. But nothing prevents one from testing and comparing countless techniques for enhanced linear system solving, focusing on special situations. For example, it might be interesting to check the performance of the modification of the conjugate gradient algorithm proposed in [START_REF] Ehrel | An augmented conjugate gradient method for solving consecutive symmetric positive definite linear systems[END_REF] for consecutive linear systems. Eventually this technique could further reduce CPU time, in the framework of the iterative solution procedure of the type (12) experimented here, as long as the problem to solve is self-adjoint and positive definite.

Figure 2 :

 2 Figure 2: P := intersection with Γ of the line joining vertex O T to the centroid M of F T ⊂ Γ h

Figure 3 :

 3 Figure 3: Q := intersection with Γ ∩ δ e of the plane orthogonal to e containing its mid-point N

Table 3 :

 3 CPU time for solving (3): iterations (12) using Cholesky's method vs. Crout's method

	p	-→	4	6	8	12	16
	Iterative solution (Cholesky's method) -→ 4.36s	17.35s	420.78s	15,786.39s 92,061.36s
	Direct solution (Crout's method)	-→ 19.00s 136.92s 3,911.48s 230,274.28s	--

Table 4 :

 4

CPU time for solving (3): iterations (12) using the CG method vs. GMRES method

Table 5 :

 5 Errors, nr.of iterations & DOF count for the nonconforming FEM + Cholesky's method

	p	-→	2	4	8	16

Table 6 :

 6 Errors, nr. of iterations & DOF count for the conforming FEM + Cholesky's method

Table 7 :

 7 Errors, nr. of iterations & DOF count for the nonconforming FEM + the CG method

	p	-→	3	6	12	24

Table 8 :

 8 Errors, nr. of iterations & DOF count for the conforming FEM + the CG method Tables

Table 9 :

 9 CPU time for solving[START_REF] Adams | Sobolev Spaces[END_REF] with the conforming FEM via direct solvers

	p	-→	4	6	8	12	16
	Iterative solution (Cholesky's method) -→ 0.26s 1.64s 21.18s	235.38s	3,013.57s
	Direct solution (Crout's method)	-→ 0.51s 6.54s 212.45s 2,617.61s 40,226.10s
	p	-→	6	8	12	16	24
	Iterative solution (CG method)	-→ 1.09s 4.70s 38.77s 399.29s 2,189.46s
	Direct solution (GMRES method) -→ 2.76s 8.16s 78.18s 426.18s 4,234.10s

Table 10 :

 10 CPU time for solving[START_REF] Adams | Sobolev Spaces[END_REF] with the conforming FEM via iterative solvers

  . An iterative-solver counterpart in terms of CPU time is supplied in Table12.

	p	-→	4	6	8	12	16
	Cholesky's method	-→ 0.20s 1.69s 14.87s	289.66s	3,552.16s
	Crout's method (with ∂ pivoting) -→ 0.58s 6.55s 202.04s 2,248.96s 39,128.77s

Table 11 :

 11 CPU time for solving[START_REF] Adams | Sobolev Spaces[END_REF] with the isoparametric quadratic FEM via direct solvers

	p	-→	6	8	12	16	24
	CG method	-→ 0.99s 4.09s 42.96s 230.94s 2,560.18s
	GMRES method -→ 1.87s 8.96s 82.78s 500.15s 4,829.93s

Table 12 :

 12 CPU time for solving (1) with the isoparametric quadratic FEM via iterative solversIt is no surprise that Tables

Table 14 :

 14 Key storage data for the symmetric band matrices handled by Cholesky's method

	p	-→	3	6	12	24
	NC FE: TNE/NU (≃ λ N C ) -→ 4, 563/513 41, 526/4, 104 353, 268/32, 832 2, 912, 328/262, 656

C FE: TNE/NU (≃ λ C ) -→ 2, 125/216 21, 052/1, 728 186, 400/13, 824 1, 566, 856/110, 592

Table 15 :

 15 Key storage data for the symmetric sparse matrices handled by the CG method and ∂T for the boundary of T ∈ T h and denoting by ∂(•)/∂n T the normal derivative on ∂T oriented outwards T we obtain: 

Acknowledgment: The first author gratefully acknowledges the financial support provided by CNPq through grant 307996/2008-5. The authors are thankful to their colleague J. A. Cuminato for helpful discussions.

† This work was partially supported by CNPq, the National Research Council of Brazil 1

Error estimates

In this section we establish error estimates for problem [START_REF] Brenner | The Mathematical Theory of Finite Element Methods[END_REF]. Akin to [START_REF] Ruas | Optimal-rate finite-element solution of Dirichlet problems in curved domains with straight-edged tetrahedra[END_REF] we distinguish the convex case from the non-convex case. First we have: Theorem 6.1 Assume that f ∈ H 1 (Ω) and g ≡ 0. As long as h is sufficiently small, if Ω is a convex domain smooth enough for the solution u of (1) to belong to H 3 (Ω), there exists a constant C(f ) depending only on f such that the solution u h of (3) satisfies :

Proof. According to [START_REF] Cuminato | Unification of distance inequalities for linear variational problems[END_REF], using Proposition 3.2 we can write:

Proof. Taking into account ( 2), all we have to do is to estimate the sup term on the right hand side of (23). As a matter of fact such an issue was basically addressed in [START_REF] Ruas | Finite element solution of 3D viscous flow problems using non standard degrees of freedom[END_REF]. More precisely the required estimate is a consequence of the fact that the L 2 -projection of the trace on a face F of the mesh of any function v ∈ V h onto the space P 1 (F ), is a linear combination of the values µ F (v) and ν e (v), where e here generically represents the edges of F . This property implies the existence of a mesh-independent constant C R such that,

Then ( 22) directly follows from ( 23), ( 2) and ( 24).

Before pursuing we introduce Ω ′ as a smooth domain of ℜ 3 close to Ω but strictly containing both Ω and Ω h for all h small enough to conform to our assumptions on the meshes. According to Stein et al. [START_REF] Stein | Immersed boundary smooth extension: A high-order method for solving PDE on arbitrary smooth domains using Fourier spectral methods[END_REF] there exists an extension u ′ of u to Ω ′ such that u ′ ∈ H 3 (Ω ′ ) and u ′ ≡ u in Ω. Now we prove Theorem 6.2 Assume that u ∈ H 3 (Ω). Provided h is sufficiently small, there exists a meshindependent constant C such that the unique solution u h to (3) satisfies:

being the regular extension of u to Ω ′ constructed in accordance to Stein et al. [START_REF] Stein | Immersed boundary smooth extension: A high-order method for solving PDE on arbitrary smooth domains using Fourier spectral methods[END_REF].

Proof. First of all combining (3) with Proposition 3.2 we can write:

The first term in the numerator of (26) can be estimated in the following manner. Following the same steps as in Theorem 5.9 of [START_REF] Ruas | Optimal-rate finite-element solution of Dirichlet problems in curved domains with straight-edged tetrahedra[END_REF], we denote by Q h the subset of O h consisting of elements T such that T ̸ = T . Next we apply First Green's identity to a h (u ′ , v). Noticing that v is not continuous across the inter-element boundaries, and recalling the notations