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