Actual source code: plexinterpolate.c

petsc-3.11.4 2019-09-28
Report Typos and Errors
  1:  #include <petsc/private/dmpleximpl.h>
  2:  #include <petsc/private/hashmapi.h>
  3:  #include <petsc/private/hashmapij.h>

  5: /* HashIJKL */

  7:  #include <petsc/private/hashmap.h>

  9: typedef struct _PetscHashIJKLKey { PetscInt i, j, k, l; } PetscHashIJKLKey;

 11: #define PetscHashIJKLKeyHash(key) \
 12:   PetscHashCombine(PetscHashCombine(PetscHashInt((key).i),PetscHashInt((key).j)), \
 13:                    PetscHashCombine(PetscHashInt((key).k),PetscHashInt((key).l)))

 15: #define PetscHashIJKLKeyEqual(k1,k2) \
 16:   (((k1).i==(k2).i) ? ((k1).j==(k2).j) ? ((k1).k==(k2).k) ? ((k1).l==(k2).l) : 0 : 0 : 0)

 18: PETSC_HASH_MAP(HashIJKL, PetscHashIJKLKey, PetscInt, PetscHashIJKLKeyHash, PetscHashIJKLKeyEqual, -1)


 21: /*
 22:   DMPlexGetFaces_Internal - Gets groups of vertices that correspond to faces for the given cell
 23:   This assumes that the mesh is not interpolated from the depth of point p to the vertices
 24: */
 25: PetscErrorCode DMPlexGetFaces_Internal(DM dm, PetscInt dim, PetscInt p, PetscInt *numFaces, PetscInt *faceSize, const PetscInt *faces[])
 26: {
 27:   const PetscInt *cone = NULL;
 28:   PetscInt        coneSize;
 29:   PetscErrorCode  ierr;

 33:   DMPlexGetConeSize(dm, p, &coneSize);
 34:   DMPlexGetCone(dm, p, &cone);
 35:   DMPlexGetRawFaces_Internal(dm, dim, coneSize, cone, numFaces, faceSize, faces);
 36:   return(0);
 37: }

 39: /*
 40:   DMPlexRestoreFaces_Internal - Restores the array
 41: */
 42: PetscErrorCode DMPlexRestoreFaces_Internal(DM dm, PetscInt dim, PetscInt p, PetscInt *numFaces, PetscInt *faceSize, const PetscInt *faces[])
 43: {
 44:   PetscErrorCode  ierr;

 47:   if (faces) { DMRestoreWorkArray(dm, 0, MPIU_INT, (void *) faces); }
 48:   return(0);
 49: }

 51: /*
 52:   DMPlexGetRawFaces_Internal - Gets groups of vertices that correspond to faces for the given cone
 53: */
 54: PetscErrorCode DMPlexGetRawFaces_Internal(DM dm, PetscInt dim, PetscInt coneSize, const PetscInt cone[], PetscInt *numFaces, PetscInt *faceSize, const PetscInt *faces[])
 55: {
 56:   PetscInt       *facesTmp;
 57:   PetscInt        maxConeSize, maxSupportSize;
 58:   PetscErrorCode  ierr;

 63:   DMPlexGetMaxSizes(dm, &maxConeSize, &maxSupportSize);
 64:   if (faces) {DMGetWorkArray(dm, PetscSqr(PetscMax(maxConeSize, maxSupportSize)), MPIU_INT, &facesTmp);}
 65:   switch (dim) {
 66:   case 1:
 67:     switch (coneSize) {
 68:     case 2:
 69:       if (faces) {
 70:         facesTmp[0] = cone[0]; facesTmp[1] = cone[1];
 71:         *faces = facesTmp;
 72:       }
 73:       if (numFaces) *numFaces = 2;
 74:       if (faceSize) *faceSize = 1;
 75:       break;
 76:     default:
 77:       SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Cone size %D not supported for dimension %D", coneSize, dim);
 78:     }
 79:     break;
 80:   case 2:
 81:     switch (coneSize) {
 82:     case 3:
 83:       if (faces) {
 84:         facesTmp[0] = cone[0]; facesTmp[1] = cone[1];
 85:         facesTmp[2] = cone[1]; facesTmp[3] = cone[2];
 86:         facesTmp[4] = cone[2]; facesTmp[5] = cone[0];
 87:         *faces = facesTmp;
 88:       }
 89:       if (numFaces) *numFaces = 3;
 90:       if (faceSize) *faceSize = 2;
 91:       break;
 92:     case 4:
 93:       /* Vertices follow right hand rule */
 94:       if (faces) {
 95:         facesTmp[0] = cone[0]; facesTmp[1] = cone[1];
 96:         facesTmp[2] = cone[1]; facesTmp[3] = cone[2];
 97:         facesTmp[4] = cone[2]; facesTmp[5] = cone[3];
 98:         facesTmp[6] = cone[3]; facesTmp[7] = cone[0];
 99:         *faces = facesTmp;
100:       }
101:       if (numFaces) *numFaces = 4;
102:       if (faceSize) *faceSize = 2;
103:       break;
104:     default:
105:       SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Cone size %D not supported for dimension %D", coneSize, dim);
106:     }
107:     break;
108:   case 3:
109:     switch (coneSize) {
110:     case 3:
111:       if (faces) {
112:         facesTmp[0] = cone[0]; facesTmp[1] = cone[1];
113:         facesTmp[2] = cone[1]; facesTmp[3] = cone[2];
114:         facesTmp[4] = cone[2]; facesTmp[5] = cone[0];
115:         *faces = facesTmp;
116:       }
117:       if (numFaces) *numFaces = 3;
118:       if (faceSize) *faceSize = 2;
119:       break;
120:     case 4:
121:       /* Vertices of first face follow right hand rule and normal points away from last vertex */
122:       if (faces) {
123:         facesTmp[0] = cone[0]; facesTmp[1]  = cone[1]; facesTmp[2]  = cone[2];
124:         facesTmp[3] = cone[0]; facesTmp[4]  = cone[3]; facesTmp[5]  = cone[1];
125:         facesTmp[6] = cone[0]; facesTmp[7]  = cone[2]; facesTmp[8]  = cone[3];
126:         facesTmp[9] = cone[2]; facesTmp[10] = cone[1]; facesTmp[11] = cone[3];
127:         *faces = facesTmp;
128:       }
129:       if (numFaces) *numFaces = 4;
130:       if (faceSize) *faceSize = 3;
131:       break;
132:     case 8:
133:       /*  7--------6
134:          /|       /|
135:         / |      / |
136:        4--------5  |
137:        |  |     |  |
138:        |  |     |  |
139:        |  1--------2
140:        | /      | /
141:        |/       |/
142:        0--------3
143:        */
144:       if (faces) {
145:         facesTmp[0]  = cone[0]; facesTmp[1]  = cone[1]; facesTmp[2]  = cone[2]; facesTmp[3]  = cone[3]; /* Bottom */
146:         facesTmp[4]  = cone[4]; facesTmp[5]  = cone[5]; facesTmp[6]  = cone[6]; facesTmp[7]  = cone[7]; /* Top */
147:         facesTmp[8]  = cone[0]; facesTmp[9]  = cone[3]; facesTmp[10] = cone[5]; facesTmp[11] = cone[4]; /* Front */
148:         facesTmp[12] = cone[2]; facesTmp[13] = cone[1]; facesTmp[14] = cone[7]; facesTmp[15] = cone[6]; /* Back */
149:         facesTmp[16] = cone[3]; facesTmp[17] = cone[2]; facesTmp[18] = cone[6]; facesTmp[19] = cone[5]; /* Right */
150:         facesTmp[20] = cone[0]; facesTmp[21] = cone[4]; facesTmp[22] = cone[7]; facesTmp[23] = cone[1]; /* Left */
151:         *faces = facesTmp;
152:       }
153:       if (numFaces) *numFaces = 6;
154:       if (faceSize) *faceSize = 4;
155:       break;
156:     default:
157:       SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Cone size %D not supported for dimension %D", coneSize, dim);
158:     }
159:     break;
160:   default:
161:     SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Dimension %D not supported", dim);
162:   }
163:   return(0);
164: }

166: /*
167:   DMPlexGetRawFacesHybrid_Internal - Gets groups of vertices that correspond to faces for the given cone using hybrid ordering (prisms)
168: */
169: static PetscErrorCode DMPlexGetRawFacesHybrid_Internal(DM dm, PetscInt dim, PetscInt coneSize, const PetscInt cone[], PetscInt *numFaces, PetscInt *numFacesNotH, PetscInt *faceSize, const PetscInt *faces[])
170: {
171:   PetscInt       *facesTmp;
172:   PetscInt        maxConeSize, maxSupportSize;
173:   PetscErrorCode  ierr;

178:   DMPlexGetMaxSizes(dm, &maxConeSize, &maxSupportSize);
179:   if (faces) {DMGetWorkArray(dm, PetscSqr(PetscMax(maxConeSize, maxSupportSize)), MPIU_INT, &facesTmp);}
180:   switch (dim) {
181:   case 1:
182:     switch (coneSize) {
183:     case 2:
184:       if (faces) {
185:         facesTmp[0] = cone[0]; facesTmp[1] = cone[1];
186:         *faces = facesTmp;
187:       }
188:       if (numFaces)     *numFaces = 2;
189:       if (numFacesNotH) *numFacesNotH = 2;
190:       if (faceSize)     *faceSize = 1;
191:       break;
192:     default:
193:       SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Cone size %D not supported for dimension %D", coneSize, dim);
194:     }
195:     break;
196:   case 2:
197:     switch (coneSize) {
198:     case 4:
199:       if (faces) {
200:         facesTmp[0] = cone[0]; facesTmp[1] = cone[1];
201:         facesTmp[2] = cone[2]; facesTmp[3] = cone[3];
202:         facesTmp[4] = cone[0]; facesTmp[5] = cone[2];
203:         facesTmp[6] = cone[1]; facesTmp[7] = cone[3];
204:         *faces = facesTmp;
205:       }
206:       if (numFaces)     *numFaces = 4;
207:       if (numFacesNotH) *numFacesNotH = 2;
208:       if (faceSize)     *faceSize = 2;
209:       break;
210:     default:
211:       SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Cone size %D not supported for dimension %D", coneSize, dim);
212:     }
213:     break;
214:   case 3:
215:     switch (coneSize) {
216:     case 6: /* triangular prism */
217:       if (faces) {
218:         facesTmp[0]  = cone[0]; facesTmp[1]  = cone[1]; facesTmp[2]  = cone[2]; facesTmp[3]  = -1;      /* Bottom */
219:         facesTmp[4]  = cone[3]; facesTmp[5]  = cone[4]; facesTmp[6]  = cone[5]; facesTmp[7]  = -1;      /* Top */
220:         facesTmp[8]  = cone[0]; facesTmp[9]  = cone[1]; facesTmp[10] = cone[3]; facesTmp[11] = cone[4]; /* Back left */
221:         facesTmp[12] = cone[1]; facesTmp[13] = cone[2]; facesTmp[14] = cone[4]; facesTmp[15] = cone[5]; /* Back right */
222:         facesTmp[16] = cone[2]; facesTmp[17] = cone[0]; facesTmp[18] = cone[5]; facesTmp[19] = cone[3]; /* Front */
223:         *faces = facesTmp;
224:       }
225:       if (numFaces)     *numFaces = 5;
226:       if (numFacesNotH) *numFacesNotH = 2;
227:       if (faceSize)     *faceSize = -4;
228:       break;
229:     default:
230:       SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Cone size %D not supported for dimension %D", coneSize, dim);
231:     }
232:     break;
233:   default:
234:     SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Dimension %D not supported", dim);
235:   }
236:   return(0);
237: }

239: static PetscErrorCode DMPlexRestoreRawFacesHybrid_Internal(DM dm, PetscInt dim, PetscInt coneSize, const PetscInt cone[], PetscInt *numFaces, PetscInt *numFacesNotH, PetscInt *faceSize, const PetscInt *faces[])
240: {
241:   PetscErrorCode  ierr;

244:   if (faces) { DMRestoreWorkArray(dm, 0, MPIU_INT, (void *) faces); }
245:   return(0);
246: }

248: static PetscErrorCode DMPlexGetFacesHybrid_Internal(DM dm, PetscInt dim, PetscInt p, PetscInt *numFaces, PetscInt *numFacesNotH, PetscInt *faceSize, const PetscInt *faces[])
249: {
250:   const PetscInt *cone = NULL;
251:   PetscInt        coneSize;
252:   PetscErrorCode  ierr;

256:   DMPlexGetConeSize(dm, p, &coneSize);
257:   DMPlexGetCone(dm, p, &cone);
258:   DMPlexGetRawFacesHybrid_Internal(dm, dim, coneSize, cone, numFaces, numFacesNotH, faceSize, faces);
259:   return(0);
260: }

262: /* This interpolates faces for cells at some stratum */
263: static PetscErrorCode DMPlexInterpolateFaces_Internal(DM dm, PetscInt cellDepth, DM idm)
264: {
265:   DMLabel        subpointMap;
266:   PetscHashIJKL  faceTable;
267:   PetscInt      *pStart, *pEnd;
268:   PetscInt       cellDim, depth, faceDepth = cellDepth, numPoints = 0, faceSizeAll = 0, face, c, d;
269:   PetscInt       coneSizeH = 0, faceSizeAllH = 0, numCellFacesH = 0, faceH, pMax = -1, dim, outerloop;
270:   PetscInt       cMax, fMax, eMax, vMax;

274:   DMGetDimension(dm, &cellDim);
275:   /* HACK: I need a better way to determine face dimension, or an alternative to GetFaces() */
276:   DMPlexGetSubpointMap(dm, &subpointMap);
277:   if (subpointMap) ++cellDim;
278:   DMPlexGetDepth(dm, &depth);
279:   ++depth;
280:   ++cellDepth;
281:   cellDim -= depth - cellDepth;
282:   PetscMalloc2(depth+1,&pStart,depth+1,&pEnd);
283:   for (d = depth-1; d >= faceDepth; --d) {
284:     DMPlexGetDepthStratum(dm, d, &pStart[d+1], &pEnd[d+1]);
285:   }
286:   DMPlexGetDepthStratum(dm, -1, NULL, &pStart[faceDepth]);
287:   pEnd[faceDepth] = pStart[faceDepth];
288:   for (d = faceDepth-1; d >= 0; --d) {
289:     DMPlexGetDepthStratum(dm, d, &pStart[d], &pEnd[d]);
290:   }
291:   cMax = fMax = eMax = vMax = PETSC_DETERMINE;
292:   DMGetDimension(dm, &dim);
293:   if (cellDim == dim) {
294:     DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);
295:     pMax = cMax;
296:   } else if (cellDim == dim -1) {
297:     DMPlexGetHybridBounds(dm, &cMax, &fMax, NULL, NULL);
298:     pMax = fMax;
299:   }
300:   pMax = pMax < 0 ? pEnd[cellDepth] : pMax;
301:   if (pMax < pEnd[cellDepth]) {
302:     const PetscInt *cellFaces, *cone;
303:     PetscInt        numCellFacesT, faceSize, cf;

305:     DMPlexGetConeSize(dm, pMax, &coneSizeH);
306:     DMPlexGetCone(dm, pMax, &cone);
307:     DMPlexGetRawFacesHybrid_Internal(dm, cellDim, coneSizeH, cone, &numCellFacesH, &numCellFacesT, &faceSize, &cellFaces);
308:     if (faceSize < 0) {
309:       PetscInt *sizes, minv, maxv;

311:       /* count vertices of hybrid and non-hybrid faces */
312:       PetscCalloc1(numCellFacesH, &sizes);
313:       for (cf = 0; cf < numCellFacesT; ++cf) { /* These are the non-hybrid faces */
314:         const PetscInt *cellFace = &cellFaces[-cf*faceSize];
315:         PetscInt       f;

317:         for (f = 0; f < -faceSize; ++f) sizes[cf] += (cellFace[f] >= 0 ? 1 : 0);
318:       }
319:       PetscSortInt(numCellFacesT, sizes);
320:       minv = sizes[0];
321:       maxv = sizes[PetscMax(numCellFacesT-1, 0)];
322:       if (minv != maxv) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_SUP, "Different number of vertices for non-hybrid face %D != %D", minv, maxv);
323:       faceSizeAll = minv;
324:       PetscMemzero(sizes, numCellFacesH*sizeof(PetscInt));
325:       for (cf = numCellFacesT; cf < numCellFacesH; ++cf) { /* These are the hybrid faces */
326:         const PetscInt *cellFace = &cellFaces[-cf*faceSize];
327:         PetscInt       f;

329:         for (f = 0; f < -faceSize; ++f) sizes[cf-numCellFacesT] += (cellFace[f] >= 0 ? 1 : 0);
330:       }
331:       PetscSortInt(numCellFacesH - numCellFacesT, sizes);
332:       minv = sizes[0];
333:       maxv = sizes[PetscMax(numCellFacesH - numCellFacesT-1, 0)];
334:       PetscFree(sizes);
335:       if (minv != maxv) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_SUP, "Different number of vertices for hybrid face %D != %D", minv, maxv);
336:       faceSizeAllH = minv;
337:     } else { /* the size of the faces in hybrid cells is the same */
338:       faceSizeAll = faceSizeAllH = faceSize;
339:     }
340:     DMPlexRestoreRawFacesHybrid_Internal(dm, cellDim, coneSizeH, cone, &numCellFacesH, &numCellFacesT, &faceSize, &cellFaces);
341:   } else if (pEnd[cellDepth] > pStart[cellDepth]) {
342:     DMPlexGetFaces_Internal(dm, cellDim, pStart[cellDepth], NULL, &faceSizeAll, NULL);
343:   }
344:   if (faceSizeAll > 4) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Do not support interpolation of meshes with faces of %D vertices", faceSizeAll);

346:   /* With hybrid grids, we first iterate on hybrid cells and start numbering the non-hybrid faces
347:      Then, faces for non-hybrid cells are numbered.
348:      This is to guarantee consistent orientations (all 0) of all the points in the cone of the hybrid cells */
349:   PetscHashIJKLCreate(&faceTable);
350:   for (outerloop = 0, face = pStart[faceDepth]; outerloop < 2; outerloop++) {
351:     PetscInt start, end;

353:     start = outerloop == 0 ? pMax : pStart[cellDepth];
354:     end = outerloop == 0 ? pEnd[cellDepth] : pMax;
355:     for (c = start; c < end; ++c) {
356:       const PetscInt *cellFaces;
357:       PetscInt        numCellFaces, faceSize, faceSizeInc, cf;

359:       if (c < pMax) {
360:         DMPlexGetFaces_Internal(dm, cellDim, c, &numCellFaces, &faceSize, &cellFaces);
361:         if (faceSize != faceSizeAll) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Inconsistent face for cell %D of size %D != %D", c, faceSize, faceSizeAll);
362:       } else { /* Hybrid cell */
363:         const PetscInt *cone;
364:         PetscInt        numCellFacesN, coneSize;

366:         DMPlexGetConeSize(dm, c, &coneSize);
367:         if (coneSize != coneSizeH) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unexpected hybrid coneSize %D != %D", coneSize, coneSizeH);
368:         DMPlexGetCone(dm, c, &cone);
369:         DMPlexGetRawFacesHybrid_Internal(dm, cellDim, coneSize, cone, &numCellFaces, &numCellFacesN, &faceSize, &cellFaces);
370:         if (numCellFaces != numCellFacesH) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unexpected numCellFaces %D != %D for hybrid cell %D", numCellFaces, numCellFacesH, c);
371:         faceSize = PetscMax(faceSize, -faceSize);
372:         if (faceSize > 4) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Do not support interpolation of meshes with faces of %D vertices", faceSize);
373:         numCellFaces = numCellFacesN; /* process only non-hybrid faces */
374:       }
375:       faceSizeInc = faceSize;
376:       for (cf = 0; cf < numCellFaces; ++cf) {
377:         const PetscInt   *cellFace = &cellFaces[cf*faceSizeInc];
378:         PetscInt          faceSizeH = faceSize;
379:         PetscHashIJKLKey  key;
380:         PetscHashIter     iter;
381:         PetscBool         missing;

383:         if (faceSizeInc == 2) {
384:           key.i = PetscMin(cellFace[0], cellFace[1]);
385:           key.j = PetscMax(cellFace[0], cellFace[1]);
386:           key.k = PETSC_MAX_INT;
387:           key.l = PETSC_MAX_INT;
388:         } else {
389:           key.i = cellFace[0];
390:           key.j = cellFace[1];
391:           key.k = cellFace[2];
392:           key.l = faceSize > 3 ? (cellFace[3] < 0 ? faceSizeH = 3, PETSC_MAX_INT : cellFace[3]) : PETSC_MAX_INT;
393:           PetscSortInt(faceSize, (PetscInt *) &key);
394:         }
395:         /* this check is redundant for non-hybrid meshes */
396:         if (faceSizeH != faceSizeAll) SETERRQ4(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unexpected number of vertices for face %D of point %D -> %D != %D", cf, c, faceSizeH, faceSizeAll);
397:         PetscHashIJKLPut(faceTable, key, &iter, &missing);
398:         if (missing) {PetscHashIJKLIterSet(faceTable, iter, face++);}
399:       }
400:       if (c < pMax) {
401:         DMPlexRestoreFaces_Internal(dm, cellDim, c, &numCellFaces, &faceSize, &cellFaces);
402:       } else {
403:         DMPlexRestoreRawFacesHybrid_Internal(dm, cellDim, coneSizeH, NULL, NULL, NULL, NULL, &cellFaces);
404:       }
405:     }
406:   }
407:   pEnd[faceDepth] = face;

409:   /* Second pass for hybrid meshes: number hybrid faces */
410:   for (c = pMax; c < pEnd[cellDepth]; ++c) {
411:     const PetscInt *cellFaces, *cone;
412:     PetscInt        numCellFaces, numCellFacesN, faceSize, cf, coneSize;

414:     DMPlexGetConeSize(dm, c, &coneSize);
415:     DMPlexGetCone(dm, c, &cone);
416:     DMPlexGetRawFacesHybrid_Internal(dm, cellDim, coneSize, cone, &numCellFaces, &numCellFacesN, &faceSize, &cellFaces);
417:     if (numCellFaces != numCellFacesH) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unexpected hybrid numCellFaces %D != %D", numCellFaces, numCellFacesH);
418:     faceSize = PetscMax(faceSize, -faceSize);
419:     for (cf = numCellFacesN; cf < numCellFaces; ++cf) { /* These are the hybrid faces */
420:       const PetscInt   *cellFace = &cellFaces[cf*faceSize];
421:       PetscHashIJKLKey  key;
422:       PetscHashIter     iter;
423:       PetscBool         missing;
424:       PetscInt          faceSizeH = faceSize;

426:       if (faceSize == 2) {
427:         key.i = PetscMin(cellFace[0], cellFace[1]);
428:         key.j = PetscMax(cellFace[0], cellFace[1]);
429:         key.k = PETSC_MAX_INT;
430:         key.l = PETSC_MAX_INT;
431:       } else {
432:         key.i = cellFace[0];
433:         key.j = cellFace[1];
434:         key.k = cellFace[2];
435:         key.l = faceSize > 3 ? (cellFace[3] < 0 ? faceSizeH = 3, PETSC_MAX_INT : cellFace[3]) : PETSC_MAX_INT;
436:         PetscSortInt(faceSize, (PetscInt *) &key);
437:       }
438:       if (faceSizeH != faceSizeAllH) SETERRQ4(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unexpected number of vertices for hybrid face %D of point %D -> %D != %D", cf, c, faceSizeH, faceSizeAllH);
439:       PetscHashIJKLPut(faceTable, key, &iter, &missing);
440:       if (missing) {PetscHashIJKLIterSet(faceTable, iter, face++);}
441:     }
442:     DMPlexRestoreRawFacesHybrid_Internal(dm, cellDim, coneSize, cone, &numCellFaces, &numCellFacesN, &faceSize, &cellFaces);
443:   }
444:   faceH = face - pEnd[faceDepth];
445:   if (faceH) {
446:     if (fMax == PETSC_DETERMINE) fMax = pEnd[faceDepth];
447:     else if (eMax == PETSC_DETERMINE) eMax = pEnd[faceDepth];
448:     else SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of unassigned hybrid facets %D for cellDim %D and dimension %D", faceH, cellDim, dim);
449:   }
450:   pEnd[faceDepth] = face;
451:   PetscHashIJKLDestroy(&faceTable);
452:   /* Count new points */
453:   for (d = 0; d <= depth; ++d) {
454:     numPoints += pEnd[d]-pStart[d];
455:   }
456:   DMPlexSetChart(idm, 0, numPoints);
457:   /* Set cone sizes */
458:   for (d = 0; d <= depth; ++d) {
459:     PetscInt coneSize, p;

461:     if (d == faceDepth) {
462:       /* I see no way to do this if we admit faces of different shapes */
463:       for (p = pStart[d]; p < pEnd[d]-faceH; ++p) {
464:         DMPlexSetConeSize(idm, p, faceSizeAll);
465:       }
466:       for (p = pEnd[d]-faceH; p < pEnd[d]; ++p) {
467:         DMPlexSetConeSize(idm, p, faceSizeAllH);
468:       }
469:     } else if (d == cellDepth) {
470:       for (p = pStart[d]; p < pEnd[d]; ++p) {
471:         /* Number of cell faces may be different from number of cell vertices*/
472:         if (p < pMax) {
473:           DMPlexGetFaces_Internal(dm, cellDim, p, &coneSize, NULL, NULL);
474:         } else {
475:           DMPlexGetFacesHybrid_Internal(dm, cellDim, p, &coneSize, NULL, NULL, NULL);
476:         }
477:         DMPlexSetConeSize(idm, p, coneSize);
478:       }
479:     } else {
480:       for (p = pStart[d]; p < pEnd[d]; ++p) {
481:         DMPlexGetConeSize(dm, p, &coneSize);
482:         DMPlexSetConeSize(idm, p, coneSize);
483:       }
484:     }
485:   }
486:   DMSetUp(idm);
487:   /* Get face cones from subsets of cell vertices */
488:   if (faceSizeAll > 4) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Do not support interpolation of meshes with faces of %D vertices", faceSizeAll);
489:   PetscHashIJKLCreate(&faceTable);
490:   for (d = depth; d > cellDepth; --d) {
491:     const PetscInt *cone;
492:     PetscInt        p;

494:     for (p = pStart[d]; p < pEnd[d]; ++p) {
495:       DMPlexGetCone(dm, p, &cone);
496:       DMPlexSetCone(idm, p, cone);
497:       DMPlexGetConeOrientation(dm, p, &cone);
498:       DMPlexSetConeOrientation(idm, p, cone);
499:     }
500:   }
501:   for (outerloop = 0, face = pStart[faceDepth]; outerloop < 2; outerloop++) {
502:     PetscInt start, end;

504:     start = outerloop == 0 ? pMax : pStart[cellDepth];
505:     end = outerloop == 0 ? pEnd[cellDepth] : pMax;
506:     for (c = start; c < end; ++c) {
507:       const PetscInt *cellFaces;
508:       PetscInt        numCellFaces, faceSize, faceSizeInc, cf;

510:       if (c < pMax) {
511:         DMPlexGetFaces_Internal(dm, cellDim, c, &numCellFaces, &faceSize, &cellFaces);
512:         if (faceSize != faceSizeAll) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Inconsistent face for cell %D of size %D != %D", c, faceSize, faceSizeAll);
513:       } else {
514:         const PetscInt *cone;
515:         PetscInt        numCellFacesN, coneSize;

517:         DMPlexGetConeSize(dm, c, &coneSize);
518:         if (coneSize != coneSizeH) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unexpected hybrid coneSize %D != %D", coneSize, coneSizeH);
519:         DMPlexGetCone(dm, c, &cone);
520:         DMPlexGetRawFacesHybrid_Internal(dm, cellDim, coneSize, cone, &numCellFaces, &numCellFacesN, &faceSize, &cellFaces);
521:         if (numCellFaces != numCellFacesH) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unexpected numCellFaces %D != %D for hybrid cell %D", numCellFaces, numCellFacesH, c);
522:         faceSize = PetscMax(faceSize, -faceSize);
523:         if (faceSize > 4) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Do not support interpolation of meshes with faces of %D vertices", faceSize);
524:         numCellFaces = numCellFacesN; /* process only non-hybrid faces */
525:       }
526:       faceSizeInc = faceSize;
527:       for (cf = 0; cf < numCellFaces; ++cf) {
528:         const PetscInt  *cellFace = &cellFaces[cf*faceSizeInc];
529:         PetscHashIJKLKey key;
530:         PetscHashIter    iter;
531:         PetscBool        missing;

533:         if (faceSizeInc == 2) {
534:           key.i = PetscMin(cellFace[0], cellFace[1]);
535:           key.j = PetscMax(cellFace[0], cellFace[1]);
536:           key.k = PETSC_MAX_INT;
537:           key.l = PETSC_MAX_INT;
538:         } else {
539:           key.i = cellFace[0];
540:           key.j = cellFace[1];
541:           key.k = cellFace[2];
542:           key.l = faceSizeInc > 3 ? (cellFace[3] < 0 ? faceSize = 3, PETSC_MAX_INT : cellFace[3]) : PETSC_MAX_INT;
543:           PetscSortInt(faceSizeInc, (PetscInt *) &key);
544:         }
545:         PetscHashIJKLPut(faceTable, key, &iter, &missing);
546:         if (missing) {
547:           DMPlexSetCone(idm, face, cellFace);
548:           PetscHashIJKLIterSet(faceTable, iter, face);
549:           DMPlexInsertCone(idm, c, cf, face++);
550:         } else {
551:           const PetscInt *cone;
552:           PetscInt        coneSize, ornt, i, j, f;

554:           PetscHashIJKLIterGet(faceTable, iter, &f);
555:           DMPlexInsertCone(idm, c, cf, f);
556:           /* Orient face: Do not allow reverse orientation at the first vertex */
557:           DMPlexGetConeSize(idm, f, &coneSize);
558:           DMPlexGetCone(idm, f, &cone);
559:           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);
560:           /* - First find the initial vertex */
561:           for (i = 0; i < faceSize; ++i) if (cellFace[0] == cone[i]) break;
562:           /* - Try forward comparison */
563:           for (j = 0; j < faceSize; ++j) if (cellFace[j] != cone[(i+j)%faceSize]) break;
564:           if (j == faceSize) {
565:             if ((faceSize == 2) && (i == 1)) ornt = -2;
566:             else                             ornt = i;
567:           } else {
568:             /* - Try backward comparison */
569:             for (j = 0; j < faceSize; ++j) if (cellFace[j] != cone[(i+faceSize-j)%faceSize]) break;
570:             if (j == faceSize) {
571:               if (i == 0) ornt = -faceSize;
572:               else        ornt = -i;
573:             } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Could not determine face orientation");
574:           }
575:           DMPlexInsertConeOrientation(idm, c, cf, ornt);
576:         }
577:       }
578:       if (c < pMax) {
579:         DMPlexRestoreFaces_Internal(dm, cellDim, c, &numCellFaces, &faceSize, &cellFaces);
580:       } else {
581:         DMPlexRestoreRawFacesHybrid_Internal(dm, cellDim, coneSizeH, NULL, NULL, NULL, NULL, &cellFaces);
582:       }
583:     }
584:   }
585:   /* Second pass for hybrid meshes: orient hybrid faces */
586:   for (c = pMax; c < pEnd[cellDepth]; ++c) {
587:     const PetscInt *cellFaces, *cone;
588:     PetscInt        numCellFaces, numCellFacesN, faceSize, cf, coneSize;

590:     DMPlexGetConeSize(dm, c, &coneSize);
591:     DMPlexGetCone(dm, c, &cone);
592:     DMPlexGetRawFacesHybrid_Internal(dm, cellDim, coneSize, cone, &numCellFaces, &numCellFacesN, &faceSize, &cellFaces);
593:     if (numCellFaces != numCellFacesH) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unexpected hybrid numCellFaces %D != %D", numCellFaces, numCellFacesH);
594:     faceSize = PetscMax(faceSize, -faceSize);
595:     for (cf = numCellFacesN; cf < numCellFaces; ++cf) { /* These are the hybrid faces */
596:       const PetscInt   *cellFace = &cellFaces[cf*faceSize];
597:       PetscHashIJKLKey key;
598:       PetscHashIter    iter;
599:       PetscBool        missing;
600:       PetscInt         faceSizeH = faceSize;

602:       if (faceSize == 2) {
603:         key.i = PetscMin(cellFace[0], cellFace[1]);
604:         key.j = PetscMax(cellFace[0], cellFace[1]);
605:         key.k = PETSC_MAX_INT;
606:         key.l = PETSC_MAX_INT;
607:       } else {
608:         key.i = cellFace[0];
609:         key.j = cellFace[1];
610:         key.k = cellFace[2];
611:         key.l = faceSize > 3 ? (cellFace[3] < 0 ? faceSizeH = 3, PETSC_MAX_INT : cellFace[3]) : PETSC_MAX_INT;
612:         PetscSortInt(faceSize, (PetscInt *) &key);
613:       }
614:       if (faceSizeH != faceSizeAllH) SETERRQ4(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unexpected number of vertices for hybrid face %D of point %D -> %D != %D", cf, c, faceSizeH, faceSizeAllH);
615:       PetscHashIJKLPut(faceTable, key, &iter, &missing);
616:       if (missing) {
617:         DMPlexSetCone(idm, face, cellFace);
618:         PetscHashIJKLIterSet(faceTable, iter, face);
619:         DMPlexInsertCone(idm, c, cf, face++);
620:       } else {
621:         const PetscInt *cone;
622:         PetscInt        coneSize, ornt, i, j, f;

624:         PetscHashIJKLIterGet(faceTable, iter, &f);
625:         DMPlexInsertCone(idm, c, cf, f);
626:         /* Orient face: Do not allow reverse orientation at the first vertex */
627:         DMPlexGetConeSize(idm, f, &coneSize);
628:         DMPlexGetCone(idm, f, &cone);
629:         if (coneSize != faceSizeH) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid number of face vertices %D for face %D should be %D", coneSize, f, faceSizeH);
630:         /* - First find the initial vertex */
631:         for (i = 0; i < faceSizeH; ++i) if (cellFace[0] == cone[i]) break;
632:         /* - Try forward comparison */
633:         for (j = 0; j < faceSizeH; ++j) if (cellFace[j] != cone[(i+j)%faceSizeH]) break;
634:         if (j == faceSizeH) {
635:           if ((faceSizeH == 2) && (i == 1)) ornt = -2;
636:           else                             ornt = i;
637:         } else {
638:           /* - Try backward comparison */
639:           for (j = 0; j < faceSizeH; ++j) if (cellFace[j] != cone[(i+faceSizeH-j)%faceSizeH]) break;
640:           if (j == faceSizeH) {
641:             if (i == 0) ornt = -faceSizeH;
642:             else        ornt = -i;
643:           } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Could not determine face orientation");
644:         }
645:         DMPlexInsertConeOrientation(idm, c, cf, ornt);
646:       }
647:     }
648:     DMPlexRestoreRawFacesHybrid_Internal(dm, cellDim, coneSize, cone, &numCellFaces, &numCellFacesN, &faceSize, &cellFaces);
649:   }
650:   if (face != pEnd[faceDepth]) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid number of faces %D should be %D", face-pStart[faceDepth], pEnd[faceDepth]-pStart[faceDepth]);
651:   PetscFree2(pStart,pEnd);
652:   PetscHashIJKLDestroy(&faceTable);
653:   PetscFree2(pStart,pEnd);
654:   DMPlexSetHybridBounds(idm, cMax, fMax, eMax, vMax);
655:   DMPlexSymmetrize(idm);
656:   DMPlexStratify(idm);
657:   return(0);
658: }

660: PetscErrorCode DMPlexOrientCell(DM dm, PetscInt p, PetscInt masterConeSize, const PetscInt masterCone[])
661: {
662:   PetscInt coneSize;
663:   PetscInt start1=0;
664:   PetscBool reverse1=PETSC_FALSE;
665:   const PetscInt *cone=NULL;

669:   DMPlexGetConeSize(dm, p, &coneSize);
670:   if (!coneSize) return(0); /* do nothing for points with no cone */
671:   DMPlexGetCone(dm, p, &cone);
672:   DMPlexFixFaceOrientations_Orient_Private(coneSize, masterConeSize, masterCone, cone, &start1, &reverse1);
673: #if defined(PETSC_USE_DEBUG)
674:   if (PetscUnlikely(cone[start1] != masterCone[0])) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "The algorithm above is wrong as cone[%d] = %d != %d = masterCone[0]", start1, cone[start1], masterCone[0]);
675: #endif
676:   DMPlexOrientCell_Internal(dm, p, start1, reverse1);
677: #if defined(PETSC_USE_DEBUG)
678:   {
679:     PetscInt c;
680:     DMPlexGetCone(dm, p, &cone);
681:     for (c = 0; c < 2; c++) {
682:       if (PetscUnlikely(cone[c] != masterCone[c])) SETERRQ4(PETSC_COMM_SELF, PETSC_ERR_PLIB, "The algorithm above is wrong as cone[%d] = %d != %d = masterCone[%d]", c, cone[c], masterCone[c], c);
683:     }
684:   }
685: #endif
686:   return(0);
687: }

689: PetscErrorCode DMPlexOrientCell_Internal(DM dm, PetscInt p, PetscInt start1, PetscBool reverse1)
690: {
691:   PetscInt i, j, k, maxConeSize, coneSize, coneConeSize, supportSize, supportConeSize;
692:   PetscInt start0, start;
693:   PetscBool reverse0, reverse;
694:   PetscInt newornt;
695:   const PetscInt *cone=NULL, *support=NULL, *supportCone=NULL, *ornts=NULL;
696:   PetscInt *newcone=NULL, *newornts=NULL;

700:   if (!start1 && !reverse1) return(0);
701:   DMPlexGetConeSize(dm, p, &coneSize);
702:   if (!coneSize) return(0); /* do nothing for points with no cone */
703:   DMPlexGetCone(dm, p, &cone);
704:   DMPlexGetMaxSizes(dm, &maxConeSize, NULL);
705:   /* permute p's cone and orientations */
706:   DMPlexGetConeOrientation(dm, p, &ornts);
707:   DMGetWorkArray(dm, maxConeSize, MPIU_INT, &newcone);
708:   DMGetWorkArray(dm, maxConeSize, MPIU_INT, &newornts);
709:   DMPlexFixFaceOrientations_Permute_Private(coneSize, cone, start1, reverse1, newcone);
710:   DMPlexFixFaceOrientations_Permute_Private(coneSize, ornts, start1, reverse1, newornts);
711:   /* if direction of p (face) is flipped, flip also p's cone points (edges) */
712:   if (reverse1) {
713:     for (i=0; i<coneSize; i++) {
714:       DMPlexGetConeSize(dm, cone[i], &coneConeSize);
715:       DMPlexFixFaceOrientations_Translate_Private(newornts[i], &start0, &reverse0);
716:       DMPlexFixFaceOrientations_Combine_Private(coneConeSize, start0, reverse0, 1, PETSC_FALSE, &start, &reverse);
717:       DMPlexFixFaceOrientations_TranslateBack_Private(coneConeSize, start, reverse, &newornts[i]);
718:     }
719:   }
720:   DMPlexSetConeOrientation(dm, p, newornts);
721:   /* fix oriention of p within cones of p's support points */
722:   DMPlexGetSupport(dm, p, &support);
723:   DMPlexGetSupportSize(dm, p, &supportSize);
724:   for (j=0; j<supportSize; j++) {
725:     DMPlexGetCone(dm, support[j], &supportCone);
726:     DMPlexGetConeSize(dm, support[j], &supportConeSize);
727:     DMPlexGetConeOrientation(dm, support[j], &ornts);
728:     for (k=0; k<supportConeSize; k++) {
729:       if (supportCone[k] != p) continue;
730:       DMPlexFixFaceOrientations_Translate_Private(ornts[k], &start0, &reverse0);
731:       DMPlexFixFaceOrientations_Combine_Private(coneSize, start0, reverse0, start1, reverse1, &start, &reverse);
732:       DMPlexFixFaceOrientations_TranslateBack_Private(coneSize, start, reverse, &newornt);
733:       DMPlexInsertConeOrientation(dm, support[j], k, newornt);
734:     }
735:   }
736:   /* rewrite cone */
737:   DMPlexSetCone(dm, p, newcone);
738:   DMRestoreWorkArray(dm, maxConeSize, MPIU_INT, &newcone);
739:   DMRestoreWorkArray(dm, maxConeSize, MPIU_INT, &newornts);
740:   return(0);
741: }

743: static PetscErrorCode SortRmineRremoteByRemote_Private(PetscSF sf, PetscInt *rmine1[], PetscInt *rremote1[])
744: {
745:   PetscInt            nleaves;
746:   PetscInt            nranks;
747:   const PetscMPIInt  *ranks=NULL;
748:   const PetscInt     *roffset=NULL, *rmine=NULL, *rremote=NULL;
749:   PetscInt            n, o, r;
750:   PetscErrorCode      ierr;

753:   PetscSFGetRanks(sf, &nranks, &ranks, &roffset, &rmine, &rremote);
754:   nleaves = roffset[nranks];
755:   PetscMalloc2(nleaves, rmine1, nleaves, rremote1);
756:   for (r=0; r<nranks; r++) {
757:     /* simultaneously sort rank-wise portions of rmine & rremote by values in rremote
758:        - to unify order with the other side */
759:     o = roffset[r];
760:     n = roffset[r+1] - o;
761:     PetscMemcpy(&(*rmine1)[o], &rmine[o], n*sizeof(PetscInt));
762:     PetscMemcpy(&(*rremote1)[o], &rremote[o], n*sizeof(PetscInt));
763:     PetscSortIntWithArray(n, &(*rremote1)[o], &(*rmine1)[o]);
764:   }
765:   return(0);
766: }

768: PetscErrorCode DMPlexOrientInterface(DM dm)
769: {
770:   PetscSF           sf=NULL;
771:   PetscInt          (*roots)[2], (*leaves)[2];
772:   PetscMPIInt       (*rootsRanks)[2], (*leavesRanks)[2];
773:   const PetscInt    *locals=NULL;
774:   const PetscSFNode *remotes=NULL;
775:   PetscInt           nroots, nleaves, p, c;
776:   PetscInt           nranks, n, o, r;
777:   const PetscMPIInt *ranks=NULL;
778:   const PetscInt    *roffset=NULL;
779:   PetscInt          *rmine1=NULL, *rremote1=NULL; /* rmine and rremote copies simultaneously sorted by rank and rremote */
780:   const PetscInt    *cone=NULL;
781:   PetscInt           coneSize, ind0;
782:   MPI_Comm           comm;
783:   PetscMPIInt        rank;
784:   PetscInt           debug = 0;
785:   PetscErrorCode     ierr;

788:   DMGetPointSF(dm, &sf);
789:   PetscSFGetGraph(sf, &nroots, &nleaves, &locals, &remotes);
790:   if (nroots < 0) return(0);
791:   PetscSFSetUp(sf);
792:   PetscSFGetRanks(sf, &nranks, &ranks, &roffset, NULL, NULL);
793: #if defined(PETSC_USE_DEBUG)
794:   DMViewFromOptions(dm, NULL, "-before_fix_dm_view");
795:   DMPlexCheckPointSF(dm);
796: #endif
797:   SortRmineRremoteByRemote_Private(sf, &rmine1, &rremote1);
798:   PetscMalloc4(nroots, &roots, nroots, &leaves, nroots, &rootsRanks, nroots, &leavesRanks);
799:   PetscObjectGetComm((PetscObject) dm, &comm);
800:   MPI_Comm_rank(comm, &rank);
801:   if (debug && rank == 0) {PetscSynchronizedPrintf(comm, "Roots\n");}
802:   for (p = 0; p < nroots; ++p) {
803:     DMPlexGetConeSize(dm, p, &coneSize);
804:     DMPlexGetCone(dm, p, &cone);
805:     /* Translate all points to root numbering */
806:     for (c = 0; c < 2; c++) {
807:       if (coneSize > 1) {
808:         PetscFindInt(cone[c], nleaves, locals, &ind0);
809:         if (ind0 < 0) {
810:           roots[p][c] = cone[c];
811:           rootsRanks[p][c] = rank;
812:         } else {
813:           roots[p][c] = remotes[ind0].index;
814:           rootsRanks[p][c] = remotes[ind0].rank;
815:         }
816:       } else {
817:         roots[p][c] = -1;
818:         rootsRanks[p][c] = -1;
819:       }
820:     }
821:   }
822:   if (debug) {
823:     for (p = 0; p < nroots; ++p) {
824:       DMPlexGetConeSize(dm, p, &coneSize);
825:       DMPlexGetCone(dm, p, &cone);
826:       if (coneSize > 1) {
827:         PetscSynchronizedPrintf(comm, "[%d]  %D: cone=[%D %D] roots=[%D %D] rootsRanks=[%D %D]\n", rank, p, cone[0], cone[1], roots[p][0], roots[p][1], rootsRanks[p][0], rootsRanks[p][1]);
828:       }
829:     }
830:     PetscSynchronizedFlush(comm, NULL);
831:   }
832:   for (p = 0; p < nroots; ++p) {
833:     for (c = 0; c < 2; c++) {
834:       leaves[p][c] = -2;
835:       leavesRanks[p][c] = -2;
836:     }
837:   }
838:   PetscSFBcastBegin(sf, MPIU_2INT, roots, leaves);
839:   PetscSFBcastBegin(sf, MPI_2INT, rootsRanks, leavesRanks);
840:   PetscSFBcastEnd(sf, MPIU_2INT, roots, leaves);
841:   PetscSFBcastEnd(sf, MPI_2INT, rootsRanks, leavesRanks);
842:   if (debug) {PetscSynchronizedFlush(comm, NULL);}
843:   if (debug && rank == 0) {PetscSynchronizedPrintf(comm, "Referred leaves\n");}
844:   for (p = 0; p < nroots; ++p) {
845:     if (leaves[p][0] < 0) continue;
846:     DMPlexGetConeSize(dm, p, &coneSize);
847:     DMPlexGetCone(dm, p, &cone);
848:     if (debug) {PetscSynchronizedPrintf(comm, "[%d]  %D: cone=[%D %D] leaves=[%D %D] roots=[%D %D] leavesRanks=[%D %D] rootsRanks=[%D %D]\n", rank, p, cone[0], cone[1], leaves[p][0], leaves[p][1], roots[p][0], roots[p][1], leavesRanks[p][0], leavesRanks[p][1], rootsRanks[p][0], rootsRanks[p][1]);}
849:     if ((leaves[p][0] != roots[p][0]) || (leaves[p][1] != roots[p][1]) || (leavesRanks[p][0] != rootsRanks[p][0]) || (leavesRanks[p][0] != rootsRanks[p][0])) {
850:       PetscInt masterCone[2];
851:       /* Translate these two cone points back to leave numbering */
852:       for (c = 0; c < 2; c++) {
853:         if (leavesRanks[p][c] == rank) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "this should never happen - remote rank of point %D is the same rank",leavesRanks[p][c]);
854:         /* Find index of rank leavesRanks[p][c] among remote ranks */
855:         /* No need for PetscMPIIntCast because these integers were originally cast from PetscMPIInt. */
856:         PetscFindMPIInt((PetscMPIInt)leavesRanks[p][c], nranks, ranks, &r);
857:         if (PetscUnlikely(r < 0)) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "this should never happen - rank %D not found among remote ranks",leavesRanks[p][c]);
858:         /* Find point leaves[p][c] among remote points aimed at rank leavesRanks[p][c] */
859:         o = roffset[r];
860:         n = roffset[r+1] - o;
861:         PetscFindInt(leaves[p][c], n, &rremote1[o], &ind0);
862:         if (PetscUnlikely(ind0 < 0)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "No cone point of %D is connected to (%D, %D) - it seems there is missing connection in point SF!",p,ranks[r],leaves[p][c]);
863:         /* Get the corresponding local point */
864:         masterCone[c] = rmine1[o+ind0];
865:       }
866:       if (debug) {PetscSynchronizedPrintf(comm, "[%d]  %D: masterCone=[%D %D]\n", rank, p, masterCone[0], masterCone[1]);}
867:       /* Vaclav's note: Here we only compare first 2 points of the cone. Full cone size would lead to stronger self-checking. */
868:       DMPlexOrientCell(dm, p, 2, masterCone);
869:     }
870:   }
871: #if defined(PETSC_USE_DEBUG)
872:   DMViewFromOptions(dm, NULL, "-after_fix_dm_view");
873:   for (r = 0; r < nleaves; ++r) {
874:     p = locals[r];
875:     DMPlexGetConeSize(dm, p, &coneSize);
876:     if (!coneSize) continue;
877:     DMPlexGetCone(dm, p, &cone);
878:     for (c = 0; c < 2; c++) {
879:       PetscFindInt(cone[c], nleaves, locals, &ind0);
880:       if (ind0 < 0) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Point SF contains %D but is missing its cone point cone[%D] = %D!", p, c, cone[c]);
881:       if (leaves[p][c] != remotes[ind0].index || leavesRanks[p][c] != remotes[ind0].rank) {
882:         if (leavesRanks[p][c] == rank) {
883:           PetscInt ind1;
884:           PetscFindInt(leaves[p][c], nleaves, locals, &ind1);
885:           if (ind1 < 0) {
886:             SETERRQ8(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Point %D = locals[%d]: cone[%D]=%D --> (%D, %D) differs from the enforced (%D, %D). The latter was not even found among the local SF points - it is probably broken!", p, r, c, cone[c], remotes[ind0].rank, remotes[ind0].index, leavesRanks[p][c], leaves[p][c]);
887:           } else {
888:             SETERRQ9(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Point %D = locals[%d]: cone[%D]=%D --> (%D, %D) differs from the enforced %D --> (%D, %D). Is the algorithm above or the point SF broken?", p, r, c, cone[c], remotes[ind0].rank, remotes[ind0].index, leaves[p][c], remotes[ind1].rank, remotes[ind1].index);
889:           }
890:         } else {
891:           SETERRQ8(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Point %D = locals[%d]: cone[%D]=%D --> (%D, %D) differs from the enforced (%D, %D). Is the algorithm above or the point SF broken?", p, r, c, cone[c], remotes[ind0].rank, remotes[ind0].index, leavesRanks[p][c], leaves[p][c]);
892:         }
893:       }
894:     }
895:   }
896: #endif
897:   if (debug) {PetscSynchronizedFlush(comm, NULL);}
898:   PetscFree4(roots, leaves, rootsRanks, leavesRanks);
899:   PetscFree2(rmine1, rremote1);
900:   return(0);
901: }

903: static PetscErrorCode IntArrayViewFromOptions(MPI_Comm comm, const char opt[], const char name[], const char idxname[], const char valname[], PetscInt n, const PetscInt a[])
904: {
905:   PetscInt       idx;
906:   PetscMPIInt    rank;
907:   PetscBool      flg;

911:   PetscOptionsHasName(NULL, NULL, opt, &flg);
912:   if (!flg) return(0);
913:   MPI_Comm_rank(comm, &rank);
914:   PetscSynchronizedPrintf(comm, "[%d]%s:\n", rank, name);
915:   for (idx = 0; idx < n; ++idx) {PetscSynchronizedPrintf(comm, "[%d]%s %D %s %D\n", rank, idxname, idx, valname, a[idx]);}
916:   PetscSynchronizedFlush(comm, NULL);
917:   return(0);
918: }

920: static PetscErrorCode SFNodeArrayViewFromOptions(MPI_Comm comm, const char opt[], const char name[], const char idxname[], PetscInt n, const PetscSFNode a[])
921: {
922:   PetscInt       idx;
923:   PetscMPIInt    rank;
924:   PetscBool      flg;

928:   PetscOptionsHasName(NULL, NULL, opt, &flg);
929:   if (!flg) return(0);
930:   MPI_Comm_rank(comm, &rank);
931:   PetscSynchronizedPrintf(comm, "[%d]%s:\n", rank, name);
932:   if (idxname) {
933:     for (idx = 0; idx < n; ++idx) {PetscSynchronizedPrintf(comm, "[%d]%s %D rank %D index %D\n", rank, idxname, idx, a[idx].rank, a[idx].index);}
934:   } else {
935:     for (idx = 0; idx < n; ++idx) {PetscSynchronizedPrintf(comm, "[%d]rank %D index %D\n", rank, a[idx].rank, a[idx].index);}
936:   }
937:   PetscSynchronizedFlush(comm, NULL);
938:   return(0);
939: }

941: static PetscErrorCode DMPlexMapToLocalPoint(PetscHMapIJ roothash, const PetscInt localPoints[], PetscMPIInt rank, PetscSFNode remotePoint, PetscInt *localPoint)
942: {

946:   if (remotePoint.rank == rank) {
947:     *localPoint = remotePoint.index;
948:   } else {
949:     PetscHashIJKey key;
950:     PetscInt       root;

952:     key.i = remotePoint.index;
953:     key.j = remotePoint.rank;
954:     PetscHMapIJGet(roothash, key, &root);
955:     if (root >= 0) {
956:       *localPoint = localPoints[root];
957:     } else PetscFunctionReturn(1);
958:   }
959:   return(0);
960: }

962: /*@
963:   DMPlexInterpolatePointSF - Insert interpolated points in the overlap into the PointSF in parallel, following local interpolation

965:   Collective on DM

967:   Input Parameters:
968: + dm      - The interpolated DM
969: - pointSF - The initial SF without interpolated points

971:   Output Parameter:
972: . pointSF - The SF including interpolated points

974:   Level: intermediate

976:    Note: All debugging for this process can be turned on with the options: -dm_interp_pre_view -petscsf_interp_pre_view -petscsection_interp_candidate_view -petscsection_interp_candidate_remote_view -petscsection_interp_claim_view -petscsf_interp_pre_view -dmplex_interp_debug

978: .keywords: mesh
979: .seealso: DMPlexInterpolate(), DMPlexUninterpolate()
980: @*/
981: PetscErrorCode DMPlexInterpolatePointSF(DM dm, PetscSF pointSF)
982: {
983:   /*
984:        Okay, the algorithm is:
985:          - Take each point in the overlap (root)
986:          - Look at the neighboring points in the overlap (candidates)
987:          - Send these candidate points to neighbors
988:          - Neighbor checks for edge between root and candidate
989:          - If edge is found, it replaces candidate point with edge point
990:          - Send back the overwritten candidates (claims)
991:          - Original guy checks for edges, different from original candidate, and gets its own edge
992:          - This pair is put into SF

994:        We need a new algorithm that tolerates groups larger than 2.
995:          - Take each point in the overlap (root)
996:          - Find all collections of points in the overlap which make faces (do early join)
997:          - Send collections as candidates (add size as first number)
998:            - Make sure to send collection to all owners of all overlap points in collection
999:          - Neighbor check for face in collections
1000:          - If face is found, it replaces candidate point with face point
1001:          - Send back the overwritten candidates (claims)
1002:          - Original guy checks for faces, different from original candidate, and gets its own face
1003:          - This pair is put into SF
1004:   */
1005:   PetscHMapI         leafhash;
1006:   PetscHMapIJ        roothash;
1007:   const PetscInt    *localPoints, *rootdegree;
1008:   const PetscSFNode *remotePoints;
1009:   PetscSFNode       *candidates, *candidatesRemote, *claims;
1010:   PetscSection       candidateSection, candidateSectionRemote, claimSection;
1011:   PetscInt           numLeaves, l, numRoots, r, candidatesSize, candidatesRemoteSize;
1012:   PetscMPIInt        size, rank;
1013:   PetscHashIJKey     key;
1014:   PetscBool          debug = PETSC_FALSE;
1015:   PetscErrorCode     ierr;

1018:   PetscOptionsHasName(NULL, ((PetscObject) dm)->prefix, "-dmplex_interp_debug", &debug);
1019:   MPI_Comm_size(PetscObjectComm((PetscObject) dm), &size);
1020:   MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank);
1021:   PetscSFGetGraph(pointSF, &numRoots, &numLeaves, &localPoints, &remotePoints);
1022:   if (size < 2 || numRoots < 0) return(0);
1023:   PetscObjectViewFromOptions((PetscObject) dm, NULL, "-dm_interp_pre_view");
1024:   PetscObjectViewFromOptions((PetscObject) pointSF, NULL, "-petscsf_interp_pre_view");
1025:   PetscLogEventBegin(DMPLEX_InterpolateSF,dm,0,0,0);
1026:   /* Build hashes of points in the SF for efficient lookup */
1027:   PetscHMapICreate(&leafhash);
1028:   PetscHMapIJCreate(&roothash);
1029:   for (l = 0; l < numLeaves; ++l) {
1030:     key.i = remotePoints[l].index;
1031:     key.j = remotePoints[l].rank;
1032:     PetscHMapISet(leafhash, localPoints[l], l);
1033:     PetscHMapIJSet(roothash, key, l);
1034:   }
1035:   /* Compute root degree to identify shared points */
1036:   PetscSFComputeDegreeBegin(pointSF, &rootdegree);
1037:   PetscSFComputeDegreeEnd(pointSF, &rootdegree);
1038:   IntArrayViewFromOptions(PetscObjectComm((PetscObject) dm), "-interp_root_degree_view", "Root degree", "point", "degree", numRoots, rootdegree);
1039:   /* Build a section / SFNode array of candidate points (face bd points) in the cone(support(leaf)),
1040:      where each candidate is defined by a set of remote points (roots) for the other points that define the face. */
1041:   PetscSectionCreate(PetscObjectComm((PetscObject) dm), &candidateSection);
1042:   PetscSectionSetChart(candidateSection, 0, numRoots);
1043:   {
1044:     PetscHMapIJ facehash;

1046:     PetscHMapIJCreate(&facehash);
1047:     for (l = 0; l < numLeaves; ++l) {
1048:       const PetscInt    localPoint = localPoints[l];
1049:       const PetscInt   *support;
1050:       PetscInt          supportSize, s;

1052:       if (debug) {PetscSynchronizedPrintf(PetscObjectComm((PetscObject) dm), "[%d]  Checking local point %D\n", rank, localPoint);}
1053:       DMPlexGetSupportSize(dm, localPoint, &supportSize);
1054:       DMPlexGetSupport(dm, localPoint, &support);
1055:       for (s = 0; s < supportSize; ++s) {
1056:         const PetscInt  face = support[s];
1057:         const PetscInt *cone;
1058:         PetscInt        coneSize, c, f, root;
1059:         PetscBool       isFace = PETSC_TRUE;

1061:         /* Only add face once */
1062:         if (debug) {PetscSynchronizedPrintf(PetscObjectComm((PetscObject) dm), "[%d]    Support point %D\n", rank, face);}
1063:         key.i = localPoint;
1064:         key.j = face;
1065:         PetscHMapIJGet(facehash, key, &f);
1066:         if (f >= 0) continue;
1067:         DMPlexGetConeSize(dm, face, &coneSize);
1068:         DMPlexGetCone(dm, face, &cone);
1069:         /* If a cone point does not map to leaves on any proc, then do not put face in SF */
1070:         for (c = 0; c < coneSize; ++c) {
1071:           if (debug) {PetscSynchronizedPrintf(PetscObjectComm((PetscObject) dm), "[%d]      Cone point %D\n", rank, cone[c]);}
1072:           PetscHMapIGet(leafhash, cone[c], &root);
1073:           if (!rootdegree[cone[c]] && (root < 0)) {isFace = PETSC_FALSE; break;}
1074:         }
1075:         if (isFace) {
1076:           if (debug) {PetscSynchronizedPrintf(PetscObjectComm((PetscObject) dm), "[%d]    Found shared face %D\n", rank, face);}
1077:           PetscHMapIJSet(facehash, key, l);
1078:           PetscSectionAddDof(candidateSection, localPoint, coneSize);
1079:         }
1080:       }
1081:     }
1082:     if (debug) {PetscSynchronizedFlush(PetscObjectComm((PetscObject) dm), NULL);}
1083:     PetscHMapIJClear(facehash);
1084:     PetscSectionSetUp(candidateSection);
1085:     PetscSectionGetStorageSize(candidateSection, &candidatesSize);
1086:     PetscMalloc1(candidatesSize, &candidates);
1087:     for (l = 0; l < numLeaves; ++l) {
1088:       const PetscInt    localPoint = localPoints[l];
1089:       const PetscInt   *support;
1090:       PetscInt          supportSize, s, offset, idx = 0;

1092:       if (debug) {PetscSynchronizedPrintf(PetscObjectComm((PetscObject) dm), "[%d]  Checking local point %D\n", rank, localPoint);}
1093:       PetscSectionGetOffset(candidateSection, localPoint, &offset);
1094:       DMPlexGetSupportSize(dm, localPoint, &supportSize);
1095:       DMPlexGetSupport(dm, localPoint, &support);
1096:       for (s = 0; s < supportSize; ++s) {
1097:         const PetscInt  face = support[s];
1098:         const PetscInt *cone;
1099:         PetscInt        coneSize, c, f, root;
1100:         PetscBool       isFace = PETSC_TRUE;

1102:         /* Only add face once */
1103:         if (debug) {PetscSynchronizedPrintf(PetscObjectComm((PetscObject) dm), "[%d]    Support point %D\n", rank, face);}
1104:         key.i = localPoint;
1105:         key.j = face;
1106:         PetscHMapIJGet(facehash, key, &f);
1107:         if (f >= 0) continue;
1108:         DMPlexGetConeSize(dm, face, &coneSize);
1109:         DMPlexGetCone(dm, face, &cone);
1110:         /* If a cone point does not map to leaves on any proc, then do not put face in SF */
1111:         for (c = 0; c < coneSize; ++c) {
1112:           if (debug) {PetscSynchronizedPrintf(PetscObjectComm((PetscObject) dm), "[%d]      Cone point %D\n", rank, cone[c]);}
1113:           PetscHMapIGet(leafhash, cone[c], &root);
1114:           if (!rootdegree[cone[c]] && (root < 0)) {isFace = PETSC_FALSE; break;}
1115:         }
1116:         if (isFace) {
1117:           if (debug) {PetscSynchronizedPrintf(PetscObjectComm((PetscObject) dm), "[%d]    Adding shared face %D at idx %D\n", rank, face, idx);}
1118:           PetscHMapIJSet(facehash, key, l);
1119:           candidates[offset+idx].rank    = -1;
1120:           candidates[offset+idx++].index = coneSize-1;
1121:           for (c = 0; c < coneSize; ++c) {
1122:             if (cone[c] == localPoint) continue;
1123:             if (rootdegree[cone[c]]) {
1124:               candidates[offset+idx].rank    = rank;
1125:               candidates[offset+idx++].index = cone[c];
1126:             } else {
1127:               PetscHMapIGet(leafhash, cone[c], &root);
1128:               if (root < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Cannot locate local point %D in SF", cone[c]);
1129:               candidates[offset+idx++] = remotePoints[root];
1130:             }
1131:           }
1132:         }
1133:       }
1134:     }
1135:     if (debug) {PetscSynchronizedFlush(PetscObjectComm((PetscObject) dm), NULL);}
1136:     PetscHMapIJDestroy(&facehash);
1137:     PetscObjectViewFromOptions((PetscObject) candidateSection, NULL, "-petscsection_interp_candidate_view");
1138:     SFNodeArrayViewFromOptions(PetscObjectComm((PetscObject) dm), "-petscsection_interp_candidate_view", "Candidates", NULL, candidatesSize, candidates);
1139:   }
1140:   /* Gather candidate section / array pair into the root partition via inverse(multi(pointSF)). */
1141:   /*   Note that this section is indexed by offsets into leaves, not by point number */
1142:   {
1143:     PetscSF   sfMulti, sfInverse, sfCandidates;
1144:     PetscInt *remoteOffsets;

1146:     PetscSFGetMultiSF(pointSF, &sfMulti);
1147:     PetscSFCreateInverseSF(sfMulti, &sfInverse);
1148:     PetscSectionCreate(PetscObjectComm((PetscObject) dm), &candidateSectionRemote);
1149:     PetscSFDistributeSection(sfInverse, candidateSection, &remoteOffsets, candidateSectionRemote);
1150:     PetscSFCreateSectionSF(sfInverse, candidateSection, remoteOffsets, candidateSectionRemote, &sfCandidates);
1151:     PetscSectionGetStorageSize(candidateSectionRemote, &candidatesRemoteSize);
1152:     PetscMalloc1(candidatesRemoteSize, &candidatesRemote);
1153:     PetscSFBcastBegin(sfCandidates, MPIU_2INT, candidates, candidatesRemote);
1154:     PetscSFBcastEnd(sfCandidates, MPIU_2INT, candidates, candidatesRemote);
1155:     PetscSFDestroy(&sfInverse);
1156:     PetscSFDestroy(&sfCandidates);
1157:     PetscFree(remoteOffsets);

1159:     PetscObjectViewFromOptions((PetscObject) candidateSectionRemote, NULL, "-petscsection_interp_candidate_remote_view");
1160:     SFNodeArrayViewFromOptions(PetscObjectComm((PetscObject) dm), "-petscsection_interp_candidate_remote_view", "Remote Candidates", NULL, candidatesRemoteSize, candidatesRemote);
1161:   }
1162:   /* */
1163:   {
1164:     PetscInt idx;
1165:     /* There is a section point for every leaf attached to a given root point */
1166:     for (r = 0, idx = 0; r < numRoots; ++r) {
1167:       PetscInt deg;
1168:       for (deg = 0; deg < rootdegree[r]; ++deg, ++idx) {
1169:         PetscInt offset, dof, d;

1171:         PetscSectionGetDof(candidateSectionRemote, idx, &dof);
1172:         PetscSectionGetOffset(candidateSectionRemote, idx, &offset);
1173:         for (d = 0; d < dof; ++d) {
1174:           const PetscInt  sizeInd   = offset+d;
1175:           const PetscInt  numPoints = candidatesRemote[sizeInd].index;
1176:           const PetscInt *join      = NULL;
1177:           PetscInt        points[1024], p, joinSize;

1179:           points[0] = r;
1180:           for (p = 0; p < numPoints; ++p) {
1181:             DMPlexMapToLocalPoint(roothash, localPoints, rank, candidatesRemote[offset+(++d)], &points[p+1]);
1182:             if (ierr) {d += numPoints-1 - p; break;} /* We got a point not in our overlap */
1183:             if (debug) {PetscSynchronizedPrintf(PetscObjectComm((PetscObject) dm), "[%d]  Checking local candidate %D\n", rank, points[p+1]);}
1184:           }
1185:           if (ierr) continue;
1186:           DMPlexGetJoin(dm, numPoints+1, points, &joinSize, &join);
1187:           if (joinSize == 1) {
1188:             if (debug) {PetscSynchronizedPrintf(PetscObjectComm((PetscObject) dm), "[%d]    Adding face %D at idx %D\n", rank, join[0], sizeInd);}
1189:             candidatesRemote[sizeInd].rank  = rank;
1190:             candidatesRemote[sizeInd].index = join[0];
1191:           }
1192:           DMPlexRestoreJoin(dm, numPoints+1, points, &joinSize, &join);
1193:         }
1194:       }
1195:     }
1196:     if (debug) {PetscSynchronizedFlush(PetscObjectComm((PetscObject) dm), NULL);}
1197:   }
1198:   /* Push claims back to receiver via the MultiSF and derive new pointSF mapping on receiver */
1199:   {
1200:     PetscSF         sfMulti, sfClaims, sfPointNew;
1201:     PetscSFNode    *remotePointsNew;
1202:     PetscHMapI      claimshash;
1203:     PetscInt       *remoteOffsets, *localPointsNew;
1204:     PetscInt        claimsSize, pStart, pEnd, root, numLocalNew, p, d;

1206:     PetscSFGetMultiSF(pointSF, &sfMulti);
1207:     PetscSectionCreate(PetscObjectComm((PetscObject) dm), &claimSection);
1208:     PetscSFDistributeSection(sfMulti, candidateSectionRemote, &remoteOffsets, claimSection);
1209:     PetscSFCreateSectionSF(sfMulti, candidateSectionRemote, remoteOffsets, claimSection, &sfClaims);
1210:     PetscSectionGetStorageSize(claimSection, &claimsSize);
1211:     PetscMalloc1(claimsSize, &claims);
1212:     PetscSFBcastBegin(sfClaims, MPIU_2INT, candidatesRemote, claims);
1213:     PetscSFBcastEnd(sfClaims, MPIU_2INT, candidatesRemote, claims);
1214:     PetscSFDestroy(&sfClaims);
1215:     PetscFree(remoteOffsets);
1216:     PetscObjectViewFromOptions((PetscObject) claimSection, NULL, "-petscsection_interp_claim_view");
1217:     SFNodeArrayViewFromOptions(PetscObjectComm((PetscObject) dm), "-petscsection_interp_claim_view", "Claims", NULL, claimsSize, claims);
1218:     /* Walk the original section of local supports and add an SF entry for each updated item */
1219:     PetscHMapICreate(&claimshash);
1220:     for (p = 0; p < numRoots; ++p) {
1221:       PetscInt dof, offset;

1223:       if (debug) {PetscSynchronizedPrintf(PetscObjectComm((PetscObject) dm), "[%d]  Checking root for claims %D\n", rank, p);}
1224:       PetscSectionGetDof(candidateSection, p, &dof);
1225:       PetscSectionGetOffset(candidateSection, p, &offset);
1226:       for (d = 0; d < dof;) {
1227:         if (claims[offset+d].rank >= 0) {
1228:           const PetscInt  faceInd   = offset+d;
1229:           const PetscInt  numPoints = candidates[faceInd].index;
1230:           const PetscInt *join      = NULL;
1231:           PetscInt        joinSize, points[1024], c;

1233:           if (debug) {PetscSynchronizedPrintf(PetscObjectComm((PetscObject) dm), "[%d]    Found claim for remote point (%D, %D)\n", rank, claims[faceInd].rank, claims[faceInd].index);}
1234:           points[0] = p;
1235:           if (debug) {PetscSynchronizedPrintf(PetscObjectComm((PetscObject) dm), "[%d]      point %D\n", rank, points[0]);}
1236:           for (c = 0, ++d; c < numPoints; ++c, ++d) {
1237:             key.i = candidates[offset+d].index;
1238:             key.j = candidates[offset+d].rank;
1239:             PetscHMapIJGet(roothash, key, &root);
1240:             points[c+1] = localPoints[root];
1241:             if (debug) {PetscSynchronizedPrintf(PetscObjectComm((PetscObject) dm), "[%d]      point %D\n", rank, points[c+1]);}
1242:           }
1243:           DMPlexGetJoin(dm, numPoints+1, points, &joinSize, &join);
1244:           if (joinSize == 1) {
1245:             if (debug) {PetscSynchronizedPrintf(PetscObjectComm((PetscObject) dm), "[%d]    Found local face %D\n", rank, join[0]);}
1246:             PetscHMapISet(claimshash, join[0], faceInd);
1247:           }
1248:           DMPlexRestoreJoin(dm, numPoints+1, points, &joinSize, &join);
1249:         } else d += claims[offset+d].index+1;
1250:       }
1251:     }
1252:     if (debug) {PetscSynchronizedFlush(PetscObjectComm((PetscObject) dm), NULL);}
1253:     /* Create new pointSF from hashed claims */
1254:     PetscHMapIGetSize(claimshash, &numLocalNew);
1255:     DMPlexGetChart(dm, &pStart, &pEnd);
1256:     PetscMalloc1(numLeaves + numLocalNew, &localPointsNew);
1257:     PetscMalloc1(numLeaves + numLocalNew, &remotePointsNew);
1258:     for (p = 0; p < numLeaves; ++p) {
1259:       localPointsNew[p] = localPoints[p];
1260:       remotePointsNew[p].index = remotePoints[p].index;
1261:       remotePointsNew[p].rank  = remotePoints[p].rank;
1262:     }
1263:     p = numLeaves;
1264:     PetscHMapIGetKeys(claimshash, &p, localPointsNew);
1265:     PetscSortInt(numLocalNew, &localPointsNew[numLeaves]);
1266:     for (p = numLeaves; p < numLeaves + numLocalNew; ++p) {
1267:       PetscInt offset;
1268:       PetscHMapIGet(claimshash, localPointsNew[p], &offset);
1269:       remotePointsNew[p] = claims[offset];
1270:     }
1271:     PetscSFCreate(PetscObjectComm((PetscObject) dm), &sfPointNew);
1272:     PetscSFSetGraph(sfPointNew, pEnd-pStart, numLeaves+numLocalNew, localPointsNew, PETSC_OWN_POINTER, remotePointsNew, PETSC_OWN_POINTER);
1273:     DMSetPointSF(dm, sfPointNew);
1274:     PetscSFDestroy(&sfPointNew);
1275:     PetscHMapIDestroy(&claimshash);
1276:   }
1277:   PetscHMapIDestroy(&leafhash);
1278:   PetscHMapIJDestroy(&roothash);
1279:   PetscSectionDestroy(&candidateSection);
1280:   PetscSectionDestroy(&candidateSectionRemote);
1281:   PetscSectionDestroy(&claimSection);
1282:   PetscFree(candidates);
1283:   PetscFree(candidatesRemote);
1284:   PetscFree(claims);
1285:   PetscLogEventEnd(DMPLEX_InterpolateSF,dm,0,0,0);
1286:   return(0);
1287: }

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

1292:   Collective on DM

1294:   Input Parameters:
1295: + dm - The DMPlex object with only cells and vertices
1296: - dmInt - The interpolated DM

1298:   Output Parameter:
1299: . dmInt - The complete DMPlex object

1301:   Level: intermediate

1303:   Notes:
1304:     It does not copy over the coordinates.

1306: .keywords: mesh
1307: .seealso: DMPlexUninterpolate(), DMPlexCreateFromCellList(), DMPlexCopyCoordinates()
1308: @*/
1309: PetscErrorCode DMPlexInterpolate(DM dm, DM *dmInt)
1310: {
1311:   DM             idm, odm = dm;
1312:   PetscSF        sfPoint;
1313:   PetscInt       depth, dim, d;
1314:   const char    *name;
1315:   PetscBool      flg=PETSC_TRUE;

1321:   PetscLogEventBegin(DMPLEX_Interpolate,dm,0,0,0);
1322:   DMPlexGetDepth(dm, &depth);
1323:   DMGetDimension(dm, &dim);
1324:   if ((depth == dim) || (dim <= 1)) {
1325:     PetscObjectReference((PetscObject) dm);
1326:     idm  = dm;
1327:   } else {
1328:     for (d = 1; d < dim; ++d) {
1329:       /* Create interpolated mesh */
1330:       DMCreate(PetscObjectComm((PetscObject)dm), &idm);
1331:       DMSetType(idm, DMPLEX);
1332:       DMSetDimension(idm, dim);
1333:       if (depth > 0) {
1334:         DMPlexInterpolateFaces_Internal(odm, 1, idm);
1335:         DMGetPointSF(odm, &sfPoint);
1336:         DMPlexInterpolatePointSF(idm, sfPoint);
1337:       }
1338:       if (odm != dm) {DMDestroy(&odm);}
1339:       odm = idm;
1340:     }
1341:     PetscObjectGetName((PetscObject) dm,  &name);
1342:     PetscObjectSetName((PetscObject) idm,  name);
1343:     DMPlexCopyCoordinates(dm, idm);
1344:     DMCopyLabels(dm, idm);
1345:     PetscOptionsGetBool(((PetscObject)dm)->options, ((PetscObject)dm)->prefix, "-dm_plex_interpolate_orient_interfaces", &flg, NULL);
1346:     if (flg) {DMPlexOrientInterface(idm);}
1347:   }
1348:   {
1349:     PetscBool            isper;
1350:     const PetscReal      *maxCell, *L;
1351:     const DMBoundaryType *bd;

1353:     DMGetPeriodicity(dm,&isper,&maxCell,&L,&bd);
1354:     DMSetPeriodicity(idm,isper,maxCell,L,bd);
1355:   }
1356:   *dmInt = idm;
1357:   PetscLogEventEnd(DMPLEX_Interpolate,dm,0,0,0);
1358:   return(0);
1359: }

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

1364:   Collective on DM

1366:   Input Parameter:
1367: . dmA - The DMPlex object with initial coordinates

1369:   Output Parameter:
1370: . dmB - The DMPlex object with copied coordinates

1372:   Level: intermediate

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

1376: .keywords: mesh
1377: .seealso: DMCopyLabels(), DMGetCoordinates(), DMGetCoordinatesLocal(), DMGetCoordinateDM(), DMGetCoordinateSection()
1378: @*/
1379: PetscErrorCode DMPlexCopyCoordinates(DM dmA, DM dmB)
1380: {
1381:   Vec            coordinatesA, coordinatesB;
1382:   VecType        vtype;
1383:   PetscSection   coordSectionA, coordSectionB;
1384:   PetscScalar   *coordsA, *coordsB;
1385:   PetscInt       spaceDim, Nf, vStartA, vStartB, vEndA, vEndB, coordSizeB, v, d;
1386:   PetscInt       cStartA, cEndA, cStartB, cEndB, cS, cE;
1387:   PetscBool      lc = PETSC_FALSE;

1393:   if (dmA == dmB) return(0);
1394:   DMPlexGetDepthStratum(dmA, 0, &vStartA, &vEndA);
1395:   DMPlexGetDepthStratum(dmB, 0, &vStartB, &vEndB);
1396:   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);
1397:   DMPlexGetHeightStratum(dmA, 0, &cStartA, &cEndA);
1398:   DMPlexGetHeightStratum(dmB, 0, &cStartB, &cEndB);
1399:   DMGetCoordinateSection(dmA, &coordSectionA);
1400:   DMGetCoordinateSection(dmB, &coordSectionB);
1401:   if (coordSectionA == coordSectionB) return(0);
1402:   PetscSectionGetNumFields(coordSectionA, &Nf);
1403:   if (!Nf) return(0);
1404:   if (Nf > 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "The number of coordinate fields must be 1, not %D", Nf);
1405:   if (!coordSectionB) {
1406:     PetscInt dim;

1408:     PetscSectionCreate(PetscObjectComm((PetscObject) coordSectionA), &coordSectionB);
1409:     DMGetCoordinateDim(dmA, &dim);
1410:     DMSetCoordinateSection(dmB, dim, coordSectionB);
1411:     PetscObjectDereference((PetscObject) coordSectionB);
1412:   }
1413:   PetscSectionSetNumFields(coordSectionB, 1);
1414:   PetscSectionGetFieldComponents(coordSectionA, 0, &spaceDim);
1415:   PetscSectionSetFieldComponents(coordSectionB, 0, spaceDim);
1416:   PetscSectionGetChart(coordSectionA, &cS, &cE);
1417:   if (cStartA <= cS && cS < cEndA) { /* localized coordinates */
1418:     if ((cEndA-cStartA) != (cEndB-cStartB)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "The number of cells in first DM %D != %D in the second DM", cEndA-cStartA, cEndB-cStartB);
1419:     cS = cS - cStartA + cStartB;
1420:     cE = vEndB;
1421:     lc = PETSC_TRUE;
1422:   } else {
1423:     cS = vStartB;
1424:     cE = vEndB;
1425:   }
1426:   PetscSectionSetChart(coordSectionB, cS, cE);
1427:   for (v = vStartB; v < vEndB; ++v) {
1428:     PetscSectionSetDof(coordSectionB, v, spaceDim);
1429:     PetscSectionSetFieldDof(coordSectionB, v, 0, spaceDim);
1430:   }
1431:   if (lc) { /* localized coordinates */
1432:     PetscInt c;

1434:     for (c = cS-cStartB; c < cEndB-cStartB; c++) {
1435:       PetscInt dof;

1437:       PetscSectionGetDof(coordSectionA, c + cStartA, &dof);
1438:       PetscSectionSetDof(coordSectionB, c + cStartB, dof);
1439:       PetscSectionSetFieldDof(coordSectionB, c + cStartB, 0, dof);
1440:     }
1441:   }
1442:   PetscSectionSetUp(coordSectionB);
1443:   PetscSectionGetStorageSize(coordSectionB, &coordSizeB);
1444:   DMGetCoordinatesLocal(dmA, &coordinatesA);
1445:   VecCreate(PETSC_COMM_SELF, &coordinatesB);
1446:   PetscObjectSetName((PetscObject) coordinatesB, "coordinates");
1447:   VecSetSizes(coordinatesB, coordSizeB, PETSC_DETERMINE);
1448:   VecGetBlockSize(coordinatesA, &d);
1449:   VecSetBlockSize(coordinatesB, d);
1450:   VecGetType(coordinatesA, &vtype);
1451:   VecSetType(coordinatesB, vtype);
1452:   VecGetArray(coordinatesA, &coordsA);
1453:   VecGetArray(coordinatesB, &coordsB);
1454:   for (v = 0; v < vEndB-vStartB; ++v) {
1455:     PetscInt offA, offB;

1457:     PetscSectionGetOffset(coordSectionA, v + vStartA, &offA);
1458:     PetscSectionGetOffset(coordSectionB, v + vStartB, &offB);
1459:     for (d = 0; d < spaceDim; ++d) {
1460:       coordsB[offB+d] = coordsA[offA+d];
1461:     }
1462:   }
1463:   if (lc) { /* localized coordinates */
1464:     PetscInt c;

1466:     for (c = cS-cStartB; c < cEndB-cStartB; c++) {
1467:       PetscInt dof, offA, offB;

1469:       PetscSectionGetOffset(coordSectionA, c + cStartA, &offA);
1470:       PetscSectionGetOffset(coordSectionB, c + cStartB, &offB);
1471:       PetscSectionGetDof(coordSectionA, c + cStartA, &dof);
1472:       PetscMemcpy(coordsB + offB,coordsA + offA,dof*sizeof(*coordsB));
1473:     }
1474:   }
1475:   VecRestoreArray(coordinatesA, &coordsA);
1476:   VecRestoreArray(coordinatesB, &coordsB);
1477:   DMSetCoordinatesLocal(dmB, coordinatesB);
1478:   VecDestroy(&coordinatesB);
1479:   return(0);
1480: }

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

1485:   Collective on DM

1487:   Input Parameter:
1488: . dm - The complete DMPlex object

1490:   Output Parameter:
1491: . dmUnint - The DMPlex object with only cells and vertices

1493:   Level: intermediate

1495:   Notes:
1496:     It does not copy over the coordinates.

1498: .keywords: mesh
1499: .seealso: DMPlexInterpolate(), DMPlexCreateFromCellList(), DMPlexCopyCoordinates()
1500: @*/
1501: PetscErrorCode DMPlexUninterpolate(DM dm, DM *dmUnint)
1502: {
1503:   DM             udm;
1504:   PetscInt       dim, vStart, vEnd, cStart, cEnd, cMax, c, maxConeSize = 0, *cone;

1510:   DMGetDimension(dm, &dim);
1511:   if (dim <= 1) {
1512:     PetscObjectReference((PetscObject) dm);
1513:     *dmUnint = dm;
1514:     return(0);
1515:   }
1516:   DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
1517:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
1518:   DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);
1519:   DMCreate(PetscObjectComm((PetscObject) dm), &udm);
1520:   DMSetType(udm, DMPLEX);
1521:   DMSetDimension(udm, dim);
1522:   DMPlexSetChart(udm, cStart, vEnd);
1523:   for (c = cStart; c < cEnd; ++c) {
1524:     PetscInt *closure = NULL, closureSize, cl, coneSize = 0;

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

1530:       if ((p >= vStart) && (p < vEnd)) ++coneSize;
1531:     }
1532:     DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
1533:     DMPlexSetConeSize(udm, c, coneSize);
1534:     maxConeSize = PetscMax(maxConeSize, coneSize);
1535:   }
1536:   DMSetUp(udm);
1537:   PetscMalloc1(maxConeSize, &cone);
1538:   for (c = cStart; c < cEnd; ++c) {
1539:     PetscInt *closure = NULL, closureSize, cl, coneSize = 0;

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

1545:       if ((p >= vStart) && (p < vEnd)) cone[coneSize++] = p;
1546:     }
1547:     DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
1548:     DMPlexSetCone(udm, c, cone);
1549:   }
1550:   PetscFree(cone);
1551:   DMPlexSetHybridBounds(udm, cMax, PETSC_DETERMINE, PETSC_DETERMINE, PETSC_DETERMINE);
1552:   DMPlexSymmetrize(udm);
1553:   DMPlexStratify(udm);
1554:   /* Reduce SF */
1555:   {
1556:     PetscSF            sfPoint, sfPointUn;
1557:     const PetscSFNode *remotePoints;
1558:     const PetscInt    *localPoints;
1559:     PetscSFNode       *remotePointsUn;
1560:     PetscInt          *localPointsUn;
1561:     PetscInt           vEnd, numRoots, numLeaves, l;
1562:     PetscInt           numLeavesUn = 0, n = 0;
1563:     PetscErrorCode     ierr;

1565:     /* Get original SF information */
1566:     DMGetPointSF(dm, &sfPoint);
1567:     DMGetPointSF(udm, &sfPointUn);
1568:     DMPlexGetDepthStratum(dm, 0, NULL, &vEnd);
1569:     PetscSFGetGraph(sfPoint, &numRoots, &numLeaves, &localPoints, &remotePoints);
1570:     /* Allocate space for cells and vertices */
1571:     for (l = 0; l < numLeaves; ++l) if (localPoints[l] < vEnd) numLeavesUn++;
1572:     /* Fill in leaves */
1573:     if (vEnd >= 0) {
1574:       PetscMalloc1(numLeavesUn, &remotePointsUn);
1575:       PetscMalloc1(numLeavesUn, &localPointsUn);
1576:       for (l = 0; l < numLeaves; l++) {
1577:         if (localPoints[l] < vEnd) {
1578:           localPointsUn[n]        = localPoints[l];
1579:           remotePointsUn[n].rank  = remotePoints[l].rank;
1580:           remotePointsUn[n].index = remotePoints[l].index;
1581:           ++n;
1582:         }
1583:       }
1584:       if (n != numLeavesUn) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Inconsistent number of leaves %d != %d", n, numLeavesUn);
1585:       PetscSFSetGraph(sfPointUn, vEnd, numLeavesUn, localPointsUn, PETSC_OWN_POINTER, remotePointsUn, PETSC_OWN_POINTER);
1586:     }
1587:   }
1588:   {
1589:     PetscBool            isper;
1590:     const PetscReal      *maxCell, *L;
1591:     const DMBoundaryType *bd;

1593:     DMGetPeriodicity(dm,&isper,&maxCell,&L,&bd);
1594:     DMSetPeriodicity(udm,isper,maxCell,L,bd);
1595:   }

1597:   *dmUnint = udm;
1598:   return(0);
1599: }