message: Permute quadrature points when evaluating velocity coefficient
diff: @@ -182,6 +182,7 @@ void DGTraceIntegrator::SetupPA(const FiniteElementSpace &fes, FaceType type)
  for (int q = 0; q < nq; ++q)
  {
  C(q,f_ind) = rho->Eval(T, ir->IntPoint(q));
  + // TODO
  }
  f_ind++;
  }
  @@ -206,16 +207,40 @@ void DGTraceIntegrator::SetupPA(const FiniteElementSpace &fes, FaceType type)
  int inf1, inf2;
  fes.GetMesh()->GetFaceElements(f, &e1, &e2);
  fes.GetMesh()->GetFaceInfos(f, &inf1, &inf2);
  + int face_id = inf1 / 64;
  if ((type==FaceType::Interior && (e2>=0 || (e2<0 && inf2>=0))) ||
  (type==FaceType::Boundary && e2<0 && inf2<0) )
  {
  ElementTransformation& T = *fes.GetMesh()->GetFaceTransformation(f);
  for (int q = 0; q < nq; ++q)
  {
  + // Convert to lexicographic ordering
  + int iq;
  + if (dim == 2)
  + {
  + iq = (face_id == 2 || face_id == 3) ? nq-1-q : q;
  + }
  + else if (dim == 3)
  + {
  + int q_i = q%quad1D;
  + int q_j = q/quad1D;
  + if (face_id==2 || face_id==1 || face_id==5)
  + {
  + iq = q_i + q_j*quad1D;
  + }
  + else if (face_id==3 || face_id==4)
  + {
  + iq = (quad1D-1-q_i) + q_j*quad1D;
  + }
  + else//face_id==0
  + {
  + iq = q_i + (quad1D-1-q_j)*quad1D;
  + }
  + }
  u->Eval(Vq, T, ir->IntPoint(q));
  for (int i = 0; i < dim; ++i)
  {
  - C(i,q,f_ind) = Vq(i);
  + C(i,iq,f_ind) = Vq(i);
  }
  }
  f_ind++;
  