High-Order Parametric Finite Element Methods

• Temporally high-order parametric FEMs are discretization techniques that approximate both evolving geometries and dynamic PDEs with high-order accuracy using advanced time integrators.
• They couple p-th order spatial discretization with schemes like BDF and Runge–Kutta to achieve optimal error estimates and robust mesh evolution for complex, time-dependent problems.
• These methods enable precise simulation of convection-dominated flows, multiphase processes, and evolving interfaces while preserving conservation laws and structural integrity.

Temporally high-order parametric finite element methods (FEMs) are advanced discretization schemes for time-dependent partial differential equations in which both the solution and evolving domain parameters—such as moving boundaries, interfaces, or geometric features—are approximated with high-order accuracy in time. These techniques enable robust and precise simulation of problems characterized by strong temporal dynamics and geometric evolution, including convection-dominated transport, phase change, surface diffusion, multiphase flows, and geometric PDEs on evolving curves and surfaces. High-order temporal accuracy is typically realized through advanced time integrators (e.g., multi-step, Runge–Kutta, BDF, or predictor–corrector schemes), geometric updates based on continuous or piecewise high-order parametric representations, and tightly coupled mesh evolution and solution strategies. The following sections provide a comprehensive overview of foundational methodologies, error analysis, practical algorithms, and structural properties across contemporary research (Bank et al., 2013, Kim, 2013, Gawlik et al., 2014, Lou et al., 2021, Ma et al., 2021, Jiang et al., 6 Feb 2024, Jiang et al., 14 Dec 2024, Gan et al., 18 Oct 2025), focusing on key advances and rigorous results.

1. High-Order Space-Time Discretizations and Error Analysis

High-order space-time finite element methods achieve balanced spatial and temporal accuracy by employing tensor-product finite element spaces or temporally polynomial mesh trajectories. In convection-dominated problems, the essential approach is to couple a p-th order spatial discretization (e.g., via $C^0$ continuous piecewise polynomials) with a temporally high-order variational or collocation scheme. Typical formulations (e.g., method of lines or space-time slabs) yield a priori energy-norm error estimates of the form

$$\|u - u_h\|_e \leq C \inf_{v \in \mathcal{V}_h^p} \|u-v\|_e,$$

where $\| \cdot \|_e$ is a mesh-dependent norm combining $L^2$ and $H^1$ (or semi-norm) contributions in time and space (Bank et al., 2013). The optimal ("symmetric") nature of these error estimates is ensured through Galerkin orthogonality, the shift lemma, and discrete Grönwall inequalities, with additional care required for dynamic meshes and characteristic-based derivatives that absorb convection.

Mesh evolution—crucial in parametric methods—is handled through continuous mesh motion (nodes tracing piecewise polynomials in time) and discrete "rezone" steps (permitting mesh topology change at time-slab boundaries), with enforced jump orthogonality to guarantee solution continuity and preserve error bounds. This combination allows for spatial and temporal refinement primarily near evolving solution features (e.g., sharp fronts), yielding high accuracy and significant computational savings.

2. Mixed and Variational Temporal High-Order Schemes

Mixed variational principles, such as the extended Hamilton’s principle (EHP) and the mixed convolved action principle (MCAP), serve as a rigorous foundation for temporally high-order FEMs in broad classes of initial boundary value problems (Kim, 2013). These frameworks treat both primal and "mixed" variables (such as impulse or internal force) as primary unknowns, leading to coupled systems naturally amenable to higher-order temporal interpolation.

Quadratic (and higher-degree) temporal approximations in the mixed variables yield update schemes of the type

$X_n = A_p X_{n-1} + A_i f_n,$

preserving unconditional stability and, in the undamped Hamiltonian case, symplectic structure ($|\lambda|=1$ for all eigenvalues $\lambda$ of $A$). These methods outperform classic integrators (e.g., Newmark) in conserving phase-space features, reducing period elongation, and suppressing spurious numerical dissipation for large time steps. When applied to damped and forced systems, convergence and robustness are maintained, as confirmed by tests on nonlinear oscillators and seismic loading scenarios.

3. Mesh Update Strategies and Isoparametric Transformation

A central challenge in parametric FEMs is the consistent, high-order representation of evolving geometry or interfaces, especially in moving boundary or domain problems (Gawlik et al., 2014, Lou et al., 2021, Heimann, 15 Jan 2024). Several strategies are employed:

• Universal Meshes: Fixed background triangulations allow for per-time-slab submeshing and adaptive mappings ($\Phi^{(n,t)}_h$) to exactly match the moving boundary/arbitrary geometry at each timestep, leveraging boundary tracking, relaxation, and blending maps. This ensures high geometric fidelity and prevents mesh distortion under large deformations (Gawlik et al., 2014).
• Isoparametric Unfitted FEM: Level set representations are mapped with isoparametric finite element transformations ($\Theta^n$), so that even if the initial computational geometry has only low-order approximation, the mapped mesh achieves geometric errors of $O(h^{q+1})$ (Lou et al., 2021). The time derivative is discretized using high-order BDF schemes, with transfer operators ensuring that solutions from previous deformed (active) finite element spaces are reliably projected onto the current one.
• Space-Time Tensor-Product Discretizations: For problems coupling bulk and evolving surfaces, high-order isoparametric mappings are constructed in space-time, leading to tensor-product finite element spaces $V_h^{(k_s)} \otimes \mathcal{P}^{(k_t)}(I_n)$ and discrete domains updated via $\Theta_h^{st}$ (Heimann, 15 Jan 2024). Stabilization (upwind, ghost penalty) is incorporated to control temporal and geometric consistency.

These approaches facilitate high-order convergence in both space and time, with error estimates (in $L^2$ or energy norm) of $O(h^r + \Delta t^q)$ up to half-order losses due to projection or interface tracking at remeshing events.

4. High-Order Time Integration: BDF, Predictor–Corrector, and Structure Preservation

Sophisticated time-integration algorithms underpin temporally high-order accuracy. Key methodologies include:

• Backward Differentiation Formulae (BDF): Multi-step BDF schemes of order two or higher (BDF2, BDF3, BDF4) provide temporal derivatives as

$$\partial_t u^{n+1} \approx \frac{a u^{n+1} - \widehat{u}^n}{\tau},$$

where $a$ and $\widehat{u}^n$ depend on the BDF order (Jiang et al., 6 Feb 2024, Jiang et al., 14 Dec 2024, Gan et al., 18 Oct 2025). For geometric flows, BDF schemes are incorporated into weak formulations, ensuring optimal order (e.g., second to fourth order) in time while maintaining good mesh quality via predictor curves or surfaces computed by lower-order variants (often classical BGN schemes).

• Predictor–Corrector (PC) Methods: To circumvent mesh distortion and nonphysical oscillations, PC strategies utilize a predictor (e.g., via first-order scheme with $\tau/2$ time step) followed by a centered correction step for the full update. This yields second-order accurate, robust algorithms for surface diffusion and general geometric evolution, preserving mesh equidistribution without mesh regularization even under complex deformations (Jiang et al., 14 Dec 2024, Gan et al., 18 Oct 2025).
• Structure-Preserving Discretization: For conserved quantities (energy, mass, perimeter, area, volume), Lagrange multipliers and auxiliary variables are introduced to enforce discrete analogues of conservation and dissipation laws. For example, geometric evolution equations are augmented with perimeter- and area-related constraints, leading to fully implicit nonlinear schemes solved via Newton iterations and decoupling techniques (Garcke et al., 24 Aug 2024, Garcke et al., 17 Aug 2025). Scalar auxiliary variables (SAV) further enable the design of unconditionally energy-stable, high-order schemes (Li et al., 12 Jul 2024).
• Stabilized Series Expansions: In high-order time series expansion methods (e.g., Borel–Padé–Laplace), artificial diffusion is systematically added and calibrated (via recurrence relations for coefficients $\alpha_k$) to minimize error amplification and achieve discrete maximum principles, stabilizing high-order FEM approximations of transient PDEs (Deeb et al., 16 May 2024).

5. Flexibility, Generalization, and Applications

Temporally high-order parametric FEMs are broadly applicable and flexible, accommodating a wide class of models:

• Moving Boundary and Interface Problems: Universal mesh and isoparametric mapping approaches enable handling prescribed or dynamically coupled boundary evolution, including Stefan-type problems, multiphase flows, solid-state dewetting, and more (Gawlik et al., 2014, Ma et al., 2021, Gan et al., 18 Oct 2025).
• Evolving Surface and Bulk–Surface Coupled Systems: Space-time methods and level set-based isoparametric updates facilitate simulation of convection–diffusion and reaction problems on domains with evolving interfaces, as well as coupled surface–bulk transport (Heimann, 15 Jan 2024).
• Geometric Flows and Structure Preservation: High-order BDF, PC, and structure-preserving parametric FEMs have been demonstrated for mean curvature flow, Willmore flow, and surface diffusion, in both isotropic and anisotropic settings, as well as area/volume-preserving and energy-dissipating Stokes flows (Jiang et al., 6 Feb 2024, Garcke et al., 24 Aug 2024, Garcke et al., 17 Aug 2025, Gan et al., 18 Oct 2025).
• Adaptive and Exascale Simulations: High-order parametric FEMs mesh naturally with adaptive refinement strategies (both in $h$ and $p$), matrix-free implementations, and highly parallel architectures. These aspects are actively developed in exascale projects such as CEED, promoting efficient, matrix-free, partially assembled approaches for large-scale unstructured PDE discretizations (Kolev et al., 2021).

6. Summary Table: Core Techniques and Features in Temporally High-Order Parametric FEMs

Methodology	Temporal Order	Geometry Handling	Notable Features
Space-time moving FEM (Bank et al., 2013)	$p$ (arbitrary)	Continuous + discrete mesh motion	Symmetric energy-norm error, adaptive feature tracking
Mixed EHP/MCAP (Kim, 2013)	Quadratic+	Single / ODE state	Symplectic, unconditionally stable, mixed variable systems
Universal mesh (Gawlik et al., 2014)	Arbitrary	Reference domain mapping	High-fidelity geometry, weak reformulation, adaptive remesh
Isoparametric BDF (Lou et al., 2021)	$k$ (BDF$k$)	Level set + mapping	Projection operators, support for moving domains, robustness
PC/BDF for surface diffusion (Jiang et al., 14 Dec 2024, Gan et al., 18 Oct 2025)	$2$–$4$	PFEM, predictor-corrector	Mesh equidistribution, structure-preserving, semi-implicit
SAV/structure-preserving (Li et al., 12 Jul 2024, Garcke et al., 24 Aug 2024, Garcke et al., 17 Aug 2025)	$1$–$2$+	Lagrange multipliers, SAV	Area/energy/volume conservation, unconditional stability

7. Implications, Performance, and Future Directions

The development of temporally high-order parametric FEMs has enabled the simulation of challenging problems exhibiting rapid geometric evolution, strong nonlinearity, or multiscale temporal phenomena. Rigorous error analysis and structure-preserving discretization guarantee long-time stability, conservation, and the avoidance of nonphysical artifacts. Mesh evolution strategies, including continuous motion, discrete "rezone," universal mesh techniques, and isoparametric mappings, provide robust frameworks for complex geometries and topological changes.

Computationally, these methods support adaptive mesh refinement, matrix-free implementations, and parallel algorithms essential for exascale scalability (Kolev et al., 2021). Modular combinations of spatial and temporal approximations (tensor-product, mixed variational, and splitting schemes), high-order time-integrators (BDF$k$, Runge–Kutta–Chebyshev, predictor–corrector), and stabilized formulations (artificial diffusion, SAV, Lagrange multipliers) yield flexible toolkits for real-world applications encompassing multiphase flows, solidification, thin film dewetting, evolving biomembranes, and geometric flows.

Open challenges remain in balancing computational cost, maintaining geometric fidelity under extreme deformations, further extending structure-preservation to nonlinear and fully coupled multiphysics, and automated adaptive procedures that simultaneously optimize spatial, temporal, and geometric accuracy. The field continues to evolve rapidly, with ongoing work integrating high-order methods into large-scale, multiphysics simulation infrastructure.