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