Matrix-oriented FEM formulation for stationary and time-dependent PDEs on x-normal domains

TL;DR

The paper introduces MO-FEM, a matrix-oriented finite element method that reformulates spatial discretisations of elliptic and parabolic PDEs on 2D domains into (multiterm) Sylvester equations, solvable efficiently in spectral space via MO-PCG with single-term preconditioners. It develops Kronecker-based decompositions for square and $x$-normal domains, including lumped $\mathbb{P}_1$ variants, and extends to time-dependent problems through IMEX Euler, producing a sequence of multiterm Sylvester equations. The framework delivers substantial memory and computational savings over traditional vector/Kronecker solvers while preserving or improving accuracy, demonstrated on Poisson, semilinear heat, and reaction-diffusion models, with concrete applications to Turing pattern formation in battery modeling on cap and jar geometries and curved surfaces. The results indicate MO-FEM’s potential for high-resolution spatial simulations on complex geometries, while highlighting the need for more effective preconditioners for multiterm Sylvester equations to broaden performance gains across all problem classes.

Abstract

When numerical solution of elliptic and parabolic partial differential equations is required to be highly accurate in space, the discrete problem usually takes the form of large-scale and sparse linear systems. In this work, as an alternative, for spatial discretization we provide a Matrix-Oriented formulation of the classical Finite Element Method, called MO-FEM, of arbitrary order $k\in\mathbb{N}$. On structured 2D domains (e.g. squares or rectangles) the discrete problem is then reformulated as a Sylvester matrix equation, that we solve by the reduced approach in the associated spectral space. On a quite general class of domains, namely normal domains, and even on special surfaces, the MO-FEM yields a multiterm Sylvester matrix equation where the additional terms account for the geometric contribution of the domain shape. In particular, we obtain a sequence of these equations after time discretization of parabolic problems by the IMEX Euler method. We apply the matrix-oriented form of the Preconditioned Conjugate Gradient (MO-PCG) method to solve each multiterm Sylvester equation for MO-FEM of degree $k=1,\dots,4$ and for the lumped $\mathbb{P}_1$ case. We choose a matrix-oriented preconditioner with a single-term form that captures the spectral properties of the whole multiterm Sylvester operator. For several numerical examples, we show a gain in computational time and memory occupation wrt the classical vector approach solving large sparse linear systems by a direct method or by the vector PCG with same preconditioning. As an application, we show the advantages of the MO-FEM-PCG to approximate Turing patterns with high spatial resolution in a reaction-diffusion PDE system for battery modeling.

Figures (8)

• Figure 1: Pictorial illustration of the classes of spatial domains considered in this work, together with their curvilinear meshes.
• Figure 2: An $x$-normal domain $\Omega^S$ and its mesh $\Omega_h^S$ are transformed into the curvilinear cylinder $\Gamma$ and its mesh $\Gamma_h$, respectively. The transformation $\sigma$ in \ref{['transformation_cylinder']} joins the top- and bottom edges of $\Omega_S$, highlighted in red.
• Figure 3: Poisson equation \ref{['example_laplace_square_dirichlet']} on the square $\Omega=[0,1]^2$. FEM of order $k=1,2,3,4$: comparison between the Kronecker (vector) approach \ref{['poisson_square_FEM']} (dashed lines), the MO-FEM solved by the reduced approach \ref{['neumann_sylvester_2']} (continuous lines) and, only for $k=1$, the reduced approach in closed form \ref{['dirichlet_sylvester_2']} (thick continuous line). Left plot: computational times. Right plot: convergence behavior. Dotted lines indicate slopes 2,4,5.
• Figure 4: Poisson equation \ref{['example_laplace_square_dirichlet']} on the square $\Omega=[0,1]^2$. FEM of order $k=1,2,3,4$: comparison between the "vector" formulation \ref{['poisson_square_FEM']} (dashed lines) solved by the direct method and the MO-FEM solved by the MO-PCG (continuous lines). Left plot: computational times. Right plot: Convergence behaviour of the two approaches for all $k$. Dotted lines indicate slopes 2,4 and 5. For $k=1,3,4$, the convergence is optimal of order $k+1$, for $k=2$, we observe superconvergence of order 4. PCG is cheaper in execution time and for all $k$ and for all $N$ stops after two iterations, see also Table 1.
• Figure 5: Poisson equation \ref{['example_laplace_cap-shaped-domain_dirichlet']} on the cap-shaped domain $\Omega^S$ defined in \ref{['cap-shaped-domain']}. Continuous lines indicate MO-PCG, while dashed lines indicate the vector (Kronecker) formulations solved by the direct method. In terms of computational time, MO-PCG is competitive only for $k=1$, without lumping (upper left panel). Both approaches exhibit optimal convergence for all FEM order $k$, with the case $k=2$ being superconvergent, round-off error becomes dominant on the finest mesh for $k=4$ (upper-right and lower-left panels). Number of iterations required by the MO-PCG increases with $k$ and $N$, lower-right panel.

Theorems & Definitions (6)

• Remark 1: Symmetric $x$-normal domains
• Remark 2: Curvilinear cylindrical surfaces
• Remark 3: More general domains
• Remark 4: Homogeneous Dirichlet boundary conditions
• Remark 5: Recap on lumped $\mathbb{P}_1$ matrix properties
• Remark 6: Memory performance of the reduced approach