2021-5-19 有限元方法从0开始学习笔记 第三天 (有限元方法计算处理简介)

又是新的一天, 昨天上午有课,今天还要备课, 明天要讲两次课, 因此这三天可能只能写一点点了, 反正能写多少是多少把. 继续给出参考资料
Brenner S, Scott R. The mathematical theory of finite element methods[M]. Springer Science & Business Media, 2007.
从第十页.

昨日内容简单回忆

2021-5-18 有限元方法从0开始学习笔记 第二天 (一维线性元简介, 与差分格式的联系)

• 给出了第一个具体的有限元空间, 一维的线性分片多项式空间, 在该空间下利用给出了 Risz-Galerkin 框架的具体刚度矩阵, 以及给出了相应的插值误差估计.
• 讲述了一维线性元与差分格式的一些联系.

0.6 有限元方法的计算处理

本节主要介绍对前面一维线性元的编程处理, 有限元处理的一个重要的观点即用局部刚度矩阵来组装整体刚度矩阵, 为此定义指标 $i(e,j)$, 其中各个记号的解释如下(这里翻译过来感觉有点别扭, 默认 j 表示局部节点编号吧, 因为书上写的 node number, 反正这里只考虑线性元, 暂时先不去引入自由度的概念.):

• $e$:表示单元的整体编号, $I_e=[x_{e-1},x_e]$.
• $j$ : 表示单元的局部节点编号, $j=0,1$ . (用 $0$表示左端点, $1$表示右端点)
• $i(e,j):$ 表示节点的全局编号,

显然, 节点的全局编号是关于局部编号以及单元整体编号的函数, 并且这里可以看出, 有
$$i(e,j)=e+j-1,e=1,2,\cdots ,n,j=0,1,$$
下面介绍有限元基本的编程思路, 基本的思想是:
参考单元上的计算 $\Rightarrow$ 局部单元上计算 $\Rightarrow$ 整体组装.
网格剖分的过程可以参考上一小节的内容, 并且记号也将沿用上一小节的记号, 因此这里不作过多的赘述, 假设 $f$ 为连续函数, 考虑 $f$ 在分片多项式空间上的插值函数:
$$f_I(x)=\sum _{i=1}^{n}f(x_i){\varphi }_i(x),$$
其中, ${\varphi }_i(x)$ 是前面所定义的帽子函数, 很显然, 上面这个式子是不利于编程的,(因为这些帽子函数需要分段地定义出来) 我们将上式子改写如下:
$$f_I(x)=\sum _e\sum _{j=0}^{1}f(x_{i(e,j)}){\psi }_j^e(x),$$
其中, ${\psi }_j^e(x),j=0,1$ 分别定义如下:
$${\psi }_0^e(x)=\{\begin{array}{ll}\frac{x_e-x}{x_e-x_{e-1}}& x\in I_e\\ 0& \text{otherwise}\end{array},{\psi }_1^e(x)=\{\begin{array}{ll}\frac{x-x_{e-1}}{x_e-x_{e-1}}& x\in I_e\\ 0& \text{otherwise}\end{array}.$$
这两个函数合起来, 称为 单元 $I_e$ 上的局部基函数.
注意, 后面改写的插值在网格节点上不再有 $f(x_i)=f_I(x_i)$, 但是这不影响积分的计算, 我的理解这是看成在 $L^2$的意义下, 新的插值函数和原来的插值函数相等? 这里只是为了方便编程时积分的计算.
此外, 一般的, 也不会直接计算 ${\psi }_j^e(x)$ 的函数值, 而是定义如下参考单元上的局部基函数:
$${\hat{\psi }}_0(\hat{x})=\{\begin{array}{ll}1-\hat{x}& \hat{x}\in \left[0,1\right]\\ 0& \text{otherwise}\end{array},{\hat{\psi }}_1(\hat{x})=\{\begin{array}{ll}\hat{x}& \hat{x}\in \left[0,1\right]\\ 0& \text{otherwise}\end{array}.$$
再通过定义仿射变换将 $\hat{x}\in \left[0,1\right]$ 映射到 $x\in \left[x_{e-1},x_e\right]$ 上如下:
$$x=x_{e-1}+\left(x_e-x_{e-1}\right)\hat{x},$$
从而
$$\hat{x}=\frac{x-x_{e-1}}{x_e-x_{e-1}},$$
因此, 可以改写公式
$$f_I(x)=\sum _e\sum _{j=0}^{1}f(x_{i(e,j)}){\psi }_j^e(x),$$
为:
$$f_I(x)=\sum _e\sum _{j=1}^{1}f(x_{i(e,j)}){\hat{\psi }}_j\left(\frac{x-x_{e-1}}{x_e-x_{e-1}}\right):=\sum _e\sum _{j=1}^{1}f(x_{i(e,j)}){\hat{\psi }}_j\left({\hat{x}}_e\right),$$
利用上述操作, 可以容易地计算如下的双线性函数, 定义 $a_e(v,w)={\int }_{x_{e-1}}^{x_e}{v}^{{}^{\prime }}(x)w^{{}^{\prime }}(x)dx$ ,
$$a(v,w)=\sum _e{a}_e(v,w)=\sum _e{\int }_{x_{e-1}}^{x_e}{v}^{\prime }(x)w^{\prime }(x)dx,$$
下面计算 ${\int }_{x_{e-1}}^{x_e}{v}^{\prime }(x)w^{\prime }(x)dx$:
$$\begin{array}{rl}{\int }_{x_{e-1}}^{x_e}{v}^{{}^{\prime }}(x)w^{{}^{\prime }}(x)dx& ={\int }_{x_{e-1}}^{x_e}\left(\sum _j{v}_{i(e,j)}{\psi }_j^{e^{\prime }}(x)\right)\left(\sum _k{w}_{i(e,k)}{\psi }_k^{e^{\prime }}(x)\right)dx\\ & =\frac{1}{x_e-x_{e-1}}{\int }_0^1\left(\sum _j{v}_{i(e,j)}{\hat{\psi }}_j^{{}^{\prime }}({\hat{x}}_e)\right)\left(\sum _k{w}_{i(e,k)}{\hat{\psi }}_k^{{}^{\prime }}({\hat{x}}_e)\right)d{\hat{x}}_e\\ & =\frac{1}{x_e-x_{e-1}}\sum _j\sum _k{v}_{i(e,j)}{w}_{i(e,k)}{\int }_0^1{\hat{\psi }}_j^{{}^{\prime }}({\hat{x}}_e){\hat{\psi }}_k^{{}^{\prime }}({\hat{x}}_e)d{\hat{x}}_e\end{array}$$
引入记号:
$$\begin{array}{r}{\mathbf{v}}_e=\left(\begin{array}{c}{v}_{i(e,0)}\\ v_{i(e,1)}\end{array}\right),{\mathbf{w}}_e=\left(\begin{array}{c}{w}_{i(e,0)}\\ w_{i(e,1)}\end{array}\right),{\mathbf{K}}_e=\left(\begin{array}{cc}{\int }_0^1{\hat{\psi }}_0^{{}^{\prime }}({\hat{x}}_e){\hat{\psi }}_0^{{}^{\prime }}({\hat{x}}_e)d{\hat{x}}_e& {\int }_0^1{\hat{\psi }}_0^{{}^{\prime }}({\hat{x}}_e){\hat{\psi }}_1^{{}^{\prime }}({\hat{x}}_e)d{\hat{x}}_e\\ {\int }_0^1{\hat{\psi }}_1^{{}^{\prime }}({\hat{x}}_e){\hat{\psi }}_0^{{}^{\prime }}({\hat{x}}_e)d{\hat{x}}_e& {\int }_0^1{\hat{\psi }}_1^{{}^{\prime }}({\hat{x}}_e){\hat{\psi }}_1^{{}^{\prime }}({\hat{x}}_e)d{\hat{x}}_e\end{array}\right)\end{array}$$
从而有:
$${\int }_{x_{e-1}}^{x_e}{v}^{\prime }(x)w^{\prime }(x)dx=\frac{1}{x_i-x_{i-1}}{\mathbf{v}}_e^{⊺}{\mathbf{K}}_e{\mathbf{w}}_e$$
上面过程是一个向量化的过程, 可以方便地编程计算, 对于整个线性元的计算, 完全可以按照上面的方法进行推导并且编程,并且上述编程思路可以推广到高次元以及高维问题, 当然, 这里暂时先不去讨论. 有兴趣的可以按照上面过程自行推导整个一维线性元的计算流程, 并编写程序实现, 后面如果有时间也会去实现一下, 并在这里附上链接.
赶紧溜了, 这两天就写到这里, 后面还要备课.