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