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