Actual source code: plexinterpolate.c
1: #include <petsc/private/dmpleximpl.h>
2: #include <petsc/private/hashmapi.h>
3: #include <petsc/private/hashmapij.h>
5: const char *const DMPlexInterpolatedFlags[] = {"none", "partial", "mixed", "full", "DMPlexInterpolatedFlag", "DMPLEX_INTERPOLATED_", NULL};
7: /* HMapIJKL */
9: #include <petsc/private/hashmapijkl.h>
11: static PetscSFNode _PetscInvalidSFNode = {-1, -1};
13: typedef struct _PetscHMapIJKLRemoteKey {
14: PetscSFNode i, j, k, l;
15: } PetscHMapIJKLRemoteKey;
17: #define PetscHMapIJKLRemoteKeyHash(key) \
18: PetscHashCombine(PetscHashCombine(PetscHashInt((key).i.rank + (key).i.index), PetscHashInt((key).j.rank + (key).j.index)), PetscHashCombine(PetscHashInt((key).k.rank + (key).k.index), PetscHashInt((key).l.rank + (key).l.index)))
20: #define PetscHMapIJKLRemoteKeyEqual(k1, k2) \
21: (((k1).i.rank == (k2).i.rank) ? ((k1).i.index == (k2).i.index) ? ((k1).j.rank == (k2).j.rank) ? ((k1).j.index == (k2).j.index) ? ((k1).k.rank == (k2).k.rank) ? ((k1).k.index == (k2).k.index) ? ((k1).l.rank == (k2).l.rank) ? ((k1).l.index == (k2).l.index) : 0 : 0 : 0 : 0 : 0 : 0 : 0)
23: PetscDisableStaticAnalyzerForExpressionUnderstandingThatThisIsDangerousAndBugprone(PETSC_HASH_MAP(HMapIJKLRemote, PetscHMapIJKLRemoteKey, PetscSFNode, PetscHMapIJKLRemoteKeyHash, PetscHMapIJKLRemoteKeyEqual, _PetscInvalidSFNode))
25: static PetscErrorCode PetscSortSFNode(PetscInt n, PetscSFNode A[])
26: {
27: PetscInt i;
29: PetscFunctionBegin;
30: for (i = 1; i < n; ++i) {
31: PetscSFNode x = A[i];
32: PetscInt j;
34: for (j = i - 1; j >= 0; --j) {
35: if ((A[j].rank > x.rank) || (A[j].rank == x.rank && A[j].index > x.index)) break;
36: A[j + 1] = A[j];
37: }
38: A[j + 1] = x;
39: }
40: PetscFunctionReturn(PETSC_SUCCESS);
41: }
43: /*
44: DMPlexGetRawFaces_Internal - Gets groups of vertices that correspond to faces for the given cone
45: */
46: PetscErrorCode DMPlexGetRawFaces_Internal(DM dm, DMPolytopeType ct, const PetscInt cone[], PetscInt *numFaces, const DMPolytopeType *faceTypes[], const PetscInt *faceSizes[], const PetscInt *faces[])
47: {
48: DMPolytopeType *typesTmp = NULL;
49: PetscInt *sizesTmp = NULL, *facesTmp = NULL;
50: PetscInt *tmp;
51: PetscInt maxConeSize, maxSupportSize, maxSize;
52: PetscInt getSize = 0;
54: PetscFunctionBegin;
56: if (cone) PetscAssertPointer(cone, 3);
57: PetscCall(DMPlexGetMaxSizes(dm, &maxConeSize, &maxSupportSize));
58: maxSize = PetscMax(maxConeSize, maxSupportSize);
59: if (faceTypes) getSize += maxSize;
60: if (faceSizes) getSize += maxSize;
61: if (faces) getSize += PetscSqr(maxSize);
62: PetscCall(DMGetWorkArray(dm, getSize, MPIU_INT, &tmp));
63: if (faceTypes) {
64: typesTmp = (DMPolytopeType *)tmp;
65: tmp += maxSize;
66: }
67: if (faceSizes) {
68: sizesTmp = tmp;
69: tmp += maxSize;
70: }
71: if (faces) facesTmp = tmp;
72: switch (ct) {
73: case DM_POLYTOPE_POINT:
74: if (numFaces) *numFaces = 0;
75: break;
76: case DM_POLYTOPE_SEGMENT:
77: if (numFaces) *numFaces = 2;
78: if (faceTypes) {
79: typesTmp[0] = DM_POLYTOPE_POINT;
80: typesTmp[1] = DM_POLYTOPE_POINT;
81: *faceTypes = typesTmp;
82: }
83: if (faceSizes) {
84: sizesTmp[0] = 1;
85: sizesTmp[1] = 1;
86: *faceSizes = sizesTmp;
87: }
88: if (faces) {
89: facesTmp[0] = cone[0];
90: facesTmp[1] = cone[1];
91: *faces = facesTmp;
92: }
93: break;
94: case DM_POLYTOPE_POINT_PRISM_TENSOR:
95: if (numFaces) *numFaces = 2;
96: if (faceTypes) {
97: typesTmp[0] = DM_POLYTOPE_POINT;
98: typesTmp[1] = DM_POLYTOPE_POINT;
99: *faceTypes = typesTmp;
100: }
101: if (faceSizes) {
102: sizesTmp[0] = 1;
103: sizesTmp[1] = 1;
104: *faceSizes = sizesTmp;
105: }
106: if (faces) {
107: facesTmp[0] = cone[0];
108: facesTmp[1] = cone[1];
109: *faces = facesTmp;
110: }
111: break;
112: case DM_POLYTOPE_TRIANGLE:
113: if (numFaces) *numFaces = 3;
114: if (faceTypes) {
115: typesTmp[0] = DM_POLYTOPE_SEGMENT;
116: typesTmp[1] = DM_POLYTOPE_SEGMENT;
117: typesTmp[2] = DM_POLYTOPE_SEGMENT;
118: *faceTypes = typesTmp;
119: }
120: if (faceSizes) {
121: sizesTmp[0] = 2;
122: sizesTmp[1] = 2;
123: sizesTmp[2] = 2;
124: *faceSizes = sizesTmp;
125: }
126: if (faces) {
127: facesTmp[0] = cone[0];
128: facesTmp[1] = cone[1];
129: facesTmp[2] = cone[1];
130: facesTmp[3] = cone[2];
131: facesTmp[4] = cone[2];
132: facesTmp[5] = cone[0];
133: *faces = facesTmp;
134: }
135: break;
136: case DM_POLYTOPE_QUADRILATERAL:
137: /* Vertices follow right hand rule */
138: if (numFaces) *numFaces = 4;
139: if (faceTypes) {
140: typesTmp[0] = DM_POLYTOPE_SEGMENT;
141: typesTmp[1] = DM_POLYTOPE_SEGMENT;
142: typesTmp[2] = DM_POLYTOPE_SEGMENT;
143: typesTmp[3] = DM_POLYTOPE_SEGMENT;
144: *faceTypes = typesTmp;
145: }
146: if (faceSizes) {
147: sizesTmp[0] = 2;
148: sizesTmp[1] = 2;
149: sizesTmp[2] = 2;
150: sizesTmp[3] = 2;
151: *faceSizes = sizesTmp;
152: }
153: if (faces) {
154: facesTmp[0] = cone[0];
155: facesTmp[1] = cone[1];
156: facesTmp[2] = cone[1];
157: facesTmp[3] = cone[2];
158: facesTmp[4] = cone[2];
159: facesTmp[5] = cone[3];
160: facesTmp[6] = cone[3];
161: facesTmp[7] = cone[0];
162: *faces = facesTmp;
163: }
164: break;
165: case DM_POLYTOPE_SEG_PRISM_TENSOR:
166: if (numFaces) *numFaces = 4;
167: if (faceTypes) {
168: typesTmp[0] = DM_POLYTOPE_SEGMENT;
169: typesTmp[1] = DM_POLYTOPE_SEGMENT;
170: typesTmp[2] = DM_POLYTOPE_POINT_PRISM_TENSOR;
171: typesTmp[3] = DM_POLYTOPE_POINT_PRISM_TENSOR;
172: *faceTypes = typesTmp;
173: }
174: if (faceSizes) {
175: sizesTmp[0] = 2;
176: sizesTmp[1] = 2;
177: sizesTmp[2] = 2;
178: sizesTmp[3] = 2;
179: *faceSizes = sizesTmp;
180: }
181: if (faces) {
182: facesTmp[0] = cone[0];
183: facesTmp[1] = cone[1];
184: facesTmp[2] = cone[2];
185: facesTmp[3] = cone[3];
186: facesTmp[4] = cone[0];
187: facesTmp[5] = cone[2];
188: facesTmp[6] = cone[1];
189: facesTmp[7] = cone[3];
190: *faces = facesTmp;
191: }
192: break;
193: case DM_POLYTOPE_TETRAHEDRON:
194: /* Vertices of first face follow right hand rule and normal points away from last vertex */
195: if (numFaces) *numFaces = 4;
196: if (faceTypes) {
197: typesTmp[0] = DM_POLYTOPE_TRIANGLE;
198: typesTmp[1] = DM_POLYTOPE_TRIANGLE;
199: typesTmp[2] = DM_POLYTOPE_TRIANGLE;
200: typesTmp[3] = DM_POLYTOPE_TRIANGLE;
201: *faceTypes = typesTmp;
202: }
203: if (faceSizes) {
204: sizesTmp[0] = 3;
205: sizesTmp[1] = 3;
206: sizesTmp[2] = 3;
207: sizesTmp[3] = 3;
208: *faceSizes = sizesTmp;
209: }
210: if (faces) {
211: facesTmp[0] = cone[0];
212: facesTmp[1] = cone[1];
213: facesTmp[2] = cone[2];
214: facesTmp[3] = cone[0];
215: facesTmp[4] = cone[3];
216: facesTmp[5] = cone[1];
217: facesTmp[6] = cone[0];
218: facesTmp[7] = cone[2];
219: facesTmp[8] = cone[3];
220: facesTmp[9] = cone[2];
221: facesTmp[10] = cone[1];
222: facesTmp[11] = cone[3];
223: *faces = facesTmp;
224: }
225: break;
226: case DM_POLYTOPE_HEXAHEDRON:
227: /* 7--------6
228: /| /|
229: / | / |
230: 4--------5 |
231: | | | |
232: | | | |
233: | 1--------2
234: | / | /
235: |/ |/
236: 0--------3
237: */
238: if (numFaces) *numFaces = 6;
239: if (faceTypes) {
240: typesTmp[0] = DM_POLYTOPE_QUADRILATERAL;
241: typesTmp[1] = DM_POLYTOPE_QUADRILATERAL;
242: typesTmp[2] = DM_POLYTOPE_QUADRILATERAL;
243: typesTmp[3] = DM_POLYTOPE_QUADRILATERAL;
244: typesTmp[4] = DM_POLYTOPE_QUADRILATERAL;
245: typesTmp[5] = DM_POLYTOPE_QUADRILATERAL;
246: *faceTypes = typesTmp;
247: }
248: if (faceSizes) {
249: sizesTmp[0] = 4;
250: sizesTmp[1] = 4;
251: sizesTmp[2] = 4;
252: sizesTmp[3] = 4;
253: sizesTmp[4] = 4;
254: sizesTmp[5] = 4;
255: *faceSizes = sizesTmp;
256: }
257: if (faces) {
258: facesTmp[0] = cone[0];
259: facesTmp[1] = cone[1];
260: facesTmp[2] = cone[2];
261: facesTmp[3] = cone[3]; /* Bottom */
262: facesTmp[4] = cone[4];
263: facesTmp[5] = cone[5];
264: facesTmp[6] = cone[6];
265: facesTmp[7] = cone[7]; /* Top */
266: facesTmp[8] = cone[0];
267: facesTmp[9] = cone[3];
268: facesTmp[10] = cone[5];
269: facesTmp[11] = cone[4]; /* Front */
270: facesTmp[12] = cone[2];
271: facesTmp[13] = cone[1];
272: facesTmp[14] = cone[7];
273: facesTmp[15] = cone[6]; /* Back */
274: facesTmp[16] = cone[3];
275: facesTmp[17] = cone[2];
276: facesTmp[18] = cone[6];
277: facesTmp[19] = cone[5]; /* Right */
278: facesTmp[20] = cone[0];
279: facesTmp[21] = cone[4];
280: facesTmp[22] = cone[7];
281: facesTmp[23] = cone[1]; /* Left */
282: *faces = facesTmp;
283: }
284: break;
285: case DM_POLYTOPE_TRI_PRISM:
286: if (numFaces) *numFaces = 5;
287: if (faceTypes) {
288: typesTmp[0] = DM_POLYTOPE_TRIANGLE;
289: typesTmp[1] = DM_POLYTOPE_TRIANGLE;
290: typesTmp[2] = DM_POLYTOPE_QUADRILATERAL;
291: typesTmp[3] = DM_POLYTOPE_QUADRILATERAL;
292: typesTmp[4] = DM_POLYTOPE_QUADRILATERAL;
293: *faceTypes = typesTmp;
294: }
295: if (faceSizes) {
296: sizesTmp[0] = 3;
297: sizesTmp[1] = 3;
298: sizesTmp[2] = 4;
299: sizesTmp[3] = 4;
300: sizesTmp[4] = 4;
301: *faceSizes = sizesTmp;
302: }
303: if (faces) {
304: facesTmp[0] = cone[0];
305: facesTmp[1] = cone[1];
306: facesTmp[2] = cone[2]; /* Bottom */
307: facesTmp[3] = cone[3];
308: facesTmp[4] = cone[4];
309: facesTmp[5] = cone[5]; /* Top */
310: facesTmp[6] = cone[0];
311: facesTmp[7] = cone[2];
312: facesTmp[8] = cone[4];
313: facesTmp[9] = cone[3]; /* Back left */
314: facesTmp[10] = cone[2];
315: facesTmp[11] = cone[1];
316: facesTmp[12] = cone[5];
317: facesTmp[13] = cone[4]; /* Front */
318: facesTmp[14] = cone[1];
319: facesTmp[15] = cone[0];
320: facesTmp[16] = cone[3];
321: facesTmp[17] = cone[5]; /* Back right */
322: *faces = facesTmp;
323: }
324: break;
325: case DM_POLYTOPE_TRI_PRISM_TENSOR:
326: if (numFaces) *numFaces = 5;
327: if (faceTypes) {
328: typesTmp[0] = DM_POLYTOPE_TRIANGLE;
329: typesTmp[1] = DM_POLYTOPE_TRIANGLE;
330: typesTmp[2] = DM_POLYTOPE_SEG_PRISM_TENSOR;
331: typesTmp[3] = DM_POLYTOPE_SEG_PRISM_TENSOR;
332: typesTmp[4] = DM_POLYTOPE_SEG_PRISM_TENSOR;
333: *faceTypes = typesTmp;
334: }
335: if (faceSizes) {
336: sizesTmp[0] = 3;
337: sizesTmp[1] = 3;
338: sizesTmp[2] = 4;
339: sizesTmp[3] = 4;
340: sizesTmp[4] = 4;
341: *faceSizes = sizesTmp;
342: }
343: if (faces) {
344: facesTmp[0] = cone[0];
345: facesTmp[1] = cone[1];
346: facesTmp[2] = cone[2]; /* Bottom */
347: facesTmp[3] = cone[3];
348: facesTmp[4] = cone[4];
349: facesTmp[5] = cone[5]; /* Top */
350: facesTmp[6] = cone[0];
351: facesTmp[7] = cone[1];
352: facesTmp[8] = cone[3];
353: facesTmp[9] = cone[4]; /* Back left */
354: facesTmp[10] = cone[1];
355: facesTmp[11] = cone[2];
356: facesTmp[12] = cone[4];
357: facesTmp[13] = cone[5]; /* Back right */
358: facesTmp[14] = cone[2];
359: facesTmp[15] = cone[0];
360: facesTmp[16] = cone[5];
361: facesTmp[17] = cone[3]; /* Front */
362: *faces = facesTmp;
363: }
364: break;
365: case DM_POLYTOPE_QUAD_PRISM_TENSOR:
366: /* 7--------6
367: /| /|
368: / | / |
369: 4--------5 |
370: | | | |
371: | | | |
372: | 3--------2
373: | / | /
374: |/ |/
375: 0--------1
376: */
377: if (numFaces) *numFaces = 6;
378: if (faceTypes) {
379: typesTmp[0] = DM_POLYTOPE_QUADRILATERAL;
380: typesTmp[1] = DM_POLYTOPE_QUADRILATERAL;
381: typesTmp[2] = DM_POLYTOPE_SEG_PRISM_TENSOR;
382: typesTmp[3] = DM_POLYTOPE_SEG_PRISM_TENSOR;
383: typesTmp[4] = DM_POLYTOPE_SEG_PRISM_TENSOR;
384: typesTmp[5] = DM_POLYTOPE_SEG_PRISM_TENSOR;
385: *faceTypes = typesTmp;
386: }
387: if (faceSizes) {
388: sizesTmp[0] = 4;
389: sizesTmp[1] = 4;
390: sizesTmp[2] = 4;
391: sizesTmp[3] = 4;
392: sizesTmp[4] = 4;
393: sizesTmp[5] = 4;
394: *faceSizes = sizesTmp;
395: }
396: if (faces) {
397: facesTmp[0] = cone[0];
398: facesTmp[1] = cone[1];
399: facesTmp[2] = cone[2];
400: facesTmp[3] = cone[3]; /* Bottom */
401: facesTmp[4] = cone[4];
402: facesTmp[5] = cone[5];
403: facesTmp[6] = cone[6];
404: facesTmp[7] = cone[7]; /* Top */
405: facesTmp[8] = cone[0];
406: facesTmp[9] = cone[1];
407: facesTmp[10] = cone[4];
408: facesTmp[11] = cone[5]; /* Front */
409: facesTmp[12] = cone[1];
410: facesTmp[13] = cone[2];
411: facesTmp[14] = cone[5];
412: facesTmp[15] = cone[6]; /* Right */
413: facesTmp[16] = cone[2];
414: facesTmp[17] = cone[3];
415: facesTmp[18] = cone[6];
416: facesTmp[19] = cone[7]; /* Back */
417: facesTmp[20] = cone[3];
418: facesTmp[21] = cone[0];
419: facesTmp[22] = cone[7];
420: facesTmp[23] = cone[4]; /* Left */
421: *faces = facesTmp;
422: }
423: break;
424: case DM_POLYTOPE_PYRAMID:
425: /*
426: 4----
427: |\-\ \-----
428: | \ -\ \
429: | 1--\-----2
430: | / \ /
431: |/ \ /
432: 0--------3
433: */
434: if (numFaces) *numFaces = 5;
435: if (faceTypes) {
436: typesTmp[0] = DM_POLYTOPE_QUADRILATERAL;
437: typesTmp[1] = DM_POLYTOPE_TRIANGLE;
438: typesTmp[2] = DM_POLYTOPE_TRIANGLE;
439: typesTmp[3] = DM_POLYTOPE_TRIANGLE;
440: typesTmp[4] = DM_POLYTOPE_TRIANGLE;
441: *faceTypes = typesTmp;
442: }
443: if (faceSizes) {
444: sizesTmp[0] = 4;
445: sizesTmp[1] = 3;
446: sizesTmp[2] = 3;
447: sizesTmp[3] = 3;
448: sizesTmp[4] = 3;
449: *faceSizes = sizesTmp;
450: }
451: if (faces) {
452: facesTmp[0] = cone[0];
453: facesTmp[1] = cone[1];
454: facesTmp[2] = cone[2];
455: facesTmp[3] = cone[3]; /* Bottom */
456: facesTmp[4] = cone[0];
457: facesTmp[5] = cone[3];
458: facesTmp[6] = cone[4]; /* Front */
459: facesTmp[7] = cone[3];
460: facesTmp[8] = cone[2];
461: facesTmp[9] = cone[4]; /* Right */
462: facesTmp[10] = cone[2];
463: facesTmp[11] = cone[1];
464: facesTmp[12] = cone[4]; /* Back */
465: facesTmp[13] = cone[1];
466: facesTmp[14] = cone[0];
467: facesTmp[15] = cone[4]; /* Left */
468: *faces = facesTmp;
469: }
470: break;
471: default:
472: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "No face description for cell type %s", DMPolytopeTypes[ct]);
473: }
474: PetscFunctionReturn(PETSC_SUCCESS);
475: }
477: PetscErrorCode DMPlexRestoreRawFaces_Internal(DM dm, DMPolytopeType ct, const PetscInt cone[], PetscInt *numFaces, const DMPolytopeType *faceTypes[], const PetscInt *faceSizes[], const PetscInt *faces[])
478: {
479: PetscFunctionBegin;
480: if (faceTypes) PetscCall(DMRestoreWorkArray(dm, 0, MPIU_INT, (void *)faceTypes));
481: else if (faceSizes) PetscCall(DMRestoreWorkArray(dm, 0, MPIU_INT, (void *)faceSizes));
482: else if (faces) PetscCall(DMRestoreWorkArray(dm, 0, MPIU_INT, (void *)faces));
483: if (faceTypes) *faceTypes = NULL;
484: if (faceSizes) *faceSizes = NULL;
485: if (faces) *faces = NULL;
486: PetscFunctionReturn(PETSC_SUCCESS);
487: }
489: /* This interpolates faces for cells at some stratum */
490: static PetscErrorCode DMPlexInterpolateFaces_Internal(DM dm, PetscInt cellDepth, DM idm)
491: {
492: DMLabel ctLabel;
493: PetscHMapIJKL faceTable;
494: PetscInt faceTypeNum[DM_NUM_POLYTOPES];
495: PetscInt depth, d, pStart, Np, cStart, cEnd, c, fStart, fEnd;
496: PetscInt cntFaces, *facesId, minCone;
498: PetscFunctionBegin;
499: PetscCall(DMPlexGetDepth(dm, &depth));
500: PetscCall(PetscHMapIJKLCreate(&faceTable));
501: PetscCall(PetscArrayzero(faceTypeNum, DM_NUM_POLYTOPES));
502: PetscCall(DMPlexGetDepthStratum(dm, cellDepth, &cStart, &cEnd));
503: /* Number new faces and save face vertices in hash table */
504: PetscCall(DMPlexGetDepthStratum(dm, depth > cellDepth ? cellDepth : 0, NULL, &fStart));
505: fEnd = fStart;
507: minCone = PETSC_MAX_INT;
508: for (c = cStart, cntFaces = 0; c < cEnd; ++c) {
509: const PetscInt *cone;
510: DMPolytopeType ct;
511: PetscInt numFaces = 0, coneSize;
513: PetscCall(DMPlexGetCellType(dm, c, &ct));
514: PetscCall(DMPlexGetCone(dm, c, &cone));
515: PetscCall(DMPlexGetConeSize(dm, c, &coneSize));
516: for (PetscInt j = 0; j < coneSize; j++) minCone = PetscMin(cone[j], minCone);
517: PetscCall(DMPlexGetRawFaces_Internal(dm, ct, cone, &numFaces, NULL, NULL, NULL));
518: cntFaces += numFaces;
519: }
520: // Encode so that we can use 0 as an excluded value, instead of PETSC_MAX_INT
521: minCone = -(minCone - 1);
523: PetscCall(PetscMalloc1(cntFaces, &facesId));
525: for (c = cStart, cntFaces = 0; c < cEnd; ++c) {
526: const PetscInt *cone, *faceSizes, *faces;
527: const DMPolytopeType *faceTypes;
528: DMPolytopeType ct;
529: PetscInt numFaces, cf, foff = 0;
531: PetscCall(DMPlexGetCellType(dm, c, &ct));
532: PetscCall(DMPlexGetCone(dm, c, &cone));
533: PetscCall(DMPlexGetRawFaces_Internal(dm, ct, cone, &numFaces, &faceTypes, &faceSizes, &faces));
534: for (cf = 0; cf < numFaces; foff += faceSizes[cf], ++cf) {
535: const PetscInt faceSize = faceSizes[cf];
536: const DMPolytopeType faceType = faceTypes[cf];
537: const PetscInt *face = &faces[foff];
538: PetscHashIJKLKey key;
539: PetscHashIter iter;
540: PetscBool missing;
542: PetscCheck(faceSize <= 4, PETSC_COMM_SELF, PETSC_ERR_SUP, "Do not support faces of size %" PetscInt_FMT " > 4", faceSize);
543: key.i = face[0] + minCone;
544: key.j = faceSize > 1 ? face[1] + minCone : 0;
545: key.k = faceSize > 2 ? face[2] + minCone : 0;
546: key.l = faceSize > 3 ? face[3] + minCone : 0;
547: PetscCall(PetscSortInt(faceSize, (PetscInt *)&key));
548: PetscCall(PetscHMapIJKLPut(faceTable, key, &iter, &missing));
549: if (missing) {
550: facesId[cntFaces] = fEnd;
551: PetscCall(PetscHMapIJKLIterSet(faceTable, iter, fEnd++));
552: ++faceTypeNum[faceType];
553: } else PetscCall(PetscHMapIJKLIterGet(faceTable, iter, &facesId[cntFaces]));
554: cntFaces++;
555: }
556: PetscCall(DMPlexRestoreRawFaces_Internal(dm, ct, cone, &numFaces, &faceTypes, &faceSizes, &faces));
557: }
558: /* We need to number faces contiguously among types */
559: {
560: PetscInt faceTypeStart[DM_NUM_POLYTOPES], ct, numFT = 0;
562: for (ct = 0; ct < DM_NUM_POLYTOPES; ++ct) {
563: if (faceTypeNum[ct]) ++numFT;
564: faceTypeStart[ct] = 0;
565: }
566: if (numFT > 1) {
567: PetscCall(PetscHMapIJKLClear(faceTable));
568: faceTypeStart[0] = fStart;
569: for (ct = 1; ct < DM_NUM_POLYTOPES; ++ct) faceTypeStart[ct] = faceTypeStart[ct - 1] + faceTypeNum[ct - 1];
570: for (c = cStart, cntFaces = 0; c < cEnd; ++c) {
571: const PetscInt *cone, *faceSizes, *faces;
572: const DMPolytopeType *faceTypes;
573: DMPolytopeType ct;
574: PetscInt numFaces, cf, foff = 0;
576: PetscCall(DMPlexGetCellType(dm, c, &ct));
577: PetscCall(DMPlexGetCone(dm, c, &cone));
578: PetscCall(DMPlexGetRawFaces_Internal(dm, ct, cone, &numFaces, &faceTypes, &faceSizes, &faces));
579: for (cf = 0; cf < numFaces; foff += faceSizes[cf], ++cf) {
580: const PetscInt faceSize = faceSizes[cf];
581: const DMPolytopeType faceType = faceTypes[cf];
582: const PetscInt *face = &faces[foff];
583: PetscHashIJKLKey key;
584: PetscHashIter iter;
585: PetscBool missing;
587: key.i = face[0] + minCone;
588: key.j = faceSize > 1 ? face[1] + minCone : 0;
589: key.k = faceSize > 2 ? face[2] + minCone : 0;
590: key.l = faceSize > 3 ? face[3] + minCone : 0;
591: PetscCall(PetscSortInt(faceSize, (PetscInt *)&key));
592: PetscCall(PetscHMapIJKLPut(faceTable, key, &iter, &missing));
593: if (missing) {
594: facesId[cntFaces] = faceTypeStart[faceType];
595: PetscCall(PetscHMapIJKLIterSet(faceTable, iter, faceTypeStart[faceType]++));
596: } else PetscCall(PetscHMapIJKLIterGet(faceTable, iter, &facesId[cntFaces]));
597: cntFaces++;
598: }
599: PetscCall(DMPlexRestoreRawFaces_Internal(dm, ct, cone, &numFaces, &faceTypes, &faceSizes, &faces));
600: }
601: for (ct = 1; ct < DM_NUM_POLYTOPES; ++ct) {
602: PetscCheck(faceTypeStart[ct] == faceTypeStart[ct - 1] + faceTypeNum[ct], PETSC_COMM_SELF, PETSC_ERR_PLIB, "Inconsistent numbering for cell type %s, %" PetscInt_FMT " != %" PetscInt_FMT " + %" PetscInt_FMT, DMPolytopeTypes[ct], faceTypeStart[ct], faceTypeStart[ct - 1], faceTypeNum[ct]);
603: }
604: }
605: }
606: PetscCall(PetscHMapIJKLDestroy(&faceTable));
608: /* Add new points, always at the end of the numbering */
609: PetscCall(DMPlexGetChart(dm, &pStart, &Np));
610: PetscCall(DMPlexSetChart(idm, pStart, Np + (fEnd - fStart)));
611: /* Set cone sizes */
612: /* Must create the celltype label here so that we do not automatically try to compute the types */
613: PetscCall(DMCreateLabel(idm, "celltype"));
614: PetscCall(DMPlexGetCellTypeLabel(idm, &ctLabel));
615: for (d = 0; d <= depth; ++d) {
616: DMPolytopeType ct;
617: PetscInt coneSize, pStart, pEnd, p;
619: if (d == cellDepth) continue;
620: PetscCall(DMPlexGetDepthStratum(dm, d, &pStart, &pEnd));
621: for (p = pStart; p < pEnd; ++p) {
622: PetscCall(DMPlexGetConeSize(dm, p, &coneSize));
623: PetscCall(DMPlexSetConeSize(idm, p, coneSize));
624: PetscCall(DMPlexGetCellType(dm, p, &ct));
625: PetscCall(DMPlexSetCellType(idm, p, ct));
626: }
627: }
628: for (c = cStart, cntFaces = 0; c < cEnd; ++c) {
629: const PetscInt *cone, *faceSizes;
630: const DMPolytopeType *faceTypes;
631: DMPolytopeType ct;
632: PetscInt numFaces, cf;
634: PetscCall(DMPlexGetCellType(dm, c, &ct));
635: PetscCall(DMPlexGetCone(dm, c, &cone));
636: PetscCall(DMPlexGetRawFaces_Internal(dm, ct, cone, &numFaces, &faceTypes, &faceSizes, NULL));
637: PetscCall(DMPlexSetCellType(idm, c, ct));
638: PetscCall(DMPlexSetConeSize(idm, c, numFaces));
639: for (cf = 0; cf < numFaces; ++cf) {
640: const PetscInt f = facesId[cntFaces];
641: DMPolytopeType faceType = faceTypes[cf];
642: const PetscInt faceSize = faceSizes[cf];
643: PetscCall(DMPlexSetConeSize(idm, f, faceSize));
644: PetscCall(DMPlexSetCellType(idm, f, faceType));
645: cntFaces++;
646: }
647: PetscCall(DMPlexRestoreRawFaces_Internal(dm, ct, cone, &numFaces, &faceTypes, &faceSizes, NULL));
648: }
649: PetscCall(DMSetUp(idm));
650: /* Initialize cones so we do not need the bash table to tell us that a cone has been set */
651: {
652: PetscSection cs;
653: PetscInt *cones, csize;
655: PetscCall(DMPlexGetConeSection(idm, &cs));
656: PetscCall(DMPlexGetCones(idm, &cones));
657: PetscCall(PetscSectionGetStorageSize(cs, &csize));
658: for (c = 0; c < csize; ++c) cones[c] = -1;
659: }
660: /* Set cones */
661: for (d = 0; d <= depth; ++d) {
662: const PetscInt *cone;
663: PetscInt pStart, pEnd, p;
665: if (d == cellDepth) continue;
666: PetscCall(DMPlexGetDepthStratum(dm, d, &pStart, &pEnd));
667: for (p = pStart; p < pEnd; ++p) {
668: PetscCall(DMPlexGetCone(dm, p, &cone));
669: PetscCall(DMPlexSetCone(idm, p, cone));
670: PetscCall(DMPlexGetConeOrientation(dm, p, &cone));
671: PetscCall(DMPlexSetConeOrientation(idm, p, cone));
672: }
673: }
674: for (c = cStart, cntFaces = 0; c < cEnd; ++c) {
675: const PetscInt *cone, *faceSizes, *faces;
676: const DMPolytopeType *faceTypes;
677: DMPolytopeType ct;
678: PetscInt numFaces, cf, foff = 0;
680: PetscCall(DMPlexGetCellType(dm, c, &ct));
681: PetscCall(DMPlexGetCone(dm, c, &cone));
682: PetscCall(DMPlexGetRawFaces_Internal(dm, ct, cone, &numFaces, &faceTypes, &faceSizes, &faces));
683: for (cf = 0; cf < numFaces; foff += faceSizes[cf], ++cf) {
684: DMPolytopeType faceType = faceTypes[cf];
685: const PetscInt faceSize = faceSizes[cf];
686: const PetscInt f = facesId[cntFaces];
687: const PetscInt *face = &faces[foff];
688: const PetscInt *fcone;
690: PetscCall(DMPlexInsertCone(idm, c, cf, f));
691: PetscCall(DMPlexGetCone(idm, f, &fcone));
692: if (fcone[0] < 0) PetscCall(DMPlexSetCone(idm, f, face));
693: {
694: const PetscInt *cone;
695: PetscInt coneSize, ornt;
697: PetscCall(DMPlexGetConeSize(idm, f, &coneSize));
698: PetscCall(DMPlexGetCone(idm, f, &cone));
699: PetscCheck(coneSize == faceSize, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid number of face vertices %" PetscInt_FMT " for face %" PetscInt_FMT " should be %" PetscInt_FMT, coneSize, f, faceSize);
700: /* Notice that we have to use vertices here because the lower dimensional faces have not been created yet */
701: PetscCall(DMPolytopeGetVertexOrientation(faceType, cone, face, &ornt));
702: PetscCall(DMPlexInsertConeOrientation(idm, c, cf, ornt));
703: }
704: cntFaces++;
705: }
706: PetscCall(DMPlexRestoreRawFaces_Internal(dm, ct, cone, &numFaces, &faceTypes, &faceSizes, &faces));
707: }
708: PetscCall(PetscFree(facesId));
709: PetscCall(DMPlexSymmetrize(idm));
710: PetscCall(DMPlexStratify(idm));
711: PetscFunctionReturn(PETSC_SUCCESS);
712: }
714: static PetscErrorCode SortRmineRremoteByRemote_Private(PetscSF sf, PetscInt *rmine1[], PetscInt *rremote1[])
715: {
716: PetscInt nleaves;
717: PetscInt nranks;
718: const PetscMPIInt *ranks = NULL;
719: const PetscInt *roffset = NULL, *rmine = NULL, *rremote = NULL;
720: PetscInt n, o, r;
722: PetscFunctionBegin;
723: PetscCall(PetscSFGetRootRanks(sf, &nranks, &ranks, &roffset, &rmine, &rremote));
724: nleaves = roffset[nranks];
725: PetscCall(PetscMalloc2(nleaves, rmine1, nleaves, rremote1));
726: for (r = 0; r < nranks; r++) {
727: /* simultaneously sort rank-wise portions of rmine & rremote by values in rremote
728: - to unify order with the other side */
729: o = roffset[r];
730: n = roffset[r + 1] - o;
731: PetscCall(PetscArraycpy(&(*rmine1)[o], &rmine[o], n));
732: PetscCall(PetscArraycpy(&(*rremote1)[o], &rremote[o], n));
733: PetscCall(PetscSortIntWithArray(n, &(*rremote1)[o], &(*rmine1)[o]));
734: }
735: PetscFunctionReturn(PETSC_SUCCESS);
736: }
738: PetscErrorCode DMPlexOrientInterface_Internal(DM dm)
739: {
740: PetscSF sf;
741: const PetscInt *locals;
742: const PetscSFNode *remotes;
743: const PetscMPIInt *ranks;
744: const PetscInt *roffset;
745: PetscInt *rmine1, *rremote1; /* rmine and rremote copies simultaneously sorted by rank and rremote */
746: PetscInt nroots, p, nleaves, nranks, r, maxConeSize = 0;
747: PetscInt(*roots)[4], (*leaves)[4], mainCone[4];
748: PetscMPIInt(*rootsRanks)[4], (*leavesRanks)[4];
749: MPI_Comm comm;
750: PetscMPIInt rank, size;
751: PetscInt debug = 0;
753: PetscFunctionBegin;
754: PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
755: PetscCallMPI(MPI_Comm_rank(comm, &rank));
756: PetscCallMPI(MPI_Comm_size(comm, &size));
757: PetscCall(DMGetPointSF(dm, &sf));
758: PetscCall(DMViewFromOptions(dm, NULL, "-before_orient_interface_dm_view"));
759: if (PetscDefined(USE_DEBUG)) PetscCall(DMPlexCheckPointSF(dm, sf, PETSC_FALSE));
760: PetscCall(PetscSFGetGraph(sf, &nroots, &nleaves, &locals, &remotes));
761: if (nroots < 0) PetscFunctionReturn(PETSC_SUCCESS);
762: PetscCall(PetscSFSetUp(sf));
763: PetscCall(SortRmineRremoteByRemote_Private(sf, &rmine1, &rremote1));
764: for (p = 0; p < nleaves; ++p) {
765: PetscInt coneSize;
766: PetscCall(DMPlexGetConeSize(dm, locals[p], &coneSize));
767: maxConeSize = PetscMax(maxConeSize, coneSize);
768: }
769: PetscCheck(maxConeSize <= 4, comm, PETSC_ERR_SUP, "This method does not support cones of size %" PetscInt_FMT, maxConeSize);
770: PetscCall(PetscMalloc4(nroots, &roots, nroots, &leaves, nroots, &rootsRanks, nroots, &leavesRanks));
771: for (p = 0; p < nroots; ++p) {
772: const PetscInt *cone;
773: PetscInt coneSize, c, ind0;
775: PetscCall(DMPlexGetConeSize(dm, p, &coneSize));
776: PetscCall(DMPlexGetCone(dm, p, &cone));
777: /* Ignore vertices */
778: if (coneSize < 2) {
779: for (c = 0; c < 4; c++) {
780: roots[p][c] = -1;
781: rootsRanks[p][c] = -1;
782: }
783: continue;
784: }
785: /* Translate all points to root numbering */
786: for (c = 0; c < PetscMin(coneSize, 4); c++) {
787: PetscCall(PetscFindInt(cone[c], nleaves, locals, &ind0));
788: if (ind0 < 0) {
789: roots[p][c] = cone[c];
790: rootsRanks[p][c] = rank;
791: } else {
792: roots[p][c] = remotes[ind0].index;
793: rootsRanks[p][c] = remotes[ind0].rank;
794: }
795: }
796: for (c = coneSize; c < 4; c++) {
797: roots[p][c] = -1;
798: rootsRanks[p][c] = -1;
799: }
800: }
801: for (p = 0; p < nroots; ++p) {
802: PetscInt c;
803: for (c = 0; c < 4; c++) {
804: leaves[p][c] = -2;
805: leavesRanks[p][c] = -2;
806: }
807: }
808: PetscCall(PetscSFBcastBegin(sf, MPIU_4INT, roots, leaves, MPI_REPLACE));
809: PetscCall(PetscSFBcastBegin(sf, MPI_4INT, rootsRanks, leavesRanks, MPI_REPLACE));
810: PetscCall(PetscSFBcastEnd(sf, MPIU_4INT, roots, leaves, MPI_REPLACE));
811: PetscCall(PetscSFBcastEnd(sf, MPI_4INT, rootsRanks, leavesRanks, MPI_REPLACE));
812: if (debug) {
813: PetscCall(PetscSynchronizedFlush(comm, NULL));
814: if (rank == 0) PetscCall(PetscSynchronizedPrintf(comm, "Referenced roots\n"));
815: }
816: PetscCall(PetscSFGetRootRanks(sf, &nranks, &ranks, &roffset, NULL, NULL));
817: for (p = 0; p < nroots; ++p) {
818: DMPolytopeType ct;
819: const PetscInt *cone;
820: PetscInt coneSize, c, ind0, o;
822: if (leaves[p][0] < 0) continue; /* Ignore vertices */
823: PetscCall(DMPlexGetCellType(dm, p, &ct));
824: PetscCall(DMPlexGetConeSize(dm, p, &coneSize));
825: PetscCall(DMPlexGetCone(dm, p, &cone));
826: if (debug) {
827: PetscCall(PetscSynchronizedPrintf(comm, "[%d] %4" PetscInt_FMT ": cone=[%4" PetscInt_FMT " %4" PetscInt_FMT " %4" PetscInt_FMT " %4" PetscInt_FMT "] roots=[(%d,%4" PetscInt_FMT ") (%d,%4" PetscInt_FMT ") (%d,%4" PetscInt_FMT ") (%d,%4" PetscInt_FMT ")] leaves=[(%d,%4" PetscInt_FMT ") (%d,%4" PetscInt_FMT ") (%d,%4" PetscInt_FMT ") (%d,%4" PetscInt_FMT ")]", rank, p, cone[0], cone[1], cone[2], cone[3], rootsRanks[p][0], roots[p][0], rootsRanks[p][1], roots[p][1], rootsRanks[p][2], roots[p][2], rootsRanks[p][3], roots[p][3], leavesRanks[p][0], leaves[p][0], leavesRanks[p][1], leaves[p][1], leavesRanks[p][2], leaves[p][2], leavesRanks[p][3], leaves[p][3]));
828: }
829: if (leavesRanks[p][0] != rootsRanks[p][0] || leaves[p][0] != roots[p][0] || leavesRanks[p][1] != rootsRanks[p][1] || leaves[p][1] != roots[p][1] || leavesRanks[p][2] != rootsRanks[p][2] || leaves[p][2] != roots[p][2] || leavesRanks[p][3] != rootsRanks[p][3] || leaves[p][3] != roots[p][3]) {
830: /* Translate these leaves to my cone points; mainCone means desired order p's cone points */
831: for (c = 0; c < PetscMin(coneSize, 4); ++c) {
832: PetscInt rS, rN;
834: if (leavesRanks[p][c] == rank) {
835: /* A local leaf is just taken as it is */
836: mainCone[c] = leaves[p][c];
837: continue;
838: }
839: /* Find index of rank leavesRanks[p][c] among remote ranks */
840: /* No need for PetscMPIIntCast because these integers were originally cast from PetscMPIInt. */
841: PetscCall(PetscFindMPIInt((PetscMPIInt)leavesRanks[p][c], nranks, ranks, &r));
842: PetscCheck(r >= 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Point %" PetscInt_FMT " cone[%" PetscInt_FMT "]=%" PetscInt_FMT " root (%d,%" PetscInt_FMT ") leaf (%d,%" PetscInt_FMT "): leaf rank not found among remote ranks", p, c, cone[c], rootsRanks[p][c], roots[p][c], leavesRanks[p][c], leaves[p][c]);
843: PetscCheck(ranks[r] >= 0 && ranks[r] < size, PETSC_COMM_SELF, PETSC_ERR_PLIB, "p=%" PetscInt_FMT " c=%" PetscInt_FMT " commsize=%d: ranks[%" PetscInt_FMT "] = %d makes no sense", p, c, size, r, ranks[r]);
844: /* Find point leaves[p][c] among remote points aimed at rank leavesRanks[p][c] */
845: rS = roffset[r];
846: rN = roffset[r + 1] - rS;
847: PetscCall(PetscFindInt(leaves[p][c], rN, &rremote1[rS], &ind0));
848: PetscCheck(ind0 >= 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Point %" PetscInt_FMT " cone[%" PetscInt_FMT "]=%" PetscInt_FMT " root (%d,%" PetscInt_FMT ") leave (%d,%" PetscInt_FMT "): corresponding remote point not found - it seems there is missing connection in point SF!", p, c, cone[c], rootsRanks[p][c], roots[p][c], leavesRanks[p][c], leaves[p][c]);
849: /* Get the corresponding local point */
850: mainCone[c] = rmine1[rS + ind0];
851: }
852: if (debug) PetscCall(PetscSynchronizedPrintf(comm, " mainCone=[%4" PetscInt_FMT " %4" PetscInt_FMT " %4" PetscInt_FMT " %4" PetscInt_FMT "]\n", mainCone[0], mainCone[1], mainCone[2], mainCone[3]));
853: /* Set the desired order of p's cone points and fix orientations accordingly */
854: PetscCall(DMPolytopeGetOrientation(ct, cone, mainCone, &o));
855: PetscCall(DMPlexOrientPoint(dm, p, o));
856: } else if (debug) PetscCall(PetscSynchronizedPrintf(comm, " ==\n"));
857: }
858: if (debug) {
859: PetscCall(PetscSynchronizedFlush(comm, NULL));
860: PetscCallMPI(MPI_Barrier(comm));
861: }
862: PetscCall(DMViewFromOptions(dm, NULL, "-after_orient_interface_dm_view"));
863: PetscCall(PetscFree4(roots, leaves, rootsRanks, leavesRanks));
864: PetscCall(PetscFree2(rmine1, rremote1));
865: PetscFunctionReturn(PETSC_SUCCESS);
866: }
868: static PetscErrorCode IntArrayViewFromOptions(MPI_Comm comm, const char opt[], const char name[], const char idxname[], const char valname[], PetscInt n, const PetscInt a[])
869: {
870: PetscInt idx;
871: PetscMPIInt rank;
872: PetscBool flg;
874: PetscFunctionBegin;
875: PetscCall(PetscOptionsHasName(NULL, NULL, opt, &flg));
876: if (!flg) PetscFunctionReturn(PETSC_SUCCESS);
877: PetscCallMPI(MPI_Comm_rank(comm, &rank));
878: PetscCall(PetscSynchronizedPrintf(comm, "[%d]%s:\n", rank, name));
879: for (idx = 0; idx < n; ++idx) PetscCall(PetscSynchronizedPrintf(comm, "[%d]%s %" PetscInt_FMT " %s %" PetscInt_FMT "\n", rank, idxname, idx, valname, a[idx]));
880: PetscCall(PetscSynchronizedFlush(comm, NULL));
881: PetscFunctionReturn(PETSC_SUCCESS);
882: }
884: static PetscErrorCode SFNodeArrayViewFromOptions(MPI_Comm comm, const char opt[], const char name[], const char idxname[], PetscInt n, const PetscSFNode a[])
885: {
886: PetscInt idx;
887: PetscMPIInt rank;
888: PetscBool flg;
890: PetscFunctionBegin;
891: PetscCall(PetscOptionsHasName(NULL, NULL, opt, &flg));
892: if (!flg) PetscFunctionReturn(PETSC_SUCCESS);
893: PetscCallMPI(MPI_Comm_rank(comm, &rank));
894: PetscCall(PetscSynchronizedPrintf(comm, "[%d]%s:\n", rank, name));
895: if (idxname) {
896: for (idx = 0; idx < n; ++idx) PetscCall(PetscSynchronizedPrintf(comm, "[%d]%s %" PetscInt_FMT " rank %" PetscInt_FMT " index %" PetscInt_FMT "\n", rank, idxname, idx, a[idx].rank, a[idx].index));
897: } else {
898: for (idx = 0; idx < n; ++idx) PetscCall(PetscSynchronizedPrintf(comm, "[%d]rank %" PetscInt_FMT " index %" PetscInt_FMT "\n", rank, a[idx].rank, a[idx].index));
899: }
900: PetscCall(PetscSynchronizedFlush(comm, NULL));
901: PetscFunctionReturn(PETSC_SUCCESS);
902: }
904: static PetscErrorCode DMPlexMapToLocalPoint(DM dm, PetscHMapIJ remotehash, PetscSFNode remotePoint, PetscInt *localPoint, PetscBool *mapFailed)
905: {
906: PetscSF sf;
907: const PetscInt *locals;
908: PetscMPIInt rank;
910: PetscFunctionBegin;
911: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank));
912: PetscCall(DMGetPointSF(dm, &sf));
913: PetscCall(PetscSFGetGraph(sf, NULL, NULL, &locals, NULL));
914: if (mapFailed) *mapFailed = PETSC_FALSE;
915: if (remotePoint.rank == rank) {
916: *localPoint = remotePoint.index;
917: } else {
918: PetscHashIJKey key;
919: PetscInt l;
921: key.i = remotePoint.index;
922: key.j = remotePoint.rank;
923: PetscCall(PetscHMapIJGet(remotehash, key, &l));
924: if (l >= 0) {
925: *localPoint = locals[l];
926: } else if (mapFailed) *mapFailed = PETSC_TRUE;
927: }
928: PetscFunctionReturn(PETSC_SUCCESS);
929: }
931: static PetscErrorCode DMPlexMapToGlobalPoint(DM dm, PetscInt localPoint, PetscSFNode *remotePoint, PetscBool *mapFailed)
932: {
933: PetscSF sf;
934: const PetscInt *locals, *rootdegree;
935: const PetscSFNode *remotes;
936: PetscInt Nl, l;
937: PetscMPIInt rank;
939: PetscFunctionBegin;
940: if (mapFailed) *mapFailed = PETSC_FALSE;
941: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank));
942: PetscCall(DMGetPointSF(dm, &sf));
943: PetscCall(PetscSFGetGraph(sf, NULL, &Nl, &locals, &remotes));
944: if (Nl < 0) goto owned;
945: PetscCall(PetscSFComputeDegreeBegin(sf, &rootdegree));
946: PetscCall(PetscSFComputeDegreeEnd(sf, &rootdegree));
947: if (rootdegree[localPoint]) goto owned;
948: PetscCall(PetscFindInt(localPoint, Nl, locals, &l));
949: if (l < 0) {
950: if (mapFailed) *mapFailed = PETSC_TRUE;
951: } else *remotePoint = remotes[l];
952: PetscFunctionReturn(PETSC_SUCCESS);
953: owned:
954: remotePoint->rank = rank;
955: remotePoint->index = localPoint;
956: PetscFunctionReturn(PETSC_SUCCESS);
957: }
959: static PetscErrorCode DMPlexPointIsShared(DM dm, PetscInt p, PetscBool *isShared)
960: {
961: PetscSF sf;
962: const PetscInt *locals, *rootdegree;
963: PetscInt Nl, idx;
965: PetscFunctionBegin;
966: *isShared = PETSC_FALSE;
967: PetscCall(DMGetPointSF(dm, &sf));
968: PetscCall(PetscSFGetGraph(sf, NULL, &Nl, &locals, NULL));
969: if (Nl < 0) PetscFunctionReturn(PETSC_SUCCESS);
970: PetscCall(PetscFindInt(p, Nl, locals, &idx));
971: if (idx >= 0) {
972: *isShared = PETSC_TRUE;
973: PetscFunctionReturn(PETSC_SUCCESS);
974: }
975: PetscCall(PetscSFComputeDegreeBegin(sf, &rootdegree));
976: PetscCall(PetscSFComputeDegreeEnd(sf, &rootdegree));
977: if (rootdegree[p] > 0) *isShared = PETSC_TRUE;
978: PetscFunctionReturn(PETSC_SUCCESS);
979: }
981: static PetscErrorCode DMPlexConeIsShared(DM dm, PetscInt p, PetscBool *isShared)
982: {
983: const PetscInt *cone;
984: PetscInt coneSize, c;
985: PetscBool cShared = PETSC_TRUE;
987: PetscFunctionBegin;
988: PetscCall(DMPlexGetConeSize(dm, p, &coneSize));
989: PetscCall(DMPlexGetCone(dm, p, &cone));
990: for (c = 0; c < coneSize; ++c) {
991: PetscBool pointShared;
993: PetscCall(DMPlexPointIsShared(dm, cone[c], &pointShared));
994: cShared = (PetscBool)(cShared && pointShared);
995: }
996: *isShared = coneSize ? cShared : PETSC_FALSE;
997: PetscFunctionReturn(PETSC_SUCCESS);
998: }
1000: static PetscErrorCode DMPlexGetConeMinimum(DM dm, PetscInt p, PetscSFNode *cpmin)
1001: {
1002: const PetscInt *cone;
1003: PetscInt coneSize, c;
1004: PetscSFNode cmin = {PETSC_MAX_INT, PETSC_MAX_INT}, missing = {-1, -1};
1006: PetscFunctionBegin;
1007: PetscCall(DMPlexGetConeSize(dm, p, &coneSize));
1008: PetscCall(DMPlexGetCone(dm, p, &cone));
1009: for (c = 0; c < coneSize; ++c) {
1010: PetscSFNode rcp;
1011: PetscBool mapFailed;
1013: PetscCall(DMPlexMapToGlobalPoint(dm, cone[c], &rcp, &mapFailed));
1014: if (mapFailed) {
1015: cmin = missing;
1016: } else {
1017: cmin = (rcp.rank < cmin.rank) || (rcp.rank == cmin.rank && rcp.index < cmin.index) ? rcp : cmin;
1018: }
1019: }
1020: *cpmin = coneSize ? cmin : missing;
1021: PetscFunctionReturn(PETSC_SUCCESS);
1022: }
1024: /*
1025: Each shared face has an entry in the candidates array:
1026: (-1, coneSize-1), {(global cone point)}
1027: where the set is missing the point p which we use as the key for the face
1028: */
1029: static PetscErrorCode DMPlexAddSharedFace_Private(DM dm, PetscSection candidateSection, PetscSFNode candidates[], PetscHMapIJ faceHash, PetscInt p, PetscBool debug)
1030: {
1031: MPI_Comm comm;
1032: const PetscInt *support;
1033: PetscInt supportSize, s, off = 0, idx = 0, overlap, cellHeight, height;
1034: PetscMPIInt rank;
1036: PetscFunctionBegin;
1037: PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
1038: PetscCallMPI(MPI_Comm_rank(comm, &rank));
1039: PetscCall(DMPlexGetOverlap(dm, &overlap));
1040: PetscCall(DMPlexGetVTKCellHeight(dm, &cellHeight));
1041: PetscCall(DMPlexGetPointHeight(dm, p, &height));
1042: if (!overlap && height <= cellHeight + 1) {
1043: /* cells can't be shared for non-overlapping meshes */
1044: if (debug) PetscCall(PetscSynchronizedPrintf(comm, "[%d] Skipping face %" PetscInt_FMT " to avoid adding cell to hashmap since this is nonoverlapping mesh\n", rank, p));
1045: PetscFunctionReturn(PETSC_SUCCESS);
1046: }
1047: PetscCall(DMPlexGetSupportSize(dm, p, &supportSize));
1048: PetscCall(DMPlexGetSupport(dm, p, &support));
1049: if (candidates) PetscCall(PetscSectionGetOffset(candidateSection, p, &off));
1050: for (s = 0; s < supportSize; ++s) {
1051: const PetscInt face = support[s];
1052: const PetscInt *cone;
1053: PetscSFNode cpmin = {-1, -1}, rp = {-1, -1};
1054: PetscInt coneSize, c, f;
1055: PetscBool isShared = PETSC_FALSE;
1056: PetscHashIJKey key;
1058: /* Only add point once */
1059: if (debug) PetscCall(PetscSynchronizedPrintf(comm, "[%d] Support face %" PetscInt_FMT "\n", rank, face));
1060: key.i = p;
1061: key.j = face;
1062: PetscCall(PetscHMapIJGet(faceHash, key, &f));
1063: if (f >= 0) continue;
1064: PetscCall(DMPlexConeIsShared(dm, face, &isShared));
1065: PetscCall(DMPlexGetConeMinimum(dm, face, &cpmin));
1066: PetscCall(DMPlexMapToGlobalPoint(dm, p, &rp, NULL));
1067: if (debug) {
1068: PetscCall(PetscSynchronizedPrintf(comm, "[%d] Face point %" PetscInt_FMT " is shared: %d\n", rank, face, (int)isShared));
1069: PetscCall(PetscSynchronizedPrintf(comm, "[%d] Global point (%" PetscInt_FMT ", %" PetscInt_FMT ") Min Cone Point (%" PetscInt_FMT ", %" PetscInt_FMT ")\n", rank, rp.rank, rp.index, cpmin.rank, cpmin.index));
1070: }
1071: if (isShared && (rp.rank == cpmin.rank && rp.index == cpmin.index)) {
1072: PetscCall(PetscHMapIJSet(faceHash, key, p));
1073: if (candidates) {
1074: if (debug) PetscCall(PetscSynchronizedPrintf(comm, "[%d] Adding shared face %" PetscInt_FMT " at idx %" PetscInt_FMT "\n[%d] ", rank, face, idx, rank));
1075: PetscCall(DMPlexGetConeSize(dm, face, &coneSize));
1076: PetscCall(DMPlexGetCone(dm, face, &cone));
1077: candidates[off + idx].rank = -1;
1078: candidates[off + idx++].index = coneSize - 1;
1079: candidates[off + idx].rank = rank;
1080: candidates[off + idx++].index = face;
1081: for (c = 0; c < coneSize; ++c) {
1082: const PetscInt cp = cone[c];
1084: if (cp == p) continue;
1085: PetscCall(DMPlexMapToGlobalPoint(dm, cp, &candidates[off + idx], NULL));
1086: if (debug) PetscCall(PetscSynchronizedPrintf(comm, " (%" PetscInt_FMT ",%" PetscInt_FMT ")", candidates[off + idx].rank, candidates[off + idx].index));
1087: ++idx;
1088: }
1089: if (debug) PetscCall(PetscSynchronizedPrintf(comm, "\n"));
1090: } else {
1091: /* Add cone size to section */
1092: if (debug) PetscCall(PetscSynchronizedPrintf(comm, "[%d] Scheduling shared face %" PetscInt_FMT "\n", rank, face));
1093: PetscCall(DMPlexGetConeSize(dm, face, &coneSize));
1094: PetscCall(PetscHMapIJSet(faceHash, key, p));
1095: PetscCall(PetscSectionAddDof(candidateSection, p, coneSize + 1));
1096: }
1097: }
1098: }
1099: PetscFunctionReturn(PETSC_SUCCESS);
1100: }
1102: /*@
1103: DMPlexInterpolatePointSF - Insert interpolated points in the overlap into the `PointSF` in parallel, following local interpolation
1105: Collective
1107: Input Parameters:
1108: + dm - The interpolated `DMPLEX`
1109: - pointSF - The initial `PetscSF` without interpolated points
1111: Level: developer
1113: Note:
1114: 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`
1116: .seealso: `DMPLEX`, `DMPlexInterpolate()`, `DMPlexUninterpolate()`
1117: @*/
1118: PetscErrorCode DMPlexInterpolatePointSF(DM dm, PetscSF pointSF)
1119: {
1120: MPI_Comm comm;
1121: PetscHMapIJ remoteHash;
1122: PetscHMapI claimshash;
1123: PetscSection candidateSection, candidateRemoteSection, claimSection;
1124: PetscSFNode *candidates, *candidatesRemote, *claims;
1125: const PetscInt *localPoints, *rootdegree;
1126: const PetscSFNode *remotePoints;
1127: PetscInt ov, Nr, r, Nl, l;
1128: PetscInt candidatesSize, candidatesRemoteSize, claimsSize;
1129: PetscBool flg, debug = PETSC_FALSE;
1130: PetscMPIInt rank;
1132: PetscFunctionBegin;
1135: PetscCall(DMPlexIsDistributed(dm, &flg));
1136: if (!flg) PetscFunctionReturn(PETSC_SUCCESS);
1137: /* Set initial SF so that lower level queries work */
1138: PetscCall(DMSetPointSF(dm, pointSF));
1139: PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
1140: PetscCallMPI(MPI_Comm_rank(comm, &rank));
1141: PetscCall(DMPlexGetOverlap(dm, &ov));
1142: PetscCheck(!ov, comm, PETSC_ERR_SUP, "Interpolation of overlapped DMPlex not implemented yet");
1143: PetscCall(PetscOptionsHasName(NULL, ((PetscObject)dm)->prefix, "-dmplex_interp_debug", &debug));
1144: PetscCall(PetscObjectViewFromOptions((PetscObject)dm, NULL, "-dm_interp_pre_view"));
1145: PetscCall(PetscObjectViewFromOptions((PetscObject)pointSF, NULL, "-petscsf_interp_pre_view"));
1146: PetscCall(PetscLogEventBegin(DMPLEX_InterpolateSF, dm, 0, 0, 0));
1147: /* Step 0: Precalculations */
1148: PetscCall(PetscSFGetGraph(pointSF, &Nr, &Nl, &localPoints, &remotePoints));
1149: PetscCheck(Nr >= 0, comm, PETSC_ERR_ARG_WRONGSTATE, "This DMPlex is distributed but input PointSF has no graph set");
1150: PetscCall(PetscHMapIJCreate(&remoteHash));
1151: for (l = 0; l < Nl; ++l) {
1152: PetscHashIJKey key;
1153: key.i = remotePoints[l].index;
1154: key.j = remotePoints[l].rank;
1155: PetscCall(PetscHMapIJSet(remoteHash, key, l));
1156: }
1157: /* Compute root degree to identify shared points */
1158: PetscCall(PetscSFComputeDegreeBegin(pointSF, &rootdegree));
1159: PetscCall(PetscSFComputeDegreeEnd(pointSF, &rootdegree));
1160: PetscCall(IntArrayViewFromOptions(comm, "-interp_root_degree_view", "Root degree", "point", "degree", Nr, rootdegree));
1161: /*
1162: 1) Loop over each leaf point $p$ at depth $d$ in the SF
1163: \item Get set $F(p)$ of faces $f$ in the support of $p$ for which
1164: \begin{itemize}
1165: \item all cone points of $f$ are shared
1166: \item $p$ is the cone point with smallest canonical number
1167: \end{itemize}
1168: \item Send $F(p)$ and the cone of each face to the active root point $r(p)$
1169: \item At the root, if at least two faces with a given cone are present, including a local face, mark the face as shared \label{alg:rootStep} and choose the root face
1170: \item Send the root face from the root back to all leaf process
1171: \item Leaf processes add the shared face to the SF
1172: */
1173: /* Step 1: Construct section+SFNode array
1174: The section has entries for all shared faces for which we have a leaf point in the cone
1175: The array holds candidate shared faces, each face is referred to by the leaf point */
1176: PetscCall(PetscSectionCreate(comm, &candidateSection));
1177: PetscCall(PetscSectionSetChart(candidateSection, 0, Nr));
1178: {
1179: PetscHMapIJ faceHash;
1181: PetscCall(PetscHMapIJCreate(&faceHash));
1182: for (l = 0; l < Nl; ++l) {
1183: const PetscInt p = localPoints[l];
1185: if (debug) PetscCall(PetscSynchronizedPrintf(comm, "[%d] First pass leaf point %" PetscInt_FMT "\n", rank, p));
1186: PetscCall(DMPlexAddSharedFace_Private(dm, candidateSection, NULL, faceHash, p, debug));
1187: }
1188: PetscCall(PetscHMapIJClear(faceHash));
1189: PetscCall(PetscSectionSetUp(candidateSection));
1190: PetscCall(PetscSectionGetStorageSize(candidateSection, &candidatesSize));
1191: PetscCall(PetscMalloc1(candidatesSize, &candidates));
1192: for (l = 0; l < Nl; ++l) {
1193: const PetscInt p = localPoints[l];
1195: if (debug) PetscCall(PetscSynchronizedPrintf(comm, "[%d] Second pass leaf point %" PetscInt_FMT "\n", rank, p));
1196: PetscCall(DMPlexAddSharedFace_Private(dm, candidateSection, candidates, faceHash, p, debug));
1197: }
1198: PetscCall(PetscHMapIJDestroy(&faceHash));
1199: if (debug) PetscCall(PetscSynchronizedFlush(comm, NULL));
1200: }
1201: PetscCall(PetscObjectSetName((PetscObject)candidateSection, "Candidate Section"));
1202: PetscCall(PetscObjectViewFromOptions((PetscObject)candidateSection, NULL, "-petscsection_interp_candidate_view"));
1203: PetscCall(SFNodeArrayViewFromOptions(comm, "-petscsection_interp_candidate_view", "Candidates", NULL, candidatesSize, candidates));
1204: /* Step 2: Gather candidate section / array pair into the root partition via inverse(multi(pointSF)). */
1205: /* Note that this section is indexed by offsets into leaves, not by point number */
1206: {
1207: PetscSF sfMulti, sfInverse, sfCandidates;
1208: PetscInt *remoteOffsets;
1210: PetscCall(PetscSFGetMultiSF(pointSF, &sfMulti));
1211: PetscCall(PetscSFCreateInverseSF(sfMulti, &sfInverse));
1212: PetscCall(PetscSectionCreate(comm, &candidateRemoteSection));
1213: PetscCall(PetscSFDistributeSection(sfInverse, candidateSection, &remoteOffsets, candidateRemoteSection));
1214: PetscCall(PetscSFCreateSectionSF(sfInverse, candidateSection, remoteOffsets, candidateRemoteSection, &sfCandidates));
1215: PetscCall(PetscSectionGetStorageSize(candidateRemoteSection, &candidatesRemoteSize));
1216: PetscCall(PetscMalloc1(candidatesRemoteSize, &candidatesRemote));
1217: PetscCall(PetscSFBcastBegin(sfCandidates, MPIU_2INT, candidates, candidatesRemote, MPI_REPLACE));
1218: PetscCall(PetscSFBcastEnd(sfCandidates, MPIU_2INT, candidates, candidatesRemote, MPI_REPLACE));
1219: PetscCall(PetscSFDestroy(&sfInverse));
1220: PetscCall(PetscSFDestroy(&sfCandidates));
1221: PetscCall(PetscFree(remoteOffsets));
1223: PetscCall(PetscObjectSetName((PetscObject)candidateRemoteSection, "Remote Candidate Section"));
1224: PetscCall(PetscObjectViewFromOptions((PetscObject)candidateRemoteSection, NULL, "-petscsection_interp_candidate_remote_view"));
1225: PetscCall(SFNodeArrayViewFromOptions(comm, "-petscsection_interp_candidate_remote_view", "Remote Candidates", NULL, candidatesRemoteSize, candidatesRemote));
1226: }
1227: /* Step 3: At the root, if at least two faces with a given cone are present, including a local face, mark the face as shared and choose the root face */
1228: {
1229: PetscHMapIJKLRemote faceTable;
1230: PetscInt idx, idx2;
1232: PetscCall(PetscHMapIJKLRemoteCreate(&faceTable));
1233: /* There is a section point for every leaf attached to a given root point */
1234: for (r = 0, idx = 0, idx2 = 0; r < Nr; ++r) {
1235: PetscInt deg;
1237: for (deg = 0; deg < rootdegree[r]; ++deg, ++idx) {
1238: PetscInt offset, dof, d;
1240: PetscCall(PetscSectionGetDof(candidateRemoteSection, idx, &dof));
1241: PetscCall(PetscSectionGetOffset(candidateRemoteSection, idx, &offset));
1242: /* dof may include many faces from the remote process */
1243: for (d = 0; d < dof; ++d) {
1244: const PetscInt hidx = offset + d;
1245: const PetscInt Np = candidatesRemote[hidx].index + 1;
1246: const PetscSFNode rface = candidatesRemote[hidx + 1];
1247: const PetscSFNode *fcone = &candidatesRemote[hidx + 2];
1248: PetscSFNode fcp0;
1249: const PetscSFNode pmax = {PETSC_MAX_INT, PETSC_MAX_INT};
1250: const PetscInt *join = NULL;
1251: PetscHMapIJKLRemoteKey key;
1252: PetscHashIter iter;
1253: PetscBool missing, mapToLocalPointFailed = PETSC_FALSE;
1254: PetscInt points[1024], p, joinSize;
1256: if (debug)
1257: PetscCall(PetscSynchronizedPrintf(PetscObjectComm((PetscObject)dm), "[%d] Checking face (%" PetscInt_FMT ", %" PetscInt_FMT ") at (%" PetscInt_FMT ", %" PetscInt_FMT ", %" PetscInt_FMT ") with cone size %" PetscInt_FMT "\n", rank, rface.rank,
1258: rface.index, r, idx, d, Np));
1259: PetscCheck(Np <= 4, PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot handle face (%" PetscInt_FMT ", %" PetscInt_FMT ") at (%" PetscInt_FMT ", %" PetscInt_FMT ", %" PetscInt_FMT ") with %" PetscInt_FMT " cone points", rface.rank, rface.index, r, idx, d, Np);
1260: fcp0.rank = rank;
1261: fcp0.index = r;
1262: d += Np;
1263: /* Put remote face in hash table */
1264: key.i = fcp0;
1265: key.j = fcone[0];
1266: key.k = Np > 2 ? fcone[1] : pmax;
1267: key.l = Np > 3 ? fcone[2] : pmax;
1268: PetscCall(PetscSortSFNode(Np, (PetscSFNode *)&key));
1269: PetscCall(PetscHMapIJKLRemotePut(faceTable, key, &iter, &missing));
1270: if (missing) {
1271: if (debug) PetscCall(PetscSynchronizedPrintf(PetscObjectComm((PetscObject)dm), "[%d] Setting remote face (%" PetscInt_FMT ", %" PetscInt_FMT ")\n", rank, rface.index, rface.rank));
1272: PetscCall(PetscHMapIJKLRemoteIterSet(faceTable, iter, rface));
1273: } else {
1274: PetscSFNode oface;
1276: PetscCall(PetscHMapIJKLRemoteIterGet(faceTable, iter, &oface));
1277: if ((rface.rank < oface.rank) || (rface.rank == oface.rank && rface.index < oface.index)) {
1278: if (debug) PetscCall(PetscSynchronizedPrintf(PetscObjectComm((PetscObject)dm), "[%d] Replacing with remote face (%" PetscInt_FMT ", %" PetscInt_FMT ")\n", rank, rface.index, rface.rank));
1279: PetscCall(PetscHMapIJKLRemoteIterSet(faceTable, iter, rface));
1280: }
1281: }
1282: /* Check for local face */
1283: points[0] = r;
1284: for (p = 1; p < Np; ++p) {
1285: PetscCall(DMPlexMapToLocalPoint(dm, remoteHash, fcone[p - 1], &points[p], &mapToLocalPointFailed));
1286: if (mapToLocalPointFailed) break; /* We got a point not in our overlap */
1287: if (debug) PetscCall(PetscSynchronizedPrintf(PetscObjectComm((PetscObject)dm), "[%d] Checking local candidate %" PetscInt_FMT "\n", rank, points[p]));
1288: }
1289: if (mapToLocalPointFailed) continue;
1290: PetscCall(DMPlexGetJoin(dm, Np, points, &joinSize, &join));
1291: if (joinSize == 1) {
1292: PetscSFNode lface;
1293: PetscSFNode oface;
1295: /* Always replace with local face */
1296: lface.rank = rank;
1297: lface.index = join[0];
1298: PetscCall(PetscHMapIJKLRemoteIterGet(faceTable, iter, &oface));
1299: if (debug)
1300: PetscCall(PetscSynchronizedPrintf(PetscObjectComm((PetscObject)dm), "[%d] Replacing (%" PetscInt_FMT ", %" PetscInt_FMT ") with local face (%" PetscInt_FMT ", %" PetscInt_FMT ")\n", rank, oface.index, oface.rank, lface.index, lface.rank));
1301: PetscCall(PetscHMapIJKLRemoteIterSet(faceTable, iter, lface));
1302: }
1303: PetscCall(DMPlexRestoreJoin(dm, Np, points, &joinSize, &join));
1304: }
1305: }
1306: /* Put back faces for this root */
1307: for (deg = 0; deg < rootdegree[r]; ++deg, ++idx2) {
1308: PetscInt offset, dof, d;
1310: PetscCall(PetscSectionGetDof(candidateRemoteSection, idx2, &dof));
1311: PetscCall(PetscSectionGetOffset(candidateRemoteSection, idx2, &offset));
1312: /* dof may include many faces from the remote process */
1313: for (d = 0; d < dof; ++d) {
1314: const PetscInt hidx = offset + d;
1315: const PetscInt Np = candidatesRemote[hidx].index + 1;
1316: const PetscSFNode *fcone = &candidatesRemote[hidx + 2];
1317: PetscSFNode fcp0;
1318: const PetscSFNode pmax = {PETSC_MAX_INT, PETSC_MAX_INT};
1319: PetscHMapIJKLRemoteKey key;
1320: PetscHashIter iter;
1321: PetscBool missing;
1323: if (debug) PetscCall(PetscSynchronizedPrintf(PetscObjectComm((PetscObject)dm), "[%d] Entering face at (%" PetscInt_FMT ", %" PetscInt_FMT ")\n", rank, r, idx));
1324: PetscCheck(Np <= 4, PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot handle faces with %" PetscInt_FMT " cone points", Np);
1325: fcp0.rank = rank;
1326: fcp0.index = r;
1327: d += Np;
1328: /* Find remote face in hash table */
1329: key.i = fcp0;
1330: key.j = fcone[0];
1331: key.k = Np > 2 ? fcone[1] : pmax;
1332: key.l = Np > 3 ? fcone[2] : pmax;
1333: PetscCall(PetscSortSFNode(Np, (PetscSFNode *)&key));
1334: if (debug)
1335: PetscCall(PetscSynchronizedPrintf(PetscObjectComm((PetscObject)dm), "[%d] key (%" PetscInt_FMT ", %" PetscInt_FMT ") (%" PetscInt_FMT ", %" PetscInt_FMT ") (%" PetscInt_FMT ", %" PetscInt_FMT ") (%" PetscInt_FMT ", %" PetscInt_FMT ")\n", rank,
1336: key.i.rank, key.i.index, key.j.rank, key.j.index, key.k.rank, key.k.index, key.l.rank, key.l.index));
1337: PetscCall(PetscHMapIJKLRemotePut(faceTable, key, &iter, &missing));
1338: PetscCheck(!missing, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Root %" PetscInt_FMT " Idx %" PetscInt_FMT " ought to have an associated face", r, idx2);
1339: PetscCall(PetscHMapIJKLRemoteIterGet(faceTable, iter, &candidatesRemote[hidx]));
1340: }
1341: }
1342: }
1343: if (debug) PetscCall(PetscSynchronizedFlush(PetscObjectComm((PetscObject)dm), NULL));
1344: PetscCall(PetscHMapIJKLRemoteDestroy(&faceTable));
1345: }
1346: /* Step 4: Push back owned faces */
1347: {
1348: PetscSF sfMulti, sfClaims, sfPointNew;
1349: PetscSFNode *remotePointsNew;
1350: PetscInt *remoteOffsets, *localPointsNew;
1351: PetscInt pStart, pEnd, r, NlNew, p;
1353: /* 4) Push claims back to receiver via the MultiSF and derive new pointSF mapping on receiver */
1354: PetscCall(PetscSFGetMultiSF(pointSF, &sfMulti));
1355: PetscCall(PetscSectionCreate(comm, &claimSection));
1356: PetscCall(PetscSFDistributeSection(sfMulti, candidateRemoteSection, &remoteOffsets, claimSection));
1357: PetscCall(PetscSFCreateSectionSF(sfMulti, candidateRemoteSection, remoteOffsets, claimSection, &sfClaims));
1358: PetscCall(PetscSectionGetStorageSize(claimSection, &claimsSize));
1359: PetscCall(PetscMalloc1(claimsSize, &claims));
1360: for (p = 0; p < claimsSize; ++p) claims[p].rank = -1;
1361: PetscCall(PetscSFBcastBegin(sfClaims, MPIU_2INT, candidatesRemote, claims, MPI_REPLACE));
1362: PetscCall(PetscSFBcastEnd(sfClaims, MPIU_2INT, candidatesRemote, claims, MPI_REPLACE));
1363: PetscCall(PetscSFDestroy(&sfClaims));
1364: PetscCall(PetscFree(remoteOffsets));
1365: PetscCall(PetscObjectSetName((PetscObject)claimSection, "Claim Section"));
1366: PetscCall(PetscObjectViewFromOptions((PetscObject)claimSection, NULL, "-petscsection_interp_claim_view"));
1367: PetscCall(SFNodeArrayViewFromOptions(comm, "-petscsection_interp_claim_view", "Claims", NULL, claimsSize, claims));
1368: /* Step 5) Walk the original section of local supports and add an SF entry for each updated item */
1369: /* TODO I should not have to do a join here since I already put the face and its cone in the candidate section */
1370: PetscCall(PetscHMapICreate(&claimshash));
1371: for (r = 0; r < Nr; ++r) {
1372: PetscInt dof, off, d;
1374: if (debug) PetscCall(PetscSynchronizedPrintf(comm, "[%d] Checking root for claims %" PetscInt_FMT "\n", rank, r));
1375: PetscCall(PetscSectionGetDof(candidateSection, r, &dof));
1376: PetscCall(PetscSectionGetOffset(candidateSection, r, &off));
1377: for (d = 0; d < dof;) {
1378: if (claims[off + d].rank >= 0) {
1379: const PetscInt faceInd = off + d;
1380: const PetscInt Np = candidates[off + d].index;
1381: const PetscInt *join = NULL;
1382: PetscInt joinSize, points[1024], c;
1384: if (debug) PetscCall(PetscSynchronizedPrintf(comm, "[%d] Found claim for remote point (%" PetscInt_FMT ", %" PetscInt_FMT ")\n", rank, claims[faceInd].rank, claims[faceInd].index));
1385: points[0] = r;
1386: if (debug) PetscCall(PetscSynchronizedPrintf(comm, "[%d] point %" PetscInt_FMT "\n", rank, points[0]));
1387: for (c = 0, d += 2; c < Np; ++c, ++d) {
1388: PetscCall(DMPlexMapToLocalPoint(dm, remoteHash, candidates[off + d], &points[c + 1], NULL));
1389: if (debug) PetscCall(PetscSynchronizedPrintf(comm, "[%d] point %" PetscInt_FMT "\n", rank, points[c + 1]));
1390: }
1391: PetscCall(DMPlexGetJoin(dm, Np + 1, points, &joinSize, &join));
1392: if (joinSize == 1) {
1393: if (claims[faceInd].rank == rank) {
1394: if (debug) PetscCall(PetscSynchronizedPrintf(comm, "[%d] Ignoring local face %" PetscInt_FMT " for non-remote partner\n", rank, join[0]));
1395: } else {
1396: if (debug) PetscCall(PetscSynchronizedPrintf(comm, "[%d] Found local face %" PetscInt_FMT "\n", rank, join[0]));
1397: PetscCall(PetscHMapISet(claimshash, join[0], faceInd));
1398: }
1399: } else {
1400: if (debug) PetscCall(PetscSynchronizedPrintf(comm, "[%d] Failed to find face\n", rank));
1401: }
1402: PetscCall(DMPlexRestoreJoin(dm, Np + 1, points, &joinSize, &join));
1403: } else {
1404: if (debug) PetscCall(PetscSynchronizedPrintf(comm, "[%d] No claim for point %" PetscInt_FMT "\n", rank, r));
1405: d += claims[off + d].index + 1;
1406: }
1407: }
1408: }
1409: if (debug) PetscCall(PetscSynchronizedFlush(comm, NULL));
1410: /* Step 6) Create new pointSF from hashed claims */
1411: PetscCall(PetscHMapIGetSize(claimshash, &NlNew));
1412: PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
1413: PetscCall(PetscMalloc1(Nl + NlNew, &localPointsNew));
1414: PetscCall(PetscMalloc1(Nl + NlNew, &remotePointsNew));
1415: for (l = 0; l < Nl; ++l) {
1416: localPointsNew[l] = localPoints[l];
1417: remotePointsNew[l].index = remotePoints[l].index;
1418: remotePointsNew[l].rank = remotePoints[l].rank;
1419: }
1420: p = Nl;
1421: PetscCall(PetscHMapIGetKeys(claimshash, &p, localPointsNew));
1422: /* We sort new points, and assume they are numbered after all existing points */
1423: PetscCall(PetscSortInt(NlNew, &localPointsNew[Nl]));
1424: for (p = Nl; p < Nl + NlNew; ++p) {
1425: PetscInt off;
1426: PetscCall(PetscHMapIGet(claimshash, localPointsNew[p], &off));
1427: PetscCheck(claims[off].rank >= 0 && claims[off].index >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid claim for local point %" PetscInt_FMT ", (%" PetscInt_FMT ", %" PetscInt_FMT ")", localPointsNew[p], claims[off].rank, claims[off].index);
1428: remotePointsNew[p] = claims[off];
1429: }
1430: PetscCall(PetscSFCreate(comm, &sfPointNew));
1431: PetscCall(PetscSFSetGraph(sfPointNew, pEnd - pStart, Nl + NlNew, localPointsNew, PETSC_OWN_POINTER, remotePointsNew, PETSC_OWN_POINTER));
1432: PetscCall(PetscSFSetUp(sfPointNew));
1433: PetscCall(DMSetPointSF(dm, sfPointNew));
1434: PetscCall(PetscObjectViewFromOptions((PetscObject)sfPointNew, NULL, "-petscsf_interp_view"));
1435: if (PetscDefined(USE_DEBUG)) PetscCall(DMPlexCheckPointSF(dm, sfPointNew, PETSC_FALSE));
1436: PetscCall(PetscSFDestroy(&sfPointNew));
1437: PetscCall(PetscHMapIDestroy(&claimshash));
1438: }
1439: PetscCall(PetscHMapIJDestroy(&remoteHash));
1440: PetscCall(PetscSectionDestroy(&candidateSection));
1441: PetscCall(PetscSectionDestroy(&candidateRemoteSection));
1442: PetscCall(PetscSectionDestroy(&claimSection));
1443: PetscCall(PetscFree(candidates));
1444: PetscCall(PetscFree(candidatesRemote));
1445: PetscCall(PetscFree(claims));
1446: PetscCall(PetscLogEventEnd(DMPLEX_InterpolateSF, dm, 0, 0, 0));
1447: PetscFunctionReturn(PETSC_SUCCESS);
1448: }
1450: /*@
1451: DMPlexInterpolate - Take in a cell-vertex mesh and return one with all intermediate faces, edges, etc.
1453: Collective
1455: Input Parameter:
1456: . dm - The `DMPLEX` object with only cells and vertices
1458: Output Parameter:
1459: . dmInt - The complete `DMPLEX` object
1461: Level: intermediate
1463: Note:
1464: Labels and coordinates are copied.
1466: Developer Notes:
1467: It sets plex->interpolated = `DMPLEX_INTERPOLATED_FULL`.
1469: .seealso: `DMPLEX`, `DMPlexUninterpolate()`, `DMPlexCreateFromCellListPetsc()`, `DMPlexCopyCoordinates()`
1470: @*/
1471: PetscErrorCode DMPlexInterpolate(DM dm, DM *dmInt)
1472: {
1473: DMPlexInterpolatedFlag interpolated;
1474: DM idm, odm = dm;
1475: PetscSF sfPoint;
1476: PetscInt depth, dim, d;
1477: const char *name;
1478: PetscBool flg = PETSC_TRUE;
1480: PetscFunctionBegin;
1482: PetscAssertPointer(dmInt, 2);
1483: PetscCall(PetscLogEventBegin(DMPLEX_Interpolate, dm, 0, 0, 0));
1484: PetscCall(DMPlexGetDepth(dm, &depth));
1485: PetscCall(DMGetDimension(dm, &dim));
1486: PetscCall(DMPlexIsInterpolated(dm, &interpolated));
1487: PetscCheck(interpolated != DMPLEX_INTERPOLATED_PARTIAL, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Not for partially interpolated meshes");
1488: if (interpolated == DMPLEX_INTERPOLATED_FULL) {
1489: PetscCall(PetscObjectReference((PetscObject)dm));
1490: idm = dm;
1491: } else {
1492: for (d = 1; d < dim; ++d) {
1493: /* Create interpolated mesh */
1494: PetscCall(DMCreate(PetscObjectComm((PetscObject)dm), &idm));
1495: PetscCall(DMSetType(idm, DMPLEX));
1496: PetscCall(DMSetDimension(idm, dim));
1497: if (depth > 0) {
1498: PetscCall(DMPlexInterpolateFaces_Internal(odm, 1, idm));
1499: PetscCall(DMGetPointSF(odm, &sfPoint));
1500: if (PetscDefined(USE_DEBUG)) PetscCall(DMPlexCheckPointSF(odm, sfPoint, PETSC_FALSE));
1501: {
1502: /* TODO: We need to systematically fix cases of distributed Plexes with no graph set */
1503: PetscInt nroots;
1504: PetscCall(PetscSFGetGraph(sfPoint, &nroots, NULL, NULL, NULL));
1505: if (nroots >= 0) PetscCall(DMPlexInterpolatePointSF(idm, sfPoint));
1506: }
1507: }
1508: if (odm != dm) PetscCall(DMDestroy(&odm));
1509: odm = idm;
1510: }
1511: PetscCall(PetscObjectGetName((PetscObject)dm, &name));
1512: PetscCall(PetscObjectSetName((PetscObject)idm, name));
1513: PetscCall(DMPlexCopyCoordinates(dm, idm));
1514: PetscCall(DMCopyLabels(dm, idm, PETSC_COPY_VALUES, PETSC_FALSE, DM_COPY_LABELS_FAIL));
1515: PetscCall(PetscOptionsGetBool(((PetscObject)dm)->options, ((PetscObject)dm)->prefix, "-dm_plex_interpolate_orient_interfaces", &flg, NULL));
1516: if (flg) PetscCall(DMPlexOrientInterface_Internal(idm));
1517: }
1518: /* This function makes the mesh fully interpolated on all ranks */
1519: {
1520: DM_Plex *plex = (DM_Plex *)idm->data;
1521: plex->interpolated = plex->interpolatedCollective = DMPLEX_INTERPOLATED_FULL;
1522: }
1523: PetscCall(DMPlexCopy_Internal(dm, PETSC_TRUE, PETSC_TRUE, idm));
1524: *dmInt = idm;
1525: PetscCall(PetscLogEventEnd(DMPLEX_Interpolate, dm, 0, 0, 0));
1526: PetscFunctionReturn(PETSC_SUCCESS);
1527: }
1529: /*@
1530: DMPlexCopyCoordinates - Copy coordinates from one mesh to another with the same vertices
1532: Collective
1534: Input Parameter:
1535: . dmA - The `DMPLEX` object with initial coordinates
1537: Output Parameter:
1538: . dmB - The `DMPLEX` object with copied coordinates
1540: Level: intermediate
1542: Note:
1543: This is typically used when adding pieces other than vertices to a mesh
1545: .seealso: `DMPLEX`, `DMCopyLabels()`, `DMGetCoordinates()`, `DMGetCoordinatesLocal()`, `DMGetCoordinateDM()`, `DMGetCoordinateSection()`
1546: @*/
1547: PetscErrorCode DMPlexCopyCoordinates(DM dmA, DM dmB)
1548: {
1549: Vec coordinatesA, coordinatesB;
1550: VecType vtype;
1551: PetscSection coordSectionA, coordSectionB;
1552: PetscScalar *coordsA, *coordsB;
1553: PetscInt spaceDim, Nf, vStartA, vStartB, vEndA, vEndB, coordSizeB, v, d;
1554: PetscInt cStartA, cEndA, cStartB, cEndB, cS, cE, cdim;
1555: PetscBool lc = PETSC_FALSE;
1557: PetscFunctionBegin;
1560: if (dmA == dmB) PetscFunctionReturn(PETSC_SUCCESS);
1561: PetscCall(DMGetCoordinateDim(dmA, &cdim));
1562: PetscCall(DMSetCoordinateDim(dmB, cdim));
1563: PetscCall(DMPlexGetDepthStratum(dmA, 0, &vStartA, &vEndA));
1564: PetscCall(DMPlexGetDepthStratum(dmB, 0, &vStartB, &vEndB));
1565: PetscCheck((vEndA - vStartA) == (vEndB - vStartB), PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "The number of vertices in first DM %" PetscInt_FMT " != %" PetscInt_FMT " in the second DM", vEndA - vStartA, vEndB - vStartB);
1566: /* Copy over discretization if it exists */
1567: {
1568: DM cdmA, cdmB;
1569: PetscDS dsA, dsB;
1570: PetscObject objA, objB;
1571: PetscClassId idA, idB;
1572: const PetscScalar *constants;
1573: PetscInt cdim, Nc;
1575: PetscCall(DMGetCoordinateDM(dmA, &cdmA));
1576: PetscCall(DMGetCoordinateDM(dmB, &cdmB));
1577: PetscCall(DMGetField(cdmA, 0, NULL, &objA));
1578: PetscCall(DMGetField(cdmB, 0, NULL, &objB));
1579: PetscCall(PetscObjectGetClassId(objA, &idA));
1580: PetscCall(PetscObjectGetClassId(objB, &idB));
1581: if ((idA == PETSCFE_CLASSID) && (idA != idB)) {
1582: PetscCall(DMSetField(cdmB, 0, NULL, objA));
1583: PetscCall(DMCreateDS(cdmB));
1584: PetscCall(DMGetDS(cdmA, &dsA));
1585: PetscCall(DMGetDS(cdmB, &dsB));
1586: PetscCall(PetscDSGetCoordinateDimension(dsA, &cdim));
1587: PetscCall(PetscDSSetCoordinateDimension(dsB, cdim));
1588: PetscCall(PetscDSGetConstants(dsA, &Nc, &constants));
1589: PetscCall(PetscDSSetConstants(dsB, Nc, (PetscScalar *)constants));
1590: }
1591: }
1592: PetscCall(DMPlexGetHeightStratum(dmA, 0, &cStartA, &cEndA));
1593: PetscCall(DMPlexGetHeightStratum(dmB, 0, &cStartB, &cEndB));
1594: PetscCall(DMGetCoordinateSection(dmA, &coordSectionA));
1595: PetscCall(DMGetCoordinateSection(dmB, &coordSectionB));
1596: if (coordSectionA == coordSectionB) PetscFunctionReturn(PETSC_SUCCESS);
1597: PetscCall(PetscSectionGetNumFields(coordSectionA, &Nf));
1598: if (!Nf) PetscFunctionReturn(PETSC_SUCCESS);
1599: PetscCheck(Nf <= 1, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "The number of coordinate fields must be 1, not %" PetscInt_FMT, Nf);
1600: if (!coordSectionB) {
1601: PetscInt dim;
1603: PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)coordSectionA), &coordSectionB));
1604: PetscCall(DMGetCoordinateDim(dmA, &dim));
1605: PetscCall(DMSetCoordinateSection(dmB, dim, coordSectionB));
1606: PetscCall(PetscObjectDereference((PetscObject)coordSectionB));
1607: }
1608: PetscCall(PetscSectionSetNumFields(coordSectionB, 1));
1609: PetscCall(PetscSectionGetFieldComponents(coordSectionA, 0, &spaceDim));
1610: PetscCall(PetscSectionSetFieldComponents(coordSectionB, 0, spaceDim));
1611: PetscCall(PetscSectionGetChart(coordSectionA, &cS, &cE));
1612: if (cStartA <= cS && cS < cEndA) { /* localized coordinates */
1613: PetscCheck((cEndA - cStartA) == (cEndB - cStartB), PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "The number of cells in first DM %" PetscInt_FMT " != %" PetscInt_FMT " in the second DM", cEndA - cStartA, cEndB - cStartB);
1614: cS = cS - cStartA + cStartB;
1615: cE = vEndB;
1616: lc = PETSC_TRUE;
1617: } else {
1618: cS = vStartB;
1619: cE = vEndB;
1620: }
1621: PetscCall(PetscSectionSetChart(coordSectionB, cS, cE));
1622: for (v = vStartB; v < vEndB; ++v) {
1623: PetscCall(PetscSectionSetDof(coordSectionB, v, spaceDim));
1624: PetscCall(PetscSectionSetFieldDof(coordSectionB, v, 0, spaceDim));
1625: }
1626: if (lc) { /* localized coordinates */
1627: PetscInt c;
1629: for (c = cS - cStartB; c < cEndB - cStartB; c++) {
1630: PetscInt dof;
1632: PetscCall(PetscSectionGetDof(coordSectionA, c + cStartA, &dof));
1633: PetscCall(PetscSectionSetDof(coordSectionB, c + cStartB, dof));
1634: PetscCall(PetscSectionSetFieldDof(coordSectionB, c + cStartB, 0, dof));
1635: }
1636: }
1637: PetscCall(PetscSectionSetUp(coordSectionB));
1638: PetscCall(PetscSectionGetStorageSize(coordSectionB, &coordSizeB));
1639: PetscCall(DMGetCoordinatesLocal(dmA, &coordinatesA));
1640: PetscCall(VecCreate(PETSC_COMM_SELF, &coordinatesB));
1641: PetscCall(PetscObjectSetName((PetscObject)coordinatesB, "coordinates"));
1642: PetscCall(VecSetSizes(coordinatesB, coordSizeB, PETSC_DETERMINE));
1643: PetscCall(VecGetBlockSize(coordinatesA, &d));
1644: PetscCall(VecSetBlockSize(coordinatesB, d));
1645: PetscCall(VecGetType(coordinatesA, &vtype));
1646: PetscCall(VecSetType(coordinatesB, vtype));
1647: PetscCall(VecGetArray(coordinatesA, &coordsA));
1648: PetscCall(VecGetArray(coordinatesB, &coordsB));
1649: for (v = 0; v < vEndB - vStartB; ++v) {
1650: PetscInt offA, offB;
1652: PetscCall(PetscSectionGetOffset(coordSectionA, v + vStartA, &offA));
1653: PetscCall(PetscSectionGetOffset(coordSectionB, v + vStartB, &offB));
1654: for (d = 0; d < spaceDim; ++d) coordsB[offB + d] = coordsA[offA + d];
1655: }
1656: if (lc) { /* localized coordinates */
1657: PetscInt c;
1659: for (c = cS - cStartB; c < cEndB - cStartB; c++) {
1660: PetscInt dof, offA, offB;
1662: PetscCall(PetscSectionGetOffset(coordSectionA, c + cStartA, &offA));
1663: PetscCall(PetscSectionGetOffset(coordSectionB, c + cStartB, &offB));
1664: PetscCall(PetscSectionGetDof(coordSectionA, c + cStartA, &dof));
1665: PetscCall(PetscArraycpy(coordsB + offB, coordsA + offA, dof));
1666: }
1667: }
1668: PetscCall(VecRestoreArray(coordinatesA, &coordsA));
1669: PetscCall(VecRestoreArray(coordinatesB, &coordsB));
1670: PetscCall(DMSetCoordinatesLocal(dmB, coordinatesB));
1671: PetscCall(VecDestroy(&coordinatesB));
1672: PetscFunctionReturn(PETSC_SUCCESS);
1673: }
1675: /*@
1676: DMPlexUninterpolate - Take in a mesh with all intermediate faces, edges, etc. and return a cell-vertex mesh
1678: Collective
1680: Input Parameter:
1681: . dm - The complete `DMPLEX` object
1683: Output Parameter:
1684: . dmUnint - The `DMPLEX` object with only cells and vertices
1686: Level: intermediate
1688: Note:
1689: It does not copy over the coordinates.
1691: Developer Notes:
1692: Sets plex->interpolated = `DMPLEX_INTERPOLATED_NONE`.
1694: .seealso: `DMPLEX`, `DMPlexInterpolate()`, `DMPlexCreateFromCellListPetsc()`, `DMPlexCopyCoordinates()`
1695: @*/
1696: PetscErrorCode DMPlexUninterpolate(DM dm, DM *dmUnint)
1697: {
1698: DMPlexInterpolatedFlag interpolated;
1699: DM udm;
1700: PetscInt dim, vStart, vEnd, cStart, cEnd, c, maxConeSize = 0, *cone;
1702: PetscFunctionBegin;
1704: PetscAssertPointer(dmUnint, 2);
1705: PetscCall(PetscLogEventBegin(DMPLEX_Uninterpolate, dm, 0, 0, 0));
1706: PetscCall(DMGetDimension(dm, &dim));
1707: PetscCall(DMPlexIsInterpolated(dm, &interpolated));
1708: PetscCheck(interpolated != DMPLEX_INTERPOLATED_PARTIAL, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Not for partially interpolated meshes");
1709: if (interpolated == DMPLEX_INTERPOLATED_NONE || dim <= 1) {
1710: /* in case dim <= 1 just keep the DMPLEX_INTERPOLATED_FULL flag */
1711: PetscCall(PetscObjectReference((PetscObject)dm));
1712: *dmUnint = dm;
1713: PetscFunctionReturn(PETSC_SUCCESS);
1714: }
1715: PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd));
1716: PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
1717: PetscCall(DMCreate(PetscObjectComm((PetscObject)dm), &udm));
1718: PetscCall(DMSetType(udm, DMPLEX));
1719: PetscCall(DMSetDimension(udm, dim));
1720: PetscCall(DMPlexSetChart(udm, cStart, vEnd));
1721: for (c = cStart; c < cEnd; ++c) {
1722: PetscInt *closure = NULL, closureSize, cl, coneSize = 0;
1724: PetscCall(DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure));
1725: for (cl = 0; cl < closureSize * 2; cl += 2) {
1726: const PetscInt p = closure[cl];
1728: if ((p >= vStart) && (p < vEnd)) ++coneSize;
1729: }
1730: PetscCall(DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure));
1731: PetscCall(DMPlexSetConeSize(udm, c, coneSize));
1732: maxConeSize = PetscMax(maxConeSize, coneSize);
1733: }
1734: PetscCall(DMSetUp(udm));
1735: PetscCall(PetscMalloc1(maxConeSize, &cone));
1736: for (c = cStart; c < cEnd; ++c) {
1737: PetscInt *closure = NULL, closureSize, cl, coneSize = 0;
1739: PetscCall(DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure));
1740: for (cl = 0; cl < closureSize * 2; cl += 2) {
1741: const PetscInt p = closure[cl];
1743: if ((p >= vStart) && (p < vEnd)) cone[coneSize++] = p;
1744: }
1745: PetscCall(DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure));
1746: PetscCall(DMPlexSetCone(udm, c, cone));
1747: }
1748: PetscCall(PetscFree(cone));
1749: PetscCall(DMPlexSymmetrize(udm));
1750: PetscCall(DMPlexStratify(udm));
1751: /* Reduce SF */
1752: {
1753: PetscSF sfPoint, sfPointUn;
1754: const PetscSFNode *remotePoints;
1755: const PetscInt *localPoints;
1756: PetscSFNode *remotePointsUn;
1757: PetscInt *localPointsUn;
1758: PetscInt numRoots, numLeaves, l;
1759: PetscInt numLeavesUn = 0, n = 0;
1761: /* Get original SF information */
1762: PetscCall(DMGetPointSF(dm, &sfPoint));
1763: if (PetscDefined(USE_DEBUG)) PetscCall(DMPlexCheckPointSF(dm, sfPoint, PETSC_FALSE));
1764: PetscCall(DMGetPointSF(udm, &sfPointUn));
1765: PetscCall(PetscSFGetGraph(sfPoint, &numRoots, &numLeaves, &localPoints, &remotePoints));
1766: if (numRoots >= 0) {
1767: /* Allocate space for cells and vertices */
1768: for (l = 0; l < numLeaves; ++l) {
1769: const PetscInt p = localPoints[l];
1771: if ((vStart <= p && p < vEnd) || (cStart <= p && p < cEnd)) numLeavesUn++;
1772: }
1773: /* Fill in leaves */
1774: PetscCall(PetscMalloc1(numLeavesUn, &remotePointsUn));
1775: PetscCall(PetscMalloc1(numLeavesUn, &localPointsUn));
1776: for (l = 0; l < numLeaves; l++) {
1777: const PetscInt p = localPoints[l];
1779: if ((vStart <= p && p < vEnd) || (cStart <= p && p < cEnd)) {
1780: localPointsUn[n] = p;
1781: remotePointsUn[n].rank = remotePoints[l].rank;
1782: remotePointsUn[n].index = remotePoints[l].index;
1783: ++n;
1784: }
1785: }
1786: PetscCheck(n == numLeavesUn, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Inconsistent number of leaves %" PetscInt_FMT " != %" PetscInt_FMT, n, numLeavesUn);
1787: PetscCall(PetscSFSetGraph(sfPointUn, cEnd - cStart + vEnd - vStart, numLeavesUn, localPointsUn, PETSC_OWN_POINTER, remotePointsUn, PETSC_OWN_POINTER));
1788: }
1789: }
1790: /* This function makes the mesh fully uninterpolated on all ranks */
1791: {
1792: DM_Plex *plex = (DM_Plex *)udm->data;
1793: plex->interpolated = plex->interpolatedCollective = DMPLEX_INTERPOLATED_NONE;
1794: }
1795: PetscCall(DMPlexCopy_Internal(dm, PETSC_TRUE, PETSC_TRUE, udm));
1796: if (PetscDefined(USE_DEBUG)) PetscCall(DMPlexCheckPointSF(udm, NULL, PETSC_FALSE));
1797: *dmUnint = udm;
1798: PetscCall(PetscLogEventEnd(DMPLEX_Uninterpolate, dm, 0, 0, 0));
1799: PetscFunctionReturn(PETSC_SUCCESS);
1800: }
1802: static PetscErrorCode DMPlexIsInterpolated_Internal(DM dm, DMPlexInterpolatedFlag *interpolated)
1803: {
1804: PetscInt coneSize, depth, dim, h, p, pStart, pEnd;
1805: MPI_Comm comm;
1807: PetscFunctionBegin;
1808: PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
1809: PetscCall(DMPlexGetDepth(dm, &depth));
1810: PetscCall(DMGetDimension(dm, &dim));
1812: if (depth == dim) {
1813: *interpolated = DMPLEX_INTERPOLATED_FULL;
1814: if (!dim) goto finish;
1816: /* Check points at height = dim are vertices (have no cones) */
1817: PetscCall(DMPlexGetHeightStratum(dm, dim, &pStart, &pEnd));
1818: for (p = pStart; p < pEnd; p++) {
1819: PetscCall(DMPlexGetConeSize(dm, p, &coneSize));
1820: if (coneSize) {
1821: *interpolated = DMPLEX_INTERPOLATED_PARTIAL;
1822: goto finish;
1823: }
1824: }
1826: /* Check points at height < dim have cones */
1827: for (h = 0; h < dim; h++) {
1828: PetscCall(DMPlexGetHeightStratum(dm, h, &pStart, &pEnd));
1829: for (p = pStart; p < pEnd; p++) {
1830: PetscCall(DMPlexGetConeSize(dm, p, &coneSize));
1831: if (!coneSize) {
1832: *interpolated = DMPLEX_INTERPOLATED_PARTIAL;
1833: goto finish;
1834: }
1835: }
1836: }
1837: } else if (depth == 1) {
1838: *interpolated = DMPLEX_INTERPOLATED_NONE;
1839: } else {
1840: *interpolated = DMPLEX_INTERPOLATED_PARTIAL;
1841: }
1842: finish:
1843: PetscFunctionReturn(PETSC_SUCCESS);
1844: }
1846: /*@
1847: DMPlexIsInterpolated - Find out to what extent the `DMPLEX` is topologically interpolated.
1849: Not Collective
1851: Input Parameter:
1852: . dm - The `DMPLEX` object
1854: Output Parameter:
1855: . interpolated - Flag whether the `DM` is interpolated
1857: Level: intermediate
1859: Notes:
1860: Unlike `DMPlexIsInterpolatedCollective()`, this is NOT collective
1861: so the results can be different on different ranks in special cases.
1862: However, `DMPlexInterpolate()` guarantees the result is the same on all.
1864: Unlike `DMPlexIsInterpolatedCollective()`, this cannot return `DMPLEX_INTERPOLATED_MIXED`.
1866: Developer Notes:
1867: Initially, plex->interpolated = `DMPLEX_INTERPOLATED_INVALID`.
1869: If plex->interpolated == `DMPLEX_INTERPOLATED_INVALID`, `DMPlexIsInterpolated_Internal()` is called.
1870: It checks the actual topology and sets plex->interpolated on each rank separately to one of
1871: `DMPLEX_INTERPOLATED_NONE`, `DMPLEX_INTERPOLATED_PARTIAL` or `DMPLEX_INTERPOLATED_FULL`.
1873: If plex->interpolated != `DMPLEX_INTERPOLATED_INVALID`, this function just returns plex->interpolated.
1875: `DMPlexInterpolate()` sets plex->interpolated = `DMPLEX_INTERPOLATED_FULL`,
1876: and DMPlexUninterpolate() sets plex->interpolated = `DMPLEX_INTERPOLATED_NONE`.
1878: .seealso: `DMPLEX`, `DMPlexInterpolate()`, `DMPlexIsInterpolatedCollective()`
1879: @*/
1880: PetscErrorCode DMPlexIsInterpolated(DM dm, DMPlexInterpolatedFlag *interpolated)
1881: {
1882: DM_Plex *plex = (DM_Plex *)dm->data;
1884: PetscFunctionBegin;
1886: PetscAssertPointer(interpolated, 2);
1887: if (plex->interpolated < 0) {
1888: PetscCall(DMPlexIsInterpolated_Internal(dm, &plex->interpolated));
1889: } else if (PetscDefined(USE_DEBUG)) {
1890: DMPlexInterpolatedFlag flg;
1892: PetscCall(DMPlexIsInterpolated_Internal(dm, &flg));
1893: PetscCheck(plex->tr || flg == plex->interpolated, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Stashed DMPlexInterpolatedFlag %s is inconsistent with current %s", DMPlexInterpolatedFlags[plex->interpolated], DMPlexInterpolatedFlags[flg]);
1894: }
1895: *interpolated = plex->interpolated;
1896: PetscFunctionReturn(PETSC_SUCCESS);
1897: }
1899: /*@
1900: DMPlexIsInterpolatedCollective - Find out to what extent the `DMPLEX` is topologically interpolated (in collective manner).
1902: Collective
1904: Input Parameter:
1905: . dm - The `DMPLEX` object
1907: Output Parameter:
1908: . interpolated - Flag whether the `DM` is interpolated
1910: Level: intermediate
1912: Notes:
1913: Unlike `DMPlexIsInterpolated()`, this is collective so the results are guaranteed to be the same on all ranks.
1915: This function will return `DMPLEX_INTERPOLATED_MIXED` if the results of `DMPlexIsInterpolated()` are different on different ranks.
1917: Developer Notes:
1918: Initially, plex->interpolatedCollective = `DMPLEX_INTERPOLATED_INVALID`.
1920: If plex->interpolatedCollective == `DMPLEX_INTERPOLATED_INVALID`, this function calls `DMPlexIsInterpolated()` which sets plex->interpolated.
1921: `MPI_Allreduce()` is then called and collectively consistent flag plex->interpolatedCollective is set and returned;
1922: if plex->interpolated varies on different ranks, plex->interpolatedCollective = `DMPLEX_INTERPOLATED_MIXED`,
1923: otherwise sets plex->interpolatedCollective = plex->interpolated.
1925: If plex->interpolatedCollective != `DMPLEX_INTERPOLATED_INVALID`, this function just returns plex->interpolatedCollective.
1927: .seealso: `DMPLEX`, `DMPlexInterpolate()`, `DMPlexIsInterpolated()`
1928: @*/
1929: PetscErrorCode DMPlexIsInterpolatedCollective(DM dm, DMPlexInterpolatedFlag *interpolated)
1930: {
1931: DM_Plex *plex = (DM_Plex *)dm->data;
1932: PetscBool debug = PETSC_FALSE;
1934: PetscFunctionBegin;
1936: PetscAssertPointer(interpolated, 2);
1937: PetscCall(PetscOptionsGetBool(((PetscObject)dm)->options, ((PetscObject)dm)->prefix, "-dm_plex_is_interpolated_collective_debug", &debug, NULL));
1938: if (plex->interpolatedCollective < 0) {
1939: DMPlexInterpolatedFlag min, max;
1940: MPI_Comm comm;
1942: PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
1943: PetscCall(DMPlexIsInterpolated(dm, &plex->interpolatedCollective));
1944: PetscCall(MPIU_Allreduce(&plex->interpolatedCollective, &min, 1, MPIU_ENUM, MPI_MIN, comm));
1945: PetscCall(MPIU_Allreduce(&plex->interpolatedCollective, &max, 1, MPIU_ENUM, MPI_MAX, comm));
1946: if (min != max) plex->interpolatedCollective = DMPLEX_INTERPOLATED_MIXED;
1947: if (debug) {
1948: PetscMPIInt rank;
1950: PetscCallMPI(MPI_Comm_rank(comm, &rank));
1951: PetscCall(PetscSynchronizedPrintf(comm, "[%d] interpolated=%s interpolatedCollective=%s\n", rank, DMPlexInterpolatedFlags[plex->interpolated], DMPlexInterpolatedFlags[plex->interpolatedCollective]));
1952: PetscCall(PetscSynchronizedFlush(comm, PETSC_STDOUT));
1953: }
1954: }
1955: *interpolated = plex->interpolatedCollective;
1956: PetscFunctionReturn(PETSC_SUCCESS);
1957: }