function [eps_x_z1,eps_x_z2,eps_y_z1,eps_y_z2,g_xy_z1,g_xy_z2,g_xz_z1,g_xz_z2,g_yz_z1,g_yz_z2,...
    eps_n_z1,eps_n_z2,eps_glo_x_z1,eps_glo_x_z2,eps_glo_y_z1,eps_glo_y_z2,eps_glo_z_z1,...
    eps_glo_z_z2,g_glo_xy_z1,g_glo_xy_z2,g_glo_xz_z1,g_glo_xz_z2,g_glo_yz_z1,g_glo_yz_z2,except,n_nodi]=post_strain(u,v,w,tx,ty,tz,nEle,nNode,...
                                                         Cne,Pnn,gcoord_index,Off,mcoord,t,desired)

%-------------------------------------------------------------------------%
%    FUNCTION FOR COMPUTING THE STRAINS FROM THE DISPLACEMENTS FOR        %
%    TRIA AND QUAD SHELL ELEMENTS                                         %
%-------------------------------------------------------------------------%
% INPUTS
% u,v,w = diplacement vector along x, y, z in th eglobal coord frame;
% tx,ty,tz = rotations vectors along x,y,z in the global ccord frame;
% nEle = number of elements;
% nNode = number of nodes;
% Cne = connectivity matrix;
% Pnn = nodal coordinates matrix in the global coord frame;
% gcoord_index = mapping vector to the indexes of the nodes;
% Off = vector of the offsets of the elements;
% mcoord = 3x3xnEle matrix of the direction cosines of the local coord of 
%          the elements;
% t = vector of the thicknesses of the elements;
% desired = matrix of the direction cosines of the coord system for the
%           computation of the nodal strains;

% OUTPUTS
% eps_x_z1,eps_x_z2,eps_y_z1,eps_y_z2,g_xy_z1,g_xy_z2,g_xz_z1,g_xz_z2,...
% g_yz_z1,g_yz_z2 = strain comonents eps_x, eps_y, gamma_xy, gamma_xz,
%                   gamma_yz expressed in the local coord frame of the
%                   elemets, computed on the top (z2) and bottom (z1)
%                   surface of the element. For each varaible each row of
%                   the cell contains the relative strain computed in the 
%                   following locations:
%                       TRIA ELEMENTS
%                   eps = {eps(CENTROID)}
%                       QUAD elements
%                   eps = {eps(CENTROID) eps(node_1) eps(node_2) ...
%                          eps(node_3) eps(node_4)}
% 
% eps_n_z1,eps_n_z2 = [nEle x 6]matrices of the nodal strain componets, 
%                     relative to the bottom (z1) and top (z2) surface of 
%                     the elements. They are expressed in the reference 
%                     frame obtained projecting the x axis of the frame 
%                     defined by the variable 'desired' on the plane of 
%                     the element. Each row of the matrices includes 
%                     the following components:
%                     [eps_x eps_y eps_z gamma_xy gamma_xz gamma_yz]
%
% eps_glo_x_z1,eps_glo_x_z2,..
% eps_glo_y_z1,eps_glo_y_z2,..
% eps_glo_z_z1,eps_glo_z_z2,..
% g_glo_xy_z1,g_glo_xy_z2,...
% g_glo_xz_z1,g_glo_xz_z2,..
% g_glo_yz_z1,g_glo_yz_z2 = strain comonents eps_x, eps_y, eps_z, gamma_xy, 
%                           gamma_xz, gamma_yz expressed in the coord frame 
%                           obtained by projecting the x axis of the frame 
%                           defined by the variable 'desired' on the plane
%                           of the element. For each varaible each row of
%                           the cell contains the relative strain computed 
%                           in the following locations:
%                               TRIA ELEMENTS
%                           eps = {eps(CENTROID)}
%                               QUAD elements
%                           eps = {eps(CENTROID) eps(node_1) eps(node_2)...
%                                 eps(node_3) eps(node_4)}
%
% THE FOLLOWING OUTPUTS ARE MAINLY FOR THE EVENTAUL COMPUTATION OF THE 
% STRESSES TROUGH THE FUNCTION 'post_stress'
%
% except = [nEle x 1] vector conating 1 for the elements where the x axis of the
%          'desired' coord frame is perpendicular to the surface of the 
%           element and 0 when that doesn't occur.
% 
% n_nodi = [nEle x 1] vector containing the dividend for each element for
%          the computation of the avarage strain value in the nodes


% INIZIALIZING ALL THE VARIABLES
eps_x_z2=cell(nEle,1);
eps_y_z2=cell(nEle,1);
g_xy_z2	=cell(nEle,1);
g_xz_z2	=cell(nEle,1);
g_yz_z2	=cell(nEle,1);
eps_x_z1=cell(nEle,1);
eps_y_z1=cell(nEle,1);
g_xy_z1	=cell(nEle,1);
g_xz_z1	=cell(nEle,1);
g_yz_z1	=cell(nEle,1);

eps_glo_x_z2=cell(nEle,1);
eps_glo_y_z2=cell(nEle,1);
eps_glo_z_z2=cell(nEle,1);
g_glo_xy_z2 =cell(nEle,1);
g_glo_xz_z2 =cell(nEle,1);
g_glo_yz_z2 =cell(nEle,1);


eps_glo_x_z1=cell(nEle,1);
eps_glo_y_z1=cell(nEle,1);
eps_glo_z_z1=cell(nEle,1);
g_glo_xy_z1 =cell(nEle,1);
g_glo_xz_z1 =cell(nEle,1);
g_glo_yz_z1 =cell(nEle,1);

eps_n_x_z1=zeros(nNode,1);
eps_n_y_z1=zeros(nNode,1);
eps_n_z_z1=zeros(nNode,1);
g_n_xy_z1=zeros(nNode,1);
g_n_xz_z1=zeros(nNode,1);
g_n_yz_z1=zeros(nNode,1);
eps_n_x_z2=zeros(nNode,1);
eps_n_y_z2=zeros(nNode,1);
eps_n_z_z2=zeros(nNode,1);
g_n_xy_z2=zeros(nNode,1);
g_n_xz_z2=zeros(nNode,1);
g_n_yz_z2=zeros(nNode,1);

n_nodi=zeros(nNode,1);

except=zeros(nEle,1);

% COMPUTING THE MATRIX FOR THE BILINEAR EXTRAPOLATION OF THE STRAINS FROM 
% THE GAUSS POINTS TO THE NODES OF THE ELEMENTS

coef=0.5*sqrt(3);

extrap=[1+coef -0.5 1-coef -0.5
        -0.5 1+coef -0.5 1-coef
        1-coef -0.5 1+coef -0.5
        -0.5 1-coef -0.5 1+coef];

% ATTRIBUTE THE DEFAULT VALUE TO desired IN CASE IT IS EMPTY
if isempty(desired)
        desired=[1 0 0
         0 1 0
         0 0 1];
end

%##### COMPUTING THE STRAINS FOR EACH ELEMENT
for e=1:nEle

    % DEFINING THE #NODES AND OF STRAIN VALUES FOR EACH ELEMENT DEPENDING
    % ON IF IT IS A TRIA OR A QUAD
    if Cne(e,end)==0
        elm_nodes=3;
        elm_sect=1;
        pen=1;
    else
        elm_nodes=4;
        elm_sect=5;
        pen=0;
    end
    
    % DEFINING THE VECTOR OF THE DOFs FOR EACH ELEMENT
    q(1:6:elm_nodes*6,1)=u(gcoord_index(Cne(e,2:elm_nodes+1)));
    q(2:6:elm_nodes*6,1)=v(gcoord_index(Cne(e,2:elm_nodes+1))); 
    q(3:6:elm_nodes*6,1)=w(gcoord_index(Cne(e,2:elm_nodes+1))); 
    q(4:6:elm_nodes*6,1)=tx(gcoord_index(Cne(e,2:elm_nodes+1))); 
    q(5:6:elm_nodes*6,1)=ty(gcoord_index(Cne(e,2:elm_nodes+1))); 
    q(6:6:elm_nodes*6,1)=tz(gcoord_index(Cne(e,2:elm_nodes+1)));
    
    
    %DEFINING THE GLOBAL COORDINATES OF THE NODES AND THE ROTATION MATIX
    % FOR EACH ELEMENT
    glocoord = Pnn(gcoord_index(Cne(e,2:elm_nodes+1)),:);
    Rot = mcoord(:,:,e);
    
    % DEFINING THE TOP (Z2) AND BOTTOM (Z1) SURFACES OF THE ELEMENT
    z2=t(e);
    z1=-t(e);

    % TRANSFORMING THE GLOBAL COORDINATES TO THE LOCAL ONES
    locCoord = glocoord * Rot';

    X = locCoord(:,1);
    Y = locCoord(:,2);

    % TRANSFORMING THE GLOBAL DOFs INTO THE  LOCAL DOFs
    Rot_q=zeros(elm_nodes*6,elm_nodes*6);
    Rot_q(1:3,1:3)=Rot;
    Rot_q(4:6,4:6)=Rot;
    Rot_q(7:9,7:9)=Rot;
    Rot_q(10:12,10:12)=Rot;
    Rot_q(13:15,13:15)=Rot;
    Rot_q(16:18,16:18)=Rot;
    Rot_q(19:elm_nodes*6-3,19:elm_nodes*6-3)=Rot(1:elm_nodes*6-21,1:elm_nodes*6-21);
    Rot_q(22:elm_nodes*6,22:elm_nodes*6)=Rot(1:elm_nodes*6-21,1:elm_nodes*6-21);

    q_loc=Rot_q*q;
    
    % TRANSFORMING THE LOCAL DISPLACEMENTS ACCORDING TO THE OFFSET OF THE
    % ELEMENT
    q_loc(1:6:end,1)=q_loc(1:6:end,1)+Off(e)*q_loc(5:6:end,1);
    q_loc(2:6:end,1)=q_loc(2:6:end,1)-Off(e)*q_loc(4:6:end,1);
    
    %## COMPUTING THE STRAINS IN THE LOCAL COORD SYSTEM
    if Cne(e,end)==0

        %## COMPUTING THE STRAINS FOR TRIA ELEMENTS

        % COMPUTING THE DERIVATIVES OF THE SHAPE FUNCTIONS
        [Bm,Bb,Bs]=compute_B_tria(X,Y,1,1,1/3,1/3);
        
        % COMPUTING THE MEMBRANE AND BENDING STRAINS ON THE TOP AND BOOTOM
        % SURFACE
        memb_centr=Bm*q_loc([1:5 7:11 13:17]);
        bend_centr=Bb*q_loc([1:5 7:11 13:17]);
        epsb_centr_z2=memb_centr+z2*bend_centr;
        epsb_centr_z1=memb_centr+z1*bend_centr;
        epss_centr_z2=Bs*q_loc([1:5 7:11 13:17]);
        eps_centr_z2=[epsb_centr_z2
                  epss_centr_z2];
        eps_centr_z1=[epsb_centr_z1
                  epss_centr_z2];
        
        % DEFINING THE STRAIN COMPONENTS IN THE CENTROID OF THE TRIA
        % ON THE TOP AND BOOTOM SURFACE
        eps_x_z2{e}=[eps_centr_z2(1)];
        eps_y_z2{e}=[eps_centr_z2(2)];
        g_xy_z2{e}=[eps_centr_z2(3) ];
        g_xz_z2{e}=[eps_centr_z2(4) ];
        g_yz_z2{e}=[eps_centr_z2(5) ];
        eps_x_z1{e}=[eps_centr_z1(1)];
        eps_y_z1{e}=[eps_centr_z1(2)];
        g_xy_z1{e}=[eps_centr_z1(3) ];
        g_xz_z1{e}=[eps_centr_z1(4) ];
        g_yz_z1{e}=[eps_centr_z1(5) ];
    
    else
        %## COMPUTING THE STRAINS FOR THE QUAD ELEMENTS

        % COMPUTING THE DERIVATIVES OF THE SHAPE FUNCTIONS
        [Bm,Bb,Bs]=compute_B_quad(X,Y,3,2,-1,1);

        % COMPUTING THE MEMBRANE AND BENDING STRAINS ON THE TOP (Z2) AND 
        % BOOTOM (Z1) SURFACE
        memb_g1=Bm(1:5:15,:)*q_loc;
        bend_g1=Bb(1:5:15,:)*q_loc;
        epsb_g1_z2=memb_g1+z2*bend_g1;
        epsb_g1_z1=memb_g1+z1*bend_g1;
        epss_g1_z2=Bs(1:5:10,:)*q_loc;
        eps_g1_z2=[epsb_g1_z2
                  epss_g1_z2];
        eps_g1_z1=[epsb_g1_z1
                  epss_g1_z2];
        memb_g2=Bm(3:5:15,:)*q_loc;
        bend_g2=Bb(3:5:15,:)*q_loc;
        epsb_g2_z2=memb_g2+z2*bend_g2;
        epsb_g2_z1=memb_g2+z1*bend_g2;
        epss_g2_z2=Bs(3:5:10,:)*q_loc;
        eps_g2_z2=[epsb_g2_z2
                  epss_g2_z2];
        eps_g2_z1=[epsb_g2_z1
                  epss_g2_z2];
        memb_g3=Bm(4:5:15,:)*q_loc;
        bend_g3=Bb(4:5:15,:)*q_loc;
        epsb_g3_z2=memb_g3+z2*bend_g3;
        epsb_g3_z1=memb_g3+z1*bend_g3;
        epss_g3_z2=Bs(4:5:10,:)*q_loc;
        eps_g3_z2=[epsb_g3_z2
                  epss_g3_z2];
        eps_g3_z1=[epsb_g3_z1
                  epss_g3_z2];
        memb_g4=Bm(2:5:15,:)*q_loc;
        bend_g4=Bb(2:5:15,:)*q_loc;
        epsb_g4_z2=memb_g4+z2*bend_g4;
        epsb_g4_z1=memb_g4+z1*bend_g4;
        epss_g4_z2=Bs(2:5:10,:)*q_loc;
        eps_g4_z2=[epsb_g4_z2
                  epss_g4_z2];
        eps_g4_z1=[epsb_g4_z1
                  epss_g4_z2];
        memb_centr=Bm(5:5:15,:)*q_loc;
        bend_centr=Bb(5:5:15,:)*q_loc;
        epsb_centr_z2=memb_centr+z2*bend_centr;
        epsb_centr_z1=memb_centr+z1*bend_centr;
        epss_centr_z2=Bs(5:5:10,:)*q_loc;
        eps_centr_z2=[epsb_centr_z2
                  epss_centr_z2];
        eps_centr_z1=[epsb_centr_z1
                  epss_centr_z2];
    
        % DEFINING THE STRAIN COMPONENTS IN THE CENTROID AND IN THE NODES
        % OF THE QUAD, ON THE TOP (Z2) AND BOOTOM (Z1) SURFACE
        eps_x_z2{e}=[eps_centr_z2(1) (extrap*[eps_g1_z2(1);eps_g2_z2(1);eps_g3_z2(1);eps_g4_z2(1)])'];
        eps_y_z2{e}=[eps_centr_z2(2) (extrap*[eps_g1_z2(2);eps_g2_z2(2);eps_g3_z2(2);eps_g4_z2(2)])'];
        g_xy_z2{e}=[eps_centr_z2(3) (extrap*[eps_g1_z2(3);eps_g2_z2(3);eps_g3_z2(3);eps_g4_z2(3)])'];
        g_xz_z2{e}=[eps_centr_z2(4) (extrap*[eps_g1_z2(4);eps_g2_z2(4);eps_g3_z2(4);eps_g4_z2(4)])'];
        g_yz_z2{e}=[eps_centr_z2(5) (extrap*[eps_g1_z2(5);eps_g2_z2(5);eps_g3_z2(5);eps_g4_z2(5)])'];
        eps_x_z1{e}=[eps_centr_z1(1) (extrap*[eps_g1_z1(1);eps_g2_z1(1);eps_g3_z1(1);eps_g4_z1(1)])'];
        eps_y_z1{e}=[eps_centr_z1(2) (extrap*[eps_g1_z1(2);eps_g2_z1(2);eps_g3_z1(2);eps_g4_z1(2)])'];
        g_xy_z1{e}=[eps_centr_z1(3) (extrap*[eps_g1_z1(3);eps_g2_z1(3);eps_g3_z1(3);eps_g4_z1(3)])'];
        g_xz_z1{e}=[eps_centr_z1(4) (extrap*[eps_g1_z1(4);eps_g2_z1(4);eps_g3_z1(4);eps_g4_z1(4)])'];
        g_yz_z1{e}=[eps_centr_z1(5) (extrap*[eps_g1_z1(5);eps_g2_z1(5);eps_g3_z1(5);eps_g4_z1(5)])'];
    end

    

    % TRANSFORMING THE STRAINS FROM THE LOCAL COORD FRAME TO THE DESIRED
    % ONE
    for j=1:elm_sect  

        [transf,except(e)]=project_coord_v2(Rot,desired);

        eps_glo_z2=transf*[eps_x_z2{e}(j) g_xy_z2{e}(j)/2 g_xz_z2{e}(j)/2
                         g_xy_z2{e}(j)/2 eps_y_z2{e}(j) g_yz_z2{e}(j)/2
                         g_xz_z2{e}(j)/2 g_yz_z2{e}(j)/2 0]*transf';
        eps_glo_z1=transf*[eps_x_z1{e}(j) g_xy_z1{e}(j)/2 g_xz_z1{e}(j)/2
                         g_xy_z1{e}(j)/2 eps_y_z1{e}(j) g_yz_z1{e}(j)/2
                         g_xz_z1{e}(j)/2 g_yz_z1{e}(j)/2 0]*transf';

        eps_glo_x_z2{e}(j)=eps_glo_z2(1,1);
        eps_glo_y_z2{e}(j)=eps_glo_z2(2,2);
        eps_glo_z_z2{e}(j)=eps_glo_z2(3,3);
        g_glo_xy_z2{e}(j)=eps_glo_z2(1,2)*2;
        g_glo_xz_z2{e}(j)=eps_glo_z2(1,3)*2;
        g_glo_yz_z2{e}(j)=eps_glo_z2(2,3)*2;
        
         eps_glo_x_z1{e}(j)=eps_glo_z1(1,1);
         eps_glo_y_z1{e}(j)=eps_glo_z1(2,2);
         eps_glo_z_z1{e}(j)=eps_glo_z1(3,3);
         g_glo_xy_z1{e}(j)=eps_glo_z1(1,2)*2;
         g_glo_xz_z1{e}(j)=eps_glo_z1(1,3)*2;
         g_glo_yz_z1{e}(j)=eps_glo_z1(2,3)*2;        
        
    end
    
    % COMPUTING THE NODAL VALUES OF THE STRAINS IN THE DESIRED REFERENCE
    % FRAME
    eps_n_x_z1(gcoord_index(Cne(e,2:elm_nodes+1),:),1)=eps_n_x_z1(gcoord_index(Cne(e,2:elm_nodes+1),:),1)+eps_glo_x_z1{e}(2-pen:elm_sect)';
    eps_n_y_z1(gcoord_index(Cne(e,2:elm_nodes+1),:),1)=eps_n_y_z1(gcoord_index(Cne(e,2:elm_nodes+1),:),1)+eps_glo_y_z1{e}(2-pen:elm_sect)';
    eps_n_z_z1(gcoord_index(Cne(e,2:elm_nodes+1),:),1)=eps_n_z_z1(gcoord_index(Cne(e,2:elm_nodes+1),:),1)+eps_glo_z_z1{e}(2-pen:elm_sect)';
    g_n_xy_z1(gcoord_index(Cne(e,2:elm_nodes+1),:),1)=g_n_xy_z1(gcoord_index(Cne(e,2:elm_nodes+1),:),1)+g_glo_xy_z1{e}(2-pen:elm_sect)';
    g_n_xz_z1(gcoord_index(Cne(e,2:elm_nodes+1),:),1)=g_n_xz_z1(gcoord_index(Cne(e,2:elm_nodes+1),:),1)+g_glo_xz_z1{e}(2-pen:elm_sect)';
    g_n_yz_z1(gcoord_index(Cne(e,2:elm_nodes+1),:),1)=g_n_yz_z1(gcoord_index(Cne(e,2:elm_nodes+1),:),1)+g_glo_yz_z1{e}(2-pen:elm_sect)';

    eps_n_x_z2(gcoord_index(Cne(e,2:elm_nodes+1),:),1)=eps_n_x_z2(gcoord_index(Cne(e,2:elm_nodes+1),:),1)+eps_glo_x_z2{e}(2-pen:elm_sect)';
    eps_n_y_z2(gcoord_index(Cne(e,2:elm_nodes+1),:),1)=eps_n_y_z2(gcoord_index(Cne(e,2:elm_nodes+1),:),1)+eps_glo_y_z2{e}(2-pen:elm_sect)';
    eps_n_z_z2(gcoord_index(Cne(e,2:elm_nodes+1),:),1)=eps_n_z_z2(gcoord_index(Cne(e,2:elm_nodes+1),:),1)+eps_glo_z_z2{e}(2-pen:elm_sect)';
    g_n_xy_z2(gcoord_index(Cne(e,2:elm_nodes+1),:),1)=g_n_xy_z2(gcoord_index(Cne(e,2:elm_nodes+1),:),1)+g_glo_xy_z2{e}(2-pen:elm_sect)';
    g_n_xz_z2(gcoord_index(Cne(e,2:elm_nodes+1),:),1)=g_n_xz_z2(gcoord_index(Cne(e,2:elm_nodes+1),:),1)+g_glo_xz_z2{e}(2-pen:elm_sect)';
    g_n_yz_z2(gcoord_index(Cne(e,2:elm_nodes+1),:),1)=g_n_yz_z2(gcoord_index(Cne(e,2:elm_nodes+1),:),1)+g_glo_yz_z2{e}(2-pen:elm_sect)';

    
    n_nodi(gcoord_index(Cne(e,2:elm_nodes+1),:),1)=n_nodi(gcoord_index(Cne(e,2:elm_nodes+1),:),1)+1;

end

% AVAREGING THE NODAL VALUES OF THE STRAINS 
eps_n_x_z1=eps_n_x_z1./n_nodi;
eps_n_y_z1=eps_n_y_z1./n_nodi;
eps_n_z_z1=eps_n_z_z1./n_nodi;
g_n_xy_z1=g_n_xy_z1./n_nodi;
g_n_xz_z1=g_n_xz_z1./n_nodi;
g_n_yz_z1=g_n_yz_z1./n_nodi;

eps_n_z1=[eps_n_x_z1 eps_n_y_z1 eps_n_z_z1 g_n_xy_z1 g_n_xz_z1 g_n_yz_z1];

eps_n_x_z2=eps_n_x_z2./n_nodi;
eps_n_y_z2=eps_n_y_z2./n_nodi;
eps_n_z_z2=eps_n_z_z2./n_nodi;
g_n_xy_z2=g_n_xy_z2./n_nodi;
g_n_xz_z2=g_n_xz_z2./n_nodi;
g_n_yz_z2=g_n_yz_z2./n_nodi;

eps_n_z2=[eps_n_x_z2 eps_n_y_z2 eps_n_z_z2 g_n_xy_z2 g_n_xz_z2 g_n_yz_z2];
end