clear all           % memoria torlese
clc                 % Command Window torlese
format shortG       % Command window megjelenítési formátuma
pkg load signal     % signal csomag betoltese

%% PREPROCESSING %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% tulajdonsagok
A=706.86;           % keresztmetszet terulete
E=200000;           % rugalmassagi modulus
% csomoponti koordinatak (x y (z))
crd=[0 0;0 10000
     5000 0;5000 10000
     2*5000 0;2*5000 10000
     3*5000 0;3*5000 10000
     4*5000 0;4*5000 10000
     5*5000 0;5*5000 10000
     6*5000 0;6*5000 10000];
% kapcsolati matrix (kezdopont, vegpont)
con=[1 2;3 4;5 6;7 8;9 10;11 12;13 14   % fuggoleges rudak
     1 3;3 5;5 7;7 9;9 11;11 13         % lenti rudak
     2 4;4 6;6 8;8 10;10 12;12 14       % fenti rudak
     1 4;3 6;5 8;9 8;11 10;13 12];      % atlos rudak
% megfogasok
kbc=[1  1 
     1  2 
     13 2];
% dinamikai peremfeltetelek
dbc=[2 1 20000
     7 2 -20000
     9 2 -40000];
% megjelenitesi beallitasok
sc=200;     % nagyitas merteke

%% SOLVING %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
dim=size(crd,2);            % dimenzio (2D vagy 3D)
non=size(crd,1);            % csomopontok szama
noe=size(con,1);            % elemek (=rudak) szama
dof=non*dim;                % szerkezet szabadsagfokszama

K=zeros(dof,dof);           % globalis merevsegi matrix inicializalasa
F=zeros(dof,1);             % csomoponti tehervektor inicializalasa
F((dbc(:,1)-1).*dim+dbc(:,2))=dbc(:,3); % adott erok tehervektorba rakasa
q=zeros(dof,1);             % csomoponti elmozdulasvektor inicializalasa

Pset=(kbc(:,1)-1).*dim+kbc(:,2); % ismert elmozdulaskoordinatak sorszama
Nset=setdiff([1:dof]',Pset);     % szabad elmozdulaskoordinatak sorszama

N1=con(:,1);                % elemek kezdopontjainak sorszamai
N2=con(:,2);                % elemek vegpontjainak sorszamai
v=crd(N2,:)-crd(N1,:);      % elemek iranyvektorai
L=rssq(v,2);                % elemek hosszai
uv=v./L;                    % elemek iranyegysegvektorai

% merevsegi matrix eloallitasa:
for e=1:noe
    T=blkdiag(uv(e,:),uv(e,:));     % elmozdulás transzformációs mátrixa   
    K1=A*E/L(e);                    % elem merevsegi matrixanak elemei
    Ke=T'*[K1 -K1;-K1 K1]*T;        % elem transzformalt merevsegi matrixa
    d1=(N1(e)-1)*dim+1:N1(e)*dim;   % elem kezdopontjanak elmozdulas sorszamai
    d2=(N2(e)-1)*dim+1:N2(e)*dim;   % elem vegpontjanak elmozdulas sorszamai
    d=[d1,d2];                      % elem csomopontjainak elmozdulas sorszamai
    K(d,d)=K(d,d)+Ke;               % elhelyezes a globalis merevsegi matrixba
end
% kinematikai peremfeltetelek figyelembe vetele:
Knn=K(Nset,Nset);       % merevsegi matrix szabad resze
Fn=F(Nset);             % tehervektor szabad resze
% elmozdulasok kiszamitasa
q(Nset)=Knn\Fn;         % Kq=F egyenlet megoldasa

%% POSTPROCESSING %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% elmozdulasbol szarmaztatott mennyisegek kiszamitasa
Fall=round(K*q.*1e8)./1e8;             % csomoponti erok kiszamítasa
qmat=reshape(q,dim,non)';              % csomoponti elmozdulasmatrix
qabs=rssq(qmat,2)                      % elmozdulasok abszolut erteke
Fmat=reshape(Fall,dim,non)';           % csomoponti eromatrix
DL=sum(uv.*(qmat(N2,:)-qmat(N1,:)),2); % rudak hosszvaltozasa
DL=round(DL.*1e8)./1e8;                % hosszvaltozasok kerekitese
EPS=DL./L;                             % rudak fajlagos nyulasa
S=E.*EPS                               % rudakban fellepo normalfeszultseg
N=A.*S;                                % ruderok

% koordinatak elmozdulas utan elemenkent abrazolashoz:
dcrd=crd+sc*qmat; 
def(1:3:3*noe-2,:)=dcrd(N1,:);
def(2:3:3*noe-1,:)=dcrd(N2,:);
def(3:3:3*noe,:)=NaN;
if dim==2 % sikbeli feladat esete
    plot(def(:,1),def(:,2),'LineWidth',3);
else % trebeli feladat esete
    plot3(def(:,1),def(:,2),def(:,3),'LineWidth',3);
end