tic;
Initial_info=[20?10?20?10];
disp（[‘該程序計算的是‘num2str（Initial_info（3）+1）‘X‘num2str（Initial_info（4）+1）‘=‘...
????num2str（（Initial_info（3）+1）*（Initial_info（4）+1））‘個結點的有限元模型‘]）;
LX=Initial_info（1）;?%%%矩形的x方向邊長
LY=Initial_info（2）;?%%%矩形的y方向邊長
nx=Initial_info（3）;%%%%%x邊劃分的段數
ny=Initial_info（4）;%%%%%y邊劃分的段數
ne=2*nx*ny;
np=（nx+1）*（ny+1）;
for?i=1:nx+1;
????j=1:ny+1;
????Np（ij）=j+（i-1）*（ny+1）;
end?
%%------生成節點編號矩陣Np
for?i=1:nx+1;
????j=1:ny+1;
????XX（ij）=（i-1）*LX/nx;
????YY（ij）=（j-1）*LY/ny;
end
XY=[reshape（XX‘np1）reshape（YY‘np1）];%各節點坐標
nx2=nx/2;
Np1=Np（1:nx2+1:）;
Np2=Np（nx2+1:end:）;
for?i=1:nx2*ny;
????if?rem（inx2）==0
???????xp=nx2;
???????yp=i/nx2;
????else
????????xp=rem（inx2）;
????????yp=fix（i/nx2）+1;
????end
????Dof1（i:）=[Np1（xpyp）Np1（xp+1yp）Np1（xpyp+1）];
????Dof1（i+nx2*ny:）=[Np1（xp+1yp）Np1（xp+1yp+1）Np1（xpyp+1）];%單元的自由度
????Dof2（i:）=[Np2（xpyp）Np2（xp+1yp）Np2（xp+1yp+1）];
????Dof2（i+nx2*ny:）=[Np2（xpyp）Np2（xp+1yp+1）Np2（xpyp+1）];%單元的自由度
end
Dof=[Dof1;Dof2];
for?i=1:ne
????unit（i:）=[XY（Dof（i1）1）XY（Dof（i2）1）XY（Dof（i3）1）...
????????XY（Dof（i1）2）XY（Dof（i2）2）XY（Dof（i3）2）];%三角形單元的三角點坐標
end
disp（‘前處理完成‘）;
%--------?前處理完成

E=1*10^7;u=0.2;%平面應力問題
D=E/（1-u^2）*[1?u?0;u?1?0;0?0?（1-u）/2];
for?i=1:ne
????xi=unit（i1）;
????yi=unit（i4）;
????xj=unit（i2）;
????yj=unit（i5）;
????xm=unit（i3）;
????ym=unit（i6）;
????ai=xj*ym-xm*yj;
????aj=xm*yi-xi*ym;
????am=xi*yj-xj*yi;
????bi=yj-ym;
????bj=ym-yi;
????bm=yi-yj;
????ci=-（xj-xm）;
????cj=-（xm-xi）;
????cm=-（xi-xj）;
????area=abs（（ai+aj+am）/2）;
????B?=?[bi??0?bj??0?bm??0
??????????0?ci??0?cj??0?cm
?????????ci?bi?cj?bj?cm?bm];
????Be{i1}?=?B/2/area?;
????ke{i1}=[Be{i1}]‘*D*Be{i1}*area;
end

%--------單元剛度矩陣計算完畢

KK=sparse（2*np2*np）;
for?ie=1:ne
????a=Dof（ie1）;
????b=Dof（ie2）;
????c=Dof（ie3）;????
????DOF（1）=2*a-1;
????DOF（2）=2*a;
????DOF（3）=2*b-1;
????DOF（4）=2*b;
????DOF（5）=2*c-1;
????DOF（6）=2*c;
????for?n1=1:6
????????for?n2=1:6
????????????KK（DOF（n1）DOF（n2））=?KK（DOF（n1）DOF（n2））+ke{ie1}（n1n2）;
????????end
????end
end

%------總剛度矩陣疊加

P=500000;%只受上部均布荷載
Re=sparse（ne6）;
for?i=1:ne;
????switch?i
????????case?num2cell（1:ne/2-nx2）
????????????Pe=[0?0?0?0?0?0];
????????case?num2cell（ne/2-nx2+1:ne/2）
????????????Pe=-LX*P*[000101]/nx/2;
????????case?num2cell（ne/2+1:ne-nx2）
????????????Pe=[0?0?0?0?0?0];
????????otherwise
????????????Pe=-LX*P*[000101]/nx/2;
????end
????Re（i:）=Pe;?
end??

%----單元等效節點荷載

Rr=sparse（12*np）;
for?i=1:ne
????a=Dof（i1）;
????b=Dof（i2）;
????c=Dof（i3）;????
????DOF（1）=2*a-1;
????DOF（2）=2*a;
????DOF（3）=2*b-1;
????DOF（4）=2*b;
????DOF（5）=2*c-1;
????DOF（6）=2*c;
????for?n1=1:6
????????Rr（DOF（n1））=?Rr（DOF（n1））+Re（in1）;
????end
end

%---------荷載疊加

cp=[1:nx+1];
ctype=ones（1length（cp））;
ctype（nx2+1）=2;%底部中間固定住
cp_all=（cp-1）*（ny+1）+1;%%%%有約束的節點對應的整體節點編號
p_stake=zeros（12*length（cp））;
for?i=1:length（cp）
?

?屬性????????????大小?????日期????時間???名稱
-----------?---------??----------?-----??----

?????文件???????7708??2010-10-28?21:31??new.m

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

?????????????????7708????????????????????1