%This file to compute the static response of a flexible structure.
% Define a vector of indices for each spring
% Indices represent node numbers
clear;
close all;
%Define springs by node ends
Spring = [1 2;3 1;1 4;4 2;2 5;4 5;2 5;3 4;3 6;6 7;6 4;4 7;7 5];
% Define nodes by (x,y) coordinates
a = 1/2;
b = 3^0.5/2;

nodes=[-a b;a b;-1 0;0 0;1 0;-a -b;a -b];

[u,v] = size(Spring);
[w,x] = size(nodes);


    

SpringV = zeros(u,v);
%Plot nodes, and define spring vectors
for i = 1:u,
    SpringV(i,:) = nodes(Spring(i,2),:)-nodes(Spring(i,1),:)
    plot([nodes(Spring(i,1),1);nodes(Spring(i,1),1)+SpringV(i,1)],[nodes(Spring(i,1),2);nodes(Spring(i,1),2)+SpringV(i,2)]);hold on;
end;

% Elongation coefficients
% for each spring
% Arranged as a matrix

SpringE = zeros(u,2*w);
for i = 1:u,
    %elongations along x axis
    SpringE(i,Spring(i,1))=-[1 0]*SpringV(i,:)';
    SpringE(i,Spring(i,2))=[1 0]*SpringV(i,:)';
    %elongations along y axis
    SpringE(i,w+Spring(i,1))=-[0 1]*SpringV(i,:)';
    SpringE(i,w+Spring(i,2))=[0 1]*SpringV(i,:)';
end;

% Now implement force balance equation for each node;
% Force terms due to Springs

FB = zeros(2*w,2*w);
for j = 1:w,
    for k=1:u,
        if Spring(k,1)==j,
            FB(j,:)=FB(j,:)+SpringE(k,:)*SpringV(k,1);
            FB(j+w,:)=FB(j+w,:)+SpringE(k,:)*SpringV(k,2);
        end;
        if Spring(k,2)==j,
            FB(j,:)=FB(j,:)-SpringE(k,:)*SpringV(k,1);
            FB(j+w,:)=FB(j+w,:)-SpringE(k,:)*SpringV(k,2);
        end;
    end;
end;
FB = [FB(:,3:7) FB(:,10:14)];
%Forces now;
FB(1,11)=1;
FB(2,12)=1;
FB(8,13)=1;
FB(9,14)=1;

%Now Force vector
F = 0.1*zeros(2*w,1);
F(5) = 2;
F(13) = -2;
F(14) = -1;
D=-FB\F*0.05;

%Total displacements
displ = [zeros(2,2);[D(1:5) D(6:10)]]

%Plot deformed structure
for i=1:u,
    plot([nodes(Spring(i,1),1)+displ(Spring(i,1),1);nodes(Spring(i,1),1)+displ(Spring(i,2),1)+SpringV(i,1)],[nodes(Spring(i,1),2)+displ(Spring(i,1),2);nodes(Spring(i,1),2)+SpringV(i,2)+displ(Spring(i,2),2)],'r');hold on;
end;