# A Brief Review on Polygonal/Polyhedral Finite Element Methods

Authors: Logah Perumal
Journal: Mathematical Problems in Engineering
 (2018)

---

## Abstract

This paper provides brief review on polygonal/polyhedral finite elements. Various techniques, together with their advantages and disadvantages, are listed. A comparison of various techniques with the recently proposed Virtual Node Polyhedral Element (VPHE) is also provided. This review would help the readers to understand the various techniques used in formation of polygonal/polyhedral finite elements.


---

## Body

## 1. Introduction

Element equations are obtained by incorporating nodal conditions of the element geometry into the shape functions. One of the requirements is that the field variable obtained from the element equations should be linear on the element boundaries. This requirement is met for triangular and quadrilateral elements, by selecting suitable linear (or bilinear) shape functions from the Pascal triangle [1]. It is also noted that the variation can be of higher order when a higher number of nodes are used on each side of the element.However, suitable first-order shape functions were not available for element geometries with more than four sides until around the 1970s. Wachspress [2, 3] introduced a new type of shape functions based on principles of perspective geometry known as Wachspress shape functions. Linear relations within shape functions for elements with more than four nodes are obtained by using rational functions. It can be seen that the shape functions consist of complex rational functions, which requires special integration techniques to solve. Wachspress method was revisited and gained more attention around the year 2000. Meanwhile, various methods have been proposed over the years to form polygonal/polyhedral finite elements and to solve problems within polygonal/polyhedral meshes. These methods are as follows:(1)
Voronoi cell finite element method (VCFEM) and polygonal finite element based on parametric variational principle and the parametric quadratic programming method(2)
Hybrid polygonal element (HPE)(3)
Conforming polygonal finite element method based on barycentric coordinates (conforming PFEM, or PFEM)(4)
n-Sided polygonal smoothed finite element method (nSFEM)(5)
Polygonal scaled boundary finite element method (PSBFEM)(6)
Mimetic finite difference (MFD) and virtual element method (VEM)(7)
Virtual node method (VNM)(8)
Discontinuous Galerkin finite element method (DGFEM)(9)
Trefftz/Hybrid Trefftz polygonal finite element (T-FEM or HT-FEM) and Boundary element based FEM (BEM-based FEM)(10)
Hybrid stress-function (HS-F) polygonal element(11)
Base forces element method (BFEM)(12)
Other recent techniques/schemesThe methods above are briefly highlighted in the following sections.


## 2. Review on Various Techniques of Polygonal/Polyhedral Finite Elements

### 2.1. Voronoi Cell Finite Element Method (VCFEM)

Around the 1990s, Ghosh and Mukhopadhyay [4] and Ghosh and Moorthy [5] proposed a technique to model and simulate polycrystalline ferroelectrics by using polygonal elements, which are formulated by using Voronoi cell. This method is known as Voronoi cell finite element method (VCFEM). The grains in the microstructure are represented by Voronoi cells, which are generated through Voronoi tessellation. Each of these Voronoi cells (in the form of polygons with arbitrary number of sides) contains heterogeneity (in the form of void or inclusions) and is treated as a single finite element. Figure 1 shows a Voronoi cell element with heterogeneity (within the element) and boundary conditions.Figure 1
A Voronoi cell element.The matrix phase is denoted asΩm and heterogeneity phase as Ωh. The heterogeneity surface and element surface are denoted as ∂Ωh and ∂Ωe, respectively. The element boundaries consist of three types, which are prescribed traction/displacement boundary, free boundary, and interelement boundary.VCFEM combines assumptions from micromechanics theories and adaptive enhancements. This technique is computationally efficient compared to the conventional FEM (displacement based triangular or quadrilateral elements), since each polygonal grain is represented by a single finite element and no further subdivision of the domain is required. Stress functions for the interior of the element are obtained in terms of polynomial expansions of the global coordinates and these polynomial functions are formulated in such a way that they satisfy equilibrium within the element.An example of VCFEM formulation for analysis of Cosserat materials based on the parametric minimum complementary energy principle is given as [6](1)∏Ωe​=∫Ω12dσTSdσdΩ+∫ΩλTQdσdΩ-∫∂ΩdTTduds,whereΩ is the computational region of the element∂ Ω represents boundary of the elementσ is a 6-by-1 matrix containing equivalent stress componentsS is a 6-by-6 inverse of the material property matrix (elastic compliance matrix)λ is a matrix containing plastic flow parametersQ is a matrix of partial derivative of flow potential function with respect to the stressu and T are the matrices for displacement and traction components along the element boundary, respectivelys represents boundary surfaceStress functions for the interior of this element is defined by using Airy’s stress function,Φ and the Mindlin stress function, Ψ in terms of polynomial expansions of the global coordinates [6]:(2)Φ=β1x2+β2xy+β3y2+β6x3+β7x2y+β8xy2+β9y3+β13x4+β14x3y+β15x2y2+β16xy3+β17y4+⋯(3)Ψ=β4x+β5y+β10x2+β11xy+β12y2+β18x3+β19x2y+β20xy2+β21y3+⋯,where βi (i=1,2,…n) are the undetermined coefficients.VCFEM has been found to perform poorly when the heterogeneity is in the form of voids (but works well for inclusions), due to the poorly defined stress functions within the interior of the element. This problem is solved by taking into account the geometry effects through conformal mapping [7].Later, Sze and Sheng [8] included other characteristic features to the VCFEM by incorporating the electromechanical Hellinger–Reissner principle. Since then, the VCFEM has been revisited, extended, and implemented in different applications such as in analysis of heterogeneous materials [5, 9, 10], determination of the effective elastic properties/crack analysis of functionally graded materials [11, 12], analysis of microstructural Representative Volume Element (RVE) [13], simulation of the crack and analysis/failure analysis of composite materials [7, 14–19], and multiscale simulations for microstructural modelling [9]. The method is also extended to 3D [14, 20, 21].Hybrid of VCFEM with other methods such as numerical conformal mapping (NCM) method can be seen in the work by Tiwary, Hu, and Ghosh [22]. The hybrid technique is known as NCM-VCFEM and it is developed so that real micrographs of heterogeneous materials with irregular shapes can be analyzed effectively. The effectiveness is attained by using NCM such as Schwarz–Christoffel mapping to convert arbitrary/irregular shapes of micrographs into a unit circle. Another example is coupling of VCFEM with parametric variational principle and the parametric quadratic programming method [6, 23]. This approach is implemented in order to incorporate the constitutive relations of the physical phenomenon and to generalize the variational principle. The new formulations are found to be able to produce good solutions, competitive with that of conventional FEM package in ANSYS and with fewer nodes (which reduces the computational cost). Simultaneously, Zhang et al. [23] developed a polygonal finite element by using similar approach and named it as parametric variational principle based polygonal finite element method (PFEM). PFEM is similar to conventional displacement based finite elements in the sense that PFEM is applicable for macroanalysis and compatible interpolation functions are used for the entire element domain (boundaries as well as interior of the elements).VCFEM is generally suitable for micromechanical analysis and multiscale modelling [24].


### 2.2. Hybrid Polygonal Element (HPE)

On the other hand, Zhang and Katsube [25] proposed a different method to analyze micromechanical properties of heterogeneous materials. The authors used the hybrid stress element method [26] with Muskhelishvili’s complex analysis approach to formulate polygonal elements. This method is known as the hybrid polygonal element (HPE) method [27] and this technique is pursued, since the stress variation around the heterogeneity (inclusions) is not well defined in the VCFEM [28]. Knowing this, HPE has later been incorporated into VCFEM as well [20]. In the HPE method, the compatible interpolation functions are formulated along the boundaries only, while the interior of the element is represented by self-equilibrating stress field.Figure2 shows a HPE in global x-y coordinate system and its equivalent mapping to a standard ellipse in reference η-ξ coordinate system.Figure 2
Mapping of a Voronoi cell element.An example of a hybrid functional∏e(H)for a HPE with hole (void) is given as [25](4)∏eH​=12∮C1TTuds-∮C2TTuds-∮C1TTu~1ds,whereu and T are the matrices for displacement and traction components along the elementC 1 and C2 represent the element’s outer boundary with adjacent elements and the inner interface between the matrix and heterogeneity, respectivelyu ~ ( 1 ) represents the components of the specified displacements along C1s represents boundary surfaceThe interior stresses are obtained by using interpolation functions (in the form of trigonometric functions) through the following relation [25]:(5)
σ
11
σ
22
σ
12
=
1
2
μ
∑
k
=
m
b
m
u
2
C
1
-
C
2
+
C
3
-
2
D
1
+
D
2
-
D
3
-
C
1
D
1
2
C
1
+
C
2
-
C
3
-
2
D
1
-
D
2
+
D
3
C
1
-
D
1
D
2
-
D
3
C
2
-
C
3
D
1
C
1
a
~
k
a
^
k
b
~
k
b
^
k
,where C1 to C3 and D1 and D3 are the functions of k, r, and θ. mu and mb are the upper and low limits of the series. r and θ are the parameters that resulted due to the mapping of the computational/matrix region Ω to the elliptical region in reference system. Expression for C1 is shown here as an example:(6)C1k,r,θ=krk-1f1cos⁡k-1θ-f2sin⁡k-1θ,where f1 and f2 are functions in terms of r, θ, a, and b. a and b represent the major and minor semiaxes of the mapped ellipse.Similar functions are used to determine the displacements within the element.Application of HPEs for the analysis of heterogeneous media in 2D as well as 3D can be seen in works by Zhang and Katsube [25] and Kaliappan and Andreas [29]. Wachspress shape functions were not utilized in these elements (in HPE and VCFEM), due to the difficulties in integrating the complex rational functions [28]. A limitation of the VCFEM and HPE methods is that the resulting polygonal elements can contain only one irregular phase (void/inclusion) within the element. Due to this, the extended multiscale finite element method [30] has been developed to analyze mechanical behaviors of heterogeneous materials with randomly distributed polygonal microstructure. Another version of HPE for plane linear elasticity problems with better performance compared to conventional FEM is presented in [31]. The polygonal meshes are generated based on the MATLAB code PolyMesher which operates based on Voronoi diagrams.


### 2.3. Conforming Polygonal Finite Element Method Based on Barycentric Coordinates (Conforming PFEM)

Around the year of 2000, the Wachspress method gained more attention and was revisited alongside with other techniques to formulate interpolation or shape functions for polygonal elements. These techniques are known as barycentric coordinates method [32, 33] and they yield complex shape functions consisting of rational, logarithmic, and trigonometric functions. Recently, polynomial spline functions (Bernstein Bezier functions) have been proposed to be included in the barycentric coordinates method [34]. Polygonal elements which are formed by these methods are implemented in FEM and known as conforming polygonal finite element method (conforming PFEM, or simply PFEM). Some of the barycentric coordinates methods used in the formulation of conforming PFEM are inverse bilinear coordinates [35], Wachspress [32, 36–40], mean value coordinates [41, 42], harmonic coordinates [43], maximum entropy coordinates [44–46], metric coordinates [47], and natural neighbor-based coordinates (Laplace shape functions) [33]. Some of the methods such as Wachspress, mean value, and harmonic coordinates have been extended to 3D [48–52]. The above mentioned barycentric coordinate methods are described below.Inverse bilinear coordinates were developed for quadrilaterals, based on bilinear mapping of a unit square to convex quadrilaterals. Rational functions are used for the mapping and their inverses were studied to develop the inverse bilinear coordinates for quadrilaterals. Wachspress developed rational polynomial functions which can be used to produce conforming shape functions for arbitrary polygons. Meyer et al. [32] modified the Wachspress coordinates by replacing the adjoint with triangle areas and rewrote the barycentric coordinates in a simpler form. Similarly, other simplifications were carried out onto the Wachspress coordinates such as representing the Wachspress coordinates in terms of perpendicular distance between two points [32], redefining the adjoint polynomials by other means [48], and so on. Advantages of Wachspress coordinates over inverse bilinear coordinates are that the Wachspress method is applicable for arbitrary polygons and do not contain square root terms [35]. However, Wachspress’ rational shape functions do not perform well for concave polygonal elements.This shortage is avoided in mean value coordinates, which are written in terms of trigonometric functions. Mean value coordinates can be adapted to complex arbitrary polygons, especially star-shaped geometries (concave polygons). It is noted that shape functions for concave polygonal elements cannot be represented by rational polynomial functions, since convex shapes (can be represented by rational polynomial functions) cannot be mapped to concave shapes [47]. Mean value coordinates are useful and vastly applied in parameterizing triangular meshes and surface fitting [35]. These barycentric coordinates are found to be robust and applicable for concave polygons as well, even though the method does not guarantee positive functions for all the cases, since it is bounded only for star-shaped geometries [47]. On the other hand, apart from being linearly precise, harmonic and maximum-entropy coordinates are guaranteed positive for both convex and concave polygons [46]. Another suitable barycentric coordinate which satisfies the boundedness requirement for both convex and concave polygonal elements is the metric coordinate method [53]. A general framework to construct barycentric coordinates was proposed by Floater, Hormann, and Kos [54]. This framework reproduces the barycentric coordinates under various values of coefficient c in the formulation.Motivated by mesh-free method, authors of references [33, 47] developed natural neighbor-based coordinates (also known as Laplace shape functions) for arbitrary polygonal elements. The development is based on natural neighbor-based schemes within a Voronoi cell. The method was later tested for utilization in FEM, together with Wachspress, metric, and mean value coordinates. Simulation results showed that the Laplace interpolant is simpler and computationally attractive and yields more accurate results compared to the rest for convex and weakly convex polygonal elements. Nonetheless, the Laplace interpolant is not suitable for concave elements and best results for convex elements are attained through mapping of parent element [33, 47].Polygonal elements based on barycentric coordinates have been implemented in various areas such as computer graphics, animation and geometric modelling [55, 56], topology optimization [57], surface parameterization [58], geometric modelling [41], analysis of a plate with a circular hole [59], crack growth modelling [60], contact-impact problems [61], mesh generation, material fracture [62], finite elasticity problems [63], modelling of rock materials [64], and so on. Recently barycentric coordinates have been implemented in static and free vibration analyses of laminated composite plates [65], multimaterial topology optimization [66], Reissner-Mindlin plate problems [67], and transient heat conduction problems [68]. Other barycentric coordinates methods have been investigated such as Poisson coordinates [69], Green coordinates, reconstructions of Green coordinates by using Cauchy’s theorem, moving least squares coordinates, and attempt to design new methods through complex representation of real-valued barycentric coordinates [35].However, evaluation of barycentric coordinates is neither simple nor efficient compared to the conventional displacement based FEM, due to the complex functions which arise in the former techniques. Furthermore, barycentric coordinates are not efficient for assembling the stiffness matrices associated with weak solutions of Poisson equations [34].Construction of shape functions∅j based on Wachspress [47] for a polygonal element in the x-y coordinate system according to Figure 3 is given as (7)∅jx,y=wjx,y∑k=1nwkx,y,wjx,y=Api,pj,pkAp,pi,pjAp,pj,pk,whereFigure 3
Construction of shape functions based on Wachspress.A represents area enclosed by the three nodes within the bracketP i, Pj, and Pk represent a particular external/surface node of the polygonal elementP represents inner node of the polygonal elementThe expression forwjx,y in (7) can be rewritten by using the angles formed (φ and Ψ) between the nodes [47] as(8)wjx,y=2cotφj+cotΨjx-xj2+y-yj2The expression for wjx,y in (7) based on mean value coordinates [47] for a polygonal element is given as(9)wjx,y=tan⁡∝i/2+tan⁡∝j/2x-xj2+y-yj2,where ∝ is the angle that is formed within the triangular partitions at the inner node, as shown in Figure 4.Figure 4
Construction of shape functions based on mean value coordinates.The expression forwjx,y in (7) based on the concept of natural neighbors (Laplace shape functions) [47] for a polygonal element according to Figure 5 is given as(10)wjx,y=sjx,yhjx,ywhereFigure 5
Construction of shape functions based on the concept of natural neighbors.s j represents length of the Voronoi edge and hj=(x-xj)2+(y-yj)2.


### 2.4.n-Sided Polygonal Smoothed Finite Element Method (nSFEM)

Another attempt to form polygonal finite element method can be seen within the smoothed finite element method (SFEM). SFEM is formed by merging conventional FEM with meshless methods. SFEM was initially formed for quadrilateral elements. Later, Dai, Liu, and Nguyen [70] extended the four-node quadrilateral smoothed elements to arbitrary sides termed asn-sided polygonal smoothed finite elements (nSFEM) and implemented the method in solid mechanics (macrolevel).InnSFEM, the polygonal element is divided into several smoothing cells in the form of triangles, which share a common node at the center of the polygon. These triangles known as smoothing cells are then subjected to smoothing techniques onto the strain components.An example ofnSFEM element is shown in Figure 6.Figure 6
Partitioning of anSFEM element into subtriangles.Point0 in Figure 6 represents the center of the polygonal element. The displacement at this point, d-0 is taken as the average of displacement of all the external nodes, given by the following equation [71]:(11)d-0=1n∑p=1nd-p,where d-p represents displacement at a particular node p and n represents total number of nodes of the polygon. The displacement within a particular subtriangle 0-1-6 (subtriangle 1) can then be represented by [71]:(12)u-i,1=N1d-1+N2d-2+N3d-0,where Nj (j = 1, 2, 3) represents the conventional shape functions for a 3 nodes’ triangular element. Substituting (12) into (11) and simplifying gives(13)u-i,1=N1+1nN3d-1+N2+1nN3d-2+1nN3d-3+⋯+1nN3d-n-1+1nN3d-nThe shape function matrix for subtriangle 0-1-6 then becomes(14)N=N1+1nN3N2+1nN31nN3⋯1nN31nN3The strain components are smoothed according to [71](15)ε-x=∑I=1nB-Ix,yd-I(16)B-Ix,y=1Aks∫ΩksB~Ix,ydΩ,where Aks represents area of smoothing domain, B~I represents the conventional compatible strain displacement matrix, and Ωks represents the smoothing domain.There are three types of smoothing techniques applicable for these smoothing cells, which are cell, node, and edge based. These techniques are known asn-sided polygonal cell-based smoothed FEM (nCS-FEM) [70],n-sided polygonal edge-based smoothed FEM (nES-FEM) [72], andn-sided polygonal node-based smoothed FEM (nNS-FEM) [73, 74], respectively.In case ofnCS-FEM, the smoothing is carried out by integrating the gradient of displacement over the particular smoothing cell’s triangular area (for 2D as shown in Figure 6), individually [71]. When combined with the conventional FEM to obtain the displacements, the area integration is reduced to surface integration for the case of a constant smoothing function. The smoothed stiffness matrix for an element is then obtained by summing up individual stiffness matrices of all the smoothing cells within the polygonal element. Shape functions are derived for each side/edge of the smoothing cells by using the two nodes that make up the particular side/edge. These shape functions for the sides/edges should be compatible, since these sides coincide with each other to form the polygonal element. Shape functions for the interior of the smoothing cells are obtained by using other methods such as PFEM or mesh-free techniques [70]. In case ofnES-FEM, the smoothing is done based on particular edge/side (instead of cell as innCS-FEM) of the smoothing cells within the polygonal element as shown in Figure 7. The smoothing is done by performing the integration over several cell domains which are linked to the particular edge or side. Therefore, the integration domain may extend to smoothing cells of adjacent elements as well (since a particular side/edge can be shared by adjacent elements) [71]. In case ofnNS-FEM, the smoothing is done based on a particular node of the smoothing cells within the polygonal element as shown in Figure 8. Similarly, the smoothing is done by performing the integration over several cell domains which are linked to the particular node and therefore, the integration domain may also extend to smoothing cells of adjacent elements (since a particular node can be shared by adjacent elements) [71].Figure 7
Smoothing domainΩks (OCDKO) fornES-FEM.Figure 8
Smoothing domainΩks (ABCDEFA) fornNS-FEM.nCS-FEM has many advantages over the conventional FEM such as the fact that stability provides accurate results as compared to the conventional FEM. This is because the stiffness matrices ofnCS-FEM tend to be less stiff and can be applied for nearly incompressible materials by using selective integration schemes to avoid volumetric locking phenomena [70]. However, for solid mechanics,nCS-FEM is proposed to be used only for regions near the boundary or very irregular parts. This is because use of these elements for interior regions would increase the number of nodes and eventually increases the computational cost [70]. Advantage ofnNS-FEM is that it is immune from volumetric locking phenomena. Disadvantage ofnNS-FEM is that the computational time is longer compared to conventional FEM for the same number of global nodes, due to larger bandwidth of stiffness matrices. Disadvantage ofnES-FEM is that there is a tendency to overestimate or underestimate the strain energy of the model for some cases. Apart from that, similarly tonNS-FEM,nES-FEM requires more computational time compared to conventional 3-node triangular elements due to the larger bandwidth [72]. Comparison between the three types ofnSFEM is provided by Nguyen et al. [72], for solid mechanics problem. It is shown thatnES-FEM provides most accurate solution compared to the others and the stiffness/softness of the model ofnES-FEM is in between the other two techniques. Combination ofnES-FEM andnNS-FEM (termed asnES/NS-FEM) to avoid volumetric locking and to achieve faster convergence can be seen in the literature [72, 75]. Applications ofnSFEM can be seen in determination of upper bound solutions to solid mechanics problems [73], fluid-solid interaction problems [75], and new application in analysis of elastic solids subjected to torsion [76]. Recently,nSFEM has been implemented for the analysis of fluid-solid interaction (FSI) problems in viscous incompressible flows together with sliding mesh [77]. Simulation results showed that the method performs better compared to the conventional finite elements. Major advantage of thenSFEM in FSI is that it is capable of performing independent domain discretization.Generally, thenSFEM is advantageous over the conventional FEM since it produces more accurate solutions, able to tolerate volumetric locking and does not require isoparametric mapping (which enables the elements to take arbitrary shape: concave and convex forms). Apart from that, the shape functions consist of polynomials, which are easier to evaluate compared to the conforming PFEM based on barycentric coordinates. Extension of the method to 3D can be seen in literatures [78–80].


### 2.5. Polygonal Scaled Boundary FEM (PSBFEM)

Scaled boundary FEM is a semianalytical method which combines the boundary element method (BEM) and FEM. It was first introduced by Song and Wolf [81], and it was soon extended to polygonal elements and named as Polygonal Scaled Boundary FEM (PSBFEM). Example of a PSBFEM element is shown in Figure 9. PSBFEM works based on scaling center, which is located at the center of the polygonal element. The scaling center is located within the polygonal element (usually at the center) in such a way that all the boundaries/sides of the polygonal element are visible from this scaling center. Radial lines are formed from the scaling center to the outer nodes of the polygonal element, and these lines are assigned value of zero at the center (scaling node) and reach value of 1 at each node. This is accomplished through implementation of radial coordinate system. Similarly, each boundary/side of the polygonal element is assigned value of 1 to -1, through implementation of local coordinate system. The boundaries are represented by conventional numerical line elements of FEM.Figure 9
A PSBFEM element.Solution along the radial direction is obtained by analytical expressions by usingm number of shape functions, wherem represents number of nodes of the polygonal element. Transformation between the Cartesian coordinates and scaled boundary coordinates (radial coordinate system and local coordinate system) is accomplished through isoparametric mapping similar to the conventional FEM. The mapping which describes scaling of the boundary has led to the name of the method. Equation (17) [82] shows transformation of coordinate system (transformation between the scaled boundary coordinate system ξ,s and Cartesian coordinate systemx,y) for a point within the triangular subdomain 0-1-2 (as shown in Figure 9). Similar transformation is done for any point within any particular triangular subdomain.(17)x=ξNsxby=ξNsyb,where Ns=[N1sN2s]=[(1/2)1-s(1/2)1+s] represents shape functions of the line element and xb=x1x2 and yb=y1y2 represent the coordinates of the boundary nodes enclosed by the specific triangular subdomain (subdomain 0-1-2).The displacementuξ,s within a particular triangular subdomain of the element is interpolated through utilization of similar shape functions as shown in (18) below [82]:(18)uξ,s=Nsuξ,where uξ=u1(ξ)u2(ξ) represent the displacements on the lines (radial lines) passing through the scaling center,0, and nodes 1 and 2, respectively.Various applications of PSBFEM can be seen in literature such as in linear elasticity [83], crack propagation [84–88], applications within geotechnical structures [89], dynamic fracture simulation [90], analysis of cracked functionally graded materials [91], polygonal mesh creation (Song et al., 2017), elastoplastic analysis of structures [92], prediction of structural responses with randomly distributed material properties [93], simulation of crack surface contact problems [94], and analysis of mesoscale concrete samples [95]. Extension of the method to 3D can be seen in the literature [96–98]. Advantages of PSBFEM compared to BEM and conventional FEM are that, in PSBFEM, analytical solutions are achieved inside the domain, discretization of free and fixed boundaries and interfaces between different materials are not required, and the calculation of stress concentrations and intensity factors based on their definition is straightforward [82]. PSBFEM exhibits some disadvantages [82]. PSBFEM is not directly applicable for unbounded domains with strongly inclined interfaces, due to the difficulty in selecting a scaling center within the body which is visible to all sides of the domain. Some modifications are needed for these cases such as the introduction of redundant nodes for subdomain creation or by moving the boundaries upward in order to create a single scaling center (alteration of the actual physical problem). In case of time dependent problems, PSBFEM cannot be directly used to process transient excitation as opposed to BEM. Unit impulse response matrices are required and there is a need for convolution integrals which increases the computational effort. It is also found that this method is not as efficient as conventional FEM or BEM when solving problems involving smooth stress variations within bounded/enclosed domain. However, PSBFEM yields highly accurate solutions for problems involving stress singularities. The stress singularity refers to a particular point within the domain in which the stress does not converge to a specific value.PSBFEM has been found to be superior to other techniques such asnSFEM and conforming PFEM within the context of linear elasticity and the linear elastic fracture mechanics [83].


### 2.6. Mimetic Finite Difference (MFD) Method and Virtual Element Method (VEM)

One of the difficulties faced in the construction of polygonal finite elements is the development of interpolation functions which extends to the interior of the element. Beir ao da Veiga, Gyrya, Lipnikov, and Manzini [99] implemented a mimetic finite difference (MFD) method to polygonal mesh and showed that the method is efficient in solving problems involving polygonal meshes, since the method uses only the surface representation of discrete unknowns and therefore the formulations are simpler.For example, consider a heat conduction phenomenon which is governed by the following governing equation (19) [100]:(19)-divK∇u=q,whereK represents the full conductivity tensor,u represents the temperature, andq is the forcing term indicating the source of heat. The crucial step in MFD method is to mimic the essential properties of the physical and mathematical model above, which can be achieved through Green formula:(20)∫ΩK-1K∇u·F→dx=-∫ΩudivF→dx+∫∂Ωu0F→·n→dx,where F→ represents the heat flux. Subsequent step is implementation of the finite difference discretization, which involves discretization of scalar functions, vector functions, and the differential operators div and K∇. Discretization of scalar and vector functions is accomplished through introduction of degree of freedom. Example of degrees of freedom for low order MFD method is shown in Figure 10.Figure 10
Low order degrees of freedom for MFD method. The arrows represent fluxes and the center represents the temperature.An additional step is necessary, which is the discretization of the integrals. This step is required in order to approximate div andK∇ (differential operators) according to mimetic approach.The method is also shown to be useful for meshes with degenerate and nonconvex polygonal elements. Since then, MFD method has been implemented in various problems (diffusion/convection-diffusion [101–104], electromagnetic field problems [105, 106], elasticity problems, Lagrangian hydrodynamics problems [107], and solving wave equations [108]) and for modelling fluid flows [109, 110]. MFD has been extended to higher order [100, 102] and also to 3D [102, 105, 111, 112].Beir ao Da Veiga, Brezzi, Cangiani, Manzini, Marini, and Russo [113] mentioned that it was quite difficult to present MFD due to nonexistence of trial functions for the interior of the element. Therefore, the MFD has been generalized and reintroduced as virtual element method (VEM). In VEM, the unknown degrees of freedom are attached to trial functions within the interior of the polygonal domain (which do not exist earlier in MFD). Three variants of VEM are presented by Russo [114], with different number of these internal degrees of freedom. This approach is now similar to conventional PFEM. This similarity opened the possibility of coupling VEM with the conforming PFEM. The resulting hybrid method (VEM-PFEM) has been found to possess high accuracy and, most importantly, the difficulties faced in the integration of complex functions resulting from barycentric coordinates in PFEM are entirely avoided [62]. Applications of VEM [115] can be seen in plate bending problems, elasticity problems, Stokes problems, Steklov eigenvalue problems, finite strain plasticity problems [116], hyperbolic problems [117], and topology optimization problems [118]. Extension of VEM to 3D can be seen in literatures [119, 120].The virtual element space on a polygonal domain (that discretizes the problem domain) is defined as [113](21)VK,k=v∈H1K:v∂K∈Bk∂K,∆vK∈Pk-2K,where VK,k represents the finite dimensional space,K represents a generic polygonal element within the domain,k represents polynomial degree of the virtual element scheme, v represents displacement field, H1(K) represents the space containingK, ∂K represents a generic edge on the polygon, Bk∂K represents the set of polynomials of degree less than or equal tok on ∂K, and Pk-2(K) denotes the space of polynomials of degree less than or equal tok-2 onK. Fork = 1, the trial functions are linear on each edge and the inside of the element is represented by harmonic functions. Fork = 2, the trial functions on each edge are made of polynomials of degree less than or equal to 2. The inside of the element is composed of polynomials with constant Laplacian.A similarity between VCFEM, HPE, MFD, and VEM is that these methods do not require the extension of compatible interpolation functions to the interior of the element. The compatible interpolation functions are required for the element boundaries only, which simplifies the formulation and enables the formulation of elements with arbitrary number of sides/nodes. Disadvantage of MFD and VEM is that they involve complex procedures and therefore require high computational effort [121].


### 2.7. Virtual Node Method (VNM)

Another attempt to overcome the difficulty in forming compatible shape functions for polygonal FEM can be seen in the literature [122]. The authors presented a novel polygonal FEM which uses a virtual node at the center of the element to formulate shape functions consisting of simple polynomials. The method was proposed as an alternative to the conforming PFEM which suffered from inaccuracy resulting from the integration of complex functions in the stiffness matrix. First, a polygonal element with arbitrary number of nodes is formed. The polygonal domain is then divided into several triangular regions, which share a common node at the center of the element (denoted as the virtual node), as shown in Figure 11.Figure 11
Partitioning of a polygon according to VNM.Compatibility among each triangular region is fulfilled by using the conventional FEM shape functions of a 3 nodes’ triangular element. These triangular regions are then coupled together (to form the polygonal element) by representing the virtual node at the center in terms of mean least square shape functions, as shown by (22) [122, 123]:(22)φTv=ϕiTvφi+ϕjTvφj+ϕkTvφ^k,where φ represents unknown field variable at a particular point/node, φ^ represents the field variable at the virtual node, ϕi, ϕj, and ϕk are the conventional FEM shape functions of a 3 nodes’ triangular element, and Tv represents a particular triangular element within the polygonal domain. For triangle T1 (an example), Tv=T1, i=1,j=2, and k=7. Equation (22) above can be written as(23)φ=ϕiφ1+ϕjφ2+ϕkφ^7Field variable at the internal node (virtual node) φ^7 can be represented by the least mean square [123]:(24)φ^7=∑m=1nNmφm=NTφ,where n represents total number of nodes of the polygonal element, NT=pT(x,y)A-1B, p is a set of basis functions, and A and B are basis matrices. Equation (24) can be rewritten in terms of shape function for a particular node:(25)ϕmTv=ϕiTv+ϕjTv·ϕiTvφi′+ϕjTvφj′+ϕkTvNmxk,yk+ϕkTvNmφi′=1,ifm=i,0,ifm≠i;φj′=1,ifm=j,0,ifm≠j.Stiffness matrix for a particular polygonal element is obtained by summing up stiffness matrices of the individual triangular regions. Next, global stiffness matrix is obtained by summing up all the stiffness matrices of the individual polygonal elements. These are shown by (26) and (27) below:(26)ke=∑v=1nkTv(27)K=∑e=1lke,wherel represents total number of elements in the mesh. It can be seen that this method is quite similar tonCS-FEM in terms of the partitioning of the domain into triangular regions/cells, integration is performed on each triangular region/cell individually, and summing of stiffness matrices of each triangular region/cell to obtain the element stiffness matrix. The difference is that strain smoothing is not carried out in VNM. The method is later extended to higher order by Oh and Lee [124].The method is advantageous compared to compatible PFEM due to the simple polynomial shape functions (which is easier to work with). VNM is found to be efficient in adaptive computation, by using quadtree or octree mesh. The method has been extended to 3D polyhedral (VPHE) and hexahedral forms and implemented in adaptive computations as can be seen in the literature [123, 125, 126]. Recently, the method has been coupled with extended FEM (XFEM) [127, 128].


### 2.8. Discontinuous Galerkin FEM (DGFEM)

This method was proposed due to the difficulties faced in executing other polygonal FEMs such as compatible PFEM, MFD, and VEM. The difficulties arise due to the complicated shape functions in compatible PFEM and complex procedures involved in MFD and VEM. These complicated entities demand high computational effort as well [121]. DGFEM, on the other hand, does not require any conforming shape function and the method is simple. Application of DGFEM in polygonal meshes can be seen in the literature [129–135]. Extension of the method to 3D can be seen in the literature [136, 137].In DGFEM, the problem domain is discretized into several polygonal cells which represent the polygonal elements. The interpolation for the elements is carried out based on a set of monomial functions which are totally independent of the element. These interpolation functions do not comply to shape function requirements of conventional FEM and therefore they are not continuous across different elements (not compatible). Due to this, integrations are carried out on the boundaries of each cell. This step is an addition compared to the conventional FEM. The degrees of freedom are obtained based on the coefficients of the linear combination of the monomial functions [138]. Additional advantages of the method are that the adjacent elements can be of different order and the elements do not need to be conforming.Interpolation of field variablesΘ within the element is achieved through [130](28)Θ=NΘ,where Θ=Θ1Θ2⋮ΘNe represents matrix of the local degrees of freedom,Ne represents the number of polygonal elements, N=s1b1s2b2…sNebNe represents matrix containing the incompatible interpolation functionsb, and s represents the element support [130]:(29)s=1forpointwithinthedomain0forpointoutsidethedomainThe interpolation functions can be made of different kinds of functions. For instance, the functions can be polynomials which are obtained from Taylor expansion as shown below for various degrees [130]:First order:b=1xySecond order:b=1xyxyx2y23rd order:b=1xyxyx2y2x3x2yxy2y3wherex and y are the local coordinates. The interpolation functions can also be made of radial basis functions [130]:(30)b=e-r12e-r22⋯e-rn2,where ri2=x-xi2+y-yi2 and xi,yi represents coordinates of the polygon vertices as well as the mass center.The stiffness matrixK for a heat transfer phenomenon is obtained through the formula [130](31)K=∫ΩkBTBdΩ-∫ΓskNTnsTBdΓ+κ∫ΓskBTnsNdΓ+∫ΓsσNTNdΓ+∫ΓΘσNTNdΓ,where N represents discontinuity jump, B represents mean value of B, Ω represents the domain space,B represents matrix containing the differentiation of interpolation functions,k represents thermal conductivity, ns represents unit vector which is tangent to the boundary, ΓΘ represents boundary of the polygon, ΓΘ represents boundary in which the unknown variable (temperature) is prescribed, σ represents discontinuity penalization parameter, and parameter κ has the value +1, -1, or 0, depending on the scheme.


### 2.9. Trefftz/Hybrid Trefftz Polygonal Finite Element (T-FEM or HT-FEM) and Boundary Element Based FEM (BEM-Based FEM)

T-FEM utilizes two different sets of functions to approximate the solutions, one for the boundary and the other is for the interior domain. For the interior of the element, a series of homogeneous solutions of the governing equation (problem equation to be solved) is used as basis functions. These basis functions are known as T-complete set and they are not conforming across the element boundaries [139]. The boundary is represented by other independent sets of conforming basis functions. An example of T-FEM element is shown in Figure 12.Figure 12
An example of T-FEM with its shape functions utilized for the case of solid mechanics.The displacementu within a domain (interior) can be interpolated by using (32)u=uˇ+∑i=1mNici,where uˇ represents the known function on the boundary, Ni is the trial function, ci is the unknown coefficient, andm represents number of trial functions. The trial functions Ni can be generated from Muskhelishvili’s complex variable formulation [139]. The displacement u~ on the boundary of the domain can be interpolated by using (33) below:(33)u~=N~xd,where d represents vector of nodal displacements and N~x=(1-ξ)/2(1+ξ)/2 represents matrix of the corresponding one-dimensional shape functions for the boundary in terms of the local coordinate system, ξ.Continuity or boundary conditions are incorporated into the interior domain by various ways, in which one of the ways is by hybrid method known as hybrid Trefftz finite element method (HT-FEM). This method uses the conforming functions of the element boundary (also known as frame) to link the interior of the elements together [140].HT-FEM has been successfully applied in linear elasticity problems [140] and found to be able to produce more accurate results and higher convergence rates compared to the conforming PFEM with Laplace/Wachspress shape functions. One of the advantages of the HT-FEM is that elements with embedded cracks or voids can be constructed. This leads to the development of new VCFEMs for microstructural analysis. The new method is known as T-Trefftz Voronoi Cell Finite Elements (VCFEM-TTs) [141] and was soon extended to 3D [142]. Recently another novel hybrid FEM has been formulated, by coupling HT-FEM with the idea of the Method of Fundamental Solution (MFS) known as hybrid fundamental solution based FEM (HFS-FEM) [143].The advantage of HT-FEM compared to conventional FEM [144] is that this method is able to handle geometry induced singularities and stress/force concentrations efficiently without mesh adjustment. This is achieved by employing special purpose Trefftz functions that satisfy both the governing equations and boundary conditions associated with the singularities. Apart from that, general polygonal elements with curved sides can be generated and the elements are tolerant to mesh distortions. Advantages of HT-FEM compared to BEM [144] are that this method is applicable for problems involving different and heterogeneous materials, and boundary integration can be avoided when field variables are to be computed inside an element and the calculation of coefficient matrices are simpler. Disadvantage of the method is that the T-complete set for some problems are either complex or difficult to formulate. Application of the method to 3D is described by Copeland, Langer, and Pusch [139]. Soon, another version of polygonal FEM emerged known as Boundary element-based FEM (BEM-based FEM) [139, 145–148]. This method is developed based on the Trefftz function.


### 2.10. Hybrid Stress-Function (HS-F) Polygonal Element

Another method has been proposed by combining the principle of minimum complementary energy (similar principal used in VCFEM) with the Airy stress function which is known as hybrid stress-function (HS-F) element method. It was developed for quadrilaterals and triangles [149–154]. These elements are found to possess excellent performance compared to the conventional elements and especially independent of the element geometry (immune to mesh distortion). Later Zhou and Cen [155] expanded the method to polygonal elements. An example of HS-F polygonal element is shown in Figure 13. Nodes 1-6 are corner nodes and nodes 7-12 are the midside nodes. Local coordinate system ξ is shown in Figure 13 for the edge 6-7-1.Figure 13
An example of a HS-F element.The displacement along a particular element edged is given by [155](34)d=Nq,where N=NiNjNk represents the vector of shape functions for the three nodes (i,j,k) along a particular edge of the element and q=uiviujvjukvk represents the displacement vector for the nodes. The shape functions are given as [155](35)N1=-0.5ξ1-ξN2=1-ξ2N3=0.5ξ1+ξ,where ξ represents local coordinate system along each element boundary. The element stress fields are given as [155](36)σ=Sβ+σ∗,whereS is the stress solution matrix (with dimensions 3 byk) which is derived fromk number of analytical solutions of the Airy stress functions, β represents matrix (with dimensionsk by 1) of unknown stress parameters and σ∗ represents particular solution corresponding to body forces.


### 2.11. Base Forces Element Method (BFEM)

Stress based FEMs such as HPE and HS-F are not well desired for most of the engineering applications, due to the difficulties in obtaining suitable/compatible stress functions. Furthermore, it is difficult to acquire nodal displacements in stress based FEMs [156]. BFEM on the other hand was formulated based on the concept of “base forces” which was introduced by Gao [157]. This method replaces the stress functions in the stress based FEMs with base forces which are easier to obtain (obtained directly from strain energy).Conventional FEM exhibits major drawback for nonlinear analysis. The conventional FEM is not able to approximate the strain and force fields accurately since these terms are dependent on the interpolation of the displacement field (displacement shape functions). This problem is avoided in BFEM which is directly based on interpolation of the internal force fields (force shape functions) [158]. Performance of BFEM is also found to be superior than conventional FEM when analyzing large strain contact problems and the nonlinear problems. This is because the deformation field in conventional FEM is complex and discontinuous near nonlinear regions and, in the case of large strain problems, the steep gradients of deformation cannot be represented accurately by the conventional shape functions [159]. BFEM has been later extended to polygonal elements based on the complementary energy principle [160] and potential energy principle [161]. Figure 14 shows an example of a BFEM element. I,J,K,L,M,N represent the edges of the element and TI,TJ,TK,TL,TM,TN represent the force vectors which act on the element edges.Figure 14
An example of a BFEM element.The stresses corresponding to a pointI on the element can be represented as [160](37)σ=1A∑I=14TI1PI1∑I=14TI1PI2∑I=14TI2PI1∑I=14TI2PI2,where TI1,TI2 are the components of force vectors which act on the center of edgeI, PI1,PI2 are the components of position vector of pointI,andA represents area of the element.


### 2.12. Other Recent Techniques/Schemes

Recently, more new schemes have been proposed for the development of polygonal/polyhedral finite elements. They are the Compatible Discrete Operator (CDO) scheme [162, 163], Hybrid High-Order (HHO) scheme [164–166], Weak Galerkin (WG) scheme [167–171], gradient correction scheme [172], and vertex-based schemes [173].Other recent techniques/methods include analysis of polygonal carbon nanotubes reinforced composite plates by using the first-order shear deformation theory (FSDT) and the element-free IMLS-Ritz method [174], an adaptive polygonal finite element method using the techniques of cut-cell and quadtree refinement [175], new adaptive mesh generation for polygonal element [176], and ultraweak formulations for high-order polygonal finite element methods [177]. New technique for 3-dimensional polyhedral elements can be seen within the framework of the finite volume method [178]. Another new approach to form polyhedral elements is by cutting a regular hexahedral element with CAD surfaces [179].


## 2.1. Voronoi Cell Finite Element Method (VCFEM)

Around the 1990s, Ghosh and Mukhopadhyay [4] and Ghosh and Moorthy [5] proposed a technique to model and simulate polycrystalline ferroelectrics by using polygonal elements, which are formulated by using Voronoi cell. This method is known as Voronoi cell finite element method (VCFEM). The grains in the microstructure are represented by Voronoi cells, which are generated through Voronoi tessellation. Each of these Voronoi cells (in the form of polygons with arbitrary number of sides) contains heterogeneity (in the form of void or inclusions) and is treated as a single finite element. Figure 1 shows a Voronoi cell element with heterogeneity (within the element) and boundary conditions.Figure 1
A Voronoi cell element.The matrix phase is denoted asΩm and heterogeneity phase as Ωh. The heterogeneity surface and element surface are denoted as ∂Ωh and ∂Ωe, respectively. The element boundaries consist of three types, which are prescribed traction/displacement boundary, free boundary, and interelement boundary.VCFEM combines assumptions from micromechanics theories and adaptive enhancements. This technique is computationally efficient compared to the conventional FEM (displacement based triangular or quadrilateral elements), since each polygonal grain is represented by a single finite element and no further subdivision of the domain is required. Stress functions for the interior of the element are obtained in terms of polynomial expansions of the global coordinates and these polynomial functions are formulated in such a way that they satisfy equilibrium within the element.An example of VCFEM formulation for analysis of Cosserat materials based on the parametric minimum complementary energy principle is given as [6](1)∏Ωe​=∫Ω12dσTSdσdΩ+∫ΩλTQdσdΩ-∫∂ΩdTTduds,whereΩ is the computational region of the element∂ Ω represents boundary of the elementσ is a 6-by-1 matrix containing equivalent stress componentsS is a 6-by-6 inverse of the material property matrix (elastic compliance matrix)λ is a matrix containing plastic flow parametersQ is a matrix of partial derivative of flow potential function with respect to the stressu and T are the matrices for displacement and traction components along the element boundary, respectivelys represents boundary surfaceStress functions for the interior of this element is defined by using Airy’s stress function,Φ and the Mindlin stress function, Ψ in terms of polynomial expansions of the global coordinates [6]:(2)Φ=β1x2+β2xy+β3y2+β6x3+β7x2y+β8xy2+β9y3+β13x4+β14x3y+β15x2y2+β16xy3+β17y4+⋯(3)Ψ=β4x+β5y+β10x2+β11xy+β12y2+β18x3+β19x2y+β20xy2+β21y3+⋯,where βi (i=1,2,…n) are the undetermined coefficients.VCFEM has been found to perform poorly when the heterogeneity is in the form of voids (but works well for inclusions), due to the poorly defined stress functions within the interior of the element. This problem is solved by taking into account the geometry effects through conformal mapping [7].Later, Sze and Sheng [8] included other characteristic features to the VCFEM by incorporating the electromechanical Hellinger–Reissner principle. Since then, the VCFEM has been revisited, extended, and implemented in different applications such as in analysis of heterogeneous materials [5, 9, 10], determination of the effective elastic properties/crack analysis of functionally graded materials [11, 12], analysis of microstructural Representative Volume Element (RVE) [13], simulation of the crack and analysis/failure analysis of composite materials [7, 14–19], and multiscale simulations for microstructural modelling [9]. The method is also extended to 3D [14, 20, 21].Hybrid of VCFEM with other methods such as numerical conformal mapping (NCM) method can be seen in the work by Tiwary, Hu, and Ghosh [22]. The hybrid technique is known as NCM-VCFEM and it is developed so that real micrographs of heterogeneous materials with irregular shapes can be analyzed effectively. The effectiveness is attained by using NCM such as Schwarz–Christoffel mapping to convert arbitrary/irregular shapes of micrographs into a unit circle. Another example is coupling of VCFEM with parametric variational principle and the parametric quadratic programming method [6, 23]. This approach is implemented in order to incorporate the constitutive relations of the physical phenomenon and to generalize the variational principle. The new formulations are found to be able to produce good solutions, competitive with that of conventional FEM package in ANSYS and with fewer nodes (which reduces the computational cost). Simultaneously, Zhang et al. [23] developed a polygonal finite element by using similar approach and named it as parametric variational principle based polygonal finite element method (PFEM). PFEM is similar to conventional displacement based finite elements in the sense that PFEM is applicable for macroanalysis and compatible interpolation functions are used for the entire element domain (boundaries as well as interior of the elements).VCFEM is generally suitable for micromechanical analysis and multiscale modelling [24].


## 2.2. Hybrid Polygonal Element (HPE)

On the other hand, Zhang and Katsube [25] proposed a different method to analyze micromechanical properties of heterogeneous materials. The authors used the hybrid stress element method [26] with Muskhelishvili’s complex analysis approach to formulate polygonal elements. This method is known as the hybrid polygonal element (HPE) method [27] and this technique is pursued, since the stress variation around the heterogeneity (inclusions) is not well defined in the VCFEM [28]. Knowing this, HPE has later been incorporated into VCFEM as well [20]. In the HPE method, the compatible interpolation functions are formulated along the boundaries only, while the interior of the element is represented by self-equilibrating stress field.Figure2 shows a HPE in global x-y coordinate system and its equivalent mapping to a standard ellipse in reference η-ξ coordinate system.Figure 2
Mapping of a Voronoi cell element.An example of a hybrid functional∏e(H)for a HPE with hole (void) is given as [25](4)∏eH​=12∮C1TTuds-∮C2TTuds-∮C1TTu~1ds,whereu and T are the matrices for displacement and traction components along the elementC 1 and C2 represent the element’s outer boundary with adjacent elements and the inner interface between the matrix and heterogeneity, respectivelyu ~ ( 1 ) represents the components of the specified displacements along C1s represents boundary surfaceThe interior stresses are obtained by using interpolation functions (in the form of trigonometric functions) through the following relation [25]:(5)
σ
11
σ
22
σ
12
=
1
2
μ
∑
k
=
m
b
m
u
2
C
1
-
C
2
+
C
3
-
2
D
1
+
D
2
-
D
3
-
C
1
D
1
2
C
1
+
C
2
-
C
3
-
2
D
1
-
D
2
+
D
3
C
1
-
D
1
D
2
-
D
3
C
2
-
C
3
D
1
C
1
a
~
k
a
^
k
b
~
k
b
^
k
,where C1 to C3 and D1 and D3 are the functions of k, r, and θ. mu and mb are the upper and low limits of the series. r and θ are the parameters that resulted due to the mapping of the computational/matrix region Ω to the elliptical region in reference system. Expression for C1 is shown here as an example:(6)C1k,r,θ=krk-1f1cos⁡k-1θ-f2sin⁡k-1θ,where f1 and f2 are functions in terms of r, θ, a, and b. a and b represent the major and minor semiaxes of the mapped ellipse.Similar functions are used to determine the displacements within the element.Application of HPEs for the analysis of heterogeneous media in 2D as well as 3D can be seen in works by Zhang and Katsube [25] and Kaliappan and Andreas [29]. Wachspress shape functions were not utilized in these elements (in HPE and VCFEM), due to the difficulties in integrating the complex rational functions [28]. A limitation of the VCFEM and HPE methods is that the resulting polygonal elements can contain only one irregular phase (void/inclusion) within the element. Due to this, the extended multiscale finite element method [30] has been developed to analyze mechanical behaviors of heterogeneous materials with randomly distributed polygonal microstructure. Another version of HPE for plane linear elasticity problems with better performance compared to conventional FEM is presented in [31]. The polygonal meshes are generated based on the MATLAB code PolyMesher which operates based on Voronoi diagrams.


## 2.3. Conforming Polygonal Finite Element Method Based on Barycentric Coordinates (Conforming PFEM)

Around the year of 2000, the Wachspress method gained more attention and was revisited alongside with other techniques to formulate interpolation or shape functions for polygonal elements. These techniques are known as barycentric coordinates method [32, 33] and they yield complex shape functions consisting of rational, logarithmic, and trigonometric functions. Recently, polynomial spline functions (Bernstein Bezier functions) have been proposed to be included in the barycentric coordinates method [34]. Polygonal elements which are formed by these methods are implemented in FEM and known as conforming polygonal finite element method (conforming PFEM, or simply PFEM). Some of the barycentric coordinates methods used in the formulation of conforming PFEM are inverse bilinear coordinates [35], Wachspress [32, 36–40], mean value coordinates [41, 42], harmonic coordinates [43], maximum entropy coordinates [44–46], metric coordinates [47], and natural neighbor-based coordinates (Laplace shape functions) [33]. Some of the methods such as Wachspress, mean value, and harmonic coordinates have been extended to 3D [48–52]. The above mentioned barycentric coordinate methods are described below.Inverse bilinear coordinates were developed for quadrilaterals, based on bilinear mapping of a unit square to convex quadrilaterals. Rational functions are used for the mapping and their inverses were studied to develop the inverse bilinear coordinates for quadrilaterals. Wachspress developed rational polynomial functions which can be used to produce conforming shape functions for arbitrary polygons. Meyer et al. [32] modified the Wachspress coordinates by replacing the adjoint with triangle areas and rewrote the barycentric coordinates in a simpler form. Similarly, other simplifications were carried out onto the Wachspress coordinates such as representing the Wachspress coordinates in terms of perpendicular distance between two points [32], redefining the adjoint polynomials by other means [48], and so on. Advantages of Wachspress coordinates over inverse bilinear coordinates are that the Wachspress method is applicable for arbitrary polygons and do not contain square root terms [35]. However, Wachspress’ rational shape functions do not perform well for concave polygonal elements.This shortage is avoided in mean value coordinates, which are written in terms of trigonometric functions. Mean value coordinates can be adapted to complex arbitrary polygons, especially star-shaped geometries (concave polygons). It is noted that shape functions for concave polygonal elements cannot be represented by rational polynomial functions, since convex shapes (can be represented by rational polynomial functions) cannot be mapped to concave shapes [47]. Mean value coordinates are useful and vastly applied in parameterizing triangular meshes and surface fitting [35]. These barycentric coordinates are found to be robust and applicable for concave polygons as well, even though the method does not guarantee positive functions for all the cases, since it is bounded only for star-shaped geometries [47]. On the other hand, apart from being linearly precise, harmonic and maximum-entropy coordinates are guaranteed positive for both convex and concave polygons [46]. Another suitable barycentric coordinate which satisfies the boundedness requirement for both convex and concave polygonal elements is the metric coordinate method [53]. A general framework to construct barycentric coordinates was proposed by Floater, Hormann, and Kos [54]. This framework reproduces the barycentric coordinates under various values of coefficient c in the formulation.Motivated by mesh-free method, authors of references [33, 47] developed natural neighbor-based coordinates (also known as Laplace shape functions) for arbitrary polygonal elements. The development is based on natural neighbor-based schemes within a Voronoi cell. The method was later tested for utilization in FEM, together with Wachspress, metric, and mean value coordinates. Simulation results showed that the Laplace interpolant is simpler and computationally attractive and yields more accurate results compared to the rest for convex and weakly convex polygonal elements. Nonetheless, the Laplace interpolant is not suitable for concave elements and best results for convex elements are attained through mapping of parent element [33, 47].Polygonal elements based on barycentric coordinates have been implemented in various areas such as computer graphics, animation and geometric modelling [55, 56], topology optimization [57], surface parameterization [58], geometric modelling [41], analysis of a plate with a circular hole [59], crack growth modelling [60], contact-impact problems [61], mesh generation, material fracture [62], finite elasticity problems [63], modelling of rock materials [64], and so on. Recently barycentric coordinates have been implemented in static and free vibration analyses of laminated composite plates [65], multimaterial topology optimization [66], Reissner-Mindlin plate problems [67], and transient heat conduction problems [68]. Other barycentric coordinates methods have been investigated such as Poisson coordinates [69], Green coordinates, reconstructions of Green coordinates by using Cauchy’s theorem, moving least squares coordinates, and attempt to design new methods through complex representation of real-valued barycentric coordinates [35].However, evaluation of barycentric coordinates is neither simple nor efficient compared to the conventional displacement based FEM, due to the complex functions which arise in the former techniques. Furthermore, barycentric coordinates are not efficient for assembling the stiffness matrices associated with weak solutions of Poisson equations [34].Construction of shape functions∅j based on Wachspress [47] for a polygonal element in the x-y coordinate system according to Figure 3 is given as (7)∅jx,y=wjx,y∑k=1nwkx,y,wjx,y=Api,pj,pkAp,pi,pjAp,pj,pk,whereFigure 3
Construction of shape functions based on Wachspress.A represents area enclosed by the three nodes within the bracketP i, Pj, and Pk represent a particular external/surface node of the polygonal elementP represents inner node of the polygonal elementThe expression forwjx,y in (7) can be rewritten by using the angles formed (φ and Ψ) between the nodes [47] as(8)wjx,y=2cotφj+cotΨjx-xj2+y-yj2The expression for wjx,y in (7) based on mean value coordinates [47] for a polygonal element is given as(9)wjx,y=tan⁡∝i/2+tan⁡∝j/2x-xj2+y-yj2,where ∝ is the angle that is formed within the triangular partitions at the inner node, as shown in Figure 4.Figure 4
Construction of shape functions based on mean value coordinates.The expression forwjx,y in (7) based on the concept of natural neighbors (Laplace shape functions) [47] for a polygonal element according to Figure 5 is given as(10)wjx,y=sjx,yhjx,ywhereFigure 5
Construction of shape functions based on the concept of natural neighbors.s j represents length of the Voronoi edge and hj=(x-xj)2+(y-yj)2.


## 2.4.n-Sided Polygonal Smoothed Finite Element Method (nSFEM)

Another attempt to form polygonal finite element method can be seen within the smoothed finite element method (SFEM). SFEM is formed by merging conventional FEM with meshless methods. SFEM was initially formed for quadrilateral elements. Later, Dai, Liu, and Nguyen [70] extended the four-node quadrilateral smoothed elements to arbitrary sides termed asn-sided polygonal smoothed finite elements (nSFEM) and implemented the method in solid mechanics (macrolevel).InnSFEM, the polygonal element is divided into several smoothing cells in the form of triangles, which share a common node at the center of the polygon. These triangles known as smoothing cells are then subjected to smoothing techniques onto the strain components.An example ofnSFEM element is shown in Figure 6.Figure 6
Partitioning of anSFEM element into subtriangles.Point0 in Figure 6 represents the center of the polygonal element. The displacement at this point, d-0 is taken as the average of displacement of all the external nodes, given by the following equation [71]:(11)d-0=1n∑p=1nd-p,where d-p represents displacement at a particular node p and n represents total number of nodes of the polygon. The displacement within a particular subtriangle 0-1-6 (subtriangle 1) can then be represented by [71]:(12)u-i,1=N1d-1+N2d-2+N3d-0,where Nj (j = 1, 2, 3) represents the conventional shape functions for a 3 nodes’ triangular element. Substituting (12) into (11) and simplifying gives(13)u-i,1=N1+1nN3d-1+N2+1nN3d-2+1nN3d-3+⋯+1nN3d-n-1+1nN3d-nThe shape function matrix for subtriangle 0-1-6 then becomes(14)N=N1+1nN3N2+1nN31nN3⋯1nN31nN3The strain components are smoothed according to [71](15)ε-x=∑I=1nB-Ix,yd-I(16)B-Ix,y=1Aks∫ΩksB~Ix,ydΩ,where Aks represents area of smoothing domain, B~I represents the conventional compatible strain displacement matrix, and Ωks represents the smoothing domain.There are three types of smoothing techniques applicable for these smoothing cells, which are cell, node, and edge based. These techniques are known asn-sided polygonal cell-based smoothed FEM (nCS-FEM) [70],n-sided polygonal edge-based smoothed FEM (nES-FEM) [72], andn-sided polygonal node-based smoothed FEM (nNS-FEM) [73, 74], respectively.In case ofnCS-FEM, the smoothing is carried out by integrating the gradient of displacement over the particular smoothing cell’s triangular area (for 2D as shown in Figure 6), individually [71]. When combined with the conventional FEM to obtain the displacements, the area integration is reduced to surface integration for the case of a constant smoothing function. The smoothed stiffness matrix for an element is then obtained by summing up individual stiffness matrices of all the smoothing cells within the polygonal element. Shape functions are derived for each side/edge of the smoothing cells by using the two nodes that make up the particular side/edge. These shape functions for the sides/edges should be compatible, since these sides coincide with each other to form the polygonal element. Shape functions for the interior of the smoothing cells are obtained by using other methods such as PFEM or mesh-free techniques [70]. In case ofnES-FEM, the smoothing is done based on particular edge/side (instead of cell as innCS-FEM) of the smoothing cells within the polygonal element as shown in Figure 7. The smoothing is done by performing the integration over several cell domains which are linked to the particular edge or side. Therefore, the integration domain may extend to smoothing cells of adjacent elements as well (since a particular side/edge can be shared by adjacent elements) [71]. In case ofnNS-FEM, the smoothing is done based on a particular node of the smoothing cells within the polygonal element as shown in Figure 8. Similarly, the smoothing is done by performing the integration over several cell domains which are linked to the particular node and therefore, the integration domain may also extend to smoothing cells of adjacent elements (since a particular node can be shared by adjacent elements) [71].Figure 7
Smoothing domainΩks (OCDKO) fornES-FEM.Figure 8
Smoothing domainΩks (ABCDEFA) fornNS-FEM.nCS-FEM has many advantages over the conventional FEM such as the fact that stability provides accurate results as compared to the conventional FEM. This is because the stiffness matrices ofnCS-FEM tend to be less stiff and can be applied for nearly incompressible materials by using selective integration schemes to avoid volumetric locking phenomena [70]. However, for solid mechanics,nCS-FEM is proposed to be used only for regions near the boundary or very irregular parts. This is because use of these elements for interior regions would increase the number of nodes and eventually increases the computational cost [70]. Advantage ofnNS-FEM is that it is immune from volumetric locking phenomena. Disadvantage ofnNS-FEM is that the computational time is longer compared to conventional FEM for the same number of global nodes, due to larger bandwidth of stiffness matrices. Disadvantage ofnES-FEM is that there is a tendency to overestimate or underestimate the strain energy of the model for some cases. Apart from that, similarly tonNS-FEM,nES-FEM requires more computational time compared to conventional 3-node triangular elements due to the larger bandwidth [72]. Comparison between the three types ofnSFEM is provided by Nguyen et al. [72], for solid mechanics problem. It is shown thatnES-FEM provides most accurate solution compared to the others and the stiffness/softness of the model ofnES-FEM is in between the other two techniques. Combination ofnES-FEM andnNS-FEM (termed asnES/NS-FEM) to avoid volumetric locking and to achieve faster convergence can be seen in the literature [72, 75]. Applications ofnSFEM can be seen in determination of upper bound solutions to solid mechanics problems [73], fluid-solid interaction problems [75], and new application in analysis of elastic solids subjected to torsion [76]. Recently,nSFEM has been implemented for the analysis of fluid-solid interaction (FSI) problems in viscous incompressible flows together with sliding mesh [77]. Simulation results showed that the method performs better compared to the conventional finite elements. Major advantage of thenSFEM in FSI is that it is capable of performing independent domain discretization.Generally, thenSFEM is advantageous over the conventional FEM since it produces more accurate solutions, able to tolerate volumetric locking and does not require isoparametric mapping (which enables the elements to take arbitrary shape: concave and convex forms). Apart from that, the shape functions consist of polynomials, which are easier to evaluate compared to the conforming PFEM based on barycentric coordinates. Extension of the method to 3D can be seen in literatures [78–80].


## 2.5. Polygonal Scaled Boundary FEM (PSBFEM)

Scaled boundary FEM is a semianalytical method which combines the boundary element method (BEM) and FEM. It was first introduced by Song and Wolf [81], and it was soon extended to polygonal elements and named as Polygonal Scaled Boundary FEM (PSBFEM). Example of a PSBFEM element is shown in Figure 9. PSBFEM works based on scaling center, which is located at the center of the polygonal element. The scaling center is located within the polygonal element (usually at the center) in such a way that all the boundaries/sides of the polygonal element are visible from this scaling center. Radial lines are formed from the scaling center to the outer nodes of the polygonal element, and these lines are assigned value of zero at the center (scaling node) and reach value of 1 at each node. This is accomplished through implementation of radial coordinate system. Similarly, each boundary/side of the polygonal element is assigned value of 1 to -1, through implementation of local coordinate system. The boundaries are represented by conventional numerical line elements of FEM.Figure 9
A PSBFEM element.Solution along the radial direction is obtained by analytical expressions by usingm number of shape functions, wherem represents number of nodes of the polygonal element. Transformation between the Cartesian coordinates and scaled boundary coordinates (radial coordinate system and local coordinate system) is accomplished through isoparametric mapping similar to the conventional FEM. The mapping which describes scaling of the boundary has led to the name of the method. Equation (17) [82] shows transformation of coordinate system (transformation between the scaled boundary coordinate system ξ,s and Cartesian coordinate systemx,y) for a point within the triangular subdomain 0-1-2 (as shown in Figure 9). Similar transformation is done for any point within any particular triangular subdomain.(17)x=ξNsxby=ξNsyb,where Ns=[N1sN2s]=[(1/2)1-s(1/2)1+s] represents shape functions of the line element and xb=x1x2 and yb=y1y2 represent the coordinates of the boundary nodes enclosed by the specific triangular subdomain (subdomain 0-1-2).The displacementuξ,s within a particular triangular subdomain of the element is interpolated through utilization of similar shape functions as shown in (18) below [82]:(18)uξ,s=Nsuξ,where uξ=u1(ξ)u2(ξ) represent the displacements on the lines (radial lines) passing through the scaling center,0, and nodes 1 and 2, respectively.Various applications of PSBFEM can be seen in literature such as in linear elasticity [83], crack propagation [84–88], applications within geotechnical structures [89], dynamic fracture simulation [90], analysis of cracked functionally graded materials [91], polygonal mesh creation (Song et al., 2017), elastoplastic analysis of structures [92], prediction of structural responses with randomly distributed material properties [93], simulation of crack surface contact problems [94], and analysis of mesoscale concrete samples [95]. Extension of the method to 3D can be seen in the literature [96–98]. Advantages of PSBFEM compared to BEM and conventional FEM are that, in PSBFEM, analytical solutions are achieved inside the domain, discretization of free and fixed boundaries and interfaces between different materials are not required, and the calculation of stress concentrations and intensity factors based on their definition is straightforward [82]. PSBFEM exhibits some disadvantages [82]. PSBFEM is not directly applicable for unbounded domains with strongly inclined interfaces, due to the difficulty in selecting a scaling center within the body which is visible to all sides of the domain. Some modifications are needed for these cases such as the introduction of redundant nodes for subdomain creation or by moving the boundaries upward in order to create a single scaling center (alteration of the actual physical problem). In case of time dependent problems, PSBFEM cannot be directly used to process transient excitation as opposed to BEM. Unit impulse response matrices are required and there is a need for convolution integrals which increases the computational effort. It is also found that this method is not as efficient as conventional FEM or BEM when solving problems involving smooth stress variations within bounded/enclosed domain. However, PSBFEM yields highly accurate solutions for problems involving stress singularities. The stress singularity refers to a particular point within the domain in which the stress does not converge to a specific value.PSBFEM has been found to be superior to other techniques such asnSFEM and conforming PFEM within the context of linear elasticity and the linear elastic fracture mechanics [83].


## 2.6. Mimetic Finite Difference (MFD) Method and Virtual Element Method (VEM)

One of the difficulties faced in the construction of polygonal finite elements is the development of interpolation functions which extends to the interior of the element. Beir ao da Veiga, Gyrya, Lipnikov, and Manzini [99] implemented a mimetic finite difference (MFD) method to polygonal mesh and showed that the method is efficient in solving problems involving polygonal meshes, since the method uses only the surface representation of discrete unknowns and therefore the formulations are simpler.For example, consider a heat conduction phenomenon which is governed by the following governing equation (19) [100]:(19)-divK∇u=q,whereK represents the full conductivity tensor,u represents the temperature, andq is the forcing term indicating the source of heat. The crucial step in MFD method is to mimic the essential properties of the physical and mathematical model above, which can be achieved through Green formula:(20)∫ΩK-1K∇u·F→dx=-∫ΩudivF→dx+∫∂Ωu0F→·n→dx,where F→ represents the heat flux. Subsequent step is implementation of the finite difference discretization, which involves discretization of scalar functions, vector functions, and the differential operators div and K∇. Discretization of scalar and vector functions is accomplished through introduction of degree of freedom. Example of degrees of freedom for low order MFD method is shown in Figure 10.Figure 10
Low order degrees of freedom for MFD method. The arrows represent fluxes and the center represents the temperature.An additional step is necessary, which is the discretization of the integrals. This step is required in order to approximate div andK∇ (differential operators) according to mimetic approach.The method is also shown to be useful for meshes with degenerate and nonconvex polygonal elements. Since then, MFD method has been implemented in various problems (diffusion/convection-diffusion [101–104], electromagnetic field problems [105, 106], elasticity problems, Lagrangian hydrodynamics problems [107], and solving wave equations [108]) and for modelling fluid flows [109, 110]. MFD has been extended to higher order [100, 102] and also to 3D [102, 105, 111, 112].Beir ao Da Veiga, Brezzi, Cangiani, Manzini, Marini, and Russo [113] mentioned that it was quite difficult to present MFD due to nonexistence of trial functions for the interior of the element. Therefore, the MFD has been generalized and reintroduced as virtual element method (VEM). In VEM, the unknown degrees of freedom are attached to trial functions within the interior of the polygonal domain (which do not exist earlier in MFD). Three variants of VEM are presented by Russo [114], with different number of these internal degrees of freedom. This approach is now similar to conventional PFEM. This similarity opened the possibility of coupling VEM with the conforming PFEM. The resulting hybrid method (VEM-PFEM) has been found to possess high accuracy and, most importantly, the difficulties faced in the integration of complex functions resulting from barycentric coordinates in PFEM are entirely avoided [62]. Applications of VEM [115] can be seen in plate bending problems, elasticity problems, Stokes problems, Steklov eigenvalue problems, finite strain plasticity problems [116], hyperbolic problems [117], and topology optimization problems [118]. Extension of VEM to 3D can be seen in literatures [119, 120].The virtual element space on a polygonal domain (that discretizes the problem domain) is defined as [113](21)VK,k=v∈H1K:v∂K∈Bk∂K,∆vK∈Pk-2K,where VK,k represents the finite dimensional space,K represents a generic polygonal element within the domain,k represents polynomial degree of the virtual element scheme, v represents displacement field, H1(K) represents the space containingK, ∂K represents a generic edge on the polygon, Bk∂K represents the set of polynomials of degree less than or equal tok on ∂K, and Pk-2(K) denotes the space of polynomials of degree less than or equal tok-2 onK. Fork = 1, the trial functions are linear on each edge and the inside of the element is represented by harmonic functions. Fork = 2, the trial functions on each edge are made of polynomials of degree less than or equal to 2. The inside of the element is composed of polynomials with constant Laplacian.A similarity between VCFEM, HPE, MFD, and VEM is that these methods do not require the extension of compatible interpolation functions to the interior of the element. The compatible interpolation functions are required for the element boundaries only, which simplifies the formulation and enables the formulation of elements with arbitrary number of sides/nodes. Disadvantage of MFD and VEM is that they involve complex procedures and therefore require high computational effort [121].


## 2.7. Virtual Node Method (VNM)

Another attempt to overcome the difficulty in forming compatible shape functions for polygonal FEM can be seen in the literature [122]. The authors presented a novel polygonal FEM which uses a virtual node at the center of the element to formulate shape functions consisting of simple polynomials. The method was proposed as an alternative to the conforming PFEM which suffered from inaccuracy resulting from the integration of complex functions in the stiffness matrix. First, a polygonal element with arbitrary number of nodes is formed. The polygonal domain is then divided into several triangular regions, which share a common node at the center of the element (denoted as the virtual node), as shown in Figure 11.Figure 11
Partitioning of a polygon according to VNM.Compatibility among each triangular region is fulfilled by using the conventional FEM shape functions of a 3 nodes’ triangular element. These triangular regions are then coupled together (to form the polygonal element) by representing the virtual node at the center in terms of mean least square shape functions, as shown by (22) [122, 123]:(22)φTv=ϕiTvφi+ϕjTvφj+ϕkTvφ^k,where φ represents unknown field variable at a particular point/node, φ^ represents the field variable at the virtual node, ϕi, ϕj, and ϕk are the conventional FEM shape functions of a 3 nodes’ triangular element, and Tv represents a particular triangular element within the polygonal domain. For triangle T1 (an example), Tv=T1, i=1,j=2, and k=7. Equation (22) above can be written as(23)φ=ϕiφ1+ϕjφ2+ϕkφ^7Field variable at the internal node (virtual node) φ^7 can be represented by the least mean square [123]:(24)φ^7=∑m=1nNmφm=NTφ,where n represents total number of nodes of the polygonal element, NT=pT(x,y)A-1B, p is a set of basis functions, and A and B are basis matrices. Equation (24) can be rewritten in terms of shape function for a particular node:(25)ϕmTv=ϕiTv+ϕjTv·ϕiTvφi′+ϕjTvφj′+ϕkTvNmxk,yk+ϕkTvNmφi′=1,ifm=i,0,ifm≠i;φj′=1,ifm=j,0,ifm≠j.Stiffness matrix for a particular polygonal element is obtained by summing up stiffness matrices of the individual triangular regions. Next, global stiffness matrix is obtained by summing up all the stiffness matrices of the individual polygonal elements. These are shown by (26) and (27) below:(26)ke=∑v=1nkTv(27)K=∑e=1lke,wherel represents total number of elements in the mesh. It can be seen that this method is quite similar tonCS-FEM in terms of the partitioning of the domain into triangular regions/cells, integration is performed on each triangular region/cell individually, and summing of stiffness matrices of each triangular region/cell to obtain the element stiffness matrix. The difference is that strain smoothing is not carried out in VNM. The method is later extended to higher order by Oh and Lee [124].The method is advantageous compared to compatible PFEM due to the simple polynomial shape functions (which is easier to work with). VNM is found to be efficient in adaptive computation, by using quadtree or octree mesh. The method has been extended to 3D polyhedral (VPHE) and hexahedral forms and implemented in adaptive computations as can be seen in the literature [123, 125, 126]. Recently, the method has been coupled with extended FEM (XFEM) [127, 128].


## 2.8. Discontinuous Galerkin FEM (DGFEM)

This method was proposed due to the difficulties faced in executing other polygonal FEMs such as compatible PFEM, MFD, and VEM. The difficulties arise due to the complicated shape functions in compatible PFEM and complex procedures involved in MFD and VEM. These complicated entities demand high computational effort as well [121]. DGFEM, on the other hand, does not require any conforming shape function and the method is simple. Application of DGFEM in polygonal meshes can be seen in the literature [129–135]. Extension of the method to 3D can be seen in the literature [136, 137].In DGFEM, the problem domain is discretized into several polygonal cells which represent the polygonal elements. The interpolation for the elements is carried out based on a set of monomial functions which are totally independent of the element. These interpolation functions do not comply to shape function requirements of conventional FEM and therefore they are not continuous across different elements (not compatible). Due to this, integrations are carried out on the boundaries of each cell. This step is an addition compared to the conventional FEM. The degrees of freedom are obtained based on the coefficients of the linear combination of the monomial functions [138]. Additional advantages of the method are that the adjacent elements can be of different order and the elements do not need to be conforming.Interpolation of field variablesΘ within the element is achieved through [130](28)Θ=NΘ,where Θ=Θ1Θ2⋮ΘNe represents matrix of the local degrees of freedom,Ne represents the number of polygonal elements, N=s1b1s2b2…sNebNe represents matrix containing the incompatible interpolation functionsb, and s represents the element support [130]:(29)s=1forpointwithinthedomain0forpointoutsidethedomainThe interpolation functions can be made of different kinds of functions. For instance, the functions can be polynomials which are obtained from Taylor expansion as shown below for various degrees [130]:First order:b=1xySecond order:b=1xyxyx2y23rd order:b=1xyxyx2y2x3x2yxy2y3wherex and y are the local coordinates. The interpolation functions can also be made of radial basis functions [130]:(30)b=e-r12e-r22⋯e-rn2,where ri2=x-xi2+y-yi2 and xi,yi represents coordinates of the polygon vertices as well as the mass center.The stiffness matrixK for a heat transfer phenomenon is obtained through the formula [130](31)K=∫ΩkBTBdΩ-∫ΓskNTnsTBdΓ+κ∫ΓskBTnsNdΓ+∫ΓsσNTNdΓ+∫ΓΘσNTNdΓ,where N represents discontinuity jump, B represents mean value of B, Ω represents the domain space,B represents matrix containing the differentiation of interpolation functions,k represents thermal conductivity, ns represents unit vector which is tangent to the boundary, ΓΘ represents boundary of the polygon, ΓΘ represents boundary in which the unknown variable (temperature) is prescribed, σ represents discontinuity penalization parameter, and parameter κ has the value +1, -1, or 0, depending on the scheme.


## 2.9. Trefftz/Hybrid Trefftz Polygonal Finite Element (T-FEM or HT-FEM) and Boundary Element Based FEM (BEM-Based FEM)

T-FEM utilizes two different sets of functions to approximate the solutions, one for the boundary and the other is for the interior domain. For the interior of the element, a series of homogeneous solutions of the governing equation (problem equation to be solved) is used as basis functions. These basis functions are known as T-complete set and they are not conforming across the element boundaries [139]. The boundary is represented by other independent sets of conforming basis functions. An example of T-FEM element is shown in Figure 12.Figure 12
An example of T-FEM with its shape functions utilized for the case of solid mechanics.The displacementu within a domain (interior) can be interpolated by using (32)u=uˇ+∑i=1mNici,where uˇ represents the known function on the boundary, Ni is the trial function, ci is the unknown coefficient, andm represents number of trial functions. The trial functions Ni can be generated from Muskhelishvili’s complex variable formulation [139]. The displacement u~ on the boundary of the domain can be interpolated by using (33) below:(33)u~=N~xd,where d represents vector of nodal displacements and N~x=(1-ξ)/2(1+ξ)/2 represents matrix of the corresponding one-dimensional shape functions for the boundary in terms of the local coordinate system, ξ.Continuity or boundary conditions are incorporated into the interior domain by various ways, in which one of the ways is by hybrid method known as hybrid Trefftz finite element method (HT-FEM). This method uses the conforming functions of the element boundary (also known as frame) to link the interior of the elements together [140].HT-FEM has been successfully applied in linear elasticity problems [140] and found to be able to produce more accurate results and higher convergence rates compared to the conforming PFEM with Laplace/Wachspress shape functions. One of the advantages of the HT-FEM is that elements with embedded cracks or voids can be constructed. This leads to the development of new VCFEMs for microstructural analysis. The new method is known as T-Trefftz Voronoi Cell Finite Elements (VCFEM-TTs) [141] and was soon extended to 3D [142]. Recently another novel hybrid FEM has been formulated, by coupling HT-FEM with the idea of the Method of Fundamental Solution (MFS) known as hybrid fundamental solution based FEM (HFS-FEM) [143].The advantage of HT-FEM compared to conventional FEM [144] is that this method is able to handle geometry induced singularities and stress/force concentrations efficiently without mesh adjustment. This is achieved by employing special purpose Trefftz functions that satisfy both the governing equations and boundary conditions associated with the singularities. Apart from that, general polygonal elements with curved sides can be generated and the elements are tolerant to mesh distortions. Advantages of HT-FEM compared to BEM [144] are that this method is applicable for problems involving different and heterogeneous materials, and boundary integration can be avoided when field variables are to be computed inside an element and the calculation of coefficient matrices are simpler. Disadvantage of the method is that the T-complete set for some problems are either complex or difficult to formulate. Application of the method to 3D is described by Copeland, Langer, and Pusch [139]. Soon, another version of polygonal FEM emerged known as Boundary element-based FEM (BEM-based FEM) [139, 145–148]. This method is developed based on the Trefftz function.


## 2.10. Hybrid Stress-Function (HS-F) Polygonal Element

Another method has been proposed by combining the principle of minimum complementary energy (similar principal used in VCFEM) with the Airy stress function which is known as hybrid stress-function (HS-F) element method. It was developed for quadrilaterals and triangles [149–154]. These elements are found to possess excellent performance compared to the conventional elements and especially independent of the element geometry (immune to mesh distortion). Later Zhou and Cen [155] expanded the method to polygonal elements. An example of HS-F polygonal element is shown in Figure 13. Nodes 1-6 are corner nodes and nodes 7-12 are the midside nodes. Local coordinate system ξ is shown in Figure 13 for the edge 6-7-1.Figure 13
An example of a HS-F element.The displacement along a particular element edged is given by [155](34)d=Nq,where N=NiNjNk represents the vector of shape functions for the three nodes (i,j,k) along a particular edge of the element and q=uiviujvjukvk represents the displacement vector for the nodes. The shape functions are given as [155](35)N1=-0.5ξ1-ξN2=1-ξ2N3=0.5ξ1+ξ,where ξ represents local coordinate system along each element boundary. The element stress fields are given as [155](36)σ=Sβ+σ∗,whereS is the stress solution matrix (with dimensions 3 byk) which is derived fromk number of analytical solutions of the Airy stress functions, β represents matrix (with dimensionsk by 1) of unknown stress parameters and σ∗ represents particular solution corresponding to body forces.


## 2.11. Base Forces Element Method (BFEM)

Stress based FEMs such as HPE and HS-F are not well desired for most of the engineering applications, due to the difficulties in obtaining suitable/compatible stress functions. Furthermore, it is difficult to acquire nodal displacements in stress based FEMs [156]. BFEM on the other hand was formulated based on the concept of “base forces” which was introduced by Gao [157]. This method replaces the stress functions in the stress based FEMs with base forces which are easier to obtain (obtained directly from strain energy).Conventional FEM exhibits major drawback for nonlinear analysis. The conventional FEM is not able to approximate the strain and force fields accurately since these terms are dependent on the interpolation of the displacement field (displacement shape functions). This problem is avoided in BFEM which is directly based on interpolation of the internal force fields (force shape functions) [158]. Performance of BFEM is also found to be superior than conventional FEM when analyzing large strain contact problems and the nonlinear problems. This is because the deformation field in conventional FEM is complex and discontinuous near nonlinear regions and, in the case of large strain problems, the steep gradients of deformation cannot be represented accurately by the conventional shape functions [159]. BFEM has been later extended to polygonal elements based on the complementary energy principle [160] and potential energy principle [161]. Figure 14 shows an example of a BFEM element. I,J,K,L,M,N represent the edges of the element and TI,TJ,TK,TL,TM,TN represent the force vectors which act on the element edges.Figure 14
An example of a BFEM element.The stresses corresponding to a pointI on the element can be represented as [160](37)σ=1A∑I=14TI1PI1∑I=14TI1PI2∑I=14TI2PI1∑I=14TI2PI2,where TI1,TI2 are the components of force vectors which act on the center of edgeI, PI1,PI2 are the components of position vector of pointI,andA represents area of the element.


## 2.12. Other Recent Techniques/Schemes

Recently, more new schemes have been proposed for the development of polygonal/polyhedral finite elements. They are the Compatible Discrete Operator (CDO) scheme [162, 163], Hybrid High-Order (HHO) scheme [164–166], Weak Galerkin (WG) scheme [167–171], gradient correction scheme [172], and vertex-based schemes [173].Other recent techniques/methods include analysis of polygonal carbon nanotubes reinforced composite plates by using the first-order shear deformation theory (FSDT) and the element-free IMLS-Ritz method [174], an adaptive polygonal finite element method using the techniques of cut-cell and quadtree refinement [175], new adaptive mesh generation for polygonal element [176], and ultraweak formulations for high-order polygonal finite element methods [177]. New technique for 3-dimensional polyhedral elements can be seen within the framework of the finite volume method [178]. Another new approach to form polyhedral elements is by cutting a regular hexahedral element with CAD surfaces [179].


## 3. Comparison of the Various Methods

Various techniques described in Section2 above are compared with the recently proposed polyhedral element, known as Virtual Node Polyhedral Element, VPHE (three-dimensional version of Virtual Node Method, as mentioned in Section 2.7 above). The comparison is given in Table 1.Table 1
Comparison of the existing methods with the proposed/present element.
Method Element Formulation Advantages Disadvantages Specialty of VPHE element compared to other techniques Application fields and typology VCFEM Principle of minimum complementary energy. 1. Computationally efficient compared to the conventional FEM. 1. Perform poorly when the heterogeneity is in the form of voids.2. Poorly defined stress functions within the interior of the element. Stress functions within the element are well defined by monomials. Simulation of microstructures (grains) and multiscale modelling. Applicable for both 2D and 3D problems. Applicable for both linear and nonlinear analyses.  NCM-VCFEM Hybrid of VCFEM with other method such as numerical conformal mapping. 1. The variational principle is generalized.2. Solution accuracy is competitive with that of conventional FEM package in ANSYS.3. Reduces the computational cost when compared to FEM package in ANSYS. 1. The NCM-based stress function construction is expensive in comparison with conformal mapping of regular shapes such as ellipses and circles.2. NCM-based stress functions introduce singularity. Special technique (such as divergence theorem) is needed to reduce the order of these singularities. Special techniques are not needed to handle singularities in the shape functions. Real micrographs of heterogeneous materials with irregular shapes can be analyzed effectively. Applicable for 2D problems for time being. Applicable for both linear and nonlinear analyses.  HPE Hybrid stress element method together with Muskhelishvili’s complex analysis. 1. Stress functions within the interior of the element are defined by self-equilibrating stress field.2. Better performance compared to conventional FEM for plane linear elasticity problems. 1. Can contain only one irregular phase (void/inclusion) within the element. - Simulation of microstructures (grains) with heterogeneity. Applicable for both 2D and 3D problems. Applicable for both linear and nonlinear analyses.  PFEM Barycentric Coordinates. 1. Able to take arbitrary form, with arbitrary number of sides and nodes. 1. Evaluation of barycentric coordinates (complex rational functions) is neither simple nor efficient compared to the conventional displacement based FEM.2. Not efficient for assembling the stiffness matrices associated with weak solutions of Poisson equations. The shape functions consist of simple monomials irrespective on number of planes/sides. Solid mechanics and heat transfer phenomena. Computer graphics, animation and geometric modelling. Quadtree/Octree mesh generation. Applicable for both 2D and 3D problems. Applicable for both linear and nonlinear analyses.  nSFEM Coupling of conventional FEM with meshless method. 1.nSFEM is advantageous over the conventional FEM since it produces more accurate solutions, able to tolerate volumetric locking and do not require isoparametric mapping. 1.nCS-FEM increases the computational cost for solid mechanics.2. Computational time of nNS-FEM and nES-FEM is longer compared to conventional FEM for the same number of global nodes, due to larger bandwidth of stiffness matrices.3. Disadvantage of nES-FEM is that there is tendency to overestimate or underestimate the strain energy of the model for some cases. Formulation of the proposed/present element is similar tonSFEM. The current technique can be improved by carrying out smoothing technique, which will be then similar to nCS-FEM. Solid mechanics and heat transfer phenomena. Fluid-solid interaction (FSI) problems. Applicable for both 2D and 3D problems. Applicable for both linear and nonlinear analyses.  PSBFEM Semi analytical method which combines boundary element method (BEM) and FEM. 1. Analytical solutions are achieved inside the domain, discretization of free and fixed boundaries and interfaces between different materials are not required, and the calculation of stress concentrations and intensity factors based on their definition is straightforward.2. Yields highly accurate solutions for problems involving stress singularities.3. Superior to other techniques such as nSFEM and conforming PFEM within the context of linear elasticity and the linear elastic fracture mechanics. 1. Not directly applicable for unbounded domains with strongly inclined interfaces.2. PSBFEM cannot be directly used to process transient excitation as opposed to BEM.3. Not as efficient as conventional FEM or BEM when solving problems involving smooth stress variations within bounded/enclosed domain. - Solid mechanics and polygonal mesh creation. Applicable for both 2D and 3D problems. Applicable for both linear and nonlinear analyses.  MFD andVEM Surface representation of discrete unknowns (MFD) and unknown degrees of freedom are attached to trial functions within interior of the polygonal domain (VEM). 1. Efficient in solving problems involving polygonal meshes.2. Difficulties faced in integration of complex functions resulting from barycentric coordinates in PFEM are entirely avoided.3. Does not require extension of compatible interpolation functions to the interior of the element. 1. Quite difficult to present MFD due to nonexistence of trial functions for the interior of the element.2. Involve complex procedures and therefore require high computational effort. Easier to execute due to simpler element formulation. Electromagnetic field problems, convection-diffusion problems, fluid flows problems, hydrodynamics problems, eigenvalue problems, solid mechanics, heat transfer, and topology optimization. Applicable for both 2D and 3D problems. Applicable for both linear and nonlinear analyses.  VNM The polygonal domain is divided into several triangular regions which use the conventional FEM shape functions. These triangular regions are then coupled together by using mean least square shape functions. 1. Do not require formulation of complex stress functions for the element (which could be difficult for some cases, as reported for stress based FEMs such as HPE and HS-F2. Numerical integration for the elements is simple and exact, as opposed to compatible PFEM.3. Simpler and easier for computer applications compared to MFD and VEM. 1. Integration within each tetrahedron can be simplified by mapping, but the mapping procedure imposes restriction to element geometry due to high aspect ratio (limited tolerance towards mesh distortion).2. Prone to element locking. - Adaptive computation, solid mechanics, and heat transfer phenomena. Applicable for both 2D and 3D problems. Applicable for both linear and nonlinear analyses.  DGFEM Problem domain is discretized into several polygonal cells which represent the polygonal elements. The interpolation for the elements is carried out based on a set of monomial functions which are totally independent of the element. 1. Does not require any conforming shape function and the method is simple.2. Adjacent elements can be of different order and the elements do not need to be conforming. 1. Interpolation functions do not comply with shape function requirements of conventional FEM (Not compatible).2. Integrations are carried out on the boundaries of each cell. This step is an addition compared to the conventional FEM. The element fulfills all the requirements of traditional FEM. Solid mechanics, heat transfer, and eigenvalue problems on polygonal meshes. Applicable for both 2D and 3D problems. Applicable for both linear and nonlinear analyses.  T-FEM or HT-FEM Utilizes two different set of functions to approximate the solutions, one for the boundary and the other is for the interior domain. 1. Able to produce more accurate results and higher convergence rates compare to the conforming PFEM with Laplace/Wachspress shape functions.2. Elements with embedded cracks or voids can be constructed.3. Able to handle geometry induced singularities and stress/force concentrations efficiently without mesh adjustment.4. General polygonal elements with curved sides can be generated and the elements are tolerant to mesh distortions. 1. T-complete sets for some problems are either complex or difficult to formulate. - Solid mechanics and heat transfer phenomena. Applicable for 2D problems for the time being. Applicable for both linear and nonlinear analyses.  BEM-based FEM Trefftz-likebasis functions are defined implicitly and treated locally by means of BoundaryElement Methods (BEMs). 1. Applicable to general polygonal meshes (immune to severe mesh distortion).2. Computational effort is reduced, since only one subspace is needed to approximate the pricewise harmonic functions. Large linear systems of equations are generated - Adaptive mesh generation, time dependent problems, and boundary value problems. Applicable for both 2D and 3D problems. Applicable for both linear and nonlinear analyses.  HS-F Combination of principle of minimum complementary energy (similar principal used in VCFEM) with Airy stress function. 1. Possess excellent performance compared to the conventional elements and especially independent of the element geometry. 1. Not well desired for most of the engineering applications, due to the difficulties in obtaining suitable/compatible stress functions.2. Difficult to acquire nodal displacements in stress based FEMs. The element fulfills all the requirements of traditional FEM. Solid mechanics phenomena. Applicable for 2D problems for the time being. Applicable for both linear and nonlinear analyses.  BFEM Replaces the stress functions in the stress based FEMs with base forces which are easier to obtain. 1. Able to approximate the strain and force fields accurately.2. Found to be superior than conventional FEM when analyzing large strain contact problems and the nonlinear problems. 1. The Lagrange multiplier method is used to deal with the equilibrium equation. So, the stiffness matrix is a full matrix. - Solid mechanics, bonding damage detection and damage mechanics. Applicable for 2D problems for the time being. Applicable for both linear and nonlinear analyses.


## 4. Software Packages

Each method described above has been tested and analyzed by using computer programs such as MATLAB/Abaqus. These computer programs have been developed specifically for the purpose of testing and analyzing the proposed techniques. However, some of the methods have been well developed and made available as commercial software. This section described some of the software packages (either commercially available or for intended use only) which have been developed for polygonal/polyhedral FEM.VCFEM has been incorporated into a software package known as Palmyra [180]. This software can be used to design composite materials and also to determine physical properties of heterogeneous materials. Three-dimensional Voronoi cell software library (an open source software) is available in the form of MATLAB code [181]. PFEM techniques have been incorporated into computer codes by using Fortran and MATLAB as well as Java [182]. Abaqus package is available fornSFEM [183]. VEM has been developed and tested in MATLAB and Abaqus packages [184, 185]. MATLAB code on PSBFEM is used in [186].It is seen that currently there are few commercial software packages which are available for polygonal/polyhedral finite elements. However, software packages for other methods can be easily developed by incorporating the source codes developed by the researchers mentioned above with the available commercial polygonal mesh generators. Some of the software packages for polygonal/polyhedral mesh generation are Platypus (MATLAB based code) [187], ReALE [188], PolyMesher [189], PolyTop [190], OpenMesh [191], and more, which can be found in [192].


## 5. Summary

It can be seen that various finite elements have been proposed for engineering analysis. These elements have been proposed to facilitate meshing of the problem domain, to facilitate the analysis of physical phenomena, and to overcome drawbacks or limitations in the existing methods. This review enables the readers to identify advantages, disadvantages, and a comparison between the various techniques used in formation of polygonal/polyhedral finite elements.



---

*Source: 5792372-2018-10-04.xml*