[Paper Review] An efficient way to perform the assembly of finite element matrices in Matlab and Octave

François Cuvelier, Caroline Japhet+1 more|arXiv (Cornell University)|May 14, 2013

Manufacturing Process and Optimization8 references17 citations

TL;DR

This paper presents a fully vectorized, loop-free algorithm in MATLAB and Octave for efficiently assembling finite element matrices—specifically mass and stiffness matrices for P1 Lagrange elements in 2D and linear elasticity—achieving performance comparable to dedicated software like FreeFEM++. The method avoids low-level coding and leverages optimized sparse matrix handling and array-based indexing to outperform standard and prior vectorized approaches.

ABSTRACT

We describe different optimization techniques to perform the assembly of finite element matrices in Matlab and Octave, from the standard approach to recent vectorized ones, without any low level language used. We finally obtain a simple and efficient vectorized algorithm able to compete in performance with dedicated software such as FreeFEM++. The principle of this assembly algorithm is general, we present it for different matrices in the P1 finite elements case and in linear elasticity. We present numerical results which illustrate the computational costs of the different approaches

Motivation & Objective

• Address the inefficiency of standard finite element matrix assembly in MATLAB and Octave due to loop-based implementations.
• Overcome performance degradation in newer MATLAB releases (R2012a/b) caused by sparse matrix reallocation in the sparse() function.
• Develop a vectorized assembly algorithm that eliminates loops and uses array operations to match or exceed performance of specialized FEM software.
• Demonstrate the method's generality and extendability to Pk elements and 3D problems, with application to linear elasticity.
• Provide a fully compatible, open-source implementation for MATLAB and Octave that is both efficient and easy to adapt.

Proposed method

• Replace standard loop-based assembly with a fully vectorized approach using precomputed index arrays (Ig, Jg, Kg) for global matrix entries.
• Use matrix-array operations to compute element matrices (e.g., mass and stiffness) for all elements simultaneously, avoiding explicit loops.
• Optimize sparse matrix initialization by replacing sparse(nq,nq) with spalloc(nq,nq,9*nme) to prevent performance regression in R2012a/b.
• Apply affine transformation and numerical integration via element-specific functions (e.g., ElemStiffMatP1) to compute local matrices.
• Construct global matrices using sparse(Ig, Jg, Kg, nq, nq) to efficiently assemble all element contributions in one call.
• Ensure compatibility with both MATLAB and GNU Octave by avoiding language-specific or low-level constructs.

Experimental results

Research questions

• RQ1 Can a fully vectorized, loop-free algorithm in MATLAB/Octave achieve performance competitive with dedicated FEM software like FreeFEM++?
• RQ2 Why does the standard sparse matrix initialization in MATLAB R2012a/b cause significant performance degradation compared to earlier versions?
• RQ3 How can vectorization be applied to finite element matrix assembly without using quadrature formulas or low-level languages?
• RQ4 To what extent can the proposed method be generalized to other finite element types (e.g., Pk, 3D) and physical problems (e.g., linear elasticity)?
• RQ5 What role does efficient sparse matrix initialization play in overall assembly performance, and how can it be optimized?

Key findings

• The proposed vectorized algorithm (OptV2) achieves performance comparable to FreeFEM++ for P1 finite elements in 2D, with computation times within 1.5× of FreeFEM++ on the same hardware.
• The performance regression in MATLAB R2012a/b for sparse matrix initialization was confirmed by the authors and reported to MathWorks, with spalloc identified as the fix.
• Using spalloc instead of sparse reduced computation time by up to 44% for 1600 degrees of freedom and up to 2.95× speedup for 40000 degrees of freedom compared to R2011a.
• The vectorized approach (OptV2) outperforms the classical loop-based method (OptV0) and prior vectorized methods (e.g., [2, 3, 12, 17]) in terms of speed and simplicity.
• The method extends naturally to linear elasticity, demonstrating its generality and reusability across different finite element problems.
• The entire implementation is fully compatible with both MATLAB and GNU Octave, enabling broad accessibility and reproducibility.