JOHANNES KEPLER UNIVERSITY LINZ Institute of Computational Mathematics IETI-DP Methods for Discontinuous Galerkin Multi-Patch Isogeometric Analysis with T-junctions

We study Dual-Primal Isogeometric Tearing and Interconnecting (IETI-DP) solvers for non-conforming multi-patch discretizations of a generalized Poisson problem. We realize the coupling between the patches using a symmetric interior penalty discontinuous Galerkin (SIPG) approach. Previously, we have assumed that the interfaces between patches always consist of whole edges. In this paper, we drop this requirement and allow T-junctions. This extension is vital for the consideration of sliding interfaces, for example between the rotor and the stator of an electrical motor. One critical part for the handling of T-junctions in IETI-DP solvers is the choice of the primal degrees of freedom. We propose to add all basis functions that are non-zero at any of the vertices to the primal space. Since there are several such basis functions at any T-junction, we call this concept “fat vertices”. For this choice, we show a condition number bound that coincides with the bound for the conforming case.


Introduction
Isogeometric Analysis (IgA), [17], is an approach to discretize partial differential equations (PDEs) that has been designed in order to overcome difficulties related to meshing of the computational domain.In IgA, the computational domain is parameterized by geometry functions, which are commonly represented in terms of B-splines or nonuniform rational B-splines (NURBS).Such a representation is also used in state-of-theart computer aided-design (CAD) software.Usually, one considers multiple patches, each parameterized with its own geometry function (multi-patch IgA).We consider the case of non-overlapping patches.

2
Conforming discretizations require that both the geometry functions and the grids agree on each interface between two patches.If this is not the case, discontinuous Galerkin (dG) methods, particularly the symmetric interior penalty discontinuous Galerkin (SIPG) approach [3], are an appropriate option.For its adaptation to IgA, see [13,15,16,26] and others.In these publications, it is assumed that the interfaces between two patches consist (in the two-dimensional case) of whole edges, which excludes the case of T-junctions between patches.Now, we include the case of Tjunctions, which allows greater flexibility for the geometry modeling.This is of vital interest for the simulation of objects with sliding interfaces, like the interface between the rotor and the stator of an electrical motor, which serves as computational domain in our model problem.The PDE in the model problem is the Poisson equation, which can be motivated as a model for the magnetostatic potential.In general, sliding interfaces lead to a non-matching decomposition of the computational domain into patches.In two dimensions, this means that T-junctions between patches occur for most rotational angles.So far, several approaches have been considered to handle such types of problems, including the more classical locked-step methods, cf.[24], the moving band technique, cf.[7], the Lagrange multiplier method, cf.[20], interpolation approaches, cf.[23].More recently, mortar techniques, cf.[4,9], domain interface methods, cf.[5] or discontinuous Galerkin methods, cf.[2] have been considered.Also combinations of different approaches have been proposed, see, e.g., [18].
We focus on the SIPG approach, the contribution of this paper is a fast solver for the linear system obtained from the proposed discretization.A canonical choice for domains with many non-overlapping patches are domain decomposition (DD) solvers, like the Dual-Primal Finite Element Tearing and Interconnecting (FETI-DP) method, cf.[10,11].These solvers have been adapted to IgA in [19] and are sometimes referred to as Dual-Primal Isogeometric Tearing and Interconnecting (IETI-DP) solvers.Recently, the authors have developed an analysis that is also explicit in the spline degree, see [25].
The IETI-DP solvers have been extended to discontinuous Galerkin discretizations in [13,14,15,26], however T-junctions have not been covered by the analysis so far.Most components of the IETI-DP framework can easily be extended to domains with T-junctions, see [27] for a numerical study.One of the critical questions is the choice of the primal degrees of freedom.We propose to add all basis functions that are non-zero at a vertex to the primal space.This yields a number of basis functions for each T-junction that grows linearly with the spline degree (fat vertices).If a vertex is the corner of the respective patch, there is only one single non-zero basis function.This means that our choice coincides with the standard choice of corner values.We introduce a scaled Dirichlet preconditioner for the Schur complement formulation of the IETI-DP system.We show that the condition number of the preconditioned system is bounded by where the constant C > 0 is independent of the grid sizes h k , the patch sizes H k , the spline degree p, the smoothness of the splines within the patches Ω (k) , and coefficient jumps between the patches.C depends on the geometry functions, the maximum number of patches that meet on any vertex, the minimal interface length and the quasi-uniformity of the grids within each patch.
The remainder of the paper is structured as follows.We introduce the model problem in Section 2 and its SIPG discretization in Section 3. In Section 4, we propose the IETI-DP solver.Numerical examples are presented in the subsequent Section 5.The paper is concluded with some final remarks in Section 6.The proof of the condition number bound is given in an Appendix.

The model problem
In this section, we introduce the model problem which we consider in this paper.We use the same notation as in [26].To keep the paper self-contained, we introduce the notation in the following.
First, we introduce the computational domain.Ω ⊂ R 2 is a simply connected and bounded open Lipschitz domain, which is composed of K non-overlapping, simply connected open patches Ω (k) , i.e., where T denotes the closure of the set T .We assume that every patch Ω (k) is parameterized by a geometry function that has a continuous extension to the closure of Ω.In IgA, the geometry functions G k are commonly represented in terms of B-splines or NURBS.For the presented analysis, it is not necessary to restrict ourselves to these representations, as long as the Jacobian of G k and its inverse are uniformly bounded.
We use the common notation for the Lebesgue and Sobolev spaces L 2 (Ω) and H s (Ω), s ∈ R, respectively.Those function spaces are equipped with the standard norms denotes the subspace of functions vanishing on ∂Ω.
The boundary value problem of interest reads as follows.Find u ∈ H 1 0 (Ω) such that with a given source function f ∈ H −1 (Ω) and a uniformly positive and bounded diffusion coefficient α, which is constant on each patch, i.e., we have with α k > 0 for all k = 1, . . ., K. Since, for simplicity, we represent the Dirichlet boundary conditions in a strong sense, we assume that the pre-images Γ (k) ) of the Dirichlet boundary Γ D = ∂Ω consist of whole edges of the parameter domain Ω.An alternative, where this restriction would not be necessary, would be a fully floating IETI-DP discretization.

The discretization using IgA and SIPG
In this section, we first introduce the patch-local discretization spaces, then we discuss the overall discretization.For the local discretization spaces, we restrict ourselves for simplicity to B-splines.Let p ∈ N := {1, 2, 3, . . .} be the spline degree, where we assume for simplicity that the degree is uniform for all patches.For each patch k ∈ {1, . . ., K} and each spatial dimension δ ∈ {1, 2}, we introduce a p-open knot vector Ξ (k,δ) = (ξ , where each inner knot might be repeated up to p times.Depending on the p-open knot vector, we introduce a B-spline basis (B[p, Ξ (k,δ) , i]) n (k,δ) i=1 via the Cox-de Boor formula, cf.[6, Eq. (2.1) and We use the standard tensor-product B-spline space V (k) over Ω, which is obtained from the tensor-product space of the two univariate spline spaces.The corresponding physical space V (k) of V (k) is defined by the pull-back principle, i.e., where v| T denotes the restriction of v to T (trace operator).The grid size h k on the parameter domain and h k on the physical domain are defined by Figure 1 shows the basis functions for the univariate case.Here and in what follows, we identify each basis function with its Greville point.We observe that there is only one active basis function on each end point of the interval; its Greville point is located on that end point.
Consequently, the standard tensor-product basis is represented as a grid of Greville points, cf. Figure 2, where an example with five patches is depicted.Note that, in the physical domain, the patches adjoin directly.Since we employ a discontinuous Galerkin discretization, the basis functions at the interfaces do not agree.Therefore, we separate the patches visually.The patches Ω (1) , Ω (3) , Ω (4) and Ω (5) meet in a regular corner.Certainly, it is also possible that only three or more than four patches meet in a regular corner.Such junctions have been previously considered.Here, we additionally allow T-junctions, like the junction between the patches Ω (1) , Ω (2) and Ω (3) .Certainly, it is possible that more than three patches meet in a T-junction.Note that in any case, the T-junction constitutes a corner of every involved patch but one patch, like patch Ω (1) in the example.
Having all patchwise discretization spaces defined, we obtain the global approximation space by On this discontinuous discretization space, we introduce a variational formulation of the model problem (2) following the symmetric interior penalty discontinuous Galerkin (SIPG) method.Find u = (u (1) , Ω ( 1) Ω ( 4) V (1) V (4) V (5)   Figure 2: Schematic representation of local spaces where and δ > 0 is some suitably chosen penalty parameter and n k is the unit normal vector pointing outwards of the patch Ω (k) and Γ (k, ) := ∂Ω (k) ∩ ∂Ω ( ) is the interface between the two patches.N Γ (k) contains the indices of the neighboring patches Ω ( ) , sharing with Ω (k) more than just a corner.
The penalty parameter δ ensures that the bilinear form a h (•, •) is bounded and coercive in the dG-norm A suitable δ can always be chosen, independently of the spline degree p and grid sizes h k , however it might depend on the geometry functions and on the quasi-uniformity of the grids, i.e., the ratios h k / h min,k , see [28,Theorem 8].The Theorem of Lax-Milgram guarantees the existence and the uniqueness of a solution to (5).If the solution u of the continuous problem is sufficiently smooth, the solution of ( 5) is an approximation to the solution of the original problem (2), cf.[28,Theorems 12 and 13].
4 The dG IETI-DP solver In this section, we introduce a IETI-DP solver for the discontinuous Galerkin discretization (5).The first step is the introduction of the local subspaces required for the domain decomposition method and the local assembling of the problem, see Subsection 4.1.In Subsection 4.2, we discuss the introduction of the primal degrees of freedom.Then, in Subsection 4.3, we discuss the coupling of the remaining degrees of freedom using Lagrange multipliers.The setup of the IETI system is discussed in Subsection 4.4, its solution is discussed in Subsection 4.5, and the corresponding convergence result is stated in Subsection 4.6.

Local subspaces and local problem
As it has been done in the seminal paper [19] and in follow-up publications on IETI-DP, the local spaces are constructed on a per-patch basis.For variational problems that are discretized using dG approaches, the setup of local spaces is not obvious.We follow the approach that has been first introduced in [8] and then adapted to IgA in [14,13,26,27]: The introduction of artificial interfaces.
The local function space for a patch Ω (k) is composed of the original function space V (k) of the patch and of the traces of the function spaces V ( ) , when restricted to the common interface Γ (k, ) .Formally, we have The discretization space is visualized in Figure 3. Again, we depict the the interfaces and artificial interfaces separately since there live different function spaces.We again represent every basis function with its Greville point.Basis functions that origin from the same basis are denoted by the same symbol.The function spaces for the artificial interfaces, like V (2,1) , are the traces of the corresponding spaces, here V (1) .Their basis just consists of the traces of those basis functions of the basis of V (1) that are active on the interface Γ (1,2) .While this is rather obvious for the case of interfaces that Ω ( 1) Ω ( 4) V (1,2)  V (1,3)   V (1,4)   V V (2,1) V (3,5)   V (4)   V (4,5)   V (4,1) V (5)   V (5,4)   V (5,3)   Figure 3: Schematic representation of local spaces with artificial interfaces span over a whole edge, it needs some more elaboration in the context of T-junctions.
As one can see in Figure 4, basis functions on the bottom side of the patch Ω (1) are chosen to be part of the artificial interface if they do not vanish on Γ (1,2) .This includes basis functions whose Greville point is not located on Γ (1,2) .The corresponding basis functions form a part of both, the basis for V (2,1) and the basis for V (3,1) .In Figure 3, we depict these degrees of freedom on extensions of the artificial interfaces.
We introduce local bilinear forms a = where we write with a slight abuse of notation gives the local linear system 4.2 Primal degrees of freedom Ω ( 4) Concerning the choice of the primal degrees of freedom, we follow the idea of vertex values.If a vertex x happens to be located on a corner of a patch Ω (k) , on each of the corresponding patches, there is only one basis function active.In this case, the primal constraint enforces for all neighbors Ω ( ) that share the the vertex x.The primal constraint is enforced by requiring that the coefficients for the corresponding basis function agree.
If a vertex happens to be a T-junction x, we select all basis functions which do not vanish on the T-junction.All corresponding degrees of freedom are then treated as primal degrees of freedom.The primal constraint is again enforced by requiring that the coefficients for the corresponding basis functions agree.
The setup of the primal degrees of freedom is visualized in Figure 5.

Jump matrices and Lagrange multipliers
Ω ( 4) Ω (5)   Figure 6: Action of the Lagrange multipliers In the following, we introduce constraints that ensure the continuity of the solution between the interface on one patch and the corresponding artificial interfaces of the neighboring patches, i.e., between the function in V (k) and the functions in V ( ,k) .Note that the function spaces k) and the corresponding bases agree.This means that we obtain continuity if the corresponding coefficients agree.Since the primal degrees of freedom are enforced strongly, there are no constraints corresponding to the primal degrees of freedom, see Figure 6.
The constraints are represented by a matrix where each row corresponds to one constraint enforcing the agreement of one of the coefficients in the usual way, i.e., such that every row has two non-zero entries: +1 and −1.The choice Bv = 0, corresponds to a function v that satisfies the constraints.Note that the proposed choice of primal degrees of freedom guarantees that there is no degree of freedom which is affected by more than one constraint.This is a property, which we use in the condition number analysis.

IETI system
Before we are able to setup the overall IETI system, we partition the degrees of freedom into the primal degrees of freedom (index Π), the degrees of freedom, which are subject to Lagrange multipliers, (index ∆), and the remaining, i.e., interior, degrees of freedom (index I).Using this partitioning, the matrices A (k) , B (k) and the vector f (k) e have the following form: Here, we make use of the fact that there are no Lagrange multipliers that are acting on the interior degrees of freedom (and thus B (k) I = 0) or on the primal degrees of freedom (and thus B (k) Π = 0).On each patch we eliminate the primal degrees of freedom.So, we define , and f These local matrices are collected to global matrices A = diag(A (1) , . . ., A (K) ), A = diag( A (1) , . . ., A (K) ) and B = ( B (1) , . . ., B (K) ), and the local vectors into global vectors ) , . . ., ( f ) ) and f = ((f e (1) ) , . . ., (f e (K) ) ) .
For the setup of the primal problem, we introduce an A-orthogonal basis.To do so, we first introduce for each patch an A (k) -orthogonal basis via Let R (k) be a binary matrix that restricts a global coefficient vector of the primal degrees of freedom to the primal degrees that are associated to space e .Then, we obtain the matrix representing the global A-orthogonal basis for the primal degrees of freedom via R (1)  . . .
The overall IETI-DP system reads as follows.Find ( u , u Π , λ ) such that This problem is equivalent to the original problem (5), cf.[21].
Remark 4.1.In this paper, we follow the approach to eliminate the primal degrees of freedom, which is a commonly used approach for handling the corner values in actual implementations.Alternatively, one can incorporate the primal constraints using Lagrange multipliers, which is a common approach if edge averages are used as primal degrees of freedom.Certainly, this approach is also possible in the framework of this paper.Here, we would obtain the formulation (8), however with the choice The matrix (0 0 I) in the definition of A (k) here is often called C (k) .
In theory papers, cf.[21], a FETI-DP or IETI-DP system is often written down in the equivalent skeleton formulation, which corresponds to the choice

Solving the IETI system
By applying a block-Gaussian elimination to (8), we obtain the Schur complement equation for the Lagrange multipliers λ, where We solve (9) with a preconditioned conjugate gradient (PCG) solver.Let us define Γ , . . ., B Γ ).The preconditioner for the PCG method is the scaled Dirichlet preconditioner M sD defined by where S = diag(S (1) , . . ., S (K) ) with is the restriction of the overall operator A to the skeleton and D = diag(D (1) , . . ., D (K) ) is a diagonal matrix defined based on the principle of coefficient scaling: Each coefficient d for degree of freedom i associated to the interface Γ (k, ) .If i corresponds to a primal degree of freedom, then can be chosen arbitrarily among the indices of the neighboring patches.
After solving the system (9), the solution vectors u and u Π and, finally, u are computed from u by means of simple patch-local postprocessing steps.
The execution of the IETI-DP method for the dG discretization described above requires basically the same computational steps as the IETI-DP method for dG discretizations on conforming patch decompositions, see [26,Section 3] for the detailed outline of the algorithm.

Condition number estimate
The following theorem allows to estimate the maximum number of iterations that the PCG solver with scaled Dirichlet preconditioner needs to reach a desired error tolerance.The condition number (and thus the number of iterations) depends on the patch sizes, the grid size and on the spline degrees as explicitly stated in the theorem (the constant C does not depend on these quantities).The dependence on the grid sizes and patch sizes is as expected for FETI-like methods.Moreover, the dependence on gird and patch sizes and spline degree is the same as for the continuous case in IgA, see [25].The condition number bound is independent of the diffusion parameters α k , of the number of patches K, of the continuity of the spline spaces, and of the choice of the penalty parameter δ (provided δ is large enough such that the overall bilinear form is coercive).The constant C (and thus the condition number) also depends on the bounds for the geometry function, on the number of patches that meet in a vertex and the quasi uniformity of the grids.
Theorem 4.2.Provided that the IETI-DP solver is set up as outlined in the previous sections, • there is a constant C 1 > 0 such that • there is a constant holds for all vertices x, • there is a constant C 3 > 0 such that holds for all k = 1, . . ., K and all ∈ N Γ (k), • and the grids are quasi-uniform, i.e., there is a constant C 4 > 0 such that holds for all k = 1, . . ., K, then the condition number of the preconditioned system satisfies the constant C only depends on the constants C 1 , C 2 , C 3 and C 4 .5 Numerical results In this section, we apply the IETI-DP method to a simple magnetostatic problem.We consider the computational domain shown in Figure 7 consisting of 272 patches representing a simplified cross section of an interior permanent magnet electric motor (IPMEM).The different colors in Fig. 7 denote different materials.The redbrown patches denote ferromagnetic material, e.g., iron, the yellow patches are the permanent magnets and the blue patches represent air regions and coils made of copper (for which we use the same material parameters).
The considered boundary reads formally as follows.Find u such that −div(ν(x, y)∇u(x, y)) = div(ν(x, y)M (x, y)) for (x, y) ∈ Ω u = 0 on ∂Ω, where ν denotes the magnetic reluctivity and M denotes the magnetization.The magnetic reluctivity is ν ferro = 1 204π 10 5 on the ferromagnetic parts, ν mag = 1 4.344π 10 7 on the permanent magnets and ν air = 1 4π 10 7 on the air and copper regions.This means that we have a jump in the order of approximately 10 4 .On each of the permanent magnets, the magnetization M is given by where ρ mag := 1.28 is the magnetic remanence, and n is unit the normal vector in positive or negative radial direction (measured from the center of the magnet), where a positive sign is used for every second magnet and a negative sign for every other second magnet.The magnetization M vanishes on the remainder of the domain (ferromagnetic parts, air, copper).
The cross section of the motor is modeled with NURBS and B-splines.For the coarsest discretization space, i.e., r = 0, we use B-splines that are global polynomials and we use only splines of maximum smoothness within the patches.The subsequent refinements  9) with the M sD preconditioner that arises from the magnetostatic model problem with a PCG solver and start the iterations with zero initial vector.We stop the iteration if the 2 -norm of the residual has been decreased by a factor of 10 −6 compared to the 2 -norm of the right-hand side.We use the penalty parameter δ = 12 for all the numerical experiments which are carried out on the Radon11 cluster located in Linz and we used the C++ library G+Smo [22].Table 1 shows the condition numbers and the iteration counts for the magnetostatic model problem.We see in this table a condition number growth with respect to h as expected and we also see that the condition numbers decrease for higher polynomial degrees from refinement level 4 on.3: Iteration counts (it) and condition numbers (κ); rotation dependence reluctivity 10 j for r = 5 refinement steps.We see in the table that the condition number is almost independent of the coefficient jumps.
In Table 3, we see the dependence of the iteration and condition numbers with respect to the angle of rotation of the motor.We choose three different angle positions ϕ at 5  36 π, 6 36 π and 7 36 π.We observe from this table that the iterations and condition numbers decrease for a larger angle ϕ.
The Table 4 shows the solving times in seconds (sec.)required to solve the IETI-DP system with the number of computing cores given in the column proc.We increase the number of cores from 2 to 16.We see in this table a very good scaling behavior of the algorithm.

Conclusions
In this paper, we have constructed a IETI-DP algorithm for computational domains with a non-matching decomposition into patches.We have adapted the idea of using corner values as primal degrees of freedom (Alg.A) from [26] according to our requirements.In this paper, we have generalized this idea to T-junctions: We add basis functions that are supported on a vertex to the primal space.For this choice, we obtain the same h and p-explicit condition number bounds as in [26].

A Appendix
In the appendix, we give a proof of Theorem 4.2.Throughout this appendix, we use the notation a b if there is a constant c > 0 that only depends on the constants C 1 , C 2 , C 3 and C 4 from Theorem 4.2 such that a ≤ cb.Moreover, we write a b if a b a.
When it is clear from the context, we do not denote the restriction of a function to an interface explicitly, so we write for example u ) .Following the usual approach, for the analysis, we need to introduce the skeleton representation of the solution which is obtained by eliminating the interior degrees of freedom.By eliminating the interior degrees of freedom from the spaces V (k) , we obtain the space W (k) := {v| ∂Ω (k) : v ∈ V (k) }.The introduction of the skeleton representation has no influence on the function spaces on the artificial interfaces, thus we define W (k, ) := V (k, ) .Based on these choices, we define analogously to V and V (k) e the function spaces Analogously to (6), a function w has the form w k) , where w (k) ∈ W (k) and w (k, ) ∈ W (k, ) .A basis for W (k) e is canonically defined by choosing the traces of the basis functions of the basis of V (k) (for which the trace does not vanish) and the basis functions of the bases of V (k, ) = W (k, ) .Finally W ⊆ W is the subspace of functions where the primal constraints are satisfied, i.e., the coefficients for the vertex basis functions agree.
As in [25], we define the seminorm for a continuous function v over the set T ⊂ R 2 .Moreover, we use the standard seminorm for T being the boundary or an edge of a patch.
The following lemma allows to estimate the action of the matrix B D B Γ , where we define e , • • • , u e ) = ((u (1) , (u (1, ) ) ∈N Γ (1) ), . ..) ∈ W with coefficient vector u and let w = (w e ) = ((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 interface Γ (k, ) that Proof.As in the proof of [26,Lemma 4.3], we have where ϕ (k) i denotes a basis function of the basis of W (k) and ϕ (k, ) i denotes a basis function of the basis of W (k, ) .The set B T (k, ) contains the pairs of indices of those basis functions in the bases for W (k) and W ( ,k) which are subject to a primal constraint.
We obtain the representation (13) since the coefficients corresponding to all basis functions on the common edge Γ (k, ) are equal to ±α /(α k + α ), except to the basis functions that correspond to the primal degrees of freedom.For the latter, the corresponding coefficients are 0 since the primal degrees of freedoms are not subject to the jump matrix.Thus, we subtract the latter.
Note that u ∈ W , which means that it satisfies the primal constraints.Hence the sum over the indices in B T (k, ) vanishes.Therefore, we immediately obtain the desired result.
Lemma A.1 and the triangle inequality immediately yield Here and in what follows, we write h k := min{ h k , h } and h k := min{h k , h }.
We estimate contributions from the artificial interfaces in the H 1/2 -and L 0 ∞ -seminorms.We start with the H 1/2 estimate.Lemma A.2. Let u ∈ W .Then, the estimate holds.
Proof.Using the triangle inequality, we obtain We apply [26, Lemma 5] and the norm equivalence, [26, Lemma 1], to the difference We use (17) to obtain An application of the norm equivalence, [26, Lemma 1], yields the estimate . Lemma A.2 and (18) finish the proof of this lemma.
Before we give a proof of the main theorem, we estimate the sum of the corresponding seminorms over all patches.

S
The combination of this estimate and (19) finishes the proof.
, •) and the local linear functional f (k) e , • that live on the spaces V (k) e .They can be seen as the local counterparts to a h (•, •), d(•, •), and f, • and we define them by

Figure 4 :
Figure 4: Basis functions selected for artificial interfaces; the dashed basis functions are selected for both bases

Figure 7 :
Figure 7: Decomposition of cross section into patches (left), and their materials (right)

Figure 8 :
Figure 8: Solution of the linear model problem

Table 2
reports on the robustness of the IETI-DP solver with respect to coefficient jumps.For the numerical tests we assign to the ferromagnetic patches the hypothetical