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 { PetscInt i, j, k, l; } PetscHashIJKLKey;

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

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

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

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

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

 26: #define PetscHashIJKLRemoteKeyHash(key) \
 27:   PetscHashCombine(PetscHashCombine(PetscHashInt((key).i.rank + (key).i.index),PetscHashInt((key).j.rank + (key).j.index)), \
 28:                    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;

 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:   return(0);
 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;
 61:   PetscErrorCode  ierr;

 66:   DMPlexGetMaxSizes(dm, &maxConeSize, &maxSupportSize);
 67:   if (faceTypes) {DMGetWorkArray(dm, PetscMax(maxConeSize, maxSupportSize),           MPIU_INT, &typesTmp);}
 68:   if (faceSizes) {DMGetWorkArray(dm, PetscMax(maxConeSize, maxSupportSize),           MPIU_INT, &sizesTmp);}
 69:   if (faces)     {DMGetWorkArray(dm, PetscSqr(PetscMax(maxConeSize, maxSupportSize)), MPIU_INT, &facesTmp);}
 70:   switch (ct) {
 71:     case DM_POLYTOPE_POINT:
 72:       if (numFaces) *numFaces = 0;
 73:       break;
 74:     case DM_POLYTOPE_SEGMENT:
 75:       if (numFaces) *numFaces = 2;
 76:       if (faceTypes) {
 77:         typesTmp[0] = DM_POLYTOPE_POINT; typesTmp[1] = DM_POLYTOPE_POINT;
 78:         *faceTypes = typesTmp;
 79:       }
 80:       if (faceSizes) {
 81:         sizesTmp[0] = 1; sizesTmp[1] = 1;
 82:         *faceSizes = sizesTmp;
 83:       }
 84:       if (faces) {
 85:         facesTmp[0] = cone[0]; facesTmp[1] = cone[1];
 86:         *faces = facesTmp;
 87:       }
 88:       break;
 89:     case DM_POLYTOPE_POINT_PRISM_TENSOR:
 90:       if (numFaces) *numFaces = 2;
 91:       if (faceTypes) {
 92:         typesTmp[0] = DM_POLYTOPE_POINT; typesTmp[1] = DM_POLYTOPE_POINT;
 93:         *faceTypes = typesTmp;
 94:       }
 95:       if (faceSizes) {
 96:         sizesTmp[0] = 1; sizesTmp[1] = 1;
 97:         *faceSizes = sizesTmp;
 98:       }
 99:       if (faces) {
100:         facesTmp[0] = cone[0]; facesTmp[1] = cone[1];
101:         *faces = facesTmp;
102:       }
103:       break;
104:     case DM_POLYTOPE_TRIANGLE:
105:       if (numFaces) *numFaces = 3;
106:       if (faceTypes) {
107:         typesTmp[0] = DM_POLYTOPE_SEGMENT; typesTmp[1] = DM_POLYTOPE_SEGMENT; typesTmp[2] = DM_POLYTOPE_SEGMENT;
108:         *faceTypes = typesTmp;
109:       }
110:       if (faceSizes) {
111:         sizesTmp[0] = 2; sizesTmp[1] = 2; sizesTmp[2] = 2;
112:         *faceSizes = sizesTmp;
113:       }
114:       if (faces) {
115:         facesTmp[0] = cone[0]; facesTmp[1] = cone[1];
116:         facesTmp[2] = cone[1]; facesTmp[3] = cone[2];
117:         facesTmp[4] = cone[2]; facesTmp[5] = cone[0];
118:         *faces = facesTmp;
119:       }
120:       break;
121:     case DM_POLYTOPE_QUADRILATERAL:
122:       /* Vertices follow right hand rule */
123:       if (numFaces) *numFaces = 4;
124:       if (faceTypes) {
125:         typesTmp[0] = DM_POLYTOPE_SEGMENT; typesTmp[1] = DM_POLYTOPE_SEGMENT; typesTmp[2] = DM_POLYTOPE_SEGMENT; typesTmp[3] = DM_POLYTOPE_SEGMENT;
126:         *faceTypes = typesTmp;
127:       }
128:       if (faceSizes) {
129:         sizesTmp[0] = 2; sizesTmp[1] = 2; sizesTmp[2] = 2; sizesTmp[3] = 2;
130:         *faceSizes = sizesTmp;
131:       }
132:       if (faces) {
133:         facesTmp[0] = cone[0]; facesTmp[1] = cone[1];
134:         facesTmp[2] = cone[1]; facesTmp[3] = cone[2];
135:         facesTmp[4] = cone[2]; facesTmp[5] = cone[3];
136:         facesTmp[6] = cone[3]; facesTmp[7] = cone[0];
137:         *faces = facesTmp;
138:       }
139:       break;
140:     case DM_POLYTOPE_SEG_PRISM_TENSOR:
141:       if (numFaces) *numFaces = 4;
142:       if (faceTypes) {
143:         typesTmp[0] = DM_POLYTOPE_SEGMENT; typesTmp[1] = DM_POLYTOPE_SEGMENT; typesTmp[2] = DM_POLYTOPE_POINT_PRISM_TENSOR; typesTmp[3] = DM_POLYTOPE_POINT_PRISM_TENSOR;
144:         *faceTypes = typesTmp;
145:       }
146:       if (faceSizes) {
147:         sizesTmp[0] = 2; sizesTmp[1] = 2; sizesTmp[2] = 2; sizesTmp[3] = 2;
148:         *faceSizes = sizesTmp;
149:       }
150:       if (faces) {
151:         facesTmp[0] = cone[0]; facesTmp[1] = cone[1];
152:         facesTmp[2] = cone[2]; facesTmp[3] = cone[3];
153:         facesTmp[4] = cone[0]; facesTmp[5] = cone[2];
154:         facesTmp[6] = cone[1]; facesTmp[7] = cone[3];
155:         *faces = facesTmp;
156:       }
157:       break;
158:     case DM_POLYTOPE_TETRAHEDRON:
159:       /* Vertices of first face follow right hand rule and normal points away from last vertex */
160:       if (numFaces) *numFaces = 4;
161:       if (faceTypes) {
162:         typesTmp[0] = DM_POLYTOPE_TRIANGLE; typesTmp[1] = DM_POLYTOPE_TRIANGLE; typesTmp[2] = DM_POLYTOPE_TRIANGLE; typesTmp[3] = DM_POLYTOPE_TRIANGLE;
163:         *faceTypes = typesTmp;
164:       }
165:       if (faceSizes) {
166:         sizesTmp[0] = 3; sizesTmp[1] = 3; sizesTmp[2] = 3; sizesTmp[3] = 3;
167:         *faceSizes = sizesTmp;
168:       }
169:       if (faces) {
170:         facesTmp[0] = cone[0]; facesTmp[1]  = cone[1]; facesTmp[2]  = cone[2];
171:         facesTmp[3] = cone[0]; facesTmp[4]  = cone[3]; facesTmp[5]  = cone[1];
172:         facesTmp[6] = cone[0]; facesTmp[7]  = cone[2]; facesTmp[8]  = cone[3];
173:         facesTmp[9] = cone[2]; facesTmp[10] = cone[1]; facesTmp[11] = cone[3];
174:         *faces = facesTmp;
175:       }
176:       break;
177:     case DM_POLYTOPE_HEXAHEDRON:
178:       /*  7--------6
179:          /|       /|
180:         / |      / |
181:        4--------5  |
182:        |  |     |  |
183:        |  |     |  |
184:        |  1--------2
185:        | /      | /
186:        |/       |/
187:        0--------3
188:        */
189:       if (numFaces) *numFaces = 6;
190:       if (faceTypes) {
191:         typesTmp[0] = DM_POLYTOPE_QUADRILATERAL; typesTmp[1] = DM_POLYTOPE_QUADRILATERAL; typesTmp[2] = DM_POLYTOPE_QUADRILATERAL;
192:         typesTmp[3] = DM_POLYTOPE_QUADRILATERAL; typesTmp[4] = DM_POLYTOPE_QUADRILATERAL; typesTmp[5] = DM_POLYTOPE_QUADRILATERAL;
193:         *faceTypes = typesTmp;
194:       }
195:       if (faceSizes) {
196:         sizesTmp[0] = 4; sizesTmp[1] = 4; sizesTmp[2] = 4; sizesTmp[3] = 4; sizesTmp[4] = 4; sizesTmp[5] = 4;
197:         *faceSizes = sizesTmp;
198:       }
199:       if (faces) {
200:         facesTmp[0]  = cone[0]; facesTmp[1]  = cone[1]; facesTmp[2]  = cone[2]; facesTmp[3]  = cone[3]; /* Bottom */
201:         facesTmp[4]  = cone[4]; facesTmp[5]  = cone[5]; facesTmp[6]  = cone[6]; facesTmp[7]  = cone[7]; /* Top */
202:         facesTmp[8]  = cone[0]; facesTmp[9]  = cone[3]; facesTmp[10] = cone[5]; facesTmp[11] = cone[4]; /* Front */
203:         facesTmp[12] = cone[2]; facesTmp[13] = cone[1]; facesTmp[14] = cone[7]; facesTmp[15] = cone[6]; /* Back */
204:         facesTmp[16] = cone[3]; facesTmp[17] = cone[2]; facesTmp[18] = cone[6]; facesTmp[19] = cone[5]; /* Right */
205:         facesTmp[20] = cone[0]; facesTmp[21] = cone[4]; facesTmp[22] = cone[7]; facesTmp[23] = cone[1]; /* Left */
206:         *faces = facesTmp;
207:       }
208:       break;
209:     case DM_POLYTOPE_TRI_PRISM:
210:       if (numFaces) *numFaces = 5;
211:       if (faceTypes) {
212:         typesTmp[0] = DM_POLYTOPE_TRIANGLE; typesTmp[1] = DM_POLYTOPE_TRIANGLE;
213:         typesTmp[2] = DM_POLYTOPE_QUADRILATERAL; typesTmp[3] = DM_POLYTOPE_QUADRILATERAL; typesTmp[4] = DM_POLYTOPE_QUADRILATERAL;
214:         *faceTypes = typesTmp;
215:       }
216:       if (faceSizes) {
217:         sizesTmp[0] = 3; sizesTmp[1] = 3;
218:         sizesTmp[2] = 4; sizesTmp[3] = 4; sizesTmp[4] = 4;
219:         *faceSizes = sizesTmp;
220:       }
221:       if (faces) {
222:         facesTmp[0]  = cone[0]; facesTmp[1]  = cone[1]; facesTmp[2]  = cone[2];                         /* Bottom */
223:         facesTmp[3]  = cone[3]; facesTmp[4]  = cone[4]; facesTmp[5]  = cone[5];                         /* Top */
224:         facesTmp[6]  = cone[0]; facesTmp[7]  = cone[2]; facesTmp[8]  = cone[4]; facesTmp[9]  = cone[3]; /* Back left */
225:         facesTmp[10] = cone[2]; facesTmp[11] = cone[1]; facesTmp[12] = cone[5]; facesTmp[13] = cone[4]; /* Front */
226:         facesTmp[14] = cone[1]; facesTmp[15] = cone[0]; facesTmp[16] = cone[3]; facesTmp[17] = cone[5]; /* Back right */
227:         *faces = facesTmp;
228:       }
229:       break;
230:     case DM_POLYTOPE_TRI_PRISM_TENSOR:
231:       if (numFaces)     *numFaces = 5;
232:       if (faceTypes) {
233:         typesTmp[0] = DM_POLYTOPE_TRIANGLE; typesTmp[1] = DM_POLYTOPE_TRIANGLE;
234:         typesTmp[2] = DM_POLYTOPE_SEG_PRISM_TENSOR; typesTmp[3] = DM_POLYTOPE_SEG_PRISM_TENSOR; typesTmp[4] = DM_POLYTOPE_SEG_PRISM_TENSOR;
235:         *faceTypes = typesTmp;
236:       }
237:       if (faceSizes) {
238:         sizesTmp[0] = 3; sizesTmp[1] = 3;
239:         sizesTmp[2] = 4; sizesTmp[3] = 4; sizesTmp[4] = 4;
240:         *faceSizes = sizesTmp;
241:       }
242:       if (faces) {
243:         facesTmp[0]  = cone[0]; facesTmp[1]  = cone[1]; facesTmp[2]  = cone[2];                         /* Bottom */
244:         facesTmp[3]  = cone[3]; facesTmp[4]  = cone[4]; facesTmp[5]  = cone[5];                         /* Top */
245:         facesTmp[6]  = cone[0]; facesTmp[7]  = cone[1]; facesTmp[8]  = cone[3]; facesTmp[9]  = cone[4]; /* Back left */
246:         facesTmp[10] = cone[1]; facesTmp[11] = cone[2]; facesTmp[12] = cone[4]; facesTmp[13] = cone[5]; /* Back right */
247:         facesTmp[14] = cone[2]; facesTmp[15] = cone[0]; facesTmp[16] = cone[5]; facesTmp[17] = cone[3]; /* Front */
248:         *faces = facesTmp;
249:       }
250:       break;
251:     case DM_POLYTOPE_QUAD_PRISM_TENSOR:
252:       /*  7--------6
253:          /|       /|
254:         / |      / |
255:        4--------5  |
256:        |  |     |  |
257:        |  |     |  |
258:        |  3--------2
259:        | /      | /
260:        |/       |/
261:        0--------1
262:        */
263:       if (numFaces) *numFaces = 6;
264:       if (faceTypes) {
265:         typesTmp[0] = DM_POLYTOPE_QUADRILATERAL;    typesTmp[1] = DM_POLYTOPE_QUADRILATERAL;    typesTmp[2] = DM_POLYTOPE_SEG_PRISM_TENSOR;
266:         typesTmp[3] = DM_POLYTOPE_SEG_PRISM_TENSOR; typesTmp[4] = DM_POLYTOPE_SEG_PRISM_TENSOR; typesTmp[5] = DM_POLYTOPE_SEG_PRISM_TENSOR;
267:         *faceTypes = typesTmp;
268:       }
269:       if (faceSizes) {
270:         sizesTmp[0] = 4; sizesTmp[1] = 4; sizesTmp[2] = 4; sizesTmp[3] = 4; sizesTmp[4] = 4; sizesTmp[5] = 4;
271:         *faceSizes = sizesTmp;
272:       }
273:       if (faces) {
274:         facesTmp[0]  = cone[0]; facesTmp[1]  = cone[1]; facesTmp[2]  = cone[2]; facesTmp[3]  = cone[3]; /* Bottom */
275:         facesTmp[4]  = cone[4]; facesTmp[5]  = cone[5]; facesTmp[6]  = cone[6]; facesTmp[7]  = cone[7]; /* Top */
276:         facesTmp[8]  = cone[0]; facesTmp[9]  = cone[1]; facesTmp[10] = cone[4]; facesTmp[11] = cone[5]; /* Front */
277:         facesTmp[12] = cone[1]; facesTmp[13] = cone[2]; facesTmp[14] = cone[5]; facesTmp[15] = cone[6]; /* Right */
278:         facesTmp[16] = cone[2]; facesTmp[17] = cone[3]; facesTmp[18] = cone[6]; facesTmp[19] = cone[7]; /* Back */
279:         facesTmp[20] = cone[3]; facesTmp[21] = cone[0]; facesTmp[22] = cone[7]; facesTmp[23] = cone[4]; /* Left */
280:         *faces = facesTmp;
281:       }
282:       break;
283:     case DM_POLYTOPE_PYRAMID:
284:       /*
285:        4----
286:        |\-\ \-----
287:        | \ -\     \
288:        |  1--\-----2
289:        | /    \   /
290:        |/      \ /
291:        0--------3
292:        */
293:       if (numFaces) *numFaces = 5;
294:       if (faceTypes) {
295:         typesTmp[0] = DM_POLYTOPE_QUADRILATERAL;
296:         typesTmp[1] = DM_POLYTOPE_TRIANGLE; typesTmp[2] = DM_POLYTOPE_TRIANGLE; typesTmp[3] = DM_POLYTOPE_TRIANGLE; typesTmp[4] = DM_POLYTOPE_TRIANGLE;
297:         *faceTypes = typesTmp;
298:       }
299:       if (faceSizes) {
300:         sizesTmp[0] = 4;
301:         sizesTmp[1] = 3; sizesTmp[2] = 3; sizesTmp[3] = 3; sizesTmp[4] = 3;
302:         *faceSizes = sizesTmp;
303:       }
304:       if (faces) {
305:         facesTmp[0]  = cone[0]; facesTmp[1]  = cone[1]; facesTmp[2]  = cone[2]; facesTmp[3]  = cone[3]; /* Bottom */
306:         facesTmp[4]  = cone[0]; facesTmp[5]  = cone[3]; facesTmp[6]  = cone[4];                         /* Front */
307:         facesTmp[7]  = cone[3]; facesTmp[8]  = cone[2]; facesTmp[9]  = cone[4];                         /* Right */
308:         facesTmp[10] = cone[2]; facesTmp[11] = cone[1]; facesTmp[12] = cone[4];                         /* Back */
309:         facesTmp[13] = cone[1]; facesTmp[14] = cone[0]; facesTmp[15] = cone[4];                         /* Left */
310:         *faces = facesTmp;
311:       }
312:       break;
313:     default: SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_OUTOFRANGE, "No face description for cell type %s", DMPolytopeTypes[ct]);
314:   }
315:   return(0);
316: }

318: PetscErrorCode DMPlexRestoreRawFaces_Internal(DM dm, DMPolytopeType ct, const PetscInt cone[], PetscInt *numFaces, const DMPolytopeType *faceTypes[], const PetscInt *faceSizes[], const PetscInt *faces[])
319: {
320:   PetscErrorCode  ierr;

323:   if (faceTypes) {DMRestoreWorkArray(dm, 0, MPIU_INT, (void *) faceTypes);}
324:   if (faceSizes) {DMRestoreWorkArray(dm, 0, MPIU_INT, (void *) faceSizes);}
325:   if (faces)     {DMRestoreWorkArray(dm, 0, MPIU_INT, (void *) faces);}
326:   return(0);
327: }

329: /* This interpolates faces for cells at some stratum */
330: static PetscErrorCode DMPlexInterpolateFaces_Internal(DM dm, PetscInt cellDepth, DM idm)
331: {
332:   DMLabel        ctLabel;
333:   PetscHashIJKL  faceTable;
334:   PetscInt       faceTypeNum[DM_NUM_POLYTOPES];
335:   PetscInt       depth, d, pStart, Np, cStart, cEnd, c, fStart, fEnd;

339:   DMPlexGetDepth(dm, &depth);
340:   PetscHashIJKLCreate(&faceTable);
341:   PetscArrayzero(faceTypeNum, DM_NUM_POLYTOPES);
342:   DMPlexGetDepthStratum(dm, cellDepth, &cStart, &cEnd);
343:   /* Number new faces and save face vertices in hash table */
344:   DMPlexGetDepthStratum(dm, depth > cellDepth ? cellDepth : 0, NULL, &fStart);
345:   fEnd = fStart;
346:   for (c = cStart; c < cEnd; ++c) {
347:     const PetscInt       *cone, *faceSizes, *faces;
348:     const DMPolytopeType *faceTypes;
349:     DMPolytopeType        ct;
350:     PetscInt              numFaces, cf, foff = 0;

352:     DMPlexGetCellType(dm, c, &ct);
353:     DMPlexGetCone(dm, c, &cone);
354:     DMPlexGetRawFaces_Internal(dm, ct, cone, &numFaces, &faceTypes, &faceSizes, &faces);
355:     for (cf = 0; cf < numFaces; foff += faceSizes[cf], ++cf) {
356:       const PetscInt       faceSize = faceSizes[cf];
357:       const DMPolytopeType faceType = faceTypes[cf];
358:       const PetscInt      *face     = &faces[foff];
359:       PetscHashIJKLKey     key;
360:       PetscHashIter        iter;
361:       PetscBool            missing;

363:       if (faceSize > 4) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "Do not support faces of size %D > 4", faceSize);
364:       key.i = face[0];
365:       key.j = faceSize > 1 ? face[1] : PETSC_MAX_INT;
366:       key.k = faceSize > 2 ? face[2] : PETSC_MAX_INT;
367:       key.l = faceSize > 3 ? face[3] : PETSC_MAX_INT;
368:       PetscSortInt(faceSize, (PetscInt *) &key);
369:       PetscHashIJKLPut(faceTable, key, &iter, &missing);
370:       if (missing) {
371:         PetscHashIJKLIterSet(faceTable, iter, fEnd++);
372:         ++faceTypeNum[faceType];
373:       }
374:     }
375:     DMPlexRestoreRawFaces_Internal(dm, ct, cone, &numFaces, &faceTypes, &faceSizes, &faces);
376:   }
377:   /* We need to number faces contiguously among types */
378:   {
379:     PetscInt faceTypeStart[DM_NUM_POLYTOPES], ct, numFT = 0;

381:     for (ct = 0; ct < DM_NUM_POLYTOPES; ++ct) {if (faceTypeNum[ct]) ++numFT; faceTypeStart[ct] = 0;}
382:     if (numFT > 1) {
383:       PetscHashIJKLClear(faceTable);
384:       faceTypeStart[0] = fStart;
385:       for (ct = 1; ct < DM_NUM_POLYTOPES; ++ct) faceTypeStart[ct] = faceTypeStart[ct-1] + faceTypeNum[ct-1];
386:       for (c = cStart; c < cEnd; ++c) {
387:         const PetscInt       *cone, *faceSizes, *faces;
388:         const DMPolytopeType *faceTypes;
389:         DMPolytopeType        ct;
390:         PetscInt              numFaces, cf, foff = 0;

392:         DMPlexGetCellType(dm, c, &ct);
393:         DMPlexGetCone(dm, c, &cone);
394:         DMPlexGetRawFaces_Internal(dm, ct, cone, &numFaces, &faceTypes, &faceSizes, &faces);
395:         for (cf = 0; cf < numFaces; foff += faceSizes[cf], ++cf) {
396:           const PetscInt       faceSize = faceSizes[cf];
397:           const DMPolytopeType faceType = faceTypes[cf];
398:           const PetscInt      *face     = &faces[foff];
399:           PetscHashIJKLKey     key;
400:           PetscHashIter        iter;
401:           PetscBool            missing;

403:           if (faceSize > 4) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "Do not support faces of size %D > 4", faceSize);
404:           key.i = face[0];
405:           key.j = faceSize > 1 ? face[1] : PETSC_MAX_INT;
406:           key.k = faceSize > 2 ? face[2] : PETSC_MAX_INT;
407:           key.l = faceSize > 3 ? face[3] : PETSC_MAX_INT;
408:           PetscSortInt(faceSize, (PetscInt *) &key);
409:           PetscHashIJKLPut(faceTable, key, &iter, &missing);
410:           if (missing) {PetscHashIJKLIterSet(faceTable, iter, faceTypeStart[faceType]++);}
411:         }
412:         DMPlexRestoreRawFaces_Internal(dm, ct, cone, &numFaces, &faceTypes, &faceSizes, &faces);
413:       }
414:       for (ct = 1; ct < DM_NUM_POLYTOPES; ++ct) {
415:         if (faceTypeStart[ct] != faceTypeStart[ct-1] + faceTypeNum[ct]) SETERRQ4(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Inconsistent numbering for cell type %s, %D != %D + %D", DMPolytopeTypes[ct], faceTypeStart[ct], faceTypeStart[ct-1], faceTypeNum[ct]);
416:       }
417:     }
418:   }
419:   /* Add new points, always at the end of the numbering */
420:   DMPlexGetChart(dm, &pStart, &Np);
421:   DMPlexSetChart(idm, pStart, Np + (fEnd - fStart));
422:   /* Set cone sizes */
423:   /*   Must create the celltype label here so that we do not automatically try to compute the types */
424:   DMCreateLabel(idm, "celltype");
425:   DMPlexGetCellTypeLabel(idm, &ctLabel);
426:   for (d = 0; d <= depth; ++d) {
427:     DMPolytopeType ct;
428:     PetscInt       coneSize, pStart, pEnd, p;

430:     if (d == cellDepth) continue;
431:     DMPlexGetDepthStratum(dm, d, &pStart, &pEnd);
432:     for (p = pStart; p < pEnd; ++p) {
433:       DMPlexGetConeSize(dm, p, &coneSize);
434:       DMPlexSetConeSize(idm, p, coneSize);
435:       DMPlexGetCellType(dm, p, &ct);
436:       DMPlexSetCellType(idm, p, ct);
437:     }
438:   }
439:   for (c = cStart; c < cEnd; ++c) {
440:     const PetscInt       *cone, *faceSizes, *faces;
441:     const DMPolytopeType *faceTypes;
442:     DMPolytopeType        ct;
443:     PetscInt              numFaces, cf, foff = 0;

445:     DMPlexGetCellType(dm, c, &ct);
446:     DMPlexGetCone(dm, c, &cone);
447:     DMPlexGetRawFaces_Internal(dm, ct, cone, &numFaces, &faceTypes, &faceSizes, &faces);
448:     DMPlexSetCellType(idm, c, ct);
449:     DMPlexSetConeSize(idm, c, numFaces);
450:     for (cf = 0; cf < numFaces; foff += faceSizes[cf], ++cf) {
451:       const PetscInt       faceSize = faceSizes[cf];
452:       const DMPolytopeType faceType = faceTypes[cf];
453:       const PetscInt      *face     = &faces[foff];
454:       PetscHashIJKLKey     key;
455:       PetscHashIter        iter;
456:       PetscBool            missing;
457:       PetscInt             f;

459:       if (faceSize > 4) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "Do not support faces of size %D > 4", faceSize);
460:       key.i = face[0];
461:       key.j = faceSize > 1 ? face[1] : PETSC_MAX_INT;
462:       key.k = faceSize > 2 ? face[2] : PETSC_MAX_INT;
463:       key.l = faceSize > 3 ? face[3] : PETSC_MAX_INT;
464:       PetscSortInt(faceSize, (PetscInt *) &key);
465:       PetscHashIJKLPut(faceTable, key, &iter, &missing);
466:       if (missing) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Missing face (cell %D, lf %D)", c, cf);
467:       PetscHashIJKLIterGet(faceTable, iter, &f);
468:       DMPlexSetConeSize(idm, f, faceSize);
469:       DMPlexSetCellType(idm, f, faceType);
470:     }
471:     DMPlexRestoreRawFaces_Internal(dm, ct, cone, &numFaces, &faceTypes, &faceSizes, &faces);
472:   }
473:   DMSetUp(idm);
474:   /* Initialize cones so we do not need the bash table to tell us that a cone has been set */
475:   {
476:     PetscSection cs;
477:     PetscInt    *cones, csize;

479:     DMPlexGetConeSection(idm, &cs);
480:     DMPlexGetCones(idm, &cones);
481:     PetscSectionGetStorageSize(cs, &csize);
482:     for (c = 0; c < csize; ++c) cones[c] = -1;
483:   }
484:   /* Set cones */
485:   for (d = 0; d <= depth; ++d) {
486:     const PetscInt *cone;
487:     PetscInt        pStart, pEnd, p;

489:     if (d == cellDepth) continue;
490:     DMPlexGetDepthStratum(dm, d, &pStart, &pEnd);
491:     for (p = pStart; p < pEnd; ++p) {
492:       DMPlexGetCone(dm, p, &cone);
493:       DMPlexSetCone(idm, p, cone);
494:       DMPlexGetConeOrientation(dm, p, &cone);
495:       DMPlexSetConeOrientation(idm, p, cone);
496:     }
497:   }
498:   for (c = cStart; c < cEnd; ++c) {
499:     const PetscInt       *cone, *faceSizes, *faces;
500:     const DMPolytopeType *faceTypes;
501:     DMPolytopeType        ct;
502:     PetscInt              numFaces, cf, foff = 0;

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

517:       if (faceSize > 4) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "Do not support faces of size %D > 4", faceSize);
518:       key.i = face[0];
519:       key.j = faceSize > 1 ? face[1] : PETSC_MAX_INT;
520:       key.k = faceSize > 2 ? face[2] : PETSC_MAX_INT;
521:       key.l = faceSize > 3 ? face[3] : PETSC_MAX_INT;
522:       PetscSortInt(faceSize, (PetscInt *) &key);
523:       PetscHashIJKLPut(faceTable, key, &iter, &missing);
524:       PetscHashIJKLIterGet(faceTable, iter, &f);
525:       DMPlexInsertCone(idm, c, cf, f);
526:       DMPlexGetCone(idm, f, &fcone);
527:       if (fcone[0] < 0) {DMPlexSetCone(idm, f, face);}
528:       {
529:         const PetscInt *cone;
530:         PetscInt        coneSize, ornt;

532:         DMPlexGetConeSize(idm, f, &coneSize);
533:         DMPlexGetCone(idm, f, &cone);
534:         if (coneSize != faceSize) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid number of face vertices %D for face %D should be %D", coneSize, f, faceSize);
535:         /* Notice that we have to use vertices here because the lower dimensional faces have not been created yet */
536:         DMPolytopeGetVertexOrientation(faceType, cone, face, &ornt);
537:         DMPlexInsertConeOrientation(idm, c, cf, ornt);
538:       }
539:     }
540:     DMPlexRestoreRawFaces_Internal(dm, ct, cone, &numFaces, &faceTypes, &faceSizes, &faces);
541:   }
542:   PetscHashIJKLDestroy(&faceTable);
543:   DMPlexSymmetrize(idm);
544:   DMPlexStratify(idm);
545:   return(0);
546: }

548: static PetscErrorCode SortRmineRremoteByRemote_Private(PetscSF sf, PetscInt *rmine1[], PetscInt *rremote1[])
549: {
550:   PetscInt            nleaves;
551:   PetscInt            nranks;
552:   const PetscMPIInt  *ranks=NULL;
553:   const PetscInt     *roffset=NULL, *rmine=NULL, *rremote=NULL;
554:   PetscInt            n, o, r;
555:   PetscErrorCode      ierr;

558:   PetscSFGetRootRanks(sf, &nranks, &ranks, &roffset, &rmine, &rremote);
559:   nleaves = roffset[nranks];
560:   PetscMalloc2(nleaves, rmine1, nleaves, rremote1);
561:   for (r=0; r<nranks; r++) {
562:     /* simultaneously sort rank-wise portions of rmine & rremote by values in rremote
563:        - to unify order with the other side */
564:     o = roffset[r];
565:     n = roffset[r+1] - o;
566:     PetscArraycpy(&(*rmine1)[o], &rmine[o], n);
567:     PetscArraycpy(&(*rremote1)[o], &rremote[o], n);
568:     PetscSortIntWithArray(n, &(*rremote1)[o], &(*rmine1)[o]);
569:   }
570:   return(0);
571: }

573: PetscErrorCode DMPlexOrientInterface_Internal(DM dm)
574: {
575:   PetscSF            sf;
576:   const PetscInt    *locals;
577:   const PetscSFNode *remotes;
578:   const PetscMPIInt *ranks;
579:   const PetscInt    *roffset;
580:   PetscInt          *rmine1, *rremote1; /* rmine and rremote copies simultaneously sorted by rank and rremote */
581:   PetscInt           nroots, p, nleaves, nranks, r, maxConeSize = 0;
582:   PetscInt         (*roots)[4],      (*leaves)[4], mainCone[4];
583:   PetscMPIInt      (*rootsRanks)[4], (*leavesRanks)[4];
584:   MPI_Comm           comm;
585:   PetscMPIInt        rank, size;
586:   PetscInt           debug = 0;
587:   PetscErrorCode     ierr;

590:   PetscObjectGetComm((PetscObject) dm, &comm);
591:   MPI_Comm_rank(comm, &rank);
592:   MPI_Comm_size(comm, &size);
593:   DMGetPointSF(dm, &sf);
594:   DMViewFromOptions(dm, NULL, "-before_orient_interface_dm_view");
595:   if (PetscDefined(USE_DEBUG)) {DMPlexCheckPointSF(dm);}
596:   PetscSFGetGraph(sf, &nroots, &nleaves, &locals, &remotes);
597:   if (nroots < 0) return(0);
598:   PetscSFSetUp(sf);
599:   SortRmineRremoteByRemote_Private(sf, &rmine1, &rremote1);
600:   for (p = 0; p < nleaves; ++p) {
601:     PetscInt coneSize;
602:     DMPlexGetConeSize(dm, locals[p], &coneSize);
603:     maxConeSize = PetscMax(maxConeSize, coneSize);
604:   }
605:   if (maxConeSize > 4) SETERRQ1(comm, PETSC_ERR_SUP, "This method does not support cones of size %D", maxConeSize);
606:   PetscMalloc4(nroots, &roots, nroots, &leaves, nroots, &rootsRanks, nroots, &leavesRanks);
607:   for (p = 0; p < nroots; ++p) {
608:     const PetscInt *cone;
609:     PetscInt        coneSize, c, ind0;

611:     DMPlexGetConeSize(dm, p, &coneSize);
612:     DMPlexGetCone(dm, p, &cone);
613:     /* Ignore vertices */
614:     if (coneSize < 2) {
615:       for (c = 0; c < 4; c++) {
616:         roots[p][c]      = -1;
617:         rootsRanks[p][c] = -1;
618:       }
619:       continue;
620:     }
621:     /* Translate all points to root numbering */
622:     for (c = 0; c < PetscMin(coneSize, 4); c++) {
623:       PetscFindInt(cone[c], nleaves, locals, &ind0);
624:       if (ind0 < 0) {
625:         roots[p][c]      = cone[c];
626:         rootsRanks[p][c] = rank;
627:       } else {
628:         roots[p][c]      = remotes[ind0].index;
629:         rootsRanks[p][c] = remotes[ind0].rank;
630:       }
631:     }
632:     for (c = coneSize; c < 4; c++) {
633:       roots[p][c]      = -1;
634:       rootsRanks[p][c] = -1;
635:     }
636:   }
637:   for (p = 0; p < nroots; ++p) {
638:     PetscInt c;
639:     for (c = 0; c < 4; c++) {
640:       leaves[p][c]      = -2;
641:       leavesRanks[p][c] = -2;
642:     }
643:   }
644:   PetscSFBcastBegin(sf, MPIU_4INT, roots, leaves, MPI_REPLACE);
645:   PetscSFBcastBegin(sf, MPI_4INT, rootsRanks, leavesRanks, MPI_REPLACE);
646:   PetscSFBcastEnd(sf, MPIU_4INT, roots, leaves, MPI_REPLACE);
647:   PetscSFBcastEnd(sf, MPI_4INT, rootsRanks, leavesRanks, MPI_REPLACE);
648:   if (debug) {
649:     PetscSynchronizedFlush(comm, NULL);
650:     if (!rank) {PetscSynchronizedPrintf(comm, "Referenced roots\n");}
651:   }
652:   PetscSFGetRootRanks(sf, &nranks, &ranks, &roffset, NULL, NULL);
653:   for (p = 0; p < nroots; ++p) {
654:     DMPolytopeType  ct;
655:     const PetscInt *cone;
656:     PetscInt        coneSize, c, ind0, o;

658:     if (leaves[p][0] < 0) continue; /* Ignore vertices */
659:     DMPlexGetCellType(dm, p, &ct);
660:     DMPlexGetConeSize(dm, p, &coneSize);
661:     DMPlexGetCone(dm, p, &cone);
662:     if (debug) {
663:       PetscSynchronizedPrintf(comm, "[%d]  %4D: cone=[%4D %4D %4D %4D] roots=[(%d,%4D) (%d,%4D) (%d,%4D) (%d,%4D)] leaves=[(%d,%4D) (%d,%4D) (%d,%4D) (%d,%4D)]",
664:        rank, p, cone[0], cone[1], cone[2], cone[3],
665:        rootsRanks[p][0], roots[p][0], rootsRanks[p][1], roots[p][1], rootsRanks[p][2], roots[p][2], rootsRanks[p][3], roots[p][3],
666:        leavesRanks[p][0], leaves[p][0], leavesRanks[p][1], leaves[p][1], leavesRanks[p][2], leaves[p][2], leavesRanks[p][3], leaves[p][3]);
667:     }
668:     if (leavesRanks[p][0] != rootsRanks[p][0] || leaves[p][0] != roots[p][0] ||
669:         leavesRanks[p][1] != rootsRanks[p][1] || leaves[p][1] != roots[p][1] ||
670:         leavesRanks[p][2] != rootsRanks[p][2] || leaves[p][2] != roots[p][2] ||
671:         leavesRanks[p][3] != rootsRanks[p][3] || leaves[p][3] != roots[p][3]) {
672:       /* Translate these leaves to my cone points; mainCone means desired order p's cone points */
673:       for (c = 0; c < PetscMin(coneSize, 4); ++c) {
674:         PetscInt rS, rN;

676:         if (leavesRanks[p][c] == rank) {
677:           /* A local leaf is just taken as it is */
678:           mainCone[c] = leaves[p][c];
679:           continue;
680:         }
681:         /* Find index of rank leavesRanks[p][c] among remote ranks */
682:         /* No need for PetscMPIIntCast because these integers were originally cast from PetscMPIInt. */
683:         PetscFindMPIInt((PetscMPIInt) leavesRanks[p][c], nranks, ranks, &r);
684:         if (PetscUnlikely(r < 0)) SETERRQ7(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Point %D cone[%D]=%D root (%d,%D) leaf (%d,%D): leaf rank not found among remote ranks", p, c, cone[c], rootsRanks[p][c], roots[p][c], leavesRanks[p][c], leaves[p][c]);
685:         if (PetscUnlikely(ranks[r] < 0 || ranks[r] >= size)) SETERRQ5(PETSC_COMM_SELF, PETSC_ERR_PLIB, "p=%D c=%D commsize=%d: ranks[%D] = %d makes no sense", p, c, size, r, ranks[r]);
686:         /* Find point leaves[p][c] among remote points aimed at rank leavesRanks[p][c] */
687:         rS = roffset[r];
688:         rN = roffset[r+1] - rS;
689:         PetscFindInt(leaves[p][c], rN, &rremote1[rS], &ind0);
690:         if (PetscUnlikely(ind0 < 0)) SETERRQ7(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Point %D cone[%D]=%D root (%d,%D) leave (%d,%D): 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]);
691:         /* Get the corresponding local point */
692:         mainCone[c] = rmine1[rS + ind0];
693:       }
694:       if (debug) {PetscSynchronizedPrintf(comm, " mainCone=[%4D %4D %4D %4D]\n", mainCone[0], mainCone[1], mainCone[2], mainCone[3]);}
695:       /* Set the desired order of p's cone points and fix orientations accordingly */
696:       DMPolytopeGetOrientation(ct, cone, mainCone, &o);
697:       DMPlexOrientPoint(dm, p, o);
698:     } else if (debug) {PetscSynchronizedPrintf(comm, " ==\n");}
699:   }
700:   if (debug) {
701:     PetscSynchronizedFlush(comm, NULL);
702:     MPI_Barrier(comm);
703:   }
704:   DMViewFromOptions(dm, NULL, "-after_orient_interface_dm_view");
705:   PetscFree4(roots, leaves, rootsRanks, leavesRanks);
706:   PetscFree2(rmine1, rremote1);
707:   return(0);
708: }

710: static PetscErrorCode IntArrayViewFromOptions(MPI_Comm comm, const char opt[], const char name[], const char idxname[], const char valname[], PetscInt n, const PetscInt a[])
711: {
712:   PetscInt       idx;
713:   PetscMPIInt    rank;
714:   PetscBool      flg;

718:   PetscOptionsHasName(NULL, NULL, opt, &flg);
719:   if (!flg) return(0);
720:   MPI_Comm_rank(comm, &rank);
721:   PetscSynchronizedPrintf(comm, "[%d]%s:\n", rank, name);
722:   for (idx = 0; idx < n; ++idx) {PetscSynchronizedPrintf(comm, "[%d]%s %D %s %D\n", rank, idxname, idx, valname, a[idx]);}
723:   PetscSynchronizedFlush(comm, NULL);
724:   return(0);
725: }

727: static PetscErrorCode SFNodeArrayViewFromOptions(MPI_Comm comm, const char opt[], const char name[], const char idxname[], PetscInt n, const PetscSFNode a[])
728: {
729:   PetscInt       idx;
730:   PetscMPIInt    rank;
731:   PetscBool      flg;

735:   PetscOptionsHasName(NULL, NULL, opt, &flg);
736:   if (!flg) return(0);
737:   MPI_Comm_rank(comm, &rank);
738:   PetscSynchronizedPrintf(comm, "[%d]%s:\n", rank, name);
739:   if (idxname) {
740:     for (idx = 0; idx < n; ++idx) {PetscSynchronizedPrintf(comm, "[%d]%s %D rank %D index %D\n", rank, idxname, idx, a[idx].rank, a[idx].index);}
741:   } else {
742:     for (idx = 0; idx < n; ++idx) {PetscSynchronizedPrintf(comm, "[%d]rank %D index %D\n", rank, a[idx].rank, a[idx].index);}
743:   }
744:   PetscSynchronizedFlush(comm, NULL);
745:   return(0);
746: }

748: static PetscErrorCode DMPlexMapToLocalPoint(DM dm, PetscHMapIJ remotehash, PetscSFNode remotePoint, PetscInt *localPoint)
749: {
750:   PetscSF         sf;
751:   const PetscInt *locals;
752:   PetscMPIInt     rank;
753:   PetscErrorCode  ierr;

756:   MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank);
757:   DMGetPointSF(dm, &sf);
758:   PetscSFGetGraph(sf, NULL, NULL, &locals, NULL);
759:   if (remotePoint.rank == rank) {
760:     *localPoint = remotePoint.index;
761:   } else {
762:     PetscHashIJKey key;
763:     PetscInt       l;

765:     key.i = remotePoint.index;
766:     key.j = remotePoint.rank;
767:     PetscHMapIJGet(remotehash, key, &l);
768:     if (l >= 0) {
769:       *localPoint = locals[l];
770:     } else PetscFunctionReturn(1);
771:   }
772:   return(0);
773: }

775: static PetscErrorCode DMPlexMapToGlobalPoint(DM dm, PetscInt localPoint, PetscSFNode *remotePoint)
776: {
777:   PetscSF            sf;
778:   const PetscInt    *locals, *rootdegree;
779:   const PetscSFNode *remotes;
780:   PetscInt           Nl, l;
781:   PetscMPIInt        rank;
782:   PetscErrorCode     ierr;

785:   MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank);
786:   DMGetPointSF(dm, &sf);
787:   PetscSFGetGraph(sf, NULL, &Nl, &locals, &remotes);
788:   if (Nl < 0) goto owned;
789:   PetscSFComputeDegreeBegin(sf, &rootdegree);
790:   PetscSFComputeDegreeEnd(sf, &rootdegree);
791:   if (rootdegree[localPoint]) goto owned;
792:   PetscFindInt(localPoint, Nl, locals, &l);
793:   if (l < 0) PetscFunctionReturn(1);
794:   *remotePoint = remotes[l];
795:   return(0);
796:   owned:
797:   remotePoint->rank  = rank;
798:   remotePoint->index = localPoint;
799:   return(0);
800: }

802: static PetscErrorCode DMPlexPointIsShared(DM dm, PetscInt p, PetscBool *isShared)
803: {
804:   PetscSF         sf;
805:   const PetscInt *locals, *rootdegree;
806:   PetscInt        Nl, idx;
807:   PetscErrorCode  ierr;

810:   *isShared = PETSC_FALSE;
811:   DMGetPointSF(dm, &sf);
812:   PetscSFGetGraph(sf, NULL, &Nl, &locals, NULL);
813:   if (Nl < 0) return(0);
814:   PetscFindInt(p, Nl, locals, &idx);
815:   if (idx >= 0) {*isShared = PETSC_TRUE; return(0);}
816:   PetscSFComputeDegreeBegin(sf, &rootdegree);
817:   PetscSFComputeDegreeEnd(sf, &rootdegree);
818:   if (rootdegree[p] > 0) *isShared = PETSC_TRUE;
819:   return(0);
820: }

822: static PetscErrorCode DMPlexConeIsShared(DM dm, PetscInt p, PetscBool *isShared)
823: {
824:   const PetscInt *cone;
825:   PetscInt        coneSize, c;
826:   PetscBool       cShared = PETSC_TRUE;
827:   PetscErrorCode  ierr;

830:   DMPlexGetConeSize(dm, p, &coneSize);
831:   DMPlexGetCone(dm, p, &cone);
832:   for (c = 0; c < coneSize; ++c) {
833:     PetscBool pointShared;

835:     DMPlexPointIsShared(dm, cone[c], &pointShared);
836:     cShared = (PetscBool) (cShared && pointShared);
837:   }
838:   *isShared = coneSize ? cShared : PETSC_FALSE;
839:   return(0);
840: }

842: static PetscErrorCode DMPlexGetConeMinimum(DM dm, PetscInt p, PetscSFNode *cpmin)
843: {
844:   const PetscInt *cone;
845:   PetscInt        coneSize, c;
846:   PetscSFNode     cmin = {PETSC_MAX_INT, PETSC_MAX_INT}, missing = {-1, -1};
847:   PetscErrorCode  ierr;

850:   DMPlexGetConeSize(dm, p, &coneSize);
851:   DMPlexGetCone(dm, p, &cone);
852:   for (c = 0; c < coneSize; ++c) {
853:     PetscSFNode rcp;

855:     DMPlexMapToGlobalPoint(dm, cone[c], &rcp);
856:     if (ierr) {
857:       cmin = missing;
858:     } else {
859:       cmin = (rcp.rank < cmin.rank) || (rcp.rank == cmin.rank && rcp.index < cmin.index) ? rcp : cmin;
860:     }
861:   }
862:   *cpmin = coneSize ? cmin : missing;
863:   return(0);
864: }

866: /*
867:   Each shared face has an entry in the candidates array:
868:     (-1, coneSize-1), {(global cone point)}
869:   where the set is missing the point p which we use as the key for the face
870: */
871: static PetscErrorCode DMPlexAddSharedFace_Private(DM dm, PetscSection candidateSection, PetscSFNode candidates[], PetscHMapIJ faceHash, PetscInt p, PetscBool debug)
872: {
873:   MPI_Comm        comm;
874:   const PetscInt *support;
875:   PetscInt        supportSize, s, off = 0, idx = 0, overlap, cellHeight, height;
876:   PetscMPIInt     rank;
877:   PetscErrorCode  ierr;

880:   PetscObjectGetComm((PetscObject) dm, &comm);
881:   MPI_Comm_rank(comm, &rank);
882:   DMPlexGetOverlap(dm, &overlap);
883:   DMPlexGetVTKCellHeight(dm, &cellHeight);
884:   DMPlexGetPointHeight(dm, p, &height);
885:   if (!overlap && height <= cellHeight+1) {
886:     /* cells can't be shared for non-overlapping meshes */
887:     if (debug) {PetscSynchronizedPrintf(comm, "[%d]    Skipping face %D to avoid adding cell to hashmap since this is nonoverlapping mesh\n", rank, p);}
888:     return(0);
889:   }
890:   DMPlexGetSupportSize(dm, p, &supportSize);
891:   DMPlexGetSupport(dm, p, &support);
892:   if (candidates) {PetscSectionGetOffset(candidateSection, p, &off);}
893:   for (s = 0; s < supportSize; ++s) {
894:     const PetscInt  face = support[s];
895:     const PetscInt *cone;
896:     PetscSFNode     cpmin={-1,-1}, rp={-1,-1};
897:     PetscInt        coneSize, c, f;
898:     PetscBool       isShared = PETSC_FALSE;
899:     PetscHashIJKey  key;

901:     /* Only add point once */
902:     if (debug) {PetscSynchronizedPrintf(comm, "[%d]    Support face %D\n", rank, face);}
903:     key.i = p;
904:     key.j = face;
905:     PetscHMapIJGet(faceHash, key, &f);
906:     if (f >= 0) continue;
907:     DMPlexConeIsShared(dm, face, &isShared);
908:     DMPlexGetConeMinimum(dm, face, &cpmin);
909:     DMPlexMapToGlobalPoint(dm, p, &rp);
910:     if (debug) {
911:       PetscSynchronizedPrintf(comm, "[%d]      Face point %D is shared: %d\n", rank, face, (int) isShared);
912:       PetscSynchronizedPrintf(comm, "[%d]      Global point (%D, %D) Min Cone Point (%D, %D)\n", rank, rp.rank, rp.index, cpmin.rank, cpmin.index);
913:     }
914:     if (isShared && (rp.rank == cpmin.rank && rp.index == cpmin.index)) {
915:       PetscHMapIJSet(faceHash, key, p);
916:       if (candidates) {
917:         if (debug) {PetscSynchronizedPrintf(comm, "[%d]    Adding shared face %D at idx %D\n[%d]     ", rank, face, idx, rank);}
918:         DMPlexGetConeSize(dm, face, &coneSize);
919:         DMPlexGetCone(dm, face, &cone);
920:         candidates[off+idx].rank    = -1;
921:         candidates[off+idx++].index = coneSize-1;
922:         candidates[off+idx].rank    = rank;
923:         candidates[off+idx++].index = face;
924:         for (c = 0; c < coneSize; ++c) {
925:           const PetscInt cp = cone[c];

927:           if (cp == p) continue;
928:           DMPlexMapToGlobalPoint(dm, cp, &candidates[off+idx]);
929:           if (debug) {PetscSynchronizedPrintf(comm, " (%D,%D)", candidates[off+idx].rank, candidates[off+idx].index);}
930:           ++idx;
931:         }
932:         if (debug) {PetscSynchronizedPrintf(comm, "\n");}
933:       } else {
934:         /* Add cone size to section */
935:         if (debug) {PetscSynchronizedPrintf(comm, "[%d]    Scheduling shared face %D\n", rank, face);}
936:         DMPlexGetConeSize(dm, face, &coneSize);
937:         PetscHMapIJSet(faceHash, key, p);
938:         PetscSectionAddDof(candidateSection, p, coneSize+1);
939:       }
940:     }
941:   }
942:   return(0);
943: }

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

948:   Collective on dm

950:   Input Parameters:
951: + dm      - The interpolated DM
952: - pointSF - The initial SF without interpolated points

954:   Output Parameter:
955: . pointSF - The SF including interpolated points

957:   Level: developer

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

961: .seealso: DMPlexInterpolate(), DMPlexUninterpolate()
962: @*/
963: PetscErrorCode DMPlexInterpolatePointSF(DM dm, PetscSF pointSF)
964: {
965:   MPI_Comm           comm;
966:   PetscHMapIJ        remoteHash;
967:   PetscHMapI         claimshash;
968:   PetscSection       candidateSection, candidateRemoteSection, claimSection;
969:   PetscSFNode       *candidates, *candidatesRemote, *claims;
970:   const PetscInt    *localPoints, *rootdegree;
971:   const PetscSFNode *remotePoints;
972:   PetscInt           ov, Nr, r, Nl, l;
973:   PetscInt           candidatesSize, candidatesRemoteSize, claimsSize;
974:   PetscBool          flg, debug = PETSC_FALSE;
975:   PetscMPIInt        rank;
976:   PetscErrorCode     ierr;

981:   DMPlexIsDistributed(dm, &flg);
982:   if (!flg) return(0);
983:   /* Set initial SF so that lower level queries work */
984:   DMSetPointSF(dm, pointSF);
985:   PetscObjectGetComm((PetscObject) dm, &comm);
986:   MPI_Comm_rank(comm, &rank);
987:   DMPlexGetOverlap(dm, &ov);
988:   if (ov) SETERRQ(comm, PETSC_ERR_SUP, "Interpolation of overlapped DMPlex not implemented yet");
989:   PetscOptionsHasName(NULL, ((PetscObject) dm)->prefix, "-dmplex_interp_debug", &debug);
990:   PetscObjectViewFromOptions((PetscObject) dm, NULL, "-dm_interp_pre_view");
991:   PetscObjectViewFromOptions((PetscObject) pointSF, NULL, "-petscsf_interp_pre_view");
992:   PetscLogEventBegin(DMPLEX_InterpolateSF,dm,0,0,0);
993:   /* Step 0: Precalculations */
994:   PetscSFGetGraph(pointSF, &Nr, &Nl, &localPoints, &remotePoints);
995:   if (Nr < 0) SETERRQ(comm, PETSC_ERR_ARG_WRONGSTATE, "This DMPlex is distributed but input PointSF has no graph set");
996:   PetscHMapIJCreate(&remoteHash);
997:   for (l = 0; l < Nl; ++l) {
998:     PetscHashIJKey key;
999:     key.i = remotePoints[l].index;
1000:     key.j = remotePoints[l].rank;
1001:     PetscHMapIJSet(remoteHash, key, l);
1002:   }
1003:   /*   Compute root degree to identify shared points */
1004:   PetscSFComputeDegreeBegin(pointSF, &rootdegree);
1005:   PetscSFComputeDegreeEnd(pointSF, &rootdegree);
1006:   IntArrayViewFromOptions(comm, "-interp_root_degree_view", "Root degree", "point", "degree", Nr, rootdegree);
1007:   /*
1008:   1) Loop over each leaf point $p$ at depth $d$ in the SF
1009:   \item Get set $F(p)$ of faces $f$ in the support of $p$ for which
1010:   \begin{itemize}
1011:     \item all cone points of $f$ are shared
1012:     \item $p$ is the cone point with smallest canonical number
1013:   \end{itemize}
1014:   \item Send $F(p)$ and the cone of each face to the active root point $r(p)$
1015:   \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
1016:   \item Send the root face from the root back to all leaf process
1017:   \item Leaf processes add the shared face to the SF
1018:   */
1019:   /* Step 1: Construct section+SFNode array
1020:        The section has entries for all shared faces for which we have a leaf point in the cone
1021:        The array holds candidate shared faces, each face is refered to by the leaf point */
1022:   PetscSectionCreate(comm, &candidateSection);
1023:   PetscSectionSetChart(candidateSection, 0, Nr);
1024:   {
1025:     PetscHMapIJ faceHash;

1027:     PetscHMapIJCreate(&faceHash);
1028:     for (l = 0; l < Nl; ++l) {
1029:       const PetscInt p = localPoints[l];

1031:       if (debug) {PetscSynchronizedPrintf(comm, "[%d]  First pass leaf point %D\n", rank, p);}
1032:       DMPlexAddSharedFace_Private(dm, candidateSection, NULL, faceHash, p, debug);
1033:     }
1034:     PetscHMapIJClear(faceHash);
1035:     PetscSectionSetUp(candidateSection);
1036:     PetscSectionGetStorageSize(candidateSection, &candidatesSize);
1037:     PetscMalloc1(candidatesSize, &candidates);
1038:     for (l = 0; l < Nl; ++l) {
1039:       const PetscInt p = localPoints[l];

1041:       if (debug) {PetscSynchronizedPrintf(comm, "[%d]  Second pass leaf point %D\n", rank, p);}
1042:       DMPlexAddSharedFace_Private(dm, candidateSection, candidates, faceHash, p, debug);
1043:     }
1044:     PetscHMapIJDestroy(&faceHash);
1045:     if (debug) {PetscSynchronizedFlush(comm, NULL);}
1046:   }
1047:   PetscObjectSetName((PetscObject) candidateSection, "Candidate Section");
1048:   PetscObjectViewFromOptions((PetscObject) candidateSection, NULL, "-petscsection_interp_candidate_view");
1049:   SFNodeArrayViewFromOptions(comm, "-petscsection_interp_candidate_view", "Candidates", NULL, candidatesSize, candidates);
1050:   /* Step 2: Gather candidate section / array pair into the root partition via inverse(multi(pointSF)). */
1051:   /*   Note that this section is indexed by offsets into leaves, not by point number */
1052:   {
1053:     PetscSF   sfMulti, sfInverse, sfCandidates;
1054:     PetscInt *remoteOffsets;

1056:     PetscSFGetMultiSF(pointSF, &sfMulti);
1057:     PetscSFCreateInverseSF(sfMulti, &sfInverse);
1058:     PetscSectionCreate(comm, &candidateRemoteSection);
1059:     PetscSFDistributeSection(sfInverse, candidateSection, &remoteOffsets, candidateRemoteSection);
1060:     PetscSFCreateSectionSF(sfInverse, candidateSection, remoteOffsets, candidateRemoteSection, &sfCandidates);
1061:     PetscSectionGetStorageSize(candidateRemoteSection, &candidatesRemoteSize);
1062:     PetscMalloc1(candidatesRemoteSize, &candidatesRemote);
1063:     PetscSFBcastBegin(sfCandidates, MPIU_2INT, candidates, candidatesRemote,MPI_REPLACE);
1064:     PetscSFBcastEnd(sfCandidates, MPIU_2INT, candidates, candidatesRemote,MPI_REPLACE);
1065:     PetscSFDestroy(&sfInverse);
1066:     PetscSFDestroy(&sfCandidates);
1067:     PetscFree(remoteOffsets);

1069:     PetscObjectSetName((PetscObject) candidateRemoteSection, "Remote Candidate Section");
1070:     PetscObjectViewFromOptions((PetscObject) candidateRemoteSection, NULL, "-petscsection_interp_candidate_remote_view");
1071:     SFNodeArrayViewFromOptions(comm, "-petscsection_interp_candidate_remote_view", "Remote Candidates", NULL, candidatesRemoteSize, candidatesRemote);
1072:   }
1073:   /* 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 */
1074:   {
1075:     PetscHashIJKLRemote faceTable;
1076:     PetscInt            idx, idx2;

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

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

1086:         PetscSectionGetDof(candidateRemoteSection, idx, &dof);
1087:         PetscSectionGetOffset(candidateRemoteSection, idx, &offset);
1088:         /* dof may include many faces from the remote process */
1089:         for (d = 0; d < dof; ++d) {
1090:           const PetscInt         hidx  = offset+d;
1091:           const PetscInt         Np    = candidatesRemote[hidx].index+1;
1092:           const PetscSFNode      rface = candidatesRemote[hidx+1];
1093:           const PetscSFNode     *fcone = &candidatesRemote[hidx+2];
1094:           PetscSFNode            fcp0;
1095:           const PetscSFNode      pmax  = {PETSC_MAX_INT, PETSC_MAX_INT};
1096:           const PetscInt        *join  = NULL;
1097:           PetscHashIJKLRemoteKey key;
1098:           PetscHashIter          iter;
1099:           PetscBool              missing;
1100:           PetscInt               points[1024], p, joinSize;

1102:           if (debug) {PetscSynchronizedPrintf(PetscObjectComm((PetscObject) dm), "[%d]  Checking face (%D, %D) at (%D, %D, %D) with cone size %D\n", rank, rface.rank, rface.index, r, idx, d, Np);}
1103:           if (Np > 4) SETERRQ6(PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot handle face (%D, %D) at (%D, %D, %D) with %D cone points", rface.rank, rface.index, r, idx, d, Np);
1104:           fcp0.rank  = rank;
1105:           fcp0.index = r;
1106:           d += Np;
1107:           /* Put remote face in hash table */
1108:           key.i = fcp0;
1109:           key.j = fcone[0];
1110:           key.k = Np > 2 ? fcone[1] : pmax;
1111:           key.l = Np > 3 ? fcone[2] : pmax;
1112:           PetscSortSFNode(Np, (PetscSFNode *) &key);
1113:           PetscHashIJKLRemotePut(faceTable, key, &iter, &missing);
1114:           if (missing) {
1115:             if (debug) {PetscSynchronizedPrintf(PetscObjectComm((PetscObject) dm), "[%d]  Setting remote face (%D, %D)\n", rank, rface.index, rface.rank);}
1116:             PetscHashIJKLRemoteIterSet(faceTable, iter, rface);
1117:           } else {
1118:             PetscSFNode oface;

1120:             PetscHashIJKLRemoteIterGet(faceTable, iter, &oface);
1121:             if ((rface.rank < oface.rank) || (rface.rank == oface.rank && rface.index < oface.index)) {
1122:               if (debug) {PetscSynchronizedPrintf(PetscObjectComm((PetscObject) dm), "[%d]  Replacing with remote face (%D, %D)\n", rank, rface.index, rface.rank);}
1123:               PetscHashIJKLRemoteIterSet(faceTable, iter, rface);
1124:             }
1125:           }
1126:           /* Check for local face */
1127:           points[0] = r;
1128:           for (p = 1; p < Np; ++p) {
1129:             DMPlexMapToLocalPoint(dm, remoteHash, fcone[p-1], &points[p]);
1130:             if (ierr) break; /* We got a point not in our overlap */
1131:             if (debug) {PetscSynchronizedPrintf(PetscObjectComm((PetscObject) dm), "[%d]  Checking local candidate %D\n", rank, points[p]);}
1132:           }
1133:           if (ierr) continue;
1134:           DMPlexGetJoin(dm, Np, points, &joinSize, &join);
1135:           if (joinSize == 1) {
1136:             PetscSFNode lface;
1137:             PetscSFNode oface;

1139:             /* Always replace with local face */
1140:             lface.rank  = rank;
1141:             lface.index = join[0];
1142:             PetscHashIJKLRemoteIterGet(faceTable, iter, &oface);
1143:             if (debug) {PetscSynchronizedPrintf(PetscObjectComm((PetscObject) dm), "[%d]  Replacing (%D, %D) with local face (%D, %D)\n", rank, oface.index, oface.rank, lface.index, lface.rank);}
1144:             PetscHashIJKLRemoteIterSet(faceTable, iter, lface);
1145:           }
1146:           DMPlexRestoreJoin(dm, Np, points, &joinSize, &join);
1147:         }
1148:       }
1149:       /* Put back faces for this root */
1150:       for (deg = 0; deg < rootdegree[r]; ++deg, ++idx2) {
1151:         PetscInt offset, dof, d;

1153:         PetscSectionGetDof(candidateRemoteSection, idx2, &dof);
1154:         PetscSectionGetOffset(candidateRemoteSection, idx2, &offset);
1155:         /* dof may include many faces from the remote process */
1156:         for (d = 0; d < dof; ++d) {
1157:           const PetscInt         hidx  = offset+d;
1158:           const PetscInt         Np    = candidatesRemote[hidx].index+1;
1159:           const PetscSFNode     *fcone = &candidatesRemote[hidx+2];
1160:           PetscSFNode            fcp0;
1161:           const PetscSFNode      pmax  = {PETSC_MAX_INT, PETSC_MAX_INT};
1162:           PetscHashIJKLRemoteKey key;
1163:           PetscHashIter          iter;
1164:           PetscBool              missing;

1166:           if (debug) {PetscSynchronizedPrintf(PetscObjectComm((PetscObject) dm), "[%d]  Entering face at (%D, %D)\n", rank, r, idx);}
1167:           if (Np > 4) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot handle faces with %D cone points", Np);
1168:           fcp0.rank  = rank;
1169:           fcp0.index = r;
1170:           d += Np;
1171:           /* Find remote face in hash table */
1172:           key.i = fcp0;
1173:           key.j = fcone[0];
1174:           key.k = Np > 2 ? fcone[1] : pmax;
1175:           key.l = Np > 3 ? fcone[2] : pmax;
1176:           PetscSortSFNode(Np, (PetscSFNode *) &key);
1177:           if (debug) {PetscSynchronizedPrintf(PetscObjectComm((PetscObject) dm), "[%d]    key (%D, %D) (%D, %D) (%D, %D) (%D, %D)\n", rank, key.i.rank, key.i.index, key.j.rank, key.j.index, key.k.rank, key.k.index, key.l.rank, key.l.index);}
1178:           PetscHashIJKLRemotePut(faceTable, key, &iter, &missing);
1179:           if (missing) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Root %D Idx %D ought to have an associated face", r, idx2);
1180:           else        {PetscHashIJKLRemoteIterGet(faceTable, iter, &candidatesRemote[hidx]);}
1181:         }
1182:       }
1183:     }
1184:     if (debug) {PetscSynchronizedFlush(PetscObjectComm((PetscObject) dm), NULL);}
1185:     PetscHashIJKLRemoteDestroy(&faceTable);
1186:   }
1187:   /* Step 4: Push back owned faces */
1188:   {
1189:     PetscSF      sfMulti, sfClaims, sfPointNew;
1190:     PetscSFNode *remotePointsNew;
1191:     PetscInt    *remoteOffsets, *localPointsNew;
1192:     PetscInt     pStart, pEnd, r, NlNew, p;

1194:     /* 4) Push claims back to receiver via the MultiSF and derive new pointSF mapping on receiver */
1195:     PetscSFGetMultiSF(pointSF, &sfMulti);
1196:     PetscSectionCreate(comm, &claimSection);
1197:     PetscSFDistributeSection(sfMulti, candidateRemoteSection, &remoteOffsets, claimSection);
1198:     PetscSFCreateSectionSF(sfMulti, candidateRemoteSection, remoteOffsets, claimSection, &sfClaims);
1199:     PetscSectionGetStorageSize(claimSection, &claimsSize);
1200:     PetscMalloc1(claimsSize, &claims);
1201:     for (p = 0; p < claimsSize; ++p) claims[p].rank = -1;
1202:     PetscSFBcastBegin(sfClaims, MPIU_2INT, candidatesRemote, claims,MPI_REPLACE);
1203:     PetscSFBcastEnd(sfClaims, MPIU_2INT, candidatesRemote, claims,MPI_REPLACE);
1204:     PetscSFDestroy(&sfClaims);
1205:     PetscFree(remoteOffsets);
1206:     PetscObjectSetName((PetscObject) claimSection, "Claim Section");
1207:     PetscObjectViewFromOptions((PetscObject) claimSection, NULL, "-petscsection_interp_claim_view");
1208:     SFNodeArrayViewFromOptions(comm, "-petscsection_interp_claim_view", "Claims", NULL, claimsSize, claims);
1209:     /* Step 5) Walk the original section of local supports and add an SF entry for each updated item */
1210:     /* TODO I should not have to do a join here since I already put the face and its cone in the candidate section */
1211:     PetscHMapICreate(&claimshash);
1212:     for (r = 0; r < Nr; ++r) {
1213:       PetscInt dof, off, d;

1215:       if (debug) {PetscSynchronizedPrintf(comm, "[%d]  Checking root for claims %D\n", rank, r);}
1216:       PetscSectionGetDof(candidateSection, r, &dof);
1217:       PetscSectionGetOffset(candidateSection, r, &off);
1218:       for (d = 0; d < dof;) {
1219:         if (claims[off+d].rank >= 0) {
1220:           const PetscInt  faceInd = off+d;
1221:           const PetscInt  Np      = candidates[off+d].index;
1222:           const PetscInt *join    = NULL;
1223:           PetscInt        joinSize, points[1024], c;

1225:           if (debug) {PetscSynchronizedPrintf(comm, "[%d]    Found claim for remote point (%D, %D)\n", rank, claims[faceInd].rank, claims[faceInd].index);}
1226:           points[0] = r;
1227:           if (debug) {PetscSynchronizedPrintf(comm, "[%d]      point %D\n", rank, points[0]);}
1228:           for (c = 0, d += 2; c < Np; ++c, ++d) {
1229:             DMPlexMapToLocalPoint(dm, remoteHash, candidates[off+d], &points[c+1]);
1230:             if (debug) {PetscSynchronizedPrintf(comm, "[%d]      point %D\n", rank, points[c+1]);}
1231:           }
1232:           DMPlexGetJoin(dm, Np+1, points, &joinSize, &join);
1233:           if (joinSize == 1) {
1234:             if (claims[faceInd].rank == rank) {
1235:               if (debug) {PetscSynchronizedPrintf(comm, "[%d]    Ignoring local face %D for non-remote partner\n", rank, join[0]);}
1236:             } else {
1237:               if (debug) {PetscSynchronizedPrintf(comm, "[%d]    Found local face %D\n", rank, join[0]);}
1238:               PetscHMapISet(claimshash, join[0], faceInd);
1239:             }
1240:           } else {
1241:             if (debug) {PetscSynchronizedPrintf(comm, "[%d]    Failed to find face\n", rank);}
1242:           }
1243:           DMPlexRestoreJoin(dm, Np+1, points, &joinSize, &join);
1244:         } else {
1245:           if (debug) {PetscSynchronizedPrintf(comm, "[%d]    No claim for point %D\n", rank, r);}
1246:           d += claims[off+d].index+1;
1247:         }
1248:       }
1249:     }
1250:     if (debug) {PetscSynchronizedFlush(comm, NULL);}
1251:     /* Step 6) Create new pointSF from hashed claims */
1252:     PetscHMapIGetSize(claimshash, &NlNew);
1253:     DMPlexGetChart(dm, &pStart, &pEnd);
1254:     PetscMalloc1(Nl + NlNew, &localPointsNew);
1255:     PetscMalloc1(Nl + NlNew, &remotePointsNew);
1256:     for (l = 0; l < Nl; ++l) {
1257:       localPointsNew[l] = localPoints[l];
1258:       remotePointsNew[l].index = remotePoints[l].index;
1259:       remotePointsNew[l].rank  = remotePoints[l].rank;
1260:     }
1261:     p = Nl;
1262:     PetscHMapIGetKeys(claimshash, &p, localPointsNew);
1263:     /* We sort new points, and assume they are numbered after all existing points */
1264:     PetscSortInt(NlNew, &localPointsNew[Nl]);
1265:     for (p = Nl; p < Nl + NlNew; ++p) {
1266:       PetscInt off;
1267:       PetscHMapIGet(claimshash, localPointsNew[p], &off);
1268:       if (claims[off].rank < 0 || claims[off].index < 0) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid claim for local point %D, (%D, %D)", localPointsNew[p], claims[off].rank, claims[off].index);
1269:       remotePointsNew[p] = claims[off];
1270:     }
1271:     PetscSFCreate(comm, &sfPointNew);
1272:     PetscSFSetGraph(sfPointNew, pEnd-pStart, Nl+NlNew, localPointsNew, PETSC_OWN_POINTER, remotePointsNew, PETSC_OWN_POINTER);
1273:     PetscSFSetUp(sfPointNew);
1274:     DMSetPointSF(dm, sfPointNew);
1275:     PetscObjectViewFromOptions((PetscObject) sfPointNew, NULL, "-petscsf_interp_view");
1276:     PetscSFDestroy(&sfPointNew);
1277:     PetscHMapIDestroy(&claimshash);
1278:   }
1279:   PetscHMapIJDestroy(&remoteHash);
1280:   PetscSectionDestroy(&candidateSection);
1281:   PetscSectionDestroy(&candidateRemoteSection);
1282:   PetscSectionDestroy(&claimSection);
1283:   PetscFree(candidates);
1284:   PetscFree(candidatesRemote);
1285:   PetscFree(claims);
1286:   PetscLogEventEnd(DMPLEX_InterpolateSF,dm,0,0,0);
1287:   return(0);
1288: }

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

1293:   Collective on dm

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

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

1302:   Level: intermediate

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

1307:   Developer Notes:
1308:     It sets plex->interpolated = DMPLEX_INTERPOLATED_FULL.

1310: .seealso: DMPlexUninterpolate(), DMPlexCreateFromCellListPetsc(), DMPlexCopyCoordinates()
1311: @*/
1312: PetscErrorCode DMPlexInterpolate(DM dm, DM *dmInt)
1313: {
1314:   DMPlexInterpolatedFlag interpolated;
1315:   DM             idm, odm = dm;
1316:   PetscSF        sfPoint;
1317:   PetscInt       depth, dim, d;
1318:   const char    *name;
1319:   PetscBool      flg=PETSC_TRUE;

1325:   PetscLogEventBegin(DMPLEX_Interpolate,dm,0,0,0);
1326:   DMPlexGetDepth(dm, &depth);
1327:   DMGetDimension(dm, &dim);
1328:   DMPlexIsInterpolated(dm, &interpolated);
1329:   if (interpolated == DMPLEX_INTERPOLATED_PARTIAL) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Not for partially interpolated meshes");
1330:   if (interpolated == DMPLEX_INTERPOLATED_FULL) {
1331:     PetscObjectReference((PetscObject) dm);
1332:     idm  = dm;
1333:   } else {
1334:     for (d = 1; d < dim; ++d) {
1335:       /* Create interpolated mesh */
1336:       DMCreate(PetscObjectComm((PetscObject)dm), &idm);
1337:       DMSetType(idm, DMPLEX);
1338:       DMSetDimension(idm, dim);
1339:       if (depth > 0) {
1340:         DMPlexInterpolateFaces_Internal(odm, 1, idm);
1341:         DMGetPointSF(odm, &sfPoint);
1342:         {
1343:           /* TODO: We need to systematically fix cases of distributed Plexes with no graph set */
1344:           PetscInt nroots;
1345:           PetscSFGetGraph(sfPoint, &nroots, NULL, NULL, NULL);
1346:           if (nroots >= 0) {DMPlexInterpolatePointSF(idm, sfPoint);}
1347:         }
1348:       }
1349:       if (odm != dm) {DMDestroy(&odm);}
1350:       odm = idm;
1351:     }
1352:     PetscObjectGetName((PetscObject) dm,  &name);
1353:     PetscObjectSetName((PetscObject) idm,  name);
1354:     DMPlexCopyCoordinates(dm, idm);
1355:     DMCopyLabels(dm, idm, PETSC_COPY_VALUES, PETSC_FALSE);
1356:     PetscOptionsGetBool(((PetscObject)dm)->options, ((PetscObject)dm)->prefix, "-dm_plex_interpolate_orient_interfaces", &flg, NULL);
1357:     if (flg) {DMPlexOrientInterface_Internal(idm);}
1358:   }
1359:   {
1360:     PetscBool            isper;
1361:     const PetscReal      *maxCell, *L;
1362:     const DMBoundaryType *bd;

1364:     DMGetPeriodicity(dm,&isper,&maxCell,&L,&bd);
1365:     DMSetPeriodicity(idm,isper,maxCell,L,bd);
1366:   }
1367:   /* This function makes the mesh fully interpolated on all ranks */
1368:   {
1369:     DM_Plex *plex = (DM_Plex *) idm->data;
1370:     plex->interpolated = plex->interpolatedCollective = DMPLEX_INTERPOLATED_FULL;
1371:   }
1372:   *dmInt = idm;
1373:   PetscLogEventEnd(DMPLEX_Interpolate,dm,0,0,0);
1374:   return(0);
1375: }

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

1380:   Collective on dmA

1382:   Input Parameter:
1383: . dmA - The DMPlex object with initial coordinates

1385:   Output Parameter:
1386: . dmB - The DMPlex object with copied coordinates

1388:   Level: intermediate

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

1392: .seealso: DMCopyLabels(), DMGetCoordinates(), DMGetCoordinatesLocal(), DMGetCoordinateDM(), DMGetCoordinateSection()
1393: @*/
1394: PetscErrorCode DMPlexCopyCoordinates(DM dmA, DM dmB)
1395: {
1396:   Vec            coordinatesA, coordinatesB;
1397:   VecType        vtype;
1398:   PetscSection   coordSectionA, coordSectionB;
1399:   PetscScalar   *coordsA, *coordsB;
1400:   PetscInt       spaceDim, Nf, vStartA, vStartB, vEndA, vEndB, coordSizeB, v, d;
1401:   PetscInt       cStartA, cEndA, cStartB, cEndB, cS, cE, cdim;
1402:   PetscBool      lc = PETSC_FALSE;

1408:   if (dmA == dmB) return(0);
1409:   DMGetCoordinateDim(dmA, &cdim);
1410:   DMSetCoordinateDim(dmB, cdim);
1411:   DMPlexGetDepthStratum(dmA, 0, &vStartA, &vEndA);
1412:   DMPlexGetDepthStratum(dmB, 0, &vStartB, &vEndB);
1413:   if ((vEndA-vStartA) != (vEndB-vStartB)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "The number of vertices in first DM %d != %d in the second DM", vEndA-vStartA, vEndB-vStartB);
1414:   /* Copy over discretization if it exists */
1415:   {
1416:     DM                 cdmA, cdmB;
1417:     PetscDS            dsA, dsB;
1418:     PetscObject        objA, objB;
1419:     PetscClassId       idA, idB;
1420:     const PetscScalar *constants;
1421:     PetscInt            cdim, Nc;

1423:     DMGetCoordinateDM(dmA, &cdmA);
1424:     DMGetCoordinateDM(dmB, &cdmB);
1425:     DMGetField(cdmA, 0, NULL, &objA);
1426:     DMGetField(cdmB, 0, NULL, &objB);
1427:     PetscObjectGetClassId(objA, &idA);
1428:     PetscObjectGetClassId(objB, &idB);
1429:     if ((idA == PETSCFE_CLASSID) && (idA != idB)) {
1430:       DMSetField(cdmB, 0, NULL, objA);
1431:       DMCreateDS(cdmB);
1432:       DMGetDS(cdmA, &dsA);
1433:       DMGetDS(cdmB, &dsB);
1434:       PetscDSGetCoordinateDimension(dsA, &cdim);
1435:       PetscDSSetCoordinateDimension(dsB, cdim);
1436:       PetscDSGetConstants(dsA, &Nc, &constants);
1437:       PetscDSSetConstants(dsB, Nc, (PetscScalar *) constants);
1438:     }
1439:   }
1440:   DMPlexGetHeightStratum(dmA, 0, &cStartA, &cEndA);
1441:   DMPlexGetHeightStratum(dmB, 0, &cStartB, &cEndB);
1442:   DMGetCoordinateSection(dmA, &coordSectionA);
1443:   DMGetCoordinateSection(dmB, &coordSectionB);
1444:   if (coordSectionA == coordSectionB) return(0);
1445:   PetscSectionGetNumFields(coordSectionA, &Nf);
1446:   if (!Nf) return(0);
1447:   if (Nf > 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "The number of coordinate fields must be 1, not %D", Nf);
1448:   if (!coordSectionB) {
1449:     PetscInt dim;

1451:     PetscSectionCreate(PetscObjectComm((PetscObject) coordSectionA), &coordSectionB);
1452:     DMGetCoordinateDim(dmA, &dim);
1453:     DMSetCoordinateSection(dmB, dim, coordSectionB);
1454:     PetscObjectDereference((PetscObject) coordSectionB);
1455:   }
1456:   PetscSectionSetNumFields(coordSectionB, 1);
1457:   PetscSectionGetFieldComponents(coordSectionA, 0, &spaceDim);
1458:   PetscSectionSetFieldComponents(coordSectionB, 0, spaceDim);
1459:   PetscSectionGetChart(coordSectionA, &cS, &cE);
1460:   if (cStartA <= cS && cS < cEndA) { /* localized coordinates */
1461:     if ((cEndA-cStartA) != (cEndB-cStartB)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "The number of cells in first DM %D != %D in the second DM", cEndA-cStartA, cEndB-cStartB);
1462:     cS = cS - cStartA + cStartB;
1463:     cE = vEndB;
1464:     lc = PETSC_TRUE;
1465:   } else {
1466:     cS = vStartB;
1467:     cE = vEndB;
1468:   }
1469:   PetscSectionSetChart(coordSectionB, cS, cE);
1470:   for (v = vStartB; v < vEndB; ++v) {
1471:     PetscSectionSetDof(coordSectionB, v, spaceDim);
1472:     PetscSectionSetFieldDof(coordSectionB, v, 0, spaceDim);
1473:   }
1474:   if (lc) { /* localized coordinates */
1475:     PetscInt c;

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

1480:       PetscSectionGetDof(coordSectionA, c + cStartA, &dof);
1481:       PetscSectionSetDof(coordSectionB, c + cStartB, dof);
1482:       PetscSectionSetFieldDof(coordSectionB, c + cStartB, 0, dof);
1483:     }
1484:   }
1485:   PetscSectionSetUp(coordSectionB);
1486:   PetscSectionGetStorageSize(coordSectionB, &coordSizeB);
1487:   DMGetCoordinatesLocal(dmA, &coordinatesA);
1488:   VecCreate(PETSC_COMM_SELF, &coordinatesB);
1489:   PetscObjectSetName((PetscObject) coordinatesB, "coordinates");
1490:   VecSetSizes(coordinatesB, coordSizeB, PETSC_DETERMINE);
1491:   VecGetBlockSize(coordinatesA, &d);
1492:   VecSetBlockSize(coordinatesB, d);
1493:   VecGetType(coordinatesA, &vtype);
1494:   VecSetType(coordinatesB, vtype);
1495:   VecGetArray(coordinatesA, &coordsA);
1496:   VecGetArray(coordinatesB, &coordsB);
1497:   for (v = 0; v < vEndB-vStartB; ++v) {
1498:     PetscInt offA, offB;

1500:     PetscSectionGetOffset(coordSectionA, v + vStartA, &offA);
1501:     PetscSectionGetOffset(coordSectionB, v + vStartB, &offB);
1502:     for (d = 0; d < spaceDim; ++d) {
1503:       coordsB[offB+d] = coordsA[offA+d];
1504:     }
1505:   }
1506:   if (lc) { /* localized coordinates */
1507:     PetscInt c;

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

1512:       PetscSectionGetOffset(coordSectionA, c + cStartA, &offA);
1513:       PetscSectionGetOffset(coordSectionB, c + cStartB, &offB);
1514:       PetscSectionGetDof(coordSectionA, c + cStartA, &dof);
1515:       PetscArraycpy(coordsB + offB,coordsA + offA,dof);
1516:     }
1517:   }
1518:   VecRestoreArray(coordinatesA, &coordsA);
1519:   VecRestoreArray(coordinatesB, &coordsB);
1520:   DMSetCoordinatesLocal(dmB, coordinatesB);
1521:   VecDestroy(&coordinatesB);
1522:   return(0);
1523: }

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

1528:   Collective on dm

1530:   Input Parameter:
1531: . dm - The complete DMPlex object

1533:   Output Parameter:
1534: . dmUnint - The DMPlex object with only cells and vertices

1536:   Level: intermediate

1538:   Notes:
1539:     It does not copy over the coordinates.

1541:   Developer Notes:
1542:     It sets plex->interpolated = DMPLEX_INTERPOLATED_NONE.

1544: .seealso: DMPlexInterpolate(), DMPlexCreateFromCellListPetsc(), DMPlexCopyCoordinates()
1545: @*/
1546: PetscErrorCode DMPlexUninterpolate(DM dm, DM *dmUnint)
1547: {
1548:   DMPlexInterpolatedFlag interpolated;
1549:   DM             udm;
1550:   PetscInt       dim, vStart, vEnd, cStart, cEnd, c, maxConeSize = 0, *cone;

1556:   DMGetDimension(dm, &dim);
1557:   DMPlexIsInterpolated(dm, &interpolated);
1558:   if (interpolated == DMPLEX_INTERPOLATED_PARTIAL) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Not for partially interpolated meshes");
1559:   if (interpolated == DMPLEX_INTERPOLATED_NONE || dim <= 1) {
1560:     /* in case dim <= 1 just keep the DMPLEX_INTERPOLATED_FULL flag */
1561:     PetscObjectReference((PetscObject) dm);
1562:     *dmUnint = dm;
1563:     return(0);
1564:   }
1565:   DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
1566:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
1567:   DMCreate(PetscObjectComm((PetscObject) dm), &udm);
1568:   DMSetType(udm, DMPLEX);
1569:   DMSetDimension(udm, dim);
1570:   DMPlexSetChart(udm, cStart, vEnd);
1571:   for (c = cStart; c < cEnd; ++c) {
1572:     PetscInt *closure = NULL, closureSize, cl, coneSize = 0;

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

1578:       if ((p >= vStart) && (p < vEnd)) ++coneSize;
1579:     }
1580:     DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
1581:     DMPlexSetConeSize(udm, c, coneSize);
1582:     maxConeSize = PetscMax(maxConeSize, coneSize);
1583:   }
1584:   DMSetUp(udm);
1585:   PetscMalloc1(maxConeSize, &cone);
1586:   for (c = cStart; c < cEnd; ++c) {
1587:     PetscInt *closure = NULL, closureSize, cl, coneSize = 0;

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

1593:       if ((p >= vStart) && (p < vEnd)) cone[coneSize++] = p;
1594:     }
1595:     DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
1596:     DMPlexSetCone(udm, c, cone);
1597:   }
1598:   PetscFree(cone);
1599:   DMPlexSymmetrize(udm);
1600:   DMPlexStratify(udm);
1601:   /* Reduce SF */
1602:   {
1603:     PetscSF            sfPoint, sfPointUn;
1604:     const PetscSFNode *remotePoints;
1605:     const PetscInt    *localPoints;
1606:     PetscSFNode       *remotePointsUn;
1607:     PetscInt          *localPointsUn;
1608:     PetscInt           vEnd, numRoots, numLeaves, l;
1609:     PetscInt           numLeavesUn = 0, n = 0;
1610:     PetscErrorCode     ierr;

1612:     /* Get original SF information */
1613:     DMGetPointSF(dm, &sfPoint);
1614:     DMGetPointSF(udm, &sfPointUn);
1615:     DMPlexGetDepthStratum(dm, 0, NULL, &vEnd);
1616:     PetscSFGetGraph(sfPoint, &numRoots, &numLeaves, &localPoints, &remotePoints);
1617:     /* Allocate space for cells and vertices */
1618:     for (l = 0; l < numLeaves; ++l) if (localPoints[l] < vEnd) numLeavesUn++;
1619:     /* Fill in leaves */
1620:     if (vEnd >= 0) {
1621:       PetscMalloc1(numLeavesUn, &remotePointsUn);
1622:       PetscMalloc1(numLeavesUn, &localPointsUn);
1623:       for (l = 0; l < numLeaves; l++) {
1624:         if (localPoints[l] < vEnd) {
1625:           localPointsUn[n]        = localPoints[l];
1626:           remotePointsUn[n].rank  = remotePoints[l].rank;
1627:           remotePointsUn[n].index = remotePoints[l].index;
1628:           ++n;
1629:         }
1630:       }
1631:       if (n != numLeavesUn) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Inconsistent number of leaves %d != %d", n, numLeavesUn);
1632:       PetscSFSetGraph(sfPointUn, vEnd, numLeavesUn, localPointsUn, PETSC_OWN_POINTER, remotePointsUn, PETSC_OWN_POINTER);
1633:     }
1634:   }
1635:   {
1636:     PetscBool            isper;
1637:     const PetscReal      *maxCell, *L;
1638:     const DMBoundaryType *bd;

1640:     DMGetPeriodicity(dm,&isper,&maxCell,&L,&bd);
1641:     DMSetPeriodicity(udm,isper,maxCell,L,bd);
1642:   }
1643:   /* This function makes the mesh fully uninterpolated on all ranks */
1644:   {
1645:     DM_Plex *plex = (DM_Plex *) udm->data;
1646:     plex->interpolated = plex->interpolatedCollective = DMPLEX_INTERPOLATED_NONE;
1647:   }
1648:   *dmUnint = udm;
1649:   return(0);
1650: }

1652: static PetscErrorCode DMPlexIsInterpolated_Internal(DM dm, DMPlexInterpolatedFlag *interpolated)
1653: {
1654:   PetscInt       coneSize, depth, dim, h, p, pStart, pEnd;
1655:   MPI_Comm       comm;

1659:   PetscObjectGetComm((PetscObject)dm, &comm);
1660:   DMPlexGetDepth(dm, &depth);
1661:   DMGetDimension(dm, &dim);

1663:   if (depth == dim) {
1664:     *interpolated = DMPLEX_INTERPOLATED_FULL;
1665:     if (!dim) goto finish;

1667:     /* Check points at height = dim are vertices (have no cones) */
1668:     DMPlexGetHeightStratum(dm, dim, &pStart, &pEnd);
1669:     for (p=pStart; p<pEnd; p++) {
1670:       DMPlexGetConeSize(dm, p, &coneSize);
1671:       if (coneSize) {
1672:         *interpolated = DMPLEX_INTERPOLATED_PARTIAL;
1673:         goto finish;
1674:       }
1675:     }

1677:     /* Check points at height < dim have cones */
1678:     for (h=0; h<dim; h++) {
1679:       DMPlexGetHeightStratum(dm, h, &pStart, &pEnd);
1680:       for (p=pStart; p<pEnd; p++) {
1681:         DMPlexGetConeSize(dm, p, &coneSize);
1682:         if (!coneSize) {
1683:           *interpolated = DMPLEX_INTERPOLATED_PARTIAL;
1684:           goto finish;
1685:         }
1686:       }
1687:     }
1688:   } else if (depth == 1) {
1689:     *interpolated = DMPLEX_INTERPOLATED_NONE;
1690:   } else {
1691:     *interpolated = DMPLEX_INTERPOLATED_PARTIAL;
1692:   }
1693: finish:
1694:   return(0);
1695: }

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

1700:   Not Collective

1702:   Input Parameter:
1703: . dm      - The DM object

1705:   Output Parameter:
1706: . interpolated - Flag whether the DM is interpolated

1708:   Level: intermediate

1710:   Notes:
1711:   Unlike DMPlexIsInterpolatedCollective(), this is NOT collective
1712:   so the results can be different on different ranks in special cases.
1713:   However, DMPlexInterpolate() guarantees the result is the same on all.

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

1717:   Developer Notes:
1718:   Initially, plex->interpolated = DMPLEX_INTERPOLATED_INVALID.

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

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

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

1729: .seealso: DMPlexInterpolate(), DMPlexIsInterpolatedCollective()
1730: @*/
1731: PetscErrorCode DMPlexIsInterpolated(DM dm, DMPlexInterpolatedFlag *interpolated)
1732: {
1733:   DM_Plex        *plex = (DM_Plex *) dm->data;
1734:   PetscErrorCode  ierr;

1739:   if (plex->interpolated < 0) {
1740:     DMPlexIsInterpolated_Internal(dm, &plex->interpolated);
1741:   } else if (PetscDefined (USE_DEBUG)) {
1742:     DMPlexInterpolatedFlag flg;

1744:     DMPlexIsInterpolated_Internal(dm, &flg);
1745:     if (flg != plex->interpolated) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Stashed DMPlexInterpolatedFlag %s is inconsistent with current %s", DMPlexInterpolatedFlags[plex->interpolated], DMPlexInterpolatedFlags[flg]);
1746:   }
1747:   *interpolated = plex->interpolated;
1748:   return(0);
1749: }

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

1754:   Collective

1756:   Input Parameter:
1757: . dm      - The DM object

1759:   Output Parameter:
1760: . interpolated - Flag whether the DM is interpolated

1762:   Level: intermediate

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

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

1769:   Developer Notes:
1770:   Initially, plex->interpolatedCollective = DMPLEX_INTERPOLATED_INVALID.

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

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

1779: .seealso: DMPlexInterpolate(), DMPlexIsInterpolated()
1780: @*/
1781: PetscErrorCode DMPlexIsInterpolatedCollective(DM dm, DMPlexInterpolatedFlag *interpolated)
1782: {
1783:   DM_Plex        *plex = (DM_Plex *) dm->data;
1784:   PetscBool       debug=PETSC_FALSE;
1785:   PetscErrorCode  ierr;

1790:   PetscOptionsGetBool(((PetscObject) dm)->options, ((PetscObject) dm)->prefix, "-dm_plex_is_interpolated_collective_debug", &debug, NULL);
1791:   if (plex->interpolatedCollective < 0) {
1792:     DMPlexInterpolatedFlag  min, max;
1793:     MPI_Comm                comm;

1795:     PetscObjectGetComm((PetscObject)dm, &comm);
1796:     DMPlexIsInterpolated(dm, &plex->interpolatedCollective);
1797:     MPI_Allreduce(&plex->interpolatedCollective, &min, 1, MPIU_ENUM, MPI_MIN, comm);
1798:     MPI_Allreduce(&plex->interpolatedCollective, &max, 1, MPIU_ENUM, MPI_MAX, comm);
1799:     if (min != max) plex->interpolatedCollective = DMPLEX_INTERPOLATED_MIXED;
1800:     if (debug) {
1801:       PetscMPIInt rank;

1803:       MPI_Comm_rank(comm, &rank);
1804:       PetscSynchronizedPrintf(comm, "[%d] interpolated=%s interpolatedCollective=%s\n", rank, DMPlexInterpolatedFlags[plex->interpolated], DMPlexInterpolatedFlags[plex->interpolatedCollective]);
1805:       PetscSynchronizedFlush(comm, PETSC_STDOUT);
1806:     }
1807:   }
1808:   *interpolated = plex->interpolatedCollective;
1809:   return(0);
1810: }