Actual source code: plexinterpolate.c
petsc-3.12.5 2020-03-29
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, cdim;
1437: PetscBool lc = PETSC_FALSE;
1443: if (dmA == dmB) return(0);
1444: DMGetCoordinateDim(dmA, &cdim);
1445: DMSetCoordinateDim(dmB, cdim);
1446: DMPlexGetDepthStratum(dmA, 0, &vStartA, &vEndA);
1447: DMPlexGetDepthStratum(dmB, 0, &vStartB, &vEndB);
1448: 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);
1449: DMPlexGetHeightStratum(dmA, 0, &cStartA, &cEndA);
1450: DMPlexGetHeightStratum(dmB, 0, &cStartB, &cEndB);
1451: DMGetCoordinateSection(dmA, &coordSectionA);
1452: DMGetCoordinateSection(dmB, &coordSectionB);
1453: if (coordSectionA == coordSectionB) return(0);
1454: PetscSectionGetNumFields(coordSectionA, &Nf);
1455: if (!Nf) return(0);
1456: if (Nf > 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "The number of coordinate fields must be 1, not %D", Nf);
1457: if (!coordSectionB) {
1458: PetscInt dim;
1460: PetscSectionCreate(PetscObjectComm((PetscObject) coordSectionA), &coordSectionB);
1461: DMGetCoordinateDim(dmA, &dim);
1462: DMSetCoordinateSection(dmB, dim, coordSectionB);
1463: PetscObjectDereference((PetscObject) coordSectionB);
1464: }
1465: PetscSectionSetNumFields(coordSectionB, 1);
1466: PetscSectionGetFieldComponents(coordSectionA, 0, &spaceDim);
1467: PetscSectionSetFieldComponents(coordSectionB, 0, spaceDim);
1468: PetscSectionGetChart(coordSectionA, &cS, &cE);
1469: if (cStartA <= cS && cS < cEndA) { /* localized coordinates */
1470: 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);
1471: cS = cS - cStartA + cStartB;
1472: cE = vEndB;
1473: lc = PETSC_TRUE;
1474: } else {
1475: cS = vStartB;
1476: cE = vEndB;
1477: }
1478: PetscSectionSetChart(coordSectionB, cS, cE);
1479: for (v = vStartB; v < vEndB; ++v) {
1480: PetscSectionSetDof(coordSectionB, v, spaceDim);
1481: PetscSectionSetFieldDof(coordSectionB, v, 0, spaceDim);
1482: }
1483: if (lc) { /* localized coordinates */
1484: PetscInt c;
1486: for (c = cS-cStartB; c < cEndB-cStartB; c++) {
1487: PetscInt dof;
1489: PetscSectionGetDof(coordSectionA, c + cStartA, &dof);
1490: PetscSectionSetDof(coordSectionB, c + cStartB, dof);
1491: PetscSectionSetFieldDof(coordSectionB, c + cStartB, 0, dof);
1492: }
1493: }
1494: PetscSectionSetUp(coordSectionB);
1495: PetscSectionGetStorageSize(coordSectionB, &coordSizeB);
1496: DMGetCoordinatesLocal(dmA, &coordinatesA);
1497: VecCreate(PETSC_COMM_SELF, &coordinatesB);
1498: PetscObjectSetName((PetscObject) coordinatesB, "coordinates");
1499: VecSetSizes(coordinatesB, coordSizeB, PETSC_DETERMINE);
1500: VecGetBlockSize(coordinatesA, &d);
1501: VecSetBlockSize(coordinatesB, d);
1502: VecGetType(coordinatesA, &vtype);
1503: VecSetType(coordinatesB, vtype);
1504: VecGetArray(coordinatesA, &coordsA);
1505: VecGetArray(coordinatesB, &coordsB);
1506: for (v = 0; v < vEndB-vStartB; ++v) {
1507: PetscInt offA, offB;
1509: PetscSectionGetOffset(coordSectionA, v + vStartA, &offA);
1510: PetscSectionGetOffset(coordSectionB, v + vStartB, &offB);
1511: for (d = 0; d < spaceDim; ++d) {
1512: coordsB[offB+d] = coordsA[offA+d];
1513: }
1514: }
1515: if (lc) { /* localized coordinates */
1516: PetscInt c;
1518: for (c = cS-cStartB; c < cEndB-cStartB; c++) {
1519: PetscInt dof, offA, offB;
1521: PetscSectionGetOffset(coordSectionA, c + cStartA, &offA);
1522: PetscSectionGetOffset(coordSectionB, c + cStartB, &offB);
1523: PetscSectionGetDof(coordSectionA, c + cStartA, &dof);
1524: PetscArraycpy(coordsB + offB,coordsA + offA,dof);
1525: }
1526: }
1527: VecRestoreArray(coordinatesA, &coordsA);
1528: VecRestoreArray(coordinatesB, &coordsB);
1529: DMSetCoordinatesLocal(dmB, coordinatesB);
1530: VecDestroy(&coordinatesB);
1531: return(0);
1532: }
1534: /*@
1535: DMPlexUninterpolate - Take in a mesh with all intermediate faces, edges, etc. and return a cell-vertex mesh
1537: Collective on dm
1539: Input Parameter:
1540: . dm - The complete DMPlex object
1542: Output Parameter:
1543: . dmUnint - The DMPlex object with only cells and vertices
1545: Level: intermediate
1547: Notes:
1548: It does not copy over the coordinates.
1550: .seealso: DMPlexInterpolate(), DMPlexCreateFromCellList(), DMPlexCopyCoordinates()
1551: @*/
1552: PetscErrorCode DMPlexUninterpolate(DM dm, DM *dmUnint)
1553: {
1554: DM udm;
1555: PetscInt dim, vStart, vEnd, cStart, cEnd, cMax, c, maxConeSize = 0, *cone;
1561: DMGetDimension(dm, &dim);
1562: if (dim <= 1) {
1563: PetscObjectReference((PetscObject) dm);
1564: *dmUnint = dm;
1565: return(0);
1566: }
1567: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
1568: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
1569: DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);
1570: DMCreate(PetscObjectComm((PetscObject) dm), &udm);
1571: DMSetType(udm, DMPLEX);
1572: DMSetDimension(udm, dim);
1573: DMPlexSetChart(udm, cStart, vEnd);
1574: for (c = cStart; c < cEnd; ++c) {
1575: PetscInt *closure = NULL, closureSize, cl, coneSize = 0;
1577: DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
1578: for (cl = 0; cl < closureSize*2; cl += 2) {
1579: const PetscInt p = closure[cl];
1581: if ((p >= vStart) && (p < vEnd)) ++coneSize;
1582: }
1583: DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
1584: DMPlexSetConeSize(udm, c, coneSize);
1585: maxConeSize = PetscMax(maxConeSize, coneSize);
1586: }
1587: DMSetUp(udm);
1588: PetscMalloc1(maxConeSize, &cone);
1589: for (c = cStart; c < cEnd; ++c) {
1590: PetscInt *closure = NULL, closureSize, cl, coneSize = 0;
1592: DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
1593: for (cl = 0; cl < closureSize*2; cl += 2) {
1594: const PetscInt p = closure[cl];
1596: if ((p >= vStart) && (p < vEnd)) cone[coneSize++] = p;
1597: }
1598: DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
1599: DMPlexSetCone(udm, c, cone);
1600: }
1601: PetscFree(cone);
1602: DMPlexSetHybridBounds(udm, cMax, PETSC_DETERMINE, PETSC_DETERMINE, PETSC_DETERMINE);
1603: DMPlexSymmetrize(udm);
1604: DMPlexStratify(udm);
1605: /* Reduce SF */
1606: {
1607: PetscSF sfPoint, sfPointUn;
1608: const PetscSFNode *remotePoints;
1609: const PetscInt *localPoints;
1610: PetscSFNode *remotePointsUn;
1611: PetscInt *localPointsUn;
1612: PetscInt vEnd, numRoots, numLeaves, l;
1613: PetscInt numLeavesUn = 0, n = 0;
1614: PetscErrorCode ierr;
1616: /* Get original SF information */
1617: DMGetPointSF(dm, &sfPoint);
1618: DMGetPointSF(udm, &sfPointUn);
1619: DMPlexGetDepthStratum(dm, 0, NULL, &vEnd);
1620: PetscSFGetGraph(sfPoint, &numRoots, &numLeaves, &localPoints, &remotePoints);
1621: /* Allocate space for cells and vertices */
1622: for (l = 0; l < numLeaves; ++l) if (localPoints[l] < vEnd) numLeavesUn++;
1623: /* Fill in leaves */
1624: if (vEnd >= 0) {
1625: PetscMalloc1(numLeavesUn, &remotePointsUn);
1626: PetscMalloc1(numLeavesUn, &localPointsUn);
1627: for (l = 0; l < numLeaves; l++) {
1628: if (localPoints[l] < vEnd) {
1629: localPointsUn[n] = localPoints[l];
1630: remotePointsUn[n].rank = remotePoints[l].rank;
1631: remotePointsUn[n].index = remotePoints[l].index;
1632: ++n;
1633: }
1634: }
1635: if (n != numLeavesUn) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Inconsistent number of leaves %d != %d", n, numLeavesUn);
1636: PetscSFSetGraph(sfPointUn, vEnd, numLeavesUn, localPointsUn, PETSC_OWN_POINTER, remotePointsUn, PETSC_OWN_POINTER);
1637: }
1638: }
1639: {
1640: PetscBool isper;
1641: const PetscReal *maxCell, *L;
1642: const DMBoundaryType *bd;
1644: DMGetPeriodicity(dm,&isper,&maxCell,&L,&bd);
1645: DMSetPeriodicity(udm,isper,maxCell,L,bd);
1646: }
1648: *dmUnint = udm;
1649: return(0);
1650: }