Discretization error estimates for discontinuous Galerkin isogeometric analysis

Isogeometric analysis is a spline-based discretization method to partial differential equations which show the approximation power of a high-order method. The number of degrees of freedom, however, is as small as the number of degrees of freedom of a low-order method. This does not come for free as the original formulation of isogeometric analysis requires a global geometry function. Since this is too restrictive for many kinds of applications, the domain is usually decomposed into patches, where each patch is parameterized with its own geometry function. In simpler cases, the patches can be combined in a conforming way. However, for non-matching discretizations or for varying coefficients, a non-conforming discretization is desired. An symmetric interior penalty discontinuous Galerkin method for isogeometric analysis has been previously introduced. In the present paper, we give error estimates that are explicit in the spline degree. This opens the door towards the construction and the analysis of fast linear solvers, particularly multigrid solvers for non-conforming multipatch isogeometric analysis.


Introduction
The original design goal of isogeometric analysis (IgA) [1] was to unite the world of computer aided design (CAD) and the world of finite element (FEM) simulation. In IgA, both the computational domain and the solution of the partial differential equation are represented by spline functions, like tensor-product B-splines or non-uniform rational B-splines (NURBS). This follows the design goal since such spline functions are also used in standard CAD systems to represent the geometric objects of interest.
The parameterization of the computational domain using just one tensor-product spline function is possible only in simple cases. A necessary condition for this to be possible is that the computational domain is topologically equivalent to the unit square or the unit cube. This might not be the case for more complicated computational domains. Such domains are typically decomposed into subdomains, in IgA called patches, where each of them is parameterized by its own geometry function. The standard approach is to set up a conforming discretization. For a standard Poisson problem, this means that the overall discretization needs to be continuous. For higher order problems, like the biharmonic problem, even more regularity is required, conforming discretizations in this case are rather hard to construct, cf. Ref. [2] and references therein.
Even for the Poisson problem, a conforming discretization requires the discretizations to agree on the interfaces. This excludes many cases of practical interest, like having different grid sizes or different spline degrees on the patches. Since such cases might be of interest, alternatives to conforming discretizations are of interest. One promising alternative is discontinuous Galerkin approaches, cf. Refs. [3,4], particularly the symmetric interior penalty discontinuous Galerkin (SIPG) method [5]. The idea of applying this technique to couple patches in IgA has been previously discussed in Refs. [6][7][8].
Concerning the approximation error, in early IgA literature, only its dependence on the grid size has been studied, cf. Refs. [1,9]. In recent publications [10][11][12][13] also the dependence on the spline degree has been investigated. These error estimates are restricted to the single-patch case. In Ref. [14], the results from Ref. [13] on approximation errors for B-splines of maximum smoothness have been extended to the conforming multipatch case.
For the case of discontinuous Galerkin discretizations, only error estimates in the grid size are known, cf. Ref. [8]. The goal of the present paper is to present an analysis that is explicit both in the grid size h and the spline degree p. We show that the penalty parameter has to grow like p 2 for the SIPG method to be well posed. If the solution is sufficiently smooth, the error in the energy norm decays like h p . For the analysis of linear solvers, we need also estimates for less smooth functions. We give a bound that degrades only poly-logarithmically with p, cf. (21). These error estimates have recently been used to analyze multigrid solvers for discontinuous Galerkin multipatch discretizations, see Ref. [15]. They might also be helpful for the analysis of Finite Element Tearing and Interconnecting (FETI) solvers, cf. Refs. [6,16], Balancing Domain Decomposition by Constraints (BDDC) solvers, cf. Refs. [17][18][19] and references therein, Schwarz type solvers, cf. Ref. [20] and references therein, and similar methods.
The remainder of the paper is organized as follows. In Section 2, we introduce the model problem and give a detailed description of its discretization. A discussion of the existence of a unique solution and the discretization and the approximation error is provided in Section 3. The proof of the approximation error estimate is given in Section 4. In Section 5, we provide numerical experiments that depict our estimates. In Section 6, we summarize the findings and discuss possible extensions.

The model problem and its discretization
We consider the following Poisson model problem. Let ⊂ R 2 be an open and simply connected Lipschitz domain. For any given source function f ∈ L 2 ( ), we are interested in the function u ∈ Here and in what follows, for any r ∈ N := {1, 2, 3, . . .}, L 2 ( ) and H r ( ) are the standard Lebesgue and Sobolev spaces with standard scalar products (·, ·) L 2 ( ) , (·, ·) H r ( ) := (∇ r ·, ∇ r ·) L 2 ( ) , norms · L 2 ( ) and · H r ( ) , and seminorms | · | H r ( ) . The Lebesgue space of function with zero mean is given by holds, where T denotes the closure of T. Each patch k is represented by a bijective geometry function which can be continuously extended to the closure of such that G k ( ) = k . We use the notation v k := v| k and v k := v k • G k for any function v on . If v ∈ H 1 ( ), we can use standard trace theorems to extend v k to k and to extend v k to .
We assume that the mesh induced by the interfaces between the patches does not have any T-junctions, i.e. we assume as follows. Assumption 2.1: For any two patches k and l with k = l, the intersection ∂ k ∩ ∂ l is either (a) empty, (b) a common vertex or (c) a common edge I k,l = I l,k such that Note that the pre-images I k,l and I l,k do not necessarily agree. We define We assume that the geometry functions agree on the interface (up to the orientation); this does not require any smoothness of the overall geometry function normal to the interface.

Remark 2.1:
For any domain satisfying Assumption 2.1, we can reparameterize each patch such that this condition is satisfied. Assume to have two patches k and l , sharing the patch we obtain a reparameterization of G k , which (a) matches the parameterization of l at the interface, (b) is unchanged on the other interfaces, and (c) keeps the patch k unchanged. By iteratively applying this approach to all patches, we obtain a discretization satisfying Assumption 2.2.
We assume that the geometry function is sufficiently smooth such that the following assumption holds.

Assumption 2.3:
There is a constant C G > 0 such that the geometry functions G k satisfy the estimates for r = 1, 2.
We assume full elliptic regularity.
For domains with a sufficiently smooth boundary, cf. Ref. [21], and for convex polygonal domains , cf. Refs. [22,23], we have u ∈ H 2 ( ) and thus also Assumption 2.4. In case of varying diffusion conditions (which are uniform on each patch), we might have u ∈ H 2 ( ), but Assumption 2.4 might still be satisfied, cf. Refs. [24,25] and others. The theory of this paper can be extended to cases where we only know u k ∈ H 3/2+ ( k ) for some > 0. For simplicity, we restrict ourselves to the case of full elliptic regularity, i.e. Assumption 2.4.
Having a representation of the domain, we introduce the isogeometric function space. Following Refs. [7,8], we use a conforming isogeometric discretization for each patch and couple the contributions for the patches using a SIPG method, cf. Ref. [5], as follows.
On the parameter domain := (0, 1) 2 , we introduce tensor-product B-spline functions The multipatch function space V h is given by Note that the grid sizes h k and the spline degrees p k can be different for each of the patches. We define for some constant C h > 0, where h k,δ,min refers to the smallest knot span. Following the assumption that u h is a patchwise function, we define for each r ∈ N, a broken Sobolev space with associated norms and scalar products For each patch, we define on its boundary ∂ k , the outer normal vector n k . On each interface I k,l , we define the jump operator [ The discretization of the variational problem using the SIPG method reads as follows. Find u h ∈ V h such that where is chosen sufficiently large. Using a basis for the space V h , we obtain a standard matrix-vector problem: find u h ∈ R N such that Here and in what follows, is the coefficient vector obtained by testing the right-hand-side functional with the basis functions.
As the dependence on the geometry function is not in the focus of this paper, unspecified constants might depend on C G , C I , and C h . Before we proceed, we introduce a convenient notation.

Definition 2.5:
Any generic constant c > 0 used within this paper is understood to be independent of the grid size h, the spline degree p, and the number of patches K, but it might depend on the constants C G , C I , and C h .
We use the notation a b if there is a generic constant c > 0 such that a ≤ cb and the notation a b if a b and b a.
For symmetric positive definite matrices A and B, we write The notations A B and A B are defined analogously.

A discretization error estimate
In Ref. [7], it has been shown that the bilinear form (·, ·) A h is coercive and bounded in the dG-norm. For our further analysis, it is vital to know these conditions to be satisfied with constants that are independent of the spline degree p. Thus, we define the dG-norm via for all u, v ∈ H 2,• ( ). Note that we define the norm differently to Ref. [7], where the dG-norm was independent of p.
Before we proceed, we give some estimates on the geometry functions.
where r = 0, 1, 2. For ease of notation, here and in what follows, we define H 0 := L 2 . If (5) holds for r = 1, 2, . . . , s with some s > r, then (12) also holds for those choices of r. Moreover, we have Proof: The statements follow directly from the chain rule for differentiation, the substitution rule for integration, and Assumption 2.3.

Lemma 3.2: We have
∇v L 2 ( I k,l ) follows directly from the chain rule for differentiation, the substitution rule for integration, and Assumption 2.3.
For σ sufficiently large, the symmetric bilinear form (·, ·) A h is coercive and bounded, i.e. a scalar product.

Theorem 3.3 (Coercivity and boundedness):
There is some σ 0 > 0 that only depends on C G , C I , and C h such that for all v h ∈ V h , k = 1, . . . , K, and l ∈ N (k).
The Cauchy-Schwarz inequality, the triangle inequality, (13), and |N (k)| ≤ 4 yield for all u h , v h ∈ V h . Let c 0 1 be the hidden constant, i.e. such that For σ ≥ 16 c 0 p 2 , we obtain i.e. coercivity. Using (14) and the Cauchy-Schwarz inequality, we obtain further i.e. boundedness.  The following theorem shows that the solution of the original problem also satisfies the discretized bilinear form.
If boundedness of the bilinear form (·, ·) A h was also satisfied for u ∈ H 2,• ( ), Ceá's Lemma (see, e.g. Ref. [26,Theorem 2.19.iii]) would allow to bound the discretization error. However, the bilinear form is not bounded in the norm · Q h , but only in the stronger norm · Q + h , given by Theorem 3.6: If σ is chosen as in Theorem 3.3, Proof: Let u ∈ H 2,• ( ) and v h ∈ V h be arbitrarily but fixed. Note that the arguments from (14) also hold if the first parameter of the bilinear form , and all β > 1. Using this estimate, and |N (k)| ≤ 4, we obtain for β : Using these estimates, we obtain which finishes the proof.
Using consistency (Theorem 3.5), coercivity, and boundedness (Theorems 3.3 and 3.6), we can bound the discretization error using the approximation error.
where u is the solution of the original problem (1) and u h is the solution of the discrete problem (8).
Theorem 3.5 and Galerkin orthogonality yield So, we obtain using Theorems 3.3 and 3.6 that Since this holds for all v h ∈ V h , this finishes the proof.
A proof of this theorem is given at the end of the next section. Assuming p p min , then we have for the case q = p min that For the analysis of linear solvers, like multigrid solvers [15], we also need low-order approximation error estimates. In the convergence proofs, we usually have to estimate errors of the iterative scheme.
Even if we know that the true solution satisfies certain regularity assumptions, this does not extend to the errors. For them, we can only rely on the regularity statements arising from the domain, which means that H 2 -regularity is usually the best we can hope for. For this case, we obtain This means that we obtain a quadratic increase in the spline degree p. Using a refined analysis, we obtain as follows.
The proof is given at the end of the next section. Assuming again p p min ≥ 2, we obtain i.e. an only poly-logarithmic increase in the spline degree p.

Proof of the approximation error estimates
Before we can give the proof, we give some auxiliary results. This section is organized as follows. In Section 4.1, we give patchwise projectors and estimates for them. We introduce a mollifying operator and give estimates for that operator in Section 4.2. Finally, in Section 4.3, we give the proof for the approximation error estimate.

Patchwise projectors
As first step, we recall the projection operators from Ref. [ In what follows, we also write S p,h (0, 1) and p,h if we refer to a uniform grid of size h. (1), we obtain for p ≥ 2 and u ∈ H 1 (0, 1) that The next step is to consider the multivariate case, more precisely the parameter domain = (0, 1) 2 . Let x k : H 2 ( ) → H 2 ( ) and y k : H 2 ( ) → H 2 ( ) be given by and let k : For the physical domain, define : Observe that we obtain using (22) that for all c ∈ R. The projectors k satisfy robust error estimates and are almost stable in H 2 . hold for all u ∈ H r+1 (0, 1).

Proof:
The identity (22) implies that the projector p,Z coincides with the projector Q 1 p from Ref.  1) holds for all u ∈ H r+1 (0, 1).

Proof:
The proof is analogous to the proof of Ref. [27,Theorem 4]. Let Q 2 p be the H 2 -orthogonal projector into S p,Z (0, 1) as introduced in Ref. [12]. Using the triangle inequality and a standard inverse estimate [26,Corollary 3.94], we obtain hold for all u ∈ H r+1 ( ).

Proof:
The proof is based on the univariate estimates given in Lemmas 4.1 and 4.2 (p min ≥ 2 and the quasi-uniformity of the grids have been required in (7)) and follows the standard construction that can be found, e.g. in the proof of Ref. [14,Theorem 3.3].
On the interfaces, we have the following approximation error estimate.

A mollifying operator
A second step of the proof is the introduction of a particular mollification operator for the interfaces. For (k, l) ∈ N , let ϒ k,l be given by ϒ k,l v : k (1 − t)), cf. Assumption 2.2. For all cases, ϒ k,l is a bijective function H s (0, 1) → H s ( I k,l ) and holds for all s. For v ∈ H s ( ), we define the abbreviated notation ϒ −1 k,l v := ϒ −1 k,l (v| I k,l ) and observe For (k, l) ∈ N ∪ N * , we define extension operators k,l : H s ( I k,l ) → H s ( ) by where φ(x) := max{0, 1 − η −1 x} and η ∈ (0, 1).
Before we proceed, we give a certain trace like estimate.

Proof: A trace theorem [14, Lemma 4.4] yields
Case 1. Assume θ < 1. In this case, we choose v to be the H 1 -orthogonal projection of u into S 3, θ −1 −1 ( ). Since the spline degree of that space is fixed, we obtain using a standard inverse inequality [26,Corollary 3.94] and a standard approximation error estimate (like from Ref. [13]) that Assume θ ≥ 1. In this case, we choose v := (u, 1) L 2 ( ) and obtain from (30) directly In this case, the Poincaré inequality finishes the proof.
As a next step, we show that the mollifier constructs functions that are very smooth on the interfaces.

Lemma 4.8: The estimate
Proof: We have using (25) and Lemma 4.6

Now, a standard inverse estimate [26, Corollary 3.94] yields
where ψ := 2 √ 3r 2 η −1 . Lemma 4.2 and the H 1 -stability of r,η yield | r,η w| 2 Using (25), we obtain By applying Lemma 4.7 to the derivative of u k , we obtain the desired result. Proof: Let u ∈ H 1,• ( ) ∩ H 2 ( ) be arbitrary but fixed. We obtain using the definition of ϒ k,l and ϒ l,k and Lemma 3.1 that holds, where w k := w k • G k and w l := w l • G l . Since u ∈ H 1 ( ), a standard trace theorem yields  Proof: Using the definition of M k and of the H 1,• -norm, we obtain We estimate the terms x,l , y,l , and •,l separately. Let without loss of generality I k,l = {0} × (0, 1).
The definition of w and (25) yield 2 x,l ≤ˆ 2 Equation ( x,l ≤ˆ 2 Step 2. Using (23) and the H 1 -stability of the H 1,D -orthogonal projection and w := ϒ k,l (I − r,η )ϒ −1 k,l u, we obtain 2 y,l = ∂ ∂y x k y k k,l w 2 L 2 ( ) ≤ ∂ ∂y x k k,l w 2 Using Using (25), we obtain further Case 2. Assume η 0 > 1. Define i.e. the set of all globally continuous functions which are locally just linear. Observe that W ⊆ V h . Using u and w being continuous, we obtain For the choice we further obtain using standard approximation error estimates and Lemma 3.1 which shows (38) also for the second case. Finally, we show that is such that the desired bound (20) follows. We again consider two cases.
Using the triangle inequality and Lemma 3.1, we have further where w k := (I − k ) u k . Using a trace theorem [14,Lemma 4.4] and AB ≤ γ A 2 + γ −1 B 2 , we obtain Using this, h 1, and 2 ≤ p 2 σ , we obtain The desired result is then a consequence of Lemma 4.3 and the assumed equivalence of the norms on the physical domain and the parameter domain.

Numerical experiments
We depict the results of this paper with numerical results. where u * (x, y) := sin(10πx) sin(10πy) is the exact solution. On each patch, we introduce a coarse discretization space V h 0 for = 0, which only consists of global polynomials of degree p. We then refine all grids = 1, 2, . . . times uniformly. Next, we modify the discretization spaces in order to obtain non-matching discretizations at the interfaces (since fully matching discretizations would be a special case that would allow a conforming discretization that would not be of interest in a discontinuous Galerkin setting): we refine the grid one additional time for one-third of the patches and we increase the spline degree to p + 1 for another third of the patches. In Table 1 and Figures 1 and 2, we depict the discretization errors e ,p in the H 1 -norm relative to that of the solution and the rates r ,p , given by for the case that the smoothness is C p−1 on the patches with spline degree p and C p on the patches with spline degree p + 1 (maximum smoothness). The numerical experiments show that the error decreases like h 2 − or even faster. Whether or not the error bound depends on p 2 or ln p cannot be seen in this experiment. The almost p-robust convergence of multigrid solvers whose analysis follows from the presented results can be seen in Ref. [15].   In Table 2 and Figure 3, we choose C 0 smoothness on all patches (minimum smoothness). The grid sizes and the degrees are chosen as for the first experiment. The configurations for which the computation could not be completed due to a lack of memory are indicated by 'OoM '. We obtain qualitatively the same behavior as for splines of maximum smoothness. We observe that the number of degrees of freedom in case of minimum smoothness is much larger than in the case of minimum  smoothness, which outweighs the fact that for the minimally smooth approximation, slightly smaller errors are obtained.
In Figure 4, one can observe the dependence of the discretization error on the number of patches. Here, we consider a decomposition of the initial domain into 2 ι × 2 ι patches, where ι = 1, 2, . . . , 6. We again consider a heterogeneous discretization (different spline degrees, different grid sizes), where we apply = 7 − ι uniform refinement steps. So, the grid sizes are constant. We observe that also the discretization errors are rather constant.
The discretization spaces on the patches are obtained by = 0, 1, 2, . . . uniform refinement steps, where we use splines of maximum smoothness with degree p + 1 on the patch 1 and with degree p on the remaining patches. This means that, again, for each interface, the intersection of the traces of the local function spaces of the involved spaces only contains global polynomials, so this is again a non-conforming discretization. We present the corresponding results in Table 3. As for the first experiment, we obtain the convergence rates that are slightly better than the expected ones.

Conclusions and possible extensions
In this paper, we have introduced an SIPG discretization of the Poisson equation in two dimensions which is well posed in the dG-norm, see Theorem 3.3. We have seen that the well-posedness is robust in the grid size and the spline degree. The approximation error estimate presented in Theorem 3.8 shows that the proposed method satisfies the expected approximation order. The approximation error estimate presented in Theorem 3.9 gives an error estimate in the H 2 -norm, which only grows logarithmically in the spline degree and which is of particular interest for the analysis of iterative solvers.

Remark 6.1 (Splines of reduced smoothness):
The proofs of Theorems 3.3 and 3.6 (coercivity and boundedness) do not use that the considered spline space is of maximum smoothness. Thus, Theorems 3.4 (existence and uniqueness), 3.5 ( consistency), and 3.7 (discretization error estimate) also apply if splines of reduced smoothness are used. Since the space of splines of maximum smoothness forms a subspace of the spaces of splines of reduced smoothness, the approximation error estimates given in Theorems 3.8 and 3.9 are also valid for spline spaces of reduced smoothness.

Remark 6.2 (Three-dimensional domains):
The extension of the proposed discretization technique to three-dimensional domains is completely straight-forward. Here, the interfaces I k,l are two-dimensional faces between patches, so the set N contains the indices of patches that share such faces. The proofs of Lemmas 3.1 and 3.2 and Theorems 3.3 and 3.6 (coercivity and boundedness) apply almost verbatim also to three-dimensional domains. Thus, Theorems 3.4 (existence and uniqueness), 3.5 ( consistency), and 3.7 (discretization error estimate) apply as well. The extension of the approximation error estimates to three dimensions is more involved since the extension of (23) to k := x k y k z k yields a projector that is well defined only in H 3 . Based on this observation, Theorem 3.8 (approximation error estimate) can be extended to three dimensions, provided that q ≥ 2. For the extension of Theorem 3.8 (low-order approximation error estimate), one would need to construct an H 2 -stable projectors and appropriate mollifiers, which goes beyond the scope of this paper.

Remark 6.3 (Other differential equations):
The theory presented in this paper can be extended to other elliptic differential equations whose variational formulation lives in the Sobolev space H 1 . If the differential equation is parameter-dependent, like (39), a parameter-dependent SIPG formulation can be introduced, cf. (40)-(41). For the problem (40)-(41), analogous versions to Theorems 3.3 and 3.6 ( coercivity and boundedness) can be shown for the corresponding parameter-dependent dG-norms, i.e. the norms are given by the combination of (11), (16), and (41). Thus, Theorems 3.4 (existence and uniqueness), 3.5 ( consistency), and 3.7 (discretization error estimate) apply as well. Approximation error estimates analogous to Theorems 3.8 ( approximation error estimate) and Theorem 3.8 (loworder approximation error estimate) hold, where the constants depend on max (k,l)∈N α k /α l and the H r -norms are parameter-dependent norms as defined in (41).

Disclosure statement
No potential conflict of interest was reported by the author.