A Simple, Accurate and Semi-Analytical Meshless Method for Solving Laplace and Helmholtz Equations in Complex Two-Dimensional Geometries

: A localized virtual boundary element–meshless collocation method (LVBE-MCM) is proposed to solve Laplace and Helmholtz equations in complex two-dimensional (2D) geometries. “Localized” refers to employing the moving least square method to locally approximate the physical quantities of the computational domain after introducing the traditional virtual boundary element method. The LVBE-MCM is a semi-analytical and domain-type meshless collocation method that is based on the fundamental solution of the governing equation, which is different from the traditional virtual boundary element method. When it comes to 2D problems, the LVBE-MCM only needs to calculate the numerical integration on the circular virtual boundary. It avoids the evaluation of singular/strong singular/hypersingular integrals seen in the boundary element method. Compared to the difﬁculty of selecting the virtual boundary and evaluating singular integrals, the LVBE-MCM is simple and straightforward. Numerical experiments, including irregular and doubly connected domains, demonstrate that the LVBE-MCM is accurate, stable, and convergent for solving both Laplace and Helmholtz equations.


Introduction
The boundary element method (BEM) [1,2] is a well-known numerical method that has become an alternative to domain methods such as the finite element method (FEM) [3,4] for the simulation of certain physical problems. The core of this method is to accurately solve singular integrals, especially nearly singular, strongly singular, and hypersingular integrals, among others. Substantial efforts have been devoted to developing and applying efficient estimation techniques for such integrals. Lutz [5] proposed a special Gaussian-type numerical integral to calculate the singular and nearly singular integrals. Johnston and Elliott [6] proposed a sinh transformation to evaluate nearly singular integrals. Niu and Zhou [7,8] suggested that the asymptotic expansion of the kernel function with respect to the local co-ordinates should be employed to address singular integrals. Besides these methods, there are other techniques that can be used to deal with various singular integrals. Although they were proven to be the effective strategies, these methods are often time consuming, tedious, and expensive.
In recent years, the virtual boundary element method (VBEM) [9 -13] has been proposed to overcome the above shortcomings. The VBEM introduces the virtual boundary to avoid the calculation of singular integrals and inherits the semi-analytical and highaccuracy features of the BEM. Sun and Yao [14] used the VBEM to successfully solve thin plate elastic obstacle problems. Yao et al. [10,15] used the VBEM to simulate magnetoelectroelastic and piezoelectric problems. Yang et al. [16] and Liu et al. [17] resolved three-dimensional inverse heat conduction problems using the VBEM. As a boundarytype scheme with global discretization, however, the VBEM encounters challenges when simulating large-scale and/or high-dimensional problems.
More recently, the localization of boundary-type meshless methods have received considerable attention, and various localized meshless methods [18][19][20][21][22][23] have been proposed to solve mathematical and mechanical problems, such as the generalized finite difference method (GFDM) [24,25], the localized method of fundamental solutions (LMFS) [26,27], the local knot method (LKM) [28,29], and the localized singular boundary method (LSBM) [30,31]. Unlike traditional boundary-type methods, these methods are not only simple, accurate, and easy-to-program, but also suitable for large-scale simulations in complicated domains. On the other hand, boundary-type meshless methods encounter many difficult issues. Similar to the fundamental solution method, the VBEM uses fundamental solutions as the basis functions and requires a virtual boundary outside of the physical domain to avoid source singularity. The selection of this artificial boundary is still a well-known tricky issue in spite of the great deal of effort that has been made to address this problem [32,33], especially in terms of complex geometries.
Motivated by the above works using localized methods, we establish a localized numerical framework for the VBEM in this paper, which we called the localized virtual boundary element-meshless collocation method (LVBE-MCM). The accuracy and effectiveness of the LVBE-MCM was verified via its ability to solve Laplace and Helmholtz equations in complex 2D domains. In the traditional VBEM, it is difficult to determine the position and shape of virtual boundaries in the complex domain because these boundaries have a certain impact on the calculation accuracy. On the contrary, the LVBE-MCM only uses the circular virtual boundary during the local approximation, and it is insensitive to the location of the boundary. Furthermore, the resulting LVBE-MCM system is sparse and can thus be easily solved using an ordinary computer. This also means that the method has certain application prospects for solving large-scale problems.
The rest of the paper is organized as follows: In Section 2, the considered problem is briefly introduced. Section 3 describes the detailed numerical procedure for the LVBE-MCM. Section 4 develops an augmented moving least squares approximation using the fundamental solutions. In Section 5, two numerical examples are provided to confirm the effectiveness and applicability of the proposed method. The conclusions are summarized in Section 6.

Preliminaries
Let Ω ∈ R 2 be an open bounded domain surrounded by the boundary Γ = ∂Ω, which is assumed to be piecewise smooth, and consider the following boundary value problem: where L is the Laplace (L = ∇ 2 ) or Helmholtz (L = ∇ 2 + λ 2 ) operator, λ is the wave number, n is the unit outward normal vector, α and β are constants, and f (x), g(x), h(x) are the provided smooth functions on the boundaries. Here, Γ D , Γ N , and Γ R represent the Dirichlet, Neumann, and Robin boundaries, respectively. The fundamental solutions for the Laplace and Helmholtz operators are determined by [34] u * (r) = − 1 2π ln(r), for Laplace operator (5) 0 (λr), for Helmholtz operator (6) where r denotes the Euclidean distance between the field point and the source point, and 0 is a zero-order Hankel function of the first kind.

Localized Virtual Boundary Element-Meshless Collocation Method
First of all, the N = n i + n b1 + n b2 + n b3 discrete nodes x (i) , i = 1, 2, . . . , N are placed over the computational domain Ω, where n i is the number of nodes inside the domain, and n b1 , n b2 , and n b3 indicate the number of nodes along the Dirichlet, Neumann, and Robin boundary, respectively. Considering an arbitrary node x (i) , which is also known as the central node, its m supporting nodes x (i) j , j = 1, 2, . . . , m can be determined based on the nearest nodes. At the same time, the local supporting domain Ω s covering m+1 nodes can also be determined, and its virtual boundary Γ (i) can be specified at a certain distance from the boundary of the supporting domain. For 2D problems, this boundary is a circle. Figure 1 shows the schematic diagram of the LVBE-MCM.
x are the provided smooth functions on the boundaries. Here, D Γ , N Γ , and R Γ represent the Dirichlet, Neumann, and Robin boundaries, respectively.
The fundamental solutions for the Laplace and Helmholtz operators are determined by [34] ( ) for Helmholtz operator (6) where r denotes the Euclidean distance between the field point and the source point, H is a zero-order Hankel function of the first kind.

Localized Virtual Boundary Element-Meshless Collocation Method
First of all, the  In the present study, the virtual boundary is discretized by M exact geometrical elements, and the physical quantity is approximated by the constant element. Using x , the unknowns at nodes ( ) , 0,1,..., Γ . Equation (7) can be rewritten in the following matrix form: In the present study, the virtual boundary is discretized by M exact geometrical elements, and the physical quantity is approximated by the constant element. Using x where is the fundamental solution of the governing equation, and ϕ(ξ) is the distribution density function associated with the virtual boundary Γ (i) . Equation (7) can be rewritten in the following matrix form: For the above local approximation, the moving least squares (MLS) can be employed to obtain the unknown coefficient vector Replacing x (i) j in Equation (7) with x (i) , the following formula is yielded: Then, substituting Equation (9) into Equation (10), we obtain If x (i) is a node on the boundary, the normal derivative can be calculated by where N (i) = n 1 σ In the above equations, n 1 , . . . , n d denote the components of the vector n, and x denote the coordinate components of the node x (i) . Taking all of the nodes x (i) , i = 1, 2, . . . , N and the boundary data provided in Equations (2)-(4) into account, the following overdetermined equations can be obtained: where u = u(x (1) ), u(x (2) ), . . . , u(x (N) ) T , b N×1 is a vector composed of zero elements and boundary data, and A N×N is a sparse matrix. Equation (14) is a well-conditioned system, and in this work, it is solved by MATLAB routine "A\b".

Augmented Moving Least Squares Approximation
The moving least squares approximation is a widely used technique in various meshless/meshfree methods. In this study, the fundamental solutions are introduced into the traditional moving least squares method, and we then developed the augmented moving least squares approximation, which is similar to the one outlined in [35]. According to its basic idea, the vector α (i) is deduced by minimizing the following functional equation: where in Equation (16), d j = x Hence, we have Mathematics 2022, 10, 833 5 of 9 By calculating and reorganizing Equation (17), we can obtain a system equation in matrix form: where It should be pointed out that we used MATLAB's mldivide (matrix left divide) function (P (i) \Q (i) ) to obtain H (i) instead of the matrix inversion.

Numerical Examples
Two numerical examples are provided to demonstrate the effectiveness and accuracy of the proposed method. To evaluate the numerical errors, we adopt the maximum absolute error (MAE) and the root-mean-square error (RMSE), which are defined as follows: where u n and u e represent the numerical and analytical solution at node x j , respectively. All computations were performed using MATLAB 2018b on a desktop PC (Intel ® Core TMi7-6700 CPU at 3.4 GHz, 16G RAM, and Hard Disk-500G).

Example 1. Consider a Laplace equation on an irregular domain with mixed boundary conditions.
The geometry and boundary conditions are shown in Figure 2. For the Robin boundary condition, α = 1 and β = 5. The analytical solution is obtained by u(x 1 , x 2 ) = cos(x 1 ) cosh(x 2 ) + sin(x 1 )sinh(x 2 ) + e x 1 cos(x 2 ) + e x 2 sin(x 1 ) + x 1 2 − x 2 2 + 2x 1 + 3x 2 + 1. First of all, N = 6784 nodes are chosen, and 8 Gaussian points are used. It can be seen from Figure 3 that the numerical error first decreases and then increases as the number of virtual elements increases, meaning that high computational accuracy has been achieved. Then, M = 15 is fixed. From Figure 4, we can observe that the number of Gaussian points has little effect on the calculation accuracy, and therefore, fewer Gaussian points can be used in the calculation.  First of all, N = 6784 nodes are chosen, and 8 Gaussian points are used. It can be seen from Figure 3 that the numerical error first decreases and then increases as the number of virtual elements increases, meaning that high computational accuracy has been achieved. Then, M = 15 is fixed. From Figure 4, we can observe that the number of Gaussian points has little effect on the calculation accuracy, and therefore, fewer Gaussian points can be used in the calculation.   The LMFS and GFDM are recently developed meshless approaches that are very similar to the present LVBE-MCM. In Table 1, these two methods are compared to the proposed approach. It can be observed that all methods are convergent. Although the LMFS and LVBE-MCM have similar numerical accuracy, the latter is slightly better than the former.   The LMFS and GFDM are recently developed meshless approaches that are very similar to the present LVBE-MCM. In Table 1, these two methods are compared to the proposed approach. It can be observed that all methods are convergent. Although the LMFS and LVBE-MCM have similar numerical accuracy, the latter is slightly better than the former.  Figure 5) is considered. The boundary conditions are specified by the analytical solution u(x 1 , x 2 ) = cos(x 1 /2+ √ 3x 2 /2) with λ = 1.   Figure 6. The maximum absolute error and the root-mean-square error are 3.3238 × 10 -8 and 6.4048 × 10 -9 , respectively. This indicates the high-accuracy of the proposed method. Furthermore, it can be observed from Table 2 that the LVBE-MCM has higher numerical accuracy than the LMFS when the same number of sources and elements is adopted.  Figure 6. The maximum absolute error and the root-mean-square error are 3.3238×10 -8 and 6.4048×10 -9 , respectively. This indicates the high-accuracy of the proposed method. Furthermore, it can be observed from Table 2 that the LVBE-MCM has higher numerical accuracy than the LMFS when the same number of sources and elements is adopted.

Conclusions
The localized virtual boundary element-meshless collocation method (LVBE-MCM) was proposed as a novel domain-type meshless method that could be used to solve Laplace and Helmholtz equations in complex 2D geometries. In this work, the traditional virtual boundary element method with a global approximation was modified to a local approximation approach by introducing the moving least square method and local approximation theory. Numerical integrations are only required on the circular virtual

Conclusions
The localized virtual boundary element-meshless collocation method (LVBE-MCM) was proposed as a novel domain-type meshless method that could be used to solve Laplace and Helmholtz equations in complex 2D geometries. In this work, the traditional virtual boundary element method with a global approximation was modified to a local approximation approach by introducing the moving least square method and local approximation theory. Numerical integrations are only required on the circular virtual boundary; thus, the exact geometry elements are convenient to use. The proposed LVBE-MCM avoids the need to evaluate the singular/strong singular/hypersingular integral in the boundary element method and has a higher calculation accuracy than the LMFS.
Two examples involving irregular geometries and doubly connected domains were investigated in detail. The numerical results indicate that the LVBE-MCM is accurate and effective for solving Laplace and Helmholtz equations in complex two-dimensional geometries. The number of Gaussian points has a little effect on the calculation accuracy, and therefore, fewer Gaussian points can be used in the calculation. Moreover, the scheme is convergent with respect to increasing the number of total nodes. It is worth noting that the proposed method can be directly extended to other partial differential equations with known fundamental solutions, such as diffusion equations, Stokes equations, and biharmonic equations.
It should also be pointed out that this paper investigates the accuracy and convergence of the LVBE-MCM numerically. Unlike the difference method and Taylor expansion, it is not an easy work to formally prove the convergence and stability of the LVBE-MCM since there are few related assumptions and theorems on approximation techniques that use the fundamental solution and the augmented moving least squares scheme. Consequently, a theoretical analysis of the LVBE-MCM will be the key issue in our subsequent work.