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;

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

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

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

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

315: PetscErrorCode DMPlexRestoreRawFaces_Internal(DM dm, DMPolytopeType ct, const PetscInt cone[], PetscInt *numFaces, const DMPolytopeType *faceTypes[], const PetscInt *faceSizes[], const PetscInt *faces[])
316: {
317:   if (faceTypes) DMRestoreWorkArray(dm, 0, MPIU_INT, (void *) faceTypes);
318:   if (faceSizes) DMRestoreWorkArray(dm, 0, MPIU_INT, (void *) faceSizes);
319:   if (faces)     DMRestoreWorkArray(dm, 0, MPIU_INT, (void *) faces);
320:   return 0;
321: }

323: /* This interpolates faces for cells at some stratum */
324: static PetscErrorCode DMPlexInterpolateFaces_Internal(DM dm, PetscInt cellDepth, DM idm)
325: {
326:   DMLabel        ctLabel;
327:   PetscHashIJKL  faceTable;
328:   PetscInt       faceTypeNum[DM_NUM_POLYTOPES];
329:   PetscInt       depth, d, pStart, Np, cStart, cEnd, c, fStart, fEnd;

331:   DMPlexGetDepth(dm, &depth);
332:   PetscHashIJKLCreate(&faceTable);
333:   PetscArrayzero(faceTypeNum, DM_NUM_POLYTOPES);
334:   DMPlexGetDepthStratum(dm, cellDepth, &cStart, &cEnd);
335:   /* Number new faces and save face vertices in hash table */
336:   DMPlexGetDepthStratum(dm, depth > cellDepth ? cellDepth : 0, NULL, &fStart);
337:   fEnd = fStart;
338:   for (c = cStart; c < cEnd; ++c) {
339:     const PetscInt       *cone, *faceSizes, *faces;
340:     const DMPolytopeType *faceTypes;
341:     DMPolytopeType        ct;
342:     PetscInt              numFaces, cf, foff = 0;

344:     DMPlexGetCellType(dm, c, &ct);
345:     DMPlexGetCone(dm, c, &cone);
346:     DMPlexGetRawFaces_Internal(dm, ct, cone, &numFaces, &faceTypes, &faceSizes, &faces);
347:     for (cf = 0; cf < numFaces; foff += faceSizes[cf], ++cf) {
348:       const PetscInt       faceSize = faceSizes[cf];
349:       const DMPolytopeType faceType = faceTypes[cf];
350:       const PetscInt      *face     = &faces[foff];
351:       PetscHashIJKLKey     key;
352:       PetscHashIter        iter;
353:       PetscBool            missing;

356:       key.i = face[0];
357:       key.j = faceSize > 1 ? face[1] : PETSC_MAX_INT;
358:       key.k = faceSize > 2 ? face[2] : PETSC_MAX_INT;
359:       key.l = faceSize > 3 ? face[3] : PETSC_MAX_INT;
360:       PetscSortInt(faceSize, (PetscInt *) &key);
361:       PetscHashIJKLPut(faceTable, key, &iter, &missing);
362:       if (missing) {
363:         PetscHashIJKLIterSet(faceTable, iter, fEnd++);
364:         ++faceTypeNum[faceType];
365:       }
366:     }
367:     DMPlexRestoreRawFaces_Internal(dm, ct, cone, &numFaces, &faceTypes, &faceSizes, &faces);
368:   }
369:   /* We need to number faces contiguously among types */
370:   {
371:     PetscInt faceTypeStart[DM_NUM_POLYTOPES], ct, numFT = 0;

373:     for (ct = 0; ct < DM_NUM_POLYTOPES; ++ct) {if (faceTypeNum[ct]) ++numFT; faceTypeStart[ct] = 0;}
374:     if (numFT > 1) {
375:       PetscHashIJKLClear(faceTable);
376:       faceTypeStart[0] = fStart;
377:       for (ct = 1; ct < DM_NUM_POLYTOPES; ++ct) faceTypeStart[ct] = faceTypeStart[ct-1] + faceTypeNum[ct-1];
378:       for (c = cStart; c < cEnd; ++c) {
379:         const PetscInt       *cone, *faceSizes, *faces;
380:         const DMPolytopeType *faceTypes;
381:         DMPolytopeType        ct;
382:         PetscInt              numFaces, cf, foff = 0;

384:         DMPlexGetCellType(dm, c, &ct);
385:         DMPlexGetCone(dm, c, &cone);
386:         DMPlexGetRawFaces_Internal(dm, ct, cone, &numFaces, &faceTypes, &faceSizes, &faces);
387:         for (cf = 0; cf < numFaces; foff += faceSizes[cf], ++cf) {
388:           const PetscInt       faceSize = faceSizes[cf];
389:           const DMPolytopeType faceType = faceTypes[cf];
390:           const PetscInt      *face     = &faces[foff];
391:           PetscHashIJKLKey     key;
392:           PetscHashIter        iter;
393:           PetscBool            missing;

396:           key.i = face[0];
397:           key.j = faceSize > 1 ? face[1] : PETSC_MAX_INT;
398:           key.k = faceSize > 2 ? face[2] : PETSC_MAX_INT;
399:           key.l = faceSize > 3 ? face[3] : PETSC_MAX_INT;
400:           PetscSortInt(faceSize, (PetscInt *) &key);
401:           PetscHashIJKLPut(faceTable, key, &iter, &missing);
402:           if (missing) PetscHashIJKLIterSet(faceTable, iter, faceTypeStart[faceType]++);
403:         }
404:         DMPlexRestoreRawFaces_Internal(dm, ct, cone, &numFaces, &faceTypes, &faceSizes, &faces);
405:       }
406:       for (ct = 1; ct < DM_NUM_POLYTOPES; ++ct) {
408:       }
409:     }
410:   }
411:   /* Add new points, always at the end of the numbering */
412:   DMPlexGetChart(dm, &pStart, &Np);
413:   DMPlexSetChart(idm, pStart, Np + (fEnd - fStart));
414:   /* Set cone sizes */
415:   /*   Must create the celltype label here so that we do not automatically try to compute the types */
416:   DMCreateLabel(idm, "celltype");
417:   DMPlexGetCellTypeLabel(idm, &ctLabel);
418:   for (d = 0; d <= depth; ++d) {
419:     DMPolytopeType ct;
420:     PetscInt       coneSize, pStart, pEnd, p;

422:     if (d == cellDepth) continue;
423:     DMPlexGetDepthStratum(dm, d, &pStart, &pEnd);
424:     for (p = pStart; p < pEnd; ++p) {
425:       DMPlexGetConeSize(dm, p, &coneSize);
426:       DMPlexSetConeSize(idm, p, coneSize);
427:       DMPlexGetCellType(dm, p, &ct);
428:       DMPlexSetCellType(idm, p, ct);
429:     }
430:   }
431:   for (c = cStart; c < cEnd; ++c) {
432:     const PetscInt       *cone, *faceSizes, *faces;
433:     const DMPolytopeType *faceTypes;
434:     DMPolytopeType        ct;
435:     PetscInt              numFaces, cf, foff = 0;

437:     DMPlexGetCellType(dm, c, &ct);
438:     DMPlexGetCone(dm, c, &cone);
439:     DMPlexGetRawFaces_Internal(dm, ct, cone, &numFaces, &faceTypes, &faceSizes, &faces);
440:     DMPlexSetCellType(idm, c, ct);
441:     DMPlexSetConeSize(idm, c, numFaces);
442:     for (cf = 0; cf < numFaces; foff += faceSizes[cf], ++cf) {
443:       const PetscInt       faceSize = faceSizes[cf];
444:       const DMPolytopeType faceType = faceTypes[cf];
445:       const PetscInt      *face     = &faces[foff];
446:       PetscHashIJKLKey     key;
447:       PetscHashIter        iter;
448:       PetscBool            missing;
449:       PetscInt             f;

452:       key.i = face[0];
453:       key.j = faceSize > 1 ? face[1] : PETSC_MAX_INT;
454:       key.k = faceSize > 2 ? face[2] : PETSC_MAX_INT;
455:       key.l = faceSize > 3 ? face[3] : PETSC_MAX_INT;
456:       PetscSortInt(faceSize, (PetscInt *) &key);
457:       PetscHashIJKLPut(faceTable, key, &iter, &missing);
459:       PetscHashIJKLIterGet(faceTable, iter, &f);
460:       DMPlexSetConeSize(idm, f, faceSize);
461:       DMPlexSetCellType(idm, f, faceType);
462:     }
463:     DMPlexRestoreRawFaces_Internal(dm, ct, cone, &numFaces, &faceTypes, &faceSizes, &faces);
464:   }
465:   DMSetUp(idm);
466:   /* Initialize cones so we do not need the bash table to tell us that a cone has been set */
467:   {
468:     PetscSection cs;
469:     PetscInt    *cones, csize;

471:     DMPlexGetConeSection(idm, &cs);
472:     DMPlexGetCones(idm, &cones);
473:     PetscSectionGetStorageSize(cs, &csize);
474:     for (c = 0; c < csize; ++c) cones[c] = -1;
475:   }
476:   /* Set cones */
477:   for (d = 0; d <= depth; ++d) {
478:     const PetscInt *cone;
479:     PetscInt        pStart, pEnd, p;

481:     if (d == cellDepth) continue;
482:     DMPlexGetDepthStratum(dm, d, &pStart, &pEnd);
483:     for (p = pStart; p < pEnd; ++p) {
484:       DMPlexGetCone(dm, p, &cone);
485:       DMPlexSetCone(idm, p, cone);
486:       DMPlexGetConeOrientation(dm, p, &cone);
487:       DMPlexSetConeOrientation(idm, p, cone);
488:     }
489:   }
490:   for (c = cStart; c < cEnd; ++c) {
491:     const PetscInt       *cone, *faceSizes, *faces;
492:     const DMPolytopeType *faceTypes;
493:     DMPolytopeType        ct;
494:     PetscInt              numFaces, cf, foff = 0;

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

510:       key.i = face[0];
511:       key.j = faceSize > 1 ? face[1] : PETSC_MAX_INT;
512:       key.k = faceSize > 2 ? face[2] : PETSC_MAX_INT;
513:       key.l = faceSize > 3 ? face[3] : PETSC_MAX_INT;
514:       PetscSortInt(faceSize, (PetscInt *) &key);
515:       PetscHashIJKLPut(faceTable, key, &iter, &missing);
516:       PetscHashIJKLIterGet(faceTable, iter, &f);
517:       DMPlexInsertCone(idm, c, cf, f);
518:       DMPlexGetCone(idm, f, &fcone);
519:       if (fcone[0] < 0) DMPlexSetCone(idm, f, face);
520:       {
521:         const PetscInt *cone;
522:         PetscInt        coneSize, ornt;

524:         DMPlexGetConeSize(idm, f, &coneSize);
525:         DMPlexGetCone(idm, f, &cone);
527:         /* Notice that we have to use vertices here because the lower dimensional faces have not been created yet */
528:         DMPolytopeGetVertexOrientation(faceType, cone, face, &ornt);
529:         DMPlexInsertConeOrientation(idm, c, cf, ornt);
530:       }
531:     }
532:     DMPlexRestoreRawFaces_Internal(dm, ct, cone, &numFaces, &faceTypes, &faceSizes, &faces);
533:   }
534:   PetscHashIJKLDestroy(&faceTable);
535:   DMPlexSymmetrize(idm);
536:   DMPlexStratify(idm);
537:   return 0;
538: }

540: static PetscErrorCode SortRmineRremoteByRemote_Private(PetscSF sf, PetscInt *rmine1[], PetscInt *rremote1[])
541: {
542:   PetscInt            nleaves;
543:   PetscInt            nranks;
544:   const PetscMPIInt  *ranks=NULL;
545:   const PetscInt     *roffset=NULL, *rmine=NULL, *rremote=NULL;
546:   PetscInt            n, o, r;

548:   PetscSFGetRootRanks(sf, &nranks, &ranks, &roffset, &rmine, &rremote);
549:   nleaves = roffset[nranks];
550:   PetscMalloc2(nleaves, rmine1, nleaves, rremote1);
551:   for (r=0; r<nranks; r++) {
552:     /* simultaneously sort rank-wise portions of rmine & rremote by values in rremote
553:        - to unify order with the other side */
554:     o = roffset[r];
555:     n = roffset[r+1] - o;
556:     PetscArraycpy(&(*rmine1)[o], &rmine[o], n);
557:     PetscArraycpy(&(*rremote1)[o], &rremote[o], n);
558:     PetscSortIntWithArray(n, &(*rremote1)[o], &(*rmine1)[o]);
559:   }
560:   return 0;
561: }

563: PetscErrorCode DMPlexOrientInterface_Internal(DM dm)
564: {
565:   PetscSF            sf;
566:   const PetscInt    *locals;
567:   const PetscSFNode *remotes;
568:   const PetscMPIInt *ranks;
569:   const PetscInt    *roffset;
570:   PetscInt          *rmine1, *rremote1; /* rmine and rremote copies simultaneously sorted by rank and rremote */
571:   PetscInt           nroots, p, nleaves, nranks, r, maxConeSize = 0;
572:   PetscInt         (*roots)[4],      (*leaves)[4], mainCone[4];
573:   PetscMPIInt      (*rootsRanks)[4], (*leavesRanks)[4];
574:   MPI_Comm           comm;
575:   PetscMPIInt        rank, size;
576:   PetscInt           debug = 0;

578:   PetscObjectGetComm((PetscObject) dm, &comm);
579:   MPI_Comm_rank(comm, &rank);
580:   MPI_Comm_size(comm, &size);
581:   DMGetPointSF(dm, &sf);
582:   DMViewFromOptions(dm, NULL, "-before_orient_interface_dm_view");
583:   if (PetscDefined(USE_DEBUG)) DMPlexCheckPointSF(dm);
584:   PetscSFGetGraph(sf, &nroots, &nleaves, &locals, &remotes);
585:   if (nroots < 0) return 0;
586:   PetscSFSetUp(sf);
587:   SortRmineRremoteByRemote_Private(sf, &rmine1, &rremote1);
588:   for (p = 0; p < nleaves; ++p) {
589:     PetscInt coneSize;
590:     DMPlexGetConeSize(dm, locals[p], &coneSize);
591:     maxConeSize = PetscMax(maxConeSize, coneSize);
592:   }
594:   PetscMalloc4(nroots, &roots, nroots, &leaves, nroots, &rootsRanks, nroots, &leavesRanks);
595:   for (p = 0; p < nroots; ++p) {
596:     const PetscInt *cone;
597:     PetscInt        coneSize, c, ind0;

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

646:     if (leaves[p][0] < 0) continue; /* Ignore vertices */
647:     DMPlexGetCellType(dm, p, &ct);
648:     DMPlexGetConeSize(dm, p, &coneSize);
649:     DMPlexGetCone(dm, p, &cone);
650:     if (debug) {
651:       PetscCall(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)]",
652:                                       rank, p, cone[0], cone[1], cone[2], cone[3],
653:                                       rootsRanks[p][0], roots[p][0], rootsRanks[p][1], roots[p][1], rootsRanks[p][2], roots[p][2], rootsRanks[p][3], roots[p][3],
654:                                       leavesRanks[p][0], leaves[p][0], leavesRanks[p][1], leaves[p][1], leavesRanks[p][2], leaves[p][2], leavesRanks[p][3], leaves[p][3]));
655:     }
656:     if (leavesRanks[p][0] != rootsRanks[p][0] || leaves[p][0] != roots[p][0] ||
657:         leavesRanks[p][1] != rootsRanks[p][1] || leaves[p][1] != roots[p][1] ||
658:         leavesRanks[p][2] != rootsRanks[p][2] || leaves[p][2] != roots[p][2] ||
659:         leavesRanks[p][3] != rootsRanks[p][3] || leaves[p][3] != roots[p][3]) {
660:       /* Translate these leaves to my cone points; mainCone means desired order p's cone points */
661:       for (c = 0; c < PetscMin(coneSize, 4); ++c) {
662:         PetscInt rS, rN;

664:         if (leavesRanks[p][c] == rank) {
665:           /* A local leaf is just taken as it is */
666:           mainCone[c] = leaves[p][c];
667:           continue;
668:         }
669:         /* Find index of rank leavesRanks[p][c] among remote ranks */
670:         /* No need for PetscMPIIntCast because these integers were originally cast from PetscMPIInt. */
671:         PetscFindMPIInt((PetscMPIInt) leavesRanks[p][c], nranks, ranks, &r);
674:         /* Find point leaves[p][c] among remote points aimed at rank leavesRanks[p][c] */
675:         rS = roffset[r];
676:         rN = roffset[r+1] - rS;
677:         PetscFindInt(leaves[p][c], rN, &rremote1[rS], &ind0);
679:         /* Get the corresponding local point */
680:         mainCone[c] = rmine1[rS + ind0];
681:       }
682:       if (debug) PetscSynchronizedPrintf(comm, " mainCone=[%4D %4D %4D %4D]\n", mainCone[0], mainCone[1], mainCone[2], mainCone[3]);
683:       /* Set the desired order of p's cone points and fix orientations accordingly */
684:       DMPolytopeGetOrientation(ct, cone, mainCone, &o);
685:       DMPlexOrientPoint(dm, p, o);
686:     } else if (debug) PetscSynchronizedPrintf(comm, " ==\n");
687:   }
688:   if (debug) {
689:     PetscSynchronizedFlush(comm, NULL);
690:     MPI_Barrier(comm);
691:   }
692:   DMViewFromOptions(dm, NULL, "-after_orient_interface_dm_view");
693:   PetscFree4(roots, leaves, rootsRanks, leavesRanks);
694:   PetscFree2(rmine1, rremote1);
695:   return 0;
696: }

698: static PetscErrorCode IntArrayViewFromOptions(MPI_Comm comm, const char opt[], const char name[], const char idxname[], const char valname[], PetscInt n, const PetscInt a[])
699: {
700:   PetscInt       idx;
701:   PetscMPIInt    rank;
702:   PetscBool      flg;

704:   PetscOptionsHasName(NULL, NULL, opt, &flg);
705:   if (!flg) return 0;
706:   MPI_Comm_rank(comm, &rank);
707:   PetscSynchronizedPrintf(comm, "[%d]%s:\n", rank, name);
708:   for (idx = 0; idx < n; ++idx) PetscSynchronizedPrintf(comm, "[%d]%s %D %s %D\n", rank, idxname, idx, valname, a[idx]);
709:   PetscSynchronizedFlush(comm, NULL);
710:   return 0;
711: }

713: static PetscErrorCode SFNodeArrayViewFromOptions(MPI_Comm comm, const char opt[], const char name[], const char idxname[], PetscInt n, const PetscSFNode a[])
714: {
715:   PetscInt       idx;
716:   PetscMPIInt    rank;
717:   PetscBool      flg;

719:   PetscOptionsHasName(NULL, NULL, opt, &flg);
720:   if (!flg) return 0;
721:   MPI_Comm_rank(comm, &rank);
722:   PetscSynchronizedPrintf(comm, "[%d]%s:\n", rank, name);
723:   if (idxname) {
724:     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);
725:   } else {
726:     for (idx = 0; idx < n; ++idx) PetscSynchronizedPrintf(comm, "[%d]rank %D index %D\n", rank, a[idx].rank, a[idx].index);
727:   }
728:   PetscSynchronizedFlush(comm, NULL);
729:   return 0;
730: }

732: static PetscErrorCode DMPlexMapToLocalPoint(DM dm, PetscHMapIJ remotehash, PetscSFNode remotePoint, PetscInt *localPoint, PetscBool *mapFailed)
733: {
734:   PetscSF         sf;
735:   const PetscInt *locals;
736:   PetscMPIInt     rank;

738:   MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank);
739:   DMGetPointSF(dm, &sf);
740:   PetscSFGetGraph(sf, NULL, NULL, &locals, NULL);
741:   if (mapFailed) *mapFailed = PETSC_FALSE;
742:   if (remotePoint.rank == rank) {
743:     *localPoint = remotePoint.index;
744:   } else {
745:     PetscHashIJKey key;
746:     PetscInt       l;

748:     key.i = remotePoint.index;
749:     key.j = remotePoint.rank;
750:     PetscHMapIJGet(remotehash, key, &l);
751:     if (l >= 0) {
752:       *localPoint = locals[l];
753:     } else if (mapFailed) *mapFailed = PETSC_TRUE;
754:   }
755:   return 0;
756: }

758: static PetscErrorCode DMPlexMapToGlobalPoint(DM dm, PetscInt localPoint, PetscSFNode *remotePoint, PetscBool *mapFailed)
759: {
760:   PetscSF            sf;
761:   const PetscInt    *locals, *rootdegree;
762:   const PetscSFNode *remotes;
763:   PetscInt           Nl, l;
764:   PetscMPIInt        rank;

766:   if (mapFailed) *mapFailed = PETSC_FALSE;
767:   MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank);
768:   DMGetPointSF(dm, &sf);
769:   PetscSFGetGraph(sf, NULL, &Nl, &locals, &remotes);
770:   if (Nl < 0) goto owned;
771:   PetscSFComputeDegreeBegin(sf, &rootdegree);
772:   PetscSFComputeDegreeEnd(sf, &rootdegree);
773:   if (rootdegree[localPoint]) goto owned;
774:   PetscFindInt(localPoint, Nl, locals, &l);
775:   if (l < 0) {if (mapFailed) *mapFailed = PETSC_TRUE;}
776:   else *remotePoint = remotes[l];
777:   return 0;
778:   owned:
779:   remotePoint->rank  = rank;
780:   remotePoint->index = localPoint;
781:   return 0;
782: }

784: static PetscErrorCode DMPlexPointIsShared(DM dm, PetscInt p, PetscBool *isShared)
785: {
786:   PetscSF         sf;
787:   const PetscInt *locals, *rootdegree;
788:   PetscInt        Nl, idx;

790:   *isShared = PETSC_FALSE;
791:   DMGetPointSF(dm, &sf);
792:   PetscSFGetGraph(sf, NULL, &Nl, &locals, NULL);
793:   if (Nl < 0) return 0;
794:   PetscFindInt(p, Nl, locals, &idx);
795:   if (idx >= 0) {*isShared = PETSC_TRUE; return 0;}
796:   PetscSFComputeDegreeBegin(sf, &rootdegree);
797:   PetscSFComputeDegreeEnd(sf, &rootdegree);
798:   if (rootdegree[p] > 0) *isShared = PETSC_TRUE;
799:   return 0;
800: }

802: static PetscErrorCode DMPlexConeIsShared(DM dm, PetscInt p, PetscBool *isShared)
803: {
804:   const PetscInt *cone;
805:   PetscInt        coneSize, c;
806:   PetscBool       cShared = PETSC_TRUE;

808:   DMPlexGetConeSize(dm, p, &coneSize);
809:   DMPlexGetCone(dm, p, &cone);
810:   for (c = 0; c < coneSize; ++c) {
811:     PetscBool pointShared;

813:     DMPlexPointIsShared(dm, cone[c], &pointShared);
814:     cShared = (PetscBool) (cShared && pointShared);
815:   }
816:   *isShared = coneSize ? cShared : PETSC_FALSE;
817:   return 0;
818: }

820: static PetscErrorCode DMPlexGetConeMinimum(DM dm, PetscInt p, PetscSFNode *cpmin)
821: {
822:   const PetscInt *cone;
823:   PetscInt        coneSize, c;
824:   PetscSFNode     cmin = {PETSC_MAX_INT, PETSC_MAX_INT}, missing = {-1, -1};

826:   DMPlexGetConeSize(dm, p, &coneSize);
827:   DMPlexGetCone(dm, p, &cone);
828:   for (c = 0; c < coneSize; ++c) {
829:     PetscSFNode rcp;
830:     PetscBool   mapFailed;

832:     DMPlexMapToGlobalPoint(dm, cone[c], &rcp, &mapFailed);
833:     if (mapFailed) {
834:       cmin = missing;
835:     } else {
836:       cmin = (rcp.rank < cmin.rank) || (rcp.rank == cmin.rank && rcp.index < cmin.index) ? rcp : cmin;
837:     }
838:   }
839:   *cpmin = coneSize ? cmin : missing;
840:   return 0;
841: }

843: /*
844:   Each shared face has an entry in the candidates array:
845:     (-1, coneSize-1), {(global cone point)}
846:   where the set is missing the point p which we use as the key for the face
847: */
848: static PetscErrorCode DMPlexAddSharedFace_Private(DM dm, PetscSection candidateSection, PetscSFNode candidates[], PetscHMapIJ faceHash, PetscInt p, PetscBool debug)
849: {
850:   MPI_Comm        comm;
851:   const PetscInt *support;
852:   PetscInt        supportSize, s, off = 0, idx = 0, overlap, cellHeight, height;
853:   PetscMPIInt     rank;

855:   PetscObjectGetComm((PetscObject) dm, &comm);
856:   MPI_Comm_rank(comm, &rank);
857:   DMPlexGetOverlap(dm, &overlap);
858:   DMPlexGetVTKCellHeight(dm, &cellHeight);
859:   DMPlexGetPointHeight(dm, p, &height);
860:   if (!overlap && height <= cellHeight+1) {
861:     /* cells can't be shared for non-overlapping meshes */
862:     if (debug) PetscSynchronizedPrintf(comm, "[%d]    Skipping face %D to avoid adding cell to hashmap since this is nonoverlapping mesh\n", rank, p);
863:     return 0;
864:   }
865:   DMPlexGetSupportSize(dm, p, &supportSize);
866:   DMPlexGetSupport(dm, p, &support);
867:   if (candidates) PetscSectionGetOffset(candidateSection, p, &off);
868:   for (s = 0; s < supportSize; ++s) {
869:     const PetscInt  face = support[s];
870:     const PetscInt *cone;
871:     PetscSFNode     cpmin={-1,-1}, rp={-1,-1};
872:     PetscInt        coneSize, c, f;
873:     PetscBool       isShared = PETSC_FALSE;
874:     PetscHashIJKey  key;

876:     /* Only add point once */
877:     if (debug) PetscSynchronizedPrintf(comm, "[%d]    Support face %D\n", rank, face);
878:     key.i = p;
879:     key.j = face;
880:     PetscHMapIJGet(faceHash, key, &f);
881:     if (f >= 0) continue;
882:     DMPlexConeIsShared(dm, face, &isShared);
883:     DMPlexGetConeMinimum(dm, face, &cpmin);
884:     DMPlexMapToGlobalPoint(dm, p, &rp, NULL);
885:     if (debug) {
886:       PetscSynchronizedPrintf(comm, "[%d]      Face point %D is shared: %d\n", rank, face, (int) isShared);
887:       PetscSynchronizedPrintf(comm, "[%d]      Global point (%D, %D) Min Cone Point (%D, %D)\n", rank, rp.rank, rp.index, cpmin.rank, cpmin.index);
888:     }
889:     if (isShared && (rp.rank == cpmin.rank && rp.index == cpmin.index)) {
890:       PetscHMapIJSet(faceHash, key, p);
891:       if (candidates) {
892:         if (debug) PetscSynchronizedPrintf(comm, "[%d]    Adding shared face %D at idx %D\n[%d]     ", rank, face, idx, rank);
893:         DMPlexGetConeSize(dm, face, &coneSize);
894:         DMPlexGetCone(dm, face, &cone);
895:         candidates[off+idx].rank    = -1;
896:         candidates[off+idx++].index = coneSize-1;
897:         candidates[off+idx].rank    = rank;
898:         candidates[off+idx++].index = face;
899:         for (c = 0; c < coneSize; ++c) {
900:           const PetscInt cp = cone[c];

902:           if (cp == p) continue;
903:           DMPlexMapToGlobalPoint(dm, cp, &candidates[off+idx], NULL);
904:           if (debug) PetscSynchronizedPrintf(comm, " (%D,%D)", candidates[off+idx].rank, candidates[off+idx].index);
905:           ++idx;
906:         }
907:         if (debug) PetscSynchronizedPrintf(comm, "\n");
908:       } else {
909:         /* Add cone size to section */
910:         if (debug) PetscSynchronizedPrintf(comm, "[%d]    Scheduling shared face %D\n", rank, face);
911:         DMPlexGetConeSize(dm, face, &coneSize);
912:         PetscHMapIJSet(faceHash, key, p);
913:         PetscSectionAddDof(candidateSection, p, coneSize+1);
914:       }
915:     }
916:   }
917:   return 0;
918: }

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

923:   Collective on dm

925:   Input Parameters:
926: + dm      - The interpolated DM
927: - pointSF - The initial SF without interpolated points

929:   Output Parameter:
930: . pointSF - The SF including interpolated points

932:   Level: developer

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

936: .seealso: DMPlexInterpolate(), DMPlexUninterpolate()
937: @*/
938: PetscErrorCode DMPlexInterpolatePointSF(DM dm, PetscSF pointSF)
939: {
940:   MPI_Comm           comm;
941:   PetscHMapIJ        remoteHash;
942:   PetscHMapI         claimshash;
943:   PetscSection       candidateSection, candidateRemoteSection, claimSection;
944:   PetscSFNode       *candidates, *candidatesRemote, *claims;
945:   const PetscInt    *localPoints, *rootdegree;
946:   const PetscSFNode *remotePoints;
947:   PetscInt           ov, Nr, r, Nl, l;
948:   PetscInt           candidatesSize, candidatesRemoteSize, claimsSize;
949:   PetscBool          flg, debug = PETSC_FALSE;
950:   PetscMPIInt        rank;

954:   DMPlexIsDistributed(dm, &flg);
955:   if (!flg) return 0;
956:   /* Set initial SF so that lower level queries work */
957:   DMSetPointSF(dm, pointSF);
958:   PetscObjectGetComm((PetscObject) dm, &comm);
959:   MPI_Comm_rank(comm, &rank);
960:   DMPlexGetOverlap(dm, &ov);
962:   PetscOptionsHasName(NULL, ((PetscObject) dm)->prefix, "-dmplex_interp_debug", &debug);
963:   PetscObjectViewFromOptions((PetscObject) dm, NULL, "-dm_interp_pre_view");
964:   PetscObjectViewFromOptions((PetscObject) pointSF, NULL, "-petscsf_interp_pre_view");
965:   PetscLogEventBegin(DMPLEX_InterpolateSF,dm,0,0,0);
966:   /* Step 0: Precalculations */
967:   PetscSFGetGraph(pointSF, &Nr, &Nl, &localPoints, &remotePoints);
969:   PetscHMapIJCreate(&remoteHash);
970:   for (l = 0; l < Nl; ++l) {
971:     PetscHashIJKey key;
972:     key.i = remotePoints[l].index;
973:     key.j = remotePoints[l].rank;
974:     PetscHMapIJSet(remoteHash, key, l);
975:   }
976:   /*   Compute root degree to identify shared points */
977:   PetscSFComputeDegreeBegin(pointSF, &rootdegree);
978:   PetscSFComputeDegreeEnd(pointSF, &rootdegree);
979:   IntArrayViewFromOptions(comm, "-interp_root_degree_view", "Root degree", "point", "degree", Nr, rootdegree);
980:   /*
981:   1) Loop over each leaf point $p$ at depth $d$ in the SF
982:   \item Get set $F(p)$ of faces $f$ in the support of $p$ for which
983:   \begin{itemize}
984:     \item all cone points of $f$ are shared
985:     \item $p$ is the cone point with smallest canonical number
986:   \end{itemize}
987:   \item Send $F(p)$ and the cone of each face to the active root point $r(p)$
988:   \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
989:   \item Send the root face from the root back to all leaf process
990:   \item Leaf processes add the shared face to the SF
991:   */
992:   /* Step 1: Construct section+SFNode array
993:        The section has entries for all shared faces for which we have a leaf point in the cone
994:        The array holds candidate shared faces, each face is refered to by the leaf point */
995:   PetscSectionCreate(comm, &candidateSection);
996:   PetscSectionSetChart(candidateSection, 0, Nr);
997:   {
998:     PetscHMapIJ faceHash;

1000:     PetscHMapIJCreate(&faceHash);
1001:     for (l = 0; l < Nl; ++l) {
1002:       const PetscInt p = localPoints[l];

1004:       if (debug) PetscSynchronizedPrintf(comm, "[%d]  First pass leaf point %D\n", rank, p);
1005:       DMPlexAddSharedFace_Private(dm, candidateSection, NULL, faceHash, p, debug);
1006:     }
1007:     PetscHMapIJClear(faceHash);
1008:     PetscSectionSetUp(candidateSection);
1009:     PetscSectionGetStorageSize(candidateSection, &candidatesSize);
1010:     PetscMalloc1(candidatesSize, &candidates);
1011:     for (l = 0; l < Nl; ++l) {
1012:       const PetscInt p = localPoints[l];

1014:       if (debug) PetscSynchronizedPrintf(comm, "[%d]  Second pass leaf point %D\n", rank, p);
1015:       DMPlexAddSharedFace_Private(dm, candidateSection, candidates, faceHash, p, debug);
1016:     }
1017:     PetscHMapIJDestroy(&faceHash);
1018:     if (debug) PetscSynchronizedFlush(comm, NULL);
1019:   }
1020:   PetscObjectSetName((PetscObject) candidateSection, "Candidate Section");
1021:   PetscObjectViewFromOptions((PetscObject) candidateSection, NULL, "-petscsection_interp_candidate_view");
1022:   SFNodeArrayViewFromOptions(comm, "-petscsection_interp_candidate_view", "Candidates", NULL, candidatesSize, candidates);
1023:   /* Step 2: Gather candidate section / array pair into the root partition via inverse(multi(pointSF)). */
1024:   /*   Note that this section is indexed by offsets into leaves, not by point number */
1025:   {
1026:     PetscSF   sfMulti, sfInverse, sfCandidates;
1027:     PetscInt *remoteOffsets;

1029:     PetscSFGetMultiSF(pointSF, &sfMulti);
1030:     PetscSFCreateInverseSF(sfMulti, &sfInverse);
1031:     PetscSectionCreate(comm, &candidateRemoteSection);
1032:     PetscSFDistributeSection(sfInverse, candidateSection, &remoteOffsets, candidateRemoteSection);
1033:     PetscSFCreateSectionSF(sfInverse, candidateSection, remoteOffsets, candidateRemoteSection, &sfCandidates);
1034:     PetscSectionGetStorageSize(candidateRemoteSection, &candidatesRemoteSize);
1035:     PetscMalloc1(candidatesRemoteSize, &candidatesRemote);
1036:     PetscSFBcastBegin(sfCandidates, MPIU_2INT, candidates, candidatesRemote,MPI_REPLACE);
1037:     PetscSFBcastEnd(sfCandidates, MPIU_2INT, candidates, candidatesRemote,MPI_REPLACE);
1038:     PetscSFDestroy(&sfInverse);
1039:     PetscSFDestroy(&sfCandidates);
1040:     PetscFree(remoteOffsets);

1042:     PetscObjectSetName((PetscObject) candidateRemoteSection, "Remote Candidate Section");
1043:     PetscObjectViewFromOptions((PetscObject) candidateRemoteSection, NULL, "-petscsection_interp_candidate_remote_view");
1044:     SFNodeArrayViewFromOptions(comm, "-petscsection_interp_candidate_remote_view", "Remote Candidates", NULL, candidatesRemoteSize, candidatesRemote);
1045:   }
1046:   /* 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 */
1047:   {
1048:     PetscHashIJKLRemote faceTable;
1049:     PetscInt            idx, idx2;

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

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

1059:         PetscSectionGetDof(candidateRemoteSection, idx, &dof);
1060:         PetscSectionGetOffset(candidateRemoteSection, idx, &offset);
1061:         /* dof may include many faces from the remote process */
1062:         for (d = 0; d < dof; ++d) {
1063:           const PetscInt         hidx  = offset+d;
1064:           const PetscInt         Np    = candidatesRemote[hidx].index+1;
1065:           const PetscSFNode      rface = candidatesRemote[hidx+1];
1066:           const PetscSFNode     *fcone = &candidatesRemote[hidx+2];
1067:           PetscSFNode            fcp0;
1068:           const PetscSFNode      pmax  = {PETSC_MAX_INT, PETSC_MAX_INT};
1069:           const PetscInt        *join  = NULL;
1070:           PetscHashIJKLRemoteKey key;
1071:           PetscHashIter          iter;
1072:           PetscBool              missing,mapToLocalPointFailed = PETSC_FALSE;
1073:           PetscInt               points[1024], p, joinSize;

1075:           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);
1077:           fcp0.rank  = rank;
1078:           fcp0.index = r;
1079:           d += Np;
1080:           /* Put remote face in hash table */
1081:           key.i = fcp0;
1082:           key.j = fcone[0];
1083:           key.k = Np > 2 ? fcone[1] : pmax;
1084:           key.l = Np > 3 ? fcone[2] : pmax;
1085:           PetscSortSFNode(Np, (PetscSFNode *) &key);
1086:           PetscHashIJKLRemotePut(faceTable, key, &iter, &missing);
1087:           if (missing) {
1088:             if (debug) PetscSynchronizedPrintf(PetscObjectComm((PetscObject) dm), "[%d]  Setting remote face (%D, %D)\n", rank, rface.index, rface.rank);
1089:             PetscHashIJKLRemoteIterSet(faceTable, iter, rface);
1090:           } else {
1091:             PetscSFNode oface;

1093:             PetscHashIJKLRemoteIterGet(faceTable, iter, &oface);
1094:             if ((rface.rank < oface.rank) || (rface.rank == oface.rank && rface.index < oface.index)) {
1095:               if (debug) PetscSynchronizedPrintf(PetscObjectComm((PetscObject) dm), "[%d]  Replacing with remote face (%D, %D)\n", rank, rface.index, rface.rank);
1096:               PetscHashIJKLRemoteIterSet(faceTable, iter, rface);
1097:             }
1098:           }
1099:           /* Check for local face */
1100:           points[0] = r;
1101:           for (p = 1; p < Np; ++p) {
1102:             DMPlexMapToLocalPoint(dm, remoteHash, fcone[p-1], &points[p], &mapToLocalPointFailed);
1103:             if (mapToLocalPointFailed) break; /* We got a point not in our overlap */
1104:             if (debug) PetscSynchronizedPrintf(PetscObjectComm((PetscObject) dm), "[%d]  Checking local candidate %D\n", rank, points[p]);
1105:           }
1106:           if (mapToLocalPointFailed) continue;
1107:           DMPlexGetJoin(dm, Np, points, &joinSize, &join);
1108:           if (joinSize == 1) {
1109:             PetscSFNode lface;
1110:             PetscSFNode oface;

1112:             /* Always replace with local face */
1113:             lface.rank  = rank;
1114:             lface.index = join[0];
1115:             PetscHashIJKLRemoteIterGet(faceTable, iter, &oface);
1116:             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);
1117:             PetscHashIJKLRemoteIterSet(faceTable, iter, lface);
1118:           }
1119:           DMPlexRestoreJoin(dm, Np, points, &joinSize, &join);
1120:         }
1121:       }
1122:       /* Put back faces for this root */
1123:       for (deg = 0; deg < rootdegree[r]; ++deg, ++idx2) {
1124:         PetscInt offset, dof, d;

1126:         PetscSectionGetDof(candidateRemoteSection, idx2, &dof);
1127:         PetscSectionGetOffset(candidateRemoteSection, idx2, &offset);
1128:         /* dof may include many faces from the remote process */
1129:         for (d = 0; d < dof; ++d) {
1130:           const PetscInt         hidx  = offset+d;
1131:           const PetscInt         Np    = candidatesRemote[hidx].index+1;
1132:           const PetscSFNode     *fcone = &candidatesRemote[hidx+2];
1133:           PetscSFNode            fcp0;
1134:           const PetscSFNode      pmax  = {PETSC_MAX_INT, PETSC_MAX_INT};
1135:           PetscHashIJKLRemoteKey key;
1136:           PetscHashIter          iter;
1137:           PetscBool              missing;

1139:           if (debug) PetscSynchronizedPrintf(PetscObjectComm((PetscObject) dm), "[%d]  Entering face at (%D, %D)\n", rank, r, idx);
1141:           fcp0.rank  = rank;
1142:           fcp0.index = r;
1143:           d += Np;
1144:           /* Find remote face in hash table */
1145:           key.i = fcp0;
1146:           key.j = fcone[0];
1147:           key.k = Np > 2 ? fcone[1] : pmax;
1148:           key.l = Np > 3 ? fcone[2] : pmax;
1149:           PetscSortSFNode(Np, (PetscSFNode *) &key);
1150:           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);
1151:           PetscHashIJKLRemotePut(faceTable, key, &iter, &missing);
1153:           else        PetscHashIJKLRemoteIterGet(faceTable, iter, &candidatesRemote[hidx]);
1154:         }
1155:       }
1156:     }
1157:     if (debug) PetscSynchronizedFlush(PetscObjectComm((PetscObject) dm), NULL);
1158:     PetscHashIJKLRemoteDestroy(&faceTable);
1159:   }
1160:   /* Step 4: Push back owned faces */
1161:   {
1162:     PetscSF      sfMulti, sfClaims, sfPointNew;
1163:     PetscSFNode *remotePointsNew;
1164:     PetscInt    *remoteOffsets, *localPointsNew;
1165:     PetscInt     pStart, pEnd, r, NlNew, p;

1167:     /* 4) Push claims back to receiver via the MultiSF and derive new pointSF mapping on receiver */
1168:     PetscSFGetMultiSF(pointSF, &sfMulti);
1169:     PetscSectionCreate(comm, &claimSection);
1170:     PetscSFDistributeSection(sfMulti, candidateRemoteSection, &remoteOffsets, claimSection);
1171:     PetscSFCreateSectionSF(sfMulti, candidateRemoteSection, remoteOffsets, claimSection, &sfClaims);
1172:     PetscSectionGetStorageSize(claimSection, &claimsSize);
1173:     PetscMalloc1(claimsSize, &claims);
1174:     for (p = 0; p < claimsSize; ++p) claims[p].rank = -1;
1175:     PetscSFBcastBegin(sfClaims, MPIU_2INT, candidatesRemote, claims,MPI_REPLACE);
1176:     PetscSFBcastEnd(sfClaims, MPIU_2INT, candidatesRemote, claims,MPI_REPLACE);
1177:     PetscSFDestroy(&sfClaims);
1178:     PetscFree(remoteOffsets);
1179:     PetscObjectSetName((PetscObject) claimSection, "Claim Section");
1180:     PetscObjectViewFromOptions((PetscObject) claimSection, NULL, "-petscsection_interp_claim_view");
1181:     SFNodeArrayViewFromOptions(comm, "-petscsection_interp_claim_view", "Claims", NULL, claimsSize, claims);
1182:     /* Step 5) Walk the original section of local supports and add an SF entry for each updated item */
1183:     /* TODO I should not have to do a join here since I already put the face and its cone in the candidate section */
1184:     PetscHMapICreate(&claimshash);
1185:     for (r = 0; r < Nr; ++r) {
1186:       PetscInt dof, off, d;

1188:       if (debug) PetscSynchronizedPrintf(comm, "[%d]  Checking root for claims %D\n", rank, r);
1189:       PetscSectionGetDof(candidateSection, r, &dof);
1190:       PetscSectionGetOffset(candidateSection, r, &off);
1191:       for (d = 0; d < dof;) {
1192:         if (claims[off+d].rank >= 0) {
1193:           const PetscInt  faceInd = off+d;
1194:           const PetscInt  Np      = candidates[off+d].index;
1195:           const PetscInt *join    = NULL;
1196:           PetscInt        joinSize, points[1024], c;

1198:           if (debug) PetscSynchronizedPrintf(comm, "[%d]    Found claim for remote point (%D, %D)\n", rank, claims[faceInd].rank, claims[faceInd].index);
1199:           points[0] = r;
1200:           if (debug) PetscSynchronizedPrintf(comm, "[%d]      point %D\n", rank, points[0]);
1201:           for (c = 0, d += 2; c < Np; ++c, ++d) {
1202:             DMPlexMapToLocalPoint(dm, remoteHash, candidates[off+d], &points[c+1], NULL);
1203:             if (debug) PetscSynchronizedPrintf(comm, "[%d]      point %D\n", rank, points[c+1]);
1204:           }
1205:           DMPlexGetJoin(dm, Np+1, points, &joinSize, &join);
1206:           if (joinSize == 1) {
1207:             if (claims[faceInd].rank == rank) {
1208:               if (debug) PetscSynchronizedPrintf(comm, "[%d]    Ignoring local face %D for non-remote partner\n", rank, join[0]);
1209:             } else {
1210:               if (debug) PetscSynchronizedPrintf(comm, "[%d]    Found local face %D\n", rank, join[0]);
1211:               PetscHMapISet(claimshash, join[0], faceInd);
1212:             }
1213:           } else {
1214:             if (debug) PetscSynchronizedPrintf(comm, "[%d]    Failed to find face\n", rank);
1215:           }
1216:           DMPlexRestoreJoin(dm, Np+1, points, &joinSize, &join);
1217:         } else {
1218:           if (debug) PetscSynchronizedPrintf(comm, "[%d]    No claim for point %D\n", rank, r);
1219:           d += claims[off+d].index+1;
1220:         }
1221:       }
1222:     }
1223:     if (debug) PetscSynchronizedFlush(comm, NULL);
1224:     /* Step 6) Create new pointSF from hashed claims */
1225:     PetscHMapIGetSize(claimshash, &NlNew);
1226:     DMPlexGetChart(dm, &pStart, &pEnd);
1227:     PetscMalloc1(Nl + NlNew, &localPointsNew);
1228:     PetscMalloc1(Nl + NlNew, &remotePointsNew);
1229:     for (l = 0; l < Nl; ++l) {
1230:       localPointsNew[l] = localPoints[l];
1231:       remotePointsNew[l].index = remotePoints[l].index;
1232:       remotePointsNew[l].rank  = remotePoints[l].rank;
1233:     }
1234:     p = Nl;
1235:     PetscHMapIGetKeys(claimshash, &p, localPointsNew);
1236:     /* We sort new points, and assume they are numbered after all existing points */
1237:     PetscSortInt(NlNew, &localPointsNew[Nl]);
1238:     for (p = Nl; p < Nl + NlNew; ++p) {
1239:       PetscInt off;
1240:       PetscHMapIGet(claimshash, localPointsNew[p], &off);
1242:       remotePointsNew[p] = claims[off];
1243:     }
1244:     PetscSFCreate(comm, &sfPointNew);
1245:     PetscSFSetGraph(sfPointNew, pEnd-pStart, Nl+NlNew, localPointsNew, PETSC_OWN_POINTER, remotePointsNew, PETSC_OWN_POINTER);
1246:     PetscSFSetUp(sfPointNew);
1247:     DMSetPointSF(dm, sfPointNew);
1248:     PetscObjectViewFromOptions((PetscObject) sfPointNew, NULL, "-petscsf_interp_view");
1249:     PetscSFDestroy(&sfPointNew);
1250:     PetscHMapIDestroy(&claimshash);
1251:   }
1252:   PetscHMapIJDestroy(&remoteHash);
1253:   PetscSectionDestroy(&candidateSection);
1254:   PetscSectionDestroy(&candidateRemoteSection);
1255:   PetscSectionDestroy(&claimSection);
1256:   PetscFree(candidates);
1257:   PetscFree(candidatesRemote);
1258:   PetscFree(claims);
1259:   PetscLogEventEnd(DMPLEX_InterpolateSF,dm,0,0,0);
1260:   return 0;
1261: }

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

1266:   Collective on dm

1268:   Input Parameters:
1269: + dm - The DMPlex object with only cells and vertices
1270: - dmInt - The interpolated DM

1272:   Output Parameter:
1273: . dmInt - The complete DMPlex object

1275:   Level: intermediate

1277:   Notes:
1278:     It does not copy over the coordinates.

1280:   Developer Notes:
1281:     It sets plex->interpolated = DMPLEX_INTERPOLATED_FULL.

1283: .seealso: DMPlexUninterpolate(), DMPlexCreateFromCellListPetsc(), DMPlexCopyCoordinates()
1284: @*/
1285: PetscErrorCode DMPlexInterpolate(DM dm, DM *dmInt)
1286: {
1287:   DMPlexInterpolatedFlag interpolated;
1288:   DM             idm, odm = dm;
1289:   PetscSF        sfPoint;
1290:   PetscInt       depth, dim, d;
1291:   const char    *name;
1292:   PetscBool      flg=PETSC_TRUE;

1296:   PetscLogEventBegin(DMPLEX_Interpolate,dm,0,0,0);
1297:   DMPlexGetDepth(dm, &depth);
1298:   DMGetDimension(dm, &dim);
1299:   DMPlexIsInterpolated(dm, &interpolated);
1301:   if (interpolated == DMPLEX_INTERPOLATED_FULL) {
1302:     PetscObjectReference((PetscObject) dm);
1303:     idm  = dm;
1304:   } else {
1305:     for (d = 1; d < dim; ++d) {
1306:       /* Create interpolated mesh */
1307:       DMCreate(PetscObjectComm((PetscObject)dm), &idm);
1308:       DMSetType(idm, DMPLEX);
1309:       DMSetDimension(idm, dim);
1310:       if (depth > 0) {
1311:         DMPlexInterpolateFaces_Internal(odm, 1, idm);
1312:         DMGetPointSF(odm, &sfPoint);
1313:         {
1314:           /* TODO: We need to systematically fix cases of distributed Plexes with no graph set */
1315:           PetscInt nroots;
1316:           PetscSFGetGraph(sfPoint, &nroots, NULL, NULL, NULL);
1317:           if (nroots >= 0) DMPlexInterpolatePointSF(idm, sfPoint);
1318:         }
1319:       }
1320:       if (odm != dm) DMDestroy(&odm);
1321:       odm = idm;
1322:     }
1323:     PetscObjectGetName((PetscObject) dm,  &name);
1324:     PetscObjectSetName((PetscObject) idm,  name);
1325:     DMPlexCopyCoordinates(dm, idm);
1326:     DMCopyLabels(dm, idm, PETSC_COPY_VALUES, PETSC_FALSE, DM_COPY_LABELS_FAIL);
1327:     PetscOptionsGetBool(((PetscObject)dm)->options, ((PetscObject)dm)->prefix, "-dm_plex_interpolate_orient_interfaces", &flg, NULL);
1328:     if (flg) DMPlexOrientInterface_Internal(idm);
1329:   }
1330:   /* This function makes the mesh fully interpolated on all ranks */
1331:   {
1332:     DM_Plex *plex = (DM_Plex *) idm->data;
1333:     plex->interpolated = plex->interpolatedCollective = DMPLEX_INTERPOLATED_FULL;
1334:   }
1335:   DMPlexCopy_Internal(dm, PETSC_TRUE, idm);
1336:   *dmInt = idm;
1337:   PetscLogEventEnd(DMPLEX_Interpolate,dm,0,0,0);
1338:   return 0;
1339: }

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

1344:   Collective on dmA

1346:   Input Parameter:
1347: . dmA - The DMPlex object with initial coordinates

1349:   Output Parameter:
1350: . dmB - The DMPlex object with copied coordinates

1352:   Level: intermediate

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

1356: .seealso: DMCopyLabels(), DMGetCoordinates(), DMGetCoordinatesLocal(), DMGetCoordinateDM(), DMGetCoordinateSection()
1357: @*/
1358: PetscErrorCode DMPlexCopyCoordinates(DM dmA, DM dmB)
1359: {
1360:   Vec            coordinatesA, coordinatesB;
1361:   VecType        vtype;
1362:   PetscSection   coordSectionA, coordSectionB;
1363:   PetscScalar   *coordsA, *coordsB;
1364:   PetscInt       spaceDim, Nf, vStartA, vStartB, vEndA, vEndB, coordSizeB, v, d;
1365:   PetscInt       cStartA, cEndA, cStartB, cEndB, cS, cE, cdim;
1366:   PetscBool      lc = PETSC_FALSE;

1370:   if (dmA == dmB) return 0;
1371:   DMGetCoordinateDim(dmA, &cdim);
1372:   DMSetCoordinateDim(dmB, cdim);
1373:   DMPlexGetDepthStratum(dmA, 0, &vStartA, &vEndA);
1374:   DMPlexGetDepthStratum(dmB, 0, &vStartB, &vEndB);
1376:   /* Copy over discretization if it exists */
1377:   {
1378:     DM                 cdmA, cdmB;
1379:     PetscDS            dsA, dsB;
1380:     PetscObject        objA, objB;
1381:     PetscClassId       idA, idB;
1382:     const PetscScalar *constants;
1383:     PetscInt            cdim, Nc;

1385:     DMGetCoordinateDM(dmA, &cdmA);
1386:     DMGetCoordinateDM(dmB, &cdmB);
1387:     DMGetField(cdmA, 0, NULL, &objA);
1388:     DMGetField(cdmB, 0, NULL, &objB);
1389:     PetscObjectGetClassId(objA, &idA);
1390:     PetscObjectGetClassId(objB, &idB);
1391:     if ((idA == PETSCFE_CLASSID) && (idA != idB)) {
1392:       DMSetField(cdmB, 0, NULL, objA);
1393:       DMCreateDS(cdmB);
1394:       DMGetDS(cdmA, &dsA);
1395:       DMGetDS(cdmB, &dsB);
1396:       PetscDSGetCoordinateDimension(dsA, &cdim);
1397:       PetscDSSetCoordinateDimension(dsB, cdim);
1398:       PetscDSGetConstants(dsA, &Nc, &constants);
1399:       PetscDSSetConstants(dsB, Nc, (PetscScalar *) constants);
1400:     }
1401:   }
1402:   DMPlexGetHeightStratum(dmA, 0, &cStartA, &cEndA);
1403:   DMPlexGetHeightStratum(dmB, 0, &cStartB, &cEndB);
1404:   DMGetCoordinateSection(dmA, &coordSectionA);
1405:   DMGetCoordinateSection(dmB, &coordSectionB);
1406:   if (coordSectionA == coordSectionB) return 0;
1407:   PetscSectionGetNumFields(coordSectionA, &Nf);
1408:   if (!Nf) return 0;
1410:   if (!coordSectionB) {
1411:     PetscInt dim;

1413:     PetscSectionCreate(PetscObjectComm((PetscObject) coordSectionA), &coordSectionB);
1414:     DMGetCoordinateDim(dmA, &dim);
1415:     DMSetCoordinateSection(dmB, dim, coordSectionB);
1416:     PetscObjectDereference((PetscObject) coordSectionB);
1417:   }
1418:   PetscSectionSetNumFields(coordSectionB, 1);
1419:   PetscSectionGetFieldComponents(coordSectionA, 0, &spaceDim);
1420:   PetscSectionSetFieldComponents(coordSectionB, 0, spaceDim);
1421:   PetscSectionGetChart(coordSectionA, &cS, &cE);
1422:   if (cStartA <= cS && cS < cEndA) { /* localized coordinates */
1424:     cS = cS - cStartA + cStartB;
1425:     cE = vEndB;
1426:     lc = PETSC_TRUE;
1427:   } else {
1428:     cS = vStartB;
1429:     cE = vEndB;
1430:   }
1431:   PetscSectionSetChart(coordSectionB, cS, cE);
1432:   for (v = vStartB; v < vEndB; ++v) {
1433:     PetscSectionSetDof(coordSectionB, v, spaceDim);
1434:     PetscSectionSetFieldDof(coordSectionB, v, 0, spaceDim);
1435:   }
1436:   if (lc) { /* localized coordinates */
1437:     PetscInt c;

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

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

1462:     PetscSectionGetOffset(coordSectionA, v + vStartA, &offA);
1463:     PetscSectionGetOffset(coordSectionB, v + vStartB, &offB);
1464:     for (d = 0; d < spaceDim; ++d) {
1465:       coordsB[offB+d] = coordsA[offA+d];
1466:     }
1467:   }
1468:   if (lc) { /* localized coordinates */
1469:     PetscInt c;

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

1474:       PetscSectionGetOffset(coordSectionA, c + cStartA, &offA);
1475:       PetscSectionGetOffset(coordSectionB, c + cStartB, &offB);
1476:       PetscSectionGetDof(coordSectionA, c + cStartA, &dof);
1477:       PetscArraycpy(coordsB + offB,coordsA + offA,dof);
1478:     }
1479:   }
1480:   VecRestoreArray(coordinatesA, &coordsA);
1481:   VecRestoreArray(coordinatesB, &coordsB);
1482:   DMSetCoordinatesLocal(dmB, coordinatesB);
1483:   VecDestroy(&coordinatesB);
1484:   return 0;
1485: }

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

1490:   Collective on dm

1492:   Input Parameter:
1493: . dm - The complete DMPlex object

1495:   Output Parameter:
1496: . dmUnint - The DMPlex object with only cells and vertices

1498:   Level: intermediate

1500:   Notes:
1501:     It does not copy over the coordinates.

1503:   Developer Notes:
1504:     It sets plex->interpolated = DMPLEX_INTERPOLATED_NONE.

1506: .seealso: DMPlexInterpolate(), DMPlexCreateFromCellListPetsc(), DMPlexCopyCoordinates()
1507: @*/
1508: PetscErrorCode DMPlexUninterpolate(DM dm, DM *dmUnint)
1509: {
1510:   DMPlexInterpolatedFlag interpolated;
1511:   DM             udm;
1512:   PetscInt       dim, vStart, vEnd, cStart, cEnd, c, maxConeSize = 0, *cone;

1516:   DMGetDimension(dm, &dim);
1517:   DMPlexIsInterpolated(dm, &interpolated);
1519:   if (interpolated == DMPLEX_INTERPOLATED_NONE || dim <= 1) {
1520:     /* in case dim <= 1 just keep the DMPLEX_INTERPOLATED_FULL flag */
1521:     PetscObjectReference((PetscObject) dm);
1522:     *dmUnint = dm;
1523:     return 0;
1524:   }
1525:   DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
1526:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
1527:   DMCreate(PetscObjectComm((PetscObject) dm), &udm);
1528:   DMSetType(udm, DMPLEX);
1529:   DMSetDimension(udm, dim);
1530:   DMPlexSetChart(udm, cStart, vEnd);
1531:   for (c = cStart; c < cEnd; ++c) {
1532:     PetscInt *closure = NULL, closureSize, cl, coneSize = 0;

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

1538:       if ((p >= vStart) && (p < vEnd)) ++coneSize;
1539:     }
1540:     DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
1541:     DMPlexSetConeSize(udm, c, coneSize);
1542:     maxConeSize = PetscMax(maxConeSize, coneSize);
1543:   }
1544:   DMSetUp(udm);
1545:   PetscMalloc1(maxConeSize, &cone);
1546:   for (c = cStart; c < cEnd; ++c) {
1547:     PetscInt *closure = NULL, closureSize, cl, coneSize = 0;

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

1553:       if ((p >= vStart) && (p < vEnd)) cone[coneSize++] = p;
1554:     }
1555:     DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
1556:     DMPlexSetCone(udm, c, cone);
1557:   }
1558:   PetscFree(cone);
1559:   DMPlexSymmetrize(udm);
1560:   DMPlexStratify(udm);
1561:   /* Reduce SF */
1562:   {
1563:     PetscSF            sfPoint, sfPointUn;
1564:     const PetscSFNode *remotePoints;
1565:     const PetscInt    *localPoints;
1566:     PetscSFNode       *remotePointsUn;
1567:     PetscInt          *localPointsUn;
1568:     PetscInt           vEnd, numRoots, numLeaves, l;
1569:     PetscInt           numLeavesUn = 0, n = 0;

1571:     /* Get original SF information */
1572:     DMGetPointSF(dm, &sfPoint);
1573:     DMGetPointSF(udm, &sfPointUn);
1574:     DMPlexGetDepthStratum(dm, 0, NULL, &vEnd);
1575:     PetscSFGetGraph(sfPoint, &numRoots, &numLeaves, &localPoints, &remotePoints);
1576:     /* Allocate space for cells and vertices */
1577:     for (l = 0; l < numLeaves; ++l) if (localPoints[l] < vEnd) numLeavesUn++;
1578:     /* Fill in leaves */
1579:     if (vEnd >= 0) {
1580:       PetscMalloc1(numLeavesUn, &remotePointsUn);
1581:       PetscMalloc1(numLeavesUn, &localPointsUn);
1582:       for (l = 0; l < numLeaves; l++) {
1583:         if (localPoints[l] < vEnd) {
1584:           localPointsUn[n]        = localPoints[l];
1585:           remotePointsUn[n].rank  = remotePoints[l].rank;
1586:           remotePointsUn[n].index = remotePoints[l].index;
1587:           ++n;
1588:         }
1589:       }
1591:       PetscSFSetGraph(sfPointUn, vEnd, numLeavesUn, localPointsUn, PETSC_OWN_POINTER, remotePointsUn, PETSC_OWN_POINTER);
1592:     }
1593:   }
1594:   /* This function makes the mesh fully uninterpolated on all ranks */
1595:   {
1596:     DM_Plex *plex = (DM_Plex *) udm->data;
1597:     plex->interpolated = plex->interpolatedCollective = DMPLEX_INTERPOLATED_NONE;
1598:   }
1599:   DMPlexCopy_Internal(dm, PETSC_TRUE, udm);
1600:   *dmUnint = udm;
1601:   return 0;
1602: }

1604: static PetscErrorCode DMPlexIsInterpolated_Internal(DM dm, DMPlexInterpolatedFlag *interpolated)
1605: {
1606:   PetscInt       coneSize, depth, dim, h, p, pStart, pEnd;
1607:   MPI_Comm       comm;

1609:   PetscObjectGetComm((PetscObject)dm, &comm);
1610:   DMPlexGetDepth(dm, &depth);
1611:   DMGetDimension(dm, &dim);

1613:   if (depth == dim) {
1614:     *interpolated = DMPLEX_INTERPOLATED_FULL;
1615:     if (!dim) goto finish;

1617:     /* Check points at height = dim are vertices (have no cones) */
1618:     DMPlexGetHeightStratum(dm, dim, &pStart, &pEnd);
1619:     for (p=pStart; p<pEnd; p++) {
1620:       DMPlexGetConeSize(dm, p, &coneSize);
1621:       if (coneSize) {
1622:         *interpolated = DMPLEX_INTERPOLATED_PARTIAL;
1623:         goto finish;
1624:       }
1625:     }

1627:     /* Check points at height < dim have cones */
1628:     for (h=0; h<dim; h++) {
1629:       DMPlexGetHeightStratum(dm, h, &pStart, &pEnd);
1630:       for (p=pStart; p<pEnd; p++) {
1631:         DMPlexGetConeSize(dm, p, &coneSize);
1632:         if (!coneSize) {
1633:           *interpolated = DMPLEX_INTERPOLATED_PARTIAL;
1634:           goto finish;
1635:         }
1636:       }
1637:     }
1638:   } else if (depth == 1) {
1639:     *interpolated = DMPLEX_INTERPOLATED_NONE;
1640:   } else {
1641:     *interpolated = DMPLEX_INTERPOLATED_PARTIAL;
1642:   }
1643: finish:
1644:   return 0;
1645: }

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

1650:   Not Collective

1652:   Input Parameter:
1653: . dm      - The DM object

1655:   Output Parameter:
1656: . interpolated - Flag whether the DM is interpolated

1658:   Level: intermediate

1660:   Notes:
1661:   Unlike DMPlexIsInterpolatedCollective(), this is NOT collective
1662:   so the results can be different on different ranks in special cases.
1663:   However, DMPlexInterpolate() guarantees the result is the same on all.

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

1667:   Developer Notes:
1668:   Initially, plex->interpolated = DMPLEX_INTERPOLATED_INVALID.

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

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

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

1679: .seealso: DMPlexInterpolate(), DMPlexIsInterpolatedCollective()
1680: @*/
1681: PetscErrorCode DMPlexIsInterpolated(DM dm, DMPlexInterpolatedFlag *interpolated)
1682: {
1683:   DM_Plex        *plex = (DM_Plex *) dm->data;

1687:   if (plex->interpolated < 0) {
1688:     DMPlexIsInterpolated_Internal(dm, &plex->interpolated);
1689:   } else if (PetscDefined (USE_DEBUG)) {
1690:     DMPlexInterpolatedFlag flg;

1692:     DMPlexIsInterpolated_Internal(dm, &flg);
1694:   }
1695:   *interpolated = plex->interpolated;
1696:   return 0;
1697: }

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

1702:   Collective

1704:   Input Parameter:
1705: . dm      - The DM object

1707:   Output Parameter:
1708: . interpolated - Flag whether the DM is interpolated

1710:   Level: intermediate

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

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

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

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

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

1727: .seealso: DMPlexInterpolate(), DMPlexIsInterpolated()
1728: @*/
1729: PetscErrorCode DMPlexIsInterpolatedCollective(DM dm, DMPlexInterpolatedFlag *interpolated)
1730: {
1731:   DM_Plex        *plex = (DM_Plex *) dm->data;
1732:   PetscBool       debug=PETSC_FALSE;

1736:   PetscOptionsGetBool(((PetscObject) dm)->options, ((PetscObject) dm)->prefix, "-dm_plex_is_interpolated_collective_debug", &debug, NULL);
1737:   if (plex->interpolatedCollective < 0) {
1738:     DMPlexInterpolatedFlag  min, max;
1739:     MPI_Comm                comm;

1741:     PetscObjectGetComm((PetscObject)dm, &comm);
1742:     DMPlexIsInterpolated(dm, &plex->interpolatedCollective);
1743:     MPI_Allreduce(&plex->interpolatedCollective, &min, 1, MPIU_ENUM, MPI_MIN, comm);
1744:     MPI_Allreduce(&plex->interpolatedCollective, &max, 1, MPIU_ENUM, MPI_MAX, comm);
1745:     if (min != max) plex->interpolatedCollective = DMPLEX_INTERPOLATED_MIXED;
1746:     if (debug) {
1747:       PetscMPIInt rank;

1749:       MPI_Comm_rank(comm, &rank);
1750:       PetscSynchronizedPrintf(comm, "[%d] interpolated=%s interpolatedCollective=%s\n", rank, DMPlexInterpolatedFlags[plex->interpolated], DMPlexInterpolatedFlags[plex->interpolatedCollective]);
1751:       PetscSynchronizedFlush(comm, PETSC_STDOUT);
1752:     }
1753:   }
1754:   *interpolated = plex->interpolatedCollective;
1755:   return 0;
1756: }