%************************************************************************** 
%                                                                         
%                     Truss    
%             Finite Element Analysis                                                             
%           
%           Variable: 
%                     displacements  (nodes)
%                     Compliance (element - global)
%                     Weight  (element - global)
%                     ECC      (element - global)
%                     Stress    (element)
%
%**************************************************************************
%   Authors:Renato Pavanello
%	Copyright (c) Outubro 2023  by The LTM/FEM/UNICAMP.	
%   $ Revision: 1.0 $					
%**************************************************************************

clear all ;close all;clc

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

Lx = 200;           % Total X Lenght in [m]
Ly = 100;           % Total y lenght in [m]
nodesx = 3;         % number of nodes in X direction
nodesy = 2;         % number of nodes in Y direction
carga = -100000;    %  Total Force in [N] (aplied in one r two nodes.) 
Ain=1;              % Inicial area in [m^2]

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

 
 %      mesh plot 
 %tit=[' Finite Element Mesh '];
 %plotMalha2Mat(coord,inci,1,tit,Tgeo,Tmat);
 
%% ************* 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);


%% ************************************************************************
% Assemblage of global Matriz and solve 
%  
[kg]=montaKg(coord,inci,Tgeo,Tmat,nnos,nel);

%% ******************* Solve Linear System  ******************************* 
ug=zeros(2*nnos,1);
ug(FreeDofs)=kg(FreeDofs,FreeDofs)\fg(FreeDofs);


%%            Pos processing 
%    * Compliance and dCdA
%    * Weight  and dWdA
%    * ECC     dECCdA
%

 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;
      ECCe(e)=A*l*rho*Tmat(inci(e,3),5);
      Sxx(e)=(E/l)*( (c*(ug(2*noj-1)-ug(2*noi-1))) + (s*(ug(2*noj)-ug(2*noi))) );
 end
 
 C=sum(dCdA);        % Compliance
 W=sum(We);          % Weight   
 ECC=sum(ECCe);       % Embodied Carbon - Carbon integré
 
 disp('     Compliance      Weight      ECC');
 format short e
 disp( [C  W  ECC]);
 
%%                    Ploting results    
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]);

Lsize=2;     % Control of line thickness

tit='';
subplot(4,1,2);plotMalha2Mat(coord,inci,1,tit,Tgeo,Tmat,Lsize);   hold on;
tit=[' Finite Element Mesh &  Max. Displamcement [m] = ',num2str((max(abs(ug))))];
subplot(4,1,2);plotMalha2Mat(a,inci,0,tit,Tgeo,Tmat,Lsize);  
axis([xmin xmax ymin ymax]);
hold off; 

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

tit=[' Stress [Pa] '];
subplot(4,1,4);plotSxx(coord,inci,Sxx,tit);   
axis([xmin xmax ymin ymax]);
hold on;