Discontinuous Galerkin Isogeometric Analysis for segmentations generating overlapping regions
Christoph Hofer Ioannis Toulopoulos


 
 a4paper,
 left=20mm,
 right=18mm,
 top=18mm,
 bottom=20mm,
 







 












airforceblue









  

  


fminipage


thmTheorem[section]
assumeAssumption
defn[thm]Definition
rem[thm]Remark
equationsection
    |ioannis.toulopoulos@ricam.oeaw.ac.at|
    |christoph.hofer@jku.at|
    
     
     


DG-IGA on multipatch representations with  overlapping regions. 

C. Hofer,   I. Toulopoulos 

 1 Johannes Kepler University (JKU),

				Altenbergerstr. 69, A-4040 Linz, Austria,

	    2 Johann Radon Institute for Computational and Applied Mathematics (RICAM),

				Austrian Academy of Sciences

				Altenbergerstr. 69, A-4040 Linz, Austria,




 

  Discontinuous Galerkin  Isogeometric Analysis for 
segmentations generating  overlapping regions.
          
    
Christoph Hofer^1
Ioannis Toulopoulos^2

    December 30, 2023
==============================================================================================================



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 <cit.> 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 to several fields have been discussed 
in several papers, see, e.g.,   the monograph   <cit.> and the references within, as well as
the survey paper <cit.>. 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 <cit.>. The approximation properties of the 
 mapped B-spline (or NURBS) spaces are discussed in several papers, 
see e.g.,  <cit.>, <cit.>, <cit.>, <cit.>.
 
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 discretise 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  ∪_i=1^N Ω_i of Ω suitable for IGA.  
Several  segmentation algorithms  and 
associated parametrization procedures have been discussed in the literature, see, e.g., <cit.>, <cit.>,<cit.>, <cit.>. 
Furthermore, we refer to  <cit.>, <cit.> and   <cit.> for different approaches for constructing   IGA planar parametrizations without 
utilizing segmentation algorithms. 
We mention the segmentation approach presented in <cit.>, and  <cit.>, 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 <cit.>. 
Note that the previous consideration is 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 <cit.>, <cit.> and <cit.>. 
A consequence of this is a 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  Fig. <ref>(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 Ω     
    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 this is an 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  (<ref>).
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.  <cit.>. 


Due to the non-matching interior patch interfaces, a direct application of the classical DG numerical fluxes
proposed in literature, see e.g.   <cit.>, <cit.>,  is not possible, as these fluxes are only applicable for matching interface parametrizations.  
In our recent papers, <cit.> and <cit.>, we  developed
DG-IGA  schemes for  multipatch unions that include  only gap regions.
In particular,  we  considered the  PDE model given  in (<ref>) 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 normal fluxes  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., 𝒪(h^r)+𝒪(d_g), where r depends on
 the B-spline degree p and  the regularity of the solution.   
  In  <cit.> and <cit.>, we have shown that,
 if d_g=𝒪(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  <cit.> and <cit.>, the whole geometric error does not include the  component coming from  the co-existence 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  (<ref>)
  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 
  <cit.> and <cit.>. 
     
   In a future work,  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 <cit.>.
   
  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., 
   <cit.>,<cit.>,<cit.>,<cit.> 
   and  references therein. We mention also that in <cit.>, 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 B-spline 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 MODEL PROBLEM



 §.§ Preliminaries

Let Ω be a bounded Lipschitz domain in ℝ^d, d=2,3, 
and let α=(α_1,…,α_d) be a multi-index of non-negative integers 
α_1,…,α_d
with degree |α| = ∑_j=1^dα_j. 
For any  α,  we define the differential operator
D^α=D_1^α_1… D_d^α_d,
with D_j = ∂ / ∂ x_j, j=1,…,d, and D^(0,…,0)ϕ=ϕ.
For a  non-negative integer m, let C^m(Ω) denote the space
of all functions ϕ:Ω→ℝ, whose 
partial derivatives
D^αϕ of all orders |α| ≤ m are continuous in Ω. 
Let  ℓ be a non-negative integer. 
 As usual, 
 L^2(Ω) denotes the Sobolev space for which 
 ∫_Ω|ϕ(x)|^2 dx < ∞, endowed with the norm
 ϕ_L^2(Ω) = (∫_Ω|ϕ(x)|^2 dx)^1/2, 
 and  L^∞(Ω) denotes the functions that are essentially bounded. Also 
 
    H^ℓ(Ω)={ϕ∈ L^2(Ω): D^αϕ∈ L^2(Ω), for all |α| ≤ℓ},

denote the standard Sobolev spaces endowed with the following  norms 
 
    ϕ_H^ℓ(Ω) = (∑_0≤ |α| ≤ℓD^αϕ_L^2(Ω)^2)^1/2.

 We identify L^2 and H^0 and also define the subspace H^1_0(Ω) and H^1_Γ(Ω) of H^1(Ω) 
 
    -0.2cm 
     	H^1_0(Ω) ={ϕ∈ H^1(Ω):ϕ=0 on ∂Ω},  
     	H^1_Γ(Ω)={ϕ∈ H^1(Ω):ϕ=0 on Γ⊂∂Ω,  |Γ| >0.}.


We recall Hölder's and Young's inequalities 

    |  ∫_Ωϕ_1 ϕ_2 dx | ≤ϕ_1_L^2(Ω)ϕ_2_L^2(Ω)    |  ∫_Ωϕ_1 ϕ_2 dx |≤ϵ/2ϕ_1^2_L^2(Ω)+  
     1/2ϵϕ_2^2_L^2(Ω),

that hold for all ϕ_1∈ L^2(Ω) and ϕ_2∈ L^2(Ω) 
and for any fixed ϵ∈ (0,∞). In addition, we  recall  trace and Poincare's inequalities, <cit.>,

    ϕ^2_L^2(∂Ω)≤    C_trϕ_L^2(Ω)ϕ_H^1( Ω),
    ϕ_L^2(Ω)≤   meas_ℝ^d(Ω)  ∇ϕ_L^2(Ω),  for  ϕ∈ H^1_Γ(Ω).




 §.§ The elliptic diffusion problem

The weak formulation of the boundary value problem (<ref>) reads as follows:
for given source function f ∈ L^2(Ω) 
find a function u∈ H_0^1(Ω) such that 
 the variational identity

    a(u,ϕ)=l_f(ϕ), ∀ϕ∈ 
    H^1_0(Ω),

is satisfied,
where the bilinear form a(·,·) and the linear form l_f(·)
are defined by 

    a(u,ϕ)=∫_Ωρ∇ u ·∇ϕ dx
        
    l_f(ϕ)=∫_Ωfϕ dx,

respectively. 

The given diffusion coefficient ρ∈ L^∞(Ω) is assumed to be
uniformly positive
and piece-wise (patch-wise, see below) constant. 
These assumptions ensure 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 which 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 which do not depend on the mesh-size h. In many cases,  
we will indicate on what may the constants depend on.  
Frequently, we will write a∼ b meaning that c   a≤ b ≤ C   a. 




 §.§ 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 <cit.>, 
<cit.> and <cit.> 
for a more detailed presentation. 
Consider,   𝒵={0=z_1 < z_2 < … <z_M=1} to be a partition 
 of I̅=[0,1] with I̅_j=[z_j,z_j+1],  j=1,…,M-1 to be the intervals
 of the partition.  
Let the integers p and n_1 denote the p spline degree and the number of the B-spline basis functions.
 
Based on 𝒵, we introduce the open knot vector  Ξ={0=ξ_1,  ξ_2, …,  ξ_n_1+p+1=1}, 
and the associated vector ℳ={m_1,…,m_M}  of knot multiplicities with m_1=m_M=p+1, i.e., 
 
    Ξ={0=ξ_1,…  ,ξ_m_1_=z_1, 
      ξ_m_1+1=…=ξ_m_1+m_2_=z_2,…,
     ξ_n_1+p+1-m_M,…,ξ_n_1+p+1=1_=z_M}.

The B-spline basis functions are defined by the Cox-de Boor formula, see, e.g.,  <cit.> and <cit.>,  

    B_i,p=   x-ξ_i/ξ_i+p-ξ_iB_i,p-1(x)+
    ξ_i+p+1-x/ξ_i+p+1-ξ_i+1B_i+1,p-1(x),
    with  B_i,0(x)=   
                     	1, if ξ_i ≤ x≤ξ_i+1,
    
                     	0, otherwise

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⊂ℝ^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 
 Ξ^k={0=ξ^k_1,  ξ^k_2, …, ξ^k_n_k+p+1=1},  k = 1,…,d, . 


 Given the knot vector Ξ^k in every direction k=1,…,d, 
 we construct the associated univariate B-spline basis functions, 
 𝔹̂_Ξ^k,p={B̂_1,k(x̂_k),…,B̂_n_k,k(x̂_k)}, 
 see, e.g., <cit.> for more details.
Accordingly,  the  B-spline basis functions of 𝔹̂_Ξ,k are defined by the tensor-product of the univariate B-spline basis functions, that is 

    𝔹̂_Ξ,p=⊗_k=1^d𝔹̂_Ξ^k,p
    =span{B̂_j(x̂)}_j=1^n=n_1·…· n_k·…· n_d,

where each  B̂_j(x̂) has the form

    B̂_j(x̂)=   B̂_j_1(x̂_1)·…·B̂_j_k(x̂_k)·…·B̂_j_d(x̂_d),  with B̂_j_k(x̂_k) ∈𝔹̂_Ξ^k,k.

 In the IGA framework, the computational domain  Ω is described as the image of Ω 
 under a B-spline, NURBS, etc.,  parametrization mapping of the form

    Φ: Ω→Ω,   
     x=Φ(x̂) = ∑_j=1^n 𝐂_jB̂_j(x̂)∈Ω,


 where 𝐂_j,  j=1,…,n are the control points and   x̂ =Φ^-1(x), 
 see Fig. <ref>(a). 
 Following the IGA methodology, <cit.>, <cit.>,
 the B-spline spaces for discretizing the PDE problem are defined by using the mapping given in (<ref>), i. e., 
we define the  B-spline space in Ω by 

    𝔹_Ξ,p:=span{B_j|_Ω: B_j(x)=
      B̂_j∘Φ^-1(x), for B̂_j∈𝔹̂_Ξ,p}.



  §.§.§  Multipatch representations and B-spline spaces
 
 
Our contribution here aims at developing a DG-IGA method appropriate for discretizing PDE models on 
non-conforming multipatch partitions  of the domain Ω. 
Let us suppose that the domain Ω is described as a union of N-subdomains  

    Ω=∪_i=1^N Ω_i,   with Ω_i∩Ω_j=∅, for  i≠ j,

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 space 𝔹̂_Ξ_i,p and the corresponding 
control points 𝐂_j^(i),  see (<ref>). Here  Ξ_i denotes the knot-vector related to Ω_i.   An illustration for N=2 is given 
in  Fig. <ref>(a). The subdomains Ω_i are referred to as patches. In an analogous way as in (<ref>), we define the physical patch-wise B-spline spaces 𝔹_Ξ_i,p for i=1,…,N. 
 We  define the global discontinuous B-spline space V_𝔹 with components on every 
 𝔹_Ξ_i,p

    V_𝔹:={ϕ_h∈ L^2(Ω): ϕ_h|_Ω_i∈𝔹_Ξ_i,p}.



	Assume that every Φ_i,  i=1,...,N is sufficiently smooth  and 
	there exist constants 0 < c < C  such that c ≤ | J_Φ_i | ≤ C, where J_Φ_i is the
        Jacobian matrix of Φ_i. 

 The components   of Ξ_i  form 
a mesh  T^(i)_h_i,Ω={Ê_m}_m=1^M_i in Ω,
where Ê_m are the micro-elements and h_i is the mesh size, which is
defined as follows. Given an element  Ê_m∈ T^(i)_h_i,Ω, 
we set h_Ê_m=diameter(Ê_m)
and the mesh size h_i is defined to be  h_i= max{h_Ê_m}.
We set
h=max_i=1,…,N{h_i}.
For every Ω_i, we construct a mesh T^(i)_h_i,Ω_i ={E_m}_m=1^M_i, 
whose vertices are the images of the vertices
 of the corresponding parametric mesh T^(i)_h_i,Ω
 under Φ_i. 
 
	The  meshes T^(i)_h_i,Ω are quasi-uniform, i.e.,
	there exist a constant θ≥ 1 such that 
	θ^-1≤h_Ê_m/h_Ê_m+1≤θ.
	Also, we assume that h_i ∼ h_j for 1 ≤ i≠ j ≤ N.

 
 



 §.§ Multipatch description of the computational domain



In many practical applications, the parametrization of a boundary represented domain Ω by a single B-spline  (NURBS) patch may not be posible.
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 
 <cit.>, 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 (<ref>), 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  Fig. <ref>(a). In this case 
the domain  Ω is  described as a union of 
two  non overlapping patches, see (<ref>),   i.e., 

 Ω = Ω_1 ∪Ω_2,      Ω_1 ∩Ω_2 = ∅,    
  with   F_12=∂Ω_1 ∩∂Ω_2,

where the interface F_12 coincides with  the  physical interface. 
We use the notation  
T_H(Ω):={Ω_1, Ω_2} for the  union (<ref>).  
For each Ω_i, i=1,2, there exists a  matching parametrization mapping such that Φ_i:Ω→Ω_i 
with Ω_i=Φ_i(Ω).
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 (<ref>), i.e.,

    u  :=u_1-u_2=0 on F_12,  and  ρ∇ u· n_F_12 :=(ρ_1∇ u_1-ρ_2∇ u_2)· n_F_12 = 0 on F_i12,

 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 (<ref>) can be ensured by  considering appropriate regularity assumptions on the solution u. 
We note that these type of multipatch representations have been considered in <cit.> and DG-IGA methods have been proposed for discretizing the problem (<ref>). 


Anyway,  for simplicity we develop  our analysis  based on Fig. <ref>. 
We introduce the appropriate spaces.  
Let  ℓ≥ 2 be  an integer, we define the broken Sobolev space

    H^ℓ(T_H(Ω))={u∈ L^2(Ω): u_i=u|_Ω_i∈ H^ℓ(Ω_i), for i=1,2}.


 We assume that the  solution u of (<ref>) belongs to  
 
 V= H^1_0(Ω)∩ H^2(Ω) ∩ H^ℓ(T_H(Ω)) with ℓ≥  2. 






 §.§ Problem statement



  §.§.§ Non-matching parametrized interfaces


Typically, the segmentation procedure will generate multipatch  representations that have
possibly non-matching interface parametrizations, <cit.>.
The result  is the existence of 
gap and overlapping regions in the multipatch representation of the domain Ω. 
In  <cit.> and <cit.>, 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  Fig. <ref>(b) and Figs. <ref>(a), (b).
Due to the  non-matching parametrization of the interior patch interfaces,
 a direct application of interface conditions, as those given in  
 (<ref>) 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 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.

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 
Fig. <ref>.  
Let  suppose that 

 Ω = Ω_1^* ∪Ω_2^*, 

where each patch has its own parametrization  Φ_1^*:Ω→Ω_1^* and 
 Φ_2^*:Ω→Ω_2^*, as it is shown in  Figs. <ref>(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^*. 
Next, 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 in a better way  our ideas, 
and to keep the notation to a minimum, e.g., the form of Jacobians, the form of face integrals etc. 
In Section <ref>,  we give details of implementing the proposed method to more complicated overlapping regions.  

	(a) Ω_1:= Ω^*_1. The face F_o12    is a an elementary face in the plane, 
	   and  it coincides with the physical interface,                i.e., F_o12=F_12, see (<ref>).

        (b) The face F_o21  can be described  as the set of points (x,y,z) 
    satisfying 
         
    0≤ x  ≤ x_M_o,  0  ≤ y ≤ y_M_o,  z=ζ_0(x,y),

           where x_Mo and y_M_o  are  real numbers, and   ζ_0(x,y) 
             is a  given  smooth  function, see Fig. <ref>.     

We note  that   we will discretize the PDE problem using the B-spline spaces defined in  
 Ω_1^* and Ω_2^*. 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 need to assign the points located on F_o12 to the diametrically opposite points located on F_o21. 
 Based on Assumption <ref>, 
we can construct   a parametrization for the face F_o21, i.e., a  mapping 
Φ_o12:F_o12→ F_o21, of the form

    x_o1∈ F_o12→Φ_o12(x_o1):=x_o2∈ F_o21,  with   Φ_o12(x_o1)=x_o1+ ζ_o(x_o1)n_F_o12,

where   n_F_o12 is the unit normal vector on F_o12 and ζ_o has the same form as in (<ref>),
and it  is a B-spline function with the same degree as the mapping Φ_2^*. More precisely, the face F_o21 
is the image of a face of ∂Ω̂ under  the mapping Φ_2^*. 
For the schematic illustration in Figs. <ref>(c),(d), we have F_o21=Φ_2^*(F̂_3). 
  Utilizing the mapping Φ_o12 given in (<ref>), we   consider 
each point
x_o2∈ F_o21 as an image  of a point x_o1∈ F_o12 under the  Φ_o12, 
see Figs. <ref>(a),(c).
Finally, we  introduce a parameter d_o, which  quantifies the width of the overlapping region Ω_o21, i.e., 

    d_o=   max_x_o1∈ F_o12 |x_o1- Φ_o12(x_o1)|.

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.,

    d_o ≤    C  h^λ,  with C>0 and  some  λ≥  1.

Based on this, we  assume that  n_F_o12≈ -n_F_o21, and 
define the mapping Φ_o21:F_o21→ F_o12  as
 
    Φ_o21(x_o2)= x_o1, with  Φ_o12(x_o1)=x_o2,

 where   Φ_o21 is  the inverse of Φ_o12. 

Our methodology can also be applied to the case where the interior faces of ∂Ω_o21 do not touch the boundary ∂Ω.


As we previously said, the face 
F_o12 is the image of a face of ∂Ω under  the mapping Φ_1^*, for example 
in Fig. <ref>(c) we have  F_o12=Φ_1^*(F̂_1). 
On the other hand, the face F_o21 is an interior curve for  Ω_1^*, see Figs. <ref>(a),(c). 
Thus, one could try to see F_o21 as an image of a curve F̂_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 (<ref>). 
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 
<ref>. 


In  Section <ref>, we present examples where the 
normal vector n_F_o12 is not constant across the face F_o12. 




§ 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 say   𝔹^*_Ξ_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 (<ref>).  
  The idea in our approach  is to
 introduce local (patch-wise)   problems a_i^*(u_i^*,ϕ_i^*)=l_i,f(ϕ_i^*) in every Ω_i^*, 
 with appropriate  bilinear forms a_i^*(·,·).  Using the triangle inequality, we split the error as  u-u^*_h_DG≤u-u^*_DG + u^*-u^*_h_DG.   Then we estimate every term separately. 


 

 §.§ The patch-wise variational problems

 Denote  T_H^*(Ω):={Ω_1^*,Ω_2^*},  let ℓ≥ 1 be an integer and let the B-spline spaces 𝔹^*_Ξ_i,p defined in Ω_i^*,  i=1,2. 
 Accordingly to the spaces  (<ref>) and (<ref>), we introduce the spaces 
 
    H^ℓ(T^*_H(Ω)):=   {u^*=(u^*_1,u^*_2):u_i^*∈ H^ℓ(Ω^*_i),  
     	                  u^*_i|_∂Ω_i^*∩∂Ω=0,for i=1,2},
    
     	H^ℓ_0(T^*_H(Ω)):=   {u^*=(u^*_1,u^*_2):u_i^*∈ H^ℓ_0(Ω^*_i), for i=1,2}, 
    
     	V^*_𝔹:=   {ϕ^*_h=(ϕ_1,h,ϕ_2,h): ϕ_i,h∈𝔹^*_Ξ_i,p, for i=1,2  }.

In order to proceed, we first define    
the DG-norm ._DG associated with  𝒯^*_H(Ω). 
For all v∈ V_h^*:=H^ℓ(T^*_H(Ω))+
V^*_𝔹,

    -0.0cmv^2_DG =  ∑_i=1^2(ρ_i∇ v_i^2_L^2(Ω^*_i) + 
                          ρ_i/hv_i^2_L^2(∂Ω^*_i∩∂Ω)
                           +  ∑_F_oij⊂∂Ω_i^*{ρ}/hv_i  ^2_L^2(F_oij)),  for 1≤ i≠ j≤ 2,
 
where  F_oij  are the interior faces related to  overlapping regions,
see Fig. <ref>(a), and {ρ}=1/2(ρ_i+ρ_j). 
We recall Assumption <ref>. 
On each Ω_i^*, i=1,2, we consider  the auxiliary problems:
 2

    
    (P1) 
                -div(ρ_1 ∇ u_1^*)     = f, in Ω_1^*, 
    
                                                u_1^*    = u_D    on ∂Ω_1^*∩∂Ω,
     
                                                u_1^*    = u     on    F_o12,

    
    (P2) 
                -div(ρ_2 ∇ u_2^*)     = f, in Ω_2^*, 
    
                                                u_2^*    = u_D    on ∂Ω_2^*∩∂Ω,
     
                                                u_2^*    = u     on   F_o21,

    

 and furthermore,  we  consider the corresponding  variational problems,  


    (P1V)   find u_1^*∈ H^1_∂Ω_1^*∩∂Ω(Ω_1^*) such that    
     u_1^*=u on F_o12, and       a_1^*(u_1^*,ϕ_1) = l^*_1, f(ϕ_1),  for  ϕ_1∈ H^1_0(Ω_1^*),

 where
 
    a_1^*(u_1^*,ϕ_1) =   ∫_Ω^*_1ρ_1∇ u_1^*·∇ϕ_1 dx, 
      and   l^*_1,f(ϕ_1)        = ∫_Ω^*_1 fϕ_1 dx,

 
    (P2V)   find u_2^*∈ H^1_∂Ω_2^*∩∂Ω(Ω_2^*) such that    
    
     u_2^*=u on F_o21, and        	a_2^*(u_2^*,ϕ_2) = l^*_2, f(ϕ_2),  for  ϕ_2∈ H^1_0(Ω_2^*),

 where
 
    a_2^*(u_2^*,ϕ_2) =   ∫_Ω^*_2ρ_2∇ u_2^*·∇ϕ_2 dx, 
      and     l^*_2,f(ϕ_2)        = ∫_Ω^*_2 fϕ_2 dx,

 
 
 By Assumption  <ref> and the definition of problem (P1), we can imply that  
 the solution u of (<ref>) satisfies the problem (P1V). The definition of (P2) and the fact that ρ_2≠ρ_1 on Ω_o21 imply  that u of (<ref>)   does not satisfy  the problem   (P2V). 
 
 According to Assumption <ref>, we make the following  assumption.
 
The solutions u_i^*, i=1,2  in  (<ref>)  belong to  H^ℓ(T^*_H(Ω)) 
with ℓ≥ 2.

In Appendix, see Subsection <ref>, we give an estimate for the distance of the solutions u and u^*. 


 §.§ The non-consistent   terms.

We multiply the problem (<ref>) by ϕ_2,h∈𝔹^*_Ξ_2,p, 
integrate over Ω_2^* and apply integration by parts, then after few calculations 
we find that 

    ∫_Ω^*_2   ρ_2∇ u_2^*·∇ϕ_2,h dx 
                        -   ∫_∂Ω^*_2ρ_2∇ u_2^*· n_∂Ω^*_2ϕ_2,h dσ
    
     =   ∫_Ω_o21ρ_1∇ u_2^*·∇ϕ_2,h dx 
                       + ∫_Ω_2ρ_2∇ u_2^*·∇ϕ_2,h dx  
      - ∫_∂Ω_2^*∩∂Ωρ_2∇ u_2^*· n_∂Ω_2ϕ_2,h dσ
     -0.73cm
     -   ∫_F_o21ρ_1∇ u_2^*· n_F_o21ϕ_2,h dσ
        - ∫_F_o12ρ_1∇ u_2^*· n_F_o12ϕ_2,h dσ  
        - ∫_F_o12ρ_2∇ u_2^*· (-n_F_o12)ϕ_2,h dσ
      
     +∫_Ω_o21(ρ_2-ρ_1)   ∇ u_2^*·∇ϕ_2,h dx - 
     ∫_F_o12 (ρ_2-ρ_1)∇ u_2^*· n_F_o12ϕ_2,h dσ 
      - ∫_F_o21 (ρ_2-ρ_1)∇ u_2^*· n_F_o21ϕ_2,h dσ
    
    = l^*_2,f(ϕ_2,h),

  and in a similar way,   multiplying 
  the problem (<ref>) by ϕ_1,h∈𝔹^*_Ξ_1,p, 
  we  have
  
  ∫_Ω^*_1ρ_1∇u_1^*·∇ϕ_h dx 
 - ∫_F_o12 ρ_1∇u_1^*·n_F_o12ϕ_h dσ- 
 ∫_∂Ω^*_1∩∂Ωρ_1∇u_1^*·n_∂Ω_1 ϕ_1,h dσ=
 l^*_1,f(ϕ_1,h).
  
  We define the forms
  

    -0.8cm
    a_2,h^*(u_2^*,ϕ_2,h):=∫_Ω^*_2   ρ_2∇ u_2^*·∇ϕ_2,h dx 
     - ∫_F_o21ρ_2∇ u_2^*· n_F_o21ϕ_2,h dσ - 
     ∫_∂Ω^*_2∩∂Ωρ_2∇ u_2^*· n_∂Ω_2ϕ_2,h dσ, 
     -0.8cm
    a_1,h^*(u_1^*,ϕ_1,h):=∫_Ω^*_1   ρ_1∇ u_1^*·∇ϕ_1,h dx 
     - ∫_F_o12ρ_1∇ u_1^*· n_F_o12ϕ_1,h dσ - 
     ∫_∂Ω^*_1∩∂Ωρ_1∇ u_1^*· n_∂Ω_1ϕ_1,h dσ

 
  and also
   
  
    -0.85cm  a_o,2(u_2^*,ϕ_2,h) =   ∫_Ω_o21ρ_1∇ u_2^*·∇ϕ_2,h dx 
                       + ∫_Ω_2ρ_2∇ u_2^*·∇ϕ_2,h dx  - 
       ∫_∂Ω_2^* ∩∂Ωρ_2∇ u_2^*· n_∂Ω_2^*ϕ_2,h dσ
    
     -    ∫_F_o21ρ_1∇ u_2^*· n_F_o21ϕ_2,h dσ
        - ∫_F_o12ρ_1∇ u_2^*· n_F_o12ϕ_2,h dσ  
        - ∫_F_o12ρ_2∇ u_2^*· (-n_F_o12)ϕ_2,h dσ
      
      a_res(u_2^*,ϕ_2,h) =    ∫_Ω_o21(ρ_2-ρ_1)∇ u_2^*·∇ϕ_2,h dx - 
      ∫_F_o12 (ρ_2-ρ_1)∇ u_2^*· n_F_o12 ϕ_2,h dσ
    
      -   ∫_F_o21 (ρ_2-ρ_1)∇ u_2^*· n_F_o21ϕ_2,h dσ.

  
By (<ref>), (<ref>) and  (<ref>), we get  that

a_2,h^*(u_2^*,ϕ_h) = a_o,2(u_2^*,ϕ_2,h)    +a_res(u_2^*,ϕ_2,h) 
= l^*_2,f(ϕ_2,h),

Also  for the solution u of (<ref>) we have that 

    ∫_Ω_o21ρ_1∇ u·∇ϕ_2,h dx 
                       + ∫_Ω_2ρ_2∇ u·∇ϕ_2,h dx  
     - ∫_F_o21ρ_1∇ u· n_F_o21ϕ_2,h dσ
    
    -    ∫_∂Ω_2^* ∩∂Ωρ_2∇ u· n_∂Ω_2^*ϕ_1,h dσ
     =
     l^*_2,f(ϕ_2,h).

From the conditions (<ref>),  the 
forms defined in (<ref>), (<ref>) and the relations  (<ref>) and (<ref>), we   derive that 
 

    a_o,2(u_2^*,ϕ_2,h) +a_res(u_2^*,ϕ_2,h) = a_o,2(u,ϕ_2,h)=l^*_2,f(ϕ_2,h)
    and    a_2,h^*(u,ϕ_2,h) - a_res(u,ϕ_2,h)  = l^*_2,f(ϕ_2,h).

 
 By a simple application of divergence theorem, we get


  a_res(u,ϕ_2,h) =    
  ∫_Ω_o21-div((ρ_2-ρ_1)∇u)ϕ_2,h dx   = 
  ∫_Ω_o21(ρ_2-ρ_1)/ρ_1f ϕ_2,h dx.
  
 
Finally, by   (<ref>) and (<ref>), we deduce that

    a_2,h^*(u,ϕ_2,h)
    + 
    ∫_Ω_o21(ρ_1-ρ_2)/ρ_1f ϕ_2,h dx=l^*_2,f(ϕ_2,h).


 Let ϕ_2,h∈𝔹^*_Ξ_2,p. There is a c>0 dependent on ρ but
 independent of u and Ω_o21 such that
 
    ϕ_2,h^2_L^2(Ω_o21)≤ c d_o h(∫_Ω_2^*|∇ϕ_2,h|^2 dx+
      {ρ}/h∫_F_o21ϕ^2_2,h dσ.)



 Let 𝐯=(0,yϕ_2,h^2). The divergence theorem for 𝐯 on Ω_o21 
  yields, 
 
    ∫_Ω_o21ϕ_2,h^2 dx+
     ∫_Ω_o212yϕ_2,h ∂_yϕ_2,h dx=∫_F_o21 yϕ^2_2,h dσ.

 Using that y≤ d_o and applying (<ref>) in (<ref>) we obtain 
 
    ϕ_2,h^2_L^2(Ω_o21)≤(ϵ^2∫_Ω_o21ϕ_2,h^2 dx+
      4/ϵ^2∫_Ω_o21d_o^2 |∇ϕ_2,h|^2 dx+
      d_o h 1/h∫_F_o21ϕ^2_2,h dσ)

Gathering similar terms and choosing ϵ appropriately small,  we get

    c_1,ϵϕ_2,h^2_L^2(Ω_o21)≤ c_2,ϵc_ρ d_oh(
      ∫_Ω_2^*ρ_2|∇ϕ_2,h|^2 dx+
       {ρ}/h∫_F_o21ϕ^2_2,h dσ),

where we  used  that  d_o^2≤ d_oh. Rearranging appropriately the constants in (<ref>) yields
(<ref>).
 


	Let f∈ L^∞(Ω), ϕ_2,h∈𝔹^*_Ξ_2,p and let u_2^* and u be the solutions of  (<ref>) and  (<ref>) respectively. 
	There are  constants c_1, c_ρ>0 dependent on F_o21 but independent
	of h such that
	
	
    ∫_Ω_o21fϕ_2,h dx    ≤ c_1 d_o f_L^∞(Ω_o21)ϕ_2,h_DG, 
    
    		|a_res(u,ϕ_2,h)|    ≤ c_ρ d_o f_L^∞(Ω_o21)ϕ_2,h_DG,

	


	It follows by from  the Cauchy-Schwartz inequality that
	
    ∫_Ω_o21fϕ_2,h dx ≤f_L^2(Ω_o21)ϕ_2,h_L^2(Ω_o21)≤ 
    	c_F_o21d_o^1/2f_L^∞(Ω_o21)ϕ_2,h_L^2(Ω_o21).

   Using (<ref>) in (<ref>), the required assertion follows easily. 

   Inequality (<ref>) follows immediately from (<ref>) and (<ref>).




 §.§ The discrete problem

In this section, we use the bilinear forms given in (<ref>) to 
define the patch-wise discrete problems. 
Using the  conditions on F_o21 and F_o12, which are given in (P1) and (P2), the 
Assumption <ref>, we imply  the following  interface conditions 
 
    ρ_1∇ u_2^*· n_F_o12 = ρ_2∇ u^*_2· n_F_o12 on F_o12,    and  
       u_1^*-u^*_2=0 on F_o21.

Next, using Taylor expansions, 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σ appearing in (<ref>). 


  §.§.§ Taylor expansions  

Let x,y∈Ω^*_2 and let  f∈ C^m≥ 2(Ω^*_2).
We recall 
Taylor's formula with integral remainder


    f(y)=   f(x) +∇ f(x)·(y-x) + R^2f(y +s(x-y)),
    
     f(x)=   f(y) -∇ f(y)· (y-x)  + R^2 f(x +s(y-x)),


where R^2 f(y +s(x-y)) and R^2 f(x +s(y-x)) are the second order remainder terms defined by


    R^2 f(y +s(x-y))=   ∑_|α|=2(y-x)^α2/α !∫_0^1sD^α  f(y +s(x-y)) d s, 
    
        R^2 f(x +s(y-x))=   ∑_|α|=2(x-y)^α2/α !∫_0^1sD^α  f(x +s(y-x)) d s.


By (<ref>) it follows that  


    ∇ f(y) · (y-x)  =    ∇ f(x) · (y-x)  +(  R^2 f(x +s(y-x))+ R^2f(y +s(x-y))  ), 
    
    	  -f(x) =    -f(y)+∇ f(x)· (y-x)  + R^2f(y +s(x-y)).




  §.§.§ Modifications of the   fluxes on ∂Ω_o21 

To illustrate  the use of (<ref>) to (<ref>) in our analysis, we consider the simple case 
of Fig. <ref>(a). 
 Let the points  x_o1∈ F_o12 and x_o2∈ F_o21 be such that 
 x_o2=Φ_o12(x_o1) as in Fig. <ref>(a). 
 These points  play the role of the points x and y in (<ref>). Then
  for a smooth function f we have
 

	f(x_o1) =    f(x_o2)   +∇f(x_o2)·(x_o1-x_o2)  + R^2 f(x_o1 +s(x_o2-x_o1)) 

	         =     f(Φ_o12(x_o1)) 
	             +∇f(Φ_o12(x_o1)) ·(x_o1-x_o2)  + R^2 f(x_o1 +s(x_o2-x_o1)).


 Now denoting r_o12=x_o1-x_o2 and using the assumption that r_o12=-r_o21, 
 see Section <ref>,  we  obtain that
 n_F_o12=r_o12/|r_o12|=-n_F_o21.
For keeping notation simple,  we  denote the Taylor's residuals  as 
R^2u^*_x_o1:=R^2 u^*(x_o1 +s(x_o2-x_o1)) and 
R^2u^*_x_o2:=R^2 u^*(x_o2 +s(x_o1-x_o2)).
Using   (<ref>) and interface conditions (<ref>),  
we  modify the  fluxes in (<ref>) as follows




∫_F_o21    ρ_2∇u_2^*(x_o2)·n_F_o21  ϕ_2,h dσ- 
{ρ}/h∫_F_o21 (u_2^*(x_o2)-u_1^*(x_o2))ϕ_2,h dσ

                          =   
  ∫_F_o21 1/2ρ_2∇u_2^*(x_o2) ϕ_2,h   dσ-{ρ}/h∫_F_o21(u_2^*(x_o2)-u_1^*(Φ_o21(x_o2))ϕ_2,h  dσ

 +   
 ∫_F_o21{ρ}/h(|r_o12|∇u_2^*(x_o2)·n_F_o21 +
 R^2u^*_x_o2 ) ϕ_2,h   dσ,
 
 
 where   {ρ}=1/2(ρ_1+ρ_2). Similarly,  we have


∫_F_o12 ρ_1∇    u_1^*(x_o1)·n_F_o12  ϕ_1,h dσ=
∫_F_o12   ρ_1∇u_1^*(x_o1)·n_F_o12 ϕ_1,h dσ-  
 {ρ}/h∫_F_o12 (u_1^*(x_o1)-u_1^*(x_o1))ϕ_1,h dσ

 =
 ∫_F_o12      ρ_1∇u_1^*(x_o1)·n_F_o12 ϕ_1,h dσ-  
 {ρ}/h∫_F_o12 (u_1^*(x_o1)-u_2^*(Φ_o12(x_o1)))ϕ_1,h dσ

+
 ∫_F_o12   {ρ}/h(|r_o21|∇u_1^*(x_o1)·n_F_o12 +
 R^2u^*_x_o1 ) ϕ_1,h   dσ.
  





  §.§.§ The global modified form

We consider the global bilinear form a^*(·,·): V_h^*× V_𝔹^* →ℝ, which is 
formed by the  contributions of a_i,h^*(·,·), i=1,2
given in (<ref>) and  the flux forms  given in (<ref>), that is 

    a^*(u^*,ϕ_h)=a_2,h^*(u_2^*,ϕ_2,h)+a_1,h^*(u_1^*,ϕ_1,h)=
    ∫_Ω_1^*ρ_1∇ u_1^*·∇ϕ_1,h dx  
    +∫_Ω_2^*ρ_2∇ u_2^*·∇ϕ_2,h dx 
    
    -   ∫_∂Ω_1^* ∩∂Ωρ_1∇ u_1^*· n_∂Ω_1^*ϕ_1,h dσ  - 
       ∫_∂Ω_2^* ∩∂Ωρ_2∇ u_2^*· n_∂Ω_2^*ϕ_2,h dσ
    
       + ρ_1/h∫_∂Ω_1^* ∩∂Ω(u_1^*-u_D)ϕ_1,h dσ 
       + ρ_2/h∫_∂Ω_2^* ∩∂Ω(u_2^*-u_D)ϕ_2,h dσ
    
       -∫_F_o12( ρ_1∇ u_1^*(x_o1)· n_F_o12       
      +{ρ}/h(u_1^*(x_o1)-u_2^*(Φ_o12(x_o1)))ϕ_1,h dσ
     
       -∫_F_o21(ρ_2∇ u_2^*(x_o2)· n_F_o21  
     +{ρ}/h(u_2^*(x_o2)-u_1^*(Φ_o21(x_o2))) ϕ_2,h dσ
      
      -∫_F_o21{ρ}/h(|r_o12|∇ u_2^*(x_o2)· n_F_o21 +
     R^2u^*_x_o2 )ϕ_2,h dσ
    
     -∫_F_o12{ρ}/h(|r_o21|∇ u_1^*(x_o1)· n_F_o12 +
             R^2u^*_x_o1 ) ϕ_1,h dσ.


	Note that  the exact 
	solution u has similar  regularity properties to the solution u^*, 
	see Assumption <ref>, and thus  we can derive for u 
	an analogous formulation as this  in (<ref>).   



  §.§.§ The DG-IGA scheme.

In view of (<ref>),  we    define the  forms
A_Ω_i^*(·,·):V_h^*× V_𝔹^* →ℝ, 
R_Ω_o21(·,·):V_h^*× V_𝔹^*  →ℝ,
and the linear functional
l_f, Ω_i^*:V_𝔹^* →ℝ by
-0.4cm
 

    A_Ω_i^*(u^*,ϕ_h) =   ∑_i=1^2(∫_Ω^*_iρ_i∇ u_i^*·∇ϕ_i,h,dx - 
       ∫_∂Ω^*_i ∩∂Ωρ_i∇ u_i^*· n_∂Ω^*_iϕ_i,h dσ
    
     -   ∑_F_oij⊂∂Ω_i^*∫_F_oijρ_i∇ u_i^*· n_F_oijϕ_i,h  -
      η{ρ}/h(u_i^*-u_j^*)ϕ_i,h  dσ) 
    
      -   ∑_i=1^2 ∫_∂Ω_i^* ∩∂Ωρ_i∇ u_i^*· n_∂Ω_i^*ϕ_i,h dσ,    1≤ i≠ j≤ 2,

 -0.2cm
  

R_Ω_o21(u^*,ϕ_h) =   ∫_F_o21  
                           -{ρ}/h(|r_o12|∇u_2^*(x_o2)·n_F_o21 +
                            R^2u^*_x_o2 )ϕ_2,h dσ

                           -∫_F_o12      
                        {ρ}/h(|r_o21|∇u_1^*(x_o1)·n_F_o12 +
                          R^2u^*_x_o1 ) ϕ_1,h   dσ

 l_f,Ω_i^*(ϕ_h) =     ∑_i=1^2 ∫_Ω_i^*fϕ_i,h dx,
 


where η >0 is a parameter that is going to be determined later. 
Based on  the forms defined  in (<ref>), 
we introduce the discrete  bilinear form
A_h(·,·):V_𝔹^*× V_𝔹^* →ℝ
 and the linear form F_h: V_𝔹^* →ℝ
as follows

    A_h(u_h^*,ϕ_h) = A_Ω_i^*(u_h^*,ϕ_h) +
    	∑_i=1^2ηρ_i/h∫_∂Ω_i^* ∩∂Ωu_i,h^*ϕ_i,h dσ,


    F_h(ϕ_h) = l_f,Ω_i^*(ϕ_h)+
    	∑_i=1^2ηρ_i/h∫_∂Ω_i^*∩∂Ωu_D ϕ_i,h dσ.
 
Finally, the DG-IGA scheme reads as follows:
find u_h^*∈ V_𝔹^* such that

    A_h(u_h^*,ϕ_h) = F_h(ϕ_h),   for all ϕ_h∈ V_𝔹^*.



From the relations (<ref>), (<ref>), the Remark <ref> and 
the forms given in (<ref>) and in (<ref>),  we can derive that 
	 
    a_2,h^*   (u_2^*,ϕ_2,h)+a_1,h^*(u_1^*,ϕ_1,h) 
      
    
    	  =     
    	  a_o,2(u_2^*,ϕ_2,h) +a_res(u_2^*,ϕ_2,h) +a_1,h^*(u_1^*,ϕ_1,h) 
    
    
     =   	 A_h(u^*,ϕ_h)+	R_Ω_o21(u^*,ϕ_h)  
    
     =    a_o,2(u,ϕ_2,h) + a_1,h^*(u,ϕ_1,h) 
    
    
    =    a_2,h(u,ϕ_2,h) -a_res(u,ϕ_2,h) +a_1,h^*(u,ϕ_1,h) 
     
    =     A_h(u,ϕ_h)+	R_Ω_o21(u,ϕ_h)  -a_res(u,ϕ_2,h) 
    
    =     F_h(ϕ_h),  8 cm      for ϕ_h:=(ϕ_1,h,ϕ_2,h)∈ V_𝔹^*.
	


Below,  we quote  few  results  that  are useful for our  error analysis. 
For the proofs we refer to
<cit.>, <cit.> and <cit.>. 

Under the assumption  (<ref>),  there exist  positive constants 
C_1 and C_2 independent of h such that the estimates
 
  | R_Ω_o21(u,ϕ_h)|≤C_1 𝒦_o(u)
  ϕ_h_DG   h^λ-0.5,      | R_Ω_o21(u^*,ϕ_h)|≤C_2 𝒦_o(u^*)
  ϕ_h_DG   h^λ-0.5,   
   hold    for the solutions u^* and u,
   and ϕ_h ∈ V_𝔹^*, 
   where   𝒦_o(v)= ∇ v_L^2(∂Ω_o21)+
   ∑_|α|=2|D^α v|_L^2(Ω_o21).
  
  
  
	The bilinear form A_h(·,·) in (<ref>) is bounded and elliptic on V_𝔹^*, i.e., 
	 there are positive constants C_M and C_m such that the estimates
	
    A_h(v_h,ϕ_h) ≤ C_M v_h_DGϕ_h_DG  and    A_h(v_h,v_h) ≥ C_m v_h^2_DG,

	hold for all v_h, ϕ_h∈ V_𝔹^* 
	 provided that η is sufficiently large.

   Let the  assumption  (<ref>) and let β=λ-1/2. Then there is a constant 
 C_* > 0 depending on the parametrization mappings but  independent of h such that the inequality 

    A_h(v,ϕ_h)≤    C_*( (v^2_DG + ∑_i=1^2
       h ρ_i^1/2∇ v^2_L^2(∂Ω^*_i))^1/2)
       ϕ_h_DG,

 holds  for all   (v,ϕ_h) ∈ V_h^* × V_𝔹^* and (v,ϕ_h)∈ (V+ V_𝔹^*)× V_𝔹^*. 

 
	Recall the definition of the pair function spaces in (<ref>). 
	In view of the form of A_h(·,·)  and applying (<ref>), we have
	
    |∑_i=1^2(∫_Ω^*_iρ_i∇ v_i·∇ϕ_i,h dx | ≤(∑_i=1^2 ρ_i^1/2∇ v_i^2_L^2(Ω_i^*))^1/2(∑_i=1^2 ρ_i^1/2∇ϕ_i,h^2_L^2(Ω_i^*))^1/2.

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 (<ref>), we obtain
      
    |   ∫_F_oijρ_i∇ v_i· n_F_oijϕ_i,h dσ|  ≤ C_i
    	   ∫_F_oijh^1/2|ρ_i^1/2∇ v_i|  ({ρ}/h)^1/2|ϕ_i,h| dσ
    ≤   
    	C_i ( h^1/2ρ_i^1/2∇ v_i_L^2(F_oij))
    	(η{ρ}/hϕ_i,h^2_L^2(F_oij))^1/2≤ 
    	C_i ( h^1/2ρ_i^1/2∇ v_i_L^2(F_oij))
    	ϕ_h_DG,

	for 1≤ i≠ j ≤ 2.
	Also, we have  
	
    -0.5cm 
    	|η{ρ}/h∫_F_o12(v_1-v_2(Φ_o12))ϕ_1,h  dσ| 
    	 ≤   
    	 2  (   η{ρ}/h∫_F_o12 v^2_1+ v^2_2(Φ_o12 )|J_Φ_o12|/|J_Φ_o12| dσ)^1/2(   η{ρ}/hϕ_1,h^2_L^2(F_o12))^1/2 2cm
    ≤ C_J_Φ_o12(   η{ρ}/hv_1  ^2_L^2(F_o12) +    η{ρ}/hv_2 ^2_L^2(F_o21)))^1/2(   η{ρ}/hϕ_1,h^2_L^2(F_o12))^1/2≤  C_J_Φ_o12v_DG ϕ_h_DG,
    where |J_Φ_o12| 
    	 is the measure of the Jacobian of Φ_o12. In the same way, we show|η{ρ}/h   ∫_F_o21(v_2-v_1(Φ_o21))ϕ_h  dσ| 
    	 ≤ C_J_Φ_o21v_DG ϕ_h_DG.

Gathering together the above bounds, we  show (<ref>). For the case where (v,ϕ_h)∈ (V+ V_𝔹^*) × V_𝔹^* 
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 (<ref>). Let 
 v ∈ H^ℓ(T^*_H(Ω))  with ℓ≥ 2.
Under Assumptions <ref>, and using the results of
<cit.> and <cit.>, we can construct a quasi-interpolant 
Π^*_h v :=(Π^*_1,h v_1,(Π^*_2,h v_2)  ∈ V_𝔹^* such that
the 
estimates 

    ∑_i=1,2|v-Π^*_h v|_H^1(Ω_i^*)≤      h^s∑_i=1,2C_1,iv_H^ℓ(Ω_i^*), 
    ∑_i=1,2|v-Π^*_h v|_L^2(∂Ω_i^*)≤     h^s-1/2∑_i=1,2C_2,iv_H^ℓ(Ω_i^*),

hold, where s=min(ℓ-1,p) and  the C_1,i, C_2,i depend on  p, Φ^*_i, θ   but not on h.  
 
 Let v ∈ H^ℓ(T^*_H(Ω)) with ℓ≥ 2  
 and let Π^*_h v be  as in (<ref>).   
 Then  there exist  constants C_i>0, i=1,2, depending on p, Φ^*_i, i=1,2 and 
 the quasi-uniformity of the meshes  but not on   h  such that 
  
    (v-Π^*_h v^2_DG +∑_i=1^2hρ_i^1/2∇ (v-Π^*_h v)^2_L^2(∂Ω_i^*))^1/2   ≤∑_i=1^2 C_i h^sv_H^ℓ(Ω_i^*),

  where s=min(ℓ-1,p).
 
  
  The   estimate (<ref>) can be shown using trace inequality and the  estimates (<ref>), 
  see  details in   Lemma 10 in <cit.>.
   See also  <cit.> and <cit.>.
    
   
   
  

Let β=λ -1/2 and   d_o = h^λ with λ≥ 1. Let u^*∈ H^ℓ(T^*_H(Ω)) with ℓ≥ 2 be the solution of the  problems in  (<ref>), and
let u_h^*∈ V_𝔹^* be the corresponding DG-IGA solution of (<ref>). 
Then the  error estimate

    u^*-u_h^*_DG≲  h^r( ∑_i=1^2u^*_H^ℓ(Ω_i^*)),

holds,  where r=min(s,β) with  s=min(ℓ-1,p).


Let z_h∈ V_𝔹^*. We set u^*_h-z_h=ϕ_h. The properties 
(<ref>),  (<ref>) 
of A_h(·,·) and (<ref>) imply 

    c_m   u^*_h-z_h^2_DG≤ 	 A_h(u^*_h-z_h,ϕ_h)=
               A_h(u^*,ϕ_h)+	R_Ω_o21(u^*,ϕ_h)
    	  - A_h(z_h,ϕ_h)
             
    =    A_h(u^*-z_h,ϕ_h)+R_Ω_o21(u^*,ϕ_h) 
    ≤    C_*( (u^*-z_h^2_DG + ∑_i=1^N
       h ρ_i^1/2∇(u^*-z_h)^2_L^2(∂Ω^*_i))^1/2)
       ϕ_h_DG
    
       +   C_2𝒦_o(u^*)
      ϕ_h_DG   h^λ-0.5,

where the bound (<ref>) has been used previously. 
Setting in (<ref>) z_h=Π^*_h u^*, and then  
using the triangle inequality c_mu^*_h-u^*_DG≤ c_mu^*_h-Π^*_h u^*_DG +
                               c_mu^*-Π^*_h u^*_DG together with the estimate in  (<ref>), we 
                               derive (<ref>).
    
   


  §.§.§ Main error estimate

	The estimate given in (<ref>) concerns the distance between the DG-IGA solution 
	u_h^*∈ V_𝔹 and the solution
	u^*∈ H^ℓ(T^*_H(Ω))  of
	the  problems  in (<ref>).  Below we give an estimate between 
	the solution u of (<ref>) and the DG-IGA solution u_h^*. 
	In the proof of this result we need the following interpolation estimate for v∈ V
		
    (v-Π^*_h v^2_DG +∑_i=1^2hρ_i^1/2∇ (v-Π^*_i,h v)^2_L^2(∂Ω_i^*))^1/2   ≤∑_i=1^2 C_i h^sv_H^ℓ(Ω_i^*),

  where the quasi-interpolant Π_h^*v=(Π_1,hv, Π_2,hv) is defined in (<ref>) and  s=min(ℓ-1,p).
The proof of (<ref>) is provided in the Appendix. 


	 Let  u be the solution of (<ref>) and let  Assumption <ref> with ℓ≥ 2.  We suppose further that d_o=h^λ, λ≥ 1 is the width of Ω_o21. 
	 The following error estimate holds
	
    u-u^*_h_DG≤C(h^s∑_i=1^2
    	(u_H^ℓ(Ω_i^*)+u^*_i_H^ℓ(Ω_i^*))
    	       +d_of_L^2(Ω) +h^β(𝒦_o(u)+𝒦_o(u^*))),

where β=λ -1/2,  s=min(ℓ-1,p), the constant C 
depends on the constants  in (<ref>),  (<ref>) and (<ref>), and 
 𝒦_o has the form given in  Lemma  <ref>. 

 
          Let z_h∈ V^*_𝔹  and let ϕ_h=u_h^*-z_h.
          By the definition of the discrete DG-IGA scheme in (<ref>),  the properties of A_h(·,·) and the Remark <ref>  we have
         
    c_m   u_h^*-z_h^2_DG≤           A_h(u^*_h - z_h,ϕ_h) - A_h(u^*,ϕ_h) - R_Ω_o21(u^*,ϕ_h) + F_h(ϕ_h)
              -A_h(Π_h^*u^*,ϕ_h)+A_h(Π_h^*u^*,ϕ_h)
     
              =    
              A_h(u_h^*-Π_h^* u^*,ϕ_h) +  A_h(u^*-Π_h^* u^*,ϕ_h) 
    
              +      A_h(-z_h,ϕ_h ) + 
              A_h(u,ϕ_h) -a_res(u,ϕ_2,h)+R_Ω_o21(u,ϕ_h)-R_Ω_o21(u^*,ϕ_h)
    
              =      A_h(u_h^*-Π_h^* u^*,ϕ_h) +  A_h(u^*-Π_h^* u^*,ϕ_h)+ A_h(u-z_h,ϕ_h)
    
              -     a_res(u,ϕ_2,h)+R_Ω_o21(u,ϕ_h)-R_Ω_o21(u^*,ϕ_h) 
    ≤     C_Mu^*_h-Π_h^*u^*_DGϕ_h_DG 6cm            
    by (<ref>),(<ref>), (<ref>), (<ref>)
    
       +     C_*( (u^*-Π_h^*u^*^2_DG + ∑_i=1^N
       h ρ_i^1/2∇(u^*-Π_h^*u^*)^2_L^2(∂Ω^*_i))^1/2)
       ϕ_h_DG
    
     +      C_*( (u-z_h^2_DG + ∑_i=1^N
       h ρ_i^1/2∇(u-z_h)^2_L^2(∂Ω^*_i))^1/2)
       ϕ_h_DG
    
     +     c_2d_of_L^2(Ω)ϕ_h_DG +C_2(𝒦_o(u^*) +𝒦_o(u))
      ϕ_h_DG   h^β

   Setting  z_h=Π^*_h u into (<ref>),  using  (<ref>), (<ref>), 
   and (<ref>) and gathering together the similar terms
   we deduce that
   
    c_mu_h^*-Π^*_h u_DG≤   ∑_i=1^2 C_i h^su_H^ℓ(Ω_i^*) +
       ∑_i=1^2 C_i h^su^*_i_H^ℓ(Ω_i^*)
    
       +    
       c_2d_of_L^2(Ω) +C_2(𝒦_o(u^*) +𝒦_o(u))   h^β

 Applying the triangle inequality 

    u-u^*_h_DG≤u-Π^*_h u_DG+Π^*_h u-u^*_h_DG,

the desired estimate follows.




§ IMPLEMENTATION AND NUMERICAL TESTS




 §.§ 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.
 

 
 2
  []
    < g r a p h i c s >

  []
    < g r a p h i c s >


  (a) Configuration of the faces and the edges on ∂Ω_o12 and their
 corresponding edges on ∂Ω which are used to compute the interface integrals,  (b) an example of an overlapping region with more than two faces. The relative edges on the opposite faces must again much.   
         
 

 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 Fig. <ref>. 
 

 For the computation of the numerical flux terms of the 
  DG-IGA scheme given  in (<ref>), 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 (<ref>) and (<ref>) 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 Fig. <ref>(a), the edge e_2^1 of F_o21 is assigned to e_2^2 of F_o12.  
  In Fig. <ref>(a) the Gauss point are denoted by 
  x_q and y_q correspondingly. 
   The edge  e_2^1 is the image of the 
  edge ê_2^1 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  Fig.  <ref>(a), the other associated edge is located on face F̂_1 and the corresponding Gauss point is denoted by x̂_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 nonmatching interfaces in a similar way. Note also that the previous approach can be easily combined with the 
   adaptive  numerical  quadrature methods presented in <cit.>, 
   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 Fig. <ref>(b).   We consider again the case where the maximum number of the overlapping patches is two.  For the  example shown in Fig. <ref>(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 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., <cit.>.  In an ongoing work we are extending the present methodology to treat these cases. We also  are  constructing  domain-decomposition methods, <cit.>,  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. The first results of this work are included in  <cit.>.  

   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_o121/2(ρ_1 ∇ u_1,h +ρ_2 ∇ u_2,h(Φ_o12) )
   · n_F_o12ϕ_1,h +η{ρ}/h(u_1,h -u_2,h(Φ_o12) ) ϕ_1,h ds., see <cit.>, <cit.>.
   


 §.§ 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 (<ref>).  All examples have been performed using  second degree (p=2)  B-spline spaces. 
We present the asymptotic behavior of  the error convergence rates  for  widths d_o=h^λ  
 with  λ∈{1,2,2.5,3}. Every example has been solved
 applying several mesh refinement steps with …,h_i,h_i+1,…, satisfying Assumption <ref>.
 The numerical convergence rates r have been  computed by the ratio 
 r =ln(e_i/e_i+1)/ln (h_i /h_i+1 ),  i=1,2,…, where the error e_i:=u-u^*_h_DG is always computed  
on the  meshes ∪_i=1^2 T^(i)_h_i,Ω_i^*.
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  (<ref>) and (<ref>) 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 <ref>. 
    In any test case, the   overlap regions  are artificially created by moving 
the control points, which are related to the interfaces F_ij, 
in the direction of n_F_ij or of -n_F_ij.  
All tests have been performed in G+SMO <cit.>,
which is a generic object-oriented C++ library for IGA computations,
<cit.>. 
In Section 3, we developed and provided a rigorous analysis for the DG-IGA method (<ref>) which includes a 
 non-symmetric numerical flux. 
In the materialization of the method, we utilized the associated symmetrized version the numerical flux, <cit.>.
For solving the resulting linear system, we use the DG-IETI-DP method presented  in  <cit.>, see also <cit.> for an analysis of the method and <cit.> for results on parallel scalability. 
Although in the analysis, we consider  meshes with similar quasi-uniform patch-wise properties,
it is known that the introduction of DG techniques
on the subdomain interfaces makes the use of non-matching and non-uniform meshes  easier, see <cit.>.
Keeping a constant linear relation between the sizes of the different patch meshes, 
the approximation properties
of the method are not affected,  <cit.>. 
In the examples  below, we exploit this advantage of the DG methods and 
first solve two-dimensional problems considering
non-matching meshes.  The   convergence rates are expected  to be the same as those displayed in Table <ref>. 




 §.§ Two-dimensional numerical examples

The control points with the corresponding knot vectors of the domains given in Example 1-3 are available under the names ,  and  as  files in 
G+SMO[G+SMO: https://www.gs.jku.at/trac/gismo].



  
Example  1: uniform diffusion coefficient ρ_i=1,  i=1,…,N. 
The first numerical example is a simple test case demonstrating
the  applicability of the proposed technique for constructing the 
DG-IGA scheme on segmentations 
including   overlaps  with general shape. 
The domain  Ω with the N=21 subdomains Ω_i^* 
 and the initial mesh  are shown in Fig. <ref>(a). We note that 
 we consider non-matching meshes across the interior  interfaces. 
The Dirichlet boundary condition and the right hand side f are determined by the exact solution 
u(x,y)=sin(π(x+0.4)/6)sin(π(y+0.3)/3)+x+y. 
 In this example, we consider the homogeneous diffusion case, i.e., 
ρ_i=1 for all  Ω_i^*, i=1,…,N.
We performed four groups of  computations, where for every group 
 the   maximum size of  d_o was defined to be 𝒪(h^λ), 
 with λ∈{1,2, 2.5,3}. 
 In Fig. <ref>(b) we present the discrete solution for d_0 = h.
 Since we are using second-order (p=2) B-spline space, based on Table 
 <ref>, we expect optimal convergence rates for λ =2.5 and λ = 3.  
 The numerical convergence rates for several levels of mesh refinement are plotted in Fig. <ref>(c). They are in very good	agreement with the theoretically predicted estimates given in Theorem <ref>, 
 see also Table <ref>. We observe that we have optimal rates r for the cases where λ≥ 2.5 and
 sub-optimal for the rest values of λ.
   
 

  
Example  2: different diffusion coefficients ρ_1≠ρ_2.  
 In the second example,  we consider a rectangular domain
 Ω, that is  described as a union of  N=12 patches, see Fig. <ref>(a).
 Here,  we study the case of having smooth solutions in each Ω_i^*
 but discontinuous coefficient,  i.e.,  we set ρ_i = 3π/2 for the patches belonging to
 half plane x≤ 0 and we set ρ_i= 2 for the rest patches  
 according to the pattern in Fig. <ref>(a). 
 By this example,  we numerically validate the 
 predicted convergence rates on 𝒯_H^* with  overlaps,  
for the case of having  smooth solutions and discontinuous coefficient ρ. 
 The exact solution is given by the formula

    u(x,y)=
        sin(π(2x+y))     if  x<0 
    sin(π(3π/2 x+y))     otherwise.
 
The  boundary conditions and the source function f are  determined by (<ref>).  
Note that, we have  u  |_F_ij=0 as 
well ρ∇ u  |_F_ij· n_F_ij=0 for 
 all the interior physical  interfaces F_ij.  
The problem has been solved   on a sequence of meshes  with h_0,...,h_i, h_i+1,...,
   following a sequential refinement process, i.e., h_i+1=h_i/2,
   where we set  d_o=h_i^λ, with λ∈{1,2,2.5,3}. 
   For the numerical tests, we use    B-splines of the degree p=2.   
   Hence, we expect optimal rates for λ≥  2.5.  
   In Fig. <ref>(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  Fig. <ref>(c). For all test cases, we can observe  that our theoretical results presented in Table <ref> are confirmed. 

  


  
Example  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 Fig. <ref>(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 (<ref>). 
We solved the problem using B-splines of degree p=2.  
In Fig. <ref>(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 Fig. <ref>(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 Fig. <ref>(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 Fig. <ref>(a).
  The knot vector in z-direction  is simply Ξ_i^3={0 , 0 ,0 , 0.5 ,1, 1, 1} with i=1,2.
  We solve the problem  using  matching meshes, as depicted in  Fig. <ref>(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 Ξ_i=1,2^k=1,2,3  together with the associated control nets can be found in
   G+SMO library in the file . 
  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 of the domain, we have to provide cuts though the domain in order to visualize them, cf.  Fig. <ref>(b).
  The Dirichlet boundary conditions u_D and the right hand side f, see (<ref>), are chosen such that the exact solution is 	
  	
    u(x,y,z)=
     	                  sin(π/2(x+y))    if  (x,y)∈Ω_1, 
    
                                  e^sin(x+y)   if (x,y)∈Ω_2.
 
  	with diffusion coefficient ρ={1,π/2}. Note that the interfaces conditions (<ref>) are satisfied. 
  	The two physical subdomains, the initial matching meshes and the 
  	exact solution are illustrated in Fig. <ref>(a).   We construct an overlap region with d_o=0.5 and solve the
  	problem using p=2 B-spline functions. In Fig. <ref>(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 Fig. <ref>(c). We observe from the plots that the  rates r are  in agreement with the   rates  predicted by the theory, see estimate (<ref>) and Table <ref>. 
  



§ CONCLUSIONS

In this article,  we have proposed  and analyzed a DG-IGA scheme for discretizing linear,
second-order, diffusion problems 
on IGA  multipatch representations with small   overlapping regions. 
This type of dmultipatch representations   lead   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 were 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. 


 

§ ACKNOWLEDGMENTS

The authors wish to thank Prof. Ulrich Langer, Prof. Bert Jüttler and  Prof.  Dirk Pauly for many interesting discussions. 
This work was supported by the Austrian Science Fund (FWF) under the grant NFN S117-03 and  W1214-N15, project DK4.



§ APPENDIX.



 §.§ A bound for the extra non-consistent term. 



Comparing the relations given in  (<ref>) and (<ref>)
we can see  that there is  an extra  term  -(ρ_2-ρ_1)∇ u_2^* in Ω_o21, 
which is a non consistent  term.  We derive   below a bound  for this term. 

Let ϕ∈ H^1_0(Ω_2^*). By a simple computations on the forms in (<ref>), we have that

a_2^*(u_2^*,ϕ_h)=   ∫_Ω_o21ρ_1∇u_2^*·∇ϕ dx 
                   + ∫_Ω_2ρ_2∇u_2^*·∇ϕ dx  
    -∫_∂Ω_2∩∂Ωρ_2∇u_2^*·n_∂Ω_2 ϕ dσ

 -    ∫_F_o21 ρ_2∇u_2^*·n_F_o21ϕ dσ=∫_Ω_o21(ρ_1-ρ_2)∇u_2^*·∇ϕ dx + l_2,f^*(ϕ). 
  
 On the other hand, under the Assumption <ref>, we immediately have that
  
  a_o,2(u,ϕ_2) =     ∫_Ω_o21ρ_1∇u·∇ϕ dx +
                       ∫_Ω_2ρ_2∇u·∇ϕ dx 

                       -   ∫_F_o21  ρ_1 ∇u  ·n_F_o21ϕ dσ-
        ∫_∂Ω_2^*∩∂Ω  ρ_2 ∇u  ·n_∂Ω_2ϕ dσ= l_2,f^*(ϕ).
 
 Subtracting (<ref>) from (<ref>) and using ϕ|_∂Ω_2^*=0
 we  obtain
  
  	
  	∫_Ω_o21ρ_1∇(u_2^*-u)·∇ϕ dx+ 
  	∫_Ω_2ρ_2∇(u_2^*-u)·∇ϕ dx = 
  	∫_Ω_o21(ρ_1-ρ_2)∇u_2^*·∇ϕ dx.
  	
  
  Applying integration by parts on the right 
  hand side in (<ref>) and then setting ϕ=u_2^*-u, we derive that

  	
  	∫_Ω_2^*ρ   |∇(u_2^*-u)|^2 dx = 
  	c_ρ(-∫_Ω_o21ρ_2Δu_2^*(u_2^*-u)  dx +
  	∫_F_o12ρ_2 ∇u_2^*·n_F_o12(u_2^*-u) dσ) 

  	  	≤   
  	c_ρ(∫_Ω_o21f (u_2^*-u)  dx +
  	∫_F_o12ρ_2 ∇u_2^*·n_F_o12(u_2^*-u) dσ) 

  	  	 (<ref>)≤    
  	 	c_ρf_L^2(Ω_o21)u_2^*-u_L^2(Ω_o21) +
  	 	 ρ_2 ∇u_2^*_L^2(F_o12)u_2^*-u_L^2(F_o12) 

          (<ref>)≤    
  	 	c_ρ f_L^2(Ω_o21) u_2^*-u_L^2(Ω_o21) +
  	 	 ρ_2 ∇u_2^*_L^2(F_o12)
  	 	 u_2^*-u^1/2_L^2(Ω_o21)u_2^*-u^1/2_H^1(Ω_o21)

        (<ref>)≤    
  	 	c_1( f_L^2(Ω_o21) d_o∇(u_2^*-u)_L^2(Ω_o21) 

  	 	   + ρ_2 ∇u_2^*_L^2(F_o12) 
  	 	d_o^1/2∇(u_2^*-u)^1/2_L^2(Ω_o21)
  	 	(d_o+1)∇(u_2^*-u)^1/2_L^2(Ω_o21)
   
  	 ≤   c_2 (f_L^2(Ω_o21)+	ρ_2 ∇u_2^*_L^2(F_o12)) 
  	 d_o^1/2∇(u_2^*-u)_L^2(Ω_o21),
  	
  
 where we have used that 0<d_o<1. By (<ref>), we can easily obtain that 
 
    ρ∇ (u_2^*-u)_L^2(Ω_2^*)≤ c_2 d_o^1/2(f_L^2(Ω_o21)+	ρ_2 ∇ u_2^*_L^2(F_o12)),
 
 and this gives an estimate of the difference between the physical solution u and the perturbed solution u^*.
 


 §.§ Proof of the  interpolation estimate (<ref>) 


Note that by Assumption <ref> and the definition of (PV1) we can conclude that 
Ω_1=Ω_1^* and u|_Ω^*_1=u^*_Ω^*_1. Hence we can construct an interpolant 
Π_1,h^*u such that 

    (∇(u_1-Π_1,h^*u_1)^2_L^2(Ω^*_1) + 
    	h∇(u_1-Π_1,h^*u_1)^2_L^2(F_o12)  
      + 
    	1/h(u_1-Π_1,h^*u_1)^2_L^2(F_o21))^1/2≤ 
    	C_1 h^min(ℓ-1,p)u_H^ℓ(Ω_1^*).

Next we show an interpolation estimate for u on Ω_2^*. Let us  denote D_1=Ω_o21 and D_2=Ω_2. 
Let the extension operator E_i:H^ℓ(D_i)→ H^ℓ(Ω^*_i), i=1,2, 
such that for each v∈ H^ℓ(D_i) it holds
(i) (E_iv)|_D_i = v and (ii) E_i v_H^ℓ(Ω_i^*)≤ C_E_iv_H^ℓ(D_i), where the constant 
C_E_i depending only on D_i and Ω_i^*, see <cit.>. 
We recall  the B-spline interpolation operatorΠ_h^* given in (<ref>) and define 

    Π_h^* ṽ :=(Π_1,h^*ṽ_1, Π_2,h^* ṽ_2),
    	  where Π_i,h^*  ṽ_i:=Π_i,h^* E_i v.

Recalling u_i=u|_Ω_i and using  the properties of the extension operator and (<ref>) we have


    ∇(u_1 - Π_1,h^* ũ_1)_L^2(Ω_o21)≤∇(E_1u - Π_1,h^* ũ_1)_L^2(Ω_1^*)≤ C_intp C_E_1  h^s u_1_H^ℓ(Ω^*_1),

and

    ∇(u_2 - Π_2,h^* ũ_2)_L^2(Ω_2)≤∇(E_2 u - Π_2,h^* ũ_2)_L^2(Ω_2^*)≤ C_intp C_E_2  h^s u_2_Ω_2,


where s=min(p, ℓ-1). 

Using the trace inequality, <cit.>,  
v^2_L^2(F_21)≤ C (h^-1v^2_L^2(Ω_o21) +h|∇ v^2_L^2(Ω_o21))
and proceeding as in (<ref>) we can show


    (h∇(u_1-Π_1,h^*ũ_1)^2_L^2(F_o21))^1/2≤ C_intp C_E_1  h^s u_1_H^ℓ(Ω^*_1), 
    (1/h(u_1-Π_1,h^*ũ_1)^2_L^2(F_o21))^1/2≤
    C_intp C_E_1   h^s u_1_H^ℓ(Ω^*_1).


Gathering the inequalities (<ref>), (<ref>) and (<ref>)  we can derive  (<ref>). 









plain



































