
%% opracowala Ewa Pabisek

function belka_symbolicznie()
format compact
clear,clc,clf

%  /|1     2     3     4     5
%  /|o=====o=====o=====o=====o  ----->x
%  /||
%    |
%    V y

syms x real
syms l EIs real  

 H1= 1 -3*(x/l)^2 + 2*(x/l)^3;
 H2= x*(1-x/l).^2;
 H3= 3*(x/l)^2 -2*(x/l)^3;
 H4= x*((x/l)^2 - x/l);

 N= [H1, H2, H3, H4];   % macierz funkcji ksztaltu 

%%-----------  B= diff(diff(N,x), x) 
%%----------- Ke= EI int_0^l B' B dx
%% macierz sztywnosci obliczona symbolicznie

 B= diff(diff(N, x), x)     
 Ke_sym= EIs * int(B'*B, 0, l) 

%% wektor sl zastepczych od obciazenia rozlozonego  dla q(x)=1;
 py_e_sym= int(N', 0, l)   % bedzie mnozone przez wartosci aktualna py

%% ============ PRE-PROCESSING ======================================
%
%%  dane dla belki
%% B25  b=30 h=50 cm 
%%  Jednostki:  kN, m

 E= 3.0e7;   I= 3.125e-3;   EI= E*I;  
 L= 4;           
 py= 3;           


%% dyskretyzacja  
 Lel = 4;         %  liczba elementow      
 Lwez= Lel+1;     %  liczba wezlow           
 le= L/Lel;       %  dlugosc elementu skonczonego

% v, fi   stopnie swobody wezla
  lssWez= 2;         % liczba stopni swobody wezla
  LssUk= Lwez*lssWez;  %liczba stopni swobody ukladu

%%  wspolrzedne wezlow
 xWez= linspace(0, L, Lwez); 

%%  macierz stopni swobody elementow
%       el  wz1    wz2
 Edof= [ 1  1  2   3  4
         2  3  4   5  6
         3  5  6   7  8
         4  7  8   9 10
        ];

 Ke= double(subs(Ke_sym, {l,EIs}, {le,EI}));
 pe= double(subs(py_e_sym, l, le));    % pe bedzie mnozone przez py_e
                                       % podczas agregacji 

%% wektor obciazen elementowych
 py_e(1:Lel)= py;

%% ============= OBLICZENIA==========================================
 
%% agregacja  Ke i ze        
 K= zeros(LssUk, LssUk);  % globalna macierz sztywnosci
 z= zeros(LssUk, 1);      % wektor obciazen wezlowych

 for ne=1:Lel
     ze = pe* py_e(ne);       % wektor od obciazenia elementowego
     edof= Edof(ne, :);        
     [K,z]= assem(edof, K, Ke, z, ze);       % funkcja  CALFEMa
 end

%% warunki brzegowe                 
%      nrss  wart.
 bc= [  1     0;  
        2     0; 
    ];

%% Obliczenie wektora przemieszczen d i reakcji r
 [d, r]= solveq(K, z, bc);  

 fprintf('\n\n W Y N I K I:\n');
 fprintf(' ============\n\n');
 fprintf('      Przemieszczenia               Reakcje\n');
 fprintf(' -----------------------       -----------------\n');
 fprintf(' wz      v           fi            rv       rfi\n');
 for wz=1:Lwez
  fprintf(' %2d %11.4e %11.4e    %8.3f %8.3f\n',...
  wz, d((wz-1)*2+1:wz*2), r((wz-1)*2+1: wz*2) );
 end

%% ============= POST-PROCESSING=====================================
figure(1);
 plot(xWez, zeros(Lwez,1),'k.-','LineWidth',3, 'MarkerSize',25);
 axis ij, hold on; grid on;
 mm= max(abs(d(1:2:end))); 
 xlim([-0.05*L 1.05*L])
 ylim([-0.1*mm mm+0.1*mm])
 title('v(x)')

 for ne=1:Lel
    % budowa interpolanta 
    n1= ne;   x1= xWez(n1);
    n2= n1+1; x2= xWez(n2);
    
    nrssel= Edof(ne,2:5);
    de= d(nrssel);
    v= subs(N,l,(x2-x1))*de; 
    v_m= matlabFunction(v);
    
    xe= linspace(x1,x2,11);
    v= v_m(xe-x1); 
    
    plot(xe,v,'b', 'LineWidth',2);
    plot(x2,v(end),'m.', 'MarkerSize',25);

    pause(1)

 end

%%-- powrot do elementu
%%-- sily przywezlowe
 fprintf('\n\n Sily przywezlowe \n');
 fprintf(' ----------------------------------------- -\n');
 fprintf(' el       Qi       Qj    |     Mi       Mj\n');

 for ne=1:Lel
    nrssel= Edof(ne,2:5);
    de= d(nrssel);
    ze = pe* py_e(ne);
%    ze= py_e(ne)*[le/2; le^2/12; le/2; -le^2/12];
    fe([1 3 2 4], ne)= Ke*de - ze;
    fel(1:4,ne)= Ke*de - ze;
    fprintf(' %2d   %8.3f %8.3f  | %8.3f %8.3f \n', ne, fe(:, ne));
 end
fprintf('\n\n');

 ryssil(xWez,fel,py_e)

end
%==================================================================

function ryssil(xnode,fe,pe)
dokl= 8;
X=[]; Q=[]; M=[];
for i= 1:length(xnode)-1
    xL= xnode(i);
    xP= xnode(i+1);
    r= fe(:,i);
    p= pe(i);
    x= linspace(xL,xP,dokl+1);
    X= [X x];
    q= r(1)+p*(x-xL);
    Q= [Q q];
    m= -r(1)*(x-xL) + r(2)-p*(x-xL).^2/2;
    M= [M m];
end

figure (11);
clf
subplot(211)
axis ij;
%plot([xnode(1),xnode(1)],1.0*[min(Q),max(Q)],'k')
area(X,Q,'FaceColor', [146/255 202/255 93/255]);
hold on
plot(xnode, zeros(length(xnode),1),'k.-','linewidth',3,'MarkerSize',25);
title('Sila poprzeczna')
xlim([-0.05*xnode(end) 1.05*xnode(end)])
grid on
if Q(1)>0
    text(X(1),Q(1),num2str(Q(1)),'HorizontalAlignment','left',...
        'VerticalAlignment','top')
else
    text(X(1),Q(1),num2str(-Q(1)),'HorizontalAlignment','left',...
        'VerticalAlignment','bottom')
end
if Q(end)>0
    text(X(end),Q(end),num2str(Q(end)),'HorizontalAlignment','right',...
        'VerticalAlignment','top')
else
    text(X(end),Q(end),num2str(-Q(end)),'HorizontalAlignment','right',...
        'VerticalAlignment','bottom')
end


subplot(212)
axis ij;
area(X,M,'FaceColor', [1 212/255 103/255]);
hold on
plot(xnode, zeros(length(xnode),1),'k.-','linewidth',3,'MarkerSize',25);
title('Moment zginajacy')
xlim([-0.05*xnode(end) 1.05*xnode(end)])
grid on
if M(1)>0
    text(X(1),M(1),num2str(M(1)),'HorizontalAlignment','left',...
        'VerticalAlignment','top')
else
    text(X(1),M(1),num2str(-M(1)),'HorizontalAlignment','left',...
        'VerticalAlignment','bottom')
end
if M(end)>0
    text(X(end),M(end),num2str(M(end)),'HorizontalAlignment','right',...
        'VerticalAlignment','top')
else
    text(X(end),M(end),num2str(-M(end)),'HorizontalAlignment','right',...
        'VerticalAlignment','bottom')
end

[wmin,nmin]=min(M);
if nmin~=1 && nmin~=length(M)
    text(X(nmin),wmin,num2str(-wmin),'HorizontalAlignment','center',...
        'VerticalAlignment','bottom')
end

[wmax,nmax]=max(M);
if nmax~=1 && nmax~=length(M)
    text(X(nmax),wmax,num2str(wmax),'HorizontalAlignment','center',...
        'VerticalAlignment','top')
end

end
