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