Actual source code: plexinterpolate.c

petsc-master 2020-02-19
Report Typos and Errors
  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_", 0};

  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: 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: 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:   DMPlexGetFaces_Internal - Gets groups of vertices that correspond to faces for the given cell
 55:   This assumes that the mesh is not interpolated from the depth of point p to the vertices
 56: */
 57: PetscErrorCode DMPlexGetFaces_Internal(DM dm, PetscInt p, PetscInt *numFaces, PetscInt *faceSize, const PetscInt *faces[])
 58: {
 59:   const PetscInt *cone = NULL;
 60:   DMPolytopeType  ct;
 61:   PetscErrorCode  ierr;

 65:   DMPlexGetCellType(dm, p, &ct);
 66:   DMPlexGetCone(dm, p, &cone);
 67:   DMPlexGetRawFaces_Internal(dm, ct, cone, numFaces, faceSize, faces);
 68:   return(0);
 69: }

 71: /*
 72:   DMPlexRestoreFaces_Internal - Restores the array
 73: */
 74: PetscErrorCode DMPlexRestoreFaces_Internal(DM dm, PetscInt p, PetscInt *numFaces, PetscInt *faceSize, const PetscInt *faces[])
 75: {
 76:   PetscErrorCode  ierr;

 79:   if (faces) { DMRestoreWorkArray(dm, 0, MPIU_INT, (void *) faces); }
 80:   return(0);
 81: }

 83: /*
 84:   DMPlexGetRawFaces_Internal - Gets groups of vertices that correspond to faces for the given cone
 85: */
 86: PetscErrorCode DMPlexGetRawFaces_Internal(DM dm, DMPolytopeType ct, const PetscInt cone[], PetscInt *numFaces, PetscInt *faceSize, const PetscInt *faces[])
 87: {
 88:   PetscInt       *facesTmp;
 89:   PetscInt        maxConeSize, maxSupportSize;
 90:   PetscErrorCode  ierr;

 95:   DMPlexGetMaxSizes(dm, &maxConeSize, &maxSupportSize);
 96:   if (faces) {DMGetWorkArray(dm, PetscSqr(PetscMax(maxConeSize, maxSupportSize)), MPIU_INT, &facesTmp);}
 97:   switch (ct) {
 98:     case DM_POLYTOPE_POINT:
 99:       if (numFaces) *numFaces = 0;
100:       if (faceSize) *faceSize = 0;
101:       break;
102:     case DM_POLYTOPE_SEGMENT:
103:       if (faces) {
104:         facesTmp[0] = cone[0]; facesTmp[1] = cone[1];
105:         *faces = facesTmp;
106:       }
107:       if (numFaces) *numFaces = 2;
108:       if (faceSize) *faceSize = 1;
109:       break;
110:     case DM_POLYTOPE_TRIANGLE:
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:       if (numFaces) *numFaces = 3;
118:       if (faceSize) *faceSize = 2;
119:       break;
120:     case DM_POLYTOPE_QUADRILATERAL:
121:       /* Vertices follow right hand rule */
122:       if (faces) {
123:         facesTmp[0] = cone[0]; facesTmp[1] = cone[1];
124:         facesTmp[2] = cone[1]; facesTmp[3] = cone[2];
125:         facesTmp[4] = cone[2]; facesTmp[5] = cone[3];
126:         facesTmp[6] = cone[3]; facesTmp[7] = cone[0];
127:         *faces = facesTmp;
128:       }
129:       if (numFaces) *numFaces = 4;
130:       if (faceSize) *faceSize = 2;
131:       break;
132:     case DM_POLYTOPE_TETRAHEDRON:
133:       /* Vertices of first face follow right hand rule and normal points away from last vertex */
134:       if (faces) {
135:         facesTmp[0] = cone[0]; facesTmp[1]  = cone[1]; facesTmp[2]  = cone[2];
136:         facesTmp[3] = cone[0]; facesTmp[4]  = cone[3]; facesTmp[5]  = cone[1];
137:         facesTmp[6] = cone[0]; facesTmp[7]  = cone[2]; facesTmp[8]  = cone[3];
138:         facesTmp[9] = cone[2]; facesTmp[10] = cone[1]; facesTmp[11] = cone[3];
139:         *faces = facesTmp;
140:       }
141:       if (numFaces) *numFaces = 4;
142:       if (faceSize) *faceSize = 3;
143:       break;
144:     case DM_POLYTOPE_HEXAHEDRON:
145:       /*  7--------6
146:          /|       /|
147:         / |      / |
148:        4--------5  |
149:        |  |     |  |
150:        |  |     |  |
151:        |  1--------2
152:        | /      | /
153:        |/       |/
154:        0--------3
155:        */
156:       if (faces) {
157:         facesTmp[0]  = cone[0]; facesTmp[1]  = cone[1]; facesTmp[2]  = cone[2]; facesTmp[3]  = cone[3]; /* Bottom */
158:         facesTmp[4]  = cone[4]; facesTmp[5]  = cone[5]; facesTmp[6]  = cone[6]; facesTmp[7]  = cone[7]; /* Top */
159:         facesTmp[8]  = cone[0]; facesTmp[9]  = cone[3]; facesTmp[10] = cone[5]; facesTmp[11] = cone[4]; /* Front */
160:         facesTmp[12] = cone[2]; facesTmp[13] = cone[1]; facesTmp[14] = cone[7]; facesTmp[15] = cone[6]; /* Back */
161:         facesTmp[16] = cone[3]; facesTmp[17] = cone[2]; facesTmp[18] = cone[6]; facesTmp[19] = cone[5]; /* Right */
162:         facesTmp[20] = cone[0]; facesTmp[21] = cone[4]; facesTmp[22] = cone[7]; facesTmp[23] = cone[1]; /* Left */
163:         *faces = facesTmp;
164:       }
165:       if (numFaces) *numFaces = 6;
166:       if (faceSize) *faceSize = 4;
167:       break;
168:     default: SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_OUTOFRANGE, "No face description for cell type %s", DMPolytopeTypes[ct]);
169:   }
170:   return(0);
171: }

173: /*
174:   DMPlexGetRawFacesHybrid_Internal - Gets groups of vertices that correspond to faces for the given cone using hybrid ordering (prisms)
175: */
176: static PetscErrorCode DMPlexGetRawFacesHybrid_Internal(DM dm, DMPolytopeType ct, const PetscInt cone[], PetscInt *numFaces, PetscInt *numFacesNotH, PetscInt *faceSize, const PetscInt *faces[])
177: {
178:   PetscInt       *facesTmp;
179:   PetscInt        maxConeSize, maxSupportSize;
180:   PetscErrorCode  ierr;

185:   DMPlexGetMaxSizes(dm, &maxConeSize, &maxSupportSize);
186:   if (faces) {DMGetWorkArray(dm, PetscSqr(PetscMax(maxConeSize, maxSupportSize)), MPIU_INT, &facesTmp);}
187:   switch (ct) {
188:     case DM_POLYTOPE_SEGMENT:
189:       if (faces) {
190:         facesTmp[0] = cone[0]; facesTmp[1] = cone[1];
191:         *faces = facesTmp;
192:       }
193:       if (numFaces)     *numFaces = 2;
194:       if (numFacesNotH) *numFacesNotH = 2;
195:       if (faceSize)     *faceSize = 1;
196:       break;
197:     case DM_POLYTOPE_QUADRILATERAL:
198:     case DM_POLYTOPE_SEG_PRISM_TENSOR:
199:       if (faces) {
200:         facesTmp[0] = cone[0]; facesTmp[1] = cone[1];
201:         facesTmp[2] = cone[2]; facesTmp[3] = cone[3];
202:         facesTmp[4] = cone[0]; facesTmp[5] = cone[2];
203:         facesTmp[6] = cone[1]; facesTmp[7] = cone[3];
204:         *faces = facesTmp;
205:       }
206:       if (numFaces)     *numFaces = 4;
207:       if (numFacesNotH) *numFacesNotH = 2;
208:       if (faceSize)     *faceSize = 2;
209:       break;
210:     case DM_POLYTOPE_TRI_PRISM:
211:       if (faces) {
212:         facesTmp[0]  = cone[0]; facesTmp[1]  = cone[2]; facesTmp[2]  = cone[1]; facesTmp[3]  = -1;      /* Bottom */
213:         facesTmp[4]  = cone[3]; facesTmp[5]  = cone[4]; facesTmp[6]  = cone[5]; facesTmp[7]  = -1;      /* Top */
214:         facesTmp[8]  = cone[0]; facesTmp[9]  = cone[1]; facesTmp[10] = cone[4]; facesTmp[11] = cone[3]; /* Back left */
215:         facesTmp[12] = cone[1]; facesTmp[13] = cone[2]; facesTmp[14] = cone[5]; facesTmp[15] = cone[4]; /* Back right */
216:         facesTmp[16] = cone[2]; facesTmp[17] = cone[0]; facesTmp[18] = cone[3]; facesTmp[19] = cone[5]; /* Front */
217:         *faces = facesTmp;
218:       }
219:       if (numFaces)     *numFaces = 5;
220:       if (numFacesNotH) *numFacesNotH = 2;
221:       if (faceSize)     *faceSize = -4;
222:       break;
223:     case DM_POLYTOPE_TRI_PRISM_TENSOR:
224:       if (faces) {
225:         facesTmp[0]  = cone[0]; facesTmp[1]  = cone[1]; facesTmp[2]  = cone[2]; facesTmp[3]  = -1;      /* Bottom */
226:         facesTmp[4]  = cone[3]; facesTmp[5]  = cone[4]; facesTmp[6]  = cone[5]; facesTmp[7]  = -1;      /* Top */
227:         facesTmp[8]  = cone[0]; facesTmp[9]  = cone[1]; facesTmp[10] = cone[3]; facesTmp[11] = cone[4]; /* Back left */
228:         facesTmp[12] = cone[1]; facesTmp[13] = cone[2]; facesTmp[14] = cone[4]; facesTmp[15] = cone[5]; /* Back right */
229:         facesTmp[16] = cone[2]; facesTmp[17] = cone[0]; facesTmp[18] = cone[5]; facesTmp[19] = cone[3]; /* Front */
230:         *faces = facesTmp;
231:       }
232:       if (numFaces)     *numFaces = 5;
233:       if (numFacesNotH) *numFacesNotH = 2;
234:       if (faceSize)     *faceSize = -4;
235:       break;
236:     default: SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_OUTOFRANGE, "No face description for cell type %s", DMPolytopeTypes[ct]);
237:   }
238:   return(0);
239: }

241: static PetscErrorCode DMPlexRestoreRawFacesHybrid_Internal(DM dm, DMPolytopeType ct, const PetscInt cone[], PetscInt *numFaces, PetscInt *numFacesNotH, PetscInt *faceSize, const PetscInt *faces[])
242: {
243:   PetscErrorCode  ierr;

246:   if (faces) { DMRestoreWorkArray(dm, 0, MPIU_INT, (void *) faces); }
247:   return(0);
248: }

250: static PetscErrorCode DMPlexGetFacesHybrid_Internal(DM dm, PetscInt p, PetscInt *numFaces, PetscInt *numFacesNotH, PetscInt *faceSize, const PetscInt *faces[])
251: {
252:   const PetscInt *cone = NULL;
253:   DMPolytopeType  ct;
254:   PetscErrorCode  ierr;

258:   DMPlexGetCellType(dm, p, &ct);
259:   DMPlexGetCone(dm, p, &cone);
260:   DMPlexGetRawFacesHybrid_Internal(dm, ct, cone, numFaces, numFacesNotH, faceSize, faces);
261:   return(0);
262: }

264: /* This interpolates faces for cells at some stratum */
265: static PetscErrorCode DMPlexInterpolateFaces_Internal(DM dm, PetscInt cellDepth, DM idm)
266: {
267:   DMLabel        subpointMap;
268:   PetscHashIJKL  faceTable;
269:   PetscInt      *pStart, *pEnd;
270:   PetscInt       cellDim, depth, faceDepth = cellDepth, numPoints = 0, faceSizeAll = 0, face, c, d;
271:   PetscInt       coneSizeH = 0, faceSizeAllH = 0, faceSizeAllT = 0, numCellFacesH = 0, faceT = 0, faceH, pMax = -1, dim, outerloop;
272:   PetscInt       cMax, fMax, eMax, vMax;

276:   DMGetDimension(dm, &cellDim);
277:   /* HACK: I need a better way to determine face dimension, or an alternative to GetFaces() */
278:   DMPlexGetSubpointMap(dm, &subpointMap);
279:   if (subpointMap) ++cellDim;
280:   DMPlexGetDepth(dm, &depth);
281:   ++depth;
282:   ++cellDepth;
283:   cellDim -= depth - cellDepth;
284:   PetscMalloc2(depth+1,&pStart,depth+1,&pEnd);
285:   for (d = depth-1; d >= faceDepth; --d) {
286:     DMPlexGetDepthStratum(dm, d, &pStart[d+1], &pEnd[d+1]);
287:   }
288:   DMPlexGetDepthStratum(dm, -1, NULL, &pStart[faceDepth]);
289:   pEnd[faceDepth] = pStart[faceDepth];
290:   for (d = faceDepth-1; d >= 0; --d) {
291:     DMPlexGetDepthStratum(dm, d, &pStart[d], &pEnd[d]);
292:   }
293:   cMax = fMax = eMax = vMax = PETSC_DETERMINE;
294:   DMGetDimension(dm, &dim);
295:   if (cellDim == dim) {
296:     DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);
297:     pMax = cMax;
298:   } else if (cellDim == dim -1) {
299:     DMPlexGetHybridBounds(dm, &cMax, &fMax, NULL, NULL);
300:     pMax = fMax;
301:   }
302:   pMax = pMax < 0 ? pEnd[cellDepth] : pMax;
303:   if (pMax < pEnd[cellDepth]) {
304:     const PetscInt *cellFaces, *cone;
305:     DMPolytopeType  ct;
306:     PetscInt        numCellFacesT, faceSize, cf;

308:     /* First get normal cell face size (we now allow hybrid cells to meet normal cells on either hybrid or normal faces */
309:     if (pStart[cellDepth] < pMax) {DMPlexGetFaces_Internal(dm, pStart[cellDepth], NULL, &faceSizeAll, NULL);}

311:     DMPlexGetCellType(dm, pMax, &ct);
312:     DMPlexGetConeSize(dm, pMax, &coneSizeH);
313:     DMPlexGetCone(dm, pMax, &cone);
314:     DMPlexGetRawFacesHybrid_Internal(dm, ct, cone, &numCellFacesH, &numCellFacesT, &faceSize, &cellFaces);
315:     if (faceSize < 0) {
316:       PetscInt *sizes, minv, maxv;

318:       /* count vertices of hybrid and non-hybrid faces */
319:       PetscCalloc1(numCellFacesH, &sizes);
320:       for (cf = 0; cf < numCellFacesT; ++cf) { /* These are the non-hybrid faces */
321:         const PetscInt *cellFace = &cellFaces[-cf*faceSize];
322:         PetscInt       f;

324:         for (f = 0; f < -faceSize; ++f) sizes[cf] += (cellFace[f] >= 0 ? 1 : 0);
325:       }
326:       PetscSortInt(numCellFacesT, sizes);
327:       minv = sizes[0];
328:       maxv = sizes[PetscMax(numCellFacesT-1, 0)];
329:       if (minv != maxv) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_SUP, "Different number of vertices for non-hybrid face %D != %D", minv, maxv);
330:       faceSizeAllT = minv;
331:       PetscArrayzero(sizes, numCellFacesH);
332:       for (cf = numCellFacesT; cf < numCellFacesH; ++cf) { /* These are the hybrid faces */
333:         const PetscInt *cellFace = &cellFaces[-cf*faceSize];
334:         PetscInt       f;

336:         for (f = 0; f < -faceSize; ++f) sizes[cf-numCellFacesT] += (cellFace[f] >= 0 ? 1 : 0);
337:       }
338:       PetscSortInt(numCellFacesH - numCellFacesT, sizes);
339:       minv = sizes[0];
340:       maxv = sizes[PetscMax(numCellFacesH - numCellFacesT-1, 0)];
341:       PetscFree(sizes);
342:       if (minv != maxv) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_SUP, "Different number of vertices for hybrid face %D != %D", minv, maxv);
343:       faceSizeAllH = minv;
344:       if (!faceSizeAll) faceSizeAll = faceSizeAllT;
345:     } else { /* the size of the faces in hybrid cells is the same */
346:       faceSizeAll = faceSizeAllH = faceSizeAllT = faceSize;
347:     }
348:     DMPlexRestoreRawFacesHybrid_Internal(dm, ct, cone, &numCellFacesH, &numCellFacesT, &faceSize, &cellFaces);
349:   } else if (pEnd[cellDepth] > pStart[cellDepth]) {
350:     DMPlexGetFaces_Internal(dm, pStart[cellDepth], NULL, &faceSizeAll, NULL);
351:     faceSizeAllH = faceSizeAllT = faceSizeAll;
352:   }
353:   if (faceSizeAll > 4) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Do not support interpolation of meshes with faces of %D vertices", faceSizeAll);

355:   /* With hybrid grids, we first iterate on hybrid cells and start numbering the non-hybrid faces
356:      Then, faces for non-hybrid cells are numbered.
357:      This is to guarantee consistent orientations (all 0) of all the points in the cone of the hybrid cells */
358:   PetscHashIJKLCreate(&faceTable);
359:   for (outerloop = 0, face = pStart[faceDepth]; outerloop < 2; outerloop++) {
360:     PetscInt start, end;

362:     start = outerloop == 0 ? pMax : pStart[cellDepth];
363:     end = outerloop == 0 ? pEnd[cellDepth] : pMax;
364:     for (c = start; c < end; ++c) {
365:       const PetscInt *cellFaces;
366:       PetscInt        numCellFaces, faceSize, faceSizeInc, faceSizeCheck, cf;

368:       if (c < pMax) {
369:         DMPlexGetFaces_Internal(dm, c, &numCellFaces, &faceSize, &cellFaces);
370:         if (faceSize != faceSizeAll) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Inconsistent face for cell %D of size %D != %D", c, faceSize, faceSizeAll);
371:         faceSizeCheck = faceSizeAll;
372:       } else { /* Hybrid cell */
373:         const PetscInt *cone;
374:         DMPolytopeType  ct;
375:         PetscInt        numCellFacesN, coneSize;

377:         DMPlexGetCellType(dm, c, &ct);
378:         DMPlexGetConeSize(dm, c, &coneSize);
379:         if (coneSize != coneSizeH) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unexpected hybrid coneSize %D != %D", coneSize, coneSizeH);
380:         DMPlexGetCone(dm, c, &cone);
381:         DMPlexGetRawFacesHybrid_Internal(dm, ct, cone, &numCellFaces, &numCellFacesN, &faceSize, &cellFaces);
382:         if (numCellFaces != numCellFacesH) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unexpected numCellFaces %D != %D for hybrid cell %D", numCellFaces, numCellFacesH, c);
383:         faceSize = PetscMax(faceSize, -faceSize);
384:         if (faceSize > 4) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Do not support interpolation of meshes with faces of %D vertices", faceSize);
385:         numCellFaces = numCellFacesN; /* process only non-hybrid faces */
386:         faceSizeCheck = faceSizeAllT;
387:       }
388:       faceSizeInc = faceSize;
389:       for (cf = 0; cf < numCellFaces; ++cf) {
390:         const PetscInt   *cellFace = &cellFaces[cf*faceSizeInc];
391:         PetscInt          faceSizeH = faceSize;
392:         PetscHashIJKLKey  key;
393:         PetscHashIter     iter;
394:         PetscBool         missing;

396:         if (faceSizeInc == 2) {
397:           key.i = PetscMin(cellFace[0], cellFace[1]);
398:           key.j = PetscMax(cellFace[0], cellFace[1]);
399:           key.k = PETSC_MAX_INT;
400:           key.l = PETSC_MAX_INT;
401:         } else {
402:           key.i = cellFace[0];
403:           key.j = cellFace[1];
404:           key.k = cellFace[2];
405:           key.l = faceSize > 3 ? (cellFace[3] < 0 ? faceSizeH = 3, PETSC_MAX_INT : cellFace[3]) : PETSC_MAX_INT;
406:           PetscSortInt(faceSize, (PetscInt *) &key);
407:         }
408:         /* this check is redundant for non-hybrid meshes */
409:         if (faceSizeH != faceSizeCheck) SETERRQ4(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unexpected number of vertices for face %D of point %D -> %D != %D", cf, c, faceSizeH, faceSizeCheck);
410:         PetscHashIJKLPut(faceTable, key, &iter, &missing);
411:         if (missing) {
412:           PetscHashIJKLIterSet(faceTable, iter, face++);
413:           if (c >= pMax) ++faceT;
414:         }
415:       }
416:       if (c < pMax) {
417:         DMPlexRestoreFaces_Internal(dm, c, &numCellFaces, &faceSize, &cellFaces);
418:       } else {
419:         DMPlexRestoreRawFacesHybrid_Internal(dm, DM_POLYTOPE_UNKNOWN, NULL, NULL, NULL, NULL, &cellFaces);
420:       }
421:     }
422:   }
423:   pEnd[faceDepth] = face;

425:   /* Second pass for hybrid meshes: number hybrid faces */
426:   for (c = pMax; c < pEnd[cellDepth]; ++c) {
427:     const PetscInt *cellFaces, *cone;
428:     DMPolytopeType  ct;
429:     PetscInt        numCellFaces, numCellFacesN, faceSize, cf;

431:     DMPlexGetCellType(dm, c, &ct);
432:     DMPlexGetCone(dm, c, &cone);
433:     DMPlexGetRawFacesHybrid_Internal(dm, ct, cone, &numCellFaces, &numCellFacesN, &faceSize, &cellFaces);
434:     if (numCellFaces != numCellFacesH) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unexpected hybrid numCellFaces %D != %D", numCellFaces, numCellFacesH);
435:     faceSize = PetscMax(faceSize, -faceSize);
436:     for (cf = numCellFacesN; cf < numCellFaces; ++cf) { /* These are the hybrid faces */
437:       const PetscInt   *cellFace = &cellFaces[cf*faceSize];
438:       PetscHashIJKLKey  key;
439:       PetscHashIter     iter;
440:       PetscBool         missing;
441:       PetscInt          faceSizeH = faceSize;

443:       if (faceSize == 2) {
444:         key.i = PetscMin(cellFace[0], cellFace[1]);
445:         key.j = PetscMax(cellFace[0], cellFace[1]);
446:         key.k = PETSC_MAX_INT;
447:         key.l = PETSC_MAX_INT;
448:       } else {
449:         key.i = cellFace[0];
450:         key.j = cellFace[1];
451:         key.k = cellFace[2];
452:         key.l = faceSize > 3 ? (cellFace[3] < 0 ? faceSizeH = 3, PETSC_MAX_INT : cellFace[3]) : PETSC_MAX_INT;
453:         PetscSortInt(faceSize, (PetscInt *) &key);
454:       }
455:       if (faceSizeH != faceSizeAllH) SETERRQ4(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unexpected number of vertices for hybrid face %D of point %D -> %D != %D", cf, c, faceSizeH, faceSizeAllH);
456:       PetscHashIJKLPut(faceTable, key, &iter, &missing);
457:       if (missing) {PetscHashIJKLIterSet(faceTable, iter, face++);}
458:     }
459:     DMPlexRestoreRawFacesHybrid_Internal(dm, ct, cone, &numCellFaces, &numCellFacesN, &faceSize, &cellFaces);
460:   }
461:   faceH = face - pEnd[faceDepth];
462:   if (faceH) {
463:     if (fMax == PETSC_DETERMINE) fMax = pEnd[faceDepth];
464:     else if (eMax == PETSC_DETERMINE) eMax = pEnd[faceDepth];
465:     else SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of unassigned hybrid facets %D for cellDim %D and dimension %D", faceH, cellDim, dim);
466:   }
467:   pEnd[faceDepth] = face;
468:   PetscHashIJKLDestroy(&faceTable);
469:   /* Count new points */
470:   for (d = 0; d <= depth; ++d) {
471:     numPoints += pEnd[d]-pStart[d];
472:   }
473:   DMPlexSetChart(idm, 0, numPoints);
474:   /* Set cone sizes */
475:   for (d = 0; d <= depth; ++d) {
476:     PetscInt coneSize, p;

478:     if (d == faceDepth) {
479:       /* Now we have two cases: */
480:       if (faceSizeAll == faceSizeAllT) {
481:         /* I see no way to do this if we admit faces of different shapes */
482:         for (p = pStart[d]; p < pEnd[d]-faceH; ++p) {
483:           DMPlexSetConeSize(idm, p, faceSizeAll);
484:         }
485:         for (p = pEnd[d]-faceH; p < pEnd[d]; ++p) {
486:           DMPlexSetConeSize(idm, p, faceSizeAllH);
487:         }
488:       } else if (faceSizeAll == faceSizeAllH) {
489:         for (p = pStart[d]; p < pStart[d]+faceT; ++p) {
490:           DMPlexSetConeSize(idm, p, faceSizeAllT);
491:         }
492:         for (p = pStart[d]+faceT; p < pEnd[d]-faceH; ++p) {
493:           DMPlexSetConeSize(idm, p, faceSizeAll);
494:         }
495:         for (p = pEnd[d]-faceH; p < pEnd[d]; ++p) {
496:           DMPlexSetConeSize(idm, p, faceSizeAllH);
497:         }
498:       } else SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Inconsistent faces sizes N: %D T: %D H: %D", faceSizeAll, faceSizeAllT, faceSizeAllH);
499:     } else if (d == cellDepth) {
500:       for (p = pStart[d]; p < pEnd[d]; ++p) {
501:         /* Number of cell faces may be different from number of cell vertices*/
502:         if (p < pMax) {
503:           DMPlexGetFaces_Internal(dm, p, &coneSize, NULL, NULL);
504:         } else {
505:           DMPlexGetFacesHybrid_Internal(dm, p, &coneSize, NULL, NULL, NULL);
506:         }
507:         DMPlexSetConeSize(idm, p, coneSize);
508:       }
509:     } else {
510:       for (p = pStart[d]; p < pEnd[d]; ++p) {
511:         DMPlexGetConeSize(dm, p, &coneSize);
512:         DMPlexSetConeSize(idm, p, coneSize);
513:       }
514:     }
515:   }
516:   DMSetUp(idm);
517:   /* Get face cones from subsets of cell vertices */
518:   if (faceSizeAll > 4) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Do not support interpolation of meshes with faces of %D vertices", faceSizeAll);
519:   PetscHashIJKLCreate(&faceTable);
520:   for (d = depth; d > cellDepth; --d) {
521:     const PetscInt *cone;
522:     PetscInt        p;

524:     for (p = pStart[d]; p < pEnd[d]; ++p) {
525:       DMPlexGetCone(dm, p, &cone);
526:       DMPlexSetCone(idm, p, cone);
527:       DMPlexGetConeOrientation(dm, p, &cone);
528:       DMPlexSetConeOrientation(idm, p, cone);
529:     }
530:   }
531:   for (outerloop = 0, face = pStart[faceDepth]; outerloop < 2; outerloop++) {
532:     PetscInt start, end;

534:     start = outerloop == 0 ? pMax : pStart[cellDepth];
535:     end = outerloop == 0 ? pEnd[cellDepth] : pMax;
536:     for (c = start; c < end; ++c) {
537:       const PetscInt *cellFaces;
538:       PetscInt        numCellFaces, faceSize, faceSizeInc, cf;

540:       if (c < pMax) {
541:         DMPlexGetFaces_Internal(dm, c, &numCellFaces, &faceSize, &cellFaces);
542:         if (faceSize != faceSizeAll) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Inconsistent face for cell %D of size %D != %D", c, faceSize, faceSizeAll);
543:       } else {
544:         const PetscInt *cone;
545:         DMPolytopeType  ct;
546:         PetscInt        numCellFacesN, coneSize;

548:         DMPlexGetCellType(dm, c, &ct);
549:         DMPlexGetConeSize(dm, c, &coneSize);
550:         if (coneSize != coneSizeH) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unexpected hybrid coneSize %D != %D", coneSize, coneSizeH);
551:         DMPlexGetCone(dm, c, &cone);
552:         DMPlexGetRawFacesHybrid_Internal(dm, ct, cone, &numCellFaces, &numCellFacesN, &faceSize, &cellFaces);
553:         if (numCellFaces != numCellFacesH) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unexpected numCellFaces %D != %D for hybrid cell %D", numCellFaces, numCellFacesH, c);
554:         faceSize = PetscMax(faceSize, -faceSize);
555:         if (faceSize > 4) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Do not support interpolation of meshes with faces of %D vertices", faceSize);
556:         numCellFaces = numCellFacesN; /* process only non-hybrid faces */
557:       }
558:       faceSizeInc = faceSize;
559:       for (cf = 0; cf < numCellFaces; ++cf) {
560:         const PetscInt  *cellFace = &cellFaces[cf*faceSizeInc];
561:         PetscHashIJKLKey key;
562:         PetscHashIter    iter;
563:         PetscBool        missing;

565:         if (faceSizeInc == 2) {
566:           key.i = PetscMin(cellFace[0], cellFace[1]);
567:           key.j = PetscMax(cellFace[0], cellFace[1]);
568:           key.k = PETSC_MAX_INT;
569:           key.l = PETSC_MAX_INT;
570:         } else {
571:           key.i = cellFace[0];
572:           key.j = cellFace[1];
573:           key.k = cellFace[2];
574:           key.l = faceSizeInc > 3 ? (cellFace[3] < 0 ? faceSize = 3, PETSC_MAX_INT : cellFace[3]) : PETSC_MAX_INT;
575:           PetscSortInt(faceSizeInc, (PetscInt *) &key);
576:         }
577:         PetscHashIJKLPut(faceTable, key, &iter, &missing);
578:         if (missing) {
579:           DMPlexSetCone(idm, face, cellFace);
580:           PetscHashIJKLIterSet(faceTable, iter, face);
581:           DMPlexInsertCone(idm, c, cf, face++);
582:         } else {
583:           const PetscInt *cone;
584:           PetscInt        coneSize, ornt, i, j, f;

586:           PetscHashIJKLIterGet(faceTable, iter, &f);
587:           DMPlexInsertCone(idm, c, cf, f);
588:           /* Orient face: Do not allow reverse orientation at the first vertex */
589:           DMPlexGetConeSize(idm, f, &coneSize);
590:           DMPlexGetCone(idm, f, &cone);
591:           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);
592:           /* - First find the initial vertex */
593:           for (i = 0; i < faceSize; ++i) if (cellFace[0] == cone[i]) break;
594:           /* - Try forward comparison */
595:           for (j = 0; j < faceSize; ++j) if (cellFace[j] != cone[(i+j)%faceSize]) break;
596:           if (j == faceSize) {
597:             if ((faceSize == 2) && (i == 1)) ornt = -2;
598:             else                             ornt = i;
599:           } else {
600:             /* - Try backward comparison */
601:             for (j = 0; j < faceSize; ++j) if (cellFace[j] != cone[(i+faceSize-j)%faceSize]) break;
602:             if (j == faceSize) {
603:               if (i == 0) ornt = -faceSize;
604:               else        ornt = -i;
605:             } else SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Could not determine orientation of face %D in cell %D", f, c);
606:           }
607:           DMPlexInsertConeOrientation(idm, c, cf, ornt);
608:         }
609:       }
610:       if (c < pMax) {
611:         DMPlexRestoreFaces_Internal(dm, c, &numCellFaces, &faceSize, &cellFaces);
612:       } else {
613:         DMPlexRestoreRawFacesHybrid_Internal(dm, DM_POLYTOPE_UNKNOWN, NULL, NULL, NULL, NULL, &cellFaces);
614:       }
615:     }
616:   }
617:   /* Second pass for hybrid meshes: orient hybrid faces */
618:   for (c = pMax; c < pEnd[cellDepth]; ++c) {
619:     const PetscInt *cellFaces, *cone;
620:     DMPolytopeType  ct;
621:     PetscInt        numCellFaces, numCellFacesN, faceSize, coneSize, cf;

623:     DMPlexGetCellType(dm, c, &ct);
624:     DMPlexGetConeSize(dm, c, &coneSize);
625:     DMPlexGetCone(dm, c, &cone);
626:     DMPlexGetRawFacesHybrid_Internal(dm, ct, cone, &numCellFaces, &numCellFacesN, &faceSize, &cellFaces);
627:     if (numCellFaces != numCellFacesH) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unexpected hybrid numCellFaces %D != %D", numCellFaces, numCellFacesH);
628:     faceSize = PetscMax(faceSize, -faceSize);
629:     for (cf = numCellFacesN; cf < numCellFaces; ++cf) { /* These are the hybrid faces */
630:       const PetscInt   *cellFace = &cellFaces[cf*faceSize];
631:       PetscHashIJKLKey key;
632:       PetscHashIter    iter;
633:       PetscBool        missing;
634:       PetscInt         faceSizeH = faceSize;

636:       if (faceSize == 2) {
637:         key.i = PetscMin(cellFace[0], cellFace[1]);
638:         key.j = PetscMax(cellFace[0], cellFace[1]);
639:         key.k = PETSC_MAX_INT;
640:         key.l = PETSC_MAX_INT;
641:       } else {
642:         key.i = cellFace[0];
643:         key.j = cellFace[1];
644:         key.k = cellFace[2];
645:         key.l = faceSize > 3 ? (cellFace[3] < 0 ? faceSizeH = 3, PETSC_MAX_INT : cellFace[3]) : PETSC_MAX_INT;
646:         PetscSortInt(faceSize, (PetscInt *) &key);
647:       }
648:       if (faceSizeH != faceSizeAllH) SETERRQ4(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unexpected number of vertices for hybrid face %D of point %D -> %D != %D", cf, c, faceSizeH, faceSizeAllH);
649:       PetscHashIJKLPut(faceTable, key, &iter, &missing);
650:       if (missing) {
651:         DMPlexSetCone(idm, face, cellFace);
652:         PetscHashIJKLIterSet(faceTable, iter, face);
653:         DMPlexInsertCone(idm, c, cf, face++);
654:       } else {
655:         PetscInt        fv[4] = {0, 1, 2, 3};
656:         const PetscInt *cone;
657:         PetscInt        coneSize, ornt, i, j, f;
658:         PetscBool       q2h = PETSC_FALSE;

660:         PetscHashIJKLIterGet(faceTable, iter, &f);
661:         DMPlexInsertCone(idm, c, cf, f);
662:         /* Orient face: Do not allow reverse orientation at the first vertex */
663:         DMPlexGetConeSize(idm, f, &coneSize);
664:         DMPlexGetCone(idm, f, &cone);
665:         if (coneSize != faceSizeH) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid number of face vertices %D for face %D should be %D", coneSize, f, faceSizeH);
666:         /* Hybrid faces are stored as tensor products of edges, so to compare them to normal faces, we have to flip */
667:         if (faceSize == 4 && c >= pMax && faceSizeAll != faceSizeAllT && f < pEnd[faceDepth] - faceH) {q2h = PETSC_TRUE; fv[2] = 3; fv[3] = 2;}
668:         /* - First find the initial vertex */
669:         for (i = 0; i < faceSizeH; ++i) if (cellFace[fv[0]] == cone[i]) break;
670:         if (q2h) { /* Matt's case: hybrid faces meeting with non-hybrid faces. This is a case that is not (and will not be) supported in general by the refinements */
671:           /* - Try forward comparison */
672:           for (j = 0; j < faceSizeH; ++j) if (cellFace[fv[j]] != cone[(i+j)%faceSizeH]) break;
673:           if (j == faceSizeH) {
674:             if ((faceSizeH == 2) && (i == 1)) ornt = -2;
675:             else                              ornt = i;
676:           } else {
677:             /* - Try backward comparison */
678:             for (j = 0; j < faceSizeH; ++j) if (cellFace[fv[j]] != cone[(i+faceSizeH-j)%faceSizeH]) break;
679:             if (j == faceSizeH) {
680:               if (i == 0) ornt = -faceSizeH;
681:               else        ornt = -i;
682:             } else SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Could not determine orientation of face %D in cell %D", f, c);
683:           }
684:         } else {
685:           /* when matching hybrid faces in 3D, only few cases are possible.
686:              Face traversal however can no longer follow the usual convention, this seems a serious issue to me */
687:           PetscInt tquad_map[4][4] = { {PETSC_MIN_INT,            0,PETSC_MIN_INT,PETSC_MIN_INT},
688:                                        {           -1,PETSC_MIN_INT,PETSC_MIN_INT,PETSC_MIN_INT},
689:                                        {PETSC_MIN_INT,PETSC_MIN_INT,PETSC_MIN_INT,            1},
690:                                        {PETSC_MIN_INT,PETSC_MIN_INT,           -2,PETSC_MIN_INT} };
691:           PetscInt i2;

693:           /* find the second vertex */
694:           for (i2 = 0; i2 < faceSizeH; ++i2) if (cellFace[fv[1]] == cone[i2]) break;
695:           switch (faceSizeH) {
696:           case 2:
697:             ornt = i ? -2 : 0;
698:             break;
699:           case 4:
700:             ornt = tquad_map[i][i2];
701:             break;
702:           default:
703:             SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Unhandled face size %D for face %D in cell %D", faceSizeH, f, c);

705:           }
706:         }
707:         DMPlexInsertConeOrientation(idm, c, cf, ornt);
708:       }
709:     }
710:     DMPlexRestoreRawFacesHybrid_Internal(dm, ct, cone, &numCellFaces, &numCellFacesN, &faceSize, &cellFaces);
711:   }
712:   if (face != pEnd[faceDepth]) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid number of faces %D should be %D", face-pStart[faceDepth], pEnd[faceDepth]-pStart[faceDepth]);
713:   PetscFree2(pStart,pEnd);
714:   PetscHashIJKLDestroy(&faceTable);
715:   PetscFree2(pStart,pEnd);
716:   DMPlexSetHybridBounds(idm, cMax, fMax, eMax, vMax);
717:   DMPlexSymmetrize(idm);
718:   DMPlexStratify(idm);
719:   return(0);
720: }

722: static PetscErrorCode SortRmineRremoteByRemote_Private(PetscSF sf, PetscInt *rmine1[], PetscInt *rremote1[])
723: {
724:   PetscInt            nleaves;
725:   PetscInt            nranks;
726:   const PetscMPIInt  *ranks=NULL;
727:   const PetscInt     *roffset=NULL, *rmine=NULL, *rremote=NULL;
728:   PetscInt            n, o, r;
729:   PetscErrorCode      ierr;

732:   PetscSFGetRootRanks(sf, &nranks, &ranks, &roffset, &rmine, &rremote);
733:   nleaves = roffset[nranks];
734:   PetscMalloc2(nleaves, rmine1, nleaves, rremote1);
735:   for (r=0; r<nranks; r++) {
736:     /* simultaneously sort rank-wise portions of rmine & rremote by values in rremote
737:        - to unify order with the other side */
738:     o = roffset[r];
739:     n = roffset[r+1] - o;
740:     PetscArraycpy(&(*rmine1)[o], &rmine[o], n);
741:     PetscArraycpy(&(*rremote1)[o], &rremote[o], n);
742:     PetscSortIntWithArray(n, &(*rremote1)[o], &(*rmine1)[o]);
743:   }
744:   return(0);
745: }

747: PetscErrorCode DMPlexOrientInterface_Internal(DM dm)
748: {
749:   /* Here we only compare first 2 points of the cone. Full cone size would lead to stronger self-checking. */
750:   PetscInt          masterCone[2];
751:   PetscInt          (*roots)[2], (*leaves)[2];
752:   PetscMPIInt       (*rootsRanks)[2], (*leavesRanks)[2];

754:   PetscSF           sf=NULL;
755:   const PetscInt    *locals=NULL;
756:   const PetscSFNode *remotes=NULL;
757:   PetscInt           nroots, nleaves, p, c;
758:   PetscInt           nranks, n, o, r;
759:   const PetscMPIInt *ranks=NULL;
760:   const PetscInt    *roffset=NULL;
761:   PetscInt          *rmine1=NULL, *rremote1=NULL; /* rmine and rremote copies simultaneously sorted by rank and rremote */
762:   const PetscInt    *cone=NULL;
763:   PetscInt           coneSize, ind0;
764:   MPI_Comm           comm;
765:   PetscMPIInt        rank,size;
766:   PetscInt           debug = 0;
767:   PetscErrorCode     ierr;

770:   DMGetPointSF(dm, &sf);
771:   PetscSFGetGraph(sf, &nroots, &nleaves, &locals, &remotes);
772:   if (nroots < 0) return(0);
773:   PetscSFSetUp(sf);
774:   PetscSFGetRootRanks(sf, &nranks, &ranks, &roffset, NULL, NULL);
775:   DMViewFromOptions(dm, NULL, "-before_fix_dm_view");
776: #if defined(PETSC_USE_DEBUG)
777:   DMPlexCheckPointSF(dm);
778: #endif
779:   SortRmineRremoteByRemote_Private(sf, &rmine1, &rremote1);
780:   PetscMalloc4(nroots, &roots, nroots, &leaves, nroots, &rootsRanks, nroots, &leavesRanks);
781:   PetscObjectGetComm((PetscObject) dm, &comm);
782:   MPI_Comm_rank(comm, &rank);
783:   MPI_Comm_size(comm, &size);
784:   for (p = 0; p < nroots; ++p) {
785:     DMPlexGetConeSize(dm, p, &coneSize);
786:     DMPlexGetCone(dm, p, &cone);
787:     if (coneSize < 2) {
788:       for (c = 0; c < 2; c++) {
789:         roots[p][c] = -1;
790:         rootsRanks[p][c] = -1;
791:       }
792:       continue;
793:     }
794:     /* Translate all points to root numbering */
795:     for (c = 0; c < 2; c++) {
796:       PetscFindInt(cone[c], nleaves, locals, &ind0);
797:       if (ind0 < 0) {
798:         roots[p][c] = cone[c];
799:         rootsRanks[p][c] = rank;
800:       } else {
801:         roots[p][c] = remotes[ind0].index;
802:         rootsRanks[p][c] = remotes[ind0].rank;
803:       }
804:     }
805:   }
806:   for (p = 0; p < nroots; ++p) {
807:     for (c = 0; c < 2; c++) {
808:       leaves[p][c] = -2;
809:       leavesRanks[p][c] = -2;
810:     }
811:   }
812:   PetscSFBcastBegin(sf, MPIU_2INT, roots, leaves);
813:   PetscSFBcastBegin(sf, MPI_2INT, rootsRanks, leavesRanks);
814:   PetscSFBcastEnd(sf, MPIU_2INT, roots, leaves);
815:   PetscSFBcastEnd(sf, MPI_2INT, rootsRanks, leavesRanks);
816:   if (debug) {PetscSynchronizedFlush(comm, NULL);}
817:   if (debug && rank == 0) {PetscSynchronizedPrintf(comm, "Referenced roots\n");}
818:   for (p = 0; p < nroots; ++p) {
819:     if (leaves[p][0] < 0) continue;
820:     DMPlexGetConeSize(dm, p, &coneSize);
821:     DMPlexGetCone(dm, p, &cone);
822:     if (debug) {PetscSynchronizedPrintf(comm, "[%d]  %4D: cone=[%4D %4D] roots=[(%d,%4D) (%d,%4D)] leaves=[(%d,%4D) (%d,%4D)]", rank, p, cone[0], cone[1], rootsRanks[p][0], roots[p][0], rootsRanks[p][1], roots[p][1], leavesRanks[p][0], leaves[p][0], leavesRanks[p][1], leaves[p][1]);}
823:     if ((leaves[p][0] != roots[p][0]) || (leaves[p][1] != roots[p][1]) || (leavesRanks[p][0] != rootsRanks[p][0]) || (leavesRanks[p][1] != rootsRanks[p][1])) {
824:       /* Translate these two leaves to my cone points; masterCone means desired order p's cone points */
825:       for (c = 0; c < 2; c++) {
826:         if (leavesRanks[p][c] == rank) {
827:           /* A local leave is just taken as it is */
828:           masterCone[c] = leaves[p][c];
829:           continue;
830:         }
831:         /* Find index of rank leavesRanks[p][c] among remote ranks */
832:         /* No need for PetscMPIIntCast because these integers were originally cast from PetscMPIInt. */
833:         PetscFindMPIInt((PetscMPIInt)leavesRanks[p][c], nranks, ranks, &r);
834:         if (PetscUnlikely(r < 0)) SETERRQ7(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Point %D cone[%D]=%D root (%d,%D) leave (%d,%D): leave rank not found among remote ranks",p,c,cone[c],rootsRanks[p][c],roots[p][c],leavesRanks[p][c],leaves[p][c]);
835:         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]);
836:         /* Find point leaves[p][c] among remote points aimed at rank leavesRanks[p][c] */
837:         o = roffset[r];
838:         n = roffset[r+1] - o;
839:         PetscFindInt(leaves[p][c], n, &rremote1[o], &ind0);
840:         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]);
841:         /* Get the corresponding local point */
842:         masterCone[c] = rmine1[o+ind0];
843:       }
844:       if (debug) {PetscSynchronizedPrintf(comm, " masterCone=[%4D %4D]\n", masterCone[0], masterCone[1]);}
845:       /* Set the desired order of p's cone points and fix orientations accordingly */
846:       /* Vaclav's note: Here we only compare first 2 points of the cone. Full cone size would lead to stronger self-checking. */
847:       DMPlexOrientCell(dm, p, 2, masterCone);
848:     } else if (debug) {PetscSynchronizedPrintf(comm, " ==\n");}
849:   }
850:   if (debug) {
851:     PetscSynchronizedFlush(comm, NULL);
852:     MPI_Barrier(comm);
853:   }
854:   DMViewFromOptions(dm, NULL, "-after_fix_dm_view");
855:   PetscFree4(roots, leaves, rootsRanks, leavesRanks);
856:   PetscFree2(rmine1, rremote1);
857:   return(0);
858: }

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

868:   PetscOptionsHasName(NULL, NULL, opt, &flg);
869:   if (!flg) return(0);
870:   MPI_Comm_rank(comm, &rank);
871:   PetscSynchronizedPrintf(comm, "[%d]%s:\n", rank, name);
872:   for (idx = 0; idx < n; ++idx) {PetscSynchronizedPrintf(comm, "[%d]%s %D %s %D\n", rank, idxname, idx, valname, a[idx]);}
873:   PetscSynchronizedFlush(comm, NULL);
874:   return(0);
875: }

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

885:   PetscOptionsHasName(NULL, NULL, opt, &flg);
886:   if (!flg) return(0);
887:   MPI_Comm_rank(comm, &rank);
888:   PetscSynchronizedPrintf(comm, "[%d]%s:\n", rank, name);
889:   if (idxname) {
890:     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);}
891:   } else {
892:     for (idx = 0; idx < n; ++idx) {PetscSynchronizedPrintf(comm, "[%d]rank %D index %D\n", rank, a[idx].rank, a[idx].index);}
893:   }
894:   PetscSynchronizedFlush(comm, NULL);
895:   return(0);
896: }

898: static PetscErrorCode DMPlexMapToLocalPoint(DM dm, PetscHMapIJ remotehash, PetscSFNode remotePoint, PetscInt *localPoint)
899: {
900:   PetscSF         sf;
901:   const PetscInt *locals;
902:   PetscMPIInt     rank;
903:   PetscErrorCode  ierr;

906:   MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank);
907:   DMGetPointSF(dm, &sf);
908:   PetscSFGetGraph(sf, NULL, NULL, &locals, NULL);
909:   if (remotePoint.rank == rank) {
910:     *localPoint = remotePoint.index;
911:   } else {
912:     PetscHashIJKey key;
913:     PetscInt       l;

915:     key.i = remotePoint.index;
916:     key.j = remotePoint.rank;
917:     PetscHMapIJGet(remotehash, key, &l);
918:     if (l >= 0) {
919:       *localPoint = locals[l];
920:     } else PetscFunctionReturn(1);
921:   }
922:   return(0);
923: }

925: static PetscErrorCode DMPlexMapToGlobalPoint(DM dm, PetscInt localPoint, PetscSFNode *remotePoint)
926: {
927:   PetscSF            sf;
928:   const PetscInt    *locals, *rootdegree;
929:   const PetscSFNode *remotes;
930:   PetscInt           Nl, l;
931:   PetscMPIInt        rank;
932:   PetscErrorCode     ierr;

935:   MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank);
936:   DMGetPointSF(dm, &sf);
937:   PetscSFGetGraph(sf, NULL, &Nl, &locals, &remotes);
938:   if (Nl < 0) goto owned;
939:   PetscSFComputeDegreeBegin(sf, &rootdegree);
940:   PetscSFComputeDegreeEnd(sf, &rootdegree);
941:   if (rootdegree[localPoint]) goto owned;
942:   PetscFindInt(localPoint, Nl, locals, &l);
943:   if (l < 0) PetscFunctionReturn(1);
944:   *remotePoint = remotes[l];
945:   return(0);
946:   owned:
947:   remotePoint->rank  = rank;
948:   remotePoint->index = localPoint;
949:   return(0);
950: }


953: static PetscErrorCode DMPlexPointIsShared(DM dm, PetscInt p, PetscBool *isShared)
954: {
955:   PetscSF         sf;
956:   const PetscInt *locals, *rootdegree;
957:   PetscInt        Nl, idx;
958:   PetscErrorCode  ierr;

961:   *isShared = PETSC_FALSE;
962:   DMGetPointSF(dm, &sf);
963:   PetscSFGetGraph(sf, NULL, &Nl, &locals, NULL);
964:   if (Nl < 0) return(0);
965:   PetscFindInt(p, Nl, locals, &idx);
966:   if (idx >= 0) {*isShared = PETSC_TRUE; return(0);}
967:   PetscSFComputeDegreeBegin(sf, &rootdegree);
968:   PetscSFComputeDegreeEnd(sf, &rootdegree);
969:   if (rootdegree[p] > 0) *isShared = PETSC_TRUE;
970:   return(0);
971: }

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

981:   DMPlexGetConeSize(dm, p, &coneSize);
982:   DMPlexGetCone(dm, p, &cone);
983:   for (c = 0; c < coneSize; ++c) {
984:     PetscBool pointShared;

986:     DMPlexPointIsShared(dm, cone[c], &pointShared);
987:     cShared = (PetscBool) (cShared && pointShared);
988:   }
989:   *isShared = coneSize ? cShared : PETSC_FALSE;
990:   return(0);
991: }

993: static PetscErrorCode DMPlexGetConeMinimum(DM dm, PetscInt p, PetscSFNode *cpmin)
994: {
995:   const PetscInt *cone;
996:   PetscInt        coneSize, c;
997:   PetscSFNode     cmin = {PETSC_MAX_INT, PETSC_MAX_INT}, missing = {-1, -1};
998:   PetscErrorCode  ierr;

1001:   DMPlexGetConeSize(dm, p, &coneSize);
1002:   DMPlexGetCone(dm, p, &cone);
1003:   for (c = 0; c < coneSize; ++c) {
1004:     PetscSFNode rcp;

1006:     DMPlexMapToGlobalPoint(dm, cone[c], &rcp);
1007:     if (ierr) {
1008:       cmin = missing;
1009:     } else {
1010:       cmin = (rcp.rank < cmin.rank) || (rcp.rank == cmin.rank && rcp.index < cmin.index) ? rcp : cmin;
1011:     }
1012:   }
1013:   *cpmin = coneSize ? cmin : missing;
1014:   return(0);
1015: }

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

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

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

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

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

1099:   Collective on dm

1101:   Input Parameters:
1102: + dm      - The interpolated DM
1103: - pointSF - The initial SF without interpolated points

1105:   Output Parameter:
1106: . pointSF - The SF including interpolated points

1108:   Level: developer

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

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

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

1178:     PetscHMapIJCreate(&faceHash);
1179:     for (l = 0; l < Nl; ++l) {
1180:       const PetscInt p = localPoints[l];

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

1444:   Collective on dm

1446:   Input Parameters:
1447: + dm - The DMPlex object with only cells and vertices
1448: - dmInt - The interpolated DM

1450:   Output Parameter:
1451: . dmInt - The complete DMPlex object

1453:   Level: intermediate

1455:   Notes:
1456:     It does not copy over the coordinates.

1458:   Developer Notes:
1459:     It sets plex->interpolated = DMPLEX_INTERPOLATED_FULL.

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

1476:   PetscLogEventBegin(DMPLEX_Interpolate,dm,0,0,0);
1477:   DMPlexGetDepth(dm, &depth);
1478:   DMGetDimension(dm, &dim);
1479:   DMPlexIsInterpolated(dm, &interpolated);
1480:   if (interpolated == DMPLEX_INTERPOLATED_PARTIAL) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Not for partially interpolated meshes");
1481:   if (interpolated == DMPLEX_INTERPOLATED_FULL) {
1482:     PetscObjectReference((PetscObject) dm);
1483:     idm  = dm;
1484:   } else {
1485:     for (d = 1; d < dim; ++d) {
1486:       /* Create interpolated mesh */
1487:       DMCreate(PetscObjectComm((PetscObject)dm), &idm);
1488:       DMSetType(idm, DMPLEX);
1489:       DMSetDimension(idm, dim);
1490:       if (depth > 0) {
1491:         DMPlexInterpolateFaces_Internal(odm, 1, idm);
1492:         DMGetPointSF(odm, &sfPoint);
1493:         {
1494:           /* TODO: We need to systematically fix cases of distributed Plexes with no graph set */
1495:           PetscInt nroots;
1496:           PetscSFGetGraph(sfPoint, &nroots, NULL, NULL, NULL);
1497:           if (nroots >= 0) {DMPlexInterpolatePointSF(idm, sfPoint);}
1498:         }
1499:       }
1500:       if (odm != dm) {DMDestroy(&odm);}
1501:       odm = idm;
1502:     }
1503:     PetscObjectGetName((PetscObject) dm,  &name);
1504:     PetscObjectSetName((PetscObject) idm,  name);
1505:     DMPlexCopyCoordinates(dm, idm);
1506:     DMCopyLabels(dm, idm, PETSC_COPY_VALUES, PETSC_FALSE);
1507:     PetscOptionsGetBool(((PetscObject)dm)->options, ((PetscObject)dm)->prefix, "-dm_plex_interpolate_orient_interfaces", &flg, NULL);
1508:     if (flg) {DMPlexOrientInterface_Internal(idm);}
1509:   }
1510:   {
1511:     PetscBool            isper;
1512:     const PetscReal      *maxCell, *L;
1513:     const DMBoundaryType *bd;

1515:     DMGetPeriodicity(dm,&isper,&maxCell,&L,&bd);
1516:     DMSetPeriodicity(idm,isper,maxCell,L,bd);
1517:   }
1518:   /* This function makes the mesh fully interpolated on all ranks */
1519:   {
1520:     DM_Plex *plex = (DM_Plex *) idm->data;
1521:     plex->interpolated = plex->interpolatedCollective = DMPLEX_INTERPOLATED_FULL;
1522:   }
1523:   *dmInt = idm;
1524:   PetscLogEventEnd(DMPLEX_Interpolate,dm,0,0,0);
1525:   return(0);
1526: }

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

1531:   Collective on dmA

1533:   Input Parameter:
1534: . dmA - The DMPlex object with initial coordinates

1536:   Output Parameter:
1537: . dmB - The DMPlex object with copied coordinates

1539:   Level: intermediate

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

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

1559:   if (dmA == dmB) return(0);
1560:   DMGetCoordinateDim(dmA, &cdim);
1561:   DMSetCoordinateDim(dmB, cdim);
1562:   DMPlexGetDepthStratum(dmA, 0, &vStartA, &vEndA);
1563:   DMPlexGetDepthStratum(dmB, 0, &vStartB, &vEndB);
1564:   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);
1565:   DMPlexGetHeightStratum(dmA, 0, &cStartA, &cEndA);
1566:   DMPlexGetHeightStratum(dmB, 0, &cStartB, &cEndB);
1567:   DMGetCoordinateSection(dmA, &coordSectionA);
1568:   DMGetCoordinateSection(dmB, &coordSectionB);
1569:   if (coordSectionA == coordSectionB) return(0);
1570:   PetscSectionGetNumFields(coordSectionA, &Nf);
1571:   if (!Nf) return(0);
1572:   if (Nf > 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "The number of coordinate fields must be 1, not %D", Nf);
1573:   if (!coordSectionB) {
1574:     PetscInt dim;

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

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

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

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

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

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

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

1653:   Collective on dm

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

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

1661:   Level: intermediate

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

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

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

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

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

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

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

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

1739:     /* Get original SF information */
1740:     DMGetPointSF(dm, &sfPoint);
1741:     DMGetPointSF(udm, &sfPointUn);
1742:     DMPlexGetDepthStratum(dm, 0, NULL, &vEnd);
1743:     PetscSFGetGraph(sfPoint, &numRoots, &numLeaves, &localPoints, &remotePoints);
1744:     /* Allocate space for cells and vertices */
1745:     for (l = 0; l < numLeaves; ++l) if (localPoints[l] < vEnd) numLeavesUn++;
1746:     /* Fill in leaves */
1747:     if (vEnd >= 0) {
1748:       PetscMalloc1(numLeavesUn, &remotePointsUn);
1749:       PetscMalloc1(numLeavesUn, &localPointsUn);
1750:       for (l = 0; l < numLeaves; l++) {
1751:         if (localPoints[l] < vEnd) {
1752:           localPointsUn[n]        = localPoints[l];
1753:           remotePointsUn[n].rank  = remotePoints[l].rank;
1754:           remotePointsUn[n].index = remotePoints[l].index;
1755:           ++n;
1756:         }
1757:       }
1758:       if (n != numLeavesUn) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Inconsistent number of leaves %d != %d", n, numLeavesUn);
1759:       PetscSFSetGraph(sfPointUn, vEnd, numLeavesUn, localPointsUn, PETSC_OWN_POINTER, remotePointsUn, PETSC_OWN_POINTER);
1760:     }
1761:   }
1762:   {
1763:     PetscBool            isper;
1764:     const PetscReal      *maxCell, *L;
1765:     const DMBoundaryType *bd;

1767:     DMGetPeriodicity(dm,&isper,&maxCell,&L,&bd);
1768:     DMSetPeriodicity(udm,isper,maxCell,L,bd);
1769:   }
1770:   /* This function makes the mesh fully uninterpolated on all ranks */
1771:   {
1772:     DM_Plex *plex = (DM_Plex *) udm->data;
1773:     plex->interpolated = plex->interpolatedCollective = DMPLEX_INTERPOLATED_NONE;
1774:   }
1775:   *dmUnint = udm;
1776:   return(0);
1777: }

1779: static PetscErrorCode DMPlexIsInterpolated_Internal(DM dm, DMPlexInterpolatedFlag *interpolated)
1780: {
1781:   PetscInt       coneSize, depth, dim, h, p, pStart, pEnd;
1782:   MPI_Comm       comm;

1786:   PetscObjectGetComm((PetscObject)dm, &comm);
1787:   DMPlexGetDepth(dm, &depth);
1788:   DMGetDimension(dm, &dim);

1790:   if (depth == dim) {
1791:     *interpolated = DMPLEX_INTERPOLATED_FULL;
1792:     if (!dim) goto finish;

1794:     /* Check points at height = dim are vertices (have no cones) */
1795:     DMPlexGetHeightStratum(dm, dim, &pStart, &pEnd);
1796:     for (p=pStart; p<pEnd; p++) {
1797:       DMPlexGetConeSize(dm, p, &coneSize);
1798:       if (coneSize) {
1799:         *interpolated = DMPLEX_INTERPOLATED_PARTIAL;
1800:         goto finish;
1801:       }
1802:     }

1804:     /* Check points at height < dim have cones */
1805:     for (h=0; h<dim; h++) {
1806:       DMPlexGetHeightStratum(dm, h, &pStart, &pEnd);
1807:       for (p=pStart; p<pEnd; p++) {
1808:         DMPlexGetConeSize(dm, p, &coneSize);
1809:         if (!coneSize) {
1810:           *interpolated = DMPLEX_INTERPOLATED_PARTIAL;
1811:           goto finish;
1812:         }
1813:       }
1814:     }
1815:   } else if (depth == 1) {
1816:     *interpolated = DMPLEX_INTERPOLATED_NONE;
1817:   } else {
1818:     *interpolated = DMPLEX_INTERPOLATED_PARTIAL;
1819:   }
1820: finish:
1821:   return(0);
1822: }

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

1827:   Not Collective

1829:   Input Parameter:
1830: . dm      - The DM object

1832:   Output Parameter:
1833: . interpolated - Flag whether the DM is interpolated

1835:   Level: intermediate

1837:   Notes:
1838:   Unlike DMPlexIsInterpolatedCollective(), this is NOT collective
1839:   so the results can be different on different ranks in special cases.
1840:   However, DMPlexInterpolate() guarantees the result is the same on all.

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

1844:   Developer Notes:
1845:   Initially, plex->interpolated = DMPLEX_INTERPOLATED_INVALID.

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

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

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

1856: .seealso: DMPlexInterpolate(), DMPlexIsInterpolatedCollective()
1857: @*/
1858: PetscErrorCode DMPlexIsInterpolated(DM dm, DMPlexInterpolatedFlag *interpolated)
1859: {
1860:   DM_Plex        *plex = (DM_Plex *) dm->data;
1861:   PetscErrorCode  ierr;

1866:   if (plex->interpolated < 0) {
1867:     DMPlexIsInterpolated_Internal(dm, &plex->interpolated);
1868:   } else {
1869: #if defined (PETSC_USE_DEBUG)
1870:     DMPlexInterpolatedFlag flg;

1872:     DMPlexIsInterpolated_Internal(dm, &flg);
1873:     if (flg != plex->interpolated) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Stashed DMPlexInterpolatedFlag %s is inconsistent with current %s", DMPlexInterpolatedFlags[plex->interpolated], DMPlexInterpolatedFlags[flg]);
1874: #endif
1875:   }
1876:   *interpolated = plex->interpolated;
1877:   return(0);
1878: }

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

1883:   Collective

1885:   Input Parameter:
1886: . dm      - The DM object

1888:   Output Parameter:
1889: . interpolated - Flag whether the DM is interpolated

1891:   Level: intermediate

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

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

1898:   Developer Notes:
1899:   Initially, plex->interpolatedCollective = DMPLEX_INTERPOLATED_INVALID.

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

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

1908: .seealso: DMPlexInterpolate(), DMPlexIsInterpolated()
1909: @*/
1910: PetscErrorCode DMPlexIsInterpolatedCollective(DM dm, DMPlexInterpolatedFlag *interpolated)
1911: {
1912:   DM_Plex        *plex = (DM_Plex *) dm->data;
1913:   PetscBool       debug=PETSC_FALSE;
1914:   PetscErrorCode  ierr;

1919:   PetscOptionsGetBool(((PetscObject) dm)->options, ((PetscObject) dm)->prefix, "-dm_plex_is_interpolated_collective_debug", &debug, NULL);
1920:   if (plex->interpolatedCollective < 0) {
1921:     DMPlexInterpolatedFlag  min, max;
1922:     MPI_Comm                comm;

1924:     PetscObjectGetComm((PetscObject)dm, &comm);
1925:     DMPlexIsInterpolated(dm, &plex->interpolatedCollective);
1926:     MPI_Allreduce(&plex->interpolatedCollective, &min, 1, MPIU_ENUM, MPI_MIN, comm);
1927:     MPI_Allreduce(&plex->interpolatedCollective, &max, 1, MPIU_ENUM, MPI_MAX, comm);
1928:     if (min != max) plex->interpolatedCollective = DMPLEX_INTERPOLATED_MIXED;
1929:     if (debug) {
1930:       PetscMPIInt rank;

1932:       MPI_Comm_rank(comm, &rank);
1933:       PetscSynchronizedPrintf(comm, "[%d] interpolated=%s interpolatedCollective=%s\n", rank, DMPlexInterpolatedFlags[plex->interpolated], DMPlexInterpolatedFlags[plex->interpolatedCollective]);
1934:       PetscSynchronizedFlush(comm, PETSC_STDOUT);
1935:     }
1936:   }
1937:   *interpolated = plex->interpolatedCollective;
1938:   return(0);
1939: }