%************************************************************************** 
%                                                                         
%                     Truss dCdA         
%                                                                            
%                      Minimize compliance C
%                      subject to volume constraint 
%
%                     Using fmincon()
%
%**************************************************************************
%   Authors:Renato Pavanello
%	Copyright (c) Outubro 2023  by The LTM/FEM/UNICAMP.	
%   $ Revision: 1.0 $					
%**************************************************************************
clear all ;close all;clc

global inci coord nel nnos Tmat fg FreeDofs ug ite vComp vVol Wlim

%% *****************   DATA BLOCK    ***************
Lx = 200;
Ly = 100;
nodesx = 5;
nodesy = 3;
carga = -4400;  %  [N]
Ain=0.065;      % [m^2]
Wlim=100;       % Wheight constraint  

%    struc_type -- type of Boundary conditions for Truss Structure - 
%     Ground  structure  model 
%           1 Cantilever Beam 
%           2 Simply supported beam
struc_type = 1;
%      mesh_type 
%         1 full mesh - all nodes connected 
%         2 local mesh X 
mesh_type = 2;
 
[inci,coord,BC,Loads,nel,nnos] = meshGen(Lx,Ly,nodesx,nodesy,carga,struc_type,mesh_type);

%  TMat = [ E Sigma_yc Sigma_yt  rho ECC]  Steel /  Timber
Tmat=[200E9   345E6   345E6  1  1.45;
       11E9   -8.6E6    6.6E6   570  0.42];

 %   TGeo = [ A  ] 
Tgeo=ones(nel,1); Tgeo=Ain*Tgeo;
 
% Modify areas and material 
%Tgeo(18:29,1)=2*Tgeo(18:29,1);
%inci(1:38,3)=2;    % material timber
 
%% **************** Boundary conditions   ******************
% ****************  Load Vector *************************** 
fg=zeros(2*nnos,1);
nLoads=size(Loads,1);
for i=1:nLoads;
fg( 2*Loads(i,1) - (2-Loads(i,2))  , 1) = fg(2*Loads(i,1)-(2-Loads(i,2)) ,1) + Loads(i,3) ;
end
% **************** Boundary conditions ********************
nBC=size(BC,1);
for i=1:nBC;
FixedDof(i)= 2*BC(i,1) - (2-BC(i,2));
end
allDof=[1:1:2*nnos];
FreeDofs=setdiff(allDof,FixedDof);


% %Plot da malha
% xmin = min(coord(:,1));
% xmax = max(coord(:,1));
% ymin = min(coord(:,2));
% ymax = max(coord(:,2));
% Lx = xmax - xmin;
% Ly = ymax - ymin;
% xmin = xmin - 0.2*Lx;
% xmax = xmax + 0.2*Lx;
% ymin = ymin - 0.2*Ly;
% ymax = ymax + 0.2*Ly;
%  scrsz = get(groot,'ScreenSize');
%  figure('Position',[ scrsz(4)*.5 scrsz(4)*.25 scrsz(3)*.35 scrsz(4)*.65])
% tit=[' Finite Element Mesh '];
% plotMalha2Mat(coord,inci,1,tit,Tgeo,Tmat);   hold on;
% axis([xmin xmax ymin ymax]);
%  

%%****************  Variable initialization ******************
ug=zeros(2*nnos,1);
ite=0;
vComp=[];
vVol=[];
A = [];
b = [];
Aeq = [];
beq = [];




%% **************** Constraints functions   ******************
[c,ceq]=nlcon(Tgeo);
con=@nlcon;


%% ********************     Box Constraints  ****************** 
%  *************   lb ub =Lower and uper bound   **************

lb=0.1*Tgeo;
ub=3*Tgeo;
x0=Tgeo; 


%% ********************   Objective functions   ****************** 
%  ********************      Compliance         ******************

[Compliance]=objective(Tgeo);

Tgeo = fmincon(@objective,x0,A,b,Aeq,beq,lb,ub,@nlcon)



%% ********************   Sensitivity analysis  ****************** 
%  *******        Sensitivuity Analysis dCda = - ue' (ke') ue  
   
 for e=1:nel
     noi=inci(e,5);    
     noj=inci(e,6);
     xi=coord(noi,1);
     yi=coord(noi,2);
     xj=coord(noj,1);
     yj=coord(noj,2);
     l=sqrt((xj-xi)^2+ (yj-yi)^2);
     c=(xj-xi)/l;
     s=(yj-yi)/l;
     E=Tmat(inci(e,3) ,1);
     A=Tgeo(inci(e,4) ,1);
          
     ke= (E/l)*[ c^2      c*s    -c^2   -c*s;
                   c*s      s^2    -c*s   -s^2;
                  -c^2    -c*s     c^2    c*s;
                  -c*s    -s^2     c*s     s^2];
     
         ue=[ ug(2*noi-1) ug(2*noi) ug(2*noj-1) ug(2*noj)]';    
      dCdA(e)=ue'*ke*ue;
      vVol(e)=A*l;
 end
 
 


xmin = min(coord(:,1));
xmax = max(coord(:,1));
ymin = min(coord(:,2));
ymax = max(coord(:,2));
Lx = xmax - xmin;
Ly = ymax - ymin;
xmin = xmin - 0.2*Lx;
xmax = xmax + 0.2*Lx;
ymin = ymin - 0.2*Ly;
ymax = ymax + 0.2*Ly;
escala = (0.1*max(Lx,Ly)/(max(abs(ug))));
a=zeros(nnos,2);
a(:,1)=coord(:,1)+escala*ug(1:2:2*nnos-1);
a(:,2)=coord(:,2)+escala*ug(2:2:2*nnos);
scrsz = get(groot,'ScreenSize');
figure('Position',[ scrsz(4)*.5 scrsz(4)*.25 scrsz(3)*.35 scrsz(4)*.65])
tit='';
subplot(3,1,1);plotMalha(coord,inci,1,tit);   hold on;
tit=[' Finite Element Mesh &  Max. Displamcement [m] = ',num2str((max(abs(ug))))];
subplot(3,1,1);plotMalha(a,inci,0,tit);  
axis([xmin xmax ymin ymax]);
hold off; 

tit=[' Finite Element Mesh '];
subplot(3,1,2);plotMalha2Mat(coord,inci,0,tit,Tgeo,Tmat);   hold on;
axis([xmin xmax ymin ymax]);
 
tit=[' Sensitivity dC/da '];
subplot(3,1,3);plotSxx(coord,inci,dCdA,tit);   
axis([xmin xmax ymin ymax]);

figure('Position',[ scrsz(4)*.5 scrsz(4)*.25 scrsz(3)*.35 scrsz(4)*.65])
subplot(2,1,1);plot(vComp); grid on;
title('Compliance along Iterations')
subplot(2,1,2);plot(vVol);grid on;  
title(' Total Weight along Iterations')
axis([xmin xmax ymin ymax]);