%MATLAB M-file pro prvky s linearni a kvadratickou aproximaci posunu
%23.10.2005 vit.smilauer@fsv.cvut.cz
%modified 17.10.2006 vit.smilauer@fsv.cvut.cz

%//////////
%----------
%    |    | x
%    |    v  
%    |
%    |
%    |
%    |
%    |
%    |
%    V F

clear all; 
echo off;

%input data=======================================================

E = 210.e+6;  %ocel [kPa]
l = 600.; %delka lana
A = 0.01; %plocha rezu [m2]
gama = 78.5; %merna tiha [kN/m3]
pocet_prvku = 1
konc_sila = 210; %kN

%end of input data================================================

L = l/pocet_prvku; %delka prvku
f = A*gama; %zatizeni [kN/m']

%presne reseni
for i = 0:1:2*pocet_prvku
    vzdalenost= i * L/2.;
    teor_reseni(i+1,1) = vzdalenost;
    teor_reseni(i+1,2) = (konc_sila + f*l)/E/A * vzdalenost - f /(2*E*A) * vzdalenost^2;
end;

%matice tuhosti linearni aproximace
k_lin = E*A/L*...
[1 -1;
 -1 1];

f_lin = f*L/2*...
[1;
1];


%matice tuhosti kvadraticka aproximace
k_quad = E*A/L*...
[7./3. -8./3. 1./3.;
-8./3. 16./3.  -8./3.;
1./3.  -8./3. 7./3.];

%f_quad = f*L/6*...
%[2;
%2;
%2];

f_quad = f*L/6*...
[1;
4;
1];


%INICIACE
K_lin_g=zeros(pocet_prvku+1);
F_lin_g=zeros(pocet_prvku+1,1);

K_quad_g=zeros(2*pocet_prvku+1);
F_quad_g=zeros(2*pocet_prvku+1,1);

%LOKALIZACE
for i=1:1:(pocet_prvku)
  K_lin_g(i:i+1,i:i+1)=K_lin_g(i:i+1,i:i+1) + k_lin(1:2,1:2);
  F_lin_g(i:i+1,1)=F_lin_g(i:i+1,1) + f_lin(1:2,1);
end;
%pridej F do posledniho uzlu v lin. aproximacich
F_lin_g(pocet_prvku+1,1)=F_lin_g(pocet_prvku+1,1) + konc_sila;

%v kvadratickem prvku jsou koncove posuny 1,3 (2 je uprostred)
for i=1:2:(2*pocet_prvku)
  K_quad_g(i:i+2,i:i+2)=K_quad_g(i:i+2,i:i+2) + k_quad(1:3,1:3);
  F_quad_g(i:i+2,1)=F_quad_g(i:i+2,1)+f_quad(1:3,1);
end;
%pridej F do posledniho uzlu v kvadr. aproximacich
F_quad_g(2*pocet_prvku+1,1) = F_quad_g(2*pocet_prvku+1,1) + konc_sila;

%redukce globalni matice tuhosti vlivem vetknuti (zrusit prvni radek a sloupec)
K_lin_red=zeros(pocet_prvku);
K_lin_red(1:pocet_prvku,1:pocet_prvku) = K_lin_g(2:pocet_prvku+1,2:pocet_prvku+1);
F_lin_red(1:pocet_prvku,1) = F_lin_g(2:pocet_prvku+1,1);

K_quad_red=zeros(2*pocet_prvku);
K_quad_red(1:2*pocet_prvku,1:2*pocet_prvku) = K_quad_g(2:2*pocet_prvku+1,2:2*pocet_prvku+1);
F_quad_red(1:2*pocet_prvku,1) = F_quad_g(2:2*pocet_prvku+1,1);

%reseni soustavy
r_lin = inv(K_lin_red) * F_lin_red;
r_quad = inv(K_quad_red) * F_quad_red;

%potencialni energie linearnich prvku dle "energie pruzin"
E_pot_lin = 0.;
for i=1:1:pocet_prvku 
   if (i==1)
       temp = 0;
   else
       temp = r_lin(i-1);
   end;
    
   E_pot_lin = E_pot_lin + 0.5*E*A/L*(r_lin(i) - (temp))^2 - F_lin_red(i) * r_lin(i);
end;
    
%klasicka definice potencialni energie MKP pro linearni prvek
E_pot_lin_MKP = 0.5*r_lin'*K_lin_red*r_lin - r_lin' * F_lin_red;


%klasicka definice potencialni energie MKP pro kvadraticky prvk
E_pot_quad_MKP = 0.5*r_quad'*K_quad_red*r_quad - r_quad' * F_quad_red;

%tisk posunu a energii
r_lin
E_pot_lin
E_pot_lin_MKP

r_quad
E_pot_quad_MKP

%graficky vystup
clear presne_reseni

for vzdalenost=0:l/20:l
   presne_reseni(vzdalenost*20/l+1,1) = (konc_sila + f*l)/E/A * vzdalenost - f /(2*E*A) * vzdalenost^2;
end;

y=zeros(20,1);
y = (0:-l/20:-l)';
a1=plot(presne_reseni,y); hold on;
clear y;
y = (-L:-L:-l)';
a2=plot(r_lin,y,'rO');
clear y;
y = (-L/2:-L/2:-l)';
a3=plot(r_quad,y,'kX');hold off;legend('presne','linearni','kvadraticke');



echo on

