%                                                        
%   BLOCK2D_FD1     2-D Heat Conduction in a Rectangular Block
%                                                        
%   Numerical Soln using the Finite Difference Method for the following heat 
%   conduction problem:
%     uxx(x,y) + uyy(x,y)  = 0    with ux(0,y) = 0   u(Lx/2,y) = 0            
%                                       u(x,0) = 0  -kuy(x,Ly) = h(u - Te)
%     where  u(x,y) = T(x,y) - Tw   and   Te = Tinf - Tw  
%
%   Note that the PDE is not modelled explicitly  --  instead, discrete energy 
%   balances are performed for each unique nodal volume in the geometry of 
%   interest.
%   
%   The goal here is to compute and plot the discrete temperature distribution 
%   We also want to compute the heat loss or gain across each face of the 
%   rectangular block.  Only the right half of the block is modeled  --  a 
%   symmetry condition at x = 0 handles this (note that the heat flows given 
%   are only for the portion of the block modeled -- need to multiply by two to 
%   account for the whole geometry)).
%
%   This same problem is solved using the Separation of Variables method in 
%   file BLOCK2D_SOV1.M.
%
%   File prepared by J. R. White, UMass-Lowell (Dec. 2003)
%

%
%   getting started
      clear all,   close all,   nfig = 0;
%                                                                     
%   problem data                                                     
      Tw = 50;          % wall temperature on left, bottom, and right sides (C)
      Tinf = 100;       % air temperature on top of block (C)
      h = 100;          % heat transfer coeff (W/m^2-C)
      k = 1;            % thermal conductivity of block (W/m-C)
      Lx = 0.06;        % width of block in x direction (m)
      Ly = 0.09;        % height of block in y direction (m)
%
%   determine spatial grid for the half -block geometry
      Nx = 15;   dx = Lx/2/Nx;  dx2 = dx*dx;  x = 0:dx:Lx/2-dx;
      Ny = 36;   dy = Ly/Ny;    dy2 = dy*dy;  y = dy:dy:Ly;
%
%   allocate storage space  
      N = Nx*Ny;   A = zeros(N,N);   B = zeros(N,1);
%
%   setup matrix coefficients for all Central Nodes
      for j = 2:Ny-1
        for i = 2:Nx-1
          n = i+(j-1)*Nx;
          A(n,n-1) = dy2;   A(n,n+1) = dy2;          % left and right neighbors
          A(n,n) = -2*(dy2 + dx2);                   % self coupling
          A(n,n-Nx) = dx2;   A(n,n+Nx) = dx2;        % bottom and top neighbors
        end
      end
%
%   setup matrix coefficients for Left BC
      i = 1;
      for j = 2:Ny-1
        n = i+(j-1)*Nx;                       
        A(n,n+1) = dy2;                              % right neighbor 
        A(n,n) = -(dy2 + dx2);                       % self coupling
        A(n,n-Nx) = dx2/2;   A(n,n+Nx) = dx2/2;      % bottom and top neighbors
      end
%
%   setup matrix coefficients for Bottom BC
      j = 1;
      for i = 2:Nx-1
        n = i+(j-1)*Nx;
        A(n,n-1) = dy2;   A(n,n+1) = dy2;            % left and right neighbors
        A(n,n) = -2*(dy2 + dx2);                     % self coupling
        B(n) = -dx2*Tw;   A(n,n+Nx) = dx2;           % bottom and top neighbors
      end
%
%   setup matrix coefficients for Right BC
      i = Nx;
      for j = 2:Ny-1
        n = i+(j-1)*Nx;
        A(n,n-1) = dy2;   B(n) = -dy2*Tw;            % left and right neighbors
        A(n,n) = -2*(dy2 + dx2);                     % self coupling
        A(n,n-Nx) = dx2;   A(n,n+Nx) = dx2;          % bottom and top neighbors
      end
%
%   setup matrix coefficients for Top BC
      j = Ny;
      for i = 2:Nx-1
        n = i+(j-1)*Nx;
        A(n,n-1) = dy2/2;   A(n,n+1) = dy2/2;        % left and right neighbors
        A(n,n) = -(dy2 + dx2 + h*dx2*dy/k);          % self coupling
        A(n,n-Nx) = dx2;   B(n) = -h*dx2*dy*Tinf/k;  % bottom and top neighbors
      end
%
%   now set up the coefficients for the four corners
%   Lower Left
      i = 1;   j = 1;   n = i+(j-1)*Nx;  
      A(n,n+1) = dy2;                                % right neighbor
      A(n,n) = -(dy2 + dx2);                         % self coupling
      B(n) = -dx2*Tw/2;   A(n,n+Nx) = dx2/2;         % bottom and top neighbors
%   Lower Right
      i = Nx;   j = 1;   n = i+(j-1)*Nx;  
      A(n,n-1) = dy2;                                % left neighbor
      A(n,n) = -2*(dy2 + dx2);                       % self coupling
      A(n,n+Nx) = dx2;                               % top neighbor
      B(n) = -(dy2 + dx2)*Tw;                        % bottom and right neighbors
%   Upper Right
      i = Nx;   j = Ny;   n = i+(j-1)*Nx;  
      A(n,n-1) = dy2/2;                              % left neighbor
      A(n,n) = -(dy2 + dx2 + h*dx2*dy/k);            % self coupling
      A(n,n-Nx) = dx2;                               % bottom neighbor
      B(n) = -dy2*Tw/2 - h*dx2*dy*Tinf/k;            % right and top neighbors
%   Upper Left
      i = 1;   j = Ny;   n = i+(j-1)*Nx;  
      A(n,n+1) = dy2/2;                              % right neighbor
      A(n,n) = -(dy2 + dx2 + h*dx2*dy/k)/2;          % self coupling
      A(n,n-Nx) = dx2/2;  B(n) = -h*dx2*dy*Tinf/k/2; % bottom and top neighbors
%
%   solve for the desired temperature profile
      T = A\B;     
%
%   let's show the 5-striped nature of the matrix using Matlab's spy function
      nfig = nfig+1;  figure(nfig)
      spy(A), axis tight
      title('Block2D\_FD1:  Structure of A matrix')
      xlabel('n value');   ylabel('n value')
%
      nfig = nfig+1;  figure(nfig)
      spy(A), axis([250 290 250 290])
      title('Block2D\_FD1:  Structure of A matrix (magnified)')
      xlabel('n value');   ylabel('n value')
%
%   let's store the temperature profile in a 2-D matrix (for ease in visualization)
      u = zeros(Ny,Nx); 
      for j = 1:Ny
        n1 = (j-1)*Nx+1;   n2 = n1+Nx-1;    u(j,:) = T(n1:n2)';
      end
%       
%   also we should add the boundary points to the final solution for plotting
      Nx = Nx + 1;  Ny = Ny+1;             % increase number of points
      xp = [x Lx/2];   yp = [0 y];         % and boundary points
      [xx,yy] = meshgrid(xp,yp);           % create mesh for 2-D plotting                          
      up = zeros(Ny,Nx);                   % allocate space for up(x,y)
      up(1,:) = Tw*ones(size(xp));         % sets bottom fixed temperatures
      for j = 1:Ny-1                       % adds right fixed temp to up(x,y) array
        up(j+1,:) = [u(j,:) Tw];
      end
%
%   now form actual temperature distribution of interest  --  T(x,y)
      T = up;  % we have used the actual temperature, so no change is necessary here
%
%   plot T(x,y) for a few y values              
      nfig = nfig+1;  figure(nfig)
      plot(xp*100,T(Ny:-8:13,:),'LineWidth',2),grid
      title(['Block2D\_FD1:  T(x) profiles for various y values (Nx = ', ...
             num2str(Nx-1),', Ny = ',num2str(Ny-1),')'])
      xlabel('x location from block centerline (cm)');    
      ylabel('Temperature (C)')                                 
      for j = Ny:-8:13,  gtext(['y = ',num2str(yp(j)*100),' cm']);   end   
%
%   various 3-D plots                           
      nfig = nfig+1;  figure(nfig)
      surf(xx*100,yy*100,T), shading interp, colorbar
      axis([0 Lx*100/2 0 Ly*100 50 90])
      title(['Block2D\_FD1:  T(x,y) in 2-D block (Nx = ', ...
             num2str(Nx-1),', Ny = ',num2str(Ny-1),')'])
      xlabel('x location (cm)');   ylabel('y location (cm)')
      zlabel('Temperature (C)')                                 
%
      nfig = nfig+1;  figure(nfig)
      surf(xx*100,yy*100,T), shading interp, colorbar, hold on
      [cc,hh] = contour3(xp*100,yp*100,T, [80 70 60 55 52.5]);
      clabel(cc), set(hh,'EdgeColor','k')
      axis([0 Lx*100/2 0 Ly*100 50 90])
      title(['Block2D\_FD1:  T(x,y) in 2-D block (Nx = ', ...
             num2str(Nx-1),', Ny = ',num2str(Ny-1),')'])
      xlabel('x location (cm)');   ylabel('y location (cm)')
      zlabel('Temperature (C)'),   hold off                                 
%
      nfig = nfig+1;  figure(nfig)
      surf(xx*100,yy*100,T), shading interp, view(2), colorbar, hold on
      [cc,hh] = contour3(xp*100,yp*100,T, [80 70 60 55 52.5]);
      clabel(cc), set(hh,'EdgeColor','k')
      axis([0 Lx*100/2 0 Ly*100 50 90])
      title(['Block2D\_FD1:  T(x,y) in 2-D block (Nx = ', ...
             num2str(Nx-1),', Ny = ',num2str(Ny-1),')'])
      xlabel('x location (cm)');   ylabel('y location (cm)')
      zlabel('Temperature (C)'),   hold off                                 
%
      nfig = nfig+1;  figure(nfig)
      surf(xx*100,yy*100,T), shading interp, view(2), colorbar, hold on
      [cc,hh] = contour3(xp*100,yp*100,T, [80 70 60 55 52.5]);
      clabel(cc), set(hh,'EdgeColor','k')
      axis equal
      axis([0 Lx*100/2 0 Ly*100 50 90])
      title(['Block2D\_FD1:  T(x,y) in 2-D block (Nx = ', ...
             num2str(Nx-1),', Ny = ',num2str(Ny-1),')'])
      xlabel('x location (cm)');   ylabel('y location (cm)')
      zlabel('Temperature (C)'),   hold off                                 
%
%   now compute the heat flows across each surface of the 2-D block
%   (here Nx and Ny are the total number of temperature points in each direction)
%   Top Surface
      dTdy_Ny = (3*T(Ny,:) - 4*T(Ny-1,:) + T(Ny-2,:))/(2*dy);
      qtop1 = -k*(dTdy_Ny(1)*dx/2 + sum(dTdy_Ny(2:Nx-1))*dx + dTdy_Ny(Nx)*dx/2);
      qtop2 = h*(T(Ny,1)*dx/2 + sum(T(Ny,2:Nx-1))*dx + T(Ny,Nx)*dx/2) - h*Tinf*Lx/2;
%   Right Side
      dTdx_Nx = (3*T(:,Nx) - 4*T(:,Nx-1) + T(:,Nx-2))/(2*dx);
      qrt = -k*(dTdx_Nx(1)*dy/2 + sum(dTdx_Nx(2:Ny-1))*dy + dTdx_Nx(Ny)*dy/2);
%   Bottom Surface
      dTdy_1 = (-3*T(1,:) + 4*T(2,:) - T(3,:))/(2*dy);
      qbot = -k*(dTdy_1(1)*dx/2 + sum(dTdy_1(2:Nx-1))*dx + dTdy_1(Nx)*dx/2);
%   Left Side
      dTdx_1 = (-3*T(:,1) + 4*T(:,2) - T(:,3))/(2*dx);
      qlt = -k*(dTdx_1(1)*dy/2 + sum(dTdx_1(2:Ny-1))*dy + dTdx_1(Ny)*dy/2);
%
      fprintf(1,'\n  BLOCK2D_FD1 -- Heat transfer (W) across boundaries for the right half of a 2-D block \n')
      fprintf(1,'                         (Finite Difference Grid:  Nx = %3i   and   Ny = %3i) \n',Nx-1,Ny-1)
      fprintf(1,'    in top boundary (conduction):   %8.2f  \n',-qtop1)
      fprintf(1,'    in top boundary (convection):   %8.2f  \n',-qtop2)
      fprintf(1,'    out right boundary:             %8.2f  \n',qrt)
      fprintf(1,'    out bottom boundary:            %8.2f  \n',-qbot)
      fprintf(1,'    out left boundary:              %8.2f  \n',-qlt)
%
%   end of problem
