Discontinuous Galerkin Isogeometric Analysis for segmentations generating overlapping regions

In the Isogeometric Analysis (IGA) framework, the computational domain has very often a multipatch representation. The multipatch domain can be obtained by a volume segmentation of a boundary represented domain, e.g., provided by a Computer Aided Design (CAD) model. Typically, small gap and overlapping regions can appear at the patch interfaces of such multipatch representations. In the current work we consider multipatch representations having only small overlapping regions between the patches. We develop a Discontinuous Galerkin (DG)- IGA method which can be immediately applied to these representations. Our method appropriately connects the fluxes of the one face of the overlapping region with the flux of the opposite face. We provide a theoretical justification of our approach by splitting the whole error into two components: the first is related to the incorrect representation of the patches (consistency error) and the second to the approximation properties of the IGA space. We show bounds for both components of the error. We verify the theoretical error estimates in a series of numerical examples.


Introduction
Isogeometric Analysis (IGA) has been introduced in [1] as a new methodology for solving numerically Partial Differential Equations (PDE). The key idea of the IGA concept is to use the superior finite dimensional spaces, which are used in Computer Aided Design (CAD), e.g. B-splines, NURBS, for both the exact representation of the computational domain and discretizing the PDE problem. Since this work, many applications of the IGA methodology in several fields have been discussed in several papers, see, e.g. the monograph [2] and the references therein, as well as the survey paper [3]. From a computational point of view, we can say that the numerical algorithm for constructing the B-spline (or NURBS) basis functions is quite simple. This helps to produce high-order approximate solutions. From the theoretical point of view, the fundamental approximation properties of the B-spline spaces on a reference domain are discussed in [4]. The approximation properties of the mapped B-spline (or NURBS) spaces are discussed in several papers, see e.g. [3,[5][6][7].
Let us consider a complex domain where its boundary is prescribed by CAD models. The CAD models can not be directly used in IGA in order to discretize the PDE problems. We need to create volumetric patch parametrizations from the CAD models. The boundary represented domain is first segmented into a collection of suitable blocks and consequently, a parametrization procedure is applied to each block. This produces the volumetric multipatch representation N i=1 i of suitable for IGA. Several segmentation algorithms and associated parametrization procedures have been discussed in the literature, see, e.g. [8][9][10][11][12]. Furthermore, we refer to [13][14][15][16] for different approaches for constructing IGA planar parametrizations without utilizing segmentation algorithms, and to [17,18] for constructing parametrizations using Bézier triangular meshes. We mention the segmentation approach presented in [19,20], from which, we have been motivated to present the current work. The main idea is to split the given boundary represented domain, using a spline curve (or face in 3d case) with the following properties: (i) must have the end points on the boundary and the tangents to be specified, (ii) the curve is reasonably regular and does not intersect the boundary of the domain, (iii) the curve cuts the domain into new subdomains with good shapes. Consequently, tensor product B-spline spaces are fitted in the collection of the subdomains for defining the tensor product B-spline surfaces or volumes [14]. Note that the previous consideration also concerns CAD models that are connected along a non-matching interface. It is important to obtain a curve that splits into new simple domains with good shapes being suitable for IGA. During the computation of the multipatch representation, errors can occur when defining the corresponding control points, see [10,14,19]. A consequence of this is non-conforming parametrizations of the patches in the sense that the images of the patch interfaces under the parametrizations are not identical. This, in turn, leads to the existence of gap and/or overlapping regions between the adjoining patches, see a schematic illustration in Figure 1(b).
This paper considers the case where there are only overlapping regions between the patches. If we apply an IGA methodology to this multipatch representation, a direct consequence is that the whole discretization error will include two (main) parts: the first naturally comes from the approximation properties of the B-spline spaces (for the purposes of this work we use B-spline spaces) and, the second comes from the geometric error. The later is due to the incorrect parametrization of the patch interfaces. Furthermore, the geometric error can be characterized as a consistency error, which consists of two error components. The first error component is related to the approximation of the jumps of the flux of the solution on the non-matching interfaces. The second component is related to the existence of more than one numerical solution in the overlapping regions.
The contribution of this paper is to develop a DG-IGA method which can be applied on volumetric patch representations with non-matching interface parametrizations. We present our methodology for discretizing the following elliptic Dirichlet boundary value problem − div(ρ∇u) = f in and u := u D = 0 on ∂ , where the diffusion coefficient ρ(x) can be discontinuous across a smooth internal interface. We derive bounds for the two main parts of the whole error. In our analysis, we derive separate bounds for the two components of the geometric error. To the best of our knowledge, we believe that this is a new area of analysis to be investigated. Our current work is the first step in the analysis, where we are developing our methodology for the numerical solution of the simple stationary diffusion problem (1). Our intention for future works is to extend the current methodology to more complicated time dependent problems, where the interface can move with time, cf. [21]. Due to the non-matching interior patch interfaces, a direct application of the classical DG numerical fluxes proposed in the literature, see e.g. [6,22], is not possible, as these fluxes are only applicable for matching interface parametrizations. In our recent papers, [23,24], we developed DG-IGA schemes for multipatch unions that include only gap regions. In particular, we considered the PDE model given in (1) and we denoted by d g the maximum distance between the diametrically opposite points located on the gap boundary. We applied Taylor expansions using the diametrically opposite points of the gap, in order to give estimates for the jumps of the solution with respect to d g . Finally, we used the same Taylor expansions in the DG-IGA scheme for constructing suitable DG numerical fluxes across the gap boundary that help on the weakly coupling of the local patch-wise discrete problems. We developed a discretization error analysis and showed a priori estimates in the DG-norm, expressed in terms of the mesh size and the gap width, i.e. O(h r ) + O(d g ), where r depends on the B-spline degree p and the regularity of the solution. In [23,24], we have shown that, if d g = O(h p+1/2 ), the proposed DG-IGA scheme has optimal approximation properties.
In this paper, we extend the previous work to multipatch unions with overlapping regions. In the analysis presented in [23,24], the whole geometric error does not include the component coming from the coexistence of different IGA solutions in the overlapping regions. Here, the new approach is to introduce local (patch-wise) auxiliary variational problems, which are compatible with the overlapping nature of the multipatch representation of . We denote the solutions of the new variational problems by u * . These problems are not consistent, in the sense that the original solution u of (1) does not satisfy them. Following the IGA concept, the B-spline spaces used for the parametrization of the patches are also used for discretizing the local auxiliary problems. We denote by u * h the produced IGA solutions. Under some regularity assumptions on u * , we can expect (see Section 3) that the IGA solution u * h has optimal approximation properties associated with u * . However, we can not directly infer that u * h can approximate in an optimal way the solution u of the original problem. In our analysis, we provide an estimate for the consistency error u − u * and consequently using the triangle inequality u − u * h DG ≤ u * − u * h DG + u − u * DG , we can derive an estimate for the error between the exact solution u and the IGA solution u * h . The mesh-dependent norm · DG is defined in Section 2. We give error estimates for both terms u * − u * h DG and u − u * DG expressed in terms of the mesh size h and the quantity d o , which is introduced in our analysis in order to quantify the width of the overlapping regions. In particular, we show that under appropriate assumptions on the data and for the case where d o is of order h λ , λ ≥ p + 1 2 , the proposed DG-IGA scheme has optimal convergence properties. This convergence result is similar to the result in [23,24].
In another work, which is under preparation, we apply the same approach to solve problems on multipatch partitions, which can include gap and overlapping regions. We present numerical solutions in multipatch unions with more complicated gaps and overlapping regions. We also provide details related to the implementation of the proposed DG-IGA scheme. In the same work, we also discuss issues related to the construction of domain-decomposition methods on these multipatch representations and provide several numerical tests for evaluating their performance. The first results in this direction can be found in [25].
We note that IGA multipatch representations with non-matching interfaces meshes, overlapping regions and trimmed patches have been considered in many publications. For the communication of the discrete patch-wise problems, several Nitsche's type coupling methods involving normal flux terms have been applied across the interfaces, see e.g. [22,[26][27][28] and references therein. We mention also that in [29], DG-IGA methods have been presented to discretize Laplace problems on multipatch unions with large overlapping regions. The proposed strategy follows the additive Schwartz methodology. To the knowledge of the authors, there are no works that analytically discuss estimates for the error, which is caused by the incorrect representation of the shape of the patches. The purpose of this work is to present such an error analysis.
The structure of the paper is as follows: Section 2 presents the PDE model, briefly reviews the Bspline spaces and describes the case of having non-matching parametrized interfaces with overlapping regions. Section 3, presents in detail the perturbation problems, the bounds for the consistency error, the proposed DG-IGA scheme and the error analysis. Section 4, includes several numerical examples that confirm the theoretical estimates. The paper closes with the Conclusions.

The elliptic diffusion problem
The weak formulation of the boundary value problem (1) reads as follows: for given source function f ∈ L 2 ( ) find a function u ∈ H 1 0 ( ) such that the variational identity is satisfied, where the bilinear form a(·, ·) and the linear form l f (·) are defined by respectively. The given diffusion coefficient ρ ∈ L ∞ ( ) is assumed to be uniformly positive and piece-wise (patch-wise, see below) constant. These assumptions ensure the existence and uniqueness of the solution due to Lax-Milgram's lemma. For simplicity, we only consider pure Dirichlet boundary conditions on ∂ . However, the analysis presented in our paper can easily be generalized to other constellations of boundary conditions that ensure existence and uniqueness such as Robin or mixed boundary conditions. In what follows, positive constants c and C appearing in inequalities are generic constants that do not depend on the mesh-size h. In many cases, we will indicate on what may the constants depend on.

B-spline spaces
In this section, we briefly present the B-spline spaces and the form of the B-spline parametrizations for the physical subdomains. For a better presentation of the B-spline spaces, we start our discussion for the one-dimensional case. Then we proceed to higher dimensions. We refer to [2,4,31] for a more detailed presentation. Consider The B-spline basis functions are defined by the Cox-de Boor formula, see, e.g. [2,31], We assume that m j ≤ p for all internal knots, which in turn gives that, at z j the B-spline basis functions have κ j = p − m j continuous derivatives. Let us now consider the unit cube = (0, 1) d ⊂ R d , which we will refer to as the parametric domain. Let the integers p and n k denote the given B-spline degree and the number of basis functions of the B-spline space that will be constructed in x k -direction with k = 1, . . . , d. We introduce the d-dimensional vector of knots = ( 1 , . . . , k , . . . , d ), with the particular components given by Given the knot vector k in every direction k = 1, . . . , d and using (7), we construct the associated univariate B-spline basis functions, {B 1,k (x k ), . . . ,B n k ,k (x k )} of the spaceB k ,p , see, e.g. [31] for more details. Accordingly, the tensor product B-spline space is defined where eachB j (x) has the form In the IGA framework, the computational domain is described as the image of under a B-spline parametrization mapping of the form where C j , j = 1, . . . , n are the control points andx = −1 (x), see Figure 1(a). Following the IGA methodology, [1,2], the B-spline spaces for discretizing the PDE problem are defined by using the mapping given in (10), i. e. we define the B-spline space in by

B-spline spaces on multipatch representations
Let us suppose that the domain can be described as a union of N-subdomains with interior interfaces F ij = ∂ i ∩ ∂ j , for 1 ≤ i = j ≤ N. We further suppose that every subdomain i has its own parametrization i , which is defined by the corresponding B-spline spaceB i ,p and the corresponding control points C j , see (10). Here i denotes the knot vector related to i . An illustration for N = 2 is given in Figure 1(a). The subdomains i are referred to as patches. In an analogous way as in (11), we define the physical patch-wise B-spline spaces B i ,p for i = 1, . . . , N. We define the global discontinuous B-spline space V B with components on every B i ,p Assumption 2.1: Assume that every i , i = 1, . . . , N is sufficiently smooth and there exist constants in , whereÊ m are the micro-elements and h i is the mesh size, which is defined as follows. Given an elementÊ m ∈ T , whose vertices are the images of the vertices of the corresponding parametric mesh T

Multipatch representation of the computational domain
In many practical applications, the parametrization of a boundary represented domain by a single B-spline patch may not be possible. In order to discretize a PDE problem following the IGA framework in this situation, we represent the domain as a multipatch. Following the methodology presented in [9,19], the initial domain is firstly segmented into a collection of simple subdomains, e.g. topological hexahedra. Consequently, a suitable parametrization mapping is constructed for each subdomain for obtaining the multipatch representation of . The final parametrization mappings of the adjoining patches must provide identical images for the common interfaces. In particular, for a DG-IGA discretization of the model (1), it would be preferable to produce a multipatch partition of compatible with the variations of the coefficient ρ, i.e. the patches to be coincided with the parts of where the coefficient ρ is constant. For example, let us consider Figure 1(a). In this case, the domain is described as a union of two non overlapping patches, see (12), i.e.
where the interface F 12 coincides with the physical interface. We use the notation T H ( ) := { 1 , 2 } for the union (14). For each i , i = 1, 2, there exists a matching parametrization mapping such that The control points, which are related to the patch interface F 12 , are appropriately matched in order for the parametrizations 1 and 2 of the neighboring patches to give the same image for the parametrized interface F 12 . Based on T H ( ), we can independently discretize the problem on the different patches i , i = 1, 2, using interface conditions across F 12 for coupling the local problems. Typically, the interface conditions across F 12 concern continuity requirements of the solution u of (1), i.e. (15) where n F 12 is the unit normal vector on F 12 with direction towards 2 , and ρ i , u i , i = 1, 2 denote the restrictions of ρ and u to i correspondingly. The conditions (15) can be ensured by considering appropriate regularity assumptions on the solution u. We note that these types of multipatch representations have been considered in [6] and DG-IGA methods have been proposed for discretizing the problem (1).
Anyway, for simplicity, we develop our analysis based on Figure 1. We introduce the appropriate spaces. Let ≥ 2 be an integer, we define the broken Sobolev space

Assumption 2.3:
We assume that the solution u of (4) belongs to

Non-matching parametrized interfaces
Typically, the segmentation procedure will generate multipatch representations that have possibly non-matching interface parametrizations, [10]. The result is the existence of gap and overlapping regions in the multipatch representation of the domain . In [23,24], we developed DG-IGA schemes for multipatch unions that only include gap regions. In this work, we focus on multipatch representations with small overlapping regions, see Figures 1(b) and 2(a,b). Due to the non-matching parametrization of the interior patch interfaces, a direct application of the interface conditions (15) for deriving DG-IGA methods, is not possible. The purpose of this paper is to investigate the construction of auxiliary interface conditions on the boundary of the overlapping regions; which can be used for constructing DG-IGA schemes. We present a discretization error analysis separating the whole discretization error into two parts: the first naturally comes from the approximation properties of the B-spline spaces and the second, is the geometric error coming from the incorrect parametrization of the patches. The geometric error is considered as a consistency error and it is further separated into two components. The first error component is related to the approximation of the auxiliary flux terms across the non-matching interfaces and the second component is related to the existence of more than one numerical solution in the overlapping regions. Remark 2.1: Alternatively, one can perform additional post-processing steps after the segmentation procedure to obtain matching interfaces. However, this procedure may increase the number of patches and the number of control points. Moreover, the newly obtained patch interfaces may not coincide with the original interface of the PDE problem, and thus the geometrical consistency error will still exist.

The overlapping regions
As we mentioned above, for the sake of simplicity, we restrict our investigation to the case where the multipatch representation of has two overlapping patches, see Figure 2. Let us suppose that where each patch has its own parametrization * 1 : → * 1 and * 2 : → * 2 , as it is shown in Figure 2(c,d).
We denote the overlapping region by o21 , i.e. o21 = * 1 ∩ * 2 ⊂ . We denote the interior boundary faces of the overlapping region by F o12 = ∂ * 1 ∩ * 2 and F o21 = ∂ * 2 ∩ * 1 , which implies that ∂ o21 = F o12 ∪ F o21 . Finally, let n F oij denote the unit exterior normal vector to F oij , for 1 ≤ i = j ≤ 2. For functions u * i defined in * i , i = 1, 2 we identify their pair (u * 1 , u * 2 ) by u * , which is equal to u * i on * i . We develop our analysis for the case where the overlapping region is, at least locally, a convex region. We introduce an assumption related to the form of the faces F o21 and F o12 . This assumption will help us to simplify the analysis, to explain our ideas in a better way, and to keep the notation to a minimum, e.g. the form of Jacobians, the form of face integrals, etc.
, and it coincides with the physical interface, i.e. F o12 = F 12 , see (14). (b) The face F o21 can be described as the set of points (x, y, z) satisfying where x Mo and y M o are real numbers, and ζ 0 (x, y) is a given smooth function, see Figure 2.
We note that we will discretize the PDE problem using the B-spline spaces defined in * 1 and * 2 , see (11). We will couple the resulting discrete problems in * 1 and in * 2 following discontinuous Galerkin techniques, this means by introducing appropriate numerical fluxes on F o12 and on F o21 . In order to construct these fluxes, we assign the points located on F o12 to the diametrically opposite points located on F o21 . Based on Assumption 2.4, we can construct a parametrization for the face F o21 , i.e. a mapping o12 : F o12 → F o21 , of the form where n F o12 is the unit normal vector on F o12 and ζ o has the same form as in (18), and it is a B-spline function with the same degree as the mapping * 2 , since the face F o21 is the image of a face of ∂ˆ under the mapping * 2 . For the schematic illustration in Figure 2(c,d), we have F o21 = * 2 (F 3 ). Utilizing the mapping o12 given in (19), we consider each point x o2 ∈ F o21 as an image of a point x o1 ∈ F o12 under the o12 , see Figure 2(a,c). Finally, we introduce a parameter d o , which quantifies the width of the overlapping region o21 , i.e.
In the present work, we are interested in overlapping regions with small size, and in particular for regions where their width d o decreases polynomially in h, i.e.
Based on this, we assume that n F o12 ≈ −n F o21 , and define the mapping o21 : where o21 is the inverse of o12 .

Remark 2.2:
By introducing Assumption 2.4, the face F o21 can be considered to be the graph of ζ 0 . This helps us to determine a parametrization for F o21 using the function ζ 0 , and through this, we can relate to every point x o2 ∈ F o21 a point x o1 ∈ F o12 . This has been achieved by the mapping o12 in (19). The mapping o12 here has a simple form and simplifies the analysis. Other mappings (even for more complex overlapping regions), for relating the points x o2 ∈ F o21 to the points x o1 ∈ F o12 can be constructed. For example, one can follow the ideas of minimum distance problem presented in [32] and can relate the points x o2 ∈ F o21 to the points x o1 ∈ F o12 by introducing the condition (x o1 − x o2 )//n F o12 . In any case, the parametrization mappings must satisfy the maximum width condition given in (21).

Remark 2.3:
As we previously said, the face F o12 is the image of a face of ∂ under the mapping * 1 , for example in Figure 2(c) we have F o12 = * 1 (F 1 ). On the other hand, the face F o21 is an interior curve for * 1 , see Figure 2(a,c). Thus, one could try to see F o21 as an image of a curveF o21 ⊂ under the mapping * 1 , i.e. F o21 = * 1 (F o21 ). In that way, it would be advantageous to have a parametric description of ∂ o21 using the mapping * 1 , which in turn would help to link the diametrically opposite points x o1 and x o2 , see (19). This approach requires the computation of the inverse ( * 1 ) −1 , which in general is very costly and demands the use of a Newton approach for solving many nonlinear systems. We are thus led to see the faces of ∂ o21 as images of both mappings * 1 and * 2 . We note also that the mappings o12 and o21 are introduced and used only for deriving the discretization error analysis. They are not used in the computation of the entries of the system matrix of the discrete DG-IGA scheme, see also discussion in Subsection 4.1.

Remark 2.4:
In Section 4, we give details of implementing the proposed method to more complicated overlapping regions. We present examples where F o12 is not an elementary face in the plane and does not coincide with the physical phase F 12 . Also, we note that for multipatch representations with large overlapping regions, one needs to work in a different direction and to use ideas coming from Schwarz domain-decomposition methods, see [33].

Remark 2.5:
The methodology can also be applied to the case where the interior faces of ∂ o21 do not touch the boundary ∂ .

The patch-wise problems and the fluxes
We compute a numerical solution in each * i , i = 1, 2 using the corresponding diffusion coefficient ρ i and the corresponding B-spline spaces defined in * i , lets us say B * i ,p . Therefore on o21 we will have the coexistence of two different numerical solutions and this makes the computation of the bounds for the error u − u * h DG more complicated. The norm · DG is defined in (24). The idea in our approach is to introduce local (patch-wise) problems a ·). Using the triangle inequality, we split the error as Then we estimate every term separately.
According to Assumption 2.3, we make the following assumption. In Appendix, see Subsection A.1, we give an estimate for the distance of the solutions u and u * .

The non-consistent terms
We multiply (25b) by φ 2,h ∈ B * 2 ,p , integrate over * 2 and apply integration by parts, then after a few calculations we find that * 2 Working in a similar way, we multiply (25a) by φ 1,h ∈ B * 1 ,p , and we have * 1 as well as By (27), (29) and (30), we get that Also for the solution u of (4) we have that From the conditions (15), the forms defined in (29), (30) and the relations (31) and (32), we derive that By a simple application of divergence theorem, we get Finally, by (33b) and (34), we deduce that Proposition 3.1: Let φ 2,h ∈ B * 2 ,p . There is a c > 0 dependent on ρ but independent of u and o21 such that Proof: Let v = (0, yφ 2 2,h ). The divergence theorem for v on o21 yields, Using that y ≤ d o and applying (2) in (37) we obtain Gathering similar terms and choosing such that 1 − 2 > 0, we get where c ρ := max{1/ρ 2 , 1/{ρ}} and we used that d 2 o ≤ d o h, see (21). Rearranging appropriately the constants in (39) yields (36).
2 ,p and let u * 2 and u be the solutions of (26d) and (4) respectively. There are constants c 1 , c ρ > 0 dependent on F o21 but independent of h such that Proof: The Cauchy-Schwartz inequality implies that Using (36) in (41), the required assertion follows easily. Inequality (40b) follows immediately from (34) and (40a).

The discrete problem
In this section, we use the bilinear forms given in (29) to define the patch-wise discrete problems. Based on Remark 3.1, and using the interface conditions on F o21 and F o12 , which are introduced in (25a) and (25b), we imply the following interface condition Next, we appropriately modify the flux terms F o21 ρ 2 ∇u * 2 · n F o21 φ 2,h dσ and F o12 ρ 1 ∇u * 1 · n F o12 φ 1,h dσ in (29) using Taylor expansions.

Taylor expansions
Let x, y ∈ * 2 and let f ∈ C m≥2 ( * 2 ). We recall Taylor's formula with integral remainder where R 2 f (y + s(x − y)) and R 2 f (x + s(y − x)) are the second order remainder terms defined by

Modifications of the fluxes on ∂ o21
To illustrate the use of (43) in our analysis, we consider the simple case of Figure 2(a). Let the points x o1 ∈ F o12 and x o2 ∈ F o21 be such that x o2 = o12 (x o1 ) as in Figure 2(a). We apply (43) using the points x o1 and x o2 and we have and setting In the same way we can get Now denoting r oij = x oi − x oj we have that n F oij = r oij /|r oij |. For keeping notation simple, we denote the Taylor's residuals as R 2 u * x o1 := R 2 u * (x o1 + s(x o2 − x o1 )) and R 2 u * x o2 := R 2 u * (x o2 + s(x o1 − x o2 )). Note that by the interface conditions (42), we have that ({ρ}/h) F o21 (u * 2 (x o2 ) − u * 1 (x o2 ))φ 2,h dσ = 0, and using also (45) we modify the fluxes in (29) as follows F o21 where {ρ} = 1 2 (ρ 1 + ρ 2 ). Similarly, we have F o12

The global modified form
We consider the global bilinear form a * (·, ·) : V * h × V * B → R, which is formed by the contributions of a * i,h (·, ·), i = 1, 2 given in (29) and the flux forms given in (46), that is

Remark 3.2:
Note that the exact solution u has similar regularity properties to the solution u * , see Assumption 2.3, and thus we can derive for u an analogous formulation as this in (47).

The DG-IGA scheme
In view of (47), we define the forms where η > 0 is a parameter that is going to be determined later. Based on the forms defined in (48), we introduce the discrete bilinear form A h (·, ·) : V * B × V * B → R and the linear form F h : V * B → R as follows Finally, the DG-IGA scheme reads as follows:

Remark 3.3:
From the relations (31), (33), the Remark 3.2 and the forms given in (49) and in (50), we can derive that Below, we quote few results that are useful for our error analysis. For the proofs we refer to [23][24][25]. (21), there exist positive constants C 1 and C 2 independent of h such that the estimates

Lemma 3.1: Under the assumption
hold for the solutions u * and u, and

Lemma 3.2:
The bilinear form A h (·, ·) in (49) is bounded and elliptic on V * B , i.e. there are positive constants C M and C m such that the estimates hold for all v h , φ h ∈ V * B provided that η is sufficiently large, see [25].

Lemma 3.3:
Let the assumption (21) and let β = λ − 1 2 . Then there is a constant C * > 0 depending on the parametrization mappings but independent of h such that the inequality Proof: Recall the definition of the pair function spaces in (23). In view of the form of A h (·, ·) and applying (2), we have Now, let us first show an estimate for the normal fluxes on F oij . Since v ∈ V * h the normal traces on the interfaces are well defined. Using again (2), we obtain where |J o12 | is the measure of the Jacobian of o12 . In the same way, we show Gathering together the above bounds, we show (55). For the case where (v, φ h ) ∈ (V + V * B ) × V * B we work similarly.

Discretization error analysis
Next, we discuss interpolation estimates that we will use to bound the discretization error. We recall the definition of the pair function spaces in (23). Let v ∈ H (T * H ( )) with ≥ 2. Under Assumptions 2.1, and using the results of [3,5], we can construct a quasi-interpolant * hold, where s = min( − 1, p) and the C 1,i , C 2,i depend on p, * i , θ but not on h.
The properties (54), (55) of A h (·, ·) and (52) imply where the bound (53) has been used previously. Setting in (61) z h = * h u * , and then using the triangle inequality c m u * h − u * DG ≤ c m u * h − * h u * DG + c m u * − * h u * DG together with the estimate in (59), we derive (60).

Main error estimate
The estimate given in (60) concerns the distance between the DG-IGA solution u * h ∈ V B and the solution u * ∈ H (T * H ( )) of the problems in (26). Below we give an estimate between the solution u of (4) and the DG-IGA solution u * h . In the proof of this result we need the following interpolation where the quasi-interpolant * is defined in (58) and s = min ( − 1, p). The proof of (62) is provided in the Appendix. We suppose further that d o = h λ , λ ≥ 1 is the width of o21 . The following error estimate holds where β = λ − 1 2 , s = min ( − 1, p), the constant C depends on the constants in (59), (55) and (54), and K o has the form given in Lemma 3.1.

Proof:
By the definition of the discrete DG-IGA scheme in (51), the properties of A h (·, ·) and the Remark 3.3 we have Setting z h = * h u into (64), using (61), (59), and (62) and gathering together the similar terms we deduce that Applying the triangle inequality the desired estimate follows.

Remark 3.4:
In the description of the problem and in the derivation of the DG-IGA scheme, we focused on using B-spline spaces. The same derivation can be applied for the case of NURBS spaces.

Implementation remarks
In this paragraph we focus on the implementation of the proposed scheme for both two and three dimensional problems. For simplicity of the presentation, we first discuss the case of having two patches. Afterwards, we explain how the same ideas can be generalized to the multipatch case. Initially, we consider interfaces with matching meshes, i.e. the number of edge elements on F o21 is the same as the number on F o12 , as shown in Figure 3.
For the computation of the numerical flux terms of the DG-IGA scheme given in (48a), a Gauss quadrature rule is applied on every edge. The first term of the numerical flux can be directly computed by using the Gauss rule and the related Jacobian term. For the computation of the jump terms, we must know the diametrically opposite edge and the associated quadrature point that are located on the other interface. We could proceed to this direction by constructing and using the mappings o21 and o12 given in (19) and (22) respectively. For the practical implementation, it would be preferable to proceed without the construction of these mappings.
We first assign the edges belonging to F o21 to the edges belonging to F o12 , for the example given in Figure 3(a), the edge e 1 2 of F o21 is assigned to e 2 2 of F o12 . In Figure 3(a) the Gauss point is denoted by x q and y q correspondingly. The edge e 1 2 is the image of the edgeê 1 2 under the parametrization * 2 , and also the edge e 2 2 is the image of the edgeê 2 2 under the parametrization * 1 . Hence, the Gauss rule is transformed back to boundary edges of the parametric domain, and for every Gauss pointŷ q there is always a corresponding Gauss point on the other associated edge to perform the numerical integration. For the configuration given in Figure 3(a), the other associated edge is located on faceF 1 and the corresponding Gauss point is denoted byx q . Thus, having defined the quadrature points on the boundary edges of , we can compute the interface terms of the numerical flux of the DG-IGA scheme.
Note that the above approach is quite simple and it follows the same ideas that we use for computing the numerical fluxes in the case of matching parametrized interfaces. It can be also applied for the case of having gap regions between the patches. The advantage of implementing this approach is that we can develop a flexible DG-IGA code which can treat patch unions with matching and non-matching interfaces in a similar way. Note also that the previous approach can be easily combined with the adaptive numerical quadrature methods presented in [34], in order to discretize the problem using non-matching structured meshes on the overlapping faces.
Overlapping regions with boundary consisting of more than two faces are shown in Figure 3(b). We consider again the case where the maximum number of the overlapping patches is two. For the example shown in Figure 3(b) the domain has four patches and the boundary of the overlapping region is compromised of the four faces F oi , i = 1, . . . , 4. Anyway, the evaluation of the interface numerical fluxes, in this case, needs more work. We first find the faces that form the boundary of the overlapping regions. Then between these faces, we determine those that are diametrically opposite, and we continue following the procedure described in the previous paragraph. This type of overlapping regions are discussed in the numerical Example 4.3.
It is clear that through a segmentation and parametrization procedure, overlapping regions with more complicated shapes than the shapes in the examples shown here can exist, e.g. more than two overlapping patches, T-joint faces on the boundary, see, e.g. [10]. In an ongoing work, we are extending the present methodology to treat these cases. We also are constructing domain-decomposition methods, [35], on these type of multipatch representations and we are discussing the influence of the size of the overlapping region on the performance of the proposed methods. We treat these cases by extending the ideas presented here. Again, we first find the interior faces that form the boundary of the overlapping regions and then we construct the numerical fluxes between the opposite faces in the same way as we presented in the previous sections. The first results of this work are included in [25]. We point out that when there are more than two overlapping patches, it is possible to have more than two (overlapped) numerical solutions on the overlapping regions. Consequently, in the error analysis, we will have more than one non-consistency terms to estimate, see Subsection 3.2. Also, we add that for multipatch unions with large overlapping regions, we can apply ideas coming from Schwarz domain-decomposition methods in order to treat the whole problem, see [29,33].
Finally, we mention that during the investigation of the proposed methodology in Section 3, we considered simple interior penalty fluxes on ∂ o21 . For the performance of the numerical examples below, we have implemented the corresponding symmetric numerical fluxes, i.e. − F o12 See [6,23]. During the derivation of the discretization error analysis in [25], estimates and bounds for the values of the parameter η are introduced. Anyway for the numerical examples here we set η = 2((p + 1)(p + d)/d) 1/2 where d is the dimension of .

Numerical examples
In this section, we perform several numerical tests with different shapes of overlapping regions as well as combinations with non-homogeneous diffusion coefficients for two-and three-dimensional problems. We investigate the order of accuracy of the DG-IGA scheme proposed in (49 . We mention that, in the test cases, we use highly smooth solutions in each patch, i.e. p + 1 ≤ , and therefore the order s in (60) and (63) becomes s = p. The predicted values of power β, the order s and the expected convergence rate r, for several values of λ, are displayed in Table 1. In any test case, the overlap regions are Table 1. The values of the expected rates r as they result from estimate (63).  solution is given by the formula The boundary conditions and the source function f are determined by (67 . For the numerical tests, we use B-splines of the degree p = 2. Hence, we expect optimal rates for λ ≥ 2.5. In Figure 5(b) the approximate solution u * h is presented on a relative coarse mesh with d o = 0.06. The results of the computed rates are presented in Figure 5(c). For all test cases, we can observe that our theoretical results presented in Table 1 are confirmed.

Example 4.3 (Overlapping regions with more than two faces):
The proposed method is now applied to a more complicated overlapping boundary with multiple faces. The geometric description of the problem in shown in Figure 6(a), the domain is decomposed into four patches and the overlapping region is defined by four interfaces. The exact solution is given by u(x, y) = sin(π(x + 0.4)) sin(2π(y + 0.3)) + x + y.
The diffusion coefficient is globally constant, i.e. ρ = 1, the right-hand side f and the Dirichlet boundary conditions are manufactured by the solution (68). We solved the problem using B-splines of degree p = 2. In Figure 6(b), we present the contours of the DG-IGA solution u * h computed on the second mesh in a sequence. The corresponding error convergence results for the four values of λ, i.e. λ ∈ {1, 2, 2.5, 3}, are given in Figure 6(c). We can observe the suboptimal behavior of the rate for λ = 1 and λ = 2 as we move to the last mesh levels. On the other hand, we have optimal rates for the rest values of λ. The numerical rates for all λ cases are in agreement with the theoretical results.

Three-dimensional numerical examples
As a final example, we consider a three-dimensional test. The domain has been constructed by a straight prolongation to the z-direction of a two-dimensional (curved) domain, see Figure 7(a). The two physical domains 1 and 2 have the physical interface F 12 consisting of all points (x, y, z) such that −1 ≤ x ≤ 0, x + y = 0 and 0 ≤ z ≤ 1, see Figure 7(a). The knot vector in the z-direction is simply 3 i = {0, 0, 0, 0.5, 1, 1, 1} with i = 1, 2. We solve the problem using matching meshes, as depicted in Figure 7(a). The B-spline parametrizations of these domains are constructed by adding a third component to the control points with the following values {0, 0.5, 1}. The completed knot vectors k=1,2,3 i=1,2 together with the associated control nets can be found in G+SMO library in the file bumper.xml. The overlap region is artificially constructed by moving only the interior control points located at the interface into the normal direction n F 12 of the related interface F 12 . Due to the fact that the overlap has to be inside the domain, we have to provide cuts through the domain in order to visualize them, cf. Figure 7 with diffusion coefficient ρ = {1, π/2}. Note that the interface conditions (15) are satisfied. The two physical subdomains, the initial matching meshes and the exact solution are illustrated in Figure 7(a).
We construct an overlap region with d o = 0.5 and solve the problem using p = 2 B-spline functions. In Figure 7(b), we show the domain meshes T (i) h i , * i , i = 1, 2, the overlapped meshes in o12 and we plot the contours of the produced solution u * h for the interior plane z = 0.5. We can see that, both faces of ∂ o12 are not parallel to the Cartesian axes. Moreover, we point out that the problem has been solved using non-matching meshes on the overlapping interfaces. We have computed the convergence rates for four different values λ ∈ {1, 2, 2.5, 3} related to the overlapping region width d o = h λ . The results of the computed rates are plotted in Figure 7(c). We observe from the plots that the rates r are in agreement with the rates predicted by the theory, see estimate (63) and Table 1.

Conclusions
In this article, we have proposed and analyzed a DG-IGA scheme for discretizing linear, secondorder, diffusion problems on IGA multipatch representations with small overlapping regions. This type of multipatch representation leads to the use of different diffusion coefficients on the overlapping patches. Auxiliary problems were introduced in every patch and DG-IGA methodology applied for discretizing these problems. The normal fluxes on the overlapped interior faces were appropriately modified using Taylor expansions, and these fluxes were further used to construct numerical fluxes in order to couple the associated discrete DG-IGA problems. The method was successfully applied to the discretization of the diffusion problem in cases with complex overlaps. A priori error estimates in the DG-norm were shown in terms of the mesh-size h and the maximum width d o of the overlapping regions. The estimates were confirmed by solving several two-and three-dimensional test problems with known exact solutions. The theoretical estimates were also confirmed by performing numerical tests using non-matching grids on the overlapping faces. Recalling u i = u| i and using the properties of the extension operator and (58) we have and ∇(u 2 − * 2,hũ 2 ) L 2 ( 2) ≤ ∇(E 2 u − * 2,hũ 2 ) L 2 ( * 2 ) ≤ C intp C E2 h s u 2 2 , ( A 8 b ) where s = min(p, − 1). Using the trace inequality, [6], v 2 ) and proceeding as in (A8) we can show h ∇(u 1 − * 1,hũ 1 ) 2 Gathering the inequalities (A6), (A8) and (A9b) we can derive (62).