A Comparative Study of Several Classical , Discrete Differential and Isogeometric Methods for Solving Poisson ’ s Equation on the Disk

This paper outlines and qualitatively compares the implementations of seven different methods for solving Poisson’s equation on the disk. The methods include two classical finite elements, a cotan formula-based discrete differential geometry approach and four isogeometric constructions. The comparison reveals numerical convergence rates and, particularly for isogeometric constructions based on Catmull–Clark elements, the need to carefully choose quadrature formulas. The seven methods include two that are new to isogeometric analysis. Both new methods yield O(h) convergence in the L norm, also when points are included where n 6= 4 pieces meet. One construction is based on a polar, singular parameterization; the other is a G tensor-product construction.


Introduction
In classical analysis, physical laws are described by a set of (partial) differential equations from which the qualitative behavior of physical systems is deduced.The differential operators used in these equations are continuous in the sense that they are based on infinitesimal change.To obtain quantitative information, one has to typically rely on computational methods.Computational methods may discretize the operator or discretize the underlying solution space of the equations.An alternative is the theory of discrete differential geometry, short DDG, which starts with a discrete description and tries to preserve key properties of the underlying continuous systems in the form of important invariants.
This paper compares convergence rates of implementations of seven different approaches for solving Poisson's equation on the disk.Besides the DDG approach [1,2], the comparison includes two classical C 0 and C 1 finite (Hsieh-Clough-Tocher) elements and four flavors of isogeometric analysis (IgA).IgA is a form of iso-parametric analysis (see Section 3) using higher-order elements such as tensor-product B-splines both to describe the domain and the approximate solution of a partial differential equation.IgA currently has some limitations, foremost sub-optimal numerical convergence rate where the spline elements are not laid out regularly, i.e.where they are not arranged as quad-grids or a hierarchical refinement thereof [3][4][5][6].Choosing Poisson's equation on the disk as the model problem, forces the introduction of irregular mesh points.
The paper's contributions to the state-of-the-art are • a first qualitative comparison between classical finite elements, a DDG approach and four isogeometric constructions; the key outcome is presented in Fig. 10, page 12; • an investigation of quadrature formulas for subdivision IgA finite elements; • implementation of an IgA method for C 1 functions on complex domains that is based on G 1 constructions and yields O(h 3 ) convergence also at irregular points; this improved convergence is confirmed for an L-shaped domain and for an elastic plate with circular hole.
• implementation of an IgA method with singular parameterization at irregular points that yields O(h 3 ) convergence also at irregular points.
Overview Section 2 gives an overview of the two classical finite element spaces and the canonical DDG approach to Poisson's equation.Section 3 describes four, partly new approaches to constructing C 1 functions over complex domains by using singular, respectively geometrically continuous splines.Section 4 succinctly reviews the classical variational framework for the Poisson equation.Section 5 compares numerical convergence rates, as summarized in Fig. 10.

Classical finite elements and DDG
This section briefly reviews standard non-linear finite elements and the cotan-formula-based DDG approach.

C 0 quadratic triangular elements
Also known as Linear Strain Triangle (LST) or Veubeke triangle, the C 0 quadratic triangular element was developed by B. M. Fraeijs de Veubeke [7].The six degrees of freedom of a polynomial of total degree 2 in two variables can be expressed as the coefficients c ijk of the polynomial in total degree Bernstein-Bézier form (BB-form, see e.g.[8]): Fig. 1 shows the two types of C 0 quadratic elements, one associated with a vertex, the other with the mid-edge of a triangle.

Hsieh-Clough-Tocher Elements
The Hsieh-Clough-Tocher (HCT for short) element is a classical C 1 finite element (see e.g.[9]).The HCT-element is piecewise polynomial of degree There are other methods of building smooth finite elements on split triangles, for example the Powell-Sabin construction [10].The general theory is presented by Lai and Schumaker [11].Recently, these classical elements have been applied to IgA, for example [12] or [13].

The discrete differential geometry approach
The theory of discrete differential geometry, short DDG, starts with a discrete formulation that strives to preserve key properties of the underlying continuous systems in the form of important invariants.An example of a DDG operator is the cotangent formula for modeling the Laplace-Beltrami operator (see e.g.Pinkall and Polthier [1]).In applications, DDG generalizes the principles underlying the continuous operator to make methods directly applicable to the data and to improve robustness over just discretizing the continuous operators.
The cotangent formula discretizes the Laplace-Beltrami operator on a triangular mesh.Among the desirable properties for discrete Laplacians enumerated in [14], we are mainly concerned with convergence in the sense that the discretization solves, in the limit under refinement, the PDEs correctly.In [15], Desbrun et al. (see also [16]) define the cotan operator, for a function f at a vertex v i of a triangular mesh M , as where A(v) is the area of all the triangles of the 1-ring neighbors of v i , N 1 (i) is the set of the vertex indices of 1-ring neighbors, and α ij and β ij are the two angles opposite to the edge in the two triangles having the edge e ij in common (see Fig. 3(a)).In [17], G. Xu proved that (2) converges to second order to the continuous operator if each v i has valence six and v i and v j lie on a sufficiently smooth surface.More general convergence guarantees appeared in [18].In [19], K. Crane et al. derive the cotangent formula from linear finite element methods whose a basis function is shown in Fig. 3(b) and, alternatively, via discrete exterior calculus.For higher order PDEs, such as thin shell simulation [20], energy methods have been adopted.DDG has found use in computer graphics and computing for architectural geometry [20,21] and is at the heart of discrete exterior calculus (see e.g.[22] ).Formula (2) has been successfully applied to geometry processing and simulation on mesh models (see e.g.[23]).

The iso-geometric approach
The next four methods model both the physical domain Ω (the disk) and the PDE discretization space by tensor-product spline-like functions where b i := N j N k is the ith tensor-product basis function defined on a standard domain such as the unit square T .We will use polynomial splines, except at the domain boundary, where the circle is exactly expressed in rational Bernstein-Bézier form.
The four methods will be used in the framework of IgA.IgA is a special case of the classical isoparametric analysis.The term iso-geometric analysis, short IgA, was coined by T. Hughes et al. [24] in an effort to eliminate the representation gap between computer aided design and engineering analysis.In particular, IgA proponents have advocated the B-spline representation [25] both for modeling the geometry of Ω and for presenting the bases b i of the differential geometric analysis.
Then the space of functions on Ω(T ) is obtained by composing test functions, also defined on T , with the inverse of x α (see Fig. 4b).In IgA, we use test functions (b i ) α : T → R, where (b i ) α is the part of the ith basis function b i on the domain piece defined by Ω α .That is, the test functions are drawn from the same space as x α .Then the discretization space on Ω is the span of the functions (see Fig. 4b).
Figure 5.The non-smooth C 0 bi-3 basis function: (a, left) A quad mesh of points associated with B-spline-like functions, for n = 5. (a, right) the coefficients of the patches in tensor-product BB-form defined by C 2 extension of the regular spline complex towards the extraordinary point (see [29]).Note the n = 5 points of valence 3 surrounding the central n-valent point.A main open challenge of IgA are extraordinary points where more or less than four tensor-product patches join.(The analogue in the case of three-sided patches is to have more or less than six elements meet at a point.)Such points occur for topological reasons, by Euler's count and are often inserted to better adjust the mesh to the local geometry.This is illustrated in Fig. 5(b) where a central n = 5-valent (5-neighbor) point is surrounded by n 3-valent points.Without a sophisticated treatment of extraordinary points, the advantage of high-order methods in IgA may be nullified by slow convergence near these points.There is an ongoing, vigorous discussion of the proper choice of refinement space for hierarchical adaptive modeling [3][4][5][6] but this does not address the modeling at extraordinary points.
Our four constructions represent alternative approaches to deal with the extraordinary points.We focus on higher-order spline-like representations that mimic bi-cubic (bi-3) tensor-product splines.Our first choice is the space of polynomials of degree bi-3 that are C 2 except near the extraordinary points where they are only continuous.We will see that this simple space has good convergence except at the extraordinary point.This flaw motivates and makes the case for our other three constructions.Our second construction leverages Catmull-Clark subdivision.This is inspired by the seminal work of F. Cirak et al. [26,27] who used subdivision surface functions over triangulations for thin shell analysis.The challenge, not emphasized in the original work, is the choice of integration rules (see Section 5.1).Our remaining constructions come from geometric modeling and are new to IgA.The third construction leverages everywhere G 1 functions (that are C 1 when considered over the physical domain).Here G 1 refers to geometrically continuity, i.e. matching of derivatives after reparameterization of one or both of the adjacent function pieces (see [28]).Just as for Catmull-Clark subdivision, the space of generating functions consists of C 2 polynomial splines of degree bi-3 away from extraordinary points.Finally, we introduce C 1 functions with polar layout, i.e. with a central pole or singularity.The last three constructions allow us to address high-order PDEs, such as Kirchhoff-Love shell model or buckling analysis that are not a focus of the present exposition.At extraordinary points this interpretation of the quad mesh breaks down.Assuming that extraordinary points are isolated, in the sense that no two extraordinary points share a quad, we can still construct a bi-3 C 2 tensor-product spline complex with 'holes' where n = 4 quads meet.A simple way to complete the spline complex is to extend the existing spline patches C 2 into the holes, as n patches in bi-3 tensor-product BB-form (see Fig. 5, (a,right)).These patches are defined up to just one BB-coefficient, corresponding to the position at the center of the hole and common to all n patches.This coefficient is trivially set to the average of the surrounding coefficients.The result is bi-3 elements that form standard C 2 bi-3 spline complex away from the extraordinary point, and that join C 0 at the extraordinary point (see Fig. 5(b)).

Catmull-Clark Elements
Subdivision splines are piecewise polynomial splines with singularities at the extraordinary points [30].The neighborhood of an extraordinary point is an infinite sequence of nested spline rings (where 'ring' indicates the connectivity, not an algebraic ring).Subdivision splines have been used for finite element analysis well before the advent of IgA (see [26]), but did not receive the attention from the engineering community that IgA is currently generating.A more complicated framework for adaptive simulation with subdivision splines was introduced by E. Grinspun et al. in [31].Subdivision-based functions for IgA on solid models were presented in [32,33].Catmull-Clark subdivision has been used in [34], the similar bi-2 spline-based Doo-Sabin subdivision in [35] and Loop's subdivision in [36], for large deformation and anisotropic growth.
Among a myriad of subdivision schemes, Catmull-Clark dominates in industrial implementations.Catmull-Clark subdivision refines a mesh by binary split in each direction (see Fig. 6).The basis function associated with a vertex not on the boundary has support on 2-ring neighbors.(For spline surfaces with boundary, we apply 'natural end conditions', i.e. do not evaluate under-defined outer quadrilaterals, after extrapolating the existing mesh.)Recently, Barendrecht performed experiments of numerical convergence of IgA with Catmull-Clark surfaces and observed poor convergence near extraordinary points [37].He conjectured that this is due to the well-known unbounded Gaussian curvature of Catmull-Clark subdivision at these points.Based on our experiments in Section 5, Table 1, we think that the poor numerical convergence is primarily due to the application of Gauss quadrature rules with respect to the original quads, rather than choosing quadrature points for each sub-polynomial of sufficiently many levels of refinement (see Fig. 6a,b).

Higher-order G 1 elements
A technique from the applied mathematics area of geometric design allows glueing-together function pieces with 'geometric continuity'.The result is a C k manifold.The designation G k is used to emphasize that derivatives of adjacent patches match only after reparameterization [28].Geometric continuity then allows smoothly joing n = 4 tensor-product pieces (often called patches) in the sense of parametric surfaces.Composition of a G 1 construction with the inverse of a G 1 construction that share the same reparameterization yields a C 1 function.There are many G 1 constructions in the literature.Some naturally complete a bi-3 C 2 tensor-product spline complex with bi-3 patches to a G 1 structure.To avoid splitting quadrilateral domains, we chose [38], a simplified version of [29] that deploys n bi-5 patches at the extraordinary point.The corresponding basis functions are shown in Fig. 7.The resulting G 1 elements are C 2 at the extraordinary point.The additional smoothness at the extraordinary points guarantees high polynomial reproduction and hence high numerical convergence also at the extraordinary points.

Polar C 1 elements
The polar parametric surface construction of [39] provides a simple element that is smooth and particularly well-behaved at points where many surfaces join in a triangle fan at the center of the disk.To match the tensor-product standard, the triangle of the fan can be interpreted as quadrilaterals that have one edge collapsed.Analogous to the G 1 edges in the previous construction, the central singularity is no cause of concern for either shape or numerical convergence.Note that this observation matches the recent results of Takacs and Jüttler [40] who analyze singularities at domain boundaries, where the test functions by themselves are not well-defined.

Solving the Poisson equation
We are solving Poisson's equation on the domain Ω, subject to zero boundary conditions on the boundary ∂Ω of Ω: The DDG method discretizes this formulation directly as −∆ M u = f where the operator ∆ M has been defined in (2).For all six methods other than DDG, we solve the equation numerically by considering its weak form: We seek an approximate solution in terms of the generating functions b i : Ω → R defined in ( 5) by determining the coefficients c i ∈ R in Using Galerkin's method, we set v = b i in (7) and obtain the constraints

*)
This yields a system of linear equations and the vector of unknown coefficients is c : For all iso-geometric methods, we define the physical domain Ω as in ( 4) and write the integrals in (9) as a sum of integrals restricted to some x α (T ).Using (5) and dropping the subscript α, we can, by the change of variables, express the local integrals with respect to each parameter domain T as where J is the transpose of the Jacobian of the mapping x : (s, t) ∈ T → [x(s, v) y(s, v)] t .For implementation, we collect Similarly, for the right hand side term,

Numerical results and comparison
Before we compare the convergence rates of the methods for Poisson's equation on the disk, we need to look in more detail at the quadrature rules that are used for Catmull-Clark functions in the IgA setting to compute (10) and (12).
Table 1.Error (scaled by 10 −5 ) of the computed solution of Poisson's equation by Catmull-Clark subdivision on Disk 1 (see Fig. 6), for different levels of subdivision when applying Gauss quadrature.The subdivision is localized to not refine the overall mesh.The p-point Gaussian quadrature rule is known to exactly calculate the integral of polynomials of degree up to 2p − 1.For piece-wise polynomials, Gaussian quadrature only gives approximate results.Table 1 shows that in order to obtain good integration results at irregular points, one needs to apply exact Gauss quadrature on many subdivision layers to obtain convergence.We found that subdivision of depth 7 was necessary for results to stabilize.A more principled approach is to take advantage of the recursive nature of subdivision and compute the quadrature rules via eigendecomposition as in Halstead et al. in [41, App B].

Convergence rates
Fig. 9 shows the three types of meshes that one might naturally associate with the methods.For C 0 quadratic, HCT and DDG elements, we optimized the aspect ratio of the triangles to guarantee numerical stability.For bi-3 C 0 , Catmull-Clark and G 1 bi-3/bi-5 elements, we chose a central extraordinary point with valence n = 5, surrounded by n satellites of valence 3. Other n can be tested or the singularities can be distributed to the boundary as in Takacs and Jüttler's approach [40].Finally, the polar configuration is natural for the polar C 1 elements.The convergence of polar elements is remarkably unaffected by the valence of the central point.We choose f := 1 in (6).Then the exact solution is u := (1 − x 2 − y 2 )/4 and we can display the exact error.Fig. 10 confirms at least an O(h 2 ) convergence for all higher-order methods as well as for DDG.The graphs are in log-scale with smaller mesh spacing displayed to the left.That is, the entries to the left correspond to more elements.Fig. 10 displays the hoped-for O(h 3 ) L 2 -norm convergence of the iso-geometric approaches.(The theory of Bazilevs et al. [42] shows that O(h 3 ) L 2 -norm convergence is optimal for regular meshes and degree bi-3 splines.)The spread factor between the polar and the other three methods is five.That is, in the L 2 sense, the easily implemented C 1 polar iso-geometric approach (that is natural for the disk) is superior.Remarkably, though, the more general G 1 construction excels in minimizing the L ∞ .The spread in the error between the four O(h 3 )-convergent methods is more than a factor of 8.The higher L ∞ error of Catmull-Clark elements Fig. 11a and C 0 bi-3 elements Fig. 11c is Figure 10.Convergence comparison between methods.Note that the graphs are in log-scale (the triangle indicates the convergence exponent in log-scale) and that higher mesh density is to the left, as the mesh spacing on the abscissa decreases.

Complexity
We do not compare execution times since implementation details, such as memory management on the GPU, pre-tabulation of basis functions or sparsity (taking advantage of finite support), etc. strongly influence the performance.However, we can state the size of the matrix K in (9) for the IgA methods.
For Disk 1, the mesh Fig. 9(e) used by bi-3 C 0 , Catmull-Clark and G 1 bi-3/bi-5 elements.K is a matrix of size 151 × 151 for 120 patches.For Disk 1, the mesh Fig. 9(f) used by C 1 polar elements K is a matrix of size 101 × 101 for 80 patches.The relative times for solving equation ( 9) by our unoptimized implementations showed a ratio of roughly 4:6:8:26 for C 1 polar, bi-3 C 0 , G 1 bi-3/bi-5 and Catmull-Clark elements respectively.We surmise that, in the natural disk setting, C 1 polar elements can achieve good results with fewer elements and fast computation.Remarkably, the quality of the polar method does not depend on the valence of the center point: the polar method's L 2 error decreases with order of O(h 3 ) even though the valence of the center point is doubled with each refinement step.The three non-IgA-methods, and the DDG method in particular, have a lower memory requirements.This allowed us to add very fine meshes for the comparison in Fig. 10.For Disk 3, the mesh of type Fig. 9(d) has 6144 triangles, and 3169, 12481 and 18819 degrees of freedom for DDG, C 0 elements and HCT elements respectively.We observed solution times with ratio 2:10:37, making DDG attractive in comparison to the classical HCT elements.We confirm the high-order convergence of the new G 1 bi-3/bi-5 elements by computing two additional well-known benchmark problems.The first is Laplace's equation on the L-shaped domain: The exact solution is u(x, y) = r 2 3 sin(2a/3 + π/3) where r(x, y) := x 2 + y 2 and a(x, y) := atan(x/y).We use this exact solution to provide non-homogeneous Dirichlet boundary conditions for the numerical problem formulation.Fig. 12(b) shows the difference between the known exact and the computed solution when solving the Laplace problem on the three mesh resolutions of Fig. 12(a).Predictably, the largest errors occur at the corner singularity.Fig. 12c,d show the convergence in the L 2 and in the L ∞ norm respectively.

Conclusion
When starting our work on IgA methods, we found reports on many individual implementations and applications.To get a sense of how IgA methods stack up against each other as well as against some of the more classical and the DDG methods, we implemented these methods.Our goal here was to confirm qualitative behavior since performance comparisons would likely depend on implementation details.Moreover, methods with the O(h 3 ) convergence at the extraordinary point are currently missing in the IgA literature.By introducing two methods with improved polynomial reproduction and smoothness at the irregular points, we were able to improve L ∞ convergence at the extraordinary point.The purpose of the paper is to share our experience concerning the qualitative behavior of implementations of the seven methods.
While neither of the two classical approaches can directly be applied to surfaces as physical domains, four of the remaining five generalize directly.Of the four IgA constructions, the G 1 bi-3/bi-5 construction needs additional work to guarantee compatible surface representations.We applied the methods to generate geodesics on surfaces by solving the heat equation.These applications confirmed the convergence characterization of Fig. 10.
Three of the IgA approaches, subdivision, G 1 bi-3/bi-5 and C 1 polar, as well as some extensions of the DDG approach, span the correct space to solve thin shell and biharmonic equations.Here we are collecting further comparative data.

3 .
It is constructed over a triangle domain split into three sub-triangles by connecting the vertices to the barycenter.The resulting 3-piece C 1 function has 12 degrees of freedom that one may choose as the value and first derivatives at each vertex, plus normal derivatives on the midpoint of each edge.The twelve basis functions overlapping a triangle are constructed by setting one of the degrees of freedom to 1 and all others to zero.This is most conveniently expressed in total-degree 3 Bernstein-Bézier basis functions b k .We associate three basis functions b 3i+j , j = 0, 1, 2 with value and the partial derivatives at each vertex v i , i = 1 . . .n and one basis function b 3n+k with each edge e k .The functions b 3i+j have support on the triangles with common vertex v i and the functions b 3n+k have support on the triangles sharing e k .

Figure 2 .
Figure 2. HCT basis functions (top view with height scale)
Notation of (2) (b) Top view of the linear 'hat' function

Figure 4 .
Figure 4.A basis function b i is the composition of the basis function b i on the tensor-product parameter domain T and the inverse of the geometry mapping x.(a) shows the union of 4×4 domains T and its image Ω under 4 × 4 maps x α .
(a) control net and C 2 extension in BB-form (b) C 0 bi-3 basis function

Figure 7 .
Figure 7. G 1 element at an extraordinary point.
Two G 1 bi-3/bi-5 basis functions (b) G 1 bi-3/bi-5 basis function with on patch in BB-form lifted up

Figure 8 .
Figure 8. Polar elements for polar configurations

Figure 9 .
Figure 9. Three types of meshes.Columns a, b, c show refinement by halving h, hence quadrupling the number of elements.Rows d, e, f show the partitions specific to each of the three classes of methods.
Error in L 2

1 (
b) Error in L ∞ concentrated at the central point, as large spikes.This is to be expected since neither method reproduces all quadratic expansions at the central point, a fact that is also reflected in their lack of C 1 , respectively C 2 smoothness at the central point.

Figure 11 .
Figure 11.Poisson's equation on Disk 1: difference graphs between the exact solution and the computed solution.
h-refinement of L-shape (b) Difference between the exact solution and the computed solution.