2D Poisson 方程与 Finite Element Method

上周上 Advanced Scientific Computing 的时候老板讲了 2D Poisson 的有限元法,用的是四边形等参数单元。他布置的作业是写第一类边界条件的求解器。由于我印象中 FEM 用三角网格比较多,因此自己又对照着写了三角形等参数单元的版本,以及把第三类边界条件的情况也做了出来。下面是一点笔记,主要包括计算用到的公式和方法。

求定解问题:
$$ -\Delta u = f, \ (x,y) \in \Omega \\ \frac{\partial u }{\partial \vec{n}} + \alpha u = g, \ (x,y) \in \partial \Omega $$
这是一个第三边值问题(Robin CB,或者称为 general Neumann BC,因为 $\alpha = 0$ 时即为 Neumann BC)。记
$$ S^1 = \{ v | \iint_{\Omega} [v^2 + (\nabla^2 v)] dxdy < \infty \} $$
则上述定解问题有变分形式:求 $u \in S^1$ s.t.
$$ \iint_{\Omega} \nabla u \cdot \nabla v dx dy + \int_{\partial \Omega} \alpha u v ds = \iint_{\Omega} fv dxdy + \int_{\partial \Omega} gv ds \ \ \ \ \ (1) $$

其中 $u$ 可表示为 $u = \sum c_i \phi_i$.

等参数单元:三角形与四边形

为了处理方便,我将三角形单元和四边形单元都变换到标准单元上进行计算,并且只讨论最简单的情形:对三角形单元使用线性插值基函数,对四边形单元使用双线性插值基函数。

任意四边形变换到 $[-1,1] \times [-1, 1]$ 上计算,变量替换为:
$$ x(\xi, \eta) = \sum_{k=1}^4 x_k^e N_k(\xi, \eta), \ y(\xi, \eta) = \sum_{k=1}^4 y_k^e N_k(\xi, \eta) \\ N_1(\xi, \eta) = \frac{1}{4}(1 - \xi)(1 - \eta) \\ N_2(\xi, \eta) = \frac{1}{4}(1 + \xi)(1 - \eta) \\ N_3(\xi, \eta) = \frac{1}{4}(1 + \xi)(1 + \eta) \\ N_4(\xi, \eta) = \frac{1}{4}(1 - \xi)(1 + \eta) \\ N(\xi, \eta) = \phi(x(\xi, \eta), y(\xi, \eta)) $$
其中 $(x_k^e, y_k^e)$ 是四边形元中按逆时针顺序(下面所有的对单元结点的编号都是按逆时针顺序的)数第 $k$ 个顶点在原求解域中的坐标,$N_{k}(\xi, \eta)$ 是标准单元上的四个插值基函数,$\phi$ 是试探空间的基函数。示意图如下:

QuadBasis

任意三角形变换到 $(0,0)-(1,0)-(0,1)$ 三角形上计算,同理变量替换为:
$$ x(\xi, \eta) = \sum_{k=1}^3 x_k^e N_k(\xi, \eta), \ y(\xi, \eta) = \sum_{k=1}^3 y_k^e N_k(\xi, \eta) \\ N_1(\xi, \eta) = 1 - \xi - \eta \\ N_2(\xi, \eta) = \xi \\ N_3(\xi, \eta) = \eta \\ N(\xi, \eta) = \phi(x(\xi, \eta), y(\xi, \eta)) $$
示意图如下:

TriBasis

变量替换以后,我们还需要处理一下二重积分以及偏导数的关系:
$$ \iint_{\Omega^e} (\frac{\partial \phi_i(x,y)}{\partial x} \frac{\partial \phi_j(x,y)}{\partial x} + \frac{\partial \phi_i(x,y)}{\partial y} \frac{\partial \phi_j(x,y)}{\partial y}) dxdy \\ = \iint_{\Omega^{std}} (\frac{\partial N_a(\xi,\eta)}{\partial x} \frac{\partial N_b(\xi,\eta)}{\partial x} + \frac{\partial N_a(\xi,\eta)}{\partial y} \frac{\partial N_b(\xi,\eta)}{\partial y}) J(\xi, \eta) d \xi d \eta \ \ (2)\\ \iint_{\Omega^e} f(x, y) \phi_i(x, y) dx dy = \iint_{\Omega^{std}} f(x(\xi, \eta), y(\xi, \eta)) N_a(\xi, \eta) J(\xi, \eta) d \xi d \eta \ \ (3) \\ J(\xi, \eta) = det \left[ \begin{matrix} \frac{\partial x}{\partial \xi} & \frac{\partial x}{\partial \eta} \\ \frac{\partial y}{\partial \xi} & \frac{\partial y}{\partial \eta} \\ \end{matrix} \right] \\ \frac{\partial \xi}{\partial x} = \ \ \frac{\partial y}{\partial \eta} / J(\xi, \eta) \\ \frac{\partial \eta}{\partial x} = -\frac{\partial y}{\partial \xi} / J(\xi, \eta) \\ \frac{\partial \xi}{\partial y} = -\frac{\partial x}{\partial \eta} / J(\xi, \eta) \\ \frac{\partial \eta}{\partial y} = \ \ \frac{\partial x}{\partial \xi} / J(\xi, \eta) $$
其中 $\Omega^{std}$ 表示变换到的标准单元。

有限元方程的形成

计算单元刚度矩阵和单元荷载向量

我们考虑在每一个单元上的积分,即对于 (1) 式,$\Omega = \Omega_1 \cup \Omega_2 \cup \cdots \cup \Omega_n$ 考虑在每一个 $\Omega_e$ 上的积分。在每个单元上进行积分时,我们计算所有 $(\phi_i^{‘}, \phi_j^{‘}), (\alpha \phi_i, \phi_j), (f, \phi_j), (g, \phi_j)$ 对的积分,等号左端形成的是 $4 \times 4$ 或者 $3 \times 3$ 的矩阵,称为单元刚度矩阵;等号右端形成的是 $4 \times 1$ 或者 $3 \times 1$ 的向量,称为单元负载向量。

首先考虑 (1) 式中等号两边的第一项,都是二重积分。用高斯求积计算二重数值积分的方法如下:对于 $[-1, 1] \times [-1, 1]$ 区域,选定高斯勒让德求积节点数 $q$, 二重积分为
$$ \iint f(x,y) dx dy \approx \sum_{ix=0}^{p} \sum_{iy=0}^{p} f(gp(ix), gp(iy)) \cdot w(ix) \cdot w(iy) \ \ (4) $$

其中 $gp(), w()$ 分别是一维高斯勒让德节点和权重。

对于单位三角形,则使用更通用的形式(上式也可写成下式的形式):

$$
\iint f(x,y) dx dy \approx \sum_{iq=0}^{p} f(gpx(iq), gpy(iq)) \cdot w(iq) \ \ (5)
$$

其中 $(gpx(iq), gpy(iq))$ 是一个二维上的高斯节点,$w(iq)$ 是与之对应的权值。我使用四个高斯点:

(x, y) w
(1/3, 1/3) -27/96
(0.6, 0.2) 25/96
(0.2, 0.2) 25/96
(0.2, 0.6) 25/96

对于等号左侧第一项,给定一个高斯点 $(gpx(iq), gpy(iq))$, 我们可以计算这样一个矩阵(以四边形单元为例):
$$ S =\left[ \begin{matrix} \frac{\partial N_1}{\partial x} & \frac{\partial N_2}{\partial x} & \frac{\partial N_3}{\partial x} & \frac{\partial N_4}{\partial x} \\ \frac{\partial N_1}{\partial y} & \frac{\partial N_2}{\partial y} & \frac{\partial N_3}{\partial y} & \frac{\partial N_4}{\partial y} \end{matrix} \right] $$
与此同时,$det = J(gpx(iq), gpy(iq))$ 也可以被计算出来。随后,我们只需要计算 $det \cdot S^T S$, 便可得到一个 $4 \times 4$ 的矩阵,即为方程 (2) 在 $(gpx(iq), gpy(iq))$ 这一点所有的 $(\frac{\partial N_a}{\partial x} \frac{\partial N_b}{\partial x} + \frac{\partial N_a}{\partial y} \frac{\partial N_b}{\partial y}) J(\xi, \eta)$, $1 \le a, b \le 4$. 三角形单元同理,只是 $S$ 变成 $2 \times 3$. 将每个高斯点得到的矩阵加起来,即为所有的 $(\frac{\partial \phi_i}{\partial x} \frac{\partial \phi_j}{\partial x} + \frac{\partial \phi_i}{\partial y} \frac{\partial \phi_j}{\partial y}) $ 在对应的单元上的积分。

对于等号右侧第一项,我们只需先将 $f(x,y)$ 变换到单位单元上,再进行二重数值积分即可。

然后我们考虑 (1) 式等号左右两边的第二项。这两项都是在边界上的线积分。对单元上任意一条在 $\partial \Omega$ 上的边 $\overline{P_i P_j}$, 记其长度为 $l$, 我们在上面引入参数 $t$ 为线段弧长,对应 $P_i$ 点有 $t = 0$, $P_j$ 点有 $t = l$. 如此,单位单元上的插值基函数都可以在 $\overline{P_i P_j}$ 上化为关于 $t$ 的一次函数。因此我们有:
$$ N_i |_{\overline{P_i P_j}} = 1 - \frac{t}{l} \\ N_j |_{\overline{P_i P_j}} = \frac{t}{l} \\ N_m |_{\overline{P_i P_j}} = 0, m \neq i, m \neq j $$
对于 $\alpha(x,y), g(x,y)$,在 $\overline{P_i P_j}$ 上有
$$ \alpha(t) = \alpha(x_i + \frac{t}{l} \Delta x, y_i + \frac{t}{l} \Delta y), \Delta x = x_j - x_i, \Delta y = y_j - y_i $$
因此,对每一条在 $\partial \Omega$ 上的边,(1) 式左端第二项在给定单元上将得到一个 $4 \times 4$ 或者 $3 \times 3$ 的矩阵,其中只有4个非零元素;(2) 式右端第二项在给定单元上将得到一个 $4 \times 1$ 或者 $3 \times 1$ 的矩阵,其中只有2个非零元素。将所有在 $\partial \Omega$ 上的边的计算得到的矩阵或向量加起来,即为所有的 $(a \phi_i, \phi_j)$ 和 $(g, \phi_j)$ 在 $\partial \Omega \cap \Omega^e$ 上的线积分。

组装全局刚度矩阵和全局荷载向量

回顾 Galerkin 法,可知每一对 $(\phi^{‘}_i, \phi^{‘}_j)$, $(f, \phi_j)$, $(a \phi_i, \phi_j)$ 和 $(g, \phi_j)$ 的结果都是属于第 $j$ 条方程的;试探空间的每一个基函数 $\phi_j$ 都可以写成若干个单元上的插值基函数之和。因此,对于一个单元,若其按逆时针方向第 $i$ 个顶点在整个网格中的编号为 $gvid_i$,那么其单元刚度矩阵的第 $(i,j)$ 个元素应该被累加到全局刚度矩阵的第 $(gvid_i, gvid_j)$ 个元素上;其单元荷载向量的第 $i$ 个元素应该被累加到全局荷载向量的第 $gvid_i$ 个元素上。

计算第一类边界问题

对于第一类边界问题(Dirichlet BC),只需将 (1) 中的 $\alpha$ 和 $g$ 设为 0,即可消去 Robin BC 的限制。要添加 Dirichlet 边值的限制,即令 $c_i = u_i$。一个简单的方法是,将全局刚度矩阵的第 $i$ 行第 j 个元素设为1,第 $i$ 行其余元素设为0,并将全局荷载向量的第 $i$ 个元素设为 $u_i$。如果将边界节点集中编号,那么可以直接划去对应的行和列,并且相应地修改右端项。

数值算例

我的代码在这里

第一类边界条件

$$ \nabla^2 u(x,y) = 3x - 6y, (x, y) \in (0, 1) \times (0, 1) \\ u(x, 0) = x, u(x, 1) = 1 + x, u(0, y) = y, u(1, y) = 1 + y $$

Mathematica 11 参考结果:

MathematicaDiri

FEM 2D 四边形网格:

FEM2DQuadDiri

FEM 2D 三角形网格:

FEM2DQuadDiri

第三类边界条件

$$ \nabla^2 u(x,y) = 5xy, (x, y) \in (0, 1) \times (0, 1) \\ \frac{\partial u(x,y)}{\partial \vec{n}} + u(x,y) = x + y $$

Mathematica 11 参考结果:

MathematicaRobin

FEM 2D 四边形网格:

FEM2DQuadRobin

FEM 2D 三角形网格:

FEM2DQuadRobin

Reference

  1. 陆金甫,关治,《偏微分方程数值解法(第二版)》,清华大学出版社
  2. Finite Elements: Basis Functions
  3. Formulation of FEM for Two-Dimensional Problems