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