A Parametric Study of Mesh Free Interpolation Based Recovery Techniques in Finite Element Elastic Analysis

The paper presents a parametric study on interpolation techniques based postprocessed error estimation in finite element elastic analysis by varying important parameters of recovery, interpolation scheme and type of patch construction. The quality of error estimation with recovery parameters is compared in terms of local and global effectivity of error estimation, rate of error convergence, and adaptively refined meshes. A mesh free moving least square interpolation technique with proven reliability and effectivity is introduced for improving the recovery of finite element solution errors. The post-processed finite element solutions of elastic problems are presented for performance study under different parameters of recovery technique. The study concludes that recovery parameters of interpolation method have pronounced effect on the recovery of finite element solution error and analysis in adaptive environment.


Introduction
The finite element method yields approximate solutions for the problems and the errors are inherent in the method. The finite element error is of two types namely, discretization and elemental error [Dow (1999)]. Elemental errors arise from modeling deficiencies in individual elements. They depend upon the algorithm employed for solving system equations and the method employed for computing derivatives of state variables. Elemental error can produce serious error in the computed solution unless special care is taken. The discretization errors are produced because the sub regions are not capable of representing the full range of behavior of the continuum. The discretization errors manifest themselves as discontinuities or jumps in components of stress or strain between elements and as misrepresentation of the actual boundary stresses. The magnitude of the jumps is indicative of the size of the local discretization error. The latter can be identified by whether the finite element solution satisfies equilibrium conditions at every point of the problem domain. The discretization error can be reduced by increasing the number of elements. The discretization error can be assessed by various developed estimation techniques. However, interpolation techniques based error estimation is most popular that rely upon the so-called "recovery" of a higher order approximation of the computed variable of interest. The development of recovery of solution errors, a-priori and aposteriori-error estimation techniques are underway for the last several years. An extensive survey of a-posteriori error estimation is due to Ainsworth et al. [Ainsworth and Oden (1997); Gratsch and Bathe (2005); Ainsworth and Oden (1997); Gratsch and Bathe (2005)]. Mirzaei [Mirzaei (2015)] provides a comprehensive analysis of error estimation considering moving least square approximation. An error estimation method based on side curvature has been presented by Xing et al. [Xing, Wang and Makinouchi (1999)]. It was used to estimate the accuracy of finite element discretization of a deformed sheet, especially one with sharp or quasi-sharp corners. Liu et al. [Liu and Elmaraghy (1992)] discussed discretized errors, and their estimation in conjunction with quadrilateral finite element meshes which are generated by an intelligent mesh generator. Zienkiewicz et al. [Zienkiewicz and Zhu (1987)] proposed an error estimator which is not only accurate but whose evaluation is computationally so simple that it can be readily implemented in existing finite element codes. The estimator computes energy norm of error on global as well as local levels. Liu et al. [Liu, Nguyen and Lam (2008)] have proposed technique to obtain the exact finite solution to mechanics problem using the alpha finite element method (αFEM). Kim et al. [Kim and Lee (2010)] applied the equilibrated residual method with set of Neumann data for nonconforming elements to improve the finite element linear elasticity solution. Ullah et al. [Ullah, Coombs and Augarde (2013)] proposed a coupled finite element method (FEM) and the element-free Galerkin method (EFGM) error estimation based on local maximum entropy shape functions, for linear and nonlinear problems. Many different post-processing techniques can be employed to improve the quality of the derivatives of the finite element approximations such as averaging [Bramble and Schatz (1977)], local or global projections [Hinton and Campbell (1974)], and those exploring the super-convergence phenomenon. The nodal average method is a simple technique that leads to optimal convergence rates when only linear elements are considered. A general recovery technique has been developed by Zienkiewicz et al. [Zienkiewicz and Zhu (1992)] for determining the derivatives. The implementation of the recovery technique is simple and cost effective. The technique has been tested for a group of linear, quadratic and cubic elements for both one and two dimensional problems. Numerical experiments demonstrate that the recovered nodal values of the derivatives with linear and cubic elements are super-convergent. One order higher accuracy is achieved by the procedure with linear and cubic elements but two order higher accuracy is achieved for derivatives with quadratic elements. In particular, an O(h4) convergence of the nodal values of the derivatives for a quadratic triangular element is reported. Zienkiewicz et al. [Zienkiewicz and Zhu (1992)] derived a theorem on dependence of the effectivity index for the Zienkiewicz-Zhu error estimator and the convergence rate of a recovered solution. They found that with super-convergent recovery, the effectivity index tends asymptotically to unity. Numerical tests on various element types illustrate the effectivity of their error estimator in the energy norm and pointwise gradient (stress) error estimation. Li et al. [Li and Wiberg (1994)] presented a post-processing technique for obtaining a-posteriori error estimators both in the energy and L2 norms, an element patch that represents the union of the considered element and its neighbours is introduced. The post-processing for determining more accurate solutions is made by fitting a higher order polynomial expansion to the computed solution at super-convergent points in the patch. The element error estimate norms are calculated directly from the improved solutions. The recovery by equilibrium in patches (REP) and the recovery by compatibility in patches (RCP) to obtain the improved stresses to the computed stress are due to Ubertini [Ubertini (2004)]. Parret et al. [Parret-Fré aud, Rey, Gosselet et al. (2016)] propose a moving least squares (MLS) recovery-based procedure to obtain post-processed smoothed stresses field in which the continuity of the recovered field is provided by the shape functions of the underlying mesh. Sharma et al. [Sharma, Zhang, Langelaar et al. (2018)] propose a stress recovery procedure for low-order finite elements in 3D in which the recovered stress field is obtained by satisfying equilibrium in an average sense and by projecting the directly calculated stress field onto a conveniently chosen space. More recently, Ahmed et al. [Ahmed, Singh and Deshmukh (2018)] has effectively implemented the element free Galerkin based recovery technique in finite element analysis of elastic problems with super-convergent properties. The researchers also present the improvement in the finite element method. Bathe et al. [Bathe and Zhang (2017)] has proposed a new finite element solution scheme with overlapping elements and present its complete formulation for twodimensional problems. Aragon et al. [Aragon and Simone (2017)] proposes a Discontinuity-Enriched Finite Element Method (DEFEM) for modeling problems with both weak and strong discontinuities independently of the finite element discretization and which adds enriched degrees of freedom only to nodes created at the intersection between a discontinuity and edges of elements in the mesh. The finite element analysis coupled with post-processing technique is influenced by several recovery parameters (i.e., parameters affecting the post-processed results) including meshing scheme, post-processed field variable, error norm, interpolation type for post processing, shape and size of influence domain, dilation parameter, order of polynomial expansion in basis function, weighting function, etc. There are only a few published works on the role of recovery parameters in interpolation based post-processed finite element procedures. Rajendran et al. [Rajendran and Liew (2003)] have proposed best-fit approach for predicting optimal stress sampling points for the patch recovery of nodal stresses. Wang et al. [Wang and Liu (2002); Kanber, Bozkurt and Erklig (2013); Wang and Liu (2002); Kanber, Bozkurt and Erklig (2013)] have propose optimal shape parameters in uniform grids while Perko et al. [Perko and Sarler (2007)] proposes shape parameters optimization in non-uniform grids. Neil et al. [Neil, Atluri and Zuo (2006)] propose a practical mathematics model for obtaining the optimal radius of support of radial weights used in MLS methods. Onate et al. [Onate, Perazzo and Miquel (2001)] have discussed the completeness, robustness, continuity and accuracy aspect of moving least squares (MLS) interpolation procedure. The present work deals with the study of influence of parameters on interpolation recovery based finite element elastic analysis. The moving least squares (MLS) interpolation technique is introduced very recently for post-processing of finite element solutions, it will be interesting to quantify the effect of various factors of MLS interpolation on the error recovery process, Moreover, presently post-processed techniques use mesh reliant patches for recovery of discretization error in the finite element solution. The meshes less patches for interpolation are investigated in the present study. The discretization error at element and global levels in the finite element solution are presented in energy norm and elemental error in the solution is assumed constant with insignificant effect. The linear and quadratic triangular elements schemes are used for problem domain discretization. The four recoveres parameter namely, interpolation type for post-processing, shape of support domain, dilation parameter and order of polynomial expansion in basis function are selected for the parametric investigation. The effect of application of error estimation techniques under different recovery parameters are also explored for finite element analysis in adaptive environment.

Least square interpolation technique (Zienkiewicz-Zhu (ZZ) Recovery)
In ZZ super-convergent patch recovery technique [Zienkiewicz and Zhu (1992)], the nodal values of stress belong to a polynomial expansion of the same complete order as that present in the basis function and is valid over a patch of nodes surrounding the particular given node. Following polynomial expansiosn may be used for each component of stress.  * (X) = P(X). a in which P(X) is the basis of the assumed polynomial, X = (x1, x2), the spatial variable and a the vector of the unknown parameters. np is the total number of sampling points.

Moving Least Squares (MLS) interpolation recovery technique
The interpolation recovery technique derived the approximate functions using a set of nodes distributed over a domain, without depending on the meshing scheme, in least square sense. The Moving Least Square (MLS) approximation, u h (x), can be defined as,u h (x) = ∑ ( ) =  ( ) =0 where  (x) is shape function and m is the total number of terms in the basis function P(x). In the present study, m is taken as six for linear and nine for quadratic elements. P T (x) = [1, x, y, x 2 , xy, y 2 ] is taken for linear element and P T (x) = [1, x, y, x 2 , xy, y 2 , x 2 y, xy 2 , x 2 y 2 ] is taken for for quadratic element.
The 1 (x) is given as,  ( ) = ( ) −1 ( ) ( ) ( − ). The vector of coefficients a(x) can be obtained by minimizing a weighted residual as, = ∑ ( − =0 ) [ ( ) ( ) − ℎ ] 2 where n is the number of nodes i and w(x−xi) is a weight function in 2-D associated to each node (domain of influence of that node) which is usually built in such a way that it takes a unit value in the vicinity of the point where the function and its derivatives are to be computed and vanishes outside a region Ωi surrounding the point xi. The support of the shape function (x) is equivalent to the support of the weight function. The following cubic spline weight function with circular or rectangular domain of influence is considered in the present study. . The solution is acceptable if   allow where allow is the allowable accuracy. If not, h-refinement i.e., reduction of element size is needed.

Parametric study of interpolation recovery technique
Numerical experiments on plate problems using finite element method were carried out to study the effect of the recovery parameters (i.e., parameters affecting the post-processed results) in terms of effectivity, rate of convergence and adaptively refined mesh. The plate domain is discretized using linear and quadratic triangular elements. The error at element and global levels in the finite element solution are computed in energy norm. The effect on the performance of interpolation recovery based post processed technique is studied by varying four parameters namely, interpolation type for post-processing, support domain shape, dilation parameter and number of terms in basis function. The mesh free circular and rectangular support domain shape is used in MLS interpolation recovery technique. The mesh dependent node patch is also used for least square and MLS interpolation. Six different dilation parameters namely, 2.4, 3.0, 4.0, 4.5, 5.5 and 6.5, are selected for the moving least square interpolation formulation with consideration to order of element. Two set of number of basis function terms namely, 3, 6, 9 and 6,9 10 are selected for fitting the derivatives of field variable over the nodes in support domain to recover the solution error in linear and quadratic element respectively.

Elastic analysis of plate problem
The square elastic plate problem subjected to body loads (bx, by) is tested to study the effect of recovery parameters on post-processed finite element analysis of the problem. The body force of the plate is in the form of polynomials and close form displacements solution (u, v) due to body load are given in Eqs.

Recovery parameter: interpolation technique type and patch construction
Two types of interpolations formulation namely, least square (ZZ recovery) interpolation and moving least square (MLS) interpolation, for fitting the field variable or its derivatives over the mesh reliant and mesh free patch of nodes, are used in recovery based finite element parametric study. The mesh dependent node patch for interpolation is the patch of nodes surrounding the particular given node of the mesh as shown in Fig. 2. The mesh dependent node patch is considered for MLS and least square interpolation. The MLS interpolation technique uses mesh free patch, i.e., support domain in circular and rectangular shape. The interpolation formulation effect on post processed results is investigated keeping dilation parameter constant with a value of 3.0 and 4.5 for linear elements and quadratic elements respectively. The number of terms in polynomial basis is taken as 6 for linear elements while the number of terms is 9 for quadratic elements in MLS interpolation technique. Tabs. 1-4 give the error convergence and effectivity obtained in post processed finite element elastic analysis using least square (ZZ) interpolation and MLS interpolation based recovery for linear elements. The finite element elastic analysis results with node patch and support domain constructed from the meshes of quadratic elements are given in Tabs. 5 to 8. From the tables, Tabs. 1-8, it is clear that recovery parameter, interpolation techniques and, types of element and patch construction significantly affect the performance of postprocessing based finite element analysis. The tables depict that order of error obtained employing recovery procedures is much lower than FEM solution error and error convergence is also manifold as compared to the FEM solution error with reduction of mesh size. The effect on error order and convergence is further enhanced for higher order elements. The MLS interpolation recovery performs better over least square interpolation recovery with different elements and patch types giving higher error convergence and effectivity converging nearly to one with increase in mesh fineness. As far as the performance of element and patch type is concerned, the overall performance of lower order element and mesh free node patch (support domain) is superior. When comparing the support domain shape effect on the post-processed analysis, the rectangular support domain construction has better performance as compared to circular shape of support domain for both lower and higher order element. The convergence rate is higher and the effectivity is closer to one in case of rectangular support domain. Though, the order of error in energy norm obtained is higher in rectangular support domain.         The element effectivity frequency and stress error distribution in plate problem domain obtained with different recovery parameters is also plotted. The element effectivity frequency for linear and quadratic elements for support domain shape used in MLS interpolation based recovery technique is plotted in Fig. 3 at last refinement level considered. The effectivity frequency plots using mesh reliant patch ZZ recovery, mesh reliant patch MLS recovery and mesh free patch (support domain) based MLS recovery technique in recovery based finite element analysis are shown in Fig. 4. The Fig. 3 and Fig. 4 depicts that unstructured discretization with higher order elements is more sensitive to patch construction type in comparison to low order elements as the effectivity results for the analysis is similar for lower order elements. The least square interpolation recovery using mesh dependent patch using higher order element perform better than MLS interpolation recovery. The stress error distribution obtained with interpolation based recovery in plate domain using different patch construction type and structured/unstructured mesh are plotted in Fig. 5 and Fig. 6. It is evident that the error pattern over problem domain is clearly distinct for different patch construction type and structured/unstructured mesh. The error distribution using recovery considering unstructured mesh seems to be more evenly distributed than error distribution using structured mesh. Also, error recovery using higher order element become further smooth and problem domain is divided two separate zone of lower and higher order error. The error distributions plots for different recovery parameters show that the error distribution patterns for exact and recovered error distributions in plate domain are in good agreement. However, the order of errors is different with different recovery parameters.
(a) Linear Triangular Element (b) Quadratic Triangular Element  Adaptive analysis of the plate problem guided through the error estimation using interpolation based post processing technique for specified accuracy limit is also carried out. Adaptive results are obtained for mesh based patch ZZ recovery, mesh based patch MLS recovery and support domain based MLS recovery technique. Tab. 9 and Tab. 10 show global errors, the number of element and DOF in refined meshes, obtained from the linear and quadratic elements initial mesh, to bring the accuracy level at target accuracy. The global error in energy norm of higher order element is much lower than global error of lower order elements. The target error of 10% is selected for adaptive finite element analyses employing linear elements discretization and a target error of 1% is fixed for adaptive analyses employing quadratic elements. The computed global error for linear elements is highest for ZZ patch recovery and minimum is for MLS patch recovery technique. For quadratic elements, global error is highest for recovery technique using MLS rectangular support domain. The adaptively refined meshes with target error using least square and MLS interpolation based error estimation are given in Figs. 7-10. The initial meshes are adaptively refined to bring the solution error within the prescribed limit.
The number of elements required to achieve target accuracy is lesser in MLS interpolation based error estimation as compared to least square based error estimation for linear elements discretization. However, for quadratic element discretization, number of elements required to achieve target accuracy are more in MLS interpolation. It infers that MLS based recovery technique is more efficient than the ZZ recovery technique for linear elements discretization. The refined mesh plots show that the meshes become finer in areas of high local error and become coarser in area of lesser error to get a uniform accuracy throughout the domain. The error distribution patterns and zone of fine mesh in plate domains show good agreement. It concludes that the interpolation techniques based adaptive analysis may predict the high errors zones and control the solution error effectively.

Recovery parameter: dilation parameter of MLS interpolation technique
The dilation parameter regulates the density of the particle distribution in support domain. The suitable dilation parameter supports the MLS interpolation function to automatically satisfy its completeness condition. The dilation parameters values with linear element meshing considered in the study are 2.5, 3.0 and 4.0. For quadratic element meshing, the dilation parameters values are 4.5, 5.5 and 6.5. The convergence rate of error and effectivity of error estimation with different dilation parameter values with increasing fineness order of mesh in MLS interpolation based post-processed finite element analysis are given in Tab. 5 and Tab. 6. The tables show that the order of error in energy norm increases with the increase of dilation parameter and effectivity of error estimation decreases with increase of the dilation parameter and, the rate of convergence is not significantly affected by the higher value of dilation parameter in case linear element meshing. For higher order element, the order of error decreases and convergence rates increase with increase of dilation parameter and the finite element analysis with quadratic element shows improved performance under different dilation parameters. However, it is observed that there will be a particular value of dilation parameter to provide an optimal performance. The local (element) effectivity frequency using linear and quadratic elements for dilation parameter is plotted in Fig. 5 at last refinement level considered, to study the effect on recovery of finite element solution. The figure depicts that the errors using linear elements is not much affected by increasing the dilation parameter value as the number of element effectivity with value about one, is more or less same with different dilation parameter. However, number of elements of effectivity with values nearly to one are shown to increase in analysis using quadratic elements mesh.

Recovery parameter: number of basis function terms
In the post-processing technique, field variable or its derivative is assumed to belong to polynomial expansion of the same complete order or higher order polynomial as that present in the basis function. According to so-called ZZ super-convergent patch recovery technique [Zienkiewicz and Zhu (1992)], the nodal values of stress belongs to a polynomial expansion of the same complete order as that present in the basis function. According to Ahmed and Singh [Ahmed and Singh (2008)] for element patch, the nodal values of field variables belong to a higher order polynomial expansion. For order of polynomial expansion in the basis function effect study, two set of number of basis function terms namely, 3, 6, 9 and 6, 9, 10 are used, for fitting the derivatives of field variable over the mesh independent node patches (support domain), for linear and quadratic elements respectively. The convergence rate of error and effectivity of error estimation with set of number of basis function terms at different element fineness in MLS interpolation based post-processed finite element analysis are given in Tab. 7 and Tab. 8. From the tables, it can be concluded that moving least square based postprocessed analysis, for the field variable derivative interpolation requires polynomial expansion of the one order higher as that present in the basis function. The local (element) effectivity frequency using linear and quadratic elements for various number of basis function terms is shown in Fig. 6 at last refinement level considered. The figure depicts that performance of recovery based finite element analysis will not improve by taking higher number of basis function terms for lower order elements as the number of element of effectivity with value nearly one are not significantly increased by using more number of basis function terms.

Elastic plate problem with central hole
The elastic square plate with central circular opening, a typical stress concentration problem, is also selected for error recovery quality study of interpolation based recovery method under different interpolations and patches construction type. The side of the square plate is 5 unit with hole radius of one unit. The normal and vertical displacement components are zero along the circular arc boundary. The normal displacement component and shear stress are zero over the line of symmetry. Static boundary conditions are imposed from the given traction equations. The plate stress field are given in Eqs. (5) to (7) [Zienkiewicz and Zhu (1992)]. The discretized domain for plate problem with circular opening is shown in Fig. 13  (7) where r = y 2 + x 2 and σ∞ is the uniaxial traction applied at infinity.

Recovery parameter: interpolation technique and patch construction
The recovery parametric study of interpolation based recovery method for the interpolation technique and patch construction considers dilation parameter value as 3.5 for linear elements and 4.0 for quadratic elements. The number of terms in polynomial basis is taken as 6 and 9 respectively for linear and quadratic elements. The quality of error estimation of the recovery analysis using least square and moving least square techniques, linear and quadratic element meshing to form circular and rectangular support domain error are presented in terms of convergence and effectivity in Tab. 15-Tab. 18. The tables show that least square interpolation method give better quality of error estimation as compared to MLS interpolation method for plate problem with central opening. The lower order elements perform better as compared to higher order element for analysed stress concentration problem with post-processed based finite element analysis using interpolation based recovery techniques. The rectangular support domain has better performance as compared to circular shape of support domain in MLS interpolation based recovery analysis. The order of error in energy norm obtained is higher in rectangular support domain as compared to the other patch construction type used in interpolation based recovery analysis.     Tab. 19 shows global errors obtained with least square interpolation recovery (ZZ) and moving least square interpolation recovery considering node patches and support domains for central hole plate problem finite element analysis using the linear and quadratic elements. The computed global error for linear elements is highest for ZZ patch recovery and minimum is for MLS patch recovery technique. For quadratic elements, global error is highest for recovery technique using MLS rectangular support domain similar to the global error obtained for plate problem without hole.

Conclusions
The study presents discretization errors estimation under varying parameters of finite element error recovery techniques. The quality of error estimation of finite element numerical experiments, using the interpolation based recovery technique parameters, on elastic plate problems is investigated in terms of global and local effectivity, rate of error convergence and adaptively refined mesh at specified error. The linear and quadratic triangular elements have been used for the discretization of the problem domain. The effect on the performance of interpolation based recovery technique is studied by considering four recovery parameters namely, type of interpolation technique and patch construction for interpolation, dilation parameter and number of terms in basis function.
The study considers two interpolation technique namely, least square (ZZ recovery) and moving least square interpolation techniques employing mesh reliant as well as mesh free patches (support domains). The main conclusions of the study are summarized below. 1. The basic recovery parameters of interpolation based recovery technique have pronounced effect on the finite element analysis post processed results. 2. The finite element analysis (FEA) employing the error recovery procedures provides much lower order of errors as compared to FEA solution errors and error convergence is also much higher than the convergence of FEA solution errors with increasing order of meshing. 3. The moving least square interpolation based recovery technique is suitable to enhance the quality of error estimation for linear element discretization while for the higher order elements, least square (ZZ recovery) technique may be used for enhanced performance. 4. The quality of error estimation of mesh free patches (support domains) used interpolation based recovery technique is superior to mesh dependent patches based error estimation. 5. The rectangular support domain construction has better convergence performance as compared to circular shape of support domain for both lower and higher element mesh in mesh free recovery based finite element analysis. The convergence rate achieved is about two orders higher in case of rectangular support domain. 6. The computational analysis results obtained for error estimation quality vary with the variation of dilation parameters. A suitable value of the dilation parameter will provide optimal performance for the problem. 7. The improved quality of error estimation in moving least square based post-processed analysis is obtained when the field variable derivative is considered belong to a higher order polynomial expansion as that present in the basis function.
8 The error distribution patterns and zone of fine mesh in problem domain show good agreement and the interpolation techniques based adaptive analysis may be utilize to predict the high errors zones and control the solution error effectively.

Conflict of interest:
The author(s) declare(s) that there is no conflict of interest regarding the publication of this paper.