fem_heat.m

this fem program is written in matlab to simulate a room temperature of a museum under the source of

tic
%    the time dependent heat equation:
%
%      dU(x,y,t)/dt = Uxx(x,y,t) + Uyy(x,y,t) + F(x,y,t) in Omega x [0,T]
%
%    initial conditions
%
%      U(x,y,0) = U_INIT(x,y,0)
%
%     Dirichlet boundary conditions
%
%      U(x,y,t) = U_D(x,y,t) on Gamma_D
%
%    Neumann boundary conditions on the outward normal derivative:
%
%      Un(x,y) = G(x,y,t) on Gamma_N

  clear

  load -mat geom.m geom
[p,e,t]=initmesh(geom,'hmax',inf); 
[p,e,t]=refinemesh(geom,p,e,t); 
[p,e,t]=refinemesh(geom,p,e,t); 
[p,e,t]=refinemesh(geom,p,e,t); 
pdemesh(p,e,t);%subplot(3,2,2),
coordinates = p';
elements3 = t([1:3],:)';

[dirichlet,neumann]= boundaryedges(coordinates, e);

k = 2.6*1e-2; %conduction coefficient windows :units = W/k
 

%  Determine the bound and free nodes.
%
  BoundNodes = unique ( dirichlet );
  FreeNodes = setdiff ( 1:size(coordinates,1), BoundNodes );

  A = sparse ( size(coordinates,1), size(coordinates,1) );
  B = sparse ( size(coordinates,1), size(coordinates,1) );

  t_start = 0.0;
  t_final = 160;
  nt = 50;

  t = t_start;
  dt = ( t_final - t_start ) / nt;

  U = zeros ( size(coordinates,1), nt+1 );
%
%  Assembly.
%
  for j = 1 : size(elements3,1)
    A(elements3(j,:),elements3(j,:)) = A(elements3(j,:),elements3(j,:)) ...
      + stima3(coordinates(elements3(j,:),:));
  end

  for j = 1 : size(elements3,1)
    B(elements3(j,:),elements3(j,:)) = B(elements3(j,:),elements3(j,:)) ...
      + det ( [1,1,1; coordinates(elements3(j,:),:)' ] )...
      * [ 2, 1, 1; 1, 2, 1; 1, 1, 2 ] / 24;
  end
%
%  Set the initial condition.
%
  U(:,1) = u_init ( coordinates, t );
%
%  Given the solution at step I-1, compute the solution at step I.
%
 i = 1;
a = 10;
%   for i = 1 : nt
       while (a<19.9)

    t = ( ( nt - i ) * t_start   ...
        + (      i ) * t_final ) ...
        /   nt;

    b = sparse ( size(coordinates,1), 1 );

    u = sparse ( size(coordinates,1), 1 );
%
%  Account for the volume forces, evaluated at the new time T.
%
    for j = 1 : size(elements3,1)
      b(elements3(j,:)) = b(elements3(j,:)) ...
        + det( [1,1,1; coordinates(elements3(j,:),:)'] ) * ...
        k*dt * f ( sum(coordinates(elements3(j,:),:))/3, t ) / 6;
    end
%
%  Account for the Neumann conditions, evaluated at the new time T.
%
    for j = 1 : size(neumann,1)
      b(neumann(j,:)) = b(neumann(j,:)) + ...
        norm(coordinates(neumann(j,1),:) - coordinates(neumann(j,2),:)) * ...
        k*dt * g (U,i, sum(coordinates(neumann(j,:),:))/2, t ) / 2;
    end
%
%  Account for terms that involve the solution at the previous timestep.
%
    b = b + B * U(:,i);
%
%  Use the Dirichlet conditions, evaluated at the new time T, to eliminate the
%  known state variables.
%
    u(BoundNodes) = u_d ( coordinates(BoundNodes,:), t );

    b = b - ( dt *k* A  + B ) * u;
%
%  Compute the remaining unknowns by solving ( dt * A + B ) * U = b.
%
    u(FreeNodes) = ( dt * k*A(FreeNodes,FreeNodes) + B(FreeNodes,FreeNodes) ) ...
      \ b(FreeNodes);

    U(:,i+1) = u;
    a = min(U(:,i+1));
   i = i + 1;
       end
  ntt = i;
% end
  figure

  show ( elements3, coordinates, U, nt, t_start, t_final,ntt-1 ,dt);
colorbar


figure
a = find(coordinates(:,2)==0);
coord = sort(coordinates(a,1));
for i=1:5:ntt+1
    T=U(a,i);
    hold all
    plot(coord,T);
end
hold off
% a = [];
% for i = 1:nt+1
% a(i) = find(min(U(:,i)))
% end
%     
fprintf ( 'Time required to heat the room to a minimum temprature of %f degree celcius =  %f min\n',a, ntt*dt );