摘译-使用有限元方法求解泊松方程（一）
本内容翻译自《Plasma Simulations by Example》(Lubos Brieda，2019），用于在三维空间内求解泊松方程。
1 强形式
泊松方程（稳态热传导方程的形式）可以写为 \[\begin{equation} \nabla^2u\equiv\frac{\partial^2u}{\partial x^2}+\frac{\partial^2u}{\partial y^2}+\frac{\partial^2u}{\partial z^2}\equiv u_{,ii}=f \end{equation}\]
在强形式中，我们通常将方程和边界条件写成积分形式，其中下标跟在逗号后面的表示导数，重复的下标意味着求和。 这种记号方法被称为指标记号或爱因斯坦标记法。我们感兴趣的是在某个空间域\(\Omega\)上求解这个方程。 为了获得唯一解，我们还需要指定边界条件。我们假设我们知道在边界的一个子集\(\Gamma_g\)上的解， 即满足狄利克雷边界条件。我们还假设在另一个子集\(\Gamma_h\)上我们知道解的导数， 即满足纽曼边界条件。因此，按照爱因斯坦标记法，上述强形式表述为：
给定\(f\)：\(\Omega\rightarrow \R\)，\(g:\Gamma_g\rightarrow \R\)，且\(h:\Gamma_h\rightarrow \R\)，寻找\(u:\Omega\rightarrow \R\)，使得\(-ku_{,ii}=f(\Omega)\),\(u=g(\Gamma_g)\)
针对泊松方程求解问题，\(u=\phi，k=1\)且\(f=\rho/\epsilon_0\)。需要注意的是，\(\Gamma_g\notin\emptyset\)，但是\(\Gamma_h\)可以是一个空集。狄利克雷边界条件是获得唯一解所必须的边界条件，但是纽曼边界条件可选。两类边界条件不需要连续，我们可以拥有多个非连续的\(\Gamma_g\)和\(\Gamma_h\)区域。
2 弱形式
弱形式使用两种函数定义解：
- 试解
- 权重函数（变量）
解需要满足边界条件，即在\(\Gamma_g\)上\(u=g\)。同时要求解的导数平方可积：
\[\begin{equation} \int_\Omega(u_{,i})^2d\Omega<\infty \end{equation}\]
满足上述条件的函数称作\(H^1\)函数，因此\(u\in H^1\)。试解的集合由下式给出： \[\begin{equation} S=\{u|u\in H^1,u=g | \Gamma_g\} \end{equation}\] 权重函数与解相似，不同的是他们只满足\(g\)边界条件的齐次形式： \[\begin{equation} V=\{w|w\in H^1,w=0|\Gamma_g\} \end{equation}\] 于是弱形式由下述形式给出：
给定\(f:\Omega\rightarrow \R, g:\Gamma_g\rightarrow \R, h:\Gamma_h\rightarrow\R, 找到u\in S,使得对所有的w\in V,满足\int_\Omega w_{,i}(ku_{,i})d\Omega=\int_\Omega wfd\Omega+\int_{\Omega_h}whd\Gamma\)
与强形式不同，弱形式通过\(u\)的导数来定义解。这两个解是相同的，如果\(u\)是\((S)\)的解，那么它也是\((W)\)的解，如果\(u\)是\((W)\)的解，那么它也是\((S)\)的解。我们现在定义几个新的算子以简化记号: \[\begin{equation}a(w,u)=\int_\Omega w_{,i}ku_{,i}d\Omega\end{equation}\] \[\begin{equation}(w,f)=\int_\Omega wfd\Omega\end{equation}\] \[\begin{equation}(w,f)_\Gamma=\int_{\Gamma_h}wfd\Gamma\end{equation}\] 式（6）和式（7）相似，只是最后是对立体表面进行积分。
3 伽辽金形式
接下来在离散化域\(\Omega^h\)上近似弱形式。第一步是构建\(S\)和\(V\)的有限维近似。代替连续函数，我们现在有了子集\(S^h\)和\(V^h\)，其中\(S^h\in S\)，\(V^h\in V\)。下标\(h\)表示与网格相关联的属性。因此，如果\(u^h\in S^h\)，则\(u^h\in S\)。边界条件也成立，即\(u^h = g\)在\(\Gamma_g\)上成立。同样，\(V^h\)的所有成员在\(\Gamma_g\)上消失。然后，我们可以让\(S^h\)的每个成员由下式表述： \[\begin{equation}u^h=v^h+g^h\end{equation}\] 式中\(v^h\in V^h\)，\(g^h\)为在\(\Gamma_g\)上满足\(u=g\)的集合，于是有了伽辽金形式的描述。
给定\(f\)，\(g\)和\(h\)，寻找\(u^h=v^h+g^h\in V^h\)，使得对于所有的\(w^h\in V^h\)，\(a(w^h,v^h)=(w^h,f)+(w^h,h)_\Gamma-a(w^h,g^h)\)
4 形状函数
假设我们的域由\(\eta = \{0, 1, …, n_{np−1}\}\)个节点组成。在这些节点中，狄利克雷\(g\)边界条件被指定的节点由\(\eta_g\)给出。那么，\(\eta - \eta_g\)就是解\(u^h\)需要被确定的节点子集。接下来，我们让\(V^h\)的成员采取以下形式： \[\begin{equation}w^h(\vec{x})=\sum_{A\in\eta-\eta_g}N_A(\vec{x})c_A\end{equation}\] \[\begin{equation}v^h(\vec{x})=\sum_{A\in\eta-\eta_g}N_A(\vec{x})d_A\end{equation}\] 式中\(d_A\)为要求解的未知变量。\(g\)的节点值使用形函数插值得到： \[\begin{equation}g^h(\vec{x})=\sum_{A\in\eta_g}N_A(\vec{x})g_A,g_A=g(\vec{x_A})\end{equation}\] 将这些插值函数代入到伽辽金表达式中，得到 \[\begin{equation}\sum_{B\in\eta-\eta_g}a(N_A,N_B)d_B=(N_A,f)+(N_A,h)_\Gamma-\sum_{B\in\eta_g}a(N_A,N_B)g_B,A\in\eta-\eta_g,A\in\eta-\eta_g\end{equation}\]
5 矩阵形式
定义 \[\begin{equation}K_{AB}=a(N_A,N_B)\end{equation}\] \[\begin{equation}F_A=(N_A,f)+(N_A,h)_\Gamma-\sum_{B\in\eta_g}a(N_A,N_B)g_B,A\in\eta-\eta_g\end{equation}\] 则式（12）变为 \[\begin{equation}\sum_{B\in\eta-\eta_g}K_{AB}d_B=F_A\end{equation}\] 进一步使用矩阵形式表述为 \[\begin{equation}\boldsymbol{K}=[K_{AB}]= \left[ \begin{matrix} K_{0,0}&K_{0,1}&\cdots&K_{0,n-1}\\ K_{1,0}&K_{1,1}&\cdots&K_{1,n-1}\\ \cdots&\cdots&\ddots&\cdots\\ K_{n-1,0}&K_{n-1,1}&\cdots&K_{n-1,n-1} \end{matrix} \right]\end{equation}\] \[\begin{equation}\vec{F}=\{F_A\}= \left\{ \begin{matrix} F_0\\ F_1\\ ...\\ F_{n-1} \end{matrix} \right\}\end{equation}\] \[\begin{equation}\vec{d}=\{d_B\}= \left\{ \begin{matrix} d_0\\ d_1\\ ...\\ d_{n-1} \end{matrix} \right\} \end{equation}\] 上述矩阵形式表述为 \[\begin{equation}\boldsymbol{K}\vec{d}=\vec{F}\end{equation}\] 式中\(K\)为刚度矩阵，\(\vec{F}\)为受力项，\(\vec{d}\)为位移矢量。矩阵形式表述为
给定系数矩阵\(\boldsymbol{K}\)和向量\(\vec{F}\)，寻找\(\vec{d}\)满足\(\boldsymbol{K}\vec{d}=\vec{F}\),其中\(\boldsymbol{K}=[K_{PQ}]\)，\(\vec{d}=\{d_q\}\)，\(\vec{F}=\{F_p\}\)，\(0\le P,Q<n_{eq}，K_{PQ}=a(N_A,N_B)，P=ID(A),Q=ID(B)，F_P=(N_A,f)+(N_A,h)_\Gamma-\sum_{B\in\eta_{eq}}a(N_A,N_B)g_B\)
需要注意的是，我们只对未知节点\(n_{eq}=n_{np}-n_g\)应用上述等式，我们需要找到节点标号与未知方程标号的对应关系，我们根据Hughes标记法将其称为ID： \[\begin{equation}ID(A)= \begin{cases} P&A\in\eta-\eta_g\\ 0&A\in\eta_g \end{cases} \end{equation}\]
6 单元视角
上述表述使用了全局视角来表述矩阵问题。在这个视角中，空间位置和插值函数\(N_A\)在物理\(\vec{x}\)坐标中给出。这种表述方式并不实用，因此我们改为使用单元视角。我们让每个单元覆盖一些逻辑域，其坐标为\(\vec{\xi} = \{\xi,\eta,\zeta\}\)。现在，插值函数\(N_{0,1,\cdots}\)以这些逻辑或自然坐标表示。整个域的积分现在可以重写为对所有单元的求和。 \[\begin{equation}\boldsymbol{K}=\sum^{n_{el}-1}_{e=0}\boldsymbol{K}^e,\boldsymbol{K}^e=[K^e_{AB}]\end{equation}\] \[\begin{equation}\vec{F}=\sum^{n_{el}-1}_{e=0}\vec{F}_e,\vec{F}_e=\{F^e_A\}\end{equation}\] 由上标\(e\)标记的\(K_{AB}\)和\(F_A\)这两个定义的主要区别是在单元视角中，只在单元内进行积分 \[\begin{equation}K^e_{AB}=a(N_A,N_B)^b=\int_{\Omega^e}\nabla N_A k\nabla N_B d\Omega\end{equation}\] \[\begin{equation}F^e_p=(N_A,f)^e+(N_A,h)^e_\Gamma-\sum_{B\in\eta_g}a(N_A,N_B)^eg_B =\int_{\Omega^e}N_Afd\Omega+\int_{\Gamma^e_h}N_Ahd\Gamma-\sum_{B\in\eta_g}a(N_A,N_B)^eg_B\end{equation}\]
在有限元方法中，矩阵方程的单元表述涉及到刚度矩阵的各个单元。\(\Gamma^e_h\)是单元\(e\)在边界\(\Gamma_h\)上的部分，这个边界仅对于位于域边界的单元非空。此外\(P = ID(A)\)，\(Q = ID(B)\)。这里的ID表示节点编号与未知方程编号之间的映射表。
对单个单元重写上述公式为 \[\begin{equation}k^e=[k^e_{ab}],f^e=\{f^e_a\},0\le a,b<n_{en}\end{equation} \] \[\begin{equation}k^e_{ab}=a(N_a,N_b)^e=\int_{\Omega^e}(\nabla N_a)^Tk\nabla N_bd\Omega \end{equation}\] \[\begin{equation} f^e_a=\int_{\Omega^e}N_afd\Omega+\int_{\Omega^e}N_ahd\Gamma-\sum^{n_{en}}_{b=1}k^e_{ab}q^e_b\end{equation}\] 这些等式与式23的区别式这里的\(a\)和\(b\)的范围在\(n_{en}\)个单元节点。在线性四面体中，我们有\(4\)个节点，得到\(4\times 4\)个矩阵。
7 正交
现在只考虑在区域\(\Omega^e\)上对任意函数\(f(x,y,z)\)进行积分。假设存在一个从物理坐标到自然（或逻辑）坐标\(\vec{\eta}=(\xi,\eta,\zeta)\)的映射，使得： \[\begin{equation} x(\vec{\eta})=\sum^4_{a=1}N_a(\vec{\xi})\vec{x}^e_a \end{equation}\] 式中\(\vec{x}^e_a\)是组成单元\(e\)的四个节点的位置，且\(\xi，\eta，\zeta\in[-1:1]\)。这个公式描述了物理坐标可以作为逻辑坐标的函数，可用于对节点位置的插值计算。接下来将积分转换为\(\vec{\xi}\)空间下的积分： \[\begin{equation} \int_{\Omega^e}f(\vec{x})dxdydz=\int^1_{-1}\int^1_{-1}\int^1_{-1}f(x(\vec{\xi}),y(\vec{\xi}),z(\vec{\xi}))j(\vec{\xi})d\xi d\eta d\zeta\end{equation}\] 式中\(j(\vec{\xi})=\partial x_i/\partial\xi_i\)称为雅可比矩阵。
现在只考虑积分\(\int_{-1}^1g(\xi)d\xi\)，在这一维度上的积分非常简单，即使用梯形公式： \[ \begin{equation} \int_{-1}^1g(\xi)d\xi=\frac{1}{2}\left(g(-1)+g(1)\right)(1+1)+R=g(-1)+g(1)+R \end{equation}\]
这是数值积分算法， 我们可以写出上式的普遍形式 \[\begin{equation} \int_{-1}^1g(\xi)d\xi=\sum_{l=1}^{n_{int}}g(\widetilde{\xi}W_l)+R\end{equation}\] 对于有两个积分点\(n_{int}=2\)使用梯形法则，我们有\(\widetilde{\xi}_1=-1\)，\(\widetilde{\xi}_2=-1\)，\(W_1=1\)，\(W_2=1\)。残差项为\(R=-(2/3)g_{,\xi\xi}(\vec{\xi})\)，当这项为0时，我们得到准确的数值积分值，如果\(g\)线性变化，则\(g_{,\xi\xi}\)将消失。线性形状函数有一个主要的问题，是他们在每个单元上的导数都是常数，因此建议考虑高阶单元。通过使用更多的控制点可以改善计算。一种优化策略是通过使用高斯正交的办法采样。使用高斯插值系数，我们可以计算如下的高维积分。 \[ \begin{equation} \int_{-1}^1\int_{-1}^1\int_{-1}^1g(\xi,\eta,\zeta)d\xi d\eta d\zeta\approx\sum_{n=1}^{n_{int}}\sum_{m=1}^{n_{int}}\sum_{l=1}^{n_{int}}g(\vec{\xi_l},\vec{xi_m},\vec{xi_n})W_lW_mW_n \end{equation}\] 对任意的极限，则有 \[\begin{equation} \int_a^bg(\xi)d\xi\approx\frac{a-b}{2}\sum_{l=1}^{n_{int}}g(\frac{b-a}{2}\vec{\xi_{l}}+\frac{a+b}{2})W_l\end{equation} \]
8 形状函数及坐标转换
下面求解 \[\begin{equation}k^e_{ab}=a(N_a,N_b)^e=\int_{\Omega^e}N_{a,i}kN_{b,i}d\Omega\end{equation}\]
我们可以在逻辑坐标空间\(\vec{\xi}\)中计算。使用链式法则， \[\begin{equation}\frac{\partial N_a}{\partial
x}=\frac{\partial N_a}{\partial\xi}\frac{\partial\xi}{\partial
x}+\frac{\partial N_a}{\partial\eta}\frac{\partial\eta}{\partial
x}\end{equation}\] 或使用爱因斯坦标记法表示为 \[\begin{equation}N_{a,x}=N_{a,\xi}\xi_{,x}+N_{a,\eta}\eta_{,x}\end{equation}\]
相似地， \[\begin{equation}N_{a,y}=N_{a,\xi}\xi_{,y}+N_{a,\eta}\eta_{,y}\end{equation}\]
从上式可以得到，我们需要获得\(N_a\)在逻辑空间中的形状函数，\(N_a=N_a(\xi)\)。我们也需要映射\(\vec{\xi}_{,\vec{x}}\)的表达式，于是有
\[\begin{equation}\int_\Omega
f(\vec{x})d\Omega=\int_\xi
f(\vec{\xi})|\vec{x}_{,\vec{\xi}}|d\vec{\xi}\end{equation}\]
现在将三角形映射到四面体上，想象在\(\{x,y\}\)空间上的某一三角形，我们需要一个坐标转换来将每个点映射到\(\{\xi,\eta\}\)空间上的每个点。我们的目标是将三角形映射到\(\{\xi,\eta\}\in[0:1]\)。
一种方式是使用三角坐标，这些坐标是基于面积形状因子\(L_i\)。为了保持已有的标记，我们将这些函数定义为\(N_i\)。考虑\(N_0\)，我们有\(N_0(\vec{x_0}=1)\)，\(N_1(\vec{x_1})=N_1(\vec{x_2})=0\)。一般地，
\[
\begin{equation}
N_i(\vec{x_i})=\begin{cases}
1:i=j\\
0:i\neq j
\end{cases}
\end{equation}
\] 由于\(N_0+N_1+N_2=1\)，我们可以定义 \[
\begin{equation}
N_0=\xi
\end{equation}\] \[
\begin{equation}
N_1=\eta
\end{equation}\] \[
\begin{equation}
N_2=1-\xi-\eta
\end{equation}\] 回到式（36）和（37），我们可以得到\(N_{a,\xi}\)和\(N_{a,\eta}\)的表达式： \[
\begin{equation}
\begin{matrix}
N_{0,\xi}=1&N_{0,\eta}=0\\
N_{1,\xi}=0&N_{1,\eta}=1\\
N_{2,\xi}=-1&N_{2,\eta}=-1
\end{matrix}
\end{equation}\] 将式（36）和式（37）重写为矩阵形式 \[\begin{equation}\left\{
\begin{matrix}
N_{a,x}\\
N_{a,y}
\end{matrix}
\right\}=
\left\{
\begin{matrix}
N_{a,\xi}\\
N_{a,\eta}
\end{matrix}
\right\}
\left[
\begin{matrix}
\xi_x&\xi_y\\
\eta_x&\eta_y
\end{matrix}
\right]\end{equation}
\] 根据\(x(\xi,\eta)=\sum^2_{a=0}N_a(\xi,\eta)x^e_a，y(\xi,\eta)=\sum^2_{a=0}N_a(\xi,\eta)y^e_a\)，计算得到
\[\begin{equation}\vec{x}_{,\vec{\xi}}=
\left[
\begin{matrix}
x_{,\xi}&x_{,\eta}\\
y_{,\xi}&y_{,\eta}
\end{matrix}
\right]\end{equation}
\] 矩阵中的每一项按下式计算 \[
\begin{equation}
x_{,\xi}=\sum_{a=0}^2N_{a,\xi}x_a^ex_{,\eta}=\sum_{a=0}^2N_{a,\eta}x_a^e
\end{equation}
\] \[
\begin{equation}
y_{,\xi}=\sum_{a=0}^2N_{y,\xi}y_a^ey_{,\eta}=\sum_{a=0}^2N_{a,\eta}y_a^e
\end{equation}
\] 注意到我们已有得到了形状函数的积分形式，即\(N_{a,\xi}\)和\(N_{a,\eta}\) 而\(\vec{\xi_{,\vec{x}}}=(\vec{x}_{,\xi})^{-1}\)。
根据以上内容，计算任意单元内\(N_{a,x}\)和\(N_{b,x}\)的步骤如下：
1.计算正交点\(l_{\xi}\)处的\(N_{a,\xi}\)；
2.使用\(N_{a,\xi}\)计算\(\vec{x}_{,\vec{\xi}}\);
3.计算\(\vec{\xi_{,\vec{x}}}=(\vec{x}_{,\xi})^{-1}\)；
4.使用\(N_{a,\xi}\xi_{,x}\)计算\(N_{a,x}\)
5.对b重复1-4。
在PIC算法中，我们通常对解的微分感兴趣，即\(\vec{E}=-\nabla\phi\)，使用有限元的形式，我们可以得到 \[ \begin{equation} E(\vec{x})=-\sum_{n=0}^{n_{en}-1}N_{a,x}(\vec{x})\phi_a^e \end{equation} \] 在本章中的形状函数均为线性，因此每个单元处的电场都为常数，但实际过程中应该使用高阶形式。 将该算法转移到三维四面体单元比较直接，我们需要有额外一个维度，并映射坐标\(\{x,y,z\}\rightarrow\{\xi,\eta,\zeta\}\)，我们仍使用归一化体作为自然四面体坐标，即 \[\begin{equation} N_0=\xi,N_1=\eta,N2=\zeta,N_3=1-\xi-\eta-\zeta \end{equation}\] 其中 \[ \begin{equation} N_0(x)\equiv L_0(x)=\frac{V_{x123}}{V_{0123}} \end{equation}\] 于是有 \[ \begin{equation} \begin{matrix} N_{0,\xi}=1&N_{0,\eta}=0&N_{0,\zeta}=0\\ N_{1,\xi}=0&N_{1,\eta}=1&N_{1,\zeta}=0\\ N_{2,\xi}=0&N_{2,\eta}=0&N_{2,\zeta}=1\\ N_{3,\xi}=-1&N_{3,\eta}=-1&N_{3,\zeta}=-1\\ \end{matrix} \end{equation}\]
