

function rama()

clear;clf 

Le=3; Lw=4;  

Lssu=Lw*3;

% stale materialowe
E= 2e7;  I=1e-4;  A=0.03;   

% parametry elementu 
pare= [E, A, I];

% wspolrzedne wezlow  
%         x   y    wez
xywez=[  0.0 0.0  % 1
         1.5 2.0  % 2
         3.0 4.0  % 3
         6.0 0.0  % 3
        ];

% macierz stopni swobody wezlow 
for i=1:Lw
  NrSSwz(i,1:3)= [(i-1)*3+1:i*3];
end

% macierz topologii 
%    elem wz.p wz.k
Top=[ 1    1    2
      2    2    3
      3    3    4
     ];

% macierz stopni swobody ES 
%                elem        p                k 
for e =1:Le    
  NrSSel(e, 1:7)= [e, NrSSwz(Top(e,2),:), NrSSwz(Top(e,3),:)];
end

% intensywnosc obciazenia ciaglego na elemencie 
%      px py  elem
Eq= [   0  0  % 1
        0  0  % 2
        4.8  -3.6  % 3
      ];


% Wektory wspolrzednych wezlow elementow
[Ex, Ey]= coordxtr( NrSSel, xywez, NrSSwz, 2 );

%  rysunek ramy
figure(1)
clf
eldraw2(Ex, Ey, [1,1,0]); hold on;  grid on;
axis([-1 8 -1 5]); 
%axis ij %xy


% Zerowanie macierzy 
F= zeros(Lssu,1);  K= zeros(Lssu,Lssu);

% sily wezlowe   
%F(4)= 100; 
%F(9)= 50; 

% macierze sztywnosci elementow i agregacja macierzy 
for i=1:Le
    [Ke,Fe]= beam2e( Ex(i,:), Ey(i,:), pare, Eq(i,:) );
    [K, F]= assem( NrSSel(i,:), K, Ke,  F, Fe );
end

F
% warunki brzegowe 
%    nrss  wart
bc= [ 1    0;  % u wz1
      2    0;  % v wz1
      3    0;  % fi wz1
     10    0
     11    0;  % v wz3
     12    0;  % v wz3
     ];

% obliczenie wektora przemieszczen i reakcji
[Q, R]= solveq( K, F, bc );
% wydruki wektora przemieszcze? i rekcji
%QiR=[Q, R]

 fprintf('\n\n W Y N I K I:\n');
 fprintf(' ============\n\n');
 fprintf('      Przemieszczenia                             Reakcje\n');
 fprintf(' -----------------------                    ----------------------\n');
 fprintf(' wz      u           v          fi            Ru       Rv       Rfi\n');
 for wz=1:Lw
  fprintf(' %2d %11.4e %11.4e %11.4e   %8.3f %8.3f %8.3f\n',...
  wz, Q((wz-1)*3+1:wz*3), R((wz-1)*3+1: wz*3) );
 end


% przemieszczenia wezlow dla poszczegolnych elementow
Qe= extract( NrSSel, Q );

N=5;
% powrot do elementu
for i=1:Le
  fe(:, :, i)=beam2s( Ex(i,:), Ey(i,:), pare, Qe(i,:), Eq(i,:), N);
end

% rysunek przemieszczeń 
figure(1)
eldisp2( Ex, Ey, Qe, [1,4,1] );
%axis([-1 8 -1 5]);  
%axis ij
title('Przemieszczenia');
hold off

plotpar=[2  1];

% Wykresy sil przekrojowych
% Sily podluzne
figure(2)
clf
scal=0.025;
for i=1:Le
     eldia2(Ex(i,:), Ey(i,:), fe(:,1,i), plotpar, scal);
     hold on
end
grid on;
title('Sily podluzne')
%%N=fe(:,1,:)

% Sily poprzeczne
figure(3)
clf
scal=0.025;
for i=1:Le
     eldia2(Ex(i,:), Ey(i,:), fe(:,2,i), plotpar, scal);
     hold on
end
grid on;
title('Sily poprzeczne')
%T=fe(:,2,:)

% Momenty
figure(4)
clf
scal=0.025;
for i=1:Le
     eldia2(Ex(i,:), Ey(i,:), fe(:,3,i), plotpar, scal);
     hold on
end
grid on;
title('Momenty zginajace')
%M=fe(:,3,:)

%fe 
 fprintf('\n      Sily przekrojowe\n');
 fprintf('      -----------------\n');
for e =1:Le
  fprintf('\n element nr%2d\n', e);
  fprintf('nrp      N         T         M\n');
  for i=1:N
    fprintf('%2d  %8.3f  %8.3f  %8.3f\n',i, fe(i,:,e));
  end
end 

%% =========================================================
% powrót do elementow - obliczenie sil przywezlowych
 fprintf('\n      Sily przywezlowe\n');
 fprintf('      -----------------\n');
 for i = 1:Le
    
  b=[ Ex(i,2)-Ex(i,1); Ey(i,2)-Ey(i,1) ];
  L=sqrt(b'*b);
  c= b(1)/L; s=b(2)/L;
 
  T=[c  s  0   0  0  0;
    -s  c  0   0  0  0;    
     0  0  1   0  0  0;
     0  0  0   c  s  0;
     0  0  0  -s  c  0;
     0  0  0   0  0  1];

% macierz sztywnosci ES w ukladzie lokalnym  
  ke=[E*A/L   0            0      -E*A/L      0          0 ;
          0   12*E*I/L^3   6*E*I/L^2  0   -12*E*I/L^3  6*E*I/L^2;
          0   6*E*I/L^2    4*E*I/L    0   -6*E*I/L^2   2*E*I/L;
        -E*A/L  0            0       E*A/L      0          0 ;
          0   -12*E*I/L^3 -6*E*I/L^2  0   12*E*I/L^3  -6*E*I/L^2;
          0   6*E*I/L^2    2*E*I/L    0   -6*E*I/L^2   4*E*I/L];
% wektor sil zastepczych od obciazenia elementowego w ukladzie lokalnym   
  qx= Eq(i,1); qy=Eq(i,2);
  ze=L*[qx/2 qy/2 qy*L/12 qx/2 qy/2 -qy*L/12]';
 
% transformacja przemieszczen wezlowych do ukladu lokalnego 
  qe= T*Qe(i,:)'; 

% sily przywezlowe 
  fe= ke*qe-ze;

  fprintf('\n element nr%2d\n', i);
  fprintf('nrw    fe(1)    fe(2)    fe(3)\n');
  fprintf('%2d  %8.3f %8.3f %8.3f \n%2d  %8.3f %8.3f %8.3f\n', ...
    Top(i,2), fe(1:3), Top(i,3), fe(4:6));
 
 end

%=================================================

end

