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