%************************************************************************** 
%                                                                         
%                     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 ECC ECClim nMax

%% *****************   DATA BLOCK    ***************
Lx = 200;
Ly = 100;
nodesx = 5;
nodesy = 3;
carga = -4400;  %  [N]
Ain=(.0254)^2;      % [m^2]
Wlim=1000;       % Wheight constraint  
ECClim=.6E2;

%    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    7800  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);


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



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

lb=[0.01*Tgeo ; zeros(nel,1)];
ub=[3*Tgeo;   ones(nel,1)];

%%  ******************** Design Variables and initial guess  ***
X=[Tgeo; 0.5*ones(nel,1)];
X0=[Tgeo; 0.5*ones(nel,1)]; 

%% ********************   Objective functions   ****************** 
%  ********************      Compliance         ******************
ite=0;

[Compliance]=objective(X);



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

%% *****************  Solver *********************************
X = 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=X(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;
   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,X(1:nel),Tmat,X(nel+1:2*nel));   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(3,1,1);plot(vComp); grid on;
title('Compliance along Iterations')
subplot(3,1,2);plot(vVol);grid on;  
title(' Total Weight along Iterations')

subplot(3,1,3);plot(ECC);grid on;  
title(' Total ECC along Iterations')