B-spline Finite Element Exterior Calculus
1. 📝 Table of contents
- 1. 📝 Table of contents
- 2. 🧮 A brief overview of the underlying mathematics
- 3. 💡 Code introduction and examples
- 3.1. One-dimensional B-splines
- 3.2. Tensor product B-spline Finite Element Exterior Calculus
- 3.2.1. Mapped domains (physical space)
- 3.3. $L^2$ projection
- 3.3.1. Solver options
- 3.3.2. Quadrature accuracy
- 3.4. Poisson problem
- 3.4.1. Adjoint differential operators
- 3.4.2. Variable coefficient Poisson problem
- 3.4.3. Custom constraints
- 3.4.4. Custom domain decomposition
- 3.4.5. Visualisation
- 3.5. Curl-curl problem
- 3.6. Eigenvalue problems
- 3.7. Lower-dimensional problems
- 3.7.1. Reduced de Rham sequence
- 3.7.2. Full de Rham sequence
- 4. 📚 References
2. 🧮 A brief overview of the underlying mathematics
3. 💡 Code introduction and examples
We give a brief overview of the most important aspects of the code. A more elaborate documentation based on Doxygen can be found here and directly compilable examples can be found under [examples/](examples) (the examples are built if the -Dbspline_feec_EXAMPLE=on CMake flag is used).
3.1. One-dimensional B-splines
The BSplineSpace type from the m_bspline module encapsulates a one-dimensional B-spline space. The space of clamped cubic B-splines on 32 equidistant subintervals can be created as follows
Using this space, we can create a B-spline function as follows
which by default has vanishing coefficients (i.e., the zero function). A B-spline function can be evaluated on the interval $[0, 1]$
or its derivative can be evaluated
More interestingly, we can compute the $L^2$ projection of any real-valued scalar function (scalar_fun) on the interval $[0, 1]$ (the coefficients of bfun now represent the B-spline approximation of scalar_fun)
and compute the $L^2$ error of this approximation
Please note that
- One-dimensional B-splines are always defined on a uniform grid on the unit interval, using either clamped or periodic boundary conditions. Non-uniformity of the grid is possible only in physical space with $m$-forms via a coordinate transformation
- One-dimensional B-splines are not MPI-aware
An example is found in examples/one_dimensional.f90.
Expected example output
3.2. Tensor product B-spline Finite Element Exterior Calculus
The finite dimensional spaces of logical $m$-forms (for $m = 0, 1, 2, 3$), that make up the de Rham sequence on the logical domain $\hat{\Omega} = [0, 1]^3$, are encapsulated in the MFormSpace type from the m_mform module, and can be created as follows
where the three BSplineSpace objects from the m_bspline module correspond to the one-dimensional B-splines used in the logical $\hat{x}, \hat{y}, \hat{z}$-directions (in this example using $24, 16, 32$ intervals and B-splines of degree $4, 6, 3$ per direction).
Using the $m$-form spaces stored in the derham array, we can create a scalar-valued $0$-form and vector-valued $1$-form as follows
which are initialised to the trivial zero-valued zero- and one-forms by default. The $0$-form can be evaluated (at $\hat{\bm{x}} = \begin{bmatrix}0.1 & 0.2 & 0.3 \end{bmatrix}$) as follows
whereas the second component of the vector-valued $1$-form is evaluated as
Similar functions exist for the evaluation of the gradient of a $0$-form, the curl of a $1$-form, and the divergence of a $2$-form; for example
yields the third component of the curl of the $1$-form form1 Alternatively, one can compute the gradient of a $0$-form (which is a $1$-form), the curl of a $1$-form (which is a $2$-form), and the divergence of a $2$-form (which is a $3$-form); for example
yields the same as the directly evaluated curl of the $1$-form.
3.2.1. Mapped domains (physical space)
Mapped domains are available by means of coordinate transformation objects, such as e.g. Cylindertransform, ToroidalTransform, FrenetSerretTransform from the m_coord_transform module. Specific user-defined coordinate transformations require the implementation of a derived object
for which user-defined functions should be implemented that evaluate the transformation (userdefined_transform), the inverse transformation (userdefined_transform_inv), the Jacobian (userdefined_jacobian), the Jacobian matrix (userdefined_jacobian_matrix), and the metric tensor (userdefined_G_matrix), as well as the inverses of the matrices. A concrete example of how to implement a user-defined coordinate transformation can be found in the examples/user_defined_transform directory.
A solid toroidal domain $\Omega$ with minor radius $0.2$ and major radius $1$ can be obtained as follows
Evaluation of physical $m$-forms (as well as the gradient/curl/divergence) requires the pushforward operator, which in turn is defined in terms of the Jacobian matrix of the coordinate transformation. This is automatically taken care of when evaluating $m$-forms when the coordinate transformation is provided as input to the evaluate function, for example
evaluates the second component of the physical $1$-form form1. The input coordinates are always assumed to be in logical space.
3.3. $L^2$ projection
Several linear PDEs can be solved 'out-of-the-box' (on a mapped domain), see also the description of the PDEs in About the project.
For each $m$-form space the corresponding $L^2$ projection problem can be solved. The right-hand side of this equation (in fact, for each of the three PDE types) involves the computation of the $L^2$ inner product of each basis function from the Finite Element space with some user-specified $m$-form (in the example below it is a scalar-valued $0$-form scalar_fun), the result of which can be stored in an $m$-form as follows
where the coordinates must be provided if the inner product is computed on a mapped domain. The left-hand side of the equation is given by the so-called mass matrix whose elements contain the inner product of each pair of two basis functions from the Finite Element space, and can be obtained as follows (in this example the $0$-form mass matrix is computed)
where we again note that the coordinates must be provided if the mass matrix is computed on a mapped domain. The Finite Element solution to the PDE (i.e., the $L^2$ projection of scalar_fun onto the space of $0$-forms) can now be obtained by solving the linear system of equations as follows
and the $L^2$ error of the approximation can be computed
An example for $1$-forms is found in examples/oneform_l2_projection.f90. This example also supports MPI (in fact all examples do, except for the one-dimensional ones):
3.3.1. Solver options
The solver is a GenericSolver object which is a high-level wrapper around PETSc, which is used for solving all of the linear systems discussed here. Several PETSc and/or hypre solvers are supported by passing an optional SolverInitOptions object to the init procedure (a clever choice is made based on the matrix type and parallelisation if no options are provided). The option ksp_type selects the method used for the outer iteration (a list of PETSc Krylov subspace methods and direct solvers), the option pc_type selects the preconditioner (a list of PETSc preconditioners), and hypre_type specifies the hypre preconditioner (a list of hypre preconditioners). For example, the BoomerAMG preconditioner of hypre can be selected in combination with the MinRes Krylov method as follows
where we note that the presence of the hypre_type option automatically sets pc_type=hypre.
Furthermore, the solve procedure supports an optional SolverSolveOptions object that specifies the absolute and relative residual tolerance (atol and rtol), the maximum number of iterations (max_iter), and verbosity (verbosity). For example
ensures that the solver always displays information about the solver's success or failure at the end of the iteration:
Supported verbosity strings are: ‘'warn never’,'warn on failure','warn always','error on failure'`.
3.3.2. Quadrature accuracy
Both the left- and right-hand sides of the PDEs involve the $L^2$ inner product, which is approximated by means of a composite Gauss-Legendre (GL) quadrature rule.
Integration of a single one-dimensional B-spline of degree $r$ is exact when using GL quadrature with $1 + \lfloor r / 2\rfloor$ nodes. For the right-hand side, we compute the inner product of a single B-spline (per logical direction) multiplied by a user-specified function as well as the Jacobian of the coordinate transformation. Depending on the specific user-specified function and Jacobian, one might require additional quadrature nodes (i.e., on top of those required for the exact integration of the single B-spline) to obtain an accurate approximation. By default, two additional quadrature nodes are used, but if so desired this number can be changed using the optional flag n_quad_extra
where n_quad_extra=0 means that $1 + \lfloor r / 2\rfloor$ nodes are used.
Integration of a product of two one-dimensional B-splines of degree $r$ is exact when using GL quadrature with $1 + r$ nodes. For the left-hand side, we compute the inner product of a product of two B-splines (per logical direction) multiplied by the Jacobian of the coordinate transformation. Depending on this Jacobian, one might require additional quadrature nodes here as well. The number of additional quadrature nodes can be specified by the optional flag n_quad_extra
where n_quad_extra=0 means that $1 + r$ nodes are used.
3.4. Poisson problem
The Poisson problem is formulated for the $0$-forms. The computation of the right-hand side is exactly as shown for the $L^2$ projection problem. For the left-hand side, however, we require the stiffness matrix rather than the $0$-form mass matrix; the stiffness matrix can be constructed as follows
Well-posedness of the Poisson problem requires the inclusion of boundary conditions, which are imposed via the so-called MFormDirichlet object, which is an object that imposes a Dirichlet boundary condition on the non-periodic boundaries for obtaining a well-posed linear system of equations, based on the provided $m$-form space. At the initialisation stage of solver, the stiffness matrix and the constraint object(s) are combined into a single constrained matrix
The constrained linear system can be solved as usual
where the $0$-form sol is defined on the entire unconstrained tensor product B-spline space, but is guaranteed to satisfy the provided constraint(s).
Note that constraints can be used with any matrix type, and thus also with a mass matrix. This results in a constrained $L^2$ projection.
An example in logical space is found in examples/divgrad_logical.f90.
3.4.1. Adjoint differential operators
We now consider a slightly more specific Poisson problem, where we want to find the longitudinal (i.e., curl-free) part of a $1$-form $\bm{f}$: find $\varphi \in V^0$ such that $$ \langle \nabla \varphi, \nabla \phi\rangle = \langle \bm{f}, \nabla \phi\rangle $$ for all $\phi \in V^0$. The right-hand side can now most easily be obtained by first computing the $L^2$ inner product of $\bm{f}$ with all $1$-forms from the space $V^1$ (where f_x, f_y, f_z are use-defined functions implementing the components of $\bm{f}$)
followed by computing the gradient adjoint (i.e., the negative divergence) of this result
An example is found in examples/divgrad_physical.f90.
3.4.2. Variable coefficient Poisson problem
Thus far, the optional left-hand side coefficient (matrix) was ignored and thus reduced to the default identity operator $\mathbf{M} = \mathbf{I}$. We extend the previous example by including a non-trivial coefficient matrix $\mathbf{M} = \mathbf{\Pi}_\perp$, where $\mathbf{\Pi}_\perp$ is an orthogonal projection operator that projects orthogonal to the third physical coordinate direction: $$ \mathbf{\Pi}_\perp \bm{f} = \begin{bmatrix}f_x\f_y\0\end{bmatrix}. $$
In order to proceed, we need to understand how the stiffness matrix is constructed. Let $\mathbf{L}^0$ denote the $0$-form stiffness matrix, $\mathbf{G}$ denote the gradient matrix (a sparse matrix with exactly one $+1$ and $-1$ per row), and $\mathbf{M}^1$ denote the $1$-form mass matrix. The stiffness matrix is constructed as follows (this is what is done in the init procedure of DiffDiffMatrix): $$ \mathbf{L} = \mathbf{G}^T \mathbf{M}^1 \mathbf{G}. $$ The $1$-form mass matrix can optionally be passed to the init procedure
which is efficient if the mass matrix is available already, but this functionality is required for the use of variable coefficients, as these coefficients are part of the mass matrix. Note that in this case the coordinate transformation need not be passed to the init procedure of DiffDiffMatrix, as the coordinate transformation is already taken into account in the assembly of the $1$-form mass matrix. In this example, the $1$-form mass matrix should be constructed as follows
where perp_projector is a user-defined function with the following signature
where s1, s2 denote the matrix indices.
Note that scalar-valued coefficients $\mathbf{M}(\hat{x},\hat{y},\hat{z}) = m(\hat{x},\hat{y},\hat{z}) \mathbf{I}$ are supported using userfun=scalar_fun.
An example is found in examples/divgrad_perp.f90.
3.4.3. Custom constraints
Consider the following example, where we combine a Dirichlet boundary condition, as well as a first-order regularity constraint for disc-like coordinate singularities (by default, the regularity degree is taken maximally)
Please note that both the Dirichlet and regularity constraints are available for each of the $m$-form spaces using the same syntax. The following homogeneous Dirichlet boundary conditions are provided by the MFormDirichlet object on the non-periodic part of the boundary $\partial \hat{\Omega}$:
- For $0$-forms: $\hat{\phi} = 0$
- For $1$-forms: $\hat{\bm{\phi}} \times \hat{\bm{n}} = \bm{0}$
- For $2$-forms: $\hat{\bm{\phi}} \cdot \hat{\bm{n}} = 0$
- For $3$-forms no boundary conditions are imposed
3.4.4. Custom domain decomposition
Thus far, the MPI parallelisation has been completely hidden behind the high-level interface, and is taken care of automatically. The parallelisation is based on a domain decomposition of the logical domain $\hat{\Omega}$ into subdomains, where each MPI rank is responsible for one subdomain. Several options are available to customise the domain decomposition. By default a hybrid parallelisation is used (using pure MPI), where a distinction is made between the ranks that share the same memory space (i.e., the ranks that run on the same compute node) and those that do not. Per logical direction one can decide how the work should be distributed: either sequentially (i.e., a single rank is responsible for the entire direction), using shared memory parallelism (i.e., all ranks that share the same memory space are responsible for the direction), or using distributed memory parallelism (i.e., compute nodes are responsible for the direction).
To customise the domain decomposition, one has to pass a DomainDecomp object to the construction function of the de Rham sequence.
We explain the input parameters by means of several examples (assuming a 3D problem).
- For a multi-node cluster, the default behaviour is to use shared memory parallelism in the logical $\hat{y}$-direction and distributed memory parallelism in the logical $\hat{z}$-direction, which can be set manually as follows: Hence, if $8$ nodes are used, where each node has $16$ shared memory ranks, then the logical $\hat{y}$-direction is decomposed into $16$ subintervals (one per shared memory rank) and the logical $\hat{z}$-direction is decomposed into $8$ subintervals (one per node), resulting in a total of $16 \times 8 = 128$ MPI ranks.domain = domaindecomp(memory_layout_y='shared', memory_layout_z='distributed')
- If so desired, the hybrid parallelisation can be disabled and pure distributed memory parallelism can be used instead (so called flat MPI): which defaults to distributed memory parallelism in the logical $\hat{z}$-direction. This incurs a memory overhead due to the use of ghost layers and also results in more MPI communication.domain = domaindecomp(flat_mpi=.true.)
- A 3D domain decomposition can be specified as follows: where a balanced factorisation of the number of shared memory ranks is used to determine the number of subintervals in each direction. Hence, if $8$ nodes are used, where each node has $16$ shared memory ranks, then the logical $\hat{x}$- and $\hat{y}$-directions are both decomposed into $4$ subintervals (resulting in $4 \times 4 = 16$ subdomains for the $16$ shared memory ranks) and the logical $\hat{z}$-direction is decomposed into $8$ subintervals (one per node), resulting in a total of $4 \times 4 \times 8 = 128$ MPI ranks.domain = domaindecomp(memory_layout_x='shared', memory_layout_y='shared', memory_layout_z='distributed')
- If so desired, the number of subintervals in each direction can be specified manually: This can be combined with the memory layout options as well. Here, it must be ensured that the product of the number of subintervals per memory layout matches the number of ranks in that memory layout (i.e.,domain = domaindecomp(nx=2, ny=4, nz=8, memory_layout_x='shared', memory_layout_y='shared', memory_layout_z='distributed')
Nx * Nymust equal the number of shared memory ranks per node, andNzmust equal the number of nodes). - The construction function of the
DomainDecompallows for the specification of the MPI communicator:which allows for the MPI parallelisation to take place on a subset of the available ranks (the default value isdomain = domaindecomp(nx=nr_shmem_ranks, nz=nr_nodes, comm=petsc_comm_world)PETSC_COMM_WORLD).
Many of the input parameters of the DomainDecomp can be combined, consistency is checked at runtime. Any missing parameters are determined automatically: the outer dimensions are prefered for distributed memory parallelism, whereas the inner dimensions are prefered for shared memory parallelism. If multiple directions have the same memory layout, then a balanced factorisation is used to determine the number of subintervals in each direction.
Please note that the domain decomposition in each direction is constrained by conditions found in the check_decomposition_1d function. The following simplified conditions are sufficient (per direction, assuming default CMake flags are used)
$$M = 1 \quad \lor \quad r+1 < M \le \frac{N}{r+2}$$ where $r$ is the B-spline degree, $N$ denotes the number of intervals and $M$ is the number of 1D subdomains in some direction (resulting in approximately $N/M$ intervals per subdomain, where we note that $N/M$ need not be integer). If the regularity constraint is used, then there should be no distributed domain decomposition in the logical $\hat{y}$-direction and the first 1D subdomain in the logical $\hat{x}$-direction should have at least $1 + r$ B-splines.
Advanced options
As mentioned in the advanced installation options, one can alter the MPI communication options by means of CMake flags. We let $\mu$ denote the size of the guard layer (1 by default) and $\eta$ denotes the number of MPI neighbours (1 by default), then the constraint on the number of 1D subdomain intervals becomes:
$$M = 1 \quad \lor \quad r+1 < M \le \frac{N \eta}{r+\mu+\eta}$$
3.4.5. Visualisation
For visualisation the VTKFortran library is required, which is automatically downloaded and compiled by means of the -Dbspline_feec_PARAVIEW=on CMake flag. The ParaView object allows for the initialisation of a sequence of .vtu files to be written to file, with a user-provided base name (in this case divgrad).
Then, MFormFun objects can be passed to the write procedure (which also accepts arrays of $m$-forms)
resulting in the divgrad00000.vtu file which can be visualised using ParaView.
3.5. Curl-curl problem
The curl-curl problem is formulated for $1$-forms. Here, we consider the specific problem of finding the transversal part (i.e., the divergence free part) of some $2$-form $\bm{f}$: find $\bm{\varphi} \in V^1$ such that $$ \langle \nabla \times \bm{\varphi}, \nabla \times \bm{\phi}\rangle = \langle \bm{f}, \nabla \times \bm{\phi}\rangle $$ for all $\bm{\phi} \in V^1$. Analogous to the way the right-hand side of the Poisson problem for the longitudinal part of a $1$-form was constructed, we construct the right-hand side in two steps
The curl-curl stiffness matrix again can be obtained by making use of a DiffDiffMatrix object
which (as usual) can be combined with a constraint object upon initialising the linear solver
The curl-curl problem is a notoriously difficult problem to solve due to the presence of the large null-space of the curl operator. Hence, most preconditioners provided by PETSc will fail or be ineffective. However, hypre provides a special linear solver4 for this class of linear problems: the Auxilliary-space Maxwell Solver (AMS). For the DiffDiffMatrix of the $1$-form space, the init procedure of the solver automatically selects this preconditioner, otherwise it can be selected manually using the optional hypre_type=ams flag.
After solving the linear system of equations, the transversal part of $\bm{f}$ is obtained by computing the curl of the resulting solution:
An example is found in examples/curlcurl.f90.
3.6. Eigenvalue problems
Eigenvalue problems can be solved using the EigenSolver object from the m_eigensolver module. As with the GenericSolver, the EigenSolver is a high-level wrapper around SLEPc. Any matrix of AbstractMatrix type (MassMatrix or DiffDiffMatrix) can be used for defining the eigenvalue problem, and similarly any two such matrices can be combined to form a generalized eigenvalue problem: $$ \mathbf{A} \bm{u} = \lambda \mathbf{B} \bm{u}. $$ In particular, the Laplace eigenvalue problem as stated in the README.md can be set up as follows
which initilizes the Arnoldi method to compute the eigenvalues with smallest real part of the generalized eigenvalue problem defined by the stiffness and mass matrices of the $0$-form space. Note that constraints can be used here as well, for example to impose regularity at coordinate singularities. The eigenvalue problem can then be solved as follows
which results in an array eigenvals_r containing the real parts of the computed eigenvalues. Depending on the solver convergence, fewer than nr_pairs eigenvalues might have been computed. The corresponding eigenfunctions can be obtained as follows
A complete example is found in examples/divgrad_eigen.f90.
3.7. Lower-dimensional problems
The functionality described thus far is also available for one- and two-dimensional problems.
3.7.1. Reduced de Rham sequence
We consider the two-dimensional de Rham sequence resulting in three $m$-form spaces, where the first and last space are scalar-valued and the middle space is vector-valued with two components. For example, the two-dimensional de Rham sequence on the unit disk can be created as follows
on which we can solve the Poisson problem for $0$-forms or the curl-curl problem for $1$-forms as before. Examples are found in examples/divgrad_2d.f90 and examples/divgrad_1d.f90.
There is no explicit implementation for two-dimensional problems, instead the three-dimensional implementation is used with one constant B-spline in the logical $\hat{z}$-direction for the $0$-form space. Several things to note here:
- The default behaviour for the MPI domain decomposition is to decompose in the logical $\hat{y}$-direction only, but this can be customised using a
DomainDecompobject as described previously, where it must be ensured thatNz=1. - User-provided functions for the right-hand side (and variable coefficients) must have the correct signature for three-dimensional functions, i.e., they must accept three logical coordinates as input arguments (the logical $\hat{z}$-coordinate can simply be ignored).
- There is no specific implementation for a two-dimensional
CoordTransform, also here the three-dimensional implementation is used where the third logical and physical coordinates are ignored. - The $1$-form space in two dimensions has only two components, but the implementation still uses three components where the third component is an empty space.
- The $2$-form space in two dimensions is a scalar-valued space, but the implementation still uses three components where the first two components are empty spaces, hence the scalar-valued 2D $2$-form should be evaluated as follows at $\hat{x} = 0.1, \hat{y} = 0.2$ val = evaluate(form2, 3, .1_wp, .2_wp, .0_wp, coord_transform=coords)
3.7.2. Full de Rham sequence
Alternatively, the full de Rham sequence can be obtained where each of the components varies only in the logical $\hat{x}$- and $\hat{y}$-directions (or in the one-dimensional case, only in the logical $\hat{x}$-direction). This functionality is available by means of the DeRhamSequence3C functions (3C stands for 'three components'). For example, the de Rham sequence with components varying only in the logical $\hat{x}$-direction can be obtained as follows
As with the reduced de Rham sequence, the implementation is unaware of the dimensionality of the problem, hence any user-provided functions, coordinate transformations, and evaluation of $m$-forms should be done as if the problem is three-dimensional.
4. 📚 References
1 Arnold DN, Falk RS, Winther R. Finite element exterior calculus, homological techniques, and applications. Acta Numerica. 2006;15:1-155. doi:10.1017/S0962492906210018
2 Remmerswaal R (2025). Higher-order FEEC discretization on a solid toroidal domain. Not yet available
3 De Boor C (1978). A practical guide to splines. Springer-Verlag New York
4 Hiptmair R, Xu J (2007). Nodal auxiliary space preconditioning in $H(\text{curl})$ and $H(\text{div})$ spaces. SIAM Journal on Numerical Analysis, 45(6), 2483-2509. doi:10.1137/060660588
