createMatrices.m

function [Mass, C] = createMatrices(truss)
% berechnet die Massen- und Steifigkeitsmatrix des gegebenen Fachwerks
% Indizes i, j von 1 .. M
% Index k von 1 .. N
% Indizes a, b  Werte 1, 2  (= x, y)
 
M = size(truss.A,1);
N = truss.N;

% Verbindungsvektoren e(i,j,a)
e0 = zeros(M,M,2);
for I=1:M
  for J = 1:(I-1)
    % oben rechts
    e0(I,J,:) = truss.x0(:,I) - truss.x0(:,J);  
    e0(I,J,:) = e0(I,J,:)/norm(squeeze(e0(I,J,:)));
    % unten links durch Antisymmetrie
    e0(J,I,:) = -e0(I,J,:);
  end
end

%
% Steifigkeitsmatrix C
%
% Doppelindex j a in Einfachindex m
%   m = 2*(k-1) + a
% also: 1x, 1y, 2x, 2y, 3x, 3y, ...
C = zeros(2*N,2*N);

% Diagonalelemente
for K = 1:N
  for J = 1:M
    C(2*K-1,2*K-1) = C(2*K-1,2*K-1) + truss.A(K,J)*e0(K,J,1)^2;  % C(kx,kx)
    C(2*K-1,2*K  ) = C(2*K-1,2*K  ) + truss.A(K,J)*e0(K,J,1)*e0(K,J,2); % C(kx,ky)
    C(2*K  ,2*K  ) = C(2*K  ,2*K  ) + truss.A(K,J)*e0(K,J,2)^2;  % C(ky,ky)
  end
  C(2*K,2*K-1) = C(2*K-1,2*K);  % C(ky,kx) = C(kx,ky)
end

% Nicht-Diagonalelemente
for K = 1:N
  for J = 1:(K-1)
     % oben rechts
     C(2*K-1,2*J-1) = -truss.A(K,J)*e0(K,J,1)^2;  % C(kx,jx)
     C(2*K-1,2*J  ) = -truss.A(K,J)*e0(K,J,1)*e0(K,J,2); % C(kx,jy)
     C(2*K  ,2*J-1) = C(2*K-1,2*J  );   % C(ky,jx)
     C(2*K  ,2*J  ) = -truss.A(K,J)*e0(K,J,2)^2;  % C(ky,jy)
     % unten links durch Symmetrie
     C(2*J-1,2*K-1) = C(2*K-1,2*J-1);
     C(2*J  ,2*K-1) = C(2*K-1,2*J  );
     C(2*J-1,2*K  ) = C(2*K  ,2*J-1);
     C(2*J  ,2*K  ) = C(2*K  ,2*J  );
  end
end
C = truss.c * C;

Mass = truss.m*eye(2*N,2*N);