krata2.m

function krata2()
 
% rozwiazanie MES
% bez uzycia pakietu CALFEM 
%
% opracował: Piotr Pluciński
 
% współrzędne węzłów
wez=[0,0 ; 3,0 ; 0 4];
 
% macierz topologii 
top=[1 2; 2 3; 1 3] ;
 
% stałe materiałowe
EA=1e4; 
 
obc_ciag=[2,0] % obciążenie ciągłe [nr_elementu wartość]
 
 
KG=zeros(2*size(wez,1));
Z=zeros(2*size(wez,1),1);
W=zeros(2*size(wez,1),1);
 
W(6)=-10; 
 
wb=[1 0 ; 2 0 ; 4 -0.001 ; 5 0 ] ;  % waunki brzegowe
 
for i=1:size(top,1)
    [Te,Le]=trans(wez,top(i,:)); % wyliczenie macierzy transformacji i długości elementu
    ke=mac_sztyw(EA,Le); % wyliczenie macierzy sztywności dla elementu
    Ke=Te'*ke*Te; % transformacja macierzy sztywności 
    KG=agre(KG,Ke,top(i,:));  % agregacja 
    if obc_ciag(1)==i     % obliczenie równoważnika obciążenia od obciążenia prostokątnego
        ze=zast(obc_ciag(2),Le); % obliczenie równoważnika 
        Ze=Te'*ze;
        Z=agre(Z,Ze,top(i,:));
    end
end
 
P=W+Z;
 
% obliczenie układów równań z uwzględnieniem warunków brzegowych
 
ss=[1:2*size(wez,1)]';
d=zeros(size(ss));
wbb=wb(:,1);
dp=wb(:,2);
ss(wbb)=[];
s=KG(ss,ss)\(P(ss)-KG(ss,wbb)*dp);
d(wbb)=dp;
d(ss)=s;
 
d                % wektor przemieszczeń (stopni swobody)
r=KG*d-P         % wektor reakcji 
 
% powrót do elementu
 
for i=1:size(top,1)
    [Te,Le]=trans(wez,top(i,:));
    de=wyciag(d,top(i,:));
    ke=mac_sztyw(EA,Le);
    de_lok=Te*de;
    fe=ke*de_lok
    if obc_ciag(1)==i     % obliczenie równoważnika obciążenia od obciążenia prostokątnego
        ze=zast(obc_ciag(2),Le); % obliczenie równoważnika 
        fe=fe-ze;
    end
    fel_kom(:,i)=fe;
end
fel_kom
 
clf 
for i=1:size(fel_kom,2)
    wez_el=[wez(top(i,1),1) wez(top(i,1),2); wez(top(i,2),1) wez(top(i,2),2)];
    rys_wykr(wez_el,fel_kom(:,i))
end
 
 
function rys_wykr(x,p)
 
ps=p*0.2; % skalowanie rysunku 
 
hold on
dx=x(2,1)-x(1,1);
dy=x(2,2)-x(1,2);
L=sqrt(dx^2+dy^2);
cos=dx/L;
sin=dy/L;
axis equal
axis off
% axis([-1.5 3.5 -0.5  4.5]);
plot([x(2,1) x(2,1)-ps(2)*sin x(1,1)+ps(1)*sin x(1,1)], [ x(2,2) x(2,2)+ps(2)*cos x(1,2)-ps(1)*cos x(1,2)],'-r')
plot([x(1,1) x(2,1)], [x(1,2) x(2,2)],'-k','LineWidth',2)
 
ppx=x(1,1)+ps(1)*sin;
ppy=x(1,2)-ps(1)*cos;
pkx=x(2,1)-ps(2)*sin;
pky=x(2,2)+ps(2)*cos;
 
    if dx<=0 
        if dy <0
 
            text(ppx,ppy,num2str(-p(1)),'HorizontalAlignment','left','VerticalAlignment','bottom')
            text(pkx,pky,num2str(p(end)),'HorizontalAlignment','right','VerticalAlignment','top')
        else
            text(ppx,ppy,num2str(-p(1)),'HorizontalAlignment','left','VerticalAlignment','top')
            text(pkx,pky,num2str(p(end)),'HorizontalAlignment','right','VerticalAlignment','bottom')
        end
    else
        if dy <0
            text(ppx,ppy,num2str(-p(1)),'HorizontalAlignment','right','VerticalAlignment','bottom')
            text(pkx,pky,num2str(p(end)),'HorizontalAlignment','left','VerticalAlignment','top')
        else
            text(ppx,ppy,num2str(-p(1)),'HorizontalAlignment','right','VerticalAlignment','top')
            text(pkx,pky,num2str(p(end)),'HorizontalAlignment','left','VerticalAlignment','bottom')
        end
    end
 
 
function de=wyciag(d,wz_e)
    idx=[wz_e(1)*2-1:wz_e(1)*2,wz_e(2)*2-1:wz_e(2)*2];
    de=d(idx);
 
function X=agre(X,xe,wz_e)
idx=[wz_e(1)*2-1:wz_e(1)*2,wz_e(2)*2-1:wz_e(2)*2];
if size(xe,2)>1 
    X(idx,idx)=X(idx,idx)+xe;
else
    X(idx)=X(idx)+xe;
end
 
 
function K=mac_sztyw(EA,L)
K=EA/L*[ 1 -1 ; -1 1];
 
function z=zast(p,L)
z=[p*L/2; p*L/2];
 
function [T,L]=trans(wez,wz_e)
dx=wez(wz_e(2),1)-wez(wz_e(1),1);
dy=wez(wz_e(2),2)-wez(wz_e(1),2);
L=sqrt(dx^2+dy^2);
cos=dx/L;
sin=dy/L;
T=[cos sin 0 0; 0 0 cos sin];