An intrinsic finite element framework for scalar- and vector-valued partial differential equations on general manifolds

1. Introduction

• Extrinsic methods incur geometric consistency error unless geometry order increases with FE order, while vector-valued methods cannot generally enforce tangency and H1-conformity simultaneously.The literature is identified as lacking an intrinsic compatible FE framework for vector-valued unknowns on multi-chart atlases.Source: §1 · pp. 1–2
• The framework formulates weak problems and discretisation in atlas parametric space, with exact geometry entering through the metric tensor at quadrature points.The presentation uses vector calculus notation to align with the data structures common in FE codes.Source: §1 · p. 2
• Arbitrary-order spaces for the whole discrete de Rham complex are constructed on multi-chart atlases, with topological gluing preserving the complex’s discrete vector-calculus structure.Compatible FE spaces inherit discrete identities associated with local mass conservation, steady geostrophic modes, and suppression of spurious modes.Source: §1 · pp. 1–2
• A nodal change of basis couples vector-valued degrees of freedom across incompatible chart coordinates, yielding exact tangency and interface jumps of approximation-error order.Inter-chart continuity holds at the nodes; rigorous convergence analysis remains open.Source: §1 · p. 2
• The intrinsic pipeline supplies the metric analytically at quadrature points and computes in parametric space; the framework is implemented in GridapGeosciences.jl.The paper compares this pipeline with an extrinsic ambient-space method and notes support for distributed-memory parallelism.Source: §1 · p. 2
• Cubed-sphere tests capture the metric to machine precision with sufficient quadrature, conserve shallow-water mass exactly and energy up to time-discretisation error, and show optimal observed surface Stokes convergence rates.The nodal vector construction has no rigorous convergence analysis.Source: §1 · pp. 2–3

2.1. Manifolds, charts and fields.

• A manifold is covered by smooth chart maps, and tangent fields are represented in each chart by contravariant components that transform across overlaps.The metric converts those components to covariant components; the chart transition map identifies overlapping parametric points.Source: §2.1 · pp. 3–4
• The metric encodes tangent-field inner products and supplies the geometric quantities used in intrinsic formulations and integration on parametric domains.The framework uses g, g^-1, and √g; geometry can be evaluated at quadrature points without forming the ambient differential dσ.Source: §2.1 · p. 4

2.2. Intrinsic differential operators.

• Intrinsic gradient, divergence, curl, and Laplace–Beltrami operators act on parametric fields, while their metric-dependent terms cancel in the stated annihilation identities.Those identities hold pointwise for any metric and are inherited exactly by the discrete spaces.Source: §2.2 · p. 5

2.3. Pullbacks, function spaces and weak identities.

• Pullbacks relate surface fields and operators to parametric representations, yielding weak identities from Euclidean integration by parts and metric-weighted pairings.The metric-free identities underpinning the compatible discretisations can be evaluated exactly; on closed manifolds, partition-interface boundary contributions cancel.Source: §2.3 · pp. 6–7, §2.4 · p. 8

2.4. The rotated de Rham complex in two dimensions.

• The two-dimensional rotated de Rham complex supplies the definitions and identities relevant to the rotating shallow water equations.Its mimetic structure is preserved exactly; non-polynomial metrics can introduce integration error in metric-weighted pairings, which the paper reports reduces to machine precision.Source: §2.4 · pp. 7–8

3. Intrinsic finite element spaces

• The framework uses standard finite element technology within each chart; the de Rham spaces require only topological inter-chart gluing.Grad-conforming vector-valued spaces additionally require transmission maps at chart interfaces.Source: §3 · p. 8

3.1. Parametric triangulations and transmission maps.

• Compatible chart triangulations match interface entities, and transmission maps transform tangent-field components between overlapping charts.The maps are needed for grad-conforming vector-valued spaces, not for the discrete de Rham complexes.Source: §3 · p. 8, §3.1 · p. 8

3.2. Discrete de Rham complexes.

• Chart-wise de Rham spaces glue through matching interface degrees of freedom to form conforming global subcomplexes at arbitrary order.The resulting algebraic identities hold exactly at the discrete level, independently of the metric and numerical quadrature.Source: §3.2 · p. 10

3.3. Grad-conforming vector-valued spaces.

• Because nodal component values are not chart-invariant, the grad-conforming construction expresses interface degrees of freedom in a common chart frame.The transition map’s transverse action is not determined by entity correspondence and orientations alone.Source: §3.3 · p. 11

3.3.1. Nodal construction.

• Owner-frame nodal degrees of freedom are transformed into local shape functions, producing fields that are exactly tangent and continuous within each chart.At every interface node, all adjacent charts define the same tangent vector; only cells touching chart interfaces are affected by the change of basis.Source: §3.3.1 · pp. 11–12
• Whole-interface single-valuedness follows when transmission maps are constant along interfaces; varying maps can leave the glued space not fully H1-conforming.The equiangular cubed sphere is an example where the transmission map varies along a chart interface.Source: §3.3.1 · p. 12

3.4. Assembly.

• Intrinsic assembly evaluates metric quantities at quadrature points on parametric meshes and does not assemble ambient-dimension objects.Quadrature accuracy affects metric-weighted integrals, not the mimetic structure of the discrete complex.Source: §3.4 · p. 13

• The section derives intrinsic finite element formulations for Hodge Laplacian problems associated with the de Rham complex and for the rotating shallow water equations.Source: §4 · p. 13
• The formulations use the glued global spaces, with chart-interface terms cancelling through the stated integration-by-parts identity.Source: §4 · p. 13
• The surface Stokes problem is treated using the grad-conforming vector-valued spaces developed earlier.Source: §4 · p. 13

4.1. Hodge Laplacian problems.

• For n = 3, the mixed scalar and vector Poisson problems use the H(div)–L2 and H1–H(curl) pairs, respectively.Source: §4.1 · p. 13

4.1.1. Scalar Poisson problem.

• The scalar Poisson discretisation uses the flux and scalar pullbacks, with metric dependence confined to the mass term and metric-free divergence terms.Source: §4.1.1 · p. 13

4.1.2. Vector Poisson problem.

• The vector Poisson formulation assumes S has no tunnels and builds its gradient and curl–curl terms from intrinsic pairings and identities.Source: §4.1.2 · pp. 13–14
• Boundary data enter the vector formulation as parametric densities; conversion from ambient data is available in closed form for the manufactured solutions.Source: §4.1.2 · p. 14

4.2. Rotating shallow water equations.

• On a closed 2D manifold, the equations evolve tangent velocity and fluid depth, with mass flux, Bernoulli potential, and potential vorticity among their associated quantities.Source: §4.2 · p. 14

4.2.1. Semi-discrete formulation.

• The semi-discrete formulation represents prognostic variables in the intrinsic finite element spaces and upwinds potential vorticity with SUPG.Source: §4.2.1 · p. 14

4.2.2. Conservation properties.

• The 2D semi-discrete system conserves mass and energy exactly regardless of metric quadrature; the 3D extension also conserves both.Source: §4.2.2 · p. 15, §4.2.3 · p. 16

4.2.3. Extension to three dimensions.

• For a 3D atmospheric shell discretised with one radial element, the formulation carries over and preserves mass and energy.Source: §4.2.3 · pp. 15–16

4.2.4. Time integration.

• Time integration uses third-order SSPRK, with a time step proportional to mesh size and stabilisation parameter ϑ = ∆t/2.Source: §4.2.4 · p. 16

4.3. Surface Stokes problem.

• The grad-conforming velocity space loses H1-conformity only along chart interfaces, while the exact intrinsic geometry ensures tangency by construction.Source: §4.3 · p. 16
• The intrinsic formulation gives exactly tangent velocity, with inter-chart continuity properties supplied by the grad-conforming space.Source: §4.3 · p. 16
• The viscous term is the Hodge Laplacian, differing by zeroth-order terms from Bochner and Boussinesq–Scriven variants, which are recovered on the sphere for divergence-free fields by adjusting α.Source: §4.3 · p. 16

5.1. Geometric accuracy under numerical quadrature.

• At sufficient quadrature degree and refinement, the scalar Poisson error settles at O(10^-14), confirming geometric accuracy to machine precision without geometric consistency error.At the coarsest level, q = 5 gives O(10^-4) error, while q = 15 gives O(10^-14); accuracy can also improve by refining at fixed q.Source: §5.1 · p. 17, Figure 2 · p. 18

5.2. Convergence of the Hodge Laplacian problems.

• The expected hp-convergence is observed for both scalar and vector Poisson problems on an atmospheric shell.The mixed formulations are assessed with manufactured solutions in the n = 3 setting.Source: §5.2 · pp. 17–18, Figure 3 · p. 18

5.3. Computational cost.

• For q = 11, the extrinsic pipeline requires nine times more operations than the intrinsic pipeline, with the difference attributed to ambient-space mappings at quadrature points.The comparison concerns the Laplace–Beltrami operator; the authors report the same conclusion for other formulation terms.Source: §5.3 · p. 19, Table 2 · p. 19

5.4. Rotating shallow water equations.

• The rotating shallow-water systems are time-integrated in two and three dimensions, with prognostic variables and potential vorticity assigned to spaces in the discrete complex.The time step scales with nominal horizontal resolution and a CFL-dependent factor.Source: §5.4 · p. 19

5.4.1. Williamson 2: steady zonal flow.

• For Williamson 2, p = 1 and p = 2 show expected second- and third-order convergence, while the p = 1 mass error reaches round-off at the finest resolutions.Semi-discrete mass conservation is exact; coarse-resolution residual error comes from the algebraic solver tolerance of 10^-8.Source: §5.4.1 · p. 19, Figure 4 · p. 20

5.4.2. Williamson 5: zonal flow over an isolated mountain.

• In Williamson 5, the two-dimensional and three-dimensional vorticity solutions are comparable; mass is conserved to machine precision and discrete energy to O(10^-4) in the two-dimensional run.The energy deviation is attributed to explicit SSPRK3 time discretisation; an implicit conservative integrator is suggested when fully discrete energy conservation is required.Source: §5.4.2 · p. 20, Figure 5 · p. 21, Figure 6 · p. 22

5.5. Surface Stokes problem.

• The surface Stokes velocity converges at rate hp+2, while its gradient and scaled interface-jump norms converge at hp+1; the jump does not degrade the measured rate.A detailed stability and convergence analysis for the surface Stokes and H1 vector Laplacian problems is left for future work.Source: §5.5 · pp. 21–22

6. Conclusions and future work

• The intrinsic framework uses the metric tensor for geometry, builds arbitrary-order conforming de Rham spaces through topological gluing, and constructs exactly tangent nodal vector fields on multi-chart atlases.Inter-chart continuity for the nodal vector fields holds at nodes, with jumps of the order of the approximation error.Source: §6 · p. 22
• Cubed-sphere experiments show Hodge Laplacian and surface Stokes convergence at expected rates, exact mass conservation, and energy conservation up to time-discretisation error for shallow water discretisations with orography.The framework captures non-polynomial geometry to machine precision at sufficient quadrature degree, and intrinsic assembly uses an order of magnitude fewer operations than its extrinsic counterpart.Source: §6 · p. 22
• Further work includes applying the atlas-agnostic framework to non-spherical and non-zero-genus manifolds, and extending it to terrain-following three-dimensional geometries.The authors also leave a priori analysis of metric quadrature and rigorous convergence analysis of the nodal vector construction for future work.Source: §6 · p. 22
• Figure 7 reports velocity and pressure L2 errors, velocity gradient error, and velocity jump error, with velocity and pressure space orders p+1 and p, respectively.The caption specifies 6 × 2^nℓ cells and parameters n = 2, R = 1, ν = 1, and α = 1.Source: Figure 7 · p. 23
• Stokes complexes whose leading spaces require C1 elements are outside the scope of this work.Source: §6 · p. 23

Declaration of competing interest.

• The authors declare that they have no competing interests.Source: Declaration of competing interest · p. 23

Data availability.

• The article contains no additional data, and the study's source code is available on Zenodo.Source: Data availability · p. 23

Appendix A. The cubed sphere atlas

• The cubed sphere is formed by projecting the six faces of a cube onto a sphere of radius R, with six atlas charts called panels.Source: Appendix A · p. 24
• Each panel uses an enlarged open square chart domain and a closed square partition domain; the equiangular projection gives the same metric on all six panels.Other geometrical maps are also possible.Source: Appendix A · p. 24
• A constant-thickness radial extrusion of the cubed sphere represents an atmospheric shell when its thickness T is much smaller than R.Source: Appendix A · p. 24
• In the orography configuration, the topography displaces the bottom of a variable-thickness shell while the top remains flat, producing metric terms proportional to the topography gradient.Source: Appendix A · p. 24