Commits

Matt Knepley committed 51dc34e Merge

Merge branch 'knepley/fix-plex-ghost-cells'

* knepley/fix-plex-ghost-cells:
DMLabel: Fixed index for ranges with nonzero lower limit
DMPlex: Fixed DMPlexConstructGhostCells() to filter out non-faces from the label
DMPlex: Added DMLabelFilter()

Comments (0)

Files changed (3)

include/petscdmplex.h

 PETSC_EXTERN PetscErrorCode DMLabelCreateIndex(DMLabel, PetscInt, PetscInt);
 PETSC_EXTERN PetscErrorCode DMLabelDestroyIndex(DMLabel);
 PETSC_EXTERN PetscErrorCode DMLabelHasPoint(DMLabel, PetscInt, PetscBool *);
+PETSC_EXTERN PetscErrorCode DMLabelFilter(DMLabel, PetscInt, PetscInt);
 
 PETSC_EXTERN PetscErrorCode DMPlexCreateLabel(DM, const char []);
 PETSC_EXTERN PetscErrorCode DMPlexGetLabelValue(DM, const char[], PetscInt, PetscInt *);

src/dm/impls/plex/plexlabel.c

       const PetscInt point = label->points[label->stratumOffsets[v]+i];
 
       if ((point < pStart) || (point >= pEnd)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Label point %d is not in [%d, %d)", point, pStart, pEnd);
-      ierr = PetscBTSet(label->bt, point);CHKERRQ(ierr);
+      ierr = PetscBTSet(label->bt, point - pStart);CHKERRQ(ierr);
     }
   }
   PetscFunctionReturn(0);
     ++label->stratumSizes[v];
     if (label->bt) {
       if ((point < label->pStart) || (point >= label->pEnd)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Label point %d is not in [%d, %d)", point, label->pStart, label->pEnd);
-      ierr = PetscBTSet(label->bt, point);CHKERRQ(ierr);
+      ierr = PetscBTSet(label->bt, point - label->pStart);CHKERRQ(ierr);
     }
   }
   PetscFunctionReturn(0);
       --label->stratumSizes[v];
       if (label->bt) {
         if ((point < label->pStart) || (point >= label->pEnd)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Label point %d is not in [%d, %d)", point, label->pStart, label->pEnd);
-        ierr = PetscBTClear(label->bt, point);CHKERRQ(ierr);
+        ierr = PetscBTClear(label->bt, point - label->pStart);CHKERRQ(ierr);
       }
       break;
     }
       const PetscInt point = label->points[label->stratumOffsets[v]+i];
 
       if ((point < label->pStart) || (point >= label->pEnd)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Label point %d is not in [%d, %d)", point, label->pStart, label->pEnd);
-      ierr = PetscBTClear(label->bt, point);CHKERRQ(ierr);
+      ierr = PetscBTClear(label->bt, point - label->pStart);CHKERRQ(ierr);
     }
   }
   label->stratumSizes[v] = 0;
   PetscFunctionReturn(0);
 }
 
+#undef __FUNCT__
+#define __FUNCT__ "DMLabelFilter"
+PetscErrorCode DMLabelFilter(DMLabel label, PetscInt start, PetscInt end)
+{
+  PetscInt       v;
+  PetscErrorCode ierr;
+
+  PetscFunctionBegin;
+  label->pStart = start;
+  label->pEnd   = end;
+  if (label->bt) {ierr = PetscBTDestroy(&label->bt);CHKERRQ(ierr);}
+  /* Could squish offsets, but would only make sense if I reallocate the storage */
+  for (v = 0; v < label->numStrata; ++v) {
+    const PetscInt offset = label->stratumOffsets[v];
+    const PetscInt size   = label->stratumSizes[v];
+    PetscInt       off    = offset, q;
+
+    for (q = offset; q < offset+size; ++q) {
+      const PetscInt point = label->points[q];
+
+      if ((point < start) || (point >= end)) continue;
+      label->points[off++] = point;
+    }
+    label->stratumSizes[v] = off-offset;
+  }
+  ierr = DMLabelCreateIndex(label, start, end);CHKERRQ(ierr);
+  PetscFunctionReturn(0);
+}
+
 
 
 #undef __FUNCT__

src/dm/impls/plex/plexsubmesh.c

   IS              valueIS;
   const PetscInt *values;
   PetscInt       *depthShift;
-  PetscInt        depth = 0, numFS, fs, ghostCell, cEnd, c;
+  PetscInt        depth = 0, numFS, fs, fStart, fEnd, ghostCell, cEnd, c;
   PetscErrorCode  ierr;
 
   PetscFunctionBegin;
+  ierr = DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);CHKERRQ(ierr);
   /* Count ghost cells */
   ierr = DMLabelGetValueIS(label, &valueIS);CHKERRQ(ierr);
   ierr = ISGetLocalSize(valueIS, &numFS);CHKERRQ(ierr);
   ierr = ISGetIndices(valueIS, &values);CHKERRQ(ierr);
-
   *numGhostCells = 0;
   for (fs = 0; fs < numFS; ++fs) {
-    PetscInt numBdFaces;
-
-    ierr = DMLabelGetStratumSize(label, values[fs], &numBdFaces);CHKERRQ(ierr);
+    IS              faceIS;
+    const PetscInt *faces;
+    PetscInt        numFaces, f, numBdFaces = 0;
 
+    ierr = DMLabelGetStratumIS(label, values[fs], &faceIS);CHKERRQ(ierr);
+    ierr = ISGetLocalSize(faceIS, &numFaces);CHKERRQ(ierr);
+    ierr = ISGetIndices(faceIS, &faces);CHKERRQ(ierr);
+    for (f = 0; f < numFaces; ++f) {
+      if ((faces[f] >= fStart) && (faces[f] < fEnd)) ++numBdFaces;
+    }
     *numGhostCells += numBdFaces;
+    ierr = ISDestroy(&faceIS);CHKERRQ(ierr);
   }
   ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr);
   ierr = PetscMalloc((depth+1) * sizeof(PetscInt), &depthShift);CHKERRQ(ierr);
     for (f = 0; f < numFaces; ++f) {
       PetscInt size;
 
+      if ((faces[f] < fStart) || (faces[f] >= fEnd)) continue;
       ierr = DMPlexGetSupportSize(dm, faces[f], &size);CHKERRQ(ierr);
       if (size != 1) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "DM has boundary face %d with %d support cells", faces[f], size);
       ierr = DMPlexSetSupportSize(gdm, faces[f] + *numGhostCells, 2);CHKERRQ(ierr);
     ierr = DMLabelGetStratumIS(label, values[fs], &faceIS);CHKERRQ(ierr);
     ierr = ISGetLocalSize(faceIS, &numFaces);CHKERRQ(ierr);
     ierr = ISGetIndices(faceIS, &faces);CHKERRQ(ierr);
-    for (f = 0; f < numFaces; ++f, ++ghostCell) {
+    for (f = 0; f < numFaces; ++f) {
       PetscInt newFace = faces[f] + *numGhostCells;
 
+      if ((faces[f] < fStart) || (faces[f] >= fEnd)) continue;
       ierr = DMPlexSetCone(gdm, ghostCell, &newFace);CHKERRQ(ierr);
       ierr = DMPlexInsertSupport(gdm, newFace, 1, ghostCell);CHKERRQ(ierr);
+      ++ghostCell;
     }
     ierr = ISRestoreIndices(faceIS, &faces);CHKERRQ(ierr);
     ierr = ISDestroy(&faceIS);CHKERRQ(ierr);