BY 4.0 license Open Access Published online by De Gruyter August 7, 2021

Convergence Theory for IETI-DP Solvers for Discontinuous Galerkin Isogeometric Analysis that is Explicit in ℎ and 𝑝

Rainer Schneckenleitner and Stefan Takacs


In this paper, we develop a convergence theory for Dual-Primal Isogeometric Tearing and Interconnecting (IETI-DP) solvers for isogeometric multi-patch discretizations of the Poisson problem, where the patches are coupled using discontinuous Galerkin. The presented theory provides condition number bounds that are explicit in the grid sizes ℎ and in the spline degrees 𝑝. We give an analysis that holds for various choices for the primal degrees of freedom: vertex values, edge averages, and a combination of both. If only the vertex values or both vertex values and edge averages are taken as primal degrees of freedom, the condition number bound is the same as for the conforming case. If only the edge averages are taken, both the convergence theory and the experiments show that the condition number of the preconditioned system grows with the ratio of the grid sizes on neighboring patches.

MSC 2010: 65D07; 65F08; 65N22; 65N55

1 Introduction

Isogeometric Analysis (IgA), see [20, 9], is an approach for discretizing partial differential equations (PDEs) that has been developed in order to improve the compatibility between computer aided design (CAD) and simulation in comparison to the standard finite element method (FEM). In CAD systems, the geometry is usually parameterized in terms of geometry functions which are usually spanned by B-splines or non-uniform rational B-splines (NURBS). In IgA, we do not mesh the computational domain, but we use the geometry function to transfer the integrals that occur in the variational formulation of the boundary value problem to the parameter domain, where we set up B-splines or NURBS bases to discretize the problem. Since only simple domains can be represented by just one geometry function (single-patch case), the overall computational domain is usually decomposed into multiple patches, each of which is parameterized by its own geometry function (multi-patch IgA). We focus on non-overlapping patches.

The patches can be coupled either in a conforming way or by means of discontinuous Galerkin methods. For a conforming discretization, both the geometry function and the discretization have to agree on the interfaces between the patches. One promising alternative to overcome these restrictions are discontinuous Galerkin (dG) approaches, cf. [31, 6], where the continuity of the solution at the patch boundaries is only enforced weakly. We consider in this paper in particular the symmetric interior penalty discontinuous Galerkin (SIPG) method, cf. [5]. The idea of applying this technique to couple patches in IgA has been previously discussed in [24, 25, 38]. Such approaches also allow the handling of inexact representations of the geometry that yield gaps between or overlaps of neighboring patches, cf. [19]. For discontinuous Galerkin methods, the discretization error decays in the grid size with the same order as for a conforming discretization, cf. [38]. However, dG methods offer much greater flexibility for the setup of the geometry.

After the discretization of the PDE, we obtain a large-scale linear system, whose condition number grows exponentially in the spline degree 𝑝 and (like other FEM like methods) like h-2 with the grid size ℎ. We are interested in fast iterative solvers for such systems. Since we are in a multi-patch framework, domain decomposition solvers are a canonical choice. One of the most popular domain decomposition solvers for large-scale systems of finite element equations in standard FEM is the Finite Element Tearing and Interconnecting method (FETI), originally proposed in [14]. Since the invention of FETI, various FETI-type methods have been developed, cf. [30, 39, 23]. In [22], it was proposed to use a FETI-type method, namely the Dual-Primal FETI method (FETI-DP), in the context of IgA. This method was called Dual-Primal Isogeometric Tearing and Interconnecting (IETI-DP) method. Later, this approach has been further analyzed, particularly in [18] and more recently in [32]. The original results, cf. [18], are only explicit in the ratio between grid and patch sizes. In [32], the authors of the paper at hand have analyzed the dependence of the condition number of the preconditioned Schur complement system also on the spline degree 𝑝. In [33], the authors did a numerical study of 3D problems. The experiments indicate a similar condition number bound as established in [34]. However, the extension of the analysis to 3D problems is known to be a non-trivial task and is ongoing research.

For conforming FEM discretizations, condition number bounds that are explicit in the polynomial degree 𝑝 have been worked out previously for FETI-DP type, cf. [28, 21], other Schwarz type, cf. [36, 15], and iterative substructuring methods, cf. [2, 29].

The extension of FETI methods or IETI methods to dG discretizations is not straight-forward. We follow the approaches that have been proposed in [12] and adapted for IETI in [16], particularly the idea of using artificial interfaces. For dG discretizations, there are only a few results concerning a 𝑝-analysis for domain decomposition approaches, like [3, 35] for two-level Schwarz type approaches or [11, 8], where the estimates for BDDC and FETI-DP type methods for spectral element methods and h⁢p-FEM are given. The publication [4] shows a bound for general non-overlapping Schwarz preconditioners for dG h⁢p-FEM that is in the order of p2. A polylogarithmic bound for a BDDC preconditioner for a hybridizable discontinuous Galerkin discretization has been developed in [10].

In this paper, we consider the Poisson problem on planar domains. We consider an IETI-DP solver with a scaled Dirichlet preconditioner. The proof follows the abstract framework from [26] and is a continuation of the paper [32], where we have analyzed IETI-DP methods for conforming discretizations. The presented theory covers three different choices for the primal degrees of freedom: vertex values only (Algorithm A), edge averages only (Algorithm B), and the combination of both (Algorithm C). We prove that the condition of the preconditioned IETI-DP system is bounded by


for Algorithms A and C and by


for Algorithm B, where 𝑝 is the spline degree, hk are the grid sizes and Hk are the patch sizes, NΓ⁢(k) contains the indices of the patches that share an edge with the 𝑘-th patch, and δ≥δ* is a suitably chosen penalty parameter. Both the constant 𝐶 and the optimal penalty parameter δ*>0 are independent of the grid sizes, the patch sizes, the spline degree, and the smoothness of the splines. Hence, the theory covers all discretizations where the smoothness within the patches is between C0 and Cp-1. Up to the knowledge of the authors, so far, there is no condition number estimate that is explicit in the spline degree 𝑝. Moreover, we are not aware of an analysis that covers Algorithm B. Such approaches might be of interest if the overall computational domain is decomposed into patches in a way that allows T-junctions or in the case of moving patches like in the case of rotating electrical machines. Note that the constant 𝐶 depends on the chosen parameterization of the geometry. Numerical experiments indicate that the convergence of the IETI-DP solver only suffers mildly from distorted geometries, see, e.g., [34]. A theoretical analysis of the dependence on the geometry is not known to the authors.

The structure of this paper is as follows. In Section 2, we introduce our model problem and its discretization using SIPG. Then, in Section 3, we present the IETI-DP solver. Section 4 is devoted to the proof of the condition number bounds. The numerical results are shown in Section 5. In Section 6, we summarize our findings and give some further remarks.

2 The Model Problem and Its Discretization

We consider the discretization of a homogeneous Poisson problem using multi-patch Isogeometric Analysis, where the coupling of the numerical solution between the individual patches is realized using an SIPG approach. In Section 2.1, we introduce the model problem. The representation of the computational domain is discussed in Section 2.2. The function spaces are introduced in Section 2.3; the variational formulation of the discretized problem is then given in Section 2.4.

Since this paper extends the results of the previous paper [32], which covered conforming discretizations, we aim to use the same notation as in the aforementioned paper. To keep the paper self-contained, we briefly reintroduce the notation.

2.1 The Model Problem

We consider an open, bounded, and simply connected domain Ω⊂R2 with Lipschitz boundary ∂⁡Ω. We are interested in solving the homogeneous Poisson problem: given f∈L2⁢(Ω), find u∈H01⁢(Ω) such that

(2.1)∫Ω∇⁡u⋅∇⁡v⁢d⁢x=∫Ωf⁢v⁢dx for all⁢v∈H01⁢(Ω),

Throughout the paper, we denote by L2⁢(Ω), L∞⁢(Ω), and Hs⁢(Ω), s∈N, the usual Lebesgue and Sobolev spaces, respectively. Moreover, H01⁢(Ω)⊂H1⁢(Ω) is the subspace of functions whose trace vanishes on ∂⁡Ω. We use the standard scalar products (⋅,⋅)L2⁢(Ω) and (⋅,⋅)H1⁢(Ω):=(∇⋅,∇⋅)L2⁢(Ω), norms ∥⋅∥L2⁢(Ω), ∥⋅∥L∞⁢(Ω), ∥⋅∥H1⁢(Ω), and seminorm |⋅|H1⁢(Ω).

2.2 Representation of the Computational Domain

The computational domain Ω is given as the composition of 𝐾 non-overlapping subdomains Ω(1),…,Ω(K), i.e., each Ω(k) is open and bounded with Lipschitz boundary,

(2.2)Ω¯=⋃k=1KΩ(k)¯ and Ω(k)∩Ω(ℓ)=∅ for all⁢k≠ℓ,

where T¯ denotes the closure of the set 𝑇. We refer to the subdomains Ω(k) as patches. We assume that every patch Ω(k) is the image of a geometry function


that can be continuously extended to the closure of Ω^.

Although in IgA commonly B-spline or NURBS functions are used for the representation of the geometry, this is not required for the proposed method. We only assume that the geometry function is sufficiently smooth such that the following assumption holds.

Assumption 1

There is a constant C1>0 such that

∥∇⁡Gk∥L∞⁢(T)≤C1⁢Hk and ∥(∇⁡Gk)-1∥L∞⁢(T)≤C1⁢1Hk

for T∈{Ω^,∂⁡Ω^} and all k=1,…,K, where Hk>0 is the diameter of the patch Ω(k).

This assumption imposes quality criteria on the representation of the geometry. It might be possible to reduce the constant C1 by a reparameterization of the geometry, i.e., by finding another function G~k such that G~k⁢(Ω^)=Gk⁢(Ω^). Moreover, also subdividing a patch into one or more patches might help in reducing the constant C1.

Note that equation (2.2) guarantees that every interface between two patches is exactly parameterized by the geometry functions of both involved patches. This, however, does not imply that the parameterizations have to agree. For the setup of the solver, the interfaces between two patches have to consist of whole edges, which means that there are no T-junctions. The next assumption guarantees this.

Assumption 2

For any two patch indices k≠ℓ, the intersection Ω(k)¯∩Ω(ℓ)¯ is either a common edge (with the neighboring vertices), a common vertex, or empty.

This assumption also ensures that the pre-images Γ^D(k):=Gk-1⁢(∂⁡Ω∩∂⁡Ω(k)) of the (Dirichlet) boundary ΓD=∂⁡Ω consist of whole edges.

For any two neighboring patches Ω(k) and Ω(ℓ), the common edge is denoted by Γ(k,ℓ)=Γ(ℓ,k), and its pre-images by Γ^(k,ℓ):=Gk-1⁢(Γ(k,ℓ)) and Γ^(ℓ,k):=Gℓ-1⁢(Γ(ℓ,k)). We denote the two end-points of Γ(ℓ,k) by x(k,ℓ,1) and x(k,ℓ,2). We collect the indices of patches sharing an edge in a set

NΓ⁢(k):={ℓ≠k:Ω(k)⁢and⁢Ω(ℓ)⁢share one edge}.

The indices of the patches that contain a certain corner are collected in the set


Moreover, we require that the number of neighbors of a patch is uniformly bounded.

Assumption 3

There is a constant C2 such that |P⁢(x)|≤C2 for all corners 𝐱.

2.3 Local Function Spaces

After the introduction of the computational domain, we establish the isogeometric function spaces. In IgA, those function spaces are either B-spline or NURBS functions. In this paper, we focus on B-splines. Let p∈N:={1,2,3,…} be a given spline degree. For ease of notation, we assume that it is the same for all patches. For n∈N, a 𝑝-open knot vector Ξ=(ξ1,…,ξn+p+1) is the building block for the B-spline basis (B⁢[p,Ξ,i])i=1n defined via the Cox–de Boor formula, cf. [9, equations (2.1) and (2.2)]. This basis spans the univariate spline space


Let Ξ(k,1) and Ξ(k,2) be two 𝑝-open knot vectors over (0,1). To get a multivariate spline space V^(k) over Ω^, we tensorize the two univariate spline spaces. The transformation of V^(k) to the physical domain is defined by the pull-back principle. We denote the resulting space by V(k). So we define

V^(k):={v∈S⁢[p,Ξ(k,1)]⊗S⁢[p,Ξ(k,2)]:v|Γ^D(k)=0} and V(k):=V^(k)∘Gk-1,

where v|T denotes the restriction of 𝑣 to 𝑇 (trace operator). The basis for the space V^(k) consists of only those tensor-product basis functions over Ω(k) that vanish on the pre-image of the Dirichlet boundary Γ^D(k). Say, the total number of basis functions of V^(k) is N(k)=NI(k)+NΓ(k). The number NI(k) denotes the number of basis function that are supported only in the interior of the patch, whereas NΓ(k) accounts for the number of basis functions that contribute to the boundary of the patch. These definitions give rise to an ordered basis (2.3)

wherej1∈{1,…,n(k,1)}andj2∈{1,…,n(k,2)}such thatϕ^|Γ^D(k)=0},
ϕ^i(k)|∂⁡Ω^=0⇔i∈{1,…,NI(k)} and ϕ^i(k)|∂⁡Ω^≠0⇔i∈NI(k)+{1,…,NΓ(k)}.
The pull-back principle gives a basis for V(k),

Φ(k):=(ϕi(k))i=1N(k) and ϕi(k):=ϕ^i(k)∘Gk-1.

We assume in the following that the grids on each of the patches are quasi-uniform.

Assumption 4

There are grid sizes h^k>0 for k=1,…,K and a constant C3>0 such that


holds for all i=1,…,n(k,δ)+p and all δ=1,2 red with ξi+1(k,δ)≠ξi(k,δ).

The corresponding grid size on the physical domain is defined via hk:=h^k⁢Hk. For any two patches Ω(k) and Ω(ℓ) sharing an edge, we define

h^k⁢ℓ:=min⁡{h^k,h^ℓ} and hk⁢ℓ:=min⁡{hk,hℓ}.

2.4 The Isogeometric SIPG Problem

Since we do not impose any strong coupling of the solution at this point, the global approximation space is the product space of the local spaces V(k), i.e.,


Using these function spaces, the dG discretization of the model problem (2.1) is given by the following: find u=(u(1),…,u(K))∈V such that

(2.4)ah⁢(u,v)=∫Ωf⁢v⁢dx for all⁢v∈V,


and δ>0 is a penalty parameter and nk is the outer normal vector of ∂⁡Ω(k).

Due to [38, Theorem 8], the parameter 𝛿 can always be chosen independently of the spline degree 𝑝 and mesh sizes hk such that the bilinear form ah⁢(⋅,⋅) is bounded and coercive in the dG-norm


Note that 𝛿 depends on the constant C1 from Assumption 1. Similar results have been shown previously, but up to the knowledge of the authors, the dependence on 𝑝 has only been addressed in [38]. Since we have coercivity and boundedness, the Lax–Milgram lemma ensures the existence of a solution of (2.4) and its uniqueness. If the solution 𝑢 to the continuous problem (2.1) satisfies u|Ω(k)∈Hq+1⁢(Ω(k)) for all k=1,…,K with 1≤q≤p, then the discretization error satisfies


for some constant C>0 independent of ℎ, see [38, Theorems 11 and 12], where more error estimates, including results on their dependence on the spline degree, are available.

3 The dG IETI-DP Solver

In this section, we propose an IETI-DP solver for (2.4). As for any tearing and interconnecting method, we have to choose local spaces, cf., e.g., [14]. As it has been done previously for IETI methods, the local spaces are chosen on a per-patch basis, cf. [22] and later publications. In the case of dG, the choice is not completely straight-forward. We follow the approach that has been proposed in [12] and that has been used in an IgA context in [17, 16] and others: for each patch, we introduce a local subspace that consists of the degrees of freedom associated to that patch and of the degrees of freedom associated to artificial interfaces. The basis functions on the artificial interface coincide with the traces of basis functions of the neighboring patch, when restricted to the interface. We introduce these local function spaces and the corresponding bilinear form in Section 3.1. In Section 3.2, we restrict the problem to the skeleton. As a next step, we introduce jump matrices that are used for the weak coupling of the degrees of freedom on an interface and on the artificial interface of the neighboring subdomain, see Section 3.3. Following the concept of dual-primal methods, in Section 3.4, we introduce the primal degrees of freedom, where we consider vertex values (Algorithm A), edge averages (Algorithm B), and a combination of both (Algorithm C) as options. Having the required ingredients, we pose the IETI-DP system and the proposed preconditioner in Section 3.5. In Section 3.6, we summarize all steps that are required for a computational realization of the solver.

3.1 Local Function Spaces and Bilinear Forms

The patch-local subspace consists of the local functions on the particular patch and of those functions from the neighboring patches required to realize the bilinear forms m(k)⁢(⋅,⋅) and r(k)⁢(⋅,⋅). We define the enriched function space


where V(k,ℓ):={v(ℓ)|Γ(k,ℓ):v(ℓ)∈V(ℓ)} is the trace of V(ℓ). We write


where ve(k)∈Ve(k), v(k)∈V(k) and v(k,ℓ)∈V(k,ℓ).

Note that


forms a basis of V(k,ℓ), which is uniquely defined up to the ordering of the basis functions. We denote the basis by Φ(k,ℓ):=(ϕi(k,ℓ))i=1N(k,ℓ).

A visualization is given in Figure 1, where all basis functions are represented using a symbol located at their maximum. The basis functions for each of the patches are represented via different symbols. The symbols that are located on the edges and corners represent the only functions that are supported on that edge or corner, respectively. Besides the patches themselves, we have artificial interfaces. While they are co-located with the (standard) interfaces, they are treated as separate entities on which the corresponding trace spaces V(k,ℓ) live. The basis functions living on these trace spaces are marked with the same kind of symbol as the basis function from the original space. Examples of basis functions in the interior, on the interfaces, and on the artificial interfaces are visualized in Figure 2.

Figure 1 Local spaces with artificial interfaces

Figure 1

Local spaces with artificial interfaces

Figure 2 Basis functions in the interior (yellow), on an interface (red), and on an artificial interface (blue)

Figure 2

Basis functions in the interior (yellow), on an interface (red), and on an artificial interface (blue)

Following the structure of (3.1), we define a basis

Φe(k):=(ϕe,i(k))i=1Ne(k) with⁢Ne(k):=dim⁡Ve(k)

for the space Ve(k) as follows. We start with the basis functions from Φ(k),

ϕe,i(k):=(ϕi(k),(0,…,0)) for⁢i=1,…,N(k).

Then the basis functions follow from the neighboring patches. Let NΓ⁢(k)={ℓ1,ℓ2,…,ℓL}. We define(3.2)

From this construction and from (2.3), we know that the first NI(k) basis functions of the basis Φe(k) live in the interior of the patch Ω(k). The following NΓ(k) basis functions live on the interfaces. The remaining basis functions live on the artificial interfaces.

On the spaces Ve(k), we define the bilinear forms ae(k)⁢(⋅,⋅) and de(k)⁢(⋅,⋅) and the linear functional ⟨fe(k),⋅⟩ analogous to the global objects ah⁢(⋅,⋅) and d⁢(⋅,⋅) by

where we write with a slight abuse of notation
The bilinear form de(k)⁢(⋅,⋅) induces the local dG-norm ∥⋅∥de(k)2:=de(k)⁢(⋅,⋅), the bilinear form ae(k)⁢(⋅,⋅) the local energy norm ∥⋅∥ae(k)2:=ae(k)⁢(⋅,⋅).

By discretizing, we obtain the local linear system



A(k)=[ae(k)⁢(ϕe,j(k),ϕe,i(k))]i,j=1Ne(k) and f¯e(k)=[⟨fe(k),ϕe,i(k)⟩]i=1Ne(k).

The vector u¯e(k)=[ue,i(k)]i=1Ne(k) is the coefficient vector representing the function ue(k)=∑i=1Ne(k)ue,i(k)⁢ϕe,i(k).

We subdivide the stiffness matrix and the load vector into blocks

A(k)=(AII(k)AI⁢Γ(k)AΓ⁢I(k)AΓ⁢Γ(k)) and f¯e(k)=(f¯I(k)f¯Γ(k)),

where the first row and the first column correspond to the first NI(k) basis functions, i.e., to those supported in the interior of the patch Ω(k). So the remainder accounts for both the standard interfaces and the artificial interfaces.

3.2 The Problem on the Skeleton

In this section, we eliminate the interior degrees of freedom by building a Schur complement system, which is equivalent to (3.3),



(3.5)S(k):=AΓ⁢Γ(k)-AΓ⁢I(k)⁢(AII(k))-1⁢AI⁢Γ(k) and g¯e(k):=f¯Γ(k)-AΓ⁢I(k)⁢(AII(k))-1⁢f¯I(k).

Once w¯e(k) has been computed, we get the solution u¯e(k) of system (3.3) by


The collection of the systems (3.4) for all patches yields a block diagonal linear system


Next, we introduce the building blocks in order to rewrite the Schur complement system in a variational setting. First, we define function spaces on the skeleton,

W(k):={v|∂⁡Ω(k):v∈V(k)} and W(k,ℓ):=V(k,ℓ). Analogously to (3.1), an extended function we(k)∈We(k) has the form

where w(k)∈W(k) and w(k,ℓ)∈W(k,ℓ). A basis for We(k) is given by

Φ˘e(k):=(ϕe,NI(k)+i(k))i=1N˘e(k) with⁢N˘e(k):=Ne(k)-NI(k).

As above, the coefficient vector w¯e(k)=[we,i(k)]i=1N˘e(k) represents the function


The Schur complements realize the discrete harmonic extension HA(k):We(k)→Ve(k) with respect to the energy norm ∥⋅∥ae(k), which is defined as follows. For any we(k)=(w(k),(w(k,ℓ))ℓ∈NΓ⁢(k)), we have


where u(k) is such that

ae(k)⁢((u(k),(w(k,ℓ))ℓ∈NΓ⁢(k)),(v(k),0))=0for all⁢v(k)∈V0(k),
where V0(k):=V(k)∩H01⁢(Ω(k)). Note that the standard discrete harmonic extension Hh(k):W(k)→V(k) is defined by Hh(k)⁢w(k)=u(k), where u(k) is such that
a(k)⁢(u(k),v(k))=0for all⁢v(k)∈V0(k).
In variational form, problem (3.7) reads as follows: find w=(we(1),…,we(K))∈W such that

∑k=1Kae(k)⁢(HA(k)⁢we(k),HA(k)⁢qe(k))⏟=⁣:s⁢(w,q)=∑k=1K⟨f,HA(k)⁢qe(k)⟩⏟=⁣:⟨g,q⟩ for all⁢q=(qe(1),…,qe(K))∈W.

3.3 The Jump Matrices

The next step is the introduction of constraints that yield the continuity of the solution. For any two patches Ω(k) and Ω(ℓ) sharing an edge Γ(k,ℓ) and any basis function ϕi(k) supported on that edge, we introduce a constraint

(3.9)we,i*(k)-we,j*(ℓ)=0 with⁢i*=i-NI(k)⁢and⁢j*=Υi(ℓ,k)-NI(ℓ),

where Υi(ℓ,k) is as in (3.2) and we,i*(k) and we,j*(k) are the coefficients of the functions we(k) and we(ℓ). The index shifts by NI(k) and NI(ℓ) account for the restriction to We(k) and We(ℓ), cf. (3.8). This constraint ensures that the function values of the solution on the artificial interfaces coincide with the function values of the solution on the corresponding standard interfaces.

If the corner values are chosen as primal degrees of freedom (Algorithms A and C), we omit the constraints for the basis functions ϕi(k) that are supported on a corner of Ω(k). A visualization of the constraints is given in Figure 3.

If only the edge averages are chosen as primal degrees of freedom (Algorithm B), we realize the constraints at the corners in a fully redundant way. This means that, besides the constraints (3.9), for every basis function ϕi(k) that is supported on a corner of Ω(k), we additionally introduce constraints of the form

(3.10)we,j1*(ℓ1)-we,j2*(ℓ2)=0 with⁢j1*=Υi(ℓ1,k)-NI(ℓ1)⁢and⁢j2*=Υi(ℓ2,k)-NI(ℓ2),

where ℓ1<ℓ2 such that the said corner lies between the edges Γ(k,ℓ1) and Γ(k,ℓ2). These additional constraints also ensure the continuity between the different artificial interfaces. A visualization of this case is given in Figure 4.

Figure 3 Omitting vertices (Algorithms A and C)

Figure 3

Omitting vertices (Algorithms A and C)

Figure 4 Fully redundant (Algorithm B)

Figure 4

Fully redundant (Algorithm B)

We collect the constraints of the form (3.9) and (3.10) in a jump matrix 𝐵 such that the constraints are equivalent to B⁢w¯=0. This is done such that each row of the matrix corresponds to one constraint, i.e., each row has only two non-zero entries, one with value 1 and one with value −1. The matrix 𝐵 can also be written as B=(B(1)⁢⋯⁢B(K)).

We are now able to write down a saddle point problem in matrix-vector form that is equivalent to problem (2.4): find (w¯,λ¯) such that


Note that the matrices A(k) and S(k) are non-singular only if the corresponding patch contributes to the Dirichlet boundary ∂⁡Ω. Since in general at least one of the matrices S(k) is singular, the matrix 𝑆 is usually singular. The remedy for this problem are primal degrees of freedom which form a global but small linear system.

3.4 The Primal Degrees of Freedom

As mentioned in the introduction, we consider three options for the choice of primal degrees of freedom: Algorithms A, B, and C. Each of them corresponds to a choice of spaces W~ and W~Δ.

Algorithm A

Algorithm A (Vertex Values)

The space W~ is the subspace of functions where the vertex values agree, i.e.,

W~:={w∈W:w(k)⁢(x(k,ℓ,i))=w(ℓ,k)⁢(x(k,ℓ,i))⁢for all⁢k,all⁢ℓ∈NΓ⁢(k),and all⁢i},

where x(k,ℓ,1) and x(k,ℓ,2) are the two end-points of Γ(k,ℓ). Furthermore, W~Δ satisfies these conditions homogeneously, i.e., W~Δ:=∏k=1KW~Δ(k) and

W~Δ(k):={we(k)∈We(k):w(k)⁢(x(k,ℓ,i))=w(ℓ,k)⁢(x(k,ℓ,i))=0⁢for all⁢ℓ∈NΓ⁢(k)⁢and all⁢i}.
Algorithm B

Algorithm B (Edge Averages)

The space W~ is the subspace of functions where the averages of the function values over the edges agree, i.e.,

W~:={w∈W:∫Γ(k,ℓ)w(k)⁢dx=∫Γ(k,ℓ)w(ℓ,k)⁢dx⁢for all⁢k⁢and all⁢ℓ∈NΓ⁢(k)}.

Furthermore, W~Δ satisfies these conditions homogeneously, i.e., W~Δ:=∏k=1KW~Δ(k) and

W~Δ(k):={we(k)∈We(k):∫Γ(k,ℓ)w(k)⁢dx=∫Γ(k,ℓ)w(ℓ,k)⁢dx=0⁢for all⁢ℓ∈NΓ⁢(k)}.
Algorithm C

Algorithm C (Vertex Values and Edge Everages)

The spaces W~ and W~Δ(k) and W~Δ are the intersections of the corresponding spaces obtained by Algorithms A and B.

Figures 5 and 6 visualize the constraints for Algorithms A and B, respectively. The arrows in Figure 5 indicate which corner degrees of freedom are coupled in Algorithm A. The different lines in Figure 6 indicate the coupling of the edge averages. Interfaces drawn in the same style have to have the same average.

Figure 5 Constraints for Algorithm A

Figure 5

Constraints for Algorithm A

Figure 6 Constraints for Algorithm B

Figure 6

Constraints for Algorithm B

The primal space W~Π is the subspace of energy minimizing functions, cf. [30]. This means that W~Π is 𝑆-orthogonal to W~Δ, i.e.,

(3.11)W~Π:={w∈W~:s⁢(w,q)=0⁢for all⁢q∈W~Δ}.

Let ψ(1),…,ψ(NΠ) be a basis of W~Π. In practice, one usually chooses a nodal basis, where the vertex values and/or the edge averages form the nodal values. The matrix Ψ represents the basis in terms of the basis for the space 𝑊, i.e.,


A representation for the spaces W~Δ(k) can be obtained by full-rank matrices C(k), C(k)⁢w¯e(k)=0⇔we(k)∈W~Δ(k). The block-diagonal collection of the matrices C(k) forms a matrix 𝐶.

3.5 The IETI-DP System and the Scaled Dirichlet Preconditioner

As in [26], the following problem is an equivalent reformulation of (2.4): find (w¯Δ,μ¯,w¯Π,λ¯) such that


The solution for the original problem is obtained by w¯=w¯Δ+Ψ⁢w¯Π. Following the DP approach, this system can be reduced to



F:=(B⁢ 0⁢B⁢Ψ)⁢(SC⊤0C0000Ψ⊤⁢S⁢Ψ)-1⁢(I0Ψ⊤)⏟=⁣:F0⁢B⊤ and d¯:=F0⁢g¯.

System (3.13) is solved with a preconditioned conjugate gradient (PCG) solver. In order to obtain an effective solver, we need a proper preconditioner. We use the scaled Dirichlet preconditioner MsD, defined by


where D∈RNΓ×NΓ is a diagonal matrix defined based on the principle of multiplicity scaling: each coefficient di,i of 𝐷 is assigned the number of constraints for the corresponding degree of freedom plus one.

3.6 The IETI-DP Algorithm

The dG IETI-DP method requires the execution of the following steps.

  • Compute the right-hand sides g¯e(k) according to (3.5).

  • Compute a basis Ψ for the primal space in accordance with (3.11) and (3.12).

  • Compute the primal Schur complement SΠ:=Ψ⊤⁢S⁢Ψ.

  • Solve (3.13) for λ¯ using a PCG solver. This requires the computation of the residual and the application of the preconditioner.

    The computation of the residual w¯^:=F⁢λ¯-d¯ requires the following steps:

    1. Compute q¯^=((q¯^e(1))⊤⁢⋯⁢(q¯^e(K))⊤)⊤:=B⊤⁢λ¯-g¯.

    2. Solve the linear system


      for all k=1,…,K. For Algorithms A and C, one usually reduces the local systems by eliminating the degrees of freedom corresponding to the vertex values and the corresponding Lagrange multipliers.

    3. Solve the (usually small) global linear system SΠ⁢w¯^Π=Ψ⊤⁢q¯^.

    4. Compute the residual w¯^:=B⁢((w¯^Δ(1))⊤⁢⋯⁢(w¯^Δ(K))⊤)⊤+B⁢Ψ⁢w¯^Π.

    Then apply the preconditioner, i.e., compute MsD⁢w¯^.

  • The vector w¯, representing the solution on the skeleton, is computed analogously to steps (2), (3), and (4) based on q=((q(1))⊤)⋯(q(K))⊤)⊤:=B⊤λ¯. Finally, the solution vector u¯ is derived by (3.6).

4 The Condition Number Estimate

In this section, we prove the following condition number estimates.

Theorem 1

Provided that the IETI-DP solver is set up as outlined in the previous sections, the condition number of the preconditioned system satisfies


for Algorithms A and C and


for Algorithm B, where in both cases the constant 𝐶 only depends on the constants from Assumptions 1, 3, and 4.

If some Lagrange multipliers are redundant to each other or to primal constraints, the matrices MsD and 𝐹 are singular. Standard iteration schemes live in the corresponding factor space. Here and in what follows, the condition number is estimated for the restriction to the factor space, cf. [26, Remark 23].

In Section 4.1, we introduce the notation and show some results concerning the well-posedness of the local systems. Then, in Section 4.2, we estimate the norm of BD⊤⁢B⁢u¯, where BD:=B⁢D-1. We give inequalities that allow to estimate terms living on the artificial interface in Section 4.3. Finally, in Section 4.4, we prove Theorem 1.

4.1 Preliminaries

First, we introduce some required notation. We write a≲b if there is a constant c>0 that only depends on the constants from Assumptions 1, 3, and 4 such that a≤c⁢b. Moreover, we write a≂b if a≲b≲a.

The space H1/2 is the standard Sobolev–Slobodeckij space with a norm defined as follows. For T⊂R2 open, bounded, and Lipschitz and Γ⊆∂⁡T, we define the seminorm |⋅|H1/2⁢(Γ)2 by


and the corresponding norm by ∥⋅∥H1/2⁢(Γ)2:=∥⋅∥L2⁢(Γ)2+|⋅|H1/2⁢(Γ)2. As in [32], we define the seminorm


Before we start with the condition number estimate, we show that it is sufficient to get the necessary estimates on the parameter domain.

Lemma 1

For all patches k=1,…,K, the equivalences

  • |u|Hs⁢(Ω(k))≂(Hk)1-s⁢|u∘Gk|Hs⁢(Ω^) for all u∈Hs⁢(Ω(k)) and s∈{0,1}, and

  • |w|Hs⁢(∂⁡Ω(k))≂(Hk)1/2-s⁢|w∘Gk|Hs⁢(∂⁡Ω^) for all w∈Hs⁢(∂⁡Ω(k)) and s∈{0,1/2,1}

hold, where we use H0⁢(⋅):=L2⁢(⋅). The equivalences for the boundary do also hold if ∂⁡Ω(k) and ∂⁡Ω^ and are replaced by Γ(k,ℓ) and Γ^(k,ℓ), respectively.


The results for s∈{0,1} are obtained by Assumption 1, the chain rule and the substitution rule for differentiation, see e.g., [7, Lemma 3.5]. The results for s=1/2 follow by Hilbert space interpolation theory, cf. [1, Theorem 7.23, Theorem 7.31]. ∎

The next lemma states that the bilinear forms ae(k)⁢(⋅,⋅) and de(k)⁢(⋅,⋅) are patchwise equivalent.

Lemma 2

For every patch Ω(k), the relation ae(k)⁢(ue(k),ue(k))≲de(k)⁢(ue(k),ue(k))≲ae(k)⁢(ue(k),ue(k)) holds for all ue(k)∈Ve(k).


Let ue(k)=(u(k),(u(k,ℓ))ℓ∈NΓ⁢(k)) and ve(k)=(v(k)⁢(v(k,ℓ))ℓ∈NΓ⁢(k))∈Ve(k). By Lemma 1, the fact that nk has length 1, the discrete trace inequality [13, Lemma 4.3], Assumption 4, and again Lemma 1, we obtain


Using the Cauchy–Schwarz inequality, estimate (4.1), and NΓ⁢(k)≲1, we get


for some constant c>0, cf. [38] for details. With the choice δ≥16⁢c2, we obtain


i.e., coercivity. Since we obtain boundedness using the Cauchy–Schwarz inequality and (4.2), this finishes the proof. ∎

4.2 Estimates for the Action of BD⊤⁢B

Analogous to [32, Lemma 4.16], the following lemma estimates the effect of BD⊤⁢B.

Lemma 3



with coefficient vector u¯, and let


with coefficient vector w¯ be such that w¯=BD⊤⁢B⁢u¯. Then we have, for each patch Ω(k) and each edge Γ(k,ℓ) connecting the vertices x(k,ℓ,1) and x(k,ℓ,2),


(4.3)Δ(k,ℓ,i):={0for Algorithms A and C,∑j∈P⁢(x(k,ℓ,i))∩NΓ⁢(k)|u(k)⁢(x(k,ℓ,i))-u(j,k)⁢(x(k,ℓ,i))|2for Algorithm B.


We start by recalling how the scaling matrix 𝐷 looks like. Remember 𝐷 is a diagonal matrix. We denote the coefficients by di,i for i=1,…,NΓ. The entries di,i are defined to be one plus the number of Lagrange multipliers of the corresponding degree of freedom 𝑖. Note that we have di,i=2 unless the corresponding degree of freedom is a corner degree of freedom. In that case, di,i takes a value, which is bounded from below by 1 and bounded from above due to Assumption 3. Thus, di,i≂1 in any case.

We will start with Algorithms A and C: simple calculations reveal for Algorithms A and C that

where θ(k,ℓ,i) is the basis function in Φ(k) such that θ(k,ℓ,i)⁢(x(k,ℓ,i))=1, and θ(ℓ,k,i) is defined analogously. Since 𝑢 satisfies the primal constraints, we have u(k)⁢(x(k,ℓ,i))=u(ℓ,k)⁢(x(k,ℓ,i)) and u(k,ℓ)⁢(x(k,ℓ,i))=u(ℓ)⁢(x(k,ℓ,i)) and thus
Therefore, we obtain, by using the triangle inequality,


for the L2-norm. For the H1/2-seminorm, we get using the triangle inequality that


Analogously, we obtain, using the triangle inequality,


which finishes the proof for Algorithms A and C.

For Algorithm B, we have that w(k) on Γ(k,ℓ) can be expanded as

where dni,ni and dmi,mi denote the entries of the matrix 𝐷 for the corresponding dofs and we use u(ℓ,ℓ):=u(ℓ). Note that θ(k,ℓ,i) behaves like max⁢{0,1-|x-x(k,ℓ,i)|/h(k)}p. Hence,

∥θ(k,ℓ,i)∥L2⁢(Γ(k,ℓ))2≂hkp and |θ(k,ℓ,i)|H1/2⁢(Γ(k,ℓ))2≂1.

Certainly, we also have ∥θ(k,ℓ,i)∥L∞⁢(Γ(k,ℓ))2=1. An application of the triangle inequality for the corner expressions yields the stated result for Algorithm B. ∎

4.3 Estimates for Contributions from Artificial Interfaces

The following lemma allows to estimate the H1/2-seminorm on the artificial edges.

Lemma 4

For any two patches Ω(k) and Ω(ℓ) that share an edge Γ(k,ℓ),


holds for all ue(k)=(u(k),(u(k,ℓ))ℓ∈NΓ⁢(k))∈Ve(k).


Let c∈R be arbitrary but fixed. Using the triangle inequality, we have


Using [1, Theorem 5.2, equation (3)] and the triangle inequality, we obtain further


Using the equivalence of the norms on the parameter domain and the physical domain, cf. Lemma 1, and an inverse inequality, cf. [32, Lemma 4.3], and using hk⁢ℓ=min⁡{hk,hℓ}, we obtain further


where C>0 is an appropriately chosen constant (that only depends on the constant from Assumption 1). Using a⁢(b+c)≤a2+(b+c)2/4≤a2+b2/2+c2/2, we have


By subtracting 12⁢A, we immediately obtain


for all ue(k)∈Ve(k). Since |u(k,ℓ)|H1/2⁢(Γ(k,ℓ))2≤A and since the bound for 𝒜 holds for all c∈R, the Poincaré inequality yields the desired result. ∎

The next step is to show that a similar estimate holds for the seminorm |⋅|L∞0⁢(T). Before we can prove that result, we need the following auxiliary result.

Lemma 5

∥u∥L∞⁢(0,1)2≤2⁢∥u∥L2⁢(0,1)⁢∥u∥H1⁢(0,1) holds for all u∈H1⁢(0,1).


Since 𝑢 is continuous, 𝑢 takes its maximum for some z∈[0,1]. Then, by the fundamental theorem of differential and integral calculus, we can write |u⁢(z)|2=|u⁢(t)|2-∫ztu⁢(s)⁢u′⁢(s)⁢ds. The desired result is obtained by integrating over the unit integral and by utilizing the Cauchy–Schwarz inequality. ∎

Lemma 6

For any two patches Ω(k) and Ω(ℓ) that share an edge Γ(k,ℓ),


holds for all ue(k)=(u(k),(u(k,ℓ))ℓ∈NΓ⁢(k))∈Ve(k).


Using the triangle inequality, we obtain


As a next step, we apply Lemma 5 to the difference |u(k)-u(k,ℓ)|L∞0⁢(Γ(k,ℓ))2. Since the L∞-norm does not change if we are on the physical or the parameter domain, we can apply Lemma 5 and subsequently utilize the equivalence of the norms on the parameter and the physical domain, cf. Lemma 1, to get


The triangle inequality allows to estimate


for all c∈R. The equivalence of the norms on the physical domain and the parameter domain, cf. Lemma 1, an inverse estimate, cf. [32, Lemma 4.3], and a⁢(b+c)≲a2+b2+c2 give


Since this holds for all c∈R, the Poincaré inequality gives further


Lemma 4 yields the final result. ∎

Lemma 4.17 in [32] states that the term Δ(k,ℓ,i) can be estimated by expressions which only involve patches sharing an edge. Now we prove a variant of [32, Lemma 4.18] that fits our needs.

Lemma 7

For any two patches Ω(k) and Ω(ℓ), sharing an edge Γ(k,ℓ) which connects the vertices x(k,ℓ,1) and x(k,ℓ,2) and ∫Γ(k,ℓ)u(k)⁢(s)-u(ℓ,k)⁢(s)⁢d⁢s=0 and i=1,2, we have


for all u=(ue(1),…,ue(K))∈W~, where Λ=1+log⁢p+maxj=1,…,K⁢log⁢Hjhj and Δ(k,ℓ,i) is as in (4.3).


Using the triangle inequality, we have




and observe that the Cauchy–Schwarz inequality and Assumption 1 yield


The triangle inequality, and then (4.5) and [32, Lemma 4.18] yield


Now, we estimate |u(ℓ)⁢(x(k,ℓ,i))-u(ℓ,k)⁢(x(k,ℓ,i))|2. Thus, we set u^(ℓ)=u(ℓ)∘Gℓ and u^(ℓ,k)=u(ℓ,k)∘Gk. With unitary transformations, i.e., rotations and reflections, we can transform the patches such that the pre-image of Γ(k,ℓ) is (0,1)×{0} and that the pre-image of x(k,ℓ,i) is 0. Let u~(ℓ):=u^(ℓ)⁢(⋅,0) and u~(ℓ,k):=u^(ℓ,k)⁢(⋅,0). Let 𝜂 be the largest value such that u~(ℓ) and u~(ℓ,k) are polynomial on (0,η). Using the quasi-uniformity assumption, cf. Assumption 4, we obtain η≂h^k⁢ℓ. Arguments that are analogous to those used in the proof of Lemma 5, together with the inverse inequality [37, Theorem 4.76, equation (4.6.5)] yield


Since the L2-norm on (0,1) is larger than the L2-norm on (0,η) and η≂h^k⁢ℓ,


An application of Lemma 1 finally yields


The combination of (4.4), (4.6), and (4.7) finishes the proof. ∎

4.4 Proof of the Main Theorem

Before we give a proof of the main theorem, we estimate the sum of we(k)-seminorms over all patches.

Lemma 8

Let 𝑢 and 𝑤 be as in Lemma 3. Then we have


where Λ:=1+log⁡p+maxk=1,…,K⁡log⁡Hkhk.


Lemma 3 yields


Lemmas 4 and 6 allow to estimate the first sum. This finishes the proof for Algorithms A and C since the second sum vanishes in this case. For Algorithm B, the second sum is estimated using Lemma 7 and |NΓ⁢(k)|≤4. ∎

Lemma 9

Let 𝑢 and 𝑤 be as in Lemma 3. Then we have



(4.8)Σ:={1for Algorithms A and C,δ⁢p⁢maxk=1,…,K⁡maxℓ∈NΓ⁢(k)⁡hkhℓfor Algorithm B.


Lemma 3 and the observation that ℓ∈NΓ⁢(k)⇔k∈NΓ⁢(ℓ) yield


Since Δ(k,ℓ,i)=0 for Algorithms A and C, this finishes the proof for this case. Now, consider Algorithm B. Lemma 7 and the observation that ℓ∈NΓ⁢(k)⇔k∈NΓ⁢(ℓ) yield


Then hk⁢ℓ=min⁡{hk,hℓ} and ℓ∈NΓ⁢(k)⇔k∈NΓ⁢(ℓ) finish the proof. ∎

Now we are able to prove the bound for the condition number of the preconditioned dG IETI method as stated in Theorem 1.

Proof of Theorem 1

Theorem 22 in [26] states that


where u¯ is the coefficient vector associated to the function u=(ue(1),…,ue(K))=((u(1),(u(1,ℓ))ℓ∈NΓ⁢(1)),…). So let 𝑢 be arbitrary but fixed, and let the function w=(we(1),…,we(K))=((w(1),(w(1,ℓ))ℓ∈NΓ⁢(1)),…) with coefficient vector w¯ be such that w¯=BD⊤⁢B⁢u¯. Let ve(k)∈Ve(k) be an arbitrary extension of we(k). Lemma 2 yields


Note that the second sum does not depend on w0(k). Thus, the infimum refers only to the H1-seminorm, which means that the seminorm coincides with the seminorm of the (standard) discrete harmonic extension. Hence, we have


First, we estimate 𝒜. Due to [32, Theorem 4.2] and [32, Lemma 4.15], we get


Using Λ≥1 and Lemma 8, we obtain further


Using [32, Lemma 4.15 and Theorem 4.2] and |NΓ⁢(k)|≤4, we further estimate


We estimate the second sum using [32, Lemma 4.14], the Poincaré inequality, Lemma 1, and |NΓ⁢(k)|≲1 and obtain


To estimate the second sum in (4.10), we use Lemma 9 to get


where Σ is as in (4.8). Then δ≳1, Λ≥1, Σ≥1, (4.10), (4.11), and (4.12) yield


The combination of this estimate and (4.9) finishes the proof. ∎

5 Numerical results

In this section, we present the results of our numerical experiments that illustrate the presented convergence theory. We consider the Poisson problem

where the computational domain Ω is one of the domains depicted in Figure 7.

Figure 7

Computational domains and the decomposition into patches

(a) Circular ring


Circular ring

(b) Yeti-footprint



The first domain (Figure 7 (a)) is a circular ring consisting of 12 patches. Each of them is the image of a NURBS mapping of degree 2. The second domain (Figure 7 (b)) is the Yeti-footprint which is composed of 84 patches. In this case, all of the patches are parameterized using B-spline functions, again of degree 2.

The numerical experiments are carried out with B-spline discretization spaces of maximum smoothness on grids that are constructed as follows. For the circular ring, the coarsest grid on each patch only consists of one element, i.e., the discretization space (on the parameter domain) for each patch consists only of polynomial functions. For the Yeti-footprint, the coarsest grid for the 20 patches in the bottom of the domain consists of two elements per patch. The grid is constructed by adding an edge that connects the midpoints of the two longer sides of the patch. The other patches of the Yeti-footprint only consist of one element. For all three domains, the finer grids are constructed by refinement. For the first refinement step, i.e., for r=1, we insert one knot into each knot span. The new knot is not located in the center, but at 4/9 times the length of the knot span if the patch index 𝑘 is even and at 6/11 times that length if 𝑘 is odd. The subsequent refinement steps r=2,3,… are done uniformly. This refinement procedure yields discretizations that are non-matching at the interfaces. In the numerical experiments underlying Tables 1–8, the patches are not modified when the grid is refined, i.e., the patch diameters Hk remain the same for the respective numerical tests.

For these discretization spaces, we set up a dG IETI discretization as introduced in Section 3. The penalty term is chosen such that hk⁢ℓ is an estimate for the size of the smaller side of the smallest element of the physical patches Ω(k) and Ω(ℓ), and δ=3 for Algorithms A and C and δ=1 for Algorithm B. To solve system (3.13), we use a preconditioned conjugate gradient (PCG) method with the proposed scaled Dirichlet preconditioner MsD. The implementation is done using G+Smo [27]. The local subproblems are solved with the sparse direct solver from the PARDISO project[1].

We start the iteration with a randomly sampled vector with entries in the interval [-1,1]. The stopping criterion is chosen as follows. The iteration stops if the l2-norm of the residual vector drops below 10-6 times the residual of the right-hand side. For each experiment, we show the number of iterations (it) required to reach the stopping criterion and an estimate of the condition number (𝜅) of the preconditioned system matrix MsD⁢F that has been derived using the PCG algorithm.

5.1 Tests on the Circular Ring

We start with the numerical experiments for the circular ring (Figure 7 (a)). Table 1 shows the results for Algorithm A (vertex values). We observe that the condition number grows not more than log2⁡H/h. Moreover, we also observe a weak growth in the spline degree, which might be smaller than the linear growth predicted by the theory. Table 2 shows the results for Algorithm B (edge averages). Here, the condition numbers are much larger than for Algorithm A, which can be explained as follows. The condition number bound in Theorem 1 depends on 𝛿 only for Algorithm B. Estimating hk⁢ℓ conservatively, which is necessary in order to obtain a coercive bilinear form, can be interpreted as increasing 𝛿. While this affects the condition number significantly, the iteration counts are not affected as much since the eigenvalues seem to be slightly clustered. The dependence on the grid size and the spline degree seem to be the same as for Algorithm A. Table 3 comprises the data collected for Algorithm C (vertex values + edge averages). As expected, this approach yields the smallest values for the condition numbers. The condition number seems to grow only like log⁡H/h in the grid size.

Table 1

Iteration counts (it) and condition numbers (𝜅), Algorithm A, triple ring


Table 2

Iteration counts (it) and condition numbers (𝜅), Algorithm B, triple ring


Table 3

Iteration counts (it) and condition numbers (𝜅), Algorithm C, triple ring


Table 4

Iteration counts (it) and condition numbers (𝜅), Algorithm A, Yeti-footprint


Table 5

Iteration counts (it) and condition numbers (𝜅), Algorithm B, Yeti-footprint

Figure 8 Condition numbers, Algorithm A, Yeti-footprint

Figure 8

Condition numbers, Algorithm A, Yeti-footprint

Figure 9 Condition numbers; Algorithm B, Yeti-footprint

Figure 9

Condition numbers; Algorithm B, Yeti-footprint

Table 6

Solving times (sec.) for Algorithms A and C, single core, Yeti-footprint


5.2 Test on the Yeti-Footprint

Next, we take a look at numerical results on the Yeti-footprint (Figure 7 (b)) if the right-hand side function is less smooth than in the model problem introduced at the beginning of the section. For the Yeti-footprint, we consider a patchwise defined right-hand side that is the zero function on patches with odd index and 2⁢π2⁢sin⁡(π⁢x)⁢sin⁡(π⁢y) on patches with even index. Table 4 reports on the results for Algorithm A. We can again see that the increase in the condition number is like log2⁡H/h, and the increase in 𝑝 is sub-linear, almost like p. Table 5 gives the condition number estimates for Algorithm B, where we observe that the condition number grows almost linearly in the spline degree 𝑝 (which indicates that the convergence theory might be sharp in this respect). Again, Algorithm B yields the largest condition numbers. The graphs in Figures 8 and 9 visualize the condition number growths for Algorithms A and B on the Yeti-footprint, respectively. Compared to the circular ring, the condition numbers for the Yeti-footprint are smaller, which might be connected to a more regular geometry mapping. Lastly, we compare the solving times between Algorithms A and C in Table 6. The solving times are given in seconds and are measured on a single core on the Radon1 cluster[2] located in Linz. We see that Algorithm C is considerably faster than Algorithm A despite the fact that the primal problem is larger for Algorithm C (392 dofs) than for Algorithm A (192 dofs).

5.3 Dependence on hk/hℓ

The next experiment shows that the dependence of the condition number on the ratio of the grid sizes of neighboring patches is only present for Algorithm B. For this purpose, we compare Algorithm A with Algorithm B on both computational domains for particular grids that are constructed as follows. Starting from the initial grid introduced above, we applied 5 refinement steps for the ring domain and 4 refinement steps for the Yeti-footprint in the same way mentioned above. Then we uniformly refine the grids on all patches Ω(k), where 𝑘 is even, additionally e∈{1,2,3} times. Thus, the grid sizes between patches with even and odd degrees vary by a factor of 2e. We observe in Table 7 for the circular ring and in Table 8 for the Yeti-footprint that the condition number is almost independent of 𝑒 if Algorithm A is used, while it increases like 2e if Algorithm B is chosen. This means that the condition number grows linearly in the ratio of the grid sizes, which coincides with the prediction of the convergence theory.

Table 7

Iteration counts (it) and condition numbers (𝜅), triple ring


Table 8

Iteration counts (it) and condition numbers (𝜅), Yeti-footprint


5.4 Discretization Errors

In this section, we report the discretization errors and the corresponding rates. Here, we focus on the circular ring depicted in Figure 7 (a). We slightly modify the problem in order to know the exact solution. We consider inhomogeneous Dirichlet boundary conditions u=u* on ∂⁡Ω, where


and observe that u* is the exact solution of the boundary value problem. The refinement procedure and the solvers are as outlined above. We stop the iteration if the ℓ2-norm of the residual vector drops below 10-8 times the residual of the right-hand side vector. In Table 9 and Figure 10, we present the discretization error ϵr,p and the corresponding rate ιr,p, which are defined by

ϵr,p:=|uhr,p-u*|H1⁢(Ω)|u*|H1⁢(Ω) and ιr,p:=log2⁡(|uhr-1,p-u*|H1⁢(Ω)|uhr,p-u*|H1⁢(Ω)).

We give the results for Algorithm A. We observe that the error rates are better than predicted by (2.5). They have the tendency to decrease as the grid is refined. Since the discretization error is almost independent of the method used to solve the linear system, we obtain for Algorithm B and C almost the same results.

Table 9

Discretization errors (ϵr,p) and rates (ιr,p), Algorithm A, triple ring

Figure 10 Discretization errors (ϵr,p)(\epsilon_{r,p}), Algorithm A, triple ring

Figure 10

Discretization errors (ϵr,p), Algorithm A, triple ring

5.5 Patch Scaling Tests

Finally, we examine the dependence of the condition number on the number of patches in practice. We use the circular ring with its K=12 patches as initial configuration. This serves as the computational domain with the fewest number of patches. The decompositions into K=48, 192, and 768 patches are obtained by splitting the patches uniformly. On each patch, we consider spline spaces that are obtained by r=5 refinement steps that are done in the same way as mentioned above. Table 10 shows that the condition numbers for Algorithms A and B increase only mildly. In the first step from 12 patches to 48 patches, one obtains the largest increase since in the initial configuration with 12 patches still most of the patches contribute to the Dirichlet boundary.

Table 10

Condition numbers (𝜅) depending on the number of patches 𝐾, triple ring


6 Conclusions

We have extended the theory from [32], where we established 𝑝-explicit condition number estimates for continuous Galerkin IETI-DP solvers, to symmetric interior penalty discontinuous Galerkin discretizations. These discretizations are of particular interest in the area of Isogeometric Analysis since they allow geometry functions and spline spaces that do not coincide on the interfaces. Again, we have analyzed both the dependence on the grid size and the spline degree. If the vertex values are chosen as primal degrees of freedom (Algorithms A and C), the results are the same as for conforming discretizations. If we use only the edge averages (Algorithm B), the condition number estimate additionally depends on the ratio between the grid sizes of neighboring patches. This can also be observed in the numerical experiments.

Algorithm B does not perform as good as the other options. However, this approach seems to be beneficial if a non-conforming decomposition of the overall domain into patches is considered, like a decomposition with T-junctions. Although the IETI-DP methods also perform well on domains with non-trivial geometry functions, analyzing the dependence of the condition number on the geometry function is an interesting topic for future research.

Funding source: Austrian Science Fund

Award Identifier / Grant number: S117

Award Identifier / Grant number: W1214-04

Award Identifier / Grant number: P31048

Funding source: Österreichische Austauschdienst

Award Identifier / Grant number: WTZ BG 03/2019

Funding source: Bulgarian National Science Fund

Award Identifier / Grant number: KP-06-Austria/8/2019

Funding statement: The first author was supported by the Austrian Science Fund (FWF): S117 and W1214-04. Also, the second author has received support from the Austrian Science Fund (FWF): P31048. This work was also supported by the bilateral project WTZ BG 03/2019 (KP-06-Austria/8/2019), funded by OeAD (Austria) and Bulgarian National Science Fund.


Moreover, the authors want to thank Ulrich Langer for fruitful discussions and help with the study of existing literature. Finally, the authors thank the anonymous reviewers. Their comments have helped to significantly improve the presentation of this paper.


[1] R. A. Adams and J. J. F. Fournier, Sobolev Spaces, Elsevier/Academic, Amsterdam, 2003. Search in Google Scholar

[2] M. Ainsworth, A preconditioner based on domain decomposition for ℎ-𝑝 finite-element approximation on quasi-uniform meshes, SIAM J. Numer. Anal. 33 (1996), no. 4, 1358–1376. Search in Google Scholar

[3] P. F. Antonietti and B. Ayuso, Schwarz domain decomposition preconditioners for discontinuous Galerkin approximations of elliptic problems: Non-overlapping case, ESAIM Math. Model. Numer. Anal. 41 (2007), no. 1, 21–54. Search in Google Scholar

[4] P. F. Antonietti and P. Houston, A class of domain decomposition preconditioners for h⁢p-discontinuous Galerkin finite element methods, J. Sci. Comput. 46 (2011), no. 1, 124–149. Search in Google Scholar

[5] D. N. Arnold, An interior penalty finite element method with discontinuous elements, SIAM J. Numer. Anal. 19 (1982), no. 4, 742–760. Search in Google Scholar

[6] D. N. Arnold, F. Brezzi, B. Cockburn and L. D. Marini, Unified analysis of discontinuous Galerkin methods for elliptic problems, SIAM J. Numer. Anal. 39 (2001/02), no. 5, 1749–1779. Search in Google Scholar

[7] Y. Bazilevs, L. Beirão da Veiga, J. A. Cottrell, T. J. R. Hughes and G. Sangalli, Isogeometric analysis: Approximation, stability and error estimates for ℎ-refined meshes, Math. Models Methods Appl. Sci. 16 (2006), no. 7, 1031–1090. Search in Google Scholar

[8] C. Canuto, L. F. Pavarino and A. B. Pieri, BDDC preconditioners for continuous and discontinuous Galerkin methods using spectral/h⁢p elements with variable local polynomial degree, IMA J. Numer. Anal. 34 (2014), no. 3, 879–903. Search in Google Scholar

[9] J. A. Cottrell, T. J. R. Hughes and Y. Bazilevs, Isogeometric Analysis–Toward Integration of CAD and FEA, John Wiley & Sons, Chichester, 2009. Search in Google Scholar

[10] L. Diosady and D. Darmofal, BDDC for higher-order discontinuous Galerkin discretizations, Domain Decomposition Methods in Science and Engineering 20, Lect. Notes Comput. Sci. Eng. 91, Springer, Heidelberg (2013), 559–567. Search in Google Scholar

[11] M. Dryja, J. Galvis and M. Sarkis, BDDC methods for discontinuous Galerkin discretization of elliptic problems, J. Complexity 23 (2007), no. 4–6, 715–739. Search in Google Scholar

[12] M. Dryja, J. Galvis and M. Sarkis, A FETI-DP preconditioner for a composite finite element and discontinuous Galerkin method, SIAM J. Numer. Anal. 51 (2013), no. 1, 400–422. Search in Google Scholar

[13] J. A. Evans and T. J. R. Hughes, Explicit trace inequalities for isogeometric analysis and parametric hexahedral finite elements, Numer. Math. 123 (2013), no. 2, 259–290. Search in Google Scholar

[14] C. Farhat and F.-X. Roux, A method of finite element tearing and interconnecting and its parallel solution algorithm, Internat. J. Numer. Methods Engrg. 32 (1991), no. 6, 1205–1227. Search in Google Scholar

[15] B. Guo and W. Cao, Additive Schwarz methods for the ℎ-𝑝 version of the finite element method in two dimensions, SIAM J. Sci. Comput. 18 (1997), no. 5, 1267–1288. Search in Google Scholar

[16] C. Hofer, Analysis of discontinuous Galerkin dual-primal isogeometric tearing and interconnecting methods, Math. Models Methods Appl. Sci. 28 (2018), no. 1, 131–158. Search in Google Scholar

[17] C. Hofer and U. Langer, Dual-primal isogeometric tearing and interconnecting solvers for multipatch dG-IgA equations, Comput. Methods Appl. Mech. Engrg. 316 (2017), 2–21. Search in Google Scholar

[18] C. Hofer and U. Langer, Dual-primal isogeometric tearing and interconnecting methods, Contributions to Partial Differential Equations and Applications, Comput. Methods Appl. Sci. 47, Springer, Cham (2019), 273–296. Search in Google Scholar

[19] C. Hofer, U. Langer and I. Toulopoulos, Isogeometric analysis on non-matching segmentation: Discontinuous Galerkin techniques and efficient solvers, J. Appl. Math. Comput. 61 (2019), no. 1–2, 297–336. Search in Google Scholar

[20] T. J. R. Hughes, J. A. Cottrell and Y. Bazilevs, Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement, Comput. Methods Appl. Mech. Engrg. 194 (2005), no. 39–41, 4135–4195. Search in Google Scholar

[21] A. Klawonn, L. F. Pavarino and O. Rheinbach, Spectral element FETI-DP and BDDC preconditioners with multi-element subdomains, Comput. Methods Appl. Mech. Engrg. 198 (2008), no. 3–4, 511–523. Search in Google Scholar

[22] S. K. Kleiss, C. Pechstein, B. Jüttler and S. Tomar, IETI—isogeometric tearing and interconnecting, Comput. Methods Appl. Mech. Engrg. 247/248 (2012), 201–215. Search in Google Scholar

[23] V. G. Korneev and U. Langer, Dirichlet–Dirichlet Domain Decomposition Methods for Elliptic Problems, World Scientific, Hackensack, 2015. Search in Google Scholar

[24] U. Langer, A. Mantzaflaris, S. E. Moore and I. Toulopoulos, Multipatch discontinuous Galerkin isogeometric analysis, Isogeometric Analysis and Applications 2014, Lect. Notes Comput. Sci. Eng. 107, Springer, Cham (2015), 1–32. Search in Google Scholar

[25] U. Langer and I. Toulopoulos, Analysis of multipatch discontinuous Galerkin IgA approximations to elliptic boundary value problems, Comput. Vis. Sci. 17 (2015), no. 5, 217–233. Search in Google Scholar

[26] J. Mandel, C. R. Dohrmann and R. Tezaur, An algebraic theory for primal and dual substructuring methods by constraints, Appl. Numer. Math. 54 (2005), no. 2, 167–193. Search in Google Scholar

[27] A. Mantzaflaris, R. Schneckenleitner, S. Takacs, G+Smo (Geometry plus Simulation modules), , 2020. Search in Google Scholar

[28] L. F. Pavarino, BDDC and FETI-DP preconditioners for spectral element discretizations, Comput. Methods Appl. Mech. Engrg. 196 (2007), no. 8, 1380–1388. Search in Google Scholar

[29] L. F. Pavarino and O. B. Widlund, A polylogarithmic bound for an iterative substructuring method for spectral elements in three dimensions, SIAM J. Numer. Anal. 33 (1996), no. 4, 1303–1335. Search in Google Scholar

[30] C. Pechstein, Finite and Boundary Element Tearing and Interconnecting Solvers for Multiscale Problems, Springer, Heidelberg, 2013. Search in Google Scholar

[31] B. Rivière, Discontinuous Galerkin Methods for Solving Elliptic and Parabolic Equations, Society for Industrial and Applied Mathematics, Philadelphia, 2008. Search in Google Scholar

[32] R. Schneckenleitner and S. Takacs, Condition number bounds for IETI-DP methods that are explicit in ℎ and 𝑝, Math. Models Methods Appl. Sci. 30 (2020), no. 11, 2067–2103. Search in Google Scholar

[33] R. Schneckenleitner and S. Takacs, IETI-DP for conforming multi-patch Isogeometric Analysis in three dimensions, preprint (2021), . Search in Google Scholar

[34] R. Schneckenleitner and S. Takacs, Towards a IETI-DP solver on non-matching multi-patch domains, preprint (2021), . Search in Google Scholar

[35] J. Schöberl and C. Lehrenfeld, Domain decomposition preconditioning for high order hybrid discontinuous Galerkin methods on tetrahedral meshes, Advanced Finite Element Methods and Applications, Lect. Notes Appl. Comput. Mech. 66, Springer, Heidelberg (2013), 27–56. Search in Google Scholar

[36] J. Schöberl, J. M. Melenk, C. G. A. Pechstein and S. C. Zaglmayr, Schwarz preconditioning for high order simplicial finite elements, Domain Decomposition Methods in Science and Engineering 26, Lect. Notes Comput. Sci. Eng. 55, Springer, Berlin (2007), 139–150. Search in Google Scholar

[37] C. Schwab, 𝑝- and h⁢p-Finite Element Methods: Theory and Applications in Solid and Fluid Mechanics, The Clarendon, New York, 1998. Search in Google Scholar

[38] S. Takacs, A quasi-robust discretization error estimate for discontinuous Galerkin isogeometric analysis, preprint (2019), . Search in Google Scholar

[39] A. Toselli and O. Widlund, Domain Decomposition Methods—Algorithms and Theory, Springer, Berlin, 2005. Search in Google Scholar

Received: 2020-10-13
Revised: 2021-07-05
Accepted: 2021-07-13
Published Online: 2021-08-07

© 2021 Walter de Gruyter GmbH, Berlin/Boston

This work is licensed under the Creative Commons Attribution 4.0 International License.