On the Identification of Symmetric Quadrature Rules for Finite Element Methods

In this paper we describe a methodology for the identification of symmetric quadrature rules inside of quadrilaterals, triangles, tetrahedra, prisms, pyramids, and hexahedra. The methodology is free from manual intervention and is capable of identifying an ensemble of rules with a given strength and a given number of points. We also present polyquad which is an implementation of our methodology. Using polyquad we proceed to derive a complete set of symmetric rules on the aforementioned domains. All rules possess purely positive weights and have all points inside the domain. Many of the rules appear to be new, and an improvement over those tabulated in the literature.


Introduction
When using the finite element method to solve a system of partial differential equations it is often necessary to evaluate surface and volume integrals inside of a standardised domain Ω [1][2][3]. A popular numerical integration technique is that of Gaussian quadrature in which where f (x) is the function to be integrated, {x i } are a set of N p points, and {ω i } the set of associated weights. The points and weights are said to define a quadrature rule. A rule is said to be of strength φ if it is capable of exactly integrating any polynomial of maximal degree φ over Ω. A degree φ polynomial p(x) with x ∈ Ω can be expressed as a linear combination of basis polynomials where P φ is the set of basis polynomials of degree ≤ φ. From the linearity of integration it therefore follows that a strength φ quadrature rule is one which can exactly integrate the basis. Taking f ∈ P φ the task of obtaining an N p point quadrature rule of strength φ is hence reduced to finding a solution to a system of |P φ | nonlinear equations. This system can be seen to possess (N D + 1)N p degrees of freedom where N D ≥ 2 corresponds to the number of spatial dimensions. In the case of N p 10 the above system can often be solved analytically using a computer algebra package. However, beyond this it is usually necessary to solve the above system-or a simplification thereof-numerically. Much of the research into multidimensional quadrature over the past five decades has been directed towards the development of such numerical methods. The traditional objective when constructing quadrature rules is to obtain a rule of strength φ inside of a domain Ω using the fewest number of points. To this end efficient quadrature rules have been derived for a variety of domains: triangles [4][5][6][7][8][9][10][11][12][13], quadrilaterals [11,14,15], tetrahedra [7,9,16,17], prisms [18], pyramids [19], and hexahedra [11,[20][21][22]. For finite element applications it is desirable that (i) points are arranged symmetrically inside of the domain, (ii) all of the points are strictly inside of the domain, and (iii) all of the weights are positive. The consideration given to these criterion in the literature cited above depends strongly on the intended field of application-not all rules are derived with finite element method in mind.
Much of the existing literature is predicated on the assumption that the integrand sits in the space of P φ . Under this assumption there is little, other than the criteria listed above, to distinguish two N p rules of strength φ; both can be expected to compute the integral exactly with the same number of functional evaluations. It is therefore common practice to terminate the rule discovery process as soon as a rule is found. However, there are cases when either the integrand is inherently non-polynomial in nature, e.g. the quotient of two polynomials, or of an high degree, e.g. a polynomial raised to a high power. In these circumstances the above assumption no longer holds and it is necessary to consider the truncation term associated with each rule. Hence, within this context it is no longer clear that the traditional objective of minimising the number of points required to obtain a rule of given strength is suitable: it is possible that the addition of an extra point will permit the integration of several of the basis functions of degree φ + 1.
Over the past five or so years there has also been an increased interest in numerical schemes where the same set of points are used for both integration and interpolation. One example of such a scheme is the flux reconstruction (FR) approach introduced by Huynh [23]. In the FR approach there is a need for quadrature rules that (i) are symmetric, (ii) remain strictly inside of the domain, (iii) have a prescribed number of points, and (iv) are associated with a well conditioned nodal basis for polynomial interpolation. These last two requirements exclude many of the points tabulated in the literature. Consequently, there is a need for bespoke or designer quadrature rules with non-standard properties.
This paper describes a methodology for the derivation of symmetric quadrature rules inside of a variety of computational domains. The method accepts both the number of points and the desired quadrature strength as free parameters and-if successful-yields an ensemble of rules. Traits, such as the positivity of the weights, can then be assessed and rules binned according to their suitability for various applications. The remainder of this paper is structured as follows. In section 2 we introduce the six reference domains and enumerate their symmetries. Our methodology is presented in section 3. Based on the approach of Witherden and Vincent [12] the methodology requires no manual intervention and avoids issues relating to ill-conditioning. In section 4 we proceed to describe our open-source implementation, polyquad. Using polyquad a variety of truncation-optimised rules, many of which appear to improve over those tabulated in the literature, are obtained and presented in section 5. Finally, conclusions are drawn in section 6.

Basis polynomials
The defining property of a quadrature rule for a domain Ω is its ability to exactly integrate the set of basis polynomials, P φ . This set has an infinite number of representations the simplest of which being the monomials. In two dimensions we can express the monomials as where φ is the maximal degree. Unfortunately, at higher degrees the monomials become extremely sensitive to small perturbations in the inputs. This gives rise to polynomial systems which are poorly conditioned and hence difficult to solve numerically [9,16]. A solution to this is to switch to an orthonormal basis set defined in two dimensions as where δ iµ is the Kronecker delta. In addition to being exceptionally well conditioned orthonormal polynomial bases have other useful properties. Taking the constant mode of the basis to be ψ 00 (x) = 1/c we see that from which we conclude that all non-constant modes of the basis integrate up to zero. Following Witherden and Vincent [12] we will use this property to define the truncation error associated with an N p point rule This definition is convenient as it is free from both integrals and normalisation factors. The task of constructing an N p point quadrature rule of strength φ is synonymous with finding a set of points and weights that minimise ξ(φ).
Although the above discussion has been presented primarily in two dimensions all of the ideas and relations carry over into three dimensions.

Symmetry orbits
A symmetric arrangement of N p points inside of a reference domain can be decomposed into a linear combination of symmetry orbits. This concept is best elucidated with an example. Consider a line segment defined by [−1, 1]. The segment possesses two symmetries: an identity transformation and a reflection about the origin. For an arrangement of distinct points to be symmetric it follows that if there is a point at α where 0 < α ≤ 1 there must also be a point at −α. We can codify this by writing S 2 (α) = ±α with |S 2 | = 2. The function S 2 is an example of a symmetry orbit that takes a single orbital parameter, α, and generates two distinct points. In the limit of α → 0 the two points become degenerate. We handle this degeneracy by introducing a second orbit, S 1 = 0, with |S 1 | = 1. Having identified the symmetries we may now decompose a symmetric arrangement of points as where n 1 ∈ {0, 1} and n 2 ≥ 0 with the constraint on n 1 being necessary to ensure uniqueness. This is a constrained linear Diophantine equation; albeit one that is trivially solvable and admits only a single solution. As a concrete example we take N p = 11. Solving the above equation we find n 1 = 1 and n 2 = 5. The n 1 orbit does not take any arguments and so does not contribute any degrees of freedom. Each n 2 orbit takes a single parameter, α, and so contributes one degree of freedom for a grand total of five. This is less than half that associated with the asymmetrical case. Hence, by parameterising the problem in terms of symmetry orbits it is possible to simultaneously reduce the number of degrees of freedom while guaranteeing a symmetric distribution of points.
Symmetries also serve to reduce the number of basis polynomials that must be considered when computing ξ(φ). Consider the following two monomials and defined inside of a square domain with vertices (−1, −1) and (1, 1). We note that p 1 (x, y) = p 2 (y, x). As this is a symmetry which is expressed by the domain it is clear that any symmetric quadrature rule capable of integrating p 1 is also capable of integrating p 2 . Further, the index i is odd we have In both cases it follows that the integral of p 1 is zero over the domain. More importantly, it also follows that any set of symmetric points are also capable of obtaining this result. This is due to terms on the right hand side of Equation 1 pairing up and cancelling out. A consequence of this is that not all of equations in the system specified by Equation 1 are independent.
Having identified such polynomials for a given domain it is legitimate to exclude them from our definition of ξ(φ). Although this exclusion does change the value of ξ(φ) in the case of a non-zero truncation error the effect is not significant. We shall denote the set of basis polynomials which are included as the objective basis, and denote this byP φ .

Reference domains
In the paragraphs which follow we will takeP (α,β) i (x) to refer to a normalised Jacobi polynomial as specified in §18.3 of [24]. In two dimensions we take the coordinate axes to be x = (x, y) and x = (x, y, z) in three dimensions.
Triangle. Our reference triangle can be seen in Figure 1 and has an area given by A triangle has six symmetries: two rotations, three reflections, and the identity transformation. A simple means of realising these symmetries is to transform into barycentric coordinates which are related to Cartesian coordinates via where the columns of the matrix can be seen to be the vertices of our reference triangle. The utility of barycentric coordinates is that the symmetric counterparts to a point λ are given by its unique permutations. The number of unique permutations depends on the number of distinct components of λ and leads us to the following three symmetry orbits where α and β are suitably constrained as to ensure the validity of the resulting coordinates. It can be easily verified that the orthonormal polynomial basis inside of our reference triangle is given by where a = 2(1 + x)/(1 − y) − 1, and b = y with the objective basis being given bỹ In the asymptotic limit the cardinality of the objective basis is half that of the complete basis. However, the modes of this objective basis are known not to be completely independent. Several authors have investigated the derivation of an optimal quadrature basis on the triangle. Details can be found in the papers of Lyness [4] and Dunavant [5].
Quadrilateral. Our reference quadrilateral can be seen in Figure 1. The area is simply A square has eight symmetries: three rotations, four reflections and the identity transformation. Applying these symmetries to a point (α, β) with 0 ≤ (α, β) ≤ 1 will yield a set χ(α, β) containing its counterparts. The cardinality of χ depends on if any of the symmetries give rise to identical points. This can be seen to occur when either β = α or β = 0. Enumerating the possible combinations of the above conditions gives rise to the following four symmetry orbits Trivially, the orthonormal basis inside of our quadrilateral is given by where a = x, and b = y. The objective basis is found to bẽ with a cardinality one eighth that of the complete basis.
Tetrahedron. Our reference tetrahedron is a right-tetrahedron as depicted in Figure 2. Integrating up the volume we find tetrahedron has a total of 24 symmetries. Once again it is convenient to work in terms of barycentric coordinates which are specified for a tetrahedron as and related to Cartesian coordinates via where as with the triangle the columns of the matrix correspond to vertices of the reference tetrahedron. Similarly the symmetric counterparts of λ are given by its unique permutations. This leads us to the following five symmetry orbits where α, β, and γ are constrained to ensure that 0 ≤ λ i ≤ 1 and i λ i = 1.
With some manipulation it can be verified that the orthonormal polynomial basis inside of our reference tetrahedron is given by where a = −2(1 + x)/(y + z) − 1, b = 2(1 + y)/(1 − z), and c = z. The objective basis is given byP Prism. Extruding the reference triangle along the z-axis gives our reference prism of Figure 2. It follows that the volume is There are a total of 12 symmetries. On account of the extrusion the most natural coordinate system is a combination of barycentric and Cartesian coordinates: (λ 1 , λ 2 , λ 3 , z). Let Perm 3 generate all of the unique permutations of its first three arguments. Using this the six symmetry groups of the prism can be expressed as S 6 (α, β, γ) = Perm 3 (α, β, 1 − α − β, ±γ), where the constraints on α and β are identical to those in a triangle and 0 < γ ≤ 1.
Combining the orthonormal polynomial bases for a right-triangle and line segment yields the orthonormal prism basis where a = 2(1 + x)/(1 − y) − 1, b = y, and c = z. The objective basis is given bỹ Pyramid. Our reference pyramid can be seen in Figure 2 with a volume determined by (1−z)/2 (z−1)/2 dx dy dz = 8/3. The symmetries are identical to those of a quadrilateral. Extending the notation employed for the quadrilateral we obtain the following symmetry orbits subject to the constraints that 0 < (α, β) ≤ (1 − γ)/2 and −1 ≤ γ ≤ 1. Inside of the reference pyramid the orthonormal polynomial basis is found to be where a = 2x/(1 − z), b = 2y/(1 − z), and c = z. The objective basis is Hexahedron. Our choice of reference hexahedron can be seen in Figure 2. The volume is, trivially, A hexahedron exhibits octahedral symmetry with a symmetry number of 48. The procedure for determining the orbits similar to that used for the quadrilateral. We consider applying these symmetries to a point (α, β, γ) with 0 ≤ (α, β, γ) ≤ 1 and let the resulting set of points by given by Ξ(α, β, γ). When α, β, and γ are all distinct and greater than zero the set has a cardinality of 48, as expected. However, when one or more parameters are either identical to one another or equal to zero some symmetries give rise to equivalent points. This reduces the cardinality of the set. Enumerating the various combinations we obtain seven symmetry orbits Ξ(0, 0, 0), Trivially, the orthonormal basis inside of our reference hexahedron is given by where a = x, b = y, and c = z. The objective basis is

Methodology
Our methodology for identifying symmetric quadrature rules is a refinement of that described by Witherden and Vincent [12] for triangles. This method is, in turn, a refinement of that of Zhang et al. [9].
To derive a quadrature rule four input parameters are required: the reference domain, Ω, the number of quadrature points N p , the target rule strength, φ, and a desired runtime, t. The algorithm begins by computing all of the possible symmetric decompositions of N p . The result is a set of vectors satisfying the relation where N s is the number of symmetry orbits associated with the domain Ω, and n i j is the number of orbits of type j in the ith decomposition. Finding these involves solving the constrained linear Diophantine equation outlined in section 2. It is possible for this equation to have no solutions. As an example we consider the case when N p = 44 for a triangular domain. From the symmetry orbits we have subject to the constraint that n 1 ∈ {0, 1}. This restricts N p to be either a multiple of three or one greater. Since forty-four is neither of these we find the equation to have no solutions. Therefore, we conclude that there can be no symmetric quadrature rules inside of a triangle with forty-four points.
Given a decomposition we are interested in finding a set of orbital parameters and weights that minimise the error associated with integrating the objective basis on Ω. This is an example of a nonlinear least squares problem. A suitable method for solving such problems is the Levenberg-Marquardt algorithm (LMA). The LMA is an iterative procedure for finding a set of parameters that correspond to a local minima of a set of functions. The minimisation process is not always successful and is dependent on an initial guess of the parameters. Within the context of quadrature rule derivation minimisation can be regarded as successful if ξ(φ) ∼ where represents machine precision.
Let us denote the number of parameters associated with symmetry orbit S i as S i . Using this we can express the total number of degrees of freedom associated with decomposition i as with the second term accounting for the presence of one quadrature weight associated with each symmetry orbit. From the list of orbits given in section 2 we expect the weights contribute approximately one third of the degrees of freedom. This is not an insignificant fraction. One way of eliminating the weights is to treat them as dependent variables. When the points are prescribed the right hand side of Equation 1 becomes linear with respect to the unknowns-the weights. In general, however, the number of weights will be different from the number of polynomials in the objective basis. It is therefore necessary to obtain a least squares solution to the systen. Linear least squares problems can be solved directly through a variety of techniques. Perhaps the most robust numerical scheme is that of singular value decomposition (SVD). Thus, at the cost of solving a small linear least squares problem at each LMA iteration we are able to reduce the number of free parameters to N s j=1 n i j S i .
Such a modification has been found to greatly reduce the number of iterations required for convergence. This reduction more than offsets the marginally greater computational cost associated with each iteration. Previous works [9,10,18] have emphasised the importance of picking a 'good' initial guess to seed the LMA. To this end several methodologies for seeding orbital parameters have been proposed. The degree of complexity associated with such strategies is not insignificant. Further, it is necessary to device a separate strategy for each symmetry orbit. Our experience, however, suggests that the choice of decomposition is far more important than the initial guess in determining whether minimisation will be successful. For larger values of N p we note that many decompositions-especially those for prisms and pyramids-are pathological.
As an example of this we consider searching for an N p = 80 point rule inside of a prism where there are 2 380 distinct symmetrical decompositions. One such Algorithm 1. Procedure for generating symmetric quadrature rules of strength φ with N p points inside of a domain. if ξ ∼ then If minimisation was successful 8: save R 9: end if 10: until CurrentTime() − t 0 > t 11: end for 12: end procedure 13: function RuleResid(R) 14: for all p i ∈P φ do For each basis function 15: b i ← Ω p i (x) dx 16: for all r j ∈ R do For each orbit 17: r j ← ClampOrbit(r j ) Ensure orbital parameters are valid 18: for all x k ∈ ExpandOrbit(r j ) do 20: ω ← b/A Use SVD to determine the weights 25: return Aω − b Compute the residual 26: end function decomposition is N p = 40 |S 2 | where all points lie in a single column down the middle of the prism. Since there is no variation in either x or y it is not possible to obtain a rule of strength φ ≥ 1. Hence, the decomposition can be dismissed without further consideration.
A presentation of our method in pseudocode can be seen in Algorithm 1. When the objective basis functions inP φ are orthonormal Equation 6 states that the integrand of all non-constant modes is zero. We can exploit this to simplify the computation of b i . The purpose of ClampOrbit is to enforce the constraints associated with a given orbit to ensure that all points remain inside of the domain.

Implementation
We have implemented the algorithms outlined above in a C++11 program called polyquad. The program is built on top of the Eigen template library [25] and is parallelised using MPI. It is capable of searching for quadrature rules on triangles, quadrilaterals, tetrahedra, prisms, pyramids, and hexahedra. All rules are guaranteed to be symmetric having all points inside of the domain. Polyquad can also, optionally, filter out rules possessing negative weights. Further, functionality exists, courtesy of MPFR [26], for refining rules to an arbitrary degree of numerical precision and for evaluating the truncation error of a ruleset.
The source code for polyquad is available under the terms of the GNU General Public License v3.0 and can be downloaded from https://github.com/ vincentlab/Polyquad.

Rules
Using polyquad we have derived a set of quadrature rules for each of the reference domains in section 2. All rules are completely symmetric, possess only positive weights, and have all points inside of the domain. It is customary in the literature to refer to quadratures with the last two attributes as being "PI" rules. As polyquad attempts to find an ensemble of rules it is necessary to have a means of differentiating between otherwise equivalent formulae. In constructing this collection the truncation term ξ(φ + 1) was employed with the rule possessing the smallest such term being chosen. The number of points N p required for a rule of strength φ can be seen in Table 1. The rules themselves are provided as electronic supplementary material and have been refined to 38 decimal places.
From the table we note that several of the rules appear to improve over those in the literature. We consider a rule to be an improvement when it either requires fewer points than any symmetric rule described in literature or when existing symmetric rules of this strength are not PI. We note that many of the rules presented by Dunavant for quadrilaterals [14] and hexahedra [21] possess either negative weights or have points outside of the domain. Using polyquad in quadrilaterals we were able to identify PI rules with point counts less than or equal to those of Dunavant at strengths φ = 8, 9, 18, 19, 20. In tetrahedra we were able to reduce the number of points required for φ = 7 and φ = 9 by one and two, respectively, compared with Zhang et al. [9]. Furthermore, in prisms and pyramids rules requiring significantly fewer points than those in literature were identified. As an example, the φ = 9 rule of [18] inside of a prism requires 71 points compared with just 60 for the rule identified by polyquad. Additionally, both of the φ = 10 rules for prisms and pyramids appear to be new.

Conclusions
We have presented a methodology for identifying symmetric quadrature rules on a variety of domains in two and three dimensions. Our scheme does not require any manual intervention and is not restricted to any particular topological configuration inside of a domain. Additionally, it is also capable of generating an ensemble of rules. We have further provided an open source implementation of our method in C++11, and used it to generate a complete set of symmetric quadrature rules that are suitable for use in finite element solvers. All rules possess purely positive weights and have all points inside the domain. Many of the rules appear to be new, and an improvement over those tabulated in the literature.