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, vStart, vEnd;
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, 0, &vStart, &vEnd));
506: PetscCall(DMPlexGetDepthStratum(dm, cellDepth, &cStart, &cEnd));
507: // If the range incorporates the vertices, it means we have a non-manifold topology, so choose just cells
508: if (cStart <= vStart && cEnd >= vEnd) cEnd = vStart;
509: // Number new faces and save face vertices in hash table
510: // If depth > cellDepth, meaning we are interpolating faces, put the new (d-1)-faces after them
511: // otherwise, we are interpolating cells, so put the faces after the vertices
512: PetscCall(DMPlexGetDepthStratum(dm, depth > cellDepth ? cellDepth : 0, NULL, &fStart));
513: fEnd = fStart;
515: minCone = PETSC_MAX_INT;
516: cntFaces = 0;
517: for (PetscInt c = cStart; c < cEnd; ++c) {
518: const PetscInt *cone;
519: DMPolytopeType ct;
520: PetscInt numFaces = 0, coneSize;
522: PetscCall(DMPlexGetCellType(dm, c, &ct));
523: PetscCall(DMPlexGetCone(dm, c, &cone));
524: PetscCall(DMPlexGetConeSize(dm, c, &coneSize));
525: for (PetscInt j = 0; j < coneSize; j++) minCone = PetscMin(cone[j], minCone);
526: // Ignore faces since they are interpolated
527: if (ct != DM_POLYTOPE_SEGMENT && ct != DM_POLYTOPE_POINT_PRISM_TENSOR) PetscCall(DMPlexGetRawFaces_Internal(dm, ct, cone, &numFaces, NULL, NULL, NULL));
528: cntFaces += numFaces;
529: }
530: // Encode so that we can use 0 as an excluded value, instead of PETSC_MAX_INT
531: minCone = -(minCone - 1);
533: PetscCall(PetscMalloc1(cntFaces, &facesId));
535: cntFaces = 0;
536: for (PetscInt c = cStart; c < cEnd; ++c) {
537: const PetscInt *cone, *faceSizes, *faces;
538: const DMPolytopeType *faceTypes;
539: DMPolytopeType ct;
540: PetscInt numFaces, foff = 0;
542: PetscCall(DMPlexGetCellType(dm, c, &ct));
543: PetscCall(DMPlexGetCone(dm, c, &cone));
544: // Ignore faces since they are interpolated
545: if (ct != DM_POLYTOPE_SEGMENT && ct != DM_POLYTOPE_POINT_PRISM_TENSOR) {
546: PetscCall(DMPlexGetRawFaces_Internal(dm, ct, cone, &numFaces, &faceTypes, &faceSizes, &faces));
547: } else {
548: numFaces = 0;
549: }
550: for (PetscInt cf = 0; cf < numFaces; foff += faceSizes[cf], ++cf) {
551: const PetscInt faceSize = faceSizes[cf];
552: const DMPolytopeType faceType = faceTypes[cf];
553: const PetscInt *face = &faces[foff];
554: PetscHashIJKLKey key;
555: PetscHashIter iter;
556: PetscBool missing;
558: PetscCheck(faceSize <= 4, PETSC_COMM_SELF, PETSC_ERR_SUP, "Do not support faces of size %" PetscInt_FMT " > 4", faceSize);
559: key.i = face[0] + minCone;
560: key.j = faceSize > 1 ? face[1] + minCone : 0;
561: key.k = faceSize > 2 ? face[2] + minCone : 0;
562: key.l = faceSize > 3 ? face[3] + minCone : 0;
563: PetscCall(PetscSortInt(faceSize, (PetscInt *)&key));
564: PetscCall(PetscHMapIJKLPut(faceTable, key, &iter, &missing));
565: if (missing) {
566: facesId[cntFaces] = fEnd;
567: PetscCall(PetscHMapIJKLIterSet(faceTable, iter, fEnd++));
568: ++faceTypeNum[faceType];
569: } else PetscCall(PetscHMapIJKLIterGet(faceTable, iter, &facesId[cntFaces]));
570: cntFaces++;
571: }
572: if (ct != DM_POLYTOPE_SEGMENT && ct != DM_POLYTOPE_POINT_PRISM_TENSOR) PetscCall(DMPlexRestoreRawFaces_Internal(dm, ct, cone, &numFaces, &faceTypes, &faceSizes, &faces));
573: }
574: /* We need to number faces contiguously among types */
575: {
576: PetscInt faceTypeStart[DM_NUM_POLYTOPES], ct, numFT = 0;
578: for (ct = 0; ct < DM_NUM_POLYTOPES; ++ct) {
579: if (faceTypeNum[ct]) ++numFT;
580: faceTypeStart[ct] = 0;
581: }
582: if (numFT > 1) {
583: PetscCall(PetscHMapIJKLClear(faceTable));
584: faceTypeStart[0] = fStart;
585: for (ct = 1; ct < DM_NUM_POLYTOPES; ++ct) faceTypeStart[ct] = faceTypeStart[ct - 1] + faceTypeNum[ct - 1];
586: cntFaces = 0;
587: for (PetscInt c = cStart; c < cEnd; ++c) {
588: const PetscInt *cone, *faceSizes, *faces;
589: const DMPolytopeType *faceTypes;
590: DMPolytopeType ct;
591: PetscInt numFaces, foff = 0;
593: PetscCall(DMPlexGetCellType(dm, c, &ct));
594: PetscCall(DMPlexGetCone(dm, c, &cone));
595: if (ct != DM_POLYTOPE_SEGMENT && ct != DM_POLYTOPE_POINT_PRISM_TENSOR) {
596: PetscCall(DMPlexGetRawFaces_Internal(dm, ct, cone, &numFaces, &faceTypes, &faceSizes, &faces));
597: } else {
598: numFaces = 0;
599: }
600: for (PetscInt cf = 0; cf < numFaces; foff += faceSizes[cf], ++cf) {
601: const PetscInt faceSize = faceSizes[cf];
602: const DMPolytopeType faceType = faceTypes[cf];
603: const PetscInt *face = &faces[foff];
604: PetscHashIJKLKey key;
605: PetscHashIter iter;
606: PetscBool missing;
608: key.i = face[0] + minCone;
609: key.j = faceSize > 1 ? face[1] + minCone : 0;
610: key.k = faceSize > 2 ? face[2] + minCone : 0;
611: key.l = faceSize > 3 ? face[3] + minCone : 0;
612: PetscCall(PetscSortInt(faceSize, (PetscInt *)&key));
613: PetscCall(PetscHMapIJKLPut(faceTable, key, &iter, &missing));
614: if (missing) {
615: facesId[cntFaces] = faceTypeStart[faceType];
616: PetscCall(PetscHMapIJKLIterSet(faceTable, iter, faceTypeStart[faceType]++));
617: } else PetscCall(PetscHMapIJKLIterGet(faceTable, iter, &facesId[cntFaces]));
618: cntFaces++;
619: }
620: if (ct != DM_POLYTOPE_SEGMENT && ct != DM_POLYTOPE_POINT_PRISM_TENSOR) PetscCall(DMPlexRestoreRawFaces_Internal(dm, ct, cone, &numFaces, &faceTypes, &faceSizes, &faces));
621: }
622: for (ct = 1; ct < DM_NUM_POLYTOPES; ++ct) {
623: 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]);
624: }
625: }
626: }
627: PetscCall(PetscHMapIJKLDestroy(&faceTable));
629: // Add new points, perhaps inserting into the numbering
630: PetscCall(DMPlexGetChart(dm, &pStart, &Np));
631: PetscCall(DMPlexSetChart(idm, pStart, Np + (fEnd - fStart)));
632: // Set cone sizes
633: // Must create the celltype label here so that we do not automatically try to compute the types
634: PetscCall(DMCreateLabel(idm, "celltype"));
635: PetscCall(DMPlexGetCellTypeLabel(idm, &ctLabel));
636: for (PetscInt d = 0; d <= depth; ++d) {
637: DMPolytopeType ct;
638: PetscInt coneSize, pStart, pEnd, poff = 0;
640: PetscCall(DMPlexGetDepthStratum(dm, d, &pStart, &pEnd));
641: // Check for non-manifold condition
642: if (d == cellDepth) {
643: if (pEnd == cEnd) continue;
644: else pStart = vEnd;
645: }
646: // Account for insertion
647: if (pStart >= fStart) poff = fEnd - fStart;
648: for (PetscInt p = pStart; p < pEnd; ++p) {
649: PetscCall(DMPlexGetConeSize(dm, p, &coneSize));
650: PetscCall(DMPlexSetConeSize(idm, p + poff, coneSize));
651: PetscCall(DMPlexGetCellType(dm, p, &ct));
652: PetscCall(DMPlexSetCellType(idm, p + poff, ct));
653: }
654: }
655: cntFaces = 0;
656: for (PetscInt c = cStart; c < cEnd; ++c) {
657: const PetscInt *cone, *faceSizes;
658: const DMPolytopeType *faceTypes;
659: DMPolytopeType ct;
660: PetscInt numFaces, poff = 0;
662: PetscCall(DMPlexGetCellType(dm, c, &ct));
663: PetscCall(DMPlexGetCone(dm, c, &cone));
664: if (c >= fStart) poff = fEnd - fStart;
665: if (ct == DM_POLYTOPE_SEGMENT || ct == DM_POLYTOPE_POINT_PRISM_TENSOR) {
666: PetscCall(DMPlexSetCellType(idm, c + poff, ct));
667: PetscCall(DMPlexSetConeSize(idm, c + poff, 2));
668: continue;
669: }
670: PetscCall(DMPlexGetRawFaces_Internal(dm, ct, cone, &numFaces, &faceTypes, &faceSizes, NULL));
671: PetscCall(DMPlexSetCellType(idm, c + poff, ct));
672: PetscCall(DMPlexSetConeSize(idm, c + poff, numFaces));
673: for (PetscInt cf = 0; cf < numFaces; ++cf) {
674: const PetscInt f = facesId[cntFaces];
675: DMPolytopeType faceType = faceTypes[cf];
676: const PetscInt faceSize = faceSizes[cf];
677: PetscCall(DMPlexSetConeSize(idm, f, faceSize));
678: PetscCall(DMPlexSetCellType(idm, f, faceType));
679: cntFaces++;
680: }
681: PetscCall(DMPlexRestoreRawFaces_Internal(dm, ct, cone, &numFaces, &faceTypes, &faceSizes, NULL));
682: }
683: PetscCall(DMSetUp(idm));
684: // Initialize cones so we do not need the bash table to tell us that a cone has been set
685: {
686: PetscSection cs;
687: PetscInt *cones, csize;
689: PetscCall(DMPlexGetConeSection(idm, &cs));
690: PetscCall(DMPlexGetCones(idm, &cones));
691: PetscCall(PetscSectionGetStorageSize(cs, &csize));
692: for (PetscInt c = 0; c < csize; ++c) cones[c] = -1;
693: }
694: // Set cones
695: {
696: PetscInt *icone;
697: PetscInt maxConeSize;
699: PetscCall(DMPlexGetMaxSizes(dm, &maxConeSize, NULL));
700: PetscCall(PetscMalloc1(maxConeSize, &icone));
701: for (PetscInt d = 0; d <= depth; ++d) {
702: const PetscInt *cone;
703: PetscInt pStart, pEnd, poff = 0, coneSize;
705: PetscCall(DMPlexGetDepthStratum(dm, d, &pStart, &pEnd));
706: // Check for non-manifold condition
707: if (d == cellDepth) {
708: if (pEnd == cEnd) continue;
709: else pStart = vEnd;
710: }
711: // Account for insertion
712: if (pStart >= fStart) poff = fEnd - fStart;
713: for (PetscInt p = pStart; p < pEnd; ++p) {
714: PetscCall(DMPlexGetCone(dm, p, &cone));
715: PetscCall(DMPlexGetConeSize(dm, p, &coneSize));
716: for (PetscInt cp = 0; cp < coneSize; ++cp) icone[cp] = cone[cp] + (cone[cp] >= fStart ? fEnd - fStart : 0);
717: PetscCall(DMPlexSetCone(idm, p + poff, icone));
718: PetscCall(DMPlexGetConeOrientation(dm, p, &cone));
719: PetscCall(DMPlexSetConeOrientation(idm, p + poff, cone));
720: }
721: }
722: cntFaces = 0;
723: for (PetscInt c = cStart; c < cEnd; ++c) {
724: const PetscInt *cone, *faceSizes, *faces;
725: const DMPolytopeType *faceTypes;
726: DMPolytopeType ct;
727: PetscInt coneSize, numFaces, foff = 0, poff = 0;
729: PetscCall(DMPlexGetCellType(dm, c, &ct));
730: PetscCall(DMPlexGetCone(dm, c, &cone));
731: PetscCall(DMPlexGetConeSize(dm, c, &coneSize));
732: if (c >= fStart) poff = fEnd - fStart;
733: if (ct == DM_POLYTOPE_SEGMENT || ct == DM_POLYTOPE_POINT_PRISM_TENSOR) {
734: for (PetscInt cp = 0; cp < coneSize; ++cp) icone[cp] = cone[cp] + (cone[cp] >= fStart ? fEnd - fStart : 0);
735: PetscCall(DMPlexSetCone(idm, c + poff, icone));
736: PetscCall(DMPlexGetConeOrientation(dm, c, &cone));
737: PetscCall(DMPlexSetConeOrientation(idm, c + poff, cone));
738: continue;
739: }
740: PetscCall(DMPlexGetRawFaces_Internal(dm, ct, cone, &numFaces, &faceTypes, &faceSizes, &faces));
741: for (PetscInt cf = 0; cf < numFaces; foff += faceSizes[cf], ++cf) {
742: DMPolytopeType faceType = faceTypes[cf];
743: const PetscInt faceSize = faceSizes[cf];
744: const PetscInt f = facesId[cntFaces];
745: const PetscInt *face = &faces[foff];
746: const PetscInt *fcone;
748: PetscCall(DMPlexInsertCone(idm, c, cf, f));
749: PetscCall(DMPlexGetCone(idm, f, &fcone));
750: if (fcone[0] < 0) PetscCall(DMPlexSetCone(idm, f, face));
751: {
752: const PetscInt *fcone2;
753: PetscInt ornt;
755: PetscCall(DMPlexGetConeSize(idm, f, &coneSize));
756: PetscCall(DMPlexGetCone(idm, f, &fcone2));
757: 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);
758: /* Notice that we have to use vertices here because the lower dimensional faces have not been created yet */
759: PetscCall(DMPolytopeGetVertexOrientation(faceType, fcone2, face, &ornt));
760: PetscCall(DMPlexInsertConeOrientation(idm, c + poff, cf, ornt));
761: }
762: cntFaces++;
763: }
764: PetscCall(DMPlexRestoreRawFaces_Internal(dm, ct, cone, &numFaces, &faceTypes, &faceSizes, &faces));
765: }
766: PetscCall(PetscFree(icone));
767: }
768: PetscCall(PetscFree(facesId));
769: PetscCall(DMPlexSymmetrize(idm));
770: PetscCall(DMPlexStratify(idm));
771: PetscFunctionReturn(PETSC_SUCCESS);
772: }
774: static PetscErrorCode SortRmineRremoteByRemote_Private(PetscSF sf, PetscInt *rmine1[], PetscInt *rremote1[])
775: {
776: PetscInt nleaves;
777: PetscInt nranks;
778: const PetscMPIInt *ranks = NULL;
779: const PetscInt *roffset = NULL, *rmine = NULL, *rremote = NULL;
780: PetscInt n, o, r;
782: PetscFunctionBegin;
783: PetscCall(PetscSFGetRootRanks(sf, &nranks, &ranks, &roffset, &rmine, &rremote));
784: nleaves = roffset[nranks];
785: PetscCall(PetscMalloc2(nleaves, rmine1, nleaves, rremote1));
786: for (r = 0; r < nranks; r++) {
787: /* simultaneously sort rank-wise portions of rmine & rremote by values in rremote
788: - to unify order with the other side */
789: o = roffset[r];
790: n = roffset[r + 1] - o;
791: PetscCall(PetscArraycpy(&(*rmine1)[o], &rmine[o], n));
792: PetscCall(PetscArraycpy(&(*rremote1)[o], &rremote[o], n));
793: PetscCall(PetscSortIntWithArray(n, &(*rremote1)[o], &(*rmine1)[o]));
794: }
795: PetscFunctionReturn(PETSC_SUCCESS);
796: }
798: PetscErrorCode DMPlexOrientInterface_Internal(DM dm)
799: {
800: PetscSF sf;
801: const PetscInt *locals;
802: const PetscSFNode *remotes;
803: const PetscMPIInt *ranks;
804: const PetscInt *roffset;
805: PetscInt *rmine1, *rremote1; /* rmine and rremote copies simultaneously sorted by rank and rremote */
806: PetscInt nroots, p, nleaves, nranks, r, maxConeSize = 0;
807: PetscInt(*roots)[4], (*leaves)[4], mainCone[4];
808: PetscMPIInt(*rootsRanks)[4], (*leavesRanks)[4];
809: MPI_Comm comm;
810: PetscMPIInt rank, size;
811: PetscInt debug = 0;
813: PetscFunctionBegin;
814: PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
815: PetscCallMPI(MPI_Comm_rank(comm, &rank));
816: PetscCallMPI(MPI_Comm_size(comm, &size));
817: PetscCall(DMGetPointSF(dm, &sf));
818: PetscCall(DMViewFromOptions(dm, NULL, "-before_orient_interface_dm_view"));
819: if (PetscDefined(USE_DEBUG)) PetscCall(DMPlexCheckPointSF(dm, sf, PETSC_FALSE));
820: PetscCall(PetscSFGetGraph(sf, &nroots, &nleaves, &locals, &remotes));
821: if (nroots < 0) PetscFunctionReturn(PETSC_SUCCESS);
822: PetscCall(PetscSFSetUp(sf));
823: PetscCall(SortRmineRremoteByRemote_Private(sf, &rmine1, &rremote1));
824: for (p = 0; p < nleaves; ++p) {
825: PetscInt coneSize;
826: PetscCall(DMPlexGetConeSize(dm, locals[p], &coneSize));
827: maxConeSize = PetscMax(maxConeSize, coneSize);
828: }
829: PetscCheck(maxConeSize <= 4, comm, PETSC_ERR_SUP, "This method does not support cones of size %" PetscInt_FMT, maxConeSize);
830: PetscCall(PetscMalloc4(nroots, &roots, nroots, &leaves, nroots, &rootsRanks, nroots, &leavesRanks));
831: for (p = 0; p < nroots; ++p) {
832: const PetscInt *cone;
833: PetscInt coneSize, c, ind0;
835: PetscCall(DMPlexGetConeSize(dm, p, &coneSize));
836: PetscCall(DMPlexGetCone(dm, p, &cone));
837: /* Ignore vertices */
838: if (coneSize < 2) {
839: for (c = 0; c < 4; c++) {
840: roots[p][c] = -1;
841: rootsRanks[p][c] = -1;
842: }
843: continue;
844: }
845: /* Translate all points to root numbering */
846: for (c = 0; c < PetscMin(coneSize, 4); c++) {
847: PetscCall(PetscFindInt(cone[c], nleaves, locals, &ind0));
848: if (ind0 < 0) {
849: roots[p][c] = cone[c];
850: rootsRanks[p][c] = rank;
851: } else {
852: roots[p][c] = remotes[ind0].index;
853: rootsRanks[p][c] = remotes[ind0].rank;
854: }
855: }
856: for (c = coneSize; c < 4; c++) {
857: roots[p][c] = -1;
858: rootsRanks[p][c] = -1;
859: }
860: }
861: for (p = 0; p < nroots; ++p) {
862: PetscInt c;
863: for (c = 0; c < 4; c++) {
864: leaves[p][c] = -2;
865: leavesRanks[p][c] = -2;
866: }
867: }
868: PetscCall(PetscSFBcastBegin(sf, MPIU_4INT, roots, leaves, MPI_REPLACE));
869: PetscCall(PetscSFBcastBegin(sf, MPI_4INT, rootsRanks, leavesRanks, MPI_REPLACE));
870: PetscCall(PetscSFBcastEnd(sf, MPIU_4INT, roots, leaves, MPI_REPLACE));
871: PetscCall(PetscSFBcastEnd(sf, MPI_4INT, rootsRanks, leavesRanks, MPI_REPLACE));
872: if (debug) {
873: PetscCall(PetscSynchronizedFlush(comm, NULL));
874: if (rank == 0) PetscCall(PetscSynchronizedPrintf(comm, "Referenced roots\n"));
875: }
876: PetscCall(PetscSFGetRootRanks(sf, &nranks, &ranks, &roffset, NULL, NULL));
877: for (p = 0; p < nroots; ++p) {
878: DMPolytopeType ct;
879: const PetscInt *cone;
880: PetscInt coneSize, c, ind0, o;
882: if (leaves[p][0] < 0) continue; /* Ignore vertices */
883: PetscCall(DMPlexGetCellType(dm, p, &ct));
884: PetscCall(DMPlexGetConeSize(dm, p, &coneSize));
885: PetscCall(DMPlexGetCone(dm, p, &cone));
886: if (debug) {
887: 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]));
888: }
889: 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]) {
890: /* Translate these leaves to my cone points; mainCone means desired order p's cone points */
891: for (c = 0; c < PetscMin(coneSize, 4); ++c) {
892: PetscInt rS, rN;
894: if (leavesRanks[p][c] == rank) {
895: /* A local leaf is just taken as it is */
896: mainCone[c] = leaves[p][c];
897: continue;
898: }
899: /* Find index of rank leavesRanks[p][c] among remote ranks */
900: /* No need for PetscMPIIntCast because these integers were originally cast from PetscMPIInt. */
901: PetscCall(PetscFindMPIInt((PetscMPIInt)leavesRanks[p][c], nranks, ranks, &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[%" PetscInt_FMT "] = %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: PetscInt idx;
931: PetscMPIInt rank;
932: PetscBool flg;
934: PetscFunctionBegin;
935: PetscCall(PetscOptionsHasName(NULL, NULL, opt, &flg));
936: if (!flg) PetscFunctionReturn(PETSC_SUCCESS);
937: PetscCallMPI(MPI_Comm_rank(comm, &rank));
938: PetscCall(PetscSynchronizedPrintf(comm, "[%d]%s:\n", rank, name));
939: for (idx = 0; idx < n; ++idx) PetscCall(PetscSynchronizedPrintf(comm, "[%d]%s %" PetscInt_FMT " %s %" PetscInt_FMT "\n", rank, idxname, idx, valname, a[idx]));
940: PetscCall(PetscSynchronizedFlush(comm, NULL));
941: PetscFunctionReturn(PETSC_SUCCESS);
942: }
944: static PetscErrorCode SFNodeArrayViewFromOptions(MPI_Comm comm, const char opt[], const char name[], const char idxname[], PetscInt n, const PetscSFNode a[])
945: {
946: PetscInt idx;
947: PetscMPIInt rank;
948: PetscBool flg;
950: PetscFunctionBegin;
951: PetscCall(PetscOptionsHasName(NULL, NULL, opt, &flg));
952: if (!flg) PetscFunctionReturn(PETSC_SUCCESS);
953: PetscCallMPI(MPI_Comm_rank(comm, &rank));
954: PetscCall(PetscSynchronizedPrintf(comm, "[%d]%s:\n", rank, name));
955: if (idxname) {
956: 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));
957: } else {
958: for (idx = 0; idx < n; ++idx) PetscCall(PetscSynchronizedPrintf(comm, "[%d]rank %" PetscInt_FMT " index %" PetscInt_FMT "\n", rank, a[idx].rank, a[idx].index));
959: }
960: PetscCall(PetscSynchronizedFlush(comm, NULL));
961: PetscFunctionReturn(PETSC_SUCCESS);
962: }
964: static PetscErrorCode DMPlexMapToLocalPoint(DM dm, PetscHMapIJ remotehash, PetscSFNode remotePoint, PetscInt *localPoint, PetscBool *mapFailed)
965: {
966: PetscSF sf;
967: const PetscInt *locals;
968: PetscMPIInt rank;
970: PetscFunctionBegin;
971: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank));
972: PetscCall(DMGetPointSF(dm, &sf));
973: PetscCall(PetscSFGetGraph(sf, NULL, NULL, &locals, NULL));
974: if (mapFailed) *mapFailed = PETSC_FALSE;
975: if (remotePoint.rank == rank) {
976: *localPoint = remotePoint.index;
977: } else {
978: PetscHashIJKey key;
979: PetscInt l;
981: key.i = remotePoint.index;
982: key.j = remotePoint.rank;
983: PetscCall(PetscHMapIJGet(remotehash, key, &l));
984: if (l >= 0) {
985: *localPoint = locals[l];
986: } else if (mapFailed) *mapFailed = PETSC_TRUE;
987: }
988: PetscFunctionReturn(PETSC_SUCCESS);
989: }
991: static PetscErrorCode DMPlexMapToGlobalPoint(DM dm, PetscInt localPoint, PetscSFNode *remotePoint, PetscBool *mapFailed)
992: {
993: PetscSF sf;
994: const PetscInt *locals, *rootdegree;
995: const PetscSFNode *remotes;
996: PetscInt Nl, l;
997: PetscMPIInt rank;
999: PetscFunctionBegin;
1000: if (mapFailed) *mapFailed = PETSC_FALSE;
1001: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank));
1002: PetscCall(DMGetPointSF(dm, &sf));
1003: PetscCall(PetscSFGetGraph(sf, NULL, &Nl, &locals, &remotes));
1004: if (Nl < 0) goto owned;
1005: PetscCall(PetscSFComputeDegreeBegin(sf, &rootdegree));
1006: PetscCall(PetscSFComputeDegreeEnd(sf, &rootdegree));
1007: if (rootdegree[localPoint]) goto owned;
1008: PetscCall(PetscFindInt(localPoint, Nl, locals, &l));
1009: if (l < 0) {
1010: if (mapFailed) *mapFailed = PETSC_TRUE;
1011: } else *remotePoint = remotes[l];
1012: PetscFunctionReturn(PETSC_SUCCESS);
1013: owned:
1014: remotePoint->rank = rank;
1015: remotePoint->index = localPoint;
1016: PetscFunctionReturn(PETSC_SUCCESS);
1017: }
1019: static PetscErrorCode DMPlexPointIsShared(DM dm, PetscInt p, PetscBool *isShared)
1020: {
1021: PetscSF sf;
1022: const PetscInt *locals, *rootdegree;
1023: PetscInt Nl, idx;
1025: PetscFunctionBegin;
1026: *isShared = PETSC_FALSE;
1027: PetscCall(DMGetPointSF(dm, &sf));
1028: PetscCall(PetscSFGetGraph(sf, NULL, &Nl, &locals, NULL));
1029: if (Nl < 0) PetscFunctionReturn(PETSC_SUCCESS);
1030: PetscCall(PetscFindInt(p, Nl, locals, &idx));
1031: if (idx >= 0) {
1032: *isShared = PETSC_TRUE;
1033: PetscFunctionReturn(PETSC_SUCCESS);
1034: }
1035: PetscCall(PetscSFComputeDegreeBegin(sf, &rootdegree));
1036: PetscCall(PetscSFComputeDegreeEnd(sf, &rootdegree));
1037: if (rootdegree[p] > 0) *isShared = PETSC_TRUE;
1038: PetscFunctionReturn(PETSC_SUCCESS);
1039: }
1041: static PetscErrorCode DMPlexConeIsShared(DM dm, PetscInt p, PetscBool *isShared)
1042: {
1043: const PetscInt *cone;
1044: PetscInt coneSize, c;
1045: PetscBool cShared = PETSC_TRUE;
1047: PetscFunctionBegin;
1048: PetscCall(DMPlexGetConeSize(dm, p, &coneSize));
1049: PetscCall(DMPlexGetCone(dm, p, &cone));
1050: for (c = 0; c < coneSize; ++c) {
1051: PetscBool pointShared;
1053: PetscCall(DMPlexPointIsShared(dm, cone[c], &pointShared));
1054: cShared = (PetscBool)(cShared && pointShared);
1055: }
1056: *isShared = coneSize ? cShared : PETSC_FALSE;
1057: PetscFunctionReturn(PETSC_SUCCESS);
1058: }
1060: static PetscErrorCode DMPlexGetConeMinimum(DM dm, PetscInt p, PetscSFNode *cpmin)
1061: {
1062: const PetscInt *cone;
1063: PetscInt coneSize, c;
1064: PetscSFNode cmin = {PETSC_MAX_INT, PETSC_MAX_INT}, missing = {-1, -1};
1066: PetscFunctionBegin;
1067: PetscCall(DMPlexGetConeSize(dm, p, &coneSize));
1068: PetscCall(DMPlexGetCone(dm, p, &cone));
1069: for (c = 0; c < coneSize; ++c) {
1070: PetscSFNode rcp;
1071: PetscBool mapFailed;
1073: PetscCall(DMPlexMapToGlobalPoint(dm, cone[c], &rcp, &mapFailed));
1074: if (mapFailed) {
1075: cmin = missing;
1076: } else {
1077: cmin = (rcp.rank < cmin.rank) || (rcp.rank == cmin.rank && rcp.index < cmin.index) ? rcp : cmin;
1078: }
1079: }
1080: *cpmin = coneSize ? cmin : missing;
1081: PetscFunctionReturn(PETSC_SUCCESS);
1082: }
1084: /*
1085: Each shared face has an entry in the candidates array:
1086: (-1, coneSize-1), {(global cone point)}
1087: where the set is missing the point p which we use as the key for the face
1088: */
1089: static PetscErrorCode DMPlexAddSharedFace_Private(DM dm, PetscSection candidateSection, PetscSFNode candidates[], PetscHMapIJ faceHash, PetscInt p, PetscBool debug)
1090: {
1091: MPI_Comm comm;
1092: const PetscInt *support;
1093: PetscInt supportSize, s, off = 0, idx = 0, overlap, cellHeight, height;
1094: PetscMPIInt rank;
1096: PetscFunctionBegin;
1097: PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
1098: PetscCallMPI(MPI_Comm_rank(comm, &rank));
1099: PetscCall(DMPlexGetOverlap(dm, &overlap));
1100: PetscCall(DMPlexGetVTKCellHeight(dm, &cellHeight));
1101: PetscCall(DMPlexGetPointHeight(dm, p, &height));
1102: if (!overlap && height <= cellHeight + 1) {
1103: /* cells can't be shared for non-overlapping meshes */
1104: if (debug) PetscCall(PetscSynchronizedPrintf(comm, "[%d] Skipping face %" PetscInt_FMT " to avoid adding cell to hashmap since this is nonoverlapping mesh\n", rank, p));
1105: PetscFunctionReturn(PETSC_SUCCESS);
1106: }
1107: PetscCall(DMPlexGetSupportSize(dm, p, &supportSize));
1108: PetscCall(DMPlexGetSupport(dm, p, &support));
1109: if (candidates) PetscCall(PetscSectionGetOffset(candidateSection, p, &off));
1110: for (s = 0; s < supportSize; ++s) {
1111: const PetscInt face = support[s];
1112: const PetscInt *cone;
1113: PetscSFNode cpmin = {-1, -1}, rp = {-1, -1};
1114: PetscInt coneSize, c, f;
1115: PetscBool isShared = PETSC_FALSE;
1116: PetscHashIJKey key;
1118: /* Only add point once */
1119: if (debug) PetscCall(PetscSynchronizedPrintf(comm, "[%d] Support face %" PetscInt_FMT "\n", rank, face));
1120: key.i = p;
1121: key.j = face;
1122: PetscCall(PetscHMapIJGet(faceHash, key, &f));
1123: if (f >= 0) continue;
1124: PetscCall(DMPlexConeIsShared(dm, face, &isShared));
1125: PetscCall(DMPlexGetConeMinimum(dm, face, &cpmin));
1126: PetscCall(DMPlexMapToGlobalPoint(dm, p, &rp, NULL));
1127: if (debug) {
1128: PetscCall(PetscSynchronizedPrintf(comm, "[%d] Face point %" PetscInt_FMT " is shared: %d\n", rank, face, (int)isShared));
1129: 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));
1130: }
1131: if (isShared && (rp.rank == cpmin.rank && rp.index == cpmin.index)) {
1132: PetscCall(PetscHMapIJSet(faceHash, key, p));
1133: if (candidates) {
1134: if (debug) PetscCall(PetscSynchronizedPrintf(comm, "[%d] Adding shared face %" PetscInt_FMT " at idx %" PetscInt_FMT "\n[%d] ", rank, face, idx, rank));
1135: PetscCall(DMPlexGetConeSize(dm, face, &coneSize));
1136: PetscCall(DMPlexGetCone(dm, face, &cone));
1137: candidates[off + idx].rank = -1;
1138: candidates[off + idx++].index = coneSize - 1;
1139: candidates[off + idx].rank = rank;
1140: candidates[off + idx++].index = face;
1141: for (c = 0; c < coneSize; ++c) {
1142: const PetscInt cp = cone[c];
1144: if (cp == p) continue;
1145: PetscCall(DMPlexMapToGlobalPoint(dm, cp, &candidates[off + idx], NULL));
1146: if (debug) PetscCall(PetscSynchronizedPrintf(comm, " (%" PetscInt_FMT ",%" PetscInt_FMT ")", candidates[off + idx].rank, candidates[off + idx].index));
1147: ++idx;
1148: }
1149: if (debug) PetscCall(PetscSynchronizedPrintf(comm, "\n"));
1150: } else {
1151: /* Add cone size to section */
1152: if (debug) PetscCall(PetscSynchronizedPrintf(comm, "[%d] Scheduling shared face %" PetscInt_FMT "\n", rank, face));
1153: PetscCall(DMPlexGetConeSize(dm, face, &coneSize));
1154: PetscCall(PetscHMapIJSet(faceHash, key, p));
1155: PetscCall(PetscSectionAddDof(candidateSection, p, coneSize + 1));
1156: }
1157: }
1158: }
1159: PetscFunctionReturn(PETSC_SUCCESS);
1160: }
1162: /*@
1163: DMPlexInterpolatePointSF - Insert interpolated points in the overlap into the `PointSF` in parallel, following local interpolation
1165: Collective
1167: Input Parameters:
1168: + dm - The interpolated `DMPLEX`
1169: - pointSF - The initial `PetscSF` without interpolated points
1171: Level: developer
1173: Note:
1174: 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`
1176: .seealso: `DMPLEX`, `DMPlexInterpolate()`, `DMPlexUninterpolate()`
1177: @*/
1178: PetscErrorCode DMPlexInterpolatePointSF(DM dm, PetscSF pointSF)
1179: {
1180: MPI_Comm comm;
1181: PetscHMapIJ remoteHash;
1182: PetscHMapI claimshash;
1183: PetscSection candidateSection, candidateRemoteSection, claimSection;
1184: PetscSFNode *candidates, *candidatesRemote, *claims;
1185: const PetscInt *localPoints, *rootdegree;
1186: const PetscSFNode *remotePoints;
1187: PetscInt ov, Nr, r, Nl, l;
1188: PetscInt candidatesSize, candidatesRemoteSize, claimsSize;
1189: PetscBool flg, debug = PETSC_FALSE;
1190: PetscMPIInt rank;
1192: PetscFunctionBegin;
1195: PetscCall(DMPlexIsDistributed(dm, &flg));
1196: if (!flg) PetscFunctionReturn(PETSC_SUCCESS);
1197: /* Set initial SF so that lower level queries work */
1198: PetscCall(DMSetPointSF(dm, pointSF));
1199: PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
1200: PetscCallMPI(MPI_Comm_rank(comm, &rank));
1201: PetscCall(DMPlexGetOverlap(dm, &ov));
1202: PetscCheck(!ov, comm, PETSC_ERR_SUP, "Interpolation of overlapped DMPlex not implemented yet");
1203: PetscCall(PetscOptionsHasName(NULL, ((PetscObject)dm)->prefix, "-dmplex_interp_debug", &debug));
1204: PetscCall(PetscObjectViewFromOptions((PetscObject)dm, NULL, "-dm_interp_pre_view"));
1205: PetscCall(PetscObjectViewFromOptions((PetscObject)pointSF, NULL, "-petscsf_interp_pre_view"));
1206: PetscCall(PetscLogEventBegin(DMPLEX_InterpolateSF, dm, 0, 0, 0));
1207: /* Step 0: Precalculations */
1208: PetscCall(PetscSFGetGraph(pointSF, &Nr, &Nl, &localPoints, &remotePoints));
1209: PetscCheck(Nr >= 0, comm, PETSC_ERR_ARG_WRONGSTATE, "This DMPlex is distributed but input PointSF has no graph set");
1210: PetscCall(PetscHMapIJCreate(&remoteHash));
1211: for (l = 0; l < Nl; ++l) {
1212: PetscHashIJKey key;
1213: key.i = remotePoints[l].index;
1214: key.j = remotePoints[l].rank;
1215: PetscCall(PetscHMapIJSet(remoteHash, key, l));
1216: }
1217: /* Compute root degree to identify shared points */
1218: PetscCall(PetscSFComputeDegreeBegin(pointSF, &rootdegree));
1219: PetscCall(PetscSFComputeDegreeEnd(pointSF, &rootdegree));
1220: PetscCall(IntArrayViewFromOptions(comm, "-interp_root_degree_view", "Root degree", "point", "degree", Nr, rootdegree));
1221: /*
1222: 1) Loop over each leaf point $p$ at depth $d$ in the SF
1223: \item Get set $F(p)$ of faces $f$ in the support of $p$ for which
1224: \begin{itemize}
1225: \item all cone points of $f$ are shared
1226: \item $p$ is the cone point with smallest canonical number
1227: \end{itemize}
1228: \item Send $F(p)$ and the cone of each face to the active root point $r(p)$
1229: \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
1230: \item Send the root face from the root back to all leaf process
1231: \item Leaf processes add the shared face to the SF
1232: */
1233: /* Step 1: Construct section+SFNode array
1234: The section has entries for all shared faces for which we have a leaf point in the cone
1235: The array holds candidate shared faces, each face is referred to by the leaf point */
1236: PetscCall(PetscSectionCreate(comm, &candidateSection));
1237: PetscCall(PetscSectionSetChart(candidateSection, 0, Nr));
1238: {
1239: PetscHMapIJ faceHash;
1241: PetscCall(PetscHMapIJCreate(&faceHash));
1242: for (l = 0; l < Nl; ++l) {
1243: const PetscInt p = localPoints[l];
1245: if (debug) PetscCall(PetscSynchronizedPrintf(comm, "[%d] First pass leaf point %" PetscInt_FMT "\n", rank, p));
1246: PetscCall(DMPlexAddSharedFace_Private(dm, candidateSection, NULL, faceHash, p, debug));
1247: }
1248: PetscCall(PetscHMapIJClear(faceHash));
1249: PetscCall(PetscSectionSetUp(candidateSection));
1250: PetscCall(PetscSectionGetStorageSize(candidateSection, &candidatesSize));
1251: PetscCall(PetscMalloc1(candidatesSize, &candidates));
1252: for (l = 0; l < Nl; ++l) {
1253: const PetscInt p = localPoints[l];
1255: if (debug) PetscCall(PetscSynchronizedPrintf(comm, "[%d] Second pass leaf point %" PetscInt_FMT "\n", rank, p));
1256: PetscCall(DMPlexAddSharedFace_Private(dm, candidateSection, candidates, faceHash, p, debug));
1257: }
1258: PetscCall(PetscHMapIJDestroy(&faceHash));
1259: if (debug) PetscCall(PetscSynchronizedFlush(comm, NULL));
1260: }
1261: PetscCall(PetscObjectSetName((PetscObject)candidateSection, "Candidate Section"));
1262: PetscCall(PetscObjectViewFromOptions((PetscObject)candidateSection, NULL, "-petscsection_interp_candidate_view"));
1263: PetscCall(SFNodeArrayViewFromOptions(comm, "-petscsection_interp_candidate_view", "Candidates", NULL, candidatesSize, candidates));
1264: /* Step 2: Gather candidate section / array pair into the root partition via inverse(multi(pointSF)). */
1265: /* Note that this section is indexed by offsets into leaves, not by point number */
1266: {
1267: PetscSF sfMulti, sfInverse, sfCandidates;
1268: PetscInt *remoteOffsets;
1270: PetscCall(PetscSFGetMultiSF(pointSF, &sfMulti));
1271: PetscCall(PetscSFCreateInverseSF(sfMulti, &sfInverse));
1272: PetscCall(PetscSectionCreate(comm, &candidateRemoteSection));
1273: PetscCall(PetscSFDistributeSection(sfInverse, candidateSection, &remoteOffsets, candidateRemoteSection));
1274: PetscCall(PetscSFCreateSectionSF(sfInverse, candidateSection, remoteOffsets, candidateRemoteSection, &sfCandidates));
1275: PetscCall(PetscSectionGetStorageSize(candidateRemoteSection, &candidatesRemoteSize));
1276: PetscCall(PetscMalloc1(candidatesRemoteSize, &candidatesRemote));
1277: PetscCall(PetscSFBcastBegin(sfCandidates, MPIU_2INT, candidates, candidatesRemote, MPI_REPLACE));
1278: PetscCall(PetscSFBcastEnd(sfCandidates, MPIU_2INT, candidates, candidatesRemote, MPI_REPLACE));
1279: PetscCall(PetscSFDestroy(&sfInverse));
1280: PetscCall(PetscSFDestroy(&sfCandidates));
1281: PetscCall(PetscFree(remoteOffsets));
1283: PetscCall(PetscObjectSetName((PetscObject)candidateRemoteSection, "Remote Candidate Section"));
1284: PetscCall(PetscObjectViewFromOptions((PetscObject)candidateRemoteSection, NULL, "-petscsection_interp_candidate_remote_view"));
1285: PetscCall(SFNodeArrayViewFromOptions(comm, "-petscsection_interp_candidate_remote_view", "Remote Candidates", NULL, candidatesRemoteSize, candidatesRemote));
1286: }
1287: /* 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 */
1288: {
1289: PetscHMapIJKLRemote faceTable;
1290: PetscInt idx, idx2;
1292: PetscCall(PetscHMapIJKLRemoteCreate(&faceTable));
1293: /* There is a section point for every leaf attached to a given root point */
1294: for (r = 0, idx = 0, idx2 = 0; r < Nr; ++r) {
1295: PetscInt deg;
1297: for (deg = 0; deg < rootdegree[r]; ++deg, ++idx) {
1298: PetscInt offset, dof, d;
1300: PetscCall(PetscSectionGetDof(candidateRemoteSection, idx, &dof));
1301: PetscCall(PetscSectionGetOffset(candidateRemoteSection, idx, &offset));
1302: /* dof may include many faces from the remote process */
1303: for (d = 0; d < dof; ++d) {
1304: const PetscInt hidx = offset + d;
1305: const PetscInt Np = candidatesRemote[hidx].index + 1;
1306: const PetscSFNode rface = candidatesRemote[hidx + 1];
1307: const PetscSFNode *fcone = &candidatesRemote[hidx + 2];
1308: PetscSFNode fcp0;
1309: const PetscSFNode pmax = {PETSC_MAX_INT, PETSC_MAX_INT};
1310: const PetscInt *join = NULL;
1311: PetscHMapIJKLRemoteKey key;
1312: PetscHashIter iter;
1313: PetscBool missing, mapToLocalPointFailed = PETSC_FALSE;
1314: PetscInt points[1024], p, joinSize;
1316: if (debug)
1317: 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,
1318: rface.index, r, idx, d, Np));
1319: 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);
1320: fcp0.rank = rank;
1321: fcp0.index = r;
1322: d += Np;
1323: /* Put remote face in hash table */
1324: key.i = fcp0;
1325: key.j = fcone[0];
1326: key.k = Np > 2 ? fcone[1] : pmax;
1327: key.l = Np > 3 ? fcone[2] : pmax;
1328: PetscCall(PetscSortSFNode(Np, (PetscSFNode *)&key));
1329: PetscCall(PetscHMapIJKLRemotePut(faceTable, key, &iter, &missing));
1330: if (missing) {
1331: if (debug) PetscCall(PetscSynchronizedPrintf(PetscObjectComm((PetscObject)dm), "[%d] Setting remote face (%" PetscInt_FMT ", %" PetscInt_FMT ")\n", rank, rface.index, rface.rank));
1332: PetscCall(PetscHMapIJKLRemoteIterSet(faceTable, iter, rface));
1333: } else {
1334: PetscSFNode oface;
1336: PetscCall(PetscHMapIJKLRemoteIterGet(faceTable, iter, &oface));
1337: if ((rface.rank < oface.rank) || (rface.rank == oface.rank && rface.index < oface.index)) {
1338: if (debug) PetscCall(PetscSynchronizedPrintf(PetscObjectComm((PetscObject)dm), "[%d] Replacing with remote face (%" PetscInt_FMT ", %" PetscInt_FMT ")\n", rank, rface.index, rface.rank));
1339: PetscCall(PetscHMapIJKLRemoteIterSet(faceTable, iter, rface));
1340: }
1341: }
1342: /* Check for local face */
1343: points[0] = r;
1344: for (p = 1; p < Np; ++p) {
1345: PetscCall(DMPlexMapToLocalPoint(dm, remoteHash, fcone[p - 1], &points[p], &mapToLocalPointFailed));
1346: if (mapToLocalPointFailed) break; /* We got a point not in our overlap */
1347: if (debug) PetscCall(PetscSynchronizedPrintf(PetscObjectComm((PetscObject)dm), "[%d] Checking local candidate %" PetscInt_FMT "\n", rank, points[p]));
1348: }
1349: if (mapToLocalPointFailed) continue;
1350: PetscCall(DMPlexGetJoin(dm, Np, points, &joinSize, &join));
1351: if (joinSize == 1) {
1352: PetscSFNode lface;
1353: PetscSFNode oface;
1355: /* Always replace with local face */
1356: lface.rank = rank;
1357: lface.index = join[0];
1358: PetscCall(PetscHMapIJKLRemoteIterGet(faceTable, iter, &oface));
1359: if (debug)
1360: 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));
1361: PetscCall(PetscHMapIJKLRemoteIterSet(faceTable, iter, lface));
1362: }
1363: PetscCall(DMPlexRestoreJoin(dm, Np, points, &joinSize, &join));
1364: }
1365: }
1366: /* Put back faces for this root */
1367: for (deg = 0; deg < rootdegree[r]; ++deg, ++idx2) {
1368: PetscInt offset, dof, d;
1370: PetscCall(PetscSectionGetDof(candidateRemoteSection, idx2, &dof));
1371: PetscCall(PetscSectionGetOffset(candidateRemoteSection, idx2, &offset));
1372: /* dof may include many faces from the remote process */
1373: for (d = 0; d < dof; ++d) {
1374: const PetscInt hidx = offset + d;
1375: const PetscInt Np = candidatesRemote[hidx].index + 1;
1376: const PetscSFNode *fcone = &candidatesRemote[hidx + 2];
1377: PetscSFNode fcp0;
1378: const PetscSFNode pmax = {PETSC_MAX_INT, PETSC_MAX_INT};
1379: PetscHMapIJKLRemoteKey key;
1380: PetscHashIter iter;
1381: PetscBool missing;
1383: if (debug) PetscCall(PetscSynchronizedPrintf(PetscObjectComm((PetscObject)dm), "[%d] Entering face at (%" PetscInt_FMT ", %" PetscInt_FMT ")\n", rank, r, idx));
1384: PetscCheck(Np <= 4, PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot handle faces with %" PetscInt_FMT " cone points", Np);
1385: fcp0.rank = rank;
1386: fcp0.index = r;
1387: d += Np;
1388: /* Find remote face in hash table */
1389: key.i = fcp0;
1390: key.j = fcone[0];
1391: key.k = Np > 2 ? fcone[1] : pmax;
1392: key.l = Np > 3 ? fcone[2] : pmax;
1393: PetscCall(PetscSortSFNode(Np, (PetscSFNode *)&key));
1394: if (debug)
1395: 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,
1396: key.i.rank, key.i.index, key.j.rank, key.j.index, key.k.rank, key.k.index, key.l.rank, key.l.index));
1397: PetscCall(PetscHMapIJKLRemotePut(faceTable, key, &iter, &missing));
1398: PetscCheck(!missing, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Root %" PetscInt_FMT " Idx %" PetscInt_FMT " ought to have an associated face", r, idx2);
1399: PetscCall(PetscHMapIJKLRemoteIterGet(faceTable, iter, &candidatesRemote[hidx]));
1400: }
1401: }
1402: }
1403: if (debug) PetscCall(PetscSynchronizedFlush(PetscObjectComm((PetscObject)dm), NULL));
1404: PetscCall(PetscHMapIJKLRemoteDestroy(&faceTable));
1405: }
1406: /* Step 4: Push back owned faces */
1407: {
1408: PetscSF sfMulti, sfClaims, sfPointNew;
1409: PetscSFNode *remotePointsNew;
1410: PetscInt *remoteOffsets, *localPointsNew;
1411: PetscInt pStart, pEnd, r, NlNew, p;
1413: /* 4) Push claims back to receiver via the MultiSF and derive new pointSF mapping on receiver */
1414: PetscCall(PetscSFGetMultiSF(pointSF, &sfMulti));
1415: PetscCall(PetscSectionCreate(comm, &claimSection));
1416: PetscCall(PetscSFDistributeSection(sfMulti, candidateRemoteSection, &remoteOffsets, claimSection));
1417: PetscCall(PetscSFCreateSectionSF(sfMulti, candidateRemoteSection, remoteOffsets, claimSection, &sfClaims));
1418: PetscCall(PetscSectionGetStorageSize(claimSection, &claimsSize));
1419: PetscCall(PetscMalloc1(claimsSize, &claims));
1420: for (p = 0; p < claimsSize; ++p) claims[p].rank = -1;
1421: PetscCall(PetscSFBcastBegin(sfClaims, MPIU_2INT, candidatesRemote, claims, MPI_REPLACE));
1422: PetscCall(PetscSFBcastEnd(sfClaims, MPIU_2INT, candidatesRemote, claims, MPI_REPLACE));
1423: PetscCall(PetscSFDestroy(&sfClaims));
1424: PetscCall(PetscFree(remoteOffsets));
1425: PetscCall(PetscObjectSetName((PetscObject)claimSection, "Claim Section"));
1426: PetscCall(PetscObjectViewFromOptions((PetscObject)claimSection, NULL, "-petscsection_interp_claim_view"));
1427: PetscCall(SFNodeArrayViewFromOptions(comm, "-petscsection_interp_claim_view", "Claims", NULL, claimsSize, claims));
1428: /* Step 5) Walk the original section of local supports and add an SF entry for each updated item */
1429: /* TODO I should not have to do a join here since I already put the face and its cone in the candidate section */
1430: PetscCall(PetscHMapICreate(&claimshash));
1431: for (r = 0; r < Nr; ++r) {
1432: PetscInt dof, off, d;
1434: if (debug) PetscCall(PetscSynchronizedPrintf(comm, "[%d] Checking root for claims %" PetscInt_FMT "\n", rank, r));
1435: PetscCall(PetscSectionGetDof(candidateSection, r, &dof));
1436: PetscCall(PetscSectionGetOffset(candidateSection, r, &off));
1437: for (d = 0; d < dof;) {
1438: if (claims[off + d].rank >= 0) {
1439: const PetscInt faceInd = off + d;
1440: const PetscInt Np = candidates[off + d].index;
1441: const PetscInt *join = NULL;
1442: PetscInt joinSize, points[1024], c;
1444: if (debug) PetscCall(PetscSynchronizedPrintf(comm, "[%d] Found claim for remote point (%" PetscInt_FMT ", %" PetscInt_FMT ")\n", rank, claims[faceInd].rank, claims[faceInd].index));
1445: points[0] = r;
1446: if (debug) PetscCall(PetscSynchronizedPrintf(comm, "[%d] point %" PetscInt_FMT "\n", rank, points[0]));
1447: for (c = 0, d += 2; c < Np; ++c, ++d) {
1448: PetscCall(DMPlexMapToLocalPoint(dm, remoteHash, candidates[off + d], &points[c + 1], NULL));
1449: if (debug) PetscCall(PetscSynchronizedPrintf(comm, "[%d] point %" PetscInt_FMT "\n", rank, points[c + 1]));
1450: }
1451: PetscCall(DMPlexGetJoin(dm, Np + 1, points, &joinSize, &join));
1452: if (joinSize == 1) {
1453: if (claims[faceInd].rank == rank) {
1454: if (debug) PetscCall(PetscSynchronizedPrintf(comm, "[%d] Ignoring local face %" PetscInt_FMT " for non-remote partner\n", rank, join[0]));
1455: } else {
1456: if (debug) PetscCall(PetscSynchronizedPrintf(comm, "[%d] Found local face %" PetscInt_FMT "\n", rank, join[0]));
1457: PetscCall(PetscHMapISet(claimshash, join[0], faceInd));
1458: }
1459: } else {
1460: if (debug) PetscCall(PetscSynchronizedPrintf(comm, "[%d] Failed to find face\n", rank));
1461: }
1462: PetscCall(DMPlexRestoreJoin(dm, Np + 1, points, &joinSize, &join));
1463: } else {
1464: if (debug) PetscCall(PetscSynchronizedPrintf(comm, "[%d] No claim for point %" PetscInt_FMT "\n", rank, r));
1465: d += claims[off + d].index + 1;
1466: }
1467: }
1468: }
1469: if (debug) PetscCall(PetscSynchronizedFlush(comm, NULL));
1470: /* Step 6) Create new pointSF from hashed claims */
1471: PetscCall(PetscHMapIGetSize(claimshash, &NlNew));
1472: PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
1473: PetscCall(PetscMalloc1(Nl + NlNew, &localPointsNew));
1474: PetscCall(PetscMalloc1(Nl + NlNew, &remotePointsNew));
1475: for (l = 0; l < Nl; ++l) {
1476: localPointsNew[l] = localPoints[l];
1477: remotePointsNew[l].index = remotePoints[l].index;
1478: remotePointsNew[l].rank = remotePoints[l].rank;
1479: }
1480: p = Nl;
1481: PetscCall(PetscHMapIGetKeys(claimshash, &p, localPointsNew));
1482: /* We sort new points, and assume they are numbered after all existing points */
1483: PetscCall(PetscSortInt(NlNew, PetscSafePointerPlusOffset(localPointsNew, Nl)));
1484: for (p = Nl; p < Nl + NlNew; ++p) {
1485: PetscInt off;
1486: PetscCall(PetscHMapIGet(claimshash, localPointsNew[p], &off));
1487: 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);
1488: remotePointsNew[p] = claims[off];
1489: }
1490: PetscCall(PetscSFCreate(comm, &sfPointNew));
1491: PetscCall(PetscSFSetGraph(sfPointNew, pEnd - pStart, Nl + NlNew, localPointsNew, PETSC_OWN_POINTER, remotePointsNew, PETSC_OWN_POINTER));
1492: PetscCall(PetscSFSetUp(sfPointNew));
1493: PetscCall(DMSetPointSF(dm, sfPointNew));
1494: PetscCall(PetscObjectViewFromOptions((PetscObject)sfPointNew, NULL, "-petscsf_interp_view"));
1495: if (PetscDefined(USE_DEBUG)) PetscCall(DMPlexCheckPointSF(dm, sfPointNew, PETSC_FALSE));
1496: PetscCall(PetscSFDestroy(&sfPointNew));
1497: PetscCall(PetscHMapIDestroy(&claimshash));
1498: }
1499: PetscCall(PetscHMapIJDestroy(&remoteHash));
1500: PetscCall(PetscSectionDestroy(&candidateSection));
1501: PetscCall(PetscSectionDestroy(&candidateRemoteSection));
1502: PetscCall(PetscSectionDestroy(&claimSection));
1503: PetscCall(PetscFree(candidates));
1504: PetscCall(PetscFree(candidatesRemote));
1505: PetscCall(PetscFree(claims));
1506: PetscCall(PetscLogEventEnd(DMPLEX_InterpolateSF, dm, 0, 0, 0));
1507: PetscFunctionReturn(PETSC_SUCCESS);
1508: }
1510: /*@
1511: DMPlexInterpolate - Take in a cell-vertex mesh and return one with all intermediate faces, edges, etc.
1513: Collective
1515: Input Parameter:
1516: . dm - The `DMPLEX` object with only cells and vertices
1518: Output Parameter:
1519: . dmInt - The complete `DMPLEX` object
1521: Level: intermediate
1523: Note:
1524: Labels and coordinates are copied.
1526: Developer Notes:
1527: It sets plex->interpolated = `DMPLEX_INTERPOLATED_FULL`.
1529: .seealso: `DMPLEX`, `DMPlexUninterpolate()`, `DMPlexCreateFromCellListPetsc()`, `DMPlexCopyCoordinates()`
1530: @*/
1531: PetscErrorCode DMPlexInterpolate(DM dm, DM *dmInt)
1532: {
1533: DMPlexInterpolatedFlag interpolated;
1534: DM idm, odm = dm;
1535: PetscSF sfPoint;
1536: PetscInt depth, dim, d;
1537: const char *name;
1538: PetscBool flg = PETSC_TRUE;
1540: PetscFunctionBegin;
1542: PetscAssertPointer(dmInt, 2);
1543: PetscCall(PetscLogEventBegin(DMPLEX_Interpolate, dm, 0, 0, 0));
1544: PetscCall(DMPlexGetDepth(dm, &depth));
1545: PetscCall(DMGetDimension(dm, &dim));
1546: PetscCall(DMPlexIsInterpolated(dm, &interpolated));
1547: PetscCheck(interpolated != DMPLEX_INTERPOLATED_PARTIAL, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Not for partially interpolated meshes");
1548: if (interpolated == DMPLEX_INTERPOLATED_FULL) {
1549: PetscCall(PetscObjectReference((PetscObject)dm));
1550: idm = dm;
1551: } else {
1552: PetscBool nonmanifold = PETSC_FALSE;
1554: PetscCall(PetscOptionsGetBool(NULL, dm->hdr.prefix, "-dm_plex_stratify_celltype", &nonmanifold, NULL));
1555: if (nonmanifold) {
1556: do {
1557: const char *prefix;
1558: PetscInt pStart, pEnd, pdepth;
1559: PetscBool done = PETSC_TRUE;
1561: // Find a point which is not correctly interpolated
1562: PetscCall(DMPlexGetChart(odm, &pStart, &pEnd));
1563: for (PetscInt p = pStart; p < pEnd; ++p) {
1564: DMPolytopeType ct;
1565: const PetscInt *cone;
1566: PetscInt coneSize, cdepth;
1568: PetscCall(DMPlexGetPointDepth(odm, p, &pdepth));
1569: PetscCall(DMPlexGetCellType(odm, p, &ct));
1570: // Check against celltype
1571: if (pdepth != DMPolytopeTypeGetDim(ct)) {
1572: done = PETSC_FALSE;
1573: break;
1574: }
1575: // Check against boundary
1576: PetscCall(DMPlexGetCone(odm, p, &cone));
1577: PetscCall(DMPlexGetConeSize(odm, p, &coneSize));
1578: for (PetscInt c = 0; c < coneSize; ++c) {
1579: PetscCall(DMPlexGetPointDepth(odm, cone[c], &cdepth));
1580: if (cdepth != pdepth - 1) {
1581: done = PETSC_FALSE;
1582: p = pEnd;
1583: break;
1584: }
1585: }
1586: }
1587: if (done) break;
1588: /* Create interpolated mesh */
1589: PetscCall(DMCreate(PetscObjectComm((PetscObject)dm), &idm));
1590: PetscCall(DMSetType(idm, DMPLEX));
1591: PetscCall(DMSetDimension(idm, dim));
1592: PetscCall(PetscObjectGetOptionsPrefix((PetscObject)dm, &prefix));
1593: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)idm, prefix));
1594: if (depth > 0) {
1595: PetscCall(DMPlexInterpolateFaces_Internal(odm, pdepth, idm));
1596: PetscCall(DMGetPointSF(odm, &sfPoint));
1597: if (PetscDefined(USE_DEBUG)) PetscCall(DMPlexCheckPointSF(odm, sfPoint, PETSC_FALSE));
1598: {
1599: /* TODO: We need to systematically fix cases of distributed Plexes with no graph set */
1600: PetscInt nroots;
1601: PetscCall(PetscSFGetGraph(sfPoint, &nroots, NULL, NULL, NULL));
1602: if (nroots >= 0) PetscCall(DMPlexInterpolatePointSF(idm, sfPoint));
1603: }
1604: }
1605: if (odm != dm) PetscCall(DMDestroy(&odm));
1606: odm = idm;
1607: } while (1);
1608: } else {
1609: for (d = 1; d < dim; ++d) {
1610: const char *prefix;
1612: /* Create interpolated mesh */
1613: PetscCall(DMCreate(PetscObjectComm((PetscObject)dm), &idm));
1614: PetscCall(DMSetType(idm, DMPLEX));
1615: PetscCall(DMSetDimension(idm, dim));
1616: PetscCall(PetscObjectGetOptionsPrefix((PetscObject)dm, &prefix));
1617: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)idm, prefix));
1618: if (depth > 0) {
1619: PetscCall(DMPlexInterpolateFaces_Internal(odm, 1, idm));
1620: PetscCall(DMGetPointSF(odm, &sfPoint));
1621: if (PetscDefined(USE_DEBUG)) PetscCall(DMPlexCheckPointSF(odm, sfPoint, PETSC_FALSE));
1622: {
1623: /* TODO: We need to systematically fix cases of distributed Plexes with no graph set */
1624: PetscInt nroots;
1625: PetscCall(PetscSFGetGraph(sfPoint, &nroots, NULL, NULL, NULL));
1626: if (nroots >= 0) PetscCall(DMPlexInterpolatePointSF(idm, sfPoint));
1627: }
1628: }
1629: if (odm != dm) PetscCall(DMDestroy(&odm));
1630: odm = idm;
1631: }
1632: }
1633: PetscCall(PetscObjectGetName((PetscObject)dm, &name));
1634: PetscCall(PetscObjectSetName((PetscObject)idm, name));
1635: PetscCall(DMPlexCopyCoordinates(dm, idm));
1636: PetscCall(DMCopyLabels(dm, idm, PETSC_COPY_VALUES, PETSC_FALSE, DM_COPY_LABELS_FAIL));
1637: PetscCall(PetscOptionsGetBool(((PetscObject)dm)->options, ((PetscObject)dm)->prefix, "-dm_plex_interpolate_orient_interfaces", &flg, NULL));
1638: if (flg) PetscCall(DMPlexOrientInterface_Internal(idm));
1639: }
1640: /* This function makes the mesh fully interpolated on all ranks */
1641: {
1642: DM_Plex *plex = (DM_Plex *)idm->data;
1643: plex->interpolated = plex->interpolatedCollective = DMPLEX_INTERPOLATED_FULL;
1644: }
1645: PetscCall(DMPlexCopy_Internal(dm, PETSC_TRUE, PETSC_TRUE, idm));
1646: *dmInt = idm;
1647: PetscCall(PetscLogEventEnd(DMPLEX_Interpolate, dm, 0, 0, 0));
1648: PetscFunctionReturn(PETSC_SUCCESS);
1649: }
1651: /*@
1652: DMPlexCopyCoordinates - Copy coordinates from one mesh to another with the same vertices
1654: Collective
1656: Input Parameter:
1657: . dmA - The `DMPLEX` object with initial coordinates
1659: Output Parameter:
1660: . dmB - The `DMPLEX` object with copied coordinates
1662: Level: intermediate
1664: Notes:
1665: This is typically used when adding pieces other than vertices to a mesh
1667: This function does not copy localized coordinates.
1669: .seealso: `DMPLEX`, `DMCopyLabels()`, `DMGetCoordinates()`, `DMGetCoordinatesLocal()`, `DMGetCoordinateDM()`, `DMGetCoordinateSection()`
1670: @*/
1671: PetscErrorCode DMPlexCopyCoordinates(DM dmA, DM dmB)
1672: {
1673: Vec coordinatesA, coordinatesB;
1674: VecType vtype;
1675: PetscSection coordSectionA, coordSectionB;
1676: PetscScalar *coordsA, *coordsB;
1677: PetscInt spaceDim, Nf, vStartA, vStartB, vEndA, vEndB, coordSizeB, v, d;
1678: PetscInt cStartA, cEndA, cStartB, cEndB, cS, cE, cdim;
1680: PetscFunctionBegin;
1683: if (dmA == dmB) PetscFunctionReturn(PETSC_SUCCESS);
1684: PetscCall(DMGetCoordinateDim(dmA, &cdim));
1685: PetscCall(DMSetCoordinateDim(dmB, cdim));
1686: PetscCall(DMPlexGetDepthStratum(dmA, 0, &vStartA, &vEndA));
1687: PetscCall(DMPlexGetDepthStratum(dmB, 0, &vStartB, &vEndB));
1688: 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);
1689: /* Copy over discretization if it exists */
1690: {
1691: DM cdmA, cdmB;
1692: PetscDS dsA, dsB;
1693: PetscObject objA, objB;
1694: PetscClassId idA, idB;
1695: const PetscScalar *constants;
1696: PetscInt cdim, Nc;
1698: PetscCall(DMGetCoordinateDM(dmA, &cdmA));
1699: PetscCall(DMGetCoordinateDM(dmB, &cdmB));
1700: PetscCall(DMGetField(cdmA, 0, NULL, &objA));
1701: PetscCall(DMGetField(cdmB, 0, NULL, &objB));
1702: PetscCall(PetscObjectGetClassId(objA, &idA));
1703: PetscCall(PetscObjectGetClassId(objB, &idB));
1704: if ((idA == PETSCFE_CLASSID) && (idA != idB)) {
1705: PetscCall(DMSetField(cdmB, 0, NULL, objA));
1706: PetscCall(DMCreateDS(cdmB));
1707: PetscCall(DMGetDS(cdmA, &dsA));
1708: PetscCall(DMGetDS(cdmB, &dsB));
1709: PetscCall(PetscDSGetCoordinateDimension(dsA, &cdim));
1710: PetscCall(PetscDSSetCoordinateDimension(dsB, cdim));
1711: PetscCall(PetscDSGetConstants(dsA, &Nc, &constants));
1712: PetscCall(PetscDSSetConstants(dsB, Nc, (PetscScalar *)constants));
1713: }
1714: }
1715: PetscCall(DMPlexGetHeightStratum(dmA, 0, &cStartA, &cEndA));
1716: PetscCall(DMPlexGetHeightStratum(dmB, 0, &cStartB, &cEndB));
1717: PetscCall(DMGetCoordinateSection(dmA, &coordSectionA));
1718: PetscCall(DMGetCoordinateSection(dmB, &coordSectionB));
1719: if (coordSectionA == coordSectionB) PetscFunctionReturn(PETSC_SUCCESS);
1720: PetscCall(PetscSectionGetNumFields(coordSectionA, &Nf));
1721: if (!Nf) PetscFunctionReturn(PETSC_SUCCESS);
1722: PetscCheck(Nf <= 1, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "The number of coordinate fields must be 1, not %" PetscInt_FMT, Nf);
1723: if (!coordSectionB) {
1724: PetscInt dim;
1726: PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)coordSectionA), &coordSectionB));
1727: PetscCall(DMGetCoordinateDim(dmA, &dim));
1728: PetscCall(DMSetCoordinateSection(dmB, dim, coordSectionB));
1729: PetscCall(PetscObjectDereference((PetscObject)coordSectionB));
1730: }
1731: PetscCall(PetscSectionSetNumFields(coordSectionB, 1));
1732: PetscCall(PetscSectionGetFieldComponents(coordSectionA, 0, &spaceDim));
1733: PetscCall(PetscSectionSetFieldComponents(coordSectionB, 0, spaceDim));
1734: PetscCall(PetscSectionGetChart(coordSectionA, &cS, &cE));
1735: cS = vStartB;
1736: cE = vEndB;
1737: PetscCall(PetscSectionSetChart(coordSectionB, cS, cE));
1738: for (v = vStartB; v < vEndB; ++v) {
1739: PetscCall(PetscSectionSetDof(coordSectionB, v, spaceDim));
1740: PetscCall(PetscSectionSetFieldDof(coordSectionB, v, 0, spaceDim));
1741: }
1742: PetscCall(PetscSectionSetUp(coordSectionB));
1743: PetscCall(PetscSectionGetStorageSize(coordSectionB, &coordSizeB));
1744: PetscCall(DMGetCoordinatesLocal(dmA, &coordinatesA));
1745: PetscCall(VecCreate(PETSC_COMM_SELF, &coordinatesB));
1746: PetscCall(PetscObjectSetName((PetscObject)coordinatesB, "coordinates"));
1747: PetscCall(VecSetSizes(coordinatesB, coordSizeB, PETSC_DETERMINE));
1748: PetscCall(VecGetBlockSize(coordinatesA, &d));
1749: PetscCall(VecSetBlockSize(coordinatesB, d));
1750: PetscCall(VecGetType(coordinatesA, &vtype));
1751: PetscCall(VecSetType(coordinatesB, vtype));
1752: PetscCall(VecGetArray(coordinatesA, &coordsA));
1753: PetscCall(VecGetArray(coordinatesB, &coordsB));
1754: for (v = 0; v < vEndB - vStartB; ++v) {
1755: PetscInt offA, offB;
1757: PetscCall(PetscSectionGetOffset(coordSectionA, v + vStartA, &offA));
1758: PetscCall(PetscSectionGetOffset(coordSectionB, v + vStartB, &offB));
1759: for (d = 0; d < spaceDim; ++d) coordsB[offB + d] = coordsA[offA + d];
1760: }
1761: PetscCall(VecRestoreArray(coordinatesA, &coordsA));
1762: PetscCall(VecRestoreArray(coordinatesB, &coordsB));
1763: PetscCall(DMSetCoordinatesLocal(dmB, coordinatesB));
1764: PetscCall(VecDestroy(&coordinatesB));
1765: PetscFunctionReturn(PETSC_SUCCESS);
1766: }
1768: /*@
1769: DMPlexUninterpolate - Take in a mesh with all intermediate faces, edges, etc. and return a cell-vertex mesh
1771: Collective
1773: Input Parameter:
1774: . dm - The complete `DMPLEX` object
1776: Output Parameter:
1777: . dmUnint - The `DMPLEX` object with only cells and vertices
1779: Level: intermediate
1781: Note:
1782: It does not copy over the coordinates.
1784: Developer Notes:
1785: Sets plex->interpolated = `DMPLEX_INTERPOLATED_NONE`.
1787: .seealso: `DMPLEX`, `DMPlexInterpolate()`, `DMPlexCreateFromCellListPetsc()`, `DMPlexCopyCoordinates()`
1788: @*/
1789: PetscErrorCode DMPlexUninterpolate(DM dm, DM *dmUnint)
1790: {
1791: DMPlexInterpolatedFlag interpolated;
1792: DM udm;
1793: PetscInt dim, vStart, vEnd, cStart, cEnd, c, maxConeSize = 0, *cone;
1795: PetscFunctionBegin;
1797: PetscAssertPointer(dmUnint, 2);
1798: PetscCall(PetscLogEventBegin(DMPLEX_Uninterpolate, dm, 0, 0, 0));
1799: PetscCall(DMGetDimension(dm, &dim));
1800: PetscCall(DMPlexIsInterpolated(dm, &interpolated));
1801: PetscCheck(interpolated != DMPLEX_INTERPOLATED_PARTIAL, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Not for partially interpolated meshes");
1802: if (interpolated == DMPLEX_INTERPOLATED_NONE || dim <= 1) {
1803: /* in case dim <= 1 just keep the DMPLEX_INTERPOLATED_FULL flag */
1804: PetscCall(PetscObjectReference((PetscObject)dm));
1805: *dmUnint = dm;
1806: PetscFunctionReturn(PETSC_SUCCESS);
1807: }
1808: PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd));
1809: PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
1810: PetscCall(DMCreate(PetscObjectComm((PetscObject)dm), &udm));
1811: PetscCall(DMSetType(udm, DMPLEX));
1812: PetscCall(DMSetDimension(udm, dim));
1813: PetscCall(DMPlexSetChart(udm, cStart, vEnd));
1814: for (c = cStart; c < cEnd; ++c) {
1815: PetscInt *closure = NULL, closureSize, cl, coneSize = 0;
1817: PetscCall(DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure));
1818: for (cl = 0; cl < closureSize * 2; cl += 2) {
1819: const PetscInt p = closure[cl];
1821: if ((p >= vStart) && (p < vEnd)) ++coneSize;
1822: }
1823: PetscCall(DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure));
1824: PetscCall(DMPlexSetConeSize(udm, c, coneSize));
1825: maxConeSize = PetscMax(maxConeSize, coneSize);
1826: }
1827: PetscCall(DMSetUp(udm));
1828: PetscCall(PetscMalloc1(maxConeSize, &cone));
1829: for (c = cStart; c < cEnd; ++c) {
1830: PetscInt *closure = NULL, closureSize, cl, coneSize = 0;
1832: PetscCall(DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure));
1833: for (cl = 0; cl < closureSize * 2; cl += 2) {
1834: const PetscInt p = closure[cl];
1836: if ((p >= vStart) && (p < vEnd)) cone[coneSize++] = p;
1837: }
1838: PetscCall(DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure));
1839: PetscCall(DMPlexSetCone(udm, c, cone));
1840: }
1841: PetscCall(PetscFree(cone));
1842: PetscCall(DMPlexSymmetrize(udm));
1843: PetscCall(DMPlexStratify(udm));
1844: /* Reduce SF */
1845: {
1846: PetscSF sfPoint, sfPointUn;
1847: const PetscSFNode *remotePoints;
1848: const PetscInt *localPoints;
1849: PetscSFNode *remotePointsUn;
1850: PetscInt *localPointsUn;
1851: PetscInt numRoots, numLeaves, l;
1852: PetscInt numLeavesUn = 0, n = 0;
1854: /* Get original SF information */
1855: PetscCall(DMGetPointSF(dm, &sfPoint));
1856: if (PetscDefined(USE_DEBUG)) PetscCall(DMPlexCheckPointSF(dm, sfPoint, PETSC_FALSE));
1857: PetscCall(DMGetPointSF(udm, &sfPointUn));
1858: PetscCall(PetscSFGetGraph(sfPoint, &numRoots, &numLeaves, &localPoints, &remotePoints));
1859: if (numRoots >= 0) {
1860: /* Allocate space for cells and vertices */
1861: for (l = 0; l < numLeaves; ++l) {
1862: const PetscInt p = localPoints[l];
1864: if ((vStart <= p && p < vEnd) || (cStart <= p && p < cEnd)) numLeavesUn++;
1865: }
1866: /* Fill in leaves */
1867: PetscCall(PetscMalloc1(numLeavesUn, &remotePointsUn));
1868: PetscCall(PetscMalloc1(numLeavesUn, &localPointsUn));
1869: for (l = 0; l < numLeaves; l++) {
1870: const PetscInt p = localPoints[l];
1872: if ((vStart <= p && p < vEnd) || (cStart <= p && p < cEnd)) {
1873: localPointsUn[n] = p;
1874: remotePointsUn[n].rank = remotePoints[l].rank;
1875: remotePointsUn[n].index = remotePoints[l].index;
1876: ++n;
1877: }
1878: }
1879: PetscCheck(n == numLeavesUn, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Inconsistent number of leaves %" PetscInt_FMT " != %" PetscInt_FMT, n, numLeavesUn);
1880: PetscCall(PetscSFSetGraph(sfPointUn, cEnd - cStart + vEnd - vStart, numLeavesUn, localPointsUn, PETSC_OWN_POINTER, remotePointsUn, PETSC_OWN_POINTER));
1881: }
1882: }
1883: /* This function makes the mesh fully uninterpolated on all ranks */
1884: {
1885: DM_Plex *plex = (DM_Plex *)udm->data;
1886: plex->interpolated = plex->interpolatedCollective = DMPLEX_INTERPOLATED_NONE;
1887: }
1888: PetscCall(DMPlexCopy_Internal(dm, PETSC_TRUE, PETSC_TRUE, udm));
1889: if (PetscDefined(USE_DEBUG)) PetscCall(DMPlexCheckPointSF(udm, NULL, PETSC_FALSE));
1890: *dmUnint = udm;
1891: PetscCall(PetscLogEventEnd(DMPLEX_Uninterpolate, dm, 0, 0, 0));
1892: PetscFunctionReturn(PETSC_SUCCESS);
1893: }
1895: static PetscErrorCode DMPlexIsInterpolated_Internal(DM dm, DMPlexInterpolatedFlag *interpolated)
1896: {
1897: PetscInt coneSize, depth, dim, h, p, pStart, pEnd;
1898: MPI_Comm comm;
1900: PetscFunctionBegin;
1901: PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
1902: PetscCall(DMPlexGetDepth(dm, &depth));
1903: PetscCall(DMGetDimension(dm, &dim));
1905: if (depth == dim) {
1906: *interpolated = DMPLEX_INTERPOLATED_FULL;
1907: if (!dim) goto finish;
1909: /* Check points at height = dim are vertices (have no cones) */
1910: PetscCall(DMPlexGetHeightStratum(dm, dim, &pStart, &pEnd));
1911: for (p = pStart; p < pEnd; p++) {
1912: PetscCall(DMPlexGetConeSize(dm, p, &coneSize));
1913: if (coneSize) {
1914: *interpolated = DMPLEX_INTERPOLATED_PARTIAL;
1915: goto finish;
1916: }
1917: }
1919: /* Check points at height < dim have cones */
1920: for (h = 0; h < dim; h++) {
1921: PetscCall(DMPlexGetHeightStratum(dm, h, &pStart, &pEnd));
1922: for (p = pStart; p < pEnd; p++) {
1923: PetscCall(DMPlexGetConeSize(dm, p, &coneSize));
1924: if (!coneSize) {
1925: *interpolated = DMPLEX_INTERPOLATED_PARTIAL;
1926: goto finish;
1927: }
1928: }
1929: }
1930: } else if (depth == 1) {
1931: *interpolated = DMPLEX_INTERPOLATED_NONE;
1932: } else {
1933: *interpolated = DMPLEX_INTERPOLATED_PARTIAL;
1934: }
1935: finish:
1936: PetscFunctionReturn(PETSC_SUCCESS);
1937: }
1939: /*@
1940: DMPlexIsInterpolated - Find out to what extent the `DMPLEX` is topologically interpolated.
1942: Not Collective
1944: Input Parameter:
1945: . dm - The `DMPLEX` object
1947: Output Parameter:
1948: . interpolated - Flag whether the `DM` is interpolated
1950: Level: intermediate
1952: Notes:
1953: Unlike `DMPlexIsInterpolatedCollective()`, this is NOT collective
1954: so the results can be different on different ranks in special cases.
1955: However, `DMPlexInterpolate()` guarantees the result is the same on all.
1957: Unlike `DMPlexIsInterpolatedCollective()`, this cannot return `DMPLEX_INTERPOLATED_MIXED`.
1959: Developer Notes:
1960: Initially, plex->interpolated = `DMPLEX_INTERPOLATED_INVALID`.
1962: If plex->interpolated == `DMPLEX_INTERPOLATED_INVALID`, `DMPlexIsInterpolated_Internal()` is called.
1963: It checks the actual topology and sets plex->interpolated on each rank separately to one of
1964: `DMPLEX_INTERPOLATED_NONE`, `DMPLEX_INTERPOLATED_PARTIAL` or `DMPLEX_INTERPOLATED_FULL`.
1966: If plex->interpolated != `DMPLEX_INTERPOLATED_INVALID`, this function just returns plex->interpolated.
1968: `DMPlexInterpolate()` sets plex->interpolated = `DMPLEX_INTERPOLATED_FULL`,
1969: and DMPlexUninterpolate() sets plex->interpolated = `DMPLEX_INTERPOLATED_NONE`.
1971: .seealso: `DMPLEX`, `DMPlexInterpolate()`, `DMPlexIsInterpolatedCollective()`
1972: @*/
1973: PetscErrorCode DMPlexIsInterpolated(DM dm, DMPlexInterpolatedFlag *interpolated)
1974: {
1975: DM_Plex *plex = (DM_Plex *)dm->data;
1977: PetscFunctionBegin;
1979: PetscAssertPointer(interpolated, 2);
1980: if (plex->interpolated < 0) {
1981: PetscCall(DMPlexIsInterpolated_Internal(dm, &plex->interpolated));
1982: } else if (PetscDefined(USE_DEBUG)) {
1983: DMPlexInterpolatedFlag flg;
1985: PetscCall(DMPlexIsInterpolated_Internal(dm, &flg));
1986: 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]);
1987: }
1988: *interpolated = plex->interpolated;
1989: PetscFunctionReturn(PETSC_SUCCESS);
1990: }
1992: /*@
1993: DMPlexIsInterpolatedCollective - Find out to what extent the `DMPLEX` is topologically interpolated (in collective manner).
1995: Collective
1997: Input Parameter:
1998: . dm - The `DMPLEX` object
2000: Output Parameter:
2001: . interpolated - Flag whether the `DM` is interpolated
2003: Level: intermediate
2005: Notes:
2006: Unlike `DMPlexIsInterpolated()`, this is collective so the results are guaranteed to be the same on all ranks.
2008: This function will return `DMPLEX_INTERPOLATED_MIXED` if the results of `DMPlexIsInterpolated()` are different on different ranks.
2010: Developer Notes:
2011: Initially, plex->interpolatedCollective = `DMPLEX_INTERPOLATED_INVALID`.
2013: If plex->interpolatedCollective == `DMPLEX_INTERPOLATED_INVALID`, this function calls `DMPlexIsInterpolated()` which sets plex->interpolated.
2014: `MPI_Allreduce()` is then called and collectively consistent flag plex->interpolatedCollective is set and returned;
2015: if plex->interpolated varies on different ranks, plex->interpolatedCollective = `DMPLEX_INTERPOLATED_MIXED`,
2016: otherwise sets plex->interpolatedCollective = plex->interpolated.
2018: If plex->interpolatedCollective != `DMPLEX_INTERPOLATED_INVALID`, this function just returns plex->interpolatedCollective.
2020: .seealso: `DMPLEX`, `DMPlexInterpolate()`, `DMPlexIsInterpolated()`
2021: @*/
2022: PetscErrorCode DMPlexIsInterpolatedCollective(DM dm, DMPlexInterpolatedFlag *interpolated)
2023: {
2024: DM_Plex *plex = (DM_Plex *)dm->data;
2025: PetscBool debug = PETSC_FALSE;
2027: PetscFunctionBegin;
2029: PetscAssertPointer(interpolated, 2);
2030: PetscCall(PetscOptionsGetBool(((PetscObject)dm)->options, ((PetscObject)dm)->prefix, "-dm_plex_is_interpolated_collective_debug", &debug, NULL));
2031: if (plex->interpolatedCollective < 0) {
2032: DMPlexInterpolatedFlag min, max;
2033: MPI_Comm comm;
2035: PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
2036: PetscCall(DMPlexIsInterpolated(dm, &plex->interpolatedCollective));
2037: PetscCall(MPIU_Allreduce(&plex->interpolatedCollective, &min, 1, MPIU_ENUM, MPI_MIN, comm));
2038: PetscCall(MPIU_Allreduce(&plex->interpolatedCollective, &max, 1, MPIU_ENUM, MPI_MAX, comm));
2039: if (min != max) plex->interpolatedCollective = DMPLEX_INTERPOLATED_MIXED;
2040: if (debug) {
2041: PetscMPIInt rank;
2043: PetscCallMPI(MPI_Comm_rank(comm, &rank));
2044: PetscCall(PetscSynchronizedPrintf(comm, "[%d] interpolated=%s interpolatedCollective=%s\n", rank, DMPlexInterpolatedFlags[plex->interpolated], DMPlexInterpolatedFlags[plex->interpolatedCollective]));
2045: PetscCall(PetscSynchronizedFlush(comm, PETSC_STDOUT));
2046: }
2047: }
2048: *interpolated = plex->interpolatedCollective;
2049: PetscFunctionReturn(PETSC_SUCCESS);
2050: }