Isogeometric discrete differential forms: Non-uniform degrees, Bezier extraction, polar splines and flows on surfaces Non-uniform degrees, Bézier extraction, polar splines and flows on surfaces

Spaces of discrete differential forms can be applied to numerically solve the partial differential equations that govern phenomena such as electromagnetics and fluid mechanics. Robustness of the resulting numerical methods is complemented by pointwise satisfaction of conservation laws (e.g., mass conservation) in the discrete setting. Here we present the construction of isogeometric discrete differential forms, i.e., differential form spaces built using smooth splines. We first present an algorithm for computing Bézier extraction operators for univariate spline differential forms that allow local degree elevation. Then, using tensor-products of the univariate splines, a complex of discrete differential forms is built on meshes that contain polar singularities, i.e., edges that are singularly mapped onto points. We prove that the spline complexes share the same cohomological structure as the de Rham complex. Several examples are presented to demonstrate the applicability of the proposed methodology. In particular, the splines spaces derived are used to simulate generalized Stokes flow on arbitrarily curved smooth surfaces and to numerically demonstrate (a) optimal approximation and inf–sup stability of the spline spaces; (b) pointwise incompressible flows; and (c) flows on deforming surfaces.


Introduction
Partial differential equations (PDEs) describing physical phenomena are built on a rich differential and geometric foundation of conservation laws, topological constraints, symmetries and invariants. The reliability of numerical algorithms that are used to discretize and approximately solve these equations is of the utmost importance for countless scientific and engineering applications, and this is intimately connected to the differential and geometric structure underlying the PDEs. Specifically, for physical problems such as electromagnetism and incompressible fluid flows, consistent, stable and accurate numerical methods that ensure physical fidelity of the discrete solutions can be built by mimicking the structure underlying the continuous problem (e.g., the identities div-curl= 0 and curl-grad= 0) at the discrete level. The formulation of such numerical methods is our focus in this article, with a special attention towards high-order accurate discretizations of PDEs defined on surfaces in R 3 .
The development of discretization methods that aim to mimic symmetries and invariants at the continuous level is an active area of research. Some of the significant contributions in this area have come in the form of mimetic finite difference methods [1,2], mimetic spectral element methods [3,4], discrete exterior calculus [5,6], finite element exterior calculus [7,8], and physics-compatible or structure-preserving isogeometric analysis [9][10][11][12]. These methods have one thing in common: they are driven by geometric interpretations of the solution fields. In particular, the solution fields are interpreted as differential forms [13], which are objects that are naturally associated to geometric objects of different dimensions. For example, for fluid flow on an d-dimensional domain, velocities may be interpreted as fluxes that flow through codimension-1 geometries (i.e., as differential forms of order d − 1) and their divergence as the mass being lost or produced in d-dimensional geometries (i.e., as differential forms of order d).
Differential forms provide a compact, clear and intuitive language for discussing both PDEs as well as their discretization. They are particularly helpful in identifying which parts of the PDE are conservation laws that do not depend on any notion of a metric, and which parts are constitutive laws -the discretizations are then constructed to exactly satisfy the former and accurately approximate the latter. The framework of finite element exterior calculus [7,8] is based on precisely this formalism. It has led to a unified approach for developing accurate finite element differential form spaces and analysing stability and well-posedness of the discrete problems. This article focuses on the methods that belong to the extension of finite element exterior calculus to isogeometric analysis.
Isogeometric analysis [14] relies on the use of smooth splines, i.e., smooth piecewise-polynomial functions, for building both the geometry on which the problem is defined, as well as the discrete finite element spaces used to solve the problem. The last decades have seen the extensive use of splines for numerical solutions of challenging problems such as the design and optimization of wind turbines [15], development of cardiac devices [16], and multiphase flows [17][18][19] and fracture dynamics [20,21] governed by high-order phase field theories. Several recent developments [22][23][24] have provided the theory supporting the large body of numerical evidence that smooth splines demonstrate better approximation behaviour per degree of freedom than less smooth or non-smooth spline spaces (e.g., traditional C 0 finite element spaces). The extensions of finite element exterior calculus to isogeometric analysis have come via the development of isogeometric discrete differential form spaces, i.e., discrete differential form spaces built using smooth splines. These isogeometric discrete differential forms are used to solve PDEs on domains that are also built using smooth splines. Examples include discrete differential forms on rectangular or cuboidal domains built using tensor-product splines [10,11,25,26] and adaptively-refined splines [12,27,28], with applications to electromagnetism and incompressible flows being common.
The above represents the state of the art on spline differential forms and the existing literature does not address the problem of simulating PDEs on arbitrary surfaces with smooth spline-based differential forms. While significant advances have been made in the understanding of spline differential forms on (locally refined) quadrilateral meshes, unstructured meshes are needed for building arbitrary surfaces. In particular, there are two types of unstructured meshes that need to be studied -ones where the number of quadrilaterals meeting at an interior vertex is not equal to 4 (excluding T-junctions), and ones where the quadrilaterals degenerate into triangles. One of the contributions of this article is taking the first steps that address the latter class of unstructured meshes; such meshes are called polar meshes and the splines built on them are called polar splines.
More broadly, the motivation for this article is construction of isogeometric discrete differential forms for numerical approximation of (scalar and vector) solutions to PDEs. We focus on discrete differential form spaces built using two particular classes of nonstandard spline spaces -univariate multi-degree splines and bivariate polar splines. Multi-degree splines [29][30][31] are splines that allow local polynomial degree adaptivity, and polar splines [29,32] are non-tensor-product splines that allow the construction of singularly parametrized, genus 0 surfaces that are nevertheless C k smooth. For instance, the geometry on the right is a C 1 polar spline surface that is singularly parametrized -it is built from a rectangular domain by collapsing a pair of opposite edges to points.
Solving PDEs on such surfaces has many applications; e.g., numerical weather prediction [33]. Of particular interest to us is the study of biological fluidic membranes such as lipid bilayers [34][35][36]. These membranes can be thought of as the envelopes for eukaryotic cell contents. These are versatile structures that behave as in-plane viscous fluids and out-of-plane solids. Computationally studying the behaviour of such structures requires the ability to simulate two-dimensional viscous fluid flow on a curved, evolving surface. Several recent methods have been proposed to solve such problems; e.g., using discrete exterior calculus [37], trace finite elements [38,39] and surface finite elements [40][41][42][43]. These methods are built using functions of low-regularity that are either C 0 or discontinuous. The methods in [37][38][39] are low-order methods that use piecewise polynomials of degree at most 1, and the methods in [38,39] are for surfaces defined implicitly on a background mesh. The method in [41] is high-order accurate but needs Lagrange multipliers to impose tangentiality of the fluid velocity on a curved surface; penalties are used to approximately achieve the same or to enforce approximate conservation in [38,39,43].
In this paper, we develop novel isogeometric discrete differential forms that, in particular, offer a high-order alternative to the above methods for simulation of flows on smooth surfaces without any recourse to Lagrange multipliers or penalties for enforcing tangentiality of the flow. Section 2 presents the mathematical preliminaries needed for our approach. We subsequently discuss the theoretical and algorithmic aspects behind the construction of multidegree spline differential forms ( Sections 4 and 3.2 ), and their application to building polar spline differential forms (Section 5). In particular, we show how this enables us to mimic the cohomological structure of the de Rham complex at the discrete level. We demonstrate the high-order accuracy, stability and applicability of the discrete differential form spaces by simulating, in particular, generalized Stokes flow on fixed and deforming smooth surfaces (Section 6). The spaces also allow us to simulate pointwise incompressible tangential flows on surfaces. See an example of such pointwise incompressible tangential flow on the right where streamfunction contours and tangential velocities are displayed; see Section 6.3 for details.

Mathematical preliminaries
Let us start by presenting some exterior calculus preliminaries and, in particular, introduce the L 2 de Rham complex. We follow the presentation of [8] in an abbreviated form. Moreover, since we are interested in building spline differential forms on smooth 2-manifolds, Ω , we restrict the following discussion to the two-dimensional setting. Note that we only present the most basic relations in this section; other necessary notation and material will be presented when needed.

Outlook
As mentioned in Section 1, the motivation for this article is the construction of stable and high-order accurate spline-based finite element methods for numerical approximation of (scalar and vector) solutions to PDEs on smooth 2-manifolds Ω . We do so within the conceptual framework of finite element exterior calculus [7,8] and its splinebased counterparts [10,11,26]. As shown in [8], for instance, well-posed problems can be formulated at the discrete level if the finite element spaces form a subcomplex (Section 2.2) of the de Rham complex of differential forms (Section 2.3). The scalar and vector fields that solve the desired PDEs can be thought of as proxies for differential forms, and well-posedness of the continuous problems is implied by properties of the de Rham complex. Then, a first step in the construction of stable methods is the construction of a finite element subcomplex that mimics the properties of the continuous de Rham complex. The next two subsections introduces Hilbert cochain complexes and the de Rham complex of differential forms.

Hilbert cochain complexes
Let V denote a sequence of Hilbert spaces : The connecting maps d (i) are called the differentials of the complex V. Moreover, V is called bounded if its differentials are bounded linear operators, and it is called closed if the image of d (i) is closed in V (i+1) for all i. The composition property of the differentials implies that the following containment holds, Members of V (i) in ker d (i) are called i-cocycles or closed, and the members of V (i) in im d (i−1) are called icoboundaries or exact. The i-th cohomology space associated to the complex V, H i (V), is defined as the following quotient, Note that, for defining H 0 (V) and H 1 (V), the beginning and the end of the complex are automatically augmented with zero maps d (−1) := 0 =: d (2) . The cohomology space H i (V) measures the extent to which the equality in Eq. (2) fails to hold.
Given two complexes V = (V, d V ) and W = (W, d W ), linear maps f (i) : V (i) → W (i) of degree 0 are called cochain maps if they commute with the differentials for all i, Cochain maps preserve closed and exact forms and, consequently, induce maps between cohomology spaces of the two complexes, f * ,(i) : V by restriction, then the complex W is called a subcomplex of V. In this case, the inclusion ι (i) : is a cochain map from W to V and induces a natural map between their cohomologies. If, additionally, there exists a cochain projection map from V to W, it induces a surjection of cohomologies. In particular, the dimensions of H i (W) are then bounded from above by those of H i (V) for all i.
Remark 2.1. In the following, to unburden the notation, we will drop the superscripts of all differentials as it will always be clear from the context which differential is being used.

The de Rham complex of differential forms
Given a (sufficiently) smooth 2-manifold Ω ⊂ R d , d = 2, 3, let T y Ω denote the 2-dimensional tangent space at y ∈ Ω . A smooth differential i-form f , i = 0, 1, 2, on Ω is a smooth field such that f y is a real-valued skew-symmetric i-linear form on T y Ω × · · · × T y Ω . Let Λ (i) (Ω ) denote the space of all smooth i-forms, i = 0, 1, 2.
For i = 0, 1, 2, and f ∈ Λ (i) (Ω ), the exterior derivative is a linear map of degree +1, d : In local (curvilinear) coordinates x = (x 1 , x 2 ) on Ω such that y = y(x), the differential forms and the action of d are simply where ∂ denotes the boundary operator. In other words, the exterior derivative can be thought of as the dual of the boundary operator with respect to the natural duality pairing of i-forms with i-dimensional submanifolds.

The univariate spline complex
In this section we present preliminary concepts about smooth polynomial splines defined on a partition of an interval, Ω := [a, b] ⊂ R. In particular, we will allow the spline pieces to have different polynomial degrees, thereby introducing the concept of multi-degree spline spaces. We also present a set of basis functions for such spaces called multi-degree B-splines (or MDB-splines) and list some of their properties. Classical B-splines are a special case of MDB-splines.

Definition of the complex
In this preliminary section, we recall the simplest univariate spline complex on Ω that we can consider. For p ∈ N, let P p be the vector space of polynomials of degree ≤ p. Then, the simplest spline space on Ω consists only of global polynomials, Choosing S −1 p as the space of 0-forms and S −1 p−1 as the space of 1-forms, the univariate polynomial complex is defined as It can be easily verified that the exterior derivative is a surjection from Λ (0) P onto Λ (1) P , thereby yielding H 1 (G) = 0. On the other hand, H 0 (G) or the nullspace of d is one dimensional and contains constants, i.e., H 0 (G) = R.

Basis for discrete differential forms
We choose the Bernstein-Bézier polynomials B i, p , i = 0, . . . , p, as the preferred basis for 0-forms; these are defined as and they span S −1 p . Therefore, f ∈ Λ (0) P can be represented as a linear combination of the B i, p with some coefficients Furthermore, if f ∈ Λ (0) P is a 0-form, then the 1-form g := d f ∈ Λ (1) P has the representation where the B i , i = 0, . . . , p − 1, are scaled Bernstein-Bézier polynomials of degree p − 1, These are chosen as the preferred basis for discrete 1-forms. Doing so helps us define a discrete representation of the exterior derivative d in the form of the sparse matrix D (0) k of size k × (k + 1), k ≥ 2, Indeed, following Eq. (15) and arranging the coefficients g i and f i in column vectors g and f , respectively, we see that D (0) p acts on the coefficients of the 0-form (with respect to the 0-form basis B i, p ) and yields the coefficients of its exterior derivative (with respect to the 1-form basis B i, p ),

Degree of freedom interpretation
We can give a geometric interpretation to Eq. (18) using a particular one-dimensional mesh. Let 0 = γ 0 < γ 1 < · · · < γ p = 1 partition the unit interval [0, 1], and consider the corresponding one-dimensional cell complex with vertices γ i , i = 0, . . . , p, and edges τ i = γ i γ i+1 , i = 0, . . . , p − 1. We orient this complex by choosing the oriented boundary of each edge τ i to be ∂( Then, we can interpret f ∈ Λ (0) P and g ∈ Λ (1) P as cochains, i.e., linear functionals on the edges and vertices, f : Doing so, we see that the preferred 1-and 0-form basis functions B i, p and B i, p , respectively, are cochain interpolants. Moreover, we see that our discrete representations mimic the continuous version of the Stokes theorem. Indeed, extending the maps in Eq. (19) to formal sums of edges and vertices via linearity, we see that since our choice of orientation makes ( D (0) p ) T the discrete representation of the boundary operator that maps edges to oriented sums of vertices. The correspondence with the Stokes theorem is now clear.
Thus, our choice of the preferred 1-and 0-form basis functions leads to a discrete version of the Stokes theorem. This not only makes it easy to compute the degrees of freedom of an exact 1-form using Eq. (18), but in higher dimensions it will also help us exactly enforce d • d = 0 at the discrete level by a judicious choice of the discrete exterior derivatives. Finally, we will graphically depict the action of the discrete exterior derivative on the degrees of freedom as in the figure below. That is, the 0-form degrees of freedom are associated to the oriented zerodimensional cells of the mesh (i.e., vertices), and the action of D (0) p yields new degrees of freedom associated to the oriented one-dimensional cells of the mesh (i.e., edges) (see Fig. 1).

Multi-degree splines
We start by partitioning the interval Ω into a finite number of points (called breakpoints) and subintervals (called elements); the space of polynomial splines on Ω will be defined with respect to this partition.  The m + 1 strictly increasing real numbers x i , such that a =: x 0 < x 1 < · · · < x m := b, will be called breakpoints that partition Ω . The breakpoints define the intervals Ω i called elements, Next, we pick polynomial degrees p k ∈ N, k = 1, . . . , m, associated to each element Ω k . We also choose a nonnegative order of smoothness r k ∈ Z ≥0 , k = 1, . . . , m − 1, for each breakpoint of the partition, and r 0 = r m ∈ Z ≥−1 . We will distinguish between the following two cases, Non-periodic setting: r 0 = r m = −1 , Periodic setting: All the p k and r k are arranged in vectors p and r, respectively. Before proceeding, we place the following mild compatibility assumption on the chosen degrees and orders of smoothness.

Assumption 1
Each "Assumption" introduced will hold for the entirety of the document following it.

Assumption 2: Degree-smoothness compatibility
For all 1 ≤ k ≤ m − 1, we assume that Definition 3.2 (Multi-degree Spline Space). Given degree and smoothness distributions, we define a spline space on Ω as Moreover, when r 0 = r m ≥ 0, a periodic spline space on Ω is defined as Remark 3. 3. In Sections 3.2.2 and 3.2.3 , the non-periodic setting will be discussed, and the periodic setting will be discussed in Sections 3.2.5 and 3.2. 5 .
The dimension formulas for S and S per can be derived in a multitude of ways [30,31,44]. With n and n per denoting their respective dimensions, we have where an empty-sum is taken to be zero.

Assumption 3: Periodic degree-smoothness compatibility
With n and n per as defined in Eq. (25), when r 0 ≥ 0 we assume that

The non-periodic setting: Definition of the complex
To build the multi-degree spline complex, we will use the multi-degree spline spaces S and S − := S r − p − , where we define vectors p − and r − such that Then, the spaces of 0-and 1-forms are respectively chosen as follows, yielding the multi-degree spline complex It can once again be easily verified (e.g., using (23)) that the exterior derivative is a surjection from Λ

The non-periodic setting: Basis for discrete differential forms
Here, we will choose the so-called multi-degree B-splines (MDB-splines) as the preferred basis for the 0-forms. MDB-splines are a multi-degree generalization of the classical B-splines and the properties of the former mirror those of the latter; e.g., see [30,44]. Let us denote the set of MDB-splines that span S with {N i : i = 0, . . . , n −1}. See Appendix A for a recursive definition of MDB-splines using integral relations. The following set of properties are relevant for us; proofs of the same and other properties can be found in [30,31,44], for instance.    . . , N µ(k)+ p k are supported on Ω k , and they span P p k , where (e) End-point interpolation: For any 0 ≤ r ≤ p 1 , only N 0 , . . . , N r have non-zero r -th derivatives at x = a. Similarly, for any 0 ≤ r ≤ p m , only N n−r −1 , . . . , N n−1 have non-zero r -th derivatives at x = b. In particular, N 0 (a) = N n−1 (b) = 1.
Appendix B presents an algorithmic computation of MDB-splines that is much more efficient for computations than the recursive definitions from Appendix A. The algorithm computes a multi-degree extraction H [29,30] that helps express N i on Ω k as a linear combination of Bernstein-Bézier polynomials on Ω k . We briefly explain this construction here.
For k = 1, . . . , m, we denote by B Ω k j, p k , j = 0, . . . , p k , the Bernstein-Bézier polynomials of degree p k defined on Ω k ; see Eq. (14). Then, we extend them outside of Ω k by 0 and relabel them as Next, arrange these relabelled basis functions in a single vector B of length θ (m). Then, the multi-degree extraction H output by the algorithm in Appendix B [45] helps build the MDB-splines using the following expression [30], where N is the vector containing all MDB-splines N i . In particular, for k = 1, . . . , m, we call H Ω k the element extraction; it is the square submatrix of H of size ( p k + 1) × ( p k + 1) such that The matrices H and H Ω k have properties that mirror those of MDB-splines as presented in Proposition 3.4; e.g., see [30]. With the above choice of the 0-form basis, we now outline how a preferred basis for the space of 1-forms can be constructed. Note that, in general, this preferred basis will not be the same as the MDB-splines for the space S − . The following mimics the exposition from Section 3.1. 2. Let {N i : i = 0, . . . , n − 1} denote the set of preferred basis functions that span S − ; note that Eq. (25) implies that Since this is the set of preferred basis functions, it means that for a 0- The following shows how the basis functions N i can be defined element-wise; its proof is presented in Appendix C. Proposition 3. 6.
where C is a lower-triangular matrix of size p k × ( p k + 1) with all entries equal to −1. All N i , i = 0, . . . , n − 1, are, moreover, locally linearly independent and form a basis for S − .

The non-periodic setting: Degree of freedom interpretation
Similarly to the discussion in Section 3.1.3 focused on the polynomial complex, we can give a geometric interpretation to the non-periodic multi-degree spline complex using Eq. (32). Let 0 = γ 0 < γ 1 < · · · < γ n−1 = 1 partition the unit interval [0, 1], and consider the corresponding one-dimensional cell complex with vertices γ i , i = 0, . . . , n − 1, and edges τ i = γ i γ i+1 , i = 0, . . . , n − 1. We again orient this complex by choosing the oriented boundary of each edge τ i to be ∂( Once again, we interpret 0-and 1-forms as linear functionals on the cell complex, and this leads to a discrete representation of the Stokes theorem. The discussion is exactly as in Section 3.1.3, therefore we do not repeat here. Instead, we only present graphical representation showing the action of the discrete exterior derivative on the 0-form degrees of freedom (see Fig. 2).

The periodic setting: Definition of the complex
In the periodic setting, i.e., with r 0 ≥ 0, we identify the right endpoint of Ω with its left endpoint, a ≡ b. We will denote this domain with Ω per and note that the end-point identification makes it a topological circle. Let us now build the multi-degree spline complex on this periodic domain; the developments are very similar to the exposition thus far in Sections 3.2.2 and 3.2. 3.
The multi-degree spline complex in the periodic setting is built by choosing the spaces of 0-and 1-forms as where S −,per is the periodic analogue of S − with C r 0 −1 smoothness enforced between the identified ends of Ω per . Then, the periodic multi-degree spline complex is given by Once again, H 0 (M per ) or the nullspace of the exterior derivative contains only constants in Λ (0) M . However, H 1 (M per ) is non-trivial and one-dimensional here, mirroring the non-trivial topology of the periodic domain Ω per . Indeed, constants are in Λ (1) M but are not in the image of d. Thus, both H 0 (M per ) and H 1 (M per ) are isomorphic as vector spaces to R.

The periodic setting: Basis for discrete differential forms
We choose periodic MDB-splines {N per i : i = 0, . . . , n per − 1} as the basis for the 0-form space Λ (0) M . These can be computed starting from the MDB-splines for the non-periodic space S [45]. In particular, we can compute a matrixH using Algorithm 3 from Appendix D such that Note that, when working in the periodic setting, all indices (of basis functions, elements, etc.) are treated in a circular fashion here. That is, if we write "N per i ", the subscript is to be understood as below, Periodic MDB-splines have the same set of properties (except end-point interpolation unless r 0 = 0) as their non-periodic counterparts, and these are summarized in the following result; these properties can be derived from the properties of H and the properties ofH as shown in Proposition D.1.
In particular, combining Eqs. (30), (31) and (35), we can write element-local representations of the periodic MDB-splines using element extraction operators, The extraction matrices H per and H Ω k ,per have properties that again mirror those of the periodic MDB-splines; see Proposition D.1. Since this is the set of preferred basis functions, it means that for 0- 3. The degree of freedom interpretation for the multi-degree complex can be extended to the periodic complex M per . The discrete exterior derivative D (0),per n per corresponding to the choice of the periodic preferred bases has a simple action as shown above.
where the coefficients g i are now obtained from f i by the action of the periodic discrete exterior derivative D is defined as the following matrix of size k × k, The following shows how the basis functions N per i can be defined element-wise; its proof can be found in Appendix E.
where C is a lower-triangular matrix of size p k ×( p k +1) with all entries equal to −1. All N per i , i = 0, . . . , n per −1, are, moreover, locally linearly independent and form a basis for S −,per .
Once again, we interpret 0-and 1-forms as linear functionals on the cell complex, and this leads to a discrete representation of the Stokes theorem. The corresponding graphical representation showing the action of the discrete exterior derivative on the 0-form degrees of freedom is presented in Fig. 3.

The tensor-product spline complex and mapped geometries
Using the univariate multi-degree spline complexes, we can build multivariate spline complexes via tensorproducts of the multi-degree spline spaces. Here, since we are mainly interested in surfaces in R 2 or R 3 , we only focus on bivariate spline complexes. First, we detail how these are built on a rectangular parametric domain Ω , and then we show how they can be used to build spline complexes on mapped surfaces in R 2 or R 3 .

Tensor-product splines
Let S i be some univariate multi-degree spline spaces built on Ω i and denote the corresponding sets of MDBsplines with {N i j : j = 0, . . . , n i − 1}, i = 1, 2. Unless otherwise specified, spline spaces S i are allowed to be non-periodic or periodic; we simply drop the superscript of "per" to simplify the notation whenever the context is unambiguous. Moreover, let {N i j : j = 0, . . . , n i − 1}, i = 1, 2, denote the respective sets of preferred 1-form basis functions corresponding to S 1 and S 2 ; these span the spline spaces S 1,− and S 2,− , respectively. Using the above univariate spline spaces, we define the following tensor-product bivariate spline spaces,

Definition of the complex
Using the above tensor-product spline spaces, we choose the spaces of 0-, 1-and 2-forms on Ω as follows, Then, the bivariate tensor-product spline complex on Ω is defined as Theorem 4.1. The cohomology spaces of the complex T satisfy: R , either S 1 or S 2 is periodic , R 2 , both S 1 and S 2 are periodic ; Proof. It is clear that only constants in Λ (0) T are annihilated by the exterior derivative, thus showing that H 0 (T) = R. The proof for H 1 (T) and H 2 (T) for non-periodic spline spaces can be obtained by, for instance, following the proof of [12,Theorem 4.1]. The cases with periodic spline spaces make Ω a topological cylinder or torus, and the cohomology spaces can be verified to be, , both S 1 and S 2 are periodic ; if only one of S 1 and S 2 is periodic , , both S 1 and S 2 are periodic . ■

Basis for discrete differential forms
We have already chosen the preferred basis for 0-, 1-and 2-forms in Eqs. (43)- (46). Since all basis functions are tensor-product, their element-wise computations are done by tensoring the respective element extraction matrices from Eqs. (31) and (38) for the splines N i j , i = 1, 2, and using Propositions 3.6 and 3.9 for the splines N i j , i = 1, 2. Therefore, it only remains to derive the discrete representations of the exterior derivatives akin to the univariate setting. We do so in the following. Let where f and N (0,0) are column vectors obtained by placing f i j and N 1 i N 2 j in the (i + jn 1 )-th locations, respectively. Then, using the univariate relations, we can write where N (1,0) and N (0,1) are column vectors obtained by placing N 1 i N 2 j and N 1 i N 2 j in the (i + jn 1 )-th and (i + jn 1 )-th locations, respectively, and the discrete exterior derivatives are given by Similarly where f 1 and f 2 are column vectors obtained by placing f 1 i j and f 2 i j in the (i + jn 1 )-th and (i + jn 1 )-th locations, respectively. Then, using the univariate relations, we can write where N (1,1) is a column vector obtained by placing N 1 i N 2 j in the (i + jn 1 )-th locations, and the discrete exterior derivatives can once again be derived to be the following sparse outer products,

Degree of freedom interpretation
The geometric interpretation of the tensor-product complex follows directly from those for the univariate multidegree complexes; see Sections 3.2.4 and 3.2. 7 . This time we consider a tensor-product partition of [0, 1] 2 into n 1 × n 2 quadrilaterals. The zero-dimensional, horizontal and vertical one-dimensional, and two-dimensional cells of this partition will be denoted, respectively, as Then, the degrees of freedom of the 0-, 1-and 2-forms are associated to these geometric objects. We define the oriented boundaries of the edges and the faces as Then, the discrete exterior derivatives from Eqs. (49) and (51) help us establish discrete versions of the Stokes theorem. The action of these on the spline degrees of freedom is presented in Fig. 4. Furthermore, it can be readily checked that, for f ∈ Λ (0) T , Eqs. (49) and (51) imply that d • d f = 0. Alternatively, this fact is implied by the duality of the discrete exterior derivatives with the boundary operator since as the boundary of a boundary is always empty.

Mapped geometries
Let us now transfer the spaces of spline differential forms onto a domainΩ ⊂ R d , d = 2 or 3, obtained via a geometric mapping of Ω . In particular, sticking to the isogeometric concept, we will look at geometric mappings built using tensor-product splines in S (0,0) and, moreover, assume thatΩ is a manifold. For Then, consider a 2-manifoldΩ obtained as the image of Ω under the spline map G defined as ThenΩ has local, curvilinear coordinates x 1 , x 2 , and global Cartesian coordinates y 1 , . . . , y d . Assuming thatΩ is a C ≥1 smooth manifold, the vectors ∂ x i = ∂ G ∂ x i (x), i = 1, 2, form a basis for vectors tangent toΩ at G(x). The vectors dual to ∂ x i are d x i , i = 1, 2, and they form a basis for the covectors tangent toΩ at G(x). The 0-, 1-and 2-forms onΩ can thus be denoted as f , f i d x i and f d x 1 ∧ d x 2 , respectively. Denote the associated metric tensor, its i j-th component and its matrix determinant with g, g i j and g, respectively, and let ∂ y i , i = 1, . . . , d, be the canonical basis vectors in R d . Then, The quantity √ g thus denotes the Jacobian determinant of the map G.
If the mapΩ is locally invertible, then √ g > 0 and the components of the inverse of the 2 × 2 matrix [g i j ] are denoted as g i j . Using the metric and its inverse, we can explicitly define the Hodge star ⋆ in the present setting, where ϵ i j is equal to 1 for (i, j) = (1, 2), equal to −1 for (i, j) = (2, 1), and zero otherwise. In particular, the L 2 inner product of i-forms can be expressed as (also see Appendix F) 4. The univariate degree of freedom interpretations for the multi-degree complexes directly lead to the same for the tensor-product complex T. The 0-, 1-and 2-forms are now associated to the vertices, edges and faces of a tensor-product cell complex, respectively; see the figures in the first row. These figures correspond to n 1 = n 2 = 5. Moreover, with respect to the preferred basis, the discrete exterior derivatives have D (1,0) , D (0,1) , D (2,0) and D (0,2) have a simple action as shown in the middle and bottom rows.
Finally, differential forms onΩ in the canonical basis dy i , i = 1, . . . , d, can be pulled back to Ω using the map G * as follows, The map G * is called the pullback and it commutes with both the wedge product and the exterior derivative. Moreover, using it, we can perform integration of an i-form on an i-dimensional geometry G ( Ω ) as Fig. 5. A single edge or a pair of opposite edges of a tensor-product spline patch can be collapsed for creating geometries with polar singularities. We will refer to the two modes of collapse as Type 1 and Type 2, respectively. The collapsed edges here are shown in red, and the black edges are identified with each other to enforce periodicity.
Using the pullback, we also define the spaces of spline differential forms onΩ aŝ and the corresponding spline complex onΩ is defined aŝ T .
The pullback commutes with the exterior derivative and, thus, forms a cochain map from the complexT to T.

The polar spline complex
In this section, we build a spline complex on geometriesΩ that are obtained via a map G that collapses one or two edges of Ω to one or two points, respectively, in R d ; see Fig. 5. These collapsed edges are called polar singularities or poles. Fig. 6 presents a topological representations of the tensor-product degree-of-freedom complexes following the introduction of polar singularities.
In general, the presence of poles means thatΩ will not be a C ≥1 smooth 2-manifold. However, by restricting each component of G to be a member of a suitable subspace of S (0,0) , we will be able to ensure smoothness of Ω . In Section 5.1, we build such a suitable subspace and use it to define smoothΩ ; the splines in the former will be called polar splines. Thereafter, in Section 5.2, we build spaces of polar spline differential forms onΩ and use them to define the polar spline complex in Section 5. 3.

Assumption 4: Tensor-product configuration for building polar splines
With Ω = Ω 1 × Ω 2 , the endpoints of Ω 1 have been identified. The univariate spline spaces S 1 and S 2 on Ω 1 and Ω 2 are periodic and non-periodic, respectively, and they are used to define all tensor-product spline spaces on Ω using Eqs. (43)- (46). Moreover, • with r i the smoothness vector used to define S i , i = 1, 2, • the dimension S 2 is at least 5, i.e., n 2 ≥ 5.

C 1 smooth polar B-splines
Thanks to the end-point interpolation and partition of unity properties of MDB-splines (see Proposition 3.4(d)), the required edge collapse shown in Fig. 5 can be achieved by choosing in Eq. (53) G 00 = G 10 = · · · = G (n 1 −1)0 ⇐⇒ ∀x 1 ∈ Ω 1 , G(x 1 , a 2 ) = G 00 ;  60) and (61). The middle row shows the tensor-product complex for (n 1 , n 2 ) = (4, 5). The rightmost vertex and vertical edge degrees of freedom have been plotted in grey to indicate that periodicity in the first parametric direction has been imposed; see Assumption 4. The bottom row shows the tensor-product complex following Type 1 collapse, while the top row shows the tensor-product complex following Type 2 collapse.
However, in general, this coefficient coalescing will introduce kinks at the poles and the surface representation will not be smooth. Nevertheless, it is possible to identify constraints on the remaining G i j that ensure thatΩ is a C 1 smooth 2-manifold or, equivalently, such that it has a well-defined tangent plane at all points. Such constraints were identified in [29] for C k smoothness, but for simplicity we restrict to the case of C 1 smoothness. In this section, we present the relevant constraints and their resolution. The discussion will be abbreviated and focused on practical considerations since the theory has already been elaborately addressed in [29] and, more recently, [32].
A polar surface will be smooth at a polar point if it can be locally (re)parametrized in a smooth way. Such parametrizations can be specified in a constructive manner and, for C 1 smoothness, they impose simple geometric constraints on the choice of the G i j [29,Section 3.3]; these are presented in the following result. (a) For the edge-collapse in Eq. (60),Ω has a well-defined tangent plane at G 00 if (i) the points G i j , i = 0, . . . , n 1 − 1 , j = 0, 1, are all coplanar; (ii) the vectors G i1 − G i0 , i = 0, . . . , n 1 − 1, are all distinct, non-zero, and form a clockwise or counter-clockwise fan around G 00 .
In particular, assigning Ω a counter-clockwise orientation, G preserves the orientation in a neighbourhood of the poles if the fans in (a) and (b) above are clockwise and counter-clockwise, respectively.

Proof. See [Section 3.3][29].
■ Depending on Type 1 or Type 2 edge collapse (see Fig. 5), we would want the satisfaction of either the conditions in Proposition 5.1(a), or those in both Proposition 5.1(a) and (b), respectively. Then, [29] suggest the following simple way of satisfying the above smoothness constraints at the poles. Choose triangles △ 1 and △ 2 with vertices , respectively. Next, require the following relations to hold, In other words, with regard to Proposition 5.1(a) (respectively, Proposition 5.1(b)), we force G i j to lie in the plane of △ 1 (respectively, △ 1 ), and the numbers χ k,1 i j (respectively, χ k,2 i j ), k = 1, 2, 3, are its corresponding barycentric coordinates. We will call △ 1 and △ 2 the domain triangles for the above sets of G i j .
Next, conditions (ii) of both Proposition 5.1(a) and (b) can be satisfied equally easily by choosing the barycentric coordinates as follows. Define θ i ∈ (0, 2π ) as Next, compute the required barycentric coordinates using the following two equations for i = 0, . . . , n 1 − 1, Lemma 5.2. All χ k,1 i j and χ k,2 i j , as specified by Eqs. Proof. As explained in [29], Eqs. (65) and (66) implicitly impose that the G i j are uniformly distributed on circles centred at the poles and, moreover, that these circles are contained within the domain triangles. Therefore, the barycentric coordinates are all non-negative.
■ Therefore, depending on Type K collapse, K = 1, 2, the number of coefficients needed to define the tensor-product coefficients is equal to n pol , Doing so, and recalling Eq. (43), let us define the polar extraction E pol of size n pol × n (0,0) as where E 1 and E 2 are defined as ( j, k) = (0, 1) and Type 1 or Type 2 or ( j, k) ∈ {(1, 1), (n 2 − 2, 2)} and Type 2 , where N pol is a vector containing the polar B-splines N Then, the following result holds.
where E pol kℓ is the kℓ-th entry of E pol . Then, it is clear that G i j satisfy the constraints of Proposition 5.1. The smoothness of the pushforwards of polar B-spline functions was shown in [29] and we omit the proof here for brevity. ■

Basis for discrete differential forms
Let us now describe the construction of polar spline discrete differential forms that are built with S pol chosen as the space of 0-forms. First, let us describe the motivation behind, and an overview of, our construction.

Motivation for the construction
The motivation for the construction presented herein is derived from the relations for mapped geometries presented in Section 4. 4. In particular, the introduction of edge-collapses implies that √ g = 0 at the poles. This implies that, in general, spline differential forms will not be bounded in a neighbourhood of the poles; e.g., see Eq. (55). Equipped with the C 1 smooth polar B-splines as 0-forms, we counteract this singular behaviour by imposing "local exactness" for all 1-and 2-forms in a neighbourhood of the poles. That is, in the vicinity of the poles, spline differential k-forms, k = 1, 2, will be restricted to be exact. Then, at the poles, C 1 smoothness of the 0-forms automatically translates to C 0 and C −1 smoothness of the 1-and 2-forms and, moreover, avoids the blowup. Note that, away from the poles, all differential forms are going to have the same smoothness as their tensor-product counterparts.
Following the above motivation, and by the construction of the 0-form polar splines as in Eq. (70), we now present polar analogues of the tensor-product cell complexes from Fig. 4. These are shown in Fig. 7; c.f. Fig. 6. Note the following about the polar degree-of-freedom complex in the top and bottom rows.
• There are three degrees of freedom for polar 0-forms near the poles.
• Imposing local exactness of polar 1-forms at the poles, there are two degrees of freedom for them near the poles.
• Imposing local exactness of polar 2-forms at the poles, and from the above bullet, there are no degrees of freedom for them near the poles.
Next, let us explain the different vertical and horizontal maps in Fig. 7. The left-most column of the figure corresponds to Eq. (73), i.e., the vertical maps in that column send degrees of freedom for polar 0-forms to those for tensor-product 0-forms. The horizontal maps in the middle row have been defined in Eqs. (49) and (51). It remains to define the remaining vertical (i.e., from polar degrees of freedom to tensor-product degrees of freedom) and horizontal maps (i.e., discrete exterior derivatives that act on polar degrees of freedom). The transposes of the vertical maps will help specify the basis functions for polar 1-and 2-forms as linear combinations of those for tensor-product 1-and 2-forms, respectively, similarly to Eq. (70) for polar 0-forms.
A concrete discussion in the following subsections requires a numbering of the degrees of freedom for 0-, 1and 2-forms in the Type 1 and Type 2 polar complexes from Fig. 7; these numberings are then shared by the basis functions for polar 0-, 1-and 2-forms, respectively. The total number of degrees of freedom associated to 0-, 1-and 2-forms can be found using Eqs. (43)- (46). Observing that Assumption 4 implies n 1 = n 1 , the number of degrees of freedom are computed as below for Type K collapse, K = 1, 2, n (1),pol := n (1,0) + n (0,1) n (2),pol := n (1,1) − K n 1 .
The degrees of freedom are numbered using the scheme shown in Fig. 8.  (2),pol such that both the top-two and bottom-two rows form commutative diagrams. The latter polar extraction operators specify the basis for polar 1-and 2-forms as suitable linear combinations of the tensor-product 1-and 2-form basis functions. (Note that, to simplify the notation, we use the same symbols to denote the extractions and discrete exterior derivatives for Type 1 and Type 2 polar complexes.).

Polar 0-forms
Define S (0),pol := S pol and E (0),pol := E pol . Then, using Eq. (70) and for any f ∈ S (0),pol , where f is a vector containing all coefficients f i . − 1 for 1-forms, and from n (2),pol − 1 for 2-forms. Also, keeping Fig. 7 in mind, the numbering near the top pole (which is applicable only for Type 2 collapse) assumes that we are looking down at the pole from above the degree-of-freedom complex.

Polar 1-forms
Let f ∈ S (0),pol and apply the exterior derivative to it. Then, using Eq. (48), Then, with reference to Fig. 9, let us define two maps D (0),pol and E (1),pol . We do so by defining their actions on the degrees of freedom, considering both regions away from the poles and near the poles.
• Fig. 9(a): Far away from the poles, the degree-of-freedom complexes locally look like their tensor-product counterparts, i.e., the topology is that of a structured quadrilateral grid; c.f. Figs ,pol g 1 kℓ := g j 1 , g 2 kℓ := g j 3 , g 1 k(ℓ+1) := g j 2 , g 2 (k+1)ℓ := g j 4 . The relations between the different subscripts is easily deciphered from Eq. (73) and the degree-of-freedom numbering shown in Fig. 8. More precisely, fixing k, ℓ, Fig. 8 implies the following relations between them, In the above, k, ℓ can be any values from the following ranges for Type K , K = 1, 2, collapse, • Fig. 9(b): Consider the degrees of freedom f i 1 , . . . , f i 5 for a polar 0-form. Then, the relations between the indices depend on whether the pole is the bottom one or the top one.
Bottom pole: Assuming that the pole is the bottom one, let us define the actions of D (0),pol and the transpose of E (1),pol . Given k, ℓ, the following relations between the different indices hold from Fig. 8(a), Here, ℓ = 0 and k ∈ {0, . . . , n 1 }. Recall Eq. (73) that helps map degrees of freedom for a polar 0-form onto those for a tensor-product 0-form. Using that map, we define the actions of D (0),pol and the transpose of E (1),pol as below, Top pole: For the case of the top pole, the above relations can be analogously defined by suitable indexrelabelling; we present them here for completeness. Firstly, the following relations between the indices hold from Fig. 8(b), Here, ℓ = n 2 − 3 and k ∈ {0, . . . , n 1 }. It is important to note here that, compared to Fig. 8(b), the edges associated to the degrees of freedom g j 3 and g j 4 are oppositely oriented. Since Fig. 8 provides a simple and global way of assigning orientations, with regards to the global operator definitions we are actually interested in the degrees of freedom g j 3 and g j 4 , g j 3 := −g j 3 , g j 4 := −g j 4 . (87) Then, using Eq. (73), we define the actions of D (0),pol and the transpose of E (1),pol as below, Proposition 5. 6. Eqs. (78), (84) and (88) imply that the diagrams in Fig. 9(a) and (b) commute.
Proof. The linear-independence claim follows from the full rank of the extraction operator E (1),pol .
• Similarly, Eqs. (84) and (88) imply that, near the bottom and top poles, respectively, the non-zero parts of the first two rows of E (1),pol are obtained by taking differences of the columns of the matrices from Eq. (69). The latter matrices have rank 3, and implies that the first and last two rows of E (1),pol are also linearly independent.
The smoothness of the one-form polar B-splines is implied by their local exactness at the poles (see Section 5.2.1) and Proposition 5. 5. ■

Polar 2-forms
Let f ∈ S (1),pol and apply the exterior derivative to it. Then, using Eq. (50), However, by the local exactness of f at the poles (see Section 5.2.1), the above implies that g i j = 0 if Type 1 collapse : j = 0 , Type 2 collapse : j = 0, n 2 − 1 .
Then, with reference to Fig. 10, let us define two maps D (1),pol and E (2),pol . As in the previous section, we do so by defining their actions on the degrees of freedom, considering both regions away from the poles and near the poles.

Definition of the complex
Using the above polar spline spaces, we choose the spaces of 0-, 1-and 2-forms on Ω as follows, .

Proof.
The fact that S is a cochain complex is immediate from the construction of the polar spline spaces for 0-, 1-and 2-forms. We prove the claims for each cohomology space separately.
As a consequence, d f = 0 implies that all degrees of freedom f i are equal to some α ∈ R. Then, by the partition of unity property of the 0-form polar B-splines, First cohomology. For the cochain complex S, we have the following equivalence between alternating sums of the dimensions of the cohomologies and the dimensions of vector spaces that form the complex, The right hand-side follows from Eq. (74)-(76) and is equal to K for type K collapse. Then, since the first term on the left is equal to 1, to prove that H 1 (S) is trivial, we only need to show that the last term on the left is equal to K − 1; this is proved in the following. Second cohomology. Note that d f = 0 for any f ∈ Λ (2) S . Then, let us build an h ∈ Λ (1) S such that dh = f . This will always be possible for Type 1 collapse, while for Type 2 collapse we will need to place one constraint on f . The claim will thus follow. We define such an h by defining its degrees of freedom, and we start at the bottom pole; see Fig. 8(a). In the following, unless specified otherwise, the index i runs from 0 to n 1 − 1.
This completes the definition of h and it can be verified that dh = f . For Type 2 collapse, such an h can be found only when the above constraint is satisfied, implying that the cohomology space H 2 (S) is one dimensional. is not in the image of d for Type 2 collapse. This polar 2-form is thus a representative element of the one dimensional cohomology space H 2 (S). In particular, this cohomology is isomorphic to the following vector space; This vector space is closely related to the idea of "discrete harmonic forms" and it will play an important part in the numerical tests; see Section 6. 3. Finally, the polar spline complex onΩ is defined as in Section 4. 4. That is, we defineΛ (i) S as below, and the corresponding spline complex onΩ is built using them, with the pullback again acting as a cochain map fromŜ to S.

Numerical tests
This section numerically investigates the approximation power and stability of the polar spline complexes by solving problems on smooth polar geometries in R d , d = 2, 3. In particular, we consider approximation of the Stokes flow on both fixed and deforming closed surfaces.
Our approach towards spline differential forms is well-suited for computations within the classical framework of finite element assembly loops. Indeed, starting from element-local representations of univariate splines (Section 3.2), tensor-product splines can be readily built. Subsequently, the tensor-product splines can themselves be combined using the polar extractions to build polar splines on each element of the two-dimensional parametric domain Ω . This approach is adopted for all computations presented here.

Spline spaces and geometries
For brevity, we only present numerical tests with the polar spline spaces as they already utilize the univariate and tensor-product splines defined in Sections 4 and 3.2 . Moreover, we only consider Type 2 collapse, i.e., closed polar manifolds. The numerical tests presented here show that the polar spline spaces demonstrate optimal approximation; similar results were obtained for the configurations not shown here (e.g., univariate and tensor-product spline spaces, Type 1 collapse).
Specifically, we consider three polar spline spaces built using Type 2 collapse. Each is built using tensor-product spline spaces of uniformly chosen bi-degree ( p, p), p = 2, 3, 4. That is, the univariate spline spaces used to build the tensor-product spline spaces are defined by choosing the polynomial degree on each element equal to p. The breakpoints and associated orders of smoothness are defined as below. : (x 0 , x 1 ) = (0, 1) , (r 0 , r 1 ) = (−1, −1) .
The polar manifolds built using the above spline spaces are shown in Fig. 11. The black lines delineate the Bézier elements of the mesh.

L 2 projection
As the most basic test of the approximation power of the individual spline spaces, we solve L 2 projection problems. Specifically, given f ex ∈ L 2 Λ (i) (Ω ), we find f ∈Λ (i) S such that The exact solutions are chosen to be where h(y 1 , y 2 , y 3 ) = sin ( 2π The L 2 -projection problems project the pullbacks of the above exact solutions onto the appropriate polar spline spaces; c.f. Eq. (57). With the approximation error defined as e := f ex − f , the L 2 norm of e and de displayed in Fig. 12. The error norms are plotted against the square root of the number of degrees of freedom. Note that the norm of de is omitted for 2-forms since the exterior derivative maps all 2-forms to zero. To the left of each error convergence plot, we also show the exact solutions for each polar geometry (see the online version of this article for high resolution pictures): • for 0-forms, the surface is coloured by f ex and the tangential vector field represents d f ex ; • for 1-forms, the surface is coloured by values of d f ex and the tangential vector field represents f ex ; • for 2-forms, the surface is coloured by values of f ex .
For optimal approximation, we expect the L 2 norm of e to decrease with order p+1 for 0-forms and p otherwise. The L 2 norm of de is expected to decrease with order p for both 0-and 1-forms. As can be observed in Fig. 12, the polar spline spaces demonstrate optimal approximation behaviour for all p.

Generalized Stokes flow
We now consider generalized Stokes flow on fixed and deforming polar manifoldsΩ . This problem is important when studying, for instance, fluid flow on biological membranes such as lipid bilayers, and the problem formulation can be derived from first principles; e.g., see [34,46] for the derivation of Stokes flow 1 . We express the strong form Fig. 11. The polar manifolds used in the numerical tests are shown above. The black lines delineate the Bézier elements at the coarsest refinement level. As can be seen, the coarsest meshes consist of 9, 8 and 5 elements for p = 2, 3 and 4, respectively. See Section 6.1 for details about the underlying polar spline spaces.  Fig. 11 are shown above; see Section 6.2. of the generalized Stokes problem onΩ as The first equation pertains to momentum conservation onΩ while the second one is the equation of mass conservation. Here, q is the pressure (2-form), u is the velocity (1-form) and ν represents the (instantaneous) normal velocity field of the deforming domainΩ (2-form). Moreover, f is an external force (1-form) on the system, h is a source of mass production (2-form), µ is the viscosity and α is a scalar constant. The remaining terms are related toΩ -the metric tensor, g; the second fundamental form, k = k i j d x i ⊗ d x j ; twice the mean curvature, H ; and the Gaussian curvature, κ. In particular, with n denoting the unit normal vector toΩ , Remark 6.1. The form of the generalized Stokes problem in Eq. (113) can be related to the usual vector calculus notation; see [34,46], for instance. In particular, we would like to mention that the velocity 1-form is related to the tangential fluid velocity onΩ , denoted u # = u #,i ∂ x i , as below, From Appendix F, u # is interpreted as a proxy field of −⋆u. In particular, du is proportional to the surface divergence of u # , i.e., Remark 6.2. Interestingly, even in the absence of mass production (h = 0) and external forcing ( f = 0), tangential flow onΩ can be produced if the surface has non-zero normal velocity, ν. This is directly related to the interpretation of the deforming surfaceΩ as an inextensible material. An interesting example of such materials is lipid bilayers [34], which are envelopes for eukaryotic cell contents. These behave as in-plane fluids and out-ofplane solids. In particular, as in Eq. (113), the out-of-plane velocities of these surfaces, governed by solid mechanics, lead to in-plane flow. See Section 6. 3.3 for examples of this phenomena.

Manufactured solution
In order to numerically verify optimal approximation of generalized Stokes flow, we create a smooth manufactured solution to the problem. SinceΩ is parametrically defined using piecewise polynomials, it does not have a simple implicit representation unlike other surfaces (e.g., spheres). This makes the derivation of a smooth manufactured solution a complicated task. For instance, fixing u and h, mass conservation implies ν = (du − h)/H . The derivatives of ν, needed for momentum conservation, therefore involve derivatives of the mean curvature H ; note that these derivatives will clearly have a lower regularity than that of the surfaceΩ . Keeping these difficulties in mind, a smooth manufactured solution can nevertheless be derived by either • assuming thatΩ is a Type 1, flat polar geometry so that, in particular, both k and H are trivial; • or, by assuming that ν = 0.
We adopt the second approach above so that we can demonstrate optimal approximation on arbitrarily curved polar geometries.
Therefore, choosing ν = 0, the exact solutions for the different variables are chosen to be Using the above u ex and p ex , we define h ex = du ex and f ex = d ⋆ q ex + µ (2κu ex − dd ⋆ u ex ) − αu ex . For the above choice of manufactured solution, we numerically solve the generalized Stokes problem in mixed form by introducing w = d ⋆ u, the vorticity (0-form). The corresponding weak form of the discrete problem is defined as follows.
(118) Fig. 13. Generalized Stokes flow w/ manufactured solution: The convergence rates for errors in vorticity, w, velocity, u, and pressure, q, are shown above; see Section 6. 3.1.
For simplicity, we set µ = α = 1. For this mixed problem, we also numerically compute the inf-sup constant γ S defined as This constant can be numerical computed by solving a generalized eigenvalue problem [47]. In view of Theorem 5.10 and the fact that we are looking at Type 2 polar geometries, we expect to have n (2),pol − 1 non-zero eigenvalues as the second cohomology space is one dimensional. The constant γ S is the square-root of the smallest non-zero eigenvalue. Fig. 13 show the results of the numerical approximation. The following information has been presented.
• The exact solutions for w ex , dw ex , u ex , du ex , q ex and d ⋆ q ex have been plotted on the polar geometries. The values of 0-and 2-forms are used to colour the surfaces, and the 1-forms are displayed as tangential vector fields.
• With e □ := □ − □ ex , □ ∈ {w, u, q}, the L 2 norms of e w , de w , e u , de u and e q have been plotted.
• Finally, the value of the inf-sup constant γ S at each refinement level has been labelled in the plot where e q has been shown.
For optimal approximation, and for polar splines built using tensor-product splines of bi-degree ( p, p), we expect the errors for all 0-forms to reduce with order p + 1 and for all 1-and 2-forms with order p. As shown in Fig. 13, all polar spline spaces demonstrate optimal approximation behaviour. Moreover, it can be seen that the spline spaces are inf-sup stable, i.e., the constant γ S does not deteriorate with mesh refinements. The left plot in each row displays the vorticity, w; the middle plot shows streamfunction contours, ψ, and the tangential vector field, u # ; and the right plot displays the pressure, q. The surface of the middle plot has been coloured by the value of du; this value is of the order of machine precision and thus the surface is uniformly coloured grey.

Pointwise incompressibility
SinceŜ is a cochain complex, we know that d mapsΛ (1) S intoΛ (2) S . (From Theorem 5.10, this map is a surjection only for Type 1 collapse as illustrated by the vanishing cohomology.) Then, if h ex ∈Λ (2) S in Eq. (118), then du is going to be pointwise equal to h ex . In particular, if h ex = 0, then the discrete velocity 1-form is going to be closed, i.e., du = 0 pointwise. Equivalently, the discrete tangential velocity u # is going to have pointwise zero surface divergence; see Remark 6.1. Moreover, from Theorem 5.10, du = 0 implies that there exists a streamfunction ψ (0-form)such that dψ = u.
We illustrate the above simple fact by solving the problem in Eq. (118) on the surfaces in Fig. 11 for the third refined level. It is important to note that pointwise incompressible solutions can be obtained for any refinement level, no matter how coarse or fine -the choice of the third refinement level is only to ensure accuracy of the discrete solutions. We choose h ex = 0 = ν ex and the forcing f ex is chosen to be equal to the one in Eq. (118). The results are shown in Fig. 14. The left figures in each row correspond to the computed w; the middle figures to the tangential velocity vector field, u # , and contours of the streamfunction, ψ; and the right figures correspond to the pressures, q.  15. Deforming domains: The above numerical tests show the discrete solutions to the generalized Stokes problem whenΩ is deforming with a prescribed normal velocity field; see Eq. (121). The left plot in each row displays the vorticity, w; the middle plot shows the tangential vector field, u # , in blue and the imposed normal velocity, (⋆ν ex )n, in red; and the right plot displays the pressure, q. The relative scaling of the normal velocity arrows with respect to the tangential velocity arrows is 1.4 for p = 2 and 2.4 for p = 3 and 4.

Deforming domains
As a final numerical example, we consider the case where ν ex ̸ = 0. We choose f ex = 0 = h ex and ν ex = √ g f 1 d x 1 ∧ d x 2 , where We also impose q ex ⊥ h (2) . The weak problem thus becomes to find (w, u, q, v) ∈Λ (0) S ×Λ (1) S ×Λ (2) S × h (2) such that for all (z 0 , z 1 , z 2 , z 3 ) ∈Λ (0) S ×Λ (1) S ×Λ (2) S × h (2) , Again, for simplicity we set µ = α = 1. Fig. 15 shows the results for the polar geometries in Fig. 11. The left figures in each row correspond to the computed w; the middle figures to the tangential fluid velocity, u # , and the normal velocity, (⋆ν ex )n = f 1 n; and the right figures correspond to the pressures, q.

Conclusions
We have investigated the development and applications of high-order accurate, spline differential forms in a variety of settings. In the univariate setting, we provide the construction of multi-degree spline differential forms, i.e., smooth, piecewise-polynomial differential forms that allow for local degree elevation; Section 3.2. The construction is presented within the paradigm of Bézier extractions, thus making our algorithms and approach easily implementable, particularly within element loop-based finite element software. In the bivariate setting, we first build spline differential forms using tensor-products of the univariate spline differential forms; Section 4. The properties of the univariate splines carry over to the tensor-product splines; this approach is easily extensible to higher dimensions. Next, in the bivariate setting, we focus on the case of singularly parametrized smooth surfaces (with and without boundary) called polar surfaces, and build spline differential forms on them; Section 5. Finally, the spline differential forms are used to solve L 2 projection problems and generalized Stokes flow on smooth polar surfaces in R 3 ; Section 6. The results demonstrate optimal approximation and inf-sup stability of the spline spaces for the Stokes problem, simulations of pointwise incompressible flows, and simulations of flows on deforming inextensible surfaces.
Our approach here has been constructive. There are several extensions of the results and applications presented herein that are related to topics in computational mechanics as well as several areas of mathematics, such as algebraic topology, differential geometry and numerical analysis. As such, there are several opportunities for future research into theoretical and practical aspects of spline-based exterior calculus with the goal of solving PDEs on surfaces; a few of these are itemized below.
• For efficient computations, development of smooth polar differential forms using adaptively-refined splines, such as hierarchical B-splines [28] is important.
• Similarly, the 0-, 1-and 2-form spline spaces developed here are H 2 , H 1 and L 2 conforming on polar surfaces. Extending this construction to higher orders of regularity [29] is interesting, for instance, for the variational multiscale framework for divergence-free flow simulations [48] where the 1-and 2-form spaces need to be H 2 and H 1 conforming, respectively.
• Another important research question is the development of commuting projection operators that help theoretically verify the stability of the spline spaces. Bounded cochain projections from the de Rham complex to spline complexes are needed for provable well-posedness at the discrete level. The theory on such cochain projections for adaptively refined and non-tensor-product splines is missing, even in the absence of singular parametrizations, and is beyond the scope of this paper. Instead, the examples in Section 6 provide numerical evidence of the well-posedness of the discrete problems.
• On the side of applications, we simulate flows on smoothly deforming surfaces with prescribed normal velocities. Incorporating mechanics into the equations to simulate the behaviour of fluid membranes, e.g., lipid bilayers, is a particularly interesting extension.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. In the above it is assumed that any undefined N k,q−1 with k < p − q + 2 or k > n is equal to the zero function, and that if M k,q−1 = 0 then

Appendix B. Algorithmic definition of MDB-splines
For 1 ≤ k ≤ m − 1 let K k,− be a matrix of size ( p k + 1) × (r k + 1), whose j-th column, j = 0, . . . , r k , is given by, and let K k,+ be a matrix of size ( p k+1 + 1) × (r k + 1), whose j-th column, j = 0, . . . , r k , is given by, Using these matrices, we can build the matrix K k of size θ (m) × (r k + 1) which contains all constraints required to enforce C r k smoothness at x k . This matrix is defined row-wise in the following manner: • j 1 = k 2 ≤ j 2 = k 1 and A 1 is a left-inverse of A 2 , A 1 A 2 = I j 1 ; • k 3 = k 1 , j 3 = k 4 and A 3 maps the nullspace of A T 2 to that of A 4 , Proof. Let A T 1 A T 2 = I k 1 + A 5 . Then, each column of A 5 must be in the nullspace of A T 2 since Then, using the nullspace-preserving property of A 3 , the claim follows,

Appendix D. Algorithmic definition of periodic MDB-splines
The periodic version of MDB-splines is built by starting from the non-periodic version. First, referring to Eqs. (122) and (123), build the matrices K m,− and K 0,+ . These are matrices of sizes ( p m + 1) × (r 0 + 1) and ( p 0 + 1) × (r 0 + 1), respectively. Using these matrices, we build the matrix K m of size θ (m) × (r 0 + 1) which contains all constraints required to enforce C r 0 smoothness at x 0 ≡ x m , With H the multi-degree extraction corresponding to the non-periodic spline space, let P be any row permutation matrix such that

such that
• 0 i are zero matrices with the number of rows equal to j i ≥ 0, i = 1, 2 (i.e., j i = 0 implies an empty matrix with no rows); •K has 2(r 0 + 1) rows.
We can always find such a P thanks to end-point derivative property of MDB-splines [30, Proposition 2.11(c)] and Eq. (25). Then, Algorithm 3 helps compute a matrixH of size n per whose rows are in the left-nullspace of P H K m . In addition to the above definitions, the algorithm uses the row permutation matrix Q defined as The highlighted entry of Q lies in its j 1 -th column (with the first column having the index 0). This permutation matrix ensures that the p k + 1 splines supported on Ω k have indices µ(k), . . . , µ(k) + p k ; recall the interpretation of indices in this periodic setting from Eq. (36).
Proposition D.1. With Assumptions 2 and 3 in place, the matrixH has full rank.