有限元之平面三角单元

有限元之平面三角单元

2020-11-04
| 科学计算 | | julia , 有限元法 | Comment 评论

本文是有限元法基础之一：三角单元。

用 \(\color{red}{红色}\) 标注的公式，都给了完整的符号推导（基于Julia）。

参考文献： 微分方程数值解法（第4版）-李荣华＆刘播-高等教育出版社-2009

三角单元内点的面积坐标

三角形元 \(\Delta(i,j,k)\) 内任意一点p，可用 \((x,y)\) 表示。也可用“归一化面积”，即所谓面积坐标 \((l_i,l_j,l_k)\) 表示。 \(l_i\) 就是 \(\Delta(j,k,p)\) 对应的归一化面积， \(l_j,l_k\) 类似。

给定 \((x,y)\) ，可求出“面积坐标” \((l_1,l_2,l_2)\) ：

\[ \boxed{\color{red}{\left\{\begin{aligned}l_1&=\frac{1}{2S}\left[(x_2 y_3-x_3 y_2)+(y_2-y_3)x+(x_3-x_2)y\right]\\ l_2&=\frac{1}{2S}\left[(x_3 y_1-x_1 y_3)+(y_3-y_1)x+(x_1-x_3)y\right]\\ l_3&=\frac{1}{2S}\left[(x_1 y_2-x_2 y_1)+(y_1-y_2)x+(x_2-x_1)y\right]\end{aligned} \right.}} \\ 2S=\begin{vmatrix}1 & x_1 & y_1 \\ 1 & x_2 & y_2 \\ 1 & x_3 & y_3 \end{vmatrix} \\ \quad \\ l_1+l_2+l_3=1 \]

这个关系可用下面的Julia代码算出（保存在变量LL中）：

using SymPy

# 解决Colab不显示输出数学公式的问题
using Markdown: MD, LaTeX
function latex(expr)
    expr |> sympy.latex |> LaTeX
end;

@vars x y real=true

# 三个顶点坐标（逆时针排序）
X = [symbols("x$i", real=true) for i in 1:3]
Y = [symbols("y$i", real=true) for i in 1:3]

# 面积坐标
L=[symbols("l$i", positive=true) for i in 1:3]

A = sympy.ones(3,3)
A[:,2]=X
A[:,3]=Y

# 总面积
S = sympy.det(A)/2

# 三个子面积并作归一化
LL = [S(X[i]=>x,Y[i]=>y) for i in 1:3]/S  .|> simplify
LL = LL .|> t->collect(t,[x,y])

# 显示
[L[i]⩵LL[i] for i in 1:3] .|> latex |> MD
\[ l_{1} = \frac{x \left(y_{2} - y_{3}\right) + x_{2} y_{3} - x_{3} y_{2} + y \left(- x_{2} + x_{3}\right)}{x_{1} y_{2} - x_{1} y_{3} - x_{2} y_{1} + x_{2} y_{3} + x_{3} y_{1} - x_{3} y_{2}} \] \[ l_{2} = \frac{x \left(- y_{1} + y_{3}\right) - x_{1} y_{3} + x_{3} y_{1} + y \left(x_{1} - x_{3}\right)}{x_{1} y_{2} - x_{1} y_{3} - x_{2} y_{1} + x_{2} y_{3} + x_{3} y_{1} - x_{3} y_{2}} \] \[ l_{3} = \frac{x \left(y_{1} - y_{2}\right) + x_{1} y_{2} - x_{2} y_{1} + y \left(- x_{1} + x_{2}\right)}{x_{1} y_{2} - x_{1} y_{3} - x_{2} y_{1} + x_{2} y_{3} + x_{3} y_{1} - x_{3} y_{2}} \] 
# 验证归一化条件
sum(LL) .|> simplify 

 1

作为逆变换，给定“面积坐标” \((l_1,l_2,l_3)\) ，可求出普通 \((x,y)\) ：

\[ \color{red}{\boxed{\left\{\begin{aligned}x&=l_1 x_1 + l_2 x_2 + l_3 x_3 \\ y&=l_1 y_1 + l_2 y_2 + l_3 y_3 \end{aligned} \right.}} \]

推导如下（用到了归一化条件，保存在变量XX和YY中）：

# 反过来，可用面积坐标表示(x,y)
eqs = [LL[i]⩵L[i] for i in 1:2]
sol = solve(eqs,[x,y])
(XX,YY) = [x,y] .|> (t->t(sol)) .|>   
        (t->collect(t, [X[3],Y[3]]))  .|> 
        (t->t(L[1]+L[2]=>1-L[3]))

# 显示
[x⩵XX,y⩵YY] .|> latex |> MD
\[ x = l_{1} x_{1} + l_{2} x_{2} + l_{3} x_{3} \] \[ y = l_{1} y_{1} + l_{2} y_{2} + l_{3} y_{3} \]

在面积坐标的表示下，三个顶点坐标分别为(每行代表一个顶点)：

[LL[j](x=>X[i],y=>Y[i]) for i in 1:3, j in 1:3] .|> simplify  |> latex |> MD
\[ \left[\begin{matrix}1 & 0 & 0\\0 & 1 & 0\\0 & 0 & 1\end{matrix}\right] \]

三角元上的多项式插值

\[ p_m(x,y)=\sum^m_{i+j=0}{c_{ij}x^iy^j} \]

注意，这是一个非齐次多项式。 改用面积坐标表示后，原则上是关于 \(l_1,l_2,l_3\) 的非齐次多项式。

但由于 \(1=l_1+l_2+l_3\) ，我们可以在任何位置反复插入 \((l_1+l_2+l_3)\) ，进而可将低次项提升到最高次项（即 \(m\) 次）。

所以总是可改写成关于 \(l_1,l_2,l_3\) 的齐次多项式：

\[ \boxed{p_m(l_1,l_2,l_3)=\sum_{i+j+k=m}{c_{ijk}l_1^i l_2^j l_3^k}} \]

这个多项式和原始多项式都是含有 \(\frac{1}{2}(m+1)(m+2)\) 个待定系数

一次多项式插值

\[ p_1(l_1,l_2,l_3)=c_1 l_1 + c_2 l_2 + c_3 l_3 \]

三个待定系数，需要三个节点值，恰好可用上三角形有三个点的值 \(u_1,u_2,u_3\) 。

据此可解出这6个待定系数,进而可写出用面积坐标表示的一次多项式插值：：

\[ \color{red}{\boxed{p_1(l_1,l_2,l_3)=l_1 u_1 + l_2 u_2 + l_3 u_3}} \]

推导如下（保存在变量pl中，其实这个结果一望而知）：

# 3个节点值
U = [symbols("u$i", real=true) for i in 1:3]
# 3个待定系数
C = [symbols("c$i", real=true) for i in 1:3]

# 3个顶点面积坐标
M= [LL[j](x=>X[i],y=>Y[i]) for i in 1:3, j in 1:3]  .|>  simplify

# 一次多项式插值（有待定系数）
p(r,s,t)=C[1]*r + C[2]*s + C[3]*t

# 3个节点方程
eqs = [ p(M[i,1],M[i,2],M[i,3])⩵U[i] for i in 1:3]

# 求解这个线性方程组
solC = solve(eqs,C)
pl = p(L[1],L[2],L[3])(solC) 

# 显示
[symbols("p1")⩵pl]  .|> latex |> MD
\[ p_{1} = l_{1} u_{1} + l_{2} u_{2} + l_{3} u_{3} \]

二次多项式插值

\[ p_2(l_1,l_2,l_3)=c_1 l_1^2 + c_2 l_2^2 + c_3 l_3^2 + c_4 l_1 l_2 + c_5 l_2 l_3 + c_6 l_3 l_1 \]

这有6个待定系数，除了3个顶点节点外，还需要3个边的中点作为节点。

据此可解出这6个待定系数,进而可写出用面积坐标表示的二次多项式插值：

\[ \boxed{\color{red}{p_2(l_1,l_2,l_3)=\sum^3_{i=1}\left[l_i(2l_i-1)u_i+4l_jl_k u_{3+i}\right]}} \\ l_{3+i}=l_i,\qquad j=i+1,k=i+2 \]

推导如下（用到了归一化条件，保存在变量pl2中）：

# 六个节点值
U = [symbols("u$i", real=true) for i in 1:6]
# 六个待定系数
C = [symbols("c$i", real=true) for i in 1:6]

# 3个顶点面积坐标
M= [LL[j](x=>X[i],y=>Y[i]) for i in 1:3, j in 1:3]  .|>  simplify
# 3个边中点坐标
M = vcat(M,[(M[i%3+1,j]+M[(i+1)%3+1,j])/2 for i in 1:3, j in 1:3]) 

# 二次多项式插值（有待定系数）
p2(r,s,t)=C[1]*r^2 + C[2]*s^2 + C[3]*t^2 + C[4]*r*s + C[5]*s*t + C[6]*t*r

# 6个节点方程
eqs = [ p2(M[i,1],M[i,2],M[i,3])⩵U[i] for i in 1:6]

# 求解这个线性方程组
solC = solve(eqs,C)

# 化简（反复利用了归一化条件）
pl2 = p2(L[1],L[2],L[3])(solC) 
pl2 = collect(pl2 |> expand, U) 
pl2 = sum(pl2.args .|> factor)
pl2 = pl2(L[2]+L[3]=>1-L[1])
pl2 = pl2(L[1]+L[3]=>1-L[2])
pl2 = pl2(L[1]+L[2]=>1-L[3])
pl2 = collect(pl2 |> expand, U) |>  simplify

# 显示
[symbols("p2")⩵pl2]  .|> latex |> MD
\[ p_{2} = 4 l_{1} l_{2} u_{6} + 4 l_{1} l_{3} u_{5} + l_{1} u_{1} \left(2 l_{1} - 1\right) + 4 l_{2} l_{3} u_{4} + l_{2} u_{2} \left(2 l_{2} - 1\right) + l_{3} u_{3} \left(2 l_{3} - 1\right) \]

三次多项式插值

\[ p_3(l_1,l_2,l_3)=c_1 l_1^3 + c_2 l_2^3 + c_3 l_3^3 + c_4 l_1^2 l_2 + c_5 l_1^2 l_3 + c_6 l_2^2 l_1 + c_7 l_2^2 l_3 + c_8 l_3^2 l_1 + c_9 l_3^2 l_2 + c10 l_1 l_2 l_3 \]

这有10个待定系数，1中心点，3个顶点，每个顶点还有2个偏导数，恰好可建10个方程。

由于涉及偏导数，所以要先给出偏导数函数dpx(r,s,t)和dpy(r,s,t)（见后面的代码）。

据此可解出这10个待定系数,进而可写出用面积坐标表示的三次多项式插值：

\[ \boxed{\color{red}{\left\{\begin{aligned}p_3(l_1,l_2,l_3) &= \alpha_0 u_{10} + \sum^3_{i=1}{\left[\alpha_i u_i+\beta_i u_{3+i}+\gamma_i u_{6+i}\right]}\\ \alpha_0 &= 27 l_1 l_2 l_3 \\ \alpha_i &=l_i^3+3 l_i^2(l_j-l_k)-7 l_1,l_2,l_3 \\ \beta_i &= (x_j-x_i)(l_i^2 l_j-l_1 l_2 l_3)+(x_k-x_i)(l_i^2 l_k-l_1 l_2 l_3) \\ \gamma_i &= (y_j-y_i)(l_i^2 l_j-l_1 l_2 l_3)+(y_k-y_i)(l_i^2 l_k-l_1 l_2 l_3) \end{aligned}\right.}} \\ l_{3+i}=l_i,\qquad j=i+1,k=i+2 \]

其中， \(u_{10}\) 对应中点值， \(u_1,u_2,u_3\) 对应3个顶点值， \(u_4,u_5,u_6\) 对应3个顶点沿 \(x\) 偏导数值， \(u_7,u_8,u_9\) 对应3个顶点沿 \(y\) 偏导数值

推导如下（保存在变量pl3中，对应的系数保持在数组αβγ中）：

# 10个节点值
U = [symbols("u$i", real=true) for i in 1:10]
# 10个待定系数
C = [symbols("c$i", real=true) for i in 1:10]

# 三次多项式插值（有待定系数）
p3(r,s,t)=C[1]*r^3+C[2]*s^3+C[3]*t^3+C[4]*r^2*s+C[5]*r^2*t+C[6]*s^2*r+C[7]*s^2*t+C[8]*t^2*r+C[9]*t^2*s+C[10]*r*s*t
# 偏导数函数
ex = p3(L[1],L[2],L[3])
ex1 = sum(diff(ex,L[i])*diff(LL[i],x) for i in 1:3) |>  simplify
ex2 = sum(diff(ex,L[i])*diff(LL[i],y) for i in 1:3) |>  simplify
dpx(r,s,t)=ex1(L[1]=>r,L[2]=>s,L[3]=>t)
dpy(r,s,t)=ex2(L[1]=>r,L[2]=>s,L[3]=>t);

# 3个顶点面积坐标
M= [LL[j](x=>X[i],y=>Y[i]) for i in 1:3, j in 1:3] .|>  simplify

# 3个顶点方程
eqs = [ p3(M[i,1],M[i,2],M[i,3])⩵U[i] for i in 1:3]

# 6个偏导数方程
eqs = vcat(eqs,[dpx(M[i,1],M[i,2],M[i,3])⩵U[3+i] for i in 1:3])
eqs = vcat(eqs,[dpy(M[i,1],M[i,2],M[i,3])⩵U[6+i] for i in 1:3])

# 1个中点方程
eqs = vcat(eqs,[p3(sum(M[:,1])/3,sum(M[:,2])/3,sum(M[:,3])/3)⩵U[10]])

# 解出待定系数
@time solC = solve(eqs,C)

# 三次多项式插值（结果）
pl3 = p3(L[1],L[2],L[3])(solC)
pl3 = collect(pl3 |> expand, U) 

# 列出u1,...,u10的系数
# αβγ[10]对应α0
# αβγ[1:3]对应α1,α2,α3
# αβγ[4:6]对应β1,β2,β3
# αβγ[7:9]对应γ1,γ2,γ3
αβγ = [diff(pl3,U[i]) for i in 1:10]

# 尽可能对每个系数合并同类项
ls = [L[1]^i*L[2]^j*L[3]^k for i in 0:3, j in 0:3, k in 0:3 if i+j+k==3]
αβγ = αβγ .|>  t->collect(t, ls) 

# 显示
tmp = [symbols("α0")⩵αβγ[10]] 
tmp = vcat(tmp,[symbols("α$i")⩵αβγ[i] for i in 1:3])
tmp = vcat(tmp,[symbols("β$i")⩵αβγ[3+i] for i in 1:3])
tmp = vcat(tmp,[symbols("γ$i")⩵αβγ[6+i] for i in 1:3])
tmp .|> latex |> MD

 15.244530 seconds (597 allocations: 18.047 KiB)
\[ α0 = 27 l_{1} l_{2} l_{3} \] \[ α1 = l_{1}^{3} + 3 l_{1}^{2} l_{2} + 3 l_{1}^{2} l_{3} - 7 l_{1} l_{2} l_{3} \] \[ α2 = 3 l_{1} l_{2}^{2} - 7 l_{1} l_{2} l_{3} + l_{2}^{3} + 3 l_{2}^{2} l_{3} \] \[ α3 = - 7 l_{1} l_{2} l_{3} + 3 l_{1} l_{3}^{2} + 3 l_{2} l_{3}^{2} + l_{3}^{3} \] \[ β1 = l_{1}^{2} l_{2} \left(- x_{1} + x_{2}\right) + l_{1}^{2} l_{3} \left(- x_{1} + x_{3}\right) + l_{1} l_{2} l_{3} \left(2 x_{1} - x_{2} - x_{3}\right) \] \[ β2 = l_{1} l_{2}^{2} \left(x_{1} - x_{2}\right) + l_{1} l_{2} l_{3} \left(- x_{1} + 2 x_{2} - x_{3}\right) + l_{2}^{2} l_{3} \left(- x_{2} + x_{3}\right) \] \[ β3 = l_{1} l_{2} l_{3} \left(- x_{1} - x_{2} + 2 x_{3}\right) + l_{1} l_{3}^{2} \left(x_{1} - x_{3}\right) + l_{2} l_{3}^{2} \left(x_{2} - x_{3}\right) \] \[ γ1 = l_{1}^{2} l_{2} \left(- y_{1} + y_{2}\right) + l_{1}^{2} l_{3} \left(- y_{1} + y_{3}\right) + l_{1} l_{2} l_{3} \left(2 y_{1} - y_{2} - y_{3}\right) \] \[ γ2 = l_{1} l_{2}^{2} \left(y_{1} - y_{2}\right) + l_{1} l_{2} l_{3} \left(- y_{1} + 2 y_{2} - y_{3}\right) + l_{2}^{2} l_{3} \left(- y_{2} + y_{3}\right) \] \[ γ3 = l_{1} l_{2} l_{3} \left(- y_{1} - y_{2} + 2 y_{3}\right) + l_{1} l_{3}^{2} \left(y_{1} - y_{3}\right) + l_{2} l_{3}^{2} \left(y_{2} - y_{3}\right) \]