Formulation, analysis and computation of an optimization-based local-to-nonlocal coupling method

We present an optimization-based coupling method for local and nonlocal continuum models. Our approach couches the coupling of the models into a control problem where the states are the solutions of the nonlocal and local equations, the objective is to minimize their mismatch on the overlap of the local and nonlocal problem domains, and the virtual controls are the nonlocal volume constraint and the local boundary condition. We present the method in the context of Local-to-Nonlocal diffusion coupling. Numerical examples illustrate the theoretical properties of the approach.


Introduction
Nonlocal continuum theories such as peridynamics [30], physics-based nonlocal elasticity [12], or nonlocal descriptions resulting from homogenization of nonlinear damage models [23] can incorporate strong nonlocal effects due to long-range forces at the mesoscale or microscale. As a result, for problems where these effects cannot be neglected, such descriptions are more accurate than local Partial Differential Equations (PDEs) models. However, their computational cost is also significantly higher than that of PDEs. Local-to-Nonlocal (LtN) coupling methods aim to combine the computational efficiency of PDEs with the accuracy of nonlocal models. The need for LtN couplings is especially acute when the size of the computational domain is such that the nonlocal solution becomes prohibitively expensive to compute, yet the nonlocal model is required to accurately resolve small scale features such as crack tips or dislocations that can affect the global material behavior.
LtN couplings involve two fundamentally different mathematical descriptions of the same physical phenomena. The principal challenge is the stable and accurate merging of these descriptions into a physically consistent coupled formulation. In this paper we address this challenge by couching the LtN coupling into an optimization problem. The objective is to minimize the mismatch of the local and nonlocal solutions on the overlap of their respective subdomains, the constraints are the associated governing equations, and the controls are the virtual nonlocal volume constraint and the local boundary condition. We formulate and analyze this optimization-based LtN approach in the context of local and nonlocal diffusion models [16].
Our coupling strategy differs fundamentally from other LtN approaches such as the extension of the Arlequin [11] method to LtN couplings [23], force-based couplings [29], or the morphing approach [5,25]. The first two schemes blend the energies or the forces of the two models over a dedicated "gluing" area, while the third one implements the coupling through a gradual change in the material properties characterizing the two models over a "morphing" region. In either case, resulting LtN methods treat the coupling condition as a constraint, similar to classical domain decomposition methods. In contrast, we treat this condition as an optimization objective, and keep the two models separate. This strategy brings about valuable theoretical and computational advantages. For instance, the coupled problem passes a patch test by construction and its well-posedness typically follows from the well-posedness of the constraint equations.
Our approach has its roots in non-standard optimization-based domain decomposition methods for PDEs [13,14,15,19,20,21,22,24]. It has also been applied to the coupling of discrete atomistic and continuum models in [27,28] and multiscale problems [1]. This paper continues the efforts in [8], which presented an initial optimization-based LtN formulation and in [10], which focussed on specializing the formulation to mixed boundary conditions and mixed volume constraints and its practical demonstration using Sandia's agile software components toolkit. The main contributions of this paper include (i) rigorous analysis of the LtN coupling error, (ii) formal proof of the well-posedness of the discretized LtN formulation, and (iii) rigorous convergence analysis.
We have organized the paper as follows. Section 2 introduces notation, basic notions of nonlocal vector calculus and the relevant mathematical models. We present the optimization-based LtN method and prove its well-posedness in Section 3 and study its error in Section 4. Section 5 focusses on the discrete LtN formulation, its well-posedness and numerical analysis. A collection of numerical examples in Section 6 illustrates the theoretical properties of the method using a simple one-dimensional setting.

Preliminaries
Let ω be a bounded open domain in R d , d = 2, 3, with Lipschitz-continuous boundary ∂ω. We use the standard notation H s (ω) for a Sobolev space of order s with norm and inner product · s,ω and (·, ·) s,ω , respectively. As usual, H 0 (ω) := L 2 (ω) is the space of all square integrable functions on ω. The subset of all functions in H 1 (ω) that vanish on ζ ⊂ ∂ω is H 1 ζ (ω) := {u ∈ H 1 (ω) : u| ζ = 0}. The nonlocal model in this paper requires nonlocal vector calculus operators [16, §3.2] acting on functions u(x) : R d → R and ν(x, y) : non-negative symmetric kernel and an antisymmetric function, respectively, i.e., γ(x, y) = γ(y, x) ≥ 0 and α(y, x) = −α(x, y). The nonlocal diffusion 1 of u is an operator L(u) : R d → R defined by and its nonlocal gradient is a mapping G(u) : Finally, the nonlocal divergence of ν(x, y) is a mapping D(ν) : Furthermore, given a second-order symmetric tensor Φ(x, y) = Φ(y, x), equations (1) imply that Thus, with the identification γ(x, y) := α(x, y) · Φ(x, y)α(x, y) the operator L is a composition of the nonlocal divergence and gradient operators: Lu = D ΦGu). We define the interaction domain of an open bounded region ω ∈ R d as η = {y ∈ R d \ ω : γ(x, y) = 0}, for x ∈ ω and set Ω = ω ∪ η. In this paper we consider kernels γ such that for where B ε (x) = {y ∈ R d : |x − y| ≤ ε}. Kernels that satisfy (2) are referred to as localized kernels with interaction radius ε. It is easy to see that for such kernels the interaction domain is a layer of thickness ε that surrounds ω, i.e. η = {y ∈ R d \ ω : inf x∈ω |y − x| ≤ ε}; see Fig. 1 for a two-dimensional example. For a symmetric positive definite tensor Φ we respectively define the nonlocal energy semi-norm, nonlocal energy space, and nonlocal volume-constrained energy space by We also define the volume-trace space V (η) := {v| η : v ∈ V (Ω)}, and an associated norm We refer to [16,17] for further information about the nonlocal vector calculus. In order to avoid technicalities not germane to the coupling scheme, in this paper we consider integrable kernels. Examples of applications modeled by the latter can be found in [2,3,4]. Specifically, we assume that there exists positive constants γ 0 and γ 2 such that for all x ∈ Ω. Note that this also implies that there exists a positive constant γ 1 such that for all x ∈ Ω Ω γ(x, y) dy ≤ γ 1 .
In [16, §4.2] this class of kernels (referred to as Case 2) is rigorously analyzed; we report below an important result, useful throughout the paper.
Lemma 2.1: Let the function γ satisfy (2) and (5), then, there exist positive constants C pn and C * such that Furthermore, the energy space The latter is a combination of results in Lemmas 4.6 and 4.7 and Corollary 4.8 in [16, §4.3.2]. Note that the lower bound in (7) represents a nonlocal Poincaré inequality. Even though not included in the analysis, singular kernels appear in applications such as peridynamics; numerical results, included in the paper, suggest that the coupling scheme can handle such kernels without difficulties. However, their analysis is beyond the scope of this paper.

Local-to-Nonlocal coupling setting
Consider a bounded open region ω ⊂ R d with interaction domain η. Given f n ∈ L 2 (ω) and σ n ∈ V (η) we assume that the volume-constrained 3 nonlocal diffusion equation provides an accurate description of the relevant physical processes in Ω = ω ∪ η. Let Γ = ∂Ω, we assume that the local diffusion model given by the Poisson equation with suitable boundary data σ l ∈ H 1 2 (Γ) and forcing term f l ∈ L 2 (Ω) is a good approximation of (8) whenever the latter has sufficiently "nice" solutions. In this work we define f l to be an extension of f n by 0 in η, specifically, For a symmetric positive definite Φ standard arguments of variational theory show that the weak form 4 Ω Ω of (8) is well-posed [16], i.e., (11) has a unique solution such that for some positive constant K n . In this work, for simplicity and without loss of generality, we set Φ = I. Although (11) and the nonlocal calculus [16] enable formulation and analysis of finite elements for (8), which parallel those for the Poisson equation (9), resulting methods may be computationally intractable for large domains. The root cause for this is that long-range interactions increase the density of the resulting algebraic system making it more expensive to assemble and solve.

Optimization-based LtN formulation
For clarity we consider (8) and (9) with homogeneous Dirichlet conditions on η and Γ, respectively. To describe the approach it suffices to examine a coupling scenario where these problems operate on two overlapping subdomains of Ω. Thus we consider partitioning of Ω into a nonlocal subdomain ω n with interaction volume η n and a local subdomain Ω l , with boundary Γ l , such that Ω n := ω n ∪ η n ⊂ Ω, Ω = Ω n ∪ Ω l and Ω o = Ω n ∩ Ω l = ∅; see Fig. 2. Let η D = η ∩ η n , η c = η n \ η D , Γ D = Γ ∩ Γ l , and Γ c = Γ l \ Γ D ; see Fig. 2 and Appendix B for a summary of notation and definitions. Restrictions of (8) and (9) to ω n and Ω l are given by respectively, where θ n ∈ Θ n = {v n | ηc : v n ∈ V η D (Ω n )} and θ l ∈ Θ l = {v l | Γc : v l ∈ H 1 Γ D (Ω l )} are an undetermined Dirichlet volume constraint and an undetermined Dirichlet boundary condition, respectively. The following constrained optimization problem min un,u l ,θn,θ l J(u n , u l ) subject to (13), where J(u n , u l ) = 1 2 defines the optimization-based LtN coupling. In this formulation the subdomain problems (13) are the optimization constraints, u n and u l are the states and θ n and θ l are the controls. We equip the control space Θ n × Θ l with the norm (σ n , σ l ) 2 Θn×Θ l = σ n 2 0,ηc + σ l In contrast to blending, (14) is an example of a divide-and-conquer strategy as the local and nonlocal problems operate independently in Ω n and Ω l .
Given an optimal solution (u * n , u * l , θ * n , θ * (14) we define the LtN solution u * ∈ L 2 (Ω) by splicing together the optimal states: The next section verifies that (14) is well-posed.

Well-posedeness
We show that for any pair of controls subproblems (13) have unique solutions u n (θ n ) and u l (θ l ), respectively. Elimination of the states from (14) yields the equivalent reduced space form of this problem in terms of θ n and θ l only: To show that (17) is well-posed we start as in [19,27] and split, for any given (θ n , θ l ), the solutions of the state equations into a harmonic and a homogeneous part. The harmonic components v n (θ n ) and v l (θ l ) of the states solve the equations respectively. The homogeneous components u 0 n and u 0 l solve a similar set of equations but with homogeneous volume constraint and boundary condition, respectively: x ∈ η n and −∆u 0 In terms of these components u n = v n (θ n ) + u 0 n , u l = v l (θ l ) + u 0 l , and the objective The Euler-Lagrange equation of (17) is given by: seek (σ n , σ l ) ∈ Θ n × Θ l such that where The following lemma establishes a key property of Q.
As a results, Q(·; ·) endows the control space Θ n × Θ l with the "energy" norm Note that Q and F are continuous with respect to the energy norm. However, the control space Θ n × Θ l may not be complete with respect to the energy norm. In this case, following [1,13], we consider the optimization problem (17) on the completion Θ n × Θ l of the control space.
Theorem 3.1: Let Θ n × Θ l denote a completion 5 of the control space with respect to the energy norm (21). Then, the minimization problem min has a unique minimizer ( θ * n , θ * l ) ∈ Θ n × Θ l such that Proof. Equation (23) is a necessary condition for any minimizer of (22). Assume first that the control space is complete, i.e. Θ n × Θ l = Θ n × Θ l . Then Θ n × Θ l is Hilbert and the projection theorem implies that (23) has a unique solution. When Θ n × Θ l is not complete, the continuous bilinear form Q and the continuous functional F defined on Θ n × Θ l can be uniquely extended by using the Hahn-Banach theorem to a continuous bilinear form Q and a continuous functional F in Θ n × Θ l . Then, the existence and uniqueness of the minimizer follow as before.
To avoid technical distractions, in what follows we assume that the minimizer (θ * n , θ * l ) belongs to Θ n × Θ l and hence u * n = u n (θ * n ) ∈ V η D (Ω n ) and u * l = u l (θ * l ) ∈ H 1 Γ D (Ω l ). We note that in the finite dimensional case the completeness is not an issue, as the discrete control space is Hilbert with respect to the discrete energy norm, see Section 5.

Analysis of the LtN coupling error
We define the LtN coupling error as the L 2 -norm of the difference between the global nonlocal solution u n of (8) with homogeneous volume constraints and the LtN solution u * ∈ L 2 (Ω) given by (16). This section shows that the coupling error is bounded by the modeling error on the local subdomain, i.e., the error made by replacing the "true" nonlocal diffusion operator on Ω l by the Laplacian.
We prove this result under the following assumptions.
We also need the trace operator T : and the lifting operator L : are a harmonic lifting operator and the homogenous part of the states, respectively. Our main result is the following theorem.
For clarity we break the proof into several steps.

The harmonic lifting operator is bounded from above
We prove that H is bounded in the operator norm · * * induced by the energy norm (21). We refer to Appendix A for additional notation and auxiliary results used in the proof. We introduce the space where κ is the thickness of Ω o .
Proof. To prove (26) it suffices to show that For some positive constant C, inversely proportional to κ. According to the definitions of H and (·, ·) * in (25) and (21), this is equivalent to The strong Cauchy-Schwarz inequality for the harmonic component (see Lemma A.3) yields the following lower bound for the right hand side: We now proceed to bound v n (σ n ) 0,Ωo and v l (σ l ) 0,Ωo from below by v n (σ n ) 0,Ωn and v l (σ l ) 0,Ω l \Ωo , respectively. We start with the nonlocal term. Let µ n ∈ V η D (η n ) denote the extension of σ n by zero in η D , i.e., µ n = σ n in η c and 0 in η D . By using the nonlocal Poincaré inequality, the well-posedness of the nonlocal problem and Lemma A.
To analyze the local term we derive a Caccioppoli-type inequality for the local harmonic component. We where κ is the thickness of the overlap, see Fig. 2. These properties imply that gv l and g 2 v l belong to H 1 0 (Ω l ). Next, we note that v l is the solution of the weak formulation of (9) with f l = 0. Using g 2 v l as a test function then yields the following identity We use the latter to find a bound on ∇(gv l ) 0,Ω l : Thus, we conclude that (27) and (28)

The approximation error is bounded by the modeling error
The optimal solution (θ * n , θ * l ) of the reduced space problem (17) approximates the trace of the global nonlocal solution u n on η c and Γ c , respectively. The following lemma shows that the error in (θ * n , θ * l ) is bounded by the modeling error on Ω l . Lemma 4.2: Let u n and (θ * n , θ * l ) solve (8) and (17), respectively. Then, Proof. Because (θ * n , θ * l ) satisfies the Euler-Lagrange equation (20) we have Using this identity together with the energy norm definition (21) yields The lemma follows by setting (µ n , µ l ) = T ( u n )−(θ * n , θ * l ) above and observing that

Proof of Theorem 4.1
Let u := L(T ( u n )). Definitions (24) and (25) together with the identities u n | Ωn = v n (T n ( u n )) + u 0 n and u l = v l (T l ( u n )) + u 0 l , imply that Likewise, the identities u * n = v n (θ * n ) + u 0 n and u * l = v l (θ * l ) + u 0 l imply that u * = L(θ * n , θ * l ). Adding and subtracting u to the LtN error then yields The first term in (33) is the consistency error of L; (32) implies that To estimate the second term we use (26) and (29): Combining (33) with this bound and (34) gives which completes the proof.

Convergence of the modeling error
In this section we show that u n − u l 0,Ω l vanishes as ε → 0. Proof. By definition u l solves the boundary value problem and so, it is also a solution of the weak equation Let ψ ∈ H 1 0 (Ω l ) solve the dual problem Since u l − u n = 0 on Γ l one can set w = u l − u n in (35) to obtain where the third equality follows from the fact that f l is extended to zero in η and the limit follows from the result in [17, Section 5] 7 .

Approximation of the optimization-based LtN formulation
This section presents the discretization and the error analysis of the LtN formulation (14). Throughout this section we assume that η c and Γ c are polygonal domains; this assumption is not restrictive as those are virtual domains that we can define at our discretion.

Discretization
We use a reduced-space approach to solve the optimization-based LtN problem (14) numerically, i.e., we discretize and solve the problem where u n (θ n ) ∈ V η D (Ω n ) solves the weak nonlocal equation B n (u n (θ n ); z n , κ n ) :=

Ωn Ωn
Gu n · Gz n dydx + ηc u n κ n dx = ωn f n z n dx + ηc θ n κ n dx, for all (z n , κ n ) ∈ V ηn × Θ * n , and u l (θ l ) ∈ H 1 Γ D (Ω l ) solves the weak local equation for all (z l , κ l ) ∈ H 1 0 (Ω l ) × Θ * l . Here, Θ * n and Θ * l are the duals of Θ n and Θ l , respectively. To discretize (36)-(38) we consider the following conforming finite element spaces for the nonlocal and local states, controls, and test functions 8 , respectively. In general, the finite element spaces for the nonlocal and local problems can be defined on different meshes with parameters h n > 0 and h l > 0, respectively, and can have different polynomial orders given by integers p n ≥ 1 and p l ≥ 1, respectively. Restriction of (36)-(38) to the finite element spaces (39) defines the discrete reduced-space LtN formulation min where u h n (θ h n ) ∈ V h η D solves the discrete nonlocal state equation and u h l (θ h l ) ∈ H h Γ D solves the discrete local state equation 9 Following Section 3.1, we write the solutions of (41) and (42) as where v h n and v h l are "harmonic" components solving (41) and (42) with f n = 0 and f l = 0 respectively, whereas u h0 n and u h0 l are "homogeneous" components solving (41) and (42) with θ h n = 0 and θ h l = 0, respectively. In terms of these components The Euler-Lagrange equation of (40) has the form: where To prove the positivity of Q h , the arguments of Lemma 3.1 cannot be extended, as the identity principle does not hold for v h l . We use instead the discrete strong Cauchy-Schwarz inequality in Lemma A.4.
Lemma 5.1: The form Q h (·, ·) defines an inner product on Θ h n × Θ h l .

Proof. We prove that
The discrete strong Cauchy-Schwarz inequality (see Lemma A.4) then implies Since δ < 1 the left hand side in the above inequality is nonnegative. Thus, we must have that v h n (σ h n ) = 0 and v h l (σ h l ) = 0, which implies (σ h n , σ h l ) = (0, 0). Lemma A.5 proves that Θ h n × Θ h l is Hilbert with respect to the discrete energy norm This fact, Lemma 5.1 and the projection theorem provide the following corollary.
Corollary 5.1: The reduced space problem (40) has a unique minimizer.

Convergence analysis
In this section we prove that the discrete solution (θ h * n , θ h * l ) converges to the exact solution (θ * n , θ * l ) assuming the latter belongs to the "raw" control space Θ n × Θ l . This assumption mirrors the one made in [1] and is necessary because the continuous problem is well-posed in the completion of the raw control space. We prove this result under the following assumptions.
We will estimate the discrete energy norm of the error (θ * n −θ h * n ; θ * l −θ h * l ) using Strang's second 11 Lemma; see, e.g., [18,Lemma 2.25,p.94]. Application of this lemma is contingent upon two conditions: (i) the discrete form Q h is continuous and coercive with respect to · h * , and (ii) there exists a positive real constant C such that µ n , µ l h * ≤ C µ n , µ l * ∀ (µ n , µ l ) ∈ Θ n × Θ l .
The first assumption holds trivially. To verify (47) note that Given µ l ∈ Θ l the function v h l (µ l ) solves the weak equation Additionally, similarly to [1], we assume that there exist positive constants γ * n , γ * l , γ n * , and γ l * such that for h n and h l small enough the following inequalities hold: The latter, the strong Cauchy-Schwarz inequality and the boundedness of the L 2 projection operators yield Application of Strang's second lemma then yields the following error estimate.
Lemma 5.2: Let (θ * n , θ * l ) and (θ h * n , θ h * l ) be the solutions of (20) and (44), then We use the result in Lemma 5.2 to obtain asymptotic convergence rates under the assumption that 1) the homogeneous problems (19) have solutions u 0 n ∈ H pn+t (Ω n ), for t ∈ [0, 1], and u 0 l ∈ H p l +1 (Ω l ); 2) the control variables (θ n , θ l ) are such that u n (θ n ) ∈ H pn+t (Ω n ) and u l (θ l ) ∈ H p l +1 (Ω l ). We treat the first term in (50) by using the norm-equivalence (67); we have where · Θn×Θ l is defined as in (15). We focus on the second term in (50). Adding and subtracting Q(θ * n , θ * l ; µ h n , µ h l ) and using conformity of Θ h n × Θ h l gives . Adding and subtracting the terms to the last expression and using the definitions of Q, F , Q h and F h yields the identity: Application of the Cauchy-Schwartz inequality then gives the following upper bound: Furthermore, note that v h n (µ h n ) − v h l (µ h l ) 0,Ωo = µ h n , µ h l h * = 1, and that u * n − u * l 0,Ωo = J(u * n , u * l ) is the optimal value of the objective functional, which is bounded by the modeling error. The regularity assumptions on the nonlocal solutions in (19) allow us to apply Theorem 6. where t ∈ [0, 1]. Furthermore, the regularity assumptions on the local solutions in (19) allow us to use Corollary 1.122 in [18, p.66] to conclude that According to Weyl's Lemma [31] the local harmonic liftings v l (θ * l ) and v l (µ h l ) are smooth functions and so there are positive constants C 1 and While a similar result holds for the nonlocal harmonic lifting v n (θ * n ), the treatment of v n (µ h n ) is more involved, due to the discrete nature of the Dirichlet data, and it requires an auxiliary function µ n ∈ C ∞ (η c ) such that µ n − µ h n L 2 (ηc) ≤ , for an arbitrarily small . Because v n and v h n depend continuously on the data, Since can be arbitrarily small, the last two terms in (52) are negligible. To complete the estimate we only need a uniform bound on v n ( µ n ) pn+t,Ωn . To this end, assume that for all is also nonlocal harmonic for all β ≤ k. Thus, D β [v n ( µ n )] has a uniformly bounded L 2 norm, i.e. D β [v n ] 0,Ωn ≤ C β , ∀ β ≤ k. This implies the existence of a positive constant C such that, v n ( µ n ) pn+t,Ωn ≤ C. It follows that there exist positive constants C 1 and We have just shown the following result.
We use Theorem 5.1 to estimate the Θ n × Θ l norm of the discretization error.
(54) Proof. Adding and subtracting Π n θ * n and Π l θ * l , and using the triangle inequality Using standard finite element approximation results for the first term yields We focus on the second term in (55). By the norm-equivalence (67) in the discrete control space, we have This result along with (53) and (56) implies (54). Since u * n and u * l depend continuously on the data, (54), Corollary 5.2 implies that that is, the L 2 norm error of the state variables is of the same order as the L 2 × H 1 2 norm error of the controls.

Numerical tests
We present numerical tests with the new LtN formulation in one dimension, including a patch test, a convergence study and approximation of discontinuous solutions. Though preliminary, these results show the effectiveness of the coupling method, illustrate the theoretical results, and provide the basis for realistic simulations. In our examples we use an integrable kernel, γ i , satisfying assumptions (2) and (5) to illustrate theoretical results and a singular kernel, γ s , often used in the literature as an approximation of a peridynamic model for nonlocal mechanics. These kernels are given by The patch test This test uses the linear manufactured solution u n = u l = x, u n | η D = x, u l (1.75) = 1.75, f n = f l = 0. We expect the LtN formulation to recover this solution exactly, i.e., u h * n ≡ u * n and u h * n ≡ u * l . Figure 4 shows the optimal states u h * n and u h * l , computed with ε = 0.065 and h = 2 −7 , for γ i (left) and γ s (right). The LtN method recovers the exact solution to machine precision.

Convergence tests
We examine the convergence of finite element approximations with respect to the grid size h using the following manufactured solutions: Note that, for both kernels, the associated nonlocal operator is equavalent to the classical Laplacian for polynomials up to the third order. For examples M.1 and M.2 we compute the convergence rates and the L 2 norm of the errors for the nonlocal state, e(u * n ), the local state, e(u * l ), and the nonlocal control parameter, e(θ * n ). The results are reported in Tables 1 and 2 for γ i and in Tables 3 and 4 for γ s in correspondence of different values of interaction radius ε and grid size h. In Fig. 5 we also report the optimal discrete solutions.
Results in Tables 1 and 2 show optimal convergence for state and control variables. We note that according to [16] and FE convergence theory [18] this is the same rate as for the independent discretization of the nonlocal and local equations by piecewise linear elements.
Remark 6.1: The convergence analysis in Section 5.2 establishes a suboptimal convergence rate in the L 2 norm of the discretization error of the state variables as we lose half order of convergence. We believe that the bound in (57) is not sharp, in fact, additional numerical tests (with h = 2 −8 , . . . 2 −12 ) show that there is no convergence deterioration.
For the singular kernel γ s there are no theoretical convergence results; however, there is numerical evidence that piecewise linear approximations of (11) are second-order accurate; see [6]. Our numerical experiments in Tables 3 and 4 show that the optimization-based LtN solution converges at the same rate.

Recovery of singular features
The tests in this section are motivated by nonlocal mechanics applications and demonstrate the effectiveness of the coupling method in the presence of point forces and discontinuities. We use the following two manufactured solution examples:   A.2 u n | η D = 0, u l (1.75) = 0, Tab. 4: Example M.2 with γ s : dependence on the grid size h and interaction radius ε of the error.
In Fig. 6 we report the optimal discrete solutions for h = 2 −7 and ε = 0.065. In particular, A.2 is a significant example that shows the usefulness of the coupling method in approximating the true solution with a local model where the nonlocality effects are not pronounced, i.e. the solution is smooth.

A Ancillary results
This appendix contains several results necessary for the well-posedness of the continuous and discrete reduced space problems and for the estimate of H * * . In the following results and proofs we let C and C i , i = 1, 2, . . ., be generic positive constants.
Proof. We need the following subspaces of V (η) and V (Ω) Let χ Σ be the indicator of Σ. Definition (4) and the fact that σ vanishes on η \ τ imply that To bound the energy norm of χ Σ v note that where the last two inequalities follow from (6) and (7) respectively. Thus, with C 2 = (1 + 2γ 1 C 2 pn ).

B Notation summary
In this appendix we report a summary of the notation we use for local and nonlocal domains. In Table  5 we report local entities on the left and nonlocal entities on the right (see Fig. 2 for a two-dimensional configuration).