% 
% optimal design of a cantilever, c cells wide and d cells high
% 
% the cantilever is fixed at its left face and loaded at the center
% of its right face, assuming that d is even.
%
%  usage   mich3(c,d)    d must be even
%
%  e.g.,  mich3(4,6)
%

function mich3(c,d)

n = c*(d+1);		    % nodes
m = 2*c*d + c*(d+1) + c*d;  % fibers

A = zeros(m,2*n);
L = zeros(m,1);
h = sqrt(1/d^2 + 4/c^2);   % length of each diagonal bar
t = asin(1/d/h);	   % angle of diagonal bar
ct = cos(t);
st = sin(t);

for j=1:d,		% build the adjacency matrix and bar length vector
                        % one vertical layer at a time
    rc = (j-1)*4*c;
    cc = (j-1)*2*c+1;   % index of horizontal motion of SE node of 1st layer
    uc = cc + 2*c;

    % cell against wall is not like the rest

    A(1+rc,cc) = 1;	% horizontal
    L(1+rc) = 2/c;

    A(2+rc,[uc uc+1]) = [ct st];   % NE-SW diagonal
    L(2+rc) = h;

    A(3+rc,[cc cc+1]) = [ct -st];  % NW-SE diagonal
    L(3+rc) = h;

    A(4+rc,[cc+1 uc+1]) = [-1 1];  % vertical
    L(4+rc) = 1/d;

    for i=1:c-1		% now do the rest

        cc = cc + 2;
        uc = uc + 2;
        rc2 = rc + 4*i;

        A(1+rc2,[cc cc-2])  = [1 -1];   % horizontal
        L(1+rc2) = 2/c;

        A(2+rc2,[cc-2 cc-1 uc uc+1])  = [-ct -st ct st];  % NE-SW diagonal
        L(2+rc2) = h;

        A(3+rc2,[cc cc+1 uc-2 uc-1])  = [ct -st -ct st];  % NW-SE diagonal
        L(3+rc2) = h;

        A(4+rc2,[cc+1 uc+1])  = [-1 1];   % vertical
        L(4+rc2) = 1/d;

   end

end

% finally, lay down the top most horizontal fibers

rc = m-c+1;
cc = cc+2;
A(rc,cc) = 1;
L(rc) = 2/c;

for i=1:c-1

    cc = cc+2;
    A(rc+i,[cc-2 cc]) = [-1 1];
    L(rc+i) = 2/c;

end

indf = c+c*(d+1);		% center right index where load is applied
f = zeros(2*n,1);
f(indf) = -0.5;

V = 20;				% total volume
a0 = V./L/m;			% initial design

LB = .001*ones(m,1);            % lower bounds on bar x-sectional areas
UB = V*ones(m,1); 		% upper bounds on bar x-sectional areas

% set options for fmincon and get best design and its compliance

options = optimset('display','iter','gradobj','on');
[besta,bestc] = fmincon(@(a) cmit(a,A,L,f),a0,[],[],L',V,LB,UB,[],options);

% now display its findings

fibdraw(c,d,besta)
title(['Compliance = ' num2str(bestc)],'fontsize',16)

%
%  [comp,grad] = cmit(a,A,L,f)
%
% evaluate the compliance (comp) and its gradient wrt a (grad)
% of the net with adjacency matrix A, fiber lengths L, fiber
% x-sec areas a, and load f.
%

function [comp,grad] = cmit(a,A,L,f)
S = A'*diag(a./L)*A;  % the stiffness matrix
x = S\f;	      % the displacement
comp = x'*f;          % the work done by the load
e = A*x;	      % the fiber elongations
grad = -e.*e./L;      % the gradient


%
%  fibdraw(c,d,a)
%
%  draw the c-by-d cantilever with fiber x-sec areas a
%
%  the logic behind the indices is exactly that used above
%  in the construction of the adjacency matrix
%

function fibdraw(c,d,a)

n = c*(d+1);	% nodes
m = 2*c*d + c*(d+1) + c*d;  % fibers

for j=1:d,

    rc = (j-1)*4*c;
    cc = (j-1)*2*c+1;   % index of horizontal motion of SE node of 1st layer
    uc = cc + 2*c;

    xc(1+rc,:) = [0 2/c];
    yc(1+rc,:) = [(j-1)/d (j-1)/d];

    xc(2+rc,:) = [0 2/c];
    yc(2+rc,:) = [(j-1)/d j/d];

    xc(3+rc,:) = [0 2/c];
    yc(3+rc,:) = [j/d (j-1)/d];

    xc(4+rc,:) = [2/c 2/c];
    yc(4+rc,:) = [(j-1)/d j/d];

    for i=1:c-1

        cc = cc + 2;
        uc = uc + 2;
        rc2 = rc + 4*i;

        xc(1+rc2,:) = [i*2/c (i+1)*2/c];
        yc(1+rc2,:) = [(j-1)/d (j-1)/d];

        xc(2+rc2,:) = [i*2/c (i+1)*2/c];
        yc(2+rc2,:) = [(j-1)/d j/d];

        xc(3+rc2,:) = [i*2/c (i+1)*2/c];
        yc(3+rc2,:) = [j/d (j-1)/d];

        xc(4+rc2,:) = [(i+1)*2/c (i+1)*2/c];
        yc(4+rc2,:) = [(j-1)/d j/d];

   end

end

rc = m-c+1;
cc = cc+2;
xc(rc,:) = [0 2/c];
yc(rc,:) = [1 1];

for i=1:c-1

    cc = cc+2;
    xc(rc+i,:) = [i*2/c (i+1)*2/c];
    yc(rc+i,:) = [1 1];

end

for j=1:m
    plot(xc(j,:),yc(j,:),'linewidth',ceil(5*a(j)))
    hold on
end

fill([-.25 0 0 -.25],[-.1 -.1 1.1 1.1],'w')  % draw the wall
axis equal
axis([-.3 2.2 -0.5 1.1]) 
axis off
hold off
return