Actual source code: plexinterpolate.c

  1: #include <petsc/private/dmpleximpl.h>
  2: #include <petsc/private/hashmapi.h>
  3: #include <petsc/private/hashmapij.h>

  5: const char *const DMPlexInterpolatedFlags[] = {"none", "partial", "mixed", "full", "DMPlexInterpolatedFlag", "DMPLEX_INTERPOLATED_", NULL};

  7: /* HashIJKL */

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

 11: typedef struct _PetscHashIJKLKey {
 12:   PetscInt i, j, k, l;
 13: } PetscHashIJKLKey;

 15: #define PetscHashIJKLKeyHash(key) PetscHashCombine(PetscHashCombine(PetscHashInt((key).i), PetscHashInt((key).j)), PetscHashCombine(PetscHashInt((key).k), PetscHashInt((key).l)))

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

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

 21:   static PetscSFNode _PetscInvalidSFNode = {-1, -1};

 23: typedef struct _PetscHashIJKLRemoteKey {
 24:   PetscSFNode i, j, k, l;
 25: } PetscHashIJKLRemoteKey;

 27: #define PetscHashIJKLRemoteKeyHash(key) \
 28:   PetscHashCombine(PetscHashCombine(PetscHashInt((key).i.rank + (key).i.index), PetscHashInt((key).j.rank + (key).j.index)), PetscHashCombine(PetscHashInt((key).k.rank + (key).k.index), PetscHashInt((key).l.rank + (key).l.index)))

 30: #define PetscHashIJKLRemoteKeyEqual(k1, k2) \
 31:   (((k1).i.rank == (k2).i.rank) ? ((k1).i.index == (k2).i.index) ? ((k1).j.rank == (k2).j.rank) ? ((k1).j.index == (k2).j.index) ? ((k1).k.rank == (k2).k.rank) ? ((k1).k.index == (k2).k.index) ? ((k1).l.rank == (k2).l.rank) ? ((k1).l.index == (k2).l.index) : 0 : 0 : 0 : 0 : 0 : 0 : 0)

 33: PetscDisableStaticAnalyzerForExpressionUnderstandingThatThisIsDangerousAndBugprone(PETSC_HASH_MAP(HashIJKLRemote, PetscHashIJKLRemoteKey, PetscSFNode, PetscHashIJKLRemoteKeyHash, PetscHashIJKLRemoteKeyEqual, _PetscInvalidSFNode))

 35:   static PetscErrorCode PetscSortSFNode(PetscInt n, PetscSFNode A[])
 36: {
 37:   PetscInt i;

 39:   for (i = 1; i < n; ++i) {
 40:     PetscSFNode x = A[i];
 41:     PetscInt    j;

 43:     for (j = i - 1; j >= 0; --j) {
 44:       if ((A[j].rank > x.rank) || (A[j].rank == x.rank && A[j].index > x.index)) break;
 45:       A[j + 1] = A[j];
 46:     }
 47:     A[j + 1] = x;
 48:   }
 49:   return 0;
 50: }

 52: /*
 53:   DMPlexGetRawFaces_Internal - Gets groups of vertices that correspond to faces for the given cone
 54: */
 55: PetscErrorCode DMPlexGetRawFaces_Internal(DM dm, DMPolytopeType ct, const PetscInt cone[], PetscInt *numFaces, const DMPolytopeType *faceTypes[], const PetscInt *faceSizes[], const PetscInt *faces[])
 56: {
 57:   DMPolytopeType *typesTmp;
 58:   PetscInt       *sizesTmp, *facesTmp;
 59:   PetscInt        maxConeSize, maxSupportSize;

 63:   DMPlexGetMaxSizes(dm, &maxConeSize, &maxSupportSize);
 64:   if (faceTypes) DMGetWorkArray(dm, PetscMax(maxConeSize, maxSupportSize), MPIU_INT, &typesTmp);
 65:   if (faceSizes) DMGetWorkArray(dm, PetscMax(maxConeSize, maxSupportSize), MPIU_INT, &sizesTmp);
 66:   if (faces) DMGetWorkArray(dm, PetscSqr(PetscMax(maxConeSize, maxSupportSize)), MPIU_INT, &facesTmp);
 67:   switch (ct) {
 68:   case DM_POLYTOPE_POINT:
 69:     if (numFaces) *numFaces = 0;
 70:     break;
 71:   case DM_POLYTOPE_SEGMENT:
 72:     if (numFaces) *numFaces = 2;
 73:     if (faceTypes) {
 74:       typesTmp[0] = DM_POLYTOPE_POINT;
 75:       typesTmp[1] = DM_POLYTOPE_POINT;
 76:       *faceTypes  = typesTmp;
 77:     }
 78:     if (faceSizes) {
 79:       sizesTmp[0] = 1;
 80:       sizesTmp[1] = 1;
 81:       *faceSizes  = sizesTmp;
 82:     }
 83:     if (faces) {
 84:       facesTmp[0] = cone[0];
 85:       facesTmp[1] = cone[1];
 86:       *faces      = facesTmp;
 87:     }
 88:     break;
 89:   case DM_POLYTOPE_POINT_PRISM_TENSOR:
 90:     if (numFaces) *numFaces = 2;
 91:     if (faceTypes) {
 92:       typesTmp[0] = DM_POLYTOPE_POINT;
 93:       typesTmp[1] = DM_POLYTOPE_POINT;
 94:       *faceTypes  = typesTmp;
 95:     }
 96:     if (faceSizes) {
 97:       sizesTmp[0] = 1;
 98:       sizesTmp[1] = 1;
 99:       *faceSizes  = sizesTmp;
100:     }
101:     if (faces) {
102:       facesTmp[0] = cone[0];
103:       facesTmp[1] = cone[1];
104:       *faces      = facesTmp;
105:     }
106:     break;
107:   case DM_POLYTOPE_TRIANGLE:
108:     if (numFaces) *numFaces = 3;
109:     if (faceTypes) {
110:       typesTmp[0] = DM_POLYTOPE_SEGMENT;
111:       typesTmp[1] = DM_POLYTOPE_SEGMENT;
112:       typesTmp[2] = DM_POLYTOPE_SEGMENT;
113:       *faceTypes  = typesTmp;
114:     }
115:     if (faceSizes) {
116:       sizesTmp[0] = 2;
117:       sizesTmp[1] = 2;
118:       sizesTmp[2] = 2;
119:       *faceSizes  = sizesTmp;
120:     }
121:     if (faces) {
122:       facesTmp[0] = cone[0];
123:       facesTmp[1] = cone[1];
124:       facesTmp[2] = cone[1];
125:       facesTmp[3] = cone[2];
126:       facesTmp[4] = cone[2];
127:       facesTmp[5] = cone[0];
128:       *faces      = facesTmp;
129:     }
130:     break;
131:   case DM_POLYTOPE_QUADRILATERAL:
132:     /* Vertices follow right hand rule */
133:     if (numFaces) *numFaces = 4;
134:     if (faceTypes) {
135:       typesTmp[0] = DM_POLYTOPE_SEGMENT;
136:       typesTmp[1] = DM_POLYTOPE_SEGMENT;
137:       typesTmp[2] = DM_POLYTOPE_SEGMENT;
138:       typesTmp[3] = DM_POLYTOPE_SEGMENT;
139:       *faceTypes  = typesTmp;
140:     }
141:     if (faceSizes) {
142:       sizesTmp[0] = 2;
143:       sizesTmp[1] = 2;
144:       sizesTmp[2] = 2;
145:       sizesTmp[3] = 2;
146:       *faceSizes  = sizesTmp;
147:     }
148:     if (faces) {
149:       facesTmp[0] = cone[0];
150:       facesTmp[1] = cone[1];
151:       facesTmp[2] = cone[1];
152:       facesTmp[3] = cone[2];
153:       facesTmp[4] = cone[2];
154:       facesTmp[5] = cone[3];
155:       facesTmp[6] = cone[3];
156:       facesTmp[7] = cone[0];
157:       *faces      = facesTmp;
158:     }
159:     break;
160:   case DM_POLYTOPE_SEG_PRISM_TENSOR:
161:     if (numFaces) *numFaces = 4;
162:     if (faceTypes) {
163:       typesTmp[0] = DM_POLYTOPE_SEGMENT;
164:       typesTmp[1] = DM_POLYTOPE_SEGMENT;
165:       typesTmp[2] = DM_POLYTOPE_POINT_PRISM_TENSOR;
166:       typesTmp[3] = DM_POLYTOPE_POINT_PRISM_TENSOR;
167:       *faceTypes  = typesTmp;
168:     }
169:     if (faceSizes) {
170:       sizesTmp[0] = 2;
171:       sizesTmp[1] = 2;
172:       sizesTmp[2] = 2;
173:       sizesTmp[3] = 2;
174:       *faceSizes  = sizesTmp;
175:     }
176:     if (faces) {
177:       facesTmp[0] = cone[0];
178:       facesTmp[1] = cone[1];
179:       facesTmp[2] = cone[2];
180:       facesTmp[3] = cone[3];
181:       facesTmp[4] = cone[0];
182:       facesTmp[5] = cone[2];
183:       facesTmp[6] = cone[1];
184:       facesTmp[7] = cone[3];
185:       *faces      = facesTmp;
186:     }
187:     break;
188:   case DM_POLYTOPE_TETRAHEDRON:
189:     /* Vertices of first face follow right hand rule and normal points away from last vertex */
190:     if (numFaces) *numFaces = 4;
191:     if (faceTypes) {
192:       typesTmp[0] = DM_POLYTOPE_TRIANGLE;
193:       typesTmp[1] = DM_POLYTOPE_TRIANGLE;
194:       typesTmp[2] = DM_POLYTOPE_TRIANGLE;
195:       typesTmp[3] = DM_POLYTOPE_TRIANGLE;
196:       *faceTypes  = typesTmp;
197:     }
198:     if (faceSizes) {
199:       sizesTmp[0] = 3;
200:       sizesTmp[1] = 3;
201:       sizesTmp[2] = 3;
202:       sizesTmp[3] = 3;
203:       *faceSizes  = sizesTmp;
204:     }
205:     if (faces) {
206:       facesTmp[0]  = cone[0];
207:       facesTmp[1]  = cone[1];
208:       facesTmp[2]  = cone[2];
209:       facesTmp[3]  = cone[0];
210:       facesTmp[4]  = cone[3];
211:       facesTmp[5]  = cone[1];
212:       facesTmp[6]  = cone[0];
213:       facesTmp[7]  = cone[2];
214:       facesTmp[8]  = cone[3];
215:       facesTmp[9]  = cone[2];
216:       facesTmp[10] = cone[1];
217:       facesTmp[11] = cone[3];
218:       *faces       = facesTmp;
219:     }
220:     break;
221:   case DM_POLYTOPE_HEXAHEDRON:
222:     /*  7--------6
223:          /|       /|
224:         / |      / |
225:        4--------5  |
226:        |  |     |  |
227:        |  |     |  |
228:        |  1--------2
229:        | /      | /
230:        |/       |/
231:        0--------3
232:        */
233:     if (numFaces) *numFaces = 6;
234:     if (faceTypes) {
235:       typesTmp[0] = DM_POLYTOPE_QUADRILATERAL;
236:       typesTmp[1] = DM_POLYTOPE_QUADRILATERAL;
237:       typesTmp[2] = DM_POLYTOPE_QUADRILATERAL;
238:       typesTmp[3] = DM_POLYTOPE_QUADRILATERAL;
239:       typesTmp[4] = DM_POLYTOPE_QUADRILATERAL;
240:       typesTmp[5] = DM_POLYTOPE_QUADRILATERAL;
241:       *faceTypes  = typesTmp;
242:     }
243:     if (faceSizes) {
244:       sizesTmp[0] = 4;
245:       sizesTmp[1] = 4;
246:       sizesTmp[2] = 4;
247:       sizesTmp[3] = 4;
248:       sizesTmp[4] = 4;
249:       sizesTmp[5] = 4;
250:       *faceSizes  = sizesTmp;
251:     }
252:     if (faces) {
253:       facesTmp[0]  = cone[0];
254:       facesTmp[1]  = cone[1];
255:       facesTmp[2]  = cone[2];
256:       facesTmp[3]  = cone[3]; /* Bottom */
257:       facesTmp[4]  = cone[4];
258:       facesTmp[5]  = cone[5];
259:       facesTmp[6]  = cone[6];
260:       facesTmp[7]  = cone[7]; /* Top */
261:       facesTmp[8]  = cone[0];
262:       facesTmp[9]  = cone[3];
263:       facesTmp[10] = cone[5];
264:       facesTmp[11] = cone[4]; /* Front */
265:       facesTmp[12] = cone[2];
266:       facesTmp[13] = cone[1];
267:       facesTmp[14] = cone[7];
268:       facesTmp[15] = cone[6]; /* Back */
269:       facesTmp[16] = cone[3];
270:       facesTmp[17] = cone[2];
271:       facesTmp[18] = cone[6];
272:       facesTmp[19] = cone[5]; /* Right */
273:       facesTmp[20] = cone[0];
274:       facesTmp[21] = cone[4];
275:       facesTmp[22] = cone[7];
276:       facesTmp[23] = cone[1]; /* Left */
277:       *faces       = facesTmp;
278:     }
279:     break;
280:   case DM_POLYTOPE_TRI_PRISM:
281:     if (numFaces) *numFaces = 5;
282:     if (faceTypes) {
283:       typesTmp[0] = DM_POLYTOPE_TRIANGLE;
284:       typesTmp[1] = DM_POLYTOPE_TRIANGLE;
285:       typesTmp[2] = DM_POLYTOPE_QUADRILATERAL;
286:       typesTmp[3] = DM_POLYTOPE_QUADRILATERAL;
287:       typesTmp[4] = DM_POLYTOPE_QUADRILATERAL;
288:       *faceTypes  = typesTmp;
289:     }
290:     if (faceSizes) {
291:       sizesTmp[0] = 3;
292:       sizesTmp[1] = 3;
293:       sizesTmp[2] = 4;
294:       sizesTmp[3] = 4;
295:       sizesTmp[4] = 4;
296:       *faceSizes  = sizesTmp;
297:     }
298:     if (faces) {
299:       facesTmp[0]  = cone[0];
300:       facesTmp[1]  = cone[1];
301:       facesTmp[2]  = cone[2]; /* Bottom */
302:       facesTmp[3]  = cone[3];
303:       facesTmp[4]  = cone[4];
304:       facesTmp[5]  = cone[5]; /* Top */
305:       facesTmp[6]  = cone[0];
306:       facesTmp[7]  = cone[2];
307:       facesTmp[8]  = cone[4];
308:       facesTmp[9]  = cone[3]; /* Back left */
309:       facesTmp[10] = cone[2];
310:       facesTmp[11] = cone[1];
311:       facesTmp[12] = cone[5];
312:       facesTmp[13] = cone[4]; /* Front */
313:       facesTmp[14] = cone[1];
314:       facesTmp[15] = cone[0];
315:       facesTmp[16] = cone[3];
316:       facesTmp[17] = cone[5]; /* Back right */
317:       *faces       = facesTmp;
318:     }
319:     break;
320:   case DM_POLYTOPE_TRI_PRISM_TENSOR:
321:     if (numFaces) *numFaces = 5;
322:     if (faceTypes) {
323:       typesTmp[0] = DM_POLYTOPE_TRIANGLE;
324:       typesTmp[1] = DM_POLYTOPE_TRIANGLE;
325:       typesTmp[2] = DM_POLYTOPE_SEG_PRISM_TENSOR;
326:       typesTmp[3] = DM_POLYTOPE_SEG_PRISM_TENSOR;
327:       typesTmp[4] = DM_POLYTOPE_SEG_PRISM_TENSOR;
328:       *faceTypes  = typesTmp;
329:     }
330:     if (faceSizes) {
331:       sizesTmp[0] = 3;
332:       sizesTmp[1] = 3;
333:       sizesTmp[2] = 4;
334:       sizesTmp[3] = 4;
335:       sizesTmp[4] = 4;
336:       *faceSizes  = sizesTmp;
337:     }
338:     if (faces) {
339:       facesTmp[0]  = cone[0];
340:       facesTmp[1]  = cone[1];
341:       facesTmp[2]  = cone[2]; /* Bottom */
342:       facesTmp[3]  = cone[3];
343:       facesTmp[4]  = cone[4];
344:       facesTmp[5]  = cone[5]; /* Top */
345:       facesTmp[6]  = cone[0];
346:       facesTmp[7]  = cone[1];
347:       facesTmp[8]  = cone[3];
348:       facesTmp[9]  = cone[4]; /* Back left */
349:       facesTmp[10] = cone[1];
350:       facesTmp[11] = cone[2];
351:       facesTmp[12] = cone[4];
352:       facesTmp[13] = cone[5]; /* Back right */
353:       facesTmp[14] = cone[2];
354:       facesTmp[15] = cone[0];
355:       facesTmp[16] = cone[5];
356:       facesTmp[17] = cone[3]; /* Front */
357:       *faces       = facesTmp;
358:     }
359:     break;
360:   case DM_POLYTOPE_QUAD_PRISM_TENSOR:
361:     /*  7--------6
362:          /|       /|
363:         / |      / |
364:        4--------5  |
365:        |  |     |  |
366:        |  |     |  |
367:        |  3--------2
368:        | /      | /
369:        |/       |/
370:        0--------1
371:        */
372:     if (numFaces) *numFaces = 6;
373:     if (faceTypes) {
374:       typesTmp[0] = DM_POLYTOPE_QUADRILATERAL;
375:       typesTmp[1] = DM_POLYTOPE_QUADRILATERAL;
376:       typesTmp[2] = DM_POLYTOPE_SEG_PRISM_TENSOR;
377:       typesTmp[3] = DM_POLYTOPE_SEG_PRISM_TENSOR;
378:       typesTmp[4] = DM_POLYTOPE_SEG_PRISM_TENSOR;
379:       typesTmp[5] = DM_POLYTOPE_SEG_PRISM_TENSOR;
380:       *faceTypes  = typesTmp;
381:     }
382:     if (faceSizes) {
383:       sizesTmp[0] = 4;
384:       sizesTmp[1] = 4;
385:       sizesTmp[2] = 4;
386:       sizesTmp[3] = 4;
387:       sizesTmp[4] = 4;
388:       sizesTmp[5] = 4;
389:       *faceSizes  = sizesTmp;
390:     }
391:     if (faces) {
392:       facesTmp[0]  = cone[0];
393:       facesTmp[1]  = cone[1];
394:       facesTmp[2]  = cone[2];
395:       facesTmp[3]  = cone[3]; /* Bottom */
396:       facesTmp[4]  = cone[4];
397:       facesTmp[5]  = cone[5];
398:       facesTmp[6]  = cone[6];
399:       facesTmp[7]  = cone[7]; /* Top */
400:       facesTmp[8]  = cone[0];
401:       facesTmp[9]  = cone[1];
402:       facesTmp[10] = cone[4];
403:       facesTmp[11] = cone[5]; /* Front */
404:       facesTmp[12] = cone[1];
405:       facesTmp[13] = cone[2];
406:       facesTmp[14] = cone[5];
407:       facesTmp[15] = cone[6]; /* Right */
408:       facesTmp[16] = cone[2];
409:       facesTmp[17] = cone[3];
410:       facesTmp[18] = cone[6];
411:       facesTmp[19] = cone[7]; /* Back */
412:       facesTmp[20] = cone[3];
413:       facesTmp[21] = cone[0];
414:       facesTmp[22] = cone[7];
415:       facesTmp[23] = cone[4]; /* Left */
416:       *faces       = facesTmp;
417:     }
418:     break;
419:   case DM_POLYTOPE_PYRAMID:
420:     /*
421:        4----
422:        |\-\ \-----
423:        | \ -\     \
424:        |  1--\-----2
425:        | /    \   /
426:        |/      \ /
427:        0--------3
428:        */
429:     if (numFaces) *numFaces = 5;
430:     if (faceTypes) {
431:       typesTmp[0] = DM_POLYTOPE_QUADRILATERAL;
432:       typesTmp[1] = DM_POLYTOPE_TRIANGLE;
433:       typesTmp[2] = DM_POLYTOPE_TRIANGLE;
434:       typesTmp[3] = DM_POLYTOPE_TRIANGLE;
435:       typesTmp[4] = DM_POLYTOPE_TRIANGLE;
436:       *faceTypes  = typesTmp;
437:     }
438:     if (faceSizes) {
439:       sizesTmp[0] = 4;
440:       sizesTmp[1] = 3;
441:       sizesTmp[2] = 3;
442:       sizesTmp[3] = 3;
443:       sizesTmp[4] = 3;
444:       *faceSizes  = sizesTmp;
445:     }
446:     if (faces) {
447:       facesTmp[0]  = cone[0];
448:       facesTmp[1]  = cone[1];
449:       facesTmp[2]  = cone[2];
450:       facesTmp[3]  = cone[3]; /* Bottom */
451:       facesTmp[4]  = cone[0];
452:       facesTmp[5]  = cone[3];
453:       facesTmp[6]  = cone[4]; /* Front */
454:       facesTmp[7]  = cone[3];
455:       facesTmp[8]  = cone[2];
456:       facesTmp[9]  = cone[4]; /* Right */
457:       facesTmp[10] = cone[2];
458:       facesTmp[11] = cone[1];
459:       facesTmp[12] = cone[4]; /* Back */
460:       facesTmp[13] = cone[1];
461:       facesTmp[14] = cone[0];
462:       facesTmp[15] = cone[4]; /* Left */
463:       *faces       = facesTmp;
464:     }
465:     break;
466:   default:
467:     SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "No face description for cell type %s", DMPolytopeTypes[ct]);
468:   }
469:   return 0;
470: }

472: PetscErrorCode DMPlexRestoreRawFaces_Internal(DM dm, DMPolytopeType ct, const PetscInt cone[], PetscInt *numFaces, const DMPolytopeType *faceTypes[], const PetscInt *faceSizes[], const PetscInt *faces[])
473: {
474:   if (faceTypes) DMRestoreWorkArray(dm, 0, MPIU_INT, (void *)faceTypes);
475:   if (faceSizes) DMRestoreWorkArray(dm, 0, MPIU_INT, (void *)faceSizes);
476:   if (faces) DMRestoreWorkArray(dm, 0, MPIU_INT, (void *)faces);
477:   return 0;
478: }

480: /* This interpolates faces for cells at some stratum */
481: static PetscErrorCode DMPlexInterpolateFaces_Internal(DM dm, PetscInt cellDepth, DM idm)
482: {
483:   DMLabel       ctLabel;
484:   PetscHashIJKL faceTable;
485:   PetscInt      faceTypeNum[DM_NUM_POLYTOPES];
486:   PetscInt      depth, d, pStart, Np, cStart, cEnd, c, fStart, fEnd;

488:   DMPlexGetDepth(dm, &depth);
489:   PetscHashIJKLCreate(&faceTable);
490:   PetscArrayzero(faceTypeNum, DM_NUM_POLYTOPES);
491:   DMPlexGetDepthStratum(dm, cellDepth, &cStart, &cEnd);
492:   /* Number new faces and save face vertices in hash table */
493:   DMPlexGetDepthStratum(dm, depth > cellDepth ? cellDepth : 0, NULL, &fStart);
494:   fEnd = fStart;
495:   for (c = cStart; c < cEnd; ++c) {
496:     const PetscInt       *cone, *faceSizes, *faces;
497:     const DMPolytopeType *faceTypes;
498:     DMPolytopeType        ct;
499:     PetscInt              numFaces, cf, foff = 0;

501:     DMPlexGetCellType(dm, c, &ct);
502:     DMPlexGetCone(dm, c, &cone);
503:     DMPlexGetRawFaces_Internal(dm, ct, cone, &numFaces, &faceTypes, &faceSizes, &faces);
504:     for (cf = 0; cf < numFaces; foff += faceSizes[cf], ++cf) {
505:       const PetscInt       faceSize = faceSizes[cf];
506:       const DMPolytopeType faceType = faceTypes[cf];
507:       const PetscInt      *face     = &faces[foff];
508:       PetscHashIJKLKey     key;
509:       PetscHashIter        iter;
510:       PetscBool            missing;

513:       key.i = face[0];
514:       key.j = faceSize > 1 ? face[1] : PETSC_MAX_INT;
515:       key.k = faceSize > 2 ? face[2] : PETSC_MAX_INT;
516:       key.l = faceSize > 3 ? face[3] : PETSC_MAX_INT;
517:       PetscSortInt(faceSize, (PetscInt *)&key);
518:       PetscHashIJKLPut(faceTable, key, &iter, &missing);
519:       if (missing) {
520:         PetscHashIJKLIterSet(faceTable, iter, fEnd++);
521:         ++faceTypeNum[faceType];
522:       }
523:     }
524:     DMPlexRestoreRawFaces_Internal(dm, ct, cone, &numFaces, &faceTypes, &faceSizes, &faces);
525:   }
526:   /* We need to number faces contiguously among types */
527:   {
528:     PetscInt faceTypeStart[DM_NUM_POLYTOPES], ct, numFT = 0;

530:     for (ct = 0; ct < DM_NUM_POLYTOPES; ++ct) {
531:       if (faceTypeNum[ct]) ++numFT;
532:       faceTypeStart[ct] = 0;
533:     }
534:     if (numFT > 1) {
535:       PetscHashIJKLClear(faceTable);
536:       faceTypeStart[0] = fStart;
537:       for (ct = 1; ct < DM_NUM_POLYTOPES; ++ct) faceTypeStart[ct] = faceTypeStart[ct - 1] + faceTypeNum[ct - 1];
538:       for (c = cStart; c < cEnd; ++c) {
539:         const PetscInt       *cone, *faceSizes, *faces;
540:         const DMPolytopeType *faceTypes;
541:         DMPolytopeType        ct;
542:         PetscInt              numFaces, cf, foff = 0;

544:         DMPlexGetCellType(dm, c, &ct);
545:         DMPlexGetCone(dm, c, &cone);
546:         DMPlexGetRawFaces_Internal(dm, ct, cone, &numFaces, &faceTypes, &faceSizes, &faces);
547:         for (cf = 0; cf < numFaces; foff += faceSizes[cf], ++cf) {
548:           const PetscInt       faceSize = faceSizes[cf];
549:           const DMPolytopeType faceType = faceTypes[cf];
550:           const PetscInt      *face     = &faces[foff];
551:           PetscHashIJKLKey     key;
552:           PetscHashIter        iter;
553:           PetscBool            missing;

556:           key.i = face[0];
557:           key.j = faceSize > 1 ? face[1] : PETSC_MAX_INT;
558:           key.k = faceSize > 2 ? face[2] : PETSC_MAX_INT;
559:           key.l = faceSize > 3 ? face[3] : PETSC_MAX_INT;
560:           PetscSortInt(faceSize, (PetscInt *)&key);
561:           PetscHashIJKLPut(faceTable, key, &iter, &missing);
562:           if (missing) PetscHashIJKLIterSet(faceTable, iter, faceTypeStart[faceType]++);
563:         }
564:         DMPlexRestoreRawFaces_Internal(dm, ct, cone, &numFaces, &faceTypes, &faceSizes, &faces);
565:       }
566:       for (ct = 1; ct < DM_NUM_POLYTOPES; ++ct) {
568:       }
569:     }
570:   }
571:   /* Add new points, always at the end of the numbering */
572:   DMPlexGetChart(dm, &pStart, &Np);
573:   DMPlexSetChart(idm, pStart, Np + (fEnd - fStart));
574:   /* Set cone sizes */
575:   /*   Must create the celltype label here so that we do not automatically try to compute the types */
576:   DMCreateLabel(idm, "celltype");
577:   DMPlexGetCellTypeLabel(idm, &ctLabel);
578:   for (d = 0; d <= depth; ++d) {
579:     DMPolytopeType ct;
580:     PetscInt       coneSize, pStart, pEnd, p;

582:     if (d == cellDepth) continue;
583:     DMPlexGetDepthStratum(dm, d, &pStart, &pEnd);
584:     for (p = pStart; p < pEnd; ++p) {
585:       DMPlexGetConeSize(dm, p, &coneSize);
586:       DMPlexSetConeSize(idm, p, coneSize);
587:       DMPlexGetCellType(dm, p, &ct);
588:       DMPlexSetCellType(idm, p, ct);
589:     }
590:   }
591:   for (c = cStart; c < cEnd; ++c) {
592:     const PetscInt       *cone, *faceSizes, *faces;
593:     const DMPolytopeType *faceTypes;
594:     DMPolytopeType        ct;
595:     PetscInt              numFaces, cf, foff = 0;

597:     DMPlexGetCellType(dm, c, &ct);
598:     DMPlexGetCone(dm, c, &cone);
599:     DMPlexGetRawFaces_Internal(dm, ct, cone, &numFaces, &faceTypes, &faceSizes, &faces);
600:     DMPlexSetCellType(idm, c, ct);
601:     DMPlexSetConeSize(idm, c, numFaces);
602:     for (cf = 0; cf < numFaces; foff += faceSizes[cf], ++cf) {
603:       const PetscInt       faceSize = faceSizes[cf];
604:       const DMPolytopeType faceType = faceTypes[cf];
605:       const PetscInt      *face     = &faces[foff];
606:       PetscHashIJKLKey     key;
607:       PetscHashIter        iter;
608:       PetscBool            missing;
609:       PetscInt             f;

612:       key.i = face[0];
613:       key.j = faceSize > 1 ? face[1] : PETSC_MAX_INT;
614:       key.k = faceSize > 2 ? face[2] : PETSC_MAX_INT;
615:       key.l = faceSize > 3 ? face[3] : PETSC_MAX_INT;
616:       PetscSortInt(faceSize, (PetscInt *)&key);
617:       PetscHashIJKLPut(faceTable, key, &iter, &missing);
619:       PetscHashIJKLIterGet(faceTable, iter, &f);
620:       DMPlexSetConeSize(idm, f, faceSize);
621:       DMPlexSetCellType(idm, f, faceType);
622:     }
623:     DMPlexRestoreRawFaces_Internal(dm, ct, cone, &numFaces, &faceTypes, &faceSizes, &faces);
624:   }
625:   DMSetUp(idm);
626:   /* Initialize cones so we do not need the bash table to tell us that a cone has been set */
627:   {
628:     PetscSection cs;
629:     PetscInt    *cones, csize;

631:     DMPlexGetConeSection(idm, &cs);
632:     DMPlexGetCones(idm, &cones);
633:     PetscSectionGetStorageSize(cs, &csize);
634:     for (c = 0; c < csize; ++c) cones[c] = -1;
635:   }
636:   /* Set cones */
637:   for (d = 0; d <= depth; ++d) {
638:     const PetscInt *cone;
639:     PetscInt        pStart, pEnd, p;

641:     if (d == cellDepth) continue;
642:     DMPlexGetDepthStratum(dm, d, &pStart, &pEnd);
643:     for (p = pStart; p < pEnd; ++p) {
644:       DMPlexGetCone(dm, p, &cone);
645:       DMPlexSetCone(idm, p, cone);
646:       DMPlexGetConeOrientation(dm, p, &cone);
647:       DMPlexSetConeOrientation(idm, p, cone);
648:     }
649:   }
650:   for (c = cStart; c < cEnd; ++c) {
651:     const PetscInt       *cone, *faceSizes, *faces;
652:     const DMPolytopeType *faceTypes;
653:     DMPolytopeType        ct;
654:     PetscInt              numFaces, cf, foff = 0;

656:     DMPlexGetCellType(dm, c, &ct);
657:     DMPlexGetCone(dm, c, &cone);
658:     DMPlexGetRawFaces_Internal(dm, ct, cone, &numFaces, &faceTypes, &faceSizes, &faces);
659:     for (cf = 0; cf < numFaces; foff += faceSizes[cf], ++cf) {
660:       DMPolytopeType   faceType = faceTypes[cf];
661:       const PetscInt   faceSize = faceSizes[cf];
662:       const PetscInt  *face     = &faces[foff];
663:       const PetscInt  *fcone;
664:       PetscHashIJKLKey key;
665:       PetscHashIter    iter;
666:       PetscBool        missing;
667:       PetscInt         f;

670:       key.i = face[0];
671:       key.j = faceSize > 1 ? face[1] : PETSC_MAX_INT;
672:       key.k = faceSize > 2 ? face[2] : PETSC_MAX_INT;
673:       key.l = faceSize > 3 ? face[3] : PETSC_MAX_INT;
674:       PetscSortInt(faceSize, (PetscInt *)&key);
675:       PetscHashIJKLPut(faceTable, key, &iter, &missing);
676:       PetscHashIJKLIterGet(faceTable, iter, &f);
677:       DMPlexInsertCone(idm, c, cf, f);
678:       DMPlexGetCone(idm, f, &fcone);
679:       if (fcone[0] < 0) DMPlexSetCone(idm, f, face);
680:       {
681:         const PetscInt *cone;
682:         PetscInt        coneSize, ornt;

684:         DMPlexGetConeSize(idm, f, &coneSize);
685:         DMPlexGetCone(idm, f, &cone);
687:         /* Notice that we have to use vertices here because the lower dimensional faces have not been created yet */
688:         DMPolytopeGetVertexOrientation(faceType, cone, face, &ornt);
689:         DMPlexInsertConeOrientation(idm, c, cf, ornt);
690:       }
691:     }
692:     DMPlexRestoreRawFaces_Internal(dm, ct, cone, &numFaces, &faceTypes, &faceSizes, &faces);
693:   }
694:   PetscHashIJKLDestroy(&faceTable);
695:   DMPlexSymmetrize(idm);
696:   DMPlexStratify(idm);
697:   return 0;
698: }

700: static PetscErrorCode SortRmineRremoteByRemote_Private(PetscSF sf, PetscInt *rmine1[], PetscInt *rremote1[])
701: {
702:   PetscInt           nleaves;
703:   PetscInt           nranks;
704:   const PetscMPIInt *ranks   = NULL;
705:   const PetscInt    *roffset = NULL, *rmine = NULL, *rremote = NULL;
706:   PetscInt           n, o, r;

708:   PetscSFGetRootRanks(sf, &nranks, &ranks, &roffset, &rmine, &rremote);
709:   nleaves = roffset[nranks];
710:   PetscMalloc2(nleaves, rmine1, nleaves, rremote1);
711:   for (r = 0; r < nranks; r++) {
712:     /* simultaneously sort rank-wise portions of rmine & rremote by values in rremote
713:        - to unify order with the other side */
714:     o = roffset[r];
715:     n = roffset[r + 1] - o;
716:     PetscArraycpy(&(*rmine1)[o], &rmine[o], n);
717:     PetscArraycpy(&(*rremote1)[o], &rremote[o], n);
718:     PetscSortIntWithArray(n, &(*rremote1)[o], &(*rmine1)[o]);
719:   }
720:   return 0;
721: }

723: PetscErrorCode DMPlexOrientInterface_Internal(DM dm)
724: {
725:   PetscSF            sf;
726:   const PetscInt    *locals;
727:   const PetscSFNode *remotes;
728:   const PetscMPIInt *ranks;
729:   const PetscInt    *roffset;
730:   PetscInt          *rmine1, *rremote1; /* rmine and rremote copies simultaneously sorted by rank and rremote */
731:   PetscInt           nroots, p, nleaves, nranks, r, maxConeSize = 0;
732:   PetscInt(*roots)[4], (*leaves)[4], mainCone[4];
733:   PetscMPIInt(*rootsRanks)[4], (*leavesRanks)[4];
734:   MPI_Comm    comm;
735:   PetscMPIInt rank, size;
736:   PetscInt    debug = 0;

738:   PetscObjectGetComm((PetscObject)dm, &comm);
739:   MPI_Comm_rank(comm, &rank);
740:   MPI_Comm_size(comm, &size);
741:   DMGetPointSF(dm, &sf);
742:   DMViewFromOptions(dm, NULL, "-before_orient_interface_dm_view");
743:   if (PetscDefined(USE_DEBUG)) DMPlexCheckPointSF(dm, sf, PETSC_FALSE);
744:   PetscSFGetGraph(sf, &nroots, &nleaves, &locals, &remotes);
745:   if (nroots < 0) return 0;
746:   PetscSFSetUp(sf);
747:   SortRmineRremoteByRemote_Private(sf, &rmine1, &rremote1);
748:   for (p = 0; p < nleaves; ++p) {
749:     PetscInt coneSize;
750:     DMPlexGetConeSize(dm, locals[p], &coneSize);
751:     maxConeSize = PetscMax(maxConeSize, coneSize);
752:   }
754:   PetscMalloc4(nroots, &roots, nroots, &leaves, nroots, &rootsRanks, nroots, &leavesRanks);
755:   for (p = 0; p < nroots; ++p) {
756:     const PetscInt *cone;
757:     PetscInt        coneSize, c, ind0;

759:     DMPlexGetConeSize(dm, p, &coneSize);
760:     DMPlexGetCone(dm, p, &cone);
761:     /* Ignore vertices */
762:     if (coneSize < 2) {
763:       for (c = 0; c < 4; c++) {
764:         roots[p][c]      = -1;
765:         rootsRanks[p][c] = -1;
766:       }
767:       continue;
768:     }
769:     /* Translate all points to root numbering */
770:     for (c = 0; c < PetscMin(coneSize, 4); c++) {
771:       PetscFindInt(cone[c], nleaves, locals, &ind0);
772:       if (ind0 < 0) {
773:         roots[p][c]      = cone[c];
774:         rootsRanks[p][c] = rank;
775:       } else {
776:         roots[p][c]      = remotes[ind0].index;
777:         rootsRanks[p][c] = remotes[ind0].rank;
778:       }
779:     }
780:     for (c = coneSize; c < 4; c++) {
781:       roots[p][c]      = -1;
782:       rootsRanks[p][c] = -1;
783:     }
784:   }
785:   for (p = 0; p < nroots; ++p) {
786:     PetscInt c;
787:     for (c = 0; c < 4; c++) {
788:       leaves[p][c]      = -2;
789:       leavesRanks[p][c] = -2;
790:     }
791:   }
792:   PetscSFBcastBegin(sf, MPIU_4INT, roots, leaves, MPI_REPLACE);
793:   PetscSFBcastBegin(sf, MPI_4INT, rootsRanks, leavesRanks, MPI_REPLACE);
794:   PetscSFBcastEnd(sf, MPIU_4INT, roots, leaves, MPI_REPLACE);
795:   PetscSFBcastEnd(sf, MPI_4INT, rootsRanks, leavesRanks, MPI_REPLACE);
796:   if (debug) {
797:     PetscSynchronizedFlush(comm, NULL);
798:     if (rank == 0) PetscSynchronizedPrintf(comm, "Referenced roots\n");
799:   }
800:   PetscSFGetRootRanks(sf, &nranks, &ranks, &roffset, NULL, NULL);
801:   for (p = 0; p < nroots; ++p) {
802:     DMPolytopeType  ct;
803:     const PetscInt *cone;
804:     PetscInt        coneSize, c, ind0, o;

806:     if (leaves[p][0] < 0) continue; /* Ignore vertices */
807:     DMPlexGetCellType(dm, p, &ct);
808:     DMPlexGetConeSize(dm, p, &coneSize);
809:     DMPlexGetCone(dm, p, &cone);
810:     if (debug) {
811:       PetscSynchronizedPrintf(comm, "[%d]  %4" PetscInt_FMT ": cone=[%4" PetscInt_FMT " %4" PetscInt_FMT " %4" PetscInt_FMT " %4" PetscInt_FMT "] roots=[(%d,%4" PetscInt_FMT ") (%d,%4" PetscInt_FMT ") (%d,%4" PetscInt_FMT ") (%d,%4" PetscInt_FMT ")] leaves=[(%d,%4" PetscInt_FMT ") (%d,%4" PetscInt_FMT ") (%d,%4" PetscInt_FMT ") (%d,%4" PetscInt_FMT ")]", rank, p, cone[0], cone[1], cone[2], cone[3], rootsRanks[p][0], roots[p][0], rootsRanks[p][1], roots[p][1], rootsRanks[p][2], roots[p][2], rootsRanks[p][3], roots[p][3], leavesRanks[p][0], leaves[p][0], leavesRanks[p][1], leaves[p][1], leavesRanks[p][2], leaves[p][2], leavesRanks[p][3], leaves[p][3]);
812:     }
813:     if (leavesRanks[p][0] != rootsRanks[p][0] || leaves[p][0] != roots[p][0] || leavesRanks[p][1] != rootsRanks[p][1] || leaves[p][1] != roots[p][1] || leavesRanks[p][2] != rootsRanks[p][2] || leaves[p][2] != roots[p][2] || leavesRanks[p][3] != rootsRanks[p][3] || leaves[p][3] != roots[p][3]) {
814:       /* Translate these leaves to my cone points; mainCone means desired order p's cone points */
815:       for (c = 0; c < PetscMin(coneSize, 4); ++c) {
816:         PetscInt rS, rN;

818:         if (leavesRanks[p][c] == rank) {
819:           /* A local leaf is just taken as it is */
820:           mainCone[c] = leaves[p][c];
821:           continue;
822:         }
823:         /* Find index of rank leavesRanks[p][c] among remote ranks */
824:         /* No need for PetscMPIIntCast because these integers were originally cast from PetscMPIInt. */
825:         PetscFindMPIInt((PetscMPIInt)leavesRanks[p][c], nranks, ranks, &r);
828:         /* Find point leaves[p][c] among remote points aimed at rank leavesRanks[p][c] */
829:         rS = roffset[r];
830:         rN = roffset[r + 1] - rS;
831:         PetscFindInt(leaves[p][c], rN, &rremote1[rS], &ind0);
833:         /* Get the corresponding local point */
834:         mainCone[c] = rmine1[rS + ind0];
835:       }
836:       if (debug) PetscSynchronizedPrintf(comm, " mainCone=[%4" PetscInt_FMT " %4" PetscInt_FMT " %4" PetscInt_FMT " %4" PetscInt_FMT "]\n", mainCone[0], mainCone[1], mainCone[2], mainCone[3]);
837:       /* Set the desired order of p's cone points and fix orientations accordingly */
838:       DMPolytopeGetOrientation(ct, cone, mainCone, &o);
839:       DMPlexOrientPoint(dm, p, o);
840:     } else if (debug) PetscSynchronizedPrintf(comm, " ==\n");
841:   }
842:   if (debug) {
843:     PetscSynchronizedFlush(comm, NULL);
844:     MPI_Barrier(comm);
845:   }
846:   DMViewFromOptions(dm, NULL, "-after_orient_interface_dm_view");
847:   PetscFree4(roots, leaves, rootsRanks, leavesRanks);
848:   PetscFree2(rmine1, rremote1);
849:   return 0;
850: }

852: static PetscErrorCode IntArrayViewFromOptions(MPI_Comm comm, const char opt[], const char name[], const char idxname[], const char valname[], PetscInt n, const PetscInt a[])
853: {
854:   PetscInt    idx;
855:   PetscMPIInt rank;
856:   PetscBool   flg;

858:   PetscOptionsHasName(NULL, NULL, opt, &flg);
859:   if (!flg) return 0;
860:   MPI_Comm_rank(comm, &rank);
861:   PetscSynchronizedPrintf(comm, "[%d]%s:\n", rank, name);
862:   for (idx = 0; idx < n; ++idx) PetscSynchronizedPrintf(comm, "[%d]%s %" PetscInt_FMT " %s %" PetscInt_FMT "\n", rank, idxname, idx, valname, a[idx]);
863:   PetscSynchronizedFlush(comm, NULL);
864:   return 0;
865: }

867: static PetscErrorCode SFNodeArrayViewFromOptions(MPI_Comm comm, const char opt[], const char name[], const char idxname[], PetscInt n, const PetscSFNode a[])
868: {
869:   PetscInt    idx;
870:   PetscMPIInt rank;
871:   PetscBool   flg;

873:   PetscOptionsHasName(NULL, NULL, opt, &flg);
874:   if (!flg) return 0;
875:   MPI_Comm_rank(comm, &rank);
876:   PetscSynchronizedPrintf(comm, "[%d]%s:\n", rank, name);
877:   if (idxname) {
878:     for (idx = 0; idx < n; ++idx) PetscSynchronizedPrintf(comm, "[%d]%s %" PetscInt_FMT " rank %" PetscInt_FMT " index %" PetscInt_FMT "\n", rank, idxname, idx, a[idx].rank, a[idx].index);
879:   } else {
880:     for (idx = 0; idx < n; ++idx) PetscSynchronizedPrintf(comm, "[%d]rank %" PetscInt_FMT " index %" PetscInt_FMT "\n", rank, a[idx].rank, a[idx].index);
881:   }
882:   PetscSynchronizedFlush(comm, NULL);
883:   return 0;
884: }

886: static PetscErrorCode DMPlexMapToLocalPoint(DM dm, PetscHMapIJ remotehash, PetscSFNode remotePoint, PetscInt *localPoint, PetscBool *mapFailed)
887: {
888:   PetscSF         sf;
889:   const PetscInt *locals;
890:   PetscMPIInt     rank;

892:   MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);
893:   DMGetPointSF(dm, &sf);
894:   PetscSFGetGraph(sf, NULL, NULL, &locals, NULL);
895:   if (mapFailed) *mapFailed = PETSC_FALSE;
896:   if (remotePoint.rank == rank) {
897:     *localPoint = remotePoint.index;
898:   } else {
899:     PetscHashIJKey key;
900:     PetscInt       l;

902:     key.i = remotePoint.index;
903:     key.j = remotePoint.rank;
904:     PetscHMapIJGet(remotehash, key, &l);
905:     if (l >= 0) {
906:       *localPoint = locals[l];
907:     } else if (mapFailed) *mapFailed = PETSC_TRUE;
908:   }
909:   return 0;
910: }

912: static PetscErrorCode DMPlexMapToGlobalPoint(DM dm, PetscInt localPoint, PetscSFNode *remotePoint, PetscBool *mapFailed)
913: {
914:   PetscSF            sf;
915:   const PetscInt    *locals, *rootdegree;
916:   const PetscSFNode *remotes;
917:   PetscInt           Nl, l;
918:   PetscMPIInt        rank;

920:   if (mapFailed) *mapFailed = PETSC_FALSE;
921:   MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);
922:   DMGetPointSF(dm, &sf);
923:   PetscSFGetGraph(sf, NULL, &Nl, &locals, &remotes);
924:   if (Nl < 0) goto owned;
925:   PetscSFComputeDegreeBegin(sf, &rootdegree);
926:   PetscSFComputeDegreeEnd(sf, &rootdegree);
927:   if (rootdegree[localPoint]) goto owned;
928:   PetscFindInt(localPoint, Nl, locals, &l);
929:   if (l < 0) {
930:     if (mapFailed) *mapFailed = PETSC_TRUE;
931:   } else *remotePoint = remotes[l];
932:   return 0;
933: owned:
934:   remotePoint->rank  = rank;
935:   remotePoint->index = localPoint;
936:   return 0;
937: }

939: static PetscErrorCode DMPlexPointIsShared(DM dm, PetscInt p, PetscBool *isShared)
940: {
941:   PetscSF         sf;
942:   const PetscInt *locals, *rootdegree;
943:   PetscInt        Nl, idx;

945:   *isShared = PETSC_FALSE;
946:   DMGetPointSF(dm, &sf);
947:   PetscSFGetGraph(sf, NULL, &Nl, &locals, NULL);
948:   if (Nl < 0) return 0;
949:   PetscFindInt(p, Nl, locals, &idx);
950:   if (idx >= 0) {
951:     *isShared = PETSC_TRUE;
952:     return 0;
953:   }
954:   PetscSFComputeDegreeBegin(sf, &rootdegree);
955:   PetscSFComputeDegreeEnd(sf, &rootdegree);
956:   if (rootdegree[p] > 0) *isShared = PETSC_TRUE;
957:   return 0;
958: }

960: static PetscErrorCode DMPlexConeIsShared(DM dm, PetscInt p, PetscBool *isShared)
961: {
962:   const PetscInt *cone;
963:   PetscInt        coneSize, c;
964:   PetscBool       cShared = PETSC_TRUE;

966:   DMPlexGetConeSize(dm, p, &coneSize);
967:   DMPlexGetCone(dm, p, &cone);
968:   for (c = 0; c < coneSize; ++c) {
969:     PetscBool pointShared;

971:     DMPlexPointIsShared(dm, cone[c], &pointShared);
972:     cShared = (PetscBool)(cShared && pointShared);
973:   }
974:   *isShared = coneSize ? cShared : PETSC_FALSE;
975:   return 0;
976: }

978: static PetscErrorCode DMPlexGetConeMinimum(DM dm, PetscInt p, PetscSFNode *cpmin)
979: {
980:   const PetscInt *cone;
981:   PetscInt        coneSize, c;
982:   PetscSFNode     cmin = {PETSC_MAX_INT, PETSC_MAX_INT}, missing = {-1, -1};

984:   DMPlexGetConeSize(dm, p, &coneSize);
985:   DMPlexGetCone(dm, p, &cone);
986:   for (c = 0; c < coneSize; ++c) {
987:     PetscSFNode rcp;
988:     PetscBool   mapFailed;

990:     DMPlexMapToGlobalPoint(dm, cone[c], &rcp, &mapFailed);
991:     if (mapFailed) {
992:       cmin = missing;
993:     } else {
994:       cmin = (rcp.rank < cmin.rank) || (rcp.rank == cmin.rank && rcp.index < cmin.index) ? rcp : cmin;
995:     }
996:   }
997:   *cpmin = coneSize ? cmin : missing;
998:   return 0;
999: }

1001: /*
1002:   Each shared face has an entry in the candidates array:
1003:     (-1, coneSize-1), {(global cone point)}
1004:   where the set is missing the point p which we use as the key for the face
1005: */
1006: static PetscErrorCode DMPlexAddSharedFace_Private(DM dm, PetscSection candidateSection, PetscSFNode candidates[], PetscHMapIJ faceHash, PetscInt p, PetscBool debug)
1007: {
1008:   MPI_Comm        comm;
1009:   const PetscInt *support;
1010:   PetscInt        supportSize, s, off = 0, idx = 0, overlap, cellHeight, height;
1011:   PetscMPIInt     rank;

1013:   PetscObjectGetComm((PetscObject)dm, &comm);
1014:   MPI_Comm_rank(comm, &rank);
1015:   DMPlexGetOverlap(dm, &overlap);
1016:   DMPlexGetVTKCellHeight(dm, &cellHeight);
1017:   DMPlexGetPointHeight(dm, p, &height);
1018:   if (!overlap && height <= cellHeight + 1) {
1019:     /* cells can't be shared for non-overlapping meshes */
1020:     if (debug) PetscSynchronizedPrintf(comm, "[%d]    Skipping face %" PetscInt_FMT " to avoid adding cell to hashmap since this is nonoverlapping mesh\n", rank, p);
1021:     return 0;
1022:   }
1023:   DMPlexGetSupportSize(dm, p, &supportSize);
1024:   DMPlexGetSupport(dm, p, &support);
1025:   if (candidates) PetscSectionGetOffset(candidateSection, p, &off);
1026:   for (s = 0; s < supportSize; ++s) {
1027:     const PetscInt  face = support[s];
1028:     const PetscInt *cone;
1029:     PetscSFNode     cpmin = {-1, -1}, rp = {-1, -1};
1030:     PetscInt        coneSize, c, f;
1031:     PetscBool       isShared = PETSC_FALSE;
1032:     PetscHashIJKey  key;

1034:     /* Only add point once */
1035:     if (debug) PetscSynchronizedPrintf(comm, "[%d]    Support face %" PetscInt_FMT "\n", rank, face);
1036:     key.i = p;
1037:     key.j = face;
1038:     PetscHMapIJGet(faceHash, key, &f);
1039:     if (f >= 0) continue;
1040:     DMPlexConeIsShared(dm, face, &isShared);
1041:     DMPlexGetConeMinimum(dm, face, &cpmin);
1042:     DMPlexMapToGlobalPoint(dm, p, &rp, NULL);
1043:     if (debug) {
1044:       PetscSynchronizedPrintf(comm, "[%d]      Face point %" PetscInt_FMT " is shared: %d\n", rank, face, (int)isShared);
1045:       PetscSynchronizedPrintf(comm, "[%d]      Global point (%" PetscInt_FMT ", %" PetscInt_FMT ") Min Cone Point (%" PetscInt_FMT ", %" PetscInt_FMT ")\n", rank, rp.rank, rp.index, cpmin.rank, cpmin.index);
1046:     }
1047:     if (isShared && (rp.rank == cpmin.rank && rp.index == cpmin.index)) {
1048:       PetscHMapIJSet(faceHash, key, p);
1049:       if (candidates) {
1050:         if (debug) PetscSynchronizedPrintf(comm, "[%d]    Adding shared face %" PetscInt_FMT " at idx %" PetscInt_FMT "\n[%d]     ", rank, face, idx, rank);
1051:         DMPlexGetConeSize(dm, face, &coneSize);
1052:         DMPlexGetCone(dm, face, &cone);
1053:         candidates[off + idx].rank    = -1;
1054:         candidates[off + idx++].index = coneSize - 1;
1055:         candidates[off + idx].rank    = rank;
1056:         candidates[off + idx++].index = face;
1057:         for (c = 0; c < coneSize; ++c) {
1058:           const PetscInt cp = cone[c];

1060:           if (cp == p) continue;
1061:           DMPlexMapToGlobalPoint(dm, cp, &candidates[off + idx], NULL);
1062:           if (debug) PetscSynchronizedPrintf(comm, " (%" PetscInt_FMT ",%" PetscInt_FMT ")", candidates[off + idx].rank, candidates[off + idx].index);
1063:           ++idx;
1064:         }
1065:         if (debug) PetscSynchronizedPrintf(comm, "\n");
1066:       } else {
1067:         /* Add cone size to section */
1068:         if (debug) PetscSynchronizedPrintf(comm, "[%d]    Scheduling shared face %" PetscInt_FMT "\n", rank, face);
1069:         DMPlexGetConeSize(dm, face, &coneSize);
1070:         PetscHMapIJSet(faceHash, key, p);
1071:         PetscSectionAddDof(candidateSection, p, coneSize + 1);
1072:       }
1073:     }
1074:   }
1075:   return 0;
1076: }

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

1081:   Collective on dm

1083:   Input Parameters:
1084: + dm      - The interpolated DM
1085: - pointSF - The initial SF without interpolated points

1087:   Output Parameter:
1088: . pointSF - The SF including interpolated points

1090:   Level: developer

1092:    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

1094: .seealso: `DMPlexInterpolate()`, `DMPlexUninterpolate()`
1095: @*/
1096: PetscErrorCode DMPlexInterpolatePointSF(DM dm, PetscSF pointSF)
1097: {
1098:   MPI_Comm           comm;
1099:   PetscHMapIJ        remoteHash;
1100:   PetscHMapI         claimshash;
1101:   PetscSection       candidateSection, candidateRemoteSection, claimSection;
1102:   PetscSFNode       *candidates, *candidatesRemote, *claims;
1103:   const PetscInt    *localPoints, *rootdegree;
1104:   const PetscSFNode *remotePoints;
1105:   PetscInt           ov, Nr, r, Nl, l;
1106:   PetscInt           candidatesSize, candidatesRemoteSize, claimsSize;
1107:   PetscBool          flg, debug = PETSC_FALSE;
1108:   PetscMPIInt        rank;

1112:   DMPlexIsDistributed(dm, &flg);
1113:   if (!flg) return 0;
1114:   /* Set initial SF so that lower level queries work */
1115:   DMSetPointSF(dm, pointSF);
1116:   PetscObjectGetComm((PetscObject)dm, &comm);
1117:   MPI_Comm_rank(comm, &rank);
1118:   DMPlexGetOverlap(dm, &ov);
1120:   PetscOptionsHasName(NULL, ((PetscObject)dm)->prefix, "-dmplex_interp_debug", &debug);
1121:   PetscObjectViewFromOptions((PetscObject)dm, NULL, "-dm_interp_pre_view");
1122:   PetscObjectViewFromOptions((PetscObject)pointSF, NULL, "-petscsf_interp_pre_view");
1123:   PetscLogEventBegin(DMPLEX_InterpolateSF, dm, 0, 0, 0);
1124:   /* Step 0: Precalculations */
1125:   PetscSFGetGraph(pointSF, &Nr, &Nl, &localPoints, &remotePoints);
1127:   PetscHMapIJCreate(&remoteHash);
1128:   for (l = 0; l < Nl; ++l) {
1129:     PetscHashIJKey key;
1130:     key.i = remotePoints[l].index;
1131:     key.j = remotePoints[l].rank;
1132:     PetscHMapIJSet(remoteHash, key, l);
1133:   }
1134:   /*   Compute root degree to identify shared points */
1135:   PetscSFComputeDegreeBegin(pointSF, &rootdegree);
1136:   PetscSFComputeDegreeEnd(pointSF, &rootdegree);
1137:   IntArrayViewFromOptions(comm, "-interp_root_degree_view", "Root degree", "point", "degree", Nr, rootdegree);
1138:   /*
1139:   1) Loop over each leaf point $p$ at depth $d$ in the SF
1140:   \item Get set $F(p)$ of faces $f$ in the support of $p$ for which
1141:   \begin{itemize}
1142:     \item all cone points of $f$ are shared
1143:     \item $p$ is the cone point with smallest canonical number
1144:   \end{itemize}
1145:   \item Send $F(p)$ and the cone of each face to the active root point $r(p)$
1146:   \item At the root, if at least two faces with a given cone are present, including a local face, mark the face as shared \label{alg:rootStep} and choose the root face
1147:   \item Send the root face from the root back to all leaf process
1148:   \item Leaf processes add the shared face to the SF
1149:   */
1150:   /* Step 1: Construct section+SFNode array
1151:        The section has entries for all shared faces for which we have a leaf point in the cone
1152:        The array holds candidate shared faces, each face is referred to by the leaf point */
1153:   PetscSectionCreate(comm, &candidateSection);
1154:   PetscSectionSetChart(candidateSection, 0, Nr);
1155:   {
1156:     PetscHMapIJ faceHash;

1158:     PetscHMapIJCreate(&faceHash);
1159:     for (l = 0; l < Nl; ++l) {
1160:       const PetscInt p = localPoints[l];

1162:       if (debug) PetscSynchronizedPrintf(comm, "[%d]  First pass leaf point %" PetscInt_FMT "\n", rank, p);
1163:       DMPlexAddSharedFace_Private(dm, candidateSection, NULL, faceHash, p, debug);
1164:     }
1165:     PetscHMapIJClear(faceHash);
1166:     PetscSectionSetUp(candidateSection);
1167:     PetscSectionGetStorageSize(candidateSection, &candidatesSize);
1168:     PetscMalloc1(candidatesSize, &candidates);
1169:     for (l = 0; l < Nl; ++l) {
1170:       const PetscInt p = localPoints[l];

1172:       if (debug) PetscSynchronizedPrintf(comm, "[%d]  Second pass leaf point %" PetscInt_FMT "\n", rank, p);
1173:       DMPlexAddSharedFace_Private(dm, candidateSection, candidates, faceHash, p, debug);
1174:     }
1175:     PetscHMapIJDestroy(&faceHash);
1176:     if (debug) PetscSynchronizedFlush(comm, NULL);
1177:   }
1178:   PetscObjectSetName((PetscObject)candidateSection, "Candidate Section");
1179:   PetscObjectViewFromOptions((PetscObject)candidateSection, NULL, "-petscsection_interp_candidate_view");
1180:   SFNodeArrayViewFromOptions(comm, "-petscsection_interp_candidate_view", "Candidates", NULL, candidatesSize, candidates);
1181:   /* Step 2: Gather candidate section / array pair into the root partition via inverse(multi(pointSF)). */
1182:   /*   Note that this section is indexed by offsets into leaves, not by point number */
1183:   {
1184:     PetscSF   sfMulti, sfInverse, sfCandidates;
1185:     PetscInt *remoteOffsets;

1187:     PetscSFGetMultiSF(pointSF, &sfMulti);
1188:     PetscSFCreateInverseSF(sfMulti, &sfInverse);
1189:     PetscSectionCreate(comm, &candidateRemoteSection);
1190:     PetscSFDistributeSection(sfInverse, candidateSection, &remoteOffsets, candidateRemoteSection);
1191:     PetscSFCreateSectionSF(sfInverse, candidateSection, remoteOffsets, candidateRemoteSection, &sfCandidates);
1192:     PetscSectionGetStorageSize(candidateRemoteSection, &candidatesRemoteSize);
1193:     PetscMalloc1(candidatesRemoteSize, &candidatesRemote);
1194:     PetscSFBcastBegin(sfCandidates, MPIU_2INT, candidates, candidatesRemote, MPI_REPLACE);
1195:     PetscSFBcastEnd(sfCandidates, MPIU_2INT, candidates, candidatesRemote, MPI_REPLACE);
1196:     PetscSFDestroy(&sfInverse);
1197:     PetscSFDestroy(&sfCandidates);
1198:     PetscFree(remoteOffsets);

1200:     PetscObjectSetName((PetscObject)candidateRemoteSection, "Remote Candidate Section");
1201:     PetscObjectViewFromOptions((PetscObject)candidateRemoteSection, NULL, "-petscsection_interp_candidate_remote_view");
1202:     SFNodeArrayViewFromOptions(comm, "-petscsection_interp_candidate_remote_view", "Remote Candidates", NULL, candidatesRemoteSize, candidatesRemote);
1203:   }
1204:   /* Step 3: At the root, if at least two faces with a given cone are present, including a local face, mark the face as shared and choose the root face */
1205:   {
1206:     PetscHashIJKLRemote faceTable;
1207:     PetscInt            idx, idx2;

1209:     PetscHashIJKLRemoteCreate(&faceTable);
1210:     /* There is a section point for every leaf attached to a given root point */
1211:     for (r = 0, idx = 0, idx2 = 0; r < Nr; ++r) {
1212:       PetscInt deg;

1214:       for (deg = 0; deg < rootdegree[r]; ++deg, ++idx) {
1215:         PetscInt offset, dof, d;

1217:         PetscSectionGetDof(candidateRemoteSection, idx, &dof);
1218:         PetscSectionGetOffset(candidateRemoteSection, idx, &offset);
1219:         /* dof may include many faces from the remote process */
1220:         for (d = 0; d < dof; ++d) {
1221:           const PetscInt         hidx  = offset + d;
1222:           const PetscInt         Np    = candidatesRemote[hidx].index + 1;
1223:           const PetscSFNode      rface = candidatesRemote[hidx + 1];
1224:           const PetscSFNode     *fcone = &candidatesRemote[hidx + 2];
1225:           PetscSFNode            fcp0;
1226:           const PetscSFNode      pmax = {PETSC_MAX_INT, PETSC_MAX_INT};
1227:           const PetscInt        *join = NULL;
1228:           PetscHashIJKLRemoteKey key;
1229:           PetscHashIter          iter;
1230:           PetscBool              missing, mapToLocalPointFailed = PETSC_FALSE;
1231:           PetscInt               points[1024], p, joinSize;

1233:           if (debug)
1234:             PetscCall(PetscSynchronizedPrintf(PetscObjectComm((PetscObject)dm), "[%d]  Checking face (%" PetscInt_FMT ", %" PetscInt_FMT ") at (%" PetscInt_FMT ", %" PetscInt_FMT ", %" PetscInt_FMT ") with cone size %" PetscInt_FMT "\n", rank, rface.rank,
1235:                                               rface.index, r, idx, d, Np));
1237:           fcp0.rank  = rank;
1238:           fcp0.index = r;
1239:           d += Np;
1240:           /* Put remote face in hash table */
1241:           key.i = fcp0;
1242:           key.j = fcone[0];
1243:           key.k = Np > 2 ? fcone[1] : pmax;
1244:           key.l = Np > 3 ? fcone[2] : pmax;
1245:           PetscSortSFNode(Np, (PetscSFNode *)&key);
1246:           PetscHashIJKLRemotePut(faceTable, key, &iter, &missing);
1247:           if (missing) {
1248:             if (debug) PetscSynchronizedPrintf(PetscObjectComm((PetscObject)dm), "[%d]  Setting remote face (%" PetscInt_FMT ", %" PetscInt_FMT ")\n", rank, rface.index, rface.rank);
1249:             PetscHashIJKLRemoteIterSet(faceTable, iter, rface);
1250:           } else {
1251:             PetscSFNode oface;

1253:             PetscHashIJKLRemoteIterGet(faceTable, iter, &oface);
1254:             if ((rface.rank < oface.rank) || (rface.rank == oface.rank && rface.index < oface.index)) {
1255:               if (debug) PetscSynchronizedPrintf(PetscObjectComm((PetscObject)dm), "[%d]  Replacing with remote face (%" PetscInt_FMT ", %" PetscInt_FMT ")\n", rank, rface.index, rface.rank);
1256:               PetscHashIJKLRemoteIterSet(faceTable, iter, rface);
1257:             }
1258:           }
1259:           /* Check for local face */
1260:           points[0] = r;
1261:           for (p = 1; p < Np; ++p) {
1262:             DMPlexMapToLocalPoint(dm, remoteHash, fcone[p - 1], &points[p], &mapToLocalPointFailed);
1263:             if (mapToLocalPointFailed) break; /* We got a point not in our overlap */
1264:             if (debug) PetscSynchronizedPrintf(PetscObjectComm((PetscObject)dm), "[%d]  Checking local candidate %" PetscInt_FMT "\n", rank, points[p]);
1265:           }
1266:           if (mapToLocalPointFailed) continue;
1267:           DMPlexGetJoin(dm, Np, points, &joinSize, &join);
1268:           if (joinSize == 1) {
1269:             PetscSFNode lface;
1270:             PetscSFNode oface;

1272:             /* Always replace with local face */
1273:             lface.rank  = rank;
1274:             lface.index = join[0];
1275:             PetscHashIJKLRemoteIterGet(faceTable, iter, &oface);
1276:             if (debug)
1277:               PetscSynchronizedPrintf(PetscObjectComm((PetscObject)dm), "[%d]  Replacing (%" PetscInt_FMT ", %" PetscInt_FMT ") with local face (%" PetscInt_FMT ", %" PetscInt_FMT ")\n", rank, oface.index, oface.rank, lface.index, lface.rank);
1278:             PetscHashIJKLRemoteIterSet(faceTable, iter, lface);
1279:           }
1280:           DMPlexRestoreJoin(dm, Np, points, &joinSize, &join);
1281:         }
1282:       }
1283:       /* Put back faces for this root */
1284:       for (deg = 0; deg < rootdegree[r]; ++deg, ++idx2) {
1285:         PetscInt offset, dof, d;

1287:         PetscSectionGetDof(candidateRemoteSection, idx2, &dof);
1288:         PetscSectionGetOffset(candidateRemoteSection, idx2, &offset);
1289:         /* dof may include many faces from the remote process */
1290:         for (d = 0; d < dof; ++d) {
1291:           const PetscInt         hidx  = offset + d;
1292:           const PetscInt         Np    = candidatesRemote[hidx].index + 1;
1293:           const PetscSFNode     *fcone = &candidatesRemote[hidx + 2];
1294:           PetscSFNode            fcp0;
1295:           const PetscSFNode      pmax = {PETSC_MAX_INT, PETSC_MAX_INT};
1296:           PetscHashIJKLRemoteKey key;
1297:           PetscHashIter          iter;
1298:           PetscBool              missing;

1300:           if (debug) PetscSynchronizedPrintf(PetscObjectComm((PetscObject)dm), "[%d]  Entering face at (%" PetscInt_FMT ", %" PetscInt_FMT ")\n", rank, r, idx);
1302:           fcp0.rank  = rank;
1303:           fcp0.index = r;
1304:           d += Np;
1305:           /* Find remote face in hash table */
1306:           key.i = fcp0;
1307:           key.j = fcone[0];
1308:           key.k = Np > 2 ? fcone[1] : pmax;
1309:           key.l = Np > 3 ? fcone[2] : pmax;
1310:           PetscSortSFNode(Np, (PetscSFNode *)&key);
1311:           if (debug)
1312:             PetscCall(PetscSynchronizedPrintf(PetscObjectComm((PetscObject)dm), "[%d]    key (%" PetscInt_FMT ", %" PetscInt_FMT ") (%" PetscInt_FMT ", %" PetscInt_FMT ") (%" PetscInt_FMT ", %" PetscInt_FMT ") (%" PetscInt_FMT ", %" PetscInt_FMT ")\n", rank,
1313:                                               key.i.rank, key.i.index, key.j.rank, key.j.index, key.k.rank, key.k.index, key.l.rank, key.l.index));
1314:           PetscHashIJKLRemotePut(faceTable, key, &iter, &missing);
1316:           PetscHashIJKLRemoteIterGet(faceTable, iter, &candidatesRemote[hidx]);
1317:         }
1318:       }
1319:     }
1320:     if (debug) PetscSynchronizedFlush(PetscObjectComm((PetscObject)dm), NULL);
1321:     PetscHashIJKLRemoteDestroy(&faceTable);
1322:   }
1323:   /* Step 4: Push back owned faces */
1324:   {
1325:     PetscSF      sfMulti, sfClaims, sfPointNew;
1326:     PetscSFNode *remotePointsNew;
1327:     PetscInt    *remoteOffsets, *localPointsNew;
1328:     PetscInt     pStart, pEnd, r, NlNew, p;

1330:     /* 4) Push claims back to receiver via the MultiSF and derive new pointSF mapping on receiver */
1331:     PetscSFGetMultiSF(pointSF, &sfMulti);
1332:     PetscSectionCreate(comm, &claimSection);
1333:     PetscSFDistributeSection(sfMulti, candidateRemoteSection, &remoteOffsets, claimSection);
1334:     PetscSFCreateSectionSF(sfMulti, candidateRemoteSection, remoteOffsets, claimSection, &sfClaims);
1335:     PetscSectionGetStorageSize(claimSection, &claimsSize);
1336:     PetscMalloc1(claimsSize, &claims);
1337:     for (p = 0; p < claimsSize; ++p) claims[p].rank = -1;
1338:     PetscSFBcastBegin(sfClaims, MPIU_2INT, candidatesRemote, claims, MPI_REPLACE);
1339:     PetscSFBcastEnd(sfClaims, MPIU_2INT, candidatesRemote, claims, MPI_REPLACE);
1340:     PetscSFDestroy(&sfClaims);
1341:     PetscFree(remoteOffsets);
1342:     PetscObjectSetName((PetscObject)claimSection, "Claim Section");
1343:     PetscObjectViewFromOptions((PetscObject)claimSection, NULL, "-petscsection_interp_claim_view");
1344:     SFNodeArrayViewFromOptions(comm, "-petscsection_interp_claim_view", "Claims", NULL, claimsSize, claims);
1345:     /* Step 5) Walk the original section of local supports and add an SF entry for each updated item */
1346:     /* TODO I should not have to do a join here since I already put the face and its cone in the candidate section */
1347:     PetscHMapICreate(&claimshash);
1348:     for (r = 0; r < Nr; ++r) {
1349:       PetscInt dof, off, d;

1351:       if (debug) PetscSynchronizedPrintf(comm, "[%d]  Checking root for claims %" PetscInt_FMT "\n", rank, r);
1352:       PetscSectionGetDof(candidateSection, r, &dof);
1353:       PetscSectionGetOffset(candidateSection, r, &off);
1354:       for (d = 0; d < dof;) {
1355:         if (claims[off + d].rank >= 0) {
1356:           const PetscInt  faceInd = off + d;
1357:           const PetscInt  Np      = candidates[off + d].index;
1358:           const PetscInt *join    = NULL;
1359:           PetscInt        joinSize, points[1024], c;

1361:           if (debug) PetscSynchronizedPrintf(comm, "[%d]    Found claim for remote point (%" PetscInt_FMT ", %" PetscInt_FMT ")\n", rank, claims[faceInd].rank, claims[faceInd].index);
1362:           points[0] = r;
1363:           if (debug) PetscSynchronizedPrintf(comm, "[%d]      point %" PetscInt_FMT "\n", rank, points[0]);
1364:           for (c = 0, d += 2; c < Np; ++c, ++d) {
1365:             DMPlexMapToLocalPoint(dm, remoteHash, candidates[off + d], &points[c + 1], NULL);
1366:             if (debug) PetscSynchronizedPrintf(comm, "[%d]      point %" PetscInt_FMT "\n", rank, points[c + 1]);
1367:           }
1368:           DMPlexGetJoin(dm, Np + 1, points, &joinSize, &join);
1369:           if (joinSize == 1) {
1370:             if (claims[faceInd].rank == rank) {
1371:               if (debug) PetscSynchronizedPrintf(comm, "[%d]    Ignoring local face %" PetscInt_FMT " for non-remote partner\n", rank, join[0]);
1372:             } else {
1373:               if (debug) PetscSynchronizedPrintf(comm, "[%d]    Found local face %" PetscInt_FMT "\n", rank, join[0]);
1374:               PetscHMapISet(claimshash, join[0], faceInd);
1375:             }
1376:           } else {
1377:             if (debug) PetscSynchronizedPrintf(comm, "[%d]    Failed to find face\n", rank);
1378:           }
1379:           DMPlexRestoreJoin(dm, Np + 1, points, &joinSize, &join);
1380:         } else {
1381:           if (debug) PetscSynchronizedPrintf(comm, "[%d]    No claim for point %" PetscInt_FMT "\n", rank, r);
1382:           d += claims[off + d].index + 1;
1383:         }
1384:       }
1385:     }
1386:     if (debug) PetscSynchronizedFlush(comm, NULL);
1387:     /* Step 6) Create new pointSF from hashed claims */
1388:     PetscHMapIGetSize(claimshash, &NlNew);
1389:     DMPlexGetChart(dm, &pStart, &pEnd);
1390:     PetscMalloc1(Nl + NlNew, &localPointsNew);
1391:     PetscMalloc1(Nl + NlNew, &remotePointsNew);
1392:     for (l = 0; l < Nl; ++l) {
1393:       localPointsNew[l]        = localPoints[l];
1394:       remotePointsNew[l].index = remotePoints[l].index;
1395:       remotePointsNew[l].rank  = remotePoints[l].rank;
1396:     }
1397:     p = Nl;
1398:     PetscHMapIGetKeys(claimshash, &p, localPointsNew);
1399:     /* We sort new points, and assume they are numbered after all existing points */
1400:     PetscSortInt(NlNew, &localPointsNew[Nl]);
1401:     for (p = Nl; p < Nl + NlNew; ++p) {
1402:       PetscInt off;
1403:       PetscHMapIGet(claimshash, localPointsNew[p], &off);
1405:       remotePointsNew[p] = claims[off];
1406:     }
1407:     PetscSFCreate(comm, &sfPointNew);
1408:     PetscSFSetGraph(sfPointNew, pEnd - pStart, Nl + NlNew, localPointsNew, PETSC_OWN_POINTER, remotePointsNew, PETSC_OWN_POINTER);
1409:     PetscSFSetUp(sfPointNew);
1410:     DMSetPointSF(dm, sfPointNew);
1411:     PetscObjectViewFromOptions((PetscObject)sfPointNew, NULL, "-petscsf_interp_view");
1412:     if (PetscDefined(USE_DEBUG)) DMPlexCheckPointSF(dm, sfPointNew, PETSC_FALSE);
1413:     PetscSFDestroy(&sfPointNew);
1414:     PetscHMapIDestroy(&claimshash);
1415:   }
1416:   PetscHMapIJDestroy(&remoteHash);
1417:   PetscSectionDestroy(&candidateSection);
1418:   PetscSectionDestroy(&candidateRemoteSection);
1419:   PetscSectionDestroy(&claimSection);
1420:   PetscFree(candidates);
1421:   PetscFree(candidatesRemote);
1422:   PetscFree(claims);
1423:   PetscLogEventEnd(DMPLEX_InterpolateSF, dm, 0, 0, 0);
1424:   return 0;
1425: }

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

1430:   Collective on dm

1432:   Input Parameters:
1433: + dm - The DMPlex object with only cells and vertices
1434: - dmInt - The interpolated DM

1436:   Output Parameter:
1437: . dmInt - The complete DMPlex object

1439:   Level: intermediate

1441:   Notes:
1442:     Labels and coordinates are copied.

1444:   Developer Notes:
1445:     It sets plex->interpolated = DMPLEX_INTERPOLATED_FULL.

1447: .seealso: `DMPlexUninterpolate()`, `DMPlexCreateFromCellListPetsc()`, `DMPlexCopyCoordinates()`
1448: @*/
1449: PetscErrorCode DMPlexInterpolate(DM dm, DM *dmInt)
1450: {
1451:   DMPlexInterpolatedFlag interpolated;
1452:   DM                     idm, odm = dm;
1453:   PetscSF                sfPoint;
1454:   PetscInt               depth, dim, d;
1455:   const char            *name;
1456:   PetscBool              flg = PETSC_TRUE;

1460:   PetscLogEventBegin(DMPLEX_Interpolate, dm, 0, 0, 0);
1461:   DMPlexGetDepth(dm, &depth);
1462:   DMGetDimension(dm, &dim);
1463:   DMPlexIsInterpolated(dm, &interpolated);
1465:   if (interpolated == DMPLEX_INTERPOLATED_FULL) {
1466:     PetscObjectReference((PetscObject)dm);
1467:     idm = dm;
1468:   } else {
1469:     for (d = 1; d < dim; ++d) {
1470:       /* Create interpolated mesh */
1471:       DMCreate(PetscObjectComm((PetscObject)dm), &idm);
1472:       DMSetType(idm, DMPLEX);
1473:       DMSetDimension(idm, dim);
1474:       if (depth > 0) {
1475:         DMPlexInterpolateFaces_Internal(odm, 1, idm);
1476:         DMGetPointSF(odm, &sfPoint);
1477:         if (PetscDefined(USE_DEBUG)) DMPlexCheckPointSF(odm, sfPoint, PETSC_FALSE);
1478:         {
1479:           /* TODO: We need to systematically fix cases of distributed Plexes with no graph set */
1480:           PetscInt nroots;
1481:           PetscSFGetGraph(sfPoint, &nroots, NULL, NULL, NULL);
1482:           if (nroots >= 0) DMPlexInterpolatePointSF(idm, sfPoint);
1483:         }
1484:       }
1485:       if (odm != dm) DMDestroy(&odm);
1486:       odm = idm;
1487:     }
1488:     PetscObjectGetName((PetscObject)dm, &name);
1489:     PetscObjectSetName((PetscObject)idm, name);
1490:     DMPlexCopyCoordinates(dm, idm);
1491:     DMCopyLabels(dm, idm, PETSC_COPY_VALUES, PETSC_FALSE, DM_COPY_LABELS_FAIL);
1492:     PetscOptionsGetBool(((PetscObject)dm)->options, ((PetscObject)dm)->prefix, "-dm_plex_interpolate_orient_interfaces", &flg, NULL);
1493:     if (flg) DMPlexOrientInterface_Internal(idm);
1494:   }
1495:   /* This function makes the mesh fully interpolated on all ranks */
1496:   {
1497:     DM_Plex *plex      = (DM_Plex *)idm->data;
1498:     plex->interpolated = plex->interpolatedCollective = DMPLEX_INTERPOLATED_FULL;
1499:   }
1500:   DMPlexCopy_Internal(dm, PETSC_TRUE, PETSC_TRUE, idm);
1501:   *dmInt = idm;
1502:   PetscLogEventEnd(DMPLEX_Interpolate, dm, 0, 0, 0);
1503:   return 0;
1504: }

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

1509:   Collective on dmA

1511:   Input Parameter:
1512: . dmA - The DMPlex object with initial coordinates

1514:   Output Parameter:
1515: . dmB - The DMPlex object with copied coordinates

1517:   Level: intermediate

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

1521: .seealso: `DMCopyLabels()`, `DMGetCoordinates()`, `DMGetCoordinatesLocal()`, `DMGetCoordinateDM()`, `DMGetCoordinateSection()`
1522: @*/
1523: PetscErrorCode DMPlexCopyCoordinates(DM dmA, DM dmB)
1524: {
1525:   Vec          coordinatesA, coordinatesB;
1526:   VecType      vtype;
1527:   PetscSection coordSectionA, coordSectionB;
1528:   PetscScalar *coordsA, *coordsB;
1529:   PetscInt     spaceDim, Nf, vStartA, vStartB, vEndA, vEndB, coordSizeB, v, d;
1530:   PetscInt     cStartA, cEndA, cStartB, cEndB, cS, cE, cdim;
1531:   PetscBool    lc = PETSC_FALSE;

1535:   if (dmA == dmB) return 0;
1536:   DMGetCoordinateDim(dmA, &cdim);
1537:   DMSetCoordinateDim(dmB, cdim);
1538:   DMPlexGetDepthStratum(dmA, 0, &vStartA, &vEndA);
1539:   DMPlexGetDepthStratum(dmB, 0, &vStartB, &vEndB);
1541:   /* Copy over discretization if it exists */
1542:   {
1543:     DM                 cdmA, cdmB;
1544:     PetscDS            dsA, dsB;
1545:     PetscObject        objA, objB;
1546:     PetscClassId       idA, idB;
1547:     const PetscScalar *constants;
1548:     PetscInt           cdim, Nc;

1550:     DMGetCoordinateDM(dmA, &cdmA);
1551:     DMGetCoordinateDM(dmB, &cdmB);
1552:     DMGetField(cdmA, 0, NULL, &objA);
1553:     DMGetField(cdmB, 0, NULL, &objB);
1554:     PetscObjectGetClassId(objA, &idA);
1555:     PetscObjectGetClassId(objB, &idB);
1556:     if ((idA == PETSCFE_CLASSID) && (idA != idB)) {
1557:       DMSetField(cdmB, 0, NULL, objA);
1558:       DMCreateDS(cdmB);
1559:       DMGetDS(cdmA, &dsA);
1560:       DMGetDS(cdmB, &dsB);
1561:       PetscDSGetCoordinateDimension(dsA, &cdim);
1562:       PetscDSSetCoordinateDimension(dsB, cdim);
1563:       PetscDSGetConstants(dsA, &Nc, &constants);
1564:       PetscDSSetConstants(dsB, Nc, (PetscScalar *)constants);
1565:     }
1566:   }
1567:   DMPlexGetHeightStratum(dmA, 0, &cStartA, &cEndA);
1568:   DMPlexGetHeightStratum(dmB, 0, &cStartB, &cEndB);
1569:   DMGetCoordinateSection(dmA, &coordSectionA);
1570:   DMGetCoordinateSection(dmB, &coordSectionB);
1571:   if (coordSectionA == coordSectionB) return 0;
1572:   PetscSectionGetNumFields(coordSectionA, &Nf);
1573:   if (!Nf) return 0;
1575:   if (!coordSectionB) {
1576:     PetscInt dim;

1578:     PetscSectionCreate(PetscObjectComm((PetscObject)coordSectionA), &coordSectionB);
1579:     DMGetCoordinateDim(dmA, &dim);
1580:     DMSetCoordinateSection(dmB, dim, coordSectionB);
1581:     PetscObjectDereference((PetscObject)coordSectionB);
1582:   }
1583:   PetscSectionSetNumFields(coordSectionB, 1);
1584:   PetscSectionGetFieldComponents(coordSectionA, 0, &spaceDim);
1585:   PetscSectionSetFieldComponents(coordSectionB, 0, spaceDim);
1586:   PetscSectionGetChart(coordSectionA, &cS, &cE);
1587:   if (cStartA <= cS && cS < cEndA) { /* localized coordinates */
1589:     cS = cS - cStartA + cStartB;
1590:     cE = vEndB;
1591:     lc = PETSC_TRUE;
1592:   } else {
1593:     cS = vStartB;
1594:     cE = vEndB;
1595:   }
1596:   PetscSectionSetChart(coordSectionB, cS, cE);
1597:   for (v = vStartB; v < vEndB; ++v) {
1598:     PetscSectionSetDof(coordSectionB, v, spaceDim);
1599:     PetscSectionSetFieldDof(coordSectionB, v, 0, spaceDim);
1600:   }
1601:   if (lc) { /* localized coordinates */
1602:     PetscInt c;

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

1607:       PetscSectionGetDof(coordSectionA, c + cStartA, &dof);
1608:       PetscSectionSetDof(coordSectionB, c + cStartB, dof);
1609:       PetscSectionSetFieldDof(coordSectionB, c + cStartB, 0, dof);
1610:     }
1611:   }
1612:   PetscSectionSetUp(coordSectionB);
1613:   PetscSectionGetStorageSize(coordSectionB, &coordSizeB);
1614:   DMGetCoordinatesLocal(dmA, &coordinatesA);
1615:   VecCreate(PETSC_COMM_SELF, &coordinatesB);
1616:   PetscObjectSetName((PetscObject)coordinatesB, "coordinates");
1617:   VecSetSizes(coordinatesB, coordSizeB, PETSC_DETERMINE);
1618:   VecGetBlockSize(coordinatesA, &d);
1619:   VecSetBlockSize(coordinatesB, d);
1620:   VecGetType(coordinatesA, &vtype);
1621:   VecSetType(coordinatesB, vtype);
1622:   VecGetArray(coordinatesA, &coordsA);
1623:   VecGetArray(coordinatesB, &coordsB);
1624:   for (v = 0; v < vEndB - vStartB; ++v) {
1625:     PetscInt offA, offB;

1627:     PetscSectionGetOffset(coordSectionA, v + vStartA, &offA);
1628:     PetscSectionGetOffset(coordSectionB, v + vStartB, &offB);
1629:     for (d = 0; d < spaceDim; ++d) coordsB[offB + d] = coordsA[offA + d];
1630:   }
1631:   if (lc) { /* localized coordinates */
1632:     PetscInt c;

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

1637:       PetscSectionGetOffset(coordSectionA, c + cStartA, &offA);
1638:       PetscSectionGetOffset(coordSectionB, c + cStartB, &offB);
1639:       PetscSectionGetDof(coordSectionA, c + cStartA, &dof);
1640:       PetscArraycpy(coordsB + offB, coordsA + offA, dof);
1641:     }
1642:   }
1643:   VecRestoreArray(coordinatesA, &coordsA);
1644:   VecRestoreArray(coordinatesB, &coordsB);
1645:   DMSetCoordinatesLocal(dmB, coordinatesB);
1646:   VecDestroy(&coordinatesB);
1647:   return 0;
1648: }

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

1653:   Collective on dm

1655:   Input Parameter:
1656: . dm - The complete DMPlex object

1658:   Output Parameter:
1659: . dmUnint - The DMPlex object with only cells and vertices

1661:   Level: intermediate

1663:   Notes:
1664:     It does not copy over the coordinates.

1666:   Developer Notes:
1667:     It sets plex->interpolated = DMPLEX_INTERPOLATED_NONE.

1669: .seealso: `DMPlexInterpolate()`, `DMPlexCreateFromCellListPetsc()`, `DMPlexCopyCoordinates()`
1670: @*/
1671: PetscErrorCode DMPlexUninterpolate(DM dm, DM *dmUnint)
1672: {
1673:   DMPlexInterpolatedFlag interpolated;
1674:   DM                     udm;
1675:   PetscInt               dim, vStart, vEnd, cStart, cEnd, c, maxConeSize = 0, *cone;

1679:   DMGetDimension(dm, &dim);
1680:   DMPlexIsInterpolated(dm, &interpolated);
1682:   if (interpolated == DMPLEX_INTERPOLATED_NONE || dim <= 1) {
1683:     /* in case dim <= 1 just keep the DMPLEX_INTERPOLATED_FULL flag */
1684:     PetscObjectReference((PetscObject)dm);
1685:     *dmUnint = dm;
1686:     return 0;
1687:   }
1688:   DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
1689:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
1690:   DMCreate(PetscObjectComm((PetscObject)dm), &udm);
1691:   DMSetType(udm, DMPLEX);
1692:   DMSetDimension(udm, dim);
1693:   DMPlexSetChart(udm, cStart, vEnd);
1694:   for (c = cStart; c < cEnd; ++c) {
1695:     PetscInt *closure = NULL, closureSize, cl, coneSize = 0;

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

1701:       if ((p >= vStart) && (p < vEnd)) ++coneSize;
1702:     }
1703:     DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
1704:     DMPlexSetConeSize(udm, c, coneSize);
1705:     maxConeSize = PetscMax(maxConeSize, coneSize);
1706:   }
1707:   DMSetUp(udm);
1708:   PetscMalloc1(maxConeSize, &cone);
1709:   for (c = cStart; c < cEnd; ++c) {
1710:     PetscInt *closure = NULL, closureSize, cl, coneSize = 0;

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

1716:       if ((p >= vStart) && (p < vEnd)) cone[coneSize++] = p;
1717:     }
1718:     DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
1719:     DMPlexSetCone(udm, c, cone);
1720:   }
1721:   PetscFree(cone);
1722:   DMPlexSymmetrize(udm);
1723:   DMPlexStratify(udm);
1724:   /* Reduce SF */
1725:   {
1726:     PetscSF            sfPoint, sfPointUn;
1727:     const PetscSFNode *remotePoints;
1728:     const PetscInt    *localPoints;
1729:     PetscSFNode       *remotePointsUn;
1730:     PetscInt          *localPointsUn;
1731:     PetscInt           numRoots, numLeaves, l;
1732:     PetscInt           numLeavesUn = 0, n = 0;

1734:     /* Get original SF information */
1735:     DMGetPointSF(dm, &sfPoint);
1736:     if (PetscDefined(USE_DEBUG)) DMPlexCheckPointSF(dm, sfPoint, PETSC_FALSE);
1737:     DMGetPointSF(udm, &sfPointUn);
1738:     PetscSFGetGraph(sfPoint, &numRoots, &numLeaves, &localPoints, &remotePoints);
1739:     if (numRoots >= 0) {
1740:       /* Allocate space for cells and vertices */
1741:       for (l = 0; l < numLeaves; ++l) {
1742:         const PetscInt p = localPoints[l];

1744:         if ((vStart <= p && p < vEnd) || (cStart <= p && p < cEnd)) numLeavesUn++;
1745:       }
1746:       /* Fill in leaves */
1747:       PetscMalloc1(numLeavesUn, &remotePointsUn);
1748:       PetscMalloc1(numLeavesUn, &localPointsUn);
1749:       for (l = 0; l < numLeaves; l++) {
1750:         const PetscInt p = localPoints[l];

1752:         if ((vStart <= p && p < vEnd) || (cStart <= p && p < cEnd)) {
1753:           localPointsUn[n]        = p;
1754:           remotePointsUn[n].rank  = remotePoints[l].rank;
1755:           remotePointsUn[n].index = remotePoints[l].index;
1756:           ++n;
1757:         }
1758:       }
1760:       PetscSFSetGraph(sfPointUn, cEnd - cStart + vEnd - vStart, numLeavesUn, localPointsUn, PETSC_OWN_POINTER, remotePointsUn, PETSC_OWN_POINTER);
1761:     }
1762:   }
1763:   /* This function makes the mesh fully uninterpolated on all ranks */
1764:   {
1765:     DM_Plex *plex      = (DM_Plex *)udm->data;
1766:     plex->interpolated = plex->interpolatedCollective = DMPLEX_INTERPOLATED_NONE;
1767:   }
1768:   DMPlexCopy_Internal(dm, PETSC_TRUE, PETSC_TRUE, udm);
1769:   if (PetscDefined(USE_DEBUG)) DMPlexCheckPointSF(udm, NULL, PETSC_FALSE);
1770:   *dmUnint = udm;
1771:   return 0;
1772: }

1774: static PetscErrorCode DMPlexIsInterpolated_Internal(DM dm, DMPlexInterpolatedFlag *interpolated)
1775: {
1776:   PetscInt coneSize, depth, dim, h, p, pStart, pEnd;
1777:   MPI_Comm comm;

1779:   PetscObjectGetComm((PetscObject)dm, &comm);
1780:   DMPlexGetDepth(dm, &depth);
1781:   DMGetDimension(dm, &dim);

1783:   if (depth == dim) {
1784:     *interpolated = DMPLEX_INTERPOLATED_FULL;
1785:     if (!dim) goto finish;

1787:     /* Check points at height = dim are vertices (have no cones) */
1788:     DMPlexGetHeightStratum(dm, dim, &pStart, &pEnd);
1789:     for (p = pStart; p < pEnd; p++) {
1790:       DMPlexGetConeSize(dm, p, &coneSize);
1791:       if (coneSize) {
1792:         *interpolated = DMPLEX_INTERPOLATED_PARTIAL;
1793:         goto finish;
1794:       }
1795:     }

1797:     /* Check points at height < dim have cones */
1798:     for (h = 0; h < dim; h++) {
1799:       DMPlexGetHeightStratum(dm, h, &pStart, &pEnd);
1800:       for (p = pStart; p < pEnd; p++) {
1801:         DMPlexGetConeSize(dm, p, &coneSize);
1802:         if (!coneSize) {
1803:           *interpolated = DMPLEX_INTERPOLATED_PARTIAL;
1804:           goto finish;
1805:         }
1806:       }
1807:     }
1808:   } else if (depth == 1) {
1809:     *interpolated = DMPLEX_INTERPOLATED_NONE;
1810:   } else {
1811:     *interpolated = DMPLEX_INTERPOLATED_PARTIAL;
1812:   }
1813: finish:
1814:   return 0;
1815: }

1817: /*@
1818:   DMPlexIsInterpolated - Find out to what extent the DMPlex is topologically interpolated.

1820:   Not Collective

1822:   Input Parameter:
1823: . dm      - The DM object

1825:   Output Parameter:
1826: . interpolated - Flag whether the DM is interpolated

1828:   Level: intermediate

1830:   Notes:
1831:   Unlike DMPlexIsInterpolatedCollective(), this is NOT collective
1832:   so the results can be different on different ranks in special cases.
1833:   However, DMPlexInterpolate() guarantees the result is the same on all.

1835:   Unlike DMPlexIsInterpolatedCollective(), this cannot return DMPLEX_INTERPOLATED_MIXED.

1837:   Developer Notes:
1838:   Initially, plex->interpolated = DMPLEX_INTERPOLATED_INVALID.

1840:   If plex->interpolated == DMPLEX_INTERPOLATED_INVALID, DMPlexIsInterpolated_Internal() is called.
1841:   It checks the actual topology and sets plex->interpolated on each rank separately to one of
1842:   DMPLEX_INTERPOLATED_NONE, DMPLEX_INTERPOLATED_PARTIAL or DMPLEX_INTERPOLATED_FULL.

1844:   If plex->interpolated != DMPLEX_INTERPOLATED_INVALID, this function just returns plex->interpolated.

1846:   DMPlexInterpolate() sets plex->interpolated = DMPLEX_INTERPOLATED_FULL,
1847:   and DMPlexUninterpolate() sets plex->interpolated = DMPLEX_INTERPOLATED_NONE.

1849: .seealso: `DMPlexInterpolate()`, `DMPlexIsInterpolatedCollective()`
1850: @*/
1851: PetscErrorCode DMPlexIsInterpolated(DM dm, DMPlexInterpolatedFlag *interpolated)
1852: {
1853:   DM_Plex *plex = (DM_Plex *)dm->data;

1857:   if (plex->interpolated < 0) {
1858:     DMPlexIsInterpolated_Internal(dm, &plex->interpolated);
1859:   } else if (PetscDefined(USE_DEBUG)) {
1860:     DMPlexInterpolatedFlag flg;

1862:     DMPlexIsInterpolated_Internal(dm, &flg);
1864:   }
1865:   *interpolated = plex->interpolated;
1866:   return 0;
1867: }

1869: /*@
1870:   DMPlexIsInterpolatedCollective - Find out to what extent the DMPlex is topologically interpolated (in collective manner).

1872:   Collective

1874:   Input Parameter:
1875: . dm      - The DM object

1877:   Output Parameter:
1878: . interpolated - Flag whether the DM is interpolated

1880:   Level: intermediate

1882:   Notes:
1883:   Unlike DMPlexIsInterpolated(), this is collective so the results are guaranteed to be the same on all ranks.

1885:   This function will return DMPLEX_INTERPOLATED_MIXED if the results of DMPlexIsInterpolated() are different on different ranks.

1887:   Developer Notes:
1888:   Initially, plex->interpolatedCollective = DMPLEX_INTERPOLATED_INVALID.

1890:   If plex->interpolatedCollective == DMPLEX_INTERPOLATED_INVALID, this function calls DMPlexIsInterpolated() which sets plex->interpolated.
1891:   MPI_Allreduce() is then called and collectively consistent flag plex->interpolatedCollective is set and returned;
1892:   if plex->interpolated varies on different ranks, plex->interpolatedCollective = DMPLEX_INTERPOLATED_MIXED,
1893:   otherwise sets plex->interpolatedCollective = plex->interpolated.

1895:   If plex->interpolatedCollective != DMPLEX_INTERPOLATED_INVALID, this function just returns plex->interpolatedCollective.

1897: .seealso: `DMPlexInterpolate()`, `DMPlexIsInterpolated()`
1898: @*/
1899: PetscErrorCode DMPlexIsInterpolatedCollective(DM dm, DMPlexInterpolatedFlag *interpolated)
1900: {
1901:   DM_Plex  *plex  = (DM_Plex *)dm->data;
1902:   PetscBool debug = PETSC_FALSE;

1906:   PetscOptionsGetBool(((PetscObject)dm)->options, ((PetscObject)dm)->prefix, "-dm_plex_is_interpolated_collective_debug", &debug, NULL);
1907:   if (plex->interpolatedCollective < 0) {
1908:     DMPlexInterpolatedFlag min, max;
1909:     MPI_Comm               comm;

1911:     PetscObjectGetComm((PetscObject)dm, &comm);
1912:     DMPlexIsInterpolated(dm, &plex->interpolatedCollective);
1913:     MPI_Allreduce(&plex->interpolatedCollective, &min, 1, MPIU_ENUM, MPI_MIN, comm);
1914:     MPI_Allreduce(&plex->interpolatedCollective, &max, 1, MPIU_ENUM, MPI_MAX, comm);
1915:     if (min != max) plex->interpolatedCollective = DMPLEX_INTERPOLATED_MIXED;
1916:     if (debug) {
1917:       PetscMPIInt rank;

1919:       MPI_Comm_rank(comm, &rank);
1920:       PetscSynchronizedPrintf(comm, "[%d] interpolated=%s interpolatedCollective=%s\n", rank, DMPlexInterpolatedFlags[plex->interpolated], DMPlexInterpolatedFlags[plex->interpolatedCollective]);
1921:       PetscSynchronizedFlush(comm, PETSC_STDOUT);
1922:     }
1923:   }
1924:   *interpolated = plex->interpolatedCollective;
1925:   return 0;
1926: }