clear
clc

% Solution In-class Exercise Week #10: Sectional analysis

%% Material properties
E=200000;
fy=355;
epsiY=fy/E;

%% Geometric properties
bf=200;
tf=12;
d=300;
tw=10;

Wel=(bf*d^2)/6-(bf-tw)*(d-2*tf)^3/(6*d);
Wpl=(bf*d^2)/4-(bf-tw)*(d-2*tf)^2/4;

%% Hand calculation
My=Wel*fy;
Mpl=Wpl*fy;

phiY=epsiY/(d/2);
phiPl=9999*phiY;

%% Discretizing section into fibers
N_fh=1;
N_fv=2;
N_wh=1;
N_wv=20;
[coordNf,Anf]=computeLocationNfFibers_WF(bf,tf,d,tw,N_fh,N_fv,N_wh,N_wv);
numFibers=size(coordNf,1);

%% Verify section properties
% I=0;
% for iFib=1:numFibers
%    I=I+coordNf(iFib,1)^2*Anf(iFib);
% end
% Wel_test=I/(d/2);
% verif_Wel=abs(Wel_test-Wel)/Wel*100;

%% Material law parameters
alphaS=0.03;
% alphaS=0.0;
Es=alphaS*E;
fu=100*fy;
epsiP=(fu-fy)/Es;
epsiPc=10;
epsiC=epsiP+epsiY;
epsiU=epsiPc+epsiC;
Epc=-fu/epsiPc;

%% Initialize parameters for material law
epsi_previous=0;
sigma_previous=0;
sigmaMaxPos=fy;
sigmaMaxNeg=-fy;
yieldFlag_Pos=0;
cappingFlag_Pos=0;
yieldFlag_Neg=0;
cappingFlag_Neg=0;
reversalFlag=0;
residualFlag=0;
epsi_sigma0_currentPos=0;
epsi_sigma0_currentNeg=0;
epsi_sigma0_projected=0;
Di_previous=0;

%Fill for each fiber
fiberMatProp=zeros(numFibers,14);
for iFib = 1:numFibers
    fiberMatProp(iFib, 1) = epsi_previous;
    fiberMatProp(iFib, 2) = sigma_previous;
    fiberMatProp(iFib, 3) = sigmaMaxPos;
    fiberMatProp(iFib, 4) = sigmaMaxNeg;
    fiberMatProp(iFib, 5) = yieldFlag_Pos;
    fiberMatProp(iFib, 6) = cappingFlag_Pos;
    fiberMatProp(iFib, 7) = yieldFlag_Neg;
    fiberMatProp(iFib, 8) = cappingFlag_Neg;
    fiberMatProp(iFib, 9) = reversalFlag;
    fiberMatProp(iFib, 10) = residualFlag;
    fiberMatProp(iFib, 11) = epsi_sigma0_currentPos;
    fiberMatProp(iFib, 12) = epsi_sigma0_currentNeg;
    fiberMatProp(iFib, 13) = epsi_sigma0_projected;
    fiberMatProp(iFib, 14) = Di_previous;
end

%% Test the material law
% nIncrements=500;
% strain=linspace(0,0.10,nIncrements);
% sigma_vector=zeros(size(strain));
% for i=1:nIncrements
%     epsi=strain(i);
%     [sigma,k,sigmaMaxPos,sigmaMaxNeg,yieldFlag_Pos,cappingFlag_Pos,yieldFlag_Neg,cappingFlag_Neg,reversalFlag,residualFlag,epsi_sigma0_currentPos,epsi_sigma0_currentNeg,epsi_sigma0_projected,Di]...
%         =computeTrilinearMatLaw(epsi,epsi_previous,sigma_previous,epsi_sigma0_currentPos,epsi_sigma0_currentNeg,epsi_sigma0_projected,E,Es,Epc,epsiY,epsiC,epsiU,fy,fu,sigmaMaxPos,sigmaMaxNeg,yieldFlag_Pos,cappingFlag_Pos,yieldFlag_Neg,cappingFlag_Neg,reversalFlag,residualFlag,Di_previous);
%     sigma_vector(i)=sigma;
%     
%     sigma_previous=sigma;
%     epsi_previous=epsi;
%     Di_previous=Di;
% end
% 
% % Plot the results
% h1=figure;
% plot(strain,sigma_vector)
% xlabel('$\varepsilon$ [-]','Interpreter','Latex')
% ylabel('$\sigma$ [MPa]','Interpreter','Latex')
% grid on
%
% % plot_settings_ASCE(h1)

%% Section analysis
phiMax=10*phiY;
nIncrements=1000;
allPhi=generateInputStrainHistory([0,phiMax], nIncrements);
allMoments=zeros(size(allPhi));

for n=1:nIncrements
    phi=allPhi(n);
    M=0;
    for iFib=1:numFibers
        y=coordNf(iFib,1);
        epsi=y*phi;
        
        % Unpack the material law parameters
        epsi_previous = fiberMatProp(iFib, 1);
        sigma_previous = fiberMatProp(iFib, 2);
        sigmaMaxPos = fiberMatProp(iFib, 3);
        sigmaMaxNeg = fiberMatProp(iFib, 4);
        yieldFlag_Pos = fiberMatProp(iFib, 5);
        cappingFlag_Pos = fiberMatProp(iFib, 6);
        yieldFlag_Neg = fiberMatProp(iFib, 7);
        cappingFlag_Neg = fiberMatProp(iFib, 8);
        reversalFlag = fiberMatProp(iFib, 9);
        residualFlag = fiberMatProp(iFib, 10);
        epsi_sigma0_currentPos = fiberMatProp(iFib, 11);
        epsi_sigma0_currentNeg = fiberMatProp(iFib, 12);
        epsi_sigma0_projected = fiberMatProp(iFib, 13);
        Di_previous = fiberMatProp(iFib, 14);
        
        % Constitutive law
        [sigma,k,sigmaMaxPos,sigmaMaxNeg,yieldFlag_Pos,cappingFlag_Pos,yieldFlag_Neg,cappingFlag_Neg,reversalFlag,residualFlag,epsi_sigma0_currentPos,epsi_sigma0_currentNeg,epsi_sigma0_projected,Di]...
            =computeTrilinearMatLaw(epsi,epsi_previous,sigma_previous,epsi_sigma0_currentPos,epsi_sigma0_currentNeg,epsi_sigma0_projected,E,Es,Epc,epsiY,epsiC,epsiU,fy,fu,sigmaMaxPos,sigmaMaxNeg,yieldFlag_Pos,cappingFlag_Pos,yieldFlag_Neg,cappingFlag_Neg,reversalFlag,residualFlag,Di_previous);
        
        % Store updated properties back into fiberMatProp
        fiberMatProp(iFib, 1) = epsi;
        fiberMatProp(iFib, 2) = sigma;
        fiberMatProp(iFib, 3) = sigmaMaxPos;
        fiberMatProp(iFib, 4) = sigmaMaxNeg;
        fiberMatProp(iFib, 5) = yieldFlag_Pos;
        fiberMatProp(iFib, 6) = cappingFlag_Pos;
        fiberMatProp(iFib, 7) = yieldFlag_Neg;
        fiberMatProp(iFib, 8) = cappingFlag_Neg;
        fiberMatProp(iFib, 9) = reversalFlag;
        fiberMatProp(iFib, 10) = residualFlag;
        fiberMatProp(iFib, 11) = epsi_sigma0_currentPos;
        fiberMatProp(iFib, 12) = epsi_sigma0_currentNeg;
        fiberMatProp(iFib, 13) = epsi_sigma0_projected;
        fiberMatProp(iFib, 14) = Di;
        
        % Compute the moment
        M=M+y*sigma*Anf(iFib);
    end
    allMoments(n)=M;
end


%% Export results
% writematrix([allPhi',allMoments'], 'results_week10_elasticPlastic.txt');
% writematrix([allPhi',allMoments'], 'results_week10_bilin003.txt');

%% Plot results
h1=figure;
plot([0,phiY,phiPl],[0,My,Mpl])
hold on
plot(allPhi,allMoments)
legend('Hand calc','Fibers','Location','best')
xlabel('Curvature \phi [-]')
ylabel('Moment M [kNm]')
xlim([0,phiMax])
grid on

% plot_settings_ASCE(h1)

%% Compare the results
% res_elasticPlastic=importdata('results_week10_elasticPlastic.txt');
% res_bilin003=importdata('results_week10_bilin003.txt');
% 
% h2=figure;
% plot(res_elasticPlastic(:,1),res_elasticPlastic(:,2)/10^6)
% hold on
% plot(res_bilin003(:,1),res_bilin003(:,2)/10^6)
% legend('Elastic-perfectly plastic','0.03% hardening','Location','best')
% xlabel('Curvature \phi [-]')
% ylabel('Moment M [kNm]')
% xlim([0,phiMax])
% grid on

% plot_settings_ASCE(h2)


test=1;


%% Function to compute the input strain history
function strain=generateInputStrainHistory(peaksVector, nbSteps)
%nbSteps is used  in maxmim excursion, other excursions are scaled proportionaly

% General informations
peaksVector=[0,peaksVector];
nbLoadExcursions=size(peaksVector,2)-1;

% Find maximum load excursion
strainDiffLoadExcursions=abs([0,diff(peaksVector)]);
maxStrainDifference=max(strainDiffLoadExcursions);

% Generate the points for each load excursions
strain=[];
for i=1:nbLoadExcursions
    nbStepsTemp=round(strainDiffLoadExcursions(i+1)/maxStrainDifference*nbSteps);
    startStrain=peaksVector(i);
    endStrain=peaksVector(i+1);
    epsi11Temp=linspace(startStrain,endStrain,nbStepsTemp);
    if i<nbLoadExcursions
        epsi11Temp=epsi11Temp(1:end-1); %Remove duplicate strain
    end
    strain=[strain,epsi11Temp];
end
end


%% Function to return the stress and tangent from material law
function [M,k,MmaxPos,MmaxNeg,yieldFlag_Pos,cappingFlag_Pos,yieldFlag_Neg,cappingFlag_Neg,reversalFlag,residualFlag,theta_M0_currentPos,theta_M0_currentNeg,theta_M0_projected,Di]...
    =computeTrilinearMatLaw(theta,theta_previous,M_previous,theta_M0_currentPos,theta_M0_currentNeg,theta_M0_projected,ke,ks,kpc,thetaY,thetaC,thetaU,My,Mu,MmaxPos,MmaxNeg,yieldFlag_Pos,cappingFlag_Pos,yieldFlag_Neg,cappingFlag_Neg,reversalFlag,residualFlag,Di_previous)
    % Find the direction of the current increment "Di": 1:if Ri is moving towards positive and -1: if Ri is moving negative
    if theta>=theta_previous
        Di=1;
	else
		Di=-1;
    end
    
    % Trial elastic moment
    M=M_previous+ke*(theta-theta_previous);
    
    % Quick return if failed
    if residualFlag==1
        M=0;
        k=0;
        return
    end
    
    %If inelastic and going to positive rotation
    if Di==1 && M>MmaxPos
        yieldFlag_Neg=0;
        cappingFlag_Neg=0;
        if yieldFlag_Pos==1 || cappingFlag_Pos==1
        else
            theta_M0_currentPos=theta_M0_projected;
        end
        
        if abs(theta_M0_currentPos-theta)>=thetaY && abs(theta_M0_currentPos-theta)<thetaC %hardening path
            yieldFlag_Pos=1;
            theta_M0_currentNeg=theta_M0_projected;
            M=Di*My+ks*(theta-theta_M0_currentPos-thetaY);
            MmaxNeg=-My;
        elseif abs(theta_M0_currentPos-theta)>=thetaC && abs(theta_M0_currentPos-theta)<thetaU %softening path
            yieldFlag_Pos=0;
            cappingFlag_Pos=1;
            M=Di*Mu+kpc*(theta-theta_M0_currentPos-thetaC);
        else %residual path
            residualFlag=1;
            M=0;
        end
    
        MmaxPos=M;
    end
    
    
    %If inelastic and going to negative rotation
    if Di==-1 && M<MmaxNeg
        yieldFlag_Pos=0;
        cappingFlag_Pos=0;
        if yieldFlag_Neg==1 || cappingFlag_Neg==1
        else
            theta_M0_currentNeg=theta_M0_projected;
        end
        
        if abs(theta_M0_currentNeg-theta)>=thetaY && abs(theta_M0_currentNeg-theta)<thetaC %hardening path
            yieldFlag_Neg=1;
            theta_M0_currentPos=theta_M0_projected;
            M=Di*My+ks*(theta-theta_M0_currentNeg+thetaY);
            MmaxPos=My;
        elseif abs(theta_M0_currentNeg-theta)>=thetaC && abs(theta_M0_currentNeg-theta)<thetaU
            yieldFlag_Neg=0;
            cappingFlag_Neg=1;
            M=Di*Mu+kpc*(theta-theta_M0_currentNeg+thetaC);
        else %failed
            residualFlag=1;
            M=0;
        end
    
        MmaxNeg=M;
    end
   
    
    % Check if unloading
    if Di_previous*Di<0 
       theta_M0_projected=theta_previous-M_previous/ke;
    end    
    
    % Compute the tangent stiffness
    k=(M-M_previous)/(theta-theta_previous);
    if theta==theta_previous
        k=0;
    elseif M==M_previous
        k=10^-6;
    end
end