Localized method of fundamental solutions for three-dimensional inhomogeneous elliptic problems: theory and MATLAB code

In this paper we investigate the application of the localized method of fundamental solutions (LMFS) for solving three-dimensional inhomogeneous elliptic boundary value problems. A direct Chebyshev collocation scheme (CCS) is employed for the approximation of the particular solutions of the given inhomogeneous problem. The Gauss–Lobatto collocation points are used in the CCS to ensure the pseudo-spectral convergence of the method. The resulting homogeneous equations are then calculated by using the LMFS. In the framework of the LMFS, the computational domain is divided into a set of overlapping local subdomains where the traditional MFS formulation and the moving least square method are applied. The proposed CCS-LMFS produces sparse and banded stiffness matrix which makes the method possible to perform large-scale simulations on a desktop computer. Numerical examples involving Poisson, Helmholtz as well as modified-Helmholtz equations (with up to 1,000,000 unknowns) are presented to illustrate the efficiency and accuracy of the proposed method.

Keywords: Particular solutions; Method of fundamental solutions; Meshless method; Chebyshev polynomials; Inhomogeneous elliptic problems

The method of fundamental solutions (MFS), which can be viewed as a special case of Trefftz method [[1]], has emerged as a robust and suitable meshless method for the solutions of certain boundary value problems [[2]–[5]]. The properties of this method, especially when high accuracy results are required, won the favor of many researchers in the engineering as well as in the mathematical community. During the last three decades, the method has been successfully tried for, for example, various homogeneous [[6]] and inhomogeneous elliptic boundary value problems [[7]], inverse problems [[9]] and non-linear problems [[11]]. Some recent modifications and different names of the MFS can be found in Refs. [[1], [12]–[17]]. However, similar to the boundary element method (BEM) [[18]–[22]], the classical MFS formulation produces dense and non-symmetric system matrices that, although smaller in sizes, require lot of memory and CPU-time to be computed. This bottleneck presents a serious problem of the method in large-scale engineering simulations [[23]–[27]].

In our previous studies, we introduced a localized version of the MFS (LMFS) [[28]–[30]] which could be an improvement to the classical MFS for solving large-scale applied mechanics problems. In the LMFS, the whole computational domain is firstly divided into a set of overlapping local subdomains in which their corresponding geometries can simply be described by a circle (for 2D problems) and/or a sphere (for 3D problems). In each of the local subdomain, the classical MFS formulation and a moving least square (MLS) method are applied to construct the small (local) systems of linear equations. By enforcing the satisfaction of the boundary conditions at the boundary points and the governing equation at the interior nodes, the LMFS will finally yield a sparse and banded matrix system, which significantly reduces the overall storage and computational burden and can be solved quickly by using various sparse matrix solvers. Prior to this study, the LMFS has been successfully tried for 2D Laplace and biharmonic equations [[28]], 3D heat conduction problems [[30]] and 2D elasticity problems [[29]]. A fast, reliable and self-contained MATLAB code for 3D steady-state heat conduction can be found in Ref. [[30]].

In this paper, we document the first attempt to apply the LMFS for large-scale simulations of 3D inhomogeneous elliptic boundary value problems. A standard strategy for solving inhomogeneous problems is to split the final numerical solution into a particular solution and a homogeneous solution. In this paper, a direct Chebyshev collocation scheme (CCS) is introduced for the calculation of the particular solutions of the given inhomogeneous problem [[31]–[33]]. In the CCS approach, the Gauss–Lobatto collocation points are employed to ensure the pseudo-spectral convergence of the Chebyshev interpolation. Once the particular solution has been derived, the given inhomogeneous problem can be transferred to a homogeneous problem which is solved by using the LMFS. Comprehensive numerical examples associated with Poisson, Helmholtz as well as modified-Helmholtz equations in 3D materials are well-studied to investigate the accuracy and efficiency of the method. Numerical examples with up to 1,000,000 unknowns are solved successfully on a Core-i7 laptop PC using the developed CCS-LMFS code.

The reminder of the paper is organized as follows: Sect. 2 introduces the methodology of the Chebyshev collocation scheme for calculating particular solutions. The LMFS formulation and its numerical implementation for homogeneous problems are presented in Sect. 3. In Sect. 4, three benchmark numerical examples involving smooth and piecewise smooth geometries have been studied. Finally, some conclusions and remarks are provided in Sect. 5.

We here considered the following inhomogeneous elliptic boundary value problem:

Graph

subject to the following Dirichlet and Neumann boundary conditions:

2 $$u(x_1,x_2,x_3)=\overline{u}(x_1,x_2,x_3)\text{for}\left(x_1,x_2,x_3\right)\in {\mathrm{\Gamma }}_D,$$

Graph

3 $$\begin{array}{cc}\hfill \frac{\partial u(x_1,x_2,x_3)}{\partial \mathit{n}}& =\mathrm{\nabla }u(x_1,x_2,x_3)·\mathit{n}(x_1,x_2,x_3)\hfill \\ \hfill & =\overline{q}(x_1,x_2,x_3)\text{for}\left(x_1,x_2,x_3\right)\in {\mathrm{\Gamma }}_N,\hfill \end{array}$$

Graph

where $\mathrm{\Omega }$ is a 3D domain bounded by a surface $\partial \mathrm{\Omega }$ , L is the linear differential operator on $u$ over $\mathrm{\Omega }$ , $\overline{u}(x_1,x_2,x_3)$ and $\overline{q}(x_1,x_2,x_3)$ indicate the given values specified along the boundary, ${\mathrm{\Gamma }}_D$ and ${\mathrm{\Gamma }}_N$ are parts of the boundary $\partial \mathrm{\Omega }$ with $\partial \mathrm{\Omega }={\mathrm{\Gamma }}_D\cup {\mathrm{\Gamma }}_N$ and ${\mathrm{\Gamma }}_D\cap {\mathrm{\Gamma }}_N=\mathrm{\varnothing }$ , $\mathit{n}(x_1,x_2,x_3)$ is the outward unit normal vector at boundary point $(x_1,x_2,x_3)$ .

A standard strategy for solving inhomogeneous Eqs. (1)–(3) is to split the numerical solution $u$ into a particular solution $u_p$ and a homogeneous solution $u_h$ , i.e., $u=u_p+u_h$ . It should be noted that the particular solution $u_p$ should satisfy the governing Eq. (1) but does not necessarily the boundary conditions (2) and/or (3). After the particular solutions $u_p$ have been determined, the homogeneous solution $u_h$ can be derived by solving the following homogeneous boundary value problem:

4 $$L{u}_h(x_1,x_2,x_3)=0,\left(x_1,x_2,x_3\right)\in \mathrm{\Omega },$$

Graph

5 $$\begin{array}{cc}\hfill u_h(x_1,x_2,x_3)& =\overline{u}(x_1,x_2,x_3)-u_p(x_1,x_2,x_3)\hfill \\ \hfill & \text{for}\left(x_1,x_2,x_3\right)\in {\mathrm{\Gamma }}_D,\hfill \end{array}$$

Graph

6 $$\begin{array}{cc}\hfill \frac{\partial u_h(x_1,x_2,x_3)}{\partial \mathit{n}}& =\overline{q}(x_1,x_2,x_3)-\frac{\partial u_p(x_1,x_2,x_3)}{\partial \mathit{n}}\hfill \\ \hfill & \text{for}\left(x_1,x_2,x_3\right)\in {\mathrm{\Gamma }}_N.\hfill \end{array}$$

Graph

The successful evaluation of particular solutions $u_p$ could be an important step for the overall accuracy of the numerical procedure employed. In the following, a direct Chebyshev collocation scheme will be introduced for the calculation of the particular solutions. Once particular solutions have been derived, the remaining homogeneous problem is then calculated by using the LMFS.

The Chebyshev polynomials of the first kind $T_n(x)$ are given explicitly as:

7 $$T_n(x)=cos\left[narccos\left(x\right)\right],n=0,1,2,...,x\in [-1,1],$$

Graph

If $F(\xi ,\eta ,\tau )$ is of a bounded function defined in [− 1, 1] × [− 1, 1] × [− 1, 1], then it can be approximated by the following truncated Chebyshev series as [[32]]:

8 $$F(\xi ,\eta ,\tau )=\sum _{m=0}^{N_1}\sum _{n=0}^{N_2}\sum _{l=0}^{N_3}{\varphi }_{m,n,l}{T}_m(\xi )T_n(\eta )T_l(\tau ),$$

Graph

where $\left\{{\varphi }_{m,n,l}\right\}$ are unknown coefficients, $N_1$ , $N_2$ and $N_3$ stand for degrees of Chebyshev's polynomials in the $x_1$ -, $x_2$ -, and $x_3$ -directions, respectively. If $F(\xi ,\eta ,\tau )$ is known, the coefficients $\left\{{\varphi }_{m,n,l}\right\}$ can then be calculated by collocating Eq. (8) on $M=\text{(}{N}_1+1\text{)(}{N}_2+1\text{)(}{N}_3+1\text{)}$ collocation points scatted inside the cubic [− 1, 1] × [− 1, 1] × [− 1, 1].

To avoid the Runge's phenomena and guarantee the pseudo-spectral convergence of the Chebyshev polynomials, it is suggested to choose the M collocation points as the roots of Chebyshev polynomials (Gauss–Lobatto nodes). For example, the M distinct zeros ${\left\{{\xi }_i\right\}}_{i=0}^{M-1}$ of the Chebyshev polynomial $T_M(\xi )$ take the form:

9 $${\xi }_i=cos\left[\frac{(2i+1)}{2M}\pi \right],0\le i\le (M-1).$$

Graph

We now suppose that $u_p(x_1,x_2,x_3)$ , $(x_1,x_2,x_3)\in \mathrm{\Omega }$ , is the particular solution of the given inhomogeneous problem (1)–(3), where $\mathrm{\Omega }$ can be an arbitrary domain defined in ${\mathit{R}}^{\mathbf{3}}$ . Since the particular solution is not unique and does not necessarily satisfy the given boundary conditions, we have the flexibility to extend smoothly $u_p(x_1,x_2,x_3)$ to a cubic $\stackrel{⌢}{\mathrm{\Omega }}=\left[a,b\right]\times \left[c,d\right]\times \left[e,f\right]$ which contains the original computational domain $\mathrm{\Omega }$ . To use the Chebyshev collocation method as shown aforementioned, we should introduce the following change of variable to transfer the computational domain [a, b] × [c, d] × [e, f] into [− 1, 1] × [− 1, 1] × [− 1, 1]:

10 $$\begin{array}{cc}\hfill u_p(x_1,x_2,x_3)& =\sum _{m=0}^{N_1}\sum _{n=0}^{N_2}\sum _{l=0}^{N_3}{\varphi }_{m,n,l}{T}_m\left(\frac{x_1-{\beta }_1}{{\alpha }_1}\right)\hfill \\ \hfill & T_n\left(\frac{x_2-{\beta }_2}{{\alpha }_2}\right)T_l\left(\frac{x_3-{\beta }_3}{{\alpha }_3}\right),\hfill \end{array}$$

Graph

where

11 $$\begin{array}{cc}\hfill {\alpha }_1& =(b-a)/2,{\alpha }_2=(d-c)/2,{\alpha }_3=(f-e)/2,\hfill \\ \hfill {\beta }_1& =(a+b)/2,{\beta }_2=(c+d)/2,{\beta }_3=(e+f)/2.\hfill \\ \hfill \end{array}$$

Graph

It should be noted that the Gauss–Lobatto nodes defined in Eq. (9) should be also imaged into [a, b] × [c, d] × [e, f], so that the aforementioned Eq. (10) can be used. For this, we introduce the following change of variable:

12 $$x_1={\alpha }_1\xi +{\beta }_1,x_2={\alpha }_2\eta +{\beta }_2,x_3={\alpha }_3\tau +{\beta }_3,$$

Graph

where $(\xi ,\eta ,\tau )$ and $(x_1,x_2,x_3)$ denote the original and shifted/modified Gauss–Lobatto nodes, respectively. Substituting Eq. (10) into the inhomogeneous Eq. (1), we have

13 $$\begin{array}{cc}& \sum _{m=0}^{N_1}\sum _{n=0}^{N_2}\sum _{l=0}^{N_3}{\varphi }_{m,n,l}L\left[T_m\left(\frac{x_1-{\beta }_1}{{\alpha }_1}\right)T_n\left(\frac{x_2-{\beta }_2}{{\alpha }_2}\right)T_l\left(\frac{x_3-{\beta }_3}{{\alpha }_3}\right)\right]\hfill \\ \hfill & =f(x_1,x_2,x_3),\left(x_1,x_2,x_3\right)\in \stackrel{⌢}{\mathrm{\Omega }}.\hfill \end{array}$$

Graph

Noting that the function $f(x_1,x_2,x_3)$ in the right-hand side of Eq. (13) is known, the unknown coefficients $\left\{{\varphi }_{m,n,l}\right\}$ , therefore, can be calculated by collocating Eq. (13) on M shifted Gauss–Lobatto nodes (12). Once the unknown coefficients $\left\{{\varphi }_{m,n,l}\right\}$ are calculated, the particular solutions $u_p(x_1,x_2,x_3)$ at any point inside the computational domain can then be obtained by using Eq. (10).

In order to obtain the explicit LMFS formulae, an irregular cloud of points should be scattered inside the whole computational domain $\mathrm{\Omega }$ . For each of a given node ${\mathit{x}}^{(0)}=\left(x_1^0,x_2^0,x_3^0\right)$ (named as the central node), the $N_s$ nearest nodes ${\mathit{x}}^{(i)}=\left(x_1^i,x_2^i,x_3^i\right)$ , $i=1,2,...,N_s$ , around ${\mathit{x}}^{(0)}$ should be found. The concept of the "local subdomain ( ${\mathrm{\Omega }}_s$ )" then refers to the small area that contains these $N_s+1$ nodes. For each of node inside ${\mathrm{\Omega }}_s$ , the following MFS formulation should be hold:

14 $$u\left({\mathit{x}}^{(i)}\right)=\sum _{j=1}^{M_s}{\phi }_{j}G\left({\mathit{x}}^{(i)},{\mathit{s}}^{(j)}\right)={\mathit{G}}^{(i)}\mathit{\phi },{\mathit{x}}^{(i)}\in {\mathrm{\Omega }}_s,i=0,1,...,N_s,$$

Graph

where ${\left({\mathit{x}}^{(i)}\right)}_{i=0}^{N_s}$ are the $N_s+1$ collocation point inside ${\mathrm{\Omega }}_s$ , ${\left({\mathit{s}}^{(j)}\right)}_{j=1}^{M_s}$ are Ms source points which are uniformly distributed along an artificial surface ${\stackrel{⌢}{\mathrm{\Omega }}}_s$ associated with ${\mathrm{\Omega }}_s$ , $\mathit{\phi }\text{=}{\left({\phi }_1,...,{\phi }_{M_s}\right)}^T$ are unknown coefficients, and

15 $$\begin{array}{cc}\hfill G\left(\mathit{x},\mathit{s}\right)& =\frac{1}{4\pi r},G\left(\mathit{x},\mathit{s}\right)=\frac{1}{4\pi r}exp(-ikr),\hfill \\ \hfill G\left(\mathit{x},\mathit{s}\right)& =\frac{1}{4\pi r}exp(-kr),\hfill \end{array}$$

Graph

are fundamental solutions for 3D Laplace, Helmholtz and modified Helmholtz equations, respectively. In Eq. (15), the parameter $k$ denotes the frequency of the acoustical field. It should be noted that, in the LMFS, the artificial surface ${\stackrel{⌢}{\mathrm{\Omega }}}_s$ associated with each local subdomain ${\mathrm{\Omega }}_s$ can be simply chosen as a sphere with radius $R_s$ and centered at ${\mathit{x}}^{(0)}$ , as shown in Fig. 1, where $R_s$ is a parameter which should be manually determined by the user.

Graph: Fig. 1 The schematic diagram for a the artificial sphere associated with Ωs, and bMs source points uniformly distributed along the artificial sphere surface

In each of the local subdomain ${\mathrm{\Omega }}_s$ , it is possible to define the following residual function:

16 $$B(u)=\sum _{i=0}^{N_s}{\left[\left({\mathit{G}}^{(i)}\mathit{\phi }\text{-}{u}^{(i)}\right){\omega }^{(i)}\right]}^2,$$

Graph

where $u^{(i)}=u\left({\mathit{x}}^{(i)}\right)$ , and

17 $${\omega }^{(i)}=\frac{exp[-{(d_i/h)}^2]-exp[-{(d_{\text{max}}/h)}^2]}{1-exp[-{(d_{\text{max}}/h)}^2]},i=0,1,...,N_s,$$

Graph

is the Gaussian weight function as employed in Ref. [[30]]. It is noted that, in Eq. (17), $d_i=∥{\mathit{x}}^{(i)}-{\mathit{x}}^{(0)}∥$ stands for the distance between the ith supporting node and the central node ${\mathit{x}}^{(0)}$ , $d_{\text{max}}$ is the size of the local subdomain (the distance between ${\mathit{x}}^{(0)}$ and the farthest node inside the subdomain), and $h$ is a constant controlling the shape of the weight function. According to the moving least square (MLS) method, we can minimize the residual function $B(u)$ with respect to $\mathit{\phi }\text{=}{\left({\phi }_1,...,{\phi }_{M_s}\right)}^T$ as follows:

18 $$\frac{\partial B(u)}{\partial \left\{\mathit{\phi }\right\}}=0,$$

Graph

and then the following (local) system of linear equations in ${\mathrm{\Omega }}_s$ can be formed:

19 $${\mathit{A}}_s\mathit{\phi }=\mathit{b},$$

Graph

where

20 $${\mathit{A}}_s=\left[\begin{array}{ccccc}{\sum }_{i=0}^{N_s}{G}_{i1}^2{\omega }^{(i)}& {\sum }_{i=0}^{N_s}{G}_{i1}{G}_{i2}{\omega }^{(i)}& {\sum }_{i=0}^{N_s}{G}_{i1}{G}_{i3}{\omega }^{(i)}& \cdots & {\sum }_{i=0}^{N_s}{G}_{i1}{G}_{i{M}_s}{\omega }^{(i)}\\ & {\sum }_{i=0}^{N_s}{G}_{i2}^2{\omega }^{(i)}& {\sum }_{i=0}^{N_s}{G}_{i2}{G}_{i3}{\omega }^{(i)}& \cdots & {\sum }_{i=0}^{N_s}{G}_{i2}{G}_{i{M}_s}{\omega }^{(i)}\\ & & {\sum }_{i=0}^{N_s}{G}_{i3}^2{\omega }^{(i)}& \cdots & {\sum }_{i=0}^{N_s}{G}_{i3}{G}_{i{M}_s}{\omega }^{(i)}\\ & \mathit{SYM}& & \ddots & ⋮\\ & & & & {\sum }_{i=0}^{N_s}{G}_{i{M}_s}^2{\omega }^{(i)}\end{array}\right],$$

Graph

and

21 $$\begin{array}{cc}\hfill \mathit{b}& =\left[\begin{array}{cccc}{G}_{01}{\omega }^{(0)}& G_{11}{\omega }^{(1)}& \cdots & G_{N_{s}1}{\omega }^{(N_s)}\\ G_{02}{\omega }^{(0)}& G_{12}{\omega }^{(1)}& \cdots & G_{N_{s}2}{\omega }^{(N_s)}\\ ⋮& ⋮& \ddots & ⋮\\ G_{0M_s}{\omega }^{(0)}& G_{1M_s}{\omega }^{(1)}& \cdots & G_{N_s{M}_s}{\omega }^{(N_s)}\end{array}\right]\left(\begin{array}{c}\hfill u^{(0)}\\ \hfill u^{(1)}\\ \hfill ⋮\\ \hfill u^{(N_s)}\\ \hfill \end{array}\right)\hfill \\ \hfill & ={\mathit{B}}_s\left(\begin{array}{c}\hfill u^{(0)}\\ \hfill u^{(1)}\\ \hfill ⋮\\ \hfill u^{(N_s)}\\ \hfill \end{array}\right).\hfill \end{array}$$

Graph

Multiplying the inverse of ${\mathit{A}}_s$ to each side of Eq. (19), one has:

22 $$\mathit{\phi }\text{=}\left(\begin{array}{c}\hfill {\phi }_1\\ \hfill {\phi }_2\\ \hfill ⋮\\ \hfill {\phi }_{M_s}\\ \hfill \end{array}\right)={\mathit{A}}_s^{-1}{\mathit{B}}_s\left(\begin{array}{c}\hfill u^{(0)}\\ \hfill u^{(1)}\\ \hfill ⋮\\ \hfill u^{(N_s)}\\ \hfill \end{array}\right).$$

Graph

We can observe from Eq. (22) that the explicit expressions of unknown coefficients $\mathit{\phi }\text{=}{\left({\phi }_1,...,{\phi }_{M_s}\right)}^T$ can now be expressed by a linear combination of function values at the $N_s+1$ nodes located inside subdomain ${\mathrm{\Omega }}_s$ . To ensure the regularity of the matrix ${\mathit{A}}_s^{-1}$ , the number of collocation points ( $N_s+1$ ) inside ${\mathrm{\Omega }}_s$ should be slightly larger than that of the sources ( $M_s$ ), that is $N_s+1\ge M_s$ . For simplicity, in our computations $N_s+1=2M_s$ collocation nodes are chosen in each of the local subdomain ${\mathrm{\Omega }}_s$ . Now, the numerical solution at the central node ${\mathit{x}}^{(0)}$ can be calculated by substituting $\mathit{\phi }\text{=}{\left({\phi }_1,...,{\phi }_{M_s}\right)}^T$ and the coordinate of ${\mathit{x}}^{(0)}$ into Eq. (14) as follows:

23 $$\begin{array}{cc}\hfill u^{(0)}& ={\mathit{G}}^{(0)}\mathit{\phi }={\mathit{G}}^{(0)}{\mathit{A}}_s^{-1}{\mathit{B}}_s\left(\begin{array}{c}\hfill u^{(0)}\\ \hfill ⋮\\ \hfill u^{(N_s)}\\ \hfill \end{array}\right)\hfill \\ \hfill & =\sum _{j=0}^{N_s}{c}^{(j)}{u}^{(j)},\mathrm{or}{u}^{(0)}-\sum _{j=0}^{N_s}{c}^{(j)}{u}^{(j)}=0,\hfill \end{array}$$

Graph

where ${\left(c^{(j)}\right)}_{j=0}^{N_s}$ are coefficients which can be calculated by ${\mathit{G}}^{(0)}{\mathit{A}}_s^{-1}{\mathit{B}}_s$ . Equation (23) represents the relation of $u^{(0)}$ and function values at its $N_s$ neighboring points, as in the case of the finite difference method (FDM). This is one of the key factors of the proposed LMFS scheme. Note that the above numerical procedures should be repeated for every node inside the computational domain.

Now we can form the final linear system of equations for the LMFS. Suppose the total number of $N=n_i+n_{b1}+n_{b2}$ nodes are arbitrarily distributed inside the computational domain as well as along the boundary, where $n_i$ , $n_{b1}$ and $n_{b2}$ are numbers of interior points, boundary points with Dirichlet and Neumann boundary conditions, respectively. To enforce the satisfaction of Eq. (23) at every interior node will yield the following linear system:

24 $$u^{(i)}-\sum _{j=0}^{N_s}{c}^{(i)(j)}{u}^{(i)(j)}=0,i=1,2,...,n_i.$$

Graph

where $u^{(i)}$ denotes the function value at the ith interior point ${\mathit{x}}^{(i)}$ (central node), $c^{(i)(j)}$ stands for the coefficient at the jth neighboring point around point ${\mathit{x}}^{(i)}$ . Note that in Eq. (24), the relation $u^{(i)}=u^{(i)(\text{0})}$ is hold. For node with Dirichlet boundary conditions, another linear system of equations can be obtained:

25 $$u^{(i)}={\overline{u}}^{(i)},i=n_i+1,n_i+2,...,n_i+n_{b1}.$$

Graph

Similarly, for nodes with Neumann boundary conditions, a third system of linear algebraic equations can be obtained. By combining the above systems of equations, the following spare and banded system of linear algebraic equations can be formed:

26 $$\mathit{Au}\text{=}\mathit{B}\text{,}$$

Graph

where ${\mathit{A}}_{N\times N}$ is the coefficient matrix, $\mathit{u}\text{=}{\left(u^{(1)},u^{(2)},...,u^{(N)}\right)}^T$ and ${\mathit{B}}_{N\times 1}$ is the vector of homogeneous term of the governing equation and the given boundary conditions. On solving this system of equations, numerical solutions of $u(\mathit{x})$ at points scattered inside the whole computational domain and along the boundary can be calculated. The flowchart of the proposed CCS-LMFS procedures for general 3D inhomogeneous elliptic problems is given in Fig. 2. The chart shows the main tasks for the program and the related functions. The source codes writing in MATLAB language can be found in "Appendix A".

Graph: Fig. 2 Flowchart for the developed CCS-LMFS program

In this section, three benchmark numerical examples will be tested to verify the methodologies developed above. The first example is related to the steady-state heat conduction problem [[34]] (Poisson's equation) while the next two are concerned with acoustic [[35]] (Helmholtz equation) and reaction diffusion/wave propagation problems [[36]] (modified Helmholtz equation). The accuracy and sensitivity of the numerical results with respect to the number of local source points ( $M_s$ ), the radius of the artificial sphere ( $R_s$ ), and the degree of Chebyshev's polynomials are carefully investigated. In this study, all computations were implemented on a Core-i7 laptop PC with a 1.99 GHz CPU, 24 GB RAM and 1000 GB hard drive. For the ease of comparison, we considered the following analytical solutions:

27 $$u(x_1,x_2,x_3)=cos(x_1+x_2+x_3)+2$$

Graph

for 3D steady-state heat conduction problem,

28 $$u(x_1,x_2,x_3)=cos(x_1+x_2+k{x}_3)$$

Graph

for 3D acoustic problem, and

29 $$u(x_1,x_2,x_3)=e^{x_1+x_2+k{x}_3}$$

Graph

for 3D reaction diffusion and wave propagation problem, respectively. The following maximum absolute error and the $L_2$ error norm (global error) are adopted:

30 $$E_{\text{max}}=\underset{1\le k\le N_{\mathit{total}}}{\text{max}}\left|I_{\mathit{numer}}^{}(k)-I_{\mathit{exact}}^{}(k)\right|,$$

Graph

31 $$E_{\mathit{global}}={\left[\sum _{k=1}^{N_{\mathit{total}}}{\left[I_{\mathit{numer}}^{}(k)-I_{\mathit{exact}}^{}(k)\right]}^2\right]}^{1/2}/{\left[\sum _{k=1}^{N_{\mathit{total}}}{\left[I_{\mathit{exact}}^{}(k)\right]}^2\right]}^{1/2},$$

Graph

where $I_{\mathit{numer}}^{}$ and $I_{\mathit{exact}}^{}$ stand for the numerical and analytical solutions, respectively, $N_{\mathit{total}}$ is the total number of points tested.

We first consider the following steady-state heat conduction problem [[34]] with Dirichlet boundary conditions in a cubic domain $\mathrm{\Omega }={\left(0,1\right)}^3\subset {\mathit{R}}^{\mathbf{3}}$ :

32 $$\frac{{\partial }^{2}u(x_1,x_2,x_3)}{\partial x_1^2}+\frac{{\partial }^{2}u(x_1,x_2,x_3)}{\partial x_2^2}+\frac{{\partial }^{2}u(x_1,x_2,x_3)}{\partial x_3^2}=f(x_1,x_2,x_3),$$

Graph

where $u(x_1,x_2,x_3)$ represents temperature change, the inhomogeneous term $f(x_1,x_2,x_3)$ and the corresponding boundary conditions are given such that satisfying the given exact solution (27). The problem sketch and the distribution of the Gauss–Lobatto nodes are depicted in Fig. 3a. To illustrate the numerical results, a total number of $N=15,341$ evenly distributed points (including 3174 boundary points and 12,167 interior points) are chosen inside the whole computational domain, as shown in Fig. 3b.

Graph: Fig. 3 Distribution of the Gauss–Lobatto nodes (a) and the LMFS nodes (b)

Firstly, we study the sensitivity of the numerical results with respect to the number of source points ( $M_s$ ) selected inside each of the local subdomain ${\mathrm{\Omega }}_s$ . For this purpose, the radius of the local artificial sphere and the degree of Chebyshev's polynomials are fixed at $R_s=2$ and $N_{\mathit{ch}}=N_1=N_2=N_3=10$ , respectively. As shown in Table 1, the proposed CCS-LMFS results are stable, accurate and rapidly convergent as the number of source points ( $M_s$ ) increases, until the number of sources reaches $M_s=30$ after which a further increase in the source number does not improve substantially the accuracy of the numerical results. It can be also seen from Table 1 that the proposed method can yield more accurate results with $M_s$ increases, but as a price paid for such a high accuracy, the method should take more time for the evaluation.

Maximum and global errors of calculated temperatures and heat fluxes with respect to the number of source points (Ms), with Rs=2 and Nch=10

Next, we fix $M_s=30$ and $N_{\mathit{ch}}=10$ and study the sensitivity of the method with respect to different values of local sphere radius $R_s$ . From the results listed in Table 2 it may be seen that the precision of the proposed CCS-LMFS results increase with the ensuing increase in the size ( $R_s$ ) of the artificial sphere. It can be also observed that the global/maximum errors keep decrease until $R_s$ reaches the value of $R_s=2$ , after which a further increase in $R_s$ does not improve substantially the accuracy of the numerical results. After that, we fix $R_s=2$ and $M_s=30$ and study the sensitivity of the proposed method with respect to the degree of Chebyshev's polynomials ( $N_{\mathit{ch}}$ ). As can be observed from Table 3 that numerical results predicted by using the proposed CCS-LMFS approach converge quickly as the degree of Chebyshev's polynomials increases.

Maximum and global errors of calculated temperatures and heat fluxes as functions of local sphere radius Rs, with Ms=30 and Nch=10

Maximum and global errors of calculated temperatures and heat fluxes with respect to degree of Chebyshev's polynomials N, with Rs=2 and Ms=30

To investigate the convergence rate (order) of the proposed CCS-LMFS approach, the relative error curves (in log–log scale) for both $u(x_1,x_2,x_3)$ and $\partial u(x_1,x_2,x_3)/\partial x_1$ are shown in Fig. 4, as the number of the CCS-LMFS nodes increases from 1000 to 1,000,000 (one million). In our computations, the radius of the local artificial sphere, the degree of Chebyshev's polynomials and the number of source points are taken to be $R_s=1$ , $N_{\mathit{ch}}=6$ and $M_s=12$ , respectively. It can be seen that the proposed CCS-LMFS results for both $u(x_1,x_2,x_3)$ and $\partial u(x_1,x_2,x_3)/\partial x_1$ converge quickly as the number of nodes increases. It is reported that, for the largest model with 1,000,000 points, the proposed CCS-LMFS approach uses only about 136.8 CPU-seconds.

Graph: Fig. 4 Convergence curves of the calculated temperatures and heat fluxes, as functions of various number of LMFS nodes

To study the stability of the proposed method, Fig. 5 gives a comparison of the condition numbers for numerical results using the traditional MFS and the proposed LMFS. In the traditional MFS, the fictitious boundary is taken to be a curve similar to the real boundary of the computational domain. To determine the optimal location of the source points, the leave-one-out cross validation (LOOCV) method proposed in Ref. [[37]] has been used for the MFS solutions. As can be seen from Fig. 5 that the condition number of the traditional MFS increases sharply as the number of MFS nodes increases, which leads the resulting matrix equation to be severely ill-conditioned. In a sharp contrast, the condition number of the proposed LMFS is significantly smaller than that of the classical MFS, which shows significant advantage of the proposed LMFS in the numerical stability compared with the traditional MFS. It is interesting to note that, for the largest model with 3000 nodes, the condition number of the proposed LMFS is only about 130. A self-contained MATLAB code for this example is provided in the part of Supplementary Materials of the paper.

Graph: Fig. 5 Condition numbers for numerical results using the traditional MFS and localized MFS

Next, we consider the acoustic problem [[35]] governed by the following Helmholtz equation:

33 $$Lu(x_1,x_2,x_3)\equiv \left({\mathrm{\nabla }}^2+k^2\right)u(x_1,x_2,x_3)=f(x_1,x_2,x_3),$$

Graph

where $u(x_1,x_2,x_3)$ represents acoustical field, the inhomogeneous term $f(x_1,x_2,x_3)$ and the corresponding boundary conditions are given such that satisfying the given exact solution (28). The problem sketch and the distribution of the Gauss–Lobatto nodes are depicted in Fig. 6. The principal dimensions of this problem are 3.6 m in length, 1.05 m in width and 1.36 m in height. As descripted in Sect. 2.2, the original computational domain should be extended firstly to a cubic so that the Chebyshev collocation method can be used. Here in our computations, the real computational domain is extended into a cubic [− 0.1, 3.7] × [− 0.1, 1.1] × [− 0.1, 1.4], as shown in Fig. 6. This 3D model is subjected to a mixed-type boundary condition, where the sounds are prescribed on the left-half surface $\left\{0\text{m}\le x\le 1.8\text{m}\right\}$ while the normal sound velocities are prescribed on the remaining surface. Here, the frequency of the acoustical field is taken to be $k=\sqrt{3}$ . To solve the problem numerically, a total number of 7307 irregularly distributed LMFS nodes are discretized inside the whole domain (containing 5278 boundary points and 2029 interior points). Here, these nodes are automatically generated by using the CAE software 'Hypermesh'.

Graph: Fig. 6 Geometry of the problem and the distribution of the Gauss–Lobatto nodes

Here, we fix the radius of the fictitious sphere at $R_s=1$ and investigate the sensitivity of the proposed method with respect to the degree of Chebyshev's polynomials ( $N_{\mathit{ch}}$ ) and the number of source points ( $M_s$ ). Figure 7a shows how the global error decays as a function of different degree of Chebyshev's polynomials ( $N_{\mathit{ch}}$ ), with $M_s=30$ . As can be seen in Fig. 7a, the CCS-LMFS results are stable, accurate and rapidly convergent as the degree of Chebyshev's polynomials ( $N_{\mathit{ch}}$ ) increases until the number of $N_{\mathit{ch}}$ reaches $N_{\mathit{ch}}=9$ after which a further increase in the degree of Chebyshev's polynomials does not improve substantially the accuracy of the numerical results. Figure 7b investigates the error variations with respect to the number of source points ( $M_s$ ) selected inside each of the local subdomain, with fixed number of $N_{\mathit{ch}}=9$ . It can be seen that the proposed CCS-LMFS results are rapidly convergent as the number of source points $M_s$ increases. For a Core-i7 1.99 GHz computer, the maximum CPU-time taken to solve above simulations is 13.5 s.

Graph: Fig. 7 Relative error curves for sounds u and sound derivatives ∂u/∂x1, as functions of a different degree of Chebyshev's polynomials Nch and b number of source points Ms

Next, to further illustrate the accuracy and efficiency of the proposed method, we compare the numerical results of the proposed CCS-LMFS approach with those obtained by using other similar works, for example, the BEM, traditional MFS and the generalized finite difference method (GFDM) [[39]]. In the GFDM approach, the second-order Taylor series expansion has been used for deriving the generalized finite difference formulae. The parameters involved in the CCS-LMFS approach are set to be $M_s=20$ , $N_{\mathit{ch}}=6$ and $R_s=1$ . Table 4 shows significant advantage of the proposed CCS-LMFS in the savings of CPU-times compared with the traditional MFS and the BEM.

Comparison of the accuracy (global errors) and efficiency (CPU-times) for sound fields

Finally, we consider the reaction diffusion and wave propagation problem [[36]] governed by the following modified Helmholtz equation:

34 $$Lu(x_1,x_2,x_3)\equiv \left({\mathrm{\nabla }}^2-k^2\right)u(x_1,x_2,x_3)=f(x_1,x_2,x_3).$$

Graph

The problem sketch and the distribution of the LMFS nodes are depicted in Fig. 8. The principal dimension of the problem is 2 m in length (along the x-axis), 1.2 m in width (along the y-axis), and 1 m in height (along the z-axis), as shown in Fig. 8. Here the computational domain is extended to [− 0.1, 2.2] × [− 0.1, 1.4] × [− 0.1, 1.3] so that the Chebyshev collocation method can be used. This 3D model is subjected to a mixed-type boundary condition, where the Dirichlet boundary conditions are prescribed on the right-half surface of the domain $\left\{-1\text{m}\le x\le 0\text{m}\right\}$ , while the Neumann boundary conditions are imposed on the remaining surface. A total number of 19,525 irregularly distributed LMFS nodes are discretized inside the whole domain (containing 16,730 boundary points and 2795 interior points).

Graph: Fig. 8 Geometry of the problem (a) and the nodes distribution of the LMFS model (b)

In the previous examples, we have already investigated the convergence and stability of the proposed method with respect to the number of local source points ( $M_s$ ), the radius of the artificial sphere ( $R_s$ ) as well as the degree of Chebyshev's polynomials ( $N_{\mathit{ch}}$ ). Here, the focus is to investigate the effect of the nondimensional wave number kD on the overall accuracy of the numerical results, where $k$ is the wave number and $D=2$ is the maximum diameter of the computational domain. Here in our computations, the related parameters are taken to be $M_s=30$ , $R_s=1$ and $N_{\mathit{ch}}=10$ .

It can be seen from Fig. 9 that as the values of kD decreases, the numerical accuracy steadily improves. This is because that the larger the value of kD, the more oscillation of the acoustical fields and the more difficult they are going to be solved. For the situation with larger wave number, more collocation points are required in order to obtain accurate solutions. Although not presented, it is reported that numerous other numerical experiments have been performed and the similar conclusions have been drawn. Overall, it can be concluded that the proposed CCS-LMFS scheme could be considered as an alternative for solving numerically the inhomogeneous elliptic problems, especially when high accuracy is required.

Graph: Fig. 9 Relative error curves of u and ∂u/∂x1 with respect to different values of nondimensional wave number kD

This paper documents the first attempt to apply the localized method of fundamental solutions (LMFS) for large-scale simulations of 3D inhomogeneous elliptic boundary value problems. A direct Chebyshev collocation scheme (CCS) is employed for the calculation of the particular solutions of the given inhomogeneous problem. In the CCS approach, the Gauss–Lobatto collocation points are employed to ensure the pseudo-spectral convergence of the Chebyshev interpolation. Once the particular solution has been derived, the original inhomogeneous problem has been transferred into a homogeneous one which is solved by using the LMFS approach. Preliminary numerical experiments, with up to 1,000,000 unknowns, show that the proposed method is accurate and robust for solving 3D inhomogeneous boundary value problems, particularly when high precision is desired.

The work described in this paper was supported by the National Natural Science Foundation of China (Nos. 11872220, 71571108), Projects of International (Regional) Cooperation and Exchanges of NSFC (No. 71611530712), and the Shandong Provincial Natural Science Foundation (Nos. ZR2017JL004, ZR2017BA003, ZR2017MF055).

Graph

Graph

Graph

Graph

Graph

Graph

Graph

Graph

Graph

• $L$
• Linear differential operator
• $\mathrm{\Omega }$
• Original computational domain
• $u_p$
• Particular solution
• $u_h$
• Homogeneous solution
• $\stackrel{⌢}{\mathrm{\Omega }}$
• A cubic which contains $\mathrm{\Omega }$
• $T_n$
• nth Chebyshev polynomial
• N 1, N2, N3
• Degrees of Chebyshev's polynomials

• M

• Number of roots of Chebyshev polynomial
• ${\mathrm{\Omega }}_s$
• Local subdomain
• $N_s$
• Number of collocation points inside ${\mathrm{\Omega }}_s$
• $M_s$
• Number of sources associated with ${\mathrm{\Omega }}_s$
• $R_s$
• Radius of the local artificial sphere
• $G$
• Fundamental solution
• $B(u)$
• Residual function

• $\xi$
• Root of Chebyshev polynomial
• $\varphi$
• Unknown coefficients in Chebyshev method
• $\phi$
• Unknown coefficients in MFS formulation
• $\omega$
• Weighting function

• m, n, l
• Indexes in a sum in Chebyshev method

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.