%nonorthog.m
%by Louie L. Yaw - Walla Walla University - 6-21-09
%Example nonlinear implicit analysis of a single dof system, bar in tension
%Demonstration of an incremental IMPLICIT static analysis
%Assume function f(u) with a limit point (like snap through)
%f(u)=u^3-12u^2+32u  the roots are at +4 and +8 and 0
%Therefore k(u)=df/du=3u^2-24u+32
%Equal increments are used throughout.
%This is a arc length control/orthogonal trajectory corrector type analysis
%started 6-21-09 (not working)

%Clear all variables
%clear

%%%%User Input
mk=[0 0];
%number of increments
ninc=47;
%Arc Length parameter
s=3.3;
%"Proportional" external applied load
P=1;
%Input Variables for Newton-Raphson iterations
tol=1e-2;%tolerance for equilibrium set by user
maxiter=20;%set maximum number of iterations
%%%%End of User Input

%%%%Start of Implicit Analysis
%The current value of the total external applied forces
fg=0;
%The current total global displacement
ug=0;
%Note that rdot=dr/dlambda=P (actually the negative of the load)
%P=-P;
%Storage arrays for plotting later
ugp=zeros(ninc+1,1);%total global displacements are recorded in this vector
fp=zeros(ninc+1,1);%total global forces are recorded in this vector
ugp(1)=0;
fp(1)=0;

%Loop over the load increments
for i=1:ninc;
    i %print to screen the current load increment
    k=3*ug^2-24*ug+32; %calculate the current stiffness, df/du
    %The initial incremental force in each load step
    df=s*k/(k^2+P^2)^0.5;
    du=P*s/(k^2+P^2)^0.5;%calculate the incremental displacement for load step i
    ug=ug+du;%calculate the current total displacement
    fg=fg+df;%Calculate the current total external force factor
    
    %%%%%Start of statements to include Newton-Raphson iterations
    %Initialize some variables for the Newton-Raphson iterations
    kiter=0;%iteration variable
    dug=0; %Variable of small delta correction to ug to achieve equilibrium
    dfg=0; %Variable of small delta correction to df
    Fint=ug^3-12*ug^2+32*ug; %Calculate the current internal force
    Fext=fg*P;%Set variable Fext equal to the current external force
    Resid=Fint-Fext;%Calculate the residual
    R=sqrt(Resid^2)%Calculate norm of the residual
    
    %Check the residual and max iterations limit, 
    %while above tolerance and max iterations not exceeded do iterations
    while R>tol && kiter<maxiter;
        k=3*(ug+dug)^2-24*(ug+dug)+32; %calculate current stiffness with latest ug+dug
%         dfg=dfg-Resid*P/(k^2+P^2); %correct the value of dfg        
%         dug=dug-Resid*k/(k^2+P^2); %Correct the value of dug based on the residual
        dfg=dfg-Resid*(-P)/(k^2+(-P)^2); %correct the value of dfg        
        dug=dug-Resid*k/(k^2+(-P)^2); %Correct the value of dug based on the residual
        
        Fint=(ug+dug)^3-12*(ug+dug)^2+32*(ug+dug); %Update Fint based on ug+dug
        Fext=(fg+dfg)*P;
        Resid=Fint-Fext; %Update the residual
        R=sqrt(Resid^2) %Update norm of the residual
        kiter=kiter+1; %Update the number of iterations used to achieve convergence
        %If convergence not achieved and maxiter reached, print a warning
        if kiter==maxiter && R>tol;
            display('WARNING -- Convergence criteria not met')
        end;
    end;
    
    %Output number of iterations required to achieve convergence
    kiter
    if mk(1,1)<kiter;
        mk(1,1)=kiter;
        mk(1,2)=i;
    end;        
    ug=ug+dug; 
    fg=fg+dfg;
    %End Newton-Raphson iterations stuff
    
    ugp(i+1)=ug;%Store the total displacement at the end of increment i
    fp(i+1)=fg;%Store the total external force applied at end of increment i
end;
%End of Implicit analysis

%The final displacement is
ug;
%The total external force is, (or Fext=fg)
Fext=fg;
%The internal force is
Fint=ug^3-12*ug^2+32*ug;

%Plot the explicit and exact force versus displacement solutions
ue=[0:0.05:9];
%ue=[0:0.1:ceil(max(ugp))]'; %create displacement points for plotting the exact solution
fe=ue.^3-12.*ue.^2+32*ue; %create the array of force data points for plotting exact

%Plot Explicit and exact solutions for comparison
plot(ugp,fp,'bo-',ue,fe,'k','LineWidth',2.0)
ylabel('Force, F')
xlabel('Displacement, u')
grid on
legend('Implicit','Exact',2)

mk
