matlab网格自动划分,有限元实现三维自动划分网格.m

% http://myangie.blog.sohu.com/6192144.html

% 在开发有限元程序的时候，

% 首先遇到的问题就是实现自动划分网格(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

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

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

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

clear k i

a2(m2,m2)=2;

a2((n2-1)*m2+1,(n2-1)*m2+1)=2;

a2(1,1)=3;

a2(n2*m2,n2*m2)=3;

for i=1:(n2-2)

a2(i*m2+1,i*m2+1)=4;

a2((i+1)*m2,(i+1)*m2)=4;

end

clear i

for j=2:(m2-1)

a2(j,j)=4;

a2((n2-1)*m2+j,(n2-1)*m2+j)=4;

end

clear j

for i=1:(n2-2)

for j=2:(m2-1)

a2(i*m2+j,i*m2+j)=6;

end

clear i j

a2=sparse(a2);

for i=1:(n2-1)*m2-1

b2(i,i+1)=-1;

end

clear i

for i=1:(n2-1)*m2

b2(i,i+m2)=-1;

end

clear i

for i=1:(n2-1)*m2-1

b2(i,i+m2+1)=-1;

end

clear i

for k=1:(n2-1)

for h=1:(n2-1)

b2(m2*k,m2*h+1)=0;

end

clear k h

for i=(n2-1)*m2+1:n2*m2-1

b2(i,i+1)=-1;

end

clear i

b2=sparse(b2);

b2=a2+b2'+b2;

clear a2

a2=sparse(b2);

clear b2

XY2=[x2',y2',z2'];

a3=zeros(m3*n3);

b3=zeros(m3*n3);

XY3=zeros(m3*n3,2);

x3=1.5*ones(1,m3*n3);

y3=zeros(1,m3*n3);

z3=zeros(1,m3*n3);

for k=1:n3

for i=1:m3

z3(i+(k-1)*m3)=zmin3+(k-1)*(zmax3-zmin3)/(n3-1);

y3(i+(k-1)*m3)=ymin3+(i-1)*(ymax3-ymin3)/(m3-1);

end

clear k i

a3(m3,m3)=2;

a3((n3-1)*m3+1,(n3-1)*m3+1)=2;

a3(1,1)=3;

a3(n3*m3,n3*m3)=3;

for i=1:(n3-2)

a3(i*m3+1,i*m3+1)=4;

a3((i+1)*m3,(i+1)*m3)=4;

end

clear i

for j=2:(m3-1)

a3(j,j)=4;

a3((n3-1)*m3+j,(n3-1)*m3+j)=4;

end

clear j

for i=1:(n3-2)

for j=2:(m3-1)

a3(i*m3+j,i*m3+j)=6;

end

clear i j

a3=sparse(a3);

for i=1:(n3-1)*m3-1

b3(i,i+1)=-1;

end

clear i

for i=1:(n3-1)*m3

b3(i,i+m3)=-1;

end

clear i

for i=1:(n3-1)*m3-1

b3(i,i+m3+1)=-1;

end

clear i

for k=1:(n3-1)

for h=1:(n3-1)

b3(m3*k,m3*h+1)=0;

end

clear k h

for i=(n3-1)*m3+1:n3*m3-1

b3(i,i+1)=-1;

end

clear i

b3=sparse(b3);

b3=a3+b3'+b3;

clear a3

a3=sparse(b3);

clear b3

XY3=[x3',y3',z3'];

a4=zeros(m4*n4);

b4=zeros(m4*n4);

XY4=zeros(m4*n4,2);

x4=0*ones(1,m4*n4);

y4=zeros(1,m4*n4);

z4=zeros(1,m4*n4);

for k=1:n4

for i=1:m4

z4(i+(k-1)*m4)=zmin4+(k-1)*(zmax4-zmin4)/(n4-1);

y4(i+(k-1)*m4)=ymin4+(i-1)*(ymax4-ymin4)/(m4-1);

end

clear k i

a4(m4,m4)=2;

a4((n4-1)*m4+1,(n4-1)*m4+1)=2;

a4(1,1)=3;

a4(n4*m4,n4*m4)=3;

for i=1:(n4-2)

a4(i*m4+1,i*m4+1)=4;

a4((i+1)*m4,(i+1)*m4)=4;

end

clear i

for j=2:(m4-1)

a4(j,j)=4;

a4((n4-1)*m4+j,(n4-1)*m4+j)=4;

end

clear j

for i=1:(n4-2)

for j=2:(m4-1)

a4(i*m4+j,i*m4+j)=6;

end

clear i j

a4=sparse(a4);

for i=1:(n4-1)*m4-1

b4(i,i+1)=-1;

end

clear i

for i=1:(n4-1)*m4

b4(i,i+m4)=-1;

end

clear i

for i=1:(n4-1)*m4-1

b4(i,i+m4+1)=-1;

end

clear i

for k=1:(n4-1)

for h=1:(n4-1)

b4(m4*k,m4*h+1)=0;

end

clear k h

for i=(n4-1)*m4+1:n4*m4-1

b4(i,i+1)=-1;

end

clear i

b4=sparse(b4);

b4=a4+b4'+b4;

clear a4

a4=sparse(b4);

clear b4

XY4=[x4',y4',z4'];

a5=zeros(m5*n5);

b5=zeros(m5*n5);

XY5=zeros(m5*n5,2);

x5=2*ones(1,m5*n5);

y5=zeros(1,m5*n5);

z5=zeros(1,m5*n5);

for k=1:n5

for i=1:m5

z5(i+(k-1)*m5)=zmin5+(k-1)*(zmax5-zmin5)/(n5-1);

y5(i+(k-1)*m5)=ymin5+(i-1)*(ymax5-ymin5)/(m5-1);

end

clear k i

a5(m5,m5)=2;

a5((n5-1)*m5+1,(n5-1)*m5+1)=2;

a5(1,1)=3;

a5(n5*m5,n5*m5)=3;

for i=1:(n5-2)

a5(i*m5+1,i*m5+1)=4;

a5((i+1)*m5,(i+1)*m5)=4;

end

clear i

for j=2:(m5-1)

a5(j,j)=4;

a5((n5-1)*m5+j,(n5-1)*m5+j)=4;

end

clear j

for i=1:(n5-2)

for j=2:(m5-1)

a5(i*m5+j,i*m5+j)=6;

end

clear i j

a5=sparse(a5);

for i=1:(n5-1)*m5-1

b5(i,i+1)=-1;

end

clear i

for i=1:(n5-1)*m5

b5(i,i+m5)=-1;

end

clear i

for i=1:(n5-1)*m5-1

b5(i,i+m5+1)=-1;

end

clear i

for k=1:(n5-1)

for h=1:(n5-1)

b5(m5*k,m5*h+1)=0;

end

clear k h

for i=(n5-1)*m5+1:n5*m5-1

b5(i,i+1)=-1;

end

clear i

b5=sparse(b5);

b5=a5+b5'+b5;

clear a5

a5=sparse(b5);

clear b5

XY5=[x5',y5',z5'];

A1=a1;

xy1=XY1;

A2=a2;

xy2=XY2;

A3=a3;

xy3=XY3;

A4=a4;

xy4=XY4;

A5=a5;

xy5=XY5;

[i1,j1] = find(A1);

[ignore, p1] = sort(max(i1,j1));

i1=i1(p1);

j1=j1(p1);

[i2,j2] = find(A2);

[ignore, p2] = sort(max(i2,j2));

i2=i2(p2);

j2=j2(p2);

[i3,j3] = find(A3);

[ignore, p3] = sort(max(i3,j3));

i3=i3(p3);

j3=j3(p3);

[i4,j4] = find(A4);

[ignore, p4] = sort(max(i4,j4));

i4=i4(p4);

j4=j4(p4);

[i5,j5] = find(A5);

[ignore, p5] = sort(max(i5,j5));

i5=i5(p5);

j5=j5(p5);

% Create a long, NaN-separated list of line segments,

% rather than individual segments.

X1 = [ xy1(i1,1) xy1(j1,1) repmat(NaN,size(i1))]';

Y1 = [ xy1(i1,3) xy1(j1,3) repmat(NaN,size(i1))]';

Z1 = [ xy1(i1,2) xy1(j1,2) repmat(NaN,size(i1))]';

X2 = [ xy2(i2,1) xy2(j2,1) repmat(NaN,size(i2))]';

Y2 = [ xy2(i2,3) xy2(j2,3) repmat(NaN,size(i2))]';

Z2 = [ xy2(i2,2) xy2(j2,2) repmat(NaN,size(i2))]';

X3 = [ xy3(i3,1) xy3(j3,1) repmat(NaN,size(i3))]';

Y3 = [ xy3(i3,3) xy3(j3,3) repmat(NaN,size(i3))]';

Z3 = [ xy3(i3,2) xy3(j3,2) repmat(NaN,size(i3))]';

X4 = [ xy4(i4,1) xy4(j4,1) repmat(NaN,size(i4))]';

Y4 = [ xy4(i4,3) xy4(j4,3) repmat(NaN,size(i4))]';

Z4 = [ xy4(i4,2) xy4(j4,2) repmat(NaN,size(i4))]';

X5 = [ xy5(i5,1) xy5(j5,1) repmat(NaN,size(i5))]';

Y5 = [ xy5(i5,3) xy5(j5,3) repmat(NaN,size(i5))]';

Z5 = [ xy5(i5,2) xy5(j5,2) repmat(NaN,size(i5))]';

X1 = X1(:);

Y1 = Y1(:);

Z1 = Z1(:);

X2 = X2(:);

Y2 = Y2(:);

Z2 = Z2(:);

X3 = X3(:);

Y3 = Y3(:);

Z3 = Z3(:);

X4 = X4(:);

Y4 = Y4(:);

Z4 = Z4(:);

X5 = X5(:);

Y5 = Y5(:);

Z5 = Z5(:);

% Create the line properties string

cla;

if nargout==0,

subplot('position',[0.05,0.05,0.5,0.8]);

h1=plot3(X1,Y1,Z1);

hold on;

h2=plot3(X2,Y2,Z2);

hold on;

h3=plot3(X3,Y3,Z3);

hold on;

h4=plot3(X4,Y4,Z4);

hold on;

h5=plot3(X5,Y5,Z5);

else

X1out = X1;

Y1out = Y1;

Z1out = Z1;

X2out = X2;

Y2out = Y2;

Z2out = Z2;

X3out = X3;

Y3out = Y3;

Z3out = Z3;

X4out = X4;

Y4out = Y4;

Z4out = Z4;

X5out = X5;

Y5out = Y5;

Z5out = Z5;

end

if nargin>2

for k=1:nargin-3

s=varargin{k};

try

if isnumeric(varargin{k+1})

set(h1,s,varargin{k+1});

set(h2,s,varargin{k+1});

set(h3,s,varargin{k+1});

set(h4,s,varargin{k+1});

set(h5,s,varargin{k+1});

else

set(h1,s);

set(h2,s);

set(h3,s);

set(h4,s);

set(h5,s);

end

catch

end

一键复制

编辑

Web IDE

原始数据

按行查看

历史