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 isoparametric or subparametric 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 Cowper^{4}. A detailed study of symmetric quadrature rules by formulating the problem in polar coordinates is made by Lyness and Jesperson^{5}. 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 future^{613}. Lague and Baldur^{14 }proposed the product formulae based only on the sampling points and weight coefficients of Gauss Legendre quadrature rules. According to Lague and Baldur^{14}, 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 Baldur^{14} 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 study^{1521}. Principal drawback in the symmetric quadrature scheme of Wandzura and Xiao^{22} 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 twodimensional 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:
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 2square in {(ξ, η) 1 < ξ, η < 1} space through isoparametric transformation and that results the equivalent integral I in Eq. 2:
where:

Fig. 1:  Arbitrary triangle divided into quadrilaterals 
and:
Now righthand 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 interelement continuity requirements. Therefore, integration over triangular domains is usually carried out in normalized coordinates.
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 N_{1}(s, t), N_{2}(s, t) and N_{3}(s, t):
where:
The original and the transformed triangles are shown in Fig. 2. Using the shape functions in Eq. 4, we obtain:

Fig. 2(ab): 
Original and transformed triangle, (a) Triangle in (x, y) space and (b) Transformed triangle in (s, t) space 
and hence:
Finally, Eq. 1 reduces to:
Using mathematical transformation equations:
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:
Now, using Eq. 8 and 9, we get:
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:
The integral I in Eq. 1 becomes:
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) 

where, (s_{r},t_{r}) are the new Gaussian points, G_{r }is the corresponding weights for triangles. For Eq. 12, we have:
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 straightline 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, m1 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:
where, (p_{r},q_{r}) are the new Gaussian points, L_{r }is the corresponding weights for triangles. For Eq. 12, we can write:
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., 54points 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 Xiao^{22} for triangle were considered:
Table 3:  Calculated values of the integrals I_{1}, I_{2}, I_{3} and I_{4} 

Computed values are summarized in Table 3. Some important remarks from the Table 3 are:
• 
Usual Gauss quadrature (GQT) for triangles e.g., 7point and 13point rules or the quadrature rule^{22 }cannot evaluate the integral of nonpolynomial 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 710 showed the percentage error in calculation of these integrals.
Again, this study considered the following integrals of rational functions due to Rathod and Karim^{23} to test the influences of formulae.

Fig. 7:  Percentage error in I_{1} 

Fig. 8: 
Percentage error in I_{2} 
Table 4:  Computed results of example 1 

Table 5:  Computed results of example 2 


Fig. 9:  Percentage error in I_{3} 

Fig. 10:  Percentage error in I_{4} 
Table 6:  Computed results of example 3 and 4 

Table 7: 
Absolute error in 

These integrals arise in axisymmetric 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 46. 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 46. 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 (example3, 4) the method described in section4 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 rule^{22} 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 Rule^{24 }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 x^{i}y^{j}, where, i , j are nonnegative 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 results^{19 }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 npoint Gauss quadrature formula can evaluate exactly the integral of polynomial of order up to 2n1. The total error in npoint Gauss quadrature formula to evaluate the integral of polynomial of high order is given by:
where, f^{2n}(x) is the 2nth 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 npoint Gaussian quadrature rule gives exact results for polynomials of degree at most 2n1. 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 x^{6}, x^{5 }y, x^{4} y^{2}, x^{3 }y^{3}, x^{2} y^{4}, x y^{5}, y^{6}. The ppoint rule gives exact result for polynomials of degree up to 2p1 and qpoint rule gives exact result for polynomials of degree up to 2q1.
Table 8:  Comparison of absolute error calculated in I_{E} 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 x^{2p1}y^{2q1}, 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 Npoint rule can calculate all the polynomials of order up to 2p1 or 2q1 in x and y
separately. Hence, N×N point rule can evaluate all the monomials of degree at most 2(p+q)2ie^{2N2}.
The proposed methods are tested on the integral of all monomials x^{i}y^{j} where, i, j are nonnegative 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 studies^{22,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 1114 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 onedimensional 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 axisline 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.