A projected super-penalty method for the $C^1$-coupling of multi-patch isogeometric Kirchhoff plates

This work focuses on the development of a super-penalty strategy based on the $L^2$-projection of suitable coupling terms to achieve $C^1$-continuity between non-conforming multi-patch isogeometric Kirchhoff plates. In particular, the choice of penalty parameters is driven by the underlying perturbed saddle point problem from which the Lagrange multipliers are eliminated and is performed to guarantee the optimal accuracy of the method. Moreover, by construction, the method does not suffer from locking also on very coarse meshes. We demonstrate the applicability of the proposed coupling algorithm to Kirchhoff plates by studying several benchmark examples discretized by non-conforming meshes. In all cases, we recover the optimal rates of convergence achievable by B-splines where we achieve a substantial gain in accuracy per degree-of-freedom compared to other choices of the penalty parameters.


Introduction
Isogeometric analysis (IGA), firstly introduced in [23], is a methodology used for the numerical discretization of partial differential equations (PDEs) based on the same building blocks used in Computer Aided Design (CAD). Indeed, in IGA, the same mathematical objects, such as B-splines and non-uniform rational B-splines (NURBS) [35], used for the geometrical description are employed for the numerical solution of the PDE at hand. A distinguishing feature of splines is the high regularity achievable by construction, which allows the approximation of higher-order variational problems directly in their primal, for instance Kirchhoff plates [34,30], Kirchhoff-Love shells [24,37,26,27] and the Cahn-Hilliard equation [16]. For a detailed review of the method and its recent applications, the reader is referred to [23,11,1], whereas its mathematical foundations can be found in [5,12].
Although smoothness is attained naturally within a patch, geometries of engineering relevance are in general described by multiple patches, where typically the underlying spline representations are non-conforming at the common interface. Clearly, in this scenario, a direct strong coupling between patches is not straightforward to achieve. Moreover, as in the scope of this work we are interested in the Kirchhoff plate model problem, an efficient strategy to obtain C 1 -coupling is needed since a global C 1 -continuity is required to obtain a well-defined bilinear form for the problem at hand. In the literature, three methods are predominantly used to achieve the latter coupling in a weak sense and they are summarized in the following. High-order mortar methods have been studied in [22,20] in the context of Kirchhoff plates and Kirchhoff-Love shells, respectively, and have been extended to a general C n -coupling in [13]. For a detailed review in the context of isogeometric analysis, we refer to the review article [19]. However, mortar methods leads to the formulation of a saddle point problem, where the associated Lagrange multipliers constitute additional unknowns to be solved for in the global system of equations. Nitsche method has been analyzed in [39] for coupling isogeometric Kirchhoff plates in the scope of immersed methods and in [17] for imposing weakly kinematic boundary conditions for fourth-order PDEs. Although this family of method is less sensitive to the choice of parameters compared to classical penalty approaches, their formulation requires additional consistency terms which, in the Kirchhoff problem, involve the computation of derivatives of shape functions up to order three. This adds some extra steps of complexity in the implementation and increases the overall computational cost of the coupling strategy. Finally, penalty methods are widely used in the engineering community due to their conceptual simplicity, see the seminal work [3]. Furthermore, they can be easily and efficiently incorporated into a numerical code, where we refer to [25,2,15,18] for more insights and some applications in the context of isogeometric Kirchhoff-Love shells. Nonetheless, a major drawback of this approach resides in their lack of robustness with respect to the choice of penalty parameters. Typically, the choice of penalty coefficients is problem-dependent and is based on a time-consuming, heuristic process. As noted in [18], on one hand, if the penalty factors are chosen too small the interface constraint is satisfied only loosely. On the other hand, if the coefficients are too high, the condition number of the resulting system matrix is negatively impacted and the convergence behavior is spoiled by spurious locking phenomena.
Our contribution falls into this realm. Inspired by the super-penalty method studied in [4], our goal is to introduce a simple coupling procedure for the displacement and rotation fields, respectively, for non-conforming multi-patch Kirchhoff plates, which preserves the high-order optimal convergence rates achievable by B-splines while mitigating the detrimental effects related to locking. To alleviate the over-constraint of the solution space we perform an L 2 -projection of the penalty terms onto a space of reduced degree defined on the slave side of the coupling interface where, motivated by the work in [8] for mortar methods, we select a p/p − 2 pairing, where p denotes the B-splines degree. In particular, starting from the perturbed saddle point formulation of the Kirchhoff plate model problem, we show how the corresponding Lagrange multipliers can be eliminated from the system and, more importantly, how the perturbation gives us insights into the optimal choice for the penalty coefficients. Indeed, the proposed methodology is truly parameter-free, as the penalty factors are fully determined by the given physical constants, the geometry and its discretization, i.e. mesh size and spline degree. We remark that the proposed methodology is especially advantageous for moderate degrees p = 2, 3, where locking phenomena are particularly pronounced and the L 2 -projection proves to be an effective and computationally efficient remedy. Then, we address the ill-conditioning issues stemming from our choice of super-penalty parameters. We adapt the block preconditioner based on an inexact Schur Complement Reduction (SCR) introduced in [28,29] and we combine it with a preconditioner tailored to the isogeometric discretization of the Kirchhoff plate, where we exploit the tensor product structure of B-splines and an efficient algorithm for the solution of the arising Sylvester-like system; for a detailed derivation we refer to [38,32,31]. Finally, we show through several numerical benchmarks the optimal convergence properties of the presented methodology, where our approach does not suffer from locking also on very coarse meshes. This leads to a substantial improvement in the accuracy achievable per degree-of-freedom (dof).
The structure of the paper is as follows. Section 2 provides a review of the fundamental concepts related to B-splines. Section 3 describes in details the derivation of the proposed methodology and motivates our choice of penalty parameters. Section 4 presents the ideas used in the construction of the preconditioner employed in this work. In Section 5 the method is validated on several numerical benchmarks and it is applied to the analysis of an idealized multi-patch design of an L-bracket. Finally, some conclusions are drawn in Section 6.

A brief introduction to B-splines
In this section, some definitions and fundamentals related to B-splines and NURBS are reviewed. We refer the reader to [35,11,21], and references therein, for a comprehensive review of B-splines and their role in isogeometric analysis. Starting from two integers p, n, a univariate B-spline basis function b i,p of degree p is generated starting from a non-decreasing sequence of real values referred to as knot vector, denoted in the following as Ξ = {ξ 1 , . . . , ξ n+p+1 }. It is worth mentioning that the smoothness of the obtained Bspline basis is C p−k at every knot, where k denotes the multiplicity of the considered knot, while it is C ∞ elsewhere. In the remainder of this work, we consider only splines of maximum continuity, i.e. C p−1 . The definition of multivariate B-splines B i,p (η) is achieved in a straight-forward manner using the tensor product of univariate B-splines as: where d denotes the dimension of the parameter space. Additionally, the multi-index i = i 1 , ..., i d denotes the position in the tensor product structure and p = p 1 , ..., p d indicates the vector of polynomial degrees, associated to the corresponding parametric dimension η = η 1 , . . . , η d , respectively. Then, let us define a domain Ω ∈ R d described by a B-spline parametrization F as a linear combination of multivariate B-spline basis functions and corresponding control points as follows: where the coefficients P i ∈ R d of the linear combination are the control points and d represents the dimensionality of the physical space. Although not treated here, it is straightforward to extend the notation to NURBS, for details see [11]. In the rest of the paper, without loss of generality, the degree vector p will be considered equal in each parametric direction and therefore simplified to a single scalar value p. Further, the vectors i and η will be omitted to simplify the notation. Finally, we can introduce the following discrete space formed by multivariate B-splines of degree p:

The strong form of the Kirchhoff plate problem
Let us introduce the governing PDE, characterized by the bilaplace differential operator, that describes the bending-dominated problem of a Kirchhoff plate, following the notation in [37]. Let us define an open set Ω ⊂ R 2 with a sufficiently smooth boundary ∂Ω, such that the normal vector n to the boundary is well-defined (almost) everywhere. Let us also introduce two admissible splittings of the boundary Γ = ∂Ω into Γ = Γ u ∪ Γ Q and Γ = Γ φ ∪ Γ M , such that Γ u ∩ Γ Q = ∅ and Γ φ ∩ Γ M = ∅, respectively. Consequently, the strong form of the problem reads: where u represents the deflection of the plate, D its bending stiffness, ν is the Poisson ratio, g is the load per unit area in the thickness direction, u Γ , φ Γ , M Γ and Q Γ are the prescribed deflection, rotation, bending moments and effective shear, respectively. The bending stiffness D is defined as: where E is the Young modulus and t denotes the thickness of the plate. For the sake of simplicity and without loss of generality, these are assumed to be a constant in Ω. Finally, the differential operator Ψ (·) reads:

The multi-patch formulation of the perturbed saddle point Kirchhoff problem
Here, following the notation used in [8], we introduce a decomposition of Ω into N non-overlapping subdomains Ω i such that: Now, let us define the interface γ k, between two adjacent patches Ω k , Ω , 1 ≤ k, ≤ N as the intersection of their corresponding boundaries: Then, the skeleton Γ is defined as the union of all non-empty interfaces (which we suppose to be labeled with an index = 1, . . . , L) and reads: Consequently, we can denote by u k and n k the value of the primary field and the outward normal on ∂Ω k , and u and n the value of the primary field and outward normal on the neighboring subdomain ∂Ω , see Figure 1 for an example on two patches. Then, for each interface γ k, we can write the following coupling conditions: which can be rewritten using the standard jump and normal jump operators, respectively, as: Further, given 1 ≤ s, t ≤ L, s = t, we denote the cross-points by c s,t = γ s ∩ γ t and we label them with an ordered index c s , s = 1, . . . , S. For ease of notation and without loss of generality, in the following we assume the flexural rigidity D to be constant in Ω and the Poisson ratio ν to be zero. Further, we assume that the values prescribed as natural boundary conditions are zero as well. Now, let us introduce for each subdomain Ω i the following space: from which the following broken Sobolev space can be characterized as: Then, let us also define the spaces: Lastly, we need to introduce the following dual spaces: We are now ready to formulate (1) as a saddle point problem.
We also define three continuous bilinear forms a : Under suitable regularity assumptions, we can provide an estimation of the error introduced by the perturbations ε where we have defined:ε

The projected super-penalty formulation
For each patch Ω i , we assume p ≥ 2 and we indicate with S p h (Ω i ) the space trivially obtained extending by zero the elements of S p h (Ω i ) over Ω \ Ω i . Additionally, let us define: Consequently, let us denote by V i,h ⊂ X i,h the finite-dimensional space given by the span of B-splines defined on the corresponding subdomain Ω i , where the exact characterization of V i,h depends on the chosen boundary conditions, for further details we refer to [10]. This allows us to introduce the following finite dimensional subspace of V , Moreover, for each interface γ , we denote by Ξ the knot vector on γ inherited from the slave side. Motivated by the choice of the p/p − 2 stable pairing in [8], we construct the following isogeometric space S p−2 h (γ ) on the reduced knot vector Ξ obtained by removing from Ξ the first and last two knots, where an example is depicted in Figure 3 for p = 2, 3. Similarly to before, we indicate with S p−2 h (γ ) the space obtained extending by zero over Γ \ γ the elements of S p−2 h (γ ). We can now define the discrete counterpart of the Lagrange multiplier spaces as: With these definitions at hand, the discretized version of (4) reads: where α defl and α rot are "large" parameters associated to the deflections and rotations, respectively. In general, they depend on the problem definition, e.g. the physical constant D, the mesh size and spline degree, where a full characterization of our choice will be given later in the section. We can now formally eliminate the Lagrange multipliers and recast (6) into its primal form. Indeed, we can write: . Finally, employing the previous results and the properties of the L 2 -projection, the resulting discretized bilinear form, augmented by suitable penalty terms that weakly enforce the coupling conditions (2), reads: find

Inf-sup test
The well-posedness of (6), independently of the value of the parameters α defl and α rot , relies on the well-posedness of the underlying unperturbed problem, i.e. the problem corresponding to (6) where we set α defl = α rot = +∞ , = 1, . . . , L. Although a rigorous proof of the inf-sup stability of such unperturbed problem is currently under investigation [9], we assess the behavior of the numerical infsup test for a domain Ω subdivided along a straight interface into two subdomains Ω i , i = 1, 2. As we are dealing with a double saddle point problem, we compute two different inf-sup constants C inf-sup defl and C inf-sup rot , corresponding to the deflection and rotation jumps, respectively. In the following we report the results for different discretization sizes of the interface h = 1/2 k , k = 3, . . . , 7 and Bspline degrees p = 2, . . . , 5, where h denotes the maximum mesh size associated to the interface γ . The numerical values of C inf-sup defl and C inf-sup rot are summarized in Table 1 and are depicted for clarity in Figure 2. In all cases we observe that the inf-sup constants converge to some values bounded away from zero, numerically suggesting that the method is inf-sup stable.

Coercivity test
Then, we also assess numerically the behavior of the coercivity constant on an example with four patches Ω i , i = 1, . . . , 4 separated by four straight interfaces meeting at a cross-point. In particular, we want to compute the biggest α 0 such that: where B 1 and B 2 are the linear operators associated to the bilinear forms b 1 and b 2 , respectively. The results for different discretization sizes of the interface h = 1/2 k , k = 2, . . . , 5 and B-spline degrees p = 2, . . . , 4 are presented in Table 2, from which we can numerically infer that the method is coercive on the intersection kernel.

Remark 1. The inf-sup and coercivity tests are performed on a reduced version of the knot vector Ξ ,
where also the first and last internal knots of Ξ are eliminated. This is justified by our preliminary mathematical analysis, where this choice is required. However, from a numerical standpoint, we  , respectively, on different uniformly refined meshes h = 1/2 k , k = 3, . . . , 7 and spline degrees p = 2, . . . , 5.  Table 2: Numerical estimation of the coercivity constant α 0 on different uniformly refined meshes h = 1/2 k , k = 2, . . . , 5 and spline degrees p = 2, 3, 4.
retain the optimality of the method without performing such a reduction and in all our examples we directly employ Ξ to define the projection spaces.

On the choice of penalty parameters
It is well-known that the penalized problem (7) is variationally consistent only in the limit α defl = α rot → ∞ = 1, . . . , L. On the other hand, the well-posedness of this problem is robust with respect to the choice of the parameters α defl and α rot . Therefore, the proposed methodology will not suffer from locking for any choice of penalty values. As a consequence, α defl and α rot can be chosen solely to guarantee the optimal accuracy of the method.

Remark 2.
A clear trade-off of this choice is the negative impact on the conditioning of the resulting system matrix. A possible remedy based on an ad-hoc preconditioner will be discussed in a later section. Another drawback consists in the loss of significant digits due to the (potentially big) difference in magnitude between the penalty contribution and the internal stiffness. For this reason (amongst other which will be pointed out in the rest of manuscript), we advise to use this method in combination with splines of degree p = 2, 3, as these round-off errors occur below a tolerance threshold of significance to most engineering applications.
Inspired by the method proposed in [18] in the context of Kirchhoff-Love shells, we want to develop a fully parameter-free penalty method. To this end, we scale the deflection and rotation penalty parameters by the physical constants, the local mesh size and the geometry as: where the exponent β is chosen to ensure the optimal convergence of the method with respect to the degree p of the underlying discretization. Note that all of these parameters are known and depend only on the problem definition, meaning that no user-defined factor is required. We highlight that our choice is based on the fact that the perturbations introduced in (4) cannot be "big" compared to the accuracy with which we want to solve the original problem and the estimate provided in (5) guides the choice of β. Moreover, as we want to recover optimal rates of convergence for the error, the exponent β must be a function of the underlying splines degree p. From the numerical experiments conducted thus far, the scaling factor β = p − 1 in (8) is necessary to ensure optimal convergence of the method in the H 2 norm, whereas for a scaling of β = p we observed optimality in the H 2 and H 1 norms. Finally, a factor of β = p+1 provides optimality in the H 2 , H 1 and L 2 norms. If not stated otherwise, we will use β = p + 1 in all our numerical examples.

Remark 3.
Although a rigorous mathematical proof of the method and the optimal choice of β are currently under development [9], we believe that this allows for some extra flexibility in the proposed methodology, where the suitable scaling factor can be chosen with respect to the corresponding quantity of interest.

Cross-points modification
In the literature of mortar methods, it is well-known that the treatment of cross-points requires extra considerations, see [14] and references therein for a discussion in the context of mortar coupling of isogeometric multi-patches. Analogously, our method also inherits the need for a cross-points modification. Indeed, in order to retain optimality of the method, a linear constraint must be (a) Interface splines of degree p red = p − 2 = 0 associated to a reduced knot vector on the slave side Ξ = [0 1/3 2/3 1] and corresponding intersection mesh for integration.
x y Ω 1 Interface splines of degree p red = p − 2 = 1 associated to a reduced knot vector on the slave side Ξ = [0 0 1/3 2/3 1 1] and corresponding intersection mesh for integration. Figure 3: Example of the projection setup on a coupling interface. We select the finer side (on Ω 1 in this example) to define the reduced space for the projection. Additionally, an intersection mesh at the interface is created only for integration purposes to properly compute the projected penalty terms.  imposed to the control variables meeting at the cross-point to ensure C 0 -continuity. An example with four patches is depicted in Figure 4, where in Figure 4a we depict the dofs associated to each coupling interface and in Figure 4b we visualize the imposition of the constraint. To explain the procedure, let us start from the following unconstrained system of equations: Now, the constraint can be incorporated easily into the standard linear system in a fully algebraic fashion, where a possible implementation is presented in Algorithm 1.

Algorithm 1
Algorithm for applying a C 0 constraint at a cross-point. 1: procedure Apply C 0 constraint(vector of dofs at cross-points u cp )

2:
Label one dof in u cp as master 3: Label the remaining dofs in u cp as slaves 4: Build the rectangular matrix C representing the linear master-slaves constraints (see (10)) 5:

Solve the reduced system
Recover the solution u h from u h = C u h 7: end procedure The construction of the rectangular matrix C is best explained with an example. Let us assume that the dofs at the cross-point are numbered as u cp = [u cp1 u cp2 u cp3 u cp4 ]. Now, without loss of generality, we pick u cp1 as the master control point and the rest as slave nodes. Then, the constraint can be expressed via the matrix C as follows: where ndof denotes the total number of degrees-of-freedom in the system. This procedure allows to eliminate the unknowns associated to the slave nodes from the system.

A nested preconditioner based on the Schur Complement Reduction
In this section, following the notation introduced in [36] and building upon the work presented in [28,29] in the context of elastodynamics and hemodynamics, we present an efficient way to mitigate the detrimental effects on the condition number stemming from our choice of super-penalty parameters. This preconditioner is based on the approximate solution of the block factorization of the system matrix known as Schur Complement Reduction (SCR). We remind the reader that before performing the algorithm described in the following, we apply a symmetric diagonal scaling to the system matrix.

The Schur Complement Reduction
We begin by reordering the matrix A ∈ R ndof×ndof stemming from (7) in blocks as follows: where the subscripts i and Γ refer to internal and interface dofs, respectively, where an example is depicted in Figure 5. Let us remark that A i,i is a block-diagonal matrix where every block is the matrix associated to an homogeneous Dirichlet problem (fully clamped) on the corresponding patch Ω i . Moreover, with a slight abuse of notation, we assume that, if needed, A has already been modified to account for the constraints related to the cross-points introduced in the previous section. Now, we can perform the following block factorization of A: where we have introduced the Schur complement S Γ,Γ : Multiplying with L on both sides we get: We highlight that, up to this point, this factorization is performed in exact algebra. Then, from (11), we can solve for x in a segregated fashion by exploiting Algorithm 2.
Algorithm 2 SCR algorithm 1: procedure Solution of Ax = r based on SCR

2:
Solve for an intermediate solutionx i 3: Update the interface residual r Γ = r Γ − B i,Γx i

4:
Solve for the interface solution x Γ from the Schur equation

5:
Update the internal residual r i = r i − B i,Γ x Γ

6:
Solve for the internal solution x i from 7: end procedure Clearly, the Schur complement S Γ,Γ is in practice expensive and often infeasible to compute explicitly. A way around this issue is given in Algorithm 3, where we summarize a matrix-free procedure to apply the Schur complement to a vector.

Algorithm 3 Algorithm for applying the Schur complement to a vector
1: procedure Application of S Γ,Γ to a vector x Γ

6:
Returnx Γ − x Γ 7: end procedure Remark 4. As noted in [28], the cost of the preconditioner is often dominated by the solution of the Schur system (13). To reduce the computational burden of this step, we use as preconditioner a coarse approximation of the Schur complement obtained by applying only a few iterations of GMRES to A i,i for assemblingS Γ, h and p) and scalable preconditioner for the Schur complement and, more in general, for fourth-order PDEs.

Nested block preconditioner strategy based on SCR
The main idea presented in [28] is to combine the robustness of the SCR factorization with the ease of application of a block preconditioners (such as SIMPLE or variants thereof [36]). Indeed, we can build a preconditioner P SCR based on an approximate factorization of (11), where Equations (12) to (14) are solved within a prescribed tolerance. Given that P SCR changes its algebraic definition at every iteration, following [28], we employ a flexible GMRES algorithm (FGMRES) as the iterative method for the most outer solve Ax = r. At each iteration of FGMRES, we can apply the preconditioner P SCR via Algorithm 2, where this entails the solution of the blocks A i,i and S Γ,Γ . This part of the algorithm is denoted as intermediate solver. Last, since we do not assemble the Schur complement explicitly, but we apply its action on a vector through Algorithm 3, we perform a final solve for A i,i in (15), denoted as inner solver. The final performance of the preconditioner is therefore determined by the prescribed tolerances for the outer, intermediate and inner layers, respectively, where the objective is finding a good balance between the computational cost and the robustness of the method. In the following, we denote the aforementioned tolerances by η o , η t and η n for the outer, intermediate and inner layers, respectively.

A preconditioner based of the Fast Diagonalization (FD) algorithm
Since each outer iteration of the nested preconditioner is based on the solution of three systems involving the block A i,i , an efficient and robust preconditioner for this block is required. In this work, we extend the isogeometric preconditioner studied in [38,32], based on the Fast Diagonalization algorithm, to the Kirchhoff plate problem. In the following, we focus our derivation on the singlepatch case. The extension to the multi-patch case is straightforward by construction, since the block A i,i is formed by disjoint sub-blocks associated to each patch Ω i . Now, exploiting the tensor product structure of the B-spline basis at the patch level, let us introduce the preconditioner P FD in Kronecker form as: where M k and K k with k = 1, 2 refer to the one-dimensional, parametric mass and hessian matrices associated to the k-th parametric dimension, respectively. They can be expanded as follows: where b indicates the univariate B-spline basis functions introduced in Section 2. Then, analogously to [33], we partially include the geometry and physical coefficients inside the preconditioner. In particular, let us denote by C the following function: where we recall that J F represents the jacobian of the B-spline parametrization F and D is the flexural stiffness of the plate. Now, as explained in [33, Appendix A.3], we perform a separation of variables on C such that we can write: where this matrix is evaluated at each quadrature point. With this, we can modify the preconditioner given in (16) to partially account for the geometry and coefficients information as follows: where Finally, each iteration of the iterative solver requires the solution of the following system: where r denotes the current residual. Due to the tensor structure of the preconditioner, we can rewrite (18) as a Sylvester matrix equation [40]: where s = vec(S) and r = vec(R).

Remark 5.
Let us recall that for any matrix Z ∈ R r×c the operator vec(Z) gives as output the vector z ∈ R rc formed by stacking the columns of Z.
Let us now consider the generalized eigendecomposition of the matrix pencils ( K 1 , M 1 ) and ( K 2 , M 2 ), respectively, as: Here, D 1 and D 2 are diagonal matrices containing the eigenvalues of M −1 1 K 1 and M −1 2 K 2 , respectively. Further, U 1 and U 2 are defined as: With these definitions at hand, we can rewrite (17) in Kronecker form as: where the preconditioner can be efficiently applied via Algorithm 4.

5:
Return s = (U 1 ⊗ U 2 )s 6: end procedure Remark 6. We remark that the application of the nested preconditioner P SCR combined with P F FD can be implemented in a fully matrix-free framework. Furthermore, although not investigated in this work, the patch-wise block structure of A i,i could be further exploited for parallelization.
For the sake of conciseness, we do not provide here further details of the FD algorithm, but we refer to [38,32] for a thorough theoretical and numerical investigation of the method in the scope of isogeometric analysis.

Numerical Examples
In this section we assess the performance of the proposed coupling method with several numerical examples defined on multi-patch geometries. All the numerical experiments presented in the following have been implemented in the open-source and free Octave/Matlab package GeoPDEs [41], a software designed for the solution of partial differential equations in the context of isogeometric analysis.

A four patches example with non-matching curved interfaces
In this example we consider the computational domain Ω = [0, 2] × [0, 2] depicted in Figure 6, split into four subdomains Ω i . We remark that all meshes are non-conforming at every coupling interface, as the irrational factor √ 2/100 has been used to shift the interface knots. The body source and boundary data are computed such that the exact solution is smooth and it reads: This setup is used to test the robustness of our method in the case of severe non-matching discretizations and with respect to the problem parameters. To this end, we present the convergence results for all combinations of Young's moduli E = [10 4 , . We compare our method to a classical penalty approach, where we set α defl = α rot = 10 4 E, = 1, . . . , L, and to a choice of penalty parameters scaled with respect to the physical parameters as proposed in [18]. In particular, they read: where the user-defined parameter δ = 10 3 is chosen. From the results in Figure 7, we observe that the projection strategy shows robustness with respect to the input parameters and allows for an easy treatment of locking phenomena, where optimal convergence rates are attained also for very coarse meshes.
In Figure 8 the convergence behavior of the error measured in the H 2 norm with and without the imposition of the C 0 constraint at the cross-point is plotted. We observe that the loss of accuracy hinders the convergence for p = 3, 4, whereas the expected optimal rates of convergence are recovered in all cases when the linear constraint is imposed to the system. This is further highlighted in Figure 9, where the element-wise H 2 error is depicted for a discretization of degree p = 4, without and with the constraint, respectively. On one hand, we remark how the error is concentrated and much higher in the elements around the cross-point, spoiling the optimal convergence, when the constraint is not imposed. On the other hand, with the linear constraint, we recover optimal convergence properties of the method.
Finally, for this example we also analyze the performance of the nested preconditioner. In Table 3 we report the iterations needed by the external solver and in brackets the average number of intermediate iterations, for several degrees of the discretization p = 2, 3, and we compare it with a classical diagonally preconditioned conjugate gradient (PCG), a PCG where an incomplete LU (ILU) is used as preconditioner and a GMRES preconditioned with ILU. All the results refer to a global tolerance (a) Geometry setup and physical parameters.
x y (b) Initial discretization. Figure 6: Problem setup and initial multi-patch non-conforming discretization for the curved four patches example.
η o of 10 −10 and, for the nested SCR-FGMRES strategy, the intermediate and inner tolerances η t and η n are set to 10 −6 . Further, the Schur complement is preconditioned by an approximationS Γ,Γ obtained with a maximum of 6 iterations of GMRES. For the sake of completeness, we perform the same test with the choice of penalty parameters studied in [18]. The results are summarized in Table 4, where we observe no substantial difference regarding the iterations needed to solve the system compared to the case where our choice of parameters is employed. This suggests that the proposed preconditioner is robust with respect to the penalty factors and it is also suitable to precondition systems stemming from other penalty approaches.   (12) to (14) in Algorithm 2, respectively. Iterations marked with -did not reached convergence within the prescribed 1000 maximum number of iterations.
In Table 5 we study the influence of the intermediate and inner tolerances on the number of outer iterations required by the FGMRES solver, on a fixed mesh of 4096 elements, for B-splines of              (12) to (14) in Algorithm 2, respectively. Iterations marked with -did not reached convergence within the prescribed 1000 maximum number of iterations.
degree p = 2, 3. We note that as the chosen tolerances become smaller and smaller, we recover the algebraically exact SCR method, where in the limit the algorithm converges in one iteration. We also remark that finding an optimal choice for these parameters is, to the best of the authors' knowledge, still an open question in the community.  Table 5: Influence of the intermediate and inner tolerances η t and η n (where we always set η t = η n ) on the number of outer iterations needed by the FGMRES solver, p = 2, 3, on a fixed mesh with 4096 elements.

A nine patches geometry
In this example we consider the computational domain Ω = [0, 3] × [0, 3] depicted in Figure 10, divided into nine subdomains Ω i . Similarly to the previous example, all meshes are non-conforming at every coupling interface, where again an irrational factor of √ 2/100 has been used to shift the interface knots. The body source and boundary data are derived from the following analytical exact solution: u ex = sin(πx) cos(πx) . (a) Geometry setup and physical parameters.
x y (b) Initial discretization.  Figure 11, for splines of degree p = 2, 3. In this example we test the robustness of the method with respect to: • floating patches; • the presence of multiple cross-points where a constraint must be applied. We observe again the expected asymptotic convergence rates of the error for all norms, where we remark that the method behaves optimally also for very coarse meshes, where locking phenomena are avoided. Indeed, on one hand, we notice again that a classical "vanilla" choice of the penalty parameters yield a severe overconstraint of the solution space, resulting in a loss of accuracy of several order of magnitudes compared to the projection method. On the other hand, the scaling studied in [18] leads to better results especially in the energy norm. However, for coarse meshes, we note that the method still suffers from locking, thus hindering the accuracy achievable by B-splines.

A three patches example with a geometrically non-conforming interface
In this example we consider the computational domain Ω = [0, 2] × [0, 2], split into three subdomains Ω i , see Figure 12a. The initial non-conforming discretization used in the following is depicted in Figure 12b, where the interface knots are again shifted by a factor of √ 2/100 to induce the non-conformity. The peculiarity of this example is the presence of a geometrically non-conforming interface between the patches, which is further used to assess the robustness of our method.

Remark 7.
Similarly to [8], we define an interface as geometrically conforming if the pull-back with respect to both slave and master domains is an entire edge of each parametric domain Ω i . The convergence results of the error measured in the L 2 , H 1 and H 2 are presented in Figure 13, for splines of degree p = 2, 3. Analogously to our previous results, our method attains optimal rates           Figure 13: Convergence study of the error measured in the L 2 , H 1 and H 2 norms in the non-matching case for the three patches example, B-splines of degree p = 2, 3. Comparison of a classic penalty method, the scaled version with respect to the problem parameters proposed in [18] (scaled) and our projection approach (proj).
of convergence, even in the presence of a geometrically non-conforming interface. Moreover, this numerical experiment confirms again that our method is insensitive to locking, starting from very coarse discretizations, where a substantial gain in accuracy per degree-of-freedom is observed.

A flat L-bracket
The last example we present is meant to show the applicability of the method to more complex multi-patch geometries. Analogously to the example studied in [6], we modeled a flat L-bracket with 28 patches, coupled along 34 interfaces, as depicted in Figure 14.  Figure 15, where we remark the smoothness of the obtained solution, especially across the coupling interfaces. In Figure 16 we also plot the bending stress tensor m, where its components are defined as: and where δ ij denotes the standard Kronecker delta. We obtain again a smooth stress field, where no visible spurious oscillations are introduced by the proposed coupling strategy. Finally, in Figure 17, we plot the convergence results of the stress component m 11 , evaluated at point A marked in Figure 14a, as a function of the number of dofs on a series of uniformly refined meshes. We note that for the classical penalty approach, and only for this example, we have tuned the penalty parameters to converge towards the reference value, where we have set α defl = 10 4 E , α rot = E, = 1, . . . , L. This example highlights once again the gain in accuracy achieved on coarse meshes by the proposed method, also for point-wise quantities of interest.        Figure 12a, for the flat L-bracket example for different B-splines of degree p = 2, 3. Comparison of a classic penalty method, the scaled version with respect to the problem parameters proposed in [18] (scaled) and our projection approach (proj).

Conclusions
In this work we have introduced a simple methodology for the C 1 -coupling of isogeometric patches based on the L 2 -projection of suitable super-penalty terms in the context of Kirchhoff plates. The method does not suffer from locking phenomena, even in the case of severe non-matching discretization, where optimal rates of convergence of the error measured in the L 2 , H 1 and H 2 norms have been attained also on very coarse meshes and a substantial gain in accuracy per degree-of-freedom has been observed compared to a classical penalty approach and to the scaled choice of parameters presented in [18] in the scope of Kirchhoff-Love shells. The method turns out to be particularly effective for moderate spline degrees p = 2, 3. Our choice of parameters is completely determined by the problem definition and is based upon the underlying perturbed saddle point formulation associated to the plate, from which the two Lagrange multipliers are eliminated and the magnitude of the corresponding perturbations gives us insights on how to appropriately select the penalty factors. Then, to mitigate the detrimental effects of this choice on the condition number of the system matrix, we have combined the nested block preconditioner introduced in [28] with a preconditioner based on the Fast Diagonalization algorithm tailored for isogeometric Kirchhoff plates, inspired by the strategy in [38].
To conclude, we have demonstrated numerically the applicability and robustness of the proposed projected super-penalty approach for coupling Kirchhoff plates discretized by non-conforming isogeometric patches, where the method does not show any locking also on very coarse meshes.