1 #include <petsc-private/dmpleximpl.h> /*I "petscdmplex.h" I*/ 2 #include <../src/sys/utils/hash.h> 3 4 #undef __FUNCT__ 5 #define __FUNCT__ "DMPlexGetFaces_Internal" 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 *facesTmp; 13 PetscInt maxConeSize, maxSupportSize, coneSize; 14 PetscErrorCode ierr; 15 16 PetscFunctionBegin; 17 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 18 ierr = DMPlexGetMaxSizes(dm, &maxConeSize, &maxSupportSize);CHKERRQ(ierr); 19 ierr = DMGetWorkArray(dm, PetscSqr(PetscMax(maxConeSize, maxSupportSize)), PETSC_INT, &facesTmp);CHKERRQ(ierr); 20 ierr = DMPlexGetConeSize(dm, p, &coneSize);CHKERRQ(ierr); 21 ierr = DMPlexGetCone(dm, p, &cone);CHKERRQ(ierr); 22 ierr = DMPlexGetRawFaces_Internal(dm, dim, coneSize, cone, numFaces, faceSize, faces);CHKERRQ(ierr); 23 PetscFunctionReturn(0); 24 } 25 26 #undef __FUNCT__ 27 #define __FUNCT__ "DMPlexRestoreFaces_Internal" 28 /* 29 DMPlexRestoreFaces_Internal - Restores the array 30 */ 31 PetscErrorCode DMPlexRestoreFaces_Internal(DM dm, PetscInt dim, PetscInt p, PetscInt *numFaces, PetscInt *faceSize, const PetscInt *faces[]) 32 { 33 PetscErrorCode ierr; 34 35 PetscFunctionBegin; 36 ierr = DMRestoreWorkArray(dm, 0, PETSC_INT, faces);CHKERRQ(ierr); 37 PetscFunctionReturn(0); 38 } 39 40 #undef __FUNCT__ 41 #define __FUNCT__ "DMPlexGetRawFaces_Internal" 42 /* 43 DMPlexGetRawFaces_Internal - Gets groups of vertices that correspond to faces for the given cone 44 */ 45 PetscErrorCode DMPlexGetRawFaces_Internal(DM dm, PetscInt dim, PetscInt coneSize, const PetscInt cone[], PetscInt *numFaces, PetscInt *faceSize, const PetscInt *faces[]) 46 { 47 PetscInt *facesTmp; 48 PetscInt maxConeSize, maxSupportSize; 49 PetscErrorCode ierr; 50 51 PetscFunctionBegin; 52 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 53 ierr = DMPlexGetMaxSizes(dm, &maxConeSize, &maxSupportSize);CHKERRQ(ierr); 54 ierr = DMGetWorkArray(dm, PetscSqr(PetscMax(maxConeSize, maxSupportSize)), PETSC_INT, &facesTmp);CHKERRQ(ierr); 55 switch (dim) { 56 case 2: 57 switch (coneSize) { 58 case 3: 59 if (faces) { 60 facesTmp[0] = cone[0]; facesTmp[1] = cone[1]; 61 facesTmp[2] = cone[1]; facesTmp[3] = cone[2]; 62 facesTmp[4] = cone[2]; facesTmp[5] = cone[0]; 63 *faces = facesTmp; 64 } 65 if (numFaces) *numFaces = 3; 66 if (faceSize) *faceSize = 2; 67 break; 68 case 4: 69 /* Vertices follow right hand rule */ 70 if (faces) { 71 facesTmp[0] = cone[0]; facesTmp[1] = cone[1]; 72 facesTmp[2] = cone[1]; facesTmp[3] = cone[2]; 73 facesTmp[4] = cone[2]; facesTmp[5] = cone[3]; 74 facesTmp[6] = cone[3]; facesTmp[7] = cone[0]; 75 *faces = facesTmp; 76 } 77 if (numFaces) *numFaces = 4; 78 if (faceSize) *faceSize = 2; 79 if (faces) *faces = facesTmp; 80 break; 81 default: 82 SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Cone size %D not supported for dimension %D", coneSize, dim); 83 } 84 break; 85 case 3: 86 switch (coneSize) { 87 case 3: 88 if (faces) { 89 facesTmp[0] = cone[0]; facesTmp[1] = cone[1]; 90 facesTmp[2] = cone[1]; facesTmp[3] = cone[2]; 91 facesTmp[4] = cone[2]; facesTmp[5] = cone[0]; 92 *faces = facesTmp; 93 } 94 if (numFaces) *numFaces = 3; 95 if (faceSize) *faceSize = 2; 96 if (faces) *faces = facesTmp; 97 break; 98 case 4: 99 /* Vertices of first face follow right hand rule and normal points away from last vertex */ 100 if (faces) { 101 facesTmp[0] = cone[0]; facesTmp[1] = cone[1]; facesTmp[2] = cone[2]; 102 facesTmp[3] = cone[0]; facesTmp[4] = cone[3]; facesTmp[5] = cone[1]; 103 facesTmp[6] = cone[0]; facesTmp[7] = cone[2]; facesTmp[8] = cone[3]; 104 facesTmp[9] = cone[2]; facesTmp[10] = cone[1]; facesTmp[11] = cone[3]; 105 *faces = facesTmp; 106 } 107 if (numFaces) *numFaces = 4; 108 if (faceSize) *faceSize = 3; 109 if (faces) *faces = facesTmp; 110 break; 111 case 8: 112 if (faces) { 113 facesTmp[0] = cone[0]; facesTmp[1] = cone[1]; facesTmp[2] = cone[2]; facesTmp[3] = cone[3]; 114 facesTmp[4] = cone[4]; facesTmp[5] = cone[5]; facesTmp[6] = cone[6]; facesTmp[7] = cone[7]; 115 facesTmp[8] = cone[0]; facesTmp[9] = cone[3]; facesTmp[10] = cone[5]; facesTmp[11] = cone[4]; 116 facesTmp[12] = cone[2]; facesTmp[13] = cone[1]; facesTmp[14] = cone[7]; facesTmp[15] = cone[6]; 117 facesTmp[16] = cone[3]; facesTmp[17] = cone[2]; facesTmp[18] = cone[6]; facesTmp[19] = cone[5]; 118 facesTmp[20] = cone[0]; facesTmp[21] = cone[4]; facesTmp[22] = cone[7]; facesTmp[23] = cone[1]; 119 *faces = facesTmp; 120 } 121 if (numFaces) *numFaces = 6; 122 if (faceSize) *faceSize = 4; 123 break; 124 default: 125 SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Cone size %D not supported for dimension %D", coneSize, dim); 126 } 127 break; 128 default: 129 SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Dimension %D not supported", dim); 130 } 131 PetscFunctionReturn(0); 132 } 133 134 #undef __FUNCT__ 135 #define __FUNCT__ "DMPlexInterpolateFaces_Internal" 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; 143 PetscErrorCode ierr; 144 145 PetscFunctionBegin; 146 ierr = DMPlexGetDimension(dm, &cellDim);CHKERRQ(ierr); 147 /* HACK: I need a better way to determine face dimension, or an alternative to GetFaces() */ 148 ierr = DMPlexGetSubpointMap(dm, &subpointMap);CHKERRQ(ierr); 149 if (subpointMap) ++cellDim; 150 ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr); 151 ++depth; 152 ++cellDepth; 153 cellDim -= depth - cellDepth; 154 ierr = PetscMalloc2(depth+1,PetscInt,&pStart,depth+1,PetscInt,&pEnd);CHKERRQ(ierr); 155 for (d = depth-1; d >= faceDepth; --d) { 156 ierr = DMPlexGetDepthStratum(dm, d, &pStart[d+1], &pEnd[d+1]);CHKERRQ(ierr); 157 } 158 ierr = DMPlexGetDepthStratum(dm, -1, NULL, &pStart[faceDepth]);CHKERRQ(ierr); 159 pEnd[faceDepth] = pStart[faceDepth]; 160 for (d = faceDepth-1; d >= 0; --d) { 161 ierr = DMPlexGetDepthStratum(dm, d, &pStart[d], &pEnd[d]);CHKERRQ(ierr); 162 } 163 if (pEnd[cellDepth] > pStart[cellDepth]) {ierr = DMPlexGetFaces_Internal(dm, cellDim, pStart[cellDepth], NULL, &faceSizeAll, NULL);CHKERRQ(ierr);} 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 ierr = PetscHashIJKLCreate(&faceTable);CHKERRQ(ierr); 166 ierr = PetscHashIJKLSetMultivalued(faceTable, PETSC_FALSE);CHKERRQ(ierr); 167 for (c = pStart[cellDepth], face = pStart[faceDepth]; c < pEnd[cellDepth]; ++c) { 168 const PetscInt *cellFaces; 169 PetscInt numCellFaces, faceSize, cf, f; 170 171 ierr = DMPlexGetFaces_Internal(dm, cellDim, c, &numCellFaces, &faceSize, &cellFaces);CHKERRQ(ierr); 172 if (faceSize != faceSizeAll) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Inconsistent face for cell %D of size %D != %D", c, faceSize, faceSizeAll); 173 for (cf = 0; cf < numCellFaces; ++cf) { 174 const PetscInt *cellFace = &cellFaces[cf*faceSize]; 175 PetscHashIJKLKey key; 176 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 ierr = PetscSortInt(faceSize, (PetscInt *) &key); 185 } 186 ierr = PetscHashIJKLGet(faceTable, key, &f);CHKERRQ(ierr); 187 if (f < 0) { 188 ierr = PetscHashIJKLAdd(faceTable, key, face);CHKERRQ(ierr); 189 f = face++; 190 } 191 } 192 ierr = DMPlexRestoreFaces_Internal(dm, cellDim, c, &numCellFaces, &faceSize, &cellFaces);CHKERRQ(ierr); 193 } 194 pEnd[faceDepth] = face; 195 ierr = PetscHashIJKLDestroy(&faceTable);CHKERRQ(ierr); 196 /* Count new points */ 197 for (d = 0; d <= depth; ++d) { 198 numPoints += pEnd[d]-pStart[d]; 199 } 200 ierr = DMPlexSetChart(idm, 0, numPoints);CHKERRQ(ierr); 201 /* Set cone sizes */ 202 for (d = 0; d <= depth; ++d) { 203 PetscInt coneSize, p; 204 205 if (d == faceDepth) { 206 for (p = pStart[d]; p < pEnd[d]; ++p) { 207 /* I see no way to do this if we admit faces of different shapes */ 208 ierr = DMPlexSetConeSize(idm, p, faceSizeAll);CHKERRQ(ierr); 209 } 210 } else if (d == cellDepth) { 211 for (p = pStart[d]; p < pEnd[d]; ++p) { 212 /* Number of cell faces may be different from number of cell vertices*/ 213 ierr = DMPlexGetFaces_Internal(dm, cellDim, p, &coneSize, NULL, NULL);CHKERRQ(ierr); 214 ierr = DMPlexSetConeSize(idm, p, coneSize);CHKERRQ(ierr); 215 } 216 } else { 217 for (p = pStart[d]; p < pEnd[d]; ++p) { 218 ierr = DMPlexGetConeSize(dm, p, &coneSize);CHKERRQ(ierr); 219 ierr = DMPlexSetConeSize(idm, p, coneSize);CHKERRQ(ierr); 220 } 221 } 222 } 223 ierr = DMSetUp(idm);CHKERRQ(ierr); 224 /* Get face cones from subsets of cell vertices */ 225 if (faceSizeAll > 4) SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "Do not support interpolation of meshes with faces of %D vertices", faceSizeAll); 226 ierr = PetscHashIJKLCreate(&faceTable);CHKERRQ(ierr); 227 ierr = PetscHashIJKLSetMultivalued(faceTable, PETSC_FALSE);CHKERRQ(ierr); 228 for (d = depth; d > cellDepth; --d) { 229 const PetscInt *cone; 230 PetscInt p; 231 232 for (p = pStart[d]; p < pEnd[d]; ++p) { 233 ierr = DMPlexGetCone(dm, p, &cone);CHKERRQ(ierr); 234 ierr = DMPlexSetCone(idm, p, cone);CHKERRQ(ierr); 235 ierr = DMPlexGetConeOrientation(dm, p, &cone);CHKERRQ(ierr); 236 ierr = DMPlexSetConeOrientation(idm, p, cone);CHKERRQ(ierr); 237 } 238 } 239 for (c = pStart[cellDepth], face = pStart[faceDepth]; c < pEnd[cellDepth]; ++c) { 240 const PetscInt *cellFaces; 241 PetscInt numCellFaces, faceSize, cf, f; 242 243 ierr = DMPlexGetFaces_Internal(dm, cellDim, c, &numCellFaces, &faceSize, &cellFaces);CHKERRQ(ierr); 244 if (faceSize != faceSizeAll) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Inconsistent face for cell %D of size %D != %D", c, faceSize, faceSizeAll); 245 for (cf = 0; cf < numCellFaces; ++cf) { 246 const PetscInt *cellFace = &cellFaces[cf*faceSize]; 247 PetscHashIJKLKey key; 248 249 if (faceSize == 2) { 250 key.i = PetscMin(cellFace[0], cellFace[1]); 251 key.j = PetscMax(cellFace[0], cellFace[1]); 252 key.k = 0; 253 key.l = 0; 254 } else { 255 key.i = cellFace[0]; key.j = cellFace[1]; key.k = cellFace[2]; key.l = faceSize > 3 ? cellFace[3] : 0; 256 ierr = PetscSortInt(faceSize, (PetscInt *) &key); 257 } 258 ierr = PetscHashIJKLGet(faceTable, key, &f);CHKERRQ(ierr); 259 if (f < 0) { 260 ierr = DMPlexSetCone(idm, face, cellFace);CHKERRQ(ierr); 261 ierr = PetscHashIJKLAdd(faceTable, key, face);CHKERRQ(ierr); 262 f = face++; 263 ierr = DMPlexInsertCone(idm, c, cf, f);CHKERRQ(ierr); 264 } else { 265 const PetscInt *cone; 266 PetscInt coneSize, ornt, i, j; 267 268 ierr = DMPlexInsertCone(idm, c, cf, f);CHKERRQ(ierr); 269 /* Orient face: Do not allow reverse orientation at the first vertex */ 270 ierr = DMPlexGetConeSize(idm, f, &coneSize);CHKERRQ(ierr); 271 ierr = DMPlexGetCone(idm, f, &cone);CHKERRQ(ierr); 272 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); 273 /* - First find the initial vertex */ 274 for (i = 0; i < faceSize; ++i) if (cellFace[0] == cone[i]) break; 275 /* - Try forward comparison */ 276 for (j = 0; j < faceSize; ++j) if (cellFace[j] != cone[(i+j)%faceSize]) break; 277 if (j == faceSize) { 278 if ((faceSize == 2) && (i == 1)) ornt = -2; 279 else ornt = i; 280 } else { 281 /* - Try backward comparison */ 282 for (j = 0; j < faceSize; ++j) if (cellFace[j] != cone[(i+faceSize-j)%faceSize]) break; 283 if (j == faceSize) { 284 if (i == 0) ornt = -faceSize; 285 else ornt = -(i+1); 286 } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Could not determine face orientation"); 287 } 288 ierr = DMPlexInsertConeOrientation(idm, c, cf, ornt);CHKERRQ(ierr); 289 } 290 } 291 ierr = DMPlexRestoreFaces_Internal(dm, cellDim, c, &numCellFaces, &faceSize, &cellFaces);CHKERRQ(ierr); 292 } 293 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]); 294 ierr = PetscFree2(pStart,pEnd);CHKERRQ(ierr); 295 ierr = PetscHashIJKLDestroy(&faceTable);CHKERRQ(ierr); 296 ierr = PetscFree2(pStart,pEnd);CHKERRQ(ierr); 297 ierr = DMPlexSymmetrize(idm);CHKERRQ(ierr); 298 ierr = DMPlexStratify(idm);CHKERRQ(ierr); 299 PetscFunctionReturn(0); 300 } 301 302 #undef __FUNCT__ 303 #define __FUNCT__ "DMPlexInterpolate" 304 /*@ 305 DMPlexInterpolate - Take in a cell-vertex mesh and return one with all intermediate faces, edges, etc. 306 307 Collective on DM 308 309 Input Parameter: 310 . dm - The DMPlex object with only cells and vertices 311 312 Output Parameter: 313 . dmInt - The complete DMPlex object 314 315 Level: intermediate 316 317 .keywords: mesh 318 .seealso: DMPlexUninterpolate(), DMPlexCreateFromCellList() 319 @*/ 320 PetscErrorCode DMPlexInterpolate(DM dm, DM *dmInt) 321 { 322 DM idm, odm = dm; 323 PetscInt depth, dim, d; 324 PetscErrorCode ierr; 325 326 PetscFunctionBegin; 327 ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr); 328 ierr = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr); 329 if (dim <= 1) { 330 ierr = PetscObjectReference((PetscObject) dm);CHKERRQ(ierr); 331 idm = dm; 332 } 333 for (d = 1; d < dim; ++d) { 334 /* Create interpolated mesh */ 335 ierr = DMCreate(PetscObjectComm((PetscObject)dm), &idm);CHKERRQ(ierr); 336 ierr = DMSetType(idm, DMPLEX);CHKERRQ(ierr); 337 ierr = DMPlexSetDimension(idm, dim);CHKERRQ(ierr); 338 if (depth > 0) {ierr = DMPlexInterpolateFaces_Internal(odm, 1, idm);CHKERRQ(ierr);} 339 if (odm != dm) {ierr = DMDestroy(&odm);CHKERRQ(ierr);} 340 odm = idm; 341 } 342 *dmInt = idm; 343 PetscFunctionReturn(0); 344 } 345 346 #undef __FUNCT__ 347 #define __FUNCT__ "DMPlexCopyCoordinates" 348 /*@ 349 DMPlexCopyCoordinates - Copy coordinates from one mesh to another with the same vertices 350 351 Collective on DM 352 353 Input Parameter: 354 . dmA - The DMPlex object with initial coordinates 355 356 Output Parameter: 357 . dmB - The DMPlex object with copied coordinates 358 359 Level: intermediate 360 361 Note: This is typically used when adding pieces other than vertices to a mesh 362 363 .keywords: mesh 364 .seealso: DMGetCoordinates(), DMGetCoordinatesLocal(), DMGetCoordinateDM(), DMPlexGetCoordinateSection() 365 @*/ 366 PetscErrorCode DMPlexCopyCoordinates(DM dmA, DM dmB) 367 { 368 Vec coordinatesA, coordinatesB; 369 PetscSection coordSectionA, coordSectionB; 370 PetscScalar *coordsA, *coordsB; 371 PetscInt spaceDim, vStartA, vStartB, vEndA, vEndB, coordSizeB, v, d; 372 PetscErrorCode ierr; 373 374 PetscFunctionBegin; 375 if (dmA == dmB) PetscFunctionReturn(0); 376 ierr = DMPlexGetDepthStratum(dmA, 0, &vStartA, &vEndA);CHKERRQ(ierr); 377 ierr = DMPlexGetDepthStratum(dmB, 0, &vStartB, &vEndB);CHKERRQ(ierr); 378 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); 379 ierr = DMPlexGetCoordinateSection(dmA, &coordSectionA);CHKERRQ(ierr); 380 ierr = DMPlexGetCoordinateSection(dmB, &coordSectionB);CHKERRQ(ierr); 381 ierr = PetscSectionSetNumFields(coordSectionB, 1);CHKERRQ(ierr); 382 ierr = PetscSectionGetFieldComponents(coordSectionA, 0, &spaceDim);CHKERRQ(ierr); 383 ierr = PetscSectionSetFieldComponents(coordSectionB, 0, spaceDim);CHKERRQ(ierr); 384 ierr = PetscSectionSetChart(coordSectionB, vStartB, vEndB);CHKERRQ(ierr); 385 for (v = vStartB; v < vEndB; ++v) { 386 ierr = PetscSectionSetDof(coordSectionB, v, spaceDim);CHKERRQ(ierr); 387 ierr = PetscSectionSetFieldDof(coordSectionB, v, 0, spaceDim);CHKERRQ(ierr); 388 } 389 ierr = PetscSectionSetUp(coordSectionB);CHKERRQ(ierr); 390 ierr = PetscSectionGetStorageSize(coordSectionB, &coordSizeB);CHKERRQ(ierr); 391 ierr = DMGetCoordinatesLocal(dmA, &coordinatesA);CHKERRQ(ierr); 392 ierr = VecCreate(PetscObjectComm((PetscObject) dmB), &coordinatesB);CHKERRQ(ierr); 393 ierr = PetscObjectSetName((PetscObject) coordinatesB, "coordinates");CHKERRQ(ierr); 394 ierr = VecSetSizes(coordinatesB, coordSizeB, PETSC_DETERMINE);CHKERRQ(ierr); 395 ierr = VecSetFromOptions(coordinatesB);CHKERRQ(ierr); 396 ierr = VecGetArray(coordinatesA, &coordsA);CHKERRQ(ierr); 397 ierr = VecGetArray(coordinatesB, &coordsB);CHKERRQ(ierr); 398 for (v = 0; v < vEndB-vStartB; ++v) { 399 for (d = 0; d < spaceDim; ++d) { 400 coordsB[v*spaceDim+d] = coordsA[v*spaceDim+d]; 401 } 402 } 403 ierr = VecRestoreArray(coordinatesA, &coordsA);CHKERRQ(ierr); 404 ierr = VecRestoreArray(coordinatesB, &coordsB);CHKERRQ(ierr); 405 ierr = DMSetCoordinatesLocal(dmB, coordinatesB);CHKERRQ(ierr); 406 ierr = VecDestroy(&coordinatesB);CHKERRQ(ierr); 407 PetscFunctionReturn(0); 408 } 409 410 #undef __FUNCT__ 411 #define __FUNCT__ "DMPlexUninterpolate" 412 /*@ 413 DMPlexUninterpolate - Take in a mesh with all intermediate faces, edges, etc. and return a cell-vertex mesh 414 415 Collective on DM 416 417 Input Parameter: 418 . dm - The complete DMPlex object 419 420 Output Parameter: 421 . dmUnint - The DMPlex object with only cells and vertices 422 423 Level: intermediate 424 425 .keywords: mesh 426 .seealso: DMPlexInterpolate(), DMPlexCreateFromCellList() 427 @*/ 428 PetscErrorCode DMPlexUninterpolate(DM dm, DM *dmUnint) 429 { 430 DM udm; 431 PetscInt dim, vStart, vEnd, cStart, cEnd, c, maxConeSize = 0, *cone; 432 PetscErrorCode ierr; 433 434 PetscFunctionBegin; 435 ierr = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr); 436 if (dim <= 1) { 437 ierr = PetscObjectReference((PetscObject) dm);CHKERRQ(ierr); 438 *dmUnint = dm; 439 PetscFunctionReturn(0); 440 } 441 ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr); 442 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 443 ierr = DMCreate(PetscObjectComm((PetscObject) dm), &udm);CHKERRQ(ierr); 444 ierr = DMSetType(udm, DMPLEX);CHKERRQ(ierr); 445 ierr = DMPlexSetDimension(udm, dim);CHKERRQ(ierr); 446 ierr = DMPlexSetChart(udm, cStart, vEnd);CHKERRQ(ierr); 447 for (c = cStart; c < cEnd; ++c) { 448 PetscInt *closure = NULL, closureSize, cl, coneSize = 0; 449 450 ierr = DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 451 for (cl = 0; cl < closureSize*2; cl += 2) { 452 const PetscInt p = closure[cl]; 453 454 if ((p >= vStart) && (p < vEnd)) ++coneSize; 455 } 456 ierr = DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 457 ierr = DMPlexSetConeSize(udm, c, coneSize);CHKERRQ(ierr); 458 maxConeSize = PetscMax(maxConeSize, coneSize); 459 } 460 ierr = DMSetUp(udm);CHKERRQ(ierr); 461 ierr = PetscMalloc(maxConeSize * sizeof(PetscInt), &cone);CHKERRQ(ierr); 462 for (c = cStart; c < cEnd; ++c) { 463 PetscInt *closure = NULL, closureSize, cl, coneSize = 0; 464 465 ierr = DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 466 for (cl = 0; cl < closureSize*2; cl += 2) { 467 const PetscInt p = closure[cl]; 468 469 if ((p >= vStart) && (p < vEnd)) cone[coneSize++] = p; 470 } 471 ierr = DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 472 ierr = DMPlexSetCone(udm, c, cone);CHKERRQ(ierr); 473 } 474 ierr = PetscFree(cone);CHKERRQ(ierr); 475 ierr = DMPlexSymmetrize(udm);CHKERRQ(ierr); 476 ierr = DMPlexStratify(udm);CHKERRQ(ierr); 477 *dmUnint = udm; 478 PetscFunctionReturn(0); 479 } 480