%nondc.m (displacement control)
%by Louie L. Yaw - Walla Walla University - 5-8-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 "displacement" control type analysis
%Used to demonstrate how displacement control succeeds in problems with
%limit points and snap through.

%Clear all variables
%clear

%%%%User Input
%number of increments
ninc=9;
%Max final displacement
Dmax=9.0;
%Input Variables for Newton-Raphson iterations
tol=1e-3;%tolerance for equilibrium set by user
maxiter=10;%set maximum number of iterations
%%%%End of User Input

%%%%Start of Implicit Analysis
%Initialize some variables
F=1.0; %Really just a dummy force scaled later to accomplish the displacements
Fnp1=0;%The current value of the total external applied forces
lambda=0;%The current load ratio
%The incremental displacement in each step
duq=Dmax/ninc;

%The current total global displacement
ug=0;

%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
    uhat=(1/k)*F; %Calculate displacement due to F using current stiffness
    dlambda=duq/uhat; %Calculate the incremental change needed in the load ratio
    lambda=lambda+dlambda; %update the load ratio
    dF=dlambda*F; %Calculate the load increment    
    du=dF/k;%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
    
    %%%%%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
    dlam=0; %iterated Load ratio correction
    Fint=ug^3-12*ug^2+32*ug; %Calculate the current internal force
    Fext=lambda*F;%Set variable Fext equal to the current external force
    Resid=Fext-Fint;%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
        %Steps to calculate load ratio correction dlam
        umlat=(1/k)*Resid;
        uhat=(1/k)*F;
        dlam=dlam-umlat/uhat;
        %Calculate correction to ug
        dug=dug+(1/k)*(Resid-(umlat/uhat)*F);
        
        Fint=(ug+dug)^3-12*(ug+dug)^2+32*(ug+dug); %Update Fint based on ug+dug
        Resid=(lambda+dlam)*F-Fint; %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
    lambda=lambda+dlam;
    ug=ug+dug; 
    %End Newton-Raphson iterations stuff
    
    ugp(i+1)=ug;%Store the total displacement at the end of increment i
    fp(i+1)=lambda*F;%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=lambda*F;
%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:Dmax];
%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)
