有限元矩形网格剖分

本程序采用MATLAB编写，可以自动进行矩形网格的剖分，以及图形的生成。

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