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:   PetscFunctionBegin;
 40:   for (i = 1; i < n; ++i) {
 41:     PetscSFNode x = A[i];
 42:     PetscInt    j;

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

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

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

474: PetscErrorCode DMPlexRestoreRawFaces_Internal(DM dm, DMPolytopeType ct, const PetscInt cone[], PetscInt *numFaces, const DMPolytopeType *faceTypes[], const PetscInt *faceSizes[], const PetscInt *faces[])
475: {
476:   PetscFunctionBegin;
477:   if (faceTypes) PetscCall(DMRestoreWorkArray(dm, 0, MPIU_INT, (void *)faceTypes));
478:   if (faceSizes) PetscCall(DMRestoreWorkArray(dm, 0, MPIU_INT, (void *)faceSizes));
479:   if (faces) PetscCall(DMRestoreWorkArray(dm, 0, MPIU_INT, (void *)faces));
480:   PetscFunctionReturn(PETSC_SUCCESS);
481: }

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

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

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

516:       PetscCheck(faceSize <= 4, PETSC_COMM_SELF, PETSC_ERR_SUP, "Do not support faces of size %" PetscInt_FMT " > 4", faceSize);
517:       key.i = face[0];
518:       key.j = faceSize > 1 ? face[1] : PETSC_MAX_INT;
519:       key.k = faceSize > 2 ? face[2] : PETSC_MAX_INT;
520:       key.l = faceSize > 3 ? face[3] : PETSC_MAX_INT;
521:       PetscCall(PetscSortInt(faceSize, (PetscInt *)&key));
522:       PetscCall(PetscHashIJKLPut(faceTable, key, &iter, &missing));
523:       if (missing) {
524:         PetscCall(PetscHashIJKLIterSet(faceTable, iter, fEnd++));
525:         ++faceTypeNum[faceType];
526:       }
527:     }
528:     PetscCall(DMPlexRestoreRawFaces_Internal(dm, ct, cone, &numFaces, &faceTypes, &faceSizes, &faces));
529:   }
530:   /* We need to number faces contiguously among types */
531:   {
532:     PetscInt faceTypeStart[DM_NUM_POLYTOPES], ct, numFT = 0;

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

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

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

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

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

615:       PetscCheck(faceSize <= 4, PETSC_COMM_SELF, PETSC_ERR_SUP, "Do not support faces of size %" PetscInt_FMT " > 4", faceSize);
616:       key.i = face[0];
617:       key.j = faceSize > 1 ? face[1] : PETSC_MAX_INT;
618:       key.k = faceSize > 2 ? face[2] : PETSC_MAX_INT;
619:       key.l = faceSize > 3 ? face[3] : PETSC_MAX_INT;
620:       PetscCall(PetscSortInt(faceSize, (PetscInt *)&key));
621:       PetscCall(PetscHashIJKLPut(faceTable, key, &iter, &missing));
622:       PetscCheck(!missing, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Missing face (cell %" PetscInt_FMT ", lf %" PetscInt_FMT ")", c, cf);
623:       PetscCall(PetscHashIJKLIterGet(faceTable, iter, &f));
624:       PetscCall(DMPlexSetConeSize(idm, f, faceSize));
625:       PetscCall(DMPlexSetCellType(idm, f, faceType));
626:     }
627:     PetscCall(DMPlexRestoreRawFaces_Internal(dm, ct, cone, &numFaces, &faceTypes, &faceSizes, &faces));
628:   }
629:   PetscCall(DMSetUp(idm));
630:   /* Initialize cones so we do not need the bash table to tell us that a cone has been set */
631:   {
632:     PetscSection cs;
633:     PetscInt    *cones, csize;

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

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

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

673:       PetscCheck(faceSize <= 4, PETSC_COMM_SELF, PETSC_ERR_SUP, "Do not support faces of size %" PetscInt_FMT " > 4", faceSize);
674:       key.i = face[0];
675:       key.j = faceSize > 1 ? face[1] : PETSC_MAX_INT;
676:       key.k = faceSize > 2 ? face[2] : PETSC_MAX_INT;
677:       key.l = faceSize > 3 ? face[3] : PETSC_MAX_INT;
678:       PetscCall(PetscSortInt(faceSize, (PetscInt *)&key));
679:       PetscCall(PetscHashIJKLPut(faceTable, key, &iter, &missing));
680:       PetscCall(PetscHashIJKLIterGet(faceTable, iter, &f));
681:       PetscCall(DMPlexInsertCone(idm, c, cf, f));
682:       PetscCall(DMPlexGetCone(idm, f, &fcone));
683:       if (fcone[0] < 0) PetscCall(DMPlexSetCone(idm, f, face));
684:       {
685:         const PetscInt *cone;
686:         PetscInt        coneSize, ornt;

688:         PetscCall(DMPlexGetConeSize(idm, f, &coneSize));
689:         PetscCall(DMPlexGetCone(idm, f, &cone));
690:         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);
691:         /* Notice that we have to use vertices here because the lower dimensional faces have not been created yet */
692:         PetscCall(DMPolytopeGetVertexOrientation(faceType, cone, face, &ornt));
693:         PetscCall(DMPlexInsertConeOrientation(idm, c, cf, ornt));
694:       }
695:     }
696:     PetscCall(DMPlexRestoreRawFaces_Internal(dm, ct, cone, &numFaces, &faceTypes, &faceSizes, &faces));
697:   }
698:   PetscCall(PetscHashIJKLDestroy(&faceTable));
699:   PetscCall(DMPlexSymmetrize(idm));
700:   PetscCall(DMPlexStratify(idm));
701:   PetscFunctionReturn(PETSC_SUCCESS);
702: }

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

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

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

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

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

812:     if (leaves[p][0] < 0) continue; /* Ignore vertices */
813:     PetscCall(DMPlexGetCellType(dm, p, &ct));
814:     PetscCall(DMPlexGetConeSize(dm, p, &coneSize));
815:     PetscCall(DMPlexGetCone(dm, p, &cone));
816:     if (debug) {
817:       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]));
818:     }
819:     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]) {
820:       /* Translate these leaves to my cone points; mainCone means desired order p's cone points */
821:       for (c = 0; c < PetscMin(coneSize, 4); ++c) {
822:         PetscInt rS, rN;

824:         if (leavesRanks[p][c] == rank) {
825:           /* A local leaf is just taken as it is */
826:           mainCone[c] = leaves[p][c];
827:           continue;
828:         }
829:         /* Find index of rank leavesRanks[p][c] among remote ranks */
830:         /* No need for PetscMPIIntCast because these integers were originally cast from PetscMPIInt. */
831:         PetscCall(PetscFindMPIInt((PetscMPIInt)leavesRanks[p][c], nranks, ranks, &r));
832:         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]);
833:         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]);
834:         /* Find point leaves[p][c] among remote points aimed at rank leavesRanks[p][c] */
835:         rS = roffset[r];
836:         rN = roffset[r + 1] - rS;
837:         PetscCall(PetscFindInt(leaves[p][c], rN, &rremote1[rS], &ind0));
838:         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]);
839:         /* Get the corresponding local point */
840:         mainCone[c] = rmine1[rS + ind0];
841:       }
842:       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]));
843:       /* Set the desired order of p's cone points and fix orientations accordingly */
844:       PetscCall(DMPolytopeGetOrientation(ct, cone, mainCone, &o));
845:       PetscCall(DMPlexOrientPoint(dm, p, o));
846:     } else if (debug) PetscCall(PetscSynchronizedPrintf(comm, " ==\n"));
847:   }
848:   if (debug) {
849:     PetscCall(PetscSynchronizedFlush(comm, NULL));
850:     PetscCallMPI(MPI_Barrier(comm));
851:   }
852:   PetscCall(DMViewFromOptions(dm, NULL, "-after_orient_interface_dm_view"));
853:   PetscCall(PetscFree4(roots, leaves, rootsRanks, leavesRanks));
854:   PetscCall(PetscFree2(rmine1, rremote1));
855:   PetscFunctionReturn(PETSC_SUCCESS);
856: }

858: static PetscErrorCode IntArrayViewFromOptions(MPI_Comm comm, const char opt[], const char name[], const char idxname[], const char valname[], PetscInt n, const PetscInt a[])
859: {
860:   PetscInt    idx;
861:   PetscMPIInt rank;
862:   PetscBool   flg;

864:   PetscFunctionBegin;
865:   PetscCall(PetscOptionsHasName(NULL, NULL, opt, &flg));
866:   if (!flg) PetscFunctionReturn(PETSC_SUCCESS);
867:   PetscCallMPI(MPI_Comm_rank(comm, &rank));
868:   PetscCall(PetscSynchronizedPrintf(comm, "[%d]%s:\n", rank, name));
869:   for (idx = 0; idx < n; ++idx) PetscCall(PetscSynchronizedPrintf(comm, "[%d]%s %" PetscInt_FMT " %s %" PetscInt_FMT "\n", rank, idxname, idx, valname, a[idx]));
870:   PetscCall(PetscSynchronizedFlush(comm, NULL));
871:   PetscFunctionReturn(PETSC_SUCCESS);
872: }

874: static PetscErrorCode SFNodeArrayViewFromOptions(MPI_Comm comm, const char opt[], const char name[], const char idxname[], PetscInt n, const PetscSFNode a[])
875: {
876:   PetscInt    idx;
877:   PetscMPIInt rank;
878:   PetscBool   flg;

880:   PetscFunctionBegin;
881:   PetscCall(PetscOptionsHasName(NULL, NULL, opt, &flg));
882:   if (!flg) PetscFunctionReturn(PETSC_SUCCESS);
883:   PetscCallMPI(MPI_Comm_rank(comm, &rank));
884:   PetscCall(PetscSynchronizedPrintf(comm, "[%d]%s:\n", rank, name));
885:   if (idxname) {
886:     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));
887:   } else {
888:     for (idx = 0; idx < n; ++idx) PetscCall(PetscSynchronizedPrintf(comm, "[%d]rank %" PetscInt_FMT " index %" PetscInt_FMT "\n", rank, a[idx].rank, a[idx].index));
889:   }
890:   PetscCall(PetscSynchronizedFlush(comm, NULL));
891:   PetscFunctionReturn(PETSC_SUCCESS);
892: }

894: static PetscErrorCode DMPlexMapToLocalPoint(DM dm, PetscHMapIJ remotehash, PetscSFNode remotePoint, PetscInt *localPoint, PetscBool *mapFailed)
895: {
896:   PetscSF         sf;
897:   const PetscInt *locals;
898:   PetscMPIInt     rank;

900:   PetscFunctionBegin;
901:   PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank));
902:   PetscCall(DMGetPointSF(dm, &sf));
903:   PetscCall(PetscSFGetGraph(sf, NULL, NULL, &locals, NULL));
904:   if (mapFailed) *mapFailed = PETSC_FALSE;
905:   if (remotePoint.rank == rank) {
906:     *localPoint = remotePoint.index;
907:   } else {
908:     PetscHashIJKey key;
909:     PetscInt       l;

911:     key.i = remotePoint.index;
912:     key.j = remotePoint.rank;
913:     PetscCall(PetscHMapIJGet(remotehash, key, &l));
914:     if (l >= 0) {
915:       *localPoint = locals[l];
916:     } else if (mapFailed) *mapFailed = PETSC_TRUE;
917:   }
918:   PetscFunctionReturn(PETSC_SUCCESS);
919: }

921: static PetscErrorCode DMPlexMapToGlobalPoint(DM dm, PetscInt localPoint, PetscSFNode *remotePoint, PetscBool *mapFailed)
922: {
923:   PetscSF            sf;
924:   const PetscInt    *locals, *rootdegree;
925:   const PetscSFNode *remotes;
926:   PetscInt           Nl, l;
927:   PetscMPIInt        rank;

929:   PetscFunctionBegin;
930:   if (mapFailed) *mapFailed = PETSC_FALSE;
931:   PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank));
932:   PetscCall(DMGetPointSF(dm, &sf));
933:   PetscCall(PetscSFGetGraph(sf, NULL, &Nl, &locals, &remotes));
934:   if (Nl < 0) goto owned;
935:   PetscCall(PetscSFComputeDegreeBegin(sf, &rootdegree));
936:   PetscCall(PetscSFComputeDegreeEnd(sf, &rootdegree));
937:   if (rootdegree[localPoint]) goto owned;
938:   PetscCall(PetscFindInt(localPoint, Nl, locals, &l));
939:   if (l < 0) {
940:     if (mapFailed) *mapFailed = PETSC_TRUE;
941:   } else *remotePoint = remotes[l];
942:   PetscFunctionReturn(PETSC_SUCCESS);
943: owned:
944:   remotePoint->rank  = rank;
945:   remotePoint->index = localPoint;
946:   PetscFunctionReturn(PETSC_SUCCESS);
947: }

949: static PetscErrorCode DMPlexPointIsShared(DM dm, PetscInt p, PetscBool *isShared)
950: {
951:   PetscSF         sf;
952:   const PetscInt *locals, *rootdegree;
953:   PetscInt        Nl, idx;

955:   PetscFunctionBegin;
956:   *isShared = PETSC_FALSE;
957:   PetscCall(DMGetPointSF(dm, &sf));
958:   PetscCall(PetscSFGetGraph(sf, NULL, &Nl, &locals, NULL));
959:   if (Nl < 0) PetscFunctionReturn(PETSC_SUCCESS);
960:   PetscCall(PetscFindInt(p, Nl, locals, &idx));
961:   if (idx >= 0) {
962:     *isShared = PETSC_TRUE;
963:     PetscFunctionReturn(PETSC_SUCCESS);
964:   }
965:   PetscCall(PetscSFComputeDegreeBegin(sf, &rootdegree));
966:   PetscCall(PetscSFComputeDegreeEnd(sf, &rootdegree));
967:   if (rootdegree[p] > 0) *isShared = PETSC_TRUE;
968:   PetscFunctionReturn(PETSC_SUCCESS);
969: }

971: static PetscErrorCode DMPlexConeIsShared(DM dm, PetscInt p, PetscBool *isShared)
972: {
973:   const PetscInt *cone;
974:   PetscInt        coneSize, c;
975:   PetscBool       cShared = PETSC_TRUE;

977:   PetscFunctionBegin;
978:   PetscCall(DMPlexGetConeSize(dm, p, &coneSize));
979:   PetscCall(DMPlexGetCone(dm, p, &cone));
980:   for (c = 0; c < coneSize; ++c) {
981:     PetscBool pointShared;

983:     PetscCall(DMPlexPointIsShared(dm, cone[c], &pointShared));
984:     cShared = (PetscBool)(cShared && pointShared);
985:   }
986:   *isShared = coneSize ? cShared : PETSC_FALSE;
987:   PetscFunctionReturn(PETSC_SUCCESS);
988: }

990: static PetscErrorCode DMPlexGetConeMinimum(DM dm, PetscInt p, PetscSFNode *cpmin)
991: {
992:   const PetscInt *cone;
993:   PetscInt        coneSize, c;
994:   PetscSFNode     cmin = {PETSC_MAX_INT, PETSC_MAX_INT}, missing = {-1, -1};

996:   PetscFunctionBegin;
997:   PetscCall(DMPlexGetConeSize(dm, p, &coneSize));
998:   PetscCall(DMPlexGetCone(dm, p, &cone));
999:   for (c = 0; c < coneSize; ++c) {
1000:     PetscSFNode rcp;
1001:     PetscBool   mapFailed;

1003:     PetscCall(DMPlexMapToGlobalPoint(dm, cone[c], &rcp, &mapFailed));
1004:     if (mapFailed) {
1005:       cmin = missing;
1006:     } else {
1007:       cmin = (rcp.rank < cmin.rank) || (rcp.rank == cmin.rank && rcp.index < cmin.index) ? rcp : cmin;
1008:     }
1009:   }
1010:   *cpmin = coneSize ? cmin : missing;
1011:   PetscFunctionReturn(PETSC_SUCCESS);
1012: }

1014: /*
1015:   Each shared face has an entry in the candidates array:
1016:     (-1, coneSize-1), {(global cone point)}
1017:   where the set is missing the point p which we use as the key for the face
1018: */
1019: static PetscErrorCode DMPlexAddSharedFace_Private(DM dm, PetscSection candidateSection, PetscSFNode candidates[], PetscHMapIJ faceHash, PetscInt p, PetscBool debug)
1020: {
1021:   MPI_Comm        comm;
1022:   const PetscInt *support;
1023:   PetscInt        supportSize, s, off = 0, idx = 0, overlap, cellHeight, height;
1024:   PetscMPIInt     rank;

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

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

1074:           if (cp == p) continue;
1075:           PetscCall(DMPlexMapToGlobalPoint(dm, cp, &candidates[off + idx], NULL));
1076:           if (debug) PetscCall(PetscSynchronizedPrintf(comm, " (%" PetscInt_FMT ",%" PetscInt_FMT ")", candidates[off + idx].rank, candidates[off + idx].index));
1077:           ++idx;
1078:         }
1079:         if (debug) PetscCall(PetscSynchronizedPrintf(comm, "\n"));
1080:       } else {
1081:         /* Add cone size to section */
1082:         if (debug) PetscCall(PetscSynchronizedPrintf(comm, "[%d]    Scheduling shared face %" PetscInt_FMT "\n", rank, face));
1083:         PetscCall(DMPlexGetConeSize(dm, face, &coneSize));
1084:         PetscCall(PetscHMapIJSet(faceHash, key, p));
1085:         PetscCall(PetscSectionAddDof(candidateSection, p, coneSize + 1));
1086:       }
1087:     }
1088:   }
1089:   PetscFunctionReturn(PETSC_SUCCESS);
1090: }

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

1095:   Collective

1097:   Input Parameters:
1098: + dm      - The interpolated `DMPLEX`
1099: - pointSF - The initial `PetscSF` without interpolated points

1101:   Level: developer

1103:    Note:
1104:    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`

1106: .seealso: `DMPLEX`, `DMPlexInterpolate()`, `DMPlexUninterpolate()`
1107: @*/
1108: PetscErrorCode DMPlexInterpolatePointSF(DM dm, PetscSF pointSF)
1109: {
1110:   MPI_Comm           comm;
1111:   PetscHMapIJ        remoteHash;
1112:   PetscHMapI         claimshash;
1113:   PetscSection       candidateSection, candidateRemoteSection, claimSection;
1114:   PetscSFNode       *candidates, *candidatesRemote, *claims;
1115:   const PetscInt    *localPoints, *rootdegree;
1116:   const PetscSFNode *remotePoints;
1117:   PetscInt           ov, Nr, r, Nl, l;
1118:   PetscInt           candidatesSize, candidatesRemoteSize, claimsSize;
1119:   PetscBool          flg, debug = PETSC_FALSE;
1120:   PetscMPIInt        rank;

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

1171:     PetscCall(PetscHMapIJCreate(&faceHash));
1172:     for (l = 0; l < Nl; ++l) {
1173:       const PetscInt p = localPoints[l];

1175:       if (debug) PetscCall(PetscSynchronizedPrintf(comm, "[%d]  First pass leaf point %" PetscInt_FMT "\n", rank, p));
1176:       PetscCall(DMPlexAddSharedFace_Private(dm, candidateSection, NULL, faceHash, p, debug));
1177:     }
1178:     PetscCall(PetscHMapIJClear(faceHash));
1179:     PetscCall(PetscSectionSetUp(candidateSection));
1180:     PetscCall(PetscSectionGetStorageSize(candidateSection, &candidatesSize));
1181:     PetscCall(PetscMalloc1(candidatesSize, &candidates));
1182:     for (l = 0; l < Nl; ++l) {
1183:       const PetscInt p = localPoints[l];

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

1200:     PetscCall(PetscSFGetMultiSF(pointSF, &sfMulti));
1201:     PetscCall(PetscSFCreateInverseSF(sfMulti, &sfInverse));
1202:     PetscCall(PetscSectionCreate(comm, &candidateRemoteSection));
1203:     PetscCall(PetscSFDistributeSection(sfInverse, candidateSection, &remoteOffsets, candidateRemoteSection));
1204:     PetscCall(PetscSFCreateSectionSF(sfInverse, candidateSection, remoteOffsets, candidateRemoteSection, &sfCandidates));
1205:     PetscCall(PetscSectionGetStorageSize(candidateRemoteSection, &candidatesRemoteSize));
1206:     PetscCall(PetscMalloc1(candidatesRemoteSize, &candidatesRemote));
1207:     PetscCall(PetscSFBcastBegin(sfCandidates, MPIU_2INT, candidates, candidatesRemote, MPI_REPLACE));
1208:     PetscCall(PetscSFBcastEnd(sfCandidates, MPIU_2INT, candidates, candidatesRemote, MPI_REPLACE));
1209:     PetscCall(PetscSFDestroy(&sfInverse));
1210:     PetscCall(PetscSFDestroy(&sfCandidates));
1211:     PetscCall(PetscFree(remoteOffsets));

1213:     PetscCall(PetscObjectSetName((PetscObject)candidateRemoteSection, "Remote Candidate Section"));
1214:     PetscCall(PetscObjectViewFromOptions((PetscObject)candidateRemoteSection, NULL, "-petscsection_interp_candidate_remote_view"));
1215:     PetscCall(SFNodeArrayViewFromOptions(comm, "-petscsection_interp_candidate_remote_view", "Remote Candidates", NULL, candidatesRemoteSize, candidatesRemote));
1216:   }
1217:   /* 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 */
1218:   {
1219:     PetscHashIJKLRemote faceTable;
1220:     PetscInt            idx, idx2;

1222:     PetscCall(PetscHashIJKLRemoteCreate(&faceTable));
1223:     /* There is a section point for every leaf attached to a given root point */
1224:     for (r = 0, idx = 0, idx2 = 0; r < Nr; ++r) {
1225:       PetscInt deg;

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

1230:         PetscCall(PetscSectionGetDof(candidateRemoteSection, idx, &dof));
1231:         PetscCall(PetscSectionGetOffset(candidateRemoteSection, idx, &offset));
1232:         /* dof may include many faces from the remote process */
1233:         for (d = 0; d < dof; ++d) {
1234:           const PetscInt         hidx  = offset + d;
1235:           const PetscInt         Np    = candidatesRemote[hidx].index + 1;
1236:           const PetscSFNode      rface = candidatesRemote[hidx + 1];
1237:           const PetscSFNode     *fcone = &candidatesRemote[hidx + 2];
1238:           PetscSFNode            fcp0;
1239:           const PetscSFNode      pmax = {PETSC_MAX_INT, PETSC_MAX_INT};
1240:           const PetscInt        *join = NULL;
1241:           PetscHashIJKLRemoteKey key;
1242:           PetscHashIter          iter;
1243:           PetscBool              missing, mapToLocalPointFailed = PETSC_FALSE;
1244:           PetscInt               points[1024], p, joinSize;

1246:           if (debug)
1247:             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,
1248:                                               rface.index, r, idx, d, Np));
1249:           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);
1250:           fcp0.rank  = rank;
1251:           fcp0.index = r;
1252:           d += Np;
1253:           /* Put remote face in hash table */
1254:           key.i = fcp0;
1255:           key.j = fcone[0];
1256:           key.k = Np > 2 ? fcone[1] : pmax;
1257:           key.l = Np > 3 ? fcone[2] : pmax;
1258:           PetscCall(PetscSortSFNode(Np, (PetscSFNode *)&key));
1259:           PetscCall(PetscHashIJKLRemotePut(faceTable, key, &iter, &missing));
1260:           if (missing) {
1261:             if (debug) PetscCall(PetscSynchronizedPrintf(PetscObjectComm((PetscObject)dm), "[%d]  Setting remote face (%" PetscInt_FMT ", %" PetscInt_FMT ")\n", rank, rface.index, rface.rank));
1262:             PetscCall(PetscHashIJKLRemoteIterSet(faceTable, iter, rface));
1263:           } else {
1264:             PetscSFNode oface;

1266:             PetscCall(PetscHashIJKLRemoteIterGet(faceTable, iter, &oface));
1267:             if ((rface.rank < oface.rank) || (rface.rank == oface.rank && rface.index < oface.index)) {
1268:               if (debug) PetscCall(PetscSynchronizedPrintf(PetscObjectComm((PetscObject)dm), "[%d]  Replacing with remote face (%" PetscInt_FMT ", %" PetscInt_FMT ")\n", rank, rface.index, rface.rank));
1269:               PetscCall(PetscHashIJKLRemoteIterSet(faceTable, iter, rface));
1270:             }
1271:           }
1272:           /* Check for local face */
1273:           points[0] = r;
1274:           for (p = 1; p < Np; ++p) {
1275:             PetscCall(DMPlexMapToLocalPoint(dm, remoteHash, fcone[p - 1], &points[p], &mapToLocalPointFailed));
1276:             if (mapToLocalPointFailed) break; /* We got a point not in our overlap */
1277:             if (debug) PetscCall(PetscSynchronizedPrintf(PetscObjectComm((PetscObject)dm), "[%d]  Checking local candidate %" PetscInt_FMT "\n", rank, points[p]));
1278:           }
1279:           if (mapToLocalPointFailed) continue;
1280:           PetscCall(DMPlexGetJoin(dm, Np, points, &joinSize, &join));
1281:           if (joinSize == 1) {
1282:             PetscSFNode lface;
1283:             PetscSFNode oface;

1285:             /* Always replace with local face */
1286:             lface.rank  = rank;
1287:             lface.index = join[0];
1288:             PetscCall(PetscHashIJKLRemoteIterGet(faceTable, iter, &oface));
1289:             if (debug)
1290:               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));
1291:             PetscCall(PetscHashIJKLRemoteIterSet(faceTable, iter, lface));
1292:           }
1293:           PetscCall(DMPlexRestoreJoin(dm, Np, points, &joinSize, &join));
1294:         }
1295:       }
1296:       /* Put back faces for this root */
1297:       for (deg = 0; deg < rootdegree[r]; ++deg, ++idx2) {
1298:         PetscInt offset, dof, d;

1300:         PetscCall(PetscSectionGetDof(candidateRemoteSection, idx2, &dof));
1301:         PetscCall(PetscSectionGetOffset(candidateRemoteSection, idx2, &offset));
1302:         /* dof may include many faces from the remote process */
1303:         for (d = 0; d < dof; ++d) {
1304:           const PetscInt         hidx  = offset + d;
1305:           const PetscInt         Np    = candidatesRemote[hidx].index + 1;
1306:           const PetscSFNode     *fcone = &candidatesRemote[hidx + 2];
1307:           PetscSFNode            fcp0;
1308:           const PetscSFNode      pmax = {PETSC_MAX_INT, PETSC_MAX_INT};
1309:           PetscHashIJKLRemoteKey key;
1310:           PetscHashIter          iter;
1311:           PetscBool              missing;

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

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

1364:       if (debug) PetscCall(PetscSynchronizedPrintf(comm, "[%d]  Checking root for claims %" PetscInt_FMT "\n", rank, r));
1365:       PetscCall(PetscSectionGetDof(candidateSection, r, &dof));
1366:       PetscCall(PetscSectionGetOffset(candidateSection, r, &off));
1367:       for (d = 0; d < dof;) {
1368:         if (claims[off + d].rank >= 0) {
1369:           const PetscInt  faceInd = off + d;
1370:           const PetscInt  Np      = candidates[off + d].index;
1371:           const PetscInt *join    = NULL;
1372:           PetscInt        joinSize, points[1024], c;

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

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

1443:   Collective

1445:   Input Parameter:
1446: . dm - The `DMPLEX` object with only cells and vertices

1448:   Output Parameter:
1449: . dmInt - The complete `DMPLEX` object

1451:   Level: intermediate

1453:   Note:
1454:     Labels and coordinates are copied.

1456:   Developer Note:
1457:     It sets plex->interpolated = `DMPLEX_INTERPOLATED_FULL`.

1459: .seealso: `DMPLEX`, `DMPlexUninterpolate()`, `DMPlexCreateFromCellListPetsc()`, `DMPlexCopyCoordinates()`
1460: @*/
1461: PetscErrorCode DMPlexInterpolate(DM dm, DM *dmInt)
1462: {
1463:   DMPlexInterpolatedFlag interpolated;
1464:   DM                     idm, odm = dm;
1465:   PetscSF                sfPoint;
1466:   PetscInt               depth, dim, d;
1467:   const char            *name;
1468:   PetscBool              flg = PETSC_TRUE;

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

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

1522:   Collective

1524:   Input Parameter:
1525: . dmA - The `DMPLEX` object with initial coordinates

1527:   Output Parameter:
1528: . dmB - The `DMPLEX` object with copied coordinates

1530:   Level: intermediate

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

1535: .seealso: `DMPLEX`, `DMCopyLabels()`, `DMGetCoordinates()`, `DMGetCoordinatesLocal()`, `DMGetCoordinateDM()`, `DMGetCoordinateSection()`
1536: @*/
1537: PetscErrorCode DMPlexCopyCoordinates(DM dmA, DM dmB)
1538: {
1539:   Vec          coordinatesA, coordinatesB;
1540:   VecType      vtype;
1541:   PetscSection coordSectionA, coordSectionB;
1542:   PetscScalar *coordsA, *coordsB;
1543:   PetscInt     spaceDim, Nf, vStartA, vStartB, vEndA, vEndB, coordSizeB, v, d;
1544:   PetscInt     cStartA, cEndA, cStartB, cEndB, cS, cE, cdim;
1545:   PetscBool    lc = PETSC_FALSE;

1547:   PetscFunctionBegin;
1550:   if (dmA == dmB) PetscFunctionReturn(PETSC_SUCCESS);
1551:   PetscCall(DMGetCoordinateDim(dmA, &cdim));
1552:   PetscCall(DMSetCoordinateDim(dmB, cdim));
1553:   PetscCall(DMPlexGetDepthStratum(dmA, 0, &vStartA, &vEndA));
1554:   PetscCall(DMPlexGetDepthStratum(dmB, 0, &vStartB, &vEndB));
1555:   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);
1556:   /* Copy over discretization if it exists */
1557:   {
1558:     DM                 cdmA, cdmB;
1559:     PetscDS            dsA, dsB;
1560:     PetscObject        objA, objB;
1561:     PetscClassId       idA, idB;
1562:     const PetscScalar *constants;
1563:     PetscInt           cdim, Nc;

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

1593:     PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)coordSectionA), &coordSectionB));
1594:     PetscCall(DMGetCoordinateDim(dmA, &dim));
1595:     PetscCall(DMSetCoordinateSection(dmB, dim, coordSectionB));
1596:     PetscCall(PetscObjectDereference((PetscObject)coordSectionB));
1597:   }
1598:   PetscCall(PetscSectionSetNumFields(coordSectionB, 1));
1599:   PetscCall(PetscSectionGetFieldComponents(coordSectionA, 0, &spaceDim));
1600:   PetscCall(PetscSectionSetFieldComponents(coordSectionB, 0, spaceDim));
1601:   PetscCall(PetscSectionGetChart(coordSectionA, &cS, &cE));
1602:   if (cStartA <= cS && cS < cEndA) { /* localized coordinates */
1603:     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);
1604:     cS = cS - cStartA + cStartB;
1605:     cE = vEndB;
1606:     lc = PETSC_TRUE;
1607:   } else {
1608:     cS = vStartB;
1609:     cE = vEndB;
1610:   }
1611:   PetscCall(PetscSectionSetChart(coordSectionB, cS, cE));
1612:   for (v = vStartB; v < vEndB; ++v) {
1613:     PetscCall(PetscSectionSetDof(coordSectionB, v, spaceDim));
1614:     PetscCall(PetscSectionSetFieldDof(coordSectionB, v, 0, spaceDim));
1615:   }
1616:   if (lc) { /* localized coordinates */
1617:     PetscInt c;

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

1622:       PetscCall(PetscSectionGetDof(coordSectionA, c + cStartA, &dof));
1623:       PetscCall(PetscSectionSetDof(coordSectionB, c + cStartB, dof));
1624:       PetscCall(PetscSectionSetFieldDof(coordSectionB, c + cStartB, 0, dof));
1625:     }
1626:   }
1627:   PetscCall(PetscSectionSetUp(coordSectionB));
1628:   PetscCall(PetscSectionGetStorageSize(coordSectionB, &coordSizeB));
1629:   PetscCall(DMGetCoordinatesLocal(dmA, &coordinatesA));
1630:   PetscCall(VecCreate(PETSC_COMM_SELF, &coordinatesB));
1631:   PetscCall(PetscObjectSetName((PetscObject)coordinatesB, "coordinates"));
1632:   PetscCall(VecSetSizes(coordinatesB, coordSizeB, PETSC_DETERMINE));
1633:   PetscCall(VecGetBlockSize(coordinatesA, &d));
1634:   PetscCall(VecSetBlockSize(coordinatesB, d));
1635:   PetscCall(VecGetType(coordinatesA, &vtype));
1636:   PetscCall(VecSetType(coordinatesB, vtype));
1637:   PetscCall(VecGetArray(coordinatesA, &coordsA));
1638:   PetscCall(VecGetArray(coordinatesB, &coordsB));
1639:   for (v = 0; v < vEndB - vStartB; ++v) {
1640:     PetscInt offA, offB;

1642:     PetscCall(PetscSectionGetOffset(coordSectionA, v + vStartA, &offA));
1643:     PetscCall(PetscSectionGetOffset(coordSectionB, v + vStartB, &offB));
1644:     for (d = 0; d < spaceDim; ++d) coordsB[offB + d] = coordsA[offA + d];
1645:   }
1646:   if (lc) { /* localized coordinates */
1647:     PetscInt c;

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

1652:       PetscCall(PetscSectionGetOffset(coordSectionA, c + cStartA, &offA));
1653:       PetscCall(PetscSectionGetOffset(coordSectionB, c + cStartB, &offB));
1654:       PetscCall(PetscSectionGetDof(coordSectionA, c + cStartA, &dof));
1655:       PetscCall(PetscArraycpy(coordsB + offB, coordsA + offA, dof));
1656:     }
1657:   }
1658:   PetscCall(VecRestoreArray(coordinatesA, &coordsA));
1659:   PetscCall(VecRestoreArray(coordinatesB, &coordsB));
1660:   PetscCall(DMSetCoordinatesLocal(dmB, coordinatesB));
1661:   PetscCall(VecDestroy(&coordinatesB));
1662:   PetscFunctionReturn(PETSC_SUCCESS);
1663: }

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

1668:   Collective

1670:   Input Parameter:
1671: . dm - The complete `DMPLEX` object

1673:   Output Parameter:
1674: . dmUnint - The `DMPLEX` object with only cells and vertices

1676:   Level: intermediate

1678:   Note:
1679:     It does not copy over the coordinates.

1681:   Developer Note:
1682:   Sets plex->interpolated = `DMPLEX_INTERPOLATED_NONE`.

1684: .seealso: `DMPLEX`, `DMPlexInterpolate()`, `DMPlexCreateFromCellListPetsc()`, `DMPlexCopyCoordinates()`
1685: @*/
1686: PetscErrorCode DMPlexUninterpolate(DM dm, DM *dmUnint)
1687: {
1688:   DMPlexInterpolatedFlag interpolated;
1689:   DM                     udm;
1690:   PetscInt               dim, vStart, vEnd, cStart, cEnd, c, maxConeSize = 0, *cone;

1692:   PetscFunctionBegin;
1695:   PetscCall(DMGetDimension(dm, &dim));
1696:   PetscCall(DMPlexIsInterpolated(dm, &interpolated));
1697:   PetscCheck(interpolated != DMPLEX_INTERPOLATED_PARTIAL, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Not for partially interpolated meshes");
1698:   if (interpolated == DMPLEX_INTERPOLATED_NONE || dim <= 1) {
1699:     /* in case dim <= 1 just keep the DMPLEX_INTERPOLATED_FULL flag */
1700:     PetscCall(PetscObjectReference((PetscObject)dm));
1701:     *dmUnint = dm;
1702:     PetscFunctionReturn(PETSC_SUCCESS);
1703:   }
1704:   PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd));
1705:   PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
1706:   PetscCall(DMCreate(PetscObjectComm((PetscObject)dm), &udm));
1707:   PetscCall(DMSetType(udm, DMPLEX));
1708:   PetscCall(DMSetDimension(udm, dim));
1709:   PetscCall(DMPlexSetChart(udm, cStart, vEnd));
1710:   for (c = cStart; c < cEnd; ++c) {
1711:     PetscInt *closure = NULL, closureSize, cl, coneSize = 0;

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

1717:       if ((p >= vStart) && (p < vEnd)) ++coneSize;
1718:     }
1719:     PetscCall(DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure));
1720:     PetscCall(DMPlexSetConeSize(udm, c, coneSize));
1721:     maxConeSize = PetscMax(maxConeSize, coneSize);
1722:   }
1723:   PetscCall(DMSetUp(udm));
1724:   PetscCall(PetscMalloc1(maxConeSize, &cone));
1725:   for (c = cStart; c < cEnd; ++c) {
1726:     PetscInt *closure = NULL, closureSize, cl, coneSize = 0;

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

1732:       if ((p >= vStart) && (p < vEnd)) cone[coneSize++] = p;
1733:     }
1734:     PetscCall(DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure));
1735:     PetscCall(DMPlexSetCone(udm, c, cone));
1736:   }
1737:   PetscCall(PetscFree(cone));
1738:   PetscCall(DMPlexSymmetrize(udm));
1739:   PetscCall(DMPlexStratify(udm));
1740:   /* Reduce SF */
1741:   {
1742:     PetscSF            sfPoint, sfPointUn;
1743:     const PetscSFNode *remotePoints;
1744:     const PetscInt    *localPoints;
1745:     PetscSFNode       *remotePointsUn;
1746:     PetscInt          *localPointsUn;
1747:     PetscInt           numRoots, numLeaves, l;
1748:     PetscInt           numLeavesUn = 0, n = 0;

1750:     /* Get original SF information */
1751:     PetscCall(DMGetPointSF(dm, &sfPoint));
1752:     if (PetscDefined(USE_DEBUG)) PetscCall(DMPlexCheckPointSF(dm, sfPoint, PETSC_FALSE));
1753:     PetscCall(DMGetPointSF(udm, &sfPointUn));
1754:     PetscCall(PetscSFGetGraph(sfPoint, &numRoots, &numLeaves, &localPoints, &remotePoints));
1755:     if (numRoots >= 0) {
1756:       /* Allocate space for cells and vertices */
1757:       for (l = 0; l < numLeaves; ++l) {
1758:         const PetscInt p = localPoints[l];

1760:         if ((vStart <= p && p < vEnd) || (cStart <= p && p < cEnd)) numLeavesUn++;
1761:       }
1762:       /* Fill in leaves */
1763:       PetscCall(PetscMalloc1(numLeavesUn, &remotePointsUn));
1764:       PetscCall(PetscMalloc1(numLeavesUn, &localPointsUn));
1765:       for (l = 0; l < numLeaves; l++) {
1766:         const PetscInt p = localPoints[l];

1768:         if ((vStart <= p && p < vEnd) || (cStart <= p && p < cEnd)) {
1769:           localPointsUn[n]        = p;
1770:           remotePointsUn[n].rank  = remotePoints[l].rank;
1771:           remotePointsUn[n].index = remotePoints[l].index;
1772:           ++n;
1773:         }
1774:       }
1775:       PetscCheck(n == numLeavesUn, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Inconsistent number of leaves %" PetscInt_FMT " != %" PetscInt_FMT, n, numLeavesUn);
1776:       PetscCall(PetscSFSetGraph(sfPointUn, cEnd - cStart + vEnd - vStart, numLeavesUn, localPointsUn, PETSC_OWN_POINTER, remotePointsUn, PETSC_OWN_POINTER));
1777:     }
1778:   }
1779:   /* This function makes the mesh fully uninterpolated on all ranks */
1780:   {
1781:     DM_Plex *plex      = (DM_Plex *)udm->data;
1782:     plex->interpolated = plex->interpolatedCollective = DMPLEX_INTERPOLATED_NONE;
1783:   }
1784:   PetscCall(DMPlexCopy_Internal(dm, PETSC_TRUE, PETSC_TRUE, udm));
1785:   if (PetscDefined(USE_DEBUG)) PetscCall(DMPlexCheckPointSF(udm, NULL, PETSC_FALSE));
1786:   *dmUnint = udm;
1787:   PetscFunctionReturn(PETSC_SUCCESS);
1788: }

1790: static PetscErrorCode DMPlexIsInterpolated_Internal(DM dm, DMPlexInterpolatedFlag *interpolated)
1791: {
1792:   PetscInt coneSize, depth, dim, h, p, pStart, pEnd;
1793:   MPI_Comm comm;

1795:   PetscFunctionBegin;
1796:   PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
1797:   PetscCall(DMPlexGetDepth(dm, &depth));
1798:   PetscCall(DMGetDimension(dm, &dim));

1800:   if (depth == dim) {
1801:     *interpolated = DMPLEX_INTERPOLATED_FULL;
1802:     if (!dim) goto finish;

1804:     /* Check points at height = dim are vertices (have no cones) */
1805:     PetscCall(DMPlexGetHeightStratum(dm, dim, &pStart, &pEnd));
1806:     for (p = pStart; p < pEnd; p++) {
1807:       PetscCall(DMPlexGetConeSize(dm, p, &coneSize));
1808:       if (coneSize) {
1809:         *interpolated = DMPLEX_INTERPOLATED_PARTIAL;
1810:         goto finish;
1811:       }
1812:     }

1814:     /* Check points at height < dim have cones */
1815:     for (h = 0; h < dim; h++) {
1816:       PetscCall(DMPlexGetHeightStratum(dm, h, &pStart, &pEnd));
1817:       for (p = pStart; p < pEnd; p++) {
1818:         PetscCall(DMPlexGetConeSize(dm, p, &coneSize));
1819:         if (!coneSize) {
1820:           *interpolated = DMPLEX_INTERPOLATED_PARTIAL;
1821:           goto finish;
1822:         }
1823:       }
1824:     }
1825:   } else if (depth == 1) {
1826:     *interpolated = DMPLEX_INTERPOLATED_NONE;
1827:   } else {
1828:     *interpolated = DMPLEX_INTERPOLATED_PARTIAL;
1829:   }
1830: finish:
1831:   PetscFunctionReturn(PETSC_SUCCESS);
1832: }

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

1837:   Not Collective

1839:   Input Parameter:
1840: . dm      - The `DMPLEX` object

1842:   Output Parameter:
1843: . interpolated - Flag whether the `DM` is interpolated

1845:   Level: intermediate

1847:   Notes:
1848:   Unlike `DMPlexIsInterpolatedCollective()`, this is NOT collective
1849:   so the results can be different on different ranks in special cases.
1850:   However, `DMPlexInterpolate()` guarantees the result is the same on all.

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

1854:   Developer Notes:
1855:   Initially, plex->interpolated = `DMPLEX_INTERPOLATED_INVALID`.

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

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

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

1866: .seealso: `DMPLEX`, `DMPlexInterpolate()`, `DMPlexIsInterpolatedCollective()`
1867: @*/
1868: PetscErrorCode DMPlexIsInterpolated(DM dm, DMPlexInterpolatedFlag *interpolated)
1869: {
1870:   DM_Plex *plex = (DM_Plex *)dm->data;

1872:   PetscFunctionBegin;
1875:   if (plex->interpolated < 0) {
1876:     PetscCall(DMPlexIsInterpolated_Internal(dm, &plex->interpolated));
1877:   } else if (PetscDefined(USE_DEBUG)) {
1878:     DMPlexInterpolatedFlag flg;

1880:     PetscCall(DMPlexIsInterpolated_Internal(dm, &flg));
1881:     PetscCheck(flg == plex->interpolated, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Stashed DMPlexInterpolatedFlag %s is inconsistent with current %s", DMPlexInterpolatedFlags[plex->interpolated], DMPlexInterpolatedFlags[flg]);
1882:   }
1883:   *interpolated = plex->interpolated;
1884:   PetscFunctionReturn(PETSC_SUCCESS);
1885: }

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

1890:   Collective

1892:   Input Parameter:
1893: . dm      - The `DMPLEX` object

1895:   Output Parameter:
1896: . interpolated - Flag whether the `DM` is interpolated

1898:   Level: intermediate

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

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

1905:   Developer Notes:
1906:   Initially, plex->interpolatedCollective = `DMPLEX_INTERPOLATED_INVALID`.

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

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

1915: .seealso: `DMPLEX`, `DMPlexInterpolate()`, `DMPlexIsInterpolated()`
1916: @*/
1917: PetscErrorCode DMPlexIsInterpolatedCollective(DM dm, DMPlexInterpolatedFlag *interpolated)
1918: {
1919:   DM_Plex  *plex  = (DM_Plex *)dm->data;
1920:   PetscBool debug = PETSC_FALSE;

1922:   PetscFunctionBegin;
1925:   PetscCall(PetscOptionsGetBool(((PetscObject)dm)->options, ((PetscObject)dm)->prefix, "-dm_plex_is_interpolated_collective_debug", &debug, NULL));
1926:   if (plex->interpolatedCollective < 0) {
1927:     DMPlexInterpolatedFlag min, max;
1928:     MPI_Comm               comm;

1930:     PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
1931:     PetscCall(DMPlexIsInterpolated(dm, &plex->interpolatedCollective));
1932:     PetscCall(MPIU_Allreduce(&plex->interpolatedCollective, &min, 1, MPIU_ENUM, MPI_MIN, comm));
1933:     PetscCall(MPIU_Allreduce(&plex->interpolatedCollective, &max, 1, MPIU_ENUM, MPI_MAX, comm));
1934:     if (min != max) plex->interpolatedCollective = DMPLEX_INTERPOLATED_MIXED;
1935:     if (debug) {
1936:       PetscMPIInt rank;

1938:       PetscCallMPI(MPI_Comm_rank(comm, &rank));
1939:       PetscCall(PetscSynchronizedPrintf(comm, "[%d] interpolated=%s interpolatedCollective=%s\n", rank, DMPlexInterpolatedFlags[plex->interpolated], DMPlexInterpolatedFlags[plex->interpolatedCollective]));
1940:       PetscCall(PetscSynchronizedFlush(comm, PETSC_STDOUT));
1941:     }
1942:   }
1943:   *interpolated = plex->interpolatedCollective;
1944:   PetscFunctionReturn(PETSC_SUCCESS);
1945: }