%function Heat_Conduction_transient( )%*****************************************************************************80
%
%% Applies the finite element method to the heat equation.
%
%    The user supplies datafiles that specify the geometry of the region
%    and its arrangement into triangular elements, and the location and
%    type of the boundary conditions, which can be any mixture of Neumann and
%    Dirichlet, and the initial condition for the solution.
%
%    Note that only uses triangular elements in this version.
%
%    The unknown state variable U(x,y,t) is assumed to satisfy
%    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]
%
%    with initial conditions
%
%      U(x,y,0) = U_INIT(x,y,0)
%
%    with Dirichlet boundary conditions
%
%      U(x,y,t) = U_D(x,y,t) on Gamma_D
%
%    and Neumann boundary conditions on the outward normal derivative:
%
%      Un(x,y) = G(x,y,t) on Gamma_N
%
%    If Gamma designates the boundary of the region Omega,
%    then we presume that
%
%      Gamma = Gamma_D + Gamma_N
%
%    but the user is free to determine which boundary conditions to
%    apply.  Note, however, that the problem will generally be singular
%    unless at least one Dirichlet boundary condition is specified.
%
%    The code uses piecewise linear basis functions for triangular elements.
%
%    The user is required to supply a number of data files and MATLAB
%    functions that specify the location of nodes, the grouping of nodes
%    into elements, the location and value of boundary conditions, and
%    the right hand side function in the heat equation.  Note that the
%    fact that the geometry is completely up to the user means that
%    just about any two dimensional region can be handled, with arbitrary
%    shape, including holes and islands.
%
%
%  Local Parameters:
%
%    Local, real DT, the size of a single time step.
%
%    Local, integer NT, the number of time steps to take.
%
%    Local, real T, the current time.
%
%    Local, real T_FINAL, the final time.
%
%    Local, real T_START, the initial time.
%%  Read the nodal coordinate data file.
%load coordinates.dat;
%
%  Read the triangular element data file.
%load elements3.dat;
%
%  Read the Neumann boundary condition data file.
%  I THINK the purpose of the EVAL command is to create an empty NEUMANN array
%  if no Neumann file is found.
%eval ( 'load neumann.dat;', 'neumann=[];' );
%
%  Read the Dirichlet boundary condition data file.
%  There must be at least one Dirichlet boundary condition.
%load dirichlet.dat;
%
%  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 = 1;nt = 10;% nt means number of time intervals instead of timestepst = 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,:),:));endfor 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.
%for i = 1 : ntt = ( ( 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,:),:)'] ) * ...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),:)) * ...dt * g ( 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 * A  + B ) * u;
%
%  Compute the remaining unknowns by solving ( dt * A + B ) * U = b.
%u(FreeNodes) = ( dt * A(FreeNodes,FreeNodes) + B(FreeNodes,FreeNodes) ) ...\ b(FreeNodes);U(:,i+1) = u;endshow ( elements3, coordinates, U, t_start, t_final );

function value = u_init ( xy, t )%*****************************************************************************80
%
%% U_INIT sets the initial condition for the state variable.
%
%  Discussion:
%
%    The user must supply the appropriate routine for a given problem
%
%    In many problems, the initial time is 0.  However, the value of
%    T is passed, in case the user wishes to use this same routine to
%    evaluate, for instance, the exact solution.
%
%
%  Parameters:
%
%    Input, real XY(N,M), contains the M-dimensional coordinates of N points.
%    N is probably the total number of points, and M is probably 2.
%
%    Input, real T, the initial time.
%
%    Output, VALUE(N), contains the value of the solution U(X,Y,T).
%n = size ( xy, 1 );value = zeros ( n, 1 );middle = floor ( n / 2 );value ( middle ) = 1.0;