Implementation of arbitrary polyhedral elements for automatic dynamic analyses of three-dimensional structures

The transition from the geometry to the mesh can be rather difficult, manual and time-consuming, especially for the large scale complex structures. The procedure of mesh generation needs massive human intervention making the automatic engineering analyses of structures from CAD geometry models hardly possible. This paper focuses on implementing a polyhedron element with arbitrary convex topology based on the Scaled Boundary Finite Element Method (SBFEM) in ABAQUS on the strength of the interface of UEL (the subroutine to define a general user-defined element) for automatic dynamic analyses of three-dimensional structures. This implementation empowers ABAQUS to analyze any model with arbitrary polyhedron elements. When the geometry of a structure is obtained from CAD, the dynamic analyses can be launched seamlessly and automatically. Cases of a cantilever subjected to a dynamic harmonic excitation with the traditional hexahedron element and this polyhedron element are compared to verify the accuracy of the UEL. Taking a practical example of the Soil-Structure Interaction analysis of a Nuclear Power Plant, the applicability and performance of this implementation are tested. The results of the two examples confirm that this polyhedron element based on SBFEM can be more accurate using much less degrees of freedom and its implementation in ABAQUS is robust and compatible.

The transition from the geometry to the mesh can be rather difficult, manual and time-consuming, especially for the large scale complex structures. The procedure of mesh generation needs massive human intervention making the automatic engineering analyses of structures from CAD geometry models hardly possible. This paper focuses on implementing a polyhedron element with arbitrary convex topology based on the Scaled Boundary Finite Element Method (SBFEM) in ABAQUS on the strength of the interface of UEL (the subroutine to define a general user-defined element) for automatic dynamic analyses of three-dimensional structures. This implementation empowers ABAQUS to analyze any model with arbitrary polyhedron elements. When the geometry of a structure is obtained from CAD, the dynamic analyses can be launched seamlessly and automatically. Cases of a cantilever subjected to a dynamic harmonic excitation with the traditional hexahedron element and this polyhedron element are compared to verify the accuracy of the UEL. Taking a practical example of the Soil-Structure Interaction analysis of a Nuclear Power Plant, the applicability and performance of this implementation are tested. The results of the two examples confirm that this polyhedron element based on SBFEM can be more accurate using much less degrees of freedom and its implementation in ABAQUS is robust and compatible.
As the development of the Computer Aided Design (CAD) and Computer Aided Engineering (CAE) tools in the industry, the design, safety assessment, health monitoring of some important modern infrastructures tend to need more and more refined and sophisticated geometry models and mesh models. Nevertheless, there is a gap between the CAD and the CAE holding a barrier between the geometry model for design and the mesh model for engineering analysis. Because except for some rarely used meshless methods, it is a difficult, manual and time-consuming process to generate a fine mesh from a complex geometry. During the manual mesh generation, the balance of accuracy and efficiency should be examined, especially for dynamic analysis in the time domain. A basic requirement is to assign heterogeneous distributed grid densities for the different location of the geometry model to balance the calculation accuracy and efficiency, which are case-sensitive, completely experimental and can be rather tedious. The automatic mesh generation technique with polyhedron elements can break this barrier due to the flexible topology that makes them highly adaptable to the complex shape and smooth transition of mesh density, which is the foundation of the automatic engineering analysis.
Since the arising of numerical methods for the Partial Differential Equation (PDE), the mesh generation has become a primary procedure for these methods that need spatial discretization like the Finite Element Method (FEM) and the Finite Difference Method (FDM). From a geometrical point of view, a mesh is the approximate shape representation of the domains or the media (or can be collectively referred to as configurations) of the original physical problem, which consists of many non-overlapping continuous small domains called elements. These elements have certain topology defined by geometric feature points such as vertices called nodes, where the degrees of freedom are defined on. Except that a numerical analysis cannot be carried on without a mesh, the quality, size and type of the mesh have a great influence on the reliability of the results as well.
For the last four decades, lots of mesh generation technologies have been developed from two-dimension (2D) domain to three-dimension (3D) domain on the strength of the computer graphics. During the 1950s, mesh generation techniques based on coordinate transformations firstly came up, such as the conformal mapping 1,2 technique. This technique was extended to simple 3D configurations by South et al. 3 and Baker 4 . Although the methods above can generate structured hexahedral or quadrilateral meshes with high quality, the process of these methods requires lots of time, manual guidance and very specific geometry. In the 1980s, triangulation methods emerged generating unstructured triangular meshes for 2D configurations or tetrahedral meshes for 3D configurations. After the generation of a triangulation mesh, an unstructured hexahedral mesh can be obtained by the domain partition [5][6][7][8] . Unstructured meshes can always be created by triangulation methods for any intricate geometry automatically and efficiently. But these kinds of meshes have low accuracy and additional artificial stiffness, especially the first order ones. On the other hand, for a same domain the degrees of freedom of a triangulation mesh are much more than those of a hexahedral or quadrilateral mesh. If the mesh requirement of conforming the boundary is omitted, the Cartesian methods can automatically generate octree meshes that contain particular polyhedron elements without creating surface meshes. But the mesh created by this method has a poor accuracy on boundary vicinity leading to an inaccurate boundary condition. At the moment, for complex geometries only the polyhedron mesh generation technique can be fully automatic with few human interventions.
Polyhedron meshes such as octree meshes have been widely used in fluid mechanics exhibiting many topological properties that make them highly and automatically adaptable to the complex shape and smooth transition of mesh density 9,10 . Nevertheless, the application of polyhedron elements for continuum mechanics is few because the topology of a polyhedron is inconstant. Hence it is impossible to construct analytic shape functions for arbitrary polyhedron elements with varying topologies. Some kind of interpolation schemes must be applied. Two kinds of interpolation schemes were developed. The first one was based on natural neighbor interpolations proposed by Sibson 11,12 . Traversoni 13 and Sambridge et al. 14 introduced Sibson interpolations modified by Watson 15 to the Galerkin-type natural neighbor, which was extended to continuum by Sukumar et al. 16 29 . However, the interpolation scheme generally needs considerable calculation to get adequate accuracy.
Another approach to constructing polyhedral elements is to utilize SBFEM. In mid-1990s, SBFEM also known as the consistent infinitesimal finite-element-cell method was bought up by Wolf and Song as an extension to FEM to simulate the wave propagation in the unbounded domain [30][31][32][33] . As a similarity-based semi-analytical numerical approach for PDE, SBFEM has many advantages over the traditional FEM and Boundary Element Method (BEM) exhibiting high accuracy and flexibility 34,35 . SBFEM just needs boundary discretization reducing the problem spatial dimension by one with no need for a fundamental solution and the solution along the radial direction is in a precise closed-form. For the unbounded media, SBFEM can perfectly satisfy the boundary condition of the radiation damping. The unit-impulse response matrix in the time domain and the dynamic-stiffness matrix in the frequency domain of the unbounded media can be directly obtained, which is equivalent to a time and space coupling artificial non-reflecting boundary 36 . For the bounded media, the SBFEM is derived in a wedge media. The wedge media can be assembled into an arbitrary convex polyhedron (3D) or polygon (2D) super element revolving around the similarity center. If this assembly is not occlusive and the crack tip is positioned in the similarity center, the singularity of the field function can be directly described by the analytical solution along the radical crack surface. This superiority of mesh adaptability can be combined with the automatic mesh generation technique based on polyhedron elements.
The first try was a SBFEM formulation for arbitrary polyhedral elements and a simple method to generate polyhedral meshes based on octree presented by Talebi et al. 37 . To automatically generate an octree polyhedral mesh from the Standard Tessellation Language (STL) widely used in the 3D printing and Computer Aided Design (CAD), Liu et al. 38 came up with a two-steps method. The first step is to generate an octree mesh within the domain enclosed by the surface defined in the STL. The second step is to cut the polyhedral elements with the surface. After the refinement of the second step, the poor boundary accuracy of the octree mesh is improved. Furthermore, Zhang et al. 39 applied the SBFEM polyhedral elements to the non-matching meshes by adding nodes to the elements adjacent to the interface. By now, there is a slight limitation when using SBFEM polyhedral elements. Since the boundary is discretized by FEM elements, the facets of a polyhedral element must be a triangle or quadrilateral, which may lead to a facet division process. To resolve this limitation of SBFEM, Ooi et al. 40 derived a dual SBFEM formulation over arbitrary faceted star convex polyhedron. On the other hand, Zou et al. 41 brought up a FEM polyhedral element formulation with shape functions derived from SBFEM circumventing this limitation as well. The automatic polyhedral element based on SBFEM was applied to dynamic nonlinear analysis of hydraulic engineering 42 , concrete fracture modelling 43 , damage modelling 44,45 , Soil-Structure Interaction (SSI) analysis of NPP 46 , etc. These applications of SBFEM polyhedral meshes were mostly implemented by self-developed programs. Few of them integrate the SBFEM into the universal FEM software such as ABAQUS 43,47 . And these integrations is suitable for traditional FEM elements with fixed topologies only. The integration of SBFEM and the universal FEM software can promote each other. Implementation of arbitrary polyhedral elements in ABAQUS is still in need.
In this paper, arbitrary polyhedron elements based on SBFEM for three dimensional continuum is developed and implemented in ABAQUS facilitating automatic dynamic analyses of structures. The governing equation of the bounded 3D elastic continuum in a polyhedron domain is derived from the SBFEM in "SBFEM governing equation for the elastic continuum" section. The FEM-coupled characteristic matrix of a polyhedron element can be obtained from the solution of the governing equation. The topology feature and corresponding data structure of a polyhedron is analyzed in "Topology and data structure of an arbitrary polyhedron element" section. The user element subroutine of UEL is developed associated with the Hilber-Hughes-Taylor (HHT) implicit time integration scheme in "Implementation in ABAQUS' section. In "Numerical application" section, there are two examples for the UEL. A cantilever subjected to a harmonic excitation with the traditional hexahedron element and the SBFEM-based polyhedron element are compared to confirm the accuracy of the UEL. Then taking a www.nature.com/scientificreports/ practical example of the soil-structure dynamic interaction analysis of a nuclear power plant, the applicability and performance of this implementation of polyhedron elements is tested.

SBFEM governing equation for the elastic continuum
SBFEM is a numerical method for solving PDE. The PDE is transformed into ordinary differential equations after specific spatial discretization, which is different from FEM, BEM, etc. SBFEM needs boundary discretization only decreasing the spatial dimension by 1 and remains analytical in the radial direction without the need of a foundation solution. The derivation of the SBFEM governing equation of elastic continuum is based on a wedge area defined as W-element. The governing equations of W-elements that shared the same scaling center can be assembled by degrees of freedom of adjacent W-elements (see Fig. 1) leading to a piece-wise boundary of the discretized domain defined as an S-element. The choice of the scaling center location can be pretty casual. The boundary of the discretized domain must be visible from the scaling center is the only geometry requirement that needs to be satisfied. The bounded domain can be divided into multiple S-elements to meet the requirement of the scaling center or to improve the poor geometry of W-elements (see Fig. 2). Apparently, the geometry of the S-element can be rather flexible for any convex polyhedron with only triangle or quadrilateral facets. The detailed derivation of SBFEM governing equation is well documented in many literatures 33,[48][49][50][51][52] . For the sake of integrity, only key equations are given below. The equations in this section are given in linear algebraic notation for the convenience of reading and programming. The brace symbol of {·} denotes a matrix and the bracket symbol of [·] denotes a column vector.    The equilibrium equations (Eq. (1)), the geometric equations (Eq. (2)) and the constitutive equations (Eq. (3)) can be combined into one set of partial differential equations of 3D dynamic elastic mechanics with the displacement vector filed as the unknown variable as where u x, y, z = u x x, y, z u y x, y, z u z x, y, z T denotes the displacement vector filed, [D] denotes the material property matrix, p denotes the body forces except the d' Alembert inertial force, ρ denotes the density, [L] denotes the spatial partial derivative operator. The linear elastic material property matrix [D] is given by where denotes the Lame parameter, E denotes the Young's modulus, G denotes the shear modulus, ν denotes the Poisson's ratio. Equation (7) represents a set of second order linear nonhomogeneous partial differential equations with constant coefficients in Cartesian coordinate, namely Navier equations. One way to solve the Navier equations analytically is to transform them into wave equations of potential functions by Helmholtz's theorem. The other way is to solve the Navier equations numerically by spatial and time discretization.
Two kinds of boundary conditions and one initial condition should be introduced to get the definite solution of the fundamental equation Eq. (7). One boundary condition is the geometry boundary condition on the surface S u defined as [L]  Coordinate mapping based on scaled boundary transformation. Taking one W-element as an example (see Fig. 3), the boundary is discretized with one plane four nodes linear finite element with circumferential axis of η and ζ (see Fig. 4). Other 2D finite elements such as the quadratic element with eight nodes can be used for this boundary discretization as well 51 . The origin of the global Cartesian coordinate system of the nodes needs to be translated to the scaling center by where {x e } , ỹ e , and {z e } denote the global coordinates of the nodes, x 0 e , ỹ 0 e , and z 0 e denote the coordinates of the scaling center in global Cartesian coordinate system, {x e } , y e , and {z e } denote the translated global coordinates of the nodes.
The mapping function between the global coordinates and the local coordinates of the points on the boundary is given by � n � x, y, z �� =   n x 0 0 n y 0 n z 0 n y 0 n x n z 0 0 0 n z 0 n y n x   , (13) u x, y, z, t = u x, y, z when t = 0, The new local scale boundary coordinate is built for this W-element with the origin at the scaling center. Points on the plane finite element is scaled along the radical direction from the scaling center to the boundary, which is defined as the ξ axis of the local scale boundary coordinate. The mapping function between the global coordinates and the local coordinates of the points in the W-element is given by where x(ξ , η, ζ ) , y(ξ , η, ζ ) , and z(ξ , η, ζ ) denote the local coordinates of the points in the W-element.
Equation (17) provides a mapping relation between the local and global coordinates. Based on the derivative method for compound function, the corresponding relationship of the derivative operators between two coordinates is given by where [J] is the Jacobian matrix expressed as The variable ξ can be isolated from Eq. (19) as with the Jacobian matrix on the boundary Ĵ defined as where (·) ,η denotes differentiating with respect to η , (·) ,ζ denotes differentiating with respect to ζ.
Based on the space analytical geometric theory, the geometry properties in the local coordinate system can be derived. The infinitesimal line vectors in the axis direction are given by where (·) denotes a vector.
The infinitesimal volume is given by   (7), the spatial partial derivative operator [L] in local coordinate system is given by  (2) and (3), the stress {σ } in local coordinate system is given by Considering the surface traction boundary condition in Eq. (10) to Eq. (12), the surface traction on the coordinate surface is defined as where t ξ , {t η } , and t ζ denote the the surface traction on the coordinate surface of ξ , η , and ζ respectively.
Compared with Eq. (31) to Eqs. (33), (36) to Eq. (38) can be reformulated as Spatial discretization of the displacement field on the boundary. Taking the same W-element in "Coordinate mapping based on scaled boundary transformation" section as example (see Fig. 5), the boundary is discretized with a finite element with circumferential axis of η and ζ (see Fig. 4).  w j T is defined as .

Topology and data structure of an arbitrary polyhedron element
The polyhedron geometry is a 3D shape with straight edges, sharp corners and vertices, and any number greater than 2 of plane polygonal faces with any number greater than 2 of straight edges and vertices. The tetrahedron, pentahedron, hexahedron and heptahedron that widely used in FEM are just the special cases of the convex polyhedron. Vastly different from the FEM element geometry, the first step to construct a polyhedron element is to describe the uncertainty of the topology of varied kinds of polyhedron elements in a model with a uniform data structure that can be recorded in the intermediate file for compatible purpose.
Essential information to describe the topology of an arbitrary polyhedron. For example, Fig. 6 shows a polyhedron element that may be encountered in an octree mesh, which has nine facets, twenty edges and thirteen vertices. The coordinates of the vertices are the basic information as shown in Table 1. In addition, the connectivity of the vertices of each facet needs to be defined as shown in Table 2. These two sets of information are sufficient and necessary for the description of the topology of any polyhedron. It is worth to notice that the order of the vertices and facets of a polyhedron is insignificant and when there are multiple polyhedra, there is no need to record the correspondence between the local and global vertex ID like FEM.
SBFEM discretization of the polyhedron. The geometry of the S-element can be adapted for any convex polyhedron with only triangle or quadrilateral facets which is a special case of the polyhedron defined in "Essential information to describe the topology of an arbitrary polyhedron" section. Nevertheless, the polygonal facet of a polyhedron can always be converted into multiple triangle or quadrilateral facets by adding a vertex on this facet (see Fig. 6). After this treatment of the polyhedron facets, a point inside the polyhedron that can see all the points on the facets is chosen as the scaling center and the vertices are defined as nodes of this polyhedron domain. The scaling center and the nodes of each facet constitute a W-element. All these W-elements are assembled into an S-element representing this polyhedron domain. It is obvious that the number of W-elements is equal to the number of facets of the polyhedron and the number of the nodes of the S-element is equal to the number of vertices of the polyhedron. The data structure to describe any 3D SBFEM polyhedron element is summarized in Table 3.     38 . Besides, at present OPENFOAM, FLUENT and STAR-CCM + provide robust polyhedron mesh generators, the latter two of which are the most convenient. The general procedure to generate a polyhedron mesh from geometry with an out-of-box tool is illustrated in Fig. 7. Because the common FEM software for solid media does not support the polyhedron element, the preprocessor built in them cannot generate polyhedron meshes. Furthermore, the raw polyhedron mesh needs some extra processing steps such as classification, grouping, adding vertices before applied. The flow of geometry and mesh information between various software through files in compatible formats is inevitable. The Initial Graphics Exchange Specification (IGES) for the geometry information and Visualization Toolkit (VTK) for the mesh information are recommended. It should be noted that no matter which polyhedron mesh generator above is chosen or developed, it is required to store the polyhedron mesh in the intermediate file that will be imported into the UEL after the polyhedron mesh is obtained.

Implementation in ABAQUS
The SBFEM polyhedron element is introduced in ABAQUS/Standard as a user element. The user element is implemented on the element iteration level of the solver in the form of a FORTRAN subroutine named UEL. During each increment of the computation, ABAQUS will invoke the UEL when calculating the coefficient matrix of each user element on and on until the end of the analysis. Assembly of the global coefficient matrix, nonlinear iteration, result output, imposition of contact, constrain, boundary and load conditions are performed by ABAQUS automatically.

Basics of UEL in ABAQUS.
UEL is essentially a function to compute specific variables with given variables. The subroutine name, declaration of the input and output variables are fixed (see "Appendix 1" 54 ). During the calculation of each element, the UEL will be called importing information of analysis type, time increment, element type, right-hand-side vectors, etc. to export the mass matrix, stiffness matrix, equivalent nodes force vector, solution-dependent state variables, etc. There are some conventions of user subroutines in ABAQUS that should be emphasized. All user subroutines and the subroutines created by users must be stored in a single source code file or its compiled object file. By default, the user subroutine can call another subroutine that is built in FORTRAN compiler, created by users or provided by ABAQUS named utility routines. This restriction suggests that some complicated calculations in a user subroutine have to be implemented by source codes provided by users. And if there were multiple user elements, these user elements should be implemented by the same subroutine of UEL with conditional branch statements on element types. Although user subroutines cannot call each other, the data can be transferred through the user-defined COMMON block variables between them.

Associated INP and intermediate file.
A convenient characteristic of ABAQUS is that the entire information of a FEM model is stored in an INP file, which is human-readable. This INP file is the only input required by the solver of ABAQUS. The user elements must be declared in the INP file before the calculation (see "Appendix 2" 54 ). The declaration of user elements provide the information of the element type, element set name, number of nodes, node connectivity, node degrees of freedom and element properties such as material parameters, etc. According to the format of the user element declaration in an INP file, the user elements need to be classified Table 3. Data structure of the SBFEM polyhedron element. a The node sequence should be in a certain order pointing out of the region according to the right-hand rule.

Element ID k
Coordinate of scaling center x 0 , y 0 , z 0

Number of nodes n
Local node ID Global node ID Node coordinates www.nature.com/scientificreports/ into different types based on at least the number of nodes and the material properties and the node connectivity of polygonal facets of a polyhedron cannot be represented. Due to the limitation of the INP file mentioned before, an intermediate file storing the topology information and material properties of arbitrary polyhedron elements is necessary. This file should be imported once the calculation is started and when iterating through the elements, the corresponding topology information of current user element should be introduced into UEL. With this intermediate file, any polyhedron mesh generator mentioned in "Generation of the SBFEM polyhedron mesh" section can be used and any mesh with arbitrary convex polyhedron shape can be calculated.

Algorithm of UEL for SBFEM polyhedron element. Although UEL has provided a specific interface,
it is necessary to understand the calculation process of ABAQUS/Standard to get a clear picture of what information is passed in UEL and what information UEL should return under different circumstances. As shown in Fig. 8, calculation process of ABAQUS/Standard has a hierarchy of three levels. Level 1 denotes the calculation procedure of steps in an analysis. Level 2 denotes the calculation procedure of increments of a step. Level 3 denotes the actual calculation of the polyhedron element matrix and vector in an increment by calling UEL for each element (see the green block in Fig. 8). The utility routine of UEXTERNALDB plays an important role in Level 1 and 2 (see the red blocks in Fig. 8), providing adequate interventions into the different phases of the calculation process, which give us more freedom to realize more complex functions. To coincide with UEL, the intermediate file storing the topology and material information of user polyhedron elements is imported by subroutine UEXTERNALDB at the beginning of the analysis, which is then transferred into UEL by COMMON block variables when the calculation of the polyhedron element starts.
The behavior of UEL is controlled by the basic branching structure (see Fig. 9) judged by the variable of LFLAGS. Considering only the general static and direct-integration dynamic analysis, there are five different cases (see colored blocks in Fig. 9) requiring returned values respectively according to the Newton-Raphson Method (NRM) and HHT direct time integration 55 .
The polyhedron element computation algorithm of UEL based on SBFEM is summarized in "Appendix 3". The assembly of the coefficient matrix of W-elements in STEP 2.3 and stiffness and mass matrix calculation that www.nature.com/scientificreports/ include Eigen-decomposition are the major differences from FEM. Subroutines from LAPACK 56 is appended to UEL in the code file to perform STEP 3. It is obvious that this UEL is compatible with the traditional FEM element topology built in ABAQUS, considering the tetrahedron, pentahedron and hexahedron as special cases of the polyhedron.

Numerical application
There are two examples for the test and application of the implementation of arbitrary polyhedral elements developed in "Implementation in ABAQUS" section. A cantilever subjected to a harmonic excitation with the traditional hexahedron element and the polyhedron element are compared to confirm the accuracy of the developed UEL. Then the polyhedron element based on SBFEM is applied in the soil-structure dynamic interaction analysis of a nuclear power plant, the applicability and performance of this implementation is tested.
Cantilever subjected to a harmonic excitation. A 5 m long horizontal beam subjected to concentrated forces on the vertices at one end and fixed at the other end is considered (see Fig. 10) www.nature.com/scientificreports/ where T denoting the period is assuemd to be 0.05 s, A denoting the amplitude is assumed to be 1000 N. F denotes the concentrated force. The HHT direct time integration calculation is terminated at 0.25 s with a fixed time step of 0.001 s.
Numerical model. From the geometry, two sets of meshes with varied mesh density are generated. The first set of meshes is made of structured hexahedron elements, namely C3D8 element (3D linear solid element with 8 nodes) built and verified in ABAQUS/Standard and the other is made of unstructured polyhedron elements implemented by UEL based on SBFEM. The details of the meshes are summarized in Table 4, where label P1 to P4 denote the polyhedron meshes with different mesh density (see Fig. 11) that can only be calculated by UEL developed in "Implementation in ABAQUS" section and H1 to H7 denote the hexahedron meshes (see Fig. 12) that can be calculated by C3D8 (built-in element) or the UEL (labeled as H2_UEL and H6_UEL). Model analysis for the first 50 modes of vibration and implicit dynamic analysis are performed.  www.nature.com/scientificreports/ Results and discussion. Model analysis checking the stiffness and mass matrix with specific boundary conditions is the most convenient approach to verify the user element for dynamic analysis, excluding the influence of load conditions and nonlinearity. The results of the first 50 natural frequencies and the displacement response of the dangling end (see Fig. 10) are compared to check the accuracy and efficiency of the UEL for the polyhedron element. As the increasing of the number of elements, the results of model analysis will converge to the true value of this problem. The first 50 natural frequencies of model H1 to H7 are listed in Fig. 13, from which the convergence values of the first 50 natural frequencies of the cantilever are obtained covering more than 95% effective mass in the main direction (97.1% in X and Y direction, 95.1% in Z direction). Figure 13 shows that the difference between the frequencies of models with varied mesh densities is no more than 9% for the first vibration mode and this difference is enlarged to less than 19% for the first 9 vibration modes, which tends to keep on increasing with the order. Even though when the mesh is coarse, the convergence speed is fast with the increasing of the element, the convergence speed is much slower near the convergence values leading to rather poor computation efficiency (see Fig. 14). For example, the highest frequency of model H6 and H7 are 1374.8 Hz and 1392 Hz, the computation times of which are 5 s and 1109 s. More than 220 times the computation time of model H6 brings less than 1.2% improvement of the results of model H7. This is a typical illustration of the need to pay great attention to the balance between mesh density and calculation efficiency. Since it is reasonable to assume that the results of model H7 are close enough to the convergence values, the true fundamental frequency of this cantilever is assumed to be 19.955 Hz.
To exclude the influence of the meshes, the same hexahedron meshes H2 and H6 are calculated by C3D8 and UEL respectively. Figure 15 shows the frequency comparison of the cantilever model H6_UEL calculated by UEL and H7 by C3D8. The consistency of the results of model H7 and H6_UEL confirms the accuracy of the UEL for the SBFEM-based polyhedron element with error less than 0.6% over the entire calculated frequency range. Figure 16 shows the frequency comparison of the cantilever mesh model H2 calculated by C3D8 and UEL. Compared to model H2 calculated by C3D8, model H2_UEL can get much better results that is almost equal to model H6 calculated by C3D8. In other words, polyhedron element based on SBFEM can get the same level of accuracy with less than 1% of the elements based on FEM, the most immediate advantage of which is the reduction in storage space for large-scale problems.
The displacement responses in Y direction from the implicit dynamic analysis of model P1, P2, P3, and H6 are demonstrated in Fig. 17. Except the displacement values that are close to zero, the error between the model calculated by the C3D8 element and the user-defined polyhedron element based on SBFEM is less than 0.5%. The results of model analysis and direct dynamic analysis of the cantilever have proven the accuracy of the UEL implementation of the SBFEM-based polyhedron element and the efficiency of the SBFEM for the solid medium.
Soil-structure dynamic interaction analysis of a nuclear power plant. To verify the performance and robustness of the implementation that has been primarily validated in "Cantilever subjected to a harmonic excitation" section, the application in the analysis of a practical engineering problem is necessary. A typical numerical elastic SSI analysis of a NPP with direct method in the time domain is taken as an example, which includes the simulation of the structure with traditional FEM elements, near field soil with the polyhedron elements, far field soil with an artificial non-reflecting boundary condition, external excitation input with site response analysis, and interaction between them with constrain or contact. For simplicity, the SSI analysis of a NPP is elastic. But it is enough to demonstrate the merit in automatic mesh generation and ability to be integrated seamlessly with other modules of the polyhedron element based on SBFEM and its implementation by UEL. The nuclear power plant is located on elastic rock soil with the Young's modulus of 2.63E10 Pa, Poisson's ratio of 0.25 and density of 1762.06 kg/m 3 . Two artificial seismic plane waves vibrating along the two horizontal directions with PGA of 0.9806 m/s 2 or 0.1 g (see Fig. 18) defined by the U.S. NRC RG 1.60 57 response acceleration spectrum (see Fig. 19) are introduced into the SSI system from the bottom of the rock propagating upward along the vertical direction. This earthquake intensity just satisfies the minimum limit recommended by the www.nature.com/scientificreports/ IAEA guidelines for the SL-2 ground motion hazard level 58 . The SSI numerical analysis model of NPP based on the direct method is illustrated in Fig. 20.
Numerical model. The SSI numerical analysis of NPP can be divided into three individual parts: the structure, the near-field and the far-field. The structure and the near-field are modeled by FEM. The far-field is modeled by the direct method. Benefiting from the FEM for unmatched mesh based on kinematic constraints, the structure and the near-field soil can be meshed separately.     www.nature.com/scientificreports/ For simplicity, the structures on the nuclear island including the containment and the auxiliary buildings are meshed with shell elements S4R or S3 (see Fig. 21). The containment of this NPP is a double shell structure including an outer reinforced concrete vessel called shield building (SB) and an inner steel vessel. The foundation of the nuclear island is meshed with solid elements embedded in the near-field soil. The connection between the structure and the near-field soil is modeled by the tie constrain defined on the interface.
The geometry of the near-field soil can be very complex once geological features must be considered like cracks or faults, etc. Therefore, the solid element is the most reasonable choice to model the near-field soil. Besides, in order to ensure the accuracy, the direct method requires that the size of the near-field soil should be as large as possible and no less than twice of the size of the structure consuming a lot of computing resources. Tremendous efforts for the mesh generation are needed to balance the accuracy and efficiency of the numerical model. In this example, the near-field soil is meshed with octree solid elements automatically (see Fig. 22). Although the geometry of the near-field soil in the example is rather simple, the geological complexity of the soil should not cause any difficulties in the automatic mesh generation process.
The SSI effect is considered by applying an artificial non-reflecting boundary condition (NRBC) at the truncation of the near field soil. Specifically, a set of springs and dashpots tied to the ground along the three directions are added to the nodes on the interface between the near-field and far-filed soil representing the approximation of the viscoelastic NRBC 59 . Unlike the substructure method, the corresponding input of exogenous excitations for a specific NRBC is needed. For viscoelastic boundary, the seismic load is imported by equivalent forces acting   Table 5. A transient dynamic analysis based on HHT implicit time integration is performed with a fixed time step of 0.01 s.
Results and discussion. To investigate the detailed acceleration and displacement responses of the NPP, five nodes are selected as monitoring points (see Fig. 23), the locations of which are listed in Table 6.
The time history of acceleration and displacement response of these nodes is compared. Figures 24 and 25 show the horizontal displacement and acceleration response of the nodes on the NPP foundation at the height of 33.576 m. It appears that the foundation is moving like a rigid plate without any rotation. The peak acceleration of the input movement at the foundation is about 0.09 g, which is a little less than the input movement at the substratum top (or the hard rock bottom). The site hardly changes the input movement due to the thin hard rock soil layer beneath the foundation indicating an insignificant SSI effect. The peak displacement of the nodes at the foundation is about 0.076 m, which can be recognized as a basic motion of the entire NPP.    www.nature.com/scientificreports/ structure. In Fig. 27, the same trend that the peak accelerations increase from 0.85 m/s 2 (X direction) and 0.81 m/ s 2 (Y direction) to 2.75 m/s 2 (X direction) and 3.43 m/s 2 (Y direction) to 5.47 m/s 2 (X direction) and 7.9 m/s 2 (Y direction) can be found. The geometry asymmetry is the main reason why there are differences in the response along the X and Y direction. The acceleration enlargement ratio from the foundation to the top is about 9.75, which is the result of the stiffness of the soil. Figure 28 illustrates the relative displacement response of the SB top to the SB foundation, eliminating the rigid transition movement due to the global site response. The peak values of the displacement in X direction and Y direction are 1.6 cm and 2.0 cm. The floor response spectra of monitoring points at different heights with the damping ratio of 0.05 are listed in Fig. 29. The same trend of spectrum enlargement for the two horizontal directions can be observed. The peak spectrum values of the three nodes appear at the same frequency of 2.78 Hz, which is near the fundamental   www.nature.com/scientificreports/ www.nature.com/scientificreports/ (a) X direction (b) Y direction (a) X direction ( b) Y direction  www.nature.com/scientificreports/ frequency of the SSI system. In the frequency range between 0.25 and 1.0 Hz, the enlargement ratio along the height is less than 1.5. The enlargement ratio increases starting from 1.0 Hz and reaches its peak value of 13.33 at 2.78 Hz. This value decreases to about 9.3 at 100.0 Hz that is nearly equivalent to the acceleration enlargement ratio from the foundation to the top of the shield building. The distribution of the peak acceleration (absolute values) can be seen from Fig. 30. The maximum PGA values arise at the top of the shield building. In general, the PGA values are enlarged along the vertical direction except regions near the foundation. The amplification ratio of the peak acceleration along three paths is detailed (a) X direction ( b) Y direction  www.nature.com/scientificreports/ in Fig. 31. The nodes of the three paths that are parallel to Z axis are chosen as the monitoring points (see Fig. 23) and the height of the foundation bottom is set as zero. There is an upheaval of amplification ratio at the height of around 75 m where the water tank is located.

Conclusion
The results of two examples indicate that the implementation of SBFEM-based arbitrary polyhedron elements can be adapted to the automatic FEM analysis of static and dynamic elastic problems with polyhedron meshes generated by any tool. The first example conforms the accuracy and efficiency of the polyhedron element based on SBFEM and the associated UEL implementation in ABAQUS. Through the second example, it can be found that (a) X direction (b) Y direction