%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Beam Optimization MATLAB script for
% maximizing compliance 
% Written by: Anil Das P V
% March 2009.IISC
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clf					% Clear graphics window
clear all			% Clear all variables
clc					% Clear command window
hold off			% No hold on the graphics window

Nn=100;             %Number of Nodes
Ne=Nn-1;            %Number of Elements
Vs=5e-3;            %Volume bound
d=0.05;             %Depth of beam
Ltot=2;             %Total length
Le=Ltot/(Nn-1);     %Element length(Uniform Size)
Elas=200e9;         %Youngs Modulus (Uniform)
qn=-5e4;            %Load/len for UDL

bmax=10e-2;         %Upper bound on breadth
bmin=1.0e-2;          %Lower bound on breadth

toler=1e-3;         %Tolerance setting for convergence

eta=0.1;            %Rate controlling variable
cc=Elas*d^2/12;     

bini=Vs/(Ltot*d);   %Initial breadth
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% PRE-PROCESSING
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Computing x and y coordinates of the Nodes
ny=zeros(1,Nn);
for i=1:Nn  
    nx(i)=(i-1)*Le;
end
%Creating the Length,Youngs modulus,breadth array's
E=Elas*ones(1,Nn);
L=Le*ones(1,Ne);
b1=bini*ones(1,Ne);
%Defining connectivity
for i=1:Ne
    ncon(i,1)=i;
    ncon(i,2)=i+1;
end
%Displacement BC for Fixed-Fixed beam
for j=1:4
    dispVal(j) = 0;
end
dispID = [1,2,3,3*Nn-1];
%Creating of nodal force vector for UDL
F = zeros(3*Nn,1);
for i=2:(Nn-1)
    F((i-1)*3+2)=qn*Le;
end
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ITERATION FOR OPTIMIZING
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
flag1=true;         %Loop stop control variable
%OUTER ITERATION
while flag1
    flag1=false;
    A=b1*d;         %Area of each element
    Inertia=b1*d^3/12;  %Inertia of each element
    %Using the feam code
    [U,R,Fint,Fint_local,Kglobal,SE] = fembeam(A,L,E,nx,ny,ncon,Ne,Nn,F, ...
        dispID,dispVal,Inertia);
    %Computing w'' using central finite difference
    for i=1:Ne
        w1(i)=(U((i+1)*3)-U(i*3))/Le;
    end
    %Estimation of Lambda
    bsum=0;
    for i=1:Ne
        bsum=bsum+b1(i)*(cc*w1(i)^2)^eta;
    end
    Lam=(Vs/(d*Le*bsum))^(-1/eta);
    %Computation of modified breadth of element
    bs=b1.*(cc*w1.^2/Lam).^eta;
    flag=false;         %Loop stop control variable for Inner loop
    %Checking if Computed breadth is within Bounds  
    for i=1:Ne
        if(bs(i)>bmax)
            bs(i)=bmax;
            flag=true;
        elseif(bs(i)<bmin)
            bs(i)=bmin;
            flag=true;
        end
    end
    %INNER ITERATION
    while(flag)
        A=bs*d;
        Inertia=bs*d^3/12;
        [U,R,Fint,Fint_local,Kglobal,SE] = fembeam(A,L,E,nx,ny,ncon,Ne,Nn,F, ...
            dispID,dispVal,Inertia);
        for i=1:Ne
            ws(i)=(U((i+1)*3)-U(i*3))/Le;
        end
        %Estimation of Lambda
        countmax=0;
        countmin=0;
        bsum=0;
        for i=1:Ne
            if(bs(i)==bmax)
                countmax=countmax+1;
                continue;
            elseif(bs(i)==bmin)
                countmin=countmin+1;
                continue;
            end
            bsum=bsum+bs(i)*(cc*ws(i)^2)^eta;
        end
        Lam=((Vs-(countmax*bmax+countmin*bmin)*d*Le)/(d*Le*bsum))^(-1/eta);
        %Computation of modified breadth of element
        for i=1:Ne
            if (bs(i)~=bmax)&&(bs(i)~=bmin)
                bs(i)=bs(i)*(cc*w1(i)^2/Lam)^eta;
            end
        end
        flag=false;
        for i=1:Ne
            if(bs(i)>bmax)
                bs(i)=bmax;
                flag=true;
            elseif(bs(i)<bmin)
                bs(i)=bmin;
                flag=true;
            end
        end
    end
    if(max(abs(bs-b1))>toler)   %Checking for convergence
        flag1=true;
        b1=bs;
    end
end
disp(SE);
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% POST-PROCESSING
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Plotting
xte=linspace(0,Ltot,Nn);
j=2;
x(1)=xte(1);
for i=2:Ne;
    x(j)=xte(i);
    bte(j-1)=bs(i-1);
    x(j+1)=xte(i);
    bte(j)=bs(i-1);
    j=j+2;
end
x(j)=xte(Nn);
bte(j-1)=bs(Ne);
bte(j)=bs(Ne);
plot(x,bte,'k',x,-bte,'k');
hold on;
plot([0 Ltot],[0 0],'k-.');
axis([0 Ltot -bmax-0.5*bmax bmax+0.5*bmax]);
xlabel('x (m)');
ylabel('b (m)');
title('Beam Area Profile');
hold off;