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