Actual source code: plexinterpolate.c
petsc-3.8.4 2018-03-24
1: #include <petsc/private/dmpleximpl.h>
3: /*
4: DMPlexGetFaces_Internal - Gets groups of vertices that correspond to faces for the given cell
5: */
6: PetscErrorCode DMPlexGetFaces_Internal(DM dm, PetscInt dim, PetscInt p, PetscInt *numFaces, PetscInt *faceSize, const PetscInt *faces[])
7: {
8: const PetscInt *cone = NULL;
9: PetscInt maxConeSize, maxSupportSize, coneSize;
10: PetscErrorCode ierr;
14: DMPlexGetMaxSizes(dm, &maxConeSize, &maxSupportSize);
15: DMPlexGetConeSize(dm, p, &coneSize);
16: DMPlexGetCone(dm, p, &cone);
17: DMPlexGetRawFaces_Internal(dm, dim, coneSize, cone, numFaces, faceSize, faces);
18: return(0);
19: }
21: /*
22: DMPlexRestoreFaces_Internal - Restores the array
23: */
24: PetscErrorCode DMPlexRestoreFaces_Internal(DM dm, PetscInt dim, PetscInt p, PetscInt *numFaces, PetscInt *faceSize, const PetscInt *faces[])
25: {
26: PetscErrorCode ierr;
29: DMRestoreWorkArray(dm, 0, PETSC_INT, (void *) faces);
30: return(0);
31: }
33: /*
34: DMPlexGetRawFaces_Internal - Gets groups of vertices that correspond to faces for the given cone
35: */
36: PetscErrorCode DMPlexGetRawFaces_Internal(DM dm, PetscInt dim, PetscInt coneSize, const PetscInt cone[], PetscInt *numFaces, PetscInt *faceSize, const PetscInt *faces[])
37: {
38: PetscInt *facesTmp;
39: PetscInt maxConeSize, maxSupportSize;
40: PetscErrorCode ierr;
44: DMPlexGetMaxSizes(dm, &maxConeSize, &maxSupportSize);
45: if (faces) {DMGetWorkArray(dm, PetscSqr(PetscMax(maxConeSize, maxSupportSize)), PETSC_INT, &facesTmp);}
46: switch (dim) {
47: case 1:
48: switch (coneSize) {
49: case 2:
50: if (faces) {
51: facesTmp[0] = cone[0]; facesTmp[1] = cone[1];
52: *faces = facesTmp;
53: }
54: if (numFaces) *numFaces = 2;
55: if (faceSize) *faceSize = 1;
56: break;
57: default:
58: SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Cone size %D not supported for dimension %D", coneSize, dim);
59: }
60: break;
61: case 2:
62: switch (coneSize) {
63: case 3:
64: if (faces) {
65: facesTmp[0] = cone[0]; facesTmp[1] = cone[1];
66: facesTmp[2] = cone[1]; facesTmp[3] = cone[2];
67: facesTmp[4] = cone[2]; facesTmp[5] = cone[0];
68: *faces = facesTmp;
69: }
70: if (numFaces) *numFaces = 3;
71: if (faceSize) *faceSize = 2;
72: break;
73: case 4:
74: /* Vertices follow right hand rule */
75: if (faces) {
76: facesTmp[0] = cone[0]; facesTmp[1] = cone[1];
77: facesTmp[2] = cone[1]; facesTmp[3] = cone[2];
78: facesTmp[4] = cone[2]; facesTmp[5] = cone[3];
79: facesTmp[6] = cone[3]; facesTmp[7] = cone[0];
80: *faces = facesTmp;
81: }
82: if (numFaces) *numFaces = 4;
83: if (faceSize) *faceSize = 2;
84: break;
85: default:
86: SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Cone size %D not supported for dimension %D", coneSize, dim);
87: }
88: break;
89: case 3:
90: switch (coneSize) {
91: case 3:
92: if (faces) {
93: facesTmp[0] = cone[0]; facesTmp[1] = cone[1];
94: facesTmp[2] = cone[1]; facesTmp[3] = cone[2];
95: facesTmp[4] = cone[2]; facesTmp[5] = cone[0];
96: *faces = facesTmp;
97: }
98: if (numFaces) *numFaces = 3;
99: if (faceSize) *faceSize = 2;
100: break;
101: case 4:
102: /* Vertices of first face follow right hand rule and normal points away from last vertex */
103: if (faces) {
104: facesTmp[0] = cone[0]; facesTmp[1] = cone[1]; facesTmp[2] = cone[2];
105: facesTmp[3] = cone[0]; facesTmp[4] = cone[3]; facesTmp[5] = cone[1];
106: facesTmp[6] = cone[0]; facesTmp[7] = cone[2]; facesTmp[8] = cone[3];
107: facesTmp[9] = cone[2]; facesTmp[10] = cone[1]; facesTmp[11] = cone[3];
108: *faces = facesTmp;
109: }
110: if (numFaces) *numFaces = 4;
111: if (faceSize) *faceSize = 3;
112: break;
113: case 8:
114: if (faces) {
115: facesTmp[0] = cone[0]; facesTmp[1] = cone[1]; facesTmp[2] = cone[2]; facesTmp[3] = cone[3]; /* Bottom */
116: facesTmp[4] = cone[4]; facesTmp[5] = cone[5]; facesTmp[6] = cone[6]; facesTmp[7] = cone[7]; /* Top */
117: facesTmp[8] = cone[0]; facesTmp[9] = cone[3]; facesTmp[10] = cone[5]; facesTmp[11] = cone[4]; /* Front */
118: facesTmp[12] = cone[2]; facesTmp[13] = cone[1]; facesTmp[14] = cone[7]; facesTmp[15] = cone[6]; /* Back */
119: facesTmp[16] = cone[3]; facesTmp[17] = cone[2]; facesTmp[18] = cone[6]; facesTmp[19] = cone[5]; /* Right */
120: facesTmp[20] = cone[0]; facesTmp[21] = cone[4]; facesTmp[22] = cone[7]; facesTmp[23] = cone[1]; /* Left */
121: *faces = facesTmp;
122: }
123: if (numFaces) *numFaces = 6;
124: if (faceSize) *faceSize = 4;
125: break;
126: default:
127: SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Cone size %D not supported for dimension %D", coneSize, dim);
128: }
129: break;
130: default:
131: SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Dimension %D not supported", dim);
132: }
133: return(0);
134: }
136: /* This interpolates faces for cells at some stratum */
137: static PetscErrorCode DMPlexInterpolateFaces_Internal(DM dm, PetscInt cellDepth, DM idm)
138: {
139: DMLabel subpointMap;
140: PetscHashIJKL faceTable;
141: PetscInt *pStart, *pEnd;
142: PetscInt cellDim, depth, faceDepth = cellDepth, numPoints = 0, faceSizeAll = 0, face, c, d;
146: DMGetDimension(dm, &cellDim);
147: /* HACK: I need a better way to determine face dimension, or an alternative to GetFaces() */
148: DMPlexGetSubpointMap(dm, &subpointMap);
149: if (subpointMap) ++cellDim;
150: DMPlexGetDepth(dm, &depth);
151: ++depth;
152: ++cellDepth;
153: cellDim -= depth - cellDepth;
154: PetscMalloc2(depth+1,&pStart,depth+1,&pEnd);
155: for (d = depth-1; d >= faceDepth; --d) {
156: DMPlexGetDepthStratum(dm, d, &pStart[d+1], &pEnd[d+1]);
157: }
158: DMPlexGetDepthStratum(dm, -1, NULL, &pStart[faceDepth]);
159: pEnd[faceDepth] = pStart[faceDepth];
160: for (d = faceDepth-1; d >= 0; --d) {
161: DMPlexGetDepthStratum(dm, d, &pStart[d], &pEnd[d]);
162: }
163: if (pEnd[cellDepth] > pStart[cellDepth]) {DMPlexGetFaces_Internal(dm, cellDim, pStart[cellDepth], NULL, &faceSizeAll, NULL);}
164: if (faceSizeAll > 4) SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "Do not support interpolation of meshes with faces of %D vertices", faceSizeAll);
165: PetscHashIJKLCreate(&faceTable);
166: for (c = pStart[cellDepth], face = pStart[faceDepth]; c < pEnd[cellDepth]; ++c) {
167: const PetscInt *cellFaces;
168: PetscInt numCellFaces, faceSize, cf;
170: DMPlexGetFaces_Internal(dm, cellDim, c, &numCellFaces, &faceSize, &cellFaces);
171: if (faceSize != faceSizeAll) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Inconsistent face for cell %D of size %D != %D", c, faceSize, faceSizeAll);
172: for (cf = 0; cf < numCellFaces; ++cf) {
173: const PetscInt *cellFace = &cellFaces[cf*faceSize];
174: PetscHashIJKLKey key;
175: PetscHashIJKLIter missing, iter;
177: if (faceSize == 2) {
178: key.i = PetscMin(cellFace[0], cellFace[1]);
179: key.j = PetscMax(cellFace[0], cellFace[1]);
180: key.k = 0;
181: key.l = 0;
182: } else {
183: key.i = cellFace[0]; key.j = cellFace[1]; key.k = cellFace[2]; key.l = faceSize > 3 ? cellFace[3] : 0;
184: PetscSortInt(faceSize, (PetscInt *) &key);
185: }
186: PetscHashIJKLPut(faceTable, key, &missing, &iter);
187: if (missing) {PetscHashIJKLSet(faceTable, iter, face++);}
188: }
189: DMPlexRestoreFaces_Internal(dm, cellDim, c, &numCellFaces, &faceSize, &cellFaces);
190: }
191: pEnd[faceDepth] = face;
192: PetscHashIJKLDestroy(&faceTable);
193: /* Count new points */
194: for (d = 0; d <= depth; ++d) {
195: numPoints += pEnd[d]-pStart[d];
196: }
197: DMPlexSetChart(idm, 0, numPoints);
198: /* Set cone sizes */
199: for (d = 0; d <= depth; ++d) {
200: PetscInt coneSize, p;
202: if (d == faceDepth) {
203: for (p = pStart[d]; p < pEnd[d]; ++p) {
204: /* I see no way to do this if we admit faces of different shapes */
205: DMPlexSetConeSize(idm, p, faceSizeAll);
206: }
207: } else if (d == cellDepth) {
208: for (p = pStart[d]; p < pEnd[d]; ++p) {
209: /* Number of cell faces may be different from number of cell vertices*/
210: DMPlexGetFaces_Internal(dm, cellDim, p, &coneSize, NULL, NULL);
211: DMPlexSetConeSize(idm, p, coneSize);
212: }
213: } else {
214: for (p = pStart[d]; p < pEnd[d]; ++p) {
215: DMPlexGetConeSize(dm, p, &coneSize);
216: DMPlexSetConeSize(idm, p, coneSize);
217: }
218: }
219: }
220: DMSetUp(idm);
221: /* Get face cones from subsets of cell vertices */
222: if (faceSizeAll > 4) SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "Do not support interpolation of meshes with faces of %D vertices", faceSizeAll);
223: PetscHashIJKLCreate(&faceTable);
224: for (d = depth; d > cellDepth; --d) {
225: const PetscInt *cone;
226: PetscInt p;
228: for (p = pStart[d]; p < pEnd[d]; ++p) {
229: DMPlexGetCone(dm, p, &cone);
230: DMPlexSetCone(idm, p, cone);
231: DMPlexGetConeOrientation(dm, p, &cone);
232: DMPlexSetConeOrientation(idm, p, cone);
233: }
234: }
235: for (c = pStart[cellDepth], face = pStart[faceDepth]; c < pEnd[cellDepth]; ++c) {
236: const PetscInt *cellFaces;
237: PetscInt numCellFaces, faceSize, cf;
239: DMPlexGetFaces_Internal(dm, cellDim, c, &numCellFaces, &faceSize, &cellFaces);
240: if (faceSize != faceSizeAll) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Inconsistent face for cell %D of size %D != %D", c, faceSize, faceSizeAll);
241: for (cf = 0; cf < numCellFaces; ++cf) {
242: const PetscInt *cellFace = &cellFaces[cf*faceSize];
243: PetscHashIJKLKey key;
244: PetscHashIJKLIter missing, iter;
246: if (faceSize == 2) {
247: key.i = PetscMin(cellFace[0], cellFace[1]);
248: key.j = PetscMax(cellFace[0], cellFace[1]);
249: key.k = 0;
250: key.l = 0;
251: } else {
252: key.i = cellFace[0]; key.j = cellFace[1]; key.k = cellFace[2]; key.l = faceSize > 3 ? cellFace[3] : 0;
253: PetscSortInt(faceSize, (PetscInt *) &key);
254: }
255: PetscHashIJKLPut(faceTable, key, &missing, &iter);
256: if (missing) {
257: DMPlexSetCone(idm, face, cellFace);
258: PetscHashIJKLSet(faceTable, iter, face);
259: DMPlexInsertCone(idm, c, cf, face++);
260: } else {
261: const PetscInt *cone;
262: PetscInt coneSize, ornt, i, j, f;
264: PetscHashIJKLGet(faceTable, iter, &f);
265: DMPlexInsertCone(idm, c, cf, f);
266: /* Orient face: Do not allow reverse orientation at the first vertex */
267: DMPlexGetConeSize(idm, f, &coneSize);
268: DMPlexGetCone(idm, f, &cone);
269: 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);
270: /* - First find the initial vertex */
271: for (i = 0; i < faceSize; ++i) if (cellFace[0] == cone[i]) break;
272: /* - Try forward comparison */
273: for (j = 0; j < faceSize; ++j) if (cellFace[j] != cone[(i+j)%faceSize]) break;
274: if (j == faceSize) {
275: if ((faceSize == 2) && (i == 1)) ornt = -2;
276: else ornt = i;
277: } else {
278: /* - Try backward comparison */
279: for (j = 0; j < faceSize; ++j) if (cellFace[j] != cone[(i+faceSize-j)%faceSize]) break;
280: if (j == faceSize) {
281: if (i == 0) ornt = -faceSize;
282: else ornt = -i;
283: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Could not determine face orientation");
284: }
285: DMPlexInsertConeOrientation(idm, c, cf, ornt);
286: }
287: }
288: DMPlexRestoreFaces_Internal(dm, cellDim, c, &numCellFaces, &faceSize, &cellFaces);
289: }
290: if (face != pEnd[faceDepth]) SETERRQ2(PetscObjectComm((PetscObject) dm), PETSC_ERR_PLIB, "Invalid number of faces %D should be %D", face-pStart[faceDepth], pEnd[faceDepth]-pStart[faceDepth]);
291: PetscFree2(pStart,pEnd);
292: PetscHashIJKLDestroy(&faceTable);
293: PetscFree2(pStart,pEnd);
294: DMPlexSymmetrize(idm);
295: DMPlexStratify(idm);
296: return(0);
297: }
299: /* This interpolates the PointSF in parallel following local interpolation */
300: static PetscErrorCode DMPlexInterpolatePointSF(DM dm, PetscSF pointSF, PetscInt depth)
301: {
302: PetscMPIInt size, rank;
303: PetscInt p, c, d, dof, offset;
304: PetscInt numLeaves, numRoots, candidatesSize, candidatesRemoteSize;
305: const PetscInt *localPoints;
306: const PetscSFNode *remotePoints;
307: PetscSFNode *candidates, *candidatesRemote, *claims;
308: PetscSection candidateSection, candidateSectionRemote, claimSection;
309: PetscHashI leafhash;
310: PetscHashIJ roothash;
311: PetscHashIJKey key;
312: PetscErrorCode ierr;
315: MPI_Comm_size(PetscObjectComm((PetscObject) dm), &size);
316: MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank);
317: PetscSFGetGraph(pointSF, &numRoots, &numLeaves, &localPoints, &remotePoints);
318: if (size < 2 || numRoots < 0) return(0);
319: PetscLogEventBegin(DMPLEX_InterpolateSF,dm,0,0,0);
320: /* Build hashes of points in the SF for efficient lookup */
321: PetscHashICreate(leafhash);
322: PetscHashIJCreate(&roothash);
323: PetscHashIJSetMultivalued(roothash, PETSC_FALSE);
324: for (p = 0; p < numLeaves; ++p) {
325: PetscHashIAdd(leafhash, localPoints[p], p);
326: key.i = remotePoints[p].index; key.j = remotePoints[p].rank;
327: PetscHashIJAdd(roothash, key, p);
328: }
329: /* Build a section / SFNode array of candidate points in the single-level adjacency of leaves,
330: where each candidate is defined by the root entry for the other vertex that defines the edge. */
331: PetscSectionCreate(PetscObjectComm((PetscObject) dm), &candidateSection);
332: PetscSectionSetChart(candidateSection, 0, numRoots);
333: {
334: PetscInt leaf, root, idx, a, *adj = NULL;
335: for (p = 0; p < numLeaves; ++p) {
336: PetscInt adjSize = PETSC_DETERMINE;
337: DMPlexGetAdjacency_Internal(dm, localPoints[p], PETSC_FALSE, PETSC_FALSE, PETSC_FALSE, &adjSize, &adj);
338: for (a = 0; a < adjSize; ++a) {
339: PetscHashIMap(leafhash, adj[a], leaf);
340: if (leaf >= 0) {PetscSectionAddDof(candidateSection, localPoints[p], 1);}
341: }
342: }
343: PetscSectionSetUp(candidateSection);
344: PetscSectionGetStorageSize(candidateSection, &candidatesSize);
345: PetscMalloc1(candidatesSize, &candidates);
346: for (p = 0; p < numLeaves; ++p) {
347: PetscInt adjSize = PETSC_DETERMINE;
348: PetscSectionGetOffset(candidateSection, localPoints[p], &offset);
349: DMPlexGetAdjacency_Internal(dm, localPoints[p], PETSC_FALSE, PETSC_FALSE, PETSC_FALSE, &adjSize, &adj);
350: for (idx = 0, a = 0; a < adjSize; ++a) {
351: PetscHashIMap(leafhash, adj[a], root);
352: if (root >= 0) candidates[offset+idx++] = remotePoints[root];
353: }
354: }
355: PetscFree(adj);
356: }
357: /* Gather candidate section / array pair into the root partition via inverse(multi(pointSF)). */
358: {
359: PetscSF sfMulti, sfInverse, sfCandidates;
360: PetscInt *remoteOffsets;
361: PetscSFGetMultiSF(pointSF, &sfMulti);
362: PetscSFCreateInverseSF(sfMulti, &sfInverse);
363: PetscSectionCreate(PetscObjectComm((PetscObject) dm), &candidateSectionRemote);
364: PetscSFDistributeSection(sfInverse, candidateSection, &remoteOffsets, candidateSectionRemote);
365: PetscSFCreateSectionSF(sfInverse, candidateSection, remoteOffsets, candidateSectionRemote, &sfCandidates);
366: PetscSectionGetStorageSize(candidateSectionRemote, &candidatesRemoteSize);
367: PetscMalloc1(candidatesRemoteSize, &candidatesRemote);
368: PetscSFBcastBegin(sfCandidates, MPIU_2INT, candidates, candidatesRemote);
369: PetscSFBcastEnd(sfCandidates, MPIU_2INT, candidates, candidatesRemote);
370: PetscSFDestroy(&sfInverse);
371: PetscSFDestroy(&sfCandidates);
372: PetscFree(remoteOffsets);
373: }
374: /* Walk local roots and check for each remote candidate whether we know all required points,
375: either from owning it or having a root entry in the point SF. If we do we place a claim
376: by replacing the vertex number with our edge ID. */
377: {
378: PetscInt idx, root, joinSize, vertices[2];
379: const PetscInt *rootdegree, *join = NULL;
380: PetscSFComputeDegreeBegin(pointSF, &rootdegree);
381: PetscSFComputeDegreeEnd(pointSF, &rootdegree);
382: /* Loop remote edge connections and put in a claim if both vertices are known */
383: for (idx = 0, p = 0; p < numRoots; ++p) {
384: for (d = 0; d < rootdegree[p]; ++d) {
385: PetscSectionGetDof(candidateSectionRemote, idx, &dof);
386: PetscSectionGetOffset(candidateSectionRemote, idx, &offset);
387: for (c = 0; c < dof; ++c) {
388: /* We own both vertices, so we claim the edge by replacing vertex with edge */
389: if (candidatesRemote[offset+c].rank == rank) {
390: vertices[0] = p; vertices[1] = candidatesRemote[offset+c].index;
391: DMPlexGetJoin(dm, 2, vertices, &joinSize, &join);
392: if (joinSize == 1) candidatesRemote[offset+c].index = join[0];
393: DMPlexRestoreJoin(dm, 2, vertices, &joinSize, &join);
394: continue;
395: }
396: /* If we own one vertex and share a root with the other, we claim it */
397: key.i = candidatesRemote[offset+c].index; key.j = candidatesRemote[offset+c].rank;
398: PetscHashIJGet(roothash, key, &root);
399: if (root >= 0) {
400: vertices[0] = p; vertices[1] = localPoints[root];
401: DMPlexGetJoin(dm, 2, vertices, &joinSize, &join);
402: if (joinSize == 1) {
403: candidatesRemote[offset+c].index = join[0];
404: candidatesRemote[offset+c].rank = rank;
405: }
406: DMPlexRestoreJoin(dm, 2, vertices, &joinSize, &join);
407: }
408: }
409: idx++;
410: }
411: }
412: }
413: /* Push claims back to receiver via the MultiSF and derive new pointSF mapping on receiver */
414: {
415: PetscSF sfMulti, sfClaims, sfPointNew;
416: PetscHashI claimshash;
417: PetscInt size, pStart, pEnd, root, joinSize, numLocalNew;
418: PetscInt *remoteOffsets, *localPointsNew, vertices[2];
419: const PetscInt *join = NULL;
420: PetscSFNode *remotePointsNew;
421: PetscSFGetMultiSF(pointSF, &sfMulti);
422: PetscSectionCreate(PetscObjectComm((PetscObject) dm), &claimSection);
423: PetscSFDistributeSection(sfMulti, candidateSectionRemote, &remoteOffsets, claimSection);
424: PetscSFCreateSectionSF(sfMulti, candidateSectionRemote, remoteOffsets, claimSection, &sfClaims);
425: PetscSectionGetStorageSize(claimSection, &size);
426: PetscMalloc1(size, &claims);
427: PetscSFBcastBegin(sfClaims, MPIU_2INT, candidatesRemote, claims);
428: PetscSFBcastEnd(sfClaims, MPIU_2INT, candidatesRemote, claims);
429: PetscSFDestroy(&sfClaims);
430: PetscFree(remoteOffsets);
431: /* Walk the original section of local supports and add an SF entry for each updated item */
432: PetscHashICreate(claimshash);
433: for (p = 0; p < numRoots; ++p) {
434: PetscSectionGetDof(candidateSection, p, &dof);
435: PetscSectionGetOffset(candidateSection, p, &offset);
436: for (d = 0; d < dof; ++d) {
437: if (candidates[offset+d].index != claims[offset+d].index) {
438: key.i = candidates[offset+d].index; key.j = candidates[offset+d].rank;
439: PetscHashIJGet(roothash, key, &root);
440: if (root >= 0) {
441: vertices[0] = p; vertices[1] = localPoints[root];
442: DMPlexGetJoin(dm, 2, vertices, &joinSize, &join);
443: if (joinSize == 1) PetscHashIAdd(claimshash, join[0], offset+d);
444: DMPlexRestoreJoin(dm, 2, vertices, &joinSize, &join);
445: }
446: }
447: }
448: }
449: /* Create new pointSF from hashed claims */
450: PetscHashISize(claimshash, numLocalNew);
451: DMPlexGetChart(dm, &pStart, &pEnd);
452: PetscMalloc1(numLeaves + numLocalNew, &localPointsNew);
453: PetscMalloc1(numLeaves + numLocalNew, &remotePointsNew);
454: for (p = 0; p < numLeaves; ++p) {
455: localPointsNew[p] = localPoints[p];
456: remotePointsNew[p].index = remotePoints[p].index;
457: remotePointsNew[p].rank = remotePoints[p].rank;
458: }
459: p = numLeaves;
460: PetscHashIGetKeys(claimshash, &p, localPointsNew);
461: PetscSortInt(numLocalNew,&localPointsNew[numLeaves]);
462: for (p = numLeaves; p < numLeaves + numLocalNew; ++p) {
463: PetscHashIMap(claimshash, localPointsNew[p], offset);
464: remotePointsNew[p] = claims[offset];
465: }
466: PetscSFCreate(PetscObjectComm((PetscObject) dm), &sfPointNew);
467: PetscSFSetGraph(sfPointNew, pEnd-pStart, numLeaves+numLocalNew, localPointsNew, PETSC_OWN_POINTER, remotePointsNew, PETSC_OWN_POINTER);
468: DMSetPointSF(dm, sfPointNew);
469: PetscSFDestroy(&sfPointNew);
470: PetscHashIDestroy(claimshash);
471: }
472: PetscHashIDestroy(leafhash);
473: PetscHashIJDestroy(&roothash);
474: PetscSectionDestroy(&candidateSection);
475: PetscSectionDestroy(&candidateSectionRemote);
476: PetscSectionDestroy(&claimSection);
477: PetscFree(candidates);
478: PetscFree(candidatesRemote);
479: PetscFree(claims);
480: PetscLogEventEnd(DMPLEX_InterpolateSF,dm,0,0,0);
481: return(0);
482: }
484: /*@C
485: DMPlexInterpolate - Take in a cell-vertex mesh and return one with all intermediate faces, edges, etc.
487: Collective on DM
489: Input Parameters:
490: + dm - The DMPlex object with only cells and vertices
491: - dmInt - If NULL a new DM is created, otherwise the interpolated DM is put into the given DM
493: Output Parameter:
494: . dmInt - The complete DMPlex object
496: Level: intermediate
498: .keywords: mesh
499: .seealso: DMPlexUninterpolate(), DMPlexCreateFromCellList()
500: @*/
501: PetscErrorCode DMPlexInterpolate(DM dm, DM *dmInt)
502: {
503: DM idm, odm = dm;
504: PetscSF sfPoint;
505: PetscInt depth, dim, d;
509: PetscLogEventBegin(DMPLEX_Interpolate,dm,0,0,0);
510: DMPlexGetDepth(dm, &depth);
511: DMGetDimension(dm, &dim);
512: if (dim <= 1) {
513: PetscObjectReference((PetscObject) dm);
514: idm = dm;
515: }
516: for (d = 1; d < dim; ++d) {
517: /* Create interpolated mesh */
518: if ((d == dim-1) && *dmInt) {idm = *dmInt;}
519: else {DMCreate(PetscObjectComm((PetscObject)dm), &idm);}
520: DMSetType(idm, DMPLEX);
521: DMSetDimension(idm, dim);
522: if (depth > 0) {
523: DMPlexInterpolateFaces_Internal(odm, 1, idm);
524: DMGetPointSF(odm, &sfPoint);
525: DMPlexInterpolatePointSF(idm, sfPoint, depth);
526: }
527: if (odm != dm) {DMDestroy(&odm);}
528: odm = idm;
529: }
530: *dmInt = idm;
531: PetscLogEventEnd(DMPLEX_Interpolate,dm,0,0,0);
532: return(0);
533: }
535: /*@
536: DMPlexCopyCoordinates - Copy coordinates from one mesh to another with the same vertices
538: Collective on DM
540: Input Parameter:
541: . dmA - The DMPlex object with initial coordinates
543: Output Parameter:
544: . dmB - The DMPlex object with copied coordinates
546: Level: intermediate
548: Note: This is typically used when adding pieces other than vertices to a mesh
550: .keywords: mesh
551: .seealso: DMCopyLabels(), DMGetCoordinates(), DMGetCoordinatesLocal(), DMGetCoordinateDM(), DMGetCoordinateSection()
552: @*/
553: PetscErrorCode DMPlexCopyCoordinates(DM dmA, DM dmB)
554: {
555: Vec coordinatesA, coordinatesB;
556: PetscSection coordSectionA, coordSectionB;
557: PetscScalar *coordsA, *coordsB;
558: PetscInt spaceDim, vStartA, vStartB, vEndA, vEndB, coordSizeB, v, d;
562: if (dmA == dmB) return(0);
563: DMPlexGetDepthStratum(dmA, 0, &vStartA, &vEndA);
564: DMPlexGetDepthStratum(dmB, 0, &vStartB, &vEndB);
565: 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);
566: DMGetCoordinateSection(dmA, &coordSectionA);
567: DMGetCoordinateSection(dmB, &coordSectionB);
568: if (coordSectionA == coordSectionB) return(0);
569: if (!coordSectionB) {
570: PetscInt dim;
572: PetscSectionCreate(PetscObjectComm((PetscObject) coordSectionA), &coordSectionB);
573: DMGetCoordinateDim(dmA, &dim);
574: DMSetCoordinateSection(dmB, dim, coordSectionB);
575: PetscObjectDereference((PetscObject) coordSectionB);
576: }
577: PetscSectionSetNumFields(coordSectionB, 1);
578: PetscSectionGetFieldComponents(coordSectionA, 0, &spaceDim);
579: PetscSectionSetFieldComponents(coordSectionB, 0, spaceDim);
580: PetscSectionSetChart(coordSectionB, vStartB, vEndB);
581: for (v = vStartB; v < vEndB; ++v) {
582: PetscSectionSetDof(coordSectionB, v, spaceDim);
583: PetscSectionSetFieldDof(coordSectionB, v, 0, spaceDim);
584: }
585: PetscSectionSetUp(coordSectionB);
586: PetscSectionGetStorageSize(coordSectionB, &coordSizeB);
587: DMGetCoordinatesLocal(dmA, &coordinatesA);
588: VecCreate(PETSC_COMM_SELF, &coordinatesB);
589: PetscObjectSetName((PetscObject) coordinatesB, "coordinates");
590: VecSetSizes(coordinatesB, coordSizeB, PETSC_DETERMINE);
591: VecSetType(coordinatesB,VECSTANDARD);
592: VecGetArray(coordinatesA, &coordsA);
593: VecGetArray(coordinatesB, &coordsB);
594: for (v = 0; v < vEndB-vStartB; ++v) {
595: for (d = 0; d < spaceDim; ++d) {
596: coordsB[v*spaceDim+d] = coordsA[v*spaceDim+d];
597: }
598: }
599: VecRestoreArray(coordinatesA, &coordsA);
600: VecRestoreArray(coordinatesB, &coordsB);
601: DMSetCoordinatesLocal(dmB, coordinatesB);
602: VecDestroy(&coordinatesB);
603: return(0);
604: }
606: /*@
607: DMPlexUninterpolate - Take in a mesh with all intermediate faces, edges, etc. and return a cell-vertex mesh
609: Collective on DM
611: Input Parameter:
612: . dm - The complete DMPlex object
614: Output Parameter:
615: . dmUnint - The DMPlex object with only cells and vertices
617: Level: intermediate
619: .keywords: mesh
620: .seealso: DMPlexInterpolate(), DMPlexCreateFromCellList()
621: @*/
622: PetscErrorCode DMPlexUninterpolate(DM dm, DM *dmUnint)
623: {
624: DM udm;
625: PetscInt dim, vStart, vEnd, cStart, cEnd, c, maxConeSize = 0, *cone;
629: DMGetDimension(dm, &dim);
630: if (dim <= 1) {
631: PetscObjectReference((PetscObject) dm);
632: *dmUnint = dm;
633: return(0);
634: }
635: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
636: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
637: DMCreate(PetscObjectComm((PetscObject) dm), &udm);
638: DMSetType(udm, DMPLEX);
639: DMSetDimension(udm, dim);
640: DMPlexSetChart(udm, cStart, vEnd);
641: for (c = cStart; c < cEnd; ++c) {
642: PetscInt *closure = NULL, closureSize, cl, coneSize = 0;
644: DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
645: for (cl = 0; cl < closureSize*2; cl += 2) {
646: const PetscInt p = closure[cl];
648: if ((p >= vStart) && (p < vEnd)) ++coneSize;
649: }
650: DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
651: DMPlexSetConeSize(udm, c, coneSize);
652: maxConeSize = PetscMax(maxConeSize, coneSize);
653: }
654: DMSetUp(udm);
655: PetscMalloc1(maxConeSize, &cone);
656: for (c = cStart; c < cEnd; ++c) {
657: PetscInt *closure = NULL, closureSize, cl, coneSize = 0;
659: DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
660: for (cl = 0; cl < closureSize*2; cl += 2) {
661: const PetscInt p = closure[cl];
663: if ((p >= vStart) && (p < vEnd)) cone[coneSize++] = p;
664: }
665: DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
666: DMPlexSetCone(udm, c, cone);
667: }
668: PetscFree(cone);
669: DMPlexSymmetrize(udm);
670: DMPlexStratify(udm);
671: /* Reduce SF */
672: {
673: PetscSF sfPoint, sfPointUn;
674: const PetscSFNode *remotePoints;
675: const PetscInt *localPoints;
676: PetscSFNode *remotePointsUn;
677: PetscInt *localPointsUn;
678: PetscInt vEnd, numRoots, numLeaves, l;
679: PetscInt numLeavesUn = 0, n = 0;
680: PetscErrorCode ierr;
682: /* Get original SF information */
683: DMGetPointSF(dm, &sfPoint);
684: DMGetPointSF(udm, &sfPointUn);
685: DMPlexGetDepthStratum(dm, 0, NULL, &vEnd);
686: PetscSFGetGraph(sfPoint, &numRoots, &numLeaves, &localPoints, &remotePoints);
687: /* Allocate space for cells and vertices */
688: for (l = 0; l < numLeaves; ++l) if (localPoints[l] < vEnd) numLeavesUn++;
689: /* Fill in leaves */
690: if (vEnd >= 0) {
691: PetscMalloc1(numLeavesUn, &remotePointsUn);
692: PetscMalloc1(numLeavesUn, &localPointsUn);
693: for (l = 0; l < numLeaves; l++) {
694: if (localPoints[l] < vEnd) {
695: localPointsUn[n] = localPoints[l];
696: remotePointsUn[n].rank = remotePoints[l].rank;
697: remotePointsUn[n].index = remotePoints[l].index;
698: ++n;
699: }
700: }
701: if (n != numLeavesUn) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Inconsistent number of leaves %d != %d", n, numLeavesUn);
702: PetscSFSetGraph(sfPointUn, vEnd, numLeavesUn, localPointsUn, PETSC_OWN_POINTER, remotePointsUn, PETSC_OWN_POINTER);
703: }
704: }
705: *dmUnint = udm;
706: return(0);
707: }