%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [K,R,I]=assm_sys(edof,K,Ke,R,Re,I,Ie)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% assemble local element contributions to global stiffness and force
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% input:  edof  = [ elem X1 Y1 X2 Y2 ]   ... incidence matrix
%         Ke    = [ ndof x ndof ]        ... element tangent  Ke 
%         Re    = [ ndof x 1]            ... element residual Re
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% output: K     = [ 4 x 4 ]              ... global tangent   K
%         R     = [ fx_1 fy_1 fx_2 fy_2] ... global residual  R
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[nie,n]=size(edof);

I(edof(1),:)=Ie(:);

t=edof(:,2:n);
for i = 1:nie
  K(t(i,:),t(i,:)) = K(t(i,:),t(i,:))+Ke;
  R(t(i,:))         =R(t(i,:))       +Re;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
