calfem_frame.sci

clc
clear
exec('solveq.sci'); 
exec('extract.sci'); 
exec('coordxtr.sci');
exec('beam2e.sci'); 
exec('beam2s.sci'); 
exec('assem.sci'); 
exec('eldraw2.sci'); 
exec('eldisp2.sci'); 
exec('beam2crd.sci'); 
//definicja wezlow i elmentow
nodes=[0 0; 0 2; 2 2; 2 0];
elements =[1 2; 2 3; 3 4]

// wlasnosci elementu skonczonego
E=200e9; 
b = 0.1
h=0.2
A=b*h; 
I = b*h^3/12
ep=[E A I];

//macierz stopni swobody i wspolrzedne elementow
for i=1:3
    np = elements(i,1)
    ne = elements(i,2)
    Ex(i,:) = [nodes(np,1),nodes(ne,1)];
    Ey(i,:) = [nodes(np,2),nodes(ne,2)];
    Edof(i,:) = [i, 3*np-2, 3*np-1,3*np, 3*ne-2, 3*ne-1, 3*ne];
end

//wykres 
Elnum=Edof(:,1);
eldraw2(Ex,Ey,Elnum);
set(gca(),"data_bounds",[-1,3,-1,3]);

// agregacja
LE=3;
LW=4;
LSSW=3;
LSSU=LW*LSSW;
F=zeros(LSSU,1);
K=zeros(LSSU,LSSU);
eq = [0 0; 0 1000; 0 0;]

for i=1:LE
    [Ke Fe]=beam2e(Ex(i,:),Ey(i,:),ep, eq(i,:));
    [K, F]=assem(Edof(i,:),K,Ke, F, Fe);
end

//warunki brzegowe, sila, podpory
//F(4)= 1000.0;
bc=[1 0; 2 0; 3 0; 10 0; 11 0; 12 0];
disp(F)

//rozwiazanie
[a,R]=solveq(K,F,bc);
Ed=extract(Edof,a);

//wykres zdeformowanej ramy
eldisp2(Ex,Ey,Ed,10000);

//wykres sil przekrojowych
figure(2)
kn=20
es = zeros(LE,kn,3)
edi = zeros(LE,kn,2)
eci = zeros(LE,kn)
for i=1:LE
    [ies,iedi,ieci]=beam2s(Ex(i,:),Ey(i,:),ep,Ed(i,:),eq(i,:),kn);
    es(i,:,:) = ies
    edi(i,:,:) = iedi
    eci(i,:) = ieci'
end
for i=1:LE
    eldia2(Ex(i,:),Ey(i,:),es(i,:,2)',eci(i,:)',1000);
end
set(gca(),"data_bounds",[-1,3,-1,3])