Spectral Method and Spectral Element Method for Three Dimensional Linear Elliptic System: Analysis and Application

In this paper, a least-squares spectral method and a non-conforming least-squares spectral element method for three dimensional linear elliptic system will be presented. Differentiability estimates and the main stability theorem for the proposed method are proven. Using the regularity estimate and the proposed stability estimates, we introduce a suitable preconditioner and show that the condition number of the preconditioned system is O((lnW)2)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$O((\ln W)^2)$$\end{document} (where W is degree of polynomials). Then we establish the error estimates of the proposed method. Specific numerical results based on linear elasticity problem, fourth order problem and sixth order problem are discussed to reflect the efficiency of the proposed method.


Introduction
In applied mathematics, engineering and scientific computing, there are many problems which attract much attention to compute the numerical solution of linear elliptic system. Linear elasticity problem, fourth order elliptic problem and sixth order elliptic problem etc. are some examples. In continuum mechanics, elasticity theory is widely used in continuum models. The main use of elasticity theory is that it describes the behaviour of the solid materials after deformation by external forces. In case of relatively small deformation, the theory of linear elasticity plays an important role. In linear elasticity, we get the stress-strain relation in form of the constitutive equation. One of the example of the elastic material is isotropic homogeneous where the constitutive equation is defined in any two terms of the six moduli (i.e. bulk modulus, Young's modulus, Lamé's first parameter, shear modulus, Poisson's ratio and p-wave modulus). The other type of the elastic material is the isotropic inhomogeneous in which the inhomogeneity is imported from position-dependent Lamé's parameters [28]. Some similar models can also get in the elasticity analysis of biomolecules [34][35][36].
In last few decades, spectral method is widely used to solve partial differential equations numerically. The main feature of spectral method is high accuracy. We refer to [3,16,[29][30][31] for the analysis and various applications of spectral method/spectral element method.
Dutt et al. [8][9][10] introduced a nonconforming hp/spectral element for three dimensional elliptic problems with mixed Neumann and Dirichlet boundary conditions on non-smooth domains. This method is an exponentially accurate method for elliptic problems with mixed Neumann and Dirichlet boundary conditions on non-smooth domains. Dutt et al. [8][9][10] used a geometric mesh and the auxiliary map of the form z = ln ξ to remove the singularities at the corners and the edges. In the remaining part of the domain usual Cartesian coordinate system is used. Schötzau et al. [26,27] presented hp-DGFEM for second-order elliptic problems in polyhedra and established exponential convergence for second-order elliptic problems in polyhedra.
Recently, Khan et al. [20] presented a spectral element method for scalar valued elliptic interface problem. In this paper, we introduce a least-squares spectral method and a fully nonconforming least-squares spectral element method (LSSEM) for three dimensional linear elliptic system. The proposed theory is valid for the coupled and the decoupled systems. Thus, we can use this theory to solve elasticity problem and higher order PDE such as fourth order and sixth order etc. In this paper, we derive the regularity result for 3-D linear elliptic systems which is main key ingredient to prove the stability theorem. We also discuss a priori error estimates for the proposed method. Specifically, the proposed formulation is based on minimizing the quadratic form which consists the sum of the squares of a weighted squared norm of the residuals in the partial differential equation, the sum of the residuals in the boundary conditions in fractional Sobolev norms and enforce the continuity along the interelement boundary by adding a term which measures the sum of the squares of the jump in the function and its derivatives in fractional Sobolev norms. We refer to [6][7][8][9][10][11][12]15,[17][18][19][20][21][22] for the analysis and application of LSSEM to solve various types of partial differential equation in 2-D and 3-D.
One of the popular least-square method is first order based formulation in which the higher order PDE is reduced into first order PDE system and it is solved by least-square principle. To solve the first order based least square formulation, we need to compute the stiffness matrix and the mass matrix. The proposed least-squares formulation is free from any kind of first order reformulation. For computing the solution, the normal equation is solved by preconditioned conjugate gradient method. The optimality of condition number for the preconditioned system is also shown. Integrals which are in the numerical scheme, are computed efficiently by Gauss-quadrature rule. One of the major advantage to use the proposed approach is that there is no need to compute the stiffness matrix and the mass matrix. For more details, we refer to Sect. 4. The proposed method is nonconforming in the sense that the discrete solution space of LSSEM is not subset of H 2 (Ω) (The definition of H 2 (Ω) is given in Sect. 2). This is the reason that we enforce the continuity along the inter-element boundary by adding a term which measures the sum of the square of the jump of the function and its derivatives in appropriate norms. The feature of the nonconformity in the proposed method is the main motivation to use parallelization. An error estimate in H 1 norm (if the solution is continuous across the interface) is proven. In case of analytic data, the proposed methods achieve exponential accuracy.
The rest of this paper is organized as follows. In Sect. 2, the required function spaces and the three dimensional elliptic system problem are presented. Section 3 is devoted to discuss the discretization of the domain and the stability estimate. In Sect. 4, we formulate the numerical scheme with the help of stability estimate. The preconditioning and parallelization issues are discussed in Sect. 5. Section 6 is devoted to derive error estimates. To validate the theory, specific numerical results are presented in Sect. 7.

Preliminaries and Problem Formulation
Let Ω be an open bounded, simply-connected domain with boundary ∂Ω = Γ . Assume that the boundary Γ is smooth (C 2 -regularity). We denote vectors in bold letters such as 3 we denote the usual Sobolev space of integer order m furnished with the norm || · || m,Ω defined by, Let w be the two dimensional function defined on E, where E is either S (master square) or T (master triangle). Now we define fractional Sobolev space of order σ , where 0 < σ < 1, However, if E is S then we prefer to use the equivalent norm [23], Moreover,

Linear Elliptic System
define as a point in a space, the solution vector, external force and boundary data. Consider the elliptic system problem The coefficients a r ,s , b r , c satisfy the following assumptions: 3. b and c satisfy the following condition: The problem (2.1)-(2.2) is well posed and its solution u satisfies the following regularity estimate.

Discretization and Stability Estimate
Assume that Ω is a hexahedron in R 3 with the boundary Γ = ∪ 6 j=1 Γ j . Now, the domain Ω is divided into a set of finite number of elements consisting of curvilinear hexahedrons, tetrahedrons and prisms (say) Ω 1 , Ω 2 , . . . , Ω L .
Let {ũ l i } be the spectral element functions as the tensor product of polynomials of degree W in each variable λ 1 , λ 2 and λ 3 as follows . Let {F u i } be the spectral element representation of the function u l i i.e.
Define {F u } as the the spectral element representation of whole system i.e.
The space of spectral element functions on Q is denoted by S W {F u }.
Remark 3.1 However, we discuss the theoretical results for the general polynomial. But we use the Legendre polynomial to solve the problem numerically. Thus, the spectral element functions {ũ l i }, defined on Q, are given as where L r (·) represents the Legendre polynomial of degree r .

Stability Estimate
Define Ω l as a typical element in Ω with its faces {Γ l j } 1≤ j≤6 . Then we have where L is the differential operator L in λ = (λ 1 , λ 2 , λ 3 ) coordinates and J l (λ) denotes the Jacobian of the mapping M l (λ) from Q to Ω l .
Assume that Ω m and Ω n are two adjacent elements in the domain Ω and Γ m j be a face of Ω m and Γ n k be a face of Ω n such that Γ m j = Γ n k i.e. Γ m j is a face common to both Ω m and Ω n . Let [u]| Γ m j denote the jump in u across the face Γ m j . We assume the face Γ m j is the image of λ 3 = −1 under the mapping M m which maps Q to Ω m and also the image of λ 3 = 1 under the mapping M n which maps Q to Ω n . Then [u]| Γ m j is a function of only λ 1 and λ 2 . Moreover, Γ m j is the image of the master square S = (−1, 1) 2 under these mappings. However, the analysis remains valid if it is the image of T , the master triangle, too. By chain rule, it follows: Then we define the jump along the inter element boundaries as follows: Now along the boundary Γ = ∪ 6 j=1 Γ j , let Γ s j ⊆ Γ j (for some j) be the image of λ 2 = 1 under the mapping M m which maps Q to Ω m . Then,

Remark 3.2
To define the quadratic form of spectral method, we just define the above . Then the resulting quadratic form V W ({F u }) for spectral method is as follows: Applying the regularity result stated in Theorem 2.1 implies Combining Lemmas 3.1 and 3.2 and (3.4) leads to the stated result (3.2).

Technical Estimates
. Moreover, the following estimate holds: In addition, these node corrections satisfy the following condition: at P j over all elements which have P j as a node. Assume that r l (λ) = (r l 1 (λ), r l 2 (λ), r l 3 (λ)) have the following values on nodes: for j = 1, . . . , 8. We can compute a j , b j , c j and d j from (3.6). Then, it holds: Using (7.8) of [33] and the result in Appendix C of [14] gives represents the averages of the values of y l k (k = 1, 2, 3) at t over all elements which have I j as a side. In addition, In same manner, we can define E k, j , F k, j , G k, j and H k, j for the remaining face k = 2, . . . , 6. If faces F j ⊆ Γ j , j = 1, . . . , 6, then E k, j , F k, j , G k, j and H k, j are defined to be identically zero. Now we estimate the following term Applying Theorem 4.79 and Theorem 4.82 of [25] gives Similarly, we can prove estimates for all the terms given in Theorem 5.2 of [2]. Then, it holds: . Then, combining (3.8), (3.11) and (3.13) leads to the stated result.
Applying the trace theorem of Sobolev space implies Using Lemma 3.1 in (3.16) gives the desired result.

Numerical Scheme
In same manner, we define u = (u 1 , u 2 , u 3 ) . Then the boundary condition (u 1 , u 2 , u 3 ) = u = g on Γ . Let Γ j ∩ ∂Ω m denotes as the image of the mapping M m of Q onto Ω m corresponding to the side λ 1 = 1 and

Remark 4.1
To define the numerical scheme for spectral method, we just define the above . Then the resulting quadratic form R W ({F u }) for spectral method is as follows: Our numerical scheme may be written as follows:

Symmetric Formulation
Our method is a least-squares method. Thus, we solve the normal equation using the preconditioned conjugate gradient method (PCGM). Assume that the normal equations be Let U W and U 2W be denoted as Note that we use Gauss-Lobatto-Legendre (GLL) quadrature formula to compute the integrals given in the minimization formulation. The details regarding computation of numerical scheme is discussed in "Appendix". Finally, we obtain the following representation of the minimization formulation for each element: Then we obtain Assume that γ W l be the normalizing factors which are used to compute the discrete Legendre transform as Then is less compare to store (G − AU ). Hence, the proposed method is cheap and efficient to compute the residual in comparison of h − p finite element method. For more details, we refer to [13,32]. Specifically, the comparison between h − p fem and LSSEM for 2 − D elliptic problems is discussed in Section 5 in [32]. The similar comparison also holds for 3 − D elliptic problems.

Preconditioning and Parallelization
First, we define the quadratic form with zero data. Next we estimate the condition number of preconditioner U W ({F u }). To obtain the lower bound, we use the trace theorem for Sobolev spaces. Then where, K is a constant. Using Theorem 3.1 gives Combining (5.2) and (5.3) implies where C is a positive constant. Thus, the condition number of the preconditioned system is O((log W ) 2 ). For each element, the preconditioner corresponds to the following quadratic form where, u = u(λ) = (u 1 (λ), u 2 (λ), u 3 (λ)) . Here each component u k (λ)(k = 1, 2, 3) is a polynomial of degree W in λ 1 , λ 2 and λ 3 separately.
Using the idea of [10], we can define the new quadratic form C(u) and show that C(u) is spectrally equivalent to B(u). The new quadratic form C(u) can be diagonalized using the separation of variables. Hence, the action of the matrix corresponding to the quadratic form C(u) will be easy to compute. In algorithm, each element is mapped onto a single processor (core) for ease of parallelism. Thus, during the PCGM process, the interchange of informations based on the value of function and its derivatives at inter-element boundaries confine the communication between neighbouring processors. Note that there are two global scalers namely the approximate solution and the search direction which need to compute for update in each iterations. Thus, it is clear that the communication between inter-processor is quite small. We note that PCGM take O(W ln W ) iteration to compute the numerical solution with exponential accuracy.

) is a polynomial of degree W in each variable separately. Then for W large enough, there exists constant C s such that the estimate
Using the idea of Theorem 6.1 from [20] and Theorem 5.1 from [21], we can conclude that Combining (6.2), (6.3) and (6.4) implies Proof First, the error is divided into two parts as follows: Next we estimate the first term of R.H.S. of (6.5).

Remark 6.2
To compute the non-conforming spectral element solution, we solve the normal equations using PCGM. After computing the numerical solution, there is a possibility to get the conforming solution by a correction. Hence, we can show that the error between the conforming numerical solution and the exact solution in the H 1 norm is also exponentially small in W . We refer to [9] for more information about these corrections. Moreover, it also holds: where C and b > 0 are constants and z is the corrected solution. . In next lemma, we present the bound of error estimate in terms of degrees of freedom.

Numerical Results
Let u ap and u ex denote the spectral element solution and the exact solution, respectively. Next we define the relative error as follows: All numerical results given in this section are computed using a FORTRAN-90 code. We use a Intel(R) Xeon(R) CPU E7-8870 @ 2.40GHz based machine for our computations. More details are as follows: Number of CPU (Physically)-8, Cores per CPU (Physically)-10 and Threads per CPU -20. Recall that each element is mapped onto a single processor (core) for

Remark 7.1
It is well known that singularities based on corners, edges and edge-corner arise in 3-D cubic domain. But, in our computations, we take our data in a way that the solution does not have any singularities.
The computed results based on spectral method and spectral element method are presented in Tables 1 and 2. Tables 1 and 2 represent the relative error in H 1 −norm with the number of iterations and CPU time against W . From Tables 1, 2 and Fig. 2, it is clear that the error decays exponentially with respect to W . In Table 3, the number of iterations and CPU time with preconditioner and with out preconditioner are provided with respect to W . It shows the effectiveness of the proposed preconditioner in terms of time and iteration count. In Fig. 3, the graph of iteration count grows with O(W ln W ). Example 2 (Biharmonic problem with simply supported boundary condition) Consider the problem Let v = −Δu, then The exact solution u is given by: u = cos(π x) sin(π y) sin(π z).
Here, we first convert the fourth order problem to second order elliptic system. Then, we solve the problem. Tables 4 and 5 display the computed results based on spectral method and spectral element method. Specifically, Tables 4 and 5 represent the relative error in H 1 −norm with the number of iterations and CPU time against W . Note that, in Tables 4, 5 and Fig. 4, the error decays exponentially with respect to W . It is easy to see that the convergence rates of u and v = Δu are similar. In Fig. 5, the graph of iteration count grows with O(W ln W ).
Example 3 (Linear elasticity problem with non-homogeneous boundary condition) Consider the following three dimensional steady state linear elasticity problem where T = λtr(σ )I + 2μσ. (7.1) Here (λ, μ), I and σ are respectively the Lamé's parameters, 3 × 3 identity matrix and stress tensor. The stress tensor and the Lamé's parameters can also be written as follows: where E and ν are Young modulus and Poisson ratio, respectively ( Table 6). The exact solution u = (u 1 , u 2 , u 3 ) is given by: In this example, we fix E = 1 and the values of ν are as follows: Here, we are presenting the numerical results of linear elasticity equation with three different choices of Poisson ratio. Numerical results based on spectral method and spectral element method for all three cases are presented in Tables 7 and 8. From Tables 7, 8 and Fig. 6,  it is clear that the error decays exponentially with respect to W . Moreover, the convergence rate of the error in all three choices are similar.
Let v = Δ 2 u, and w = Δu then The exact solution u is given by: We first introduce two new variables v = Δ 2 u and w = Δu, and convert six order problem into second order elliptic system. Numerical results based on spectral method and spectral element method are presented in Tables 9 and 10. From Tables 9, 10 and Fig. 7, it is clear that the error decays exponentially with respect to W . To convert the six order elliptic equation to second order elliptic system, we obtain the same convergence order for each components. Hence, higher order derivative Δ 2 u = v and Δu = w also obtain same convergence order as the convergence order of u.  where A = exp(x + y + z)I and I is the 3 × 3 identity matrix. The exact solution u = (u 1 , u 2 , u 3 ) is given by:

Example 5 Consider the problem
In Table 11, we present the computed results based on spectral method and spectral element method. Specifically, Table 11 represents the relative error in H 1 −norm. It is clear from Table 11 and Fig. 8 that the error decays exponentially with respect to W . Thus, it confirms that the proposed theory is also valid for the problem with variable coefficients.

Remark 7.2
In this remark, we present the numerical results for the problems discussed in Examples 3 and 5. Specifically, we solve the problems in the complex domain (unit sphere) Ω = {(x, y, z) ∈ R 3 : x 2 + y 2 + z 2 < 1}. We divided our domain into 7 curvilinear hexahedrons, for more details see also [20,Example 7.3]. We present the convergence results for Example 3 in Fig. 9a and for Example 5 in Fig. 9b using spectral element method. Figure 9a, b reflect that the error decays exponentially in W .     The exact solution u = (u 1 , u 2 , u 3 ) is given by: where α ≥ 2 is a positive real number.
In this example, we want to show the convergence rate of the proposed methods for the solutions that are not smooth enough. The exact solution u has two properties. First, It is smooth for the positive integer α. Second, it is not analytic for the positive non-integer α. We denoteα by the greatest integer less than or equal to α. If α is not integer, then ∂α x , ∂α y and ∂α z have singularities near to zero. In Tables 12 and 13, we can easily see that the error is equal to machine accuracy for the polynomial order W ≥ 3 and α = 3. In case of α = 3.5, the error decays with the low convergence rate. We plot the convergence rate for the different choices of α in Fig. 10a, b. The convergence rate of the error increases from α = 2.5 to α = 6.5 in Fig 10a, b.

Conclusions and Future Work
In this article, a least-squares spectral method and a fully non-conforming least-squares spectral element method are studied for three dimensional elliptic system problems. In both cases (spectral method and spectral element method), it is shown that the error in the computed solution decays exponentially in polynomial degree W . Computational results for specific test problems confirm the estimates obtained for the error and computational complexity. Our algorithm is quite simple and easy to implement on parallel computers. We plan to develop numerical schemes for three dimensional elliptic system with corner singularity and edge singularity in future work.

Thus, we have
Qṽ Similarly, we can write the expression for the remaining terms of R term . To arrangeũ = (ũ 1 ,ũ 2 ,ũ 3 ) in lexicographic order, we define In same manner, we define Thus, the final form is as follows: Here each R k i (i, k = 1, 2, 3) is a matrix such that R k i Z 2W i (i, k = 1, 2, 3) is easily computed.
Similarly, we have Now we compute the boundary terms in norm H 1/2 (S). Letl be the polynomial of degree 2W in λ 1 and λ 2 separately.