%MATLAB M-file pro simulaci 1D vedeni vlhkosti - nestacionarni uloha
%linearni aproximace relativni vlhkosti
%3.4.2007 vit.smilauer@fsv.cvut.cz

clear all; 
echo off;

%input data=======================================================

lambda = 0.1; %soucinitel vodivosti
l = .1; %delka konstrukce [m]
pocet_prvku = 10;
kapacita = 5000; %kapacita vlhkosti
Tau = 0.5;%0-explicitni, 1-implicitni, 0.5-Crank-Nicolson
delta_t = 1.;
t_max = 10;


%end of input data================================================

L = l/pocet_prvku; %delka prvku

%matice vodivosti
K_elem = lambda/L*...
[1 -1;
 -1 1];

C_elem = kapacita*L*...
[1./3. 1./6.;
 1./6. 1./3.];

B_elem = [1; 1];%integral pres B po prvku
N_elem = [1./2.; 1./2.];%integral pres B po prvku

%iniciace 
C_g=zeros(pocet_prvku+1);globalnich matic
K_g=zeros(pocet_prvku+1);
R_g=zeros(pocet_prvku+1);

%pocatecni podminky
h=ones(pocet_prvku+1,1);
h(pocet_prvku+1,1) = 0.;

t = 0;

while t < t_max

%lokalizace
for i=1:1:(pocet_prvku)
  K_g(i:i+1,i:i+1)=K_g(i:i+1,i:i+1)+K_elem(1:2,1:2);
  C_g(i:i+1,i:i+1)=C_g(i:i+1,i:i+1)+C_elem(1:2,1:2);
  R_g(i:i+1,1) = R_g(i:i+1,1) - lambda * (h(i+1) - h(i))/L* B_elem;%gradienty z predchoziho kroku
  R_g(i:i+1,1) = R_g(i:i+1,1) + kapacita * (h(i+1) + h(i))/2* N_elem;%hodnoty z predchoziho kroku
end;

%redukce soustavy
K_red=K_g(1:pocet_prvku,1:pocet_prvku);
C_red=C_g(1:pocet_prvku,1:pocet_prvku);
R_red=R_g(1:pocet_prvku,1);

%reseni soustavy
h_red = inv(Tau*K_red + 1/delta_t*C_red) * R_red;

t = t + delta_t;

h(1:pocet_prvku)=h_red;
h(pocet_prvku+1) = 0.0;

%graficky vystup

y=zeros(pocet_prvku+1,1);
y = (0:L:l)';
a2=plot(y,h);

pause(2);

end;

echo on

