二维问题的矩形元

一、问题与变分形式

二维 Poisson 方程（ Dirichlet 边界）：

$\begin{cases} -\Delta u = f(x,y), & (x,y) \in \Omega = [a,b] \times [c,d] \\ u|_{\partial\Omega} = g(x,y) \end{cases}$

变分形式（虚功原理）：

找 $u \in H_0^1(\Omega)$（满足边界条件的 Sobolev 空间），使

$\int_{\Omega} \nabla u \cdot \nabla v \, dxdy = \int_{\Omega} fv \, dxdy$

对任意 $v \in H_0^1(\Omega)$ 成立（$\Delta$ 是 Laplace 算子，$\nabla$ 是梯度）

二、Galerkin 离散（有限元方法）

有限元空间 $U_h$：

取 $U_h$ 为分片双线性函数空间（Lagrange 型，在矩形剖分单元上线性），满足 $u_h|_{\partial\Omega} = g_h$（$g_h$ 是 $g$ 的边界插值）

离散方程：

找 $u_h \in U_h$ 使

$\int_{\Omega} \nabla u_h \cdot \nabla v_h \, dxdy = \int_{\Omega} f v_h \, dxdy$

对任意 $v_h \in U_h$（$v_h|_{\partial\Omega} = 0$ ）成立，代入基函数展开得线性方程组

三、矩形剖分与基函数

区域剖分：

将 $\Omega$ 剖分为矩形单元 $R_{ij} = [x_i, x_{i+1}] \times [y_j, y_{j+1}]$ ，步长 $\Delta x = x_{i+1}-x_i$，$\Delta y = y_{j+1}-y_j$ 。

双线性基函数：

单元上基函数为 $1, x, y, xy$ 的线性组合（双线性插值），通过节点值 $u_{ij}, u_{i+1,j}, u_{i,j+1}, u_{i+1,j+1}$ 构造，满足分片双线性与 Kronecker delta 性质（节点处取 1，其余为 0 ）

标准化与坐标变换

对矩形单元 $R_{ij} = [x_i, x_{i+\Delta x}] \times [y_j, y_{j+\Delta y}]$，通过仿射变换标准化到参考单元 $[0,1] \times [0,1]$：

$\xi = \frac{x - x_i}{\Delta x}, \quad \eta = \frac{y - y_j}{\Delta y}$

双线性插值函数构造

设单元顶点值为 $u_{00}, u_{10}, u_{01}, u_{11}$（对应 $(\xi,\eta)=(0,0),(1,0),(0,1),(1,1)$ ），双线性插值函数 $p(\xi,\eta)$ 满足：

$p(\xi,\eta) = c_0 + c_1 \xi + c_2 \eta + c_3 \xi\eta$

通过顶点插值条件

$p(0,0)=u_{00},\quad p(1,0)=u_{10}\\ p(0,1)=u_{01},\quad p(1,1)=u_{11}$

解得 Lagrange 插值多项式：

$p(\xi, \eta) = u_{00}(1 - \xi)(1 - \eta) + u_{10} \xi(1 - \eta) + u_{01}(1 - \xi)\eta + u_{11} \xi \eta$

二维 Lagrange 基函数：

$\begin{align*} p^{(0,0)}(\xi, \eta) &= (1 - \xi)(1 - \eta) = N_0(\xi)N_0(\eta) \\ p^{(1,0)}(\xi, \eta) &= \xi(1 - \eta) = N_1(\xi)N_0(\eta) \\ p^{(0,1)}(\xi, \eta) &= (1 - \xi)\eta = N_0(\xi)N_1(\eta) \\ p^{(1,1)}(\xi, \eta) &= \xi \eta = N_1(\xi)N_1(\eta) \end{align*}$

表示为Lagrange 基函数的乘积形式：

$p(\xi,\eta) = u_{00} N_0(\xi) N_0(\eta) + u_{10} N_1(\xi) N_0(\eta) + u_{01} N_0(\xi) N_1(\eta) + u_{11} N_1(\xi) N_1(\eta)$

基函数为：

$N_0(t) = 1 - t, \quad N_1(t) = t$

四、单元插值函数（$ (x,y) \in R_{ij} $ ）

$u_h(x,y) = u_{ij} p_{ij}^{(0,0)}(\xi,\eta) + u_{i+1,j} p_{ij}^{(1,0)}(\xi,\eta) + u_{i,j+1} p_{ij}^{(0,1)}(\xi,\eta) + u_{i+1,j+1} p_{ij}^{(1,1)}(\xi,\eta)$

给矩形单元排序（从左下至上依次记为 $R_{11},R_{12},\cdots,R_{MN}$）

$u_h(x,y) = \begin{cases} u_{11} p_{11}^{(0,0)}(\xi,\eta) + u_{21} p_{11}^{(1,0)}(\xi,\eta) + u_{12} p_{11}^{(0,1)}(\xi,\eta) + u_{22} p_{11}^{(1,1)}(\xi,\eta), \quad R_{11}\\ u_{12} p_{12}^{(0,0)}(\xi,\eta) + u_{22} p_{12}^{(1,0)}(\xi,\eta) + u_{13} p_{12}^{(0,1)}(\xi,\eta) + u_{23} p_{11}^{(1,1)}(\xi,\eta), \quad R_{12} \\ \quad \quad \vdots \\ u_{1,N} p_{1N}^{(0,0)}(\xi,\eta) + u_{2,N} p_{1N}^{(1,0)}(\xi,\eta) + u_{1,N+1} p_{1N}^{(0,1)}(\xi,\eta) + u_{2,N+1} p_{1N}^{(1,1)}(\xi,\eta) , \quad R_{1N}\\ \quad \quad \vdots \\ u_{M,N} p_{MN}^{(0,0)}(\xi,\eta) + u_{M+1,N} p_{MN}^{(1,0)}(\xi,\eta) + u_{M,N+1} p_{MN}^{(0,1)}(\xi,\eta) + u_{M+1,N+1} p_{MN}^{(1,1)}(\xi,\eta), \quad R_{MN} \end{cases}$

五、整体双线性基函数

$\varphi_{ij}(x,y) = \begin{cases} \left(1 - \frac{|x - x_i|}{\Delta x}\right)\left(1 - \frac{|y - y_j|}{\Delta y}\right), & (x,y) \in R_{ij} \cup R_{i+1,j} \cup R_{i,j+1} \cup R_{i+1,j+1} \\ 0, & \text{其他} \end{cases}$

六、自由度与空间维数

若把 $u_h$ 在所有节点的函数值都当作自变量，节点数 $ (M+1) \times (N+1) $，故：

$\text{$\varphi_{i,j}$ 的个数为} (M+1)(N+1), \quad \dim(U_h) = (M+1)(N+1)$

七、网格剖分与数据结构

一、网格剖分与标号规则

节点标号：

矩形区域剖分为 $M \times N$ 单元，节点总数 $K = (M+1)(N+1)$，按顺序标号 $1, 2, \dots, K$ 。

单元标号：
单元总数 $M \cdot N$，按顺序标号 $1, 2, \dots, M \cdot N$ ，每个单元对应 4 个顶点（左下→右下→右上→左上）。

二、数据结构定义

节点数组（node）：

存储节点坐标，维度 $K \times 2$ ，每行对应一个节点的 $(x,y)$ 坐标：

$\text{node} = \begin{bmatrix} x^{(1)} & y^{(1)} \\ \vdots & \vdots \\ x^{(K)} & y^{(K)} \end{bmatrix}$

单元数组（elem）：
存储单元顶点的节点标号，维度 $M \cdot N \times 4$ ，每行对应一个单元的 4 个顶点标号（按左下→右下→右上→左上顺序）：

$\text{elem} = \begin{bmatrix} i_{j,1} & i_{j,2} & i_{j,3} & i_{j,4} \\ \vdots & \vdots & \vdots & \vdots \end{bmatrix}$ $R_{ij}(i,j)=[u_{ij},u_{i+1,j},u_{i,j+1},u_{i+1,j+1}]$

八、单元刚度矩阵积分

$\iint_{R_{ij}} \nabla u_h \cdot \nabla v_h \, dxdy = \iint_{R_{ij}} \left( \frac{\partial u_h}{\partial x} \frac{\partial v_h}{\partial x} + \frac{\partial u_h}{\partial y} \frac{\partial v_h}{\partial y} \right) dxdy$

九、插值函数与导数（局部坐标 $ \xi, \eta $ ）

单元内 $ u_h(x,y) $ 展开：

$u_h(x,y) = u_{ij} N_0(\xi)N_0(\eta) + u_{i+1,j} N_1(\xi)N_0(\eta) + u_{i,j+1} N_0(\xi)N_1(\eta) + u_{i+1,j+1} N_1(\xi)N_1(\eta)$

$ x $ 方向导数（$ \displaystyle \frac{\partial \xi}{\partial x} = \frac{1}{\Delta x} $ ）：

$\frac{\partial u_h}{\partial x} = \frac{1}{\Delta x} \left[ u_{ij} (- (1 - \eta)) + u_{i+1,j} (1 - \eta) + u_{i,j+1} (- \eta) + u_{i+1,j+1} \eta \right]$

$v_h$ 同理

十、单元刚度矩阵计算（$ x $ 方向示例）

$\begin{align} \iint_{R_{ij}} \frac{\partial u_h}{\partial x} \frac{\partial v_h}{\partial x}\,dxdy&= \int_{R_{ij}} [u_{ij}[-(1-\eta)]+\cdots]\cdot[v_{ij}[-(1-\eta)]+\cdots]\cdot\frac{1}{\Delta x_i^2}\,dxdy\\ &=\int_0^1 \int_0^1 [u_{ij}[-(1-\eta)]+\cdots]\cdot[v_{ij}[-(1-\eta)]+\cdots]\cdot\frac{1}{\Delta x_i^2}\cdot \Delta x\cdot\Delta y\,d\xi d\eta\\ &=(v_{ij},v_{i+1,j}, v_{i,j+1},v_{i+1,j+1}) \begin{bmatrix} \displaystyle \int_0^1 (1-\eta)^2 d\eta &\displaystyle \int_0^1 -(1-\eta)^2 d\eta & \times &\times \\ \times & \times & \times &\times \\ \times & \times & \times &\times \\ \times & \times & \times &\times \end{bmatrix} \frac{\Delta y}{\Delta x} \begin{bmatrix} u_{i j} \\ u_{i+1,j} \\ u_{i,j+1}\\u_{i+1,j+1} \end{bmatrix}\\ &=(v_{ij},\cdots,v_{i+1,j+1})K_ex\begin{bmatrix} u_{i j} \\ \vdots\\u_{i+1,j+1} \end{bmatrix} \end{align}$

$K_ey$ 同理

十一、整体单元刚度矩阵

$K_e = K_{e}x + K_{e}y$

$ K_{e}y $ 为 $ y $ 方向单元刚度矩阵，对称推导

十二、总刚度矩阵 $K$ 的组装

初始化：

总刚度矩阵 $K$ 维度为 $(M+1)(N+1) \times (M+1)(N+1)$，初始化为零矩阵：

$K = \text{zeros}(((M+1)(N+1), (M+1)(N+1)))$

单元循环组装：

遍历每个单元（标号 $k=1$ 到 $ne = M \cdot N$ ），计算单元刚度矩阵并叠加到总矩阵对应位置（通过节点标号索引）

伪代码示意：

$\text{for } k = 1 : n_e \\ \quad \text{将单元刚度矩阵按节点索引叠加到 } K \\ \text{end}$

十三、载荷向量的积分计算

对载荷项 $\displaystyle \int_{\Omega} f v_h \, dxdy$，在矩形单元 $R_{ij}$ 上通过双线性基函数展开，利用变量替换（标准化到参考单元 $[0,1] \times [0,1]$ ），将积分转化为：

$\begin{align} \int_{R_{ij}} f v_h \, dxdy &=\int_0^1\int_0^1f(x_i+\Delta x\cdot\xi,y_j+\Delta y\cdot\eta)[v_{ij}N_0(\xi)N_0(\eta)+\cdots+v_{i+1,j+1}N_1(\xi)N_1(\eta)]\cdot \Delta x\cdot\Delta y\,d\xi d\eta\\ &= (v_{ij}, \dots, v_{i+1,j+1}) \begin{bmatrix} \displaystyle \int_0^1 \int_0^1 f N_0(\xi)N_0(\eta) d\xi d\eta \\ \vdots \\ \displaystyle\int_0^1 \int_0^1 f N_1(\xi)N_1(\eta) d\xi d\eta \end{bmatrix} \Delta x \Delta y\\ &\approx(v_{ij}, \dots, v_{i+1,j+1})\begin{bmatrix} \frac{1}{4}f_{ij} \\ \vdots \\\frac{1}{4}f_{i+1,j+1}\end{bmatrix}\Delta x \Delta y \end{align}$

$\Delta x, \Delta y$ 是单元步长，$v_j$ 是节点载荷系数

十四、最终变分形式

离散后的变分方程表示为：

$\mathbf{V}^T A \mathbf{U} = \mathbf{V}^T \mathbf{F}\\ A=K$

其中，$A = K$ 是总刚度矩阵，$\mathbf{U}$ 是节点值向量，$\mathbf{V}$ 是检验函数向量，$\mathbf{F}$ 是载荷向量，体现有限元离散的线性方程组形式（$K \mathbf{U} = \mathbf{F}$ ）

十五、边界条件处理

节点分类：

边界节点 $ \text{bdnode} $（已知 $ U(\text{bdnode}) $ ）、自由节点 $ \text{freenode} $（待求 $ U(\text{freenode}) $ ）

变分约束：

检验函数 $ V(\text{bdnode}) = 0 $，自由节点 $ V(\text{freenode}) $ 任意，由此从整体系统 $ A U = F $ 中提取自由节点子系统：

$A(\text{freenode}, :) U - F(\text{freenode}) = 0$

子系统分解：
展开为自由节点与边界节点的矩阵块形式：

$A(\text{freenode}, \text{freenode}) U(\text{freenode}) + A(\text{freenode}, \text{bdnode}) U(\text{bdnode}) = F(\text{freenode})$

求解自由节点：
移项得自由节点的代数方程：

$A(\text{freenode}, \text{freenode}) U(\text{freenode}) = F(\text{freenode}) - A(\text{freenode}, \text{bdnode}) U(\text{bdnode})$

最终解为：

$U(\text{freenode}) = A(\text{freenode}, \text{freenode})^{-1} \left( F(\text{freenode}) - A(\text{freenode}, \text{bdnode}) U(\text{bdnode}) \right)$

核心是有限元中Dirichlet 边界条件的离散处理方法，通过节点分类和矩阵分块，将已知边界值代入整体系统，求解自由节点未知量

三角元法

一、问题与变分形式

二维 Poisson 方程（ Dirichlet 边界）：

$\begin{cases} -\Delta u = f(x,y), & (x,y) \in \Omega \\ u|_{\partial\Omega} = g(x,y) \end{cases}$

变分形式（虚功原理）：

找 $u \in H_0^1(\Omega)$（满足边界条件的 Sobolev 空间），使

$\int_{\Omega} \nabla u \cdot \nabla v \, dxdy = \int_{\Omega} fv \, dxdy$

对任意 $v \in H_0^1(\Omega)$ 成立（$\Delta$ 是 Laplace 算子，$\nabla$ 是梯度）

二、三角元法的核心优势

三角形剖分可逼近任意复杂曲线边界（通过多边形域近似 $\Omega$ ），适配工程中不规则区域（如飞机翼型、地质构造）的数值求解

三、三角剖分规则

将多边形域 $\Omega$ 剖分为有限个三角形单元 $\Delta_j$，需满足：

单元不重叠（$\Delta_i^o \cap \Delta_j^o = \emptyset, i \neq j$ ，$^o$ 表示内部）

无悬挂节点（三角形顶点不落在另一三角形的边上）

进一步定义正则剖分（最小角有下界）和拟一致剖分（单元尺寸均匀），保证有限元解的收敛性

四、Lagrange 型三角元空间

有限元空间 $U_h$ 定义为：

$U_h = \left\{ u \in H_0^1(\Omega) \mid u|_{\Delta_j} \in P^m(\Delta_j) \right\}$

其中，$P^m(\Delta_j)$ 是三角形单元上的 $m$ 次 Lagrange 多项式，形式为：

$P_m(x,y) = \sum_{i+j \leq m} c_{ij} x^i y^j$

$m$ 次多项式的系数个数为 $\frac{1}{2}(m+1)(m+2)$，需取对应数量的插值节点保证唯一性

五、单元与插值函数推导

对二维 Poisson 方程的三角剖分单元 $\Delta_{(123)}$（节点按逆时针排序为 $(x_1,y_1),(x_2,y_2),(x_3,y_3)$ ），构造线性插值函数 $p(x,y) = ax + by + c$ ，满足顶点插值条件 $p(x_i,y_i) = u_i$（$i=1,2,3$ ）

通过顶点条件建立线性方程组：

$\begin{cases} a x_1 + b y_1 + c = u_1 \\ a x_2 + b y_2 + c = u_2 \\ a x_3 + b y_3 + c = u_3 \end{cases}$

写成矩阵形式 $A \begin{bmatrix} a \ b \ c \end{bmatrix} = \begin{bmatrix} u_1 \ u_2 \ u_3 \end{bmatrix}$ ，其中 $A = \begin{bmatrix} x_1 & y_1 & 1 \ x_2 & y_2 & 1 \ x_3 & y_3 & 1 \end{bmatrix}$

解得系数：

$\begin{bmatrix} a \\ b \\ c \end{bmatrix} = A^{-1} \begin{bmatrix} u_1 \\ u_2 \\ u_3 \end{bmatrix}$

因此，插值函数可表示为：

$p(x,y) = (x, y, 1) A^{-1} \begin{bmatrix} u_1 \\ u_2 \\ u_3 \end{bmatrix}$

六、面积坐标（Lagrange 基函数）

定义面积坐标 $L_1,L_2,L_3$：

$(L_1, L_2, L_3) = (x, y, 1) A^{-1}$

则插值函数简化为：

$p(x,y) = u_1 L_1(x,y) + u_2 L_2(x,y) + u_3 L_3(x,y)$

其中，$L_i(x,y)$ 是线性 Lagrange 基函数，满足 $L_i(x_j,y_j) = \delta_{ij}$（Kronecker delta ），且 $L_1 + L_2 + L_3 = 1$

$\begin{align} L_1&=\frac{1}{\det(A^T)}\begin{vmatrix} x &x_2 &x_3 \\ y & y_2 & y_3 \\ 1 & 1 & 1 \end{vmatrix}\\ L_2&=\frac{1}{\det(A^T)}\begin{vmatrix} x_1 &x &x_3 \\ y_1 & y & y_3 \\ 1 & 1 & 1 \end{vmatrix}\\ L_3&=\frac{1}{\det(A^T)}\begin{vmatrix} x_1 &x_2 &x \\ y_1 & y_2 & y \\ 1 & 1 & 1 \end{vmatrix} \end{align}$

七、几何意义（Cramer 法则）

通过 Cramer 法则计算 $L_i$，发现分母 $ \det(A^T) = \det(A) $ 与三角形面积相关：

$\det(A) = 2 S_{\Delta_{(123)}}$

其中

$S_{\Delta_{(123)}} = \frac{1}{2} \left| \det \begin{bmatrix} x_1 & y_1 & 1 \\ x_2 & y_2 & 1 \\ x_3 & y_3 & 1 \end{bmatrix} \right|$

记 $(x_1,y_1)\rightarrow(x_2,y_2)$ 的向量为 $\vec{a}$，$(x_1,y_1)\rightarrow(x_3,y_3)$ 的向量为 $\vec{b}$，则

$\vec a \times \vec b=\begin{vmatrix} i &j & k \\ x_2-x_1 & y_2-y_1 & 0 \\ x_3-x_1 & y_3-y_1 & 0 \end{vmatrix}$ $S_{\Delta_{(123)}} = \frac{1}{2} \begin{vmatrix}\vec a \times \vec b\end{vmatrix}$

八、面积坐标与直角坐标的关系

$\begin{align} L_1&=\frac{2S_1}{2S}=\frac{S_1}{S}=\frac{(y_2 - y_3)(x - x_3) + (x_3 - x_2)(y - y_3)}{(y_2 - y_3)(x_1 - x_3) + (x_3 - x_2)(y_1 - y_3)}\\ L_2&=\frac{2S_2}{2S}=\frac{S_2}{S}=\frac{(y_3 - y_1)(x - x_3) + (x_1 - x_3)(y - y_3)}{(y_2 - y_3)(x_1 - x_3) + (x_3 - x_2)(y_1 - y_3)}\\ L_3&=\frac{2S_3}{2S}=\frac{S_3}{S}=\frac{(y_1 - y_2)(x - x_2) + (x_2 - x_1)(y - y_2)}{(y_2 - y_3)(x_1 - x_3) + (x_3 - x_2)(y_1 - y_3)}\\ \end{align}$ $x = L_1 x_1 + L_2 x_2 + L_3 x_3\\ y = L_1 y_1 + L_2 y_2 + L_3 y_3$

九、单元插值函数

三角形单元上的线性插值函数表示为：

$u_h(x,y) = u_1 L_1(x,y) + u_2 L_2(x,y) + u_3 L_3(x,y)$

$u_i$ 是顶点 $i$ 的函数值，$L_i$ 是面积坐标基函数

十、全域基函数与自由度

通过 “单元拼接” 构造全域基函数 $\varphi_j(x,y)$（在节点 $P_j$ 处取 1，其他单元为 0 ），有限元空间 $U_h$ 的维度为节点总数 $nd$（自由度）

$u_h(x,y)=\sum_{j=1}^{nd}u_j\cdot \varphi_j(x,y)$

十一、单元刚度矩阵

单元刚度矩阵 $K_e$ 由变分形式的积分推导：

$\begin{align} \int_{\Delta_j} \nabla u_h \cdot \nabla v_h \, dxdy &= \int_{\Delta_j}(u_1^{(j)}\nabla L_1^{(j)}+u_2^{(j)}\nabla L_2^{(j)}+u_3^{(j)}\nabla L_3^{(j)})(v_1^{(j)}\nabla L_1^{(j)}+v_2^{(j)}\nabla L_2^{(j)}+v_3^{(j)}\nabla L_3^{(j)})\\ &=(v_1^{(j)},v_2^{(j)}, v_3^{(j)}) \begin{bmatrix} \displaystyle \int_{\Delta j} \nabla L_1^{(j)}\nabla L_1^{(j)} dxdy &\displaystyle \int_{\Delta j} \nabla L_1^{(j)}\nabla L_2^{(j)} dxdy & \times &\\\ \times & \times & \times \\ \times & \times & \times \\ \end{bmatrix} \begin{bmatrix} u_1^j \\ u_2^j \\u_3^j \end{bmatrix}\\ &=(v_1^{(j)},v_2^{(j)}, v_3^{(j)}) \cdot K_e \cdot \begin{bmatrix} u_1^j \\ u_2^j \\u_3^j \end{bmatrix} \end{align}$

十二、面积坐标的梯度

对面积坐标 $L_i(x,y)$，通过对坐标变换矩阵求导，得梯度分量：

$\frac{\partial L_i}{\partial x} = \frac{y_j - y_k}{2S}, \quad \frac{\partial L_i}{\partial y} = \frac{x_k - x_j}{2S}$

$(i,j,k)$ 是三角形顶点的循环标号，$S$ 是三角形面积

性质：$\nabla L_i$ 是常向量（与位置无关）

十三、刚度矩阵的积分

单元刚度矩阵的积分项（梯度内积）可简化为：

$\int_{\Delta_j} \nabla L_k \cdot \nabla L_l \, dxdy = \nabla L_k \cdot \nabla L_l \cdot S$

利用 $\nabla L_i$ 是常向量，积分转化为向量点乘与面积乘积

十四、单元载荷的积分

载荷项 $\int_{\Delta_j} f u_h \, dxdy$ 分解为基函数的线性组合，积分后表示为：

$\begin{align} \int_{\Delta_j} f u_h \, dxdy &= \int_{\Delta_j} f (v_1^{(j)}L_1^{(j)}+v_2^{(j)}L_2^{(j)}+v_3^{(j)}L_3^{(j)})\, dxdy \\ &=(v_1^{(j)},v_2^{(j)},v_3^{(j)})\begin{bmatrix}\displaystyle \int_{\Delta_j} f\cdot L_1^{(j)}\,dxdy \\ \displaystyle \int_{\Delta_j} f\cdot L_2^{(j)}\,dxdy \\ \displaystyle \int_{\Delta_j} f\cdot L_3^{(j)}\,dxdy \end{bmatrix} \end{align}$