%********************************************************************
%                         ECO_TOP_2
%
%                    Truss static Analysis 
%                                2D 
%                   Gradient Based Design 
%                          Compliance 
%               ECC - Embodied Carbon Coeficient
%             
%                             
%
%   Author: Renato Pavanello
%   Copyright (c) 2023 by  LTM/UNICAMP
%********************************************************************
close all
clear all
all_fig = findall(0, 'type', 'figure');
close(all_fig)
clc


%********************************************************************
%%  DATA BLOCK MESH Generation 
% geometric data 
Lx = 200;
Ly = 100;
nodesx = 3;
nodesy = 2;
Load = -4400;  %  [N]
Area=6.5E-4;  % [m^2]

[inci,coord,BC,Loads,nel,nnos] = meshGen(Lx,Ly,nodesx,nodesy,Load,2,2);


%  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];
%inci(8:11,3)=2;    % material timber

        
 %   Tgeo = [ A  ] 
x=ones(nel,1);
%x(8:11)=4;
Tgeo=Area*x;

%********************************************************************
 
 
 %%  *************** Solver  *******************************************
 %   Passo zero  
 %   Assemblage of global  stiffness matrix 
kg=zeros(2*nnos);
ug=zeros(2*nnos,1);
fg=zeros(2*nnos,1);
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*A/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];
     
     loc=[2*noi-1   2*noi  2*noj-1   2*noj];
     kg(loc,loc)=kg(loc,loc)+ke;
end

%% ********  Vetor de cargas ********** 

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

%% ******** Solver Linear ************ 
ug(FreeDofs)=kg(FreeDofs,FreeDofs)\fg(FreeDofs);
%********************************************************************




%% ******** Pós - processamento *******
% ****** Calculo das tensoes ******
  for e=1:nel
      no1=inci(e,5);
      no2=inci(e,6);
      x1=coord(no1,1);
      y1=coord(no1,2);
      x2=coord(no2,1);
      y2=coord(no2,2);
      le=sqrt((x2-x1)^2+(y2-y1)^2);
      c=(x2-x1)/le;
      s=(y2-y1)/le;
      E=Tmat(inci(e,3),1);
      Sxx(e)=(E/le)*( (c*(ug(2*no2-1)-ug(2*no1-1))) + (s*(ug(2*no2)-ug(2*no1))) );
  end
  
  %  Sensitivuity Analysis dCda = - ue' (ke') ue  and  ECC
   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);
     h=sqrt((xj-xi)^2+ (yj-yi)^2);
     c=(xj-xi)/h;
     s=(yj-yi)/h;
     E=Tmat(inci(e,3) ,1);
     A=Tgeo(inci(e,4) ,1);
     rho= Tmat(inci(e,3) ,4);
     ke= (E/h)*[ 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;
      We(e)=rho*A*h*9.81;
      ECCe(e)=(rho*A*h*9.81)*Tmat(inci(e,3) ,5);
   end
   % calculo da compliance
   Wt=sum(We);
   ECCt=sum(ECCe);
  C=fg'*ug; 
   
%% ******** Plotting the FEM mesh, Results and display numbers  ***************
        
%figure1 = figure('Color',[1 1 1]);
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;
tit='';
plotMalha(coord,inci,1,tit);   hold on;
tit=[' Finite Element Mesh &  Max. Displamcement [m] = ',num2str((max(abs(ug))))];
plotMalha(a,inci,0,tit);  
axis([xmin xmax ymin ymax]);
hold off;

figure;
tit='Bar stress \sigma_{xx} [N/m^2]';
plotSxx(coord,inci,Sxx,tit);   
axis([xmin xmax ymin ymax]);

figure;
tit=[' Area of Elements'];
plotMalha2Mat(coord,inci,1,tit,Tgeo,Tmat);   hold on;
axis([xmin xmax ymin ymax]);
 
figure;
tit=[' Sensitivity dC/da and Compliance C= ', num2str(C), '  [N m]'];
plotSxx(coord,inci,dCdA,tit);   
axis([xmin xmax ymin ymax]);


%%  ******************    Data Viewer    **************************************
         
fig = uifigure;
tdata = table([1:nel]',inci(:,3),[Tgeo(:,1)],Sxx', dCdA',We',ECCe','VariableNames', {'Element','Material n.', 'Area','Stress','Sensitivity dCdA','weight [kg]',' ECCe [kg/kg]'});
uit = uitable(fig,'Data',tdata);
uit.Position(3) = 400;
%uit.RowName = 'numbered';

disp([' Compliance C [Nm] =', num2str(fg'*ug)]);
disp([' ECC total  [kg/kg] =', num2str(ECCt)]);
disp([' Total weight [Kg] =', num2str(Wt)]);