% Clear variables

clear all;


% Set gas temperature and wall heat transfer coefficients at
% boundaries of the blade.  Note: Tcool(i) and hcool(i) are the
% values of Tcool and hcool for the ith boundary which are numbered
% as follows:  
%
%   1 = 1st internal cooling passage (from leading edge)
%   2 = 2nd internal cooling passage 
%   3 = 3rd internal cooling passage (including trailing edge slot)

Tcool = [600, 600, 600];    % ( K )
hcool = [1500, 1500, 1500]; % ( W/(K m^2) )

% Set blade material properties
kblade   =   30; % ( W/( K m) )
rhoblade = 8000; % ( kg/m^3 )
cpblade  =  400; % ( J/(kg K) ) 

% Set blade chord
chord = 0.04; % (m)

% Set initial temperature
Tinit = 300; % Initial blade temperature (K)

% Time step and final time
dt = 0.01;
tf = 4;
Ntime = ceil(tf/dt)+1;

% Load in the grid file

fname = input('Enter gridfile name: ','s');
load(fname);

fprintf('Number of nodes      = %i\n',Nv);
fprintf('Number of elements   = %i\n',Nt);

% NOTE:  after loading a gridfile using the load(fname) command,
%        three important grid variables and data arrays exist.  These are:
%
% Nt: Number of triangles (i.e. elements) in mesh
%
% Nv: Number of nodes (i.e. vertices) in mesh
%
% Nbc: Number of edges which lie on a boundary of the computational
%      domain.
% 
% tri2nod(3,Nt):  list of the 3 node numbers which form the current
%                 triangle.  Thus, tri2nod(1,i) is the 1st node of
%                 the i'th triangle, tri2nod(2,i) is the 2nd node
%                 of the i'th triangle, etc.
%
% xy(2,Nv): list of the x and y locations of each node.  Thus,
%           xy(1,i) is the x-location of the i'th node, xy(2,i)
%           is the y-location of the i'th node, etc.
%
% bedge(3,Nbc): For each boundary edge, bedge(1,i) and bedge(2,i) 
%               are the node numbers for the nodes at the end
%               points of the i'th boundary edge.  bedge(3,i) is an
%               integer which identifies which boundary the edge is
%               on. In this solver, the third value has the
%               following meaning:
%
%               bedge(3,i) = 0: edge is on the airfoil
%               bedge(3,i) = 1: edge is on the first cooling passage
%               bedge(3,i) = 2: edge is on the second cooling passage
%               bedge(3,i) = 3: edge is on the third cooling passage
% 

% Re-scale xy location because they currently have unit chord.
xy = chord*xy;

% Find the elements containing points A-D
xyp = [ 0.015, 0.225, 0.475, 0.800; ...
        0.015, 0.050, 0.050, 0.020 ]*chord;

iel(1) = findloc(xyp(1,1),xyp(2,1),xy,tri2nod);
iel(2) = findloc(xyp(1,2),xyp(2,2),xy,tri2nod);
iel(3) = findloc(xyp(1,3),xyp(2,3),xy,tri2nod);
iel(4) = findloc(xyp(1,4),xyp(2,4),xy,tri2nod);

% Zero stiffness matrix

K = sparse(Nv, Nv);
M = sparse(Nv, Nv);
b = zeros(Nv, 1);

% Loop over elements and calculate residual and stiffness matrix

Atmp = 0;
dsmin = 0;
for ii = 1:Nt,
  
  kn(1) = tri2nod(1,ii);
  kn(2) = tri2nod(2,ii);
  kn(3) = tri2nod(3,ii);
    
  xe(1) = xy(1,kn(1));
  xe(2) = xy(1,kn(2));
  xe(3) = xy(1,kn(3));

  ye(1) = xy(2,kn(1));
  ye(2) = xy(2,kn(2));
  ye(3) = xy(2,kn(3));

  % Find minimum edge length
  ds21 = 1/sqrt((xe(2)-xe(1))^2 + (ye(2)-ye(1))^2);
  ds32 = 1/sqrt((xe(3)-xe(2))^2 + (ye(3)-ye(2))^2);
  ds13 = 1/sqrt((xe(1)-xe(3))^2 + (ye(1)-ye(3))^2);
  
  dsmin = max([dsmin, ds21, ds32, ds13]);
  
  % Calculate all of the necessary shape function derivatives, the
  % Jacobian of the element, etc.
  
  % Derivatives of node 1's interpolant 
  dNdxi(1,1) = -1.0; % with respect to xi1
  dNdxi(1,2) = -1.0; % with respect to xi2
  
  % Derivatives of node 2's interpolant
  dNdxi(2,1) =  1.0; % with respect to xi1
  dNdxi(2,2) =  0.0; % with respect to xi2

  % Derivatives of node 3's interpolant
  dNdxi(3,1) =  0.0; % with respect to xi1
  dNdxi(3,2) =  1.0; % with respect to xi2
  
  % Sum these to find dxdxi (note: these are constant within an element)
  dxdxi = zeros(2,2);
  for nn = 1:3,
    dxdxi(1,:) = dxdxi(1,:) + xe(nn)*dNdxi(nn,:);
    dxdxi(2,:) = dxdxi(2,:) + ye(nn)*dNdxi(nn,:);
  end
  
  % Calculate determinant for area weighting
  J = abs(dxdxi(1,1)*dxdxi(2,2) - dxdxi(1,2)*dxdxi(2,1));
  A = 0.5*J; % Area is half of the Jacobian
  
  Atmp = Atmp + A;
  
  % Invert dxdxi to find dxidx using inversion rule for a 2x2 matrix
  dxidx = [ dxdxi(2,2)/J, -dxdxi(1,2)/J; ...
	   -dxdxi(2,1)/J,  dxdxi(1,1)/J];
  
  % Calculate dNdx 
  dNdx = dNdxi*dxidx;

  % Add contributions to stiffness matrix for node 1 weighted residual
  K(kn(1), kn(1)) = K(kn(1), kn(1)) - kblade*(dNdx(1,1)*dNdx(1,1) + dNdx(1,2)*dNdx(1,2))*A;
  K(kn(1), kn(2)) = K(kn(1), kn(2)) - kblade*(dNdx(1,1)*dNdx(2,1) + dNdx(1,2)*dNdx(2,2))*A;
  K(kn(1), kn(3)) = K(kn(1), kn(3)) - kblade*(dNdx(1,1)*dNdx(3,1) + dNdx(1,2)*dNdx(3,2))*A;
  
  % Add contributions to stiffness matrix for node 2 weighted residual
  K(kn(2), kn(1)) = K(kn(2), kn(1)) - kblade*(dNdx(2,1)*dNdx(1,1) + dNdx(2,2)*dNdx(1,2))*A;
  K(kn(2), kn(2)) = K(kn(2), kn(2)) - kblade*(dNdx(2,1)*dNdx(2,1) + dNdx(2,2)*dNdx(2,2))*A;
  K(kn(2), kn(3)) = K(kn(2), kn(3)) - kblade*(dNdx(2,1)*dNdx(3,1) + dNdx(2,2)*dNdx(3,2))*A;
  
  % Add contributions to stiffness matrix for node 3 weighted residual
  K(kn(3), kn(1)) = K(kn(3), kn(1)) - kblade*(dNdx(3,1)*dNdx(1,1) + dNdx(3,2)*dNdx(1,2))*A;
  K(kn(3), kn(2)) = K(kn(3), kn(2)) - kblade*(dNdx(3,1)*dNdx(2,1) + dNdx(3,2)*dNdx(2,2))*A;
  K(kn(3), kn(3)) = K(kn(3), kn(3)) - kblade*(dNdx(3,1)*dNdx(3,1) + dNdx(3,2)*dNdx(3,2))*A;

  % Add contributions to mass matrix for node 1 weighted residual
  M(kn(1), kn(1)) = M(kn(1), kn(1)) + (1/12)*rhoblade*cpblade*J;
  M(kn(1), kn(2)) = M(kn(1), kn(2)) + (1/24)*rhoblade*cpblade*J;
  M(kn(1), kn(3)) = M(kn(1), kn(3)) + (1/24)*rhoblade*cpblade*J;
  
  % Add contributions to mass matrix for node 2 weighted residual
  M(kn(2), kn(1)) = M(kn(2), kn(1)) + (1/24)*rhoblade*cpblade*J; 
  M(kn(2), kn(2)) = M(kn(2), kn(2)) + (1/12)*rhoblade*cpblade*J; 
  M(kn(2), kn(3)) = M(kn(2), kn(3)) + (1/24)*rhoblade*cpblade*J; 
  
  % Add contributions to mass matrix for node 3 weighted residual
  M(kn(3), kn(1)) = M(kn(3), kn(1)) + (1/24)*rhoblade*cpblade*J; 
  M(kn(3), kn(2)) = M(kn(3), kn(2)) + (1/24)*rhoblade*cpblade*J; 
  M(kn(3), kn(3)) = M(kn(3), kn(3)) + (1/12)*rhoblade*cpblade*J; 
  
end

dsmin = 1/dsmin;
fprintf('Minimum edge length = %e\n',dsmin);

% Loop over boundary edges and account for bc's
% Note: the bc's are all convective heat transfer coefficient bc's
% so the are of 'Robin' form.  This requires modification of the
% stiffness matrix as well as impacting the right-hand side, b.
%

scool = 0;
shot = 0;
for ii = 1:Nbc,
  
  % Get node numbers on edge
  kn(1) = bedge(1,ii);
  kn(2) = bedge(2,ii);
  
  % Get node coordinates
  xe(1) = xy(1,kn(1));
  xe(2) = xy(1,kn(2));
  
  ye(1) = xy(2,kn(1));
  ye(2) = xy(2,kn(2));
  
  % Calculate edge length
  ds = sqrt((xe(1)-xe(2))^2 + (ye(1)-ye(2))^2);
  
  % Determine the boundary number
  bnum = bedge(3,ii);
  if (bnum == 0), % Blade surface
    
    [Text, hext] = Thgas( 0.5*(xe(1)+xe(2)), 0.5*(ye(1)+ye(2)), ...
                          chord );
    shot = shot + ds;
    
  else, % cooling passage
    Text = Tcool(bnum);
    hext = hcool(bnum);
    scool = scool + ds;
    
  end
  
  % Based on boundary number, set heat transfer bc
  K(kn(1), kn(1)) = K(kn(1), kn(1)) - hext*ds*(1/3);
  K(kn(1), kn(2)) = K(kn(1), kn(2)) - hext*ds*(1/6);
  b(kn(1))        = b(kn(1))        + hext*ds*0.5*Text;
  
  K(kn(2), kn(1)) = K(kn(2), kn(1)) - hext*ds*(1/6);
  K(kn(2), kn(2)) = K(kn(2), kn(2)) - hext*ds*(1/3);
  b(kn(2))        = b(kn(2))        + hext*ds*0.5*Text;
  
end


% Allocate storage for temperature histories
Tmin_hist = zeros(Ntime, 1);
Tmax_hist = zeros(Ntime, 1);
TA_hist   = zeros(Ntime, 1);
TB_hist   = zeros(Ntime, 1);
TC_hist   = zeros(Ntime, 1);
TD_hist   = zeros(Ntime, 1);
t_hist    = zeros(Ntime, 1);


% Set time
t = 0;
T = Tinit*ones(Nv, 1);
n = 1;

% Report outputs
Tmax = max(T);
Tmin = min(T);
TA = findval(xyp(1,1),xyp(2,1),iel(1),xy,tri2nod,T);
TB = findval(xyp(1,2),xyp(2,2),iel(2),xy,tri2nod,T);
TC = findval(xyp(1,3),xyp(2,3),iel(3),xy,tri2nod,T);
TD = findval(xyp(1,4),xyp(2,4),iel(4),xy,tri2nod,T);

fprintf(['t = %e, Tmin = %7.2f, Tmax = %7.2f, TA = %7.2f, TB = ', ...
         '%7.2f, TC = %7.2f, TD = %7.2f\n'],t,Tmin,Tmax,TA,TB,TC,TD);

Tmax_hist(n) = Tmax;
Tmin_hist(n) = Tmin;
TA_hist(n)   = TA;
TB_hist(n)   = TB;
TC_hist(n)   = TC;
TD_hist(n)   = TD;
t_hist(n)    = t;


% First step is backward Euler

% Calculate Kp
Kp = (1/dt)*M - K;

t = t + dt;
n = n + 1;
Tn = T;

% Set bp
bp = (1/dt)*M*Tn + b;

% Solve for temperature at n+1
T = Kp\bp;

% Plot solution
%figure(1);
%bladeplot; drawnow;

% Report outputs
Tmax = max(T);
Tmin = min(T);
TA = findval(xyp(1,1),xyp(2,1),iel(1),xy,tri2nod,T);
TB = findval(xyp(1,2),xyp(2,2),iel(2),xy,tri2nod,T);
TC = findval(xyp(1,3),xyp(2,3),iel(3),xy,tri2nod,T);
TD = findval(xyp(1,4),xyp(2,4),iel(4),xy,tri2nod,T);

fprintf(['t = %e, Tmin = %7.2f, Tmax = %7.2f, TA = %7.2f, TB = ', ...
         '%7.2f, TC = %7.2f, TD = %7.2f\n'],t,Tmin,Tmax,TA,TB,TC,TD);

Tmax_hist(n) = Tmax;
Tmin_hist(n) = Tmin;
TA_hist(n)   = TA;
TB_hist(n)   = TB;
TC_hist(n)   = TC;
TD_hist(n)   = TD;
t_hist(n)    = t;

% Time march using BDF2

% Calculate linear system for BDF2
Kp = (1/dt)*M - (2/3)*K;

while (n < Ntime),
  
  t = t + dt;
  n = n + 1;
  Tn1 = Tn;
  Tn  = T;
  
  % Set bp
  bp = (1/dt)*M*((4/3)*Tn - (1/3)*Tn1) + (2/3)*b;  

  % Solve for temperature at n+1
  T = Kp\bp;

  % Plot solution
%  bladeplot; drawnow;

  % Report outputs
  Tmax = max(T);
  Tmin = min(T);
  TA = findval(xyp(1,1),xyp(2,1),iel(1),xy,tri2nod,T);
  TB = findval(xyp(1,2),xyp(2,2),iel(2),xy,tri2nod,T);
  TC = findval(xyp(1,3),xyp(2,3),iel(3),xy,tri2nod,T);
  TD = findval(xyp(1,4),xyp(2,4),iel(4),xy,tri2nod,T);
  
  fprintf(['t = %e, Tmin = %7.2f, Tmax = %7.2f, TA = %7.2f, TB = ', ...
  '%7.2f, TC = %7.2f, TD = %7.2f\n'],t,Tmin,Tmax,TA,TB,TC,TD);
  
  Tmax_hist(n) = Tmax;
  Tmin_hist(n) = Tmin;
  TA_hist(n)   = TA;
  TB_hist(n)   = TB;
  TC_hist(n)   = TC;
  TD_hist(n)   = TD;
  t_hist(n)    = t;

end

% Plot time histories

subplot(321);
plot(t_hist,Tmin_hist,'-'); hold on;
xlabel('t (sec)');
ylabel('T_{min}');

subplot(322);
plot(t_hist,Tmax_hist,'-'); hold on;
xlabel('t (sec)');
ylabel('T_{max}');

subplot(323);
plot(t_hist,TA_hist,'-'); hold on;
xlabel('t (sec)');
ylabel('T_A');

subplot(324);
plot(t_hist,TB_hist,'-'); hold on;
xlabel('t (sec)');
ylabel('T_B');

subplot(325);
plot(t_hist,TC_hist,'-'); hold on;
xlabel('t (sec)');
ylabel('T_C');

subplot(326);
plot(t_hist,TD_hist,'-'); hold on;
xlabel('t (sec)');
ylabel('T_D');