有限元法之有限元空间的构造

一、区域Ω的剖分

二、三角形一次元

三、一次元的基函数与面积坐标

四、三角形二次元及其基函数

前两节我们介绍了有限元基本概念和变分理论的推导，本节我们继续探讨有限元空间的构造。

一、区域Ω的剖分

对矩形区域 $\Omega =[x_{a},x_{b}]\times [y_{c},y_{d}]$ 进行三角剖分，其中x方向剖分m份，y方向剖分n份，共得到 $(m+1)(n+1)$ 个节点及 $2mn$ 个三角形单元。图1是 $m=5,n=4$ 的剖分情况，节点编号用数字表示，单元用带圈的数字表示。为了实现后面的程序编写，必须明确单元上的局部编号与整体编号，如图2所示。通过设置剖分数，可以建立单元上整体编号与局部编号之间的关系，可设置二维数组 $lnd[\;][\;]$ ，第一个参数为单元编号，第二个参数为局部节点编号，如 $lnd[3][0]=8$ 等，表示第3个单元第0号局部节点的整体节点编号为8，而 $lnd[2][1]=2$ 则表示第2个单元第1号局部节点的整体节点编号为2。可以通过循环设置所有的节点。

二、三角形一次元

前面两节提到，可以选取 $V_{h}\subset V=H^{1}_{0}(\Omega)$ 为分片连续的一次多项式函数空间，也就是在每个单元e上， $V_{h}$ 中的函数都是一次多项式，且要保证整体连续。因此对于相邻的两个三角形单元，它们有一条公共边，只要保证分片一次多项式在这条公共边的两个端点（也是剖分节点）处函数值相同即可保证函数整体连续。这样，分片一次多项式在每个单元上的表达式就可以由它在3个顶点处的值唯一确定。下面，在节点 $P_{i},P_{j},P_{k}$ （对应整体编号为i，j，k）的单元e上考虑数值解 $u_{h}$ 的表达式，尝试用基函数来表示 $u_{h}(x,y)|_{e}=u_{i}\lambda _{0}(x,y)+u_{j}\lambda _{1}(x,y)+u_{k}\lambda _{2}(x,y)$ ，其中 $\lambda_{0},\lambda_{1},\lambda_{2}$ 为待定基函数，满足以下性质：

$\lambda_{0}(P_{i})=1,\lambda_{0}(P_{j})=0,\lambda_{0}(P_{k})=0 \;\;\;\;\; (1)$

$\lambda_{1}(P_{i})=0,\lambda_{1}(P_{j})=1,\lambda_{1}(P_{k})=0 \;\;\;\;\; (2)$

$\lambda_{2}(P_{i})=0,\lambda_{2}(P_{j})=0,\lambda_{2}(P_{k})=1 \;\;\;\;\; (3)$

且它们都是一次函数。这样，数值解 $u_{h}$ 在单元e上的表达式完全由它在3个顶点处 $P_{i},P_{j},P_{k}$ 处的值 $u_{i},u_{j},u_{k}$ 决定， $u_{i},u_{j},u_{k}$ 可以看作精确解u在整体编号i，j，k的节点处的近似。一旦把所有 $u_{i},i=0,1,\cdots,(m+1)(n+1)-1$ 求出来（边界点除外，因为 $u_{h}\in V_{h}$ 从而边界节点处 $u_{h}$ 的值为零），则数值解 $u_{h}$ 的表达式也就确定了。所以现在的基本问题是对离散问题式

求 $u_{h}(x,y)\in V_{h}$ ，使得 $a(u_{h},v_{h})=(f,v_{h})\;\;\;\;\forall v_{h}(x,y)\in V_{h}$

建立 $u_{i},i=0,1,\cdots,(m+1)(n+1)-1$ 的关系式。

三、一次元的基函数与面积坐标

由于基函数在单元e上是一次多项式，尝试设 $\lambda_{0}(x,y)|_{e}=ax+by+c$ ，其中a，b，c为待定系数，且单元e上s号节点 $P_{s}$ 的坐标为 $(x_{s},y_{s}),s=i,j,k$ ，则由条件公式（1）可知：

$\left\{\begin{matrix} ax_{i}+by_{i}+c=1,\\ ax_{j}+by_{j}+c=0,\\ ax_{k}+by_{k}+c=0, \end{matrix}\right.\;\;is$ $\begin{pmatrix} x_{i} & y_{i} & 1\\ x_{j} & y_{j} & 1\\ x_{k} & y_{k} & 1 \end{pmatrix}\begin{pmatrix} a\\ b\\ c \end{pmatrix}=\begin{pmatrix} 1\\ 0\\ 0 \end{pmatrix}$

从而解出

$a=\frac{\begin{vmatrix} 1 & y_{i} & 1\\ 0 & y_{j} & 1\\ 0 & y_{k} & 1 \end{vmatrix}}{\begin{vmatrix} x_{i} & y_{i} & 1\\ x_{j} & y_{j} & 1\\ x_{k} & y_{k} & 1 \end{vmatrix}}=\frac{y_{j}-y_{k}}{\begin{vmatrix} x_{i} & y_{i} & 1\\ x_{j} & y_{j} & 1\\ x_{k} & y_{k} & 1 \end{vmatrix}}$ $b=\frac{\begin{vmatrix} x_{i} & 1 & 1\\ x_{j} & 0 & 1\\ x_{k} & 0 & 1 \end{vmatrix}}{\begin{vmatrix} x_{i} & y_{i} & 1\\ x_{j} & y_{j} & 1\\ x_{k} & y_{k} & 1 \end{vmatrix}}=\frac{x_{k}-x_{j}}{\begin{vmatrix} x_{i} & y_{i} & 1\\ x_{j} & y_{j} & 1\\ x_{k} & y_{k} & 1 \end{vmatrix}}$

$c=\frac{\begin{vmatrix} x_{i} & y_{i} & 1\\ x_{j} & y_{j} & 0\\ x_{k} & y_{k} & 0 \end{vmatrix}}{\begin{vmatrix} x_{i} & y_{i} & 1\\ x_{j} & y_{j} & 1\\ x_{k} & y_{k} & 1 \end{vmatrix}}=\frac{x_{j}y_{k}-x_{k}y_{j}}{\begin{vmatrix} x_{i} & y_{i} & 1\\ x_{j} & y_{j} & 1\\ x_{k} & y_{k} & 1 \end{vmatrix}}$

代入可得

$\lambda_{0}(x,y)|_{e}=\frac{x(y_{j}-y_{k})+y(x_{k}-x_{j})+(x_{j}y_{k}-x_{k}y_{j})}{\begin{vmatrix} x_{i} & y_{i} & 1\\ x_{j} & y_{j} & 1\\ x_{k} & y_{k} & 1 \end{vmatrix}}=\frac{\begin{vmatrix} x & y & 1\\ x_{j} & y_{j} & 1\\ x_{k} & y_{k} & 1 \end{vmatrix}}{\begin{vmatrix} x_{i} & y_{i} & 1\\ x_{j} & y_{j} & 1\\ x_{k} &y_{k} &1 \end{vmatrix}}$

可以证明以 $P_{i},P_{j},P_{k}$ （逆时针排列）为顶点的三角形单元e的面积 $S_{e}=\frac{1}{2}\begin{vmatrix} x_{i} & y_{i} & 1\\ x_{j} & y_{j} & 1\\ x_{k} & y_{k} & 1 \end{vmatrix}$ 。

于是，若 $\Delta P_{i}P_{j}P_{k}$ 内有一点P的坐标为 $(x,y)$ ，如图3所示，则

$\lambda_{0}(x,y)|_{e}=\frac{\begin{vmatrix} x & y & 1\\ x_{j} & y_{j} & 1\\ x_{k} & y_{k} & 1 \end{vmatrix}}{\begin{vmatrix} x_{i} & y_{i} & 1\\ x_{j} & y_{j} & 1\\ x_{k} &y_{k} &1 \end{vmatrix}}=\frac{2S_{\Delta PP_{j}P_{k}}}{2S_{\Delta P_{i}P_{j}P_{k}}}=\frac{S_{\Delta PP_{j}P_{k}}}{S_{e}}\;\;\;(4)$

同理，

$\lambda_{1}(x,y)|_{e}=\frac{\begin{vmatrix} x_{i} & y_{i} & 1\\ x & y & 1\\ x_{k} & y_{k} & 1 \end{vmatrix}}{\begin{vmatrix} x_{i} & y_{i} & 1\\ x_{j} & y_{j} & 1\\ x_{k} &y_{k} &1 \end{vmatrix}}=\frac{S_{\Delta P_{i}PP_{k}}}{S_{e}},\lambda_{2}(x,y)|_{e}=\frac{\begin{vmatrix} x_{i} & y_{i} & 1\\ x_{j} & y_{j} & 1\\ x & y & 1 \end{vmatrix}}{\begin{vmatrix} x_{i} & y_{i} & 1\\ x_{j} & y_{j} & 1\\ x_{k} &y_{k} &1 \end{vmatrix}}=\frac{S_{\Delta P_{i}P_{j}P}}{S_{e}}\;\;\;(5)$

注意到 $S_{e}=S_{\Delta P_{i}P_{j}P_{k}}=S_{\Delta PP_{j}P_{k}}+S_{\Delta P_{i}PP_{k}}+S_{\Delta P_{i}P_{j}P}$ ，显然有

$\lambda_{0}+\lambda_{1}+\lambda_{2}=1\;\;\;\;(6)$

也就是说 $\lambda_{0},\lambda_{1},\lambda_{2}$ 不是相互独立的。换言之， $\Delta P_{i}P_{j}P_{k}$ 内任一点 $P(x,y)$ ，必然可以唯一对应一组坐标 $(\lambda_{0},\lambda_{1})$ ，基函数 $\lambda_{0},\lambda_{1},\lambda_{2}$ 被称为重心坐标。由于它们又都是三角形的面积比，所以它们也称为面积坐标。面积坐标在有限元分析中非常重要，它是从一般单元变化到标准单元的工具，也是进行Sobolev空间范数估计的有效手段。事实上，公式（4）、（5）可以反解出直角坐标 $(x,y)$ 与重心坐标之间的对应关系式：

$\left\{\begin{matrix} x=x_{i}\lambda_{0}+x_{j}\lambda_{1}+x_{k}\lambda_{2}\\ y=y_{i}\lambda_{0}+y_{j}\lambda_{1}+y_{k}\lambda_{2} \end{matrix}\right.\;\;\;\; or\;\;\;\left\{\begin{matrix} x=(x_{i}-x_{k})\lambda_{0}+(x_{j}-x_{k})\lambda_{1}+x_{k}\\ y=(y_{i}-y_{k})\lambda_{0}+(y_{j}-y_{k})\lambda_{1}+y_{k} \end{matrix}\right.\;\;\;(7)$

从而可以实现将一般的三角形单元 $\Delta P_{i}P_{j}P_{k}$ 变换成标准单元 $\widehat{e}$ ，如图4所示。

四、三角形二次元及其基函数

我们除了可以选取 $V_{h}$ 为分片连续的一次多项式函数空间外，也可以选取 $V_{h}$ 为分片连续的二次多项式函数空间，也就是在每个单元e上， $V_{h}$ 中的函数都是二次多项式，且要保证整体连续。因此在每个单元e上， $V_{h}$ 中的分片二次多项式函数 $v(x,y)$ 就形如 $v|_{e}=Ax^{2}+Bxy+Cy^{2}+Dx+Ey+F$ ，其中 $A,B,C,D,E,F$ 均为待定常数，从而需要有6个条件来唯一确定这个表达式。与一次元相似，要确定这6个常数，我们可以取三角形单元e的3个顶点及3条边的中点值作为条件（这些条件称为自由度），即分片二次多项式在每个单元上的表达式就可以由它在这个单元3个顶点和3条边的中点处的值唯一确定，这样也可以保证函数的整体连续性。事实上，在相邻的两个三角形单元上的公共边上，位置变量x和y有一个直线方程的线性约束，从而 $v(x,y)$ 在这条边上成为一个只关于自变量x的二次函数，这个函数在3个不同的点（两个顶点和一个中点）上取值相同，说明 $v(x,y)$ 在公共边上的表达式所示唯一确定的，也就是说，这个分片二次多项式在相邻两个单元上虽然整体表达式不相同，但在其公共边上表达式相同，这就保证了函数在 $\Omega$ 上整体连续，从而实现 $V_{h}\subset V=H^{1}_{0}(\Omega)$ 。

对于以上的三角形二次元，由于涉及到三角形单元的中点，所以尽管三角形剖分情况不变，即共有2mn个三角形单元，但整体节点数变为 $(2m+1)(2n+1)$ 个，且节点的编号将随之发生改变。例如，图1将变为图5。

接下来，在单元e上考虑数值解 $u_{h}\in V_{h}$ 的表达式，其中e的3个顶点为 $P_{i},P_{j},P_{k}$ （对应整体编号为i，j，k），3条边的中点为 $P_{jk},P_{ki},P_{ij}$ （对应整体编号为 $\frac{j+k}{2},\frac{k+i}{2},\frac{i+j}{2}$ ），如图6。

$u_{h}$ 在单元e上的表达式尝试用基函数表示为

$u_{h}(x,y)|_{e}=u_{i}\varphi_{0}(x,y)+u_{j}\varphi_{1}(x,y)+u_{k}\varphi_{2}(x,y)+u_{jk}\Psi_{0}(x,y)+u_{ki}\Psi_{1}(x,y)+u_{ij}\Psi_{2}(x,y)$

其中 $\varphi_{0},\varphi_{1},\varphi_{2},\Psi_{0},\Psi_{1},\Psi_{2}$ 为待定基函数，满足以下性质：

$\varphi_{0}(P_{i})=1,\varphi_{0}(P_{j})=0,\varphi_{0}(P_{k})=0,\varphi_{0}(P_{jk})=0,\varphi_{0}(P_{ki})=0,\varphi_{0}(P_{ij})=0,$

$\varphi_{1}(P_{i})=0,\varphi_{1}(P_{j})=1,\varphi_{1}(P_{k})=0,\varphi_{1}(P_{jk})=0,\varphi_{1}(P_{ki})=0,\varphi_{1}(P_{ij})=0,$

$\varphi_{2}(P_{i})=0,\varphi_{2}(P_{j})=0,\varphi_{2}(P_{k})=1,\varphi_{2}(P_{jk})=0,\varphi_{2}(P_{ki})=0,\varphi_{2}(P_{ij})=0,$

$\Psi_{0}(P_{i})=0,\Psi_{0}(P_{j})=0,\Psi_{0}(P_{k})=0,\Psi_{0}(P_{jk})=1,\Psi_{0}(P_{ki})=0,\Psi_{0}(P_{ij})=0,$

$\Psi_{1}(P_{i})=0,\Psi_{1}(P_{j})=0,\Psi_{1}(P_{k})=0,\Psi_{1}(P_{jk})=1,\Psi_{1}(P_{ki})=1,\Psi_{1}(P_{ij})=0,$

$\Psi_{2}(P_{i})=0,\Psi_{2}(P_{j})=0,\Psi_{2}(P_{k})=0,\Psi_{2}(P_{jk})=1,\Psi_{2}(P_{ki})=0,\Psi_{2}(P_{ij})=1.$

利用重心坐标，很容易将上述基函数表示出来，即有分别对应于三角形单元3个顶点 $P_{i},P_{j},P_{k}$ 的基函数：

$\varphi_{0}(x,y)=\lambda_{0}(2\lambda_{0}-1),\varphi_{1}(x,y)=\lambda_{1}(2\lambda_{1}-1),\varphi_{2}(x,y)=\lambda_{2}(2\lambda_{2}-1)$

及对应于三角形3条边中点 $P_{jk},P_{ki},P_{ij}$ 的基函数：

$\Psi_{0}=4\lambda_{1}\lambda_{2},\Psi_{1}=4\lambda_{2}\lambda_{0},\Psi_{2}=4\lambda_{0}\lambda_{1}$

至此，数值解 $u_{h}$ 在单元e上的表达式就确定为：

$u_{h}(x,y)|_{e}=u_{i}\lambda_{0}(2\lambda_{0}-1)+u_{j}\lambda_{1}(2\lambda_{1}-1)+u_{k}\lambda_{2}(2\lambda_{2}-1)+4u_{jk}\lambda_{1}\lambda_{2}+4u_{ki}\lambda_{2}\lambda_{0}+4u_{ij}\lambda_{0}\lambda_{1}$

综上，有限元空间 $X_{h}$ 由一个三元组 $(e,V_{h},\sum)$ 确定。具体的，设 $\tau_{h}$ 是区域Ω的一个剖分，e是剖分 $\tau_{h}$ 中的单元，参数h定义为所有单元的最大直径，即 $h=\underset{e\in\tau_{h}}{max}(diam(e))$ ， $V_{h}$ 是选定的分片多项式函数空间， $\sum$ 是每个e上用于唯一确定 $V_{h}$ 内的多项式函数所需要的条件。