Actual source code: plexinterpolate.c

petsc-3.8.3 2017-12-09
Report Typos and Errors
  1:  #include <petsc/private/dmpleximpl.h>

  3: /*
  4:   DMPlexGetFaces_Internal - Gets groups of vertices that correspond to faces for the given cell
  5: */
  6: PetscErrorCode DMPlexGetFaces_Internal(DM dm, PetscInt dim, PetscInt p, PetscInt *numFaces, PetscInt *faceSize, const PetscInt *faces[])
  7: {
  8:   const PetscInt *cone = NULL;
  9:   PetscInt        maxConeSize, maxSupportSize, coneSize;
 10:   PetscErrorCode  ierr;

 14:   DMPlexGetMaxSizes(dm, &maxConeSize, &maxSupportSize);
 15:   DMPlexGetConeSize(dm, p, &coneSize);
 16:   DMPlexGetCone(dm, p, &cone);
 17:   DMPlexGetRawFaces_Internal(dm, dim, coneSize, cone, numFaces, faceSize, faces);
 18:   return(0);
 19: }

 21: /*
 22:   DMPlexRestoreFaces_Internal - Restores the array
 23: */
 24: PetscErrorCode DMPlexRestoreFaces_Internal(DM dm, PetscInt dim, PetscInt p, PetscInt *numFaces, PetscInt *faceSize, const PetscInt *faces[])
 25: {
 26:   PetscErrorCode  ierr;

 29:   DMRestoreWorkArray(dm, 0, PETSC_INT, (void *) faces);
 30:   return(0);
 31: }

 33: /*
 34:   DMPlexGetRawFaces_Internal - Gets groups of vertices that correspond to faces for the given cone
 35: */
 36: PetscErrorCode DMPlexGetRawFaces_Internal(DM dm, PetscInt dim, PetscInt coneSize, const PetscInt cone[], PetscInt *numFaces, PetscInt *faceSize, const PetscInt *faces[])
 37: {
 38:   PetscInt       *facesTmp;
 39:   PetscInt        maxConeSize, maxSupportSize;
 40:   PetscErrorCode  ierr;

 44:   DMPlexGetMaxSizes(dm, &maxConeSize, &maxSupportSize);
 45:   if (faces) {DMGetWorkArray(dm, PetscSqr(PetscMax(maxConeSize, maxSupportSize)), PETSC_INT, &facesTmp);}
 46:   switch (dim) {
 47:   case 1:
 48:     switch (coneSize) {
 49:     case 2:
 50:       if (faces) {
 51:         facesTmp[0] = cone[0]; facesTmp[1] = cone[1];
 52:         *faces = facesTmp;
 53:       }
 54:       if (numFaces) *numFaces         = 2;
 55:       if (faceSize) *faceSize         = 1;
 56:       break;
 57:     default:
 58:       SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Cone size %D not supported for dimension %D", coneSize, dim);
 59:     }
 60:     break;
 61:   case 2:
 62:     switch (coneSize) {
 63:     case 3:
 64:       if (faces) {
 65:         facesTmp[0] = cone[0]; facesTmp[1] = cone[1];
 66:         facesTmp[2] = cone[1]; facesTmp[3] = cone[2];
 67:         facesTmp[4] = cone[2]; facesTmp[5] = cone[0];
 68:         *faces = facesTmp;
 69:       }
 70:       if (numFaces) *numFaces         = 3;
 71:       if (faceSize) *faceSize         = 2;
 72:       break;
 73:     case 4:
 74:       /* Vertices follow right hand rule */
 75:       if (faces) {
 76:         facesTmp[0] = cone[0]; facesTmp[1] = cone[1];
 77:         facesTmp[2] = cone[1]; facesTmp[3] = cone[2];
 78:         facesTmp[4] = cone[2]; facesTmp[5] = cone[3];
 79:         facesTmp[6] = cone[3]; facesTmp[7] = cone[0];
 80:         *faces = facesTmp;
 81:       }
 82:       if (numFaces) *numFaces         = 4;
 83:       if (faceSize) *faceSize         = 2;
 84:       break;
 85:     default:
 86:       SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Cone size %D not supported for dimension %D", coneSize, dim);
 87:     }
 88:     break;
 89:   case 3:
 90:     switch (coneSize) {
 91:     case 3:
 92:       if (faces) {
 93:         facesTmp[0] = cone[0]; facesTmp[1] = cone[1];
 94:         facesTmp[2] = cone[1]; facesTmp[3] = cone[2];
 95:         facesTmp[4] = cone[2]; facesTmp[5] = cone[0];
 96:         *faces = facesTmp;
 97:       }
 98:       if (numFaces) *numFaces         = 3;
 99:       if (faceSize) *faceSize         = 2;
100:       break;
101:     case 4:
102:       /* Vertices of first face follow right hand rule and normal points away from last vertex */
103:       if (faces) {
104:         facesTmp[0] = cone[0]; facesTmp[1]  = cone[1]; facesTmp[2]  = cone[2];
105:         facesTmp[3] = cone[0]; facesTmp[4]  = cone[3]; facesTmp[5]  = cone[1];
106:         facesTmp[6] = cone[0]; facesTmp[7]  = cone[2]; facesTmp[8]  = cone[3];
107:         facesTmp[9] = cone[2]; facesTmp[10] = cone[1]; facesTmp[11] = cone[3];
108:         *faces = facesTmp;
109:       }
110:       if (numFaces) *numFaces         = 4;
111:       if (faceSize) *faceSize         = 3;
112:       break;
113:     case 8:
114:       if (faces) {
115:         facesTmp[0]  = cone[0]; facesTmp[1]  = cone[1]; facesTmp[2]  = cone[2]; facesTmp[3]  = cone[3]; /* Bottom */
116:         facesTmp[4]  = cone[4]; facesTmp[5]  = cone[5]; facesTmp[6]  = cone[6]; facesTmp[7]  = cone[7]; /* Top */
117:         facesTmp[8]  = cone[0]; facesTmp[9]  = cone[3]; facesTmp[10] = cone[5]; facesTmp[11] = cone[4]; /* Front */
118:         facesTmp[12] = cone[2]; facesTmp[13] = cone[1]; facesTmp[14] = cone[7]; facesTmp[15] = cone[6]; /* Back */
119:         facesTmp[16] = cone[3]; facesTmp[17] = cone[2]; facesTmp[18] = cone[6]; facesTmp[19] = cone[5]; /* Right */
120:         facesTmp[20] = cone[0]; facesTmp[21] = cone[4]; facesTmp[22] = cone[7]; facesTmp[23] = cone[1]; /* Left */
121:         *faces = facesTmp;
122:       }
123:       if (numFaces) *numFaces         = 6;
124:       if (faceSize) *faceSize         = 4;
125:       break;
126:     default:
127:       SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Cone size %D not supported for dimension %D", coneSize, dim);
128:     }
129:     break;
130:   default:
131:     SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Dimension %D not supported", dim);
132:   }
133:   return(0);
134: }

136: /* This interpolates faces for cells at some stratum */
137: static PetscErrorCode DMPlexInterpolateFaces_Internal(DM dm, PetscInt cellDepth, DM idm)
138: {
139:   DMLabel        subpointMap;
140:   PetscHashIJKL  faceTable;
141:   PetscInt      *pStart, *pEnd;
142:   PetscInt       cellDim, depth, faceDepth = cellDepth, numPoints = 0, faceSizeAll = 0, face, c, d;

146:   DMGetDimension(dm, &cellDim);
147:   /* HACK: I need a better way to determine face dimension, or an alternative to GetFaces() */
148:   DMPlexGetSubpointMap(dm, &subpointMap);
149:   if (subpointMap) ++cellDim;
150:   DMPlexGetDepth(dm, &depth);
151:   ++depth;
152:   ++cellDepth;
153:   cellDim -= depth - cellDepth;
154:   PetscMalloc2(depth+1,&pStart,depth+1,&pEnd);
155:   for (d = depth-1; d >= faceDepth; --d) {
156:     DMPlexGetDepthStratum(dm, d, &pStart[d+1], &pEnd[d+1]);
157:   }
158:   DMPlexGetDepthStratum(dm, -1, NULL, &pStart[faceDepth]);
159:   pEnd[faceDepth] = pStart[faceDepth];
160:   for (d = faceDepth-1; d >= 0; --d) {
161:     DMPlexGetDepthStratum(dm, d, &pStart[d], &pEnd[d]);
162:   }
163:   if (pEnd[cellDepth] > pStart[cellDepth]) {DMPlexGetFaces_Internal(dm, cellDim, pStart[cellDepth], NULL, &faceSizeAll, NULL);}
164:   if (faceSizeAll > 4) SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "Do not support interpolation of meshes with faces of %D vertices", faceSizeAll);
165:   PetscHashIJKLCreate(&faceTable);
166:   for (c = pStart[cellDepth], face = pStart[faceDepth]; c < pEnd[cellDepth]; ++c) {
167:     const PetscInt *cellFaces;
168:     PetscInt        numCellFaces, faceSize, cf;

170:     DMPlexGetFaces_Internal(dm, cellDim, c, &numCellFaces, &faceSize, &cellFaces);
171:     if (faceSize != faceSizeAll) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Inconsistent face for cell %D of size %D != %D", c, faceSize, faceSizeAll);
172:     for (cf = 0; cf < numCellFaces; ++cf) {
173:       const PetscInt   *cellFace = &cellFaces[cf*faceSize];
174:       PetscHashIJKLKey  key;
175:       PetscHashIJKLIter missing, iter;

177:       if (faceSize == 2) {
178:         key.i = PetscMin(cellFace[0], cellFace[1]);
179:         key.j = PetscMax(cellFace[0], cellFace[1]);
180:         key.k = 0;
181:         key.l = 0;
182:       } else {
183:         key.i = cellFace[0]; key.j = cellFace[1]; key.k = cellFace[2]; key.l = faceSize > 3 ? cellFace[3] : 0;
184:         PetscSortInt(faceSize, (PetscInt *) &key);
185:       }
186:       PetscHashIJKLPut(faceTable, key, &missing, &iter);
187:       if (missing) {PetscHashIJKLSet(faceTable, iter, face++);}
188:     }
189:     DMPlexRestoreFaces_Internal(dm, cellDim, c, &numCellFaces, &faceSize, &cellFaces);
190:   }
191:   pEnd[faceDepth] = face;
192:   PetscHashIJKLDestroy(&faceTable);
193:   /* Count new points */
194:   for (d = 0; d <= depth; ++d) {
195:     numPoints += pEnd[d]-pStart[d];
196:   }
197:   DMPlexSetChart(idm, 0, numPoints);
198:   /* Set cone sizes */
199:   for (d = 0; d <= depth; ++d) {
200:     PetscInt coneSize, p;

202:     if (d == faceDepth) {
203:       for (p = pStart[d]; p < pEnd[d]; ++p) {
204:         /* I see no way to do this if we admit faces of different shapes */
205:         DMPlexSetConeSize(idm, p, faceSizeAll);
206:       }
207:     } else if (d == cellDepth) {
208:       for (p = pStart[d]; p < pEnd[d]; ++p) {
209:         /* Number of cell faces may be different from number of cell vertices*/
210:         DMPlexGetFaces_Internal(dm, cellDim, p, &coneSize, NULL, NULL);
211:         DMPlexSetConeSize(idm, p, coneSize);
212:       }
213:     } else {
214:       for (p = pStart[d]; p < pEnd[d]; ++p) {
215:         DMPlexGetConeSize(dm, p, &coneSize);
216:         DMPlexSetConeSize(idm, p, coneSize);
217:       }
218:     }
219:   }
220:   DMSetUp(idm);
221:   /* Get face cones from subsets of cell vertices */
222:   if (faceSizeAll > 4) SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "Do not support interpolation of meshes with faces of %D vertices", faceSizeAll);
223:   PetscHashIJKLCreate(&faceTable);
224:   for (d = depth; d > cellDepth; --d) {
225:     const PetscInt *cone;
226:     PetscInt        p;

228:     for (p = pStart[d]; p < pEnd[d]; ++p) {
229:       DMPlexGetCone(dm, p, &cone);
230:       DMPlexSetCone(idm, p, cone);
231:       DMPlexGetConeOrientation(dm, p, &cone);
232:       DMPlexSetConeOrientation(idm, p, cone);
233:     }
234:   }
235:   for (c = pStart[cellDepth], face = pStart[faceDepth]; c < pEnd[cellDepth]; ++c) {
236:     const PetscInt *cellFaces;
237:     PetscInt        numCellFaces, faceSize, cf;

239:     DMPlexGetFaces_Internal(dm, cellDim, c, &numCellFaces, &faceSize, &cellFaces);
240:     if (faceSize != faceSizeAll) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Inconsistent face for cell %D of size %D != %D", c, faceSize, faceSizeAll);
241:     for (cf = 0; cf < numCellFaces; ++cf) {
242:       const PetscInt  *cellFace = &cellFaces[cf*faceSize];
243:       PetscHashIJKLKey key;
244:       PetscHashIJKLIter missing, iter;

246:       if (faceSize == 2) {
247:         key.i = PetscMin(cellFace[0], cellFace[1]);
248:         key.j = PetscMax(cellFace[0], cellFace[1]);
249:         key.k = 0;
250:         key.l = 0;
251:       } else {
252:         key.i = cellFace[0]; key.j = cellFace[1]; key.k = cellFace[2]; key.l = faceSize > 3 ? cellFace[3] : 0;
253:         PetscSortInt(faceSize, (PetscInt *) &key);
254:       }
255:       PetscHashIJKLPut(faceTable, key, &missing, &iter);
256:       if (missing) {
257:         DMPlexSetCone(idm, face, cellFace);
258:         PetscHashIJKLSet(faceTable, iter, face);
259:         DMPlexInsertCone(idm, c, cf, face++);
260:       } else {
261:         const PetscInt *cone;
262:         PetscInt        coneSize, ornt, i, j, f;

264:         PetscHashIJKLGet(faceTable, iter, &f);
265:         DMPlexInsertCone(idm, c, cf, f);
266:         /* Orient face: Do not allow reverse orientation at the first vertex */
267:         DMPlexGetConeSize(idm, f, &coneSize);
268:         DMPlexGetCone(idm, f, &cone);
269:         if (coneSize != faceSize) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid number of face vertices %D for face %D should be %D", coneSize, f, faceSize);
270:         /* - First find the initial vertex */
271:         for (i = 0; i < faceSize; ++i) if (cellFace[0] == cone[i]) break;
272:         /* - Try forward comparison */
273:         for (j = 0; j < faceSize; ++j) if (cellFace[j] != cone[(i+j)%faceSize]) break;
274:         if (j == faceSize) {
275:           if ((faceSize == 2) && (i == 1)) ornt = -2;
276:           else                             ornt = i;
277:         } else {
278:           /* - Try backward comparison */
279:           for (j = 0; j < faceSize; ++j) if (cellFace[j] != cone[(i+faceSize-j)%faceSize]) break;
280:           if (j == faceSize) {
281:             if (i == 0) ornt = -faceSize;
282:             else        ornt = -i;
283:           } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Could not determine face orientation");
284:         }
285:         DMPlexInsertConeOrientation(idm, c, cf, ornt);
286:       }
287:     }
288:     DMPlexRestoreFaces_Internal(dm, cellDim, c, &numCellFaces, &faceSize, &cellFaces);
289:   }
290:   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]);
291:   PetscFree2(pStart,pEnd);
292:   PetscHashIJKLDestroy(&faceTable);
293:   PetscFree2(pStart,pEnd);
294:   DMPlexSymmetrize(idm);
295:   DMPlexStratify(idm);
296:   return(0);
297: }

299: /* This interpolates the PointSF in parallel following local interpolation */
300: static PetscErrorCode DMPlexInterpolatePointSF(DM dm, PetscSF pointSF, PetscInt depth)
301: {
302:   PetscMPIInt        size, rank;
303:   PetscInt           p, c, d, dof, offset;
304:   PetscInt           numLeaves, numRoots, candidatesSize, candidatesRemoteSize;
305:   const PetscInt    *localPoints;
306:   const PetscSFNode *remotePoints;
307:   PetscSFNode       *candidates, *candidatesRemote, *claims;
308:   PetscSection       candidateSection, candidateSectionRemote, claimSection;
309:   PetscHashI         leafhash;
310:   PetscHashIJ        roothash;
311:   PetscHashIJKey     key;
312:   PetscErrorCode     ierr;

315:   MPI_Comm_size(PetscObjectComm((PetscObject) dm), &size);
316:   MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank);
317:   PetscSFGetGraph(pointSF, &numRoots, &numLeaves, &localPoints, &remotePoints);
318:   if (size < 2 || numRoots < 0) return(0);
319:   PetscLogEventBegin(DMPLEX_InterpolateSF,dm,0,0,0);
320:   /* Build hashes of points in the SF for efficient lookup */
321:   PetscHashICreate(leafhash);
322:   PetscHashIJCreate(&roothash);
323:   PetscHashIJSetMultivalued(roothash, PETSC_FALSE);
324:   for (p = 0; p < numLeaves; ++p) {
325:     PetscHashIAdd(leafhash, localPoints[p], p);
326:     key.i = remotePoints[p].index; key.j = remotePoints[p].rank;
327:     PetscHashIJAdd(roothash, key, p);
328:   }
329:   /* Build a section / SFNode array of candidate points in the single-level adjacency of leaves,
330:      where each candidate is defined by the root entry for the other vertex that defines the edge. */
331:   PetscSectionCreate(PetscObjectComm((PetscObject) dm), &candidateSection);
332:   PetscSectionSetChart(candidateSection, 0, numRoots);
333:   {
334:     PetscInt leaf, root, idx, a, *adj = NULL;
335:     for (p = 0; p < numLeaves; ++p) {
336:       PetscInt adjSize = PETSC_DETERMINE;
337:       DMPlexGetAdjacency_Internal(dm, localPoints[p], PETSC_FALSE, PETSC_FALSE, PETSC_FALSE, &adjSize, &adj);
338:       for (a = 0; a < adjSize; ++a) {
339:         PetscHashIMap(leafhash, adj[a], leaf);
340:         if (leaf >= 0) {PetscSectionAddDof(candidateSection, localPoints[p], 1);}
341:       }
342:     }
343:     PetscSectionSetUp(candidateSection);
344:     PetscSectionGetStorageSize(candidateSection, &candidatesSize);
345:     PetscMalloc1(candidatesSize, &candidates);
346:     for (p = 0; p < numLeaves; ++p) {
347:       PetscInt adjSize = PETSC_DETERMINE;
348:       PetscSectionGetOffset(candidateSection, localPoints[p], &offset);
349:       DMPlexGetAdjacency_Internal(dm, localPoints[p], PETSC_FALSE, PETSC_FALSE, PETSC_FALSE, &adjSize, &adj);
350:       for (idx = 0, a = 0; a < adjSize; ++a) {
351:         PetscHashIMap(leafhash, adj[a], root);
352:         if (root >= 0) candidates[offset+idx++] = remotePoints[root];
353:       }
354:     }
355:     PetscFree(adj);
356:   }
357:   /* Gather candidate section / array pair into the root partition via inverse(multi(pointSF)). */
358:   {
359:     PetscSF   sfMulti, sfInverse, sfCandidates;
360:     PetscInt *remoteOffsets;
361:     PetscSFGetMultiSF(pointSF, &sfMulti);
362:     PetscSFCreateInverseSF(sfMulti, &sfInverse);
363:     PetscSectionCreate(PetscObjectComm((PetscObject) dm), &candidateSectionRemote);
364:     PetscSFDistributeSection(sfInverse, candidateSection, &remoteOffsets, candidateSectionRemote);
365:     PetscSFCreateSectionSF(sfInverse, candidateSection, remoteOffsets, candidateSectionRemote, &sfCandidates);
366:     PetscSectionGetStorageSize(candidateSectionRemote, &candidatesRemoteSize);
367:     PetscMalloc1(candidatesRemoteSize, &candidatesRemote);
368:     PetscSFBcastBegin(sfCandidates, MPIU_2INT, candidates, candidatesRemote);
369:     PetscSFBcastEnd(sfCandidates, MPIU_2INT, candidates, candidatesRemote);
370:     PetscSFDestroy(&sfInverse);
371:     PetscSFDestroy(&sfCandidates);
372:     PetscFree(remoteOffsets);
373:   }
374:   /* Walk local roots and check for each remote candidate whether we know all required points,
375:      either from owning it or having a root entry in the point SF. If we do we place a claim
376:      by replacing the vertex number with our edge ID. */
377:   {
378:     PetscInt        idx, root, joinSize, vertices[2];
379:     const PetscInt *rootdegree, *join = NULL;
380:     PetscSFComputeDegreeBegin(pointSF, &rootdegree);
381:     PetscSFComputeDegreeEnd(pointSF, &rootdegree);
382:     /* Loop remote edge connections and put in a claim if both vertices are known */
383:     for (idx = 0, p = 0; p < numRoots; ++p) {
384:       for (d = 0; d < rootdegree[p]; ++d) {
385:         PetscSectionGetDof(candidateSectionRemote, idx, &dof);
386:         PetscSectionGetOffset(candidateSectionRemote, idx, &offset);
387:         for (c = 0; c < dof; ++c) {
388:           /* We own both vertices, so we claim the edge by replacing vertex with edge */
389:           if (candidatesRemote[offset+c].rank == rank) {
390:             vertices[0] = p; vertices[1] = candidatesRemote[offset+c].index;
391:             DMPlexGetJoin(dm, 2, vertices, &joinSize, &join);
392:             if (joinSize == 1) candidatesRemote[offset+c].index = join[0];
393:             DMPlexRestoreJoin(dm, 2, vertices, &joinSize, &join);
394:             continue;
395:           }
396:           /* If we own one vertex and share a root with the other, we claim it */
397:           key.i = candidatesRemote[offset+c].index; key.j = candidatesRemote[offset+c].rank;
398:           PetscHashIJGet(roothash, key, &root);
399:           if (root >= 0) {
400:             vertices[0] = p; vertices[1] = localPoints[root];
401:             DMPlexGetJoin(dm, 2, vertices, &joinSize, &join);
402:             if (joinSize == 1) {
403:               candidatesRemote[offset+c].index = join[0];
404:               candidatesRemote[offset+c].rank = rank;
405:             }
406:             DMPlexRestoreJoin(dm, 2, vertices, &joinSize, &join);
407:           }
408:         }
409:         idx++;
410:       }
411:     }
412:   }
413:   /* Push claims back to receiver via the MultiSF and derive new pointSF mapping on receiver */
414:   {
415:     PetscSF         sfMulti, sfClaims, sfPointNew;
416:     PetscHashI      claimshash;
417:     PetscInt        size, pStart, pEnd, root, joinSize, numLocalNew;
418:     PetscInt       *remoteOffsets, *localPointsNew, vertices[2];
419:     const PetscInt *join = NULL;
420:     PetscSFNode    *remotePointsNew;
421:     PetscSFGetMultiSF(pointSF, &sfMulti);
422:     PetscSectionCreate(PetscObjectComm((PetscObject) dm), &claimSection);
423:     PetscSFDistributeSection(sfMulti, candidateSectionRemote, &remoteOffsets, claimSection);
424:     PetscSFCreateSectionSF(sfMulti, candidateSectionRemote, remoteOffsets, claimSection, &sfClaims);
425:     PetscSectionGetStorageSize(claimSection, &size);
426:     PetscMalloc1(size, &claims);
427:     PetscSFBcastBegin(sfClaims, MPIU_2INT, candidatesRemote, claims);
428:     PetscSFBcastEnd(sfClaims, MPIU_2INT, candidatesRemote, claims);
429:     PetscSFDestroy(&sfClaims);
430:     PetscFree(remoteOffsets);
431:     /* Walk the original section of local supports and add an SF entry for each updated item */
432:     PetscHashICreate(claimshash);
433:     for (p = 0; p < numRoots; ++p) {
434:       PetscSectionGetDof(candidateSection, p, &dof);
435:       PetscSectionGetOffset(candidateSection, p, &offset);
436:       for (d = 0; d < dof; ++d) {
437:         if (candidates[offset+d].index != claims[offset+d].index) {
438:           key.i = candidates[offset+d].index; key.j = candidates[offset+d].rank;
439:           PetscHashIJGet(roothash, key, &root);
440:           if (root >= 0) {
441:             vertices[0] = p; vertices[1] = localPoints[root];
442:             DMPlexGetJoin(dm, 2, vertices, &joinSize, &join);
443:             if (joinSize == 1) PetscHashIAdd(claimshash, join[0], offset+d);
444:             DMPlexRestoreJoin(dm, 2, vertices, &joinSize, &join);
445:           }
446:         }
447:       }
448:     }
449:     /* Create new pointSF from hashed claims */
450:     PetscHashISize(claimshash, numLocalNew);
451:     DMPlexGetChart(dm, &pStart, &pEnd);
452:     PetscMalloc1(numLeaves + numLocalNew, &localPointsNew);
453:     PetscMalloc1(numLeaves + numLocalNew, &remotePointsNew);
454:     for (p = 0; p < numLeaves; ++p) {
455:       localPointsNew[p] = localPoints[p];
456:       remotePointsNew[p].index = remotePoints[p].index;
457:       remotePointsNew[p].rank = remotePoints[p].rank;
458:     }
459:     p = numLeaves;
460:     PetscHashIGetKeys(claimshash, &p, localPointsNew);
461:     PetscSortInt(numLocalNew,&localPointsNew[numLeaves]);
462:     for (p = numLeaves; p < numLeaves + numLocalNew; ++p) {
463:       PetscHashIMap(claimshash, localPointsNew[p], offset);
464:       remotePointsNew[p] = claims[offset];
465:     }
466:     PetscSFCreate(PetscObjectComm((PetscObject) dm), &sfPointNew);
467:     PetscSFSetGraph(sfPointNew, pEnd-pStart, numLeaves+numLocalNew, localPointsNew, PETSC_OWN_POINTER, remotePointsNew, PETSC_OWN_POINTER);
468:     DMSetPointSF(dm, sfPointNew);
469:     PetscSFDestroy(&sfPointNew);
470:     PetscHashIDestroy(claimshash);
471:   }
472:   PetscHashIDestroy(leafhash);
473:   PetscHashIJDestroy(&roothash);
474:   PetscSectionDestroy(&candidateSection);
475:   PetscSectionDestroy(&candidateSectionRemote);
476:   PetscSectionDestroy(&claimSection);
477:   PetscFree(candidates);
478:   PetscFree(candidatesRemote);
479:   PetscFree(claims);
480:   PetscLogEventEnd(DMPLEX_InterpolateSF,dm,0,0,0);
481:   return(0);
482: }

484: /*@C
485:   DMPlexInterpolate - Take in a cell-vertex mesh and return one with all intermediate faces, edges, etc.

487:   Collective on DM

489:   Input Parameters:
490: + dm - The DMPlex object with only cells and vertices
491: - dmInt - If NULL a new DM is created, otherwise the interpolated DM is put into the given DM

493:   Output Parameter:
494: . dmInt - The complete DMPlex object

496:   Level: intermediate

498: .keywords: mesh
499: .seealso: DMPlexUninterpolate(), DMPlexCreateFromCellList()
500: @*/
501: PetscErrorCode DMPlexInterpolate(DM dm, DM *dmInt)
502: {
503:   DM             idm, odm = dm;
504:   PetscSF        sfPoint;
505:   PetscInt       depth, dim, d;

509:   PetscLogEventBegin(DMPLEX_Interpolate,dm,0,0,0);
510:   DMPlexGetDepth(dm, &depth);
511:   DMGetDimension(dm, &dim);
512:   if (dim <= 1) {
513:     PetscObjectReference((PetscObject) dm);
514:     idm  = dm;
515:   }
516:   for (d = 1; d < dim; ++d) {
517:     /* Create interpolated mesh */
518:     if ((d == dim-1) && *dmInt) {idm  = *dmInt;}
519:     else                        {DMCreate(PetscObjectComm((PetscObject)dm), &idm);}
520:     DMSetType(idm, DMPLEX);
521:     DMSetDimension(idm, dim);
522:     if (depth > 0) {
523:       DMPlexInterpolateFaces_Internal(odm, 1, idm);
524:       DMGetPointSF(odm, &sfPoint);
525:       DMPlexInterpolatePointSF(idm, sfPoint, depth);
526:     }
527:     if (odm != dm) {DMDestroy(&odm);}
528:     odm  = idm;
529:   }
530:   *dmInt = idm;
531:   PetscLogEventEnd(DMPLEX_Interpolate,dm,0,0,0);
532:   return(0);
533: }

535: /*@
536:   DMPlexCopyCoordinates - Copy coordinates from one mesh to another with the same vertices

538:   Collective on DM

540:   Input Parameter:
541: . dmA - The DMPlex object with initial coordinates

543:   Output Parameter:
544: . dmB - The DMPlex object with copied coordinates

546:   Level: intermediate

548:   Note: This is typically used when adding pieces other than vertices to a mesh

550: .keywords: mesh
551: .seealso: DMCopyLabels(), DMGetCoordinates(), DMGetCoordinatesLocal(), DMGetCoordinateDM(), DMGetCoordinateSection()
552: @*/
553: PetscErrorCode DMPlexCopyCoordinates(DM dmA, DM dmB)
554: {
555:   Vec            coordinatesA, coordinatesB;
556:   PetscSection   coordSectionA, coordSectionB;
557:   PetscScalar   *coordsA, *coordsB;
558:   PetscInt       spaceDim, vStartA, vStartB, vEndA, vEndB, coordSizeB, v, d;

562:   if (dmA == dmB) return(0);
563:   DMPlexGetDepthStratum(dmA, 0, &vStartA, &vEndA);
564:   DMPlexGetDepthStratum(dmB, 0, &vStartB, &vEndB);
565:   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);
566:   DMGetCoordinateSection(dmA, &coordSectionA);
567:   DMGetCoordinateSection(dmB, &coordSectionB);
568:   if (coordSectionA == coordSectionB) return(0);
569:   if (!coordSectionB) {
570:     PetscInt dim;

572:     PetscSectionCreate(PetscObjectComm((PetscObject) coordSectionA), &coordSectionB);
573:     DMGetCoordinateDim(dmA, &dim);
574:     DMSetCoordinateSection(dmB, dim, coordSectionB);
575:     PetscObjectDereference((PetscObject) coordSectionB);
576:   }
577:   PetscSectionSetNumFields(coordSectionB, 1);
578:   PetscSectionGetFieldComponents(coordSectionA, 0, &spaceDim);
579:   PetscSectionSetFieldComponents(coordSectionB, 0, spaceDim);
580:   PetscSectionSetChart(coordSectionB, vStartB, vEndB);
581:   for (v = vStartB; v < vEndB; ++v) {
582:     PetscSectionSetDof(coordSectionB, v, spaceDim);
583:     PetscSectionSetFieldDof(coordSectionB, v, 0, spaceDim);
584:   }
585:   PetscSectionSetUp(coordSectionB);
586:   PetscSectionGetStorageSize(coordSectionB, &coordSizeB);
587:   DMGetCoordinatesLocal(dmA, &coordinatesA);
588:   VecCreate(PETSC_COMM_SELF, &coordinatesB);
589:   PetscObjectSetName((PetscObject) coordinatesB, "coordinates");
590:   VecSetSizes(coordinatesB, coordSizeB, PETSC_DETERMINE);
591:   VecSetType(coordinatesB,VECSTANDARD);
592:   VecGetArray(coordinatesA, &coordsA);
593:   VecGetArray(coordinatesB, &coordsB);
594:   for (v = 0; v < vEndB-vStartB; ++v) {
595:     for (d = 0; d < spaceDim; ++d) {
596:       coordsB[v*spaceDim+d] = coordsA[v*spaceDim+d];
597:     }
598:   }
599:   VecRestoreArray(coordinatesA, &coordsA);
600:   VecRestoreArray(coordinatesB, &coordsB);
601:   DMSetCoordinatesLocal(dmB, coordinatesB);
602:   VecDestroy(&coordinatesB);
603:   return(0);
604: }

606: /*@
607:   DMPlexUninterpolate - Take in a mesh with all intermediate faces, edges, etc. and return a cell-vertex mesh

609:   Collective on DM

611:   Input Parameter:
612: . dm - The complete DMPlex object

614:   Output Parameter:
615: . dmUnint - The DMPlex object with only cells and vertices

617:   Level: intermediate

619: .keywords: mesh
620: .seealso: DMPlexInterpolate(), DMPlexCreateFromCellList()
621: @*/
622: PetscErrorCode DMPlexUninterpolate(DM dm, DM *dmUnint)
623: {
624:   DM             udm;
625:   PetscInt       dim, vStart, vEnd, cStart, cEnd, c, maxConeSize = 0, *cone;

629:   DMGetDimension(dm, &dim);
630:   if (dim <= 1) {
631:     PetscObjectReference((PetscObject) dm);
632:     *dmUnint = dm;
633:     return(0);
634:   }
635:   DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
636:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
637:   DMCreate(PetscObjectComm((PetscObject) dm), &udm);
638:   DMSetType(udm, DMPLEX);
639:   DMSetDimension(udm, dim);
640:   DMPlexSetChart(udm, cStart, vEnd);
641:   for (c = cStart; c < cEnd; ++c) {
642:     PetscInt *closure = NULL, closureSize, cl, coneSize = 0;

644:     DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
645:     for (cl = 0; cl < closureSize*2; cl += 2) {
646:       const PetscInt p = closure[cl];

648:       if ((p >= vStart) && (p < vEnd)) ++coneSize;
649:     }
650:     DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
651:     DMPlexSetConeSize(udm, c, coneSize);
652:     maxConeSize = PetscMax(maxConeSize, coneSize);
653:   }
654:   DMSetUp(udm);
655:   PetscMalloc1(maxConeSize, &cone);
656:   for (c = cStart; c < cEnd; ++c) {
657:     PetscInt *closure = NULL, closureSize, cl, coneSize = 0;

659:     DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
660:     for (cl = 0; cl < closureSize*2; cl += 2) {
661:       const PetscInt p = closure[cl];

663:       if ((p >= vStart) && (p < vEnd)) cone[coneSize++] = p;
664:     }
665:     DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
666:     DMPlexSetCone(udm, c, cone);
667:   }
668:   PetscFree(cone);
669:   DMPlexSymmetrize(udm);
670:   DMPlexStratify(udm);
671:   /* Reduce SF */
672:   {
673:     PetscSF            sfPoint, sfPointUn;
674:     const PetscSFNode *remotePoints;
675:     const PetscInt    *localPoints;
676:     PetscSFNode       *remotePointsUn;
677:     PetscInt          *localPointsUn;
678:     PetscInt           vEnd, numRoots, numLeaves, l;
679:     PetscInt           numLeavesUn = 0, n = 0;
680:     PetscErrorCode     ierr;

682:     /* Get original SF information */
683:     DMGetPointSF(dm, &sfPoint);
684:     DMGetPointSF(udm, &sfPointUn);
685:     DMPlexGetDepthStratum(dm, 0, NULL, &vEnd);
686:     PetscSFGetGraph(sfPoint, &numRoots, &numLeaves, &localPoints, &remotePoints);
687:     /* Allocate space for cells and vertices */
688:     for (l = 0; l < numLeaves; ++l) if (localPoints[l] < vEnd) numLeavesUn++;
689:     /* Fill in leaves */
690:     if (vEnd >= 0) {
691:       PetscMalloc1(numLeavesUn, &remotePointsUn);
692:       PetscMalloc1(numLeavesUn, &localPointsUn);
693:       for (l = 0; l < numLeaves; l++) {
694:         if (localPoints[l] < vEnd) {
695:           localPointsUn[n]        = localPoints[l];
696:           remotePointsUn[n].rank  = remotePoints[l].rank;
697:           remotePointsUn[n].index = remotePoints[l].index;
698:           ++n;
699:         }
700:       }
701:       if (n != numLeavesUn) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Inconsistent number of leaves %d != %d", n, numLeavesUn);
702:       PetscSFSetGraph(sfPointUn, vEnd, numLeavesUn, localPointsUn, PETSC_OWN_POINTER, remotePointsUn, PETSC_OWN_POINTER);
703:     }
704:   }
705:   *dmUnint = udm;
706:   return(0);
707: }