Anisotropic mesh refinement in polyhedral domains: error estimates with data in
We consider the homogeneous Dirichlet problem for the Laplace equation,
where is a polyhedral domain. Note that we could consider a more general elliptic equation of second order. But by a linear change of the independent variables the main part of the differential operator could be transformed to the Laplace operator in another polyhedral domain such that it is sufficient to consider the Laplace operator here.
The aim of the paper is to prove the discretization error estimate
for the finite element solution which is constructed by using piecewise linear and continuous functions on a family of appropriate finite element meshes . Note that we assume here not more than such that the -error estimate
follows by the Aubin–Nitsche method immediately. The generic constant may have different values on each occurrence.
If the solution of the boundary value problem was in then the finite element meshes could be chosen quasi-uniform, and the error estimates and would be standard. However, if the domain is non-convex, the solution will in general contain vertex and edge singularities, that means . In this case the convergence order is reduced in comparison with and when quasi-uniform meshes are used. As a remedy, we focus here on a priori anisotropic mesh grading techniques as they were investigated by Apel and Nicaise in . In comparison with isotropic local mesh refinement, the use of anisotropic elements avoids an unnecessary refinement along the edges.
The estimate is in general proven by using the Céa lemma (or the best approximation property of the finite element method),
and by proving an interpolation error estimate as an upper bound for the right-hand side of . The particular difficulty is that when the Lagrange interpolant is used together with anisotropic mesh grading, then the local interpolation error estimate
does not hold for but only for , see . Hence the classical proof of a finite element error estimate via
does not work. This problem was overcome by Apel and Nicaise, , by using and related estimates in weighted spaces, as well as the Hölder inequality for the prize that with has to be assumed in problem . Hence estimate cannot be proved in this way.
For prismatic domains and tensor product type meshes the problem was overcome in  by proving local estimates for a certain quasi-interpolation operator. This work cannot be easily extended to general polyhedral domains since the orthogonality of certain edges of the elements was used there. The aim of the current paper is to construct a quasi-interpolation operator such that the error estimate
can be proved for the anisotropic meshes introduced in .
Quasi-interpolants were introduced by Clément . The idea is to replace nodal values by certain averaged values such that non-smooth functions can be interpolated. This original idea has been modified by many authors since then. The contribution by Scott and Zhang  was most influential to our work.
The plan of the paper is as follows. In Section ? we introduce notation, recall regularity results for the solution of and describe the finite element discretization. The main results are proved in Section ?. The paper continues with numerical results in Section ? and ends with two sections where we describe applications which motivated us to improve the approximation result from , , to . The first one is a discretization of a distributed optimal control problem with as the state equation. The second application consists in a simpler proof of the discrete compactness property for edge elements of any order on this kind of finite element meshes.
We finish this introduction by commenting on related work. The idea to treat singularities due to a non-smooth boundary by using graded finite element meshes is old. The two-dimensional case was investigated by Oganesyan and Rukhovets , Babu ska , Raugel , and Schatz and Wahlbin . In three dimensions we can distinguish isotropic mesh grading, see the papers by Apel and Heinrich  and Apel, Sändig, and Whiteman , and anisotropic mesh grading, see the already mentioned papers  for the special case of prismatic domains, and  for general polyhedral domains. This work has been extended by Băcuţă, Nistor, and Zikatanov  to higher order finite element approximations where naturally higher regularity of the right-hand side has to be assumed. Boundary element methods with anisotropic, graded meshes have been considered by von Petersdorff and Stephan . The main alternative to mesh grading is augmenting the finite element space with singular functions, see for example Strang and Fix , Blum and Dobrowolski , or Assous, Ciarlet Jr., and Segré  for various variants. It works well in two dimensions where the coefficient in front of the singular function is constant. In the case of edge singularities this coefficient is a function which can be approximated, see Beagles and Whiteman , or it can be treated by Fourier analysis, see Lubuma and Nicaise .
2Notation, regularity, discretization
It is well known that the solution of the boundary value problem contains edge and vertex singularities which are characterized by singular exponents. For each edge , the corresponding leading (smallest) singular exponent is simply defined by where is the interior dihedral angle at the edge . For vertices of , the leading singular exponent has to be computed via the eigenvalue problem of the Laplace-Beltrami operator on the intersection of and the unit sphere centered at . Note that and . A vertex or an edge will be called singular if or , respectively. We exclude the case that is a singular exponent of any vertex. For a detailed discussion of edge and vertex singularities we refer to .
As in  we subdivide the domain into a finite number of disjoint tetrahedral subdomains, subsequently called macro-elements,
We assume that each contains at most one singular edge and at most one singular vertex. In the case that contains both a singular edge and a singular vertex, that vertex is contained in that edge. Note that the edges of are considered to have length. For , the closures of the macroelements and may be disjoint or they intersect defining a coupling face, or a coupling edge, or a coupling node. Denote by , and the sets of coupling faces, edges and nodes, respectively.
For the description of the regularity of the solution of , we set if the macro-element contains the singular vertex of . If does not contain any singular vertex we set . Moreover, we set if contains the singular edge of , otherwise we set . Furthermore, we define in each macro-element a Cartesian coordinate system such that the singular vertex, if existing, is located in the origin, and the singular edge, if existing, is contained in the -axis. We also introduce by
the distance to the -axis, the distance to the origin, the angular distance from the -axis, respectively.
For and we define the weighted Sobolev space
Here, we have used the standard multi-index notation to describe partial derivatives, and we have omitted the index in and for simplicity.
Following  we consider a triangulation of ,
made up of tetrahedra which match the initial partition: if then . Four cases are considered:
If does neither contain a singular edge nor a singular vertex then is assumed to be isotropic and quasi-uniform with element size , see Figure ?, top left.
If contains a singular vertex but no singular edges then is isotropic and has a singular vertex refinement, i.e., the mesh is graded towards the singular vertex with a grading parameter . This can be achieved by using a coordinate transformation of the vertices from Case 1, see Figure ?, top right.
If contains a singular edge but no singular vertices then is anisotropically graded towards the singular edge. The grading parameter is . To this end, we introduce a family of planes transversal to the singular edge and containing the opposite one. These planes split the macro element into strips and contain all nodes. In the planes the position of the nodes is achieved by applying a coordinate transformation to a uniform triangulation, see Figure ?, bottom left.
If contains both a singular vertex and a singular edge then is graded towards the singular edge with grading parameter and towards the singular vertex with grading parameter . The mesh is topologically equivalent to the mesh of Case 3 but the planes of do not divide the singular edge equidistantly but with a grading towards the singular vertex.
We point out that anisotropic elements can appear only in Cases 3 and 4, for which contains needle elements near the singular edge and flat elements near the opposite one, see Figure ?. We further observe that if is of type 3 or 4, the elements in do not intersect any plane of .
For each element we introduce its lengths and as follows. Let be the diameter of . If with of type 1 or 2, then . If with of type 3 or 4 then is the length of the edge of parallel to the singular edge, and where and are the edges of intersecting and each one of them is contained in some plane of .
By classical regularity theory, the solution of the boundary value problem is continuous, see e.g. , such that the Lagrange interpolant with respect to the subdivision is well defined. We consider the decomposition
It follows that the restriction has the same smoothness properties as , see Theorem ?. Furthermore, vanishes in coupling nodes and on singular edges. We construct now an interpolant which also vanishes on these nodes such that can be used to estimate the discretization error via .
To this end, let , , and be the set of all nodes of , the set of all the interior nodes, the set of coupling nodes, and the set of nodes which belong to some singular edge, respectively. The terminal points of the singular edges are included in . The piecewise linear nodal basis on is denoted by . We associate (as specified below) with each an edge with as an endpoint. Note that since with . Hence the operator with
is well defined when is the -projection operator onto the space of polynomials of degree less than or equal to one. Note that vanishes on coupling nodes and on singular edges by construction. In order to impose the boundary conditions and to be able to prove interpolation error estimates we need to select the edges in an appropriate way, compare the illustration in Figure ?.
First, we demand that
for each node , and belong to the same macroelement.
This requires in particular the following restrictions.
If lays on a boundary or coupling face, then is contained in that face.
If lays on a coupling edge, then is contained in that coupling edge.
Note that these requirements made the treatment of the coupling nodes via the interpolation on the initial necessary. Note further that this construction leads to a preservation of the homogeneous Dirichlet boundary condition.
In order to prove the stability of in the anisotropic refinement regions we also require:
If is a vertex of a tetrahedron contained in a macroelement of types 3 or 4, then is an edge contained on some plane of .
If and belong to a macroelement of types 3 or 4 and have the same orthogonal projection onto the -plane, then the same holds for and .
In order to estimate the interpolation error we need to define for each a set which should satisfy the following assumptions.
The set is a union of elements of (plus some faces) and in particular .
The set is an open connected domain, and as small as possible.
We have for all nodes of .
If , then .
If with of type 3 of 4, then is a prism where the top and bottom faces are contained in two planes of (and so they are not parallel) and the other faces are parallel to the singular edge.
The following properties follow from the definitions of the edges and the sets .
Let be contained in a macroelement of type or . If intersects two planes and of , then intersects exactly the same planes and .
If the node , , belongs to a coupling face, that means that there exist tetrahedra and with and , then but .
If is an isotropic element then all the elements in are also isotropic and of size of the same order.
The second point is essential for our proof of the approximation properties. It was the target for which we made the construction as it is.
The aim of this section is to derive error estimates for our discretization. They are based on local interpolation error estimates for our interpolant . For proving these estimates we have to distinguish several cases, see also Figure ? for an illustration:
is an isotropic element without coupling node, has full regularity,
is an isotropic element with coupling node, has full regularity,
is an isotropic element with coupling node, has reduced regularity,
is an anisotropic flat element without coupling node, has full regularity,
is an anisotropic flat element with coupling node, has full regularity,
is an anisotropic needle element without node on the singular edge, has full regularity,
is an anisotropic needle element with node on the singular edge, has reduced regularity.
In Lemma ? we present the general approach for the proof of the local interpolation error estimate by considering isotropic elements with and without coupling nodes (cases 1 and 2). We proceed with Lemmas ? where we introduce for isotropic elements how to cope with the weighted norms in the case of reduced regularity (case 3). The interpolated function is only from a weighted Sobolev space but we will see that this even simplifies some parts of the proof.
For anisotropic elements the use of an inverse inequality (as was done in the previous lemmas) has to be avoided; instead we use the structure of the meshes in the macroelements of types 3 and 4. We start with a stability estimate of which allows immediately the treatment of anisotropic flat elements (cases 4 and 5) in Lemma ?. Then we prove stability estimates for the remaining derivatives and continue with the interpolation error estimates for needle elements. Lemma ? is devoted to case 6, and Lemma ? to case 7.
All these local estimates can then be combined to prove the global interpolation error estimate, see Theorem ?, and the finite element error estimates, see Corollary ?.
where we denote by the set of nodes of without the coupling nodes. Note that
compare . (By some calculation one can even specify that .) With , the direct computation
the trace theorem
and we obtain
If does not contain a node we find that for all such that we get by using the triangle inequality and the stability estimate
We use now a Deny–Lions type argument (see e.g. ) and conclude estimate .
In the case when contains a node we do not have the property that for all but we can use that . Let be an edge contained in having as an endpoint, and let be the Lagrange basis function associated with . (Note that we deal here with nodes which are not used in the definition of . Therefore we can assume that is local in .) Consequently, we have with the previous argument that
Let be the linear Lagrange interpolation of on . Since is linear, we have have . From this fact and using – as in the derivation of (here with the specific instead of since ), we have
where we used standard estimates for the Lagrange interpolant in the last step. With and the triangle inequality we conclude estimate also in this case.
We start as in the proof of Lemma ? but use the sharper trace theorem
With , , , and we obtain
and hence via the triangle inequality
For the first two terms we just use that , hence , to get
To estimate the third term we use the Cauchy–Schwarz inequality and again , to obtain for
where is obtained by executing the integration and using that . All these estimates imply estimate .
In order to prove interpolation error estimates for the anisotropic elements we derive stability estimates for where we avoid the use of the inverse inequality. Let and be a Cartesian coordinate system with the -direction parallel to the singular edge of . We will estimate separately the -norm of the derivatives of .
Let be an anisotropic element with the characteristic lengths and . We will not use that , , in the next lemma in order to use this estimate both for the needle and the flat elements.
We observe that has an edge parallel to the singular edge, and so, parallel to the -axis. Since is linear on , we have . If is contained on the singular edge, then since and we are done. Now, consider the case that is not contained in a singular edge and denote its endpoints by and such that and . Then we have
We observe now that by our assumptions and have the same projection into the -plane and hence form two opposite edges of a plane quadrilateral which is parallel to the -axis and which we will denote by . We note further that and can be considered as the same function defined on and . With this insight we obtain
We integrate this estimate over , apply the standard trace theorem
and obtain the desired estimate.
We are now prepared to estimate the interpolation error for the flat elements occurring far away from the singular edge in cases 3 and 4.
The proof for can be done on the basis of Lemma ?. Assume for the moment that the element does not contain a coupling node. Similar to the proof of Lemma ? we obtain for any
We choose now such that the constant satisfies and such that we can conclude by using the Poincaré–Friedrichs inequality (or again a Deny–Lions type argument)
Note that the polynomial can be chosen such that it vanishes in three nodes of . It is completely described by choosing the appropriate value at one endpoint of the edge of which is parallel to the -axis. Since a possible coupling node is not an endpoint of this edge, the argument above can also be used in the case of coupling nodes.
For the other directions we can proceed as in the proof of Lemma ?. In the case of coupling nodes the interpolation error estimate is used there which does not hold for anisotropic elements. However, the estimate , , does hold, see for example .
It remains to prove interpolation error estimates for needle elements such that we will assume for the next lemmas.
For each node we denote by the top or bottom face of the prismatic domain such that . Observe that we have for all . Observe further that is isotropic with diameter of order and recall the standard trace inequality
for all . We need also the trace inequality
which can be proved by using Lemma ? from page and the facts that is a union of prisms, and is a face of .
Let be one of the short edges of and denote its endpoints by and . We use the same notation for the direction of this edge in order to denote by the directional derivative. In the following we first estimate . After that, the desired estimates easily follow as we will show.
Notice that if we have , and if then . For all we have (and here we use that the element does not contain a node )
From the trace inequality we have for each
Since the definition of implies , we have
Now we choose as the average of on and use a Poincaré type inequality on to get
Therefore we arrive at
where we used again the trace inequality .
Now, let and be two different short edges (edge vectors) of such that the determinant of the matrix made up of , and as columns is greater than a constant depending only the maximum angle of . Note that this is possible due to the maximal angle condition, see . Then, if the canonical vector , , is expressed as
it follows that , and are bounded by above by a constant depending only on the maximum angle condition. Since
we obtain from with and , Lemma ?, and recalling that .
For each node of we select one short edge with an endpoint at and contained in the same macroelement as such that we can apply Lemma ?. We have for
Now we deal with which is first estimated by