% FEAcomp. Solution of Example 5.1. 
% Ever J. Barbero (c) 2007, 2013
clear all; close all;

bx = 1;
by = 1;

%%% Define loads [N] [Nmm]
Nx = 175e3/bx; 
Ny = 0/by; 
Nxy= 0/bx/by; 
Mx = 0/bx; 
My = 0/by; 
Mxy= 0/bx/by;

n_file='0_90_s'
n_laminate='[0/0/90/90]s'

%%% Laminate Input Data for transversaly isotropic materials, as
%%% E11 (MPa), E22 (MPa), G12 (MPa), nu12, angle (º), thickness (mm)

L(1,:)=[139e3, 14.5e3, 5.86e3, 0.21, 0, 1];    
L(2,:)=[139e3, 14.5e3, 5.86e3, 0.21, 0, 1];    
L(3,:)=[139e3, 14.5e3, 5.86e3, 0.21, 90, 1];    
L(4,:)=[139e3, 14.5e3, 5.86e3, 0.21, 90, 1];    
L(5,:)=[139e3, 14.5e3, 5.86e3, 0.21, 90, 1];    
L(6,:)=[139e3, 14.5e3, 5.86e3, 0.21, 90, 1];    
L(7,:)=[139e3, 14.5e3, 5.86e3, 0.21, 0, 1];    
L(8,:)=[139e3, 14.5e3, 5.86e3, 0.21, 0, 1];


% n_file='45_0_90_s'
% n_laminate='[45/-45/0/90]s'
% 
% %%% Laminate Input Data for transversaly isotropic materials, as
% %%% E11 (MPa), E22 (MPa), G12 (MPa), nu12, angle (º), thickness (mm)
% 
% L(1,:)=[139e3, 14.5e3, 5.86e3, 0.21, 45, 1];    
% L(2,:)=[139e3, 14.5e3, 5.86e3, 0.21,-45, 1];    
% L(3,:)=[139e3, 14.5e3, 5.86e3, 0.21,  0, 1];    
% L(4,:)=[139e3, 14.5e3, 5.86e3, 0.21, 90, 1];    
% L(5,:)=[139e3, 14.5e3, 5.86e3, 0.21, 90, 1];    
% L(6,:)=[139e3, 14.5e3, 5.86e3, 0.21,  0, 1];    
% L(7,:)=[139e3, 14.5e3, 5.86e3, 0.21,-45, 1];    
% L(8,:)=[139e3, 14.5e3, 5.86e3, 0.21, 45, 1];    

%%% The first layer is the bottom layer....
%%% Define equal number of L rows as laminae number

[M,N]=size(L);

%%% Compute stiffness and compliance matrices of each lamina

%% Using material lamina axis (1 - long. fiber, 2 - transversal fiber)
Q=zeros(3,3,M);

for i=1:M;
    Q(1,1,i)=L(i,1)/(1-L(i,4)^2*L(i,2)/L(i,1));
    Q(1,2,i)=L(i,4)*L(i,2)/(1-L(i,4)^2*L(i,2)/L(i,1));
    Q(2,1,i)=Q(1,2,i);
    Q(2,2,i)=L(i,2)/(1-L(i,4)^2*L(i,2)/L(i,1));
    Q(3,3,i)=L(i,3);
end;

%% Transform to laminate global axis (x - 0 degrees , y - 90 degrees)
for i=1:M;
    theta=L(i,5)*pi/180;
    m=cos(theta);
    n=sin(theta);
   
% Transformation matrix T (for stresses)    
    T(1,1,i)=m^2;
    T(1,2,i)=n^2;
    T(1,3,i)=2*m*n;
    T(2,1,i)=n^2;
    T(2,2,i)=m^2;
    T(2,3,i)=-2*m*n;
    T(3,1,i)=-m*n;
    T(3,2,i)=m*n;
    T(3,3,i)=m^2-n^2;
   % Ttest=(T(:,:,i))
% Invert T
    Tii=inv(T(:,:,i));
    Ti(:,:,i)=Tii;
    
% Transformation matrix Tg (for engineering strains)   
    Tg(:,:,i)=Tii';
    
% Stiffness an Compliance matrix in laminate global axis   
    TTg=Tg(:,:,i);
    QQbb=Q(:,:,i);
    QQb=Tii*QQbb*TTg;
    Qb(:,:,i)=QQb;
end;

%%% Compute Middle plane location
TH=0;
z=zeros(M,1);
th=zeros(M,1);
for i=1:M;
    TH=TH+L(i,6);
end;
th_mid=TH/2;

%%% Compute location and thickness of laminae
TH=0;
for i=1:M;
    th(i)=L(i,6);
    TH=TH+th(i);
    z(i)=TH-th(i)/2-th_mid;
end;

%%% Compute A,B,D matrices 
A=zeros(3,3);
B=zeros(3,3);
D=zeros(3,3);
for i=1:M;
    for j=1:3;
        for k=1:3;
            A(j,k)=A(j,k)+Qb(j,k,i)*(th(i));
            B(j,k)=B(j,k)+Qb(j,k,i)*(th(i)*z(i));
            D(j,k)=D(j,k)+Qb(j,k,i)*(th(i)*z(i)^2+th(i)^3/12);
        end;
    end;
end;


ABD =[A B;B D]

N = [Nx ; Ny ; Nxy; Mx; My; Mxy];

% Evaluate the midsurface strains(eo) and change in curvature (ko)
Eo= inv(ABD)*N;
eo= Eo(1:3)';
ko= Eo(4:6)';

e_inter=zeros(3,M*2)
s_inter=zeros(3,M*2)
% Evaluate strain at the top and at the bottom of each layer
for i=1:M;
    for j=1:3;
        e_inter(j,i*2-1)=eo(j)+ko(j)*(z(i)-th(i)/2);
        e_inter(j,i*2)=eo(j)+ko(j)*(z(i)+th(i)/2);
    end;
end;

% Evaluate stress at the top and at the bottom of each layer
for i=1:M;
        s_inter(:,i*2-1)=Qb(:,:,i)*e_inter(:,i*2-1);
        s_inter(:,i*2)=Qb(:,:,i)*e_inter(:,i);
end;
