% FEAcomp. Solution of problem 3.x. 
%
% J.A. Mayugo, febraury 2006

%%% Compute the FC in 3D 

clear all;

%%% In this program version the material data is defined in these file directelly

%%Stress state

%           SX          SY          SZ          SXY         SYZ         SXZ     
str(1,:)=[ 26.391    -0.13551E-01  0.0000    -0.28242E-14 0.29466E-16 0.26472E-15]
str(2,:)=[-7.7621      1.3850      0.0000     0.25357E-14 -0.29088E-15 0.78807E-16]
str(3,:)=[ 9.3143     0.68572      0.0000     -1.1122     -0.21245E-15 0.47077E-15]
str(4,:)=[ 9.3143     0.68572      0.0000      1.1122     0.30339E-15 0.33374E-15]

[nlayer,N]=size(str);

%Failure criteria parameters definition
F1t=1830;       %! F1t strength
F1c=1096;       %! F1c strength
F2t=57;         %! F2t strength
F2c=228;        %! F2c strength
F3t=1e6;        %! F3t=F3c strength (large value, not compute)
F3c=F3t;        %! F3t=F3c strength (large value, not compute)
F6=71;          %! F6 strength
F4=1e6;         %! F4 strength (large value, not compute)
F5=1e6;         %! F5 strength (large value, not compute)
c4=-1;          %! Tsai-Wu coefficients
c5=-1;
c6=-1;

%% Maximum Stress Criterion
for ilay=1:nlayer
    if str(ilay,1)>=0.0
        I(ilay,1)=str(ilay,1)/F1t;
    else
        I(ilay,1)=-str(ilay,1)/F1c;
    end
    if str(ilay,2)>=0.0
        I(ilay,2)=str(ilay,2)/F2t;
    else
        I(ilay,2)=-str(ilay,2)/F2c;
    end
    if str(ilay,3)>=0.0
        I(ilay,3)=str(ilay,3)/F3t;
    else
        I(ilay,3)=-str(ilay,3)/F3c;
    end
    I(ilay,6)=abs(str(ilay,6))/F4;
    I(ilay,5)=abs(str(ilay,5))/F5;
    I(ilay,4)=abs(str(ilay,4))/F6;

    I;
%   Failure index  
    If_ms(ilay)=max(I(ilay,:));
%   Strength ratio
    R_ms(ilay)=1/If_ms(ilay);
end;

If_ms
R_ms

%% Tsai-Wu Criterion
for ilay=1:nlayer
    A(ilay)=str(ilay,1)^2/(F1t*F1c)+str(ilay,2)^2/(F2t*F2c)+str(ilay,3)^2/(F3t*F3c)+ ...
            str(ilay,4)^2/(F6)^2+str(ilay,5)^2/(F5)^2+str(ilay,6)^2/(F4)^2+ ...
            c4*str(ilay,2)*str(ilay,3)/sqrt(F2t*F2c*F3t*F3c)+ ...
            c5*str(ilay,1)*str(ilay,3)/sqrt(F1t*F1c*F3t*F3c)+ ...
            c6*str(ilay,1)*str(ilay,2)/sqrt(F1t*F1c*F2t*F2c);
    B(ilay)=(1/F1t-1/F1c)*str(ilay,1)+(1/F2t-1/F2c)*str(ilay,2)+(1/F3t-1/F3c)*str(ilay,3);
%    Iff_tw(ilay)=A(ilay)+B(ilay);  % "Strength index", not recommended 

%   Strength ratio
    R_tw(ilay)=-B(ilay)/(2*A(ilay))+sqrt((B(ilay)/2/A(ilay))^2+1/A(ilay));
%   Failure index    
    If_tw(ilay)=1/R_tw(ilay);    
end;
A
B
If_tw
R_tw