Variational Formulations Yielding High-Order Finite-Element Solutions in Smooth Domains without Curved Elements ()
1. Introduction
This work deals with a new variational formulation designed for finite-element 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
or
, with Dirichlet boundary conditions, namely,
(1)
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 2 we present our method to solve the model problem with Dirichlet boundary conditions in a smooth curved two-dimensional domain with conforming Lagrange finite elements based on meshes with straight triangles, in connection with the standard Galerkin formulation. Numerical examples illustrating technique’s potential are given. 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
of
into (closed) ordinary triangles or tetrahedra, according to the value of N, satisfying the usual compatibility conditions (see e.g. [4] ). Every
is assumed to belong to a uniformly regular family of partitions. We denote by
the set
and by
the boundary of
. Letting
be the diameter of
, we set
, as usual. We also recall that if
is convex
is a proper subset of
.
2. The Two-Dimensional Case
To begin with we describe our methodology in the case where
. In order to simplify the presentation in this section we assume that
, leaving for the next one its extension to the case of an arbitrary d.
2.1. Method Description
Here we make the very reasonable assumption on the mesh that no element in
has more than one edge on
.
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
. This property is guaranteed if the points separating such portions of
are vertices of polygon
. In doing so, let
be the subset of
consisting of triangles having one edge on
. Now for every
we denote by
the set delimited by
and the edge
of T whose end-points belong to
and set
if
is not a subset of T and
otherwise (see Figure 1).
Notice that if
lies on a convex portion of
, T is a proper subset of
, while the opposite occurs if
lies on a concave portion of
. With such a definition we can assert that there is a partition
of
associated with
consisting of non overlapping sets
for
, besides the elements in
.
For convenience henceforth we refer to the
points in a triangle T which are vertices of the
equal triangles T can be subdivided into, where
for
as the lagrangian nodes of order k (cf. [4] ).
Next we introduce two function spaces
and
associated with
.
is the standard Lagrange finite element space consisting of continuous functions v defined in
that vanish on
, whose restriction to every
is a polynomial of degree less than or equal to k for
. For
Figure 1. Skin
related to a mesh triangle T next to a convex (right) or a concave (left) portion of
.
convenience we extend by zero every function
to
.
in turn is the space of functions defined in
having the pro- perties listed below.
1) The restriction of
to every
is a polynomial of degree less than or equal to k;
2) Every
is continuous in
and vanishes at the vertices of
;
3) A function
is also defined in
in such a way that its polynomial expression in
also applies to points in
;
4)
,
for every P among the
intersections with
of the line passing through the vertex
of T not belonging to
and the points M different from vertices of T subdividing the edge opposite to
into k segments of equal length (cf. Figure 2).
Remark The construction of the nodes associated with
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
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
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
be the space of polynomials defined in
of degree less than or equal to k. Provided h is small enough
, given a set of
real values
with
, there exists a unique function
that vanishes at both vertices of T located on
and at the
points P of
defined in accordance with item 4. of the above definition of
, and takes value
respectively at the
lagrangian nodes of T not located on
.
Figure 2. Construction of nodes
for space
related to lagrangian nodes
for
.
Now let us set the problem associated with spaces
and
, whose solution is an approximation of u, that is, the solution of (1). Denoting by
a sufficiently smooth extension of f to
in case this set is not empty, and renaming f by
in
, we wish to solve,
(2)
The following result is borrowed from [5] :
Proposition 2 Provided h is sufficiently small problem (2) has a unique solution.
2.2. 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
. 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
with
for an exact solution u given by
. Thus we take
and
and owing to symmetry we consider only the quarter domain given by
and
by prescribing Neumann boundary conditions on
and
. We take
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
whose edges are parallel to the coordinate axes and to the line
.
Hereunder, and in the remainder of this work we denote by
the standard mean-square norm in
of a function or a vector field
.
In Table 1 we display the absolute errors in the energy norm, namely
and in the mean-square norm, that is
for increasing
Table 1. Absolute errors in different senses for the new approach (with ordinary triangles).
values of J; more precisely
for
. 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
, as predicted therein. The error in the mean-square norm in turn tends to decrease as
, 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
with
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
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
, while in case
is a polygon it is known to be an
for sufficiently smooth solutions (see e.g. [7] ).
3. The Three-Dimensional Case
In this section we consider the solution of (1) by our method in case
.
In order to avoid non essential difficulties we make the assumption that no element in
has more than one face on
, which is nothing but reasonable.
3.1. Method Description
First of all we need some definitions regarding the set
.
Let
be the subset of
consisting of tetrahedra having one face on
and
be the subset of
of tetrahedra having exactly one edge on
. Notice that, owing to our initial assumption, no tetrahedron in
has a non empty intersection with
.
To every edge e of
we associate a plane skin
containing e, and delimited by
and e itself. Except for the fact that each skin contains an edge of
, its plane can be arbitrarily chosen. In Figure 3 we illustrate one out of three such skins corresponding to the edges of a face
or
contained in
Figure 3. Sets
,
,
for tetrahedra
with a common edge e and a tetrahedron
.
Table 2. Absolute errors in different senses for the classical approach with ordinary triangles.
, of two tetrahedra T and
belonging to
. More precisely in Figure 3 we show the skin
, e being the edge common to
and
. Further, for every
, we define a set
delimited by
, the face
and the three plane skins associated with the edges of
, as illustrated in Figure 3. In this manner we can assert that, if
is convex,
is a proper subset of
and moreover
is the union of the disjoint sets
and
(cf. Figure 3). Otherwise
is a non empty set that equals the union of certain parts of the sets
corresponding to non convex portions of
.
Next we introduce two sets of functions
and
, both associated with
.
is the standard Lagrange finite element space consisting of continuous functions v defined in
that vanish on
, whose restriction to every
is a polynomial of degree less than or equal to k for
. For convenience we extend by zero every function
to
. We recall that a function in
is uniquely defined by its values at the points which are vertices of the partition of each mesh tetrahedron into
equal tetrahedra (see e.g. [4] ). Akin to the two-dimensional case these points will be referred to as the lagrangian nodes of order k of the mesh.
in turn is a linear manifold consisting of functions defined in
satisfying the following conditions:
1) The restriction of
to every
is a polynomial of degree less than or equal to k;
2) Every
is single-valued at all the inner lagrangian nodes of the mesh, that is all its lagrangian nodes of order k except those located on
;
3) A function
is also defined in
in such a way that its polynomial expression in
also applies to points in
;
4) A function
takes the value
at any vertex S of
;
5)
and for
,
for every point P among the
intersections with
of the line passing through the vertex
of T not belonging to
and the
points M not belonging to any edge of
among the
points of
that subdivide this face (opposite to
) into
equal triangles (see illustration in Figure 4 for
);
6)
,
for every Q among the
intersections with
of the line orthogonal to e in the skin
, passing through the points
different from vertices of T, subdividing e into k equal segments, where e represents a generic edge of T contained in
(see illustration in Figure 5 for
).
Remark The construction of the nodes associated with
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
commonly used in the isoparametric technique. The main advantage of this proposal is the determination by linearity of the coordinates of the boundary nodes
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
is a non empty set is a trivial consequence of the two following results proved in [8] , where
represents the space of polynomials defined in
of degree not greater than k.
Figure 4. Construction of node
of
related to the Lagrange node M in the interior of
.
Figure 5. Construction of nodes
of
related to the Lagrange nodes
.
Proposition 3 (cf. [8] ).
Provided h is small enough
given a set of
real values
,
with
for
and
for
, there exists a unique function
that takes the value of d at the vertices S of T located on
, at the points P of
defined in accordance with item 5. for
only, and at the points Q of
defined in accordance with item 6. of the above definition of
, and takes the value
respectively at the
lagrangian nodes of T not located on
.
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
does not necessarily consist of continuous functions. This is because of the interfaces between elements in
and
. 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
instead of the corresponding lagrangian node
, in accordance with item 6. in the definition of
. On the other hand
is necessarily continuous over all other faces common to two mesh tetrahedra.
Next we set the problem associated with the space
and the manifold
, whose solution is an approximation of u, that is, the solution of (1). Extending
by a smooth
in
if necessary, and renaming
by
in any case, we wish to solve,
(3)
3.2. 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
and major radius
. This means that the torus’ inner radius
equals
and its outer radius
equals
. Hence
is given by the
equation
. We consider a test-problem with symmetry
about the z-axis, and with respect to the plane
. For this reason we may work with a computational domain given by
. A family of meshes of this domain depending on a single even integer parameter
containing
tetrahedra is generated by the following procedure. First we generate a partition of the cube
into
equal rectangular boxes 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
. This mesh with
tetrahedra is transformed into the mesh of the quarter cylinder
, following the transformation of the mesh consisting of
equal right triangles formed by the faces of the mesh elements contained in the unit cube’s section given by
, for
. The latter transformation is based on the mapping of the cartesian coordinates
into the polar coordinates
with
, using a procedure described in [8] (cf. Figure 6). Then the resulting mesh of the quarter cylinder is transformed into the mesh with
the trahedra of the half cylinder
by symmetry with respect to the plane
. Finally this mesh is transformed into the computational mesh (of an eighth of half-torus) by first mapping the cartesian coordinates
into polar coordinates
, with
Figure 6. Trace of the intermediate mesh of 1/4 cylinder on sections
,
, for
.
and
, and then the latter coordinates into new cartesian coordinates
using the relations
and
. Notice that the faces of the final tetrahedral mesh on the sections of the torus given by
, for
, form a triangular mesh of a disk with radius equal to
, having the pattern illustrated in Figure 6 for a quarter disk, taking
,
and
.
Recalling that here
, we take
,
and
. For
the exact solution is given by
. In order to enable the calculation of mean-square norms of the error in
, we extend u to
in a neighborhood of
lying outside
, taking the same expression as above. In Table 3 we display the absolute errors in the energy norm, that is
and in the mean-square norm
, for increasing values of I, namely
for
. Now we take as a reference
.
As one can observe from Table 3, here again the quality of the approxi- mations 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
, as predicted. On the other hand here again the mean-square norm of the error function
tends to decrease as
. Likewise the two-dimensional case, Table 4 in turn shows a qualitative erosion of the solution errors if
is replaced by
in (3).
4. Final Comments
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 3. Absolute errors for the new approach (with ordinary tetrahedra) in two different norms.
Table 4. Absolute errors for the classical approach with ordinary tetrahedra in two different norms.
a paper to appear shortly.
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.
3) As for the Poisson equation with homogeneous Neumann boundary conditions
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
is modified, in such a way that boundary integrals for elements
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
or in
, the inverse of an
matrix has to be computed, where
is the dimension of
. However this extra effort is not really a problem nowadays, in view of the significant progress already accomplished in Computational Linear Algebra.
Acknowledgements
The author is thankful for the support provided by CNPq-Brazil, through grant 307996/2008-5, to carry out this work, which was first presented in August 2017 at both the MFET Conference in Bad Honnef, Germany, and the French Congress of Mechanics in Lille.
Sincere thanks are due to the managing editor Hellen XU for her precious assistance.