% A program to generate a grid for a domain consisting of
% rectangular blocks (corrected 4/17/02)

xp(1,1:4)=[0,0,1,2];xp(2,1:4)=[0.5,0.5,1,1.5];
yp(1,1:4)=[1,0,0,0.5];yp(2,1:4)=[1,0.5,0.5,1];

NumBlocks=3;
Nblock(1)=8;Nblock(2)=5;Nblock(3)=12;NTot=0;
Mblock(1)=6;Mblock(2)=6;Mblock(3)=6;MTot=6;

for l=1:NumBlocks
  N=Nblock(l);M=Mblock(l);
  x0=xp(1,l);x1=xp(1,l+1);x3=xp(2,l);x2=xp(2,l+1);
  y0=yp(1,l);y1=yp(1,l+1);y3=yp(2,l);y2=yp(2,l+1);
  Nf=2;if(l == 1),Nf=1;end
  for i=Nf:N, for j=1:M	
      ii=i+NTot;
      x(ii,j)=((M-j)/(M-1))*(((N-i)/(N-1))*x0+((i-1)/(N-1))*x1)+...
             ((j-1)/(M-1))*(((N-i)/(N-1))*x3+((i-1)/(N-1))*x2);
      y(ii,j)=((M-j)/(M-1))*(((N-i)/(N-1))*y0+((i-1)/(N-1))*y1)+...
             ((j-1)/(M-1))*(((N-i)/(N-1))*y3+((i-1)/(N-1))*y2);
    end, end;
  NTot=NTot+Nblock(l)-1;
end
NTot=NTot+1;  %corrected 4/17/02

axis('equal'), hold on
for j=1:MTot;plot(x(1:NTot,j),y(1:NTot,j),'linewidth',2);end;
for i=1:NTot;plot(x(i,1:MTot),y(i,1:MTot),'linewidth',2);end;
axis('square'),axis([-0.10, 2, 0, 2]); hold on
axis off

pause
for m=1:2
	for i=2:NTot-1
	for j=2:MTot-1
		x(i,j)=0.25*(x(i+1,j)+x(i-1,j)+x(i,j+1)+x(i,j-1));
		y(i,j)=0.25*(y(i+1,j)+y(i-1,j)+y(i,j+1)+y(i,j-1));
	end;
	end;
end;
axis('square'),axis([-0.10, 2, 0, 2]); hold on
for j=1:MTot;plot(x(1:NTot,j),y(1:NTot,j),'r','linewidth',2);end;
for i=1:NTot;plot(x(i,1:MTot),y(i,1:MTot),'r','linewidth',2);end;

axis off