在开发有限元程序的时候,首先遇到的问题就是实现自动划分网格(meshgrid),matlab可以实现这一功能,它的基本原理就是根据图论和邻接矩阵的原理利用gplot函数和sparse稀疏矩阵实现的,但是它只能实现二维!本人重新编写了一个程序,可以实现三维的网格自动划分(包括三角形单元、矩形单元和六面体单元),程序如下:
a1=zeros(m1*n1);
b1=zeros(m1*n1);
XY1=zeros(m1*n1,2);
x1=zeros(1,m1*n1);
y1=zeros(1,m1*n1);
z1=zeros(1,m1*n1);
for k=1:n1
for i=1:m1
x1(i+(k-1)*m1)=xmin1+(k-1)*(xmax1-xmin1)/(n1-1);
z1(i+(k-1)*m1)=zmin1+(i-1)*(zmax1-zmin1)/(m1-1);
end
end
clear k i
a1(m1,m1)=2;
a1((n1-1)*m1+1,(n1-1)*m1+1)=2;
a1(1,1)=3;
a1(n1*m1,n1*m1)=3;
for i=1:(n1-2)
a1(i*m1+1,i*m1+1)=4;
a1((i+1)*m1,(i+1)*m1)=4;
end
clear i
for j=2:(m1-1)
a1(j,j)=4;
a1((n1-1)*m1+j,(n1-1)*m1+j)=4;
end
clear j
for i=1:(n1-2)
for j=2:(m1-1)
a1(i*m1+j,i*m1+j)=6;
end
end
clear i j
a1=sparse(a1);
for i=1:(n1-1)*m1-1
b1(i,i+1)=-1;
end
clear i
for i=1:(n1-1)*m1
b1(i,i+m1)=-1;
end
clear i
for i=1:(n1-1)*m1-1
b1(i,i+m1+1)=-1;
end
clear i
for k=1:(n1-1)
for h=1:(n1-1)
b1(m1*k,m1*h+1)=0;
end
end
clear k h
for i=(n1-1)*m1+1:n1*m1-1
b1(i,i+1)=-1;
end
clear i
b1=sparse(b1);
b1=a1+b1'+b1;
clear a1
a1=sparse(b1);
clear b1
XY1=[x1',y1',z1'];
a2=zeros(m2*n2);
b2=zeros(m2*n2);
XY2=zeros(m2*n2,2);
x2=0.5*ones(1,m2*n2);
y2=zeros(1,m2*n2);
z2=zeros(1,m2*n2);
for k=1:n2
for i=1:m2
z2(i+(k-1)*m2)=zmin2+(k-1)*(zmax2-zmin2)/(n2-1);
y2(i+(k-1)*m2)=ymin2+(i-1)*(ymax2-ymin2)/(m2-1);
end
end
clear k i
a2(m2,m2)=2;
a2((n2-1)*m2+1,(n2-1)*m2+1)=2;
a2(1,1)=3;