Fast Numerical Integration on Polytopic Meshes with Applications to Discontinuous Galerkin Finite Element Methods

In this paper we present efficient quadrature rules for the numerical approximation of integrals of polynomial functions over general polygonal/polyhedral elements that do not require an explicit construction of a sub-tessellation into triangular/tetrahedral elements. The method is based on successive application of Stokes’ theorem; thereby, the underlying integral may be evaluated using only the values of the integrand and its derivatives at the vertices of the polytopic domain, and hence leads to an exact cubature rule whose quadrature points are the vertices of the polytope. We demonstrate the capabilities of the proposed approach by efficiently computing the stiffness and mass matrices arising from hp-version symmetric interior penalty discontinuous Galerkin discretizations of second-order elliptic partial differential equations.


Introduction
In recent years the exploitation of computational meshes composed of polygonal and polyhedral elements has become very popular in the field of numerical methods for partial differential equations.Indeed, the flexibility offered by polygonal/polyhedral elements allows for the design of efficient computational grids when the underlying problem is characterized by a strong complexity of the physical domain, such as, for example, in geophysical applications, fluid-structure interaction, or crack propagation problems.Moreover, the possibility to adopt computational meshes with hanging nodes is included in this framework by observing that, for example, a classical quadrilateral element with a hanging node on one of its edges can be treated as a pentagon with two aligned edges.Several conforming numerical discretization methods which admit polygonal/polyhedral meshes have been proposed within the current literature; here, we mention, for example, the Composite Finite Element Method [5,41,42], the Mimetic Finite Difference (MFD) method [4,[18][19][20][21]44], the Polygonal Finite Element Method [63], the Extended Finite Element Method [37,64], the Virtual Element Method (VEM) [10,11,[15][16][17] and the Hybrid High-Order (HHO) method [33][34][35].In the non-conforming setting, we mention Discontinuous Galerkin (DG) methods [1][2][3]6,9,14,[22][23][24][25]25], Hybridizable DG methods [29][30][31][32], non-conforming VEM [8,13,26], and the Gradient Schemes [36]; here the possibility of defining local polynomial discrete spaces follows naturally with the flexibility provided by polytopic meshes.
One of the key aspects concerning the development of efficient finite element discretizations with polygonal/polyhedral grids is the construction of quadrature formulae for the approximate computation of the terms appearing in the underlying weak formulation.Indeed, the design of efficient quadrature rules for the numerical computation of integrals over general shaped polytopes is far from being a trivial task.The classical and most widely employed technique for the integration over polytopes is the Sub-Tessellation method, cf.[38,51,62]; here, the domain of integration is subdivided into standard-shaped elements, such as triangular/quadrilateral elements in 2D or tetrahedral/hexahedral elements in 3D, whereby standard efficient quadrature rules are employed, cf.[50,60,70], and also [71] and [48], for an interpolation technique based on the same idea.On the one hand this technique is easy to implement, however, it is generally computationally expensive, particularly for high order polynomials, since the number of function evaluations may be very large.
For this reason, the development of quadrature rules that avoid sub-tessellation is an active research field.Several approaches have been proposed; in particular, we mention [43,54,67,68], for example.One interesting method in this direction is represented by the Moment Fitting Equation technique, firstly proposed by Lyness and Monegato in [49], for the construction of quadrature rules on polygons featuring the same symmetry as the regular hexagon.Generalizations to convex and non-convex polygons and polyhedra were proposed by Mousavi et al. in [53].Here, starting from an initial quadrature rule, given, for example, by the sub-tessellation method described above, an iterative node elimination algorithm is performed based on employing the least-squares Newton method [69] in order to minimise the number of quadrature points while retaining exact integration.Further improvements of the moment fitting equation algorithm can also be found in [52] and [61].While this method is optimal with respect to the number of function evaluations, the nodes and weights must be stored for every polygon, thus affecting memory efficiency.An alternative approach designed to overcome the limitations of the sub-tessellation approach is based on employing the generalized version of Stokes' theorem; here, the exploitation of Stokes' theorem reduces the integral over a polytope to an integration over its boundary; see [66] for details.For the two-dimensional case, in [59], Sommariva and Vianello proposed a quadrature rule based on employing Green's theorem.In particular, if an x-or y-primitive of the integrand is available (as for bivariate polynomial functions), the integral over the polygon is reduced to a sum of line integrations over its edges.When the primitive is not known, this method does not directly require a sub-tessellation of the polygon, but a careful choice of the parameters in the proposed formula leads to a cubature rule that can be viewed as a particular sub-tessellation of the polygon itself.However, it is not possible to guarantee that all of the quadrature points lie inside the domain of integration.An alternative and very efficient formula has been proposed by Lasserre in [47] for the integration of homogeneous functions over convex polytopes.This technique has been recently extended to general convex and non-convex polytopes in [27].The essential idea here is to exploit the generalized Stokes' theorem together with Euler's homogeneous function theorem, cf.[58], in order to reduce the integration over a polytope only to boundary evaluations.The main difference with respect to the work presented in [59] is the possibility to apply the same idea recursively, leading to a quadrature formula which exactly evaluates integrals over a polygon/polyhedron by employing only point-evaluations of the integrand and its derivatives at the vertices of the polytope.
In this article we extend the approach of [27] to the efficient computation of the volume/face integral terms appearing in the discrete weak formulation of second-order elliptic problems, discretized by means of high-order DG methods.We point out that our approach is completely general and can be directly applied to other discretization schemes, such as VEM, HHO, Hybridisable DG, and MFD, for example.We focus on the DG approach presented in [25], where the local polynomial discrete spaces are defined based on employing the bounded box technique [39].We show that our integration approach leads to a considerable improvement in the performance compared to classical quadrature algorithms based on sub-tessellation, in both two-and three-dimensions.The outline of this article is as follows: in Sect. 2 we recall the work introduced in [27], and outline how this approach can be utilized to efficiently compute the integral of d-variate polynomial functions over general polytopes.In Sect. 3 we introduce the interior penalty DG formulation for the numerical approximation of a secondorder diffusion-reaction equation on general polytopic meshes.In Sect. 4 we outline the exploitation of the method presented in Sect. 2 for the assembly of the mass and stiffness matrices appearing in the DG formulation.Several two-and three-dimensional numerical results are presented in Sect. 5 in order to show the efficiency of the proposed approach.Finally, in Sect.6 we summarise the work undertaken in this article and discuss future extensions.

Integrating Polynomials over General Polygons/Polyhedra
In this section we review the procedure introduced by Chin et al. in [27] for the integration of homogeneous functions over a polytopic domain.To this end, we consider the numerical computation of P g(x)dx, where Each face F i lies in a hyperplane H i identified by a vector a i ∈ R d and a scalar number b i , such that We observe that a i , i = 1, . . ., m, can be chosen as the unit outward normal vector to F i , i = 1, . . ., m, respectively, relative to P, cf.Figs. 1 and 2.
Fig. 1 Example of a two-dimensional polytope P and its face F i .The hyperplane H i is defined by the local origin x 0,i and the vector e i1 Fig. 2 The dodecahedron P with pentagonal faces and the face F i ⊂ ∂P with unit outward normal vector n i .Here, F i has five edges F i j , j = 1, . . ., 5, and five unit outward normal vectors n i j , j = 1, . . ., 5, lying on the plane H i .The hyperplane H i is identified by the local origin x 0,i and the orthonormal vectors e i1 , e i2 g : P → R is a homogeneous function of degree q ∈ R, i.e., for all λ > 0, g(λx) = λ q g(x) for all x ∈ P.
We recall that Euler's homogeneous function theorem [58] states that, if g is a homogeneous function of degree q ≥ 0, then the following identity holds: Next we introduce the generalized Stokes' theorem, which can be stated as follows, cf.[66]: given a generic vector field X : P → R d , the following identity holds where n is the unit outward normal vector to P and dσ denotes the (d − 1)-dimensional (surface) measure.Selecting X = x in (3), and employing (2), gives Equation ( 4) states that if g is homogeneous, then the integral of g over a polytope P can be evaluated by computing the integral of the same function over the boundary faces then the following identity holds: ( Before outlining the details regarding the recursive application of the Stokes' Theorem to (4), we first require the following lemma.
Lemma 1 Let F i j ⊂ ∂F i j = 1, . . ., m i , be the vertices/edges of F i , i = 1, . . ., m, for d = 2, 3, respectively, and let n i j be the unit outward normal vectors to F i j lying in H i .Moreover, let F i j ⊂ ∂ F i be the preimage of F i j with respect to the map γ , and n i j be the corresponding unit outward normal vector.Then, the following holds 1) , whose columns are the vectors {e in } d−1 n=1 , i = 1, . . ., m. Proof We first note that employing the definition of γ we have that The proof now follows immediately from simple linear algebra considerations; for full details, we refer to [7].
Given identity (5) and Lemma 1, we can prove the following result.
Proposition 1 Let F i , i = 1, . . ., m, be a face of the polytope P, and let F i j , j = 1, . . ., m i , be the planar/straight faces/edges such that ∂F i = ∪ m i j=1 F i j for some m i ∈ N.Then, for any homogeneous function g, of degree q ≥ 0, the following identity holds where d i j denotes the Euclidean distance between F i j and x 0,i , x 0,i ∈ H i , is arbitrary, i = 1, . . ., m, and dν denotes the (d − 2)-dimensional (surface) measure.

Proof If we denote by
the gradient operator on the hyperplane H i , i = 1, . . ., m, with respect to the coordinate system ( x1 , . . ., xd−1 ), then, upon application of Stokes' theorem, we have where ñ is the unit outward normal vector of F i and X is a vector field on R d−1 .Next, we transform (7) back to the original coordinate system.To this end, denoting E ∈ R d×(d−1) to be the matrix whose columns are the vectors {e in } d−1 n=1 , we observe that, if we choose X = x, then its divergence is ∇ i • X = d − 1. Exploiting (6), the term ∇ i g(γ (x)) can be written as follows: Exploiting ( 6) and ( 8), we can write 1 and 2 as 1 = (d − 1) respectively.Employing Lemma 1, together with (6), we have that We observe that the term (x − x 0,i ) • n i j is constant for any x ∈ F i j , and that it represents the Euclidean distance between F i j and x 0,i ; thereby, we define d i j = (x − x 0,i ) • n i j .From the above identities ( 9), (10) and (11) we deduce the statement of the Proposition.
Using Proposition 1, together with Eq. ( 4), we obtain the following identity (12) where we recall that ∂P = ∪ m i=1 F i and Remark 1 If d = 2, then F i j is a point and (12) states that the integral of g on P can be computed by vertex-evaluations of the integrand plus a line integration of the partial derivative of g.If d = 3 we can apply Stokes' Theorem recursively to F i j g(x)dν.Proceeding as before, we get where ∂F i j = ∪ m i j k=1 F i jk , x 0,i j is an arbitrarily chosen origin for F i j , and d i jk is the Euclidean distance between F i jk and x 0,i j .
In view of the application of Proposition 1 to finite element methods, we are interested in the integration of a particular class of homogeneous functions, namely polynomial homogeneous functions of the form where k n ∈ N 0 for n = 1, . . ., d.
In this case, g is a homogeneous function of degree q = k 1 + • • • + k d , and the general partial derivative ∂ g ∂ x n is a homogeneous function of degree q − 1.With this in mind, it is possible to recursively apply formula (12) to the terms involving the integration of the derivatives of g.To this end, we write m, will be an edge or a face, respectively; see Table 1 for details.We define the function end if 123 which returns the integral of the polynomial x k 1 1 . . .x k d d over E, where dσ N is the Ndimensional (surface) measure, N = 1, 2, . . ., d.According to Proposition 1, the recursive definition of the function Remark 2 With a slight abuse of notation, when 1 ≤ N ≤ d − 1, in Algorithm 1 (and for the purposes of the following discussion), the point x 0 = (x 0,1 , . . ., x 0,d ) denotes an arbitrarily chosen origin for the coordinate system which defines the N -polytope E and d i represents the Euclidean distance between the (N − 1)-polytopes E i , which form the boundary of E, and x 0 , i = 1, . . ., m.Furthermore, in Algorithm 1, b i , i = 1, . . ., m, is the same constant appearing in (1).Here it can be evaluated as b i = n i • v, where v is a vertex of E i and n i is the unit outward normal vector, i = 1, . . ., m.

Remark 3
We point out that in (12), cf. also (13), the shape of the underlying polytope can be general: indeed, nonconvex simply-connected domains E are admissable.

Integration of Bivariate Polynomials over Polygonal Domains
In order to test the performance of the method proposed in Algorithm 1, we consider the integration of bivariate homogeneous functions on a given polygon P ⊂ R 2 based on using the three different approaches: A. 1 Recursive algorithm described in Sect.2, based on the formula (13): P x k y l dx = I(2, P, k, l), cf.Algorithm 1. A.2 Use of the formula (4) together with numerical integration employed for the evaluation of the edge integrals with known one-dimensional Gaussian quadrature rules, as recently proposed in [28]; A. 3 Sub-tessellation technique: the domain of integration P is firstly decomposed into triangles where standard efficient quadrature rules are then employed.
We test the three different approaches for integrating bivariate polynomials of different polynomial degrees on the triangle depicted in Fig. 3 and the two irregular polygons shown in Figs. 4 and 5, cf.Table 2 for the list of coordinates for each domain; the actual values of the integrals are given in Table 3.In Table 4 we show the average CPU-time taken to evaluate the underlying integral using each method.We point out that, for each integrand and each integration domain P, the relative errors between the output of the three different approaches are of the order of machine precision; that is, all three algorithms return the exact integral up to roundoff error.For completeness, we note that the times for A.1 include the computation of b i , n i , and d i j , j = 1, . . ., m i , i = 1, . . ., m.For A.2 we take into account the evaluation of b i , n i , i = 1, . . ., m, and the one-time computation of the one-dimensional quadrature defined on (−1, 1), consisting of N nodes and weights, employed for the line integrations.
Here, we select N = k+l 2 +1, in order to guarantee the exact integration of x k y l .The times for A.3 include the one-time computation of the N 2 nodes and weights on the reference tri-  angle, where N is selected as in A.2, the time required for sub-tessellation, as well as the time needed for numerical integration on each sub-triangle.The results shown in Table 4 illustrate that the sub-tessellation approach A.3 is the slowest while the proposed method A.1 is the fastest for all of the considered cases; in particular, we highlight that, even for just a single domain of integration, the former method is between one-to two-orders of magnitude slower than the latter approach proposed in this article.Moreover, when the integration domain consists of a triangle, our algorithm A.1 still outperforms classical quadrature rules, cf.A.3, even though in this case no sub-tessellation is undertaken.When comparing A.1 and A.2, we observe that the former algorithm is again superior in terms of CPU time in comparison with the latter approach; this difference seems to grow when the exponents k and l of the integrand function x k y l are very different.This is because in A. 1 we have made an optimal selection of the points x 0,i = (x 0i,1 , x 0i,2 ) , i = 1, . . ., m, appearing in (12).Indeed, performing the geometric reduction of the edges of the domain of integration, we then choose x 0i,1 = 0 or x 0i,2 = 0, i = 1, . . ., m, if the exponents of the integrand function x k y l are k ≥ l or k < l, respectively.The choice x 0i,1 = 0 or x 0i,2 = 0, i = 1, . . ., m, allows us to avoid the recursive calls to the function I(•, •, . . ., •) related to the x-or y-partial derivatives, respectively.In this way the approach A.1 is able to exploit the form of the integrand in order to optimize the evaluation of the corresponding integral.To explore this issue further, in the following section we consider the computational complexity of A.1 in both the cases when an optimal and non-optimal selection of the points x 0,i , i = 1, . . ., m, is made.

Computational Complexity of Algorithm 1
The computational complexity of Algorithm 1, which is employed in A.1, depends in general on the number of recursive calls of the function I(•, •, . . ., •).In particular, using the shorthand notation introduced in Remark 2, the selection of the points x 0 = (x 0,1 , . . ., x 0,d ) , which are used to define the origin of the coordinate system of each N -polytope E which defines the facets of P is crucial.In general, any (d − 1)-dimensional hyperplane in R d possesses a non-empty intersection with some axis of the Cartesian reference system, which means that it is always possible to choose (d − 1) components of x 0 as zero.Without loss of generality we select x 0,r = b i/n i,r and x 0,s = 0 for s = r , where b i and n i are as defined in Remark 2, and r ∈ {1, . . ., d} is chosen so that k r = min{k 1 , . . ., k d }. -Set d = 2 and assume k 1 ≤ k 2 , so that we can choose x 0,1 = 0 and x 0,2 = 0 on each of the edges of P.Then, according to Algorithm 1 we have where we have denoted the vertices of the edge E i as v i1 and v i2 .Hence, In general, for d = 2 we deduce that -Set d = 3 and assume k 1 = min{k 1 , k 2 , k 3 }, so that we may select x 0,1 = 0 and x 0,2 = x 0,3 = 0 on each of the faces of P. Thereby, employing Algorithm 1 we deduce that where, for each i = 1, . . ., m, Here, the computational complexity of I(1, E i j , k, k 2 , k 3 ) depends on the choice of x 0 ≡ x 0,i j which defines the origin of the coordinate system for E i j , j = 1, . . ., m i , i = 1, . . ., m.According to Remark 4, two components of x 0,i j can possibly be different from zero, which implies that the complexity of Algorithm 1 increases exponentially when d = 3.However, it is possible to modify Algorithm 1 in order to avoid the double recursive calls which cause this exponential complexity.In particular, in Sect.2.3 we propose an alternative algorithm which exploits the same idea of Algorithm 1 and allows us to overcome this issue.
In order to confirm (14), we use the tool [56] to measure the number of FLOPs required to exactly compute P x k 1 y k 2 dx; moreover, comparisons will also be made with A.3.To simplify the presentation, the polygon P is selected to be the triangle with vertices (−1, 0.3), (1, −1), and (0.3, 1); thereby, A.3 does not require the computation of a sub-tessellation.In Fig. 6, we plot the number of FLOPs needed to evaluate P x k 1 y k 2 dx by fixing k 1 and varying k 2 ∈ {0, . . ., 50} employing both A.1 and A.3.In particular, Fig. 6a shows that the number of FLOPs required by the quadrature free method A.1 growths linearly with respect to k 2 when k 1 > k 2 and becomes constant as k 2 increases when k 1 ≤ k 2 .Figure 7 confirms the asymptotic behaviour of the two algorithms in the case when k 1 = k 2 ; here, the number of FLOPs required by the sub-tessellation method is reported in both the case when the cost of the evaluation of the quadrature nodes and weights employing the function gauleg, cf.[55], for example, is included/excluded.In particular, we show results both in the case when the cost of the function evaluations is excluded, cf.Fig. 7a, as well as the total number of FLOPs required by each algorithm to exactly evaluate P x k y k dx, cf.Fig. 7b.As expected, the computational complexity of A.3 grows as O(k 2 ) and O(k 3 ) in these two latter cases, respectively, while the cost of the quadrature free method is always O(k) as k increases.Thus far, the numerical results presented for the proposed quadrature free method have assumed that the points x 0 , which are used to define the origin of the coordinate system of each N -polytope E which defines the facets of P, has been chosen in an optimal manner to ensure that the number of recursive calls of I(•, •, . . ., •), cf.Algorithm 1, is minimized.Indeed, a sub-optimal choice of these points leads to an exponential growth in the number of recursive calls of the function I(•, •, . . ., •) in Algorithm 1.For example, if d = 2 the non-optimal choice of x 0 implies that each call of I(•, •, . . ., •) with N = 1 leads to a double recursive call of I(•, •, . . ., •), up to when a zero exponent k 1 or k 2 appears as input.In particular, if k 1 = k 2 = k, it is possible to show that the number of FLOPs required by the quadrature free method grows as O(2 2k−1 ), as k increases, cf.Fig. 8.In the following section, we present an alternative implementation of the quadrature free algorithm which avoids this exponential growth, irrespective of the selection of the points x 0 .

Integration of Families of Monomial Functions
In the context of employing the quadrature free approach within a finite element method, in practice we are not interested in integrating a single monomial function, but instead an entire family of monomials, which, for example, form a basis for the space of polynomials of a given degree over a given polytopic element κ which belongs to the underlying computational mesh.For example, when d = 2, let us consider the evaluation of We note that even when employing the Approach A.1 with an optimal choice of the points x 0 , the total number of FLOPs required for the computation of ( 15) is approximately O( p 3 ), as p increases.
To improve the dependence on p we propose an alternative approach, cf.Algorithm 2; this is based on the observation that, using the notation of Algorithm 1, if the values of in Algorithm 1, have already been computed, then the computation of I(N , E, k 1 , . . ., k d ) is extremely cheap.Indeed, since we must store the integrals of all the monomials on κ anyway, we can start by computing and storing κ x k 1 y k 2 dx 1 dx 2 related to the lower degrees k 1 , k 2 and N = 1, then exploit these values in order to compute the integrals with higher degrees k 1 , k 2 and higher dimension N of the integration domain E.This leads to an algorithm, whereby the number of FLOPs required to compute and store ), as p increases, irrespective of the selection of choice of the points x 0 .In Fig. 9 we now compare these two approaches for d = 2, when the underlying element is selected to be the triangular region employed in the previous section.Here, we compare Algorithm 1, with an end for end for end procedure optimal selection of the points x 0 , with Algorithm 2, where in the latter case the points x 0 are simply selected to be equal to the first vertex defining each edge; here, we clearly observe the predicted increase in FLOPs of O( p 3 ) and O( p 2 ), as p increases, for each of the two algorithms, respectively.

Application to hp-Version DG Methods
We consider the following elliptic model problem, given by: find u such that where Ω ⊂ R d , d = 2, 3, is a polygonal/polyhedral domain with boundary ∂Ω and f is a given function in L 2 (Ω).
In order to discretize problem ( 16), we introduce a partition T h of the domain Ω, which consists of disjoint (possibly non-convex) open polygonal/polyhedral elements κ of diameter h κ , such that Ω = κ∈T h κ.We denote the mesh size of T h by h = max κ∈T h h κ .Furthermore, we define the faces of the mesh T h as the planar/straight intersections of the (d − 1)-dimensional facets of neighbouring elements.This implies that, for d = 2, a face consists of a line segment, while for d = 3, the faces of T h are general shaped polygons; without loss of generality, for the definition of the proceeding DG method we assume that the faces are (d − 1)-dimensional simplices, cf.[24,25] for a discussion of this issue.In order to introduce the DG formulation, it is helpful to distinguish between boundary and interior element faces, denoted by F B h and F I h , respectively.In particular, we observe that F ⊂ ∂Ω for F ∈ F B h , while for any F ∈ F I h we assume that F ⊂ ∂κ ± , where κ ± are two adjacent elements in T h .Furthermore, we write F h = F I h ∪ F B h to denote the set of all mesh faces of T h .For simplicity of presentation we assume that each element κ ∈ T h possesses a uniformly bounded number of faces under mesh refinement, cf.[24,25].
In order to define the DG method, we introduce the jump and average operators: where v ± and τ ± denote the traces of sufficiently smooth scalar-and vector-valued functions v and τ , respectively, on F taken from the interior of κ ± , respectively, and n ± are the unit outward normal vectors to ∂κ ± , respectively, cf.[12].
We then consider the bilinear form

corresponding to the symmetric interior penalty DG method, defined by
where ∇ h denotes the broken gradient operator, defined elementwise, and α h ∈ L ∞ (F h ) denotes the interior penalty stabilization function, whose precise definition, based on the analysis introduced in [24,25], is given below.To this end, we first need the following definition.
Definition 1 Let Th be the subset of elements κ ∈ T h such that each κ ∈ Th can be covered by at most n T shape-regular simplices , and |K i | ≥ c as |κ| for all i = 1, . . ., n T , for some n T ∈ N, where C as and c as are positive constants, independent of κ and T h .
Lemma 2 Let κ ∈ T h , F ⊂ ∂κ denote one of its faces, and Th be defined as in Definition 1.
Then, for each v ∈ P p κ (κ), we have the inverse estimate where and κ F denotes a d-dimensional simplex contained in κ which shares the face F with κ ∈ T h .Furthermore, C inv is a positive constant, which if κ ∈ Th depends on the shape regularity of the covering of κ given in Definition 1, but is always independent of |κ|/ sup κ F ⊂κ |κ F |, p κ and v.
Based on Lemma 2, together with the analysis presented in [24,25], the parameter α h can be defined as follows.
Definition 2 Let α h : F h → R + be defined facewise by where C min α is a sufficiently large lower bound.
The DG discretization of the problem ( 16) is given by: find u h ∈ V h such that By fixing a basis {φ i } N h i=1 , N h denoting the dimension of the discrete space V h , ( 19) can be rewritten as: find U ∈ R N h such that where f i = Ω f φ i dx ∀i = 1, . . ., N h , A is the stiffness matrix, given by A i j = A h (φ j , φ i ) ∀i, j = 1, . . ., N h , M is the mass matrix, and U contains the expansion coefficients of u h ∈ V h with respect to the chosen basis.In order to assemble (A + M) we need to compute the following matrices: for i, j = 1, . . ., N h , where as before N h denotes the dimension of the DG space V h .In particular, the stiffness matrix related to the DG approximation of problem ( 19) is defined as A = V − I − I + S.

Shape Functions for the Discrete Space V h
To construct the discrete space V h we exploit the approach presented in [25], based on employing polynomial spaces defined over the bounding box of each element.More precisely, given an element κ ∈ T h , we first construct the Cartesian bounding box B κ , such that κ ⊂ B κ .Given B κ , κ ∈ T h , it is easy to define a linear map between B κ and the reference element B = (−1, 1) d as follow: where J κ ∈ R d×d is the Jacobi matrix of the transformation which describes the stretching in each direction, and t κ ∈ R d is the translation between the point 0 ∈ B and the baricenter of the bounded box B κ , see Fig. 10.We note that since F κ affinely maps one bounding box to another (without rotation), the Jacobi matrix J κ is diagonal.
Employing the map F κ , κ ∈ T h , we may define a standard polynomial space P p (B κ ) on B κ spanned by a set of basis functions {φ i,κ } for i = 1, . . ., N p κ = dim(P p (B κ )).More precisely, we denote by {L n (x)} ∞ n=0 the family of one-dimensional and L 2 -orthonormal Legendre polynomials, defined over L 2 (−1, 1), i.e., cf. [40,57].We then define the basis functions for the polynomial space P p ( B) as follows: writing I = (i 1 , i 2 , . . ., i d ) to denote the multi-index used to identify each basis function Then, the basis functions for the polynomial space P p κ (κ) are defined by using the map F κ , namely: The set {φ I ,κ : 0 ≤ |I | ≤ p κ , κ ∈ T h } forms a basis for the space V h .On each element κ ∈ T h we introduce a bijective relation between the set of multi-indices {I = (i 1 , . . ., i d ) : 0 ≤ |I | ≤ p κ } and the set {1, 2, . . ., N p κ }.

Volume Integrals Over Polytopic Mesh Elements
In the following we describe the application of Algorithm 2 to compute the entries in the local mass and element-based stiffness matrices respectively, for all κ ∈ T h .For simplicity of presentation, we restrict ourselves to twodimensions, though we emphasize that the three-dimensional case is analogous, cf.Sect.5.2 below.Since the basis functions are supported only on one element, employing the transformation F κ , we have where in the last integral κ = F −1 κ (κ) ⊂ B, see Fig. 10.Here, the Jacobian of the transformation F κ is given by |J κ | = (J κ ) 1,1 (J κ ) 2,2 , which is constant, due to the definition of the map.In order to employ the homogeneous function integration method described in the previous section, we need to identify the coefficients of the homogeneous polynomial expansion for the function φi ( x, ŷ) φ j ( x, ŷ).We observe that φi ( x, ŷ) = L i 1 ( x)L i 2 ( ŷ), and each one-dimensional Legendre polynomial can be expanded as Therefore, we have Here, we have written Notice that the coefficients C i, j,k can be evaluated, once and for all, independently of the polygonal element κ.We now consider the general element of the volume matrix V i, j , cf. (24).Proceeding as before, let I , J be the two multi-indices corresponding respectively to i and j, we have Proceeding as before, we apply a change of variables to the terms 1 and 2 with respect to the map F κ ; thereby, we obtain From the definition of F κ , the inverse map is given by Then, using the definition (23) of the basis functions, we have the following characterization of the partial derivatives appearing in the terms 1 and 2 : where we have used that Then, 1 can be written as: , the integrand function of term 1 is a polynomial.Thereby, we have the following relation: From the expansion (25) of the Legendre polynomials, we note that xm , for i > 0; (28) where the indices C i,m = (m + 1)C i,m+1 are the coefficients for the expansion of L i (•).We deduce that 1 = 0 if i 1 = 0 or j 1 = 0, and where C i 2 , j 2 ,l is defined in (26), and with C i,n = (n+1)C i,n+1 , C j,m = (m+1)C j,m+1 , cf. (28), is the expansion of the derivatives of the Legendre polynomials which is computable independently of the element κ, κ ∈ T h .Analogously, we deduce the following expression for the second term of Eq. ( 27):

Interface Integrals over Polytopic Mesh Elements
With regards the interface integrals appearing in Eq. ( 18), we describe the method by expanding the jump and average operators and computing each term separately, working, for simplicity, again in two-dimensions.Firstly, we discuss how to transform the integral over a physical face F ⊂ ∂κ to the corresponding integral over the face F = F −1 κ (F) ⊂ ∂ κ on the reference rectangular element κ.To this end, let F ⊂ ∂κ be a face of the polygon κ, κ ∈ T h , and let x 1 = (x 1 , y 1 ) and x 2 = (x 2 , y 2 ) denote the vertices of the face, based on counter clock-wise ordering of the polygon vertices.The face F = F −1 κ (F) is identified by the two vertices x1 = F −1 κ (x) and x2 = F −1 κ (x 2 ).For a general integrable function g : κ → R we have where dσ (F κ ( x, ŷ)) = J F d σ and J F is defined as where n F is the unit outward normal vector to F.
We next describe how to compute the interface integrals.From the definition of the jump and average operators, cf.(17), on each edge F ∈ F I h shared by the elements κ ± we need to assemble for i, j = 1, . . ., N p κ ± .Analogously, on the boundary face F ∈ F B h belonging to κ + ∈ T h we only have to compute for i, j = 1, . . ., N p κ + .We next show how to efficiently compute a term of the form where I , J are the suitable multi-indices associated to i, j = 1, . . ., N p κ + , respectively.Proceeding as before, we have Analogously, we have For the term a , we directly apply the definition of the basis function, and obtain while for the term b we have In order to obtain a homogeneous polynomial expansion for b we have to write explicitly the composite map where the matrix J is diagonal since J −1 κ − and J κ + are diagonal.We then have b = φJ ( F(x)) = φJ ( Jx + t) = φJ ( J1,1 x + t1 , J2,2 ŷ + t2 ) Combining (30) and (31), and denoting by F+ = F −1 κ + (F), cf.Fig. 11, from (29) we obtain where X and Ỹ are defined as Here, as before, C i,n are the coefficients of the homogeneous function expansion of the Legendre polynomials in (−1, 1), while X j,m and Ỹ j,m are defined by here, we have exploited the Newton-binomial expansion of the terms ( J1,1 x + t1 ) k and ( J2,2 ŷ + t2 ) l appearing in equation (31).Similar considerations allow us to compute where C i, j,k are defined as and where n + = [n + x , n + y ] is the unit outward normal vector to the physical face F from κ + .Similarly, where we have also introduced X and Ỹ defined as

Remark 5
The coefficients X and Ỹ depend on the maps F κ + and F κ − , as well as X , X , Ỹ and Ỹ ; thereby, they must be computed for each element κ in the mesh T h .

Remark 6
With regards the computation of the forcing term we point out that the quadrature method proposed in this paper allows to exactly evaluate (32) when f is a constant or a polynomial function.If f is a general function, an explicit polynomial approximation of f is required.

Numerical Experiments
We present some two-and three-dimensional numerical experiments to test the practical performance of the proposed approach.Here, the results are compared with standard assembly algorithms based on employing efficient quadrature rules on a sub-tessellation.

Two-Dimensional Test Case
We test the performance of the algorithm outlined in Sect. 4 for the computation of the elemental mass and stiffness matrices resulting from the DG discretization (19) on Voronoi decompositions as shown in Fig. 12.In particular, we compare the CPU-time needed to assemble the local and global elemental matrices using Algorithm 2, cf.Sect.4, with Quadrature Integration over polygonal domains, based on the sub-tessellation method on polygons and Gaussian line integration for the related interface terms.More precisely, given κ ∈ T h , the sub-tessellation scheme on κ is performed by constructing a non-overlapping sub-tessellation κ S = {τ κ } consisting of standard triangular elements; in particular, as, for our tests, we consider Voronoi numerical grids, we exploit the convexity of κ and define κ S by connecting the centre of mass of κ with its vertices.As an example, if we consider computing the elemental mass matrix M κ i, j , we have that where F τ κ : τ → τ κ is the mapping from the reference simplex τ to τ κ , with Jacobian |J τ κ |, and {(ξ r , ω r )} q κ r =1 denotes the quadrature rule defined on τ .The construction of quadrature rules on τ may be computed based on employing the Duffy transformation, whereby the reference tensor-product element (−1, 1) 2 is mapped to the reference simplex.As the algorithm outlined in Sect. 4 does not require the definition of quadrature nodes and weights, in the following we will refer to it as the Quadrature Free Method.Consider the problem (19) introduced in Sect. 3 with d = 2 and Ω = (0, 1) 2 , where we select the set of basis functions {φ i } N h i=1 for V h as described in Sect. 4. In order to quantify the performance of the proposed approach, we consider a series of numerical tests obtained by varying the polynomial degree p κ = p for all κ ∈ T h , between 1 and 6 and by employing a series of uniform polygonal meshes of different granularity, cf.Fig. 12.The numerical grids are constructed based on employing PolyMesher, cf.[65].Here, we are interested in the CPU time needed to assemble the matrices ( 21) and (22).
In the first test case, we consider the CPU time needed to assemble the matrices M and V.As pointed out in Sect.4, these matrices are block diagonal and each block consists of an integral over each polygonal element κ ∈ T h .In Fig. 13 we present the comparison between the CPU times needed to assemble the global matrices M and V based on employing the quadrature free method and quadrature integration (based on sub-tessellation) when varying the number of elements N e ∈ {64, 256, 1024, 4094, 16384, 65536} and the polynomial degree p ∈ {1, 2, 3} (left), and p ∈ {4, 5, 6} (right).Clearly, our approach outperforms the classical sub-tessellation method leading to substantial gains in efficiency.For a more detailed comparison, we have presented in Fig. 15a the logarithmic-scaled graphs of each computation: from the results of Fig. 15a we observe that the CPU time grows like O(N e ), as N e increases, as expected.is more expensive than the time it takes to assemble the face-based matrices S and I, when the classical Gaussian line integration method is employed.This is, of course, due to the greater number of function evaluations required to compute M and V on the underlying subtessellation; note that in two-dimensions, a sub-tesellation of the faces is not necessary, since they simply consist of line segments.However, the opposite behaviour is observed when the quadrature free method is employed; in this case, the volume integrals can be very efficiently computed since the coefficients C i, j,k and C i, j,k only need to be computed once, cf.Sect.4.2.On the other hand, computing the face integrals present in S and I requires the evaluation of the coefficients X j,m , Xi, j,k , X i, j,k , Ỹ j,m , Ỹi, j,k , and Ỹ i, j,k , cf.Sect.4.3, which must be computed for each face F ∈ F h .

Three-Dimensional Test Case
We now consider the diffusion-reaction problem (19) with d = 3 and Ω = (0, 1) 3 .The polyhedral grids employed for this test case are defined by agglomeration: starting from a fine partition T f ine of Ω consisting of N f ine disjoint tetrahedrons {κ i f } N f ine i=1 , such that Ω = ∪ N f ine i=1 κ i f , a coarse mesh T h of Ω consisting of disjoint polyhedral elements κ can be defined such that where S κ ⊂ T f ine h denotes the set of fine elements which forms κ.Here, the agglomeration of fine tetrahedral elements is performed based on employing the METIS library for graph partitioning, cf., for example, [45,46].With this definition each polyhedral element is typically non-convex.For simplicity, we have considered only the case of simply connected elements.In this particular case, the faces of the mesh T h are the triangular intersections of two-dimensional facets of neighbouring elements.Figure 16 shows three examples of the polyhedral elements resulting from agglomeration.
We perform a similar set of experiments as the ones outlined in Sect.5.2 for the twodimensional case.Again, we compare the CPU time required by the proposed quadrature free method with the quadrature integration/sub-tessellation approach to assemble the stiff- sibly polygons of arbitrary shape, a sub-tessellation is needed also for surface integrals.The proposed technique is completely general and can be extended to several numerical methods based on discrete spaces defined on polygonal/polyhedral meshes, such as Virtual Element Methods, Mimetic Finite Differences, Hybrid High-Order Methods, Hybridizable Discontinuous Galerkin Schemes, and Polygonal Finite Element Methods, for example.We stress that for moderate polynomial degrees, the proposed integration technique, which involves exact integration of bivariate and trivariate functions in two-and three-dimensions, respectively, has been observed to be numerically stable.

Fig. 3 Triangle (P 1 )Fig. 4 Fig. 5
Fig. 3 Triangle (P 1 ) be selected to be zero.In this way, the selection of k r essentially fixes the number of recursive calls of I(•, •, . . ., •) in Algorithm 1.More precisley, we write |I(N , E, k 1 , . . ., k d )| Fl to denote the number of FLOPs to perform I(N , E, k 1 , . . ., k d ), and let C N be the number of FLOPs required by I(N , E, k 1 , . . ., k d ), without considering the recursive calls of I(•, •, . . ., •) to itself.With this in mind, let us consider the following two examples:

Fig. 6 Fig. 7
Fig.6 Comparison of the number of FLOPs required to evaluate P x k 1 y k 2 dx, based on fixing k 1 and varying k 2 ∈ {0, . . ., 50}: a Quadrature free method A.1; b Sub-tessellation method A.3

Fig. 8
Fig.8 Number of FLOPs required to evaluate P x k y k dx as k increases, based on employing the quadrature free method, with a sub-optimal choice of x 0

Fig. 9
Fig. 9 Number of FLOPs required to evaluate{ κ x k 1 y k 2 dx ∀ k 1 , k 2 ≥ 0, k 1 + k 2 ≤p} based on employing Algorithm 1 (with an optimal selection of the points x 0 ), and Algorithm 2

Fig. 11
Fig.11Example of a polygonal elements κ ± ∈ T h , together with the bounded boxes B κ ± , and the local maps F κ ± : κ → κ ± for the common face F ⊂ κ ±

Fig. 16
Fig. 16 Example of polyhedral elements κ ∈ T h obtained by agglomeration of tetrahedrons.κ 1 has 18 triangular faces, κ 2 has 20 triangular faces and κ 3 has 22 triangular faces.a Element κ 1 .b Element κ 2 .c Element κ 3 Here, x 0,i ∈ H i , i = 1, ..., m, is an arbitrary point which represents the origin of the coordinate system on H i , and {e in } d−1 n=1 is an orthonormal basis on H i , i = 1, ..., m; see Figs.1 and 2for two-and three-dimensional examples, respectively.Notice that x 0,i does not have to lie inside m.By applying Stokes' theorem recursively, we can further reduce each term end, Stokes' theorem needs to be applied on the hyperplane H i , i = 1, ..., m, in which each F i , i = 1, ..., m, lies, respectively.In order to proceed, let γ : R d−1 → R d be the function which expresses a generic point x = ( x1 , ..., xd−1 ) ∈ R d−1 as a point in R d that lies on H i , i = 1, ..., m, i.e.,x −→ γ (x) = x 0,i + d−1 n=1xn e in , with e in ∈ R d , e in • e im = δ nm .

Table 1
Polytopic domains of integration E considered in Algorithm 1 as a function of the dimension d N

Table 2
Coordinates of the polygons of Figs.3, 4, and 5

Table 3
1he approximated values of the integral over the three polygons in Figs. 3, 4 and 5 obtained with approach A.1

Table 4
3PU times as a function of the integrand and the integration domain P for the three approaches A.1, A.2 and A.3 Algorithm for integrating all monomials up to order p∂E = {E 1 , . .., E m } where E i ⊂ ∂E; F = FaceI ntegrals(d − 1, E 1 , . .., E m , k 1 , . .., k d ); for a 1 = 0 : k 1 , . .., a d = 0 : k d ; k 1 + k 2 + . ..+ k d ≤ p do V (a 1 , . . ., a d ) = FaceI ntegrals(N , E 1 , . . ., E m , k 1 , . . ., k d ) F(−1 : k 1 , . .., −1 : k d , 1 : m) = 0;for i=1:m do choose x 0 as the first vertex of E i ; Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/),which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.