Skip to content
BY 4.0 license Open Access Published by De Gruyter August 7, 2021

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

  • Rainer Schneckenleitner EMAIL logo and Stefan Takacs

Abstract

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 p 2 . 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

C ⁢ p ⁢ ( 1 + log ⁡ p + max k = 1 , … , K ⁡ log ⁡ H k h k ) 2

for Algorithms A and C and by

C ⁢ δ ⁢ p ⁢ ( max k = 1 , … , K ⁡ max ℓ ∈ N Γ ⁢ ( k ) ⁡ h k h ℓ ) ⁢ ( 1 + log ⁡ p + max k = 1 , … , K ⁡ log ⁡ H k h k ) 2

for Algorithm B, where 𝑝 is the spline degree, h k are the grid sizes and H k 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 C 0 and C p - 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 Ω ⊂ R 2 with Lipschitz boundary ∂ ⁡ Ω . We are interested in solving the homogeneous Poisson problem: given f ∈ L 2 ⁢ ( Ω ) , find u ∈ H 0 1 ⁢ ( Ω ) such that

(2.1) ∫ Ω ∇ ⁡ u ⋅ ∇ ⁡ v ⁢ d ⁢ x = ∫ Ω f ⁢ v ⁢ d x   for all ⁢ v ∈ H 0 1 ⁢ ( Ω ) ,

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

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 = 1 K Ω ( 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

G k : Ω ^ := ( 0 , 1 ) 2 → Ω ( k ) := G k ⁢ ( Ω ^ ) ⊂ R 2 ,

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 C 1 > 0 such that

∥ ∇ ⁡ G k ∥ L ∞ ⁢ ( T ) ≤ C 1 ⁢ H k   and   ∥ ( ∇ ⁡ G k ) - 1 ∥ L ∞ ⁢ ( T ) ≤ C 1 ⁢ 1 H k

for T ∈ { Ω ^ , ∂ ⁡ Ω ^ } and all k = 1 , … , K , where H k > 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 C 1 by a reparameterization of the geometry, i.e., by finding another function G ~ k such that G ~ k ⁢ ( Ω ^ ) = G k ⁢ ( Ω ^ ) . Moreover, also subdividing a patch into one or more patches might help in reducing the constant C 1 .

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 ) := G k - 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 , ℓ ) := G k - 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

P ⁢ ( x ) := { k : x ∈ Ω ( k ) ¯ } .

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

Assumption 3

There is a constant C 2 such that | P ⁢ ( x ) | ≤ C 2 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 = 1 n defined via the Cox–de Boor formula, cf. [9, equations (2.1) and (2.2)]. This basis spans the univariate spline space

S ⁢ [ p , Ξ ] := span ⁡ { B ⁢ [ p , Ξ , 1 ] , … , B ⁢ [ p , Ξ , n ] } .

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 ) ∘ G k - 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 ) = N I ( k ) + N Γ ( k ) . The number N I ( 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)

Φ ^ ( k ) := ( ϕ ^ i ( k ) ) i = 1 N ( k ) ,
{ ϕ ^ i ( k ) } = { ϕ ^ ∈ H 1 ( Ω ^ ) : ϕ ^ ( x , y ) = B [ p , Ξ ( k , 1 ) , j 1 ] ( x ) B [ p , Ξ ( k , 2 ) , j 2 ] ( y ) ,
where j 1 ∈ { 1 , … , n ( k , 1 ) } and j 2 ∈ { 1 , … , n ( k , 2 ) } such that ϕ ^ | Γ ^ D ( k ) = 0 } ,
ϕ ^ i ( k ) | ∂ ⁡ Ω ^ = 0 ⇔ i ∈ { 1 , … , N I ( k ) }   and   ϕ ^ i ( k ) | ∂ ⁡ Ω ^ ≠ 0 ⇔ i ∈ N I ( k ) + { 1 , … , N Γ ( k ) } .
The pull-back principle gives a basis for V ( k ) ,

Φ ( k ) := ( ϕ i ( k ) ) i = 1 N ( k )   and   ϕ i ( k ) := ϕ ^ i ( k ) ∘ G k - 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 C 3 > 0 such that

C 3 ⁢ h ^ k ≤ ξ i + 1 ( k , δ ) - ξ i ( k , δ ) ≤ h ^ k

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 h k := h ^ k ⁢ H k . For any two patches Ω ( k ) and Ω ( ℓ ) sharing an edge, we define

h ^ k ⁢ ℓ := min ⁡ { h ^ k , h ^ ℓ }   and   h k ⁢ ℓ := min ⁡ { h k , 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.,

V := ∏ k = 1 K V ( k ) := V ( 1 ) × ⋯ × V ( K ) .

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) a h ⁢ ( u , v ) = ∫ Ω f ⁢ v ⁢ d x   for all ⁢ v ∈ V ,

where

a h ⁢ ( u , v ) := ∑ k = 1 K ( a ( k ) ⁢ ( u , v ) + m ( k ) ⁢ ( u , v ) + r ( k ) ⁢ ( u , v ) ) ,
a ( k ) ⁢ ( u , v ) := ∫ Ω ( k ) ∇ ⁡ u ( k ) ⋅ ∇ ⁡ v ( k ) ⁢ d ⁢ x ,
m ( k ) ⁢ ( u , v ) := ∑ ℓ ∈ N Γ ⁢ ( k ) ∫ Γ ( k , ℓ ) 1 2 ⁢ ( ∂ ⁡ u ( k ) ∂ ⁡ n k ⁢ ( v ( ℓ ) - v ( k ) ) + ∂ ⁡ v ( k ) ∂ ⁡ n k ⁢ ( u ( ℓ ) - u ( k ) ) ) ⁢ d s ,
r ( k ) ⁢ ( u , v ) := ∑ ℓ ∈ N Γ ⁢ ( k ) ∫ Γ ( k , ℓ ) δ ⁢ p 2 h k ⁢ ℓ ⁢ ( u ( ℓ ) - u ( k ) ) ⁢ ( v ( ℓ ) - v ( k ) ) ⁢ d s ,
and δ > 0 is a penalty parameter and n k 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 h k such that the bilinear form a h ⁢ ( ⋅ , ⋅ ) is bounded and coercive in the dG-norm

∥ v ∥ d 2 := d ⁢ ( v , v ) , where ⁢ d ⁢ ( u , v ) := ∑ k = 1 K ( a ( k ) ⁢ ( u , v ) + r ( k ) ⁢ ( u , v ) ) .

Note that 𝛿 depends on the constant C 1 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 ) ∈ H q + 1 ⁢ ( Ω ( k ) ) for all k = 1 , … , K with 1 ≤ q ≤ p , then the discretization error satisfies

(2.5) ∥ u - u h ∥ d 2 ≤ C ⁢ h 2 ⁢ q ⁢ ∑ k = 1 K ∥ u ∥ H q + 1 ⁢ ( Ω ( k ) ) 2

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

V e ( k ) := V ( k ) × ∏ ℓ ∈ N Γ ⁢ ( k ) V ( k , ℓ ) ,

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

(3.1) v e ( k ) = ( v ( k ) , ( v ( k , ℓ ) ) ℓ ∈ N Γ ⁢ ( k ) ) ,

where v e ( k ) ∈ V e ( k ) , v ( k ) ∈ V ( k ) and v ( k , ℓ ) ∈ V ( k , ℓ ) .

Note that

{ ϕ 1 ( k , ℓ ) , … , ϕ N ( k , ℓ ) ( k , ℓ ) } := { ϕ i ( ℓ ) | Γ ( k , ℓ ) : ϕ i ( ℓ ) ∈ Φ ( ℓ ) , ϕ i ( ℓ ) | Γ ( k , ℓ ) ≠ 0 }

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 = 1 N ( 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 = 1 N e ( k )   with ⁢ N e ( k ) := dim ⁡ V e ( k )

for the space V e ( 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)

ϕ e , Υ i ( k , ℓ 1 ) ( k ) := ( 0 , ( ϕ i ( k , ℓ 1 ) , 0 ⁢ … , 0 ) )
with ⁢ Υ i ( k , ℓ 1 ) := N ( k ) + i for ⁢ i = 1 , … , N ( k , ℓ 1 ) ,
ϕ e , Υ i ( k , ℓ 2 ) ( k ) := ( 0 , ( 0 , ϕ i ( k , ℓ 2 ) , 0 ⁢ … , 0 ) )
with ⁢ Υ i ( k , ℓ 2 ) := N ( k ) + N ( k , ℓ 1 ) + i for ⁢ i = 1 , … , N ( k , ℓ 2 ) ,
⋮
ϕ e , Υ i ( k , ℓ L ) ( k ) := ( 0 , ( 0 , … , 0 , ϕ i ( k , ℓ L ) ) )
with ⁢ Υ i ( k , ℓ L ) := N ( k ) + N ( k , ℓ 1 ) + ⋯ + N ( k , ℓ L - 1 ) + i for ⁢ i = 1 , … , N ( k , ℓ L ) .
From this construction and from (2.3), we know that the first N I ( 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 V e ( k ) , we define the bilinear forms a e ( k ) ⁢ ( ⋅ , ⋅ ) and d e ( k ) ⁢ ( ⋅ , ⋅ ) and the linear functional ⟨ f e ( k ) , ⋅ ⟩ analogous to the global objects a h ⁢ ( ⋅ , ⋅ ) and d ⁢ ( ⋅ , ⋅ ) by

a e ( k ) ⁢ ( u e ( k ) , v e ( k ) ) := a ( k ) ⁢ ( u e ( k ) , v e ( k ) ) + m ( k ) ⁢ ( u e ( k ) , v e ( k ) ) + r ( k ) ⁢ ( u e ( k ) , v e ( k ) ) ,
d e ( k ) ⁢ ( u e ( k ) , v e ( k ) ) := a ( k ) ⁢ ( u e ( k ) , v e ( k ) ) + r ( k ) ⁢ ( u e ( k ) , v e ( k ) ) ,
⟨ f e ( k ) , v e ( k ) ⟩ := ∫ Ω ( k ) f ⁢ v ( k ) ⁢ d x ,
where we write with a slight abuse of notation
a ( k ) ⁢ ( u e ( k ) , v e ( k ) ) := ∫ Ω ( k ) ∇ ⁡ u ( k ) ⋅ ∇ ⁡ v ( k ) ⁢ d ⁢ x ,
m ( k ) ⁢ ( u e ( k ) , v e ( k ) ) := ∑ ℓ ∈ N Γ ⁢ ( k ) ∫ Γ ( k , ℓ ) 1 2 ⁢ ( ∂ ⁡ u ( k ) ∂ ⁡ n k ⁢ ( v ( k , ℓ ) - v ( k ) ) + ∂ ⁡ v ( k ) ∂ ⁡ n k ⁢ ( u ( k , ℓ ) - u ( k ) ) ) ⁢ d s ,
r ( k ) ⁢ ( u e ( k ) , v e ( k ) ) := ∑ ℓ ∈ N Γ ⁢ ( k ) ∫ Γ ( k , ℓ ) δ ⁢ p 2 h k ⁢ ℓ ⁢ ( u ( k , ℓ ) - u ( k ) ) ⁢ ( v ( k , ℓ ) - v ( k ) ) ⁢ d s .
The bilinear form d e ( k ) ⁢ ( ⋅ , ⋅ ) induces the local dG-norm ∥ ⋅ ∥ d e ( k ) 2 := d e ( k ) ⁢ ( ⋅ , ⋅ ) , the bilinear form a e ( k ) ⁢ ( ⋅ , ⋅ ) the local energy norm ∥ ⋅ ∥ a e ( k ) 2 := a e ( k ) ⁢ ( ⋅ , ⋅ ) .

By discretizing, we obtain the local linear system

(3.3) A ( k ) ⁢ u ¯ e ( k ) = f ¯ e ( k ) ,

where

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

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

We subdivide the stiffness matrix and the load vector into blocks

A ( k ) = ( A II ( k ) A I ⁢ Γ ( 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 N I ( 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.4) S ( k ) ⁢ w ¯ e ( k ) = g ¯ e ( k ) ,

where

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

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

(3.6) u ¯ e ( k ) = ( u ¯ I ( k ) u ¯ Γ ( k ) ) = ( ( A II ( k ) ) - 1 ⁢ ( f ¯ I ( k ) - A I ⁢ Γ ( k ) ⁢ w ¯ e ( k ) ) w ¯ e ( k ) ) .

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

(3.7) S ⁢ w ¯ = g ¯ .

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 = 1 K W e ( k ) , where ⁢ W e ( k ) := W ( k ) × ∏ ℓ ∈ N Γ ⁢ ( k ) W ( k , ℓ ) ,

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

w e ( k ) = ( w ( k ) , ( w ( k , ℓ ) ) ℓ ∈ N Γ ⁢ ( k ) ) ,

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

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

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

(3.8) w e ( k ) = ∑ i = 1 N ˘ e ( k ) w e , i ( k ) ⁢ ϕ e , N I ( k ) + i ( k ) .

The Schur complements realize the discrete harmonic extension H A ( k ) : W e ( k ) → V e ( k ) with respect to the energy norm ∥ ⋅ ∥ a e ( k ) , which is defined as follows. For any w e ( k ) = ( w ( k ) , ( w ( k , ℓ ) ) ℓ ∈ N Γ ⁢ ( k ) ) , we have

H A ( k ) ⁢ w e ( k ) = ( u ( k ) , ( w ( k , ℓ ) ) ℓ ∈ N Γ ⁢ ( k ) ) ,

where u ( k ) is such that

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

∑ k = 1 K a e ( k ) ⁢ ( H A ( k ) ⁢ w e ( k ) , H A ( k ) ⁢ q e ( k ) ) ⏟ = ⁣ : s ⁢ ( w , q ) = ∑ k = 1 K ⟨ f , H A ( k ) ⁢ q e ( k ) ⟩ ⏟ = ⁣ : ⟨ g , q ⟩   for all ⁢ q = ( q e ( 1 ) , … , q e ( 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) w e , i * ( k ) - w e , j * ( ℓ ) = 0   with ⁢ i * = i - N I ( k ) ⁢ and ⁢ j * = Υ i ( ℓ , k ) - N I ( ℓ ) ,

where Υ i ( ℓ , k ) is as in (3.2) and w e , i * ( k ) and w e , j * ( k ) are the coefficients of the functions w e ( k ) and w e ( ℓ ) . The index shifts by N I ( k ) and N I ( ℓ ) account for the restriction to W e ( k ) and W e ( ℓ ) , 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) w e , j 1 * ( ℓ 1 ) - w e , j 2 * ( ℓ 2 ) = 0   with ⁢ j 1 * = Υ i ( ℓ 1 , k ) - N I ( ℓ 1 ) ⁢ and ⁢ j 2 * = Υ i ( ℓ 2 , k ) - N I ( ℓ 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

( S B ⊤ B ) ⁢ ( w ¯ λ ¯ ) = ( g ¯ 0 ) .

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 = 1 K W ~ Δ ( k ) and

W ~ Δ ( k ) := { w e ( k ) ∈ W e ( 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 ) ⁢ d x = ∫ Γ ( k , ℓ ) w ( ℓ , k ) ⁢ d x ⁢ for all ⁢ k ⁢ and all ⁢ ℓ ∈ N Γ ⁢ ( k ) } .

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

W ~ Δ ( k ) := { w e ( k ) ∈ W e ( k ) : ∫ Γ ( k , ℓ ) w ( k ) ⁢ d x = ∫ Γ ( k , ℓ ) w ( ℓ , k ) ⁢ d x = 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.,

(3.12) Ψ ⊤ = ( ( ψ ¯ ( 1 ) ) ⊤ ⁢ ⋯ ⁢ ( ψ ¯ ( N Π ) ) ⊤ ) ⊤ .

A representation for the spaces W ~ Δ ( k ) can be obtained by full-rank matrices C ( k ) , C ( k ) ⁢ w ¯ e ( k ) = 0 ⇔ w e ( 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

( S C ⊤ 0 B ⊤ C 0 0 0 0 0 Ψ ⊤ ⁢ S ⁢ Ψ Ψ ⊤ ⁢ B ⊤ B 0 B ⁢ Ψ 0 ) ⁢ ( w ¯ Δ μ ¯ w ¯ Π λ ¯ ) = ( g ¯ 0 Ψ ⊤ ⁢ g ¯ 0 ) .

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

(3.13) F ⁢ λ ¯ = d ¯ ,

where

F := ( B ⁢  0 ⁢ B ⁢ Ψ ) ⁢ ( S C ⊤ 0 C 0 0 0 0 Ψ ⊤ ⁢ S ⁢ Ψ ) - 1 ⁢ ( I 0 Ψ ⊤ ) ⏟ = ⁣ : F 0 ⁢ B ⊤   and   d ¯ := F 0 ⁢ 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 M sD , defined by

M sD := B ⁢ D - 1 ⁢ S ⁢ D - 1 ⁢ B ⊤ ,

where D ∈ R N Γ × N Γ is a diagonal matrix defined based on the principle of multiplicity scaling: each coefficient d i , 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

      ( S ( k ) ( C ( k ) ) ⊤ C ( k ) ) ⁢ ( w ¯ ^ Δ ( k ) μ ¯ ^ ( k ) ) = ( q ¯ ^ e ( k ) 0 )

      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 M sD ⁢ 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

κ ⁢ ( M sD ⁢ F ) ≤ C ⁢ p ⁢ ( 1 + log ⁡ p + max k = 1 , … , K ⁡ log ⁡ H k h k ) 2

for Algorithms A and C and

κ ⁢ ( M sD ⁢ F ) ≤ C ⁢ δ ⁢ p ⁢ ( max k = 1 , … , K ⁡ max ℓ ∈ N Γ ⁢ ( k ) ⁡ h k h ℓ ) ⁢ ( 1 + log ⁡ p + max k = 1 , … , K ⁡ log ⁡ H k h k ) 2

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 M sD 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 B D ⊤ ⁢ B ⁢ u ¯ , where B D := 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 H 1 / 2 is the standard Sobolev–Slobodeckij space with a norm defined as follows. For T ⊂ R 2 open, bounded, and Lipschitz and Γ ⊆ ∂ ⁡ T , we define the seminorm | ⋅ | H 1 / 2 ⁢ ( Γ ) 2 by

| u | H 1 / 2 ⁢ ( Γ ) 2 := ∫ Γ ∫ Γ | u ⁢ ( x ) - u ⁢ ( y ) | 2 | x - y | 2 ⁢ d x ⁢ d y

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

| v | L ∞ 0 ⁢ ( T ) := inf c ∈ R ⁡ ∥ v - c ∥ L ∞ ⁢ ( T ) .

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 | H s ⁢ ( Ω ( k ) ) ≂ ( H k ) 1 - s ⁢ | u ∘ G k | H s ⁢ ( Ω ^ ) for all u ∈ H s ⁢ ( Ω ( k ) ) and s ∈ { 0 , 1 } , and

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

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

Proof

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 a e ( k ) ⁢ ( ⋅ , ⋅ ) and d e ( k ) ⁢ ( ⋅ , ⋅ ) are patchwise equivalent.

Lemma 2

For every patch Ω ( k ) , the relation a e ( k ) ⁢ ( u e ( k ) , u e ( k ) ) ≲ d e ( k ) ⁢ ( u e ( k ) , u e ( k ) ) ≲ a e ( k ) ⁢ ( u e ( k ) , u e ( k ) ) holds for all u e ( k ) ∈ V e ( k ) .

Proof

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

(4.1) ∥ ∇ ⁡ u ( k ) ⋅ n k ∥ L 2 ⁢ ( Γ ( k , ℓ ) ) 2 ≲ ∥ ∇ ⁡ u ( k ) ∥ L 2 ⁢ ( Γ ( k , ℓ ) ) 2 ⁢ ∥ n k ∥ L ∞ ⁢ ( Γ ( k , ℓ ) ) 2 ≲ ∥ ∇ ⁡ ( u ( k ) ∘ G k ) ∥ L 2 ⁢ ( Γ ^ ( k , ℓ ) ) 2 ≲ p 2 h ^ k ⁢ ℓ ⁢ | u ( k ) ∘ G k | H 1 ⁢ ( Ω ^ ) 2 ≲ p 2 h k ⁢ ℓ ⁢ | u ( k ) | H 1 ⁢ ( Ω ( k ) ) 2 .

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

(4.2) | m ( k ) ⁢ ( u e ( k ) , v e ( k ) ) | ≤ 2 ⁢ c δ 1 / 2 ⁢ ∥ u e ( k ) ∥ d e ( k ) ⁢ ∥ v e ( k ) ∥ d e ( k )

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

a e ( k ) ⁢ ( u e ( k ) , u e ( k ) ) = ∥ u e ( k ) ∥ d e ( k ) 2 + m ( k ) ⁢ ( u e ( k ) , u e ( k ) ) ≥ ∥ u e ( k ) ∥ d e ( k ) 2 - | m ( k ) ⁢ ( u e ( k ) , u e ( k ) ) | ≥ 1 2 ⁢ ∥ u e ( k ) ∥ d e ( k ) 2 ,

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 B D ⊤ ⁢ B

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

Lemma 3

Let

u = ( u e ( 1 ) , … , u e ( K ) ) = ( ( u ( 1 ) , ( u ( 1 , ℓ ) ) ℓ ∈ N Γ ⁢ ( 1 ) ) , … ) ∈ W ~

with coefficient vector u ¯ , and let

w = ( w e ( 1 ) , … , w e ( K ) ) = ( ( w ( 1 ) , ( w ( 1 , ℓ ) ) ℓ ∈ N Γ ⁢ ( 1 ) ) , … ) ∈ W ~

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

∥ w ( k ) - w ( k , ℓ ) ∥ L 2 ⁢ ( Γ ( k , ℓ ) ) 2 ≲ ∥ u ( k ) - u ( k , ℓ ) ∥ L 2 ⁢ ( Γ ( k , ℓ ) ) 2 + ∥ u ( ℓ ) - u ( ℓ , k ) ∥ L 2 ⁢ ( Γ ( k , ℓ ) ) 2 + ∑ i = 1 2 ( h k p ⁢ Δ ( k , ℓ , i ) + h ℓ p ⁢ Δ ( ℓ , k , i ) ) ,
| w ( k ) | H 1 / 2 ⁢ ( Γ ( k , ℓ ) ) 2 ≲ | u ( k ) | H 1 / 2 ⁢ ( Γ ( k , ℓ ) ) 2 + | u ( ℓ , k ) | H 1 / 2 ⁢ ( Γ ( k , ℓ ) ) 2 + ∑ i = 1 2 Δ ( k , ℓ , i ) ,
| w ( k ) | L ∞ 0 ⁢ ( Γ ( k , ℓ ) ) 2 ≲ | u ( k ) | L ∞ 0 ⁢ ( Γ ( k , ℓ ) ) 2 + | u ( ℓ , k ) | L ∞ 0 ⁢ ( Γ ( k , ℓ ) ) 2 + ∑ i = 1 2 Δ ( k , ℓ , i ) ,
where

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

Proof

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

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

w ( k ) | Γ ( k , ℓ ) = 1 2 ⁢ ( u ( k ) | Γ ( k , ℓ ) - u ( ℓ , k ) - ∑ i = 1 2 θ ( k , ℓ , i ) ⁢ ( u ( k ) ⁢ ( x ( k , ℓ , i ) ) -