1 #include <petsc-private/pleximpl.h> /*I "petscdmplex.h" I*/ 2 3 extern PetscErrorCode DMPlexGetNumFaceVertices_Internal(DM, PetscInt, PetscInt *); 4 5 #undef __FUNCT__ 6 #define __FUNCT__ "DMPlexGetFaceOrientation" 7 PetscErrorCode DMPlexGetFaceOrientation(DM dm, PetscInt cell, PetscInt numCorners, PetscInt indices[], PetscInt oppositeVertex, PetscInt origVertices[], PetscInt faceVertices[], PetscBool *posOriented) 8 { 9 MPI_Comm comm = ((PetscObject) dm)->comm; 10 PetscBool posOrient = PETSC_FALSE; 11 const PetscInt debug = 0; 12 PetscInt cellDim, faceSize, f; 13 PetscErrorCode ierr; 14 15 ierr = DMPlexGetDimension(dm, &cellDim);CHKERRQ(ierr); 16 if (debug) {PetscPrintf(comm, "cellDim: %d numCorners: %d\n", cellDim, numCorners);CHKERRQ(ierr);} 17 18 if (cellDim == numCorners-1) { 19 /* Simplices */ 20 faceSize = numCorners-1; 21 posOrient = !(oppositeVertex%2) ? PETSC_TRUE : PETSC_FALSE; 22 } else if (cellDim == 1 && numCorners == 3) { 23 /* Quadratic line */ 24 faceSize = 1; 25 posOrient = PETSC_TRUE; 26 } else if (cellDim == 2 && numCorners == 4) { 27 /* Quads */ 28 faceSize = 2; 29 if ((indices[1] > indices[0]) && (indices[1] - indices[0] == 1)) { 30 posOrient = PETSC_TRUE; 31 } else if ((indices[0] == 3) && (indices[1] == 0)) { 32 posOrient = PETSC_TRUE; 33 } else { 34 if (((indices[0] > indices[1]) && (indices[0] - indices[1] == 1)) || ((indices[0] == 0) && (indices[1] == 3))) { 35 posOrient = PETSC_FALSE; 36 } else SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid quad crossedge"); 37 } 38 } else if (cellDim == 2 && numCorners == 6) { 39 /* Quadratic triangle (I hate this) */ 40 /* Edges are determined by the first 2 vertices (corners of edges) */ 41 const PetscInt faceSizeTri = 3; 42 PetscInt sortedIndices[3], i, iFace; 43 PetscBool found = PETSC_FALSE; 44 PetscInt faceVerticesTriSorted[9] = { 45 0, 3, 4, /* bottom */ 46 1, 4, 5, /* right */ 47 2, 3, 5, /* left */ 48 }; 49 PetscInt faceVerticesTri[9] = { 50 0, 3, 4, /* bottom */ 51 1, 4, 5, /* right */ 52 2, 5, 3, /* left */ 53 }; 54 55 faceSize = faceSizeTri; 56 for (i = 0; i < faceSizeTri; ++i) sortedIndices[i] = indices[i]; 57 ierr = PetscSortInt(faceSizeTri, sortedIndices);CHKERRQ(ierr); 58 for (iFace = 0; iFace < 3; ++iFace) { 59 const PetscInt ii = iFace*faceSizeTri; 60 PetscInt fVertex, cVertex; 61 62 if ((sortedIndices[0] == faceVerticesTriSorted[ii+0]) && 63 (sortedIndices[1] == faceVerticesTriSorted[ii+1])) { 64 for (fVertex = 0; fVertex < faceSizeTri; ++fVertex) { 65 for (cVertex = 0; cVertex < faceSizeTri; ++cVertex) { 66 if (indices[cVertex] == faceVerticesTri[ii+fVertex]) { 67 faceVertices[fVertex] = origVertices[cVertex]; 68 break; 69 } 70 } 71 } 72 found = PETSC_TRUE; 73 break; 74 } 75 } 76 if (!found) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid tri crossface"); 77 if (posOriented) *posOriented = PETSC_TRUE; 78 PetscFunctionReturn(0); 79 } else if (cellDim == 2 && numCorners == 9) { 80 /* Quadratic quad (I hate this) */ 81 /* Edges are determined by the first 2 vertices (corners of edges) */ 82 const PetscInt faceSizeQuad = 3; 83 PetscInt sortedIndices[3], i, iFace; 84 PetscBool found = PETSC_FALSE; 85 PetscInt faceVerticesQuadSorted[12] = { 86 0, 1, 4, /* bottom */ 87 1, 2, 5, /* right */ 88 2, 3, 6, /* top */ 89 0, 3, 7, /* left */ 90 }; 91 PetscInt faceVerticesQuad[12] = { 92 0, 1, 4, /* bottom */ 93 1, 2, 5, /* right */ 94 2, 3, 6, /* top */ 95 3, 0, 7, /* left */ 96 }; 97 98 faceSize = faceSizeQuad; 99 for (i = 0; i < faceSizeQuad; ++i) sortedIndices[i] = indices[i]; 100 ierr = PetscSortInt(faceSizeQuad, sortedIndices);CHKERRQ(ierr); 101 for (iFace = 0; iFace < 4; ++iFace) { 102 const PetscInt ii = iFace*faceSizeQuad; 103 PetscInt fVertex, cVertex; 104 105 if ((sortedIndices[0] == faceVerticesQuadSorted[ii+0]) && 106 (sortedIndices[1] == faceVerticesQuadSorted[ii+1])) { 107 for (fVertex = 0; fVertex < faceSizeQuad; ++fVertex) { 108 for (cVertex = 0; cVertex < faceSizeQuad; ++cVertex) { 109 if (indices[cVertex] == faceVerticesQuad[ii+fVertex]) { 110 faceVertices[fVertex] = origVertices[cVertex]; 111 break; 112 } 113 } 114 } 115 found = PETSC_TRUE; 116 break; 117 } 118 } 119 if (!found) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid quad crossface"); 120 if (posOriented) *posOriented = PETSC_TRUE; 121 PetscFunctionReturn(0); 122 } else if (cellDim == 3 && numCorners == 8) { 123 /* Hexes 124 A hex is two oriented quads with the normal of the first 125 pointing up at the second. 126 127 7---6 128 /| /| 129 4---5 | 130 | 3-|-2 131 |/ |/ 132 0---1 133 134 Faces are determined by the first 4 vertices (corners of faces) */ 135 const PetscInt faceSizeHex = 4; 136 PetscInt sortedIndices[4], i, iFace; 137 PetscBool found = PETSC_FALSE; 138 PetscInt faceVerticesHexSorted[24] = { 139 0, 1, 2, 3, /* bottom */ 140 4, 5, 6, 7, /* top */ 141 0, 1, 4, 5, /* front */ 142 1, 2, 5, 6, /* right */ 143 2, 3, 6, 7, /* back */ 144 0, 3, 4, 7, /* left */ 145 }; 146 PetscInt faceVerticesHex[24] = { 147 3, 2, 1, 0, /* bottom */ 148 4, 5, 6, 7, /* top */ 149 0, 1, 5, 4, /* front */ 150 1, 2, 6, 5, /* right */ 151 2, 3, 7, 6, /* back */ 152 3, 0, 4, 7, /* left */ 153 }; 154 155 faceSize = faceSizeHex; 156 for (i = 0; i < faceSizeHex; ++i) sortedIndices[i] = indices[i]; 157 ierr = PetscSortInt(faceSizeHex, sortedIndices);CHKERRQ(ierr); 158 for (iFace = 0; iFace < 6; ++iFace) { 159 const PetscInt ii = iFace*faceSizeHex; 160 PetscInt fVertex, cVertex; 161 162 if ((sortedIndices[0] == faceVerticesHexSorted[ii+0]) && 163 (sortedIndices[1] == faceVerticesHexSorted[ii+1]) && 164 (sortedIndices[2] == faceVerticesHexSorted[ii+2]) && 165 (sortedIndices[3] == faceVerticesHexSorted[ii+3])) { 166 for (fVertex = 0; fVertex < faceSizeHex; ++fVertex) { 167 for (cVertex = 0; cVertex < faceSizeHex; ++cVertex) { 168 if (indices[cVertex] == faceVerticesHex[ii+fVertex]) { 169 faceVertices[fVertex] = origVertices[cVertex]; 170 break; 171 } 172 } 173 } 174 found = PETSC_TRUE; 175 break; 176 } 177 } 178 if (!found) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid hex crossface"); 179 if (posOriented) *posOriented = PETSC_TRUE; 180 PetscFunctionReturn(0); 181 } else if (cellDim == 3 && numCorners == 10) { 182 /* Quadratic tet */ 183 /* Faces are determined by the first 3 vertices (corners of faces) */ 184 const PetscInt faceSizeTet = 6; 185 PetscInt sortedIndices[6], i, iFace; 186 PetscBool found = PETSC_FALSE; 187 PetscInt faceVerticesTetSorted[24] = { 188 0, 1, 2, 6, 7, 8, /* bottom */ 189 0, 3, 4, 6, 7, 9, /* front */ 190 1, 4, 5, 7, 8, 9, /* right */ 191 2, 3, 5, 6, 8, 9, /* left */ 192 }; 193 PetscInt faceVerticesTet[24] = { 194 0, 1, 2, 6, 7, 8, /* bottom */ 195 0, 4, 3, 6, 7, 9, /* front */ 196 1, 5, 4, 7, 8, 9, /* right */ 197 2, 3, 5, 8, 6, 9, /* left */ 198 }; 199 200 faceSize = faceSizeTet; 201 for (i = 0; i < faceSizeTet; ++i) sortedIndices[i] = indices[i]; 202 ierr = PetscSortInt(faceSizeTet, sortedIndices);CHKERRQ(ierr); 203 for (iFace=0; iFace < 4; ++iFace) { 204 const PetscInt ii = iFace*faceSizeTet; 205 PetscInt fVertex, cVertex; 206 207 if ((sortedIndices[0] == faceVerticesTetSorted[ii+0]) && 208 (sortedIndices[1] == faceVerticesTetSorted[ii+1]) && 209 (sortedIndices[2] == faceVerticesTetSorted[ii+2]) && 210 (sortedIndices[3] == faceVerticesTetSorted[ii+3])) { 211 for (fVertex = 0; fVertex < faceSizeTet; ++fVertex) { 212 for (cVertex = 0; cVertex < faceSizeTet; ++cVertex) { 213 if (indices[cVertex] == faceVerticesTet[ii+fVertex]) { 214 faceVertices[fVertex] = origVertices[cVertex]; 215 break; 216 } 217 } 218 } 219 found = PETSC_TRUE; 220 break; 221 } 222 } 223 if (!found) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid tet crossface"); 224 if (posOriented) *posOriented = PETSC_TRUE; 225 PetscFunctionReturn(0); 226 } else if (cellDim == 3 && numCorners == 27) { 227 /* Quadratic hexes (I hate this) 228 A hex is two oriented quads with the normal of the first 229 pointing up at the second. 230 231 7---6 232 /| /| 233 4---5 | 234 | 3-|-2 235 |/ |/ 236 0---1 237 238 Faces are determined by the first 4 vertices (corners of faces) */ 239 const PetscInt faceSizeQuadHex = 9; 240 PetscInt sortedIndices[9], i, iFace; 241 PetscBool found = PETSC_FALSE; 242 PetscInt faceVerticesQuadHexSorted[54] = { 243 0, 1, 2, 3, 8, 9, 10, 11, 24, /* bottom */ 244 4, 5, 6, 7, 12, 13, 14, 15, 25, /* top */ 245 0, 1, 4, 5, 8, 12, 16, 17, 22, /* front */ 246 1, 2, 5, 6, 9, 13, 17, 18, 21, /* right */ 247 2, 3, 6, 7, 10, 14, 18, 19, 23, /* back */ 248 0, 3, 4, 7, 11, 15, 16, 19, 20, /* left */ 249 }; 250 PetscInt faceVerticesQuadHex[54] = { 251 3, 2, 1, 0, 10, 9, 8, 11, 24, /* bottom */ 252 4, 5, 6, 7, 12, 13, 14, 15, 25, /* top */ 253 0, 1, 5, 4, 8, 17, 12, 16, 22, /* front */ 254 1, 2, 6, 5, 9, 18, 13, 17, 21, /* right */ 255 2, 3, 7, 6, 10, 19, 14, 18, 23, /* back */ 256 3, 0, 4, 7, 11, 16, 15, 19, 20 /* left */ 257 }; 258 259 faceSize = faceSizeQuadHex; 260 for (i = 0; i < faceSizeQuadHex; ++i) sortedIndices[i] = indices[i]; 261 ierr = PetscSortInt(faceSizeQuadHex, sortedIndices);CHKERRQ(ierr); 262 for (iFace = 0; iFace < 6; ++iFace) { 263 const PetscInt ii = iFace*faceSizeQuadHex; 264 PetscInt fVertex, cVertex; 265 266 if ((sortedIndices[0] == faceVerticesQuadHexSorted[ii+0]) && 267 (sortedIndices[1] == faceVerticesQuadHexSorted[ii+1]) && 268 (sortedIndices[2] == faceVerticesQuadHexSorted[ii+2]) && 269 (sortedIndices[3] == faceVerticesQuadHexSorted[ii+3])) { 270 for (fVertex = 0; fVertex < faceSizeQuadHex; ++fVertex) { 271 for (cVertex = 0; cVertex < faceSizeQuadHex; ++cVertex) { 272 if (indices[cVertex] == faceVerticesQuadHex[ii+fVertex]) { 273 faceVertices[fVertex] = origVertices[cVertex]; 274 break; 275 } 276 } 277 } 278 found = PETSC_TRUE; 279 break; 280 } 281 } 282 if (!found) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid hex crossface"); 283 if (posOriented) *posOriented = PETSC_TRUE; 284 PetscFunctionReturn(0); 285 } else SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Unknown cell type for faceOrientation()."); 286 if (!posOrient) { 287 if (debug) {ierr = PetscPrintf(comm, " Reversing initial face orientation\n");CHKERRQ(ierr);} 288 for (f = 0; f < faceSize; ++f) faceVertices[f] = origVertices[faceSize-1 - f]; 289 } else { 290 if (debug) {ierr = PetscPrintf(comm, " Keeping initial face orientation\n");CHKERRQ(ierr);} 291 for (f = 0; f < faceSize; ++f) faceVertices[f] = origVertices[f]; 292 } 293 if (posOriented) *posOriented = posOrient; 294 PetscFunctionReturn(0); 295 } 296 297 #undef __FUNCT__ 298 #define __FUNCT__ "DMPlexGetOrientedFace" 299 /* 300 Given a cell and a face, as a set of vertices, 301 return the oriented face, as a set of vertices, in faceVertices 302 The orientation is such that the face normal points out of the cell 303 */ 304 PetscErrorCode DMPlexGetOrientedFace(DM dm, PetscInt cell, PetscInt faceSize, const PetscInt face[], PetscInt numCorners, PetscInt indices[], PetscInt origVertices[], PetscInt faceVertices[], PetscBool *posOriented) 305 { 306 const PetscInt *cone = PETSC_NULL; 307 PetscInt coneSize, v, f, v2; 308 PetscInt oppositeVertex = -1; 309 PetscErrorCode ierr; 310 311 PetscFunctionBegin; 312 ierr = DMPlexGetConeSize(dm, cell, &coneSize);CHKERRQ(ierr); 313 ierr = DMPlexGetCone(dm, cell, &cone);CHKERRQ(ierr); 314 for (v = 0, v2 = 0; v < coneSize; ++v) { 315 PetscBool found = PETSC_FALSE; 316 317 for (f = 0; f < faceSize; ++f) { 318 if (face[f] == cone[v]) { 319 found = PETSC_TRUE; break; 320 } 321 } 322 if (found) { 323 indices[v2] = v; 324 origVertices[v2] = cone[v]; 325 ++v2; 326 } else { 327 oppositeVertex = v; 328 } 329 } 330 ierr = DMPlexGetFaceOrientation(dm, cell, numCorners, indices, oppositeVertex, origVertices, faceVertices, posOriented);CHKERRQ(ierr); 331 PetscFunctionReturn(0); 332 } 333 334 #undef __FUNCT__ 335 #define __FUNCT__ "DMPlexInsertFace_Private" 336 /* 337 DMPlexInsertFace_Private - Puts a face into the mesh 338 339 Not collective 340 341 Input Parameters: 342 + dm - The DMPlex 343 . numFaceVertex - The number of vertices in the face 344 . faceVertices - The vertices in the face for dm 345 . subfaceVertices - The vertices in the face for subdm 346 . numCorners - The number of vertices in the cell 347 . cell - A cell in dm containing the face 348 . subcell - A cell in subdm containing the face 349 . firstFace - First face in the mesh 350 - newFacePoint - Next face in the mesh 351 352 Output Parameters: 353 . newFacePoint - Contains next face point number on input, updated on output 354 355 Level: developer 356 */ 357 PetscErrorCode DMPlexInsertFace_Private(DM dm, DM subdm, PetscInt numFaceVertices, const PetscInt faceVertices[], const PetscInt subfaceVertices[], PetscInt numCorners, PetscInt cell, PetscInt subcell, PetscInt firstFace, PetscInt *newFacePoint) 358 { 359 MPI_Comm comm = ((PetscObject) dm)->comm; 360 DM_Plex *submesh = (DM_Plex*) subdm->data; 361 const PetscInt *faces; 362 PetscInt numFaces, coneSize; 363 PetscErrorCode ierr; 364 365 PetscFunctionBegin; 366 ierr = DMPlexGetConeSize(subdm, subcell, &coneSize);CHKERRQ(ierr); 367 if (coneSize != 1) SETERRQ2(comm, PETSC_ERR_ARG_OUTOFRANGE, "Cone size of cell %d is %d != 1", cell, coneSize); 368 #if 0 369 /* Cannot use this because support() has not been constructed yet */ 370 ierr = DMPlexGetJoin(subdm, numFaceVertices, subfaceVertices, &numFaces, &faces);CHKERRQ(ierr); 371 #else 372 { 373 PetscInt f; 374 375 numFaces = 0; 376 ierr = DMGetWorkArray(subdm, 1, PETSC_INT, (void**) &faces);CHKERRQ(ierr); 377 for (f = firstFace; f < *newFacePoint; ++f) { 378 PetscInt dof, off, d; 379 380 ierr = PetscSectionGetDof(submesh->coneSection, f, &dof);CHKERRQ(ierr); 381 ierr = PetscSectionGetOffset(submesh->coneSection, f, &off);CHKERRQ(ierr); 382 /* Yes, I know this is quadratic, but I expect the sizes to be <5 */ 383 for (d = 0; d < dof; ++d) { 384 const PetscInt p = submesh->cones[off+d]; 385 PetscInt v; 386 387 for (v = 0; v < numFaceVertices; ++v) { 388 if (subfaceVertices[v] == p) break; 389 } 390 if (v == numFaceVertices) break; 391 } 392 if (d == dof) { 393 numFaces = 1; 394 ((PetscInt*) faces)[0] = f; 395 } 396 } 397 } 398 #endif 399 if (numFaces > 1) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Vertex set had %d faces, not one", numFaces); 400 else if (numFaces == 1) { 401 /* Add the other cell neighbor for this face */ 402 ierr = DMPlexSetCone(subdm, cell, faces);CHKERRQ(ierr); 403 } else { 404 PetscInt *indices, *origVertices, *orientedVertices, *orientedSubVertices, v, ov; 405 PetscBool posOriented; 406 407 ierr = DMGetWorkArray(subdm, 4*numFaceVertices * sizeof(PetscInt), PETSC_INT, &orientedVertices);CHKERRQ(ierr); 408 origVertices = &orientedVertices[numFaceVertices]; 409 indices = &orientedVertices[numFaceVertices*2]; 410 orientedSubVertices = &orientedVertices[numFaceVertices*3]; 411 ierr = DMPlexGetOrientedFace(dm, cell, numFaceVertices, faceVertices, numCorners, indices, origVertices, orientedVertices, &posOriented);CHKERRQ(ierr); 412 /* TODO: I know that routine should return a permutation, not the indices */ 413 for (v = 0; v < numFaceVertices; ++v) { 414 const PetscInt vertex = faceVertices[v], subvertex = subfaceVertices[v]; 415 for (ov = 0; ov < numFaceVertices; ++ov) { 416 if (orientedVertices[ov] == vertex) { 417 orientedSubVertices[ov] = subvertex; 418 break; 419 } 420 } 421 if (ov == numFaceVertices) SETERRQ1(comm, PETSC_ERR_PLIB, "Could not find face vertex %d in orientated set", vertex); 422 } 423 ierr = DMPlexSetCone(subdm, *newFacePoint, orientedSubVertices);CHKERRQ(ierr); 424 ierr = DMPlexSetCone(subdm, subcell, newFacePoint);CHKERRQ(ierr); 425 ierr = DMRestoreWorkArray(subdm, 4*numFaceVertices * sizeof(PetscInt), PETSC_INT, &orientedVertices);CHKERRQ(ierr); 426 ++(*newFacePoint); 427 } 428 ierr = DMPlexRestoreJoin(subdm, numFaceVertices, subfaceVertices, &numFaces, &faces);CHKERRQ(ierr); 429 PetscFunctionReturn(0); 430 } 431 432 #undef __FUNCT__ 433 #define __FUNCT__ "DMPlexCreateSubmesh_Uninterpolated" 434 static PetscErrorCode DMPlexCreateSubmesh_Uninterpolated(DM dm, const char label[], DM subdm) 435 { 436 MPI_Comm comm = ((PetscObject) dm)->comm; 437 PetscBool boundaryFaces = PETSC_FALSE; 438 PetscSection coordSection, subCoordSection; 439 Vec coordinates, subCoordinates; 440 PetscScalar *coords, *subCoords; 441 DMLabel subpointMap; 442 IS labelIS; 443 const PetscInt *subVertices; 444 PetscInt *subVerticesActive; 445 PetscInt *subCells = PETSC_NULL; 446 PetscInt numSubVertices, numSubVerticesActive, firstSubVertex, numSubCells = 0, maxSubCells = 0, numOldSubCells; 447 PetscInt *face, *subface, maxConeSize, numSubFaces = 0, firstSubFace, newFacePoint, nFV = 0, coordSize; 448 PetscInt dim; /* Right now, do not specify dimension */ 449 PetscInt cStart, cEnd, cMax, c, vStart, vEnd, vMax, v, p, corner, i, d, f; 450 PetscErrorCode ierr; 451 452 PetscFunctionBegin; 453 ierr = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr); 454 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 455 ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr); 456 ierr = DMPlexGetMaxSizes(dm, &maxConeSize, PETSC_NULL);CHKERRQ(ierr); 457 ierr = DMPlexGetHybridBounds(dm, &cMax, PETSC_NULL, PETSC_NULL, &vMax);CHKERRQ(ierr); 458 if (cMax >= 0) cEnd = PetscMin(cEnd, cMax); 459 if (vMax >= 0) vEnd = PetscMin(vEnd, vMax); 460 ierr = DMGetWorkArray(dm, 2*maxConeSize, PETSC_INT, &face);CHKERRQ(ierr); 461 subface = &face[maxConeSize]; 462 ierr = DMPlexGetStratumIS(dm, label, 1, &labelIS);CHKERRQ(ierr); 463 if (labelIS) { 464 ierr = ISGetSize(labelIS, &numSubVertices);CHKERRQ(ierr); 465 ierr = ISGetIndices(labelIS, &subVertices);CHKERRQ(ierr); 466 } 467 maxSubCells = numSubVertices; 468 ierr = PetscMalloc(maxSubCells * sizeof(PetscInt), &subCells);CHKERRQ(ierr); 469 ierr = PetscMalloc(numSubVertices * sizeof(PetscInt), &subVerticesActive);CHKERRQ(ierr); 470 ierr = PetscMemzero(subVerticesActive, numSubVertices * sizeof(PetscInt));CHKERRQ(ierr); 471 for (v = 0; v < numSubVertices; ++v) { 472 const PetscInt vertex = subVertices[v]; 473 PetscInt *star = PETSC_NULL; 474 PetscInt starSize, numCells = 0; 475 476 ierr = DMPlexGetTransitiveClosure(dm, vertex, PETSC_FALSE, &starSize, &star);CHKERRQ(ierr); 477 for (p = 0; p < starSize*2; p += 2) { 478 const PetscInt point = star[p]; 479 if ((point >= cStart) && (point < cEnd)) star[numCells++] = point; 480 } 481 numOldSubCells = numSubCells; 482 for (c = 0; c < numCells; ++c) { 483 const PetscInt cell = star[c]; 484 PetscInt *closure = PETSC_NULL; 485 PetscInt closureSize, numCorners = 0, faceSize = 0; 486 PetscInt cellLoc; 487 488 ierr = PetscFindInt(cell, numOldSubCells, subCells, &cellLoc);CHKERRQ(ierr); 489 if (cellLoc >= 0) continue; 490 ierr = DMPlexGetTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 491 for (p = 0; p < closureSize*2; p += 2) { 492 const PetscInt point = closure[p]; 493 if ((point >= vStart) && (point < vEnd)) closure[numCorners++] = point; 494 } 495 if (!nFV) {ierr = DMPlexGetNumFaceVertices_Internal(dm, numCorners, &nFV);CHKERRQ(ierr);} 496 for (corner = 0; corner < numCorners; ++corner) { 497 const PetscInt cellVertex = closure[corner]; 498 PetscInt subVertex; 499 500 ierr = PetscFindInt(cellVertex, numSubVertices, subVertices, &subVertex);CHKERRQ(ierr); 501 if (subVertex >= 0) { /* contains submesh vertex */ 502 for (i = 0; i < faceSize; ++i) { 503 if (cellVertex == face[i]) break; 504 } 505 if (i == faceSize) { 506 if (faceSize >= maxConeSize) SETERRQ2(comm, PETSC_ERR_ARG_OUTOFRANGE, "Number of vertices in face %d should not exceed %d", faceSize+1, maxConeSize); 507 face[faceSize] = cellVertex; 508 subface[faceSize] = subVertex; 509 ++faceSize; 510 } 511 } 512 } 513 ierr = DMPlexRestoreTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 514 if (faceSize >= nFV) { 515 if (faceSize > nFV && !boundaryFaces) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Invalid submesh: Too many vertices %d of an element on the surface", faceSize); 516 if (numSubCells >= maxSubCells) { 517 PetscInt *tmpCells; 518 maxSubCells *= 2; 519 ierr = PetscMalloc(maxSubCells * sizeof(PetscInt), &tmpCells);CHKERRQ(ierr); 520 ierr = PetscMemcpy(tmpCells, subCells, numSubCells * sizeof(PetscInt));CHKERRQ(ierr); 521 ierr = PetscFree(subCells);CHKERRQ(ierr); 522 523 subCells = tmpCells; 524 } 525 /* TOOD: Maybe overestimate then squeeze out empty faces */ 526 if (faceSize > nFV) { 527 /* TODO: This is tricky. Maybe just add all faces */ 528 numSubFaces++; 529 } else { 530 numSubFaces++; 531 } 532 for (f = 0; f < faceSize; ++f) subVerticesActive[subface[f]] = 1; 533 subCells[numSubCells++] = cell; 534 } 535 } 536 ierr = DMPlexRestoreTransitiveClosure(dm, vertex, PETSC_FALSE, &starSize, &star);CHKERRQ(ierr); 537 ierr = PetscSortRemoveDupsInt(&numSubCells, subCells);CHKERRQ(ierr); 538 } 539 /* Pick out active subvertices */ 540 for (v = 0, numSubVerticesActive = 0; v < numSubVertices; ++v) { 541 if (subVerticesActive[v]) { 542 subVerticesActive[numSubVerticesActive++] = subVertices[v]; 543 } 544 } 545 ierr = DMPlexSetChart(subdm, 0, numSubCells+numSubFaces+numSubVerticesActive);CHKERRQ(ierr); 546 /* Set cone sizes */ 547 firstSubVertex = numSubCells; 548 firstSubFace = numSubCells+numSubVerticesActive; 549 newFacePoint = firstSubFace; 550 for (c = 0; c < numSubCells; ++c) { 551 ierr = DMPlexSetConeSize(subdm, c, 1);CHKERRQ(ierr); 552 } 553 for (f = firstSubFace; f < firstSubFace+numSubFaces; ++f) { 554 ierr = DMPlexSetConeSize(subdm, f, nFV);CHKERRQ(ierr); 555 } 556 ierr = DMSetUp(subdm);CHKERRQ(ierr); 557 /* Create face cones */ 558 for (c = 0; c < numSubCells; ++c) { 559 const PetscInt cell = subCells[c]; 560 PetscInt *closure = PETSC_NULL; 561 PetscInt closureSize, numCorners = 0, faceSize = 0; 562 563 ierr = DMPlexGetTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 564 for (p = 0; p < closureSize*2; p += 2) { 565 const PetscInt point = closure[p]; 566 if ((point >= vStart) && (point < vEnd)) closure[numCorners++] = point; 567 } 568 for (corner = 0; corner < numCorners; ++corner) { 569 const PetscInt cellVertex = closure[corner]; 570 PetscInt subVertex; 571 572 ierr = PetscFindInt(cellVertex, numSubVerticesActive, subVerticesActive, &subVertex);CHKERRQ(ierr); 573 if (subVertex >= 0) { /* contains submesh vertex */ 574 for (i = 0; i < faceSize; ++i) { 575 if (cellVertex == face[i]) break; 576 } 577 if (i == faceSize) { 578 face[faceSize] = cellVertex; 579 subface[faceSize] = numSubCells+subVertex; 580 ++faceSize; 581 } 582 } 583 } 584 ierr = DMPlexRestoreTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 585 if (faceSize >= nFV) { 586 if (faceSize > nFV && !boundaryFaces) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Invalid submesh: Too many vertices %d of an element on the surface", faceSize); 587 /* Here we allow a set of vertices to lie completely on a boundary cell (like a corner tetrahedron) */ 588 /* We have to take all the faces, and discard those in the interior */ 589 /* We check the join of the face vertices, which produces 2 cells if in the interior */ 590 #if 0 591 /* This object just calls insert on each face that comes from subsets() */ 592 /* In fact, we can just always acll subsets(), since when we pass a single face it is a single call */ 593 FaceInserterV<FlexMesh::sieve_type> inserter(mesh, sieve, subSieve, f, *c_iter, numCorners, indices, &origVertices, &faceVertices, &submeshCells); 594 PointArray faceVec(face->begin(), face->end()); 595 596 subsets(faceVec, nFV, inserter); 597 #endif 598 ierr = DMPlexInsertFace_Private(dm, subdm, faceSize, face, subface, numCorners, cell, c, firstSubFace, &newFacePoint);CHKERRQ(ierr); 599 } 600 } 601 ierr = DMPlexSymmetrize(subdm);CHKERRQ(ierr); 602 ierr = DMPlexStratify(subdm);CHKERRQ(ierr); 603 /* Build coordinates */ 604 ierr = DMPlexGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 605 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr); 606 ierr = DMPlexGetCoordinateSection(subdm, &subCoordSection);CHKERRQ(ierr); 607 ierr = PetscSectionSetChart(subCoordSection, firstSubVertex, firstSubVertex+numSubVerticesActive);CHKERRQ(ierr); 608 for (v = firstSubVertex; v < firstSubVertex+numSubVerticesActive; ++v) { 609 ierr = PetscSectionSetDof(subCoordSection, v, dim);CHKERRQ(ierr); 610 } 611 ierr = PetscSectionSetUp(subCoordSection);CHKERRQ(ierr); 612 ierr = PetscSectionGetStorageSize(subCoordSection, &coordSize);CHKERRQ(ierr); 613 ierr = VecCreate(comm, &subCoordinates);CHKERRQ(ierr); 614 ierr = VecSetSizes(subCoordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 615 ierr = VecSetFromOptions(subCoordinates);CHKERRQ(ierr); 616 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 617 ierr = VecGetArray(subCoordinates, &subCoords);CHKERRQ(ierr); 618 for (v = 0; v < numSubVerticesActive; ++v) { 619 const PetscInt vertex = subVerticesActive[v]; 620 const PetscInt subVertex = firstSubVertex+v; 621 PetscInt dof, off, sdof, soff; 622 623 ierr = PetscSectionGetDof(coordSection, vertex, &dof);CHKERRQ(ierr); 624 ierr = PetscSectionGetOffset(coordSection, vertex, &off);CHKERRQ(ierr); 625 ierr = PetscSectionGetDof(subCoordSection, subVertex, &sdof);CHKERRQ(ierr); 626 ierr = PetscSectionGetOffset(subCoordSection, subVertex, &soff);CHKERRQ(ierr); 627 if (dof != sdof) SETERRQ4(comm, PETSC_ERR_PLIB, "Coordinate dimension %d on subvertex %d, vertex %d should be %d", sdof, subVertex, vertex, dof); 628 for (d = 0; d < dof; ++d) subCoords[soff+d] = coords[off+d]; 629 } 630 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 631 ierr = VecRestoreArray(subCoordinates, &subCoords);CHKERRQ(ierr); 632 ierr = DMSetCoordinatesLocal(subdm, subCoordinates);CHKERRQ(ierr); 633 ierr = VecDestroy(&subCoordinates);CHKERRQ(ierr); 634 635 ierr = DMPlexSetVTKCellHeight(subdm, 1);CHKERRQ(ierr); 636 /* Create map from submesh points to original mesh points */ 637 ierr = DMLabelCreate("subpoint_map", &subpointMap);CHKERRQ(ierr); 638 for (c = 0; c < numSubCells; ++c) {ierr = DMLabelSetValue(subpointMap, subCells[c], 2);CHKERRQ(ierr);} 639 for (v = 0; v < numSubVerticesActive; ++v) {ierr = DMLabelSetValue(subpointMap, subVerticesActive[v], 0);CHKERRQ(ierr);} 640 ierr = DMPlexSetSubpointMap(subdm, subpointMap);CHKERRQ(ierr); 641 ierr = DMLabelDestroy(&subpointMap);CHKERRQ(ierr); 642 643 ierr = PetscFree(subCells);CHKERRQ(ierr); 644 ierr = PetscFree(subVerticesActive);CHKERRQ(ierr); 645 if (labelIS) { 646 ierr = ISRestoreIndices(labelIS, &subVertices);CHKERRQ(ierr); 647 ierr = ISDestroy(&labelIS);CHKERRQ(ierr); 648 } 649 ierr = DMRestoreWorkArray(dm, 2*maxConeSize, PETSC_INT, &face);CHKERRQ(ierr); 650 PetscFunctionReturn(0); 651 } 652 653 #undef __FUNCT__ 654 #define __FUNCT__ "DMPlexCreateSubmesh_Interpolated" 655 static PetscErrorCode DMPlexCreateSubmesh_Interpolated(DM dm, const char vertexLabelName[], DM subdm) 656 { 657 MPI_Comm comm = ((PetscObject) dm)->comm; 658 DMLabel subpointMap, vertexLabel; 659 IS *subpointIS, subvertexIS; 660 const PetscInt **subpoints, *subvertices; 661 PetscInt *pStart, *pEnd, *numSubPoints, *firstSubPoint, *coneNew; 662 PetscInt numSubVerticesInitial, totSubPoints = 0, maxConeSize, dim, p, d, v; 663 PetscErrorCode ierr; 664 665 PetscFunctionBegin; 666 ierr = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr); 667 ierr = PetscMalloc6(dim+1,PetscInt,&pStart,dim+1,PetscInt,&pEnd,dim+1,PetscInt,&numSubPoints,dim+1,PetscInt,&firstSubPoint,dim+1,IS,&subpointIS,dim+1,const PetscInt *,&subpoints);CHKERRQ(ierr); 668 for (d = 0; d <= dim; ++d) { 669 ierr = DMPlexGetDepthStratum(dm, d, &pStart[d], &pEnd[d]);CHKERRQ(ierr); 670 } 671 ierr = DMLabelCreate("subpoint_map", &subpointMap);CHKERRQ(ierr); 672 ierr = DMPlexGetLabel(dm, vertexLabelName, &vertexLabel);CHKERRQ(ierr); 673 ierr = DMLabelGetStratumIS(vertexLabel, 1, &subvertexIS);CHKERRQ(ierr); 674 ierr = ISGetSize(subvertexIS, &numSubVerticesInitial);CHKERRQ(ierr); 675 ierr = ISGetIndices(subvertexIS, &subvertices);CHKERRQ(ierr); 676 /* Loop over initial vertices and mark all faces in the collective star() */ 677 for (v = 0; v < numSubVerticesInitial; ++v) { 678 const PetscInt vertex = subvertices[v]; 679 PetscInt *star = PETSC_NULL; 680 PetscInt starSize, s, numFaces = 0, f; 681 682 ierr = DMPlexGetTransitiveClosure(dm, vertex, PETSC_FALSE, &starSize, &star);CHKERRQ(ierr); 683 for (s = 0; s < starSize*2; s += 2) { 684 const PetscInt point = star[s]; 685 if ((point >= pStart[dim-1]) && (point < pEnd[dim-1])) { 686 star[numFaces++] = point; 687 } 688 } 689 for (f = 0; f < numFaces; ++f) { 690 const PetscInt face = star[f]; 691 PetscInt *closure = PETSC_NULL; 692 PetscInt closureSize, c; 693 PetscInt faceLoc; 694 695 ierr = DMLabelGetValue(subpointMap, face, &faceLoc);CHKERRQ(ierr); 696 if (faceLoc == dim-1) continue; 697 if (faceLoc >= 0) SETERRQ2(comm, PETSC_ERR_PLIB, "Face %d has dimension %d in the surface label", face, faceLoc); 698 ierr = DMPlexGetTransitiveClosure(dm, face, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 699 for (c = 0; c < closureSize*2; c += 2) { 700 const PetscInt point = closure[c]; 701 PetscInt vertexLoc; 702 703 if ((point >= pStart[0]) && (point < pEnd[0])) { 704 ierr = DMLabelGetValue(vertexLabel, point, &vertexLoc);CHKERRQ(ierr); 705 if (vertexLoc < 0) break; 706 } 707 } 708 if (c == closureSize*2) { 709 const PetscInt *support; 710 PetscInt supportSize, s; 711 712 for (c = 0; c < closureSize*2; c += 2) { 713 const PetscInt point = closure[c]; 714 715 for (d = 0; d < dim; ++d) { 716 if ((point >= pStart[d]) && (point < pEnd[d])) { 717 ierr = DMLabelSetValue(subpointMap, point, d);CHKERRQ(ierr); 718 break; 719 } 720 } 721 } 722 ierr = DMPlexGetSupportSize(dm, face, &supportSize);CHKERRQ(ierr); 723 ierr = DMPlexGetSupport(dm, face, &support);CHKERRQ(ierr); 724 for (s = 0; s < supportSize; ++s) { 725 ierr = DMLabelSetValue(subpointMap, support[s], dim);CHKERRQ(ierr); 726 } 727 } 728 ierr = DMPlexRestoreTransitiveClosure(dm, face, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 729 } 730 ierr = DMPlexRestoreTransitiveClosure(dm, vertex, PETSC_FALSE, &starSize, &star);CHKERRQ(ierr); 731 } 732 ierr = ISRestoreIndices(subvertexIS, &subvertices);CHKERRQ(ierr); 733 ierr = ISDestroy(&subvertexIS);CHKERRQ(ierr); 734 for (d = 0; d <= dim; ++d) { 735 ierr = DMLabelGetStratumSize(subpointMap, d, &numSubPoints[d]);CHKERRQ(ierr); 736 totSubPoints += numSubPoints[d]; 737 } 738 ierr = DMPlexSetChart(subdm, 0, totSubPoints);CHKERRQ(ierr); 739 ierr = DMPlexSetVTKCellHeight(subdm, 1);CHKERRQ(ierr); 740 /* Set cone sizes */ 741 firstSubPoint[dim] = 0; 742 firstSubPoint[0] = firstSubPoint[dim] + numSubPoints[dim]; 743 if (dim > 1) {firstSubPoint[dim-1] = firstSubPoint[0] + numSubPoints[0];} 744 if (dim > 2) {firstSubPoint[dim-2] = firstSubPoint[dim-1] + numSubPoints[dim-1];} 745 for (d = 0; d <= dim; ++d) { 746 ierr = DMLabelGetStratumIS(subpointMap, d, &subpointIS[d]);CHKERRQ(ierr); 747 ierr = ISGetIndices(subpointIS[d], &subpoints[d]);CHKERRQ(ierr); 748 } 749 for (d = 0; d <= dim; ++d) { 750 for (p = 0; p < numSubPoints[d]; ++p) { 751 const PetscInt point = subpoints[d][p]; 752 const PetscInt subpoint = firstSubPoint[d] + p; 753 const PetscInt *cone; 754 PetscInt coneSize, coneSizeNew, c, val; 755 756 ierr = DMPlexGetConeSize(dm, point, &coneSize);CHKERRQ(ierr); 757 ierr = DMPlexSetConeSize(subdm, subpoint, coneSize);CHKERRQ(ierr); 758 if (d == dim) { 759 ierr = DMPlexGetCone(dm, point, &cone);CHKERRQ(ierr); 760 for (c = 0, coneSizeNew = 0; c < coneSize; ++c) { 761 ierr = DMLabelGetValue(subpointMap, cone[c], &val);CHKERRQ(ierr); 762 if (val >= 0) coneSizeNew++; 763 } 764 ierr = DMPlexSetConeSize(subdm, subpoint, coneSizeNew);CHKERRQ(ierr); 765 } 766 } 767 } 768 ierr = DMSetUp(subdm);CHKERRQ(ierr); 769 /* Set cones */ 770 ierr = DMPlexGetMaxSizes(dm, &maxConeSize, PETSC_NULL);CHKERRQ(ierr); 771 ierr = PetscMalloc(maxConeSize * sizeof(PetscInt), &coneNew);CHKERRQ(ierr); 772 for (d = 0; d <= dim; ++d) { 773 for (p = 0; p < numSubPoints[d]; ++p) { 774 const PetscInt point = subpoints[d][p]; 775 const PetscInt subpoint = firstSubPoint[d] + p; 776 const PetscInt *cone; 777 PetscInt coneSize, subconeSize, coneSizeNew, c, subc; 778 779 ierr = DMPlexGetConeSize(dm, point, &coneSize);CHKERRQ(ierr); 780 ierr = DMPlexGetConeSize(subdm, subpoint, &subconeSize);CHKERRQ(ierr); 781 ierr = DMPlexGetCone(dm, point, &cone);CHKERRQ(ierr); 782 for (c = 0, coneSizeNew = 0; c < coneSize; ++c) { 783 ierr = PetscFindInt(cone[c], numSubPoints[d-1], subpoints[d-1], &subc);CHKERRQ(ierr); 784 if (subc >= 0) coneNew[coneSizeNew++] = firstSubPoint[d-1] + subc; 785 } 786 if (coneSizeNew != subconeSize) SETERRQ2(comm, PETSC_ERR_PLIB, "Number of cone points located %d does not match subcone size %d", coneSizeNew, subconeSize); 787 ierr = DMPlexSetCone(subdm, subpoint, coneNew);CHKERRQ(ierr); 788 } 789 } 790 ierr = PetscFree(coneNew);CHKERRQ(ierr); 791 ierr = DMPlexSymmetrize(subdm);CHKERRQ(ierr); 792 ierr = DMPlexStratify(subdm);CHKERRQ(ierr); 793 /* Build coordinates */ 794 { 795 PetscSection coordSection, subCoordSection; 796 Vec coordinates, subCoordinates; 797 PetscScalar *coords, *subCoords; 798 PetscInt coordSize; 799 800 ierr = DMPlexGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 801 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr); 802 ierr = DMPlexGetCoordinateSection(subdm, &subCoordSection);CHKERRQ(ierr); 803 ierr = PetscSectionSetChart(subCoordSection, firstSubPoint[0], firstSubPoint[0]+numSubPoints[0]);CHKERRQ(ierr); 804 for (v = 0; v < numSubPoints[0]; ++v) { 805 const PetscInt vertex = subpoints[0][v]; 806 const PetscInt subvertex = firstSubPoint[0]+v; 807 PetscInt dof; 808 809 ierr = PetscSectionGetDof(coordSection, vertex, &dof);CHKERRQ(ierr); 810 ierr = PetscSectionSetDof(subCoordSection, subvertex, dof);CHKERRQ(ierr); 811 } 812 ierr = PetscSectionSetUp(subCoordSection);CHKERRQ(ierr); 813 ierr = PetscSectionGetStorageSize(subCoordSection, &coordSize);CHKERRQ(ierr); 814 ierr = VecCreate(comm, &subCoordinates);CHKERRQ(ierr); 815 ierr = VecSetSizes(subCoordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 816 ierr = VecSetFromOptions(subCoordinates);CHKERRQ(ierr); 817 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 818 ierr = VecGetArray(subCoordinates, &subCoords);CHKERRQ(ierr); 819 for (v = 0; v < numSubPoints[0]; ++v) { 820 const PetscInt vertex = subpoints[0][v]; 821 const PetscInt subvertex = firstSubPoint[0]+v; 822 PetscInt dof, off, sdof, soff, d; 823 824 ierr = PetscSectionGetDof(coordSection, vertex, &dof);CHKERRQ(ierr); 825 ierr = PetscSectionGetOffset(coordSection, vertex, &off);CHKERRQ(ierr); 826 ierr = PetscSectionGetDof(subCoordSection, subvertex, &sdof);CHKERRQ(ierr); 827 ierr = PetscSectionGetOffset(subCoordSection, subvertex, &soff);CHKERRQ(ierr); 828 if (dof != sdof) SETERRQ4(comm, PETSC_ERR_PLIB, "Coordinate dimension %d on subvertex %d, vertex %d should be %d", sdof, subvertex, vertex, dof); 829 for (d = 0; d < dof; ++d) { 830 subCoords[soff+d] = coords[off+d]; 831 } 832 } 833 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 834 ierr = VecRestoreArray(subCoordinates, &subCoords);CHKERRQ(ierr); 835 ierr = DMSetCoordinatesLocal(subdm, subCoordinates);CHKERRQ(ierr); 836 ierr = VecDestroy(&subCoordinates);CHKERRQ(ierr); 837 } 838 for (d = 0; d <= dim; ++d) { 839 ierr = ISRestoreIndices(subpointIS[d], &subpoints[d]);CHKERRQ(ierr); 840 ierr = ISDestroy(&subpointIS[d]);CHKERRQ(ierr); 841 } 842 ierr = DMPlexSetSubpointMap(subdm, subpointMap);CHKERRQ(ierr); 843 ierr = DMLabelDestroy(&subpointMap);CHKERRQ(ierr); 844 PetscFunctionReturn(0); 845 } 846 847 #undef __FUNCT__ 848 #define __FUNCT__ "DMPlexCreateSubmesh" 849 /* 850 DMPlexCreateSubmesh - Extract a hypersurface from the mesh using vertices defined by a label 851 852 Input Parameters: 853 + dm - The original mesh 854 - vertexLabel - The DMLabel marking vertices contained in the surface 855 856 Output Parameter: 857 . subdm - The surface mesh 858 859 Note: This function produces a DMLabel mapping original points in the submesh to their depth. This can be obtained using DMPlexGetSubpointMap(). 860 861 Level: developer 862 863 .seealso: DMPlexGetSubpointMap(), DMPlexGetLabel(), DMLabelSetValue() 864 */ 865 PetscErrorCode DMPlexCreateSubmesh(DM dm, const char vertexLabel[], DM *subdm) 866 { 867 PetscInt dim, depth; 868 PetscErrorCode ierr; 869 870 PetscFunctionBegin; 871 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 872 PetscValidCharPointer(vertexLabel, 2); 873 PetscValidPointer(subdm, 4); 874 ierr = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr); 875 ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr); 876 ierr = DMCreate(((PetscObject) dm)->comm, subdm);CHKERRQ(ierr); 877 ierr = DMSetType(*subdm, DMPLEX);CHKERRQ(ierr); 878 ierr = DMPlexSetDimension(*subdm, dim-1);CHKERRQ(ierr); 879 /* TODO Remove the dim guard */ 880 if (depth == dim) { 881 ierr = DMPlexCreateSubmesh_Interpolated(dm, vertexLabel, *subdm);CHKERRQ(ierr); 882 } else { 883 ierr = DMPlexCreateSubmesh_Uninterpolated(dm, vertexLabel, *subdm);CHKERRQ(ierr); 884 } 885 PetscFunctionReturn(0); 886 } 887 888 #undef __FUNCT__ 889 #define __FUNCT__ "DMPlexGetSubpointMap" 890 PetscErrorCode DMPlexGetSubpointMap(DM dm, DMLabel *subpointMap) 891 { 892 DM_Plex *mesh = (DM_Plex*) dm->data; 893 894 PetscFunctionBegin; 895 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 896 PetscValidPointer(subpointMap, 2); 897 *subpointMap = mesh->subpointMap; 898 PetscFunctionReturn(0); 899 } 900 901 #undef __FUNCT__ 902 #define __FUNCT__ "DMPlexSetSubpointMap" 903 /* Note: Should normally not be called by the user, since it is set in DMPlexCreateSubmesh() */ 904 PetscErrorCode DMPlexSetSubpointMap(DM dm, DMLabel subpointMap) 905 { 906 DM_Plex *mesh = (DM_Plex *) dm->data; 907 PetscErrorCode ierr; 908 909 PetscFunctionBegin; 910 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 911 ierr = DMLabelDestroy(&mesh->subpointMap);CHKERRQ(ierr); 912 mesh->subpointMap = subpointMap; 913 ++mesh->subpointMap->refct;CHKERRQ(ierr); 914 PetscFunctionReturn(0); 915 } 916 917 #undef __FUNCT__ 918 #define __FUNCT__ "DMPlexCreateSubpointIS" 919 /* 920 DMPlexCreateSubpointIS - Creates an IS covering the entire subdm chart with the original points as data 921 922 Input Parameter: 923 . dm - The submesh DM 924 925 Output Parameter: 926 . subpointIS - The IS of all the points from the original mesh in this submesh, or PETSC_NULL if this is not a submesh 927 928 Note: This is IS is guaranteed to be sorted by the construction of the submesh 929 */ 930 PetscErrorCode DMPlexCreateSubpointIS(DM dm, IS *subpointIS) 931 { 932 MPI_Comm comm = ((PetscObject) dm)->comm; 933 DMLabel subpointMap; 934 IS is; 935 const PetscInt *opoints; 936 PetscInt *points, *depths; 937 PetscInt depth, depStart, depEnd, d, pStart, pEnd, p, n, off; 938 PetscErrorCode ierr; 939 940 PetscFunctionBegin; 941 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 942 PetscValidPointer(subpointIS, 2); 943 *subpointIS = PETSC_NULL; 944 ierr = DMPlexGetSubpointMap(dm, &subpointMap);CHKERRQ(ierr); 945 if (subpointMap) { 946 ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr); 947 ierr = DMPlexGetChart(dm, &pStart, &pEnd);CHKERRQ(ierr); 948 if (pStart) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Submeshes must start the point numbering at 0, not %d", pStart); 949 ierr = DMGetWorkArray(dm, depth+1, PETSC_INT, &depths);CHKERRQ(ierr); 950 depths[0] = depth; 951 depths[1] = 0; 952 for(d = 2; d <= depth; ++d) {depths[d] = depth+1 - d;} 953 ierr = PetscMalloc(pEnd * sizeof(PetscInt), &points);CHKERRQ(ierr); 954 for(d = 0, off = 0; d <= depth; ++d) { 955 const PetscInt dep = depths[d]; 956 957 ierr = DMPlexGetDepthStratum(dm, dep, &depStart, &depEnd);CHKERRQ(ierr); 958 ierr = DMLabelGetStratumSize(subpointMap, dep, &n);CHKERRQ(ierr); 959 if (d < 2) { /* Only check vertices and cells for now since the map is broken for others */ 960 if (n != depEnd-depStart) SETERRQ3(comm, PETSC_ERR_ARG_WRONG, "The number of mapped submesh points %d at depth %d should be %d", n, dep, depEnd-depStart); 961 } else { 962 if (!n) for(p = 0; p < depEnd-depStart; ++p, ++off) points[off] = PETSC_MAX_INT; 963 } 964 if (n) { 965 ierr = DMLabelGetStratumIS(subpointMap, dep, &is);CHKERRQ(ierr); 966 ierr = ISGetIndices(is, &opoints);CHKERRQ(ierr); 967 for(p = 0; p < n; ++p, ++off) points[off] = opoints[p]; 968 ierr = ISRestoreIndices(is, &opoints);CHKERRQ(ierr); 969 ierr = ISDestroy(&is);CHKERRQ(ierr); 970 } 971 } 972 ierr = DMRestoreWorkArray(dm, depth+1, PETSC_INT, &depths);CHKERRQ(ierr); 973 if (off != pEnd) SETERRQ2(comm, PETSC_ERR_ARG_WRONG, "The number of mapped submesh points %d should be %d", off, pEnd); 974 ierr = ISCreateGeneral(comm, pEnd, points, PETSC_OWN_POINTER, subpointIS);CHKERRQ(ierr); 975 } 976 PetscFunctionReturn(0); 977 } 978