%
%  cantilever.m
%
%  solve and diplay a tip loaded cantilever
%
%  usage:  cantilever(Ea,F,nof,mov)
%
%  where:  Ea = vector of fiber modulus-x-area  prodcuts
%          F = strength of the downward force at the upper right tip
%          nof = number of frames
%          mov = 1 for movie or 0 for Duchamp
%
%  example:  cantilever(ones(16,1),.01,10,0);
%

function cantilever(Ea,F,nof,mov)

close all
m = length(Ea);
nos = m/4; 	% number of stories

s = 1/sqrt(2);
A = zeros(m);
A(1,[1 2]) = [s s];   % build 1st story of A
A(2,2) = 1;
A(3,[3 4]) = [s s];
A(4,[1 3]) = [-1 1];

for i=1:nos-1,		% build remaining stories
    hl = 4*i - 3;       % index hor-dof of bottom left node
    A(5+(i-1)*4,[hl hl+1 hl+4 hl+5]) = [-s -s s s];
    A(6+(i-1)*4,[hl+3 hl+5]) = [-1 1];
    A(7+(i-1)*4,[hl+2 hl+3 hl+6 hl+7]) = [-s -s s s];
    A(8+(i-1)*4,[hl+4 hl+6]) = [-1 1];
end

figure(1)		% visualize the adjacencies
spy(A)
xlabel('degree of freedom','fontsize',16)
ylabel('fiber','fontsize',16)
title(['Nonzero adjacencies, nos = ' num2str(nos)],'fontsize',16)

for i=0:nos-1	% construct the coordinates of the undeformed net

    xc(1+i*4:4+i*4,:) = [i-1 i; i i; i i+1; i i+1];
    yc(1+i*4:4+i*4,:) = [i i+1; i i+1; i i+1; i+1 i+1];

end

figure(2)

line(xc',yc','linewidth',2)			% draw the undeformed net
hold on
fill([-1.2 nos+2 nos+2 -1.2],[0 0 -.2 -.2],'w')    % draw the foundation
axis off
axis equal
axis([-1.5 1.75*nos -.5 nos+1])		% set a fixed frame size

if mov
   M(1) = getframe;			% get the inital frame
   hold off
   clf
end

s = 1/sqrt(2);
L = ones(size(Ea));
L(1:2:end) = 1/s;

k = Ea./L;			% the fiber stiffnesses
K = diag(k);

S = A'*K*A;			% the mother stiffness matrix

f = zeros(4*nos,1);

for j = 2:nof,
 
    f(end) = -F*(j-1)/(nof-1);  % scale f to achieve max at final fram

    x = S\f;		
    
    xh = [0; 0; x(1:2:end)];
    yh = [0; 0; x(2:2:end)];
   
    for i=0:nos-1   % get the deformed coordinates

        dx(1+i*4:4+i*4,:) = xc(1+i*4:4+i*4,:) + [xh(1+2*i) xh(3+2*i) 
                                                 xh(2+2*i) xh(3+2*i)
                                                 xh(2+2*i) xh(4+2*i)
                                                 xh(3+2*i) xh(4+2*i)];

        dy(1+4*i:4+4*i,:) = yc(1+i*4:4+i*4,:) + [yh(1+2*i) yh(3+2*i);
                                                 yh(2+2*i) yh(3+2*i)
                                                 yh(2+2*i) yh(4+2*i)
                                                 yh(3+2*i) yh(4+2*i)];


    end
    
    line(dx',dy','linewidth',2)		% draw the deformed net
    hold on
    fill([-1.2 nos+2 nos+2 -1.2],[0 0 -.2 -.2],'w')    % draw the foundation
    axis equal
    axis([-1.5 1.75*nos -.5 nos+1])	
    axis off
    if mov
       M(j) = getframe;		% get the jth frame
       hold off
       clf
    end


end

if mov
   movie(M,-4,10)	% play the movie 4 times, both back and forth, at
                        % 10 frames/second

   movie2avi(M,'flw','quality',100)  % save a high quality portable copy
                                       % for your portfolio
end