简易分析，对流扩散方程，在四边形等参网格中的有限元离散过程，以及数值积分

1、生成弱形式

a.等式两侧同时乘测试函数$\psi$

$$\psi (-\nabla \cdot \nabla u)+\psi (\vec{v}\cdot \nabla u)=0\psi ,u\in V\text{\hspace{1em}(2.1)}$$
这里$V=H^1(\Omega )$

b.对整个求解域$\Omega$积分

$$-{\int }_{\Omega }\psi (\nabla \cdot \nabla u)+{\int }_{\Omega }\psi (\vec{v}\cdot \nabla u)=0\text{\hspace{1em}(2.2)}$$

c.分部积分

由于求导的性质有$\nabla \cdot (\psi \nabla u)=\psi \nabla \cdot \nabla u+\nabla \psi \cdot \nabla u$，带入（2.2）式左侧第一项得:
$${\int }_{\Omega }\nabla \psi \cdot \nabla u-{\int }_{\Omega }\nabla \cdot (\psi \nabla u)+{\int }_{\Omega }\psi (\vec{v}\cdot \nabla u)=0\text{\hspace{1em}(2.3)}$$

d.应用散度定理

散度定理可以理解成：考虑在流动得液体中取出体积为$V$的部分，则流过该部分表面$S$的液体通量等于散度在整个体积上的积分。
$${\iiint }_V\left(\nabla \cdot \mathbf{F}\right)\mathrm{d}V=\iint s(\mathbf{F}\cdot \hat{\mathbf{n}})\mathrm{d}S$$
左侧是在整个体积$V$上的积分，右侧是在$V$的表面边界$S$上的积分,$\hat{\mathbf{n}}$是边界上外法线方向。
对（2.3）式的左侧第二项应用散度定理，（2.3）式可以写成
$${\int }_{\Omega }\nabla \psi \cdot \nabla u-{\int }_{\partial \Omega }\psi (\nabla u\cdot \hat{\mathbf{n}})+{\int }_{\Omega }\psi (\vec{v}\cdot \nabla u)=0\psi ,u\in V\text{\hspace{1em}(2.4)}$$

e.写成内积的形式

可以把（2.4）式用内积形式表示，其中每项都能在MOOSE中找到现有类进行继承
$$\underset{Kernel}{\underbrace{\left(\nabla \psi ,\nabla u\right)}}-\underset{BoundaryCondition}{\underbrace{⟨\psi ,\nabla u\cdot \hat{n}⟩}}+\underset{Kernel}{\underbrace{\left(\psi ,\vec{v}\cdot \nabla u\right)}}=0\psi ,u\in V\text{\hspace{1em}(2.5)}$$

2、Rizt—Galerkin法

现在我们有了方程的弱形式，但是考虑到真实解$u$所在的空间$V$是无线维的，在大多情况下依旧难以求解，根据Ritz方法，可以选择一个$V$的有限维子空间$W$，即$W$的基 { w 1 , w 2 , . . . w n } \left\{w_{1}, w_{2},...w_{n} \right\} {w1​,w2​,...wn​}是有限个，称为基函数，则真实解$u$可以由W的基构成的$u_h$进行逼近表示
$$u\approx u_h=\sum _i^n{U}_i{w}_i$$
为了方便记
$$a(\psi ,u)=\left(\nabla \psi ,\nabla u\right)-⟨\psi ,\nabla u\cdot \hat{n}⟩+\left(\psi ,\vec{v}\cdot \nabla u\right)$$
则逼近解替换（2.5）式中的真实解后
$$\begin{array}{rl}& \left(\nabla \psi ,\nabla u_h\right)-⟨\psi ,\nabla u_h\cdot \hat{n}⟩+\left(\psi ,\vec{v}\cdot \nabla u_h\right)\\ & =\left(\nabla \psi ,\sum _i^n{U}_i\nabla w_i\right)-⟨\psi ,\sum _i^n{U}_i\nabla w_i\cdot \hat{n}⟩+\left(\psi ,\vec{v}\cdot \sum _i^n{U}_i\nabla w_i\right)\\ & =\sum _i^n{U}_i\left(\nabla \psi ,\nabla w_i\right)-⟨\psi ,\nabla w_i\cdot \hat{n}⟩+\left(\psi ,\vec{v}\cdot \nabla w_i\right)\\ & ={\left[\begin{array}{c}\left(\nabla \psi ,\nabla w_1\right)-⟨\psi ,\nabla w_1\cdot \hat{n}⟩+\left(\psi ,\vec{v}\cdot \nabla w_1\right)\\ \left(\nabla \psi ,\nabla w_2\right)-⟨\psi ,\nabla w_2\cdot \hat{n}⟩+\left(\psi ,\vec{v}\cdot \nabla w_2\right)\\ ⋮\\ ⋮\\ \left(\nabla \psi ,\nabla w_n\right)-⟨\psi ,\nabla w_n\cdot \hat{n}⟩+\left(\psi ,\vec{v}\cdot \nabla w_n\right)\end{array}\right]}^{\mathrm{T}}\left[\begin{array}{c}{U}_1\\ U_2\\ ⋮\\ ⋮\\ U_n\end{array}\right]\\ & =\left[\begin{array}{c}a(\psi ,w_1),a(\psi ,w_2),\cdots ,a(\psi ,w_n)\end{array}\right]\left[\begin{array}{c}{U}_1\\ U_2\\ ⋮\\ ⋮\\ U_n\end{array}\right]\\ & =0\end{array}\text{\hspace{1em}(2.6)}$$
有限维子空间$W$可以自由选取，所以基函数$w_i$可以当作已知量，但是基的系数$U_i$是未知数，共有$n$个，需要我们计算，由于测试函数的任意性，我们可以选用基函数$w_i$当作试函数，我们就能得到$n$个方程组，这样$n$个方程$n$个未知数，我们就能够求解$U_i$了。
$$\left(\nabla w_i,\nabla u_h\right)-⟨w_i,\nabla u_h\cdot \hat{n}⟩+\left(w_i,\vec{v}\cdot \nabla u_h\right)=0i=1,2,\cdots ,n\text{\hspace{1em}(2.7)}$$
把（2.7）式写成矩阵的形式
$$\begin{array}{rl}& \left(\nabla \psi ,\nabla u\right)-⟨\psi ,\nabla u\cdot \hat{n}⟩+\left(\psi ,\vec{v}\cdot \nabla u\right)\\ & =\left[\begin{array}{c}a(\psi ,w_1),a(\psi ,w_2),\cdots ,a(\psi ,w_n)\end{array}\right]\left[\begin{array}{c}{U}_1\\ U_2\\ ⋮\\ ⋮\\ U_n\end{array}\right]\\ & =\left[\begin{array}{cccc}a(w_1,w_1)& a(w_1,w_2)& \cdots & a(w_1,w_n)\\ a(w_2,w_1)& a(w_2,w_2)& \cdots & a(w_3,w_n)\\ ⋮& ⋮& \ddots & ⋮\\ & & & \\ a(w_n,w_1)& a(w_n,w_2)& \cdots & a(w_n,w_n)\end{array}\right]\left[\begin{array}{c}{U}_1\\ U_2\\ ⋮\\ ⋮\\ U_n\end{array}\right]\\ & =KU\\ & =0\end{array}\text{\hspace{1em}(2.8)}$$
这样就得到了一个简单的写法
$$KU=0\text{\hspace{1em}(2.9)}$$
其中$U$称作载荷向量，$K$称为刚度矩阵,我们的目标是求解$U$

3、分段多项式基函数

这一小节的讨论都基于规则四边形网格，就是求解域和子域都是矩形。
求解为了快速有效的求解式（2.9），必须选择一个很好的逼近子空间，我们可以选择连续分段多项式函数组成的子空间作为逼近子空间。为了在域$\Omega$上定义分段多项式函数，当作逼近子空间$W$的基，将域划分为子域${\Omega }_e$，分段多项式是由每个子域上的多项式定义的函数，子域的集合称为网格。在本例中选择了四边形网格剖分，每个四边形有四个顶点，可以考虑一个具有四个自由度的多项式
$$a+bx+cy+dxy$$
这样的多项式函数，称为双线性函数，对于一个平行于坐标轴的四边形，设四个顶点坐标为$(x_1,y_1),(x_2,y_1),(x_2,y_2),(x_1,y_2)$以及相应的值$u_1,u_2,u_3,u_4$，我们可以得到四个等式
$$\begin{array}{r}a+b{x}_1+c{y}_1+d{x}_1{y}_1=u_1\\ a+b{x}_2+c{y}_1+d{x}_2{y}_1=u_2\\ a+b{x}_2+c{y}_2+d{x}_2{y}_2=u_3\\ a+b{x}_1+c{y}_2+d{x}_1{y}_2=u_4\end{array}$$
这样可以求出定义在该四边形区域上的多项式的系数$a,b,c,d$，如果网格单元是规则的矩形，双线性函数在矩形单元的任何边上都可以简化为一元线性函数。因此，在规则的矩形剖分上，网格的每个节点上都能确定一个连续分段线双线性函数，这样在剖分出的$n$个节点上我们定义出逼近子空间的基函数 { ψ 1 , ψ 2 , . . . ψ n } \left\{\psi_{1}, \psi_{2},...\psi_{n} \right\} {ψ1​,ψ2​,...ψn​}，分段多项式函数的定义方式，所以不在同一个矩形的顶点上定义的${\psi }_i,{\psi }_i$乘积为零。原本定义在整个域的方程，现在可以定义到每个子域上求解然后再求和，则（2.4）式化成
$$\begin{array}{rlr}& {\int }_{\Omega }\nabla \psi \cdot \nabla u-{\int }_{\partial \Omega }\psi (\nabla u\cdot \hat{\mathbf{n}})+{\int }_{\Omega }\psi (\vec{v}\cdot \nabla u)& \\ =& \sum _e\left[{\int }_{{\Omega }_e}\nabla {\psi }_i\cdot \nabla u_h-{\int }_{\partial \Omega }{\psi }_i(\nabla u_h\cdot \hat{\mathbf{n}})+{\int }_{{\Omega }_e}{\psi }_i(\vec{v}\cdot \nabla u_h)\right]& \end{array}$$
)
现在我们考虑一个9节点的网格，其共有4个子域，则刚度矩形$K$就是$9\times 9$的$$\begin{array}{rl}& {\int }_{\Omega }\nabla \psi \cdot \nabla u-{\int }_{\partial \Omega }\psi (\nabla u\cdot \hat{\mathbf{n}})+{\int }_{\Omega }\psi (\vec{v}\cdot \nabla u)\\ =& \sum _{e=1}^4\left[{\int }_{{\Omega }_e}\nabla {\psi }_i\cdot \nabla u_h-{\int }_{\partial \Omega }{\psi }_i(\nabla u_h\cdot \hat{\mathbf{n}})+{\int }_{{\Omega }_e}{\psi }_i(\vec{v}\cdot \nabla u_h)\right]\\ =& a_{{\Omega }_1}({\psi }_i,u)+a_{{\Omega }_2}({\psi }_i,u)+a_{{\Omega }_3}({\psi }_i,u)+a_{{\Omega }_4}({\psi }_i,u)i=1,2,\cdots ,9\end{array}$$
在区域${\Omega }_1$中，包含四个写成刚度矩阵的形式
$$\begin{array}{rl}& a_{{\Omega }_1}(\psi ,u)+a_{{\Omega }_2}(\psi ,u)+a_{{\Omega }_3}(\psi ,u)+a_{{\Omega }_4}(\psi ,u)\\ =& \sum _{e=1}^4{\left[\begin{array}{ccccccccc}{K}_{11}& K_{12}& K_{13}& K_{14}& K_{15}& K_{16}& K_{17}& K_{18}& K_{19}\\ K_{21}& K_{22}& K_{23}& K_{24}& K_{25}& K_{26}& K_{27}& K_{28}& K_{29}\\ K_{31}& K_{32}& K_{33}& K_{34}& K_{35}& K_{36}& K_{37}& K_{38}& K_{39}\\ K_{41}& K_{42}& K_{43}& K_{44}& K_{45}& K_{46}& K_{47}& K_{48}& K_{49}\\ K_{51}& K_{52}& K_{53}& K_{54}& K_{55}& K_{56}& K_{57}& K_{58}& K_{59}\\ K_{61}& K_{62}& K_{63}& K_{64}& K_{65}& K_{66}& K_{67}& K_{68}& K_{69}\\ K_{71}& K_{72}& K_{73}& K_{74}& K_{75}& K_{76}& K_{77}& K_{78}& K_{79}\\ K_{81}& K_{82}& K_{83}& K_{84}& K_{85}& K_{86}& K_{87}& K_{88}& K_{89}\\ K_{91}& K_{92}& K_{93}& K_{94}& K_{95}& K_{96}& K_{97}& K_{98}& K_{99}\end{array}\right]}_{{\Omega }_e}\left[\begin{array}{c}{U}_1\\ U_2\\ ⋮\\ ⋮\\ U_n\end{array}\right]\\ =& {\left[\begin{array}{ccccccccc}{K}_{11}& K_{12}& 0& K_{14}& K_{15}& 0& 0& 0& 0\\ K_{21}& K_{22}& 0& K_{24}& K_{25}& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 0\\ K_{41}& K_{42}& 0& K_{44}& K_{45}& 0& 0& 0& 0\\ K_{51}& K_{52}& 0& K_{54}& K_{55}& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 0\end{array}\right]}_{{\Omega }_1}\left[\begin{array}{c}{U}_1\\ U_2\\ ⋮\\ ⋮\\ U_n\end{array}\right]\\ +& {\Omega }_2+{\Omega }_3+{\Omega }_4(\Omega 234区域的刚度矩阵与1类似，不同时包含区域节点上积函数的元素为零)\end{array}\text{\hspace{1em}(2.10)}$$

4、等参单元

多数情况的求解域都不是规则的矩形，所以为了减少计算复杂度，同时让定义的基函数有效，可以引入等参单元${\hat{\Omega }}_e$，构造一个与坐标轴平行的规则矩形，把要在网格中求解的不规则四边形${\Omega }_e$通过仿射变换映射到该等参单元上求解。把等参单元的坐标用$(\xi ,\eta )$表示，那么$(\xi ,\eta )$到$(x,y)$的映射$F(\xi ,\eta )$可以表示为
$$\begin{array}{rl}x& =a_1+a_2\xi +a_3\eta +a_4\xi \eta \\ y& =b_1+b_2\xi +b_3\eta +b_4\xi \eta \end{array}\text{\hspace{1em}(2.11)}$$
等参单元${\hat{\Omega }}_e$的四个顶点坐标可以规定为$(-1,-1),(1,-1),(1,1),(-1,1)$,设${\Omega }_e$的坐标为$(x_1,y_1),(x_2,y_2),(x_3,y_3),(x_4,y_4)$，带入式（3.1）可以得到
$$\begin{array}{rl}{a}_1-a_2-a_3+a_4& =x_1\\ a_1+a_2-a_3-a_4& =x_2\\ a_1+a_2+a_3+a_4& =x_3\\ a_1-a_2+a_3-a_4& =x_4\end{array}\text{\hspace{1em}(2.12)}$$
可以发现关于$x$和$y$的方程组系数是一样的
$$M=\left[\begin{array}{cccc}1& -1& -1& 1\\ 1& 1& -1& -1\\ 1& 1& 1& 1\\ 1& -1& 1& -1\end{array}\right]$$
所以可以得到
$$Ma=xMb=y$$
这样就很容易的的到$a,b$的值，$(x,y)$和$(\xi ,\eta )$之间就能相互转化。然后就可以定义等参单元${\hat{\Omega }}_e$上的形状函数，我们可以选择山形基函数，在当前节点为1其余节点为0，对于等参单元的（-1，-1）点，根据拉格朗日插值的想法（$\xi$或$\eta$等于1的时要等于0），可以设形状函数${\gamma }_1$为$c(\xi -1)(\eta -1)$，其中$c$为常数，当$\xi$=$\eta =-1$时基函数为1，则${\gamma }_1=1$，解得$c=\frac{1}{4}$，对于等参单元${\hat{\Omega }}_e$另外三个顶点，采用相同的想法，可以构造出四个形状函数：
$$\begin{array}{r}{\gamma }_1=\frac{1}{4}（1-\xi -\eta +\xi \eta ）\\ {\gamma }_2=\frac{1}{4}（1+\xi -\eta -\xi \eta ）\\ {\gamma }_3=\frac{1}{4}（1+\xi +\eta +\xi \eta ）\\ {\gamma }_4=\frac{1}{4}（1-\xi +\eta -\xi \eta ）\end{array}\text{\hspace{1em}(2.13)}$$
这样通过坐标转换就能得到${\Omega }_e$上的形状函数，值得注意的是${\hat{\Omega }}_e$上的形状函数是双线性的，但是${\Omega }_e$上的形状函数通常不是双线性的。
这样在网格的每个节点都能计算出相应的基函数${\varphi }_i$，这时再看$KU=0$，我们可以用${\varphi }_i$当作K的逼近子空间的基，局部刚度矩阵
$$K_{ij}=\sum _e\left(w_i,w_j\right)=\sum _e{\int }_{{\Omega }_{e_k}}\nabla {\varphi }_i\cdot \nabla {\varphi }_{j}dxdy$$
根据从$(x,y)$到$(\xi ,\eta )$的坐标变换（3.1）式可以求出这个变换的雅可比矩阵
$$J(\xi ,\eta )=F^{\prime }(\xi ,\eta )=\left[\begin{array}{cc}\frac{\partial F_1}{\partial \xi }(\xi ,\eta )& \frac{\partial F_1}{\partial \eta }(\xi ,\eta )\\ \frac{\partial F_2}{\partial \xi }(\xi ,\eta )& \frac{\partial F_2}{\partial \eta }(\xi ,\eta )\end{array}\right]=\left[\begin{array}{cc}{a}_2+a_4\eta & a_3+a_4\xi \\ b_2+b_4\eta & b_3+b_4\xi \end{array}\right]$$
则根据多重积分的性质，对于任意$f(x,y)$的积分有
$${\int }_{\Omega }f(x,y)\mathrm{d}x\mathrm{d}y=\sum _e{\int }_{{\Omega }_e}f(x,y)\mathrm{d}x\mathrm{d}y=\sum _e{\int }_{{\hat{\Omega }}_e}f(\xi ,\eta )\left|{\mathcal{J}}_e\right|\mathrm{d}\xi \mathrm{d}\eta$$

其中$\left|{\mathcal{J}}_e\right|$ ,是等参单元到积分单元的雅可比行列式，再根据（2.7）式，则（2.4）可以写成
$$\begin{array}{rl}a(\psi ,u)=& \left(\nabla \psi ,\nabla u\right)-⟨\psi ,\nabla u\cdot \hat{n}⟩+\left(\psi ,\vec{v}\cdot \nabla u\right)\\ =& {\int }_{\Omega }\nabla \psi \cdot \nabla u-{\int }_{\partial \Omega }\psi (\nabla u\cdot \hat{\mathbf{n}})+{\int }_{\Omega }\psi (\vec{v}\cdot \nabla u)\\ =& \sum _e\left|{\mathcal{J}}_e\right|\left[{\int }_{{\hat{\Omega }}_e}\nabla {\psi }_i\cdot \nabla u_h-{\int }_{\partial \hat{\Omega }}{\psi }_i(\nabla u_h\cdot \hat{\mathbf{n}})+{\int }_{{\hat{\Omega }}_e}{\psi }_i(\vec{v}\cdot \nabla u_h)\right]\\ =& 0\end{array}$$
其中
$$u\approx u_h=\sum _i^n{U}_i{\psi }_{i}i=1,2,\cdots n$$
这样刚度矩阵的

5、高斯积分

进行离散后的积分求解起来很麻烦，为了方便就采用了高斯积分的形式
$$\sum _e{\int }_{{\hat{\Omega }}_e}f(\bar{\xi })\left|{\mathcal{J}}_e\right|\mathrm{d}\bar{\xi }\approx \sum _e\sum _q{w}_{q}f\left({\bar{x}}_q\right)\left|{\mathcal{J}}_e\left({\bar{x}}_q\right)\right|$$
${\bar{x}}_q$是正交点的位置，$w_q$是它的相关权重。

二维的高斯积分点和权重对应表

高斯点个数$n^2(n\times n)$	高斯点坐标	权重	精度
$1(1\times 1)$	$x_1=y_1=1$	$w_1=4$	1
$4(2\times 2)$	$x_i,y_j=\pm \sqrt{\frac{1}{3}}$	$w_{ij}=1$	3
$9(3\times 3)$	$$\begin{gathered} x_1=y_1=0 \\ {x}_{i2},y_{j2}=0,\pm \sqrt{\frac{3}{5}} \\ {x}_{i3},y_{j3}=\pm \sqrt{\frac{3}{5}},0 \\ {x}_{i4},y_{j4}=\pm \sqrt{\frac{3}{5}},\pm \sqrt{\frac{3}{5}} \end{gathered}$$	$$\begin{gathered} w_1=\frac{8}{9}\times \frac{8}{9}=\frac{64}{81} \\ {w}_2=\frac{8}{9}\times \pm \frac{5}{9}=\frac{40}{81} \\ {w}_3=\pm \frac{5}{9}\times \frac{8}{9}=\frac{40}{81} \\ {w}_4=\frac{5}{9}\times \frac{5}{9}=\frac{25}{81} \end{gathered}$$	5

对于四边形等参单元，可以选取9个高斯积分点。当$x,y=0$时权重为$\frac{8}{9}$，当$x,y=\pm \sqrt{\frac{3}{5}}$时权重为$\pm \frac{5}{9}$,那对该等参单元积分就可以写成
$$\begin{array}{rl}{\int }_{{\Omega }_e}f(x,y)\mathrm{d}x\mathrm{d}y& ={\int }_{{\hat{\Omega }}_e}f(\xi ,\eta )\left|{\mathcal{J}}_e\right|\mathrm{d}\xi \mathrm{d}\eta \\ & \approx [\frac{25}{81}f\left(-\sqrt{\frac{3}{5}},\sqrt{\frac{3}{5}}\right)+\frac{40}{81}f\left(0,\sqrt{\frac{3}{5}}\right)+\frac{25}{81}f\left(\sqrt{\frac{3}{5}},\sqrt{\frac{3}{5}}\right)\\ & +\frac{40}{81}f\left(-\sqrt{\frac{3}{5}},0\right)+\frac{64}{81}f\left(0，0)\right)+\frac{40}{81}f\left(\sqrt{\frac{3}{5}},0\right)\\ & +\frac{25}{81}f\left(-\sqrt{\frac{3}{5}},-\sqrt{\frac{3}{5}}\right)+\frac{40}{81}f\left(-\sqrt{\frac{3}{5}},0\right)+\frac{25}{81}f\left(-\sqrt{\frac{3}{5}},\sqrt{\frac{3}{5}}\right)]\left|{\mathcal{J}}_e\right|\end{array}$$

最终我们只需求解:
$$\begin{array}{rl}{\vec{R}}_i(u_h)& =\sum _e\sum _q{w}_q\left|{\mathcal{J}}_e\right|\left[\nabla {\psi }_i\cdot \nabla u_h+{\psi }_i\left(\vec{v}\cdot \nabla u_h\right)\right]({\vec{x}}_q)\\ & -\sum _f\sum _{q_{face}}{w}_{q_{face}}\left|{\mathcal{J}}_f\right|\left[{\psi }_i\nabla u_h\cdot \vec{n}\right]({\vec{x}}_{q_{face}})\}\end{array}$$
这个式子也可以按着（2.6）式的推导方式，写成KU的形式。其中$R(u_h)$是逼近解产生的残差，是逼近解求得的右端项$b_h$与真实的右端项之间的差值.

上文肯定错误百出，希望大家能够指出。我学习改正！