%?用有限元法求解三角形形區(qū)域上的Possion方程
function?Finite_element_tri（ImaxJmax）
global?ndm?nel?na
%?ndm?總節(jié)點(diǎn)數(shù)
%?nel?基元數(shù)
%?na?活動(dòng)節(jié)點(diǎn)數(shù)
Imax=30;Jmax=60;%設(shè)定網(wǎng)格數(shù)
V=0;?J=0;X0=1/Imax;Y0=X0;
domain_tri???
[XYNNNE]=setelm_tri（ImaxJmax）;??%?給節(jié)點(diǎn)和三角形元素編號(hào)，并設(shè)定節(jié)點(diǎn)坐標(biāo)
%?求解有限元方程的求系數(shù)矩陣
T=zeros（ndmndm）;
for?n=1:nel
????n1=NE（1n）;???n2=NE（2n）;???n3=NE（3n）;
????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（n1n2）=T（n1n2）+（（Y（n2）-Y（n3））*（Y（n3）-Y（n1））+（X（n3）-X（n2））*（X（n1）-X（n3）））/（4*s）;
????????T（n2n1）=T（n1n2）;
????????T（n1n1）=T（n1n1）+（（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:na1:na）;
%?求有限元方程的右端項(xiàng)
G=zeros（na1）;
for?n=1:nel
????n1=NE（1n）;???n2=NE（2n）;???n3=NE（3n）;??
????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

%??求解方程得結(jié)果
F=M\G;?
NNV=zeros（Imax+1Jmax+1）;
fi=zeros（ndm1）;
fi（1:na）=F（1:na）;
fi（na+1:ndm）=V;
for?j=0:Jmax
????for?i=0:Imax
????????n=NN（i+1j+1）;
????????if?n<=0
????????????n=na+1;
????????end
????????NNV（i+1j+1）=fi（n）;
????end
end

%?畫等電勢(shì)線
X1=zeros（1Imax+1）;
Y1=zeros（1Jmax+1）;
for?i=1:Imax+1
????X1（i）=（i-1）*X0;
end
for?i=1:Jmax+1
????Y1（i）=（i-1）*Y0;
end

%?畫解函數(shù)的曲面圖
figure（2）
surf（X1Y1NNV‘）;
fid=fopen（‘Finite_element_tri.txt‘‘w‘）;
fprintf（fid‘\n?*********有限元法求解三角形區(qū)域上Possion方程的結(jié)果**********?\n?\n‘）;
fprintf（fid‘\n?節(jié)點(diǎn)編號(hào)?\n?\n‘）;
Nna=fliplr（NN）;
fprintf（fid‘%4d%4d%4d%4d%4d%4d%4d%4d%4d\n‘Nna）;fprintf（fid‘\n?各節(jié)點(diǎn)的電勢(shì)?\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????節(jié)點(diǎn)編號(hào)????坐標(biāo)分量x???坐標(biāo)分量y??????u（xy）的值\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?[XYNNNE]=setelm_tri（ImaxJmax）
%?給節(jié)點(diǎn)和三角形元素編號(hào)，并設(shè)定節(jié)點(diǎn)坐標(biāo)
global??ndm?nel?na
%?I1?I2?J1?J2?Imax?Jmax分別描述網(wǎng)線縱向和橫向數(shù)目的變量
%?X?Y表示節(jié)點(diǎn)坐標(biāo)
%?NN描述節(jié)點(diǎn)編號(hào)
%?NE?描述各基點(diǎn)局域節(jié)點(diǎn)的矩陣
%?ndm?總節(jié)點(diǎn)數(shù)
%?nel?基元數(shù)
%?na?表示活動(dòng)節(jié)點(diǎn)數(shù)
nlm=Imax*Jmax;
dx=1/Imax;
dy=1/Jmax;
X=nlm:1;
Y=nlm:1;
NN=zeros（Imax+1Jmax+1）;
n1=0;?
%?活動(dòng)節(jié)點(diǎn)編號(hào)
for?j=3:Jmax/2
????for?i=2:j-1
????????n1=n1+1;
????????NN（ij）=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（ij）=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（1j）=n1;
????X（n1）=0;
????Y（n1）=1+（j-Jmax-1）*dy;
end
for?j=Jmax/2:-1:1
????n1=n1+1;
????NN（1j）=n1;
????X（n1）=0;
????Y（n1）=-1+（j-1）*dy;
end
for?i=2:Imax+1
????n1=n1+1;
????NN（ii）=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（ii+K）=n1;
????X（n1）=（i-

?屬性????????????大小?????日期????時(shí)間???名稱
-----------?---------??----------?-----??----

?????文件???????4717??2007-01-13?21:11??Finite_element_tri.m

-----------?---------??----------?-----??----

?????????????????4935????????????????????2