[论文解读] Programming of Finite Element Methods in MATLAB

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

Matrix Theory and Algorithms参考文献 5被引用 6

一句话总结

本文提出了一种系统性方法，用于在 MATLAB 中实现线性有限元法求解泊松方程，重点在于高效使用稀疏矩阵及三角剖分与边界条件的数据结构。它展示了如何通过稀疏矩阵运算组装刚度矩阵和质量矩阵，处理各种边界条件，并应用数值积分规则，通过高级 MATLAB 编程和稀疏矩阵技术实现最优性能。

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.

研究动机与目标

• 开发一种清晰且高效的线性有限元法在 MATLAB 中的实现，用于求解泊松方程。
• 设计最优的数据结构，使用全局和局部节点索引表示三角剖分和边界条件。
• 展示 MATLAB 中稀疏矩阵的使用，以提高有限元组装过程中的计算效率并减少内存占用。
• 提供实用指导，说明如何在单纯形单元上使用数值积分规则组装刚度矩阵和质量矩阵。
• 解决诺伊曼问题中的适定性问题，包括零空间处理和相容性条件。

提出的方法

• 使用两个矩阵表示三角剖分：'e' 用于几何坐标，'m' 用于单纯形的全局节点索引，确保方向一致（2D 中为逆时针方向，3D 中为右手定则）。
• 使用稠密矩阵 'lag' 编码边界条件，其中值 0（内部）、1（狄利克雷）、2（诺伊曼）、3（罗宾）标记边界面。
• 利用 MATLAB 内置的稀疏矩阵格式（CSC）高效存储和操作大规模稀疏的刚度矩阵和质量矩阵。
• 通过遍历单元并使用稀疏矩阵加法累加局部贡献来组装全局矩阵，避免使用稠密矩阵运算。
• 在重心坐标系中应用数值积分规则（例如，形心、中点、高斯积分），以准确且高效地计算单元矩阵。
• 在诺伊曼问题中通过解后执行 u = u - mean(u) 的方式处理刚度矩阵的零空间。

实验结果

研究问题

• RQ1 如何在 MATLAB 中使用稀疏数据结构高效表示三角剖分和边界条件？
• RQ2 使用 MATLAB 的稀疏矩阵运算组装有限元矩阵的最优方法是什么？
• RQ3 不同的数值积分规则如何影响有限元组装的精度和效率？
• RQ4 在离散设置中，如何解决带有诺伊曼边界条件的泊松方程的适定性问题？
• RQ5 在有限元实现中，使用高级编程风格与低级编程风格对性能和可读性有何影响？

主要发现

• 在 MATLAB 中使用稀疏矩阵可显著减少内存使用并加快矩阵组装，尤其在大规模问题中效果明显。
• 所提出的使用矩阵 'e' 和 'm' 的数据结构可高效支持三角剖分的局部加密与粗化。
• 在 2D 中，梯形积分规则比形心规则更高效，因为函数求值次数更少（N 对比 NT，其中 NT ≈ 2N）。
• 存在用于三角形的二阶积分规则（使用边中点），但不存在用于四面体的二阶规则，需采用其他对称规则。
• 必须在数值上强制执行离散相容性条件 ∫f dx = ∫g dS；否则诺伊曼问题可能无解。
• 诺伊曼问题的解仅在常数范围内唯一；通过强制 ∫u dx = 0 可确保唯一性并改善 L² 误差计算。