Research Article

# A Symmetric Extended Gaussian Quadrature Formula for Evaluation of Triangular Domain Integrals Farzana Hussain and M.S. Karim

ABSTRACT

This article proposed symmetrical Gaussian quadrature formulae for triangular domain integrals. As a result, it presents n×n points (for n>1) and points (for n>2) quadrature formulae in which the second one is totally free of crowding of Gaussian quadrature points and weights. By suitable transformation of a triangle in global space into its contiguous space, Gauss points and weights are computed which are symmetric about the line of symmetry. For clarity and reference, Gaussian integration points and weights for different values of n are presented in tabular form. The efficiency and accuracy of the schemes are tested through application examples. Finally, an error formula also presented to evaluate the error in monomial/polynomial integration using m×n points method successfully. The error calculated by the new error formula and the error in calculation of integrals by the proposed methods are found in good agreement.

 Services Related Articles in ASCI Similar Articles in this Journal Search in Google Scholar View Citation Report Citation Science Alert

 How to cite this article: Farzana Hussain and M.S. Karim, 2020. A Symmetric Extended Gaussian Quadrature Formula for Evaluation of Triangular Domain Integrals. Journal of Applied Sciences, 20: 145-158. DOI: 10.3923/jas.2020.145.158 URL: https://scialert.net/abstract/?doi=jas.2020.145.158

Received: December 26, 2019; Accepted: February 04, 2020; Published: March 15, 2020

Copyright: © 2020. This is an open access article distributed under the terms of the creative commons attribution License, which permits unrestricted use, distribution and reproduction in any medium, provided the original author and source are credited.

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 points (for n>2) formulae for which Gauss points are symmetrical about the line of symmetry. It is easy to observe that n×n points formulae give rise to huge crowding of Gauss points, but points formulae are totally free of such crowding. Through application examples it is demonstrated that the formulae so presented are accurate in view of accuracy and the point formula is faster as it utilizes minimum number of Gauss points and weights in the calculation process.

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 point formula (SGQTM) 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 points Gaussian quadrature formulae for Eq. 10 is given by: (15)

where, (pr,qr) are the new Gaussian points, Lr is the corresponding weights for triangles. For Eq. 12, we can write: (16)

The points Gaussian quadrature formulae (SGQTM) is now obtained which are crowding free and calculates points instead of m×m points. For clarity and reference, computed Gauss points (p, q) and weights L for different values of m are listed in Table 2. Figure 5 and 6 showed the distribution of Gaussian points for m = 10 i.e., 54-points formula. If we interchange p and q then we obtain (q, p) = (p', q') and G = G'. It is now clear from the Table 2 and Fig. 5 and 6 that the method SGQTM is totally free of crowding of Gauss points and use significantly less number of points and weights.

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 with β ≠ γ = 0 first formulae described in Eq. 13 (SGQTS) and Eq. 15 (SGQTM) are more accurate and rate of convergence is higher, but the formula in Eq. 15 requires very less computational effort • Similarly, for the integrand with γ ≠ β = 0 second formula described in Eq. 14 (SGQTS) and Eq. 16 (SGQTM) are more accurate and convergence is higher. Here, also the formula in Eq. 16 requires very less computational effort • Similarly, for the integrand (example-3, 4) the method described in section-4 are more accurate and convergence is higher, also SGQTM requires very less computational afford • 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.

REFERENCES
1:  Hammer, P.C., O.J. Marlowe and A.H. Stroud, 1956. Numerical integration over simplexes and cones. Math. Tables Other Aids Comput., 10: 130-137.
CrossRef  |  Direct Link  |

2:  Hammer, P.C. and A.H. Stroud, 1956. Numerical integration over simplexes. Math. Tables Other Aids Comput., 10: 137-139.
CrossRef  |  Direct Link  |

3:  Hammer, P.C. and A.H. Stroud, 1958. Numerical evaluation of multiple integrals II. Math. Comput., 12: 272-280.
CrossRef  |  Direct Link  |

4:  Cowper, G.R., 1973. Gaussian quadrature formulas for triangles. Int. J. Numer. Methods Eng., 7: 405-408.
CrossRef  |  Direct Link  |

5:  Lyness, J.N. and D. Jespersen, 1975. Moderate degree symmetric quadrature rules for the triangle. IMA J. Applied Math., 15: 19-32.
CrossRef  |  Direct Link  |

6:  Lannoy, F.G., 1977. Triangular finite elements and numerical integration. Comput. Struct., 7: 613-613.
CrossRef  |  Direct Link  |

7:  Laurie, D.P., 1977. Automatic numerical integration over a triangle. CSIR Spec. Rep. WISK 273, National Institute for Mathematical Science, Pretoria.

8:  Laursen, M.E. and M. Gellert, 1978. Some criteria for numerically integrated matrices and quadrature formulas for triangles. Int. J. Numer. Methods Eng., 12: 67-76.
CrossRef  |  Direct Link  |

9:  Lether, F.G., 1976. Computation of double integrals over a triangle. J. Comput. Applied Math., 2: 219-224.
CrossRef  |  Direct Link  |

10:  Hillion, P., 1977. Numerical integration on a triangle. Int. J. Numer. Methods Eng., 11: 797-815.
CrossRef  |  Direct Link  |

11:  Reddy, C.T., 1978. Improved three point integration schemes for triangular finite elements. Int. J. Numer. Methods Eng., 12: 1890-1896.
CrossRef  |  Direct Link  |

12:  Reddy, C.T. and D.J. Shippy, 1981. Alternative integration formulae for triangular finite elements. Int. J. Numer. Methods Eng., 17: 133-139.
CrossRef  |  Direct Link  |

13:  Dunavant, D.A., 1985. High degree efficient symmetrical Gaussian quadrature rules for the triangle. Int. J. Numer. Methods Eng., 21: 1129-1148.
CrossRef  |  Direct Link  |

14:  Lague, G. and R. Baldur, 1977. Extended numerical integration method for triangular surfaces. Int. J. Numer. Methods Eng., 11: 388-392.
CrossRef  |  Direct Link  |

15:  Lyness, J.N. and R. Cools, 1994. A survey of numerical cubature over triangles. Proc. Symposia Applied Math., 48: 127-150.

16:  Cools, R. and P. Rabinowitz, 1993. Monomial cubature rules since “Stroud”: A compilation. J. Comput. Applied Math., 48: 309-326.
CrossRef  |  Direct Link  |

17:  Cools, R., 2003. An encyclopaedia of cubature formulas. J. Complexity, 19: 445-453.
CrossRef  |  Direct Link  |

18:  Bernardini, F., 1991. Integration of polynomials over n-dimensional polyhedra. Comput.-Aided Design, 23: 51-58.
CrossRef  |  Direct Link  |

19:  Rathod, H.T. and H.S.G. Rao, 1996. Integration of polynomials over an arbitrary tetrahedron in Euclidean three-dimensional space. Comput. Struct., 59: 55-65.
CrossRef  |  Direct Link  |

20:  Rathod, H.T. and H.S.G. Rao, 1998. Integration of trivariate polynomials over linear polyhedra in Euclidean three-dimensional space. ANZIAM J., 39: 355-385.
CrossRef  |  Direct Link  |

21:  Rathod, H.T., K.V. Nagaraja, B. Venkatesudu and N.L. Ramesh, 2004. Gauss legendre quadrature over a triangle. J. Indian Inst. Sci., 84: 183-188.

22:  Wandzurat, S. and H. Xiao, 2003. Symmetric quadrature rules on a triangle. Comput. Math. Applic., 45: 1829-1840.
CrossRef  |  Direct Link  |

23:  Rathod, H.T. and M.S. Karim, 2002. An explicit integration scheme based on recursion for the curved triangular finite elements. Comput. Struct., 80: 43-76.
CrossRef  |  Direct Link  |

24:  Taylor, M.A., B.A. Wingate and R.E. Vincent, 2000. An algorithm for computing Fekete points in the triangle. SIAM J. Numer. Anal., 38: 1707-1720.
CrossRef  |  Direct Link  |

25:  Lewis, P.E. and J.P. Ward, 1991. The Finite Element Method. Addison-Wesley Publishing Company, New York.

26:  Karim, M.S., 2001. Integration of some bivariate polynomials with rational denominators-An application to finite element method. Ph.D. Thesis, Department of Mathematics, Bangalore University, Bangalore, India. ©  2021 Science Alert. All Rights Reserved