Commits

Jed Brown committed 48e10e4 Merge

Merge branch 'knepley/pylith'

* knepley/pylith:
DMPlex: Fixed bug in DMPlexOrient() - Must visit faces in BFS order so that all comparisons are in a connected component of the mesh
DM: Fix leak when creating global section
DMPlex: Fix allocation for embedded submesh
DMPlex: Removed debugging prints
DMPlex: Fix leak in interpolation
DMPlex: Added DMPlexOrient() which correctly orients a set of adjacent cells - You need this when extracting an embedded submesh
DMPlex: Fix memory management in submesh creation

Comments (0)

Files changed (5)

include/petscdmplex.h

 PETSC_EXTERN PetscErrorCode DMPlexGetMaxSizes(DM, PetscInt *, PetscInt *);
 PETSC_EXTERN PetscErrorCode DMPlexSymmetrize(DM);
 PETSC_EXTERN PetscErrorCode DMPlexStratify(DM);
+PETSC_EXTERN PetscErrorCode DMPlexOrient(DM);
 PETSC_EXTERN PetscErrorCode DMPlexGetCoordinateSection(DM, PetscSection *);
 PETSC_EXTERN PetscErrorCode DMPlexSetCoordinateSection(DM, PetscSection);
 PETSC_EXTERN PetscErrorCode DMPlexSetPreallocationCenterDimension(DM,PetscInt);

src/dm/impls/plex/plex.c

 }
 
 #undef __FUNCT__
+#define __FUNCT__ "DMPlexOrient"
+/* Trys to give the mesh a consistent orientation */
+PetscErrorCode DMPlexOrient(DM dm)
+{
+  PetscBT        seenCells, flippedCells, seenFaces;
+  PetscInt      *faceFIFO, fTop, fBottom;
+  PetscInt       dim, h, cStart, cEnd, c, fStart, fEnd, face, maxConeSize, *revcone, *revconeO;
+  PetscErrorCode ierr;
+
+  PetscFunctionBegin;
+  /* Truth Table
+     mismatch    flips   do action   mismatch   flipA ^ flipB   action
+         F       0 flips     no         F             F           F
+         F       1 flip      yes        F             T           T
+         F       2 flips     no         T             F           T
+         T       0 flips     yes        T             T           F
+         T       1 flip      no
+         T       2 flips     yes
+  */
+  ierr = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr);
+  ierr = DMPlexGetVTKCellHeight(dm, &h);CHKERRQ(ierr);
+  ierr = DMPlexGetHeightStratum(dm, h,   &cStart, &cEnd);CHKERRQ(ierr);
+  ierr = DMPlexGetHeightStratum(dm, h+1, &fStart, &fEnd);CHKERRQ(ierr);
+  ierr = PetscBTCreate(cEnd - cStart, &seenCells);CHKERRQ(ierr);
+  ierr = PetscBTMemzero(cEnd - cStart, seenCells);CHKERRQ(ierr);
+  ierr = PetscBTCreate(cEnd - cStart, &flippedCells);CHKERRQ(ierr);
+  ierr = PetscBTMemzero(cEnd - cStart, flippedCells);CHKERRQ(ierr);
+  ierr = PetscBTCreate(fEnd - fStart, &seenFaces);CHKERRQ(ierr);
+  ierr = PetscBTMemzero(fEnd - fStart, seenFaces);CHKERRQ(ierr);
+  ierr = PetscMalloc((fEnd - fStart) * sizeof(PetscInt), &faceFIFO);CHKERRQ(ierr);
+  fTop = fBottom = 0;
+  /* Initialize FIFO with first cell */
+  {
+    const PetscInt *cone;
+    PetscInt        coneSize;
+
+    ierr = DMPlexGetConeSize(dm, cStart, &coneSize);CHKERRQ(ierr);
+    ierr = DMPlexGetCone(dm, cStart, &cone);CHKERRQ(ierr);
+    for (c = 0; c < coneSize; ++c) {
+      faceFIFO[fBottom++] = cone[c];
+      ierr = PetscBTSet(seenFaces, cone[c]-fStart);CHKERRQ(ierr);
+    }
+  }
+  /* Consider each face in FIFO */
+  while (fTop < fBottom) {
+    const PetscInt *support, *coneA, *coneB, *coneOA, *coneOB;
+    PetscInt        supportSize, coneSizeA, coneSizeB, posA = -1, posB = -1;
+    PetscInt        seenA, flippedA, seenB, flippedB, mismatch;
+
+    face = faceFIFO[fTop++];
+    ierr = DMPlexGetSupportSize(dm, face, &supportSize);CHKERRQ(ierr);
+    ierr = DMPlexGetSupport(dm, face, &support);CHKERRQ(ierr);
+    if (supportSize < 2) continue;
+    if (supportSize != 2) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Faces should separate only two cells, not %d", supportSize);
+    seenA    = PetscBTLookup(seenCells,    support[0]-cStart);
+    flippedA = PetscBTLookup(flippedCells, support[0]-cStart);
+    seenB    = PetscBTLookup(seenCells,    support[1]-cStart);
+    flippedB = PetscBTLookup(flippedCells, support[1]-cStart);
+
+    ierr = DMPlexGetConeSize(dm, support[0], &coneSizeA);CHKERRQ(ierr);
+    ierr = DMPlexGetConeSize(dm, support[1], &coneSizeB);CHKERRQ(ierr);
+    ierr = DMPlexGetCone(dm, support[0], &coneA);CHKERRQ(ierr);
+    ierr = DMPlexGetCone(dm, support[1], &coneB);CHKERRQ(ierr);
+    ierr = DMPlexGetConeOrientation(dm, support[0], &coneOA);CHKERRQ(ierr);
+    ierr = DMPlexGetConeOrientation(dm, support[1], &coneOB);CHKERRQ(ierr);
+    for (c = 0; c < coneSizeA; ++c) {
+      if (!PetscBTLookup(seenFaces, coneA[c]-fStart)) {
+        faceFIFO[fBottom++] = coneA[c];
+        ierr = PetscBTSet(seenFaces, coneA[c]-fStart);CHKERRQ(ierr);
+      }
+      if (coneA[c] == face) posA = c;
+      if (fBottom > fEnd-fStart) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Face %d was pushed exceeding capacity %d > %d", coneA[c], fBottom, fEnd-fStart);
+    }
+    if (posA < 0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Face %d could not be located in cell %d", face, support[0]);
+    for (c = 0; c < coneSizeB; ++c) {
+      if (!PetscBTLookup(seenFaces, coneB[c]-fStart)) {
+        faceFIFO[fBottom++] = coneB[c];
+        ierr = PetscBTSet(seenFaces, coneB[c]-fStart);CHKERRQ(ierr);
+      }
+      if (coneB[c] == face) posB = c;
+      if (fBottom > fEnd-fStart) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Face %d was pushed exceeding capacity %d > %d", coneA[c], fBottom, fEnd-fStart);
+    }
+    if (posB < 0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Face %d could not be located in cell %d", face, support[1]);
+
+    if (dim == 1) {
+      mismatch = posA == posB;
+    } else {
+      mismatch = coneOA[posA] == coneOB[posB];
+    }
+
+    if (mismatch ^ (flippedA ^ flippedB)) {
+      if (seenA && seenB) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Previously seen cells %d and %d do not match: Fault mesh is non-orientable", support[0], support[1]);
+      if (!seenA && !flippedA) {
+        ierr = PetscBTSet(flippedCells, support[0]-cStart);CHKERRQ(ierr);
+      } else if (!seenB && !flippedB) {
+        ierr = PetscBTSet(flippedCells, support[1]-cStart);CHKERRQ(ierr);
+      } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Inconsistent mesh orientation: Fault mesh is non-orientable");
+    } else if (flippedA && flippedB) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Attempt to flip already flipped cell: Fault mesh is non-orientable");
+    ierr = PetscBTSet(seenCells, support[0]-cStart);CHKERRQ(ierr);
+    ierr = PetscBTSet(seenCells, support[1]-cStart);CHKERRQ(ierr);
+  }
+
+  ierr = DMPlexGetMaxSizes(dm, &maxConeSize, NULL);CHKERRQ(ierr);
+  ierr = DMGetWorkArray(dm, maxConeSize, PETSC_INT, &revcone);CHKERRQ(ierr);
+  ierr = DMGetWorkArray(dm, maxConeSize, PETSC_INT, &revconeO);CHKERRQ(ierr);
+  for (c = cStart; c < cEnd; ++c) {
+    const PetscInt *cone, *coneO;
+    PetscInt        coneSize, faceSize, cp;
+
+    if (!PetscBTLookup(flippedCells, c-cStart)) continue;
+    ierr = DMPlexGetConeSize(dm, c, &coneSize);CHKERRQ(ierr);
+    ierr = DMPlexGetCone(dm, c, &cone);CHKERRQ(ierr);
+    ierr = DMPlexGetConeOrientation(dm, c, &coneO);CHKERRQ(ierr);
+    for (cp = 0; cp < coneSize; ++cp) {
+      const PetscInt rcp = coneSize-cp-1;
+
+      ierr = DMPlexGetConeSize(dm, cone[rcp], &faceSize);CHKERRQ(ierr);
+      revcone[cp]  = cone[rcp];
+      revconeO[cp] = coneO[rcp] >= 0 ? -(faceSize-coneO[rcp]) : faceSize+coneO[rcp];
+    }
+    ierr = DMPlexSetCone(dm, c, revcone);CHKERRQ(ierr);
+    ierr = DMPlexSetConeOrientation(dm, c, revconeO);CHKERRQ(ierr);
+  }
+  ierr = DMRestoreWorkArray(dm, maxConeSize, PETSC_INT, &revcone);CHKERRQ(ierr);
+  ierr = DMRestoreWorkArray(dm, maxConeSize, PETSC_INT, &revconeO);CHKERRQ(ierr);
+  ierr = PetscBTDestroy(&seenCells);CHKERRQ(ierr);
+  ierr = PetscBTDestroy(&flippedCells);CHKERRQ(ierr);
+  ierr = PetscBTDestroy(&seenFaces);CHKERRQ(ierr);
+  ierr = PetscFree(faceFIFO);CHKERRQ(ierr);
+  PetscFunctionReturn(0);
+}
+
+#undef __FUNCT__
 #define __FUNCT__ "DMPlexGetAdjacencySingleLevel_Internal"
 static PetscErrorCode DMPlexGetAdjacencySingleLevel_Internal(DM dm, PetscInt p, PetscBool useClosure, const PetscInt *tmpClosure, PetscInt *adjSize, PetscInt adj[])
 {

src/dm/impls/plex/plexinterpolate.c

   if (face != pEnd[faceDepth]) SETERRQ2(PetscObjectComm((PetscObject) dm), PETSC_ERR_PLIB, "Invalid number of faces %D should be %D", face-pStart[faceDepth], pEnd[faceDepth]-pStart[faceDepth]);
   ierr = PetscFree2(pStart,pEnd);CHKERRQ(ierr);
   ierr = PetscHashIJKLDestroy(&faceTable);CHKERRQ(ierr);
+  ierr = PetscFree2(pStart,pEnd);CHKERRQ(ierr);
   ierr = DMPlexSymmetrize(idm);CHKERRQ(ierr);
   ierr = DMPlexStratify(idm);CHKERRQ(ierr);
   PetscFunctionReturn(0);

src/dm/impls/plex/plexsubmesh.c

 {
   IS               subvertexIS;
   const PetscInt  *subvertices;
-  PetscInt        *pStart, *pEnd, *pMax;
+  PetscInt        *pStart, *pEnd, *pMax, pSize;
   PetscInt         depth, dim, d, numSubVerticesInitial = 0, v;
   PetscErrorCode   ierr;
 
   *nFV      = 0;
   ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr);
   ierr = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr);
-  ierr = PetscMalloc3(dim+1,PetscInt,&pStart,dim+1,PetscInt,&pEnd,dim+1,PetscInt,&pMax);CHKERRQ(ierr);
+  pSize = PetscMax(depth, dim) + 1;
+  ierr = PetscMalloc3(pSize,PetscInt,&pStart,pSize,PetscInt,&pEnd,pSize,PetscInt,&pMax);CHKERRQ(ierr);
   ierr = DMPlexGetHybridBounds(dm, &pMax[depth], depth>1 ? &pMax[depth-1] : NULL, depth > 2 ? &pMax[1] : NULL, &pMax[0]);CHKERRQ(ierr);
   for (d = 0; d <= depth; ++d) {
     ierr = DMPlexGetDepthStratum(dm, d, &pStart[d], &pEnd[d]);CHKERRQ(ierr);
     PetscInt f;
 
     numFaces = 0;
-    ierr     = DMGetWorkArray(subdm, 1, PETSC_INT, (void**) &faces);CHKERRQ(ierr);
+    ierr     = DMGetWorkArray(subdm, 1, PETSC_INT, (void **) &faces);CHKERRQ(ierr);
     for (f = firstFace; f < *newFacePoint; ++f) {
       PetscInt dof, off, d;
 
     ierr = DMRestoreWorkArray(subdm, 4*numFaceVertices * sizeof(PetscInt), PETSC_INT, &orientedVertices);CHKERRQ(ierr);
     ++(*newFacePoint);
   }
+#if 0
   ierr = DMPlexRestoreJoin(subdm, numFaceVertices, subfaceVertices, &numFaces, &faces);CHKERRQ(ierr);
+#else
+  ierr = DMRestoreWorkArray(subdm, 1, PETSC_INT, (void **) &faces);CHKERRQ(ierr);
+#endif
   PetscFunctionReturn(0);
 }
 

src/dm/interface/dm.c

   if (!dm->defaultGlobalSection) {
     if (!dm->defaultSection || !dm->sf) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "DM must have a default PetscSection and PetscSF in order to create a global PetscSection");
     ierr = PetscSectionCreateGlobalSection(dm->defaultSection, dm->sf, PETSC_FALSE, &dm->defaultGlobalSection);CHKERRQ(ierr);
+    ierr = PetscLayoutDestroy(&dm->map);CHKERRQ(ierr);
     ierr = PetscSectionGetValueLayout(PetscObjectComm((PetscObject)dm),dm->defaultGlobalSection,&dm->map);CHKERRQ(ierr);
   }
   *section = dm->defaultGlobalSection;