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

Vitoriano Ruas^{1,2,3}

Show more

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 $\Omega $ with boundary $\Gamma $ , for $N=2$ or $N=3$ , with Dirichlet boundary conditions, namely,

$\{\begin{array}{l}-\Delta u=f\text{\hspace{0.05em}}\text{\hspace{0.17em}}\text{in}\text{\hspace{0.17em}}\text{\hspace{0.05em}}\Omega \hfill \\ u=d\text{\hspace{0.05em}}\text{\hspace{0.17em}}\text{on}\text{\hspace{0.17em}}\text{\hspace{0.05em}}\Gamma ,\hfill \end{array}$ (1)

where f and d are given functions defined in $\Omega $ and on $\Gamma $ , 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 ${\mathcal{T}}_{h}$ of $\Omega $ into (closed) ordinary triangles or tetrahedra, according to the value of N, satisfying the usual compatibility conditions (see e.g. [4] ). Every ${\mathcal{T}}_{h}$ is assumed to belong to a uniformly regular family of partitions. We denote by ${\Omega}_{h}$ the set ${\cup}_{T\in {\mathcal{T}}_{h}}T$ and by ${\Gamma}_{h}$ the boundary of ${\Omega}_{h}$ . Letting ${h}_{T}$ be the diameter of $T\in {\mathcal{T}}_{h}$ , we set $h:={\mathrm{max}}_{T\in {\mathcal{T}}_{h}}{h}_{T}$ , as usual. We also recall that if $\Omega $ is convex ${\Omega}_{h}$ is a proper subset of $\Omega $ .

2. The Two-Dimensional Case

To begin with we describe our methodology in the case where $N=2$ . In order to simplify the presentation in this section we assume that $d\equiv 0$ , 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 ${\mathcal{T}}_{h}$ has more than one edge on ${\Gamma}_{h}$ .

We also need some definitions regarding the skin $\left(\Omega \backslash {\Omega}_{h}\right)\cup \left({\Omega}_{h}\backslash \Omega \right)$ . 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 $\Gamma $ correspond to convex and concave portions of ${\Gamma}_{h}$ . This property is guaranteed if the points separating such portions of $\Gamma $ are vertices of polygon ${\Omega}_{h}$ . In doing so, let ${\mathcal{S}}_{h}$ be the subset of ${\mathcal{T}}_{h}$ consisting of triangles having one edge on ${\Gamma}_{h}$ . Now for every $T\in {\mathcal{S}}_{h}$ we denote by ${\Delta}_{T}$ the set delimited by $\Gamma $ and the edge ${e}_{T}$ of T whose end-points belong to $\Gamma $ and set ${T}^{\prime}\mathrm{:}=T\cup {\Delta}_{T}$ if ${\Delta}_{T}$ is not a subset of T and ${T}^{\prime}:=\stackrel{\xaf}{T\backslash {\Delta}_{T}}$ otherwise (see Figure 1).

Notice that if ${e}_{T}$ lies on a convex portion of ${\Gamma}_{h}$ , T is a proper subset of ${T}^{\prime}$ , while the opposite occurs if ${e}_{T}$ lies on a concave portion of ${\Gamma}_{h}$ . With such a definition we can assert that there is a partition ${{\mathcal{T}}^{\prime}}_{h}$ of $\Omega $ associated with ${\mathcal{T}}_{h}$ consisting of non overlapping sets ${T}^{\prime}$ for $T\in {\mathcal{S}}_{h}$ , besides the elements in ${\mathcal{T}}_{h}\backslash {\mathcal{S}}_{h}$ .

For convenience henceforth we refer to the ${n}_{k}$ points in a triangle T which are vertices of the ${k}^{2}$ equal triangles T can be subdivided into, where ${n}_{k}:=\left(k+2\right)\left(k+1\right)/2$ for $k>1$ as the lagrangian nodes of order k (cf. [4] ).

Next we introduce two function spaces ${V}_{h}$ and ${W}_{h}$ associated with ${\mathcal{T}}_{h}$ .

${V}_{h}$ is the standard Lagrange finite element space consisting of continuous functions v defined in ${\Omega}_{h}$ that vanish on ${\Gamma}_{h}$ , whose restriction to every $T\in {\mathcal{T}}_{h}$ is a polynomial of degree less than or equal to k for $k\ge 2$ . For

Figure 1. Skin ${\Delta}_{T}$ related to a mesh triangle T next to a convex (right) or a concave (left) portion of $\Gamma $ .

convenience we extend by zero every function $v\in {V}_{h}$ to $\Omega \backslash {\Omega}_{h}$ .

${W}_{h}$ in turn is the space of functions defined in ${\Omega}_{h}\cup \Omega $ having the pro- perties listed below.

1) The restriction of $w\in {W}_{h}$ to every $T\in {\mathcal{T}}_{h}$ is a polynomial of degree less than or equal to k;

2) Every $w\in {W}_{h}$ is continuous in ${\Omega}_{h}$ and vanishes at the vertices of ${\Gamma}_{h}$ ;

3) A function $w\in {W}_{h}$ is also defined in $\Omega \backslash {\Omega}_{h}$ in such a way that its polynomial expression in $T\in {\mathcal{S}}_{h}$ also applies to points in ${\Delta}_{T}$ ;

4) $\forall T\in {\mathcal{S}}_{h}$ , $w\left(P\right)=0$ for every P among the $k-1$ intersections with $\Gamma $ of the line passing through the vertex ${O}_{T}$ of T not belonging to $\Gamma $ and the points M different from vertices of T subdividing the edge opposite to ${O}_{T}$ into k segments of equal length (cf. Figure 2).

Remark The construction of the nodes associated with ${W}_{h}$ located on $\Gamma $ advocated in item 4 is not mandatory. Notice that it differs from the intuitive construction of such nodes lying on normals to edges of ${\Gamma}_{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 $\Gamma $ . 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 ${W}_{h}$ 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 ${\mathcal{P}}_{k}\left(T\right)$ be the space of polynomials defined in $T\in {\mathcal{S}}_{h}$ of degree less than or equal to k. Provided h is small enough $\forall T\in {\mathcal{S}}_{h}$ , given a set of ${m}_{k}$ real values ${b}_{i},i=1,\cdots ,{m}_{k}$ with ${m}_{k}=\left(k+1\right)k/2$ , there exists a unique function ${w}_{T}\in {\mathcal{P}}_{k}\left(T\right)$ that vanishes at both vertices of T located on $\Gamma $ and at the $k-1$ points P of $\Gamma $ defined in accordance with item 4. of the above definition of ${W}_{h}$ , and takes value ${b}_{i}$ respectively at the ${m}_{k}$ lagrangian nodes of T not located on ${\Gamma}_{h}$ .

Figure 2. Construction of nodes $P\in \Gamma $ for space ${W}_{h}$ related to lagrangian nodes $M\in {\Gamma}_{h}$ for $k=3$ .

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

$\{\begin{array}{l}\text{\hspace{0.05em}}\text{Find}\text{\hspace{0.17em}}\text{\hspace{0.05em}}{u}_{h}\in {W}_{h}\text{\hspace{0.05em}}\text{\hspace{0.17em}}\text{such}\text{\hspace{0.17em}}\text{that}\text{\hspace{0.05em}}\text{\hspace{0.17em}}{a}_{h}\left({u}_{h},v\right)={F}_{h}\left(v\right)\text{\hspace{0.17em}}\forall v\in {V}_{h}\\ \text{where}\text{\hspace{0.17em}}\text{\hspace{0.05em}}{a}_{h}\left(w,v\right):={\displaystyle {\int}_{{\Omega}_{h}}}grad\text{\hspace{0.17em}}w\cdot grad\text{\hspace{0.17em}}v\text{\hspace{0.05em}}\text{\hspace{0.17em}}\text{and}\text{\hspace{0.17em}}\text{\hspace{0.05em}}{F}_{h}\left(v\right):={\displaystyle {\int}_{{\Omega}_{h}}}{f}^{\prime}\text{\hspace{0.05em}}v\end{array}$ (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 $k=2$ . Then we solved Equation (1) for several test-cases already reported in [5] and [6] . Here we present the results for the following one:

$\Omega $ is the ellipse delimited by the curve ${\left(x/e\right)}^{2}+{y}^{2}=1$ with $e>0$ for an exact solution u given by $u=\left({e}^{2}-{e}^{2}{x}^{2}-{y}^{2}\right)\left({e}^{2}-{x}^{2}-{e}^{2}{y}^{2}\right)$ . Thus we take $f:=-\Delta u$ and $d\equiv 0$ and owing to symmetry we consider only the quarter domain given by $x>0$ and $y>0$ by prescribing Neumann boundary conditions on $x=0$ and $y=0$ . We take $e=0.5$ 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 $\left(\mathrm{0,1}\right)\times \left(\mathrm{0,1}\right)$ 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 ${\Vert \mathcal{F}\Vert}_{h}$ the standard mean-square norm in ${\Omega}_{h}$ of a function or a vector field $\mathcal{F}$ .

In Table 1 we display the absolute errors in the energy norm, namely ${\Vert grad\left(u-{u}_{h}\right)\Vert}_{h}$ and in the mean-square norm, that is ${\Vert u-{u}_{h}\Vert}_{h}$ for increasing

Table 1. Absolute errors in different senses for the new approach (with ordinary triangles).

values of J; more precisely $J={2}^{m}$ for $m=3,4,5,6$ . 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 ${\left(1/J\right)}^{2}$ , as predicted therein. The error in the mean-square norm in turn tends to decrease as ${\left(1/J\right)}^{3}$ , while the maximum absolute error at the nodes seem to behave like an $O\left({h}^{2}\right)$ .

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 ${W}_{h}$ with ${V}_{h}$ 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\left({h}^{1.5}\right)$ 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\left({h}^{2}\right)$ , while in case $\Omega $ is a polygon it is known to be an $O\left({h}^{3}\right)$ 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 $N=3$ .

In order to avoid non essential difficulties we make the assumption that no element in ${\mathcal{T}}_{h}$ has more than one face on ${\Gamma}_{h}$ , which is nothing but reasonable.

3.1. Method Description

First of all we need some definitions regarding the set $\left(\Omega \backslash {\Omega}_{h}\right)\cup \left({\Omega}_{h}\backslash \Omega \right)$ .

Let ${\mathcal{S}}_{h}$ be the subset of ${\mathcal{T}}_{h}$ consisting of tetrahedra having one face on ${\Gamma}_{h}$ and ${\mathcal{R}}_{h}$ be the subset of ${\mathcal{T}}_{h}\backslash {\mathcal{S}}_{h}$ of tetrahedra having exactly one edge on ${\Gamma}_{h}$ . Notice that, owing to our initial assumption, no tetrahedron in ${\mathcal{T}}_{h}\backslash \left[{\mathcal{S}}_{h}\cup {\mathcal{R}}_{h}\right]$ has a non empty intersection with ${\Gamma}_{h}$ .

To every edge e of ${\Gamma}_{h}$ we associate a plane skin ${\delta}_{e}$ containing e, and delimited by $\Gamma $ and e itself. Except for the fact that each skin contains an edge of ${\Gamma}_{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 ${F}_{T}$ or ${F}_{{T}^{\prime}}$ contained in

Figure 3. Sets ${\Delta}_{T}$ , ${\Delta}_{{T}^{\prime}}$ , ${\delta}_{e}$ for tetrahedra $T\mathrm{,}{T}^{\prime}\in {\mathcal{S}}_{h}$ with a common edge e and a tetrahedron ${T}^{\u2033}\in {\mathcal{R}}_{h}$ .

Table 2. Absolute errors in different senses for the classical approach with ordinary triangles.

${\Gamma}_{h}$ , of two tetrahedra T and ${T}^{\prime}$ belonging to ${\mathcal{S}}_{h}$ . More precisely in Figure 3 we show the skin ${\delta}_{e}$ , e being the edge common to ${F}_{T}$ and ${F}_{{T}^{\prime}}$ . Further, for every $T\in {\mathcal{S}}_{h}$ , we define a set ${\Delta}_{T}$ delimited by $\Gamma $ , the face ${F}_{T}$ and the three plane skins associated with the edges of ${F}_{T}$ , as illustrated in Figure 3. In this manner we can assert that, if $\Omega $ is convex, ${\Omega}_{h}$ is a proper subset of $\Omega $ and moreover $\Omega $ is the union of the disjoint sets ${\Omega}_{h}$ and $\cup}_{T\in {\mathcal{T}}_{h}}{\Delta}_{T$ (cf. Figure 3). Otherwise ${\Omega}_{h}\backslash \Omega $ is a non empty set that equals the union of certain parts of the sets ${\Delta}_{T}$ corresponding to non convex portions of $\Gamma $ .

Next we introduce two sets of functions ${V}_{h}$ and ${W}_{h}^{d}$ , both associated with ${\mathcal{T}}_{h}$ .

${V}_{h}$ is the standard Lagrange finite element space consisting of continuous functions v defined in ${\Omega}_{h}$ that vanish on ${\Gamma}_{h}$ , whose restriction to every $T\in {\mathcal{T}}_{h}$ is a polynomial of degree less than or equal to k for $k\ge 2$ . For convenience we extend by zero every function $v\in {V}_{h}$ to $\Omega \backslash {\Omega}_{h}$ . We recall that a function in ${V}_{h}$ is uniquely defined by its values at the points which are vertices of the partition of each mesh tetrahedron into ${k}^{3}$ 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.

${W}_{h}^{d}$ in turn is a linear manifold consisting of functions defined in ${\Omega}_{h}\cup \Omega $ satisfying the following conditions:

1) The restriction of $w\in {W}_{h}^{d}$ to every $T\in {\mathcal{T}}_{h}$ is a polynomial of degree less than or equal to k;

2) Every $w\in {W}_{h}^{d}$ 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 ${\Gamma}_{h}$ ;

3) A function $w\in {W}_{h}^{d}$ is also defined in $\Omega \backslash {\Omega}_{h}$ in such a way that its polynomial expression in $T\in {\mathcal{S}}_{h}$ also applies to points in ${\Delta}_{T}$ ;

4) A function $w\in {W}_{h}^{d}$ takes the value $d\left(S\right)$ at any vertex S of ${\Gamma}_{h}$ ;

5) $\forall T\in {\mathcal{S}}_{h}$ and for $k>2$ , $w\left(P\right)=d\left(P\right)$ for every point P among the $\left(k-1\right)\left(k-2\right)/2$ intersections with $\Gamma $ of the line passing through the vertex ${O}_{T}$ of T not belonging to $\Gamma $ and the $\left(k-1\right)\left(k-2\right)/2$ points M not belonging to any edge of ${F}_{T}$ among the $\left(k+2\right)\left(k+1\right)/2$ points of ${F}_{T}$ that subdivide this face (opposite to ${O}_{T}$ ) into ${k}^{2}$ equal triangles (see illustration in Figure 4 for $k=3$ );

6) $\forall T\in {\mathcal{S}}_{h}\cup {\mathcal{R}}_{h}$ , $w\left(Q\right)=d\left(Q\right)$ for every Q among the $k-1$ intersections with $\Gamma $ of the line orthogonal to e in the skin ${\delta}_{e}$ , passing through the points $M\in e$ different from vertices of T, subdividing e into k equal segments, where e represents a generic edge of T contained in ${\Gamma}_{h}$ (see illustration in Figure 5 for $k=3$ ).

Remark The construction of the nodes associated with ${W}_{h}^{d}$ located on $\Gamma $ 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 ${\Gamma}_{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 ${W}_{h}^{d}$ is a non empty set is a trivial consequence of the two following results proved in [8] , where ${\mathcal{P}}_{k}\left(T\right)$ represents the space of polynomials defined in $T\in {\mathcal{S}}_{h}\cup {\mathcal{R}}_{h}$ of degree not greater than k.

Figure 4. Construction of node $P\in \Gamma $ of ${W}_{h}^{d}$ related to the Lagrange node M in the interior of ${F}_{T}\subset {\Gamma}_{h}$ .

Figure 5. Construction of nodes $Q\in \Gamma \cap \stackrel{\xaf}{{\delta}_{e}}$ of ${W}_{h}^{d}$ related to the Lagrange nodes $M\in e\subset {\Gamma}_{h}$ .

Proposition 3 (cf. [8] ).

Provided h is small enough $\forall T\in {\mathcal{S}}_{h}\cup {\mathcal{R}}_{h}$ given a set of ${m}_{k}$ real values ${b}_{i}$ , $i=1,\cdots ,{m}_{k}$ with ${m}_{k}=k\left(k+1\right)\left(k+2\right)/6$ for $T\in {\mathcal{S}}_{h}$ and ${m}_{k}=\left(k+1\right)\left(k+2\right)\left(k+3\right)/6-\left(k+1\right)$ for $T\in {\mathcal{R}}_{h}$ , there exists a unique function ${w}_{T}\in {P}_{k}\mathrm{(}T\mathrm{)}$ that takes the value of d at the vertices S of T located on $\Gamma $ , at the points P of $\Gamma $ defined in accordance with item 5. for $T\in {\mathcal{S}}_{h}$ only, and at the points Q of $\Gamma $ defined in accordance with item 6. of the above definition of ${W}_{h}^{d}$ , and takes the value ${b}_{i}$ respectively at the ${m}_{k}$ lagrangian nodes of T not located on ${\Gamma}_{h}$ .

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 ${W}_{h}^{d}$ does not necessarily consist of continuous functions. This is because of the interfaces between elements in ${\mathcal{R}}_{h}$ and ${\mathcal{S}}_{h}$ . Indeed a function $w\in {W}_{h}^{d}$ 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\in \Gamma $ instead of the corresponding lagrangian node $M\in {\Gamma}_{h}$ , in accordance with item 6. in the definition of ${W}_{h}^{d}$ . 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 ${V}_{h}$ and the manifold ${W}_{h}^{d}$ , whose solution is an approximation of u, that is, the solution of (1). Extending $f$ by a smooth ${f}^{\prime}$ in ${\Omega}_{h}\backslash \Omega $ if necessary, and renaming $f$ by ${f}^{\prime}$ in any case, we wish to solve,

$\{\begin{array}{l}\text{\hspace{0.05em}}\text{Find}\text{\hspace{0.17em}}\text{\hspace{0.05em}}{u}_{h}\in {W}_{h}^{d}\text{\hspace{0.05em}}\text{\hspace{0.17em}}\text{such}\text{\hspace{0.17em}}\text{that}\text{\hspace{0.05em}}\text{\hspace{0.17em}}{a}_{h}\left({u}_{h},v\right)={F}_{h}\left(v\right)\text{\hspace{0.17em}}\forall v\in {V}_{h}\\ \text{\hspace{0.05em}}\text{where}\text{\hspace{0.17em}}\text{\hspace{0.05em}}{a}_{h}\left(w,v\right):={\displaystyle {\int}_{{\Omega}_{h}}}grad\text{\hspace{0.17em}}w\cdot grad\text{\hspace{0.17em}}v\text{\hspace{0.05em}}\text{\hspace{0.17em}}\text{and}\text{\hspace{0.17em}}\text{\hspace{0.05em}}{F}_{h}\left(v\right):={\displaystyle {\int}_{{\Omega}_{h}}}{f}^{\prime}\text{\hspace{0.05em}}v\end{array}$ (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 $\Omega $ with minor radius ${r}_{m}$ and major radius ${r}_{M}$ . This means that the torus’ inner radius ${r}_{i}$ equals ${r}_{M}-{r}_{m}$ and its outer radius ${r}_{e}$ equals ${r}_{M}+{r}_{m}$ . Hence $\Gamma $ is given by the

equation ${\left({r}_{M}-\sqrt{{x}^{2}+{y}^{2}}\right)}^{2}+{z}^{2}={r}_{m}^{2}$ . We consider a test-problem with symmetry

about the z-axis, and with respect to the plane $z=0$ . For this reason we may work with a computational domain given by $\left\{\left(x,y,z\right)\in \Omega |z\ge 0;\text{\hspace{0.17em}}0\le \theta \le \text{\pi}/4\text{\hspace{0.17em}}\text{with}\text{\hspace{0.17em}}\text{\hspace{0.05em}}\theta =a\mathrm{tan}\left(y/x\right)\right\}$ . A family of meshes of this domain depending on a single even integer parameter $I$ containing $6{I}^{3}$ tetrahedra is generated by the following procedure. First we generate a partition of the cube $\left(\mathrm{0,1}\right)\times \left(\mathrm{0,1}\right)\times \left(\mathrm{0,1}\right)$ into ${I}^{3}/2$ 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 $4x=y=z$ . This mesh with $3{I}^{3}$ tetrahedra is transformed into the mesh of the quarter cylinder $\left\{\left(x,y,z\right)|0\le x\le 1,\text{\hspace{0.17em}}y\ge 0,\text{\hspace{0.17em}}z\ge 0,\text{\hspace{0.17em}}{y}^{2}+{z}^{2}\le 1\right\}$ , following the transformation of the mesh consisting of ${I}^{2}/2$ equal right triangles formed by the faces of the mesh elements contained in the unit cube’s section given by $x=j/\left(2I\right)$ , for $j=0,1,\cdots ,2I$ . The latter transformation is based on the mapping of the cartesian coordinates $\left(y\mathrm{,}z\right)$ into the polar coordinates $\left(r\mathrm{,}\phi \right)$ with

$r=\sqrt{{y}^{2}+{z}^{2}}$ , using a procedure described in [8] (cf. Figure 6). Then the resulting mesh of the quarter cylinder is transformed into the mesh with $6{I}^{3}$ the trahedra of the half cylinder $\left\{\left(x,y,z\right)|0\le x\le 1,\text{\hspace{0.17em}}-1\le y\le 1,\text{\hspace{0.17em}}z\ge 0,\text{\hspace{0.17em}}{y}^{2}+{z}^{2}\le 1\right\}$ by symmetry with respect to the plane $y=0$ . Finally this mesh is transformed into the computational mesh (of an eighth of half-torus) by first mapping the cartesian coordinates $\left(x\mathrm{,}y\right)$ into polar coordinates $\left(\rho ,\theta \right)$ , with $\rho ={r}_{M}+y{r}_{m}$

Figure 6. Trace of the intermediate mesh of 1/4 cylinder on sections $x=j/\left(2I\right)$ , $0\le j\le 2I$ , for $I=4$ .

and $\theta =x\text{\pi}/4$ , and then the latter coordinates into new cartesian coordinates $\left(x\mathrm{,}y\right)$ using the relations $x=\rho \mathrm{cos}\theta $ and $y=\rho \mathrm{sin}\theta $ . Notice that the faces of the final tetrahedral mesh on the sections of the torus given by $\theta =j\text{\hspace{0.05em}}\text{\pi}/\left(8I\right)$ , for $j=0,1,\cdots ,2I$ , form a triangular mesh of a disk with radius equal to ${r}_{m}$ , having the pattern illustrated in Figure 6 for a quarter disk, taking $I=4$ , $\theta =0$ and ${r}_{m}=1$ .

Recalling that here $\rho =\sqrt{{x}^{2}+{y}^{2}}$ , we take ${r}_{M}=5/6$ , ${r}_{m}=1/6$ and ${f}^{\prime}=6-5/\left(3\rho \right)$ . For $d\equiv 0$ the exact solution is given by $u=1/36-{z}^{2}-{\left(5/6-\rho \right)}^{2}$ . In order to enable the calculation of mean-square norms of the error in $\Omega $ , we extend u to ${u}^{\prime}$ in a neighborhood of $\Gamma $ lying outside $\Omega $ , taking the same expression as above. In Table 3 we display the absolute errors in the energy norm, that is ${\Vert grad\left({u}^{\prime}-{u}_{h}\right)\Vert}_{h}$ and in the mean-square norm ${\Vert {u}^{\prime}-{u}_{h}\Vert}_{h}$ , for increasing values of I, namely $I={2}^{m}$ for $m=1,2,3,4$ . Now we take as a reference $h=\text{\pi}/\left(8I\right)$ .

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 $1/{I}^{2}$ , as predicted. On the other hand here again the mean-square norm of the error function ${u}^{\prime}-{u}_{h}$ tends to decrease as $1/{I}^{3}$ . Likewise the two-dimensional case, Table 4 in turn shows a qualitative erosion of the solution errors if ${W}_{h}^{d}$ is replaced by ${V}_{h}$ 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 $\partial u/\partial n=0$ on $\Gamma $ (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 ${F}_{h}$ is modified, in such a way that boundary integrals for elements $T\in {\mathcal{S}}_{h}$ are shifted to the curved boundary portion sufficiently close to $\Gamma $ 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 ${\mathcal{S}}_{h}$ or in ${\mathcal{R}}_{h}$ , the inverse of an ${n}_{k}\times {n}_{k}$ matrix has to be computed, where ${n}_{k}$ is the dimension of ${P}_{k}\left(T\right)$ . 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.

References

[1] Brenner, S.C. and Scott, L.R. (2008) The Mathematical Theory of Finite Element Methods. Texts in Applied Mathematics, Vol. 15, Springer, Berlin.

https://doi.org/10.1007/978-0-387-75934-0

[2] Nitsche, J.A. (1972) On Dirichlet Problems Using Subspaces with Nearly Zero Boundary Conditions. In: Aziz, A.K., Ed., The Mathematical Foundations of the Finite Element Method with Applications to Partial Differential Equations, Academic Press, Cambridge, Massachusetts.

[3] Scott, L.R. (1973) Finite Element Techniques for Curved Boundaries. PhD Thesis, MIT, Cambridge, Massachusetts,.

[4] Ciarlet, P.G. (1978) The Finite Element Method for Elliptic Problems. North Holland, Amsterdam.

[5] Ruas, V. (2017) Optimal Simplex Finite-Element Approximations of Arbitrary Order in Curved Domains Circumventing the Isoparametric Technique. arXiv Numerical Analysis.

https://arxiv.org/abs/1701.00663

[6] Ruas, V. (2017) A Simple Alternative for Accurate Finite-Element Modeling in Curved Domains. Congrès Français de Mécanique, Lille, France.

https://cfm2017.sciencesconf.org/133073/document

[7] Nitsche, J.A. (1979) L∞-Convergence of Finite Element Galerkin Approximations for Parabolic Problems. ESAIM: Mathematical Modelling and Numerical Analysis, 13, 31-54.

[8] Ruas, V. (2017) Methods of Arbitrary Optimal Order with Tetrahedral Finite-Element Meshes Forming Polyhedral Approximations of Curved Domains. arXiv Numerical Analysis.

https://arxiv.org/abs/1706.08004