A Brief Review on Polygonal/Polyhedral Finite Element Methods

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


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/schemes The methods above are briefly highlighted in the following sections.

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.The matrix phase is denoted as Ω  and heterogeneity phase as Ω ℎ .The heterogeneity surface and element surface are denoted as Ω ℎ and Ω  , 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] where Ω is the computational region of the element Ω represents boundary of the element  is a 6-by-1 matrix containing equivalent stress components  is a 6-by-6 inverse of the material property matrix (elastic compliance matrix) is a matrix containing plastic flow parameters  is a matrix of partial derivative of flow potential function with respect to the stress  and  are the matrices for displacement and traction components along the element boundary, respectively  represents boundary surface Stress 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]: where   ( = 1, 2, . . .) 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].
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].

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.
Figure 2 shows a HPE in global - coordinate system and its equivalent mapping to a standard ellipse in reference - coordinate system.
An example of a hybrid functional ∏ ()  for a HPE with hole (void) is given as [25] () where  and  are the matrices for displacement and traction components along the element  1 and  2 represent the element's outer boundary with adjacent elements and the inner interface between the matrix and heterogeneity, respectively ũ(1) represents the components of the specified displacements along  1  represents boundary surface The interior stresses are obtained by using interpolation functions (in the form of trigonometric functions) through the following relation [25]: where  1 to  3 and  1 and  3 are the functions of , , and .  and   are the upper and low limits of the series. and  are the parameters that resulted due to the mapping of the computational/matrix region Ω to the elliptical region in reference system.Expression for  1 is shown here as an example: where  1 and  2 are functions in terms of , , , and . and  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.

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][37][38][39][40], mean value coordinates [41,42], harmonic coordinates [43], maximum entropy coordinates [44][45][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][49][50][51][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  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 neighborbased 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 0  based on Wachspress [47] for a polygonal element in the - coordinate system according to Figure 3 is given as where  represents area enclosed by the three nodes within the bracket   ,   , and   represent a particular external/surface node of the polygonal element  represents inner node of the polygonal element  The expression for   (, ) in ( 7) can be rewritten by using the angles formed ( and Ψ) between the nodes [47] as The expression for   (, ) in ( 7) based on mean value coordinates [47] for a polygonal element is given as where ∝ is the angle that is formed within the triangular partitions at the inner node, as shown in Figure 4.
The expression for   (, ) in ( 7) based on the concept of natural neighbors (Laplace shape functions) [47] for a polygonal element according to Figure 5 is given as where   represents length of the Voronoi edge and

𝑛-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 as n-sided polygonal smoothed finite elements (nSFEM) and implemented the method in solid mechanics (macrolevel).
In nSFEM, 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 of nSFEM element is shown in Figure 6.Point 0 in Figure 6 represents the center of the polygonal element.The displacement at this point,  0 is taken as the average of displacement of all the external nodes, given by the following equation [71]: where   represents displacement at a particular node  and  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]: where   ( = 1, 2, 3) represents the conventional shape functions for a 3 nodes' triangular element.Substituting ( 12) into ( 11) and simplifying gives The shape function matrix for subtriangle 0-1-6 then becomes The strain components are smoothed according to [71] where    represents area of smoothing domain, B represents the conventional compatible strain displacement matrix, and Ω   represents the smoothing domain.
In case of nCS-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 of nES-FEM, the smoothing is done based on particular edge/side (instead of cell as in nCS-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 of nNS-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].
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 of nCS-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 of nNS-FEM is that it is immune from volumetric locking phenomena.Disadvantage of nNS-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 of nES-FEM is that there is a tendency to overestimate or underestimate the strain energy of the model for some cases.Apart from that, similarly to nNS-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 of nSFEM is provided by Nguyen et al. [72], for solid mechanics problem.It is shown that nES-FEM provides most accurate solution compared to the others and the stiffness/softness of the model of nES-FEM is in between the other two techniques.Combination of nES-FEM and nNS-FEM (termed as nES/NS-FEM) to avoid volumetric locking and to achieve faster convergence can be seen in the literature [72,75].Applications of nSFEM 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 fluidsolid 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 the nSFEM in FSI is that it is capable of performing independent domain discretization.
Generally, the nSFEM 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][79][80].

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.Solution along the radial direction is obtained by analytical expressions by using m number of shape functions, where m 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 ,  and Cartesian coordinate system x, 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.
The displacement (, ) within a particular triangular subdomain of the element is interpolated through utilization of similar shape functions as shown in (18) below [82]: where () = {  1 ()  2 () } 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][85][86][87][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][97][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 as nSFEM and conforming PFEM within the context of linear elasticity and the linear elastic fracture mechanics [83].

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]: where K represents the full conductivity tensor, u represents the temperature, and q 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: where  →  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 ∇.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.
An additional step is necessary, which is the discretization of the integrals.This step is required in order to approximate div and ∇ (differential operators) according to mimetic approach.
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]  , = {V ∈  1 () : V | ∈ B  () , V | ∈ P −2 ()} , (21) where  , 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,  1 () represents the space containing K,  represents a generic edge on the polygon, B  () represents the set of polynomials of degree less than or equal to k on , and P −2 () denotes the space of polynomials of degree less than or equal to k-2 on K.For k = 1, the trial functions are linear on each edge and the inside of the element is represented by harmonic functions.For k = 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].

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.
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]: where  represents unknown field variable at a particular point/node, φ represents the field variable at the virtual node,   ,   , and   are the conventional FEM shape functions of a 3 nodes' triangular element, and  V represents a particular triangular element within the polygonal domain.For triangle  22) above can be written as Field variable at the internal node (virtual node) φ7 can be represented by the least mean square [123]: where  represents total number of nodes of the polygonal element,   =   (, ) −1 ,  is a set of basis functions, and  and  are basis matrices.Equation ( 24) can be rewritten in terms of shape function for a particular node: 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: where l represents total number of elements in the mesh.It can be seen that this method is quite similar to nCS-FEM Mathematical Problems in Engineering 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].

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][130][131][132][133][134][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] where represents matrix of the local degrees of freedom, Ne represents the number of polygonal elements, {} = [ 1  1  2  2 . . .    ] represents matrix containing the incompatible interpolation functions b, and  represents the element support [130]: The 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]: where  and  are the local coordinates.The interpolation functions can also be made of radial basis functions [130]: where  2  = ( −   ) 2 + ( −   ) 2 and (  ,   ) represents coordinates of the polygon vertices as well as the mass center.
The stiffness matrix  for a heat transfer phenomenon is obtained through the formula [130] where ⟦⟧ represents discontinuity jump, ⟨⟩ represents mean value of , Ω represents the domain space, B represents matrix containing the differentiation of interpolation functions, k represents thermal conductivity,   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.

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.
The displacement  within a domain (interior) can be interpolated by using where ǔ represents the known function on the boundary,   is the trial function,   is the unknown coefficient, and m represents number of trial functions.The trial functions   can be generated from Muskhelishvili's complex variable formulation [139].The displacement ũ on the boundary of the domain can be interpolated by using (33) below: where  represents vector of nodal displacements and Ñ() = { (1−)/2 (1+)/2 represents matrix of the corresponding onedimensional 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][146][147][148]. This method is developed based on the Trefftz function.

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][150][151][152][153][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.
The displacement along a particular element edge d is given by [155] where  = {      } represents the vector of shape functions for the three nodes (i, j, k) along a particular edge of the element and  = {  V    V    V  } represents the displacement vector for the nodes.The shape functions are given as [155] where  represents local coordinate system along each element boundary.The element stress fields are given as [155] where S is the stress solution matrix (with dimensions 3 by k) which is derived from k number of analytical solutions of the Airy stress functions,  represents matrix (with dimensions k by 1) of unknown stress parameters and  * represents particular solution corresponding to body forces.

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., , , , ,  represent the edges of the element and   ,   ,   ,   ,   ,   represent the force vectors which act on the element edges.
The stresses corresponding to a point I on the element can be represented as [160] where  1 ,  2 are the components of force vectors which act on the center of edge I,  1 ,  2 are the components of position vector of point I, and A represents area of the element.
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].

Comparison of the Various Methods
Various techniques described in Section 2 above are compared with the recently proposed polyhedral element, known as Virtual Node Polyhedral Element, VPHE (threedimensional version of Virtual Node Method, as mentioned in Section 2.7 above).The comparison is given in Table 1.

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 for nSFEM [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].

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.

Figure 2 :
Figure 2: Mapping of a Voronoi cell element.

Figure 3 :
Figure 3: Construction of shape functions based on Wachspress.

Figure 4 :
Figure 4: Construction of shape functions based on mean value coordinates.

Figure 5 :Figure 6 :
Figure 5: Construction of shape functions based on the concept of natural neighbors.

Figure 10 :
Figure 10: Low order degrees of freedom for MFD method.The arrows represent fluxes and the center represents the temperature.

2 Figure 11 :
Figure 11: Partitioning of a polygon according to VNM.

Figure 12 :
Figure 12: An example of T-FEM with its shape functions utilized for the case of solid mechanics.

Figure 13 :
Figure 13: An example of a HS-F element.

Figure 14 :
Figure 14: An example of a BFEM element.

Table 1 :
Comparison of the existing methods with the proposed/present element.