%%%%%%%%%%%% Main Function for the assignment: optimality.m %%%%%%%%%%%%%%%
%%%%%%%%%%%%%%% Written by Aman P.V. on 15/03/09 for ME 256 Course in the
%%%% Indian Institute of Science, Bangalore %%%%%%%%%%%%%%%%% 
% Modify modifydata.m for different boundary conditions and loads

clear all                   % Clear all variables
clc;                        % Clear command window
clf;                        % Clear figure    
maxiter=100;                % Setting max # iterations
toler=1e-6;                 % For convergence criteria
n=80;                       % number of Finite Elements to be taken
len=1;                      % Length of the beam
h=len/n;                    % Element size
vol=0.01;                   % Volume limit
runi=sqrt(vol/len/pi);      % Initial radius that is feasible
Y=1;                        % Young's Modulus in N/cm^2
%  typ  boundary condition
%   1-> fixed-free
%   2-> pinned-pinned
typ=1;
r0=runi*ones(1,n);          % initial guess for width of elements
count=0;                    % for counting the no. of iterations
flag=1;                     % for checking if the convergence criteria are met
eta=0.01;                   % parameter for exponent and multiplying factor
rmin=0.001;
rmax=0.06;
while flag
    count=count+1;
    % Calling the fucnction for editing the data files
    modifydat4col(len,n,r0,Y,typ)
    % calling the script beam.m for the deflection by FEM analysis
    [w,eigval]=fem4buck(typ);    
    w1(1)=(-3*w(1)+4*w(2)-w(3))/2/h;
    w2(1)=(2*w(1)-5*w(2)+4*w(3)-w(4))/h^2;
    for i=2:n
        w1(i)=(w(i+1)-w(i-1))/2/h;
        w2(i)=(w(i+1)+w(i-1)-2*w(i))/h^2;
    end
    w1(n+1)=(3*w(n+1)-4*w(n)+w(n-1))/2/h;
    w2(n+1)=(2*w(n+1)-5*w(n)+4*w(n-1)-w(n-2))/h^2;
    w=(w(1:end-1)+w(2:end))/2;
    w1=(w1(1:end-1)+w1(2:end))/2;
    w2=(w2(1:end-1)+w2(2:end))/2;
    B=Y*r0.^2.*w2.^2/(2*h*sum(w1.^2));   % Calculating a term in the optimality criteria
    % Lagrange multiplier
    lambda=real(pi/vol*h*sum(r0.^2.*B.^(2*eta)))^(1/(2*eta));
    % new values of radius
    r1=r0.*(B/lambda).^eta;    
    for i=1:length(r1);
        if r1(i)<rmin             
            r1(i)=rmin;            
        end
        if r1(i)>rmax
            r1(i)=rmax;
        end
    end
    % checking the convergence of the width values
    maxerr=max(abs(r1-r0));
    if maxerr<toler || count>=maxiter
        flag=0;
    end
    r0=r1;
end
disp(['Number of Iterations required is:',num2str(count)]);
disp(['Corresponding Buckling Load is: ',num2str(eigval)]);
x=linspace(0,len,n+1);
xx=reshape([x(1:end-1);x(2:end)],1,2*n);
rr=reshape([r0;r0],1,2*n);
plot(rr,xx);
hold on;
plot(-rr,xx);
% axis equal;