Startseite A Priori Analysis of an Anisotropic Finite Element Method for Elliptic Equations in Polyhedral Domains
Artikel Öffentlich zugänglich

A Priori Analysis of an Anisotropic Finite Element Method for Elliptic Equations in Polyhedral Domains

  • Hengguang Li EMAIL logo und Serge Nicaise
Veröffentlicht/Copyright: 18. Februar 2020

Abstract

Consider the Poisson equation in a polyhedral domain with mixed boundary conditions. We establish new regularity results for the solution with possible vertex and edge singularities with interior data in usual Sobolev spaces H σ with σ [ 0 , 1 ) . We propose anisotropic finite element algorithms approximating the singular solution in the optimal convergence rate. In particular, our numerical method involves anisotropic graded meshes with fewer geometric constraints but lacking the maximum angle condition. Optimal convergence on such meshes usually requires the pure Dirichlet boundary condition. Thus, a by-product of our result is to extend the application of these anisotropic meshes to broader practical computations with the price to have “smoother” interior data. Numerical tests validate the theoretical analysis.

MSC 2010: 65N30; 35B65; 35J25

1 Introduction

Consider the elliptic problem associated with the Laplace operator in a bounded polyhedral domain Ω 3 with the mixed boundary condition

(1.1) { - Δ u = f in Ω , u = 0 on Γ Dir , n u = 0 on Γ Neu ,

where Γ Dir and Γ Neu are open subsets of the boundary Ω such that Γ Dir ¯ Γ Neu ¯ = Ω . For simplicity, we suppose that each face of Ω is included either in Γ Dir or in Γ Neu and Γ Dir . The solution of equation (1.1) is uniquely defined in H Γ Dir 1 ( Ω ) (see (2.1)) for f ( H Γ Dir 1 ( Ω ) ) (cf. [15, 23]). The solution regularity, however, is determined by the smoothness of the given function, the geometry of the domain and the boundary conditions. Let us refer to the non-smooth boundary points and the points where the boundary condition changes as singular points. Then, even for a smooth function f, the solution may possess singularities in high-order Sobolev spaces near the singular points [13, 16, 18, 19]. These singularities, often being the main theoretical concern, can also severely deteriorate the efficacy of the numerical approximation.

For equation (1.1), the singular points are in fact the non-smooth boundary points (namely, vertices and edges), provided that each face is either in Γ Dir or in Γ Neu . Then, given a sufficiently smooth function f, there are two types of solution singularities near the singular points in Ω ¯ : the vertex singularity and the anisotropic edge singularity. For such singularities, anisotropic meshes are usually designed to improve the effectiveness of the finite element method (FEM). This is different from the isotropic graded meshes in two-dimensional polygonal domains, where only corner (vertex) singularities need special numerical treatment. The development of optimal FEMs for elliptic equations in polyhedral domains is a technically challenging task due to the combination of different types of singularities and due to the complexity in the three-dimensional geometry. Meanwhile, the error analysis for the numerical scheme often demands specific anisotropic regularity estimates. Compared with estimates in isotropic Sobolev spaces, such anisotropic regularity results are limited and less known, most of which are for pure Dirichlet problems. Consequently, the developments of effective mesh algorithms are extensively centered around pure Dirichlet problems.

The existing mesh algorithms for polyhedral domains usually require restrictive geometric conditions on the mesh and on the domain. For example, the mesh in [2, 14] is based on the method of dyadic partitioning. These meshes are isotropic and optimal only for weaker singular solutions. The mesh in [1, 3, 4] is based on a coordinate transformation from a quasi-uniform mesh. It is anisotropic along the edges and requires confining angle conditions for the simplex. The mesh in [7, 8] is also anisotropic and leads to optimal convergence rate. The algorithm, however, requires extra steps for prism refinements to maintain the angle condition in the simplex. There are also tensor-product anisotropic meshes based on 2D graded meshes [5, 24] that are usually effective on a domain with simple geometry. Recently, a new anisotropic FEM has emerged [21, 22] based on explicit recursive refinements. With fewer geometric requirements on the simplex and on the domain, this algorithm leads to conforming triangulations, which however violates the maximum angle condition in simplexes [6, 20]. Nevertheless, it was shown that, for the pure Dirichlet problem, the solution has extra regularity in the edge direction to compensate for the lack of mesh shape regularity, and this algorithm gives rise to optimal FEMs approximating the anisotropic singular solutions. Equations with mixed boundary conditions can possess solutions with a singular structure different from that in pure Dirichlet problems. Especially, near an edge or a vertex that is surrounded by Neumann faces, the solution does not vanish and therefore does not belong to the same Sobolev space as in the Dirichlet case. The rigorous theoretical and numerical justification of anisotropic algorithms for problems with mixed boundary conditions, which occurs often in practical computations, remains an open investigation.

In this paper, we extend the application of anisotropic algorithms to problems with mixed boundary conditions by developing new finite element algorithms and new regularity results for equation (1.1). In particular, we study the singular expansion of the solution near singular points surrounded by Neumann faces. It turns out that part of the singular expansion resembles the singularity in the Dirichlet problem and therefore belong to a similar weighted space. For the rest of the singular expansion, a series of estimates on its fundamental properties shall reveal its directional regularities. We summarize our findings by establishing new regularity results in Theorems 6.1 and 6.2, Lemma 6.3 and Corollary 6.5 in different parts of the domain. Then we propose an optimal finite element algorithm (Algorithm 3.4) and validate it based on interpolation error analysis in anisotropic weighted spaces.

The paper is organized as follows. In Section 2, we introduce necessary notation regarding the finite element approximation of equation (1.1). We also define a domain decomposition according to the distance to the singular points. In Section 3, we first review the anisotropic mesh developed in [21]. Then we propose the anisotropic FEM for equation (1.1) with the mixed boundary condition. In Section 4, we study the regularity of the equation in a dihedron, which shall lead to the local regularity estimates near an open edge. In Section 5, we investigate the regularity of the equation in an infinite cone, which shall lead to the local regularity estimates near a vertex of the domain. In Section 6, we present the regularity results for the solution in the domain. In Section 7, we include detailed interpolation error analysis for the anisotropic finite element algorithm in weighted spaces. These optimal interpolation error estimates in turn lead to the conclusion that the proposed FEMs obtain the optimal convergence rate approximating the target problem. Numerical tests are implemented in a polyhedral prism domain for different mixed boundary conditions, and the results are reported in Section 8. These numerical results are in agreement with our theoretical prediction and hence validate our method.

Throughout the text below, we adopt the bold notation for vector fields. Let T be a triangle (resp. tetrahedron) with vertices a , b , c (resp. a , b , c , d ). Then we denote T by its vertices: 3 a b c for the triangle and 4 a b c d for the tetrahedron, where the sup-index implies the number of vertices for T. We denote by ab the open line segment with endpoints a and b. By a b (resp. a b ) we mean that there exists a constant C > 0 independent of a and b such that C - 1 a b C a (resp. a C b ). The generic constant C > 0 in our estimates may be different at different occurrences. It will depend on the computational domain, but not on the functions involved or the mesh level in the finite element algorithms. In addition, both of the terms are used to represent the same directional derivative: 1 = x , 2 = y , and 3 = z . For a bounded domain D (or its boundary), the usual norm and semi-norm of H s ( D ) ( s 0 ) are denoted by s , D and | | s , D , respectively. For s = 0 , we will drop the index 0, and for D = Ω the index Ω. For two positive parameters s and ρ, we finally introduce the norm s , D , ρ on H s ( D ) (see [13, Definition AA.17] for instance) defined by

u s , D , ρ = ( ρ 2 s u D 2 + | u | s , D 2 ) 1 2 for all u H s ( D ) ,

that is equivalent to the usual norm s , D (with constants of equivalence depending on ρ) if D has a Lipschitz boundary.

2 Preliminaries

In this section, we introduce the notation and recall some existing results regarding the solution of (1.1).

2.1 The Finite Element Approximation

By a polyhedral domain Ω 3 we mean a bounded domain with a Lipschitz boundary Ω made of plane faces (i.e., its boundary is a finite union of polygons). Thus, the boundary of Ω is smooth, except at the vertex points and along the edges. In a neighborhood of a vertex c, Ω coincides with a three-dimensional cone, while near an interior point of an edge e, Ω resembles a dihedral angle.

For a bounded domain 𝒪 of 3 , let H m ( 𝒪 ) , m 0 , be the usual Sobolev space that consists of functions defined in 𝒪 whose kth derivatives are square-integrable for 0 k m (hence L 2 ( 𝒪 ) : = H 0 ( 𝒪 ) ). Let H loc m ( Ω ) : = { v , v H m ( G ) for any open subset G with compact closure G ¯ Ω } . The trace operator from H 1 ( Ω ) into H 1 2 ( Ω ) will be denoted by γ. We define

(2.1) H Γ Dir 1 ( Ω ) : = { u H 1 ( Ω ) , γ u = 0 on Γ Dir } ,

which is clearly a closed subspace of H 1 ( Ω ) .

Then, for f L 2 ( Ω ) , the variational solution u H Γ Dir 1 ( Ω ) of problem (1.1) is defined by

(2.2) a ( u , v ) : = Ω u v d x = ( f , v ) : = Ω f v d x for all v H Γ Dir 1 ( Ω ) .

Let 𝒯 n be a triangulation of Ω with tetrahedra. Let S n H Γ Dir 1 ( Ω ) be the linear Lagrange finite element space associated with 𝒯 n . Then the finite element solution u n S n for equation (1.1) is given by

(2.3) a ( u n , v n ) = ( f , v n ) for all v n S n .

Remark 2.1.

By Poincaré’s inequality, the bilinear form a ( , ) is both continuous and coercive on

V : = H Γ Dir 1 ( Ω ) .

Then, by Céa’s lemma [9, 11], u n is the best approximation from S n in V, u - u n V inf v n S n u - v n V . It is well known that the solution u may not belong to H 2 ( Ω ) due to the presence of the non-smooth points (vertices and edges) on the boundary. On a standard quasi-uniform triangulation 𝒯 n , the limited regularity of u in the Sobolev space can result in a sub-optimal convergence rate for the finite element approximation. Namely, u - u n H 1 ( Ω ) C h s u H s + 1 ( Ω ) , where h is the mesh size in 𝒯 n and 0 < s < 1 depends on the geometry of the domain.

For equation (1.1), there are two types of singularities in the solution that may affect the convergence of numerical methods. The vertex singularity appears in the neighborhood of a vertex and concentrates at the vertex. The edge singularity occurs in the neighborhood of an edge; it is however anisotropic in the sense that the solution is smoother in the direction along the edge than toward the edge. Consequently, anisotropic graded meshes are frequently applied to improve the convergence of the finite element solution.

2.2 The Domain and the Weighted Sobolev Space

We denote by the finite set of open edges and by 𝒞 the finite set of vertices of Ω. We also denote by c the set of edges joining at c 𝒞 and by 𝒞 e 𝒞 the set of endpoints of e . We say an edge e is a Dirichlet (Neumann) edge if the Dirichlet (Neumann) boundary conditions are imposed on both adjacent faces of e. We say e is a DN edge if the Dirichlet condition is imposed on one adjacent face of e and the Neumann condition is on the other. Let ω e be the opening angle between the two adjacent faces of e. For each e , define

(2.4) ν e = { π ω e if e is a Dirichlet edge or a Neumann edge , π 2 ω e if e is a DN edge .

The edge e is called singular if ν e < 1 ; otherwise, it is called regular. Denote by Γ c the cone that coincides with the domain Ω at c 𝒞 . Let ν c be the first positive eigenvalue of the Laplace–Beltrami operator on the intersection of Γ c with the unit sphere with boundary conditions inherited from equation (1.1). Then, if - 1 2 + ( ν c + 1 4 ) 1 2 < 1 2 , c is called singular; it is regular otherwise. For e and c 𝒞 , we set

λ e = { ν e if e is singular , otherwise ;
λ c = { - 1 2 + ( ν c + 1 4 ) 1 2 if c is singular , otherwise .

To better describe the singular behavior of the solution near the non-smooth points, we further define the distance functions. For any c 𝒞 (resp. e ), we define R c ( x ) (resp. r e ( x ) ) to be the distance from x Ω to c (resp. to e). We further define θ c , e ( x ) : = r e ( x ) R c ( x ) as the angular distance from x to the edge of e near c. Then, for any vertex c 𝒞 and edge e , as in [10, 22], we define the following subsets of Ω:

(2.5) { 𝒱 c = { x Ω , R c ( x ) < ε } , 𝒱 c e = { x 𝒱 c , θ c , e ( x ) < ε } , 𝒱 c 0 = { x 𝒱 c , θ c , e ( x ) ε for all e c } , 𝒱 e 0 = { x Ω , R c ( x ) ε , θ c , e ( x ) < ε for all c 𝒞 e }

with ε > 0 small enough such that all these sets are disjoint for different vertices c and edges e. We further define

(2.6) 𝒱 0 = Ω ( ( c 𝒞 𝒱 c ) ( e 𝒱 e 0 ) ) .

It is clear that the subsets in (2.5) are neighborhoods of different non-smoothness points on the boundary. In the neighborhoods 𝒱 e 0 and 𝒱 c e , we choose a local Cartesian coordinate system in which the edge e lies on the z-axis. Let α = ( α 1 , α 2 ) consist of the first two entries of the multi-index α = ( α 1 , α 2 , α 3 ) 0 3 . Therefore, in 𝒱 e 0 and 𝒱 c e , α = x α 1 y α 2 is a partial derivative in a direction perpendicular to the edge e. Meanwhile, we define | α | : = α 1 + α 2 + α 3 and | α | : = α 1 + α 2 .

We shall need the following weighted Sobolev space of Kondratiev’s type. Let 𝒪 be a subset of 3 such that 0 belongs to its boundary. Then, for any β , k 0 , we define the space

V β k ( 𝒪 ) = { v L loc 2 ( 𝒪 ) r β + | α | - k D α v L 2 ( 𝒪 ) for all | α | k } ,

where r is the distance to 0. Define λ : = max e ( 0 , 1 - λ e ) . We also assume the given data in equation (1.1) to satisfy f H σ ( Ω ) with

(2.7) σ { ( λ , 1 ) if λ > 0 , [ 0 , 1 ) if λ = 0 .

3 Anisotropic Finite Element Algorithms

In this section, we propose new anisotropic FEMs approximating equation (1.1). In particular, we give explicit values for the associated parameters in the algorithm, with which we shall prove the proposed method achieves the optimal rate of convergence, even when the solution is singular.

3.1 Anisotropic Algorithms

Recall the vertex set 𝒞 and the edge set . Following the notation in [21], we first classify tetrahedra in the triangulation of Ω.

Definition 3.1 (Tetrahedron Types).

Let T be a tetrahedron. If an edge e T of T lies on e , we call e T a marked edge. Let c T be a vertex of T. If c T 𝒞 , or if c T is an interior point of an edge e and c T = e T ¯ , we call c T a marked vertex. Let 𝒯 be a tetrahedral triangulation of Ω such that (I) each tetrahedron contains at most one marked vertex and at most one marked edge, (II) if a tetrahedron contains both a marked vertex and a marked edge, the marked vertex is an endpoint of the marked edge. Let 𝒮 = 𝒞 . Then, for each tetrahedron T 𝒯 , according to its relation with 𝒮 , there are five possible types.

  1. o-tetrahedron: T ¯ 𝒮 = .

  2. v-tetrahedron: T ¯ 𝒮 = c 𝒞 .

  3. v e -tetrahedron: T ¯ 𝒮 is an interior point of an edge e .

  4. e-tetrahedron: T ¯ 𝒮 is a marked edge, but contains no vertex in 𝒞 .

  5. ev-tetrahedron: T ¯ 𝒮 contains a marked edge and a marked vertex.

Note that different types of tetrahedra in Definition 3.1 are associated to different sub-regions of Ω in (2.5) and (2.6). In addition, we recall the following anisotropic mesh algorithm [21].

Algorithm 3.2 (Anisotropic Refinement).

Let 𝒯 be a triangulation of Ω as in Definition 3.1. To each c 𝒞 (resp. e ), we associate a grading parameter κ c (resp. κ e ) ( 0 , 1 2 ] . Let T = 4 x 0 x 1 x 2 x 3 𝒯 be a tetrahedron with vertices x 0 , x 1 , x 2 , x 3 such that x 0 is the marked vertex if T is a v-, v e - or ev-tetrahedron, and x 0 x 1 is the marked edge if T is an e- or ev-tetrahedron. Let 𝜿 be the collection of the parameters κ c and κ e for all c 𝒞 and e . Then the refinement, denoted by 𝜿 ( 𝒯 ) , proceeds as follows. We first generate new nodes x k l , 0 k < l 3 , on each edge x k x l of T, based on its type.

  1. o-tetrahedron: x k l = x k + x l 2 .

  2. v-tetrahedron: Suppose x 0 = c 𝒞 . Define κ = κ e c : = min e c ( κ c , κ e ) . Then x k l = x k + x l 2 for 1 k < l 3 ; x 0 l = ( 1 - κ ) x 0 + κ x l for 1 l 3 .

  3. v e -tetrahedron: Suppose x 0 is an interior point of e . Let κ = κ e . Then x k l = x k + x l 2 for 1 k < l 3 ; x 0 l = ( 1 - κ ) x 0 + κ x l for 1 l 3 .

  4. e-tetrahedron: Suppose x 0 x 1 e . Then x k l = ( 1 - κ e ) x k + κ e x l for 0 k 1 and 2 l 3 ; x 01 = x 0 + x 1 2 , x 23 = x 2 + x 3 2 .

  5. ev-tetrahedron: Suppose x 0 = c 𝒞 and x 0 x 1 e c . Define κ e c : = min e c ( κ c , κ e ) . Then, for 2 l 3 , x 0 l = ( 1 - κ e c ) x 0 + κ e c x l and x 1 l = ( 1 - κ e ) x 1 + κ e x l ; x 01 = ( 1 - κ c ) x 0 + κ c x 1 , x 23 = x 2 + x 3 2 .

Connecting these nodes x k l on all the faces of T, we obtain four sub-tetrahedra and one octahedron. The octahedron then is cut into four tetrahedra using x 13 as the common vertex. Therefore, after one refinement, we obtain eight sub-tetrahedra for each T 𝒯 denoted by their vertices:

4 x 0 x 01 x 02 x 03 , 4 x 01 x 1 x 12 x 13 , 4 x 02 x 12 x 2 x 23 , 4 x 03 x 13 x 23 x 3 ,
4 x 01 x 02 x 03 x 13 , 4 x 01 x 02 x 12 x 13 , 4 x 02 x 03 x 13 x 23 , 4 x 02 x 12 x 13 x 23 .

See Figure 1 for different types of decompositions. Given an initial mesh 𝒯 0 satisfying the condition in Definition 3.1, the associated family of anisotropic meshes { 𝒯 n , n 0 } is defined recursively, 𝒯 n = 𝜿 ( 𝒯 n - 1 ) . See Figure 2 for example.

Figure 1 
                  Refinements of a tetrahedron 
                        
                           
                              
                                 
                                    △
                                    4
                                 
                                 ⁢
                                 
                                    x
                                    0
                                 
                                 ⁢
                                 
                                    x
                                    1
                                 
                                 ⁢
                                 
                                    x
                                    2
                                 
                                 ⁢
                                 
                                    x
                                    3
                                 
                              
                           
                           
                           {\triangle^{4}x_{0}x_{1}x_{2}x_{3}}
                        
                     , top (left to right): o-tetrahedron, v- or 
                        
                           
                              
                                 v
                                 e
                              
                           
                           
                           {v_{e}}
                        
                     -tetrahedron, e-tetrahedron; bottom (left to right): two ev-tetrahedra with 
                        
                           
                              
                                 
                                    κ
                                    
                                       e
                                       ⁢
                                       c
                                    
                                 
                                 =
                                 
                                    κ
                                    e
                                 
                              
                           
                           
                           {\kappa_{ec}=\kappa_{e}}
                        
                      and 
                        
                           
                              
                                 
                                    κ
                                    
                                       e
                                       ⁢
                                       c
                                    
                                 
                                 =
                                 
                                    κ
                                    c
                                 
                              
                           
                           
                           {\kappa_{ec}=\kappa_{c}}
                        
                     .
Figure 1 
                  Refinements of a tetrahedron 
                        
                           
                              
                                 
                                    △
                                    4
                                 
                                 ⁢
                                 
                                    x
                                    0
                                 
                                 ⁢
                                 
                                    x
                                    1
                                 
                                 ⁢
                                 
                                    x
                                    2
                                 
                                 ⁢
                                 
                                    x
                                    3
                                 
                              
                           
                           
                           {\triangle^{4}x_{0}x_{1}x_{2}x_{3}}
                        
                     , top (left to right): o-tetrahedron, v- or 
                        
                           
                              
                                 v
                                 e
                              
                           
                           
                           {v_{e}}
                        
                     -tetrahedron, e-tetrahedron; bottom (left to right): two ev-tetrahedra with 
                        
                           
                              
                                 
                                    κ
                                    
                                       e
                                       ⁢
                                       c
                                    
                                 
                                 =
                                 
                                    κ
                                    e
                                 
                              
                           
                           
                           {\kappa_{ec}=\kappa_{e}}
                        
                      and 
                        
                           
                              
                                 
                                    κ
                                    
                                       e
                                       ⁢
                                       c
                                    
                                 
                                 =
                                 
                                    κ
                                    c
                                 
                              
                           
                           
                           {\kappa_{ec}=\kappa_{c}}
                        
                     .
Figure 1 
                  Refinements of a tetrahedron 
                        
                           
                              
                                 
                                    △
                                    4
                                 
                                 ⁢
                                 
                                    x
                                    0
                                 
                                 ⁢
                                 
                                    x
                                    1
                                 
                                 ⁢
                                 
                                    x
                                    2
                                 
                                 ⁢
                                 
                                    x
                                    3
                                 
                              
                           
                           
                           {\triangle^{4}x_{0}x_{1}x_{2}x_{3}}
                        
                     , top (left to right): o-tetrahedron, v- or 
                        
                           
                              
                                 v
                                 e
                              
                           
                           
                           {v_{e}}
                        
                     -tetrahedron, e-tetrahedron; bottom (left to right): two ev-tetrahedra with 
                        
                           
                              
                                 
                                    κ
                                    
                                       e
                                       ⁢
                                       c
                                    
                                 
                                 =
                                 
                                    κ
                                    e
                                 
                              
                           
                           
                           {\kappa_{ec}=\kappa_{e}}
                        
                      and 
                        
                           
                              
                                 
                                    κ
                                    
                                       e
                                       ⁢
                                       c
                                    
                                 
                                 =
                                 
                                    κ
                                    c
                                 
                              
                           
                           
                           {\kappa_{ec}=\kappa_{c}}
                        
                     .
Figure 1 
                  Refinements of a tetrahedron 
                        
                           
                              
                                 
                                    △
                                    4
                                 
                                 ⁢
                                 
                                    x
                                    0
                                 
                                 ⁢
                                 
                                    x
                                    1
                                 
                                 ⁢
                                 
                                    x
                                    2
                                 
                                 ⁢
                                 
                                    x
                                    3
                                 
                              
                           
                           
                           {\triangle^{4}x_{0}x_{1}x_{2}x_{3}}
                        
                     , top (left to right): o-tetrahedron, v- or 
                        
                           
                              
                                 v
                                 e
                              
                           
                           
                           {v_{e}}
                        
                     -tetrahedron, e-tetrahedron; bottom (left to right): two ev-tetrahedra with 
                        
                           
                              
                                 
                                    κ
                                    
                                       e
                                       ⁢
                                       c
                                    
                                 
                                 =
                                 
                                    κ
                                    e
                                 
                              
                           
                           
                           {\kappa_{ec}=\kappa_{e}}
                        
                      and 
                        
                           
                              
                                 
                                    κ
                                    
                                       e
                                       ⁢
                                       c
                                    
                                 
                                 =
                                 
                                    κ
                                    c
                                 
                              
                           
                           
                           {\kappa_{ec}=\kappa_{c}}
                        
                     .
Figure 1 
                  Refinements of a tetrahedron 
                        
                           
                              
                                 
                                    △
                                    4
                                 
                                 ⁢
                                 
                                    x
                                    0
                                 
                                 ⁢
                                 
                                    x
                                    1
                                 
                                 ⁢
                                 
                                    x
                                    2
                                 
                                 ⁢
                                 
                                    x
                                    3
                                 
                              
                           
                           
                           {\triangle^{4}x_{0}x_{1}x_{2}x_{3}}
                        
                     , top (left to right): o-tetrahedron, v- or 
                        
                           
                              
                                 v
                                 e
                              
                           
                           
                           {v_{e}}
                        
                     -tetrahedron, e-tetrahedron; bottom (left to right): two ev-tetrahedra with 
                        
                           
                              
                                 
                                    κ
                                    
                                       e
                                       ⁢
                                       c
                                    
                                 
                                 =
                                 
                                    κ
                                    e
                                 
                              
                           
                           
                           {\kappa_{ec}=\kappa_{e}}
                        
                      and 
                        
                           
                              
                                 
                                    κ
                                    
                                       e
                                       ⁢
                                       c
                                    
                                 
                                 =
                                 
                                    κ
                                    c
                                 
                              
                           
                           
                           {\kappa_{ec}=\kappa_{c}}
                        
                     .
Figure 1

Refinements of a tetrahedron 4 x 0 x 1 x 2 x 3 , top (left to right): o-tetrahedron, v- or v e -tetrahedron, e-tetrahedron; bottom (left to right): two ev-tetrahedra with κ e c = κ e and κ e c = κ c .

Remark 3.3.

The anisotropic mesh in Algorithm 3.2 is explicitly determined by the grading parameters κ c and κ e that are associated to each vertex and edge. A smaller value of the parameter leads to a higher mesh density near the vertex or the edge, while the value κ c = κ e = 1 2 corresponds to a quasi-uniform refinement. In different regions of the domain, the resulting mesh may have different shape regularities. In 𝒱 0 , the mesh is isotropic and quasi-uniform. The local refinement for a v- or v e -tetrahedron in fact follows the same rule: the mesh is isotropic and graded toward the vertex x 0 based on the grading parameter κ associated to the vertex x 0 . In 𝒱 e 0 , the resulting mesh in general is anisotropic and graded toward the edge e . The refinement in 𝒱 e c depends on the parameters κ c and κ e , e c , which is also anisotropic, graded toward the edge e and the vertex c 𝒞 . We also mention that the mesh in 𝒱 e 0 and in 𝒱 e c does not satisfy the maximum angle condition [6, 20] if κ e < 1 2 , which can lead to a fair amount of difficulty in analysis. Nevertheless, it has been shown in [21, 22] that these anisotropic meshes are effective in approximating three-dimensional singular solutions provided the pure Dirichlet boundary condition is imposed. For mixed boundary conditions, the singular solution no longer belongs to the same space as in [21, 22]. The algorithm design and analysis is therefore more technically involved.

Now, we proceed to propose our finite element algorithm for equation (1.1) with f H σ ( Ω ) .

Algorithm 3.4 (Anisotropic Finite Element Algorithm).

Let 𝒯 0 be the initial triangulation of Ω that satisfy the condition in Definition 3.1. Then each parameter κ c ( resp. κ e ) ( 0 , 1 2 ] is uniquely determined by a new parameter a c ( resp. a e ) ( 0 , 1 ] such that

(3.1) κ c = 2 - 1 a c and κ e = 2 - 1 a e .

We choose a c and a e such that a c a e for any e c and

(3.2) 1 - σ a e < λ e if e is singular ;    a e = 1 if e is regular ;
(3.3) a c < λ c + 1 2 if c is singular ;    a c = 1 if c and all e c are regular .

Let 𝒯 n be the mesh obtained after n anisotropic refinements (Algorithm 3.2) from 𝒯 0 based on the parameters κ c and κ e defined by a e and a c through (3.1)–(3.3). Then the proposed linear finite element approximation u n to equation (1.1) is defined by (2.3) on the mesh 𝒯 n .

Remark 3.5.

For any c 𝒞 , recall κ e c : = min e c ( κ c , κ e ) in Algorithm 3.2. Based on the selections in (3.1)–(3.3), it is clear that, for any c 𝒞 , κ e c = κ c . Note that a e has a lower bound 1 - σ in (3.2). Condition (2.7), σ > λ max e ( 1 - λ e ) , ensures the set given in (3.2) is not empty. For 0 < a e < 1 , it is clear that refinements for an e- or ev-tetrahedron lead to anisotropic meshes toward the edge that do not preserve the maximum angle condition. Namely, the maximum edge angle in the face of the tetrahedron approaches π as the level of refinement n increases. This is a main difficulty that we shall overcome in the error analysis.

3.2 Mesh Layers

To better facilitate the error analysis for the proposed finite element algorithm (Algorithm 3.4) solving equation (1.1), for each initial tetrahedron T ( 0 ) 𝒯 0 , we introduce the mesh layers that are the collections of tetrahedra in 𝒯 n .

We first define mesh layers for a v- or v e -tetrahedron in 𝒯 0 .

Definition 3.6 (Mesh Layers in v- and v e -Tetrahedra).

Let T ( 0 ) = 4 x 0 x 1 x 2 x 3 𝒯 0 be either a v- or a v e -tetrahedron with x 0 𝒞 or x 0 e . We use a local Cartesian coordinate system such that x 0 is the origin. For 1 i n , the ith refinement on T ( 0 ) produces a small tetrahedron with x 0 as a vertex and with one face, denoted by P v , i , parallel to the face 3 x 1 x 2 x 3 of T ( 0 ) . See Figure 1 for example. Then, after n refinements, we define the ith mesh layer L v , i of T ( 0 ) , 1 i < n , as the region in T ( 0 ) between P v , i and P v , i + 1 . We denote by L v , 0 the region in T ( 0 ) between 3 x 1 x 2 x 3 and P v , 1 , and let L v , n be the small tetrahedron with x 0 as a vertex that is bounded by P v , n and three faces of T ( 0 ) . Since x 0 is the only point for the special refinement, we drop the sub-index in the grading parameter. Namely, for such T ( 0 ) , we use κ = 2 - 1 a to denote the grading parameter near x 0 ( κ = κ c if x 0 𝒞 , and κ = κ e if x 0 e ). See the second picture in Figure 2. Then, by Algorithm 3.2, the dilation matrix

(3.4) 𝐁 v , i : = ( κ - i 0 0 0 κ - i 0 0 0 κ - i )

maps L v , i to L v , 0 for 0 i < n , and maps L v , n to T ( 0 ) . We define the initial triangulation of L v , i , 0 i < n , to be the first decomposition of L v , i into tetrahedra. Thus, the initial triangulation of L v , i consists of those tetrahedra in 𝒯 i + 1 that are contained in the layer L v , i .

Figure 2 
                  Anisotropic triangulations after two consecutive refinements and mesh layers on an initial tetrahedron (left to right): o-tetrahedron, v- or 
                        
                           
                              
                                 v
                                 e
                              
                           
                           
                           {v_{e}}
                        
                     -tetrahedron (
                        
                           
                              
                                 κ
                                 =
                                 0.3
                              
                           
                           
                           {\kappa=0.3}
                        
                     ), e-tetrahedron (
                        
                           
                              
                                 
                                    κ
                                    e
                                 
                                 =
                                 0.3
                              
                           
                           
                           {\kappa_{e}=0.3}
                        
                     ); ev-tetrahedron (
                        
                           
                              
                                 
                                    κ
                                    
                                       e
                                       ⁢
                                       c
                                    
                                 
                                 =
                                 
                                    κ
                                    c
                                 
                                 =
                                 0.3
                              
                           
                           
                           {\kappa_{ec}=\kappa_{c}=0.3}
                        
                     , 
                        
                           
                              
                                 
                                    κ
                                    e
                                 
                                 =
                                 0.4
                              
                           
                           
                           {\kappa_{e}=0.4}
                        
                     ).
Figure 2 
                  Anisotropic triangulations after two consecutive refinements and mesh layers on an initial tetrahedron (left to right): o-tetrahedron, v- or 
                        
                           
                              
                                 v
                                 e
                              
                           
                           
                           {v_{e}}
                        
                     -tetrahedron (
                        
                           
                              
                                 κ
                                 =
                                 0.3
                              
                           
                           
                           {\kappa=0.3}
                        
                     ), e-tetrahedron (
                        
                           
                              
                                 
                                    κ
                                    e
                                 
                                 =
                                 0.3
                              
                           
                           
                           {\kappa_{e}=0.3}
                        
                     ); ev-tetrahedron (
                        
                           
                              
                                 
                                    κ
                                    
                                       e
                                       ⁢
                                       c
                                    
                                 
                                 =
                                 
                                    κ
                                    c
                                 
                                 =
                                 0.3
                              
                           
                           
                           {\kappa_{ec}=\kappa_{c}=0.3}
                        
                     , 
                        
                           
                              
                                 
                                    κ
                                    e
                                 
                                 =
                                 0.4
                              
                           
                           
                           {\kappa_{e}=0.4}
                        
                     ).
Figure 2 
                  Anisotropic triangulations after two consecutive refinements and mesh layers on an initial tetrahedron (left to right): o-tetrahedron, v- or 
                        
                           
                              
                                 v
                                 e
                              
                           
                           
                           {v_{e}}
                        
                     -tetrahedron (
                        
                           
                              
                                 κ
                                 =
                                 0.3
                              
                           
                           
                           {\kappa=0.3}
                        
                     ), e-tetrahedron (
                        
                           
                              
                                 
                                    κ
                                    e
                                 
                                 =
                                 0.3
                              
                           
                           
                           {\kappa_{e}=0.3}
                        
                     ); ev-tetrahedron (
                        
                           
                              
                                 
                                    κ
                                    
                                       e
                                       ⁢
                                       c
                                    
                                 
                                 =
                                 
                                    κ
                                    c
                                 
                                 =
                                 0.3
                              
                           
                           
                           {\kappa_{ec}=\kappa_{c}=0.3}
                        
                     , 
                        
                           
                              
                                 
                                    κ
                                    e
                                 
                                 =
                                 0.4
                              
                           
                           
                           {\kappa_{e}=0.4}
                        
                     ).
Figure 2 
                  Anisotropic triangulations after two consecutive refinements and mesh layers on an initial tetrahedron (left to right): o-tetrahedron, v- or 
                        
                           
                              
                                 v
                                 e
                              
                           
                           
                           {v_{e}}
                        
                     -tetrahedron (
                        
                           
                              
                                 κ
                                 =
                                 0.3
                              
                           
                           
                           {\kappa=0.3}
                        
                     ), e-tetrahedron (
                        
                           
                              
                                 
                                    κ
                                    e
                                 
                                 =
                                 0.3
                              
                           
                           
                           {\kappa_{e}=0.3}
                        
                     ); ev-tetrahedron (
                        
                           
                              
                                 
                                    κ
                                    
                                       e
                                       ⁢
                                       c
                                    
                                 
                                 =
                                 
                                    κ
                                    c
                                 
                                 =
                                 0.3
                              
                           
                           
                           {\kappa_{ec}=\kappa_{c}=0.3}
                        
                     , 
                        
                           
                              
                                 
                                    κ
                                    e
                                 
                                 =
                                 0.4
                              
                           
                           
                           {\kappa_{e}=0.4}
                        
                     ).
Figure 2

Anisotropic triangulations after two consecutive refinements and mesh layers on an initial tetrahedron (left to right): o-tetrahedron, v- or v e -tetrahedron ( κ = 0.3 ), e-tetrahedron ( κ e = 0.3 ); ev-tetrahedron ( κ e c = κ c = 0.3 , κ e = 0.4 ).

Now, we define mesh layers for an initial e-tetrahedron T ( 0 ) .

Definition 3.7 (Mesh Layers in e-Tetrahedra).

Based on Algorithm 3.2, in each refinement, an e-tetrahedron is cut by a parallelogram parallel to x 0 x 1 . For example, in the e-tetrahedron of Figure 1, the quadrilateral with vertices x 02 , x 12 , x 13 , x 03 is the aforementioned parallelogram. We denote by P e , i the parallelogram produced in the ith refinement, 1 i n . For the mesh 𝒯 n , let the ith layer L e , i on T ( 0 ) , 0 < i < n , be the region bounded by P e , i , P e , i + 1 and the faces of T ( 0 ) . Define L e , 0 to be the sub-region of T ( 0 ) away from e that is separated by P e , 1 . Define L e , n to be the sub-region of T ( 0 ) between P e , n and e. See also the third picture in Figure 2. As in Definition 3.6, the initial triangulation of the layer L e , i , 0 i < n , consists of the tetrahedra in 𝒯 i + 1 that are contained in L e , i . Therefore, r e κ e i on L e , i , 0 i < n .

In addition, we have the following anisotropic mapping that transforms a tetrahedron in L e , i to a reference element [21, Lemma 4.13].

Proposition 3.8.

Let T ( i + 1 ) T i + 1 be a tetrahedron such that T ( i + 1 ) L e , i T ( 0 ) , 0 i < n . Then T ( i + 1 ) is contained either in an e-tetrahedron in T i or in a v e -tetrahedron in T i .

Case I: T ( i + 1 ) is contained in an e-tetrahedron T ( i ) T i . Using a T ( i ) -based local coordinate system, there is a transformation

(3.5) 𝐁 e , i = ( κ e - i 0 0 0 κ e - i 0 b 1 κ e - i b 2 κ e - i 2 i )

that maps T ( i + 1 ) to one of the four o-tetrahedra in T ^ 1 (hence, we have finitely many reference elements for all T ( i + 1 ) ). Here, T ^ 1 is the triangulation on a reference tetrahedron T ^ that is obtained after one graded refinement to the edge. For an e-tetrahedron in the last layer T ( n ) L e , n T ( 0 ) , using a T ( n ) -based local coordinate system, there exists a transformation B e , n of form (3.5) with i = n that maps T ( n ) to a reference tetrahedron T ^ .

Case II: T ( i + 1 ) is contained in a v e -tetrahedron T ( i ) T i . Let T ( k ) T k , 1 k i , be the v e -tetrahedron such that T ( i ) T ( k ) and T ( k ) is contained in an e-tetrahedron T ( k - 1 ) T k - 1 . Using a T ( k - 1 ) -based local coordinate system, there is a transformation

(3.6) 𝐁 i , k = ( κ e - i + 1 0 0 0 κ e - i + 1 0 b 1 κ e - i + 1 b 2 κ e - i + 1 2 k - 1 κ e k - i )

that maps T ( i + 1 ) to one of the o-tetrahedra in T ^ 2 (as in Case I, we again have finitely many reference elements for all T ( i + 1 ) ). Here, T ^ 2 is the triangulation on a reference tetrahedron T ^ that is obtained after two graded refinements to the edge. For a v e -tetrahedron in the last layer T ( n ) L e , n T ( 0 ) , let T ( k ) T k be the v e -tetrahedron such that T ( n ) T ( k ) and T ( k ) is contained in an e-tetrahedron T ( k - 1 ) T k - 1 . Using a T ( k - 1 ) -based local coordinate system, there exists a transformation B n , k of form (3.6) with i = n that maps T ( n ) to a v e -tetrahedron in T ^ 1 .

In both cases, | b 1 | , | b 2 | C 0 for C 0 > 0 depending on T ( 0 ) but not on i, n or k.

In addition, we define the mesh layers on an initial ev-tetrahedron T ( 0 ) 𝒯 0 .

Definition 3.9 (Mesh Layers in ev-Tetrahedra).

For 1 i n , the ith refinement on T ( 0 ) produces a small tetrahedron with x 0 as a vertex. We denote by P e v , i the face of this small tetrahedron whose closure does not contain x 0 (see the last two pictures in Figure 1). Then, for the mesh 𝒯 n on T ( 0 ) , we define the ith mesh layer L e v , i , 1 i < n , as the region in T ( 0 ) between P e v , i and P e v , i + 1 . We define L e v , 0 to be the region in T ( 0 ) between 3 x 1 x 2 x 3 and P e v , 1 and let L e v , n T ( 0 ) be the small tetrahedron with x 0 as a vertex that is generated in the nth refinement.

Given the condition κ e c = κ c in Algorithm 3.4, we see that the layer L e v , i and the layer L v , i in Definition 3.6 are obtained from the same procedure. Therefore, use a local Cartesian coordinate system such that c is the origin. For 0 i n , the mapping

(3.7) 𝐁 e v , i = ( κ c - i 0 0 0 κ c - i 0 0 0 κ c - i )

is a bijection from L e v , i to L ^ , where L ^ is the reference domain for L e v , i that satisfies L ^ : = T ( 0 ) when i = n and L ^ : = L e v , 0 when 0 i < n . Recall that one graded refinement using the same parameters κ c and κ e gives rise to a triangulation on T ( 0 ) , which we denote by 𝒯 ^ 1 . We further denote by ^ the initial triangulation of L e v , 0 that consists of the seven tetrahedra in 𝒯 ^ 1 .

4 Regularity Results in a Dihedron

In this section, we develop new regularity estimates for equation (1.1), especially in the region that is close to the edges where different boundary conditions are imposed.

Let D = K × be a dihedron, with K a two-dimensional cone of center 0 and opening angle ω. In this domain, we consider u H 1 ( D ) with a support included in ( K B ( 0 , R ) ) × for some R > 0 to be the solution of

(4.1) { - Δ u = f in D , u = 0 on Γ Dir × , n u = 0 on Γ Neu × ,

where f H σ ( D ) for some σ [ 0 , 1 ) and Γ Dir Γ Neu = K such that Γ Dir (resp. Γ Neu ) is either empty or a full half-line. In that way, we consider either the pure Dirichlet, the pure Neumann or the mixed problem. In this section, to simplify the exposition, we use ( x 1 , x 2 , x 3 ) (instead of ( x , y , z ) ) to denote a point in D and suppose the edge of D is on the x 3 -axis.

The behavior of this solution is well known in the case of the pure Dirichlet problem [17, 13] for data in L 2 but is less studied for smoother data and in the two other cases of boundary conditions. Our goal is to show that this solution is decomposed into a regular part and a singular one with the appropriate behavior. For that purpose, we perform a partial Fourier transform in the x 3 -variable that allows to reduce the study to a Helmholtz equation in K.

4.1 Helmholtz Equation in a Cone

For all ξ , we consider the solution v H 1 ( K ) with a support included in K B ( 0 , R ) of

(4.2) { - Δ v + ξ 2 v = g in K , v = 0 on Γ Dir , n v = 0 on Γ Neu ,

where g H σ ( K ) for some σ [ 0 , 1 ) and show that v admits a decomposition into a regular part and a singular one. Recall that the singularities of problem (4.2) are related to the singularities of the Laplace equation, namely to the singularities of problem (4.2) with ξ = 0 . Such singularities are in the form [16, 13]

S k = r λ k φ k ( θ )

with

(4.3) λ k = k π ω for all k * = { 0 }

in the pure Dirichlet and Neumann case, while

(4.4) λ k = ( 2 k - 1 ) π 2 ω for all k *

in the mixed case. Here, r is the distance to the vertex of K. For shortness, the smallest singular exponent λ 1 is denoted by λ. The function φ k is given by φ k ( θ ) = sin ( λ k θ ) in the pure Dirichlet case and in the mixed case, while φ k ( θ ) = cos ( λ k θ ) in the pure Neumann case.

Now, we can prove the next result.

Theorem 4.1.

Let σ [ 0 , 1 ) be such that σ λ k - 1 for all k N * . Then, for all ξ R , the solution v H 1 ( K ) with a support included in K B ( 0 , R ) of (4.2) with g H σ ( D ) can be split up into

(4.5) v = v reg ( ξ ) + v sing ( ξ ) ,

with v reg ( ξ ) H 2 + σ ( K ) and v sing ( ξ ) V δ 2 ( K ) H 1 ( K ) for any δ > 1 - λ satisfying the estimates[1]

(4.6) v reg ( ξ ) 2 + σ , K , 1 + | ξ | g σ , K , 1 + | ξ | ,
(4.7) ( 1 + ξ 2 ) r - σ v reg ( ξ ) K g σ , K , 1 + | ξ | ,
(4.8) ( 1 + | ξ | δ + σ ) v sing ( ξ ) V δ 2 ( K ) g σ , K , 1 + | ξ | ,
(4.9) ( 1 + | ξ | ) | v sing ( ξ ) | 1 , K g σ , K , 1 + | ξ | .

Proof.

We distinguish the case | ξ | > 1 to the case | ξ | 1 .

(a) For | ξ | > 1 , as in the proof of [13, Theorem 16.9], we use a scaling argument; namely, by setting x ^ = | ξ | x and v ^ ( x ^ ) = v ( x ) , we see that v ^ is solution of

(4.10) { - Δ v ^ + v ^ = g ^ in K , v ^ = 0 on Γ Dir , n v ^ = 0 on Γ Neu ,

where g ^ H σ ( K ) is defined by g ^ ( x ^ ) = ξ - 2 g ( x ) . Clearly, the weak formulation of this problem is

( v ^ , w ) 1 , K = K g ^ w d x for all w H Γ Dir 1 ( K ) ,

where H Γ Dir 1 ( K ) : = { w H 1 ( K ) : w = 0 on Γ Dir } is a Hilbert space with its natural inner product

( u , w ) 1 , K = K ( u w + u w ) d x for all u , w H Γ Dir 1 ( K ) .

As a direct consequence, the above problem has a unique solution v ^ H Γ Dir 1 ( K ) with the continuous dependence

(4.11) v ^ 1 , K g ^ 0 , K .

Now, using a localization argument, by [13, Theorem 23.7] (see also [13, Remark 23.8]) near the origin and the standard shift theorem far from the origin, one deduces that v ^ admits the splitting

v ^ = v ^ reg + η ( r ^ ) k * : 0 < λ k < 1 + σ c k r ^ λ k φ k ,

where η 𝒟 ( 2 ) is a smooth cut-off function equal to 1 in a neighborhood of the origin that, without loss of generality, is assumed to have a support included in K B ( 0 , R ) , v ^ reg H 2 + σ ( K ) and c k with

v ^ reg 2 + σ , K + k * : 0 < λ k < 1 + σ | c k | g ^ σ , K + v ^ 1 , K .

Combined with (4.11), we find that

(4.12) v ^ reg 2 + σ , K + k * : 0 < λ k < 1 + σ | c k | g ^ σ , K .

By a transformation back, this yields (4.5) by setting v reg ( x ) = v ^ reg ( x ^ ) and[2]

v sing ( ξ ) = η ( | ξ | r ) k * : 0 < λ k < 1 + σ c k | ξ | λ k r λ k φ k .

Furthermore, using [13, Lemma AA.19], estimate (4.12) is equivalent to

(4.13) v reg 2 + σ , K , | ξ | + | ξ | 1 + σ k * : 0 < λ k < 1 + σ | c k | g σ , K , | ξ | .

This estimate directly leads to the first estimate (4.6) recalling that | ξ | > 1 . To prove the second one, we first notice that, the support of v sing being included in B ( 0 , R | ξ | ) B ( 0 , R ) , v reg has a compact support included in K B ( 0 , R ) . Furthermore, using the estimate v reg 2 + σ , K B ( 0 , R ) , | ξ | g σ , K , | ξ | and an interpolation inequality (see [16, Theorem 1.4.3.3]), we find that v reg σ , K B ( 0 , R ) | ξ | - 2 g σ , K , | ξ | . Since [13, Theorem AA.7] guarantees that H σ ( K B ( 0 , R ) ) = V 0 σ ( K B ( 0 , R ) ) , we deduce that

r - σ v reg 0 , K = r - σ v reg 0 , K B ( 0 , R ) v reg σ , K B ( 0 , R ) | ξ | - 2 g σ , K , | ξ | ,

which is exactly (4.7).

To prove (4.8), it suffices to check that, for all k * such that 0 < λ k < 1 + σ , one has

| ξ | δ + σ + λ k | c k | η ( | ξ | r ) r λ k φ k V δ 2 ( K ) g σ , K , | ξ | ,

which, in view of (4.13), holds as soon as

(4.14) | ξ | δ + λ k - 1 η ( | ξ | r ) r λ k φ k V δ 2 ( K ) 1 .

Now, using polar coordinates, one can show that

η ( | ξ | r ) r λ k φ k V δ 2 ( K ) 2 0 r 2 ( δ + λ k - 2 ) { | η ( | ξ | r ) | 2 + | | ξ | r η ( | ξ | r ) | 2 + | | ξ | 2 r 2 η ′′ ( | ξ | r ) | 2 } r d r ,

and by the change of variable r ^ = | ξ | r , one finds

η ( | ξ | r ) r λ k φ k V δ 2 ( K ) 2 | ξ | 2 ( δ + λ k - 1 ) 0 r ^ 2 ( δ + λ k - 2 ) { | η ( r ^ ) | 2 + | r ^ η ( r ^ ) | 2 + | r ^ 2 η ′′ ( r ^ ) | 2 } r ^ d r ^ .

The integral of this right-hand side being finite as soon as δ + λ k - 1 > 0 (which holds if δ + λ - 1 > 0 ), we have found that (4.14) is valid. The proof of (4.8) is fully similar and is left to the reader.

(b) For | ξ | 1 , we first notice that

(4.15) v 1 , K | v | 1 , K

because v has a compact support included into K B ( 0 , R ) . Since the weak formulation of problem (4.2) is

K ( v w + ξ 2 v w ) d x = K g w d x for all w H Γ Dir 1 ( K ) ,

by taking w = v in this identity and using (4.15), we find

v 1 , K 2 | v | 1 , K 2 K ( | v | 2 + ξ 2 | v | 2 ) d x = K g v d x .

Consequently, by Cauchy–Schwarz’s inequality, we get

(4.16) v 1 , K g 0 , K .

Now, v can be seen as the solution of (compare with (4.10))

{ - Δ v + v = g ~ in K , v = 0 on Γ Dir , n v = 0 on Γ Neu ,

where g ~ : = g + v - ξ 2 v H σ ( K ) that, owing to (4.16), satisfies (recalling that | ξ | 1 ) g ~ σ , K g σ , K . As in the previous point, we then get the decomposition

v = v reg + η ( r ) k * : 0 < λ k < 1 + σ c k r λ k φ k ,

where v reg H 2 + σ ( K ) and c k with

(4.17) v reg 2 + σ , K + k * : 0 < λ k < 1 + σ | c k | g ~ σ , K + v 1 , K g ~ σ , K .

This yields (4.5) with[3]

v sing ( ξ ) = η ( r ) k * : 0 < λ k < 1 + σ c k r λ k φ k ,

estimate (4.17) corresponding to (4.13) with | ξ | = 1 . Estimate (4.6) is a direct consequence of (4.17), while estimates (4.7)–(4.9) follow by using the previous arguments simply replacing | ξ | by 1. ∎

4.2 Singular Decomposition in a Dihedron

Define the weighted space γ , σ 2 ( D ) in the dihedron D = K × ,

(4.18) γ , σ 2 ( D ) : = { v H loc 2 ( D ) r - 1 - γ v , r - γ v , r - 1 3 v L 2 ( D ) , r 1 - γ 2 v , 3 v , r - σ 3 2 v L 2 ( D ) }

with the norm

v γ , σ 2 ( D ) 2 : = r - σ 3 2 v L 2 ( D ) 2 + | α | = 1 α 3 v L 2 ( D ) 2 + r - 1 3 v L 2 ( D ) 2 + | α | 2 r | α | - 1 - γ α v L 2 ( D ) 2 ,

where means the first-order derivatives in the x 1 , x 2 variables, 3 = x 3 , while α means that the third component of the multi-index is zero.

Then we have the following regularity estimates for equation (4.1).

Theorem 4.2.

Let σ [ 0 , 1 ) be such that σ λ k - 1 for all k N * . Recall λ : = λ 1 . Suppose f H σ ( D ) . Then the solution u H 1 ( D ) of (4.1) with a support included in ( K B ( 0 , R ) ) × R for some R > 0 can be split up into

(4.19) u = u reg + u sing ,

with u reg H 2 + σ ( D ) such that r - σ 3 2 u reg L 2 ( D ) and u sing H γ , σ 2 ( D ) for any γ < λ satisfying the estimates

(4.20) u reg 2 + σ , D f σ , D ,
(4.21) r - σ 3 j u reg 0 , D f σ , D for all j = 0 , 1 , 2 ,
(4.22) u sing γ , σ 2 ( D ) f σ , D .

Proof.

We perform a partial Fourier transform in x 3 . Namely, let v ( ξ ) = x 3 ξ u and g ( ξ ) = x 3 ξ f . Then we see that v is solution of (4.2). Applying Theorem 4.1 to v and performing inverse Fourier transform, we find decomposition (4.19) with u reg = x 3 ξ - 1 v reg , u sing = x 3 ξ - 1 v sing . Estimate (4.20) (resp. (4.21)) follows from (4.6) (resp. (4.7)) and [13, Proposition AA.20]. Similarly, using estimate (4.8), we get (since δ + σ > 0 )

| α | 2 r | α | - 2 + δ α u sing L 2 ( D ) 2 f σ , D 2 for all δ > 1 - λ .

By setting γ = 1 - δ , this yields

(4.23) | α | 2 r | α | - 1 - γ α u sing L 2 ( D ) 2 f σ , D 2 .

Again, applying (4.8) with δ = 2 - σ (resp. δ = 0 ) that is clearly larger than 1 - λ , we find

(4.24) r - σ 3 2 u sing L 2 ( D ) 2 f σ , D 2 ,
(4.25) r - 1 3 u sing L 2 ( D ) 2 f σ , D 2 .

Finally, applying (4.9), we clearly obtain

(4.26) 3 u sing L 2 ( D ) 2 f σ , D 2 .

Estimates (4.23)–(4.26) show that (4.22) holds. ∎

5 Regularity Results in a Dihedral Cone

In this section, we investigate the regularity of the solution of (1.1) in the region where the vertex and the edges meet. Let Γ be a dihedral cone of 3 of vertex c 𝒞 (that can be identified with 0), in the sense that

Γ = { x 3 : x | x | G } ,

with G an open subset of the unit sphere S 2 with a piecewise smooth boundary, each smooth part being included in a great circle.

Let γ c and γ e be the parameters corresponding to a vertex c 𝒞 and an edge e of Γ, respectively. Let 𝜸 be the collection of all the parameters γ c and γ e for Γ. Recall that R c and r e are the distances to the vertex c and to the edge e, respectively. Then, given γ c , γ e 0 for all edges e of Γ, we define the weighted space

(5.1) 𝜸 , σ 2 ( Γ ) : = { v H loc 2 ( Γ ) R c γ e - γ c r e - 1 - γ e v , R c γ e - γ c r e - γ e v , R c 1 - γ c r e - 1 3 v L 2 ( 𝒱 c e ) , R c γ e - γ c r e 1 - γ e 2 v , R c 1 - γ c 3 v , R c 1 + σ - γ c r e - σ 3 2 v L 2 ( 𝒱 c e ) ; r e - 1 - γ e v , r e - γ e v , r e - 1 3 v L 2 ( 𝒱 e 0 ) , r e 1 - γ e 2 v , 3 v , r e - σ 3 2 v L 2 ( 𝒱 e 0 ) ; R c - 1 - γ c v , R c - γ c v , R c - γ c 3 v L 2 ( 𝒱 c 0 ) , R c 1 - γ c 2 v , R c 1 - γ c 3 v , R c 1 - γ c 3 2 v L 2 ( 𝒱 c 0 ) } ,

with the norm

v 𝜸 , σ 2 ( Γ ) 2 : = v H 2 ( 𝒱 0 ) 2 + e c ( R c 1 - γ c θ c , e - σ 3 2 v L 2 ( 𝒱 c e ) 2 + | α | = 1 R c 1 - γ c α 3 v L 2 ( 𝒱 c e ) 2 + R c - γ c θ c , e - 1 3 v L 2 ( 𝒱 c e ) 2 + | α | 2 R c | α | - 1 - γ c θ c , e | α | - 1 - γ e α v L 2 ( 𝒱 c e ) 2 ) + c 𝒞 , | α | 2 R c | α | - 1 - γ c α v L 2 ( 𝒱 c 0 ) 2 + e ( r e - σ 3 2 v L 2 ( 𝒱 e 0 ) 2 + | α | = 1 α 3 v L 2 ( 𝒱 e 0 ) 2 + r e - 1 3 v L 2 ( 𝒱 e 0 ) 2 + | α | 2 r e | α | - 1 - γ e α v L 2 ( 𝒱 e 0 ) 2 ) ,

where 3 is the derivative in the direction of e, α = 1 α 1 2 α 2 for α = ( α 1 , α 2 ) , and α = ( α 1 , α 2 , α 3 ) .

In this domain, we consider u H 1 ( Γ ) with a support included in Γ B ( 0 , R ) for some R > 0 being the solution of

(5.2) { - Δ u = f in Γ , u = 0 on Γ Dir , n u = 0 on Γ Neu ,

where f H σ ( C ) for some σ [ 0 , 1 ) and Γ Dir Γ Neu = Γ such that Γ Dir (resp. Γ Neu ) is either empty or a finite union of plane faces. Denote by γ Dir = Γ Dir S 2 .

Since u is only in H 1 , by a solution of (5.2) we mean that u H Dir 1 ( Γ ) = { v H 1 ( Γ ) u = 0 on Γ Dir } satisfies

(5.3) Γ u v d x = Γ f v d x for all v H Dir 1 ( Γ ) .

Note that the vertex singular exponent of problem (5.3) near c (see [16, 3]) is given by - 1 2 ± ν c , k + 1 4 , where { ν c , k } k = 0 is the spectrum (enumerated in increasing order and repeated according to its multiplicity) of the non-negative Laplace–Beltrami operator L G mixed on the intersection G between Γ and the unit sphere with Dirichlet boundary condition on Γ Dir G and Neumann boundary condition on Γ Neu G . Here we are only interested in exponents larger than - 1 2 ; hence, for all k : = { 0 , 1 , 2 , } , we set λ c , k = - 1 2 + ν c , k + 1 4 , which is always non-negative. The associated singular function σ c , k is given by

(5.4) σ c , k = R c λ c , k φ c , k ,

where φ c , k is the eigenvector of L G mixed associated with ν c , k , namely L G mixed φ c , k = ν c , k φ c , k . Note that λ c , 0 = 0 if and only if Γ Dir G is empty, and in that case, σ c , 0 = 1 ; otherwise, λ c , 0 > 0 .

Note that, in Γ, u may consist of singularities from the edge and singularities from the vertex. In a first step, we subtract from u its corner (vertex) singularities. Namely, as in [13, Lemma 17.4], we have the following result.

Lemma 5.1.

Let σ [ 0 , 1 ) be such that λ c , k σ + 1 2 for all k N . Let u H 1 ( Γ ) with a support included in Γ B ( 0 , R ) for some R > 0 be the solution of equation (5.3) with f H σ ( Γ ) . Then u admits the splitting

(5.5) u = u 0 + - 1 2 < λ c , k < σ + 1 2 c k σ c , k ,

where u 0 V - ( 1 + σ ) 1 ( Γ ) , c k C and Δ u 0 H σ ( Γ ) = V 0 σ ( Γ ) .

Proof.

By [13, Proposition AA.27], the Mellin transform [ u ] ( λ ) of u exists for all λ with λ = - 1 2 and is the variational solution of ( Δ + λ ( λ + 1 ) ) [ u ] ( λ ) = [ f ] ( λ - 2 ) . Since, by assumption on the line λ = σ + 1 2 , the operator Δ + λ ( λ + 1 ) is invertible from H Dir 1 ( G ) = { v H 1 ( G ) v = 0 on γ Dir } into its dual with

( Δ + λ ( λ + 1 ) ) - 1 h 1 , G , 1 + | λ | h G ,

by the inverse Mellin transform on the line λ = σ + 1 2 , we find the result (5.5) as in [13, Lemma 17.4]. ∎

We now split up u 0 into a regular part and a singular one that contains the edge contribution.

Lemma 5.2.

Under the assumption of Lemma 5.1, u 0 admits the splitting u 0 = u reg + u sing with u reg V 0 2 + σ ( Γ ) , r e - σ 3 2 u reg L 2 ( Γ ) and u sing H 𝛄 , σ 2 ( Γ ) with γ c = 1 + σ and γ e < λ e , where λ e is the smallest exponent λ 1 determined in either (4.3) or (4.4) according to the associated boundary condition with ω e instead of ω.

Proof.

We start as in the proof of [13, Proposition 17.12] by setting

w ( t , θ ) = e η t u 0 ( e t , θ ) , h ( t , θ ) = e ( η + 2 ) t ( Δ u 0 ) ( e t , θ )

with η = - ( σ + 1 2 ) . These functions have the regularity w H 1 ( × G ) , h H σ ( × G ) , and w is the weak solution of (meaning that w satisfies the Dirichlet condition and Neumann condition in a weak sense)

( Δ + t 2 + ( 1 - 2 η ) t + η ( η - 1 ) ) w = h .

Since H 1 ( × G ) is embedded into H σ ( × G ) , this implies that

( Δ + t 2 + ( 1 - 2 η ) t ) w = h ~ = h - η ( η - 1 ) w H σ ( × G ) .

Then, as in Section 4, we apply a partial Fourier transform in t to find that W = t ξ v is the weak solution of ( Δ - ξ 2 + ( 1 - 2 η ) i ξ ) ) W = H with H = t ξ h ~ . This operator is mainly the same one as in problem (4.2), and therefore we conclude that W admits a splitting similar to (4.5). Hence, taking the Fourier transform back, we find that (see Theorem 4.2) w = w reg + w sing with w reg H 2 + σ ( × G ) such that θ c , e - σ t j w reg L 2 ( × G ) for j = 0 , 1 and  2 (recalling that θ c , e is the distance to the edges in × G and hence the angular distance in Γ), and w sing γ , σ 2 ( × G ) for any γ e < λ e for all e.

Coming back to u 0 , we find the result by setting (recalling that R c is the distance to the vertex c)

u reg ( R c , θ ) = R c - η w reg ( R c , θ ) , u sing ( r , θ ) = R c - η w sing ( R c , θ ) .

Indeed, the regularity u reg V 0 2 + σ ( Γ ) follows from w reg H 2 + σ ( × G ) by using [13, Theorem AA.3], while the property r e - σ 3 2 u reg L 2 ( Γ ) follows from the expression of 3 2 in spherical coordinates, an Euler’s change of variables and the regularities of w reg mentioned above (noticing that r e - σ R c - σ θ c , e - σ ).

The regularity of u sing is proved similarly. ∎

In summary, we have the following decomposition of the solution u of equation (5.3).

Corollary 5.3.

Under the assumption of Lemma 5.1, u H 1 ( Γ ) with a support included in Γ B ( 0 , R ) for some R > 0 being the solution of equation (5.3) with f H σ ( Γ ) for some σ [ 0 , 1 ) admits the splitting

u = u reg + u sing + - 1 2 < λ c , k < σ + 1 2 c k ψ σ c , k

with u reg V 0 2 + σ ( Γ ) , r e - σ 3 2 u reg L 2 ( Γ ) and u sing H 𝛄 , σ 2 ( Γ ) with γ c = 1 + σ and γ e < λ e , and ψ being a smooth (and radial) cut-off function with a compact support and equal to 1 on the support of u.

Proof.

Since u = ψ u , the result follows from the two previous lemmas and Leibniz’s rule. ∎

6 Regularity Analysis

In this section, we obtain anisotropic regularity results for equation (1.1) with f H σ ( Ω ) for some σ [ 0 , 1 ) . Our analysis is based on regularity estimates in different sub-regions (see (2.5)) of the domain near the vertices and edges and uses a localization argument.

6.1 Regularity Estimates in 𝒱 e 0

Let us start with an improved regularity of the solution along the edges. For that purpose, for any e and any point ξ e e , we can fix a Cartesian system of coordinates x e = ( x e , 1 , x e , 2 , x e , 3 ) such that ξ e corresponds to ( 0 , 0 , 0 ) . In such a situation, we can fix a cut-off function η ξ e in the form η ξ e ( x e ) = η 0 ( x e , 1 , x e , 2 ) η 1 ( x e , 3 ) with η 0 , η 1 two cut-off functions such that η 0 (resp. η 1 ) is equal to 1 near ( 0 , 0 ) (resp. 0) with a sufficiently small support such that the support of η ξ e is included in 𝒱 e 0 . For all k * , we further denote by λ e , k the associated singular exponents defined by (4.3) or (4.4) (according to the boundary conditions on its adjacent faces) with ω e instead of ω.

Theorem 6.1.

Recall the space H γ , σ 2 from (4.18) and λ e from Lemma 5.2. Let u H Γ Dir 1 ( Ω ) be the solution of (2.2) with f H σ ( Ω ) for some σ [ 0 , 1 ) such that σ λ e , k - 1 for all k N * . Then, for any e E and any point ξ e e , we have

(6.1) η ξ e u = u ξ e , reg + u ξ e , sing

with u ξ e , reg H 2 + σ ( V e 0 ) such that r e - σ 3 2 u ξ e , reg L 2 ( V e 0 ) and v ξ e , sing H γ e , σ 2 ( V e 0 ) for any γ e [ 0 , 1 ] and γ e < λ e satisfying the estimates

u ξ e , reg 2 + σ , 𝒱 e 0 f σ , Ω ,
r e - σ 3 j u ξ e , reg 𝒱 e 0 f σ , Ω for all j = 0 , 1 , 2 ,
(6.2) u ξ e , sing γ e , σ 2 ( 𝒱 e 0 ) f σ , Ω .

Proof.

For the sake of simplicity, we drop the indices e and ξ e . Denote by D = K × the dihedral cone that coincides with Ω near ξ, where K is a two-dimensional cone of opening angle ω.

Now, u ~ : = η u is clearly a weak solution of

(6.3) - Δ u ~ = f ~ in D ,

where f ~ is given by f ~ = η f - 2 u η - u Δ η L 2 ( D ) . But the main point is that f ~ actually belongs to H σ ( D ) . Indeed, the first term has trivially this property; the third term is even in H 1 ( D ) , hence also in H σ ( D ) . Finally, for the second term, we will show that it belongs to H 1 ( D ) as well. Indeed, we notice that

u η = η 1 i = 1 , 2 i η 0 i u + η 0 3 η 1 3 u .

Now, the first term belongs to H 1 ( D ) because η 1 i η 0 is zero in a neighborhood of the edges and corners; hence the H 2 regularity of u inside the domain suffices to get the request regularity. On the other hand, for the second term, we notice that the method of tangential differential quotients of Nirenberg (see for instance [16, pp. 87–90]) can be applied to ψ u in the x 3 direction, with a cut-off function similar to η such that ψ 1 on the support of η, and deduce that 3 ( ψ u ) H 1 ( Ω ) . This obviously leads to η 0 3 η 1 3 u H 1 ( D ) .

Once f ~ belongs to H σ ( D ) , we conclude by applying Theorem 4.2. ∎

6.2 Regularity Estimates in 𝒱 c

Now, we describe the extra regularity in a neighborhood of a vertex c 𝒞 . For that purpose, we fix a cut-off function χ c that depends only on the R c variable and such that χ c 1 in a neighborhood of c and with a support included in 𝒱 c (hence χ c 0 in a neighborhood of the other vertices of Ω). We further denote by Γ c the infinite cone that coincides with Ω near c.

Then we have the following splitting of the solution near the vertex.

Theorem 6.2.

Recall the space H 𝛄 , σ 2 from (5.1) and λ e from Lemma 5.2. Under the assumption of Lemma 5.1, let u H Γ Dir 1 ( Ω ) be the solution of (2.2) with f H σ ( Ω ) for some σ [ 0 , 1 ) . Then, for any c C , χ c u admits the splitting

(6.4) χ c u = u c , reg + u c , sing + - 1 2 < λ c , k < σ + 1 2 c k ψ σ c , k ,

with u c , reg V 0 2 + σ ( Γ c ) , r e - σ 3 2 u c , reg L 2 ( V c e ) for all e E having c as an endpoint, and u c , sing H 𝛄 , σ 2 ( Γ c ) with γ c = 1 + σ and γ e [ 0 , 1 ] and γ e < λ e , and finally ψ being a smooth (and radial) cut-off function with a compact support and equal to 1 on the support of χ c .

Proof.

We can apply Corollary 5.3 to χ u (for shortness, we drop the index c) if we show that Δ ( χ u ) H σ ( K ) . As before, using Leibniz’s rule, u ~ : = χ u is clearly a weak solution of (6.3) with f ~ = χ f - 2 u χ - u Δ χ , which actually belongs to H σ ( K ) . The first and third term have trivially this property. Hence, only the second term requires a careful inspection. Due to the choice of χ and using Cartesian coordinates centered at c, we have

u χ = χ ( R ) R i = 1 , 2 , 3 x i i u .

First we may notice that χ is zero near c; hence the regularity of u χ is related to the regularity of u far from the corners. So, as u belongs to H 2 in 𝒱 0 c 𝒞 𝒱 c , we get that

(6.5) u χ H 1 ( 𝒱 0 ) .

Now, for a fixed edge e having c as an endpoint, we can use Cartesian coordinates such that the x 3 -axis contains the edge e and can use splitting (6.1) of Theorem 6.1. The regular part contributes to an H 1 function, so let us show that this is the same for the singular part u e , sing . Indeed, we have to show that

x i i u e , sing H 1 ( 𝒱 c ) for all i = 1 , 2 , 3 .

For i = 3 , this is a direct consequence of (6.2), while for i = 1 or  2 , this is a consequence of (6.2) and of the bound | x i | θ c , e θ c , e 1 - γ e . All together, (6.5) is valid, and the proof is complete. ∎

The third term of splitting (6.4) of χ c u is not in 𝜸 , σ 2 ( Γ B ( 0 , R ) ) because (see (5.4)) σ c , k = R c λ c , k φ c , k and φ c , k is not necessarily equal to zero at the corners of G (intersection of Γ with the unit sphere), but it can be split up into a regular part and a singular one in the spirit of Theorem 6.1 with ξ = 0 . This allows us to show the next result.

Lemma 6.3.

Let λ c , k > 0 be fixed such that λ c , k < σ + 1 2 . Recall λ e from Lemma 5.2. Then v c , k : = ψ σ c , k with σ c , k given by (5.4) can be split up into v c , k = v 1 + v 2 such that v 1 satisfies

R c σ c θ c , e - σ c , e 3 2 v 1 L 2 ( 𝒱 c e ) ,
R c σ c α 3 v 1 L 2 ( 𝒱 c e ) for all | α | = 1 ,
R c σ c α v 1 L 2 ( 𝒱 c e ) for all | α | = 2 ,
R c σ c - 2 + | α | α v 1 L 2 ( 𝒱 e c ) for all | α | 1 ,

and v 2 satisfies

R c σ c θ c , e - σ c , e 3 2 v 2 L 2 ( 𝒱 c e ) ,
R c σ c α 3 v 2 L 2 ( 𝒱 c e ) for all | α | = 1 ,
R c σ c θ c , e σ ~ c , e α v 2 L 2 ( 𝒱 c e ) for all | α | = 2 ,
R c σ c - 1 θ c , e σ ~ c , e - 1 α v 2 L 2 ( 𝒱 c e ) for all | α | = 1 ,
R c σ c - 1 θ c , e - 1 3 v 2 L 2 ( 𝒱 c e ) ,
R c σ c - 2 θ c , e σ ~ c , e - 2 v 2 L 2 ( 𝒱 c e ) ,

for any σ c , e < 1 , any σ ~ c , e > 1 - λ e and any σ c > 1 2 - λ c , where λ c is the smallest positive λ c , k such that - 1 2 < λ c , k < σ + 1 2 and θ c , e ( x ) : = r e ( x ) R c ( x ) is the angular distance.

Proof.

The proof is based on simple calculations expressing the Cartesian derivatives in spherical coordinates ( R c , θ c , e , φ c , e ) (where θ c , e is the angular distance to the edge e). We use the splitting

φ c , k = κ c , e ( 0 ) χ c , e ( θ c , e ) + κ c , e θ c , e λ c , e φ c , k , e ( φ c , e ) + φ c , k , R ,

where κ c , e ( 0 ) , κ c , e , χ c , e is a cut-off function equal to 1 near θ c , e = 0 , φ c , k , e is the singular function associated with the corner singular exponent λ c , e and φ c , k , R is the regular part of φ c , k that either belongs to H 2 ( G ) or has to be split up into a singular part (similar to the second term of the above right-hands side) and a regular part in H 2 ( G ) . In this situation,

v 2 = κ c , e ψ R c λ c , k θ c , e λ c , e φ c , k , e ( φ c , e ) ,

while v 1 = v c , k - v 2 . ∎

Remark 6.4.

The function v c , k = ψ σ c , k = ψ R λ c , k φ c , k that characterizes the vertex singularity in (6.4) satisfies the following. (1) In the case λ c , k = 0 , σ c , k = constant because the eigenvector φ c , k of the Laplace–Beltrami operator corresponding to the zero eigenvalue is a constant. (2) In the case λ c , k > 0 , according to Lemma 6.3, the function v c , k admits a splitting into two functions v 1 and v 2 . In particular, let

(6.6) σ c , e ( 0 , 1 ) , σ c = 1 - a c , σ ~ c , e = 1 - a e ,

where a c and a e are parameters defined for the anisotropic mesh (Algorithm 3.4). Given the conditions in (3.2) and (3.3), we conclude that the selections in (6.6) satisfy the conditions in Lemma 6.3. Namely, σ c > 1 2 - λ c and σ ~ c , e > 1 - λ e . Recall the weighted space 𝜸 2 from (5.1). It is clear that the function v 2 𝜸 , σ 2 ( 𝒱 c e ) with γ e = a e , γ c = a c and σ = σ c , e .

Recall that 𝒱 c 0 is part of the neighborhood of c 𝒞 that excludes the edges. Based on Theorem 6.2 and Lemma 6.3, we further obtain the regularity estimate in 𝒱 c 0 .

Corollary 6.5.

Under the assumption of Lemma 5.1, let u H Γ Dir 1 ( Ω ) be the solution of (2.2) with f H σ ( Ω ) for some σ [ 0 , 1 ) . Let c C be a vertex, and let λ c be the smallest positive λ c , k . If λ c < σ + 1 2 , in V c 0 , u admits the decomposition u = u c , reg + u c , sing , where u c , reg H 2 ( V c 0 ) and u c , sing V σ c 2 ( V c 0 ) for any σ c > 1 2 - λ c . If λ c σ + 1 2 , we have u H 2 ( V c 0 ) .

Proof.

Note that the angular distance θ c , e is bounded below from 0 in 𝒱 c 0 . Taking this into account in (6.4) and in the regularity estimates (Theorem 6.2 and Lemma 6.3), we can derive the result in this corollary by straightforward calculations. ∎

7 Interpolation Error Analysis

In this section, we carry out the interpolation error analysis for the proposed anisotropic finite element scheme (Algorithm 3.4) for equation (1.1) with f H σ ( Ω ) given in (2.7). Different from the error analysis for the pure Dirichlet problem [22], the numerical analysis for the problem with the mixed boundary condition, especially in the case that Neumann boundary conditions are imposed in the adjacent faces of non-smooth points, requires new analytical tools and more involved, since the underlying solution has quite different singular behaviors near vertices and edges as elaborated in previous sections. We shall conduct the analysis on initial tetrahedra according to their types (Definition 3.1).

We first have the estimate for an o-tetrahedron in the initial mesh.

Lemma 7.1.

Let T ( 0 ) T 0 be an o-tetrahedron. For u H 2 ( T ( 0 ) ) , let I u S n be its nodal interpolation on T n . Then we have | u - I u | H 1 ( T ( 0 ) ) C h u H 2 ( T ( 0 ) ) , where h = 2 - n and C is independent of n and u.

Proof.

Based on Algorithms 3.2 and 3.4, the restriction of 𝒯 n on T ( 0 ) is a quasi-uniform mesh with size O ( 2 - n ) . By the standard interpolation error estimate, we obtain | u - I u | H 1 ( T ( 0 ) ) C 2 - n u H 2 ( T ( 0 ) ) C h u H 2 ( T ( 0 ) ) . ∎

7.1 Analysis on Initial v- and v e -Tetrahedra

For an initial v- or v e -tetrahedron T ( 0 ) = 4 x 0 x 1 x 2 x 3 in 𝒯 0 , recall the mesh layers L v , i , 0 i n in Definition 3.6. Based on the refinement, on each L v , i , the tetrahedra in 𝒯 n are isotropic with mesh size O ( κ i 2 i - n ) . In T ( 0 ) , let ρ be the distance to x 0 . Therefore,

(7.1) ρ κ i on L v , i , 0 i < n .

Namely, if T ( 0 ) is a v-tetrahedron, ρ R c for c = x 0 𝒞 , and if T ( 0 ) is a v e -tetrahedron, ρ r e , where e is the edge containing x 0 .

Recall from Theorem 6.1 that the solution admits the following decomposition on an initial v e -tetrahedron T ( 0 ) :

(7.2) u = u ξ e , reg + u ξ e , sing ,

where u ξ e , reg H 2 ( T ( 0 ) ) and u ξ e , sing γ e , σ ( T ( 0 ) ) for γ e [ 0 , 1 ] and γ e < λ e . Meanwhile, by Corollary 6.5, the solution admits the following decompositions on an initial v-tetrahedron T ( 0 ) :

(7.3) u = u c , reg + u c , sing ,

where u c , reg H 2 ( T ( 0 ) ) and u c , sing V σ c 2 ( T ( 0 ) ) for any σ c > 1 2 - λ c , where λ c < σ + 1 2 ; in the case λ c σ + 1 2 , we have u = u c , reg .

Then we have the interpolation error estimate in the layer L v , i .

Lemma 7.2.

For a continuous function v, let Iv be its nodal interpolation on T n . Let T ( 0 ) T 0 be either a v- or a v e -tetrahedron. For v H 2 ( T ( 0 ) ) , we have

| v - I v | H 1 ( L v , i ) C h v H 2 ( L v , i ) , 0 i n .

If T ( 0 ) T 0 is a v-tetrahedron and v V σ c 2 ( T ( 0 ) ) , where σ c satisfies the condition in Corollary 6.5, then

| v - I v | H 1 ( L v , i ) C h v V 1 - a c 2 ( L v , i ) , 0 i < n .

If T ( 0 ) T 0 is a v e -tetrahedron and v H γ e , σ 2 ( T ( 0 ) ) , where γ e satisfies the condition in Theorem 6.1, then

| v - I v | H 1 ( L v , i ) C h v a e , σ 2 ( L v , i ) , 0 i < n .

In all these estimates, h = 2 - n , C is independent of i and v, and a c and a e are the mesh grading parameters in (3.2) and (3.3).

Proof.

For ( x , y , z ) L v , i , let ( x ^ , y ^ , z ^ ) L ^ be its image under the dilation 𝐁 v , i in (3.4), where L ^ = L v , 0 when i < n and L ^ = T ( 0 ) when i = n . For a function v on L v , i , we define v ^ on L ^ by v ^ ( x ^ , y ^ , z ^ ) : = v ( x , y , z ) . As part of 𝒯 n , the triangulation on L v , i is mapped by 𝐁 v , i to a triangulation on L ^ with mesh size O ( 2 i - n ) . Then, by the scaling argument and the interpolation error estimate on L ^ , we have

(7.4) | v - I v | H 1 ( L v , i ) 2 = κ i | v ^ - I v ^ | H 1 ( L ^ ) 2 C κ i 2 2 ( i - n ) | v ^ | H 2 ( L ^ ) 2 C 2 2 ( i - n ) κ 2 i | v | H 2 ( L v , i ) 2 .

Thus, for v H 2 ( T ( 0 ) ) , by (7.4) and κ 1 2 , for 0 i n , we have

| v - I v | H 1 ( L v , i ) C h | v | H 2 ( L v , i ) ,

which proves the first estimate of the lemma.

If T ( 0 ) is a v-tetrahedron, for i < n , we have R c κ i = κ c i on L v , i . Note that if v V σ c 2 ( T ( 0 ) ) for any σ c > 1 2 - λ c (a condition in (7.3)), according to condition (3.3), it is clear that v V 1 - a c 2 ( T ( 0 ) ) . Therefore, by (7.4) and (3.1), we have

| v - I v | H 1 ( L v , i ) 2 C 2 2 ( i - n ) κ 2 i | v | H 2 ( L v , i ) 2 C 2 2 ( i - n ) κ c 2 i a c | α | = 2 R c 1 - a c α v L 2 ( L v , i ) 2 C 2 - 2 n v V 1 - a c 2 ( L i , v ) 2 .

This proves the second estimate of the lemma.

If T ( 0 ) is a v e -tetrahedron, for i < n , we have r e κ i = κ e i on L v , i . Note that if v γ e , σ 2 ( T ( 0 ) ) for any γ e < λ e (a condition in (7.2)), according to the condition (3.2), it is clear that v a e , σ 2 ( T ( 0 ) ) . Thus, by (7.4), (3.1) and (3.2), we have

| v - I v | H 1 ( L v , i ) 2 C 2 2 ( i - n ) κ 2 i | v | H 2 ( L v , i ) 2 C 2 2 ( i - n ) κ e 2 i a e | α | = 2 r e 1 - a e α v L 2 ( L v , i ) 2 C 2 - 2 n ( | α | = 2 r e 1 - a e α v L 2 ( L v , i ) 2 + | α | = 1 α z v L 2 ( L v , i ) 2 + r e - σ z 2 v L 2 ( L v , i ) 2 ) C 2 - 2 n v a e , σ 2 ( L i , v ) 2 .

This proves the third estimate and concludes the proof of this lemma. ∎

Then we give the error estimate on the whole initial tetrahedron T ( 0 ) .

Corollary 7.3.

Let T ( 0 ) T 0 be either a v- or a v e -tetrahedron. For the solution u of equation (1.1) with f given in (2.7), let Iu be its nodal interpolation on T n . Then we have | u - I u | H 1 ( T ( 0 ) ) C h , where h = 2 - n and C is independent of n.

Proof.

We first show the interpolation error estimate on the last layer L v , n .

For ( x , y , z ) L v , n , let ( x ^ , y ^ , z ^ ) T ( 0 ) be its image under the dilation 𝐁 v , n . For a function v on L v , n , we define v ^ on T ( 0 ) by v ^ ( x ^ , y ^ , z ^ ) : = v ( x , y , z ) . Now, let χ be a smooth cut-off function on T ( 0 ) such that χ = 0 in a neighborhood of x 0 and = 1 at every other node of T ( 0 ) . Recall the distance function ρ from (7.1). Thus, ρ ( x ^ , y ^ , z ^ ) = κ - n ρ ( x , y , z ) . Since χ v ^ = 0 in the neighborhood of x 0 , we have

(7.5) | χ v ^ | H 2 ( T ( 0 ) ) 2 C | α | 2 ρ | α | - 1 α v ^ L 2 ( T ( 0 ) ) 2 .

Define w ^ : = v ^ - χ v ^ , and note that I ^ ( χ v ^ ) = I ^ v ^ , where I ^ v ^ is the interpolation on T ( 0 ) . We have

(7.6) | v ^ - I ^ v ^ | H 1 ( T ( 0 ) ) = | w ^ + χ v ^ - I ^ v ^ | H 1 ( T ( 0 ) ) | w ^ | H 1 ( T ( 0 ) ) + | χ v ^ - I ^ v ^ | H 1 ( T ( 0 ) ) = | w ^ | H 1 ( T ( 0 ) ) + | χ v ^ - I ^ ( χ v ^ ) | H 1 ( T ( 0 ) ) C ( v ^ H 1 ( T ( 0 ) ) + | χ v ^ | H 2 ( T ( 0 ) ) ) ,

where C depends on T ( 0 ) . Then, using (7.6), (7.5), the scaling argument and κ - n ρ - 1 in L v , n , we have

| v - I v | H 1 ( L v , n ) 2 = κ n | v ^ - I ^ v ^ | H 1 ( T ( 0 ) ) 2 C κ n ( v ^ H 1 ( T ( 0 ) ) 2 + | α | 2 ρ | α | - 1 α v ^ L 2 ( T ( 0 ) ) 2 ) C | α | 2 ρ | α | - 1 α v L 2 ( L v , n ) 2 C κ 2 n a | α | 2 ρ | α | - 1 - a α v L 2 ( L v , n ) 2 .

When T ( 0 ) is a v-tetrahedron, by ρ R c , the definition of the weighted space, (3.1) and (3.3), we have

(7.7) | v - I v | H 1 ( L v , n ) 2 C 2 - 2 n | α | 2 R c | α | - 1 - a c α v L 2 ( L v , n ) 2 C h 2 v V 1 - a c 2 ( L v , n ) 2 .

When T ( 0 ) is a v e -tetrahedron, by ρ r e , the definition of the weighted space, (3.1) and (3.2), one obtains

(7.8) | v - I v | H 1 ( L v , n ) 2 C 2 - 2 n | α | 2 r e | α | - 1 - a e α v L 2 ( L v , n ) 2 C h 2 ( r e - 1 z v L 2 ( L v , n ) 2 + r e - σ z 2 v L 2 ( L v , n ) 2 + | α | = 1 α z v L 2 ( L v , n ) 2 + | α | 2 r e | α | - 1 - a e α v L 2 ( L v , n ) 2 ) C h 2 v a e , σ 2 ( L v , n ) 2 .

For a v-tetrahedron T ( 0 ) , recall decomposition (7.3), u = u c , reg + u c , sing . In Lemma 7.2, replace v by u c , reg in the first estimate, replace v by u c , sing in the second estimate, and replace v by u c , sing in (7.7). Summing up these estimates over all the layers, by the regularity estimates in Corollary 6.5 and by a c ( 0 , 1 ] and a c < λ c + 1 2 , we therefore have

| u - I u | H 1 ( T ( 0 ) ) | u c , reg - I u c , reg | H 1 ( T ( 0 ) ) + | u c , sing - I u c , sing | H 1 ( T ( 0 ) ) C h ( u c , reg H 2 ( T ( 0 ) ) + u c , sing V 1 - a c 2 ( T ( 0 ) ) ) C h .

Similarly, for a v e -tetrahedron T ( 0 ) , recall decomposition (7.2), u = u ξ e , reg + u ξ e , sing . In Lemma 7.2, replace v by u ξ e , reg in the first estimate, replace v by u ξ e , sing in the third estimate, and replace v by u ξ e , sing in (7.8). Summing up these estimates over all the layers, by the regularity estimates in Theorem 6.1 and by a e ( 0 , 1 ] and 1 - σ a e < λ e , we therefore have

| u - I u | H 1 ( T ( 0 ) ) | u ξ e , reg - I u ξ e , reg | H 1 ( T ( 0 ) ) + | u ξ e , sing - I u ξ e , sing | H 1 ( T ( 0 ) ) C h ( u ξ e , reg H 2 ( T ( 0 ) ) + u ξ e , sing a e , σ 2 ( T ( 0 ) ) ) C h .

Hence, the proof is completed. ∎

7.2 Analysis on Initial e-Tetrahedra

In the neighborhood of an edge e, according to Theorem 6.1, we write u = u reg + u sing . Recall the nodal interpolation I u S n on 𝒯 n . Then the interpolation error on an initial e-tetrahedron T ( 0 ) 𝒯 0 satisfies

(7.9) | u - I u | H 1 ( T ( 0 ) ) | u reg - I u reg | H 1 ( T ( 0 ) ) + | u sing - I u sing | H 1 ( T ( 0 ) ) .

Lemma 7.4.

Let T ( 0 ) T 0 be an e-tetrahedron, and let L e , i be the ith mesh layer (Definition 3.7). Then, for σ given in (2.7) and a e given in (3.2) (namely, a e ( 0 , 1 ] and 1 - σ a e < λ e ), we have

| u reg - I u reg | H 1 ( L e , i ) C h ( u reg H 2 ( L e , i ) + r e - σ z 2 u reg L 2 ( L e , i ) ) , 0 i n ,
| u sing - I u sing | H 1 ( L e , i ) C h u sing a e , σ 2 ( L e , i ) , 0 i < n ,

where n is the number of refinements, h = 2 - n , and C depends on T ( 0 ) , but not on i.

Proof.

Recall that the layer L e , n T ( 0 ) is the union of tetrahedra in 𝒯 n that touch the edge e; if i < n , L e , i is formed in the ( i + 1 ) th refinement and is the union of tetrahedra in 𝒯 i + 1 between P e , i and P e , i + 1 . Therefore, it suffices to study the interpolation error estimate on each tetrahedron T ( n ) L e , n and on 𝒯 i + 1 T ( i + 1 ) L e , i if i < n . Let T ( i ) 𝒯 i be the tetrahedron containing T ( i + 1 ) if i < n , and let T ( i ) = T ( n ) if i = n . Then T ( i ) is either an e-tetrahedron or a v e -tetrahedron. To simplify the notation, in what follows, we denote by v a function and by Iv its nodal interpolation on 𝒯 n . The analysis is based on T ( i ) ’s type.

Case I: T ( i ) is an e-tetrahedron. Let ( x , y , z ) T ( i + 1 ) , and let ( x ^ , y ^ , z ^ ) T ^ ( i + 1 ) : = 𝐁 e , i T ( i + 1 ) as defined in Proposition 3.8, where T ^ ( i + 1 ) is the reference tetrahedron. Let v ^ ( x ^ , y ^ , z ^ ) : = v ( x , y , z ) . Then, by the mapping in (3.5), we have

(7.10) { d x d y d z = 2 - i κ e 2 i d x ^ d y ^ d z ^ ; x v = ( κ e - i x ^ + b 1 κ e - i z ^ ) v ^ , y v = ( κ e - i y ^ + b 2 κ e - i z ^ ) v ^ , z v = 2 i z ^ v ^ ; x ^ v ^ = ( κ e i x - b 1 2 - i z ) v , y ^ v ^ = ( κ e i y - b 2 2 - i z ) v , z ^ v ^ = 2 - i z v .

Therefore, by (7.10) and the standard interpolation estimate on T ^ ( i + 1 ) , we have

(7.11) x ( v - I v ) L 2 ( T ( i + 1 ) ) 2 C 2 - i ( x ^ ( v ^ - I v ^ ) L 2 ( T ^ ( i + 1 ) ) 2 + z ^ ( v ^ - I v ^ ) L 2 ( T ^ ( i + 1 ) ) 2 ) C 2 - i 2 2 ( i - n ) | v ^ | H 2 ( T ^ ( i + 1 ) ) 2 C 2 2 ( i - n ) | α | + α 3 = 2 2 - 2 i α 3 κ e 2 i ( | α | - 1 ) α z α 3 v L 2 ( T ( i + 1 ) ) 2 .

A similar calculation for the derivative with respect to y gives

(7.12) y ( v - I v ) L 2 ( T ( i + 1 ) ) 2 C 2 2 ( i - n ) | α | + α 3 = 2 2 - 2 i α 3 κ e 2 i ( | α | - 1 ) α z α 3 v L 2 ( T ( i + 1 ) ) 2 .

In the z-direction, by (7.10) and κ e 1 2 , following the calculation in (7.11), we have

(7.13) z ( v - I v ) L 2 ( T ( i + 1 ) ) 2 C 2 i κ e 2 i z ^ ( v ^ - I v ^ ) L 2 ( T ^ ( i + 1 ) ) 2 C 2 - i ( x ^ ( v ^ - I v ^ ) L 2 ( T ^ ( i + 1 ) ) 2 + z ^ ( v ^ - I v ^ ) L 2 ( T ^ ( i + 1 ) ) 2 ) C 2 2 ( i - n ) | α | + α 3 = 2 2 - 2 i α 3 κ e 2 i ( | α | - 1 ) α z α 3 v L 2 ( T ( i + 1 ) ) 2 .

Thus, by (7.11)–(7.13), the estimate of the term

E : = 2 2 ( i - n ) | α | + α 3 = 2 2 - 2 i α 3 κ e 2 i ( | α | - 1 ) α z α 3 v L 2 ( T ( i + 1 ) ) 2

is important for the error analysis. By (3.1) and (3.2), we first have

(7.14) E C 2 - 2 n ( 2 - 2 i κ e - 2 i z 2 v L 2 ( T ( i + 1 ) ) 2 + | α | = 1 α z v L 2 ( T ( i + 1 ) ) 2 + 2 2 i κ e 2 i | α | = 2 α v L 2 ( T ( i + 1 ) ) 2 ) C h 2 ( κ e 2 i ( a e - 1 ) z 2 v L 2 ( T ( i + 1 ) ) 2 + | α | = 1 α z v L 2 ( T ( i + 1 ) ) 2 + κ e 2 i ( 1 - a e ) | α | = 2 α v L 2 ( T ( i + 1 ) ) 2 ) .

Recall κ e = 2 - 1 a e , a e 1 , r e κ e i on T ( i + 1 ) if i < n , and r e < C κ e n on T ( n ) . By 1 - σ a e < λ e , for i < n , we have

(7.15) E C h 2 ( r e a e - 1 z 2 v L 2 ( T ( i + 1 ) ) 2 + | α | = 1 α z v L 2 ( T ( i + 1 ) ) 2 + | α | = 2 r e 1 - a e α v L 2 ( T ( i + 1 ) ) 2 ) C h 2 ( r e - σ z 2 v L 2 ( T ( i + 1 ) ) 2 + | α | = 1 α z v L 2 ( T ( i + 1 ) ) 2 + | α | = 2 r e 1 - a e α v L 2 ( T ( i + 1 ) ) 2 ) .

For i = n , we have

(7.16) E C h 2 ( r e a e - 1 z 2 v L 2 ( T ( n ) ) 2 + | α | = 1 α z v L 2 ( T ( n ) ) 2 + κ e 2 i ( 1 - a e ) | α | = 2 α v L 2 ( T ( n ) ) 2 ) C h 2 ( r e - σ z 2 v L 2 ( T ( n ) ) 2 + | α | = 1 α z v L 2 ( T ( n ) ) 2 + | α | = 2 α v L 2 ( T ( n ) ) 2 ) .

Thus, for v = u reg , by (7.11)–(7.16), for i n , we have

(7.17) | u reg - I u reg | H 1 ( T ( i + 1 ) ) C h ( u reg H 2 ( T ( i + 1 ) ) + r e - σ z 2 u reg L 2 ( T ( i + 1 ) ) ) .

For v = u sing , by (7.11)–(7.15), for i < n , we have

(7.18) | u sing - I u sing | H 1 ( T ( i + 1 ) ) C h u sing a e , σ 2 ( T ( i + 1 ) ) .

Hence, we have completed the proof of Case I. Note that, for i = n , the scaling arguments in (7.11)–(7.13) hold for u reg since u ^ reg H 2 on the reference element. They hold for u sing only when i < n because u sing H 2 on the last layer of the mesh. This is why estimates (7.17) and (7.18) have different forms.

Case II: T ( i ) is a v e -tetrahedron. Let T ( k ) 𝒯 k , 1 k i , be the v e -tetrahedron such that T ( i ) T ( k ) and T ( k ) is contained in an e-tetrahedron T ( k - 1 ) 𝒯 k - 1 . For ( x , y , z ) T ( i + 1 ) , let ( x ^ , y ^ , z ^ ) T ^ ( i + 1 ) = 𝐁 i , k T ( i + 1 ) (see Proposition 3.8), where T ^ ( i + 1 ) is the reference element. Let v ^ ( x ^ , y ^ , z ^ ) : = v ( x , y , z ) . We have

(7.19) { d x d y d z = 2 1 - k κ e 3 i - k - 2 d x ^ d y ^ d z ^ ; x ^ v ^ = ( κ e i - 1 x - b 1 2 1 - k κ e i - k z ) v , y ^ v ^ = ( κ e i - 1 y - b 2 2 1 - k κ e i - k z ) v , z ^ v ^ = 2 1 - k κ e i - k z v ; x v = ( κ e 1 - i x ^ + b 1 κ e 1 - i z ^ ) v ^ , y v = ( κ e 1 - i y ^ + b 2 κ e 1 - i z ^ ) v ^ , z v = 2 k - 1 κ e k - i z ^ v ^ .

Therefore, by (7.19), κ e 1 2 and the standard interpolation estimate, we have

(7.20) x ( v - I v ) L 2 ( T ( i + 1 ) ) 2 C 2 1 - k κ e i - k ( x ^ ( v ^ - I v ^ ) L 2 ( T ^ ( i + 1 ) ) 2 + z ^ ( v ^ - I v ^ ) L 2 ( T ^ ( i + 1 ) ) 2 ) C 2 1 - k κ e i - k 2 2 ( i - n ) | v ^ | H 2 ( T ^ ( i + 1 ) ) 2 C 2 2 ( i - n ) | α | + α 3 = 2 2 2 ( 1 - k ) α 3 κ e 2 ( i - k ) α 3 κ e ( 2 i - 2 ) ( | α | - 1 ) α z α 3 v L 2 ( T ( i + 1 ) ) 2 C 2 2 ( i - n ) | α | + α 3 = 2 2 - 2 i α 3 κ e 2 i ( | α | - 1 ) α z α 3 v L 2 ( T ( i + 1 ) ) 2 .

A similar calculation for the derivative with respect to y gives

(7.21) y ( v - I v ) L 2 ( T ( i + 1 ) ) 2 C 2 2 ( i - n ) | α | + α 3 = 2 2 - 2 i α 3 κ e 2 i ( | α | - 1 ) α z α 3 v L 2 ( T ( i + 1 ) ) 2 .

In the z-direction, by Proposition 3.8, (7.19), κ e 1 2 and following the calculation in (7.20),

(7.22) z ( v - I v ) L 2 ( T ( i + 1 ) ) 2 C ( 2 1 - k κ e i - k ) κ e 2 ( i - 1 ) ( 2 k - 1 κ e k - i ) 2 z ^ ( v ^ - I v ^ ) L 2 ( T ^ ( i + 1 ) ) 2 C 2 1 - k κ e i - k ( x ^ ( v ^ - I v ^ ) L 2 ( T ^ ( i + 1 ) ) 2 + z ^ ( v ^ - I v ^ ) L 2 ( T ^ ( i + 1 ) ) 2 ) C 2 2 ( i - n ) | α | + α 3 = 2 2 - 2 i α 3 κ e 2 i ( | α | - 1 ) α z α 3 v L 2 ( T ( i + 1 ) ) 2 .

Then, by (7.20)–(7.22) and (7.15)–(7.16), we have obtained the desired estimates for Case II.

Thus, we complete the proof by summing up the estimates for all the tetrahedra T ( i + 1 ) in L e , i . ∎

Then we have the interpolation error analysis on an initial e-tetrahedron.

Theorem 7.5.

Let T ( 0 ) T 0 be an e-tetrahedron. Recall the decomposition u = u reg + u sing on T ( 0 ) from Theorem 6.1. Let Iu be its nodal interpolation on T n . Then we have

| u - I u | H 1 ( T ( 0 ) ) C h ( u sing a e , σ 2 ( T ( 0 ) ) + u reg H 2 + σ ( T ( 0 ) ) + r e - σ z 2 u reg L 2 ( T ( 0 ) ) ) ,

where h = 2 - n and C depends on T ( 0 ) but not on n.

Proof.

By (7.9) and Lemma 7.4, it suffices to show

| u sing - I u sing | H 1 ( T ( n ) ) C h u sing a e , σ 2 ( T ( n ) )

for any tetrahedron T ( n ) 𝒯 n in the last layer L e , n . Since T ( n ) is either an e- or a v e -tetrahedron, we derive this estimate in two cases. To simplify the notation, we let v = u sing .

Case I: T ( n ) is an e-tetrahedron. By Proposition 3.8, the mapping 𝐁 e , n translates T ( n ) to the reference tetrahedron T ^ = 4 x ^ 0 x ^ 1 x ^ 2 x ^ 3 . Consequently, it maps any point ( x , y , z ) T ( n ) to ( x ^ , y ^ , z ^ ) T ^ . For a function v on T ( n ) , we define v ^ on T ^ by v ^ ( x ^ , y ^ , z ^ ) : = v ( x , y , z ) . Now, let χ be a smooth cut-off function on T ^ such that χ = 0 in a neighborhood of the edge e ^ : = x ^ 0 x ^ 1 and = 1 at every other Lagrange node of T ^ . Let r e ^ be the distance to e ^ . Let I ^ v ^ be the interpolation of u ^ on the reference tetrahedron T ^ . Since χ v ^ = 0 in the neighborhood of e ^ , I ^ ( χ v ^ ) = I ^ v ^ ( v ^ = 0 on e ^ since v γ e , σ 2 ( T ( 0 ) ) (Theorem 6.1)), and we have

(7.23) | χ v ^ | H 2 ( T ^ ) 2 C ( r e ^ a e - 1 z ^ 2 v ^ L 2 ( T ^ ) 2 + | α | + α 3 2 , α 3 < 2 r e ^ | α | - 1 α z ^ α 3 v ^ L 2 ( T ^ ) 2 ) .

Define w ^ : = v ^ - χ v ^ . Then, by the usual interpolation error estimate, we have

(7.24) | v ^ - I v ^ | H 1 ( T ^ ) = | w ^ + χ v ^ - I v ^ | H 1 ( T ^ ) | w ^ | H 1 ( T ^ ) + | χ v ^ - I v ^ | H 1 ( T ^ ) = | w ^ | H 1 ( T ^ ) + | χ v ^ - I ^ ( χ v ^ ) | H 1 ( T ^ ) C ( v ^ H 1 ( T ^ ) + | χ v ^ | H 2 ( T ^ ) ) ,

where C depends on, through χ, the nodes on T ^ . Then, using the scaling argument based on (7.10), by (7.24), (7.23), the relation r e ^ ( x ^ , y ^ , z ^ ) = κ e - n r e ( x , y , z ) and (3.1), we have

(7.25) x ( v - I v ) L 2 ( T ( n ) ) 2 C 2 - n ( x ^ ( v ^ - I v ^ ) L 2 ( T ^ ) 2 + z ^ ( v ^ - I v ^ ) L 2 ( T ^ ) 2 ) C 2 - n ( r e ^ a e - 1 z ^ 2 v ^ L 2 ( T ^ ) 2 + | α | + α 3 2 , α 3 < 2 r e ^ | α | - 1 α z ^ α 3 v ^ L 2 ( T ^ ) 2 ) C ( 2 - 2 n r e a e - 1 z 2 v L 2 ( T ( n ) ) 2 + | α | + α 3 2 , α 3 < 2 2 - 2 n α 3 r e | α | - 1 α z α 3 v L 2 ( T ( n ) ) 2 ) C ( 2 - 2 n ( r e a e - 1 z 2 v L 2 ( T ( n ) ) 2 + | α | = 1 α z v L 2 ( T ( n ) ) 2 + r e - 1 z v L 2 ( T ( n ) ) 2 ) + | α | 2 r e | α | - 1 α v L 2 ( T ( n ) ) 2 ) C 2 - 2 n ( r e a e - 1 z 2 v L 2 ( T ( n ) ) 2 + | α | = 1 α z v L 2 ( T ( n ) ) 2 + r e - 1 z v L 2 ( T ( n ) ) 2 + | α | 2 r e | α | - 1 - a e α v L 2 ( T ( n ) ) 2 ) .

A similar calculation for the derivative with respect to y gives

(7.26) y ( v - I v ) L 2 ( T ( n ) ) 2 C 2 - 2 n ( r e a e - 1 z 2 v L 2 ( T ( n ) ) 2 + | α | = 1 α z v L 2 ( T ( n ) ) 2 + r e - 1 z v L 2 ( T ( n ) ) 2 + | α | 2 r e | α | - 1 - a e α v L 2 ( T ( n ) ) 2 ) .

In the z-direction, using (7.24), (7.23), (7.10), (3.1) and the calculation in (7.25), we have

(7.27) z ( v - I v ) L 2 ( T ( n ) ) 2 = 2 n κ e 2 n z ^ ( v ^ - I v ^ ) L 2 ( T ^ ) 2 2 - n ( x ^ ( v ^ - I v ^ ) L 2 ( T ^ ) 2 + z ^ ( v ^ - I v ^ ) L 2 ( T ^ ) 2 ) C 2 - 2 n ( r e a e - 1 z 2 v L 2 ( T ( n ) ) 2 + | α | = 1 α z v L 2 ( T ( n ) ) 2 + r e - 1 z v L 2 ( T ( n ) ) 2 + | α | 2 r e | α | - 1 - a e α v L 2 ( T ( n ) ) 2 ) .

Recall a e - 1 - σ and a e < λ e . Then, by (7.25)–(7.27), we have

| v - I v | H 1 ( T ( n ) ) 2 C 2 - 2 n ( r e - σ z 2 v L 2 ( T ( n ) ) 2 + | α | = 1 α z v L 2 ( T ( n ) ) 2 + r e - 1 z v L 2 ( T ( n ) ) 2 + | α | 2 r e | α | - 1 - a e α v L 2 ( T ( n ) ) 2 ) C h 2 v a e , σ 2 ( T ( n ) ) 2 ,

which proves the estimate for Case I.

Case II: T ( n ) is a v e -tetrahedron. Let T ( k ) 𝒯 k , 1 k n , be the v e -tetrahedron such that T ( n ) T ( k ) and T ( k ) is contained in an e-tetrahedron T ( k - 1 ) 𝒯 k - 1 . By Proposition 3.8, the mapping 𝐁 n , k transforms T ( n ) to a v e -tetrahedron T ^ ( n ) 𝒯 ^ 1 . Thus, 𝐁 n , k maps every point ( x , y , z ) T ( n ) to ( x ^ , y ^ , z ^ ) T ^ ( n ) . As in Case I, for a function v on T ( n ) , we define v ^ on T ^ ( n ) by v ^ ( x ^ , y ^ , z ^ ) : = v ( x , y , z ) . Now, let χ be a smooth cut-off function on T ^ ( n ) such that χ = 0 in a neighborhood of e ^ : = x ^ 0 x ^ 1 of T ^ = 4 x ^ 0 x ^ 1 x ^ 2 x ^ 3 and χ = 1 at every other Lagrange node of T ^ ( n ) . Recall the distance r e ^ to e ^ . Since χ v ^ = 0 in the neighborhood of the refined vertex, we have I ^ ( χ v ^ ) = I ^ v ^ on T ^ ( n ) and

(7.28) | χ v ^ | H 2 ( T ^ ( n ) ) 2 C ( r e ^ a e - 1 z ^ 2 v ^ L 2 ( T ^ ( n ) ) 2 + | α | + α 3 2 , α 3 < 2 r e ^ | α | - 1 α z ^ α 3 v ^ L 2 ( T ^ ( n ) ) 2 ) .

Define w ^ : = v ^ - χ v ^ . Then, by the usual interpolation error estimate, we have

(7.29) | v ^ - I v ^ | H 1 ( T ^ ( n ) ) = | w ^ + χ v ^ - I v ^ | H 1 ( T ^ ( n ) ) | w ^ | H 1 ( T ^ ( n ) ) + | χ v ^ - I v ^ | H 1 ( T ^ ( n ) ) = | w ^ | H 1 ( T ^ ( n ) ) + | χ v ^ - I ^ ( χ v ^ ) | H 1 ( T ^ ( n ) ) C ( v ^ H 1 ( T ^ ( n ) ) + | χ v ^ | H 2 ( T ^ ( n ) ) ) ,

where C depends on, through χ, the nodes in the reference element T ^ ( n ) . In L e , n , r e ( x , y , z ) = κ e n - 1 r e ^ ( x ^ , y ^ , z ^ ) . Therefore, by (7.19), (7.29), (7.28), (3.1) and κ e 1 2 , we have

(7.30) x ( v - I v ) L 2 ( T ( n ) ) 2 C 2 1 - k κ e n - k ( x ^ ( v ^ - I v ^ ) L 2 ( T ^ ( n ) ) 2 + z ^ ( v ^ - I v ^ ) L 2 ( T ^ ( n ) ) 2 ) C 2 1 - k κ e n - k ( r e ^ a e - 1 z ^ 2 v ^ L 2 ( T ^ ( n ) ) 2 + | α | + α 3 2 , α 3 < 2 r e ^ | α | - 1 α z ^ α 3 v ^ L 2 ( T ^ ( n ) ) 2 ) C ( 2 2 ( 1 - k ) κ e 2 ( n - k ) r e a e - 1 z 2 v L 2 ( T ( n ) ) 2 + | α | + α 3 2 , α 3 < 2 2 2 ( 1 - k ) α 3 κ e 2 ( n - k ) α 3 r e | α | - 1 α z α 3 v L 2 ( T ( n ) ) 2 ) C 2 - 2 n ( r e a e - 1 z 2 v L 2 ( T ( n ) ) 2 + | α | = 1 α z v L 2 ( T ( n ) ) 2 + r e - 1 z v L 2 ( T ( n ) ) 2 + | α | 2 r e | α | - 1 - a e α v L 2 ( T ( n ) ) 2 ) .

A similar calculation for the derivative with respect to y gives

(7.31) y ( v - I v ) L 2 ( T ( n ) ) 2 C 2 - 2 n ( r e a e - 1 z 2 v L 2 ( T ( n ) ) 2 + | α | = 1 α z v L 2 ( T ( n ) ) 2 + r e - 1 z v L 2 ( T ( n ) ) 2 + | α | 2 r e | α | - 1 - a e α v L 2 ( T ( n ) ) 2 ) .

In the z-direction, by (7.19), (7.29), (7.28) and the calculation in (7.30), we have

(7.32) z ( v - I v ) L 2 ( T ( n ) ) 2 = ( 2 1 - k κ e n - k ) κ e 2 ( n - 1 ) ( 2 k - 1 κ e k - n ) 2 z ^ ( v ^ - I v ^ ) L 2 ( T ^ ( n ) ) 2 2 1 - k κ e n - k ( x ^ ( v ^ - I v ^ ) L 2 ( T ^ ( n ) ) 2 + z ^ ( v ^ - I v ^ ) L 2 ( T ^ ( n ) ) 2 ) C 2 - 2 n ( r e a e - 1 z 2 v L 2 ( T ( n ) ) 2 + | α | = 1 α z v L 2 ( T ( n ) ) 2 + r e - 1 z v L 2 ( T ( n ) ) 2 + | α | 2 r e | α | - 1 - a e α v L 2 ( T ( n ) ) 2 ) .

Therefore, by a e - 1 - σ and a e < λ e and by (7.30)–(7.32), we have

| v - I v | H 1 ( T ( n ) ) 2 C h 2 v a e , σ 2 ( T ( n ) ) 2 ,

which proves the estimate for Case II. Hence, the theorem is proved by summing up the estimates in Lemma 7.4 and the estimates | u sing - I u sing | H 1 ( T ( n ) ) for all the tetrahedra T ( n ) in L e , n . ∎

Based on the regularity estimates in Theorem 6.1 and the interpolation error estimate in Theorem 7.5, it is clear that, on an initial e-tetrahedron T ( 0 ) , | u - I u | H 1 ( T ( 0 ) ) C h , where h = 2 - n and C is independent of n.

7.3 Analysis on Initial ev-Tetrahedra

In the neighborhood 𝒱 c e of an edge e and a vertex c, according to Theorem 6.2, we write u = u reg + u sing + u c , where u c = - 1 2 < λ c , k < σ + 1 2 c k ψ σ c , k . Then the interpolation error on an initial ev-tetrahedron T ( 0 ) 𝒯 0 satisfies

(7.33) | u - I u | H 1 ( T ( 0 ) ) | u reg - I u reg | H 1 ( T ( 0 ) ) + | u sing - I u sing | H 1 ( T ( 0 ) ) + | u c - I u c | H 1 ( T ( 0 ) ) .

Theorem 7.6.

Let T ( 0 ) T 0 be an ev-tetrahedron, and let L e v , i be its ith mesh layer (Definition 3.9). Recall the weighted space H 𝛄 , σ 2 specified in Theorem 6.2. Then, for σ given in (2.7), a e and a c given in (3.2) and (3.3), respectively, and 0 i n , in (7.33), we have

| u reg - I u reg | H 1 ( L e , i ) C h ( u reg V 0 2 ( L e , i ) + r e - σ z 2 u reg L 2 ( L e , i ) ) ,
| u sing - I u sing | H 1 ( T ( 0 ) ) C h u sing 𝜸 , σ 2 ( T ( 0 ) ) ,
| u c - I u c | H 1 ( T ( 0 ) ) C h ,

where n is the number of refinements, h = 2 - n ; in the first two estimates, C depends on T ( 0 ) , but not on i; in the last estimate, C depends on T ( 0 ) and f in (1.1).

Proof.

Let T ( i ) T ( 0 ) be the ev-tetrahedron in 𝒯 i . Recall the mesh layer L e v , i in Definition 3.9 and the mapping 𝐁 e v , i in (3.7) that translates L e v , i to the reference domain L ^ . For a point ( x , y , z ) L e v , i , let ( x ^ , y ^ , z ^ ) L ^ be its image under 𝐁 e v , i . For a function v on L e v , i , define the function v ^ on L ^ by

(7.34) v ^ ( x ^ , y ^ , z ^ ) : = v ( x , y , z ) .

The distance function r e to the edge e satisfies r e ( x , y , z ) = κ c i r e ( x ^ , y ^ , z ^ ) on L e v , i . Meanwhile, 𝐁 e v , i maps the triangulation 𝒯 n on L e v , i to a graded triangulation on L ^ that is obtained after i + 1 - n refinements of the initial mesh ^ . Note that the subsequent refinements on ^ are anisotropic with the parameter κ e toward e since ^ does not contain ev- or v-tetrahedra. Then, by the scaling argument using the mapping (3.7), we have

(7.35) | v - I v | H 1 ( L e v , i ) 2 C κ c i | v ^ - I v ^ | H 1 ( L ^ ) 2 .

Case I: v = u reg V 0 2 + σ ( T ( 0 ) ) and r e - σ z 2 v L 2 ( T ( 0 ) ) . Note that this implies v H 2 ( L e v , i ) for i n . Then, by (7.35) and by similar calculations as in (7.15) and (7.16), and 1 - σ a e , we have

(7.36) | v - I v | H 1 ( L e v , i ) 2 C κ c i 2 2 ( i - n ) ( r e a e - 1 z ^ 2 v ^ L 2 ( L ^ ) 2 + | α | = 1 α z ^ v ^ L 2 ( L ^ ) 2 + | α | = 2 α v ^ L 2 ( L ^ ) 2 ) C 2 - 2 n ( 2 2 i κ c 4 i - 2 i a e r e a e - 1 z 2 v L 2 ( L e v , i ) 2 + 2 2 i κ c 2 i | α | = 1 α z v L 2 ( L e v , i ) 2 + 2 2 i κ c 2 i | α | = 2 α v L 2 ( L e v , i ) 2 ) C 2 - 2 n ( r e a e - 1 z 2 v L 2 ( L e v , i ) 2 + | α | = 2 α v L 2 ( L e v , i ) 2 ) C 2 - 2 n ( r e - σ z 2 v L 2 ( L e v , i ) 2 + | α | = 2 α v L 2 ( L e v , i ) 2 ) C h 2 ( v V 0 2 ( L e , i ) 2 + r e - σ z 2 v L 2 ( L e , i ) 2 ) .

This proves the first estimate in the theorem.

Case II: v = u sing H 𝛄 , σ 2 ( T ( 0 ) ) with γ c = 1 + σ and γ e < λ e . The following estimate was obtained in [22, Corollary 5.16] for functions in 𝜸 * , μ 2 ( T ( 0 ) ) , where γ e * , γ c * [ 0 , 1 ] and γ e * < λ e , γ c * < λ c + 1 2 and 1 - μ a e :

(7.37) | v - I v | H 1 ( T ( 0 ) ) C h v 𝜸 * , μ 2 ( T ( 0 ) ) .

Since 𝜸 , σ 2 ( T ( 0 ) ) 𝜸 * , μ 2 ( T ( 0 ) ) , we have

| v - I v | H 1 ( T ( 0 ) ) C h v 𝜸 * , μ 2 ( T ( 0 ) ) C h v 𝜸 , σ 2 ( T ( 0 ) ) .

This proves the second estimate in the theorem.

Case III: v = u c in (7.33). Recall that if zero is an eigenvalue of the Laplace–Beltrami operator in the expansion of u c , namely λ c , 0 = 0 , the corresponding term is a constant function, which the finite element interpolation can completely resolve. Therefore, we proceed to consider the interpolation error in the case that λ c , 0 > 0 . Recall from Lemma 6.3 that v = v 1 + v 2 . Therefore, to obtain | v - I v | H 1 ( T 0 ) in this case, it is sufficient to analyze | v 1 - I v 1 | H 1 ( T 0 ) and | v 2 - I v 2 | H 1 ( T 0 ) , respectively.

We first estimate | v 1 - I v 1 | H 1 ( T 0 ) . For i < n , by the scaling map (3.7), R c κ c i , (3.1), and the estimates in (7.20)–(7.22) and (7.15)–(7.16), we have

| v 1 - I v 1 | H 1 ( L e v , i ) 2 C κ c i 2 2 ( i - n ) ( r e a e - 1 z ^ 2 v ^ 1 L 2 ( L ^ ) 2 + | α | = 1 α z ^ v ^ 1 L 2 ( L ^ ) 2 + | α | = 2 α v ^ 1 L 2 ( L ^ ) 2 ) C 2 - 2 n ( 2 2 i κ c 4 i - 2 i a e r e a e - 1 z 2 v 1 L 2 ( L e v , i ) 2 + 2 2 i κ c 2 i | α | = 1 α z v 1 L 2 ( L e v , i ) 2 + 2 2 i κ c 2 i | α | = 2 α v 1 L 2 ( L e v , i ) 2 ) C 2 - 2 n ( R c 2 - a c - a e r e a e - 1 z 2 v 1 L 2 ( L e v , i ) 2 + | α | = 1 R c 1 - a c α z v 1 L 2 ( L e v , i ) 2 + | α | = 2 R c 1 - a c α v 1 L 2 ( L e v , i ) 2 ) .

According to Lemma 6.3, R c σ c + σ c , e r e - σ c , e z 2 v 1 L 2 ( 𝒱 c e ) , R c σ c α z v 1 L 2 ( 𝒱 c e ) ( | α | = 1 ), and R σ c α v 1 L 2 ( 𝒱 c e ) ( | α | = 2 ), where σ c , e ( 0 , 1 ) and any σ c > 1 2 - λ c , 0 . Choose σ c , e > 1 - a e and σ c = 1 - a c . Then we have

(7.38) | v 1 - I v 1 | H 1 ( L e v , i ) 2 C 2 - 2 n ( R c 1 - a c + σ c , e r e - σ c , e z 2 v 1 L 2 ( L e v , i ) 2 + | α | = 1 R c 1 - a c α z v 1 L 2 ( L e v , i ) 2 + | α | = 2 R c 1 - a c α v 1 L 2 ( L e v , i ) 2 ) C 2 - 2 n ( R c σ c + σ c , e r e - σ c , e z 2 v 1 L 2 ( L e v , i ) 2 + | α | = 1 R c σ c α z v 1 L 2 ( L e v , i ) 2 + | α | = 2 R c σ c α v 1 L 2 ( L e v , i ) 2 ) .

For i = n , 𝐁 e v , n ( L e v , n ) = L ^ = T ( 0 ) . With λ c > 0 , recall the condition R c σ c - 2 v 1 L 2 ( 𝒱 e c ) for any σ c > 1 2 - λ c in Lemma 6.3. This implies v 1 ( c ) = 0 . For ( x , y , z ) L e v , n , let ( x ^ , y ^ , z ^ ) L ^ be its image under 𝐁 e v , n . For a function v on L e v , n , recall the scaling (7.34) v ^ on L ^ .

Now, let χ be a smooth cut-off function on L ^ such that χ = 0 in a neighborhood of the vertex c and χ = 1 at every other node of L ^ . Let I ^ v ^ be the interpolation of v ^ on the reference tetrahedron L ^ . Since χ v ^ = 0 in the neighborhood of c, I ^ ( χ v ^ ) = I ^ v ^ I = I v ^ . Therefore, by (7.35), R c C κ c n , R c ( x ^ , y ^ , z ^ ) = κ c - n R c ( x , y , z ) and r e ( x ^ , y ^ , z ^ ) = κ c - n r e ( x , y , z ) , we have

| v 1 - I v 1 | H 1 ( L e v , n ) 2 C κ c n | v ^ 1 - I v 1 ^ | H 1 ( L ^ ) 2 C κ c n ( | v ^ 1 - χ v ^ 1 | H 1 ( L ^ ) 2 + | χ v ^ 1 - I ^ ( χ v ^ 1 ) | H 1 ( L ^ ) 2 )
C κ c n ( v ^ 1 H 1 ( L ^ ) 2 + | χ v ^ 1 | H 2 ( L ^ ) 2 )
C κ c n ( v ^ 1 H 1 ( L ^ ) 2 + | α | = 2 , α 3 1 R c 1 - a c α z ^ α 3 v ^ 1 L 2 ( L ^ ) 2 + R c 1 - a c + σ c , e r e - σ c , e z ^ 2 v ^ 1 L 2 ( L ^ ) 2 )
C 2 - 2 n ( 2 2 n κ c - 2 n v 1 L 2 ( L e v , n ) 2 + 2 2 n | v 1 | H 1 ( L e v , n ) 2 + 2 2 n κ c 2 n ( κ c 2 n a c - 2 n )
       × ( | α | = 2 R c 1 - a c α v 1 L 2 ( L e v , n ) 2 + | α | = 1 R c 1 - a c α z v 1 L 2 ( L e v , n ) 2
          + R c 1 - a c + σ c , e r e - σ c , e z 2 v 1 L 2 ( L e v , n ) 2 ) )
C 2 - 2 n ( R c - 1 - a c v 1 L 2 ( L e v , n ) 2 + | α | = 1 R c - a c α v 1 L 2 ( L e v , n ) 2
    + | α | = 2 R c 1 - a c α v 1 L 2 ( L e v , n ) 2 + | α | = 1 R c 1 - a c α z v 1 L 2 ( L e v , n ) 2
    + R c 1 - a c + σ c , e r e - σ c , e z 2 v 1 L 2 ( L e v , n ) 2 )
C 2 - 2 n ( R c σ c - 2 v 1 L 2 ( L e v , n ) 2 + | α | = 1 R c σ c - 1 α v 1 L 2 ( L e v , n ) 2
    + | α | = 2 R c σ c α v 1 L 2 ( L e v , n ) 2 + | α | = 1 R c σ c α z v 1 L 2 ( L e v , n ) 2
(7.39)     + R c σ c + σ c , e r e - σ c , e z 2 v 1 L 2 ( L e v , n ) 2 ) ,

where we chose σ c = 1 - a c and any σ c , e > 1 - a e . Therefore, by (7.38), (7.39) and Lemma 6.3, we have

(7.40) | v 1 - I v 1 | H 1 ( T ( 0 ) ) 2 C h 2 .

Now, for | v 2 - I v 2 | H 1 ( T 0 ) . Recall from Remark 6.4 that v 2 satisfies v 2 𝜸 , σ 2 ( 𝒱 c e ) with γ e = a e , γ c = a c and σ = σ c , e . Choose σ c , e 1 - a e ( 0 , 1 ) . Then, based on the estimates in [22, Corollary 5.16] and in Lemma 6.3, we have

(7.41) | v 2 - I v 2 | H 1 ( T ( 0 ) ) C h v 2 𝜸 , σ 2 ( T ( 0 ) ) C h .

Thus, the theorem is proved by the estimates in (7.36), (7.37), (7.40) and (7.41). ∎

Hence, based on the interpolation error estimates for different initial tetrahedra, we obtain the global error analysis.

Corollary 7.7.

For the anisotropic finite element method proposed in Algorithm 3.4 solving equation (1.1) with f H σ ( Ω ) and σ [ 0 , 1 ) satisfying (2.7), σ λ e , k - 1 , for all k N * , e E (see Theorem 6.1), and σ λ c , k - 1 2 for all k N , c C (see Theorem 6.2), we have u - u n H 1 ( Ω ) C h , where h = 2 - n , and C depends on the initial triangulation T 0 and f, but not on n.

Proof.

We first have the optimal interpolation error estimate | u - I u | H 1 ( Ω ) C h by summing up the interpolation error estimates Lemma 7.1 (o-tetrahedra), Corollary 7.3 (v- and v e -tetrahedra), Theorem 7.5 (e-tetrahedra) and Theorem 7.6 (ev-tetrahedra), and by the regularity estimates in Section 6. The desired estimate for u - u n H 1 ( Ω ) is a consequence of the Poincaré inequality and Céa’s lemma. ∎

8 Numerical Results

In this section, we present numerical test results to verify the error analysis in Corollary 7.7 for the proposed anisotropic finite element method (Algorithm 3.4) solving elliptic equations. Our numerical tests are implemented in two typical polyhedral domains: the prism domain and the Fichera domain. We shall demonstrate three numerical examples (Examples 8.18.3) for the prism domain and one example (Example 8.4) for the Fichera domain.

The first set of tests are for the prism domain. Let Ω 1 be the square with vertices at ( 1 , 0 ) , ( 0 , 1 ) , ( - 1 , 0 ) , ( 0 , - 1 ) , and let Ω 2 be the triangle with vertices at ( 0 , 0 ) , ( - 1 , 0 ) , ( 0 , - 1 ) . Define the prism domain Ω p = ( Ω 1 Ω 2 ) × ( 0 , 1 ) . For a point ( x , y , z ) Ω , denote by ( r , θ ) the polar coordinates of its projection in the xy-plane ( r ( x , y , z ) = r ( x , y ) and θ ( x , y , z ) = θ ( x , y ) ). In the first two numerical examples (Examples 8.1 and 8.2), we are especially interested in the performance of the numerical method when the Neumann condition is imposed on both adjacent faces of the singular edge. For pure Dirichlet problems, we refer to [21, 22]. In the third example (Example 8.3), we illustrate our method with the presence of a DN singular edge.

In the first two examples, consider the elliptic equation with the mixed boundary condition

(8.1) { - Δ u = 1 in Ω p , u = r 2 3 cos 2 ( θ + π 2 ) 3 on Γ Dir , n u = 0 on Γ Neu ,

where Γ Dir and Γ Neu consist of boundary faces of the polyhedron Ω p . By imposing Dirichlet and Neumann boundary conditions on different faces, we keep the edge e = { ( 0 , 0 , z ) for  0 < z < 1 } as the singular Neumann edge. We will use the graded mesh toward the edge e and its two endpoints, as described in Algorithm 3.4, in the finite element approximation. See Figure 3 for the case κ e = κ c = 0.2 .

Figure 3 
               Graded meshes on the prism domain (left to right): the initial mesh, mesh after one refinement, mesh after three refinements (
                     
                        
                           
                              
                                 κ
                                 e
                              
                              =
                              
                                 κ
                                 c
                              
                              =
                              0.2
                           
                        
                        
                        {\kappa_{e}=\kappa_{c}=0.2}
                     
                  ).
Figure 3 
               Graded meshes on the prism domain (left to right): the initial mesh, mesh after one refinement, mesh after three refinements (
                     
                        
                           
                              
                                 κ
                                 e
                              
                              =
                              
                                 κ
                                 c
                              
                              =
                              0.2
                           
                        
                        
                        {\kappa_{e}=\kappa_{c}=0.2}
                     
                  ).
Figure 3 
               Graded meshes on the prism domain (left to right): the initial mesh, mesh after one refinement, mesh after three refinements (
                     
                        
                           
                              
                                 κ
                                 e
                              
                              =
                              
                                 κ
                                 c
                              
                              =
                              0.2
                           
                        
                        
                        {\kappa_{e}=\kappa_{c}=0.2}
                     
                  ).
Figure 3

Graded meshes on the prism domain (left to right): the initial mesh, mesh after one refinement, mesh after three refinements ( κ e = κ c = 0.2 ).

Example 8.1.

We impose the Neumann condition on the two faces adjacent to the edge e and impose the Dirichlet condition on all the other faces (including the top and bottom faces) in equation (8.1). See Figure 4. Thus, e is a Neumann edge, and other edges are either Dirichlet edges or DN edges (see the description before (2.4)).

According to (2.4) (see also [12]), the edge e is the only singular edge with λ e = 2 3 , and the two vertices c (its two endpoints) satisfy λ c + 1 2 > λ e . Other vertices and edges of Ω p are regular in this case. In addition, the right-hand side function in (8.1) belongs to H σ ( Ω p ) for any σ [ 0 , 1 ) . In fact, the solution u H 1 + s ( Ω p ) for s < 2 3 . Therefore, based on the conditions in (3.2) and (3.3), and by Corollary 7.7, it is sufficient to obtain the optimal convergence rate in the finite element method if we choose a e ( 0 , 2 3 ) for the edge e and a c ( 0 , a e ] for its two endpoints. This gives rise to the following optimal range of the grading parameters: κ c κ e = 2 - 1 a e < 2 - 3 2 0.353 near e, and quasi-uniform meshes (the associated grading parameters being 1 2 ) for all the other edges and vertices.

Figure 4 
               Example 8.1 (two Neumann faces): the top view of 
                     
                        
                           
                              Ω
                              p
                           
                        
                        
                        {\Omega_{p}}
                     
                  , top Dirichlet face marked in blue (left); the absolute value of the numerical solution (right).
Figure 4 
               Example 8.1 (two Neumann faces): the top view of 
                     
                        
                           
                              Ω
                              p
                           
                        
                        
                        {\Omega_{p}}
                     
                  , top Dirichlet face marked in blue (left); the absolute value of the numerical solution (right).
Figure 4

Example 8.1 (two Neumann faces): the top view of Ω p , top Dirichlet face marked in blue (left); the absolute value of the numerical solution (right).

In Table 1, we display the convergence rates of the finite element solution on proposed anisotropic meshes associated with different values of the grading parameter for Example 8.1. In all these meshes, we choose κ = κ e = κ c for the singular edge e and the two endpoints c (Figure 3). Here, j is the level of refinements. Denote by u j the linear finite element solution on the mesh after j refinements. Since the exact solution is not known, the convergence rate is computed using the numerical solutions for successive mesh refinements

(8.2) convergence rate = log 2 ( | u j - u j - 1 | H 1 ( Ω ) | u j + 1 - u j | H 1 ( Ω ) ) .

As j increases, the dimension of the discrete system is O ( 2 3 j ) . Therefore, the asymptotic convergence rate in (8.2) is a reasonable indicator of the actual convergence rate for the numerical solution. For example, when the numerical solution approximates the singular solution at the optimal convergence rate as described in Corollary 7.7, the convergence rate in (8.2) shall converge to 1 as the number of refinements j increases. On quasi-uniform meshes, the convergence rate (8.2) is however asymptotically bounded by 2 3 0.667 due to the fact that the solution is singular in Ω p ( u H 1 + s ( Ω p ) for s < 2 3 ).

Table 1

The H 1 convergence rates in Example 8.1.

j κ = 0.1 κ = 0.2 κ = 0.3 κ = 0.4 κ = 0.5
3 0.65 0.72 0.75 0.71 0.61
4 0.81 0.83 0.84 0.77 0.64
5 0.89 0.91 0.90 0.81 0.66
6 0.95 0.96 0.94 0.83 0.67
7 0.98 0.98 0.96 0.84 0.67
8 0.99 0.99 0.97 0.85 0.67

It is clear that the above theoretical predictions are confirmed by the numbers in Table 1. Namely, when the grading parameter to the singular edge and its two endpoints is in the optimal range

κ = κ e = κ c = 0.1 , 0.2 , 0.3 < 0.353 ,

the convergence rates in Table 1 converge to 1, which implies the optimal convergence rate in the finite element method is achieved on the anisotropic meshes proposed in Algorithm 3.4. For κ = 0.4 , 0.5 > 0.353 , the convergence is not optimal. In particular, the convergence rates for κ = 0.5 is also very close to the theoretical bound 2 3 as discussed above.

In the second example, in addition to the Neumann edge, we shall confirm the effectiveness of our numerical scheme when the domain has Neumann vertices, namely vertices surrounded by Neumann faces.

Example 8.2.

In equation (8.1), the Neumann condition is imposed on the two faces adjacent to the edge e and also on the top and bottom faces of the domain Ω p . The Dirichlet condition is imposed on all the other side faces. Therefore, e is a Neumann edge, and its two endpoints c are Neumann vertices.

Figure 5 
               Example 8.2 (four Neumann faces): the top view of 
                     
                        
                           
                              Ω
                              p
                           
                        
                        
                        {\Omega_{p}}
                     
                  , top Neumann face marked in red (left); the absolute value of the numerical solution (right).
Figure 5 
               Example 8.2 (four Neumann faces): the top view of 
                     
                        
                           
                              Ω
                              p
                           
                        
                        
                        {\Omega_{p}}
                     
                  , top Neumann face marked in red (left); the absolute value of the numerical solution (right).
Figure 5

Example 8.2 (four Neumann faces): the top view of Ω p , top Neumann face marked in red (left); the absolute value of the numerical solution (right).

According to (2.4) (see also [12]), the edge e is the only singular edge with λ e = 2 3 and the two vertices c are regular vertices. Similar to Example 8.1, all the other vertices and edges of Ω p are regular. Note that the right-hand side function in (8.1) still belongs to H σ ( Ω p ) for any σ [ 0 , 1 ) , and the solution u H 1 + s ( Ω p ) for s < 2 3 . Therefore, based on the conditions in (3.2) and (3.3), and by Corollary 7.7, it is sufficient to obtain the optimal convergence rate in the finite element method if we choose a e ( 0 , 2 3 ) for the edge e and a c ( 0 , a e ] for its two endpoints. These are the same conditions as in Example 8.1 because the error analysis in Section 7 ensures that the solution near the Neumann vertex can be approximated well by the finite element solution. Hence, the optimal range of the grading parameters are κ c κ e = 2 - 1 a e < 2 - 3 2 0.353 near e, and quasi-uniform meshes (the associated grading parameters being 1 2 ) for all the other edges and vertices.

Table 2

The H 1 convergence rates in Example 8.2.

level κ = 0.1 κ = 0.2 κ = 0.3 κ = 0.4 κ = 0.5
3 0.64 0.70 0.73 0.68 0.57
4 0.80 0.83 0.82 0.75 0.63
5 0.88 0.90 0.89 0.81 0.66
6 0.95 0.96 0.94 0.84 0.67
7 0.98 0.98 0.96 0.85 0.67
8 0.99 0.99 0.97 0.86 0.67

The numerical convergence rates (8.2) are reported in Table 2 on the anisotropic meshes associated with different values of the grading parameter. As in Example 8.1, we choose κ = κ e = κ c for the singular edge e and the two endpoints c. It is clear that the above theoretical predictions are confirmed by the numbers in Table 2. The convergence rate is optimal when κ = 0.1 , 0.2 , 0.3 < 0.353 is in the optimal range, and it slows down when κ = 0.4 , 0.5 > 0.353 . As discussed above, the anisotropic algorithm can give rise to the optimal numerical approximation for equations with mixed boundary conditions, even with Neumann edges and vertices.

In the third example for the prism domain, we test the anisotropic algorithm with the presence of a singular edge with the mixed boundary condition. Instead of equation (8.1), we shall solve another form of equation (1.1).

Example 8.3.

Consider equation (1.1) with f = 1 in the prism domain Ω p . We impose the Neumann boundary condition on one face adjacent to the edge e and impose the Dirichlet boundary condition on all the other faces (see Figure 6). Thus, e is the singular edge surrounded by faces with mixed boundary conditions.

Figure 6 
               Example 8.3 (a DN singular edge): Neumann face marked in red and Dirichlet face in blue (left); the numerical solution (right).
Figure 6 
               Example 8.3 (a DN singular edge): Neumann face marked in red and Dirichlet face in blue (left); the numerical solution (right).
Figure 6

Example 8.3 (a DN singular edge): Neumann face marked in red and Dirichlet face in blue (left); the numerical solution (right).

According to (2.4), e is the only singular edge with λ e = 1 3 . Note that the right-hand side function in (1.1) belongs to H σ ( Ω p ) for any σ [ 0 , 1 ) , and the solution u H 1 + s ( Ω p ) for s < 1 3 . Based on Algorithm 3.4, it is sufficient to obtain the optimal convergence rate in the finite element method if we choose a e ( 0 , 1 3 ) for the edge e and a c ( 0 , a e ] for its two endpoints. Hence, the optimal range of the grading parameters are: κ c κ e = 2 - 1 a e < 2 - 3 = 0.125 near e, and quasi-uniform meshes for all the other edges and vertices.

Table 3

The H 1 convergence rates in Example 8.3.

level κ = 0.1 κ = 0.2 κ = 0.3 κ = 0.4 κ = 0.5
3 0.43 0.43 0.48 0.56 0.60
4 0.72 0.73 0.72 0.70 0.68
5 0.87 0.86 0.81 0.72 0.64
6 0.93 0.90 0.80 0.66 0.54
7 0.96 0.90 0.76 0.58 0.45
8 0.97 0.89 0.70 0.52 0.39

We summarize the numerical convergence rates for Example 8.3 in Table 3, which again verifies the theoretical prediction. Namely, when the grading parameter κ = κ e = κ c < 0.125 , the optimal convergence rate is achieved, while it is not the case for κ > 0.125 .

Before we discuss the last numerical example, we define a Fichera domain as follows. Let Ω 1 = ( - 1 , 1 ) 3 and Ω 2 = ( 0 , 1 ) 3 be two cubes. Then the Fichera domain is Ω f : = Ω 1 Ω 2 . In this example, we shall show the test results for solving equation (1.1) in Ω f .

Example 8.4.

In equation (1.1), let f = 1 . We impose the Neumann boundary condition on the two faces that touch the center vertex of the Fichera domain Ω f and impose the Dirichlet boundary condition on all the other faces. See Figure 7.

Figure 7 
               Example 8.4 (three singular edges in a Fichera domain): Neumann faces marked in red and Dirichlet faces in blue (left); the numerical solution (right).
Figure 7 
               Example 8.4 (three singular edges in a Fichera domain): Neumann faces marked in red and Dirichlet faces in blue (left); the numerical solution (right).
Figure 7

Example 8.4 (three singular edges in a Fichera domain): Neumann faces marked in red and Dirichlet faces in blue (left); the numerical solution (right).

Figure 8 
               Graded meshes on the Fichera domain (left to right): the initial mesh, mesh after one refinement, mesh after three refinements (
                     
                        
                           
                              
                                 κ
                                 e
                              
                              =
                              
                                 κ
                                 c
                              
                              =
                              0.2
                           
                        
                        
                        {\kappa_{e}=\kappa_{c}=0.2}
                     
                  ).
Figure 8 
               Graded meshes on the Fichera domain (left to right): the initial mesh, mesh after one refinement, mesh after three refinements (
                     
                        
                           
                              
                                 κ
                                 e
                              
                              =
                              
                                 κ
                                 c
                              
                              =
                              0.2
                           
                        
                        
                        {\kappa_{e}=\kappa_{c}=0.2}
                     
                  ).
Figure 8 
               Graded meshes on the Fichera domain (left to right): the initial mesh, mesh after one refinement, mesh after three refinements (
                     
                        
                           
                              
                                 κ
                                 e
                              
                              =
                              
                                 κ
                                 c
                              
                              =
                              0.2
                           
                        
                        
                        {\kappa_{e}=\kappa_{c}=0.2}
                     
                  ).
Figure 8

Graded meshes on the Fichera domain (left to right): the initial mesh, mesh after one refinement, mesh after three refinements ( κ e = κ c = 0.2 ).

Therefore, according to (2.4), there are three singular edges e (the edges touching the center vertex), and the endpoints c of these singular edges are possible singular vertices. See Figure 8. For the Neumann singular edge, we have λ e = 2 3 , and for the two DN singular edges, we have λ e = 1 3 . In fact, the solution satisfies the global regularity u H 1 + s ( Ω f ) for s < 1 3 . We will use the graded mesh toward the singular edges and their endpoints, as described in Algorithm 3.4, in the finite element approximation. In particular, to simplify the presentation, we shall choose the same parameter κ e for all the singular edges, and the same parameter κ c = κ e for all their endpoints. See Figure 8 for the case κ e = κ c = 0.2 . Recall that, for any possible singular vertex c, λ c + 1 2 > 1 2 . Hence, based on (3.2) and (3.3), it is sufficient to choose κ c = κ e < 2 - 3 = 0.125 for all the singular edges and their endpoints in order to obtain the optimal convergence in the numerical approximation.

The numerical convergence rates for different values of the grading parameter in Example 8.4 are listed in Table 4. As predicted by the theory, for κ = κ c = κ e > 0.125 , the convergence is not optimal, while for κ = 0.1 < 0.125 , the numbers are increasing toward the optimal rate. We stopped at level 7 because the resources needed to extend the calculation to the next level of refinement have exceeded our computing capability. For instance, the refinement to the next level will generate more than 5 billion tetrahedra and 1 billion nodes. Nevertheless, according to Table 4, there is a clear improvement in convergence rates using the appropriate graded meshes ( κ < 0.125 ) compared with other graded meshes ( κ > 0.125 ), and it is reasonable to expect the rates for κ = 0.1 will converge to 1 when the asymptotic region is reached with further refinements.

Table 4

The H 1 convergence rates in Example 8.4.

level κ = 0.1 κ = 0.2 κ = 0.3 κ = 0.4 κ = 0.5
3 0.58 0.63 0.65 0.65 0.63
4 0.78 0.81 0.80 0.76 0.70
5 0.85 0.88 0.83 0.75 0.65
6 0.88 0.89 0.81 0.68 0.54
7 0.91 0.88 0.75 0.59 0.45

Award Identifier / Grant number: DMS-1819041

Funding statement: The first author was supported in part by the NSF Grant DMS-1819041 and by the Wayne State University Career Development Chair Award.

Acknowledgements

The authors would also like to thank Xun Lu for providing the numerical test results in this paper.

References

[1] T. Apel, Anisotropic Finite Elements: Local Estimates and Applications, Adv. Numer. Math., B. G. Teubner, Stuttgart, 1999. Suche in Google Scholar

[2] T. Apel and B. Heinrich, Mesh refinement and windowing near edges for some elliptic problem, SIAM J. Numer. Anal. 31 (1994), no. 3, 695–708. 10.1137/0731037Suche in Google Scholar

[3] T. Apel and S. Nicaise, The finite element method with anisotropic mesh grading for elliptic problems in domains with corners and edges, Math. Methods Appl. Sci. 21 (1998), no. 6, 519–549. 10.1002/(SICI)1099-1476(199804)21:6<519::AID-MMA962>3.0.CO;2-RSuche in Google Scholar

[4] T. Apel, A.-M. Sändig and J. R. Whiteman, Graded mesh refinement and error estimates for finite element solutions of elliptic boundary value problems in non-smooth domains, Math. Methods Appl. Sci. 19 (1996), no. 1, 63–85. 10.1002/(SICI)1099-1476(19960110)19:1<63::AID-MMA764>3.0.CO;2-SSuche in Google Scholar

[5] T. Apel and J. Schöberl, Multigrid methods for anisotropic edge refinement, SIAM J. Numer. Anal. 40 (2002), no. 5, 1993–2006. 10.1137/S0036142900375414Suche in Google Scholar

[6] I. Babuška and A. K. Aziz, On the angle condition in the finite element method, SIAM J. Numer. Anal. 13 (1976), no. 2, 214–226. 10.1137/0713021Suche in Google Scholar

[7] C. Bacuta, H. Li and V. Nistor, Anisotropic graded meshes and quasi-optimal rates of convergence for the FEM on polyhedral domains in 3D, European Congress on Computational Methods in Applied Sciences and Engineering—ECCOMAS 2012, ECCOMAS, 2012, 9003–9014. Suche in Google Scholar

[8] C. Bacuta, V. Nistor and L. T. Zikatanov, Improving the rate of convergence of high-order finite elements on polyhedra. II. Mesh refinements and interpolation, Numer. Funct. Anal. Optim. 28 (2007), no. 7–8, 775–824. 10.1080/01630560701493263Suche in Google Scholar

[9] S. C. Brenner and L. R. Scott, The Mathematical Theory of Finite Element Methods, 2nd ed., Texts Appl. Math. 15, Springer, New York, 2002. 10.1007/978-1-4757-3658-8Suche in Google Scholar

[10] A. Buffa, M. Costabel and M. Dauge, Anisotropic regularity results for Laplace and Maxwell operators in a polyhedron, C. R. Math. Acad. Sci. Paris 336 (2003), no. 7, 565–570. 10.1016/S1631-073X(03)00138-9Suche in Google Scholar

[11] P. G. Ciarlet, The Finite Element Method for Elliptic Problems, Stud. Math. Appl. 4, North-Holland, Amsterdam, 1978. 10.1115/1.3424474Suche in Google Scholar

[12] M. Costabel, M. Dauge and S. Nicaise, Weighted analytic regularity in polyhedra, Comput. Math. Appl. 67 (2014), no. 4, 807–817. 10.1016/j.camwa.2013.03.006Suche in Google Scholar

[13] M. Dauge, Elliptic Boundary Value Problems on Corner Domains. Smoothness and Asymptotics of Solutions, Lecture Notes in Math. 1341, Springer, Berlin, 1988. 10.1007/BFb0086682Suche in Google Scholar

[14] R. Fritzsch, Optimale Finite-Elemente-Approximationen für Funktionen mit Singularitäten, Ph.D. Thesis, TU Dresden, 1990. Suche in Google Scholar

[15] D. Gilbarg and N. S. Trudinger, Elliptic Partial Differential Equations of Second Order, Grundlehren Math. Wiss. 224, Springer, Berlin, 1977. 10.1007/978-3-642-96379-7Suche in Google Scholar

[16] P. Grisvard, Elliptic Problems in Nonsmooth Domains, Monogr. Stud. Math. 24, Pitman, Boston, 1985. Suche in Google Scholar

[17] P. Grisvard, Edge behavior of the solution of an elliptic problem, Math. Nachr. 132 (1987), 281–299. 10.1002/mana.19871320119Suche in Google Scholar

[18] V. A. Kondrat’ev, Boundary value problems for elliptic equations in domains with conical or angular points, Trudy Moskov. Mat. Obšč. 16 (1967), 209–292. Suche in Google Scholar

[19] V. A. Kozlov, V. G. Maz’ya and J. Rossmann, Elliptic Boundary Value Problems in Domains with Point Singularities, Math. Surveys Monogr. 52, American Mathematical Society, Providence, 1997. Suche in Google Scholar

[20] M. Křížek, On the maximum angle condition for linear tetrahedral elements, SIAM J. Numer. Anal. 29 (1992), no. 2, 513–520. 10.1137/0729031Suche in Google Scholar

[21] H. Li, An anisotropic finite element method on polyhedral domains: interpolation error analysis, Math. Comp. 87 (2018), no. 312, 1567–1600. 10.1090/mcom/3290Suche in Google Scholar

[22] H. Li and S. Nicaise, Regularity and a priori error analysis on anisotropic meshes of a Dirichlet problem in polyhedral domains, Numer. Math. 139 (2018), no. 1, 47–92. 10.1007/s00211-017-0936-0Suche in Google Scholar

[23] J.-L. Lions and E. Magenes, Non-homogeneous Boundary Value Problems and Applications. Vol. I, Grundlehren Math. Wiss. 181, Springer, New York, 1972. 10.1007/978-3-642-65161-8Suche in Google Scholar

[24] D. Schötzau, C. Schwab and T. P. Wihler, hp-dGFEM for second-order mixed elliptic problems in polyhedra, Math. Comp. 85 (2016), no. 299, 1051–1083. 10.1090/mcom/3062Suche in Google Scholar

Received: 2019-10-09
Accepted: 2020-01-24
Published Online: 2020-02-18
Published in Print: 2021-01-01

© 2021 Walter de Gruyter GmbH, Berlin/Boston

Heruntergeladen am 29.10.2025 von https://www.degruyterbrill.com/document/doi/10.1515/cmam-2019-0148/html
Button zum nach oben scrollen