Isogeometric analysis with $C^1$ cubic Powell-Sabin splines

1. Problem Statement

The paper addresses the challenge of bridging the gap between Computer-Aided Design (CAD) and numerical simulation (Isogeometric Analysis, or IGA) for complex geometries, specifically planar and spatial surface domains.

• Limitations of Current Methods:
  • Finite Element Analysis (FEA): Traditionally uses $C^0$ continuous elements (like cubic Lagrange) on triangular meshes. While flexible for unstructured grids, they lack the higher-order smoothness ($C^1$ or higher) required for solving fourth-order PDEs (e.g., biharmonic problems) without complex mixed formulations or penalty methods.
  • Standard IGA: Typically relies on tensor-product B-splines or NURBS. While offering high smoothness, these are restricted to four-sided (quadrilateral) domains. Handling complex, multi-patch, or unstructured geometries often requires "extraordinary points" where smoothness degrades to $C^0$, or complex multi-patch gluing techniques.
• The Goal: To develop a robust IGA framework using $C^1$ cubic Powell–Sabin splines on unstructured triangulations. This approach aims to combine the geometric flexibility of triangular meshes with the high-order smoothness ($C^1$) and rational representation capabilities of NURBS, enabling the solution of both second-order (Poisson) and fourth-order (biharmonic) boundary value problems on arbitrary surface domains.

2. Methodology

The authors propose a comprehensive framework involving three main stages:

A. Mathematical Foundation: $C^1$ Cubic Powell–Sabin Splines

• Refinement: The method starts with an unstructured triangulation $\Delta$ of a domain $\Theta$. A Powell–Sabin refinement ($\Delta_{PS}$) is applied, splitting each triangle into six smaller triangles by introducing split points (triangle centers and edge midpoints).
• Spline Space: The authors utilize the space $S_3(\Delta_{PS})$, consisting of piecewise cubic polynomials that are globally $C^1$ continuous and possess $C^2$ super-smoothness at specific split points.
• Basis Functions: The space is spanned by B-spline-like basis functions associated with vertices and edges of the original triangulation.
  • 3 basis functions per vertex ($B^v_{i,r}$).
  • 2 basis functions per edge ($B^e_{l,r}$).
  • These functions are non-negative and form a partition of unity, mimicking the desirable properties of standard B-splines.
• Rational Generalization: To handle exact geometry representation (like circles and cylinders), the polynomial space is generalized to a rational space ($S^w_3$). Weights are assigned to the basis functions, creating Rational Powell–Sabin splines analogous to NURBS but defined on unstructured triangulations.

B. Geometry Parametrization (IGA Framework)

• Mapping: A geometry mapping $F: \Theta \to \Omega$ is constructed using the rational basis functions. The control points and weights define the physical surface $\Omega$ (either 2D planar or 3D spatial).
• Handling Complex Geometries: The paper demonstrates how to parametrize:
  • Planar domains with curved boundaries.
  • Annular sectors and cylindrical shells (spatial surfaces).
  • Domains with "extraordinary points" (vertices where more than 4 edges meet) without losing $C^1$ continuity, a common failure point for tensor-product NURBS.
• Boundary Conditions: The authors detail how to impose Dirichlet conditions ($u=g_0$) and Neumann-like conditions ($\nabla u \cdot n = g_1$) on the boundary. They utilize a least-squares fitting approach to determine the coefficients of the basis functions that satisfy these boundary constraints, ensuring the resulting spline space is suitable for the weak formulation of PDEs.

C. Numerical Solution

• Weak Formulation: The authors apply the isoparametric approach. The physical domain integrals are pulled back to the parametric domain $\Theta$.
• Poisson Problem (2nd Order): Solved using the standard Ritz-Galerkin method. The stiffness matrix involves the surface gradient and the Jacobian of the mapping.
• Biharmonic Problem (4th Order): Solved directly using the $C^1$ continuity of the Powell–Sabin splines. This avoids the need for mixed formulations required by $C^0$ elements. The weak form involves the Laplace-Beltrami operator applied twice.

3. Key Contributions

1. Rational Powell–Sabin IGA Framework: The paper successfully extends the polynomial $C^1$ cubic Powell–Sabin splines to a rational setting, allowing for the exact or high-fidelity approximation of complex geometries (including circular arcs and cylinders) on unstructured triangular meshes.
2. Handling Fourth-Order PDEs: By leveraging the global $C^1$ continuity, the method solves biharmonic problems directly, a task that is computationally expensive and complex with standard $C^0$ finite elements.
3. Geometric Flexibility: The approach removes the topological restriction of tensor-product NURBS to four-sided domains. It supports arbitrary triangulations, local refinement, and complex boundary shapes while maintaining high-order smoothness.
4. Boundary Condition Treatment: A robust strategy is presented for imposing boundary conditions on unstructured grids using the specific support properties of the Powell–Sabin basis functions.

4. Results

The authors validate the method through extensive numerical experiments comparing $C^1$ Powell–Sabin splines against $C^0$ cubic Lagrange elements and bicubic NURBS.

• Convergence Rates:
  • For Poisson problems, the Powell–Sabin splines achieve optimal convergence rates ($O(h^4)$ in $L^2$ norm) for both polynomial and rational cases.
  • For Biharmonic problems, the method also demonstrates optimal convergence ($O(h^4)$), confirming the suitability of the $C^1$ space for fourth-order operators.
• Comparison with Lagrange Elements:
  • Powell–Sabin splines consistently outperform or match $C^0$ cubic Lagrange elements in terms of error magnitude for the same number of degrees of freedom (NDOF).
  • The $C^1$ continuity provides a distinct advantage in accuracy, particularly for higher-order problems.
• Comparison with NURBS:
  • On spatial domains (e.g., cylindrical shells), the Powell–Sabin results are comparable to those obtained with bicubic NURBS.
  • Crucially, Powell–Sabin splines handle unstructured meshes and local refinements more naturally than NURBS, which struggle with extraordinary points and irregular topology.
• Local Refinement: In Example 8, the method successfully handles a solution with a sharp spike using local mesh refinement, demonstrating the adaptability of the triangular framework.

5. Significance

This paper presents a powerful alternative to standard IGA and FEA methodologies:

• Unification of Benefits: It combines the geometric flexibility of triangular finite elements (handling complex, unstructured domains) with the high-order smoothness of NURBS ($C^1$ continuity).
• Versatility: It enables the solution of both second and fourth-order PDEs on general surface domains without the topological constraints of tensor-product splines.
• Future Potential: The authors suggest that this framework can be extended to higher-degree Powell–Sabin splines, potentially offering even higher accuracy and efficiency for complex engineering and scientific simulations.

In summary, the work establishes $C^1$ cubic Powell–Sabin splines as a robust, flexible, and accurate tool for isogeometric analysis, effectively overcoming the limitations of both traditional finite elements and tensor-product NURBS in the context of unstructured surface domains.