Refined Meshless Local Strong Form solution of
CauchyNavier equation on an irregular domain
Abstract
This paper considers a numerical solution of a linear elasticity problem, namely the CauchyNavier equation, using a strong form method based on a local Weighted Least Squares (WLS) approximation. The main advantage of the employed numerical approach, also referred to as a Meshless Local Strong Form method, is its generality in terms of approximation setup and positions of computational nodes. In this paper, flexibility regarding the nodal position is demonstrated through two numerical examples, i.e. a drilled cantilever beam, where an irregular domain is treated with a relatively simple nodal positioning algorithm, and a Hertzian contact problem, where again, a relatively simple hrefinement algorithm is used to extensively refine discretization under the contact area. The results are presented in terms of accuracy and convergence rates, using different approximations and refinement setups, namely Gaussian and monomial based approximations, and a comparison of execution time for each block of the solution procedure.
keywords:
Meshless Local Strong Form Method, Weighted Least Squares, Shape functions, CauchyNavier equation, Cantilever beam, Hertzian contact, hrefinement, Irregular domain1 Introduction
Linear elasticity problems, governed by the CauchyNavier equation, are typically addressed in their weak form with the Finite Elements Method (FEM) fem_solid (). However, the problem has also been addressed in its strong form in the past, e.g. componentwise iterative solution with the Finite Differences Method (FDM) volume_based_finite_difference () and with the Finite Volumes Method (FVM) solid_FVM1 (). Besides mesh based methods, meshless methods have also been employed for solving solid mechanics problems in strong and weak form meshless_solid (); strong_form_solid (). The conceptual difference between meshless methods and mesh based methods is in the treatment of relations between nodes. In the mesh based methods the nodes need to be structured into polygons (mesh) that cover the whole computational domain, while on the other hand, meshless methods fully define relations between nodes through the relative inter nodal positions trobec_kosec (), with an immediate consequence of greater generality of the meshless methods.
Strong form meshless methods can be understood as generalizations of FDM, where instead of predetermined interpolation over a local support, a more general approach with variable support and basis is used to evaluate partial differential operators Kouszeg1 (), e.g. collocation using Radial Basis Functions strong_form_solid (); Kouszeg_PN () or approximation with monomial basis Prax (). There are many other methods with more or less similar methodology introducing new variants of the strong form meshless principle meshless_review (). On the other hand, weak form meshless methods are generalizations of FEM. Probably the most known method among weak form meshless methods is the Meshless Local Petrov Galerkin Method (MLPG) Atluri (), where for each integration point a local support is used to evaluate field values and weight functions of a Moving Least Squares (MLS) approximation are used as test functions. In last few decades there have been many variants of MLPG introduced to mitigate numerical instabilities and to improve accuracy and convergence rate, etc. meshless_review ().
In general, recent developments in meshless community are vivid, ranging from analyses of computer execution on different platforms trobec_kosec (); meshless_gpu (), reducing computational cost by introducing a piecewise approximation li_piecewise_2017 () to implementation of more complex multiphase flow kovarik_meshless_2016 (), and many more.
This paper extends the spectra of published papers with a generalized formulation of a local strong form meshless method, termed Meshless Local Strong Form Method (MLSM) enriched with hrefinement Kouszeg_cmc () and ability to discretize arbitrary domains Kouszeg1 ().
The introduced meshless approach is demonstrated on a solution of a benchmark cantilever beam case slaughter2012linearized () and a Hertzian contact problem pereira2016convergence (). The results are presented in terms of displacement and stress plots, comparison against closed form solutions, convergence analyses, and execution time analyses.
The goal of this paper is to demonstrate generality of MLSM that is driven by the fact that all the building blocks of the method depend only on the relative positions between the computational nodes. This is a very useful feature, especially when dealing with problems in multidimensional spaces, complex geometries, and moving boundaries. This feature can be also exploited to write elegant generic code utils_web ().
The rest of the paper is organized as follows: in section 2 the MLSM principle is explained, in section 3 the governing problem is introduced, section 4 is focused on solution procedure, section 5 focuses on discussing the results, and finally, the paper offers some conclusions and guidelines for future work in the last section.
2 MLSM formulation
The core of the spatial discretization used in this paper is a local approximation of a considered field over the overlapping local support domains, i.e. in each node we use approximation over a small local subset of neighbouring nodes. The trial function is thus introduced as
(1) 
with , , and standing for the number of basis functions, approximation coefficients, basis functions and the position vector, respectively. In cases when the number of basis functions and the number of nodes in the support domain are the same, , the determination of coefficients simplifies to solving a system of linear equations, resulting from evaluating equation (1) in each support node and setting it equal to a true value , for from 1 to :
(2) 
where are positions of support nodes and is the actual value of considered field in the support node . The above system can be written in matrix form as
(3) 
where stands for coefficient matrix with elements . The most known method that uses such an approach is the Local Radial Basis Function Collocation Method (LRBFCM) that has been recently used in various problems strong_form_solid (); Kouszeg_PN ().
In cases when the number of support nodes is higher than the number of basis functions () a WLS approximation is chosen as a solution of equation (3), which becomes an overdetermined problem. An example of this approach is DAM Prax () that was originally formulated to solve fluid flow in porous media. DAM uses six monomials for basis and nine noded support domains to evaluate first and second derivatives of physical fields required to solve the problem at hand, namely the Navier Stokes equation. To determine the approximation coefficients , a norm
(4) 
is minimized, where is a diagonal matrix with elements with
(5) 
where stands for weight shape parameter, for centre of support domain and for the distance to the nearest support domain node. There are different computational approaches to minimizing (4). The most intuitive and also computationally effective is to simply compute the gradient of with respect to and setting it to zero, resulting in a positive definite system
(6) 
The problem of this approach is bad conditioning, as the condition number of is the square of the condition number of , unnecessarily increasing numerical instability. A more stable and more expensive approach is QR decomposition. An even more stable approach is SVD decomposition, which is of course even more expensive. Nevertheless, the solution of equation (6) can be written generally in matrix form as
(7) 
where stands for a Moore–Penrose pseudo inverse of matrix . By explicitly inserting equation (7) for into (1), the equation
(8) 
is obtained, where is called a shape function. Now, we can apply a partial differential operator to the trial function, and get
(9) 
In this paper we deal with a CauchyNavier equation and therefore following shape functions are needed, expressed explicitly as
(10)  
(11)  
(12)  
(13)  
(14) 
The shape functions depend only on the numerical setup, namely nodal distribution, shape parameter, basis and support selection, and can as such be precomputed for a specific computation.
3 Governing problem
The goal in this paper is to numerically determine the stress and displacement distributions in a solid body subjected to the applied external force. To obtain a displacement vector field throughout the domain, a CauchyNavier equation is solved, which can expressed concisely in vector form as
(15) 
where and stand for Lamé constants. In two dimensions we express and the equation reads
(16)  
(17) 
Two types of boundary conditions are commonly used when solving these types of problems, namely essential boundary conditions and traction (also called natural) boundary conditions. Essential boundary conditions specify displacements on some portion of the boundary of the domain, i.e. , while traction boundary conditions specify surface traction , where is an outside unit normal to the boundary of the domain and
(18) 
is the stress tensor. In terms of displacement vector the traction boundary conditions read
(19)  
(20) 
where and denote the Cartesian components of and .
4 Solution procedure
4.1 Discretization of the problem
The elliptic boundary value problem at hand is discretized into a linear system of algebraic equations by approximating the differential operations using MLSM, as described in section 2. A block system of linear equations for two vectors and of unknowns representing values and , respectively, is constructed. This system is a discrete analogy of PDE (15) and can symbolically be represented as
(21) 
where and stand for unknown displacements, and for values of boundary conditions and blocks contain precomputed shape functions (10–14). With standing for a list of indices of the chosen neighbours of a point , as introduced in the beginning of section 2, we can, for all indices of internal nodes, express
(22) 
(23) 
for each . Note that equation (22) represents direct discrete analogue of (16) and, likewise, (23) of (17).
Similarly, for all indices of boundary nodes with traction boundary conditions we express
(24) 
(25) 
for each , where are the Cartesian components of the outside unit normal to the boundary in node . Again, equation (24) is a direct analogue of (19) and (25) of (20). And finally, for indices of nodes with essential boundary condition, we express
(26) 
4.2 Positioning of nodes in a complex domain
Meshless methods are advertised as the methods that do not require any topological relations among nodes. That implies that even randomly distributed nodes could be used random_mesh (). However, it is wellknown that with regularly distributed nodes one achieves much better results in terms of accuracy and stability amani (). This has also been recently reported for MLSM in a solution of a NavierStokes problem Kouszeg1 (). The reason behind the sensitivity regarding the distribution of nodes lies in the generation of shape functions. To construct a stable method well balanced support domains are needed, i.e. the nodes in support domain need to be distributed evenly enough Kouszeg1 (). This condition is obviously fulfilled in regular nodal distributions, but when working with more interesting geometries, the positioning of nodes requires additional treatment. In literature one can find several algorithms for distributing the nodes within the domain of different shapes onate (); bubble (). In this paper we will use an extremely simple algorithm, introduced in Kouszeg1 () to minimize the variations in distances between nodes in the support domain. The basic idea is to “relax” the nodes based on a potential between them. Since a Gaussian function is a suitable potential and already used as weight in the shape functions, the nodes are translated simply as
(27) 
where , , and stand for the translation step of the node, position of th support node, relaxation parameter and number of support nodes, respectively (Figure 1(a)). After offsets in all nodes are computed, the nodes are repositioned as
(28) 
Presented iterative process procedure begins by positioning the boundary nodes, which are considered as the definition of the domain and are kept fixed throughout the process.
4.3 hrefinement
Besides flexibility regarding the shape of the domain, nodal refinement is often mandatory to achieve desired accuracy in cases with pronounced differences in stress within the domain. A typical example of such situation is a contact problem pereira2016convergence (). To mitigate the error in areas with high stresses the hrefinement scheme is used., which has already been introduced into different meshless solutions libre (); Bourantas (). In context of local RBF approximation the hrefinement has been used in the solution of the Burger’s equation Kouszeg_cmc (), where a quadtree based algorithm has been used to add and remove child nodes symmetrically around the parent node in transient solution of Burgers’ equation. However, the algorithm presented in Kouszeg_cmc () supported only regular nodal distribution. In this paper we generalize it also to irregular nodal distribution.
In each node to be refined, new nodes are added on the half distances between the node itself and its support nodes
(29) 
where index indicates th support node. When adding new nodes, checks are performed if the newly added node is too close to any of the existing nodes; in that case the node is not added. Moreover, if the refined node and support node are both boundary nodes, newly added node is positioned on the boundary (Figure 1(b)). This procedure can be repeated several times if an even more refined domain is desired. These subsequent refinements will be called levels of refinement and will be denoted as level for the refinement that resulted from applications of the described algorithm.
The described algorithm follows the concept of meshless methods and as such does not require any special topological relations between nodes to refine a certain part of the computation domain. It is also flexible regarding the dimensionality of the domain, i.e. there is no difference in implementation of 2D or 3D variant of the algorithm.
An example on a nontrivial refinement is demonstrated in Fig. 3, where a domain with a hole is considered. The vicinity of the hole is four times refined and then, to mitigate possible irregularities during refinement, relaxed.
4.4 Asymptotic complexity of MLSM
The asymptotic complexity analysis begins with an assumption that evaluations of basis functions, weights, linear operators and boundary conditions take time. For simple domain discretization, such as the uniform grid in a rectangle or random positioning, time is required, where stands for number of computational nodes. To find the neighbours of each point, a tree based data structure such as kdtree trobec_kosec (), taking time to construct and time to query closest nodes, is used. The relaxation of nodal positions (see section 4.2) with iterations costs additional time. Refinding the support nodes by rebuilding the tree and querying for support nodes once again, requires another time. Calculation of the shape functions requires SVD decompositions, each taking time, as well as some matrix and vector multiplication of lower complexity. Assembling the matrix takes time and assembling the right hand side takes of time. Then, the system is solved using BiCGSTAB iterative algorithm. The final time complexity is thus , where stands for the time spent by BiCGSTAB.
For comparison, the complexity of a wellknown weak form Element Free Galerking method (EFG) EFG () differs from MLSM in construction of the shape functions, whose computation requires time using EFG method, with standing for the number of Gauss integration points per node. Additionally, the number of nonzero elements in the final system of EFG is of order times higher than that of MLSM, again increasing the complexity of EFG.
5 Numerical examples
5.1 Cantilever beam
First, the standard cantilever beam test is solved to assess accuracy and stability of the method. Consider an ideal thin cantilever beam of length and height covering the area . Timoshenko beam theory offers a closed form solution for displacements and stresses in such a beam under plane stress conditions and a parabolic load on the left side. The solution is widely known and derived in e.g. slaughter2012linearized (), giving stresses in the beam as
(30) 
and displacements as
(31)  
where is the moment of inertia around the horizontal axis, is Young’s modulus, is the Poisson’s ratio and is the total load force.
In the numerical solution, traction free boundary conditions are used on the top and bottom of the domain, essential boundary conditions given by (31) are used on the right and traction boundary conditions given by (30) on the left
(32)  
(33) 
(34)  
(35)  
(36)  
(37)  
(38)  
(39) 
The problem is solved using MLSM method with or support nodes and Gaussian weight with (see (5)). Two sets of basis functions are considered, 9 monomials
(40) 
and 9 Gaussian basis functions (see (5) for definition) centred in support nodes. In the following discussions these two choices of basis functions will be referred to as M9 and G9, respectively.
System (21) is solved with BiCGSTAB iterative algorithm van1992bi () with ILUT preconditioning saad1994ilut (). Values of , , , and were chosen as physical parameters of the problem.
The acquired numerical solution of the cantilever beam problem is shown in Figure 4.
The error of the numerical approximation of stresses and displacements is measured in relative discrete norm, using
(41)  
(42) 
as error indicators, with representing the set of all nodes. Convergence with respect to the number of computational nodes is shown in Figure 5. The numerical approximations converge towards the correct solution in stress norm as well, with approximately the same convergence rate.
It can be seen that monomials converge very regularly with order 1 as expected, while Gaussian functions exhibit slightly worse convergence. Such behaviour has already been reported in solution of diffusion equation, where MLSM with Gaussian basis failed to obtain accurate solution with a high number of computational nodes. More details about the phenomenon and further reading can be found in KosecZinterhof ().
The method was compared to the standard Element Free Galerkin (EFG) method MMrev (). The EFG method used circular domains of influence with radius equal to 3.5 times internodal distance, a cubic spline
(43) 
for a weight function, Gaussian points for approximation of line integrals and points for approximating area integrals. Lagrange multipliers were used to impose essential boundary conditions.
The performance of EFG with respect to the number of nodes is much better than MLSM. However, a more fair comparison would take into account also a higher complexity of the EFG. This can be achieved by comparing error with respect to the number of MLS evaluations, which is the most time consuming part of the solution procedure. In Figure 5 it is demonstrated that although EFG provides much better results in comparison to MLSM at a given number of nodes, its accuracy becomes comparable to MLSM, when compared per number of MLS evaluations.
To asses the stability of the method regarding the nodal distribution, the following analysis was performed. A regular distribution of points as used in the solution in Figure 4 was distorted by adding a random perturbation to each internal node. Its position is altered by
(44) 
where is the distance to the closest node, and measuring the accuracy of the solution with respect to , representing magnitude of the perturbation. An example of original and perturbed node distributions are shown in Figure 6.
Accuracy of the solution with respect to the perturbation magnitude is presented in Figure 7. It is demonstrated that using monomials as a basis with 9 support nodes results in an unstable setup. On the other hand monomials with 13 support nodes are much more stable and equally accurate, while using Gaussian basis with high shape parameter is the most unstable setup. To mitigate the stability issue, a lower shape parameter can be chosen, however, at the cost of accuracy. Regardless of the setup one can expect the solution to be stable at least up to . Note that using more nodes in support domain can also increase stability. Refer to Kouszeg1 () for more details.
Time spent on each part of the solution procedure is shown in
Figure 8. All measurements were performed on a laptop
computer with an Intel(R) Core(TM) i74700MQ @2.40GHz CPU and with
16 of DDR3 RAM. MLSM is implemented in C++ utils_web ()
and compiled using g++ 7.1.7
for Linux with
std=c++14 O3 DNDEBUG flags.
It can be seen that solving the system (21) makes up for more than 50 of total time spent. Around 70 of that time is spent on computing the preconditioner. The only other significant factor is computing the shape functions taking approximately 40 of total time. Domain construction and matrix assembly take negligible amounts of time, matching the predictions made by complexity analysis in section 4.4.
To emphasize the generality of MLSM method, a “drilled” domain is considered in the next step. Arbitrarily positioned holes are added to the rectangular domain. The positioning algorithm described in section 4.2 and hrefinement algorithm described in section 4.3 are used to distribute the nodes inside the domain and refine the areas around the holes. The boundary conditions in this example are on the right, traction free on the inside of the holes and on top and bottom and uniform load of on the left. The computed solution is shown in Figure 9 along with the ordinary cantilever beam example.
Both solutions are coloured using von Mises stress , computed for the plane stress case as
(45) 
To further illustrate the generality of the method, an even more deformed domain is considered (Figure 10).
5.2 Hertzian contact
Another more interesting case arises from basic theory of contact mechanics, called Hertzian contact theory williams2001contact (). Consider two cylinders with radii and and parallel axes pressed together by a force per unit length of magnitude . The theory predicts they form a small contact surface of width , where
(46) 
and
(47)  
(48) 
Elastic modulus and Poisson’s ratio for the first material are denoted with and , and with and for the second material. The pressure distribution between the bodies along the contact surface is semielliptical, i.e. of the form
(49) 
A problem can be reduced to two dimensions using plane stress assumption. A special case of this problem is when , and , describing a contact of a cylinder and a half plane. This is the second numerical example tackled in this paper. The setup is ideal for testing the refinement, since a pronounced difference in behaviour of numerical solution near the contact in comparison to the rest of the domain is expected.
A displacement field satisfying (15) on with boundary conditions
(50)  
(51) 
is sought. Vector represents traction force on the surface and the upwards direction. Analytical solution for internal stresses in the plane in general point is calculated using the method of complex potentials mcewen1949stresses () and the stresses are given in terms of and , defined as
(52)  
(53) 
where in . The stresses are then expressed as
(54)  
(55)  
(56) 
Numerically the problem is solved by truncating the infinite domain to a rectangle for large enough and setting the essential boundary conditions everywhere but on the top boundary. The top boundary has a traction boundary condition with normal traction given by and no tangential traction. An illustration of the problem domain along with the boundary conditions is given in Figure 11.
The described contact problem is solved numerically and the error is measured between calculated and given stresses in relative norm as before, using
as an error indicator. Values , , , were chosen for the physical parameters of the problem. These values yield contact halfwidth and peak pressure . A value of for domain height is chosen, approximately 38 times greater than width of the contact surface. Convergence of the method is shown in Figure 12.
It is clear that the convergence of the method is very irregular and slow. This is to be expected, as means only approximately 30 nodes positioned within the contact surface, and that naturally leads to large changes as a change of a single node bears a relatively high influence. Another problem is that the boundary conditions are only continuous, and exhibit no higher regularity, not even Lipschitz continuity. The accuracy of the approximation may seem bad, but is in fact comparable to the cantilever beam case. Using the comparable value of nodes in the contact area it can be seen from Figure 5 that the approximation using this number of nodes in cantilever beam case achieved very similar results.
The total error of the approximation is composed of two main parts, the truncation error due to the nonexact boundary conditions and the discretization error, due to solving a discrete problem instead of the continuous one. First, we analyse the total error in terms of domain height . A graph showing the total error with respect to domain height is shown in Figure 13.
The total error decreases as domain height increases, regardless of the discretization density used. However, as soon as truncation error becomes lower than discretization error, increasing the height further yields little to no gain in total error. The higher the discretization density is, the later this happens. When convergence of a method stops or significantly decreases in order, an error limit imposed by the truncation error was reached.
It soon becomes impossible to uniformly increase discretization density due to limited resources, and the immediate solution is to refine the discretization in the contact area with the hrefinement algorithm introduced in section 4.3. A domain of height is chosen. Primary refinement is done in rectangle areas of the form
and secondary refinement around points on the surface is done in rectangle areas
The refined domain as described above is shown in Figure 14. This domain was used to solve the considered contact problem. Different levels of secondary refinement were tested to prove that refinement helps with accuracy. Convergence of the method on the refined domain is shown in Figure 15.
Comparing Figure 15 to Figure 12, it can be seen that refinement greatly improves the accuracy of the method. Using nodes without refinement yields worse results than nodes with only primary refinement. Each additional level of secondary refinement helps to decrease the error even further while keeping the same order of convergence. A solution of the problem on the final mesh is shown in Figure 16.
6 Conclusions
A MLSM solution of a linear elasticity problem on regular and irregular domains with a refined nodal distribution of two different numerical examples is presented in this paper. The method is analysed in terms of accuracy by comparison against available closed form solutions and by comparison against weak form EFG method. The convergence of the method is evaluated with respect to the number of computational nodes, selection of different basis functions, different refinement strategies and different boundary conditions. MLSM is also analysed from complexity point of view, first, theoretically, and then also experimentally by timing the computer execution time of all main blocks of the method. It is clearly demonstrated that the method is accurate and stable. Furthermore, it is demonstrated that nodal adaptivity is mandatory when solving contact problems in order to obtain accurate results and that the proposed MLSM method can handle extensive refinement with the smallest internodal distance being times smaller than the initial one. It is also demonstrated that proposed MLSM configuration can handle computations in complex domains.
In our opinion the presented meshless setup can be used, not only to solve academic cases with the sole goal to show excellent convergences, but also in more complex engineering problems. The C++ implementation of presented MLSM is freely available at utils_web ().
In future work we will continue to develop a meshless solution of a contact problem with a final goal to simulate a crack propagation due to the fretting fatigue pereira2016convergence () in a general 3D domain with added padaptivity to treat singularities near the crack tip.
Acknowledgement
The authors would like to acknowledge the financial support of the Research Foundation Flanders (FWO), The Luxembourg National Research Fund (FNR) and Slovenian Research Agency (ARRS) in the framework of the FWO Lead Agency project: G018916N Multianalysis of fretting fatigue using physical and virtual experiments and the ARRS research core funding No. P20095.
References
 [1] O. C. Zienkiewicz, R. L. Taylor, The Finite Element Method: Solid Mechanics, ButterworthHeinemann, 2000.
 [2] J. H. Hattel, P. N. Hansen, A control volumebased finite difference method for solving the equilibrium equations in terms of displacements, Applied Mathematical Modelling 19 (1994) 210–243.
 [3] Y. D. Freyer, C. Bailey, M. Cross, C.H. Lai, A control volume procedure for solving the elastic stressstrain equations on an unstructured mesh, Applied Mathematical Modelling 15 (1991) 639–645.
 [4] Y. Chen, J. D. Lee, A. Eskandarian, Meshless methods in solid mechanics, Springer, New York, NY, 2006.
 [5] B. Mavrič, B. Šarler, Local radial basis function collocation method for linear thermoelasticity in two dimensions, International Journal of Numerical Methods for Heat and Fluid Flow 25 (2015) 1488–1510.
 [6] R. Trobec, G. Kosec, Parallel scientific computing: theory, algorithms, and applications of mesh based and meshless methods, Springer, 2015.
 [7] G. Kosec, A local numerical solution of a fluidflow problem on an irregular domain, Advances in Engineering Software.
 [8] G. Kosec, R. Trobec, Simulation of semiconductor devices with a local numerical approach, Engineering Analysis with Boundary Elements 50 (2015) 69–75.
 [9] C. A. Wang, H. Sadat, C. Prax, A new meshless approach for three dimensional fluid flow and related heat transfer problems, Computers and Fluids 69 (2012) 136–146.
 [10] S. Li, K. W. Liu, Meshfree and particle methods and their applications, Applied Mechanics Reviews 55 (2002) 1–34.
 [11] S. N. Atluri, The Meshless Method (MLPG) for Domain and BIE Discretization, Tech Science Press, Forsyth, 2004.

[12]
J. Zhang, H. Chen, C. Cao,
A graphics processing
unitaccelerated meshless method for twodimensional compressible flows,
Engineering Applications of Computational Fluid Mechanics 11 (1) (2017)
526–543.
arXiv:http://dx.doi.org/10.1080/19942060.2017.1317027, doi:10.1080/19942060.2017.1317027.
URL http://dx.doi.org/10.1080/19942060.2017.1317027 
[13]
W. Li, G. Song, G. Yao,
Piecewise
moving least squares approximation, Applied Numerical Mathematics 115 (2017)
68–81.
doi:10.1016/j.apnum.2017.01.001.
URL http://linkinghub.elsevier.com/retrieve/pii/S0168927417300016 
[14]
K. Kovářík, S. Masarovičová, J. Mužík, D. Sitányiová,
A
meshless solution of two dimensional multiphase flow in porous media,
Engineering Analysis with Boundary Elements 70 (2016) 12–22.
doi:10.1016/j.enganabound.2016.05.008.
URL http://linkinghub.elsevier.com/retrieve/pii/S0955799716300911  [15] G. Kosec, B. Šarler, Hadaptive local radial basis function collocation meshless method, CMC: Computers, Materials, and Continua 26 (3) (2011) 227–253.
 [16] W. S. Slaughter, The linearized theory of elasticity, Springer Science & Business Media, 2012.
 [17] K. Pereira, S. Bordas, S. Tomar, R. Trobec, M. Depolli, G. Kosec, M. Abdel Wahab, On the convergence of stresses in fretting fatigue, Materials 9 (8) (2016) 639.

[18]
G. Kosec, M. Kolman, J. Slak,
Utilities for solving pdes
with meshless methods (2016).
URL https://gitlab.com/e62Lab/e62numcodes.git  [19] K. Reuther, B. Sarler, M. Rettenmayr, Solving diffusion problems on an unstructured, amorphous grid by a meshless method, International Journal of Thermal Sciences 51 (2012) 16–22.
 [20] J. Amani, M. Afshar, M. Naisipour, Mixed discrete least squares meshless method for planar elasticity problems using regular and irregular nodal distributions, Engineering Analysis with Boundary Elements 36 (2012) 894–902.
 [21] R. Löhner, E. Oñate, A general advancing front technique for filling space with arbitrary objects, International Journal for Numerical Methods in Engineering 61 (12) (2004) 1977–1991.
 [22] Y. Liu, Y. Nie, W. Zhang, L. Wang, Node placement method by bubble simulation and its application, CMES: Computer Modeling in Engineering & Sciences 55 (1) (2010) 89.
 [23] N. A. Libre, A. Emdadi, E. J. Kansa, M. Rahimian, M. Shekarchi, A stabilised rbf collocation scheme for neumann type boundary conditions., CMES: Computer Modeling in Engineering & Sciences 24 (2008) 61–80.
 [24] G. C. Bourantas, E. D. Skouras, G. C. Nikiforidis, Adaptive support domain implementation on the moving least squares approximation for mfree methods applied on elliptic and parabolic pde problems using strongform description, CMES: Computer Modeling in Engineering & Sciences 43 (2009) 1–25.
 [25] T. Belytschko, Y. Y. Lu, L. Gu, Elementfree galerkin methods, Int. J. Numer. Methods Eng. 37 (2) (1994) 229–256.
 [26] H. A. Van der Vorst, BiCGSTAB: A fast and smoothly converging variant of BiCG for the solution of nonsymmetric linear systems, SIAM Journal on scientific and Statistical Computing 13 (2) (1992) 631–644.
 [27] Y. Saad, ILUT: A dual threshold incomplete LU factorization, Numerical linear algebra with applications 1 (4) (1994) 387–402.
 [28] G. Kosec, P. Zinterhof, Local strong form meshless method on multiple graphics processing units, CMES: Computer Modeling in Engineering & Sciences 91 (5) (2013) 377–396.
 [29] V. P. Nguyen, T. Rabczuk, S. Bordas, M. Duflot, Meshless methods: a review and computer implementation aspects, Mathematics and computers in simulation 79 (3) (2008) 763–813.
 [30] J. A. Williams, R. S. DwyerJoyce, Contact between solid surfaces, Modern tribology handbook 1 (2001) 121–162.
 [31] E. M’Ewen, Stresses in elastic cylinders in contact along a generatrix (including the effect of tangential friction), The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science 40 (303) (1949) 454–459. doi:10.1080/14786444908521733.