Actual source code: plexinterpolate.c

petsc-master 2019-10-20
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)


 23: /*
 24:   DMPlexGetFaces_Internal - Gets groups of vertices that correspond to faces for the given cell
 25:   This assumes that the mesh is not interpolated from the depth of point p to the vertices
 26: */
 27: PetscErrorCode DMPlexGetFaces_Internal(DM dm, PetscInt dim, PetscInt p, PetscInt *numFaces, PetscInt *faceSize, const PetscInt *faces[])
 28: {
 29:   const PetscInt *cone = NULL;
 30:   PetscInt        coneSize;
 31:   PetscErrorCode  ierr;

 35:   DMPlexGetConeSize(dm, p, &coneSize);
 36:   DMPlexGetCone(dm, p, &cone);
 37:   DMPlexGetRawFaces_Internal(dm, dim, coneSize, cone, numFaces, faceSize, faces);
 38:   return(0);
 39: }

 41: /*
 42:   DMPlexRestoreFaces_Internal - Restores the array
 43: */
 44: PetscErrorCode DMPlexRestoreFaces_Internal(DM dm, PetscInt dim, PetscInt p, PetscInt *numFaces, PetscInt *faceSize, const PetscInt *faces[])
 45: {
 46:   PetscErrorCode  ierr;

 49:   if (faces) { DMRestoreWorkArray(dm, 0, MPIU_INT, (void *) faces); }
 50:   return(0);
 51: }

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

 65:   DMPlexGetMaxSizes(dm, &maxConeSize, &maxSupportSize);
 66:   if (faces) {DMGetWorkArray(dm, PetscSqr(PetscMax(maxConeSize, maxSupportSize)), MPIU_INT, &facesTmp);}
 67:   switch (dim) {
 68:   case 1:
 69:     switch (coneSize) {
 70:     case 2:
 71:       if (faces) {
 72:         facesTmp[0] = cone[0]; facesTmp[1] = cone[1];
 73:         *faces = facesTmp;
 74:       }
 75:       if (numFaces) *numFaces = 2;
 76:       if (faceSize) *faceSize = 1;
 77:       break;
 78:     default:
 79:       SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Cone size %D not supported for dimension %D", coneSize, dim);
 80:     }
 81:     break;
 82:   case 2:
 83:     switch (coneSize) {
 84:     case 3:
 85:       if (faces) {
 86:         facesTmp[0] = cone[0]; facesTmp[1] = cone[1];
 87:         facesTmp[2] = cone[1]; facesTmp[3] = cone[2];
 88:         facesTmp[4] = cone[2]; facesTmp[5] = cone[0];
 89:         *faces = facesTmp;
 90:       }
 91:       if (numFaces) *numFaces = 3;
 92:       if (faceSize) *faceSize = 2;
 93:       break;
 94:     case 4:
 95:       /* Vertices follow right hand rule */
 96:       if (faces) {
 97:         facesTmp[0] = cone[0]; facesTmp[1] = cone[1];
 98:         facesTmp[2] = cone[1]; facesTmp[3] = cone[2];
 99:         facesTmp[4] = cone[2]; facesTmp[5] = cone[3];
100:         facesTmp[6] = cone[3]; facesTmp[7] = cone[0];
101:         *faces = facesTmp;
102:       }
103:       if (numFaces) *numFaces = 4;
104:       if (faceSize) *faceSize = 2;
105:       break;
106:     default:
107:       SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Cone size %D not supported for dimension %D", coneSize, dim);
108:     }
109:     break;
110:   case 3:
111:     switch (coneSize) {
112:     case 3:
113:       if (faces) {
114:         facesTmp[0] = cone[0]; facesTmp[1] = cone[1];
115:         facesTmp[2] = cone[1]; facesTmp[3] = cone[2];
116:         facesTmp[4] = cone[2]; facesTmp[5] = cone[0];
117:         *faces = facesTmp;
118:       }
119:       if (numFaces) *numFaces = 3;
120:       if (faceSize) *faceSize = 2;
121:       break;
122:     case 4:
123:       /* Vertices of first face follow right hand rule and normal points away from last vertex */
124:       if (faces) {
125:         facesTmp[0] = cone[0]; facesTmp[1]  = cone[1]; facesTmp[2]  = cone[2];
126:         facesTmp[3] = cone[0]; facesTmp[4]  = cone[3]; facesTmp[5]  = cone[1];
127:         facesTmp[6] = cone[0]; facesTmp[7]  = cone[2]; facesTmp[8]  = cone[3];
128:         facesTmp[9] = cone[2]; facesTmp[10] = cone[1]; facesTmp[11] = cone[3];
129:         *faces = facesTmp;
130:       }
131:       if (numFaces) *numFaces = 4;
132:       if (faceSize) *faceSize = 3;
133:       break;
134:     case 8:
135:       /*  7--------6
136:          /|       /|
137:         / |      / |
138:        4--------5  |
139:        |  |     |  |
140:        |  |     |  |
141:        |  1--------2
142:        | /      | /
143:        |/       |/
144:        0--------3
145:        */
146:       if (faces) {
147:         facesTmp[0]  = cone[0]; facesTmp[1]  = cone[1]; facesTmp[2]  = cone[2]; facesTmp[3]  = cone[3]; /* Bottom */
148:         facesTmp[4]  = cone[4]; facesTmp[5]  = cone[5]; facesTmp[6]  = cone[6]; facesTmp[7]  = cone[7]; /* Top */
149:         facesTmp[8]  = cone[0]; facesTmp[9]  = cone[3]; facesTmp[10] = cone[5]; facesTmp[11] = cone[4]; /* Front */
150:         facesTmp[12] = cone[2]; facesTmp[13] = cone[1]; facesTmp[14] = cone[7]; facesTmp[15] = cone[6]; /* Back */
151:         facesTmp[16] = cone[3]; facesTmp[17] = cone[2]; facesTmp[18] = cone[6]; facesTmp[19] = cone[5]; /* Right */
152:         facesTmp[20] = cone[0]; facesTmp[21] = cone[4]; facesTmp[22] = cone[7]; facesTmp[23] = cone[1]; /* Left */
153:         *faces = facesTmp;
154:       }
155:       if (numFaces) *numFaces = 6;
156:       if (faceSize) *faceSize = 4;
157:       break;
158:     default:
159:       SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Cone size %D not supported for dimension %D", coneSize, dim);
160:     }
161:     break;
162:   default:
163:     SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Dimension %D not supported", dim);
164:   }
165:   return(0);
166: }

168: /*
169:   DMPlexGetRawFacesHybrid_Internal - Gets groups of vertices that correspond to faces for the given cone using hybrid ordering (prisms)
170: */
171: static PetscErrorCode DMPlexGetRawFacesHybrid_Internal(DM dm, PetscInt dim, PetscInt coneSize, const PetscInt cone[], PetscInt *numFaces, PetscInt *numFacesNotH, PetscInt *faceSize, const PetscInt *faces[])
172: {
173:   PetscInt       *facesTmp;
174:   PetscInt        maxConeSize, maxSupportSize;
175:   PetscErrorCode  ierr;

180:   DMPlexGetMaxSizes(dm, &maxConeSize, &maxSupportSize);
181:   if (faces) {DMGetWorkArray(dm, PetscSqr(PetscMax(maxConeSize, maxSupportSize)), MPIU_INT, &facesTmp);}
182:   switch (dim) {
183:   case 1:
184:     switch (coneSize) {
185:     case 2:
186:       if (faces) {
187:         facesTmp[0] = cone[0]; facesTmp[1] = cone[1];
188:         *faces = facesTmp;
189:       }
190:       if (numFaces)     *numFaces = 2;
191:       if (numFacesNotH) *numFacesNotH = 2;
192:       if (faceSize)     *faceSize = 1;
193:       break;
194:     default:
195:       SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Cone size %D not supported for dimension %D", coneSize, dim);
196:     }
197:     break;
198:   case 2:
199:     switch (coneSize) {
200:     case 4:
201:       if (faces) {
202:         facesTmp[0] = cone[0]; facesTmp[1] = cone[1];
203:         facesTmp[2] = cone[2]; facesTmp[3] = cone[3];
204:         facesTmp[4] = cone[0]; facesTmp[5] = cone[2];
205:         facesTmp[6] = cone[1]; facesTmp[7] = cone[3];
206:         *faces = facesTmp;
207:       }
208:       if (numFaces)     *numFaces = 4;
209:       if (numFacesNotH) *numFacesNotH = 2;
210:       if (faceSize)     *faceSize = 2;
211:       break;
212:     default:
213:       SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Cone size %D not supported for dimension %D", coneSize, dim);
214:     }
215:     break;
216:   case 3:
217:     switch (coneSize) {
218:     case 6: /* triangular prism */
219:       if (faces) {
220:         facesTmp[0]  = cone[0]; facesTmp[1]  = cone[1]; facesTmp[2]  = cone[2]; facesTmp[3]  = -1;      /* Bottom */
221:         facesTmp[4]  = cone[3]; facesTmp[5]  = cone[4]; facesTmp[6]  = cone[5]; facesTmp[7]  = -1;      /* Top */
222:         facesTmp[8]  = cone[0]; facesTmp[9]  = cone[1]; facesTmp[10] = cone[3]; facesTmp[11] = cone[4]; /* Back left */
223:         facesTmp[12] = cone[1]; facesTmp[13] = cone[2]; facesTmp[14] = cone[4]; facesTmp[15] = cone[5]; /* Back right */
224:         facesTmp[16] = cone[2]; facesTmp[17] = cone[0]; facesTmp[18] = cone[5]; facesTmp[19] = cone[3]; /* Front */
225:         *faces = facesTmp;
226:       }
227:       if (numFaces)     *numFaces = 5;
228:       if (numFacesNotH) *numFacesNotH = 2;
229:       if (faceSize)     *faceSize = -4;
230:       break;
231:     default:
232:       SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Cone size %D not supported for dimension %D", coneSize, dim);
233:     }
234:     break;
235:   default:
236:     SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Dimension %D not supported", dim);
237:   }
238:   return(0);
239: }

241: static PetscErrorCode DMPlexRestoreRawFacesHybrid_Internal(DM dm, PetscInt dim, PetscInt coneSize, 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 dim, PetscInt p, PetscInt *numFaces, PetscInt *numFacesNotH, PetscInt *faceSize, const PetscInt *faces[])
251: {
252:   const PetscInt *cone = NULL;
253:   PetscInt        coneSize;
254:   PetscErrorCode  ierr;

258:   DMPlexGetConeSize(dm, p, &coneSize);
259:   DMPlexGetCone(dm, p, &cone);
260:   DMPlexGetRawFacesHybrid_Internal(dm, dim, coneSize, 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:     PetscInt        numCellFacesT, faceSize, cf;

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

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

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

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

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

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

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

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

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

392:         if (faceSizeInc == 2) {
393:           key.i = PetscMin(cellFace[0], cellFace[1]);
394:           key.j = PetscMax(cellFace[0], cellFace[1]);
395:           key.k = PETSC_MAX_INT;
396:           key.l = PETSC_MAX_INT;
397:         } else {
398:           key.i = cellFace[0];
399:           key.j = cellFace[1];
400:           key.k = cellFace[2];
401:           key.l = faceSize > 3 ? (cellFace[3] < 0 ? faceSizeH = 3, PETSC_MAX_INT : cellFace[3]) : PETSC_MAX_INT;
402:           PetscSortInt(faceSize, (PetscInt *) &key);
403:         }
404:         /* this check is redundant for non-hybrid meshes */
405:         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);
406:         PetscHashIJKLPut(faceTable, key, &iter, &missing);
407:         if (missing) {
408:           PetscHashIJKLIterSet(faceTable, iter, face++);
409:           if (c >= pMax) ++faceT;
410:         }
411:       }
412:       if (c < pMax) {
413:         DMPlexRestoreFaces_Internal(dm, cellDim, c, &numCellFaces, &faceSize, &cellFaces);
414:       } else {
415:         DMPlexRestoreRawFacesHybrid_Internal(dm, cellDim, coneSizeH, NULL, NULL, NULL, NULL, &cellFaces);
416:       }
417:     }
418:   }
419:   pEnd[faceDepth] = face;

421:   /* Second pass for hybrid meshes: number hybrid faces */
422:   for (c = pMax; c < pEnd[cellDepth]; ++c) {
423:     const PetscInt *cellFaces, *cone;
424:     PetscInt        numCellFaces, numCellFacesN, faceSize, cf, coneSize;

426:     DMPlexGetConeSize(dm, c, &coneSize);
427:     DMPlexGetCone(dm, c, &cone);
428:     DMPlexGetRawFacesHybrid_Internal(dm, cellDim, coneSize, cone, &numCellFaces, &numCellFacesN, &faceSize, &cellFaces);
429:     if (numCellFaces != numCellFacesH) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unexpected hybrid numCellFaces %D != %D", numCellFaces, numCellFacesH);
430:     faceSize = PetscMax(faceSize, -faceSize);
431:     for (cf = numCellFacesN; cf < numCellFaces; ++cf) { /* These are the hybrid faces */
432:       const PetscInt   *cellFace = &cellFaces[cf*faceSize];
433:       PetscHashIJKLKey  key;
434:       PetscHashIter     iter;
435:       PetscBool         missing;
436:       PetscInt          faceSizeH = faceSize;

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

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

519:     for (p = pStart[d]; p < pEnd[d]; ++p) {
520:       DMPlexGetCone(dm, p, &cone);
521:       DMPlexSetCone(idm, p, cone);
522:       DMPlexGetConeOrientation(dm, p, &cone);
523:       DMPlexSetConeOrientation(idm, p, cone);
524:     }
525:   }
526:   for (outerloop = 0, face = pStart[faceDepth]; outerloop < 2; outerloop++) {
527:     PetscInt start, end;

529:     start = outerloop == 0 ? pMax : pStart[cellDepth];
530:     end = outerloop == 0 ? pEnd[cellDepth] : pMax;
531:     for (c = start; c < end; ++c) {
532:       const PetscInt *cellFaces;
533:       PetscInt        numCellFaces, faceSize, faceSizeInc, cf;

535:       if (c < pMax) {
536:         DMPlexGetFaces_Internal(dm, cellDim, c, &numCellFaces, &faceSize, &cellFaces);
537:         if (faceSize != faceSizeAll) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Inconsistent face for cell %D of size %D != %D", c, faceSize, faceSizeAll);
538:       } else {
539:         const PetscInt *cone;
540:         PetscInt        numCellFacesN, coneSize;

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

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

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

615:     DMPlexGetConeSize(dm, c, &coneSize);
616:     DMPlexGetCone(dm, c, &cone);
617:     DMPlexGetRawFacesHybrid_Internal(dm, cellDim, coneSize, cone, &numCellFaces, &numCellFacesN, &faceSize, &cellFaces);
618:     if (numCellFaces != numCellFacesH) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unexpected hybrid numCellFaces %D != %D", numCellFaces, numCellFacesH);
619:     faceSize = PetscMax(faceSize, -faceSize);
620:     for (cf = numCellFacesN; cf < numCellFaces; ++cf) { /* These are the hybrid faces */
621:       const PetscInt   *cellFace = &cellFaces[cf*faceSize];
622:       PetscHashIJKLKey key;
623:       PetscHashIter    iter;
624:       PetscBool        missing;
625:       PetscInt         faceSizeH = faceSize;

627:       if (faceSize == 2) {
628:         key.i = PetscMin(cellFace[0], cellFace[1]);
629:         key.j = PetscMax(cellFace[0], cellFace[1]);
630:         key.k = PETSC_MAX_INT;
631:         key.l = PETSC_MAX_INT;
632:       } else {
633:         key.i = cellFace[0];
634:         key.j = cellFace[1];
635:         key.k = cellFace[2];
636:         key.l = faceSize > 3 ? (cellFace[3] < 0 ? faceSizeH = 3, PETSC_MAX_INT : cellFace[3]) : PETSC_MAX_INT;
637:         PetscSortInt(faceSize, (PetscInt *) &key);
638:       }
639:       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);
640:       PetscHashIJKLPut(faceTable, key, &iter, &missing);
641:       if (missing) {
642:         DMPlexSetCone(idm, face, cellFace);
643:         PetscHashIJKLIterSet(faceTable, iter, face);
644:         DMPlexInsertCone(idm, c, cf, face++);
645:       } else {
646:         PetscInt        fv[4] = {0, 1, 2, 3};
647:         const PetscInt *cone;
648:         PetscInt        coneSize, ornt, i, j, f;
649:         PetscBool       q2h = PETSC_FALSE;

651:         PetscHashIJKLIterGet(faceTable, iter, &f);
652:         DMPlexInsertCone(idm, c, cf, f);
653:         /* Orient face: Do not allow reverse orientation at the first vertex */
654:         DMPlexGetConeSize(idm, f, &coneSize);
655:         DMPlexGetCone(idm, f, &cone);
656:         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);
657:         /* Hybrid faces are stored as tensor products of edges, so to compare them to normal faces, we have to flip */
658:         if (faceSize == 4 && c >= pMax && faceSizeAll != faceSizeAllT && f < pEnd[faceDepth] - faceH) {q2h = PETSC_TRUE; fv[2] = 3; fv[3] = 2;}
659:         /* - First find the initial vertex */
660:         for (i = 0; i < faceSizeH; ++i) if (cellFace[fv[0]] == cone[i]) break;
661:         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 */
662:           /* - Try forward comparison */
663:           for (j = 0; j < faceSizeH; ++j) if (cellFace[fv[j]] != cone[(i+j)%faceSizeH]) break;
664:           if (j == faceSizeH) {
665:             if ((faceSizeH == 2) && (i == 1)) ornt = -2;
666:             else                              ornt = i;
667:           } else {
668:             /* - Try backward comparison */
669:             for (j = 0; j < faceSizeH; ++j) if (cellFace[fv[j]] != cone[(i+faceSizeH-j)%faceSizeH]) break;
670:             if (j == faceSizeH) {
671:               if (i == 0) ornt = -faceSizeH;
672:               else        ornt = -i;
673:             } else SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Could not determine orientation of face %D in cell %D", f, c);
674:           }
675:         } else {
676:           /* when matching hybrid faces in 3D, only few cases are possible.
677:              Face traversal however can no longer follow the usual convention, this seems a serious issue to me */
678:           PetscInt tquad_map[4][4] = { {PETSC_MIN_INT,            0,PETSC_MIN_INT,PETSC_MIN_INT},
679:                                        {           -1,PETSC_MIN_INT,PETSC_MIN_INT,PETSC_MIN_INT},
680:                                        {PETSC_MIN_INT,PETSC_MIN_INT,PETSC_MIN_INT,            1},
681:                                        {PETSC_MIN_INT,PETSC_MIN_INT,           -2,PETSC_MIN_INT} };
682:           PetscInt i2;

684:           /* find the second vertex */
685:           for (i2 = 0; i2 < faceSizeH; ++i2) if (cellFace[fv[1]] == cone[i2]) break;
686:           switch (faceSizeH) {
687:           case 2:
688:             ornt = i ? -2 : 0;
689:             break;
690:           case 4:
691:             ornt = tquad_map[i][i2];
692:             break;
693:           default:
694:             SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Unhandled face size %D for face %D in cell %D", faceSizeH, f, c);

696:           }
697:         }
698:         DMPlexInsertConeOrientation(idm, c, cf, ornt);
699:       }
700:     }
701:     DMPlexRestoreRawFacesHybrid_Internal(dm, cellDim, coneSize, cone, &numCellFaces, &numCellFacesN, &faceSize, &cellFaces);
702:   }
703:   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]);
704:   PetscFree2(pStart,pEnd);
705:   PetscHashIJKLDestroy(&faceTable);
706:   PetscFree2(pStart,pEnd);
707:   DMPlexSetHybridBounds(idm, cMax, fMax, eMax, vMax);
708:   DMPlexSymmetrize(idm);
709:   DMPlexStratify(idm);
710:   return(0);
711: }

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

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

738: PetscErrorCode DMPlexOrientInterface_Internal(DM dm)
739: {
740:   /* Here we only compare first 2 points of the cone. Full cone size would lead to stronger self-checking. */
741:   PetscInt          masterCone[2];
742:   PetscInt          (*roots)[2], (*leaves)[2];
743:   PetscMPIInt       (*rootsRanks)[2], (*leavesRanks)[2];

745:   PetscSF           sf=NULL;
746:   const PetscInt    *locals=NULL;
747:   const PetscSFNode *remotes=NULL;
748:   PetscInt           nroots, nleaves, p, c;
749:   PetscInt           nranks, n, o, r;
750:   const PetscMPIInt *ranks=NULL;
751:   const PetscInt    *roffset=NULL;
752:   PetscInt          *rmine1=NULL, *rremote1=NULL; /* rmine and rremote copies simultaneously sorted by rank and rremote */
753:   const PetscInt    *cone=NULL;
754:   PetscInt           coneSize, ind0;
755:   MPI_Comm           comm;
756:   PetscMPIInt        rank,size;
757:   PetscInt           debug = 0;
758:   PetscErrorCode     ierr;

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

851: static PetscErrorCode IntArrayViewFromOptions(MPI_Comm comm, const char opt[], const char name[], const char idxname[], const char valname[], PetscInt n, const PetscInt a[])
852: {
853:   PetscInt       idx;
854:   PetscMPIInt    rank;
855:   PetscBool      flg;

859:   PetscOptionsHasName(NULL, NULL, opt, &flg);
860:   if (!flg) return(0);
861:   MPI_Comm_rank(comm, &rank);
862:   PetscSynchronizedPrintf(comm, "[%d]%s:\n", rank, name);
863:   for (idx = 0; idx < n; ++idx) {PetscSynchronizedPrintf(comm, "[%d]%s %D %s %D\n", rank, idxname, idx, valname, a[idx]);}
864:   PetscSynchronizedFlush(comm, NULL);
865:   return(0);
866: }

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

876:   PetscOptionsHasName(NULL, NULL, opt, &flg);
877:   if (!flg) return(0);
878:   MPI_Comm_rank(comm, &rank);
879:   PetscSynchronizedPrintf(comm, "[%d]%s:\n", rank, name);
880:   if (idxname) {
881:     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);}
882:   } else {
883:     for (idx = 0; idx < n; ++idx) {PetscSynchronizedPrintf(comm, "[%d]rank %D index %D\n", rank, a[idx].rank, a[idx].index);}
884:   }
885:   PetscSynchronizedFlush(comm, NULL);
886:   return(0);
887: }

889: static PetscErrorCode DMPlexMapToLocalPoint(PetscHMapIJ roothash, const PetscInt localPoints[], PetscMPIInt rank, PetscSFNode remotePoint, PetscInt *localPoint)
890: {

894:   if (remotePoint.rank == rank) {
895:     *localPoint = remotePoint.index;
896:   } else {
897:     PetscHashIJKey key;
898:     PetscInt       root;

900:     key.i = remotePoint.index;
901:     key.j = remotePoint.rank;
902:     PetscHMapIJGet(roothash, key, &root);
903:     if (root >= 0) {
904:       *localPoint = localPoints[root];
905:     } else PetscFunctionReturn(1);
906:   }
907:   return(0);
908: }

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

913:   Collective on dm

915:   Input Parameters:
916: + dm      - The interpolated DM
917: - pointSF - The initial SF without interpolated points

919:   Output Parameter:
920: . pointSF - The SF including interpolated points

922:   Level: intermediate

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

926: .seealso: DMPlexInterpolate(), DMPlexUninterpolate()
927: @*/
928: PetscErrorCode DMPlexInterpolatePointSF(DM dm, PetscSF pointSF)
929: {
930:   /*
931:        Okay, the algorithm is:
932:          - Take each point in the overlap (root)
933:          - Look at the neighboring points in the overlap (candidates)
934:          - Send these candidate points to neighbors
935:          - Neighbor checks for edge between root and candidate
936:          - If edge is found, it replaces candidate point with edge point
937:          - Send back the overwritten candidates (claims)
938:          - Original guy checks for edges, different from original candidate, and gets its own edge
939:          - This pair is put into SF

941:        We need a new algorithm that tolerates groups larger than 2.
942:          - Take each point in the overlap (root)
943:          - Find all collections of points in the overlap which make faces (do early join)
944:          - Send collections as candidates (add size as first number)
945:            - Make sure to send collection to all owners of all overlap points in collection
946:          - Neighbor check for face in collections
947:          - If face is found, it replaces candidate point with face point
948:          - Send back the overwritten candidates (claims)
949:          - Original guy checks for faces, different from original candidate, and gets its own face
950:          - This pair is put into SF
951:   */
952:   PetscHMapI         leafhash;
953:   PetscHMapIJ        roothash;
954:   const PetscInt    *localPoints, *rootdegree;
955:   const PetscSFNode *remotePoints;
956:   PetscSFNode       *candidates, *candidatesRemote, *claims;
957:   PetscSection       candidateSection, candidateSectionRemote, claimSection;
958:   PetscInt           numLeaves, l, numRoots, r, candidatesSize, candidatesRemoteSize;
959:   PetscMPIInt        size, rank;
960:   PetscHashIJKey     key;
961:   PetscBool          debug = PETSC_FALSE;
962:   PetscErrorCode     ierr;

965:   PetscOptionsHasName(NULL, ((PetscObject) dm)->prefix, "-dmplex_interp_debug", &debug);
966:   MPI_Comm_size(PetscObjectComm((PetscObject) dm), &size);
967:   MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank);
968:   PetscSFGetGraph(pointSF, &numRoots, &numLeaves, &localPoints, &remotePoints);
969:   if (size < 2 || numRoots < 0) return(0);
970:   DMPlexGetOverlap(dm, &r);
971:   if (r) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Interpolation of overlapped DMPlex not implemented yet");
972:   PetscObjectViewFromOptions((PetscObject) dm, NULL, "-dm_interp_pre_view");
973:   PetscObjectViewFromOptions((PetscObject) pointSF, NULL, "-petscsf_interp_pre_view");
974:   PetscLogEventBegin(DMPLEX_InterpolateSF,dm,0,0,0);
975:   /* Build hashes of points in the SF for efficient lookup */
976:   PetscHMapICreate(&leafhash);
977:   PetscHMapIJCreate(&roothash);
978:   for (l = 0; l < numLeaves; ++l) {
979:     key.i = remotePoints[l].index;
980:     key.j = remotePoints[l].rank;
981:     PetscHMapISet(leafhash, localPoints[l], l);
982:     PetscHMapIJSet(roothash, key, l);
983:   }
984:   /* Compute root degree to identify shared points */
985:   PetscSFComputeDegreeBegin(pointSF, &rootdegree);
986:   PetscSFComputeDegreeEnd(pointSF, &rootdegree);
987:   IntArrayViewFromOptions(PetscObjectComm((PetscObject) dm), "-interp_root_degree_view", "Root degree", "point", "degree", numRoots, rootdegree);
988:   /* Build a section / SFNode array of candidate points (face bd points) in the cone(support(leaf)),
989:      where each candidate is defined by a set of remote points (roots) for the other points that define the face. */
990:   PetscSectionCreate(PetscObjectComm((PetscObject) dm), &candidateSection);
991:   PetscSectionSetChart(candidateSection, 0, numRoots);
992:   {
993:     PetscHMapIJ facehash;

995:     PetscHMapIJCreate(&facehash);
996:     for (l = 0; l < numLeaves; ++l) {
997:       const PetscInt    localPoint = localPoints[l];
998:       const PetscInt   *support;
999:       PetscInt          supportSize, s;

1001:       if (debug) {PetscSynchronizedPrintf(PetscObjectComm((PetscObject) dm), "[%d]  Checking local point %D\n", rank, localPoint);}
1002:       DMPlexGetSupportSize(dm, localPoint, &supportSize);
1003:       DMPlexGetSupport(dm, localPoint, &support);
1004:       for (s = 0; s < supportSize; ++s) {
1005:         const PetscInt  face = support[s];
1006:         const PetscInt *cone;
1007:         PetscInt        coneSize, c, f, root;
1008:         PetscBool       isFace = PETSC_TRUE;

1010:         /* Only add face once */
1011:         if (debug) {PetscSynchronizedPrintf(PetscObjectComm((PetscObject) dm), "[%d]    Support point %D\n", rank, face);}
1012:         key.i = localPoint;
1013:         key.j = face;
1014:         PetscHMapIJGet(facehash, key, &f);
1015:         if (f >= 0) continue;
1016:         DMPlexGetConeSize(dm, face, &coneSize);
1017:         DMPlexGetCone(dm, face, &cone);
1018:         /* If a cone point does not map to leaves on any proc, then do not put face in SF */
1019:         for (c = 0; c < coneSize; ++c) {
1020:           if (debug) {PetscSynchronizedPrintf(PetscObjectComm((PetscObject) dm), "[%d]      Cone point %D\n", rank, cone[c]);}
1021:           PetscHMapIGet(leafhash, cone[c], &root);
1022:           if (!rootdegree[cone[c]] && (root < 0)) {isFace = PETSC_FALSE; break;}
1023:         }
1024:         if (isFace) {
1025:           if (debug) {PetscSynchronizedPrintf(PetscObjectComm((PetscObject) dm), "[%d]    Found shared face %D\n", rank, face);}
1026:           PetscHMapIJSet(facehash, key, l);
1027:           PetscSectionAddDof(candidateSection, localPoint, coneSize);
1028:         }
1029:       }
1030:     }
1031:     if (debug) {PetscSynchronizedFlush(PetscObjectComm((PetscObject) dm), NULL);}
1032:     PetscHMapIJClear(facehash);
1033:     PetscSectionSetUp(candidateSection);
1034:     PetscSectionGetStorageSize(candidateSection, &candidatesSize);
1035:     PetscMalloc1(candidatesSize, &candidates);
1036:     for (l = 0; l < numLeaves; ++l) {
1037:       const PetscInt    localPoint = localPoints[l];
1038:       const PetscInt   *support;
1039:       PetscInt          supportSize, s, offset, idx = 0;

1041:       if (debug) {PetscSynchronizedPrintf(PetscObjectComm((PetscObject) dm), "[%d]  Checking local point %D\n", rank, localPoint);}
1042:       PetscSectionGetOffset(candidateSection, localPoint, &offset);
1043:       DMPlexGetSupportSize(dm, localPoint, &supportSize);
1044:       DMPlexGetSupport(dm, localPoint, &support);
1045:       for (s = 0; s < supportSize; ++s) {
1046:         const PetscInt  face = support[s];
1047:         const PetscInt *cone;
1048:         PetscInt        coneSize, c, f, root;
1049:         PetscBool       isFace = PETSC_TRUE;

1051:         /* Only add face once */
1052:         if (debug) {PetscSynchronizedPrintf(PetscObjectComm((PetscObject) dm), "[%d]    Support point %D\n", rank, face);}
1053:         key.i = localPoint;
1054:         key.j = face;
1055:         PetscHMapIJGet(facehash, key, &f);
1056:         if (f >= 0) continue;
1057:         DMPlexGetConeSize(dm, face, &coneSize);
1058:         DMPlexGetCone(dm, face, &cone);
1059:         /* If a cone point does not map to leaves on any proc, then do not put face in SF */
1060:         for (c = 0; c < coneSize; ++c) {
1061:           if (debug) {PetscSynchronizedPrintf(PetscObjectComm((PetscObject) dm), "[%d]      Cone point %D\n", rank, cone[c]);}
1062:           PetscHMapIGet(leafhash, cone[c], &root);
1063:           if (!rootdegree[cone[c]] && (root < 0)) {isFace = PETSC_FALSE; break;}
1064:         }
1065:         if (isFace) {
1066:           if (debug) {PetscSynchronizedPrintf(PetscObjectComm((PetscObject) dm), "[%d]    Adding shared face %D at idx %D\n", rank, face, idx);}
1067:           PetscHMapIJSet(facehash, key, l);
1068:           candidates[offset+idx].rank    = -1;
1069:           candidates[offset+idx++].index = coneSize-1;
1070:           for (c = 0; c < coneSize; ++c) {
1071:             if (cone[c] == localPoint) continue;
1072:             if (rootdegree[cone[c]]) {
1073:               candidates[offset+idx].rank    = rank;
1074:               candidates[offset+idx++].index = cone[c];
1075:             } else {
1076:               PetscHMapIGet(leafhash, cone[c], &root);
1077:               if (root < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Cannot locate local point %D in SF", cone[c]);
1078:               candidates[offset+idx++] = remotePoints[root];
1079:             }
1080:           }
1081:         }
1082:       }
1083:     }
1084:     if (debug) {PetscSynchronizedFlush(PetscObjectComm((PetscObject) dm), NULL);}
1085:     PetscHMapIJDestroy(&facehash);
1086:     PetscObjectViewFromOptions((PetscObject) candidateSection, NULL, "-petscsection_interp_candidate_view");
1087:     SFNodeArrayViewFromOptions(PetscObjectComm((PetscObject) dm), "-petscsection_interp_candidate_view", "Candidates", NULL, candidatesSize, candidates);
1088:   }
1089:   /* Gather candidate section / array pair into the root partition via inverse(multi(pointSF)). */
1090:   /*   Note that this section is indexed by offsets into leaves, not by point number */
1091:   {
1092:     PetscSF   sfMulti, sfInverse, sfCandidates;
1093:     PetscInt *remoteOffsets;

1095:     PetscSFGetMultiSF(pointSF, &sfMulti);
1096:     PetscSFCreateInverseSF(sfMulti, &sfInverse);
1097:     PetscSectionCreate(PetscObjectComm((PetscObject) dm), &candidateSectionRemote);
1098:     PetscSFDistributeSection(sfInverse, candidateSection, &remoteOffsets, candidateSectionRemote);
1099:     PetscSFCreateSectionSF(sfInverse, candidateSection, remoteOffsets, candidateSectionRemote, &sfCandidates);
1100:     PetscSectionGetStorageSize(candidateSectionRemote, &candidatesRemoteSize);
1101:     PetscMalloc1(candidatesRemoteSize, &candidatesRemote);
1102:     PetscSFBcastBegin(sfCandidates, MPIU_2INT, candidates, candidatesRemote);
1103:     PetscSFBcastEnd(sfCandidates, MPIU_2INT, candidates, candidatesRemote);
1104:     PetscSFDestroy(&sfInverse);
1105:     PetscSFDestroy(&sfCandidates);
1106:     PetscFree(remoteOffsets);

1108:     PetscObjectViewFromOptions((PetscObject) candidateSectionRemote, NULL, "-petscsection_interp_candidate_remote_view");
1109:     SFNodeArrayViewFromOptions(PetscObjectComm((PetscObject) dm), "-petscsection_interp_candidate_remote_view", "Remote Candidates", NULL, candidatesRemoteSize, candidatesRemote);
1110:   }
1111:   /* */
1112:   {
1113:     PetscInt idx;
1114:     /* There is a section point for every leaf attached to a given root point */
1115:     for (r = 0, idx = 0; r < numRoots; ++r) {
1116:       PetscInt deg;
1117:       for (deg = 0; deg < rootdegree[r]; ++deg, ++idx) {
1118:         PetscInt offset, dof, d;

1120:         PetscSectionGetDof(candidateSectionRemote, idx, &dof);
1121:         PetscSectionGetOffset(candidateSectionRemote, idx, &offset);
1122:         for (d = 0; d < dof; ++d) {
1123:           const PetscInt  sizeInd   = offset+d;
1124:           const PetscInt  numPoints = candidatesRemote[sizeInd].index;
1125:           const PetscInt *join      = NULL;
1126:           PetscInt        points[1024], p, joinSize;

1128:           points[0] = r;
1129:           for (p = 0; p < numPoints; ++p) {
1130:             DMPlexMapToLocalPoint(roothash, localPoints, rank, candidatesRemote[offset+(++d)], &points[p+1]);
1131:             if (ierr) {d += numPoints-1 - p; break;} /* We got a point not in our overlap */
1132:             if (debug) {PetscSynchronizedPrintf(PetscObjectComm((PetscObject) dm), "[%d]  Checking local candidate %D\n", rank, points[p+1]);}
1133:           }
1134:           if (ierr) continue;
1135:           DMPlexGetJoin(dm, numPoints+1, points, &joinSize, &join);
1136:           if (joinSize == 1) {
1137:             if (debug) {PetscSynchronizedPrintf(PetscObjectComm((PetscObject) dm), "[%d]    Adding face %D at idx %D\n", rank, join[0], sizeInd);}
1138:             candidatesRemote[sizeInd].rank  = rank;
1139:             candidatesRemote[sizeInd].index = join[0];
1140:           }
1141:           DMPlexRestoreJoin(dm, numPoints+1, points, &joinSize, &join);
1142:         }
1143:       }
1144:     }
1145:     if (debug) {PetscSynchronizedFlush(PetscObjectComm((PetscObject) dm), NULL);}
1146:   }
1147:   /* Push claims back to receiver via the MultiSF and derive new pointSF mapping on receiver */
1148:   {
1149:     PetscSF         sfMulti, sfClaims, sfPointNew;
1150:     PetscSFNode    *remotePointsNew;
1151:     PetscHMapI      claimshash;
1152:     PetscInt       *remoteOffsets, *localPointsNew;
1153:     PetscInt        claimsSize, pStart, pEnd, root, numLocalNew, p, d;

1155:     PetscSFGetMultiSF(pointSF, &sfMulti);
1156:     PetscSectionCreate(PetscObjectComm((PetscObject) dm), &claimSection);
1157:     PetscSFDistributeSection(sfMulti, candidateSectionRemote, &remoteOffsets, claimSection);
1158:     PetscSFCreateSectionSF(sfMulti, candidateSectionRemote, remoteOffsets, claimSection, &sfClaims);
1159:     PetscSectionGetStorageSize(claimSection, &claimsSize);
1160:     PetscMalloc1(claimsSize, &claims);
1161:     PetscSFBcastBegin(sfClaims, MPIU_2INT, candidatesRemote, claims);
1162:     PetscSFBcastEnd(sfClaims, MPIU_2INT, candidatesRemote, claims);
1163:     PetscSFDestroy(&sfClaims);
1164:     PetscFree(remoteOffsets);
1165:     PetscObjectViewFromOptions((PetscObject) claimSection, NULL, "-petscsection_interp_claim_view");
1166:     SFNodeArrayViewFromOptions(PetscObjectComm((PetscObject) dm), "-petscsection_interp_claim_view", "Claims", NULL, claimsSize, claims);
1167:     /* Walk the original section of local supports and add an SF entry for each updated item */
1168:     PetscHMapICreate(&claimshash);
1169:     for (p = 0; p < numRoots; ++p) {
1170:       PetscInt dof, offset;

1172:       if (debug) {PetscSynchronizedPrintf(PetscObjectComm((PetscObject) dm), "[%d]  Checking root for claims %D\n", rank, p);}
1173:       PetscSectionGetDof(candidateSection, p, &dof);
1174:       PetscSectionGetOffset(candidateSection, p, &offset);
1175:       for (d = 0; d < dof;) {
1176:         if (claims[offset+d].rank >= 0) {
1177:           const PetscInt  faceInd   = offset+d;
1178:           const PetscInt  numPoints = candidates[faceInd].index;
1179:           const PetscInt *join      = NULL;
1180:           PetscInt        joinSize, points[1024], c;

1182:           if (debug) {PetscSynchronizedPrintf(PetscObjectComm((PetscObject) dm), "[%d]    Found claim for remote point (%D, %D)\n", rank, claims[faceInd].rank, claims[faceInd].index);}
1183:           points[0] = p;
1184:           if (debug) {PetscSynchronizedPrintf(PetscObjectComm((PetscObject) dm), "[%d]      point %D\n", rank, points[0]);}
1185:           for (c = 0, ++d; c < numPoints; ++c, ++d) {
1186:             key.i = candidates[offset+d].index;
1187:             key.j = candidates[offset+d].rank;
1188:             PetscHMapIJGet(roothash, key, &root);
1189:             points[c+1] = localPoints[root];
1190:             if (debug) {PetscSynchronizedPrintf(PetscObjectComm((PetscObject) dm), "[%d]      point %D\n", rank, points[c+1]);}
1191:           }
1192:           DMPlexGetJoin(dm, numPoints+1, points, &joinSize, &join);
1193:           if (joinSize == 1) {
1194:             if (debug) {PetscSynchronizedPrintf(PetscObjectComm((PetscObject) dm), "[%d]    Found local face %D\n", rank, join[0]);}
1195:             PetscHMapISet(claimshash, join[0], faceInd);
1196:           }
1197:           DMPlexRestoreJoin(dm, numPoints+1, points, &joinSize, &join);
1198:         } else d += claims[offset+d].index+1;
1199:       }
1200:     }
1201:     if (debug) {PetscSynchronizedFlush(PetscObjectComm((PetscObject) dm), NULL);}
1202:     /* Create new pointSF from hashed claims */
1203:     PetscHMapIGetSize(claimshash, &numLocalNew);
1204:     DMPlexGetChart(dm, &pStart, &pEnd);
1205:     PetscMalloc1(numLeaves + numLocalNew, &localPointsNew);
1206:     PetscMalloc1(numLeaves + numLocalNew, &remotePointsNew);
1207:     for (p = 0; p < numLeaves; ++p) {
1208:       localPointsNew[p] = localPoints[p];
1209:       remotePointsNew[p].index = remotePoints[p].index;
1210:       remotePointsNew[p].rank  = remotePoints[p].rank;
1211:     }
1212:     p = numLeaves;
1213:     PetscHMapIGetKeys(claimshash, &p, localPointsNew);
1214:     PetscSortInt(numLocalNew, &localPointsNew[numLeaves]);
1215:     for (p = numLeaves; p < numLeaves + numLocalNew; ++p) {
1216:       PetscInt offset;
1217:       PetscHMapIGet(claimshash, localPointsNew[p], &offset);
1218:       remotePointsNew[p] = claims[offset];
1219:     }
1220:     PetscSFCreate(PetscObjectComm((PetscObject) dm), &sfPointNew);
1221:     PetscSFSetGraph(sfPointNew, pEnd-pStart, numLeaves+numLocalNew, localPointsNew, PETSC_OWN_POINTER, remotePointsNew, PETSC_OWN_POINTER);
1222:     DMSetPointSF(dm, sfPointNew);
1223:     PetscSFDestroy(&sfPointNew);
1224:     PetscHMapIDestroy(&claimshash);
1225:   }
1226:   PetscHMapIDestroy(&leafhash);
1227:   PetscHMapIJDestroy(&roothash);
1228:   PetscSectionDestroy(&candidateSection);
1229:   PetscSectionDestroy(&candidateSectionRemote);
1230:   PetscSectionDestroy(&claimSection);
1231:   PetscFree(candidates);
1232:   PetscFree(candidatesRemote);
1233:   PetscFree(claims);
1234:   PetscLogEventEnd(DMPLEX_InterpolateSF,dm,0,0,0);
1235:   return(0);
1236: }

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

1241:   Collective on dm

1243:   Input Parameters:
1244: + dm - The DMPlex object with only cells and vertices
1245: - dmInt - The interpolated DM

1247:   Output Parameter:
1248: . dmInt - The complete DMPlex object

1250:   Level: intermediate

1252:   Notes:
1253:     It does not copy over the coordinates.

1255: .seealso: DMPlexUninterpolate(), DMPlexCreateFromCellList(), DMPlexCopyCoordinates()
1256: @*/
1257: PetscErrorCode DMPlexInterpolate(DM dm, DM *dmInt)
1258: {
1259:   DMPlexInterpolatedFlag interpolated;
1260:   DM             idm, odm = dm;
1261:   PetscSF        sfPoint;
1262:   PetscInt       depth, dim, d;
1263:   const char    *name;
1264:   PetscBool      flg=PETSC_TRUE;

1270:   PetscLogEventBegin(DMPLEX_Interpolate,dm,0,0,0);
1271:   DMPlexGetDepth(dm, &depth);
1272:   DMGetDimension(dm, &dim);
1273:   DMPlexIsInterpolated(dm, &interpolated);
1274:   if (interpolated == DMPLEX_INTERPOLATED_PARTIAL) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Not for partially interpolated meshes");
1275:   if (interpolated == DMPLEX_INTERPOLATED_FULL) {
1276:     PetscObjectReference((PetscObject) dm);
1277:     idm  = dm;
1278:   } else {
1279:     for (d = 1; d < dim; ++d) {
1280:       /* Create interpolated mesh */
1281:       DMCreate(PetscObjectComm((PetscObject)dm), &idm);
1282:       DMSetType(idm, DMPLEX);
1283:       DMSetDimension(idm, dim);
1284:       if (depth > 0) {
1285:         DMPlexInterpolateFaces_Internal(odm, 1, idm);
1286:         DMGetPointSF(odm, &sfPoint);
1287:         DMPlexInterpolatePointSF(idm, sfPoint);
1288:       }
1289:       if (odm != dm) {DMDestroy(&odm);}
1290:       odm = idm;
1291:     }
1292:     PetscObjectGetName((PetscObject) dm,  &name);
1293:     PetscObjectSetName((PetscObject) idm,  name);
1294:     DMPlexCopyCoordinates(dm, idm);
1295:     DMCopyLabels(dm, idm, PETSC_COPY_VALUES, PETSC_FALSE);
1296:     PetscOptionsGetBool(((PetscObject)dm)->options, ((PetscObject)dm)->prefix, "-dm_plex_interpolate_orient_interfaces", &flg, NULL);
1297:     if (flg) {DMPlexOrientInterface_Internal(idm);}
1298:   }
1299:   {
1300:     PetscBool            isper;
1301:     const PetscReal      *maxCell, *L;
1302:     const DMBoundaryType *bd;

1304:     DMGetPeriodicity(dm,&isper,&maxCell,&L,&bd);
1305:     DMSetPeriodicity(idm,isper,maxCell,L,bd);
1306:   }
1307:   /* This function makes the mesh fully interpolated on all ranks */
1308:   {
1309:     DM_Plex *plex = (DM_Plex *) dm->data;
1310:     plex->interpolated = plex->interpolatedCollective = DMPLEX_INTERPOLATED_FULL;
1311:   }
1312:   *dmInt = idm;
1313:   PetscLogEventEnd(DMPLEX_Interpolate,dm,0,0,0);
1314:   return(0);
1315: }

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

1320:   Collective on dmA

1322:   Input Parameter:
1323: . dmA - The DMPlex object with initial coordinates

1325:   Output Parameter:
1326: . dmB - The DMPlex object with copied coordinates

1328:   Level: intermediate

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

1332: .seealso: DMCopyLabels(), DMGetCoordinates(), DMGetCoordinatesLocal(), DMGetCoordinateDM(), DMGetCoordinateSection()
1333: @*/
1334: PetscErrorCode DMPlexCopyCoordinates(DM dmA, DM dmB)
1335: {
1336:   Vec            coordinatesA, coordinatesB;
1337:   VecType        vtype;
1338:   PetscSection   coordSectionA, coordSectionB;
1339:   PetscScalar   *coordsA, *coordsB;
1340:   PetscInt       spaceDim, Nf, vStartA, vStartB, vEndA, vEndB, coordSizeB, v, d;
1341:   PetscInt       cStartA, cEndA, cStartB, cEndB, cS, cE;
1342:   PetscBool      lc = PETSC_FALSE;

1348:   if (dmA == dmB) return(0);
1349:   DMPlexGetDepthStratum(dmA, 0, &vStartA, &vEndA);
1350:   DMPlexGetDepthStratum(dmB, 0, &vStartB, &vEndB);
1351:   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);
1352:   DMPlexGetHeightStratum(dmA, 0, &cStartA, &cEndA);
1353:   DMPlexGetHeightStratum(dmB, 0, &cStartB, &cEndB);
1354:   DMGetCoordinateSection(dmA, &coordSectionA);
1355:   DMGetCoordinateSection(dmB, &coordSectionB);
1356:   if (coordSectionA == coordSectionB) return(0);
1357:   PetscSectionGetNumFields(coordSectionA, &Nf);
1358:   if (!Nf) return(0);
1359:   if (Nf > 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "The number of coordinate fields must be 1, not %D", Nf);
1360:   if (!coordSectionB) {
1361:     PetscInt dim;

1363:     PetscSectionCreate(PetscObjectComm((PetscObject) coordSectionA), &coordSectionB);
1364:     DMGetCoordinateDim(dmA, &dim);
1365:     DMSetCoordinateSection(dmB, dim, coordSectionB);
1366:     PetscObjectDereference((PetscObject) coordSectionB);
1367:   }
1368:   PetscSectionSetNumFields(coordSectionB, 1);
1369:   PetscSectionGetFieldComponents(coordSectionA, 0, &spaceDim);
1370:   PetscSectionSetFieldComponents(coordSectionB, 0, spaceDim);
1371:   PetscSectionGetChart(coordSectionA, &cS, &cE);
1372:   if (cStartA <= cS && cS < cEndA) { /* localized coordinates */
1373:     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);
1374:     cS = cS - cStartA + cStartB;
1375:     cE = vEndB;
1376:     lc = PETSC_TRUE;
1377:   } else {
1378:     cS = vStartB;
1379:     cE = vEndB;
1380:   }
1381:   PetscSectionSetChart(coordSectionB, cS, cE);
1382:   for (v = vStartB; v < vEndB; ++v) {
1383:     PetscSectionSetDof(coordSectionB, v, spaceDim);
1384:     PetscSectionSetFieldDof(coordSectionB, v, 0, spaceDim);
1385:   }
1386:   if (lc) { /* localized coordinates */
1387:     PetscInt c;

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

1392:       PetscSectionGetDof(coordSectionA, c + cStartA, &dof);
1393:       PetscSectionSetDof(coordSectionB, c + cStartB, dof);
1394:       PetscSectionSetFieldDof(coordSectionB, c + cStartB, 0, dof);
1395:     }
1396:   }
1397:   PetscSectionSetUp(coordSectionB);
1398:   PetscSectionGetStorageSize(coordSectionB, &coordSizeB);
1399:   DMGetCoordinatesLocal(dmA, &coordinatesA);
1400:   VecCreate(PETSC_COMM_SELF, &coordinatesB);
1401:   PetscObjectSetName((PetscObject) coordinatesB, "coordinates");
1402:   VecSetSizes(coordinatesB, coordSizeB, PETSC_DETERMINE);
1403:   VecGetBlockSize(coordinatesA, &d);
1404:   VecSetBlockSize(coordinatesB, d);
1405:   VecGetType(coordinatesA, &vtype);
1406:   VecSetType(coordinatesB, vtype);
1407:   VecGetArray(coordinatesA, &coordsA);
1408:   VecGetArray(coordinatesB, &coordsB);
1409:   for (v = 0; v < vEndB-vStartB; ++v) {
1410:     PetscInt offA, offB;

1412:     PetscSectionGetOffset(coordSectionA, v + vStartA, &offA);
1413:     PetscSectionGetOffset(coordSectionB, v + vStartB, &offB);
1414:     for (d = 0; d < spaceDim; ++d) {
1415:       coordsB[offB+d] = coordsA[offA+d];
1416:     }
1417:   }
1418:   if (lc) { /* localized coordinates */
1419:     PetscInt c;

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

1424:       PetscSectionGetOffset(coordSectionA, c + cStartA, &offA);
1425:       PetscSectionGetOffset(coordSectionB, c + cStartB, &offB);
1426:       PetscSectionGetDof(coordSectionA, c + cStartA, &dof);
1427:       PetscArraycpy(coordsB + offB,coordsA + offA,dof);
1428:     }
1429:   }
1430:   VecRestoreArray(coordinatesA, &coordsA);
1431:   VecRestoreArray(coordinatesB, &coordsB);
1432:   DMSetCoordinatesLocal(dmB, coordinatesB);
1433:   VecDestroy(&coordinatesB);
1434:   return(0);
1435: }

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

1440:   Collective on dm

1442:   Input Parameter:
1443: . dm - The complete DMPlex object

1445:   Output Parameter:
1446: . dmUnint - The DMPlex object with only cells and vertices

1448:   Level: intermediate

1450:   Notes:
1451:     It does not copy over the coordinates.

1453: .seealso: DMPlexInterpolate(), DMPlexCreateFromCellList(), DMPlexCopyCoordinates()
1454: @*/
1455: PetscErrorCode DMPlexUninterpolate(DM dm, DM *dmUnint)
1456: {
1457:   DMPlexInterpolatedFlag interpolated;
1458:   DM             udm;
1459:   PetscInt       dim, vStart, vEnd, cStart, cEnd, cMax, c, maxConeSize = 0, *cone;

1465:   DMGetDimension(dm, &dim);
1466:   DMPlexIsInterpolated(dm, &interpolated);
1467:   if (interpolated == DMPLEX_INTERPOLATED_PARTIAL) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Not for partially interpolated meshes");
1468:   if (interpolated == DMPLEX_INTERPOLATED_NONE || dim <= 1) {
1469:     /* in case dim <= 1 just keep the DMPLEX_INTERPOLATED_FULL flag */
1470:     PetscObjectReference((PetscObject) dm);
1471:     *dmUnint = dm;
1472:     return(0);
1473:   }
1474:   DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
1475:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
1476:   DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);
1477:   DMCreate(PetscObjectComm((PetscObject) dm), &udm);
1478:   DMSetType(udm, DMPLEX);
1479:   DMSetDimension(udm, dim);
1480:   DMPlexSetChart(udm, cStart, vEnd);
1481:   for (c = cStart; c < cEnd; ++c) {
1482:     PetscInt *closure = NULL, closureSize, cl, coneSize = 0;

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

1488:       if ((p >= vStart) && (p < vEnd)) ++coneSize;
1489:     }
1490:     DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
1491:     DMPlexSetConeSize(udm, c, coneSize);
1492:     maxConeSize = PetscMax(maxConeSize, coneSize);
1493:   }
1494:   DMSetUp(udm);
1495:   PetscMalloc1(maxConeSize, &cone);
1496:   for (c = cStart; c < cEnd; ++c) {
1497:     PetscInt *closure = NULL, closureSize, cl, coneSize = 0;

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

1503:       if ((p >= vStart) && (p < vEnd)) cone[coneSize++] = p;
1504:     }
1505:     DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
1506:     DMPlexSetCone(udm, c, cone);
1507:   }
1508:   PetscFree(cone);
1509:   DMPlexSetHybridBounds(udm, cMax, PETSC_DETERMINE, PETSC_DETERMINE, PETSC_DETERMINE);
1510:   DMPlexSymmetrize(udm);
1511:   DMPlexStratify(udm);
1512:   /* Reduce SF */
1513:   {
1514:     PetscSF            sfPoint, sfPointUn;
1515:     const PetscSFNode *remotePoints;
1516:     const PetscInt    *localPoints;
1517:     PetscSFNode       *remotePointsUn;
1518:     PetscInt          *localPointsUn;
1519:     PetscInt           vEnd, numRoots, numLeaves, l;
1520:     PetscInt           numLeavesUn = 0, n = 0;
1521:     PetscErrorCode     ierr;

1523:     /* Get original SF information */
1524:     DMGetPointSF(dm, &sfPoint);
1525:     DMGetPointSF(udm, &sfPointUn);
1526:     DMPlexGetDepthStratum(dm, 0, NULL, &vEnd);
1527:     PetscSFGetGraph(sfPoint, &numRoots, &numLeaves, &localPoints, &remotePoints);
1528:     /* Allocate space for cells and vertices */
1529:     for (l = 0; l < numLeaves; ++l) if (localPoints[l] < vEnd) numLeavesUn++;
1530:     /* Fill in leaves */
1531:     if (vEnd >= 0) {
1532:       PetscMalloc1(numLeavesUn, &remotePointsUn);
1533:       PetscMalloc1(numLeavesUn, &localPointsUn);
1534:       for (l = 0; l < numLeaves; l++) {
1535:         if (localPoints[l] < vEnd) {
1536:           localPointsUn[n]        = localPoints[l];
1537:           remotePointsUn[n].rank  = remotePoints[l].rank;
1538:           remotePointsUn[n].index = remotePoints[l].index;
1539:           ++n;
1540:         }
1541:       }
1542:       if (n != numLeavesUn) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Inconsistent number of leaves %d != %d", n, numLeavesUn);
1543:       PetscSFSetGraph(sfPointUn, vEnd, numLeavesUn, localPointsUn, PETSC_OWN_POINTER, remotePointsUn, PETSC_OWN_POINTER);
1544:     }
1545:   }
1546:   {
1547:     PetscBool            isper;
1548:     const PetscReal      *maxCell, *L;
1549:     const DMBoundaryType *bd;

1551:     DMGetPeriodicity(dm,&isper,&maxCell,&L,&bd);
1552:     DMSetPeriodicity(udm,isper,maxCell,L,bd);
1553:   }
1554:   /* This function makes the mesh fully uninterpolated on all ranks */
1555:   {
1556:     DM_Plex *plex = (DM_Plex *) dm->data;
1557:     plex->interpolated = plex->interpolatedCollective = DMPLEX_INTERPOLATED_NONE;
1558:   }
1559:   *dmUnint = udm;
1560:   return(0);
1561: }

1563: static PetscErrorCode DMPlexIsInterpolated_Internal(DM dm, DMPlexInterpolatedFlag *interpolated)
1564: {
1565:   PetscInt       coneSize, depth, dim, h, p, pStart, pEnd;
1566:   MPI_Comm       comm;

1570:   PetscObjectGetComm((PetscObject)dm, &comm);
1571:   DMPlexGetDepth(dm, &depth);
1572:   DMGetDimension(dm, &dim);

1574:   if (depth == dim) {
1575:     *interpolated = DMPLEX_INTERPOLATED_FULL;
1576:     if (!dim) goto finish;

1578:     /* Check points at height = dim are vertices (have no cones) */
1579:     DMPlexGetHeightStratum(dm, dim, &pStart, &pEnd);
1580:     for (p=pStart; p<pEnd; p++) {
1581:       DMPlexGetConeSize(dm, p, &coneSize);
1582:       if (coneSize) {
1583:         *interpolated = DMPLEX_INTERPOLATED_PARTIAL;
1584:         goto finish;
1585:       }
1586:     }

1588:     /* Check points at height < dim have cones */
1589:     for (h=0; h<dim; h++) {
1590:       DMPlexGetHeightStratum(dm, h, &pStart, &pEnd);
1591:       for (p=pStart; p<pEnd; p++) {
1592:         DMPlexGetConeSize(dm, p, &coneSize);
1593:         if (!coneSize) {
1594:           *interpolated = DMPLEX_INTERPOLATED_PARTIAL;
1595:           goto finish;
1596:         }
1597:       }
1598:     }
1599:   } else if (depth == 1) {
1600:     *interpolated = DMPLEX_INTERPOLATED_NONE;
1601:   } else {
1602:     *interpolated = DMPLEX_INTERPOLATED_PARTIAL;
1603:   }
1604: finish:
1605:   return(0);
1606: }

1608: /*@
1609:   DMPlexIsInterpolated - Find out whether this DM is interpolated, i.e. number of strata is equal to dimension.

1611:   Not Collective

1613:   Input Parameter:
1614: . dm      - The DM object

1616:   Output Parameter:
1617: . interpolated - Flag whether the DM is interpolated

1619:   Level: intermediate

1621:   Notes:
1622:   This is NOT collective so the results can be different on different ranks in special cases.
1623:   However, DMPlexInterpolate() guarantees the result is the same on all.
1624:   Unlike DMPlexIsInterpolatedCollective(), this cannot return DMPLEX_INTERPOLATED_MIXED.

1626: .seealso: DMPlexInterpolate(), DMPlexIsInterpolatedCollective()
1627: @*/
1628: PetscErrorCode DMPlexIsInterpolated(DM dm, DMPlexInterpolatedFlag *interpolated)
1629: {
1630:   DM_Plex        *plex = (DM_Plex *) dm->data;
1631:   PetscErrorCode  ierr;

1636:   if (plex->interpolated < 0) {
1637:     DMPlexIsInterpolated_Internal(dm, &plex->interpolated);
1638:   } else {
1639: #if defined (PETSC_USE_DEBUG)
1640:     DMPlexInterpolatedFlag flg;

1642:     DMPlexIsInterpolated_Internal(dm, &flg);
1643:     if (flg != plex->interpolated) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "stashed DMPlexInterpolatedFlag is inconsistent");
1644: #endif
1645:   }
1646:   *interpolated = plex->interpolated;
1647:   return(0);
1648: }

1650: /*@
1651:   DMPlexIsInterpolatedCollective - Find out whether this DM is interpolated, i.e. number of strata is equal to dimension.

1653:   Collective

1655:   Input Parameter:
1656: . dm      - The DM object

1658:   Output Parameter:
1659: . interpolated - Flag whether the DM is interpolated

1661:   Level: intermediate

1663:   Notes:
1664:   This is collective so the results are always guaranteed to be the same on all ranks.
1665:   Unlike DMPlexIsInterpolated(), this will return DMPLEX_INTERPOLATED_MIXED if the results of DMPlexIsInterpolated() are different on different ranks.

1667: .seealso: DMPlexInterpolate(), DMPlexIsInterpolated()
1668: @*/
1669: PetscErrorCode DMPlexIsInterpolatedCollective(DM dm, DMPlexInterpolatedFlag *interpolated)
1670: {
1671:   DM_Plex        *plex = (DM_Plex *) dm->data;
1672:   PetscBool       debug=PETSC_FALSE;
1673:   PetscErrorCode  ierr;

1678:   PetscOptionsGetBool(((PetscObject) dm)->options, ((PetscObject) dm)->prefix, "-dm_plex_is_interpolated_collective_debug", &debug, NULL);
1679:   if (plex->interpolatedCollective < 0) {
1680:     DMPlexInterpolatedFlag  min, max;
1681:     MPI_Comm                comm;

1683:     PetscObjectGetComm((PetscObject)dm, &comm);
1684:     DMPlexIsInterpolated(dm, &plex->interpolatedCollective);
1685:     MPI_Allreduce(&plex->interpolatedCollective, &min, 1, MPIU_ENUM, MPI_MIN, comm);
1686:     MPI_Allreduce(&plex->interpolatedCollective, &max, 1, MPIU_ENUM, MPI_MAX, comm);
1687:     if (min != max) plex->interpolatedCollective = DMPLEX_INTERPOLATED_MIXED;
1688:     if (debug) {
1689:       PetscMPIInt rank;

1691:       MPI_Comm_rank(comm, &rank);
1692:       PetscSynchronizedPrintf(comm, "[%d] interpolated=%s interpolatedCollective=%s\n", rank, DMPlexInterpolatedFlags[plex->interpolated], DMPlexInterpolatedFlags[plex->interpolatedCollective]);
1693:       PetscSynchronizedFlush(comm, PETSC_STDOUT);
1694:     }
1695:   }
1696:   *interpolated = plex->interpolatedCollective;
1697:   return(0);
1698: }