1 Introduction

Cubature, i.e. multivariate numerical integration, numerically approximates definite integrals over multi-dimensional domains [1,2,3]. This is of interest across a wide range of applications in science and engineering; for example, cubature is fundamental to the calculation of element stiffness matrices in the finite element method.

The most common numerical integration method, and the one considered in this paper, is Gaussian-type cubature, which approximates the integral through a weighted sum of evaluations of the integrand at specific points. Any such set of points and corresponding weights constitutes a cubature rule. Collections of cubature rules have been published [4,5,6], but newer results are still being obtained. Indeed, theoretical results (e.g. [7]) indicate that many rules improving on existing ones are still to be calculated.

The construction of cubature rules is a challenging task because it usually requires solving a strongly nonlinear system of equations called the moment equations. Much of the relevant literature, therefore, focuses on developing better methods to solve the moment equations and obtain cubature rules with specific characteristics (e.g. integration domain or degree of accuracy), for example on the triangle [8,9,10,11,12,13] and the tetrahedron [8, 11, 14,15,16,17].

Invariance of cubature rules with respect to given transformations, i.e. symmetry, is often sought in practical applications. Full symmetry for the triangle and the tetrahedron, for example, where the rule is invariant with respect to any permutation of the vertices, allows for numerical integration that is independent of the ordering of the vertices. Invariance can also significantly simplify the moment equations, allowing them to be solved. For this reason, most known rules on the triangle and tetrahedron are fully symmetric.

Symmetry defines orbits of integration points that are invariant under the chosen transformations. Different types of orbits exist for a given symmetry, thus a rule structure is introduced which indicates the number of orbits of each type. The search for better rules of a given degree is greatly aided by the calculation of consistency conditions [18, 19], which indicate which rule structures are expected to generate solvable moment equations, and can, therefore, also be used to estimate the optimal (in terms of number of points) rule structure.

Consistency conditions for fully symmetric rules on the triangle [20] are widely known and used. For the tetrahedron, a method to derive consistency conditions, and the resulting optimal consistent rule structures, were obtained by Maeztu and Sainz de la Maza [7]. This method however is provided without a detailed implementation and does not lead to explicit expressions for the consistency conditions. This is probably the reason why consistency conditions are not used in recent papers where new rules are obtained on the tetrahedron [8, 11, 14, 15], except for [17].

In this paper, we derive for the first time explicit expressions for the consistency conditions for fully symmetric cubature rules on the tetrahedron. Significantly extending the method in [10], we introduce symmetric polynomials to rigorously generate a new fully-symmetric basis for the polynomial space so that the structure of the moment equations can be easily analyzed to obtain consistency conditions. These in turn are used to generate optimal and quasi-optimal rule structures. We also calculate several new cubature rules, including results that show the limitations in the current approach to generating consistency conditions. While the exposition focuses on fully symmetric rules on the tetrahedron, the proposed approach is general and can form the basis of similar results for different domain types and/or symmetries.

The rest of the paper is separated into five sections: Sect. 2 presents concisely several essential mathematical concepts. Section 3 shows the transformation and decomposition leading to the new basis for symmetric polynomials. In Sect. 4, consistency conditions are derived based on the new basis. An algorithm to calculate consistent rule structure is also presented and explicit formulas to estimate optimal rule structures are given. In Sect. 5, moment equations are solved analytically using Gröbner bases and several interesting results are presented. Finally, Sect. 6 summarises the main results obtained in this paper and mentions possible future research directions.

2 Theoretical background

2.1 Barycentric coordinates

Consider a tetrahedron defined by four vertices with Cartesian coordinates \((x_{|\kappa },y_{|\kappa },z_{|\kappa })\) where \(\kappa =1 \ldots 4\). A point with Cartesian coordinates (x, y, z) can be described, in relation to the tetrahedron, by the barycentric coordinates \((L_1, L_2, L_3, L_4)\) such that

$$\begin{aligned} x&= L_1 x_{|1} + L_2 x_{|2} + L_3 x_{|3} + L_4 x_{|4} \end{aligned}$$
(1a)
$$\begin{aligned} y&= L_1 y_{|1} + L_2 y_{|2} + L_3 y_{|3} + L_4 y_{|4} \end{aligned}$$
(1b)
$$\begin{aligned} z&= L_1 z_{|1} + L_2 z_{|2} + L_3 z_{|3} + L_4 z_{|4} \end{aligned}$$
(1c)
$$\begin{aligned} 1&= L_1+L_2+L_3+L_4 \end{aligned}$$
(1d)

Barycentric coordinates are not independent, as seen by Eq. (1d). Although their use increases the number of unknowns, it greatly simplifies the calculations when considering symmetry. Barycentric coordinates can also be defined without the normalisation (1d); when such a normalisation is adopted, as in this paper, the resulting coordinates are also called volume coordinates.

A polynomial of degree d in the Cartesian coordinates can be expressed as a homogeneous polynomial in the barycentric coordinates, i.e. as a linear combination of monomials \(L_1^iL_2^jL_3^kL_4^{d-i-j-k}\).

2.2 Symmetric polynomials

A symmetric polynomial is a multivariate polynomial invariant under any permutation of its variables [21]. Considering n variables \(v_1, v_2,\ldots , v_n\), we define the elementary symmetric polynomials \({\tilde{v}}_k\) as the sum of all products of k distinct variables \(v_i\), with negative sign when k is odd, that is

$$\begin{aligned} {\tilde{v}}_k=(-1)^k\sum _{1 \le i_1<i_2< \cdots <i_k \le n}v_{i_1}v_{i_2} \cdots v_{i_k} \end{aligned}$$
(2)

with \({\tilde{v}}_0=1\). The summation in (2) is taken over all ordered sets of k distinct indices in the range \(1 \ldots n\). Any symmetric polynomial in the variables \(v_i\) can be uniquely expressed as a polynomial in the elementary symmetric polynomials \({\tilde{v}}_k\). The values of \(v_i\) can be calculated from \({\tilde{v}}_k\) as the solutions for v of the polynomial equation

$$\begin{aligned} \sum _{j=0}^{n}{\tilde{v}}_{n-j}v^j=0. \end{aligned}$$
(3)

2.3 Solution of a polynomial system

Consider a generic system of m polynomial equations with n variables \(v_1,v_2,\ldots ,v_n\) and real coefficients. The system is called overdetermined if \(m>n\) and underdetermined if \(m<n\).

A solution of a polynomial system is a tuple of (possibly complex) values \((v_1,v_2,\ldots ,v_n)\) that satisfy all equations in the system. A system is called inconsistent if it has no solution, and consistent if it has at least one solution. A consistent system is zero-dimensional if it has a finite number of solutions and positive-dimensional if it has an infinite number of solutions. For more details on polynomial system solving see [22].

Numerical methods, such as Newton’s method, can obtain numerical approximations to individual solutions but cannot systematically obtain all solutions. Indeed, for positive-dimensional systems there is no clear answer to what is the complete solution. Using algebraic geometry, algebraic solutions to a system can be obtained by expressing the system in a form that is exact and easy to solve numerically, such as Gröbner bases or, for zero-dimensional systems, the rational univariate representation.

2.4 Cubature rules, moment equations and consistency conditions

For a function f over a domain \(\Omega \subset {\mathbb {R}}^n\) with n-volume V, consider the scaled integral

$$\begin{aligned} I(f)= \frac{1}{V}\int _{\Omega }fd\Omega . \end{aligned}$$
(4)

A cubature formula (or cubature rule) is an approximation of the integral I as

$$\begin{aligned} Q(f) = \sum _{i=1}^{n_K} w_i f(\varvec{x}_i) \approx I(f), \end{aligned}$$
(5)

where \(f(\varvec{x}_i)\) is the value of the function f at the point \(\varvec{x}_i \in {\mathbb {R}}^n\), \(w_i\) is the corresponding weight, and \(n_K\) is the number of integration points. Only cubature rules of (polynomial) degree d are considered here, for which Eq. (5) is exact for all polynomials of degree equal or less than d and not exact for at least one polynomial of degree \(d+1\). The accuracy of the approximation (5) is not considered in this paper, details can be found in [1].

Let \({\mathbb {P}}^n_d\) denote the vector space of all polynomials in n independent variables of degree at most d, and B be a basis of this vector space. The cubature rule Q is of degree at least d if and only if it is exact for any element of B, that is

$$\begin{aligned} Q(f)=I(f) \quad \forall f \in B. \end{aligned}$$
(6)

Equation (6) is a polynomial system of \(\dim ({\mathbb {P}}^n_d)\) equations in \((n+1)n_K\) unknowns, where \(\dim ({\mathbb {P}}^n_d)\) is the dimension of \({\mathbb {P}}^n_d\). These so-called moment equations can be solved to obtain the coordinates of the integration points and their associated weights, which define the cubature rule.

The consistency conditions are assumed conditions on the number of integration points for a given rule and degree, based on the (not necessarily true) assumption that the system of moment equations, and any of its subsystems, will be consistent if and only if it is not overdetermined. These are, therefore, consistency conditions for non-overdeterminedness of the moment equations; for conciseness, and following previous works, throughout this paper we use the simpler term “consistency conditions”.

2.5 Quality of cubature rules

A quality is assigned to each cubature rule, using two letters. The first letter refers to the integration weights. It is P if all weights are positive, N if at least one weight is negative but all weights are real, and C if there is at least one complex weight. The second letter refers to the position of the integration points with respect to the integration domain. It is I if all points are inside the integration domain, B if there is at least one point on the boundary of the domain and the remaining are inside, O if at least one point is located outside the integration space but all points have real coordinates, and C if at least one point has complex coordinates. The possible qualities, in the usually assumed order of preference, are: PI, NI, PB, NB, PO, NO, PC, NC, CC.

2.6 Invariant cubature rules

Let G be a finite group of transformations \(g:{\mathbb {R}}^n \rightarrow {\mathbb {R}}^n\). A function \(f(\varvec{x})\) is invariant with respect to G if it does not change under any transformation of the group, that is \(f(g(\varvec{x}))=f(\varvec{x}) \;\forall g \in G\) [23]. The G-orbit of a point \(\varvec{x}\in {\mathbb {R}}^n\), denoted by \(G(\varvec{x})\), is the set \(\{g(\varvec{x}): g\in G\}\). The point \(\varvec{x}\) is called a generator of \(G(\varvec{x})\).

A cubature rule is invariant with respect to the group G if the region \(\Omega\) is G-invariant (i.e. \(g(\Omega )=\Omega \;\forall g \in G\)), the set of integration points is a union of G-orbits, and all points in the same orbit have the same weight. Using Sobolev’s theorem (see [2] and references therein), invariant cubature rules can be obtained by solving the moment equations

$$\begin{aligned} Q(f)=I(f) \quad \forall f \in B_G \end{aligned}$$
(7)

where \(B_G\) is a basis of the vector space \({\mathbb {P}}^n_d(G) \subset {\mathbb {P}}^n_d\) of all G-invariant polynomials of degree d in n variables.

3 Bases for cubature on the tetrahedron

3.1 Asymmetric basis

In three dimensions, a basis of \({\mathbb {P}}^3_d\) is the set of monomials \(x^i y^j z^k\) with total degree \(i+j+k \le d\), which has dimension

$$\begin{aligned} m_{\textrm{asym}}=\frac{(d+3)(d+2)(d+1)}{6}. \end{aligned}$$
(8)

In this case, there are four unknowns for every integration point (three coordinates and one weight) so the consistency condition is

$$\begin{aligned} 4 n_K \ge m_{\textrm{asym}}. \end{aligned}$$
(9)

This consistency condition is obtained for any domain in \({\mathbb {R}}^3\), for asymmetric rules (i.e. for rules that are not necessarily invariant). We know that there exist rules that do not follow this consistency condition: for degree \(d=9\), Eq. (9) yields \(n_K \ge 55\), but there are, for example, known cubature rules with 52 points on the cube [19] and with 53 points on the tetrahedron [24].

For tetrahedra, a convenient alternative basis is the set of monomials \(L_1^i L_2^j L_3^k L_4^{d-i-j-k}\) of total degree d. The use of barycentric coordinates simplifies the expression of the integration point coordinates in a way independent of a specific tetrahedron. Even more importantly, barycentric coordinates simplify expressing invariant rules, as any symmetry of the tetrahedron can be expressed as invariance with respect to specific permutations of the vertices, or of the barycentric coordinates.

3.2 Fully symmetric basis

We consider here fully symmetric rules on tetrahedra, which are invariant with respect to the action of the symmetric group \(\textrm{S}_4\) (the group of all permutations of four elements). For given degree d, the \(\textrm{S}_4\)-invariant polynomials are the symmetric polynomials in the four barycentric coordinates. A basis of these polynomials are the products \({\tilde{L}}_1^{t}{\tilde{L}}_2^{i}{\tilde{L}}_3^{j}{\tilde{L}}_4^{k}\) with \(t+2i+3j + 4k=d\), where the elementary symmetric polynomials in the barycentric coordinates are given by Eq. (2) as

$$\begin{aligned} {\tilde{L}}_1&= -(L_1+L_2+L_3+L_4) = -1 \end{aligned}$$
(10a)
$$\begin{aligned} {\tilde{L}}_2&= L_1L_2+L_1L_3+L_1L_4+L_2L_3+L_2L_4+L_3L_4 \end{aligned}$$
(10b)
$$\begin{aligned} {\tilde{L}}_3&= -(L_1L_2L_3+L_1L_3L_4+L_1L_2L_4+L_2L_3L_4) \end{aligned}$$
(10c)
$$\begin{aligned} {\tilde{L}}_4&= L_1L_2L_3L_4. \end{aligned}$$
(10d)

Due to Eq. (10a), the \(\textrm{S}_4\)-invariant basis is actually the products \({\tilde{L}}_2^{i}{\tilde{L}}_3^{j}{\tilde{L}}_4^{k}\) with \(2i+3j+4k \le d\). The number of elements in the basis, and, therefore, also the number of moment equations \(m_{\textrm{e}}\), is the number of non-negative integer solutions to \(2i+3j+4k \le d\), given by [25]

$$\begin{aligned} m_{\textrm{e}} = \left\lfloor {\frac{(d+4)^3 + 3(d+4)^2 - 9(d+4)\bigl ((d+4) \bmod 2\bigr )}{144}}\right\rceil \end{aligned}$$
(11)

where \({\lfloor }x{\rceil }\) denotes the nearest integer to x. Comparing Eqs. (11) and (8) shows that the fully symmetric case has significantly fewer moment equations than the asymmetric case (a ratio of 1/24 as \(d \rightarrow \infty\)).

The symmetric group \(\textrm{S}_4\) generates five different types of orbits, depending on the number of repeated barycentric coordinates in the generator. These orbit types, numbered from 0 to 4, are shown in Table 1.

Table 1 Types of orbits
Full size table

Table 1 also shows the number of points for each orbit type, and the number of variables introduced at the moment equations for each orbit of a given type (this is equal to the number of variables defining the generator, plus one variable which is the weight, common to all points in the orbit).

The orbit structure of a rule is the list \([n_0,n_1,n_2,n_3,n_4]\), where \(n_i\) is the number of orbits of type i. The total number of points \(n_K\) for a rule with orbit structure \([n_0,n_1,n_2,n_3,n_4]\) is:

$$\begin{aligned} n_K=n_0+4n_1+6n_2+12n_3+24n_4, \end{aligned}$$
(12)

while the number of unknowns is

$$\begin{aligned} n_e=n_0+2n_1+2n_2+3n_3+4n_4. \end{aligned}$$
(13)

Equations (11) and (13) can be used to derive a consistency condition for fully symmetric rules. A more precise set of consistency conditions can however be obtained by adopting a different basis that has as many elements as possible that are zero for as many orbit types as possible.

3.3 A simpler fully symmetric basis

Substituting L for the generic variable v in Eq. (3), and using Eq. (10a), gives

$$\begin{aligned} L^4 - L^3 + {\tilde{L}}_2 L^2 + {\tilde{L}}_3 L + {\tilde{L}}_4 = 0. \end{aligned}$$
(14)

Studying the multiplicity of the roots of (14), considered as a quartic function in L, gives the relation between \({\tilde{L}}_2\), \({\tilde{L}}_3\) and \({\tilde{L}}_4\) for each orbit type. This study is greatly simplified by considering the depressed quartic, therefore we use the transformation

$$\begin{aligned} l_{\kappa }=L_{\kappa }-1/4 \quad \hbox { with}\quad \kappa =1 \cdots 4 \end{aligned}$$
(15)

so that the elementary symmetric polynomials are

$$\begin{aligned} {\tilde{l}}_1&= -(l_1+l_2+l_3+l_4) = 0 \end{aligned}$$
(16a)
$$\begin{aligned} {\tilde{l}}_2&= l_1 l_2 + l_1 l_3 + l_1 l_4 + l_2 l_3 + l_2 l_4 + l_3 l_4 \end{aligned}$$
(16b)
$$\begin{aligned} {\tilde{l}}_3&= -(l_1 l_2 l_3 + l_1 l_3 l_4 + l_1 l_2 l_4 + l_2 l_3 l_4) \end{aligned}$$
(16c)
$$\begin{aligned} {\tilde{l}}_4&= l_1 l_2 l_3 l_4 \end{aligned}$$
(16d)

and Eq. (3) becomes

$$\begin{aligned} l^4 + {\tilde{l}}_2 l^2 + {\tilde{l}}_3 l + {\tilde{l}}_4 = 0. \end{aligned}$$
(17)

The discriminant of (17) with respect to l is

$$\begin{aligned} \Delta =-27{\tilde{l}}_3^4-4{\tilde{l}}_2({\tilde{l}}_2^2-36{\tilde{l}}_4){\tilde{l}}_3^2+16{\tilde{l}}_4({\tilde{l}}_2^2-4{\tilde{l}}_4)^2. \end{aligned}$$
(18)

We, therefore, have the following cases [26]:

  • Type-0 orbit (four equal roots, all zero) for \({\tilde{l}}_2={\tilde{l}}_3={\tilde{l}}_4=0\)

  • Type-1 orbit (three equal roots) for \(8 {\tilde{l}}_2^3 + 27 {\tilde{l}}_3^2=0\) and \({\tilde{l}}_2^2 + 12 {\tilde{l}}_4=0\), with \({\tilde{l}}_2 \ne 0\)

  • Type-2 orbit (two pairs of equal roots) for \({\tilde{l}}_3=0\) and \({\tilde{l}}_2^2 - 4 {\tilde{l}}_4=0\), with \({\tilde{l}}_2 \ne 0\)

  • Type-3 orbit (only one pair of equal roots) for \(\Delta = 0\) but none of the previous cases holding

  • Type-4 orbit (four distinct roots) for \(\Delta \ne 0\).

Note that the conditions given for orbits of types 0, 1 and 2 ensure that \(\Delta =0\). Further simplification is achieved by introducing the quantities

$$\begin{aligned} p=-\frac{2{\tilde{l}}_2}{3}, q=-{\tilde{l}}_3, r=\frac{{\tilde{l}}_2^2 +12{\tilde{l}}_4}{9}, \end{aligned}$$
(19)

resulting in the following simpler conditions

  • Type-0 orbit for \(p=q=r=0\)

  • Type-1 orbit for \(p^3 - q^2=0\) and \(r=0\), with \(p \ne 0\)

  • Type-2 orbit for \(q=0\) and \(p^2-r=0\), with \(p \ne 0\)

  • Type-3 orbit for \(\Delta = 0\) but none of the previous cases holding

  • Type-4 orbit for \(\Delta \ne 0\)

where now

$$\begin{aligned} \Delta = 27 \bigl ( -q^4 + 2 p (p^2 - 3r) q^2 - (p^2 - 4 r) (p^2 - r)^2 \bigr ). \end{aligned}$$
(20)

We, therefore, consider the fully symmetric monomial basis \(p^i q^j r^k\) with weighted degree \(2i+3j+4k \le d\). This basis is simpler than the previous ones, as it is easier to express the conditions holding on orbits of types 0 to 3.

3.4 Fully symmetric basis for consistency conditions

As already mentioned, a more precise set of consistency conditions can be obtained by adopting a basis that has as many elements as possible that are zero for as many orbit types as possible. We create such a basis starting from the monomial basis \(p^i q^j r^k\) and then splitting and transforming (by taking linear combinations with constant coefficients) groups of elements. This process is summarised in Fig. 1 and explained in the following.

Fig. 1
figure 1

Diagram showing the derivation of the fully symmetric basis for consistency conditions

Full size image

Orbits of types 0 to 3 are identified by the condition \(\Delta =0\). Equation (20) shows that \(\Delta\) is of degree 3 in r, so we split \(p^i q^j r^k\), by degree of r, into

$$\begin{aligned}{}[p^i q^j r^k] \rightarrow [p^i q^j r^k r^3, p^i q^j r^2, p^i q^j, p^i q^j r] \end{aligned}$$
(21)

where i, j, k can take different values for each term. To simplify notation, here and in the following we implicitly consider that \(i,j,k \ge 0\) and that the monomials and polynomials shown are of total weighted degree up to d. The reason for swapping the last two terms will be seen shortly.

The terms \(p^i q^j r^k r^3\) can easily be transformed, by taking linear combinations with the other terms, into elements \(p^i q^j r^k \Delta\), which are zero for orbit types 0 to 3. The remaining terms, of degree less than 3 in r, cannot be linearly combined to give a polynomial with a factor \(\Delta\), and will, therefore, not be zero for type-3 orbits.

To proceed further, we use algebraic geometry [27] to obtain that any polynomial in p, q, r that is zero for orbit types 0 to 2 will be the sum of polynomials with factors \(p^3-pr-q^2\), \(r(p^2-r)\), or qr (as these three polynomials generate the radical of the product of the ideals generated by the polynomials that have to be zero for each orbit type from 0 to 2).

The terms \(p^i q^j r^2\) in (21), taken in linear combination with the terms \(p^i q^j r\), give terms \(p^i q^j r (p^2-r)\). This leaves from (21) the terms \(p^i q^j\) and \(p^i q^j r\), which are split as

$$\begin{aligned}{}[p^i q^j, p^i q^j r] \rightarrow [p^i q^j q^2, p^i q, p^i, p^i q^j qr, p^i r] \end{aligned}$$
(22)

The first term on the r.h.s. of (22), taken in linear combination with the second, fourth and fifth terms, yields terms \(p^i q^j (p^3-pr-q^2)\), while the fourth term is already in the form \(p^i q^j (qr)\), completing the set of terms that are zero for orbit types 0 to 2.

The third term on the r.h.s. of (22) can be further split as

$$\begin{aligned}{}[p^i] \rightarrow [p^i p^2, p, 1], \end{aligned}$$
(23)

where the terms \(p^i p^2\), taken in linear combination with the terms \(p^i r\), yields the terms \(p^i (p^2-r)\) which are zero for orbit types 0 and 2.

Collecting all the resulting terms yields a new fully symmetric basis, summarised in Table 2, which no longer contains only monomial terms. The way in which the new basis is derived from the monomial basis in p, q, r ensures that the resulting basis is indeed a basis of the same vector space of polynomials. Additionally, this basis maximises the number of elements that are zero for different orbit types, since no linear combination of elements that are zero for fewer orbit types can give a polynomial that is zero for more orbit types.

Table 2 Fully symmetric basis for consistency conditions
Full size table

4 Consistency conditions and (quasi-)optimal rules

In the previous section, we obtained a new fully symmetric basis for the vector space of \(\textrm{S}_4\)-invariant polynomials, which has as many elements as possible that are zero for as many orbit types as possible. This allows for deriving the consistency conditions, and, therefore, also the estimated lower bounds on the number of integration points for cubature rules in tetrahedra.

4.1 Number of basis elements equations

The last column of Table 2 gives the number of basis elements for each element type. This is the number of non-negative integer solutions, for the indices appearing in the weighted degree, for which the weighted degree is less or equal to the degree d. Specifically, \(m_{p3}(d)\) is given by eq. (11), extended to also cover negative values of d

$$\begin{aligned} m_{p3}(d)= & {} \left\lfloor \frac{(d+4)^3 + 3(d+4)^2 - 9(d+4)\bigl ((d+4) \bmod 2\bigr )}{144}\right\rceil \nonumber \\{} & {} \llbracket {d \ge 0}{\rrbracket } \end{aligned}$$
(24a)

where the Iverson brackets \(\llbracket {\ldots }{\rrbracket }\) are defined as [28]

$$\begin{aligned} \llbracket {S}{\rrbracket } = {\left\{ \begin{array}{ll} 0 &{} \text {if }\,S\, \text { is false} \\ 1 &{} \text {if }\, S\, \text { is true} \end{array}\right. } \end{aligned}$$
(24b)

Similarly, \(m_{p2}(d)\) is the number of non-negative integer solutions of \(2i+3j \le d\), given by [29]

$$\begin{aligned} m_{p2}(d) = \left\lfloor { \frac{(d + 3)^2}{12} }\right\rceil \llbracket {d \ge 0}{\rrbracket } \end{aligned}$$
(24c)

and \(m_{p1}(d)\) is the number of non-negative integer solutions of \(2i \le d\), that is

$$\begin{aligned} m_{p1}(d) =\left\lfloor { \frac{d+2}{2} }\right\rfloor \llbracket {d \ge 0}{\rrbracket } \end{aligned}$$
(24d)

where \({\lfloor }{x}{\rfloor }\) is the largest integer that is smaller or equal to x. Finally, \(m_{p0}(d)\) is simply given by

$$\begin{aligned} m_{p0}(d) = \llbracket {d \ge 0}{\rrbracket } \end{aligned}$$
(24e)

Table 2 shows that the basis elements can be grouped by the orbit types for which they are not necessarily zero, giving six groups: \(\{4\}, \{3,4\}, \{2,3,4\}, \{1,3,4\}, \{1,2,3,4\}, \{0,1,2,3,4\}\), having respectively \(m_4, m_3, m_2, m_1, m_{12}, m_0\) elements. From Table 2 and Eq. (24), the number of basis elements in each group is computed as

$$\begin{aligned} m_4&= m_{p3}(d-12) \end{aligned}$$
(25a)
$$\begin{aligned} m_3&= \left\lfloor {\Bigl (\frac{d}{2}-2 \Bigr )^2}\right\rfloor \llbracket {d\ge 6}{\rrbracket } \end{aligned}$$
(25b)
$$\begin{aligned} m_2&= \left\lfloor { \frac{d}{2} - 1 }\right\rfloor \llbracket {d \ge 4}{\rrbracket } \end{aligned}$$
(25c)
$$\begin{aligned} m_1&= (d - 2)\llbracket {d \ge 2}\rrbracket \end{aligned}$$
(25d)
$$\begin{aligned} m_{12}&= \llbracket {d \ge 2}\rrbracket \end{aligned}$$
(25e)
$$\begin{aligned} m_0&= 1 \end{aligned}$$
(25f)

Table 3 lists the number of basis elements in each element group, and the total number of elements, for degrees \(d \le 20\).

Table 3 Number of basis elements for each orbit group
Full size table

4.2 Consistency conditions

Every element inside the basis must satisfy Eq. (6), therefore each element group represents a subsystem of the moment equations. For each of these groups, consistency requires that the number of unknowns across all orbits in the group must be greater or equal to the number of equations that only involve the orbits of the group. Applying this to each group we get the consistency conditions

$$\begin{aligned}&n_0+2n_1+2n_2+3n_3+4n_4 \ge m_4 \nonumber \\&\quad + m_3 + m_2 + m_1 + m_{12} + m_0 \end{aligned}$$
(26a)
$$\begin{aligned}&2n_1+2n_2+3n_3+4n_4\ge m_4 + m_3 + m_2 + m_1 + m_{12} \end{aligned}$$
(26b)
$$\begin{aligned}&2n_1 +3n_3+4n_4\ge m_4 + m_3 + m_1 \end{aligned}$$
(26c)
$$\begin{aligned}&2n_2+3n_3+4n_4\ge m_4 + m_3 + m_2 \end{aligned}$$
(26d)
$$\begin{aligned}&3n_3+4n_4\ge m_4 + m_3 \end{aligned}$$
(26e)
$$\begin{aligned}&4n_4\ge m_4 \end{aligned}$$
(26f)

An additional condition is that there can be at most one orbit of type 0, i.e. \(n_0 \in \{0,1\}\). This means, since \(m_0=1\), that Eq. (26b) can be omitted as it is implied by Eq. (26a). The consistency conditions can then be written. as

$$\begin{aligned}&n_0 \in \{0,1\} \end{aligned}$$
(27a)
$$\begin{aligned}&n_0+2n_1+2n_2+3n_3+4n_4 \ge m_4 \nonumber \\&\quad + m_3 + m_2 + m_1 + m_{12} + m_0 \end{aligned}$$
(27b)
$$\begin{aligned}&2n_1 +3n_3+4n_4\ge m_4 + m_3 + m_1 \end{aligned}$$
(27c)
$$\begin{aligned}&2n_2+3n_3+4n_4\ge m_4 + m_3 + m_2 \end{aligned}$$
(27d)
$$\begin{aligned}&3n_3+4n_4\ge m_4 + m_3 \end{aligned}$$
(27e)
$$\begin{aligned}&4n_4\ge m_4 \end{aligned}$$
(27f)

4.3 Consistent rule structures

For a given degree d, a rule structure \([n_0,n_1,n_2,n_3,n_4]\) is consistent if it satisfies the consistency conditions (27). Consistent rule structures with a given maximum number of points can be easily calculated using Algorithm 1.

figure a

A consistent rule structure is optimal for a given degree if there are no other consistent structures for the same degree with fewer points. Finding optimal consistent rule structures is in general an integer linear programming problem [19], but we propose here a simpler approach. Considering the consistency conditions (27) in reverse order, we estimate an optimal consistent rule structure as

$$\begin{aligned} n_4^*&=\left\lceil {\frac{m_4}{4}}\right\rceil \end{aligned}$$
(28a)
$$\begin{aligned} n_3^*&=\left\lceil {\frac{m_4 + m_3-4 n_4^*}{3}}\right\rceil \end{aligned}$$
(28b)
$$\begin{aligned} n_2^*&=\left\lceil {\frac{m_4 + m_3 + m_2 - 3 n_3^* - 4 n_4^*}{2}}\right\rceil \end{aligned}$$
(28c)
$$\begin{aligned} n_1^*&=\left\lfloor {\frac{m_{\textrm{e}}-2n_2^*-3n_3^*-4n_4^*}{2}}\right\rfloor \end{aligned}$$
(28d)
$$\begin{aligned} n_0^*&=m_{\textrm{e}}-2n_1^*-2n_2^*-3n_3^*-4n_4^* \end{aligned}$$
(28e)

Algorithm 1, with a maximum number of points equal to the number of points of rule structure (28), shows that this rule is optimal, and is unique (at least up to an unrealistically high degree \(d=200\)).

Table 4 shows the optimal consistent rule structures for degree \(d \le 20\), and the corresponding number of points. This number of points represents an estimate of the lower bound for the number of integration points.

Table 4 Optimal consistent rule structures for tetrahedra
Full size table

It is, however, important to remember that consistency conditions only provide an estimate for which cubature rule structures will yield actual cubature rules, since the satisfaction of the consistency conditions does not guarantee that the system of moment equations is indeed consistent, and, therefore, a cubature rule with a given structure actually exists.

Additionally, in most practical applications only rules of quality PI, or at most NI, are considered acceptable. Even when rules with the optimal consistent structure exist, therefore, the quality of such rules may not be acceptable. It is then worth looking for quasi-optimal consistent rule structures, i.e. structures that follow the consistency conditions but have a few more integration points than the optimal ones. For a given degree, the search is in practice limited to rules that have fewer integration points than any known PI rule. Quasi-optimal consistent rule structures are easily computed using Algorithm 1.

5 New results for cubature rules

Consistency conditions help limit the search space when searching for cubature rules with better quality or a lower number of points than existing ones. Additionally, the choice of an appropriate basis of the polynomials simplifies the calculations needed to solve the moment equations and derive individual rules.

5.1 Summary of new results

Table 5 shows, for degrees up to 20, the (estimated) lowest number of integration points as well as the lowest number of points achieved in known cubature rules of different quality. Known results for NI or *O (i.e. PO or NO) quality, and for PI quality up to degree 6, are given in [6]. More recent results, all of PI quality, are from [8] for degree 8, [11] for degrees 7 and 9, and [17] for degrees 12 to 20. Results for degrees 10 and 11 are taken from the source code of version 0.9.7 of the PHG (Parallel Hierarchical Grid) code (http://lsec.cc.ac.cn/phg/download.htm) and were obtained using an improved version of the code described in [8].Footnote 1

For degrees 1–6 and 8, rules of PI or NI quality with the estimated lower number of points were already given in the literature. For other degrees, a gap exists between estimated and an actual lower number, which increases with a degree. Using the results obtained in the previous sections and extending the approach presented in [10], some new results were obtained, indicated with an underline in table 5.

Table 5 Estimated and known lower bounds for number of integration points in fully symmetric rules on the tetrahedron
Full size table

For degrees 7 and 9 we obtain rules with the optimal number of points, but only with complex point coordinates. While not of practical interest, these results confirm that consistency of the moment equations is correctly predicted for these rule structures. As the moment equations for both cases are zero-dimensional and all solutions are obtained, using Gröbner bases, we prove that there are no optimal rules of better quality. For degree 9, a new 55-point NI rule with structure [1, 3, 1, 3, 0] is also obtained, which improves on the 59 points of the existing PI rule (see Table 6).

Table 6 Degree 9, 55-point NI rule generators and weights
Full size table

For degree 10, the optimal consistent rule structure given in Table 4 is [0, 5, 2, 3, 0], with 68 points. The corresponding moment equations, using the non-monomial basis from Sect. 3.4, have a subsystem for the type-3 orbits with 9 equations and 9 unknowns which is however inconsistent. There are, therefore, no degree-10 rules with structure \([*,*,*,3,0]\); eliminating those from the list of quasi-optimal consistent structures results in an optimal consistent structure [1, 4, 1, 4, 0] with 71 points. Similar calculations show that there is no rule of degree 11 with structure \([*,*,*,4,0]\), so the optimal consistent structure is [0, 5, 1, 5, 0] with 86 points.

5.2 Details on the consistency conditions for degree 10

For a degree-10 rule with structure \([*,*,*,3,0]\), Table 2 shows that there are 9 basis elements that only involve type-3 orbits, namely

$$\begin{aligned} \begin{aligned}{}[g_1, \ldots , g_9]&= [ p^3 - pr - q^2, qr, p (p^3 - pr - q^2),\\&\qquad r(p^2 - r), q(p^3 - pr - q^2), \\&\qquad pqr, p^2(p^3 - pr - q^2), q^2 r, pr(p^2 - r)] \end{aligned} \end{aligned}$$
(29)

Therefore, the corresponding subsystem of the moment equations becomes

$$\begin{aligned} Q_i = I_i \quad i=1 \ldots 9 \end{aligned}$$
(30)

where

$$\begin{aligned} Q_i = 12 \sum _{j=1}^{3} w_j g_i(p_j,q_j,r_j) \end{aligned}$$
(31)

and \(I_i\) are the exact integrals (4) evaluated for the elements \(g_i\), resulting in

$$\begin{aligned} \begin{aligned}{}[I_1, \ldots , I_9] = [&1/22680, 1/151200, 19/4989600, \\&23/9979200, 1/1108800, 13/19958400, \\&1/2620800, 29/172972800, 17/74131200 ] \end{aligned} \end{aligned}$$
(32)

The quantities \(p_j,q_j,r_j\) for every orbit satisfy \(\Delta _j = \Delta (p_j,q_j,r_j)=0\), where the discriminant \(\Delta\) is given by Eq. (20). Therefore, the system (30) is essentially a system of 9 equations with nine unknowns, which we, therefore, assume to be consistent.

Computing a Gröbner basis for the system, however, shows that the system is inconsistent. In this case, therefore, the consistency conditions fail to predict the consistency of the moment equations. This is not due to an incorrect choice of the non-monomial basis in Sect. 3.4, but due to a more complex relationship between the quantities \(Q_i\). Indeed, we can calculate that in this case the \(Q_i\) are not independent, but satisfy the equation

$$\begin{aligned}{} & {} Q_1 Q_3 Q_7^2 + Q_1 Q_3 Q_7 Q_8 - 5 Q_1 Q_3 Q_7 Q_9 \\{} & {} \quad + 4 Q_1 Q_3 Q_9^2 - 4 Q_1 Q_4 Q_7^2 - 4 Q_1 Q_4 Q_7 Q_8 \\{} & {} \quad + 20 Q_1 Q_4 Q_7 Q_9 - 16 Q_1 Q_4 Q_9^2 - Q_1 Q_5^2 Q_7 \\{} & {} \quad + 8 Q_1 Q_5 Q_6 Q_7 - 8 Q_1 Q_5 Q_6 Q_9 \\{} & {} \quad - 12 Q_1 Q_6^2 Q_7 + 4 Q_1 Q_6^2 Q_8 + 12 Q_1 Q_6^2 Q_9 \\{} & {} \quad + 4 Q_2^2 Q_7^2 + 4 Q_2^2 Q_7 Q_8 - 20 Q_2^2 Q_7 Q_9 \\{} & {} \quad + 16 Q_2^2 Q_9^2 + 8 Q_2 Q_3 Q_5 Q_9 - 8 Q_2 Q_3 Q_6 Q_7 \\{} & {} \quad - 8 Q_2 Q_3 Q_6 Q_8 + 8 Q_2 Q_3 Q_6 Q_9 \\{} & {} \quad - 8 Q_2 Q_4 Q_5 Q_7 + 32 Q_2 Q_4 Q_6 Q_7 - 32 Q_2 Q_4 Q_6 Q_9 \\{} & {} \quad - Q_3^3 Q_7 - Q_3^3 Q_8 + 5 Q_3^3 Q_9 \\{} & {} \quad + 4 Q_3^2 Q_4 Q_7 + 4 Q_3^2 Q_4 Q_8 - 28 Q_3^2 Q_4 Q_9 \\{} & {} \quad + Q_3^2 Q_5^2 - 8 Q_3^2 Q_5 Q_6 + 16 Q_3^2 Q_6^2 \\{} & {} \quad + 4 Q_3 Q_4^2 Q_7 + 32 Q_3 Q_4^2 Q_9 + 8 Q_3 Q_4 Q_5 Q_6 \\{} & {} \quad - 32 Q_3 Q_4 Q_6^2 - 16 Q_4^3 Q_7 + 16 Q_4^2 Q_6^2 = 0 \end{aligned}$$
(33)

Equation (33) is non-linear, and does not hold for \(n_3 > 3\). It is, therefore, clear that deriving consistency conditions that correctly indicate that there are no degree-10 rules with structure \([*,*,*,3,0]\) requires a different approach to the one in this paper (and in general in the literature) which is based on linearly independent basis elements.

6 Conclusion

In this paper, we have rigorously developed a new non-monomial fully symmetric polynomial basis for the tetrahedron. By having as many elements as possible to be zero for as many orbit types as possible, this basis directly leads to the formulation of consistency conditions. For the first time, we are able to obtain explicit formulas for the consistency conditions, and thus for determining the optimal consistent rule structures. Additionally, an algorithm is presented that generates quasi-optimal rule structures.

The new basis is also useful in calculating specific cubature rules, since it allows decomposing the moment equations into a series of successively independent smaller subsystems, which can be exploited in designing more efficient solution methods. Solving the moment equations for specific cases, we obtained a new NI rule of degree 9 with 55 points (lower than existing rules of PI/NI quality). We also proved that there exist rules of degree 7 and 9 with the optimal structure but they are not of practical use as their point coordinates are complex numbers.

Finally, we proved that the optimal rule structures estimated by our formulas for degrees 10 and 11 lead to inconsistent moment equations. This is not due to an incorrect derivation of the consistency conditions but is a result of a more complex non-linear relationship between the moment equations, which cannot be captured by the usual assumptions that are employed to derive consistency conditions.

The quasi-optimal rule structures obtained can be used as the starting point for calculating additional cubature rules. For higher degrees the large number of equations and unknowns necessitates improved solving techniques, possibly mirroring those already developed for the triangle. The overall approach described in this paper can be further applied to obtain consistency conditions for different domains and types of symmetry.