MATERIAL DISCONTINUITY PROBLEMS SOLVED BY A MESHLESS METHOD BASED ON VARIABLY SCALED DISCONTINUOUS RADIAL FUNCTIONS 1

The paper presents a meshless method based on global interpolation with radial base functions and shows its application in solving interface problems. Such problems arise when two or more diﬀerent materials are used to construct the element under consideration. Across the interface between the materials, a discontinuity of material parameters arises. To solve the problem, the radial basis function-based collocation method is applied. To properly reﬂect the discontinuity, the base functions are modiﬁed. In this paper, the method is applied to solve problems described by elliptic equations. Using these examples, the accuracy, stability and convergence of the method are examined.


Introduction
Many problems in science and engineering are defined for elements composed of two or more different materials.Such problems can be encountered, for example, when considering heat transfer in a rod or plate composed of two different materials in contact to form an interface.Across this interface, a discontinuity of material parameters arises, which can lead to non-smooth or discontinuous solutions (LeVeque and Li, 1994).
Low-order approximation methods such as finite element method or finite difference method can generally be used to solve numerically these interface problems because they can easily approximate a discontinuous derivative of an appropriate order across the interface, leading to a non-smooth or discontinuous solution.Moreover, several more sophisticated low order methods have been developed to handle the interface problems.These include the immersed interface method (LeVeque and Li, 1994;Li, 2003) and discontinuous Galerkin immersed finite element method (Yang and Zhang, 2016) to name a few.
Recently, meshless techniques have been strongly developed (Belytschko et al., 1996;Liu, 2003).They, oppositely to the mentioned methods, need only scattered nodes to discretize a problem and therefore they are more versatile than mesh-based methods.Among the meshless techniques, one can also distinguish the methods that can be successfully employed to solve the interface problems (Yoon and Song, 2014;Stevens and Power, 2015;Martin and Fornberg, 2017).They are based on local approximation schemes.
Although all these numerical methods can solve problems with adequate accuracy, they require imposition of many degrees of freedom to ensure this accuracy.On the other hand, there are methods based on global approximation (Fornberg, 1996;Trefethen, 2000;Krowiak, 2008) -high order approximation methods -which are able to solve some problems with high accuracy using only a few nodes.Among the meshless techniques, one can also distinguish such methods.They usually use global approximation with radial basis functions (RBFs) (Ferreira and Fasshauer, 2007;Chen et al., 2014;Krowiak, 2018).It was found that the RBFs have very interesting features: they are very useful in approximation of scattered data and they are very convenient when solving multidimensional problems (Fasshauer, 2007).Unfortunately, the high order approximation methods are not well suited for problems possessing discontinuities.They are generally not able to catch some local features like, for example, discontinuities across interfaces.
Therefore, the goal of the present paper is to modify a meshless method based on a global approximation with the RBFs to accurately solve the interface problems.In this way, we want to preserve all the advantages of this kind of method, extending the possibilities of its use to problems with discontinuities of material parameters across the interfaces.The well-known Kansa (1990) approach is taken under consideration since it possesses all the advantages of the RBFs and is characterized by a simple formulation.To extend the usefulness of the Kansa method for the interface problems, the so-called variably scaled discontinuous kernel interpolation is introduced into the method.This type of interpolation was recently studied by de Marchi et al. (2020), showing its high accuracy in the approximation of discontinuous functions.

Definition of the interface problem
The examples of one-, two-and three-dimensional physical models composed of different materials that are in contact forming interfaces are presented in Fig. 1.In the figure, Ω + and Ω − denote the areas occupied by materials of different parameters and Γ is the interface.Since the method considered in the present paper is used to discretize the spatial variable, we take an elliptic type equation to show an example of mathematical model of the interface problem.The general form of such an equation is as follows Equation (2.1) with proper boundary conditions can describe boundary or stationary physical problems such as linear elasticity or steady-state heat transfer.In such problems.Eq. (2.1) contains usually the first term called the diffusion one.Then, parameter β(x) corresponds to material parameters of the physical model.Therefore, when the interface problem is considered, this parameter has discontinuity across the interface, which can be put as Moreover, all the parameters in Eq. (2.1) as well as the right hand side function can be discontinuous across the interface.
To find an unique solution of the interface problem one should enrich the mathematical model given by Eq. (2.1) and proper boundary conditions with other equations that describe the continuity or a jump of the solution as well as the derivative of the solution at the interface.These conditions have the following form for In Eq. ( 2.3) u + (x), u − (x) denote the solutions at the interface, when approaching from the Ω + or Ω − subdomain, I + , I − are differential operators associated with appropriate subdomains, imposed on the solution at the interface, and w(x), v(x) are functions that define the jump.
In mechanical engineering problems, the interface conditions usually have a homogeneous form (w(x) = 0, v(x) = 0) describing the continuity of the solution variable u and the continuity of a flux or stress dependently on the problem considered.

Solution method
To solve the problem described in Section 2, a meshless collocation method that uses global interpolation with the RBFs is applied.In the method, the sought solution is interpolated by a series composed of RBFs, which has the following form where N denotes the number of so-called source points ξ j , where the RBFs are centered.In the present paper, the three most often used types of RBFs (Table 1) are considered.

RBFs
Function form Radial functions presented in Table 1 are infinitely differentiable.They possess a parameter ε that controls their flatness influencing the accuracy and stability of the method.Some interesting issues on approximating as well as solving differential equations using the RBFs are presented by Fasshauer (2007).
Global interpolation with infinitely differentiable RBFs does not allow function (3.1) to accurately reflect the discontinuities across interfaces.Therefore, in the present paper, the RBFs are modified according to the idea presented by de Marchi et al. (2020) for kernel functions.Following this idea, the RBFs in the present work are modified by introducing a scaling function ψ(x).
Assuming that the domain of interest belongs to d-dimensional space, i.e. x ∈ Ω ⊆ R d , then Ψ : R d → R and the variably scaled RBFs are defined as where φ : R d+1 × R d+1 → R and φ ψ : In order to introduce a discontinuity to the variably scaled RBFs, the scaling function should be a piecewise continuous one, possessing discontinuity across the interface.It is sufficient to introduce this function as a piecewise constant one, as follows Taking into account the definition of variably scaled RBFs (3.2) and scaling function (3.3), we can conclude that if x and ξ j come from the same subdomain one obtains the original RBF, i.e. φ( x − ξ j 2 ) but if they are taken from different subdomains, we have φ x − ξ j 2 + (α − β) 2 .Now, the interpolant written in terms of the modified RBFs approximates the sought solution of the interface problem, what can be put as To find the solution, the domain of the problem is discretized with the source points of RBFs ξ j , j = 1, . . ., N and with the collocation points inside the domain x I i , i = 1, . . ., N I , on the boundaries x B i , i = 1, . . ., N B and at the interface x Γ i , i = 1, . . ., N Γ .Then, the function is introduced into the mathematical model of the interface problem.By collocating the equations at appropriate points, the algebraic system of equations is obtained in the following form For simplicity of the presentation, the differential operator from the governing equation is denoted by L and the one from boundary conditions is denoted by B. To properly reflect the interface conditions given by the last two equations of system (3.5), the following explanation on how to evaluate function φ + ψ and φ − ψ at x Γ i should be made.In this case, taking into account the definition of scaling function (3.3), we have Similarly, the derivative of this function at interface nodes is evaluated.Equations (3.5) can be written in a convenient matrix notation Γ , Φ Γ ] T , b = [f , g, w, v] T , and the vector c contains all the interpolation coefficients c j , j = 1, . . ., N .
The elements of the submatrices and subvectors contained in Eq. (3.7) are computed from the formulas for j = 1, . . ., N (Φ L ) ij = Lφ (x, ψ(x)), (ξ j , ψ(ξ j )) Γ ) ij = I − φ (x, α), (ξ j , ψ(ξ j )) − I + φ (x, β), (ξ j , ψ(ξ j )) and from the formulas In this paper, the discretization which leads to the square system of equations is considered.In this case, the numbers of discretization points have to fulfil the following condition: ) is conditioned by invertibility of Φ matrix.Some notes on the invertibility of this type of matrix can be found in the work of Hon and Schaback (2001).From this information, it follows that although the invertibility is not guaranteed, the cases where the matrix is singular are very rare and the method has been successfully applied to many problems in science and engineering (Chen et al., 2014).To guarantee the solvability of the problem, the so-called symmetric collocation approach should be used (Wu, 1992).
The matrices resulting from global interpolation using the RBFs are dense and prone to illconditioning.Therefore, some approaches to control this numerical phenomenon with a proper value of the shape parameter of RBFs have to be included during the calculation process.Among the algorithms designed for this purpose, it was found in the present work that a kind of cross--validation approach presented by Rippa (1999) yields good results.This algorithm is used in finding the value of the shape parameter of RBFs, when solving examples in the present paper.

Numerical examples and discussion
To examine the method, several examples from mechanical engineering and science have been solved.In the present paper, a few of them are presented.The first example presents the interface problem that can describe one-dimensional elasticity.The governing equation and boundary conditions are as follows Due to the use of different materials, the Young modulus is discontinuous across the interface at x = 5 In the calculations the following values are assumed: E + = 10 4 , E − = 10 3 and the right hand side function describing the body forces is as follows b(x) = 25x − 7.5x 2 + 0.5x 3 .In the considered problem, the interface conditions define the continuity of displacement and the continuity of stress at the interface point.They have the form According to Eq. ( 4.3), we expect a continuous but non-smooth solution and discontinuous derivative of the solution.The scaling function is assumed as follows In Fig. 2, the results of application of the method with conventional RBFs are shown.To obtain the results, the multiquadric RBFs have been used and the domain has been discretized by N = 40 equispaced source points.The method that uses classical RBFs does not catch the interface features.
The same problem has been solved by the method with variably scaled discontinuous RBFs using the same discretization parameters.The obtained results are shown in Fig. 3.
The method introduced in Section 3 accurately reflects the non-smooth solution and the discontinuity of its derivative, as can be seen in Fig. 3, where the approximate solution as well as its derivative coincide with the exact one.
In this paper, three types of RBFs listed in Table 1 have been tested for accuracy.For each type of RBFs and each type of point distribution, the shape parameter ε of RBFs has been determined.The results in the form of root-mean-square error between the approximate solution and the exact one are included in Table 2.
The presented results show that the method with variably scaled discontinuous RBFs gives very accurate results using a relatively small number of discretization points.As the second example, a problem modelled by a partial differential equation is considered.It can describe a steady-state heat transfer in a square domain with s circle interface presented in Fig. 4. The mathematical description of the problem is as follows where diffusion parameter β(x) has a jump across the interface due to the use of different materials The interface conditions have the form for where n is the normal outward vector to the interface.
In the calculation, the following values and functions have been assumed in the model To reflect the discontinuity, the scaling function has been taken as The solution obtained using classical RBFs is presented in Fig. 5a along with the exact one (Fig. 5b).The solution obtained using variably scaled discontinuous RBFs and the absolute error δ = |u h − u| are shown in Fig. 6.To obtain the results, Gaussian RBFs have been applied with following discretization parameters: N = 400, N I = 196, N B = 60, N Γ = 72.A pseudorandom distribution of the RBF source points (Fig. 4b) and the uniform one for the collocation points (Fig. 4c) have been applied in the calculation.
For a detailed assessment, the profile of the solution u h (x, y = 0) and the profile of the derivative of the solution u ′ hx (x, y = 0) along with the exact ones are shown in Fig. 7. Various types of RBFs have been tested in the present work and several point distributions.In all cases presented in this paper, the RBF source points were distributed pseudorandomly but the collocation points were evenly distributed.The detailed results are shown in Table 3.
The results show that the use of variably scaled discontinuous RBFs allows the method to accurately reflect a non-smooth solution unlike the classical method.All types of RBFs lead to similar accuracy if only the shape parameter is appropriately determined.

Conclusion
The conventional RBF collocation method has many advantages: rapid convergence, convenient use in high dimensional problems and simple implementation but it is not able to accurately catch some local features.In the present work, the method is extended to be useful in the interface problems.To this end, the scaling function is introduced into the RBFs.A piecewise constant scaling function is sufficient to make the RBFs discontinuous across the interface.The use of such RBFs in the conventional collocation procedure of Kansa method allows the method to accurately solve the interface problems while retaining all the other advantages of this method.The results show that various types of RBFs can be used giving similar accuracy, if only appropriate values of the shape parameter are determined.
In the future work, it is planned to find whether the values assumed in the scaling function affect the quality of the results and to find a method of determining these values, which lead to the best accuracy.

Fig. 1 .
Fig. 1.Examples of physical models with interfaces Solution of Eq. (3.7) gives the values of the interpolation coefficients, i.e. c = Φ −1 b (3.8) Finally the approximate, analytic solution of the problem, in the form of interpolation function (3.4) is obtained.Unique solution (3.8

Fig. 2 .
Fig. 2. The solution (a) and derivative of the solution (b) obtained by the classical RBF method; solid line -approximate solution, dash line -exact one

Fig. 3 .
Fig. 3.The solution (a) and derivative of the solution (b) obtained by the modified multiquadric RBFs; solid line -approximate solution, dash line -exact one

Fig. 4 .
Fig. 4. The domain of the problem (a) with examples of the RBF source point distribution (b) and collocation point distribution (c)

Fig. 5 .
Fig. 5.The solution obtained using classical Gaussian RBFs (a) and the exact one (b)

Fig. 6 .Fig. 7 .
Fig. 6.The solution obtained with the modified Gaussian RBFs (a) and the error of the solution (b)

Table 2 .
The accuracy (RMS error) of making use of various types of RBFs

Table 3 .
The accuracy (RMS error) of making use of various types of RBFs