Actual source code: plexinterpolate.c

petsc-3.12.1 2019-10-22
Report Typos and Errors
  1:  #include <petsc/private/dmpleximpl.h>
  2:  #include <petsc/private/hashmapi.h>
  3:  #include <petsc/private/hashmapij.h>

  5: /* HashIJKL */

  7:  #include <petsc/private/hashmap.h>

  9: typedef struct _PetscHashIJKLKey { PetscInt i, j, k, l; } PetscHashIJKLKey;

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

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

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


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

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

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

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

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

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

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

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

239: static PetscErrorCode DMPlexRestoreRawFacesHybrid_Internal(DM dm, PetscInt dim, PetscInt coneSize, const PetscInt cone[], PetscInt *numFaces, PetscInt *numFacesNotH, PetscInt *faceSize, const PetscInt *faces[])
240: {
241:   PetscErrorCode  ierr;

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

248: static PetscErrorCode DMPlexGetFacesHybrid_Internal(DM dm, PetscInt dim, PetscInt p, PetscInt *numFaces, PetscInt *numFacesNotH, PetscInt *faceSize, const PetscInt *faces[])
249: {
250:   const PetscInt *cone = NULL;
251:   PetscInt        coneSize;
252:   PetscErrorCode  ierr;

256:   DMPlexGetConeSize(dm, p, &coneSize);
257:   DMPlexGetCone(dm, p, &cone);
258:   DMPlexGetRawFacesHybrid_Internal(dm, dim, coneSize, cone, numFaces, numFacesNotH, faceSize, faces);
259:   return(0);
260: }

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

711: PetscErrorCode DMPlexOrientCell(DM dm, PetscInt p, PetscInt masterConeSize, const PetscInt masterCone[])
712: {
713:   PetscInt coneSize;
714:   PetscInt start1=0;
715:   PetscBool reverse1=PETSC_FALSE;
716:   const PetscInt *cone=NULL;

720:   DMPlexGetConeSize(dm, p, &coneSize);
721:   if (!coneSize) return(0); /* do nothing for points with no cone */
722:   DMPlexGetCone(dm, p, &cone);
723:   DMPlexFixFaceOrientations_Orient_Private(coneSize, masterConeSize, masterCone, cone, &start1, &reverse1);
724: #if defined(PETSC_USE_DEBUG)
725:   if (PetscUnlikely(cone[start1] != masterCone[0])) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "The algorithm above is wrong as cone[%d] = %d != %d = masterCone[0]", start1, cone[start1], masterCone[0]);
726: #endif
727:   DMPlexOrientCell_Internal(dm, p, start1, reverse1);
728: #if defined(PETSC_USE_DEBUG)
729:   {
730:     PetscInt c;
731:     DMPlexGetCone(dm, p, &cone);
732:     for (c = 0; c < 2; c++) {
733:       if (PetscUnlikely(cone[c] != masterCone[c])) SETERRQ4(PETSC_COMM_SELF, PETSC_ERR_PLIB, "The algorithm above is wrong as cone[%d] = %d != %d = masterCone[%d]", c, cone[c], masterCone[c], c);
734:     }
735:   }
736: #endif
737:   return(0);
738: }

740: PetscErrorCode DMPlexOrientCell_Internal(DM dm, PetscInt p, PetscInt start1, PetscBool reverse1)
741: {
742:   PetscInt i, j, k, maxConeSize, coneSize, coneConeSize, supportSize, supportConeSize;
743:   PetscInt start0, start;
744:   PetscBool reverse0, reverse;
745:   PetscInt newornt;
746:   const PetscInt *cone=NULL, *support=NULL, *supportCone=NULL, *ornts=NULL;
747:   PetscInt *newcone=NULL, *newornts=NULL;

751:   if (!start1 && !reverse1) return(0);
752:   DMPlexGetConeSize(dm, p, &coneSize);
753:   if (!coneSize) return(0); /* do nothing for points with no cone */
754:   DMPlexGetCone(dm, p, &cone);
755:   DMPlexGetMaxSizes(dm, &maxConeSize, NULL);
756:   /* permute p's cone and orientations */
757:   DMPlexGetConeOrientation(dm, p, &ornts);
758:   DMGetWorkArray(dm, maxConeSize, MPIU_INT, &newcone);
759:   DMGetWorkArray(dm, maxConeSize, MPIU_INT, &newornts);
760:   DMPlexFixFaceOrientations_Permute_Private(coneSize, cone, start1, reverse1, newcone);
761:   DMPlexFixFaceOrientations_Permute_Private(coneSize, ornts, start1, reverse1, newornts);
762:   /* if direction of p (face) is flipped, flip also p's cone points (edges) */
763:   if (reverse1) {
764:     for (i=0; i<coneSize; i++) {
765:       DMPlexGetConeSize(dm, cone[i], &coneConeSize);
766:       DMPlexFixFaceOrientations_Translate_Private(newornts[i], &start0, &reverse0);
767:       DMPlexFixFaceOrientations_Combine_Private(coneConeSize, start0, reverse0, 1, PETSC_FALSE, &start, &reverse);
768:       DMPlexFixFaceOrientations_TranslateBack_Private(coneConeSize, start, reverse, &newornts[i]);
769:     }
770:   }
771:   DMPlexSetConeOrientation(dm, p, newornts);
772:   /* fix oriention of p within cones of p's support points */
773:   DMPlexGetSupport(dm, p, &support);
774:   DMPlexGetSupportSize(dm, p, &supportSize);
775:   for (j=0; j<supportSize; j++) {
776:     DMPlexGetCone(dm, support[j], &supportCone);
777:     DMPlexGetConeSize(dm, support[j], &supportConeSize);
778:     DMPlexGetConeOrientation(dm, support[j], &ornts);
779:     for (k=0; k<supportConeSize; k++) {
780:       if (supportCone[k] != p) continue;
781:       DMPlexFixFaceOrientations_Translate_Private(ornts[k], &start0, &reverse0);
782:       DMPlexFixFaceOrientations_Combine_Private(coneSize, start0, reverse0, start1, reverse1, &start, &reverse);
783:       DMPlexFixFaceOrientations_TranslateBack_Private(coneSize, start, reverse, &newornt);
784:       DMPlexInsertConeOrientation(dm, support[j], k, newornt);
785:     }
786:   }
787:   /* rewrite cone */
788:   DMPlexSetCone(dm, p, newcone);
789:   DMRestoreWorkArray(dm, maxConeSize, MPIU_INT, &newcone);
790:   DMRestoreWorkArray(dm, maxConeSize, MPIU_INT, &newornts);
791:   return(0);
792: }

794: static PetscErrorCode SortRmineRremoteByRemote_Private(PetscSF sf, PetscInt *rmine1[], PetscInt *rremote1[])
795: {
796:   PetscInt            nleaves;
797:   PetscInt            nranks;
798:   const PetscMPIInt  *ranks=NULL;
799:   const PetscInt     *roffset=NULL, *rmine=NULL, *rremote=NULL;
800:   PetscInt            n, o, r;
801:   PetscErrorCode      ierr;

804:   PetscSFGetRootRanks(sf, &nranks, &ranks, &roffset, &rmine, &rremote);
805:   nleaves = roffset[nranks];
806:   PetscMalloc2(nleaves, rmine1, nleaves, rremote1);
807:   for (r=0; r<nranks; r++) {
808:     /* simultaneously sort rank-wise portions of rmine & rremote by values in rremote
809:        - to unify order with the other side */
810:     o = roffset[r];
811:     n = roffset[r+1] - o;
812:     PetscArraycpy(&(*rmine1)[o], &rmine[o], n);
813:     PetscArraycpy(&(*rremote1)[o], &rremote[o], n);
814:     PetscSortIntWithArray(n, &(*rremote1)[o], &(*rmine1)[o]);
815:   }
816:   return(0);
817: }

819: PetscErrorCode DMPlexOrientInterface(DM dm)
820: {
821:   PetscSF           sf=NULL;
822:   PetscInt          (*roots)[2], (*leaves)[2];
823:   PetscMPIInt       (*rootsRanks)[2], (*leavesRanks)[2];
824:   const PetscInt    *locals=NULL;
825:   const PetscSFNode *remotes=NULL;
826:   PetscInt           nroots, nleaves, p, c;
827:   PetscInt           nranks, n, o, r;
828:   const PetscMPIInt *ranks=NULL;
829:   const PetscInt    *roffset=NULL;
830:   PetscInt          *rmine1=NULL, *rremote1=NULL; /* rmine and rremote copies simultaneously sorted by rank and rremote */
831:   const PetscInt    *cone=NULL;
832:   PetscInt           coneSize, ind0;
833:   MPI_Comm           comm;
834:   PetscMPIInt        rank;
835:   PetscInt           debug = 0;
836:   PetscErrorCode     ierr;

839:   DMGetPointSF(dm, &sf);
840:   PetscSFGetGraph(sf, &nroots, &nleaves, &locals, &remotes);
841:   if (nroots < 0) return(0);
842:   PetscSFSetUp(sf);
843:   PetscSFGetRootRanks(sf, &nranks, &ranks, &roffset, NULL, NULL);
844: #if defined(PETSC_USE_DEBUG)
845:   DMViewFromOptions(dm, NULL, "-before_fix_dm_view");
846:   DMPlexCheckPointSF(dm);
847: #endif
848:   SortRmineRremoteByRemote_Private(sf, &rmine1, &rremote1);
849:   PetscMalloc4(nroots, &roots, nroots, &leaves, nroots, &rootsRanks, nroots, &leavesRanks);
850:   PetscObjectGetComm((PetscObject) dm, &comm);
851:   MPI_Comm_rank(comm, &rank);
852:   if (debug && rank == 0) {PetscSynchronizedPrintf(comm, "Roots\n");}
853:   for (p = 0; p < nroots; ++p) {
854:     DMPlexGetConeSize(dm, p, &coneSize);
855:     DMPlexGetCone(dm, p, &cone);
856:     /* Translate all points to root numbering */
857:     for (c = 0; c < 2; c++) {
858:       if (coneSize > 1) {
859:         PetscFindInt(cone[c], nleaves, locals, &ind0);
860:         if (ind0 < 0) {
861:           roots[p][c] = cone[c];
862:           rootsRanks[p][c] = rank;
863:         } else {
864:           roots[p][c] = remotes[ind0].index;
865:           rootsRanks[p][c] = remotes[ind0].rank;
866:         }
867:       } else {
868:         roots[p][c] = -1;
869:         rootsRanks[p][c] = -1;
870:       }
871:     }
872:   }
873:   if (debug) {
874:     for (p = 0; p < nroots; ++p) {
875:       DMPlexGetConeSize(dm, p, &coneSize);
876:       DMPlexGetCone(dm, p, &cone);
877:       if (coneSize > 1) {
878:         PetscSynchronizedPrintf(comm, "[%d]  %D: cone=[%D %D] roots=[%D %D] rootsRanks=[%D %D]\n", rank, p, cone[0], cone[1], roots[p][0], roots[p][1], rootsRanks[p][0], rootsRanks[p][1]);
879:       }
880:     }
881:     PetscSynchronizedFlush(comm, NULL);
882:   }
883:   for (p = 0; p < nroots; ++p) {
884:     for (c = 0; c < 2; c++) {
885:       leaves[p][c] = -2;
886:       leavesRanks[p][c] = -2;
887:     }
888:   }
889:   PetscSFBcastBegin(sf, MPIU_2INT, roots, leaves);
890:   PetscSFBcastBegin(sf, MPI_2INT, rootsRanks, leavesRanks);
891:   PetscSFBcastEnd(sf, MPIU_2INT, roots, leaves);
892:   PetscSFBcastEnd(sf, MPI_2INT, rootsRanks, leavesRanks);
893:   if (debug) {PetscSynchronizedFlush(comm, NULL);}
894:   if (debug && rank == 0) {PetscSynchronizedPrintf(comm, "Referred leaves\n");}
895:   for (p = 0; p < nroots; ++p) {
896:     if (leaves[p][0] < 0) continue;
897:     DMPlexGetConeSize(dm, p, &coneSize);
898:     DMPlexGetCone(dm, p, &cone);
899:     if (debug) {PetscSynchronizedPrintf(comm, "[%d]  %D: cone=[%D %D] leaves=[%D %D] roots=[%D %D] leavesRanks=[%D %D] rootsRanks=[%D %D]\n", rank, p, cone[0], cone[1], leaves[p][0], leaves[p][1], roots[p][0], roots[p][1], leavesRanks[p][0], leavesRanks[p][1], rootsRanks[p][0], rootsRanks[p][1]);}
900:     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])) {
901:       PetscInt masterCone[2];
902:       /* Translate these two cone points back to leave numbering */
903:       for (c = 0; c < 2; c++) {
904:         if (leavesRanks[p][c] == rank) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "this should never happen - remote rank of point %D is the same rank",leavesRanks[p][c]);
905:         /* Find index of rank leavesRanks[p][c] among remote ranks */
906:         /* No need for PetscMPIIntCast because these integers were originally cast from PetscMPIInt. */
907:         PetscFindMPIInt((PetscMPIInt)leavesRanks[p][c], nranks, ranks, &r);
908:         if (PetscUnlikely(r < 0)) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "this should never happen - rank %D not found among remote ranks",leavesRanks[p][c]);
909:         /* Find point leaves[p][c] among remote points aimed at rank leavesRanks[p][c] */
910:         o = roffset[r];
911:         n = roffset[r+1] - o;
912:         PetscFindInt(leaves[p][c], n, &rremote1[o], &ind0);
913:         if (PetscUnlikely(ind0 < 0)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "No cone point of %D is connected to (%D, %D) - it seems there is missing connection in point SF!",p,ranks[r],leaves[p][c]);
914:         /* Get the corresponding local point */
915:         masterCone[c] = rmine1[o+ind0];
916:       }
917:       if (debug) {PetscSynchronizedPrintf(comm, "[%d]  %D: masterCone=[%D %D]\n", rank, p, masterCone[0], masterCone[1]);}
918:       /* Vaclav's note: Here we only compare first 2 points of the cone. Full cone size would lead to stronger self-checking. */
919:       DMPlexOrientCell(dm, p, 2, masterCone);
920:     }
921:   }
922: #if defined(PETSC_USE_DEBUG)
923:   DMViewFromOptions(dm, NULL, "-after_fix_dm_view");
924:   for (r = 0; r < nleaves; ++r) {
925:     p = locals[r];
926:     DMPlexGetConeSize(dm, p, &coneSize);
927:     if (!coneSize) continue;
928:     DMPlexGetCone(dm, p, &cone);
929:     for (c = 0; c < 2; c++) {
930:       PetscFindInt(cone[c], nleaves, locals, &ind0);
931:       if (ind0 < 0) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Point SF contains %D but is missing its cone point cone[%D] = %D!", p, c, cone[c]);
932:       if (leaves[p][c] != remotes[ind0].index || leavesRanks[p][c] != remotes[ind0].rank) {
933:         if (leavesRanks[p][c] == rank) {
934:           PetscInt ind1;
935:           PetscFindInt(leaves[p][c], nleaves, locals, &ind1);
936:           if (ind1 < 0) {
937:             SETERRQ8(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Point %D = locals[%d]: cone[%D]=%D --> (%D, %D) differs from the enforced (%D, %D). The latter was not even found among the local SF points - it is probably broken!", p, r, c, cone[c], remotes[ind0].rank, remotes[ind0].index, leavesRanks[p][c], leaves[p][c]);
938:           } else {
939:             SETERRQ9(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Point %D = locals[%d]: cone[%D]=%D --> (%D, %D) differs from the enforced %D --> (%D, %D). Is the algorithm above or the point SF broken?", p, r, c, cone[c], remotes[ind0].rank, remotes[ind0].index, leaves[p][c], remotes[ind1].rank, remotes[ind1].index);
940:           }
941:         } else {
942:           SETERRQ8(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Point %D = locals[%d]: cone[%D]=%D --> (%D, %D) differs from the enforced (%D, %D). Is the algorithm above or the point SF broken?", p, r, c, cone[c], remotes[ind0].rank, remotes[ind0].index, leavesRanks[p][c], leaves[p][c]);
943:         }
944:       }
945:     }
946:   }
947: #endif
948:   if (debug) {PetscSynchronizedFlush(comm, NULL);}
949:   PetscFree4(roots, leaves, rootsRanks, leavesRanks);
950:   PetscFree2(rmine1, rremote1);
951:   return(0);
952: }

954: static PetscErrorCode IntArrayViewFromOptions(MPI_Comm comm, const char opt[], const char name[], const char idxname[], const char valname[], PetscInt n, const PetscInt a[])
955: {
956:   PetscInt       idx;
957:   PetscMPIInt    rank;
958:   PetscBool      flg;

962:   PetscOptionsHasName(NULL, NULL, opt, &flg);
963:   if (!flg) return(0);
964:   MPI_Comm_rank(comm, &rank);
965:   PetscSynchronizedPrintf(comm, "[%d]%s:\n", rank, name);
966:   for (idx = 0; idx < n; ++idx) {PetscSynchronizedPrintf(comm, "[%d]%s %D %s %D\n", rank, idxname, idx, valname, a[idx]);}
967:   PetscSynchronizedFlush(comm, NULL);
968:   return(0);
969: }

971: static PetscErrorCode SFNodeArrayViewFromOptions(MPI_Comm comm, const char opt[], const char name[], const char idxname[], PetscInt n, const PetscSFNode a[])
972: {
973:   PetscInt       idx;
974:   PetscMPIInt    rank;
975:   PetscBool      flg;

979:   PetscOptionsHasName(NULL, NULL, opt, &flg);
980:   if (!flg) return(0);
981:   MPI_Comm_rank(comm, &rank);
982:   PetscSynchronizedPrintf(comm, "[%d]%s:\n", rank, name);
983:   if (idxname) {
984:     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);}
985:   } else {
986:     for (idx = 0; idx < n; ++idx) {PetscSynchronizedPrintf(comm, "[%d]rank %D index %D\n", rank, a[idx].rank, a[idx].index);}
987:   }
988:   PetscSynchronizedFlush(comm, NULL);
989:   return(0);
990: }

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

997:   if (remotePoint.rank == rank) {
998:     *localPoint = remotePoint.index;
999:   } else {
1000:     PetscHashIJKey key;
1001:     PetscInt       root;

1003:     key.i = remotePoint.index;
1004:     key.j = remotePoint.rank;
1005:     PetscHMapIJGet(roothash, key, &root);
1006:     if (root >= 0) {
1007:       *localPoint = localPoints[root];
1008:     } else PetscFunctionReturn(1);
1009:   }
1010:   return(0);
1011: }

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

1016:   Collective on dm

1018:   Input Parameters:
1019: + dm      - The interpolated DM
1020: - pointSF - The initial SF without interpolated points

1022:   Output Parameter:
1023: . pointSF - The SF including interpolated points

1025:   Level: intermediate

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

1029: .seealso: DMPlexInterpolate(), DMPlexUninterpolate()
1030: @*/
1031: PetscErrorCode DMPlexInterpolatePointSF(DM dm, PetscSF pointSF)
1032: {
1033:   /*
1034:        Okay, the algorithm is:
1035:          - Take each point in the overlap (root)
1036:          - Look at the neighboring points in the overlap (candidates)
1037:          - Send these candidate points to neighbors
1038:          - Neighbor checks for edge between root and candidate
1039:          - If edge is found, it replaces candidate point with edge point
1040:          - Send back the overwritten candidates (claims)
1041:          - Original guy checks for edges, different from original candidate, and gets its own edge
1042:          - This pair is put into SF

1044:        We need a new algorithm that tolerates groups larger than 2.
1045:          - Take each point in the overlap (root)
1046:          - Find all collections of points in the overlap which make faces (do early join)
1047:          - Send collections as candidates (add size as first number)
1048:            - Make sure to send collection to all owners of all overlap points in collection
1049:          - Neighbor check for face in collections
1050:          - If face is found, it replaces candidate point with face point
1051:          - Send back the overwritten candidates (claims)
1052:          - Original guy checks for faces, different from original candidate, and gets its own face
1053:          - This pair is put into SF
1054:   */
1055:   PetscHMapI         leafhash;
1056:   PetscHMapIJ        roothash;
1057:   const PetscInt    *localPoints, *rootdegree;
1058:   const PetscSFNode *remotePoints;
1059:   PetscSFNode       *candidates, *candidatesRemote, *claims;
1060:   PetscSection       candidateSection, candidateSectionRemote, claimSection;
1061:   PetscInt           numLeaves, l, numRoots, r, candidatesSize, candidatesRemoteSize;
1062:   PetscMPIInt        size, rank;
1063:   PetscHashIJKey     key;
1064:   PetscBool          debug = PETSC_FALSE;
1065:   PetscErrorCode     ierr;

1068:   PetscOptionsHasName(NULL, ((PetscObject) dm)->prefix, "-dmplex_interp_debug", &debug);
1069:   MPI_Comm_size(PetscObjectComm((PetscObject) dm), &size);
1070:   MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank);
1071:   PetscSFGetGraph(pointSF, &numRoots, &numLeaves, &localPoints, &remotePoints);
1072:   if (size < 2 || numRoots < 0) return(0);
1073:   DMPlexGetOverlap(dm, &r);
1074:   if (r) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Interpolation of overlapped DMPlex not implemented yet");
1075:   PetscObjectViewFromOptions((PetscObject) dm, NULL, "-dm_interp_pre_view");
1076:   PetscObjectViewFromOptions((PetscObject) pointSF, NULL, "-petscsf_interp_pre_view");
1077:   PetscLogEventBegin(DMPLEX_InterpolateSF,dm,0,0,0);
1078:   /* Build hashes of points in the SF for efficient lookup */
1079:   PetscHMapICreate(&leafhash);
1080:   PetscHMapIJCreate(&roothash);
1081:   for (l = 0; l < numLeaves; ++l) {
1082:     key.i = remotePoints[l].index;
1083:     key.j = remotePoints[l].rank;
1084:     PetscHMapISet(leafhash, localPoints[l], l);
1085:     PetscHMapIJSet(roothash, key, l);
1086:   }
1087:   /* Compute root degree to identify shared points */
1088:   PetscSFComputeDegreeBegin(pointSF, &rootdegree);
1089:   PetscSFComputeDegreeEnd(pointSF, &rootdegree);
1090:   IntArrayViewFromOptions(PetscObjectComm((PetscObject) dm), "-interp_root_degree_view", "Root degree", "point", "degree", numRoots, rootdegree);
1091:   /* Build a section / SFNode array of candidate points (face bd points) in the cone(support(leaf)),
1092:      where each candidate is defined by a set of remote points (roots) for the other points that define the face. */
1093:   PetscSectionCreate(PetscObjectComm((PetscObject) dm), &candidateSection);
1094:   PetscSectionSetChart(candidateSection, 0, numRoots);
1095:   {
1096:     PetscHMapIJ facehash;

1098:     PetscHMapIJCreate(&facehash);
1099:     for (l = 0; l < numLeaves; ++l) {
1100:       const PetscInt    localPoint = localPoints[l];
1101:       const PetscInt   *support;
1102:       PetscInt          supportSize, s;

1104:       if (debug) {PetscSynchronizedPrintf(PetscObjectComm((PetscObject) dm), "[%d]  Checking local point %D\n", rank, localPoint);}
1105:       DMPlexGetSupportSize(dm, localPoint, &supportSize);
1106:       DMPlexGetSupport(dm, localPoint, &support);
1107:       for (s = 0; s < supportSize; ++s) {
1108:         const PetscInt  face = support[s];
1109:         const PetscInt *cone;
1110:         PetscInt        coneSize, c, f, root;
1111:         PetscBool       isFace = PETSC_TRUE;

1113:         /* Only add face once */
1114:         if (debug) {PetscSynchronizedPrintf(PetscObjectComm((PetscObject) dm), "[%d]    Support point %D\n", rank, face);}
1115:         key.i = localPoint;
1116:         key.j = face;
1117:         PetscHMapIJGet(facehash, key, &f);
1118:         if (f >= 0) continue;
1119:         DMPlexGetConeSize(dm, face, &coneSize);
1120:         DMPlexGetCone(dm, face, &cone);
1121:         /* If a cone point does not map to leaves on any proc, then do not put face in SF */
1122:         for (c = 0; c < coneSize; ++c) {
1123:           if (debug) {PetscSynchronizedPrintf(PetscObjectComm((PetscObject) dm), "[%d]      Cone point %D\n", rank, cone[c]);}
1124:           PetscHMapIGet(leafhash, cone[c], &root);
1125:           if (!rootdegree[cone[c]] && (root < 0)) {isFace = PETSC_FALSE; break;}
1126:         }
1127:         if (isFace) {
1128:           if (debug) {PetscSynchronizedPrintf(PetscObjectComm((PetscObject) dm), "[%d]    Found shared face %D\n", rank, face);}
1129:           PetscHMapIJSet(facehash, key, l);
1130:           PetscSectionAddDof(candidateSection, localPoint, coneSize);
1131:         }
1132:       }
1133:     }
1134:     if (debug) {PetscSynchronizedFlush(PetscObjectComm((PetscObject) dm), NULL);}
1135:     PetscHMapIJClear(facehash);
1136:     PetscSectionSetUp(candidateSection);
1137:     PetscSectionGetStorageSize(candidateSection, &candidatesSize);
1138:     PetscMalloc1(candidatesSize, &candidates);
1139:     for (l = 0; l < numLeaves; ++l) {
1140:       const PetscInt    localPoint = localPoints[l];
1141:       const PetscInt   *support;
1142:       PetscInt          supportSize, s, offset, idx = 0;

1144:       if (debug) {PetscSynchronizedPrintf(PetscObjectComm((PetscObject) dm), "[%d]  Checking local point %D\n", rank, localPoint);}
1145:       PetscSectionGetOffset(candidateSection, localPoint, &offset);
1146:       DMPlexGetSupportSize(dm, localPoint, &supportSize);
1147:       DMPlexGetSupport(dm, localPoint, &support);
1148:       for (s = 0; s < supportSize; ++s) {
1149:         const PetscInt  face = support[s];
1150:         const PetscInt *cone;
1151:         PetscInt        coneSize, c, f, root;
1152:         PetscBool       isFace = PETSC_TRUE;

1154:         /* Only add face once */
1155:         if (debug) {PetscSynchronizedPrintf(PetscObjectComm((PetscObject) dm), "[%d]    Support point %D\n", rank, face);}
1156:         key.i = localPoint;
1157:         key.j = face;
1158:         PetscHMapIJGet(facehash, key, &f);
1159:         if (f >= 0) continue;
1160:         DMPlexGetConeSize(dm, face, &coneSize);
1161:         DMPlexGetCone(dm, face, &cone);
1162:         /* If a cone point does not map to leaves on any proc, then do not put face in SF */
1163:         for (c = 0; c < coneSize; ++c) {
1164:           if (debug) {PetscSynchronizedPrintf(PetscObjectComm((PetscObject) dm), "[%d]      Cone point %D\n", rank, cone[c]);}
1165:           PetscHMapIGet(leafhash, cone[c], &root);
1166:           if (!rootdegree[cone[c]] && (root < 0)) {isFace = PETSC_FALSE; break;}
1167:         }
1168:         if (isFace) {
1169:           if (debug) {PetscSynchronizedPrintf(PetscObjectComm((PetscObject) dm), "[%d]    Adding shared face %D at idx %D\n", rank, face, idx);}
1170:           PetscHMapIJSet(facehash, key, l);
1171:           candidates[offset+idx].rank    = -1;
1172:           candidates[offset+idx++].index = coneSize-1;
1173:           for (c = 0; c < coneSize; ++c) {
1174:             if (cone[c] == localPoint) continue;
1175:             if (rootdegree[cone[c]]) {
1176:               candidates[offset+idx].rank    = rank;
1177:               candidates[offset+idx++].index = cone[c];
1178:             } else {
1179:               PetscHMapIGet(leafhash, cone[c], &root);
1180:               if (root < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Cannot locate local point %D in SF", cone[c]);
1181:               candidates[offset+idx++] = remotePoints[root];
1182:             }
1183:           }
1184:         }
1185:       }
1186:     }
1187:     if (debug) {PetscSynchronizedFlush(PetscObjectComm((PetscObject) dm), NULL);}
1188:     PetscHMapIJDestroy(&facehash);
1189:     PetscObjectViewFromOptions((PetscObject) candidateSection, NULL, "-petscsection_interp_candidate_view");
1190:     SFNodeArrayViewFromOptions(PetscObjectComm((PetscObject) dm), "-petscsection_interp_candidate_view", "Candidates", NULL, candidatesSize, candidates);
1191:   }
1192:   /* Gather candidate section / array pair into the root partition via inverse(multi(pointSF)). */
1193:   /*   Note that this section is indexed by offsets into leaves, not by point number */
1194:   {
1195:     PetscSF   sfMulti, sfInverse, sfCandidates;
1196:     PetscInt *remoteOffsets;

1198:     PetscSFGetMultiSF(pointSF, &sfMulti);
1199:     PetscSFCreateInverseSF(sfMulti, &sfInverse);
1200:     PetscSectionCreate(PetscObjectComm((PetscObject) dm), &candidateSectionRemote);
1201:     PetscSFDistributeSection(sfInverse, candidateSection, &remoteOffsets, candidateSectionRemote);
1202:     PetscSFCreateSectionSF(sfInverse, candidateSection, remoteOffsets, candidateSectionRemote, &sfCandidates);
1203:     PetscSectionGetStorageSize(candidateSectionRemote, &candidatesRemoteSize);
1204:     PetscMalloc1(candidatesRemoteSize, &candidatesRemote);
1205:     PetscSFBcastBegin(sfCandidates, MPIU_2INT, candidates, candidatesRemote);
1206:     PetscSFBcastEnd(sfCandidates, MPIU_2INT, candidates, candidatesRemote);
1207:     PetscSFDestroy(&sfInverse);
1208:     PetscSFDestroy(&sfCandidates);
1209:     PetscFree(remoteOffsets);

1211:     PetscObjectViewFromOptions((PetscObject) candidateSectionRemote, NULL, "-petscsection_interp_candidate_remote_view");
1212:     SFNodeArrayViewFromOptions(PetscObjectComm((PetscObject) dm), "-petscsection_interp_candidate_remote_view", "Remote Candidates", NULL, candidatesRemoteSize, candidatesRemote);
1213:   }
1214:   /* */
1215:   {
1216:     PetscInt idx;
1217:     /* There is a section point for every leaf attached to a given root point */
1218:     for (r = 0, idx = 0; r < numRoots; ++r) {
1219:       PetscInt deg;
1220:       for (deg = 0; deg < rootdegree[r]; ++deg, ++idx) {
1221:         PetscInt offset, dof, d;

1223:         PetscSectionGetDof(candidateSectionRemote, idx, &dof);
1224:         PetscSectionGetOffset(candidateSectionRemote, idx, &offset);
1225:         for (d = 0; d < dof; ++d) {
1226:           const PetscInt  sizeInd   = offset+d;
1227:           const PetscInt  numPoints = candidatesRemote[sizeInd].index;
1228:           const PetscInt *join      = NULL;
1229:           PetscInt        points[1024], p, joinSize;

1231:           points[0] = r;
1232:           for (p = 0; p < numPoints; ++p) {
1233:             DMPlexMapToLocalPoint(roothash, localPoints, rank, candidatesRemote[offset+(++d)], &points[p+1]);
1234:             if (ierr) {d += numPoints-1 - p; break;} /* We got a point not in our overlap */
1235:             if (debug) {PetscSynchronizedPrintf(PetscObjectComm((PetscObject) dm), "[%d]  Checking local candidate %D\n", rank, points[p+1]);}
1236:           }
1237:           if (ierr) continue;
1238:           DMPlexGetJoin(dm, numPoints+1, points, &joinSize, &join);
1239:           if (joinSize == 1) {
1240:             if (debug) {PetscSynchronizedPrintf(PetscObjectComm((PetscObject) dm), "[%d]    Adding face %D at idx %D\n", rank, join[0], sizeInd);}
1241:             candidatesRemote[sizeInd].rank  = rank;
1242:             candidatesRemote[sizeInd].index = join[0];
1243:           }
1244:           DMPlexRestoreJoin(dm, numPoints+1, points, &joinSize, &join);
1245:         }
1246:       }
1247:     }
1248:     if (debug) {PetscSynchronizedFlush(PetscObjectComm((PetscObject) dm), NULL);}
1249:   }
1250:   /* Push claims back to receiver via the MultiSF and derive new pointSF mapping on receiver */
1251:   {
1252:     PetscSF         sfMulti, sfClaims, sfPointNew;
1253:     PetscSFNode    *remotePointsNew;
1254:     PetscHMapI      claimshash;
1255:     PetscInt       *remoteOffsets, *localPointsNew;
1256:     PetscInt        claimsSize, pStart, pEnd, root, numLocalNew, p, d;

1258:     PetscSFGetMultiSF(pointSF, &sfMulti);
1259:     PetscSectionCreate(PetscObjectComm((PetscObject) dm), &claimSection);
1260:     PetscSFDistributeSection(sfMulti, candidateSectionRemote, &remoteOffsets, claimSection);
1261:     PetscSFCreateSectionSF(sfMulti, candidateSectionRemote, remoteOffsets, claimSection, &sfClaims);
1262:     PetscSectionGetStorageSize(claimSection, &claimsSize);
1263:     PetscMalloc1(claimsSize, &claims);
1264:     PetscSFBcastBegin(sfClaims, MPIU_2INT, candidatesRemote, claims);
1265:     PetscSFBcastEnd(sfClaims, MPIU_2INT, candidatesRemote, claims);
1266:     PetscSFDestroy(&sfClaims);
1267:     PetscFree(remoteOffsets);
1268:     PetscObjectViewFromOptions((PetscObject) claimSection, NULL, "-petscsection_interp_claim_view");
1269:     SFNodeArrayViewFromOptions(PetscObjectComm((PetscObject) dm), "-petscsection_interp_claim_view", "Claims", NULL, claimsSize, claims);
1270:     /* Walk the original section of local supports and add an SF entry for each updated item */
1271:     PetscHMapICreate(&claimshash);
1272:     for (p = 0; p < numRoots; ++p) {
1273:       PetscInt dof, offset;

1275:       if (debug) {PetscSynchronizedPrintf(PetscObjectComm((PetscObject) dm), "[%d]  Checking root for claims %D\n", rank, p);}
1276:       PetscSectionGetDof(candidateSection, p, &dof);
1277:       PetscSectionGetOffset(candidateSection, p, &offset);
1278:       for (d = 0; d < dof;) {
1279:         if (claims[offset+d].rank >= 0) {
1280:           const PetscInt  faceInd   = offset+d;
1281:           const PetscInt  numPoints = candidates[faceInd].index;
1282:           const PetscInt *join      = NULL;
1283:           PetscInt        joinSize, points[1024], c;

1285:           if (debug) {PetscSynchronizedPrintf(PetscObjectComm((PetscObject) dm), "[%d]    Found claim for remote point (%D, %D)\n", rank, claims[faceInd].rank, claims[faceInd].index);}
1286:           points[0] = p;
1287:           if (debug) {PetscSynchronizedPrintf(PetscObjectComm((PetscObject) dm), "[%d]      point %D\n", rank, points[0]);}
1288:           for (c = 0, ++d; c < numPoints; ++c, ++d) {
1289:             key.i = candidates[offset+d].index;
1290:             key.j = candidates[offset+d].rank;
1291:             PetscHMapIJGet(roothash, key, &root);
1292:             points[c+1] = localPoints[root];
1293:             if (debug) {PetscSynchronizedPrintf(PetscObjectComm((PetscObject) dm), "[%d]      point %D\n", rank, points[c+1]);}
1294:           }
1295:           DMPlexGetJoin(dm, numPoints+1, points, &joinSize, &join);
1296:           if (joinSize == 1) {
1297:             if (debug) {PetscSynchronizedPrintf(PetscObjectComm((PetscObject) dm), "[%d]    Found local face %D\n", rank, join[0]);}
1298:             PetscHMapISet(claimshash, join[0], faceInd);
1299:           }
1300:           DMPlexRestoreJoin(dm, numPoints+1, points, &joinSize, &join);
1301:         } else d += claims[offset+d].index+1;
1302:       }
1303:     }
1304:     if (debug) {PetscSynchronizedFlush(PetscObjectComm((PetscObject) dm), NULL);}
1305:     /* Create new pointSF from hashed claims */
1306:     PetscHMapIGetSize(claimshash, &numLocalNew);
1307:     DMPlexGetChart(dm, &pStart, &pEnd);
1308:     PetscMalloc1(numLeaves + numLocalNew, &localPointsNew);
1309:     PetscMalloc1(numLeaves + numLocalNew, &remotePointsNew);
1310:     for (p = 0; p < numLeaves; ++p) {
1311:       localPointsNew[p] = localPoints[p];
1312:       remotePointsNew[p].index = remotePoints[p].index;
1313:       remotePointsNew[p].rank  = remotePoints[p].rank;
1314:     }
1315:     p = numLeaves;
1316:     PetscHMapIGetKeys(claimshash, &p, localPointsNew);
1317:     PetscSortInt(numLocalNew, &localPointsNew[numLeaves]);
1318:     for (p = numLeaves; p < numLeaves + numLocalNew; ++p) {
1319:       PetscInt offset;
1320:       PetscHMapIGet(claimshash, localPointsNew[p], &offset);
1321:       remotePointsNew[p] = claims[offset];
1322:     }
1323:     PetscSFCreate(PetscObjectComm((PetscObject) dm), &sfPointNew);
1324:     PetscSFSetGraph(sfPointNew, pEnd-pStart, numLeaves+numLocalNew, localPointsNew, PETSC_OWN_POINTER, remotePointsNew, PETSC_OWN_POINTER);
1325:     DMSetPointSF(dm, sfPointNew);
1326:     PetscSFDestroy(&sfPointNew);
1327:     PetscHMapIDestroy(&claimshash);
1328:   }
1329:   PetscHMapIDestroy(&leafhash);
1330:   PetscHMapIJDestroy(&roothash);
1331:   PetscSectionDestroy(&candidateSection);
1332:   PetscSectionDestroy(&candidateSectionRemote);
1333:   PetscSectionDestroy(&claimSection);
1334:   PetscFree(candidates);
1335:   PetscFree(candidatesRemote);
1336:   PetscFree(claims);
1337:   PetscLogEventEnd(DMPLEX_InterpolateSF,dm,0,0,0);
1338:   return(0);
1339: }

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

1344:   Collective on dm

1346:   Input Parameters:
1347: + dm - The DMPlex object with only cells and vertices
1348: - dmInt - The interpolated DM

1350:   Output Parameter:
1351: . dmInt - The complete DMPlex object

1353:   Level: intermediate

1355:   Notes:
1356:     It does not copy over the coordinates.

1358: .seealso: DMPlexUninterpolate(), DMPlexCreateFromCellList(), DMPlexCopyCoordinates()
1359: @*/
1360: PetscErrorCode DMPlexInterpolate(DM dm, DM *dmInt)
1361: {
1362:   DM             idm, odm = dm;
1363:   PetscSF        sfPoint;
1364:   PetscInt       depth, dim, d;
1365:   const char    *name;
1366:   PetscBool      flg=PETSC_TRUE;

1372:   PetscLogEventBegin(DMPLEX_Interpolate,dm,0,0,0);
1373:   DMPlexGetDepth(dm, &depth);
1374:   DMGetDimension(dm, &dim);
1375:   if ((depth == dim) || (dim <= 1)) {
1376:     PetscObjectReference((PetscObject) dm);
1377:     idm  = dm;
1378:   } else {
1379:     for (d = 1; d < dim; ++d) {
1380:       /* Create interpolated mesh */
1381:       DMCreate(PetscObjectComm((PetscObject)dm), &idm);
1382:       DMSetType(idm, DMPLEX);
1383:       DMSetDimension(idm, dim);
1384:       if (depth > 0) {
1385:         DMPlexInterpolateFaces_Internal(odm, 1, idm);
1386:         DMGetPointSF(odm, &sfPoint);
1387:         DMPlexInterpolatePointSF(idm, sfPoint);
1388:       }
1389:       if (odm != dm) {DMDestroy(&odm);}
1390:       odm = idm;
1391:     }
1392:     PetscObjectGetName((PetscObject) dm,  &name);
1393:     PetscObjectSetName((PetscObject) idm,  name);
1394:     DMPlexCopyCoordinates(dm, idm);
1395:     DMCopyLabels(dm, idm);
1396:     PetscOptionsGetBool(((PetscObject)dm)->options, ((PetscObject)dm)->prefix, "-dm_plex_interpolate_orient_interfaces", &flg, NULL);
1397:     if (flg) {DMPlexOrientInterface(idm);}
1398:   }
1399:   {
1400:     PetscBool            isper;
1401:     const PetscReal      *maxCell, *L;
1402:     const DMBoundaryType *bd;

1404:     DMGetPeriodicity(dm,&isper,&maxCell,&L,&bd);
1405:     DMSetPeriodicity(idm,isper,maxCell,L,bd);
1406:   }
1407:   *dmInt = idm;
1408:   PetscLogEventEnd(DMPLEX_Interpolate,dm,0,0,0);
1409:   return(0);
1410: }

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

1415:   Collective on dmA

1417:   Input Parameter:
1418: . dmA - The DMPlex object with initial coordinates

1420:   Output Parameter:
1421: . dmB - The DMPlex object with copied coordinates

1423:   Level: intermediate

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

1427: .seealso: DMCopyLabels(), DMGetCoordinates(), DMGetCoordinatesLocal(), DMGetCoordinateDM(), DMGetCoordinateSection()
1428: @*/
1429: PetscErrorCode DMPlexCopyCoordinates(DM dmA, DM dmB)
1430: {
1431:   Vec            coordinatesA, coordinatesB;
1432:   VecType        vtype;
1433:   PetscSection   coordSectionA, coordSectionB;
1434:   PetscScalar   *coordsA, *coordsB;
1435:   PetscInt       spaceDim, Nf, vStartA, vStartB, vEndA, vEndB, coordSizeB, v, d;
1436:   PetscInt       cStartA, cEndA, cStartB, cEndB, cS, cE;
1437:   PetscBool      lc = PETSC_FALSE;

1443:   if (dmA == dmB) return(0);
1444:   DMPlexGetDepthStratum(dmA, 0, &vStartA, &vEndA);
1445:   DMPlexGetDepthStratum(dmB, 0, &vStartB, &vEndB);
1446:   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);
1447:   DMPlexGetHeightStratum(dmA, 0, &cStartA, &cEndA);
1448:   DMPlexGetHeightStratum(dmB, 0, &cStartB, &cEndB);
1449:   DMGetCoordinateSection(dmA, &coordSectionA);
1450:   DMGetCoordinateSection(dmB, &coordSectionB);
1451:   if (coordSectionA == coordSectionB) return(0);
1452:   PetscSectionGetNumFields(coordSectionA, &Nf);
1453:   if (!Nf) return(0);
1454:   if (Nf > 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "The number of coordinate fields must be 1, not %D", Nf);
1455:   if (!coordSectionB) {
1456:     PetscInt dim;

1458:     PetscSectionCreate(PetscObjectComm((PetscObject) coordSectionA), &coordSectionB);
1459:     DMGetCoordinateDim(dmA, &dim);
1460:     DMSetCoordinateSection(dmB, dim, coordSectionB);
1461:     PetscObjectDereference((PetscObject) coordSectionB);
1462:   }
1463:   PetscSectionSetNumFields(coordSectionB, 1);
1464:   PetscSectionGetFieldComponents(coordSectionA, 0, &spaceDim);
1465:   PetscSectionSetFieldComponents(coordSectionB, 0, spaceDim);
1466:   PetscSectionGetChart(coordSectionA, &cS, &cE);
1467:   if (cStartA <= cS && cS < cEndA) { /* localized coordinates */
1468:     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);
1469:     cS = cS - cStartA + cStartB;
1470:     cE = vEndB;
1471:     lc = PETSC_TRUE;
1472:   } else {
1473:     cS = vStartB;
1474:     cE = vEndB;
1475:   }
1476:   PetscSectionSetChart(coordSectionB, cS, cE);
1477:   for (v = vStartB; v < vEndB; ++v) {
1478:     PetscSectionSetDof(coordSectionB, v, spaceDim);
1479:     PetscSectionSetFieldDof(coordSectionB, v, 0, spaceDim);
1480:   }
1481:   if (lc) { /* localized coordinates */
1482:     PetscInt c;

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

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

1507:     PetscSectionGetOffset(coordSectionA, v + vStartA, &offA);
1508:     PetscSectionGetOffset(coordSectionB, v + vStartB, &offB);
1509:     for (d = 0; d < spaceDim; ++d) {
1510:       coordsB[offB+d] = coordsA[offA+d];
1511:     }
1512:   }
1513:   if (lc) { /* localized coordinates */
1514:     PetscInt c;

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

1519:       PetscSectionGetOffset(coordSectionA, c + cStartA, &offA);
1520:       PetscSectionGetOffset(coordSectionB, c + cStartB, &offB);
1521:       PetscSectionGetDof(coordSectionA, c + cStartA, &dof);
1522:       PetscArraycpy(coordsB + offB,coordsA + offA,dof);
1523:     }
1524:   }
1525:   VecRestoreArray(coordinatesA, &coordsA);
1526:   VecRestoreArray(coordinatesB, &coordsB);
1527:   DMSetCoordinatesLocal(dmB, coordinatesB);
1528:   VecDestroy(&coordinatesB);
1529:   return(0);
1530: }

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

1535:   Collective on dm

1537:   Input Parameter:
1538: . dm - The complete DMPlex object

1540:   Output Parameter:
1541: . dmUnint - The DMPlex object with only cells and vertices

1543:   Level: intermediate

1545:   Notes:
1546:     It does not copy over the coordinates.

1548: .seealso: DMPlexInterpolate(), DMPlexCreateFromCellList(), DMPlexCopyCoordinates()
1549: @*/
1550: PetscErrorCode DMPlexUninterpolate(DM dm, DM *dmUnint)
1551: {
1552:   DM             udm;
1553:   PetscInt       dim, vStart, vEnd, cStart, cEnd, cMax, c, maxConeSize = 0, *cone;

1559:   DMGetDimension(dm, &dim);
1560:   if (dim <= 1) {
1561:     PetscObjectReference((PetscObject) dm);
1562:     *dmUnint = dm;
1563:     return(0);
1564:   }
1565:   DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
1566:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
1567:   DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);
1568:   DMCreate(PetscObjectComm((PetscObject) dm), &udm);
1569:   DMSetType(udm, DMPLEX);
1570:   DMSetDimension(udm, dim);
1571:   DMPlexSetChart(udm, cStart, vEnd);
1572:   for (c = cStart; c < cEnd; ++c) {
1573:     PetscInt *closure = NULL, closureSize, cl, coneSize = 0;

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

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

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

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

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

1642:     DMGetPeriodicity(dm,&isper,&maxCell,&L,&bd);
1643:     DMSetPeriodicity(udm,isper,maxCell,L,bd);
1644:   }

1646:   *dmUnint = udm;
1647:   return(0);
1648: }