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  for two-dimensional problems. Although the latter technique is very intuitive and has been known since the seventies (cf.   ), 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,
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.  ), 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.   ), 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.  ). 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.  ).
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.  ).
The fact that is a non empty finite-dimensional space was established in  . Furthermore the following result was also proved in the same reference:
Proposition 1 (cf.  ).
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,
The following result is borrowed from  :
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  , we implemented it taking . Then we solved Equation (1) for several test-cases already reported in  and  . 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  . 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  . 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.  ). 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.  ).
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.  ). 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  , 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.  ).
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  , namely,
Proposition 4 (cf.  )
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.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  (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  , 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.
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.