function kratownicaZwykladu()


%      ^ y
%      |
%      |   e3
%  wz3 o + + + + + o wz4
%      + +
%    e1+   +e2
%      +     +
%      +       +
%      +         +
%      o wz1       o wz2---------------->x


 LssWez = 2;                % liczba stopni swobody elentu
 LwezEl = 2;                % liczbawezlow elementu
 Lel =  3;                  % liczba elementów skoczonych
 Lwez= 4;                   % liczba wêz³ów  kratownicy 
 LssUk = LssWez*Lwez;       % liczba stopni swobody ukladu

% wspó³rzêdne wezlow
%         x  y    wez
 xyWez = [0  0   % 1
          3  0   % 2
          0  4   % 3
          3  4]; % 4
% relacje przylegania
%           wz1  wz2  elem
 RelPrzyl = [1     3   % 1
             2     3   % 2
             3     4]; % 3                  

%--- materia³ , przekrój elementu
 E = 10^4; %kN
 A = 1;
 ParEl = [E A; E A; E A];
 
%-- Wektory wspolrzednych i  stopni swobody  wezlow  elementow
 for el=1:Lel
   for wz=1:2
     wez = RelPrzyl(el,wz); 
     Ex(el,wz)= xyWez(wez,1);
     Ey(el,wz)= xyWez(wez,2);
     ss(wz)= wez*2-1;
   end
   NrSSel(el,1:5)= [el, ss(1), ss(1)+1, ss(2), ss(2)+1]; 
 end

figure(1)
rysunekKratownicy(xyWez, Ex, Ey);

figure(2)
eldraw2(Ex,Ey,[1,4,1]); hold on; axis equal;

%--- obci±zenie elem
 qel=[ 30, 0, -50];  %kN/m

 K=zeros(LssUk); P=zeros(LssUk,1);  P(5)= 35; %kN

%--  agregacja
 for el=1:Lel
    [Kel] = bar2e(Ex(el,:), Ey(el,:), ParEl(el,:));
    Zel(:,el)= zastepnikiWezloweObciazeniaEl(Ex(el,:), Ey(el,:), qel(el) );
    [K, P] = assem(NrSSel(el,:), K, Kel, P, Zel(:,el));
 end


%-- warunki brzegowe
%    ss  wart.
 bc= [1    0; 
      2    0; 
      3    0; 
      4   -0.01;
      7    0; 
      8    0];
 [Q, R]=solveq(K,P,bc);

 fprintf('\n\n W Y N I K I:\n');
 fprintf(' ============\n\n');
 fprintf('      Przemieszczenia          Reakcje\n');
 fprintf(' ----------------------------------------\n'); 
 fprintf(' wz      u       v           ru      rv\n'); 
 k=1;
 for i=1:Lwez
  fprintf(' %2d %8.4f %8.4f    %8.3f %8.3f\n', i, Q(k),Q(k+1), R(k), R(k+1));
  k=k+2;
 end

 figure(2)
 Qel= extract(NrSSel,Q);   % przemieszczenia w wezlach elementu skonczonego
 rysunekKratownicy(xyWez, Ex, Ey, Qel);

%-- rysunek 
% figure(3)
% grid on;
% axis([ -0.5 3.5  -0.5 4.5])
% eldraw2(Ex, Ey, [2,4,1] ); hold on; %axis equal;
% Qel= extract(NrSSel,Q);  % przemieszczenia w wezlach elementu skonczonego
% eldisp2(Ex, Ey, Qel, [1,2,1]); hold off;
% %set(gca,'XDir','reverse')


%--- powrot do elmentu
%--  obliczenie si³ przywêz³owych w elementach (uk³ad lokalny)
 fprintf('\n\n Si³y przywêzlowe\n');
 fprintf(' ---------------------\n'); 
 fprintf(' el     w1       w2   \n'); 
 for el=1:Lel
  fel(el,:) = silyPrzywezloweEl(Ex(el,:),Ey(el,:),Qel(el,:),Zel(:,el),ParEl(el,:))';
  fprintf(' %2d  %8.3f %8.3f  \n', el, fel(el,:));
 end
  fprintf('\n\n');

figure(4)
grid on;
axis([ -0.5 3.5  -0.5 4.5])
for i=1:Lel
 eldia2a(Ex(i,:),Ey(i,:),fel(i,:)',[3,1], 0.01);
end
title('sily w pretach')



%---------------- Koniec zadania ---------------------

function rysunekKratownicy(xyWez, Ex, Ey, Qe)
 if ~((nargin==3)|(nargin==4))
    disp('??? Nieprawidlowe wolanie funkcji: rysunekKratownicy!')
    exit
 end

% grid
 
switch(nargin)
 case 3
   mx(2) =max(max(Ex)); my(2) =max(max(Ey));
   mx(1) =min(min(Ex)); my(1) =min(min(Ey));
   for el = 1:length(Ex)
     plot (Ex(el,:),Ey(el,:),'b-', Ex(el,:),Ey(el,:),'ro','LineWidth',3); 
     hold on;
     x = Ex(el,1)+Ex(el,2); y = Ey(el,1)+Ey(el,2);
     d = 0.07;
     h = text(x/2+d,y/2+d,int2str(el));
     set(h,'fontsize',10);
     axis([ mx(1)-0.5 mx(2)+0.5  my(1)-0.5 my(2)+0.5])
     grid on
   end
 case 4
   dxmx = max(max(Ex')-min(Ex')); dymx = max(max(Ey')-min(Ey'));
   dlmx = max(dxmx,dymx); Qemx = max(max(abs(Qe))); krel = 0.1;
   sk = krel*dlmx/Qemx;
   qwx = Ex+sk*Qe(:,[1,3]); qwy = Ey+sk*Qe(:,[2,4]);
   mx(2) = max(max(qwx));   my(2) = max(max(qwy));
   mx(1) = min(min(qwx));   my(1) = min(min(qwy));
   for el = 1:length(Ex)
    plot(Ex(el,:),Ey(el,:),'b--', Ex(el,:),Ey(el,:),'bo','LineWidth',3);
    hold on;
    plot(qwx(el,:),qwy(el,:),'r-',qwx(el,:),qwy(el,:),'ro','LineWidth',3); 
    axis([ mx(1)-0.5 mx(2)+0.5  my(1)-0.5 my(2)+0.5])
    grid on
   end
 otherwise
   disp('!! Nieznany rysunek');
 end
%---koniec funkcji rysunekKratownicy--------------



function Zel= zastepnikiWezloweObciazeniaEl(Ex, Ey, qel )
%-------------------------------------------------------
% INP: 
%      Ex(x_1, x_2), Ey(y_1, y_2)  - wspó³rzedne wezlow elementu
%      qel -  obci±zenie rozlozone 
%
% OUT: Zel(1:4,1)   
%-------------------------------------------------------
    l = [ Ex(2)-Ex(1); Ey(2)-Ey(1) ];
    le = sqrt(l'*l);
    pel(1:2,1) = qel*[le-0.5*le; 0.5*le];
    s = l(2)/le; c = l(1)/le;
    t = [c s 0 0
        0 0 c s];
    Zel = t'*pel;
%---koniec funkcji zastepnikiObciazenia-----------


function fel = silyPrzywezloweEl(Ex, Ey, Qel, Zel, ParEl)
%-------------------------------------------------------
% INP: 
%      Ex(x_1, x_2), Ey(y_1, y_2)  - wspó³rzedne wezlow elementu
%      Qel(1:4, 1)  - uogólnione przemieszczenia wezlow elementu       
%      Zel((1:4, 1) - zastepniki wezlowe obciazenia elementu
%      ParEl(E,A)   - parametry elementu
%
% OUT: fel(1:2, 1) (w ukladzie lokalnym elementu)  
%-------------------------------------------------------

%-- I-sposob na wykladzie
    [Kel] = bar2e(Ex, Ey, ParEl);
    Fel = Kel*Qel' - Zel;
    l = [Ex(2)-Ex(1); Ey(2)-Ey(1)];
    le = sqrt(l'*l);
    s = l(2)/le; c= l(1)/le;
    t =[c s 0 0
        0 0 c s];
    fel = t*Fel;
%---koniec funkcji  silyPrzywezlowe---------------





