发布时间 : 星期五 文章matlaB程序的有限元法解泊松方程更新完毕开始阅读
[X,Y,NN,NE]=setelm_tri(Imax,Jmax); % 给节点和三角形元素编号，并设定节点坐标 % 以下求解有限元方程的求系数矩阵 T=zeros(ndm,ndm); for n=1:nel
n1=NE(1,n); n2=NE(2,n); n3=NE(3,n);%整体编号
s=abs((X(n2)-X(n1))*(Y(n3)-Y(n1))-(X(n3)-X(n1))*(Y(n2)-Y(n1)))/2;%三角形面积 for k=1:3
if n1<=na|n2<=na
T(n1,n2)=T(n1,n2)+((Y(n2)-Y(n3))*(Y(n3)-Y(n1))+(X(n3)-X(n2))*(X(n1)-X(n3)))/(4*s); T(n2,n1)=T(n1,n2);
T(n1,n1)=T(n1,n1)+((Y(n2)-Y(n3))^2+(X(n3)-X(n2))^2)/(4*s);%V=0则边界积分为零，非零时积分编程类似，再加边界积分。 end
k=n1;n1=n2;n2=n3;n3=k; % 轮换坐标将值赋入3阶主子矩阵中 end end
M=T(1:na,1:na);
% 求有限元方程的右端项 f=X;%场源函数 G=zeros(na,1); for n=1:nel
n1=NE(1,n); n2=NE(2,n); n3=NE(3,n);
s=abs((X(n2)-X(n1))*(Y(n3)-Y(n1))-(X(n3)-X(n1))*(Y(n2)-Y(n1)))/2; for k=1:3
if n1<=na
G(n1)=G(n1)+(2*f(n1)+f(n2)+f(n3))*s/12;%f在单元上为线性差值时场域单元的积分公式 end
n4=n1; n1=n2; n2=n3; n3=n4; % 轮换坐标标 end end
F=M\\G; % 求解方程得结果 NNV=zeros(Imax+1,Jmax+1); fi=zeros(ndm,1); fi(1:na)=F(1:na); fi(na+1:ndm)=V; for j=0:Jmax for i=0:Imax
n=NN(i+1,j+1); if n<=0
n=na+1; end
NNV(i+1,j+1)=fi(n);
end end
figure(2)
imagesc(NNV);%画解函数的平面图 X1=zeros(1,Imax+1); Y1=zeros(1,Jmax+1); for i=1:Imax+1
X1(i)=(i-1)*X0; end
for i=1:Jmax+1
Y1(i)=(i-1)*Y0; end
figure(3)
surf(X1,Y1,NNV');% 画解函数的曲面图 % 以下是结果的输出
fid=fopen('Finite_element_tri.txt','w');
fprintf(fid,'\\n *********有限元法求解三角形区域上Possion方程的结果********** \\n \\n'); L=[1:ndm]';
fprintf(fid,'\\n\\n 节点编号 坐标分量x 坐标分量y u（x,y）的值\\n\\n'); for i=1:ndm
fprintf(fid,'?.5f.5f.5f\\n',L(i),X(i),Y(i),fi(i)); end
fclose(fid); End
2、domain_tri.m文件
function domain_tri % 画求解区域 xy=[0 1;0 -1;1 0]; A=zeros(3,3);
A(1,1)=2; A(1,2)=-1;A(1,3)=-1; A(2,2)=2; A(2,1)=-1;A(2,3)=-1; A(3,3)=2; A(3,2)=-1;A(3,1)=-1; A=sparse(A); figure(1); gplot(A,xy);
3、setelm_tri.m文件
function [X,Y,NN,NE]=setelm_tri(Imax,Jmax)
% 给节点和三角形单元编号，并设定节点坐标 % 定义一些全局变量 global ndm nel na
% I1 I2 J1 J2 Imax Jmax分别描述网线纵向和横向数目的变量 % X Y表示节点坐标
% NN描述节点编号
% NE 描述各单元局部节点编号与总体编号对应的矩阵 % ndm 总节点数 % nel 单元数
% na 表示不含边界的节点数 nlm=Imax*Jmax; dx=1/Imax; dy=1/Jmax; X=nlm:1; Y=nlm:1;
NN=zeros(Imax+1,Jmax+1); n1=0;
for j=3:Jmax/2 for i=2:j-1 n1=n1+1;
NN(i,j)=n1; %X=i列，Y=j行处节点 X(n1)=(i-1)*dx; Y(n1)=-1+(j-1)*dy; end end
k=Jmax/2+1;
for j=Jmax/2+1:Jmax-1 %三角形区域上下两部分节点坐标分别求 k=k-1; for i=2:k
n1=n1+1; NN(i,j)=n1; X(n1)=(i-1)*dx;
Y(n1)=1+(j-Jmax-1)*dy; end end
na=n1;%不含边界节点数
for j=Jmax+1:-1:Jmax/2+1 %降序 n1=n1+1; NN(1,j)=n1; X(n1)=0;
Y(n1)=1+(j-Jmax-1)*dy; end
for j=Jmax/2:-1:1 n1=n1+1; NN(1,j)=n1; X(n1)=0;
Y(n1)=-1+(j-1)*dy; end %
for i=2:Imax+1
n1=n1+1; NN(i,i)=n1; X(n1)=(i-1)*dx; Y(n1)=-1+(i-1)*dy; end K=0;
for i=Imax:-1:2 K=K+2; n1=n1+1;
NN(i,i+K)=n1; X(n1)=(i-1)*dx;
Y(n1)=1+(i+K-Jmax-1)*dy; end
% 以上四个循环为对边界节点进行编号 ndm=n1;
NE=zeros(3,2*ndm); n1=0; for j=3:Jmax/2 for i=2:j-1 n1=n1+1;
NE(1,n1)=NN(i,j); NE(2,n1)=NN(i-1,j+1); NE(3,n1)=NN(i-1,j); n1=n1+1;
NE(1,n1)=NN(i,j); NE(2,n1)=NN(i,j+1); NE(3,n1)=NN(i-1,j+1); end end
k=Jmax/2+1;
for j=Jmax/2+1:Jmax-1 k=k-1; for i=2:k
n1=n1+1;
NE(1,n1)=NN(i,j); NE(2,n1)=NN(i-1,j+1); NE(3,n1)=NN(i-1,j); n1=n1+1;
NE(1,n1)=NN(i,j); NE(2,n1)=NN(i,j+1); NE(3,n1)=NN(i-1,j+1); end
end %内部节点对应左上角正方形的两个三角形单元，上左，左上 for i=2:Imax n1=n1+1;