% Calculate error indicator
ErrInd = zeros(Nt,1);

for i = 1:Ne,
  iL = edge2tri(3,i);
  iR = edge2tri(4,i);
  Mdiff = abs(M(iL) - M(iR));
  ErrInd(iL) = ErrInd(iL) + Mdiff*elength(i);
  ErrInd(iR) = ErrInd(iR) + Mdiff*elength(i);
end
  
for i = 1:Nbe,
  iL = bedge2tri(3,i);
  if (bedge2tri(4,i) == 0), % Cylinder surface
    rL = U(1,iL);
    uL = U(2,iL)/rL;
    vL = U(3,iL)/rL;
    unL = uL*bnormal(1,i) + vL*bnormal(2,i);
    qL2 = uL^2 + vL^2;
    pL = (gamma-1)*(U(4,iL) - 0.5*rL*qL2);
    cL = sqrt(gamma*pL/rL);
    MnL = unL/cL;
    Mdiff = abs(MnL);
  else, % Outer boundary
    Mdiff = 0;
  end
  ErrInd(iL) = ErrInd(iL) + Mdiff*blength(i);
end

% Sort cells by error indicator
[ErrIndSort, TSort] = sort(-ErrInd);

% Refine the top eta cells (approximately)
eta = 0.25;
Nrefine = floor(eta*Nt);
RefineList = TSort(1:Nrefine);