% Coded by Safvan P (11711,Phd Student)for ME256 Homework #7
% QUESTION 1 :
% 18-3-2015
% ==========================================================
%Modify the Matlab code given for bar-optimization to solvin-
%g beam-optimization problems.
% ==========================================================
close all; clc;clear;
%======================================
% Constants
%======================================
nElem = 300;                                 %no of elements
length=20;                               %length of the beam
vStar = 50;                       %maximum volume constraint
Amax =6;                            %maximum area constraint
Amin = 1;                           %minimum area constraint
tol = 1e-3;        %tolerance for optimality criteria method

depth=1;                   %depth of the constant depth beam
alpha=depth^2/12;                        %Inertia=alpha*area
width=ones(nElem,1);

%======================================
% Inputs to fembeam()
%======================================
areaElem=width*depth;           %area of individual elements

lengthElem = length/nElem*ones(nElem,1);%length of individua
%l elements

youngsElem = 210*ones(nElem,1);              %youngs modulus

xElem = zeros(nElem,1);%preallocation of Element X-cordinates
for i=1:nElem
    xElem(i+1) = xElem(i)+lengthElem(i);
end

yElem= zeros(nElem+1,1);%preallocation of Element Y-cordinates

con = zeros(nElem,2);          %connectivty matrix for fembeam
for i=1:nElem
    con(i,1)=i;
    con(i,2)=i+1;
end

nNodes=nElem+1;                               %number of nodes

F = zeros(3*nNodes,1);
for q=1:nNodes
    F(3*q-1)=1;
end
%force applied

dispID = [1,2,3,3*nNodes-2,3*nNodes];%node ID's where displacement is constrained
boundaryText='fixed-free';

inertia = alpha*areaElem;                %intertia of elements


eta=.2;                    %eta for optimality criteria method
its=50;                                      %no of iterations

historyArea=zeros(nElem,its+1); %to store values of Area after
%every iteration
historyArea(:,1)=areaElem(:,1);

ux=zeros(1,nNodes);              %x diaplcement from fembeam()
uy=zeros(1,nNodes);              %y diaplcement from fembeam()
uDash=zeros(1,nNodes);       %theta diaplcement from fembeam()
uDashDash=zeros(1,nElem);                 %derivtaive of uDash
%=========================================
% Iteration for optimality criteria method
%=========================================
for i=1:its
    lambda1=0;
    inertia = alpha*areaElem;          %updated inertia matrix
    [u,~] = fembeam(areaElem, lengthElem, youngsElem, xElem, yElem, con,...
        nElem,nNodes,F,dispID,inertia);
    for q=1:nNodes
        ux(q)=u(3*q-2);
        uy(q)=u(3*q-1);
        uDash(q)=u(3*q);
    end
    for a=1:nElem
        uDashDash(a)=(uDash(a+1)-uDash(a))/lengthElem(a);
    end
    for j=1:nElem
        lambda1=lambda1+areaElem(j)*lengthElem(j)*(alpha*youngsElem(j)*...
            uDashDash(j)^2)^eta;
    end
    lambda=(lambda1/vStar)^(1/eta);       %lagrange multiplier
    for j=1:nElem
        areaElem(j)=areaElem(j)*(alpha*youngsElem(j)*uDashDash(j)^2/...
            lambda)^eta;                         %updated area
    end
    flag=false;       %flag to check area constraint violation
    %======================================
    % Inner loop to check area constraints
    %======================================
    for j=1:nElem
        if (areaElem(j)>Amax)
            areaElem(j)=Amax;
            flag=true;
        end
        if (areaElem(j)<Amin)
            areaElem(j)=Amin;
            flag=true;
        end
    end
    
    while flag
        inertia = alpha*areaElem;      %updated inertia matrix
        [u,~] = fembeam(areaElem, lengthElem, youngsElem, xElem, yElem,...
            con,nElem,nNodes,F,dispID,inertia);
        for q=1:nNodes
            ux(q)=u(3*q-2);
            uy(q)=u(3*q-1);
            uDash(q)=u(3*q);
        end
        for a=1:nElem
            uDashDash(a)=(uDash(a+1)-uDash(a))/lengthElem(a);
        end
        upcount=0;lowcount=0;lambda1=0;
        for j=1:nElem
            if (areaElem(j)==Amax)
                upcount=upcount+1;
                continue;
            end
            if (areaElem(j)==Amin)
                lowcount=lowcount+1;
                continue;
            end
            
            lambda1=lambda1+areaElem(j)*lengthElem(j)*(alpha*...
                youngsElem(j)*uDashDash(j)^2)^eta;
        end
        vStar1=vStar-(upcount*Amax+lowcount*Amin)*length/nElem;
        %updated volume constraint for the
        %inner loop
        lambda=(lambda1/vStar1)^(1/eta);
        for j=1:nElem
            if areaElem(j)~=Amax && areaElem(j)~=Amin
                areaElem(j)=areaElem(j)*(alpha*youngsElem(j)*...
                    uDashDash(j)^2/lambda)^eta;
            end
        end
        flag=false;
        for j=1:nElem
            if (areaElem(j)>Amax)
                areaElem(j)=Amax;
                flag=true;
            end
            if (areaElem(j)<Amin)
                areaElem(j)=Amin;
                flag=true;
            end
        end
    end
    historyArea(:,i+1)=areaElem(:,1);
    %area after every iteration stored
    %to history area
    %========================
    % check for convergence
    %========================
    tolcheck=max(abs(historyArea(:,i+1)-historyArea(:,i)));
%     disp(tolcheck);
    k=i;
    if tolcheck<=tol
        break;
    end
%     disp(tolcheck);
end
j=k;

%=========================================
% Plot 3-D beam
%=========================================

for f=1:nElem
    width(f)=historyArea(f,j)/depth;
    x1=-width(f)/2;
    x2=+width(f)/2;
    x3=+width(f)/2;
    x4=-width(f)/2;
    y1=-depth/2;
    y2=-depth/2;
    y3=depth/2;
    y4=depth/2;
    fill3(xElem(f)*ones(1,4),[x1,x2,x3,x4],[y1,y2,y3,y4],'r')
    hold on
end
axis equal
title(['3D view, Amax=',num2str(Amax),' Amin=',num2str(Amin),...
    ' for boundary condition',boundaryText,',uniform load']);
%==========================================
% 2D Plot of variation in area crosssection
%==========================================
historyArea(:,:)=historyArea(:,:)/depth;
historyArea(:,:)=historyArea(:,:)/2;
figure
for f=1:nElem
    plot([xElem(f);xElem(f+1)],[historyArea(f,j);historyArea(f,j)],...
        'Linewidth',3);
    hold on
    plot([xElem(f);xElem(f+1)],[-historyArea(f,j);-historyArea(f,j)],...
        'Linewidth',3);
    hold on
end
axis equal
grid on
title(['Amax=',num2str(Amax),' Amin=',num2str(Amin),...
    ' for boundary condition',boundaryText,',uniform load']);



