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_INT_MAX;
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_INT_MAX
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: PetscMPIInt nranks;
778: const PetscMPIInt *ranks = NULL;
779: const PetscInt *roffset = NULL, *rmine = NULL, *rremote = NULL;
780: PetscInt n, o;
782: PetscFunctionBegin;
783: PetscCall(PetscSFGetRootRanks(sf, &nranks, &ranks, &roffset, &rmine, &rremote));
784: nleaves = roffset[nranks];
785: PetscCall(PetscMalloc2(nleaves, rmine1, nleaves, rremote1));
786: for (PetscMPIInt 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, maxConeSize = 0;
807: PetscMPIInt nranks, r;
808: PetscInt(*roots)[4], (*leaves)[4], mainCone[4];
809: PetscMPIInt(*rootsRanks)[4], (*leavesRanks)[4];
810: MPI_Comm comm;
811: PetscMPIInt rank, size;
812: PetscInt debug = 0;
814: PetscFunctionBegin;
815: PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
816: PetscCallMPI(MPI_Comm_rank(comm, &rank));
817: PetscCallMPI(MPI_Comm_size(comm, &size));
818: PetscCall(DMGetPointSF(dm, &sf));
819: PetscCall(DMViewFromOptions(dm, NULL, "-before_orient_interface_dm_view"));
820: if (PetscDefined(USE_DEBUG)) PetscCall(DMPlexCheckPointSF(dm, sf, PETSC_FALSE));
821: PetscCall(PetscSFGetGraph(sf, &nroots, &nleaves, &locals, &remotes));
822: if (nroots < 0) PetscFunctionReturn(PETSC_SUCCESS);
823: PetscCall(PetscSFSetUp(sf));
824: PetscCall(SortRmineRremoteByRemote_Private(sf, &rmine1, &rremote1));
825: for (p = 0; p < nleaves; ++p) {
826: PetscInt coneSize;
827: PetscCall(DMPlexGetConeSize(dm, locals[p], &coneSize));
828: maxConeSize = PetscMax(maxConeSize, coneSize);
829: }
830: PetscCheck(maxConeSize <= 4, comm, PETSC_ERR_SUP, "This method does not support cones of size %" PetscInt_FMT, maxConeSize);
831: PetscCall(PetscMalloc4(nroots, &roots, nroots, &leaves, nroots, &rootsRanks, nroots, &leavesRanks));
832: for (p = 0; p < nroots; ++p) {
833: const PetscInt *cone;
834: PetscInt coneSize, c, ind0;
836: PetscCall(DMPlexGetConeSize(dm, p, &coneSize));
837: PetscCall(DMPlexGetCone(dm, p, &cone));
838: /* Ignore vertices */
839: if (coneSize < 2) {
840: for (c = 0; c < 4; c++) {
841: roots[p][c] = -1;
842: rootsRanks[p][c] = -1;
843: }
844: continue;
845: }
846: /* Translate all points to root numbering */
847: for (c = 0; c < PetscMin(coneSize, 4); c++) {
848: PetscCall(PetscFindInt(cone[c], nleaves, locals, &ind0));
849: if (ind0 < 0) {
850: roots[p][c] = cone[c];
851: rootsRanks[p][c] = rank;
852: } else {
853: roots[p][c] = remotes[ind0].index;
854: rootsRanks[p][c] = (PetscMPIInt)remotes[ind0].rank;
855: }
856: }
857: for (c = coneSize; c < 4; c++) {
858: roots[p][c] = -1;
859: rootsRanks[p][c] = -1;
860: }
861: }
862: for (p = 0; p < nroots; ++p) {
863: PetscInt c;
864: for (c = 0; c < 4; c++) {
865: leaves[p][c] = -2;
866: leavesRanks[p][c] = -2;
867: }
868: }
869: PetscCall(PetscSFBcastBegin(sf, MPIU_4INT, roots, leaves, MPI_REPLACE));
870: PetscCall(PetscSFBcastBegin(sf, MPI_4INT, rootsRanks, leavesRanks, MPI_REPLACE));
871: PetscCall(PetscSFBcastEnd(sf, MPIU_4INT, roots, leaves, MPI_REPLACE));
872: PetscCall(PetscSFBcastEnd(sf, MPI_4INT, rootsRanks, leavesRanks, MPI_REPLACE));
873: if (debug) {
874: PetscCall(PetscSynchronizedFlush(comm, NULL));
875: if (rank == 0) PetscCall(PetscSynchronizedPrintf(comm, "Referenced roots\n"));
876: }
877: PetscCall(PetscSFGetRootRanks(sf, &nranks, &ranks, &roffset, NULL, NULL));
878: for (p = 0; p < nroots; ++p) {
879: DMPolytopeType ct;
880: const PetscInt *cone;
881: PetscInt coneSize, c, ind0, o, ir;
883: if (leaves[p][0] < 0) continue; /* Ignore vertices */
884: PetscCall(DMPlexGetCellType(dm, p, &ct));
885: PetscCall(DMPlexGetConeSize(dm, p, &coneSize));
886: PetscCall(DMPlexGetCone(dm, p, &cone));
887: if (debug) {
888: 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]));
889: }
890: 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]) {
891: /* Translate these leaves to my cone points; mainCone means desired order p's cone points */
892: for (c = 0; c < PetscMin(coneSize, 4); ++c) {
893: PetscInt rS, rN;
895: if (leavesRanks[p][c] == rank) {
896: /* A local leaf is just taken as it is */
897: mainCone[c] = leaves[p][c];
898: continue;
899: }
900: /* Find index of rank leavesRanks[p][c] among remote ranks */
901: /* No need for PetscMPIIntCast because these integers were originally cast from PetscMPIInt. */
902: PetscCall(PetscFindMPIInt((PetscMPIInt)leavesRanks[p][c], nranks, ranks, &ir));
903: PetscCall(PetscMPIIntCast(ir, &r));
904: 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]);
905: 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]);
906: /* Find point leaves[p][c] among remote points aimed at rank leavesRanks[p][c] */
907: rS = roffset[r];
908: rN = roffset[r + 1] - rS;
909: PetscCall(PetscFindInt(leaves[p][c], rN, &rremote1[rS], &ind0));
910: 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]);
911: /* Get the corresponding local point */
912: mainCone[c] = rmine1[rS + ind0];
913: }
914: 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]));
915: /* Set the desired order of p's cone points and fix orientations accordingly */
916: PetscCall(DMPolytopeGetOrientation(ct, cone, mainCone, &o));
917: PetscCall(DMPlexOrientPoint(dm, p, o));
918: } else if (debug) PetscCall(PetscSynchronizedPrintf(comm, " ==\n"));
919: }
920: if (debug) {
921: PetscCall(PetscSynchronizedFlush(comm, NULL));
922: PetscCallMPI(MPI_Barrier(comm));
923: }
924: PetscCall(DMViewFromOptions(dm, NULL, "-after_orient_interface_dm_view"));
925: PetscCall(PetscFree4(roots, leaves, rootsRanks, leavesRanks));
926: PetscCall(PetscFree2(rmine1, rremote1));
927: PetscFunctionReturn(PETSC_SUCCESS);
928: }
930: static PetscErrorCode IntArrayViewFromOptions(MPI_Comm comm, const char opt[], const char name[], const char idxname[], const char valname[], PetscInt n, const PetscInt a[])
931: {
932: PetscInt idx;
933: PetscMPIInt rank;
934: PetscBool flg;
936: PetscFunctionBegin;
937: PetscCall(PetscOptionsHasName(NULL, NULL, opt, &flg));
938: if (!flg) PetscFunctionReturn(PETSC_SUCCESS);
939: PetscCallMPI(MPI_Comm_rank(comm, &rank));
940: PetscCall(PetscSynchronizedPrintf(comm, "[%d]%s:\n", rank, name));
941: for (idx = 0; idx < n; ++idx) PetscCall(PetscSynchronizedPrintf(comm, "[%d]%s %" PetscInt_FMT " %s %" PetscInt_FMT "\n", rank, idxname, idx, valname, a[idx]));
942: PetscCall(PetscSynchronizedFlush(comm, NULL));
943: PetscFunctionReturn(PETSC_SUCCESS);
944: }
946: static PetscErrorCode SFNodeArrayViewFromOptions(MPI_Comm comm, const char opt[], const char name[], const char idxname[], PetscInt n, const PetscSFNode a[])
947: {
948: PetscInt idx;
949: PetscMPIInt rank;
950: PetscBool flg;
952: PetscFunctionBegin;
953: PetscCall(PetscOptionsHasName(NULL, NULL, opt, &flg));
954: if (!flg) PetscFunctionReturn(PETSC_SUCCESS);
955: PetscCallMPI(MPI_Comm_rank(comm, &rank));
956: PetscCall(PetscSynchronizedPrintf(comm, "[%d]%s:\n", rank, name));
957: if (idxname) {
958: for (idx = 0; idx < n; ++idx) PetscCall(PetscSynchronizedPrintf(comm, "[%d]%s %" PetscInt_FMT " rank %d index %" PetscInt_FMT "\n", rank, idxname, idx, (PetscMPIInt)a[idx].rank, a[idx].index));
959: } else {
960: for (idx = 0; idx < n; ++idx) PetscCall(PetscSynchronizedPrintf(comm, "[%d]rank %d index %" PetscInt_FMT "\n", rank, (PetscMPIInt)a[idx].rank, a[idx].index));
961: }
962: PetscCall(PetscSynchronizedFlush(comm, NULL));
963: PetscFunctionReturn(PETSC_SUCCESS);
964: }
966: static PetscErrorCode DMPlexMapToLocalPoint(DM dm, PetscHMapIJ remotehash, PetscSFNode remotePoint, PetscInt *localPoint, PetscBool *mapFailed)
967: {
968: PetscSF sf;
969: const PetscInt *locals;
970: PetscMPIInt rank;
972: PetscFunctionBegin;
973: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank));
974: PetscCall(DMGetPointSF(dm, &sf));
975: PetscCall(PetscSFGetGraph(sf, NULL, NULL, &locals, NULL));
976: if (mapFailed) *mapFailed = PETSC_FALSE;
977: if (remotePoint.rank == rank) {
978: *localPoint = remotePoint.index;
979: } else {
980: PetscHashIJKey key;
981: PetscInt l;
983: key.i = remotePoint.index;
984: key.j = remotePoint.rank;
985: PetscCall(PetscHMapIJGet(remotehash, key, &l));
986: if (l >= 0) {
987: *localPoint = locals[l];
988: } else if (mapFailed) *mapFailed = PETSC_TRUE;
989: }
990: PetscFunctionReturn(PETSC_SUCCESS);
991: }
993: static PetscErrorCode DMPlexMapToGlobalPoint(DM dm, PetscInt localPoint, PetscSFNode *remotePoint, PetscBool *mapFailed)
994: {
995: PetscSF sf;
996: const PetscInt *locals, *rootdegree;
997: const PetscSFNode *remotes;
998: PetscInt Nl, l;
999: PetscMPIInt rank;
1001: PetscFunctionBegin;
1002: if (mapFailed) *mapFailed = PETSC_FALSE;
1003: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank));
1004: PetscCall(DMGetPointSF(dm, &sf));
1005: PetscCall(PetscSFGetGraph(sf, NULL, &Nl, &locals, &remotes));
1006: if (Nl < 0) goto owned;
1007: PetscCall(PetscSFComputeDegreeBegin(sf, &rootdegree));
1008: PetscCall(PetscSFComputeDegreeEnd(sf, &rootdegree));
1009: if (rootdegree[localPoint]) goto owned;
1010: PetscCall(PetscFindInt(localPoint, Nl, locals, &l));
1011: if (l < 0) {
1012: if (mapFailed) *mapFailed = PETSC_TRUE;
1013: } else *remotePoint = remotes[l];
1014: PetscFunctionReturn(PETSC_SUCCESS);
1015: owned:
1016: remotePoint->rank = rank;
1017: remotePoint->index = localPoint;
1018: PetscFunctionReturn(PETSC_SUCCESS);
1019: }
1021: static PetscErrorCode DMPlexPointIsShared(DM dm, PetscInt p, PetscBool *isShared)
1022: {
1023: PetscSF sf;
1024: const PetscInt *locals, *rootdegree;
1025: PetscInt Nl, idx;
1027: PetscFunctionBegin;
1028: *isShared = PETSC_FALSE;
1029: PetscCall(DMGetPointSF(dm, &sf));
1030: PetscCall(PetscSFGetGraph(sf, NULL, &Nl, &locals, NULL));
1031: if (Nl < 0) PetscFunctionReturn(PETSC_SUCCESS);
1032: PetscCall(PetscFindInt(p, Nl, locals, &idx));
1033: if (idx >= 0) {
1034: *isShared = PETSC_TRUE;
1035: PetscFunctionReturn(PETSC_SUCCESS);
1036: }
1037: PetscCall(PetscSFComputeDegreeBegin(sf, &rootdegree));
1038: PetscCall(PetscSFComputeDegreeEnd(sf, &rootdegree));
1039: if (rootdegree[p] > 0) *isShared = PETSC_TRUE;
1040: PetscFunctionReturn(PETSC_SUCCESS);
1041: }
1043: static PetscErrorCode DMPlexConeIsShared(DM dm, PetscInt p, PetscBool *isShared)
1044: {
1045: const PetscInt *cone;
1046: PetscInt coneSize, c;
1047: PetscBool cShared = PETSC_TRUE;
1049: PetscFunctionBegin;
1050: PetscCall(DMPlexGetConeSize(dm, p, &coneSize));
1051: PetscCall(DMPlexGetCone(dm, p, &cone));
1052: for (c = 0; c < coneSize; ++c) {
1053: PetscBool pointShared;
1055: PetscCall(DMPlexPointIsShared(dm, cone[c], &pointShared));
1056: cShared = (PetscBool)(cShared && pointShared);
1057: }
1058: *isShared = coneSize ? cShared : PETSC_FALSE;
1059: PetscFunctionReturn(PETSC_SUCCESS);
1060: }
1062: static PetscErrorCode DMPlexGetConeMinimum(DM dm, PetscInt p, PetscSFNode *cpmin)
1063: {
1064: const PetscInt *cone;
1065: PetscInt coneSize, c;
1066: PetscSFNode cmin = {PETSC_INT_MAX, PETSC_MPI_INT_MAX}, missing = {-1, -1};
1068: PetscFunctionBegin;
1069: PetscCall(DMPlexGetConeSize(dm, p, &coneSize));
1070: PetscCall(DMPlexGetCone(dm, p, &cone));
1071: for (c = 0; c < coneSize; ++c) {
1072: PetscSFNode rcp;
1073: PetscBool mapFailed;
1075: PetscCall(DMPlexMapToGlobalPoint(dm, cone[c], &rcp, &mapFailed));
1076: if (mapFailed) {
1077: cmin = missing;
1078: } else {
1079: cmin = (rcp.rank < cmin.rank) || (rcp.rank == cmin.rank && rcp.index < cmin.index) ? rcp : cmin;
1080: }
1081: }
1082: *cpmin = coneSize ? cmin : missing;
1083: PetscFunctionReturn(PETSC_SUCCESS);
1084: }
1086: /*
1087: Each shared face has an entry in the candidates array:
1088: (-1, coneSize-1), {(global cone point)}
1089: where the set is missing the point p which we use as the key for the face
1090: */
1091: static PetscErrorCode DMPlexAddSharedFace_Private(DM dm, PetscSection candidateSection, PetscSFNode candidates[], PetscHMapIJ faceHash, PetscInt p, PetscBool debug)
1092: {
1093: MPI_Comm comm;
1094: const PetscInt *support;
1095: PetscInt supportSize, s, off = 0, idx = 0, overlap, cellHeight, height;
1096: PetscMPIInt rank;
1098: PetscFunctionBegin;
1099: PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
1100: PetscCallMPI(MPI_Comm_rank(comm, &rank));
1101: PetscCall(DMPlexGetOverlap(dm, &overlap));
1102: PetscCall(DMPlexGetVTKCellHeight(dm, &cellHeight));
1103: PetscCall(DMPlexGetPointHeight(dm, p, &height));
1104: if (!overlap && height <= cellHeight + 1) {
1105: /* cells can't be shared for non-overlapping meshes */
1106: if (debug) PetscCall(PetscSynchronizedPrintf(comm, "[%d] Skipping face %" PetscInt_FMT " to avoid adding cell to hashmap since this is nonoverlapping mesh\n", rank, p));
1107: PetscFunctionReturn(PETSC_SUCCESS);
1108: }
1109: PetscCall(DMPlexGetSupportSize(dm, p, &supportSize));
1110: PetscCall(DMPlexGetSupport(dm, p, &support));
1111: if (candidates) PetscCall(PetscSectionGetOffset(candidateSection, p, &off));
1112: for (s = 0; s < supportSize; ++s) {
1113: const PetscInt face = support[s];
1114: const PetscInt *cone;
1115: PetscSFNode cpmin = {-1, -1}, rp = {-1, -1};
1116: PetscInt coneSize, c, f;
1117: PetscBool isShared = PETSC_FALSE;
1118: PetscHashIJKey key;
1120: /* Only add point once */
1121: if (debug) PetscCall(PetscSynchronizedPrintf(comm, "[%d] Support face %" PetscInt_FMT "\n", rank, face));
1122: key.i = p;
1123: key.j = face;
1124: PetscCall(PetscHMapIJGet(faceHash, key, &f));
1125: if (f >= 0) continue;
1126: PetscCall(DMPlexConeIsShared(dm, face, &isShared));
1127: PetscCall(DMPlexGetConeMinimum(dm, face, &cpmin));
1128: PetscCall(DMPlexMapToGlobalPoint(dm, p, &rp, NULL));
1129: if (debug) {
1130: PetscCall(PetscSynchronizedPrintf(comm, "[%d] Face point %" PetscInt_FMT " is shared: %d\n", rank, face, (int)isShared));
1131: PetscCall(PetscSynchronizedPrintf(comm, "[%d] Global point (%d, %" PetscInt_FMT ") Min Cone Point (%d, %" PetscInt_FMT ")\n", rank, (PetscMPIInt)rp.rank, rp.index, (PetscMPIInt)cpmin.rank, cpmin.index));
1132: }
1133: if (isShared && (rp.rank == cpmin.rank && rp.index == cpmin.index)) {
1134: PetscCall(PetscHMapIJSet(faceHash, key, p));
1135: if (candidates) {
1136: if (debug) PetscCall(PetscSynchronizedPrintf(comm, "[%d] Adding shared face %" PetscInt_FMT " at idx %" PetscInt_FMT "\n[%d] ", rank, face, idx, rank));
1137: PetscCall(DMPlexGetConeSize(dm, face, &coneSize));
1138: PetscCall(DMPlexGetCone(dm, face, &cone));
1139: candidates[off + idx].rank = -1;
1140: candidates[off + idx++].index = coneSize - 1;
1141: candidates[off + idx].rank = rank;
1142: candidates[off + idx++].index = face;
1143: for (c = 0; c < coneSize; ++c) {
1144: const PetscInt cp = cone[c];
1146: if (cp == p) continue;
1147: PetscCall(DMPlexMapToGlobalPoint(dm, cp, &candidates[off + idx], NULL));
1148: if (debug) PetscCall(PetscSynchronizedPrintf(comm, " (%d,%" PetscInt_FMT ")", (PetscMPIInt)candidates[off + idx].rank, candidates[off + idx].index));
1149: ++idx;
1150: }
1151: if (debug) PetscCall(PetscSynchronizedPrintf(comm, "\n"));
1152: } else {
1153: /* Add cone size to section */
1154: if (debug) PetscCall(PetscSynchronizedPrintf(comm, "[%d] Scheduling shared face %" PetscInt_FMT "\n", rank, face));
1155: PetscCall(DMPlexGetConeSize(dm, face, &coneSize));
1156: PetscCall(PetscHMapIJSet(faceHash, key, p));
1157: PetscCall(PetscSectionAddDof(candidateSection, p, coneSize + 1));
1158: }
1159: }
1160: }
1161: PetscFunctionReturn(PETSC_SUCCESS);
1162: }
1164: /*@
1165: DMPlexInterpolatePointSF - Insert interpolated points in the overlap into the `PointSF` in parallel, following local interpolation
1167: Collective
1169: Input Parameters:
1170: + dm - The interpolated `DMPLEX`
1171: - pointSF - The initial `PetscSF` without interpolated points
1173: Level: developer
1175: Note:
1176: 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`
1178: .seealso: `DMPLEX`, `DMPlexInterpolate()`, `DMPlexUninterpolate()`
1179: @*/
1180: PetscErrorCode DMPlexInterpolatePointSF(DM dm, PetscSF pointSF)
1181: {
1182: MPI_Comm comm;
1183: PetscHMapIJ remoteHash;
1184: PetscHMapI claimshash;
1185: PetscSection candidateSection, candidateRemoteSection, claimSection;
1186: PetscSFNode *candidates, *candidatesRemote, *claims;
1187: const PetscInt *localPoints, *rootdegree;
1188: const PetscSFNode *remotePoints;
1189: PetscInt ov, Nr, r, Nl, l;
1190: PetscInt candidatesSize, candidatesRemoteSize, claimsSize;
1191: PetscBool flg, debug = PETSC_FALSE;
1192: PetscMPIInt rank;
1194: PetscFunctionBegin;
1197: PetscCall(DMPlexIsDistributed(dm, &flg));
1198: if (!flg) PetscFunctionReturn(PETSC_SUCCESS);
1199: /* Set initial SF so that lower level queries work */
1200: PetscCall(DMSetPointSF(dm, pointSF));
1201: PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
1202: PetscCallMPI(MPI_Comm_rank(comm, &rank));
1203: PetscCall(DMPlexGetOverlap(dm, &ov));
1204: PetscCheck(!ov, comm, PETSC_ERR_SUP, "Interpolation of overlapped DMPlex not implemented yet");
1205: PetscCall(PetscOptionsHasName(NULL, ((PetscObject)dm)->prefix, "-dmplex_interp_debug", &debug));
1206: PetscCall(PetscObjectViewFromOptions((PetscObject)dm, NULL, "-dm_interp_pre_view"));
1207: PetscCall(PetscObjectViewFromOptions((PetscObject)pointSF, NULL, "-petscsf_interp_pre_view"));
1208: PetscCall(PetscLogEventBegin(DMPLEX_InterpolateSF, dm, 0, 0, 0));
1209: /* Step 0: Precalculations */
1210: PetscCall(PetscSFGetGraph(pointSF, &Nr, &Nl, &localPoints, &remotePoints));
1211: PetscCheck(Nr >= 0, comm, PETSC_ERR_ARG_WRONGSTATE, "This DMPlex is distributed but input PointSF has no graph set");
1212: PetscCall(PetscHMapIJCreate(&remoteHash));
1213: for (l = 0; l < Nl; ++l) {
1214: PetscHashIJKey key;
1216: key.i = remotePoints[l].index;
1217: key.j = remotePoints[l].rank;
1218: PetscCall(PetscHMapIJSet(remoteHash, key, l));
1219: }
1220: /* Compute root degree to identify shared points */
1221: PetscCall(PetscSFComputeDegreeBegin(pointSF, &rootdegree));
1222: PetscCall(PetscSFComputeDegreeEnd(pointSF, &rootdegree));
1223: PetscCall(IntArrayViewFromOptions(comm, "-interp_root_degree_view", "Root degree", "point", "degree", Nr, rootdegree));
1224: /*
1225: 1) Loop over each leaf point $p$ at depth $d$ in the SF
1226: \item Get set $F(p)$ of faces $f$ in the support of $p$ for which
1227: \begin{itemize}
1228: \item all cone points of $f$ are shared
1229: \item $p$ is the cone point with smallest canonical number
1230: \end{itemize}
1231: \item Send $F(p)$ and the cone of each face to the active root point $r(p)$
1232: \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
1233: \item Send the root face from the root back to all leaf process
1234: \item Leaf processes add the shared face to the SF
1235: */
1236: /* Step 1: Construct section+SFNode array
1237: The section has entries for all shared faces for which we have a leaf point in the cone
1238: The array holds candidate shared faces, each face is referred to by the leaf point */
1239: PetscCall(PetscSectionCreate(comm, &candidateSection));
1240: PetscCall(PetscSectionSetChart(candidateSection, 0, Nr));
1241: {
1242: PetscHMapIJ faceHash;
1244: PetscCall(PetscHMapIJCreate(&faceHash));
1245: for (l = 0; l < Nl; ++l) {
1246: const PetscInt p = localPoints[l];
1248: if (debug) PetscCall(PetscSynchronizedPrintf(comm, "[%d] First pass leaf point %" PetscInt_FMT "\n", rank, p));
1249: PetscCall(DMPlexAddSharedFace_Private(dm, candidateSection, NULL, faceHash, p, debug));
1250: }
1251: PetscCall(PetscHMapIJClear(faceHash));
1252: PetscCall(PetscSectionSetUp(candidateSection));
1253: PetscCall(PetscSectionGetStorageSize(candidateSection, &candidatesSize));
1254: PetscCall(PetscMalloc1(candidatesSize, &candidates));
1255: for (l = 0; l < Nl; ++l) {
1256: const PetscInt p = localPoints[l];
1258: if (debug) PetscCall(PetscSynchronizedPrintf(comm, "[%d] Second pass leaf point %" PetscInt_FMT "\n", rank, p));
1259: PetscCall(DMPlexAddSharedFace_Private(dm, candidateSection, candidates, faceHash, p, debug));
1260: }
1261: PetscCall(PetscHMapIJDestroy(&faceHash));
1262: if (debug) PetscCall(PetscSynchronizedFlush(comm, NULL));
1263: }
1264: PetscCall(PetscObjectSetName((PetscObject)candidateSection, "Candidate Section"));
1265: PetscCall(PetscObjectViewFromOptions((PetscObject)candidateSection, NULL, "-petscsection_interp_candidate_view"));
1266: PetscCall(SFNodeArrayViewFromOptions(comm, "-petscsection_interp_candidate_view", "Candidates", NULL, candidatesSize, candidates));
1267: /* Step 2: Gather candidate section / array pair into the root partition via inverse(multi(pointSF)). */
1268: /* Note that this section is indexed by offsets into leaves, not by point number */
1269: {
1270: PetscSF sfMulti, sfInverse, sfCandidates;
1271: PetscInt *remoteOffsets;
1273: PetscCall(PetscSFGetMultiSF(pointSF, &sfMulti));
1274: PetscCall(PetscSFCreateInverseSF(sfMulti, &sfInverse));
1275: PetscCall(PetscSectionCreate(comm, &candidateRemoteSection));
1276: PetscCall(PetscSFDistributeSection(sfInverse, candidateSection, &remoteOffsets, candidateRemoteSection));
1277: PetscCall(PetscSFCreateSectionSF(sfInverse, candidateSection, remoteOffsets, candidateRemoteSection, &sfCandidates));
1278: PetscCall(PetscSectionGetStorageSize(candidateRemoteSection, &candidatesRemoteSize));
1279: PetscCall(PetscMalloc1(candidatesRemoteSize, &candidatesRemote));
1280: PetscCall(PetscSFBcastBegin(sfCandidates, MPIU_SF_NODE, candidates, candidatesRemote, MPI_REPLACE));
1281: PetscCall(PetscSFBcastEnd(sfCandidates, MPIU_SF_NODE, candidates, candidatesRemote, MPI_REPLACE));
1282: PetscCall(PetscSFDestroy(&sfInverse));
1283: PetscCall(PetscSFDestroy(&sfCandidates));
1284: PetscCall(PetscFree(remoteOffsets));
1286: PetscCall(PetscObjectSetName((PetscObject)candidateRemoteSection, "Remote Candidate Section"));
1287: PetscCall(PetscObjectViewFromOptions((PetscObject)candidateRemoteSection, NULL, "-petscsection_interp_candidate_remote_view"));
1288: PetscCall(SFNodeArrayViewFromOptions(comm, "-petscsection_interp_candidate_remote_view", "Remote Candidates", NULL, candidatesRemoteSize, candidatesRemote));
1289: }
1290: /* 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 */
1291: {
1292: PetscHMapIJKLRemote faceTable;
1293: PetscInt idx, idx2;
1295: PetscCall(PetscHMapIJKLRemoteCreate(&faceTable));
1296: /* There is a section point for every leaf attached to a given root point */
1297: for (r = 0, idx = 0, idx2 = 0; r < Nr; ++r) {
1298: PetscInt deg;
1300: for (deg = 0; deg < rootdegree[r]; ++deg, ++idx) {
1301: PetscInt offset, dof, d;
1303: PetscCall(PetscSectionGetDof(candidateRemoteSection, idx, &dof));
1304: PetscCall(PetscSectionGetOffset(candidateRemoteSection, idx, &offset));
1305: /* dof may include many faces from the remote process */
1306: for (d = 0; d < dof; ++d) {
1307: const PetscInt hidx = offset + d;
1308: const PetscInt Np = candidatesRemote[hidx].index + 1;
1309: const PetscSFNode rface = candidatesRemote[hidx + 1];
1310: const PetscSFNode *fcone = &candidatesRemote[hidx + 2];
1311: PetscSFNode fcp0;
1312: const PetscSFNode pmax = {-1, -1};
1313: const PetscInt *join = NULL;
1314: PetscHMapIJKLRemoteKey key;
1315: PetscHashIter iter;
1316: PetscBool missing, mapToLocalPointFailed = PETSC_FALSE;
1317: PetscInt points[1024], p, joinSize;
1319: if (debug)
1320: PetscCall(PetscSynchronizedPrintf(PetscObjectComm((PetscObject)dm), "[%d] Checking face (%d, %" PetscInt_FMT ") at (%" PetscInt_FMT ", %" PetscInt_FMT ", %" PetscInt_FMT ") with cone size %" PetscInt_FMT "\n", rank, (PetscMPIInt)rface.rank,
1321: rface.index, r, idx, d, Np));
1323: PetscCheck(Np <= 4, PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot handle face (%d, %" PetscInt_FMT ") at (%" PetscInt_FMT ", %" PetscInt_FMT ", %" PetscInt_FMT ") with %" PetscInt_FMT " cone points", (PetscMPIInt)rface.rank, rface.index, r, idx, d, Np);
1324: fcp0.rank = rank;
1325: fcp0.index = r;
1326: d += Np;
1327: /* Put remote face in hash table */
1328: key.i = fcp0;
1329: key.j = fcone[0];
1330: key.k = Np > 2 ? fcone[1] : pmax;
1331: key.l = Np > 3 ? fcone[2] : pmax;
1332: PetscCall(PetscSortSFNode(Np, (PetscSFNode *)&key));
1333: PetscCall(PetscHMapIJKLRemotePut(faceTable, key, &iter, &missing));
1334: if (missing) {
1335: if (debug) PetscCall(PetscSynchronizedPrintf(PetscObjectComm((PetscObject)dm), "[%d] Setting remote face (%" PetscInt_FMT ", %d)\n", rank, rface.index, (PetscMPIInt)rface.rank));
1336: PetscCall(PetscHMapIJKLRemoteIterSet(faceTable, iter, rface));
1337: } else {
1338: PetscSFNode oface;
1340: PetscCall(PetscHMapIJKLRemoteIterGet(faceTable, iter, &oface));
1341: if ((rface.rank < oface.rank) || (rface.rank == oface.rank && rface.index < oface.index)) {
1342: if (debug) PetscCall(PetscSynchronizedPrintf(PetscObjectComm((PetscObject)dm), "[%d] Replacing with remote face (%" PetscInt_FMT ", %d)\n", rank, rface.index, (PetscMPIInt)rface.rank));
1343: PetscCall(PetscHMapIJKLRemoteIterSet(faceTable, iter, rface));
1344: }
1345: }
1346: /* Check for local face */
1347: points[0] = r;
1348: for (p = 1; p < Np; ++p) {
1349: PetscCall(DMPlexMapToLocalPoint(dm, remoteHash, fcone[p - 1], &points[p], &mapToLocalPointFailed));
1350: if (mapToLocalPointFailed) break; /* We got a point not in our overlap */
1351: if (debug) PetscCall(PetscSynchronizedPrintf(PetscObjectComm((PetscObject)dm), "[%d] Checking local candidate %" PetscInt_FMT "\n", rank, points[p]));
1352: }
1353: if (mapToLocalPointFailed) continue;
1354: PetscCall(DMPlexGetJoin(dm, Np, points, &joinSize, &join));
1355: if (joinSize == 1) {
1356: PetscSFNode lface;
1357: PetscSFNode oface;
1359: /* Always replace with local face */
1360: lface.rank = rank;
1361: lface.index = join[0];
1362: PetscCall(PetscHMapIJKLRemoteIterGet(faceTable, iter, &oface));
1363: if (debug)
1364: PetscCall(PetscSynchronizedPrintf(PetscObjectComm((PetscObject)dm), "[%d] Replacing (%" PetscInt_FMT ", %d) with local face (%" PetscInt_FMT ", %d)\n", rank, oface.index, (PetscMPIInt)oface.rank, lface.index, (PetscMPIInt)lface.rank));
1365: PetscCall(PetscHMapIJKLRemoteIterSet(faceTable, iter, lface));
1366: }
1367: PetscCall(DMPlexRestoreJoin(dm, Np, points, &joinSize, &join));
1368: }
1369: }
1370: /* Put back faces for this root */
1371: for (deg = 0; deg < rootdegree[r]; ++deg, ++idx2) {
1372: PetscInt offset, dof, d;
1374: PetscCall(PetscSectionGetDof(candidateRemoteSection, idx2, &dof));
1375: PetscCall(PetscSectionGetOffset(candidateRemoteSection, idx2, &offset));
1376: /* dof may include many faces from the remote process */
1377: for (d = 0; d < dof; ++d) {
1378: const PetscInt hidx = offset + d;
1379: const PetscInt Np = candidatesRemote[hidx].index + 1;
1380: const PetscSFNode *fcone = &candidatesRemote[hidx + 2];
1381: PetscSFNode fcp0;
1382: const PetscSFNode pmax = {-1, -1};
1383: PetscHMapIJKLRemoteKey key;
1384: PetscHashIter iter;
1385: PetscBool missing;
1387: if (debug) PetscCall(PetscSynchronizedPrintf(PetscObjectComm((PetscObject)dm), "[%d] Entering face at (%" PetscInt_FMT ", %" PetscInt_FMT ")\n", rank, r, idx));
1388: PetscCheck(Np <= 4, PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot handle faces with %" PetscInt_FMT " cone points", Np);
1389: fcp0.rank = rank;
1390: fcp0.index = r;
1391: d += Np;
1392: /* Find remote face in hash table */
1393: key.i = fcp0;
1394: key.j = fcone[0];
1395: key.k = Np > 2 ? fcone[1] : pmax;
1396: key.l = Np > 3 ? fcone[2] : pmax;
1397: PetscCall(PetscSortSFNode(Np, (PetscSFNode *)&key));
1398: if (debug)
1399: PetscCall(PetscSynchronizedPrintf(PetscObjectComm((PetscObject)dm), "[%d] key (%d, %" PetscInt_FMT ") (%d, %" PetscInt_FMT ") (%d, %" PetscInt_FMT ") (%d, %" PetscInt_FMT ")\n", rank, (PetscMPIInt)key.i.rank, key.i.index,
1400: (PetscMPIInt)key.j.rank, key.j.index, (PetscMPIInt)key.k.rank, key.k.index, (PetscMPIInt)key.l.rank, key.l.index));
1401: PetscCall(PetscHMapIJKLRemotePut(faceTable, key, &iter, &missing));
1402: PetscCheck(!missing, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Root %" PetscInt_FMT " Idx %" PetscInt_FMT " ought to have an associated face", r, idx2);
1403: PetscCall(PetscHMapIJKLRemoteIterGet(faceTable, iter, &candidatesRemote[hidx]));
1404: }
1405: }
1406: }
1407: if (debug) PetscCall(PetscSynchronizedFlush(PetscObjectComm((PetscObject)dm), NULL));
1408: PetscCall(PetscHMapIJKLRemoteDestroy(&faceTable));
1409: }
1410: /* Step 4: Push back owned faces */
1411: {
1412: PetscSF sfMulti, sfClaims, sfPointNew;
1413: PetscSFNode *remotePointsNew;
1414: PetscInt *remoteOffsets, *localPointsNew;
1415: PetscInt pStart, pEnd, r, NlNew, p;
1417: /* 4) Push claims back to receiver via the MultiSF and derive new pointSF mapping on receiver */
1418: PetscCall(PetscSFGetMultiSF(pointSF, &sfMulti));
1419: PetscCall(PetscSectionCreate(comm, &claimSection));
1420: PetscCall(PetscSFDistributeSection(sfMulti, candidateRemoteSection, &remoteOffsets, claimSection));
1421: PetscCall(PetscSFCreateSectionSF(sfMulti, candidateRemoteSection, remoteOffsets, claimSection, &sfClaims));
1422: PetscCall(PetscSectionGetStorageSize(claimSection, &claimsSize));
1423: PetscCall(PetscMalloc1(claimsSize, &claims));
1424: for (p = 0; p < claimsSize; ++p) claims[p].rank = -1;
1425: PetscCall(PetscSFBcastBegin(sfClaims, MPIU_SF_NODE, candidatesRemote, claims, MPI_REPLACE));
1426: PetscCall(PetscSFBcastEnd(sfClaims, MPIU_SF_NODE, candidatesRemote, claims, MPI_REPLACE));
1427: PetscCall(PetscSFDestroy(&sfClaims));
1428: PetscCall(PetscFree(remoteOffsets));
1429: PetscCall(PetscObjectSetName((PetscObject)claimSection, "Claim Section"));
1430: PetscCall(PetscObjectViewFromOptions((PetscObject)claimSection, NULL, "-petscsection_interp_claim_view"));
1431: PetscCall(SFNodeArrayViewFromOptions(comm, "-petscsection_interp_claim_view", "Claims", NULL, claimsSize, claims));
1432: /* Step 5) Walk the original section of local supports and add an SF entry for each updated item */
1433: /* TODO I should not have to do a join here since I already put the face and its cone in the candidate section */
1434: PetscCall(PetscHMapICreate(&claimshash));
1435: for (r = 0; r < Nr; ++r) {
1436: PetscInt dof, off, d;
1438: if (debug) PetscCall(PetscSynchronizedPrintf(comm, "[%d] Checking root for claims %" PetscInt_FMT "\n", rank, r));
1439: PetscCall(PetscSectionGetDof(candidateSection, r, &dof));
1440: PetscCall(PetscSectionGetOffset(candidateSection, r, &off));
1441: for (d = 0; d < dof;) {
1442: if (claims[off + d].rank >= 0) {
1443: const PetscInt faceInd = off + d;
1444: const PetscInt Np = candidates[off + d].index;
1445: const PetscInt *join = NULL;
1446: PetscInt joinSize, points[1024], c;
1448: if (debug) PetscCall(PetscSynchronizedPrintf(comm, "[%d] Found claim for remote point (%d, %" PetscInt_FMT ")\n", rank, (PetscMPIInt)claims[faceInd].rank, claims[faceInd].index));
1449: points[0] = r;
1450: if (debug) PetscCall(PetscSynchronizedPrintf(comm, "[%d] point %" PetscInt_FMT "\n", rank, points[0]));
1451: for (c = 0, d += 2; c < Np; ++c, ++d) {
1452: PetscCall(DMPlexMapToLocalPoint(dm, remoteHash, candidates[off + d], &points[c + 1], NULL));
1453: if (debug) PetscCall(PetscSynchronizedPrintf(comm, "[%d] point %" PetscInt_FMT "\n", rank, points[c + 1]));
1454: }
1455: PetscCall(DMPlexGetJoin(dm, Np + 1, points, &joinSize, &join));
1456: if (joinSize == 1) {
1457: if (claims[faceInd].rank == rank) {
1458: if (debug) PetscCall(PetscSynchronizedPrintf(comm, "[%d] Ignoring local face %" PetscInt_FMT " for non-remote partner\n", rank, join[0]));
1459: } else {
1460: if (debug) PetscCall(PetscSynchronizedPrintf(comm, "[%d] Found local face %" PetscInt_FMT "\n", rank, join[0]));
1461: PetscCall(PetscHMapISet(claimshash, join[0], faceInd));
1462: }
1463: } else {
1464: if (debug) PetscCall(PetscSynchronizedPrintf(comm, "[%d] Failed to find face\n", rank));
1465: }
1466: PetscCall(DMPlexRestoreJoin(dm, Np + 1, points, &joinSize, &join));
1467: } else {
1468: if (debug) PetscCall(PetscSynchronizedPrintf(comm, "[%d] No claim for point %" PetscInt_FMT "\n", rank, r));
1469: d += claims[off + d].index + 1;
1470: }
1471: }
1472: }
1473: if (debug) PetscCall(PetscSynchronizedFlush(comm, NULL));
1474: /* Step 6) Create new pointSF from hashed claims */
1475: PetscCall(PetscHMapIGetSize(claimshash, &NlNew));
1476: PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
1477: PetscCall(PetscMalloc1(Nl + NlNew, &localPointsNew));
1478: PetscCall(PetscMalloc1(Nl + NlNew, &remotePointsNew));
1479: for (l = 0; l < Nl; ++l) {
1480: localPointsNew[l] = localPoints[l];
1481: remotePointsNew[l].index = remotePoints[l].index;
1482: remotePointsNew[l].rank = remotePoints[l].rank;
1483: }
1484: p = Nl;
1485: PetscCall(PetscHMapIGetKeys(claimshash, &p, localPointsNew));
1486: /* We sort new points, and assume they are numbered after all existing points */
1487: PetscCall(PetscSortInt(NlNew, PetscSafePointerPlusOffset(localPointsNew, Nl)));
1488: for (p = Nl; p < Nl + NlNew; ++p) {
1489: PetscInt off;
1490: PetscCall(PetscHMapIGet(claimshash, localPointsNew[p], &off));
1491: PetscCheck(claims[off].rank >= 0 && claims[off].index >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid claim for local point %" PetscInt_FMT ", (%d, %" PetscInt_FMT ")", localPointsNew[p], (PetscMPIInt)claims[off].rank, claims[off].index);
1492: remotePointsNew[p] = claims[off];
1493: }
1494: PetscCall(PetscSFCreate(comm, &sfPointNew));
1495: PetscCall(PetscSFSetGraph(sfPointNew, pEnd - pStart, Nl + NlNew, localPointsNew, PETSC_OWN_POINTER, remotePointsNew, PETSC_OWN_POINTER));
1496: PetscCall(PetscSFSetUp(sfPointNew));
1497: PetscCall(DMSetPointSF(dm, sfPointNew));
1498: PetscCall(PetscObjectViewFromOptions((PetscObject)sfPointNew, NULL, "-petscsf_interp_view"));
1499: if (PetscDefined(USE_DEBUG)) PetscCall(DMPlexCheckPointSF(dm, sfPointNew, PETSC_FALSE));
1500: PetscCall(PetscSFDestroy(&sfPointNew));
1501: PetscCall(PetscHMapIDestroy(&claimshash));
1502: }
1503: PetscCall(PetscHMapIJDestroy(&remoteHash));
1504: PetscCall(PetscSectionDestroy(&candidateSection));
1505: PetscCall(PetscSectionDestroy(&candidateRemoteSection));
1506: PetscCall(PetscSectionDestroy(&claimSection));
1507: PetscCall(PetscFree(candidates));
1508: PetscCall(PetscFree(candidatesRemote));
1509: PetscCall(PetscFree(claims));
1510: PetscCall(PetscLogEventEnd(DMPLEX_InterpolateSF, dm, 0, 0, 0));
1511: PetscFunctionReturn(PETSC_SUCCESS);
1512: }
1514: /*@
1515: DMPlexInterpolate - Take in a cell-vertex mesh and return one with all intermediate faces, edges, etc.
1517: Collective
1519: Input Parameter:
1520: . dm - The `DMPLEX` object with only cells and vertices
1522: Output Parameter:
1523: . dmInt - The complete `DMPLEX` object
1525: Level: intermediate
1527: Note:
1528: Labels and coordinates are copied.
1530: Developer Notes:
1531: It sets plex->interpolated = `DMPLEX_INTERPOLATED_FULL`.
1533: .seealso: `DMPLEX`, `DMPlexUninterpolate()`, `DMPlexCreateFromCellListPetsc()`, `DMPlexCopyCoordinates()`
1534: @*/
1535: PetscErrorCode DMPlexInterpolate(DM dm, DM *dmInt)
1536: {
1537: DMPlexInterpolatedFlag interpolated;
1538: DM idm, odm = dm;
1539: PetscSF sfPoint;
1540: PetscInt depth, dim, d;
1541: const char *name;
1542: PetscBool flg = PETSC_TRUE;
1544: PetscFunctionBegin;
1546: PetscAssertPointer(dmInt, 2);
1547: PetscCall(PetscLogEventBegin(DMPLEX_Interpolate, dm, 0, 0, 0));
1548: PetscCall(DMPlexGetDepth(dm, &depth));
1549: PetscCall(DMGetDimension(dm, &dim));
1550: PetscCall(DMPlexIsInterpolated(dm, &interpolated));
1551: PetscCheck(interpolated != DMPLEX_INTERPOLATED_PARTIAL, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Not for partially interpolated meshes");
1552: if (interpolated == DMPLEX_INTERPOLATED_FULL) {
1553: PetscCall(PetscObjectReference((PetscObject)dm));
1554: idm = dm;
1555: } else {
1556: PetscBool nonmanifold = PETSC_FALSE;
1558: PetscCall(PetscOptionsGetBool(NULL, dm->hdr.prefix, "-dm_plex_stratify_celltype", &nonmanifold, NULL));
1559: if (nonmanifold) {
1560: do {
1561: const char *prefix;
1562: PetscInt pStart, pEnd, pdepth;
1563: PetscBool done = PETSC_TRUE;
1565: // Find a point which is not correctly interpolated
1566: PetscCall(DMPlexGetChart(odm, &pStart, &pEnd));
1567: for (PetscInt p = pStart; p < pEnd; ++p) {
1568: DMPolytopeType ct;
1569: const PetscInt *cone;
1570: PetscInt coneSize, cdepth;
1572: PetscCall(DMPlexGetPointDepth(odm, p, &pdepth));
1573: PetscCall(DMPlexGetCellType(odm, p, &ct));
1574: // Check against celltype
1575: if (pdepth != DMPolytopeTypeGetDim(ct)) {
1576: done = PETSC_FALSE;
1577: break;
1578: }
1579: // Check against boundary
1580: PetscCall(DMPlexGetCone(odm, p, &cone));
1581: PetscCall(DMPlexGetConeSize(odm, p, &coneSize));
1582: for (PetscInt c = 0; c < coneSize; ++c) {
1583: PetscCall(DMPlexGetPointDepth(odm, cone[c], &cdepth));
1584: if (cdepth != pdepth - 1) {
1585: done = PETSC_FALSE;
1586: p = pEnd;
1587: break;
1588: }
1589: }
1590: }
1591: if (done) break;
1592: /* Create interpolated mesh */
1593: PetscCall(DMCreate(PetscObjectComm((PetscObject)dm), &idm));
1594: PetscCall(DMSetType(idm, DMPLEX));
1595: PetscCall(DMSetDimension(idm, dim));
1596: PetscCall(PetscObjectGetOptionsPrefix((PetscObject)dm, &prefix));
1597: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)idm, prefix));
1598: if (depth > 0) {
1599: PetscCall(DMPlexInterpolateFaces_Internal(odm, pdepth, idm));
1600: PetscCall(DMGetPointSF(odm, &sfPoint));
1601: if (PetscDefined(USE_DEBUG)) PetscCall(DMPlexCheckPointSF(odm, sfPoint, PETSC_FALSE));
1602: {
1603: /* TODO: We need to systematically fix cases of distributed Plexes with no graph set */
1604: PetscInt nroots;
1605: PetscCall(PetscSFGetGraph(sfPoint, &nroots, NULL, NULL, NULL));
1606: if (nroots >= 0) PetscCall(DMPlexInterpolatePointSF(idm, sfPoint));
1607: }
1608: }
1609: if (odm != dm) PetscCall(DMDestroy(&odm));
1610: odm = idm;
1611: } while (1);
1612: } else {
1613: for (d = 1; d < dim; ++d) {
1614: const char *prefix;
1616: /* Create interpolated mesh */
1617: PetscCall(DMCreate(PetscObjectComm((PetscObject)dm), &idm));
1618: PetscCall(DMSetType(idm, DMPLEX));
1619: PetscCall(DMSetDimension(idm, dim));
1620: PetscCall(PetscObjectGetOptionsPrefix((PetscObject)dm, &prefix));
1621: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)idm, prefix));
1622: if (depth > 0) {
1623: PetscCall(DMPlexInterpolateFaces_Internal(odm, 1, idm));
1624: PetscCall(DMGetPointSF(odm, &sfPoint));
1625: if (PetscDefined(USE_DEBUG)) PetscCall(DMPlexCheckPointSF(odm, sfPoint, PETSC_FALSE));
1626: {
1627: /* TODO: We need to systematically fix cases of distributed Plexes with no graph set */
1628: PetscInt nroots;
1629: PetscCall(PetscSFGetGraph(sfPoint, &nroots, NULL, NULL, NULL));
1630: if (nroots >= 0) PetscCall(DMPlexInterpolatePointSF(idm, sfPoint));
1631: }
1632: }
1633: if (odm != dm) PetscCall(DMDestroy(&odm));
1634: odm = idm;
1635: }
1636: }
1637: PetscCall(PetscObjectGetName((PetscObject)dm, &name));
1638: PetscCall(PetscObjectSetName((PetscObject)idm, name));
1639: PetscCall(DMPlexCopyCoordinates(dm, idm));
1640: PetscCall(DMCopyLabels(dm, idm, PETSC_COPY_VALUES, PETSC_FALSE, DM_COPY_LABELS_FAIL));
1641: PetscCall(PetscOptionsGetBool(((PetscObject)dm)->options, ((PetscObject)dm)->prefix, "-dm_plex_interpolate_orient_interfaces", &flg, NULL));
1642: if (flg) PetscCall(DMPlexOrientInterface_Internal(idm));
1643: }
1644: /* This function makes the mesh fully interpolated on all ranks */
1645: {
1646: DM_Plex *plex = (DM_Plex *)idm->data;
1647: plex->interpolated = plex->interpolatedCollective = DMPLEX_INTERPOLATED_FULL;
1648: }
1649: PetscCall(DMPlexCopy_Internal(dm, PETSC_TRUE, PETSC_TRUE, idm));
1650: *dmInt = idm;
1651: PetscCall(PetscLogEventEnd(DMPLEX_Interpolate, dm, 0, 0, 0));
1652: PetscFunctionReturn(PETSC_SUCCESS);
1653: }
1655: /*@
1656: DMPlexCopyCoordinates - Copy coordinates from one mesh to another with the same vertices
1658: Collective
1660: Input Parameter:
1661: . dmA - The `DMPLEX` object with initial coordinates
1663: Output Parameter:
1664: . dmB - The `DMPLEX` object with copied coordinates
1666: Level: intermediate
1668: Notes:
1669: This is typically used when adding pieces other than vertices to a mesh
1671: This function does not copy localized coordinates.
1673: .seealso: `DMPLEX`, `DMCopyLabels()`, `DMGetCoordinates()`, `DMGetCoordinatesLocal()`, `DMGetCoordinateDM()`, `DMGetCoordinateSection()`
1674: @*/
1675: PetscErrorCode DMPlexCopyCoordinates(DM dmA, DM dmB)
1676: {
1677: Vec coordinatesA, coordinatesB;
1678: VecType vtype;
1679: PetscSection coordSectionA, coordSectionB;
1680: PetscScalar *coordsA, *coordsB;
1681: PetscInt spaceDim, Nf, vStartA, vStartB, vEndA, vEndB, coordSizeB, v, d;
1682: PetscInt cStartA, cEndA, cStartB, cEndB, cS, cE, cdim;
1684: PetscFunctionBegin;
1687: if (dmA == dmB) PetscFunctionReturn(PETSC_SUCCESS);
1688: PetscCall(DMGetCoordinateDim(dmA, &cdim));
1689: PetscCall(DMSetCoordinateDim(dmB, cdim));
1690: PetscCall(DMPlexGetDepthStratum(dmA, 0, &vStartA, &vEndA));
1691: PetscCall(DMPlexGetDepthStratum(dmB, 0, &vStartB, &vEndB));
1692: 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);
1693: /* Copy over discretization if it exists */
1694: {
1695: DM cdmA, cdmB;
1696: PetscDS dsA, dsB;
1697: PetscObject objA, objB;
1698: PetscClassId idA, idB;
1699: const PetscScalar *constants;
1700: PetscInt cdim, Nc;
1702: PetscCall(DMGetCoordinateDM(dmA, &cdmA));
1703: PetscCall(DMGetCoordinateDM(dmB, &cdmB));
1704: PetscCall(DMGetField(cdmA, 0, NULL, &objA));
1705: PetscCall(DMGetField(cdmB, 0, NULL, &objB));
1706: PetscCall(PetscObjectGetClassId(objA, &idA));
1707: PetscCall(PetscObjectGetClassId(objB, &idB));
1708: if ((idA == PETSCFE_CLASSID) && (idA != idB)) {
1709: PetscCall(DMSetField(cdmB, 0, NULL, objA));
1710: PetscCall(DMCreateDS(cdmB));
1711: PetscCall(DMGetDS(cdmA, &dsA));
1712: PetscCall(DMGetDS(cdmB, &dsB));
1713: PetscCall(PetscDSGetCoordinateDimension(dsA, &cdim));
1714: PetscCall(PetscDSSetCoordinateDimension(dsB, cdim));
1715: PetscCall(PetscDSGetConstants(dsA, &Nc, &constants));
1716: PetscCall(PetscDSSetConstants(dsB, Nc, (PetscScalar *)constants));
1717: }
1718: }
1719: PetscCall(DMPlexGetHeightStratum(dmA, 0, &cStartA, &cEndA));
1720: PetscCall(DMPlexGetHeightStratum(dmB, 0, &cStartB, &cEndB));
1721: PetscCall(DMGetCoordinateSection(dmA, &coordSectionA));
1722: PetscCall(DMGetCoordinateSection(dmB, &coordSectionB));
1723: if (coordSectionA == coordSectionB) PetscFunctionReturn(PETSC_SUCCESS);
1724: PetscCall(PetscSectionGetNumFields(coordSectionA, &Nf));
1725: if (!Nf) PetscFunctionReturn(PETSC_SUCCESS);
1726: PetscCheck(Nf <= 1, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "The number of coordinate fields must be 1, not %" PetscInt_FMT, Nf);
1727: if (!coordSectionB) {
1728: PetscInt dim;
1730: PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)coordSectionA), &coordSectionB));
1731: PetscCall(DMGetCoordinateDim(dmA, &dim));
1732: PetscCall(DMSetCoordinateSection(dmB, dim, coordSectionB));
1733: PetscCall(PetscObjectDereference((PetscObject)coordSectionB));
1734: }
1735: PetscCall(PetscSectionSetNumFields(coordSectionB, 1));
1736: PetscCall(PetscSectionGetFieldComponents(coordSectionA, 0, &spaceDim));
1737: PetscCall(PetscSectionSetFieldComponents(coordSectionB, 0, spaceDim));
1738: PetscCall(PetscSectionGetChart(coordSectionA, &cS, &cE));
1739: cS = vStartB;
1740: cE = vEndB;
1741: PetscCall(PetscSectionSetChart(coordSectionB, cS, cE));
1742: for (v = vStartB; v < vEndB; ++v) {
1743: PetscCall(PetscSectionSetDof(coordSectionB, v, spaceDim));
1744: PetscCall(PetscSectionSetFieldDof(coordSectionB, v, 0, spaceDim));
1745: }
1746: PetscCall(PetscSectionSetUp(coordSectionB));
1747: PetscCall(PetscSectionGetStorageSize(coordSectionB, &coordSizeB));
1748: PetscCall(DMGetCoordinatesLocal(dmA, &coordinatesA));
1749: PetscCall(VecCreate(PETSC_COMM_SELF, &coordinatesB));
1750: PetscCall(PetscObjectSetName((PetscObject)coordinatesB, "coordinates"));
1751: PetscCall(VecSetSizes(coordinatesB, coordSizeB, PETSC_DETERMINE));
1752: PetscCall(VecGetBlockSize(coordinatesA, &d));
1753: PetscCall(VecSetBlockSize(coordinatesB, d));
1754: PetscCall(VecGetType(coordinatesA, &vtype));
1755: PetscCall(VecSetType(coordinatesB, vtype));
1756: PetscCall(VecGetArray(coordinatesA, &coordsA));
1757: PetscCall(VecGetArray(coordinatesB, &coordsB));
1758: for (v = 0; v < vEndB - vStartB; ++v) {
1759: PetscInt offA, offB;
1761: PetscCall(PetscSectionGetOffset(coordSectionA, v + vStartA, &offA));
1762: PetscCall(PetscSectionGetOffset(coordSectionB, v + vStartB, &offB));
1763: for (d = 0; d < spaceDim; ++d) coordsB[offB + d] = coordsA[offA + d];
1764: }
1765: PetscCall(VecRestoreArray(coordinatesA, &coordsA));
1766: PetscCall(VecRestoreArray(coordinatesB, &coordsB));
1767: PetscCall(DMSetCoordinatesLocal(dmB, coordinatesB));
1768: PetscCall(VecDestroy(&coordinatesB));
1769: PetscFunctionReturn(PETSC_SUCCESS);
1770: }
1772: /*@
1773: DMPlexUninterpolate - Take in a mesh with all intermediate faces, edges, etc. and return a cell-vertex mesh
1775: Collective
1777: Input Parameter:
1778: . dm - The complete `DMPLEX` object
1780: Output Parameter:
1781: . dmUnint - The `DMPLEX` object with only cells and vertices
1783: Level: intermediate
1785: Note:
1786: It does not copy over the coordinates.
1788: Developer Notes:
1789: Sets plex->interpolated = `DMPLEX_INTERPOLATED_NONE`.
1791: .seealso: `DMPLEX`, `DMPlexInterpolate()`, `DMPlexCreateFromCellListPetsc()`, `DMPlexCopyCoordinates()`
1792: @*/
1793: PetscErrorCode DMPlexUninterpolate(DM dm, DM *dmUnint)
1794: {
1795: DMPlexInterpolatedFlag interpolated;
1796: DM udm;
1797: PetscInt dim, vStart, vEnd, cStart, cEnd, c, maxConeSize = 0, *cone;
1799: PetscFunctionBegin;
1801: PetscAssertPointer(dmUnint, 2);
1802: PetscCall(PetscLogEventBegin(DMPLEX_Uninterpolate, dm, 0, 0, 0));
1803: PetscCall(DMGetDimension(dm, &dim));
1804: PetscCall(DMPlexIsInterpolated(dm, &interpolated));
1805: PetscCheck(interpolated != DMPLEX_INTERPOLATED_PARTIAL, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Not for partially interpolated meshes");
1806: if (interpolated == DMPLEX_INTERPOLATED_NONE || dim <= 1) {
1807: /* in case dim <= 1 just keep the DMPLEX_INTERPOLATED_FULL flag */
1808: PetscCall(PetscObjectReference((PetscObject)dm));
1809: *dmUnint = dm;
1810: PetscFunctionReturn(PETSC_SUCCESS);
1811: }
1812: PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd));
1813: PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
1814: PetscCall(DMCreate(PetscObjectComm((PetscObject)dm), &udm));
1815: PetscCall(DMSetType(udm, DMPLEX));
1816: PetscCall(DMSetDimension(udm, dim));
1817: PetscCall(DMPlexSetChart(udm, cStart, vEnd));
1818: for (c = cStart; c < cEnd; ++c) {
1819: PetscInt *closure = NULL, closureSize, cl, coneSize = 0;
1821: PetscCall(DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure));
1822: for (cl = 0; cl < closureSize * 2; cl += 2) {
1823: const PetscInt p = closure[cl];
1825: if ((p >= vStart) && (p < vEnd)) ++coneSize;
1826: }
1827: PetscCall(DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure));
1828: PetscCall(DMPlexSetConeSize(udm, c, coneSize));
1829: maxConeSize = PetscMax(maxConeSize, coneSize);
1830: }
1831: PetscCall(DMSetUp(udm));
1832: PetscCall(PetscMalloc1(maxConeSize, &cone));
1833: for (c = cStart; c < cEnd; ++c) {
1834: PetscInt *closure = NULL, closureSize, cl, coneSize = 0;
1836: PetscCall(DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure));
1837: for (cl = 0; cl < closureSize * 2; cl += 2) {
1838: const PetscInt p = closure[cl];
1840: if ((p >= vStart) && (p < vEnd)) cone[coneSize++] = p;
1841: }
1842: PetscCall(DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure));
1843: PetscCall(DMPlexSetCone(udm, c, cone));
1844: }
1845: PetscCall(PetscFree(cone));
1846: PetscCall(DMPlexSymmetrize(udm));
1847: PetscCall(DMPlexStratify(udm));
1848: /* Reduce SF */
1849: {
1850: PetscSF sfPoint, sfPointUn;
1851: const PetscSFNode *remotePoints;
1852: const PetscInt *localPoints;
1853: PetscSFNode *remotePointsUn;
1854: PetscInt *localPointsUn;
1855: PetscInt numRoots, numLeaves, l;
1856: PetscInt numLeavesUn = 0, n = 0;
1858: /* Get original SF information */
1859: PetscCall(DMGetPointSF(dm, &sfPoint));
1860: if (PetscDefined(USE_DEBUG)) PetscCall(DMPlexCheckPointSF(dm, sfPoint, PETSC_FALSE));
1861: PetscCall(DMGetPointSF(udm, &sfPointUn));
1862: PetscCall(PetscSFGetGraph(sfPoint, &numRoots, &numLeaves, &localPoints, &remotePoints));
1863: if (numRoots >= 0) {
1864: /* Allocate space for cells and vertices */
1865: for (l = 0; l < numLeaves; ++l) {
1866: const PetscInt p = localPoints[l];
1868: if ((vStart <= p && p < vEnd) || (cStart <= p && p < cEnd)) numLeavesUn++;
1869: }
1870: /* Fill in leaves */
1871: PetscCall(PetscMalloc1(numLeavesUn, &remotePointsUn));
1872: PetscCall(PetscMalloc1(numLeavesUn, &localPointsUn));
1873: for (l = 0; l < numLeaves; l++) {
1874: const PetscInt p = localPoints[l];
1876: if ((vStart <= p && p < vEnd) || (cStart <= p && p < cEnd)) {
1877: localPointsUn[n] = p;
1878: remotePointsUn[n].rank = remotePoints[l].rank;
1879: remotePointsUn[n].index = remotePoints[l].index;
1880: ++n;
1881: }
1882: }
1883: PetscCheck(n == numLeavesUn, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Inconsistent number of leaves %" PetscInt_FMT " != %" PetscInt_FMT, n, numLeavesUn);
1884: PetscCall(PetscSFSetGraph(sfPointUn, cEnd - cStart + vEnd - vStart, numLeavesUn, localPointsUn, PETSC_OWN_POINTER, remotePointsUn, PETSC_OWN_POINTER));
1885: }
1886: }
1887: /* This function makes the mesh fully uninterpolated on all ranks */
1888: {
1889: DM_Plex *plex = (DM_Plex *)udm->data;
1890: plex->interpolated = plex->interpolatedCollective = DMPLEX_INTERPOLATED_NONE;
1891: }
1892: PetscCall(DMPlexCopy_Internal(dm, PETSC_TRUE, PETSC_TRUE, udm));
1893: if (PetscDefined(USE_DEBUG)) PetscCall(DMPlexCheckPointSF(udm, NULL, PETSC_FALSE));
1894: *dmUnint = udm;
1895: PetscCall(PetscLogEventEnd(DMPLEX_Uninterpolate, dm, 0, 0, 0));
1896: PetscFunctionReturn(PETSC_SUCCESS);
1897: }
1899: static PetscErrorCode DMPlexIsInterpolated_Internal(DM dm, DMPlexInterpolatedFlag *interpolated)
1900: {
1901: PetscInt coneSize, depth, dim, h, p, pStart, pEnd;
1902: MPI_Comm comm;
1904: PetscFunctionBegin;
1905: PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
1906: PetscCall(DMPlexGetDepth(dm, &depth));
1907: PetscCall(DMGetDimension(dm, &dim));
1909: if (depth == dim) {
1910: *interpolated = DMPLEX_INTERPOLATED_FULL;
1911: if (!dim) goto finish;
1913: /* Check points at height = dim are vertices (have no cones) */
1914: PetscCall(DMPlexGetHeightStratum(dm, dim, &pStart, &pEnd));
1915: for (p = pStart; p < pEnd; p++) {
1916: PetscCall(DMPlexGetConeSize(dm, p, &coneSize));
1917: if (coneSize) {
1918: *interpolated = DMPLEX_INTERPOLATED_PARTIAL;
1919: goto finish;
1920: }
1921: }
1923: /* Check points at height < dim have cones */
1924: for (h = 0; h < dim; h++) {
1925: PetscCall(DMPlexGetHeightStratum(dm, h, &pStart, &pEnd));
1926: for (p = pStart; p < pEnd; p++) {
1927: PetscCall(DMPlexGetConeSize(dm, p, &coneSize));
1928: if (!coneSize) {
1929: *interpolated = DMPLEX_INTERPOLATED_PARTIAL;
1930: goto finish;
1931: }
1932: }
1933: }
1934: } else if (depth == 1) {
1935: *interpolated = DMPLEX_INTERPOLATED_NONE;
1936: } else {
1937: *interpolated = DMPLEX_INTERPOLATED_PARTIAL;
1938: }
1939: finish:
1940: PetscFunctionReturn(PETSC_SUCCESS);
1941: }
1943: /*@
1944: DMPlexIsInterpolated - Find out to what extent the `DMPLEX` is topologically interpolated.
1946: Not Collective
1948: Input Parameter:
1949: . dm - The `DMPLEX` object
1951: Output Parameter:
1952: . interpolated - Flag whether the `DM` is interpolated
1954: Level: intermediate
1956: Notes:
1957: Unlike `DMPlexIsInterpolatedCollective()`, this is NOT collective
1958: so the results can be different on different ranks in special cases.
1959: However, `DMPlexInterpolate()` guarantees the result is the same on all.
1961: Unlike `DMPlexIsInterpolatedCollective()`, this cannot return `DMPLEX_INTERPOLATED_MIXED`.
1963: Developer Notes:
1964: Initially, plex->interpolated = `DMPLEX_INTERPOLATED_INVALID`.
1966: If plex->interpolated == `DMPLEX_INTERPOLATED_INVALID`, `DMPlexIsInterpolated_Internal()` is called.
1967: It checks the actual topology and sets plex->interpolated on each rank separately to one of
1968: `DMPLEX_INTERPOLATED_NONE`, `DMPLEX_INTERPOLATED_PARTIAL` or `DMPLEX_INTERPOLATED_FULL`.
1970: If plex->interpolated != `DMPLEX_INTERPOLATED_INVALID`, this function just returns plex->interpolated.
1972: `DMPlexInterpolate()` sets plex->interpolated = `DMPLEX_INTERPOLATED_FULL`,
1973: and DMPlexUninterpolate() sets plex->interpolated = `DMPLEX_INTERPOLATED_NONE`.
1975: .seealso: `DMPLEX`, `DMPlexInterpolate()`, `DMPlexIsInterpolatedCollective()`
1976: @*/
1977: PetscErrorCode DMPlexIsInterpolated(DM dm, DMPlexInterpolatedFlag *interpolated)
1978: {
1979: DM_Plex *plex = (DM_Plex *)dm->data;
1981: PetscFunctionBegin;
1983: PetscAssertPointer(interpolated, 2);
1984: if (plex->interpolated < 0) {
1985: PetscCall(DMPlexIsInterpolated_Internal(dm, &plex->interpolated));
1986: } else if (PetscDefined(USE_DEBUG)) {
1987: DMPlexInterpolatedFlag flg;
1989: PetscCall(DMPlexIsInterpolated_Internal(dm, &flg));
1990: 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]);
1991: }
1992: *interpolated = plex->interpolated;
1993: PetscFunctionReturn(PETSC_SUCCESS);
1994: }
1996: /*@
1997: DMPlexIsInterpolatedCollective - Find out to what extent the `DMPLEX` is topologically interpolated (in collective manner).
1999: Collective
2001: Input Parameter:
2002: . dm - The `DMPLEX` object
2004: Output Parameter:
2005: . interpolated - Flag whether the `DM` is interpolated
2007: Level: intermediate
2009: Notes:
2010: Unlike `DMPlexIsInterpolated()`, this is collective so the results are guaranteed to be the same on all ranks.
2012: This function will return `DMPLEX_INTERPOLATED_MIXED` if the results of `DMPlexIsInterpolated()` are different on different ranks.
2014: Developer Notes:
2015: Initially, plex->interpolatedCollective = `DMPLEX_INTERPOLATED_INVALID`.
2017: If plex->interpolatedCollective == `DMPLEX_INTERPOLATED_INVALID`, this function calls `DMPlexIsInterpolated()` which sets plex->interpolated.
2018: `MPI_Allreduce()` is then called and collectively consistent flag plex->interpolatedCollective is set and returned;
2019: if plex->interpolated varies on different ranks, plex->interpolatedCollective = `DMPLEX_INTERPOLATED_MIXED`,
2020: otherwise sets plex->interpolatedCollective = plex->interpolated.
2022: If plex->interpolatedCollective != `DMPLEX_INTERPOLATED_INVALID`, this function just returns plex->interpolatedCollective.
2024: .seealso: `DMPLEX`, `DMPlexInterpolate()`, `DMPlexIsInterpolated()`
2025: @*/
2026: PetscErrorCode DMPlexIsInterpolatedCollective(DM dm, DMPlexInterpolatedFlag *interpolated)
2027: {
2028: DM_Plex *plex = (DM_Plex *)dm->data;
2029: PetscBool debug = PETSC_FALSE;
2031: PetscFunctionBegin;
2033: PetscAssertPointer(interpolated, 2);
2034: PetscCall(PetscOptionsGetBool(((PetscObject)dm)->options, ((PetscObject)dm)->prefix, "-dm_plex_is_interpolated_collective_debug", &debug, NULL));
2035: if (plex->interpolatedCollective < 0) {
2036: DMPlexInterpolatedFlag min, max;
2037: MPI_Comm comm;
2039: PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
2040: PetscCall(DMPlexIsInterpolated(dm, &plex->interpolatedCollective));
2041: PetscCallMPI(MPIU_Allreduce(&plex->interpolatedCollective, &min, 1, MPIU_ENUM, MPI_MIN, comm));
2042: PetscCallMPI(MPIU_Allreduce(&plex->interpolatedCollective, &max, 1, MPIU_ENUM, MPI_MAX, comm));
2043: if (min != max) plex->interpolatedCollective = DMPLEX_INTERPOLATED_MIXED;
2044: if (debug) {
2045: PetscMPIInt rank;
2047: PetscCallMPI(MPI_Comm_rank(comm, &rank));
2048: PetscCall(PetscSynchronizedPrintf(comm, "[%d] interpolated=%s interpolatedCollective=%s\n", rank, DMPlexInterpolatedFlags[plex->interpolated], DMPlexInterpolatedFlags[plex->interpolatedCollective]));
2049: PetscCall(PetscSynchronizedFlush(comm, PETSC_STDOUT));
2050: }
2051: }
2052: *interpolated = plex->interpolatedCollective;
2053: PetscFunctionReturn(PETSC_SUCCESS);
2054: }