function out=element_truss(task,varargin)
% stiff=element_truss('stiff',coef,nodes)
% input
%   coef    the number A*E
%   nodes   columns are node coordinates, any dimension
% output
%   stiff     the element stiffness matrix
%   

switch task
    case 'stiff'
        % compute local stiffness matrix
        out=element_truss_stiff(varargin{:});     
    otherwise
        error('bad function')
end

function stiff=element_truss_stiff(coef,nodes)

if size(nodes,2) ~= 2,
    error('truss has 2 nodes')
end
dim=size(nodes,1);  % the dimension
d=nodes(:,2)-nodes(:,1); % element vector
% Hookes law:
% deltaL = L*F/A*E
% deltaL is elongation in axial direction
% L is length of the element
% F is force in the axial direction
% A is the crossection
% E is the Young's modulus
% here: 
% A*E=coef
% L = norm(d) where d=node(:,2)-node(:,1) is element vector
% deltaU = U(:,2)-U(:,1) is relative displacement vector
% deltaL = (projection of deltaU on direction d) = (d*d')/(d'*d)*deltaU
% force F = deltaL*A*E/L acts on both ends in direction of d
%         -F <- node(1) xxxxxxxxxxxxxxxxxxxxxxx node(2) -> +F
% denote f = [F    
%            -F] 
% The local stiffness matrix returns
% f = stiff*u
% assume the dofs are ordered all node1, all node2
% u = [U(:,1)   % displacements at node 1
%      U(:,2)]; % displacements at node 2
dim=size(nodes,1);
d = nodes(:,2)-nodes(:,1);
delta=[eye(dim),-eye(dim)];  % deltaU = delta*u, f=delta'*F
% L = norm(d);
% stiff=(coef/L)*delta'*d*d'/(L*L)*delta
% a bit more efficient
L1=1/norm(d);
dd=(d*L1)'*delta;
stiff=(coef*L1)*dd'*dd;