basix::polyset Namespace Reference

Polynomial expansion sets. More...

Enumerations
enum class	type { standard = 0 , macroedge = 1 }
	Cell type.

Functions
template<std::floating_point T>
std::pair< std::vector< T >, std::array< std::size_t, 3 > >	tabulate (cell::type celltype, polyset::type ptype, int d, int n, MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan< const T, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents< std::size_t, 2 >> x)
	Tabulate the orthonormal polynomial basis, and derivatives, at points on the reference cell. More...
template<std::floating_point T>
void	tabulate (MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan< T, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents< std::size_t, 3 >> P, cell::type celltype, polyset::type ptype, int d, int n, MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan< const T, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents< std::size_t, 2 >> x)
	Tabulate the orthonormal polynomial basis, and derivatives, at points on the reference cell. More...
int	dim (cell::type cell, polyset::type ptype, int d)
	Dimension of a polynomial space. More...
int	nderivs (cell::type cell, int d)
	Number of derivatives that the orthonormal basis will have on the given cell. More...
polyset::type	superset (cell::type cell, polyset::type type1, polyset::type type2)
	Get the polyset types that is a superset of two types on the given cell. More...
polyset::type	restriction (polyset::type ptype, cell::type cell, cell::type restriction_cell)
	Get the polyset type that represents the restrictions of a type on a subentity. More...

Detailed Description

Polynomial expansion sets.

In Basix, the bases of the polynomial sets spanned by finite elements are represented by a set of coefficients of orthonormal polynomials on the cell. By orthogonality, we mean that

\[\int_R p_ip_j=\begin{cases}1&i=j\\0&i\not=j\end{cases}\]

Interval

Legendre polynomials form a set of orthogonal polynomials on an interval. Legrendre polynomials are a special case of the Jacobi polynomials, \(P_n^{\alpha, \beta}\), with the weights \((\alpha, \beta) = (0, 0)\). Using the Jacobi recurrence relation, we have

\[ n P^{0,0}_n(x) = (2n - 1) x P^{0,0}_{n - 1}(x) - (n-1) P^{0,0}_{n-2}(x),\]

starting with \(P^{0,0}_0 = 1\). Legendre polynomials are orthogonal when integrated over the interval [-1, 1].

In Basix, the interval [0, 1] is used as the reference, so we apply the transformation \(x\mapsto 2x-1\) to the standard Legendre polynomials. The orthogonal polynomials we use are therefore given by the recurrence relation

\[n p_n(x)=(2n-1)(2x-1)p_{n-1}(x)-(n-1)p_{n-2}(x).\]

For a given set of points, the recurrence relation can be used to compute the polynomial set up to arbitrary order.

Triangle (Dubiner's basis)

See Sherwin & Karniadakis (1995).

An orthogonal set over a triangle can be obtained with a change of variable to \(\zeta = 2\frac{1+x}{1 - y} - 1\), and a modified polynomial. The polynomial basis becomes:

\[ Q_{p,q} = P^{0,0}_p(\zeta) \left(\frac{1-y}{2}\right)^p\ P_q^{2p+1, 0}(y), \]

with a Legendre Polynomial of \(\zeta\) and a weighted Jacobi Polynomial of \(y\). In order to calculate these without actually multiplying together the two polynomials directly, we can first compute \(Q_{p, 0}\) using a recurrence relation (since \(P_0^{2p + 1, 0} = 1\)):

\[ p Q_{p,0} = (2p-1)(2x + y+1)Q_{p-1,0} - (p-1)(1-y)^2 Q_{p-2, 0}.\]

Subsequently, we can calculate \(Q_{p,q}\) by building up from \(Q_{p,0}\) with another recurrence relation:

\[ Q_{p,q} = (a_0 y + a_1) Q_{p, q-1} + a_2 Q_{p, q-2}, \]

where \(a_0, a_1, a_2\) are the coefficients in the Jacobi recurrence relation with \((\alpha,\beta) = (2p+1, 0)\), see Wikipedia. Note that \(a_2 = 0\) when \(q < 2\).

Tetrahedron

See Sherwin & Karniadakis (1995).

Let \(\zeta = 2\frac{1+x}{y+z} + 1\) and \(\xi = 2\frac{1+y}{1-z} - 1\). Orthogonal polynomials on a tetrahedron are then given by

\[ Q_{p, q, r} = P^{0,0}_p(\zeta)\left(\frac{y+z}{2}\right)^p\ P_q^{2p+1,0}(\xi)\left(\frac{1-z}{2}\right)^q\ P_r^{2(p+q+1), 0}(z).\]

This can similarly be built up by first computing \(Q_{p,0,0}\), then \(Q_{p, q, 0}\) and finally \(Q_{p, q, r}\) with recurrence relations on each axis.

Quadrilateral and hexahedron

The polynomial sets for quadrilateral and hexahedral elements can be formed by applying the recurrence relation for intervals in each direction.

Prism

The polynomial set for a prism can be computed by using using the recurrence relation for the triangle to get the polynomials that are constant with respect to \(z\), then applying the recurrence relation for the interval to compute further polynomials with \(z\)s.

Pyramid

Orthogonal polynomials on the pyramid element are best calculated in the same way as the tetrahedron, using recurrence relations on each axis. Let \(\zeta_x = 2\frac{1+x}{1-z} - 1\), \(\zeta_y = 2\frac{1+y}{1-z} - 1\). The recurrence relation is then

\[Q_{p, q, r} = P^{0,0}_p(\zeta_x) P^{0,0}_q(\zeta_y) \left(\frac{1-z}{2}\right)^{(p+q)} P_r^{2(p+q+1), 0}(z).\]

Normalisation

For each cell type, we obtain an orthonormal set of polynomials by scaling each of the orthogonal polynomials.

Derivatives

Recurrence relations can also be used to find the derivatives of the polynomials at given points. For example, on the interval, the first derivatives are given by

\[ n P'_n(x) = (2n - 1) \left(P_{n-1}(x) + x P'_{n - 1}(x)\right) + (n-1)P'_{n-2}(x). \]

More generally, the \(k\)-th derivative is given by

\[ n P^k_n(x) = (2n - 1) \left(k P^{k-1}_{n-1}(x) + x P^k_{n - 1}(x)\right)+ (n-1) P^k_{n-2}(x) \]

This is now a recurrence relation in both \(n\) and \(k\). Similar recurrence relations can be obtained for the derivatives of all the polynomial sets on the other shapes. Care must be taken with quadratic terms, and cross-terms in two and three dimensions.

Function Documentation

tabulate() [1/2]

template<std::floating_point T>

std::pair< std::vector< T >, std::array< std::size_t, 3 > > basix::polyset::tabulate	(	cell::type	celltype,
		polyset::type	ptype,
		int	d,
		int	n,
		MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan< const T, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents< std::size_t, 2 >>	x
	)

Tabulate the orthonormal polynomial basis, and derivatives, at points on the reference cell.

All derivatives up to the given order are computed. If derivatives are not required, use n = 0. For example, order n = 2 for a 2D cell, will compute the basis \(\psi, d\psi/dx, d\psi/dy, d^2 \psi/dx^2, d^2\psi/dxdy, d^2\psi/dy^2\) in that order (0, 0), (1, 0), (0, 1), (2, 0), (1, 1), (0 ,2).

For an interval cell there are nderiv + 1 derivatives, for a 2D cell, there are (nderiv + 1)(nderiv + 2)/2 derivatives, and in 3D, there are (nderiv + 1)(nderiv + 2)(nderiv + 3)/6. The ordering is 'triangular' with the lower derivatives appearing first.

Parameters

[in]	celltype	Cell type
[in]	ptype	The polynomial type
[in]	d	Polynomial degree
[in]	n	Maximum derivative order. Use n = 0 for the basis only.
[in]	x	Points at which to evaluate the basis. The shape is (number of points, geometric dimension).

Returns

Polynomial sets, for each derivative, tabulated at points. The shape is (number of derivatives computed, number of points, basis index).

• The first index is the derivative. The first entry is the basis itself. Derivatives are stored in triangular (2D) or tetrahedral (3D) ordering, eg if (p, q) denotes p order derivative with respect to x and q order derivative with respect to y, [0] -> (0, 0), [1] -> (1, 0), [2] -> (0, 1), [3] -> (2, 0), [4] -> (1, 1), [5] -> (0, 2), [6] -> (3, 0),... The function basix::indexing::idx maps tuples (p, q, r) to the array index.
• The second index is the point, with index i corresponding to the point in row i of x.
• The third index is the basis function index.

  Todo:

  Does the order for the third index need to be documented?

tabulate() [2/2]

template<std::floating_point T>

void basix::polyset::tabulate	(	MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan< T, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents< std::size_t, 3 >>	P,
		cell::type	celltype,
		polyset::type	ptype,
		int	d,
		int	n,
		MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan< const T, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents< std::size_t, 2 >>	x
	)

Tabulate the orthonormal polynomial basis, and derivatives, at points on the reference cell.

All derivatives up to the given order are computed. If derivatives are not required, use n = 0. For example, order n = 2 for a 2D cell, will compute the basis \(\psi, d\psi/dx, d\psi/dy, d^2 \psi/dx^2, d^2\psi/dxdy, d^2\psi/dy^2\) in that order (0, 0), (1, 0), (0, 1), (2, 0), (1, 1), (0 ,2).

For an interval cell there are nderiv + 1 derivatives, for a 2D cell, there are (nderiv + 1)(nderiv + 2)/2 derivatives, and in 3D, there are (nderiv + 1)(nderiv + 2)(nderiv + 3)/6. The ordering is 'triangular' with the lower derivatives appearing first.

Note

This function will be called at runtime when tabulating a finite element, so its performance is critical.

Parameters

[in,out]

P

Polynomial sets, for each derivative, tabulated at points. The shape is (number of derivatives computed, number of points, basis index). The first index is the derivative. The first entry is the basis itself. Derivatives are stored in triangular (2D) or tetrahedral (3D) ordering, eg if (p, q) denotes p order derivative with respect to x and q order derivative with respect to y, [0] -> (0, 0), [1] -> (1, 0), [2] -> (0, 1), [3] -> (2, 0), [4] -> (1, 1), [5] -> (0, 2), [6] -> (3, 0),... The function basix::indexing::idx maps tuples (p, q, r) to the array index. The second index is the point, with index i corresponding to the point in row i of x. The third index is the basis function index.

Todo:

Does the order for the third index need to be documented?

Parameters

[in]	celltype	Cell type
[in]	ptype	The polynomial type
[in]	d	Polynomial degree
[in]	n	Maximum derivative order. Use n=0 for the basis only.
[in]	x	Points at which to evaluate the basis. The shape is (number of points, geometric dimension).

dim()

int basix::polyset::dim	(	cell::type	cell,
		polyset::type	ptype,
		int	d
	)

Dimension of a polynomial space.

Parameters

[in]	cell	Cell type
[in]	ptype	Polynomial type
[in]	d	Polynomial degree

Returns

The number of terms in the basis spanning a space of polynomial degree d.

nderivs()

int basix::polyset::nderivs	(	cell::type	cell,
		int	d
	)

Number of derivatives that the orthonormal basis will have on the given cell.

Parameters

[in]	cell	Cell type
[in]	d	Highest derivative order

Returns

Number of derivatives

superset()

polyset::type basix::polyset::superset	(	cell::type	cell,
		polyset::type	type1,
		polyset::type	type2
	)

Get the polyset types that is a superset of two types on the given cell.

Parameters

[in]	cell	Cell type
[in]	type1	First polyset type
[in]	type2	Decond polyset type

Returns

Superset type

restriction()

polyset::type basix::polyset::restriction	(	polyset::type	ptype,
		cell::type	cell,
		cell::type	restriction_cell
	)

Get the polyset type that represents the restrictions of a type on a subentity.

Parameters

[in]	ptype	Polyset type
[in]	cell	Cell type
[in]	restriction_cell	Cell type of the subentity

Returns

Restricted polyset type