Mesh grading in isogeometric analysis

Computers and Mathematics with Applications (

)

–

Contents lists available at ScienceDirect

Computers and Mathematics with Applications journal homepage: www.elsevier.com/locate/camwa

Mesh grading in isogeometric analysis Ulrich Langer, Angelos Mantzaflaris, Stephen E. Moore, Ioannis Toulopoulos ∗ Johann Radon Institute for Computational and Applied Mathematics, Austrian Academy of Sciences, Altenbergerstr. 69, A-4040 Linz, Austria

article

info

Article history: Available online xxxx Keywords: Elliptic boundary value problems Domains with geometric singular points or edges Discontinuous coefficients Isogeometric analysis Mesh grading Recovering optimal convergence rates

abstract This paper is concerned with the construction of graded meshes for approximating socalled singular solutions of elliptic boundary value problems by means of multipatch discontinuous Galerkin Isogeometric Analysis schemes. Such solutions appear, for instance, in domains with re-entrant corners on the boundary of the computational domain, in problems with changing boundary conditions, in interface problems, or in problems with singular source terms. Making use of the analytic behavior of the solution, we construct the graded meshes in the neighborhoods of such singular points following a multipatch approach. We prove that appropriately graded meshes lead to the same convergence rates as in the case of smooth solutions with approximately the same number of degrees of freedom. Representative numerical examples are studied in order to confirm the theoretical convergence rates and to demonstrate the efficiency of the mesh grading technology in Isogeometric Analysis. © 2015 Elsevier Ltd. All rights reserved.

1. Introduction The gradient of the solution of elliptic boundary value problems can exhibit singularities in the vicinities of re-entrant corners or edges. The same is true in case of changing boundary conditions or interface problems. This singular behavior of the gradients was discovered and analyzed in the famous work by Kondrat’ev [1]. We refer the reader to the monographs [2–4] for a more recent and comprehensive presentation of related results. It is well known that these singularities may cause loss in the approximation order of the standard discretization methods like the finite element method, see the classical monograph [5] or the more recent paper [6]. In the case of two dimensional problems with singular boundary points, grading mesh techniques have been developed for finite element methods in order to recover the full approximation order, see the classical textbook [7] and the more recent publications [6,8–10] for three-dimensional problems. Here, we devise graded meshes for solving elliptic problems with singular solutions by means of discontinuous Galerkin Isogeometric Analysis method (dG IgA). In the frame of IgA, the use of B-splines or NURBS basis functions allows complicated CAD (Compute-Aided Design) geometries to be exactly represented. The key point of Hughes et al. [11] was to make use of the same basis to approximate the solution of the problem under consideration. Since this pioneer paper, applications of the IgA approach have been considered in many fields, see [12]. Here, we apply a multipatch symmetric dG IgA method which has been extensively studied for diffusion problems in volumetric computational domains and on surfaces in [13,14], respectively, see also [15] for a comprehensive presentation. The solution of the problem is independently approximated in every subdomain by IgA without imposing any matching grid conditions and without any continuity requirements for the discrete solution across

∗

Corresponding author. E-mail addresses: [email protected] (U. Langer), [email protected] (A. Mantzaflaris), [email protected] (S.E. Moore), [email protected] (I. Toulopoulos). http://dx.doi.org/10.1016/j.camwa.2015.03.011 0898-1221/© 2015 Elsevier Ltd. All rights reserved.

2

U. Langer et al. / Computers and Mathematics with Applications (

)

–

the subdomain interfaces. Symmetrized numerical fluxes with interior penalty jump terms, see, e.g., [16–18], are introduced on the interfaces in order to treat the discontinuities of the discrete solution and to interchange information between the non-matching grids. As we will see later, considering a numerical scheme in this general context is a flexible approach, notably when applied on zone-type subdivisions of Ω , which are convenient for analysing elliptic boundary value problems in domains with singular boundary points. This paper aims at the construction of graded dG IgA meshes in the zones located near the singular points in order to recover full convergence rates like in the case of smooth solutions on uniform meshes. The grading of the mesh is mainly determined by the analytic behavior of the solution u around the singular points and follows the spirit of grading mesh techniques using layers, which have been proposed for finite element methods in [7,6,8]. According to this, having an a priori knowledge about the location of the singular point, e.g. the re-entrant corner, the domain Ω is subdivided into zones, called layers in [6,8], and then a further subdivision of Ω into subdomains (also called patches in IgA), say TH (Ω ) := {Ωi }Ni=1 , is performed in such a way that TH (Ω ) is in correspondence with the initial zone partition. On the other hand, the solution can be split into a sum of a regular part ur ∈ W l≥2,2 (Ω ) and a singular part us ∈ W 1+ε,2 (Ω ), with known ε ∈ (0, 1), i.e., u = ur + us , see, e.g., [2]. The analytical form of us contains terms with singular exponents in the radial direction. We use this information and construct appropriately graded meshes in the zones around the singular points. The resulting graded meshes have a ‘‘zone-wise character’’, this means that the grid size of the graded mesh in every zone determines the mesh of every subdomain which belongs to this zone, where we assume that every subdomain belongs to only one zone (the ideal situation is every zone to be a subdomain). We mention that the mesh grading methodology is developed and is analyzed for the classical two dimensional problem with a re-entrant corner. The proposed methodology can be generalized and applied to other situations. This is shown by the numerical examples presented in Section 4. The particular properties of the produced graded meshes help us show optimal error estimates for the dG IgA method, which exhibit optimal convergence rates. The error estimates for the proposed method are proved by using a variation of Céa’s lemma and using B-spline quasi-interpolation estimates for u ∈ W 1,2 (Ω ) ∩ W l≥2,p∈(1,2] (TH (Ω )), which have been proved in [13]. More precisely, these interpolation estimates are expressed with respect to the mesh size hi of the corresponding subdomain Ωi . For the domains away from the singular point, the solution is smooth (see ur part in previous splitting), and we can derive the usual interpolation estimates. Conversely, for the subdomains Ωi , for which the boundary ∂ Ωi touches the singular point, the singular part us of the solution u can be considered as a function from the Sobolev space W 2,p∈(1,2) (Ωi ). Now the estimates given in [13] enable us to derive error estimates for the singular part us . This makes the whole error analysis easier in comparison with the techniques earlier developed for the finite element method, e.g., in [6,19,9]. We mention that, in the literature, other IgA techniques have been proposed for solving problems with singularities. In [20,21], the mapping technique has been developed, where the original B-spline finite dimensional space has been enriched by generating singular functions which resemble the types of the singularities of the problem. The mappings constructed on this enriched space describe the geometry singularities explicitly. Also in [22], by studying the anisotropic character of the singularities of the problem, the one-dimensional approximation properties of the B-splines are generalized for two-dimensional problems, in order to produce anisotropic refined meshes in the regions of the singular points. Furthermore, we note the extended isogeometric element formulation, which has been developed for solving fracture mechanics problems in [23–26]. Utilizing knot insertion, the original NURBS basis is enriched by crack-discontinuous functions, which are able to reproduce the singular field near the crack. In [27–29], the flexibility of the PHT-spline and T-spline basis on local adaptive refinement has been used for solving elasticity problems on non-smooth domains. Recently, other types of adaptivity refinement procedures have been presented by using hierarchical B-splines, see [30]. The authors introduce different levels of resolution in an adaptive framework, keeping simultaneously the good properties of standard B-splines. The same methodology has been generalized for truncated hierarchical B-splines and has been applied in CAD modeling as well in [31]. Lastly, new adaptive refinement techniques have been developed in the IgA framework on the basis of functional a posteriori error estimates in [32]. The rest of the paper is organized as follows. The problem description, the weak formulation and the dG IgA discrete analogue are presented in Section 2. Section 3 discusses the construction of the appropriately graded IgA meshes, and provides the proof for obtaining the full approximation order of the dG IgA method on the graded meshes. Several two and three dimensional examples are presented in Section 4. Finally, we draw our conclusions.

2. Problem description and dG IgA discretization First, let us introduce some notation. We define the differential operator α

α

Da = D1 1 · · · Dd d ,

with Dj =

∂ , D(0,...,0) u = u, ∂ xj

where α = (α1 , . . . , αd ), with αj ≥ 0, j = 1, . . . , d, denotes a multi-index of the degree |α| =

(2.1)

d

j =1

αj . For a bounded

Lipschitz domain Ω ⊂ Rd , d = 2, 3 we denote by W l,p (Ω ), with l ≥ 1 and 1 ≤ p ≤ ∞, the usual Sobolev function spaces

U. Langer et al. / Computers and Mathematics with Applications (

)

–

3

endowed with the norms

∥u∥W l,p (Ω ) =

 

∥Dα u∥pLp (Ω )

 1p

,

(2.2a)

0≤|α|≤m

∥u∥W l,∞ (Ω ) = max ∥Dα u∥∞ .

(2.2b)

0≤|α|≤m

More details about Sobolev’s function spaces can be found in [33]. We often write a ∼ b, meaning that Cm a ≤ b ≤ CM a, with Cm and CM are positive constants independent of the discretization parameters. 2.1. The model problem Let us assume that the boundary of ΓD = ∂ Ω of Ω contains geometric singular parts. In particular, for d = 2, we consider domains which have corner boundary points with internal angles greater than π . For d = 3, we consider that case where the domain Ω can be described in the form Ω = Ω2 × Z , where Ω2 ⊂ R2 and Z = [0, zM ] is an interval. The cross section of Ω has only one corner with an interior angle ω ∈ (π , 2π ). This means that the ∂ Ω has only one singular edge which is Γs := {(0, 0, z ), 0 ≤ z ≤ zM }. The remaining parts of ΓD are assumed to be smooth, see Fig. 2(a) and (b) for an illustration of the domains. For simplicity, we restrict our study to the following model problem:

− div(α∇ u) = f

in Ω , u = uD on ∂ Ω ,

(2.3)

where the coefficient α(x) ∈ L∞ (Ω ) is a piecewise constant function, bounded from above and below by positive constants, 1

and f ∈ L2 (Ω ), uD ∈ H 2 (∂ Ω ) are given data. The variational formulation of (2.3) reads as follows: find u ∈ W 1,2 (Ω ) such that u = uD on ΓD = ∂ Ω and a(u, v) = l(v),

∀v ∈ W01,2 (Ω ),

(2.4a)

where a(u, v) =

 Ω

α∇ u · ∇v dx and l(v) =

 Ω

f v dx.

(2.4b)

It is clear that, under the assumptions made above, there exists a unique solution of the variational problem (2.4) due to Lax–Milgram’s lemma, [34]. We follow the theoretical analysis of the regularity of solution presented in [3] and we introduce few theoretical tools which are necessary. We consider the two-dimensional case. Suppose that the ΓD has only one singular corner, say Ps , with internal angle ω ∈ (π , 2π ), and that the boundary parts from the one and the other side of Ps are smooth curves, see Fig. 2(a). We consider the local cylindrical coordinates (r , θ ) with origin Ps , and define the cone (a circular sector with angular point Ps )

C = {(x, y) ∈ Ω : x = r cos(θ ), y = r sin(θ ), 0 < r < R, 0 < θ < ω}.

(2.5)

We construct a highly smooth cut-off function ξ in C , such that ξ ∈ C , and it is supported inside the cone C . It has been shown in [3], that the solution u of the problem (2.4) can be written as a sum of a regular function ur ∈ W l≥2,2 (Ω ) and a singular function us , ∞

u = ur + us ,

(2.6)

us = ξ (r )γ r λ sin(λθ ),

(2.7)

with

where γ is the stress intensity factor (for the two-dimensional problems is a real number depending only on f ) and λ = πω ∈ (0, 1) is an exponent which determines the strength of the singularity. Since λ < 1, based on the theory and the examples of Ch. 5, pp. 244–245 of [34], we can show that the singular function us does not belong to W 2,2 (Ω ) but to W l=2,p (Ω ) with p = 2/(2 − λ). Consequently, the regularity properties of u in C are mainly determined by the regularity properties of us , and we can assume that u ∈ W 1,2 (Ω ) ∩ W l,p (TH (Ω )) (see below details for the TH (Ω )). Remark 2.1. For the expression (2.7), we admit that the computational domain has only one non-convex corner and that only Dirichlet boundary conditions are prescribed on ∂ Ω . A similar expression can be derived if there are more non-convex corners or other types of boundary conditions, see details in [3].

4

U. Langer et al. / Computers and Mathematics with Applications (

)

–

2.2. The dG IgA discrete scheme 2.2.1. Isogeometric analysis spaces

¯ = i=1 Ω ¯i We assume a non-overlapping subdivision TH (Ω ) := {Ωi }Ni=1 of the computational domain Ω such that Ω with Ωi ∩ Ωj = ∅ for i ̸= j. The subdivision TH (Ω ) is considered to be compatible with the discontinuities of the coefficient α , i.e., the jumps can only appear on the interfaces Fij = ∂ Ωi ∩ ∂ Ωj between the subdomains. For the sake of brevity in our notations, the set of common interior faces are denoted by FI . The collection of the faces that belong to ∂ Ω are denoted by FB , i.e., F ∈ FB , if there is a Ωi ∈ TH (Ω ) such that F = ∂ Ωi ∩ ∂ Ω . We denote the set of all subdomain faces by F = FI ∪ FB . In the multi-patch IgA context, each subdomain is represented by a B-spline (or NURBS) mapping. To accomplish this, we associate each Ωi with a vector of knots 4di = (Ξi1 , . . . , Ξiι , . . . , Ξid ), with Ξiι = {ξ1ι , ξ2ι , . . . , ξnι }, ι = 1, . . . , d, which are N

i ˆ   = (0, 1)d . The interior knots of 4di form a mesh T (i) = {Eˆ m }M set on the parametric domain Ω m=1 in Ω , where Em are the h i ,Ω

(i)

ˆ micro elements, see Fig. 1. Given a micro element Eˆ m ∈ Th ,Ω  , we denote by hEˆ m the diameter of Em , and the local grid size hi i is defined to be hi = max{hEˆm }. Finally, we define the maximum grid size to be h = max{hi }. We refer the reader to [12] for more information about the meaning of the knot vectors in CAD and IgA. (i)

d Assumption 1. The meshes Th ,Ω  defined by the knots 4i are quasi-uniform, i.e., there exists a constant σ ≥ 1 such that

σ −1 ≤

hˆ

Em

hˆ

Em+1

i

≤ σ. (i)

(i)

ˆ On each Th ,Ω  , we derive the finite dimensional space Bhi spanned by B-spline (or NURBS) basis functions of degree k, see i

more details in [12,35,36], (i)

(i)

ˆ (i) ) dim(B h

ˆ h = span{Bˆ j (ˆx)}j=1 B i

i

.

(2.8a)

(i)

ˆ in (2.8a) is the product of one-dimensional B-spline basis functions, i.e. Every function, Bˆ j (ˆx), xˆ ∈ Ω (i) (i) (i) Bˆ j (ˆx) = Bˆ j1 (ˆx1 ) · · · Bˆ jd (ˆxd ),

(2.8b)

where the index jι corresponds to a univariate B-spline basis function in direction xι . In the following, we suppose that the one-dimensional B-splines in (2.8b) have the same degree k. Finally, having the B-spline spaces, we can represent each subdomain Ωi by the parametric mapping

 → Ωi , 8i : Ω

8i (ˆx) =



(i) (i)

Cj Bˆ j (ˆx) := x ∈ Ωi ,

(2.9a)

j 1 with xˆ = 8− i (x),

(2.9b)

(i)

where Cj are the B-spline control points, i = 1, . . . , N, cf. [12]. (i)

M

i We construct a mesh Thi ,Ωi = {Em }m= 1 for every Ωi , whose vertices are the images of the vertices of the corresponding

(i)

parametric mesh Th ,Ω  through 8i , see Fig. 1. Notice that, the above subdomain mesh construction can result in non-matching i

meshes along the patch interfaces. Further, by taking advantage of the properties of 8i , we define the global finite dimensional B-spline (dG) space (i)

(N )

Bh (TH ) := Bhi (Ωi ) × · · · × BhN (ΩN ), (i)

(2.10a)

(i)

where every Bhi (Ωi ) is defined on Thi ,Ωi as follows: (i) (i) (i) (i) 1 Bhi (Ωi ) := {Bj |Ωi : Bj (x) = Bˆ j ◦ 8− i (x),

(i)

(i)

ˆ h }. for Bˆ j ∈ B i

(2.10b)

Later, the solution u of the problem (2.4) will be approximated by the discrete solution uh ∈ Bh (TH ). 2.2.2. Discrete problem The problem (2.4) is independently discretized in every Ωi using the spaces (2.10b), without imposing continuity requirements for the B-spline basis functions on the interfaces Fij = ∂ Ωi ∩ ∂ Ωj and also non-matching grids may exist. (i)

Using the notation φh := φh |Ωi , we define the average and the jump of φh ∈ Bh (TH ) on Fij ∈ FI by

{φh } :=

1 2

(φh(i) + φh(j) ),

and

[[φh ]] := φh(i) − φh(j) ,

(2.11a)

U. Langer et al. / Computers and Mathematics with Applications (

)

–

5

Fig. 1. The parametric domain and two adjacent subdomains with different underlying meshes red and blue. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 2. The domains, (a) two-dimensional with corner singularity, (b) three-dimensional with re-entrant edge, and (c) subdivision of Ω into zones and subdomains.

and, for Fi ∈ FB ,

{φh } := φh(i) ,

(i)

and [[φh ]] := φh .

(2.11b)

Since the functions in Bh (TH ) are not continuous across the subdomain interfaces, we apply the discontinuous Galerkin technology. We multiply (2.3) by φh ∈ Bh (TH ), integrate over the computational domain Ω , rewrite the integral over Ω as a sum of integrals over the subdomains Ωi , and then integrate by parts over every Ωi . After few computations, we arrive at N   i =1

Ωi

α∇ u∇φh dx −

 

  {α∇ u} · nFij [[φh ]] ds = f φh dx,

(2.12)

Ω

Fij

Fij ⊂∂ Ωi

where the unit normal vector nFij is oriented from Ωi towards the interior of Ωj . Exploiting the fact that [[u]] = 0 on the solution of our variational problem (2.4), we add two terms to (2.12) which are zero on the solution and which make the dG bilinear form symmetric and stable, see [13,18]. This standard dG procedure leads to the following symmetric dG IgA variational problem: find uh ∈ Bh (TH ) such that ah (uh , φh ) = l(φh ) + pD (uD , φh ),

∀φh ∈ Bh (TH ),

(2.13a)

where the dG bilinear form is given by ah (uh , φh ) =

N  

 1

ai (uh , φh ) −

Fij ⊂∂ Ωi

i =1

2

si (uh , φh ) + pi (uh , φh )



(2.13b)

with the bilinear forms (cf. also [16]): ai (uh , φh ) = si (uh , φh ) =

 Ωi



α∇ uh ∇φh dx,

{α∇ uh } · nFij [[φh ]] + {α∇φh } · nFij [[uh ]] ds, Fij

(2.13c)

(2.13d)

6

U. Langer et al. / Computers and Mathematics with Applications (

pi (uh , φh ) =

  (j) σ α (i)  σα   + [[uh ]][[φh ]] ds,   hj

Fij

    pD (uD , φh ) =

σ α (i)



hi

Fi

σ α (i)

 Fi

hi

hi

[[uh ]][[φh ]] ds,

uD φh ds,

)

–

if Fij ∈ FI , (2.13e) if Fij ∈ FB ,

F i ∈ FB .

(2.13f)

We again mention that the term {α∇φh } · nFij [[uh ]] that appears in (2.13d) ensures the symmetry of the dG bilinear form, and the jump terms in (2.13e) help to achieve the discrete coercivity of the dG bilinear form ah (., .). The penalty parameter

σ > 0 must be chosen large enough, e.g., σ >

 (k+1)(k+2)  21 2

, in order to ensure the stability of the dG IgA method [13].

3. IgA on graded meshes In many realistic applications, we very often have to solve problems similar to (2.4) in domains with non-smooth boundary parts, that possess geometric singularities, for instance, non-convex corners, see Fig. 2. It is well-known that the numerical methods loose accuracy when they are applied to this type of problems. This occurs as a result of the reduced regularity of the solutions in the vicinity of the non-smooth parts [2]. When finite element methods are used, graded meshes have been utilized around the singular boundary parts in order to obtain optimal convergence rates, see e.g. [6,10], see also [9] for dG methods. The basic idea of this grading mesh technique is to use the a priori knowledge of the singular behavior of the solution around the singular boundary points, cf. (2.7) and (2.6), and consequently adjust accordingly the size of the elements. In this section, we extend the grading mesh techniques of the finite element method to dG IgA framework for solving boundary value problems like (2.4) in the presence of singular points. Our approach is inspired by the grading mesh methodology using layers, cf. [6,8]. In the next section, we construct the graded mesh and show that the proposed dG IgA method exhibits optimal convergence rates as for problems with high regularity solutions. We present the technique and the corresponding analysis for two-dimensional problems. However, in Section 4, we also apply our methodology to some three-dimensional examples and discuss the numerical results. 3.1. A priori mesh grading The grading of the meshes around the singular points is guided by the exponent λ, which specifies the regularity of the function us , see (2.7), as well as by the location of the singular boundary point. Next, we discuss the construction of the mesh for the case of one singular geometric point on ∂ Ω . Let Ps be the singular point and let Us := {x ∈ Ω : |Ps − x| ≤ R = LU h, with LU ≥ 2} be an area around Ps in Ω , which is 1

further subdivided into ζM ring-type zones Zζ , ζ = 0, . . . , ζM , such that the distance from Ps is D(Zζ ,Ps ) := C (nζ h) µ , where C =R

1 1− µ

and 0 ≤ nζ < LU . By µ ∈ (0, 1], we denote the grading control parameter. The radius of every zone is defined to 1

1

be RZζ := D(Zζ +1 ,Ps ) − D(Zζ ,Ps ) = C (nζ +1 h) µ − C (nζ h) µ , where we suppose that there is a ν > 0 such that nζ +1 = nζ + ν with 1 ≤ ν < LU − 1. In particular, we set RZM = R − D(ZM −1 ,Ps ) . For convenience, we assume that the initial subdivision TH (Ω ) fits to the Zζ ring zone partition in order to fulfill the following conditions, for an illustration, see Fig. 2(c) with ζM = 3:

• The subdomains can be grouped into those which belong (entirely) to the area Us and those which belong (entirely) to Ω \ Us . This means that there is no Ωi , i = 1, . . . , N such that Us ∩ Ωi ̸= ∅ and (Ω \ Us ) ∩ Ωi ̸= ∅. • Every ring zone Zζ is partitioned into ‘‘circular’’ subdomains Ωiζ , which have radius RΩiζ equal to the radius of the zone, that is RΩi = RZζ . For computational efficiency reasons, we prefer, if it is possible, every zone to be only represented by ζ one subdomain. This essentially depends on the characteristics of the problem, i.e., the shape of Ω and the coefficient α . (i ) • The zone Z0 is represented by one subdomain, say Ωi0 , and the mesh Thi 0 (Ωi0 ) includes all the micro-elements E such 0 that ∂ E ∩ Ps ̸= ∅. (iζ )

In the sequel, we will construct the meshes Thi (Ωiζ ) such that they satisfy the following property: for ζ ̸= 0, the mesh size ζ

1

1−µ

hiζ is of order hiζ = O (hRΩi ), where D(Zζ ,Ps ) is the distance between Ωiζ and Ps , and hi0 is of order O (h µ ). Thus, we have ζ

the following relations: 1

1

Cm h µ ≤ hi0 ≤ CM h µ , 1−µ

Cm hRΩi

ζ

if Ω iζ ∩ Ps ̸= ∅, 1−µ

≤ hiζ ≤ CM hD(Zζ ,Ps ) ,

if Ω iζ ∩ Ps = ∅.

(3.1a) (3.1b)

U. Langer et al. / Computers and Mathematics with Applications (

)

–

7

(iζ )

We need to specify the mesh size for every Thi (Ωiζ ) in order to satisfy inequalities (3.1). Having predefined the maximum ζ

grid size h, see Section 2.2.1, we set 1

1

((nζ + ν)h) µ − (nζ h) µ

hiζ = C

,

1

[ν µ ]

where C = R

(3.2) (iζ )

1 1− µ

and [ν −µ ] denotes the nearest integer to ν −µ . Consequently, the mesh size of Thi (Ωiζ ) is of order ζ

hiζ = O (RZζ ν −(1/µ) ). Notice that the grading has ‘‘a subdomain character’’ and is mainly determined by the parameter µ ∈ (0, 1]. For µ = 1, we have hiζ = h, that is, we obtain quasi-uniform meshes. Using the fact that µ ≤ 1 and the inequality (a + b)γ ≤ 2γ −1 (aγ + bγ ), for γ ≥ 1 and a, b ≥ 0 (which can be shown easily since the function t γ is convex in (0, ∞)) we arrive at the estimates 1

hi ζ = C

1

((nζ + ν)h) µ − (nζ h) µ 1

1

≤C

≤ (C1,R,µ,ν nζ )

1

1

≤

1 (C1,R,µ,ν nζ ) µ  1−µ

nζ 1

1

1

[ν µ ] µ hh µ −1

1

Cµ (nζ h) µ + Cµ (ν h) µ − (nζ h) µ

h

µ

nζ h

[ν µ ]  1 1−µ µ

1−µ

≤ (C2,µ,ν nζ ) µ hD(Zζ ,Ps ) ,

(3.3)

which gives the right inequality in (3.1b). Since 1 > 1 − µ ≥ 0, by the definition of the subdomain radius RΩi , we can easily ζ

show that 1

1−1+µ  1  RΩ ≥ C m h, iζ µ ν

(3.4a)

and immediately we get 1

 1  RΩi ≥ Cm R1Ω−µ h, iζ ζ νµ with Cm =

1 2

1

(3.4b)

1

((nζ + ν) µ − nζµ ). Inserting the definition of the grid size hiζ = RΩiζ /[ν −µ ] in (3.4b), we derive the left

inequality in (3.1b). Remark 3.1. There exist several possible strategies for building graded meshes in finite elements. The interested reader is referred to [8,6] for an overview of different possibilities. We have chosen the construction presented in this section due to both its simplicity and suitability for IgA. 3.2. Quasi-interpolant, error estimates Next, we study the error estimates of the method (2.13). For the purposes of our analysis, we consider the enlarged space l ,p

Wh := W 1,2 (Ω ) ∩ W l≥2,p (TH (Ω )) + Bh (TH (Ω )),

(3.5)

where p ∈ (max{1, d+22d }, 2]. Let us mention that we allow different l and p in different subdomains Ωi . In particular, for (l−1) 2 subdomains with Ωi ∩ Us = ∅ we set in (3.5) p = 2, while for subdomains with Ωi ∩ Us ̸= ∅ we set 1 < p = 2−λ < 2. The 1,2

space Wh

is equipped with the broken dG-norm

∥u∥2dG(Ω ) =

N  

 α (i) ∥∇ u(i) ∥2L2 (Ω ) + pi (u(i) , u(i) ) , i

1,2

u ∈ Wh .

(3.6)

i =1

Let f ∈ W 1,2 (Ω ) ∩ W l≥2,p (TH (Ω )) with p ∈ (max{1, d+22d }, 2], then we can construct a quasi-interpolant Πh f ∈ Bh (TH ) (l−1) such that Πh f = f for all f ∈ Bh (TH ). We refer the reader to [36], see also [35], for more details about the construction of Πh f . We have the following approximation estimate.

8

U. Langer et al. / Computers and Mathematics with Applications (

)

–

(i)

Lemma 3.2. Let u ∈ W 1,2 (Ω ) ∩ W l,p (TH (Ω )) with p ∈ (max{1, d+22d }, 2] and l ≥ 2, and let E = 8i (Eˆ ), Eˆ ∈ Th ,Ωˆ . Then, (l−1) i   for 0 ≤ m ≤ l ≤ k + 1, there exist an quasi-interpolant Πh u ∈ Bh (TH ) and constants Ci := Ci maxl0 ≤l (∥Dl0 8i ∥L∞ (Ωi ) ) such that

|u − Πh u|pW m,p (E ) ≤ Ci ∥u∥W l,p (Ωi ) hpi (l−m) .

 (i) E ∈Th ,Ω i

(3.7)

i

Furthermore, we have the following estimates in the ∥.∥dG(Ω ) norm

∥u − Πh u∥dG(Ω ) ≤

N  

δ(l,p,d)

C i hi

N     hi  δ(l,p,d) ∥u∥W l,p (Ωi ) + Ci α (j) hi ∥u∥W l,p (Ωi ) , i=1 Fij ⊂∂ Ωi

i=1

∥uh − Πh u∥dG(Ω ) ≤ ∥u − Πh u∥dG(Ω ) +

N 

δ(l,p,d)

C i hi

hj

∥u∥W l,p (Ωi ) ,

(3.8a)

(3.8b)

i =1

where δ(l, p, d) = l + ( 2d −

d p

− 1).

Proof. The proof is given in [13].



Remark 3.3. If hi and hj are the grid sizes of two adjacent subdomains Ωi and Ωj , then relations (3.1) immediately yield the two-side estimate σm,ζ ≤ hi /hj ≤ σM ,ζ , where the positive constants σm,ζ and σM ,ζ only depend on the quantities which specify the initial zone partition Zζ . Hence, in what follows, the estimate (3.8a) will be used in the form ∥u − Πh u∥dG(Ω ) ≤

N

i=1

δ(l,p,d)

C i hi

. We mention that the analysis presented here can easily be extended to non-matching grids, see [13]. We (iζ )

also note that the meshes Thi (Ωiζ ) satisfy Assumption 1. ζ

We emphasize that Lemma 3.2 provides local estimates which hold in every subdomain Ωi . These help us investigate the accuracy of the method in every zone Zζ of Us separately. Since us ∈ W

2 l=2,p= 2−λ

estimate (3.8a) for l = 2, d = 2 and p = 2−λ , we obtain

(Ωiζ ) for Ωiζ ∈ Zζ , by the local interpolation

2

∥us − Πh us ∥dG(Us ) ≤



hλiζ Ciζ .

(3.9)

iζ

Next, we prove the discretization error estimate under the assumption that ur ∈ W l,2 (Ω ),

with l ≥ k + 1.

(3.10)

Theorem 3.4. Suppose that (3.10) holds and let Zζ be a zone partition of Ω with the properties listed in the previous section. The (i)

meshes of the subdomains Thi (Ωi ) are as described in Section 3.1. Then, for the solution u of (2.4a), we have the error estimate

∥u − uh ∥dG(Ω ) ≤ Chr ,

with r = min{k, λ/µ},

(3.11)

where the constant C > 0 depends only on the quasi-uniform mesh properties, see (3.1), and the constants Ci of Lemma 3.2. Proof. Let Πh u ∈ Bh (TH ) be the quasi-interpolant of Lemma 3.2. Using the triangle inequality, we obtain

∥u − uh ∥dG(Ω ) ≤ ∥uh − Πh u∥dG(Ω ) + ∥u − Πh u∥dG(Ω ) .

(3.12)

Moreover, representation (2.6) yields

∥u − Πh u∥dG(Ω ) ≤ ∥us − Πh us ∥dG(Ω ) + ∥ur − Πh ur ∥dG(Ω ) .

(3.13)

Using (3.10), Lemma 3.2, the mesh properties (3.1) and 0 < µ ≤ 1, we have k

∥ur − Πh ur ∥dG(Ω ) ≤ C1 h µ + C2 hk ≤ Chk ,

(3.14)

where the constants C1 and C2 are determined by the constants that appear in (3.1) and (3.8a). Therefore, it remains to estimate the first term in (3.13). By (3.9) and (3.1a), we obtain the estimate λ

∥us − Πh us ∥dG(Z0 ) ≤ Chλi0 ≤ Ch µ

(3.15)

U. Langer et al. / Computers and Mathematics with Applications (

)

–

9

(iζ )

Fig. 3. Basis functions on the graded mesh Th ,Ω  using µ = 0.6: (a) 1D case and (b) 2D case. iζ

in Z0 , where the constant C is determined by the constants in (3.1) and (3.8a). For the subdomains belonging to the remaining zones of Zζ , ζ ̸= 0, (3.9) and (3.1b) yield the estimates

 1−µ λ ∥us − Πh us ∥dG(Zζ ) ≤ Cζ ̸=0 hλiζ ≤ Cζ ̸=0 hD(Zζ ,Ps )  1−µ λ λ ≤ Cζ ̸=0 hh µ ≤ Cζ ̸=0 h µ .

(3.16)

Collecting (3.14)–(3.16), we arrive at the interpolation error estimate

∥u − Πh u∥dG(Ω ) ≤ Chr ,

(3.17)

with r = min{k, λ/µ}. Now, inserting estimate (3.17) into (3.12) and recalling estimate (3.8b), we derive the error estimate (3.11).  4. Numerical examples In this section, we present a series of numerical examples in order to confirm the theoretical results and to assess the effectiveness of the proposed grading mesh technique. First we look at two-dimensional problems with boundary point singularities and with highly discontinuous coefficients. In the last examples, we consider applications of the method to three-dimensional problems with an interior singularity and in domains with singular edges having ω = 3π /2 interior angle. The numerical examples have been performed in G + Smo.1 4.1. Implementation details All the numerical tests presented below, are performed using the fewest possible subdomains, and, according to the spirit of IgA, the grading of the mesh is done in the parameter domain. The underlying assumption is that the given parameterization of the domain has uniform speed along the patch. An ideal situation is to have an arc-length parameterization. Nevertheless, a well-behaving parameterization is one whose speed is within a constant factor of the arc-length parameterization. This is a reasonable assumption, also because CAD software typically try to adhere to such a requirement, since it is desirable for CAD operations as well. For constructing the graded parameter mesh, we choose a number of interior knots in each parametric direction and we place the knots according to the grading parameter and the location of the singular point in parameter space. This means that, in the parametric mesh the knots are located in such a way that the size of the edges of the micro elements which are close to singular point satisfies (3.1a), and similarly the edges of the micro elements which are away from the singular point satisfy (3.1b). In Fig. 3(a), we show a one-dimensional B-spline basis on a graded mesh, and similarly in Fig. 3(b), we present a two-dimensional one on the corresponding graded mesh. If the location of the singular point is given in physical coordinates, we invert the point to parameter space with a Newton iteration, to map it back to parameter space. Under the (iζ )

(i)

assumption of a well-behaved B-spline geometry map, the parameter mesh is transformed from Th ,Ω  to Thi ,Ωi in such a iζ

1 Geometry + Simulation Modules, http://www.gs.jku.at/gismo, see also [37].

10

U. Langer et al. / Computers and Mathematics with Applications (

)

–

Fig. 4. Heart shaped problem: (a) the computational domain and the subdomains, (b) the graded meshes of the two subdomains, and (c) the contours of uh solution. Table 1 The control points for the two B-spline surfaces each with degree k = 2 depicted in Fig. 4. These points where constructed by interactive design, using the G + Smo library and the user interface of Axel geometric modeler (http://axel.inria.fr). j1

j2

B1j (j1 , j2 )

B2j (j1 , j2 )

1 2 3 1 2 3 1 2 3

1 1 1 2 2 2 3 3 3

(0.00, 0.00) (0.49, 0.49) (0.97, 0.97) (0.00, −0.81) (0.46, −0.16) (1.00, 0.94) (0.35, −0.84) (0.71, −0.84) (0.85, 0.042)

(0.00, 0.00) (0.49, 0.49) (0.97, 0.97) (−0.81, 0.00) (−0.16, 0.46) (0.94, 1.00) (−0.84, 0.35) (−0.84, 0.71) (0.04, 0.85)

way that the size of the physical elements are proportional to their (pre-image) parametric elements. This property ensures that our theoretical analysis applies in the experiments that we conducted. For efficiency reasons, the mesh TH (Ω ) is created using the same number of knots as an equivalent uniform mesh with grid size hi for each subdomain Ωi , where each knot is relocated towards the singularity Ps according to the grading parameter µ. This strategy satisfies Assumption 1 without increasing the total number of knots. This approach also mimics the zones construction introduced in Section 3.1, since it corresponds to a zone partition that shrinks towards Ps at every refinement step. In our experiments, we consider a mapping 8i produced on an initial knot vector 4di , which exactly represents the subdomain Ωi , as the isogeometric paradigm suggests. During the relocation of the knots the shape and the parameterization of the subdomains Ωi remain unaltered. Nevertheless, in our implementation, we have the freedom to add additional knots in the graded mesh in order to preserve the initial shape of the subdomain. In the previous analysis, we assumed that each mapping from the parametric domain to the physical subdomains is regular, see (2.9). However, as we will see in the heart example, singularities in the mapping do not affect the procedure, provided that they are not near to the singularities of the solution. 4.2. Numerical examples in two dimensions 4.2.1. Heart shaped domain To illustrate the efficiency of the proposed mesh grading methodology and to validate the estimates of Section 2, we consider the problem (2.4) in a curved domain (heart shaped) having a singular point Ps (re-entrant corner) with internal angle ω = 3π /2, see Fig. 4(a). The computational domain Ω consists of two subdomains shown in Fig. 4(a), where the corresponding knot vectors are 42i = (Ξi1 , Ξi2 ), i = 1, 2 with Ξi1 = Ξi2 = {0, 0, 0, 1, 1, 1} and are parametrized by k = 2 B-spline basis with the control π

points given in Table 1. The exact solution is given by u = r ω sin(θ π /ω) with f and uD in (2.4) are specified by the exact solution. We set α = 1 in the entire Ω . Note that u ∈ W 1.5,2 (Ω ). The problem is solved using first (k = 1) and second (k = 2) order B-spline spaces with grading parameters µ = 0.6 and µ = 0.3, respectively, see Fig. 4(b). We plot the contours of the solution uh computed using quadratic B-splines in Fig. 4(c). In Table 2, the relative error ∥u − uh ∥dG /∥u∥dG computed on every mesh is presented. We observe that the relative error associated with the graded meshes is much more reduced (specifically for the k = 2 case) in comparison with the relative error of the non-graded meshes. In particular, when quadratics are used, two levels of graded mesh refinement resulted in the same magnitude of the relative error as five levels

U. Langer et al. / Computers and Mathematics with Applications (

)

–

11

Table 2 Heart shaped problem: the relative error ∥uh − u∥dG(Ω ) /∥u∥dG(Ω ) with and without grading. h/2s

Without grading k=1

With grading k=2

k = 1, µ = 0.6

k = 2, µ = 0.3

0.0795555 0.049578 0.0310918 0.0195497 0.0123059 0.0077496

0.16685 0.0932192 0.0501646 0.0264986 0.0138272 0.0071509

0.0654141 0.0215324 0.0061074 0.0016288 0.0004300 0.0001128

Relative error s s s s s s

=0 =1 =2 =3 =4 =5

0.190142 0.119593 0.074713 0.046712 0.029257 0.018355

of uniform refinement. This is due to the fact that the graded mesh assists in increasing the resolution of the approximate solution around the singular corners. In Table 3, we display the convergence rate of the error. In the left column (without grading), we present the rates using quasi-uniform meshes. In the right column of the table, we show the rates in the case of using graded meshes. As the theory predicts, the convergence rates in left columns of both cases k = 1 and k = 2 are mainly determined by the regularity of the solution. On the other hand, the rates which correspond to the graded meshes tend to be optimal with respect to the B-spline degree k. This shows the adequacy of the proposed graded mesh for solving this problem. The present example has a few boundary points with vanishing Jacobian determinant of the parameterization. For instance, two such points are close to (0, 0.5) and (0.5, 0), see Fig. 4(b), where (curved) triangular elements can be observed. These singularities in the mapping did not affect the convergence rates. However, if the singularity of the mapping was located near to the re-entrant corner, we would expect numerical problems. In particular, the quadrature points would be located too close to the singularity of the mapping, and thus we might be lead to division by (almost) zero while performing the numerical integration. In any case, if the required grading is too strong, then after several steps of refinement we may arrive at having degenerate elements close to singular boundary points, which would have to be excluded or treated differently. For example, in [24,26], mappings with singularities are used to solve elliptic problems with singular solutions. However, all these cases are beyond the scope of the present paper. 4.2.2. Kellogg’s problem It is known that the solutions of problem (2.4) with rough diffusion coefficients may not be very smooth. Thus, standard numerical methods cannot provide an (optimal) accurate approximation [38]. We examine such a case by solving the socalled Kellogg test problem [39]. We consider the computational domain Ω = (−1, 1)2 . Also, we choose a piecewise constant diffusion coefficient α in (2.3), having the same value, say α := α13 , in the first and third quadrants and, similarly, α := α24 in the second and fourth quadrants. This choice of the diffusion coefficient leads to the subdivision of Ω in the four subdomains as it is shown in Fig. 5(a). The exact solution of the problem for f = 0 is given in polar coordinates by u(r , θ ) = r λ ϕ(θ ), where

 cos((π /2 − σ )λ) cos((θ − π /2 + ρ)λ),   cos(ρλ) cos((θ − π + σ )λ), ϕ(θ ) = cos(σ λ) cos((θ − π − ρ)λ),  cos((π /2 − ρ)λ) cos((θ − 3π /2 − σ )λ),

if 0 ≤ θ < π /2, if π /2 ≤ θ < π , if π ≤ θ < 3π /2, if 3π /2 ≤ θ ≤ 2π ,

and the numbers λ, ρ, σ satisfy the nonlinear relations

− tan((π /2 − σ )λ) cot(ρλ) = R,    − tan(ρλ) cot(σ λ) = 1/R,   − tan(σ λ) cot((π /2 − ρ)λ) = R, 0 < λ < 2 ,    max{0, π λ − π } < 2λρ < min{π λ, π}, max{0, π − π λ} < −2λσ < min{π , 2π − λπ},

(4.1)

where R := α13 /α24 is the ratio of the diffusion coefficient values. The above system of equations admits several solutions. For this example, we set λ = 0.4 and ρ = π /4. A Newton iteration recovers one solution for the rest of the parameters, namely σ = 4.7123889 and R = 9.4721359. Therefore, we set the diffusion coefficients to α13 = R and α24 = 1. On the one hand, we obtain a problem with an exact solution u ∈ W 1.4,2 (Ω ) and discontinuous derivatives across the interfaces. On the other hand, u ∈ W 2,1.25 (Ω ), and the estimates presented in Section 3 can be applied. We solved the problem using B-spline spaces with degrees k = 1 and k = 2 on uniform meshes. We performed again the test using graded meshes with grading parameter chosen such that λ/µ = k, see (3.11). In Fig. 5(b), we can see the graded meshes of the subdomains. Fig. 5(c) shows the contours of the dG solution uh computed for linear B-splines. In Table 4, we display the convergence rates of the error. We observe that, in the case of uniform meshes, the experimental order of convergence of the method approaches 0.4, which agrees with the regularity

12

U. Langer et al. / Computers and Mathematics with Applications (

)

–

Table 3 Heart shaped problem: the convergence rates in ∥.∥dG(Ω ) with and without grading. h/2s

Without grading k=1

With grading k=2

k = 1, µ = 0.6

k = 2, µ = 0.3

– 0.68221 0.67322 0.669385 0.667797 0.667156

– 0.843026 0.894636 0.921219 0.938709 0.951475

– 1.49519 1.85785 2.02913 2.02562 2.00987

Convergence rates s s s s s s

=0 =1 =2 =3 =4 =5

– 0.671469 0.678694 0.677558 0.675018 0.672622

Table 4 Kellogg’s test: the convergence rates in the dG-norm ∥.∥dg (Ω ) . h/2s

Without grading k=1

With grading k=2

k = 1, µ = 0.40

k = 2, µ = 0.20

– 0.591165 0.355586 0.378796 0.385672 0.390348

– 0.830217 0.858329 0.879976 0.895984 0.906179

– 0.477657 1.15442 1.78696 1.84425 1.95223

Convergence rates s s s s s s

=0 =1 =2 =3 =4 =5

– 0.655814 0.354865 0.368103 0.378375 0.385464

Fig. 5. Kellogg’s test problem: (a) the computational domain and the four subdomains, (b) the graded meshes of the four subdomains, and (c) the contours of uh solution.

of the solution. Conversely, the rates in the right columns which correspond to the results using graded meshes tend to be optimal with respect the order of the B-spline space. A glimpse of the discrete solution uh is given in Fig. 6. The grading required for recovering a quadratic rate of convergence in this example is quite strong, and reaches the limits of our implementation. Smaller grading parameters lead to undesirable effects. One effect of this ‘‘over-grading’’, i.e. µ < 0.2 in our case, is that elements close to the singular points tend to become degenerate (i.e. have area close to zero). Numerical quadrature rules can be problematic in this case, due to division by a very small number in the denominator of the integrand. In addition, the linear solver needs a careful preconditioning. Nevertheless, our simple ILU preconditioned Conjugate Gradient solver was sufficient for tackling this example even for µ = 0.2. 4.3. Three dimensional examples 4.3.1. Cube with interior point singularity This test case is inspired by [40]. The computational domain is Ω = (−1, 1)3 and is decomposed into eight subdomains, see Fig. 7(a). We choose the diffusion coefficient α = 1 in the whole computational domain Ω . The solution of the problem has a singular point at the origin of the axis and is given by u(x) = |x|λ with λ = 0.85. It is easy to show that u ∈ W l=2,2 (Ω ). We solved the problem using k = 1, k = 2 and k = 3 B-spline spaces on quasi-uniform meshes. In Fig. 7(b), we plot the contours of solution uh computed using linear B-splines. The convergence rates of the error corresponding to non-graded meshes are shown in left columns of Table 6.

U. Langer et al. / Computers and Mathematics with Applications (

)

–

13

Fig. 6. Discrete solution of Kellogg’s problem plotted over the graded mesh.

Fig. 7. Cube with interior singularity: (a) the decomposition of Ω into 8 subdomains with the graded meshes of the subdomain, (b) the contours of the solution uh , and (c) the variance of the uh contours around the singular point. Table 5 Cube with interior singularity: the relative error ∥uh − u∥dG(Ω ) /∥u∥dG(Ω ) . h/2s

Without grading k=1

With grading

k=2

k=3

k = 2, µ = 0.6

k = 3, µ = 0.4

0.11820 0.05402 0.02175 0.00860 0.00338 0.00132

0.04164 0.02289 0.00958 0.00377 0.00147 0.00058

0.11820 0.04346 0.01293 0.00341 0.00087 0.00022

0.0416433 0.0208749 0.00565056 0.00077189 9.24677e−05 1.11311e−05

Relative error s s s s s s

=0 =1 =2 =3 =4 =5

0.33124 0.21410 0.11989 0.06361 0.03289 0.01677

The rates are optimal for the linear B-spline space, but not for the two other B-spline spaces. This was expected, according to the regularity of the solution u. Note that the rates presented on the left columns in Table 6 are in agreement with the estimate given in (3.8a). We have performed again the test using graded meshes for the last two B-spline spaces. The grading parameter µ has been chosen to be δ(l = 2, p = 2, d = 3)/µ = k, see Lemma 3.2 and (3.11). In Fig. 7(c), the variate of the uh contours around the singular point is shown. ∥u−u ∥ The variations of the relative error ∥u∥h dG computed on the successively refined meshes are given in Table 5. From dG the comparison of the results in Table 5, we can see that the reduction of the error computed on the graded meshes is higher. This indicates that higher accuracy solutions can be achieved using the proposed graded meshes. The corresponding convergence rates obtained on the graded meshes are displayed in the right columns of Table 6. We can observe that these rates approach the optimal rate for both higher-order B-spline spaces. This numerical example demonstrates that the dG IgA method applied on the proposed graded meshes can exhibit optimal convergence rates for interior singularity type problems as well. 4.3.2. Three-dimensional L-shape domain   Now the computational domain Ω has the form of a 3D L-shape and is given by (−1, 1)2 \ (−1, 0)2 ×[0, 1]. Even though the ‘‘L-shape’’ example has been mostly studied in the literature in its two-dimensional set up (see for example anisotropic 2D meshes for IgA discretizations in [22]) we believe that it is an interesting test case, because we will see that the graded

14

U. Langer et al. / Computers and Mathematics with Applications (

)

–

Table 6 Cube with interior singularity: the convergence rate of the error on uniform and graded meshes. h/2s

Without grading k=1

With grading

k=2

k=3

k = 2, µ = 0.6

k = 3, µ = 0.4

– 1.066 1.306 1.340 1.346 1.348

– 0.687 1.234 1.343 1.350 1.350

– 1.393 1.766 1.928 1.959 1.974

– 0.791 1.870 2.942 3.080 3.066

Convergence rates s s s s s s

=0 =1 =2 =3 =4 =5

– 0.593 0.839 0.917 0.953 0.972

Table 7 3D L-shape: the convergence rates of the error with respect to the dG norm on uniform and graded meshes. h/2s

Without grading k=1

With grading k=2

k = 1, µ = 0.6

k = 2, µ = 0.3

– 0.477178 0.639951 0.670841 0.669949 0.668371

– 0.629909 0.869128 0.883655 0.902467 0.920065

– 0.387338 1.11198 1.80531 1.96533 2.00296

Convergence rates s s s s s s

=0 =1 =2 =3 =4 =5

– 0.645078 0.650805 0.642971 0.644107 0.648100

Fig. 8. 3D L-shape test: (a) The domain Ω with the corners and the edge boundary singularities, (b) The graded meshes of the two subdomain, and (c) The contours of uh .

mesh of the plane can be prolonged in a direction perpendicular to the singular edge for treating the boundary singularities. Note that in this three-dimensional setting, the domain includes both corner and edge singularities, see Fig. 8(a). We consider an exact solution given by u = r λ sin( θωπ ), where λ = π /ω and ω = 3π /2. We set ΓD = ∂ Ω . The data f and uD of (2.3) are given by the exact solution. The computational domain Ω consists of two subdomains. We have solved the problem using B-spline spaces of order k = 1 and k = 2 using quasi-uniform and graded meshes in both subdomains. The grading parameter is defined by the relation δ(l, p, d)/µ = k. In Fig. 8(b), we can see the graded meshes for µ = 0.6. The contours of the corresponding approximate solution uh computed for k = 1 are presented in Fig. 8(c). Table 7 displays the convergence rates of the error. We observe the same behavior of the rates as in the previous examples. The rates of the uniform meshes are determined by the regularity of the solution (u ∈ W 1+λ,p=2 (Ω )) for both B-spline spaces. The convergence rates corresponding to graded meshes approach the optimal value. We remark here that similar grading approaches have been used in finite element methods for approximating solutions of elliptic problems in three-dimensional domains with edges, see [6,8,41]. 4.3.3. Three-dimensional heart shaped domain In this example, we consider an exact solution given by u = r λ sin(θ π /ω), where λ = π /ω and ω = 3π /2. We again set ΓD = ∂ Ω , and the data f and uD of (2.3) are specified by the given exact solution. The computational domain Ω consists of two subdomains. The problem is solved with B-spline spaces of order k = 1 and k = 2 using quasi-uniform and graded meshes in both subdomains. The grading parameter is defined by the relation λ/µ = k. In Fig. 9(b), we can see the graded meshes for µ = 0.6. The contours of the corresponding approximate solution uh computed with degree k = 1 is presented in Fig. 9(c).

U. Langer et al. / Computers and Mathematics with Applications (

)

–

15

Fig. 9. 3D heart test: (a) the domain Ω with the corners and edge boundary singularities, (b) the graded meshes of the two subdomain, and (c) the contours of uh . Table 8 3D heart: the convergence rates of the error with respect to the dG norm on uniform and graded meshes. h/2s

Without grading k=1

With grading k=2

k = 1, µ = 0.6

k = 2, µ = 0.3

– 0.675611 0.685756 0.674337 0.669817 0.667968

– 0.686633 0.846524 0.902119 0.925762 0.94134

– 0.964287 1.55143 1.91781 2.10561 2.09457

Convergence rates s s s s s s

=0 =1 =2 =3 =4 =5

– 0.650805 0.642971 0.644107 0.6481 0.65251

The convergence rates of the error corresponding to the quasi-uniform meshes are shown in the left columns (without grading) of Table 8, and the rates corresponding to the graded meshes are shown in the right columns of Table 8. 5. Conclusion We have presented mesh grading techniques for dG IgA discretizations of elliptic boundary value problems in the presence of so-called singular points. Based on the a priori or a posteriori knowledge of the behavior of the exact solution around the singular points, we pre-defined the grading of the mesh by performing a relocation of the knots, instead of increasing their number. The grading has been devised in a patch-wise fashion in order to fit well in the IgA framework. Optimal error estimates of the multipatch dG IgA method have been shown, when used in conjunction with the proposed graded meshes. The theoretical results have been confirmed by a number of two- and three-dimensional test problems with known exact solutions. At the first glance, mesh grading is an a priori adaptive mesh generation technique. It requires a certain knowledge about the singularities of the solution. However, this knowledge can be obtained after a few uniform refinement steps with linear B-splines, and via computing the local convergence rates in the neighborhood of critical geometrical points, cf., also the left part ‘‘without grading’’ of the tables presented in Section 4. This adaptive mesh grading strategy was used in the finite element context in [42] in connection with multigrid solvers, see also [43,44]. On the basis of these ideas and our results, it is certainly possible to construct a fully adaptive mesh grading IgA multilevel method. Acknowledgment This research was supported by the National Research Network NFN S117-03 ‘‘Geometry + Simulation’’ of the Austrian Science Fund (FWF). References [1] V.A. Kondrat’ev, Boundary value problems for elliptic equations in domains with conical or angular points, Transl. Moscow Math. Soc. 16 (1967) 227–313. [2] P. Grisvard, Elliptic Problems in Nonsmooth Domains, in: Monographs and Studies in Mathematics, Pitman Advanced Pub. Program, 1985. [3] P. Grisvard, Singularities in Boundary Value Problems, in: Recherches en Mathématiques Appliquées, Masson, 1992. [4] V.A. Kozlov, V.G. Maz’ya, J. Rossmann, Spectral Problems Associated with Corner Singularities of Solutions to Elliptic Equations, in: Mathematical Surveys and Monographs, vol. 85, American Mathematical Society, Rhode Issland, USA, 2001. [5] G. Strang, G. Fix, An Analysis of the Finite Element Method, Prentice-Hall, Englewood Cliffs, NJ, 1973. [6] T. Apel, A.-M. Sändig, 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 (30) (1996) 63–85.

16

U. Langer et al. / Computers and Mathematics with Applications (

)

–

[7] L.A. Oganesjan, L.A. Ruchovetz, Variational Difference Methods for the Solution of Elliptic Equations, Isdatelstvo Akademi Nank Armjanskoj SSR, Erevan, 1979, (in Russian). [8] T. Apel, F. Milde, Comparison of several mesh refinement strategies near edges, Comm. Numer. Methods Engrg. 12 (1996) 373–381. [9] M. Feistauer, A.-M. Sändig, Graded mesh refinement and error estimates of higher order for DGFE solutions of elliptic boundary value problems in polygons, Numer. Methods Partial Differential Equations 28 (4) (2012) 1124–1151. [10] T. Apel, B. Heinrich, Mesh refinement and windowing near edges for some elliptic problem, SIAM J. Numer. Anal. 31 (3) (1994) 695–708. [11] T.J.R. Hughes, J.A. Cottrell, Y. Bazilevs, Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement, Comput. Methods Appl. Mech. Engrg. 194 (2005) 4135–4195. [12] J.A. Cotrell, T.J.R. Hughes, Y. Bazilevs, Isogeometric Analysis, Toward Integration of CAD and FEA, John Wiley and Sons, 2009. [13] U. Langer, I. Toulopoulos, Analysis of multipatch discontinuous Galerkin IgA approximations to elliptic boundary value problems. RICAM Reports 201408, Johann Radon Institute for Computational and Applied Mathematics, Austrian Academy of Sciences, Linz, 2014, http://arxiv.org/abs/1408.0182. [14] U. Langer, S.E. Moore, Discontinuous Galerkin isogeometric analysis of elliptic PDEs on surfaces. NFN Technical Report 12, Johannes Kepler University Linz, NFN Geometry and Simulation, Linz, 2014, http://arxiv.org/abs/1402.1185 and (in press) in the DD22 proceedings. [15] U. Langer, A. Mantzaflaris, S.E. Moore, I. Toulopoulos, Multipatch discontinuous Galerkin isogeometric analysis. RICAM Reports 2014-18, Johann Radon Institute for Computational and Applied Mathematics, Austrian Academy of Sciences, Linz, 2014, http://arxiv.org/abs/1411.2478v1. [16] M. Dryja, On discontinuous Galerkin methods for elliptic problems with discontinuous coeffcients, Comput. Methods Appl. Math. 3 (2003) 76–85. [17] B. Rivière, Discontinuous Galerkin Methods for Solving Elliptic and Parabolic Equations: Theory and Implementation, SIAM, Philadelphia, 2008. [18] Daniele A. Di Pietro, Alexandre Ern, Mathematical Aspects of Discontinuous Galerkin Methods, in: Mathématiques et Applications, vol. 69, SpringerVerlag, Heidelberg, Dordrecht, London, New York, 2012. [19] T. Apel, Interpolation of non-smooth functions on anisotropic finite element meshes, M2AN 33 (6) (1999) 1149–1185. [20] H.S. Oh, H. Kim, J.W. Jeong, Enriched isogeometric analysis of elliptic boundary value problems in domains with cracks and/or corners, Internat. J. Numer. Methods Engrg. 97 (3) (2014). [21] J.W. Jeong, H.S. Oh, S. Kang, H. Kim, Mapping techniques for isogeometric analysis of elliptic boundary value problems containing singularities, Comput. Methods Appl. Mech. Engrg. 254 (0) (2013) 334–352. [22] L. Beirão da Veiga, D. Cho, G. Sangalli, Anisotropic NURBS approximation in isogeometric analysis, Comput. Methods Appl. Mech. Engrg. 209–212 (2012) 1–11. [23] E. De Luycker, D.J. Benson, T. Belytschko, Y. Bazilevs, M.C. Hsu, X-FEM in isogeometric analysis for linear fracture mechanics, Internat. J. Numer. Methods Engrg. 87 (6) (2011) 541–565. [24] S.Sh. Ghorashi, N. Valizadeh, S. Mohammadi, Extended isogeometric analysis for simulation of stationary and propagating cracks, Internat. J. Numer. Methods Engrg. 89 (9) (2012) 1069–1101. [25] S.Sh. Ghorashi, N. Valizadeh, S. Mohammadi, T. Rabczuk, T-spline based XIGA for fracture analysis of orthotropic media, Comput. Struct. 147 (0) (2015) 138–146. [26] N. Nguyen-Thanh, N. Valizadeh, M.N. Nguyen, H. Nguyen-Xuan, X. Zhuang, P. Areias, G. Zi, Y. Bazilevs, L. De Lorenzis, T. Rabczuk, An extended isogeometric thin shell analysis based on Kirchhoff-love theory, Comput. Methods Appl. Mech. Engrg. 284 (2015) 265–291. Isogeometric Analysis Special Issue. [27] N. Nguyen-Thanh, H. Nguyen-Xuan, S.P.A. Bordas, T. Rabczuk, Isogeometric analysis using polynomial splines over hierarchical t-meshes for twodimensional elastic solids, Comput. Methods Appl. Mech. Engrg. 200 (21–22) (2011) 1892–1908. [28] Y. Bazilevs, V.M. Calo, J.A. Cottrell, J.A. Evans, T.J.R. Hughes, S. Lipton, M.A. Scott, T.W. Sederberg, Isogeometric analysis using T-splines, Comput. Methods Appl. Mech. Engrg. 199 (2010) 229–263. [29] N. Nguyen-Thanh, J. Kiendl, H. Nguyen-Xuan, R. Wüchner, K.U. Bletzinger, Y. Bazilevs, T. Rabczuk, Rotation free isogeometric thin shell analysis using PHT-splines, Comput. Methods Appl. Mech. Engrg. 200 (47–48) (2011) 3410–3424. [30] A.-V. Vuong, C. Giannelli, B. Jüttler, B. Simeon, A hierarchical approach to adaptive local refinement in isogeometric analysis, Comput. Methods Appl. Mech. Engrg. 200 (49–52) (2011) 3554–3567. [31] G. Kiss, C. Giannelli, U. Zore, B. Jüttler, D. Grossmann, J. Barner, Adaptive {CAD} model (re-)construction with THB-splines, Graph. Models 76 (5) (2014) 273–288. [32] S.K. Kleiss, S.K. Tomar, Guaranteed and sharp a posteriori error estimates in isogeometric analysis, 2014. arXiv:1304.7712. [33] R.A. Adams, J.J.F. Fournier, Sobolev Spaces, second ed., in: Pure and Applied Mathematics, vol. 140, ACADEMIC PRESS-imprint Elsevier Science, 2003. [34] L.C. Evans, Partial Differential Equestions, first ed., in: Graduate Studies in Mathematics, vol. 19, American Mathematical Society, 1998. [35] Y. Bazilevs, L. Beirão da Veiga, J.A. Cottrell, T.J.R. Hughes, G. Sangalli, Isogeometric analysis: approximation, stability and error estimates for h-refined meshes, M3AS 16 (07) (2006) 1031–1090. [36] L.L. Schumaker, Spline Functions: Basic Theory, third ed., University Press, Cambridge, 2007. [37] B. Jüttler, U. Langer, A. Mantzaflaris, S.E. Moore, W. Zulehner, Geometry + simulation modules: implementing isogeometric analysis, PAMM 14 (1) (2014) 961–962. [38] S.R. Falkand, J.E. Osborn, Remarks on mixed finite element methods for problems with rough coefficients, Math. Comp. 62 (205) (1994) 1–19. [39] R.B. Kellogg, On the Poisson equation with intersecting interfaces, Appl. Anal. 4 (1975) 101–129. [40] D. Kröner, M. Růžička, I. Toulopoulos, Numerical solutions of systems with (p, δ)-structure using local discontinuous Galerkin finite element methods, Internat. J. Numer. Methods Fluids 76 (11) (2014) 855–874. [41] T. Apel, B. Heinrich, The finite element method with anisotropic mesh grading for elliptic problems in domains with corner and edges, SIAM J. Numer. Anal. 31 (3) (1998) 695–708. [42] M. Jung, On adaptive grids in multilevel methods, in: S. Hengst (Ed.), GAMM–Seminar on Multigrid-Methods, Gosen, Germany, 21–25 September 1992, in: WIAS Report, vol. 5, WIAS, 1993, pp. 67–80. [43] V.V. Scheidurov, Multigrid Methods for Finite Elements, Nauka, Moscow, 1989, (in Russian). [44] M. Jung, S. Nicaise, J. Tabka, Some multilevel methods on graded meshes, J. Comput. Appl. Math. 138 (1) (2002) 151–171.