finite_element_tri.m

有限元法求解泊松方程

% 用有限元法求解三角形形区域上的Possion方程
function Finite_element_tri(Imax,Jmax)
global ndm nel na
% ndm 总节点数
% nel 基元数
% na 活动节点数
Imax=30;Jmax=60;%设定网格数
V=0; J=0;X0=1/Imax;Y0=X0;
domain_tri   
[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);
        end
        k=n1;n1=n2;n2=n3;n3=k; 
    end
end
M=T(1:na,1:na);
% 求有限元方程的右端项
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*X(n1)+X(n2)+X(n3))*s/12;
        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

% 画等电势线
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(2)
surf(X1,Y1,NNV');
fid=fopen('Finite_element_tri.txt','w');
fprintf(fid,'\n *********有限元法求解三角形区域上Possion方程的结果********** \n \n');
fprintf(fid,'\n 节点编号 \n \n');
Nna=fliplr(NN);
fprintf(fid,'%4d%4d%4d%4d%4d%4d%4d%4d%4d\n',Nna);fprintf(fid,'\n 各节点的电势 \n \n');
NNV=fliplr(NNV);
fprintf(fid,'%10.6f%10.6f%10.6f%10.6f%10.6f%10.6f%10.6f%10.6f%10.6f\n',NNV);
L=[1:ndm]';
fprintf(fid,'\n\n    节点编号    坐标分量x   坐标分量y      u（x,y）的值\n\n');
for i=1:ndm
    fprintf(fid,'%8d%14.5f%14.5f%14.5f\n',L(i),X(i),Y(i),fi(i));
end
fclose(fid);
end 
 
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(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;
    NE(1,n1)=NN(i,i);
    NE(2,n1)=NN(i-1,i);
    NE(3,n1)=NN(i-1,i-1);
    n1=n1+1;
    NE(1,n1)=NN(i,i);
    NE(2,n1)=NN(i-1,i+1);
    NE(3,n1)=NN(i-1,i);
    n1=n1+1;
    NE(1,n1)=NN(i,i);
    NE(2,n1)=NN(i,i+1);
    NE(3,n1)=NN(i-1,i+1);
end

n1=n1+1;
NE(1,n1)=NN(Imax+1,Imax+1);
NE(2,n1)=NN(Imax,Imax+1);
NE(3,n1)=NN(Imax,Imax);
n1=n1+1;
NE(1,n1)=NN(Imax+1,Imax+1);
NE(2,n1)=NN(Imax,Imax+2);
NE(3,n1)=NN(Imax,Imax+1);
K=0;
for i=Imax:-1:2
    K=K+2;
    n1=n1+1;
    NE(1,n1)=NN(i,i+K);
    NE(2,n1)=NN(i-1,i+K+1);
    NE(3,n1)=NN(i-1,i+K);
end
nel=n1;
end
   
function domain_tri
% 定义求解区域
xy=[0 -1;-1 1;1 1;];
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);