Bernstein bases algorithms

Barycentric coordinates

A $D$-dimensional simplex $T$ is defined by $N=D+1$ vertices $\{v_1, v_2, …, v_N\}=\{v_i\}_{i∈1:N}$. The barycentric coordinates $λ(\bm{x})=\{λ^j(\bm{x})\}_{1 ≤ j ≤ N}$ are uniquely defined by:

\[\bm{x} = ∑_{1 ≤ j ≤ N} λ^j(\bm{x})v_j \quad\text{and}\quad ∑_{1≤ j≤ N} λ^j(\bm{x}) = 1,\]

as long as the simplex is non-degenerate (vertices are not all in one hyperplane).

Assuming the simplex polytopal (has flat faces), this change of coordinates is affine, and is implemented using:

\[λ(\bm{x}) = M\left(\begin{array}{c} 1\\ x^1\\ ⋮\\ x^D \end{array}\right) \quad\text{with}\quad M = \left(\begin{array}{cccc} 1 & 1 & ⋯ & 1 \\ (v_1)^1 & (v_2)^1 & ⋯ & (v_N)^1 \\ ⋮ & ⋮ & ⋯ & ⋮ \\ (v_1)^D & (v_2)^D & ⋯ & (v_N)^D \\ \end{array}\right)^{-1}\]

where the inverse exists because $T$ is non-degenerate [1], cf. functions _cart_to_bary and _compute_cart_to_bary_matrix. Additionally, we have $∂_{x^i} λ^j(\bm{x}) = M^j_{i+1}$, so

\[∇ λ^j = M^j_{2:N}.\]

The matrix $M$ is all we need that depends on $T$ in order to compute Bernstein polynomials and their derivatives, it is stored in the field x_to_λ of BernsteinBasisOnSimplex.

On the reference simplex defined by the vertices get_vertex_coordinates(SEGMENT / TRI / TET⋯):

\[\begin{aligned} v_1 & = (0\ 0\ ⋯\ 0), \\ v_2 & = (1\ 0\ ⋯\ 0), \\ ⋮ & \\ v_N & = (0\ ⋯\ 0\ 1), \end{aligned}\]

the matrix $M$ reduces to

\[λ(\bm{x}) = \Big(1-∑_{1≤ i≤ D} x^i, x^1, x^2, ⋯, x^D\Big) \quad\text{and}\quad ∂_{x^i} λ^j = δ^j_{i+1} - δ^j_1 = M^j_{i+1}.\]

Bernstein polynomials definition

The univariate Bernstein polynomials forming a basis of $\mathcal{P}_K$ are defined by

\[B_{K}^{n}(x) = \binom{K}{n} x^n (1-x)^{K-n}\qquad\text{ for } 0≤ n≤ K.\]

The $D$-multivariate Bernstein polynomials of degree $K$ relative to a simplex $T$ are defined by

\[B_{D,K}^{α}(\bm{x}) = \binom{K}{α} λ(\bm{x})^α\qquad\text{for all }α ∈\mathcal{I}_K^D\]

where

• $\mathcal{I}_K^D = \{\ α∈(\mathbb{Z}_+)^{D+1} \quad|\quad |α|=K\ \}$
• $|α|=∑_{1≤ i≤ N} α_i$
• $\binom{K}{α} = \frac{K!}{α_1 !α_2 !… α_N!}$
• $λ$ are the barycentric coordinates relative to $T$ (defined above)

The subscripts $D$ and $K$ in $B_{D,K}^{α}(x)$ are parameters, not indices, and can be omitted because they are always determined by $α$ using ${D=\#(α)-1}$ and $K=|α|$. The set $\{B^α\}_{α∈\mathcal{I}_K^D}$ is a basis of $\mathcal{P}^D_K$, implemented by BernsteinBasisOnSimplex.

Bernstein indices and indexing

Working with Bernstein polynomials requires dealing with several quantities indexed by some $α ∈ \mathcal{I}_K^D$, the polynomials themselves but also the coefficients $c_α$ of a polynomial in the basis, the domain points ${\bm{x}_α = \underset{1≤i≤N}{∑} α_i v_i}$ and the intermediate coefficients used in the de Casteljau algorithm.

These indices are returned by bernstein_terms(K,D). When storing such quantities in arrays, $∙_α$ is stored at index bernstein_term_id(α), which is the index of α in bernstein_terms(sum(α),length(α)-1).

We adopt the convention that a quantity indexed by a $α ∉ ℤ_+^N$ is equal to zero (to simplify the definition of algorithms where $α=β-e_i$ appears).

The de Casteljau algorithms

A polynomial $p ∈ \mathcal{P}^D_K$ in Bernstein form $p = ∑_{α∈\mathcal{I}^D_K}\, p_α B^α$ can be evaluated at $\bm{x}$ using the de Casteljau algorithms [1, Algo. 2.9] by iteratively computing

\[\qquad p_β^{(l)} = \underset{1 ≤ i ≤ N}{∑} λ^i\, p_{β+e_i}^{(l-1)} \qquad ∀β ∈ \mathcal{I}^D_{K-l},\]

for $l=1, 2, …, K$ where $p_α^{(0)}=p_α$, $λ=λ(\bm{x})$ and the result is $p(\bm{x})=p_𝟎^{(K)}$. This algorithm is implemented (in place) by _de_Casteljau_nD!.

But Gridap implements the polynomial bases themselves instead of individual polynomials in a basis. To compute all $B^α$ at $\bm{x}$, one can use the de Casteljau algorithm going "downwards" (from the tip of the pyramid to the base). The idea is to use the relation

\[B^α = ∑_{1 ≤ i ≤ N} λ^i B^{α-e_i}\qquad ∀α ∈ ℤ_+^N,\ |α|≥1.\]

Starting from $b^{𝟎,(0)}=B^𝟎(\bm{x})=1$, compute iteratively

\[\qquad b^{β,(l)} = \underset{1 ≤ i ≤ N}{∑} λ^i\, b^{β-e_i,(l-1)} \qquad ∀β ∈ \mathcal{I}^D_{l},\]

for $l=1,2, …, K$, where again $λ=λ(\bm{x})$ and the result is $B^α(\bm{x})=b^{α,(K)}$ for all $α$ in $\mathcal{I}^D_K$. This algorithm is implemented (in place) by _downwards_de_Casteljau_nD!. The implementation is a bit tricky, because the iterations must be done in reverse order to avoid erasing coefficients needed later, and a lot of summands disappear (when $(β-e_i)_i < 0$).

The gradient and hessian of the BernsteinBasisOnSimplex are also implemented. They rely on the following

\[∂_q B^α(\bm{x}) = K\!∑_{1 ≤ i ≤ N} ∂_qλ^i\, B^{α-e_i}(\bm{x}),\qquad ∂_t ∂_q B^α(\bm{x}) = K\!∑_{1 ≤ i,j ≤ N} ∂_tλ^j\, ∂_qλ^i\, B^{α-e_i-e_j}(\bm{x}).\]

The gradient formula comes from [1, Eq. (2.28)], and the second is derived from the first using the fact that $∂_qλ$ is homogeneous. The implementation of the gradient and hessian compute the $B^β$ using _downwards_de_Casteljau_nD! up to order $K-1$ and $K-2$ respectively, and then the results are assembled by _grad_Bα_from_Bαm! and _hess_Bα_from_Bαmm! respectively. The implementation makes sure to only access each relevant $B^β$ once per $(∇/H)B^α$ computed.

Bernstein basis generalization for $\mathcal{P}Λ$ spaces

The BarycentricPmΛBasis and BarycentricPΛBasis bases respectively implement the polynomial bases for the spaces $\mathcal{P}_r^-Λ^k(T^D)$ and $\mathcal{P}_rΛ^k(T^D)$ (we write $\mathcal{P}_r^{(-)}Λ^k$ for either one of them) derived in [2] on simplices of any dimension, for any form degree $k$ and polynomial degree $r$. These spaces include and generalize several standard FE polynomial spaces, see the Periodic Table of the Finite Elements [3].

The following notes explain the implementation in detail. For the moment, only the space with form order ${k = 0,1,D-1}$ and $D$ are available, because the forms are translated into their vector calculus proxy.

Face and form coefficients indexing

Again, a $D$-dimensional simplex $T$ is defined by $N=D+1$ vertices $\{v_1, v_2, ..., v_N\}=\{v_i\}_{i∈1:N}$. We uniquely identify a $d$-dimensional face $F$ of $T$ by the set of the $d+1$ increasing indices of its vertices:

\[F = \{F_1, F_2, ..., F_{d+1}\} \qquad\text{such that } 1≤ F_1 < F_2 < ... <F_{d+1}≤ N .\]

In particular, $T\sim \{1:N\}$. We write $F⊆ T$ for any face of $T$, including $T$ itself or its vertices. $T$ has $\binom{N}{d+1}$ $d$-dimensional faces, indexed $\forall\,1≤ F_1 < F_2 < ... < F_{d+1} ≤ N$. The dimension of a face $F$ is $\#F\;$ (length(F)), and we write ${"∀\,\#J=d+1"}$ for all the increasing index sets of the $d$-dimensional faces of $T$. We will sometimes write $J(i)$ instead of $J_i$ for readability purpose when it appears as a sub- or superscript.

Using Einstein's convention of summation on repeated indices, a degree-$k$ dimension-$D$ form $ω$ can be written in the canonical Cartesian basis as $ω = ω_I\,\text{d}x^I$, where the basis is

\[\big\{ \text{d}x^I = \underset{i∈I}{⋀}\text{d}x^{i} =\text{d}x^{I_1}∧ ...∧ \text{d}x^{I_k} \quad\big|\quad I=\{I_1, ..., I_k\} \text{ for }1≤ I_1 < ... < I_k ≤ D\big\},\]

$\{ω_I\}_I∈\mathbb{R}^\binom{D}{k}$ is the vector of coefficients of $ω$, and $\{\text{d}x^i\}_{1≤ i≤ D}$ is the canonical covector basis (basis of $\text{T}_x T$) such that $\text{d}x^i(∂_{x^j})=δ^i_j$.

These sets of indices $I,J,F$ are $k$-combinations of ${1:D/N}$, stored in Vector{Int}. A generator sorted_combinations returns a vector containing all the $D$-dimensional $k$-combinations in lexicographic order, e.g. for $D=4$, $k=2$

\[\{1,2\},\ \{1,3\},\ \{1,4\},\ \{2,3\},\ \{2,4\},\ \{3,4\},\]

and combination_index computes the index of a combination in this vector.

This order define the linear indices of basis k-forms such as $I$ and $J$. The order depends on the dimension D.: $\{2,3\}$ is the third length-$2$ combination for $D=3$ but the fourth for $D=4$.

The faces indices $F$ are ordered with right-to-left lexicographic, obtain with kwarg right_to_left=true of sorted_combinations and combination_index. The combination indices are compared from last to first, so above $\{2,3\}$ would swaps with $\{1,4\}$. This order is consistent with the vertices ordering of simplices faces, e.g. get_faces(TET,k,0)

Translation between forms and vectors

By default, the polynomial forms of order $k = 0,1,D-1$ and $D$ are translated into their equivalents in the standard vector calculus framework (assuming the simplex is Euclidean).

This change of coordinate is implemented by _basis_forms_components, the indices of a basis b::BarycentricP(m)ΛBasis are stored in b._indices.components. For ${D=2}$ and ${k=1}$, the default proxy is $ω♯$. The user may choose the $(⋆ω)♯$ proxy (for div-conforming spaces) using the kwarg rotate_90=true.

Geometric decomposition

The main feature of the BarycentricP(m)ΛBasis bases is that each basis polynomial $ω^{α,J}$ is associated with a face $F$ of $T$ via ${F=⟦α⟧∪J}$ with $J$ a face of $F$ and $α$ a Bernstein index whose associated domain point $\boldsymbol{x}_α$ is geometrically inside $F$. Importantly, the trace of $ω^{α,J}$ on another face $G⊆ T$ is zero when $G$ does not contain $F$:

\[F\not⊆ G\ \rightarrow\ \text{tr}_G\, ω^{α,J} = 0, \quad\forall F,G ⊆ T,\ \forall α,J \text{ s.t. }\llbracket α\rrbracket\cup J = F,\]

including any face $G\neq F$ of dimension less or equal that of $F$.

These basis polynomials $ω^{α,J}$ are called bubble functions associated to $F$, the space they span is called $\mathring{\mathcal{P}}_r^{(-)}Λ^k(T,F)$. There are no bubble functions of degree $k$ on faces of dimension $<k$, so the spaces $\mathcal{P}_r^{(-)}Λ^k(T)$ admit the geometric decomposition:

\[\mathcal{P}_r^{(-)}Λ^k(T) = \underset{F⊆ T}{\oplus}\ \mathring{\mathcal{P}}_r^{(-)}Λ^k(F) = \underset{k≤d≤D}{\oplus}\underset{\quad F=1≤ F_1 < ... < F_{d+1} ≤ N}{\oplus}\ \mathring{\mathcal{P}}_r^{(-)}Λ^k(T,F).\]

Bubble functions $\mathring{\mathcal{P}}_r^-Λ^k$

The $\mathcal{P}^-Λ$ type bubble basis polynomials associated to a face $F⊆T$ defined by [2, Th. 6.1-4] are

\[\mathring{\mathcal{P}}_r^-Λ^k(T,F) = \text{span}\big\{ ω̄^{α,J} = B^α φ^J \ \big| \ α∈\mathcal{I}_{r-1}^D,\ \#J=k\!+\!1,\ ⟦α⟧∪J=F,\ α_i=0 \text{ if } i< \text{min}(J) \big\}\]

where $B^α$ are the scalar Bernstein polynomials implemented by BernsteinBasisOnSimplex, and $φ^J$ [2, Eq. (6.3)] are the Whitney forms:

\[φ^J = \sum_{1≤l≤k+1} (-1)^{l+1} λ^{J(l)} \, \text{d}λ^{J\backslash l} \quad\text{where}\quad \text{d}λ^{J\backslash l} = \underset{j∈J\backslash \{J_l\} }{⋀}\text{d}λ^{j},\]

$φ^J$ is a $k$-form of polynomial order $1$.

Given $k,r$ and $D$, the function PmΛ_bubbles(r,k,D) computes, for each $d$-face$F$ "owning" a bubble space, the indices necessary to compute its bubble polynomials. PmΛ_bubbles is used as follows:

for (F, bubble_functions) in PΛ_bubbles(r,k,D)  # d = length(F)
    for (w, α, α_id, J, sub_J_ids, sup_α_ids) in bubble_functions
        # do stuff for ω̄^{α,J}
    end
end

where

• w is the index of $ω̄^{α,J}$ in the whole BarycentricPmΛBasis,
• α is a Vector{Int},
• α_id is bernstein_term_id(α), the index of Bα in the scalar BernsteinBasisOnSimplex,
• J is a Vector{Int},
• sub_J_ids is a ::Vector{Int} are the combination_index of each $J\backslash \{J(l)\}$ for $1\leq l\leq \#J$, taken among the combinations of $1\!:\!N$,
• sup_α_ids is a ::Vector{Int} are the bernstein_term_id of each $α+e_i$ for $1\leq i\leq \#α$.

The implementation is flexible enough to select a subset of the bubble spaces, the bubbles of a b::BarycentricPmΛBasis are obtained via get_bubbles(b) (do NOT modify them).

We now need to express $\text{d}λ^{J\backslash l}$ in the Cartesian basis ${\text{d}x^I}$. In a polytopal simplex $T$ (flat faces), the 1-forms $\text{d}λ^j:=\text{d}(λ^j)$ is homogeneous, its coefficients in the canonical basis are derived by

\[\text{d}λ^j = (\nabla(λ^j))^♭ = δ_{ki}∂_{k}λ^j\,\text{d}x^i = ∂_{i}λ^j\,\text{d}x^i = M^j_{i+1}\,\text{d}x^i\]

where ${}^♭$ is the flat map, the metric $g_{ki}=δ_{ki}$ is trivial and $M^j_{i+1}$ are components of the barycentric change of coordinate matrix $M$ introduced in the Barycentric coordinates section above.

So the exterior products $\text{d}λ^{J\backslash l}$ are expressed using the $k$-minors $m_I^{J\backslash l}$ of $M^\intercal$ as follows:

\[\text{d}λ^{J\backslash l} = m_I^{J\backslash l}\text{d}x^I \quad\text{where}\quad m_I^J = \text{det}\big( (∂_{I(i)}λ^{J(j)})_{1≤ i,j≤ k} \big) = \text{det}\big( (M^{J(j)}_{I(i)+1})_{1≤ i,j≤ k} \big),\]

and we obtain the components of $ω̄^{α,J}=B^α φ^J$ in the basis $\mathrm{d}x^I$

\[ω̄_{I}^{α,J} = B^α \sum_{1≤l≤k+1} (-1)^{l+1} λ^{J(l)} \, m_I^{J\backslash l}.\]

For $k \leq 1$, BarycentricPmΛBasis also offers the flavor=:BMM variant, which substitutes the bare monomial $λ^α = \binom{|α|}{α}^{-1}B^α$ for $B^α$.

The $\binom{D}{k}\binom{N}{k}$ coefficients $\{m_I^{J}\}_{I,J}$ are constant in $T$ and are pre-computed from $M$ in _compute_PmΛ_basis_coefficients! at the creation of BarycentricPmΛBasis and stored in its field m.

Finally, the pseudocode to evaluate our basis $ω̄$ of $\mathcal{P}_r^-Λ^k(T)$ at $\boldsymbol{x}$ is (spelled with the implementation's names, so its indices are Julia identifiers rather than the notation above)

compute λ(x)
compute B(x) = { Bα(λ(x)) } for all |α|=r-1

for (F, bubble_functions) in get_bubbles(b)
    for (w, α, α_id, J, sub_J_ids) in bubble_functions

        ω̄_w = 0 # ω̄^{α,J}
        for (l, J_sub_Jl_id) in enumerate(sub_J_ids)
            λ_j = λ[J[l]]
            m_J_l = m[J_sub_Jl_id] # is a coordinate vector for all I
            ω̄_w += -(-1)^l * λ_j * m_J_l
        end

        Bα = B[α_id]
        ω̄[w] = Bα * ω̄_w
    end
end

Bubble functions $\mathring{\mathcal{P}}_rΛ^k$

The $\mathcal{P}Λ$ type bubble basis polynomials associated to a face $F⊆T$ defined by [2, Th. 6.1-2] – where the basis function Eq. (8.3) replace Eq. (8.1) – are

\[\mathring{\mathcal{P}}_rΛ^k(T,F) = \text{span}\big\{ ω^{α,J}=B^α Ψ^{α,J} \quad\big|\quad \ α∈\mathcal{I}_{r}^D,\ \#J=k,\ ⟦α⟧∪J=F,\ α_i=0 \text{ if } i< \text{min}(F \backslash J) \big\},\]

where $Ψ^{α,J}$ [2, Eq. (8.3)] are defined by

\[Ψ^{α,J} = \underset{j∈J}{⋀} Ψ^{α,F(α,J),j} \quad\text{and}\quad Ψ^{α,F,j} = \mathrm{d}λ^j - \frac{α_j}{|α|}\sum_{l∈F}\mathrm{d}λ^l,\]

where $F(α,J)=⟦α⟧∪J$. get_bubbles(b::BarycentricPΛBasis) provides the bubbles of b, their bubble function indices are (w, α, α_id, J) only (do NOT modify them).

Again, we need their components in the Cartesian basis $\mathrm{d}x^I$:

\[Ψ^{α,F,j} = M^j_{i+1}\mathrm{d}x^i - \frac{α_j}{|α|}\sum_{l∈F}M^l_{i+1}\mathrm{d}x^i = \big(M^j_{i+1} - \frac{α_j}{|α|}\sum_{l∈F}M^l_{i+1}\big)\mathrm{d}x^i\]

so

\[Ψ^{α,F,j} = ψ_{i}^{α,F,j} \mathrm{d}x^i \quad\text{where}\quad ψ_{i}^{α,F,j} = M^j_{i+1} - \frac{α_j}{|α|}\sum_{l∈F}M^l_{i+1}\]

and

\[Ψ^{α,J} = ψ_I^{α,J} \mathrm{d}x^I \quad\text{where}\quad ψ_I^{α,J} = \text{det}\big( (ψ_{i}^{α,F,j})_{i∈I,\,j∈J} \big).\]

For $k \leq 1$, BarycentricPΛBasis also offers the flavor=:BMM direction forms, which substitute the support indicator $s(α)_j = 1_{α_j>0}$ for $α_j$ and $|\mathrm{supp}(α)|$ for $|α|$ in the formula above. All the $α$ of a given support then share their direction forms.

Finally, the $\binom{D}{k}$ components of $ω^{α,J}=B^α Ψ^{α,J}$ in the basis $\mathrm{d}x^I$ are

\[ω_{I}^{α,J} = B^α\, ψ_I^{α,J},\]

where the $\binom{D+r}{k+r}\binom{r+k}{k}\binom{D}{k} =\mathrm{dim}(\mathcal{P}_rΛ^k(T^D))\times\# (\{\mathrm{d}x^I\}_I)$ coefficients $ψ_I^{α,J}$ depend only on $T$ and are pre-computed in _compute_PΛ_basis_form_coefficient! at the construction of BarycentricPΛBasis and stored in its field Ψ.

The pseudocode to evaluate our basis $ω$ of $\mathcal{P}_rΛ^k(T)$ at $\boldsymbol{x}$ is

compute λ(x)
compute B(x) = { Bα(λ(x)) } for all |α|=r

for (F, bubble_functions) in get_bubbles(b)
    for (w, α, α_id) in bubble_functions
        Bα = B[α_id]
        ω[w] = Bα * Ψ[w]
    end
end

Gradient and Hessian of the coefficient vectors

Let us derive the formula for the gradient and hessian of the basis forms coefficient vectors $\{ω̄_{I}^{α,J}\}_I$ and $\{ω_{I}^{α,J}\}_I$. We will express them in function of the scalar Bernstein polynomial derivatives already implemented by BernsteinBasisOnSimplex. They are only supported for scalar or VectorValue'd bases (vector calculus style).

Coefficient vector $\{ω_{I}^{α,J}\}_I$

Recall $ω_{I}^{α,J} = B^α\, ψ_I^{α,J}$. The derivatives are easy to compute because only $B^α$ depends on $\boldsymbol{x}$, leading to

\[∂_q\, ω_{I}^{α,J} = ψ_I^{α,J}\; ∂_q B^α,\qquad\text{or}\qquad ∇ω^{α,J} = ∇B^α ⊗ ψ^{α,J}\\ ∂_t∂_q\, ω_{I}^{α,J} = ψ_I^{α,J}\; ∂_t∂_q B^α\qquad\text{or}\qquad ∇∇ω^{α,J} = ∇∇B^α ⊗ ψ^{α,J}.\]

where $∇$ and $∇∇$ are the standard gradient and hessian operators and $ψ^{α,J}$ is seen as a length-$\binom{D}{k}$ vector with no variance (not as a $k$-form, which is an order $k$ covariant tensor).

Coefficient vector $\{ω̄_{I}^{α,J}\}_I$

Recall $ω̄_{I}^{α,J} = B^α \sum_{1≤l≤k+1} (-1)^{l+1} λ^{J(l)} \, m_I^{J\backslash l}$, the derivatives are not immediate to compute because both $B^α$ and $λ^{J(l)}$ depend on $\boldsymbol{x}$, let us first use $B^α λ^{J(l)} = \frac{α_{J(l)} + 1}{|α|+1}B^{α+e(J,l)}$ where $e(J,l) = \big(δ_i^{J_l}\big)_{1≤ i≤ N}$ to write the coefficients in Bernstein form as follows

\[ω̄_{I}^{α,J} = B^α \sum_{1≤l≤k+1} (-1)^{l+1} λ^{J(l)} \, m_I^{J\backslash l} = \frac{1}{r}\sum_{1≤l≤k+1} (-1)^{l+1} (α_{J(l)} +1)\ B^{α+e(J,l)}\, m_I^{J\backslash l},\]

where $|α|+1$ was replaced with $r$, the polynomial degree of $ω̄_{I}^{α,J}$. As a consequence, for any Cartesian coordinate indices $1≤ p,q≤ D$, we get

\[∂_q ω̄_{I}^{α,J} = \frac{1}{r}\sum_{1≤l≤k+1} (-1)^{l+1} (α_{J(l)} +1) \ ∂_q B^{α+e(J,l)}\, m_I^{J\backslash l},\\ ∂_t∂_q ω̄_{I}^{α,J} = \frac{1}{r}\sum_{1≤l≤k+1} (-1)^{l+1} (α_{J(l)} +1) \ ∂_t∂_q B^{α+e(J,l)}\, m_I^{J\backslash l}.\]

In tensor form, this is

\[\mathrm{D}ω^{α,J} = \frac{1}{r}\sum_{1≤l≤k+1} (-1)^{l+1} (α_{J(l)} +1)\ \mathrm{D}\!B^{α+e(J,l)} ⊗ m^{J\backslash l}\\\]

where $\mathrm{D}$ is $∇$ or $∇∇$, the standard gradient and hessian operators, and $ω̄^{α,J}$ is again seen as a length-$\binom{D}{k}$ vector with no variance.

Exterior derivative of the basis forms (to be implemented)

The exterior derivative of a $k$-form $ω=ω_{\tilde{I}}\,\mathrm{d}x^{\tilde{I}}$ is the $k\!+\!1$-form

\[\mathrm{d}ω = ∂_i ω_{\tilde{I}}\, \mathrm{d}x^i ∧ \mathrm{d}x^{\tilde{I}} \qquad\text{ where }\qquad \#\tilde{I} = k\]

We need to express $\mathrm{d}ω$ in the basis of $k+1$ forms. Let $I$ such that $\#I=k\!+\!1$ with $k<D$ (otherwise $\mathrm{d}ω=0$). Because the exterior product is alternating, the coefficients that contribute to $(\mathrm{d}ω)_{I}$ are $∂_i ω_{\tilde{I}}$ for which $i=I_q$ and $\tilde{I} = I\backslash \{I_q\}$ with $1≤ q≤ k+1$, so one can deduce

\[(\mathrm{d}ω)_I = \underset{1≤ q≤ k+1}{\sum} (-1)^{q-1}\ ∂_{I(q)} ω_{I\backslash q}.\]

Polynomial forms $\mathrm{d}\,ω̄^{α,J}$

For all $|α|=r\!-\!1$, $\,\#J=k\!+\!1$ and $\#I = k\!+\!1$ (with $k\!<\!D$):

\[(\mathrm{d}\,ω̄^{α,J})_I = \frac{1}{r}\underset{1≤ l≤ k+1}{\sum} (-1)^{l+1}(α_{J(l)}+1) \underset{1≤ q≤ k+1}{\sum} (-1)^{q-1}\ m_{I\backslash q}^{J\backslash l}\ ∂_{I(q)} B^{α+e(J,l)}\;\]

Polynomial forms $\mathrm{d}\,ω^{α,J}$

For all $|α|=r$, $\,\#J=k$ and $\#I = k\!+\!1$ (with $k\!<\!D$):

\[(\mathrm{d}\,ω^{α,J})_I = \underset{1≤ q≤ k+1}{\sum} (-1)^{q-1}\ ψ_{I\backslash q}^{α,J}\ ∂_{I(q)} B^α.\]

Hodge operator of the basis forms

The Hodge operator of the canonical basis forms $\mathrm{d}x^I$ in an Euclidean (Riemannian) space is

\[\star \mathrm{d}x^I = \mathrm{sgn}\left(I\!*\!\bar{I}\,\right)\mathrm{d}x^{\bar{I}}\]

where $\bar{I}$ is the complement of $I$ in $1\!:\!D$, that is the only combination such that the concatenated permutation $I\!*\!\bar{I}$ is a permutation of $1\!:\!D$. $\bar{I}$ is implemented by _complement(I,D). _combination_sign(I) computes the sign of $I\!*\!\bar{I}$.

Analytical formulas in the reference simplex

In the reference simplex $\hat{T}$, the vertices and thus coefficients of $M$ are known at compile time, so the coefficients $m_I^J$ and $ψ_I^{α,J}$ in $\hat{T}$, denoted by $\hat{m}_I^J$ and $\hat{ψ}_I^{α,J}$ respectively, could be hard-coded at compile time in @generated functions to avoid storing them in the basis and accessing them at runtime.

In the following, analytical formulas are derived for them. They are implemented and used in tests to validate the numerical computation of the coefficients.

Coefficients $\hat{m}_I^J$

It was shown in the Barycentric coordinates section above that $M^j_{i+1} = δ^j_{i+1} - δ^j_1$. Let $\#I=k$ and $\#J=k$. We need to compute the determinant of the matrix

\[\hat{M}^J_I=(δ^{J(j)}_{I(i)+1}-δ^{J(j)}_{1})_{1≤ i,j≤ k}.\]

Let us define:

• $s=δ_1^{J_1}$, that indicates if $\hat{M}^J_I$ contains a column of $-1$,
• $p = \text{min } \{j\,|\, I_j+1 ∉J\}$ where $\text{min}\,∅=0$, the index of the first row of $\hat{M}^J_I$ containing no $1$. $p=0$ if and only if $\hat{M}^J_I$ is the identity matrix,
• $n =\# \{\ i\ |\ i>s,\, J_i-1∉I\}$, the number of columns of zeros of $\hat{M}^J_I$.

Then it can be shown that $k-n$ is the rank of $\hat{M}^J_I$, and that

\[\hat{m}_I^J = \mathrm{det}(\hat{M}^J_I) = (-1)^{p(I,J)}δ_0^{n(I,J)},\]

where the dependency of $p,n$ on $I,J$ is made explicit, so in $\hat{T}$, there is

\[ω̄_{I}^{α,J} = B^α \sum_{1≤l≤k+1} (-1)^{l+1} λ^{J(l)} \, \hat{m}_I^{J\backslash l}.\]

Coefficients $\hat{ψ}_I^{α,J}$

The expression of $ψ_{i}^{α,F,j}$ in $\hat{T}$ is

\[\hat{ψ}_{i}^{α,F,j} = M^j_{i+1} - \frac{α_j}{|α|}\sum_{l∈F}\big(δ^l_{i+1} - δ^l_1\big) = M^j_{i+1} + \sum_{l∈F}\big( δ^l_1-δ^l_{i+1} \big)\frac{α_j}{|α|}\]

leading to

\[\hat{ψ}_I^{α,J} = \mathrm{det}\Big(\hat{M}^J_I + u\,v^{\intercal}\Big) \quad\text{where}\quad u^i = \sum_{l∈F}\big(δ^l_1-δ^l_{I(i)+1}\big), \qquad v^j = \frac{α_{J(j)}}{|α|}.\]

We can use the following matrix determinant lemma:

\[\mathrm{det}(\hat{M}^J_I + uv^\intercal) = \mathrm{det}(\hat{M}^J_I) + v^\intercal\mathrm{adj}(\hat{M}^J_I)u.\]

The determinant $\mathrm{det}(\hat{M}^J_I)=\hat{m}_I^{J}$ was computed above, but $\mathrm{adj}(\hat{M}^J_I)$, the transpose of the cofactor matrix of $\hat{M}^J_I$, is also needed. Let $s=δ_1^{J_1}$, $n$ and $p$ be defined as above, and additionally define

• $q = \text{min } \{j\,|\,j>p,\ I_j+1 ∉J\}$, the index of the second row of $\hat{M}^J_I$ containing no $1$ ($q=0$ if there isn't any),
• $m = \text{min } \{i\,|\,i>s,\ J_i-1 ∉I\}$, the index of the first column of $\hat{M}^J_I$ containing only zeros ($m=0$ if there isn't any).

Then the following table gives the required information to apply the matrix determinant lemma and formulas for $\hat{ψ}_I^{α,J}$

\[\begin{array}{|c|c|c|c|c|} \hline s & n & \mathrm{rank}\hat{M}^J_I & \mathrm{adj}\hat{M}^J_I & \hat{ψ}_I^{α,J} \\ \hline \hline 0 & 0 & k & δ^i_j & 1 + u \cdot v\\ \hline 0 & 1 & k-1 & (-1)^{m+p}δ^i_m δ^p_j & (-1)^{m+p}v^m u^p \\ \hline 1 & 0 & k & (-1)^p(δ^{J(i)}_{I(j)+1}-δ^{J(j)}_{p})& (-1)^p(1-u^p|v|+\underset{1≤ l<p}{\sum}v^{l+1}u^l + \underset{p<l≤ k}{\sum}v^{l}u^l) \\ \hline\hspace{1mm} 1 & 1 & k-1 & (-1)^{m+p+q}δ^i_m(δ^q_j-δ^p_j) & (-1)^{m+p+q}v^m(u^q-u^p) \\ \hline 0/1 & \geq 2 & ≤ k-2 & 0 & 0 \\ \hline \end{array}\]

In this table, $m$, $p$ and $q$ depend on $I$ and $J$, $u$ depends on $F$ and $I$, and $v$ depends on $α$ and $J$.

Low level docstrings

_compute_cart_to_bary_matrix(vertices, ::Val{N})

_cart_to_bary(x::Point{D,T}, x_to_λ)

_de_Casteljau_nD!(c, λ,::Val{K},::Val{D},::Val{Kf}=Val(0))

_downwards_de_Casteljau_nD!(c, λ,::Val{K},::Val{D},::Val{K0}=Val(1))

_basis_forms_components(D,k,DG_style,rotate_90)

v[I_proxy_id] = I_proxy_sgn * ω[I_id]

_update_φ_αF!(φ_αF, b, α, F, r, flavor)

φ^{α,F,j} = dλʲ - (sⱼ/|s|) Σ_{l∈F} dλˡ,     1 ≤ j ≤ N.

_combination_sign(I)

_complement(I, D)

References

[1] M.J. Lai & L.L. Schumaker, Spline Functions on Triangulations, Chapter 2 - Bernstein–Bézier Methods for Bivariate Polynomials, pp. 18 - 61.

[2] D.N. Arnold, R.S. Falk & R. Winther, Geometric decompositions and local bases for spaces of finite element differential forms, Computer Methods in Applied Mechanics and Engineering

[3] D.N. Arnold and A. Logg, Periodic Table of the Finite Elements, SIAM News, vol. 47 no. 9, November 2014.