有限元法 FEM

有限元法 (Finite Element Method, FEM) 的数学过程：

二维泊松方程

以二维泊松方程 \(-\nabla^2 u=f\) 在区域 \(\Omega\) 上为例：

1. 网格剖分（meshing）。将区域 \(\Omega\) 剖分为 \(N_e\) 个单元（如三角形）。可以用 p 次多项式估计每个三角形内部的 \(u\) 的值，比如 1 次多项式，那么设三角形顶点上的函数值分别为 \(u_{e,1},u_{e,2},u_{e,3}\)，那么就可以插值出来三角形内部的函数值（一次函数近似）\(u_e(x,y)=c_{e,1}(x,y)u_{e,1}+c_{e,2}(x,y)u_{e,2}+c_{e,3}(x,y)u_{e,3}\)，其中 \(c_n(x,y)=\sum_{n=node_{e,i}} c_{e,i}(x,y)\) 称作形函数（Shape Functions），是一个只与第 \(e\) 个单元位置相关的已知的 p 次函数. 注意未知元是顶点上的函数值，而非单元（三角形）的什么什么值。
2. 弱形式与离散化：将原泊松方程转化为弱形式（乘测试函数 \(v\) 后积分）（涉及到一点张量符号之间的变换）：\(\int_\Omega \nabla u \cdot \nabla v ~ d\Omega = \int_\Omega fv ~ d\Omega\)，然后把 \(u\) 用形函数展开：（假设一共有 \(n\) 个顶点）

\[\sum_{n=1}^N u_n \int_\Omega \nabla c_n \cdot \nabla v ~ d\Omega = \int_\Omega f v ~ d\Omega \]

分别取 \(v\) 为 \(c_k, ~ k=1,2,\dots,N\) 可以得到 N 个方程，且第 \(k\) 个方程中的 \(u_n\) 的系数只有当第 \(k,n\) 顶点同处于同一个单元上的时候系数才非零。因此实际上只有相邻点之间才有边，边集正好就是三角形剖分中的边集。

可以证明（应该吧）矩阵系数对称正定。因此离散后的能量泛函 \(J(u)=\frac{1}{2} u^T A u - f^T u\) 有唯一极小值点，可以用 CG。

三维 Stokes 方程

参考论文：Kohl, Nils, and Ulrich Rüde. "Textbook efficiency: massively parallel matrix-free multigrid for the Stokes system." SIAM Journal on Scientific Computing 44.2 (2022): C124-C155.

以三维 Stokes 方程为例。包括动量方程（粘性力 + 压力梯度 = 外力）\(-\Delta u + \nabla p = f\) 和连续性方程（不可压缩条件） \(\nabla \cdot u = 0\)。其中 \(u\) 为速度向量，\(p\) 为压力向量。

这里有一个FEM的关键问题——inf-sup条件（LBB 条件）。如果速度和压力用相同的线性多项式（P1-P1）近似，离散得到的线性系统会不稳定，压力会出现虚假震荡。因此，要么使用更精细的速度空间 + 较粗糙的压力空间，例如 P2-P1（Taylor-Hood 元），要么对 P1-P1 做稳定化处理（PSPG 方法）

例如 P2-P1 方法，对速度 \(u\) 做二次多项式近似，压力 \(p\) 做一次多项式近似，则\(u\) 每个单元（四面体）的每个（速度方向）分量需要 10 个自由度来插值（\(a + b*x + c*y + d*z + e*x^2 + f*y^2 + g*z^2 + h*x*y + i*x*z + j*y*z\)），一共 30 个自由度；压力 \(p\) 每个单元需要 4 个自由度来插值（\(a + b*x + c*y + d*z\)）。

因此，可以设速度 \(u\) 的每个方向的分量在四面体的4个顶点、6个边中点的值，和压力 \(p\) 在四面体的4个顶点的值，作为未知数，得到一共 \(4\#(vertex)+\#(edge)\) 个自由度。

HHG