%
% skin.m    Steve Cox, Linear Algebra in situ
%
% Given a data array of node numbers and fiber angles generate
% the incidence matrix, A, for a small patch of skin. Assuming unit
% fiber stiffnesses, K, build A'KA. Then apply uniform traction, f, at the
% outer nodes, solve A^TKA x = f for x and plot the deformed and 
% undeformed nodes.
%

data = [			% one row of data for each fiber, the
1 	4	pi/2		% first two columns are end-node numbers
1	5	pi/4		% while the third is the fiber angle
1	2	0
2	4	3*pi/4
2	5	pi/2
2	6	pi/4
2	3	0
3	5	3*pi/4
3	6	pi/2
4	7	pi/2
4	8	pi/4
4	5	0
5	7	3*pi/4
5	8	pi/2
5	9	pi/4
5	6	0
6	8	3*pi/4
6	9	pi/2
7	8	0
8	9	0];

A = zeros(20,18);		% initialize the incidence matrix

for j=1:20,

    m = data(j,1);
    n = data(j,2);
    theta = data(j,3);

    A(j,2*n-1) = cos(theta);	% direct implementation of general formula
    A(j,2*m-1) = -cos(theta);
    A(j,2*n) = sin(theta);
    A(j,2*m) = -sin(theta);

end

figure(1)
spy(A)
 
K = eye(20);    % unit fiber stiffnesses
 
S = A'*K*A;     % the composite stiffness matrix

s = 1/sqrt(2);

f = randn(18,1); % [-s -s 0 -1 s -s -1 0 0 0 1 0 -s s 0 1 s s]';   % uniform traction

pS = pinv(S);   % the pseudoinverse if S

x = pS*f;       % the solution to A'KAx=f (with rigid body motions)

% draw the nodes of the original and deformed specimens

figure(2)

for h=1:3,    % plot the reference circle and its square translate

    plot(0,(h-1)*6,'ko','markersize',12)	% original
    hold on
    plot(x(1+(h-1)*6),(h-1)*6 + x(2+(h-1)*6),'rs','markersize',12') % deformed
    plot(6,(h-1)*6,'ko','markersize',12')
    plot(6+x(3+(h-1)*6),(h-1)*6 + x(4+(h-1)*6),'rs','markersize',12')
    plot(12,(h-1)*6,'ko','markersize',12')
    plot(12+x(5+(h-1)*6),(h-1)*6 + x(6+(h-1)*6),'rs','markersize',12')

end

axis equal
axis off
hold off

