# Coupling of polynomial approximations with application to a boundary meshless method

Yendoubouam Tampango, Michel Potier-Ferry, Yao Koutsawa, Sonnou Tiem

► To cite this version:

Yendoubouam Tampango, Michel Potier-Ferry, Yao Koutsawa, Sonnou Tiem. Coupling of polynomial approximations with application to a boundary meshless method. International Journal for Numerical Methods in Engineering, 2013, 95 (13), pp.1094-1112. 10.1002/nme.4549 . hal-01500816

HAL Id: hal-01500816

https://hal.univ-lorraine.fr/hal-01500816

Submitted on 4 May 2017

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Coupling of polynomial approximations with application to a boundary meshless method

Yendoubouam Tampango<fnref n="1" />,<fnref n="3" />, Michel Potier-Ferry<fnref n="1" />,<fnref n="2" />, Yao Koutsawa<fnref n="3" /> and
Sonnou Tiem<fnref n="4" />

New bridging techniques are introduced to match high degree polynomials. This permits piecewise resolutions of elliptic partial differential equations (PDEs) in the framework of a boundary meshless method introduced recently. This new meshless method relies on the computation of Taylor series approximations deduced from the PDE, the shape functions being high degree polynomials. In this way, the PDE is solved quasi-exactly inside the subdomains so that only discretization of the boundary and the interfaces are needed, which leads to small size matricial problems. The bridging techniques are based on the introduction of Lagrange multipliers and a set of collocation points on the boundary and the interfaces. Several numerical applications establish that the method is robust and permits an exponential convergence with the degree.

KEY WORDS: Taylor series; meshless; domain decomposition; Arlequin; bridging

## 1. INTRODUCTION

In this paper, we improve a high-order boundary meshless method that has been recently introduced [1, 2]. The basic idea of this method is to reduce the number of unknowns by analytically solving the partial differential equation (PDE) using the technique of a Taylor series. More precisely, the solution is sought in the form of a high-order polynomial, and the number of its coefficients is reduced by vanishing the Taylor series of the PDE at a given expansion point. For instance, in the case of a second-order elliptic PDE in a 2D domain, the number of coefficients is reduced from 231 to 41 for a degree $p = 20$, which is much easier to handle numerically. Next, the boundary conditions are accounted for using a least-square collocation method [3], because the pure collocation proves to be numerically unstable. For convenience, this association of a Taylor series and boundary least-square collocation will be named Taylor meshless method (TMM). It has several advantages. First, the approximate solution often converges exponentially with the degree ($p$-convergence), as in the $p$-version of the finite element method [4, 5], which can lead to very accurate solutions. Second, because of the quasi-exact solving of the interior problem, one needs only to discretize the boundary, which leads to a true boundary meshless technique. Lastly, the computation cost is much smaller

<footnote n="1">¹ LEM3, Laboratoire d' Etude des Microstructures et de Mécanique des Matériaux, UMR CNRS 7239, Université de Lorraine, Ile du Saulcy, 57045 Metz, Cedex 01, France</footnote>

<footnote n="2">² Laboratory of Excellence on Design of Alloy Metals for low-mAss Structures (DAMAS), Université de Lorraine, Metz, France</footnote>

<footnote n="3">³ CRP Henri Tudor, 66, rue de Luxembourg, L-4221 Esch-sur-Alzette, Luxembourg</footnote>

<footnote n="4">⁴ Ecole Nationale Supérieure d'Ingénieurs de Lomé, Université de Lomé, BP 1515 Lomé, Togo</footnote>
than with the p-version of the finite element method [1] because one avoids the cost of integrations, and there are many fewer unknowns.

The TMM differs significantly from other meshless methods by the type of shape functions. The most popular ones are based on moving least-square approximations [6, 7], radial functions [8, 9], or fundamental solutions [10–12]. As in the method of fundamental solutions (MFS), the shape functions of TMM are built up from the differential operator, but in MFS, they result from an analytical solution solving, which limits the choice of the considered operators. In TMM, the use of a Taylor series makes it possible to obtain quasi-exactly the particular solution and the general solution of the associated homogeneous equation. Contrary to MFS and because the PDE is solved quasi-exactly through a Taylor series, the present method remains truly a boundary technique even for non-homogeneous and nonlinear equations [13].

Taylor series may diverge, and it is not always possible to solve a boundary value problem with a single Taylor series. That is why a piecewise approach was proposed in the original paper [2]. Thus, there is a need to develop numerical methods for connecting high degree polynomials by keeping the property of exponential convergence with the degree: this is the aim of this paper.

There are many numerical bridging methods in the literature, many of them being based on Lagrange multipliers. For instance, in domain decomposition methods such as finite element tearing and interconnecting (FETI) [14–16], the exact continuity of the main unknown is enforced: when the convergence is reached, the approximation is in conformity. Over the last few decades, several non-conforming bridging procedures have been proposed, especially mortar [17] and Arlequin methods [18, 19], which rely on the introduction and the discretization of Lagrange multipliers defined on a matching zone. Within the mortar method, the continuity is enforced only in some nodes of the interface, and this continuity is understood in a mean sense in the Arlequin case. Intuitively, such weak continuity requirements seem more relevant to couple quite different functions. Nowadays, it is possible to match for instance non-conforming meshes [20, 21], 3D continua and beams or shell models [21, 22], and molecular and continuum models [23, 24]. Apparently, the coupling between several meshless approximations has not yet been studied, and the exponential convergence property with the degree (p-convergence) has not been discussed. For this purpose, we propose some variants of mortar and Arlequin techniques that are consistent with the framework of collocation meshless methods.

The paper is organized as follows. In Section 2, we present the boundary meshless method using Taylor series approximations and some applications showing its accuracy and cases where the method diverges in a single domain. In Sections 3–5, some new coupling techniques are developed in order to build a piecewise resolution within the presented meshless method. The coupling techniques are developed so that continuity and consistency are preserved, in some sense, along the interfaces between subdomains. The proposed piecewise techniques are based on Lagrange multipliers.

## 2. THE TAYLOR MESHLESS METHOD IN A SINGLE DOMAIN

### 2.1. Statement of the TMM

Let $\Omega$ be a bounded domain in $\mathbb{R}^2$. We consider a PDE in $\Omega$ with Dirichlet or mixed boundary conditions. The aim of the TMM is to build some Taylor series approximations of the solution of the PDE. For this purpose, the PDE is solved in its strong form by using a polynomial approximation. This operation leads to a family of equations, which makes it possible to reduce the number of unknown coefficients. To illustrate this procedure, let us consider the Laplace equation, with Dirichlet boundary condition:

$$ \begin{cases} \Delta u = 0 & \text{in } \Omega \\ u = u^d & \text{on } \partial\Omega \end{cases} \tag{1} $$

The unknown function $u$ will be sought as a Taylor series $u^p$ developed at a given point $c = [c_x, c_y]$ of $\Omega$:
$$u^p(x, y) = \sum_{n=0}^{p} \sum_{i=0}^{n} u_{n-i, i} (x - c_x)^{n-i} (y - c_y)^i \quad (2)$$

To illustrate this resolution technique, let us consider an order $p = 3$ and a development point $c = [0,0]$. Then the basic approximation solution is a polynomial of degree 3 involving ten coefficients ($u_{ij}$):

$$u^p(x, y) = u_{00} + u_{10}x + u_{01}y + u_{20}x^2 + u_{11}xy + u_{02}y^2 + u_{30}x^3 + u_{21}x^2y + u_{12}xy^2 + u_{03}y^3 \quad (3)$$

Next, the residual of the Laplace operator is

$$R(x, y) = 2(u_{20} + u_{02}) + 2(3u_{30} + u_{12})x + 2(u_{21} + 3u_{03})y \quad (4)$$

and it vanishes if the following conditions are satisfied:

$$\left\{ \begin{array}{l} u_{20} + u_{02} = 0 \\ 3u_{30} + u_{12} = 0 \\ u_{21} + 3u_{03} = 0 \end{array} \right. \quad (5)$$

From these last equations, we can find three of the unknown coefficients as a function of the others. Therefore, the previous approximated form (3) can be rewritten in the following form:

$$u^p(x, y) = u_{00} + u_{10}x + u_{01}y + u_{20}(x^2 - y^2) + u_{11}xy + u_{30}(x^3 - 3xy^2) + u_{21}\left(x^2y - \frac{1}{3}y^3\right) \quad (6)$$

The number of unknown coefficients is then reduced to 7. This methodology can be easily extended to any degree $p$ and to PDE with variables coefficients: in this case, one vanishes the Taylor series of the residual up to order $p-2$. For any degree $p$, the number of unknown coefficients goes from $(p+1)(p+2)/2$ to $2p+1$, after resolution of the PDE in its strong form.

The remaining coefficients, $\{v\} = \{u_{00}, u_{10}, u_{01}, u_{20}, u_{11}, u_{30}, u_{21}\}$ in the case of the approximation (6), will be found by applying boundary conditions. For this purpose, a least-square method combined with collocation method is used [2]. This technique has been suggested by Zhang et al. [3]. A set of $M$ collocation points $\underline{x}_j$ is chosen on the boundary of the domain, and one minimizes the error between the approximate value and the given value of $u$ in these points. It comes to minimize the function

$$J(v) = \frac{1}{2} \sum_{j=1}^{M} \left| u^p(\underline{x}_j) - u^d(\underline{x}_j) \right|^2 \quad (7)$$

This minimization leads to a linear system $[K]\{v\} = \{F\}$, where $[K]$ is an invertible matrix. Solving this system gives the vector $\{v\}$ and therefore the approximate solution of the problem (1).

## 2.2. An application with p-convergence

Let us consider a Laplace problem in a disk:

$$\left\{ \begin{array}{ll} -\Delta u = 0 & \text{in } \Omega = \{(x, y)/x^2 + y^2 \le 1\} \\ u(\underline{x}) = \frac{x-x_0}{(x-x_0)^2+(y-y_0)^2} & \text{on } \partial\Omega. \end{array} \right. \quad (8)$$

The exact solution of this problem is $u(\underline{x}) = (x - x_0)/[(x - x_0)^2 + (y - y_0)^2]$. This solution has a singularity at $X_0 = [x_0, y_0]$.

The solution of this problem by using TMM has been previously discussed in past papers [2, 25]. For instance, in [25], it has been compared with the finite element method: while a given accuracy
is obtained by the TMM with 90 degrees of freedom (DOFs), the same quality of solution requires at least 5000 DOFs with the best finite elements.

The numerical solution obtained by TMM depends on two parameters: the degree $p$ of the polynomials and the number of collocation points $M$. It has been established in [1, 2] that this procedure is robust with respect to the number of collocation points provided that it is large enough. For completeness, the influence of the number of collocation points has been illustrated in Figure 1 in the case $X_0 = [2, 0.3]$ for three values of the degree $p = 5, p = 10$, and $p = 20$. The maximal error decreases with $M$ until an optimal number where it becomes stable, as represented by the plateau on the graph. This optimal number is a little higher than $2p$. To be safe, we shall choose $M = 4p$ as recommended in [2, 25].

The $p$-convergence is illustrated in Figure 1 and Table I. It can be seen that it is possible to obtain a very accurate result by increasing the degree and that one converges exponentially by increasing the degree. This means that the method naturally satisfies the consistency condition with an arbitrarily large degree of convergence, in the same way as in the $p$-version of finite elements [4]. If the solution is a polynomial, the method provides the exact values with an infinitesimal margin of error (see, for instance, Figure 1-a of [26]). This is not surprising for a method based on Taylor series. Of course, a little care is needed to obtain a robust procedure within this meshless framework. The rate of convergence depends on the studied problem: in this case, the convergence is faster if the singularity is far away from the domain.

## 2.3. An application with boundary layers

Let us now consider the following problem in a rectangle:

$$ egin{cases} -\Delta u + u = 0 & \text{in } \Omega = \{(x, y) / -30 \le x \le 30 \text{ and } 0 \le y \le 4\} \\ u(x, 0) = u(x, 4) = 0 \\ u(\pm 30, y) = \sin(\frac{\pi}{4}y) \end{cases} \quad (9) $$

<figure><img src="image_2.png" /><figcaption>Figure 1. The influence of the number of boundary collocation points on the convergence.</figcaption></figure>

<table><caption>Table I. Logarithm of the relative error for the Laplace problem (8) with two different positions of the singularity.</caption><thead><tr><th></th><th>p = 5</th><th>p = 10</th><th>p = 15</th><th>p = 20</th><th>p = 40</th><th>p = 60</th></tr></thead><tbody><tr><td>X<sub>0</sub> = [2, 0.3]</td><td>-1.8781</td><td>-3.3457</td><td>-4.9492</td><td>-6.4344</td><td>-12.5522</td><td>-14.9132</td></tr><tr><td>X<sub>0</sub> = [1.2, 0.3]</td><td>0.0529</td><td>-0.5184</td><td>-0.9204</td><td>-1.3623</td><td>-3.2096</td><td>-5.0606</td></tr></tbody></table>

Throughout the paper, the error is the difference between exact and approximate solutions, divided by the maximum value of the solution.
<figure><img src="image_3.png" /><figcaption>Figure 2. Exact and approximate solutions of problem (9) along the line y = 1.5, degree 15.</figcaption></figure>

The analytical solution of this problem is given by

$$u(x, y) = \frac{\cosh\left(x\sqrt{1 + \frac{2}{16}}\right)}{\cosh\left(30\sqrt{1 + \frac{2}{16}}\right)} \sin\left(\frac{1}{4}y\right) \quad (10)$$

This problem is solved by the TMM with a development point at $c = [0, 0]$ and various approximation degrees $p$. Figure 2 presents the analytical and exact solutions along the line $y = 1.5$ for the degree $p = 15$. The graph shows that the approximate solution coincides with the exact solution in the center of the domain, but when approaching the edge, the error increases. This is due to the effect of boundary layers that are pronounced in this case.

Nevertheless, one can obtain an accurate solution in this way (error less than $10^{-4}$) with a degree $p = 60$, but one may be asked whether a better approach would be a splitting of the domain into several subdomains.

## 2.4. An application with a divergence in the whole domain

We now propose to solve the following the Poisson problem:

$$\left\{
\begin{array}{ll}
-\Delta u = -\frac{4}{(x^2+y^2)^2} & \text{in } \Omega \\
u(\underline{x}) = u^d(\underline{x}) & \text{on } \Gamma_d \\
\frac{\partial u}{\partial n}(\underline{x}) = g(\underline{x}) & \text{on } \Gamma_n
\end{array}
\right.
\quad (11)$$

The domain is a crown ($\Omega = \{(x, y) | r_1^2 \le x^2 + y^2 \le r_2^2\}$, with $r_1 = 0.8$ and $r_2 = 1$) or a sector of this crown. The discretization for the application of the boundary conditions is made as shown in Figure 3.

The Dirichlet and Neumann data are chosen in such a way that the exact solution is

$$u(x, y) = \frac{1}{x^2 + y^2} \quad (12)$$

This solution presents a singularity at the center of the crown ($X_0 = [0, 0]$). Contrary to the previous problems, the singularity is included in the convex hull of the domain. For this problem, we have mixed boundary conditions. The boundary least-square cost function is rewritten in the form

$$J(v) = \frac{1}{2} \sum_{j=1}^{M} \alpha(|u^p(\underline{x}_j) - u^d(\underline{x}_j)|^2) + \beta \left( |\frac{\partial u^p}{\partial n}(\underline{x}_j) - g(\underline{x}_j)|^2 \right) \quad (13)$$
<figure><img src="image_4.png" /><figcaption>Figure 3. Discretization of the crown and the sector for the Poisson problem (11).</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Figure 4. Representation of the domain of convergence of series calculated by the TMM for a development point chosen at $[0, 1]$.</figcaption></figure>

This new formulation presents a contribution for Dirichlet-type conditions and a contribution for Neumann-type conditions. The weight of each contribution is managed by the constants $\alpha$ and $\beta$. Many investigations have been made to discuss the influence of these constants on the convergence. The results showed that they have little influence except for very large or very small values. So in our work, we fixed $\alpha = \beta = 1$.

This problem has been solved in the whole crown with TMM by choosing various development points and various degrees: the results were always completely wrong. This observed divergence is not a surprise. Indeed, suppose that the domain of convergence is a disk limited by the first singularity. Then it is not possible to cover the whole domain with a single disk avoiding the singularity $X_0 = [0, 0]$ (Figure 4).

Next, the problem (11) has been solved in a sector with Dirichlet boundary conditions on the arcs and Neumann conditions along the radii (Figure 3). The development point is chosen at the center of the sector. The obtained results are presented in Table II: one sees that the algorithm converges very rapidly with the degree if the angle is lower than $\pi/2$.

In summary, these three applications (Laplace problem, problem (9), and Poisson problem) have shown that the TMM is a very accurate method, but in some cases, the resolution in the whole domain converges slowly (problem (8)) or diverges (Poisson problem (11)). However, a good degree of accuracy is easily obtained on a part of the domain. Then in the case of a complex domain or of a solution with singularities near the domain, it is necessary to split the domain into several subdomains and to apply the TMM in each subdomain. A question then arises: how can the continuity of the built series at the interfaces between the subdomains be ensured?
<table><caption>Table II. The maximal error of the Poisson problem (11) with different domains (whole crown or a part of the crown).</caption><thead><tr><th>Domain Ω</th><th>log(maximal error) p = 10</th><th>log(maximal error) p = 20</th></tr></thead><tbody><tr><td>Crown</td><td>3.8028</td><td>6.6711</td></tr><tr><td>1/2 Crown</td><td>0.8691</td><td>2.3435</td></tr><tr><td>1/4 Crown</td><td>-2.1389</td><td>-3.0064</td></tr><tr><td>1/6 Crown</td><td>-3.6026</td><td>-6.5307</td></tr><tr><td>1/8 Crown</td><td>-4.6345</td><td>-9.3923</td></tr></tbody></table>

### 3. COUPLINGS BY COLLOCATION

#### 3.1. Generalities

Let us consider an elliptic equation in a domain $\Omega$

$$ \begin{cases} Lu = f & \text{in } \Omega \\ u(\underline{x}) = u^d(\underline{x}) & \text{on } \partial\Omega \end{cases} \tag{14} $$

where $L$ is a second-order differential operator. In order to make a piecewise resolution, an interface ($\Gamma'$) is introduced by splitting the domain $\Omega$ into two subdomains $\Omega^1$ and $\Omega^2$. For simplicity, we only present Dirichlet conditions, but other boundary conditions have been considered in previous papers [25, 26]. A scheme of the coupling method is shown in Figure 5. The coupling techniques proposed here can be easily extended for any number of subdomains.

Therefore, solving problem (14) is equivalent to solving the following problems:

$$ \begin{cases} Lu^1 = f & \text{in } \Omega^1 \\ u^1(\underline{x}) = u^d(\underline{x}) & \text{on } \Gamma^1 \end{cases} \quad \text{and} \quad \begin{cases} Lu^2 = f & \text{in } \Omega^2 \\ u^2(\underline{x}) = u^d(\underline{x}) & \text{on } \Gamma^2 \end{cases} \tag{15} $$

with transmission conditions along the interface:

$$ u^1(\underline{x}) = u^2(\underline{x}) \quad \text{and} \quad \frac{\partial u^1}{\partial n}(\underline{x}) = \frac{\partial u^2}{\partial n}(\underline{x}) \quad \text{on } \Gamma^r \tag{16} $$

where $\Gamma^1 = \partial\Omega^1 \cap \partial\Omega$, $\Gamma^2 = \partial\Omega^2 \cap \partial\Omega$ and $\Gamma^r = \partial\Omega^1 \cap \partial\Omega^2$.

With the initial formulation of the TMM, it is possible to build polynomial shape functions for problems (15). Indeed, vanishing the coefficients of the Taylor series of the residual of these problems leads to two families of shape functions $P_k^1(\underline{x})$ and $P_k^2(\underline{x})$ such that

$$ u^1(\underline{x}) = P_0^1 + \sum_k P_k^1(\underline{x}) v_k^{(1)} \quad \text{and} \quad u^2(\underline{x}) = P_0^2 + \sum_k P_k^2(\underline{x}) v_k^{(2)} \tag{17} $$

<figure><img src="image_6.png" /><figcaption>Figure 5. Representation of the continuous and discretized problems for the resolution with two subdomains $\Omega^1$ and $\Omega^2$. The discretization introduces two types of nodes: the nodes on boundary $\Gamma^1$ and $\Gamma^2$ where the boundary conditions are accounted for and the collocation nodes on interface $\Gamma^r$ where the transmission conditions are ensured.</figcaption></figure>
The variables {$v^{(1)}$} and {$v^{(2)}$} will be found by applying the boundary conditions and the transmission conditions at the interface (16).

The discretization of the domain then generates two types of nodes (Figure 5):

* nodes on the external boundaries $\Gamma^1$ and $\Gamma^2$ where boundary conditions will be applied;

* nodes at the interface $\Gamma^r$ where the transmission conditions will be applied by using a coupling technique.

The boundary conditions are applied by using a least-square method as shown in the previous section. Thus, a set of nodes is chosen on each sub-boundary:

* A set of $M_1$ collocation nodes is chosen on the boundary $\Gamma^1$. On these nodes, the boundary conditions on $\Gamma^1$ are applied.

* A set of $M_2$ collocation nodes is chosen on the boundary $\Gamma^2$. On these nodes, the boundary conditions on $\Gamma^2$ are applied.

Then the functions accounting for boundary conditions are

$$ 
\begin{aligned}
J_1(v^{(1)}) &= \frac{1}{2} \sum_{j_1=1}^{M_1} \left| P_0^1(\underline{x}_{j_1}) + \sum_{i=1}^{2p_1+1} P_i^1(\underline{x}_{j_1}) v_i^1 - u^d(\underline{x}_{j_1}) \right|^2 \\
J_2(v^{(2)}) &= \frac{1}{2} \sum_{j_2=1}^{M_2} \left| P_0^2(\underline{x}_{j_2}) + \sum_{i=1}^{2p_2+1} P_i^2(\underline{x}_{j_2}) v_i^2 - u^d(\underline{x}_{j_2}) \right|^2
\end{aligned}
\quad (18) $$

In this paper, we are interested in the formulation of numerical techniques to ensure transmission conditions.

## 3.2. Discretization of the interface problem by least-square collocation

In [1], it was proposed that transmission conditions be ensured using the least-square method. For this purpose, a set of $M_r$ collocation points is chosen on the interface, and the application of transmission conditions would be satisfied in a mean sense via the following coupling function:

$$ c(v^{(1)}, v^{(2)}) = \frac{1}{2} \sum_{j_r=1}^{M_r} \left[ \alpha |u^1(\underline{x}_{j_r}) - u^2(\underline{x}_{j_r})|^2 + \beta \left| \frac{\partial u^1}{\partial n}(\underline{x}_{j_r}) - \frac{\partial u^2}{\partial n}(\underline{x}_{j_r}) \right|^2 \right] \quad (19) $$

As for the boundary least-square function for mixed conditions, in the function $c$, two constants $\alpha$ and $\beta$ have been introduced. Studies have been carried out on the influence of these constants. The results have shown that we can fix $\alpha = \beta = 1$.

Then for the whole problem (14), search the variables {$v^{(1)}$} and {$v^{(2)}$} equivalent to minimizing the following function:

$$ J(v^{(1)}, v^{(2)}) = J_1(v^{(1)}) + J_2(v^{(2)}) + c(v^{(1)}, v^{(2)}) \quad (20) $$

The minimization of the function $J$ leads to a linear symmetric and invertible system $[K]\{V\} = \{F\}$ to solve for the vectors {$v^{(1)}$} and {$v^{(2)}$}.

Nevertheless, according to several numerical applications in [1, 2], this procedure does not seem very robust, and it seems difficult to obtain a very high accuracy by increasing the degree. That is why two alternative techniques will be presented in the next sections.

# 4. COUPLING BY LAGRANGE MULTIPLIERS

## 4.1. Statement of a coupling method by Lagrange multipliers

The basic idea of this new procedure is to consider the transmission conditions as constraints in an optimization problem. We minimize the function ($L = J_1 + J_2$) introduced to account for the boundary conditions. Two sets of collocation points were introduced as in Section 3:
* $M_1$ (respectively $M_2$) collocation nodes on $\Gamma^1$ (respectively $\Gamma^2$) for the boundary conditions.

* $M_r$ collocation nodes on $\Gamma^r$ for the transmission conditions.

Then the discretized version of the optimization problem can be formulated as

$$ \begin{array}{ll} \text{Minimize} & L(v^{(1)}, v^{(2)}) = J_1(v^{(1)}) + J_2(v^{(2)}) \\ \text{with constraints} & \left\{ \begin{array}{l} \begin{pmatrix} 1 \\ j \end{pmatrix} (v^{(1)}, v^{(2)}) = u^1(\underline{x}_j) - u^2(\underline{x}_j) = 0 \\ \begin{pmatrix} 2 \\ j \end{pmatrix} (v^{(1)}, v^{(2)}) = \frac{\partial u^1}{\partial n}(\underline{x}_j) - \frac{\partial u^2}{\partial n}(\underline{x}_j) = 0 \end{array} \right. \quad \forall \underline{x}_j \in \Gamma^r \end{array} \quad (21) $$

According to the 'Lagrange multiplier technique', $(v^{(1)}, v^{(2)})$ is the solution of the problem (21) if and only if there are a couple of vectors $(\lambda^1, \lambda^2)$ such that the 4-uplet $(v^{(1)}, v^{(2)}, \lambda^1, \lambda^2)$ is a stationary point of the function:

$$ J(v^{(1)}, v^{(2)}, \lambda^1, \lambda^2) = J_1(v^{(1)}) + J_2(v^{(2)}) + c(v^{(1)}, v^{(2)}, \lambda^1, \lambda^2) \quad (22) $$

where $c$ is the coupling function defined by

$$ c(v^{(1)}, v^{(2)}, \lambda^1, \lambda^2) = \sum_{j=1}^{M_r} \left[ \begin{pmatrix} \lambda_j^1 & 1 \\ 1 & j \end{pmatrix} (v^{(1)}, v^{(2)}) + \lambda_j^2 \begin{pmatrix} 2 \\ j \end{pmatrix} (v^{(1)}, v^{(2)}) \right] \quad (23) $$

Setting $\Lambda = \langle \lambda_1, \lambda_2 \rangle$, the stationary points of the function $J$ are solutions of a linear system in the form

$$ \begin{bmatrix} K_1 & 0 & C_1^T \\ 0 & K_2 & -C_2^T \\ C_1 & -C_2 & 0 \end{bmatrix} \begin{Bmatrix} v^1 \\ v^2 \\ \Lambda \end{Bmatrix} = \begin{Bmatrix} F_1 \\ F_2 \\ 0 \end{Bmatrix} \quad (24) $$

The continuity at the interface is enforced only at the collocation points: as the approximation is not conforming, there is no perfect continuity along the interface.

## 4.2. Application to Laplace equation in a disk

Let us consider the Laplace problem (8), with a singularity at $[1.5, 0]$. The domain is subdivided in two parts as shown in Figure 6. Each collocation point $(\underline{x}_j)$ on the interface is associated to two Lagrange multipliers $(\lambda_j^1$ and $\lambda_j^2)$.

The solving parameters are

* the approximation degrees in the subdomains $(p_1$ and $p_2)$,

* the number of collocation nodes on the boundary of each subdomain $(M_1$ and $M_2)$,

* the number of collocation nodes on the interface $(M_r)$.

We present the error on the inner circle of radius $r = 0.8$. This circle encounters the interface at $\theta = \pi/2$ and $\theta = 3\pi/2$. Figure 7 presents the analytical solution and the exact solution for approximation degree $p_1 = p_2 = 10$ and $M_r = 10$ collocation points, that is, 20 Lagrange multipliers. This figure shows a good correspondence between the exact and approximate solutions. Thus, the transmission conditions have been well accounted for. It is also possible to see that the piecewise resolution does not introduce additional errors: we find an error of $3 \times 10^{-3}$ as in the resolution in a single domain.

We are now interested in the influence of the number of coupling nodes $(M_r)$ on the quality of the solution. It defines the quality of the information transmitted along the interface. Indeed, it is
<figure><img src="image_7.png" /><figcaption>Figure 6. Domain decomposition of the Laplace problem (8) for a piecewise resolution in two subdomains.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Figure 7. Laplace equation in the unit disk (8). Solutions and error on the inner circle of radius 0.8, resolution made by the Lagrange multiplier technique, $p_1 = p_2 = 10$. There is a good accuracy at the interface ($\theta = \pi/2$ and $\theta = 3\pi/2$).</figcaption></figure>

directly related to the number of continuity equations written at the interface (21). Moreover, this number defines the number of DOFs ($N_\lambda$) introduced by the coupling.

Table III presents the error for several approximation degrees and for different collocation points $M_r$. This table shows that a minimum number of points of collocation are needed to obtain a good convergence. An optimal convergence is obtained for a number of Lagrange multipliers around $N_\lambda = 2p$. From this value, the convergence is stabilized by increasing the number of Lagrange multipliers. This proves the robustness of the Lagrange multiplier technique as a coupling method. Furthermore, the table shows a $p$-convergence as in the case of resolution in a single domain: the convergence is rapidly improved by increasing the degree.

The optimal number of boundary collocation points had been previously discussed [2]: a choice $M_1 = 2 \times p$, $M_2 = 2 \times p$ will be sufficient. Beyond this limit, the error is stable, which establishes the robustness of the procedure.

### 4.3. Coupling polynomials with different degrees

Here, we are interested in the behavior of the coupling technique in the case where the polynomials on both sides of the interface have different degrees. We consider the preceding problem by choosing two different degrees $p_1, p_2$ in the subdomains $\Omega^1$ and $\Omega^2$. The development points of the
<table><caption>Table III. Laplace problem (8) split in two subdomains.</caption><thead><tr><th>Degree</th><th><i>M</i><sub>1</sub> − <i>M</i><sub>r</sub> − <i>M</i><sub>2</sub></th><th><i>N</i><sub>λ</sub></th><th>log(maximal error)</th></tr></thead><tbody><tr><td rowspan="3"><i>p</i> = 4</td><td>10–2–10</td><td>4</td><td>−0.1592</td></tr><tr><td>10–4–10</td><td>8</td><td>−1.1002</td></tr><tr><td>10–6–10</td><td>12</td><td>−0.9847</td></tr><tr><td rowspan="6"><i>p</i> = 10</td><td>20–6–20</td><td>12</td><td>−1.7369</td></tr><tr><td>20–8–20</td><td>16</td><td>−2.4695</td></tr><tr><td>20–10–20</td><td>20</td><td>−2.6641</td></tr><tr><td>20–12–20</td><td>24</td><td>−2.5460</td></tr><tr><td>20–14–20</td><td>28</td><td>−2.5460</td></tr><tr><td>20–20–20</td><td>40</td><td>−2.5460</td></tr><tr><td rowspan="9"><i>p</i> = 20</td><td>40–12–40</td><td>24</td><td>−2.9640</td></tr><tr><td>20–16–20</td><td>32</td><td>−4.4692</td></tr><tr><td>20–18–20</td><td>36</td><td>−5.2678</td></tr><tr><td>40–20–40</td><td>40</td><td>−5.1621</td></tr><tr><td>40–22–40</td><td>44</td><td>−5.1009</td></tr><tr><td>40–24–40</td><td>48</td><td>−5.1008</td></tr><tr><td>20–30–20</td><td>60</td><td>−5.0993</td></tr></tbody></table>

Influence of the number of the Lagrange multipliers $N_λ$.

<table><caption>Table IV. Laplace problem (8) split in two subdomains.</caption><thead><tr><th><i>p</i><sub>1</sub></th><th><i>p</i><sub>2</sub></th><th><i>N</i><sub>λ</sub></th><th><i>M</i><sub>r</sub></th><th>log(maximal error)</th></tr></thead><tbody><tr><td>5</td><td>15</td><td>12</td><td>6</td><td>−2.17</td></tr><tr><td>7</td><td>15</td><td>16</td><td>8</td><td>−3.47</td></tr><tr><td>8</td><td>15</td><td>16</td><td>8</td><td>−3.29</td></tr><tr><td rowspan="3">9</td><td rowspan="3">15</td><td>16</td><td>8</td><td>−2.88</td></tr><tr><td>20</td><td>10</td><td>−3.79</td></tr><tr><td>24</td><td>12</td><td>−3.67</td></tr><tr><td rowspan="2">12</td><td rowspan="2">15</td><td>24</td><td>12</td><td>−4.55</td></tr><tr><td>28</td><td>14</td><td>−4.32</td></tr><tr><td rowspan="2">15</td><td rowspan="2">15</td><td>28</td><td>14</td><td>−4.42</td></tr><tr><td>32</td><td>16</td><td>−4.18</td></tr><tr><td colspan="4">Single zone <i>p</i> = 15</td><td>−2.70</td></tr></tbody></table>

Coupling polynomials with different degrees by Lagrange multipliers method.

approximate series are chosen at $c_1 = [-0.5, 0]$ for the subdomain on the left side and at $c_2 = [0, 0]$ for the right one (Figure 6). The second development point has been chosen far from the singularity.

The results obtained after this resolution are presented in Table IV. This shows that the convergence is good even in cases where the degrees of the polynomials on both sides of the interface are very different. In any case, a correct accuracy has been obtained, which establishes the robustness of the method. Especially, good results have been found with various degrees in the left part of the domain ($5 \le p_1 \le 15$). It is also not very sensitive with respect to the number of multipliers $N_λ$, provided that $N_λ$ is larger than $2p_1$.

A comparison with the solution obtained with a single domain resolution shows that the piecewise resolutions generally improve the convergence.

However, the subdomain resolutions introduce more DOFs than the single domain resolution. In a single domain resolution, for an approximation degree $p$, the number of DOFs is $2p + 1$. For the subdomain resolutions, there are (1) the DOFs of the TMM in each subdomain ($2p_i + 1$) and (2) the DOFs introduced by the matching, which, in this case, is the number of Lagrange multipliers $N_λ$.
### 4.4. Application to Poisson equation in a crown

To illustrate the interest of the piecewise resolution, we are now looking at the Poisson problem in a crown ($0.8 \le r \le 1$) (11) using TMM coupled with the Lagrange multiplier technique. To obtain a convergent algorithm, we split the crown into several sectors (Figure 8), and a resolution of the Poisson equation is made in each sector. In each sector, the development point of the approximate series ($c_i$) has been chosen in the middle of the outer arc.

Figure 9 presents the maximal error with respect to the number of subdomains. On each interface, we have chosen $p/2$ collocation points. This is equivalent to $p$ Lagrange multipliers, that is, $p$ transmission equations. This result shows that the convergence is improved by increasing the number of subdomains. The error is about $10^{-2}$ for 4 subdomains, and it decreases to $10^{-7}$ for 10 subdomains. This confirms that the speed of convergence is related to the size of the subdomains and to the distance between point of development and first singularity. Indeed, by splitting the domain, one builds several Taylor series with domains of convergence, which widely cover the subdomains.

The convergence with the degree is now investigated by fixing the number of subdomains to 6 and increasing the approximation degree. The results are presented in Figure 10. This shows that the convergence is improved by increasing the degree until a value $p = 20$, where it becomes stationary, but with a very good degree of accuracy. Such a plateau is due to numerical noise related to the accuracy limit of the computer.

In summary, we have built a coupling technique, which makes it possible to maintain the exponential convergence in the case of piecewise resolution. Nevertheless, in the proposed technique, the

<figure><img src="image_9.png" /><figcaption>Figure 8. Discretization of the crown for a piecewise resolution.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Figure 9. A piecewise resolution of the Poisson problem in the crown (11) with the Lagrange multiplier technique. Influence of the number of subdomains on the convergence with $p = 10$ in each subdomain.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>Figure 10. A piecewise resolution of the Poisson problem in the crown (11) with the Lagrange multiplier technique. p-convergence for six subdomains.</figcaption></figure>

number of DOFs introduced by the coupling is directly related to the number of collocation nodes chosen at the interface. For instance, with two subdomains and a common degree p, the numerical experiments suggest a number of multipliers about 2p. In the next part, an alternative coupling technique is presented with a view to reducing the number of multipliers.

## 5. COUPLING BY AN ARLEQUIN-TYPE METHOD

### 5.1. *Outline of Arlequin method*

The Arlequin method is a multi-scale technique that consists in the superposition of numerical models of different nature. It was initiated by Ben Dhia [18] in 1998. Since then, many studies have been carried out on this method to prove its accuracy and efficiency [19, 20, 22, 24].

The original formulation of the Arlequin method uses the weak form of the PDE. To keep the spirit of the TMM, we propose adapting the Arlequin formulation to a discrete distribution of nodes. First, we recall the original formulation of the Arlequin method, and next, we will present how it has been adapted to the TMM.

The Arlequin method can be applied to a problem with several subdomains. Here, for simplicity, let us consider a domain $\Omega$ split into two subdomains $\Omega^1$ and $\Omega^2$ (Figure 5). Let us assume that these two subdomains have a nonempty intersection $S = \Omega^1 \cap \Omega^2$. As the Arlequin method is a multi-scale method, it is assumed that the two areas are discretized differently. To illustrate, let us consider a Poisson problem:

$$ \begin{cases} -\Delta u = f & \text{in } \Omega \\ u = u^d & \text{on } \partial\Omega \end{cases} \qquad (25) $$

Let $W_1$ and $W_2$ be the space of kinematically admissible displacement fields associated respectively to $\Omega^1$ and $\Omega^2$. From the weak formulation of the PDE, the Arlequin method is formulated as follows:

$$ \begin{gathered} \text{Find}(u_1, u_2, \lambda) \in W_1 \times W_2 \times M \\ \int_{\Omega^1} \alpha_1 \nabla u_1 \cdot \nabla v_1 \, d\Omega + c(\lambda, v_1) = \int_{\Omega^1} \alpha_1 f \cdot v_1 \, d\Omega, \quad \forall v_1 \in W_1 \\ \int_{\Omega^2} \alpha_2 \nabla u_2 \cdot \nabla v_2 \, d\Omega - c(\lambda, v_2) = \int_{\Omega^2} \alpha_2 f \cdot v_2 \, d\Omega, \quad \forall v_2 \in W_2 \\ c(\mu, u_1 - u_2) = 0 \quad \forall \mu \in M \end{gathered} \qquad (26) $$

where the non-negative weighting coefficients $\alpha_1$ and $\alpha_2$ satisfy

$$ \begin{aligned} \alpha_1(x) + \alpha_2(x) &= 1 \quad \forall x \text{ in } \Omega \\ \alpha_i &= 0 \quad \text{in } \Omega \setminus \Omega^i \end{aligned} \qquad (27) $$
The matching operator $c(\lambda, \nu)$ is defined on the matching area $S$ and is often chosen as

$$c(\lambda, \mu) = \int_S (\lambda \cdot \mu + l^2 \nabla \lambda \cdot \nabla \mu) d\Omega, \quad \forall (\lambda, \mu) \in M^2 \quad (28)$$

where $M$ is the set of admissible functions on $S$ and $l$ is a characteristic length. This matching operator is referred as $H^1$ coupling. There are other ways to define this operator, for instance, the $L^2$ type operator. Bauman et al. [23] have presented a detailed description of the mathematical properties of the Arlequin method associated with various types of coupling operators.

## 5.2. A new coupling method by collocation

To keep the spirit of Arlequin method, the matching zone, which was a curve in Sections 3 and 4, is replaced by a surface as shown in Figure 11. Let us choose $M_r$ collocation nodes in this area. The new matching function inspired by the Arlequin method with respect to the discrete distribution of nodes is defined by

$$c(\lambda, u) = \sum_{i=1}^{M_r} \left( \lambda(\underline{x}_i) u(\underline{x}_i) + l^2 \nabla \lambda(\underline{x}_i) \nabla u(\underline{x}_i) \right) \quad (29)$$

where $\underline{x}_i$ are the collocation points chosen on the matching zone.

The function $\lambda(\underline{x})$ is discretized by radial functions, and this introduces discrete Lagrange multipliers $\lambda_i^1; \lambda_i^2$:

$$\lambda(\underline{x}) = \sum_{j=1}^{N_{\lambda}} (\lambda_j^1 \Phi_j^1(\underline{x}) + \lambda_j^2 \Phi_j^2(\underline{x}))$$

$$\Phi_j^1(\underline{x}) = \exp\left(\frac{-|\underline{x} - \underline{x}_j|^2}{d^2}\right) \quad \text{and} \quad \Phi_j^2(\underline{x}) = (\underline{x} - \underline{x}_j) \cdot \underline{n} \exp\left(\frac{-|\underline{x} - \underline{x}_j|^2}{d^2}\right) \quad (30)$$

where $d$ represents the radius of influence of the shape functions $\Phi_i$.

Thus, in the matching zone, there are two types of nodes (Figure 11):

* $M_r$ collocation points, where the transmission conditions are satisfied in a mean sense.

<figure><img src="image_12.png" /><figcaption>Figure 11. Discretization of the matching area. This zone is a surface intersecting the two subdomains. It is discretized by two types of points: the collocation points and the centers of the radial functions.</figcaption></figure>
* $N_λ$ points, which are related to discrete Lagrange multipliers (center of the radial functions). Around each of these nodes, the transmission conditions are ensured in an average way, on all the collocation points inside the disk of radius $d$.

Hence, solving the Laplace problem in the domain $\Omega$ is equivalent to finding $(v^{(1)}, v^{(2)}, \Lambda^1, \Lambda^1)$ such that the following function $L$ is stationary:

$$L(v^{(1)}, v^{(2)}, \Lambda^1, \Lambda^1) = J_1(v^{(1)}) + J_2(v^{(2)}) + c(\Lambda^1, \Lambda^2, u^1 - u^2) \quad (31)$$

where $u^i(\underline{x})$ is the approximate solution in the area $\Omega^i$:

$$u^i(\underline{x}) = P_0^i + \sum_k P_k^i(\underline{x}) v_k^{(i)} \quad (32)$$

and $J_i$ is the least-square function for boundary conditions:

$$J_i(v^i) = \frac{1}{2} \sum_{\underline{x}_j \in \partial\Omega^i} \left| P_0^i(\underline{x}_j) + \sum_{k=1}^{2p_i+1} P_k^i(\underline{x}_j) v_k^i - u^d(\underline{x}_j) \right|^2 \quad (33)$$

Setting $\Lambda = \langle \Lambda_1, \Lambda_2 \rangle$, the minimization of the function $L$ leads to a linear system in the same form as in Section 4 (24).

### 5.3. Application to Laplace equation in a disk

Let us consider the Laplace problem (8). Here, we solve this problem in a piecewise way by using the coupled technique TMM–Arlequin. The domain is split into two areas separated by the line $x = 0$. To conform with the Arlequin spirit, this line interface has been extended to a surface in which the transmission conditions will be applied. The chosen matching zone is a rectangle of length $2R$ and of width $2h$ (Figure 12), and we shall discuss the influence of this length $h$, as well as of the degree $p$ and the radius of influence $d$. In each subdomain, the TMM is applied, and the boundary is discretized as shown in Figure 12.

We first investigated the influence of the area of the coupling zone on the quality of convergence. The results are presented in Table V. These results are obtained with the following parameters: $p_1 = p_2 = p$, $M_1 = M_2 = 4p$, $M_r = 4p$, and $N_λ = p$. This problem leads to accurate results with an error about $10^{-3}$ for degree $p = 10$. The results show that the value of the parameter $h$ (width of the matching zone) has no effect on the accuracy of the method. The error is the same if the matching zone is a surface or a curve. Therefore, in the following, the matching zone will be reduced to a

<figure><img src="image_13.png" /><figcaption>Figure 12. Discretization of the unit disk in two subdomains by Arlequin-collocation method.</figcaption></figure>
<table><caption>Table V. Laplace problem (8) in two subdomains solved by the Arlequin-collocation technique.</caption><thead><tr><th><i>h</i></th><th><i>d</i></th><th>log(Emax)<i>p</i> = 10</th><th>log(Emax)<i>p</i> = 20</th></tr></thead><tbody><tr><td rowspan="4">0.3</td><td>0.5</td><td>-2.8173</td><td>-4.7552</td></tr><tr><td>0.4</td><td>-2.8339</td><td>-4.9001</td></tr><tr><td>0.2</td><td>-2.7004</td><td>-4.7628</td></tr><tr><td>0.1</td><td>-2.6117</td><td>-4.7522</td></tr><tr><td rowspan="4">0.2</td><td>0.5</td><td>-2.8173</td><td>-4.7552</td></tr><tr><td>0.4</td><td>-2.8339</td><td>-4.9001</td></tr><tr><td>0.2</td><td>-2.7004</td><td>-4.7628</td></tr><tr><td>0.1</td><td>-2.6117</td><td>-4.7522</td></tr><tr><td rowspan="4">0</td><td>0.5</td><td>-2.8173</td><td>-4.7552</td></tr><tr><td>0.4</td><td>-2.8339</td><td>-4.9001</td></tr><tr><td>0.2</td><td>-2.7004</td><td>-4.7628</td></tr><tr><td>0.1</td><td>-2.6117</td><td>-4.7522</td></tr></tbody></table>

Sensitivity of the algorithm to parameters *h*, *d*, and *p*.

line ($h = 0$). Table V shows also the effect of the parameter *d* on the convergence. Let us remember that the parameter *d* is the radius of zones of influence around each Lagrange multiplier node. The results show that this parameter does not have a significant influence on the convergence.

## 5.4. Coupling polynomials with different degrees by Arlequin method

Let us consider again the Laplace problem with a singularity at [1.5, 0]. Here, we solve this problem by using the coupling technique TMM–Arlequin. The results are presented in Table VI. These results are similar to those of Table IV obtained by using the Lagrange multiplier coupling technique. A comparison with the solution obtained by a single domain resolution shows that the piecewise resolution improves the convergence. Furthermore, the subdomain resolution makes it possible to use different approximations on each subdomain. The accuracy obtained in single resolution for degree *p* = 15 is recovered with two subdomains by making an approximation of degree *p* = 5 in the subdomain at the left side and an approximation of degree *p* = 15 in the subdomain near the singularity.

<table><caption>Table VI. Laplace problem (8).</caption><thead><tr><th><i>p</i><sub>1</sub></th><th><i>p</i><sub>2</sub></th><th><i>N</i><sub>λ</sub></th><th><i>M</i><sub>r</sub></th><th>log(maximal error)</th></tr></thead><tbody><tr><td>3</td><td>15</td><td>6</td><td>10</td><td>-1.76</td></tr><tr><td>4</td><td>15</td><td>6</td><td>10</td><td>-1.95</td></tr><tr><td>5</td><td>15</td><td>6</td><td>10</td><td>-2.13</td></tr><tr><td>6</td><td>15</td><td>6</td><td>10</td><td>-2.22</td></tr><tr><td>7</td><td>15</td><td>6</td><td>10</td><td>-2.20</td></tr><tr><td></td><td></td><td>7</td><td>15</td><td>-2.45</td></tr><tr><td>8</td><td>15</td><td>7</td><td>15</td><td>-2.35</td></tr><tr><td></td><td></td><td>8</td><td>15</td><td>-2.32</td></tr><tr><td>9</td><td>15</td><td>7</td><td>15</td><td>-2.10</td></tr><tr><td></td><td></td><td>9</td><td>15</td><td>-3.63</td></tr><tr><td>10</td><td>15</td><td>10</td><td>18</td><td>-3.97</td></tr><tr><td>12</td><td>15</td><td>12</td><td>20</td><td>-3.42</td></tr><tr><td></td><td></td><td>15</td><td>15</td><td>-4.37</td></tr><tr><td>15</td><td>15</td><td>15</td><td>20</td><td>-4.37</td></tr><tr><td></td><td></td><td>15</td><td>30</td><td>-4.37</td></tr><tr><td colspan="4">Single zone <i>p</i> = 15</td><td>-2.70</td></tr></tbody></table>

Coupling two different polynomials by Arlequin-collocation method.
However, a comparison of Tables IV and VI shows that, in the Arlequin method, the number of needed Lagrange multipliers $N_{\lambda}$ is much smaller than in the Lagrange multipliers technique to obtain about the same accuracy. Thus, the Arlequin method reduces the number of DOFs while keeping a good convergence.

### 5.5. Application of Arlequin method to the Poisson problem

The last application of the Arlequin method will concern the Poisson equation in a crown with a singularity at the center of the crown (11). In the second section, it was shown that this problem diverges when solving in the whole domain and that a good convergence is obtained on a part of the crown. Here, the crown is split into several subdomains (Figure 8), and the Arlequin-collocation method is used to ensure transmission conditions at the interfaces.

The results are similar to the previous results when using the Lagrange multipliers. Figure 13 presents the error with respect to the number of subdomains considered for a degree $p = 10$ (with $N_{\lambda} = M_r = 10$), and Figure 14 presents the error with respect to the approximation degree for six subdomains. From these two figures, one notices that the Arlequin-collocation method ensures the transmission conditions at the interface and that the $p$-convergence is obtained as for the same problem in a sector.

In summary, from these last sections about the application of the Arlequin-collocation method, we can say that the Arlequin-collocation technique ensures transmission conditions as well as the Lagrange multipliers method. However, the Arlequin-collocation method has the advantage to need fewer DOFs than the Lagrange multiplier method. With the Lagrange multipliers and a degree $p = 20$, we need about $N_{\lambda} = 2 \times 20$ DOFs for the interface, while $N_{\lambda} = 20$ will be sufficient with the Arlequin-collocation method. Thus, the latter method permits us to reduce the number of

<figure><img src="image_14.png" /><figcaption>Figure 13. A piecewise resolution of the Poisson problem in the crown (11) with the Arlequin-collocation technique. Influence of the number of subdomains on the convergence with $p = 10$ in each subdomain.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>Figure 14. A piecewise resolution of the Poisson problem in the crown (11) with the Arlequin-collocation technique. $p$-convergence for six subdomains.</figcaption></figure>
DOFs from about $p = 20$ (i.e., 102 instead of 122). We have just established a first bridging technique of Arlequin type in a meshless framework. It is efficient and permits the reduction of the number of DOFs for the multiplier, as compared with the method in Section 4. There are likely many alternatives to discretize the multiplier, especially by radial functions (cubic splines, etc).

## 6. CONCLUSION

Two bridging techniques have been proposed to link high-order polynomial approximations, in the framework of a new discretization method named TMM, in which the shape functions are deduced from the PDE via Taylor series expansions. Within TMM, the solution converges exponentially with the degree, and it was a challenge to recover the same exponential convergence in a piece-wise approach. The proposed bridging methods are based on Lagrange multipliers. In the first one, the transmission conditions are enforced at some collocation points, and in the second one, they are satisfied in a mean sense, in the same spirit as in the Arlequin method [18]. This second technique seems to require fewer DOFs. Meshless bridging techniques are associated with a meshless discretization of a PDE for the first time, to our knowledge.

With these bridging techniques, a significant progress has been made in view of a wide application of TMM, because in practice, one cannot hope to approximate the solution of any elliptic PDE with a single polynomial. The examples in this paper are quite simple, but one can refer to [1] for an application in a complex domain, and to [13] for an application to nonlinear elasticity. No case in a 3D-domain has yet been studied, but we do not anticipate problems specific to 3D. Another difficulty of TMM is the calculation of Taylor series that is easy in the simple examples presented in this paper but can be intricate for most of the PDEs of practical interest. In this respect, automatic differentiation [27] is a well-known tool to calculate derivatives, and we plan to use it in order to establish a user-friendly implementation procedure, the goal being to solve generic PDEs without explicitly implementing the calculations of series. Finally, let us underline the strong potential of TMM that is a true boundary meshless method and converges exponentially by avoiding any integration procedure.

## ACKNOWLEDGEMENT

The authors acknowledge the financial support from the Luxembourg AFR grant PHD 09-017.

## REFERENCES

1. Zézé DS. Calcul de fonctions de formes de haut degré par une technique de perturbation. Ph.D. Thesis, Université Paul Verlaine Metz, 2009.

2. Zézé DS, Potier-Ferry M, Damil N. A boundary meshless method with shape functions computed from the PDE. *Engineering Analysis with Boundary Elements* 2010; **34**:747–754.

3. Zhang X, Liu XH, Song KZ, Lu MW. Least-squares collocation meshless method. *International Journal for Numerical Methods in Engineering* 2001; **51**:1089–1100.

4. Babuska I, Szabo BA, Katz LN. The p-version of finite element method. *SIAM Journal of Numerical Analysis* 1981; **18**:515–545.

5. Schwab C, Suri M. The p and hp versions of the finite element method for problem with boundary layers. *Mathematics of Computation* 1996; **65**:1403–1429.

6. Belytschko T, Lu YY, Gu L. Element free Galerkin methods. *International Journal for Numerical Methods in Engineering* 1994; **37**:229–256.

7. Nayroles B, Touzot G, Villon P. Generalizing the finite element method: diffuse approximation and diffuse elements. *Computational Mechanics* 1992; **10**:307–318.

8. Kansa EJ. Multiquadrics—a scattered data approximation scheme with applications to computational fluid-dynamics. I. Surface approximations and partial derivative estimates. *Computers Mathematics with Applications* 1990; **19**: 127–145.

9. Kansa EJ. Multiquadrics—a scattered data approximation scheme with applications to computational fluid-dynamics. II. Solutions to parabolic, hyperbolic and elliptic partial differential equations. *Computers Mathematics with Applications* 1990; **19**:147–161.

10. Fairweather G, Karageorghis A. The method of fundamental solutions for elliptic boundary value problems. *Advances in Computational Mathematics* 1998; **9**:69–95.
11. Golberg MA. The method of fundamental solutions for Poisson's equation. *Engineering Analysis with Boundary Elements* 1995; **16**:205-213.

12. Golberg MA, Chen CS. The method of fundamental solutions for potential, Helmholtz and diffusion problems. In *Boundary Integral Methods—Numerical and mathematical Aspects*. WIT Press: Southampton, 1998; 103–176.

13. Tampango Y. Développement d'une méthode sans maillage utilisant les approximations de Taylor. Ph.D. Thesis, Université de Lorraine, 2012.

14. Farhat C, Roux FX. A method of finite element tearing and interconnecting and its parallel solution algorithm. *International Journal for Numerical Methods in Engineering* 1991; **32**:1205-1227.

15. Gosselet P, Rey C. Non-overlapping domain decomposition methods in structural mechanics. *Archives of Computational Methods in Engineering* 2006; **13**:515-572.

16. Stefanica D. Parallel FETI algorithms for mortars. *Applied Numerical Mathematics* 2005; **54**(2):266-279.

17. Bernardi C, Maday Y, Patera AT. A new nonconforming approach to domain decomposition: the mortar element method. In *Nonlinear "Partial Differential Equations and Their Applications"*, H. Brezis & J.-L. Lions (ed.). Pitman: Harlow, 1994; 13–51.

18. Ben Dhia H. Problèmes mécaniques multi-échelles: la méthode Arlequin. *Comptes Rendus de l'Académie des Sciences - Series IIB - Mechanics-Physics-Astronomy* 1998; **326**:899-904.

19. Ben Dhia H, Rateau G. Analyse mathématique de la méthode Arlequin mixte. *Comptes Rendus de l'Académie des Sciences - Series I - Mathematics* 2001; **332**:649-654.

20. Ben Dhia H, Jamond O. On the use of XFEM within the Arlequin framework for the simulation of crack propagation. *Computer Methods in Applied Mechanics and Engineering* 2010; **199**:1403–1414.

21. Prudhomme S, Chamoin L, Ben Dhia H, Bauman PT. An adaptive strategy for the control of modeling error in two-dimensional atomic-to-continuum coupling simulations. *Computer Methods in Applied Mechanics and Engineering* 2009; **198**:1887–1901.

22. Hu H, Belouettar S, Potier-Ferry M, Daya EM. Multi-scale modelling of sandwich structures using the Arlequin method Part I: linear modelling. *Finite Elements in Analysis and Design* 2008; **45**:37-51.

23. Bauman P, Ben Dhia H, Elkhodja N, Oden J, Prudhomme S. On the application of the Arlequin method to the coupling of particle and continuum models. *Computational Mechanics* 2008; **42**:511–530.

24. Prudhomme S, Ben Dhia H, Bauman PT, Elkhodja N, Oden JT. Computational analysis of modeling error for the coupling of particle and continuum models by the Arlequin method. *Computer Methods in Applied Mechanics and Engineering* 2008; **197**:3399–3409.

25. Tampango Y, Potier-Ferry M, Koutsawa Y, Belouettar S. Convergence analysis and detection of singularities within a boundary meshless method based on Taylor series. *Engineering Analysis with Boundary Elements* 2012; **36**: 1465–1472.

26. Tampango Y, Potier-Ferry M, Koutsawa Y, Belouettar S. A new meshless method using Taylor series to solve elasticity problems. *European Journal of Computational Mechanics* 2012; **21**(3-6):365–373.

27. Griewank A. *Evaluating derivatives: principles and techniques of algorithmic differentiation*. Frontiers in applied mathematics, Vol. 19. SIAM: Philadelphia, 2000.