Actual source code: plexinterpolate.c
petsc-3.4.5 2014-06-29
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: static PetscErrorCode DMPlexGetFaces_Internal(DM dm, PetscInt dim, PetscInt p, PetscInt *numFaces, PetscInt *faceSize, const PetscInt *faces[])
10: {
11: const PetscInt *cone = NULL;
12: PetscInt *facesTmp;
13: PetscInt maxConeSize, maxSupportSize, coneSize;
14: PetscErrorCode ierr;
18: DMPlexGetMaxSizes(dm, &maxConeSize, &maxSupportSize);
19: DMGetWorkArray(dm, PetscSqr(PetscMax(maxConeSize, maxSupportSize)), PETSC_INT, &facesTmp);
20: DMPlexGetConeSize(dm, p, &coneSize);
21: DMPlexGetCone(dm, p, &cone);
22: switch (dim) {
23: case 2:
24: switch (coneSize) {
25: case 3:
26: if (faces) {
27: facesTmp[0] = cone[0]; facesTmp[1] = cone[1];
28: facesTmp[2] = cone[1]; facesTmp[3] = cone[2];
29: facesTmp[4] = cone[2]; facesTmp[5] = cone[0];
30: *faces = facesTmp;
31: }
32: if (numFaces) *numFaces = 3;
33: if (faceSize) *faceSize = 2;
34: break;
35: case 4:
36: /* Vertices follow right hand rule */
37: if (faces) {
38: facesTmp[0] = cone[0]; facesTmp[1] = cone[1];
39: facesTmp[2] = cone[1]; facesTmp[3] = cone[2];
40: facesTmp[4] = cone[2]; facesTmp[5] = cone[3];
41: facesTmp[6] = cone[3]; facesTmp[7] = cone[0];
42: *faces = facesTmp;
43: }
44: if (numFaces) *numFaces = 4;
45: if (faceSize) *faceSize = 2;
46: if (faces) *faces = facesTmp;
47: break;
48: default:
49: SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Cone size %D not supported for dimension %D", coneSize, dim);
50: }
51: break;
52: case 3:
53: switch (coneSize) {
54: case 3:
55: if (faces) {
56: facesTmp[0] = cone[0]; facesTmp[1] = cone[1];
57: facesTmp[2] = cone[1]; facesTmp[3] = cone[2];
58: facesTmp[4] = cone[2]; facesTmp[5] = cone[0];
59: *faces = facesTmp;
60: }
61: if (numFaces) *numFaces = 3;
62: if (faceSize) *faceSize = 2;
63: if (faces) *faces = facesTmp;
64: break;
65: case 4:
66: /* Vertices of first face follow right hand rule and normal points away from last vertex */
67: if (faces) {
68: facesTmp[0] = cone[0]; facesTmp[1] = cone[1]; facesTmp[2] = cone[2];
69: facesTmp[3] = cone[0]; facesTmp[4] = cone[3]; facesTmp[5] = cone[1];
70: facesTmp[6] = cone[0]; facesTmp[7] = cone[2]; facesTmp[8] = cone[3];
71: facesTmp[9] = cone[2]; facesTmp[10] = cone[1]; facesTmp[11] = cone[3];
72: *faces = facesTmp;
73: }
74: if (numFaces) *numFaces = 4;
75: if (faceSize) *faceSize = 3;
76: if (faces) *faces = facesTmp;
77: break;
78: case 8:
79: if (faces) {
80: facesTmp[0] = cone[0]; facesTmp[1] = cone[1]; facesTmp[2] = cone[2]; facesTmp[3] = cone[3];
81: facesTmp[4] = cone[4]; facesTmp[5] = cone[5]; facesTmp[6] = cone[6]; facesTmp[7] = cone[7];
82: facesTmp[8] = cone[0]; facesTmp[9] = cone[3]; facesTmp[10] = cone[5]; facesTmp[11] = cone[4];
83: facesTmp[12] = cone[2]; facesTmp[13] = cone[1]; facesTmp[14] = cone[7]; facesTmp[15] = cone[6];
84: facesTmp[16] = cone[3]; facesTmp[17] = cone[2]; facesTmp[18] = cone[6]; facesTmp[19] = cone[5];
85: facesTmp[20] = cone[0]; facesTmp[21] = cone[4]; facesTmp[22] = cone[7]; facesTmp[23] = cone[1];
86: *faces = facesTmp;
87: }
88: if (numFaces) *numFaces = 6;
89: if (faceSize) *faceSize = 4;
90: break;
91: default:
92: SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Cone size %D not supported for dimension %D", coneSize, dim);
93: }
94: break;
95: default:
96: SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Dimension %D not supported", dim);
97: }
98: DMRestoreWorkArray(dm, 0, PETSC_INT, &facesTmp);
99: return(0);
100: }
104: /* This interpolates faces for cells at some stratum */
105: static PetscErrorCode DMPlexInterpolateFaces_Internal(DM dm, PetscInt cellDepth, DM idm)
106: {
107: DMLabel subpointMap;
108: PetscHashIJKL faceTable;
109: PetscInt *pStart, *pEnd;
110: PetscInt cellDim, depth, faceDepth = cellDepth, numPoints = 0, faceSizeAll = 0, face, c, d;
114: DMPlexGetDimension(dm, &cellDim);
115: /* HACK: I need a better way to determine face dimension, or an alternative to GetFaces() */
116: DMPlexGetSubpointMap(dm, &subpointMap);
117: if (subpointMap) ++cellDim;
118: DMPlexGetDepth(dm, &depth);
119: ++depth;
120: ++cellDepth;
121: cellDim -= depth - cellDepth;
122: PetscMalloc2(depth+1,PetscInt,&pStart,depth+1,PetscInt,&pEnd);
123: for (d = depth-1; d >= faceDepth; --d) {
124: DMPlexGetDepthStratum(dm, d, &pStart[d+1], &pEnd[d+1]);
125: }
126: DMPlexGetDepthStratum(dm, -1, NULL, &pStart[faceDepth]);
127: pEnd[faceDepth] = pStart[faceDepth];
128: for (d = faceDepth-1; d >= 0; --d) {
129: DMPlexGetDepthStratum(dm, d, &pStart[d], &pEnd[d]);
130: }
131: if (pEnd[cellDepth] > pStart[cellDepth]) {DMPlexGetFaces_Internal(dm, cellDim, pStart[cellDepth], NULL, &faceSizeAll, NULL);}
132: if (faceSizeAll > 4) SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "Do not support interpolation of meshes with faces of %D vertices", faceSizeAll);
133: PetscHashIJKLCreate(&faceTable);
134: PetscHashIJKLSetMultivalued(faceTable, PETSC_FALSE);
135: for (c = pStart[cellDepth], face = pStart[faceDepth]; c < pEnd[cellDepth]; ++c) {
136: const PetscInt *cellFaces;
137: PetscInt numCellFaces, faceSize, cf, f;
139: DMPlexGetFaces_Internal(dm, cellDim, c, &numCellFaces, &faceSize, &cellFaces);
140: if (faceSize != faceSizeAll) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Inconsistent face for cell %D of size %D != %D", c, faceSize, faceSizeAll);
141: for (cf = 0; cf < numCellFaces; ++cf) {
142: const PetscInt *cellFace = &cellFaces[cf*faceSize];
143: PetscHashIJKLKey key;
145: if (faceSize == 2) {
146: key.i = PetscMin(cellFace[0], cellFace[1]);
147: key.j = PetscMax(cellFace[0], cellFace[1]);
148: key.k = 0;
149: key.l = 0;
150: } else {
151: key.i = cellFace[0]; key.j = cellFace[1]; key.k = cellFace[2]; key.l = faceSize > 3 ? cellFace[3] : 0;
152: PetscSortInt(faceSize, (PetscInt *) &key);
153: }
154: PetscHashIJKLGet(faceTable, key, &f);
155: if (f < 0) {
156: PetscHashIJKLAdd(faceTable, key, face);
157: f = face++;
158: }
159: }
160: }
161: pEnd[faceDepth] = face;
162: PetscHashIJKLDestroy(&faceTable);
163: /* Count new points */
164: for (d = 0; d <= depth; ++d) {
165: numPoints += pEnd[d]-pStart[d];
166: }
167: DMPlexSetChart(idm, 0, numPoints);
168: /* Set cone sizes */
169: for (d = 0; d <= depth; ++d) {
170: PetscInt coneSize, p;
172: if (d == faceDepth) {
173: for (p = pStart[d]; p < pEnd[d]; ++p) {
174: /* I see no way to do this if we admit faces of different shapes */
175: DMPlexSetConeSize(idm, p, faceSizeAll);
176: }
177: } else if (d == cellDepth) {
178: for (p = pStart[d]; p < pEnd[d]; ++p) {
179: /* Number of cell faces may be different from number of cell vertices*/
180: DMPlexGetFaces_Internal(dm, cellDim, p, &coneSize, NULL, NULL);
181: DMPlexSetConeSize(idm, p, coneSize);
182: }
183: } else {
184: for (p = pStart[d]; p < pEnd[d]; ++p) {
185: DMPlexGetConeSize(dm, p, &coneSize);
186: DMPlexSetConeSize(idm, p, coneSize);
187: }
188: }
189: }
190: DMSetUp(idm);
191: /* Get face cones from subsets of cell vertices */
192: if (faceSizeAll > 4) SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "Do not support interpolation of meshes with faces of %D vertices", faceSizeAll);
193: PetscHashIJKLCreate(&faceTable);
194: PetscHashIJKLSetMultivalued(faceTable, PETSC_FALSE);
195: for (d = depth; d > cellDepth; --d) {
196: const PetscInt *cone;
197: PetscInt p;
199: for (p = pStart[d]; p < pEnd[d]; ++p) {
200: DMPlexGetCone(dm, p, &cone);
201: DMPlexSetCone(idm, p, cone);
202: DMPlexGetConeOrientation(dm, p, &cone);
203: DMPlexSetConeOrientation(idm, p, cone);
204: }
205: }
206: for (c = pStart[cellDepth], face = pStart[faceDepth]; c < pEnd[cellDepth]; ++c) {
207: const PetscInt *cellFaces;
208: PetscInt numCellFaces, faceSize, cf, f;
210: DMPlexGetFaces_Internal(dm, cellDim, c, &numCellFaces, &faceSize, &cellFaces);
211: if (faceSize != faceSizeAll) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Inconsistent face for cell %D of size %D != %D", c, faceSize, faceSizeAll);
212: for (cf = 0; cf < numCellFaces; ++cf) {
213: const PetscInt *cellFace = &cellFaces[cf*faceSize];
214: PetscHashIJKLKey key;
216: if (faceSize == 2) {
217: key.i = PetscMin(cellFace[0], cellFace[1]);
218: key.j = PetscMax(cellFace[0], cellFace[1]);
219: key.k = 0;
220: key.l = 0;
221: } else {
222: key.i = cellFace[0]; key.j = cellFace[1]; key.k = cellFace[2]; key.l = faceSize > 3 ? cellFace[3] : 0;
223: PetscSortInt(faceSize, (PetscInt *) &key);
224: }
225: PetscHashIJKLGet(faceTable, key, &f);
226: if (f < 0) {
227: DMPlexSetCone(idm, face, cellFace);
228: PetscHashIJKLAdd(faceTable, key, face);
229: f = face++;
230: DMPlexInsertCone(idm, c, cf, f);
231: } else {
232: const PetscInt *cone;
233: PetscInt coneSize, ornt, i, j;
235: DMPlexInsertCone(idm, c, cf, f);
236: /* Orient face: Do not allow reverse orientation at the first vertex */
237: DMPlexGetConeSize(idm, f, &coneSize);
238: DMPlexGetCone(idm, f, &cone);
239: 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);
240: /* - First find the initial vertex */
241: for (i = 0; i < faceSize; ++i) if (cellFace[0] == cone[i]) break;
242: /* - Try forward comparison */
243: for (j = 0; j < faceSize; ++j) if (cellFace[j] != cone[(i+j)%faceSize]) break;
244: if (j == faceSize) {
245: if ((faceSize == 2) && (i == 1)) ornt = -2;
246: else ornt = i;
247: } else {
248: /* - Try backward comparison */
249: for (j = 0; j < faceSize; ++j) if (cellFace[j] != cone[(i+faceSize-j)%faceSize]) break;
250: if (j == faceSize) {
251: if (i == 0) ornt = -faceSize;
252: else ornt = -(i+1);
253: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Could not determine face orientation");
254: }
255: DMPlexInsertConeOrientation(idm, c, cf, ornt);
256: }
257: }
258: }
259: 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]);
260: PetscFree2(pStart,pEnd);
261: PetscHashIJKLDestroy(&faceTable);
262: PetscFree2(pStart,pEnd);
263: DMPlexSymmetrize(idm);
264: DMPlexStratify(idm);
265: return(0);
266: }
270: /*@
271: DMPlexInterpolate - Take in a cell-vertex mesh and return one with all intermediate faces, edges, etc.
273: Collective on DM
275: Input Parameter:
276: . dmA - The DMPlex object with only cells and vertices
278: Output Parameter:
279: . dmB - The complete DMPlex object
281: Level: intermediate
283: .keywords: mesh
284: .seealso: DMPlexCreateFromCellList()
285: @*/
286: PetscErrorCode DMPlexInterpolate(DM dm, DM *dmInt)
287: {
288: DM idm, odm = dm;
289: PetscInt depth, dim, d;
293: DMPlexGetDepth(dm, &depth);
294: DMPlexGetDimension(dm, &dim);
295: if (dim <= 1) {
296: PetscObjectReference((PetscObject) dm);
297: idm = dm;
298: }
299: for (d = 1; d < dim; ++d) {
300: /* Create interpolated mesh */
301: DMCreate(PetscObjectComm((PetscObject)dm), &idm);
302: DMSetType(idm, DMPLEX);
303: DMPlexSetDimension(idm, dim);
304: if (depth > 0) {DMPlexInterpolateFaces_Internal(odm, 1, idm);}
305: if (odm != dm) {DMDestroy(&odm);}
306: odm = idm;
307: }
308: *dmInt = idm;
309: return(0);
310: }
314: /*@
315: DMPlexCopyCoordinates - Copy coordinates from one mesh to another with the same vertices
317: Collective on DM
319: Input Parameter:
320: . dmA - The DMPlex object with initial coordinates
322: Output Parameter:
323: . dmB - The DMPlex object with copied coordinates
325: Level: intermediate
327: Note: This is typically used when adding pieces other than vertices to a mesh
329: .keywords: mesh
330: .seealso: DMGetCoordinates(), DMGetCoordinatesLocal(), DMGetCoordinateDM(), DMPlexGetCoordinateSection()
331: @*/
332: PetscErrorCode DMPlexCopyCoordinates(DM dmA, DM dmB)
333: {
334: Vec coordinatesA, coordinatesB;
335: PetscSection coordSectionA, coordSectionB;
336: PetscScalar *coordsA, *coordsB;
337: PetscInt spaceDim, vStartA, vStartB, vEndA, vEndB, coordSizeB, v, d;
341: if (dmA == dmB) return(0);
342: DMPlexGetDepthStratum(dmA, 0, &vStartA, &vEndA);
343: DMPlexGetDepthStratum(dmB, 0, &vStartB, &vEndB);
344: 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);
345: DMPlexGetCoordinateSection(dmA, &coordSectionA);
346: DMPlexGetCoordinateSection(dmB, &coordSectionB);
347: PetscSectionSetNumFields(coordSectionB, 1);
348: PetscSectionGetFieldComponents(coordSectionA, 0, &spaceDim);
349: PetscSectionSetFieldComponents(coordSectionB, 0, spaceDim);
350: PetscSectionSetChart(coordSectionB, vStartB, vEndB);
351: for (v = vStartB; v < vEndB; ++v) {
352: PetscSectionSetDof(coordSectionB, v, spaceDim);
353: PetscSectionSetFieldDof(coordSectionB, v, 0, spaceDim);
354: }
355: PetscSectionSetUp(coordSectionB);
356: PetscSectionGetStorageSize(coordSectionB, &coordSizeB);
357: DMGetCoordinatesLocal(dmA, &coordinatesA);
358: VecCreate(PetscObjectComm((PetscObject) dmB), &coordinatesB);
359: PetscObjectSetName((PetscObject) coordinatesB, "coordinates");
360: VecSetSizes(coordinatesB, coordSizeB, PETSC_DETERMINE);
361: VecSetFromOptions(coordinatesB);
362: VecGetArray(coordinatesA, &coordsA);
363: VecGetArray(coordinatesB, &coordsB);
364: for (v = 0; v < vEndB-vStartB; ++v) {
365: for (d = 0; d < spaceDim; ++d) {
366: coordsB[v*spaceDim+d] = coordsA[v*spaceDim+d];
367: }
368: }
369: VecRestoreArray(coordinatesA, &coordsA);
370: VecRestoreArray(coordinatesB, &coordsB);
371: DMSetCoordinatesLocal(dmB, coordinatesB);
372: VecDestroy(&coordinatesB);
373: return(0);
374: }