Hybrid Fundamental Solution Based Finite Element Method: Theory and Applications

An overview on the development of hybrid fundamental solution based finite element method (HFS-FEM) and its application in engineering problems is presented in this paper. The framework and formulations of HFS-FEM for potential problem, plane elasticity, three-dimensional elasticity, thermoelasticity, anisotropic elasticity, and plane piezoelectricity are presented. In this method, two independent assumed fields (intraelement filed and auxiliary frame field) are employed.The formulations for all cases are derived from the modified variational functionals and the fundamental solutions to a given problem. Generation of elemental stiffness equations from the modified variational principle is also described. Typical numerical examples are given to demonstrate the validity and performance of the HFS-FEM. Finally, a brief summary of the approach is provided and future trends in this field are identified.


Introduction
A novel hybrid finite element formulation, called the hybrid fundamental solution based FEM (HFS-FEM), was recently developed based on the framework of hybrid Trefftz finite element method (HT-FEM) and the idea of the method of fundamental solution (MFS) [1][2][3][4][5].In this method, two independent assumed fields (intraelement filed and auxiliary frame field) are employed and the domain integrals in the variational functional can be directly converted to boundary integrals without any appreciable increase in computational effort as in HT-FEM [6][7][8].It should be mentioned that the intraelement field of HFS-FEM is approximated by the linear combination of fundamental solutions analytically satisfying the related governing equation, instead of -complete functions as in HT-FEM.The resulting system of equations from the modified variational functional is expressed in terms of symmetric stiffness matrix and nodal displacements only, which is easy to implement into the standard FEM.It is noted that no singular integrals are involved in the HFS-FEM by locating the source point outside the element of interest and do not overlap with field point during the computation [9].
The HFS-FEM mentioned above inherits all the advantages of HT-FEM over the traditional FEM and the boundary element method (BEM), namely, domain decomposition and boundary integral expressions, while avoiding the major weaknesses of BEM [10][11][12], that is, singular element boundary integral and loss of symmetry and sparsity [13].The employment of two independent fields also makes the HFS-FEM easier to generate arbitrary polygonal or even curvesided elements.It also obviates the difficulties that occur in HT-FEM [14,15] in deriving -complete functions for certain complex or new physical problems [16].The HFS-FEM has simpler expressions of interpolation functions for intraelement fields (fundamental solutions) and avoids the coordinate transformation procedure required in the HT-FEM to keep the matrix inversion stable.Moreover, this approach also has the potential to achieve high accuracy using coarse meshes of high-degree elements, to enhance insensitivity to mesh distortion, to give great liberty in element shape, and to accurately represent various local effects (such as hole, crack, and inclusions) without troublesome mesh adjustment [17][18][19][20].Additionally, HFS-FEM makes it possible for a more 2 Advances in Mathematical Physics flexible element material definition which is important in dealing with multimaterial problems, rather than the material definition being the same in the entire domain in BEM.However, we noticed that there are also some limitations of HFS-FEM, for example, determining the positions of source points used for approximation interpolations.It is also known that fundamental solution based approximations can perform remarkably well in smooth problems but tend to deteriorate when high-gradient stress fields are presented.
This paper is organized as follows: in Section 2, the basic idea and formulations of the HFS-FEM are presented through a simple potential problem.Then, plane elasticity problems are described in Section 3. Section 4 extends the 2D formulations of the HFS-FEM to general three-dimensional (3D) elasticity problems.The method of particular solution and radial basis function approximation are shown to deal with body force in this Section.In Section 5 we extend the HFS-FEM to thermoelastic problems with arbitrary body force and temperature change.In Section 6, the HFS-FEM for 2D anisotropic elastic materials is described based on the powerful Stroh formalism.Plane piezoelectric problem is discussed in Section 7. Finally, typical numerical examples are presented in Section 8 to illustrate applications and performance of the HFS-FEM.Concluding remarks and future development are discussed at the end of this paper.

Basic Equations of Potential
Problems.The Laplace equation of a well-posed potential problem (e.g., heat conduction) in a general plane domain Ω can be expressed as [21,22] ∇ 2  (x) = 0 ∀x ∈ Ω (1) with the boundary conditions where  is the unknown field variable and  represents the boundary flux,   is the th component of outward normal vector to the boundary Γ = Γ  ∪ Γ  , and  and  are specified functions on the related boundaries, respectively.The space derivatives are indicated by a subscript comma, that is,  , = /  , and the subscript index  takes values (1, 2) for twodimensional and (1, 2, 3) for three-dimensional problems.Additionally, the repeated subscript indices imply summation convention.
For convenience, (3) can be rewritten in matrix form as with 2.2.Assumed Independent Fields.In this section, the procedure for developing a hybrid finite element model with fundamental solution as interior trial function is described based on the boundary value problem defined by ( 1)- (3).Similar to the conventional FEM and HT-FEM, the domain under consideration is divided into a series of elements [15,16,21,[23][24][25][26][27][28][29][30].In each element, two independent fields are assumed in the way as described in [31] and are given in Section 2.2.
2.2.1.Intraelement Field.Similar to the method of fundamental solution (MFS) in removing singularities of fundamental solution, for a particular element  occupying subdomain Ω  , we assume that the field variable defined in the element domain is extracted from a linear combination of fundamental solutions centered at different source points located outside the element (see Figure 1): where   is undetermined coefficients,   is the number of virtual sources, and   (x, y  ) is the fundamental solution to the partial differential equation: as Evidently, (5) analytically satisfies (1) due to the solution property of   (x, y  ).
In implementation, the number of source points is taken to be the same as the number of element nodes, which is free of spurious energy modes and can keep the stiffness equations in full rank, as indicated in [21].The source point y  ( = 1, 2, . . .,   ) can be generated by means of the method employed in the MFS [32][33][34][35] where  is a dimensionless coefficient, x 0 is the point on the element boundary (the nodal point in this work), and x  is the geometrical centroid of the element (see Figure 1).Determination of  was discussed in [31,36] and  = 5-10 is usually used in practice.
The corresponding outward normal derivative of   on Γ  is where with

Auxiliary Frame Field.
In order to enforce the conformity on the field variable , for instance,   =   on Γ  ∩ Γ  of any two neighboring elements  and , an auxiliary interelement frame field ũ is used and expressed in terms of the same degrees of freedom (DOF) d as those used in the conventional finite elements.In this case, ũ is confined to the whole element boundary as ũ (x) = Ñ (x) d  (12) which is independently assumed along the element boundary in terms of nodal DOF d  , where Ñ (x) represents the conventional finite element interpolating functions.For example, a simple interpolation of the frame field on a side with three nodes of a particular element can be given in the form where Ñ ( = 1, 2, 3) stands for shape functions in terms of natural coordinate  defined in Figure 2.

Modified Variational Principle.
For the boundary value problem defined in (1)-( 3) and (5), since the stationary conditions of the traditional potential or complementary variational functional cannot guarantee the interelement continuity condition required in the proposed HFS FE model, as in the HT FEM [21,26], a variational functional corresponding to the new trial functions should be constructed to assure the additional continuity across the common boundaries Γ Ief between intraelement fields of element "" and element "" (see Figure 3) [36,37]: Advances in Mathematical Physics A modified variational functional is developed as follows: where in which the governing equation ( 1) is assumed to be satisfied, a priori, for deriving the HFS FE model.The boundary Γ  of a particular element consists of the following parts: where Γ Ie represents the interelement boundary of the element "" shown in Figure 3.
To show that the stationary condition of the functional (15) leads to the governing equation (Euler equation), boundary conditions, and continuity conditions, invoking (16) and (15) gives the following functional for the problem domain: (18) from which the first-order variational yields Using divergence theorem we can obtain For the displacement-based method, the potential conformity should be satisfied in advance then ( 21) can be rewritten as The Euler equation and boundary conditions can be obtained as using the stationary condition Π  = 0.
As for the continuous requirement between two adjacent elements "" and "" given in (14), we can obtain it in the following way.When assembling elements "" and ", " we have from which the vanishing variation of Π (+) leads to the reciprocity condition   +   = 0 on the interelement boundary Γ Ief .
If the following expression: is uniformly positive (or negative) in the neighborhood of {} 0 , where the displacement {} 0 has such a value that Π  ({} 0 ) = (Π  ) 0 and where (Π  ) 0 stands for the stationary value of Π  , we have in which the relation that {ũ}  = {ũ}  is identical on Γ  ∩Γ  has been used.This is due to the definition in ( 14) in Section 2.3.
Proof.For the proof of the theorem on the existence of extremum, we may complete it by the so-called "second variational approach" [38].In doing this, performing variation of Π  and using the constrained conditions (26), we find Therefore the theorem has been proved from the sufficient condition of the existence of a local extreme of a functional [38].

Element Stiffness Equation.
With the intraelement field and frame field independently defined in a particular element (see Figure 1), we can generate element stiffness equation by the variational functional derived in Section 2.3.Following the approach described in [21], the variational functional Π  corresponding to a particular element  of the present problem can be written as Appling the divergence theorem (20) to the functional (29), we have the final functional for the HFS-FE model Then, substituting (5), (9), and (12) into the functional (30) produces in which The symmetry of H  is obvious from the scalar definition (31) of variational functional Π  .
To enforce interelement continuity on the common element boundary, the unknown vector c  should be expressed in terms of nodal DOF d  .The minimization of the functional Π  with respect to c  and d  , respectively, yields from which the optional relationship between c  and d  and the stiffness equation can be produced where K  = G   H −1  G  stands for the element stiffness matrix.

Numerical
Integral for H and G Matrix.Generally, it is difficult to obtain the analytical expression of the integral in (32) and numerical integration along the element boundary is required.Herein the widely used Gaussian integration is employed [22].
For the H  matrix, one can express it as by introducing the matrix function Equation ( 36) can be further rewritten as where and  is the Jacobean expressed as where Thus, the Gaussian numerical integration for H matrix can be calculated by where   is the number of edges of the element and   is the Gaussian sampling points employed in the Gaussian numerical integration.Similarly, we can calculate the G  matrix using

Advances in Mathematical Physics
The calculation of vector g  in (32) is the same as that in the conventional FEM, so it is convenient to incorporate the proposed HFS-FEM into the standard FEM program.Besides, the flux is directly computed from (9).The boundary DOF can be directly computed from (12) while the unknown variable at interior points of the element can be determined from (5) plus the recovered rigid-body modes in each element, which is discussed in the following section.

Recovery of Rigid-Body Motion.
Considering the physical definition of the fundamental solution, it is necessary to recover the missing rigid-body motion modes from the above results.Following the method presented in [21], the missing rigid-body motion can be recovered by writing the internal potential field of a particular element  as where the undetermined rigid-body motion parameter  0 can be calculated using the least square matching of   and ũ at element nodes which finally gives in which Δ  = (ũ  − N  c  )| node and  is the number of element nodes.Once the nodal field is determined by solving the final stiffness equation, the coefficient vector c  can be evaluated from (34), and then  0 is evaluated from (45).Finally, the potential field  at any internal point in an element can be obtained by means of (43).

Plane Elasticity Problems
3.1.Linear Theory of Plane Elasticity.In linear elastic theory, the strain displacement relations can be used and equilibrium equations refer to the undeformed geometry [39].In the rectangular Cartesian coordinates ( 1 ,  2 ), the governing equations of a plane elastic body can be expressed as If written as matrix form, it can be presented as where is a body force vector, and the differential operator matrix L is given as where  = [ 11  22  12 ]  is a strain vector and u = [ 1 ,  2 ]  is a displacement vector.The constitutive equations for the linear elasticity are given in matrix form as where D is the material coefficient matrix with constant components for isotropic materials, which can be expressed as follows: where for plane stress. ( The two different kinds of boundary conditions can be expressed as where t = [ 1  2 ]  denotes the traction vector and A is a transformation matrix related to the direction cosine of the outward normal Substituting ( 49) and ( 50) into (47) yields the well-known Navier partial differential equations in terms of displacements: 3.2.Assumed Independent Fields.For elasticity problem, two different assumed fields are employed as in potential problems: intraelement and frame field [1,31,36].The intraelement continuity is enforced on nonconforming internal displacement field chosen as the fundamental solution of the problem [36].The intraelement displacement fields are approximated in terms of a linear combination of fundamental solutions of the problem of interest: where   is again the number of source points outside the element domain, which is equal to the number of nodes of an element in the present research based on the generation approach of the source points [31]., in which x and y  are, respectively, the field point and source point in the local coordinate system ( 1 ,  2 ).The components  *  (x, y  ) are the fundamental solution, that is, induced displacement component in -direction at the field point x due to a unit point load applied in -direction at the source point y  , which are given by [40,41] where The virtual source points for elasticity problems are generated in the same manner as that in potential problems described in Section 2.
With the assumption of intraelement field in (56), the corresponding stress fields can be obtained by the constitutive equation (50): where As a consequence, the traction is written as in which The components  *  (x, y) for plane strain problems are given as The unknown c e in ( 56) is calculated by a hybrid technique [31], in which the elements are linked through an auxiliary conforming displacement frame which has the same form as that in conventional FEM (see Figure 1).This means that, in the HFS-FEM, a conforming displacement field should be independently defined on the element boundary to enforce the field continuity between elements and also to link the unknown c with nodal displacement d  .Thus, the frame is defined as where the symbol "∼" is used to specify that the field is defined on the element boundary only, Ñ is the matrix of shape functions, and d  is the nodal displacements of elements.
Taking the side 3-4-5 of a particular 8-node quadrilateral element (see Figure 1) as an example, Ñ and d  can be expressed as , where Ñ1 , Ñ2 , and Ñ3 can be expressed by natural coordinate as 3.3.Modified Functional for the Hybrid FEM.As in Section 2, HFS-FE formulation for a plane elastic problem can also be established by the variational approach [36].In the absence of body forces, the hybrid functional Π  used for deriving the present HFS-FEM can be constructed as [22] where ũ and   are the intraelement displacement field defined within the element and the frame displacement field defined on the element boundary, respectively.Ω  and Γ  are the element domain and element boundary, respectively.Γ  , Γ  , and Γ  stand, respectively, for the specified traction boundary, specified displacement boundary, and interelement boundary (Γ  = Γ  + Γ  + Γ  ).Compared to the functional employed in the conventional FEM, the present variational functional is constructed by adding a hybrid integral term related to the intraelement and element frame displacement fields to guarantee the satisfaction of displacement and traction continuity conditions on the common boundary of two adjacent elements.

Advances in Mathematical Physics
By applying the Gaussian theorem, (67) can be simplified to Due to the satisfaction of the equilibrium equation with the constructed intraelement field, we have the following expression for HFS finite element model: The variational functional in (69) contains boundary integrals only and will be used to derive HFS-FEM formulation for the plane isotropic elastic problem.

Element Stiffness Matrix.
As in Section 2, the element stiffness equation can be generated by setting Π  = 0. Substituting ( 56), (64), and (61) into the functional of (69), we have where To enforce interelement continuity on the common element boundary, the unknown vector c  should be expressed in terms of nodal DOF d  .The stationary condition of the functional Π  with respect to c  and d  yields, respectively, (33) and (34).

Recovery of Rigid-Body Motion.
For the same reason stated in Section 2.6, it is necessary to reintroduce the discarded rigid-body motion terms after we have obtained the internal field of an element.The least squares method can be employed for this purpose and the missing terms can be recovered easily by setting for the augmented internal field [22] where the undetermined rigid-body motion parameter  0 can be calculated using the least square matching of   and ũ at element nodes which finally gives where ] , in which Δ  = (ũ  − û ) ( = 1, 2) and  is the number of element nodes.Once the nodal field is determined by solving the final stiffness equation, the coefficient vector c  can be evaluated from ( 56), and then  0 is evaluated from (74).
Finally, the displacement field u  at any internal point in an element can be obtained by (72).

Three-Dimensional Elastic Problems
In this section, the HFS-FEM approach is extended to threedimensional (3D) elastic problem with/without body force.
The detailed 3D formulations of HFS-FEM are firstly derived for elastic problems by ignoring body forces, and then a procedure based on the method of particular solution and radial basis function approximation are introduced to deal with the body force [42].As a consequence, the homogeneous solution is obtained by using the HFS-FEM and the particular solution associated with body force is approximated by using the strong form of basis function interpolation.

Governing Equations and Boundary
Conditions.Let ( 1 ,  2 ,  3 ) denote the coordinates in a Cartesian coordinate system and consider a finite isotropic body occupying the domain Ω, as shown in Figure 4.The equilibrium equation for this finite body with body force can be expressed as The constitutive equations for linear elasticity and the kinematical relation are given as where   is the stress tensor,   is the strain tensor, and   is the Kronecker delta.Substituting ( 77) into ( 76), the equilibrium equation is rewritten as For a well-posed boundary value problem, the boundary conditions are prescribed as follows: where Γ  ∪Γ  = Γ is the boundary of the solution domain Ω,   and   are the prescribed boundary values.In the following parts, we will present the procedure for handling the body force appearing in (78).

The Method of Particular Solution.
The inhomogeneous term   associated with the body force in (78) can be effectively handled by means of the method of particular solution presented in [22].In this approach, the displacement   is decomposed into two parts, a homogeneous solution  ℎ  and a particular solution where the particular solution    should satisfy the governing equation: without any restriction of boundary condition.However, the homogeneous solution should satisfy with the modified boundary conditions From the above equations it can be seen that once the particular solution    is known, the homogeneous solution  ℎ  in ( 82) and ( 83) can be obtained using HFS-FEM.The final solution can then be given by ( 80).In the next section, radial basis function approximation is introduced to obtain the particular solution, and the HFS-FEM is presented for solving (82) and (83).

Radial Basis Function Approximation.
For body force   , it is generally impossible to find an analytical solution which enables us to convert the domain integral into a boundary one.So we must approximate it by a combination of basis (trial) functions or other methods with the HFS-FEM.Here, we use radial basis function (RBF) [33,43] to interpolate the body force.We assume where  is the number of interpolation points,   are the RBFs, and    are the coefficients to be determined.Subsequently, the particular solution can be approximated by where Φ   is the approximated particular solution kernel of displacement satisfying (86) below.Once the RBFs are selected, the problem of finding a particular solution is reduced to solve the following equation: To solve this equation, the displacement is expressed in terms of the Galerkin-Papkovich vectors Substituting ( 87) into (86), we can obtain the following biharmonic equation: Taking the Spline Type RBF  =  2−1 as an example, we get the following solutions: where , Advances in Mathematical Physics and   represents the Euclidean distance between a field point ( 1 ,  2 ,  3 ) and a given point ( 1 ,  2 ,  3 ) in the domain of interest.The corresponding particular solution of stresses can be obtained as where  = (2V/(1−2V)).Substituting ( 90) into (92), we have where , (94)

HFS-FEM for Homogeneous Solution.
After obtaining the particular solution in Sections 4.2 and 4.3, we can determine the modified boundary conditions (83).Finally, we can treat the 3D problem as a homogeneous problem governed by ( 82) and ( 83) by using the HFS-FEM presented below.It is clear that once the particular and homogeneous solutions for displacement and stress components at nodal points are determined, the distribution of displacement and stress fields at any point in the domain can be calculated using the element interpolation function.However, for 3D elasticity problems in the absent of body force, that is,   = 0, the procedures in Sections 4.2 and 4.3 will become unnecessary and we can employ the following procedures to find the solution directly.

Assumed Intraelement and Auxiliary Frame Fields.
The intraelement displacement fields are approximated in terms of a linear combination of fundamental solutions of the problem as where the matrix N e and unknown vector c e can be further written as ] in which x and y  are, respectively, the field point and source point in the local coordinate system ( 1 ,  2 ).The fundamental solution  *  (x, y  ) is given by [40] where ,   is the number of source points.The source point y  ( = 1, 2, . . .,   ) can also be generated by means of the following method [36] as in two-dimensional cases where  is a dimensionless coefficient, x 0 is the point on the element boundary (the nodal point in this work), and x  is the geometrical centroid of the element (see Figure 5).According to ( 50) and ( 49), the corresponding stress fields can be expressed as where The components  *  (x, y) are given by As a consequence, the traction can be written in the form in which To link the unknown c  and the nodal displacement d  , the frame is defined as where the symbol "∼" is used to specify that the field is defined on the element boundary only, Ñ is the matrix of shape functions, and d  is the nodal displacements of elements.Taking the surface 2-3-7-6 of a particular 8-node brick element (see Figure 5) as an example, matrix Ñ and vector d  can be expressed as where the shape functions are expressed as where and (  ,   ) is the natural coordinate of the -node of the element (Figure 6).

Modified Functional for Hybrid Finite Element Method.
In the absence of body forces, the hybrid functional Π  used for deriving the present HFS-FEM is written as [22] By applying the Gaussian theorem, (108) can be simplified as Due to the satisfaction of the equilibrium equation with the constructed intraelement fields, we have the following expression for HFS finite element model: The functional (110) contains only boundary integrals of the element and will be used to derive HFS-FEM formulation for the three-dimensional elastic problem in the following section.

Element Stiffness
Matrix.Substituting (95), (102), and (104) into the functional (110), we have where To enforce interelement continuity on the common element boundary, the unknown vector c  should be expressed in terms of nodal DOF d  .The stationary condition of the functional Π  with respect to c  and d  yields again, respectively, ( 33) and (34).

Numerical Integral over Element.
Considering a surface of the 3D hexahedron element, as shown in Figure 6, the vector normal to the surface can be obtained by where V  and V  are the tangential vectors in the -direction and -direction, respectively, calculated by where   is the number of nodes of the surface and (  ,   ,   ) are the nodal coordinates.Thus the unit normal vector is given by where is the Jacobian of the transformation from Cartesian coordinates (, ) to natural coordinates (, ).

Advances in Mathematical Physics
For the  matrix, we introduce the matrix function Then we can get and we rewrite it to the component form as Using the relationship d =  (, ) d d (120) and the Gaussian numerical integration, we can obtain where   and   are, respectively, the number of surface of the 3D element and the number of Gaussian integral points in each direction of the element surface.Similarly, we can calculate the G  matrix by F [ (, ) ,  (, )]  (, ) d d It should be mentioned that the calculation of vector g  in equation is the same as that in the conventional FEM, so it is convenient to incorporate the proposed HFS-FEM into the standard FEM program.Besides, the stress and traction estimations are directly computed from (99) and (100), respectively.The boundary displacements can be directly computed from (104) while the displacements at interior points of element can be determined from (95) plus the recovered rigid-body modes in each element, which is introduced in the following section.

Recovery of Rigid-Body Motion Terms.
As in Section 2, the least square method is employed to recover the missing terms of rigid-body motions.The missing terms can be recovered by setting for the augmented internal field and using a least-square procedure to match  ℎ and ũℎ at the nodes of the element boundary The above equation finally yields where ] ,

Thermoelasticity Problems
Thermoelasticity problems arise in many practical designs such as steam and gas turbines, jet engines, rocket motors, and nuclear reactors.Thermal stress induced in these structures is one of the major concerns in product design and analysis.The general thermoelasticity is governed by two time-dependent coupled differential equations: the heat conduction equation and the Navier equation with thermal body force [44].In many engineering applications, the coupling term of the heat equation and the inertia term in Navier equation are generally negligible [44].As a consequence, most of the analyses are employing the uncoupled thermoelasticity theory which is adopted in this topic [45][46][47][48][49][50][51][52].In this section, the HFS-FEM is presented to solve 2D and 3D thermoelastic problems with considering arbitrary body forces and temperature changes [53].The method used herein is similar to that in Section 4.

Basic Equations for Thermoelasticity.
Consider an isotropic material in a finite domain Ω and let ( 1 ,  2 ,  3 ) denote the coordinates in Cartesian coordinate system.The equilibrium governing equations of the thermoelasticity with the body force are expressed as where   is the stress tensor,   is the body force vector, and ,  = 1, 2, 3.The generalized thermoelastic stress-strain relations and kinematical relation are given as where   is the strain tensor,   is the displacement vector,  is the temperature change,  is the shear modulus, ] is the Poisson's ratio,   is the Kronecker delta, and is the thermal constant with  being the coefficient of linear thermal expansion.Substituting (128) into (127), the equilibrium equations may be rewritten as For a well-posed boundary value problem, the following boundary conditions, either displacement or traction boundary condition, should be prescribed as where Γ  ∪ Γ  = Γ is the boundary of the solution domain Ω,   and   are the prescribed boundary values, and is the boundary traction, in which   denotes the boundary outward normal.

The Method of Particular Solution.
For the governing equation (130) given in the previous section, the inhomogeneous term  , −   can be eliminated by employing the method of particular solution [16,36,44].Using superposition principle, we decompose the displacement   into two parts, the homogeneous solution  ℎ  and the particular solution    as follows: but does not necessarily satisfy any boundary condition.It should be pointed out that its solution is not unique and can be obtained by various numerical techniques.However, the homogeneous solution should satisfy with modified boundary conditions From above equations, it can be seen that once the particular solution is known, the homogeneous solution  ℎ  in ( 135)-(137) can be solved by (135).In the following section, RBF approximation is described to illustrate the particular solution procedure, and the HFS-FEM is presented which can be used for solving (135)-(137).

Radial Basis Function Approximation.
RBF is to be used to approximate the body force   and the temperature field  in order to obtain the particular solution.To implement this approximation, we may consider two different ways: one is to treat body force   and the temperature field  separately as in Tsai [54].The other is to treat  , −   as a whole [53].
Here we demonstrated that the performance of the latter one is usually better than the former one.

Interpolating Temperature and Body Force Separately.
The body force   and temperature  are assumed to be by the following two equations: where  is the number of interpolation points,   are the basis functions, and    and   are the coefficients to be determined by collocation.Subsequently, the approximate particular solution can be written as follows: where Φ   and Ψ   are the approximated particular solution kernels.Once the RBF is selected, the problem of finding a particular solution will be reduced to solve the following equations: To solve (140), the displacement is expressed in terms of the Galerkin-Papkovich vectors [43,[55][56][57]] Substituting ( 142) into (140), one can obtain the following biharmonic equation: If taking the Spline Type RBF  =  2−1 , one can get the following solutions: where , for two-dimensional problem and where , for three-dimensional problem, where   represents the Euclidean distance of the given point ( 1 ,  2 ) from a fixed point ( 1 ,  2 ) in the domain of interest.The corresponding stress particular solution can be obtained by where  = (2V/(1 − 2V)).Substituting (147) into (149), one can obtain where for two-dimensional problem, and substituting (147) into (149), one can obtain where , for three-dimensional problem.
To solve (141), one can treat Ψ  as the gradient of a scalar function Substituting ( 154) into (141) obtains the Poisson's equation Thus, taking  =  2−1 , its particular solution can be obtained [55] as follows: Then from (154) we can get Ψ  as follows: The corresponding stress particular solution can be obtained by substituting (147) into Then we have Advances in Mathematical Physics

Interpolating Temperature and Body Force Together.
Considering the temperature gradient plays the role of body force, we can approximate  , −   together by the following equation: Thus the approximate particular solution can be written as Consequently, we can follow the same way as that for body force in Method 1 and employing ( 140) and ( 142)-( 150) to obtain the desired Φ   and   , which are the same as those for body force case only.It is noted that Method 2 has a relatively smaller number of equations to solve for the coefficients and the condition number of the coefficient matrix is smaller as well, which will be beneficial to the solution.
Once we have obtained the particular solutions of (133), we can use them to get the modified boundary conditions in (136) to obtain the homogeneous solution by considering (135).Then, we can employ the HFS-FEM described in Section 3 for 2D problem and Section 4.4 for 3D problem to obtain the homogeneous solutions.

Anisotropic Composite Materials
In materials science, composite laminates are usually assemblies of layers of fibrous composite materials which can be joined together to provide required engineering properties, such as specified in-plane stiffness, bending stiffness, strength, and coefficient of thermal expansion [58].Individual layers (or laminas) of the laminates consist of highmodulus, high-strength fibers in a polymeric, metallic, or ceramic matrix material.From the viewpoint of micromechanics, the fiber and matrix in each lamina can be treated as the inclusion and matrix, respectively.On the other hand, from the viewpoint of macromechanics, both a lamina and the whole laminate can be viewed as a general anisotropic body by classical lamination theory.Hence, the analysis of anisotropic bodies is important for understanding of the micro-or macromechanical behavior of composites [58].
In the literature, there are two main approaches dealing with generalized two-dimensional anisotropic elastic problems.One is Lekhnitskii formalism [59,60] which begins with stresses as basic variables and the other is Stroh formalism [61,62] which starts with displacements as basic variables.Both of them are formulated in terms of complex variable functions.The Stroh formalism, which has been shown to be elegant and powerful, is used to find the analytical solutions for the corresponding infinite bodies [61,63].The formalism is also widely employed in the derivation of the inclusion or crack problems of anisotropic materials [64,65].Because of the limitations of the analytical solutions which are only available for some problems with simple geometry and boundary conditions [66,67], numerical methods such as finite element method (FEM), boundary element method (BEM), mesh free method (MFM), and Hybrid-Trefftz (HT) FEM are usually resorted to solve more complex problems with complicated boundary constraints and loading conditions [68][69][70].In this section, we presented the HFS-FEM for analyzing anisotropic composite materials based on the associated fundamental solutions in terms of Stroh formalism [71].

Basic Equations and Stroh
Formalism.In the Cartesian coordinate system ( 1 ,  2 ,  3 ), if we neglect the body force   , the equilibrium equations, stress-strain laws, and straindisplacement equations for anisotropic elasticity are [61]  , = 0, (162) where ,  = 1, 2, 3,   is the fourth-rank anisotropic elasticity tensor.The equilibrium equations can be rewritten in terms of displacements by substituting ( 163) and ( 164) into (162) as The boundary conditions of the boundary value problem (163)-(165) are where   and   are the prescribed boundary displacement and traction vector, respectively.In addition,   is the unit outward normal to the boundary and Γ = Γ  + Γ  is the boundary of the solution domain Ω.
For the generalized two-dimensional deformation of anisotropic elasticity   is assumed to depend on  1 and  2 only.Based on this assumption, the general solution to (165) can be written as [61,62] where u = ( 1 ,  2 ,  3 )  is the displacement vector,  = ( 1 ,  2 ,  3 )  is the stress function vector, and f() = [ 1 ( 1 ),  2 ( 2 ),  3 ( 3 )]  is a function vector composed of three holomorphic complex functions   (  ),  = 1, 2, 3, which is an arbitrary function with argument   =  1 +    2 and will be determined by satisfying the boundary and loading conditions of a given problem.In (167), Re stands for the real part of a complex number,   are the material eigenvalues with positive imaginary part, and A = [a 1 , a 2 , a 3 ] and B = [b 1 , b 2 , b 3 ] are 3 × 3 complex matrices formed by the material eigenvector associated with   , which can be obtained by the following Eigen relations [61]: where N is a 6 × 6 foundational elasticity matrix and  is a 6 × 1 column vector defined by where and the matrices Q, R and T are 3 × 3 matrices extracted from   as follows: The stresses can be obtained from the derivation of stress functions  as follows: where Thus, the Green's function satisfying the above boundary conditions is found to be [72] Therefore, fundamental solutions of the problem can be expressed as The corresponding stress components can be obtained from stress function  as where p are chosen to be (1, 0, 0)  , (0, 1, 0)  , and (0, 0, 1)  , respectively, ⟨⋅⟩ stands for the diagonal matrix corresponding to subscript , Im denotes the imagery part of a complex number, and superscript  denotes the matrix transpose.For the two coordinate systems mentioned in Figure 7, the angle between the axis- 1 and axis- is denoted by , which is positive along the anticlockwise direction, and the relationship for transformation of stress and strain components between the local material coordinates and global coordinates is given by where the transformation matrix T and its inverse matrix are defined as where the matrix N e and unknown vector c e can be further written as ] , in which   is the number of source points, x and y  are, respectively, the field point and source point in the coordinate system ( 1 ,  2 ) local to the element under consideration.The fundamental solution  *  (x, y  ) is given by (175) for general elements. is a dimensionless coefficient for determining source points as described in previous sections.
The corresponding stress fields can be expressed as where The components  *  (x, y) are given by (176) when p is selected to be (1, 0, 0)  , (0, 1, 0)  , and (0, 0, 1)  , respectively.As a consequence, the traction can be written as in which The unknown c  in ( 180) and (182) may be calculated using a hybrid technique [31], in which the elements are linked through an auxiliary conforming displacement frame which has the same form as in conventional FEM (see Figure 1).Thus, the frame field is defined as where the symbol "∼" is used to specify that the field is defined on the element boundary only, Ñ is the matrix of shape functions, and d  is the nodal displacements of elements.
Taking the side 3-4-5 of a particular 8-node quadrilateral element (see Figure 1) as an example, Ñ and d  can be expressed as where the shape functions are expressed as and Ñ1 , Ñ2 , and Ñ3 are expressed by natural coordinate  ∈ (189)

Modified Functional for HFS-FEM.
With the assumption of two distinct intraelement field and frame field for elements, we can establish the modified variational principle based on ( 165) and (166) for the hybrid finite element method of anisotropic materials [21,73].In the absence of the body forces, the hybrid variational functional Π  for a particular element  is constructed as where the boundary Γ  of the element  is and Γ  is the interelement boundary of element .Performing a variation of Π  , one obtains Applying Gaussian theorem and the definitions of traction force we obtain Then, substituting (195) into (192) gives Considering the fact that we can finally obtain the following form: Therefore, the Euler equations for (190) result in ( 165) and (166) because the quantities   ,   , and ũ  may be arbitrary.
As for the continuity condition between elements, it can be easily seen from the following variational of two adjacent elements such as  and This indicates that the stationary condition of the functional satisfies both the required boundary and interelement continuity equations.In addition, the existence of extremum of functional (190) can be easily proved by the so-called "second variational approach" as well, which indicates functional (190) has a local extreme.Therefore, we can conclude that the variational functional (190) can be used for deriving hybrid finite element formulations.

Element Stiffness Equation.
Using Gaussian theorem and equilibrium equations, all domain integrals in (190) can be converted into boundary integrals as follows: Substituting ( 180), (184), and (186) into the functional (200) yields the formulation as where The stationary condition of the functional Π  with respect to c  and d  , respectively, yields ( 33) and (34).

Piezoelectric Materials
Piezoelectric materials have the property of converting electrical energy into mechanical energy and vice versa.This reciprocity in the energy conversion makes them very attractive for using in electromechanical devices, such as sensors, actuators, transducers, and frequency generators.To enhance understanding of the electromechanical coupling mechanism in piezoelectric materials and to explore their potential applications in practical engineering, numerous investigations, either analytically or numerically, have been reported over the past decades [73][74][75][76][77][78][79][80][81].In this section, the HFS-FEM is developed for modeling two-dimensional piezoelectric material [82,83].The detailed formulations and procedures are given in the following sections.

Basic Equations for Piezoelectric Materials.
For a linear piezoelectric material in absence of body forces and electric charge density, the differential governing equations in the Cartesian coordinate system   ( = 1, 2, 3) are given by where   is the stress tensor,   is the electric displacement vector, and Ω is the solution domain.With strain and electric field as the independent variables, the constitutive equations are written as where   is the elasticity tensor measured under zero electric field,   and   are, respectively, the piezoelectric tensor and dielectric tensor measured under zero strain.  and   are the elastic strain tensor and the electric field vector, respectively.The relation between the strain tensor   and the displacement   is given by and the electric field component   is related to the electric potential  by The boundary conditions of the boundary value problem ( 203)-( 206) can be defined by where   ,   ,   , and  are, respectively, the prescribed boundary displacement, the prescribed traction vector, the prescribed surface charge, and the prescribed electric potential.In addition,   is the unit outward normal to the boundary and is the boundary of the solution domain Ω.
For the transversely isotropic material, if  1 - 2 is taken as the isotropic plane, one can employ either  1 - 3 or  2 - 3 plane to study the plane electromechanical phenomenon.Thus, choosing the former and considering the plane strain conditions ( 22 =  32 =  12 = 0 and  2 = 0) (204) can be reduced to For the plane stress piezoelectric problem ( (213) 7.2.Assumed Independent Fields.For the piezoelectric problems, HFS-FEM is based on assuming two distinct displacement and electric potential (DEP) fields: intraelement DEP field u and an independent DEP frame field ũ along element boundaries [21,29].

Intraelement Field.
The intraelement DEP field u identically fulfills the governing differential equations ( 203) and is approximated by a linear combination of foundational solutions at different source points located outside of the element domain where the fundamental solution matrix N e is now given by , in which x and y  are, respectively, the field point (i.e., the nodal points of the element in this work) and source point in the local coordinate system ( 1 ,  2 ).The component  *  (x, y  ) is the induced displacement component ( = 1, 2) or electric potential ( = 3) in -direction at the field point  due to a unit point load ( = 1,2) or point charge ( = 3) applied in -direction at the source point y  .The fundamental solution  *  (x, y  ) is given as [76,82,84] where 2 and   is the three different roots of the characteristic equation  6   −  4  +  3  −  = 0.The source point y  ( = 1, 2, . . .,   ) can be generated by the following method [36]: Making use of (205) and the expression of intraelement DEP field u in (214), the corresponding stress and electric displacement in (212) can be written as where  = [ 11  22  12  1  2 ]  and in which  *  (x, y  ) denotes the corresponding stress components ( = 1, 2, 3) or electric displacement ( = 4, 5) along -direction at the field point  due to a unit point load ( = 1, 2) or a unit point charge ( = 3) applied in -direction at the source point  and can be derived from (216) and are listed below, which are derived by substituting Advances in Mathematical Physics the fundamental solutions into constitutive equations [85]: in which the coefficients   ,   ,   ,   ,  11 ,  12 , and  13 are defined as in literature [84].
From (204), (208), and (209), the generalized traction forces and electric displacement are given as Advances in Mathematical Physics

23
where 7.2.2.Auxiliary Frame Field.For the two-dimensional piezoelectric problem under consideration, the frame field is assumed as where Ñ is a matrix of the corresponding shape functions.
For the side 3-4-5 of a particular quadratic element as shown in Figure 1, the shape function matrix Ñ and nodal vector d e can be given in the form , where the shape functions Ñ are expressed by natural coordinate (225)

Variational Principles.
Based on the assumption of two distinct DEP fields, the Euler equations of the proposed variational functional should also satisfy the following interelement continuity requirements in addition to (203)-(210) [82,83]: where "" and "" stand for any two neighboring elements.Equations ( 203)-(210) together with ( 226) and ( 227) can now be taken as the basis to establish the modified variational principle for the hybrid finite element method of piezoelectric materials [21,29].
Since the stationary conditions of the traditional potential or complementary variational functional cannot satisfy the interelement continuity condition required in the proposed HFS-FEM, new modified variational functional should be developed.In the absence of the body forces and electric charge density, the hybrid functional Π  for a particular element  is constructed as where and the boundary Γ  of the element  is and Γ  is the interelement boundary of element .Compared to the functional employed in the conventional FEM, the present hybrid functional is constructed by adding two integral terms related to the intraelement and element frame DEP fields to guarantee the satisfaction of displacement and electrical potential continuity condition on the common boundary of two adjacent elements.
It can be proved that the stationary conditions of the above functional (228) lead to (203)-(210).To this end, performing a variation of Π  , one obtains [86] in which the first term is Advances in Mathematical Physics and the definitions of traction force and electrical displacement we obtain Then, substituting (235) into (231) gives Considering the fact that we finally get the following form: Therefore, the Euler equations for (238) result in (203)-( 210) and ( 226) because the quantities   ,   , ,   , ũ  , and  φ may be arbitrary.As for the continuity condition of (227), it can easily be seen from the following variation of two adjacent elements such as  and : which indicates that the stationary condition of the functional satisfies both the required boundary and interelement continuity equations.In addition, the existence of extremum of functional (228) can be easily proved by the "second variational approach" as well, which indicates functional (228) has a local extreme..The element stiffness equation can be generated by setting Π  = 0. To simplify the derivation, we first transform all domain integrals in (228) into boundary ones.With Gaussian theorem, the functional in (228) may be simplified as

Element Stiffness Equation
Due to the satisfaction of the equilibrium equation with the constructed intraelement fields, we have the following expression for the HT finite element model: Substituting ( 214), (223), and (221) into the above functional (241) yields the formulation as where  The stationary condition of the functional Π  with respect to c  and d  , respectively, yields (33) and (34).Following the procedure described in [21,22], the missing rigid-body motion can be recovered by setting the augmented internal field of a particular element  as where the undetermined rigid-body motion parameter c 0 can be calculated using the least square matching of u e and ũe at element nodes min which finally gives in which Δu  = (ũ  − û )| nodei and  is the number of element nodes.As a consequence, c 0 can be calculated by (246) once  the nodal DEP fields d  and the interpolation coefficients c  are, respectively, determined by ( 214) and (223).Then the complete DEP fields u  can be obtained from (244).

Normalization.
The order of magnitudes of the material constants and the corresponding field variables in piezoelectricity have a wide spectrum as large as 10 19 in SI unit.This will lead to ill-conditioned matrix of the system [72,87].To resolve this problem, normalization of each quantity by its reference value is employed.The reference values for the stiffness, piezoelectric stress constant, dielectric constants, and strain are selected to be  0 = 10 11 (N/m 2 ),  0 = 10 1 (N/mV),  0 = 10 −9 (C/mV), and  0 = 10 −3 (V/m), respectively.The reference values of other quantities, as Advances in Mathematical Physics shown in Table 1, are determined in terms of these four fundamental reference variables and the characteristic length  0 = 10 0 (m) of the problem so that the normalized governing equations remain in exactly the same form as the original equations.

Numerical Examples
Several numerical examples are presented in this section to illustrate the application of the HFS-FEM and to demonstrate its effectiveness and accuracy.Unless otherwise indicated, mesh convergence tests were conducted for the reference solutions obtained from ABAQUS in the following examples.

Cubic Block under Uniform Tension and Body Force.
An isotropic cubic block, with dimension 10 × 10 × 10 and subject to a uniform tension as shown in Figure 8, is considered in this example.A constant body force of 10 Mpa and uniform distributed tension of 100 MPa are applied to the cube.Figure 9 presents the displacement component  1 and the stress component  11 at point  of the block, which are calculated by the HFS-FEM on the three meshes with distorted 8-node brick elements (Figure 10).The results calculated by ABAQUS with a very fine mesh (with 40 × 40 × 40 C3D8 and EAS element with 68921 nodes) are taken as a reference value for comparison.The results from C3D8 and EAS elements in ABAQUS are also presented in Figure 9 for Table 1: Reference values for material constants and field variables in piezoelectricity derived from basic reference variables:  0 ,  0 ,  0 ,  0 , and  0 .
comparison.It can be seen from these figures that the results obtained from both the HFS-FEM and ABAQUS converge to the benchmark value with the number of degree of freedom (DOF) increasing.For Mesh 1, the hybrid EAS element has the best performance while for Mesh 2 and Mesh 3 it can be seen that HFS-FEM with 8-node brick elements exhibits better accuracy in both displacement and stress compared with EAS element in traditional FEM.Contour plots of  1 and  11 obtained by HFS-FEM on Mesh 3 are also presented in Figure 11.It should be noted that for problems involving body forces the accuracy of the RBF interpolation has to be considered for a satisfactory solution.The details on the RBF interpolation can be found in previous literatures [43,[88][89][90].Figure 13 presents the variation of the radial and the circumferential thermal stresses with the cylinder radius, respectively, in which the theoretical values are given for comparison [39].It is seen from Figure 13 that the results from Method 2 are much better than those obtained from Method 1 for both radial and circumferential stress.It can be inferred that the error may be to a large extent due to the RBF interpolation, for which the number of interpolation points has a significant influence on its accuracy.
Figure 14 displays the contour plots of (a) radial and (b) circumferential thermal stresses (the meshes used for contour plot are different from that for calculation due to using quadratic elements).It demonstrates that treating temperature gradient and body force together is more superior to dealing with them separately.However, we have to rely on Method 1 when the temperature change is discrete distribution or the gradient of the temperature field is not available.

3D Cube under Arbitrary Temperature and Body Force.
As shown in Figure 15, a 3D cube of 1 × 1 × 1 with center located at (0.5, 0.5, 0.5) is considered in this example.The material properties of the cube are Young's modulus  = 5000, Poisson's ratio ] = 0.3, and linear thermal expansion Because there is no analytical solution available, the results from ABAQUS herein are employed for comparison.The meshes used by HFS-FEM and ABAQUS are given in Figure 15.
Figure 16 presents the displacement  3 and stress  33 along one edge of the cube which is coinciding with  3 axis.It can be seen that the results from HFS-FEM again agree very well with those by ABAQUS.It is demonstrated that the proposed HFS-FEM is able to predict the response of 3D thermoelastic problems under arbitrary temperature and body force.It is also shown that HFS-FEM with RBF interpolation can give satisfactory results using coarse meshes.

Orthotropic Composite Plate with an Elliptic Hole under
Tension.A finite composite plate containing an elliptical hole (Figure 17) is investigated in this example.A uniform The length and width of the plate are  = 20 mm and  = 20 mm; the major and minor lengths of the ellipses  and  are, respectively, 2 mm and 1 mm.In the computation, plane stress condition is used.The mesh configurations used by HFS-FEM and ABAQUS are given in Figure 17.
Figure 18 shows the corresponding variation of the hoop stress along the rim of the elliptical hole when the orientation angle  of reinforced fibers is equal to 0 ∘ , 45 ∘ , and 90 ∘ , respectively.It is found from Figure 18 that the results from HFS-FEM have a good agreement with the reference solutions from ABAQUS.This indicates that the proposed method is able to capture the variations of the hoop stress induced by the elliptical hole in the plate.The contour plots of stress components  11 ,  22 , and  12 around the elliptic hole in the composite plate for several different fiber angles are shown in Figure 19.
Figure 20 shows the stress concentration factor (SCF) along with the inclined angle  of the reinforced fibers, which exhibits a good agreement with the solutions from ABAQUS.It is obvious that the SCF of the punched plate rises with the increasing fiber angle .It is found from Figure 20 that the largest SCF occurs at  = 90 ∘ , whereas the smallest appears at  = 0 ∘ .It indicates the effectiveness of the proposed method in predicting the SCF for anisotropic composites as well.

Isotropic Plate with Multianisotropic Inclusions.
In this example, a multi-inclusion problem is investigated to show the capability of the HFS-FEM to deal with both isotropic and anisotropic materials in a unified way.As shown in Figure 21, an isotropic plate containing multianisotropic inclusions of square geometry (edge length  = 2) is considered.The distance between any two inclusions is assumed to be  = 3.The material parameters for the inclusions are chosen as  1 = 134.45GPa,  2 =  3 = 11.03GPa,  23 = 2.98 GPa,  31 =  12 = 28.5 GPa, and V 23 = 0.49, V 31 = V 12 = 0.301.The material parameters for isotropic matrix are elastic modules  = 2.8 Gpa and poison's ratio ] = 0.3.The mesh configuration of the plate for HFS-FEM is shown in Figure 21, which uses 272 quadratic general elements.
In general, the Stroh formalism is suitable for the anisotropic material with distinct material eigenvalues, and it fails for the degenerated materials like isotropic material with repeated eigenvalues   = , ( = 1, 2, 3) [59].However, a small perturbation of the material constants, such as  1 = (1 − 0.004),  2 = , and  3 = (1 + 0.004), can be applied to make the eigenvalues distinct and the results can be applied conveniently.
Table 2 shows the displacement and stresses at points A, B, and C as indicated in Figure 21.It is observed that there is a good agreement between the results by the HFS-FEM and those from ABAQUS using very fine mesh, in which the maximum relative error for displacement and stress by HFS-FEM occur at Point B (i.e.,  2 = 0) and are 0.7% and 1.3%, respectively.Additionally, it can be found that the results from HFS-FEM are better than those from ABAQUS using the same mesh.Although the stress at the vicinity of the inclusions (Point C) has a little degradation due to the high stress concentration and stress contrast at the adjacent elements, the displacement still agrees well with the reference solution.The variations of displacement components  1 and  2 along the right edge ( = 8) by HFS-FEM are shown in Figure 22.
8.6.Infinite Piezoelectric Medium with Hole.Consider an infinite piezoelectric plane with a circular hole as shown in Figure 23.The material parameters are given in Table 3. Suppose that mechanical load  ∞  =  0 = 10 parallel to the  axis is imposed at infinity with traction and electric chargefree at the boundary of the hole.In our calculation, the radius of the hole is set to be  = 1 and / = 20 is employed in the analysis.
Figure 24 shows the distribution of hoop stress   and radial stress   along the line  = 0 for remote loading  ∞  and along the line  = 0 for remote loading  ∞  , respectively.Figure 25 presents the variations of the normalized stress   / 0 and the normalized electric displacement   / 0 × 10 10 along the hole edge under remote mechanical loading.It is obvious that the results obtained from HFS-FEM agree well with the results from ABAQUS and Sosa [74].
It can be seen from Figure 24 that hoop stress   has maximum value on the rim of the hole and it decreases dramatically with the increase of the distance from the hole edge.It is also shown that   tends to be equal to the remote applied load  0 when  increases toward infinity.Compared with the hoop stress   , it is obvious that the radial stress   is much smaller and usually does not need to be considered.It is obvious that loading along the poling direction will produce smaller stress concentration due to coupling effect.The maximum values of   appear at  = 90 ∘ for case of  ∞  and at  = 0 ∘ and  = 180 ∘ for case of loading  ∞  , both of which agree well with the analytical solution from Sosa [74].The electric displacement   / 0 × 10 10 produced by  ∞  and  ∞  is nearly the same and is symmetrical with respect to the  axis.It is found that the maximum values of   appear at  = 65 ∘ and  = 115 ∘ , which also agrees well with the analytical solution.

Conclusions
In this paper, we have reviewed the HFS-FEM and its application in engineering applications.The HFS-FEM is a promising numerical method for solving complex engineering problems.The main advantages of this method include integration along the element boundaries only, easily adopting arbitrary polygonal or even curve-sided elements and symmetric and sparse stiffness matrix, and avoiding the singularity integral problem as encountered in BEM.Moreover, as in HT-FEM, this method offers the attractive possibility to develop accurate crack singular, corner, or perforated elements, simply by using appropriate special fundamental solutions as the trial functions of the intraelement displacements.It is noted that the HFS-FEM has attracted more attention of researchers in computational mechanics in the past few years, and good progress has been made in the field of potential problems, plane elasticity, piezoelectric problems, and so on.However, there are still many possible extensions and areas in need of further development in the future: (1) to develop various special-purpose elements to effectively handle singularities attributable to local geometrical or load effects (holes, cracks, inclusions, interface, corner, and load singularities), with the special-purpose functions warranting that excellent results are obtained at minimal computational cost and without local mesh refinement, (2) to extend the HFS-FEM to elastodynamics, fluid flow, thin and thick plate bending, and fracture mechanics, (3) to develop efficient schemes for complex engineering structures and improve the related general purpose computer codes with good preprocessing and postprocessing capabilities, (4) to extend this method to the case of multifield problems such as thermoelastic-piezoelectric materials and thermomagnetic-electric-mechanical materials, and to develop multiscale framework across from continuum to micro-and nanoscales for modeling heterogeneous materials.
Frame field ũ(x) = e d e u = N e c e

Figure 1 :
Figure 1: Intraelement field and frame field of a HFS-FEM element for 2D potential problems.

Figure 2 :
Figure 2: Typical quadratic interpolation for frame field.

Figure 4 :
Figure 4: Geometrical definitions and boundary conditions for a general 3D solid.

Figure 5 :
Figure 5: Intraelement field and frame field of a hexahedron HFS-FEM element for 3D elastic problem (the source points and centroid of 20-node element are omitted in the figure for clarity and clear view, which is similar to that of the 8-node element).

Figure 6 :
Figure 6: Typical linear interpolation for the frame fields of 3D brick elements.

) 6 . 1 . 2 .
Foundational Solutions.To find the fundamental solution needed in our analysis, we have to first derive the Green's function of the problem: an infinite homogeneous anisotropic elastic medium loaded by a concentrated point force (or line force for two-dimensional problems) p = ( p1 , p2 , p3 ) applied at an internal point x = (x 1 , x2 ) far from the boundary.The boundary conditions of this problem can be written as ∫  d = p for any closed curve  enclosing x, ∫  du = p for any closed curve , lim x → ∞   = 0.

Figure 8 :Figure 9 :Figure 10 :
Figure 8: Cubic block under uniform tension and body force: geometry, boundary condition, and loading.

Figure 13 :
Figure 13: (a) Radial thermal stresses and (b) circumferential thermal stresses with the cylinder radius.

Figure 14 :
Figure 14: Contour plots of (a) radial and (b) circumferential thermal stresses (the background mesh in (a) and (b) is used for plots of the calculated results in postprocessing).

Figure 16 :
Figure 16: (a) Displacement  3 and (b) stress  33 along one cube edge when subjected to arbitrary temperature and body force.

8. 2 .
Circular Cylinder with Axisymmetric Temperature Change.In this example, a long circular cylinder with axisymmetric temperature change in domain is considered.Both inside and outside surfaces of the cylinder are assumed to be free from traction.The temperature  changes logarithmically along the radial direction.With the symmetry condition of the problem, only one quarter of the cylinder is modeled.The configurations of geometry and the boundary conditions are shown in Figure 12.In our computation, the parameters  = 5,  = 20,  = 1000, ] = 0.3,  = 0.001, and  0 = 10.The two approaches listed in Section 4 to approximate the body force and temperature are discussed and analyzed in this example.

Figure 17 :
Figure 17: (a) Schematic of an orthotropic composite plate with an elliptic hole under uniform tension, and its mesh configurations for (b) HFS-FEM, 1515 quadratic elements, and (c) ABAQUS, 9498 quadratic elements.

Figure 18 :
Figure 18: Variation of hoop stresses along the rim of the elliptical hole for different fiber orientation .

Figure 22 :
Figure 22: The variation of displacement component (a)  1 and (b)  2 along the right edge of the plate ( = 8) by HFS-FEM and ABAQUS.

Figure 23 :Figure 24 :Figure 25 :
Figure 23: (a) An infinite piezoelectric plate with a circular hole subjected to remote stress.(b) Mesh used by HFS-FEM.(c) Mesh used by ABAQUS.

Table 2 :
Comparison of displacement and stress at points A and B.

Table 3 :
Properties of the material PZT-4 used for the model.