[论文解读] 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 Optimization参考文献 8被引用 17

一句话总结

本文提出了一种在 MATLAB 和 Octave 中完全向量化、无循环的算法，用于高效组装有限元矩阵——特别是二维 P1 拉格朗日元和线弹性问题的质点矩阵与刚度矩阵——其性能可与专用软件（如 FreeFEM++）相媲美。该方法避免了低级编程，利用优化的稀疏矩阵处理和基于数组的索引，优于标准方法和以往的向量化方法。

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

研究动机与目标

• 解决由于基于循环的实现导致 MATLAB 和 Octave 中标准有限元矩阵组装效率低下的问题。
• 克服在新版 MATLAB 发行版（R2012a/b）中因 sparse() 函数内稀疏矩阵重新分配导致的性能下降问题。
• 开发一种无循环的向量化组装算法，利用数组运算，实现与专用 FEM 软件相当或更优的性能。
• 证明该方法在 Pk 单元和三维问题上的通用性与可扩展性，并应用于线弹性问题。
• 提供一个完全兼容 MATLAB 和 Octave 的开源实现，兼具高效性与可修改性。

提出的方法

• 用预计算的索引数组（Ig, Jg, Kg）替代标准的基于循环的组装方式，用于全局矩阵条目。
• 通过矩阵-数组运算同时计算所有单元的单元矩阵（如质量矩阵和刚度矩阵），避免显式循环。
• 通过将 sparse(nq,nq) 替换为 spalloc(nq,nq,9*nme) 优化稀疏矩阵初始化，防止在 R2012a/b 版本中出现性能退化。
• 通过单元特定函数（如 ElemStiffMatP1）应用仿射变换和数值积分，计算局部矩阵。
• 通过 sparse(Ig, Jg, Kg, nq, nq) 一次性高效组装所有单元贡献，构建全局矩阵。
• 通过避免使用语言特定或低级构造，确保与 MATLAB 和 GNU Octave 的兼容性。

实验结果

研究问题

• RQ1 在 MATLAB/Octave 中，是否可以实现完全向量化、无循环的算法，使其性能与专用 FEM 软件（如 FreeFEM++）相当？
• RQ2 为何 MATLAB R2012a/b 中标准的稀疏矩阵初始化会导致与早期版本相比显著的性能下降？
• RQ3 如何在不使用积分公式或低级语言的情况下，将向量化应用于有限元矩阵组装？
• RQ4 所提出的算法在多大程度上可推广至其他有限元类型（如 Pk、三维问题）和物理问题（如线弹性）？
• RQ5 高效的稀疏矩阵初始化在整体组装性能中起什么作用？如何对其进行优化？

主要发现

• 所提出的向量化算法（OptV2）在二维 P1 有限元上的性能与 FreeFEM++ 相当，相同硬件条件下计算时间在 FreeFEM++ 的 1.5 倍以内。
• 作者确认了 MATLAB R2012a/b 中稀疏矩阵初始化导致的性能退化问题，并已向 MathWorks 报告，确认使用 spalloc 可解决该问题。
• 与 R2011a 相比，使用 spalloc 替代 sparse 可将 1600 自由度问题的计算时间减少最多 44%，40000 自由度问题的加速比最高达 2.95 倍。
• 与经典循环方法（OptV0）及以往的向量化方法（如 [2, 3, 12, 17]）相比，该向量化方法在速度和简洁性方面均表现更优。
• 该方法可自然推广至线弹性问题，证明了其在不同类型有限元问题中的通用性与可重用性。
• 整个实现完全兼容 MATLAB 和 GNU Octave，确保了广泛可用性与可复现性。