Cross-points in the Dirichlet-Neumann method I: well-posedness and convergence issues

Cross-points in domain decomposition, i.e., points where more than two subdomains meet, have received substantial attention over the past years, since domain decomposition methods often need special attention in their definition at cross-points, in particular if the transmission conditions of the domain decomposition method contain derivatives, like in the Dirichlet-Neumann method. We study here for the first time the convergence of the Dirichlet-Neumann method at the continuous level in the presence of cross-points. We show that its iterates can be uniquely decomposed into two parts, an even symmetric part that converges geometrically, like when there are no cross-points present, and an odd symmetric part, which generates a singularity at the cross-point and is not convergent. We illustrate our analysis with numerical experiments.


Introduction
Domain decomposition methods for partial differential equations are naturally defined and analyzed at the continuous level, like the original overlapping Schwarz method from 1869 [1], and also the original Dirichlet-Neumann [2] and the Neumann-Neumann [3] method. Even FETI (Finite Element Tearing and Interconnect) was first presented at the continuous level in the original publication [4], as a minimization problem, before the authors proceeded to the finite element discretization that led to its name. For symmetric and positive definite problems, early domain decomposition research focused then however on condition number estimates for such methods used as preconditioners at the discrete level, leading to groundbreaking results (see, e.g., [5] and references therein). The methods were thus intimately linked with the conjugate gradient method and not considered as standalone solvers, in contrast to multigrid methods, for example [6,7]. For more general problems which are non-symmetric and/or indefinite, condition numbers are not the key quantities anymore for understanding their convergence when used as preconditioners, and directly studying preconditioning properties for more general Krylov methods becomes difficult. There has therefore been an effort to also investigate the underlying iterative versions of these methods, and the study of their convergence properties at the continuous level (see, for example, [8] and references therein). This reveals many interesting properties of domain decomposition methods which are masked by the Krylov method that can correct convergence problems when the domain decomposition method is used as preconditioner. An interesting example is the Additive Schwarz method, which needs Krylov acceleration to be used, while Restricted Additive Schwarz does not, since it corresponds directly to the discretization of the parallel Schwarz method of Lions [9], and thus converges as a standalone iterative method [10,11].
We are interested here in understanding the convergence properties of the Dirichlet-Neumann method in the presence of cross-points. Cross-points in domain decomposition methods have become a focus of attention over the past years because of an increasing interest in the domain decomposition research community to better understand the discretization of domain decomposition methods at cross-points, and the influence on the convergence of the iterative solution. For the Helmholtz equation and Després' seminal non-overlapping Schwarz method with Robin transmission conditions, cross-points do not hamper convergence [12], and the same holds more generally for non-overlapping optimized Schwarz methods when studied at the continuous level (see [13]). Care needs to be taken however when discretizing such methods (see, for example, [14][15][16][17][18]), and this is even more important when higher order transmission conditions like Ventcell conditions are used [19][20][21][22][23]. For the Neumann-Neumann method, a well-posedness issue has been identified in the presence of cross-points (see [24]), and the authors present a modification of the method to get around this difficulty. Also multitrace formulations pose problems at cross-points, and a solution involving specific non-local operators has been proposed in [25] (see also [26] for a purely algebraic formulation). There was even a dedicated mini-symposium on cross-points in domain decomposition methods at the last international conference on domain decomposition methods, DD27, in Hong Kong [27].
For the Dirichlet-Neumann method, the well-posedness issue in the presence of cross points was already mentioned in early work [28], but has so far not been fully analyzed for the iterative version of the algorithm. Most research works using this method or variants of this method do not encounter this problem. Indeed, many authors consider the case of domain decompositions with two subdomains [2,29,30] or many subdomains in stripes [31,32], which excludes the presence of cross-points.
Others use the Dirichlet-Neumann method as a preconditioner for a Krylov method (see, for example, [8,33,34]), which hides the problematic behavior of the Dirichlet-Neumann method. We focus here on a simple but instructive Laplace problem on a square divided into four squared subdomains of equal area. This allows us to develop a complete analytical understanding of the convergence of the Dirichlet-Neumann method in the presence of cross-points at the continuous level. We will show that the even symmetric part of the iteration converges like when no cross-points are present, but the odd symmetric part of the iteration generates singularities at cross-points and is thus not convergent. Even though our analysis is limited to the case of four subdomains, it reveals the behavior of the iterates locally at the cross-point. Therefore, it is representative of the behavior of the iterates near cross-points for more general domain decompositions of grid-shape. The study on how to correct the convergence problem of the odd symmetric part will be addressed in a further research paper.

Geometry and model problem
As shown in Fig. 1, we use as domain ⊂ R 2 for our Laplace model problem the square (−1, 1) × (−1, 1), divided into four non-overlapping square subdomains i , i ∈ I := {1, 2, 3, 4}, of equal area. With such a partition, there is one interior cross-point (red dot), and there are also four boundary cross-points (black dots). The presence of boundary cross-points has been identified as an obstacle for multitrace formulations, which requires some additional specific treatment (see, for example, [35,Section 2.4]). Conversely, these points are not an issue here for the Dirichlet-Neumann method. We denote the interfaces between adjacent subdomains by ij := int(∂ i ∩ ∂ j ), and the skeleton of the partition by := i,j ij , and we have ij = ji for all i, j . We further denote the interior of the intersection between ∂ i and the boundary ∂ by ∂ 0 i := int(∂ i ∩ ∂ ), and the interior of the left, right, bottom, and top sides of by ∂ l , ∂ r , ∂ b , and ∂ t . Thus, an arbitrary side of is denoted by ∂ σ , where σ is in the set of indices S := {l, r, b, t}. We use the same notation also for an arbitrary side of a subdomain i , namely ∂ i,σ with σ ∈ S.

Even and odd symmetric functions
We recall now the definition of an even/odd symmetric function in multivariate calculus, and a useful decomposition result, which we express here in terms of functions in L p , for p ∈ [1, +∞].
Definition 1 (Symmetric set) Let n ≥ 1 be an integer. A subset U of R n is said to be symmetric if, for any (x 1 , · · · , x n ) ∈ R n , Definition 2 (Even/odd symmetric functions) Let U ⊂ R n be open and symmetric, and p ∈ [1, +∞] Similarly, the function h is called odd symmetric if for almost all ( Proof Taking This implies that h e −h e =h o − h o a.e. in U . Since the left hand-side is even symmetric and the right hand-side is odd symmetric, we must have h e −h e =h o − h o = 0 a.e. in U , which contradicts the initial assumption. Definition 3 (Symmetry preservation of operators) Let X and Y be two Banach spaces such that X ⊂ L p (U ) and Y ⊂ L q (V ) for some open symmetric sets U , V and some p, q ∈ [1, +∞]. An operator T : X → Y is said to preserve symmetry if it satisfies the following properties: • for all h ∈ X such that h is even symmetric, T h is even symmetric, • for all h ∈ X such that h is odd symmetric, T h is odd symmetric.
The Laplace operator : H 2 ( ) → L 2 ( ) preserves symmetry, which can be proved using the standard chain rule for C 2 functions together with the density of C ∞ ( ) in H 2 ( ). In the same way, if p denotes an even symmetric function in W 1 2 ,∞ (∂ ), one can prove that the operator (∂ n + p) : H 3 2 (∂ ) → H 1 2 (∂ ) also preserves symmetry because for each x ∈ , n(−x) = −n(x).

Laplace model problem
We consider the Laplace problem with Dirichlet boundary condition on , that is: where f ∈ H −1 ( ) and g ∈ H 1 2 (∂ ). Of course, since is Lipschitz, it is known that (1) admits a unique solution u ∈ H 1 ( ). We also consider the Laplace problem with Robin boundary condition: find u solution to where p ∈ L ∞ (∂ ) (p ≥ 0 a.e. on ∂ ) is even symmetric, f ∈ H −1 ( ) and g ∈ H − 1 2 (∂ ). To ensure that the problem is well-posed in H 1 ( ), we assume that p is strictly positive on a subset of ∂ of non-zero measure.

Regularity results
In this section, we briefly recall some important results about the theory of elliptic boundary value problems in nonsmooth two-dimensional domains, adapted to our specific context of a square. The general results in arbitrary polygons, together with their proofs, can be found in [36,Chapter 4], [37], and [38]. For a brief review on the subject, the reader is also referred to the lecture notes [39,Chapter 2]. H 1 ( ) regularity As it has just been mentioned, the standard variational approach to study problems (1) and (2) leads to existence and uniqueness of solutions in H 1 ( ). However, for this result to hold, the regularity required on the boundary data is g ∈ H 1 2 (∂ ) for the Dirichlet case and g ∈ H − 1 2 (∂ ) for the Neumann (Robin) case. This means that g can be identified as the trace on ∂ of a function in H 1 ( ) in the Dirichlet case. Similarly, in the Neumann (Robin) case, it can be identified as the normal derivative on ∂ of a function in H 1 ( ). When is a polygon, the regularity on ∂ can be expressed by means of the regularity on each side of the polygon associated with so-called compatibility relations at the corners (see [37,Chapter 1]). Let us introduce the set of corners C which consists in four pairs of indices in S: For the reader's convenience, the restriction of g to a part of the boundary ∂ σ will be denoted by g σ := g | ∂ σ , for each σ ∈ S. With these notations, we have the following useful characterization for the Dirichlet case: g ∈ H 1 2 (∂ ) iff for some small ε > 0. There exists a similar (and simpler) characterization for the Robin case: Therefore, replacing the regularity conditions on the boundary data g in problems (1) and (2) by conditions (3) and (4) leads to well-posed formulations in H 1 ( ).

H 2 ( ) regularity
The existence and uniqueness results mentioned above only give us H 1 ( ) regularity for the solution, which means u might not even be continuous since ⊂ R 2 . Since we are interested in pointwise properties, especially what happens at the cross-point, we need more regular functions.

Theorem 2 (Dirichlet case) If in addition to the previous assumptions, we have
Theorem 3 (Robin case) If in addition to the previous assumptions, we have f ∈ L 2 ( ), p | ∂ σ ∈ C ∞ (∂ σ ), and g σ ∈ H 1 2 (∂ σ ) for all σ ∈ S, then the solution u to (2) is in H 2 ( ).
We will also need regularity results for the mixed (Dirichlet/Neumann) problem for subproblems generated by the Dirichlet-Neumann method: let D and N be two subsets of S corresponding to the indices for Dirichlet and Neumann boundary conditions. In order to avoid well-posedness issues, we assume that D = ∅. The mixed problem reads: find u solution to The set of corners C can be split into three subsets C = C D ∪ C N ∪ C M corresponding to Dirichlet corners, Neumann corners, or mixed corners, such that for each (σ, σ ) ∈ C,

Theorem 4 (Dirichlet/Neumann case) If in addition to the previous assumptions, we
Finally, let us state one last useful result about the regularity of such functions, which is a direct consequence of the Sobolev embedding theorem (see, for example, [40]).

Proposition 5 Let U be a bounded Lipschitz open subset of R 2 . Then, any function in
Proof From the Sobolev embedding theorem, since 2 > 3 2 and U ⊂ R 2 , we know that H 2 (U ) is (continuously) embedded in the Hölder space C 0,μ b (U ) for some μ ∈ (0, 1). In addition, any function in C 0,μ b (U ) is uniformly continuous on U ; therefore, it can be extended to a continuous function on U .
From now on, we will assume that the data satisfy the regularity required by the assumptions in Theorems 2 and 3. This ensures that the regularity of the solutions we deal with is at least H 2 ( ).

Even/odd symmetric decomposition
We now decompose the problems (1) and (2) into two subproblems, in order to analyze the subproblems separately. With Theorem 1, we can uniquely decompose the data, and thus the associated problem, into the even symmetric part and the odd symmetric part. For problem (1), this leads to: find u e and u o solutions to − u e = f e in , These subproblems are still well-posed and the solutions still have H 2 ( ) regularity. Using the symmetry preserving property of the operator : H 2 ( ) → L 2 ( ), we see that u e is even symmetric and u o is odd symmetric. Moreover, since (− ) is linear and the boundary conditions are linear as well, u e + u o solves (1). Therefore, u e + u o = u, and since the decomposition is unique, the solution to (6) u e coincides with u e , the even symmetric part of u, and similarly u o = u o . As we will see, the convergence analysis of the Dirichlet-Neumann method reveals different behaviors for the even and odd symmetric subproblems. Note that the geometric domain decomposition itself is also symmetric with respect to the origin (0, 0). As for the Dirichlet problem, we define the even and odd symmetric parts of (2), and still denote by u e and u o their solutions. When p is even symmetric, in addition to the operator , we have seen that the operator (∂ n + p) : H 3 2 (∂ ) → H 1 2 (∂ ) also preserves symmetry. Hence, the even and odd symmetry properties of the solutions u e and u o . Thus, the linearity of (2) yields u e + u o = u, and the uniqueness of the decomposition finally leads to u e = u e and u o = u o like in the Dirichlet case.

Analysis of the Dirichlet-Neumann method
We use as in [8, Section 1.4] a gray and white coloring (see Fig. 1) and define the sets of indices The transmission conditions are indicated in Fig. 1, where "D" stands for Dirichlet and "N" for Neumann. Given an initial guess u 0 and a relaxation parameter θ ∈ R, each iteration k ≥ 1 of the method can be split into two steps: To start the Dirichlet-Neumann method, we need an initial guess u 0 , or equivalently λ 0 := u 0 | , since only the traces of u 0 on the interfaces ij are used in the initialization step k = 1. This initial guess needs to satisfy the following compatibility condition.
Definition 4 (Compatible initial guess) An initial guess u 0 (or equivalently λ 0 ) is said to be compatible with the Dirichlet boundary condition if it satisfies u 0 | ∂ ∩ = g | (or equivalently λ 0 = g a.e. on ∂ ∩ representing the set of boundary crosspoints).
In what follows, we fix an initial guess λ 0 such that λ 0 ∈ C 0 ( ) ∩ H 3 2 ( ij ) for all (i, j ) such that ij = ∅, and λ 0 is compatible with the boundary condition.

Case of the even symmetric part
We begin with applying the Dirichlet-Neumann method to the even symmetric part of problem (1).

Theorem 6
Taking λ 0 e as the initial guess for the Dirichlet-Neumann method applied to the even symmetric part of problem (1) produces a sequence {u k e } k that converges geometrically to the solution u e in the L 2 -norm and the broken H 1 -norm, for any θ ∈ (0, 1). Moreover, the convergence factor is given by | 1 − 2θ |, which also proves that this method becomes a direct solver for the specific choice θ = 1 2 . 1 Proof We perform the first two iterations of the Dirichlet-Neumann method in terms of the local errors e k i := u i − u k i where u i := u | i is the restriction of the original solution to the i-th subdomain (see Section 4 (Example 1) for a numerical illustration).
Iteration k = 1, Dirichlet step: e is compatible with the even part of the Dirichlet boundary condition and u e | ∂ 1,σ ∈ H 3 2 (∂ 1,σ ) for all σ ∈ S, by Theorem 2, e 1 e,1 ∈ H 2 ( 1 ) exists and is unique. Moreover, we know that it is also continuous in 1 due to Proposition 5.
In the same way, we solve in 3 3 = u e − λ 0 e on 23 ∪ 34 . Since u e − λ 0 e is even symmetric, it follows that the only solution e 1 e,3 ∈ H 2 ( 3 ) to this problem verifies e 1 e,3 (x, y) = e 1 e,1 (−x, −y), for all (x, y) ∈ 3 . Iteration k = 1, Neumann step: Now, in 2 , we solve 23 . This problem is well-posed in H 1 ( 2 ). In addition, it is clear that the function defined in 2 by (x, y) → −e 1 e,1 (−x, y) solves the problem. By uniqueness, one deduces that the solution e 1 e,2 is given by e 1 e,2 (x, y) = −e 1 e,1 (−x, y), for all (x, y) ∈ 2 . Note that this equality extends to the whole 2 . Again, in the exact same way, we get that the solution e 1 e,4 in 4 is given by e 1 e,4 (x, y) = −e 1 e,1 (x, −y), for all (x, y) ∈ 4 . We are left with a recombined error e 1 e (defined in \ ) that is discontinuous across all parts of the skeleton where e 1 e,1 = 0 (see Fig. 3c). This may lead to discontinuities for the recombined solution u 1 e = u e + e 1 e . Also note that e 1 e is even symmetric in \ .
Iteration k = 2, Dirichlet step: The unique solution to this problem is e 2 e,1 = (1 − 2θ)e 1 e,1 . In the exact same way, we get in 3 3 .
The unique solution to this problem is e 2 e,2 = (1 − 2θ)e 1 e,2 . Again, in the exact same way, we get in 4 4 . For all θ ∈ (0, 1), the recombined solution e 2 e is exactly (1 − 2θ)e 1 e . Iterations k ≥ 3: From the analysis of the first two iterations, it follows by induction that, at iteration k, one has Fig. 3d). Therefore, the proposed domain decomposition method converges geometrically to the solution u e both in the L 2 -norm and the broken H 1 -norm for all θ ∈ (0, 1), Note that the recombined solution u k e is in general not continuous across the skeleton .

Remark 1
It is actually not difficult to find an initial guess λ 0 satisfying the assumptions of Theorem 6: the simplest way is to build a piecewise linear function on . If we denote P k , k ∈ {1, · · · , 4}, the four boundary cross points, and set λ 0 (P k ) := g(P k ) for each k, then λ 0 is compatible with the Dirichlet boundary condition. Next, we set λ 0 (0, 0) := 1 Of course, any other choice for λ 0 (0, 0) would work as well, but this specific choice leads to a function λ 0 with minimal slopes as it satisfies

Case of the odd symmetric part
Let us now turn to the odd symmetric part of problem (1). Given our initial guess, we obtain the following result. The previous result shows that in this case, at some point, the Dirichlet-Neumann method is not applicable anymore. In practice, since we are able to find a solution to this ill-posed problem, one may wonder if it is possible to recover a nice behavior of the method if we go past this iteration k 0 . The next theorem provides a negative answer. Before giving the proofs of these theorems, we need four technical lemmas which provide solutions to the Laplace problem on the two-dimensional cone C := R * + × (− π 4 , π 4 ) with different types of Dirichlet and Neumann boundary conditions on ∂C − := R * + × {− π 4 } and ∂C + := R * + × { π 4 } (see Fig. 2). Induction step. Let q ≥ 1 be fixed. We assume that the result of the lemma holds for all ≤ q. In order to find a solution to (8) for q + 1, we begin with the initial guess (ln r) q+1 , which obviously satisfies the boundary condition. First, we compute a correction r q+1 D such that (ln r) q+1 + r q+1 D is harmonic in C. Then, since the boundary conditions are no longer satisfied due to this first correction, we compute a second correctionr q+1 D such that (ln r) q+1 + r q+1 D +r q+1 D solves (8). For the first step, we start by computing the Laplacian in polar coordinates. Let us recall its expression for some smooth function f depending on (r, φ)

Lemma 9 For any integer
For our intial guess (ln r) q+1 , this leads to We can now cancel the term remaining on the right by adding on the left Again, to remove the new term on the right, we add on the left and continuing until the exponent on the (ln r) term in the right hand side reaches 0 or 1 gives ⎛ Therefore, introducing β q+1 j := (−1) j q+1 2j for all j , and defining the correction For the second step, let us note that this function verifies we obtain a function v Proof Let us proceed by induction, in the same way as in the proof of Lemma 9. Base case. It is not difficult to check that v 0 D := 4 π φ and v 1 D := 4 π φ ln r solve (9) for q = 0 and q = 1, respectively.
Induction step. Let q ≥ 1 be fixed. We assume that the result of the lemma holds for all ≤ q. As for Lemma 9, in order to find a solution to (9) for q + 1, we proceed in two steps, this time starting from 4 π φ(ln r) q+1 as initial guess. For the first step, following a similar iterative approach, we obtain a correction for all j . This enables us to get a function 4 π φ(ln r) q+1 + r q+1 D that is harmonic and is in on ∂C ± . This notation means that the sign depends on which part of the boundary is considered: + on ∂C + and − on ∂C − . Therefore, using the induction hypothesis, we introduce a second correctioñ Proof Again, we proceed by induction, in the same spirit as in the proofs of the previous lemmas. Base case Obviously, v 0 N := φ and v 1 N := φ ln r solve (10) for q = 0 and q = 1, respectively.
Induction step Let q ≥ 1 be fixed. The result of the lemma is assumed to hold for all ≤ q. As before, we aim at finding a solution to (10) for q + 1 by proceeding in two steps, this time starting from φ(ln r) q+1 as initial guess.
This initial guess is the same as in the proof of Lemma 10 (up to a factor 4 π ); therefore, we have for the first correction where β Proof We keep proceeding by induction. Base case Take q = 0. Then we have v 0 N := 2 π φ 2 − (ln r) 2 , which is indeed solution to (11) for q = 0.
Induction step Let q ≥ 1 be fixed. The result of the lemma is assumed to hold for all ≤ q. In order to find a solution to (11) for q + 1, we follow the usual two steps starting from the initial guess 2 π φ 2 (ln r) q+1 , which satisfies the boundary conditions. For the first step, we reuse the computations performed in the proof of Lemma 9, replacing q + 1 by q + 3. Rewriting the function (ln r) q+3 + r q+3 D , we get that is a harmonic function. We can easily deduce from this a first correction for our initial guess, q+3 2j for all j . As desired, the new function 2 π φ 2 (ln r) q+1 + r q+1 N is harmonic and belongs to L 2 (C) ∩ C ∞ C \ {(0, 0)} . In addition, it satisfies where the coefficients α q+1 N ,j,m depend on the β q+1 j and the α N ,j,m , for ≤ q −1.
We can now prove Theorem 7 and Theorem 8.

Proof of Theorem 7
We want to prove that the algorithm has a nice behavior up to some iteration k 0 , where the iterate becomes non-unique and singular near the crosspoint, with a leading singularity of type (ln r) 2 . The idea is to show that, in general, after two iterations only, the approximate solution given by the method is singular near the cross-point. In order to do so, let us apply the Dirichlet-Neumann method step by step to (7), and write the local subproblems in terms of the local errors e k o,i . Iteration k = 1, Dirichlet step: As for the even symmetric case, by Theorem 2, e 1 o,1 ∈ H 2 ( 1 ) ⊂ C 0 ( 1 ) exists and is unique. In the same way, we solve in 3 Since Here again, we know the problem is well-posed in H 1 ( 2 ). However, in contrast to the even symmetric case, we cannot argue that ±e 1 o,1 (−x, y) solves this problem. There is a sign incompatibility in the boundary condition: minus sign on 12 Since the boundary conditions are continuous on 12  In other words, whenever δ 1 = 0, the method enforces a discontinuous Dirichlet boundary condition at this step, which may lead to a singular solution since the compatibility relation (3) is no longer satisfied. Especially, the problem is not necessarily well-posed so the Dirichlet-Neumann method might not even be valid in this case.
Note that there is no other discontinuity enforced since the boundary conditions on 12 and 41 extend to 0 at 12 ∩ 1 and 41 ∩ 1 . Case δ 1 = 0 In what follows, we prove existence and uniqueness of e 2 o,1 , exhibiting the type of singularity induced by this nonsmooth boundary condition. In order to do so, we try to decompose e 2 o,1 as the sum of a regular part v 2 1 and a singular part w 2 1 , in the same spirit as in [37] for more regular problems. First, for each i ∈ I, let us introduce the angle φ i ∈ (0, 2π) such that the rotation R i of angle −φ i , given in polar coordinates by maps the quadrant containing i onto the cone C. More specifically, we have Using these notations, we define w 2 1 := (θ δ 1 ) · v 0 D • R 1 , whose expression in polar coordinates (r, φ) reads We know from Lemma 10 that w 2 1 is of class . Then, since we would like v 2 1 + w 2 1 to solve (12), we must define v 2 1 such that Note that the Dirichlet boundary condition enforced here is in C 0 (∂ 1 )∩H 1 2 (∂ 1,σ ) for all σ ∈ S. Therefore, v 2 1 ∈ H 1 ( 1 ) exists and is unique. Therefore, we have built (in a unique way) a solution v 2 1 + w 2 1 to (12) which is in Now, in order to conclude that e 2 o,1 = v 2 1 + w 2 1 , we must have uniqueness of the solution to (12) in L 2 ( 1 ). This uniqueness property is indeed guaranteed since we know from [37,Theorem 4.4.3.3] that the subspace of all solutions z ∈ L 2 ( 1 ) to − z = 0 in 1 , is of dimension 0. Hence, e 2 o,1 := v 2 1 + w 2 1 exists and is unique. In the same way, one can conclude that, in 3 , e 2 o,3 := v 2 3 + w 2 3 exists and is unique, with v 2 3 ∈ H 1 ( 3 ) and w 2 3 ∈ L 2 ( 3 ) \ H 1 ( 3 ). Of course, v 2 3 and w 2 3 can be obtained immediately from v 2 1 and w 2 1 using symmetry arguments, which gives for It is now clear that the algorithm generates a singular solution at this step. In order to estimate how the singularity propagates in the Neumann step, let us keep going and see what happens. Iteration k = 2, Neumann step: From Lemma 12, we know that w 2 2 ∈ C ∞ (R + × R − \ {(0, 0)}), and that it satisfies This time, in order for v 2 2 + w 2 2 to solve the problem for e 2 o,2 , we define v 2 2 such that Given the regularities of w 2 2 , v 2 1 , and v 2 3 , we deduce that v 2 2 exists and is unique in H 1 ( 2 ). Again, we have built (in a unique way) a solution v 2 2 + w 2 2 to (15), where v 2 2 ∈ H 1 ( 2 ) and w 2 2 ∈ L 2 ( 2 ) \ H 1 ( 2 ). But this is not enough to conclude that e 2 o,2 = v 2 2 + w 2 2 . As for 1 , we know that this last equality holds provided that the subspace of all solutions z ∈ L 2 ( 2 ) to is of dimension 0. Unfortunately, it follows from [37,Theorem 4.4.3.3] that its dimension is 1, and that it is spanned by a function in L 2 ( 2 ), say z 2 , that admits a singularity of type ln r at the cross-point. Therefore, one has that e 2 o,2 is not unique, and it can be written as e 2 o,2 = v 2 2 + w 2 2 + C 2 2 z 2 , for some constant C 2 2 ∈ R. However, no matter the value of C 2 2 , one may always deduce that, in a neighborhood V 0 of (0, 0), Note that there are actually three singular terms in e 2 o,2 : (φ − φ 2 ) 2 , (ln r) 2 and ln r. So the leading singularity is indeed (ln r) 2 .
Besides, one gets a similar result for e 2 o,4 . That is e 2 o,4 is not unique and can be written as and z 4 ∈ L 2 ( 4 ) admits a singularity of type ln r near the cross-point, Of course, one obtains a similar asymptotic result in the neighborhood V 0 of (0, 0), i.e., −θδ 1 8 π 2 (ln r) 2 in 4 ∩ V 0 . Note that in this case, the integer k 0 in the statement of the theorem equals 2.
Case δ 1 = 0 In this case, well-posedness is guaranteed and no singularity is generated by the method at the current iteration k = 2. Indeed, the method behaves exactly as in the first iteration and all e 2 o,i are well defined. Thus, we can introduce δ 2 ∈ R such that e 2 o,2 (0, 0) = −e 2 o,4 (0, 0) =: δ 2 . Then we are again facing two possible situations: δ 2 = 0, in which case well-posedness is lost and a (ln r) 2 singularity is generated at the next iteration k = 3 (i.e., k 0 = 3), or δ 2 = 0, in which case we still have the same behavior as in the first iteration. If we keep going with the same reasoning, we end up with two possible cases: either there exists some integer k 0 > 1 such that all iterates are uniquely defined and regular up to k 0 and both regularity and well-posedness are lost at k = k 0 , or δ k = 0 for all k ≥ 1, in which case all iterates are well defined and regular. This last case, which has never been encountered in practice, is not treated here as it seems difficult to study convergence properties in this specific situation.

Proof of Theorem 8
Here, we study how the singularity exhibited in Theorem 7 propagates through the next iterations k ≥ k 0 . To begin, we assume that the algorithm is capable of finding one of the solutions to a problem for which uniqueness is not guaranteed (typically problem (15)). To simplify notations, we introduce the integer p ≥ 0 such that k = k 0 + p. Then, we claim that at iteration k, for each i ∈ I, there exists a regular function v k i ∈ H 1 ( i ) and real coefficients γ k i,j,m j,m such that the local error e k o,i is given by In addition, for i ∈ I W , we have γ k i,2p,m = 0 if m > 1, which means that the leading singularity in i is of type (φ − φ i )(ln r) 2p . And for i ∈ I G , we have γ k i,2p+2,m = 0 if m = 0; thus, the leading singularity in i is of type (ln r) 2p+2 .
To prove this, we proceed by induction and use the results of the four lemmas stated earlier.
Base case The case k = k 0 , or equivalently p = 0, has already been seen in the proof of Theorem 7. Indeed, we have shown that the Dirichlet step in i (i ∈ I W ) led to a local error with a singular part w k 0 i that matches the one in (19) for k = k 0 (see expression (13) or (14)). In addition, the Neumann step in i (i ∈ I G ) led to a local error with a singular part w  (16) and (18) that this matches the singular part in (20) for k = k 0 .
Induction step Let k > k 0 , or equivalently p > 0, be fixed. Assuming our statement holds for any integer ≤ k, let us prove that it still holds for k + 1. Given e k o,i for each i ∈ I, the only way to get information about e k+1 o,i is to perform the Dirichlet and Neumann steps of the algorithm.
Dirichlet step: In 1 we solve where we have introduced for each j μ k 2,j := Note that we have performed some simplifications using the equality γ k 2,j,m = −γ k 4,j,m which holds for all k, j , m, due to the odd symmetry of the problem. Now, we consider each part of (22), and express the general form of its solution. First, since v k 2 and v k 4 are smooth and verify v k 2 (x, y) = −v k 4 (−x, −y) for almost all (x, y) ∈ 2 , we have already seen (in problem (12)) that there exists a solutioñ to part (a), where the coefficient in the linear combination depends on the jump v k 2 (0, 0) at the cross-point. Then, for each j , we get from Lemma 10 that there exists a solutioñ to part (b) j . In addition, we also get from Lemma 9 that there exists a solutioñ to part (c) j . Summing up all these contributions over the values of j , we deduce that there exists a functionṽ k+1  (21) that is in L 2 ( 1 ), which also proves that it is unique in L 2 ( 1 ). Moreover, using the previous decomposition forw k+1 1 and the induction hypothesis for e k o,1 , we get that there exists a function v k+1 with γ k+1 1,2p+2,m = θγ k+1 1,2p+2,m = 0 for m > 1. Obviously, the same result (existence, uniqueness, and decomposition) holds for e k+1 o,3 in 3 due to the odd symmetry property of the problem.
Neumann step: In 2 we solve We have already seen (again from [37,Theorem 4.4.3.3]) that, if there exists a solution to (23) that is in L 2 ( 2 ), then the set of all solutions is an affine subspace of dimension 1. As in the Dirichlet step, in order to prove existence, we decompose e k+1 o,2 using the decompositions obtained for e k+1 o,1 and e k+1 o, 3 . More specifically, we consider for 1 ≤ j ≤ 2p + 2 the sum of the solutions to the boundary value problems 1,j (ln r) j on 23 . . This time we have used the equality γ k+1 1,j,m = −γ k+1 3,j,m to simplify the formulations. Let us now study separately each part of (24), and express the general form of its solution. To begin with, due to the regularities of v k+1 1 and v k+1 3 , we know from standard results that there exists a (unique) solution to part (a). Then, for each j , we get from Lemma 12 that there exists a solution to part (b) j . Finally, we get from Lemma 10 that there exists a solution to part (c) j . Combining these results, we end up with a function v k+1 with C k+1 2 ∈ R. Since the singularity in z 2 is of type ln r, it follows that every solution can be written as (see (19) and (20)). The type of singularity generated by the domain decomposition method at iteration k is directly given by those expressions, so the proof of the result is complete. Remark 2 Following the same computations as in the previous proof while taking care of the coefficient in front of the leading singularity at each step enables us to end up with the following approximations of the errors in subdomains 1 , . . . , 4 where the correct sign is + for i ∈ {1, 2} and − for i ∈ {3, 4}.
The results given by Theorem 6, Theorem 7, and Theorem 8 can be easily extended to the Laplace problem with Robin boundary conditions. In this case, no compatibility condition is required, we only impose that the initial guess λ 0 satisfies the regularity assumption that λ 0 ∈ C 0 ( ) ∩ H 3 2 ( ij ) for all (i, j ) such that ij = ∅.

Theorem 13
Taking λ 0 e as the initial guess for the Dirichlet-Neumann method applied to the even symmetric part of problem (2) produces a sequence {u k e } k that converges geometrically to the solution u e with respect to the L 2 -norm and the broken H 1 -norm. Moreover, the convergence factor is given by | 1 − 2θ |, which also proves that this method becomes a direct solver for the specific choice θ = 1 2 . Proof The proofs can be obtained by following the same steps as in the proofs of Theorem 6, Theorem 7 and Theorem 8.

Remark 3
Note that the formulas given by the asymptotic analysis near the origin remain valid in the Robin case. Indeed, as this study has been conducted in the neighborhood of (0, 0), what happens "far away" from this point (e.g., on the boundary ∂ ) has no influence on the results.

Numerical experiments
We now illustrate the theoretical results obtained in the previous section. The square domain is discretized using a regular grid of size h, and our numerical method is based on a standard five-point finite difference scheme. In problems with mixed boundary conditions (Dirichlet and Robin), the Dirichlet boundary condition is enforced weakly using a penalty parameter ε of order 10 −10 . Unless otherwise stated, the mesh size will be set to h = 2 · 10 −2 in all experiments. In addition, each convergence analysis is performed taking u ex (exact solution) as the discrete solution obtained when solving on the whole domain with a direct solver.

Example 1
In order to illustrate the result of Theorem 6, i.e., convergence for the even symmetric part, we take the source term f = f e = 1 in , and set the Dirichlet boundary condition to g = g e = 0 on ∂ . A simple initial guess compatible with the Dirichlet boundary condition is in this case λ 0 = λ 0 e = 0 on . The results are displayed in Fig. 3. As expected, when θ = 1 2 , the DN method becomes a direct solver; thus, the error at iteration 2 is "zero" (here it cannot be smaller than the order of magnitude of ε) (see Fig. 3b). For θ = 1 2 , we see from Fig. 3c and d that the error on is multiplied by a constant from one iteration to the next. Moreover, the plot of the L 2 and broken H 1 norms of the error in Fig. 3e confirms that the DN method converges geometrically in this case. We also see that, as predicted by the proof of Theorem 6, at each iteration, the error remains continuous at the cross-point from 1 to 3 , and from 2 to 4 , and it has the expected symmetry property.

Example 2
We illustrate now the result of Theorem 13 for the even symmetric part with Robin conditions: the source term is f = f e = 1 in , and the Robin boundary condition is defined by p = p e = 1 and g = g e on ∂ , where g e is such that g e = 1 on ∂ b ∪ ∂ t (bottom and top sides) and g e = y 2 on ∂ l ∪ ∂ r (left and right sides). The same initial guess λ 0 = λ 0 e = 0 on is considered. As shown in Fig. 4, we observe the same convergence properties as for the Dirichlet problem (Example 1). Especially, the DN method becomes a direct solver when θ is set to 1 2 (see Fig. 4b), and for other choices of θ, it converges geometrically (see Fig. 4e), with the expected common ratio (1 − 2θ).

Example 3
Finally, we give an illustration of the problematic case described in Theorem 7 and Theorem 8. We consider a Dirichlet problem with the odd symmetric data: 2 4 sin(π x) cos( π 2 y) in and g = g o = 0 on ∂ . As in the previous examples, the initial guess is set to λ 0 = λ 0 o = 0. The results displayed in Fig. 5 show that, as expected, the DN method applied to this odd symmetric problem does not converge. More specifically, we see that after the first iteration (see Fig. 5b), continuity across the cross-point from 2 to 4 is already lost. This jump (referred to as δ 1 = 0 in the proof of Theorem 7) generates a "singularity" at the next iteration (see Fig. 5c). This "singularity" keeps propagating in the following iterations so that the method diverges, as predicted by Theorem 8. Moreover, the graph in Fig. 5d reveals a geometric (divergent) behavior of the errors in L 2 and broken H 1 norms. Fig. 4 Results for the DN method applied to the even symmetric part of (2) (Example 2) Remark 4 Since we use a standard finite difference scheme, it is not possible to enforce a discontinuous Dirichlet boundary condition. In practice, when two Dirichlet boundary values do not match at a corner, the average is computed and imposed at this corner. Thus, every numerical solution in i necessarily belongs to a finite dimensional subspace of C 0 ( i ), which explains the quotation marks for singularity in the previous paragraph. More generally, solving this problem using standard discretization methods (such as the finite difference method or the finite element method) involves a regularization step, namely the projection of the discontinuous boundary condition onto some finite dimensional subspace of C 0 ( i ).
Given our choice of projection (computing the average), the boundary data are regularized in the disk D h of radius h centered at the origin. Outside this disk, they are not modified. Consequently, one may argue that, in each subdomain i , the local numerical solution should be "not too far" from the local real solution in i \ D h . In order to verify this in the numerical results, we have plotted (red marks) the value of the error at the cross-point in subdomain 2 with respect to the mesh size h (see Fig. 6). We see that these points follow a curve of (ln h) 2 type. An interpretation of that result is that the discretization process (which acts as a regularization here) turns the singularity of type (ln r) 2 from Theorem 7 (see formula (17)) into a pseudo "singularity" of type (ln h) 2 . Now, we would like to analyze how the "singularity" propagates through the numerical iterates. In other words, given the "singularity" of type (ln h) 2 obtained at iteration k = k 0 , do we observe a "singularity" of type (ln h) 2p+2 at the following iterations k = k 0 + p with p > 0, as predicted by Theorem 8 ? To answer this, let us first note that the error at the cross-point | e k o,2 (0, 0) | seems to grow geometrically with respect to k (see Fig. 5d). Thus, for each h, we are able to compute constants α 2 , β 2 such that, for k ≥ 2, ln | e k o,2 (0, 0) | α 2 k + β 2 .
In addition, computing the logarithm of formula (25) (with i = 2), we get in the neighborhood of the cross-point ln | e k o,2 | ln θ 4 π ln r 2 k +β 2 , whereβ 2 depends on θ and k (only logarithmically). In Fig. 7, we have plotted the value computed for the coefficient α 2 as a function of the mesh size h, and tried to make it fit with a curve of type ln(c 1 (ln h) 2 ) (drawn in orange). As shown in the figure, this fitting was not successful and it appears that the appropriate fitting is a curve of type c 2 ln(c 3 (ln h) 2 ) (drawn in blue), with a constant c 2 0.87. This suggests that, in the numerical experiments, the "singularity" at iteration k = k 0 + p is of type (ln h) 1.74(p+1) rather than (ln h) 2(p+1) . One possible explanation for this is Fig. 7 Slope of the curve k → ln | e k o,2 (0, 0) | with respect to h, for the DN method applied to the odd symmetric problem described in Example 3, with θ = 0.45 the regularizing effect of the discretization, which may slow down the propagation of the "singularity."

Conclusion
We presented a complete analysis of the Dirichlet-Neumann method at the continuous level in a specific configuration involving one cross-point. Based on an even/odd symmetric decomposition of the data, we proved that the even symmetric part of the iterates converges geometrically to the right solution, while the boundary value problems associated to the odd symmetric part are not well-posed, which generates singular iterates. We also exhibited the type of singularity generated, and showed how this singularity propagates through the iterations. Finally, we studied the impact of our theoretical findings on numerical experiments.
A natural extension of this work would be to conduct a similar analysis for the Neumann-Neumann method, which is also known to pose problems in configurations with cross-points (see [41]). Another direction of future work, which will be the subject of the second part of this paper, is to build a modified (and convergent) Dirichlet-Neumann method taking advantage of this even/odd symmetric decomposition of the data.
Funding Open access funding provided by University of Geneva Data availability Data sharing not applicable to this article as no datasets were generated or analyzed during the current study.

Conflict of interest The authors declare no competing interests.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4. 0/.