Commits

Matt Knepley  committed a8ecc83

DMPlex: Fixed refinement of hybrid hex cells
- Copied code from the tet case

  • Participants
  • Parent commits 61311cf

Comments (0)

Files changed (1)

File src/dm/impls/plex/plexrefine.c

     for (c = cMax; c < cEnd; ++c) {
       const PetscInt  newp = (cMax - cStart)*8 + (c - cMax)*4;
       const PetscInt *cone, *ornt, *fornt;
-      PetscInt        coneNew[6], orntNew[6];
+      PetscInt        coneNew[6], orntNew[6], o, of, i;
 
       ierr = DMPlexGetCone(dm, c, &cone);CHKERRQ(ierr);
       ierr = DMPlexGetConeOrientation(dm, c, &ornt);CHKERRQ(ierr);
       ierr = DMPlexGetConeOrientation(dm, cone[0], &fornt);CHKERRQ(ierr);
+      o = ornt[0] < 0 ? -1 : 1;
       for (r = 0; r < 4; ++r) {
         PetscInt subfA = GetQuadSubface_Static(ornt[0], r);
         PetscInt edgeA = GetQuadEdge_Static(ornt[0], r);
-        PetscInt edgeB = (edgeA+3)%4;
+        PetscInt edgeB = GetQuadEdge_Static(ornt[0], (r+3)%4);
         if (ornt[0] != ornt[1]) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Inconsistent ordering for matching ends of hybrid cell %d: %d != %d", c, ornt[0], ornt[1]);
         coneNew[0]         = fStartNew + (cone[0] - fStart)*4 + subfA;
         orntNew[0]         = ornt[0];
         coneNew[1]         = fStartNew + (cone[1] - fStart)*4 + subfA;
         orntNew[1]         = ornt[0];
-        coneNew[(r+0)%4+2] = fStartNew + (fMax - fStart)*4 + (cMax - cStart)*12 + (cone[edgeA+2] - fMax)*2 + (fornt[edgeA] < 0 ? 1 : 0);
-        orntNew[(r+0)%4+2] = ornt[edgeA];
-        coneNew[(r+1)%4+2] = fStartNew + (fMax - fStart)*4 + (cMax - cStart)*12 + (fEnd          - fMax)*2 + (c - cMax)*4 + edgeA;
-        orntNew[(r+1)%4+2] = 0;
-        coneNew[(r+2)%4+2] = fStartNew + (fMax - fStart)*4 + (cMax - cStart)*12 + (fEnd          - fMax)*2 + (c - cMax)*4 + edgeB;
-        orntNew[(r+2)%4+2] = -2;
-        coneNew[(r+3)%4+2] = fStartNew + (fMax - fStart)*4 + (cMax - cStart)*12 + (cone[edgeB+2] - fMax)*2 + (fornt[edgeB] < 0 ? 0 : 1);
-        orntNew[(r+3)%4+2] = ornt[edgeB];
+        of = fornt[edgeA] < 0 ? -1 : 1;
+        i  = GetQuadEdgeInverse_Static(ornt[0], r) + 2;
+        coneNew[i] = fStartNew + (fMax - fStart)*4 + (cMax - cStart)*12 + (cone[2+edgeA] - fMax)*2 + (o*of < 0 ? 1 : 0);
+        orntNew[i] = ornt[edgeA];
+        i  = GetQuadEdgeInverse_Static(ornt[0], (r+1)%4) + 2;
+        coneNew[i] = fStartNew + (fMax - fStart)*4 + (cMax - cStart)*12 + (fEnd          - fMax)*2 + (c - cMax)*4 + edgeA;
+        orntNew[i] = 0;
+        i  = GetQuadEdgeInverse_Static(ornt[0], (r+2)%4) + 2;
+        coneNew[i] = fStartNew + (fMax - fStart)*4 + (cMax - cStart)*12 + (fEnd          - fMax)*2 + (c - cMax)*4 + edgeB;
+        orntNew[i] = -2;
+        of = fornt[edgeB] < 0 ? -1 : 1;
+        i  = GetQuadEdgeInverse_Static(ornt[0], (r+3)%4) + 2;
+        coneNew[i] = fStartNew + (fMax - fStart)*4 + (cMax - cStart)*12 + (cone[2+edgeB] - fMax)*2 + (o*of < 0 ? 0 : 1);
+        orntNew[i] = ornt[edgeB];
         ierr       = DMPlexSetCone(rdm, newp+r, coneNew);CHKERRQ(ierr);
         ierr       = DMPlexSetConeOrientation(rdm, newp+r, orntNew);CHKERRQ(ierr);
 #if 1
 #endif
       }
     }
-    /* Interior faces have 4 edges and 2 cells */
+    /* Interior cell faces have 4 edges and 2 cells */
     for (c = cStart; c < cMax; ++c) {
       const PetscInt  newCells[24] = {0, 3,  2, 3,  1, 2,  0, 1,  4, 5,  5, 6,  6, 7,  4, 7,  0, 4,  3, 5,  2, 6,  1, 7};
       const PetscInt *cone, *ornt;
 #endif
         for (s = 0; s < size; ++s) {
           const PetscInt *coneCell, *orntCell, *fornt;
+          PetscInt        o, of;
 
           ierr = DMPlexGetCone(dm, support[s], &coneCell);CHKERRQ(ierr);
           ierr = DMPlexGetConeOrientation(dm, support[s], &orntCell);CHKERRQ(ierr);
+          o = orntCell[0] < 0 ? -1 : 1;
           for (c = 2; c < 6; ++c) if (coneCell[c] == f) break;
           if (c >= 6) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Could not find face %d in cone of cell %d", f, support[s]);
           ierr = DMPlexGetConeOrientation(dm, coneCell[0], &fornt);CHKERRQ(ierr);
-          supportNew[s] = cStartNew + (cMax - cStart)*8 + (support[s] - cMax)*4 + (GetQuadEdgeInverse_Static(orntCell[0], c-2) + (fornt[c-2] < 0 ? 1-r : r))%4;
+          of = fornt[c-2] < 0 ? -1 : 1;
+          supportNew[s] = cStartNew + (cMax - cStart)*8 + (support[s] - cMax)*4 + (GetQuadEdgeInverse_Static(orntCell[0], c-2) + (o*of < 0 ? 1-r : r))%4;
         }
         ierr = DMPlexSetSupport(rdm, newp, supportNew);CHKERRQ(ierr);
 #if 1
       ierr = DMPlexGetCone(dm, c, &cone);CHKERRQ(ierr);
       ierr = DMPlexGetConeOrientation(dm, c, &ornt);CHKERRQ(ierr);
       for (r = 0; r < 4; ++r) {
+#if 0
         coneNew[0] = eStartNew + (eMax - eStart)*2 + (cone[0] - fStart)*4 + GetQuadSubface_Static(ornt[0], r);
         orntNew[0] = 0;
         coneNew[1] = eStartNew + (eMax - eStart)*2 + (cone[1] - fStart)*4 + GetQuadSubface_Static(ornt[1], r);
         orntNew[2] = 0;
         coneNew[3] = eStartNew + (eMax - eStart)*2 + (fMax - fStart)*4 + (cMax - cStart)*6 + (eEnd - eMax) + (fEnd                                   - fMax) + (c - cMax);
         orntNew[3] = 0;
+#else
+        coneNew[0] = eStartNew + (eMax - eStart)*2 + (cone[0] - fStart)*4 + r;
+        orntNew[0] = 0;
+        coneNew[1] = eStartNew + (eMax - eStart)*2 + (cone[1] - fStart)*4 + r;
+        orntNew[1] = 0;
+        coneNew[2] = eStartNew + (eMax - eStart)*2 + (fMax - fStart)*4 + (cMax - cStart)*6 + (eEnd - eMax) + (cone[2+r] - fMax);
+        orntNew[2] = 0;
+        coneNew[3] = eStartNew + (eMax - eStart)*2 + (fMax - fStart)*4 + (cMax - cStart)*6 + (eEnd - eMax) + (fEnd      - fMax) + (c - cMax);
+        orntNew[3] = 0;
+#endif
         ierr = DMPlexSetCone(rdm, newp+r, coneNew);CHKERRQ(ierr);
         ierr = DMPlexSetConeOrientation(rdm, newp+r, orntNew);CHKERRQ(ierr);
 #if 1
           if (support[s] < cMax) {
             supportRef[2+s] = fStartNew + (fMax - fStart)*4 + (support[s] - cStart)*12 + newFaces[c*4 + GetQuadEdgeInverse_Static(orntCell[c], r)];
           } else {
-            supportRef[2+s] = fStartNew + (fMax - fStart)*4 + (cMax - cStart)*12 + (fEnd - fMax)*2 + (support[s] - cMax)*4 + GetQuadEdgeInverse_Static(orntCell[c], r);
+            supportRef[2+s] = fStartNew + (fMax - fStart)*4 + (cMax - cStart)*12 + (fEnd - fMax)*2 + (support[s] - cMax)*4 + r;
           }
         }
         ierr = DMPlexSetSupport(rdm, newp, supportRef);CHKERRQ(ierr);
         ierr = DMPlexGetConeOrientation(dm, support[s], &cornt);CHKERRQ(ierr);
         for (c = 0; c < csize; ++c) if (ccone[c] == f) break;
         if ((c < 2) || (c >= csize)) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Hybrid face %d is not in cone of hybrid cell %d", f, support[s]);
-        supportRef[2+s] = fStartNew + (fMax - fStart)*4 + (cMax - cStart)*12 + (fEnd - fMax)*2 + (support[s] - cMax)*4 + GetQuadSubfaceInverse_Static(cornt[0], c-2);
+        supportRef[2+s] = fStartNew + (fMax - fStart)*4 + (cMax - cStart)*12 + (fEnd - fMax)*2 + (support[s] - cMax)*4 + c-2;
       }
       ierr = DMPlexSetSupport(rdm, newp, supportRef);CHKERRQ(ierr);
 #if 1