The Implementation of Spectral Element Method in a CAE System for the Solution of Elasticity Problems on Hybrid Curvilinear Meshes

Modern high-performance computing systems allow us to explore and implement new technologies and mathematical modeling algorithms into industrial software systems of engineering analysis. For a long time the finite elementmethod (FEM)was considered as the basic approach to mathematical simulation of elasticity theory problems; it provided the problems solution within an engineering error. However, modern high-tech equipment allows us to implement design solutions with a high enough accuracy, which requires more sophisticated approaches within the mathematical simulation of elasticity problems in industrial packages of engineering analysis. One of such approaches is the spectral element method (SEM).The implementation of SEM in a CAE system for the solution of elasticity problems is considered. An important feature of the proposed variant of SEM implementation is a support of hybrid curvilinear meshes. The main advantages of SEM over the FEM are discussed. The shape functions for different classes of spectral elements are written. Some results of computations are given for model problems that have analytical solutions. The results show the better accuracy of SEM in comparison with FEM for the same meshes.


Introduction
The finite element method (FEM) [1,2] was considered as the main approach for solution of the problems of elasticity theory taking into account finite deformations.The desire to find the stress-strain state of structures with high accuracy is still a hot topic at the moment and it forces us to seek new and effective methods for solving engineering problems.One of such approaches is the spectral element method (SEM) [3][4][5].The SEM was firstly applied for the modeling of liquid flows [6].These problems require a high accuracy and high rate of computations.Later the SEM was successfully adopted for seismic problems [3,7,8].A special quadrature formula was constructed for integration over space.This quadrature formula permits one to develop a fully explicit scheme for the integration over time, which is an important advantage of the SEM.
Among the main advantages of the SEM over the FEM is the high accuracy of approximation of the solution at a substantially smaller number of mesh elements required.The error of the numerical solution of SEM decreases exponentially with the order of elements [9].When working with a model, the user does not need to rebuild and refine the mesh to verify the mesh convergence [9] of the obtained solution, as it has been when using the FEM, since with the use of the SEM the mesh can remain initial with only changing the order of elements.The possibility of effective paralleling of computing systems with shared and distributed memory using OpenMP and MPI technology makes SEM attractive for industrial applications in various software systems [10].In particular, a spectral element finite element scheme that efficiently solves elliptic problems on unstructured hexahedral meshes is developed in [11].It is demonstrated that problems with over 50 million degrees of freedom can be solved in a few seconds on an off-the-shelf GPU.
The industrial application of the SEM is hindered by the problem of mesh generation [12,13].Typically, for a multibody model geometry it is quite difficult to build a conformal finite element mesh consisting of hexahedral elements only for which the SEM [3] was initially developed.In general, the three-dimensional finite element mesh contains inmost dominantly hexahedral and tetrahedral elements.The prismatic and pyramidal elements are used in order to connect triangular and quadrilateral surfaces of elements.A two-dimensional mesh may contain rectangular and triangular elements.Such finite element meshes are called hybrid (mixed) meshes [14]; the SEM was adopted for them within the present work.
This article discusses the variant of implementation of spectral element method on hybrid curvilinear meshes for three-dimensional problems of elasticity theory and its industrial application in CAE Fidesys [15].The shape functions for different classes of spectral elements are written.Some results of computations are given for model problems that have analytical solutions.

Materials and Methods
Let the domain Ω ∈ R 2 (R 3 ) be divided into n e elements Ω e .The classes of quadrilaterals and triangles serve as elements for two-dimensional case.The classes of hexahedrons, tetrahedrons, pyramids, and prisms serve as elements for threedimensional case.Each element Ω e is defined by reference points.The number of reference points for a particular element is indicated in its name.
Each reference point is determined by the index a varying from 1 to d, where d is a number of reference points of the finite element.A nondegenerate mapping F e : Λ → Ω e of the base (reference) element Λ in Ω e is built for each element.The coordinates of points in the global coordinate system x e = (x e , y e , z e ) and in the reference coordinate system  = (, , ) are related as follows: where x e = (x e , y e , z e ) is the vector of global coordinates,   ∈ Ω  , x e a the global coordinates of the reference point a of the element Ω e ,  = (, , ) the reference coordinates, and N a () the ath shape function of the finite element,  ∈ Λ.
Generally, the Jacobian [16] of the mapping F e () is required for the calculation of the integral of an arbitrary function h(x) over the element Ω e through reference coordinates : where h e (x) = h(x)| Ω e and J e () = det J e () is the Jacobian in the reference point .The Jacobian matrix J e () can be calculated in a standard way: In the three-dimensional case: J e () = ( In the two-dimensional case, we assume that the finite elements lie in a plane (x, y) and reference points lie in a plane (, ):

QUAD:
The Class of Quadrilaterals.The class of quadrilaterals includes the following types of finite elements, named by number of reference points: four-node, QUAD4 ( = 4); eight-node, QUAD8 ( = 8); nine-node, QUAD9 ( = 9).The reference element Λ for the class of quadrilaterals The nodes of a quadrilateral spectral element are the Gauss-Lobatto-Legendre (GLL) nodes.Let  be an order of the spectral element; then the number of nodes in the element  QUAD  = ( + 1) 2 .In one-dimensional case the GLL-nodes are calculated as roots   of the derivative of the Legendre polynomial   of an order , which can be defined as   () = (1/2  !)(  /  )( 2 − 1)  .In two-dimensional case the coordinates of GLL-nodes are defined as a direct product of one-dimensional coordinates.
The shape functions for the element Ω  are built based on the direct product of one-dimensional Lagrange polynomials: is the th shape function of the order .The one-dimensional Lagrange polynomials ℓ   of the order  are defined as follows [3]: Each polynomial ℓ   is equal to 1 at its node  and is equal to 0 at the remaining nodes of the element Ω  : ℓ   (  ) =   , where   is the Kronecker symbol.
The approximation of solution of the problem u e = u| Ω  on a quadrilateral spectral element Ω  will be as follows ( =  + ): In order to integrate an arbitrary function ℎ(x) over the element Ω  , the Gauss-Lobatto-Legendre (GLL) quadrature formula is used: where w  are the GLL-weights and  l are the coordinates of GLL-nodes.GLL-weights are calculated as follows: The coordinates of GLL-nodes and GLL-weights (indexes ,  vary from 0 to ) can be written as follows: One of the most important features of the SEM is that the computation of the Lagrange polynomials that are used for the approximation of solution is based on the same GLLnodes that are necessary to calculate integrals over the region Ω  using the Gauss-Lobatto-Legendre quadrature formula.

TRI:
The Class of Triangles.The class of triangles contains the following types of finite elements named by the number of reference points: three-node, TRI3; six-node, TRI6.The reference element Λ for the class of triangles The nodes of a triangular spectral element are the points  j , which are obtained from the solution of the electrostatic problem described by Hesthaven and Teng [17] and Taylor et al. [18].Let  be an order of the spectral element; then the number of nodes in the element  TRI  = ( + 1)( + 2)/2.Shape functions for the element Ω  are built based on an orthogonal basis ( is the order of a spectral element): where  ,  () are the Jacobi polynomials of the order  on the interval [−1, 1], orthogonal with the weight function (1 − )  (1+)  .The Vandermonde matrix elements are calculated in the nodes  j of the element  TRI , =  TRI  ( j ),1 ≤ ,  ≤  TRI  .Then the shape functions will be as follows: The approximation of solution of the problem u e = u| Ω  on a triangular spectral element Ω  will be as follows: To integrate an arbitrary function ℎ(x) over the element Ω  , the symmetric quadrature formulas, described by Zhang et al. [19], are used: where w  are positive weights,  l are the coordinates of points, and  max is the number of points in the quadrature formula corresponding to the spectral element of the order .

HEX:
The Class of Hexahedrons.The class of hexahedrons contains the following types of finite elements named by the number of reference points: eight-node, HEX8; twenty-node, HEX20; twenty-seven-node, HEX27.The reference element Λ for the class of hexahedrons The nodes of quadrilateral spectral element are the GLLnodes.Let  be an order of the spectral element; then the number of nodes in the element  HEX  = ( + 1) 3 .The shape functions for the element Ω  are built based on the direct product of one-dimensional Lagrange polynomials:  HEX , () = ℓ   ()ℓ   ()ℓ   (), th shape function of the order .
The approximation of solution of the problem u e = u| Ω  on hexahedral spectral element Ω  will be as follows ( =  2 +  + ): To integrate an arbitrary function ℎ(x) over the element Ω  , the Gauss-Lobatto-Legendre (GLL) quadrature formula is used: where w  are the GLL-weights and  l are the coordinates of GLL-nodes.
Coordinates of the GLL-nodes and GLL-weights (indexes , ,  vary from 0 to ) The nodes of a tetrahedron spectral element are the points  j , which are obtained from the solution of the electrostatic problem described Hesthaven and Teng [17].Let  be an order of the spectral element; then the number of nodes in the element  TET  = ( + 1)( + 2)( + 3)/6.The shape functions for the element Ω  are built based on the following orthogonal basis ( is the order of a spectral element): where  ,  () are the Jacobi polynomials of the order  on the interval [−1, 1], orthogonal with the weight function (1 − )  (1 + )  .The Vandermonde matrix elements are calculated in the nodes  j of the element  TET , =  TET  ( j ), 1 ≤ ,  ≤  TET  .Then the shape functions will be as follows: The approximation of solution of the problem u e = u| Ω  on a tetrahedral spectral element Ω  will be as follows: In order to integrate an arbitrary function ℎ(x) over the element Ω  , the symmetric quadrature formulas, described by Zhang et al. [19], are used: where w  are positive weights,  l are coordinates of points, and  max is the number of points in the quadrature formula corresponding to the spectral element of the order .
The nodes of a pyramidal spectral element are the points  j , matching the GLL-nodes on a square pyramid base and nodes for triangular spectral element on the side surfaces of the pyramid.Internal points are located in the planes parallel to the square base of the pyramid in the scaled GLL-nodes.Let  be an order of the spectral element; then the number of nodes in the element The shape functions for the element Ω  are built based on the following orthogonal basis ( is the order of a spectral element): where  ,  () are the Jacobi polynomials of the order  on the interval [−1, 1], orthogonal with the weight function (1 − )  (1+)  .The Vandermonde matrix elements are calculated in the nodes  j of the element  PYR , =  PYR  ( j ) 1 ≤ ,  ≤  PYR  .Then the shape functions will be as follows: The approximation of solution of the problem u e = u| Ω  on pyramidal spectral element Ω  will be as follows: To integrate an arbitrary function ℎ(x) over the element Ω  , the symmetric conical quadrature formulas, described by Felippa [20], are used: where w  are positive weights,  l are coordinates of points, and  max is the number of points in the quadrature formula corresponding to the spectral element of the order .

WEDGE:
The Class of Prisms.The class of prisms contains the following types of finite elements named by the number of reference points: six-node: WEDGE6; fifteennode, WEDGE15.The reference element Λ for the class of prisms The nodes of a prismatic spectral element are the points  j arising from the direct product of the nodes for the triangular spectral element and one-dimensional GLL-points.Let  be an order of the spectral element; then the number of nodes in the element  WEDGE  =  TRI  ( + 1).The shape functions for the element Ω  are built based on the product of the shape functions for the triangular spectral element and the Lagrange polynomials: The approximation of solution of the problem u e = u| Ω  on a triangular spectral element Ω  will be as follows ( =  TRI  + ): To integrate an arbitrary function ℎ(x) over the element Ω  , the quadrature formulas are used, which are the unions of quadrature formulas for triangular spectral elements and the Gauss-Lobatto-Legendre quadrature formula: where w  are weights,  l are coordinates of points, and  max is the number of points in the quadrature formula corresponding to the spectral element of the order .

Features of Implementation.
The main steps of application of the spectral element method to the problems of the theory of elasticity are similar to the steps of problems solving using the finite element method, such as discretization of equilibrium equations in the integral form; selection of quadrature for calculation of integrals; building of local matrices of stiffness, mass, and damping for each element; assembling global matrices of stiffness, mass, and damping.At present, most of mesh generators for geometric models build only the finite element meshes of the first and second order, which forces us to rebuild the mesh model for calculation using the spectral element method.The most long-lasting operation at this step is to construct a graph of mesh connectivity, so one can use the connectivity graph of the initial finite element mesh model in order to speed up the algorithm.Increasing the speed of the algorithm is associated with an increase in memory consumption.In particular, for high order elements it is necessary to store the points of integration, quadrature weights, and computed values of derivatives of functions at these points.An important feature of the algorithm is the ability to use the spectral element method on any hybrid finite element mesh of the first and second order.In fact, uncoated classes of reference mesh elements currently are beam and shell elements.

Results
This approach to the implementation of spectral element method on hybrid curvilinear meshes for elasticity problems was industrially implemented in CAE Fidesys.Let us give the spectral element solution of the problem of stress analysis for a structural element and the results of spectral element analysis of wave propagation [13] in the elastic medium in comparison with the finite element and analytical solutions.Numerical experiments clearly demonstrate that the mesh convergence in the model is achieved much faster when using the spectral element method as compared with the finite element method.

Calculation of Wave Propagation in Unbounded Elastic Medium under the Action of a Point Source of Disturbances.
The wave propagation in unbounded elastic medium was analysed using the spectral element method, and the results were compared with the solution obtained via analytical calculation of displacement vector components, depending on time at fixed points of the body [1].In Figure 1, one can see the time dependencies of the displacement vector components   ,   and the relative errors  SEM  ,  SEM  of calculation results for these components obtained via the SEM against the analytical solution.The results are given for one of the receivers within the medium;  denotes time, and  is the order of approximating polynomials.As it can be seen in the graphs, the maximum error of the calculation via SEM on a mesh of about 10 thousands of elements does not exceed 2.0% at the 8th order of approximating polynomials and 0.4% at the 10th order of approximating polynomials.

Analysis of Stress-Strain State of a Bridge Span under the
Action of Pressure on Its Base.The calculation was performed for the finite elements of the first and second orders, as well as for the spectral elements of the orders from 1 to 4. An error of vector norm of maximum displacement of the bridge span was estimated.As a reference value for the assessment of the accuracy of the obtained solution, a value corresponding to the model mesh convergence was taken.Results are provided for the following cases: case 1, 10 thousand elements (the character element size is 0.8); case 2, 28 thousand elements (the character element size is 0.4); case 3, 104 thousand elements (the character element size is 0.2); case 4, 849 thousand elements (the character element size is 0.1).The case number is laid off as abscissa.
As it can be seen in Figure 2 on the graphs, the spectral element method on the coarsest mesh just on the 4th order demonstrates the accuracy of the order of 1%, whereas the finite element method on a mesh of 849 thousand elements still gives an error of more than 10% for the finite elements of the first order and about 2% for the finite elements of the second order.

Discussion
This article discusses the variant of spectral element method implementation on the hybrid curvilinear meshes for problems of elasticity theory and its industrial application in CAE Fidesys.The comparison with the finite element method was conducted, which allowed us to draw a conclusion about the high accuracy of the method and the correctness of algorithms and the program developed.In the future, it is planned to expand the implementation of the method for the cases of shell and beam elements within the static and dynamic elasticity problems.