Symbolic Evaluation of hp‐FEM Element Matrices

In this paper we consider higher order shape functions for finite elements on a triangle. On the reference element Dubiner‐like ansatz functions based on suitable integrated Jacobi polynomials are chosen. It can be proved that the corresponding mass and stiffness matrices are sparse for all polynomial degree p. Due to the orthogonal relations between Jacobi polynomials the exact values of the entries of mass and stiffness matrix can be determined. Using symbolic computation, we can find simple recurrence relations which allow us to compute the remaining nonzero entries in optimal arithmetic complexity.

In this paper we consider higher order shape functions for finite elements on a triangle. On the reference element Dubiner-like ansatz functions based on suitable integrated Jacobi polynomials are chosen. It can be proved that the corresponding mass and stiffness matrices are sparse for all polynomial degree p. Due to the orthogonal relations between Jacobi polynomials the exact values of the entries of mass and stiffness matrix can be determined. Using symbolic computation, we can find simple recurrence relations which allow us to compute the remaining nonzero entries in optimal arithmetic complexity.
with an elementwise constant diffusion matrix D ∈ R 2×2 and coefficient κ. For a simple notation we assume a Dirichlet boundary condition. Problem (1) is discretized by means of the hp-version of finite elements, see e.g. [1]. Denote by Φ = [φ 1 , . . . , φ N ] a basis for the finite element space contained in H 1 0 (Ω). Following [2] the shape functions can be seperated into three groups φ v , φ e and φ I using a vertex (v), edge (e), interior (I) ansatz.

Discretization -Shape functions
For the definition of our ansatz functions, Jacobi polynomials are required. They are defined by the Rodriguez formula, see e.g. [3]. More precisely, let be the n-th Jacobi polynomial with respect to the weight function (1 − x) α (1 + x) 0 . Let△ be the reference triangle with vertices (−1, −1), (1, −1) and (0, 1). Following [4], see also [5], we choose the interior bubbles on△ as i≥2,j≥1 be the vector of interior bubbles. In addition the vertex and edge bubbles are defined accordingly, see [4], [5]. The entries of the interior block of the element mass matrix M and element stiffness matrix K on△ are then defined by It has been shown by [4] that the matrices M and K in (2) are sparse on the reference triangle, i.e.
This result can easily be extended to any non degenerate triangle, since the transformation is affine linear. If D is a diagonal matrix the stronger result K ij,kl = 0 holds for |i − k| / ∈ {0, 2}. Using orthogonality relations of Jacobi polynomials, see e.g. [3], the exact value of the integrals can be computed. This leads to big broken rational functions, see [6]. Alternatively the non zero entries can be computed by using sum factorization in O(p 3 ) operations in 2D, [7]. An alternative approach is to determine recurrence relation between the mass or stiffness matrix entries. Out of simplicity we will consider only the stiffness matrix.

Symbolic Evaluation
The Mathematica package Guess [8] can guess a recurrence relation between different values of broken rational functions, if we provide enough data. Using the algorithms from [6], we compute the exact values of the stiffness matrix entries and use Guess to find a multivariate relation between this entries, which is much simpler than the explicit entries. Using the exact formulas, these relations can be easily verified. It holds that (2). Then all non zero entries can be computed in optimal arithmetical complexity O(p 2 ).

Outlook
This recursion formulas yield the possibility of an efficient way to assemble element matrices of high polynomial degree. Similar recursion formulas to the shown formula can be found for the interior block of the mass matrix as well. By choosing efficient edge bubbles an extension to the other blocks is possible. We refer to [6] for generalizations to the three dimensional case and to [9] for other PDE's like the Maxwell equations. Finally the approach can be used for preconditioning in case of an arbitrary diffusion matrix D.