A quasi-robust discretization error estimate for discontinuous Galerkin Isogeometric Analysis

Isogeometric Analysis is a spline-based discretization method to partial differential equations which shows 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 (SIPG) method for Isogeometric Analysis has been previously introduced. In the present paper, we give error estimates that only depend poly-logarithmically on 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.

Authors

• 13 publications
02/05/2019

Fast multigrid solvers for conforming and non-conforming multi-patch Isogeometric Analysis

Isogeometric Analysis is a high-order discretization method for boundary...
05/19/2020

Convergence theory for IETI-DP solvers for discontinuous Galerkin Isogeometric Analysis that is explicit in h and p

In this paper, we develop a convergence theory for Dual-Primal Isogeomet...
12/12/2018

Efficient discontinuous Galerkin implementations and preconditioners for implicit unsteady compressible flow simulations

This work presents and compares efficient implementations of high-order ...
09/03/2020

A two level method for isogeometric discretizations

Isogeometric Analysis (IGA) is a computational technique for the numeric...
07/23/2020

Linearizing the hybridizable discontinuous Galerkin method: A linearly scaling operator

This paper proposes a matrix-free residual evaluation technique for the ...
12/01/2018

Recovering missing CFD data for high-order discretizations using deep neural networks and dynamics learning

Data I/O poses a significant bottleneck in large-scale CFD simulations; ...
06/11/2020

Frontiers in Mortar Methods for Isogeometric Analysis

Complex geometries as common in industrial applications consist of multi...
This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1 Introduction

The original design goal of Isogeometric Analysis (IgA), Hughes:2005

, 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 (PDE) 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. Kapl:Sangalli:Takacs:2018 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 are discontinuous Galerkin approaches, cf. Riviere:2008 ; Arnold:Brezzi:Cockburn:Marini:2002 , particularly the symmetric interior penalty discontinuous Galerkin (SIPG) method Arnold:1982 . The idea of applying this technique to couple patches in IgA, has been previously discussed in LMMT:2015 ; LT:2015 .

Concerning the approximation error, in early IgA literature, only its dependence on the grid size has been studied, cf. Hughes:2005 ; Bazilevs:BeiraoDaVeiga:Cottrell:Hughes:Sangalli:2006 . In recent publications BeiraoDaVeiga:Buffa:Rivas:Sangalli:2011 ; Takacs:Takacs:2015 ; Floater:Sande:2017 ; Sande:Manni:Speleers:2018 also the dependence on the spline degree has been investigated. These error estimates are restricted to the single-patch case. In Takacs:2018 , the results from Takacs:Takacs:2015 on approximation errors for B-splines of maximum smoothness have been extended to the conforming multi-patch case.

For the case of discontinuous Galerkin discretizations, only error estimates in the grid size are known, cf. LT:2015 . The goal of the present paper, is to present an error analysis in the grid size , the spline degree and patchwise constant diffusion coefficients . We observe that under reasonable assumptions, the approximation error drops like , while being robust in the coefficients . The dependence on the spline degree is only poly-logarithmically, cf. (15). This might be surprising as the penalization parameter has to grow like for the SIPG method to be well-posed.

The robust error estimates presented of this paper can be used to analyze multigrid solvers for discontinuous Galerkin multipatch discretizations, see Takacs:2019a for solvers that show robust convergence behavior in numerical experiments.

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. We provide numerical experiments that depict our estimates, in Section 5.

2 The model problem and its discretization

We consider the following Poisson model problem. Let be an open and simply connected Lipschitz domain. For any given source function , we are interested in the function solving

 (α∇u,∇v)L2(Ω)=(f,v)L2(Ω)for% all v∈H1,∘(Ω), (1)

where is piecewise constant. Here and in what follows, for any , and are the standard Lebesgue and Sobolev spaces with standard scalar products , , norms and , and seminorms . The Lebesgue space of function with zero mean is given by .

The computational domain is the union of non-overlapping open patches , i.e.,

 ¯¯¯¯Ω=K⋃k=1¯¯¯¯¯¯ΩkandΩk∩Ωl=∅ for any k≠l (2)

holds, where denotes the closure of . We assume that the patches are constructed such that the coefficient function is constant on each patch, i.e.,

 α=αkonΩk,whereαk∈R+.

Each patch is represented by a bijective geometry function

 Gk:ˆΩ:=(0,1)2→Ωk:=Gk(ˆΩ)⊂R2,

which can be continuously extended to the closure of such that . We use the notation

 vk:=v|Ωkandˆvk:=vk∘Gk

for any function on . If , we can use standard trace theorems to extend to and to extend 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 1

For any two patches and with , the intersection is either (a) empty, (b) a common vertex, or (c) a common edge such that

 ˆIk,l:=Ik,l∘Gk∈ˆI:={{0}×(0,1),{1}×(0,1),(0,1)×{0},(0,1)×{1}}. (3)

Note that the pre-images and do not necessarily agree. We define

 N(k) :={l∈{1,…,K}:Ωk and Ωl% have common edge}, N :={(k,l)∈{1,…,K}2:kl and l∈N(k)},

and the parameterization via

 γk,l(t):={(t,s)ifˆIk,l=(0,1)×{s},s∈{0,1}(s,t)ifˆIk,l={s}×(0,1),s∈{0,1}.. (4)

We assume that the geometry functions agree on the interface; this does not require any smoothness of the overall geometry function normal to the interface.

Assumption 2

For all and , we have

 γk,l(t)=G−1k∘Gl∘γl,k(t)orγk,l(t)=G−1k∘Gl∘γl,k(1−t).

We can reparameterize each patch such that this condition is satisfied. Assume to have two patches and , sharing the patch . Using

 ˜Gk(x,y):=Gk(yx+(1−y)ρ(x),y),where(ρ(t),0):=G−1k∘Gl(t,0),

we obtain a reparameterization of , which (a) matches the parameterization of at the interface, (b) is unchanged on the other interfaces, and (c) keeps the patch unchanged. By iteratively applying this approach to all patches, we obtain a discretization satisfying Assumption 2.

We assume that the geometry function is sufficiently smooth such that the following assumption holds.

Assumption 3

There is a constant such that the geometry functions satisfy the estimate

 supx∈ˆΩ∥∇rGk(x)∥ℓ2≤CGandsupx∈ˆΩ∥(∇rGk(x))−1∥ℓ2≤CGforr∈{1,2}.

We assume full elliptic regularity.

Assumption 4

The solution of the model problem (1) is patch-wise , i.e.,

 uk∈H2(Ωk)

holds for all .

If all are equal, we obtain (and thus also Assumption 4) for domains with a sufficiently smooth boundary, cf. Necas:1967 , and for convex polygonal domains , cf. Dauge:1988 ; Dauge:1992 . This case is of interest if the dG discretization is used to obtain a flexible combination of the patches. If not all values of agree, in general , but Assumption 4 might be satisfied under certain circumstances, cf. Petzoldt:2001 ; Petzoldt:2002 and others. The theory of this paper can be extended to cases where we only know for some . For simplicity, we restrict ourselves to the case of full elliptic regularity (Assumption 4).

Having a representation of the domain, we introduce the isogeometric function space. Following LMMT:2015 ; LT:2015 , we use a conforming isogeometric discretization for each patch and couple the contributions for the patches using a symmetric interior penalty discontinuous Galerkin (SIPG) method, cf. Arnold:1982 , as follows.

For the univariate case, the space of spline functions of degree and size with is given by

 Sp,h(0,1):={v∈Hp(0,1):v|(ih,(i+1)h]∈Pp% for all j=1,…,n−1},

where is the space of polynomials of degree . On the parameter domain , we introduce tensor-product B-spline functions

 Sp,h(ˆΩ):=Sp,h(0,1)⊗Sp,h(0,1).

The multi-patch function space is given by

 Vh:={uh∈L∘2(Ω):uh∘Gk∈Spk,hk(ˆΩ) for k=1,…,K,}. (5)

Note that the grid sizes and the spline degrees can be different for each of the patches. We define

 p:=maxk∈{1,…,K}pk,pmin:=mink∈{1,…,K}pk,andh:=maxk∈{1,…,K}hk

to be the largest spline degree, the smallest spline degree and the largest grid size and assume there to be a constant such that for all .

Following the assumption that is a patchwise function, we define for each a broken Sobolev space

 Hr(Ω):={v∈L2(Ω):vk∈Hr(Ω)},

with associated norms and scalar products

 ∥v∥Hr(Ω):=(v,v)1/2Hr(Ω)and(u,v)Hr(Ω):=K∑k=1(u,v)Hr(Ωk)

and weighted norms and scalar products

 ∥v∥Hrα(Ω):=(v,v)1/2Hrα(Ω)and(u,v)Hrα(Ω):=K∑k=1αk(u,v)Hr(Ωk).

For each patch, we define on its boundary

the outer normal vector

. On each interface , we define the jump operator by

 ⟦v⟧:=vk−vlon Ik,l=Il,k where% (k,l)∈N

and the average operator by

 {v}:=12(vk+vl)on Ik,l=Il,k where (k,l)∈N.

The discretization of the variational problem using the symmetric interior penalty discontinuous Galerkin method reads as follows. Find such that

 (uh,vh)Ah=(f,vh)L2(Ω)for all vh∈Vh, (6)

where

 (u,v)Ah :=(u,v)H1α(Ω)−(u,v)Bh−(v,u)Bh+(u,v)Ch, (u,v)Bh :=∑(k,l)∈N(⟦u⟧,{α∇v}⋅{n}k)L2(Ik,l), (u,v)Ch :=σh∑(k,l)∈Nαk,l(⟦u⟧,⟦v⟧)L2(Ik,l)

for all ,

 αk,l:=max{αk,αl}

and the penalty parameter is chosen sufficiently large.

Using a basis for the space , we obtain a standard matrix-vector problem: Find such that

 Ahu––h=f––h. (7)

Here and in what follows, is the coefficient vector representing with respect to the chosen basis, i.e., , and 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 , and . Before we proceed, we introduce a convenient notation.

Definition 5

Any generic constant used within this paper is understood to be independent of the grid size , the spline degree and the number of patches , but it might depend on the constants , and .

We use the notation if there is a generic constant such that and the notation if and .

For symmetric positive definite matrices and , we write

 A≤Bifv–⊤hAv–h≤v–⊤hBv–h for all vectors v–h.

The notations and are defined analogously.

3 A discretization error estimate

In LMMT:2015 , it has been shown that the bilinear form 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 . Thus, we define the dG-norm via

 ∥u∥2Qh:=(u,u)Qh,where(u,v)Qh:=(u,v)H1α(Ω)+(u,v)Ch

for all . Note that we define the norm differently to LMMT:2015 , where the dG-norm was independent of .

Before we proceed, we give some estimates on the geometry functions.

Lemma 6

The geometry functions satisfy

 ∥v∘G−1k∥Hr(Ωk)≂∥v∥Hr(ˆΩ) for all v∈Hr(ˆΩ),r∈{0,1,2},and ∥v∘G−1k∥L2(Ik,l)≂∥v∥L2(ˆIk,l) for all v∈H1(ˆΩ).

For ease of notation, here and in what follows, we define .

The statements follow directly from the chain rule for differentiation, the substitution rule for integration and Assumption

3. ∎

Lemma 7

The geometry functions satisfy

 ∥(∇v∘G−1k)⋅{n}k∥L2(Ik,l)≲∥∇v∥L2(ˆIk,l) for all v∈H2(ˆΩ).

We have

 ∥(∇v∘G−1k)⋅{n}k∥L2(Ik,l)≤∥∇v∘G−1k∥L2(Ik,l)∥{n}k∥L∞(Ik,l),

where certainly because the length of is always . The estimate follows directly from the chain rule for differentiation, the substitution rule for integration and Assumption 3. ∎

For sufficiently large, the symmetric bilinear form is coercive and bounded, i.e., a scalar product.

Theorem 8 (Coercivity and boundedness)

There is some that only depends on and such that

 (uh,uh)Ah≳∥uh∥2Qhand(uh,vh)Ah≲∥uh∥Qh∥v∥Qh

holds for all and all .

Note that Using Lemma 7, (Takacs:2018, , Lemma 4.4), (Schwab:1998, , Corollary 3.94) and Lemma 6, we obtain

 ∥∇vh⋅{n}k∥2L2(Ik,l)≲∥∇(vh∘Gk)⋅{n}k∥2L2(ˆIk,l)≲∥vh∘Gk∥H2(ˆΩ)∥vh∘Gk∥H1(ˆΩ) ≲p2h∥vh∘Gk∥2H1(ˆΩ)≲p2h∥vh∥2H1(Ωk) (8)

for all , and . As , the Poincaré inequality (see, e.g., (Schwab:1998, , Theorem A.25)) yields also

 ∥∇vh⋅{n}k∥2L2(Ik,l)≲p2h|vh|2H1(Ωk).

The Cauchy-Schwarz inequality, the triangle inequality, (8), , and yield

 |(uh,vh)Bh| (9) ≲⎛⎝∑(k,l)∈Nαk,l∥⟦uh⟧∥2L2(Ik,l)⎞⎠1/2⎛⎝K∑k=1∑j∈N(k)α−1k,lα2k∥∇vh⋅{n}k∥2L2(Ik,l)⎞⎠1/2 ≲(p2h)1/2⎛⎝∑(k,l)∈Nαk,l∥⟦uh⟧∥2L2(Ik,l)⎞⎠1/2(K∑k=1αk|vh|2H1(Ωk))1/2 ≤pσ−1/2∥uh∥Qh∥vh∥Qh

for all . Let be the hidden constant, i.e., such that

 |(uh,vh)Bh|≤c0pσ−1/2∥uh∥Qh∥vh∥Qh. (10)

For , we obtain

 (uh,uh)Ah=∥uh∥2Qh−2(uh,uh)Bh≥12∥uh∥2Qh,

i.e., coercivity. Using (9) and Cauchy-Schwarz inequality, we obtain further

 (uh,vh)Ah=(uh,vh)Qh−(uh,vh)Bh−(vh,uh)Bh≤32∥uh∥Qh∥vh∥Qh,

i.e., boundedness. ∎

As we have boundedness and coercivity (Theorem 8), the Lax Milgram theorem (see, e.g., (Schwab:1998, , Theorem 1.24)) yields states existence and uniqueness of a solution, i.e., the following statement.

Theorem 9 (Existence and uniqueness)

If is chosen as in Theorem 8, the problem (6) has exactly one solution .

The following theorem shows that the solution of the original problem also satisfies the discretized bilinear form.

Theorem 10 (Consistency)

The solution of the original problem (1) satisfies

 (u,vh)Ah=(f,vh)L2(Ω)for all vh∈Vh.

For a proof, see, e.g., (Riviere:2008, , Proposition 2.9); the proof requires elliptic regularity (cf. Assumption 4).

If boundedness of the bilinear form was also satisfied for , Ceá’s Lemma (see, e.g., (Schwab:1998, , Theorem 2.19.iii)) would allow to bound the discretization error. However, the bilinear form is not bounded in the norm , but only in the stronger norm , given by

 (11)
Theorem 11

There is some that only depends on and such that

 (u,vh)Ah≲∥u∥Q+h∥vh∥Qh

holds for all , all all .

Let and be arbitrarily but fixed and assume , where is as in (10). Note that the arguments from (9) also hold if the first parameter of the bilinear form is not in . So, we obtain

 |(u,vh)Bh|≤14∥u∥Qh∥vh∥Qh.

Using Lemma 7, (Takacs:2018, , Lemma 4.4), Lemma 6 and the Poincaré inequality, we obtain

 ∥∇v⋅{n}k∥2L2(Ik,l)≲∥∇(v∘Gk)⋅{n}k∥2L2(ˆIk,l)≲∥v∘Gk∥H2(ˆΩ)∥v∘Gk∥H1(ˆΩ) ≲∥v∥H2(Ωk)∥v∥H1(Ωk)≤1β∥v∥2H2(Ωk)+β∥v∥2H1(Ωk)≤1β|v|2H2(Ωk)+β|v|2H1(Ωk) (12)

for all , all , all and all . Using this estimate, , and , we obtain for

 |(vh,u)Bh| ≤⎛⎝σh∑(k,l)∈Nαk,l∥⟦vh⟧∥2L2(Ik,l)⎞⎠1/2⎛⎝hσK∑k=1∑l∈N(k)α2kαk,l∥∇u⋅{n}k∥2L2(Ik,l)⎞⎠1/2 ≤∥vh∥Qh(K∑k=1αk|u|2H1(Ωk)+h2σ2K∑k=1αk|u|2H2(Ωk))1/2 ≲∥vh∥Qh∥u∥Q+h.

Using these estimates, we obtain

 (u,vh)Ah =(u,vh)Qh−(u,vh)Bh−(vh,u)Bh≲∥u∥Q+h∥vh∥Qh,

which finishes the proof. ∎

Using consistency (Theorem 10), coercivity and boundedness (Theorems 8 and 11), we can bound the discretization error using a the approximation error.

Theorem 12 (Discretization error estimate)

Provided the assumptions of Theorems 8 and 10, the estimate

 ∥u−uh∥Qh≲infvh∈Vh∥u−vh∥Q+h

holds, where is the solution of the original problem (1) and is the solution of the discrete problem (6).

For any , the triangle inequality yields

 ∥u−uh∥Qh≤∥u−vh∥Qh+∥vh−uh∥Qh. (13)

Theorem 10 and Galerkin orthogonality yield for all . So, we obtain using Theorems 8 and 11 that

 ∥vh−uh∥2Qh≲(vh−uh,vh−uh)Ah=(vh−u,vh−uh)Ah≲∥vh−u∥Q+h∥vh−uh∥Qh,

which shows . Together with (13), this shows . Since this holds for all , this finishes the proof. ∎

It is rather straight forward to give approximation error estimates that bound the approximation error as follows:

 infvh∈Vh∥u−vh∥Q+h≲σ1/2h|u|H2(Ω)

for all . If is chosen in an optimal way, this yields a result of the form

 infvh∈Vh∥u−vh∥Q+h≲σ1/20p2h|u|H2(Ω),

i.e., a quadratic increase in the spline degree . Using a refined analysis, we obtain as follows.

Theorem 13 (Approximation error estimate)

Provided , and , the estimate

 infvh∈Vh∥u−vh∥Q+h≲(lnσ)2σ1/(2pmin−1)h|u|H2α(Ω) (14)

holds for all .

The proof is given at the end of the next section.

If we consider the case and if we do not consider over-penalization, i.e., we assume , we obtain using Theorem 12 the estimate

 ∥u−uh∥Q+h≲infvh∈Vh∥u−vh∥Q+h≲(lnσ)2h|u|H2α(Ω), (15)

where is the solution of the original problem and is the solution of the problem discretized with the proposed SIPG approach.

4 Proof of the approximation error estimate

Before we can give the proof, we give some auxiliary results. This section is organized as follows. In Section 4.1, we give patch-wise 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.

4.1 Patch-wise projectors

As first step, we recall the projection operators from (Takacs:2018, , Sections 3.1 and 3.2). Let be the -orthogonal projection into , where

 (u,v)H1D(0,1)=(u′,v′)L2(0,1)+u(0)v(0).

(Takacs:2018, , Lemma 3.1) states that and . Using , we obtain for and that

 (u−Πp,hu,1)L2(0,1)=0. (16)

The next step is to consider the multivariate case, more precisely the parameter domain . Let and be given by

 (Πxp,hu)(x,y)=(Πp,hu(⋅,y))(x)and(Πyp,hu)(x,y)=(Πp,hu(x,⋅))(y)

and let be such that

 ˆΠk=Πxpk,hkΠypk,hk. (17)

For the physical domain, define to be such that

 (Πv)|Ωk=(ˆΠk(v∘Gk))∘G−1kfor all v∈H2(Ω) and k=1,…,K.

Observe that we obtain using (16) that

 (u−ˆΠku,1)L2(ˆΩ)=0,ˆΠkc=candΠc=c. (18)

for all .

The projectors satisfy robust error estimates and are almost stable in .

Lemma 14

holds for all .

This result follows directly from (Takacs:2018, , Theorem 3.3).

Lemma 15

holds for all .

The proof is analogous to the proof of (Hofreither:Takacs:2017, , Theorem 4). Let