{'fraction_non_alphanumeric': 0.0469390368852459, 'fraction_numerical': 0.012743340163934426, 'mean_word_length': 4.069469242006168, 'pattern_counts': {'":': 0, '<': 7, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Mesh Grading in Isogeometric Analysis

This paper is concerned with the construction of graded meshes for approximating so-called 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.


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 [19]. We refer the reader to the monographs [14,15,20] 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 [29] 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 [25] and the more recent publications [6,5,13], and [3] 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 IgA frame, the use of B-splines or NURBS basis functions allow complicated CAD geometries to be exactly represented, and the key point of Hughes et al. [16] was to make use of the same basis to approximate the solution of the problem under consideration. Since this pioneer paper, applications of IgA method have been considered in many fields, see [9]. 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 [24] and [23], respectively, see also [22] for 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 the subdomain interfaces. Symmetrized numerical fluxes with interior penalty jump terms, see, e.g., [11,27,10], 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, the consideration of the numerical scheme in this general context makes it more flexible to be applied on zone-type subdivisions of Ω, which have been found to be quite convenient for treating 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 [25,6,5]. 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,5], and then a further subdivision of Ω into subdomains (also called patches in IgA), say T H (Ω) := {Ω i } N i=1 , is performed in such way that T H (Ω) is in correspondence with the initial zone partition. On the other hand, the solution can be split into a sum of a regular part u r ∈ W l≥2,2 (Ω) and a singular part u s ∈ W 1+ε,2 (Ω), with known ε ∈ (0, 1), i.e., u = u r + u s , see, e.g., [14]. The analytical form of u s 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 into 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 to 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 quasiinterpolation estimates for u ∈ W 1,2 (Ω) ∩ W l≥2,p∈ (1,2] (T H (Ω)), which have been proved in [24]. More precisely, these interpolation estimates have subdomain character and are expressed with respect to the mesh size h i of the corresponding subdomain Ω i . For the domains away from the singular point, the solution is smooth (see u r 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 u s 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 [24] enable us to derive error estimates for the singular part u s . This makes the whole error analysis easier in comparison with the techniques earlier developed for the finite element method, e.g., in [6,2,13].
We mention that, in the literature, other IgA techniques have been proposed for solving two-dimensional problems with singularities very recently. In [26] and [17], 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 [8], by studying the anisotropic character of the singularities of the problem, the onedimensional 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.
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 some conclusion.

Problem description and dG IgA discretization
First, let us introduce some notation. We define the differential operator where α = (α 1 , ..., α d ), with α j ≥ 0, j = 1, ..., d, denotes a multi-index of the degree |α| = d j=1 α j . For a bounded Lipschitz domain Ω ⊂ R d , d = 2, 3 we denote by W l,p (Ω), with l ≥ 1 and 1 ≤ p ≤ ∞, the usual Sobolev function spaces endowed with the norms More details about Sobolev's function spaces can be found in [1]. We often write a ∼ b, meaning that C m a ≤ b ≤ C M a, with C m and C M are positive constants independent of the discretization parameters.

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 ⊂ R 2 and Z = [0, z M ] 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 The remaining parts of Γ D are considered as smooth, see Fig. 2(a) and Fig. 2(b) for an illustration of the domains. For simplicity, we restrict our study to the following model problem where the coefficient α(x) ∈ L ∞ (Ω) is a piecewise constant function, bounded from above and below by positive constants, f ∈ L 2 (Ω) and u D ∈ H 1 2 (∂Ω) are given data. The variational formulation of (2.3) reads as follows: find u ∈ W 1,2 (Ω) such that u = u D on Γ D = ∂Ω and 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. We follow the theoretical analysis of the regularity of solution presented in [15]. We consider the two-dimensional case. Suppose that the Γ D has only one singular corner, say P s , with internal angle ω ∈ (π, 2π), and that the boundary parts from the one and the other side of P s are straight lines, see Fig. 2(a). We consider the local cylindrical coordinates (r, θ) with origin P s , and define the cone (a circular sector with angular point P s ). 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 [15], that the solution u of the problem (2.4) can be written as a sum of a regular function u r ∈ W l≥2,2 (Ω) and a singular function u s , with u s = ξ(r)γr λ sin(λθ), 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, by an easy computation, we can show that the singular function u s 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 u s , and we can assume that u ∈ W 1,2 (Ω) ∩ W l,p (T H (Ω)), (see below details for the T H (Ω)). Remark 1. For the expression (2.7), we admit that the computational domain has only one non-convex corner and only Dirichlet boundary conditions are prescribed on ∂Ω. Similar expression can be derived if there are more non-convex corners and if there are other type of boundary conditions, see details in [15].
The subdivision T H (Ω) is considered to be compatible with the discontinuities of the coefficient α, i.e., the jumps can only appear on the interfaces F ij = ∂Ω i ∩ ∂Ω j between the subdomains. For the sake of brevity in our notations, the set of common interior faces are denoted by F I . The collection of the faces that belong to ∂Ω are denoted by In the multi-patch (multi-subdomain) 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 .., ξ ι n }, ι = 1, . . . , d, which are set on the parametric domain Ω = (0, 1) d . The interior knots of Ξ d i are considered without repetitions and form a mesh T in Ω, whereÊ m are the micro elements. Given a micro elementÊ m ∈ T (i) h i , Ω , we denote by hÊ m = diameter(Ê m ), and the local grid size h i is defined to be the maximum diameter of allÊ m ∈ T We refer the reader to [9] for more information about the meaning of the knot vectors in CAD and IgA.

Assumption 1 The meshes T
h i spanned by B-spline (or NURBS) basis functions of degree k, see more details in [9,7,28],B (i) (2.8a) j (x) function in (2.8a) is derived by means of tensor products of one-dimensional B-spline basis functions, i.e. (2.8b) 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 j are the B-spline control points, i = 1, ..., N , cf. [9]. We construct a mesh T for every Ω i , whose vertices are the images of the vertices of the corresponding parametric mesh T Notice that, the above subdomain mesh construction can result in non-matching meshes along the patch interfaces.
Further, by taking advantage of the properties of Φ i , we define the global finite dimensional B-spline (dG) space Later, the solution u of the problem (2.4) will be approximated by the discrete (dG) solution u h ∈ B h (T H ).

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 F ij = ∂Ω i ∩ ∂Ω j and also non-matching grids may exist. Using the notation φ The parametric domain and two adjacent subdomains with different underlying meshes red and blue.

and, for
The discrete problem is specified by the symmetric dG IgA method, see [24], and reads as follows: where the dG bilinear form is given by with the bilinear forms (cf. also [11]): Here the unit normal vector n F ij is oriented from Ω i towards the interior of Ω j . The penalty parameter µ > 0 must be chosen large enough in order to ensure the stability of the dG IgA method [24].

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 [14]. 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,3], see also [13] 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. The purpose of this paper is to extend the grading mesh techniques from the finite element method to dG IgA framework for solving boundary value problems like (2.4) in the presence of singular points. We develop a mesh grading algorithm around the singular boundary parts inspired by the grading mesh methodology using layers, therefore, extending the approach used in finite element methods, cf. [6,5], to isogeometric analysis.
Next, we construct the graded mesh and show that the proposed dG IgA method exhibits optimal convergence rates as for the problems with high regularity solutions. We present our mesh grading technique and the corresponding analysis for two-dimensional problems. In Section 4, we also apply our methodology to some three-dimensional examples and discuss the numerical results.

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 u s , see (2.7), and by the location of the singular boundary point too. Next, we discuss the construction of the mesh for the case of one singular geometric point on ∂Ω.
Let P s be the singular point and let U s := {x ∈ Ω : |P s − x| ≤ R = L U h, with L U ≥ 2} be an area around P s in Ω, which is further subdivided into ζ M ring-type zones Z ζ , ζ = 0, .., ζ M , such that the distance from µ and 0 ≤ n ζ < L U . By µ ∈ (0, 1], we denote the grading control parameter. The radius of every zone is defined to be where we suppose that there is a ν > 0 such that n ζ+1 = n ζ + ν with 1 ≤ ν < L U − 1. In particular, we set For convenience, we assume that the initial subdivision T H (Ω) 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) into the area U s and those that belong (entirely) into Ω \ U s . This means that there is no Ω i , i = 1, ..., N such that U s ∩ Ω i = ∅ and (Ω \ U s ) ∩ Ω 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 ζ = R Z ζ . 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 α. -The zone Z 0 is represented by one subdomain, say Ω i 0 , and the mesh T We construct the meshes T (Ω i ζ ) (we will explain later how we can choose the grid size) in order to satisfy the following properties: for . Thus, we have the following relations: We need to specify the mesh size for every T (Ω i ζ ) in order to satisfy inequalities (3.1). We set the mesh size of T , and, for a uniform subdomain mesh, we can set µ and int(ν −µ ) denotes the nearest integer to ν −µ . Notice that the grading has "a subdomain character" and is mainly determined by the parameter µ ∈ (0, 1]. For µ = 1, we get h i ζ = h, i.e., means we get quasi-uniform meshes. Using inequality µ ≤ 1 and inequality (a + b) γ ≤ 2 γ−1 (a γ + b γ ), which can easily be shown since the function t γ is convex in (0, ∞), we arrive at the estimates which gives the right inequality in (3.1b). By the initial choice of h i ζ , we with C m = 1 2 ((n ζ + ν)

Remark 2.
It is possible to apply other techniques, see for example [6,5], of constructing graded meshes, where we could prove optimal rates for the dG IgA method. We prefer the way that is described above for its simplicity and because it suits to the spirit of the dG IgA methodology.

Remark 3.
If h i and h j are the grid sizes of two adjacent subdomains Ω i and Ω j , then relations (3.1) immediately yield the two-side estimate σ m,ζ ≤ h i /h j ≤ σ M,ζ , where the positive constants σ m,ζ and σ M,ζ only depend on the quantities which specify the initial zone partition Z ζ . Hence, in what follows, estimate (3.7a) will be used in the form . We mention that the analysis presented here can easily be extended to non-matching grids, see [24]. We also note that the meshes T (Ω i ζ ) satisfy the Assumption 1.
We emphasize that Lemma 1 provides local estimates which hold in every subdomain Ω i . This help us to investigate the accuracy of the method in every zone Z ζ of U s separately. We give an approximation estimate for the case where u r ∈ W l,2 (Ω) with l ≥ k + 1, see (2.6).
Theorem 1. Let Z ζ be a zone partition of Ω with the properties listed in the previous section, and let T (i) h i (Ω i ) be the meshes of the subdomains as described in Section 3.1. Then, for the solution u of (2.4a), we have the error estimate u − u h dG(Ω) ≤ h r C, with r = min{k, λ/µ}, (3.9) where the constant C > 0 is determined by the quasi uniform mesh properties, see (3.1), and the constants C i of Lemma 1.

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. The first examples concern 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 .

Implementation details
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. In Fig. 3(a), we show the one-dimensional B-spline basis on a graded mesh, and similarly in In Fig. 3(b), we present the two-dimensional B-spline basis 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 assumption of a well-behaved B-spline geometry map, the parameter mesh is transformed from T in such a 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 T H (Ω) is created by the grading function with the same number of knots as an equivalent uniform mesh with grid size h i for each subdomain Ω i , but pulled towards the singularity P s using the grading parameter µ. This strategy satisfies Assumption 1 with a minimal number of knots. This approach also mimics the zones construction introduced in Subsection 3.1, since it corresponds to a zone partition that shrinks towards P s at every refinement step.
In our experiments, we consider a mapping Φ i produced on an initial knot vector Ξ d i , which exactly represents the subdomain Ω i , as the isogeometric paradigm suggests. The knots are relocated during the grading procedure but without changing the shape or the parameterization of the subdomains Ω i . If needed, the original coarse knots are inserted in the discretization basis such that the exact representation of the original shape is feasible. Nevertheless, in our implementation, we have the freedom to use a different sequence of knots for the discretization space without refining the initial geometry in this basis.   Table 2. Heart shape problem: The convergence rates in the dG-norm . dg(Ω) with and without grading.
The computational domain Ω consists of two subdomains shown in Fig. 4(a), where the corresponding knot vectors are , 0, 0, 1, 1, 1} and are parametrized by k = 2 Bspline basis with the control points given in Table 1. The exact solution is given by u = r π ω sin(θπ/ω) with f and u D in (2.4) are specified by the exact solution. We set α = 1 in the entire Ω. Note that u ∈ W 1.5,2 (Ω). The problem has been solved using first (k = 1) and second (k = 2) order Bspline spaces with grading parameter µ = 0.6 and µ = 0.3, respectively, see Fig. 4(b). We plot the contours of the solution u h computed using k = 2 B-splines in Figure 4(c). In Table 2, 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. One 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.

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 method can not provide an (optimal) accurate approximation [12]. We examine such a case by solving the so-called Kellogg test problem [18]. We consider the computational domain Ω = (−1, 1) 2 . The diffusion coefficient α in (2.3) is supposed to be piecewise constant taking the same value, say α := α 13 , in the first and third quadrants, and, similarly, for the second and fourth quadrants, α := α 24 . 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 where the numbers λ, ρ, σ satisfy the nonlinear relations For λ = 0.4, the solution u ∈ W 1.4,2 (Ω), and has discontinuous derivatives across the interfaces. On the other hand, u ∈ W 2,1.25 (Ω), and the estimates presented in Section 3.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.9). In Fig. 5(b), we can see the graded meshes of the subdomains. Fig. 5(c) shows the plot of the contours of the dG solution u h computed for degree k = 1 B-splines. In Table 3, we display the convergence rates of the solution. We observe that, in the case of uniform meshes, the experimental order of convergence of the method is 0.4 which is determined by the regularity of the solution. Conversely, the rates in the right columns which correspond to the results using mesh grading tend to be optimal with respect the order of the B-spline space. A glimpse of the discrete solution u h is given in Fig. 6.  Table 3. Kellogg's test: The convergence rates in the dG-norm . dg(Ω) .

Three dimensional examples
4.3.1 Cube with interior point singularity. This test case is inspired by [21]. The computational domain is Ω = (−1, 1) 3 which is decom- posed into 8 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 u h computed by k = 1 B-spline space. The convergence rates of the error corresponding to non-graded meshes are shown in left columns of Table 4. The rates are optimal for k = 1 B-spline space and sub-optimal for the two other B-spline spaces, as it was expected according to the regularity of the solution u. Note that the rates presented on the left columns in Table  4 are in agreement with the estimate given in (3.7a). We have performed again the test using grading 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 1 and (3.9). In Fig. 7(c), the variate of the u h contours around the singular point is shown. The rates obtained on graded meshes are displayed in right columns of Table 4.  Table 4. Cube with interior singularity: The convergence rate of the error on uniform and graded meshes.
We can observe that the rates approach the optimal rate for both high-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 3d L-shape form 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 [8]), we believe that it is an interesting test case, because we will see that the graded 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 u D 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 u h computed for k = 1 are pre- sented in Fig. 8(c). Table 5 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 the same type of graded meshes have also been used in finite element methods for approximating solutions of elliptic problems in three-dimensional domains with edges, see [6,5,4].  Table 5. 3d L-shape : The convergence rates of the error with respect to the dG norm on uniform and graded meshes.

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 u D 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 u h computed with degree k = 1 is presented in Fig. 9(c).
The convergence rates of the error corresponding to the quasi-uniform meshes are shown in left columns (without grading) of Table 6, and the rates corresponding to the graded meshes are shown in the right columns of Table 6.  Table 6. 3d Heart : The convergence rates of the error with respect to the dG norm on uniform and graded meshes.

Conclusion
We have presented mesh grading techniques for dG IgA discretizions of elliptic boundary value problems in the presence of so-called singular points. Based on the a priori or a posteriori knowledge of the behaviour of the exact solution around the singular points, we pre-defined the grading of the mesh without increasing the knots but performing a relocation. The grading refinement has a subdomain (patch) character in order to fit well into the IgA framework. Optimal error estimates of the multipatch dG IgA method have been shown when it is used on the graded meshes proposed. The theoretical results have been confirmed by a number of two-and three-dimensional test problems with known exact solutions.