% Stiffest bar for given volume and axial load
% Solved using the optimality criteria method.
% Code written by Gaurav Nair in 2011, IISc-Bangalore
% Commented and modified thereafter for ME 256 course in IISc-Bangalore
%
close all; clear all; clc; hold off; % Housekeeping commands
% Data
length = 10; n = 100; vstar = 15;  %data (length, no of elements,volume)
Ae = 6*ones(n,1); % Initial approx of area
Le = length/n*ones(n,1);
Ee = 210*ones(n,1); %elastic modulus
nx = zeros(n,1);
Amax = 3; %upper bound on area
Amin = 0.5; %lower bound on area
tol = 1e-3; %tolerance to stop iteration
for i=1:n
    nx(i+1) = nx(i)+Le(i); %x co-ordinates of the nodes along the axis of the bar
end
ny = zeros(n+1,1); %y co-ordinates of the nodes along the axis of the bar
F = ones(n+1,1); %Uniformly distributed force
%F(n/2,1)=50;%F(n+1,1)=100;
ncon = zeros(n,2);
for i=1:n
    ncon(i,1)=i;
    ncon(i,2)=i+1;
end
dispID = [1; n+1]%;n+1]; %for Fixed free condition , for fixed fixed dispID=[1;n+1] and for free fixed dispID=[n+1]
dispVal = [0; 0]; %;0]; %for Fixed free condition and free fixed condition, for fixed fixed condition dispVal=[0;0]
NELEM=n; %Number of elements
NNODE=n+1; %Number of nodes
beta = 0.1;
its=50; % Max number of iterations
history_A=zeros(n,its+1); % Natrix to store values of Area after every iteration
history_A(:,1)=Ae(:,1);
for i=1:its
    lambda1=0;
    [u,Rdisp,P,Ksing,SE] = fembar(Ae, Le, Ee, nx, ny, ncon,NELEM,NNODE,F,dispID,dispVal);
    for a=1:n
        udash(a)=(u(a+1)-u(a))/Le(a); 
    end
    for j=1:n
        lambda1=lambda1+Ae(j)*Le(j)*(Ee(j)*udash(j)^2)^beta;
    end
    lambda=(lambda1/vstar)^(1/beta);
    for j=1:n
        Ae(j)=Ae(j)*(Ee(j)*udash(j)^2/lambda)^beta;
    end
    flag=false;
    for j=1:n
        if (Ae(j)>Amax)
            Ae(j)=Amax;
            flag=true;
        end
        if (Ae(j)<Amin)
            Ae(j)=Amin;
            flag=true;
        end
    end
    while flag
        [u,Rdisp,P,Ksing,SE] = fembar(Ae, Le, Ee, nx, ny, ncon,NELEM,NNODE,F,dispID,dispVal);
        for a=1:n
            udash(a)=(u(a+1)-u(a))/Le(a);
        end
        upcount=0;lowcount=0;lambda1=0;
        for j=1:n
            if (Ae(j)==Amax)
                upcount=upcount+1;
                continue;
            end
            if (Ae(j)==Amin)
                lowcount=lowcount+1;
                continue;
            end
            lambda1=lambda1+Ae(j)*Le(j)*(Ee(j)*udash(j)^2)^beta;
        end
        vstar1=vstar-(upcount*Amax+lowcount*Amin)*length/n;
        lambda=(lambda1/vstar1)^(1/beta);
        for j=1:n
            if Ae(j)~=Amax && Ae(j)~=Amin
                Ae(j)=Ae(j)*((Ee(j)*udash(j)^2)/lambda)^beta;
            end
        end
        flag=false;
        for j=1:n
            if (Ae(j)>Amax)
                Ae(j)=Amax;
                flag=true;
            end
            if (Ae(j)<Amin)
                Ae(j)=Amin;
                flag=true;
            end
        end
    end
    history_A(:,i+1)=Ae(:,1);
    % To check convergence
    tolcheck=max(abs(history_A(:,i+1)-history_A(:,i)));
    k=i;
    if tolcheck<=tol
        break;
    end
end
history_A(n+1,:)=history_A(n,:);
iter = k;
% Plotting changing area during iteration
% for j = 1:iter,
%     clf
% for i=1:n
%     X = length/n * [i-1 i i i-1 i-1];
%     Y = history_A(i,j) * [-0.5 -0.5 0.5 0.5 -0.5];
%     fill(X,Y,'c');
%     hold on
% end
% pause
% end
% Animating the iterative progress during optimization
for j = 1:iter,
    clf
    axis([0 10 -Amax/2 Amax/2]);
    hold on
    for i=1:n
        X = length/n * [i-1 i i i-1 i-1];
        Y = history_A(i,j) * [-0.5 -0.5 0.5 0.5 -0.5];
        fill(X,Y,'c');
    end
    pause(1)
end
xlabel('x');ylabel('Area of cross-section, A(x)');
grid on
% Plotting the result of the final iteration
for i=1:n
    X = length/n * [i-1 i i i-1 i-1];
    Y = history_A(i,iter) * [-0.5 -0.5 0.5 0.5 -0.5];
    fill(X,Y,'c');
    hold on
end
xlabel('x');ylabel('Area of cross-section, A(x)');
grid on

