Variational Formulations Yielding High-Order Finite-Element Solutions in Smooth Domains without Curved Elements

One of the reasons for the great success of the finite element method is its versatility to deal with different types of geometries. This is particularly true of problems posed in curved domains. Nevertheless it is well-known that, for standard variational formulations, the optimal approximation properties known to hold for polytopic domains are lost, if meshes consisting of ordinary elements are still used in the case of curved domains. That is why method’s isoparametric version for meshes consisting of curved triangles or tetrahedra has been widely employed, especially in case Dirichlet boundary conditions are prescribed all over a curved boundary. However, besides geometric inconveniences, the isoparametric technique helplessly requires the manipulation of rational functions and the use of numerical integration. In this work we consider a simple alternative that bypasses these drawbacks, without eroding qualitative approximation properties. More specifically we work with a variational formulation leading to high order finite element methods based only on polynomial algebra, since they do not require the use of curved elements. Application of the new approach to Lagrange methods of arbitrary order illustrates its potential to take the best advantage of finite-element discretizations in the solution of wide classes of problems posed in curved domains.


Introduction solution methods of boundary value problems with Dirichlet boundary conditions,
posed in a two-or three-dimensional domain having a smooth curved boundary of arbitrary shape.The principle it is based upon is related to the technique called interpolated boundary conditions studied in [1] for two-dimensional problems.
Although the latter technique is very intuitive and has been known since the seventies (cf.[2] [3]), it has been of limited use so far.Among the reasons for this we could quote its difficult implementation, the lack of an extension to three-dimensional problems, and most of all, restrictions on the choice of boundary nodal points to reach optimal convergence rates.In contrast our method is simple to implement in both in two-and three-dimensional geometries.
Moreover optimality is attained very naturally in both cases for various choices of boundary nodal points.
In order to allow an easier description of our methodology, thereby avoiding non essential technicalities, we consider as a model the Poisson equation in an N-dimensional smooth domain Ω with boundary Γ , for 2 N = or 3 N = , with Dirichlet boundary conditions, namely, in on , where f and d are given functions defined in Ω and on Γ , having suitable regularity properties.
Here ( 1) is supposed to be solved by different N-simplex based finite element methods, incorporating degrees of freedom other than function values at the mesh vertices.For instance, if standard quadratic Lagrange finite elements are employed, it is well-known that approximations of an order not greater than 1.5 in the energy norm are generated (cf.[4]), in contrast to the second order ones that apply to the case of a polygonal or polyhedral domain, assuming that the solution is sufficiently smooth.If we are to recover the optimal second order approximation property something different has to be done.Since long the isoparametric version of the finite element method for meshes consisting of curved triangles or tetrahedra (cf.[3] [4]), has been considered as the ideal way to achieve this.It turns out that, besides a more elaborated description of the mesh, the isoparametric technique inevitably leads to the integration of rational functions to compute the system matrix, which raises the delicate question on how to choose the right numerical quadrature formula in the master element.In contrast, in the technique to be introduced in this paper exact numerical integration can always be used for this purpose, since we only have to deal with polynomial integrands.Moreover the element geometry remains the same as in the case of polygonal or polyhedral domains.It is noteworthy that both advantages are conjugated with the fact that no erosion of qualitative approximation properties results from the application of our technique, as compared to the equivalent isoparametric one.
An outline of the paper is as follows.In Section In Section 3 we extend the approach adopted in Section 2 to the three-dimensional case including also numerical experimentation.We conclude in Section 4 with some comments on possible extensions of the methodology studied in this work.
In the remainder of this paper we will be given partitions h  of Ω into (closed) ordinary triangles or tetrahedra, according to the value of N, satisfying the usual compatibility conditions (see e.g.[4]).Every h  is assumed to belong to a uniformly regular family of partitions.We denote by and by h Γ the boundary of h Ω .Letting T h be the diameter of h T ∈  , we set

The Two-Dimensional Case
To begin with we describe our methodology in the case where 2 N = .In order to simplify the presentation in this section we assume that 0 d ≡ , leaving for the next one its extension to the case of an arbitrary d.

Method Description
Here we make the very reasonable assumption on the mesh that no element in h  has more than one edge on h Γ .
We also need some definitions regarding the skin ( ) ( ) . First of all, in order to avoid non essential difficulties, we assume that the mesh is constructed in such a way that convex and concave portions of Γ correspond to convex and concave portions of h Γ .This property is guaranteed if the points separating such portions of Γ are vertices of polygon h Ω .In doing so, let h  be the subset of h  consisting of triangles having one edge on h Γ .Now for every h T ∈  we denote by T ∆ the set delimited by Γ and the edge T e of T whose end-points belong to Γ and set : k > as the lagrangian nodes of order k (cf.[4]).
Next we introduce two function spaces h V and h W associated with h  .h V is the standard Lagrange finite element space consisting of continuous functions v defined in h Ω that vanish on h Γ , whose restriction to every h T ∈  is a polynomial of degree less than or equal to k for for every P among the of the line passing through the vertex T O of T not belonging to Γ and the points M different from vertices of T subdividing the edge opposite to T O into k segments of equal length (cf. Figure 2).
Remark The construction of the nodes associated with h W located on Γ advocated in item 4 is not mandatory.Notice that it differs from the intuitive construction of such nodes lying on normals to edges of h Γ commonly used in the isoparametric technique.The main advantage of this proposal is an easy determination of boundary node coordinates by linearity, using a supposedly available analytical expression of Γ .Nonetheless the choice of boundary nodes ensuring our method's optimality is really wide, in contrast to the restrictions inherent to the interpolated boundary condition method (cf.[1]).
The fact that h W is a non empty finite-dimensional space was established in [5].Furthermore the following result was also proved in the same reference: Proposition 1 (cf.[5]).

Let ( )
k T  be the space of polynomials defined in h T ∈  of degree less than or equal to k. Provided h is small enough , there exists a unique function

Now let us set the problem associated with spaces h
V and h W , whose solution is an approximation of u, that is, the solution of (1).Denoting by f ′ a sufficiently smooth extension of f to \ h Ω Ω in case this set is not empty, and renaming f by f ′ in Ω , we wish to solve,

Method Assessment
In order to illustrate the accuracy and the optimal order of the method described in the previous subsection rigorously demonstrated in [5], we implemented it taking 2 k = .Then we solved Equation (1) for several test-cases already reported in [5] and [6].Here we present the results for the following one: Ω is the ellipse delimited by the curve ( ) and compute with quasi-uniform meshes defined by a single integer parameter J, constructed by the procedure described in [5].Roughly speaking the mesh of the quarter domain is the polar coordinate counterpart of the standard uniform mesh of the unit square ( ) ( ) 0,1 0,1 × whose edges are parallel to the coordinate axes and to the line x y = .
Hereunder, and in the remainder of this work we denote by h  the standard mean-square norm in h Ω of a function or a vector field  .
In Table 1  .We also show the evolution of the maximum absolute errors at the mesh nodes.
As one infers from Table 1, the approximations obtained with our method perfectly conform to the theoretical estimate given in [5].Indeed as J increases the errors in the energy norm decrease roughly as ( ) 1 J , as predicted therein.
The error in the mean-square norm in turn tends to decrease as ( ) 1 J , while the maximum absolute error at the nodes seem to behave like an ( ) Now in order to rule out any particularity inherent to the above test-problem, we also solved it using the classical approach, that is, by replacing h W with h V in (2).In Table 2 we display the same kind of results as in Table 1 for this approach.
Table 2 confirms the error decrease in the energy norm like an ( ) O h as predicted in classical texts (cf.[4]).The behavior of the classical approach deteriorates even more, as compared to the new approach, when the errors are measured in the mean-square norm, whose order seem to decrease from three to two.The quality of the maximum nodal absolute errors in turn are not affected at all by the way boundary conditions are handled.Actually in both cases this error is roughly an ( ) O h , while in case Ω is a polygon it is known to be an ( )

The Three-Dimensional Case
In this section we consider the solution of (1) by our method in case 3 N = .In order to avoid non essential difficulties we make the assumption that no element in h  has more than one face on h Γ , which is nothing but reasonable.

Method Description
First of all we need some definitions regarding the set ( ) ( )  be the subset of h  consisting of tetrahedra having one face on h Γ and h  be the subset of \ h h   of tetrahedra having exactly one edge on h Γ .
Notice that, owing to our initial assumption, no tetrahedron in has a non empty intersection with h Γ .
To every edge e of h Γ we associate a plane skin e δ containing e, and delimited by Γ and e itself.Except for the fact that each skin contains an edge of h Γ , its plane can be arbitrarily chosen.In Figure 3 we illustrate one out of three such skins corresponding to the edges of a face T F or T F ′ contained in Table 2. Absolute errors in different senses for the classical approach with ordinary triangles.δ , e being the edge common to T F and T F ′ .Further, for every h T ∈  , we define a set T ∆ delimited by Γ , the face T F and the three plane skins associated with the edges of T F , as illustrated in Figure 3.In this manner we can assert that, if Ω is convex, h Ω is a proper subset of Ω and moreover Ω is the union of the disjoint sets h Ω and h T T∈ ∆   (cf. Figure 3).Otherwise \ h Ω Ω is a non empty set that equals the union of certain parts of the sets T ∆ corresponding to non convex portions of Γ .V is the standard Lagrange finite element space consisting of continuous functions v defined in h Ω that vanish on h Γ , whose restriction to every h T ∈  is a polynomial of degree less than or equal to k for 2 k ≥ .For convenience we extend by zero every function

Next we introduce two sets of functions h V and
V is uniquely defined by its values at the points which are vertices of the partition of each mesh tetrahedron into 3  k equal tetrahedra (see e.g.[4]).for with Γ of the line orthogonal to e in the skin e δ , passing through the points

M e
∈ different from vertices of T, subdividing e into k equal segments, where e represents a generic edge of T contained in h Γ (see illustration in Figure 5 for Remark The construction of the nodes associated with d h W located on Γ advocated in items 5. and 6. is not mandatory.Notice that it differs from the intuitive construction of such nodes lying on normals to faces of h Γ commonly used in the isoparametric technique.The main advantage of this proposal is the determination by linearity of the coordinates of the boundary nodes P in the case of item 5. Nonetheless, akin to the two-dimensional case, the choice of boundary nodes ensuring our method's optimality is absolutely very wide.

The fact that d h
W is a non empty set is a trivial consequence of the two following results proved in [8], where   Proposition 3 (cf.[8]).Provided h is small enough that takes the value of d at the vertices of T located on Γ at the points P of Γ defined in accordance with item 5. for h T ∈  only, and at the points Q of Γ defined in accordance with item 6. of the above definition of A well-posedness result analogous to Proposition 2.2 holds for problem (3), according to [8], namely, Proposition 4 (cf.[8]) As long as h is sufficiently small problem (3) has a unique solution.Remark It is important to stress that, in contrast to its two-dimensional counterpart, the set d h W does not necessarily consist of continuous functions.This is because of the interfaces between elements in h  and h  .Indeed a function is not forcibly single-valued at all the lagrangian nodes located on one such an interface, owing to the enforcement of the boundary condition at the points Q ∈ Γ instead of the corresponding lagrangian node h M ∈ Γ , in accordance with item 6. in the definition of d h W .On the other hand w is necessarily continuous over all other faces common to two mesh tetrahedra.
Next we set the problem associated with the space h V and the manifold whose solution is an approximation of u, that is, the solution of (1).Extending f by a smooth f ′ in \ h Ω Ω if necessary, and renaming f by f ′ in any case, we wish to solve,

Method Assessment
In this sub-section we assess the behavior of the new method, by solving the Poisson equation in a non convex domain.Now we consider a non polynomial exact solution.More precisely (1) is solved in the torus Ω with minor radius m r and major radius M r .This means that the torus' inner radius i r equals M m r r − and its outer radius e r equals M m r r + .Hence Γ is given by the equation ( ) r y z = + , using a procedure described in [8] (cf. Figure 6).Then the resulting mesh of the quarter cylinder is transformed into the mesh with  Recalling that here .In order to enable the calculation of mean-square norms of the error in Ω , we extend u to u′ in a neighborhood of Γ lying outside Ω , taking the same expression as above.In Table 3  As one can observe from Table 3, here again the quality of the approximations obtained with the new method is in very good agreement with the theoretical result in [8], for as I increases the errors in the energy norm decrease roughly as  4 in turn shows a qualitative erosion of the solution errors if  2) The technique studied in this paper is also particularly handy, to treat problems posed in curved domains in terms of vector fields, such as the linear elasticity system and Maxwell's equations of electromagnetism.The same remark applies to multi-field systems such as the Navier-Stokes equations, and more generally mixed formulations of several types with curved boundaries, to be approximated by the finite element method.

Final Comments
3) As for the Poisson equation with homogeneous Neumann boundary conditions 0 u n ∂ ∂ = on Γ (provided f satisfies the underlying scalar condition) our method coincides with the standard Lagrange finite element method.Notice that if inhomogeneous Neumann boundary conditions are prescribed, optimality can only be recovered if the linear form h F is modified, in such a way that boundary integrals for elements h T ∈  are shifted to the curved boundary portion sufficiently close to Γ of an extension or reduction of T. But this is an issue that has nothing to do with our method, which is basically aimed at resolving those related to the prescription of degrees of freedom in the case of Dirichlet boundary conditions.
4) Finally we note that our method leads to linear systems of equations with a non symmetric matrix, even when the original problem is symmetric.Moreover in order to compute the element matrix and right side vector for an element T in h  or in h  , the inverse of an k k n n × matrix has to be computed, where k n is the dimension of ( ) k P T .However this extra effort is not really a problem nowadays, in view of the significant progress already accomplished in Computational Linear Algebra.

Figure 1 .
Figure 1.Skin T ∆ related to a mesh triangle T next to a convex (right) or a concave (left) h both vertices of T located on Γ and at the 1 k − points P of Γ defined in accordance with item 4. of the above definition of h W , and takes value i b respectively at the k m lagrangian nodes of T not located on h Γ .

Figure 2 .
Figure 2. Construction of nodes P ∈ Γ for space h W related to lagrangian nodes h M ∈ Γ for 3 k = .
and 0 d ≡ and owing to symmetry we consider only the quarter domain given by 0

4 hΓ
, of two tetrahedra T and T ′ belonging to h  .More precisely in Figure3we show the skin e h   of degree not greater than k.

Figure 4 .
Figure 4. Construction of node P ∈ Γ of d h W related to the

Figure 5 .
Figure 5. Construction of nodes

W
, and takes the value i b respectively at the k m lagrangian nodes of T not located on h Γ .

3 6I 2 I.
We consider a test-problem with symmetry about the z-axis, and with respect to the plane 0 z = .For this reason we may work with a computational domain given by of meshes of this domain depending on a single even integer parameter I containing tetrahedra is generated by the following procedure.First we generate a partition of the cube ( ) ( ) ( ) by subdividing the edges parallel to the x-axis, the y-axis and the z-axis into 2I, I/2 and I/2 equal segments, respectively.Then each box is subdivided into six tetrahedra having an edge parallel to the line 4x y z = = .This mesh with3 3Itetrahedra is transformed into the mesh of the quarter cylinder following the transformation of the mesh consisting of 2 equal right triangles formed by the faces of the mesh elements contained in the unit cube's section given by The latter transformation is based on the mapping of the cartesian coordinates ( ) , y z into the polar coordinates ( )
the faces of the final tetrahedral mesh on the sections of the torus given by triangular mesh of a disk with radius equal to m r , having the pattern illustrated in Figure6for a quarter disk, taking 4

2 1 I 1 I
, as predicted.On the other hand here again the mean-square norm of the error function .Likewise the two-dimensional case, Table E−3 0.143425 E−3 0.343823 E−4 0.834136 E−5 a paper to appear shortly.
in turn is the space of functions defined in h

Table 1 .
Absolute errors in different senses for the new approach (with ordinary triangles).
we display the absolute errors in the energy norm, namely ( )

Table 3 .
Absolute errors for the new approach (with ordinary tetrahedra) in two different norms.
1) The method addressed in this work to solve the Poisson equation with Dirichlet boundary conditions in curved domains with classical Lagrange finite elements provides a simple and reliable manner to overcome technical difficulties brought about by more complicated problems and interpolations.For example, Hermite finite element methods to solve fourth order problems in curved domains with normal derivative degrees of freedom can also be dealt with very easily by means of our new method.The author intends to show this in

Table 4 .
Absolute errors for the classical approach with ordinary tetrahedra in two different norms.