Geometric finite elements (FEEC)

Basics

In conforming finite element (FE) methods, the general principle of approximating a function \(u \in V\) in a space \(V\) is

\[u \approx u_h = \sum_{i=1}^{M} u_i \Lambda_i \quad \in V_h\,,\]

where \(\Lambda_i \in V\) are linearly independent basis functions, \((u_i)_i \in \mathbb R^M\) are called coefficients and \(M \in \mathbb N\) is the dimension of the subspace \(V_h \subset V\) spanned by the basis functions. The thing the differentiates FE methods from spectral methods is the fact that the \(\Lambda_i\) have local support around some grid point \(x_i\). A nice feature of FE methods is that \(u_h \in V\), which leads often to an elegant analysis of the method. The implementation of a FE method consists in writing down a system of equations for the coefficients \((u_i)_i \in \mathbb R^M\), leaving the basis functions fixed.

We denote by “geometric finite elements” the FE method based on finite element exterior calculus (FEEC). In 3D, it consists of a sequence of discrete (=finite dimensional) spaces \(V_h^n,\,0\leq n \leq 3\) that satisfy the 3d Derham diagram:

Fig. 1: 3d Derham diagram.

The upper row is a sequence of continuous function spaces, where the last space \(L^2(\Omega)\) is the space of square-integrable functions on \(\Omega\). The other spaces are:

\(H^1(\Omega)\): functions \(f \in L^2(\Omega)\) for which \(\nabla f \in (L^2(\Omega))^3\)

\(H(\textnormal{curl}, \Omega)\): functions \(\mathbf u \in (L^2(\Omega))^3\) for which \(\nabla \times \mathbf u \in (L^2(\Omega))^3\)

\(H(\textnormal{div}, \Omega)\): functions \(\mathbf u \in (L^2(\Omega))^3\) for which \(\nabla \cdot \mathbf u \in L^2(\Omega)\)

More information on these spaces can be found in many textbooks on FE methods, for instance in Brezzi, Fortin, “Mixed and Hybrid Finite Element Methods”.

Struphy uses a conforming” FE method,

\[V_h^0 \subset H^1(\Omega),\qquad V_h^1 \subset H(\textnormal{curl}, \Omega),\qquad V_h^2 \subset H(\textnormal{div}, \Omega),\qquad V_h^3 \subset L^2(\Omega)\]

The operators \(\Pi_n,\,0\leq n \leq 3\) project (\(\Pi_n^2 = \Pi_n\)) into the subspaces :

\[\Pi_0: \, H^1(\Omega) \to V_h^0,\qquad \Pi_1:\, H(\textnormal{curl}, \Omega) \to V_h^1,\qquad \Pi_2:\, H(\textnormal{div}, \Omega) \to V_h^2,\qquad \Pi_3:\, L^2(\Omega) \to V_h^3\]

These spaces and the associated operators of the Derham diagram have been implemented in the open-source package Psydac. Struphy interfaces to this library by means of the class struphy.feec.psydac_derham.Derham.

Both the continuous and the discrete spaces form a complex, which means

\[\textnormal{curl}\,\textnormal{grad} = 0\,,\qquad \textnormal{div}\,\textnormal{curl} = 0\]

holds on both levels. Moreover, the above diagram can be viewed as composed of three commuting diagrams, namely

(1)\[\Pi_1\textnormal{grad} = \textnormal{grad}\,\Pi_0\,,\qquad \Pi_2\textnormal{curl} = \textnormal{curl}\,\Pi_1\,,\qquad \Pi_3\textnormal{div} = \textnormal{div}\,\Pi_2\]

In Struphy the discrete spaces \(V_h^n,\,0\leq n \leq 3\) are spanned by tensor-product B-spline basis functions. Building blocks are the Uni-variate spline spaces of degree \(p\) with basis functions \(N^p(\eta)\) and \(D^{p -1}(\eta)\) defined on the unit interval \(\eta \in [0, 1]\). Hence, in Struphy the simualtion domain is always the unit cube, \((\eta_1,\eta_2,\eta_3) \in \Omega = [0, 1]^3\). The tensor-product basis functions are denoted as follows:

\begin{align} \Lambda^0_{ijk}(\eta_1,\eta_2,\eta_3) &= N_i^{p_1}(\eta_1) N_j^{p_2}(\eta_2) N_k^{p_3}(\eta_3)\,,\qquad &&\mathbf \Lambda^0 = (\Lambda^0_{ijk})_{i=1, j=1, k=1}^{n_1, n_2, n_3} \qquad \\\\ \vec{\Lambda}^1_{1,ijk}(\eta_1,\eta_2,\eta_3) &= \begin{pmatrix} D_i^{p_1 - 1}(\eta_1) N_j^{p_2}(\eta_2) N_k^{p_3}(\eta_3) \\ 0 \\ 0 \end{pmatrix}\,,\qquad &&\vec{\mathbf \Lambda}^1_1 = (\vec \Lambda^1_{1,ijk})_{i=1, j=1, k=1}^{d_1, n_2, n_3} \\\\ \vec{\Lambda}^1_{2,ijk}(\eta_1,\eta_2,\eta_3) &= \begin{pmatrix} 0 \\ N_i^{p_1}(\eta_1) D_j^{p_2 - 1}(\eta_2) N_k^{p_3}(\eta_3) \\ 0 \end{pmatrix}\,,\qquad &&\vec{\mathbf \Lambda}^1_2 = (\vec \Lambda^1_{2,ijk})_{i=1, j=1, k=1}^{n_1, d_2, n_3} \\\\ \vec{\Lambda}^1_{3,ijk}(\eta_1,\eta_2,\eta_3) &= \begin{pmatrix} 0 \\ 0 \\ N_i^{p_1}(\eta_1) N_j^{p_2}(\eta_2) D_k^{p_3 - 1}(\eta_3) \end{pmatrix}\,,\qquad &&\vec{\mathbf \Lambda}^1_3 = (\vec \Lambda^1_{3,ijk})_{i=1, j=1, k=1}^{n_1, n_2, d_3} \\\\ \vec{\Lambda}^2_{1,ijk}(\eta_1,\eta_2,\eta_3) &= \begin{pmatrix} N_i^{p_1}(\eta_1) D_j^{p_2 - 1}(\eta_2) D_k^{p_3 - 1}(\eta_3) \\ 0 \\ 0 \end{pmatrix}\,,\qquad &&\vec{\mathbf \Lambda}^2_1 = (\vec \Lambda^2_{1,ijk})_{i=1, j=1, k=1}^{n_1, d_2, d_3} \\\\ \vec{\Lambda}^2_{2,ijk}(\eta_1,\eta_2,\eta_3) &= \begin{pmatrix} 0 \\ D_i^{p_1 - 1}(\eta_1) N_j^{p_2}(\eta_2) D_k^{p_3 - 1}(\eta_3) \\ 0 \end{pmatrix}\,,\qquad &&\vec{\mathbf \Lambda}^2_2 = (\vec \Lambda^2_{2,ijk})_{i=1, j=1, k=1}^{d_1, n_2, d_3} \\\\ \vec{\Lambda}^2_{3,ijk}(\eta_1,\eta_2,\eta_3) &= \begin{pmatrix} 0 \\ 0 \\ D_i^{p_1 - 1}(\eta_1) D_j^{p_2 - 1}(\eta_2) N_k^{p_3}(\eta_3) \end{pmatrix}\,,\qquad &&\vec{\mathbf \Lambda}^2_3 = (\vec \Lambda^2_{3,ijk})_{i=1, j=1, k=1}^{d_1, d_2, n_3} \\\\ \Lambda^3_{ijk}(\eta_1,\eta_2,\eta_3) &= D_i^{p_1 - 1}(\eta_1) D_j^{p_2 - 1}(\eta_2) D_k^{p_3 - 1}(\eta_3)\,,\qquad &&\mathbf \Lambda^3 = (\Lambda^3_{ijk})_{i=1, j=1, k=1}^{d_1, d_2, d_3} \qquad \end{align}

Elements of the discrete spaces are linear combinations of the respective basis functions:

\begin{align} V_h^0 &= \textnormal{span}(\mathbf \Lambda^0)\,,\qquad p_h^0(\eta_1,\eta_2,\eta_3) = \sum_{ijk} p_{ijk}\, \Lambda^0_{ijk}(\eta_1,\eta_2,\eta_3) = \mathbf p^\top \mathbf \Lambda^0 \quad \in \,V_h^0 \\ V_h^1 &= \textnormal{span}(\vec{\mathbf \Lambda}^1_1, \vec{\mathbf \Lambda}^1_2, \vec{\mathbf \Lambda}^1_3)\,, \\ &\qquad\qquad \mathbf E_h^1(\eta_1,\eta_2,\eta_3) = \sum_{\mu=1}^3 \sum_{ijk} e_{\mu, ijk}\,\vec{\Lambda}^1_{\mu, ijk}(\eta_1,\eta_2,\eta_3) = \sum_{\mu=1}^3 \mathbf e_\mu^\top \vec{\mathbf \Lambda}^1_\mu = \mathbf e^\top \vec{\mathbf \Lambda}^1 \quad \in \,V_h^1 \\ V_h^2 &= \textnormal{span}(\vec{\mathbf \Lambda}^2_1, \vec{\mathbf \Lambda}^2_2, \vec{\mathbf \Lambda}^2_3)\,, \\ &\qquad \qquad \mathbf B_h^2(\eta_1,\eta_2,\eta_3) = \sum_{\mu=1}^3 \sum_{ijk} b_{\mu, ijk}\,\vec{\Lambda}^2_{\mu, ijk}(\eta_1,\eta_2,\eta_3) = \sum_{\mu=1}^3 \mathbf b_\mu^\top \vec{\mathbf \Lambda}^2_\mu = \mathbf b^\top \vec{\mathbf \Lambda}^2 \quad \in \,V_h^2 \\ V_h^3 &= \textnormal{span}(\mathbf \Lambda^3)\,,\qquad n_h^3(\eta_1,\eta_2,\eta_3) = \sum_{ijk}^{\phantom{3}} n_{ijk}\, \Lambda^3_{ijk}(\eta_1,\eta_2,\eta_3) = \mathbf n^\top \mathbf \Lambda^3 \quad \in \,V_h^3 \end{align}

The discrete spline functions are represented by their FE coefficients \(\mathbf p \in \mathbb R^{N_0}\), \(\mathbf e \in \mathbb R^{N_1}\), \(\mathbf b \in \mathbb R^{N_2}\) and \(\mathbf n \in \mathbb R^{N_3}\). The class for callable spline functions in Struphy is struphy.feec.derham.SplineFuntion. In particular, each SplineFuntion object has

1. the property vector holding the FE coefficients, along with a setter method,

2. a __call__() method for evaluating the spline at (an array of) points \((\eta_{ijk})_{ijk}\).

A splien function can be created via the factory function struphy.feec.derham.Derham.create_spline_function(). The space dimensions are products of the uni-variate dimensions:

\begin{align} \textnormal{dim} V_h^0 &= N_0 = n_1 n_2 n_3 \\\\ \textnormal{dim} V_h^1 &= N_1 = d_1 n_2 n_3 + n_1 d_2 n_3 + n_1 n_2 d_3 \\\\ \textnormal{dim} V_h^2 &= N_2 = n_1 d_2 d_3 + d_1 n_2 d_3 + d_1 d_2 n_3 \\\\ \textnormal{dim} V_h^3 &= N_3 = d_1 d_2 d_3 \end{align}

The derivatives act as follows on the FE coefficients:

\begin{align} \textnormal{grad}\, p_h^0 &= (\mathbb G \mathbf p)^\top \vec{\mathbf \Lambda}^1 \quad \in \,V_h^1 \\\\ \textnormal{curl}\, \mathbf E_h^1 &= (\mathbb C \mathbf e)^\top \vec{\mathbf \Lambda}^2 \quad \in \,V_h^2 \\\\ \textnormal{div}\, \mathbf B_h^2 &= (\mathbb D \mathbf b)^\top \mathbf \Lambda^3 \quad \in \,V_h^3 \end{align}

Here, we introduced the discrete linear operators \(\mathbb G: \mathbb R^{N_0} \to \mathbb R^{N_1}\), \(\mathbb C: \mathbb R^{N_1} \to \mathbb R^{N_2}\) and \(\mathbb D: \mathbb R^{N_2} \to \mathbb R^{N_3}\), which satisfy the complex property

\[\mathbb C\, \mathbb G = 0\,,\qquad \mathbb D\, \mathbb C = 0\]

Note

A struphy userguide for the operators \(\mathbb G\), \(\mathbb C\) and \(\mathbb D\) and for the projection operators \(\Pi_n\) is given in the Tutorials.

The projectors \(\Pi_n,\, 0\leq n \leq 3\) into the discrete spaces \(V_h^n,\, 0\leq n \leq 3\) are constructed such that the commuting relations (1) hold. They are defined by Degrees of freedom (DOFs) \(\sigma^n,\, 0\leq n \leq 3\) which are functionals that map continuous functions to real numbers. The number of DOFs for each space is equal to its dimension. A common DOF is the evaluation of a function at a certain point in \(\Omega\); if we choose \(N_0\) such distinct points and call the corresponding DOFs \((\sigma^0_{ijk})_{i=1,j=1, k=1}^{n_1,n_2,n_3}\), we can define a projector \(\Pi_0:H^1 \to V_h^0\) by

(2)\[\sigma^0_{ijk}(p_h^0) = \sigma^0_{ijk}(p^0)\,,\qquad p^0 \in H^1 \mapsto p_h^0 \in V_h^0\]

These are \(N_0\) equations for the coefficients \(\mathbf p \in \mathbb R^{N_0}\) of \(p_h^0 \in V_h^0\). The above formula is nothing else than classical interpolation and is used for the \(\Pi_0\)-projector in Struphy. The other projectors are a mix of interpolation and “histopolation” (integration between two interpolation points), as depicted in Fig. 2. For instance, to project into the space \(V_h^1\) (1-form), the first component is “histopolated” in \(\eta_1\)-direction (orange lines) and interpolated in \((\eta_2, \eta_3)\) (blue points); the second component is “histopolated” in \(\eta_2\)-direction and interpolated in \((\eta_1, \eta_3)\); the third component is “histopolated” in \(\eta_3\)-direction and interpolated in \((\eta_1, \eta_2)\). The DOFs for \(V_h^2\) contain histopolation over surface elements (green) and for \(V_h^3\) one has 3d histopolation over volume elements (violet). This ultimately leads to the commuting property.

Fig. 2: Degrees of freedom (DOFs) of the commuting projectors.

Uni-variate spline spaces

Finite element spaces in Struphy are based on uni-variate B-spline spaces on the interval \(\eta \in [0, 1]\), as introduced for instance in Section 2.2 of Hughes, Cottrell, Bazilevs, “Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement”, Comput. Methods Appl. Mech. Engrg. 194 (2005). B-spline basis functions of degree \(p\) are denoted by \(N_i^p(\eta)\) (thus also called N-splines in the following). They are piece-wise polynomials with a regularity of at most \(p-1\) (they are piece-wise constants for \(p=0\)). Linear combinations af these basis functions give a spline function

\[f_h(\eta) = \sum_{i=0}^{n-1} f_i N^p_i(\eta)\,,\qquad f_i \in \mathbb R\,,\]

where the \(\mathbf f = (f_i)_{i=0}^{n-1} \in \mathbb R^n\) are the spline coefficients. The dimension of the spline space is \(n\). In Struphy two different boundary conditions can be chosen: periodic or clamped. For the latter, the first and last basis functions are “interpolatory” at \(\eta=0\) and \(\eta=1\), respectively, meaning that \(f_h(0) = f_0\) and \(f_h(1) = f_{n-1}\). Assuming that the interval \([0, 1]\) is split into \(num_elements\) cells (or elements), the dimension of the periodic space is \(n=num_elements\), whereas the dimension of the clamped space is \(n=num_elements + p\). Spline basis functions of degree \(p=2\) are potted in Fig. 3 for \(num_elements=8\) cells.

Fig. 3: Uni-variate spline spaces of degree \(p=2\).

In Fig. 3 the black crosses are the break points (cell interfaces) and the red dots are the Greville points, which are the interpolation points used to define the DOFs for the projectors in (2). Moreover, we encounter the so-called D-splines which are related to the derivatives of N-splines via

\[\frac{\partial }{\partial \eta} N_i^p = D^{p-1}_{i-1} - D^{p-1}_i\]

D-splines are rescaled B-splines of degree \(p-1\). In the clamped case the dimension \(d\) of the D-spline space is thus one less than that of N-splines; in this case \(D^{p-1}_{-1} = D^{p-1}_{n-1} = 0\), where \(n\) is the dimension of the N-spline space.

Fig. 4 shows the uni-variate spline spaces for degree \(p=1\). Note that in this case the D-Splines are piece-wise constants and a spline function is therefore generally not continuous.

Fig. 4: Uni-variate spline spaces of degree \(p=1\).

Higher degree leads to higher smoothness in Struphy, where spline functions have the maximal regularity \(p-1\). Moreover, the support of each basis function is the union of \(p+1\) adjacent cells and thus increases with the degree. This allows for high-order methods in Struphy. See Fig. 5 for an example with \(p=4\).

Fig. 4: Uni-variate spline spaces of degree \(p=4\).

Polar splines

For domains with a polar singularity at \(\eta_1 = 0\) (e.g. IGAPolarCylinder, IGAPolarTorus), DerhamOptions(polar_splines=True) builds the C1 polar spline spaces of Holderied and Possanner (JCP 2022). The first two rings of tensor-product splines around the pole are replaced by a few polar basis functions, via extraction operators \(\mathbb E\) (basis) and \(\mathbb P\) (degrees of freedom) that act in the poloidal plane \((\eta_1, \eta_2)\) only.

Hence, the inter-/histopolation matrices of the commuting projectors, \(\mathbb P \mathcal I \mathbb E^T\), and the mass matrices with weights independent of \(\eta_3\) have a 2D x 1D Kronecker structure \(\mathbb A_\textnormal{pol} \otimes \mathbb A_3\) (struphy.polar.kron_2dx1d). They are solved with a 1d direct solver in \(\eta_3\) and one of the following solvers for the sparse poloidal matrix \(\mathbb A_\textnormal{pol}\):

• "schur" (default): block elimination of the pole block (polar coefficients and first tensor-product ring), with an exact (Kronecker or sparse LU) solve of the outer block and a dense LU of the small Schur complement,

• "lu": sparse LU factorization,

• "ilu": incomplete LU factorization, used as preconditioner of an iterative solver.

The solver of the projectors is set with DerhamOptions(polar_poloidal_solver=...). For weighted mass matrices, use PolarMassMatrixPreconditioner (precond="PolarMassMatrixPreconditioner" in the propagator options); it is exact for axisymmetric weights and mappings.

Stiffness operators \(\mathbb d^T \mathbb M_{k+1} \mathbb d + \sigma \mathbb M_k\) (gradient and curl) are a sum of two Kronecker products, \((\mathbb K_\textnormal{pol} + \sigma \mathbb M_\textnormal{pol}) \otimes \mathbb M_3 + \mathbb W_\textnormal{pol} \otimes \mathbb S_3\), where \(\mathbb W_\textnormal{pol}\) is the poloidal mass matrix of the \(\eta_3\)-derivative terms and \(\mathbb S_3\) the 1d stiffness matrix in \(\eta_3\). PolarStiffnessPreconditioner inverts this sum by fast diagonalization in \(\eta_3\), with one poloidal solve per \(\eta_3\) mode. It is used automatically for precond="StiffnessPreconditioner" with polar splines. For the gradient it is exact on axisymmetric domains (without straight field line coordinates). For the curl, the components \((E_1, E_2)\) and \(E_3\) are decoupled (block Jacobi) and a correction on the gradients is added. Weighted averages over \(\eta_2\) and/or \(\eta_3\) (e.g. the flux-surface average of the gyrokinetic Poisson equation, stab_average) are available on polar 0-forms through PolarWeightedAverageProjection, and are included in the polar stiffness preconditioner (the average over \(\eta_2\) by a low-rank correction).

On the CuPy backend, extraction, derivatives, Kronecker operators and these solves keep their coefficients on the device. Factorization and mass-preconditioner assembly run once on the host; cached device factors apply subsequent solves. Identical factors are shared across polar spaces and projector/preconditioner instances of the same Derham. Matrix contents and solver settings are matched exactly; updated weights use new factors, and unused entries are released. Polar MPI collectives support CUDA-aware device buffers or host staging through cunumpy.mpi.mpi_buffer. The feectools halo exchanges require CUDA-aware MPI on real GPUs. SplineFunction.__call__ dispatches CUDA kernels for marker and dense/sparse grid evaluation. The fixed local-grid helper eval_tp_fixed_loc remains host-only; distributed field-evaluation MPI handling and coverage are separate follow-ups.