%%
%%  This small was program developed by Luzhong Yin, March 18,2002
%%  University of Pennsylvania
%%
clear all;
%% filename.yin  ------>the file name containing problem specifications
%%% important note: you need to change the input file name for your own problem
%%% below. It should be inserted as below: "!pennsynF filename.yin".
!pennsynF iiscS.yin
% Deformation display scale factor; change to small or large value as 
% needed.
scale = 10;
%% bou.dat     ------> fixed boundary information
%% ouf.dat     ------> output displacement information
%% inf.dat     ------> input load information 
%% ijm.dat     ------> element connetion information
%% xy0.dat     ------> node definition data
%% y00.dis     ------> initial node displacements  
%% y00.num     ------> initial design material densities
%% y01.dis     ------> final node displacements  
%% y01.num     ------> final design material densities
load bou.dat;
load ouf.dat;
load inf.dat;

load ijm.dat;
load xy0.dat;

load y00.dis
dis0=y00;
load y00.num
um0=y00;

load y01.dis
dis=y01;
load y01.num
um=y01;
[ne,nh]=size(ijm);

j=[1:1:4];

for i=1:ne
   x0(j,i)=xy0(ijm(i,j),1);
   x0(5,i)=x0(1,i);
   y0(j,i)=xy0(ijm(i,j),2);
   y0(5,i)=y0(1,i);
end
for i=1:ne
      x(j,i)=xy0(ijm(i,j),1)+scale*dis(ijm(i,j),1);
      x(5,i)=x(1,i);
      y(j,i)=xy0(ijm(i,j),2)+scale*dis(ijm(i,j),2);
      y(5,i)=y(1,i);
end

figure;
hold on;
for i=1:ne
        c=[1-um0(i) 1-um0(i) 1-um0(i)];
%%%%%%%initial density distribution
        fill(x0(:,i),y0(:,i),c) 
end
%%%%%%plot boundary condition
[nb,i]=size(bou);
i=1;
while i < nb
    hold on;
    switch bou(i,1)==bou(i+1,1)
    case 1 
        plot(xy0(bou(i,1),1),xy0(bou(i,1),2),'s');
        i=i+2;
    case 0
        plot(xy0(bou(i,1),1),xy0(bou(i,1),2),'O');
        i=i+1;
    end
end
hold on;
len=abs(max(max(xy0)-min(xy0)))/10.;
%%%%%%plot input information
[next,i]=size(inf);
for i=1:1:next
plot([inf(i,2) inf(i,2)+(2-inf(i,5))*len]',[inf(i,3) inf(i,3)+(inf(i,5)-1)*len]','g-');
hold on;
if inf(i,6)>0
    switch inf(i,5)
    case 1
      plot(inf(i,2)+(2-inf(i,5))*len,inf(i,3)+(inf(i,5)-1)*len,'g>')
    case 2
      plot(inf(i,2)+(2-inf(i,5))*len,inf(i,3)+(inf(i,5)-1)*len,'g^')
    end
end
if inf(i,6)<0
    switch inf(i,5)
    case 1
      plot(inf(i,2),inf(i,3),'g<')
    case 2
      plot(inf(i,2),inf(i,3),'gv')
    end
end
end
hold on;
%%%%%plot output information
[nout,i]=size(ouf);
for i=1:1:nout
plot([ouf(i,1) ouf(i,1)+(2-ouf(i,3))*len]',[ouf(i,2) ouf(i,2)+(ouf(i,3)-1)*len]','r-');
hold on;
if ouf(i,4)>0
    switch ouf(i,3)
    case 1
      plot(ouf(i,1)+(2-ouf(i,3))*len,ouf(i,2)+(ouf(i,3)-1)*len,'r>')
    case 2
      plot(ouf(i,1)+(2-ouf(i,3))*len,ouf(i,2)+(ouf(i,3)-1)*len,'r^')
    end
end
if ouf(i,4)<0
    switch ouf(i,3)
    case 1
      plot(ouf(i,1),ouf(i,2),'r<')
    case 2 
      plot(ouf(i,1),ouf(i,2),'rv')
    end
end
end
    
title('Initial Design')
axis equal
axis off

getaxis = axis;

%%%%%%%% plot the final design
figure;
hold on;
for i=1:ne
   if(um(i)>0.1)  
        c=[1-um(i) 1-um(i) 1-um(i)];
%%%%%%%final density distribution
        fill(x0(:,i),y0(:,i),c) 
%%%%%%%deformed mesh
        %plot(x(:,i),y(:,i),'r')
    end
end
title('Final Design')
axis(getaxis);
axis equal
axis off
%break % Added by Suresh
%%%%%%plot boundary condition
[nb,i]=size(bou);
i=1;
while i < nb
    hold on;
    switch bou(i,1)==bou(i+1,1)
    case 1 
        plot(xy0(bou(i,1),1),xy0(bou(i,1),2),'s');
        i=i+2;
    case 0
        plot(xy0(bou(i,1),1),xy0(bou(i,1),2),'O');
        i=i+1;
    end
end
%%%%%%plot input information
[next,i]=size(inf);
for i=1:1:next
plot([inf(i,2) inf(i,2)+(2-inf(i,5))*len]',[inf(i,3) inf(i,3)+(inf(i,5)-1)*len]','g-');
hold on;
if inf(i,6)>0
    switch inf(i,5)
    case 1
      plot(inf(i,2)+(2-inf(i,5))*len,inf(i,3)+(inf(i,5)-1)*len,'g>')
    case 2
      plot(inf(i,2)+(2-inf(i,5))*len,inf(i,3)+(inf(i,5)-1)*len,'g^')
    end
end
if inf(i,6)<0
    switch inf(i,5)
    case 1
      plot(inf(i,2),inf(i,3),'g<')
    case 2
      plot(inf(i,2),inf(i,3),'gv')
    end
end
end
hold on;
%%%%%plot output information
[nout,i]=size(ouf);
for i=1:1:nout
plot([ouf(i,1) ouf(i,1)+(2-ouf(i,3))*len]',[ouf(i,2) ouf(i,2)+(ouf(i,3)-1)*len]','r-');
hold on;
if ouf(i,4)>0
    switch ouf(i,3)
    case 1
      plot(ouf(i,1)+(2-ouf(i,3))*len,ouf(i,2)+(ouf(i,3)-1)*len,'r>')
    case 2
      plot(ouf(i,1)+(2-ouf(i,3))*len,ouf(i,2)+(ouf(i,3)-1)*len,'r^')
    end
end
if ouf(i,4)<0
    switch ouf(i,3)
    case 1
      plot(ouf(i,1),ouf(i,2),'r<')
    case 2 
      plot(ouf(i,1),ouf(i,2),'rv')
    end
end
end


%clear all;



