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