message: Use lexicographic vertex ordering for refined tetrahedron
diff: @@ -1165,24 +1165,31 @@ RefinedGeometry * GeometryRefiner::Refine(Geometry::Type Geom,
  Array<int> vi((n+1)*(n+1)*(n+1));
  vi = -1;
  m = 0;
  - for (k = 0; k <= n; k++)
  - for (j = 0; j <= k; j++)
  - for (i = 0; i <= j; i++)
  +
  + // vertices are given in lexicographic ordering on the reference
  + // element
  + for (int kk = 0; kk <= n; kk++)
  + for (int jj = 0; jj <= n-kk; jj++)
  + for (int ii = 0; ii <= n-jj-kk; ii++)
  {
  IntegrationPoint &ip = RG->RefPts.IntPoint(m);
  - // map the coordinates to the reference tetrahedron
  - // (0,0,0) -> (0,0,0)
  - // (0,0,1) -> (1,0,0)
  - // (1,1,1) -> (0,1,0)
  - // (0,1,1) -> (0,0,1)
  - double w = cp[k-j] + cp[i] + cp[j-i] + cp[Times-k];
  - ip.x = cp[k-j]/w;
  - ip.y = cp[i]/w;
  - ip.z = cp[j-i]/w;
  + double w = cp[ii] + cp[jj] + cp[kk] + cp[Times-ii-jj-kk];
  + ip.x = cp[ii]/w;
  + ip.y = cp[jj]/w;
  + ip.z = cp[kk]/w;
  + // (ii,jj,kk) are coordinates in the refernce tetrahedron,
  + // transform to coordinates (i,j,k) in the auxiliary
  + // tetrahedron defined by (0,0,0), (0,0,1), (1,1,1), (0,1,1)
  + int i = jj;
  + int j = jj+kk;
  + int k = ii+jj+kk;
  l = i + (j + k * (n+1)) * (n+1);
  + // map from linear Cartesian hex index in the auxiliary tet
  + // to lexicographic in the reference tet
  vi[l] = m;
  m++;
  }
  +
  if (m != (n+3)*(n+2)*(n+1)/6)
  {
  mfem_error("GeometryRefiner::Refine() for TETRAHEDRON #1");
  