Actual source code: plexinterpolate.c
petsc-3.5.4 2015-05-23
1: #include <petsc-private/dmpleximpl.h> /*I "petscdmplex.h" I*/
2: #include <../src/sys/utils/hash.h>
6: /*
7: DMPlexGetFaces_Internal - Gets groups of vertices that correspond to faces for the given cell
8: */
9: PetscErrorCode DMPlexGetFaces_Internal(DM dm, PetscInt dim, PetscInt p, PetscInt *numFaces, PetscInt *faceSize, const PetscInt *faces[])
10: {
11: const PetscInt *cone = NULL;
12: PetscInt maxConeSize, maxSupportSize, coneSize;
13: PetscErrorCode ierr;
17: DMPlexGetMaxSizes(dm, &maxConeSize, &maxSupportSize);
18: DMPlexGetConeSize(dm, p, &coneSize);
19: DMPlexGetCone(dm, p, &cone);
20: DMPlexGetRawFaces_Internal(dm, dim, coneSize, cone, numFaces, faceSize, faces);
21: return(0);
22: }
26: /*
27: DMPlexRestoreFaces_Internal - Restores the array
28: */
29: PetscErrorCode DMPlexRestoreFaces_Internal(DM dm, PetscInt dim, PetscInt p, PetscInt *numFaces, PetscInt *faceSize, const PetscInt *faces[])
30: {
31: PetscErrorCode ierr;
34: DMRestoreWorkArray(dm, 0, PETSC_INT, (void *) faces);
35: return(0);
36: }
40: /*
41: DMPlexGetRawFaces_Internal - Gets groups of vertices that correspond to faces for the given cone
42: */
43: PetscErrorCode DMPlexGetRawFaces_Internal(DM dm, PetscInt dim, PetscInt coneSize, const PetscInt cone[], PetscInt *numFaces, PetscInt *faceSize, const PetscInt *faces[])
44: {
45: PetscInt *facesTmp;
46: PetscInt maxConeSize, maxSupportSize;
47: PetscErrorCode ierr;
51: DMPlexGetMaxSizes(dm, &maxConeSize, &maxSupportSize);
52: if (faces) {DMGetWorkArray(dm, PetscSqr(PetscMax(maxConeSize, maxSupportSize)), PETSC_INT, &facesTmp);}
53: switch (dim) {
54: case 1:
55: switch (coneSize) {
56: case 2:
57: if (faces) {
58: facesTmp[0] = cone[0]; facesTmp[1] = cone[1];
59: *faces = facesTmp;
60: }
61: if (numFaces) *numFaces = 2;
62: if (faceSize) *faceSize = 1;
63: break;
64: default:
65: SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Cone size %D not supported for dimension %D", coneSize, dim);
66: }
67: break;
68: case 2:
69: switch (coneSize) {
70: case 3:
71: if (faces) {
72: facesTmp[0] = cone[0]; facesTmp[1] = cone[1];
73: facesTmp[2] = cone[1]; facesTmp[3] = cone[2];
74: facesTmp[4] = cone[2]; facesTmp[5] = cone[0];
75: *faces = facesTmp;
76: }
77: if (numFaces) *numFaces = 3;
78: if (faceSize) *faceSize = 2;
79: break;
80: case 4:
81: /* Vertices follow right hand rule */
82: if (faces) {
83: facesTmp[0] = cone[0]; facesTmp[1] = cone[1];
84: facesTmp[2] = cone[1]; facesTmp[3] = cone[2];
85: facesTmp[4] = cone[2]; facesTmp[5] = cone[3];
86: facesTmp[6] = cone[3]; facesTmp[7] = cone[0];
87: *faces = facesTmp;
88: }
89: if (numFaces) *numFaces = 4;
90: if (faceSize) *faceSize = 2;
91: break;
92: default:
93: SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Cone size %D not supported for dimension %D", coneSize, dim);
94: }
95: break;
96: case 3:
97: switch (coneSize) {
98: case 3:
99: if (faces) {
100: facesTmp[0] = cone[0]; facesTmp[1] = cone[1];
101: facesTmp[2] = cone[1]; facesTmp[3] = cone[2];
102: facesTmp[4] = cone[2]; facesTmp[5] = cone[0];
103: *faces = facesTmp;
104: }
105: if (numFaces) *numFaces = 3;
106: if (faceSize) *faceSize = 2;
107: break;
108: case 4:
109: /* Vertices of first face follow right hand rule and normal points away from last vertex */
110: if (faces) {
111: facesTmp[0] = cone[0]; facesTmp[1] = cone[1]; facesTmp[2] = cone[2];
112: facesTmp[3] = cone[0]; facesTmp[4] = cone[3]; facesTmp[5] = cone[1];
113: facesTmp[6] = cone[0]; facesTmp[7] = cone[2]; facesTmp[8] = cone[3];
114: facesTmp[9] = cone[2]; facesTmp[10] = cone[1]; facesTmp[11] = cone[3];
115: *faces = facesTmp;
116: }
117: if (numFaces) *numFaces = 4;
118: if (faceSize) *faceSize = 3;
119: break;
120: case 8:
121: if (faces) {
122: facesTmp[0] = cone[0]; facesTmp[1] = cone[1]; facesTmp[2] = cone[2]; facesTmp[3] = cone[3]; /* Bottom */
123: facesTmp[4] = cone[4]; facesTmp[5] = cone[5]; facesTmp[6] = cone[6]; facesTmp[7] = cone[7]; /* Top */
124: facesTmp[8] = cone[0]; facesTmp[9] = cone[3]; facesTmp[10] = cone[5]; facesTmp[11] = cone[4]; /* Front */
125: facesTmp[12] = cone[2]; facesTmp[13] = cone[1]; facesTmp[14] = cone[7]; facesTmp[15] = cone[6]; /* Back */
126: facesTmp[16] = cone[3]; facesTmp[17] = cone[2]; facesTmp[18] = cone[6]; facesTmp[19] = cone[5]; /* Right */
127: facesTmp[20] = cone[0]; facesTmp[21] = cone[4]; facesTmp[22] = cone[7]; facesTmp[23] = cone[1]; /* Left */
128: *faces = facesTmp;
129: }
130: if (numFaces) *numFaces = 6;
131: if (faceSize) *faceSize = 4;
132: break;
133: default:
134: SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Cone size %D not supported for dimension %D", coneSize, dim);
135: }
136: break;
137: default:
138: SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Dimension %D not supported", dim);
139: }
140: return(0);
141: }
145: /* This interpolates faces for cells at some stratum */
146: static PetscErrorCode DMPlexInterpolateFaces_Internal(DM dm, PetscInt cellDepth, DM idm)
147: {
148: DMLabel subpointMap;
149: PetscHashIJKL faceTable;
150: PetscInt *pStart, *pEnd;
151: PetscInt cellDim, depth, faceDepth = cellDepth, numPoints = 0, faceSizeAll = 0, face, c, d;
155: DMPlexGetDimension(dm, &cellDim);
156: /* HACK: I need a better way to determine face dimension, or an alternative to GetFaces() */
157: DMPlexGetSubpointMap(dm, &subpointMap);
158: if (subpointMap) ++cellDim;
159: DMPlexGetDepth(dm, &depth);
160: ++depth;
161: ++cellDepth;
162: cellDim -= depth - cellDepth;
163: PetscMalloc2(depth+1,&pStart,depth+1,&pEnd);
164: for (d = depth-1; d >= faceDepth; --d) {
165: DMPlexGetDepthStratum(dm, d, &pStart[d+1], &pEnd[d+1]);
166: }
167: DMPlexGetDepthStratum(dm, -1, NULL, &pStart[faceDepth]);
168: pEnd[faceDepth] = pStart[faceDepth];
169: for (d = faceDepth-1; d >= 0; --d) {
170: DMPlexGetDepthStratum(dm, d, &pStart[d], &pEnd[d]);
171: }
172: if (pEnd[cellDepth] > pStart[cellDepth]) {DMPlexGetFaces_Internal(dm, cellDim, pStart[cellDepth], NULL, &faceSizeAll, NULL);}
173: if (faceSizeAll > 4) SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "Do not support interpolation of meshes with faces of %D vertices", faceSizeAll);
174: PetscHashIJKLCreate(&faceTable);
175: for (c = pStart[cellDepth], face = pStart[faceDepth]; c < pEnd[cellDepth]; ++c) {
176: const PetscInt *cellFaces;
177: PetscInt numCellFaces, faceSize, cf;
179: DMPlexGetFaces_Internal(dm, cellDim, c, &numCellFaces, &faceSize, &cellFaces);
180: if (faceSize != faceSizeAll) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Inconsistent face for cell %D of size %D != %D", c, faceSize, faceSizeAll);
181: for (cf = 0; cf < numCellFaces; ++cf) {
182: const PetscInt *cellFace = &cellFaces[cf*faceSize];
183: PetscHashIJKLKey key;
184: PetscHashIJKLIter missing, iter;
186: if (faceSize == 2) {
187: key.i = PetscMin(cellFace[0], cellFace[1]);
188: key.j = PetscMax(cellFace[0], cellFace[1]);
189: key.k = 0;
190: key.l = 0;
191: } else {
192: key.i = cellFace[0]; key.j = cellFace[1]; key.k = cellFace[2]; key.l = faceSize > 3 ? cellFace[3] : 0;
193: PetscSortInt(faceSize, (PetscInt *) &key);
194: }
195: PetscHashIJKLPut(faceTable, key, &missing, &iter);
196: if (missing) {PetscHashIJKLSet(faceTable, iter, face++);}
197: }
198: DMPlexRestoreFaces_Internal(dm, cellDim, c, &numCellFaces, &faceSize, &cellFaces);
199: }
200: pEnd[faceDepth] = face;
201: PetscHashIJKLDestroy(&faceTable);
202: /* Count new points */
203: for (d = 0; d <= depth; ++d) {
204: numPoints += pEnd[d]-pStart[d];
205: }
206: DMPlexSetChart(idm, 0, numPoints);
207: /* Set cone sizes */
208: for (d = 0; d <= depth; ++d) {
209: PetscInt coneSize, p;
211: if (d == faceDepth) {
212: for (p = pStart[d]; p < pEnd[d]; ++p) {
213: /* I see no way to do this if we admit faces of different shapes */
214: DMPlexSetConeSize(idm, p, faceSizeAll);
215: }
216: } else if (d == cellDepth) {
217: for (p = pStart[d]; p < pEnd[d]; ++p) {
218: /* Number of cell faces may be different from number of cell vertices*/
219: DMPlexGetFaces_Internal(dm, cellDim, p, &coneSize, NULL, NULL);
220: DMPlexSetConeSize(idm, p, coneSize);
221: }
222: } else {
223: for (p = pStart[d]; p < pEnd[d]; ++p) {
224: DMPlexGetConeSize(dm, p, &coneSize);
225: DMPlexSetConeSize(idm, p, coneSize);
226: }
227: }
228: }
229: DMSetUp(idm);
230: /* Get face cones from subsets of cell vertices */
231: if (faceSizeAll > 4) SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "Do not support interpolation of meshes with faces of %D vertices", faceSizeAll);
232: PetscHashIJKLCreate(&faceTable);
233: for (d = depth; d > cellDepth; --d) {
234: const PetscInt *cone;
235: PetscInt p;
237: for (p = pStart[d]; p < pEnd[d]; ++p) {
238: DMPlexGetCone(dm, p, &cone);
239: DMPlexSetCone(idm, p, cone);
240: DMPlexGetConeOrientation(dm, p, &cone);
241: DMPlexSetConeOrientation(idm, p, cone);
242: }
243: }
244: for (c = pStart[cellDepth], face = pStart[faceDepth]; c < pEnd[cellDepth]; ++c) {
245: const PetscInt *cellFaces;
246: PetscInt numCellFaces, faceSize, cf;
248: DMPlexGetFaces_Internal(dm, cellDim, c, &numCellFaces, &faceSize, &cellFaces);
249: if (faceSize != faceSizeAll) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Inconsistent face for cell %D of size %D != %D", c, faceSize, faceSizeAll);
250: for (cf = 0; cf < numCellFaces; ++cf) {
251: const PetscInt *cellFace = &cellFaces[cf*faceSize];
252: PetscHashIJKLKey key;
253: PetscHashIJKLIter missing, iter;
255: if (faceSize == 2) {
256: key.i = PetscMin(cellFace[0], cellFace[1]);
257: key.j = PetscMax(cellFace[0], cellFace[1]);
258: key.k = 0;
259: key.l = 0;
260: } else {
261: key.i = cellFace[0]; key.j = cellFace[1]; key.k = cellFace[2]; key.l = faceSize > 3 ? cellFace[3] : 0;
262: PetscSortInt(faceSize, (PetscInt *) &key);
263: }
264: PetscHashIJKLPut(faceTable, key, &missing, &iter);
265: if (missing) {
266: DMPlexSetCone(idm, face, cellFace);
267: PetscHashIJKLSet(faceTable, iter, face);
268: DMPlexInsertCone(idm, c, cf, face++);
269: } else {
270: const PetscInt *cone;
271: PetscInt coneSize, ornt, i, j, f;
273: PetscHashIJKLGet(faceTable, iter, &f);
274: DMPlexInsertCone(idm, c, cf, f);
275: /* Orient face: Do not allow reverse orientation at the first vertex */
276: DMPlexGetConeSize(idm, f, &coneSize);
277: DMPlexGetCone(idm, f, &cone);
278: 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);
279: /* - First find the initial vertex */
280: for (i = 0; i < faceSize; ++i) if (cellFace[0] == cone[i]) break;
281: /* - Try forward comparison */
282: for (j = 0; j < faceSize; ++j) if (cellFace[j] != cone[(i+j)%faceSize]) break;
283: if (j == faceSize) {
284: if ((faceSize == 2) && (i == 1)) ornt = -2;
285: else ornt = i;
286: } else {
287: /* - Try backward comparison */
288: for (j = 0; j < faceSize; ++j) if (cellFace[j] != cone[(i+faceSize-j)%faceSize]) break;
289: if (j == faceSize) {
290: if (i == 0) ornt = -faceSize;
291: else ornt = -i;
292: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Could not determine face orientation");
293: }
294: DMPlexInsertConeOrientation(idm, c, cf, ornt);
295: }
296: }
297: DMPlexRestoreFaces_Internal(dm, cellDim, c, &numCellFaces, &faceSize, &cellFaces);
298: }
299: 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]);
300: PetscFree2(pStart,pEnd);
301: PetscHashIJKLDestroy(&faceTable);
302: PetscFree2(pStart,pEnd);
303: DMPlexSymmetrize(idm);
304: DMPlexStratify(idm);
305: return(0);
306: }
310: /*@C
311: DMPlexInterpolate - Take in a cell-vertex mesh and return one with all intermediate faces, edges, etc.
313: Collective on DM
315: Input Parameters:
316: + dm - The DMPlex object with only cells and vertices
317: - dmInt - If NULL a new DM is created, otherwise the interpolated DM is put into the given DM
319: Output Parameter:
320: . dmInt - The complete DMPlex object
322: Level: intermediate
324: .keywords: mesh
325: .seealso: DMPlexUninterpolate(), DMPlexCreateFromCellList()
326: @*/
327: PetscErrorCode DMPlexInterpolate(DM dm, DM *dmInt)
328: {
329: DM idm, odm = dm;
330: PetscInt depth, dim, d;
334: PetscLogEventBegin(DMPLEX_Interpolate,dm,0,0,0);
335: DMPlexGetDepth(dm, &depth);
336: DMPlexGetDimension(dm, &dim);
337: if (dim <= 1) {
338: PetscObjectReference((PetscObject) dm);
339: idm = dm;
340: }
341: for (d = 1; d < dim; ++d) {
342: /* Create interpolated mesh */
343: if ((d == dim-1) && *dmInt) {idm = *dmInt;}
344: else {DMCreate(PetscObjectComm((PetscObject)dm), &idm);}
345: DMSetType(idm, DMPLEX);
346: DMPlexSetDimension(idm, dim);
347: if (depth > 0) {DMPlexInterpolateFaces_Internal(odm, 1, idm);}
348: if (odm != dm) {DMDestroy(&odm);}
349: odm = idm;
350: }
351: *dmInt = idm;
352: PetscLogEventEnd(DMPLEX_Interpolate,dm,0,0,0);
353: return(0);
354: }
358: /*@
359: DMPlexCopyCoordinates - Copy coordinates from one mesh to another with the same vertices
361: Collective on DM
363: Input Parameter:
364: . dmA - The DMPlex object with initial coordinates
366: Output Parameter:
367: . dmB - The DMPlex object with copied coordinates
369: Level: intermediate
371: Note: This is typically used when adding pieces other than vertices to a mesh
373: .keywords: mesh
374: .seealso: DMCopyLabels(), DMGetCoordinates(), DMGetCoordinatesLocal(), DMGetCoordinateDM(), DMGetCoordinateSection()
375: @*/
376: PetscErrorCode DMPlexCopyCoordinates(DM dmA, DM dmB)
377: {
378: Vec coordinatesA, coordinatesB;
379: PetscSection coordSectionA, coordSectionB;
380: PetscScalar *coordsA, *coordsB;
381: PetscInt spaceDim, vStartA, vStartB, vEndA, vEndB, coordSizeB, v, d;
385: if (dmA == dmB) return(0);
386: DMPlexGetDepthStratum(dmA, 0, &vStartA, &vEndA);
387: DMPlexGetDepthStratum(dmB, 0, &vStartB, &vEndB);
388: 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);
389: DMGetCoordinateSection(dmA, &coordSectionA);
390: DMGetCoordinateSection(dmB, &coordSectionB);
391: PetscSectionSetNumFields(coordSectionB, 1);
392: PetscSectionGetFieldComponents(coordSectionA, 0, &spaceDim);
393: PetscSectionSetFieldComponents(coordSectionB, 0, spaceDim);
394: PetscSectionSetChart(coordSectionB, vStartB, vEndB);
395: for (v = vStartB; v < vEndB; ++v) {
396: PetscSectionSetDof(coordSectionB, v, spaceDim);
397: PetscSectionSetFieldDof(coordSectionB, v, 0, spaceDim);
398: }
399: PetscSectionSetUp(coordSectionB);
400: PetscSectionGetStorageSize(coordSectionB, &coordSizeB);
401: DMGetCoordinatesLocal(dmA, &coordinatesA);
402: VecCreate(PetscObjectComm((PetscObject) dmB), &coordinatesB);
403: PetscObjectSetName((PetscObject) coordinatesB, "coordinates");
404: VecSetSizes(coordinatesB, coordSizeB, PETSC_DETERMINE);
405: VecSetType(coordinatesB,VECSTANDARD);
406: VecGetArray(coordinatesA, &coordsA);
407: VecGetArray(coordinatesB, &coordsB);
408: for (v = 0; v < vEndB-vStartB; ++v) {
409: for (d = 0; d < spaceDim; ++d) {
410: coordsB[v*spaceDim+d] = coordsA[v*spaceDim+d];
411: }
412: }
413: VecRestoreArray(coordinatesA, &coordsA);
414: VecRestoreArray(coordinatesB, &coordsB);
415: DMSetCoordinatesLocal(dmB, coordinatesB);
416: VecDestroy(&coordinatesB);
417: return(0);
418: }
422: /*@
423: DMPlexCopyLabels - Copy labels from one mesh to another with a superset of the points
425: Collective on DM
427: Input Parameter:
428: . dmA - The DMPlex object with initial labels
430: Output Parameter:
431: . dmB - The DMPlex object with copied labels
433: Level: intermediate
435: Note: This is typically used when interpolating or otherwise adding to a mesh
437: .keywords: mesh
438: .seealso: DMCopyCoordinates(), DMGetCoordinates(), DMGetCoordinatesLocal(), DMGetCoordinateDM(), DMGetCoordinateSection()
439: @*/
440: PetscErrorCode DMPlexCopyLabels(DM dmA, DM dmB)
441: {
442: PetscInt numLabels, l;
446: if (dmA == dmB) return(0);
447: DMPlexGetNumLabels(dmA, &numLabels);
448: for (l = 0; l < numLabels; ++l) {
449: DMLabel label, labelNew;
450: const char *name;
451: PetscBool flg;
453: DMPlexGetLabelName(dmA, l, &name);
454: PetscStrcmp(name, "depth", &flg);
455: if (flg) continue;
456: DMPlexGetLabel(dmA, name, &label);
457: DMLabelDuplicate(label, &labelNew);
458: DMPlexAddLabel(dmB, labelNew);
459: }
460: return(0);
461: }
465: /*@
466: DMPlexUninterpolate - Take in a mesh with all intermediate faces, edges, etc. and return a cell-vertex mesh
468: Collective on DM
470: Input Parameter:
471: . dm - The complete DMPlex object
473: Output Parameter:
474: . dmUnint - The DMPlex object with only cells and vertices
476: Level: intermediate
478: .keywords: mesh
479: .seealso: DMPlexInterpolate(), DMPlexCreateFromCellList()
480: @*/
481: PetscErrorCode DMPlexUninterpolate(DM dm, DM *dmUnint)
482: {
483: DM udm;
484: PetscInt dim, vStart, vEnd, cStart, cEnd, c, maxConeSize = 0, *cone;
488: DMPlexGetDimension(dm, &dim);
489: if (dim <= 1) {
490: PetscObjectReference((PetscObject) dm);
491: *dmUnint = dm;
492: return(0);
493: }
494: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
495: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
496: DMCreate(PetscObjectComm((PetscObject) dm), &udm);
497: DMSetType(udm, DMPLEX);
498: DMPlexSetDimension(udm, dim);
499: DMPlexSetChart(udm, cStart, vEnd);
500: for (c = cStart; c < cEnd; ++c) {
501: PetscInt *closure = NULL, closureSize, cl, coneSize = 0;
503: DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
504: for (cl = 0; cl < closureSize*2; cl += 2) {
505: const PetscInt p = closure[cl];
507: if ((p >= vStart) && (p < vEnd)) ++coneSize;
508: }
509: DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
510: DMPlexSetConeSize(udm, c, coneSize);
511: maxConeSize = PetscMax(maxConeSize, coneSize);
512: }
513: DMSetUp(udm);
514: PetscMalloc1(maxConeSize, &cone);
515: for (c = cStart; c < cEnd; ++c) {
516: PetscInt *closure = NULL, closureSize, cl, coneSize = 0;
518: DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
519: for (cl = 0; cl < closureSize*2; cl += 2) {
520: const PetscInt p = closure[cl];
522: if ((p >= vStart) && (p < vEnd)) cone[coneSize++] = p;
523: }
524: DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
525: DMPlexSetCone(udm, c, cone);
526: }
527: PetscFree(cone);
528: DMPlexSymmetrize(udm);
529: DMPlexStratify(udm);
530: *dmUnint = udm;
531: return(0);
532: }