Bizarre Love Triangle | 31.01.1823:59 Uhr Tina York - Ich Darf Das! (Remastered) Pop-Schlager320 kbit/s 0 / 04.067 Hits VID P2P DDL 0 Kommentare | One Week Friends

spectral element method for elliptic problems

spectral element method for elliptic problems

Applied Numerical Mathematics 60 (2010) 38–54 Contents lists available at ScienceDirect Applied Numerical Mathematics www.elsevier.com/locate/apnum ...

416KB Sizes 3 Downloads 42 Views

Applied Numerical Mathematics 60 (2010) 38–54

Contents lists available at ScienceDirect

Applied Numerical Mathematics www.elsevier.com/locate/apnum

Least-squares hp/spectral element method for elliptic problems N. Kishore Kumar a,∗ , G. Naga Raju b a b

Max Planck Institute for Mathematics in the Sciences, Insel Strasse, 22, Leipzig, Germany Department of Mathematics and Statistics, Indian Institute of Technology, Kanpur, India

a r t i c l e

i n f o

Article history: Received 20 April 2009 Received in revised form 29 August 2009 Accepted 31 August 2009 Available online 8 September 2009 MSC: 65M70 65N35 65Y05 74B05

a b s t r a c t The solution of elliptic boundary value problems often leads to singularities due to nonsmoothness of the domains on which the problem is posed. This paper studies the performance of the nonconforming hp/spectral element method for elliptic problems on nonsmooth domains. This paper deals with monotone singularities of type r α and r α logδ r as well as the oscillating singularities of type r α sin(ε log r ). © 2009 IMACS. Published by Elsevier B.V. All rights reserved.

Keywords: Geometric mesh Least-squares solution Preconditioner Auxiliary mapping Exponential accuracy

1. Introduction In [1,2] Babuska and Guo proposed an exponentially accurate method in the frame work of hp finite element method to deal the singularities in the solution. They were able to resolve the singularities which arise at the corners by using a geometric mesh. In [3,14] Babuska and H.S. Oh have introduced the method of auxiliary mapping (MAM). With this method exponential rates of convergence was recovered for Laplace equation with corner singularities, in the context of p version of finite element method. In [13] Lucas and Oh extended this method for Helmholtz equations. The method of auxiliary mapping (MAM) introduced by Babuska and Oh in [3] was proven to be successful in dealing with r α singularities. However, the effectiveness of MAM is reduced in handling r α sin( log r ) type singularities. In [15] H.S. Oh et al., introduced the power auxiliary mapping (PAM) and the exponential auxiliary mapping (EAM) and shown that the method is highly accurate in dealing the singularities of type r α , r α logδ r and r α sin( log r ), where 0 < α < 1. They presented numerical results for various test problems and compared the results with the results obtained by the hp finite element method. In the latest book by Pavel and Gunzburger [4] the least-squares finite element method (LSFEM) for elliptic problems have been summarized. The standard techniques used to convert the second order elliptic equations into first order system, for example like div/curl systems shows poor rates of convergence in the presence of singularities in the solution. A weighted norm first order system least-squares (FOSLS) for problems with corner singularities have been proposed in [5,11]. The

*

Corresponding author. E-mail address: [email protected] (N.K. Kumar).

0168-9274/$30.00 © 2009 IMACS. Published by Elsevier B.V. All rights reserved. doi:10.1016/j.apnum.2009.08.008

N.K. Kumar, G.N. Raju / Applied Numerical Mathematics 60 (2010) 38–54

39

method recovers optimal order accuracy in the weighted L 2 and H 1 norms and retains optimal L 2 convergence near the singularities. The error estimates shows only an algebraic convergence. The elliptic problem on nonsmooth domains have been examined by Pathria and Karniadakis in [16] and Karniadakis and Spencer in [8], in the framework of spectral/hp element methods. In [19,20] hp version of the mortar finite element method (nonconforming) have been studied. The method is optimal for both h and p version and exponential accuracy is obtained with hp version. Mathematical frame work for different discontinuous Galerkin methods for elliptic problems have been discussed in [8]. In [23] spectral/hp element methods for solving elliptic boundary value problems on polygonal domains using parallel computers were proposed. For problems with Dirichlet boundary conditions the spectral element functions were nonconforming. For problems with Neumann and mixed boundary conditions the spectral element functions have to be continuous at the vertices of the elements only and nonconforming otherwise. In [7] P. Dutt et al. proposed an exponentially accurate nonconforming hp/spectral element method to solve general elliptic boundary value problems with mixed Neumann and Dirichlet boundary conditions on nonsmooth domains. The stability and error estimates have been proved. In this paper we briefly describe the method and present the energy norm performance of elliptic boundary value problems containing the singularities of the type r α , r α logδ r and r α sin( log r ). A geometric mesh is used in the neighborhood of the corners and the auxiliary map of the form z = ln ξ is introduced to remove the singularities at the corners, which was first introduced by Kondratiev in [9]. In the remaining part of the domain usual Cartesian coordinate system is used. The spectral element functions are nonconforming. The method is essentially a least-squares method and the solution can be obtained by solving the normal equations using the preconditioned conjugate gradient method (PCGM) without computing the mass and stiffness matrices [7,10,23]. A novel preconditioner is proposed for the method which is a block diagonal matrix, where each diagonal block corresponds to an element [6]. The condition number of the preconditioner is O (ln W )2 , where W is the degree of the approximating polynomial. Let N denote the number of layers in the geometric mesh such that W is proportional to N. Then the method requires O ( W ln W ) iterations of the PCGM to obtain the solution to exponential accuracy. The contents of this paper are organized as follows: In Section 2 the problem is stated and the numerical scheme is described briefly. In Section 3 computational results are provided for various test problems. Appendix A contains details of discretization of the domain and the numerical formulation in brief. 2. Numerical scheme The numerical method is briefly described in this section. The complete details of the stability estimate, error estimates and numerical scheme were given in [7]. The stability estimate has been proved for strongly elliptic operator which satisfy the Babuska–Brezzi inf-sup condition on curvilinear polygons whose sides are piecewise analytic. 2.1. Elliptic equation on a polygonal domain Let Ω be a polygonal domain in R2 with boundary ∂Ω = Γ as shown in Fig. 1. Let the vertices of Ω be given by E 1 , E 2 , . . . , E p and the corresponding sides by the segments Γ1 , Γ2 , . . . , Γ p , where Γi joins the points E i −1 and E i . Let the   angle subtended at E j be ω j . Further, let Γ = Γ [0] ∪ Γ [1] , Γ [0] = i ∈D Γ¯i , Γ [1] = i ∈N Γ¯i where D is a subset of the set {i | i = 1, . . . , p } and N = {i | i = 1, . . . , p } \ D .

Fig. 1. Polygonal domain.

40

N.K. Kumar, G.N. Raju / Applied Numerical Mathematics 60 (2010) 38–54

Denote by H m (Ω) the Sobolev space of functions with square integrable derivatives of integer order  m on Ω furnished with the norm    D α u 2 u 2 = . H m (Ω)

L 2 (Ω)

|α |m

Further, let

 



u 2s, I =

u 2 (x) dx + I

I

I

|u (x) − u (x )|2 dx dx |x − x |1+2s

denote the fractional Sobolev norm of order s, where 0 < s < 1. Here I denotes an interval contained in R. Consider a two dimensional elliptic boundary value problem

Lu = f u=g



∂u ∂N

in Ω,

[0]

on Γ [0] ,



= g [1] on Γ [1] . A

2

Here L = −



i , j =1 ∂ xi

(ai j (x) ∂∂x j ) +

(2.1)

2



i =1 b i (x) ∂ xi

+ c (x) is a strongly elliptic operator which satisfy the Babuska–Brezzi inf-

¯ , and ( ∂ u ) A denotes the usual conormal sup condition on Ω and the coefficients ai , j (x) = a j ,i (x), b i (x), c (x) are analytic on Ω ∂N derivative which is defined as follows. Let N = ( N 1 , N 2 ) denote the outward normal to the curve Γi for i ∈ N . Then



∂u ∂N

 (x) = A

2 

N r ar , s

r , s =1

∂u . ∂ xs

(2.2)

¯ and g [l] , l = 0, 1 is analytic on every closed arc Γ¯i and g [0] is continuous Assume that the given data f is analytic on Ω on Γ [0] . 2.2. Discretization and local transformation Discretize the polygonal domain Ω into p nonoverlapping polygonal subdomains S 1 , S 2 , . . . , S p , where S k denotes a subdomain which contains the vertex E k only. Let S k = {Ωik, j : j = 1, 2, . . . , J k , i = 1, 2, . . . , I k } be a partition of S k , where J k and I k are integers (Fig. 1). I k is bounded for all k. Let (rk , θk ) denote polar coordinates with center at E k . Choose ρ so that the sector Ω k with sides Γk and Γk+1 bounded by the circular arc B kρ centered at E k with radius ρ , is such that

Ω k ⊆ S k . Then Ω k can be represented as

Ω k = (x1 , x2 ) ∈ Ω : 0 < rk < ρ .

Let {ψik }i =1,..., I k +1 be an increasing sequence of points such that ψ1k = ψlk and ψ Ik +1 = ψuk . Let ψik = ψik+1 − ψik . Choose k these points so that







max max ψik  λ min min ψik k

i

k



i

for some constant λ. Now choose a geometric mesh with N layers in Ω k with a geometric ratio qk (0 < qk < 1). Let j  N + 1 and Let

σ jk = ρ (qk )N +1− j for 2 

σ = 0. k 1

Ωik, j = (x1 , x2 ): σ jk < rk < σ jk+1 , ψik < θk < ψik+1 , for 1  i  I k , 1  j  N. In the remaining part of S k , for 1  k  p, we retain the Cartesian coordinate system (x1 , x2 ) i.e., in Ωik, j for 1  i  I k , N < j  Jk. Let

Ω p +1 = Ωik, j : 1  i  I k , N < j  J k , 1  k  p . Relabel the elements of Ω p +1 and write

p +1 Ω p +1 = Ωl , 1  l  L , where L denotes the cardinality of Ω p +1 .

N.K. Kumar, G.N. Raju / Applied Numerical Mathematics 60 (2010) 38–54

Fig. 2. Quasi-uniform mesh in

41

τk and θk coordinates.

τk = ln rk in the sector Ω k for 1  k  p. Define ζ kj = ln σ jk for 1  j  N + 1. Here ζ1k = −∞. Define

Ω˜ ik, j = (τk , θk ): ζ kj < τk < ζ kj+1 , ψik < θk < ψik+1 ,

Now let

for 1  i  I k , 1  j  N. Hence the geometric mesh Ωik, j , 2  j  N becomes a quasi-uniform mesh in modified polar

˜ k is a semi-infinite strip. coordinates (Fig. 2). However, Ω i ,1 2.3. Approximation

The nonconforming spectral element functions are sum of tensor products of polynomials of degree W j , 1  W j  W ˜ k for 1  k  p, 1  i  I k , 2  j  N. In the infinite sector i.e., in in their respective modified polar coordinates (A.1) in Ω i, j

Ω˜ ik,1 , the solution is approximated by a constant which is the value of the function u at the corresponding vertex E k . The

constant value is computed by treating it as a common boundary value during the numerical computation. The quadrilateral elements of Ω p +1 are mapped onto the square S = (−1, 1) × (−1, 1) and the element function is represented as a sum of tensor product of polynomials of degree W in ξ and η , the transformed variables (A.2). 2.4. The numerical formulation

We seek a solution which minimizes the sum of the squares of a weighted squared norm of the residuals in the partial differential equation and the sum of the squares of the residuals in the boundary conditions in fractional Sobolev norms and enforce continuity by adding a term which measures the sum of the squares of the jump in the function and its derivatives in fractional Sobolev norms (see (A.5)). The method is essentially a least-squares method and the normal equations can be solved using the preconditioned conjugate gradient method (PCGM). Let the normal equations be

AU = h.

(2.3)

The vector U composed of the values of the spectral element functions at Gauss–Legendre–Lobatto points is divided into two sub vectors one consisting of the values of the spectral element functions at the vertices of the domain constitute the set of common boundary values U B (corresponding to the constant approximation in the semi-infinite strip) and the other

U  consisting of the remaining values which we denote by U I . Now corresponding to the decomposition of U = U I , A and h B

has the forms



A=

AII

AIB

ABI

ABB



 and h =

hI hB

 (2.4)

.

To solve the matrix Eq. (2.3) we use the block L-U factorization of A, viz.

 A=

I

0

1 A TIB A − II

I



AII

0

0

S



I

1 A− II AIB

0

I

where the Schur complement S is defined as 1 S = A B B − A TIB A − II AIB.

 ,

(2.5)

42

N.K. Kumar, G.N. Raju / Applied Numerical Mathematics 60 (2010) 38–54

To solve Eq. (2.3) based on the L-U factorization of A as given in (2.5) reduces to first solving the system of equations

SU B = h˜ B ,

(2.6)

1 h˜ B = h B − A TIB A − I I hI .

(2.7)

where

Once solved for U B using (2.6), U I can be obtained by solving

A I I U I = hI − A I B U B . The feasibility of such a process depends on the ability to compute A I B U B , A I I U I and A B B U B for any U I , U B efficiently and this can always be done if A V can be computed inexpensively for any vector V . It has been shown in [23] that A V can be computed inexpensively without computing the mass and stiffness matrices. However in addition to this it is imperative that we should be able to construct effective preconditioner for the matrix A so that the condition number of the preconditioned system is as small as possible. If this can be done then it will be possible to compute A −1 V efficiently using the PCGM for any vector V . It has been shown in [6] that a block diagonal matrix can be constructed as a preconditioner matrix for the matrix in the normal equations, where each diagonal block corresponds to a particular element which is mapped onto the master square S. The condition number of the preconditioned system is O ((ln W )2 ). Hence to compute A −1 V to an accuracy of O (e −bW ) would require O ( W ln W ) iterations of the PCGM. The solution vector U is obtained to an accuracy of O (e −bW ) using O ( W ln W ) iterations of the PCGM. After obtaining the nonconforming solution at the Gauss–Legendre–Lobatto points, a set of corrections are performed so that the solution is conforming and belongs to H 1 (Ω). These corrections are similar to Lemma 4.57 of [18]. Then for W large enough the error estimate

u ex − z1,Ω  C e −bW

(2.8)

holds, where C and b are constants and z is the corrected solution. The proof for this estimate follows similar to the proof of Theorem 3.1 of [22]. 3. Numerical results To show the effectiveness of the method we considered the Laplace and the Poisson equations for which the exact solution is in one of the singular form r α , r α logδ r and r α sin( log r ). Further to check the performance of the method we have presented the results for homogeneous Helmholtz equation and for the Motz problem. Elasticity equations are also considered as an another example for showing the applicability of the method for system of equations. Taken W j = W for all j (A.1) and the number of layers N in the geometric mesh to be equal to W . The computations are carried out on a single processor. As explained in the previous section, after obtaining the nonconforming solution a set of correction are performed so that the solution belongs to H 1 . In this section we mean by Iters the total number of iterations required to compute the Schur complement matrix, solve for the common boundary value problems and finally to obtain the solution. e  E The relative error e  E R is defined as e  E R = , where . E denotes energy norm (H 1 -norm). Since the total number u  E of degrees of freedom (DOF) M is proportional to W 3 the relative error e  E R in the energy norm satisfies the estimate

e  E R  C e −bM

1/3

(3.1)

.

Example 1. Consider the Laplace equation on a circular domain of radius 1 as shown in Fig. 3(a) with Dirichlet boundary conditions on crack panels O A, O C and on the boundary of the circular region O AC . Let (r , θ) denote the polar coordinates

Fig. 3. The scheme of cracked domain.

N.K. Kumar, G.N. Raju / Applied Numerical Mathematics 60 (2010) 38–54

43

with origin at the vertex O = (0, 0). Let us choose the data so that the solution u has the form of the leading singularity 1

r 2 sin( 12 θ) at the origin. Considering the symmetry in the domain, the problem is solved on semi-circular disk as shown in Fig. 3(b). A geometric mesh is used in the neighborhood of O , with a geometric ratio q = 0.15. Choose λ = 1/5 (in (A.3)). In Fig. 4 the sequence of meshes for the domain Fig. 3(b) are shown, with geometric mesh refinement.

Fig. 4. Geometric mesh with different mesh refinements in θ direction.

The values of relative error e  E R (in percentage) are reported for different values of W in Table 1. The relative error is almost same for all the three different meshes. It is noted that the DOF is less for the discretization shown in Mesh 3 in comparison with Mesh 2 and Mesh 1. The DOF can be further reduced by using variable degree of polynomial approximation (here a uniform degree of polynomial approximation is used for computational simplicity). Table 1 Relative error in percent against W for different meshes. W

Mesh 1

Mesh 2

Mesh 3

2 3 4 5 6 7 8 9

10.0939 4.00533 1.55574 0.60104 0.23338 0.09020 0.03500 0.01353

10.100974 4.005878 1.556295 0.600855 0.233532 0.090178 0.034998 0.013532

10.113946 4.0100051 1.5560239 0.6000372 0.2333148 0.0900233 0.0349554 0.0135070

Table 2 contains the iteration count, number of degrees of freedom, the relative error in percentage and the constants b, C in (3.1) for the Mesh 3, for different values of the polynomial order W . Table 2 Values of Iterations, M 1/3 , relative error (%), b and C against W . W

Iters

M

e  E R %

b

C

2 3 4 5 6 7 8 9

22 40 49 66 76 87 101 115

31 89 191 349 575 881 1279 1781

10.113946 4.0100051 1.5560239 0.6000372 0.2333148 0.0900233 0.0349554 0.0135070

0.741966 0.745101 0.747182 0.747927 0.748825 0.749115 0.749944 0.749972

1.090843 1.130016 1.157822 1.168322 1.181549 1.186005 1.199242 1.199703

In Fig. 5 a graph is plotted for log e  E R against M 1/3 , for Mesh 3. The graph is shown to be a straight line which shows the exponential rate of convergence as it obeys exactly the error estimate (3.1). Example 2. Consider the Laplace equation on a domain O A B, as shown in Fig. 6. Neumann boundary conditions are taken on sides O A, O B and Dirichlet boundary condition on the rest of the boundary. Let r and θ denote the polar coordinates with origin at the vertex O . Chosen the data such that the exact solution u has the form r 2/3 cos( 23 θ) + r 4/3 cos( 43 θ). Mesh 1 and Mesh 2 in Fig. 6 are two different discretization of the domain with a geometric ratio q = 0.15. Let us choose λ = 1/4. In Table 3 the values of degrees of freedom ( M ), relative error in percentage and constant values b, C for Mesh 1 and Mesh 2 are tabulated against W . Since the size of the elements is large in Mesh 2, the difference in the relative errors can be noted for smaller values of W , compared to the relative error for Mesh 1. But as W increases the error is same for both Mesh 1 and Mesh 2. The number of degrees of freedom is less for Mesh 2 and it can be further reduced by considering variable polynomial approximation. Fig. 7 shows the graph in M 1/3 × log e  E R scale. The graph is a straight line for Mesh 1 and for Mesh 2 it has some jumps initially and as W increases it becomes a straight line.

44

N.K. Kumar, G.N. Raju / Applied Numerical Mathematics 60 (2010) 38–54

Fig. 5. Log of relative error against M 1/3 .

Fig. 6. Sectoral domain.

Table 3 Computational values for Mesh 1 and Mesh 2. Mesh 1

Mesh 2

W

M

e  E R %

M

e  E R %

b

C

Iters

2 3 4 5 6 7 8 9

46 133 286 523 862 1321 1918 2671

22.7853 3.1268 0.5927 0.1055 0.0265 0.0073 0.0020 0.0005

31 89 191 349 575 881 1279 1781

55.4689 3.5228 2.1518 0.1151 0.0494 0.0073 0.0021 0.0005

1.19044 1.14055 1.14816 1.05222 1.07138 1.00165 1.00501 1.00027

10.994 6.2712 6.8545 2.1441 2.7263 1.1041 1.1548 1.0816

20 38 52 71 87 110 127 145

Example 3. Consider a homogeneous Helmholtz equation on a domain Ω as shown in Fig. 8(a),

−u + u = 0 in Ω, ∂u = 0 on Γ1 , ∂n  sinh √ (r ) cos(θ/2) on Γ2 , r u= 0 on Γ3 . Discretization of the domain, with a geometric ratio q = 0.15 is shown in Fig. 8(b). The problem has a singularity at sinh(r ) O = (0, 0) and the exact solution of the problem is √ cos(θ/2). r

N.K. Kumar, G.N. Raju / Applied Numerical Mathematics 60 (2010) 38–54

45

Fig. 7. Log of relative error against M 1/3 .

Fig. 8. (a) The domain. (b) Discretization.

Table 4 Values of Iteration count, DOF ( M ), the relative error (%), b and C against W . W

Iters

M

e  E R %

2 3 4 5 6 7 8 9

18 43 55 74 89 102 114 127

31 89 191 349 575 881 1279 1781

16.923085 6.76869 2.037648 0.661117 0.235673 0.088490 0.033923 0.013098

b

C 0.793474 0.790778 0.774419 0.761727 0.754788 0.751745 0.750621 0.750326

1.9648 1.9061 1.5744 1.3500 1.2375 1.1896 1.1719 1.1671

Let us choose λ = 1/5. In Table 4 the values of iteration count, degrees of freedom, relative error (in percentage) and constants b, C for different values of W are tabulated. Fig. 9 shows the graph in M 1/3 × log e  E R scale. The graph is a straight line and this shows the exponential rate of convergence. Example 4 (Motz problem). Consider the Laplace’s equation −u = 0 in a rectangular domain Ω1 = {(x, y ) | −1 < x < 1, 0 < y < 1} as shown in Fig. 10, satisfying the following boundary conditions

u |x<0, y =0 = 0,

u |x=1 = 500,

u y | y =1 = u y | y =0,x>0 = u x |x=−1 = 0. This problem has been widely studied by many researchers. Rosser and Papamichael in 1975 [17] succeeded in finding the closed form solution for this problem and represented it in the following form

46

N.K. Kumar, G.N. Raju / Applied Numerical Mathematics 60 (2010) 38–54

Fig. 9. Log of relative error against M 1/3 .

u=

P 

Fig. 10. The domain Ω1 .



 1 θ bl r l+1/2 cos l +

(3.2)

2

l =0

for P = 19. Here (r , θ) denotes the polar coordinates. For the convenience of the readers the coefficients bl , l = 0, 2, . . . , 19 are provided in Appendix A3. More accurate solution for this problem has also been given in [12], that is for P = 33. Here we have considered the Motz problem on a domain Ω which is a semi-circle of radius 1 as shown in Fig. 8(a). Chosen (3.2) as the exact solution and restricted the expansion of u on Γ2 i.e., applied Dirichlet boundary condition on semi-circular arc. Further, chosen the following conditions on the other part of the boundaries

u |Γ3 = 0,

 ∂ u  = 0. ∂ n Γ1

The discretization of the domain, with a geometric ratio q = 0.15 is as shown in Fig. 8(b). Choose λ = 1/5. In Table 5 the values of iteration count, relative error (in percentage) for different values of W are tabulated. Table 5 Values of Iteration count, relative error (%) against W . W

Iters

M

e  E R %

2 3 4 5 6 7 8 9

29 41 60 76 92 107 119 135

31 89 191 349 575 881 1279 1781

17.4791 4.7940 1.8151 0.6587 0.2618 0.0969 0.0375 0.0145

N.K. Kumar, G.N. Raju / Applied Numerical Mathematics 60 (2010) 38–54

47

Fig. 11. Log of relative error against M 1/3 .

Fig. 11 shows the graph in M 1/3 × log e  E R scale which is a straight line and this indicates the exponential rate of convergence of the method. Example 5. Consider the Poisson equation

− u = f in Ω = (r , θ): r  r0 , 0  θ  π . In Fig. 12, the discretization of the domain Ω is shown. Chosen the data such that exact solution u is of the form r 0.5 log2 r cos θ . Dirichlet boundary condition is imposed on r = r0 and Neumann boundary condition is imposed on the remaining part of the boundary.

Fig. 12. Discretization of the domain.

As mentioned in [15], for the geometric mesh refinement the ratio q = e −1.5π gives the better results (but not optimal), here too a geometric ratio q = e −1.5π is taken for the geometric mesh and r0 = 2. Auxiliary map is used over the whole domain for simplicity. Choose λ = 1/5. Table 6 contains the values of the iterations, degrees of freedom and the relative error in percentage against the degree of polynomial approximation W . Fig. 13 shows the log of relative error against degrees of freedom. Results shows that our approach gives a better results than the hp finite element method and yields the same rates of convergence as EAM comparing with the values given in [15]. Example 6. Consider the Poisson equation on the domain as shown in Fig. 12 with r0 = 2. Choosing f such that the exact solution has the oscillatory singularity of the form r α sin( log r ) cos θ with respect to various sizes of oscillating factor  . As the oscillating factor becomes smaller, the singular function is less oscillating. Here we consider the weak singular function with  = 0.1 and next the highly oscillating function with  = 3.0.

48

N.K. Kumar, G.N. Raju / Applied Numerical Mathematics 60 (2010) 38–54

Table 6 Values of Iterations, DOF and relative error (%) against W . W

Iters

M

e  E R %

2 3 4 5 6 7 8 9

25 44 64 78 94 103 115 123

31 89 191 349 575 881 1279 1781

37.566 10.212 3.726 0.765 0.106 0.012 0.001 0.0001

Fig. 13. Log of relative error against M 1/3 .

Consider the case with  = 0.1 and a geometric mesh with a geometric ratio q = e −1.5π . Auxiliary map is used over the whole domain for simplicity. Choose λ = 1/5. The computational values such as iterations, degrees of freedom and relative error in percentage are provided in Table 7 for different values of W . In Fig. 14 the log of relative error against degrees of freedom is drawn. Table 7 Values of Iterations, DOF and relative error (%) against W . W

Iters

M

e  E R %

2 3 4 5 6 7 8 9

17 38 47 62 70 86 103 111

31 89 191 349 575 881 1279 1781

23.1583 4.59133 1.25379 0.19646 0.02104 0.00183 0.00013 0.00001

Let u = r 0.5 sin( log r ) cos θ with  = 3.0, it is a highly oscillatory singular function. Consider r0 = 2 and geometric mesh with geometric ratio q = 0.15. Choose λ = 1/5. The geometric ratio q = 0.15 gives better results among the four ratios q = 0.15, q = e −π , q = e −1.5π and q = e −2π . Table 8 contains the values of iteration count, degrees of freedom M and relative error in percentage against W . Fig. 15 shows graph of the log of relative error in percentage against degrees of freedom. Results suggests that our approach gives a better results than the hp finite element method and MAM by comparing the results shown in [15].

N.K. Kumar, G.N. Raju / Applied Numerical Mathematics 60 (2010) 38–54

Fig. 14. Log of relative error against M 1/3 .

Table 8 Values of Iterations, DOF and relative error (%) against W . W

Iters

M

e  E R %

2 3 4 5 6 7 8 9

22 36 52 64 77 85 94 100

31 89 191 349 575 881 1279 1781

78.166 54.840 12.04 6.464 0.68 0.31 0.036 0.010

Fig. 15. Log of relative error against M 1/3 .

49

50

N.K. Kumar, G.N. Raju / Applied Numerical Mathematics 60 (2010) 38–54

Example 7. Let u = (u 1 , u 2 ) T be a displacement vector. Consider the following plane strain linear elasticity problem on the domain with a re-entrant crack as shown in Fig. 16(a) when the body forces are not present,

     ∂ ∂ ∂ u1 ∂ u2 ∂ u1 ∂ u2 − c 11 + c 12 c 66 + − = 0, ∂ x1 ∂ x1 ∂ x2 ∂ x2 ∂ x2 ∂ x1      ∂ ∂ u1 ∂ u2 ∂ ∂ u1 ∂ u2 c 66 + c 12 + c 22 − − = 0. ∂ x1 ∂ x2 ∂ x1 ∂ x2 ∂ x1 ∂ x2

(3.3)

Fig. 16. (a) Crack domain. (b) Discretization.

Let E denote the modulus of elasticity and ν (0 < ν < 0.5) denote the Poisson’s ratio, then coefficients c 11 , c 12 , c 22 and E (1−ν ) Eν ν . c 66 are given by c 11 = c 22 = (1+ν )(1−2ν ) , c 12 = (1+ν )( , and c 66 = 2(1E+ 1−2ν ) ν) On the crack panels i.e., on the boundary Γ1 and Γ2 traction boundary conditions T u = ( T 1 u , T 2 u ) T are taken and Dirichlet boundary conditions are taken on the rest of the boundaries Γi , i = 3, . . . , 7. Let n = (n1 , n2 ) be the unit outward normal on the boundary, then the traction components T 1 u and T 2 u are given by

   ∂ u1 ∂ u2 ∂ u1 ∂ u2 + c 12 n1 + c 66 + n2 = g 1 , ∂ x1 ∂ x2 ∂ x2 ∂ x1     ∂ u1 ∂ u2 ∂ u1 ∂ u2 T 2 u = c 66 + n1 + c 12 + c 22 n2 = g 2 . ∂ x2 ∂ x1 ∂ x1 ∂ x2 

T 1 u = c 11

Let (r , θ) denote the polar coordinates with origin at E 1 . The data chosen such that the solution is of the form of the leading singularity at E 1 . Thus we choose the Mode 1 displacement components [21]

u1 =

1











κ − Q (α + 1) cos α θ − α cos (α − 2)θ ,

2c 66    1 α  u2 = r κ + Q (α + 1) sin α θ + α sin (α − 2)θ , 2c 66

where κ = 3 − 4ν , α = 0.5 and Q = 0.333. Further, we choose E = 1, ν = 0.3 and the geometric ratio q = 0.15. Let us choose λ = 1/7. Table 9 shows the values of iteration count, relative error e  E R in percentage against W . Table 9 Values of Iterations and relative error (%) against W . W

Iters

M

e  E R %

2 3 4 5 6 7 8

66 234 406 577 767 1086 1678

386 962 1922 3362 5378 8066 11522

0.1331E+02 0.3720E+01 0.1050E+01 0.3916E+00 0.1439E+00 0.5005E−01 0.1954E−01

N.K. Kumar, G.N. Raju / Applied Numerical Mathematics 60 (2010) 38–54

51

Fig. 17. Log of relative error against M 1/3 .

In Fig. 17 a graph is plotted for log e  E R against M 1/3 . Note. In Example 7 the geometric mesh is used in the neighborhood of the crack tip, shown in Fig. 16(b) as circular rings. Modified polar coordinate system is used in this region for numerical computation and away from this region usual Cartesian coordinate system is used. For all other examples the numerical computations are done by only using modified polar coordinates as the geometry of the domains in these examples supports polar coordinates easily. 4. Conclusions The proposed method is exponentially accurate. The dimension of the Schur complement matrix is small, since the cardinality of the common boundary values is small. So it is easy to construct a nearly exact approximation to the Schur complement. The preconditioner is a block diagonal matrix where each diagonal block corresponds to an element and it’s inverse is trivial with almost optimal condition number. The algorithm for preconditioner is quite easy to implement with minimum extra effort. The residuals in the normal equation can be obtained efficiently without computing the mass and stiffness matrices. The method can also be implement on parallel computers more efficiently. Acknowledgements Authors gratefully acknowledge the help and guidance from Prof. P. Dutt and Prof. C.S. Upadhyay (I.I.T. Kanpur), without them this work would not have been possible. Appendix A A.1. Spectral element function Let uki,1 (τk , θk ) = hk , a constant on Ω˜ ik,1 (in the semi-infinite strip). Define the spectral element function

uki, j (τk , θk ) =

Wj Wj  

gr ,s τkr θks ,

(A.1)

r =0 s =0

˜ k for 1  i  I k , 2  j  N , 1  k  p. Here 1  W j  W . on Ω i, j p +1 p +1 Moreover there is an analytic mapping M l from the master square S = (−1, 1)2 to Ωl . Define p +1 

ul

p +1

Ml

W W    (ξ, η) = g r ,s ξ r η s . r =0 s =0

(A.2)

52

N.K. Kumar, G.N. Raju / Applied Numerical Mathematics 60 (2010) 38–54

A.2. Numerical scheme As described in Section 2, Ω˜ ik, j is the image of Ωik, j in (τk , θk ) coordinates. Let Lk be the operator defined by Lk u = rk2 Lu. Let y 1 = τk and y 2 = θk then

L˜k u = −

   2 2  ∂ ∂u ∂u a˜ ki, j b˜ ki + c˜ k u . + ∂ yi ∂yj ∂ yi

i , j =1

Let

 Ok =

i =1

cos θk

− sin θk

sin θk

cos θk

 and

˜k = A

 ˜k

a1,1

a˜ k1,2

a˜ k2,1

a˜ k2,2

 .

˜ k = ( O k )T A O k , where A is the matrix ( A )r ,s = ar ,s . Then A Next, let the vertex E k = (xk1 , xk2 ) and



F ik, j (τk , θk ) = e 2τk f xk1 + e τk cos θk , xk2 + e τk sin θk



in Ω˜ ik, j for 1  k  p, 2  j  N, 1  i  I k (since Lk u = rk2 f ). Let γs ⊆ Γ [1] ∩ ∂Ω k for 1  k  p, and γ˜s denote the image of point P˜ on γ˜s can be written as n = (n1 , n2 ). Then



∂ uk ∂n

 ˜k A

=

2 

ni a˜ ki, j

i , j =1

γs in (τk , θk ) coordinates. Now the unit normal n at a

∂ uk . ∂yj

Fig. 18. Edge Γk common to Ω k−1 and Ω k .

Let

Consider the boundary conditions u = gk on Γk ∩ ∂Ω k for k ∈ D , and ( ∂∂Nu ) A = gk on Γk ∩ ∂Ω k for k ∈ N . (See Fig. 18.)

 lk1 ( k )

τ =

u = gk (xk1 + e τk cos(ψlk ), xk2 + e τk sin(ψlk )),

( ∂∂ un ) A˜ k

for k ∈ D,

= e τk gk (xk1 + e τk cos(ψlk ), xk2 + e τk sin(ψlk )),

for k ∈ N .

Consider the boundary condition u = gk for k ∈ D , and ( ∂∂Nu ) A = gk for k ∈ N on Γk ∩ ∂Ω k−1 . Define



lk2 ( k−1 )

τ

=

u = gk (xk1−1 + e τk−1 cos(ψuk−1 ), xk2−1 + e τk−1 sin(ψuk−1 )),

( ∂∂ un ) A˜ k

p +1

Now we consider the elements in Ωl



 p +1  2 L u  dx1 dx2 =



l

p +1

Ωl

p +1

Here J l

for k ∈ D,

= e τk−1 gk (xk1−1 + e τk−1 cos(ψuk−1 ), xk2−1 + e τk−1 sin(ψuk−1 )) p +1

. In Ωl

for k ∈ N .

for 1  l  L,

 p +1  p +1 Lu J dξ dη . l l

S p +1

is the Jacobian of the mapping Ml

p +1

from S to Ωl

p +1

. Define Ll

=



p +1

Jl

L.

N.K. Kumar, G.N. Raju / Applied Numerical Mathematics 60 (2010) 38–54

53

 p +1 p +1 p +1 p +1 p +1 Let f l (ξ, η) = f ( Ml (ξ, η)) for 1  l  L and define F l (ξ, η) = f l (ξ, η) J l (ξ, η).

γs we shall denote a side common to the elements Ωmp+1 and Ωnp+1 . It may be assumed that γs is the image of η = −1 under the mapping Mmp+1 which maps S to Ωmp+1 and also the image of η = 1 under the mapping Mnp+1 which By

p +1

maps S to Ωn

 and



p +1 

um

x1

 p +1   p +1  = um ξ ξx1 + um η ηx1 ,

x2

 p +1   p +1  = um ξ ξx2 + um η ηx2 .

p +1 

um

. By the chain rule

Then

 p +1 2  p +1 2 p +1  u = um (ξ, −1) − un (ξ, 1)0, I , 0,γs   p +1 2   p +1  2  p +1   ux  =  um x (ξ, −1) − un (ξ, 1)1/2, I , i 1/2,γ x s

i

i

i = 1, 2.

Here I is the interval (−1, 1). Next, let γs ⊆ Γ [0] ∩ ∂Ω p +1 and let

γs be the image of η = −1 under the mapping Mmp+1 which maps S to Ωmp+1 . Then  p +1 2   p +1  2  ∂u     p +1 2  p +1     u  um (ξ, −1)2 +  ∂ um + = 1 ) . (ξ, −     , I 0,γs 0 ∂T ∂T 1/2,γs 1/2, I

In the same way if p +1

p +1

γs ⊆ Γ [1] ∩ ∂Ω p+1 , ( ∂ u∂ N ) A 21/2,γs can be defined. p +1

Let Γk ∩ ∂Ωm

p +1

k = Cm be the image of the mapping M m of S onto Ωm corresponding to the side η = −1, and p +1 k om (η) = gk ( M m (ξ, −1)), where −1  ξ  1. Let γs ⊆ Ω¯ k and d( E k , γs ) = infx∈γs {distance( E k , x)}. Choose λk < αk where αk is defined as in [7]. Let Fu = p +1 k {{u i , j (τk , θk )}i , j,k , {ul (ξ, η)}l } ∈ Π N ,W , the space of spectral element functions. Define ak = u ( E k ).

Define the functional N ,W rvertices (Fu ) =

Ik p N    

ρ qkN +1− j

k =1 j =2 i =1

+

p 



m 



i, j

2 d( E k , γs )−2λk  uk 0,γ˜ +  ukτk s





m∈D k=m−1 γs ⊆∂Ω k ∩Γm , μ(γ˜s )<∞ m  

+



2 

1/2,γ˜s

  2  +  ukθk 1/2,γ˜ s







 



m 2 2  k d( E k , γs )−2λk  uk − hk − lm m−k+1 − ak 0,γ˜s + u τk − lm−k+1 τk 1/2,γ˜s

(hk − ak )2 +

m∈D k=m−1

In the above Define

 

 

k=1 γs ⊆Ω k ∪ B kρ , μ(γ˜s )<∞

+

 −2λk  k  k  L˜ u (τk , θk ) − F k (τk , θk )2 ˜ k i, j i, j 0,Ω

m  



m∈N k=m−1 γs ⊆∂Ω k ∩Γm , μ(γ˜s )<∞

 k   ∂u

d( E k , γs )−2λk  

∂n

˜k A

2   − lm m−k+1 

.



(A.3)

1/2,γ˜s

μ(γ˜s ) denotes the measure of γ˜s .

N ,W rinterior (Fu ) =

L    p +1  p +1 2 p +1  L u (ξ, η) − F (ξ, η) l

l =1

+

l

    u p +1 2

0,γs

γs ⊆Ω p +1

+





l∈D γs ⊆∂Ω p +1 ∩Γl

+

l





l∈N

γs ⊆∂Ω p +1 ∩Γl

0, S

  p +1 2   p +1 2    +  u x1 +  u x2 1/2,γ 1/2,γ s

s

 p +1   l 2   ∂u ∂o    − + 0,γs ∂T ∂ T 1/2,γs

 2  p +1 u − ol 

 p +1  2  ∂u  l  −o  .  ∂N A 1/2,γs

(A.4)

54

N.K. Kumar, G.N. Raju / Applied Numerical Mathematics 60 (2010) 38–54

Let N ,W N ,W r N ,W (Fu ) = rvertices (Fu ) + rinterior (Fu ).

(A.5)

We choose as our approximate solution the unique Fz ∈ Π N , W , the space of spectral element functions, which minimizes the functional r N , W (Fu ) over all Fu . The numerical scheme presented is based on the stability estimate, Theorem 3.2 of [7]. The stability estimate in addition with the trace theorems for Sobolev spaces ensures the norm equivalence of residual norms and the solution norm. A.3. The coefficients in the solution of the Motz problem For the coefficients bl see Table 10. Table 10 Coefficients bl . l 0 1 2 3 4 5 6 7 8 9

bl 401.1624537452 87.6559201951 17.2379150794 −8.0712152597 1.4402727170 0.3310548859 0.2754373445 −0.0869329945 0.0336048784 0.0153843745

l

bl

10 11 12 13 14 15 16 17 18 19

0.0073023017 −0.0031841139 0.0012206461 0.0005309655 0.0002715122 −0.0001200463 0.0000505400 0.000023167 0.000011535 −0.000005295

References [1] I. Babuska, B.Q. Guo, Regularity of the solution of elliptic problems with piecewise analytic data, Part – I, SIAM J. Math. Anal. 19 (1988) 172–203. [2] I. Babuska, B.Q. Guo, The h − p version of the finite element method on domains with curved boundaries, SIAM J. Numer. Anal. 25 (1988) 837–861. [3] I. Babuska, H.S. Oh, The p version of the finite element method for domains with corners and infinite domains, Numer. Methods Partial Differential Equations 6 (4) (1990) 371–392. [4] P.B. Bochev, M.D. Gunzburger, Least-Squares Finite Element Methods, Springer, 2009. [5] Z. Cai, C. Westphal, A weighted H (div) least-squares method for second order elliptic problems, SIAM J. Numer. Anal. 46 (3) (2008) 1640–1651. [6] P. Dutt, P. Biswas, G. Naga Raju, Preconditioners for spectral element methods for elliptic and parabolic problems, J. Comput. Appl. Math. 215 (1) (2008) 152–166. [7] P. Dutt, N. Kishore Kumar, C.S. Upadhyay, Nonconforming h − p spectral element methods for elliptic problems, Proc. Indian Acad. Sci. Math. Sci. 117 (2007) 109–145. [8] G. Karniadakis, S.J. Spencer, Spectral/hp Element Methods for CFD, Oxford University Press, 1999. [9] V.A. Kondratiev, The smoothness of a solution of Dirichlet’s problem for second order elliptic equations in a region with a piecewise smooth boundary, Differ. Uravn. 6 (10) (1970) 1831–1843, also: Differ. Equ. 6 (1970) 1392–1401. [10] N. Kishore Kumar, P. Dutt, C.S. Upadhyay, Nonconforming spectral/hp element methods for elliptic systems, J. Numer. Math. 17 (2) (2009) 119–142. [11] E. Lee, T.A. Manteuffel, C.R. Westphal, Weighted-norm first order system least-squares (FOSLS) for problems with corner singularities, SIAM J. Numer. Anal. 44 (5) (2006) 1974–1996. [12] Z.C. Li, R. Mathon, P. Sermer, Boundary methods for solving elliptic problems with singularities and interfaces, SIAM J. Numer. Anal. 24 (3) (1987) 487–498. [13] T.R. Lucas, H.S. Oh, The method of auxiliary mapping for the finite element solutions of elliptic problems containing singularities, J. Comput. Phys. 108 (1993) 327–342. [14] H.S. Oh, I. Babuska, The method of auxiliary mapping for the finite element solutions of elasticity problems containing singularities, J. Comput. Phys. 121 (1995) 193–212. [15] H.S. Oh, H. Kim, S.J. Lee, The numerical methods for oscillating singularities in elliptic boundary value problems, J. Comput. Phys. 170 (2001) 742–763. [16] D. Pathria, G.E. Karniadakis, Spectral element methods for elliptic problems in nonsmooth domains, J. Comput. Phys. 122 (1995) 83–95. [17] J.B. Rosser, N. Papamichael, MRC Technical Summary Report, No. 1405. 1975, University of Wisconsin. [18] Ch. Schwab, p and h − p Finite Element Methods, Clarendon Press, Oxford, 1998. [19] P. Seshaiyer, Non-conforming hp finite element methods, PhD dissertation, University of Maryland Baltimore Country, 1998. [20] P. Seshaiyer, M. Suri, Convergence results for non-conforming hp methods: The mortar finite element method, Contemp. Math. 218 (1998). [21] B. Szabo, I. Babuska, Finite Element Analysis, John Wiley & Sons, 1991. [22] S.K. Tomar, h − p spectral element methods for elliptic problems on non-smooth domains using parallel computers, PhD thesis, IIT Kanpur, India, 2001. Reprint available as Tec. Rep. No. 1631, Department of Applied Mathematics, University of Twente, The Netherlands, http://www.math.utwente.nl/ publications. [23] S.K. Tomar, h − p spectral element method for elliptic problems on non-smooth domains using parallel computers, Computing 78 (2006) 117–143.