calfem_heat_transfer.sci

clc
clear
exec('assem.sci'); 
exec('coordxtr.sci');
exec('eldraw2.sci');
exec('extract.sci');
exec('solveq.sci'); 
exec('flw2qs.sci'); 
exec('statcon.sci'); 
exec('flw2qe.sci');
exec('flw2te.sci');
exec('flw2ts.sci');
exec('elflux2.sci');
exec('eliso2.sci');

wn = 8; x0 = 0.0; xn = 1.0;
hn = 10; y0 = 0.0; yn = 3.0;
xx = linspace(x0, xn, wn)
yy = linspace(y0, yn, hn)
it = 1
LW = wn*hn
nodes = zeros(LW,2)
for i=1:wn
    for j=1:hn
        nodes(it,1) = xx(i)
        nodes(it,2) = yy(j)
        it = it + 1
    end
end
LE = (wn-1)*(hn-1)

elements = zeros(LE,4)
it = 1
for i=1:wn-1
    for j=1:hn-1
        elements(it,:) = [(i-1)*hn+j, i*hn+j, i*hn+j+1, (i-1)*hn+j+1]
        it = it + 1
    end
end

for i=1:LE
    n1 = elements(i,1)
    n2 = elements(i,2)
    n3 = elements(i,3)
    n4 = elements(i,4)
    Ex(i,:) = [nodes(n1,1),nodes(n2,1),nodes(n3,1),nodes(n4,1)];
    Ey(i,:) = [nodes(n1,2),nodes(n2,2),nodes(n3,2),nodes(n4,2)];
    Edof(i,:)=[i, n1, n2, n3, n4];
end

eldraw2(Ex,Ey);

t0 = 20.0
t1 = 100.0
it0 = 1:hn:LW
it1 = hn:hn:LW
n = length(it0)
m = length(it1)

bc = zeros(n+m, 2)
bc(1:n,1) = it0
bc(1:n,2) = t0
bc(n+1:$,1) = it1
bc(n+1:$,2) = t1

disp(bc)

LSSW = 1
LSSU = LW*LSSW
K=zeros(LSSU,LSSU)
F=zeros(LSSU,1)

ep=1; D=[1 0;0 1];

for i=1:LE
    Ke=flw2qe(Ex(i,:),Ey(i,:),ep,D);
    K=assem(Edof(i,:),K,Ke);
end;

[a,Q]=solveq(K,F,bc)

Ed=extract(Edof,a);
for i=1:LE
    Es(i,1:2)=flw2qs(Ex(i,:),Ey(i,:),ep,D,Ed(i,:));
end

elflux2(Ex,Ey,Es,0.01);