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: PetscDisableStaticAnalyzerForExpressionUnderstandingThatThisIsDangerousAndBugprone(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: PetscDisableStaticAnalyzerForExpressionUnderstandingThatThisIsDangerousAndBugprone(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, pStart, 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, &pStart, &Np);
421: DMPlexSetChart(idm, pStart, 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: {
529: const PetscInt *cone;
530: PetscInt coneSize, ornt;
532: DMPlexGetConeSize(idm, f, &coneSize);
533: DMPlexGetCone(idm, f, &cone);
534: 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);
535: /* Notice that we have to use vertices here because the lower dimensional faces have not been created yet */
536: DMPolytopeGetVertexOrientation(faceType, cone, face, &ornt);
537: DMPlexInsertConeOrientation(idm, c, cf, ornt);
538: }
539: }
540: DMPlexRestoreRawFaces_Internal(dm, ct, cone, &numFaces, &faceTypes, &faceSizes, &faces);
541: }
542: PetscHashIJKLDestroy(&faceTable);
543: DMPlexSymmetrize(idm);
544: DMPlexStratify(idm);
545: return(0);
546: }
548: static PetscErrorCode SortRmineRremoteByRemote_Private(PetscSF sf, PetscInt *rmine1[], PetscInt *rremote1[])
549: {
550: PetscInt nleaves;
551: PetscInt nranks;
552: const PetscMPIInt *ranks=NULL;
553: const PetscInt *roffset=NULL, *rmine=NULL, *rremote=NULL;
554: PetscInt n, o, r;
555: PetscErrorCode ierr;
558: PetscSFGetRootRanks(sf, &nranks, &ranks, &roffset, &rmine, &rremote);
559: nleaves = roffset[nranks];
560: PetscMalloc2(nleaves, rmine1, nleaves, rremote1);
561: for (r=0; r<nranks; r++) {
562: /* simultaneously sort rank-wise portions of rmine & rremote by values in rremote
563: - to unify order with the other side */
564: o = roffset[r];
565: n = roffset[r+1] - o;
566: PetscArraycpy(&(*rmine1)[o], &rmine[o], n);
567: PetscArraycpy(&(*rremote1)[o], &rremote[o], n);
568: PetscSortIntWithArray(n, &(*rremote1)[o], &(*rmine1)[o]);
569: }
570: return(0);
571: }
573: PetscErrorCode DMPlexOrientInterface_Internal(DM dm)
574: {
575: PetscSF sf;
576: const PetscInt *locals;
577: const PetscSFNode *remotes;
578: const PetscMPIInt *ranks;
579: const PetscInt *roffset;
580: PetscInt *rmine1, *rremote1; /* rmine and rremote copies simultaneously sorted by rank and rremote */
581: PetscInt nroots, p, nleaves, nranks, r, maxConeSize = 0;
582: PetscInt (*roots)[4], (*leaves)[4], mainCone[4];
583: PetscMPIInt (*rootsRanks)[4], (*leavesRanks)[4];
584: MPI_Comm comm;
585: PetscMPIInt rank, size;
586: PetscInt debug = 0;
587: PetscErrorCode ierr;
590: PetscObjectGetComm((PetscObject) dm, &comm);
591: MPI_Comm_rank(comm, &rank);
592: MPI_Comm_size(comm, &size);
593: DMGetPointSF(dm, &sf);
594: DMViewFromOptions(dm, NULL, "-before_orient_interface_dm_view");
595: if (PetscDefined(USE_DEBUG)) {DMPlexCheckPointSF(dm);}
596: PetscSFGetGraph(sf, &nroots, &nleaves, &locals, &remotes);
597: if (nroots < 0) return(0);
598: PetscSFSetUp(sf);
599: SortRmineRremoteByRemote_Private(sf, &rmine1, &rremote1);
600: for (p = 0; p < nleaves; ++p) {
601: PetscInt coneSize;
602: DMPlexGetConeSize(dm, locals[p], &coneSize);
603: maxConeSize = PetscMax(maxConeSize, coneSize);
604: }
605: if (maxConeSize > 4) SETERRQ1(comm, PETSC_ERR_SUP, "This method does not support cones of size %D", maxConeSize);
606: PetscMalloc4(nroots, &roots, nroots, &leaves, nroots, &rootsRanks, nroots, &leavesRanks);
607: for (p = 0; p < nroots; ++p) {
608: const PetscInt *cone;
609: PetscInt coneSize, c, ind0;
611: DMPlexGetConeSize(dm, p, &coneSize);
612: DMPlexGetCone(dm, p, &cone);
613: /* Ignore vertices */
614: if (coneSize < 2) {
615: for (c = 0; c < 4; c++) {
616: roots[p][c] = -1;
617: rootsRanks[p][c] = -1;
618: }
619: continue;
620: }
621: /* Translate all points to root numbering */
622: for (c = 0; c < PetscMin(coneSize, 4); c++) {
623: PetscFindInt(cone[c], nleaves, locals, &ind0);
624: if (ind0 < 0) {
625: roots[p][c] = cone[c];
626: rootsRanks[p][c] = rank;
627: } else {
628: roots[p][c] = remotes[ind0].index;
629: rootsRanks[p][c] = remotes[ind0].rank;
630: }
631: }
632: for (c = coneSize; c < 4; c++) {
633: roots[p][c] = -1;
634: rootsRanks[p][c] = -1;
635: }
636: }
637: for (p = 0; p < nroots; ++p) {
638: PetscInt c;
639: for (c = 0; c < 4; c++) {
640: leaves[p][c] = -2;
641: leavesRanks[p][c] = -2;
642: }
643: }
644: PetscSFBcastBegin(sf, MPIU_4INT, roots, leaves, MPI_REPLACE);
645: PetscSFBcastBegin(sf, MPI_4INT, rootsRanks, leavesRanks, MPI_REPLACE);
646: PetscSFBcastEnd(sf, MPIU_4INT, roots, leaves, MPI_REPLACE);
647: PetscSFBcastEnd(sf, MPI_4INT, rootsRanks, leavesRanks, MPI_REPLACE);
648: if (debug) {
649: PetscSynchronizedFlush(comm, NULL);
650: if (!rank) {PetscSynchronizedPrintf(comm, "Referenced roots\n");}
651: }
652: PetscSFGetRootRanks(sf, &nranks, &ranks, &roffset, NULL, NULL);
653: for (p = 0; p < nroots; ++p) {
654: DMPolytopeType ct;
655: const PetscInt *cone;
656: PetscInt coneSize, c, ind0, o;
658: if (leaves[p][0] < 0) continue; /* Ignore vertices */
659: DMPlexGetCellType(dm, p, &ct);
660: DMPlexGetConeSize(dm, p, &coneSize);
661: DMPlexGetCone(dm, p, &cone);
662: if (debug) {
663: PetscSynchronizedPrintf(comm, "[%d] %4D: cone=[%4D %4D %4D %4D] roots=[(%d,%4D) (%d,%4D) (%d,%4D) (%d,%4D)] leaves=[(%d,%4D) (%d,%4D) (%d,%4D) (%d,%4D)]",
664: rank, p, cone[0], cone[1], cone[2], cone[3],
665: rootsRanks[p][0], roots[p][0], rootsRanks[p][1], roots[p][1], rootsRanks[p][2], roots[p][2], rootsRanks[p][3], roots[p][3],
666: leavesRanks[p][0], leaves[p][0], leavesRanks[p][1], leaves[p][1], leavesRanks[p][2], leaves[p][2], leavesRanks[p][3], leaves[p][3]);
667: }
668: if (leavesRanks[p][0] != rootsRanks[p][0] || leaves[p][0] != roots[p][0] ||
669: leavesRanks[p][1] != rootsRanks[p][1] || leaves[p][1] != roots[p][1] ||
670: leavesRanks[p][2] != rootsRanks[p][2] || leaves[p][2] != roots[p][2] ||
671: leavesRanks[p][3] != rootsRanks[p][3] || leaves[p][3] != roots[p][3]) {
672: /* Translate these leaves to my cone points; mainCone means desired order p's cone points */
673: for (c = 0; c < PetscMin(coneSize, 4); ++c) {
674: PetscInt rS, rN;
676: if (leavesRanks[p][c] == rank) {
677: /* A local leaf is just taken as it is */
678: mainCone[c] = leaves[p][c];
679: continue;
680: }
681: /* Find index of rank leavesRanks[p][c] among remote ranks */
682: /* No need for PetscMPIIntCast because these integers were originally cast from PetscMPIInt. */
683: PetscFindMPIInt((PetscMPIInt) leavesRanks[p][c], nranks, ranks, &r);
684: if (PetscUnlikely(r < 0)) SETERRQ7(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Point %D cone[%D]=%D root (%d,%D) leaf (%d,%D): leaf rank not found among remote ranks", p, c, cone[c], rootsRanks[p][c], roots[p][c], leavesRanks[p][c], leaves[p][c]);
685: 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]);
686: /* Find point leaves[p][c] among remote points aimed at rank leavesRanks[p][c] */
687: rS = roffset[r];
688: rN = roffset[r+1] - rS;
689: PetscFindInt(leaves[p][c], rN, &rremote1[rS], &ind0);
690: 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]);
691: /* Get the corresponding local point */
692: mainCone[c] = rmine1[rS + ind0];
693: }
694: if (debug) {PetscSynchronizedPrintf(comm, " mainCone=[%4D %4D %4D %4D]\n", mainCone[0], mainCone[1], mainCone[2], mainCone[3]);}
695: /* Set the desired order of p's cone points and fix orientations accordingly */
696: DMPolytopeGetOrientation(ct, cone, mainCone, &o);
697: DMPlexOrientPoint(dm, p, o);
698: } else if (debug) {PetscSynchronizedPrintf(comm, " ==\n");}
699: }
700: if (debug) {
701: PetscSynchronizedFlush(comm, NULL);
702: MPI_Barrier(comm);
703: }
704: DMViewFromOptions(dm, NULL, "-after_orient_interface_dm_view");
705: PetscFree4(roots, leaves, rootsRanks, leavesRanks);
706: PetscFree2(rmine1, rremote1);
707: return(0);
708: }
710: static PetscErrorCode IntArrayViewFromOptions(MPI_Comm comm, const char opt[], const char name[], const char idxname[], const char valname[], PetscInt n, const PetscInt a[])
711: {
712: PetscInt idx;
713: PetscMPIInt rank;
714: PetscBool flg;
718: PetscOptionsHasName(NULL, NULL, opt, &flg);
719: if (!flg) return(0);
720: MPI_Comm_rank(comm, &rank);
721: PetscSynchronizedPrintf(comm, "[%d]%s:\n", rank, name);
722: for (idx = 0; idx < n; ++idx) {PetscSynchronizedPrintf(comm, "[%d]%s %D %s %D\n", rank, idxname, idx, valname, a[idx]);}
723: PetscSynchronizedFlush(comm, NULL);
724: return(0);
725: }
727: static PetscErrorCode SFNodeArrayViewFromOptions(MPI_Comm comm, const char opt[], const char name[], const char idxname[], PetscInt n, const PetscSFNode a[])
728: {
729: PetscInt idx;
730: PetscMPIInt rank;
731: PetscBool flg;
735: PetscOptionsHasName(NULL, NULL, opt, &flg);
736: if (!flg) return(0);
737: MPI_Comm_rank(comm, &rank);
738: PetscSynchronizedPrintf(comm, "[%d]%s:\n", rank, name);
739: if (idxname) {
740: 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);}
741: } else {
742: for (idx = 0; idx < n; ++idx) {PetscSynchronizedPrintf(comm, "[%d]rank %D index %D\n", rank, a[idx].rank, a[idx].index);}
743: }
744: PetscSynchronizedFlush(comm, NULL);
745: return(0);
746: }
748: static PetscErrorCode DMPlexMapToLocalPoint(DM dm, PetscHMapIJ remotehash, PetscSFNode remotePoint, PetscInt *localPoint)
749: {
750: PetscSF sf;
751: const PetscInt *locals;
752: PetscMPIInt rank;
753: PetscErrorCode ierr;
756: MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank);
757: DMGetPointSF(dm, &sf);
758: PetscSFGetGraph(sf, NULL, NULL, &locals, NULL);
759: if (remotePoint.rank == rank) {
760: *localPoint = remotePoint.index;
761: } else {
762: PetscHashIJKey key;
763: PetscInt l;
765: key.i = remotePoint.index;
766: key.j = remotePoint.rank;
767: PetscHMapIJGet(remotehash, key, &l);
768: if (l >= 0) {
769: *localPoint = locals[l];
770: } else PetscFunctionReturn(1);
771: }
772: return(0);
773: }
775: static PetscErrorCode DMPlexMapToGlobalPoint(DM dm, PetscInt localPoint, PetscSFNode *remotePoint)
776: {
777: PetscSF sf;
778: const PetscInt *locals, *rootdegree;
779: const PetscSFNode *remotes;
780: PetscInt Nl, l;
781: PetscMPIInt rank;
782: PetscErrorCode ierr;
785: MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank);
786: DMGetPointSF(dm, &sf);
787: PetscSFGetGraph(sf, NULL, &Nl, &locals, &remotes);
788: if (Nl < 0) goto owned;
789: PetscSFComputeDegreeBegin(sf, &rootdegree);
790: PetscSFComputeDegreeEnd(sf, &rootdegree);
791: if (rootdegree[localPoint]) goto owned;
792: PetscFindInt(localPoint, Nl, locals, &l);
793: if (l < 0) PetscFunctionReturn(1);
794: *remotePoint = remotes[l];
795: return(0);
796: owned:
797: remotePoint->rank = rank;
798: remotePoint->index = localPoint;
799: return(0);
800: }
802: static PetscErrorCode DMPlexPointIsShared(DM dm, PetscInt p, PetscBool *isShared)
803: {
804: PetscSF sf;
805: const PetscInt *locals, *rootdegree;
806: PetscInt Nl, idx;
807: PetscErrorCode ierr;
810: *isShared = PETSC_FALSE;
811: DMGetPointSF(dm, &sf);
812: PetscSFGetGraph(sf, NULL, &Nl, &locals, NULL);
813: if (Nl < 0) return(0);
814: PetscFindInt(p, Nl, locals, &idx);
815: if (idx >= 0) {*isShared = PETSC_TRUE; return(0);}
816: PetscSFComputeDegreeBegin(sf, &rootdegree);
817: PetscSFComputeDegreeEnd(sf, &rootdegree);
818: if (rootdegree[p] > 0) *isShared = PETSC_TRUE;
819: return(0);
820: }
822: static PetscErrorCode DMPlexConeIsShared(DM dm, PetscInt p, PetscBool *isShared)
823: {
824: const PetscInt *cone;
825: PetscInt coneSize, c;
826: PetscBool cShared = PETSC_TRUE;
827: PetscErrorCode ierr;
830: DMPlexGetConeSize(dm, p, &coneSize);
831: DMPlexGetCone(dm, p, &cone);
832: for (c = 0; c < coneSize; ++c) {
833: PetscBool pointShared;
835: DMPlexPointIsShared(dm, cone[c], &pointShared);
836: cShared = (PetscBool) (cShared && pointShared);
837: }
838: *isShared = coneSize ? cShared : PETSC_FALSE;
839: return(0);
840: }
842: static PetscErrorCode DMPlexGetConeMinimum(DM dm, PetscInt p, PetscSFNode *cpmin)
843: {
844: const PetscInt *cone;
845: PetscInt coneSize, c;
846: PetscSFNode cmin = {PETSC_MAX_INT, PETSC_MAX_INT}, missing = {-1, -1};
847: PetscErrorCode ierr;
850: DMPlexGetConeSize(dm, p, &coneSize);
851: DMPlexGetCone(dm, p, &cone);
852: for (c = 0; c < coneSize; ++c) {
853: PetscSFNode rcp;
855: DMPlexMapToGlobalPoint(dm, cone[c], &rcp);
856: if (ierr) {
857: cmin = missing;
858: } else {
859: cmin = (rcp.rank < cmin.rank) || (rcp.rank == cmin.rank && rcp.index < cmin.index) ? rcp : cmin;
860: }
861: }
862: *cpmin = coneSize ? cmin : missing;
863: return(0);
864: }
866: /*
867: Each shared face has an entry in the candidates array:
868: (-1, coneSize-1), {(global cone point)}
869: where the set is missing the point p which we use as the key for the face
870: */
871: static PetscErrorCode DMPlexAddSharedFace_Private(DM dm, PetscSection candidateSection, PetscSFNode candidates[], PetscHMapIJ faceHash, PetscInt p, PetscBool debug)
872: {
873: MPI_Comm comm;
874: const PetscInt *support;
875: PetscInt supportSize, s, off = 0, idx = 0, overlap, cellHeight, height;
876: PetscMPIInt rank;
877: PetscErrorCode ierr;
880: PetscObjectGetComm((PetscObject) dm, &comm);
881: MPI_Comm_rank(comm, &rank);
882: DMPlexGetOverlap(dm, &overlap);
883: DMPlexGetVTKCellHeight(dm, &cellHeight);
884: DMPlexGetPointHeight(dm, p, &height);
885: if (!overlap && height <= cellHeight+1) {
886: /* cells can't be shared for non-overlapping meshes */
887: if (debug) {PetscSynchronizedPrintf(comm, "[%d] Skipping face %D to avoid adding cell to hashmap since this is nonoverlapping mesh\n", rank, p);}
888: return(0);
889: }
890: DMPlexGetSupportSize(dm, p, &supportSize);
891: DMPlexGetSupport(dm, p, &support);
892: if (candidates) {PetscSectionGetOffset(candidateSection, p, &off);}
893: for (s = 0; s < supportSize; ++s) {
894: const PetscInt face = support[s];
895: const PetscInt *cone;
896: PetscSFNode cpmin={-1,-1}, rp={-1,-1};
897: PetscInt coneSize, c, f;
898: PetscBool isShared = PETSC_FALSE;
899: PetscHashIJKey key;
901: /* Only add point once */
902: if (debug) {PetscSynchronizedPrintf(comm, "[%d] Support face %D\n", rank, face);}
903: key.i = p;
904: key.j = face;
905: PetscHMapIJGet(faceHash, key, &f);
906: if (f >= 0) continue;
907: DMPlexConeIsShared(dm, face, &isShared);
908: DMPlexGetConeMinimum(dm, face, &cpmin);
909: DMPlexMapToGlobalPoint(dm, p, &rp);
910: if (debug) {
911: PetscSynchronizedPrintf(comm, "[%d] Face point %D is shared: %d\n", rank, face, (int) isShared);
912: PetscSynchronizedPrintf(comm, "[%d] Global point (%D, %D) Min Cone Point (%D, %D)\n", rank, rp.rank, rp.index, cpmin.rank, cpmin.index);
913: }
914: if (isShared && (rp.rank == cpmin.rank && rp.index == cpmin.index)) {
915: PetscHMapIJSet(faceHash, key, p);
916: if (candidates) {
917: if (debug) {PetscSynchronizedPrintf(comm, "[%d] Adding shared face %D at idx %D\n[%d] ", rank, face, idx, rank);}
918: DMPlexGetConeSize(dm, face, &coneSize);
919: DMPlexGetCone(dm, face, &cone);
920: candidates[off+idx].rank = -1;
921: candidates[off+idx++].index = coneSize-1;
922: candidates[off+idx].rank = rank;
923: candidates[off+idx++].index = face;
924: for (c = 0; c < coneSize; ++c) {
925: const PetscInt cp = cone[c];
927: if (cp == p) continue;
928: DMPlexMapToGlobalPoint(dm, cp, &candidates[off+idx]);
929: if (debug) {PetscSynchronizedPrintf(comm, " (%D,%D)", candidates[off+idx].rank, candidates[off+idx].index);}
930: ++idx;
931: }
932: if (debug) {PetscSynchronizedPrintf(comm, "\n");}
933: } else {
934: /* Add cone size to section */
935: if (debug) {PetscSynchronizedPrintf(comm, "[%d] Scheduling shared face %D\n", rank, face);}
936: DMPlexGetConeSize(dm, face, &coneSize);
937: PetscHMapIJSet(faceHash, key, p);
938: PetscSectionAddDof(candidateSection, p, coneSize+1);
939: }
940: }
941: }
942: return(0);
943: }
945: /*@
946: DMPlexInterpolatePointSF - Insert interpolated points in the overlap into the PointSF in parallel, following local interpolation
948: Collective on dm
950: Input Parameters:
951: + dm - The interpolated DM
952: - pointSF - The initial SF without interpolated points
954: Output Parameter:
955: . pointSF - The SF including interpolated points
957: Level: developer
959: 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
961: .seealso: DMPlexInterpolate(), DMPlexUninterpolate()
962: @*/
963: PetscErrorCode DMPlexInterpolatePointSF(DM dm, PetscSF pointSF)
964: {
965: MPI_Comm comm;
966: PetscHMapIJ remoteHash;
967: PetscHMapI claimshash;
968: PetscSection candidateSection, candidateRemoteSection, claimSection;
969: PetscSFNode *candidates, *candidatesRemote, *claims;
970: const PetscInt *localPoints, *rootdegree;
971: const PetscSFNode *remotePoints;
972: PetscInt ov, Nr, r, Nl, l;
973: PetscInt candidatesSize, candidatesRemoteSize, claimsSize;
974: PetscBool flg, debug = PETSC_FALSE;
975: PetscMPIInt rank;
976: PetscErrorCode ierr;
981: DMPlexIsDistributed(dm, &flg);
982: if (!flg) return(0);
983: /* Set initial SF so that lower level queries work */
984: DMSetPointSF(dm, pointSF);
985: PetscObjectGetComm((PetscObject) dm, &comm);
986: MPI_Comm_rank(comm, &rank);
987: DMPlexGetOverlap(dm, &ov);
988: if (ov) SETERRQ(comm, PETSC_ERR_SUP, "Interpolation of overlapped DMPlex not implemented yet");
989: PetscOptionsHasName(NULL, ((PetscObject) dm)->prefix, "-dmplex_interp_debug", &debug);
990: PetscObjectViewFromOptions((PetscObject) dm, NULL, "-dm_interp_pre_view");
991: PetscObjectViewFromOptions((PetscObject) pointSF, NULL, "-petscsf_interp_pre_view");
992: PetscLogEventBegin(DMPLEX_InterpolateSF,dm,0,0,0);
993: /* Step 0: Precalculations */
994: PetscSFGetGraph(pointSF, &Nr, &Nl, &localPoints, &remotePoints);
995: if (Nr < 0) SETERRQ(comm, PETSC_ERR_ARG_WRONGSTATE, "This DMPlex is distributed but input PointSF has no graph set");
996: PetscHMapIJCreate(&remoteHash);
997: for (l = 0; l < Nl; ++l) {
998: PetscHashIJKey key;
999: key.i = remotePoints[l].index;
1000: key.j = remotePoints[l].rank;
1001: PetscHMapIJSet(remoteHash, key, l);
1002: }
1003: /* Compute root degree to identify shared points */
1004: PetscSFComputeDegreeBegin(pointSF, &rootdegree);
1005: PetscSFComputeDegreeEnd(pointSF, &rootdegree);
1006: IntArrayViewFromOptions(comm, "-interp_root_degree_view", "Root degree", "point", "degree", Nr, rootdegree);
1007: /*
1008: 1) Loop over each leaf point $p$ at depth $d$ in the SF
1009: \item Get set $F(p)$ of faces $f$ in the support of $p$ for which
1010: \begin{itemize}
1011: \item all cone points of $f$ are shared
1012: \item $p$ is the cone point with smallest canonical number
1013: \end{itemize}
1014: \item Send $F(p)$ and the cone of each face to the active root point $r(p)$
1015: \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
1016: \item Send the root face from the root back to all leaf process
1017: \item Leaf processes add the shared face to the SF
1018: */
1019: /* Step 1: Construct section+SFNode array
1020: The section has entries for all shared faces for which we have a leaf point in the cone
1021: The array holds candidate shared faces, each face is refered to by the leaf point */
1022: PetscSectionCreate(comm, &candidateSection);
1023: PetscSectionSetChart(candidateSection, 0, Nr);
1024: {
1025: PetscHMapIJ faceHash;
1027: PetscHMapIJCreate(&faceHash);
1028: for (l = 0; l < Nl; ++l) {
1029: const PetscInt p = localPoints[l];
1031: if (debug) {PetscSynchronizedPrintf(comm, "[%d] First pass leaf point %D\n", rank, p);}
1032: DMPlexAddSharedFace_Private(dm, candidateSection, NULL, faceHash, p, debug);
1033: }
1034: PetscHMapIJClear(faceHash);
1035: PetscSectionSetUp(candidateSection);
1036: PetscSectionGetStorageSize(candidateSection, &candidatesSize);
1037: PetscMalloc1(candidatesSize, &candidates);
1038: for (l = 0; l < Nl; ++l) {
1039: const PetscInt p = localPoints[l];
1041: if (debug) {PetscSynchronizedPrintf(comm, "[%d] Second pass leaf point %D\n", rank, p);}
1042: DMPlexAddSharedFace_Private(dm, candidateSection, candidates, faceHash, p, debug);
1043: }
1044: PetscHMapIJDestroy(&faceHash);
1045: if (debug) {PetscSynchronizedFlush(comm, NULL);}
1046: }
1047: PetscObjectSetName((PetscObject) candidateSection, "Candidate Section");
1048: PetscObjectViewFromOptions((PetscObject) candidateSection, NULL, "-petscsection_interp_candidate_view");
1049: SFNodeArrayViewFromOptions(comm, "-petscsection_interp_candidate_view", "Candidates", NULL, candidatesSize, candidates);
1050: /* Step 2: Gather candidate section / array pair into the root partition via inverse(multi(pointSF)). */
1051: /* Note that this section is indexed by offsets into leaves, not by point number */
1052: {
1053: PetscSF sfMulti, sfInverse, sfCandidates;
1054: PetscInt *remoteOffsets;
1056: PetscSFGetMultiSF(pointSF, &sfMulti);
1057: PetscSFCreateInverseSF(sfMulti, &sfInverse);
1058: PetscSectionCreate(comm, &candidateRemoteSection);
1059: PetscSFDistributeSection(sfInverse, candidateSection, &remoteOffsets, candidateRemoteSection);
1060: PetscSFCreateSectionSF(sfInverse, candidateSection, remoteOffsets, candidateRemoteSection, &sfCandidates);
1061: PetscSectionGetStorageSize(candidateRemoteSection, &candidatesRemoteSize);
1062: PetscMalloc1(candidatesRemoteSize, &candidatesRemote);
1063: PetscSFBcastBegin(sfCandidates, MPIU_2INT, candidates, candidatesRemote,MPI_REPLACE);
1064: PetscSFBcastEnd(sfCandidates, MPIU_2INT, candidates, candidatesRemote,MPI_REPLACE);
1065: PetscSFDestroy(&sfInverse);
1066: PetscSFDestroy(&sfCandidates);
1067: PetscFree(remoteOffsets);
1069: PetscObjectSetName((PetscObject) candidateRemoteSection, "Remote Candidate Section");
1070: PetscObjectViewFromOptions((PetscObject) candidateRemoteSection, NULL, "-petscsection_interp_candidate_remote_view");
1071: SFNodeArrayViewFromOptions(comm, "-petscsection_interp_candidate_remote_view", "Remote Candidates", NULL, candidatesRemoteSize, candidatesRemote);
1072: }
1073: /* 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 */
1074: {
1075: PetscHashIJKLRemote faceTable;
1076: PetscInt idx, idx2;
1078: PetscHashIJKLRemoteCreate(&faceTable);
1079: /* There is a section point for every leaf attached to a given root point */
1080: for (r = 0, idx = 0, idx2 = 0; r < Nr; ++r) {
1081: PetscInt deg;
1083: for (deg = 0; deg < rootdegree[r]; ++deg, ++idx) {
1084: PetscInt offset, dof, d;
1086: PetscSectionGetDof(candidateRemoteSection, idx, &dof);
1087: PetscSectionGetOffset(candidateRemoteSection, idx, &offset);
1088: /* dof may include many faces from the remote process */
1089: for (d = 0; d < dof; ++d) {
1090: const PetscInt hidx = offset+d;
1091: const PetscInt Np = candidatesRemote[hidx].index+1;
1092: const PetscSFNode rface = candidatesRemote[hidx+1];
1093: const PetscSFNode *fcone = &candidatesRemote[hidx+2];
1094: PetscSFNode fcp0;
1095: const PetscSFNode pmax = {PETSC_MAX_INT, PETSC_MAX_INT};
1096: const PetscInt *join = NULL;
1097: PetscHashIJKLRemoteKey key;
1098: PetscHashIter iter;
1099: PetscBool missing;
1100: PetscInt points[1024], p, joinSize;
1102: 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);}
1103: 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);
1104: fcp0.rank = rank;
1105: fcp0.index = r;
1106: d += Np;
1107: /* Put remote face in hash table */
1108: key.i = fcp0;
1109: key.j = fcone[0];
1110: key.k = Np > 2 ? fcone[1] : pmax;
1111: key.l = Np > 3 ? fcone[2] : pmax;
1112: PetscSortSFNode(Np, (PetscSFNode *) &key);
1113: PetscHashIJKLRemotePut(faceTable, key, &iter, &missing);
1114: if (missing) {
1115: if (debug) {PetscSynchronizedPrintf(PetscObjectComm((PetscObject) dm), "[%d] Setting remote face (%D, %D)\n", rank, rface.index, rface.rank);}
1116: PetscHashIJKLRemoteIterSet(faceTable, iter, rface);
1117: } else {
1118: PetscSFNode oface;
1120: PetscHashIJKLRemoteIterGet(faceTable, iter, &oface);
1121: if ((rface.rank < oface.rank) || (rface.rank == oface.rank && rface.index < oface.index)) {
1122: if (debug) {PetscSynchronizedPrintf(PetscObjectComm((PetscObject) dm), "[%d] Replacing with remote face (%D, %D)\n", rank, rface.index, rface.rank);}
1123: PetscHashIJKLRemoteIterSet(faceTable, iter, rface);
1124: }
1125: }
1126: /* Check for local face */
1127: points[0] = r;
1128: for (p = 1; p < Np; ++p) {
1129: DMPlexMapToLocalPoint(dm, remoteHash, fcone[p-1], &points[p]);
1130: if (ierr) break; /* We got a point not in our overlap */
1131: if (debug) {PetscSynchronizedPrintf(PetscObjectComm((PetscObject) dm), "[%d] Checking local candidate %D\n", rank, points[p]);}
1132: }
1133: if (ierr) continue;
1134: DMPlexGetJoin(dm, Np, points, &joinSize, &join);
1135: if (joinSize == 1) {
1136: PetscSFNode lface;
1137: PetscSFNode oface;
1139: /* Always replace with local face */
1140: lface.rank = rank;
1141: lface.index = join[0];
1142: PetscHashIJKLRemoteIterGet(faceTable, iter, &oface);
1143: 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);}
1144: PetscHashIJKLRemoteIterSet(faceTable, iter, lface);
1145: }
1146: DMPlexRestoreJoin(dm, Np, points, &joinSize, &join);
1147: }
1148: }
1149: /* Put back faces for this root */
1150: for (deg = 0; deg < rootdegree[r]; ++deg, ++idx2) {
1151: PetscInt offset, dof, d;
1153: PetscSectionGetDof(candidateRemoteSection, idx2, &dof);
1154: PetscSectionGetOffset(candidateRemoteSection, idx2, &offset);
1155: /* dof may include many faces from the remote process */
1156: for (d = 0; d < dof; ++d) {
1157: const PetscInt hidx = offset+d;
1158: const PetscInt Np = candidatesRemote[hidx].index+1;
1159: const PetscSFNode *fcone = &candidatesRemote[hidx+2];
1160: PetscSFNode fcp0;
1161: const PetscSFNode pmax = {PETSC_MAX_INT, PETSC_MAX_INT};
1162: PetscHashIJKLRemoteKey key;
1163: PetscHashIter iter;
1164: PetscBool missing;
1166: if (debug) {PetscSynchronizedPrintf(PetscObjectComm((PetscObject) dm), "[%d] Entering face at (%D, %D)\n", rank, r, idx);}
1167: if (Np > 4) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot handle faces with %D cone points", Np);
1168: fcp0.rank = rank;
1169: fcp0.index = r;
1170: d += Np;
1171: /* Find remote face in hash table */
1172: key.i = fcp0;
1173: key.j = fcone[0];
1174: key.k = Np > 2 ? fcone[1] : pmax;
1175: key.l = Np > 3 ? fcone[2] : pmax;
1176: PetscSortSFNode(Np, (PetscSFNode *) &key);
1177: 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);}
1178: PetscHashIJKLRemotePut(faceTable, key, &iter, &missing);
1179: if (missing) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Root %D Idx %D ought to have an associated face", r, idx2);
1180: else {PetscHashIJKLRemoteIterGet(faceTable, iter, &candidatesRemote[hidx]);}
1181: }
1182: }
1183: }
1184: if (debug) {PetscSynchronizedFlush(PetscObjectComm((PetscObject) dm), NULL);}
1185: PetscHashIJKLRemoteDestroy(&faceTable);
1186: }
1187: /* Step 4: Push back owned faces */
1188: {
1189: PetscSF sfMulti, sfClaims, sfPointNew;
1190: PetscSFNode *remotePointsNew;
1191: PetscInt *remoteOffsets, *localPointsNew;
1192: PetscInt pStart, pEnd, r, NlNew, p;
1194: /* 4) Push claims back to receiver via the MultiSF and derive new pointSF mapping on receiver */
1195: PetscSFGetMultiSF(pointSF, &sfMulti);
1196: PetscSectionCreate(comm, &claimSection);
1197: PetscSFDistributeSection(sfMulti, candidateRemoteSection, &remoteOffsets, claimSection);
1198: PetscSFCreateSectionSF(sfMulti, candidateRemoteSection, remoteOffsets, claimSection, &sfClaims);
1199: PetscSectionGetStorageSize(claimSection, &claimsSize);
1200: PetscMalloc1(claimsSize, &claims);
1201: for (p = 0; p < claimsSize; ++p) claims[p].rank = -1;
1202: PetscSFBcastBegin(sfClaims, MPIU_2INT, candidatesRemote, claims,MPI_REPLACE);
1203: PetscSFBcastEnd(sfClaims, MPIU_2INT, candidatesRemote, claims,MPI_REPLACE);
1204: PetscSFDestroy(&sfClaims);
1205: PetscFree(remoteOffsets);
1206: PetscObjectSetName((PetscObject) claimSection, "Claim Section");
1207: PetscObjectViewFromOptions((PetscObject) claimSection, NULL, "-petscsection_interp_claim_view");
1208: SFNodeArrayViewFromOptions(comm, "-petscsection_interp_claim_view", "Claims", NULL, claimsSize, claims);
1209: /* Step 5) Walk the original section of local supports and add an SF entry for each updated item */
1210: /* TODO I should not have to do a join here since I already put the face and its cone in the candidate section */
1211: PetscHMapICreate(&claimshash);
1212: for (r = 0; r < Nr; ++r) {
1213: PetscInt dof, off, d;
1215: if (debug) {PetscSynchronizedPrintf(comm, "[%d] Checking root for claims %D\n", rank, r);}
1216: PetscSectionGetDof(candidateSection, r, &dof);
1217: PetscSectionGetOffset(candidateSection, r, &off);
1218: for (d = 0; d < dof;) {
1219: if (claims[off+d].rank >= 0) {
1220: const PetscInt faceInd = off+d;
1221: const PetscInt Np = candidates[off+d].index;
1222: const PetscInt *join = NULL;
1223: PetscInt joinSize, points[1024], c;
1225: if (debug) {PetscSynchronizedPrintf(comm, "[%d] Found claim for remote point (%D, %D)\n", rank, claims[faceInd].rank, claims[faceInd].index);}
1226: points[0] = r;
1227: if (debug) {PetscSynchronizedPrintf(comm, "[%d] point %D\n", rank, points[0]);}
1228: for (c = 0, d += 2; c < Np; ++c, ++d) {
1229: DMPlexMapToLocalPoint(dm, remoteHash, candidates[off+d], &points[c+1]);
1230: if (debug) {PetscSynchronizedPrintf(comm, "[%d] point %D\n", rank, points[c+1]);}
1231: }
1232: DMPlexGetJoin(dm, Np+1, points, &joinSize, &join);
1233: if (joinSize == 1) {
1234: if (claims[faceInd].rank == rank) {
1235: if (debug) {PetscSynchronizedPrintf(comm, "[%d] Ignoring local face %D for non-remote partner\n", rank, join[0]);}
1236: } else {
1237: if (debug) {PetscSynchronizedPrintf(comm, "[%d] Found local face %D\n", rank, join[0]);}
1238: PetscHMapISet(claimshash, join[0], faceInd);
1239: }
1240: } else {
1241: if (debug) {PetscSynchronizedPrintf(comm, "[%d] Failed to find face\n", rank);}
1242: }
1243: DMPlexRestoreJoin(dm, Np+1, points, &joinSize, &join);
1244: } else {
1245: if (debug) {PetscSynchronizedPrintf(comm, "[%d] No claim for point %D\n", rank, r);}
1246: d += claims[off+d].index+1;
1247: }
1248: }
1249: }
1250: if (debug) {PetscSynchronizedFlush(comm, NULL);}
1251: /* Step 6) Create new pointSF from hashed claims */
1252: PetscHMapIGetSize(claimshash, &NlNew);
1253: DMPlexGetChart(dm, &pStart, &pEnd);
1254: PetscMalloc1(Nl + NlNew, &localPointsNew);
1255: PetscMalloc1(Nl + NlNew, &remotePointsNew);
1256: for (l = 0; l < Nl; ++l) {
1257: localPointsNew[l] = localPoints[l];
1258: remotePointsNew[l].index = remotePoints[l].index;
1259: remotePointsNew[l].rank = remotePoints[l].rank;
1260: }
1261: p = Nl;
1262: PetscHMapIGetKeys(claimshash, &p, localPointsNew);
1263: /* We sort new points, and assume they are numbered after all existing points */
1264: PetscSortInt(NlNew, &localPointsNew[Nl]);
1265: for (p = Nl; p < Nl + NlNew; ++p) {
1266: PetscInt off;
1267: PetscHMapIGet(claimshash, localPointsNew[p], &off);
1268: 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);
1269: remotePointsNew[p] = claims[off];
1270: }
1271: PetscSFCreate(comm, &sfPointNew);
1272: PetscSFSetGraph(sfPointNew, pEnd-pStart, Nl+NlNew, localPointsNew, PETSC_OWN_POINTER, remotePointsNew, PETSC_OWN_POINTER);
1273: PetscSFSetUp(sfPointNew);
1274: DMSetPointSF(dm, sfPointNew);
1275: PetscObjectViewFromOptions((PetscObject) sfPointNew, NULL, "-petscsf_interp_view");
1276: PetscSFDestroy(&sfPointNew);
1277: PetscHMapIDestroy(&claimshash);
1278: }
1279: PetscHMapIJDestroy(&remoteHash);
1280: PetscSectionDestroy(&candidateSection);
1281: PetscSectionDestroy(&candidateRemoteSection);
1282: PetscSectionDestroy(&claimSection);
1283: PetscFree(candidates);
1284: PetscFree(candidatesRemote);
1285: PetscFree(claims);
1286: PetscLogEventEnd(DMPLEX_InterpolateSF,dm,0,0,0);
1287: return(0);
1288: }
1290: /*@
1291: DMPlexInterpolate - Take in a cell-vertex mesh and return one with all intermediate faces, edges, etc.
1293: Collective on dm
1295: Input Parameters:
1296: + dm - The DMPlex object with only cells and vertices
1297: - dmInt - The interpolated DM
1299: Output Parameter:
1300: . dmInt - The complete DMPlex object
1302: Level: intermediate
1304: Notes:
1305: It does not copy over the coordinates.
1307: Developer Notes:
1308: It sets plex->interpolated = DMPLEX_INTERPOLATED_FULL.
1310: .seealso: DMPlexUninterpolate(), DMPlexCreateFromCellListPetsc(), DMPlexCopyCoordinates()
1311: @*/
1312: PetscErrorCode DMPlexInterpolate(DM dm, DM *dmInt)
1313: {
1314: DMPlexInterpolatedFlag interpolated;
1315: DM idm, odm = dm;
1316: PetscSF sfPoint;
1317: PetscInt depth, dim, d;
1318: const char *name;
1319: PetscBool flg=PETSC_TRUE;
1325: PetscLogEventBegin(DMPLEX_Interpolate,dm,0,0,0);
1326: DMPlexGetDepth(dm, &depth);
1327: DMGetDimension(dm, &dim);
1328: DMPlexIsInterpolated(dm, &interpolated);
1329: if (interpolated == DMPLEX_INTERPOLATED_PARTIAL) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Not for partially interpolated meshes");
1330: if (interpolated == DMPLEX_INTERPOLATED_FULL) {
1331: PetscObjectReference((PetscObject) dm);
1332: idm = dm;
1333: } else {
1334: for (d = 1; d < dim; ++d) {
1335: /* Create interpolated mesh */
1336: DMCreate(PetscObjectComm((PetscObject)dm), &idm);
1337: DMSetType(idm, DMPLEX);
1338: DMSetDimension(idm, dim);
1339: if (depth > 0) {
1340: DMPlexInterpolateFaces_Internal(odm, 1, idm);
1341: DMGetPointSF(odm, &sfPoint);
1342: {
1343: /* TODO: We need to systematically fix cases of distributed Plexes with no graph set */
1344: PetscInt nroots;
1345: PetscSFGetGraph(sfPoint, &nroots, NULL, NULL, NULL);
1346: if (nroots >= 0) {DMPlexInterpolatePointSF(idm, sfPoint);}
1347: }
1348: }
1349: if (odm != dm) {DMDestroy(&odm);}
1350: odm = idm;
1351: }
1352: PetscObjectGetName((PetscObject) dm, &name);
1353: PetscObjectSetName((PetscObject) idm, name);
1354: DMPlexCopyCoordinates(dm, idm);
1355: DMCopyLabels(dm, idm, PETSC_COPY_VALUES, PETSC_FALSE);
1356: PetscOptionsGetBool(((PetscObject)dm)->options, ((PetscObject)dm)->prefix, "-dm_plex_interpolate_orient_interfaces", &flg, NULL);
1357: if (flg) {DMPlexOrientInterface_Internal(idm);}
1358: }
1359: {
1360: PetscBool isper;
1361: const PetscReal *maxCell, *L;
1362: const DMBoundaryType *bd;
1364: DMGetPeriodicity(dm,&isper,&maxCell,&L,&bd);
1365: DMSetPeriodicity(idm,isper,maxCell,L,bd);
1366: }
1367: /* This function makes the mesh fully interpolated on all ranks */
1368: {
1369: DM_Plex *plex = (DM_Plex *) idm->data;
1370: plex->interpolated = plex->interpolatedCollective = DMPLEX_INTERPOLATED_FULL;
1371: }
1372: *dmInt = idm;
1373: PetscLogEventEnd(DMPLEX_Interpolate,dm,0,0,0);
1374: return(0);
1375: }
1377: /*@
1378: DMPlexCopyCoordinates - Copy coordinates from one mesh to another with the same vertices
1380: Collective on dmA
1382: Input Parameter:
1383: . dmA - The DMPlex object with initial coordinates
1385: Output Parameter:
1386: . dmB - The DMPlex object with copied coordinates
1388: Level: intermediate
1390: Note: This is typically used when adding pieces other than vertices to a mesh
1392: .seealso: DMCopyLabels(), DMGetCoordinates(), DMGetCoordinatesLocal(), DMGetCoordinateDM(), DMGetCoordinateSection()
1393: @*/
1394: PetscErrorCode DMPlexCopyCoordinates(DM dmA, DM dmB)
1395: {
1396: Vec coordinatesA, coordinatesB;
1397: VecType vtype;
1398: PetscSection coordSectionA, coordSectionB;
1399: PetscScalar *coordsA, *coordsB;
1400: PetscInt spaceDim, Nf, vStartA, vStartB, vEndA, vEndB, coordSizeB, v, d;
1401: PetscInt cStartA, cEndA, cStartB, cEndB, cS, cE, cdim;
1402: PetscBool lc = PETSC_FALSE;
1408: if (dmA == dmB) return(0);
1409: DMGetCoordinateDim(dmA, &cdim);
1410: DMSetCoordinateDim(dmB, cdim);
1411: DMPlexGetDepthStratum(dmA, 0, &vStartA, &vEndA);
1412: DMPlexGetDepthStratum(dmB, 0, &vStartB, &vEndB);
1413: 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);
1414: /* Copy over discretization if it exists */
1415: {
1416: DM cdmA, cdmB;
1417: PetscDS dsA, dsB;
1418: PetscObject objA, objB;
1419: PetscClassId idA, idB;
1420: const PetscScalar *constants;
1421: PetscInt cdim, Nc;
1423: DMGetCoordinateDM(dmA, &cdmA);
1424: DMGetCoordinateDM(dmB, &cdmB);
1425: DMGetField(cdmA, 0, NULL, &objA);
1426: DMGetField(cdmB, 0, NULL, &objB);
1427: PetscObjectGetClassId(objA, &idA);
1428: PetscObjectGetClassId(objB, &idB);
1429: if ((idA == PETSCFE_CLASSID) && (idA != idB)) {
1430: DMSetField(cdmB, 0, NULL, objA);
1431: DMCreateDS(cdmB);
1432: DMGetDS(cdmA, &dsA);
1433: DMGetDS(cdmB, &dsB);
1434: PetscDSGetCoordinateDimension(dsA, &cdim);
1435: PetscDSSetCoordinateDimension(dsB, cdim);
1436: PetscDSGetConstants(dsA, &Nc, &constants);
1437: PetscDSSetConstants(dsB, Nc, (PetscScalar *) constants);
1438: }
1439: }
1440: DMPlexGetHeightStratum(dmA, 0, &cStartA, &cEndA);
1441: DMPlexGetHeightStratum(dmB, 0, &cStartB, &cEndB);
1442: DMGetCoordinateSection(dmA, &coordSectionA);
1443: DMGetCoordinateSection(dmB, &coordSectionB);
1444: if (coordSectionA == coordSectionB) return(0);
1445: PetscSectionGetNumFields(coordSectionA, &Nf);
1446: if (!Nf) return(0);
1447: if (Nf > 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "The number of coordinate fields must be 1, not %D", Nf);
1448: if (!coordSectionB) {
1449: PetscInt dim;
1451: PetscSectionCreate(PetscObjectComm((PetscObject) coordSectionA), &coordSectionB);
1452: DMGetCoordinateDim(dmA, &dim);
1453: DMSetCoordinateSection(dmB, dim, coordSectionB);
1454: PetscObjectDereference((PetscObject) coordSectionB);
1455: }
1456: PetscSectionSetNumFields(coordSectionB, 1);
1457: PetscSectionGetFieldComponents(coordSectionA, 0, &spaceDim);
1458: PetscSectionSetFieldComponents(coordSectionB, 0, spaceDim);
1459: PetscSectionGetChart(coordSectionA, &cS, &cE);
1460: if (cStartA <= cS && cS < cEndA) { /* localized coordinates */
1461: 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);
1462: cS = cS - cStartA + cStartB;
1463: cE = vEndB;
1464: lc = PETSC_TRUE;
1465: } else {
1466: cS = vStartB;
1467: cE = vEndB;
1468: }
1469: PetscSectionSetChart(coordSectionB, cS, cE);
1470: for (v = vStartB; v < vEndB; ++v) {
1471: PetscSectionSetDof(coordSectionB, v, spaceDim);
1472: PetscSectionSetFieldDof(coordSectionB, v, 0, spaceDim);
1473: }
1474: if (lc) { /* localized coordinates */
1475: PetscInt c;
1477: for (c = cS-cStartB; c < cEndB-cStartB; c++) {
1478: PetscInt dof;
1480: PetscSectionGetDof(coordSectionA, c + cStartA, &dof);
1481: PetscSectionSetDof(coordSectionB, c + cStartB, dof);
1482: PetscSectionSetFieldDof(coordSectionB, c + cStartB, 0, dof);
1483: }
1484: }
1485: PetscSectionSetUp(coordSectionB);
1486: PetscSectionGetStorageSize(coordSectionB, &coordSizeB);
1487: DMGetCoordinatesLocal(dmA, &coordinatesA);
1488: VecCreate(PETSC_COMM_SELF, &coordinatesB);
1489: PetscObjectSetName((PetscObject) coordinatesB, "coordinates");
1490: VecSetSizes(coordinatesB, coordSizeB, PETSC_DETERMINE);
1491: VecGetBlockSize(coordinatesA, &d);
1492: VecSetBlockSize(coordinatesB, d);
1493: VecGetType(coordinatesA, &vtype);
1494: VecSetType(coordinatesB, vtype);
1495: VecGetArray(coordinatesA, &coordsA);
1496: VecGetArray(coordinatesB, &coordsB);
1497: for (v = 0; v < vEndB-vStartB; ++v) {
1498: PetscInt offA, offB;
1500: PetscSectionGetOffset(coordSectionA, v + vStartA, &offA);
1501: PetscSectionGetOffset(coordSectionB, v + vStartB, &offB);
1502: for (d = 0; d < spaceDim; ++d) {
1503: coordsB[offB+d] = coordsA[offA+d];
1504: }
1505: }
1506: if (lc) { /* localized coordinates */
1507: PetscInt c;
1509: for (c = cS-cStartB; c < cEndB-cStartB; c++) {
1510: PetscInt dof, offA, offB;
1512: PetscSectionGetOffset(coordSectionA, c + cStartA, &offA);
1513: PetscSectionGetOffset(coordSectionB, c + cStartB, &offB);
1514: PetscSectionGetDof(coordSectionA, c + cStartA, &dof);
1515: PetscArraycpy(coordsB + offB,coordsA + offA,dof);
1516: }
1517: }
1518: VecRestoreArray(coordinatesA, &coordsA);
1519: VecRestoreArray(coordinatesB, &coordsB);
1520: DMSetCoordinatesLocal(dmB, coordinatesB);
1521: VecDestroy(&coordinatesB);
1522: return(0);
1523: }
1525: /*@
1526: DMPlexUninterpolate - Take in a mesh with all intermediate faces, edges, etc. and return a cell-vertex mesh
1528: Collective on dm
1530: Input Parameter:
1531: . dm - The complete DMPlex object
1533: Output Parameter:
1534: . dmUnint - The DMPlex object with only cells and vertices
1536: Level: intermediate
1538: Notes:
1539: It does not copy over the coordinates.
1541: Developer Notes:
1542: It sets plex->interpolated = DMPLEX_INTERPOLATED_NONE.
1544: .seealso: DMPlexInterpolate(), DMPlexCreateFromCellListPetsc(), DMPlexCopyCoordinates()
1545: @*/
1546: PetscErrorCode DMPlexUninterpolate(DM dm, DM *dmUnint)
1547: {
1548: DMPlexInterpolatedFlag interpolated;
1549: DM udm;
1550: PetscInt dim, vStart, vEnd, cStart, cEnd, c, maxConeSize = 0, *cone;
1556: DMGetDimension(dm, &dim);
1557: DMPlexIsInterpolated(dm, &interpolated);
1558: if (interpolated == DMPLEX_INTERPOLATED_PARTIAL) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Not for partially interpolated meshes");
1559: if (interpolated == DMPLEX_INTERPOLATED_NONE || dim <= 1) {
1560: /* in case dim <= 1 just keep the DMPLEX_INTERPOLATED_FULL flag */
1561: PetscObjectReference((PetscObject) dm);
1562: *dmUnint = dm;
1563: return(0);
1564: }
1565: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
1566: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
1567: DMCreate(PetscObjectComm((PetscObject) dm), &udm);
1568: DMSetType(udm, DMPLEX);
1569: DMSetDimension(udm, dim);
1570: DMPlexSetChart(udm, cStart, vEnd);
1571: for (c = cStart; c < cEnd; ++c) {
1572: PetscInt *closure = NULL, closureSize, cl, coneSize = 0;
1574: DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
1575: for (cl = 0; cl < closureSize*2; cl += 2) {
1576: const PetscInt p = closure[cl];
1578: if ((p >= vStart) && (p < vEnd)) ++coneSize;
1579: }
1580: DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
1581: DMPlexSetConeSize(udm, c, coneSize);
1582: maxConeSize = PetscMax(maxConeSize, coneSize);
1583: }
1584: DMSetUp(udm);
1585: PetscMalloc1(maxConeSize, &cone);
1586: for (c = cStart; c < cEnd; ++c) {
1587: PetscInt *closure = NULL, closureSize, cl, coneSize = 0;
1589: DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
1590: for (cl = 0; cl < closureSize*2; cl += 2) {
1591: const PetscInt p = closure[cl];
1593: if ((p >= vStart) && (p < vEnd)) cone[coneSize++] = p;
1594: }
1595: DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
1596: DMPlexSetCone(udm, c, cone);
1597: }
1598: PetscFree(cone);
1599: DMPlexSymmetrize(udm);
1600: DMPlexStratify(udm);
1601: /* Reduce SF */
1602: {
1603: PetscSF sfPoint, sfPointUn;
1604: const PetscSFNode *remotePoints;
1605: const PetscInt *localPoints;
1606: PetscSFNode *remotePointsUn;
1607: PetscInt *localPointsUn;
1608: PetscInt vEnd, numRoots, numLeaves, l;
1609: PetscInt numLeavesUn = 0, n = 0;
1610: PetscErrorCode ierr;
1612: /* Get original SF information */
1613: DMGetPointSF(dm, &sfPoint);
1614: DMGetPointSF(udm, &sfPointUn);
1615: DMPlexGetDepthStratum(dm, 0, NULL, &vEnd);
1616: PetscSFGetGraph(sfPoint, &numRoots, &numLeaves, &localPoints, &remotePoints);
1617: /* Allocate space for cells and vertices */
1618: for (l = 0; l < numLeaves; ++l) if (localPoints[l] < vEnd) numLeavesUn++;
1619: /* Fill in leaves */
1620: if (vEnd >= 0) {
1621: PetscMalloc1(numLeavesUn, &remotePointsUn);
1622: PetscMalloc1(numLeavesUn, &localPointsUn);
1623: for (l = 0; l < numLeaves; l++) {
1624: if (localPoints[l] < vEnd) {
1625: localPointsUn[n] = localPoints[l];
1626: remotePointsUn[n].rank = remotePoints[l].rank;
1627: remotePointsUn[n].index = remotePoints[l].index;
1628: ++n;
1629: }
1630: }
1631: if (n != numLeavesUn) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Inconsistent number of leaves %d != %d", n, numLeavesUn);
1632: PetscSFSetGraph(sfPointUn, vEnd, numLeavesUn, localPointsUn, PETSC_OWN_POINTER, remotePointsUn, PETSC_OWN_POINTER);
1633: }
1634: }
1635: {
1636: PetscBool isper;
1637: const PetscReal *maxCell, *L;
1638: const DMBoundaryType *bd;
1640: DMGetPeriodicity(dm,&isper,&maxCell,&L,&bd);
1641: DMSetPeriodicity(udm,isper,maxCell,L,bd);
1642: }
1643: /* This function makes the mesh fully uninterpolated on all ranks */
1644: {
1645: DM_Plex *plex = (DM_Plex *) udm->data;
1646: plex->interpolated = plex->interpolatedCollective = DMPLEX_INTERPOLATED_NONE;
1647: }
1648: *dmUnint = udm;
1649: return(0);
1650: }
1652: static PetscErrorCode DMPlexIsInterpolated_Internal(DM dm, DMPlexInterpolatedFlag *interpolated)
1653: {
1654: PetscInt coneSize, depth, dim, h, p, pStart, pEnd;
1655: MPI_Comm comm;
1659: PetscObjectGetComm((PetscObject)dm, &comm);
1660: DMPlexGetDepth(dm, &depth);
1661: DMGetDimension(dm, &dim);
1663: if (depth == dim) {
1664: *interpolated = DMPLEX_INTERPOLATED_FULL;
1665: if (!dim) goto finish;
1667: /* Check points at height = dim are vertices (have no cones) */
1668: DMPlexGetHeightStratum(dm, dim, &pStart, &pEnd);
1669: for (p=pStart; p<pEnd; p++) {
1670: DMPlexGetConeSize(dm, p, &coneSize);
1671: if (coneSize) {
1672: *interpolated = DMPLEX_INTERPOLATED_PARTIAL;
1673: goto finish;
1674: }
1675: }
1677: /* Check points at height < dim have cones */
1678: for (h=0; h<dim; h++) {
1679: DMPlexGetHeightStratum(dm, h, &pStart, &pEnd);
1680: for (p=pStart; p<pEnd; p++) {
1681: DMPlexGetConeSize(dm, p, &coneSize);
1682: if (!coneSize) {
1683: *interpolated = DMPLEX_INTERPOLATED_PARTIAL;
1684: goto finish;
1685: }
1686: }
1687: }
1688: } else if (depth == 1) {
1689: *interpolated = DMPLEX_INTERPOLATED_NONE;
1690: } else {
1691: *interpolated = DMPLEX_INTERPOLATED_PARTIAL;
1692: }
1693: finish:
1694: return(0);
1695: }
1697: /*@
1698: DMPlexIsInterpolated - Find out to what extent the DMPlex is topologically interpolated.
1700: Not Collective
1702: Input Parameter:
1703: . dm - The DM object
1705: Output Parameter:
1706: . interpolated - Flag whether the DM is interpolated
1708: Level: intermediate
1710: Notes:
1711: Unlike DMPlexIsInterpolatedCollective(), this is NOT collective
1712: so the results can be different on different ranks in special cases.
1713: However, DMPlexInterpolate() guarantees the result is the same on all.
1715: Unlike DMPlexIsInterpolatedCollective(), this cannot return DMPLEX_INTERPOLATED_MIXED.
1717: Developer Notes:
1718: Initially, plex->interpolated = DMPLEX_INTERPOLATED_INVALID.
1720: If plex->interpolated == DMPLEX_INTERPOLATED_INVALID, DMPlexIsInterpolated_Internal() is called.
1721: It checks the actual topology and sets plex->interpolated on each rank separately to one of
1722: DMPLEX_INTERPOLATED_NONE, DMPLEX_INTERPOLATED_PARTIAL or DMPLEX_INTERPOLATED_FULL.
1724: If plex->interpolated != DMPLEX_INTERPOLATED_INVALID, this function just returns plex->interpolated.
1726: DMPlexInterpolate() sets plex->interpolated = DMPLEX_INTERPOLATED_FULL,
1727: and DMPlexUninterpolate() sets plex->interpolated = DMPLEX_INTERPOLATED_NONE.
1729: .seealso: DMPlexInterpolate(), DMPlexIsInterpolatedCollective()
1730: @*/
1731: PetscErrorCode DMPlexIsInterpolated(DM dm, DMPlexInterpolatedFlag *interpolated)
1732: {
1733: DM_Plex *plex = (DM_Plex *) dm->data;
1734: PetscErrorCode ierr;
1739: if (plex->interpolated < 0) {
1740: DMPlexIsInterpolated_Internal(dm, &plex->interpolated);
1741: } else if (PetscDefined (USE_DEBUG)) {
1742: DMPlexInterpolatedFlag flg;
1744: DMPlexIsInterpolated_Internal(dm, &flg);
1745: if (flg != plex->interpolated) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Stashed DMPlexInterpolatedFlag %s is inconsistent with current %s", DMPlexInterpolatedFlags[plex->interpolated], DMPlexInterpolatedFlags[flg]);
1746: }
1747: *interpolated = plex->interpolated;
1748: return(0);
1749: }
1751: /*@
1752: DMPlexIsInterpolatedCollective - Find out to what extent the DMPlex is topologically interpolated (in collective manner).
1754: Collective
1756: Input Parameter:
1757: . dm - The DM object
1759: Output Parameter:
1760: . interpolated - Flag whether the DM is interpolated
1762: Level: intermediate
1764: Notes:
1765: Unlike DMPlexIsInterpolated(), this is collective so the results are guaranteed to be the same on all ranks.
1767: This function will return DMPLEX_INTERPOLATED_MIXED if the results of DMPlexIsInterpolated() are different on different ranks.
1769: Developer Notes:
1770: Initially, plex->interpolatedCollective = DMPLEX_INTERPOLATED_INVALID.
1772: If plex->interpolatedCollective == DMPLEX_INTERPOLATED_INVALID, this function calls DMPlexIsInterpolated() which sets plex->interpolated.
1773: MPI_Allreduce() is then called and collectively consistent flag plex->interpolatedCollective is set and returned;
1774: if plex->interpolated varies on different ranks, plex->interpolatedCollective = DMPLEX_INTERPOLATED_MIXED,
1775: otherwise sets plex->interpolatedCollective = plex->interpolated.
1777: If plex->interpolatedCollective != DMPLEX_INTERPOLATED_INVALID, this function just returns plex->interpolatedCollective.
1779: .seealso: DMPlexInterpolate(), DMPlexIsInterpolated()
1780: @*/
1781: PetscErrorCode DMPlexIsInterpolatedCollective(DM dm, DMPlexInterpolatedFlag *interpolated)
1782: {
1783: DM_Plex *plex = (DM_Plex *) dm->data;
1784: PetscBool debug=PETSC_FALSE;
1785: PetscErrorCode ierr;
1790: PetscOptionsGetBool(((PetscObject) dm)->options, ((PetscObject) dm)->prefix, "-dm_plex_is_interpolated_collective_debug", &debug, NULL);
1791: if (plex->interpolatedCollective < 0) {
1792: DMPlexInterpolatedFlag min, max;
1793: MPI_Comm comm;
1795: PetscObjectGetComm((PetscObject)dm, &comm);
1796: DMPlexIsInterpolated(dm, &plex->interpolatedCollective);
1797: MPI_Allreduce(&plex->interpolatedCollective, &min, 1, MPIU_ENUM, MPI_MIN, comm);
1798: MPI_Allreduce(&plex->interpolatedCollective, &max, 1, MPIU_ENUM, MPI_MAX, comm);
1799: if (min != max) plex->interpolatedCollective = DMPLEX_INTERPOLATED_MIXED;
1800: if (debug) {
1801: PetscMPIInt rank;
1803: MPI_Comm_rank(comm, &rank);
1804: PetscSynchronizedPrintf(comm, "[%d] interpolated=%s interpolatedCollective=%s\n", rank, DMPlexInterpolatedFlags[plex->interpolated], DMPlexInterpolatedFlags[plex->interpolatedCollective]);
1805: PetscSynchronizedFlush(comm, PETSC_STDOUT);
1806: }
1807: }
1808: *interpolated = plex->interpolatedCollective;
1809: return(0);
1810: }