Abstract:
This article proposed symmetrical Gaussian quadrature formulae for triangular domain integrals. As a result, it presents n×n points (for n>1) and
INTRODUCTION
For numerical solution of continuum mechanics problem using Finite Element Method (FEM), an extra and important consideration is involved in the stiffness matrix calculation which involves numerical integration over the corresponding domain. Among all the numerical rules, Gaussian quadrature rule occupies a central role for such calculations. Complication arise from two main sources, firstly the large number of integrations that need to be performed and secondly in methods which use iso-parametric or sub-parametric elements, the presence of the determinant of Jacobean matrix in the denominator of the stiffness matrix for which the integrands are rational functions. Most of the domain integrals encountered in several areas of science and engineering are not amenable to evaluate analytically or tedious in calculation. Such integrals encountered for employing linear elements in the discretization are simple and may be evaluated analytically but, large numbers of integrals are needed to be evaluated. However, encountered integrals for employing higher order elements or for some distorted elements are too complicated. It is highly expected that the expressions for the exact values of the integrals must be evaluated with care and hence, the numerical integration techniques are the best choice. Finite Element Method (FEM) got importance due to the most obvious reason that it can provide solutions to many complicated problems that would be intractable by other numerical methods. The crucial problem of integrating arbitrary functions of two variables over the surface of the triangle were first described by Hammer et al.1 and Hammer and Stroud 2,3.
A table of Gaussian quadrature formulae with symmetrically placed integration points is provided by Cowper4. A detailed study of symmetric quadrature rules by formulating the problem in polar co-ordinates is made by Lyness and Jesperson5. Different researchers described some numerical integration formulae for triangles with precision limited up to 10° and it is not likely that the techniques can be extended much further to give a greater accuracy which may be demanded in future6-13. Lague and Baldur14 proposed the product formulae based only on the sampling points and weight coefficients of Gauss Legendre quadrature rules. According to Lague and Baldur14, one can obtain numerical integration rules of very high degree of precision as the derivation rely on standard Gauss Legendre quadrature rules. However, Lague and Baldur14 have not worked out on explicit weight coefficients and sampling points for application to integrals over a triangular surface. Different reports provided the information about the quadrature for triangle in a systematic manner in their study15-21. Principal drawback in the symmetric quadrature scheme of Wandzura and Xiao22 was that, one must manually adjust the annealing parameters several times, before the process yields a satisfactory initial approximation of weights and abscissae, also, it provides only 6 types of quadrature rules of order up to 30 over triangles.
The versatility of the popular triangular elements can be further enhanced by improved numerical integration schemes and hence, evaluation of the triangular domain integrals with desired accuracy by other technique is preferable. It is notable that the high order Gaussian quadrature formulae available only for the square domain integrals and the same are demanded for the triangular domain integrals. But, the derivation of the higher order Gaussian quadrature for triangular domain integrals is not so easy and indeed very difficult task.
The aim of this article was to present symmetrical extended Gaussian quadrature formulae avoiding the crowd of Gaussian integration points and weights in the calculation process.
MATERIALS AND METHODS
The study of this article was carried out at Mathematics Department Lab, Huston Tillotson University, Texas, USA, from August, 2016 to July, 2019. The results are tested and verified at Mathematics Department Lab, Shahjalal University of Science and Technology, Bangladesh. First, it proposed a numerical integration scheme to evaluate the triangular domain integral employing Gaussian quadrature schemes for square domain integrals. This scheme is used as a tool for testing the accuracy for the derived numerical integration formulae for triangular domain integrals. Secondly, it presented 2 types of extended quadrature n×n points (for n>1) and
An error formula is described to calculate the error in two-dimensional domain integral. The error calculated by the new error formulae and the error in the resultant integrals of the proposed methods are found in good agreement. Therefore, a proper balance between accuracy and efficiency is ensured for the presented quadrature schemes.
Consider the triangular domain integral:
| (1) |
The aim of this article is to derive a suitable, highly accurate, efficient method to evaluate the integral I.
INTEGRATION OVER ARBITRARY TRIANGLE (IOAT)
Integration over any triangle can be calculated as a sum of integrals evaluated over three quadrilaterals (Fig. 1). Each quadrilateral in Fig. 1 is transformed into 2-square in {(ξ, η) |-1 < ξ, η < 1} space through iso-parametric transformation and that results the equivalent integral I in Eq. 2:
| (2) |
where:
| Fig. 1: | Arbitrary triangle divided into quadrilaterals |
and:
Now right-hand side of Eq. 2 can be evaluated by use of available higher order Gaussian quadrature for square domain integrals.
TRANSFORMATION OF TRIANGULAR DOMAIN
The simple shapes of the elements restrict severely their applications in the analysis of practical problems, where often quite complex geometrical boundaries needed to be modelled. In FEM solution process, this restriction can be removed by mapping a simple element in the local coordinates into a more complex shape in the global coordinate system. Once a particular form of mapping is adopted and the coordinates are chosen for every element so that these map into contiguous space, then shape functions written in the local element space can be used to represent the function variation over the element in the global space without upsetting the inter-element continuity requirements. Therefore, integration over triangular domains is usually carried out in normalized co-ordinates.
To perform the integration, first map one vertex (vertex 1) to the (1,0), the second vertex (vertex 2) to point (-1,1) and the third vertex (vertex 3) to point (-1, -1)) (Fig. 2). This geometrical transformation is most easily accomplished by use of shape functions N1(s, t), N2(s, t) and N3(s, t):
| (3) |
where:
| (4) |
The original and the transformed triangles are shown in Fig. 2. Using the shape functions in Eq. 4, we obtain:
| Fig. 2(a-b): | Original and transformed triangle, (a) Triangle in (x, y) space and (b) Transformed triangle in (s, t) space |
| (5) |
and hence:
| (6) |
Finally, Eq. 1 reduces to:
| (7) |
Using mathematical transformation equations:
| (8) |
Then integral I of Eq. 7 is transformed into an integral over the surface of the standard square {(ξ, η)| -1< ξ, η <1}.
Now, the determinant of the Jacobean and the differential area are:
and:
| (9) |
Now, using Eq. 8 and 9, we get:
| (10) |
Thus, the triangular domain integral in Eq. 1 is finally converted to square domain integral which can be evaluated by using available higher order Gaussian quadrature for square domain integrals.
Again, if we use the shape functions:
then using mathematical transformations, the original triangle can be transformed to a triangle with vertices 1(1, -1), 2(0, 1), 3(-1, -1). Now consider the transformation:
| (11) |
The integral I in Eq. 1 becomes:
| (12) |
SYMMETRIC GAUSS QUADRATURE FOR TRIANGLE (SGQTS)
The Gauss points are calculated simply for i = 1, m and j = 1, n. Thus, the m×n points Gaussian quadrature formula for Eq. 10 is:
| Table 1: | Computed Gauss points (s, t) and corresponding weights G for m×m point formula (SGQTS) |
| (13) |
where, (sr,tr) are the new Gaussian points, Gr is the corresponding weights for triangles. For Eq. 12, we have:
| (14) |
Computed Gauss points and weights for different values of n are listed in Table 1 and Fig. 3 and 4 showed the distribution of Gaussian points for n = 10. The Gauss points in Fig. 4 can be obtained by simply interchanging s and t shown in Fig. 3 and the corresponding weights are the same. Also, the distribution of Gauss points is symmetric about the straight-line y = 0 or x = 0, which can significantly minimize the computational effort to calculate the quadrature points and weights. But, in both cases it is seen that there are crowding of gauss points at one side within the triangle and that is one of the major causes of error germane in the calculation. To get rid of these crowding, further modification is obtained in the next section.
SYMMETRIC GAUSS QUADRATURE FOR TRIANGLE (SGQTM)
It is clearly noticed in the Eq. 13 and 14 that for each i (i = 1, 2, 3, ...., m), j varies from 1 to n and hence, at the terminal value i = m there are n crowding points as shown in Table 1 and Fig. 3 and 4.
| Fig. 3: | Distribution of Gauss points (s, t) for n = 10 SGQTS (100 points) |
| Fig. 4: | Distribution of Gauss points (s’, t’) for n = 10 SGQTS (100 points) |
| Fig. 5: | Distribution of Gaussian points (p, q) for n = 10 SGQTM (54 points) |
| Fig. 6: | Distribution of Gaussian points (p', q') for n = 10 SGQTM (54 points) |
| Table 2: | Computed Gauss points (p,q) and Corresponding weights L for |
To overcome the crowding of points, considering the geometry of the elements and using algebraic manipulation, we are taking j dependent on i for the calculation of quadrature points and corresponding weights. Gauss points and weights are calculated for i = 1, m-1 and j = 1, n + 1 – i, where, n>m, that is for n = m; total
| (15) |
where, (pr,qr) are the new Gaussian points, Lr is the corresponding weights for triangles. For Eq. 12, we can write:
| (16) |
The
APPLICATION EXAMPLES
To show the accuracy and efficiency of the derived formulae some examples with known results are considered. To compare the results, the results obtained by using available for Gauss 7×7 points and 13×13 points methods and the quadrature rule of Wandzura and Xiao22 for triangle were considered:
| Table 3: | Calculated values of the integrals I1, I2, I3 and I4 |
Computed values are summarized in Table 3. Some important remarks from the Table 3 are:
| • | Usual Gauss quadrature (GQT) for triangles e.g., 7-point and 13-point rules or the quadrature rule22 cannot evaluate the integral of non-polynomial functions accurately |
| • | Splitting any triangle into quadrilaterals (IOAT) provides the way of using Gaussian quadrature for square and the convergence rate is slow, but satisfactory in view of accuracy |
| • | The new Gaussian quadrature formula for triangle (SGQTS and SGQTM) are exact in view of accuracy, efficiency and rate of convergence is high also SGQTM used less computational effort |
Figure 7-10 showed the percentage error in calculation of these integrals.
Again, this study considered the following integrals of rational functions due to Rathod and Karim23 to test the influences of formulae.
| Fig. 7: | Percentage error in I1 |
| Fig. 8: | Percentage error in I2 |
| Table 4: | Computed results of example 1 |
| Table 5: | Computed results of example 2 |
| Fig. 9: | Percentage error in I3 |
| Fig. 10: | Percentage error in I4 |
| Table 6: | Computed results of example 3 and 4 |
| Table 7: | Absolute error in |
These integrals arise in axi-symmetric finite element method with linear triangular element as well as in finite element formulations of second order linear differential equations by use of triangular element with two straight sides and one curved side. Consider:
Results are summarized in Table 4-6. These data strongly substantiated the influences of numerical evaluation of the integrals described in this study. Some important comments may be drawn from the Table 4-6. In these tables formula 1 stands for Eq. 13 and 15, whereas, formula 2 stands for Eq. 14 and 16:
| • | For the integrand |
| • | Similarly, for the integrand |
| • | Similarly, for the integrand |
| • | Similar influences of these formulae may be observed for different conditions on β, γ |
| • | The existing 7 - point, 13 - point GQT or quadrature rule22 cannot evaluate the rational integrals accurately |
It is evident that the new formulae e.g., SGQTS and SGQTM are accurate in view of accuracy and equally applicable for any geometry that is for different values of β and γ. The Fekete quadrature Rule24 is also tested for all the examples and for all problems the present methods are found the best.
The proposed methods are also tested on the integral of all monomials xiyj, where, i , j are non-negative integers such that i+j<30. Table 7 presented the absolute error over corresponding monomials integrals for each quadrature of order between 1 and 30. The results are compared with the previous study results19 and it is observed that the new method SGQTM is always accurate in view of both accuracy and efficiency and hence a proper balance is observed.
ERROR ANALYSIS
The n-point Gauss quadrature formula can evaluate exactly the integral of polynomial of order up to 2n-1. The total error in n-point Gauss quadrature formula to evaluate the integral of polynomial of high order is given by:
where, f2n(x) is the 2n-th derivative of the function at a point x in the interval [-1, 1]. In this article the triangular domain integral is evaluated by converting it to a square domain integral. Consider the square domain integral:
Integrating first with respect to x keeping y fixed, then integrating with respect to y, using m points along x direction and n points along y direction, the error in this method is found to be of the form:
| Fig. 11: | Polynomial of order 15 |
| Fig. 12: | Polynomial of order 20 |
where, is the 2m th partial derivative of the function with respect to x, is the 2n th partial derivative of the function with respect to y and and are points somewhere in the domain [-1, 1] ×[-1, 1].
The n-point Gaussian quadrature rule gives exact results for polynomials of degree at most 2n-1. Thus, for n = 2 we have a rule with 4 nodes which is exact for any polynomial of degrees at most 3 in x and y separately, so the total degree of this monomial is at most 6. But this rule is not exact for all monomials of degree at most 6, which includes x6, x5 y, x4 y2, x3 y3, x2 y4, x y5, y6. The p-point rule gives exact result for polynomials of degree up to 2p-1 and q-point rule gives exact result for polynomials of degree up to 2q-1.
| Table 8: | Comparison of absolute error calculated in IE using the new error formula and error in SGQTS |
| Fig. 13: | Polynomial of order 25 |
| Fig. 14: | Polynomial of order 30 |
But p×q -point rule cannot give exact result for x2p-1y2q-1, which is a monomial of order 2(p + q)-2. Let, N be the maximum value of p+q for each term in the monomial. Then N-point rule can calculate all the polynomials of order up to 2p-1 or 2q-1 in x and y
separately. Hence, N×N- point rule can evaluate all the monomials of degree at most 2(p+q)-2ie2N-2.
The proposed methods are tested on the integral of all monomials xiyj where, i, j are non-negative integers such that i+j<N. Table 8 presented the absolute error over corresponding monomial integral of order up to 30. Table 8 is in good agreement with the above statement. The results are compared with the results of previous studies22,25,26 and it is observed that the new method SGQTS is accurate in view of both accuracy and efficiency and hence a proper balance is observed. Figure 11-14 showed the absolute error in the integral of polynomial of order 15, 20, 25 and 30.
CONCLUSION
This paper showed first the integral over the triangular domain can be computed as the sum of three integrals over the square domain, then the readily available quadrature formulae for the square can be used for the desired accuracy. Secondly, it presented new techniques to derive quadrature formulae utilizing the one-dimensional Gaussian quadrature formulae and that overcomes all the difficulties pertinent to the higher order formulae. The symmetric distribution of Gauss points derived in these methods can minimize the computational afford of numerical evaluation of integrals in triangular domains. It is observed that the scheme SGQTM is appropriate for the triangular domain integrals as it requires less computational effort for desired accuracy and efficiency.
SIGNIFICANCE STATEMENT
The current work is on a new triangle element type with vertices (0, -1), (0, 1), (1, 0). Gauss points are symmetrically distributed about the axis-line y = 0, reducing the calculation effort significantly. Several techniques are discussed for the domain integrals with their drawbacks, then the drawbacks are also removed by updated findings. An error formula is described for the integration over triangular element. The model is verified by comparing with other existing methods. Analyzed the error on each step for each example types. The findings of this study are accurate, efficient, faster with less computational effort. Thus, they are believed to earn a better place in the field of triangular domain integrals.