Recursion relations for hp‐FEM Element Matrices on quadrilaterals

In this paper we consider higher order shape functions for finite elements on quadrilaterals. Using tensor products of suitable Jacobi polynomials, it can be proved that the corresponding mass and stiffness matrices are sparse with respect to 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. Besides the H1 case also the conformal basis functions for H(Div) and H(Curl) are investigated.


Formulation of the problem in H 1
For a polygonal Lipschitz domain Ω ∈ R 2 and given f ∈ L 2 (Ω), we consider the model problem − div(D grad u)) + κu = f in Ω, u = 0 on ∂Ω, with an elementwise constant diffusion matrix D ∈ R 2×2 and coefficient κ. For 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 separated into three groups ϕ v , ϕ e and ϕ I using a vertex (v), edge (e), interior (I) ansatz.

Discretization -Shape functions for H 1
For the definition of our ansatz functions, Legendre polynomials and their integrals are required. They are defined by the Rodriguez formula, see e.g. [3]. More precisely, they and their integrals are given by respectively. Let C = (−1, 1) 2 be the reference square, see [4] for the triangular case. The H 1 interior basis functions on C are given by , n, k = 2, . . . , p.
In addition the vertex and edge bubbles are defined accordingly, see [1], [5]. By using a Kronecker product formulation, we only need to evaluate the integrals I (n,m) 1 An extension to the cube is straightforward. Since the H 1 basis functions are given by P n,m,k (x, y, z) = P (0,0) we can again factor the integrals and use a Kronecker product formulation. Thus the recursion formulas stay the same.
In this section, the differential equations −∇∇ · u + κu = f and −∇ × ∇ × u + κu = f for vector valued functions are discretized. These are H(Div) and H(Curl) model problems, respectively. For ease of notation, only the interior bubble functions on the reference element are presented.
P r o o f. The formula can be found using the Mathematica package Guess [7] and proven using the closed form.
For the stiffness matrix one needs to assemble C Div(P (I/II) n1,k1 (x, y)) Div(P (I/II) n1,k1 )(x, y) dx dy. Although the divergence of the basis functions is rather simple, i.e. Div P (I) n,k (x, y) = 0 and Div P (II) n,k (x, y) = L n−1 (x)L k−1 (y) one can still reuse (3) for a fast assembly.

H(Curl) conformal basis functions
The two types of basis functions for H(Curl) on the reference square are given by , n, k = 2, . . . , p.
Again the recursion formulas (3) come in handy. Additionally recursion formula (4) and similar formulas can be used as well.
For the stiffness part, we need to solve integrals of the form C curl(P (I/II) n1,k1 (x, y)) curl(P (I/II) n2,k2 )(x, y) dx dy. Therefore we first compute the curl of the basis functions, i.e. curl P (I) n,k (x, y) = 0 and curl P (II) n,k (x, y) = −L n−1 (x)L k−1 (y). Up to sign, these are the same functions as in the divergence for the H(Div) case. Thus we can compute the H(Curl) in the same way as the H(Div) case.