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.
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. . 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. . For discontinuous Galerkin methods, the discretization error decays in the grid size with the same order as for a conforming discretization, cf. . 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 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 . Since the invention of FETI, various FETI-type methods have been developed, cf. [30, 39, 23]. In , 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  and more recently in . The original results, cf. , are only explicit in the ratio between grid and patch sizes. In , 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 , the authors did a numerical study of 3D problems. The experiments indicate a similar condition number bound as established in . 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  and adapted for IETI in , 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 -FEM are given. The publication  shows a bound for general non-overlapping Schwarz preconditioners for dG -FEM that is in the order of . A polylogarithmic bound for a BDDC preconditioner for a hybridizable discontinuous Galerkin discretization has been developed in .
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  and is a continuation of the paper , 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 Algorithm B, where 𝑝 is the spline degree, are the grid sizes and are the patch sizes, 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 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 and . 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., . 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 , 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 with Lipschitz boundary . We are interested in solving the homogeneous Poisson problem: given , find such that
Throughout the paper, we denote by , , and , , the usual Lebesgue and Sobolev spaces, respectively. Moreover, is the subspace of functions whose trace vanishes on . We use the standard scalar products and , norms , , , and seminorm .
2.2 Representation of the Computational Domain
The computational domain Ω is given as the composition of 𝐾 non-overlapping subdomains , i.e., each is open and bounded with Lipschitz boundary,
where denotes the closure of the set 𝑇. We refer to the subdomains as patches. We assume that every patch 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.
There is a constant such that
for and all , where is the diameter of the patch .
This assumption imposes quality criteria on the representation of the geometry. It might be possible to reduce the constant by a reparameterization of the geometry, i.e., by finding another function such that . Moreover, also subdividing a patch into one or more patches might help in reducing the constant .
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.
For any two patch indices , the intersection is either a common edge (with the neighboring vertices), a common vertex, or empty.
This assumption also ensures that the pre-images of the (Dirichlet) boundary consist of whole edges.
For any two neighboring patches and , the common edge is denoted by , and its pre-images by and . We denote the two end-points of by and . We collect the indices of patches sharing an edge in a set
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.
There is a constant such that 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 be a given spline degree. For ease of notation, we assume that it is the same for all patches. For , a 𝑝-open knot vector is the building block for the B-spline basis defined via the Cox–de Boor formula, cf. [9, equations (2.1) and (2.2)]. This basis spans the univariate spline space
Let and be two 𝑝-open knot vectors over . To get a multivariate spline space over , we tensorize the two univariate spline spaces. The transformation of to the physical domain is defined by the pull-back principle. We denote the resulting space by . So we define
where denotes the restriction of 𝑣 to 𝑇 (trace operator).
The basis for the space consists of only those tensor-product basis functions over that vanish on the pre-image of the Dirichlet boundary .
Say, the total number of basis functions of is .
The number denotes the number of basis function that are supported only in the interior of the patch, whereas accounts for the number of basis functions that contribute to the boundary of the patch.
These definitions give rise to an ordered basis
We assume in the following that the grids on each of the patches are quasi-uniform.
There are grid sizes for and a constant such that
holds for all and all red with .
The corresponding grid size on the physical domain is defined via . For any two patches and sharing an edge, we define
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 , i.e.,
Using these function spaces, the dG discretization of the model problem (2.1) is given by the following: find such that
Due to [38, Theorem 8], the parameter 𝛿 can always be chosen independently of the spline degree 𝑝 and mesh sizes such that the bilinear form is bounded and coercive in the dG-norm
Note that 𝛿 depends on the constant 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 . 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 for all with , then the discretization error satisfies
for some constant 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., . As it has been done previously for IETI methods, the local spaces are chosen on a per-patch basis, cf.  and later publications. In the case of dG, the choice is not completely straight-forward. We follow the approach that has been proposed in  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 and . We define the enriched function space
where is the trace of . We write
where , and .
forms a basis of , which is uniquely defined up to the ordering of the basis functions. We denote the basis by .
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 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.
Following the structure of (3.1), we define a basis
for the space as follows. We start with the basis functions from ,
Then the basis functions follow from the neighboring patches.
On the spaces , we define the bilinear forms and and the linear functional analogous to the global objects and by
By discretizing, we obtain the local linear system
The vector is the coefficient vector representing the function .
We subdivide the stiffness matrix and the load vector into blocks
where the first row and the first column correspond to the first basis functions, i.e., to those supported in the interior of the patch . 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),
Once has been computed, we get the solution 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,
where and . A basis for is given by
As above, the coefficient vector represents the function
The Schur complements realize the discrete harmonic extension with respect to the energy norm , which is defined as follows. For any , we have
where is such that
3.3 The Jump Matrices
The next step is the introduction of constraints that yield the continuity of the solution. For any two patches and sharing an edge and any basis function supported on that edge, we introduce a constraint
where is as in (3.2) and and are the coefficients of the functions and . The index shifts by and account for the restriction to and , 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 that are supported on a corner of . 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 that is supported on a corner of , we additionally introduce constraints of the form
where such that the said corner lies between the edges and . These additional constraints also ensure the continuity between the different artificial interfaces. A visualization of this case is given in Figure 4.
We collect the constraints of the form (3.9) and (3.10) in a jump matrix 𝐵 such that the constraints are equivalent to . 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 .
We are now able to write down a saddle point problem in matrix-vector form that is equivalent to problem (2.4): find such that
Note that the matrices and are non-singular only if the corresponding patch contributes to the Dirichlet boundary . Since in general at least one of the matrices 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
Algorithm A (Vertex Values)
The space is the subspace of functions where the vertex values agree, i.e.,
where and are the two end-points of . Furthermore, satisfies these conditions homogeneously, i.e., and
Algorithm B (Edge Averages)
The space is the subspace of functions where the averages of the function values over the edges agree, i.e.,
Furthermore, satisfies these conditions homogeneously, i.e., and
Algorithm C (Vertex Values and Edge Everages)
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.
The primal space is the subspace of energy minimizing functions, cf. . This means that is 𝑆-orthogonal to , i.e.,
Let be a basis of . 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 can be obtained by full-rank matrices , . The block-diagonal collection of the matrices forms a matrix 𝐶.
3.5 The IETI-DP System and the Scaled Dirichlet Preconditioner
The solution for the original problem is obtained by . Following the DP approach, this system can be reduced to
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 , defined by
where is a diagonal matrix defined based on the principle of multiplicity scaling: each coefficient 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 according to (3.5).
Compute the primal Schur complement .
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 requires the following steps:
Solve the (usually small) global linear system .
Compute the residual .
The vector , representing the solution on the skeleton, is computed analogously to steps (2), (3), and (4) based on . Finally, the solution vector is derived by (3.6).
4 The Condition Number Estimate
In this section, we prove the following condition number estimates.
Provided that the IETI-DP solver is set up as outlined in the previous sections, the condition number of the preconditioned system satisfies
If some Lagrange multipliers are redundant to each other or to primal constraints, the matrices 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 , where . 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.
The space is the standard Sobolev–Slobodeckij space with a norm defined as follows. For open, bounded, and Lipschitz and , we define the seminorm by
and the corresponding norm by . As in , 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.
For all patches , the equivalences
for all and , and
for all and
The results for are obtained by Assumption 1, the chain rule and the substitution rule for differentiation, see e.g., [7, Lemma 3.5]. The results for follow by Hilbert space interpolation theory, cf. [1, Theorem 7.23, Theorem 7.31]. ∎
The next lemma states that the bilinear forms and are patchwise equivalent.
For every patch , the relation holds for all .
Using the Cauchy–Schwarz inequality, estimate (4.1), and , we get
for some constant , cf.  for details. With the choice , 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
Analogous to [32, Lemma 4.16], the following lemma estimates the effect of .
with coefficient vector , and let
with coefficient vector be such that .
Then we have, for each patch and each edge connecting the vertices and ,
We start by recalling how the scaling matrix 𝐷 looks like. Remember 𝐷 is a diagonal matrix. We denote the coefficients by for . The entries are defined to be one plus the number of Lagrange multipliers of the corresponding degree of freedom 𝑖. Note that we have unless the corresponding degree of freedom is a corner degree of freedom. In that case, takes a value, which is bounded from below by 1 and bounded from above due to Assumption 3. Thus, in any case.
We will start with Algorithms A and C: simple calculations reveal for Algorithms A and C that
for the -norm. For the -seminorm, we get using the triangle inequality that
Analogously, we obtain, using the triangle inequality,
For Algorithm B, we have that on can be expanded as
Certainly, we also have . 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 -seminorm on the artificial edges.
For any two patches and that share an edge ,
holds for all .
Let be arbitrary but fixed. Using the triangle inequality, we have
Using [1, Theorem 5.2, equation (3)] and the triangle inequality, we obtain further
where is an appropriately chosen constant (that only depends on the constant from Assumption 1). Using , we have
By subtracting , we immediately obtain
for all . Since and since the bound for 𝒜 holds for all , the Poincaré inequality yields the desired result. ∎
The next step is to show that a similar estimate holds for the seminorm . Before we can prove that result, we need the following auxiliary result.
holds for all .
Since 𝑢 is continuous, 𝑢 takes its maximum for some . Then, by the fundamental theorem of differential and integral calculus, we can write . The desired result is obtained by integrating over the unit integral and by utilizing the Cauchy–Schwarz inequality. ∎
For any two patches and that share an edge ,
holds for all .
Using the triangle inequality, we obtain
As a next step, we apply Lemma 5 to the difference . Since the -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
Since this holds for all , the Poincaré inequality gives further
Lemma 4 yields the final result. ∎
For any two patches and , sharing an edge which connects the vertices and and and , we have
for all , where and is as in (4.3).
Using the triangle inequality, we have
and observe that the Cauchy–Schwarz inequality and Assumption 1 yield
Now, we estimate . Thus, we set and . With unitary transformations, i.e., rotations and reflections, we can transform the patches such that the pre-image of is and that the pre-image of is 0. Let and . Let 𝜂 be the largest value such that and are polynomial on . Using the quasi-uniformity assumption, cf. Assumption 4, we obtain . 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 -norm on is larger than the -norm on and ,
An application of Lemma 1 finally yields
4.4 Proof of the Main Theorem
Before we give a proof of the main theorem, we estimate the sum of -seminorms over all patches.
Let 𝑢 and 𝑤 be as in Lemma 3. Then we have
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 . ∎
Let 𝑢 and 𝑤 be as in Lemma 3. Then we have
Lemma 3 and the observation that yield
Then and 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  states that
where is the coefficient vector associated to the function . So let 𝑢 be arbitrary but fixed, and let the function with coefficient vector be such that . Let be an arbitrary extension of . Lemma 2 yields
Note that the second sum does not depend on . Thus, the infimum refers only to the -seminorm, which means that the seminorm coincides with the seminorm of the (standard) discrete harmonic extension. Hence, we have
Using and Lemma 8, we obtain further
Using [32, Lemma 4.15 and Theorem 4.2] and , we further estimate
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
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 , we insert one knot into each knot span. The new knot is not located in the center, but at times the length of the knot span if the patch index 𝑘 is even and at times that length if 𝑘 is odd. The subsequent refinement steps 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 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 is an estimate for the size of the smaller side of the smallest element of the physical patches and , and for Algorithms A and C and for Algorithm B. To solve system (3.13), we use a preconditioned conjugate gradient (PCG) method with the proposed scaled Dirichlet preconditioner . The implementation is done using G+Smo . The local subproblems are solved with the sparse direct solver from the PARDISO project.
We start the iteration with a randomly sampled vector with entries in the interval . The stopping criterion is chosen as follows. The iteration stops if the -norm of the residual vector drops below 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 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 . 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 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 in the grid size.
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 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 , and the increase in 𝑝 is sub-linear, almost like . 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 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
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 , where 𝑘 is even, additionally times. Thus, the grid sizes between patches with even and odd degrees vary by a factor of . 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 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.
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 on , where
and observe that 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 -norm of the residual vector drops below times the residual of the right-hand side vector. In Table 9 and Figure 10, we present the discretization error and the corresponding rate , which are defined by
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.
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 patches as initial configuration. This serves as the computational domain with the fewest number of patches. The decompositions into , 192, and 768 patches are obtained by splitting the patches uniformly. On each patch, we consider spline spaces that are obtained by 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.
We have extended the theory from , 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.
 R. A. Adams and J. J. F. Fournier, Sobolev Spaces, Elsevier/Academic, Amsterdam, 2003. Search in Google Scholar
 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
 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
 P. F. Antonietti and P. Houston, A class of domain decomposition preconditioners for -discontinuous Galerkin finite element methods, J. Sci. Comput. 46 (2011), no. 1, 124–149. Search in Google Scholar
 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
 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
 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
 C. Canuto, L. F. Pavarino and A. B. Pieri, BDDC preconditioners for continuous and discontinuous Galerkin methods using spectral/ elements with variable local polynomial degree, IMA J. Numer. Anal. 34 (2014), no. 3, 879–903. Search in Google Scholar
 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
 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
 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
 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
 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
 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
 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
 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
 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
 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
 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
 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
 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
 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
 V. G. Korneev and U. Langer, Dirichlet–Dirichlet Domain Decomposition Methods for Elliptic Problems, World Scientific, Hackensack, 2015. Search in Google Scholar
 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
 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
 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
 A. Mantzaflaris, R. Schneckenleitner, S. Takacs, G+Smo (Geometry plus Simulation modules), , 2020. Search in Google Scholar
 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
 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
 C. Pechstein, Finite and Boundary Element Tearing and Interconnecting Solvers for Multiscale Problems, Springer, Heidelberg, 2013. Search in Google Scholar
 B. Rivière, Discontinuous Galerkin Methods for Solving Elliptic and Parabolic Equations, Society for Industrial and Applied Mathematics, Philadelphia, 2008. Search in Google Scholar
 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
 R. Schneckenleitner and S. Takacs, IETI-DP for conforming multi-patch Isogeometric Analysis in three dimensions, preprint (2021), . Search in Google Scholar
 R. Schneckenleitner and S. Takacs, Towards a IETI-DP solver on non-matching multi-patch domains, preprint (2021), . Search in Google Scholar
 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
 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
 C. Schwab, 𝑝- and -Finite Element Methods: Theory and Applications in Solid and Fluid Mechanics, The Clarendon, New York, 1998. Search in Google Scholar
 S. Takacs, A quasi-robust discretization error estimate for discontinuous Galerkin isogeometric analysis, preprint (2019), . Search in Google Scholar
 A. Toselli and O. Widlund, Domain Decomposition Methods—Algorithms and Theory, Springer, Berlin, 2005. Search in Google Scholar
© 2021 Walter de Gruyter GmbH, Berlin/Boston
This work is licensed under the Creative Commons Attribution 4.0 International License.