有限元方法简单的二维算例（三角形剖分）

文章目录

有限元方法简单的二维算例（三角形剖分）
算例描述
变分问题
有限元离散
- 问题转化
- 有限元三要素
参考单元与一般单元
- 一般单元上的形函数
- 一般单元上的积分计算
问题求解
- 局部基函数到全局基函数
- 单元扫描计算刚度矩阵和载荷向量
- 边界条件处理
数值实验和收敛阶
- 程序代码
- 误差和收敛阶分析

算例描述

我们对下述椭圆边值问题 \label{eq1}（破编辑器不容易给公式编号，只能这样写）{−Δu=fu∣∂Ω=0\left \{ \begin{aligned} & -\Delta u = f\\ & u|_{\partial \Omega} = 0 & \end{aligned} \right.{−Δu=fu∣∂Ω=0
其中，Ω=(0,1)2\Omega = (0,1)^2Ω=(0,1)2且f=2π2sinπxsinπyf=2\pi^2\text{sin}\pi x\text{sin}\pi yf=2π2sinπxsinπy。考虑其变分问题，对其变分问题有限元离散并求解，并验证其为一阶收敛。
注：问题的真解为u(x)=sinπxsinπyu(x) = \text{sin}\pi x\text{sin}\pi yu(x)=sinπxsinπy。

变分问题

∀v∈H01\forall v \in H_0^1∀v∈H01，乘([eq1])式两边，使用格林公式，利用边界条件，易得：\label{eq2}
∫Ω∇u∇vdx=∫Ωfvdx\int _\Omega \nabla u \nabla v dx = \int _\Omega fv dx∫Ω∇u∇vdx=∫Ωfvdx
其中fff为方程中的右端项。令\label{eq3}

a(u,v)=∫Ω∇u∇vdxf(v)=∫Ωfvdx\begin{aligned} &a(u,v) = \int _\Omega \nabla u \nabla v dx & \\ &f(v) = \int _\Omega fv dx & \end{aligned}a(u,v)=∫Ω∇u∇vdxf(v)=∫Ωfvdx

容易证明，原问题([eq1])等价于变分问题: \label{eq4}
{求u∈H01，使得a(u,v)=f(v)∀v∈H01\left \{ \begin{aligned} & \text{求}u \in H_0^1 \text{，使得}\\ & a(u,v) = f(v)&\forall v \in H_0^1 \end{aligned} \right. {求u∈H01，使得a(u,v)=f(v)∀v∈H01
事实上，在一定连续性的要求下，强解为弱解，弱解也是强解，二者等价。故求解问题([eq1])变为了求解问题([eq4])。更一般的变分问题，描述为：\label{eq5}
{求u∈V,s.t.a(u,v)=f(v)∀v∈V\left \{ \begin{aligned} & \text{求}u \in V ,s.t.&\\ & a(u,v) = f(v)&\forall v \in V& \end{aligned} \right.{求u∈V,s.t.a(u,v)=f(v)∀v∈V

有限元离散

问题转化

我们来考虑上述变分问题的有限维逼近，即构造VVV的有限维子空间Vh⊂VV_h \subset VVh⊂V，考虑如下的离散问题：\label{eq6}
{求uh∈Vh,s.t.a(uh,vh)=f(vh)∀vh∈Vh\left \{ \begin{aligned} & \text{求}u_h \in V_h ,s.t.&\\ & a(u_h,v_h) = f(v_h)&\forall v_h \in V_h& \end{aligned} \right.{求uh∈Vh,s.t.a(uh,vh)=f(vh)∀vh∈Vh

我们用问题([eq6])近似问题([eq5])，后者的解逼近前者的解。所以我们可以通过求([eq6])的解作为近似解。
设VhV_hVh空间中的一组基为{ϕi}i=1N\{\phi_i\}_{i=1}^N{ϕi}i=1N，若uh=Σuiϕiu_h = \Sigma u_i \phi_iuh=Σuiϕi（为了书写方便，不加说明，求和指标都为1到N），我们需要求的是组合系数uiu_iui，将uh=Σuiϕiu_h = \Sigma u_i \phi_iuh=Σuiϕi代入，并依次分别取vhv_hvh 为每个基函数，我们可以得到： \label{eq7}Σuia(ϕi,ϕj)=f(ϕj),j=1,2,...N.\Sigma u_i a(\phi_i,\phi_j)=f(\phi_j),j=1,2,...N.Σuia(ϕi,ϕj)=f(ϕj),j=1,2,...N.
这事实上就是一个线性方程组AU=FAU=FAU=F，其中A=(a(ϕj,ϕi)i,j=1N),U=(u1,...,uN)T,F=(f(ϕ1,...,ϕN))TA = (a(\phi_j,\phi_i)_{i,j=1}^N),U=(u_1,...,u_N)^T,F=(f(\phi_1,...,\phi_N))^TA=(a(ϕj,ϕi)i,j=1N),U=(u1,...,uN)T,F=(f(ϕ1,...,ϕN))T。那么，求解问题就转为了在VhV_hVh的约束下，求解这个线性方程组。

由我们选取vhv_hvh分别为ϕi\phi_iϕi可知，这个方程组的解对于原问题是必要而不充分的。但是在合适的条件之下，由Lax-Milgram定理可知，离散变分问题的解是存在唯一的。一般来说，我们在初边值条件的约束下，方程组的解存在唯一，那么这个唯一解必然是原问题的解。故而，对于一般的变分问题，我们通过离散化，最后可以通过求解([eq7])来近似。那么问题本质上就转化为了找一个合适的空间逼近，在这个空间中求基函数，和基函数之间的某种相互作用以及基函数和fff之间的作用，构建出方程组，最后求解方程组，得到逼近阶关于基函数的组合系数，最后得到原问题的一个逼近的数值解。

需要注意的是，由于Vh∈VV_h\in VVh∈V，我们对ϕi\phi_iϕi是有一定的要求的，或者说我们不管基函数是否属于VVV,我们最后对组合系数uiu_iui施以一定的要求（为满足边界条件），使得uh=Σuiϕi∈Vu_h = \Sigma u_i \phi_i \in Vuh=Σuiϕi∈V，具有相同的效力。

有限元三要素

我们现在来考虑单元上的插值问题。对于问题([eq1])，我们考虑其三要素(T,PT,ΣT)(T,P_T,\Sigma_T)(T,PT,ΣT)，其中，TTT为三角形分割，如图[pic1]所示，为了方便，我这里假设的是两个方向上是等距分割，也就是三角形是等边三角形。

我们先来考虑如图所示的参考单元，即采用面积坐标。

分割单元上的形函数空间为多项式集合PT=span{λ1,λ2,λ3}P_T = \text{span}\{\lambda_1,\lambda_2,\lambda_3\}PT=span{λ1,λ2,λ3}，其中λ1+λ2+λ3=1\lambda_1+ \lambda_2+ \lambda_3=1λ1+λ2+λ3=1。单元自由度ΣT={u(a1),u(a2),u(a3)}\Sigma_T = \{u(a_1),u(a_2),u(a_3)\}ΣT={u(a1),u(a2),u(a3)}
表示三个端点的值。(注意，自由度是个集合，而不是个数，很多人在口头表述的时候，总喜欢表述“自由度是几”，这个本身容易引起误导，也侧面反应了这些人基础概念的混淆。)

适定性
显然，若Σ={0}\Sigma = \{0\}Σ={0}，那么u≡0u\equiv0u≡0。
求基函数
显然，三个面积坐标λ1,λ2,λ3\lambda_1,\lambda_2,\lambda_3λ1,λ2,λ3即为单元上的形函数。
连续性（连续才能保证Vh∈H1V_h \in H^1Vh∈H1，只有这样才能和前面的基础理论相容。）
考虑到λ1+λ2+λ3=1\lambda_1+ \lambda_2+ \lambda_3=1λ1+λ2+λ3=1，相邻插值函数在共用边上都为某个变量的一次函数，两端点值固定，一次函数也就确定了，故而，在边界处是连续的。

参考单元与一般单元

一般单元上的形函数

我们知道在参考单元上的一组基函数为λ1,λ2,λ3\lambda_1,\lambda_2,\lambda_3λ1,λ2,λ3，那么一般单元上的形函数是什么呢？

首先我们注意到，两者转换满足以下公式： x=∑i=13xiλi=x3+λ1(x1−x3)+λ2(x2−x3)y=∑i=13yiλi=y3+λ1(y1−y3)+λ2(y2−y3)\begin{aligned}x = \sum\limits_{i=1}^3 x_i\lambda_i=x_3+\lambda_1(x_1-x_3)+\lambda_2(x_2-x_3) \\ y = \sum\limits_{i=1}^3 y_i\lambda_i=y_3+\lambda_1(y_1-y_3)+\lambda_2(y_2-y_3)\end{aligned}x=i=1∑3xiλi=x3+λ1(x1−x3)+λ2(x2−x3)y=i=1∑3yiλi=y3+λ1(y1−y3)+λ2(y2−y3)

反过来，令三角形节点标号为逆时针，记为i,j,ki,j,ki,j,k，且引入记号ξi=xj−xk,ηi=yj−yk,ωi=xjyk−xkyj\xi_i = x_j -x_k,\eta_i = y_j-y_k,\omega_i = x_jy_k-x_ky_jξi=xj−xk,ηi=yj−yk,ωi=xjyk−xkyj时，面积坐标可表示为（很容易 check）：
λi=ηix−ξiy+ωi2∣T∣,i=1,2,3\lambda_i = \frac{\eta_ix-\xi_iy+\omega_i}{2|T|},i=1,2,3λi=2∣T∣ηix−ξiy+ωi,i=1,2,3

这里的∣T∣|T|∣T∣表示三角形面积，这里的i,j,ki,j,ki,j,k标号一定要以逆时针顺序编号，且ξi\xi_iξi和ηi\eta_iηi的表达式为前面一个坐标减去后面一个坐标，否则λi\lambda_iλi的表达式会差一个符号，所以，为了不出错，我们提倡节点编号都以逆时针方向来编。

由此，我们得到了λi\lambda_iλi到一般单元上的映射，即知道了一般单元上的形函数。

一般单元上的积分计算

我们知道，搞清楚了参考单元上的情况，并不意味着一般单元上的情况就清楚了。我们利用需要两者之间的转换和变量替换关系，搞清楚一般单元上的积分和参考单元上的积分的关系。
放到参考单元上去做，表达起来会更加简洁有力美观，同时，在计算多个单元上的积分时，一些中间量如雅克比行列式和雅克比矩阵等可被复用，一定程度上减少了计算量。当然，你不做 scaling，直接在原单元上求基函数，及其相互之间以及右端项之间的作用的积分，也是可以的，表达式复杂一点点而已。

由此，我们很容易求得雅克比矩阵（红色部分修正，对雅克比矩阵的定义做了统一，避免前后不一致）。
J=∂(x,y)∂(λ1,λ2)=[−ξ2,ξ1;−η2,η1]TJ = \frac{\partial(x,y)}{\partial(\lambda_1,\lambda_2)}=[ -\xi_2,\xi_1;-\eta_2,\eta_1 ]^TJ=∂(λ1,λ2)∂(x,y)=[−ξ2,ξ1;−η2,η1]T
J=∂(x,y)∂(λ1,λ2)=[xλ1,xλ2;yλ1,yλ2]=[−ξ2,ξ1;−η2,η1]（修正后）J = \frac{\partial(x,y)}{\partial(\lambda_1,\lambda_2)}=[x_{\lambda_1},x_{\lambda_2}; y_{\lambda_1},y_{\lambda_2}]=[ -\xi_2,\xi_1;-\eta_2,\eta_1 ]（修正后）J=∂(λ1,λ2)∂(x,y)=[xλ1,xλ2;yλ1,yλ2]=[−ξ2,ξ1;−η2,η1]（修正后）

雅克比行列式∣J∣=η2ξ1−η1ξ2|J| = \eta_2\xi_1 - \eta_1\xi_2∣J∣=η2ξ1−η1ξ2。

另外，我们也需要雅克比矩阵的逆J−1=∂(λ1,λ2)∂(x,y)=[η1/(2T),−ξ1/(2T);η2/(2T),−ξ2/(2T)]J^{-1} = \frac{\partial(\lambda_1,\lambda_2)}{\partial(x,y)}= [\eta_1/(2T), -\xi_1/(2T);\eta_2/(2T), -\xi_2/(2T)] J−1=∂(x,y)∂(λ1,λ2)=[η1/(2T),−ξ1/(2T);η2/(2T),−ξ2/(2T)]。
（这里少了个转置符号，需要加上，具体看评论）

那么，在一般单元上的积分计算就可以和参考单元上的积分计算联系上，如果不涉及求导，那么之差雅克比行列式，涉及求导，两者的积分就不仅仅是差一个常数，如下所示：
∫T∇u(x,y)T∇v(x,y)dxdy=∫T^∇u(λ1,λ2)TJ−TJ−1∇v(λ1,λ2)∣J∣dλ1dλ2\int_T {\nabla u(x,y)}^T\nabla v(x,y) dxdy = \int _{\hat T} {\nabla u(\lambda_1,\lambda_2)}^T J^{-T}J^{-1}\nabla v(\lambda_1,\lambda_2)|J|d\lambda_1d \lambda_2 ∫T∇u(x,y)T∇v(x,y)dxdy=∫T^∇u(λ1,λ2)TJ−TJ−1∇v(λ1,λ2)∣J∣dλ1dλ2
∫T∇u(x,y)T∇v(x,y)dxdy=∫T^∇u(λ1,λ2)TJ−1J−T∇v(λ1,λ2)∣J∣dλ1dλ2（修正后）\int_T {\nabla u(x,y)}^T\nabla v(x,y) dxdy = \int _{\hat T} {\nabla u(\lambda_1,\lambda_2)}^T J^{-1}J^{-T}\nabla v(\lambda_1,\lambda_2)|J|d\lambda_1d \lambda_2 （修正后）∫T∇u(x,y)T∇v(x,y)dxdy=∫T^∇u(λ1,λ2)TJ−1J−T∇v(λ1,λ2)∣J∣dλ1dλ2（修正后）

∫Tf(x,y)v(x,y)dxdy=∫T^f(λ1,λ2)v(λ1,λ2)∣J∣dλ1dλ2\int_T f(x,y)v(x,y) dxdy = \int _{\hat T} f(\lambda_1,\lambda_2)v(\lambda_1,\lambda_2)|J|d\lambda_1 d\lambda_2∫Tf(x,y)v(x,y)dxdy=∫T^f(λ1,λ2)v(λ1,λ2)∣J∣dλ1dλ2

问题求解

局部基函数到全局基函数

求得了在每个单元上的插值的节点基函数，对于每个节点，我们将其相关单元的相关基函数并成一个关于这个节点的全局基函数ϕi\phi_iϕi，有多少个节点，就有多少个基函数。

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

有了基函数，理论上就可以求解变分问题离散后转化成的方程组。我们需要计算a(ϕi,ϕj)a(\phi_i,\phi_j)a(ϕi,ϕj)和f(ϕj)f(\phi_j)f(ϕj)。
这里对于积分的计算(a(ϕi,ϕj)a(\phi_i,\phi_j)a(ϕi,ϕj)和f(Φj)f(\Phi_j)f(Φj))可以采用扫描单元的方式，即在每个单元上计算基函数之间的相互作用并分配到相对应的全局基函数的相互作用上，计算每个单元上的基函数和右端项的作用，分配到与之相关的全局基函数和右端项的作用上。我们就需要计算单元上基函数作用a(fi,fj)a(f_i,f_j)a(fi,fj)和单元上基函数和右端项的作用f(fi)f(f_i)f(fi)。这里为了区分局部基函数和全局基函数，我用fif_ifi表示单元上的基函数，ϕi\phi_iϕi表示全局基函数。

由前可知，我们事实上已经能够通过将一般单元变换到参考单元来计算一般三角形单元上两个形函数之间的相互作用
a(fi,fj)=∫T∇fi(x,y)T∇fj(x,y)dxdy=∫T^∇λ1TJ−TJ−1∇λ2∣J∣dλ1dλ2a(f_i,f_j) = \int_T {\nabla f_i(x,y)}^T\nabla f_j(x,y) dxdy = \int _{\hat T} {\nabla \lambda_1}^T J^{-T}J^{-1}\nabla \lambda_2|J|d\lambda_1d \lambda_2a(fi,fj)=∫T∇fi(x,y)T∇fj(x,y)dxdy=∫T^∇λ1TJ−TJ−1∇λ2∣J∣dλ1dλ2

a(fi,fj)=∫T∇fi(x,y)T∇fj(x,y)dxdy=∫T^∇λ1TJ−1J−T∇λ2∣J∣dλ1dλ2（修正后）a(f_i,f_j) = \int_T {\nabla f_i(x,y)}^T\nabla f_j(x,y) dxdy = \int _{\hat T} {\nabla \lambda_1}^T J^{-1}J^{-T}\nabla \lambda_2|J|d\lambda_1d \lambda_2（修正后）a(fi,fj)=∫T∇fi(x,y)T∇fj(x,y)dxdy=∫T^∇λ1TJ−1J−T∇λ2∣J∣dλ1dλ2（修正后）

以及单元上fff和基函数之间的作用f(fj)=∫Tf(x,y)fj(x,y)dxdy=∫T^f(λ1,λ2)λj∣J∣dλ1dλ2f(f_j)=\int_T f(x,y)f_j(x,y) dxdy = \int _{\hat T} f(\lambda_1,\lambda_2)\lambda_j|J|d\lambda_1 d\lambda_2f(fj)=∫Tf(x,y)fj(x,y)dxdy=∫T^f(λ1,λ2)λj∣J∣dλ1dλ2。

对于f(fj)=∫T^f(λ1,λ2)λj∣J∣dλ1dλ2f(f_j) = \int _{\hat T} f(\lambda_1,\lambda_2)\lambda_j|J|d\lambda_1 d\lambda_2f(fj)=∫T^f(λ1,λ2)λj∣J∣dλ1dλ2的计算，由于在一般问题中，fff不一定就是这个问题的fff，所以，对于具体问题的具体的fff我们可以用 MATLAB 内置的求积分函数，也可以自己编写程序实现。

边界条件处理

但是在边值条件问题中，粗算出来的方程组的刚度矩阵KKK，就不是满秩的，那么解就不唯一。这是因为所求解的uh∈Vh⊂Vu_h \in V_h \subset Vuh∈Vh⊂V这个条件没满足，所以，一个思路是根据边界条件选取一个合适的基函数子集，另外，也可以最后根据边值条件来调整刚度矩阵KKK，我采取第二个思路。

调整的思路就是调整KU=FKU=FKU=F，将其等价转换，使得最后的解满足UUU的对应位置为固定的值，所谓对应位置，就是边界所对应的位置，所谓固定的值，就是边界值。简单地说，我们可以这样操作，直接令UUU的对应位置强制为边界值，然后将UUU中新的不知道的数视为未知数，构建新的方程组求解，本质上无非就是右端项减去边界条件乘以KKK中与其位置对应的列……假若，原来的KKK秩比起满秩少了m，那么补上m个边界条件，方程组的解应该是存在唯一的，认识到这一点，再加上一些线性代数的基本知识，事情就很明了了……就不再赘述……

数值实验和收敛阶

程序代码

基于以上的思想，编写了程序代码，直接粘贴如下，为了方便，直接将函数文件直接粘贴在主文件末尾了。程序的细节可以看注释，就不再详述。我基本上将每一句代码都打上了注释。

%% 这是一个二维的有限元程序
%% 这是一个二维的有限元程序
clc;  % 清空命令行窗口
clear; %清除工作空间
close all; %关闭所有图像
%% 参数设置
error_L2s = [];
for k = [2,4,8,16,32,64]
Lx = 1; %定义单元右边界（左边界为0，如果不是，可以平移为0）
Ly = 1;%定义单元上边界
N = k;%分割的一个方向的单元数目
numelx = N;%定义分割的x方向单元数目（按矩形计算）
numely = N;%定义分割的y方向单元数目（按矩形计算）
hx = Lx/numelx;%x方向上的单元长度（对应着三角形两条直角边之一）
hy = Ly/numely;%y方向上的单元长度（对应着三角形两条直角边之一）
numel = numelx*numely*2;%小单元的数目,每个矩形分成两个三角形
u_b = zeros(2*(numelx+numely),1); %定义第一类边界条件，一圈过来都是0
numnodx = numelx + 1; % x方向节点个数比单元个数多1
numnody = numely + 1; % y方向节点个数比单元个数多1
numnod = numnodx*numnody;%总的节点个数
nel = 3;%每个单元的节点数目,即每个单元上有几个形函数参与作用，单元自由度
coordx = linspace(0,Lx,numnodx)'; %等分节点的坐标（为了方便，我这里采用等分的方式，事实上单元长度可以不一致，非均匀网格）
coordy = linspace(0,Ly,numnody)'; %等分节点的坐标（为了方便，我这里采用等分的方式，事实上单元长度可以不一致）
[X, Y] = meshgrid(coordx,coordy);%张成网格，X和Y分别表示对应位置的横纵坐标
X = X';Y = Y';coord = [X(:) Y(:)];%把网格一行一行扯开，coord的每一行是对应节点的坐标，按顺序排
connect = connect_mat(numnodx,numnody,nel);%连接矩阵，表示每个单元周围的节点编号，也就是涉及的形函数编号
ebcdof = unique([1:numnodx numnodx*numely+1:numnodx*numely+numnodx ...numnodx+1:numnodx:numnodx*(numely-1)+1 2*numnodx:numnodx:numely*numnodx]); % 强制性边界点的编号，本例子中是四条边，下上左右边
ebcval = u_b; %假设边界值都为u_b
bigk = sparse(numnod,numnod); % 刚度矩阵[K]，初始化为0，使用稀疏矩阵存储
fext = sparse(numnod,1);      % 载荷向量{f},初始化为0%%  计算系数矩阵K和右端项f，单刚组装总刚
for e = 1:numel %同一维的情况，依然按单元来扫描ke = elemstiff2d(e,nel,hx,hy,coord,connect);%计算单元刚度矩阵fe = elemforce2d(e,nel,hx,hy,coord,connect);%计算单元载荷向量sctr = connect(e,:);bigk(sctr,sctr) = bigk(sctr,sctr) + ke;fext(sctr) = fext(sctr) + fe;%full(bigk)%full(fext)
end
for i = 1:length(ebcdof)n = ebcdof(i);for j = 1:numnodif (isempty(find(ebcdof == j, 1))) % 第j个点若不是固定点fext(j) = fext(j) - bigk(j,n)*ebcval(i);endendbigk(n,:) = 0.0;bigk(:,n) = 0.0;bigk(n,n) = 1.0;fext(n) = ebcval(i);
end
u_coeff = bigk\fext;%求出系数，事实上也是函数在对应点上的值
u_cal = u_coeff;
u_cal_re = reshape(u_coeff,numnodx,numnody);
u_cal_re = full(u_cal_re);%% 求精确解
L =Lx;
nsamp = 1001;
xsamp = linspace(0,L,nsamp);%100等分区间中间有100个数
[X,Y] = meshgrid(xsamp,xsamp);
uexact = exactsolution2d(X(:),Y(:));
uexact_re = reshape(uexact,nsamp,nsamp);
mesh(xsamp,xsamp,uexact_re)
%% 绘图，可视化
h = mesh(coordx,coordy,u_cal_re);
title(' FE Solutions');%标题
saveas(h,['报告\pics\k' num2str(k) '.eps'],'epsc2')
%% 求精确解，计算误差，没有计算单元准确积分，依然使用近似误差u_ex = exactsolution2d(coord(:,1),coord(:,2));
u_ex_re = reshape(u_ex,numnodx,numnody);
error_L1 = sum(abs(u_cal - u_ex))*hx*hy
error_L2 = sqrt(sum((u_cal - u_ex).^2)*hx*hy)
error_Linf = max(abs(u_cal - u_ex))
error_L2s(end+1) = error_L2;end
es = log2(error_L2s(1:end-1)./error_L2s(2:end))function [ke] = elemstiff2d(e,nel,hx,hy,coord,connect)
T = hx*hy/2;
ke    = zeros(nel,nel); %单刚初始化
nodes = connect(e,:);%相关形函数（节点）编号
xe    = coord(nodes,:); %相关节点的坐标
xi1 = xe(2,1) - xe(3,1);xi2 = xe(3,1) - xe(1,1);%xi3 = xe(1,1) - xe(2,1);
eta1 = xe(2,2) - xe(3,2);eta2 = xe(3,2) - xe(1,2);%eta3 = xe(1,2) - xe(2,2);
%以下内容可修改，具体看评论
auv11 = -((eta1*xi2 - eta2*xi1)*(eta1^2 + eta2^2))/(8*T^2);
auv12 = ((eta1*xi1 + eta2*xi2)*(eta1*xi2 - eta2*xi1))/(8*T^2);
auv13 = -((eta1*xi2 - eta2*xi1)*(- eta1^2 + xi1*eta1 - eta2^2 + xi2*eta2))/(8*T^2);
auv22 = -((xi1^2 + xi2^2)*(eta1*xi2 - eta2*xi1))/(8*T^2);
auv23 = -((eta1*xi2 - eta2*xi1)*(- xi1^2 + eta1*xi1 - xi2^2 + eta2*xi2))/(8*T^2);
auv33 = -((eta1*xi2 - eta2*xi1)*(eta1^2 - 2*eta1*xi1 + eta2^2 - 2*eta2*xi2 + xi1^2 + xi2^2))/(8*T^2);
ke = ke + [auv11 auv12 auv13;auv12 auv22 auv23;auv13 auv23 auv33];
return
endfunction [fe] = elemforce2d(e,nel,hx,hy,coord,connect)
%输入单元编号e，单元自由度nel数目，单元长度hx，单元宽度hy，单元节点坐标coord，和连接矩阵connect
%输出单元载荷fe
%考虑以一般单元为基准，可求f在这个一般单元上的取值，再在这个一般单元上积分
%比较规范的划分，求积分不妨直接在一般单元上求，没必要变换到标准单元上，一来若求梯度更加麻烦了，标准单元上的依然需要求不同的积分
%二来因为f的值未定，依然需要求积分
%fe = zeros(nel,1); %初始化载荷向量
nodes = connect(e,:);%自由度编号
xe = coord(nodes,:); % 单元自由节点坐标
xi1 = xe(2,1) - xe(3,1);xi2 = xe(3,1) - xe(1,1);%xi3 = xe(1,1) - xe(2,1);
eta1 = xe(2,2) - xe(3,2);eta2 = xe(3,2) - xe(1,2);%eta3 = xe(1,2) - xe(2,2);
%syms lam1 lam2;%x =  xe(3,1) + lam1*(-xi2)+lam2*(xi1);
%y =  xe(3,2) + lam1*(-eta2)+lam2*(eta1);detJ = eta2*xi1 - eta1*xi2;
g1 = @(lam1,lam2) fun(xe(3,1) + lam1*(xi2)+lam2*(xi1),xe(3,2) + lam1*(-eta2)+lam2*(eta1)).*lam1*detJ;%g1 = matlabFunction(g1);
g2 = @(lam1,lam2) fun(xe(3,1) + lam1*(xi2)+lam2*(xi1),xe(3,2) + lam1*(-eta2)+lam2*(eta1)).*lam2*detJ;%g2 = matlabFunction(g2);
g3 = @(lam1,lam2) fun(xe(3,1) + lam1*(xi2)+lam2*(xi1),xe(3,2) + lam1*(-eta2)+lam2*(eta1)).*(1-lam1-lam2)*detJ;%g3 = matlabFunction(g3);lammax = @(lam1) 1 - lam1;
gx(1) = integral2(g1,0,1,0,lammax);
gx(2) = integral2(g2,0,1,0,lammax);
gx(3) = integral2(g3,0,1,0,lammax);
fe = [gx(1);gx(2);gx(3)];
endfunction bx = fun(x,y)
bx = (2*pi^2)*sin(pi*x).*sin(pi*y);
endfunction connect_mat = connect_mat( numnodx,numnody,nel)
%输入横纵坐标的节点数目，和单元自由度
%输出连接矩阵，每个单元涉及的节点的编号
xn = 1:(numnodx*numnody);%拉成一条编号
A = reshape(xn,numnodx,numnody);%同形状编号
for i = 1:(numnodx-1)*(numnody-1)xg = rem(i,numnodx-1);%xg表示单元为左边界数起第几个if xg == 0xg = numnodx-1;endyg = ceil(i/(numnodx-1));%下边界其数第几个a = A(xg:xg+1,yg:yg+1);%这个小矩阵，拉直了就是连接矩阵a_vec = a(:);connect_mat(2*i-1:2*i,1:nel) = [a_vec([1 4 3])';a_vec([4 1 2])'];
end
endfunction uexact = exactsolution2d(xg,yg)
uexact = sin(pi*xg).*sin(pi*yg);return
end

美中不足的是，代码中求解方程组部分，我直接使用matlab的右除算子\\backslash\，对于这种稀疏矩阵，其实可以考虑使用一些快速的迭代方法。我这里单元内部的节点编号，使用了一个连接矩阵来存储，即按照一定的顺序将每个单元上的节点全局编号按行存到连接矩阵中，需要的时候再取出来。

误差和收敛阶分析

对于这个二维问题，精确解和数值解不好画在一张图上做对比。我们可以看到，该问题的精确解如图[pic3]所示。

定义单元长度为1，不妨假设两个方向分割的单元数目相同，记为k，结果如下图所示。收敛阶的计算是基于L2L_2L2误差的，使用其它度量结果差不多。分别设置k的不同，使用程序求解，最后的结果如图下图所示。

计算其收敛阶，结果如下所示。

由此可见，二维线性三角形元的在LpL_pLp空间的意义下收敛阶是2。

有限元方法入门：有限元方法简单的二维算例（三角形剖分）相关推荐

有限元方法入门：有限元方法简单的二维算例（矩形剖分）
#有限元方法简单的二维算例(矩形剖分) 算例描述我们对下述椭圆边值问题 \label{eq1}{−Δu=fu∣∂Ω=0\left \{ \begin{aligned} & -\Delta u ...
autocad2007二维图画法_cad怎样绘制简单的二维图形
CAD绘制二维图形非常的简单,大家经常用它来画图,下面是学习啦小编带来关于cad怎样绘制简单的二维图形的内容,希望可以让大家有所收获! cad绘制简单二维图形的方法 1.绘图菜单绘图菜单是绘制图形最基 ...
如何用Python制作一个简单的二维码生成器
目录前言 1.安装第三方库 2.QRCode参数解释 3.自定义二维码生成器 4.给二维码加图片 5.全部代码 6.结果前言二维码又称二维条码,常见的二维码为QR Code,QR全称Quick ...
使用C++实现简单的二维细胞自动机
使用C++实现简单的二维细胞自动机细胞自动机(cellular automata)是为模拟包括自组织结构在内的复杂现象提供的一个强有力的方法,也称为元胞自动机(Cellular Automaton) ...
如何使用JAVA代码生成一个简单的二维码
提示:文章写完后,目录可以自动生成,如何生成可参考右边的帮助文档文章目录前言一.二维码是什么? 二.使用步骤 1.引入依赖 2.开始操作 3.二维码容错 4.结果展示总结前言二维码大家应该 ...
Android 简单生成二维码名片
二维码名片是现在很常见的,这里只是一个简单生成二维码,如果对二维码名片的内容没有过多的要求,可以借鉴一下.生成二维码用的是谷歌的Zxing库,关于扫码,这里就不多说了. zxing下载地址:http: ...
【Matlab】一种超简单的二维矩阵降维方法
1.Introduction Matlab里图像处理时,经常会把一维数组转二维数组,二维数组转一维,如下图所示: 一般经常使用的函数是 reshape ,可以在不同维度之间进行转换,不过需要事先计算数 ...
二维码制作方法有哪些？教你简单的二维码制作方法
二维码是怎么制作的呢?二维码是用某种特定的几何图形按照一定规律在平面(二维方向)分布的黑白相间的图形记录数据符号信息的.现如今,随着智能手机的广泛普及和技术的不断改进,二维码已经被广泛应用于商业领域中 ...
php 数组降维,php 数组去重的方法参考(一维数组去重、二维数组去重)
本文介绍下,对php数组去除重复的方法,包括一维数组的去重.二维数组的去重.有需要的朋友参考下. 首先,来看一维数组重复项的去除方法. 使用array_unique函数,例如: 输出结果: Array ...

有限元方法入门：有限元方法简单的二维算例（三角形剖分）