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);