%M-file pro demonstraci pruhybu konce konzoly (shear locking)

clear all; 
echo off;

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

E = 1.e+4;
b = 1.;
h = 1.;
l = 25.; %delka konzoly
k = 5./6.;
G = 1.e+6;
P = 10.;
pocet_prvku = 10

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

L = l/pocet_prvku; %delka prvku
A = b*h;
I = 1./12.*b*h^3;

k_b = E*I/L*...
[0  0 0 0;
 0  1 0 -1;
 0  0 0 0;
 0 -1 0 1];

%full solution of FE
k_s = k*G*A/L*...
[1  L/2   -1   L/2;
L/2 L^2/3 -L/2 L^2/6;
-1  -L/2  1    -L/2;
L/2 L^2/6 -L/2 L^2/3];

%reduced integration of FE
k_s_red = k*G*A/L*...
[1  L/2   -1   L/2;
L/2 L^2/4 -L/2 L^2/4;
-1  -L/2  1    -L/2;
L/2 L^2/4 -L/2 L^2/4];

%k_s=k_s_red; %uncomment for reduced integration

staticky_max_pruhyb = P*l^3/3/E/I;

K_g=zeros(2*pocet_prvku+2);

%lokalizace
for i=1:2:(2*pocet_prvku)
  K_g(i:i+3,i:i+3)=K_g(i:i+3,i:i+3)+k_b(1:4,1:4)+k_s(1:4,1:4);
end;

%redukce globalni matice tuhosti vlivem vetknuti (zrusit posledni 2 radky a
%2 sloupce)
K_red=zeros(2*pocet_prvku);
K_red(1:2*pocet_prvku,1:2*pocet_prvku) = K_g(1:2*pocet_prvku,1:2*pocet_prvku);

%vektor zatizeni
F = zeros(2*pocet_prvku,1);
F(1,1) = P;

%reseni soustavy
r = inv(K_red) * F;

%tisk posunu konce a sil na koncovem prvku
r(1,1)
k_elem=k_b(1:4,1:4)+k_s(1:4,1:4);

if pocet_prvku ~= 1
    k_elem*r(1:4,1)
end

echo on

