%************************************************************************** 
%                                                                         
%                     Truss          
%       Minize     C= alfa_1 *fg^T*ug +alfa_2* ECC 
%         A_e
%      Subject to   Kg(A_e) * ug(A_e) = fg
%                   g_w= \SUM^nel( We) <= Wf
%                   Amin<A(ite)<Amax
%
%
%**************************************************************************
%   Authors:Renato Pavanello
%	Copyright (c) Outubro 2023  by The LTM/FEM/UNICAMP.	
%   $ Revision: 1.0 $					
%**************************************************************************


%-------------------------------------------------------------------------%
% Evolutionary Compliance Optimization - Based on BESO method             %
%-------------------------------------------------------------------------%
clear all ;close all;clc

%% ********************  DATA BLOCK ***************************************

Lx = 200;           % Total X Lenght in [m]
Ly = 100;           % Total y lenght in [m]
nodesx = 5;         % number of nodes in X direction
nodesy = 3;         % number of nodes in Y direction
carga = -100000;    %  Total Force in [N] (aplied in one r two nodes.) 
Lsize=6;
Wf=1850;
%                  variávei de estado

Ain=1;              % Inicial area in [m^2]
Amin=0.001*Ain;       % pode ser considerado vazio
Amax=3*Ain;         % afrouchar a variávellaterla se for o caso
Ai=0.1*Ain;


%   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 X mesh 
mesh_type = 2;
 
[inci,coord,BC,Loads,nel,nnos] = meshGen(Lx,Ly,nodesx,nodesy,carga,struc_type,mesh_type);

%       Material properties 
%   TMat = [ E Sigma_yc Sigma_yt  rho ECC]   lines 1 Steel / 2 Timber
Tmat=[200E9   345E6   345E6   7800  1.45;
      11E9   -8.6E6    6.6E6   570  0.42];
  
%       Geometric properties
Tgeo=ones(nel,1); Tgeo=Ain*Tgeo;
 
%       malha com duas areas e dois materiais 
%Tgeo(8:11,1)=2*Tgeo(18:29,1);
inci(23:38,3)=2;    % material timber
 
%% ************* Loads and  Boundary conditions  **************************                    
%   Loads 
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);

%%                mesh plot

%tit=[' Finite Element Mesh '];
%plotMalha2Mat(coord,inci,1,tit,Tgeo,Tmat,Lsize);

%%  ***************        Optimization Parameter      ****************

tau=0.0001;              % Stop criteria limit - Variation of obj. function 
%  elem_list = 1:nel;        %Lista de elementos
%  no_list = 1:nnos;         %Lista de nos
N=5;                    % Iteration for error convergence
nel_met=2;              % number of elements that change for iterations

%   ******************   Stiffness Matrix Update **********************  
[kg]=montaKg(coord,inci,Tgeo,Tmat,nnos,nel);

erro=1;
ite=0;
dCdATotal=[];
WeTotal=[];
ECCTotal=[];
Vol_ite=Wf;
while erro>tau & ite<250 & Vol_ite>=Wf
     ite=ite+1;
     
%%   ****************** Solver Linear ******************** 
ug=zeros(2*nnos,1);
ug(FreeDofs)=kg(FreeDofs,FreeDofs)\fg(FreeDofs);

%%    **************    Função Objetivo    ******************

C=fg'*ug;
%%    **************    Sensitivity 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;
      rho=Tmat(inci(e,3),4);
      We(e)=A*l*rho;
      ECC(e)=A*l*rho*Tmat(inci(e,3),5);
 end
 
 dCdATotal(ite)=sum(dCdA);
 WeTotal(ite)=sum(We);
 ECCTotal(ite)=sum(ECC);

% *****************  Compliance Normalization 
   if ite==1  
           nMax=max(dCdA);
       else
       end
dCdA=(1/nMax)*dCdA;              %  Initial Normalization dCdA  

[a,b] = sort(dCdA,'descend');    %  Sort of dCdA  

%   *************   Add and Remove Area section Ai
    i=0;
    j=0;
    while j<nel_met+1
    if Tgeo(b(end-i))>Amin
    Tgeo(b(end-i))=Tgeo(b(end-i))-Ai;j=j+1;
    else 
      Tgeo(b(end-i))=Amin;
    end
        i=i+1;
  
     end
 
    i=1;
    j=0;
    while j<(nel_met-1)
    if Tgeo(b(i))<Amax
    Tgeo(b(i))=Tgeo(b(i))+Ai;
        j=j+1;
    end
        i=i+1;
    end

    for i=1:nel
        if Tgeo(i)<Amin
            Tgeo(i)=Amin;
        end
    end
    
vArea(:,ite)=Tgeo;
vC(ite)=C;
vdCdA(:,ite)= dCdA;
Vol(ite)=sum(vVol);
Vol_ite=sum(vVol);
%   *****************    Matrix Updated     *************************

[kg]=montaKg(coord,inci,Tgeo,Tmat,nnos,nel);

% *******************   Error calculation ***************************

if ite>=2*N
   erro=abs(sum(vC(ite-N+1:ite))- sum(vC(ite-N*2+1:ite-N)) )/ sum(vC(ite-N+1:ite));
end


end    % loop while 


% %   
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(2,1,1);plotMalha(coord,inci,1,tit);   hold on;
tit=[' Finite Element Mesh &  Max. Displamcement [m] = ',num2str((max(abs(ug))))];
subplot(2,1,1);plotMalha(a,inci,0,tit);  
axis([xmin xmax ymin ymax]);
hold off

 
tit=[' Finite Element Mesh '];

subplot(2,1,2);plotMalha2Mat(coord,inci,0,tit,Tgeo,Tmat,Lsize);   hold on;
axis([xmin xmax ymin ymax]);

figure;
yyaxis left
plot(1:ite,vC); grid on;hold on;
yyaxis right
plot(1:ite,Vol);grid on;
legend('Compliance along Iterations','Volume along Iterations'); 
hold off;