Actual source code: plexinterpolate.c
petsc-3.11.4 2019-09-28
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, numCellFacesH = 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: DMPlexGetConeSize(dm, pMax, &coneSizeH);
306: DMPlexGetCone(dm, pMax, &cone);
307: DMPlexGetRawFacesHybrid_Internal(dm, cellDim, coneSizeH, cone, &numCellFacesH, &numCellFacesT, &faceSize, &cellFaces);
308: if (faceSize < 0) {
309: PetscInt *sizes, minv, maxv;
311: /* count vertices of hybrid and non-hybrid faces */
312: PetscCalloc1(numCellFacesH, &sizes);
313: for (cf = 0; cf < numCellFacesT; ++cf) { /* These are the non-hybrid faces */
314: const PetscInt *cellFace = &cellFaces[-cf*faceSize];
315: PetscInt f;
317: for (f = 0; f < -faceSize; ++f) sizes[cf] += (cellFace[f] >= 0 ? 1 : 0);
318: }
319: PetscSortInt(numCellFacesT, sizes);
320: minv = sizes[0];
321: maxv = sizes[PetscMax(numCellFacesT-1, 0)];
322: if (minv != maxv) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_SUP, "Different number of vertices for non-hybrid face %D != %D", minv, maxv);
323: faceSizeAll = minv;
324: PetscMemzero(sizes, numCellFacesH*sizeof(PetscInt));
325: for (cf = numCellFacesT; cf < numCellFacesH; ++cf) { /* These are the hybrid faces */
326: const PetscInt *cellFace = &cellFaces[-cf*faceSize];
327: PetscInt f;
329: for (f = 0; f < -faceSize; ++f) sizes[cf-numCellFacesT] += (cellFace[f] >= 0 ? 1 : 0);
330: }
331: PetscSortInt(numCellFacesH - numCellFacesT, sizes);
332: minv = sizes[0];
333: maxv = sizes[PetscMax(numCellFacesH - numCellFacesT-1, 0)];
334: PetscFree(sizes);
335: if (minv != maxv) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_SUP, "Different number of vertices for hybrid face %D != %D", minv, maxv);
336: faceSizeAllH = minv;
337: } else { /* the size of the faces in hybrid cells is the same */
338: faceSizeAll = faceSizeAllH = faceSize;
339: }
340: DMPlexRestoreRawFacesHybrid_Internal(dm, cellDim, coneSizeH, cone, &numCellFacesH, &numCellFacesT, &faceSize, &cellFaces);
341: } else if (pEnd[cellDepth] > pStart[cellDepth]) {
342: DMPlexGetFaces_Internal(dm, cellDim, pStart[cellDepth], NULL, &faceSizeAll, NULL);
343: }
344: if (faceSizeAll > 4) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Do not support interpolation of meshes with faces of %D vertices", faceSizeAll);
346: /* With hybrid grids, we first iterate on hybrid cells and start numbering the non-hybrid faces
347: Then, faces for non-hybrid cells are numbered.
348: This is to guarantee consistent orientations (all 0) of all the points in the cone of the hybrid cells */
349: PetscHashIJKLCreate(&faceTable);
350: for (outerloop = 0, face = pStart[faceDepth]; outerloop < 2; outerloop++) {
351: PetscInt start, end;
353: start = outerloop == 0 ? pMax : pStart[cellDepth];
354: end = outerloop == 0 ? pEnd[cellDepth] : pMax;
355: for (c = start; c < end; ++c) {
356: const PetscInt *cellFaces;
357: PetscInt numCellFaces, faceSize, faceSizeInc, cf;
359: if (c < pMax) {
360: DMPlexGetFaces_Internal(dm, cellDim, c, &numCellFaces, &faceSize, &cellFaces);
361: if (faceSize != faceSizeAll) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Inconsistent face for cell %D of size %D != %D", c, faceSize, faceSizeAll);
362: } else { /* Hybrid cell */
363: const PetscInt *cone;
364: PetscInt numCellFacesN, coneSize;
366: DMPlexGetConeSize(dm, c, &coneSize);
367: if (coneSize != coneSizeH) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unexpected hybrid coneSize %D != %D", coneSize, coneSizeH);
368: DMPlexGetCone(dm, c, &cone);
369: DMPlexGetRawFacesHybrid_Internal(dm, cellDim, coneSize, cone, &numCellFaces, &numCellFacesN, &faceSize, &cellFaces);
370: if (numCellFaces != numCellFacesH) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unexpected numCellFaces %D != %D for hybrid cell %D", numCellFaces, numCellFacesH, c);
371: faceSize = PetscMax(faceSize, -faceSize);
372: if (faceSize > 4) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Do not support interpolation of meshes with faces of %D vertices", faceSize);
373: numCellFaces = numCellFacesN; /* process only non-hybrid faces */
374: }
375: faceSizeInc = faceSize;
376: for (cf = 0; cf < numCellFaces; ++cf) {
377: const PetscInt *cellFace = &cellFaces[cf*faceSizeInc];
378: PetscInt faceSizeH = faceSize;
379: PetscHashIJKLKey key;
380: PetscHashIter iter;
381: PetscBool missing;
383: if (faceSizeInc == 2) {
384: key.i = PetscMin(cellFace[0], cellFace[1]);
385: key.j = PetscMax(cellFace[0], cellFace[1]);
386: key.k = PETSC_MAX_INT;
387: key.l = PETSC_MAX_INT;
388: } else {
389: key.i = cellFace[0];
390: key.j = cellFace[1];
391: key.k = cellFace[2];
392: key.l = faceSize > 3 ? (cellFace[3] < 0 ? faceSizeH = 3, PETSC_MAX_INT : cellFace[3]) : PETSC_MAX_INT;
393: PetscSortInt(faceSize, (PetscInt *) &key);
394: }
395: /* this check is redundant for non-hybrid meshes */
396: if (faceSizeH != faceSizeAll) SETERRQ4(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unexpected number of vertices for face %D of point %D -> %D != %D", cf, c, faceSizeH, faceSizeAll);
397: PetscHashIJKLPut(faceTable, key, &iter, &missing);
398: if (missing) {PetscHashIJKLIterSet(faceTable, iter, face++);}
399: }
400: if (c < pMax) {
401: DMPlexRestoreFaces_Internal(dm, cellDim, c, &numCellFaces, &faceSize, &cellFaces);
402: } else {
403: DMPlexRestoreRawFacesHybrid_Internal(dm, cellDim, coneSizeH, NULL, NULL, NULL, NULL, &cellFaces);
404: }
405: }
406: }
407: pEnd[faceDepth] = face;
409: /* Second pass for hybrid meshes: number hybrid faces */
410: for (c = pMax; c < pEnd[cellDepth]; ++c) {
411: const PetscInt *cellFaces, *cone;
412: PetscInt numCellFaces, numCellFacesN, faceSize, cf, coneSize;
414: DMPlexGetConeSize(dm, c, &coneSize);
415: DMPlexGetCone(dm, c, &cone);
416: DMPlexGetRawFacesHybrid_Internal(dm, cellDim, coneSize, cone, &numCellFaces, &numCellFacesN, &faceSize, &cellFaces);
417: if (numCellFaces != numCellFacesH) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unexpected hybrid numCellFaces %D != %D", numCellFaces, numCellFacesH);
418: faceSize = PetscMax(faceSize, -faceSize);
419: for (cf = numCellFacesN; cf < numCellFaces; ++cf) { /* These are the hybrid faces */
420: const PetscInt *cellFace = &cellFaces[cf*faceSize];
421: PetscHashIJKLKey key;
422: PetscHashIter iter;
423: PetscBool missing;
424: PetscInt faceSizeH = faceSize;
426: if (faceSize == 2) {
427: key.i = PetscMin(cellFace[0], cellFace[1]);
428: key.j = PetscMax(cellFace[0], cellFace[1]);
429: key.k = PETSC_MAX_INT;
430: key.l = PETSC_MAX_INT;
431: } else {
432: key.i = cellFace[0];
433: key.j = cellFace[1];
434: key.k = cellFace[2];
435: key.l = faceSize > 3 ? (cellFace[3] < 0 ? faceSizeH = 3, PETSC_MAX_INT : cellFace[3]) : PETSC_MAX_INT;
436: PetscSortInt(faceSize, (PetscInt *) &key);
437: }
438: 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);
439: PetscHashIJKLPut(faceTable, key, &iter, &missing);
440: if (missing) {PetscHashIJKLIterSet(faceTable, iter, face++);}
441: }
442: DMPlexRestoreRawFacesHybrid_Internal(dm, cellDim, coneSize, cone, &numCellFaces, &numCellFacesN, &faceSize, &cellFaces);
443: }
444: faceH = face - pEnd[faceDepth];
445: if (faceH) {
446: if (fMax == PETSC_DETERMINE) fMax = pEnd[faceDepth];
447: else if (eMax == PETSC_DETERMINE) eMax = pEnd[faceDepth];
448: else SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of unassigned hybrid facets %D for cellDim %D and dimension %D", faceH, cellDim, dim);
449: }
450: pEnd[faceDepth] = face;
451: PetscHashIJKLDestroy(&faceTable);
452: /* Count new points */
453: for (d = 0; d <= depth; ++d) {
454: numPoints += pEnd[d]-pStart[d];
455: }
456: DMPlexSetChart(idm, 0, numPoints);
457: /* Set cone sizes */
458: for (d = 0; d <= depth; ++d) {
459: PetscInt coneSize, p;
461: if (d == faceDepth) {
462: /* I see no way to do this if we admit faces of different shapes */
463: for (p = pStart[d]; p < pEnd[d]-faceH; ++p) {
464: DMPlexSetConeSize(idm, p, faceSizeAll);
465: }
466: for (p = pEnd[d]-faceH; p < pEnd[d]; ++p) {
467: DMPlexSetConeSize(idm, p, faceSizeAllH);
468: }
469: } else if (d == cellDepth) {
470: for (p = pStart[d]; p < pEnd[d]; ++p) {
471: /* Number of cell faces may be different from number of cell vertices*/
472: if (p < pMax) {
473: DMPlexGetFaces_Internal(dm, cellDim, p, &coneSize, NULL, NULL);
474: } else {
475: DMPlexGetFacesHybrid_Internal(dm, cellDim, p, &coneSize, NULL, NULL, NULL);
476: }
477: DMPlexSetConeSize(idm, p, coneSize);
478: }
479: } else {
480: for (p = pStart[d]; p < pEnd[d]; ++p) {
481: DMPlexGetConeSize(dm, p, &coneSize);
482: DMPlexSetConeSize(idm, p, coneSize);
483: }
484: }
485: }
486: DMSetUp(idm);
487: /* Get face cones from subsets of cell vertices */
488: if (faceSizeAll > 4) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Do not support interpolation of meshes with faces of %D vertices", faceSizeAll);
489: PetscHashIJKLCreate(&faceTable);
490: for (d = depth; d > cellDepth; --d) {
491: const PetscInt *cone;
492: PetscInt p;
494: for (p = pStart[d]; p < pEnd[d]; ++p) {
495: DMPlexGetCone(dm, p, &cone);
496: DMPlexSetCone(idm, p, cone);
497: DMPlexGetConeOrientation(dm, p, &cone);
498: DMPlexSetConeOrientation(idm, p, cone);
499: }
500: }
501: for (outerloop = 0, face = pStart[faceDepth]; outerloop < 2; outerloop++) {
502: PetscInt start, end;
504: start = outerloop == 0 ? pMax : pStart[cellDepth];
505: end = outerloop == 0 ? pEnd[cellDepth] : pMax;
506: for (c = start; c < end; ++c) {
507: const PetscInt *cellFaces;
508: PetscInt numCellFaces, faceSize, faceSizeInc, cf;
510: if (c < pMax) {
511: DMPlexGetFaces_Internal(dm, cellDim, c, &numCellFaces, &faceSize, &cellFaces);
512: if (faceSize != faceSizeAll) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Inconsistent face for cell %D of size %D != %D", c, faceSize, faceSizeAll);
513: } else {
514: const PetscInt *cone;
515: PetscInt numCellFacesN, coneSize;
517: DMPlexGetConeSize(dm, c, &coneSize);
518: if (coneSize != coneSizeH) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unexpected hybrid coneSize %D != %D", coneSize, coneSizeH);
519: DMPlexGetCone(dm, c, &cone);
520: DMPlexGetRawFacesHybrid_Internal(dm, cellDim, coneSize, cone, &numCellFaces, &numCellFacesN, &faceSize, &cellFaces);
521: if (numCellFaces != numCellFacesH) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unexpected numCellFaces %D != %D for hybrid cell %D", numCellFaces, numCellFacesH, c);
522: faceSize = PetscMax(faceSize, -faceSize);
523: if (faceSize > 4) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Do not support interpolation of meshes with faces of %D vertices", faceSize);
524: numCellFaces = numCellFacesN; /* process only non-hybrid faces */
525: }
526: faceSizeInc = faceSize;
527: for (cf = 0; cf < numCellFaces; ++cf) {
528: const PetscInt *cellFace = &cellFaces[cf*faceSizeInc];
529: PetscHashIJKLKey key;
530: PetscHashIter iter;
531: PetscBool missing;
533: if (faceSizeInc == 2) {
534: key.i = PetscMin(cellFace[0], cellFace[1]);
535: key.j = PetscMax(cellFace[0], cellFace[1]);
536: key.k = PETSC_MAX_INT;
537: key.l = PETSC_MAX_INT;
538: } else {
539: key.i = cellFace[0];
540: key.j = cellFace[1];
541: key.k = cellFace[2];
542: key.l = faceSizeInc > 3 ? (cellFace[3] < 0 ? faceSize = 3, PETSC_MAX_INT : cellFace[3]) : PETSC_MAX_INT;
543: PetscSortInt(faceSizeInc, (PetscInt *) &key);
544: }
545: PetscHashIJKLPut(faceTable, key, &iter, &missing);
546: if (missing) {
547: DMPlexSetCone(idm, face, cellFace);
548: PetscHashIJKLIterSet(faceTable, iter, face);
549: DMPlexInsertCone(idm, c, cf, face++);
550: } else {
551: const PetscInt *cone;
552: PetscInt coneSize, ornt, i, j, f;
554: PetscHashIJKLIterGet(faceTable, iter, &f);
555: DMPlexInsertCone(idm, c, cf, f);
556: /* Orient face: Do not allow reverse orientation at the first vertex */
557: DMPlexGetConeSize(idm, f, &coneSize);
558: DMPlexGetCone(idm, f, &cone);
559: 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);
560: /* - First find the initial vertex */
561: for (i = 0; i < faceSize; ++i) if (cellFace[0] == cone[i]) break;
562: /* - Try forward comparison */
563: for (j = 0; j < faceSize; ++j) if (cellFace[j] != cone[(i+j)%faceSize]) break;
564: if (j == faceSize) {
565: if ((faceSize == 2) && (i == 1)) ornt = -2;
566: else ornt = i;
567: } else {
568: /* - Try backward comparison */
569: for (j = 0; j < faceSize; ++j) if (cellFace[j] != cone[(i+faceSize-j)%faceSize]) break;
570: if (j == faceSize) {
571: if (i == 0) ornt = -faceSize;
572: else ornt = -i;
573: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Could not determine face orientation");
574: }
575: DMPlexInsertConeOrientation(idm, c, cf, ornt);
576: }
577: }
578: if (c < pMax) {
579: DMPlexRestoreFaces_Internal(dm, cellDim, c, &numCellFaces, &faceSize, &cellFaces);
580: } else {
581: DMPlexRestoreRawFacesHybrid_Internal(dm, cellDim, coneSizeH, NULL, NULL, NULL, NULL, &cellFaces);
582: }
583: }
584: }
585: /* Second pass for hybrid meshes: orient hybrid faces */
586: for (c = pMax; c < pEnd[cellDepth]; ++c) {
587: const PetscInt *cellFaces, *cone;
588: PetscInt numCellFaces, numCellFacesN, faceSize, cf, coneSize;
590: DMPlexGetConeSize(dm, c, &coneSize);
591: DMPlexGetCone(dm, c, &cone);
592: DMPlexGetRawFacesHybrid_Internal(dm, cellDim, coneSize, cone, &numCellFaces, &numCellFacesN, &faceSize, &cellFaces);
593: if (numCellFaces != numCellFacesH) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unexpected hybrid numCellFaces %D != %D", numCellFaces, numCellFacesH);
594: faceSize = PetscMax(faceSize, -faceSize);
595: for (cf = numCellFacesN; cf < numCellFaces; ++cf) { /* These are the hybrid faces */
596: const PetscInt *cellFace = &cellFaces[cf*faceSize];
597: PetscHashIJKLKey key;
598: PetscHashIter iter;
599: PetscBool missing;
600: PetscInt faceSizeH = faceSize;
602: if (faceSize == 2) {
603: key.i = PetscMin(cellFace[0], cellFace[1]);
604: key.j = PetscMax(cellFace[0], cellFace[1]);
605: key.k = PETSC_MAX_INT;
606: key.l = PETSC_MAX_INT;
607: } else {
608: key.i = cellFace[0];
609: key.j = cellFace[1];
610: key.k = cellFace[2];
611: key.l = faceSize > 3 ? (cellFace[3] < 0 ? faceSizeH = 3, PETSC_MAX_INT : cellFace[3]) : PETSC_MAX_INT;
612: PetscSortInt(faceSize, (PetscInt *) &key);
613: }
614: 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);
615: PetscHashIJKLPut(faceTable, key, &iter, &missing);
616: if (missing) {
617: DMPlexSetCone(idm, face, cellFace);
618: PetscHashIJKLIterSet(faceTable, iter, face);
619: DMPlexInsertCone(idm, c, cf, face++);
620: } else {
621: const PetscInt *cone;
622: PetscInt coneSize, ornt, i, j, f;
624: PetscHashIJKLIterGet(faceTable, iter, &f);
625: DMPlexInsertCone(idm, c, cf, f);
626: /* Orient face: Do not allow reverse orientation at the first vertex */
627: DMPlexGetConeSize(idm, f, &coneSize);
628: DMPlexGetCone(idm, f, &cone);
629: 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);
630: /* - First find the initial vertex */
631: for (i = 0; i < faceSizeH; ++i) if (cellFace[0] == cone[i]) break;
632: /* - Try forward comparison */
633: for (j = 0; j < faceSizeH; ++j) if (cellFace[j] != cone[(i+j)%faceSizeH]) break;
634: if (j == faceSizeH) {
635: if ((faceSizeH == 2) && (i == 1)) ornt = -2;
636: else ornt = i;
637: } else {
638: /* - Try backward comparison */
639: for (j = 0; j < faceSizeH; ++j) if (cellFace[j] != cone[(i+faceSizeH-j)%faceSizeH]) break;
640: if (j == faceSizeH) {
641: if (i == 0) ornt = -faceSizeH;
642: else ornt = -i;
643: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Could not determine face orientation");
644: }
645: DMPlexInsertConeOrientation(idm, c, cf, ornt);
646: }
647: }
648: DMPlexRestoreRawFacesHybrid_Internal(dm, cellDim, coneSize, cone, &numCellFaces, &numCellFacesN, &faceSize, &cellFaces);
649: }
650: 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]);
651: PetscFree2(pStart,pEnd);
652: PetscHashIJKLDestroy(&faceTable);
653: PetscFree2(pStart,pEnd);
654: DMPlexSetHybridBounds(idm, cMax, fMax, eMax, vMax);
655: DMPlexSymmetrize(idm);
656: DMPlexStratify(idm);
657: return(0);
658: }
660: PetscErrorCode DMPlexOrientCell(DM dm, PetscInt p, PetscInt masterConeSize, const PetscInt masterCone[])
661: {
662: PetscInt coneSize;
663: PetscInt start1=0;
664: PetscBool reverse1=PETSC_FALSE;
665: const PetscInt *cone=NULL;
669: DMPlexGetConeSize(dm, p, &coneSize);
670: if (!coneSize) return(0); /* do nothing for points with no cone */
671: DMPlexGetCone(dm, p, &cone);
672: DMPlexFixFaceOrientations_Orient_Private(coneSize, masterConeSize, masterCone, cone, &start1, &reverse1);
673: #if defined(PETSC_USE_DEBUG)
674: 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]);
675: #endif
676: DMPlexOrientCell_Internal(dm, p, start1, reverse1);
677: #if defined(PETSC_USE_DEBUG)
678: {
679: PetscInt c;
680: DMPlexGetCone(dm, p, &cone);
681: for (c = 0; c < 2; c++) {
682: 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);
683: }
684: }
685: #endif
686: return(0);
687: }
689: PetscErrorCode DMPlexOrientCell_Internal(DM dm, PetscInt p, PetscInt start1, PetscBool reverse1)
690: {
691: PetscInt i, j, k, maxConeSize, coneSize, coneConeSize, supportSize, supportConeSize;
692: PetscInt start0, start;
693: PetscBool reverse0, reverse;
694: PetscInt newornt;
695: const PetscInt *cone=NULL, *support=NULL, *supportCone=NULL, *ornts=NULL;
696: PetscInt *newcone=NULL, *newornts=NULL;
700: if (!start1 && !reverse1) return(0);
701: DMPlexGetConeSize(dm, p, &coneSize);
702: if (!coneSize) return(0); /* do nothing for points with no cone */
703: DMPlexGetCone(dm, p, &cone);
704: DMPlexGetMaxSizes(dm, &maxConeSize, NULL);
705: /* permute p's cone and orientations */
706: DMPlexGetConeOrientation(dm, p, &ornts);
707: DMGetWorkArray(dm, maxConeSize, MPIU_INT, &newcone);
708: DMGetWorkArray(dm, maxConeSize, MPIU_INT, &newornts);
709: DMPlexFixFaceOrientations_Permute_Private(coneSize, cone, start1, reverse1, newcone);
710: DMPlexFixFaceOrientations_Permute_Private(coneSize, ornts, start1, reverse1, newornts);
711: /* if direction of p (face) is flipped, flip also p's cone points (edges) */
712: if (reverse1) {
713: for (i=0; i<coneSize; i++) {
714: DMPlexGetConeSize(dm, cone[i], &coneConeSize);
715: DMPlexFixFaceOrientations_Translate_Private(newornts[i], &start0, &reverse0);
716: DMPlexFixFaceOrientations_Combine_Private(coneConeSize, start0, reverse0, 1, PETSC_FALSE, &start, &reverse);
717: DMPlexFixFaceOrientations_TranslateBack_Private(coneConeSize, start, reverse, &newornts[i]);
718: }
719: }
720: DMPlexSetConeOrientation(dm, p, newornts);
721: /* fix oriention of p within cones of p's support points */
722: DMPlexGetSupport(dm, p, &support);
723: DMPlexGetSupportSize(dm, p, &supportSize);
724: for (j=0; j<supportSize; j++) {
725: DMPlexGetCone(dm, support[j], &supportCone);
726: DMPlexGetConeSize(dm, support[j], &supportConeSize);
727: DMPlexGetConeOrientation(dm, support[j], &ornts);
728: for (k=0; k<supportConeSize; k++) {
729: if (supportCone[k] != p) continue;
730: DMPlexFixFaceOrientations_Translate_Private(ornts[k], &start0, &reverse0);
731: DMPlexFixFaceOrientations_Combine_Private(coneSize, start0, reverse0, start1, reverse1, &start, &reverse);
732: DMPlexFixFaceOrientations_TranslateBack_Private(coneSize, start, reverse, &newornt);
733: DMPlexInsertConeOrientation(dm, support[j], k, newornt);
734: }
735: }
736: /* rewrite cone */
737: DMPlexSetCone(dm, p, newcone);
738: DMRestoreWorkArray(dm, maxConeSize, MPIU_INT, &newcone);
739: DMRestoreWorkArray(dm, maxConeSize, MPIU_INT, &newornts);
740: return(0);
741: }
743: static PetscErrorCode SortRmineRremoteByRemote_Private(PetscSF sf, PetscInt *rmine1[], PetscInt *rremote1[])
744: {
745: PetscInt nleaves;
746: PetscInt nranks;
747: const PetscMPIInt *ranks=NULL;
748: const PetscInt *roffset=NULL, *rmine=NULL, *rremote=NULL;
749: PetscInt n, o, r;
750: PetscErrorCode ierr;
753: PetscSFGetRanks(sf, &nranks, &ranks, &roffset, &rmine, &rremote);
754: nleaves = roffset[nranks];
755: PetscMalloc2(nleaves, rmine1, nleaves, rremote1);
756: for (r=0; r<nranks; r++) {
757: /* simultaneously sort rank-wise portions of rmine & rremote by values in rremote
758: - to unify order with the other side */
759: o = roffset[r];
760: n = roffset[r+1] - o;
761: PetscMemcpy(&(*rmine1)[o], &rmine[o], n*sizeof(PetscInt));
762: PetscMemcpy(&(*rremote1)[o], &rremote[o], n*sizeof(PetscInt));
763: PetscSortIntWithArray(n, &(*rremote1)[o], &(*rmine1)[o]);
764: }
765: return(0);
766: }
768: PetscErrorCode DMPlexOrientInterface(DM dm)
769: {
770: PetscSF sf=NULL;
771: PetscInt (*roots)[2], (*leaves)[2];
772: PetscMPIInt (*rootsRanks)[2], (*leavesRanks)[2];
773: const PetscInt *locals=NULL;
774: const PetscSFNode *remotes=NULL;
775: PetscInt nroots, nleaves, p, c;
776: PetscInt nranks, n, o, r;
777: const PetscMPIInt *ranks=NULL;
778: const PetscInt *roffset=NULL;
779: PetscInt *rmine1=NULL, *rremote1=NULL; /* rmine and rremote copies simultaneously sorted by rank and rremote */
780: const PetscInt *cone=NULL;
781: PetscInt coneSize, ind0;
782: MPI_Comm comm;
783: PetscMPIInt rank;
784: PetscInt debug = 0;
785: PetscErrorCode ierr;
788: DMGetPointSF(dm, &sf);
789: PetscSFGetGraph(sf, &nroots, &nleaves, &locals, &remotes);
790: if (nroots < 0) return(0);
791: PetscSFSetUp(sf);
792: PetscSFGetRanks(sf, &nranks, &ranks, &roffset, NULL, NULL);
793: #if defined(PETSC_USE_DEBUG)
794: DMViewFromOptions(dm, NULL, "-before_fix_dm_view");
795: DMPlexCheckPointSF(dm);
796: #endif
797: SortRmineRremoteByRemote_Private(sf, &rmine1, &rremote1);
798: PetscMalloc4(nroots, &roots, nroots, &leaves, nroots, &rootsRanks, nroots, &leavesRanks);
799: PetscObjectGetComm((PetscObject) dm, &comm);
800: MPI_Comm_rank(comm, &rank);
801: if (debug && rank == 0) {PetscSynchronizedPrintf(comm, "Roots\n");}
802: for (p = 0; p < nroots; ++p) {
803: DMPlexGetConeSize(dm, p, &coneSize);
804: DMPlexGetCone(dm, p, &cone);
805: /* Translate all points to root numbering */
806: for (c = 0; c < 2; c++) {
807: if (coneSize > 1) {
808: PetscFindInt(cone[c], nleaves, locals, &ind0);
809: if (ind0 < 0) {
810: roots[p][c] = cone[c];
811: rootsRanks[p][c] = rank;
812: } else {
813: roots[p][c] = remotes[ind0].index;
814: rootsRanks[p][c] = remotes[ind0].rank;
815: }
816: } else {
817: roots[p][c] = -1;
818: rootsRanks[p][c] = -1;
819: }
820: }
821: }
822: if (debug) {
823: for (p = 0; p < nroots; ++p) {
824: DMPlexGetConeSize(dm, p, &coneSize);
825: DMPlexGetCone(dm, p, &cone);
826: if (coneSize > 1) {
827: 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]);
828: }
829: }
830: PetscSynchronizedFlush(comm, NULL);
831: }
832: for (p = 0; p < nroots; ++p) {
833: for (c = 0; c < 2; c++) {
834: leaves[p][c] = -2;
835: leavesRanks[p][c] = -2;
836: }
837: }
838: PetscSFBcastBegin(sf, MPIU_2INT, roots, leaves);
839: PetscSFBcastBegin(sf, MPI_2INT, rootsRanks, leavesRanks);
840: PetscSFBcastEnd(sf, MPIU_2INT, roots, leaves);
841: PetscSFBcastEnd(sf, MPI_2INT, rootsRanks, leavesRanks);
842: if (debug) {PetscSynchronizedFlush(comm, NULL);}
843: if (debug && rank == 0) {PetscSynchronizedPrintf(comm, "Referred leaves\n");}
844: for (p = 0; p < nroots; ++p) {
845: if (leaves[p][0] < 0) continue;
846: DMPlexGetConeSize(dm, p, &coneSize);
847: DMPlexGetCone(dm, p, &cone);
848: 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]);}
849: if ((leaves[p][0] != roots[p][0]) || (leaves[p][1] != roots[p][1]) || (leavesRanks[p][0] != rootsRanks[p][0]) || (leavesRanks[p][0] != rootsRanks[p][0])) {
850: PetscInt masterCone[2];
851: /* Translate these two cone points back to leave numbering */
852: for (c = 0; c < 2; c++) {
853: 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]);
854: /* Find index of rank leavesRanks[p][c] among remote ranks */
855: /* No need for PetscMPIIntCast because these integers were originally cast from PetscMPIInt. */
856: PetscFindMPIInt((PetscMPIInt)leavesRanks[p][c], nranks, ranks, &r);
857: 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]);
858: /* Find point leaves[p][c] among remote points aimed at rank leavesRanks[p][c] */
859: o = roffset[r];
860: n = roffset[r+1] - o;
861: PetscFindInt(leaves[p][c], n, &rremote1[o], &ind0);
862: 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]);
863: /* Get the corresponding local point */
864: masterCone[c] = rmine1[o+ind0];
865: }
866: if (debug) {PetscSynchronizedPrintf(comm, "[%d] %D: masterCone=[%D %D]\n", rank, p, masterCone[0], masterCone[1]);}
867: /* Vaclav's note: Here we only compare first 2 points of the cone. Full cone size would lead to stronger self-checking. */
868: DMPlexOrientCell(dm, p, 2, masterCone);
869: }
870: }
871: #if defined(PETSC_USE_DEBUG)
872: DMViewFromOptions(dm, NULL, "-after_fix_dm_view");
873: for (r = 0; r < nleaves; ++r) {
874: p = locals[r];
875: DMPlexGetConeSize(dm, p, &coneSize);
876: if (!coneSize) continue;
877: DMPlexGetCone(dm, p, &cone);
878: for (c = 0; c < 2; c++) {
879: PetscFindInt(cone[c], nleaves, locals, &ind0);
880: 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]);
881: if (leaves[p][c] != remotes[ind0].index || leavesRanks[p][c] != remotes[ind0].rank) {
882: if (leavesRanks[p][c] == rank) {
883: PetscInt ind1;
884: PetscFindInt(leaves[p][c], nleaves, locals, &ind1);
885: if (ind1 < 0) {
886: 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]);
887: } else {
888: 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);
889: }
890: } else {
891: 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]);
892: }
893: }
894: }
895: }
896: #endif
897: if (debug) {PetscSynchronizedFlush(comm, NULL);}
898: PetscFree4(roots, leaves, rootsRanks, leavesRanks);
899: PetscFree2(rmine1, rremote1);
900: return(0);
901: }
903: static PetscErrorCode IntArrayViewFromOptions(MPI_Comm comm, const char opt[], const char name[], const char idxname[], const char valname[], PetscInt n, const PetscInt a[])
904: {
905: PetscInt idx;
906: PetscMPIInt rank;
907: PetscBool flg;
911: PetscOptionsHasName(NULL, NULL, opt, &flg);
912: if (!flg) return(0);
913: MPI_Comm_rank(comm, &rank);
914: PetscSynchronizedPrintf(comm, "[%d]%s:\n", rank, name);
915: for (idx = 0; idx < n; ++idx) {PetscSynchronizedPrintf(comm, "[%d]%s %D %s %D\n", rank, idxname, idx, valname, a[idx]);}
916: PetscSynchronizedFlush(comm, NULL);
917: return(0);
918: }
920: static PetscErrorCode SFNodeArrayViewFromOptions(MPI_Comm comm, const char opt[], const char name[], const char idxname[], PetscInt n, const PetscSFNode a[])
921: {
922: PetscInt idx;
923: PetscMPIInt rank;
924: PetscBool flg;
928: PetscOptionsHasName(NULL, NULL, opt, &flg);
929: if (!flg) return(0);
930: MPI_Comm_rank(comm, &rank);
931: PetscSynchronizedPrintf(comm, "[%d]%s:\n", rank, name);
932: if (idxname) {
933: 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);}
934: } else {
935: for (idx = 0; idx < n; ++idx) {PetscSynchronizedPrintf(comm, "[%d]rank %D index %D\n", rank, a[idx].rank, a[idx].index);}
936: }
937: PetscSynchronizedFlush(comm, NULL);
938: return(0);
939: }
941: static PetscErrorCode DMPlexMapToLocalPoint(PetscHMapIJ roothash, const PetscInt localPoints[], PetscMPIInt rank, PetscSFNode remotePoint, PetscInt *localPoint)
942: {
946: if (remotePoint.rank == rank) {
947: *localPoint = remotePoint.index;
948: } else {
949: PetscHashIJKey key;
950: PetscInt root;
952: key.i = remotePoint.index;
953: key.j = remotePoint.rank;
954: PetscHMapIJGet(roothash, key, &root);
955: if (root >= 0) {
956: *localPoint = localPoints[root];
957: } else PetscFunctionReturn(1);
958: }
959: return(0);
960: }
962: /*@
963: DMPlexInterpolatePointSF - Insert interpolated points in the overlap into the PointSF in parallel, following local interpolation
965: Collective on DM
967: Input Parameters:
968: + dm - The interpolated DM
969: - pointSF - The initial SF without interpolated points
971: Output Parameter:
972: . pointSF - The SF including interpolated points
974: Level: intermediate
976: 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
978: .keywords: mesh
979: .seealso: DMPlexInterpolate(), DMPlexUninterpolate()
980: @*/
981: PetscErrorCode DMPlexInterpolatePointSF(DM dm, PetscSF pointSF)
982: {
983: /*
984: Okay, the algorithm is:
985: - Take each point in the overlap (root)
986: - Look at the neighboring points in the overlap (candidates)
987: - Send these candidate points to neighbors
988: - Neighbor checks for edge between root and candidate
989: - If edge is found, it replaces candidate point with edge point
990: - Send back the overwritten candidates (claims)
991: - Original guy checks for edges, different from original candidate, and gets its own edge
992: - This pair is put into SF
994: We need a new algorithm that tolerates groups larger than 2.
995: - Take each point in the overlap (root)
996: - Find all collections of points in the overlap which make faces (do early join)
997: - Send collections as candidates (add size as first number)
998: - Make sure to send collection to all owners of all overlap points in collection
999: - Neighbor check for face in collections
1000: - If face is found, it replaces candidate point with face point
1001: - Send back the overwritten candidates (claims)
1002: - Original guy checks for faces, different from original candidate, and gets its own face
1003: - This pair is put into SF
1004: */
1005: PetscHMapI leafhash;
1006: PetscHMapIJ roothash;
1007: const PetscInt *localPoints, *rootdegree;
1008: const PetscSFNode *remotePoints;
1009: PetscSFNode *candidates, *candidatesRemote, *claims;
1010: PetscSection candidateSection, candidateSectionRemote, claimSection;
1011: PetscInt numLeaves, l, numRoots, r, candidatesSize, candidatesRemoteSize;
1012: PetscMPIInt size, rank;
1013: PetscHashIJKey key;
1014: PetscBool debug = PETSC_FALSE;
1015: PetscErrorCode ierr;
1018: PetscOptionsHasName(NULL, ((PetscObject) dm)->prefix, "-dmplex_interp_debug", &debug);
1019: MPI_Comm_size(PetscObjectComm((PetscObject) dm), &size);
1020: MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank);
1021: PetscSFGetGraph(pointSF, &numRoots, &numLeaves, &localPoints, &remotePoints);
1022: if (size < 2 || numRoots < 0) return(0);
1023: PetscObjectViewFromOptions((PetscObject) dm, NULL, "-dm_interp_pre_view");
1024: PetscObjectViewFromOptions((PetscObject) pointSF, NULL, "-petscsf_interp_pre_view");
1025: PetscLogEventBegin(DMPLEX_InterpolateSF,dm,0,0,0);
1026: /* Build hashes of points in the SF for efficient lookup */
1027: PetscHMapICreate(&leafhash);
1028: PetscHMapIJCreate(&roothash);
1029: for (l = 0; l < numLeaves; ++l) {
1030: key.i = remotePoints[l].index;
1031: key.j = remotePoints[l].rank;
1032: PetscHMapISet(leafhash, localPoints[l], l);
1033: PetscHMapIJSet(roothash, key, l);
1034: }
1035: /* Compute root degree to identify shared points */
1036: PetscSFComputeDegreeBegin(pointSF, &rootdegree);
1037: PetscSFComputeDegreeEnd(pointSF, &rootdegree);
1038: IntArrayViewFromOptions(PetscObjectComm((PetscObject) dm), "-interp_root_degree_view", "Root degree", "point", "degree", numRoots, rootdegree);
1039: /* Build a section / SFNode array of candidate points (face bd points) in the cone(support(leaf)),
1040: where each candidate is defined by a set of remote points (roots) for the other points that define the face. */
1041: PetscSectionCreate(PetscObjectComm((PetscObject) dm), &candidateSection);
1042: PetscSectionSetChart(candidateSection, 0, numRoots);
1043: {
1044: PetscHMapIJ facehash;
1046: PetscHMapIJCreate(&facehash);
1047: for (l = 0; l < numLeaves; ++l) {
1048: const PetscInt localPoint = localPoints[l];
1049: const PetscInt *support;
1050: PetscInt supportSize, s;
1052: if (debug) {PetscSynchronizedPrintf(PetscObjectComm((PetscObject) dm), "[%d] Checking local point %D\n", rank, localPoint);}
1053: DMPlexGetSupportSize(dm, localPoint, &supportSize);
1054: DMPlexGetSupport(dm, localPoint, &support);
1055: for (s = 0; s < supportSize; ++s) {
1056: const PetscInt face = support[s];
1057: const PetscInt *cone;
1058: PetscInt coneSize, c, f, root;
1059: PetscBool isFace = PETSC_TRUE;
1061: /* Only add face once */
1062: if (debug) {PetscSynchronizedPrintf(PetscObjectComm((PetscObject) dm), "[%d] Support point %D\n", rank, face);}
1063: key.i = localPoint;
1064: key.j = face;
1065: PetscHMapIJGet(facehash, key, &f);
1066: if (f >= 0) continue;
1067: DMPlexGetConeSize(dm, face, &coneSize);
1068: DMPlexGetCone(dm, face, &cone);
1069: /* If a cone point does not map to leaves on any proc, then do not put face in SF */
1070: for (c = 0; c < coneSize; ++c) {
1071: if (debug) {PetscSynchronizedPrintf(PetscObjectComm((PetscObject) dm), "[%d] Cone point %D\n", rank, cone[c]);}
1072: PetscHMapIGet(leafhash, cone[c], &root);
1073: if (!rootdegree[cone[c]] && (root < 0)) {isFace = PETSC_FALSE; break;}
1074: }
1075: if (isFace) {
1076: if (debug) {PetscSynchronizedPrintf(PetscObjectComm((PetscObject) dm), "[%d] Found shared face %D\n", rank, face);}
1077: PetscHMapIJSet(facehash, key, l);
1078: PetscSectionAddDof(candidateSection, localPoint, coneSize);
1079: }
1080: }
1081: }
1082: if (debug) {PetscSynchronizedFlush(PetscObjectComm((PetscObject) dm), NULL);}
1083: PetscHMapIJClear(facehash);
1084: PetscSectionSetUp(candidateSection);
1085: PetscSectionGetStorageSize(candidateSection, &candidatesSize);
1086: PetscMalloc1(candidatesSize, &candidates);
1087: for (l = 0; l < numLeaves; ++l) {
1088: const PetscInt localPoint = localPoints[l];
1089: const PetscInt *support;
1090: PetscInt supportSize, s, offset, idx = 0;
1092: if (debug) {PetscSynchronizedPrintf(PetscObjectComm((PetscObject) dm), "[%d] Checking local point %D\n", rank, localPoint);}
1093: PetscSectionGetOffset(candidateSection, localPoint, &offset);
1094: DMPlexGetSupportSize(dm, localPoint, &supportSize);
1095: DMPlexGetSupport(dm, localPoint, &support);
1096: for (s = 0; s < supportSize; ++s) {
1097: const PetscInt face = support[s];
1098: const PetscInt *cone;
1099: PetscInt coneSize, c, f, root;
1100: PetscBool isFace = PETSC_TRUE;
1102: /* Only add face once */
1103: if (debug) {PetscSynchronizedPrintf(PetscObjectComm((PetscObject) dm), "[%d] Support point %D\n", rank, face);}
1104: key.i = localPoint;
1105: key.j = face;
1106: PetscHMapIJGet(facehash, key, &f);
1107: if (f >= 0) continue;
1108: DMPlexGetConeSize(dm, face, &coneSize);
1109: DMPlexGetCone(dm, face, &cone);
1110: /* If a cone point does not map to leaves on any proc, then do not put face in SF */
1111: for (c = 0; c < coneSize; ++c) {
1112: if (debug) {PetscSynchronizedPrintf(PetscObjectComm((PetscObject) dm), "[%d] Cone point %D\n", rank, cone[c]);}
1113: PetscHMapIGet(leafhash, cone[c], &root);
1114: if (!rootdegree[cone[c]] && (root < 0)) {isFace = PETSC_FALSE; break;}
1115: }
1116: if (isFace) {
1117: if (debug) {PetscSynchronizedPrintf(PetscObjectComm((PetscObject) dm), "[%d] Adding shared face %D at idx %D\n", rank, face, idx);}
1118: PetscHMapIJSet(facehash, key, l);
1119: candidates[offset+idx].rank = -1;
1120: candidates[offset+idx++].index = coneSize-1;
1121: for (c = 0; c < coneSize; ++c) {
1122: if (cone[c] == localPoint) continue;
1123: if (rootdegree[cone[c]]) {
1124: candidates[offset+idx].rank = rank;
1125: candidates[offset+idx++].index = cone[c];
1126: } else {
1127: PetscHMapIGet(leafhash, cone[c], &root);
1128: if (root < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Cannot locate local point %D in SF", cone[c]);
1129: candidates[offset+idx++] = remotePoints[root];
1130: }
1131: }
1132: }
1133: }
1134: }
1135: if (debug) {PetscSynchronizedFlush(PetscObjectComm((PetscObject) dm), NULL);}
1136: PetscHMapIJDestroy(&facehash);
1137: PetscObjectViewFromOptions((PetscObject) candidateSection, NULL, "-petscsection_interp_candidate_view");
1138: SFNodeArrayViewFromOptions(PetscObjectComm((PetscObject) dm), "-petscsection_interp_candidate_view", "Candidates", NULL, candidatesSize, candidates);
1139: }
1140: /* Gather candidate section / array pair into the root partition via inverse(multi(pointSF)). */
1141: /* Note that this section is indexed by offsets into leaves, not by point number */
1142: {
1143: PetscSF sfMulti, sfInverse, sfCandidates;
1144: PetscInt *remoteOffsets;
1146: PetscSFGetMultiSF(pointSF, &sfMulti);
1147: PetscSFCreateInverseSF(sfMulti, &sfInverse);
1148: PetscSectionCreate(PetscObjectComm((PetscObject) dm), &candidateSectionRemote);
1149: PetscSFDistributeSection(sfInverse, candidateSection, &remoteOffsets, candidateSectionRemote);
1150: PetscSFCreateSectionSF(sfInverse, candidateSection, remoteOffsets, candidateSectionRemote, &sfCandidates);
1151: PetscSectionGetStorageSize(candidateSectionRemote, &candidatesRemoteSize);
1152: PetscMalloc1(candidatesRemoteSize, &candidatesRemote);
1153: PetscSFBcastBegin(sfCandidates, MPIU_2INT, candidates, candidatesRemote);
1154: PetscSFBcastEnd(sfCandidates, MPIU_2INT, candidates, candidatesRemote);
1155: PetscSFDestroy(&sfInverse);
1156: PetscSFDestroy(&sfCandidates);
1157: PetscFree(remoteOffsets);
1159: PetscObjectViewFromOptions((PetscObject) candidateSectionRemote, NULL, "-petscsection_interp_candidate_remote_view");
1160: SFNodeArrayViewFromOptions(PetscObjectComm((PetscObject) dm), "-petscsection_interp_candidate_remote_view", "Remote Candidates", NULL, candidatesRemoteSize, candidatesRemote);
1161: }
1162: /* */
1163: {
1164: PetscInt idx;
1165: /* There is a section point for every leaf attached to a given root point */
1166: for (r = 0, idx = 0; r < numRoots; ++r) {
1167: PetscInt deg;
1168: for (deg = 0; deg < rootdegree[r]; ++deg, ++idx) {
1169: PetscInt offset, dof, d;
1171: PetscSectionGetDof(candidateSectionRemote, idx, &dof);
1172: PetscSectionGetOffset(candidateSectionRemote, idx, &offset);
1173: for (d = 0; d < dof; ++d) {
1174: const PetscInt sizeInd = offset+d;
1175: const PetscInt numPoints = candidatesRemote[sizeInd].index;
1176: const PetscInt *join = NULL;
1177: PetscInt points[1024], p, joinSize;
1179: points[0] = r;
1180: for (p = 0; p < numPoints; ++p) {
1181: DMPlexMapToLocalPoint(roothash, localPoints, rank, candidatesRemote[offset+(++d)], &points[p+1]);
1182: if (ierr) {d += numPoints-1 - p; break;} /* We got a point not in our overlap */
1183: if (debug) {PetscSynchronizedPrintf(PetscObjectComm((PetscObject) dm), "[%d] Checking local candidate %D\n", rank, points[p+1]);}
1184: }
1185: if (ierr) continue;
1186: DMPlexGetJoin(dm, numPoints+1, points, &joinSize, &join);
1187: if (joinSize == 1) {
1188: if (debug) {PetscSynchronizedPrintf(PetscObjectComm((PetscObject) dm), "[%d] Adding face %D at idx %D\n", rank, join[0], sizeInd);}
1189: candidatesRemote[sizeInd].rank = rank;
1190: candidatesRemote[sizeInd].index = join[0];
1191: }
1192: DMPlexRestoreJoin(dm, numPoints+1, points, &joinSize, &join);
1193: }
1194: }
1195: }
1196: if (debug) {PetscSynchronizedFlush(PetscObjectComm((PetscObject) dm), NULL);}
1197: }
1198: /* Push claims back to receiver via the MultiSF and derive new pointSF mapping on receiver */
1199: {
1200: PetscSF sfMulti, sfClaims, sfPointNew;
1201: PetscSFNode *remotePointsNew;
1202: PetscHMapI claimshash;
1203: PetscInt *remoteOffsets, *localPointsNew;
1204: PetscInt claimsSize, pStart, pEnd, root, numLocalNew, p, d;
1206: PetscSFGetMultiSF(pointSF, &sfMulti);
1207: PetscSectionCreate(PetscObjectComm((PetscObject) dm), &claimSection);
1208: PetscSFDistributeSection(sfMulti, candidateSectionRemote, &remoteOffsets, claimSection);
1209: PetscSFCreateSectionSF(sfMulti, candidateSectionRemote, remoteOffsets, claimSection, &sfClaims);
1210: PetscSectionGetStorageSize(claimSection, &claimsSize);
1211: PetscMalloc1(claimsSize, &claims);
1212: PetscSFBcastBegin(sfClaims, MPIU_2INT, candidatesRemote, claims);
1213: PetscSFBcastEnd(sfClaims, MPIU_2INT, candidatesRemote, claims);
1214: PetscSFDestroy(&sfClaims);
1215: PetscFree(remoteOffsets);
1216: PetscObjectViewFromOptions((PetscObject) claimSection, NULL, "-petscsection_interp_claim_view");
1217: SFNodeArrayViewFromOptions(PetscObjectComm((PetscObject) dm), "-petscsection_interp_claim_view", "Claims", NULL, claimsSize, claims);
1218: /* Walk the original section of local supports and add an SF entry for each updated item */
1219: PetscHMapICreate(&claimshash);
1220: for (p = 0; p < numRoots; ++p) {
1221: PetscInt dof, offset;
1223: if (debug) {PetscSynchronizedPrintf(PetscObjectComm((PetscObject) dm), "[%d] Checking root for claims %D\n", rank, p);}
1224: PetscSectionGetDof(candidateSection, p, &dof);
1225: PetscSectionGetOffset(candidateSection, p, &offset);
1226: for (d = 0; d < dof;) {
1227: if (claims[offset+d].rank >= 0) {
1228: const PetscInt faceInd = offset+d;
1229: const PetscInt numPoints = candidates[faceInd].index;
1230: const PetscInt *join = NULL;
1231: PetscInt joinSize, points[1024], c;
1233: if (debug) {PetscSynchronizedPrintf(PetscObjectComm((PetscObject) dm), "[%d] Found claim for remote point (%D, %D)\n", rank, claims[faceInd].rank, claims[faceInd].index);}
1234: points[0] = p;
1235: if (debug) {PetscSynchronizedPrintf(PetscObjectComm((PetscObject) dm), "[%d] point %D\n", rank, points[0]);}
1236: for (c = 0, ++d; c < numPoints; ++c, ++d) {
1237: key.i = candidates[offset+d].index;
1238: key.j = candidates[offset+d].rank;
1239: PetscHMapIJGet(roothash, key, &root);
1240: points[c+1] = localPoints[root];
1241: if (debug) {PetscSynchronizedPrintf(PetscObjectComm((PetscObject) dm), "[%d] point %D\n", rank, points[c+1]);}
1242: }
1243: DMPlexGetJoin(dm, numPoints+1, points, &joinSize, &join);
1244: if (joinSize == 1) {
1245: if (debug) {PetscSynchronizedPrintf(PetscObjectComm((PetscObject) dm), "[%d] Found local face %D\n", rank, join[0]);}
1246: PetscHMapISet(claimshash, join[0], faceInd);
1247: }
1248: DMPlexRestoreJoin(dm, numPoints+1, points, &joinSize, &join);
1249: } else d += claims[offset+d].index+1;
1250: }
1251: }
1252: if (debug) {PetscSynchronizedFlush(PetscObjectComm((PetscObject) dm), NULL);}
1253: /* Create new pointSF from hashed claims */
1254: PetscHMapIGetSize(claimshash, &numLocalNew);
1255: DMPlexGetChart(dm, &pStart, &pEnd);
1256: PetscMalloc1(numLeaves + numLocalNew, &localPointsNew);
1257: PetscMalloc1(numLeaves + numLocalNew, &remotePointsNew);
1258: for (p = 0; p < numLeaves; ++p) {
1259: localPointsNew[p] = localPoints[p];
1260: remotePointsNew[p].index = remotePoints[p].index;
1261: remotePointsNew[p].rank = remotePoints[p].rank;
1262: }
1263: p = numLeaves;
1264: PetscHMapIGetKeys(claimshash, &p, localPointsNew);
1265: PetscSortInt(numLocalNew, &localPointsNew[numLeaves]);
1266: for (p = numLeaves; p < numLeaves + numLocalNew; ++p) {
1267: PetscInt offset;
1268: PetscHMapIGet(claimshash, localPointsNew[p], &offset);
1269: remotePointsNew[p] = claims[offset];
1270: }
1271: PetscSFCreate(PetscObjectComm((PetscObject) dm), &sfPointNew);
1272: PetscSFSetGraph(sfPointNew, pEnd-pStart, numLeaves+numLocalNew, localPointsNew, PETSC_OWN_POINTER, remotePointsNew, PETSC_OWN_POINTER);
1273: DMSetPointSF(dm, sfPointNew);
1274: PetscSFDestroy(&sfPointNew);
1275: PetscHMapIDestroy(&claimshash);
1276: }
1277: PetscHMapIDestroy(&leafhash);
1278: PetscHMapIJDestroy(&roothash);
1279: PetscSectionDestroy(&candidateSection);
1280: PetscSectionDestroy(&candidateSectionRemote);
1281: PetscSectionDestroy(&claimSection);
1282: PetscFree(candidates);
1283: PetscFree(candidatesRemote);
1284: PetscFree(claims);
1285: PetscLogEventEnd(DMPLEX_InterpolateSF,dm,0,0,0);
1286: return(0);
1287: }
1289: /*@C
1290: DMPlexInterpolate - Take in a cell-vertex mesh and return one with all intermediate faces, edges, etc.
1292: Collective on DM
1294: Input Parameters:
1295: + dm - The DMPlex object with only cells and vertices
1296: - dmInt - The interpolated DM
1298: Output Parameter:
1299: . dmInt - The complete DMPlex object
1301: Level: intermediate
1303: Notes:
1304: It does not copy over the coordinates.
1306: .keywords: mesh
1307: .seealso: DMPlexUninterpolate(), DMPlexCreateFromCellList(), DMPlexCopyCoordinates()
1308: @*/
1309: PetscErrorCode DMPlexInterpolate(DM dm, DM *dmInt)
1310: {
1311: DM idm, odm = dm;
1312: PetscSF sfPoint;
1313: PetscInt depth, dim, d;
1314: const char *name;
1315: PetscBool flg=PETSC_TRUE;
1321: PetscLogEventBegin(DMPLEX_Interpolate,dm,0,0,0);
1322: DMPlexGetDepth(dm, &depth);
1323: DMGetDimension(dm, &dim);
1324: if ((depth == dim) || (dim <= 1)) {
1325: PetscObjectReference((PetscObject) dm);
1326: idm = dm;
1327: } else {
1328: for (d = 1; d < dim; ++d) {
1329: /* Create interpolated mesh */
1330: DMCreate(PetscObjectComm((PetscObject)dm), &idm);
1331: DMSetType(idm, DMPLEX);
1332: DMSetDimension(idm, dim);
1333: if (depth > 0) {
1334: DMPlexInterpolateFaces_Internal(odm, 1, idm);
1335: DMGetPointSF(odm, &sfPoint);
1336: DMPlexInterpolatePointSF(idm, sfPoint);
1337: }
1338: if (odm != dm) {DMDestroy(&odm);}
1339: odm = idm;
1340: }
1341: PetscObjectGetName((PetscObject) dm, &name);
1342: PetscObjectSetName((PetscObject) idm, name);
1343: DMPlexCopyCoordinates(dm, idm);
1344: DMCopyLabels(dm, idm);
1345: PetscOptionsGetBool(((PetscObject)dm)->options, ((PetscObject)dm)->prefix, "-dm_plex_interpolate_orient_interfaces", &flg, NULL);
1346: if (flg) {DMPlexOrientInterface(idm);}
1347: }
1348: {
1349: PetscBool isper;
1350: const PetscReal *maxCell, *L;
1351: const DMBoundaryType *bd;
1353: DMGetPeriodicity(dm,&isper,&maxCell,&L,&bd);
1354: DMSetPeriodicity(idm,isper,maxCell,L,bd);
1355: }
1356: *dmInt = idm;
1357: PetscLogEventEnd(DMPLEX_Interpolate,dm,0,0,0);
1358: return(0);
1359: }
1361: /*@
1362: DMPlexCopyCoordinates - Copy coordinates from one mesh to another with the same vertices
1364: Collective on DM
1366: Input Parameter:
1367: . dmA - The DMPlex object with initial coordinates
1369: Output Parameter:
1370: . dmB - The DMPlex object with copied coordinates
1372: Level: intermediate
1374: Note: This is typically used when adding pieces other than vertices to a mesh
1376: .keywords: mesh
1377: .seealso: DMCopyLabels(), DMGetCoordinates(), DMGetCoordinatesLocal(), DMGetCoordinateDM(), DMGetCoordinateSection()
1378: @*/
1379: PetscErrorCode DMPlexCopyCoordinates(DM dmA, DM dmB)
1380: {
1381: Vec coordinatesA, coordinatesB;
1382: VecType vtype;
1383: PetscSection coordSectionA, coordSectionB;
1384: PetscScalar *coordsA, *coordsB;
1385: PetscInt spaceDim, Nf, vStartA, vStartB, vEndA, vEndB, coordSizeB, v, d;
1386: PetscInt cStartA, cEndA, cStartB, cEndB, cS, cE;
1387: PetscBool lc = PETSC_FALSE;
1393: if (dmA == dmB) return(0);
1394: DMPlexGetDepthStratum(dmA, 0, &vStartA, &vEndA);
1395: DMPlexGetDepthStratum(dmB, 0, &vStartB, &vEndB);
1396: 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);
1397: DMPlexGetHeightStratum(dmA, 0, &cStartA, &cEndA);
1398: DMPlexGetHeightStratum(dmB, 0, &cStartB, &cEndB);
1399: DMGetCoordinateSection(dmA, &coordSectionA);
1400: DMGetCoordinateSection(dmB, &coordSectionB);
1401: if (coordSectionA == coordSectionB) return(0);
1402: PetscSectionGetNumFields(coordSectionA, &Nf);
1403: if (!Nf) return(0);
1404: if (Nf > 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "The number of coordinate fields must be 1, not %D", Nf);
1405: if (!coordSectionB) {
1406: PetscInt dim;
1408: PetscSectionCreate(PetscObjectComm((PetscObject) coordSectionA), &coordSectionB);
1409: DMGetCoordinateDim(dmA, &dim);
1410: DMSetCoordinateSection(dmB, dim, coordSectionB);
1411: PetscObjectDereference((PetscObject) coordSectionB);
1412: }
1413: PetscSectionSetNumFields(coordSectionB, 1);
1414: PetscSectionGetFieldComponents(coordSectionA, 0, &spaceDim);
1415: PetscSectionSetFieldComponents(coordSectionB, 0, spaceDim);
1416: PetscSectionGetChart(coordSectionA, &cS, &cE);
1417: if (cStartA <= cS && cS < cEndA) { /* localized coordinates */
1418: 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);
1419: cS = cS - cStartA + cStartB;
1420: cE = vEndB;
1421: lc = PETSC_TRUE;
1422: } else {
1423: cS = vStartB;
1424: cE = vEndB;
1425: }
1426: PetscSectionSetChart(coordSectionB, cS, cE);
1427: for (v = vStartB; v < vEndB; ++v) {
1428: PetscSectionSetDof(coordSectionB, v, spaceDim);
1429: PetscSectionSetFieldDof(coordSectionB, v, 0, spaceDim);
1430: }
1431: if (lc) { /* localized coordinates */
1432: PetscInt c;
1434: for (c = cS-cStartB; c < cEndB-cStartB; c++) {
1435: PetscInt dof;
1437: PetscSectionGetDof(coordSectionA, c + cStartA, &dof);
1438: PetscSectionSetDof(coordSectionB, c + cStartB, dof);
1439: PetscSectionSetFieldDof(coordSectionB, c + cStartB, 0, dof);
1440: }
1441: }
1442: PetscSectionSetUp(coordSectionB);
1443: PetscSectionGetStorageSize(coordSectionB, &coordSizeB);
1444: DMGetCoordinatesLocal(dmA, &coordinatesA);
1445: VecCreate(PETSC_COMM_SELF, &coordinatesB);
1446: PetscObjectSetName((PetscObject) coordinatesB, "coordinates");
1447: VecSetSizes(coordinatesB, coordSizeB, PETSC_DETERMINE);
1448: VecGetBlockSize(coordinatesA, &d);
1449: VecSetBlockSize(coordinatesB, d);
1450: VecGetType(coordinatesA, &vtype);
1451: VecSetType(coordinatesB, vtype);
1452: VecGetArray(coordinatesA, &coordsA);
1453: VecGetArray(coordinatesB, &coordsB);
1454: for (v = 0; v < vEndB-vStartB; ++v) {
1455: PetscInt offA, offB;
1457: PetscSectionGetOffset(coordSectionA, v + vStartA, &offA);
1458: PetscSectionGetOffset(coordSectionB, v + vStartB, &offB);
1459: for (d = 0; d < spaceDim; ++d) {
1460: coordsB[offB+d] = coordsA[offA+d];
1461: }
1462: }
1463: if (lc) { /* localized coordinates */
1464: PetscInt c;
1466: for (c = cS-cStartB; c < cEndB-cStartB; c++) {
1467: PetscInt dof, offA, offB;
1469: PetscSectionGetOffset(coordSectionA, c + cStartA, &offA);
1470: PetscSectionGetOffset(coordSectionB, c + cStartB, &offB);
1471: PetscSectionGetDof(coordSectionA, c + cStartA, &dof);
1472: PetscMemcpy(coordsB + offB,coordsA + offA,dof*sizeof(*coordsB));
1473: }
1474: }
1475: VecRestoreArray(coordinatesA, &coordsA);
1476: VecRestoreArray(coordinatesB, &coordsB);
1477: DMSetCoordinatesLocal(dmB, coordinatesB);
1478: VecDestroy(&coordinatesB);
1479: return(0);
1480: }
1482: /*@
1483: DMPlexUninterpolate - Take in a mesh with all intermediate faces, edges, etc. and return a cell-vertex mesh
1485: Collective on DM
1487: Input Parameter:
1488: . dm - The complete DMPlex object
1490: Output Parameter:
1491: . dmUnint - The DMPlex object with only cells and vertices
1493: Level: intermediate
1495: Notes:
1496: It does not copy over the coordinates.
1498: .keywords: mesh
1499: .seealso: DMPlexInterpolate(), DMPlexCreateFromCellList(), DMPlexCopyCoordinates()
1500: @*/
1501: PetscErrorCode DMPlexUninterpolate(DM dm, DM *dmUnint)
1502: {
1503: DM udm;
1504: PetscInt dim, vStart, vEnd, cStart, cEnd, cMax, c, maxConeSize = 0, *cone;
1510: DMGetDimension(dm, &dim);
1511: if (dim <= 1) {
1512: PetscObjectReference((PetscObject) dm);
1513: *dmUnint = dm;
1514: return(0);
1515: }
1516: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
1517: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
1518: DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);
1519: DMCreate(PetscObjectComm((PetscObject) dm), &udm);
1520: DMSetType(udm, DMPLEX);
1521: DMSetDimension(udm, dim);
1522: DMPlexSetChart(udm, cStart, vEnd);
1523: for (c = cStart; c < cEnd; ++c) {
1524: PetscInt *closure = NULL, closureSize, cl, coneSize = 0;
1526: DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
1527: for (cl = 0; cl < closureSize*2; cl += 2) {
1528: const PetscInt p = closure[cl];
1530: if ((p >= vStart) && (p < vEnd)) ++coneSize;
1531: }
1532: DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
1533: DMPlexSetConeSize(udm, c, coneSize);
1534: maxConeSize = PetscMax(maxConeSize, coneSize);
1535: }
1536: DMSetUp(udm);
1537: PetscMalloc1(maxConeSize, &cone);
1538: for (c = cStart; c < cEnd; ++c) {
1539: PetscInt *closure = NULL, closureSize, cl, coneSize = 0;
1541: DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
1542: for (cl = 0; cl < closureSize*2; cl += 2) {
1543: const PetscInt p = closure[cl];
1545: if ((p >= vStart) && (p < vEnd)) cone[coneSize++] = p;
1546: }
1547: DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
1548: DMPlexSetCone(udm, c, cone);
1549: }
1550: PetscFree(cone);
1551: DMPlexSetHybridBounds(udm, cMax, PETSC_DETERMINE, PETSC_DETERMINE, PETSC_DETERMINE);
1552: DMPlexSymmetrize(udm);
1553: DMPlexStratify(udm);
1554: /* Reduce SF */
1555: {
1556: PetscSF sfPoint, sfPointUn;
1557: const PetscSFNode *remotePoints;
1558: const PetscInt *localPoints;
1559: PetscSFNode *remotePointsUn;
1560: PetscInt *localPointsUn;
1561: PetscInt vEnd, numRoots, numLeaves, l;
1562: PetscInt numLeavesUn = 0, n = 0;
1563: PetscErrorCode ierr;
1565: /* Get original SF information */
1566: DMGetPointSF(dm, &sfPoint);
1567: DMGetPointSF(udm, &sfPointUn);
1568: DMPlexGetDepthStratum(dm, 0, NULL, &vEnd);
1569: PetscSFGetGraph(sfPoint, &numRoots, &numLeaves, &localPoints, &remotePoints);
1570: /* Allocate space for cells and vertices */
1571: for (l = 0; l < numLeaves; ++l) if (localPoints[l] < vEnd) numLeavesUn++;
1572: /* Fill in leaves */
1573: if (vEnd >= 0) {
1574: PetscMalloc1(numLeavesUn, &remotePointsUn);
1575: PetscMalloc1(numLeavesUn, &localPointsUn);
1576: for (l = 0; l < numLeaves; l++) {
1577: if (localPoints[l] < vEnd) {
1578: localPointsUn[n] = localPoints[l];
1579: remotePointsUn[n].rank = remotePoints[l].rank;
1580: remotePointsUn[n].index = remotePoints[l].index;
1581: ++n;
1582: }
1583: }
1584: if (n != numLeavesUn) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Inconsistent number of leaves %d != %d", n, numLeavesUn);
1585: PetscSFSetGraph(sfPointUn, vEnd, numLeavesUn, localPointsUn, PETSC_OWN_POINTER, remotePointsUn, PETSC_OWN_POINTER);
1586: }
1587: }
1588: {
1589: PetscBool isper;
1590: const PetscReal *maxCell, *L;
1591: const DMBoundaryType *bd;
1593: DMGetPeriodicity(dm,&isper,&maxCell,&L,&bd);
1594: DMSetPeriodicity(udm,isper,maxCell,L,bd);
1595: }
1597: *dmUnint = udm;
1598: return(0);
1599: }