High-Accuracy Finite Element Method : Benchmark Calculations

We describe a new high-accuracy finite element scheme with simplex elements for solving the elliptic boundary-value problems and show its efficiency on benchmark solutions of the Helmholtz equation for the triangle membrane and hypercube.


Introduction
In Ref. [1] we present a new algorithm for the calculation of high-order Lagrange and Hermite interpolation polynomials (LIP and HIP) of the simplex in analytical form, their classification and a typical example of the triangle element for high-accuracy finite element method (FEM) schemes.In this paper we show the efficiency of the schemes using high-order accuracy LIP and HIP for benchmark calculations of the exactly solvable problems for the triangle membrane and hypercube.

Formulation of the Problem
Consider the self-adjoint boundary-value problem (BVP) for the elliptic differential equation: It is assumed that g 0 (z) > 0, g ji (z) = g i j (z) and V(z) are real-valued functions, continuous together with their generalized derivatives to a given order in the domain z ∈ Ω = Ω ∪ ∂Ω ∈ R d with the piecewise continuous boundary S = ∂Ω, which provides the existence of nontrivial solutions obeying the boundary conditions (BCs) of the first (I) or the second (II) kind e-mail: gooseff@jinr.rue-mail: vinitsky@theor.jinr.ruwhere ∂Φ m (z) ∂n D is the derivative along the conormal direction, n is the outer normal to the boundary of the domain S = ∂Ω, êi is the unit vector of z = d i=1 êi z i , and (n, êi ) is the scalar product in R d .For a discrete spectrum problem, the functions Φ m (z) from the Sobolev space H s≥1 2 (Ω), Φ m (z) ∈ H s≥1 2 (Ω), corresponding to the real eigenvalues E: The FEM solution of the BVPs (1)-( 3) reduces to the stationary points of the functional The polyhedral domain Ω = Ω h (z) = Q q=1 ∆ q is divided in finite elements in the form of simplexes ∆ q with HIPs or LIPs ϕ κ r (z), calculated using the Algorithm of [1].The piecewise polynomial functions P l (z) of the order p with continuous derivatives to the order κ are constructed by joining the polynomials ϕ κ r (z) to the finite elements ∆ q ∈ Ω h (z).The expansion of the sought solution Φ m (z) in the basis of piecewise polynomial functions P l (z), Φ h m (z) = N l=1 P l (z)Φ h lm and its substitution into the variational functional (4) leads to the generalized algebraic eigenvalue problem (AEP), (A − BE h m )Φ h m = 0, solved using the standard method of [2].The elements of the symmetric matrices of stiffness A and mass B comprise integrals like Eq. ( 4), which are calculated on the elements ∆ q in the domain Ω = Ω h (z) = Q q=1 ∆ q , recalculated into the local coordinates on the master element ∆.

The deviation of the approximate solution
where , h is the maximal size of the finite element ∆ q , p is the order of the FEM scheme, m is the number of the eigenvalue, c 1 and c 2 are coefficients independent of h.

BVP for Helmholtz equation in an equilateral triangle
As an example, we consider the solution of the BVP ( 1 assumed to be an equilateral triangles with the side 4π/3 partitioned into Q = n 2 equilateral triangles ∆ q of sides h = 4π/3n.The eigenvalues of this problem having degenerate spectrum are the integers  1 of paper [1], conserving the continuity of the first and the second derivative of the approximate solution, respectively.As seen from Fig. 1, the errors of the eigenvalue ∆E h 4 (z) of the FEM schemes of the same order p = κ max (p + 1) − 1 are nearly similar and correspond to the theoretical estimates (5), but in the FEM schemes conserving the continuity of the first and the second derivatives of the approximate solution, the matrices of smaller dimension are used that correspond to the length of the vector N smaller by 1.5-2 times than in the schemes with LIP that conserve only the continuity of the functions themselves at the boundaries of the finite elements.Computations were carried out using a dual processor desktop 2 x Xeon 3.2 GHz, 4 GB RAM, Intel Fortran 77 with quadruple precision real*16.The computing time was 3 minutes.

BVP for Helmholtz equation in a d-dimensional hypercube
For benchmark calculations we use the BVP for the HEQ with the BC (II) in a d-dimensional hypercube with the edge length π.Since the variables are separated, the eigenvalues d.Assertion (see also [3]).The hypercube is divided into d! equal simplexes.The vertices of each simplex are located on broken lines, composed of d mutually perpendicular edges, and the extreme vertices of all polygons are located on one of the diagonals of the hypercube (for d = 3 see Fig. 2a).
Algorithm.Input.A single d-dimensional hypercube with vertices the coordinates of which are either 0 or 1 in the Euclidean space R d .The chosen diagonal of the hypercube connects the vertices with the coordinates (0, . . ., 0) and (1, . . ., 1). Output.
) the coordinates of i-th simplex.Local.The coordinates of the vertices of the polygonal line are z k = (z k,1 , . . ., z k,d ), k = 0, . . ., d. 1.For all i = (i 1 , . . ., i d ), the permutations of the numbers (1, . . ., d): 3D HEQ for the cube.In Fig. 2b we show the error of eigenvalues of the 3D BVP for the HEQ at d = 3 with the BC (II) using the FEM scheme with 3D LIP of the order p = p =6.As can be seen from Fig. 2b, the errors of the eigenvalues lie on parallel lines in the double logarithmic scale which corresponds to the theoretical error estimates (5) for the eigenvalues depending on the maximal size of the finite element.In the case of the cube with the edge π divided into 4 3 cubes, each of them comprised by 6 tetrahedrons, the matrices A and B had the dimension 15625 × 15625.The matrices A and B were calculated analytically using Maple 2015, 2x 8-core Xeon E5-2667 v2 3.3 GHz, 512 GB RAM, GPU Tesla 2075, and the AEP was solved during 20 minutes using Intel Fortran.
6D HEQ for the hypercube.We solved HEQ at d = 6 with the BC (II) using FEM scheme with 6D LIP of the order p = p =3.The 6D hypercube having the edge π was divided into n = d! = 6! = 720 simplexes (the size of the finite element being equal to π).On each of them N 1 (p) = (p + d)!/(d!p!) = 84 third-order LIPs were used.The matrices A and B had the dimension 4096 × 4096.The lower part of the spectrum E m is shown in Table 1.The errors of the second, the third and the fourth degenerate eigenvalues are equal to 0.0003, 0.05 and 0.15, respectively.Note, that applying the thirdorder scheme for solving the BVPs of smaller dimension d, we obtained errors of the same order.The calculation time was 9234.46 seconds using Maple 2015, 2x 8-core Xeon E5-2667 v2 3.3 GHz, 512 GB RAM, GPU Tesla 2075.

Conclusion
We demonstrated the efficiency of the schemes using the high-order accuracy LIP and HIP for benchmark calculations of the exactly solvable Helmholtz problems for the triangle and hypercube.The detailed description of the FEM schemes with multivariate LIP and HIP, using a more economical numerical calculation of FEM integrals by means of the cubature formulas, together with benchmark calculations of the parametric BVP for the two-center Coulomb problem, that provide accuracy of the order 10 −12 for eigenvalues of a low part of the spectrum, will be presented elsewhere.

Figure 1 .
Figure 1.The error ∆E 4 of the eigenvalue E h 4 versus the number of elements n and the length of the vector N

Figure 2 .
Figure 2. a) The division of the 3D cube into 3! = 6 equal tetrahedrons (T1,. . .,T6).b) The error ∆E m = E h m − E m calculated by FEM using sixth-order LIPs versus the exact eigenvalue E m .Squares: the cube having the edge π divided into 6 tetrahedrons.Circles: the cube with having the edge π divided into 2 3 cubes, each comprised of 6 tetrahedrons.Solid circles: the cube with the edge π divided into 4 3 cubes, each comprised of 6 tetrahedrons.

Table 1 .
The lower part of the exact spectrum E m and the calculated spectrum E h m for the 6D hypercube.For all k = 0, . . ., d and s = 1, . . ., d: z