PSBFEM-Abaqus: Development of User Element Subroutine (UEL) for Polygonal Scaled Boundary Finite Element Method in Abaqus

*e polygonal scaled boundary finite element method (PSBFEM) is a novel method integrating the standard scaled boundary finite element method (SBFEM) and the polygonal mesh technique. *is work discusses developing a PSBFEM framework within the commercial finite element software Abaqus. *e PSBFEM is implemented by the User Element Subroutine (UEL) feature of the software. *e details on the main procedures to interact with Abaqus, defining the UEL element, and solving the stiffness matrix by the eigenvalue decomposition are present. Moreover, we also develop the preprocessing module and the postprocessing module using the Python script to generate meshes automatically and visualize results. Several benchmark problems from two-dimensional linear elastostatics are solved to validate the proposed implementation. *e results show that PSBFEM-UEL has significantly better than FEM convergence and accuracy rate with mesh refinement. *e implementation of PSBFEM-UEL can conveniently use arbitrary polygon elements by the polygon/ quadtree discretizations in the Abaqus. *e developed UEL and the associated input files can be downloaded from https:// github.com/hhupde/PSBFEM-Abaqus.


Introduction
e finite element method (FEM) is a reliable computational tool to solve partial differential equations (PDE) in science and engineering [1][2][3][4]. A domain of complex geometry is partitioned into a finite number of nonoverlapping subdomains of simplex shapes by introducing the concept of discretization [5]. Typically, the shape of the conventional two-dimensional finite element method is triangles or quadrilaterals. At present, the conventional FEM also faces several problems. For example, (a) the accuracy of the solution depends on the quality of the mesh and (b) requires sophisticated discretization techniques to generate highquality meshes and to capture topological changes [6]. ese elements used in conventional FEM must conform to the domain's boundary, which leads to difficulty in solving many complex problems.
To overcome the weakness above the conventional FEM, researchers proposed other alternatives methods, such as the meshfree method [7][8][9], the smoothed finite element method [6,10], the Isogeometric Analysis (IGA) [11], Deep Neural Networks (DNNs) [12], and the polygonal finite element method [10,13]. Polygonal element with more than four sides involves more nodes in their interpolation compared with a conventional FEM. Simultaneously, they generally exhibit superior solution accuracy [5]. Moreover, it is more flexible in the discretization of complex geometry. erefore, these advantages have further motivated polygonal elements as an alternative to conventional FEM using triangles or quadrilaterals. e scaled boundary finite element method (SBFEM) is an alternative method to construct polygonal elements. Song and Wolf developed the technique in the 1990s [14]. e scaled boundary finite element method is a semianalytical method that attempts to fuse the advantages and characteristics of FEM and the boundary element method (BEM) into one new approach. e SBFEM has been applied to many physical field problems, such as wave propagation [15,16], heat conduction [17,18], fracture [19,20], acoustic [21], seepage [22,23], elastoplastic [5], and fluid [24]. For these problems, the SBFEM presents more efficiency compared with the conventional FEM. e polygonal scaled boundary finite element method (PSBFEM) is a novel method integrating the standard SBFEM and the polygonal mesh technique. is method is flexible in meshing complex geometries, and the use of polygons to discretize the computational domain naturally complements the SBFEM. Recently, an alternative mesh technique has been widely used in geometric discretization. In computational mechanics, the quadtree algorithm is usually used in large-scale simulations typical in the modeling of earthquake and ground motions [25], flood [26], and tsunami [27]. e quadtree algorithm is fast, efficient, and capable of achieving rapid and smooth transitions of element sizes between mesh refinement regions [28]. Mesh generation and adaptive refinement of quadtree meshes are straightforward [29]. Due to hanging nodes between two adjacent elements of different sizes, it is problematic that quadtree meshes are directly used to simulate within the finite element method's framework. However, the SBFEM only discretizes in the boundary of geometry. Hence, each cell in a quadtree mesh, regardless of hanging nodes, is treated as a generic polygon. is enables the structure of the quadtree to be exploited for efficient computations. e ability to assume any number of sides also enables the SBFEM to discretize curved boundaries better.
Although the PSBFEM has become quite mature, these studies only exist as a few stand-alone codes. e PSBFEM has not yet formed a part of commercial software at present. Hence, it is not easy to use the method in a specific community or laboratory. e commercial software Abaqus has powerful linear or nonlinear, static, or dynamic analysis capabilities. Also, Abaqus/Standard analysis provides a User Element Subroutine (UEL) to define an element with a very available option to interface with the code. Liang et al. [30] developed a UEL for dynamic analysis of saturated porous media based FEM. Kumbhar et al. [6] implemented the element based smoothed finite element method (CSFEM) by the UEL subroutines. Molnar et al. [31] implemented an implicit, staggered elastoplastic version of the phase-field approach Abaqus through the UEL. erefore, the UEL can be used to extend the Abaqus to solve new problems conveniently.
Recently, Ya et al. [32] implemented an open-source polyhedral SBFEM element for three-dimensional and nonlinear problems in the commercial software Abaqus through the UEL. e code can use the polyhedral element to enhance the performances of ABAQUS for interfacial problems, and it significantly reduces the meshing burden encountered in the FEM. At the same time, it also implemented the technique of octree mesh generation, which is efficient and robust to mesh complex geometries and promise to integrate geometric models and numerical analysis in a fully automatic manner.
In engineering designs, to calculate models to be used as design tools, two-dimensional (2D) models are relatively easy to set up and have reasonably short computational times, which would allow sensitivity and optimization analyses [33]. ere are no available subroutines of SBFEM to solve 2D problems in the commercial software Abaqus at present. However, the implementation of two-dimension can be easily derived from three-dimensional implementation. Hence, Ya et al. [32] provided a reference for developing the UEL to solve 2D problems using the polygonal scaled boundary finite element method with the polygon/quadtree meshes.
In addition, to solve the stiffness matrix, we need to perform the eigenvalue decomposition in the SBFEM. Ya et al. [32] used the Linear Algebra Package (LAPACK) [34] to perform the eigenvalue decomposition by adding the LAPACK source codes into UEL subroutines. Hence, LAPACK needs to be compiled every time. It will take extra time when calling UEL subroutines. Abaqus provides mathematical libraries (the Intel Math Kernel Library, MKL [35]) to perform the eigenvalue decomposition. e user can directly use MKL only by modifying the Abaqus environment file, which can avoid taking extra time when calling UEL subroutines.
In this paper, we implement the PSBFEM technique with Abaqus and provide the details. For simplicity, we consider a two-dimensional linear elastostatics problem using the PSBFEM technique. is work is organized as follows: the basic principles of the SBFEM are introduced in Section 2. Section 3 describes the implementation of PSBFEM by the Abaqus UEL subroutine. Moreover, we also develop the preprocessing and postprocessing module using the Python script. Section 4 presents a detailed convergence study by comparing the convergence and accuracy rates with the conventional FEM using several numerical benchmarks in two-dimensional linear elastostatics. Finally, the concluding remarks of this work are summarized in Section 5.

e SBFEM Equations in the 2D Linear Elastostatics.
e fundamental difference between SBFEM and FEM is introducing a scaling center O [36]. As shown in Figure 1, a scaling center can be directly visible from any point on the whole boundary of the S-element. By a scaling center, we can transform the scaled boundary coordinates to polar coordinates. Hence, the S-element can be described by a radial coordinate ξ and a circumferential coordinate η.
In this section, we consider a two-dimensional isotropic linear elastostatics problem. e governing equation and the boundary conditions can be expressed as where = is the differential operator, σ is the Cauchy stress tensor, b is the body force, Ω is the computational domain. Γ u is the displacement boundary conditions, and Γ t is the surface traction boundary conditions. e displacement u and the surface traction t are imposed on the domain boundaries Γ u and Γ t with the outward unit normal n.
As illustrated in Figure 2, the SBFEM presents a local coordinate system (ξ, η). e coordinates of a point (x, y) along the radial line and inside the domain can be expressed as follows [36]: where ξ,η are the scaled boundary coordinates in two dimensions, ξ is a radial coordinate, and η is the circumferential coordinate. e differential operator can be transformed from the Cartesian coordinate system x, y into the scaled boundary coordinates system ξ, η as follows: where the Jacobian matrix at the boundary (ξ � 1) can be expressed as A comma followed by a subscript is used to denote partial differentiation to the variable in the subscript. e displacement field u(ξ, η) at any point in SBFEM coordinates is written as where u(ξ) is radial displacement functions along a line connecting the scaling center O and a node at the boundary.
[N u (η)] is the shape function matrix Bounded domains Unbounded domains e strain field ε { } is expressed in the scaled boundary coordinates as e stress field can be express as where [D] is an elasticity matrix and has the form and for plane stress cases, E is Young's modulus, and v is Poisson's ratio. According to the virtual work principle, the radial displacements function u(ξ) is the solution of the SBFEM equation in displacement [36]: where u(ξ) { } ,ξξ is the second partial derivative with respect to the variable ξ, u(ξ) { } ,ξ is the first partial derivative with respect to variable ξ, [E i ], i � 0, 1, 2, are the coefficient matrices, and F(ξ) is a load vector.
In the derivation, the governing equations of linear elasticity are weakened in the circumferential direction, while the strong form remains in the radial direction. e coefficient matrices [E i ], i � 0, 1, 2, depends only on the geometry and material properties. ese coefficient matrices can be given as follows: When F(ξ) � 0, the solution of equation (12) is secondorder ordinary differential equations that can be obtained by introducing the variable: where q(ξ) is the internal force vector. Equation (12) can be transformed into a first-order ordinary differential equation with twice the number of an unknown as where the coefficient matrix e solution of the equation (15) can be obtained by the mathematical theorem and eigenvalue decomposition technique. e eigenvalue decomposition of [Z p ] is expressed as u ] and [Φ (p) u ] are the transformation matrices corresponding to the modal displacements and forces, respectively. e general solution of equation (15) can be written as where c (n) and c (p) are the integration constants. For a super element with 0 ≤ ξ ≤ 1, the solution of equation (15) is where the integration constants c (n) can be extracted from the nodal displacements on the boundary Eliminating the integration constants c (n) at ξ � 1 is given Hence, the stiffness matrix of the S-element can be written as e displacement field u(ξ, η) inside subdomains by a line element on the S-element can be given

e Element of PSBFEM.
e PSBFEM is inherently appropriate for modeling polygons and has other promising capabilities. Because the PSBFEM is discretized only in the boundary, and an element of PSBFEM can assume more complex shapes than a finite element method, the PSBFEM element is much more flexible than the FEM element. With such an important advantage, the PSBFEM is more applicable to complex geometries than other methods. Figure 3 shows the supported element type in the PSBFEM-UEL. In the PSBFEM-UEL, we can use the Abaqus element type, such as CPS3, CPE3, CPS4, and CPE4. Besides, the PSBFEM-UEL also provides the complex element type: polygonal elements and complex quadrilateral element (quadtree discretization). Hence, it is an effective tool to solve complex elements in numerical analysis.

Automatic Meshes Generation.
is section uses a Python script to automatically generate PSBFEM elements by the Delaunay triangulation [37]. is algorithm is robust and efficient, and it has been integrated into Abaqus CAE. e algorithm can be used to triangulate any set of points on a two-dimensional plane. e triangulated mesh is then used to generate the polygon elements. e detail of the polygon generation algorithm mainly contains two steps. Firstly, we automatically generate triangular mesh by the Abaqus CAE. Secondly, considering each triangle interior node as the center of the polygon element, a polygonal element can be generated by connecting the centroids of all the triangular elements circumventing. More detail of this algorithm is presented in Figure 4. Finally, the PSBFEM-UEL can be directly used with a polygon mesh generator to analyze complex geometry problems. e quadtree decomposition is a tree data structure in which each parent has precisely four children [38]. e quadtree meshes are fast, efficient, and capable of achieving rapid and smooth transitions of element sizes between mesh refinement regions. e quadtree mesh can provide vital support in the preprocessing for the adaptive analysis of the SBFEM. We developed a quadtree mesh automatic generation code by a Python script. A simple example of quadtree discretization with three levels is illustrated in Figure 5. e quadtree algorithm is described in more detail in [28,38].

UEL Implementation of SBFEM in Abaqus
e Abaqus/Standard analysis provides a programming interface UEL to define the customized elements. In this section, we present the major implementation details of UEL for the PSBFEM in Abaqus. e Abaqus solver can be carried out with the UEL subroutine by the command: abaqus job � 〈input file name〉user � 〈UEL subroutine file〉.
(25) Figure 6 shows the framework sketch of implementing the PSBFEM within Abaqus.
is system contains three parts: preprocessing module, PSBFEM-UEL module, and postprocessing module. A Python script is used in preprocessing to generate polygon SBFEM mesh and define the loading and boundary conditions. Finally, we can obtain an Abaqus input file by preprocessing. e Abaqus CAE does not support the visualization of the UEL element. us, a Python script is provided in the postprocessing to extract the VTU format results and visualize them in the software Paraview [39].

e Implementation of PSBFEM-UEL.
e most critical work of UEL is to update the contribution of the element to the internal force vector RHS and the stiffness matrix AMATRX according to the information from ABAQUS/ Standard analysis. In this paper, the UEL code for implementing the PSBFEM is written in FORTRAN77. Figure 7 shows an overall implementation of the UEL subroutine. Based on the input file's connectivity information, the UEL computes the scaling centers and transforms the global coordinate into the local coordinate. Equation (13)  To solve the stiffness matrix, we need to employ the eigenvalue decomposition (see equation (17)). At present, many mathematical libraries to perform the eigenvalue decomposition exist. In this work, we use the Intel Math Kernel Library (MKL) [35] to decompose the eigenvalue. In the Abaqus/Standard analysis, we can directly use MKL by modifying the Abaqus environment file.

Defining the UEL Elements.
Abaqus's input file usually contains a model information section (such as defining nodes, elements, the active degrees of freedom, and materials). At present, this information cannot be set in Abaqus CAE and must be defined through an input file. e Abaqus provides the keyword * USER ELEMENT to defined a new user element. e main contents of defining the user element are as follows: (1) Assigning an element type key to a user-defined element and the number of nodes. e element type key must be of the form "Un" in Abaqus/Standard analysis, where ′ n ′ is a positive integer that identifies the element type uniquely. In this implementation of PSBFEM, the integer is equal to the number of nodes of the element.  In the case of PSBFEM, we present a simple polygonal mesh of PSBFEM (see Figure 8).
is mesh consists of three element types: triangular element (U3), quadrilateral element (U4), and Pentagon element (U5). In the input file (see listing 1), 1 ∼ 19 is the line number; the actual input file does not contain the line number. Lines 1 ∼ 7 are used to define two triangular elements (U3). Line 1 assigns the element type, the number of nodes, the number of element properties, and the number of freedom degrees per node. Line 2 sets the active degrees of freedom. Lines 3 ∼ 5 define the element sets' E3'. Lines 6 ∼ 7 set the element properties (Young's modulus and Poisson's ratio) of "E3". Similarly, other types of elements are defined using the same way.

Numerical Examples
In this section, we validate the convergence and accuracy of the implementation by solving a few benchmark problems. Moreover, the results of the PSBFEM are compared with the FEM. e FEM analysis uses the commercial finite element software Abaqus. For validation, the relative error L 2 norms in the displacement are computed as follows: where u h is the numerical solution and u is the analytical or reference solution. e relative error in the energy as given by (27) where ε h is the numerical solution and ε is the analytical or reference solution.

A Two-Dimensional Cantilever Beam.
A two-dimensional cantilever beam of height H � 1.0 m and length L � 5.0 m subjected to a uniform loading P � 10 kPa is considered, as shown in Figure 9(a). e material properties are given by Young's modulus E � 2 GPa and Poisson's ratio v � 0. e right boundaries are constrained without displacement (ΔX � 0, ΔY � 0). e domain is discretized with the quadrangle and arbitrary polygonal elements. A convergence study is performed by mesh refinement. e meshes are refined successively following the sequence n � 2, 4, 8, 16, 32. e element size is chosen as (h � H/n). A representative mesh is presented in Figures 9(b) and 9(c). In this work, the Abaqus uses the CPS4 element. e analytical solution of the stress fields can be expressed as [40][41][42] (28c) e results of the relative error in the vertical displacement u y for the point O are given in Table 1. It can be observed that the errors decrease as the mesh refinement. e errors of PSBFEM are less than the FEM at the same element size. Moreover, the polygonal element shows higher accuracy than the quadrilateral element in the PSBFEM because the polygonal element has more nodes for the same element size. Figure 10 shows that the results of mesh size sensitivity for three element types. It is clear that these elements would obtain more accurate results when mesh sizeis less than 1/8 m. e convergence of the relative error in the displacement and the energy norm with mesh refinement are presented in Figures 11 and 12. It is observed that the PSBFEM converges to an exact solution with an optimal convergence rate.
Moreover, Figure 13 shows the contours of the vertical displacement u y of Abaqus CPS4, PSBFEM-Quad, and PSBFEM polygon. It is clearly shown that the results are virtually the same for the FEM and PSBFEM. In addition, Figure 14 presents the computational cost comparison of the developed UEL and Abaqus standard elements. Due to the increment size setting affecting the solving time, we use the automatic incrementation type. e comparison is evaluated with an Intel Core i7-4710MQ CPU running at 2.50 GHz and 4.0 GB of RAM. e total CPU time is normalized. It is noted that computation costs for the PSBFEM-UEL and Abaqus standard elements are comparable and of the same trend. Moreover, the computational cost of PSBFEM-UEL is slightly higher than the Abaqus standard element. It is due to the semianalytical method of SBFEM being bound to generate relatively more computations in its concept than FEM [43].

Infinite Plate with a Circular Hole.
In this example, an infinite plate with a circular hole of radius a under remote uniaxial tension σ x is considered as shown in Figure 15(a). Considering the symmetry of geometry, we only analyze a quarter of the infinite plate, as shown in Figure 15(b). e analytical solution of stress fields in the polar coordinates (r, θ) is expressed as [36]

Mathematical Problems in Engineering
and the displacement fields can be expressed as (30b) e shear modulus G and Kolosov's constant κ can be expressed as where a is the radius of the hole. e circular hole radius is 0.4 m, the length of the quarter of the circle hole is 1.0 m, and remote tension σ x � 1 kPa. e domain is discretized using quadrilateral and polygonal meshes, respectively. e problem is also modeled by the plane stress, and the material properties are given by Young's modulus E � 1 × 10 5 Pa and Poisson's ratio v � 0.25. e vertical displacement of bottom boundaries (ΔY � 0) and the horizontal displacement of left boundaries (ΔX � 0) are constrained without displacement. e results of PSBFEM and FEM in terms of convergence behavior and accuracy are compared and demonstrated. Figures 16 and 17 show the convergence of the relative error norms decrease in displacement and the energy norm with increasing elements. It is observed that the PSBFEM element is significantly more accurate than the FEM element. Besides, PSBFEM Polygon has a more slightly fast convergence rate.
From the results, it is clear that all the methods asymptotically converge to the analytical solution with mesh refinement. It is noted that the PSBFEM requires fewer DOFs when compared to conventional FEM. Besides, Figure 18(a) provides that the contours of the vertical displacement u y of Abaqus CPS4, PSBFEM-Quad, and PSBFEM polygon in an infinite plate with a circular hole. e results show a good agreement for the FEM and PSBFEM. Moreover, Figures 18(b) and 18(c) also show that the distribution of stress component σ y and strain component ε y obtained from the PSBFEM and the FEM by Abaqus. Similarly, good agreement between the three sets of results is observed.

Pressurized
ick Cylinder Modeled by a Quarter-Annulus Model. In this example, we consider a benchmark problem of a quarter-thick cylinder subjected to internal pressure at the inner circular edge. e analytical solution of stress fields in the polar coordinates is [44] and the displacement fields can be expressed as where R a and R b are the inner and the outer radius of the cylinder, respectively, with R a � 1 m and R b � 2 m. P is the pressure exerted along the inner circular edge, with P � 1000 Pa. e material properties are given by Young's modulus E � 1 × 10 5 Pa and Poisson's ratio v � 0.30. e geometry and boundary conditions are presented in Figure 19(a). e domain is discretized with the quadrilateral and polygonal elements, as shown in Figures 19(b) and 19(c). Figures 20 and 21 show the convergence of the relative error in both the displacement and the energy norm with mesh refinement. Results show that all the techniques asymptotically converge to analytical solutions with reducing element size, and the PSBFEM Polygon shows slightly accurate results. Moreover, a good agreement between PSFEM and FEM is observed in Figure 22. For the error estimation, the relative error L 2 norms in the displacement are computed as equation (26). Moreover, a simple regression is investigated to compare the results of FEM and PSBFEM. e quadtree mesh is generated by setting the same number of mesh seeds on every hole, as shown in Figure 23(c). It is clearly present that the mesh transition between the holes of different sizes is effectively handled. e square body with multiple holes is modeled with 1632 quadtree elements, and the total nodes are 2478.
Moreover, this problem is also analyzed with a similar number of nodes (2532 nodes) using the Abaqus CPS4 element. Due to the symmetry of geometry, only the results on the right side of the top edge are provided for comparison. Results in Figure 24 show the comparison between quadtree PSBFEM and Abaqus CPS4 in the vertical displacement u y . A simple regression on the Abaqus CPS4 and quadtree PSBFEM had presented R 2 � 0.99, and the error in the L 2 norm is 0.0097. e strain energy of FEM is 85.934 kJ, and the strain energy of PSBFEM is 86.5917 kJ. e relative error of strain energy is 0.76%. Hence, these results indicate PSBFEM accuracy and reliability for the quadtree mesh. Also, the contour plots of the vertical displacement u y obtained from the PSBFEM analysis and the FEM analysis by Abaqus are shown in Figure 25. It is noted that the contour plots present a good agreement.    (ΔX � 0, ΔY � 0), and the right edge is applied a horizontal displacement U � 0.1 m, as shown in Figure 26. We chose four nodes, A, B, C, D, to compare results, as shown in Figure 26. e square body with a rabbit shape cavity is modeled with 3147 quadtree elements. Moreover, this problem is also analyzed with a similar number of elements using the Abaqus CPS4 element, and the total elements are 3358, as shown in Figure 27. Table 2 shows the relative error in the horizontal displacement for the PSBFEM and FEM. e relative errors are less than 0.5%. Moreover, the strain energy of FEM is 4.304 kJ, and the strain energy of PSBFEM is 4.299 kJ. e relative error of strain energy is 0.12%. Figure 28 presents that the contour plots present a good agreement for FEM and PSBFEM. erefore, the quadtree mesh of PSBFEM has sufficient accuracy and reliability.

Conclusions
is paper implements the PSBFEM of two-dimensional linear elastostatic problems within the Abaqus/Standard analysis by the UEL subroutine. is work mainly focuses on the main procedures to interact with Abaqus, defining the UEL element in the input file, and solving the stiffness matrix by the eigenvalue decomposition in the UEL implementation procedure. Also, we discuss the automatic mesh generation of polygon/quadtree and the visualization of results by the Paraview. e implementation of PSBFEM is validated against the FEM by solving a few benchmark problems. e results demonstrate that PSBFEM-UEL has a significantly better than FEM convergence rate. Moreover, the polygon mesh has a higher accuracy rate than the quadrangle mesh in the PSBFEM-UEL. Notably, the implementation of PSBFEM can conveniently use arbitrary polygon elements by the polygon/quadtree discretizations in the commercial finite element software Abaqus. In the future, the scope of this approach developed here can be extended to higher-order elements. e source code of the implementation can be downloaded from https://github.com/hhupde/PSBFEM-Abaqus with input files of numerical examples presented in this work.

Abbreviation
FEM: e finite element method PDE: Partial differential equations IGA: Isogeometric analysis DNNs: Deep neural networks SBFEM: e scaled boundary finite element method PSBFEM: e polygonal scaled boundary finite element method UEL: User element subroutine MKL: Intel math kernel library CPS3: 3-node linear plane stress elements CPE3: 3-node linear plane strain elements CPS4: 4-node bilinear plane stress elements CPE4: 4-node bilinear plane strain elements RAM: Random-access memory CPU: Central processing unit.

Data Availability
e data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest
e authors declare that they have no conflicts of interest.