% This Matlab script solves the two-dimensional Euler
% equations using a finite volume algorithm for the flow
% around a cylinder.
%

% If FVMrestart does not exist, set it to zero (assume not restarting)
if (~exist('FVMrestart')),
  FVMrestart = 0;
end

if (FVMrestart == 0),
  % Not a restart, need to initialize a bunch of data
  cyl_initmesh;

  % After this command, the following variables are present and
  % will be useful in the finite volume implementation.
  %
  % Nt: Number of triangles (i.e. elements) in mesh
  %
  % Nv: Number of nodes (i.e. vertices) in mesh
  %
  % Ne: Number of internal edges (between two triangles)
  %
  % Nbe: 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.
  %
  % edge2tri(1, i) = index of first  node of edge i
  % edge2tri(2, i) = index of second node of edge i
  % edge2tri(3, i) = index of triangle to left of edge i
  % edge2tri(4, i) = index of triangle to right of edge i
  %
  % The edge2tri nodes are oriented so that the triangle in 
  % edge2tri(3,i) is to the left when moving from node 1 to node 2
  %
  % bedge2tri(1, i) = index of first node for bedge i
  % bedge2tri(2, i) = index of second node for bedge i
  % bedge2tri(3, i) = index of the triangle for bedge i
  % bedge2tri(4, i) = boundary ID for bedge i 
  %                   (0 for cylinder, 1 for outer boundary)
  % 
  % The bedge2tri nodes are oriented so that the triangle is to 
  % the left when moving from node 1 to node 2

  % Enter freestream and initial Mach number
  Minf =  2.0;
  Minit = 0.0;
  
  % Specific heat ratio
  gamma = 1.4;
  
  % Initialize solution
  Uinit = [1; Minit; 0; 1/(gamma-1)/gamma + 0.5*Minit^2]; 
  U = ones(4,Nt);
  U(1,:) = Uinit(1)*U(1,:);
  U(2,:) = Uinit(2)*U(2,:);
  U(3,:) = Uinit(3)*U(3,:);
  U(4,:) = Uinit(4)*U(4,:);

  n = 0;

  % Clear figures
  figure(1); clf;
  figure(2); clf;

end

% Set tolerances on iteration and residual tolerance
if (~exist('ntol')),
  ntol = 100;
end
if (~exist('Rtol')),
  Rtol = -5;
end
Rmax = 0;

% Set CFL
if (~exist('CFL')),
  CFL = 0.9;
end

% Calculate edge normals and lengths
enormal = [xy(2,edge2tri(1,:))-xy(2,edge2tri(2,:)); ...
           xy(1,edge2tri(2,:))-xy(1,edge2tri(1,:))];
elength = sqrt(enormal(1,:).^2 + enormal(2,:).^2);
enormal = [enormal(1,:)./elength; enormal(2,:)./elength];

% Calculate bedge normals and lengths
bnormal = [xy(2,bedge2tri(1,:))-xy(2,bedge2tri(2,:)); ...
           xy(1,bedge2tri(2,:))-xy(1,bedge2tri(1,:))];
blength = sqrt(bnormal(1,:).^2 + bnormal(2,:).^2);
bnormal = [bnormal(1,:)./blength; bnormal(2,:)./blength];

% Calculate cell areas
A = zeros(1,Nt);
for i = 1:Nt,
  dx21 = xy(1,tri2nod(2,i)) - xy(1,tri2nod(1,i));
  dy21 = xy(2,tri2nod(2,i)) - xy(2,tri2nod(1,i));
  
  dx31 = xy(1,tri2nod(3,i)) - xy(1,tri2nod(1,i));
  dy31 = xy(2,tri2nod(3,i)) - xy(2,tri2nod(1,i));

  A(i) = 0.5*(dx21*dy31 - dy21*dx31);
end
if (min(A) < 0),
  fprintf('ERROR: Negative area found!\n');
end

% Calculate freestream states for use in setting boundary conditions
Uinf = [1; Minf; 0; 1/(gamma-1)/gamma + 0.5*Minf^2]; 

% Loop until one of the termination criterion are met
while ((n < ntol) & (Rmax > Rtol)),
  
  % Zero residual and dt
  R  = zeros(4,Nt);
  dt = zeros(1,Nt);
  
  % Calculate the flux at each interior edge
  for i = 1:Ne,
    iL = edge2tri(3,i);
    iR = edge2tri(4,i);
    [H, smax] = eulerflux(U(:,iL),U(:,iR),enormal(:,i),gamma);
    R(:,iL) = R(:,iL) - H*elength(i);
    R(:,iR) = R(:,iR) + H*elength(i);
    dt(iL) = dt(iL) + smax*elength(i);
    dt(iR) = dt(iR) + smax*elength(i);
  end
  
  % Calculate the flux at each boundary edge
  for i = 1:Nbe,
    iL = bedge2tri(3,i);
    if (bedge2tri(4,i) == 1), % Outer boundary
      [H, smax] = eulerflux(U(:,iL),Uinf,bnormal(:,i),gamma);
    else, % Cylinder boundary
      [H, smax] = wallflux(U(:,iL),bnormal(:,i),gamma);
    end
    R(:,iL) = R(:,iL) - H*blength(i);
    dt(iL) = dt(iL) + smax*blength(i);
  end
  
  % Local timestep calculation
  dtA = 2*CFL./dt;
  
  % Calculate maximum residual
  Rmax = log10(max(max(abs(R))));

  % Forward Euler step
  U(1,:) = U(1,:) - dtA.*R(1,:);
  U(2,:) = U(2,:) - dtA.*R(2,:);
  U(3,:) = U(3,:) - dtA.*R(3,:);
  U(4,:) = U(4,:) - dtA.*R(4,:);

  % Increment iteration
  n = n + 1;

  % Calculate some useful quantities (Mach number)
  q = sqrt(U(2,:).^2 + U(3,:).^2)./U(1,:);
  p = (gamma-1)*(U(4,:) - 0.5*U(1,:).*q.^2);
  c = sqrt(gamma*p./U(1,:));
  M = q./c;
  
  % Calculate lift and drag coefficients
  CalcForces;
  
  % Plot current solution
  figure(1);
  clf;
  patch('Vertices',xy','Faces',tri2nod','CData',M','FaceColor', ...
        'flat','EdgeColor','none');
  title(sprintf('Minf = %5.2f, CD = %4.2f, Log10 Max Res = %7.3f\n', ...
                Minf,CD,Rmax));
  xlabel('x'); ylabel('y');
  caxis([0,4]);
  colorbar;
  
  % Plot max residual and drag coefficient
  figure(2);
  subplot(211);
  hold on; 
  plot(n,Rmax,'*'); 
  hold off;
  xlabel('Iteration');
  ylabel('Log10 Rmax');

  subplot(212);
  hold on;
  plot(n,CD,'*');
  hold off;
  xlabel('Iteration');
  ylabel('C_D');

  drawnow;

end