Actual source code: plexinterpolate.c
petsc-3.10.5 2019-03-28
1: #include <petsc/private/dmpleximpl.h>
2: #include <petsc/private/hashmapi.h>
3: #include <petsc/private/hashmapij.h>
5: /* HashIJKL */
7: #include <petsc/private/hashmap.h>
9: typedef struct _PetscHashIJKLKey { PetscInt i, j, k, l; } PetscHashIJKLKey;
11: #define PetscHashIJKLKeyHash(key) \
12: PetscHashCombine(PetscHashCombine(PetscHashInt((key).i),PetscHashInt((key).j)), \
13: PetscHashCombine(PetscHashInt((key).k),PetscHashInt((key).l)))
15: #define PetscHashIJKLKeyEqual(k1,k2) \
16: (((k1).i==(k2).i) ? ((k1).j==(k2).j) ? ((k1).k==(k2).k) ? ((k1).l==(k2).l) : 0 : 0 : 0)
18: PETSC_HASH_MAP(HashIJKL, PetscHashIJKLKey, PetscInt, PetscHashIJKLKeyHash, PetscHashIJKLKeyEqual, -1)
21: /*
22: DMPlexGetFaces_Internal - Gets groups of vertices that correspond to faces for the given cell
23: This assumes that the mesh is not interpolated from the depth of point p to the vertices
24: */
25: PetscErrorCode DMPlexGetFaces_Internal(DM dm, PetscInt dim, PetscInt p, PetscInt *numFaces, PetscInt *faceSize, const PetscInt *faces[])
26: {
27: const PetscInt *cone = NULL;
28: PetscInt coneSize;
29: PetscErrorCode ierr;
33: DMPlexGetConeSize(dm, p, &coneSize);
34: DMPlexGetCone(dm, p, &cone);
35: DMPlexGetRawFaces_Internal(dm, dim, coneSize, cone, numFaces, faceSize, faces);
36: return(0);
37: }
39: /*
40: DMPlexRestoreFaces_Internal - Restores the array
41: */
42: PetscErrorCode DMPlexRestoreFaces_Internal(DM dm, PetscInt dim, PetscInt p, PetscInt *numFaces, PetscInt *faceSize, const PetscInt *faces[])
43: {
44: PetscErrorCode ierr;
47: if (faces) { DMRestoreWorkArray(dm, 0, MPIU_INT, (void *) faces); }
48: return(0);
49: }
51: /*
52: DMPlexGetRawFaces_Internal - Gets groups of vertices that correspond to faces for the given cone
53: */
54: PetscErrorCode DMPlexGetRawFaces_Internal(DM dm, PetscInt dim, PetscInt coneSize, const PetscInt cone[], PetscInt *numFaces, PetscInt *faceSize, const PetscInt *faces[])
55: {
56: PetscInt *facesTmp;
57: PetscInt maxConeSize, maxSupportSize;
58: PetscErrorCode ierr;
63: DMPlexGetMaxSizes(dm, &maxConeSize, &maxSupportSize);
64: if (faces) {DMGetWorkArray(dm, PetscSqr(PetscMax(maxConeSize, maxSupportSize)), MPIU_INT, &facesTmp);}
65: switch (dim) {
66: case 1:
67: switch (coneSize) {
68: case 2:
69: if (faces) {
70: facesTmp[0] = cone[0]; facesTmp[1] = cone[1];
71: *faces = facesTmp;
72: }
73: if (numFaces) *numFaces = 2;
74: if (faceSize) *faceSize = 1;
75: break;
76: default:
77: SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Cone size %D not supported for dimension %D", coneSize, dim);
78: }
79: break;
80: case 2:
81: switch (coneSize) {
82: case 3:
83: if (faces) {
84: facesTmp[0] = cone[0]; facesTmp[1] = cone[1];
85: facesTmp[2] = cone[1]; facesTmp[3] = cone[2];
86: facesTmp[4] = cone[2]; facesTmp[5] = cone[0];
87: *faces = facesTmp;
88: }
89: if (numFaces) *numFaces = 3;
90: if (faceSize) *faceSize = 2;
91: break;
92: case 4:
93: /* Vertices follow right hand rule */
94: if (faces) {
95: facesTmp[0] = cone[0]; facesTmp[1] = cone[1];
96: facesTmp[2] = cone[1]; facesTmp[3] = cone[2];
97: facesTmp[4] = cone[2]; facesTmp[5] = cone[3];
98: facesTmp[6] = cone[3]; facesTmp[7] = cone[0];
99: *faces = facesTmp;
100: }
101: if (numFaces) *numFaces = 4;
102: if (faceSize) *faceSize = 2;
103: break;
104: default:
105: SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Cone size %D not supported for dimension %D", coneSize, dim);
106: }
107: break;
108: case 3:
109: switch (coneSize) {
110: case 3:
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: if (numFaces) *numFaces = 3;
118: if (faceSize) *faceSize = 2;
119: break;
120: case 4:
121: /* Vertices of first face follow right hand rule and normal points away from last vertex */
122: if (faces) {
123: facesTmp[0] = cone[0]; facesTmp[1] = cone[1]; facesTmp[2] = cone[2];
124: facesTmp[3] = cone[0]; facesTmp[4] = cone[3]; facesTmp[5] = cone[1];
125: facesTmp[6] = cone[0]; facesTmp[7] = cone[2]; facesTmp[8] = cone[3];
126: facesTmp[9] = cone[2]; facesTmp[10] = cone[1]; facesTmp[11] = cone[3];
127: *faces = facesTmp;
128: }
129: if (numFaces) *numFaces = 4;
130: if (faceSize) *faceSize = 3;
131: break;
132: case 8:
133: /* 7--------6
134: /| /|
135: / | / |
136: 4--------5 |
137: | | | |
138: | | | |
139: | 1--------2
140: | / | /
141: |/ |/
142: 0--------3
143: */
144: if (faces) {
145: facesTmp[0] = cone[0]; facesTmp[1] = cone[1]; facesTmp[2] = cone[2]; facesTmp[3] = cone[3]; /* Bottom */
146: facesTmp[4] = cone[4]; facesTmp[5] = cone[5]; facesTmp[6] = cone[6]; facesTmp[7] = cone[7]; /* Top */
147: facesTmp[8] = cone[0]; facesTmp[9] = cone[3]; facesTmp[10] = cone[5]; facesTmp[11] = cone[4]; /* Front */
148: facesTmp[12] = cone[2]; facesTmp[13] = cone[1]; facesTmp[14] = cone[7]; facesTmp[15] = cone[6]; /* Back */
149: facesTmp[16] = cone[3]; facesTmp[17] = cone[2]; facesTmp[18] = cone[6]; facesTmp[19] = cone[5]; /* Right */
150: facesTmp[20] = cone[0]; facesTmp[21] = cone[4]; facesTmp[22] = cone[7]; facesTmp[23] = cone[1]; /* Left */
151: *faces = facesTmp;
152: }
153: if (numFaces) *numFaces = 6;
154: if (faceSize) *faceSize = 4;
155: break;
156: default:
157: SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Cone size %D not supported for dimension %D", coneSize, dim);
158: }
159: break;
160: default:
161: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Dimension %D not supported", dim);
162: }
163: return(0);
164: }
166: /*
167: DMPlexGetRawFacesHybrid_Internal - Gets groups of vertices that correspond to faces for the given cone using hybrid ordering (prisms)
168: */
169: static PetscErrorCode DMPlexGetRawFacesHybrid_Internal(DM dm, PetscInt dim, PetscInt coneSize, const PetscInt cone[], PetscInt *numFaces, PetscInt *numFacesNotH, PetscInt *faceSize, const PetscInt *faces[])
170: {
171: PetscInt *facesTmp;
172: PetscInt maxConeSize, maxSupportSize;
173: PetscErrorCode ierr;
178: DMPlexGetMaxSizes(dm, &maxConeSize, &maxSupportSize);
179: if (faces) {DMGetWorkArray(dm, PetscSqr(PetscMax(maxConeSize, maxSupportSize)), MPIU_INT, &facesTmp);}
180: switch (dim) {
181: case 1:
182: switch (coneSize) {
183: case 2:
184: if (faces) {
185: facesTmp[0] = cone[0]; facesTmp[1] = cone[1];
186: *faces = facesTmp;
187: }
188: if (numFaces) *numFaces = 2;
189: if (numFacesNotH) *numFacesNotH = 2;
190: if (faceSize) *faceSize = 1;
191: break;
192: default:
193: SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Cone size %D not supported for dimension %D", coneSize, dim);
194: }
195: break;
196: case 2:
197: switch (coneSize) {
198: case 4:
199: if (faces) {
200: facesTmp[0] = cone[0]; facesTmp[1] = cone[1];
201: facesTmp[2] = cone[2]; facesTmp[3] = cone[3];
202: facesTmp[4] = cone[0]; facesTmp[5] = cone[2];
203: facesTmp[6] = cone[1]; facesTmp[7] = cone[3];
204: *faces = facesTmp;
205: }
206: if (numFaces) *numFaces = 4;
207: if (numFacesNotH) *numFacesNotH = 2;
208: if (faceSize) *faceSize = 2;
209: break;
210: default:
211: SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Cone size %D not supported for dimension %D", coneSize, dim);
212: }
213: break;
214: case 3:
215: switch (coneSize) {
216: case 6: /* triangular prism */
217: if (faces) {
218: facesTmp[0] = cone[0]; facesTmp[1] = cone[1]; facesTmp[2] = cone[2]; facesTmp[3] = -1; /* Bottom */
219: facesTmp[4] = cone[3]; facesTmp[5] = cone[4]; facesTmp[6] = cone[5]; facesTmp[7] = -1; /* Top */
220: facesTmp[8] = cone[0]; facesTmp[9] = cone[1]; facesTmp[10] = cone[3]; facesTmp[11] = cone[4]; /* Back left */
221: facesTmp[12] = cone[1]; facesTmp[13] = cone[2]; facesTmp[14] = cone[4]; facesTmp[15] = cone[5]; /* Back right */
222: facesTmp[16] = cone[2]; facesTmp[17] = cone[0]; facesTmp[18] = cone[5]; facesTmp[19] = cone[3]; /* Front */
223: *faces = facesTmp;
224: }
225: if (numFaces) *numFaces = 5;
226: if (numFacesNotH) *numFacesNotH = 2;
227: if (faceSize) *faceSize = -4;
228: break;
229: default:
230: SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Cone size %D not supported for dimension %D", coneSize, dim);
231: }
232: break;
233: default:
234: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Dimension %D not supported", dim);
235: }
236: return(0);
237: }
239: static PetscErrorCode DMPlexRestoreRawFacesHybrid_Internal(DM dm, PetscInt dim, PetscInt coneSize, const PetscInt cone[], PetscInt *numFaces, PetscInt *numFacesNotH, PetscInt *faceSize, const PetscInt *faces[])
240: {
241: PetscErrorCode ierr;
244: if (faces) { DMRestoreWorkArray(dm, 0, MPIU_INT, (void *) faces); }
245: return(0);
246: }
248: static PetscErrorCode DMPlexGetFacesHybrid_Internal(DM dm, PetscInt dim, PetscInt p, PetscInt *numFaces, PetscInt *numFacesNotH, PetscInt *faceSize, const PetscInt *faces[])
249: {
250: const PetscInt *cone = NULL;
251: PetscInt coneSize;
252: PetscErrorCode ierr;
256: DMPlexGetConeSize(dm, p, &coneSize);
257: DMPlexGetCone(dm, p, &cone);
258: DMPlexGetRawFacesHybrid_Internal(dm, dim, coneSize, cone, numFaces, numFacesNotH, faceSize, faces);
259: return(0);
260: }
262: /* This interpolates faces for cells at some stratum */
263: static PetscErrorCode DMPlexInterpolateFaces_Internal(DM dm, PetscInt cellDepth, DM idm)
264: {
265: DMLabel subpointMap;
266: PetscHashIJKL faceTable;
267: PetscInt *pStart, *pEnd;
268: PetscInt cellDim, depth, faceDepth = cellDepth, numPoints = 0, faceSizeAll = 0, face, c, d;
269: PetscInt coneSizeH = 0, faceSizeAllH = 0, numCellFacesH = 0, faceH, pMax = -1, dim, outerloop;
270: PetscInt cMax, fMax, eMax, vMax;
274: DMGetDimension(dm, &cellDim);
275: /* HACK: I need a better way to determine face dimension, or an alternative to GetFaces() */
276: DMPlexGetSubpointMap(dm, &subpointMap);
277: if (subpointMap) ++cellDim;
278: DMPlexGetDepth(dm, &depth);
279: ++depth;
280: ++cellDepth;
281: cellDim -= depth - cellDepth;
282: PetscMalloc2(depth+1,&pStart,depth+1,&pEnd);
283: for (d = depth-1; d >= faceDepth; --d) {
284: DMPlexGetDepthStratum(dm, d, &pStart[d+1], &pEnd[d+1]);
285: }
286: DMPlexGetDepthStratum(dm, -1, NULL, &pStart[faceDepth]);
287: pEnd[faceDepth] = pStart[faceDepth];
288: for (d = faceDepth-1; d >= 0; --d) {
289: DMPlexGetDepthStratum(dm, d, &pStart[d], &pEnd[d]);
290: }
291: cMax = fMax = eMax = vMax = PETSC_DETERMINE;
292: DMGetDimension(dm, &dim);
293: if (cellDim == dim) {
294: DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);
295: pMax = cMax;
296: } else if (cellDim == dim -1) {
297: DMPlexGetHybridBounds(dm, &cMax, &fMax, NULL, NULL);
298: pMax = fMax;
299: }
300: pMax = pMax < 0 ? pEnd[cellDepth] : pMax;
301: if (pMax < pEnd[cellDepth]) {
302: const PetscInt *cellFaces, *cone;
303: PetscInt numCellFacesT, faceSize, cf;
305: DMPlexGetConeSize(dm, pMax, &coneSizeH);
306: DMPlexGetCone(dm, pMax, &cone);
307: DMPlexGetRawFacesHybrid_Internal(dm, cellDim, coneSizeH, cone, &numCellFacesH, &numCellFacesT, &faceSize, &cellFaces);
308: if (faceSize < 0) {
309: PetscInt *sizes, minv, maxv;
311: /* count vertices of hybrid and non-hybrid faces */
312: PetscCalloc1(numCellFacesH, &sizes);
313: for (cf = 0; cf < numCellFacesT; ++cf) { /* These are the non-hybrid faces */
314: const PetscInt *cellFace = &cellFaces[-cf*faceSize];
315: PetscInt f;
317: for (f = 0; f < -faceSize; ++f) sizes[cf] += (cellFace[f] >= 0 ? 1 : 0);
318: }
319: PetscSortInt(numCellFacesT, sizes);
320: minv = sizes[0];
321: maxv = sizes[PetscMax(numCellFacesT-1, 0)];
322: if (minv != maxv) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_SUP, "Different number of vertices for non-hybrid face %D != %D", minv, maxv);
323: faceSizeAll = minv;
324: PetscMemzero(sizes, numCellFacesH*sizeof(PetscInt));
325: for (cf = numCellFacesT; cf < numCellFacesH; ++cf) { /* These are the hybrid faces */
326: const PetscInt *cellFace = &cellFaces[-cf*faceSize];
327: PetscInt f;
329: for (f = 0; f < -faceSize; ++f) sizes[cf-numCellFacesT] += (cellFace[f] >= 0 ? 1 : 0);
330: }
331: PetscSortInt(numCellFacesH - numCellFacesT, sizes);
332: minv = sizes[0];
333: maxv = sizes[PetscMax(numCellFacesH - numCellFacesT-1, 0)];
334: PetscFree(sizes);
335: if (minv != maxv) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_SUP, "Different number of vertices for hybrid face %D != %D", minv, maxv);
336: faceSizeAllH = minv;
337: } else { /* the size of the faces in hybrid cells is the same */
338: faceSizeAll = faceSizeAllH = faceSize;
339: }
340: DMPlexRestoreRawFacesHybrid_Internal(dm, cellDim, coneSizeH, cone, &numCellFacesH, &numCellFacesT, &faceSize, &cellFaces);
341: } else if (pEnd[cellDepth] > pStart[cellDepth]) {
342: DMPlexGetFaces_Internal(dm, cellDim, pStart[cellDepth], NULL, &faceSizeAll, NULL);
343: }
344: if (faceSizeAll > 4) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Do not support interpolation of meshes with faces of %D vertices", faceSizeAll);
346: /* With hybrid grids, we first iterate on hybrid cells and start numbering the non-hybrid grids
347: Then, faces for non-hybrid cells are numbered.
348: This is to guarantee consistent orientations (all 0) of all the points in the cone of the hybrid cells */
349: PetscHashIJKLCreate(&faceTable);
350: for (outerloop = 0, face = pStart[faceDepth]; outerloop < 2; outerloop++) {
351: PetscInt start, end;
353: start = outerloop == 0 ? pMax : pStart[cellDepth];
354: end = outerloop == 0 ? pEnd[cellDepth] : pMax;
355: for (c = start; c < end; ++c) {
356: const PetscInt *cellFaces;
357: PetscInt numCellFaces, faceSize, faceSizeInc, cf;
359: if (c < pMax) {
360: DMPlexGetFaces_Internal(dm, cellDim, c, &numCellFaces, &faceSize, &cellFaces);
361: if (faceSize != faceSizeAll) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Inconsistent face for cell %D of size %D != %D", c, faceSize, faceSizeAll);
362: } else { /* Hybrid cell */
363: const PetscInt *cone;
364: PetscInt numCellFacesN, coneSize;
366: DMPlexGetConeSize(dm, c, &coneSize);
367: if (coneSize != coneSizeH) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unexpected hybrid coneSize %D != %D", coneSize, coneSizeH);
368: DMPlexGetCone(dm, c, &cone);
369: DMPlexGetRawFacesHybrid_Internal(dm, cellDim, coneSize, cone, &numCellFaces, &numCellFacesN, &faceSize, &cellFaces);
370: if (numCellFaces != numCellFacesH) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unexpected numCellFaces %D != %D for hybrid cell %D", numCellFaces, numCellFacesH, c);
371: faceSize = PetscMax(faceSize, -faceSize);
372: if (faceSize > 4) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Do not support interpolation of meshes with faces of %D vertices", faceSize);
373: numCellFaces = numCellFacesN; /* process only non-hybrid faces */
374: }
375: faceSizeInc = faceSize;
376: for (cf = 0; cf < numCellFaces; ++cf) {
377: const PetscInt *cellFace = &cellFaces[cf*faceSizeInc];
378: PetscInt faceSizeH = faceSize;
379: PetscHashIJKLKey key;
380: PetscHashIter iter;
381: PetscBool missing;
383: if (faceSizeInc == 2) {
384: key.i = PetscMin(cellFace[0], cellFace[1]);
385: key.j = PetscMax(cellFace[0], cellFace[1]);
386: key.k = PETSC_MAX_INT;
387: key.l = PETSC_MAX_INT;
388: } else {
389: key.i = cellFace[0];
390: key.j = cellFace[1];
391: key.k = cellFace[2];
392: key.l = faceSize > 3 ? (cellFace[3] < 0 ? faceSizeH = 3, PETSC_MAX_INT : cellFace[3]) : PETSC_MAX_INT;
393: PetscSortInt(faceSize, (PetscInt *) &key);
394: }
395: /* this check is redundant for non-hybrid meshes */
396: if (faceSizeH != faceSizeAll) SETERRQ4(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unexpected number of vertices for face %D of point %D -> %D != %D", cf, c, faceSizeH, faceSizeAll);
397: PetscHashIJKLPut(faceTable, key, &iter, &missing);
398: if (missing) {PetscHashIJKLIterSet(faceTable, iter, face++);}
399: }
400: if (c < pMax) {
401: DMPlexRestoreFaces_Internal(dm, cellDim, c, &numCellFaces, &faceSize, &cellFaces);
402: } else {
403: DMPlexRestoreRawFacesHybrid_Internal(dm, cellDim, coneSizeH, NULL, NULL, NULL, NULL, &cellFaces);
404: }
405: }
406: }
407: pEnd[faceDepth] = face;
409: /* Second pass for hybrid meshes: number hybrid faces */
410: for (c = pMax; c < pEnd[cellDepth]; ++c) {
411: const PetscInt *cellFaces, *cone;
412: PetscInt numCellFaces, numCellFacesN, faceSize, cf, coneSize;
414: DMPlexGetConeSize(dm, c, &coneSize);
415: DMPlexGetCone(dm, c, &cone);
416: DMPlexGetRawFacesHybrid_Internal(dm, cellDim, coneSize, cone, &numCellFaces, &numCellFacesN, &faceSize, &cellFaces);
417: if (numCellFaces != numCellFacesH) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unexpected hybrid numCellFaces %D != %D", numCellFaces, numCellFacesH);
418: faceSize = PetscMax(faceSize, -faceSize);
419: for (cf = numCellFacesN; cf < numCellFaces; ++cf) { /* These are the hybrid faces */
420: const PetscInt *cellFace = &cellFaces[cf*faceSize];
421: PetscHashIJKLKey key;
422: PetscHashIter iter;
423: PetscBool missing;
424: PetscInt faceSizeH = faceSize;
426: if (faceSize == 2) {
427: key.i = PetscMin(cellFace[0], cellFace[1]);
428: key.j = PetscMax(cellFace[0], cellFace[1]);
429: key.k = PETSC_MAX_INT;
430: key.l = PETSC_MAX_INT;
431: } else {
432: key.i = cellFace[0];
433: key.j = cellFace[1];
434: key.k = cellFace[2];
435: key.l = faceSize > 3 ? (cellFace[3] < 0 ? faceSizeH = 3, PETSC_MAX_INT : cellFace[3]) : PETSC_MAX_INT;
436: PetscSortInt(faceSize, (PetscInt *) &key);
437: }
438: if (faceSizeH != faceSizeAllH) SETERRQ4(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unexpected number of vertices for hybrid face %D of point %D -> %D != %D", cf, c, faceSizeH, faceSizeAllH);
439: PetscHashIJKLPut(faceTable, key, &iter, &missing);
440: if (missing) {PetscHashIJKLIterSet(faceTable, iter, face++);}
441: }
442: DMPlexRestoreRawFacesHybrid_Internal(dm, cellDim, coneSize, cone, &numCellFaces, &numCellFacesN, &faceSize, &cellFaces);
443: }
444: faceH = face - pEnd[faceDepth];
445: if (faceH) {
446: if (fMax == PETSC_DETERMINE) fMax = pEnd[faceDepth];
447: else if (eMax == PETSC_DETERMINE) eMax = pEnd[faceDepth];
448: else SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of unassigned hybrid facets %D for cellDim %D and dimension %D", faceH, cellDim, dim);
449: }
450: pEnd[faceDepth] = face;
451: PetscHashIJKLDestroy(&faceTable);
452: /* Count new points */
453: for (d = 0; d <= depth; ++d) {
454: numPoints += pEnd[d]-pStart[d];
455: }
456: DMPlexSetChart(idm, 0, numPoints);
457: /* Set cone sizes */
458: for (d = 0; d <= depth; ++d) {
459: PetscInt coneSize, p;
461: if (d == faceDepth) {
462: /* I see no way to do this if we admit faces of different shapes */
463: for (p = pStart[d]; p < pEnd[d]-faceH; ++p) {
464: DMPlexSetConeSize(idm, p, faceSizeAll);
465: }
466: for (p = pEnd[d]-faceH; p < pEnd[d]; ++p) {
467: DMPlexSetConeSize(idm, p, faceSizeAllH);
468: }
469: } else if (d == cellDepth) {
470: for (p = pStart[d]; p < pEnd[d]; ++p) {
471: /* Number of cell faces may be different from number of cell vertices*/
472: if (p < pMax) {
473: DMPlexGetFaces_Internal(dm, cellDim, p, &coneSize, NULL, NULL);
474: } else {
475: DMPlexGetFacesHybrid_Internal(dm, cellDim, p, &coneSize, NULL, NULL, NULL);
476: }
477: DMPlexSetConeSize(idm, p, coneSize);
478: }
479: } else {
480: for (p = pStart[d]; p < pEnd[d]; ++p) {
481: DMPlexGetConeSize(dm, p, &coneSize);
482: DMPlexSetConeSize(idm, p, coneSize);
483: }
484: }
485: }
486: DMSetUp(idm);
487: /* Get face cones from subsets of cell vertices */
488: if (faceSizeAll > 4) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Do not support interpolation of meshes with faces of %D vertices", faceSizeAll);
489: PetscHashIJKLCreate(&faceTable);
490: for (d = depth; d > cellDepth; --d) {
491: const PetscInt *cone;
492: PetscInt p;
494: for (p = pStart[d]; p < pEnd[d]; ++p) {
495: DMPlexGetCone(dm, p, &cone);
496: DMPlexSetCone(idm, p, cone);
497: DMPlexGetConeOrientation(dm, p, &cone);
498: DMPlexSetConeOrientation(idm, p, cone);
499: }
500: }
501: for (outerloop = 0, face = pStart[faceDepth]; outerloop < 2; outerloop++) {
502: PetscInt start, end;
504: start = outerloop == 0 ? pMax : pStart[cellDepth];
505: end = outerloop == 0 ? pEnd[cellDepth] : pMax;
506: for (c = start; c < end; ++c) {
507: const PetscInt *cellFaces;
508: PetscInt numCellFaces, faceSize, faceSizeInc, cf;
510: if (c < pMax) {
511: DMPlexGetFaces_Internal(dm, cellDim, c, &numCellFaces, &faceSize, &cellFaces);
512: if (faceSize != faceSizeAll) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Inconsistent face for cell %D of size %D != %D", c, faceSize, faceSizeAll);
513: } else {
514: const PetscInt *cone;
515: PetscInt numCellFacesN, coneSize;
517: DMPlexGetConeSize(dm, c, &coneSize);
518: if (coneSize != coneSizeH) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unexpected hybrid coneSize %D != %D", coneSize, coneSizeH);
519: DMPlexGetCone(dm, c, &cone);
520: DMPlexGetRawFacesHybrid_Internal(dm, cellDim, coneSize, cone, &numCellFaces, &numCellFacesN, &faceSize, &cellFaces);
521: if (numCellFaces != numCellFacesH) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unexpected numCellFaces %D != %D for hybrid cell %D", numCellFaces, numCellFacesH, c);
522: faceSize = PetscMax(faceSize, -faceSize);
523: if (faceSize > 4) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Do not support interpolation of meshes with faces of %D vertices", faceSize);
524: numCellFaces = numCellFacesN; /* process only non-hybrid faces */
525: }
526: faceSizeInc = faceSize;
527: for (cf = 0; cf < numCellFaces; ++cf) {
528: const PetscInt *cellFace = &cellFaces[cf*faceSizeInc];
529: PetscHashIJKLKey key;
530: PetscHashIter iter;
531: PetscBool missing;
533: if (faceSizeInc == 2) {
534: key.i = PetscMin(cellFace[0], cellFace[1]);
535: key.j = PetscMax(cellFace[0], cellFace[1]);
536: key.k = PETSC_MAX_INT;
537: key.l = PETSC_MAX_INT;
538: } else {
539: key.i = cellFace[0];
540: key.j = cellFace[1];
541: key.k = cellFace[2];
542: key.l = faceSizeInc > 3 ? (cellFace[3] < 0 ? faceSize = 3, PETSC_MAX_INT : cellFace[3]) : PETSC_MAX_INT;
543: PetscSortInt(faceSizeInc, (PetscInt *) &key);
544: }
545: PetscHashIJKLPut(faceTable, key, &iter, &missing);
546: if (missing) {
547: DMPlexSetCone(idm, face, cellFace);
548: PetscHashIJKLIterSet(faceTable, iter, face);
549: DMPlexInsertCone(idm, c, cf, face++);
550: } else {
551: const PetscInt *cone;
552: PetscInt coneSize, ornt, i, j, f;
554: PetscHashIJKLIterGet(faceTable, iter, &f);
555: DMPlexInsertCone(idm, c, cf, f);
556: /* Orient face: Do not allow reverse orientation at the first vertex */
557: DMPlexGetConeSize(idm, f, &coneSize);
558: DMPlexGetCone(idm, f, &cone);
559: 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);
560: /* - First find the initial vertex */
561: for (i = 0; i < faceSize; ++i) if (cellFace[0] == cone[i]) break;
562: /* - Try forward comparison */
563: for (j = 0; j < faceSize; ++j) if (cellFace[j] != cone[(i+j)%faceSize]) break;
564: if (j == faceSize) {
565: if ((faceSize == 2) && (i == 1)) ornt = -2;
566: else ornt = i;
567: } else {
568: /* - Try backward comparison */
569: for (j = 0; j < faceSize; ++j) if (cellFace[j] != cone[(i+faceSize-j)%faceSize]) break;
570: if (j == faceSize) {
571: if (i == 0) ornt = -faceSize;
572: else ornt = -i;
573: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Could not determine face orientation");
574: }
575: DMPlexInsertConeOrientation(idm, c, cf, ornt);
576: }
577: }
578: if (c < pMax) {
579: DMPlexRestoreFaces_Internal(dm, cellDim, c, &numCellFaces, &faceSize, &cellFaces);
580: } else {
581: DMPlexRestoreRawFacesHybrid_Internal(dm, cellDim, coneSizeH, NULL, NULL, NULL, NULL, &cellFaces);
582: }
583: }
584: }
585: /* Second pass for hybrid meshes: orient hybrid faces */
586: for (c = pMax; c < pEnd[cellDepth]; ++c) {
587: const PetscInt *cellFaces, *cone;
588: PetscInt numCellFaces, numCellFacesN, faceSize, cf, coneSize;
590: DMPlexGetConeSize(dm, c, &coneSize);
591: DMPlexGetCone(dm, c, &cone);
592: DMPlexGetRawFacesHybrid_Internal(dm, cellDim, coneSize, cone, &numCellFaces, &numCellFacesN, &faceSize, &cellFaces);
593: if (numCellFaces != numCellFacesH) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unexpected hybrid numCellFaces %D != %D", numCellFaces, numCellFacesH);
594: faceSize = PetscMax(faceSize, -faceSize);
595: for (cf = numCellFacesN; cf < numCellFaces; ++cf) { /* These are the hybrid faces */
596: const PetscInt *cellFace = &cellFaces[cf*faceSize];
597: PetscHashIJKLKey key;
598: PetscHashIter iter;
599: PetscBool missing;
600: PetscInt faceSizeH = faceSize;
602: if (faceSize == 2) {
603: key.i = PetscMin(cellFace[0], cellFace[1]);
604: key.j = PetscMax(cellFace[0], cellFace[1]);
605: key.k = PETSC_MAX_INT;
606: key.l = PETSC_MAX_INT;
607: } else {
608: key.i = cellFace[0];
609: key.j = cellFace[1];
610: key.k = cellFace[2];
611: key.l = faceSize > 3 ? (cellFace[3] < 0 ? faceSizeH = 3, PETSC_MAX_INT : cellFace[3]) : PETSC_MAX_INT;
612: PetscSortInt(faceSize, (PetscInt *) &key);
613: }
614: if (faceSizeH != faceSizeAllH) SETERRQ4(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unexpected number of vertices for hybrid face %D of point %D -> %D != %D", cf, c, faceSizeH, faceSizeAllH);
615: PetscHashIJKLPut(faceTable, key, &iter, &missing);
616: if (missing) {
617: DMPlexSetCone(idm, face, cellFace);
618: PetscHashIJKLIterSet(faceTable, iter, face);
619: DMPlexInsertCone(idm, c, cf, face++);
620: } else {
621: const PetscInt *cone;
622: PetscInt coneSize, ornt, i, j, f;
624: PetscHashIJKLIterGet(faceTable, iter, &f);
625: DMPlexInsertCone(idm, c, cf, f);
626: /* Orient face: Do not allow reverse orientation at the first vertex */
627: DMPlexGetConeSize(idm, f, &coneSize);
628: DMPlexGetCone(idm, f, &cone);
629: if (coneSize != faceSizeH) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid number of face vertices %D for face %D should be %D", coneSize, f, faceSizeH);
630: /* - First find the initial vertex */
631: for (i = 0; i < faceSizeH; ++i) if (cellFace[0] == cone[i]) break;
632: /* - Try forward comparison */
633: for (j = 0; j < faceSizeH; ++j) if (cellFace[j] != cone[(i+j)%faceSizeH]) break;
634: if (j == faceSizeH) {
635: if ((faceSizeH == 2) && (i == 1)) ornt = -2;
636: else ornt = i;
637: } else {
638: /* - Try backward comparison */
639: for (j = 0; j < faceSizeH; ++j) if (cellFace[j] != cone[(i+faceSizeH-j)%faceSizeH]) break;
640: if (j == faceSizeH) {
641: if (i == 0) ornt = -faceSizeH;
642: else ornt = -i;
643: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Could not determine face orientation");
644: }
645: DMPlexInsertConeOrientation(idm, c, cf, ornt);
646: }
647: }
648: DMPlexRestoreRawFacesHybrid_Internal(dm, cellDim, coneSize, cone, &numCellFaces, &numCellFacesN, &faceSize, &cellFaces);
649: }
650: if (face != pEnd[faceDepth]) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid number of faces %D should be %D", face-pStart[faceDepth], pEnd[faceDepth]-pStart[faceDepth]);
651: PetscFree2(pStart,pEnd);
652: PetscHashIJKLDestroy(&faceTable);
653: PetscFree2(pStart,pEnd);
654: DMPlexSetHybridBounds(idm, cMax, fMax, eMax, vMax);
655: DMPlexSymmetrize(idm);
656: DMPlexStratify(idm);
657: return(0);
658: }
660: /* This interpolates the PointSF in parallel following local interpolation */
661: static PetscErrorCode DMPlexInterpolatePointSF(DM dm, PetscSF pointSF, PetscInt depth)
662: {
663: PetscMPIInt size, rank;
664: PetscInt p, c, d, dof, offset;
665: PetscInt numLeaves, numRoots, candidatesSize, candidatesRemoteSize;
666: const PetscInt *localPoints;
667: const PetscSFNode *remotePoints;
668: PetscSFNode *candidates, *candidatesRemote, *claims;
669: PetscSection candidateSection, candidateSectionRemote, claimSection;
670: PetscHMapI leafhash;
671: PetscHMapIJ roothash;
672: PetscHashIJKey key;
673: PetscErrorCode ierr;
676: MPI_Comm_size(PetscObjectComm((PetscObject) dm), &size);
677: MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank);
678: PetscSFGetGraph(pointSF, &numRoots, &numLeaves, &localPoints, &remotePoints);
679: if (size < 2 || numRoots < 0) return(0);
680: PetscLogEventBegin(DMPLEX_InterpolateSF,dm,0,0,0);
681: /* Build hashes of points in the SF for efficient lookup */
682: PetscHMapICreate(&leafhash);
683: PetscHMapIJCreate(&roothash);
684: for (p = 0; p < numLeaves; ++p) {
685: PetscHMapISet(leafhash, localPoints[p], p);
686: key.i = remotePoints[p].index;
687: key.j = remotePoints[p].rank;
688: PetscHMapIJSet(roothash, key, p);
689: }
690: /* Build a section / SFNode array of candidate points in the single-level adjacency of leaves,
691: where each candidate is defined by the root entry for the other vertex that defines the edge. */
692: PetscSectionCreate(PetscObjectComm((PetscObject) dm), &candidateSection);
693: PetscSectionSetChart(candidateSection, 0, numRoots);
694: {
695: PetscInt leaf, root, idx, a, *adj = NULL;
696: for (p = 0; p < numLeaves; ++p) {
697: PetscInt adjSize = PETSC_DETERMINE;
698: DMPlexGetAdjacency_Internal(dm, localPoints[p], PETSC_FALSE, PETSC_FALSE, PETSC_FALSE, &adjSize, &adj);
699: for (a = 0; a < adjSize; ++a) {
700: PetscHMapIGet(leafhash, adj[a], &leaf);
701: if (leaf >= 0) {PetscSectionAddDof(candidateSection, localPoints[p], 1);}
702: }
703: }
704: PetscSectionSetUp(candidateSection);
705: PetscSectionGetStorageSize(candidateSection, &candidatesSize);
706: PetscMalloc1(candidatesSize, &candidates);
707: for (p = 0; p < numLeaves; ++p) {
708: PetscInt adjSize = PETSC_DETERMINE;
709: PetscSectionGetOffset(candidateSection, localPoints[p], &offset);
710: DMPlexGetAdjacency_Internal(dm, localPoints[p], PETSC_FALSE, PETSC_FALSE, PETSC_FALSE, &adjSize, &adj);
711: for (idx = 0, a = 0; a < adjSize; ++a) {
712: PetscHMapIGet(leafhash, adj[a], &root);
713: if (root >= 0) candidates[offset+idx++] = remotePoints[root];
714: }
715: }
716: PetscFree(adj);
717: }
718: /* Gather candidate section / array pair into the root partition via inverse(multi(pointSF)). */
719: {
720: PetscSF sfMulti, sfInverse, sfCandidates;
721: PetscInt *remoteOffsets;
722: PetscSFGetMultiSF(pointSF, &sfMulti);
723: PetscSFCreateInverseSF(sfMulti, &sfInverse);
724: PetscSectionCreate(PetscObjectComm((PetscObject) dm), &candidateSectionRemote);
725: PetscSFDistributeSection(sfInverse, candidateSection, &remoteOffsets, candidateSectionRemote);
726: PetscSFCreateSectionSF(sfInverse, candidateSection, remoteOffsets, candidateSectionRemote, &sfCandidates);
727: PetscSectionGetStorageSize(candidateSectionRemote, &candidatesRemoteSize);
728: PetscMalloc1(candidatesRemoteSize, &candidatesRemote);
729: PetscSFBcastBegin(sfCandidates, MPIU_2INT, candidates, candidatesRemote);
730: PetscSFBcastEnd(sfCandidates, MPIU_2INT, candidates, candidatesRemote);
731: PetscSFDestroy(&sfInverse);
732: PetscSFDestroy(&sfCandidates);
733: PetscFree(remoteOffsets);
734: }
735: /* Walk local roots and check for each remote candidate whether we know all required points,
736: either from owning it or having a root entry in the point SF. If we do we place a claim
737: by replacing the vertex number with our edge ID. */
738: {
739: PetscInt idx, root, joinSize, vertices[2];
740: const PetscInt *rootdegree, *join = NULL;
741: PetscSFComputeDegreeBegin(pointSF, &rootdegree);
742: PetscSFComputeDegreeEnd(pointSF, &rootdegree);
743: /* Loop remote edge connections and put in a claim if both vertices are known */
744: for (idx = 0, p = 0; p < numRoots; ++p) {
745: for (d = 0; d < rootdegree[p]; ++d) {
746: PetscSectionGetDof(candidateSectionRemote, idx, &dof);
747: PetscSectionGetOffset(candidateSectionRemote, idx, &offset);
748: for (c = 0; c < dof; ++c) {
749: /* We own both vertices, so we claim the edge by replacing vertex with edge */
750: if (candidatesRemote[offset+c].rank == rank) {
751: vertices[0] = p; vertices[1] = candidatesRemote[offset+c].index;
752: DMPlexGetJoin(dm, 2, vertices, &joinSize, &join);
753: if (joinSize == 1) candidatesRemote[offset+c].index = join[0];
754: DMPlexRestoreJoin(dm, 2, vertices, &joinSize, &join);
755: continue;
756: }
757: /* If we own one vertex and share a root with the other, we claim it */
758: key.i = candidatesRemote[offset+c].index;
759: key.j = candidatesRemote[offset+c].rank;
760: PetscHMapIJGet(roothash, key, &root);
761: if (root >= 0) {
762: vertices[0] = p; vertices[1] = localPoints[root];
763: DMPlexGetJoin(dm, 2, vertices, &joinSize, &join);
764: if (joinSize == 1) {
765: candidatesRemote[offset+c].index = join[0];
766: candidatesRemote[offset+c].rank = rank;
767: }
768: DMPlexRestoreJoin(dm, 2, vertices, &joinSize, &join);
769: }
770: }
771: idx++;
772: }
773: }
774: }
775: /* Push claims back to receiver via the MultiSF and derive new pointSF mapping on receiver */
776: {
777: PetscSF sfMulti, sfClaims, sfPointNew;
778: PetscHMapI claimshash;
779: PetscInt size, pStart, pEnd, root, joinSize, numLocalNew;
780: PetscInt *remoteOffsets, *localPointsNew, vertices[2];
781: const PetscInt *join = NULL;
782: PetscSFNode *remotePointsNew;
783: PetscSFGetMultiSF(pointSF, &sfMulti);
784: PetscSectionCreate(PetscObjectComm((PetscObject) dm), &claimSection);
785: PetscSFDistributeSection(sfMulti, candidateSectionRemote, &remoteOffsets, claimSection);
786: PetscSFCreateSectionSF(sfMulti, candidateSectionRemote, remoteOffsets, claimSection, &sfClaims);
787: PetscSectionGetStorageSize(claimSection, &size);
788: PetscMalloc1(size, &claims);
789: PetscSFBcastBegin(sfClaims, MPIU_2INT, candidatesRemote, claims);
790: PetscSFBcastEnd(sfClaims, MPIU_2INT, candidatesRemote, claims);
791: PetscSFDestroy(&sfClaims);
792: PetscFree(remoteOffsets);
793: /* Walk the original section of local supports and add an SF entry for each updated item */
794: PetscHMapICreate(&claimshash);
795: for (p = 0; p < numRoots; ++p) {
796: PetscSectionGetDof(candidateSection, p, &dof);
797: PetscSectionGetOffset(candidateSection, p, &offset);
798: for (d = 0; d < dof; ++d) {
799: if (candidates[offset+d].index != claims[offset+d].index) {
800: key.i = candidates[offset+d].index;
801: key.j = candidates[offset+d].rank;
802: PetscHMapIJGet(roothash, key, &root);
803: if (root >= 0) {
804: vertices[0] = p; vertices[1] = localPoints[root];
805: DMPlexGetJoin(dm, 2, vertices, &joinSize, &join);
806: if (joinSize == 1) {PetscHMapISet(claimshash, join[0], offset+d);}
807: DMPlexRestoreJoin(dm, 2, vertices, &joinSize, &join);
808: }
809: }
810: }
811: }
812: /* Create new pointSF from hashed claims */
813: PetscHMapIGetSize(claimshash, &numLocalNew);
814: DMPlexGetChart(dm, &pStart, &pEnd);
815: PetscMalloc1(numLeaves + numLocalNew, &localPointsNew);
816: PetscMalloc1(numLeaves + numLocalNew, &remotePointsNew);
817: for (p = 0; p < numLeaves; ++p) {
818: localPointsNew[p] = localPoints[p];
819: remotePointsNew[p].index = remotePoints[p].index;
820: remotePointsNew[p].rank = remotePoints[p].rank;
821: }
822: p = numLeaves;
823: PetscHMapIGetKeys(claimshash, &p, localPointsNew);
824: PetscSortInt(numLocalNew, &localPointsNew[numLeaves]);
825: for (p = numLeaves; p < numLeaves + numLocalNew; ++p) {
826: PetscHMapIGet(claimshash, localPointsNew[p], &offset);
827: remotePointsNew[p] = claims[offset];
828: }
829: PetscSFCreate(PetscObjectComm((PetscObject) dm), &sfPointNew);
830: PetscSFSetGraph(sfPointNew, pEnd-pStart, numLeaves+numLocalNew, localPointsNew, PETSC_OWN_POINTER, remotePointsNew, PETSC_OWN_POINTER);
831: DMSetPointSF(dm, sfPointNew);
832: PetscSFDestroy(&sfPointNew);
833: PetscHMapIDestroy(&claimshash);
834: }
835: PetscHMapIDestroy(&leafhash);
836: PetscHMapIJDestroy(&roothash);
837: PetscSectionDestroy(&candidateSection);
838: PetscSectionDestroy(&candidateSectionRemote);
839: PetscSectionDestroy(&claimSection);
840: PetscFree(candidates);
841: PetscFree(candidatesRemote);
842: PetscFree(claims);
843: PetscLogEventEnd(DMPLEX_InterpolateSF,dm,0,0,0);
844: return(0);
845: }
847: /*@C
848: DMPlexInterpolate - Take in a cell-vertex mesh and return one with all intermediate faces, edges, etc.
850: Collective on DM
852: Input Parameters:
853: + dm - The DMPlex object with only cells and vertices
854: - dmInt - The interpolated DM
856: Output Parameter:
857: . dmInt - The complete DMPlex object
859: Level: intermediate
861: Notes:
862: It does not copy over the coordinates.
864: .keywords: mesh
865: .seealso: DMPlexUninterpolate(), DMPlexCreateFromCellList(), DMPlexCopyCoordinates()
866: @*/
867: PetscErrorCode DMPlexInterpolate(DM dm, DM *dmInt)
868: {
869: DM idm, odm = dm;
870: PetscSF sfPoint;
871: PetscInt depth, dim, d;
872: const char *name;
878: PetscLogEventBegin(DMPLEX_Interpolate,dm,0,0,0);
879: DMPlexGetDepth(dm, &depth);
880: DMGetDimension(dm, &dim);
881: if ((depth == dim) || (dim <= 1)) {
882: PetscObjectReference((PetscObject) dm);
883: idm = dm;
884: } else {
885: for (d = 1; d < dim; ++d) {
886: /* Create interpolated mesh */
887: DMCreate(PetscObjectComm((PetscObject)dm), &idm);
888: DMSetType(idm, DMPLEX);
889: DMSetDimension(idm, dim);
890: if (depth > 0) {
891: DMPlexInterpolateFaces_Internal(odm, 1, idm);
892: DMGetPointSF(odm, &sfPoint);
893: DMPlexInterpolatePointSF(idm, sfPoint, depth);
894: }
895: if (odm != dm) {DMDestroy(&odm);}
896: odm = idm;
897: }
898: PetscObjectGetName((PetscObject) dm, &name);
899: PetscObjectSetName((PetscObject) idm, name);
900: DMPlexCopyCoordinates(dm, idm);
901: DMCopyLabels(dm, idm);
902: }
903: {
904: PetscBool isper;
905: const PetscReal *maxCell, *L;
906: const DMBoundaryType *bd;
908: DMGetPeriodicity(dm,&isper,&maxCell,&L,&bd);
909: DMSetPeriodicity(idm,isper,maxCell,L,bd);
910: }
911: *dmInt = idm;
912: PetscLogEventEnd(DMPLEX_Interpolate,dm,0,0,0);
913: return(0);
914: }
916: /*@
917: DMPlexCopyCoordinates - Copy coordinates from one mesh to another with the same vertices
919: Collective on DM
921: Input Parameter:
922: . dmA - The DMPlex object with initial coordinates
924: Output Parameter:
925: . dmB - The DMPlex object with copied coordinates
927: Level: intermediate
929: Note: This is typically used when adding pieces other than vertices to a mesh
931: .keywords: mesh
932: .seealso: DMCopyLabels(), DMGetCoordinates(), DMGetCoordinatesLocal(), DMGetCoordinateDM(), DMGetCoordinateSection()
933: @*/
934: PetscErrorCode DMPlexCopyCoordinates(DM dmA, DM dmB)
935: {
936: Vec coordinatesA, coordinatesB;
937: VecType vtype;
938: PetscSection coordSectionA, coordSectionB;
939: PetscScalar *coordsA, *coordsB;
940: PetscInt spaceDim, Nf, vStartA, vStartB, vEndA, vEndB, coordSizeB, v, d;
941: PetscInt cStartA, cEndA, cStartB, cEndB, cS, cE;
942: PetscBool lc = PETSC_FALSE;
948: if (dmA == dmB) return(0);
949: DMPlexGetDepthStratum(dmA, 0, &vStartA, &vEndA);
950: DMPlexGetDepthStratum(dmB, 0, &vStartB, &vEndB);
951: 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);
952: DMPlexGetHeightStratum(dmA, 0, &cStartA, &cEndA);
953: DMPlexGetHeightStratum(dmB, 0, &cStartB, &cEndB);
954: DMGetCoordinateSection(dmA, &coordSectionA);
955: DMGetCoordinateSection(dmB, &coordSectionB);
956: if (coordSectionA == coordSectionB) return(0);
957: PetscSectionGetNumFields(coordSectionA, &Nf);
958: if (!Nf) return(0);
959: if (Nf > 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "The number of coordinate fields must be 1, not %D", Nf);
960: if (!coordSectionB) {
961: PetscInt dim;
963: PetscSectionCreate(PetscObjectComm((PetscObject) coordSectionA), &coordSectionB);
964: DMGetCoordinateDim(dmA, &dim);
965: DMSetCoordinateSection(dmB, dim, coordSectionB);
966: PetscObjectDereference((PetscObject) coordSectionB);
967: }
968: PetscSectionSetNumFields(coordSectionB, 1);
969: PetscSectionGetFieldComponents(coordSectionA, 0, &spaceDim);
970: PetscSectionSetFieldComponents(coordSectionB, 0, spaceDim);
971: PetscSectionGetChart(coordSectionA, &cS, &cE);
972: if (cStartA <= cS && cS < cEndA) { /* localized coordinates */
973: if ((cEndA-cStartA) != (cEndB-cStartB)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "The number of cellls in first DM %d != %d in the second DM", cEndA-cStartA, cEndB-cStartB);
974: cS = cS - cStartA + cStartB;
975: cE = vEndB;
976: lc = PETSC_TRUE;
977: } else {
978: cS = vStartB;
979: cE = vEndB;
980: }
981: PetscSectionSetChart(coordSectionB, cS, cE);
982: for (v = vStartB; v < vEndB; ++v) {
983: PetscSectionSetDof(coordSectionB, v, spaceDim);
984: PetscSectionSetFieldDof(coordSectionB, v, 0, spaceDim);
985: }
986: if (lc) { /* localized coordinates */
987: PetscInt c;
989: for (c = cS-cStartB; c < cEndB-cStartB; c++) {
990: PetscInt dof;
992: PetscSectionGetDof(coordSectionA, c + cStartA, &dof);
993: PetscSectionSetDof(coordSectionB, c + cStartB, dof);
994: PetscSectionSetFieldDof(coordSectionB, c + cStartB, 0, dof);
995: }
996: }
997: PetscSectionSetUp(coordSectionB);
998: PetscSectionGetStorageSize(coordSectionB, &coordSizeB);
999: DMGetCoordinatesLocal(dmA, &coordinatesA);
1000: VecCreate(PETSC_COMM_SELF, &coordinatesB);
1001: PetscObjectSetName((PetscObject) coordinatesB, "coordinates");
1002: VecSetSizes(coordinatesB, coordSizeB, PETSC_DETERMINE);
1003: VecGetBlockSize(coordinatesA, &d);
1004: VecSetBlockSize(coordinatesB, d);
1005: VecGetType(coordinatesA, &vtype);
1006: VecSetType(coordinatesB, vtype);
1007: VecGetArray(coordinatesA, &coordsA);
1008: VecGetArray(coordinatesB, &coordsB);
1009: for (v = 0; v < vEndB-vStartB; ++v) {
1010: PetscInt offA, offB;
1012: PetscSectionGetOffset(coordSectionA, v + vStartA, &offA);
1013: PetscSectionGetOffset(coordSectionB, v + vStartB, &offB);
1014: for (d = 0; d < spaceDim; ++d) {
1015: coordsB[offB+d] = coordsA[offA+d];
1016: }
1017: }
1018: if (lc) { /* localized coordinates */
1019: PetscInt c;
1021: for (c = cS-cStartB; c < cEndB-cStartB; c++) {
1022: PetscInt dof, offA, offB;
1024: PetscSectionGetOffset(coordSectionA, c + cStartA, &offA);
1025: PetscSectionGetOffset(coordSectionB, c + cStartB, &offB);
1026: PetscSectionGetDof(coordSectionA, c + cStartA, &dof);
1027: PetscMemcpy(coordsB + offB,coordsA + offA,dof*sizeof(*coordsB));
1028: }
1029: }
1030: VecRestoreArray(coordinatesA, &coordsA);
1031: VecRestoreArray(coordinatesB, &coordsB);
1032: DMSetCoordinatesLocal(dmB, coordinatesB);
1033: VecDestroy(&coordinatesB);
1034: return(0);
1035: }
1037: /*@
1038: DMPlexUninterpolate - Take in a mesh with all intermediate faces, edges, etc. and return a cell-vertex mesh
1040: Collective on DM
1042: Input Parameter:
1043: . dm - The complete DMPlex object
1045: Output Parameter:
1046: . dmUnint - The DMPlex object with only cells and vertices
1048: Level: intermediate
1050: Notes:
1051: It does not copy over the coordinates.
1053: .keywords: mesh
1054: .seealso: DMPlexInterpolate(), DMPlexCreateFromCellList(), DMPlexCopyCoordinates()
1055: @*/
1056: PetscErrorCode DMPlexUninterpolate(DM dm, DM *dmUnint)
1057: {
1058: DM udm;
1059: PetscInt dim, vStart, vEnd, cStart, cEnd, cMax, c, maxConeSize = 0, *cone;
1065: DMGetDimension(dm, &dim);
1066: if (dim <= 1) {
1067: PetscObjectReference((PetscObject) dm);
1068: *dmUnint = dm;
1069: return(0);
1070: }
1071: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
1072: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
1073: DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);
1074: DMCreate(PetscObjectComm((PetscObject) dm), &udm);
1075: DMSetType(udm, DMPLEX);
1076: DMSetDimension(udm, dim);
1077: DMPlexSetChart(udm, cStart, vEnd);
1078: for (c = cStart; c < cEnd; ++c) {
1079: PetscInt *closure = NULL, closureSize, cl, coneSize = 0;
1081: DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
1082: for (cl = 0; cl < closureSize*2; cl += 2) {
1083: const PetscInt p = closure[cl];
1085: if ((p >= vStart) && (p < vEnd)) ++coneSize;
1086: }
1087: DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
1088: DMPlexSetConeSize(udm, c, coneSize);
1089: maxConeSize = PetscMax(maxConeSize, coneSize);
1090: }
1091: DMSetUp(udm);
1092: PetscMalloc1(maxConeSize, &cone);
1093: for (c = cStart; c < cEnd; ++c) {
1094: PetscInt *closure = NULL, closureSize, cl, coneSize = 0;
1096: DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
1097: for (cl = 0; cl < closureSize*2; cl += 2) {
1098: const PetscInt p = closure[cl];
1100: if ((p >= vStart) && (p < vEnd)) cone[coneSize++] = p;
1101: }
1102: DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
1103: DMPlexSetCone(udm, c, cone);
1104: }
1105: PetscFree(cone);
1106: DMPlexSetHybridBounds(udm, cMax, PETSC_DETERMINE, PETSC_DETERMINE, PETSC_DETERMINE);
1107: DMPlexSymmetrize(udm);
1108: DMPlexStratify(udm);
1109: /* Reduce SF */
1110: {
1111: PetscSF sfPoint, sfPointUn;
1112: const PetscSFNode *remotePoints;
1113: const PetscInt *localPoints;
1114: PetscSFNode *remotePointsUn;
1115: PetscInt *localPointsUn;
1116: PetscInt vEnd, numRoots, numLeaves, l;
1117: PetscInt numLeavesUn = 0, n = 0;
1118: PetscErrorCode ierr;
1120: /* Get original SF information */
1121: DMGetPointSF(dm, &sfPoint);
1122: DMGetPointSF(udm, &sfPointUn);
1123: DMPlexGetDepthStratum(dm, 0, NULL, &vEnd);
1124: PetscSFGetGraph(sfPoint, &numRoots, &numLeaves, &localPoints, &remotePoints);
1125: /* Allocate space for cells and vertices */
1126: for (l = 0; l < numLeaves; ++l) if (localPoints[l] < vEnd) numLeavesUn++;
1127: /* Fill in leaves */
1128: if (vEnd >= 0) {
1129: PetscMalloc1(numLeavesUn, &remotePointsUn);
1130: PetscMalloc1(numLeavesUn, &localPointsUn);
1131: for (l = 0; l < numLeaves; l++) {
1132: if (localPoints[l] < vEnd) {
1133: localPointsUn[n] = localPoints[l];
1134: remotePointsUn[n].rank = remotePoints[l].rank;
1135: remotePointsUn[n].index = remotePoints[l].index;
1136: ++n;
1137: }
1138: }
1139: if (n != numLeavesUn) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Inconsistent number of leaves %d != %d", n, numLeavesUn);
1140: PetscSFSetGraph(sfPointUn, vEnd, numLeavesUn, localPointsUn, PETSC_OWN_POINTER, remotePointsUn, PETSC_OWN_POINTER);
1141: }
1142: }
1143: {
1144: PetscBool isper;
1145: const PetscReal *maxCell, *L;
1146: const DMBoundaryType *bd;
1148: DMGetPeriodicity(dm,&isper,&maxCell,&L,&bd);
1149: DMSetPeriodicity(udm,isper,maxCell,L,bd);
1150: }
1152: *dmUnint = udm;
1153: return(0);
1154: }