% Ratkaistaan ristikkorakenteen siirtymät elementtimenetelmällä
% e viittaa aina elementtikohtaisiin matriiseihin

clear;

% Parametrien maarittelyja
d=0.01;		% Sauvan halkaisija
A=pi*(d/2)^2;	% Poikkipinnan ala
E=200e9;	% Kimmokerroin
EA=E*A;		% Aksiaalijäykkyys

% Solmupisteet
Ns=9;
X=[0 0 1 1 2 2 3 3 4	% x-koordinaatit
   0 1 0 1 0 1 0 1 0];	% y-koordinaatit

% Reunaehdot ( 1=kiinnitetty, 0=vapaa )
D=[1 1 0 0 0 0 0 0 0	% x-siirtymä
   1 1 0 0 0 0 0 0 0];	% y-siirtymä

% Ulkoiset pistekuormat
L=[0 0 0 0 0 0 0 0 0		% x-voima
   0 0 0 0 0 0 0 0 -1e5];	% y-voima

% Elementit
Ne=15;
E=[1 1 2 2 3 4 4 5 6 6 7 8 3 5 7	% Solmusta Xk
   2 3 4 3 5 6 5 7 8 7 9 9 4 6 8];	% solmuun Xl

% Piirretään rakenne kuvaan 1
figure(1); clf; hold on;
for e=1:Ne;
	k=E(1,e); l=E(2,e);	% Elementin e solmujen numerot k ja l
	Xk=X(:,k); Xl=X(:,l);	% Solmupisteet Xk ja Xl
	plot([Xk(1) Xl(1)],[Xk(2) Xl(2)]);
end;
axis image; axis off;

disp('Paina namiskuukkelia');
pause

% Muodostetaan muunnosmatriisi T
Nv=0;
for i=1:Ns;
	if D(1,i) == 0; Nv=Nv+1; T(2*i-1,Nv)=1; end;
	if D(2,i) == 0; Nv=Nv+1; T(2*i  ,Nv)=1; end;
end;

% Muodostetaan vapaan rakenteen jäykkyysmatriisi elementeittäin
K=zeros(2*Ns,2*Ns);
for e=1:Ne;
	k=E(1,e); l=E(2,e);	% Elementin solmupisteiden numerot
	Xk=X(:,k); Xl=X(:,l);	% Elementin solmupisteiden koordinaatit
	he=norm(Xk-Xl);		% Elementin pituus
	de=(Xk-Xl)/he;		% Elementin yksikkösuuntavektori
	De=[de' [0 0]		% Suuntamatriisi
	    [0 0] de'];
	ke=(EA/he)*[1 -1;-1 1];	% Aksiaalijäykkyysmatriisi
	Ke=De'*ke*De;		% Elementin jäykkyysmatriisi
	Ae=zeros(4,2*Ns);	% Asennusmatriisi
	Ae(1,2*k-1)=1; Ae(2,2*k)=1; Ae(3,2*l-1)=1; Ae(4,2*l)=1;
	K=K+Ae'*Ke*Ae;		% Summataan globaaliin jäykkyysmatriisiin
end;

% Muodostetaan vapaan rakenteen voimavektori (vain solmupistevoimia)
F=zeros(2*Ns,1);
for i=1:Ns; F(2*i-1)=L(1,i); F(2*i)=+L(2,i); end;

% Ratkaistaan vapausasteet (red=redusoitu matriisi)
Kred=T'*K*T; Fred=T'*F; Ured=inv(Kred)*Fred; U=T*Ured;

% Tulostetaan siirtymät
for e=1:Ne;  
	k=E(1,e); l=E(2,e);	% Elementin solmupisteiden numerot
	Xk=X(:,k); Xl=X(:,l);	% Elementin solmupisteiden koordinaatit
	Ae=zeros(4,2*Ns);	% Asennusmatriisi
	Ae(1,2*k-1)=1; Ae(2,2*k)=1; Ae(3,2*l-1)=1; Ae(4,2*l)=1;
	Ue=Ae*U;
	plot([Xk(1)+Ue(1) Xl(1)+Ue(3)],[Xk(2)+Ue(2) Xl(2)+Ue(4)], '--');
end;

axis image; axis off;

