CONSTRUCTING STEKLOV-TYPE CUBATURE FORMULAS FOR A FINITE ELEMENT IN THE SHAPE OF A BIPYRAMID

A n z h e l i k a M o t a i l o Corresponding author PhD Department of Natural Sciences Kherson State Maritime Academy Ushakova ave., 20, Kherson, Ukraine, 73009 E-mail: akillehzna@gmail.com G a l i n a T u l u c h e n k o Doctor of Technical Sciences Department of Higher Mathematics National Technical University «Kharkiv Polytechnic Institute» Kyrpychova str., 2, Kharkiv, Ukraine, 61002 This paper reports the construction of cubature formulas for a finite element in the form of a bipyramid, which have a second algebraic order of accuracy. The proposed formulas explicitly take into consideration the parameter of bipyramid deformation, which is important when using irregular grids. The cubature formulas were constructed by applying two schemes for the location of interpolation nodes along the polyhedron axes: symmetrical and asymmetrical. The intervals of change in the elongation (compression) parameter of a bipyramid semi-axis have been determined, within which interpolation nodes of the constructed formulas belong to the integration region, while the weight coefficients are positive, which warrants the stability of calculations based on these cubature formulas. If the deformation parameter of the bipyramid is equal to unity, then both cubature formulas hold for the octahedron and have a third algebraic order of accuracy. The resulting formulas make it possible to find elements of the local stiffness matrix on a finite element in the form of a bipyramid. When calculating with a finite number of digits, a rounding error occurs, which has the same order for each of the two cubature formulas. The intervals of change in the elongation (compression) parameter of the bipyramid semi-axis have been determined, which meet the requirements, which are employed in the ANSYS software package, for deviations in the volume of the bipyramid from the volume of the octahedron. Among the constructed cubature formulas for a bipyramid, the optimal formula in terms of the accuracy of calculations has been chosen, derived from applying a symmetrical scheme of the arrangement of nodes relative to the center of the bipyramid. This formula is invariant in relation to any affinity transformations of the local bipyramid coordinate system. The constructed cubature formulas could be included in libraries of methods for approximate integration used by those software suites that implement the finite element method


Introduction
The finite element method (FEM) is one of the widely used techniques for solving problems related to mathematical physics. The scope of its application covers all boundary-value problems that can be described by differential equations. However, the main disadvantage of a given method is cumbersome calculations that require the use of a large amount of computer memory. Therefore, in line with the process of developing high-performance computers, there is an issue of reducing the volume of FEM calculations.
The search for new opportunities to reduce the time complexity of the FEM algorithm requires a reasonable choice not only of the finite-element grid and node coordinate functions but also of numerical integration formulas. It should also be taken into account that such a choice should ensure the stability and accuracy of the resulting solutions. This task is especially difficult in 3D. If the computational domain has a complex geometric shape, then its discretization only by regular polyhedra is often impossible. In the libraries of computational software that implement FEM, more and more irregular polyhedra appear, which are used in the near-boundary layer when discretization of the region of the problem being solved.
Despite the large enough selection of software implementing FEM, the search for new formulas and methods of numerical integration in the regions of both regular and irregular polyhedra continues. Confirmation of this fact is in works [1][2][3], which address the construction of symmetrical cubature formulas of a high order of accuracy on the tetrahedron, hexahedron, prism, and pyramid. The resulting formulas are accurate for algebraic polynomials of three variables up to the order of 20 inclusive, and allow for the rapid construction of stiffness and mass matrices in finite-element analysis. However, the use of these formulas in the region of an octahedron or a bipyramid, which can be represented by a combination of tetrahedra or prisms, would have an excessively large number of nodes and increase the volume of calculations. That means that the subject of research tackling the numerical integration in the region of the octahedron and the bipyramid is relevant.

Literature review and problem statement
Paper [4] shows that an octahedron in the ensemble with a tetrahedron forms a hybrid grid, which, in terms of the accuracy of calculations, is not inferior to the tetrahedral one, and is optimal as regards the time of calculations.
To implement finite-element calculations, the authors of the cited paper built the piecewise-linear node coordinate functions corresponding to the vertices of the octahedron and its center. Computational experiments were performed for the thermal conductivity problem. Numerical integration in the octahedron region is based on its representation in the form of a combination of eight tetrahedra. In this case, a cubature formula, which is known from work [5], is used for the tetrahedron. The issue of constructing a cubature formula for the octahedron was not considered in [4]. The reason for this may relate to that the task of finding new effective methods of discretization the three-dimensional region was resolved only for the visualization of objects, which makes the relevant studies of numerical analysis impractical. At the same time, the authors of [4] note that the grids of the tetrahedral-octahedral structure can be used to improve the algorithm of parallel FEM calculations. In the case where an extremely large number of finite elements (FE) are needed to obtain a solution, parallel processing is useful at each FEM stage, including grid refinement.
Paper [6] shows that tetrahedral-octahedral grids have two advantages compared to tetrahedral ones. First, the grids containing octahedra are symmetrical, which simplifies the implementation of the algorithm for refining such a grid. Second, the grids of the tetrahedral-octahedral structure generate fewer comparison classes than the simplest tetrahedral grids, which is important in multilevel FEM calculations. The refinement algorithm [5,6] generates two comparison classes in 3D instead of three, which is important in terms of practical application.
The effectiveness of using the grid of a tetrahedral-octahedral structure was confirmed in work [7] when solving practical problems related to electromagnetism. In the cited work, the methods of dynamic load balancing are studied, in the parallel refinement of the grid in order to optimize the architecture of a parallel system of calculations and program methods. With a parallel design of the FEM algorithm, load balancing should reach the minimum execution time, distributing the task evenly among processors. The authors of work [7] devised a new approach to modeling, analysis, and design of communication circuits with the parallel improvement of the FE grid using Petri networks. A given approach has been introduced to an effective concurrent data processing pipeline strategy to improve the FEM tetrahedral grid, which employs octahedra at two stages of refinement [8]. The issue of numerical analysis was not considered in [6][7][8] since those works addressed the optimization of system resources and the balance of overall computational loads.
Work [9] investigated the possibilities of using grids of a tetrahedral-octahedral structure in solving boundary problems of mathematical physics for elliptical type equations. It has been proved that the piecewise-linear approximation of the temperature field, as well as the piecewise-linear and quadratic finite-element approximations of functions of the displacement field on the grids of the tetrahedral-octahedral structure converge to the exact solution in problems with the differential operator of the second-order of the elliptical type. Verification of the obtained theoretical conclusions was carried out in the Maple computer mathematics system for problems related to thermal conductivity and linear elasticity theory in the simplest regions that have the shape of a parallelepiped, which makes it possible to compare the accuracy of numerical solutions with those derived analytically. To find the elements of the local stiffness matrix of an octa-hedron with piecewise-linear functions, the cited work's authors built a formula for numerical integration, which is a consequence of the application of the properties of the triple integral in the region of the octahedron. To find elements of the local stiffness matrix of an octahedron with quadratic functions, the procedures for integrating by analytical methods that are built into the Maple processor were applied. Formulas of numerical integration for FE in the form of an octahedron to the seventh algebraic order of accuracy were constructed in work [10]. These formulas reduce the time complexity of the FEM algorithm when using cells in the form of an octahedron but cannot be applied to FE of irregular geometric shape.
The first step in overcoming the difficulties associated with the adaptation of a tetrahedral-octahedral grid to regions of complex geometry is work [11], in which a bipyramid was studied as FE.
In this case, the bipyramid is considered, which is obtained from the octahedron by elongating (compressing) one of its half-diagonals. In the cited work, two systems of node coordinate functions were built that correspond to the vertices and center of the polyhedron. Both systems of the node coordinate functions of such a bipyramid depend on a parameter equal to the length of the deformed half-diagonal.
To perform finite-element calculations involving bipyramids, which could be used when discretization three-dimensional regions in the near-border layer, it is necessary to derive formulas for numerical integration on a given polyhedron. That allows us to assert that it is advisable to conduct a study on the construction of cubature formulas on FE in the shape of a bipyramid.

The aim and objectives of the study
The purpose of this work is to study the possibility of constructing cubature formulas of the Steklov type for FE in the shape of a bipyramid to calculate the elements of the stiffness matrix on a given geometric carrier.
To accomplish the aim, the following tasks have been set: -to construct a cubature formula of the second algebraic order of accuracy for a bipyramid with the symmetrical localization of nodes relative to the center of the polyhedron; -to build a cubature formula of the second algebraic order of accuracy for a bipyramid with the asymmetric localization of nodes relative to the center of the polyhedron; -to define the conditions for the use of constructed cubature formulas in the algorithmization of FEM; -to assess the effect of the approximate execution of arithmetic operations in a finite bit grid on the final accuracy of the application of the proposed cubature formulas.

The study materials and methods
During the research, we employed the main provisions of a finite element method, in particular, the properties of the basis functions of finite elements; rules for the construction of local stiffness matrices (thermal conductivity) for the problems of mathematical physics in an elliptical type; the estimates of permissible deviations of the geometric shape of finite elements from regular polyhedra [12].
The construction of new cubature formulas relies on wellknown theorems in the theory of approximate integration methods. Namely, a theorem about the existence of a cubature formula with positive weight coefficients and the property of computational stability of calculations using it; a theorem about possible algebraic accuracy of the cubature formula depending on the localization of interpolation nodes [13].
The research was carried out using the tools of Maple computer mathematics software.

5.
The results of studying the possibilities of constructing cubature formulas of Steklov type and evaluation of their computational characteristics

1. The symmetrical scheme of localization of cubature formula nodes
Let Ω Ω Ω = ∪ 1 2 be the region in the shape of a bipyramid ( Fig. 1), where Ω 1 1 0 x y z x y z p z , , : , and , , p is the K 0 K 5 half-diagonal elongation (compression) coefficient that meets the conditions p>0 and p ∈ R. , , were derived by us earlier and reported in [12]: Let the function f(x,y,z) be continuous in the region Ω.
Consider the integral f x y z x y z , , , where A s are the weight coefficients, a s are the interpolation nodes, Q is the number of nodes.
The following value is considered to be an error in cubature formula (2): which is the residual term of formula (2) in the Lagrange form, where n is the highest order of the algebraic polynomial for which formula (2) is exact.
Let the triple integral in formula (2) ,z) is Ω, that is: where a ijk are the coefficients, α = α(i,j,k) is the multi-index, |α| = i+j+k. Given the symmetry in the region Ω, the nodes a s in interpolation formula (2) are to be chosen separately along the bipyramid's half-diagonals at a distance t from the center K 0. Then the nodes of the cubature formula are the vertices of the octahedron x y z t + + ≤ . Integrate polynomial (3): Since (2) must be exact for the second-power polynomials, we derive the following system of equations: .
The following is the solution to the system of equations (5) Thus, the corresponding cubature formula takes the form: where the coefficients A s s = ( ) 1 4 , that are equal correspond to the nodes a s s = ( ) 1 4 , , located in the Оху plane, where the integration region Ω is centrally symmetrical [13].
Let us analyze a given solution graphically. Fig. 2 shows the dependence of the obtained values of the coordinate t of the nodes and weight coefficients in formula (2) on the parameter p of the bipyramid. Obviously, the number of nodes in (7) can be reduced to Q = 5 by choosing p ≈ 0.42.

2. The asymmetrical scheme of localization of cubature formula nodes
Let us arrange the nodes a s in cubature formula (2) one at a time along the half-diagonals of the bipyramid Ω (Fig. 1) in such a way that the following ratio is satisfied: where |a s | is the distance from a node in formula (2) to the bipyramid center K 0 , |K 0 K i | is the distance from the top of the bipyramid to the center K 0 . Then the nodes in (2) are at the vertices of the bipyramid: where ω 1 0 x y z x y z p t z , , : , and ω 2 0 x y z x y z t z , , : , .
The requirement that the desired cubature formula should be exact for polynomials (3), and the integration of polynomials (4), leads to the following system of equations: .

A p A A t A A t A p A t p
The following is the solution to the system of equations (8) Thus, we obtain another cubature formula in form (7) with six as nodes a s and the weight coefficients A s , to be determined from formulas (9).
Our analysis of the graphical parameters of formulas (9) for cubature formula (7) indicates that the specified cubature formula has five interpolation nodes when p ≈ 0.42 (Fig. 3). In this case, the value of p, at which A 5 = 0, is the same in formulas (6) and (9)

3. Determining the conditions for using constructed cubature formulas in the algorithmization of FEM
Studying the regions of values for functions a s = a s (p) and A s = A s (p) shows that the nodes a s that satisfy formulas (6) belong to the integration region Ω at p ≥ 0.52.
The nodes a s that satisfy (8) belong to the integration Ω at p>0. In this case, the node a s , which has coordinates (0,0,pt), also belongs to the half-diagonal K 0 K 5 of the bipyramid (Fig. 1). Indeed, the condition pt ≤ p, where p, t >0, is met at t ≤ 1, which is confirmed by Fig. 3.
The weight coefficients A s (s ≠ 5), which satisfy formulas (6) and (9), are positive at p>0, and the coefficient A 5 ≥ 0 at p>0.42 regardless of the choice of formulas (6) or (9).
The condition under which the nodes of the cubature formula belong to the integration region is optional [13], unlike the condition of the positiveness of the weight coefficients, which warrants the stability of calculations based on the cubature formula [14].
Let us evaluate the resulting cubature formulas according to the Skewness asymmetry indicators used by the ANSYS software package (USA) [15,16]. High and permissible indicators of the quality of calculations are considered such at which the value of the FE volume of an irregular geometric shape deviates from the volume of the regular polyhedron by no more than 10 % and 25 %. For the bipyramid, the values 0.66 ≤ p ≤ 0.86 and 0.51 ≤ p ≤ 1.1 correspond to the high and sufficient indicator.
Thus, the cubature formulas constructed according to (7) satisfy the conditions of stability and precision of calculations at 0.51 ≤ p ≤ 1.1. Within a given interval, the resulting formulas have six interpolation nodes, which can be determined from (6) or (9). In addition, note that at p = 1, the systems of equations (5) and (8) have a common solution, which, when substituting into formula (7), is a cubature formula for the octahedron [10], which is accurate for the third-power algebraic polynomials.
The constructed cubature formulas hold when calculating the elements of the FE stiffness matrix in the shape of a bipyramid with seven nodes located at the vertices of the polyhedron K i i = ( ) 1 6 , and its center -point K 0 (Fig. 1). The corresponding node coordinate functions of the bipyramid are as follows [11]: where p>0. At the same time, the choice of FE with one or another set of node coordinate functions depends, first of all, on the type of problem being solved, the geometry of the calculated region, and the type of boundary conditions. However numerical integration formulas are part of the FEM algorithm, so it is important to get estimates of the accuracy of calculations according to these formulas.

4. Estimating the accuracy of resulting cubature formulas in a finite bit grid
When calculating the bipyramid stiffness matrix, the accuracy of the obtained formulas at different values for the parameter p is high but depends on the rounding errors that occur when performing arithmetic operations. The results of    The calculation was performed over a segment of 0.51 ≤ p ≤ 1.1, which corresponds to the studied changes in the p parameter. According to the results from calculating an error in (7) for the two options for the selection of nodes and weights, the optimal accuracy of the calculations is produced by the cubature formula, which was obtained by applying a symmetrical scheme of the nodes relative to the center of the bipyramid. In addition, this formula is invariant in relation to any affinity transformations of the local bipyramid coordinate system, which simplifies its application in the algorithmization of FEM when bypassing polyhedron nodes.

Discussion of results of studying the conditions for applying the obtained cubature formulas
Our results are based on meeting the conditions in the theorems about the existence of a cubature formula with positive weight coefficients and the algebraic accuracy of the cubature formula, depending on the localization of interpolation nodes.
The proposed cubature formulas (7) with nodes and coefficients determined from (6) or (9) have a second algebraic accuracy order. However, these cubature formulas are used in calculations with a finite number of digits, so when performing arithmetic operations, a rounding error occurs. Indeed, the results of the calculations in Fig. 4 clearly illustrate that the accuracy of the values received is not less than 7•10 -20 for FE in the form of a bipyramid with six interpolation nodes. In addition, the accuracy of calculations is not less than 2•10 -19 for a bipyramid with seven interpolation nodes. In both cases, the accuracy of calculations meets both the high and sufficient requirements put forward by the ANSYS software package. That is why the resulting cubature formula (7) with nodes and weight coefficients that satisfy formulas (6) can be considered optimal in terms of the accuracy of calculations, to be used in finite-element calculations.
The constructed formulas of numerical integration in the region of a bipyramid are endowed with the advantage that they could be used to calculate local stiffness matrices for FE in the form of a bipyramid at different values of the coefficient of elongation (compression) p.
This study was carried out with certain restrictions on the parameter of bipyramid deformation. The resulting cubature formulas can be applied at the values of the deformation parameter 0.66 ≤ p ≤ 0.86 at the high quality indicators of finite-element calculations employed by the ANSYS software package, and 0.51 ≤ p ≤ 1.1 -at the sufficient ones.
The disadvantages of the current study include the fact that only one bipyramid semi-axis is taken into consideration. At the same time, bipyramids that have two or three variable parameters are of practical interest. Such polyhedra better adapt to the boundaries of regions of complex geometric shape.
It should be noted that the proposed cubature formulas cannot be used to find a local mass matrix of bipyramid, the elements of which are the integrals from the three-dimensional algebraic polynomials of the fourth power. This fact points the direction of further research aimed at obtaining a cubature formula in the region of a bipyramid, which would make it possible to find exactly the elements of the matrix of masses on a given polyhedron. That could make it possible to obtain solutions to the dynamic boundary problems by using the FEM, provided that the estimation region is discretized with an irregular grid with cells in the shape of a bipyramid. (7) was constructed for a bipyramid with the symmetrical localization of nodes relative to the center of the polyhedron that satisfies (6). The resulting formula has a second algebraic accuracy order and makes it possible to find exactly the elements of the local bipyramid stiffness matrix. The intervals of change in the bipyramid deformation parameter have been determined, for which interpolation nodes belong to the integration region while the weight coefficients are positive, which warrants the stability of numerical calculations based on a given formula.

A cubature formula
2. A cubature formula (7) was constructed for a bipyramid with the asymmetrical localization of nodes relative to the center of the polyhedron that satisfies (9). The resulting formula has characteristics that are similar to the characteristics of the preceding formula. If the bipyramid deformation parameter p = 1, then the constructed formulas have the same nodes and weight coefficients and are also accurate for the third-order algebraic polynomials.
3. When algorithmizing FEM, the resulting cubature formulas can be used with certain limitations. According to the requirements by the ANSYS software complex for bipyramid volume deviations from the volume of the octahedron, a high level of calculation accuracy is achieved at parameter values 0.66 ≤ p ≤ 0.86, and a sufficient one -at 0.51 ≤ p ≤ 1.1.
4. The constructed cubature formulas are used in calculations with a finite number of digits, so when performing arithmetic operations, a rounding error occurs. At the same time, the accuracy of calculations for a bipyramid with six interpolation nodes is not less than 7•10 -20 , and, for a bipyramid with seven interpolation nodes, not less than 2•10 -19 .