Solving Biharmonic Equations with Tri-Cubic C 1 Splines on Unstructured Hex Meshes

: Unstructured hex meshes are partitions of three spaces into boxes that can include irregular edges, where n (cid:54) = 4 boxes meet along an edge, and irregular points, where the box arrangement is not consistent with a tensor-product grid. A new class of tri-cubic C 1 splines is evaluated as a tool for solving elliptic higher-order partial differential equations over unstructured hex meshes. Convergence rates for four levels of reﬁnement are computed for an implementation of the isogeometric Galerkin approach applied to Poisson’s equation and the biharmonic equation. The ratios of error are contrasted and superior to an implementation of Catmull-Clark solids. For the trivariate Poisson problem on irregularly partitioned domains, the reduction by 2 4 in the L 2 norm is consistent with the optimal convergence on a regular grid, whereas the convergence rate for Catmull-Clark solids is measured as O ( h 3 ). The tri-cubic splines in the isogeometric framework correctly solve the trivariate biharmonic equation, but the convergence rate in the irregular case is lower than O( h 4 ). An optimal reduction of 2 4 is observed when the functions on the C 1 geometry are relaxed to be C 0 .


Introduction
The efficient representation of volumetric C 1 fields over hexahedral meshes is of interest in areas ranging from scientific data visualization to solving higher-order differential equations. For example, to visualize a flow computed by the Discontinuous Galerkin approach currently requires substantial post-processing to extract streamlines that the theory predicts to be smooth [1]. Engineering analysis based on splines is efficient in that it a priori bakes in the smoothness required of the solution of higher-order partial differential equations. Splines can serve both to define the geometry of the physical space and supply the degrees of freedom for numerical analysis on the manifold. The case of volumetric physical space, in three variables, is of obvious high practical interest. Where symmetries can not reduce the dimension, splines need to be well-defined over unstructured partitions of the physical domain. In particular, at irregularities, where the volumetric partition into boxes does not form a topological grid layout suitable for box splines [2], the construction of smooth solutions is challenging. Unstructured partitions of the physical domain into boxes can include irregular edges, where n = 4 boxes meet along an edge, and irregular points, where n = 8 boxes meet. The challenge is that, at the irregularities, there is no consistent extension of the individual pieces' parametric derivatives to the whole neighborhood unless their cross-product vanishes [3] (Lemma 3.7).
Until recently, the literature did not offer conforming polynomial C 1 splines over irregular box complexes. One approach with two variables introduces a removable singularity into the parameterization at irregular points, see [4][5][6][7]. This approach collapses the 1-jet of the first derivatives at the irregularity and forces the 2-jet of the second derivatives of Numerical experiments with four refinement steps, i.e., with up to half a million degrees of freedom, indicate • a 2 4 (fourth-order) convergence rate for Poisson's equation on irregular box-complexes, i.e., the error between the computed and the known exact solution in the L 2 norm decreases by a factor of 2 4 under halving of the mesh interval h, by 2 3 in the H 1 error, and 2 2 in the H 2 norm. For the biharmonic (i.e., fourth-order, bi-Laplacian) equation, the observed convergence rate is also • 2 4 for the regular case and for C 0 elements on C 1 geometry.
However, although still converging to the correct solution, the • convergence rate of singular C 1 tri-3 splines on singular C 1 tri-3 spline geometry is less than 2 4 on irregular box-complexes.
We note that the available convergence estimates, e.g., of [18,19], assume higher smoothness of the space, typically as high as for splines on the regular box-complex. That is, these estimates do not apply to solving the biharmonic equation with tri-3 C 1 splines. A likely explanation is that the singular, C 1 -constrained tri-variate spline space does not have full approximation power. On the other hand, convergence is consistently better than that of Catmull-Clark solids [20,21]: • in all cases enumerated above, tri-3 C 1 splines exhibit faster convergence than Catmull-Clark solids.
Overview. After a brief literature review of trivariate smooth elements and the bivariate antecedents of tri-3 C 1 splines, Section 2 defines the tri-3 C 1 splines for unstructured boxcomplexes. The space is C 1 on regular local grids but is C 2 if initialized by knot insertion on a locally tensor-product grid. The space has zero first derivatives across irregularities but is C 1 after a change in variables and has 2 3 linear-independent B-spline-like basis functions per box. Section 4 shows the numerical convergence for Poisson's equation and compares the convergence to that of Catmull-Clark solids. Section 5 shows and discusses the convergence for the biharmonic equation and compares it to Catmull-Clark solids.

Smooth Trivariate Finite Elements
The grid points of a regular partition of 3-space into boxes can be interpreted as the control points of a tri-variate tensor-product spline with one polynomial piece per cube. The theory of such splines is well-understood, see, e.g., [22,23]. However, the complex outer shape and internal partition lead to unstructured hex-meshes, see [24][25][26][27][28][29][30]. Our improved understanding of fields via their singularity graph [31,32] has not been matched by corresponding progress to more flexible spline spaces. For box-complexes where the tensor grid gives way to an irregular arrangement of boxes, including irregular points and irregular edges, there are multiple options, none of them perfect.
Geometric continuity in three variables is, in principle, well understood as a change in variables between pieces, see [33,34]. In practice, the trivariate geometric continuity has been barely explored: [35,36] join just one pair of trilinearly parameterized face-adjacent boxes and [37] consider trilinearly parameterized multipatch volumes with exactly one inner (irregular) edge. The challenge is the complicated interaction of reparameterizations surrounding an irregular point. This complexity is particularly pronounced when the polynomial (tensor-) degree is low, below tri-5, which is important in three variables to obtain manageable spline spaces (tri-4 polynomial pieces already have 125 coefficients). However, geometric continuity requires increased polynomial degrees near irregularities and careful book keeping to adjust reparameterizations under refinement. Generalized subdivision [3,38,39] creates an infinite sequence of nested piecewise polynomial layers that complicate the analysis, e.g., computing integrals near irregularities. Trivariate subdivision rules analogous to Catmull-Clark subdivision [38] have been proposed in [20] but come without a guarantee of smoothness and approximation order. Ref. [21] pioneered the use of Catmull-Clark solids in engineering applications. More recent work can be found in [40]. Ref. [41] solves the heat equation using interpolatory Catmull-Clark solids. We compare the convergence of the tri-3 C 1 splines to (non-interpolatory) Catmull-Clark solids.
Fixed-grid immersed representations, such as web splines [42] or unstructured collections of radial basis functions [43], require careful adaptation of computations near the implicitly enforced boundaries. Penalty methods, e.g., in [44], can add smoothness constraints as part of the solution process at the cost of increasing the size of the problem. The approach requires a judicious choice of penalty parameters.

Singular Jet Collapse Constructions in Two Variables
There are three types of singular spline spaces in the bivariate case: (1) collapse of the domain as for generalized subdivision, (2) collapse of an edge or face as for polar constructions or Duffy-type [8] elements with a removable singularity, or (3) collapse of the set of derivatives (jets) at irregularities in the grid. This review focuses on the third option. Singular corner constructions that collapse the Taylor expansion (1-jet) at the irregular point by setting derivatives to zero have been proposed by [4][5][6][7]9,[45][46][47][48]. The induced singularity side-steps the vertex enclosure problem [34,49], a non-trivial algebraic requirement that arises from forcing the mixed parametric derivatives ∂ u ∂ v f and ∂ v ∂ u f to agree. Ref. [4] suggested to simply set the mixed derivatives to zero, and [47] proved that if higher partial derivatives are suitably constrained, such singularities are locally removable. That is, the parametric singularity does not result in a loss of geometric smoothness. In contrast to subdivision, the singular corner approach yields a finite number of polynomial pieces compatible with existing CAD modeling environments. The resulting C 1 surfaces typically have poor shape when used for free-form design. A recent variant, proposed in [9], removes visible shape defects but at the cost of an increased polynomial degree and overall complexity. After a gap of 20 years, [6] re-discovered the usefulness of singular parameterization to circumvent the challenges of refinement for engineering analysis and functions on irregular bivariate 2-manifolds. Combining the singular splines with PHTsplines [50] (also known as bi-cubic finite elements with hanging nodes in the finite element literature) yields a bi-cubic C 1 space with adaptive refinability. Ref. [6] demonstrates the space's effectiveness for modeling and solving thin plate challenge problems of the finite element obstacle course, such as the 'octant of a sphere' and the 'Scordelis-Lo roof'.

Constructions in Three Variables
Ref. [10] base their spline space on tri-cubics but need no jet collapse since the space is only C 0 across extraordinary edges and vertices. Their tri-cubic C012 splines have three types of degrees of freedom: C 2 spline control points (mesh vertices), C 1 spline control points (8 per box), and individual BB-coefficients near the irregularities (64 per box). Ref. [11] built a tri-3 C 1 spline space with singular parameterization. The space offers eight degrees of freedom per box.

Tri-3 C 1 Splines on Unstructured Box-Complexes
This section gives a brief summary of the tri-3 C 1 splines defined in [11], starting with the definition of an unstructured box-complex. The spline space has 2 3 basis functions per box, see Figure 2a. By default, the map x : R 3 → R 3 that defines the geometry of the physical domain, is initialized by interpreting the vertices of the box-complex, wherever possible, as B-spline coefficients [22,23]. Knot insertion (averaging) converts the C 2 spline coefficients into C 1 spline coefficients (of a C 2 function). Then, at each irregularity, a wellbehaved linear function is determined and composed with a singular local volumetric reparameterizationx consistent with the local layout of the box-complex. All first derivatives ofx are continuous, albeit zero across irregularities. However, since the inversex −1 is well defined, the local expansion of the linear function composed ofx can be reparameterized to remove the singularity. The polynomial pieces of the spline space, therefore, join not just nominally C 1 (with a singularity), but smoothly over the whole box-complex. The tri-3 C 1 spline space can reproduce linear functions and is refinable. Each box with an irregularity is dyadically split into 2 3 sub-boxes to localize the operations that make the space C 1 . Splitting allows irregular points to be in close proximity without interfering with one another and simplifies the space's use for computations: every input box, regular or irregular, uniformly contributes exactly 2 × 2 × 2 degrees of freedom. Notation and Indexing. Analogous to a simplicial complex, a box-complex (also known as a hex-mesh) in R 3 is a collection of d-dimensional boxes, 0 ≤ d ≤ 3, called d-boxes. Boxes of any dimension overlap only in complete lower-dimensional d-boxes. A 0-box is a vertex, a 1-box an edge, a 2-box a quadrilateral, and a 3-box is a quadrilateral-faced hexahedron. A box without a prefix is a 3-box. Irregularities. For d < 3, an interior d-box is regular if it is completely surrounded by 2 3−d boxes and for 3 >d > d, all incidentd-boxes are regular. For example, for a vertex to be regular, all edges incident to it must be regular. In R 3 , a regular vertex (d = 0) is surrounded by 8 boxes, a regular edge (d = 1) by 4 boxes, and a regular quadrilateral face (d = 2) by 2 boxes. Interior faces are always regular since they are shared by exactly 2 1 boxes.
Within stacked bi-variate irregularities, we define a 0-box to be a semi-regular box (blue point in Figure 2b) if (i) the box is shared by exactly two edges that are each surrounded by the same number of n e = 4 boxes and (ii) the box is fully surrounded by 2n e boxes. An edge connecting two semi-regular points is a semi-regular edge. A 1-box that is not a semi-regular edge but is surrounded by n e = 4 boxes is called an irregular edge. A point that is neither regular nor semi-regular is an irregular point.
Boxes with at least one irregular point are evenly split into 2 m sub-boxes so that the sub-boxes contain at most one irregular point each. A box is regular if all its vertices are regular. Otherwise, the box is irregular. Example 1. Two layers of five boxes share a semi-regular point (blue in Figure 2b). If the top and bottom 5-valent points are not semi-regular, the two edges forming the axis are irregular. If the stacked configuration were to continue to a third layer of five boxes, the middle edge would be a semi-regular edge. The dotted lines in Figure 2b hint at the partition of the C 1 spline into polynomial pieces near edge irregularities. Figure 2c illustrates an irregular point enclosed by four boxes.

Solving Poisson's Equation over Unstructured Hex Meshes
Ref.
[10] uses a supercomputer implementation to solve Poisson's equation to analyze mechanical models similar to those in Figure 1. Their tri-variate splines are akin to [6] and, therefore, to the tri-3 C 1 splines in [11] but differ both in that they are truncated to transition from irregular to regular regions of C 2 tri-3 splines and require no jet collapse because they are C 0 near irregular edges. For the Poisson equation, their numerical tests yield straight lines with a slope consistent with optimal convergence, e.g., O(h 4 ) for the L 2 error. (They also point out that pure C 0 splines have too many degrees of freedom and that (finite) C 1 splines over irregular box-complexes, i.e., tri-3 C 1 splines, still await study in the literature.) For Catmull-Clark solids, combined with [53], Ref. [10] reports sub-optimal convergence.
To ensure the correct implementation of the tri-3 C 1 splines and the overall computational framework, and since the solution of second-order elliptic equations is of interest in its own right, we test the convergence of the tri-3 C 1 splines to the solution of the Poisson's equation over the physical domain Ω = x( ): The weak form of Poisson's equation, projected into the C 1 space of basis functions B j , i.e., Galerkin's approach, is This can be rewritten as the matrix equation Kc = f to be solved for the coefficient vector c of u h where Here, the sum is over all pieces α where B i has support and J α := ∇ s x α . Boundary constraints are enforced following [54] (Section 1): we split the set of basis functions {B i } into two sets: I := {B 0 , · · · B n } and B := {B n+1 , · · · B n+n δ }, with where the basis functions in I vanish on the boundary while those in B are nonzero on the boundary and vanish towards the interior. We set the coefficients of the functions in B to enforce the boundary conditions, leaving the basis functions in I free to solve the system of equations.
Then the exact solution of Equation (2) is Starting with the tensor grid as a basic 'sanity' test, we refine each of the three types of meshes, solve, and measure the errors. The error plots in Figures 5 and 6

Comparison to Catmull-Clark Solids
For further comparison and to calibrate with respect to a competing, easy-to-implement (in) finite element approach, we coded Catmull-Clark solids following [20,21]. Catmull-Clark solids are assumed to be C 1 at irregular points C 1 across irregular edges and C 2 elsewhere. In the bivariate setting, [56] reported the convergence of O(h 2 ) and [57] reported an L 2 convergence of O(h 2.5 ) for the 35-mesh (one slice of the 35-extruded mesh). For trivariate interpolatory Catmull-Clark solids, [41] (Figure 7) observed an initial decay in the L 2 error of O(h 3.52 ), with the error concentrated at irregularities. Catmull-Clark solids are, therefore, a natural competing representation for functions over unstructured box-complexes. Here the vertices of the input mesh are the degrees of freedom (as opposed to the 2 3 control points per box of the tri-3 C 1 spline). On sub-complexes without irregularities, the function can be evaluated using de Boor's algorithm. Otherwise, we locally perform Catmull-Clark solid subdivision until each Gauss point is enclosed by a regular neighborhood so that we can evaluate. The coefficients of the outer boundary of the mesh are set to enforce (zero) Dirichlet boundary conditions. Figure 7

Solving the Biharmonic Equation over Unstructured Hex Meshes
One of the applications of the biharmonic equation is the stream function formulation of Stokes and Navier-Stokes equations. A classic approach, e.g., [58][59][60], is to rewrite the biharmonic fourth-order partial differential equations as a system of first-order equations. Numerical solutions then typically require post-processing to yield a proper higher-order approximation. Using a finite difference multigrid approach, [61] report solving a tensor product grid of size 512 3 . Ref. [62] model a fourth-order Cahn-Hilliard flow on a regularly partitioned cube. This and the T-spline constructions [63,64] do not discuss convergence rates. The regular tensor grid and C 2 splines [19] predict a convergence of O(h 4 ). Our implementation of tri-3 C 1 spline achieves this rate, also for a geometrically displaced grid, see Figure 8.  To accommodate irregular box-complexes, discontinuous Galerkin elements with polygonal boundaries have been applied, see, e.g., [65][66][67]. Ref. [68] states optimal numerical convergence for quad meshes but provides no data for hexahedral meshes. Their approach matches boundary data via a penalty function. Below, we will compare the results of tri-3 C 1 splines to Catmull-Clark solids on unstructured box-complexes. The authors of an efficient implementation of Catmull-Clark solids on unstructured box-complexes [40] were not aware of publications that cover Catmull-Clark solid convergence for fourth-order equations. (Recall that the estimates of [18,19] rely on higher than C 1 smoothness for fourth-order problems and, therefore, do not apply.) The biharmonic equation find u : Ω → R : has the weak form that can be rewritten as the matrix equation Kc = f to be solved for the vector of coefficients c where for ξ = x(s) ∈ Ω, Here the sum is over all pieces α where B i has support. By the chain rule and denoting the kth column of the inverse of the Jacobian J of x and s := x −1 (ξ) as J −1 (:,k) , For the convergence measurements, we chose Ω := [0, 6] 3 partitioned as shown in Figure 4, and chose f so that the solution is For Catmull-Clark solids we observe sub-optimal convergence in the first steps already on the perturbed tensor grid, see Figure 9. On the 35-extruded mesh and ball octant, the convergence drops to sub-quadratic. Figure 10a indicates tri-3 C 1 splines O(h 3 ) convergence in the L 2 norm and optimal convergence in the H 1 and H 2 norms for the 35-extruded mesh. For the ball octant, see Figure 10b, the convergence deteriorates notably in the third and fourth refinements. Since the computations do not betray any of the usual signs of a software bug and have been carefully and repeatedly checked-Indeed, the computations still converge-this observation points to a gap in the published theory of numerical convergence that, while intriguing, lies outside the scope of the present investigation into the numerical properties of tri-3 C 1 splines over box-complexes. As a further observation, the error on the tensor grid is maximal near high curvature. For 35-extruded mesh and ball octant, the error is maximal near, but not at, the irregularities. Finally, since the integrals in matrix K and the right-hand side f are computed using Gauss points interior to each polynomial domain, K and f are well-defined, regardless of smoothness. However, when the space of functions permits jumps in the derivative that the weak formulation of the biharmonic equation does not allow, the solution of linear Equation (7) may not be a solution to the Galerkin projection of the biharmonic equation. Nevertheless, using the same C 1 parameterization x of the physical domain geometry but relaxing the space of analysis functions u h to be only C 0 -by not applying the projection P-also yields near-optimal convergence rates for the irregular meshes, see Figure 10c,d.

Conclusions and Future Work
When solving the Poisson equation with tri-3 C 1 splines, the convergence rate is optimal, i.e., in line with the regular tensor-product case for the test problems. For biharmonic equations on unstructured hex meshes, the isogeometric approach using tri-3 C 1 spline converges to the correct solution. Intriguingly, the observed convergence rate in the irregular case is less than in the regular case, and an implementation error is very unlikely.
Combining tri-3 C 1 spline geometry with C 0 functions in the isogeometric setting recovers the convergence rate of the regular tensor-product case. This indicates a knowledge gap in the theory of a priori numerical convergence estimates. A likely explanation is that the singular, tri-variate C 1 spline space lacks full approximation power at irregular points, and this raises the question of whether and what cost-effective remedy exists and is needed to have a maximal convergence rate for fourth-order elliptic problems.
Author Contributions: The authors contributed equally to this article. All authors have read and agreed to the published version of the manuscript.
Funding: This work supported in part by DARPA TRADES and NIH EB018625.

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.
Appendix A. C 1 Projector P This brief summary condenses the longer exposition and motivation of the steps in [11]. We use the same notation. Applying P ensures C 1 continuity, via removable singularity, of the tri-3 C 1 splines corresponding to the irregular sub-boxes H s , s = 1, . . . , n. Let c s α ∈ R 3 be coefficients of H s with c s 0 the central irregular point. Denote by α ∈ T the labels of the direct neighbors of c s 0 and by α ∈ G the labels of the 2-neighborhood of c s 0 . The values y α for α ∈ T are collected in the n-vector y T .

1.
Compute a best-fit linear map to the BB-coefficients c s α , α ∈ T, e.g., by computing a vector as
where I 3 is the 3 × 3 identity matrix, e j (i) = 1 for i = j and zero otherwise are the (local unit) labels of the three directions emanating from the irregular corner of H s , n k is the valence of the edge with label e k and, C k ∈ R 2×3 are the first two rows of Q n k (c 2e k +e j − c 2e k ), j = k where Q n k is defined in (A5). Check that the system is well-conditioned, i.e., the sub-box is well-formed according to [11] (Def 2).