Non-symmetric isogeometric FEM-BEM couplings

We present a coupling of the Finite Element and the Boundary Element Method in an isogeometric framework to approximate either two-dimensional Laplace interface problems or boundary value problems consisting in two disjoint domains. We consider the Finite Element Method in the bounded domains to simulate possibly non-linear materials. The Boundary Element Method is applied in unbounded or thin domains where the material behavior is linear. The isogeometric framework allows to combine different design and analysis tools: first, we consider the same type of NURBS parameterizations for an exact geometry representation and second, we use the numerical analysis for the Galerkin approximation. Moreover, it facilitates to perform h- and p-refinements. For the sake of analysis, we consider the framework of strongly monotone and Lipschitz continuous operators to ensure well-posedness of the coupled system. Furthermore, we provide an a priori error estimate. We additionally show an improved convergence behavior for the errors in functionals of the solution that may double the rate under certain assumptions. Numerical examples conclude the work which illustrate the theoretical results.


Introduction and preliminaries
In the last decades, simulation gained more and more importance as a forth pillar of sciences besides theory, experiments, and observations. A successful simulation means a good imitation of some phenomena. This allows the analysis, optimization, and predictions to be made ad-hoc with a certain reliability, which depends on the application. For example, this can be achieved for problems that are formulated as boundary and initial value problems by choosing the right mathematical model, a good representation of the computational domain, and a suitable numerical method. We encounter in this work two types of model problems. First, we consider a Laplacian interface problem in Section 2. Its specificity lies in the combination of a possibly non-linear and non-homogeneous problem in a bounded domain with a linear and homogeneous problem in an unbounded domain. This type of model describes a wide class of engineering and physical applications. One example for electromagnetic scattering problems and elastostatics can be found in [Ste11]. We address the second type of model in Section 4. It is a Boundary Value Problem (BVP) with two disjoint domains, which are separated by a (thin) gap. In the domains we allow non-linear equations. However, the gap is assumed to be filled with a linear material, where the simplest form is air. For a visualization we refer to Figure 1. This model is particularly used for the simulation of electro-mechanical energy converters. An example is an electric machine discussed in [BCSDG17]. In general, the air gap is very thin. The other two domains are modeled separately for a facilitation of a possible rotation of the interior part, called rotor in case of an electric machine. This movement is induced by the interaction of electromagnetic fields in the air gap. The computation of forces and torques are therefore one central goal in this type of simulations. Formally, this can be achieved by using the so called Maxwell Stress Tensor (MST) method, see, e.g., [KFK + 97]. For this, the solution in the air gap as well as its derivatives are needed. These aspects have to be kept in mind for a suitable choice of a numerical method. The coupling of the Finite Element Method (FEM) and the Boundary Element Methods (BEM) appears to be an intuitive and straightforward choice for the above described problems. Indeed, the FEM is well established and widely used for possibly non-linear problems in bounded domains. On the other side the BEM relies on the transfer of the model problem to an integral representation. Further steps then lead to a Galerkin discretization problem on its boundary with certain integral operators. In a post-processing step, a solution can be found in every point of the underlying domain. Hence, BEM is suitable to handle problems with an unbounded domain where we do not have to truncate the domain since the discretization itself is done on the boundary. We remark that a truncation would be mandatory if we would apply FEM. Since the BEM discretization takes place on a boundary of the domain, it is also very attractive to get a solution in the thin gap described above. For a mesh-based method in the whole domain, e.g., like FEM, it is very difficult to find a mesh for such a thin gap, where the numerical method remains stable. The discretization on the boundary with BEM and the post-processing afterwards avoids this problem. However, to apply the BEM we need to know the fundamental solution of the underlying problem. Note that this restricts the application of BEM especially for non-linear problems. Therefore, we apply BEM in this work for two different applications: first in the exterior unbounded domain and second in the thin air gap. In both cases we consider the Laplace operator for the BEM part, where the fundamental solution can be given explicitly.
In the literature, we distinguish several types of FEM-BEM coupling techniques. These coupling procedures differ solely in the considered representation of the Boundary Integral Equations (BIE), which are the basis for BEM. In order to introduce briefly the considered BIEs, we envisage first the following Laplace equation (1) − ∆u = 0 in Ω κ , κ = 0, 1, where Ω 0 ⊆ R 2 is a bounded domain with Lipschitz boundary Γ, and Ω 1 = R 2 \Ω 0 is the corresponding unbounded domain. Hence, Γ = Ω 0 ∩ Ω 1 . Note that (1) is an interior problem for κ = 0 and an exterior problem for κ = 1. In the latter case, we additionally assume the radiation condition u(x) = C ∞ log |x| + O(1/|x|) for |x| → ∞ with the unknown constant C ∞ , see also Remark 2.1. For some x ∈ Ω κ , the solution u(x) is given by the representation formula where G(x, y) = − 1 2π log |x − y| denotes the fundamental solution of the Laplace operator, ν(y) is an outer normal vector on Γ pointing outward with respect to Ω 0 at y ∈ Γ, and u |Γ , φ := ∂ ν u |Γ are the unknown or partially unknown Cauchy data. Hereby, the notation u |Γ means the trace of u with respect to Γ. Note that we omit to write the trace operators in this work due to readability. Taking the trace of the representation formula yield to the following BIEs, see, e.g., [Ste07, Chapter 7] for more details, The invoked Boundary Integral Operators (BIO), the single layer operator V and the double layer operator K, are given for smooth enough inputs by and can be extended continuously to linear and bounded operators such that Theorem 1]. In particular, the boundary integral operator V is additionally symmetric, and H − 1 2 (Γ)-elliptic, if diam (Ω) < 1, see, e.g., [Ste07,Theorem 6.23]. The properties of V induce the norm equivalence .
In the previous lines, the mentioned spaces have to be understood as follows: for k > 0, H k (·) denotes the standard Sobolev space equipped with the usual norm · H k (·) . Moreover, the space H k− 1 2 (Γ) is the trace space of H k (Ω), and spaces with negative exponents H −k (Γ) are defined as dual spaces of H k (·) using the natural duality pairing ·, · Γ , which is obtained by the extended L 2 -scalar product (·, ·) Γ . Furthermore, for the unbounded domain Ω 1 we need functions with local behavior and denote them by H 1 loc (Ω 1 ) := {v : Ω 1 → R v |K ∈ H 1 (K) for all K ⊂ Ω 1 compact}. Finally, we write H 1 (Ω) ′ for the dual space of H 1 (Ω). We recall that holds for all v ∈ H 1 (Ω) and ψ ∈ H − 1 2 (Γ), where the trace inequality is encoded with the trace constant C tr > 0. For some bounded domain Ω, (·, ·) Ω denotes equivalently the standard L 2 -scalar product in Ω.
In the following, we describe a variational ansatz to get a weak form of the model problem. As mentioned above, there are several coupling strategies possible. If we describe the FEM part by the weak form of the well-known Green's first formula, a coupling with the weak form of (3) (κ = 1) leads to the so called Johnson-Nédélec coupling introduced in [JN80]. The combined weak form is non-symmetric even though the model problem itself is symmetric. Also a Galerkin discretization leads to a non-symmetric system of linear equations. Therefore, this coupling is also known as non-symmetric coupling, where the unknowns are the u of the FEM part and the conormal derivate φ of the BEM part. To symmetrize this system in case of a symmetric model problem, we first observe that taking the conormal derivative of (2) leads to another integral equation with two other integral operators. A modification of the Johnson-Nédélec coupling with this additional integral equation renders the coupled problem symmetric. This procedure appeared first in [Cos88b] and is known as Costabel's symmetric coupling. The price of the symmetry is the use of four BIOs, which is computationally more expensive. However, there are still only two unknowns involved. A coupling method with three unknowns, i.e., additionally the trace u |Γ of the BEM part is an unknown, is called a three field coupling [Era12]. A coupling procedure with the so called indirect ansatz is also possible and is called Bielak-MacCamy coupling [BM83]. With this strategy, however, one unknown of the BEM part has no physical meaning. Because of the advantages of the non-symmetric coupling, we consider in this work only this type of coupling. We will introduce it formally in Section 2. For a long time a mathematical analysis for this coupling was only available for smooth boundaries due to the use of a compactness argument of the double layer operator K; [JN80]. In particular, Lipschitz boundaries were excluded . However, a decade ago Sayas [Say09] provided in fact the first analysis also for Lipschitz boundaries. This work influenced several variations and improvements, e.g., [Ste11, AFF + 13, EOS17] to mention a few but not all. Hence, the non-symmetric coupling became a more natural choice, especially, if a part of the model problem is non-symmetric or non-linear. In this work, we use the results of [AFF + 13], where an extension to nonlinear interface problems has been addressed and combine this result with the proof shown in [EOS17]. Furthermore, we also use [OS14], which extended the proofs for the linear interface problem to certain Boundary Value Problems, i.e, to the second type of problem considered here. For the linear interface problem with a general second order problem in the interior domain, we also refer to [EOS17] for a rigorous and, to the authors knowledge, sharpest ellipticity estimate. Recently, a complete analysis of a parabolic-elliptic interface problem with a full discretization in the sense of a non-symmetric FEM-BEM coupling for spatial discretization was published in [EES18]. Note that such a system arises, for instance, in the modeling of eddy currents in the magneto-quasi-static regime [MS87]. Now, having described the weak form of the model problem with the proposed FEM and BEM parts, we still need to take two major decisions for a successful simulation: a suitable discretization technique, i.e., choosing concrete ansatz spaces for the FEM and BEM, and a good representation of the geometry. These steps are typically made independently, which complicates meshing and remeshing procedures without altering the original geometry. In order to circumvent this, design step and numerical analysis can be combined by considering the same type of basis functions. Hence the geometrical modeling is also used to design ansatz functions in the Galerkin discretization schemes for the approximation of the solution. Such a method is proposed in [HCB05,CHB09]. It is based on using Non-Uniform Rational B-Splines (NURBS) for the unification of Computer Aided Design (CAD) and Finite Element Analysis (FEA). This method is called IsoGeometric Analysis (IGA). The first isogeometric BEM simulation of collocation type can be found in [PGK + 09, SSE + 13]. Moreover, fast methods for isogeometric BEM have been successfully implemented in [HR09,MZBF15,DKSW19], which reduces the known high computational complexity of such an application due to the dense matrices produced by the BEM. This makes the method more attractive even for more realistic and complex applications, see, e.g., [CdFDGS16] and [BCSDG17]. A rigorous mathematical analysis for isogemetric FEM started in [BBadVC + 06, BadVBSV14] and for isogemetric Galerkin BEM in [Gan14,FGP15,FGHP16,Gan17,FGPS19]. For our purpose the results [BadVBSV14, BDK + 20] together with [AFF + 13, EOS17, OS14] play a central role in proving the validity and an a priori error estimate of the FEM-BEM coupling in the isogeometric context, which is done in this manuscript for the first time.
The rest of this paper is organized as follows: in Section 2, the non-linear interface problem is addressed. We consider the framework of Lipschitz continuous and strongly monotone operators such as given in [Zei86] and used in [AFF + 13]. Strong monotonicity of the non-symmetric weak form is showed equivalently to [EOS17] by adapting the setting to non-linear operators. Moreover, well-posedness of the coupling is stated. Section 3 is devoted to the Galerkin discretization of the non-symmetric coupling. Thereby, we introduce the isogeometric framework and the necessary discrete spaces. We derive some error estimates for the conforming isogeometric discretization. In Section 4, we extend the model to a Boundary Value Problem. More precisely, the model domain is split in two disjoint domains, which are separated by a thin (air) gap. First, a variational formulation of the coupled problem is derived. Then we show well-posedness and stability of the method. Furthermore, we discuss a super-convergence result for the evaluation of the solution in the BEM domain. In the last Section 5, we confirm the theoretical results by conducting one numerical example for each model problem. The work is completed by some conclusions and an outlook.

Interface problem
Let Ω ⊂ R 2 be a bounded domain with Lipschitz boundary Γ = ∂Ω and Ω e := R 2 \Ω the corresponding unbounded (exterior) domain. Furthermore, to guarantee the H − 1 2 (Γ)-ellipticity of the boundary integral operator V, we assume diam (Ω) < 1. This assumption can merely be achieved by scaling. We consider the following interface problem: in Ω e , (7b) We remind that ν denotes the outer normal vector with respect to Ω and U : R 2 → R 2 is a possibly non-linear diffusion tensor. The right-hand side is given by f ∈ H 1 (Ω) ′ , u 0 ∈ H 1 2 (Γ) is the jump in the Dirichlet data, and φ 0 ∈ H − 1 2 (Γ) the jump in the Neumann data. To ensure the right radiation condition (7e) at infinity, we have to assume the additional condition ∇u e |Γ · ν, 1 Γ = 0. This can be transformed into a compatibility condition on the data, i.e., (f, 1) Ω + φ 0 , 1 Γ = 0.
Remark 2.1. Note that the assumption to ensure the radiation condition is only needed in the two dimensional case. Alternatively, (7e) can be replaced by a logarithmic decay of the solution in two dimensions to avoid the additional assumption on ∇u e |Γ , i.e., u e = C log |x| + O |x| −1 , for |x| → ∞, with C := 1 2π ∇u e |Γ · ν, 1 Γ or equivalently C := − 1 2π ((f, 1) Ω + φ 0 , 1 Γ ), which can be easily verified in the weak formulation below.
As mentioned above, the diffusion tensor U can be a non-linear operator. To apply standard theory for non-linear operators, see, e.g., [Zei86], we assume throughout the manuscript that U is Lipschitz continuous and strongly monotone: (A1) Lipschitz continuity: . The derivation of a non-symmetric variational form follows a standard procedure: In the variational form of (7a), we replace the Neumann data by the jump condition (7d) to couple the interior problem with the conormal derivative with φ := ∂ ν u e |Γ = ∇u e |Γ · ν of the exterior problem. For the second equation we use the exterior integral equation (3) with κ = 1, and insert the jump condition (7c) to couple this with the interior trace.
Hence, the weak formulation of the non-symmetric coupling problem reads: . This variational form can be written in a compact form. For this we introduce a product space with corresponding norm, i.e., holds ∀v ∈ H with the linear form (linear in the second argument) a : H × H → R, and the linear functional ℓ on H, It is easy to check that a(v, v) is not elliptic, e.g., insert v = (1, 0). Hence, [AFF + 13] suggested an implicit stabilization where the stabilized problem is equivalent to the original one, i.e., a solution of the original problem is also a solution of the stabilized one and vice versa. Thus, the analysis is done with the aid of the stabilized form, i.e., well-posedness is inherited to the original problem. For implementation purposes, we still use the original problem. The stabilized problem reads: In order to state well-posedness for Problem 2.3, and thanks to Lemma 2.4 also for Problem 2.2, we follow standard results for monotone operatos [Zei86]. First, we note that the form a(u, v) induces a non-linear operator A : H → H ′ by where H ′ denotes the dual space of H. This allows us to prove the following lemma.
for all u = (u, φ) ∈ H, v = (v, ψ) ∈ H with the norm ψ 2 V := ψ, Vψ and with C stab = min 1, , then A is strongly monotone , i.e., there exists C ell > 0 such that Proof. The Lipschitz continuity of A follows from the Lipschitz continuity of U, and the continuity of the integral operators. The proof of the second assertion follows the lines of [EOS17, Theorem 1] for β = 1. We replace the coercivity estimate of the bilinear form (U∇u, ∇v) Ω considered in [EOS17] for a linear U by the strong monotonicity property of U, i.e, ell is a direct result of the use of a contractivity result for the double layer operator K [OS13, Lemma 2.1] with a constant C K ∈ [ 1 2 , 1), where we use the worst case of C K = 1 in the statement. For the last assertion we note the norm equivalence (5), and by a a Rellich compactness argument it can be shown [AFF + 13, Lemma 10] that This leads together with (12) to the last assertion.
The following theorem follows directly from the theoretical result [Zei86, Theorem 25.B].
Theorem 2.6 (Well-posedness, [Zei86]). Let C U ell > 1 4 . Since the induced operator A of a(·, ·) is strongly monotone and Lipschitz continuous, there exists a unique solution u : Thanks to Lemma 2.4, this is also the unique solution of Problem 2.2.
For some engineering applications, where the non-linear operator U has a special form, we can state the following stabilization result.
Lemma 2.7. Let all the assumptions of Theorem 2.6 hold and let the non-linear operator U be of the form U∇u := g(|∇u|)∇u with a non-linear function g : R → R. Then, for the solution u := (u, φ) ∈ H of Problem 2.2, we have the stability result Proof. Let v ∈ H be arbitrary. We know from the strong monotonicity of A that Without loss of generality, we choose v = (0, 0) and note that U∇v = 0. Thanks to Lemma 2.4 u := (u, φ) is also the unique solution of the Problem 2.3. Thus, we conclude that with s(u) := 1, 1 2 − K u |Γ Γ + 1, Vφ Γ . Next we use inequality (6) along with the boundedness of K and V. Then, rearranging the terms yield to where C K , C V > 0 denote the continuity constants of the boundary integral operators K and V, respectively. From this follows the assertion with a constant C > 0 that depends on C K , C V , C tr , C ell and Γ.
Remark 2.8. Because of Lemma 2.4, the results obtained for the stabilized formulation also hold true for the original non-symmetric coupling of Problem 2.2. Hence, in the next section, we only discretize the original problem using a Galerkin approximation. In fact, the stabilized version is only used for analysis purposes.

Galerkin discretization
Let V ℓ ⊂ H 1 (Ω) and X ℓ ⊂ H − 1 2 (Γ) be some finite dimensional subspaces, where the index ℓ expresses a refinement level, e.g., in a sequence of mesh refinements. We assume that: (A3) The discrete space X ℓ contains the constants, i.e., We consider a conforming Galerkin discretization of the Problem 2.2. Replacing the spaces H 1 (Ω) and H − 1 2 (Γ) with V ℓ and X ℓ , respectively, yields to the following discrete problem: The compact form in the product space H ℓ reads: The linear form a(·, ·) and the linear functional ℓ are defined in (9) and (10), respectively.
Provided that Assumption (A3) is satisfied, the analysis for Problem 3.1 is done analogously to the continuous Problem 2.2 since the discrete spaces are conform. In other words, all the above results including the introduction of a stabilized form and Lemma 2.4 also hold for the subspaces. In particular, due to Theorem 2.6 the discrete solution u ℓ = (u ℓ , φ ℓ ) ∈ H ℓ = V ℓ × X ℓ of Problem 3.1 exists and is unique. The following quasi-optimality result in the sense of the Céa-type Lemma is a standard but central result, which will be needed in Subsection 3.2 for the a priori error estimate of the non-symmetric coupling.
Theorem 3.2 (Quasi-optimality). Let Assumption (A3) hold, and C U ell > 1 4 . Moreover, let u := (u, φ) ∈ H be the unique solution of Problem 2.2, and u ℓ := (u ℓ , φ ℓ ) ∈ H ℓ the solution of its discrete counterpart Problem 3.1. It holds that The assertion follows as a result of the main theorem on strongly monotone operators [Zei86,Corollary 25.7]. That means with v ℓ = (v ℓ , ψ ℓ ), thanks to Lemma 2.4, the strong monotonicity, Galerkin orthogonality, Cauchy-Schwarz inequality, and the Lipschitz continuity we get where the assertion follows directly.
3.1. Isogeometric Analysis. The basis functions that are considered for the geometry design in the isogeometric framework are used as ansatz functions for the Galerkin discretization. These functions are typically B-Splines or some extensions of B-Splines, e.g., NURBS, T-Splines etc. In the following, we introduce briefly the concept of Isogeometric Analysis and refer to [CHB09] for a more detailed introduction, and to [BadVBSV14] and [BDK + 20] for a mathematical analysis of IGA in the FEM and BEM context, respectively.
Associated to the knot vector Ξ, k B-Spline basis functions can be defined recursively for p ≥ 1 by for all i = 0 . . . k − 1, starting with piecewise constant basis functions for p = 0, namely, ..k−1 } the space of B-Splines of degree p and dimension k in the parameter domain over the knot vector Ξ.
with c i ∈ R d representing an element of a set of control points. The mapping f describes a onedimensional curve embedded in a d-dimensional Euclidian space and is called a B-Spline curve. Moreover, we call γ a patch if the mapping f is regular.
As long as the B-Spline mapping f is regular, we can define B-Spline spaces in the physical domain, i.e., over a patch γ by using the following transformation Definition 3.5. B-Spline spaces on higher dimensional domains are constructed by using tensor product relationships. For example, in 2D, we write with p 1 , p 2 ∈ N and k 1 , where p 1 and p 2 denote the degrees in each parametric direction, and k 1 k 2 is the number of the B-Splines basis functions. A B-Spline surface is thus represented by and c i1,i2 ∈ R d representing an element of a set of control points. If f is regular, we call ω a patch. Equivalently, the two-dimensional B-Spline space in the physical domain is defined over a patch ω by Note that for the sake of simplicity, if p 1 = p 2 = p, a B-Spline space of degree p should be understood as a B-Spline space of degree p in each parametric direction.
The parametrization of curves and surfaces using B-Spline functions allows an exact representation of a large spectrum of geometries. However, they fail to represent conic sections exactly, which are widely present in the design of various engineering applications. In order to circumvent this, Non-Uniform Rational B-Splines (NURBS) are used instead, see [CHB09] and [PT12], for instance.
Definition 3.6. Let p, k, p 1 , p 2 , k 1 , k 2 ∈ N as above. NURBS mappings can be considered as weighted B-Spline mappings. They can be defined as follows, Thereby, w i , w i1,i2 ∈ R are elements of a vector of dimension k and a matrix of dimension k 1 × k 2 , containing weighting coefficients of the NURBS, respectively, and c i , c i1,i2 are the control points.
Remark 3.7. Contrary to B-Splines, NURBS spaces on higher dimensional domains cannot be defined using simple tensor product relationships.
In order to guarantee the existence of a regular mapping between the parameter and the physical domain, multiple patches defined through a family of regular parameterizations may in some cases be necessary.
Definition 3.8. Let Ω be a two-dimensional Lipschitz domain with boundary Γ. The domain Ω is called a multipatch domain, if there exists a family of N Ω disjoint patches such that Ω = i Ω i and a regular parametrization r i (x) : [0, 1] 2 → Ω i for every single patch Ω i , with 0 ≤ i < N Ω . Furthermore, we require the parametrization at interfaces to coincide. Equivalently, Γ is also considered a multipatch domain with Γ = i Γ i and r i (x) : Knowing that B-Splines form a partition of unity [PT12], it is easy to see that B-Splines are a special type of NURBS, when the weightings are equal to 1.
In the following, if we refer to the geometry, we mean NURBS mappings. If we refer to the spaces used for the discretizations, we mean B-Spline mappings. The motivation for this follows from [BDK + 20], namely, the spline preserving property of B-Splines is needed for a conforming discretization of the De Rham complex.
3.2. Error estimates for an isogeometric FEM-BEM discretization. Let the assumptions of Section 2 on Ω hold. We consider the discrete Problem 3.1 with V ℓ = S 0 (Ω) and X ℓ = S 2 (Γ), where S 0 (Ω) and S 2 (Γ) are B-Spline spaces defined as in [BDK + 20] and [BadVBSV14]. Namely, Thereby, N Ω and N Γ denote the number of domain patches and boundary patches, respectively. Note that the degrees of the B-Spline spaces (13) and (14) are solely fixed by one parameter p > 0. Throughout the rest of this work, we assume the following: (A4) All knot vectors are p-open and locally quasi-uniform, i.e., for all non-empty, neighboring elements [ξ i1 , ξ i1+1 ] and [ξ i2 , ξ i2+1 ], there exists θ ≥ 1, such that The multipatch geometry of Ω is generated by a family of regular, smooth parameterizations.
Proof. From [BDK + 20] we know that S 0 (Ω) and S 2 (Γ) are closed subspaces of H 1 (Ω) and H − 1 2 (Γ), respectively. Moreover, Assumption (A3) holds true per construction of the B-Spline spaces. Hence, the usual analysis for a conforming Galerkin discretization of a non-symmetric FEM-BEM coupling can be considered also in the isogeometric context. Now, using Lemma 3.11 and the quasi-optimality stated in Theorem 3.2 yields the assertion.

Extension of the model problem
Let Ω, Ω 1 , Ω b , Ω 2 ⊂ R 2 , be bounded Lipschitz domains, see Figure 1. We denote by Γ b = Γ 1 ∪ Γ 2 the boundary of Ω b and by Γ 0,1 and Γ 0,2 the Dirichlet boundaries of Ω 1 and Ω 2 , respectively. Furthermore, we define H 1 0 (Ω i , Γ 0,i ) := {u ∈ H 1 (Ω i ) : u |Γ0,i = 0} for i = 1, 2. We consider the following boundary value problem: Find (u 1 , u 2 , u Hereby, ν i and ν b denote the outer normal vector of Ω i and Ω b , respectively, (f i , u 0,i , φ 0,i ) ∈ H 1 (Ω i ) ′ × H 1 2 (Γ i ) × H − 1 2 (Γ i ) with i = 1, 2 are some given data, and U i are possibly non-linear operators with the Assumptions (A1) and (A2). We emphasize that the model problem (16) can be used to simulate electric machines, see also the example in Section 5.2, which motivates its consideration. Next, we want to derive a weak formulation for Problem (16). We consider the weak form of the two problems in Ω 1 and Ω 2 . Hence, we multiply (16a) with test functions and apply the first Green's identity and get for i = 1, 2. Note that u i = 0 on Γ 0,i . We may transfer (16b) in Ω b to an integral equation on Γ b in order to apply BEM in the following. Hence, the (interior) representation formula (2) (κ = 0) hold if we replace u by u b . Let φ := ∂ ν b u b denote the conormal derivative of u b on Γ b , the BIE is obtained as in Section 1 where the the single layer operator V and the double layer operator K are defined in (4) over Γ b instead of Γ but of course with the same fundamental solution G(x, y). Note that the normal vector ν b points outwards with respect to Ω b since it is considered as an interior problem in our integral equation notation.
In what follows we strongly follow the work of [OS14], where a boundary value problem with hard inclusion is considered. As in [OS14], we can derive two equivalent weak formulations. It is enough to consider here only one. In what follows, the following considerations might help for a better understanding for the weak coupling formulation below. Note that for a constant it follows ( 1 2 + K)1 = 0 on Γ b . Furthermore, if K ′ is the adjoint operator of K and it holds V −1 K = K ′ V −1 . Then, with (18) we see Note that this φ together with the representation formula leads to u b in Ω b , see also [McL00,Theorem 7.5]. Therefore, we introduce the following subspace . Furthermore, similar as in Section 2 we introduce a product space with its norm, namely Remark 4.1. Instead of considering a subspace and thus eliminating the constants from the solution space, a suitable orthogonal decomposition of H − 1 2 (Γ b ) in the following proofs could also be considered, see [OS14].
In this case no stabilization is needed, since both subproblems involve a Dirichlet boundary condition. Hence, we prove directly the strong monotonicity of b(·, ·). Equivalently to (11), the form b(·, ·) induces a non-linear operator B : The next theorem states the strong monotonicity of the method for the extended BVP. It can be considered as an extension to our problem setting of the stability estimate result given in [OS14] for an interior Dirichlet BVP of a diffusion equation with a hard inclusion. The key idea therein is to estimate the energy of the bounded finite element domains with the energy of some related problem in the exterior domain. If both corresponding Steklov-Poincaré operators are H 1 2 (Γ)-elliptic, then it holds with λ > 0 the minimal eigenvalue of the related exterior problem that where S ext and S int are the Steklov-Poincaré operators of the exterior and the interior domain, respectively, c.f. [OS14].

Theorem 4.3. Let us consider the non-linear operator
. Furthermore, λ 1 , λ 2 > 0 are the eigenvalues of (21) with respect to the domains Ω 1 and Ω 2 . Then the following assertions hold.
• B is Lipschitz continuous, i.e., there exists C Lip > 0 such that for all u := (u 1 , u 2 , φ) ∈ H 0 , v := (v 1 , v 2 , ψ) ∈ H 0 with C stab = min 1, • if C Ui ell > 1 4λi for i = 1, 2, then B is strongly monotone, i.e., there exists C ell > 0 such that Proof. The Lipschitz continuity follows merely from the Lipschitz continuity of U 1 and U 2 and the continuity of the boundary integral operators. The stability estimate follows strongly the steps of the proofs of [OS14, Theorem 2.2.ii.] and in [OS14, Section 5.1]. Since we are dealing with a different BVP and non-linear material tensors, we sketch the main steps of the proof, for convenience. For ease of notation, let w : First, we start with the domain parts. Provided U i , i = 1, 2, are strongly monotone, then it holds , we now consider the splitting w i = w i + w 0,i , where w i is the harmonic extension of w i|Γ i and w 0,i ∈ H 1 0 (Ω i , ∂Ω i ) as in [EOS17], for instance. From this follows where S i , i = 1, 2, denote the interior Steklov-Poincaré operators of the bounded domains Ω 1 and Ω 2 , respectively. Hence, Next, by using the contractivity of K, as given in [OS14, Lemma 2.1] (we consider here the worst case C K = 1), as well as the invertibility of V, we obtain ξ, are the Steklov-Poincaré operators associated to the corresponding exterior eigenvalue problem, see [OS14, Section 2.2]. Similarly, we assume the following spectral equivalence where λ i , i = 1, 2 are characterized as minimal eigenvalues of the related problem. Thus, Inserting (26) and (27) in (25), ξ, Vξ Γ b = ξ 2 V and some manipulations as in the proof of [EOS17,Theorem1] lead to the assertion. To prove the last claim we consider v := (v 1 , v 2 , ψ) ∈ H 0 . Note that v 1 = 0 on Γ 0,1 and v 2 = 0 on Γ 0,2 with |Γ 0,1 |, |Γ 0,2 | > 0. Due to Friedrich's inequality and (5) it follows that ∇v 1 2 is an equivalent norm on H 0 . Thus, (24) follows directly from (23).
Problem 4.4. Find u ℓ := (u 1,ℓ , u 2,ℓ , φ ℓ ) ∈ H 0,ℓ : Analogously to the interface problem, we state in the following theorem the quasi-optimality in the sense of the Céa-type Lemma of the Galerkin discretization of Problem 4.2, as well as an a priori error estimate for the introduced B-Spline discretization. To simplify the presentation, we introduce in this section a piecewise defined product space for s ≥ 0, which is used to get convergence rates with the aid of Lemma 3.11. The corresponding norm defined in the sense of (15) is denoted by · H s pw .
Theorem 4.5. For i = 1, 2, let C Ui ell > 1 4λi , where λ i > 0 are the eigenvalues of (21) with respect to the domains Ω i . Moreover, let u ∈ H 0 be the solution of Problem 4.2 and u ℓ ∈ H 0,ℓ be the discrete solution of Problem 4.4. Then the following results hold: • Quasi-optimality: Quasi-optimality follows from the strong monotonicity and Lipschitz continuity stated in Theorem 4.3, by following the lines of Theorem 3.2. The a priori estimate follows from the quasi-optimality and Lemma 3.11, as is done in Theorem 3.12 for the interface problem.
The non-linear operators U i , i = 1, 2, are now considered to have the form U i ∇u := g i (|∇u|)∇u with non-linear functions g i : R → R. Similarly to the interface problem, we state the following stability result.
Lemma 4.6. Let C Ui ell > 1 4λi , i = 1, 2, with λ i as in Theorem 4.5. Furthermore, the non-linear operators U i , i = 1, 2, shall have the form U i ∇u := g i (|∇u|)∇u with the non-linear functions g i : R → R. Moreover, let u ∈ H 0 be the unique solution of Problem 4.2 and , with i = 1, 2, be some suitable inputs. There exists C > 0 such that Proof. We know from the strong monotonicity of B that Without loss of generality, we choose v = (0, 0, 0) and note that U i ∇v i = 0, i = 1, 2, for our specific non-linearity. Since u := (u 1 , u 2 , φ) is the unique solution of the problem, we conclude that Using inequality (6) along with the boundedness of K and V, and rearranging the terms yields the assertion.
In many practical applications, one is not directly interested in the solution (u 1 , u 2 , φ) of Problem 4.2 rather than in some derived quantities. These quantities are, for example, evaluated in the exterior/air gap domain. As it can be observed for standalone BEM applications, estimating the error in functionals of the solution may lead to a so called super-convergence, i.e., linear functionals of the solution may converge better than the solution in the energy norm, see [SS10, Section 4.2.5]. With enough regularity the convergence rate doubles. In the following, this behavior is also showed for the coupled problem. For this, we use the following Aubin-Nitsche argument, similarly to [SS10, Theorem 4.2.14].
Theorem 4.7. Let F ∈ H ′ 0 be a continuous and linear functional on the solution u : of Problem 4.2 and u ℓ := (u 1,ℓ , u 2,ℓ , φ ℓ ) ∈ H 0,ℓ is the discrete solution of Problem 4.4. Furthermore, let w ∈ H 0 be the unique solution of the dual problem for all v ∈ H 0 . Then there exists a constant C 1 > 0 such that for arbitrary v ℓ ∈ H 0,ℓ , z ℓ ∈ H 0,ℓ . Furthermore, let 1 2 ≤ s, t ≤ p and remember the product space defined in (28). Provided u and w are additionally in H s pw and H t pw , respectively, there exists a constant C 2 > 0 such that Proof. The proof follows strongly the lines in [SS10, Theorem 4.2.14]. Since we allow non-linearities, we give a brief sketch. First of all, we note that Theorem 4.3 holds for arbitrary functions. Thus, wellposedness, and hence the existence of a unique solution can be established also for the dual problem (30). Furthermore, the dual problem (30), the Galerkin orthogonality b(u − u ℓ , z ℓ ) = 0 for all z ℓ ∈ H 0,ℓ , and the Lipschitz continuity of the form b(·, ·) yield to for arbitrary z ℓ = (z 1,ℓ , z 2,ℓ , ϕ ℓ ) ∈ H 0,ℓ .
Remark 4.8. In practice the functional of Theorem 4.7 may be, e.g., the representation formula of the BEM part Ω b , i.e., for u = (u 1 , u 2 , φ) ∈ H 0 there holds Next, let us assume the regularity u ∈ H p pw of the solution 4.2 and w ∈ H p pw of its dual problem (30), where the spaces are defined in (28). Then with the discrete solution u ℓ ∈ H 0,ℓ and (32) we calculate the pointwise error in Ω b as which is the maximal possible super-convergence. Since the constant C depends on u H p pw and w H p pw , a possible estimate of these norms would probably involve their right-hand sides. The right-hand side of the dual problem (30) is the functional F(u). Thus, the constant C might include a factor like . Note that this term is finite for all x ∈ R 2 \Γ b and p ≥ 0. However, because of the singularity of the kernels, its tends to infinity when approaching the boundaries. Thus, also C from (33) might tend to infinity. This effect is even more severe, if we consider functionals that involve derivatives of the kernels, e.g., for the computation of forces and torques using the Maxwell Stress Tensor. Finally, we mention that the regularity assumptions might only hold for smooth surfaces.

Numerical illustration
To illustrate the theoretical results, we consider for each model problem one example. The description of NURBS geometric entities are obtained by means of the NURBS toolbox included in GeoPDEs, which is implemented in MATLAB, see [dFRV11]. In the same spirit, the required matrices associated to the boundary integral operators are implemented by using, adapting, and supplementing some structures of GeoPDEs. The implementation of the BIOs for arbitrary ansatz functions is performed numerically using standard Gauss-Legendre quadrature for regular contributions and by means of some Duffy-type transformations with a subsequent combination of logarithmic and Gaussian quadrature for the singular parts, see, e.g., [Ban15,Chapter 4.3]. In the following, the H and H 0 -norm of (8) and (19), respectively, are computed by Gaussian quadrature. However, we replace the non-computable norm · by the equivalent norm · 2 V stated in (5). Moreover, we measure the error for the evaluated solution in the BEM-domain in the following way. First we define an evaluation path Γ e in the BEM domain. For a certain number N of evaluations points x i ∈ Γ e , i = 1, . . . , N , x i = x j , with i = j, we define the pointwise error as Here, u e ℓ and u b,ℓ are the discrete evaluations of the corresponding representation formula (2) with the Cauchy data from the corresponding discrete coupling problem. Note that for both problem types the trace has to be calculated with the aid of the jump condition (7c) and (16c), respectively.
In all our experiments, we consider uniform h-refinement, for different degrees of B-Splines, starting from the minimal degrees needed to represent the geometry exactly. Increasing the degree of basis functions is called p-refinement. Furthermore, note that the number of elements in every h-refinement 5.1. Single domain. In the first example, we consider a square domain Ω := (−0.25, 0.25) 2 and denote its boundary by Γ. We parametrize Ω as a single patch domain using linear B-Spline functions in each parametric direction. It is obvious that Assumption (A5) about the multipatch geometry is satisfied. Moreover, we consider the interface problem (7) with a linear material tensor U := Id. As in [EOS17], we prescribe the exact solutions and u e (x) = log( We calculate the jumps u 0 , φ 0 , and the right-hand side f appropriately. Solving the coupled problem using the isogeometric framework, as described in the previous section, yields a discrete solution (u ℓ , φ ℓ ) ∈ H ℓ := S 0 (Ω) × S 2 (Γ). An isogeometric approach for this example is not mandatory since the domain Ω is standard Cartesian, see, e.g., [EOS17]. However, we want to demonstrate our higher order coupling approach and in particular the super-convergence behaviour of this example. Figure 2 shows the solution u ℓ ∈ S 0 (Ω) in the interior domain, as well as the exterior solution u e ℓ in a subset of Ω e , which we call an evaluation domain Ω e e := − 1 2 , 1 2 \Ω. The exterior solution u e ℓ is obtained from the representation formula (2) (κ = 1) with the computed Cauchy data (u ℓ|Γ − u 0 , φ) from our discrete solution of the interface problem. Thereby, the degree of the considered B-Spline space for the domain discretization is p = 2 and its dimension corresponds to an h-refinement level ℓ = 20. As a first numerical experiment, we analyze the convergence of the isogeometric FEM-BEM coupling with respect to the V , which is equivalent to H-norm in Ω. Since the solution is smooth, the expected order of convergence is equal to the degree of the considered discrete space H ℓ , as given in the a priori estimate from Theorem 3.12. In Figure 3 we observe the predicted optimal convergence of the method for B-Spline spaces of degree p = 1, 2, 3.  In the second experiment, we investigate the convergence of the solution in the exterior domain. Note that our exterior solution is smooth. At a first step, we evaluate the solution on an evaluation path Γ e , which we define here as the boundary of (−0.35, 0.35) 2 . We calculate the error according to (34) with N = 20 evaluations points. In Figure 4, we observe a doubling of the convergence rates with respect to the pointwise error, which confirms the theoretical considerations in Remark 4.8, see also Remark 4.9. Furthermore, we want to investigate the dependency of the super-convergence on the position of the evaluation point for a fixed degree p = 3 of the B-Spline space. For this, we compare the convergence behavior of the exterior solution on three distinct evaluation paths. We denote the paths by Γ e,1 , Γ e,2 , and Γ e,3 , which are the boundaries of (−1, 1) 2 , (−0.35, 0.35) 2 , and (−0.26, 0.26) 2 , respectively. For each evaluation path we choose again 20 evaluation points to compute the pointwise error (34). The result is visualized in Figure 5, where we observe the expected behavior, see Remark 4.8. In particular, super-convergence is readily observed for the solution on Γ e,1 and Γ e,2 . We note that for the error in Γ e,1 we are already at machine precision. However, the related constant is larger for the solution on Γ e,2 , since the path is closer to the interface boundary Γ = ∂Ω with Ω := (−0.25, 0.25) 2 . The same behavior can be observed for the path Γ e,3 , which is even closer to Γ. However, the quality of the computation is also deteriorated in the asymptotic. Additionally, we observe saturation effects for higher refinement levels. This can be improved by increasing the number of Gaussian quadrature points N Gauss on each boundary element, as it is shown in Figure 6. However, this in turn is time consuming. With using special extraction techniques, such as the ones developed for 3-D in [SW99], this undesirable effect can be reduced. However, a further investigation is beyond the scope of this work.  ℓ (x i )| are on Γ e,1 , Γ e,2 , and Γ e,3 , which are the boundaries of (−1, 1) 2 , (−0.35, 0.35) 2 , and (−0.26, 0.26) 2 , respectively. We observe the growing constant of the super-convergence constant, which leads to an undesirable saturation for the closest path Γ e,3 with respect to Γ.  5.2. Multiple domains. In this second example, we consider the non-symmetric isogeometric FEM-BEM coupling for the extended boundary value problem (16) as described in Section 4. The topology of the model problem and the notation can be adopted from Figure 1. However, we consider here a problem domain constructed over circles, see Figure 7. In particular, if we denote by B((x 1 , x 2 ); r) a circular domain with midpoint (x 1 , x 2 ) and radius r we arrive at the following setting: Ω 1 = B((0, 0); 0.39)\B((0, 0); 0.1), Ω 2 = B((0, 0); 0.6)\B((0, 0); 0.4), and the thin air gap Ω b = B((0, 0); 0.4)\B((0, 0); 0.39), which describes in fact three rings. We prescribe the right-hand side f i in Ω i as f 1 (x 1 , x 2 ) = 0 and f 2 (x 1 , x 2 ) = 100 sin(ϕ), where ϕ is the standard angle in a polar coordinate system. The non-linear material tensor is chosen as where we choose ǫ > 0 arbitrarily such that g(t) < 1, for all 0 < t ≤ t c , and α, β such that g(t) is continuously differentiable for all t > 0 1 . In addition, we do not allow jumps, i.e., u 0,i = 0 and φ 0,i = 0.
Following the isogeometric approach, we model both domains separately and according to Definition 3.8 as multipatch domains consisting of four patches, see Figure 7. Each patch is represented exactly by a NURBS of degree p = 2 in each parametric direction. Moreover, the Assumption (A5) is obviously satisfied. Note that this model configuration with the circular geometry can be interpreted as a 2D section of a simplified 2-pole synchronous machine [Kur98, Section 5.2]. This type of applications motivates also the consideration of non-linear operators. In fact, these devices are mainly made of ferromagnetic materials, which are known to be non-linear. In particular, by neglecting anisotropies and hysteresis effects, ferromagnetic materials can be modeled by using non-linear operators of the same type as the ones we considered in Lemma 2.7 and Lemma 4.6, and for this example in (35). For more details about this topic, see [Pec04] and [Röm15], for instance. Furthermore, we refer to [BCSDG17] for electrical engineering simulations of electric machines. In this experiment, the arising non-linear problem is solved by using a standard Picard iteration method. For the stopping criterion, we consider a relative residual error of 10 −10 . In our simulation below we need an average of 35 Picard iterations to fulfill the criterion. The solutions u 1 and u 2 in the interior domains Ω 1 and Ω 2 , respectively, are visualized in Figure 8. In the context of electric machines, u i , i = 1, 2, can be interpreted as the third component of the magnetic vector potential. Note that the equipotential lines, i.e., the continuous black lines in Figure 8 are the magnetic field lines. The interaction of the magnetic fields stemming from the rotor and the stator in the air gap may induce a mechanical torque. This leads the rotor, i.e., the interior ring to move in order to reduce the (spatial) phase shift between both magnetic fields. In particular, the computation of torques and forces involves the computation of the magnetic flux density, which in turn requires the evaluation of the solution and its derivatives in the gap domain, c.f., [KFK + 97], for instance. Therefore, the super-convergence behavior in the air gap of the machine is of particular interest. In a post-processing step, we compute the magnetic reluctivity, which is defined as the reciprocal of the magnetic permeability. Formally, it is given by the function g(|∇u i |), i = 1, 2, in (35), which we evaluate by using the solutions u 1 and u 2 . Since we are considering non-linear materials, the reluctivity is not constant across the electric machine. This is depicted in Figure 9. Note that the thick black lines are the same as those in Figure 8, i.e., they represent the equipotential lines of the solution u 1 and u 2 . To verify Remark 4.8 numerically, we evaluate the solution in the BEM domain Ω b on the evaluation path given as the parametrized circle ∂B((0, 0); 0.395). Note that in this example, the BEM is applied in an interior domain Ω b . Hence, we use for the evaluation the representation formula (2) with κ = 0 and the complete Cauchy data on Γ 1 and Γ 2 , which are available after solving Problem 4.4 with the jump condition (16c). An analytical solution for our model problem is not known. Hence, to verify Figure 9. Saturation effects caused by the non-linear material tensors. The color-bar and the thin lines represent the levels of the magnetic reluctivity, which is given by g(∇u i ), for i = 1, 2, and evaluated using the derivatives of the solutions u 1 and u 2 . The thick equipotential lines show the flow direction of the magnetic field. the convergence order, we follow a standard procedure: The mesh of the current solution is successively refined three times and we calculate the corresponding discrete solutions. We apply the Aitkin's ∆ 2extrapolation to this sequence of discrete solutions and this extrapolated value is the reference solution u e (x i ) for the error = max i=1,...,N |u e (x i ) − u e ℓ (x i )| calculated from N = 20 evaluations points. This error is visualized in Figure 10 for ansatz spaces of degree p = 2 and p = 3, where we observe an amelioration of the convergence rates. Note that this amelioration depends on the quality of the numerical integration, as shown in Figure 6 for the example of Section 5.1. For this example, noticeable amelioration of the convergence rates were only observable for a high number of Gaussian quadrature points. In this case, N Gauss = 400 points were considered for the assembling of the BEM matrices, which is very time consuming. The dominance of the quadrature error for this type of evaluations can however be tackled, as mentioned in the previous section, by using special extraction techniques. Moreover, efficient assembly of the BEM matrices based on B-spline tailored quadrature rules, as given in [ACD + 18], together with suitable compression methods, see e.g., [DKSW19], would accelerate the computation considerably. However, this investigation is beyond the scope of this work.

Conclusions
The non-symmetric FEM-BEM coupling in the isogeometric context for simulating practical problems with complex geometries turns out to be a promising alternative to classical approaches. A transformation to an integral formulation allows a problem in a domain to be reduced to its boundary, where the BEM can The error = max i=1,...,N |u e (x i ) − u e ℓ (x i )| is calculated with N = 20 evaluations points on Γ e . As a replacement for the unknown analytical solution we use an Aitken ∆ 2 extrapolation of a sequence of successively refined discrete solutions. be applied. For exterior problems there is no need to truncate the unbounded domain, and for simulating thin gaps there is no need for a complicated remeshing. In both cases numerical errors can be avoided. Thanks to the definition of B-Splines, h-and p-refinements are applied in a straightforward manner. Furthermore, multiple domain modeling can be done independently. This is particularly advantageous if we consider moving or deforming geometries. A classical transmission and a multiple domain problem with parts of non-linear material are considered. Obviously, FEM is applied to the non-linear areas, whereas BEM is exclusively used for the linear problem. For both model problems, well-posedness for the continuous and discrete problem, and quasi-optimality and convergence rates for the numerical approximation are mathematically analyzed in the isogeometric framework. Furthermore, we show an improvement of the convergence behavior, if we consider the error in functionals of the solution. This is motivated by a practice-oriented application such as electric machines. Here the computation of torques are a central task and involve the evaluation of some derivatives of the solution in the BEM domain. We observe for both model applications this super-convergence, which confirms the theory. Future extensions of the method may include the consideration of parabolic-elliptic problems and a rigorous analysis of the coupling for curl curl-type equations in 3D.