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