[Paper Review] Programming of Finite Element Methods in MATLAB

Long Chen|arXiv (Cornell University)|Apr 14, 2018

Matrix Theory and Algorithms5 references6 citations

TL;DR

This paper presents a systematic approach to implementing the linear finite element method for solving the Poisson equation in MATLAB, emphasizing efficient sparse matrix usage and data structures for triangulations and boundary conditions. It demonstrates how to assemble stiffness and mass matrices using sparse matrix operations, handle various boundary conditions, and apply numerical quadrature rules, achieving optimal performance through high-level MATLAB coding and sparse matrix techniques.

ABSTRACT

We discuss how to implement the linear finite element method for solving the Poisson equation. We begin with the data structure to represent the triangulation and boundary conditions, introduce the sparse matrix, and then discuss the assembling process. We pay special attention to an efficient programming style using sparse matrices in MATLAB.

Motivation & Objective

• To develop a clear and efficient implementation of the linear finite element method in MATLAB for solving the Poisson equation.
• To design optimal data structures for representing triangulations and boundary conditions using global and local vertex indexing.
• To demonstrate the use of sparse matrices in MATLAB to improve computational efficiency and reduce memory usage in finite element assembly.
• To provide practical guidance on assembling stiffness and mass matrices using numerical quadrature on simplicial elements.
• To address well-posedness issues in Neumann problems, including kernel handling and compatibility conditions.

Proposed method

• Represent triangulations using two matrices: 'e' for geometric coordinates and 'm' for global vertex indices of simplices, ensuring consistent orientation (counter-clockwise in 2D, right-hand rule in 3D).
• Encode boundary conditions using a dense matrix 'lag' where values 0 (interior), 1 (Dirichlet), 2 (Neumann), 3 (Robin) label boundary faces.
• Utilize MATLAB’s built-in sparse matrix format (CSC) to store and manipulate large, sparse stiffness and mass matrices efficiently.
• Assemble global matrices by looping over elements and adding local contributions using sparse matrix addition, avoiding dense matrix operations.
• Apply numerical quadrature rules (e.g., centroid, midpoint, Gauss) in barycentric coordinates to compute element matrices accurately and efficiently.
• Handle the nullspace of the stiffness matrix in Neumann problems by enforcing a zero-mean constraint via u = u - mean(u) after solution.

Experimental results

Research questions

• RQ1 How can triangulations and boundary conditions be efficiently represented in MATLAB using sparse data structures?
• RQ2 What is the optimal way to assemble finite element matrices using MATLAB’s sparse matrix operations?
• RQ3 How do different numerical quadrature rules affect the accuracy and efficiency of finite element assembly?
• RQ4 How can the well-posedness issues of the Poisson equation with Neumann boundary conditions be addressed in the discrete setting?
• RQ5 What is the impact of using high-level versus low-level coding styles on performance and readability in finite element implementations?

Key findings

• The use of sparse matrices in MATLAB significantly reduces memory usage and accelerates matrix assembly, especially for large-scale problems.
• The proposed data structure with matrices 'e' and 'm' enables efficient local refinement and coarsening of triangulations.
• The trapezoidal rule is more efficient than the centroid rule in 2D due to fewer function evaluations (N vs. NT, with NT ≈ 2N).
• Second-order quadrature rules exist for triangles (using edge midpoints) but not for tetrahedra, requiring alternative symmetric rules.
• The discrete compatibility condition ∫f dx = ∫g dS must be enforced numerically; otherwise, the Neumann problem may be unsolvable.
• Solutions to Neumann problems are unique only up to a constant; enforcing ∫u dx = 0 ensures uniqueness and improves L² error computation.