Patch coupling in isogeometric analysis of solids in boundary representation using a mortar approach

This contribution is concerned with a coupling approach for nonconforming NURBS patches in the framework of an isogeometric formulation for solids in boundary representation. The boundary representation modeling technique in CAD is the starting point of this approach. We parameterize the solid according to the scaled boundary finite element method and employ NURBS basis functions for the approximation of the solution. Therefore, solid surfaces consist of several sections, which can be regarded as patches and discretized independently. The main objective of this study is to derive an approach for the connection of independent sections in order to allow for local refinement and thus an accurate and efficient discretization of the computational domain. Nonconforming sections are coupled with a mortar approach within a master‐slave framework. The coupling of adjacent sections ensures the equality of mutual deformations along the interface in a weak sense and is enforced by constraining the NURBS basis functions on the interface. We apply this approach to nonlinear problems in two dimensions and compare the results with conforming discretizations.

propagate throughout the entire patch. Due to these issues, approaches for the coupling of patches are essential in the context of isogeometric analysis for solids. In regard to the discretization in the interior of the solid, mutual refinement could result in unnecessary computational cost. Local refinement is, in principle, possible for the special case of hierarchical meshes, where an exact relation exists and the control points of the refined side are a function of the control points on the coarse side. This approach has been elaborated by Kagan et al 8 for B-splines and Cottrell et al 9 for NURBS. More general approaches are, however, required in order to handle arbitrary, nonconforming meshes, and more general constraints.
The most commonly employed methods for the treatment of nonconforming discretizations are the penalty method, Nitsche's method, and the Lagrange multiplier method. For isogeometric analysis, a comprehensive comparison of these methods is given by Apostolatos et al. 10 The main advantage of the penalty method is its ease of implementation. However, the resulting accuracy is highly dependent on the choice of the penalty parameter and the condition of the stiffness matrix can be impaired severely. Recent advances focus, for example, on higher-order penalty formulations for multi-patch shells 11 and plates. 12 On the other hand, Nitsche's method is independent of empirical parameters and requires no additional degrees of freedom. The stabilization term of the Nitsche formulation requires a suitable parameter that can be bounded by the solution of a small eigenvalue problem. [13][14][15] However, the variational formulation is altered, which imposes challenges from the implementation viewpoint. Further works focus, for example, on the use of the Nitsche's method for NURBS patches 16,17 and shell structures. [18][19][20] Moreover, the Lagrange multiplier approach has been widely applied in the context of mortar methods. [21][22][23][24][25] Due to its saddle point nature, this approach is highly dependent on the choice of appropriate Lagrange multiplier spaces. Different choices for isogeometric analysis are discussed by Brivadis et al. 26 Alternatively, the function spaces can be constrained, which leads to a positive definite problem as referred to by Bernardi et al. 27 This approach is fully numerical and does not alter the variational formulation, which renders its implementation straightforward. Moreover, it forms the basis of the weak substitution method, which is proposed for the coupling of NURBS patches in isogeometric analysis. 28 The use of dual basis functions and approximate dual basis functions makes this approach competitive to conforming discretizations in terms of the computational cost. 29,30 This contribution deals with the coupling of patches for isogeometric analysis of solids in boundary representation. An isogeometric framework for the analysis of solids in boundary representation is presented in our previous works for linear 31,32 and nonlinear problems. 33,34 In the latter, a general solution approach based on the principle of virtual work is adopted instead of the scaled boundary finite element equation for linear elasticity. 35,36 This enables the treatment of nonlinearities. For the sake of conciseness, we consider here the 2D case. Due to the scaled boundary parameterization, the domain is partitioned into sections. In order to allow local refinement, adjacent sections may be nonconforming. The coupling of the individual sections is realized by a mortar approach following the idea of the weak substitution method, 28 which is applicable to the general case of nonconforming discretizations. The main difference to the weak substitution method is that here, we directly constrain the basis functions on the interface instead of coupling the global system of equations by means of static condensation. Since the applicability of this approach concerns a wide class of constraints, for example higher continuity, 37 it is the approach of choice in this work. Therefore, the coupling is realized by decomposing the NURBS basis functions on the interface based on a master-slave relation. This means that the NURBS basis functions on the slave side are expressed as a function of the NURBS basis functions on the master side in order to fulfill the equality of mutual deformations in a weak sense. Thus, the constraint is imposed in a weak manner and ensures C 0 continuity of the displacements on the interface. This approach constitutes an application of the mortar method by Bernardi et al 27 to isogeometric analysis of solids in boundary representation. We study several numerical problems in nonlinear solid mechanics. Both nonconforming and hierarchical discretizations are investigated. The main features of the proposed approach are summarized: • A coupling approach for nonconforming discretizations is derived and applied to isogeometric analysis in boundary representation. It ensures the equality of mutual deformations in a weak sense.
• The coupling is realized in a master-slave framework by constraining the NURBS basis functions on the interface. A relation is attained numerically based on the weak form of the displacement continuity and is used to decompose the slave interface basis functions.
• The approach is fully numerical and does not alter the variational formulation. All slave interface degrees of freedom are related to master interface degrees of freedom based on the decomposed basis functions.
• It enables local refinement. The application to nonlinear problems is straightforward and is demonstrated by means of several numerical examples.
The article is structured as follows: Section 2 provides a brief overview of the nonlinear isogeometric formulation in boundary representation. Section 3 presents the coupling approach based on the weak substitution method and Section 4 its extension to multiple interfaces. Furthermore, we study several nonlinear problems in order to evaluate the behavior of the connection method in Section 5. For comparison, we employ computations with conforming isogeometric discretizations in boundary representation. In Section 6 we summarize the main conclusions of our study and provide an outlook for future work.

NONLINEAR ISOGEOMETRIC FORMULATION IN BOUNDARY REPRESENTATION
Hereafter, the nonlinear isogeometric formulation in boundary representation is briefly presented. For the sake of conciseness, the formulation is presented for a 2D domain with nonlinear kinematics and elastic material behavior. A more detailed description for material and geometrical nonlinearities can be found in previous works. 33,34 The formulation employs the weak form of the boundary value problem where u is the virtual displacement, b is the body force, and P is the first Piola-Kirchhoff stress tensor in the reference configuration of the domain Ω.
The parameterization of the domain is inspired by the scaled boundary finite element method. 35,36 The domain is partitioned into sections Ω s in respect to the scaling center C. Each section is described by a circumferential parameter on the boundary and a scaling parameter , which runs from the scaling center to the boundary as depicted in Figure 1 (left). For an arbitrary point on the boundaryX or in the domain X, it holdsX = N b ( ) X s on Ω s and X =X + (X( ) −X) in Ω s , respectively.X represents the coordinates of the scaling center, N b the NURBS basis functions describing the geometry of the boundary, and X s the coordinates of the boundary control points. After transforming the geometry with the scaled boundary coordinates, the Jacobian matrix reads The differential operator in scaled boundary coordinates is only a function of the boundary coordinate and reads where J = det J( ). The submatrices b 1 and b 2 contain the components of the Jacobian matrix and are only a function of the boundary parameter . 33,34 For the numerical approximation, we employ NURBS basis functions in both parametric directions following the isogeometric paradigm. 1 An exemplary illustration for a section Ω s of the domain is given in Figure 1 (right). Besides the description of the geometry with NURBS, the displacements u are also approximated with NURBS basis functions as u( , ) = N b ( )U s ( ) for each section Ω s . Note that the displacement U s is only a function of the scaling parameter . For better illustration, the approximation in matrix notation reads Here, n bc denotes the number of boundary control points and p the polynomial degree of boundary NURBS. The dimension of U s,i for i = 1, … , n bc is 2. Moreover, NURBS are employed to approximate the displacement in scaling direction in order to allow nonlinear problems. Note that due to the straight lines in scaling direction, B-splines apply here. These are termed as N s and can be expressed for one radial scaling line in matrix form as Here, n cp denotes the number of control points per radial scaling line and q the polynomial degree in scaling direction. The vector U r,i contains the degrees of freedom along a single scaling line associated with the i-th control point on the boundary. After rearranging all control point vectors U r,i into the vector U, the displacement in scaling direction can be rewritten as where n bs = n cp ⋅ n bc is the total number of control points in Ω s . The numbering of control points in the vector U is depicted in Figure 1 (middle). The dimension of N s is (2 ⋅ n bc , 2 ⋅ n bs ). Inserting Equation (6) into Equation (4) leads to a multiplicative decomposition of the displacements, where it holds Moreover N( , ) can be expressed in matrix form as The dimension of N is 2 ⋅ n bs . It is remarked that the displacement vector U is associated to the control points of a section in the direction of . The same approximation is also adopted for the virtual displacements u. Based on the NURBS approximation, we define the discretized strain-displacement transformation matrix as where N , = N b N s , and N , = N b , N s . The discretized weak form of equilibrium can be derived and linearized as elaborated in previous works. 33,34

COUPLING APPROACH
In this section, we present the coupling of sections based on a mortar approach. 27,28 For the sake of conciseness, we derive the approach for a single interface. Let us consider the case, where two adjacent sections are nonconforming. Each section can be regarded as a patch. The interface between sections is denoted with Ω s . The deformations of each section Ω s are split into interface and domain deformations. Thus, it holds The deformation vectors u (k) c and u (k) d contain all control points on the interface and in the domain, respectively. The superscript k denotes the section number. The coupling conditions along the interface read where the traction t c along the interface can be regarded as a virtual load. The equality of mutual deformations on the interface is enforced in a weak sense and can be expressed as The interface deformation vector in Equation (12) is interpolated according to Equation (4) and can be rewritten as where the basis functions associated to interface degrees of freedom are rearranged in Ñ s , Ñ (k) . A master-slave relation is employed for the coupling as, for example, illustrated in Figure 2. The superscripts (ma) and (sl) denote the master and slave patch, accordingly. Note that the choice of master and slave patch is, in principle, flexible. In analogy to the interface deformations, for the traction vector it holds where Ñ (sl) concerns the basis functions associated to slave interface degrees of freedom and the vector t contains the interface tractions related to slave degrees of freedom along the radial scaling lines. These values have only virtual character and will drop out of the formulation in the following. Note that we have nonlocal support of the interface traction basis functions on the slave side. Local support can be in principle recovered with dual basis functions. 29,30 Moreover, it is remarked that the classification of patches in slave and master is done according to the criteria of the weak substitution method. 28 The discretized form of the equality of mutual deformations in Equation (12) reads The basis functions Ñ are defined as shown in Equation (13). Note that U c denotes the degrees of freedom along the interface in scaling direction, see also  1] of the master side. Standard Gauss quadrature is employed for the integration. The computation of the mortar matrices follows the approach of the weak substitution method. 28 Several projections are required, as the integrals of Equation (15) need to be evaluated in the parametric space of each patch. It can be observed in Equation (15) that the master mortar matrix D (ma) requires integration in the parameter space of the master patch, however, the integral D (ma) contains the basis functions Ñ (sl) , which are defined on Ω (sl) s . Therefore, every master interface point with parametric coordinates ( (ma) , (ma) ) has to be projected onto its counterpart on the slave side with parametric coordinates ( (sl) , (sl) ) in order to evaluate Ñ (sl) . The projection is performed in the physical space by determining the physical location of X (ma) for each master interface point ( (ma) , (ma) ) and then computing the physical point X (sl) on the slave side with closest point projection. The parametric coordinates ( (sl) , (sl) ) of the point X (sl) are then determined using the point inversion algorithm. Note that watertightness is assumed here, therefore it should hold that X (sl) = X (ma) . More details on the employed algorithms are given by Piegl and Tiller. 38 After employing the mortar matrices and dropping the virtual traction from the equation, the relation on the interface can be finally rewritten as The relation matrix T is obtained from the weak form of the displacement continuity in Equation (15) and can be used to enforce the constraint in a weak sense. Considering Equation (16), the slave interface degrees of freedom can be expressed as a function of the master interface degrees of freedom and substituted based on this relation. Thus, the interface deformation vector on the slave side can be rewritten considering Equation (13) and Equation (16) as In this work we impose the constraint on the basis functions. Therefore, the slave basis functions having influence on the interface are decomposed using the relation matrix T. The modified basis functions are assigned to the adjacent master degrees of freedom on the interface. This yields a relation between the basis functions of the master and slave side in order to fulfill the constraint in a weak sense. As a result, the slave interface degrees of freedom are related to the master interface degrees of freedom by means of the decomposed basis functions. Thereafter, the deformation vector of each section can be expressed considering Equation (13) and Equation (17) as where N d contains the basis functions associated to the domain degrees of freedom U d of each section Ω s . For better illustration, the basis function matrix is expressed as bs is the total number of control points in the master section and n (ma) cp denotes the number of master interface control points. The basis functions in Equation (19) are computed according to Equation (8). This results in the combined contributions of N (ma) d and Ñ (ma) , which are contained in the left block of the matrix. Note that the master interface basis functions Ñ (ma) remain unchanged, therefore 2 ⋅ n (ma) cp entries are set to zero. For the slave patch, an example is given for the case where only one degree of freedom is substituted. The basis function matrix is then modified as For better illustration, it is assumed here that only N 1 is associated with a slave degree of freedom on the interface. The dimension of N (sl) is 2 ⋅ n (sl) bs + 2 ⋅ n (ma) cp , where n (sl) bs is the total number of control points in the slave section. The slave interface basis functions are decomposed by Ñ (sl) T and assigned to n (ma) cp master interface control points. The first block of the matrix contains only the contributions of N (sl) d , whereas slave interface basis functions are set to zero. The transformation matrix for the interface can be expressed as The dimension of T is n (sl) cp ⋅ n (ma) cp , where n (sl) cp is the number of slave interface control points. The derivatives of the basis functions N (sl) , , N (ma) are derived in a similar fashion in order to form the displacement gradient and the stiffness matrix. The derivatives of the slave interface control points, which are required to compute the B matrix of Equation (9), are decomposed considering Equation (17) as This approach causes the slave interface degrees of freedom U (sl) c to vanish from Equation (18) and their recovery is possible with the help of Equation (16). It is remarked that Equation (18) holds also for the position vector X and the virtual displacements u. Moreover, the basis functions approximating the body force and traction forces are decomposed likewise following Equation (17). In this way, a natural connection is established between sections of the computational domain.

EXTENSION TO MULTIPLE INTERFACES
In general, multiple interfaces are present in the computational domain as a result of the scaled boundary parameterization, which partitions the domain into an arbitrary number of sections. The coupling approach can be extended in a straightforward manner for multiple intersections within star-shaped domains. In this case, Equation (15) is rewritten as where the summation is performed over all interfaces n int until the domain is covered. Thereafter, Equation (16) can be rewritten as All interface constraints can be assembled in a global constraint matrix T G , the dimensions of which depend on the total number of interface control points. It is remarked that each section can be discretized independently, that is, the knot vector can vary between different interfaces of the computational domain. This means that the number of master and slave interface control points may differ between interfaces. Therefore, the dimension of the global constraint matrix T G is defined by the maximum number of master interface control points and the total number of slave interface control points (max(n (ma),(i) whereas entries which do not correspond to constraints are set to zero. Moreover, any combination of conforming, hierarchical and nonconforming interfaces is possible within the domain as long as a master-slave relation is established between adjacent sections. As a result of the scaled boundary parameterization, all interfaces meet at a common point, the scaling center. This constitutes a cross-point with the same physical coordinates for all interfaces. In general, cross-points can induce circular dependencies, as they can act as both master and slaves for different interfaces. This is illustrated with a multi-patch example consisting of four patches in Figure 3.
In order to circumvent this, we introduce additional degrees of freedom at the scaling center, which are then coupled with its adjacent counterparts. In this way, each control point is considered only once as master or slave, respectively, and interface constraints are linearly independent. This is illustrated in Figure 4 for a domain with four nonconforming interfaces. The repeated control points at the scaling center of each section can be linked and considered as shared degrees of freedom in the analysis. In this way, multiple interfaces can be connected successively in a natural way.
In case of non star-shaped domains, further points where multiple interfaces meet may be present. This would require an approach for the combination of interface constraints. 29,39 The treatment of linearly dependent interface constraints seems promising for complex geometries, this approach is, however, not further followed in this work.

NUMERICAL EXAMPLES
In this section we demonstrate four numerical benchmarks in order to investigate the accuracy of the coupling approach and its capability to enable local refinement. We employ different polynomial degrees of NURBS and discretization schemes in order to evaluate the accuracy of the connection method. The focus is placed on nonlinear problems with The objective here is to investigate the behavior of local refinement for strongly nonlinear behavior and more complex geometries, which consist of curved and arbitrary boundaries accordingly. For this purpose, hierarchical discretizations are employed. The results of the computations using the coupling approach are compared to conforming computations using the isogeometric formulation in boundary representation, where adjacent sections are directly connected with shared degrees of freedom. In cases where no reference solution is available, we employ precise isogeometric computations as reference. In order to evaluate the accuracy of the connection method, we employ the L 2 -error norm for the error of the Dirichlet interface condition ‖ ‖ u (ma) − u (sl) ‖ ‖ Ω s and the Neumann interface condition ‖ ‖ t (ma) + t (sl)‖ ‖ Ω s . The error is evaluated by integration based on the mutual deformations and tractions along the interface. Conforming computations are denoted as conf. in the legend. Computations with the proposed approach are denoted as SF and MF depending on the choice of the master patch as the patch with less or more interface control points, respectively. The presented approach has been implemented into a customized isogeometric version of the Finite Element Program FEAP. 40 Table 1 summarizes the declarations for the description of the computational models.

Cook's membrane
The first example is the quasi incompressible Cook's membrane as originally presented by Simo and Armero. 41 The aim of this test is to validate the coupling approach for different discretizations and polynomial degrees of NURBS. Therefore, we evaluate the convergence behavior in terms of point-wise displacements and the L 2 -error norm of the mutual deformations along the interface. Furthermore, the capability of the approach to enable local refinement is investigated. Figure 5 depicts the geometry and boundary conditions. The plate has a thickness of t = 1 mm and dimensions as shown in Figure 5.  The system is discretized by four sections as illustrated in Figure 5 (middle). In order to validate the coupling approach, a finer discretization is chosen in respective sections as shown in Figure 5 (right). The areas with higher stress concentrations, which are depicted in Figure 5 (middle), are discretized finer. This leads to a total of four nonconforming interfaces in the domain which are numbered as shown in Figure 5 (right). Note that a singularity is present at the top of the clamped edge, therefore refinement of the clamped section is not employed here. We employ two different discretization schemes. In the first scheme, uneven sections are discretized with n (i) b = n (i) c = j and even sections are discretized finer with We remark that k-refinement is employed for all computations. The finer discretized sections are considered as slaves (SF) and also as masters (MF). To validate the connection method, we perform the analysis while monitoring the vertical displacement v A at point A as shown in Figure 5. Reference results with conforming computations are available. 34 In order to evaluate the relative error in the vertical displacement, we compute a reference solution with IGA using sixth-order basis functions and a 64 × 64 element mesh. This yields a system with 9660 equations and v A = 6.92192 mm. 34 We observe a quadratic convergence of the Newton-Raphson scheme, which is in all cases comparable to conforming computations and thus omitted here. The deformation results are illustrated in Figure 6 for the polynomial degrees p = 2 and p = 4 of NURBS basis functions. We observe that the relative error decreases with mesh refinement in all cases. For the first discretization scheme and quadratic basis functions, the results of SF are comparable to conforming computations whereas the results of MF are marginally more accurate for coarse discretizations. The second discretization scheme requires also slightly less elements to produce as accurate results as conforming computations for the same refinement level in the coarse range of discretizations, which indicates that although the number of elements in scaling direction is reduced this does not affect significantly the accuracy. For higher polynomial degree only the results of the first discretization scheme are illustrated here. The results of SF for quartic basis functions nearly coincide with conforming computations while MF is marginally more accurate for coarse discretizations. This is due to the fact that the master patch is discretized finer in MF, which may lead to more accurate results with less elements in total compared to conforming computations. For finer discretizations the results are nearly identical for all respective discretizations and polynomial degrees. This indicates that the accuracy of the solution is not impaired by the coupling and the convergence behavior is in general comparable to conforming computations.
Moreover, the error of mutual deformations is depicted for each interface in Figure 7 in order to compare the solution of SF and MF. The error for conforming computations is not depicted here as it is in the range of numerical precision. The results of the second discretization scheme are displayed for quadratic NURBS. In all cases the error decreases as the mesh   is refined, reaching a magnitude of around 10 −5 for all interfaces. Note that SF reproduces the exact substitution relation for the coarsest mesh, which is a hierarchical mesh. 8,9 Therefore, the results of SF for the coarsest mesh are in the range of numerical precision for all interfaces. The error of SF is slightly lower than MF for interfaces 1 and 3. Overall, the results of both SF and MF for finer discretizations are nearly identical. Moreover, the slope of MF is in general comparable to SF if the same interface is concerned. The results indicate that SF and MF are in general comparable, with SF being slightly more accurate in terms of the mutual deformations along the interface. This is also consistent with previous observations of the weak substitution method. 28 Furthermore, the capability for local refinement is investigated. For comparison, we employ conforming computations with uniform refinement. 34 Moreover, an isogeometric solution computed using a single patch with p = 2 and p = 4 is employed as reference. The system is refined locally in the area where the displacement response is higher, see Figure 8. The locally refined section is discretized with n b = n c = 2i and all other sections with n b = 2i, n c = i, where i = 1, 2, 4, 8, 16, 32 indicates the refinement level. A representative example of this discretization is depicted in Figure 8 for i = 2. The polynomial degree is chosen on the boundary as p b = 2 and in scaling direction as p c = 4. For better illustration the results of uniform refinement with p b = p c = 2 and the isogeometric solution with p = 2 are also depicted. Note that the chosen discretization rule produces a hierarchical mesh, as the knot vector of the locally refined section is a knot inserted version of its adjacent counterpart while the polynomial degree remains the same. In this case, an exact substitution relation exists 8,9 and can be reproduced with SF. Thereafter, we compare the response between SF and MF. The results are depicted in Figure 9. Furthermore, increasing the polynomial degree in scaling direction improves the accuracy for both uniform and local refinement. Local refinement leads to more accurate results for all respective discretizations. Increasing the polynomial degree in both directions leads to more accurate results for all discretizations. Uniform refinement is marginally more accurate than local refinement in the fine range of respective discretizations, whereas in both cases the results are slightly more accurate than the isogeometric reference solution. We observe that in terms of the global behavior, SF and MF are identical. To sum up, the proposed approach has been validated for different discretizations and can be used for local refinement.

Perforated square plate
The perforated square plate is presented here in order to investigate the performance of local refinement for arbitrary nonconforming meshes using the presented coupling approach. A reference solution with nonlinear isogeometric computations is employed for the evaluation. 28 We employ computations with conforming discretizations for comparison. 34 Figure 10 depicts the geometry and boundary conditions. The thickness of the plate is t = 1 m. All dimensions are in [m] and depicted in Figure 10. The traction is applied uniformly on the left side of the plate, see also Figure 10. We employ the St. Venant-Kirchhoff material model with the following properties: Young's modulus E = 100 Pa and Poisson's ratio = 0.3. Plane stress condition is assumed. The domain is partitioned into five sections in relation to the scaling center, which results in a total of five interfaces depicted with dashed lines in Figure 10 (left). We apply local refinement and employ a finer discretization in the areas around the hole where the targeted deformation u A is evaluated. Note that the higher stress concentrations are also expected in this area, see Figure 10 (right). The discretization results in two nonconforming interfaces, which are denoted with blue color and numbered according to Figure 10 (right). The remaining sections are discretized conforming, that is, three from a total of five interfaces are conforming. The conforming sections are discretized with n b = n c = j and their adjacent sections are discretized finer with n b = j, n c = j + 1, where j = 2, 4, 8, 16, 32 indicates the refinement level. A representative example of the discretization is depicted in Figure 10 (middle) for j = 4. We perform nonlinear analysis and monitor the horizontal displacement u A at point A (see Figure 10 (left)). Note that point A is located on the physical structure and not on the control net. Thus, it is evaluated by substituting the known control variables in Equation (7). The relative error is computed with the reference solution u A,ref = 0.342816274990393. 28 The convergence of the displacements at point A can be seen in Figure 11 for polynomial degree p = 2 and p = 4. We observe that the error decreases as the mesh is refined. SF is slightly more accurate for coarse quadratic discretizations than conforming computations. MF is more accurate than SF with the exception of a kink where the results seem to be unstable. For finer discretizations, both cases coincide very well with conforming computations. For higher polynomial degree, SF and MF are more accurate for coarse discretizations and less accurate than conforming computations for finer discretizations. The convergence rate of SF is in most cases comparable to conforming computations, whereas MF converges slower than conforming computations.
Furthermore, we evaluate the global energy norm L E = ||S h ∶ E h || in order to assess the global stress behavior. The relative error is computed with the reference solution L E,ref = 3.896722804022. 28 Figure 12 depicts the convergence of the energy norm. We observe that the locally refined quadratic discretizations produce more accurate results than conforming computations with uniform refinement, as we have locally refined the areas with higher stress concentrations. Note that the relative error for MF is lower, which is also consistent with previous observations of the weak substitution method. 28 For discretizations in the fine range both cases agree very well with conforming computations. For higher polynomial degree, SF is more accurate for coarse discretizations and less accurate than conforming computations for fine discretizations. Conforming computations are more accurate than MF. The convergence rate of SF is nearly quadratic  and comparable to conforming computations for low polynomial degrees, whereas MF converges slightly slower. For higher polynomial degree and fine discretizations, conforming computations converge faster. It is remarked that the stress behavior of high-order discretizations might be additionally affected by the presence of C 0 -continuity between sections as a result of the imposition of the displacement continuity constraint. The equivalent stress distribution is depicted in Figure 13. Overall, a good agreement is observed between SF and MF. In both cases, jumps are present along interface 1, whereas MF seems to capture the stress response along interface 2 more smoothly than SF. Moreover, we further compare SF and MF for p = 2 and p = 4 by investigating the L 2 -error norm of the Dirichlet and Neumann interface conditions in Figures 14 and 15, respectively. For conciseness only the error of nonconforming interfaces is depicted here. We observe that in all cases the error decreases as the mesh is refined, reaching a magnitude of up to 10 −10 for the Dirichlet interface condition and up to 10 −6 for the Neumann interface condition for the second interface. For p = 2, the error for interface 1 has the same order of magnitude for both SF and MF, while SF is slightly more accurate than MF. The slope is also similar for both solutions. The error of SF in interface 2 is one order of magnitude higher than MF. SF seems to perform better for interface 1, whereas MF is more accurate for interface 2. This is also the case for p = 4, where it can be observed that changing the choice of the master section leads to more accurate results than increasing the polynomial degree. Note that the stresses are in general higher for interface 2 due to transverse deformation and MF leads to a smoother stress distribution in this area (cf. Figure 13).
Overall, it can be said that the proposed approach can be used for the local refinement with arbitrary, nonconforming interfaces. Both conforming and nonconforming interfaces can be coupled inside the computational domain. Regarding the classification of patches, it cannot be properly said which choice is preferable. The results of both the SF and the MF choice are as good as the results of conforming computations. There is no significant negative impact from the nonconforming discretization on the global stress error behavior for low polynomial degrees, whereas high-order discretizations might benefit from higher continuity measures between adjacent sections.

Propeller blade
The objective of this example is to apply the coupling approach on a complex geometry with curved boundaries. The employed system and boundary conditions are depicted in Figure 16. The system is fixed in radial direction as shown in Figure 16. The thickness is t = 1 m. More details on the geometry of the system can be found in previous studies. 33 We assume here a constant rotational velocity of = 600 rpm and density = 1 kg/m 2 . The components of the body load are defined according to Barber 42 as This results in body forces with the value of p x = −100 N/m 2 and p y = −100 N/m 2 . We employ the St. Venant-Kirchhoff material model under plane stress condition with the following material properties: Young's modulus E = 1 MPa and the Poisson's ratio = 0.3. Due to symmetry, only one blade of the propeller system is modeled. The blade is partitioned into five sections. We apply local refinement in the areas where the radial deformations are higher, see Figure 16 (right). It is remarked that the deformation is depicted by evaluating the NURBS basis functions and plotting on the physical structure. Our aim is to investigate the convergence of deformations using local refinement and the presented coupling approach. For comparison, we employ computations of conforming discretizations with p = 2, 3. The locally refined section is discretized with n b = n c = 2i and all other sections with n b = 2i, n c = i, where i = 1, 2, 4, 8, 16, 32 indicates the refinement level. A representative example of the resulting hierarchical discretization is depicted in Figure 16 (right) for i = 1. It is remarked that the basis functions interpolating the body forces are obtained as shown in Equations (19) and (20). We perform the analysis and monitor the radial deformation at point A, see Figure 16 (left). To the authors' knowledge, there is no analytical solution for nonlinear computations of this geometry. Therefore, we perform precise computations with p = 6 and 20 480 elements using the isogeometric formulation in boundary representation. This yields a radial displacement of u r,A = 0.01769696835 m. Figure 17 depicts the convergence of the radial deformation at point A. We observe that the error decreases as the mesh is refined. SF results in lower error for all discretizations. Note that SF reproduces the exact substitution relation and the error of the Dirichlet interface condition is in the range of numerical precision for all discretizations, therefore it is omitted here. Moreover, the stress distribution is compared between local and uniform refinement for a mesh with p = 2 and i = 4 ( Figure 18). Due to the C 0 -continuity between sections, jumps can be observed in both cases. For the locally refined mesh, the jumps are slightly higher due to the weakly imposed displacement continuity constraint. Overall, the stresses agree very well. It can be concluded that the global behavior is not impaired by the coupling approach.

Flower
The last example concerns a smooth geometry with flower-like boundaries. The aim of this benchmark is to investigate the behavior of local refinement on an arbitrary domain with strong nonlinearities. Figure 19 depicts the geometry and boundary conditions. The geometry is adopted from the literature. 43 The system is fixed along the Dirichlet boundary u Ω and loaded with a constant traction along the Neumann boundary t Ω as shown in Figure 19. The thickness is t = 1 m. In order to account for large strains, we employ the Neo-Hookean material model with shear modulus = 38.4615 Pa and Lamé parameter Λ = 57.6923 Pa. Figure 20 (left) depicts the distorted mesh for the conforming case. Local refinement is applied in the area with higher distortions, see Figure 20 (right). Similarly to the previous example, the deformation is depicted by evaluating the NURBS basis functions and plotting on the physical structure. The locally refined area is discretized with n b = n c = 2i and all other sections with n b = 2i, n c = i, where i = 1, 2, 4, 8, 16 indicates the refinement level. This results in hierarchical meshes in analogy to the previous example. A representative example of the local refinement for i = 2 is depicted in Figure 16 (right). The polynomial degree of NURBS is chosen p = 3 for all computations. The deformation u is monitored, see also Figure 19. Note that the point under consideration is located on the curve and not on the control net. Therefore, after obtaining the solution of the control variables, it has to be further evaluated using Equation (7). To the authors' knowledge, there is no analytical solution for nonlinear computations of this geometry. Therefore, we perform precise computations with p = 6 and 9344 elements using the isogeometric formulation in boundary representation. This yields a horizontal displacement of u x,ref = 2.67859352106265 m and vertical displacement of u y,ref = 1.02328660586063 m . Figure 21 depicts the relative error in displacements. We observe that the results of SF and MF are identical in terms of the global behavior. It is noted that SF reproduces the exact substitution relation for all discretizations whereas for MF no exact relation exists. 8 This means that an error in the Dirichlet interface condition may be present in MF, which however does not seem to impair the global behavior. The results in Figure 21 indicate that local refinement achieves the same accuracy as uniform refinement with less elements in total for all respective discretizations. Overall, the proposed approach enables the independent discretization of sections with arbitrary boundaries.

CONCLUSIONS
This contribution deals with an approach to couple nonconforming patches for isogeometric analysis of solids in boundary representation. The motivation of this work is to perform local refinement in the interior of the computational domain. Due to the adopted scaled boundary parameterization, the domain is described by a circumferential parameter on the boundary and a scaling parameter in the interior of the domain. Thus, the domain is decomposed into sections, which can be understood as patches and discretized independently, that is, nonconforming to each other. For this reason we propose a coupling approach based on the weak substitution method. It is remarked that this approach is chosen here mainly due to its straightforward treatment of nonlinearities and multiple interfaces. Also other domain decomposition methods could in principle be used to enable local refinement. A master-slave relation is established where the degrees of freedom on the slave side can be expressed as a function of the degrees of freedom on the master side. The coupling is realized in a weak manner by constraining the basis functions on the interface. The slave basis functions on the interface are decomposed using the relation matrix and assigned to the corresponding master interface control points. In this way C 0 -continuity of the displacements is obtained along the interface. This procedure is fully numerical and applicable to nonlinear problems. We have applied this approach to several nonlinear problems in solid mechanics. The results indicate that the coupling does not impair the global solution and produces in general results of comparable accuracy with conforming computations. Moreover, local refinement produces more efficient results in areas with high distortions.
Regarding the master-slave relation, choosing the coarser side as master seems in general to be more accurate for the Dirichlet interface condition whereas choosing the finer side as master seems more accurate for the global stress behavior. This choice is in principle flexible and depends on the targeted computational accuracy and efficiency. The main advantages of the presented approach are summarized: • The coupling approach does not alter the variational formulation and can be applied to nonlinear problems without further measures.
• It enables local refinement of the computational domain and reduces unnecessary computational cost caused by mutual refinement of adjacent sections.
• It can be broadly applied in its current form for different types of problems as it is easily extendable to more general constraints.
Future work concerns the extension of this approach to ensure higher order continuity along the interface and thus smooth solutions. 25,44 Currently, only 2D problems have been considered. Further work needs to be done regarding the application to the 3D case. Then, the presented approach could be used in the future to couple nonconforming surface patches which have been designed in CAD with the boundary representation modeling technique.