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: /* HMapIJKL */

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

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

 13: typedef struct _PetscHMapIJKLRemoteKey {
 14:   PetscSFNode i, j, k, l;
 15: } PetscHMapIJKLRemoteKey;

 17: #define PetscHMapIJKLRemoteKeyHash(key) \
 18:   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)))

 20: #define PetscHMapIJKLRemoteKeyEqual(k1, k2) \
 21:   (((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)

 23: PetscDisableStaticAnalyzerForExpressionUnderstandingThatThisIsDangerousAndBugprone(PETSC_HASH_MAP(HMapIJKLRemote, PetscHMapIJKLRemoteKey, PetscSFNode, PetscHMapIJKLRemoteKeyHash, PetscHMapIJKLRemoteKeyEqual, _PetscInvalidSFNode))

 25:   static PetscErrorCode PetscSortSFNode(PetscInt n, PetscSFNode A[])
 26: {
 27:   PetscInt i;

 29:   PetscFunctionBegin;
 30:   for (i = 1; i < n; ++i) {
 31:     PetscSFNode x = A[i];
 32:     PetscInt    j;

 34:     for (j = i - 1; j >= 0; --j) {
 35:       if ((A[j].rank > x.rank) || (A[j].rank == x.rank && A[j].index > x.index)) break;
 36:       A[j + 1] = A[j];
 37:     }
 38:     A[j + 1] = x;
 39:   }
 40:   PetscFunctionReturn(PETSC_SUCCESS);
 41: }

 43: /*
 44:   DMPlexGetRawFaces_Internal - Gets groups of vertices that correspond to faces for the given cone
 45: */
 46: PetscErrorCode DMPlexGetRawFaces_Internal(DM dm, DMPolytopeType ct, const PetscInt cone[], PetscInt *numFaces, const DMPolytopeType *faceTypes[], const PetscInt *faceSizes[], const PetscInt *faces[])
 47: {
 48:   DMPolytopeType *typesTmp = NULL;
 49:   PetscInt       *sizesTmp = NULL, *facesTmp = NULL;
 50:   PetscInt       *tmp;
 51:   PetscInt        maxConeSize, maxSupportSize, maxSize;
 52:   PetscInt        getSize = 0;

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

477: PetscErrorCode DMPlexRestoreRawFaces_Internal(DM dm, DMPolytopeType ct, const PetscInt cone[], PetscInt *numFaces, const DMPolytopeType *faceTypes[], const PetscInt *faceSizes[], const PetscInt *faces[])
478: {
479:   PetscFunctionBegin;
480:   if (faceTypes) PetscCall(DMRestoreWorkArray(dm, 0, MPIU_INT, (void *)faceTypes));
481:   else if (faceSizes) PetscCall(DMRestoreWorkArray(dm, 0, MPIU_INT, (void *)faceSizes));
482:   else if (faces) PetscCall(DMRestoreWorkArray(dm, 0, MPIU_INT, (void *)faces));
483:   if (faceTypes) *faceTypes = NULL;
484:   if (faceSizes) *faceSizes = NULL;
485:   if (faces) *faces = NULL;
486:   PetscFunctionReturn(PETSC_SUCCESS);
487: }

489: /* This interpolates faces for cells at some stratum */
490: static PetscErrorCode DMPlexInterpolateFaces_Internal(DM dm, PetscInt cellDepth, DM idm)
491: {
492:   DMLabel       ctLabel;
493:   PetscHMapIJKL faceTable;
494:   PetscInt      faceTypeNum[DM_NUM_POLYTOPES];
495:   PetscInt      depth, d, pStart, Np, cStart, cEnd, c, fStart, fEnd;
496:   PetscInt      cntFaces, *facesId, minCone;

498:   PetscFunctionBegin;
499:   PetscCall(DMPlexGetDepth(dm, &depth));
500:   PetscCall(PetscHMapIJKLCreate(&faceTable));
501:   PetscCall(PetscArrayzero(faceTypeNum, DM_NUM_POLYTOPES));
502:   PetscCall(DMPlexGetDepthStratum(dm, cellDepth, &cStart, &cEnd));
503:   /* Number new faces and save face vertices in hash table */
504:   PetscCall(DMPlexGetDepthStratum(dm, depth > cellDepth ? cellDepth : 0, NULL, &fStart));
505:   fEnd = fStart;

507:   minCone = PETSC_MAX_INT;
508:   for (c = cStart, cntFaces = 0; c < cEnd; ++c) {
509:     const PetscInt *cone;
510:     DMPolytopeType  ct;
511:     PetscInt        numFaces = 0, coneSize;

513:     PetscCall(DMPlexGetCellType(dm, c, &ct));
514:     PetscCall(DMPlexGetCone(dm, c, &cone));
515:     PetscCall(DMPlexGetConeSize(dm, c, &coneSize));
516:     for (PetscInt j = 0; j < coneSize; j++) minCone = PetscMin(cone[j], minCone);
517:     PetscCall(DMPlexGetRawFaces_Internal(dm, ct, cone, &numFaces, NULL, NULL, NULL));
518:     cntFaces += numFaces;
519:   }
520:   // Encode so that we can use 0 as an excluded value, instead of PETSC_MAX_INT
521:   minCone = -(minCone - 1);

523:   PetscCall(PetscMalloc1(cntFaces, &facesId));

525:   for (c = cStart, cntFaces = 0; c < cEnd; ++c) {
526:     const PetscInt       *cone, *faceSizes, *faces;
527:     const DMPolytopeType *faceTypes;
528:     DMPolytopeType        ct;
529:     PetscInt              numFaces, cf, foff = 0;

531:     PetscCall(DMPlexGetCellType(dm, c, &ct));
532:     PetscCall(DMPlexGetCone(dm, c, &cone));
533:     PetscCall(DMPlexGetRawFaces_Internal(dm, ct, cone, &numFaces, &faceTypes, &faceSizes, &faces));
534:     for (cf = 0; cf < numFaces; foff += faceSizes[cf], ++cf) {
535:       const PetscInt       faceSize = faceSizes[cf];
536:       const DMPolytopeType faceType = faceTypes[cf];
537:       const PetscInt      *face     = &faces[foff];
538:       PetscHashIJKLKey     key;
539:       PetscHashIter        iter;
540:       PetscBool            missing;

542:       PetscCheck(faceSize <= 4, PETSC_COMM_SELF, PETSC_ERR_SUP, "Do not support faces of size %" PetscInt_FMT " > 4", faceSize);
543:       key.i = face[0] + minCone;
544:       key.j = faceSize > 1 ? face[1] + minCone : 0;
545:       key.k = faceSize > 2 ? face[2] + minCone : 0;
546:       key.l = faceSize > 3 ? face[3] + minCone : 0;
547:       PetscCall(PetscSortInt(faceSize, (PetscInt *)&key));
548:       PetscCall(PetscHMapIJKLPut(faceTable, key, &iter, &missing));
549:       if (missing) {
550:         facesId[cntFaces] = fEnd;
551:         PetscCall(PetscHMapIJKLIterSet(faceTable, iter, fEnd++));
552:         ++faceTypeNum[faceType];
553:       } else PetscCall(PetscHMapIJKLIterGet(faceTable, iter, &facesId[cntFaces]));
554:       cntFaces++;
555:     }
556:     PetscCall(DMPlexRestoreRawFaces_Internal(dm, ct, cone, &numFaces, &faceTypes, &faceSizes, &faces));
557:   }
558:   /* We need to number faces contiguously among types */
559:   {
560:     PetscInt faceTypeStart[DM_NUM_POLYTOPES], ct, numFT = 0;

562:     for (ct = 0; ct < DM_NUM_POLYTOPES; ++ct) {
563:       if (faceTypeNum[ct]) ++numFT;
564:       faceTypeStart[ct] = 0;
565:     }
566:     if (numFT > 1) {
567:       PetscCall(PetscHMapIJKLClear(faceTable));
568:       faceTypeStart[0] = fStart;
569:       for (ct = 1; ct < DM_NUM_POLYTOPES; ++ct) faceTypeStart[ct] = faceTypeStart[ct - 1] + faceTypeNum[ct - 1];
570:       for (c = cStart, cntFaces = 0; c < cEnd; ++c) {
571:         const PetscInt       *cone, *faceSizes, *faces;
572:         const DMPolytopeType *faceTypes;
573:         DMPolytopeType        ct;
574:         PetscInt              numFaces, cf, foff = 0;

576:         PetscCall(DMPlexGetCellType(dm, c, &ct));
577:         PetscCall(DMPlexGetCone(dm, c, &cone));
578:         PetscCall(DMPlexGetRawFaces_Internal(dm, ct, cone, &numFaces, &faceTypes, &faceSizes, &faces));
579:         for (cf = 0; cf < numFaces; foff += faceSizes[cf], ++cf) {
580:           const PetscInt       faceSize = faceSizes[cf];
581:           const DMPolytopeType faceType = faceTypes[cf];
582:           const PetscInt      *face     = &faces[foff];
583:           PetscHashIJKLKey     key;
584:           PetscHashIter        iter;
585:           PetscBool            missing;

587:           key.i = face[0] + minCone;
588:           key.j = faceSize > 1 ? face[1] + minCone : 0;
589:           key.k = faceSize > 2 ? face[2] + minCone : 0;
590:           key.l = faceSize > 3 ? face[3] + minCone : 0;
591:           PetscCall(PetscSortInt(faceSize, (PetscInt *)&key));
592:           PetscCall(PetscHMapIJKLPut(faceTable, key, &iter, &missing));
593:           if (missing) {
594:             facesId[cntFaces] = faceTypeStart[faceType];
595:             PetscCall(PetscHMapIJKLIterSet(faceTable, iter, faceTypeStart[faceType]++));
596:           } else PetscCall(PetscHMapIJKLIterGet(faceTable, iter, &facesId[cntFaces]));
597:           cntFaces++;
598:         }
599:         PetscCall(DMPlexRestoreRawFaces_Internal(dm, ct, cone, &numFaces, &faceTypes, &faceSizes, &faces));
600:       }
601:       for (ct = 1; ct < DM_NUM_POLYTOPES; ++ct) {
602:         PetscCheck(faceTypeStart[ct] == faceTypeStart[ct - 1] + faceTypeNum[ct], PETSC_COMM_SELF, PETSC_ERR_PLIB, "Inconsistent numbering for cell type %s, %" PetscInt_FMT " != %" PetscInt_FMT " + %" PetscInt_FMT, DMPolytopeTypes[ct], faceTypeStart[ct], faceTypeStart[ct - 1], faceTypeNum[ct]);
603:       }
604:     }
605:   }
606:   PetscCall(PetscHMapIJKLDestroy(&faceTable));

608:   /* Add new points, always at the end of the numbering */
609:   PetscCall(DMPlexGetChart(dm, &pStart, &Np));
610:   PetscCall(DMPlexSetChart(idm, pStart, Np + (fEnd - fStart)));
611:   /* Set cone sizes */
612:   /*   Must create the celltype label here so that we do not automatically try to compute the types */
613:   PetscCall(DMCreateLabel(idm, "celltype"));
614:   PetscCall(DMPlexGetCellTypeLabel(idm, &ctLabel));
615:   for (d = 0; d <= depth; ++d) {
616:     DMPolytopeType ct;
617:     PetscInt       coneSize, pStart, pEnd, p;

619:     if (d == cellDepth) continue;
620:     PetscCall(DMPlexGetDepthStratum(dm, d, &pStart, &pEnd));
621:     for (p = pStart; p < pEnd; ++p) {
622:       PetscCall(DMPlexGetConeSize(dm, p, &coneSize));
623:       PetscCall(DMPlexSetConeSize(idm, p, coneSize));
624:       PetscCall(DMPlexGetCellType(dm, p, &ct));
625:       PetscCall(DMPlexSetCellType(idm, p, ct));
626:     }
627:   }
628:   for (c = cStart, cntFaces = 0; c < cEnd; ++c) {
629:     const PetscInt       *cone, *faceSizes;
630:     const DMPolytopeType *faceTypes;
631:     DMPolytopeType        ct;
632:     PetscInt              numFaces, cf;

634:     PetscCall(DMPlexGetCellType(dm, c, &ct));
635:     PetscCall(DMPlexGetCone(dm, c, &cone));
636:     PetscCall(DMPlexGetRawFaces_Internal(dm, ct, cone, &numFaces, &faceTypes, &faceSizes, NULL));
637:     PetscCall(DMPlexSetCellType(idm, c, ct));
638:     PetscCall(DMPlexSetConeSize(idm, c, numFaces));
639:     for (cf = 0; cf < numFaces; ++cf) {
640:       const PetscInt f        = facesId[cntFaces];
641:       DMPolytopeType faceType = faceTypes[cf];
642:       const PetscInt faceSize = faceSizes[cf];
643:       PetscCall(DMPlexSetConeSize(idm, f, faceSize));
644:       PetscCall(DMPlexSetCellType(idm, f, faceType));
645:       cntFaces++;
646:     }
647:     PetscCall(DMPlexRestoreRawFaces_Internal(dm, ct, cone, &numFaces, &faceTypes, &faceSizes, NULL));
648:   }
649:   PetscCall(DMSetUp(idm));
650:   /* Initialize cones so we do not need the bash table to tell us that a cone has been set */
651:   {
652:     PetscSection cs;
653:     PetscInt    *cones, csize;

655:     PetscCall(DMPlexGetConeSection(idm, &cs));
656:     PetscCall(DMPlexGetCones(idm, &cones));
657:     PetscCall(PetscSectionGetStorageSize(cs, &csize));
658:     for (c = 0; c < csize; ++c) cones[c] = -1;
659:   }
660:   /* Set cones */
661:   for (d = 0; d <= depth; ++d) {
662:     const PetscInt *cone;
663:     PetscInt        pStart, pEnd, p;

665:     if (d == cellDepth) continue;
666:     PetscCall(DMPlexGetDepthStratum(dm, d, &pStart, &pEnd));
667:     for (p = pStart; p < pEnd; ++p) {
668:       PetscCall(DMPlexGetCone(dm, p, &cone));
669:       PetscCall(DMPlexSetCone(idm, p, cone));
670:       PetscCall(DMPlexGetConeOrientation(dm, p, &cone));
671:       PetscCall(DMPlexSetConeOrientation(idm, p, cone));
672:     }
673:   }
674:   for (c = cStart, cntFaces = 0; c < cEnd; ++c) {
675:     const PetscInt       *cone, *faceSizes, *faces;
676:     const DMPolytopeType *faceTypes;
677:     DMPolytopeType        ct;
678:     PetscInt              numFaces, cf, foff = 0;

680:     PetscCall(DMPlexGetCellType(dm, c, &ct));
681:     PetscCall(DMPlexGetCone(dm, c, &cone));
682:     PetscCall(DMPlexGetRawFaces_Internal(dm, ct, cone, &numFaces, &faceTypes, &faceSizes, &faces));
683:     for (cf = 0; cf < numFaces; foff += faceSizes[cf], ++cf) {
684:       DMPolytopeType  faceType = faceTypes[cf];
685:       const PetscInt  faceSize = faceSizes[cf];
686:       const PetscInt  f        = facesId[cntFaces];
687:       const PetscInt *face     = &faces[foff];
688:       const PetscInt *fcone;

690:       PetscCall(DMPlexInsertCone(idm, c, cf, f));
691:       PetscCall(DMPlexGetCone(idm, f, &fcone));
692:       if (fcone[0] < 0) PetscCall(DMPlexSetCone(idm, f, face));
693:       {
694:         const PetscInt *cone;
695:         PetscInt        coneSize, ornt;

697:         PetscCall(DMPlexGetConeSize(idm, f, &coneSize));
698:         PetscCall(DMPlexGetCone(idm, f, &cone));
699:         PetscCheck(coneSize == faceSize, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid number of face vertices %" PetscInt_FMT " for face %" PetscInt_FMT " should be %" PetscInt_FMT, coneSize, f, faceSize);
700:         /* Notice that we have to use vertices here because the lower dimensional faces have not been created yet */
701:         PetscCall(DMPolytopeGetVertexOrientation(faceType, cone, face, &ornt));
702:         PetscCall(DMPlexInsertConeOrientation(idm, c, cf, ornt));
703:       }
704:       cntFaces++;
705:     }
706:     PetscCall(DMPlexRestoreRawFaces_Internal(dm, ct, cone, &numFaces, &faceTypes, &faceSizes, &faces));
707:   }
708:   PetscCall(PetscFree(facesId));
709:   PetscCall(DMPlexSymmetrize(idm));
710:   PetscCall(DMPlexStratify(idm));
711:   PetscFunctionReturn(PETSC_SUCCESS);
712: }

714: static PetscErrorCode SortRmineRremoteByRemote_Private(PetscSF sf, PetscInt *rmine1[], PetscInt *rremote1[])
715: {
716:   PetscInt           nleaves;
717:   PetscInt           nranks;
718:   const PetscMPIInt *ranks   = NULL;
719:   const PetscInt    *roffset = NULL, *rmine = NULL, *rremote = NULL;
720:   PetscInt           n, o, r;

722:   PetscFunctionBegin;
723:   PetscCall(PetscSFGetRootRanks(sf, &nranks, &ranks, &roffset, &rmine, &rremote));
724:   nleaves = roffset[nranks];
725:   PetscCall(PetscMalloc2(nleaves, rmine1, nleaves, rremote1));
726:   for (r = 0; r < nranks; r++) {
727:     /* simultaneously sort rank-wise portions of rmine & rremote by values in rremote
728:        - to unify order with the other side */
729:     o = roffset[r];
730:     n = roffset[r + 1] - o;
731:     PetscCall(PetscArraycpy(&(*rmine1)[o], &rmine[o], n));
732:     PetscCall(PetscArraycpy(&(*rremote1)[o], &rremote[o], n));
733:     PetscCall(PetscSortIntWithArray(n, &(*rremote1)[o], &(*rmine1)[o]));
734:   }
735:   PetscFunctionReturn(PETSC_SUCCESS);
736: }

738: PetscErrorCode DMPlexOrientInterface_Internal(DM dm)
739: {
740:   PetscSF            sf;
741:   const PetscInt    *locals;
742:   const PetscSFNode *remotes;
743:   const PetscMPIInt *ranks;
744:   const PetscInt    *roffset;
745:   PetscInt          *rmine1, *rremote1; /* rmine and rremote copies simultaneously sorted by rank and rremote */
746:   PetscInt           nroots, p, nleaves, nranks, r, maxConeSize = 0;
747:   PetscInt(*roots)[4], (*leaves)[4], mainCone[4];
748:   PetscMPIInt(*rootsRanks)[4], (*leavesRanks)[4];
749:   MPI_Comm    comm;
750:   PetscMPIInt rank, size;
751:   PetscInt    debug = 0;

753:   PetscFunctionBegin;
754:   PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
755:   PetscCallMPI(MPI_Comm_rank(comm, &rank));
756:   PetscCallMPI(MPI_Comm_size(comm, &size));
757:   PetscCall(DMGetPointSF(dm, &sf));
758:   PetscCall(DMViewFromOptions(dm, NULL, "-before_orient_interface_dm_view"));
759:   if (PetscDefined(USE_DEBUG)) PetscCall(DMPlexCheckPointSF(dm, sf, PETSC_FALSE));
760:   PetscCall(PetscSFGetGraph(sf, &nroots, &nleaves, &locals, &remotes));
761:   if (nroots < 0) PetscFunctionReturn(PETSC_SUCCESS);
762:   PetscCall(PetscSFSetUp(sf));
763:   PetscCall(SortRmineRremoteByRemote_Private(sf, &rmine1, &rremote1));
764:   for (p = 0; p < nleaves; ++p) {
765:     PetscInt coneSize;
766:     PetscCall(DMPlexGetConeSize(dm, locals[p], &coneSize));
767:     maxConeSize = PetscMax(maxConeSize, coneSize);
768:   }
769:   PetscCheck(maxConeSize <= 4, comm, PETSC_ERR_SUP, "This method does not support cones of size %" PetscInt_FMT, maxConeSize);
770:   PetscCall(PetscMalloc4(nroots, &roots, nroots, &leaves, nroots, &rootsRanks, nroots, &leavesRanks));
771:   for (p = 0; p < nroots; ++p) {
772:     const PetscInt *cone;
773:     PetscInt        coneSize, c, ind0;

775:     PetscCall(DMPlexGetConeSize(dm, p, &coneSize));
776:     PetscCall(DMPlexGetCone(dm, p, &cone));
777:     /* Ignore vertices */
778:     if (coneSize < 2) {
779:       for (c = 0; c < 4; c++) {
780:         roots[p][c]      = -1;
781:         rootsRanks[p][c] = -1;
782:       }
783:       continue;
784:     }
785:     /* Translate all points to root numbering */
786:     for (c = 0; c < PetscMin(coneSize, 4); c++) {
787:       PetscCall(PetscFindInt(cone[c], nleaves, locals, &ind0));
788:       if (ind0 < 0) {
789:         roots[p][c]      = cone[c];
790:         rootsRanks[p][c] = rank;
791:       } else {
792:         roots[p][c]      = remotes[ind0].index;
793:         rootsRanks[p][c] = remotes[ind0].rank;
794:       }
795:     }
796:     for (c = coneSize; c < 4; c++) {
797:       roots[p][c]      = -1;
798:       rootsRanks[p][c] = -1;
799:     }
800:   }
801:   for (p = 0; p < nroots; ++p) {
802:     PetscInt c;
803:     for (c = 0; c < 4; c++) {
804:       leaves[p][c]      = -2;
805:       leavesRanks[p][c] = -2;
806:     }
807:   }
808:   PetscCall(PetscSFBcastBegin(sf, MPIU_4INT, roots, leaves, MPI_REPLACE));
809:   PetscCall(PetscSFBcastBegin(sf, MPI_4INT, rootsRanks, leavesRanks, MPI_REPLACE));
810:   PetscCall(PetscSFBcastEnd(sf, MPIU_4INT, roots, leaves, MPI_REPLACE));
811:   PetscCall(PetscSFBcastEnd(sf, MPI_4INT, rootsRanks, leavesRanks, MPI_REPLACE));
812:   if (debug) {
813:     PetscCall(PetscSynchronizedFlush(comm, NULL));
814:     if (rank == 0) PetscCall(PetscSynchronizedPrintf(comm, "Referenced roots\n"));
815:   }
816:   PetscCall(PetscSFGetRootRanks(sf, &nranks, &ranks, &roffset, NULL, NULL));
817:   for (p = 0; p < nroots; ++p) {
818:     DMPolytopeType  ct;
819:     const PetscInt *cone;
820:     PetscInt        coneSize, c, ind0, o;

822:     if (leaves[p][0] < 0) continue; /* Ignore vertices */
823:     PetscCall(DMPlexGetCellType(dm, p, &ct));
824:     PetscCall(DMPlexGetConeSize(dm, p, &coneSize));
825:     PetscCall(DMPlexGetCone(dm, p, &cone));
826:     if (debug) {
827:       PetscCall(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]));
828:     }
829:     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]) {
830:       /* Translate these leaves to my cone points; mainCone means desired order p's cone points */
831:       for (c = 0; c < PetscMin(coneSize, 4); ++c) {
832:         PetscInt rS, rN;

834:         if (leavesRanks[p][c] == rank) {
835:           /* A local leaf is just taken as it is */
836:           mainCone[c] = leaves[p][c];
837:           continue;
838:         }
839:         /* Find index of rank leavesRanks[p][c] among remote ranks */
840:         /* No need for PetscMPIIntCast because these integers were originally cast from PetscMPIInt. */
841:         PetscCall(PetscFindMPIInt((PetscMPIInt)leavesRanks[p][c], nranks, ranks, &r));
842:         PetscCheck(r >= 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Point %" PetscInt_FMT " cone[%" PetscInt_FMT "]=%" PetscInt_FMT " root (%d,%" PetscInt_FMT ") leaf (%d,%" PetscInt_FMT "): leaf rank not found among remote ranks", p, c, cone[c], rootsRanks[p][c], roots[p][c], leavesRanks[p][c], leaves[p][c]);
843:         PetscCheck(ranks[r] >= 0 && ranks[r] < size, PETSC_COMM_SELF, PETSC_ERR_PLIB, "p=%" PetscInt_FMT " c=%" PetscInt_FMT " commsize=%d: ranks[%" PetscInt_FMT "] = %d makes no sense", p, c, size, r, ranks[r]);
844:         /* Find point leaves[p][c] among remote points aimed at rank leavesRanks[p][c] */
845:         rS = roffset[r];
846:         rN = roffset[r + 1] - rS;
847:         PetscCall(PetscFindInt(leaves[p][c], rN, &rremote1[rS], &ind0));
848:         PetscCheck(ind0 >= 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Point %" PetscInt_FMT " cone[%" PetscInt_FMT "]=%" PetscInt_FMT " root (%d,%" PetscInt_FMT ") leave (%d,%" PetscInt_FMT "): corresponding remote point not found - it seems there is missing connection in point SF!", p, c, cone[c], rootsRanks[p][c], roots[p][c], leavesRanks[p][c], leaves[p][c]);
849:         /* Get the corresponding local point */
850:         mainCone[c] = rmine1[rS + ind0];
851:       }
852:       if (debug) PetscCall(PetscSynchronizedPrintf(comm, " mainCone=[%4" PetscInt_FMT " %4" PetscInt_FMT " %4" PetscInt_FMT " %4" PetscInt_FMT "]\n", mainCone[0], mainCone[1], mainCone[2], mainCone[3]));
853:       /* Set the desired order of p's cone points and fix orientations accordingly */
854:       PetscCall(DMPolytopeGetOrientation(ct, cone, mainCone, &o));
855:       PetscCall(DMPlexOrientPoint(dm, p, o));
856:     } else if (debug) PetscCall(PetscSynchronizedPrintf(comm, " ==\n"));
857:   }
858:   if (debug) {
859:     PetscCall(PetscSynchronizedFlush(comm, NULL));
860:     PetscCallMPI(MPI_Barrier(comm));
861:   }
862:   PetscCall(DMViewFromOptions(dm, NULL, "-after_orient_interface_dm_view"));
863:   PetscCall(PetscFree4(roots, leaves, rootsRanks, leavesRanks));
864:   PetscCall(PetscFree2(rmine1, rremote1));
865:   PetscFunctionReturn(PETSC_SUCCESS);
866: }

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

874:   PetscFunctionBegin;
875:   PetscCall(PetscOptionsHasName(NULL, NULL, opt, &flg));
876:   if (!flg) PetscFunctionReturn(PETSC_SUCCESS);
877:   PetscCallMPI(MPI_Comm_rank(comm, &rank));
878:   PetscCall(PetscSynchronizedPrintf(comm, "[%d]%s:\n", rank, name));
879:   for (idx = 0; idx < n; ++idx) PetscCall(PetscSynchronizedPrintf(comm, "[%d]%s %" PetscInt_FMT " %s %" PetscInt_FMT "\n", rank, idxname, idx, valname, a[idx]));
880:   PetscCall(PetscSynchronizedFlush(comm, NULL));
881:   PetscFunctionReturn(PETSC_SUCCESS);
882: }

884: static PetscErrorCode SFNodeArrayViewFromOptions(MPI_Comm comm, const char opt[], const char name[], const char idxname[], PetscInt n, const PetscSFNode a[])
885: {
886:   PetscInt    idx;
887:   PetscMPIInt rank;
888:   PetscBool   flg;

890:   PetscFunctionBegin;
891:   PetscCall(PetscOptionsHasName(NULL, NULL, opt, &flg));
892:   if (!flg) PetscFunctionReturn(PETSC_SUCCESS);
893:   PetscCallMPI(MPI_Comm_rank(comm, &rank));
894:   PetscCall(PetscSynchronizedPrintf(comm, "[%d]%s:\n", rank, name));
895:   if (idxname) {
896:     for (idx = 0; idx < n; ++idx) PetscCall(PetscSynchronizedPrintf(comm, "[%d]%s %" PetscInt_FMT " rank %" PetscInt_FMT " index %" PetscInt_FMT "\n", rank, idxname, idx, a[idx].rank, a[idx].index));
897:   } else {
898:     for (idx = 0; idx < n; ++idx) PetscCall(PetscSynchronizedPrintf(comm, "[%d]rank %" PetscInt_FMT " index %" PetscInt_FMT "\n", rank, a[idx].rank, a[idx].index));
899:   }
900:   PetscCall(PetscSynchronizedFlush(comm, NULL));
901:   PetscFunctionReturn(PETSC_SUCCESS);
902: }

904: static PetscErrorCode DMPlexMapToLocalPoint(DM dm, PetscHMapIJ remotehash, PetscSFNode remotePoint, PetscInt *localPoint, PetscBool *mapFailed)
905: {
906:   PetscSF         sf;
907:   const PetscInt *locals;
908:   PetscMPIInt     rank;

910:   PetscFunctionBegin;
911:   PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank));
912:   PetscCall(DMGetPointSF(dm, &sf));
913:   PetscCall(PetscSFGetGraph(sf, NULL, NULL, &locals, NULL));
914:   if (mapFailed) *mapFailed = PETSC_FALSE;
915:   if (remotePoint.rank == rank) {
916:     *localPoint = remotePoint.index;
917:   } else {
918:     PetscHashIJKey key;
919:     PetscInt       l;

921:     key.i = remotePoint.index;
922:     key.j = remotePoint.rank;
923:     PetscCall(PetscHMapIJGet(remotehash, key, &l));
924:     if (l >= 0) {
925:       *localPoint = locals[l];
926:     } else if (mapFailed) *mapFailed = PETSC_TRUE;
927:   }
928:   PetscFunctionReturn(PETSC_SUCCESS);
929: }

931: static PetscErrorCode DMPlexMapToGlobalPoint(DM dm, PetscInt localPoint, PetscSFNode *remotePoint, PetscBool *mapFailed)
932: {
933:   PetscSF            sf;
934:   const PetscInt    *locals, *rootdegree;
935:   const PetscSFNode *remotes;
936:   PetscInt           Nl, l;
937:   PetscMPIInt        rank;

939:   PetscFunctionBegin;
940:   if (mapFailed) *mapFailed = PETSC_FALSE;
941:   PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank));
942:   PetscCall(DMGetPointSF(dm, &sf));
943:   PetscCall(PetscSFGetGraph(sf, NULL, &Nl, &locals, &remotes));
944:   if (Nl < 0) goto owned;
945:   PetscCall(PetscSFComputeDegreeBegin(sf, &rootdegree));
946:   PetscCall(PetscSFComputeDegreeEnd(sf, &rootdegree));
947:   if (rootdegree[localPoint]) goto owned;
948:   PetscCall(PetscFindInt(localPoint, Nl, locals, &l));
949:   if (l < 0) {
950:     if (mapFailed) *mapFailed = PETSC_TRUE;
951:   } else *remotePoint = remotes[l];
952:   PetscFunctionReturn(PETSC_SUCCESS);
953: owned:
954:   remotePoint->rank  = rank;
955:   remotePoint->index = localPoint;
956:   PetscFunctionReturn(PETSC_SUCCESS);
957: }

959: static PetscErrorCode DMPlexPointIsShared(DM dm, PetscInt p, PetscBool *isShared)
960: {
961:   PetscSF         sf;
962:   const PetscInt *locals, *rootdegree;
963:   PetscInt        Nl, idx;

965:   PetscFunctionBegin;
966:   *isShared = PETSC_FALSE;
967:   PetscCall(DMGetPointSF(dm, &sf));
968:   PetscCall(PetscSFGetGraph(sf, NULL, &Nl, &locals, NULL));
969:   if (Nl < 0) PetscFunctionReturn(PETSC_SUCCESS);
970:   PetscCall(PetscFindInt(p, Nl, locals, &idx));
971:   if (idx >= 0) {
972:     *isShared = PETSC_TRUE;
973:     PetscFunctionReturn(PETSC_SUCCESS);
974:   }
975:   PetscCall(PetscSFComputeDegreeBegin(sf, &rootdegree));
976:   PetscCall(PetscSFComputeDegreeEnd(sf, &rootdegree));
977:   if (rootdegree[p] > 0) *isShared = PETSC_TRUE;
978:   PetscFunctionReturn(PETSC_SUCCESS);
979: }

981: static PetscErrorCode DMPlexConeIsShared(DM dm, PetscInt p, PetscBool *isShared)
982: {
983:   const PetscInt *cone;
984:   PetscInt        coneSize, c;
985:   PetscBool       cShared = PETSC_TRUE;

987:   PetscFunctionBegin;
988:   PetscCall(DMPlexGetConeSize(dm, p, &coneSize));
989:   PetscCall(DMPlexGetCone(dm, p, &cone));
990:   for (c = 0; c < coneSize; ++c) {
991:     PetscBool pointShared;

993:     PetscCall(DMPlexPointIsShared(dm, cone[c], &pointShared));
994:     cShared = (PetscBool)(cShared && pointShared);
995:   }
996:   *isShared = coneSize ? cShared : PETSC_FALSE;
997:   PetscFunctionReturn(PETSC_SUCCESS);
998: }

1000: static PetscErrorCode DMPlexGetConeMinimum(DM dm, PetscInt p, PetscSFNode *cpmin)
1001: {
1002:   const PetscInt *cone;
1003:   PetscInt        coneSize, c;
1004:   PetscSFNode     cmin = {PETSC_MAX_INT, PETSC_MAX_INT}, missing = {-1, -1};

1006:   PetscFunctionBegin;
1007:   PetscCall(DMPlexGetConeSize(dm, p, &coneSize));
1008:   PetscCall(DMPlexGetCone(dm, p, &cone));
1009:   for (c = 0; c < coneSize; ++c) {
1010:     PetscSFNode rcp;
1011:     PetscBool   mapFailed;

1013:     PetscCall(DMPlexMapToGlobalPoint(dm, cone[c], &rcp, &mapFailed));
1014:     if (mapFailed) {
1015:       cmin = missing;
1016:     } else {
1017:       cmin = (rcp.rank < cmin.rank) || (rcp.rank == cmin.rank && rcp.index < cmin.index) ? rcp : cmin;
1018:     }
1019:   }
1020:   *cpmin = coneSize ? cmin : missing;
1021:   PetscFunctionReturn(PETSC_SUCCESS);
1022: }

1024: /*
1025:   Each shared face has an entry in the candidates array:
1026:     (-1, coneSize-1), {(global cone point)}
1027:   where the set is missing the point p which we use as the key for the face
1028: */
1029: static PetscErrorCode DMPlexAddSharedFace_Private(DM dm, PetscSection candidateSection, PetscSFNode candidates[], PetscHMapIJ faceHash, PetscInt p, PetscBool debug)
1030: {
1031:   MPI_Comm        comm;
1032:   const PetscInt *support;
1033:   PetscInt        supportSize, s, off = 0, idx = 0, overlap, cellHeight, height;
1034:   PetscMPIInt     rank;

1036:   PetscFunctionBegin;
1037:   PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
1038:   PetscCallMPI(MPI_Comm_rank(comm, &rank));
1039:   PetscCall(DMPlexGetOverlap(dm, &overlap));
1040:   PetscCall(DMPlexGetVTKCellHeight(dm, &cellHeight));
1041:   PetscCall(DMPlexGetPointHeight(dm, p, &height));
1042:   if (!overlap && height <= cellHeight + 1) {
1043:     /* cells can't be shared for non-overlapping meshes */
1044:     if (debug) PetscCall(PetscSynchronizedPrintf(comm, "[%d]    Skipping face %" PetscInt_FMT " to avoid adding cell to hashmap since this is nonoverlapping mesh\n", rank, p));
1045:     PetscFunctionReturn(PETSC_SUCCESS);
1046:   }
1047:   PetscCall(DMPlexGetSupportSize(dm, p, &supportSize));
1048:   PetscCall(DMPlexGetSupport(dm, p, &support));
1049:   if (candidates) PetscCall(PetscSectionGetOffset(candidateSection, p, &off));
1050:   for (s = 0; s < supportSize; ++s) {
1051:     const PetscInt  face = support[s];
1052:     const PetscInt *cone;
1053:     PetscSFNode     cpmin = {-1, -1}, rp = {-1, -1};
1054:     PetscInt        coneSize, c, f;
1055:     PetscBool       isShared = PETSC_FALSE;
1056:     PetscHashIJKey  key;

1058:     /* Only add point once */
1059:     if (debug) PetscCall(PetscSynchronizedPrintf(comm, "[%d]    Support face %" PetscInt_FMT "\n", rank, face));
1060:     key.i = p;
1061:     key.j = face;
1062:     PetscCall(PetscHMapIJGet(faceHash, key, &f));
1063:     if (f >= 0) continue;
1064:     PetscCall(DMPlexConeIsShared(dm, face, &isShared));
1065:     PetscCall(DMPlexGetConeMinimum(dm, face, &cpmin));
1066:     PetscCall(DMPlexMapToGlobalPoint(dm, p, &rp, NULL));
1067:     if (debug) {
1068:       PetscCall(PetscSynchronizedPrintf(comm, "[%d]      Face point %" PetscInt_FMT " is shared: %d\n", rank, face, (int)isShared));
1069:       PetscCall(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));
1070:     }
1071:     if (isShared && (rp.rank == cpmin.rank && rp.index == cpmin.index)) {
1072:       PetscCall(PetscHMapIJSet(faceHash, key, p));
1073:       if (candidates) {
1074:         if (debug) PetscCall(PetscSynchronizedPrintf(comm, "[%d]    Adding shared face %" PetscInt_FMT " at idx %" PetscInt_FMT "\n[%d]     ", rank, face, idx, rank));
1075:         PetscCall(DMPlexGetConeSize(dm, face, &coneSize));
1076:         PetscCall(DMPlexGetCone(dm, face, &cone));
1077:         candidates[off + idx].rank    = -1;
1078:         candidates[off + idx++].index = coneSize - 1;
1079:         candidates[off + idx].rank    = rank;
1080:         candidates[off + idx++].index = face;
1081:         for (c = 0; c < coneSize; ++c) {
1082:           const PetscInt cp = cone[c];

1084:           if (cp == p) continue;
1085:           PetscCall(DMPlexMapToGlobalPoint(dm, cp, &candidates[off + idx], NULL));
1086:           if (debug) PetscCall(PetscSynchronizedPrintf(comm, " (%" PetscInt_FMT ",%" PetscInt_FMT ")", candidates[off + idx].rank, candidates[off + idx].index));
1087:           ++idx;
1088:         }
1089:         if (debug) PetscCall(PetscSynchronizedPrintf(comm, "\n"));
1090:       } else {
1091:         /* Add cone size to section */
1092:         if (debug) PetscCall(PetscSynchronizedPrintf(comm, "[%d]    Scheduling shared face %" PetscInt_FMT "\n", rank, face));
1093:         PetscCall(DMPlexGetConeSize(dm, face, &coneSize));
1094:         PetscCall(PetscHMapIJSet(faceHash, key, p));
1095:         PetscCall(PetscSectionAddDof(candidateSection, p, coneSize + 1));
1096:       }
1097:     }
1098:   }
1099:   PetscFunctionReturn(PETSC_SUCCESS);
1100: }

1102: /*@
1103:   DMPlexInterpolatePointSF - Insert interpolated points in the overlap into the `PointSF` in parallel, following local interpolation

1105:   Collective

1107:   Input Parameters:
1108: + dm      - The interpolated `DMPLEX`
1109: - pointSF - The initial `PetscSF` without interpolated points

1111:   Level: developer

1113:   Note:
1114:   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`

1116: .seealso: `DMPLEX`, `DMPlexInterpolate()`, `DMPlexUninterpolate()`
1117: @*/
1118: PetscErrorCode DMPlexInterpolatePointSF(DM dm, PetscSF pointSF)
1119: {
1120:   MPI_Comm           comm;
1121:   PetscHMapIJ        remoteHash;
1122:   PetscHMapI         claimshash;
1123:   PetscSection       candidateSection, candidateRemoteSection, claimSection;
1124:   PetscSFNode       *candidates, *candidatesRemote, *claims;
1125:   const PetscInt    *localPoints, *rootdegree;
1126:   const PetscSFNode *remotePoints;
1127:   PetscInt           ov, Nr, r, Nl, l;
1128:   PetscInt           candidatesSize, candidatesRemoteSize, claimsSize;
1129:   PetscBool          flg, debug = PETSC_FALSE;
1130:   PetscMPIInt        rank;

1132:   PetscFunctionBegin;
1135:   PetscCall(DMPlexIsDistributed(dm, &flg));
1136:   if (!flg) PetscFunctionReturn(PETSC_SUCCESS);
1137:   /* Set initial SF so that lower level queries work */
1138:   PetscCall(DMSetPointSF(dm, pointSF));
1139:   PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
1140:   PetscCallMPI(MPI_Comm_rank(comm, &rank));
1141:   PetscCall(DMPlexGetOverlap(dm, &ov));
1142:   PetscCheck(!ov, comm, PETSC_ERR_SUP, "Interpolation of overlapped DMPlex not implemented yet");
1143:   PetscCall(PetscOptionsHasName(NULL, ((PetscObject)dm)->prefix, "-dmplex_interp_debug", &debug));
1144:   PetscCall(PetscObjectViewFromOptions((PetscObject)dm, NULL, "-dm_interp_pre_view"));
1145:   PetscCall(PetscObjectViewFromOptions((PetscObject)pointSF, NULL, "-petscsf_interp_pre_view"));
1146:   PetscCall(PetscLogEventBegin(DMPLEX_InterpolateSF, dm, 0, 0, 0));
1147:   /* Step 0: Precalculations */
1148:   PetscCall(PetscSFGetGraph(pointSF, &Nr, &Nl, &localPoints, &remotePoints));
1149:   PetscCheck(Nr >= 0, comm, PETSC_ERR_ARG_WRONGSTATE, "This DMPlex is distributed but input PointSF has no graph set");
1150:   PetscCall(PetscHMapIJCreate(&remoteHash));
1151:   for (l = 0; l < Nl; ++l) {
1152:     PetscHashIJKey key;
1153:     key.i = remotePoints[l].index;
1154:     key.j = remotePoints[l].rank;
1155:     PetscCall(PetscHMapIJSet(remoteHash, key, l));
1156:   }
1157:   /*   Compute root degree to identify shared points */
1158:   PetscCall(PetscSFComputeDegreeBegin(pointSF, &rootdegree));
1159:   PetscCall(PetscSFComputeDegreeEnd(pointSF, &rootdegree));
1160:   PetscCall(IntArrayViewFromOptions(comm, "-interp_root_degree_view", "Root degree", "point", "degree", Nr, rootdegree));
1161:   /*
1162:   1) Loop over each leaf point $p$ at depth $d$ in the SF
1163:   \item Get set $F(p)$ of faces $f$ in the support of $p$ for which
1164:   \begin{itemize}
1165:     \item all cone points of $f$ are shared
1166:     \item $p$ is the cone point with smallest canonical number
1167:   \end{itemize}
1168:   \item Send $F(p)$ and the cone of each face to the active root point $r(p)$
1169:   \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
1170:   \item Send the root face from the root back to all leaf process
1171:   \item Leaf processes add the shared face to the SF
1172:   */
1173:   /* Step 1: Construct section+SFNode array
1174:        The section has entries for all shared faces for which we have a leaf point in the cone
1175:        The array holds candidate shared faces, each face is referred to by the leaf point */
1176:   PetscCall(PetscSectionCreate(comm, &candidateSection));
1177:   PetscCall(PetscSectionSetChart(candidateSection, 0, Nr));
1178:   {
1179:     PetscHMapIJ faceHash;

1181:     PetscCall(PetscHMapIJCreate(&faceHash));
1182:     for (l = 0; l < Nl; ++l) {
1183:       const PetscInt p = localPoints[l];

1185:       if (debug) PetscCall(PetscSynchronizedPrintf(comm, "[%d]  First pass leaf point %" PetscInt_FMT "\n", rank, p));
1186:       PetscCall(DMPlexAddSharedFace_Private(dm, candidateSection, NULL, faceHash, p, debug));
1187:     }
1188:     PetscCall(PetscHMapIJClear(faceHash));
1189:     PetscCall(PetscSectionSetUp(candidateSection));
1190:     PetscCall(PetscSectionGetStorageSize(candidateSection, &candidatesSize));
1191:     PetscCall(PetscMalloc1(candidatesSize, &candidates));
1192:     for (l = 0; l < Nl; ++l) {
1193:       const PetscInt p = localPoints[l];

1195:       if (debug) PetscCall(PetscSynchronizedPrintf(comm, "[%d]  Second pass leaf point %" PetscInt_FMT "\n", rank, p));
1196:       PetscCall(DMPlexAddSharedFace_Private(dm, candidateSection, candidates, faceHash, p, debug));
1197:     }
1198:     PetscCall(PetscHMapIJDestroy(&faceHash));
1199:     if (debug) PetscCall(PetscSynchronizedFlush(comm, NULL));
1200:   }
1201:   PetscCall(PetscObjectSetName((PetscObject)candidateSection, "Candidate Section"));
1202:   PetscCall(PetscObjectViewFromOptions((PetscObject)candidateSection, NULL, "-petscsection_interp_candidate_view"));
1203:   PetscCall(SFNodeArrayViewFromOptions(comm, "-petscsection_interp_candidate_view", "Candidates", NULL, candidatesSize, candidates));
1204:   /* Step 2: Gather candidate section / array pair into the root partition via inverse(multi(pointSF)). */
1205:   /*   Note that this section is indexed by offsets into leaves, not by point number */
1206:   {
1207:     PetscSF   sfMulti, sfInverse, sfCandidates;
1208:     PetscInt *remoteOffsets;

1210:     PetscCall(PetscSFGetMultiSF(pointSF, &sfMulti));
1211:     PetscCall(PetscSFCreateInverseSF(sfMulti, &sfInverse));
1212:     PetscCall(PetscSectionCreate(comm, &candidateRemoteSection));
1213:     PetscCall(PetscSFDistributeSection(sfInverse, candidateSection, &remoteOffsets, candidateRemoteSection));
1214:     PetscCall(PetscSFCreateSectionSF(sfInverse, candidateSection, remoteOffsets, candidateRemoteSection, &sfCandidates));
1215:     PetscCall(PetscSectionGetStorageSize(candidateRemoteSection, &candidatesRemoteSize));
1216:     PetscCall(PetscMalloc1(candidatesRemoteSize, &candidatesRemote));
1217:     PetscCall(PetscSFBcastBegin(sfCandidates, MPIU_2INT, candidates, candidatesRemote, MPI_REPLACE));
1218:     PetscCall(PetscSFBcastEnd(sfCandidates, MPIU_2INT, candidates, candidatesRemote, MPI_REPLACE));
1219:     PetscCall(PetscSFDestroy(&sfInverse));
1220:     PetscCall(PetscSFDestroy(&sfCandidates));
1221:     PetscCall(PetscFree(remoteOffsets));

1223:     PetscCall(PetscObjectSetName((PetscObject)candidateRemoteSection, "Remote Candidate Section"));
1224:     PetscCall(PetscObjectViewFromOptions((PetscObject)candidateRemoteSection, NULL, "-petscsection_interp_candidate_remote_view"));
1225:     PetscCall(SFNodeArrayViewFromOptions(comm, "-petscsection_interp_candidate_remote_view", "Remote Candidates", NULL, candidatesRemoteSize, candidatesRemote));
1226:   }
1227:   /* 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 */
1228:   {
1229:     PetscHMapIJKLRemote faceTable;
1230:     PetscInt            idx, idx2;

1232:     PetscCall(PetscHMapIJKLRemoteCreate(&faceTable));
1233:     /* There is a section point for every leaf attached to a given root point */
1234:     for (r = 0, idx = 0, idx2 = 0; r < Nr; ++r) {
1235:       PetscInt deg;

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

1240:         PetscCall(PetscSectionGetDof(candidateRemoteSection, idx, &dof));
1241:         PetscCall(PetscSectionGetOffset(candidateRemoteSection, idx, &offset));
1242:         /* dof may include many faces from the remote process */
1243:         for (d = 0; d < dof; ++d) {
1244:           const PetscInt         hidx  = offset + d;
1245:           const PetscInt         Np    = candidatesRemote[hidx].index + 1;
1246:           const PetscSFNode      rface = candidatesRemote[hidx + 1];
1247:           const PetscSFNode     *fcone = &candidatesRemote[hidx + 2];
1248:           PetscSFNode            fcp0;
1249:           const PetscSFNode      pmax = {PETSC_MAX_INT, PETSC_MAX_INT};
1250:           const PetscInt        *join = NULL;
1251:           PetscHMapIJKLRemoteKey key;
1252:           PetscHashIter          iter;
1253:           PetscBool              missing, mapToLocalPointFailed = PETSC_FALSE;
1254:           PetscInt               points[1024], p, joinSize;

1256:           if (debug)
1257:             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,
1258:                                               rface.index, r, idx, d, Np));
1259:           PetscCheck(Np <= 4, PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot handle face (%" PetscInt_FMT ", %" PetscInt_FMT ") at (%" PetscInt_FMT ", %" PetscInt_FMT ", %" PetscInt_FMT ") with %" PetscInt_FMT " cone points", rface.rank, rface.index, r, idx, d, Np);
1260:           fcp0.rank  = rank;
1261:           fcp0.index = r;
1262:           d += Np;
1263:           /* Put remote face in hash table */
1264:           key.i = fcp0;
1265:           key.j = fcone[0];
1266:           key.k = Np > 2 ? fcone[1] : pmax;
1267:           key.l = Np > 3 ? fcone[2] : pmax;
1268:           PetscCall(PetscSortSFNode(Np, (PetscSFNode *)&key));
1269:           PetscCall(PetscHMapIJKLRemotePut(faceTable, key, &iter, &missing));
1270:           if (missing) {
1271:             if (debug) PetscCall(PetscSynchronizedPrintf(PetscObjectComm((PetscObject)dm), "[%d]  Setting remote face (%" PetscInt_FMT ", %" PetscInt_FMT ")\n", rank, rface.index, rface.rank));
1272:             PetscCall(PetscHMapIJKLRemoteIterSet(faceTable, iter, rface));
1273:           } else {
1274:             PetscSFNode oface;

1276:             PetscCall(PetscHMapIJKLRemoteIterGet(faceTable, iter, &oface));
1277:             if ((rface.rank < oface.rank) || (rface.rank == oface.rank && rface.index < oface.index)) {
1278:               if (debug) PetscCall(PetscSynchronizedPrintf(PetscObjectComm((PetscObject)dm), "[%d]  Replacing with remote face (%" PetscInt_FMT ", %" PetscInt_FMT ")\n", rank, rface.index, rface.rank));
1279:               PetscCall(PetscHMapIJKLRemoteIterSet(faceTable, iter, rface));
1280:             }
1281:           }
1282:           /* Check for local face */
1283:           points[0] = r;
1284:           for (p = 1; p < Np; ++p) {
1285:             PetscCall(DMPlexMapToLocalPoint(dm, remoteHash, fcone[p - 1], &points[p], &mapToLocalPointFailed));
1286:             if (mapToLocalPointFailed) break; /* We got a point not in our overlap */
1287:             if (debug) PetscCall(PetscSynchronizedPrintf(PetscObjectComm((PetscObject)dm), "[%d]  Checking local candidate %" PetscInt_FMT "\n", rank, points[p]));
1288:           }
1289:           if (mapToLocalPointFailed) continue;
1290:           PetscCall(DMPlexGetJoin(dm, Np, points, &joinSize, &join));
1291:           if (joinSize == 1) {
1292:             PetscSFNode lface;
1293:             PetscSFNode oface;

1295:             /* Always replace with local face */
1296:             lface.rank  = rank;
1297:             lface.index = join[0];
1298:             PetscCall(PetscHMapIJKLRemoteIterGet(faceTable, iter, &oface));
1299:             if (debug)
1300:               PetscCall(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));
1301:             PetscCall(PetscHMapIJKLRemoteIterSet(faceTable, iter, lface));
1302:           }
1303:           PetscCall(DMPlexRestoreJoin(dm, Np, points, &joinSize, &join));
1304:         }
1305:       }
1306:       /* Put back faces for this root */
1307:       for (deg = 0; deg < rootdegree[r]; ++deg, ++idx2) {
1308:         PetscInt offset, dof, d;

1310:         PetscCall(PetscSectionGetDof(candidateRemoteSection, idx2, &dof));
1311:         PetscCall(PetscSectionGetOffset(candidateRemoteSection, idx2, &offset));
1312:         /* dof may include many faces from the remote process */
1313:         for (d = 0; d < dof; ++d) {
1314:           const PetscInt         hidx  = offset + d;
1315:           const PetscInt         Np    = candidatesRemote[hidx].index + 1;
1316:           const PetscSFNode     *fcone = &candidatesRemote[hidx + 2];
1317:           PetscSFNode            fcp0;
1318:           const PetscSFNode      pmax = {PETSC_MAX_INT, PETSC_MAX_INT};
1319:           PetscHMapIJKLRemoteKey key;
1320:           PetscHashIter          iter;
1321:           PetscBool              missing;

1323:           if (debug) PetscCall(PetscSynchronizedPrintf(PetscObjectComm((PetscObject)dm), "[%d]  Entering face at (%" PetscInt_FMT ", %" PetscInt_FMT ")\n", rank, r, idx));
1324:           PetscCheck(Np <= 4, PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot handle faces with %" PetscInt_FMT " cone points", Np);
1325:           fcp0.rank  = rank;
1326:           fcp0.index = r;
1327:           d += Np;
1328:           /* Find remote face in hash table */
1329:           key.i = fcp0;
1330:           key.j = fcone[0];
1331:           key.k = Np > 2 ? fcone[1] : pmax;
1332:           key.l = Np > 3 ? fcone[2] : pmax;
1333:           PetscCall(PetscSortSFNode(Np, (PetscSFNode *)&key));
1334:           if (debug)
1335:             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,
1336:                                               key.i.rank, key.i.index, key.j.rank, key.j.index, key.k.rank, key.k.index, key.l.rank, key.l.index));
1337:           PetscCall(PetscHMapIJKLRemotePut(faceTable, key, &iter, &missing));
1338:           PetscCheck(!missing, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Root %" PetscInt_FMT " Idx %" PetscInt_FMT " ought to have an associated face", r, idx2);
1339:           PetscCall(PetscHMapIJKLRemoteIterGet(faceTable, iter, &candidatesRemote[hidx]));
1340:         }
1341:       }
1342:     }
1343:     if (debug) PetscCall(PetscSynchronizedFlush(PetscObjectComm((PetscObject)dm), NULL));
1344:     PetscCall(PetscHMapIJKLRemoteDestroy(&faceTable));
1345:   }
1346:   /* Step 4: Push back owned faces */
1347:   {
1348:     PetscSF      sfMulti, sfClaims, sfPointNew;
1349:     PetscSFNode *remotePointsNew;
1350:     PetscInt    *remoteOffsets, *localPointsNew;
1351:     PetscInt     pStart, pEnd, r, NlNew, p;

1353:     /* 4) Push claims back to receiver via the MultiSF and derive new pointSF mapping on receiver */
1354:     PetscCall(PetscSFGetMultiSF(pointSF, &sfMulti));
1355:     PetscCall(PetscSectionCreate(comm, &claimSection));
1356:     PetscCall(PetscSFDistributeSection(sfMulti, candidateRemoteSection, &remoteOffsets, claimSection));
1357:     PetscCall(PetscSFCreateSectionSF(sfMulti, candidateRemoteSection, remoteOffsets, claimSection, &sfClaims));
1358:     PetscCall(PetscSectionGetStorageSize(claimSection, &claimsSize));
1359:     PetscCall(PetscMalloc1(claimsSize, &claims));
1360:     for (p = 0; p < claimsSize; ++p) claims[p].rank = -1;
1361:     PetscCall(PetscSFBcastBegin(sfClaims, MPIU_2INT, candidatesRemote, claims, MPI_REPLACE));
1362:     PetscCall(PetscSFBcastEnd(sfClaims, MPIU_2INT, candidatesRemote, claims, MPI_REPLACE));
1363:     PetscCall(PetscSFDestroy(&sfClaims));
1364:     PetscCall(PetscFree(remoteOffsets));
1365:     PetscCall(PetscObjectSetName((PetscObject)claimSection, "Claim Section"));
1366:     PetscCall(PetscObjectViewFromOptions((PetscObject)claimSection, NULL, "-petscsection_interp_claim_view"));
1367:     PetscCall(SFNodeArrayViewFromOptions(comm, "-petscsection_interp_claim_view", "Claims", NULL, claimsSize, claims));
1368:     /* Step 5) Walk the original section of local supports and add an SF entry for each updated item */
1369:     /* TODO I should not have to do a join here since I already put the face and its cone in the candidate section */
1370:     PetscCall(PetscHMapICreate(&claimshash));
1371:     for (r = 0; r < Nr; ++r) {
1372:       PetscInt dof, off, d;

1374:       if (debug) PetscCall(PetscSynchronizedPrintf(comm, "[%d]  Checking root for claims %" PetscInt_FMT "\n", rank, r));
1375:       PetscCall(PetscSectionGetDof(candidateSection, r, &dof));
1376:       PetscCall(PetscSectionGetOffset(candidateSection, r, &off));
1377:       for (d = 0; d < dof;) {
1378:         if (claims[off + d].rank >= 0) {
1379:           const PetscInt  faceInd = off + d;
1380:           const PetscInt  Np      = candidates[off + d].index;
1381:           const PetscInt *join    = NULL;
1382:           PetscInt        joinSize, points[1024], c;

1384:           if (debug) PetscCall(PetscSynchronizedPrintf(comm, "[%d]    Found claim for remote point (%" PetscInt_FMT ", %" PetscInt_FMT ")\n", rank, claims[faceInd].rank, claims[faceInd].index));
1385:           points[0] = r;
1386:           if (debug) PetscCall(PetscSynchronizedPrintf(comm, "[%d]      point %" PetscInt_FMT "\n", rank, points[0]));
1387:           for (c = 0, d += 2; c < Np; ++c, ++d) {
1388:             PetscCall(DMPlexMapToLocalPoint(dm, remoteHash, candidates[off + d], &points[c + 1], NULL));
1389:             if (debug) PetscCall(PetscSynchronizedPrintf(comm, "[%d]      point %" PetscInt_FMT "\n", rank, points[c + 1]));
1390:           }
1391:           PetscCall(DMPlexGetJoin(dm, Np + 1, points, &joinSize, &join));
1392:           if (joinSize == 1) {
1393:             if (claims[faceInd].rank == rank) {
1394:               if (debug) PetscCall(PetscSynchronizedPrintf(comm, "[%d]    Ignoring local face %" PetscInt_FMT " for non-remote partner\n", rank, join[0]));
1395:             } else {
1396:               if (debug) PetscCall(PetscSynchronizedPrintf(comm, "[%d]    Found local face %" PetscInt_FMT "\n", rank, join[0]));
1397:               PetscCall(PetscHMapISet(claimshash, join[0], faceInd));
1398:             }
1399:           } else {
1400:             if (debug) PetscCall(PetscSynchronizedPrintf(comm, "[%d]    Failed to find face\n", rank));
1401:           }
1402:           PetscCall(DMPlexRestoreJoin(dm, Np + 1, points, &joinSize, &join));
1403:         } else {
1404:           if (debug) PetscCall(PetscSynchronizedPrintf(comm, "[%d]    No claim for point %" PetscInt_FMT "\n", rank, r));
1405:           d += claims[off + d].index + 1;
1406:         }
1407:       }
1408:     }
1409:     if (debug) PetscCall(PetscSynchronizedFlush(comm, NULL));
1410:     /* Step 6) Create new pointSF from hashed claims */
1411:     PetscCall(PetscHMapIGetSize(claimshash, &NlNew));
1412:     PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
1413:     PetscCall(PetscMalloc1(Nl + NlNew, &localPointsNew));
1414:     PetscCall(PetscMalloc1(Nl + NlNew, &remotePointsNew));
1415:     for (l = 0; l < Nl; ++l) {
1416:       localPointsNew[l]        = localPoints[l];
1417:       remotePointsNew[l].index = remotePoints[l].index;
1418:       remotePointsNew[l].rank  = remotePoints[l].rank;
1419:     }
1420:     p = Nl;
1421:     PetscCall(PetscHMapIGetKeys(claimshash, &p, localPointsNew));
1422:     /* We sort new points, and assume they are numbered after all existing points */
1423:     PetscCall(PetscSortInt(NlNew, &localPointsNew[Nl]));
1424:     for (p = Nl; p < Nl + NlNew; ++p) {
1425:       PetscInt off;
1426:       PetscCall(PetscHMapIGet(claimshash, localPointsNew[p], &off));
1427:       PetscCheck(claims[off].rank >= 0 && claims[off].index >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid claim for local point %" PetscInt_FMT ", (%" PetscInt_FMT ", %" PetscInt_FMT ")", localPointsNew[p], claims[off].rank, claims[off].index);
1428:       remotePointsNew[p] = claims[off];
1429:     }
1430:     PetscCall(PetscSFCreate(comm, &sfPointNew));
1431:     PetscCall(PetscSFSetGraph(sfPointNew, pEnd - pStart, Nl + NlNew, localPointsNew, PETSC_OWN_POINTER, remotePointsNew, PETSC_OWN_POINTER));
1432:     PetscCall(PetscSFSetUp(sfPointNew));
1433:     PetscCall(DMSetPointSF(dm, sfPointNew));
1434:     PetscCall(PetscObjectViewFromOptions((PetscObject)sfPointNew, NULL, "-petscsf_interp_view"));
1435:     if (PetscDefined(USE_DEBUG)) PetscCall(DMPlexCheckPointSF(dm, sfPointNew, PETSC_FALSE));
1436:     PetscCall(PetscSFDestroy(&sfPointNew));
1437:     PetscCall(PetscHMapIDestroy(&claimshash));
1438:   }
1439:   PetscCall(PetscHMapIJDestroy(&remoteHash));
1440:   PetscCall(PetscSectionDestroy(&candidateSection));
1441:   PetscCall(PetscSectionDestroy(&candidateRemoteSection));
1442:   PetscCall(PetscSectionDestroy(&claimSection));
1443:   PetscCall(PetscFree(candidates));
1444:   PetscCall(PetscFree(candidatesRemote));
1445:   PetscCall(PetscFree(claims));
1446:   PetscCall(PetscLogEventEnd(DMPLEX_InterpolateSF, dm, 0, 0, 0));
1447:   PetscFunctionReturn(PETSC_SUCCESS);
1448: }

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

1453:   Collective

1455:   Input Parameter:
1456: . dm - The `DMPLEX` object with only cells and vertices

1458:   Output Parameter:
1459: . dmInt - The complete `DMPLEX` object

1461:   Level: intermediate

1463:   Note:
1464:   Labels and coordinates are copied.

1466:   Developer Notes:
1467:   It sets plex->interpolated = `DMPLEX_INTERPOLATED_FULL`.

1469: .seealso: `DMPLEX`, `DMPlexUninterpolate()`, `DMPlexCreateFromCellListPetsc()`, `DMPlexCopyCoordinates()`
1470: @*/
1471: PetscErrorCode DMPlexInterpolate(DM dm, DM *dmInt)
1472: {
1473:   DMPlexInterpolatedFlag interpolated;
1474:   DM                     idm, odm = dm;
1475:   PetscSF                sfPoint;
1476:   PetscInt               depth, dim, d;
1477:   const char            *name;
1478:   PetscBool              flg = PETSC_TRUE;

1480:   PetscFunctionBegin;
1482:   PetscAssertPointer(dmInt, 2);
1483:   PetscCall(PetscLogEventBegin(DMPLEX_Interpolate, dm, 0, 0, 0));
1484:   PetscCall(DMPlexGetDepth(dm, &depth));
1485:   PetscCall(DMGetDimension(dm, &dim));
1486:   PetscCall(DMPlexIsInterpolated(dm, &interpolated));
1487:   PetscCheck(interpolated != DMPLEX_INTERPOLATED_PARTIAL, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Not for partially interpolated meshes");
1488:   if (interpolated == DMPLEX_INTERPOLATED_FULL) {
1489:     PetscCall(PetscObjectReference((PetscObject)dm));
1490:     idm = dm;
1491:   } else {
1492:     for (d = 1; d < dim; ++d) {
1493:       /* Create interpolated mesh */
1494:       PetscCall(DMCreate(PetscObjectComm((PetscObject)dm), &idm));
1495:       PetscCall(DMSetType(idm, DMPLEX));
1496:       PetscCall(DMSetDimension(idm, dim));
1497:       if (depth > 0) {
1498:         PetscCall(DMPlexInterpolateFaces_Internal(odm, 1, idm));
1499:         PetscCall(DMGetPointSF(odm, &sfPoint));
1500:         if (PetscDefined(USE_DEBUG)) PetscCall(DMPlexCheckPointSF(odm, sfPoint, PETSC_FALSE));
1501:         {
1502:           /* TODO: We need to systematically fix cases of distributed Plexes with no graph set */
1503:           PetscInt nroots;
1504:           PetscCall(PetscSFGetGraph(sfPoint, &nroots, NULL, NULL, NULL));
1505:           if (nroots >= 0) PetscCall(DMPlexInterpolatePointSF(idm, sfPoint));
1506:         }
1507:       }
1508:       if (odm != dm) PetscCall(DMDestroy(&odm));
1509:       odm = idm;
1510:     }
1511:     PetscCall(PetscObjectGetName((PetscObject)dm, &name));
1512:     PetscCall(PetscObjectSetName((PetscObject)idm, name));
1513:     PetscCall(DMPlexCopyCoordinates(dm, idm));
1514:     PetscCall(DMCopyLabels(dm, idm, PETSC_COPY_VALUES, PETSC_FALSE, DM_COPY_LABELS_FAIL));
1515:     PetscCall(PetscOptionsGetBool(((PetscObject)dm)->options, ((PetscObject)dm)->prefix, "-dm_plex_interpolate_orient_interfaces", &flg, NULL));
1516:     if (flg) PetscCall(DMPlexOrientInterface_Internal(idm));
1517:   }
1518:   /* This function makes the mesh fully interpolated on all ranks */
1519:   {
1520:     DM_Plex *plex      = (DM_Plex *)idm->data;
1521:     plex->interpolated = plex->interpolatedCollective = DMPLEX_INTERPOLATED_FULL;
1522:   }
1523:   PetscCall(DMPlexCopy_Internal(dm, PETSC_TRUE, PETSC_TRUE, idm));
1524:   *dmInt = idm;
1525:   PetscCall(PetscLogEventEnd(DMPLEX_Interpolate, dm, 0, 0, 0));
1526:   PetscFunctionReturn(PETSC_SUCCESS);
1527: }

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

1532:   Collective

1534:   Input Parameter:
1535: . dmA - The `DMPLEX` object with initial coordinates

1537:   Output Parameter:
1538: . dmB - The `DMPLEX` object with copied coordinates

1540:   Level: intermediate

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

1545: .seealso: `DMPLEX`, `DMCopyLabels()`, `DMGetCoordinates()`, `DMGetCoordinatesLocal()`, `DMGetCoordinateDM()`, `DMGetCoordinateSection()`
1546: @*/
1547: PetscErrorCode DMPlexCopyCoordinates(DM dmA, DM dmB)
1548: {
1549:   Vec          coordinatesA, coordinatesB;
1550:   VecType      vtype;
1551:   PetscSection coordSectionA, coordSectionB;
1552:   PetscScalar *coordsA, *coordsB;
1553:   PetscInt     spaceDim, Nf, vStartA, vStartB, vEndA, vEndB, coordSizeB, v, d;
1554:   PetscInt     cStartA, cEndA, cStartB, cEndB, cS, cE, cdim;
1555:   PetscBool    lc = PETSC_FALSE;

1557:   PetscFunctionBegin;
1560:   if (dmA == dmB) PetscFunctionReturn(PETSC_SUCCESS);
1561:   PetscCall(DMGetCoordinateDim(dmA, &cdim));
1562:   PetscCall(DMSetCoordinateDim(dmB, cdim));
1563:   PetscCall(DMPlexGetDepthStratum(dmA, 0, &vStartA, &vEndA));
1564:   PetscCall(DMPlexGetDepthStratum(dmB, 0, &vStartB, &vEndB));
1565:   PetscCheck((vEndA - vStartA) == (vEndB - vStartB), PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "The number of vertices in first DM %" PetscInt_FMT " != %" PetscInt_FMT " in the second DM", vEndA - vStartA, vEndB - vStartB);
1566:   /* Copy over discretization if it exists */
1567:   {
1568:     DM                 cdmA, cdmB;
1569:     PetscDS            dsA, dsB;
1570:     PetscObject        objA, objB;
1571:     PetscClassId       idA, idB;
1572:     const PetscScalar *constants;
1573:     PetscInt           cdim, Nc;

1575:     PetscCall(DMGetCoordinateDM(dmA, &cdmA));
1576:     PetscCall(DMGetCoordinateDM(dmB, &cdmB));
1577:     PetscCall(DMGetField(cdmA, 0, NULL, &objA));
1578:     PetscCall(DMGetField(cdmB, 0, NULL, &objB));
1579:     PetscCall(PetscObjectGetClassId(objA, &idA));
1580:     PetscCall(PetscObjectGetClassId(objB, &idB));
1581:     if ((idA == PETSCFE_CLASSID) && (idA != idB)) {
1582:       PetscCall(DMSetField(cdmB, 0, NULL, objA));
1583:       PetscCall(DMCreateDS(cdmB));
1584:       PetscCall(DMGetDS(cdmA, &dsA));
1585:       PetscCall(DMGetDS(cdmB, &dsB));
1586:       PetscCall(PetscDSGetCoordinateDimension(dsA, &cdim));
1587:       PetscCall(PetscDSSetCoordinateDimension(dsB, cdim));
1588:       PetscCall(PetscDSGetConstants(dsA, &Nc, &constants));
1589:       PetscCall(PetscDSSetConstants(dsB, Nc, (PetscScalar *)constants));
1590:     }
1591:   }
1592:   PetscCall(DMPlexGetHeightStratum(dmA, 0, &cStartA, &cEndA));
1593:   PetscCall(DMPlexGetHeightStratum(dmB, 0, &cStartB, &cEndB));
1594:   PetscCall(DMGetCoordinateSection(dmA, &coordSectionA));
1595:   PetscCall(DMGetCoordinateSection(dmB, &coordSectionB));
1596:   if (coordSectionA == coordSectionB) PetscFunctionReturn(PETSC_SUCCESS);
1597:   PetscCall(PetscSectionGetNumFields(coordSectionA, &Nf));
1598:   if (!Nf) PetscFunctionReturn(PETSC_SUCCESS);
1599:   PetscCheck(Nf <= 1, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "The number of coordinate fields must be 1, not %" PetscInt_FMT, Nf);
1600:   if (!coordSectionB) {
1601:     PetscInt dim;

1603:     PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)coordSectionA), &coordSectionB));
1604:     PetscCall(DMGetCoordinateDim(dmA, &dim));
1605:     PetscCall(DMSetCoordinateSection(dmB, dim, coordSectionB));
1606:     PetscCall(PetscObjectDereference((PetscObject)coordSectionB));
1607:   }
1608:   PetscCall(PetscSectionSetNumFields(coordSectionB, 1));
1609:   PetscCall(PetscSectionGetFieldComponents(coordSectionA, 0, &spaceDim));
1610:   PetscCall(PetscSectionSetFieldComponents(coordSectionB, 0, spaceDim));
1611:   PetscCall(PetscSectionGetChart(coordSectionA, &cS, &cE));
1612:   if (cStartA <= cS && cS < cEndA) { /* localized coordinates */
1613:     PetscCheck((cEndA - cStartA) == (cEndB - cStartB), PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "The number of cells in first DM %" PetscInt_FMT " != %" PetscInt_FMT " in the second DM", cEndA - cStartA, cEndB - cStartB);
1614:     cS = cS - cStartA + cStartB;
1615:     cE = vEndB;
1616:     lc = PETSC_TRUE;
1617:   } else {
1618:     cS = vStartB;
1619:     cE = vEndB;
1620:   }
1621:   PetscCall(PetscSectionSetChart(coordSectionB, cS, cE));
1622:   for (v = vStartB; v < vEndB; ++v) {
1623:     PetscCall(PetscSectionSetDof(coordSectionB, v, spaceDim));
1624:     PetscCall(PetscSectionSetFieldDof(coordSectionB, v, 0, spaceDim));
1625:   }
1626:   if (lc) { /* localized coordinates */
1627:     PetscInt c;

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

1632:       PetscCall(PetscSectionGetDof(coordSectionA, c + cStartA, &dof));
1633:       PetscCall(PetscSectionSetDof(coordSectionB, c + cStartB, dof));
1634:       PetscCall(PetscSectionSetFieldDof(coordSectionB, c + cStartB, 0, dof));
1635:     }
1636:   }
1637:   PetscCall(PetscSectionSetUp(coordSectionB));
1638:   PetscCall(PetscSectionGetStorageSize(coordSectionB, &coordSizeB));
1639:   PetscCall(DMGetCoordinatesLocal(dmA, &coordinatesA));
1640:   PetscCall(VecCreate(PETSC_COMM_SELF, &coordinatesB));
1641:   PetscCall(PetscObjectSetName((PetscObject)coordinatesB, "coordinates"));
1642:   PetscCall(VecSetSizes(coordinatesB, coordSizeB, PETSC_DETERMINE));
1643:   PetscCall(VecGetBlockSize(coordinatesA, &d));
1644:   PetscCall(VecSetBlockSize(coordinatesB, d));
1645:   PetscCall(VecGetType(coordinatesA, &vtype));
1646:   PetscCall(VecSetType(coordinatesB, vtype));
1647:   PetscCall(VecGetArray(coordinatesA, &coordsA));
1648:   PetscCall(VecGetArray(coordinatesB, &coordsB));
1649:   for (v = 0; v < vEndB - vStartB; ++v) {
1650:     PetscInt offA, offB;

1652:     PetscCall(PetscSectionGetOffset(coordSectionA, v + vStartA, &offA));
1653:     PetscCall(PetscSectionGetOffset(coordSectionB, v + vStartB, &offB));
1654:     for (d = 0; d < spaceDim; ++d) coordsB[offB + d] = coordsA[offA + d];
1655:   }
1656:   if (lc) { /* localized coordinates */
1657:     PetscInt c;

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

1662:       PetscCall(PetscSectionGetOffset(coordSectionA, c + cStartA, &offA));
1663:       PetscCall(PetscSectionGetOffset(coordSectionB, c + cStartB, &offB));
1664:       PetscCall(PetscSectionGetDof(coordSectionA, c + cStartA, &dof));
1665:       PetscCall(PetscArraycpy(coordsB + offB, coordsA + offA, dof));
1666:     }
1667:   }
1668:   PetscCall(VecRestoreArray(coordinatesA, &coordsA));
1669:   PetscCall(VecRestoreArray(coordinatesB, &coordsB));
1670:   PetscCall(DMSetCoordinatesLocal(dmB, coordinatesB));
1671:   PetscCall(VecDestroy(&coordinatesB));
1672:   PetscFunctionReturn(PETSC_SUCCESS);
1673: }

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

1678:   Collective

1680:   Input Parameter:
1681: . dm - The complete `DMPLEX` object

1683:   Output Parameter:
1684: . dmUnint - The `DMPLEX` object with only cells and vertices

1686:   Level: intermediate

1688:   Note:
1689:   It does not copy over the coordinates.

1691:   Developer Notes:
1692:   Sets plex->interpolated = `DMPLEX_INTERPOLATED_NONE`.

1694: .seealso: `DMPLEX`, `DMPlexInterpolate()`, `DMPlexCreateFromCellListPetsc()`, `DMPlexCopyCoordinates()`
1695: @*/
1696: PetscErrorCode DMPlexUninterpolate(DM dm, DM *dmUnint)
1697: {
1698:   DMPlexInterpolatedFlag interpolated;
1699:   DM                     udm;
1700:   PetscInt               dim, vStart, vEnd, cStart, cEnd, c, maxConeSize = 0, *cone;

1702:   PetscFunctionBegin;
1704:   PetscAssertPointer(dmUnint, 2);
1705:   PetscCall(PetscLogEventBegin(DMPLEX_Uninterpolate, dm, 0, 0, 0));
1706:   PetscCall(DMGetDimension(dm, &dim));
1707:   PetscCall(DMPlexIsInterpolated(dm, &interpolated));
1708:   PetscCheck(interpolated != DMPLEX_INTERPOLATED_PARTIAL, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Not for partially interpolated meshes");
1709:   if (interpolated == DMPLEX_INTERPOLATED_NONE || dim <= 1) {
1710:     /* in case dim <= 1 just keep the DMPLEX_INTERPOLATED_FULL flag */
1711:     PetscCall(PetscObjectReference((PetscObject)dm));
1712:     *dmUnint = dm;
1713:     PetscFunctionReturn(PETSC_SUCCESS);
1714:   }
1715:   PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd));
1716:   PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
1717:   PetscCall(DMCreate(PetscObjectComm((PetscObject)dm), &udm));
1718:   PetscCall(DMSetType(udm, DMPLEX));
1719:   PetscCall(DMSetDimension(udm, dim));
1720:   PetscCall(DMPlexSetChart(udm, cStart, vEnd));
1721:   for (c = cStart; c < cEnd; ++c) {
1722:     PetscInt *closure = NULL, closureSize, cl, coneSize = 0;

1724:     PetscCall(DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure));
1725:     for (cl = 0; cl < closureSize * 2; cl += 2) {
1726:       const PetscInt p = closure[cl];

1728:       if ((p >= vStart) && (p < vEnd)) ++coneSize;
1729:     }
1730:     PetscCall(DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure));
1731:     PetscCall(DMPlexSetConeSize(udm, c, coneSize));
1732:     maxConeSize = PetscMax(maxConeSize, coneSize);
1733:   }
1734:   PetscCall(DMSetUp(udm));
1735:   PetscCall(PetscMalloc1(maxConeSize, &cone));
1736:   for (c = cStart; c < cEnd; ++c) {
1737:     PetscInt *closure = NULL, closureSize, cl, coneSize = 0;

1739:     PetscCall(DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure));
1740:     for (cl = 0; cl < closureSize * 2; cl += 2) {
1741:       const PetscInt p = closure[cl];

1743:       if ((p >= vStart) && (p < vEnd)) cone[coneSize++] = p;
1744:     }
1745:     PetscCall(DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure));
1746:     PetscCall(DMPlexSetCone(udm, c, cone));
1747:   }
1748:   PetscCall(PetscFree(cone));
1749:   PetscCall(DMPlexSymmetrize(udm));
1750:   PetscCall(DMPlexStratify(udm));
1751:   /* Reduce SF */
1752:   {
1753:     PetscSF            sfPoint, sfPointUn;
1754:     const PetscSFNode *remotePoints;
1755:     const PetscInt    *localPoints;
1756:     PetscSFNode       *remotePointsUn;
1757:     PetscInt          *localPointsUn;
1758:     PetscInt           numRoots, numLeaves, l;
1759:     PetscInt           numLeavesUn = 0, n = 0;

1761:     /* Get original SF information */
1762:     PetscCall(DMGetPointSF(dm, &sfPoint));
1763:     if (PetscDefined(USE_DEBUG)) PetscCall(DMPlexCheckPointSF(dm, sfPoint, PETSC_FALSE));
1764:     PetscCall(DMGetPointSF(udm, &sfPointUn));
1765:     PetscCall(PetscSFGetGraph(sfPoint, &numRoots, &numLeaves, &localPoints, &remotePoints));
1766:     if (numRoots >= 0) {
1767:       /* Allocate space for cells and vertices */
1768:       for (l = 0; l < numLeaves; ++l) {
1769:         const PetscInt p = localPoints[l];

1771:         if ((vStart <= p && p < vEnd) || (cStart <= p && p < cEnd)) numLeavesUn++;
1772:       }
1773:       /* Fill in leaves */
1774:       PetscCall(PetscMalloc1(numLeavesUn, &remotePointsUn));
1775:       PetscCall(PetscMalloc1(numLeavesUn, &localPointsUn));
1776:       for (l = 0; l < numLeaves; l++) {
1777:         const PetscInt p = localPoints[l];

1779:         if ((vStart <= p && p < vEnd) || (cStart <= p && p < cEnd)) {
1780:           localPointsUn[n]        = p;
1781:           remotePointsUn[n].rank  = remotePoints[l].rank;
1782:           remotePointsUn[n].index = remotePoints[l].index;
1783:           ++n;
1784:         }
1785:       }
1786:       PetscCheck(n == numLeavesUn, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Inconsistent number of leaves %" PetscInt_FMT " != %" PetscInt_FMT, n, numLeavesUn);
1787:       PetscCall(PetscSFSetGraph(sfPointUn, cEnd - cStart + vEnd - vStart, numLeavesUn, localPointsUn, PETSC_OWN_POINTER, remotePointsUn, PETSC_OWN_POINTER));
1788:     }
1789:   }
1790:   /* This function makes the mesh fully uninterpolated on all ranks */
1791:   {
1792:     DM_Plex *plex      = (DM_Plex *)udm->data;
1793:     plex->interpolated = plex->interpolatedCollective = DMPLEX_INTERPOLATED_NONE;
1794:   }
1795:   PetscCall(DMPlexCopy_Internal(dm, PETSC_TRUE, PETSC_TRUE, udm));
1796:   if (PetscDefined(USE_DEBUG)) PetscCall(DMPlexCheckPointSF(udm, NULL, PETSC_FALSE));
1797:   *dmUnint = udm;
1798:   PetscCall(PetscLogEventEnd(DMPLEX_Uninterpolate, dm, 0, 0, 0));
1799:   PetscFunctionReturn(PETSC_SUCCESS);
1800: }

1802: static PetscErrorCode DMPlexIsInterpolated_Internal(DM dm, DMPlexInterpolatedFlag *interpolated)
1803: {
1804:   PetscInt coneSize, depth, dim, h, p, pStart, pEnd;
1805:   MPI_Comm comm;

1807:   PetscFunctionBegin;
1808:   PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
1809:   PetscCall(DMPlexGetDepth(dm, &depth));
1810:   PetscCall(DMGetDimension(dm, &dim));

1812:   if (depth == dim) {
1813:     *interpolated = DMPLEX_INTERPOLATED_FULL;
1814:     if (!dim) goto finish;

1816:     /* Check points at height = dim are vertices (have no cones) */
1817:     PetscCall(DMPlexGetHeightStratum(dm, dim, &pStart, &pEnd));
1818:     for (p = pStart; p < pEnd; p++) {
1819:       PetscCall(DMPlexGetConeSize(dm, p, &coneSize));
1820:       if (coneSize) {
1821:         *interpolated = DMPLEX_INTERPOLATED_PARTIAL;
1822:         goto finish;
1823:       }
1824:     }

1826:     /* Check points at height < dim have cones */
1827:     for (h = 0; h < dim; h++) {
1828:       PetscCall(DMPlexGetHeightStratum(dm, h, &pStart, &pEnd));
1829:       for (p = pStart; p < pEnd; p++) {
1830:         PetscCall(DMPlexGetConeSize(dm, p, &coneSize));
1831:         if (!coneSize) {
1832:           *interpolated = DMPLEX_INTERPOLATED_PARTIAL;
1833:           goto finish;
1834:         }
1835:       }
1836:     }
1837:   } else if (depth == 1) {
1838:     *interpolated = DMPLEX_INTERPOLATED_NONE;
1839:   } else {
1840:     *interpolated = DMPLEX_INTERPOLATED_PARTIAL;
1841:   }
1842: finish:
1843:   PetscFunctionReturn(PETSC_SUCCESS);
1844: }

1846: /*@
1847:   DMPlexIsInterpolated - Find out to what extent the `DMPLEX` is topologically interpolated.

1849:   Not Collective

1851:   Input Parameter:
1852: . dm - The `DMPLEX` object

1854:   Output Parameter:
1855: . interpolated - Flag whether the `DM` is interpolated

1857:   Level: intermediate

1859:   Notes:
1860:   Unlike `DMPlexIsInterpolatedCollective()`, this is NOT collective
1861:   so the results can be different on different ranks in special cases.
1862:   However, `DMPlexInterpolate()` guarantees the result is the same on all.

1864:   Unlike `DMPlexIsInterpolatedCollective()`, this cannot return `DMPLEX_INTERPOLATED_MIXED`.

1866:   Developer Notes:
1867:   Initially, plex->interpolated = `DMPLEX_INTERPOLATED_INVALID`.

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

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

1875:   `DMPlexInterpolate()` sets plex->interpolated = `DMPLEX_INTERPOLATED_FULL`,
1876:   and DMPlexUninterpolate() sets plex->interpolated = `DMPLEX_INTERPOLATED_NONE`.

1878: .seealso: `DMPLEX`, `DMPlexInterpolate()`, `DMPlexIsInterpolatedCollective()`
1879: @*/
1880: PetscErrorCode DMPlexIsInterpolated(DM dm, DMPlexInterpolatedFlag *interpolated)
1881: {
1882:   DM_Plex *plex = (DM_Plex *)dm->data;

1884:   PetscFunctionBegin;
1886:   PetscAssertPointer(interpolated, 2);
1887:   if (plex->interpolated < 0) {
1888:     PetscCall(DMPlexIsInterpolated_Internal(dm, &plex->interpolated));
1889:   } else if (PetscDefined(USE_DEBUG)) {
1890:     DMPlexInterpolatedFlag flg;

1892:     PetscCall(DMPlexIsInterpolated_Internal(dm, &flg));
1893:     PetscCheck(plex->tr || flg == plex->interpolated, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Stashed DMPlexInterpolatedFlag %s is inconsistent with current %s", DMPlexInterpolatedFlags[plex->interpolated], DMPlexInterpolatedFlags[flg]);
1894:   }
1895:   *interpolated = plex->interpolated;
1896:   PetscFunctionReturn(PETSC_SUCCESS);
1897: }

1899: /*@
1900:   DMPlexIsInterpolatedCollective - Find out to what extent the `DMPLEX` is topologically interpolated (in collective manner).

1902:   Collective

1904:   Input Parameter:
1905: . dm - The `DMPLEX` object

1907:   Output Parameter:
1908: . interpolated - Flag whether the `DM` is interpolated

1910:   Level: intermediate

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

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

1917:   Developer Notes:
1918:   Initially, plex->interpolatedCollective = `DMPLEX_INTERPOLATED_INVALID`.

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

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

1927: .seealso: `DMPLEX`, `DMPlexInterpolate()`, `DMPlexIsInterpolated()`
1928: @*/
1929: PetscErrorCode DMPlexIsInterpolatedCollective(DM dm, DMPlexInterpolatedFlag *interpolated)
1930: {
1931:   DM_Plex  *plex  = (DM_Plex *)dm->data;
1932:   PetscBool debug = PETSC_FALSE;

1934:   PetscFunctionBegin;
1936:   PetscAssertPointer(interpolated, 2);
1937:   PetscCall(PetscOptionsGetBool(((PetscObject)dm)->options, ((PetscObject)dm)->prefix, "-dm_plex_is_interpolated_collective_debug", &debug, NULL));
1938:   if (plex->interpolatedCollective < 0) {
1939:     DMPlexInterpolatedFlag min, max;
1940:     MPI_Comm               comm;

1942:     PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
1943:     PetscCall(DMPlexIsInterpolated(dm, &plex->interpolatedCollective));
1944:     PetscCall(MPIU_Allreduce(&plex->interpolatedCollective, &min, 1, MPIU_ENUM, MPI_MIN, comm));
1945:     PetscCall(MPIU_Allreduce(&plex->interpolatedCollective, &max, 1, MPIU_ENUM, MPI_MAX, comm));
1946:     if (min != max) plex->interpolatedCollective = DMPLEX_INTERPOLATED_MIXED;
1947:     if (debug) {
1948:       PetscMPIInt rank;

1950:       PetscCallMPI(MPI_Comm_rank(comm, &rank));
1951:       PetscCall(PetscSynchronizedPrintf(comm, "[%d] interpolated=%s interpolatedCollective=%s\n", rank, DMPlexInterpolatedFlags[plex->interpolated], DMPlexInterpolatedFlags[plex->interpolatedCollective]));
1952:       PetscCall(PetscSynchronizedFlush(comm, PETSC_STDOUT));
1953:     }
1954:   }
1955:   *interpolated = plex->interpolatedCollective;
1956:   PetscFunctionReturn(PETSC_SUCCESS);
1957: }