title: A general framework for substructuring-based domain decomposition methods for models having nonlocal interactions
authors: Giacomo Capodaglio Marta D'Elia Max Gunzburger Pavel Bochev Manuel Klar Christian Vollmann
text: 
  Accurate Deep Potential model for the Al-Cu-Mg alloy in the full concentration space The work of H. W.
  is supported by the National Science Foundation of China under Grant No. 11871110, the National Key Research and Development Program of China under Grants No. 2016YFB0201200 and No. 2016YFB0201203, and Beijing Academy of Artificial Intelligence (BAAI). The work of L. Z. was supported in part by  the Center of Chemistry in Solution and at Interfaces (CSI) funded by the DOE Award de-sc0019394.
      Wanrun Jiang^1,2,  Yuzhi Zhang^3,4, 
   Linfeng Zhang^5,Corresponding author. E-mail: linfengz@princeton.edu ,  and  Han Wang^6,Corresponding author. E-mail: wang_han@iapcm.ac.cn 
  
  ^1Songshan Lake Materials Laboratory,
  Dongguan, Guangdong 523808, People's Republic of China
   ^2Institute of Physics, Chinese Academy of Sciences,
   Beijing 100190, People's Republic of China
    
  ^3Beijing Institute of Big Data Research,
   Beijing 100871, People's Republic of China
  
  ^4Yuanpei College of Peking University,
  
  Beijing 100871, People's Republic of China
  
   
  ^5Program in Applied and Computational Mathematics,
  
  Princeton University, Princeton, NJ 08544, USA
  
  ^6Laboratory of Computational Physics,
  
  Institute of Applied Physics and Computational Mathematics,
  
  Huayuan Road 6, Beijing 100088, People's Republic of China
  
  
      December 30, 2023
  =============================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================
  
  
  
  
  A rigorous mathematical framework is provided for a substructuring-based domain-decomposition approach for nonlocal problems that feature  interactions between points separated by a finite distance. Here, by  substructuring it is meant that a traditional geometric configuration for local partial differential equation problems is used in which a computational domain is subdivided into non-overlapping subdomains. In the nonlocal setting, this approach is substructuring-based in the sense that those subdomains interact with neighboring domains over interface regions having finite volume, in contrast to the local PDE setting in which  interfaces are lower dimensional manifolds separating abutting subdomains. Key results include the equivalence between the global, single-domain nonlocal problem and its multi-domain reformulation, both at the continuous and discrete levels. These results provide the rigorous foundation necessary for the development of efficient solution strategies for nonlocal domain-decomposition methods.
   
  
  
  Nonlocal models, domain decomposition, finite element methods
  
  
  
  34B10, 65M60, 45P05, 45A99, 65R99
  
  
  
  myheadings
  plain
  G. Capodaglio, M. D'Elia, M. Gunzburger, P. Bochev, M. Klar and C. Vollman
  A General Framework for Nonlocal Domain Decomposition
  
  
  
  
  
  
  § INTRODUCTION
  
  
  Nonlocal models have become a popular alternative to partial differential equation (PDE) models due to their ability to describe effects that PDEs fail to capture. In particular, a nonlocal model can describe multiscale and anomalous behavior for applications that exhibit hierarchical features that cannot be reproduced by a classical model. These applications include, among others, subsurface transport <cit.>,
  image processing <cit.>,
  multiscale and multiphysics systems <cit.>,
  magnetohydrodynamic <cit.>,
  finance <cit.>, and stochastic processes <cit.>.
  
  The general class of nonlocal models we consider are characterized by integral operators having the form
  
  
      u()= ∫_B_δ() (u()-u())(,) d,
  
  
  where B_δ() denotes the ball (usually Euclidean) centered at  with radius δ (usually referred to as the horizon or interaction radius) and (,) is an application-dependent kernel function (usually symmetric in its arguments and nonnegative) that determines the regularity properties of the solution. The nonlocality inherent in (<ref>) is clear: points  interact with all points  located within a distance δ. Compared to that for the local PDE setting, the integral form clearly reduces regularity requirements on the solution and allows for the capture of long-range interactions. 
  
  However, the utilization of nonlocal models in applications that could benefit from their improved predictive capabilities is hindered by several modeling and numerical challenges. These include the unresolved treatment of nonlocal interfaces <cit.>, the nontrivial prescription of nonlocal volume constraints (the nonlocal counterpart of boundary conditions) <cit.>, and the fact that computational costs attendant to the use of nonlocal problems may become prohibitive as the extent of the nonlocal interactions increases; see, e.g., <cit.> for variational methods and <cit.> for mesh-free methods. Other critical challenges are related to the uncertain nature of model parameters; in fact, modeling parameters such as δ and those characterizing the kernel, applied forces, and/or sources can be non-measurable, sparse, and/or subject to noise. Research on such topics is very active (see, e.g., <cit.>) but further consideration of them is beyond the scope of this work. 
  
  Here, we focus on the treatment of nonlocal interfaces and, notably on the design of nonlocal domain-decomposition (DD) formulations with the aim of reducing computational costs by increasing the parallel concurrency in the numerical solution of nonlocal problems. Specifically, the goal is to address the high computational cost associated with nonlocal models by providing a foundational algorithmic framework for their parallel solution, mirroring that of successful parallel DD algorithms for PDEs such as, e.g., Finite Element Tearing and Interconnecting (FETI) <cit.> and other approaches <cit.>.
  
  This work is part of a comprehensive effort by the authors to fill the theoretical and practical gaps in the current understanding of nonlocal interfaces (both physical ones and those created by DD solution algorithms) by developing a rigorous nonlocal interface theory for nonlocal diffusion (see the preliminary work <cit.>), including pure fractional diffusion, and nonlocal mechanics. Our ultimate goal is to design efficient and scalable DD solvers to unlock the full potential of nonlocal models. To this end, as is often done in nonlocal modeling, we draw inspiration from the existing work on classical DD methods for PDEs. 
  
  Notwithstanding the vast knowledge accumulated by the DD community, extension of local DD methods to nonlocal models remains a nontrivial task: nonlocality introduces many challenges and limitations.
  
  Some of these challenges are shared with local DD methods. An example of a shared challenge is the proper treatment of floating subdomains[In a multi-domain setting, floating domains are subdomains of the decomposition that are either internal (they do not share boundaries with the physical boundary) or share boundaries only with parts of the physical boundary at which Neumann-type conditions are prescribed. In this context “floating” refers to the fact that these domains either do not have volume constraints because they are internal, or have partial constraints of Neumann type.] which require special attention due to the singularity of the discretized equations on such subdomains.
  
  However, other challenges are unique to the nonlocal setting and require new approaches that have no analogues in the local setting. 
  For example, in the local PDE setting, many DD methods typically start by breaking up the computational domain into non-overlapping subdomains, a process we refer to as substructuring of the domain. These subdomains interact only through their shared boundaries on which one usually imposes some appropriate continuity conditions. 
  
  Although one starts from the same initial geometric configuration in which the domain is substructured into non-overlapping subdomains, inherent nonlocal interactions between the subdomains force one to expand these subdomains to include parts of neighboring subdomains having nonzero volume, causing an overlap of a thickness determined by the interaction radius δ. This overlap is required because, in the nonlocal setting, it is not possible to define subdomain problems simply by restricting the global operator to the subdomains. For this reason we refer to our approach as being substructuring-based. 
  
  It is important to point out that the thickness of the overlap regions depends solely on the modeling parameter δ and is unrelated to the discretization method employed and, in particular, to the grid size. 
  As a result, discretization of the decomposed nonlocal problem requires special care because the overlapping regions induced by the decomposition do not, in general, match the underlying mesh.
  
  The current state of the art of nonlocal DD methods is very limited, with <cit.> being perhaps the most relevant work. In that paper, the authors consider a simple two-domain configuration and develop a variational approach to DD based on adding an interface equation and a new variable that lives on the overlap between the subdomains, while using test functions that vanish on the interface for each subdomain. Decoupling is achieved by solving a Schur-complement problem for the interface variable, similar to a conventional FETI scheme. The subproblem definition in <cit.> does not consider  multi-domain configurations nor does it consider floating subdomains. 
  
  As a result, extension of this approach to the general multi-domain case, which is a prerequisite for efficient parallel nonlocal DD algorithms, remains an open question.
  
  In this paper we formulate a general framework for nonlocal DD problems that has the following equivalence property:
  
  
      
   
  
  While the significance of this property is self-evident, its fulfillment is far from trivial, and is one of the key challenges addressed in the paper. 
  
  As already mentioned, we refer to our approach as “substructuring-based” because, much like as it is in standard non-overlapping DD, the subdomains interact only through their shared interfaces. Of course, the key difference is that in the nonlocal setting these shared interfaces are nonlocal, i.e., they  are  regions having finite volumes as opposed to the local case in which interfaces are lower-dimensional manifolds.
  
  The main contributions of this paper are as follows.
  
      
    * We introduce a systematic way to decompose the domain given an existing mesh for the single-domain problem and we discuss ways to make the decomposed domains compatible with the given mesh. Specifically, we provide a recipe for decomposing the domain that prevents integration over partial (cut) finite elements by using approximate neighborhoods in a manner such that the equivalence of the decomposed solution and the single-problem one is not compromised.
      
    * We formulate a continuous nonlocal DD system of subdomain problems and prove that it is equivalent to the single-domain problem, i.e., we provide a solution for (<ref>). The key ingredient is the appropriate definition of indicator functions that keep track of the number of overlapping subdomains. 
      
    * We define a Galerkin finite element discretization of the nonlocal DD problem and prove that it is equivalent to the discretization of the single-domain domain problem effected using the same type of finite element functions. This equivalence holds for the finite-dimensional variational formulation and for the corresponding matrix form.
  
  
  Our nonlocal DD formulation provides a mathematical foundation for the development of a range of efficient numerical algorithms for the parallel solution of nonlocal problems that mirror existing approaches for local problems. For example, treating the nonlocal interface equations as constraints and using Lagrange multipliers to enforce them lends itself to the development of nonlocal FETI <cit.> or Arlequin-like <cit.>  algorithms.  Alternatively, one can choose to view these coupling conditions as an optimization objective and treat the subdomain equations as constraints. Such an approach would lead to nonlocal optimization-based DD methods that are nonlocal counterparts of the methods in <cit.>.
  
  It should be noted though that realizing the potential of our DD framework to reduce the computational burden of solving nonlocal problems requires the ratio between δ and the diameter of the subdomains to be smaller than 1, i.e., we target problems for which the extent of the nonlocal interactions is much smaller than the diameter of the domain. Such problems arise in several engineering applications such as, e.g., nonlocal mechanics, and are the main motivator for this work. In contrast, for applications described by nonlocal operators with infinite interactions, a DD approach may not be as effective because  the interaction regions would span a portion of the domain that is of the same size (or even larger) than the domain itself.
  
  Finally, we mention that often one may be given a decomposition of Ω into a few subdomains which are constructed to follow well-defined geometric entities, e.g., a wing and a fuselage, or different media properties, e.g., different diffusion coefficients, within Ω. Such “physically”-motivated DDs typically arise in the context of mesh tying <cit.> in which a complex geometric entity is broken into smaller parts to enable efficient mesh generation. In contrast, here we focus on DD as a means for faster and more efficient parallel solution methods for nonlocal problems in which case the number of the subdomains is very large and they do not generally follow any “physics”-motivated interfaces. 
  
  
  
  
  
    
  Outline of the paper The paper is organized as follows. In Section <ref>, we recall the variational formulation of a single-domain volume-constrained nonlocal Poisson problem and briefly describe its discretization via finite element methods. In Section <ref>, we introduce the continuous formulation of a multi-domain nonlocal DD method and prove its equivalence to the single-domain problem presented in Section <ref>. Section <ref> explains how we address the decomposition problem by formulating rules for the construction of the subdomains and their interaction regions that fulfill (<ref>), i.e.,  guarantee the equivalence of the DD and the single-domain problems. In the same section we also introduce the discretized subproblems and their matrix forms and show their equivalence  to the underlying single-domain formulation. Concluding remarks are provided in Section <ref>. 
  
  In the technical report  <cit.>, we further elucidate the equivalence between the multi-domain formulation and the single-domain problem and illustrate the application of the nonlocal DD framework.  Specifically, we use the framework developed in this work along with a FETI solution approach to obtain, for a very simplified setting, illustrative numerical examples of nonlocal DD problems.
  
  
  
  
  § A NONLOCAL (SINGLE-DOMAIN) VOLUME-CONSTRAINED PROBLEM AND ITS FINITE ELEMENT DISCRETIZATION
  
  
  For simplicity, in this work, we consider the two-dimensional case. Let Ω denote a bounded, open subset of ^2. For any δ > 0, often referred to as the horizon or interaction radius, we define the associated interaction domain as the closed region
  
      Γ_interaction= {∈^2  ∖Ω : ∃ ∈Ω  such that   | - | ≤δ } .
  
  
  The interaction domain is split into two disjoint parts Γ and Γ_Neumann, where Γ is a nonempty closed domain, whereas Γ_Neumann is allowed to be empty. Thus, we have that Γ∪Γ_Neumann=Γ_interaction and Γ∩Γ_Neumann=∅, where Γ_Neumann is open along its common boundary with Γ. Also, note that Γ_interaction and therefore also Γ and Γ_Neumann depend on δ, even though that dependence is not explicitly indicated. Figure <ref>-left illustrates this geometric configuration.[Domains such as Γ and Γ_Neumann in Figure <ref>-left and therefore also in subsequent figures are stylized versions of their true shapes. For example, because points  interact only with points ∈ B_δ(), where B_δ() denotes the Euclidean ball of radius δ centered at , those domains have rounded corners. However, in practice, one can keep the stylized domains because the points outside the true interaction domains are not accessed during a finite element assembly process.] 
  
  
  
  
  The strong formulation of a nonlocal volume-constrained Poisson problem is given by <cit.>[The problem (<ref>) is a nonlocal analogue of the PDE Poisson problem -∇·(κ∇ u)=f_Ω in Ω, u=g on a nonempty part of the boundary of Ω, and κ∇ u· n= f_Neumann on the rest of that boundary.]
  
  
      -2∫_Ω∪Γ∪Γ_Neumann(u() -u())γ(,)d 
      = f_Ω()    ∈Ω          (a)
      
      
      u() = g()    ∈Γ           (b)
      
      
       2∫_Ω∪Γ∪Γ_Neumann(u() -u())γ(,)d 
      = f_Neumann()    ∈Γ_Neumann,    (c)
  
  
  where f_Ω(), f_Neumann(), and g() are given functions and γ(,) is a  given symmetric nonnegative kernel, i.e., γ(,) = γ(,) for , ∈Ω∪Γ∪Γ_Neumann.
  
  Equations (<ref>b) and (<ref>c) are nonlocal analogues of Dirichlet and Neumann boundary conditions, respectively, for PDEs; specifically, equation (<ref>b) is a Dirichlet volume constraint imposed on a set Γ having nonzero measure in ^2 and, if Γ_Neumann also has nonzero measure in ^2, (<ref>c) is a Neumann volume constraint on that set.
  
  The nonlocal operators in (<ref>a) and (<ref>c) are identical up to a sign so that these equations can be combined to obtain an equivalent, more compact,  strong form
  
  
      -2∫_Ω∪Γ(u() -u())γ(,)d 
      = f()    ∈Ω
      
      u() = g()    ∈Γ,
  
  
  where Ω=Ω∪Γ_Neumann, f()|_Ω=f_Ω(), and f()|_Γ_Neumann=-f_Neumann(). The strong form (<ref>) corresponds to the configuration in Figure <ref>-right and is used in the remainder of the paper.
  
  We assume that Ω does not include the boundary portion ∂Ω∩∂Γ, but does include the boundary portion ∂Ω∖(∂Ω∩∂Γ), where ∂Ω and ∂Γ denote the boundaries of Ω and Γ, respectively.
  
  We define the function spaces
  
      {    W= {w ∈ L^2(Ω∪Γ)   :   |||w||| < ∞}
                
        |||w|||^2 =  ∫_Ω∪Γ∫_Ω∪Γ
        |w() - w()|^2γ(,)d  d +
        w^2_L^2(Ω∪Γ)
          W^0 = {w∈ W  :   w=0   }.
  
  
  and, for u,v∈ W, we define the bilinear form 𝒜(·,·) and linear functional (·) as
  
  
      {   𝒜(u,v) = 
      ∫_Ω∪Γ∫_Ω∪Γ(u() - u())(v() - v())γ(,)d  d
         (v) = ∫_Ω v()f()d. 
      .
  
  
  Then, a weak formulation of the nonlocal volume-constrained problem  (<ref>) can be stated as <cit.>
  
  
      𝒜(u,v) = (v)   ∀  v ∈ W^0
              .
  
  Here, W' denotes the dual space of bounded linear functionals on W^0 with respect to the standard L^2 duality pairing and W_Γ denotes the nonlocal “trace” space defined as W_Γ={w|_Γ  :   w∈ W}.
  If Γ has nonzero measure in ^2, the well posedness of the problem (<ref>) is proved in, e.g., <cit.>.   
  
  
  For v∈ W, define the energy functional
  
      ℰ_single(v) = 1/2𝒜(v,v) - (v).
  
  Then, (<ref>) is equivalent to the  minimization problem <cit.>
  
      ℰ_single(u) = inf_v∈ Wℰ_single(v) 
               .
  
  
  
  The functional setting and the well posedness of the nonlocal problem depend on the kernel  γ(,). For example, if the kernel is square integrable (i.e., ∫_Ω∪Γ(γ(,))^2 d<∞ for all ∈Ω∪Γ) or if the kernel is integrable and translationally invariant (i.e., ∫_Ω∪Γγ(,) d<∞ for all ∈Ω∪Γ and γ(,))=γ(-)), it is known that W=L^2(Ω∪Γ); see, e.g., <cit.>. On the other hand, for the fractional kernel γ(,)∝ |-|^-d-2s with d denoting the space dimension and 0<s<1, it is known that W=H^s(Ω∪Γ), i.e., a fractional Sobolev space; again, see, e.g., <cit.>.
  
  However, the technical aspects of this work are largely independent of the choice of the kernel as long as the nonlocal problem remains well posed. Moreover, some conditions such as the symmetry of  γ(,) can be further relaxed under some additional assumptions <cit.> that ensure the well posedness of the nonlocal problem. Likewise, one can also relax the condition that γ(,) is nonnegative everywhere; see, e.g., <cit.>. Although for such kernels, an energy minimization characterization of the nonlocal problem may not be available, 
  well-defined weak formulations of the strong form nonlocal equations still exist for more general cases. As a result, the algorithms developed in this work can be extended in a straightforward manner to problems that cannot be characterized in terms of an energy minimization setting.
      
  
  
  
  
  
  
   §.§ Finite element discretization of the nonlocal volume-constrained problem
  
  The weak formulation (<ref>) of the nonlocal volume-constrained problem can be discretized using a finite element method as follows. Let 𝒯^h denote a finite element triangulation of Ω∪Γ parameterized by a grid size parameter h. We assume that 𝒯^h conforms to the boundary of Ω, i.e., ∂Ω consists of finite element edges. This requirement can be satisfied by first constructing a grid in Ω after which a grid is constructed in Γ that shares element vertices with those of the grid in Ω along their common boundary. For simplicity, in the sequel, we restrict the discussion to Lagrangian finite element spaces.
  
  Let W^h⊂ W and W^0,h⊂ W^0 denote finite element (FE) subspaces. Then, a FE approximation u^h()∈ W^h of the solution u∈ W of (<ref>) is defined to be the solution of the discretized weak formulation
  
  
      , 
      
         
             𝒜(u^h,v^h) = (v^h)   ∀  v^h ∈ W^0,h    ,
  
  where g^h() denotes an approximation of g() that is usually chosen to be the FE interpolant[When g() is not of class C^0, g^h() can be defined by, e.g., least-squares approximation or Clement interpolation <cit.>.] of  g(). As long as Γ has nonzero measure in ^2, the well posedness of problem (<ref>) is guaranteed; see, e.g., <cit.>. 
  
  Let N^h denote the number of degrees of freedom corresponding to the nodes in Ω∪Γ and let N_h denote the number of degrees of freedom corresponding to the (possibly semi-) open domain Ω, with the remaining and N^h-N^h degrees of freedom corresponding to the closed domain Γ.
  Note that because Ω is an open domain with respect to its common boundary with Γ and Γ itself is a closed domain, nodes and degrees of freedom along their common boundary are assigned to Γ. We then define the finite element subspace W^h⊂ W as the span of a nodal finite element basis {ϕ_i()}_i=1^N^h so that a finite element approximation u^h() of the solution u() of (<ref>) can be expressed as
  
  
      u^h() = ∑_i=1^N^h_i ϕ_i()
      + ∑_i=N^h+1^N^h_i ϕ_i(),
  
  where  denotes an N^h-vector of unknown coefficients and  denotes an (N^h-N^h)-vector of nodal values of the approximation g^h() of g(). 
  
  We recall that the support of each nodal basis function ϕ_i() comprises all elements sharing the node _i. As a result, all basis functions corresponding to nodes in Ω vanish on Γ and  span{ϕ_i()}_i=1^N^h⊂ W^0,h. 
  
  Let 𝔸_single  and b_single denote the N^h× N^h matrix and the N^h-vector with elements
  
  
      {
      (𝔸_single)_ij = 𝒜(ϕ_j,ϕ_i)
                                    i,j=1,…,N^h
      
      
      (b_single)_i = (ϕ_i)
      -
      ∑_j=N^h+1^N^h 𝒜(ϕ_j,ϕ_i) _j
               
      i=1,…,N^h,
      .
  
  respectively. Then, the discrete FE problem (<ref>) is equivalent to the linear algebraic system
  
  
      𝔸_single = _single
  
  
  for the unknown nodal coefficient vector u. The matrix 𝔸_single is symmetric and, owing to the fact that Γ has nonzero measure in ^2, it is also positive definite <cit.>.
  
  
  For v^h∈ W^h, we define the discrete energy functional  
  
      ℰ_single^h(v^h) = 1/2𝒜(v^h,v^h) - (v^h).
  
  Then, the discrete nonlocal volume-constrained problem (<ref>), respectively (<ref>), can be expressed in terms of the equivalent minimization problem <cit.>
  
      ℰ^h_single(u^h()) = min_v^h()∈ W^hℰ^h_single( v^h())       v^h()=g^h() .
  
  
  
  As implied by (<ref>), any point ∈Ω∪Γ interacts only with points in the ball B_δ(). This raises a serious issue in FE methods for nonlocal problems because the intersection of these balls with the finite element grid produces cut elements, i.e., partial elements, within the ball. As a result, one either has to deal with cut elements or, if one wants to only deal with uncut elements, one is faced with discontinuous integrands that vanish outside the ball. This issue is glossed over in many FE papers for nonlocal problems, especially those that only provide one-dimensional numerical results. However, we address this issue in Section <ref>. A comprehensive discussion of how to effectively handle cut elements can be found in <cit.>.
  
      
  
  
  
  
  
  
  
  § NONLOCAL DOMAIN DECOMPOSITION IN THE CONTINUOUS SETTING
  
  
  In this section we first describe, in the continuous setting, how to define a nonlocal decomposition of the domain and then introduce the formulation of the multi-domain system. The central result of this section proves the equivalence of the solution of the single-domain system and the one corresponding to the multi-domain system. The significance of this result is that it establishes the consistency and the well posedness of our multi-domain formulation.
  
  
  
  
   §.§ Construction of the geometric domain decomposition
  
  
  In a standard PDE domain decomposition setting, one can partition Ω into non-overlapping subdomains and then simply define the subdomain problems by restricting the global operator to each subdomain. Such a construction is impossible in the nonlocal setting due to the inherent nonlocal interactions which require any two adjacent subdomains to share an interface having nonzero volume. As a result, our substructuring-based domain decomposition starts from a non-overlapping, covering subdivision of Ω into N_s subdomains {Ω_n}_n=1^N_s, as illustrated in Figure <ref>-left for N_s=6, and then adds the overlaps necessary for the nonlocal interactions. 
  
  Note that some of the domains Ω_n include part of the boundary ∂Ω. For example, in Figure <ref>-left, we have that this is the case for Ω_2, Ω_3, and Ω_5.
  
  
  
  For each subdomain Ω_n, n=1,…,N_s, we define the (possibly semi-) open smaller subdomain 
  
      Ω_n =
         {∈Ω_n   
         :    |-|>δ/2  ∀ ∈Ω∖Ω_n  
      };
  
  see Figure <ref>-right for an illustration. Note that {Ω_n}_n=1^N_s is, by construction, a set of non-overlapping domains. We also subdivide the interaction domain Γ into a set of overlapping, covering subdomains 
  
  
      Γ_n = {∈Γ 
          :  ∃ ∈Ω_n   such that  |-| ≤δ }
  
  and also define the set of overlapping subdomains
  
      Γ_n = {∈Ω∖Ω_n 
          :  ∃ ∈Ω_n   such that  |-| ≤δ }    .
  
  
  In Figure <ref>-right, for each n, Γ_n∪Γ_n consists of all the strips of thickness δ that surround Ω_n, including in some instances a portion of Γ. Because each subdomain Γ_n overlaps with at least one other subdomain Γ_n', n' n, and some subdomain Γ_n may overlap with another subdomain Γ_n', n'≠ n, the above construction results in the overlapping domain decomposition of Ω∪Γ given by
  
      Ω∪Γ = ∪_n=1^N_sΩ_n∪Γ_n∪Γ_n.
  
  
  Following conventional DD nomenclature, we subdivide the set of subdomains {Ω_n∪Γ_n∪Γ}_n=1^N_s into two classes:
  
      
   
  For example, in Figure <ref>-right, Ω_3 and Ω_4 are floating subdomains, whereas Ω_1, Ω_2, Ω_5, and Ω_6 are non-floating. 
  
  Analogous to the conventional local DD setting, a floating domain is endowed with a purely Neumann nonlocal volume constraint so that its associated nonlocal problem has a non-trivial null space. In a typical local or nonlocal DD configuration, the number of subdomains is large and most of them are of the floating type. 
    
  
  
  As alluded to in Section <ref>, in a local PDE setting one can consider both non-overlapping and overlapping DD algorithms because the latter offer some computational conveniences and may allow for a faster convergence of iterative solution methods. Typically, the overlap depends on the grid size and its thickness goes to zero as the mesh size is reduced. In contrast, the nonlocal setting requires any two adjacent subdomains to overlap in order to compute the necessary nonlocal interactions between them. As a result, the size of this mandatory overlap is determined not by the mesh size but by the interaction radius δ, i.e., its thickness is independent of the underlying discretization mesh.
      
  
  
  
  
  
   §.§ The domain decomposition (multi-domain) system
  
  
  A multi-domain system is a system of N_s equations, each of which holds for ∈Ω_n∪Γ_n∪Γ_n, n=1,…,N_s. In constructing the multi-domain system we have to deal with overlapping domains, i.e., although the three domains ∪_n=1^N_sΩ_n, ∪_n=1^N_sΓ_n, and ∪_n=1^N_sΓ_n are mutually disjoint as are the N_s domains Ω_n, there are overlaps among the N_s domains in ∪_n=1^N_sΓ_n and there could be overlaps among the N_s domains in ∪_n=1^N_sΓ_n.
  
  
  
  
  To properly deal with the consequences of having overlapping domains, we define the generalized characteristic  functions
  
      ζ_(,)=
      ∑_n=1^N_s_Ω_n∪Γ_n∪Γ_n()
      _Ω_n∪Γ_n∪Γ_n()   
          ζ_()= ∑_n=1^N_s_Ω_n∪Γ_n().
  
  Note that ζ_(,) is a symmetric function, i.e., ζ_(,)=ζ_(,) and is a non-negative piecewise integer-valued function and  ζ_() is a positive piecewise integer-valued function. 
  
  
  It is possible for ζ_(,)=0. For example, this is the case if ∈Ω_n and  ∈Ω_n' with n' n. However, this does not present a problem for (<ref>) below, which involves ζ_(,)^-1, because points  interact only with points  such that |-|≤δ and for such pairs of points, ζ_(,)>0.
  
      
  
  
  
  
  
    §.§.§ The subdomain system
  
  
  For n=1,…,N_s and any pair of functions u_n() and v_n() defined on Ω_n ∪Γ_n∪Γ_n, 
  we define the subdomain bilinear form
  
      𝒜_n(u_n,v_n)
      
         
      =∫_Ω_n∪Γ_n∪Γ_n^  ∫_Ω_n∪Γ_n∪Γ_n^ζ_(,)^-1(u_n() - u_n())(v_n() - v_n())γ(,)d  d
  
  and the associated subdomain linear functional
  
      _n(v_n) = ∫_Ω_n∪Γ_nζ_()^-1v_n()f()d,
  
  where, of course, for floating domains, i.e., if Γ_n=∅, the integrals over Γ_n vanish.
  
  For n=1,…,N_s, we define the function spaces
  
      {   W_n = {w∈ L^2 (Ω_n∪Γ_n∪Γ_n)   :   |||w|||_n < ∞}
         W_n^0 = {w∈ W_n  :   w=0   },
      .
      
                            |||w|||_n^2 =  𝒜_n(w,w) + w^2_L^2(Ω_n∪Γ_n∪Γ_n).
  
  Let
  
      W_n = {   W_n^0   ,     
         W_n   ,    .
  
  We then define the domain-decomposition or multi-domain system of equations as
  
      {   
         
              𝒜_n(u_n,v_n) = _n(v_n) 
        ∀  v_n∈W_n
      (a)
      
         
           
       u_n() = u_n'()     (b)
      
         
            u_n() = g()     ∀ ∈Γ_n  ,   .
           (c)
      .
  
  
  
  Equation (<ref>b) can be thought of as a nonlocal version of the standard continuity constraint in non-overlapping local DD methods. 
  The constraints in that equation are needed
  because the solution u() of (<ref>) is a single-valued function on Ω∪Γ, and in particular on ∪_n=1^N_sΓ_n. However, by construction, we have that for all n' such that Γ_n∩Γ_n'∅, both u_n() and u_n'() are defined on Γ_n∩Γ_n'∅. Clearly, we have that (<ref>b) must be imposed on that domain. 
  
  Of course, this equation automatically holds on Γ_n∩Γ_n' because both u_n()=g() and u_n'()=g() there.
  
  
  As in the standard (local) DD case, the constraints in (<ref>b) ensure single-valued solutions of the multi-domain system, and are appropriate  when the global nonlocal solution is also continuous. However, one of the principal advantages of nonlocal models is that, for some kernels in common use, they admit solutions with jump discontinuities  <cit.>. To allow for such solutions, one may choose to enforce (<ref>b) weakly, i.e., for n=1,…,N_s and n'=n+1,…,N_s,
  
  
      ∫_Γ_n∩Γ_n'( u_n() - u_n'() )v() d = 0
          ∀  v∈ W|_Γ_n∩Γ_n',
  
  
  Note that (<ref>) also arises when (<ref>b) is enforced using Lagrange multipliers; such a treatment of (<ref>b) would result in finite element formulations of the constraints (<ref>b) involving mass matrices.
  
      
  
  
  
  The constraints in (<ref>b) are not independent. For example, consider a point 
  ∈Γ_1∩Γ_4∩Γ_6 near the bottom right corner  of Ω_1 in Figure <ref>-right. Then, (<ref>b) would include the constraints u_1()= u_4(), u_1()= u_6(), and u_4()= u_6(), only two of which are independent. These redundancies in (<ref>b) have implications in the design of discretization algorithms as is discussed in Section <ref>.
  
      
  
  
  
  
  
  
    §.§.§ Equivalence of the single-domain and multi-domain problems
  
  
  Our next task is to show that the solution u() for ∈Ω∪Γ of the single-domain system (<ref>) and the solutions u_n(), n=1,…,N_s, of the multi-domain system (<ref>) are the same or, more precisely, that u_n()=u() for ∈Ω_n∪Γ_n∪Γ_n. 
  The first step towards that end is the following lemma. We refer to Appendix <ref> for a proof.
  
  
  Given functions u() and v() for ∈Ω∪Γ, define the functions u_n() and v_n() for ∈Ω_n∪Γ_n∪Γ_n as
  
      u_n() = u()|_Ω_n∪Γ_n∪Γ_n    
        v_n() = v()|_Ω_n∪Γ_n∪Γ_n  .
  
  Then, for n=1,…,N_s and n'=n+1,…,N_s, we have that, if Γ_n∩Γ_n'∅ and Γ_n∩Γ_n'∅,
  
      u_n'() = u_n() 
                    v_n'() = v_n() 
                 ∈Γ_n∩Γ_n'        ∈Γ_n∩Γ_n'.
  
  
  Then, for the bilinear forms and linear functionals defined in (<ref>), (<ref>), and (<ref>), we have that
  
      𝒜(u,v) =∑_i=1^N_s𝒜_n(u_n,v_n) 
              (v) =  ∑_n=1^N_s_n(v_n).
  
  
  
  
  Using Lemma <ref> and the fact that W = ⊗_n=1^N_s W_n and W^0 = ⊗_n=1^N_sW_n one can easily prove the following equivalence result.
  
  
  
  Let u()∈ W denote the solution of the single-domain system (<ref>) and, for n=1,…,N_s, let u_n()∈ W_n denote the solution of the n-th subproblem in multi-domain system (<ref>). Then, the multi-domain system (<ref>) and the single-domain system (<ref>) are equivalent and their respective solutions coincide, i.e, u_n() = u()|_Ω_n∪Γ_n∪Γ_n.  
  
  
  
  The equivalence of the multi-domain and the global problems and the fact that the latter is well posed imply that (<ref>) is also well posed.
  
  We note that the equivalence between the multi-domain weak formulation (<ref>) and the single-domain system (<ref>) can also be established by noting that the single-domain energy functional (<ref>) can be written as a sum of subdomain energy functionals defining a multi-domain energy functional. The proof is straightforward and is omitted.
  
  
  
  Define the subdomain energy functionals ℰ_n(u_n) by
  
      ℰ_n(u_n):= 1/2𝒜_n(u_n,u_n) - _n(u_n)     ,
  
  where the bilinear form _n(·,·) and linear functional _n(·) are defined in (<ref>) and (<ref>), respectively. Then,
  
      ∑_n=1^N_sℰ_n(u_n) = ℰ_single(u).
  
  Furthermore, the multi-domain weak formulation (<ref>) is the Euler-Lagrange equation corresponding to the minimization problem
  
      inf_v_n∈ W_n,  n=1,…,N_s∑_n=1^N_sℰ_n(v_n)
           {   
         
         .
    
  
  
  
  
  
  
  
  § FINITE ELEMENT DISCRETIZATION OF THE SUBDOMAINS SYSTEMS
  
  
  In this section we describe a finite element discretization of the multi-domain system (<ref>). We follow a standard FE assembly procedure specialized to our needs. 
  
  Specifically,  for each n=1,…,N_s, we proceed as follows.
  
   1. Subdomain meshing. This step introduces a suitable finite element partition on  each subdomain Ω_n∪Γ_n∪Γ_n.
  
  2. Selection of a subdomain finite element space.
  We endow each subdomain Ω_n∪Γ_n∪Γ_n  with a finite element space W_n^h defined as the span of a set of basis functions, usually chosen to be piecewise polynomials with respect to the underlying grid.  
  
  3. Subdomain finite element ansatz. We seek an approximation of the solution u_n() of (<ref>) out of the finite element space W_n^h.
  
  4. Subdomain finite element assembly. We construct the stiffness matrix and right-hand side vector for each subdomain in the same manner as for the global problem, i.e., by inserting the finite element basis functions into the respective weak subdomain equations.
  
  As is the case with any DD method, the goal is to perform these steps in such a way so as to ensure that
  
  
      
  
  
  A standard way to achieve this goal in DD methods is to use a partition based on the underlying mesh. In this approach the subdomains are defined as unions of adjacent elements using a graph partitioning tool. The associated finite element spaces are then simply restrictions of the global space to the resulting subdomain. The first step, i.e., the subdomain meshing in this case is trivial because it uses an already defined global mesh. 
  
  This construction is not restrictive as we are interested in DD as means to devise fast parallel solution algorithms for a given global nonlocal problem. In such a case a global mesh and a finite element space can be assumed to have been already constructed. This should be contrasted with mortar methods and mesh tying where from the onset one assumes that each subdomain is meshed separately.
  
  
  
  
  
   §.§ Definition of subdomain grids
  
  
  We assume the setting of Section <ref>, i.e., that we are given a finite element mesh ^h of Ω∪Γ that conforms to the interface shared by Ω and Γ. As before, N^h denotes the total number of nodal degrees of freedom on this mesh and N^h are the nodal values associated with Ω. The remaining N^h-N^h degrees of freedom live on the closed domain Γ.
  
  In local PDE settings, for both non-overlapping and overlapping DD methods, the next step is to subdivide the domain Ω into subdomains that contain only whole finite elements. If one is going to invoke a non-overlapping DD method, the construction of the subdomains is complete. If instead an overlapping DD method is to be used, one adds, to each non-overlapping subdomain, whole elements in neighboring subdomains that are within a certain distance from the common boundary between the two domains; the distance used is usually related to some multiple of the local grid size, although other criteria are also in use <cit.>. 
   
  In the nonlocal case, the practical construction of a subdivision of Ω is similar to that for overlapping DD in the local case. We again start by subdividing Ω into the set {Ω_n}_n=1^N_s of subdomains with each Ω_n consisting of whole finite elements. We now want to add to and subtract from each subdomain Ω_n strips of thickness δ/2 to create the subdomains Ω_n and Γ_n. Unfortunately, because we are given a global grid over Ω∪Γ to work with, in general, we will not be able to (see Remark <ref>) define grids consisting of whole finite elements that respect the boundaries between Ω_n and Γ_n, i.e., those domains would also contain partial (cut) elements which is something we want to avoid. 
  
  We thus see that there is a big difference between the definitions of overlapping grids for the local and nonlocal cases. To recapitulate, in the local case we are free to add whole elements to effect an overlap with neighboring elements. In the nonlocal case, we do not have this freedom because the strips to be created have thickness δ irrespective of the given global finite element grid so that, in general, that strip will not consist of whole elements. Thus we have the choice of truncating elements so that the δ thickness of the strip is respected or instead approximate the strip by a strip consisting of whole elements which is tantamount to approximating the common boundaries of Ω_n and Γ_n by element edges. We use the latter choice because it is substantially easier to implement and, as shown below, does not compromise achieving the goal (<ref>).
  
  The above discussion motivates the following procedure for the construction, in the discretized setting, of a subdivision of Ω into subdomains that is analogous, but not the same, as that in Section <ref> for the continuous problem. We begin by assuming, as is done in Section <ref> for the single-domain setting, that 
  
   – we are given an integer N_s>1 and a finite element meshing ^h of Ω∪Γ which respects their common boundary Ω∩Γ. 
  
  We denote by ^h_Ω and ^h_Γ the sets of finite elements in Ω and Γ, respectively, and we denote by  a typical element in ^h. Then, 
  
   – we subdivide Ω into N_s non-overlapping, covering subdomains Ω_n, n=1,…,N_s, such that each subdomain Ω_n consists entirely of whole finite elements.
  
  This step is effected exactly in the same manner as for the local PDE non-overlapping DD setting so that no further comments are needed. 
  
  Note that the boundary ∂Ω_n of Ω_n consists of two or three disjoint, covering parts. First, we have for all n, 
  
   Type 1. ⋃_n'=1, n' n^N_s∂Ω_n∩∂Ω_n', i.e., the common boundary shared by Ω_n and subdomains Ω_n' that abut to Ω_n.
  
  We also have either one or both of
  
   Type 2. ∂Ω_n∩∂Γ, i.e., the common boundary shared by Ω_n and Γ.
  
   Type 3. ∂Ω_n∖[(⋃_n'=1, n' n^N_s∂Ω_n∩∂Ω_n')] ∪(∂Ω_n∩∂Γ), i.e., the part of ∂Ω_n that is not shared with the boundary of Γ or with any of the boundaries of other subdomains Ω_n'. 
  
  For example, referring to Figure <ref>-left, we have that the boundaries of floating domains such as Ω_3 consist of only Type 1 and 3 parts, the boundaries of domains such as Ω_6 consist of only Type 1 and 2 parts, and the boundaries of domains such as Ω_5 consists of all three parts.
  
  Recall that, by construction, the domains Ω_n consist of whole FE triangles. However, in general, the subdomains Ω_n, Γ_n, and Γ_n consist of whole FE triangles and additionally partial (cut) FE triangles. To obtain the equivalence of the single- and multi-domain FE solutions, we necessarily have to work with subdomains that consist of only whole FE triangles because only such triangles, i.e., triangles T∈𝒯_h, are used in the single-domain FE method. For this reason, we define an approach for the construction of subdomains Ω_n^h, Γ_n^h, and Γ_n^h of Ω∪Γ in such a way that all subdomains consists of only whole FE triangles. Additionally, the construction process is required to account for all interactions that occur between two subdomains Ω_n and Ω_n'. Meeting this requirement is guaranteed if all triangles that overlap with the Γ_n are included in Γ_n^h.
  
  
  The specific geometric domain decomposition we use is defined as follows.
  We denote by ^vertex a typical vertex on ∂Ω_n and by ^barycenter_ the barycenter of a typical finite element ∈^h. Then, for n=1,…,N_s, we define the subdomains
  
  
      {Γ_n^h    =
      {∈^h_Ω : ∃ ^vertex∈
                         such that  
      | ^vertex - _^barycenter| ≤δ/2+h}
      Ω_n^h    = {∈Ω_n∖ 
      (Ω_n∩Γ_n^h )   }
      Γ_n^h     = {∈^h_Γ : ∃ ^vertex∈
                         such that  
      | ^vertex - _^barycenter | ≤δ+h }.
      .
  
  
  We have that Γ_n^h consists of all elements ∈^h_Ω whose barycenters are within a distance δ/2+h of some element vertex on the Type 1 part of the boundary of Ω_n. Also, Γ_n^h consists of all elements ∈^h_Γ whose barycenters are within a distance δ+h of some element vertex on the Type 2 part of the boundary of Ω_n. Note that this procedure guarantees that the true interface region Γ_n is fully contained in the approximate interface region Γ_n^h. 
  
  
  We illustrate the above discussion in Figure <ref>. Note that all of that discussion applies even to the case of the single domain Ω being a rectangle and a single-domain FE grid that is Cartesian and uniform. Figure <ref>a depicts a portion of the  grid in the single domain Ω that respects the common boundary (depicted by the thick line segment) between two subdomains Ω_n and Ω_n'. For Figure <ref>b, we have that the orange subdomains depict portions of the subdomains Ω_n and Ω_n' and the blue domain depicts a portion of Γ_n∩Γ_n'. Note that the common boundaries of both Ω_n and Ω_n' with Γ_n∩Γ_n' do not respect the  grid so that Ω_n, Ω_n', and Γ_n∩Γ_n' all contain some partial (cut) triangles. Figure <ref>c illustrates the need to make changes to the single-domain FE grid so that the new grid does respect those common boundaries. Of course, if we define the subdomain FE discretization using the new grid of Figure <ref>c, there is no hope for the solution of the FE discretization of (<ref>) to be the same as the solution of single-domain FE discretization (<ref>), i.e., the goal (<ref>) cannot be achieved. Note also that the re-meshing of Figure <ref> is relatively easy to effect for Cartesian grids, but becomes a much more complex task for general grids, especially in three dimensions. Figure <ref>d illustrates the process defined in (<ref>). Now the orange-shaded regions depict portions of the subdomains Ω_n^h and Ω_n'^h and the magenta region depicts a portion of Γ_n^h∩Γ_n'^h. Note that, in Figure <ref>d, those three domains all contain only whole FE triangles. Also, the blue region in Figure <ref>b, i.e. Γ_n∩Γ_n', is fully contained in the set of magenta triangles.
  Figure <ref>e illustrates why, e.g., in the first equation in (<ref>), we used the criteria | ^vertex - _^barycenter| ≤δ/2+h and not simply the criteria | ^vertex - _^barycenter| ≤δ/2. The barycenters of the white triangles are such that | ^vertex - _^barycenter| > δ/2 so that using the latter criteria means that the white triangles are not included in Γ_n^h even though it is obvious from Figure <ref>b that those triangles overlap with Γ_n^h∩Γ_n'^h. On the other hand, Figure <ref>f illustrates that using the criteria | ^vertex - _^barycenter| ≤δ/2+h results in the yellow triangles which are in Γ_n^h∩Γ_n'^h but do not overlap with Γ_n∩Γ_n' so that those elements do not interact with elements on the other side of Γ_n^h∩Γ_n'^h. However, the equivalence between single- and multi-domain solutions is not compromised because in assembling the FE stiffness matrix, any additional entries in that matrix corresponding to points in the yellow elements are taken care of by weighting the integrand with the inverse of ζ_(,).
  
  
  
  
  Comparing the definitions of Γ_n^h, Γ_n^h, and Ω_n^h with the definitions of Γ_n, Γ_n, and Ω_n given in Section <ref>, one can certainly view the first trio as approximations to the second trio. However, this view does not intrude on any aspect of the developments that follow. For example, the accuracy of these domain approximations is, as is made evident in Section <ref>, irrelevant with respect to the goal stated in (<ref>).    
   
  
  
  
  
  
  
  
   §.§ Multi-domain finite element system
  
  
  Let W^h denote the finite element space spanned by the set of basis functions {ϕ_i()}_i=1^N^h used in Section (<ref>) to define the single-domain FE system (<ref>). Then,  let  
  
      W_n^h = {
      ϕ_i()   :  _i∈Ω_n^h∪Γ_n^h∪Γ_n^h}    ,
  
  
  i.e., W_n^h is spanned by the basis functions ϕ_i such that the associated node _i belongs to  Ω_n^h∪Γ_n^h∪Γ_n^h.
  The subspaces W_n^h ⊂ W^h are well defined because of (<ref>), i.e., because Ω_n^h, Γ_n^h, and Γ_n^h all consist of whole elements from the triangulation ^h. We also define the spaces
  
  
      W_n^0,h = {w^h()    ∈ W_n^h   :   
      w^h() = 0   }
             .
  
  Note that although we have that W_1^h⊗⋯⊗  W_N_s^h = W^h⊂ W, in general, W_n^h⊄W_n.
  
  We define
  
      ζ_^h(,)=
      ∑_n=1^N_sχ_Ω_n^h∪Γ_n^h∪Γ_n^h()
      χ_Ω_n^h∪Γ_n^h∪Γ_n^h()   
          ζ_^h()= ∑_n=1^N_sχ_Ω_n^h∪Γ_n^h().
  
  For n=1,…,N_s and for u_n^h(), v_n^h()∈ W_n^h, we define the discretized subdomain bilinear form 
  
  
      𝒜_n^h(u_n^h,v_n^h) =
        
         ∫_Ω_n^h∪Γ_n^h∪Γ_n^h∫_Ω_n^h∪Γ_n^h∪Γ_n^hζ^h_(,)^-1(v_n^h() - v_n^h())(u_n^h() - u_n^h())γ(,)d  d
  
  and the discretized subdomain linear functional
  
      _n^h(v_n^h) = ∫_Ω_n^h∪Γ_n^hζ^h_()^-1 v_n^h()f()d
  
  with the tacit understanding that for floating domains, i.e., if Γ_n^h=∅, the domains of integration of both integrals in (<ref>) reduce to Ω_n^h ∪Γ_n^h. 
  
  As done for the continuous problem, we introduce the function space
  
  
      W^h_n = {    W_n^0,h  ,     
         W_n^h   ,    .
  
  
  Then, we define the system of equations
  
      {   
         
              𝒜_n^h(u_n^h,v_n^h) = _n^h(v_n^h)     ∀  v_n^h∈W_n^h
      (a)
      
         
           
       u_n^h() = u_n'^h()     (b)
      
         
            u_n^h() = g^h()       ∀ ∈Γ_n^h
             ,   .
           (c)
      .
  
  
  
  The system (<ref>) is not a discretization of the system (<ref>) because  𝒜_n^h(·,·)𝒜_n(·,·), i.e., the former is defined with respect to Ω_n^h∪Γ_n^h∪Γ_n^h whereas the latter is defined in terms of Ω_n∪Γ_n∪Γ_n. However, this observation is unimportant because what is true is that
  
  
      u^h_n()=u^dd,h()|_Ω_n^h∪Γ_n^h∪Γ_n^h,
  
  
  which is, after all, the goal (<ref>) we want to achieve. The truthfulness of (<ref>) is verified following the same steps as those used to prove Proposition <ref> with, of course, Ω_n∪Γ_n∪Γ_n replaced by Ω_n^h∪Γ_n^h∪Γ_n^h. 
  
  When we define, in Section <ref>, the matrix form of (<ref>), it is useful to  differentiate between the bilinear forms in the two cases in (<ref>a). First, because Γ_n=∅ for floating domains, we have 
  
      𝒜_n^h(u_n^h,v_n^h)=𝒜_n^fl,h(u_n^h,v_n^h) = 
       
           ∫_Ω_n^h∪Γ_n^h∫_Ω_n^h∪Γ_n^hζ_^h(,)^-1(v_n^h() - v_n^h())(u_n^h() - u_n^h())γ(,)d  d.
  
  For non-floating domains, the bilinear form involves integrals with respect to Γ_n^h. For such integrals, we have that either v_n^h(·)=0 [because v_n^h∈ W_n^0,h] or u_n^h(·)=g^h(·) [because of the constraint (<ref>c)], so that then
  
      𝒜_n^h(u_n^h,v_n^h)=𝒜^nfl,h_n(u_n^h,v_n^h) - 𝒜^g,h_n(g^h,v_n^h),
   
  where
  
      𝒜^nfl,h_n(u_n^h,v_n^h) =
      
             ∫_Ω_n^h∪Γ_n^h∫_Ω_n^h∪Γ_n^hζ_^h(,)^-1(v_n^h() - v_n^h())(u_n^h() - u_n^h())γ(,)d  d
                 +
      2∫_Ω_n^h∪Γ_n^hu_n^h()v_n^h()∫_Γ_n^hζ^h_(,)^-1γ(,) d  d
  
  and
  
      𝒜^g,h_n(g^h,v_n^h)= 
           ∫_Ω_n^h∪Γ_n^h^∫_Γ_n^hζ_^h(,)^-1
      g^h()v_n^h()
      γ(,) d  d
             
      +∫_Γ_n^h∫_Ω_n^h∪Γ_n^hζ^h_(,)^-1g^h()v_n^h() 
       γ(,)  d  d.
  
  Note that in (<ref>), (<ref>), and (<ref>), both u_n^h(·) and v_n^h(·) are evaluated only at points in Ω_n^h∪Γ_n^h and g_n^h(·) is evaluated only at points in Γ_n^h.
  
  In light of (<ref>), (<ref>), and (<ref>),
  (<ref>) can be rewritten as, for floating domains,
  
      {𝒜_n^fl,h(u_n^h,v_n^h)    = _n^h(v_n^h)
      (a)
      
      
       u_n^h()    = u_n'^h()              (b)
      .
  
   and for non-floating domains
  
      {𝒜_n^nfl,h(u_n^h,v_n^h)    = _n^h(v_n^h) + 𝒜^g,h_n(g^h,v_n^h)
      (a)
      
      
       u_n^h()    = u_n'^h()           (b)
      .
  
  
  
  For any domain indexed by n, each of the subdomain problems in (<ref>) and (<ref>) is coupled, through (<ref>)b or (<ref>)b, to other domains indexed by n' with n' n. Of course, this defeats the goal of domain decomposition which is to construct uncoupled subdomain problems so that, e.g., parallelization can be realized. This becomes the task for algorithms of obtaining solutions of the subdomain problems. Further comments in this regard are provided in Section <ref>.
      
  
  
  
  If instead of (<ref>a) we consider
  
      -2∫_Ω∪Γ∪Γ_Neumann(u() -u())γ(,)d 
      +c()u()
      = f_Ω()
  
  with c()>0 so that the bilinear form 𝒜(u,v) in (<ref>) has the additional term ∫_Ω c()u()v()d, then that bilinear form is coercive even for floating subdomains. In this case, the design of solution methods for (<ref>), and in particular for (<ref>), becomes substantially simpler.
      
  
  
  
  The discussion that includes (<ref>)–(<ref>) as well as the comments made in Remarks <ref> and <ref> about the discrete bilinear form 𝒜_n^h(·,·) and the discrete subdomain system (<ref>) also hold for the continuous bilinear form 𝒜_n(·,·) defined in (<ref>) and continuous subdomain system (<ref>). 
      
  
  
  
  
  
    §.§.§ Matrix form of the multi-domain finite element system
  
  
  Based on the numbering introduced in Section <ref>, we let
  
  X_Ω∪Γ={_i}_i=1^N^h denote the set of nodes in the grid used for the FE discretization of the single-domain Ω∪Γ
  
  X_Ω={_i}_i=1^N^h 
  and
  X_Γ={_i}_i=N^h+1^N^h denote the set of nodes in Ω and  Γ, respectively. 
  
  The sets X_Ω⊂ X_Ω∪Γ and X_Γ⊂ X_Ω∪Γ are disjoint and X_Ω∪ X_Γ=X_Ω∪Γ. For n=1,…,N_s, let 
  
  X_n^h={_i^n}_i=1^N_n^h denote a local numbering of the set of nodes in Ω_n^h∪Γ_n^h∪Γ_n^h; clearly, by construction, X_n^h⊂ X_Ω∪Γ
  
  X_n^h={_i^n}_i=1^N_n^h and  X_Γ_n^h^h={_i^n}_N_n^h+1^N_n^h denote the nodes in X_n^h located in Ω_n^h∪Γ_n^h and Γ_n^h, respectively.
  
  The sets X_n^h⊂X_n^h and X_Γ_n^h^h⊂X_n^h are disjoint and X_n^h∪ X_Γ_n^h^h=X_n^h. Note that for n' n, the sets X_n^h and X_n'^h overlap whenever Γ_n^h∩Γ_n'^h∅ and similarly for the sets X_Γ_n^h and X_Γ_n'^h. Also note that if Ω_n^h is a floating domain, then the set X_Γ_n^h is vacuous.
  
  Bilinear forms in matrix notation. We first consider the conversion of (<ref>a) and (<ref>a) to matrix notation. Corresponding to the nodes in X_Ω∪Γ, we have the set of basis functions {ϕ_i()}_i=1^N^h whose span is used to define the finite element space W^h for the single-domain finite element system. We introduce the set of basis functions corresponding to each subdomain. For n=1,…,N_s, let
  
  {ϕ_i^n()}_i=1^N_n^h denote a local numbering of the basis functions in the spanning set for W^h which correspond to the nodes in X_n^h. 
  
  By construction, {ϕ_i^n()}_i=1^N_n^h spans the finite element space W_n^h. We then 
  
  divide {ϕ_i^n()}_i=1^N_n^h into the sets {ϕ_i^n()}_i=1^N_n^h and {ϕ_i^n()}_i=N_n^h+1^N_n^h that correspond to nodes in X_n^h and X_Γ_n^h^h, respectively.
  
  Note that if Ω_n is a floating domain, then N_n^h=N_n^h so that the set {ϕ_i^n()}_i=N_n^h+1^N_n^h is vacuous.
  
  
  
  Let  denote an N^h_n-vector of nodal values of a function u^h_n() defined for the nodes in X_n^h (i.e., nodes in Ω^h_n∪Γ^h_n) and let  denote the (N^h_n-N^h_n)-vector of nodal values of g^h() defined for the nodes in X_Γ_n^h^h (i.e., nodes in Γ^h_n). Then, we have that
  
  
      u^h_n() = 
          {∑_j=1^N^h_n_j ϕ_j^n()
          + ∑_j=N^h_n+1^N^h_n_j ϕ_j^n()    for non-floating domains
      ∑_j=1^N^h_n_j ϕ_j^n()    for floating domains..
  
  Note that components of  are ordered according the local indexing of nodes.
  
  We then define, for n=1,…,N_s, the entries of the N_n^h× N_n^h matrix 𝔸_n as
  
  
      (𝔸_n)_ij =
         {𝒜^nfl,h_n(ϕ_j^n,ϕ_i^n)    for non-floating domains
      𝒜^fl,h_n(ϕ_j^n,ϕ_i^n)    for floating domains.
  
  
  and, for i=1,…,N_n^h, the components of the N_n^h-vector  as
  
  
      _i =
         {ℱ^h_n(ϕ_i^n)-∑_j=N^h+1^N^h𝒜^g,h_n(ϕ_j^n,ϕ_i^n)_j    for non-floating domains
      ℱ^h_n(ϕ_i^n)    for floating domains.
          .
  
  Then, the finite element problems (<ref>a) and (<ref>a) both have the matrix notation equivalent
  
      𝔸_n =     .
  
  
  Constraints in matrix notation. We next turn to the conversion of the constraints in (<ref>b) [or equivalently (<ref>b) and (<ref>b)] to matrix notation. 
  
  We first transform those constraints to vector notation. We have global and local indices of nodes. Thus, if the node _j^n∈ X_n corresponds to the node _i∈ X, i.e., if we have _i = _j^n, then i and j are the global and local indices, respectively, for the same node. We define a mapping from global to local indices, specifically, for a globally indexed node _i∈Γ_n^h, we let
  
      I_ni=
  
  Now suppose that Γ_n^h∩Γ_n'^h∅. Then, for a globally indexed node _i∈Γ_n^h∩Γ_n'^h, we have that
  
      _i = _j^n = _j'^n'    .
  
  Because the FE approximation u_n^h() is uniquely determined by its nodal values, the constraints in (<ref>b) can be equivalently expressed as
  
      _j = (_n')_j'  ,  
          n'= n+1,…,N_s.
  
  We keep in mind that, as was the case for the continuous multi-domain system (see Remark <ref>), the constraints in (<ref>) are not independent. For example, if Γ_1^h ∩Γ_2^h∩Γ_3^h∅, then, for a node _i in that domain, we have from (<ref>) that
  (_1)_j = (_2)_j', (_1)_j = (_3)_j”, and (_2)_j' = (_3)_j”, where j=I_1i, j'=I_2i, and j”=I_3i.
  Clearly, these three equations are not independent.
  
  Our task is then reduced to expressing the constraints in (<ref>) in an economical matrix form. Here, we mimic the process given in <cit.> for DD in the local PDE case, In fact, we construct an M×∑_n=1^N_sN_n^h matrix 𝕄, where M is specified below, of the form
  
      𝕄 = ( 𝕄_1 ⋯𝕄_N_s )
  
  such that the constraints in (<ref>) can be equivalently expressed as
  
      ∑_n=1^N_s𝕄_n  = 0,
  
  where 𝕄_n, n=1,…,N_s, are M× N_n^h matrices. We construct two such matrices, one for the constraints in (<ref>), the other for an equivalent non-redundant set of constraints. To this end, letting Γ^h=∪_n=1^N_sΓ_n^h,
  
      θ(_i)= {n  :  _i∈Γ_n^h}
  
  so that θ(_i) consists of the indices of all the subdomains Γ_n^h that contain the globally indexed node _i and m(_i) denotes the number of distinct subdomains which the node _i belongs to.
  
  Non-redundant constraints – full-rank matrix 𝕄.
  For each node _i∈Γ^h,
  
      - arrange the indices in θ(_i) in increasing order
  
      - for each consecutive pair of indices, impose one constraint.
  
  This results in, with N_Γ^h denoting the number of nodes in Γ^h,
  
      M = ∑_i=1^N_Γ^h( m(_i) - 1 )
  
  non-redundant constraints. For example, if Γ_1^h ∩Γ_2^h∩Γ_3^h∅, then, for a node _i in that domain, we now have that
  (_1)_j = (_2)_j' and (_2)_j' = (_3)_j”, where j=I_1i, j'=I_2i, and j”=I_3i. Clearly, these two equations are not redundant and together imply the constraint (_1)_j  = (_3)_j” from (<ref>) that is now missing.
  
  The entries of the matrix 𝕄 can be determined as follows: set k=0 and then,
  
      - for i=1,…, N_Γ^h
  
      - for each pair (n,n'), n<n', of consecutive indices in θ(_i)
  
      - set k⟵ k+1
  
      - for j and j' such that ^n_j=^n'_j'=_i, 
  
      (𝕄_n)_kj = 1      
      (𝕄_n')_kj' = -1
  
  and all other entries in the k-th row of 𝕄 to zero.
  
  Redundant constraints
  We proceed as we do above for the non-redundant set of constraints, except that now we do not require that n<n' be consecutive indices in θ(_i), i.e., we impose a constraint for every distinct pair n<n' of indices in θ(_i). This approach lends itself better for parallelization compared to the use of non-redundant constraints; see, e.g., <cit.> for a discussion in the local DD setting. Thus, we impose all the constraints in (<ref>). Note that redundant constraints are caused only for m(_i)≥ 3. Because there are m(_i) distinct indices in θ(_i), we have that the number of rows in the matrix 𝕄 is now given by
  
      M = ∑_i=1^N_Γ^h1/2 m(_i)( m(_i) - 1 ).
  
  
  
  
  
  
  
  § CONCLUDING REMARKS
  
  
  We have defined and analyzed a general framework for the construction of domain decomposition methods for nonlocal problems that achieves the goal stated in (<ref>) or, more precisely, in (<ref>). However, there is still work to be done because we have not met a second goal which is that the nonlocal DD method is amenable to parallelization. What we have so far are N_s systems (<ref>)–(<ref>) in which the subdomain systems are coupled through the constraints (<ref>b) and (<ref>b) or, equivalently, the subdomain matrix systems in (<ref>) that are coupled through the constraints in (<ref>). Such couplings prevent the direct use of (<ref>)–(<ref>) (or equivalently (<ref>)–(<ref>)) for achieving the second goal.
  
  At a similar stage in the development DD algorithms in the local PDE setting, one is faced with the analogous situation; for example, in the non-overlapping DD setting, there is coupling of the subdomain problems, say in the matrix formulation, at the nodes located along the common boundaries between subdomains. The uncoupling between subdomain problems is then effected through the design of solution methods in which the bulk of the computational effort is borne by steps in those methods that are parallelizable and for which the non-paralellizable steps and the communications between subdomain problems, i.e., between processors, is kept to a minimum. 
  
  In the nonlocal DD setting, solution methods have to be designed to meet the same criterion: the bulk of the computational costs has to be borne by parallelizable steps. In a follow-up paper, we will develop, analyze, and implement such methods. Parallelizable solution methods for local PDE non-overlapping DD will be generalized to the nonlocal setting. For example, Lagrange multiplier methods, e.g., FETI <cit.>, Arlequin methods <cit.>, and optimization-based DD methods <cit.>, all of which are in use for local PDE non-overlapping DD and which are also all good candidates for generalization to the nonlocal DD setting. Other local DD solution methods could also be considered for generalization. As is the case for the framework developed in this paper, generalizations of solution methods will pose challenges because of nonlocality. One thing to keep in mind is that a preferred solution method in the local setting may or may not remain so when generalized to the nonlocal setting.
  
  
  
  
  
  
  
  § ACKNOWLEDGMENTS
  
  Sandia National Laboratories is a multimission laboratory managed and operated by National Technology and Engineering Solutions of Sandia, LLC., a wholly owned subsidiary of Honeywell International, Inc., for the U.S. Department of Energy's National Nuclear Security Administration under contract DE-NA-0003525. This paper describes objective technical results and analysis. Any subjective views or opinions that might be expressed in this paper (SAND2020-8735) do not necessarily represent the views of the U.S. Department of Energy or the United States Government.
  
  This work was supported by the Sandia National Laboratories (SNL) Laboratory-directed Research and Development (LDRD) program, and the U.S. Department of Energy, Office of Science, Office of Advanced Scientific Computing Research under Award Number DE-SC-0000230927 and under the Collaboratory on Mathematics and Physics-Informed Learning Machines for Multiscale and Multiphysics Problems (PhILMs) project.
  
  
  
  
  
  
  
  
  
  
  
  
  § PROOF OF LEMMA <REF>
  
  It is convenient for what follows to introduce, for each n=1,…,N_s such that Γ_n∅, the splitting
  
      Γ_n = Γ_n^*∪Γ_n^†    {   Γ_n^*⊂Γ_n
            Γ_n^*∩(∪_n'=1, n' n^N_sΓ_n')=∅
         Γ_n^†=Γ_n∖Γ_n^*⊂Γ_n
      .
  
  so that Γ_n^* (resp. Γ_n^†) are the disjoint parts of Γ_n that do not (resp. do) overlap with any other Γ_n' with n' n. The blue regions in Figure <ref>-right illustrate examples of the sets Γ_n^†. Note that, by definition,
  Γ_n^*∩Γ_n^†=∅. 
  With these definitions in hand, we have that, for n=1,…,N_s, 
  Ω_n∪Γ_n∪Γ_n =Ω_n∪Γ_n∪Γ_n^*∪Γ_n^†.
  We can then express (<ref>) as
  
      ζ_𝒜(,)
      =
      ∑_n=1^N_s_Ω_n∪Γ_n^*∪Γ_n∪Γ_n^†()
      _Ω_n∪Γ_n^*∪Γ_n∪Γ_n^†().
  
  
  Note that among the sets Ω_n, Γ_n^*, Γ_n, and Γ_n^† the only two that may possibly overlap with other sets are Γ_n^† and Γ_n. For this reason and for ease of notation, we further introduce the set Γ_n=Γ_n^†∪Γ_n. Based on this consideration, we split the outer integral of the bilinear form _n into a set that does not overlap with any other sets (i.e. Ω_n∪Γ_n^*) and Γ_n. We have:
  
  
      𝒜_n(u_n,v_n) 
              = ∫_Ω_n∪Γ_n^*∪Γ_n∫_Ω_n∪Γ_n^*∪Γ_nζ_𝒜(𝐱,𝐲)^-1(v_n(𝐲) - v_n(𝐱))(u_n(𝐲) - u_n(𝐱))γ(𝐱,𝐲)d𝐲  d𝐱
              = 𝒜_n^disjoint(u_n,v_n) + 𝒜_n^overlap(u_n,v_n),
  
  where
  
      𝒜_n^disjoint(u_n,v_n) := 
             ∫_Ω_n∪Γ_n^*∫_Ω_n∪Γ_n^*∪Γ_nζ_𝒜(𝐱,𝐲)^-1(v_n(𝐲) - v_n(𝐱))(u_n(𝐲) - u_n(𝐱))γ(𝐱,𝐲)d𝐲  d𝐱,
  
  
  
      𝒜_n^overlap(u_n,v_n) := 
             ∫_Γ_n∫_Ω_n∪Γ_n^*∪Γ_nζ_𝒜(𝐱,𝐲)^-1(v_n(𝐲) - v_n(𝐱))(u_n(𝐲) - u_n(𝐱))γ(𝐱,𝐲)d𝐲  d𝐱.
  
  
  To simplify the notation, we let w_n(,)= (v_n(𝐲) - v_n(𝐱))(u_n(𝐲) - u_n(𝐱))γ(𝐱,𝐲) and w(,) = (v(𝐲) - v(𝐱))(u(𝐲) - u(𝐱))γ(𝐱,𝐲). Note that whenever  or  belong to a set that overlaps with other sets, the definitions of u_n and v_n in (<ref>) guarantee that w_n(,)=w(,). We first analyze _n^disjoint; we have that
  
  
      𝒜_n^disjoint(u_n,v_n) 
         =∫_Ω_n∪Γ_n^*∫_Ω_n∪Γ_n^*∪Γ_n w(𝐱,𝐲)d d
         = ∫_Ω_n∪Γ_n^*∫_Ω_n∪Γ_n^*∪Γ_n w(𝐱,𝐲)d d
         +  ∫_Ω_n∪Γ_n^*∫_(Ω∪Γ)∖(Ω_n∪Γ_n^*∪Γ_n)w(𝐱,𝐲)d𝐲  d𝐱
  
  
  where the first equality follows from the fact that ζ_𝒜(𝐱,𝐲)=1 for 𝐱∈Ω_n ∪Γ_n^* and ∈Ω_n∪Γ_n^*∪Γ_n and the second inequality from the fact that γ(𝐱,𝐲)=0 for 𝐱∈Ω_n ∪Γ_n^* and ∈(Ω∪Γ)∖(Ω_n∪Γ_n^*∪Γ_n). Hence,
  
  
      ∑_n=1^N_s𝒜_n^disjoint(u_n,v_n) = ∫_⋃_n=1^N_s(Ω_n∪Γ_n^*)∫_Ω∪Γ w(𝐱,𝐲)d𝐲  d𝐱.
  
  
  By definition of ζ_𝒜, for any 𝐲∈ (Ω_n∪Γ_n^*) and for any 𝐱∈Γ_n we have ζ_𝒜(𝐱,𝐲)=1, hence 
  
  
      𝒜_n^overlap(u_n,v_n)    = ∫_Γ_n∫_Ω_n∪Γ_n^* w_n(𝐱,𝐲)  d𝐲  d𝐱
         + ∫_Γ_n∫_Γ_nζ_𝒜(𝐱,𝐲)^-1
        w_n(𝐱,𝐲) d𝐲  d𝐱
         = 𝒜_n^overlap,I(u_n,v_n) + 𝒜_n^overlap,II(u_n,v_n).
  
  
  We introduce the partition {Λ_i}_i=1^N_λ of ⋃_n=1^N_sΓ_n, i.e. 
  
      ⋃_i=1^N_λΛ_i = ⋃_i=1^N_sΓ_n,    and  Λ_i∩Λ_j=∅,   for  i≠ j,
  
  such that there exists an index set I_n, n=1, … N_s, for which 
  
      ⋃_i∈ I_nΛ_i = Γ_n.
  
  In practice, this means that for every n=1,… N_s there exists a subset of the partition that provides a partition of Γ_n. Note that such a partition always exists. Then we have
  
  
      ∑_n=1^N_s𝒜_n^overlap,I(u_n,v_n) 
         = ∑_n=1^N_s∫_Γ_n∫_Ω_n∪Γ_n^*w_n(𝐱,𝐲) d𝐲  d𝐱
         = ∑_n=1^N_s∑_i∈ I_n∫_Λ_i∫_Ω_n∪Γ_n^*w_n(𝐱,𝐲) d𝐲  d𝐱
         = ∑_n=1^N_s∑_i=1^N_λ∫_Λ_i∫_Ω_n∪Γ_n^*w_n(𝐱,𝐲) d𝐲d
         = ∑_i=1^N_λ∫_Λ_i[∑_n=1^N_s∫_Ω_n∪Γ_n^*w_n(𝐱,𝐲) d]d
         =∫_⋃_n=1^N_sΓ_n∫_⋃_n=1^N_s(Ω_n∪Γ_n^*)w(𝐱,𝐲) d𝐲  d𝐱,
  
  
  where the second equality follows from the fact that {Λ_i}_i∈ I_n form a disjoint covering of Γ_n, the third from the fact that γ(,)=0 for all the extra terms, the fourth from the fact that the sums are independent and hence can be switched, and the fifth from the fact that the sets are disjoint.
  
  For 𝒜_n^overlap,II we have 
  
  
      ∑_n=1^N_s𝒜_n^overlap,II(u_n,v_n) 
         = ∑_n=1^N_s∫_Γ_n∫_Γ_nζ_A(,)^-1w(,) d𝐲  d𝐱
         = ∑_n=1^N_s∑_i∈ I_n∑_j∈ I_n∫_Λ_i∫_Λ_jζ_A(,)^-1 w(,) d𝐲  d𝐱,
  
  
  where, again, the second equality follows from the fact that the sets {Λ_i}_i∈ I_n are disjoint and form a covering of Γ_n.
  
  Next, we introduce the index set _ij(,) that contains all indices n such that (,)∈Λ_i×Λ_j and i,j∈ I_n. Formally,
  
  
      _ij(,)={n∈{1,…,N_s}  :  ∈Λ_i,  ∈Λ_j,   and  i, j∈ I_n}.
  
  
  
  
  
  
  
  
  Note that for ∈Λ_i and ∈Λ_j 
  
      ζ_(,) 
         = ∑_n=1^N_s_Ω_n∪Γ_n^*∪Γ_n()
         _Ω_n∪Γ_n^*∪Γ_n() 
         = ∑_n=1^N_s__ij(,)(n)
         _Λ_i()_Λ_j()
      =  ∑_n=1^N_s__ij(,)(n).
  
  
  Thus,
  
  
      ∑_n=1^N_s   𝒜_n^overlap,II(u_n,v_n) 
      = ∑_n=1^N_s∑_i∈ I_n∑_j∈ I_n∫_Λ_i∫_Λ_jζ_A(,)^-1 w(,) d d
         = ∑_n=1^N_s∑_i∈ I_n∑_j∈ I_n∫_Λ_i∫_Λ_j__ij(,)(n) _Λ_i()_Λ_j()ζ_A(,)^-1 w(,) d d
         = ∑_n=1^N_s∑_i=1^N_λ∑_j=1^N_λ∫_Λ_i∫_Λ_j__ij(,)(n) _Λ_i()_Λ_j()ζ_A(,)^-1 w(,) d d
         = ∑_i=1^N_λ∑_j=1^N_λ∫_Λ_i∫_Λ_jζ_A(,)ζ_A(,)^-1w(,) d d
          = ∫_⋃_n=1^N^sΓ_n∫_⋃_n=1^N^sΓ_n w(,)d d.
  
  
  Here, in the second equality we only added terms that are equal to 1 because all the indicator functions are active. This allows us to extend the sums over i and j in the third equality because all the extra terms are zero. Then, the fourth equality follows from (<ref>) and the fifth from the fact that all sets Λ_i and Λ_j are disjoint.
  
  It then follows that 
  
      ∑_n=1^N_s𝒜_n^overlap(u_n,v_n)  =  ∑_n=1^N_s(𝒜_n^overlap,I(u_n,v_n)  
      +𝒜_n^overlap,II(u_n,v_n))
         =∫_⋃_n=1^N_sΓ_n( ∫_⋃_n=1^N_s(Ω_n∪Γ_n^*)w(𝐱,𝐲) d𝐲  +∫_⋃_n=1^N_sΓ_nw(𝐱,𝐲) d𝐲)d𝐱
         = ∫_⋃_n=1^N_sΓ_n∫_Ω∪Γ w(𝐱,𝐲) d𝐲 d𝐱.
  
    The last equality follows from the fact that ⋃_n=1^N_s(Ω_n∪Γ_n^*) and ⋃_n=1^N_s (Γ_n^†∪Γ_n) are two disjoint sets. For the same reason, we have
  
      ∑_n=1^N_s𝒜_n(u_n,v_n) =  ∑_n=1^N_s(𝒜_n^disjoint(u_n,v_n)  
      +𝒜_n^overlap(u_n,v_n))
         ∫_⋃_n=1^N_s(Ω_n∪Γ_n^*)∫_Ω∪Γ w(𝐱,𝐲)d𝐲  d𝐱 + ∫_⋃_n=1^N_sΓ_n∫_Ω∪Γ w(𝐱,𝐲) d𝐲 d𝐱
         = ∫_Ω∪Γ∫_Ω∪Γ w(𝐱,𝐲) d𝐲 d𝐱 = 𝒜(u,v).
  
  
  The proof for  directly follows from the same arguments and it is not reported.
  
  
  
  
  
  
  
  