Actual source code: plexinterpolate.c
petsc-3.9.4 2018-09-11
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, MPIU_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)), MPIU_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 - The interpolated DM
493: Output Parameter:
494: . dmInt - The complete DMPlex object
496: Level: intermediate
498: Notes: It does not copy over the coordinates.
500: .keywords: mesh
501: .seealso: DMPlexUninterpolate(), DMPlexCreateFromCellList(), DMPlexCopyCoordinates()
502: @*/
503: PetscErrorCode DMPlexInterpolate(DM dm, DM *dmInt)
504: {
505: DM idm, odm = dm;
506: PetscSF sfPoint;
507: PetscInt depth, dim, d;
508: const char *name;
514: PetscLogEventBegin(DMPLEX_Interpolate,dm,0,0,0);
515: DMPlexGetDepth(dm, &depth);
516: DMGetDimension(dm, &dim);
517: if ((depth == dim) || (dim <= 1)) {
518: PetscObjectReference((PetscObject) dm);
519: idm = dm;
520: } else {
521: for (d = 1; d < dim; ++d) {
522: /* Create interpolated mesh */
523: DMCreate(PetscObjectComm((PetscObject)dm), &idm);
524: DMSetType(idm, DMPLEX);
525: DMSetDimension(idm, dim);
526: if (depth > 0) {
527: DMPlexInterpolateFaces_Internal(odm, 1, idm);
528: DMGetPointSF(odm, &sfPoint);
529: DMPlexInterpolatePointSF(idm, sfPoint, depth);
530: }
531: if (odm != dm) {DMDestroy(&odm);}
532: odm = idm;
533: }
534: PetscObjectGetName((PetscObject) dm, &name);
535: PetscObjectSetName((PetscObject) idm, name);
536: DMPlexCopyCoordinates(dm, idm);
537: DMCopyLabels(dm, idm);
538: }
539: {
540: PetscBool isper;
541: const PetscReal *maxCell, *L;
542: const DMBoundaryType *bd;
544: DMGetPeriodicity(dm,&isper,&maxCell,&L,&bd);
545: DMSetPeriodicity(idm,isper,maxCell,L,bd);
546: }
547: *dmInt = idm;
548: PetscLogEventEnd(DMPLEX_Interpolate,dm,0,0,0);
549: return(0);
550: }
552: /*@
553: DMPlexCopyCoordinates - Copy coordinates from one mesh to another with the same vertices
555: Collective on DM
557: Input Parameter:
558: . dmA - The DMPlex object with initial coordinates
560: Output Parameter:
561: . dmB - The DMPlex object with copied coordinates
563: Level: intermediate
565: Note: This is typically used when adding pieces other than vertices to a mesh
567: .keywords: mesh
568: .seealso: DMCopyLabels(), DMGetCoordinates(), DMGetCoordinatesLocal(), DMGetCoordinateDM(), DMGetCoordinateSection()
569: @*/
570: PetscErrorCode DMPlexCopyCoordinates(DM dmA, DM dmB)
571: {
572: Vec coordinatesA, coordinatesB;
573: VecType vtype;
574: PetscSection coordSectionA, coordSectionB;
575: PetscScalar *coordsA, *coordsB;
576: PetscInt spaceDim, Nf, vStartA, vStartB, vEndA, vEndB, coordSizeB, v, d;
577: PetscInt cStartA, cEndA, cStartB, cEndB, cS, cE;
578: PetscBool lc = PETSC_FALSE;
584: if (dmA == dmB) return(0);
585: DMPlexGetDepthStratum(dmA, 0, &vStartA, &vEndA);
586: DMPlexGetDepthStratum(dmB, 0, &vStartB, &vEndB);
587: 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);
588: DMPlexGetHeightStratum(dmA, 0, &cStartA, &cEndA);
589: DMPlexGetHeightStratum(dmB, 0, &cStartB, &cEndB);
590: DMGetCoordinateSection(dmA, &coordSectionA);
591: DMGetCoordinateSection(dmB, &coordSectionB);
592: if (coordSectionA == coordSectionB) return(0);
593: PetscSectionGetNumFields(coordSectionA, &Nf);
594: if (!Nf) return(0);
595: if (Nf > 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "The number of coordinate fields must be 1, not %D", Nf);
596: if (!coordSectionB) {
597: PetscInt dim;
599: PetscSectionCreate(PetscObjectComm((PetscObject) coordSectionA), &coordSectionB);
600: DMGetCoordinateDim(dmA, &dim);
601: DMSetCoordinateSection(dmB, dim, coordSectionB);
602: PetscObjectDereference((PetscObject) coordSectionB);
603: }
604: PetscSectionSetNumFields(coordSectionB, 1);
605: PetscSectionGetFieldComponents(coordSectionA, 0, &spaceDim);
606: PetscSectionSetFieldComponents(coordSectionB, 0, spaceDim);
607: PetscSectionGetChart(coordSectionA, &cS, &cE);
608: if (cStartA <= cS && cS < cEndA) { /* localized coordinates */
609: 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);
610: cS = cS - cStartA + cStartB;
611: cE = vEndB;
612: lc = PETSC_TRUE;
613: } else {
614: cS = vStartB;
615: cE = vEndB;
616: }
617: PetscSectionSetChart(coordSectionB, cS, cE);
618: for (v = vStartB; v < vEndB; ++v) {
619: PetscSectionSetDof(coordSectionB, v, spaceDim);
620: PetscSectionSetFieldDof(coordSectionB, v, 0, spaceDim);
621: }
622: if (lc) { /* localized coordinates */
623: PetscInt c;
625: for (c = cS-cStartB; c < cEndB-cStartB; c++) {
626: PetscInt dof;
628: PetscSectionGetDof(coordSectionA, c + cStartA, &dof);
629: PetscSectionSetDof(coordSectionB, c + cStartB, dof);
630: PetscSectionSetFieldDof(coordSectionB, c + cStartB, 0, dof);
631: }
632: }
633: PetscSectionSetUp(coordSectionB);
634: PetscSectionGetStorageSize(coordSectionB, &coordSizeB);
635: DMGetCoordinatesLocal(dmA, &coordinatesA);
636: VecCreate(PETSC_COMM_SELF, &coordinatesB);
637: PetscObjectSetName((PetscObject) coordinatesB, "coordinates");
638: VecSetSizes(coordinatesB, coordSizeB, PETSC_DETERMINE);
639: VecGetBlockSize(coordinatesA, &d);
640: VecSetBlockSize(coordinatesB, d);
641: VecGetType(coordinatesA, &vtype);
642: VecSetType(coordinatesB, vtype);
643: VecGetArray(coordinatesA, &coordsA);
644: VecGetArray(coordinatesB, &coordsB);
645: for (v = 0; v < vEndB-vStartB; ++v) {
646: PetscInt offA, offB;
648: PetscSectionGetOffset(coordSectionA, v + vStartA, &offA);
649: PetscSectionGetOffset(coordSectionB, v + vStartB, &offB);
650: for (d = 0; d < spaceDim; ++d) {
651: coordsB[offB+d] = coordsA[offA+d];
652: }
653: }
654: if (lc) { /* localized coordinates */
655: PetscInt c;
657: for (c = cS-cStartB; c < cEndB-cStartB; c++) {
658: PetscInt dof, offA, offB;
660: PetscSectionGetOffset(coordSectionA, c + cStartA, &offA);
661: PetscSectionGetOffset(coordSectionB, c + cStartB, &offB);
662: PetscSectionGetDof(coordSectionA, c + cStartA, &dof);
663: PetscMemcpy(coordsB + offB,coordsA + offA,dof*sizeof(*coordsB));
664: }
665: }
666: VecRestoreArray(coordinatesA, &coordsA);
667: VecRestoreArray(coordinatesB, &coordsB);
668: DMSetCoordinatesLocal(dmB, coordinatesB);
669: VecDestroy(&coordinatesB);
670: return(0);
671: }
673: /*@
674: DMPlexUninterpolate - Take in a mesh with all intermediate faces, edges, etc. and return a cell-vertex mesh
676: Collective on DM
678: Input Parameter:
679: . dm - The complete DMPlex object
681: Output Parameter:
682: . dmUnint - The DMPlex object with only cells and vertices
684: Level: intermediate
686: Notes: It does not copy over the coordinates.
688: .keywords: mesh
689: .seealso: DMPlexInterpolate(), DMPlexCreateFromCellList(), DMPlexCopyCoordinates()
690: @*/
691: PetscErrorCode DMPlexUninterpolate(DM dm, DM *dmUnint)
692: {
693: DM udm;
694: PetscInt dim, vStart, vEnd, cStart, cEnd, c, maxConeSize = 0, *cone;
700: DMGetDimension(dm, &dim);
701: if (dim <= 1) {
702: PetscObjectReference((PetscObject) dm);
703: *dmUnint = dm;
704: return(0);
705: }
706: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
707: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
708: DMCreate(PetscObjectComm((PetscObject) dm), &udm);
709: DMSetType(udm, DMPLEX);
710: DMSetDimension(udm, dim);
711: DMPlexSetChart(udm, cStart, vEnd);
712: for (c = cStart; c < cEnd; ++c) {
713: PetscInt *closure = NULL, closureSize, cl, coneSize = 0;
715: DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
716: for (cl = 0; cl < closureSize*2; cl += 2) {
717: const PetscInt p = closure[cl];
719: if ((p >= vStart) && (p < vEnd)) ++coneSize;
720: }
721: DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
722: DMPlexSetConeSize(udm, c, coneSize);
723: maxConeSize = PetscMax(maxConeSize, coneSize);
724: }
725: DMSetUp(udm);
726: PetscMalloc1(maxConeSize, &cone);
727: for (c = cStart; c < cEnd; ++c) {
728: PetscInt *closure = NULL, closureSize, cl, coneSize = 0;
730: DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
731: for (cl = 0; cl < closureSize*2; cl += 2) {
732: const PetscInt p = closure[cl];
734: if ((p >= vStart) && (p < vEnd)) cone[coneSize++] = p;
735: }
736: DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
737: DMPlexSetCone(udm, c, cone);
738: }
739: PetscFree(cone);
740: DMPlexSymmetrize(udm);
741: DMPlexStratify(udm);
742: /* Reduce SF */
743: {
744: PetscSF sfPoint, sfPointUn;
745: const PetscSFNode *remotePoints;
746: const PetscInt *localPoints;
747: PetscSFNode *remotePointsUn;
748: PetscInt *localPointsUn;
749: PetscInt vEnd, numRoots, numLeaves, l;
750: PetscInt numLeavesUn = 0, n = 0;
751: PetscErrorCode ierr;
753: /* Get original SF information */
754: DMGetPointSF(dm, &sfPoint);
755: DMGetPointSF(udm, &sfPointUn);
756: DMPlexGetDepthStratum(dm, 0, NULL, &vEnd);
757: PetscSFGetGraph(sfPoint, &numRoots, &numLeaves, &localPoints, &remotePoints);
758: /* Allocate space for cells and vertices */
759: for (l = 0; l < numLeaves; ++l) if (localPoints[l] < vEnd) numLeavesUn++;
760: /* Fill in leaves */
761: if (vEnd >= 0) {
762: PetscMalloc1(numLeavesUn, &remotePointsUn);
763: PetscMalloc1(numLeavesUn, &localPointsUn);
764: for (l = 0; l < numLeaves; l++) {
765: if (localPoints[l] < vEnd) {
766: localPointsUn[n] = localPoints[l];
767: remotePointsUn[n].rank = remotePoints[l].rank;
768: remotePointsUn[n].index = remotePoints[l].index;
769: ++n;
770: }
771: }
772: if (n != numLeavesUn) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Inconsistent number of leaves %d != %d", n, numLeavesUn);
773: PetscSFSetGraph(sfPointUn, vEnd, numLeavesUn, localPointsUn, PETSC_OWN_POINTER, remotePointsUn, PETSC_OWN_POINTER);
774: }
775: }
776: {
777: PetscBool isper;
778: const PetscReal *maxCell, *L;
779: const DMBoundaryType *bd;
781: DMGetPeriodicity(dm,&isper,&maxCell,&L,&bd);
782: DMSetPeriodicity(udm,isper,maxCell,L,bd);
783: }
785: *dmUnint = udm;
786: return(0);
787: }