Commits

Jed Brown committed b30c3ef Merge

Merge branch 'knepley/pylith'

* knepley/pylith:
DMPlex: Add missing declaration - This should merge alright with the other branch which fixes it
DMPlex: Submesh fixes - Coordinates must use a field - Allow same subpointMap to be specified again
DMPlex: Fix submesh marking for uninterpolated meshes - This now works for embedded submeshes
DMPlex: Interpolation corner cases - Do not interpolate a 1D mesh - No need to copy coordinates of the same DM
DMPlex: Interpolation of a submesh now starts at the faces - Do not know a better way to do this right now
DMPlex: Fix error in work space allocation
ISieve: Coordinate conversion should add a field

  • Participants
  • Parent commits bc22db3, 930319c

Comments (0)

Files changed (5)

File include/petscdmplex.h

 PETSC_EXTERN PetscErrorCode DMPlexGetRefinementUniform(DM, PetscBool *);
 PETSC_EXTERN PetscErrorCode DMPlexSetRefinementUniform(DM, PetscBool);
 PETSC_EXTERN PetscErrorCode DMPlexInterpolate(DM, DM *);
+PETSC_EXTERN PetscErrorCode DMPlexCopyCoordinates(DM, DM);
 PETSC_EXTERN PetscErrorCode DMPlexDistribute(DM, const char[], PetscInt, DM*);
 PETSC_EXTERN PetscErrorCode DMPlexLoad(PetscViewer, DM);
 PETSC_EXTERN PetscErrorCode DMPlexGetSubpointMap(DM, DMLabel*);

File include/sieve/ISieve.hh

       PetscInt                            n;
       PetscErrorCode                      ierr;
 
-      if (!chart.size()) return;
+      ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRXX(ierr);
+      if (!chart.size()) {
+        ierr = PetscSectionSetFieldComponents(coordSection, 0, 1);CHKERRXX(ierr);
+        return;
+      }
       for(typename Section::chart_type::const_iterator p_iter = chart.begin(); p_iter != chart.end(); ++p_iter) {
         min = std::min(min, renumbering[*p_iter]);
         max = std::max(max, renumbering[*p_iter]);
       }
+      ierr = PetscSectionSetFieldComponents(coordSection, 0, coordinates.getFiberDimension(*chart.begin()));CHKERRXX(ierr);
       ierr = PetscSectionSetChart(coordSection, min, max+1);CHKERRXX(ierr);
       for(typename Section::chart_type::const_iterator p_iter = chart.begin(); p_iter != chart.end(); ++p_iter) {
         ierr = PetscSectionSetDof(coordSection, renumbering[*p_iter], coordinates.getFiberDimension(*p_iter));CHKERRXX(ierr);
+        ierr = PetscSectionSetFieldDof(coordSection, renumbering[*p_iter], 0, coordinates.getFiberDimension(*p_iter));CHKERRXX(ierr);
       }
       ierr = PetscSectionSetUp(coordSection);CHKERRXX(ierr);
       ierr = PetscSectionGetStorageSize(coordSection, &n);CHKERRXX(ierr);

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

   PetscFunctionBegin;
   PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
   ierr    = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr);
-  maxSize = 2*PetscMax(PetscMax(PetscPowInt(mesh->maxConeSize,depth+1),PetscPowInt(mesh->maxSupportSize,depth+1)),depth);
+  maxSize = 2*PetscMax(PetscMax(PetscPowInt(mesh->maxConeSize,depth+1),PetscPowInt(mesh->maxSupportSize,depth+1)),depth+1);
   ierr    = DMGetWorkArray(dm, maxSize, PETSC_INT, &fifo);CHKERRQ(ierr);
   if (*points) {
     closure = *points;

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

 /* This interpolates faces for cells at some stratum */
 static PetscErrorCode DMPlexInterpolateFaces_Internal(DM dm, PetscInt cellDepth, DM idm)
 {
+  DMLabel        subpointMap;
   PetscHashIJKL  faceTable;
   PetscInt      *pStart, *pEnd;
   PetscInt       cellDim, depth, faceDepth = cellDepth, numPoints = 0, faceSizeAll = 0, face, c, d;
 
   PetscFunctionBegin;
   ierr = DMPlexGetDimension(dm, &cellDim);CHKERRQ(ierr);
+  /* HACK: I need a better way to determine face dimension, or an alternative to GetFaces() */
+  ierr = DMPlexGetSubpointMap(dm, &subpointMap);CHKERRQ(ierr);
+  if (subpointMap) ++cellDim;
   ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr);
   ++depth;
   ++cellDepth;
 
   PetscFunctionBegin;
   ierr = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr);
+  if (dim <= 1) {
+    ierr = PetscObjectReference((PetscObject) dm);CHKERRQ(ierr);
+    idm  = dm;
+  }
   for (d = 1; d < dim; ++d) {
     /* Create interpolated mesh */
     ierr = DMCreate(PetscObjectComm((PetscObject)dm), &idm);CHKERRQ(ierr);
   PetscErrorCode ierr;
 
   PetscFunctionBegin;
+  if (dmA == dmB) PetscFunctionReturn(0);
   ierr = DMPlexGetDepthStratum(dmA, 0, &vStartA, &vEndA);CHKERRQ(ierr);
   ierr = DMPlexGetDepthStratum(dmB, 0, &vStartB, &vEndB);CHKERRQ(ierr);
   if ((vEndA-vStartA) != (vEndB-vStartB)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "The number of vertices in first DM %d != %d in the second DM", vEndA-vStartA, vEndB-vStartB);

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

       if (!(*nFV)) {ierr = DMPlexGetNumFaceVertices(dm, dim, numCorners, nFV);CHKERRQ(ierr);}
       if (faceSize > *nFV) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Invalid submesh: Too many vertices %d of an element on the surface", faceSize);
       if (faceSize == *nFV) {
+        const PetscInt *cells = NULL;
+        PetscInt        numCells, nc;
+
         ++(*numFaces);
         for (cl = 0; cl < faceSize; ++cl) {
           ierr = DMLabelSetValue(subpointMap, closure[cl], 0);CHKERRQ(ierr);
         }
-        ierr = DMLabelSetValue(subpointMap, cell, 2);CHKERRQ(ierr);
+        ierr = DMPlexGetJoin(dm, faceSize, closure, &numCells, &cells);CHKERRQ(ierr);
+        for (nc = 0; nc < numCells; ++nc) {
+          ierr = DMLabelSetValue(subpointMap, cells[nc], 2);CHKERRQ(ierr);
+        }
+        ierr = DMPlexRestoreJoin(dm, faceSize, closure, &numCells, &cells);CHKERRQ(ierr);
       }
       ierr = DMPlexRestoreTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr);
     }
     PetscSection coordSection, subCoordSection;
     Vec          coordinates, subCoordinates;
     PetscScalar *coords, *subCoords;
-    PetscInt     coordSize, v;
+    PetscInt     numComp, coordSize, v;
 
     ierr = DMPlexGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr);
     ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr);
     ierr = DMPlexGetCoordinateSection(subdm, &subCoordSection);CHKERRQ(ierr);
+    ierr = PetscSectionSetNumFields(subCoordSection, 1);CHKERRQ(ierr);
+    ierr = PetscSectionGetFieldComponents(coordSection, 0, &numComp);CHKERRQ(ierr);
+    ierr = PetscSectionSetFieldComponents(subCoordSection, 0, numComp);CHKERRQ(ierr);
     ierr = PetscSectionSetChart(subCoordSection, firstSubVertex, firstSubVertex+numSubVertices);CHKERRQ(ierr);
     for (v = 0; v < numSubVertices; ++v) {
       const PetscInt vertex    = subVertices[v];
 
       ierr = PetscSectionGetDof(coordSection, vertex, &dof);CHKERRQ(ierr);
       ierr = PetscSectionSetDof(subCoordSection, subvertex, dof);CHKERRQ(ierr);
+      ierr = PetscSectionSetFieldDof(subCoordSection, subvertex, 0, dof);CHKERRQ(ierr);
     }
     ierr = PetscSectionSetUp(subCoordSection);CHKERRQ(ierr);
     ierr = PetscSectionGetStorageSize(subCoordSection, &coordSize);CHKERRQ(ierr);
     PetscSection coordSection, subCoordSection;
     Vec          coordinates, subCoordinates;
     PetscScalar *coords, *subCoords;
-    PetscInt     coordSize;
+    PetscInt     numComp, coordSize;
 
     ierr = DMPlexGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr);
     ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr);
     ierr = DMPlexGetCoordinateSection(subdm, &subCoordSection);CHKERRQ(ierr);
+    ierr = PetscSectionSetNumFields(subCoordSection, 1);CHKERRQ(ierr);
+    ierr = PetscSectionGetFieldComponents(coordSection, 0, &numComp);CHKERRQ(ierr);
+    ierr = PetscSectionSetFieldComponents(subCoordSection, 0, numComp);CHKERRQ(ierr);
     ierr = PetscSectionSetChart(subCoordSection, firstSubPoint[0], firstSubPoint[0]+numSubPoints[0]);CHKERRQ(ierr);
     for (v = 0; v < numSubPoints[0]; ++v) {
       const PetscInt vertex    = subpoints[0][v];
 
       ierr = PetscSectionGetDof(coordSection, vertex, &dof);CHKERRQ(ierr);
       ierr = PetscSectionSetDof(subCoordSection, subvertex, dof);CHKERRQ(ierr);
+      ierr = PetscSectionSetFieldDof(subCoordSection, subvertex, 0, dof);CHKERRQ(ierr);
     }
     ierr = PetscSectionSetUp(subCoordSection);CHKERRQ(ierr);
     ierr = PetscSectionGetStorageSize(subCoordSection, &coordSize);CHKERRQ(ierr);
     PetscSection coordSection, subCoordSection;
     Vec          coordinates, subCoordinates;
     PetscScalar *coords, *subCoords;
-    PetscInt     coordSize, v;
+    PetscInt     numComp, coordSize, v;
 
     ierr = DMPlexGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr);
     ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr);
     ierr = DMPlexGetCoordinateSection(subdm, &subCoordSection);CHKERRQ(ierr);
+    ierr = PetscSectionSetNumFields(subCoordSection, 1);CHKERRQ(ierr);
+    ierr = PetscSectionGetFieldComponents(coordSection, 0, &numComp);CHKERRQ(ierr);
+    ierr = PetscSectionSetFieldComponents(subCoordSection, 0, numComp);CHKERRQ(ierr);
     ierr = PetscSectionSetChart(subCoordSection, firstSubVertex, firstSubVertex+numSubVertices);CHKERRQ(ierr);
     for (v = 0; v < numSubVertices; ++v) {
       const PetscInt vertex    = subVertices[v];
 
       ierr = PetscSectionGetDof(coordSection, vertex, &dof);CHKERRQ(ierr);
       ierr = PetscSectionSetDof(subCoordSection, subvertex, dof);CHKERRQ(ierr);
+      ierr = PetscSectionSetFieldDof(subCoordSection, subvertex, 0, dof);CHKERRQ(ierr);
     }
     ierr = PetscSectionSetUp(subCoordSection);CHKERRQ(ierr);
     ierr = PetscSectionGetStorageSize(subCoordSection, &coordSize);CHKERRQ(ierr);
     PetscSection coordSection, subCoordSection;
     Vec          coordinates, subCoordinates;
     PetscScalar *coords, *subCoords;
-    PetscInt     coordSize;
+    PetscInt     numComp, coordSize;
 
     ierr = DMPlexGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr);
     ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr);
     ierr = DMPlexGetCoordinateSection(subdm, &subCoordSection);CHKERRQ(ierr);
+    ierr = PetscSectionSetNumFields(subCoordSection, 1);CHKERRQ(ierr);
+    ierr = PetscSectionGetFieldComponents(coordSection, 0, &numComp);CHKERRQ(ierr);
+    ierr = PetscSectionSetFieldComponents(subCoordSection, 0, numComp);CHKERRQ(ierr);
     ierr = PetscSectionSetChart(subCoordSection, firstSubPoint[0], firstSubPoint[0]+numSubPoints[0]);CHKERRQ(ierr);
     for (v = 0; v < numSubPoints[0]; ++v) {
       const PetscInt vertex    = subpoints[0][v];
 
       ierr = PetscSectionGetDof(coordSection, vertex, &dof);CHKERRQ(ierr);
       ierr = PetscSectionSetDof(subCoordSection, subvertex, dof);CHKERRQ(ierr);
+      ierr = PetscSectionSetFieldDof(subCoordSection, subvertex, 0, dof);CHKERRQ(ierr);
     }
     ierr = PetscSectionSetUp(subCoordSection);CHKERRQ(ierr);
     ierr = PetscSectionGetStorageSize(subCoordSection, &coordSize);CHKERRQ(ierr);
 PetscErrorCode DMPlexSetSubpointMap(DM dm, DMLabel subpointMap)
 {
   DM_Plex       *mesh = (DM_Plex *) dm->data;
+  DMLabel        tmp;
   PetscErrorCode ierr;
 
   PetscFunctionBegin;
   PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
-  ierr = DMLabelDestroy(&mesh->subpointMap);CHKERRQ(ierr);
+  tmp  = mesh->subpointMap;
   mesh->subpointMap = subpointMap;
-  ++mesh->subpointMap->refct;CHKERRQ(ierr);
+  ++mesh->subpointMap->refct;
+  ierr = DMLabelDestroy(&tmp);CHKERRQ(ierr);
   PetscFunctionReturn(0);
 }