Piecewise resolution

Keywords: Série de Taylor, méthode sans maillage, résolution par morceaux, différenciation automatique, fonctions de forme singulières, équations aux Dérivées partielles Taylor series, meshless method, boundary collocation, coupling techniques in piecewise resolution, automatic differentiation, singular shape function, partial differential equation -V - Taylor series, Meshless, Least-square, Lagrange multiplier, Piecewise resolution Taylor series, PDE, large-scale problem, Round-off errors. Chapter Taylor series, Multivariate algorithmic differentiation, Meshless, Newton Method

Based on Taylor Meshless Method (TMM), the aim of this thesis is to develop a simple, robust, efficient and accurate numerical method which is capable of solving large scale engineering problems and to provide a new idea for the follow-up study on meshless methods. To this end, the influence of the key factors in TMM has been studied by solving three-dimensional and non-linear Partial Differential Equations (PDEs). The main idea of TMM is to use high order polynomials as shape functions which are approximated solutions of the PDE and the discretization concerns only the boundary. To solve the unknown coefficients, boundary conditions are accounted by collocation procedures associated with least-square method. TMM that needs only boundary collocation without integration process, is a true meshless method.

The main contributions of this thesis are as following: 1) Based on TMM, a general and efficient algorithm has been developed for solving three-dimensional PDEs; 2) Three coupling techniques in piecewise resolutions have been discussed and tested in cases of large-scale problems, including least-square collocation method and two coupling methods based on Lagrange multipliers; 3) A general numerical method for solving non-linear PDEs has been proposed by combining Newton Method, TMM and Automatic Differentiation technique; 4) To apply TMM for solving problems with singularities, the singular solutions satisfying the control equation are introduced as complementary shape functions, which provides a theoretical basis for solving singular problems.

Résumé

Le but de cette thèse est de développer une méthode numérique simple, robuste, efficace et précise pour résoudre des problèmes d'ingénierie de grande taille à partir de la méthode Taylor Meshless (TMM) et fournir de nouvelles idée principale de TMM est d'utiliser comme fonctions de forme des polynômes d'ordre élevé qui sont des solutions approchées de l'EDP. Ainsi la discrétisation ne concerne que la frontière. Les coefficients de ces fonctions de forme sont obtenus en discrétisant les conditions aux limites par des procédures de collocation associées à la méthode des moindres carrés. TMM est alors une véritable méthode sans maillage sans processus d'intégration, les conditions aux limites étant obtenues par collocation.

Les principales contributions de cette thèse sont les suivantes: 1) Basé sur TMM, un algorithme général et efficace a été développé pour résoudre des EDP elliptiques tridimensionnelles; 2) Trois techniques de couplage pour des résolutions par morceaux ont été discutées dans des cas de problèmes à grande échelle: la méthode de collocation par les moindres carrés et deux méthodes de couplage basées sur les multiplicateurs de Lagrange; 3) Une méthode numérique générale pour résoudre les EDP non-linéaires a été proposée en combinant la méthode de Newton, la TMM et la technique de différentiation automatique. 4) Pour résoudre des problèmes avec un bord non régulier, des solutions singulières satisfaisant l'équation de contrôle sont introduites comme des fonctions de forme complémentaires, ce qui fournit une base théorique pour la résolution de problèmes singuliers.

Research background and significance

The physical quantity of the objective world is generally changed over time and space, and its intrinsic laws to be presented in the form of differential or partial differential equations (PDEs). There are many methods to solve these PDEs, for instance finite difference method (FDM), finite element method (FEM), boundary element method (BEM) and meshless methods.

FEM [START_REF] Zienkiewicz | The finite element method (Fifth edition) Volume 1: The Basis[END_REF]2] has been widely used in the field of engineering due to its robustness and universality. It contains the following advantages: 1) the use of equivalent integral weak form equations reduces the requirement of continuity of the interpolation functions; 2) the use of localized interpolation form a narrow bandwidth sparse stiffness matrix and this improves the stability of the calculation; 3) since the Kronecker delta property of the interpolation functions, it is convenient to account for essential boundary conditions.

There also exist the following drawbacks in FEM: 1) creation of a mesh for a complicated domain could be time-consuming; 2) when handling large deformation, considerable accuracy is lost because of the elements distortion; 3) in stress calculations, the stresses obtained by using FEM are discontinuous and less accurate. 4) the re-meshing technique may solve the problem of mesh distortion, however this will reduce the computational efficiency and accuracy; 5) it is very difficult to simulate the breakage of material into a large number of fragments as FEM is essentially based on continuum mechanics.

The main idea of BEM [START_REF] Brebbia | The boundary element method for engineers[END_REF] is to transform the PDE into boundary integral equation by employing the fundamental solutions and weighted residual approach. Since the boundary integral equation concerns the boundary, only the boundary needs to be discretized. The main advantages of BEM are as follows: 1) since no integration inside the domain is needed, the dimension and the number of degrees of freedom is strongly reduced; 2) as the fundamental solutions adopted in BEM satisfying the boundary condition at infinite, it is convenient to solve problems with unbounded domain.

There also exist some drawbacks when using BEM: 1) it is difficult to obtain the fundamental solutions for a complicated PDE and to deal with integration of singular fields; 2) when computing the crack propagation, one needs to repeatedly update the boundary grids; 3) normally the global matrix is unsymmetrical and full, which reduces the accuracy and robustness; 4) when dealing with non-linear cases, the integration inside the domain is required.

The main idea of FDM [4] is to use finite differences instead of derivatives. As the advantages of simplicity, flexibility and versatility of FDM, it is been widely used in fields of solid and fluid mechanics. It suffers from a major disadvantage in that it relies on regularly distributed nodes.

These difficulties associated with FEM, BEM and FDM mainly come from the mesh or grid in which a predefined connection between neighbor points is required. Thus the idea of eliminating the elements has evolved naturally. The concept of meshless or mesh free methods has been proposed, in which the domain of the problem is represented by a set of arbitrarily distributed nodes. The meshless framework not only provides a great convenience of pre-processing work, but also can avoid problems with mesh, such as mesh distortion, crack propagation, high velocity impact or explosive mechanics. It can effectively compensate the drawbacks of methods based on a mesh.

The meshless methods consists of two main steps: the approximation of unknown functions and the discretization of the PDE. The latter step has two main categories: Galerkin-based technique and collocation approach. Since background mesh and integration are required in Galerkin-based meshless method, it leads to expensive computational cost. The collocation-based meshless method is sometimes efficient since no integration is needed, but it is difficult to solve large-scale problems due to ill-conditioned matrices.

This thesis aims to discuss a newly proposed collocation-based meshless technique, named Taylor Meshless Method (TMM), for solving threedimensional non-linear PDEs. The final goal is to develop a simple, robust, efficient and accurate numerical method which is capable of solving large scale engineering problems and to provide a new idea for the follow-up study on meshless methods. To this end, the effect of matrix ill-conditioning and the propagation of round-off errors are analyzed carefully.

Review of meshless method

Mesh free or meshless method is defined with respect to the word "mesh". It is a common name of discretization that is different in different numerical methods, for instance, it is called grids in FDM, volumes or cells in FVM and named elements in FEM. The grids, volumes or cells and elements are collectively referred to as meshes since they are used to predefine the connection between nodes. To get rid of the tedious work of meshing and at the same time to avoid the computational difficulty caused by mesh distortion, a class of numerical methods without mesh based on interpolation of scatter points or based on shape functions fitting came into being, is collectively referred to as the meshless or mesh free methods.

Classification

The Smoothed Particle Hydrodynamics (SPH) proposed by Gingold and Monaghan [5] and Lucy et al. [START_REF] Lucy | A numerical approach to the testing of the fission hypothesis[END_REF] in 1977 is considered to be one of the earliest meshfree methods in literature. In the 1990s, a new class of meshfree methods emerged based on Galerkin method. Among which, the first one called diffuse element method (DEM) was proposed by Nayroles et al. [START_REF] Nayroles | Generalizing the finite element method: Diffuse approximation and diffuse elements[END_REF]. Thereafter in the framework of DEM, Belytschko made some improvements and proposed the Element Free Galerkin Method [START_REF] Belytschko | Element-free Galerkin methods[END_REF] (EFGM). Then in the following few decades, a variety of new meshless methods have sprung up.

According to the difference of discrete areas, the meshless method can be divided into two categories: the domain-type meshless methods and the boundary-type meshless methods.

Domain-type meshless methods

There are two main steps in domain-type meshless methods: 1) approximation of unknown functions; 2) discretization of control equation.

The first step is realized by using the interpolation of the arbitrary and irregular scattered points in the whole domain. There exist the following approaches to approximate the unknown functions in the literature: kernel particle approximation [5,[START_REF] Lucy | A numerical approach to the testing of the fission hypothesis[END_REF], reproducing kernel particle [START_REF] Liu | Reproducing kernel particle methods[END_REF], moving leastsquare [START_REF] Lancaster | Surfaces generated by moving least squares methods[END_REF][START_REF] Mclain | Drawing contours from arbitrary data points[END_REF][START_REF] Gordon | Shepard's method of "Metric Interpolatio" to bivariate and multivariate interpolation[END_REF], partition of unity [START_REF] Duarte | Hp clouds: A meshless method to solve boundary-value problems[END_REF], radial basis function [START_REF] Hardy | Multiquadric equations of topography and other irregular surfaces[END_REF][START_REF] Franke | Scattered data interpolation: tests of some methods[END_REF][START_REF] Micchelli | Interpolation of scattered data: distance matrices and conditionally positive definite functions[END_REF][START_REF] Buhmann | Multivariable interpolation using radial basis functions[END_REF][START_REF] Hardy | Theory and applications of the multiquadric-biharmonic method 20 years of discovery[END_REF], point interpolation method [START_REF] Liu | A local radial point interpolation method (lrpim) for free vibration analyses of 2-d solids[END_REF], etc.

Concerning the way to discretize the control equation, there are two major techniques: Galerkin method and collocation method. As for Galerkin-based methods, there are two main types of Galerkin methods. The first one is based on background integration in the whole domain while the second one named local Petrov-Galerkin method takes account of the integration in a rather small local sub-domain and no background mesh is required. Therefore meshless methods based on local Petrov-Galerkin integration are considered to be truly mesh free methods. Galerkin-based meshless methods benefit from their robustness and versatility that ensure the ability of solving large-scale problems. However the existence of integration may lead to expansive computational costs.

An alternative way to discretize the control equation is collocation technique. No background mesh and no integration is required, that makes it very efficient. Generally we believe that the collocation-based meshless methods are truly integration-free meshless method. However, when increasing the scale of considered problems, the ill-conditioned matrices in collocation-based methods may lead to numerical instability and low accuracy. By combining the above-mentioned ways to approximate the unknown functions with the approaches to discretize control equation, one can obtain a variety of domain-type meshless methods. The main domain-type meshless methods have been collected in Table 1.1. Additional information about meshless methods can be found in several review papers and books, see for instance [START_REF] Belytschko | Meshless methods: an overview and recent developments[END_REF][START_REF] Liu | Meshfree methods: moving beyond the finite element method[END_REF]. Table 1.1: Main domain-type meshless methods.

Galerkin

Local Petrov-Galerkin Collocation KP / / SPH [5,[START_REF] Lucy | A numerical approach to the testing of the fission hypothesis[END_REF] RKP RKPM [9] , MRKPM [START_REF] Liu | Wavelet and multiple scale reproducing kernel methods[END_REF] MLPG [START_REF] Shen | The meshless local petrov-galerkin (mlpg) method: A simple & less-costly alternative to the finite element and boundary element methods[END_REF] PCM [START_REF] Aluru | A point collocation method based on reproducing kernel approximations[END_REF] MLS DEM [7] , EFGM [START_REF] Belytschko | Element-free Galerkin methods[END_REF] MLPG [START_REF] Shen | The meshless local petrov-galerkin (mlpg) method: A simple & less-costly alternative to the finite element and boundary element methods[END_REF][START_REF] Atluri | A new meshless local Petrov-Galerkin (MLPG) approach[END_REF] FPM [26] , LSCM [START_REF] Zhang | Least-squares collocation meshless method[END_REF] RBF MG-RBF [START_REF] Wendland | Meshless Galerkin methods using radial basis functions[END_REF] MLPG [START_REF] Shen | The meshless local petrov-galerkin (mlpg) method: A simple & less-costly alternative to the finite element and boundary element methods[END_REF] RBF [START_REF] Kansa | Multiquadrics-a scattered data approximation scheme with applications to computational fluid-dynamics-I surface approximations and partial derivative estimates[END_REF][START_REF] Kansa | Multiquadrics-a scattered data approximation scheme with applications to computational fluid-dynamics-II solutions to parabolic, hyperbolic and elliptic partial differential equations[END_REF] , BKM [START_REF] Chen | A meshless, integration-free, and boundary-only RBF technique[END_REF] PU Hp Clouds [START_REF] Duarte | Hp clouds: A meshless method to solve boundary-value problems[END_REF] MLPG [START_REF] Shen | The meshless local petrov-galerkin (mlpg) method: A simple & less-costly alternative to the finite element and boundary element methods[END_REF] Hp-meshless clouds [START_REF] Liszka | hp-meshless cloud method[END_REF] PI PIM [START_REF] Liu | Meshfree methods: moving beyond the finite element method[END_REF] LPIM [START_REF] Liu | A local point interpolation method for stress analysis of two-dimensional solids[END_REF] /

Boundary-type meshless methods

The boundary-type meshless methods are based on complete families of shape functions that are fully exact solutions of considered problems. They are also known as Trefftz methods and many information about Trefftz methods can be found in several review papers, see for instance [START_REF] Herrera | Trefftz Method[END_REF][START_REF] Kita | Trefftz method: an overview[END_REF][START_REF] Jirousek | Survey of trefftz-type element formulations[END_REF][START_REF] Qin | Formulation of hybrid trefftz finite element method for elastoplasticity[END_REF] or in some books [START_REF] Li | Trefftz and collocation methods[END_REF][START_REF] Kołodziej | Boundary collocation techniques and their application in engineering[END_REF]. Since the PDEs are automatically satisfied, only the discritization of boundary is needed. There are two types of boundary-type meshless method: boundary integration and boundary collocation. The first one is more or less similar to BEM that is based on the boundary integral equation. Many such kinds of meshless methods can be found in the literature, see for instance Boundary Node Method [START_REF] Mukherjee | The boundary node method for potential problems[END_REF] (BNM), Boundary Element-Free Method [START_REF] Liew | Boundary element-free method (BEFM) and its application to two-dimensional elasticity problems[END_REF] (BEFM), Boundary Point Interpolation Method [START_REF] Liu | Boundary meshfree methods based on the boundary point interpolation methods[END_REF] (BPIM), Local Boundary Integral Equation Method [START_REF] Zhu | A local boundary integral equation (LBIE) method in computational mechanics, and a meshless discretization approach[END_REF] (LBIEM), etc. Although the boundary integral type meshless method benefits from the reduction of dimensions and the number of degrees of freedom, it is difficult to obtain the fundamental solutions of complex PDEs and to compute singular boundary integration. Mover, the requirement of integration leads to a low computational efficiency.

Due to the drawbacks of the first kind of boundary-type meshless methods, another boundary-type meshless methods based on fundamental solutions and collocation technique have been proposed. The computational efficiency is usually high since the discretization concerns only the boundary and no mesh or integration is required in this kind of algorithms. For instance, the Method of Fundamental Solution [START_REF] Kupradze | The method of functional equations for the approximate solution of certain boundary value problems[END_REF] (MFS) initialized by Kupradze and Aleksidze, Boundary Knot Method [START_REF] Chen | A meshless, integration-free, and boundary-only RBF technique[END_REF] (BKM) proposed by Chen et al., Regularized Meshless Method [START_REF] Young | Novel meshless method for solving the potential problems with arbitrary domain[END_REF] (RMM), Singular Boundary Method [START_REF] Chen | Singular boundary method: a novel, simple, meshfree, boundary collocation numerical method[END_REF] (SBM), etc. The main idea of the previously mentioned methods is to approximate the unknown function by a linear combination of the singular fundamental solutions and then the unknown coefficients are determined by applying the boundary conditions.

The present thesis focuses on a new boundary integration-free meshless method proposed by Zézé et al. [START_REF] Zézé | A boundary meshless method with shape functions computed from the PDE[END_REF], named Taylor Meshless Method (TMM), that relies on approximated solutions of the PDEs in the sense of Taylor series. When applied to a linear homogeneous PDE with constant coefficients, it coincides with Trefftz method associated with harmonic polynomials. In the following, few two dimensional applications of TMM are considered to briefly recall the basic properties of TMM, especially exponential convergence, robustness and efficiency.

Taylor Meshless Method

State of the art

The basic idea of TMM is to use the high order polynomial shape functions that are approximated solutions of the considered PDE. First we consider a simplest Laplace equation: ∆u(x, y) = 0 (1.1)

The harmonic polynomial solutions of the Laplace equation can be denoted by the real and imaginary parts of (x + iy) n , Re(x + iy) n and Im(x + iy) n . The general solution from order zero to p are collected in Table 1

P i • v i = P v (1.2)
where v i represents the unknown coefficients. Since Laplace equation has been fully satisfied by the approximation u p , only boundary conditions need to be considered to determine the unknown vector v. Here, we consider a mixed boundary conditions as follows:

u(x) = u d x ∈ Γ d , Tu(x) = t n x ∈ Γ n . (1.3)
Here, a collocation technique combined with the least-square method [START_REF] Zhang | Least-squares collocation meshless method[END_REF][START_REF] Fairweather | The method of fundamental solutions for elliptic boundary value problems[END_REF] is used to apply the boundary conditions. One chooses a set of nodes x i on Γ d and another set of nodes x j on Γ n , see Fig. Then one minimizes the error between the approximate value u p , Tu p and the given value of u d , t n at these points. It comes to minimize the following function:

T (v) = 1 2 x i ∈Γ d u p (x i ) -u d (x i ) 2 + w • 1 2
x j ∈Γn |Tu p (x j )t n (x j )| 2 (1.4)

The minimization leads to a linear system and solving this system gives the vector v and then the numerical solution of the considered problem.

The influence of the number of collocation points

Here we consider a Helmholtz equation in a rectangular domain:

    
-∆u + u = 0 in Ω u| y=0,4 = 0 u| x=±2.5 = sin(πy/4)

(1.5)
The exact solution is as following:

u(x, y) = cosh(x 1 + π 2 /16) cosh(2.5 1 + π 2 /16) sin(πy/4) (1.6)

The influence of the number of collocation points has been illustrated in Fig. 1.2. Three values of the degree are considered, p = 10, 15 and 20 respectively. The maximum accuracy is obtained for about M = 2p + 1 and it remains constant beyond this threshold. Almost the same behavior has been observed in the other 2D cases. In general, one chooses 4p collocation points to ensure the convergence.

Compared with FEM

Now we compare the accuracy of the proposed meshless method with classical finite element discretization. Here we consider Laplace equation in a rectangular domain (Ω = {(x, y)|0 x 10, 0 y π}), see Fig. In this discussion, we focus on the number of degrees of freedom necessary to get the same accuracy. The results are presented in Table 1.3. First one can see that the TMM converges exponentially with the degree of Taylor series while the FEM (with Q8 elements) converges slowly with the refinement of mesh. To obtain an accuracy of 10 -6 , only 31 DOFs are necessary with TMM while 9341 DOFs are required in FEM. The presented results show that a significant reduction in the number of DOFs has been obtained.

Convergence analysis of linear system

Another property of TMM we care about is the performance on convergence when solving a linear system. Here we consider Stokes equations in a square domain (see Fig. where f x and f y denote:

f x = (2µπ 2 -2π) sin(πx) sin(πy) f y = (2µπ 2 + 2π) cos(πx) cos(πy)
The boundary conditions are presented in Eq. (1.10).

       u = q = 0; v = -cos (πx); Γ 1 = {(x, y)|y = 1} u = 0; v = cos (πx); Γ 2 = {(x, y)|y = 0} u = 0; v = cos (πy); Γ 3 = {(x, y)|x = 0} u = 0; v = -cos (πy); Γ 4 = {(x, y)|x = 1} (1.10)
The exact solution of linear system Eq. (1.9) reads: One should note that there are three unknown functions in the considered case and the number of DOFs is 4p + 2. The convergence with the degree of Taylor series has been presented in Table 1.4. The maximum relative error decreases exponentially with the degree of Taylor series.

     u = sin (πx) sin (πy) v = cos (πx) cos (πy) q = 2 cos (πx) sin (πy) (1.11) * x c x i y x Γ 1 Γ 2 Γ 4 Γ 3 ℓ x ℓ y

Comments on TMM

In the previous papers [START_REF] Zézé | A boundary meshless method with shape functions computed from the PDE[END_REF][START_REF] Tampango | A new meshless method using Taylor series to solve elasticity problems[END_REF][START_REF] Tampango | Convergence analysis and detection of singularities within a boundary meshless method based on Taylor series[END_REF][START_REF] Tampango | Coupling of polynomial approximations with application to a boundary meshless method[END_REF], it has been validated that TMM is robust and efficient for solving two-dimensional elliptic PDEs. One can also use TMM to handle problems with any shape of domain and any boundary condition.

-10 - For instance, TMM works well for a problem with an amoeba-like boundary shape in [START_REF] Zézé | A boundary meshless method with shape functions computed from the PDE[END_REF] and for elasticity problems with displacement and stress boundary conditions in [START_REF] Tampango | Convergence analysis and detection of singularities within a boundary meshless method based on Taylor series[END_REF]. In practice, TMM needs much less DOFs as compared with other discretization techniques since the PDE is solved analytically.

The main content of this thesis

In this thesis, a boundary integration-free meshless method based on Taylor series, named Taylor Meshless Method (TMM), is discussed. The main idea of TMM is to compute a family of shape functions that are solutions of the PDE in the sense of Taylor series. And the specificity of TMM is the quasiexact solution inside the domain of the considered problems, which leads to a reduction of the dimensions and a strong reduction of the number of degrees of freedom. Some one-dimensional and two-dimensional applications of TMM have been proposed in the PhD thesis of Zézé [START_REF] Zézé | Calcul de fonctions de forme de haut degré par une technique de perturbation (Computing high degree shape functions by a perturbation technique[END_REF] and PhD thesis of Tampango [START_REF] Tampango | Développement d'une méthode sans maillage utilisant les approximations de Taylor (Development of a meshless method using Taylor approximations[END_REF], but also to two-dimensional non-linear elasticity in Master thesis of Akpama [START_REF] Akpama | Amélioration de la "Taylor meshless method" et application à des structures élastiques en grandes transformations[END_REF]. TMM is confirmed to be robust and efficient for at least two-dimensional cases with smooth solutions. The purpose of this thesis is to extend the applications of TMM to three dimensions, large-scale and non-linear range, but also try to initiate the problems with singularities. Along with continued research work, a variety of serious problems arise, for instance the performance of bridging technique to connect several sub-domains, the ill-conditioned matrix and the propagation of round-off errors for large-scale applications, the efficiency for non-linear applications and the handling of singular problems involving cracks, corners or notches. The detailed content of this thesis is as the follows:

• In chapter 2, in the framework of piecewise resolutions of TMM, Lagrange multipliers are applied to account for boundary conditions and least-square collocation is revisited to account for transmission conditions.

• In chapter 3, the applications of TMM are extended to three dimensions. The computation time and the influence of ill-conditioning of matrices on the accuracy are discussed carefully.

• In chapter 4, TMM is combined with Newton method and automatic differentiation to solve non-linear PDEs in one dimension and three dimensions.

• In chapter 5, few known singular solutions are introduced in TMM as singular shape functions to solve problems with singularities.

Introduction

Most of the numerical methods for PDEs rely on a priori chosen shape functions. Low degree polynomials are most wide-spread in the Finite Element Method. High-order polynomials permit to get very accurate solutions without re-meshing in the p-version of finite elements [START_REF] Babuska | The p-version of the finite element method[END_REF][START_REF] Schwab | The p and hp versions of the finite element method for problems with boundary layers[END_REF]. There were many works in the last two decades about radial functions [START_REF] Kansa | Multiquadrics-a scattered data approximation scheme with applications to computational fluid-dynamics-I surface approximations and partial derivative estimates[END_REF][START_REF] Kansa | Multiquadrics-a scattered data approximation scheme with applications to computational fluid-dynamics-II solutions to parabolic, hyperbolic and elliptic partial differential equations[END_REF], moving leastsquares approximations [START_REF] Belytschko | Element-free Galerkin methods[END_REF][START_REF] Nayroles | Generalizing the finite element method: diffuse approximation and diffuse elements[END_REF][START_REF] Nguyen | Meshless methods: a review and computer implementation aspects[END_REF] and Method of Fundamental Solution [START_REF] Fairweather | The method of fundamental solutions for elliptic boundary value problems[END_REF][START_REF] Golberg | The method of fundamental solutions for Poisson's equation[END_REF][START_REF] Golberg | The method of fundamental solutions for potential, helmholtz and diffusion problems[END_REF], among others. Nevertheless, except for small degree polynomials, it is not easy to manage either the integration procedures [START_REF] Frangi | A direct approach for boundary integral equations with high-order singularities[END_REF][START_REF] Frangi | 3d fracture analysis by the symmetric Galerkin BEM[END_REF] in case of a discretisation based on the weak form, or ill-conditioning matrices [START_REF] Kansa | Circumventing the ill-conditioning problem with multiquadric radial basis functions: applications to elliptic partial differential equations[END_REF] in case of collocation-based discretisation.

In this chapter, we are interested in the discretisation of PDEs by the method of Taylor series. Taylor series are used quite often for solving ordinary differential equations, see for instance [START_REF] Rentrop | A Taylor series method for the numerical solution of twopoint boundary value problems[END_REF][START_REF] Corliss | Solving ordinary differential equations using Taylor series[END_REF][START_REF] López | Two-point Taylor expansions of analytic functions[END_REF][START_REF] López | Multi-point Taylor approximations in one-dimensional linear boundary value problems[END_REF], but they have been introduced only recently for PDE's in a meshless and integration-less framework [START_REF] Zézé | A boundary meshless method with shape functions computed from the PDE[END_REF][START_REF] Tampango | A new meshless method using Taylor series to solve elasticity problems[END_REF][START_REF] Tampango | Convergence analysis and detection of singularities within a boundary meshless method based on Taylor series[END_REF][START_REF] Tampango | Coupling of polynomial approximations with application to a boundary meshless method[END_REF]. The solution is sought in the form of a polynomial with an arbitrary large degree and the size of the polynomial basis is reduced by considering the Taylor series of the PDE. So only the boundary of the domain has to be discretized. In this framework of large degree polynomials, the point-collocation is highly unstable, while the least-squares collocation [START_REF] Fairweather | The method of fundamental solutions for elliptic boundary value problems[END_REF][START_REF] Zhang | Least-squares collocation meshless method[END_REF] leads to a reliable numerical technique when applied together with Taylor series [START_REF] Zézé | A boundary meshless method with shape functions computed from the PDE[END_REF][START_REF] Zézé | Calcul de fonctions de forme de haut degré par une technique de perturbation (Computing high degree shape functions by a perturbation technique[END_REF]. Often the procedure converges exponentially with the degree (p-convergence), as in the p-version of the finite element method, which can lead to very accurate solutions. For instance, numerical tests in [START_REF] Tampango | Convergence analysis and detection of singularities within a boundary meshless method based on Taylor series[END_REF] have established that it requires much less degrees of freedom (DOFs) than finite element methods. Likely it could be applied to any linear or nonlinear elliptic system, as shown by the solution of hyperelastic boundary value problems in [START_REF] Akpama | Amélioration de la "Taylor meshless method" et application à des structures élastiques en grandes transformations[END_REF]. Nevertheless this convergence depends on the radius of convergence of the series and it is not possible to solve any boundary value problem with a single series.

Thus, numerical methods have to be introduced to connect several highorder polynomial approximations. A least-square approach was proposed in [START_REF] Zézé | A boundary meshless method with shape functions computed from the PDE[END_REF], but this piecewise resolution did not yield exponential convergence as with a single Taylor series. Such p-convergence properties have been established in [START_REF] Tampango | Coupling of polynomial approximations with application to a boundary meshless method[END_REF] by using two methods based on Lagrange multipliers. The first one is more or less similar with FETI [START_REF] Farhat | A method of finite element tearing and interconnecting and its parallel solution algorithm[END_REF][START_REF] Gosselet | Non-overlapping domain decomposition methods in structural mechanics[END_REF][START_REF] Stefanica | Parallel FETI algorithms for mortars[END_REF] or mortar methods [START_REF] Bernardi | A new nonconforming approach to domain decomposition: the mortar element method[END_REF], the continuity being enforced at some nodes of the interface. The second one defines the connection in a weak sense, which is a discrete version of Arlequin method [START_REF] Lions | Problèmes mécaniques multi-échelles: la méthode Arlequin[END_REF][START_REF] Ben Dhia | The Arlequin method as a flexible engineering design tool[END_REF]. In the two cases, the bridging techniques are not perfectly conforming and there are small discontinuities between the two polynomials, but of course these discontinuities decrease exponentially. With these recent results, an efficient and reliable method is established: it associates Taylor series inside each sub-domain, least-squares collocation for the boundary conditions and Lagrange multipliers at the interfaces between sub-domains. As in [START_REF] Tampango | Coupling of polynomial approximations with application to a boundary meshless method[END_REF], it will be called TMM. Thus we have two families of numerical techniques to account for boundary and interface conditions that have been validated in only one case. One may wonder if Lagrange multipliers may also be used for boundary conditions and least-squares at the interfaces. The aim of the chapter is to discuss these questions.

The chapter is organized as follows. In section 2.2, we present the boundary meshless method using Taylor series approximations. In section 5.3, two methods are presented to account for boundary conditions and some applications are considered to assess the validity of the two presented techniques: least-square collocation and Lagrange multiplier method based on radial functions. In section 2.4, least-square collocation is applied for the interface conditions in piecewise resolutions and the boundary conditions are ensured by the two presented methods respectively. In section 2.5, 2D elasticity problems are considered to discuss the behavior of TMM for mechanical problems. Finally in section 2.6, the robustness of the sub-domain techniques is assessed by considering a 3D benchmark that requires more than one million of unknowns within the finite element method.

State of the art

In this section, a meshless method based on Taylor series approximation is briefly recalled. More detailed descriptions of this method can be found in Zézé et al. [START_REF] Zézé | A boundary meshless method with shape functions computed from the PDE[END_REF][START_REF] Zézé | Calcul de fonctions de forme de haut degré par une technique de perturbation (Computing high degree shape functions by a perturbation technique[END_REF]. To illustrate the technique, we consider a Dirichlet problem:

-∆u + c • u = 0 in Ω u(x) = u d (x) on Γ (2.1)
where c is a constant. The main idea of this technique is to introduce high degree polynomial shape functions. This technique involves two steps. The first step consists in determining the shape functions by a quasi-exact resolution of the PDE in the domain. The approximate solution of Eq. (2.1) is sought in the form of a polynomial of degree p:

u p (x, y) = p k=0 k i=0 u i,k-i x i y k-i = p k=0 X k { u k } (2.2)
where

X k = x k , x k-1 y, • • • , y k , { u k } t = u k,0 , u k-1,1 , • • • , u 0,k .
For each degree k, the unknown is the vector { u k } ∈ R k+1 . For the complete polynomial Eq. (2.2), there are (p + 1)(p + 2)/2 coefficients to be found. Then if Eq. (2.1) is satisfied, the polynomial -∆u + c • u vanishes and this leads to a linear system of equations:

-[∆ k ]{ u k+2 } + c{ u k } = {0} ∀ k, 0 ≤ k ≤ p -2 (2.3)
For each k, there are k+1 equations in Eq. ( 2.3) for k+3 unknowns ({ u k+2 } ∈ R k+3 ). Each vector { u k+2 } can then be written as a function of its two first components:

v k = u k,0 u k-1,1 ∈ R 2 (2.4)
The account of the PDE reduces the general polynomial Eq. (2.2) to a family of (2p + 1) polynomials instead of (p + 1)(p + 2)/2:

u p (x, y) = 2p+1 i=0 P i (x, y)v i (2.5)
For instance in the case of Laplace equation, one recovers the well known result: P 2n = Re(x+iy) n and P 2n+1 = Im(x+iy) n . A similar procedure can be applied to nonlinear equations or to equations with variable coefficients [START_REF] Zézé | Calcul de fonctions de forme de haut degré par une technique de perturbation (Computing high degree shape functions by a perturbation technique[END_REF][START_REF] Akpama | Amélioration de la "Taylor meshless method" et application à des structures élastiques en grandes transformations[END_REF], but we will not give more details in this paper that focuses on the treatment of boundary and transmission conditions. At this level, the coefficients {v} are unknown. So the approximate solution of the problem is completely obtained by determining these 2p + 1 variables. That is the goal of the second step. It consists in the application of boundary conditions. A least-square method combined with collocation technique can be used. This technique is rather common within meshless methods and avoids numerical instabilities occurring with pure collocation [START_REF] Fairweather | The method of fundamental solutions for elliptic boundary value problems[END_REF][START_REF] Zhang | Least-squares collocation meshless method[END_REF]. One chooses a set of nodes x j on the boundary of the domain and one minimizes the error between the approximate values and the given values of u in these points. It comes to minimize the function:

T (v) = 1 2 M j=1 u p (x j ) -u d (x j ) 2 (2.6)
This minimization leads to a linear system [K]{v} = {b} where [K] is an invertible matrix. Solving this system gives the vector {v} and so one gets the numerical solution of the problem Eq. (2.1).

Methods for setting boundary conditions

Two methods to account for boundary conditions

Previous papers about Taylor Meshless Method apply least-square collocation to account for the boundary conditions. Here a new approach based on discrete Lagrange multipliers is presented that looks like the techniques used in [START_REF] Tampango | Coupling of polynomial approximations with application to a boundary meshless method[END_REF] for the transmission conditions.

Least-square collocation

The first method to account for the boundary conditions is the least-square collocation as shown in previous section. After solving the PDE inside the domain, see Eq. (2.5), one chooses M collocation points x j on the boundary of the domain and minimizes the error between the approximate values and the given values of u in these points, according to Eq. (2.34). Thus this first procedure is completely defined, first by the point where the Taylor series are developed (x = [0, 0] in section 2.2) that will be called development point, second by those boundary collocation points. Previous results have established that the procedure is quite robust and converges provided that the number of collocation points is large enough. This question of robustness will be re-discussed in this paper. Other boundary conditions can be studied by functions similar as Eq. (2.34); an example will be discussed in section 2.5 within linear elasticity.

Lagrange multipliers

The second method to account for the boundary conditions is based on Lagrange multipliers and it has many common characteristics with the method introduced in [START_REF] Tampango | Coupling of polynomial approximations with application to a boundary meshless method[END_REF] for the discretization of the transmission conditions. One starts from the boundary conditions in a discrete weak form and the Lagrange multiplier is discretized by radial functions that are wide-spread within meshless methods [START_REF] Kansa | Multiquadrics-a scattered data approximation scheme with applications to computational fluid-dynamics-I surface approximations and partial derivative estimates[END_REF][START_REF] Kansa | Multiquadrics-a scattered data approximation scheme with applications to computational fluid-dynamics-II solutions to parabolic, hyperbolic and elliptic partial differential equations[END_REF]. For simplicity, let us consider only the Dirichlet boundary condition:

u(x) = u d (x) on Γ (2.7)
The weak form of Eq. (2.7) involves a Lagrange multiplier µ(x) defined on the boundary:

Γ µ(x)[u(x) -u d (x)]dΓ = 0 ∀ µ(x) defined on Γ (2.8)
Two sets of discretization points have to be chosen on the boundary: a cloud of M collocation points (x k ∈ Γ) and a cloud of N points (x i ∈ Γ) to define the discrete form of the Lagrange multiplier µ(x). Then the continuous weak form Eq. (2.8) is replaced by a discrete one:

M k=1 µ(x k )[u(x k ) -u d (x k )] = 0 ∀ µ(x) defined on Γ.
(2.9)

As in [START_REF] Tampango | Coupling of polynomial approximations with application to a boundary meshless method[END_REF], we introduce a family of Gaussians, associated to the centers x i and depending on a characteristic distance d:

Φ i (x k ) = exp( -x k -x i 2 d 2 ) (2.10)
These Gaussians are chosen as weighting functions in the variational equation (2.9):

µ(x k ) = N i=1 µ i Φ i (x k ) or µ(x k ) = Φ i (x k ) (2.11)
The number of discrete Lagrange multipliers µ i must be the same as the number of the remaining degrees of freedom in Eq. (2.5) in order to obtain a square matrix N = 2p + 1. Hence the proposed method can be seen as a boundary-collocation method, where the collocation conditions are satisfied in a mean sense. By combining Eqs. (2.5), (2.9) and (2.11), we arrive at a linear problem in R N :

[K]{v} = {b} (2.12)

where

K ij = M k=1 Φ i (x k )P j (x k ) and b i = M k=1 Φ i (x k )u d (x k ) (2.13)

Application in a rectangular domain

In this part, an application will be discussed to assess the validity of the two presented techniques accounting for boundary condition. The first goal is to show the accuracy of the solutions. The second goal is to discuss the influence of each parameter, check the robustness of the technique and get the optimal values of the parameters. Here we consider the Dirichlet problem in a rectangular domain, 0 x 10, 0 y π.

     -∆u = 0 in Ω u(x, 0) = u(x, π) = u(10, y) = 0 u(0, y) = sin(y) (2.14)
The exact solution of Eq. (2.14) is

u(x, y) = sin(y) • sinh(10 -x) sinh(10) (2.15)
This problem is solved by TMM with a development point at c = [0, 0] in the following. Throughout this paper, the relative error is the difference between exact and approximate solutions, divided by the maximum value of the solution:

E = max|u(x) -u h (x)| max|u(x)| .
(2.16)

TMM with least-square collocation

The Taylor meshless method with least-square collocation has been previously discussed in past papers [START_REF] Zézé | A boundary meshless method with shape functions computed from the PDE[END_REF][START_REF] Tampango | Convergence analysis and detection of singularities within a boundary meshless method based on Taylor series[END_REF][START_REF] Tampango | Coupling of polynomial approximations with application to a boundary meshless method[END_REF]. The numerical solution obtained by TMM with least-square collocation depends on two parameters: the degree p of the polynomials and the number of collocation points M . It has been established in [START_REF] Zézé | A boundary meshless method with shape functions computed from the PDE[END_REF][START_REF] Tampango | Coupling of polynomial approximations with application to a boundary meshless method[END_REF][START_REF] Zézé | Calcul de fonctions de forme de haut degré par une technique de perturbation (Computing high degree shape functions by a perturbation technique[END_REF] that the minimal number of collocation points is a little higher than 2p. This will be re-discussed here, but generally we shall choose M = 4p as recommended in [START_REF] Zézé | A boundary meshless method with shape functions computed from the PDE[END_REF][START_REF] Tampango | Convergence analysis and detection of singularities within a boundary meshless method based on Taylor series[END_REF]. The p-convergence illustrated in Fig. 2.1 shows that one can get accurate results by increasing the degree. The error decreases more or less exponentially with the degree (p-convergence). The accuracy is improved until p = 25 where it becomes stationary. Next, due to the numerical noise related to the accuracy limit of the computer, a small infection is observed, but the accuracy stands around 10 -8 .

TMM with Lagrange multipliers

The numerical solution obtained by TMM with Lagrange multipliers mainly depends on three parameters: the degree p of the polynomials, the number of collocation points M and the effective radius d of the radial functions. The p-convergence is illustrated in Fig. 2.1. The same convergence is obtained as in TMM with least-square collocation with a very good accuracy.

The influence of the number of collocation points in TMM with Lagrange multipliers has been illustrated in Fig. 2.2 for three values of degree p = 5, p = 10 and p = 15. The maximal error decreases with M until an optimal number where it becomes stable. In this case, the optimal number is a little higher than 3p. It can be seen in Fig. 2.2 that 4p is large enough to ensure the best convergence. The parameter d is the radius of the influence zone around each Lagrange multiplier node. Four collocation points will be located in the influence zone of each Lagrange multiplier node when d = 2 * δd is chosen, where the δd represents the minimum distance between two adjacent collocation points. The results are presented in Fig. 2.3. These results are obtained with M = 4p collocation points and different degrees from p = 8 to p = 44. The results show that the influence of parameter d is increasing with the degree of the polynomial. One can get always an optimal convergence when the chosen parameter d is less than 2 * δd.

Application to Laplace equation in an unit circle

Here, we consider problem Eq. (2.1) in a unit disk x 2 + y 2 1 with c = 0:

-∆u(x) = 0 x ∈ Ω u(x) = u d (x) x ∈ Γ (2.17)
where u d (x, y) = (xx 0 )/[(xx 0 ) 2 + (yy 0 ) 2 ] and the singular point is set to be x 0 = [START_REF] Zienkiewicz | The finite element method (Fifth edition) Volume 1: The Basis[END_REF][START_REF] Zienkiewicz | The finite element method (Fifth edition) Volume 1: The Basis[END_REF]. The exact solution of Eq. (2.17) reads:

u(x) = (x -x 0 ) (x -x 0 ) 2 + (y -y 0 ) 2 (2.18)
The same case solved by TMM with least-square collocation has been presented in Tampango et al. [START_REF] Tampango | Coupling of polynomial approximations with application to a boundary meshless method[END_REF]. The influence of the degree of the polynomials shape functions p and the number of collocation points M have been well discussed. Here, we will only analyze the influence of the parameters on the algorithm with Lagrange multipliers. The p-convergence is illustrated in Tab. 2.1. It can be seen that one converges exponentially by increasing the degree.

The influence of the number of collocation points is illustrated in Fig. 2.4 for three values of degree p = 10, 20 and 30. First one checks that the method fails if the number of collocation points is too small. Next the maximal error decreases with M until an optimal number where it becomes stable. The optimal number is a little higher than 2p. It can be seen in Fig. 2.4 that 4p is large enough to ensure the best convergence. Here, one gets the same conclusion as in section 2.3.2.

The influence of the parameter d is illustrated in Fig. 2.5. The convergence keeps stable with a good accuracy by increasing the value of the radius d of the radial functions until a value d = 6 * δd where the maximal error increases rapidly. The influence of d becomes sensitive with high degree polynomials. Anyway, TMM with Lagrange multipliers always works well with small values of parameter d.

Comments of applying boundary conditions

A new technique based on Lagrange multiplier and radial functions was presented to account for boundary conditions and assessed by several numerical tests. Contrarily to the pure collocation method, it converges, is robust and the error decreases exponentially with the degree in all the considered examples. It has about the same efficiency as the least-square collocation previously used within TMM. It depends on three parameters: the degree p, the number of collocation points M and d the characteristic distance appearing in the gaussian functions. A safe carrying out is obtained for M = 4p and d/δd between 1 and 4, δd being the distance between neighbor collocation points. With these choices, the algorithm depends only on the degree p that controls the convergence. As in previous versions of TMM, very high degrees (say 20, 40) can be chosen, which can yield very accurate solutions without huge computational cost unlike high degree finite elements. Nevertheless the previously used least-square collocation method is as efficient as the new one and has the advantage to need only one family of discretization points.

Using least-square collocation to connect various sub-domains

It has not been proved yet that the least-square collocation method can be applied to the transmission conditions by keeping the property of exponential convergence with the degree. By comparison with Lagrange multipliers, it has the drawback to lead to full matrices, but the number of degrees of freedom is smaller and it is more easily defined. Only a set of collocation points is needed and their number and position can be chosen very freely. Of course this can be associated with a treatment of the boundary conditions and the two methods discussed in section 5.3 will be investigated.

Full least-square method

To complete the solution of the PDE by Taylor series obtained in section 2.2, a numerical procedure requires a treatment of boundary and transmission conditions. Let us begin by the simplest procedure where least-square collocation is applied in the two cases.

Statement

Let us consider an elliptic equation in a domain Ω

Lu = f in Ω u(x) = u d (x) on Γ (2.19)
where L is a second order differential operator.

In order to make a piecewise resolution, an interface (Γ r ) is introduced by splitting the domain Ω into two sub-domains Ω 1 and Ω 2 (see Fig. The discretization of the domain is generated by two types of nodes:

• nodes on the external boundaries Γ 1 and Γ 2 for boundary conditions.

• nodes at the interface Γ r for transmission conditions. The two conditions will be accounted by using a least-square method as shown in the previous section. Thus a set of nodes is chosen on each sub-boundary and on the interface: • a set of M 1 collocation nodes on the boundary Γ 1 .

• a set of M 2 collocation nodes on the boundary Γ 2 .

• a set of M r collocation points on the interface Γ r .

With the initial formulation of the TMM, it is possible to build polynomial shape functions for problems (2.19) in each sub-domain:

u 1 (x) = P 1 0 (x) + 2p 1 +1 i=1 P 1 i (x)v 1 i (2.20)
-25 -Chapter 2. Piecewise resolution of Taylor Meshless Method

u 2 (x) = P 2 0 (x) + 2p 2 +1 i=1 P 2 i (x)v 2 i (2.21)
where the P 1 0 (x) and P 2 0 (x) are used to balance the right-hand side analytical function f (x) of Eq. (2.19).

The function accounting for boundary conditions on Γ 1 and Γ 2 is exactly the same as in previous papers:

T 1 (v 1 ) = 1 2 M 1 j 1 =1 P 1 0 (x j1 ) + 2p 1 +1 i=1 P 1 i (x j 1 )v 1 i -u d (x j 1 ) 2 (2.22) T 2 (v 2 ) = 1 2 M 2 j 2 =1 P 2 0 (x j2 ) + 2p 2 +1 i=1 P 2 i (x j 2 )v 2 i -u d (x j 2 ) 2 (2.23)
In the same way, the two transmission conditions will be ensured via the following coupling function:

C (v 1 , v 2 ) = 1 2 Mr jr=1 u 1 (x jr ) -u 2 (x jr ) 2 + ∂u 1 ∂n (x jr ) - ∂u 2 ∂n (x jr ) 2 (2.24)
Then for the whole problem, searching the variables {v (1) } and {v (2) } is equivalent to minimizing the following function:

T (v 1 , v 2 ) = T 1 (v 1 ) + T 2 (v 2 ) + C (v 1 , v 2 ) (2.25)
The minimization of the function T (v 1 , v 2 ) leads to a linear symmetric and invertible system as follows:

K 11 K 12 K 21 K 22 v 1 v 2 = b 1 b 2 (2.26)

Application to Laplace equation in a disk

In this section, the numerical solutions obtained by this piecewise approach are presented to study the validity of the TMM proposed in previous section. The Laplace problem Eq. (2.17 The parameters of the discretization are the following ones:

• approximation degrees in the sub-domains (p 1 and p 2 ).

• number of collocation on the boundary of each sub-domain (M 1 and M 2 ).

• center of the Taylor expansion in each sub-domain (x c 1 and x c 2 ).

• number of collocation nodes on the interface (M r ).

The development points of the approximate series are chosen at x c 1 = [-0.5, 0] for the sub-domain on the left-hand side and at x c 2 = [0, 0] on the right-hand side. The influence of the other parameters will be discussed in details.

* * x c 1 x c 2 Ω 1 Ω 2 Figure 2
.7: Two sub-domains for the Laplace problem (2.17).

We present the error on the inner circle of radius r = 0.8 in Fig. 2.8. This circle encounters the interface at θ = π/2 and θ = 3π/2. The Fig. 2.8 presents the approximate solution and the exact solution for p 1 = p 2 = p = 10,

M 1 = M 2 = 2p, M r = p.
The figure shows a good correspondence between the exact and the approximate solution. Thus the transmission conditions have been well accounted for. The optimal number of boundary collocation points was previously discussed [START_REF] Zézé | A boundary meshless method with shape functions computed from the PDE[END_REF]: a choice M 1 = 2p, M 2 = 2p is sufficient to ensure the convergence. Beyond this limit, the error is stable. The method works well with more points, but this does not improve the convergence.

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. The Tab. 2.2 presents the error of the whole domain for several approximation degrees and for different number of collocation points M r . As for the boundary collocation, the error decreases when M r increases and becomes constant from some threshold, around M r = p. From this value, the convergence is stabilized for larger M r . Furthermore, the convergence is rapidly improved by increasing the degree. The above results prove the robustness of the least-square technique applied for both boundary and interface condition. In the previous case, the polynomials on both sides of the interface had the same degree. Now we are interested in the case with different degrees in each sub-domain, always in order to test the robustness of the coupling technique. The results, presented in the Tab. 2.3, are similar to the previous approximation with a common degree. The error is stabilized from some threshold around M r = (p 1 + p 2 )/2. One sees that a coarse discretization is sufficient in the left part, in which the gradient is smaller: indeed one gets about the same accuracy with (p 1 , p 2 ) = (7, 15) as with p 1 = p 2 = 15. This also establishes the robustness of the least-square coupling technique that works well with different degrees.

     -u(x, y) = - 4 (x 2 + y 2 ) 2 in Ω u d (x, y) = 1 x 2 + y 2 on Γ (2.27)
The domain is a crown (Ω = {(x, y)/r 2 1 ≤ x 2 + y 2 ≤ r 2 2 }, with r 1 = 0.8 and r 2 = 1). The exact solution of Eq. (2.27) is u(x, y) = 1/x 2 + y 2 . This problem has been solved in [START_REF] Tampango | Coupling of polynomial approximations with application to a boundary meshless method[END_REF]. The result obtained from TMM with one single domain were always completely wrong. That is due to the singularity at x 0 = [0, 0] in the middle of crown. There is no development point that can cover the whole domain by avoiding the singularity (see Fig. 2.9).

*

x 0

x c

Domain of convergence

Figure 2.9: Discretization of the crown for a piecewise resolution.

In the paper [START_REF] Tampango | Coupling of polynomial approximations with application to a boundary meshless method[END_REF], this problem has been solved by dividing the whole domain into several sub-domains coupled by using Lagrange multipliers and Arlequin method. Here we are interested in the behavior of the least-square collocation approach used both for boundary and interface conditions.

To validate TMM coupled with the least-square technique, we split the crown into several sectors (Fig. 2.10). In each sector, the development point of the approximate series (x c i ) has been chosen in the middle of the outer arc. The Fig. 2.11 presents the maximal error of the TMM model with different numbers of sub-domains. We choose p collocation points on each interface and 2p collocation points on the boundary. This result shows that the convergence is improved by increasing the number of sub-domains. The main reason is that each Taylor series can widely cover its sub-domain by avoiding the singularity.

* x 0 Ω 2 Ω 1 Ω 3 Ω 4 Ω 5 Ω 6 x c 1 x c 2 x c 3 x c 4 x c 5 x c 6
The convergence with the degree is illustrated in Fig. 2.12. Here the crown is split into six sub-domains. The results show that the convergence is improved by increasing the degree until a value p = 20 where it becomes stable, with a very good accuracy. Likely such a plateau is due to numerical noise related to the accuracy limit of the computer. The obtained accuracy (10 -4.19 for p = 10, 10 -7.49 for p = 20) is almost the same as those of Reference [START_REF] Tampango | Coupling of polynomial approximations with application to a boundary meshless method[END_REF] where various Lagrange multipliers techniques were applied for the interface conditions.

Outcome of full least-square method

The full least-square method works well. It is robust and does not depend on the number of collocation points. It converges exponentially with the degree, about as the previous technique least-square/Lagrange multipliers [START_REF] Tampango | Coupling of polynomial approximations with application to a boundary meshless method[END_REF]. Thus it provides an interesting alternative.

A mixed Lagrange/least-square method

We now check that transmission least-square collocation is also compatible with a Lagrange multiplier discretization of boundary conditions. This will be complementary to the method in [START_REF] Tampango | Coupling of polynomial approximations with application to a boundary meshless method[END_REF] where the interface condition was treated by Lagrange multipliers and the boundary condition by least-square collocation.

Statement

In this section, Lagrange multiplier method is applied for boundary condition in a piecewise resolution and least-square collocation for interface condition.

Let us consider a Poisson problem with Dirichlet boundary condition:

-u(x) = q(x) in Ω u(x) = g(x)
on Γ (2.28)

An interface (Γ r ) is introduced by splitting the domain Ω into two subdomains Ω 1 and Ω 2 (see Fig. Here, transmission conditions will be considered as the target function and discretized by the least-square collocation. For this purpose a set of M r collocation points is chosen on the interface and transmission conditions would be satisfied in a mean sense via Eq. (2.24). The weak form of the boundary condition involving Lagrange multipliers µ s (x)(s ∈ [START_REF] Zienkiewicz | The finite element method (Fifth edition) Volume 1: The Basis[END_REF]2]) reads:

Γs µ s (x)[u s (x) -g s (x)]dΓ = 0 ∀ µ s (x) defined on Γ s .
(2.29)

Then for the whole problem, searching the variables v and µ is equivalent to vanishing the gradient of the following function:

L (v, µ) = C (v 1 , v 2 ) + 2 s=1 Γs µ s (x)[u s (x) -g s (x)]dΓ (2.30)
where C (v 1 , v 2 ) is the same as in Eq. (2.24).

As in Section 5.3, two sets of discretization points have to be chosen on the boundary in each sub-domain:

• M s collocation points (x m ∈ Γ s ) to discrete the boundary condition.

• N s points (x n ∈ Γ s ) to define the Lagrange multiplier.

Then the discretization form of Eq. (2.30) reads:

L (v, µ) = 1 2 t {v}[K r ]{v} -t {b r }{v} + t {µ}[C]{v} -t {µ}{b} (2.31)
where

[K r ] = K r11 -K r12 -K r21 K r22 , [C] = C 1 0 0 C 2 [K rst ] ij = Mr m=1 P s i (x m )P t j (x m ) + ∂P s i (x m ) ∂n ∂P t j (x m ) ∂n . [C s ] ij = Ms n=1 Φ s i (x n )P s j (x n ). [b] = b 1 b 2 , [b r ] = b 1 r b 2 r . [b s ] i = Ms n=1 Φ s i (x n ) (g s (x n ) -P s 0 (x n )) . [b s r ] i = Mr m=1 P s i (x m ) P 2 0 (x m ) -P 1 0 (x m ) + ∂P s i (x m ) ∂n ∂P 2 0 (x m ) ∂n - ∂P 1 0 (x m ) ∂n .
(Note that: s, t = 1, 2 and x m ∈ Γ r and x n ∈ Γ s .)

The stationarity of the function L (v, µ) leads to a linear symmetric and invertible system as follows:

K r t C C 0 v µ = b r b (2.32)

Application to Laplace equation in a disk

Let us consider the Laplace problem (2.17) with a singularity at [START_REF] Zienkiewicz | The finite element method (Fifth edition) Volume 1: The Basis[END_REF][START_REF] Zienkiewicz | The finite element method (Fifth edition) Volume 1: The Basis[END_REF]. We solve this problem in a piecewise way by the method described in section 2.4.2.1. The domain is split into two sub-domains. The development points of the Taylor series are chosen at x c 1 = [-0.5, 0] for the sub-domain on the left-hand side and at x c 2 = [0, 0] on the right-hand side. Sufficiently large numbers of collocation points and of Lagrange multipliers have been chosen in such a way that the accuracy is optimal. The radial functions needed to discretize the Lagrange multipliers depend on a length d that is compared with the minimal distance between neighbor collocation points, called δd. In Tab. 2.4, the maximal error is reported as a function of the degree and of the ratio d/δd that was varied in the range [START_REF] Zienkiewicz | The finite element method (Fifth edition) Volume 1: The Basis[END_REF][START_REF] Lucy | A numerical approach to the testing of the fission hypothesis[END_REF]. Clearly the efficiency of the method depends very little on this ratio that can be chosen rather freely.

Table 2.4: The p-convergence for Laplace Eq. (2.17) solved by using two subdomains. δd represents the minimum distance between to adjacent collocation points. These results can be compared with the single domain solution presented in section 2.3.3, to the pure least-square method of section 2.4.1 or the method of reference [START_REF] Tampango | Coupling of polynomial approximations with application to a boundary meshless method[END_REF] that combines Lagrange multipliers and least-square collocation, but in the opposite way, see Tab. 2.5. Clearly all these methods give almost the same results and all converge exponentially with the degree. This means that the discretization of boundary and transmission conditions can be done by either method, without significant influence on the effectiveness.

M 1 = M 2 = 2p, N 1 = N 2 = p and M r = 2p.

Application to Poisson equation in a crown

Poisson Eq. (2.27) is reconsidered in this section and solved by the method of section 2.4.2.1. The crown has been split in six sub-domains with the same development points as in section 2.4.1.3. The degree varies and the numbers of discretization nodes are as follows: M r = p, M s = 2p and N s = 2p. As for the previous example, the influence of the ratio d/δd is studied in the range [START_REF] Zienkiewicz | The finite element method (Fifth edition) Volume 1: The Basis[END_REF][START_REF] Lucy | A numerical approach to the testing of the fission hypothesis[END_REF]. Roughly, this influence is very weak except sometimes for d = 6 * δd, where the accuracy is slightly lower. The accuracy obtained by this Lagrange/leastsquare technique has been compared with the same alternative methods as for the previous example. The results illustrated in Tab. 2.7 confirm that all these methods are more or less equivalent. Let us mention the quasi-perfect solution obtained by the present method with a degree p = 30 and d = δd or d = 2 * δd: the error is less than 10 -12 .

Table 2.6: The maximal relative error of the solutions for Poisson problem (2.27) split into six sub-domains. δd represents the minimum distance between two adjacent collocation points. M r = p, M s = 2p and N s = 2p, where M r , M s and N s represent the number of collocation points on the interface, the number of collocation points on the boundary of each sub-domain and the number of discretization points of Lagrange multipliers µ s on the boundary of each sub-domain.

Comments about applying interface conditions

In this section, the discretization of the transmission condition by least-square collocation has been discussed and compared with previous techniques based on Lagrange multipliers. From numerical tests, we can conclude that all these methods have about the same efficiency and robustness. All converge exponentially with the degree and are few sensitive to the discretization parameters.

Application in 2D elasticity

The previous examples concern a single elliptic equation of Laplace or Poisson type. We now discuss examples in 2D elasticity. In the literature, one can find few applications of TMM in linear elasticity [START_REF] Tampango | A new meshless method using Taylor series to solve elasticity problems[END_REF][START_REF] Tampango | Convergence analysis and detection of singularities within a boundary meshless method based on Taylor series[END_REF] or even in nonlinear elasticity [START_REF] Akpama | Amélioration de la "Taylor meshless method" et application à des structures élastiques en grandes transformations[END_REF]. It seems that the algorithm works about in the same way for a system as for a single PDE, but this could be no longer the case if the sought solution is not a perfectly smooth function. In real life, one needs to represent point-wise forces, corners or cracks, in which cases the solution is more or less singular. Thus the representation of non-smooth solutions by Taylor series is a serious challenge and this research will be initiated here by an example with a square domain.

The reliability is a central point in the assessment of a new numerical method. In the previous examples presented in this chapter, the convergence is exponential with a limit at a very high accuracy of 10 -7 or 10 -8 . This limit may be greater, for instance if the development point is not well chosen, see for instance Fig. 4 in [START_REF] Tampango | A new meshless method using Taylor series to solve elasticity problems[END_REF]. Thus for the sake of reliability, it is important that this limit of accuracy remain sufficiently low.

In this part we apply least-square collocation whose reliability and efficiency were proved in the front part of this paper. Here, we consider 2D elasticity problems in the generic form:

     Gu(x) = 0 in Ω u(x) = u d (x) on Γ u T u(x) = t d (x) on Γ t (2.33)
where G and T are differential operators:

G = E 1 -ν 2     ∂ 2 ∂x 2 + 1 -ν 2 ∂ 2 ∂y 2 1 + ν 2 ∂ 2 ∂x∂y 1 + ν 2 ∂ 2 ∂x∂y ∂ 2 ∂y 2 + 1 -ν 2 ∂ 2 ∂x 2     T = E 1 -ν 2    n x ∂ ∂x + n y 1 -ν 2 ∂ ∂y n x ν ∂ ∂y + n y 1 -ν 2 ∂ ∂x n y ν ∂ ∂x + n x 1 -ν 2 ∂ ∂y n y ∂ ∂y + n x 1 -ν 2 ∂ ∂x   
in which ν is the Poisson ratio and E is the Young modulus. The n x and n y are the components of the normal to the boundary. Assuming u h (x) has satisfied the control equation Gu(x) = 0, the only thing need to be considered is the boundary conditions. It comes to minimize the following function:

J(v) = 1 2 x j ∈Γu u h (x j ) -u d (x j ) 2 + w • 1 2 x j ∈Γt T u h (x j ) -t d (x j ) 2
(2.34)

A 2D elasticity problem without singularity

In this section, we consider an elasticity problem (2.33) with a uniformly distributed force σ x = q 1 along x = L and σ y = q 2 along y = 0 and y = B (see Fig. The constants are as follows: E = 1000, ν = 0.3, L = 10, B = 10, q 1 = 100 and q 2 = ν • q 1 . The development point is chosen at the center X 0 = [5,5]. The analytical solution in this case varies linearly:

u(x, y) = x E (q 1 -νq 2 ) v(x, y) = 0 (2.35)
where the u and v represent the displacements along x-axis direction and yaxis direction respectively. The numerical results presented in Tab. 2.8 show that TMM satisfies naturally the consistency condition. It is able to recover a spatially linear function with a very high accuracy. Nevertheless this accuracy slightly decreases for large degrees due to numerical noise, but it remains lower than 10 -10 .

A 2D elasticity problem with singularity

In this section, we consider a 2D elasticity problem (2.33) with a linear distribution of force σ x along the end (see Fig. 2.15). In this case, there will be a singularity in the top-left corner of the square. The data are as follows: E = 1000, ν = 0.3, L = 10, B = 10 and q = 100. This example has no analytical solution so that the reference solution has to be defined numerically. Several finite element solutions have been computed with linear and quadratic elements and with several meshes. We consider that a 500 × 500 mesh and four node quadrilateral elements (Q4) leads to a sort of reference solution.

Indeed with respect to a 250 × 250 mesh of Q4 elements or a 100 × 100 mesh of Q8 elements, the difference is less than 6.8 × 10 -4 .

TMM with one domain

A single Taylor series is first considered with a development point at the center of the square X 0 = [5,5]. The numerical results presented in Tab.

2.9 indicate that the convergence is improved slowly by increasing the degree of the polynomials. There is an accuracy limit that is never less than 1%.

According to the distribution of the error presented in Fig. 2.16, this maximal error is mainly governed by the singularity located at the left corner of the square [0, 10]. The relative errors at the point X = [START_REF] Lancaster | Surfaces generated by moving least squares methods[END_REF][START_REF] Lancaster | Surfaces generated by moving least squares methods[END_REF] that is far away from the singularity point is illustrated in Fig. 2.17. One observes also a Table 2.9: TMM with a single domain for the 2D elasticity problem shown in Fig. 2.15. The maximal error in the whole domain.

Relative error p = 10 p = 20 p = 30 p = 40

Max error of u 0.009 0.009 0.0055 0.0035 Max error of v 0.080 0.030 0.018 0.012

TMM with piecewise resolution

In this section, we split the whole domain into two or three subdomains, as pictured in Fig. 2.18. Least-square method is applied both for boundary and interface conditions. The comparison with the solution obtained by a single domain resolution (Tab. 2.9) shows that the piecewise resolution improves only slightly the convergence, the order of magnitude of the error remaining about 0.8%. Several other tests have been performed, but this did not lead to significant improvements. Hence it seems difficult to represent very accurately more or less singular solutions by Taylor series, even with piecewise resolutions. Other tests involving cracks will confirm this result, see chapter 5. Thus one of the next challenges for TMM is to find new procedures for a better account of singular solutions. The results presented in Tab. 2.10 and Tab. 2.11 are obtained by using TMM with sub-domains. Least-square method is applied for boundary and

Degree p

Log10(relative error)

The convergence at the point [START_REF] Lancaster | Surfaces generated by moving least squares methods[END_REF][START_REF] Lancaster | Surfaces generated by moving least squares methods[END_REF] error of u error of v interface conditions. The comparison with the solution obtained by a single domain resolution (Tab. 2.9) shows that the piecewise resolution slightly improves the convergence.

6,5] [2,8] [2,4] [2,5] [2,8] Ω 1 Ω 2 Ω 1 Ω 2 Ω 3

A very large scale test

The previous examples concern only two dimensional partial differential equations and have been solved by TMM with a single domain or a few number of sub-domains. It was established that the Lagrange multiplier and least-square collocation are both robust and efficient when accounting for boundary and interface conditions. However the number of degrees of freedom of the previous cases are really small. Hence one considers a three dimensional case divided into many sub-domains to check the robustness and efficiency of TMM when solving large scale problems. This benchmark is motivated by wrinkling instabilities of film-substrate systems [START_REF] Chen | Herringbone buckling patterns of compressed thin films on compliant substrates[END_REF][START_REF] Xu | 3D finite element modeling for instabilities in thin films on soft substrates[END_REF], where one observes many spatial oscillations. Here one accounts only for the substrate and the film is represented by a prescribed sinusoidal displacement. The domain is a rectangular hexahedron with sides l x × l y × d, see Fig. 2.19. One considers the Laplace equation with Dirichlet boundary conditions:

     u(x, y, z) = 0 in Ω u(x, y, 0) = 0 u(x, y, d) = sin(πx)sin(πy) (2.36)
The exact solution is u(x, y, z) = sin(πx)sin(πy)sinh( √ 2πz)/sinh( √ 2πd). The domain is split in cubes of size 1 × 1 × 1 that will be called "elements" for simplicity, see Fig. 2.20. This corresponds to one element through the thickness and two elements along each period. The boundary conditions are accounted by using the least-square method and three bridging techniques will be compared to account for the interface conditions: full least-square described in section 2.4.1, Lagrange multipliers and Arlequin methods as introduced in [START_REF] Tampango | Coupling of polynomial approximations with application to a boundary meshless method[END_REF]. First one considers the domain l x × l y × d = 10 × 10 × 1 that is discretized by 100 TMM-elements. The convergence of the algorithm with the degree (p-convergence) is presented in Fig. 2.21. The three bridging techniques lead about to the same results. The error decreases exponentially with the degree up to a very high accuracy, even for this large scale problem. This establishes the robustness of all these procedures with respect to the number of sub-domains. Next these TMM-calculations have been compared the finite element method. In this respect the code FreeFem++ has been chosen [START_REF] Hecht | New development in FreeFem++[END_REF].

It is an open source finite element code to solve systems of Partial Differential Equations. Its efficiency relies especially on very rapid multi-frontal linear solvers like UMFPACK [START_REF] Davis | Algorithm 832: UMFPACK V4.3 -An unsymmetricpattern multifrontal method[END_REF]. The domain has been discretized by tetrahedral elements with quadratic shape functions (the best 3D-element available in FreeFem++). Various discretizations have been tested and we present the results in Tab. 2.12, the best calculations providing a relative error less than 10 -3 . In this case, the finite element mesh contains 417501 degrees of freedom, what is 20 times more than with the three TMM-procedures. Next we have studied a larger domain l x × l y × d = 20 × 20 × 1 that is split in 400 cubes. In this case the FreeFem++ mesh contains 1656441 degrees of freedom to get an error lower than 10 -3 and the corresponding computation time is 46.44s. With the full least-square TMM and 400 sub-domains, the same accuracy was obtained with 48400 degrees of freedom (degree p = 10) and a computation time of only 25.67s. The performance of full least-square TMM with increase of number of TMM-elements has been collected in Table 2.13. The condition numbers of the final linear system remain nearly constant with the increase of the number of sub-domains. This shows the ability of the present techniques to solve large scale problems with at least the same efficiency as a high-performance finite element code. Up to now, we did not find any limit regarding the number of sub-domains.

Conclusion

New numerical techniques were assessed in this paper to account for boundary conditions and least-square collocation was reconsidered to account for transmission conditions within Taylor Meshless Method. Because TMM solves quasi-exactly the partial differential equations by the technique of Taylor series, the treatment of boundary and interface conditions becomes a central question in view of a wide application of this method. Sometimes a single Taylor series is sufficient to solve the boundary value problem, but in other cases the domain of convergence of the series imposes a splitting in several sub-domains and therefore transmission conditions. The simplest method is the least-square collocation, previously validated for boundary conditions [START_REF] Zézé | A boundary meshless method with shape functions computed from the PDE[END_REF]. In a second technique validated in [START_REF] Tampango | Coupling of polynomial approximations with application to a boundary meshless method[END_REF] for transmission conditions, Lagrange multipliers are introduced and discretized by radial functions. This Lagrange multiplier method has been first applied in this paper to account for boundary conditions. These two techniques are rather simple to be implemented and applied and few discretization parameters are needed: the degrees of the polynomials and a set of collocation points on boundary and interfaces; in the case of Lagrange multipliers, the radial functions are also to be defined from a cloud of points and a characteristic distance.

In most of the cases that were tested in this paper, least-square method is very efficient and exponential convergence with the degree has been obtained. Very few DOFs are needed as compared with classical discretization methods, the ratio being typically in the range 10 ∼ 100. It is also quite robust with respect to the parameters of the algorithm and its only drawback is to lead to full matrices when used with few sub-domains. As for Lagrange multiplier method, it converges about in the same way and is quite robust. Its main drawback is to require additional parameters to define radial functions. In our last 2D elasticity test, we did not find an exponential convergence with the degree, what is very probably due to the fact that the solution is not a smooth function near a corner. The error in this case was much greater than in the previous cases, about 10 -2 instead of 10 -7 , 10 -8 . Thus the method needs improvements to better capture non-smooth solutions which will be discussed in the chapter 5. Further work is also necessary to compute Taylor series for more complex models like nonlinear elasticity [START_REF] Akpama | Amélioration de la "Taylor meshless method" et application à des structures élastiques en grandes transformations[END_REF].

Finally all these methods have been assessed in a 3D case, whose discretization requires hundreds of sub-domains and more than one million of degrees of freedom in an equivalent finite element simulation. Similar performances were found, whatever be the bridging technique (least-square or Lagrange multipliers). All these TMM-variants seem very efficient and robust according to all the benchmarks discussed in this paper.

Chapter 3

Taylor Meshless Method for large-scale problems

Abstract

A true integration-free meshless method based on Taylor series named Taylor Meshless Method (TMM) has been proposed to solve two dimensional Partial Differential Equations (PDE). In this framework, the shape functions are approximated solutions of the PDE and the discretization concerns only the boundary. In this chapter, the applicability of TMM to solve large-scale problems is discussed under two aspects. First, as in some other meshless methods, ill-conditioned matrices and round-off error propagation could lead to a loss of accuracy when the number of unknowns increases. This point will be investigated in the case of large-scale problems. Second the computation time and its distribution are analyzed from numerical experiments for PDE's in a 3D domain. It is established that the TMM method is efficient and robust, even in the case of large-scale problems while the finite element numerical model involves more than three millions degrees of freedom.

Introduction

This chapter analyzes the behavior of Taylor Meshless Method (TMM) in the case of large-scale problems, i.e. for Partial Differential Equation (PDE) whose discretization would require millions of degrees of freedom (DOFs) with a standard technique like Finite Element Method (FEM) or Finite Difference Method (FDM). TMM is based on a quasi-exact resolution of the PDE's by the technique of Taylor series so that the resulting matricial problem comes from the discretization of the boundary conditions [START_REF] Zézé | A boundary meshless method with shape functions computed from the PDE[END_REF]. Two main techniques have been proposed for boundary conditions: least-square collocation [START_REF] Zézé | A boundary meshless method with shape functions computed from the PDE[END_REF] and Lagrange multipliers [START_REF] Yang | Least-square collocation and Lagrange multipliers for Taylor meshless method[END_REF]. The same methods can be also applied for the transmission conditions in piecewise resolutions [START_REF] Tampango | Coupling of polynomial approximations with application to a boundary meshless method[END_REF][START_REF] Yang | Least-square collocation and Lagrange multipliers for Taylor meshless method[END_REF]. The numerical solution is obtained in the form of a single high order polynomial or by a family of polynomials, each one being valid in a sub-domain. An important property of TMM is the exponential convergence with the degree, which permits to obtain very accurate approximations.

There are two challenges in view of applications to large-scale problems: control the computation time, control the condition number of the matrix or at least the effect of ill-conditioning. For the first challenge, controlling computation time, two aspects need to be considered: the way of discretization, the number of DOFs. Two types of discretization techniques, Galerkin method and collocation method are mainly applied in meshless methods. The Galerkin-based meshless methods [START_REF] Belytschko | Element-free Galerkin methods[END_REF][START_REF] Li | A novel complex variable elementfree Galerkin method for two-dimensional large deformation problems[END_REF][START_REF] Zhang | The improved complex variable element-free galerkin method for two-dimensional schrödinger equation[END_REF] always need background grid and integration procedures which lead to expensive computation costs, while the collocation-based methods [START_REF] Kansa | Multiquadrics-a scattered data approximation scheme with applications to computational fluid-dynamics-I surface approximations and partial derivative estimates[END_REF][START_REF] Fairweather | The method of fundamental solutions for elliptic boundary value problems[END_REF][START_REF] Zhang | Least-squares collocation meshless method[END_REF][START_REF] Chen | Numerical convergence of boundary knot method in the analysis of helmholtz, modified helmholtz, and convectiondiffusion problems[END_REF][START_REF] Ferreira | A formulation of the multiquadric radial basis function method for the analysis of laminated composite plates[END_REF], are truly meshless, sometimes very efficient, but difficult to apply for large-scale problems. Nevertheless TMM is quite different from previous meshless techniques due to the convergence properties of Taylor series and to an original analytical solving inside the domain. On the other hand, according to previous numerical experiments, TMM seems competitive in terms of computation time because it needs few DOFs: Tampango et al. [START_REF] Tampango | Convergence analysis and detection of singularities within a boundary meshless method based on Taylor series[END_REF] were able to solve a Laplace problem with less than 100 DOFs while about 5000 were necessary with quadratic interpolation (P2) finite elements and more than 40000 with linear interpolation (P1) triangles at the same level of accuracy. With account of this ratio 5000/100, a large-scale problem involving 500000 DOFs by finite elements corresponds about to 10000 DOFs within TMM. The numerical cost for this class of problems will be analyzed in this chapter for 3D Laplace and elasticity problems, in order to assess the performance of TMM as compared to other discretization techniques and to locate which steps of the algorithm are time consuming.

The effect of a too large condition number is a challenging problem, but it is well documented in the case of mesh-free methods based on a single packet of shape functions, i.e. without splitting into sub-domains. It has been early recognized that the Method of Fundamental Solution converges better when the source points are away from the domain, which means that the method works very well with rather flat shape functions [START_REF] Bogomolny | Fundamental solutions method for elliptic boundary value problems[END_REF]. The same behaviour is observed with Radial Basis Functions, especially with multi-quadrics [START_REF] Kansa | Multiquadrics-a scattered data approximation scheme with applications to computational fluid-dynamics-I surface approximations and partial derivative estimates[END_REF][START_REF] Ferreira | Free vibration analysis of symmetric laminated composite plates by FSDT and radial basis functions[END_REF][START_REF] Ferreira | Computation of natural frequencies of shear deformable beams and plates by an RBF-pseudospectral method[END_REF][START_REF] Driscoll | Interpolation in the limit of increasingly flat radial basis functions[END_REF][START_REF] Cheng | Exponential convergence and H-c multiquadric collocation method for partial differential equations[END_REF] that can be made arbitrarily flat. Even though the use of flat shape functions may decrease the error of approximate solutions, the flatness of shape functions makes them indistinguishable which ultimately leads to the ill-conditioning of the matrices [START_REF] Driscoll | Interpolation in the limit of increasingly flat radial basis functions[END_REF]. This phenomenon is known as "Schaback's uncertainty principle" [START_REF] Schaback | Error estimates and condition numbers for radial basis function interpolation[END_REF]: "the error and the condition number can not both be kept small. Either one goes for a small error and gets a bad sensitivity, or one wants a stable algorithm and has to take a comparably large error." A comprehensive study of this phenomenon can be found in [START_REF] Cheng | Exponential convergence and H-c multiquadric collocation method for partial differential equations[END_REF] by Cheng et al. A similar behaviour has been also observed with Taylor Meshless Method: for small degrees p, the method converges exponentially, generally up to a very high accuracy; then a plateau is reached and beyond this level, the accuracy is progressively affected [START_REF] Tampango | Coupling of polynomial approximations with application to a boundary meshless method[END_REF][START_REF] Yang | Least-square collocation and Lagrange multipliers for Taylor meshless method[END_REF]. Many proposals can be found in the literature to control the effect of ill-conditioned matrices by using localized shape functions [START_REF] Wendland | Piecewise polynomial, positive definite and compactly supported radial functions of minimal degree[END_REF] or sub-domains techniques [START_REF] Kansa | Circumventing the ill-conditioning problem with multiquadric radial basis functions: applications to elliptic partial differential equations[END_REF][START_REF]Domain decomposition and local refinement for multiquadric approximations. I: second-order equations in one dimension[END_REF][START_REF] Wong | Multizone decomposition for simulation of time-dependent problems using the multiquadric scheme[END_REF]. For instance according to recent numerical experiments [START_REF] Chu | Finite subdomain radial basis collocation method[END_REF], it appears that the sub-domain technique improves the conditioning, reduces the computation time and delays the loss of stability due to an ill-conditioned matrix, but can not overcome Schaback's uncertainty principle. Such a sensitivity to conditioning is rarely mentioned with other meshless methods based on more localized approximations [START_REF] Belytschko | Element-free Galerkin methods[END_REF][START_REF] Nayroles | Generalizing the finite element method: diffuse approximation and diffuse elements[END_REF][START_REF] Nguyen | Meshless methods: a review and computer implementation aspects[END_REF][START_REF] Liew | Nonlinear analysis of corrugated plates using a FSDT and a meshfree method[END_REF][START_REF] Li | A novel complex variable elementfree Galerkin method for two-dimensional large deformation problems[END_REF]. In this chapter, the effect of matrix conditioning on the stability of TMM will be discussed carefully from numerical experiments in the case of large-scale problems.

The present study does not fully exhaust the question of applicability of TMM for practical engineering problems. The treatment of non-linear problems has been sketched in [START_REF] Zézé | Calcul de fonctions de forme de haut degré par une technique de perturbation (Computing high degree shape functions by a perturbation technique[END_REF][START_REF] Tampango | Développement d'une méthode sans maillage utilisant les approximations de Taylor (Development of a meshless method using Taylor approximations[END_REF][START_REF] Akpama | A Taylor meshless method for hyperelasticity[END_REF] by associating Newton-like algorithms with TMM-discretization and Automatic Differentiation [START_REF] Griewank | Evaluating higher derivative tensors by forward propagation of univariate Taylor series[END_REF]. Another question is posed by the performance of TMM for non-smooth solutions that are common in engineering computations because of corners, cracks, concentrated forces and so on. These two questions will be re-discussed in details in the next chapters. On the contrary, previous papers have established that the method works well for any boundary condition and any domain geometry with an easy implementation.

This chapter is organized as follows: section 3.2 states the algorithm for TMM and its piecewise resolutions. Our numerical tests will be presented in section 3.3. We analyze the influence of the number of collocation points in section 3.4.1 and the robustness of convergence within a single domain in section 3.4.2 and piecewise resolutions in section 3.4.3. The effect of a large condition number of matrix will be discussed in section 3.4.4. The computation time for each step of the algorithm and its comparison with FEM are presented in section 3.5, especially in section 3.5.3 where one discusses a case needing three millions DOFs with FEM.

Algorithm for Taylor Meshless Method

In this section, the Taylor Meshless Method is adapted in a 3D framework. It is based on a resolution of the PDE's from Taylor series. One seeks the unknowns as rather high degree polynomials by vanishing the residual of the PDE's up to a given order. This leads to a general solution of the PDE's that depends on many parameters, but these unknowns are much less numerous than with other discretization techniques. Typically the general solution of 2D Laplace equation ∆u = 0 is a linear combination of the polynomials Re(x + iy) n and Im(x + iy) n that are the shape functions of the method. Here we extend in the 3D case the method of [START_REF] Zézé | A boundary meshless method with shape functions computed from the PDE[END_REF] that brings together all the shape functions having a same degree. Next the boundary value problem is solved by leastsquare method as proposed in [START_REF] Zézé | A boundary meshless method with shape functions computed from the PDE[END_REF]. If necessary, a piecewise resolution can be done either by introducing Lagrange multipliers [START_REF] Tampango | Coupling of polynomial approximations with application to a boundary meshless method[END_REF] or also by least-square method [START_REF] Yang | Least-square collocation and Lagrange multipliers for Taylor meshless method[END_REF]. All these techniques will be recalled shortly for completeness.

Algorithm to compute the shape functions

In this section, we limit ourselves to the case of a linear equation with constant coefficients and a general right hand side f (x):

-∆u(x) + u(x) = f (x), x ∈ Ω (3.1)
where Ω is a bounded domain in R 3 . The unknown function u(x) is expanded in Taylor series with a development point c = [0, 0, 0]:

u(x, y, z) = p k=0 k i=0 k-i j=0 u j,i,k-i-j x j y i z k-i-j (3.2)
The idea is to build the affine subspace of the polynomials that solves the PDE (3.1) in the sense of Taylor series. For convenience, each polynomial is split into homogeneous polynomials. A homogeneous polynomial of degree k combines (k + 2)(k + 1)/2 monomials and it is written in a matricial form

x k { u k }, where u k =              { u k 0 } { u k 1 } { u k 2 } . . . { u k k }              , u k i =              u k-i,i,0 u k-i-1,i,1 u k-i-2,i,2 . . . u 0,i,k-i              . x k = x k 0 , x k 1 , x k 2 , • • • , x k k x k i = x k-i y i , x k-i-1 y i z 1 , x k-i-2 y i z 2 , • • • , y i z k-i
In this way, all the polynomials whose degree is lower or equal to p can be written in the following form:

u(x, y, z) = p k=0 x k { u k } = p k=0 u k (x, y, z) (3.3)
The basis of monomials within a homogeneous polynomial of degree k and the corresponding coefficients are organized as follows:

u k,0,0 u k-1,1,0 u k-1,0,1 u k-2,2,0 u k-2,1,1 u k-2,0,2 u k-3,3,0 u k-3,2,1 u k-3,1,2 u k-3,0,3 . . . . . . . . . . . . . . . u 0,k,0 • • • u 0,3,k-3 u 0,2,k-2 u 0,1,k-1 u 0,0,k x k x k-1 y 1 x k-1 z 1 x k-2 y 2 x k-2 y 1 z 1 x k-2 z 2 x k-3 y 3 x k-3 y 2 z 1 x k-3 y 1 z 2 x k-3 z 3 . . . . . . . . . . . . . . . y k • • • y 3 z k-3 y 2 z k-2 y 1 z k-1 z k
The red color terms represent the terms {x k 0 } and { u k 0 } that do not contain y (i = 0). The green color items represent the terms {x k 1 } and { u k 1 } that contain y 1 (i = 1). If the PDE (3.1) is considered as a differential equation with respect to y, its solution can be expressed in terms of the "initial values" u(x, 0, z) and ∂u ∂y (x, 0, z). That is why we have distinguished the terms i = 0 and i = 1. The corresponding coefficients of these terms i = 0 and i = 1 are chosen as independent variables and collected in the following vector:

{ v k } = { u k 0 } { u k 1 } ∈ R 2k+1 (3.4)          ∈ R k(k-1)/2 (3.5)
Two matrices [S vk ] and [S wk ] can be constructed (see Appendix A) to split the coefficients vector { u k } into two parts:

{ u k } = [S vk ]{ v k } + [S wk ]{ w k } (3.6)
Next, if u(x, y, z) is a homogeneous polynomial of degree k + 2, ∆u is a homogeneous polynomial of degree k. Hence, there exists a matrix

[L k ] such that ∆u k+2 (x) = x k [L k ]{ u k+2 } (3.7)
The matrix [L k ] represents the action of the Laplacian on the coefficients of homogeneous polynomials. It is defined in Appendix A Eq. ( 1). The principle is to vanish the residual of ∆uu + f , up to the degree p -2. With account of Eq. (3.7), this leads to the following equation:

x k [L k ]{ u k+2 } -{ u k } + { f k } = 0, k ∈ [0, p -2] (3.8)
By inserting Eq. (3.6) into Eq. (3.8), we get

[L k ][S vk+2 ]{ v k+2 }+[L k ][S wk+2 ]{ w k+2 } = [S vk ]{ v k } + [S wk ]{ w k } -{ f k } (3.9)
One checks easily that the square matrix

[L k ][S wk+2 ] given in the Appendix A, is invertible. Assuming that the vector { v k }(k ∈ [0, p]
) is known, the Eq. (3.9) becomes the relation of recurrence between { w k+2 } and { w k }. Then one can get the expression of

{ w k }(k ∈ [0, p]) as a function of { v k }(k ∈ [0, p]
). Finally, one can get the approximated general solution of the Laplace equation as follows:

u(x) = P 0 (x) + (p+1) 2 i=1 P i (x)v i = P 0 (x) + P (x) { v} (3.10)
where the P 0 (x, y, z) is used to balance the right-hand side of Eq. (3.1) and the unknown vector { v} is defined by

{ v} =          { v 0 } { v 1 } . . . { v p }          ∈ R (p+1) 2 (3.11)

Boundary least-square collocation

Next one has to apply the boundary conditions. It is well established that a pure collocation does not converge. As proposed in [START_REF] Zézé | A boundary meshless method with shape functions computed from the PDE[END_REF], the boundary conditions will be accounted by a least-square collocation procedure that was previously used in other meshless context, for instance by Fairweather and Karageorghis [START_REF] Fairweather | The method of fundamental solutions for elliptic boundary value problems[END_REF] and by Zhang et al. [START_REF] Zhang | Least-squares collocation meshless method[END_REF]. One chooses a set of nodes x j , j ∈ [1, M ], on the boundary of the 3D domain and one minimizes the error between the approximate value and the given value of u in these points. Here we consider a mixed boundary conditions as follows:

u(x) = u d , x ∈ Γ d ; Tu(x) = t n , x ∈ Γ n . (3.12)
It comes to minimize the following function with respect to the components of the vector {v}:

T (v) = 1 2 x i ∈Γ d u p (x i ) -u d (x i ) 2 + w • 1 2 x j ∈Γn |Tu p (x j ) -t n (x j )| 2 (3.13)
where Γ d represents the Dirichlet boundary condition, and Γ n represents the Neumann boundary. This minimization leads to a linear system [K]{v} = {b} where [K] is an invertible matrix. Solving this system gives the vector {v} and therefore the approximate solution of the boundary value problem Eq. (3.1) with boundary Eq. (3.12). Note the simple form of the matrix (α, β ∈ [1, (p + 1) 2 ]):

[K] αβ = M d i=1 P α (x i ) • P β (x i ) + w • Mn j=1 TP α (x j ) • TP β (x j ), (3.14
)

{b} α = M d i=1 P α (x i ) u d (x i ) -P 0 (x i ) + w • Mn j=1 TP α (x j ) [t n (x j ) -TP 0 (x j )] . (3.15)

Piecewise resolution

It is not possible to solve any boundary value problem with a single series. Thus, numerical methods have to be introduced to connect several high-order polynomial approximations. Within this aim, several techniques have been presented and well investigated in the previous works for two dimensional problems: discrete Lagrange multipliers or discretized continuous Lagrange multipliers [START_REF] Tampango | Coupling of polynomial approximations with application to a boundary meshless method[END_REF] and least-square collocation method [START_REF] Yang | Least-square collocation and Lagrange multipliers for Taylor meshless method[END_REF]. The aim of this section is to briefly recall the main procedures of these bridging techniques in 3D piecewise resolutions.

Here we consider a domain as shown in Fig. 3.1 which has been divided into two parts Ω 1 and Ω 2 . For each sub-domain, an independent high-order polynomial approximation based on Taylor series, u 1 (x) and u 2 (x), will be obtained as stated in section 3.2.1. The functions accounting for boundary conditions of each sub-domain, T 1 (v 1 ) and T 2 (v 2 ), are as in Eq. (3.13). In order to connect the two approximations, a set of M r collocation points x j on the interface Γ r are needed for the discretization of the transmission conditions. The constraints on the interface are as the following:

ψ 1 j (v 1 , v 2 ) = u 1 (x j ) -u 2 (x j ) = 0 ψ 2 j (v 1 , v 2 ) = ∂u 1 ∂n (x j ) -∂u 2 ∂n (x j ) = 0 ∀x j ∈ Γ r (3.16)
In the following, three methods for the transmission conditions are illustrated respectively.

Least-square collocation method

By applying the least-square collocation method, the transmission conditions (3.16) will be satisfied in a mean sense via the following coupling function:

C ls (v 1 , v 2 ) = 1 2 Mr j=1 u 1 (x j ) -u 2 (x j ) 2 + ∂u 1 ∂n (x j ) - ∂u 2 ∂n (x j ) 2 (3.17)
Then for the whole problem, searching the variables {v 1 } and {v 2 } is equivalent to minimizing the following function:

T (v 1 , v 2 ) = T 1 (v 1 ) + T 2 (v 2 ) + C ls (v 1 , v 2 ) (3.18)
The minimization of this function T (v 1 , v 2 ) leads to a linear system. Solving this system gives the vectors {v 1 } and {v 2 } and then the approximate solution for each sub-domain.

Discrete Lagrange multipliers

Accounting for the discrete transmission conditions (3.16) by Lagrange multipliers, the constraints will be satisfied in a strong form by using two Lagrange multipliers on each collocation point on the interface via the following coupling function:

C lm (v 1 , v 2 , λ 1 , λ 2 ) = Mr j=1 [λ 1 j ψ 1 j (v 1 , v 2 ) + λ 2 j ψ 2 j (v 1 , v 2 )] (3.19)
Then for the whole problem, searching the variables {v 1 } and {v 2 } is equivalent to setting to zero the gradient of the following function:

T (v 1 , v 2 , λ 1 , λ 2 ) = T 1 (v 1 ) + T 2 (v 2 ) + C lm (v 1 , v 2 , λ 1 , λ 2 )
(3.20)

Continuous Lagrange multiplier and its discretization

Another way is to define the Lagrange multiplier at the continuous level. A bilinear form is introduced to account for the transmission conditions:

C ar (λ, u) = Mr i=1 (λ(x i )µ(x i ) + l 2 ∇λ(x i )∇µ(x i )) (3.21)
where l is a characteristic length and λ(x) is the continuous version of the Lagrange multiplier. In Eq. (3.21) a set of M r collocation points x i on the interface has always to be chosen. This bilinear form looks like the bilinear form of Arlequin method [START_REF] Lions | Problèmes mécaniques multi-échelles: la méthode Arlequin[END_REF][START_REF] Hu | A bridging technique to analyze the influence of boundary conditions on instability patterns[END_REF], but for consistency with the meshless framework, it is defined from points and not from an integral. Numerical experiments in [START_REF] Tampango | Coupling of polynomial approximations with application to a boundary meshless method[END_REF] have shown that a distribution of collocation points in a volume is not necessary. So the points are located along the interface between the sub-domains. Next the continuous Lagrange multiplier λ(x) has to be discretized. As in [START_REF] Tampango | Coupling of polynomial approximations with application to a boundary meshless method[END_REF], this discretization is done from two families of radial functions depending of another set of N λ discretization nodes x j (that are the centers of the radial functions):

λ(x) = N λ j=1 (λ 1 j Φ 1 j (x) + λ 2 j Φ 2 j (x)) (3.22) Φ 1 j (x) = exp( -x -x j 2 d 2 ) and Φ 2 j (x) = (x -x j ) • n exp( -x -x j 2 d 2 )
(3.23) where d represents the radius of influence of the radial functions and n is the unit normal to the interface. Then solving the whole problem in the domain is equivalent to taking the stationary values of the following function:

T (v 1 , v 2 , Λ 1 , Λ 2 ) = T 1 (v 1 ) + T 2 (v 2 ) + C ar (Λ 1 , Λ 2 , u 1 -u 2 ) (3.24)
As compared with the previous discrete Lagrange multipliers method, it permits to use less degrees of freedom (N λ < M r ).

Numerical examples

This paper focuses on the behavior of Taylor Meshless Method for large-scale problems. That is why numerical tests in 3D domains are considered and discussed carefully, what was not done before. Only cases having an exact solution are studied in this chapter, what permits to evaluate the error even it is very small. Throughout this paper, the relative error is the difference between exact and approximate solutions, divided by the maximum value of the solution:

E = max |u(x ) -u ex (x )| max |u(x )| . ( 3.25)
Three boundary value problems will be tested in this paper, two concerning the Laplace equation and the last one the equations of linear isotropic elasticity. In the first case, the exact solution is a polynomial of small degree. This permits to check the consistency of the method: the exact solution lying in the subspace of our shape functions, is it easily recovered within TMM ? In the next two tests, the exact solution is a fundamental solution. According to the location of the singularity, this function can be more or less flat or strongly varying. In the last test, the solution contains many waves so that the discretization involves a very large number of DOFs.

Laplace equation with polynomial solution

The first example is the Laplace equation in a spherical domain (Ω = {x, y, z|x 2

+ y 2 + z 2 ≤ 1}) -∆u = 0 in Ω u(x, y, z) = x 3 yz + xy 2 z on ∂Ω (3.26)
The exact solution is a polynomial of degree 5:

u ex = x 3 yz + xy 2 z + ( 1 3 xyz + 1 7 xz)(1 -x 2 -y 2 -z 2 ) (3.27)

Laplace equation with singular solution

The second example is always a Laplace equation with alternative Dirichlet boundary conditions in the unit sphere (Ω = {x, y, z|x

2 + y 2 + z 2 ≤ 1})    -∆u = 0 in Ω u(x, y, z) = 1 r on ∂Ω (3.28)
where

r 2 = (x -x 0 ) 2 + (y -y 0 ) 2 + (z -z 0 ) 2 .
The exact solution is

u ex = 1 (x -x 0 ) 2 + (y -y 0 ) 2 + (z -z 0 ) 2 (3.29)
The solution has a singularity at X 0 = [x 0 , y 0 , z 0 ] that must lie outside the domain. If this singularity is far away from the domain, the solution is rather flat. If it is closer, the solution involves large gradients, what shall require finer discretizations.

3D elasticity

The previous examples concern a single elliptic equation. In view of practical applications, it is important to discuss the efficiency of the method for systems of PDEs. Linear isotropic 3D elasticity in the absence of body forces is considered in the unit sphere. The unknown fields are the components u 1 , u 2 and u 3 of the displacements. In Cartesian coordinates, the model is based on the Cauchy-Navier equations:

2 -2ν 1 -2ν ∂ 2 u 1 ∂x 2 1 + ∂ 2 u 1 ∂x 2 2 + ∂ 2 u 1 ∂x 2 3 + 1 1 -2ν ∂ 2 u 2 ∂x 1 ∂x 2 + 1 1 -2ν ∂ 2 u 3 ∂x 1 ∂x 3 = 0 (3.30) 2 -2ν 1 -2ν ∂ 2 u 2 ∂x 2 2 + ∂ 2 u 2 ∂x 2 1 + ∂ 2 u 2 ∂x 2 3 + 1 1 -2ν ∂ 2 u 1 ∂x 1 ∂x 2 + 1 1 -2ν ∂ 2 u 3 ∂x 2 ∂x 3 = 0 (3.31) 2 -2ν 1 -2ν ∂ 2 u 3 ∂x 2 3 + ∂ 2 u 3 ∂x 2 1 + ∂ 2 u 3 ∂x 2 2 + 1 1 -2ν ∂ 2 u 1 ∂x 1 ∂x 3 + 1 1 -2ν ∂ 2 u 2 ∂x 2 ∂x 3 = 0 (3.32)
The studied benchmark is the Dirichlet problem, where the exact solution is a fundamental solution. For a source located at a point Q acting at a point P , the fundamental solutions of the system (3.30, 3.31, 3.32) are (see, e.g., Refs. [START_REF] Poullikkas | The method of fundamental solutions for three-dimensional elastostatics problems[END_REF]):

G 11 (P, Q) = 1 16πµ(1 -ν) (3 -4ν)r 2 P Q + (x 1P -x 1Q ) 2 r 3 P Q (3.33) G 12 (P, Q) = 1 16πµ(1 -ν) (x 1P -x 1Q )(x 2P -x 2Q ) r 3 P Q = G 21 (P, Q) (3.34) G 13 (P, Q) = 1 16πµ(1 -ν) (x 1P -x 1Q )(x 3P -x 3Q ) r 3 P Q = G 31 (P, Q) (3.35) G 22 (P, Q) = 1 16πµ(1 -ν) (3 -4ν)r 2 P Q + (x 2P -x 2Q ) 2 r 3 P Q (3.36) G 23 (P, Q) = 1 16πµ(1 -ν) (x 2P -x 2Q )(x 3P -x 3Q ) r 3 P Q = G 32 (P, Q) (3.37) G 33 (P, Q) = 1 16πµ(1 -ν) (3 -4ν)r 2 P Q + (x 3P -x 3Q ) 2 r 3 P Q (3.38)
where

r P Q = (x 1P -x 1Q ) 2 + (x 2P -x 2Q ) 2 + (x 3P -x 3Q ) 2 (3.39)

A very large-scale test

This benchmark is issued from the study of instabilities of film-substrate systems [START_REF] Chen | Herringbone buckling patterns of compressed thin films on compliant substrates[END_REF][START_REF] Xu | 3D finite element modeling for instabilities in thin films on soft substrates[END_REF], where one has to describe many spatial oscillations. Here one accounts only for the substrate and the film is represented by a prescribed sinusoidal displacement. The domain is a cuboid with sides l x × l y × d (d = 1 in any case), see Fig. 3.2. For simplicity we consider only the Laplace problem ∆u = 0 with Dirichlet boundary conditions (u(x, y, d) = sin(πx)sin(πy) on the top face, u(x, y, 0) = 0 on the bottom), the exact solution being u(x, y, z) = sin(πx)sin(πy)sinh( √ 2πz)/sinh( √ 2πd). The domain is split in sub-domains that are rectangular hexahedra of size 1 × 1 × 1 that will be called "elements" for simplicity, see Fig. 3.3. This corresponds to one element through the thickness and two elements along each period. The connection between sub-domains is achieved by the least-square method and a rather small degree (p = 10) is used.

Convergence and conditioning

In previous papers [START_REF] Zézé | A boundary meshless method with shape functions computed from the PDE[END_REF][START_REF] Tampango | A new meshless method using Taylor series to solve elasticity problems[END_REF][START_REF] Tampango | Convergence analysis and detection of singularities within a boundary meshless method based on Taylor series[END_REF][START_REF] Tampango | Coupling of polynomial approximations with application to a boundary meshless method[END_REF] and in chapters 1 and 2, many results were presented, especially the influence of the number of collocation points, the convergence with the degree, the splitting into sub-domains, the account of various boundary conditions, of various domain shapes and globally the robustness of the procedure. For consistency of the paper, some points will be re-discussed here, but not the last one's: we limit ourselves to Dirichlet boundary conditions and to very simple geometries, the sought conclusion being the ability of the method to deal with large systems. In addition to the issues mentioned above, there exists another key point: the influence of round-off errors and matrix conditioning. Indeed looking for very accurate solutions leads to ill-conditioned matrix and finally to loss of accuracy [START_REF] Cheng | Exponential convergence and H-c multiquadric collocation method for partial differential equations[END_REF]. Thus it is an important issue to check whether these problems prevent the use of TMM to solve large-scale problems.

Influence of the number of collocation points

The influence of the number of collocation points has been illustrated in Fig. 3.4, Fig. 3.5 and Fig. 3.6 for the polynomial solution of the Laplace problem Eq. (3.26), the singular solution of the Laplace problem Eq. (3.28) and for the elasticity test Eq. (3.30, 3.31, 3.32) respectively, and for three values of the degree p = 10, 20 and 30. Let us recall that the number of unknowns is (p + 1) 2 for a single equation and 3(p + 1) 2 in 3D elasticity. In the three cases, the method diverges for the pure collocation or when the number of collocation points is only slightly larger than the number of unknowns. Next a maximal accuracy is obtained for about M = 2(p+1) 2 and it remains constant beyond this threshold. A similar behavior has been previously observed in 2D cases [START_REF] Zézé | A boundary meshless method with shape functions computed from the PDE[END_REF]. So the method is remarkably stable with respect to the number of collocation points. In what follows, we choose a number of collocation points M = 2(p+2) 2 that is slightly larger than two times of the number of unknowns (p + 1) 2 necessary for a single equation, what turns out to be sufficient in the considered applications.

Of course the obtained maximal accuracy depends on the problem and of the degree of the shape functions. In the case of a polynomial solution of degree 5, see Fig. 3.4, this accuracy is very large, up to 10 -15 , what is very close to the unit round-off error 2 -53 ≈ 1.11e -16 for a single real number within the double-precision floating-point format [START_REF] Muller | Handbook of Floating-Point Arithmetic[END_REF]. Because the exact solution lies in the subspace generated by the shape functions, this very small error is only due to the propagation of these round-off errors. Hence it is natural that the final error increases with the degree because this induces a larger number of elementary operations. In the two other cases, see Fig. 3.5 and Fig. 3.6, the exact solution is not a polynomial. The maximal accuracy increases with the degree, from about 10 -3 for p = 10 to 10 -8 for p = 30, what corresponds to the exponential convergence of any Taylor series. In those cases, the influence of the round-off errors is too small to be observed in the range p ∈ [START_REF] Lancaster | Surfaces generated by moving least squares methods[END_REF][START_REF] Kansa | Multiquadrics-a scattered data approximation scheme with applications to computational fluid-dynamics-II solutions to parabolic, hyperbolic and elliptic partial differential equations[END_REF]. This balance between round-off errors and Taylor series errors is the heart of the present discussion.

Exponential convergence

The convergence with the degree is summarized in Fig. 3.7 for the singular solution of Laplace equation (3.28). Two cases are presented, a first one where the singularity is away from the domain, X 0 = [2, 2, 2], and therefore the solution is rather flat, a second one where the singularity X 0 = [1, 1, 1] is closer and the gradient is larger. The rate of convergence is better if the singularity is far away from the domain, i.e. the convergence is faster for the flat solution, but in the two cases the convergence is exponential for a small degree. A second stage can be seen from a degree p=25, from which the accuracy is no longer improved but it deteriorates slightly. It seems that the convergence is significantly affected by the propagation of round-off errors for large values of the degree, while before this level the behavior is dominated by the convergence of the Taylor series. This scheme is not new. Within TMM, it had been observed in [START_REF] Tampango | Convergence analysis and detection of singularities within a boundary meshless method based on Taylor series[END_REF] Figure 3, or in [START_REF] Tampango | Coupling of polynomial approximations with application to a boundary meshless method[END_REF] Figures 10 and14, even in cases of piecewise resolutions. It is also similar to "Schaback's uncertainty principle" discussed in [START_REF] Schaback | Error estimates and condition numbers for radial basis function interpolation[END_REF] for meshless method based on radial basis functions. In the latter paper, this influence of round-off errors becomes also significant beyond some level of accuracy and this had been attributed to the ill-conditioning of the matrix. This point will be discussed below.

Next the elasticity problem with the nearest source point X 0 = [1, 1, 1] is discussed and some features are reported in Table 3.1. Once again, a reversal of the accuracy curve is found for a degree p = 40, which corresponds to a matrix condition number of about 10 28 . Nevertheless the accuracy remains very high after the reversal point (p = 50, 60) as well as before (p = 30), with an error always lower than 10 -6.6 . This illustrates the robustness and the high precision of the method, whether before or after the accuracy reversal.

Piecewise resolutions

Piecewise resolutions are necessary because Taylor series often diverge [START_REF] Tampango | Coupling of polynomial approximations with application to a boundary meshless method[END_REF]. So it is important that the method keeps the same properties of convergence and of robustness in cases of piecewise resolutions. The Laplace problem Eq. (3.28) with singular point X 0 = [1, 1, 1] is solved by the three techniques mentioned in section 3.2.3. The one-domain method works in this case, see Fig. 3.7, and the test will permit to check the influence of the matching technique on the features of the discretization. The sphere is split in two sub-domains by the xy-plane. In the following, we choose M i = 2(p + 1) 2 collocation points on the boundary of each sub-domain and M r = 2p 2 + 1 collocation points on the interface. For piecewise resolutions coupled by Lagrange multipliers, the number of Lagrange multipliers is two times the number of the collocation points on the interface (N λ = 4p 2 + 2). For the weak sense coupling technique (Arlequin type, section 3.2.3.3), there are N λ = 4(p/2 -1) 2 + 2 Lagrange multipliers on the interface and the radius of influence of the radial functions is chosen to be d = 0.1.

The three matching techniques work well, with about the same characteristics as the one-domain method, see Table 3.2, 3.3, 3.4: first there is always a turning point in the accuracy curve for p ≈ 30, what corresponds to an accuracy in the range [10 -7 , 10 -6 ]; second the condition number of the matrix [K] has about the same value 10 27 at this accuracy reversal so that the equilibrium point between round-off errors and Taylor-series errors seems very strongly related to the matrix conditioning. They are minor differences between the three matching techniques but they are not important by comparison with other criteria, as the simplicity of least-square and Lagrange multiplier techniques, the lower number of dof's with least-square method or the sparsity of the matrix with the weak sense coupling technique. Anyway the three matching methods are available and none degrades the convergence properties of the one-domain TMM.

More about conditioning

To go further in the analysis of conditioning, many calculations have been done, for the three benchmarks described in section 3.3, for many degrees in the range [5,[START_REF] Golberg | The method of fundamental solutions for potential, helmholtz and diffusion problems[END_REF]. In the cases of singular solutions (Eq. (3.28) and Eq. (3.30, 3.31, 3.32)), many locations of the source points were considered, from very flat solutions for X 0 = [100, 100, 100] up to source points rather close to the domain for X 0 = [1, 1, 1]. When the solution is pushed to higher accuracy by increasing the degree p of Taylor series, the matrix condition number increases. For instance in the case of Eq. (3.26), the condition number changes from 10 3.5 to 10 37.6 , (Fig. 3.8). Because of the propagation of roundoff errors, the accuracy decreases with the degree by starting from a very high accuracy (error = 10 -15.1 for p = 5). A remarkable feature is the quasilinearity of the curve in the Log-Log plot that can be approximately given as y ≈ 0.206x -16.4, where the notations y and x represent the Log 10 (maximal error) and Log 10 (condition number) respectively, the value at the origin being very close to the unit round-off error in double precision (1.1×10 -16 ). The very high accuracy for a low degree p = 5 is due to fact that the exact solution is a polynomial and the number of elementary calculations is small at this level. By increasing the degree beyond p = 5, the number of operations increases and this induces further errors without any compensation due to Taylor series that can be quasi-exact from p = 5.

The behavior is different in the general case of non-polynomial solutions, except for very flat solutions (X 0 = [100, 100, 100]), in which case the accuracy is also very high from p = 5. In the other cases (see Fig. 3.9 and 3.10), one recovers the behavior described in Sections 3.4.1 and 3.4.2. On the left part, one obtains an exponential convergence with the degree and a more or less affine relation between error and conditioning. This corresponds to a convergence dominated by the Taylor series. On the right part of these figures, convergence is deteriorated by the increase of the degree and therefore dominated by the propagation of round-off errors. There is another remarkable feature: all these accuracy curves tend asymptotically to straight lines whose equations are reported in the figures. Moreover these straight lines are very close in the three considered boundary value problems (Laplace equation with a polynomial solution, Laplace equation with source points and 3D elasticity) and this could be approximated by y = 0.2x -16, see Fig. 3.11. Of course the value at the origin (-16) is connected to unit round-off error in double precision, but the slope (0.2) seems also more or less universal to characterize the loss of accuracy due to ill-conditioning. According to this prediction y = 0.2x -16, the accuracy remains fabulous (y ≈ 10 -8 ∼ 10 -6 ) for condition number of matrix up to 10 40 ∼ 10 50 (x = 40 ∼ 50) as degree up to p = 65 ∼ 80. Nevertheless it is not necessary to choose a so large degree, an optimal degree will be at the reversal or before.

In summary, the accuracy of TMM is subjected to two antagonistic effects: the convergence of Taylor series that decreases the error when the degree p increases and the propagation of round-off errors that grows with p. The second effect yields a maximal value of the degree that is very large. The combination of the two effects leads to an accuracy reversal so that an optimal computation requires a degree before this reversal.

Computation time

The computation time of TMM has not yet been seriously analyzed in the literature. In 2D cases, a comparison with FEM has been done in [START_REF] Tampango | Convergence analysis and detection of singularities within a boundary meshless method based on Taylor series[END_REF] and it has established that TMM permits a very strong reduction of the number of degrees of freedom as compared with FEM, typically less than 100 DOFs instead of 5000 with FEM at the same level of accuracy. In this section, the origin of the computational cost within TMM is discussed in details. A comparison with the FEM code FreeFem++ [START_REF] Hecht | New development in FreeFem++[END_REF] will be carried out, mainly to locate the time consuming steps of the TMM computations. This section can be considered as a first attempt to analyze the computational cost of TMM. It has been done with a commercial software (MATLAB R2012b) in a personal computer under the system of Windows 7 with intel(R) core(TM) i5-2300 type CPU and 8 GB computer memory.

Analysis of the computation time

We considered the Laplace problem Eq. (3.28) with a source point X 0 = [1, 1, 1]. This problem is solved by TMM with a development point at c = [0, 0, 0] and a sufficiently large number of collocation points. The convergence analysis is similar with the previous ones, there is an accuracy reversal at p = 40 and a high optimal accuracy (10 -10.2 ), in Table 3.5. Globally with this new case, one recovers about the same features as with the elasticity problem of Table 3.1 and the multi-domain calculations of Table 3.2, 3.3 and 3.4: same order of magnitude of the errors and of the CPU time, similar accuracy curves with a reversal around p = 30 or 40, same growth of the CPU time with the degree. The most significant difference is a larger computation time for the 2 sub-domains cases because they use more DOFs. In the present case, a sufficient accuracy (10 -5 ) is obtained with a rather low degree p = 20, few degrees of freedom (441) and a very small computation time (0.1s). Few details about the CPU time are given in Table 3.6. It turns out that a large part of CPU time (71% to 86%) is due to the calculation of the polynomial shape functions, more precisely the numerical values of these polynomials at the cloud of collocation points P i (x j ). A small part of the time is devoted to the calculation of the matrix K, i.e. of products of these quantities P i (x j ), see Eq. (3.14-3.15). The inversion of the matrix requires few time, less than 6%, even with the largest matrix (p = 60, 3721 DOFs). Note that the case p=10 is different because the computer needs a preparing work that requires about 21.5% of the time, but this small size case has little significance in the present discussion. This time distribution is quite different with classical discretization methods where the main cost is due to matrix inversion, what is due to a much larger number of unknowns. A first consequence of these features is an evaluation of the computation time as a function of the degree: the number of polynomials corresponds to the number of unknowns (p + 1) 2 and it has been seen that the number of collocation points is also proportional to this number of DOFs. Hence the number of quantities to be computed P i (x j ) is of the order of p 4 . In Fig. 3.12, we plotted the experimental CPU time as a function of the degree and we found that is about O(p 4.04 ), which has a very good agreement compared with the predicted order O(p 4 ). There is a second interesting consequence of the time distribution: the calculation of all the quantities P i (x j ) does not require a sequential treatment. Hence this main task in TMM can be achieved in parallel. The used processor has four cores so that first parallel computations have been achieved, the results being reported in Table 3.1 and Table 3.6. Significant time saving have been obtained for the most expansive calculations (40 ≤ p ≤ 60), with a gain between 35 and 45% for the results of Table 3.6 and between 15 and 25% in Table 3.1.

First comparison with FEM

An assessment of the efficiency of TMM requires a first comparison with a classical discretization method. In this respect, we have chosen FreeFem++ that is an open source finite element code to solve systems of Partial Differential Equations [START_REF] Hecht | New development in FreeFem++[END_REF]. It has a lot of users and offers a large variety of triangular elements, an automatic mesh generator and a large variety of linear solvers.

Here tetrahedrons with linear (T1) or quadratic (T2) interpolation will be used to solve the Laplace problem Eq. (3.28) in a spherical domain corresponding to singularities X 0 = [1, 1, 1] and X 0 = [2, 2, 2]. To compare TMM and FreeFem++, we look at the number of degrees of freedom and the computation time to get a given accuracy, namely a relative error less than 10 -2 , 10 -3 or 10 -4 . In each case, several finite element meshes have been tested and the most significant two have been reported in Table 3.7. With the linear tetrahedron, we were not able to reach the highest accuracies because of a lack of memory in our computer.

In any case the meshless method is more efficient than the FEM, of course in terms of number of unknowns, but also of computation time. The best results with FreeFem++ correspond to the cases with a weak accuracy requirement (error less than 1%), with a time ratio 32/1 and 100/1 and a DOFs ratio of 6/1 and 13/1. The results with a strong accuracy (error less than 10 -3 or than 10 -4 ) are more in favor of the meshless method that converges more rapidly with an increase of the degree, the corresponding ratios being respectively 50/1 and 140/1 (time) and 22/1 and 70/1 (DOFs).

Hence Taylor Meshless Method seems a competitive approach to reduce significantly the number of unknowns and the computation time as compared with well established discretization methods. Of course these first results would deserve to be deepened with larger scale problems, other tests with several subdomains, more complex PDEs • • • , a first attempt is presented in the next section.

Large boxes submitted to sinusoidal loading

In most of the previous examples, the focus was on the effect of matrix conditioning, on the robustness of the algorithm and the evaluation of computation cost. It was established that the algorithm still give rather accurate results with very large degrees [START_REF] Tampango | Convergence analysis and detection of singularities within a boundary meshless method based on Taylor series[END_REF][START_REF] Golberg | The method of fundamental solutions for potential, helmholtz and diffusion problems[END_REF], but the optimal degree lies before the accuracy reversal presented in Fig. 3.8, 3.9 and Fig. 3.10. Moreover large degrees induce increasing computational costs of the order O(p 4 ). At the opposite, it does not seem a good idea to settle for a small degree (p = 1, 2, • • • ) as in the FEM. Indeed, first, the solution by Taylor series and the absence of numerical Table 3.7: Comparison between the proposed method and FEM: Tetrahedrons with linear and quadratic shape functions (T1 and T2). The table presents the number of NDOFs necessary to get a relative error lower than 10 -2 , 10 -3 , 10 -4 , for problem Eq. (3.28) with two source points X 0 . In parentheses, the actually obtained error.

Element Type

Singularity at integration permits to work with moderately large degrees (p = 10, 20) without overgrowth of the computation time like in the p-version of the FEM [START_REF] Szabo | Finite element analysis[END_REF].

X 0 = [1, 1, 1] Singularity at X 0 = [2,
Second the accuracy with p = 1, 2 should be very poor, which should require a too large number of sub-domains. Hence a good compromise is to choose a degree in the range p ∈ [START_REF] Lucy | A numerical approach to the testing of the fission hypothesis[END_REF][START_REF] Belytschko | Meshless methods: an overview and recent developments[END_REF].

The TMM-computations are compared with FreeFem++ calculations using tetrahedral elements and quadratic shape functions. Linear interpolation is available with FreeFem++, but of course the convergence is lower. Elements with higher degrees are not available in 3D. One should mention that FreeFem++ is associated with several multifrontal solvers that are efficient for large-scale problems and we use the default sparse solver UMFPACK [START_REF] Davis | Algorithm 832: UMFPACK V4.3 -An unsymmetricpattern multifrontal method[END_REF].

First we concern the domain l x × l y × d = 10 × 10 × 1, discretized by 180000 elements and 255025 DOFs within FreeFem++ and by 12100 DOFs with TMM (degree p = 10, 100 "elements"), these two discretizations being designed to get an error less than 10 -3 . The approximate and exact solutions and the relative error on the plane z = 0.8 inside the domain obtained by TMM have been illustrated in Fig. 3.13, 3.14 and 3.15. The corresponding CPU-time are 6.39s with FreeFem++ and 5.08s with TMM. A larger domain (l x × l y × d = 25 × 25 × 1) has been considered (second test) by requiring the same accuracy. Several finite element meshes have been tested and the one which is able to reach this accuracy involves more than three millions DOFs (3213051), the TMM-mesh containing 625 "elements" and 75625 DOFs for a degree p = 10. The obtained CPU-times are 98.73s for FreeFem++ and 71.93s for TMM. Note that the condition number remains relatively small (in the range 10 9 ∼ 10 11 ) for all these large-scale TMM-tests, while it is much larger (about 10 32 ) within FreeFem++ due to a specific procedure to apply the boundary conditions.

In summary the present TMM procedure works well for problems re- quiring more than one million of unknowns by classical finite elements, even with hundreds of TMM-elements. Especially the method converges rapidly with the degree, even with a large number of "elements", which assesses the piecewise resolution techniques presented in section 3.2.3. Of course these first results deserve to be deepened, especially by a consistent use of parallel computations and by an optimal memory management.

Conclusion

A truly integration-free meshless method was assessed in this chapter for the purpose of application to large-scale problems involving three-dimensional domains. The specificity lies in the shape functions that are approximated solutions of the PDEs. So the discretization concerns only the boundary and it is accounted by least-square collocation. Then the control of the computation time and of the ill-conditioning of the matrix become central questions in view of practical applications of this method. It has been found that the accuracy of the method is governed by two antagonistic effects: first the classical convergence of Taylor series that leads to an exponential convergence with the degree, second the propagation of round-off errors, what grows with the degree and is facilitated by matrix illconditioning. According to our few numerical tests, the effect of round-off errors yields a lower bound of the error that seems more or less universal within the double-precision floating-point format. When the two effects are balanced, there is an accuracy reversal point beyond which the error increases. In practice, the method yields very accurate solutions up to a large degree (p = 60) and with many degrees of freedom: more than 70000, what would correspond to more than three millions of DOFs with a classical discretization method. This establishes the ability of the method to solve rather large-scale problems. Likely this result is due to the analytical solving inside the domain that limits the number of unknowns. In practice, the optimal degrees seem smaller (p = 10, 20) and the last test (section 3.5.3) shows that one can combine a rather large degree and many sub-domains. This allows to account for very large-scale problem, which is difficult to achieve within several other meshless techniques.

The computation time of the method has been discussed to locate the time consuming steps. In the present state of the algorithm, the largest part of the computation time is spent on calculating the value of the shape functions at the collocation points. The time distribution is quite different from usual methods (FEM, FDM) and this is an advantage in view of parallel computing. Nevertheless by comparing with a well established finite element code, one finds that TMM needs much less degrees of freedom and less computation time.

The numerical tests were based on 3D Laplace equation and linear isotropic elasticity. The application to non-linear problems has been sketched in previous works [START_REF] Zézé | Calcul de fonctions de forme de haut degré par une technique de perturbation (Computing high degree shape functions by a perturbation technique[END_REF][START_REF] Akpama | Amélioration de la "Taylor meshless method" et application à des structures élastiques en grandes transformations[END_REF][START_REF] Tampango | Développement d'une méthode sans maillage utilisant les approximations de Taylor (Development of a meshless method using Taylor approximations[END_REF] and a more complete study will be published elsewhere, the only difficulty being the calculation of shape functions for more generic PDEs.

Introduction

The non-linear Partial Differential Equations (PDEs) play an important role in science and engineering. To solve these non-linear PDEs, plenty of numerical methods have been developed over the past decades, such as finite difference [START_REF] Liszka | The finite difference method at arbitrary irregular grids and its application in applied mechanics[END_REF][START_REF] Chan | Generalized finite difference method for solving two-dimensional non-linear obstacle problems[END_REF], finite element, finite volume [START_REF] Shu | High-order finite difference and finite volume WENO schemes and discontinuous galerkin methods for CFD[END_REF][START_REF] Boyer | Finite volume method for 2D linear and nonlinear elliptic problems with discontinuities[END_REF] and spectral methods [START_REF] Tadmor | Convergence of spectral methods for nonlinear conservation laws[END_REF][START_REF] Tadmor | A review of numerical methods for nonlinear partial differential equations[END_REF][START_REF] Hussaini | Spectral methods in fluid dynamics[END_REF]. Meshless methods based on fundamental solution [START_REF] Balakrishnan | Osculatory interpolation in the method of fundamental solution for nonlinear Poisson problems[END_REF][START_REF] Tri | Perturbation technique and method of fundamental solution to solve nonlinear Poisson problems[END_REF][START_REF] Tri | High order continuation algorithm and meshless procedures to solve nonlinear Poisson problems[END_REF][START_REF] Gu | Singular boundary method for solving plane strain elastostatic problems[END_REF], radial basis function [START_REF] Kansa | Multiquadrics-a scattered data approximation scheme with applications to computational fluid-dynamics-I surface approximations and partial derivative estimates[END_REF][START_REF] Kansa | Multiquadrics-a scattered data approximation scheme with applications to computational fluid-dynamics-II solutions to parabolic, hyperbolic and elliptic partial differential equations[END_REF][START_REF] Chen | A meshless, integration-free, and boundary-only RBF technique[END_REF] and moving least-squares approximations [START_REF] Belytschko | Element-free Galerkin methods[END_REF][START_REF] Nayroles | Generalizing the finite element method: diffuse approximation and diffuse elements[END_REF][START_REF] Zhang | Least-squares collocation meshless method[END_REF][START_REF] Liu | Meshfree methods: moving beyond the finite element method[END_REF][START_REF] Timesli | A new algorithm based on moving least square method to simulate material mixing in friction stir welding[END_REF] were also applied to address non-linear PDEs. In these numerical methods, two main classes of discretization techniques were applied: Galerkin procedures and pointwise methods. Galerkin-based methods may work for large-scale problems, but their requirement for integration procedures could definitely lead to high computational costs. Pointwise-based techniques avoid the drawback of integration, but they could cause ill-conditioned matrices and numerical instabilities, which make it difficult to solve large-scale problems.

This paper proposes a new numerical method to solve non-linear elliptic PDEs. The original point of the method is to compute shape functions from Taylor series as proposed in [START_REF] Zézé | A boundary meshless method with shape functions computed from the PDE[END_REF]. Since the PDE is solved analytically, a significant reduction in the number of degrees of freedom (DOFs) has been obtained comparing with other discretization techniques, typically a reduction ratio in the range [START_REF] Belytschko | Meshless methods: an overview and recent developments[END_REF][START_REF] Poullikkas | The method of fundamental solutions for three-dimensional elastostatics problems[END_REF], see [START_REF] Tampango | Convergence analysis and detection of singularities within a boundary meshless method based on Taylor series[END_REF][START_REF] Yang | Solving large-scale problems by taylor meshless method[END_REF]. This strong reduction has two origins. First, because the PDE is solved analytically inside the domain, the discretization concerns only the boundary. Second, the Taylor series shows exponential convergence properties with respect to its degree and can characterize the unknown field in a large domain using a single series.

In practice, the method is confirmed to be efficient for polynomials with a large degree, generally in the range [START_REF] Belytschko | Element-free Galerkin methods[END_REF][START_REF] Kansa | Multiquadrics-a scattered data approximation scheme with applications to computational fluid-dynamics-II solutions to parabolic, hyperbolic and elliptic partial differential equations[END_REF], but the method also works for larger degrees. It is never easy to define an efficient integration procedure with such a class of polynomials. For this reason, a collocation-based discretization method was combined with Taylor series. As the pure collocation may show the numerical instabilities, a least-square minimization proposed in [START_REF] Fairweather | The method of fundamental solutions for elliptic boundary value problems[END_REF][START_REF] Zhang | Least-squares collocation meshless method[END_REF] was introduced to overcome this drawback. In this way, the method also enjoys many advantages of the collocation-based meshless methods: the discretization is fully defined from clouds of points, what's more, there is no required connection with the neighboring and no integration as well.

As a counterpart, collocation-based methods may be less reliable and lead to very ill-conditioned matrices, making them inapplicable to large-scale problems. This aspect has been well analyzed in the chapter 3. It has been found that the accuracy of the method is dominated by two antagonistic effects. The refinement of discretization results in exponential convergence, while the propagation of round-off errors governed by matrix ill-conditioning grows with the degree of polynomial shape functions. In practice, the proposed method yields very accurate solutions up to a large degree and with many DOFs.

In the previous papers [47, 49-51, 116, 117], it has been well established that the method always converges exponentially with the degree and can be applied to any boundary condition and any shape of domain. Bridging techniques based on Lagrange multipliers [START_REF] Tampango | Coupling of polynomial approximations with application to a boundary meshless method[END_REF] or least-square collocation [START_REF] Yang | Least-square collocation and Lagrange multipliers for Taylor meshless method[END_REF] have been presented to connect several sub-domains. It has been proven that the method also works well (p-convergence) in piecewise resolutions even with hundreds of sub-domains and more than 70, 000 DOFs, while the corresponding finite element numerical model requires over 3 millions DOFs.

The applications of TMM to non-linear PDEs have been attempted before [START_REF] Zézé | Calcul de fonctions de forme de haut degré par une technique de perturbation (Computing high degree shape functions by a perturbation technique[END_REF][START_REF] Tampango | Développement d'une méthode sans maillage utilisant les approximations de Taylor (Development of a meshless method using Taylor approximations[END_REF][START_REF] Akpama | Amélioration de la "Taylor meshless method" et application à des structures élastiques en grandes transformations[END_REF] but they are limited to simple or specific equations such as ODEs in [START_REF] Zézé | Calcul de fonctions de forme de haut degré par une technique de perturbation (Computing high degree shape functions by a perturbation technique[END_REF][START_REF] Tampango | Développement d'une méthode sans maillage utilisant les approximations de Taylor (Development of a meshless method using Taylor approximations[END_REF] and 2D-hyperelasticity in [START_REF] Akpama | Amélioration de la "Taylor meshless method" et application à des structures élastiques en grandes transformations[END_REF]. The purpose of this chapter is to define a general procedure for solving non-linear elliptic PDEs. Normally a numerical method for non-linear PDEs combines two steps: a linearization procedure and a discretization procedure. Some linearization techniques are available in literature such as Asymptotic Numerical Method (ANM) [START_REF] Damil | A new method to compute perturbed bifurcations: Application to the buckling of imperfect elastic structures[END_REF][START_REF] Yang | A new family of finite elements for wrinkling analysis of thin films on compliant substrates[END_REF], Newton method, Newton-Raphson method and Picard iteration [START_REF] Chen | Evaluation of thin plate spline based particular solutions for Helmholtz-type operators for the DRM[END_REF][START_REF] Golberg | An efficient mesh-free method for nonlinear reaction-diffusion equations[END_REF]. In this paper we use Newton method since it is the simplest and the most commonly used technique to achieve this linearization. The application of Newton method to a non-linear PDE leads to an inhomogeneous linear equation at each iteration. The TMM-solution of the latter linear PDE relies on the computation of Taylor series of a function composed by two analytic functions. This calculation is achieved by Automatic Differentiation (AD) of multivariate Taylor series [START_REF] Griewank | Evaluating higher derivative tensors by forward propagation of univariate Taylor series[END_REF]. Finally a general procedure will be presented to compute the polynomial shape functions for a large class of elliptic equations. This does not necessarily imply an efficient convergence of the non-linear algorithm nor does the convergence of the discretization. For instance it should be non-obvious to capture very localized patterns occurring in contact mechanics [START_REF] Johnson | One hundred years of hertz contact[END_REF] or in convection-dominated problems [START_REF] Liu | Solving singular convection-diffusion equation by exponentially-fitted trial functions and adjoint Trefftz test functions[END_REF].

The layout of this chapter is as the following. In section 4.2.1, the procedure of solving a PDE by the method of Taylor series is shortly recalled. section section 4.2.2 gives a brief description of Newton method. In section 4.2.3, a procedure based on AD is presented to compute the Taylor series of a function composed by two analytic functions. In section 4.2.4, one recalls the main features of TMM including its ability to solve large-scale problems. In section 4.3, a few one-dimensional and three-dimensional applications are discussed to assess the efficiency and robustness of the proposed technique.

Description of the method

The aim of this section is to present a general computational method based on Taylor series for solving non-linear elliptic PDEs. First we recall how to build the general solution of a linear PDE with variable coefficients, the latter being given by their Taylor series. Roughly, this corresponds to a method presented in previous papers ( [START_REF] Zézé | A boundary meshless method with shape functions computed from the PDE[END_REF][START_REF] Tampango | A new meshless method using Taylor series to solve elasticity problems[END_REF][START_REF] Tampango | Convergence analysis and detection of singularities within a boundary meshless method based on Taylor series[END_REF][START_REF] Tampango | Coupling of polynomial approximations with application to a boundary meshless method[END_REF]). In general, a non-linear equation is solved by a linearization procedure (Newton-Raphson, Newton, etc.) that transforms the non-linear systems into a sequence of linear problems. This leads to compose analytic functions by the method of Taylor series (section 4.2.2). Hence the key point is to establish a numerical technique giving the Taylor series of a function that is the composition of two analytic functions. Such techniques are available in literature concerning algorithmic differentiation (or automatic differentiation) [START_REF] Griewank | Evaluating higher derivative tensors by forward propagation of univariate Taylor series[END_REF]. The chosen algorithm is presented in section 4.2.3. The combination of Newton method and algorithmic differentiation provides a general procedure to solve elliptic PDEs from Taylor series. The main characteristics of TMM are the rapid convergence with the degree (p-convergence), the robustness of the algorithm with respect to several parameters and the ability to deal with rather large degrees. For the sake of completeness, these points are discussed briefly in section 4.2.4.

From PDE to Taylor series

Computing the shape functions

The aim of this section is to recall the procedures to build a general solution of a PDE. Here, we consider a linear PDE in Ω, a bounded domain in R 3 , with mixed boundary conditions as follows:

     -∆u(x) + g(x)u(x) = f (x) x ∈ Ω, u(x) = u d x ∈ Γ d , T u(x) = t n x ∈ Γ n . (4.1)
where the function u(x) is unknown while g(x) and f (x) are known, and T is a linear differential operator. The functions u(x), g(x) and f (x) are represented by truncated Taylor series with a development point

x c = [x c , y c , z c ]. For convenience, the vector x -x c is denoted by [ξ, η, γ]. u(x, y, z) = p k=0 k i=0 k-i j=0 u(i, j, k -i -j)ξ i η j γ k-i-j (4.2) g(x, y, z) = p-2 k=0 k i=0 k-i j=0 g(i, j, k -i -j)ξ i η j γ k-i-j (4.3) f (x, y, z) = p-2 k=0 k i=0 k-i j=0 f (i, j, k -i -j)ξ i η j γ k-i-j (4.4)
where u( , , ), g( , , ) and f ( , , ) represent the coefficients in the Taylor series of u(x), g(x) and f (x), respectively. Note that only the coefficients u( , , ) are unknown. The expansion of u(x) contains (p + 1)(p + 2)(p + 3)/6 coefficients which will be a large value when a large degree p is chosen.

The principle of TMM is to vanish the Taylor series of the residual R(x) = ∆u(x)g(x)u(x) + f (x) up to order p -2. The resulting equations then permit to reduce the number of unknown coefficients. In this respect, we shall distinguish (p -1)p(p + 1)/6 coefficients that will be expressed in terms of the (p + 1) 2 other coefficients. After this elimination, (p + 1) 2 independent unknowns will remain, which will be collected in a vector {v}.

First ∆u(x) is a polynomial of order p -2. The coefficient of the term containing ξ m η n γ s in ∆u(x) reads:

∆u(m, n, s) = (m + 2)(m + 1) u(m + 2, n, s) +(n + 2)(n + 1) u(m, n + 2, s) +(s + 2)(s + 1) u(m, n, s + 2). (4.5)
where the expression ∆u( , , ) represents the Taylor coefficients of ∆u(x).

Next we come to the second term g(x)u(x) in Eq. (4.1). The product of these two Taylor series leads to a function h(x). The coefficient h(m, n, s) of the term containing ξ m η n γ s in h(x) denotes:

h(m, n, s) = m i=0 n j=0 s k=0 g(i, j, k) u(m -i, n -j, s -k) (4.6)
Then the expression of the residual R(x) is as following:

R(x) = ∆u(m, n, s) -h(m, n, s) + f (m, n, s) ξ m η n γ s (4.7)
where ((m, n, s) ∈ [0, p -2], m + n + s ≤ p -2). The residual vanishes if the following conditions are satisfied:

∆u(m, n, s) -h(m, n, s) + f (m, n, s) = 0 (4.8)
There exist (p -1)p(p + 1)/6 linear equations in Eq. (4.8) containing (p + 1)(p + 2)(p + 3)/6 unknown coefficients u( , , ). We choose the coefficients u(0, , ) and u(1, , ) as independent unknowns. If the PDE (4.1) is consider as a differential equation with respect to x, its solution can be expressed in terms of the "initial values" u(0, y, z) and ∂u/∂x(0, y, z). Then the non-independent coefficients can be expressed in terms of the independent ones, which finally leads to a reduced expression u p containing (p + 1) 2 unknowns.

This procedure can be easily applied to a large degree of polynomials p and one can get the approximated solution as following:

u p =P 0 (x, y, z) + (p+1) 2 i=1 P i (x, y, z)v i =P 0 (x, y, z) + P {v}. (4.9)
where the first polynomial P 0 (x, y, z) balances the right-hand side f (x) of Eq. (4.1).

Boundary least-square collocation

Since the PDE in Eq. (4.1) has been approximately satisfied by the reduced approximation u p in Eq. (4.9), only boundary conditions need to be considered to determine the unknown vector {v}. Here, a collocation technique combined with the least-square method [START_REF] Fairweather | The method of fundamental solutions for elliptic boundary value problems[END_REF][START_REF] Zhang | Least-squares collocation meshless method[END_REF] is used. One chooses a set of nodes x i on Γ d and another set of nodes x j on Γ n . Then one minimizes the error between the approximate value u p , T u p and the given value of u d , t n at these points. It comes to minimize the following function:

T (v) = 1 2 x i ∈Γ d u p (x i ) -u d (x i ) 2 + w • 1 2 x j ∈Γn |T u p (x j ) -t n (x j )| 2 (4.10)
where the weight coefficient w is used to balance the two kinds of boundary conditions. This minimization leads to a linear system [K]{v} = {b} where [K] is an invertible matrix. Solving this system gives the vector {v} and therefore the approximate solution of the problem Eq. (4.1). Note the simple form of the matrix [K] and the vector {b} writes:

[K] αβ = M d i=1 P α (x i ) • P β (x i ) + w • Mn j=1 T P α (x j ) • T P β (x j ), α, β ∈ [1, (p + 1) 2 ] (4.11) {b} α = M d i=1 P α (x i )(u d (x i ) -P 0 (x i )) + w • Mn j=1 T P α (x j )(t n (x j ) -T P 0 (x j )) (4.12)
where M d and M n represent the number of collocation points on the displacement boundary Γ d and the Neumann boundary Γ n respectively.

It is impossible to solve any boundary value problem with a single Taylor series. Thus, numerical methods are introduced to connect several high-order polynomial approximations. Within this aim, several techniques have been presented and well investigated in the previous works: discrete Lagrange multipliers or discretized continuous Lagrange multipliers [START_REF] Tampango | Coupling of polynomial approximations with application to a boundary meshless method[END_REF][START_REF] Hu | A bridging technique to analyze the influence of boundary conditions on instability patterns[END_REF] and least-square collocation method [START_REF] Yang | Least-square collocation and Lagrange multipliers for Taylor meshless method[END_REF]. In this paper, we focus on the least-square collocation for the interface conditions in piecewise resolution which will be briefly recalled in the following. Here we consider a domain as shown in Fig. 4.1. The independent highorder general solution for each sub-domain, u 1 p and u 2 p , is deduced from Taylor series respectively as stated in section 4.2.1.1. The consideration of the boundary conditions for each sub-domain yields two functions, T 1 (v 1 ) and T 2 (v 2 ), as described in section 4.2.1.2. Note that the number of unknowns for each sub-domain is (p + 1) 2 . The constraints on the interface are as follows:

ψ 1 j (v 1 , v 2 ) = u 1 p (x j ) -u 2 p (x j ) = 0 ψ 2 j (v 1 , v 2 ) = ∂u 1 p ∂n (x j ) - ∂u 2 p ∂n (x j ) = 0 ∀x j ∈ Γ r (4.13)
where x j represents the collocation points on the interface Γ r . By applying the least-square collocation method, the transmission conditions Eq. (4.13) are satisfied in a mean square sense via the following coupling function:

C ls (v 1 , v 2 ) = 1 2 Mr j=1 u 1 p (x j ) -u 2 p (x j ) 2 + w • ∂u 1 p ∂n (x j ) - ∂u 2 p ∂n (x j ) 2 (4.14)
where the weight coefficient w is used to balance the two kinds of constraints on the interface. For the whole problem, searching the variables {v 1 } and {v 2 } is equivalent to minimizing the following function:

T (v 1 , v 2 ) = T 1 (v 1 ) + T 2 (v 2 ) + C ls (v 1 , v 2 ) (4.15)

Newton Method

As the method recalled in section 4.2.1 relies on a superposition principle (see Eq. (4.9)), it can be only applied to a linear equation. In this respect, a linearization procedure is required if one wants to solve a non-linear equation. In this paper, we limit ourselves to Newton method that is the simplest way to achieve this linearization. To show this linearization procedure in details, we consider the following non-linear equation:

     -∆u(x) + h(u(x)) = q(x), x ∈ Ω, u(x) = u d x ∈ Γ d , T u(x) = t n x ∈ Γ n . (4.16)
where Ω is a bounded domain and h(•) represents a non-linear function.

The initial solution and the correction of the solution are represented by u 0 and δu respectively. Usually the initial solution is set to be u 0 = 0. The (i + 1) th Newton iteration is as follows:

           -∆(δu) + h (u i )δu + R(u i ) = 0, x ∈ Ω, u i+1 (x) = u d x ∈ Γ d , T u i+1 (x) = t n x ∈ Γ n , u i+1 = u i + δu. (4.17)
where h (•) represents the derivative of function h(u) with respect to u and R(u i ) represents the residual of the solution obtained at the previous iteration.

One stops the iteration when one of the following conditions is satisfied:

• |log 10 R n -log 10 R n-1 | ≤ 10 -2 ,
• The maximum number of iterations is equal to 20, where R n represents the maximum value of the function R(u n ).

At each iteration of the algorithm, one gets a linear PDE which will be solved by TMM. According to the computing procedure of TMM in section 4.2.1, one needs to expand the known functions h (u i ) and R(u i ) into Taylor series. This can be realized by AD and will be described in the next section.

Automatic Differentiation

At each step of the non-linear algorithm, one needs to solve a linear PDE with variable coefficients as Eq. (4.1), where the two functions g(x) and f (x), x ∈ Ω, depend on the approximation of the unknown function u(x) at the previous iteration:

g(x) = h (u i (x)) f (x) = ∆(u i (x)) -h(u i (x)) + q(x) (4.18)
The computation of TMM-shape functions P (x) requires the derivatives of composed functions G(x) = F (U (x)), where F (•) is a given analytic function and U (x) is another analytic function approximated by a truncated Taylor series. The numerical composition of a Taylor series with a given analytic function can be achieved with the standard tools of AD.

First let us recall the statement of AD for univariate functions. The aim is to compute the numerical values of Taylor coefficients (equivalent of the derivatives) of G(x) at a point x c . There are well known recurrence formulae if the function F (•) is simple: product, fraction, exponential, etc. [START_REF] Griewank | Evaluating derivatives: principles and techniques of algorithmic differentiation[END_REF]. These elementary functions are included in a small library. As illustrated in Table 4.1, such a recurrence can be extended to more intricate functions by combining these basic operations. AD works for any function that is the composition of elementary functions, which is not restrictive in practice. This can be implemented within Matlab by building a class for Taylor series, including the type and the overloaded operations (i.e. the rules in right column of Table 4.1). Table 4.1: Successive elementary operations to compute F (U (x)) = U 2 /(1 + e U ) and its derivatives.

1 V = U • U V n = n k=0 U k • U n-k 2 W = e U W n = 1 n n-1 k=0 (n -k)W k • U n-k 3 X = 1 X 0 = 1, X n = 0 (n ≥ 1) 4 Y = X + W Y n = X n + W n 5 G = V /Y G n = (V n - n-1 k=0 G k • Y n-k )/Y 0
The extension to multivariate Taylor series has been addressed since a long time [START_REF] Griewank | Evaluating higher derivative tensors by forward propagation of univariate Taylor series[END_REF][START_REF] Bischof | Structured second-and higher-order derivatives through univariate Taylor series[END_REF][START_REF] Neidinger | Directions for computing truncated multivariate Taylor series[END_REF]. The best way seems to compute Taylor series of univariate functions t → G(x c + tr) in a sufficiently large number of directions r and to deduce all the derivatives of the multivariate function G(x) up to a given order from all these univariate Taylor series. We refer to the AD literature for the details of the algorithm, including [START_REF] Neidinger | Efficient recurrence relations for univariate and multivariate Taylor series coefficients[END_REF].

Neidinger has done an open source Matlab implementation of this theory in his work on "Directions for computing truncated multivariate Taylor series" [START_REF] Neidinger | Directions for computing truncated multivariate Taylor series[END_REF]. The code sources can be found on his web page [START_REF] Neidinger | Directions for computing truncated multivariate Taylor series: MATLAB implementation[END_REF]. However for our needs in the TMM, we propose an implementation of the truncated power series algebra (TPSA) based on the works of [START_REF] Kaltchev | Building truncated Taylor maps with mathematica and applications to FFAG[END_REF][START_REF] Kaltchev | Implementation of TPSA in the mathematica code LieMath[END_REF][START_REF] Kaltchev | Mathematica program for extracting one-turn lie generator map. Application of TPSA[END_REF]. The Mathematica codes of Kaltchev have been translated and adapted in Matlab for the automatic computation of the truncated multivariate Taylor series of any Matlab functions. The TMM Matlab Toolbox contains two classes: "pindexes" and "ndseries". The class "pindexes" is a singleton class and contains all the "boxes" described in [START_REF] Kaltchev | Building truncated Taylor maps with mathematica and applications to FFAG[END_REF][START_REF] Kaltchev | Implementation of TPSA in the mathematica code LieMath[END_REF][START_REF] Kaltchev | Mathematica program for extracting one-turn lie generator map. Application of TPSA[END_REF]. Once the space dimension of the problem (i.e. three dimensional problems) and the order of the truncated series are chosen, an instance of "pindexes" is created with all the data to manage the series computation. The class "ndseries" is the main class where the overloading of the operators and elementary functions have been implemented. An instance of "ndseries" has a member called "coeffs" (1D Matlab array) which stores the multivariate Taylor coefficients of the object. This member is manipulated by the operators and elementary functions using the "boxes" defined in the singleton object "pindexes". The following operators and functions are currently overloaded in the TMM Toolbox: acos, asin,

. Thanks to this Matlab toolbox, the Taylor series of the functions (4.18) can always be computed.

Recalling the basic properties of TMM

In this section, two 3D linear problems are considered to briefly recall the basic properties of TMM, especially exponential convergence, robustness, efficiency and ability to solve large-scale problems. The first case is the Laplace equation in a sphere domain with singular solution and this case will be solved by a single Taylor series. The next one is issued from the study of instabilities of film-substrate systems [START_REF] Chen | Herringbone buckling patterns of compressed thin films on compliant substrates[END_REF][START_REF] Xu | 3D finite element modeling for instabilities in thin films on soft substrates[END_REF][START_REF] Yu | A novel two-dimensional finite element to study the instability phenomena of sandwich plates[END_REF], where one has to describe many spatial oscillations. This second case will be solved with hundreds of sub-domains.

Laplace equation with singular solution

The problem considered here reads:

   ∆u(x, y, z) = 0 in Ω, u d (x, y, z) = 1 r on ∂Ω. (4.19)
where Ω = {x, y, z|x 2 + y 2 + z 2 ≤ 1} and r = xx 0 . The exact solution is as following with a singularity at x 0 = [x 0 , y 0 , z 0 ]:

u ex (x, y, z) = 1 x -x 0 (4.20)
Throughout this chapter, the technique will be evaluated through the maximum relative error which is defined by:

E = max |u p (x) -u ex (x)| max |u ex (x)| (4.21)
The numerical solution obtained by TMM depends on two parameters: the degree p of the polynomials and the number of collocation points M . The influence of the number of collocation points is illustrated in Fig. 4.2 for three values of the degree p = 10, 20 and 30 with a singularity point at

X 0 = [1, 1, 1]
. The maximum error decreases with M until an optimal number M = 2(p + 1) 2 where it becomes stable. The divergence for small value of M is easily understood by considering the equivalent 2D problem: in this case, the minimization of Eq. (4.10) is simply an identification of Taylor series that cannot work for a small M . For larger value of M , it is observed that the identified solution does not depend on the collocation points and the remaining error is only due to the truncation of Taylor series. The convergence with the polynomial degree p is presented in Fig. 4.3 for different singular points X 0 = [1, 1, 1] and X 0 = [2, 2, 2] with 2(p + 2) 2 collocation points to ensure the convergence. The results illustrated in Fig. 4.3 show that this technique converges exponentially with the degree and the rate of convergence depends on each specific problem. This 3D example illustrates the rapid convergence with the degree, as well as the robustness of the method with respect to the number of discretization points and to the degree. For more assessments, we refer to [START_REF] Tampango | Coupling of polynomial approximations with application to a boundary meshless method[END_REF][START_REF] Yang | Solving large-scale problems by taylor meshless method[END_REF].

A large-scale test

This benchmark is motivated by wrinkling instabilities of film-substrate systems [START_REF] Chen | Herringbone buckling patterns of compressed thin films on compliant substrates[END_REF][START_REF] Xu | 3D finite element modeling for instabilities in thin films on soft substrates[END_REF][START_REF] Yu | A novel two-dimensional finite element to study the instability phenomena of sandwich plates[END_REF], where one observes many spatial oscillations. Here one accounts only for the substrate, and the film is represented by a prescribed sinusoidal displacement. The domain is a rectangular hexahedron with sides x × y ×d = 10×10×1. We split this domain into cubes of size 1×1×1 that will be called "elements" for simplicity. This results in one element through the thickness and two elements along each period. In this case, Dirichlet boundary conditions are applied and the exact solution is u(x, y, z) = sin(πx)sin(πy)sinh( √ 2πz)/sinh( √ 2πd). On each surface of the cube "element", p 2 collocation points for the boundary conditions and p 2 /4 collocation points for the transmission conditions are chosen to ensure the convergence. The results presented in Table 4.2 show that the error decreases rapidly with the degree (p-convergence). This demonstrates the robustness and rapid convergence of the piecewise resolution even with a large number of subdomains. A comparison between the TMM and FEM has been done for this example. The efficient open source code FreeFem++ [START_REF] Hecht | New development in FreeFem++[END_REF] is used to perform the finite element simulation. It contains rapid multi-frontal linear solvers such as UMFPACK [START_REF] Davis | Algorithm 832: UMFPACK V4.3 -An unsymmetricpattern multifrontal method[END_REF]. The domain is discretized by tetrahedral elements with quadratic shape functions (elements T2). To obtain a relative error less than 10 -3.2 , the finite element mesh requires 417501 degrees of freedom. The corresponding CPU-time is 10.31s with FreeFem++ and 5.08s with TMM. We have also implemented the standard Finite Difference Method (FDM) and tested various regular meshes, especially with steps h = 0.0625, h = 0.05 and h = 0.03875. The FDM converges rather slowly, with respective errors of 10 -2.64 , 10 -2.86 and 10 -3.21 . The last case requires a huge computational time because of the large number of DOFs (2338588). Finally our most significant FDM result corresponds to a step h = 0.05, 752419 DOFs, an error 10 -2.86 , and a CPU-time of 77.3s, as compared with FEM (417501 DOFs, 10.31s) and TMM (12100 DOFs, 5.08s). Remark that such comparisons are never perfectly fair because they can be computer-dependent, implementation-dependent and even author-dependent, FreeFem++ benefits from an optimized implementation and from elements T2 (with linear interpolation elements (T1), it converged also more slowly). A larger domain ( x × y × d = 20 × 20 × 1) has been considered (second test) by requiring the same accuracy. Several finite element meshes have been tested and the one which is able to reach this accuracy involves more than one million DOFs (1656441), the TMM-mesh containing 400 "elements" and 48400 DOFs for a degree p = 10. The obtained CPU-times are 46.44s for FreeFem++ and 25.67s for TMM. The results about CPU-time are obtained by using a personal computer with inter(R) core(TM) i5-2300 type CPU and 8 GB computer memory.

     u(x, y, z) = 0 in Ω u(x, y, 0) = 0 u(x,

Comments for TMM

In the previous papers [START_REF] Zézé | A boundary meshless method with shape functions computed from the PDE[END_REF][START_REF] Tampango | A new meshless method using Taylor series to solve elasticity problems[END_REF][START_REF] Tampango | Convergence analysis and detection of singularities within a boundary meshless method based on Taylor series[END_REF][START_REF] Tampango | Coupling of polynomial approximations with application to a boundary meshless method[END_REF], it has been validated that TMM is robust and efficient for solving linear PDEs, even with large numbers of DOFs or hundreds of sub-domains. One can also use TMM to handle problems with any shape of domain and any boundary condition. For instance, TMM works well for a problem with an amoeba-like boundary shape in [START_REF] Zézé | A boundary meshless method with shape functions computed from the PDE[END_REF] and for elasticity problems with displacement and stress boundary conditions in [START_REF] Tampango | A new meshless method using Taylor series to solve elasticity problems[END_REF]. In practice, TMM needs less DOFs as compared with other discretization techniques since the PDE is solved analytically. This small number of unknowns is a typical feature of TMM. One observes a significant reduction of the number of DOFs with a ratio in the range [START_REF] Belytschko | Meshless methods: an overview and recent developments[END_REF][START_REF] Poullikkas | The method of fundamental solutions for three-dimensional elastostatics problems[END_REF] with respect to the FEM.

Numerical applications

In this section, several non-linear problems in one-dimension and three-dimension are considered to validate the robustness and convergence of the proposed technique which combines Taylor meshless method with algorithmic differentiation. The proposed method for non-linear PDEs includes two main steps: firstly, one linearizes the non-linear PDE by Newton method and gets a linear PDE at each iteration as in Eq. (4.17); secondly, this linear PDE is solved by combining Taylor series inside the domain and least-square collocation at boundaries and interfaces. The TMM procedure for one-dimensional linearized ODEs is recalled in Appendix B and for three-dimensional cases it has been presented in section 4.2.1, see also [START_REF] Yang | Solving large-scale problems by taylor meshless method[END_REF] for more details.

One-dimensional non-linear problems

An one-dimensional example

Here, we consider a non-linear ODE as following:

-u (x) + u 3 (x) -1 = 0, x ∈ [-1, 1] u(±1) = 0. (4.23)
One should note that only 2 DOFs and 2 points are necessary in this considered problem. The maximum residual of the ODE is considered as an indicator for the measurement of the accuracy of the approximation. The results will be discussed under two aspects: first the convergence of the Newton iterations for a fixed degree of polynomials p; second the convergence with the refinement of the discretization (p increases).

The convergence with the number of iterations and the degree of polynomials is presented in Table 4.3. For each fixed degree of polynomials p, the maximum residual decreases rapidly with the number of iterations until a maximum accuracy is reached and then it remains constant. Only 3 iterations are necessary to get the maximum accuracy for p = 10 and 20, and 4 iterations are sufficient for p = 30 and 40.

Hence with a sufficient number of Newton iterations, the proposed technique converges exponentially with the degree of polynomials p and up to a very high accuracy.

Bratu equation

The classical one-dimensional Bratu's problem is considered as following:

u (x) + λe u(x) = 0, x ∈ [0, 1] u(0) = u(1) = 0. (4.24)
where λ is a constant. The analytical solution of Bratu's problem is in the following form:

u(x) = -2ln cosh((x -1 2 ) θ 2 ) cosh( θ 4 ) (4.25)
where θ is the solution of θ = √ 2λcosh( θ 4 ). The accuracy of the discrete problem depends on two parameters: the number of iterations and the degree of polynomial shape functions p. Here we focus on the influence of the degree of polynomials on the accuracy. Three cases are considered here with different values of λ: λ = 1, 1.5 and 2 respectively. The numerical solutions of these cases obtained by TMM with p = 20 are shown in Fig. 4.4. It appears that the approximated solutions agree well with the analytical solutions. The convergence with the degree of polynomials (p ∈ [START_REF] Lancaster | Surfaces generated by moving least squares methods[END_REF][START_REF] Mukherjee | The boundary node method for potential problems[END_REF]) is summarized in Table 4. 4. In all these cases, the results show that the maximum relative error decreases rapidly when increasing the degree of polynomials, as in the linear cases. The relationship between the maximum residual and the maximum relative error of the approximated solution has been illustrated in Fig. 4.5. It can be seen that the maximum relative error decreases with the maximum residual when increasing the degree of polynomials (p ∈ [START_REF] Lancaster | Surfaces generated by moving least squares methods[END_REF][START_REF] Mukherjee | The boundary node method for potential problems[END_REF]). This means that the maximum residual can be a good indicator of accuracy of the approximated solution when there is no reference solution of the non-linear problem.

Three-dimensional non-linear problems

In the previous section, we have tested the efficiency and robustness of the proposed technique for non-linear ODEs. Since our final goal is to apply the TMM associated with Newton method for large-scale problems, several three-dimensional non-linear PDEs are discussed in this section.

Problem one

We consider the following non-linear problem in a unit sphere:

-∆u(x, y, z) + u 3 (x, y, z) = f (x, y, z) in Ω, u(x, y, z) = sin(x) sin(y)e -z on ∂Ω. (4.26) where the right-hand function f (x, y, z) denotes:

f (x, y, z) = sin(x) sin(y)e -z + sin 3 (x) sin 3 (y)e -3z

The exact solution is u ex = sin(x) sin(y)e -z . The accuracy of the solution depends on three factors: the number of collocation points on the boundary; the degree of polynomials and the number of Newton iterations, among which the influence of the number of collocation points for three-dimensional linear PDEs has been well investigated in Yang et al. [START_REF] Yang | Solving large-scale problems by taylor meshless method[END_REF]. This has been also briefly recalled in section 4.2.4. The maximum accuracy is obtained for about 2(p + 1) 2 collocation points and it remains constant beyond this threshold. Here 2(p + 2) 2 collocation points have been chosen to ensure the convergence. The convergence with Newton iterations (within 10 iterations) and the discretization of TMM (p ∈ [5,[START_REF] Kansa | Multiquadrics-a scattered data approximation scheme with applications to computational fluid-dynamics-II solutions to parabolic, hyperbolic and elliptic partial differential equations[END_REF]) are presented in Fig. 4.6. For each degree p, Newton iterations converge rapidly until obtaining a maximum accuracy and then it remains constant beyond this threshold. The number of iterations necessary for the convergence is rather small. For instance, only 3 Newton iterations are sufficient to get the maximum accuracy for p = 20 and 30.

The convergence with the discretization can be also seen in Fig. 4.6 and detailed information is presented in Table 4.5. One should note that the number of DOFs is (p + 1) 2 . The maximum relative error and the maximum residual decrease exponentially with the degree of polynomials. From a degree p = 20, the accuracy is no longer improved because one is closed to the unit round-off error within the double precision floating-point format. where r = (xx 0 ) 2 + (yy 0 ) 2 + (zz 0 ) 2 . The exact solution is as follows:

u ex (x, y, z) = 1 (x -x 0 ) 2 + (y -y 0 ) 2 + (z -z 0 ) 2 (4.28)
This solution has a singularity at X 0 = [x 0 , y 0 , z 0 ]. The critical number of collocation points on the boundary for the 3D case with p order polynomial shape functions is about 2(p + 1) 2 and the accuracy remains constant beyond this threshold. Here 2(p + 2) 2 collocation points have been chosen to ensure the convergence. The p-convergence has been reported in Table 4.6 for the problem Eq. (4.27) with the nearest singularity point X 0 = [1.0, 1.0, 1.0]. The maximum relative error and maximum residual decrease exponentially with the degree in several iterations. The influence of the singular point on the rate of the convergence has been shown in Fig. 4.7. Two different singular points X 0 = [1, 1, 1] and X 0 = [2, 2, 2] have been chosen. The results show that the quality of convergence depends on the considered problem, especially on the position of the singular point X 0 . The rate of convergence is better if the singularity is far from the domain which means that TMM is easier to solve problem with a flat solution. In the case with singularity at X 0 = [2, 2, 2], from a degree p = 30, the accuracy is no longer improved because one is closed to the unit round-off error within the double precision floating-point format. The CPU-time presented in Table 4.6 increases rapidly with the degree of polynomials p, the CPU-time is 0.09s for p = 10 while 434.47s for p = 50. As analyzed in a previous paper [START_REF] Yang | Solving large-scale problems by taylor meshless method[END_REF], the computation of shape functions is the most time-consuming part in TMM and it relates to the number of unknowns (p + 1) 2 and the number of collocation points 2(p + 2) 2 as shown in Eq. (5.9). Hence the number of quantities to be computed P α (x j ) is of the order of p 4 . Nevertheless it is not necessary to choose a so large degree and in practice the optimal degree seems in the range p ∈ [START_REF] Lancaster | Surfaces generated by moving least squares methods[END_REF][START_REF] Belytschko | Meshless methods: an overview and recent developments[END_REF].

Problem three

We consider the third non-linear problem in a spherical domain (Ω = {x, y, z|x 2 + y 2 + z 2 1}):

-∆u(x, y, z) + u 3 (x, y, z) = 1 in Ω, u(x, y, z) = 0 on ∂Ω. (4.29)
The number of DOFs of TMM with any degree of polynomials p is (p + 1) 2 and the critical number of collocation points is about 2(p + 1) 2 . In the following 2(p + 2) 2 collocation points will be chosen to ensure the convergence. The maximum residual of the PDE is considered here as an indicator for the measurement of the accuracy of the approximation. The convergence with the refinement of discretization is presented in Table 4.7. Only 2 iterations are necessary to get the maximum accuracy for p = 10 and 3 iterations for The exact solution is u(x, y, z) = e x sin(πy/2) sin(πz/2). The critical number of collocation points on the boundary is about 2(p + 1) 2 . Here 2(p + 2) 2 collocation points have been chosen to ensure the convergence. Then the accuracy of the approximated solution depends on two factors: the discretization of TMM (the degree of Taylor series p) and the number of Newton iterations. The convergence with the previous two factors has been illustrated in Table 4.8. The number of Newton iterations necessary to ensure the convergence is still small. Only 4 iterations are needed to obtain the best accuracy in the case with p ≤ 30. With a sufficient number of Newton iterations, the maximum relative error and the maximum residual decrease exponentially with the degree of polynomial shape functions.

A large-scale test

One considers a non-linear PDE with Dirichlet boundary conditions, the linear version being discussed in section 4.2.4.2:

     -∆u(x, y, z) + u 3 (x, y, z) = f 3 (x, y, z)
in Ω u(x, y, 0) = 0 u(x, y, d) = sin(πx) sin(πy) We split the domain into 100 cubes of size 1 × 1 × 1 that will be called "elements" for simplicity, see Fig. 4.9. This corresponds to one element through the thickness and two elements along each period. The Dirichlet boundary conditions and connections between sub-domains are achieved by the leastsquare collocation technique. On each surface of the cube "element", p 2 collocation points for the boundary conditions and p 2 /4 collocation points for the transmission conditions are used to ensure the convergence.

The p-convergence of the proposed technique for the problem Eq. (4.32) has been shown in Table 4.9. The error decreases exponentially with the degree when taking enough Newton iterations. The number of iterations necessary for convergence is still small, only 2 Newton iterations for a degree p ≤ 10, while 3 iterations are enough for a degree p in the range [START_REF] Mclain | Drawing contours from arbitrary data points[END_REF][START_REF] Belytschko | Meshless methods: an overview and recent developments[END_REF]. This demonstrates the robustness of TMM associated with Newton method for solving non-linear elliptic PDEs, even with large numbers of TMM-elements and large numbers of DOFs. The CPU-time for this case Eq. (4.32) with 100 elements, p = 10 and 2 Newton iterations, is 18.25s. The CPU-time for a similar linear problem Eq. (4.22) recalled in section 4.2.4.2, with the same domain and boundary conditions, is 5.08s. The CPU-time of the non-linear case for each iteration (about 9s) is a little larger than that of the linear case (5.08s). Likely this extra time is due to the calculation of the product of h (u i )δu (in Eq. (4.17)) by AD. This point deserves further discussion, but in the present state, the automatization has already a moderate counterpart in CPU-time.

Conclusion

A general numerical method associating Taylor meshless method with Automatic Differentiation, called TMM, is proposed in this chapter to solve a class of non-linear elliptic PDEs. First Newton method is used to linearize the non-linear PDE, but the present method could be easily adapted to any linearization technique. Then the obtained linear problem is solved by TMM. The specificity of TMM is the numerical computation of a family of shape functions that are solutions of the PDE in the sense of Taylor series. Solving the linearized equations needs to compute the Taylor series of the functions composed by two analytic functions. This key point is realized by using the Automatic Differentiation.

The numerical tests demonstrate that non-linear TMM keeps the same advantages as in the linear case. First the discretization technique converges exponentially with the degree and it requires much less degrees of freedom than other techniques like the finite element method. Next it is quite robust with respect to the number of collocation points and to the degree of Taylor series that can be very large. In the considered numerical example, the Newton method converges rapidly with one or two additional iterations for large degrees. TMM does not require a mesh, but merely requires a cloud of collocation points on the boundary. The computation time has been estimated in a large-scale test (12, 100 DOFs for TMM, but more than 400,000 for the equivalent FEM mesh) and it remains of the same order as that required in the linear case, which demonstrates that the implementation via Automatic Differentiation does not alter the effectiveness of the method. The proposed technique also works well for a very large degree, for instance p = 50, but leads to high computational costs of the order O(p 4 ). Nevertheless it's not necessary to choose such a large degree, in practice an optimal degree seems smaller and usually in the range [START_REF] Lancaster | Surfaces generated by moving least squares methods[END_REF][START_REF] Kansa | Multiquadrics-a scattered data approximation scheme with applications to computational fluid-dynamics-II solutions to parabolic, hyperbolic and elliptic partial differential equations[END_REF]. The numerical tests in this chapter are limited to a single non-linear elliptic PDE. Some applications to practical engineering problems in the fields of solid and fluid mechanics are in progress. Accounting for singularities is very important in solid mechanics because of cracks, corners or concentrated forces, and an ongoing research will include singular shape functions in the present computational framework.

Introduction

This chapter deals with the numerical solution of elliptic Partial Differential Equations (PDEs) by using Taylor series. The main interest of a Taylor series is the property of exponential convergence with the degree, which leads to strong reductions of the number of degrees of freedom (DOFs) as compared, for instance, with the finite element method [START_REF] Tampango | Convergence analysis and detection of singularities within a boundary meshless method based on Taylor series[END_REF][START_REF] Yang | Solving large-scale problems by taylor meshless method[END_REF]. Nevertheless the radius of convergence of a Taylor series is limited by the nearest singularity and singularities are common in engineering mechanics because of corners, cracks or pointwise forces. So it is expected that such singularities will slow the convergence. In this chapter, the relation between singular solutions and convergence of Taylor series is addressed.

Polynomial solutions of some PDEs have been used since a long time to build large finite elements, for instance in plane elasticity [START_REF] Jirousek | Large finite elements method for the solution of problems in the theory of elasticity[END_REF] or for plate bending [START_REF] Jirousek | Hybrid-Trefftz plate bending elements with p-method capabilities[END_REF], even with relatively large degrees. When associated with a specific treatment of boundary and interface conditions, this method has been called hybrid-Trefftz method, but there are also many other ways to account for these boundary-interfaces conditions [START_REF] Kita | Trefftz method: an overview[END_REF]. Trefftz methods reduce dramatically the number of DOFs because only the boundary and interfaces have to be discretized, but it does not work as well for non-homogeneous and nonlinear problems, in which cases exact solutions are not known. To extend the method in the non-linear range, it has been proposed to compute approximate solutions in the sense of Taylor series [START_REF] Zézé | A boundary meshless method with shape functions computed from the PDE[END_REF] and the latter procedure is efficient for a large class of PDEs, even in non-linear cases [START_REF] Yang | Taylor meshless method for solving non-linear partial differential equations[END_REF], see chapter 3 and 4. When coupled with appropriate collocation-based subdomain methods [START_REF] Tampango | Coupling of polynomial approximations with application to a boundary meshless method[END_REF][START_REF] Yang | Solving large-scale problems by taylor meshless method[END_REF], it leads always to exponential convergence results. In this form, it has been called Taylor Meshless Method (TMM). These excellent convergence properties are illustrated in Fig. 5.1, where three problems are studied, one with a very flat response and two others with larger gradients. In the flat case, the convergence is very rapid and a small degree is sufficient to get a high precision ( 10 -14 ) that is close to the limit for a real number in double-precision floating-point format. In the two other cases, there are two types of response. For a small degree, the response is governed by the exponential convergence of Taylor series, even for large domains. For larger degrees, the accuracy deteriorates slightly due to the propagation of round-off errors, what is favored by matrix ill-conditioning. The accuracy reversal observed in Fig. 5.1 is obtained for degrees p = 25 and p = 45, which means that TMM works well with large degrees. How to keep these excellent convergence properties when the solution is not analytic?

The treatment of singular solutions is a key issue in computational fracture mechanics. The bibliography on the topics is huge and we just refer for One tries to recover fundamental solutions within 3D linear elasticity. Three cases were considered, from [START_REF] Yang | Solving large-scale problems by taylor meshless method[END_REF].

instance to [START_REF] Liebowitz | Finite element methods in fracture mechanics[END_REF] or [START_REF] Belytschko | A review of extended/generalized finite element methods for material modeling[END_REF]. Standard finite elements, which ignore the singular fields, could be used to analyse elastic cracked bodies, but this "should require extremely fine meshes" [START_REF] Liebowitz | Finite element methods in fracture mechanics[END_REF], so that various specific elements were introduced to better capture the singular behavior near crack tip. The first one [START_REF] Wilson | Combined mode fracture mechanics[END_REF] simply consists in including the analytic singular functions and minimizing with respect to the stress intensity factors, but there are many variants in the literature of the seventies, see for instance [START_REF] Tracey | Finite elements for determination of crack tip elastic stress intensity factors[END_REF] or [START_REF] Benzley | Representation of singularities with isoparametric finite elements[END_REF]. The modern approaches, called extended finite element method or generalized finite element method, are based on the partition of unity [START_REF] Melenk | The partition of unity finite element method: basic theory and applications[END_REF] and can include analytic singular fields and discontinuities of the displacement [START_REF] Moës | A finite element method for crack growth without remeshing[END_REF]. Note that some problems with singular solutions were solved in a meshless framework by using a full basis of singular solutions [START_REF] Li | Combined methods for elliptic equations with singularities, interfaces and infinities[END_REF][START_REF] Li | Trefftz, collocation, and other boundary methods: a comparison[END_REF]. The convergence properties of Taylor Meshless Method in the presence of singularities are adressed in this chapter, especially for domains with corners and for classical crack problems. As expected, it is not possible to recover exponential convergence up to a high accuracy, in the same way as in Fig. 5.1. Building on the significant achievements within fracture mechanics, we shall combine the high degree polynomials computed by Taylor series with analytically known singular solutions. Indeed the ability to solve problems with singular solutions is an important issue in the development of Taylor-based numerical methods.

Combining Taylor series and singular solution

The issue of the present chapter is the behavior of Taylor methods for elliptic systems having a singular solution. In this part 5.2, it will be shown that the convergence of Taylor meshless method, as illustrated in Fig. 5.1, deteriorates for corner-domains (section 5.2.3). Next a variant of Taylor Meshless Method is proposed simply by including few singular solutions in the basis of shape functions (section 5.2.4). For completeness, two main features of TMM are recalled: first the algorithm to compute polynomial approximate solutions of the PDE by the method of Taylor series (section 5.2.1), next a procedure to build the global problem by a boundary least-square collocation method (section 5.2.2).

Compute shape functions from Taylor series

The Taylor Meshless Method can be split in two parts. First the differential equation is solved in a quasi-exact manner by using the concept of Taylor series, next the boundary and interface conditions are accounted by a collocation-based technique. Let us begin by the solving inside the domain and describe the method for a simple two-dimensional linear PDE with mixed boundary conditions:

     -∆u(x) + g(x)u(x) = f (x) x ∈ Ω, u(x) = u d x ∈ Γ d , Tu(x) = t n x ∈ Γ n .
(5.1)

where the function u(x) is unknown while g(x) and f (x) are given analytical functions, and T is a linear differential operator. Firstly, one expands u(x), g(x) and f (x) by using truncated Taylor series with a development point

x c = [x c , y c ]. For convenience, the vector x -x c is denoted by [ξ, η]. u(x, y) = p j=0 j i=0 u(i, j -i)ξ i η j-i (5.2) g(x, y) = p-2 j=0 j i=0 g(i, j -i)ξ i η j-i (5.3) f (x, y) = p-2 j=0 j i=0 f (i, j -i)ξ i η j-i (5.4)
where u( , ), g( , ) and f ( , ) represent the coefficients of the Taylor series of u(x), g(x) and f (x) respectively. Note that only the coefficients u( , ) are unknown.

The principle of TMM is to vanish the residual R(x) of the PDE up to the order (p -2). Substituting Eqs. (5.2-5.4) into the PDE in Eq. ( 5.1), one can express R(x) in the form of Taylor series:

R(x) =∆u(x) -g(x)u(x) + f (x) ≈ p-2 m=0 m n=0 R(m, n)ξ m η n (5.5)
where R(m, n) represents the coefficients of the Taylor series of the residual R(x). The principle is to vanish these Taylor coefficients of the residual:

R(m, n) =(m + 2)(m + 1) u(m + 2, n) + (n + 2)(n + 1) u(m, n + 2) - m i=0 n j=0 g(i, j) u(m -i, n -j) + f (m, n) = 0 (5.6)
The resulting equations can be considered as recurrence formulae permitting to compute a family of polynomials that are quasi-exact solutions of the PDE. In this way, one reduces the number of unknowns from (p + 1)(p + 2)/2 to 2p + 1. The later independent unknowns are collected in a vector {v}. Then the reduced expression u p containing 2p + 1 unknowns writes:

u p =P 0 (x, y, z) + 2p+1 i=1 P i (x, y, z)v i =P 0 (x, y, z) + P {v}. (5.7)
where the first polynomial P 0 (x, y, z) balances the right-hand side f (x) of Eq. (5.1) and P collects the complete family of solutions of the associated homogeneous problems, see [START_REF] Zézé | A boundary meshless method with shape functions computed from the PDE[END_REF][START_REF] Yang | Solving large-scale problems by taylor meshless method[END_REF] for more details. As compared with Trefftz methods [START_REF] Kita | Trefftz method: an overview[END_REF], the application to non-homogeneous and non-linear equations is straightforward. chosen on Γ d and Γ n respectively (see Fig. 5.2). Then searching the variables {v} is equivalent to minimize the following function:

Boundary least-square collocation

T (v) = 1 2 x i ∈Γ d u p (x i ) -u d (x i ) 2 + ω • 1 2 x j ∈Γn |Tu p (x j ) -t n (x j )| 2 (5.8)
where ω permits to balance the two kinds of boundary conditions. This minimization leads to a linear system [K]{v} = {b} where [K] is an invertible matrix. Solving this system gives the vector {v} and therefore the approximate solution of the problem Eq. (5.1). The detailed expressions of matrix [K] and the vector {b} write:

[K] αβ = M d i=1 P α (x i ) • P β (x i ) + ω • Mn j=1 TP α (x j ) • TP β (x j ), α, β ∈ [1, 2p + 1]
(5.9)

{b} α = M d i=1 P α (x i )(u d (x i ) -P 0 (x i )) + ω • Mn j=1 TP α (x j )(t n (x j ) -TP 0 (x j )) (5.10)
where M d and M n represent the number of collocation points on the Dirichlet boundary Γ d and the Neumann boundary Γ n respectively. This simple least-square collocation procedure will be applied in this chapter, a single polynomial solution Eq. (5.7) being valid in the whole domain. Nevertheless the same collocation method works well also to discretize interface conditions in a multi-domain approach [START_REF] Yang | Solving large-scale problems by taylor meshless method[END_REF] and alternative techniques are available for the discretization of the interface [START_REF] Tampango | Coupling of polynomial approximations with application to a boundary meshless method[END_REF]. There are also finite element-based alternatives within the Trefftz community [START_REF] Kita | Trefftz method: an overview[END_REF][START_REF] Jirousek | Large finite elements method for the solution of problems in the theory of elasticity[END_REF].

Convergence when the domain has a corner

What is the behavior of the latter method when the boundary is not smooth? Let us consider a problem of plane stress isotropic elasticity in a square domain with Dirichlet-Neumann boundary conditions as pictured in Fig. 5.3. The parameters are as follows: L = B = 10, E = 1000 (Young's modulus), ν = 0.3 (Poisson's ratio) and q = 100. This problem has no analytical solution and we tried to define a reference solution by the finite element method. We consider that a 500 × 500 mesh with Q4 elements leads to a sort of reference solution. Indeed with respect to a 250×250 mesh with Q4 elements or a 100×100 mesh with Q8 elements, the difference is less than 10 -3.17 . Throughout this chapter, the effectiveness of the proposed technique will be evaluated by the maximum relative error between the approximated solution u p and the reference solution u ref which is defined by:

ε = max |u p (x) -u ref (x)| max |u ref (x)| (5.11)
This elasticity problem has been solved by the Taylor method presented in the previous sections. In this case, it was not possible to get a TMMsolution as accurate as in the Fig. 5.1. The error decreases less quickly and there is a plateau of accuracy at about 1%, what is much larger than with problems having smooth solutions and larger than the difference between the two best FEM-results. More precisely the maximal error is about 4% for a degree p = 10, 1.5% for p = 30 and only 1% for p = 50. Note also that the accuracy is better in the center of the domain ( 0.03%). Clearly the worse convergence is due to the presence of the corners and this will be confirmed by forthcoming results.

A new TMM including singular shape functions

From the numerical test of section 5.2.3, one sees the difficulty to recover the property of exponential convergence when the boundary is not smooth. To establish undoubtedly the connection between presence of corner and poor convergence, the Taylor series method of Sections 5.2.1 and 5.2.2 will be modified by including some singular shape functions. Indeed it has been mathematically proved a long time ago [START_REF] Grisvard | Elliptic problems in nonsmooth domains[END_REF][START_REF] Dauge | Elliptic boundary value problems on corner domains: smoothness and asymptotics of solutions[END_REF] that the solution of an elliptic PDE in a domain with angular points can be split in two parts: a first part combining singular functions and a second part lying in a vectorial space of smoother functions. Thus it seems consistent to approximate this second part by a truncated Taylor series.

u = 0 ϕ ∂ u ∂ n = 0 θ = 0 θ = ϕ
one gets a countable family of possible exponents α:

α j ϕ = π 2 + (j -1)π, j = 1, 2, 3, • • • (5.13)
and the expressions of the corresponding singular functions as follows

Q j (r, θ) = r α j sin(α j θ) (5.14)
Such singular solutions can be built analytically for a number of PDEs and we refer to Appendix C in the case of 2D elasticity. Next the simple way to take them into account is to combine linearly the non-singular solutions in Eq. (5.7) and the singular ones in Eq. (5.14). In the case of a homogeneous problem, this leads to the following approximation:

u p (x) = 2p+1 i=1 v i P i (x) + n j=1 v 2p+1+j Q j (r, θ) (5.15)
Next we have to determine the 2p + 1 + n unknown coefficients v i from the boundary conditions, what will be achieved by the same least-square collocation method described in section 5.2.2. Transmission conditions could be accounted in the same way in case of multi-domain discretization and there are alternative methods based on Lagrange multipliers [START_REF] Tampango | Coupling of polynomial approximations with application to a boundary meshless method[END_REF][START_REF] Yang | Solving large-scale problems by taylor meshless method[END_REF]. Note that the coefficients {v j , j ≥ 2p + 2} are often considered as quantities of interest. In fracture mechanics, the stress intensity factors belong to this set of coefficients.

Numerical applications

The modified Taylor method described in section 5.2.4 is now analyzed by considering four numerical benchmarks related to Laplace equation and 2D elasticity. Two of them have exact solutions, which permits to discuss the convergence up to a very high accuracy. These latter tests correspond to basic solutions of fracture mechanics. The two other cases are simple, but more generic, examples without exact solutions, which avoids some bias due to a too specific choice of the example. For these four simple cases a multidomain method is not necessary: the approximation Eq. (5.15) holds in the whole domain.

Laplace equation with singularity

Here, we consider the Laplace equation in a domain that is a circular segment, as shown in Fig. 5.5. The left corner induces a rather weak singularity without infinite gradient at the corner (α 1 = 3/2). The geometric parameters are as follows: Since no analytical solution is available, a reference solution has been defined numerically by using the efficient open source code FreeFem++ [START_REF] Hecht | New development in FreeFem++[END_REF]. Three different meshes with 6 nodes triangular elements were used as shown in Fig. 5.6. The finite element solution obtained by the most dense mesh (mesh three with 20576 P2 elements) is identified as a reference solution. Comparing mesh one (with 3810 P2 elements) and mesh two (with 6078 P2 elements) with mesh three respectively, one obtains an error equal to 10 -4.7 and 10 -4.62 respectively. The optimal number of collocation points has been carefully discussed in previous papers [START_REF] Tampango | Coupling of polynomial approximations with application to a boundary meshless method[END_REF][START_REF] Yang | Solving large-scale problems by taylor meshless method[END_REF]. Without singular functions, the number of degrees of freedom is 2p + 1 in 2D and (p + 1) 2 in 3D and generally, the recommended number of collocation points is about the double. The influence of the number of collocation points is rediscussed here in the presence of singular shape functions and this is illustrated in Fig. 5.7. It appears that the error remains always stable from a critical number of collocation points that is a little higher than 4(p+n). In the following, we choose 6(p+n) collocation points to ensure the convergence. When enough collocation points are adopted, the accuracy is mainly dominated by the degree of polynomial shape functions p and the number of singular shape functions n. The convergence with p and n is presented in Table 5.1. Without singular shape functions (n = 0), one converges slowly and reaches an accuracy of 10 -3.8 for a large degree p = 50. Just one or two singular shape functions permit to gain two orders of magnitude: for instance an accuracy of 10 -4 require 101 DOFs without singular functions (n = 0, p = 50) while only 22 are sufficient with one singular shape function (n = 1, p = 10). Of course such results are consistent with the state of art within fracture mechanics [START_REF] Liebowitz | Finite element methods in fracture mechanics[END_REF] or with alternative results in a meshfree framework [START_REF] Li | Combined methods for elliptic equations with singularities, interfaces and infinities[END_REF]. Table 5.1: Convergence with the number of singular shape functions and the degree of non-singular polynomial shape functions for the Laplace problem defined in Fig. 5.5.

L x = √ 3, L y = 1, c = [ √ 3

Two tests from linear elastic fracture mechanics

The two next examples are basic benchmarks of fracture mechanics. The behavior near a crack tip is more singular than in the example of section 5.3.1, which means that the exponent α is smaller (α 1 = α 2 = 1/2) and the stress is infinite at the crack tip (σ ∼ r -1/2 ). The two first coefficients in the singular part, called stress intensity factors (SIFs), are very important in fracture mechanics and the accuracy with which these coefficients are obtained will be discussed. There is an ample literature about the computation of SIFs, but also of the full set of singular coefficients, see for instance [START_REF] Costabel | A quasi-dual function method for extracting edge stress intensity functions[END_REF]. In this paper, we are focused on the connection between singular shape functions and convergence of Taylor series, but the accuracy of the SIFs will be also shortly analyzed. Some informations about the analytical calculation of the singular functions are recalled in Appendix C for completeness. In this section, two fracture mechanics problems are considered. They are designed to have an analytical solution, which permits to measure very high accuracies. The solution of the first one is nothing but the famous mode I (or crack opening mode) and it is one of the singular shape functions. As for the solution of the second problem, it comes from the problem of a cracked infinite plate under uniaxial tension. Note that the maximal error in Eq. (5.11) remains consistent for the displacement field that is continuous, but not for the stress that is not bounded at the crack tip. Nevertheless this maximal error will be applied on a cloud of 50 × 50 points uniformly distributed on the domain except at the crack tip, where the control point has been moved of 0.001 (typically from x = 0.2 to x = 0.201).

A test to recover the crack opening model

The physical model is two-dimensional plane strain elasticity. The famous mode I of fracture mechanics corresponds to the following stress field: exactly given by Eqs. (5.16) and (5.17). The numerical value of the Poisson ratio is ν = 0.3 and the crack length is a = 0.2. Table 5.2: Convergence with the degree p and the number of singular functions n for the mode I problem described in Fig. 5 The Taylor method of section 5.2.4 has been applied to this problem, without and with singular shape functions. As expected, it is not possible to get accurate solutions without singular solutions, the errors in displacement being always larger than 50%, even for large degrees p = 30 and p = 50: hence it seems difficult to use only polynomial shape functions to recover a singular solution. Some results with singular shape functions are presented in Table 5.2. When two singular shape functions (α 1 = α 2 = 1/2) and three polynomial shape functions (i.e. a degree p = 1) are used, the maximum relative errors on displacements and stresses are very small: 10 -15.1 and 10 -15.3 respectively. Such an error is very close to the unit round-off error 2 -52 ≈ 2.22 × 10 -16 for a single real number within the double-precision floating-point format [START_REF] Muller | Handbook of Floating-Point Arithmetic[END_REF][START_REF] Kansa | On the ill-conditioned nature of C ∞ RBF strong collocation[END_REF]. If one increases the degree, the accuracy deteriorates slightly, because of the propagation of round-off errors. The same very specific scheme was also observed if one tries to get a polynomial solution by the pure Taylor method of section 5.2.1. This behavior is typical of problems where the exact solution lies in the vectorial space generated by the shape functions [START_REF] Yang | Solving large-scale problems by taylor meshless method[END_REF].

             σ x = K I √ 2πr

The infinite plate with a central crack under tensile stress

The second fracture benchmark emanates from the problem of an infinite medium with a horizontal central crack of length 2a and submitted to a uniaxial stress σ ∞ e y ⊗ e y at infinity, see Fig. 5.9a. The analytical solution is well known [START_REF] Sih | On the Westergaard method of crack analysis[END_REF][START_REF] Muskhelishvili | Some basic problems of the mathematical theory of elasticity[END_REF]: where the r 1 , θ 1 and r 2 , θ 2 are defined as follows:

2a x y σ ∞ σ ∞ ∞ ∞ B L (b)
   r 1 = x 2 + y 2 ; θ 1 = arctan( y x ).

r 2 = (x 2y 2a 2 ) 2 + (-2xy) 2 ; θ 2 = arctan( -2xy x 2y 2a 2 ).

(5.20)

The chosen parameters are as follows: a = 0.2, σ ∞ = 100, ν = 0.3, κ = 3 -4ν, E = 1000 (Young's modulus) and µ = E/[2(1 + µ)]. A square domain has been cut in the infinite domain (see Fig. 5.9a) and a corresponding boundary value problem has been posed in Fig. 5.9b so as to keep the same solution Eqs. (5.18) and (5.19).

The convergence with the degree p of the polynomial shape functions and the number n of singular shape functions is now discussed. The numerical solution is wrong without singular shape functions, the error being at least 50%. The Table 5.3 presents the obtained accuracy from two to eight singular functions. The introduction of two singular functions is sufficient to get an error of about 1%. Next the accuracy is easily improved by increasing the number of singular functions and the degree. For instance with a degree p = 10 (42 polynomials), one passes from an error of 10 -1.4 for n = 2 to an error of 10 -2.3 for n = 6 and finally to 10 -3.3 for n = 20. There is a plateau of accuracy at a high level (10 -5 or 10 -6 ). Moreover the error in stress and Stress Intensity Factor has about the same order of magnitude as for the displacement. The convergence with these two parameters p and n is also illustrated in Fig. 5.10, where the distribution of transverse stress σ y is plotted: one sees clearly the strong improvement obtained by increasing the number of singular functions from n = 2 to n = 8. Table 5.3: Convergence with the degree p and the number of singular functions n for the cracked plate problem described in Fig. 5

Application in two-dimensional elasticity

Let us come back to the 2D elasticity problem of Fig. 5.3. The accuracy of TMM without singular functions was not better than 1%. The accuracy of the vertical displacement v with two and four singular functions is reported in Table 5.4. Two types of singular functions are used: n 1 and n 2 singular functions located at upper-left and bottom-left corners respectively. At least one order of magnitude is gained simply by adding one singular function at each corner (n 1 = n 2 = 1), the error being reduced to 0.035% for a degree p = 30 and for n 1 = n 2 = 1. This establishes once more the connection between singularities and bad convergence of Taylor series, as well as a simple manner to get tremendous improvements of this convergence. Note that the maximal error mentioned in the Table (10 -3.67 ) may be due as well to FEM as to TMM.

Conclusion

A boundary collocation meshless method involving high order polynomial shape functions and singular shape functions was proposed in this chapter to solve problems with singularities. Since these shape functions are the approximated solutions of the considered problems, the discretization concerns only the boundary and least-square collocation is employed to determine the unknown coefficients. Based on previous research, the proposed technique without singular shape functions always works well when the solution of a considered problem is smooth, while it fails to recover the solution of problem with singularities. So we wonder the probability of solving singular problem effectively by introducing the singular shape functions. It was well established [START_REF] Grisvard | Elliptic problems in nonsmooth domains[END_REF] that the solution of a problem in an angular domain can be split in a smooth part and a linear combination of singular functions: in the present method, the smooth part is discretized by Taylor series. Few benchmark tests established that the accuracy of the approximated solution is widely improved by introducing few singular shape functions. When enough collocation points are chosen, the accuracy is dominated by two factors: the degree of polynomial shape functions and the number of singular shape functions. For a fixed number of singular shape functions, the proposed technique converges with the degree of polynomial shape functions. And for a fixed degree of polynomial shape functions, the accuracy increases with the number of singular shape functions.

Chapter 6

Conclusion and perspectives

Concluding remarks

A boundary collocation meshless method based on Taylor series, named Taylor Meshless Method (TMM), has been discussed in this thesis. The specificity of TMM is the numerical computation of a family of polynomial shape functions that are solutions of the considered PDE in the sense of Taylor series. Since the PDE has been satisfied inside the domain, the discretization concerns only the boundary. Then collocation technique associated with least-square method has been employed to determine the unknown coefficients.

In Chapter 2, the interest has been focused on the treatment of boundary and interface conditions in piecewise resolutions. Least-square collocation has been validated for boundary conditions while bridging techniques based on Lagrange multipliers have been validated for transmission conditions by Tampango et al. [START_REF] Tampango | Coupling of polynomial approximations with application to a boundary meshless method[END_REF]. In this chapter, we verified the possibility of using Lagrange multiplier method to account for boundary conditions and re-visiting least-square collocation method to account for interface conditions. Based on the numerical results, least-square collocation is confirmed to be very easy, efficient and robust for accounting for both boundary conditions and interface conditions.

In chapter 3, the application of TMM has been extended to three dimensions. The computation time and ill-conditioning have been well analyzed. The results show that the accuracy of TMM is governed by two antagonistic effects: exponential convergence of Taylor series with the degree while propagation of round-off errors facilitated by ill-conditioned matrix increases the error. As for computation time, the most time-consuming part in TMM is the calculation of shape functions at all collocation points and this tedious work can be issued by parallel computing. The considered large-scale cases show the possibility of combining a rather large degree of Taylor series and many sub-domains.

In chapter 4, TMM associated with Newton method and Automatic Differentiation is proposed to solve non-linear elliptic PDEs. The numerical tests demonstrate that non-linear TMM keeps the same exponential convergence and robustness features as the linear one. Non-linear numerical applications in one dimension and three dimensions illustrate that the maximum relative error decreases exponentially with the refinement of discretization of TMM and the number of Newton iterations.

In Chapter 5, singular shape functions are introduced in TMM to account for problems with singularities. It was well established that the solution of a problem in an angular domain can be split into a smooth part and a linear combination of singular functions. In this chapter, the polynomial shape functions deduced from Taylor series were employed to account for the smooth part solutions while singular part solutions were recovered by using a few analytically known singular solutions. Based on the numerical tests, the accuracy of the approximated solution is widely improved by introducing a few singular shape functions.

In this thesis, the applications of TMM are extended from two dimensions to three dimensions, from linear cases to non-linear cases, from one Taylor series or a few Taylor series to hundreds of sub-domains. TMM is confirmed to be efficient, fast convergent and robust.

Future perspective

The present thesis leaves a number of challenging open issues that are interesting to address in future research.

As mentioned previously, the most time-consuming part in TMM is the calculation of shape functions at all collocation points. Fortunately, since the shape functions of different collocation points are independent of each other, this time-consuming part can be realized by using parallel computing. Only small-scale parallel computing is tested in the present thesis. With the rapid development of the computing technique, larger-scale parallel computing is worthy of further study.

Although TMM has many advantages, for instance the exponential convergence, the strong reduction of the number of degrees of freedom and the advantages what meshles methods possess, it is not as versatile as the classical numerical methods, for instance finite element method. One may wonder the possibility of combining TMM with FEM and play their respective advantages. This will further extend the applications of TMM and ensure the accuracy of approximations while improving the efficiency of computing.

Although the algorithm discussed in this thesis is very simple, it is still difficult for potential users to solve practical engineering problems. If TMM is combined with graphics software and packaged into a user-oriented program software, the more value of TMM will be achieved.

Appendix A Discretization of the Laplace operator

Few formulae are presented to define the matrices in section 3.2.

[L k ] = [G k+1 x ][G k+2 x ] + [G k+1 y ][G k+2 y ] + [G k+1 z ][G k+2 z ] (1)
where g n-j-2 β j )u 1 (14) Comparing Eq.( 13) with Eq.( 14), one can get the iterative relationships for γ n , α n and β n .

[G k x ] =        G k xe 0 0 • • • 0 0 0 G k-1 xe 0 • • • 0 0 0 0 G k-2 xe 0 0 . . . . . . . . . . . . 0 0 0 G 1 xe 0        (2)
[G k xe ] =        k 0 0 • • • 0 0 0 k -1 0 • • • 0 0 0 0 k -2 0 0 . . . . . . . . . . . . 0 0 0 1 0        (3)
                               α 0 = 1, α 1 = 0, α i = 1 i(i -1) i-2 j=0
g i-j-2 α j ∀i, 2 i p; β 0 = 0, β 1 = 1,

β i = 1 i(i -1) i-2 j=0
g i-j-2 β j ∀i, 2 i p; γ 0 = 0, γ 1 = 0,

γ i = 1 i(i -1) (-f i-2 + i-2 j=0 g i-j-2 γ j ) ∀i, 2 i p. (15)
Then, u(x) can be written as u(x) = P 0 + P {v} = P 0 + P 1 u 0 + P 2 u 1 [START_REF] Micchelli | Interpolation of scattered data: distance matrices and conditionally positive definite functions[END_REF] where P 0 = p i=2 γ i x i , P 1 = 1 + p i=2 α i x i and P 2 = x + p i=2 β i x i .

So the approximate solution of the problem is completely obtained as a function of u 0 and u 1 . Last the two unknowns u 0 and u 1 can be deduced from the boundary conditions.

[S wk ] = 0 (2k+1)×k(k-1)/2 I k(k-1)/2×k(k-1)/2 [START_REF] Hardy | Theory and applications of the multiquadric-biharmonic method 20 years of discovery[END_REF]

Figure 1 . 1 :

 11 Figure 1.1: Sketch for boundary collocation.

Figure 1 . 2 :

 12 Figure 1.2: The influence of the number of collocation points in problem (1.5).

1

 1

( 1 . 7 )Figure 1 . 3 :

 1713 Figure 1.3: Sketch for collocation in a rectangular domain.

1 . 4 ,-

 14 x = y = 1) as follows: µ∆u + ∂q/∂x = f x -µ∆v + ∂q/∂y = f y ∂u/∂x + ∂v/∂y = 0 (1.9)

Figure 1 . 4 :

 14 Figure 1.4: Sketch of boundary collocation for Stokes equations.

0

  TMM using least-square method for B. C. TMM using Lagrange multipliers method for B. C.

Figure 2 . 1 :

 21 Figure 2.1: The p-convergence of TMM models for Laplace Eq. (2.14).

Figure 2 . 2 :

 22 Figure 2.2: The influence of the number of collocation points on the convergence for the Laplace problem (2.14).

Figure 2 . 3 :

 23 Figure 2.3: The influence of the effective radius of radial functions in TMM with Lagrange multipliers for problem (2.14). d = n * δd. ∆d represents the minimum distance between two adjacent collocation points. M = 4p.

Figure 2 . 4 :

 24 Figure 2.4: The influence of the number of collocation points on the convergence for Laplace problem Eq. (2.17).

Figure 2 . 5 :

 25 Figure 2.5: The influence of the effective radius d for Eq. (2.17). d = n * δd and M = 4p, where δd represents the minimum distance between two adjacent collocation points.

  2.6). A splitting in several sub-domains can be achieved in the same way.

Figure 2 . 6 :

 26 Figure 2.6: Two sub-domains for the Laplace problem (2.19).

  ) is reconsidered with a singularity at x 0 = [1.0, 1.0]. The domain is subdivided in two parts as shown on Fig.2.7.

3 θFigure 2 . 8 :

 328 Figure 2.8: Laplace equation (2.17) in the unit disk. Solutions on the inner circle of radius 0.8, p 1 = p 2 = p = 10, M 1 = M 2 = 2p, M r = p.

Figure 2 . 10 :

 210 Figure 2.10: Discretization of the crown for a piecewise resolution.

Figure 2 . 11 :

 211 Figure 2.11: The convergence with the number of sub-domains (p = 10).

Figure 2 . 12 :

 212 Figure 2.12: The p-convergence solved by TMM with six sub-domains.

  2.13).

Figure 2 . 13 :

 213 Figure 2.13: Two sub-domains for Laplace problem (2.28).

d p = 10 p

 10 = 20 p = 30 p = 40 p = 50 1 * δd -1.5550 -3.2542 -5.0940 -6.3572 -6.2538 3 * δd -1.6715 -3.4039 -5.1856 -6.3524 -7.6927 6 * δd -1.6236 -3.3400 -4.9338 -5.4374 -5.5595

d p = 10 p

 10 = 15 p = 20 p = 25 p = 30

  2.14).

Figure 2 .

 2 Figure 2.14: A first 2D elasticity problem

Figure 2 .

 2 Figure 2.15: A second 2D elasticity problem

Figure 2 . 16 :

 216 Figure 2.16: Problem of Fig. 2.15 solved with one domain. The relative error of v in the whole domain with the degree p=10.

Figure 2 . 17 :

 217 Figure 2.17: Problem of Fig. 2.15 solved with one domain. The p-convergence of the relative error at the point X = [10, 10] far away from the singularity.

[

Figure 2 . 18 :

 218 Figure 2.18: TMM model with sub-domains. The red points are the development points of each sub-domain.

Figure 2 . 19 :

 219 Figure 2.19: Sketch for the shape of domain, problem Eq. (2.36).

Figure 2 . 20 :

 220 Figure 2.20: Sketch for the plate divided into 100 sub-domains, problem Eq. (2.36).

Figure 2 . 21 :

 221 Figure 2.21: The p-convergence for problem Eq. (2.36). The size of the domain is 10 × 10 × 1 and with 100 TMM-elements coupled by three different bridging techniques.

Figure 3 . 1 :

 31 Figure 3.1: Sketch for two sub-domains.

u 1

 1 Poisson ratio and shear modulus are denoted by ν and µ respectively.The displacement values u 1 , u 2 and u 3 in three coordinate directions of point P are approximated by linear combinations of fundamental solutions: (a, b, c, Q; P ) = aG 11 (P, Q) + bG 12 (P, Q)+ cG 13 (P, Q) u 2 (a, b, c, Q; P ) = aG 21 (P, Q) + bG 22 (P, Q) + cG 23 (P, Q) u 3 (a, b, c, Q; P ) = aG 31 (P, Q) + bG 32 (P, Q) + cG 33 (P, Q) (3.40)where the coefficients a, b and c denote three coordinate components of the point load acting at the source point Q. In the following, we choose values of elastic constants appropriate for steel namely as recommended in[START_REF] Poullikkas | The method of fundamental solutions for three-dimensional elastostatics problems[END_REF], ν = 0.3 and µ = 1.15 (10 5 MPa). The coordinate components of the point load acting at the source point Q is set to be a = b = c = 1 (10 5 N) and different locations of the source points are tested in the following.

Figure 3 . 2 :

 32 Figure 3.2: Sketch for the shape of domain, problem of section 3.3.4.

Figure 3 . 3 :

 33 Figure 3.3: Sketch for the plate divided into 100 sub-domains, problem of section 3.3.4.

Figure 3 . 4 :

 34 Figure 3.4: Influence of number of collocation points M for degrees p = 10, 20, 30 to recover the polynomial solution of problem Eq. (3.26).

Figure 3 . 5 :Figure 3 . 6 :

 3536 Figure 3.5: Influence of number of collocation points M for degrees p = 10, 20, 30 to recover the singular solution of the Laplace problem Eq. (3.28) according to the source point X 0 = [1, 1, 1].

Figure 3 . 7 :

 37 Figure 3.7: Exponential convergence with the degree for Laplace problems Eq. (3.28). Two cases are considered according to the source point X 0 .

Figure 3 . 8 :

 38 Figure 3.8: The relative error versus the condition number for problem Eq. (3.26).

  Source point [2.0,2.0,2.0] Source point [3.0,3.0,3.0] Source point [4.0,4.0,4.0] Source point [5.0,5.0,5.0] Source point [10.,10.,10.] Source point [100,100,100] Linear fit of the results y = 0.2404x -16.9560

Figure 3 . 9 :

 39 Figure 3.9: The relative error versus the condition number for problem Eq. (3.28).

  Source point [2.0,2.0,2.0] Source point [3.0,3.0,3.0] Source point [4.0,4.0,4.0] Source point [5.0,5.0,5.0] Source point [10.,10.,10.] Source point [100,100,100] Linear fit of the results y = 0.2092x -15.2550

Figure 3 . 10 :

 310 Figure 3.10: The relative error versus the condition number for system (3.30, 3.31, 3.32).

Figure 3 . 11 :

 311 Figure 3.11: Comparing the limits of convergence.

Figure 3 . 12 :

 312 Figure 3.12: CPU time for shape functions versus p.

Figure 3 . 13 :

 313 Figure 3.13: Test of section 3.3.4, the approximation on the plane z = 0.8.

Figure 3 . 14 :

 314 Figure 3.14: Test of section 3.3.4, The exact solution on the plane z = 0.8.

Figure 3 . 15 :

 315 Figure 3.15: Test of section 3.3.4, The relative error on the plane z = 0.8.

Figure 4 . 1 :

 41 Figure 4.1: Sketch for two sub-domains.

  y, d) = sin(πx)sin(πy)

Figure 4 . 2 :

 42 Figure 4.2: Influence of number of collocation points for problem Eq. (4.19).

Figure 4 . 3 :

 43 Figure 4.3: p-convergence for two Laplace problems Eq. (4.19) that differ by the position of the singularity.

  Analytical solution for λ=1.0 Approximated solution for λ=1.0 Analytical solution for λ=1.5 Approximated solution for λ=1.5 Analytical solution for λ=2.0 Approximated solution for λ=2.0

Figure 4 . 4 :

 44 Figure 4.4: The values of u for Bratu's equation (4.24) obtained by TMM with p = 20.

40 λ

 40 Degree of polynomial p = 10 p = 20 p = 30 p = = 1.0 -3.1302 -6.9192 -9.7169 -12.3561 log 10 (maximum error) λ = 1.5 -2.0113 -4.5520 -7.0990 -9.5040 λ = 2.0 -1.3025 -3.1504 -5.0475 -7.1640

Figure 4 . 5 :

 45 Figure 4.5: The residual of Bratu's equation (4.24) versus the error for different degrees (p ∈ [10, 40]).

Figure 4 . 6 :

 46 Figure 4.6: Problem Eq. (4.26). Convergence with Newton iterations and with the discretization (p-convergence).

Figure 4 . 7 :

 47 Figure 4.7: p-convergence for two non-linear problems Eq. (4.27) that differ by the singularity X 0 .

( 4 . 32 )

 432 where f (x, y, z) = sin(πx) sin(πy)sinh( √ 2πz)/sinh( √ 2πd) and Ω is a rectangular hexahedron with sides x × y × d = 10 × 10 × 1, see Fig.4.8. The exact solution reads u ex = sin(πx) sin(πy)sinh(

Figure 4 . 8 :

 48 Figure 4.8: Sketch for the shape of domain, problem Eq. (4.32).

Figure 4 . 9 :

 49 Figure 4.9: Sketch for the plate divided into 100 sub-domains, problem Eq. (4.32).

Figure 5 . 1 :

 51 Figure 5.1: Relative error versus conditioning for a degree varying from p = 5 up to p = 60. One tries to recover fundamental solutions within 3D linear elasticity. Three cases were considered, from[START_REF] Yang | Solving large-scale problems by taylor meshless method[END_REF].

Figure 5 . 2 :

 52 Figure 5.2: Sketch for the boundary collocation method.

Figure 5 . 3 :

 53 Figure 5.3: A two-dimensional elasticity problem.

Figure 5 . 4 :

 54 Figure 5.4: Laplace equation with Neumann-Dirichlet boundary conditions.

, - 1 ]Figure 5 . 5 :

 155 Figure 5.5: Laplace equation in a circular segment. The solution is singular at the left corner.

1 (Figure 5 . 6 :

 156 Figure 5.6: Three meshes of the finite element models for the Laplace problem described in Fig. 5.5.

p 1 Figure 5 . 7 :

 157 Figure 5.7: The influence of the number of collocation points for the Laplace problem in Fig. 5.5, according to the degree p and the number n of singular shape functions

log 10 (

 10 ε) p = 10 p = 20 p = 30 p = 40 p = 50 n = 0 -2.3498 -3.0184 -3.4256 -3.6733 -3.8487 n = 1 -3.8493 -4.7543 -5.0174 -5.0227 -5.0246 n = 3 -4.4218 -4.6921 -5.0171 -5.0225 -5.0217 n = 5 -4.6823 -4.7546 -5.0174 -5.0221 -5.0196

  Figure 5.8: A boundary value problem to recover the crack opening mode.

Figure 5 . 9 :2- 1 /2 2 sin(θ 1 + θ 2 / 2 ) -x 2 v2- 1 /2 2 cos(θ 1 + θ 2 / 2 - 1 /2 2 •

 591212221212212 Figure 5.9: (a) The infinite plate with a central crack and submitted to a unixial stress at infinity; (b) A boundary value problem in a finite domain whose solution is the same as for the infinite plate of Fig. 5.9a.

Figure 5 . 10 :

 510 Figure 5.10: The distribution of transverse stress σ y and of the corresponding error along the horizontal line y = 0.

Table 5 . 4 :n 1 = n 2

 5412 Convergence with degree in the presence of few singular functions for the square elasticity shown in Fig. 5.3. n p Error of v (log 10 (ε))

  [START_REF] Duarte | Hp clouds: A meshless method to solve boundary-value problems[END_REF] into Eq.(12), one can get:n(n -1)u n = (-f n-2 +

[

  S vk ] = I (2k+1)×(2k+1) 0 k(k-1)/2×(2k+1)

Table 1 .

 1 .2. 2: Polynomial solutions of Laplace equation P 2p-2 • x -P 2p-1 • y P 2p+1 = P 2p-2 •y+P 2p-1 •xThen the approximated solution of Laplace equation can be expressed by -6 -1.3. Taylor Meshless Method the linear combination of the polynomial solutions shown in Table1.2.

	n	Re(x + iy) n	Im(x + iy) n
	0	P 1 = 1
	1	P 2 = x	P 3 = y
	2 . . .	P 4 = x 2 -y 2 . . .	P 5 = 2xy . . .
	p -1 p	P 2p-2 P 2p =	P 2p-1

Table 1 .

 1 3: For problem (1.7), The comparison between TMM and FEM.

	TMM Order of Taylor series p = 10	p = 15	p = 20	p = 30
		log 10 (max E)	-2.8018 -5.9848 -9.4708	-9.4532
		Degrees of freedom	21	31	41	61
	FEM	Mesh (n x ×n y )	30 × 10 60 × 20 90 × 30 120 × 40
		log 10 (max E)	-4.1210 -5.2704 -5.9561	-6.4464
		Degrees of freedom	981	3761	9341	14721

Table 1 .

 1 4: Convergence analysis with the degree of Taylor series in case of Stokes equations.

	Taylor series	p	5	10	15	20	25	30
		u -0.8196 -3.8717 -6.5825 -11.3130 -11.0942 -10.2116
	log 10 (max E)	v -0.8885 -3.6310 -6.3422 -11.0601 -11.0125 -10.3790
		q 0.4908 -2.2985 -4.8797 -9.1742	-8.7675	-7.9215
	Total DOFs		22	42	62	82	102	122

Table 2 .

 2

	1: The p-convergence of TMM with Lagrange multipliers for Eq.
	(2.17). x 0 = [1, 1], M = 4p and d = 3 * δd.
	Degree p	10	15	20	25	30	40	60
	log 10 (E) -1.609 -2.398 -3.016 -3.783 -4.639 -6.032 -8.952

Table 2 .

 2 2: Laplace problem (2.17) split into two sub-domains. The influence of the number of coupling nodes M r . p 1 = p 2 = p. The maximal error of the whole domain.

	Degree M r log 10 (maximum error)
		6	-0.8724
	p = 10	8	-1.8752
		20	-1.7609
		12	-1.9941
	p = 20	16	-3.6059
		40	-3.6191
		22	-4.7267
	p = 30	28	-5.2471
		60	-5.3946

Table 2 .

 2 3: Laplace problem (2.17) split into two sub-domains. Coupling polynomials with different degrees by least-square method. The maximal error of the whole domain. M 1 = 2p 1 and M 2 = 2p 2 . p 1 p 2 M r log 10 (maximal error) p 1 p 2 M r log 10 (maximal error)

		6	-0.1261		6	-0.0261
	5 15	8 10	-1.5376 -2.3962	12 15	10 12	-2.3342 -2.7563
		20	-2.3150		30	-2.6764
		6	-0.0368		8	-0.9729
	9 15	10 12	-2.4050 -2.7529	15 15	10 12	-2.2286 -2.7885
		28	-2.6682		30	-2.7231
	2.4.1.3 Application to Poisson equation in a crown
	Here we consider a Poisson's problem:

Table 2 .

 2 5: The maximal relative error of the solutions of Laplace problem (2.17) obtained by various methods.

	Degree p	p = 10 p = 20 p = 30
	Single domain (section 2.3.3)	-1.609 -3.016 -4.639
	Least-square/Lagrange (Method of reference [51])	-1.578 -3.081 -4.767
	Least-square/Least-square -1.893 -3.667 -5.331 (section 2.4.1.2)
	Lagrange/Least-square (Present method)	-1.672 -3.404 -5.186

Table 2 .

 2 8: 2D elasticity problem shown in Fig. 2.14. The error shown in table is log 10 (maximal error).

	Degree p	1	3	5	6	10	15

error of u -15.612 -15.135 -14.737 -12.587 -11.670 -10.167 error of v -16.286 -15.279 -14.525 -12.627 -11.870 -10.557

Table 2 .

 2

	p 1	p 2	Max error of u	Max error of v
	5	5	0.018	0.093
	10	10	0.011	0.020
	15	10	0.012	0.018
	20	15	0.0055	0.012
	20	20	0.004	0.010

10: . The p-convergence of TMM with two sub-domains for the 2D elasticity problem shown in Fig.

2

.15.

Table 2 .

 2 11: The p-convergence of TMM with three sub-domains for the 2D elasticity problem shown in Fig. 2.15.

	p 1 p 2 p 3	Max error of u	Max error of v
	5	5 5	0.016	0.045
	10 10 10	0.009	0.02
	15 15 15	0.0055	0.012
	20 15 15	0.0045	0.009
	20 10 15	0.0045	0.008
	20 15 20	0.0045	0.009

Table 2 .

 2 12: The three versions of TMM are compared with the finite element code FreeFem++. Problem (2.36), l x × l y × d = 10 × 10 × 1.

			TMM		FreeFem++
		Least-square Lagrange multiplier Arlequin
	Degree DOFs	10 12100	11 21600	10 19300	\ 417501
	log 10 (Error)	-3.1917	-3.4995	-3.1416	-3.3425
	Time (s)	5.08	13.34	14.77	10.31
	log 10 (Cond(K))	10.9	14.2	15.9	32.7

Table 2 .

 2

	10 × 10 × 1 20 × 20 × 1 30 × 30 × 1	100 400 900	12100 48400 108900	-3.24 idem idem	10.7 idem idem

13: One fixes the size and degree p = 10 of each TMM-element and increases the size of the domain. The sub-domains are connected by using least-square collocation. Domain Number of elements DOFs log 10 (Error) log 10 (Cond(K))

  Contents 3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . 3.2 Algorithm for Taylor Meshless Method . . . . . . . . 3.2.1 Algorithm to compute the shape functions . . . . . . . 3.2.2 Boundary least-square collocation . . . . . . . . . . . . 3.2.3 Piecewise resolution . . . . . . . . . . . . . . . . . . . 3.3 Numerical examples . . . . . . . . . . . . . . . . . . . . 3.3.1 Laplace equation with polynomial solution . . . . . . .

3.3.2 Laplace equation with singular solution . . . . . . . . 3.3.3 3D elasticity . . . . . . . . . . . . . . . . . . . . . . . . 3.3.4 A very large-scale test . . . . . . . . . . . . . . . . . .

3.4 Convergence and conditioning . . . . . . . . . . . . . . 3.4.1 Influence of the number of collocation points . . . . . 3.4.2 Exponential convergence . . . . . . . . . . . . . . . . . 3.4.3 Piecewise resolutions . . . . . . . . . . . . . . . . . . . 3.4.4 More about conditioning . . . . . . . . . . . . . . . . . 3.5 Computation time . . . . . . . . . . . . . . . . . . . . . 3.5.1 Analysis of the computation time . . . . . . . . . . . . 3.5.2 First comparison with FEM . . . . . . . . . . . . . . . 3.5.3 Large boxes submitted to sinusoidal loading . . . . . . 3.6 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . .

Table 3 .

 3 1: Computational features for the 3D elastic problem Eq. (3.30-3.32) with X 0 = [1, 1, 1], according to the degree.

	Taylor Meshless Method	p = 10 p = 20 p = 30 p = 40 p = 50 p = 60
	Number of degrees of freedom	363	1323	2883	5043	7803	11163
	Total CPU time (s)	0.04	0.39	2.74	14.26	51.97	144.25
	Parallel calculation 4 cores (s)	0.37	0.42	2.73	12.11	40.78	108.86
	Relative error (Log 10 )	-2.2747 -4.4167 -6.6648 -9.0538 -8.5838 -7.2927
	Condition number of K (Log 10 )	8.48	15.69	22.46	27.76	33.56	37.54

Table 3 .

 3 2: Computational features for the singular solution of the Laplace problem Eq. (3.28) with X 0 = [1, 1, 1], according to the degree. Two subdomains connected by least-square collocation method.

	Taylor Meshless Method	p=10	p=20	p=30	p=40	p=50	p=60
	Number of DOFs	242	882	1922	3362	5202	7442
	Total CPU time (s)	0.22	0.49	2.53	10.10	32.44	87.43
	Relative error (log 10 )	-2.8691 -5.6706 -7.1521 -7.5321 -6.8928 -7.1045
	Cond(K) (log 10 )	10.76	21.40	28.22	31.30	35.21	38.81

Table 3 .

 3 3: Computational features for the singular solution of the Laplace problem Eq. (3.28) with X 0 = [1, 1, 1], according to the degree. Two subdomains connected by Lagrange multipliers method.

	Taylor Meshless Method	p=10	p=20	p=30	p=40	p=50	p=60
	Number of DOFs	644	2484	5524	9764	15204	21844
	Total CPU time (s)	0.10	0.84	4.83	19.01	62.90	192.60
	Relative error (log 10 )	-2.6576 -5.2419 -7.3609 -7.5256 -6.5958 -5.5101
	Cond(K) (log 10 )	21.06	22.56	23.16	27.73	31.67	36.61

Table 3 .

 3 4: Computational features for the singular solution of the Laplace problem Eq. (3.28) with X 0 = [1, 1, 1], according to the degree. Two subdomains connected by the weak sense coupling method.

	Taylor Meshless Method	p=10	p=20	p=30	p=40	p=50	p=60
	Number of DOFs	308	1208	2708	4808	7508	10808
	Total CPU time (s)	0.15	0.70	4.01	15.69	50.11	142.08
	Relative error (log 10 )	-2.8645 -5.7040 -6.6714 -6.6013 -5.1441 -5.9329
	Cond(K) (log 10 )	12.01	19.50	24.17	27.21	30.56	33.28

Table 3 .

 3 5: Convergence analysis for problem Eq. (3.28).

	Taylor Meshless Method p = 10 p = 20 p = 30	p = 40	p = 50 p = 60
	Number of DOFs	121	441	961	1681	2601	3721
	Collocation points	288	968	2048	3528	5408	7688
	Relative error (Log 10 )	-2.6617 -5.1167 -7.6343 -10.2583 -9.6329 -7.8021
	Cond(K) (Log 10 )	6.67	12.95	19.13	25.28	31.39	37.57

Table 3 .

 3 6: Analysis of CPU time for problem Eq. (3.28).

	Taylor Meshless Method	p = 10 p = 20 p = 30 p = 40 p = 50 p = 60
	Number of DOFs	121	441	961	1681	2601	3721
	Total CPU time (s)	0.015	0.10	0.56	1.96	5.94	18.66
	of which	time for P (x) (%) 71.00	84.5	85.7	85.1	85.8	84.5
		time to K (%)	3.5	5.6	7.7	9.7	10.3	8.7
		time for K\F (%)	4.0	2.0	2.5	2.7	2.7	6.0
	Parallel calculation (4 cores) (s)	0.05	0.09	0.31	1.15	3.81	10.45

Table 4 .

 4 2: The p-convergence of problem Eq.(4.22).

	Degree of polynomial	p = 5	p = 10 p = 15 p = 18
	Number of DOFs	3600	12100	25600	36100
	log 10 (max error)	-0.9140 -3.2407 -6.5568 -8.8391

Table 4 .

 4 3:The convergence with the number of iterations and the degree of polynomials for problem Eq.(4.23).

	Degree of polynomial	Iteration p = 10 p = 20 p = 30	p = 40
		1	-0.9031 -0.9031 -0.9031 -0.9031
	log 10 (maximum residual)	2	-2.153 -2.7844 -2.7844 -2.7844
		3	-1.9376 -4.9503 -6.6041 -6.6041
		4	-1.9376 -4.9543 -8.4665 -12.7815
		8	-1.9376 -4.9543 -8.4665 -12.7817

Table 4 .

 4 4: The p-convergence for Bratu's problem (4.24) with different λ.

Table 4 .

 4 5: The maximal relative error of each iteration for the 3-D problem Eq.(4.26).

	Degree of polynomial	p = 5	p = 8	p = 10	p = 15	p = 20	p = 30
	Number of DOFs	36	81	121	256	441	961
	Number of iterations	1	2	2	3	3	3
	log 10 (max error)	-2.3103 -4.7874 -6.6143 -11.9125 -14.6874 -14.8070
	log 10 (max residual)	-1.0902 -3.0926 -4.6054 -9.1402 -14.6186 -15.5758

Table 4 .

 4 6: p-convergence of problem Eq.(4.27) with singularity at X 0 = [1, 1, 1].

	Degree of polynomial p = 10 p = 20 p = 30	p = 40	p = 50
	Number of iterations	2	3	3	4	4
	log 10 (max error)	-2.6448 -5.2524 -7.4341 -10.0566 -12.2454
	log 10 (max residual)	-2.1183 -5.9682 -6.6592 -10.6087 -12.4773
	CPU-time (s)	0.09	1.15	10.85	94.76	434.47

Table 4 .

 4 7: p-convergence of problem Eq.(4.29). We consider the non-linear convection equation in a spherical domain (Ω = {x, y, z|x 2 + y 2 + z 2 1}):

	Degree of polynomial	p=5	p=8	p=10	p=20	p=30
	Number of DOFs	36	81	121	441	961
	Number of iterations	1	2	2	3	3
	log 10 (max residual)	-2.3345 -3.9346 -4.4316 -9.5295 -15.8210
	p = 20 and 30. With a sufficient number of Newton iterations, the proposed
	technique converges exponentially with the degree of polynomials p and up to
	a very high accuracy.
	4.3.2.4 Problem four
	-∆u + ( ∂u ∂x + ∂u ∂y + ∂u ∂z )u = f (x, y, z) in u(x, y, z) = e x sin(πy/2) sin(πz/2) on ∂Ω. Ω,	(4.30)
	where the right-hand side denotes:
	f (x, y, z) = (1 + π)e 2x sin 2 (πy/2) sin 2 (πz/2) + (π 2 /2 -1)e x sin(πy/2) sin(πz/2) (4.31)

Table 4 .

 4 8: p-convergence of problem Eq.(4.30).

	Degree of polynomial	p = 5	p = 10 p = 15	p = 20	p = 30
	Number of iterations	2	3	4	4	4
	log 10 (max error)	-1.3483 -4.7346 -8.6229 -13.4719 -14.1924
	log 10 (max residual)	-0.5556 -3.5428 -7.7981 -11.1863 -13.5416

Table 4 .

 4 9: The p-convergence of problem Eq. (4.32) solved by TMM with 100 elements.

	Degree of polynomial	p=5	p=8	p=10	p=15	p=18
	Number of DOFs	3600	8100	12100	25600	36100
	Number of iterations	1	2	2	3	3
	log 10 (max error)	-0.8537 -2.0655 -3.1639 -6.5548 -8.8333
	log 10 (max residual)	-1.0988 -1.787 -2.7377 -4.8785 -6.7583

  .8.

	n	p NDOFs Error in displacement Error in stress
				log 10 (ε)	log 10 (ε)
		1	8	-15.1204	-15.2682
	n = 2	10	44	-12.8060	-12.6139
		20	84	-11.2045	-10.1424
		1	12	-15.1804	-15.3474
	n = 6	10	48	-12.5630	-12.2353
		20	88	-10.8113	-10.0309

  .10.

	n	p Error in displacements Error in stresses Error in SIFs
			log 10 (ε)	log 10 (ε)	log 10 (ε)
		10	-1.4137	-1.3121	-1.5543
	n = 2	30	-2.0785	-1.6041	-1.5914
		50	-2.2081	-1.5592	-1.7093
		10	-2.3236	-1.6326	-2.5541
	n = 6	30	-4.7309	-4.2851	-4.4160
		50	-5.8604	-5.1474	-5.0110
		10	-2.4972	-1.7541	-2.4210
	n = 8	30 40	-5.0535 -6.2975	-4.5884 -5.6968	-5.8802 -6.0031
		50	-5.6359	-4.8222	-5.8995

The polynomials u p in Eq. (5.7) are quasi-exact solutions of the PDE, in such a way that the discretization concerns only the boundary. As proposed in[START_REF] Zézé | A boundary meshless method with shape functions computed from the PDE[END_REF][START_REF] Fairweather | The method of fundamental solutions for elliptic boundary value problems[END_REF][START_REF] Zhang | Least-squares collocation meshless method[END_REF], a collocation technique combined with the least-square method can be used to account for boundary conditions. Two sets of nodes x i , x j are

Acknowledgments

This thesis has been written during my study as a joint-PhD candidate in Université de Lorraine and Wuhan University.

First I would like to express my deeply-felt gratitude to my parents who create the space of chasing my dreams freely. I also would like to extend my heartfelt appreciation to my supervisors, Prof. Michel Potier-Ferry and Prof. Heng Hu, whose patient guidance, valuable suggestions and constant encouragement make me successfully finish my PhD thesis. Their enthusiasm for scientific research and the attitude towards lift inspire me all the time both in academic study and daily life. They offer me selfless help that made me who I am today.

My heartfelt appreciate also extends to all my colleges in the team. Thank them for their help not only in the scientific research but also in the daily life. I am so grateful for the learning, companionship and encouragement from each other.

-I -

The other coefficients can be expressed in terms of the independent ones. They are collected in the following vector:

Computing singular solutions of PDEs by Taylor series

Abstract

Taylor Meshless Method (TMM) is a true meshless integration-free numerical method for solving elliptic Partial Differential Equations (PDEs). The basic idea of this method is to use high order polynomial shape functions that are approximated solutions of the PDE and are computed by the technique of Taylor series. Currently, this new method has proved robust and efficient and it has the property of exponential convergence with the degree, when solving problems with smooth solutions. This exponential convergence is no longer obtained for problems involving cracks, corners or notches. On the basis of numerical tests, this chapter establishes that the presence of a singularity leads to a worsened convergence of the Taylor series, but highly accurate solutions can be recovered by including a few singular solutions in the basis of shape functions.

Present chapter corresponds to the submitted research paper (Yang et

Appendix B TMM for ODE

Here, we consider an ODE:

In Eq.( 7), the functions g(x) and f (x) are known. We expand g(x) and f (x) by Taylor series up to the degree p -2 as follows.

where g i and f i are known coefficients.

The approximate solution of Eq.( 7) is sought in the form of a polynomial of degree p:

we have the following expression of derivative respect to u(x):

The product g(x)u(x) is truncated at the order p:

If the Eq.( 7) is satisfied, one can get the following relationship:

There are p-1 equations in Eq.( 12) for p+1 unknowns (u i , i = 0, 1, • • • , p): u 0 and u 1 are chosen as the independent unknowns. One can find the relationship between u n and u 0 , u 1 :

Appendix C

Singular solutions of opening crack

In this Appendix, we compute analytically the whole set of singular functions in the case of a crack tip. This is based on the well-known complex variable method and we refer to [START_REF] Muskhelishvili | Some basic problems of the mathematical theory of elasticity[END_REF] for the details. A similar approach can be done for the square domain problem for Fig. 5.3. The Airy stress function U satisfying the biharmonic equation ∆ 2 U = 0 reads:

where z = x + iy and both ϕ(z) and χ(z) are harmonic functions.

The complex representation of displacements denotes

where ν and µ represent the Poisson's ratio and the shear modulus respectively. The number κ is a constant related to ν and it takes the value 3 -4ν for plane strain while it takes (3ν)/(1 + ν) for generalized plane stress. The complex representations of stresses denote

The force acting on an arc AB then can be expressed by:

where ψ(z) = χ (z).

One sets the harmonic functions as follows: ϕ(z) = c 1 z α and ψ(z) = c 2 z α , where α is a real number and c 1 , c 2 are unknown complex numbers determined by the boundary conditions on the upper and lower surfaces of the crack tip. The force (X + iY ) vanishes along the crack tip, hence the Eq. ( 22) satisfies the following conditions:

Following relationships can be obtained from Eq. ( 23):

Considering Eq. ( 24), one can get

which means α = 1 2 , 1, 3 2 , 2, 5 2 , • • • . To ensure the singularity of the approximation, the integer values of α should be removed.

Considering the value of α deduced from the boundary conditions, one can get the following relationship between c 1 and c 2 by Eq. ( 24):

If c 1 is set to be 1 and i respectively, these corresponding values of c 2 are as follows:

Substituting the values of c 1 and c 2 in Eqs. ( 20) and ( 21), one can obtain the expressions of displacements and stresses as follows:

• The first set of singular solutions, c 1 = 1 and c 2 = 1α:

• The second set of singular solutions, c 1 = i and c 2 = -i(1 + α):

Selecting one value of α, there will be two sets of singular solutions. The way to deduce the non-singular shape functions from Taylor series has been introduced in section 5.2.1. The linear combination of the non-singular shape functions and singular solutions shown in Eqs. [START_REF] Wendland | Meshless Galerkin methods using radial basis functions[END_REF] and [START_REF] Kansa | Multiquadrics-a scattered data approximation scheme with applications to computational fluid-dynamics-II solutions to parabolic, hyperbolic and elliptic partial differential equations[END_REF] leads to the approximations of this considered problem. Since the approximations satisfy the control equations, only boundary conditions need to be considered to determine the unknown coefficients. This last step is realized by using leastsquare collocation technique as shown in section 5.2.2.