function rama
E=27*10^6; % Modu³ Younga
A=0.1*0.1; % Pole powierzchni 
J=(0.1*0.1^3)/12; % Moment bezwladnosci
N=3; % Liczba elementów
n=4; % Liczba wêz³ów
Nss=3*n; 
X=[0 1 2 2]; %  Wspó³rzêdne x-owe wêz³ów
Y=[0 0 1 2]; % Wspó³rzêdne y-owe wêz³ów
inc1=[1 2;2 3;3 4]; % Relacja: element - wêz³y
inc2=zeros(N,7);
Xel=zeros(N,2);
Yel=zeros(N,2);
for i=1:N
    n1=inc1(i,1);
    n2=inc1(i,2);
    inc2(i,:)=[i 3*n1-2 3*n1-1 3*n1 3*n2-2 3*n2-1 3*n2];
    Xel(i,:)=[X(n1) X(n2)];
    Yel(i,:)=[Y(n1) Y(n2)];
end
close all
figure(1)
hold on
grid on
eldraw2(Xel, Yel, [ 1 1 1], 1:N)

K=zeros(Nss, Nss);
F=zeros(Nss,1);
F(5)=-20; % Obci¹¿enie skupione - si³a
F(9)=10; % Obci¹¿enie skupione - moment
q=zeros(N,2);
q(2,1)=-5; % Obci¹¿enie ci¹g³e styczne elementu = 2
q(2,2)=5; % Obci¹¿enie ci¹g³e normalne elementu = 2
q(3,2)=10; % Obci¹¿enie ci¹g³e normalne elementu = 3
for i=1:N
    [Kel,Fel]=beam2e(Xel(i,:), Yel(i,:), [E A J],q(i,:));
    [K,F]=assem(inc2(i,:), K, Kel, F, Fel);
end
wb=[1 0; 2 0; 10 0;11 0;12 0]; % Warunki brzegowe (numer ss - 0)
[U, R]=solveq(K, F, wb);
U
R
Uel=extract(inc2, U);
m=10;
sk=scalfact2(Xel, Yel, Uel, 0.35);
eldisp2(Xel, Yel, Uel, [2 2 1], sk)
title('przemieszczenia')
Sn=zeros(N,m);
Sq=zeros(N,m);
Sm=zeros(N,m);
for i=1:N
    z=beam2s(Xel(i,:), Yel(i,:), [E A J], Uel(i,:), q(i,:), m);
    Sn(i,:)=z(:,1)';
    Sq(i,:)=z(:,2)';
    Sm(i,:)=z(:,3)';
end
figure(2)
subplot(2,2,1);
hold on
grid on
sk=scalfact2(Xel, Yel, Sn, 0.5);
for i=1:N
    eldia2(Xel(i,:), Yel(i,:), Sn(i,:)', [2,1], sk);
end
title('si³y pod³u¿ne')

subplot(2,2,2);
hold on
grid on
sk=scalfact2(Xel, Yel, Sq, 0.5);
for i=1:N
    eldia2(Xel(i,:), Yel(i,:), Sq(i,:)', [2,1], sk);
end
title('si³y poprzeczne')

subplot(2,2,3);
hold on
grid on
sk=scalfact2(Xel, Yel, Sm, 0.5);
for i=1:N
    eldia2(Xel(i,:), Yel(i,:), Sm(i,:)', [2,1], sk);
end
title('moment zginaj¹cy')


