1 #include <petsc/private/dmpleximpl.h> /*I "petscdmplex.h" I*/ 2 #include <petscsf.h> 3 4 #undef __FUNCT__ 5 #define __FUNCT__ "DMPlexMarkBoundaryFaces_Internal" 6 PetscErrorCode DMPlexMarkBoundaryFaces_Internal(DM dm, PetscInt cellHeight, DMLabel label) 7 { 8 PetscInt fStart, fEnd, f; 9 PetscErrorCode ierr; 10 11 PetscFunctionBegin; 12 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 13 ierr = DMPlexGetHeightStratum(dm, cellHeight+1, &fStart, &fEnd);CHKERRQ(ierr); 14 for (f = fStart; f < fEnd; ++f) { 15 PetscInt supportSize; 16 17 ierr = DMPlexGetSupportSize(dm, f, &supportSize);CHKERRQ(ierr); 18 if (supportSize == 1) {ierr = DMLabelSetValue(label, f, 1);CHKERRQ(ierr);} 19 } 20 PetscFunctionReturn(0); 21 } 22 23 #undef __FUNCT__ 24 #define __FUNCT__ "DMPlexMarkBoundaryFaces" 25 /*@ 26 DMPlexMarkBoundaryFaces - Mark all faces on the boundary 27 28 Not Collective 29 30 Input Parameter: 31 . dm - The original DM 32 33 Output Parameter: 34 . label - The DMLabel marking boundary faces with value 1 35 36 Level: developer 37 38 .seealso: DMLabelCreate(), DMPlexCreateLabel() 39 @*/ 40 PetscErrorCode DMPlexMarkBoundaryFaces(DM dm, DMLabel label) 41 { 42 PetscErrorCode ierr; 43 44 PetscFunctionBegin; 45 ierr = DMPlexMarkBoundaryFaces_Internal(dm, 0, label);CHKERRQ(ierr); 46 PetscFunctionReturn(0); 47 } 48 49 #undef __FUNCT__ 50 #define __FUNCT__ "DMPlexLabelComplete" 51 /*@ 52 DMPlexLabelComplete - Starting with a label marking points on a surface, we add the transitive closure to the surface 53 54 Input Parameters: 55 + dm - The DM 56 - label - A DMLabel marking the surface points 57 58 Output Parameter: 59 . label - A DMLabel marking all surface points in the transitive closure 60 61 Level: developer 62 63 .seealso: DMPlexLabelCohesiveComplete() 64 @*/ 65 PetscErrorCode DMPlexLabelComplete(DM dm, DMLabel label) 66 { 67 IS valueIS; 68 const PetscInt *values; 69 PetscInt numValues, v; 70 PetscErrorCode ierr; 71 72 PetscFunctionBegin; 73 ierr = DMLabelGetNumValues(label, &numValues);CHKERRQ(ierr); 74 ierr = DMLabelGetValueIS(label, &valueIS);CHKERRQ(ierr); 75 ierr = ISGetIndices(valueIS, &values);CHKERRQ(ierr); 76 for (v = 0; v < numValues; ++v) { 77 IS pointIS; 78 const PetscInt *points; 79 PetscInt numPoints, p; 80 81 ierr = DMLabelGetStratumSize(label, values[v], &numPoints);CHKERRQ(ierr); 82 ierr = DMLabelGetStratumIS(label, values[v], &pointIS);CHKERRQ(ierr); 83 ierr = ISGetIndices(pointIS, &points);CHKERRQ(ierr); 84 for (p = 0; p < numPoints; ++p) { 85 PetscInt *closure = NULL; 86 PetscInt closureSize, c; 87 88 ierr = DMPlexGetTransitiveClosure(dm, points[p], PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 89 for (c = 0; c < closureSize*2; c += 2) { 90 ierr = DMLabelSetValue(label, closure[c], values[v]);CHKERRQ(ierr); 91 } 92 ierr = DMPlexRestoreTransitiveClosure(dm, points[p], PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 93 } 94 ierr = ISRestoreIndices(pointIS, &points);CHKERRQ(ierr); 95 ierr = ISDestroy(&pointIS);CHKERRQ(ierr); 96 } 97 ierr = ISRestoreIndices(valueIS, &values);CHKERRQ(ierr); 98 ierr = ISDestroy(&valueIS);CHKERRQ(ierr); 99 PetscFunctionReturn(0); 100 } 101 102 #undef __FUNCT__ 103 #define __FUNCT__ "DMPlexLabelAddCells" 104 /*@ 105 DMPlexLabelAddCells - Starting with a label marking faces on a surface, we add a cell for each face 106 107 Input Parameters: 108 + dm - The DM 109 - label - A DMLabel marking the surface points 110 111 Output Parameter: 112 . label - A DMLabel incorporating cells 113 114 Level: developer 115 116 Note: The cells allow FEM boundary conditions to be applied using the cell geometry 117 118 .seealso: DMPlexLabelComplete(), DMPlexLabelCohesiveComplete() 119 @*/ 120 PetscErrorCode DMPlexLabelAddCells(DM dm, DMLabel label) 121 { 122 IS valueIS; 123 const PetscInt *values; 124 PetscInt numValues, v, cStart, cEnd, cEndInterior; 125 PetscErrorCode ierr; 126 127 PetscFunctionBegin; 128 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 129 ierr = DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 130 cEnd = cEndInterior < 0 ? cEnd : cEndInterior; 131 ierr = DMLabelGetNumValues(label, &numValues);CHKERRQ(ierr); 132 ierr = DMLabelGetValueIS(label, &valueIS);CHKERRQ(ierr); 133 ierr = ISGetIndices(valueIS, &values);CHKERRQ(ierr); 134 for (v = 0; v < numValues; ++v) { 135 IS pointIS; 136 const PetscInt *points; 137 PetscInt numPoints, p; 138 139 ierr = DMLabelGetStratumSize(label, values[v], &numPoints);CHKERRQ(ierr); 140 ierr = DMLabelGetStratumIS(label, values[v], &pointIS);CHKERRQ(ierr); 141 ierr = ISGetIndices(pointIS, &points);CHKERRQ(ierr); 142 for (p = 0; p < numPoints; ++p) { 143 PetscInt *closure = NULL; 144 PetscInt closureSize, point, cl; 145 146 ierr = DMPlexGetTransitiveClosure(dm, points[p], PETSC_FALSE, &closureSize, &closure);CHKERRQ(ierr); 147 for (cl = closureSize-1; cl > 0; --cl) { 148 point = closure[cl*2]; 149 if ((point >= cStart) && (point < cEnd)) {ierr = DMLabelSetValue(label, point, values[v]);CHKERRQ(ierr); break;} 150 } 151 ierr = DMPlexRestoreTransitiveClosure(dm, points[p], PETSC_FALSE, &closureSize, &closure);CHKERRQ(ierr); 152 } 153 ierr = ISRestoreIndices(pointIS, &points);CHKERRQ(ierr); 154 ierr = ISDestroy(&pointIS);CHKERRQ(ierr); 155 } 156 ierr = ISRestoreIndices(valueIS, &values);CHKERRQ(ierr); 157 ierr = ISDestroy(&valueIS);CHKERRQ(ierr); 158 PetscFunctionReturn(0); 159 } 160 161 #undef __FUNCT__ 162 #define __FUNCT__ "DMPlexShiftPoint_Internal" 163 PETSC_STATIC_INLINE PetscInt DMPlexShiftPoint_Internal(PetscInt p, PetscInt depth, PetscInt depthMax[], PetscInt depthEnd[], PetscInt depthShift[]) 164 { 165 if (depth < 0) return p; 166 /* Normal Cells */ if (p < depthMax[depth]) return p; 167 /* Hybrid Cells+Normal Vertices */ if (p < depthMax[0]) return p + depthShift[depth]; 168 /* Hybrid Vertices+Normal Faces */ if (depth < 2 || p < depthMax[depth-1]) return p + depthShift[depth] + depthShift[0]; 169 /* Hybrid Faces+Normal Edges */ if (depth < 3 || p < depthMax[depth-2]) return p + depthShift[depth] + depthShift[0] + depthShift[depth-1]; 170 /* Hybrid Edges */ return p + depthShift[depth] + depthShift[0] + depthShift[depth-1] + depthShift[depth-2]; 171 } 172 173 #undef __FUNCT__ 174 #define __FUNCT__ "DMPlexShiftSizes_Internal" 175 static PetscErrorCode DMPlexShiftSizes_Internal(DM dm, PetscInt depthShift[], DM dmNew) 176 { 177 PetscInt *depthMax, *depthEnd; 178 PetscInt depth = 0, d, pStart, pEnd, p; 179 PetscErrorCode ierr; 180 181 PetscFunctionBegin; 182 ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr); 183 if (depth < 0) PetscFunctionReturn(0); 184 ierr = PetscMalloc2(depth+1,&depthMax,depth+1,&depthEnd);CHKERRQ(ierr); 185 /* Step 1: Expand chart */ 186 ierr = DMPlexGetChart(dm, &pStart, &pEnd);CHKERRQ(ierr); 187 ierr = DMPlexGetHybridBounds(dm, &depthMax[depth], depth > 0 ? &depthMax[depth-1] : NULL, &depthMax[1], &depthMax[0]);CHKERRQ(ierr); 188 for (d = 0; d <= depth; ++d) { 189 pEnd += depthShift[d]; 190 ierr = DMPlexGetDepthStratum(dm, d, NULL, &depthEnd[d]);CHKERRQ(ierr); 191 depthMax[d] = depthMax[d] < 0 ? depthEnd[d] : depthMax[d]; 192 } 193 ierr = DMPlexSetChart(dmNew, pStart, pEnd);CHKERRQ(ierr); 194 /* Step 2: Set cone and support sizes */ 195 for (d = 0; d <= depth; ++d) { 196 ierr = DMPlexGetDepthStratum(dm, d, &pStart, &pEnd);CHKERRQ(ierr); 197 for (p = pStart; p < pEnd; ++p) { 198 PetscInt newp = DMPlexShiftPoint_Internal(p, depth, depthMax, depthEnd, depthShift); 199 PetscInt size; 200 201 ierr = DMPlexGetConeSize(dm, p, &size);CHKERRQ(ierr); 202 ierr = DMPlexSetConeSize(dmNew, newp, size);CHKERRQ(ierr); 203 ierr = DMPlexGetSupportSize(dm, p, &size);CHKERRQ(ierr); 204 ierr = DMPlexSetSupportSize(dmNew, newp, size);CHKERRQ(ierr); 205 } 206 } 207 ierr = PetscFree2(depthMax,depthEnd);CHKERRQ(ierr); 208 PetscFunctionReturn(0); 209 } 210 211 #undef __FUNCT__ 212 #define __FUNCT__ "DMPlexShiftPoints_Internal" 213 static PetscErrorCode DMPlexShiftPoints_Internal(DM dm, PetscInt depthShift[], DM dmNew) 214 { 215 PetscInt *depthEnd, *depthMax, *newpoints; 216 PetscInt depth = 0, d, maxConeSize, maxSupportSize, maxConeSizeNew, maxSupportSizeNew, pStart, pEnd, p; 217 PetscErrorCode ierr; 218 219 PetscFunctionBegin; 220 ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr); 221 if (depth < 0) PetscFunctionReturn(0); 222 ierr = DMPlexGetMaxSizes(dm, &maxConeSize, &maxSupportSize);CHKERRQ(ierr); 223 ierr = DMPlexGetMaxSizes(dmNew, &maxConeSizeNew, &maxSupportSizeNew);CHKERRQ(ierr); 224 ierr = PetscMalloc3(depth+1,&depthEnd,depth+1,&depthMax,PetscMax(PetscMax(maxConeSize, maxSupportSize), PetscMax(maxConeSizeNew, maxSupportSizeNew)),&newpoints);CHKERRQ(ierr); 225 ierr = DMPlexGetHybridBounds(dm, &depthMax[depth], depth > 0 ? &depthMax[depth-1] : NULL, &depthMax[1], &depthMax[0]);CHKERRQ(ierr); 226 for (d = 0; d <= depth; ++d) { 227 ierr = DMPlexGetDepthStratum(dm, d, NULL, &depthEnd[d]);CHKERRQ(ierr); 228 depthMax[d] = depthMax[d] < 0 ? depthEnd[d] : depthMax[d]; 229 } 230 /* Step 5: Set cones and supports */ 231 ierr = DMPlexGetChart(dm, &pStart, &pEnd);CHKERRQ(ierr); 232 for (p = pStart; p < pEnd; ++p) { 233 const PetscInt *points = NULL, *orientations = NULL; 234 PetscInt size,sizeNew, i, newp = DMPlexShiftPoint_Internal(p, depth, depthMax, depthEnd, depthShift); 235 236 ierr = DMPlexGetConeSize(dm, p, &size);CHKERRQ(ierr); 237 ierr = DMPlexGetCone(dm, p, &points);CHKERRQ(ierr); 238 ierr = DMPlexGetConeOrientation(dm, p, &orientations);CHKERRQ(ierr); 239 for (i = 0; i < size; ++i) { 240 newpoints[i] = DMPlexShiftPoint_Internal(points[i], depth, depthMax, depthEnd, depthShift); 241 } 242 ierr = DMPlexSetCone(dmNew, newp, newpoints);CHKERRQ(ierr); 243 ierr = DMPlexSetConeOrientation(dmNew, newp, orientations);CHKERRQ(ierr); 244 ierr = DMPlexGetSupportSize(dm, p, &size);CHKERRQ(ierr); 245 ierr = DMPlexGetSupportSize(dmNew, newp, &sizeNew);CHKERRQ(ierr); 246 ierr = DMPlexGetSupport(dm, p, &points);CHKERRQ(ierr); 247 for (i = 0; i < size; ++i) { 248 newpoints[i] = DMPlexShiftPoint_Internal(points[i], depth, depthMax, depthEnd, depthShift); 249 } 250 for (i = size; i < sizeNew; ++i) newpoints[i] = 0; 251 ierr = DMPlexSetSupport(dmNew, newp, newpoints);CHKERRQ(ierr); 252 } 253 ierr = PetscFree3(depthEnd,depthMax,newpoints);CHKERRQ(ierr); 254 PetscFunctionReturn(0); 255 } 256 257 #undef __FUNCT__ 258 #define __FUNCT__ "DMPlexShiftCoordinates_Internal" 259 static PetscErrorCode DMPlexShiftCoordinates_Internal(DM dm, PetscInt depthShift[], DM dmNew) 260 { 261 PetscSection coordSection, newCoordSection; 262 Vec coordinates, newCoordinates; 263 PetscScalar *coords, *newCoords; 264 PetscInt *depthEnd, coordSize; 265 PetscInt dim, depth = 0, d, vStart, vEnd, vStartNew, vEndNew, v; 266 PetscErrorCode ierr; 267 268 PetscFunctionBegin; 269 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 270 ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr); 271 ierr = PetscMalloc1(depth+1, &depthEnd);CHKERRQ(ierr); 272 for (d = 0; d <= depth; ++d) { 273 ierr = DMPlexGetDepthStratum(dm, d, NULL, &depthEnd[d]);CHKERRQ(ierr); 274 } 275 /* Step 8: Convert coordinates */ 276 ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr); 277 ierr = DMPlexGetDepthStratum(dmNew, 0, &vStartNew, &vEndNew);CHKERRQ(ierr); 278 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 279 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)dm), &newCoordSection);CHKERRQ(ierr); 280 ierr = PetscSectionSetNumFields(newCoordSection, 1);CHKERRQ(ierr); 281 ierr = PetscSectionSetFieldComponents(newCoordSection, 0, dim);CHKERRQ(ierr); 282 ierr = PetscSectionSetChart(newCoordSection, vStartNew, vEndNew);CHKERRQ(ierr); 283 for (v = vStartNew; v < vEndNew; ++v) { 284 ierr = PetscSectionSetDof(newCoordSection, v, dim);CHKERRQ(ierr); 285 ierr = PetscSectionSetFieldDof(newCoordSection, v, 0, dim);CHKERRQ(ierr); 286 } 287 ierr = PetscSectionSetUp(newCoordSection);CHKERRQ(ierr); 288 ierr = DMSetCoordinateSection(dmNew, PETSC_DETERMINE, newCoordSection);CHKERRQ(ierr); 289 ierr = PetscSectionGetStorageSize(newCoordSection, &coordSize);CHKERRQ(ierr); 290 ierr = VecCreate(PetscObjectComm((PetscObject)dm), &newCoordinates);CHKERRQ(ierr); 291 ierr = PetscObjectSetName((PetscObject) newCoordinates, "coordinates");CHKERRQ(ierr); 292 ierr = VecSetSizes(newCoordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 293 ierr = VecSetBlockSize(newCoordinates, dim);CHKERRQ(ierr); 294 ierr = VecSetType(newCoordinates,VECSTANDARD);CHKERRQ(ierr); 295 ierr = DMSetCoordinatesLocal(dmNew, newCoordinates);CHKERRQ(ierr); 296 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr); 297 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 298 ierr = VecGetArray(newCoordinates, &newCoords);CHKERRQ(ierr); 299 for (v = vStart; v < vEnd; ++v) { 300 PetscInt dof, off, noff, d; 301 302 ierr = PetscSectionGetDof(coordSection, v, &dof);CHKERRQ(ierr); 303 ierr = PetscSectionGetOffset(coordSection, v, &off);CHKERRQ(ierr); 304 ierr = PetscSectionGetOffset(newCoordSection, DMPlexShiftPoint_Internal(v, depth, depthEnd, depthEnd, depthShift), &noff);CHKERRQ(ierr); 305 for (d = 0; d < dof; ++d) { 306 newCoords[noff+d] = coords[off+d]; 307 } 308 } 309 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 310 ierr = VecRestoreArray(newCoordinates, &newCoords);CHKERRQ(ierr); 311 ierr = VecDestroy(&newCoordinates);CHKERRQ(ierr); 312 ierr = PetscSectionDestroy(&newCoordSection);CHKERRQ(ierr); 313 ierr = PetscFree(depthEnd);CHKERRQ(ierr); 314 PetscFunctionReturn(0); 315 } 316 317 #undef __FUNCT__ 318 #define __FUNCT__ "DMPlexShiftSF_Internal" 319 static PetscErrorCode DMPlexShiftSF_Internal(DM dm, PetscInt depthShift[], DM dmNew) 320 { 321 PetscInt *depthMax, *depthEnd; 322 PetscInt depth = 0, d; 323 PetscSF sfPoint, sfPointNew; 324 const PetscSFNode *remotePoints; 325 PetscSFNode *gremotePoints; 326 const PetscInt *localPoints; 327 PetscInt *glocalPoints, *newLocation, *newRemoteLocation; 328 PetscInt numRoots, numLeaves, l, pStart, pEnd, totShift = 0; 329 PetscErrorCode ierr; 330 331 PetscFunctionBegin; 332 ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr); 333 ierr = PetscMalloc2(depth+1,&depthMax,depth+1,&depthEnd);CHKERRQ(ierr); 334 ierr = DMPlexGetHybridBounds(dm, depth >= 0 ? &depthMax[depth] : NULL, depth>1 ? &depthMax[depth-1] : NULL, depth>2 ? &depthMax[1] : NULL, depth >= 0 ? &depthMax[0] : NULL);CHKERRQ(ierr); 335 for (d = 0; d <= depth; ++d) { 336 totShift += depthShift[d]; 337 ierr = DMPlexGetDepthStratum(dm, d, NULL, &depthEnd[d]);CHKERRQ(ierr); 338 depthMax[d] = depthMax[d] < 0 ? depthEnd[d] : depthMax[d]; 339 } 340 /* Step 9: Convert pointSF */ 341 ierr = DMGetPointSF(dm, &sfPoint);CHKERRQ(ierr); 342 ierr = DMGetPointSF(dmNew, &sfPointNew);CHKERRQ(ierr); 343 ierr = DMPlexGetChart(dm, &pStart, &pEnd);CHKERRQ(ierr); 344 ierr = PetscSFGetGraph(sfPoint, &numRoots, &numLeaves, &localPoints, &remotePoints);CHKERRQ(ierr); 345 if (numRoots >= 0) { 346 ierr = PetscMalloc2(numRoots,&newLocation,pEnd-pStart,&newRemoteLocation);CHKERRQ(ierr); 347 for (l=0; l<numRoots; l++) newLocation[l] = DMPlexShiftPoint_Internal(l, depth, depthMax, depthEnd, depthShift); 348 ierr = PetscSFBcastBegin(sfPoint, MPIU_INT, newLocation, newRemoteLocation);CHKERRQ(ierr); 349 ierr = PetscSFBcastEnd(sfPoint, MPIU_INT, newLocation, newRemoteLocation);CHKERRQ(ierr); 350 ierr = PetscMalloc1(numLeaves, &glocalPoints);CHKERRQ(ierr); 351 ierr = PetscMalloc1(numLeaves, &gremotePoints);CHKERRQ(ierr); 352 for (l = 0; l < numLeaves; ++l) { 353 glocalPoints[l] = DMPlexShiftPoint_Internal(localPoints[l], depth, depthMax, depthEnd, depthShift); 354 gremotePoints[l].rank = remotePoints[l].rank; 355 gremotePoints[l].index = newRemoteLocation[localPoints[l]]; 356 } 357 ierr = PetscFree2(newLocation,newRemoteLocation);CHKERRQ(ierr); 358 ierr = PetscSFSetGraph(sfPointNew, numRoots + totShift, numLeaves, glocalPoints, PETSC_OWN_POINTER, gremotePoints, PETSC_OWN_POINTER);CHKERRQ(ierr); 359 } 360 ierr = PetscFree2(depthMax,depthEnd);CHKERRQ(ierr); 361 PetscFunctionReturn(0); 362 } 363 364 #undef __FUNCT__ 365 #define __FUNCT__ "DMPlexShiftLabels_Internal" 366 static PetscErrorCode DMPlexShiftLabels_Internal(DM dm, PetscInt depthShift[], DM dmNew) 367 { 368 PetscSF sfPoint; 369 DMLabel vtkLabel, ghostLabel; 370 PetscInt *depthMax, *depthEnd; 371 const PetscSFNode *leafRemote; 372 const PetscInt *leafLocal; 373 PetscInt depth = 0, d, numLeaves, numLabels, l, cStart, cEnd, c, fStart, fEnd, f; 374 PetscMPIInt rank; 375 PetscErrorCode ierr; 376 377 PetscFunctionBegin; 378 ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr); 379 ierr = PetscMalloc2(depth+1,&depthMax,depth+1,&depthEnd);CHKERRQ(ierr); 380 ierr = DMPlexGetHybridBounds(dm, depth >= 0 ? &depthMax[depth] : NULL, depth>1 ? &depthMax[depth-1] : NULL, depth>2 ? &depthMax[1] : NULL, depth >= 0 ? &depthMax[0] : NULL);CHKERRQ(ierr); 381 for (d = 0; d <= depth; ++d) { 382 ierr = DMPlexGetDepthStratum(dm, d, NULL, &depthEnd[d]);CHKERRQ(ierr); 383 depthMax[d] = depthMax[d] < 0 ? depthEnd[d] : depthMax[d]; 384 } 385 /* Step 10: Convert labels */ 386 ierr = DMPlexGetNumLabels(dm, &numLabels);CHKERRQ(ierr); 387 for (l = 0; l < numLabels; ++l) { 388 DMLabel label, newlabel; 389 const char *lname; 390 PetscBool isDepth; 391 IS valueIS; 392 const PetscInt *values; 393 PetscInt numValues, val; 394 395 ierr = DMPlexGetLabelName(dm, l, &lname);CHKERRQ(ierr); 396 ierr = PetscStrcmp(lname, "depth", &isDepth);CHKERRQ(ierr); 397 if (isDepth) continue; 398 ierr = DMPlexCreateLabel(dmNew, lname);CHKERRQ(ierr); 399 ierr = DMPlexGetLabel(dm, lname, &label);CHKERRQ(ierr); 400 ierr = DMPlexGetLabel(dmNew, lname, &newlabel);CHKERRQ(ierr); 401 ierr = DMLabelGetValueIS(label, &valueIS);CHKERRQ(ierr); 402 ierr = ISGetLocalSize(valueIS, &numValues);CHKERRQ(ierr); 403 ierr = ISGetIndices(valueIS, &values);CHKERRQ(ierr); 404 for (val = 0; val < numValues; ++val) { 405 IS pointIS; 406 const PetscInt *points; 407 PetscInt numPoints, p; 408 409 ierr = DMLabelGetStratumIS(label, values[val], &pointIS);CHKERRQ(ierr); 410 ierr = ISGetLocalSize(pointIS, &numPoints);CHKERRQ(ierr); 411 ierr = ISGetIndices(pointIS, &points);CHKERRQ(ierr); 412 for (p = 0; p < numPoints; ++p) { 413 const PetscInt newpoint = DMPlexShiftPoint_Internal(points[p], depth, depthMax, depthEnd, depthShift); 414 415 ierr = DMLabelSetValue(newlabel, newpoint, values[val]);CHKERRQ(ierr); 416 } 417 ierr = ISRestoreIndices(pointIS, &points);CHKERRQ(ierr); 418 ierr = ISDestroy(&pointIS);CHKERRQ(ierr); 419 } 420 ierr = ISRestoreIndices(valueIS, &values);CHKERRQ(ierr); 421 ierr = ISDestroy(&valueIS);CHKERRQ(ierr); 422 } 423 ierr = PetscFree2(depthMax,depthEnd);CHKERRQ(ierr); 424 /* Step 11: Make label for output (vtk) and to mark ghost points (ghost) */ 425 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);CHKERRQ(ierr); 426 ierr = DMGetPointSF(dm, &sfPoint);CHKERRQ(ierr); 427 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 428 ierr = PetscSFGetGraph(sfPoint, NULL, &numLeaves, &leafLocal, &leafRemote);CHKERRQ(ierr); 429 ierr = DMPlexCreateLabel(dmNew, "vtk");CHKERRQ(ierr); 430 ierr = DMPlexCreateLabel(dmNew, "ghost");CHKERRQ(ierr); 431 ierr = DMPlexGetLabel(dmNew, "vtk", &vtkLabel);CHKERRQ(ierr); 432 ierr = DMPlexGetLabel(dmNew, "ghost", &ghostLabel);CHKERRQ(ierr); 433 for (l = 0, c = cStart; l < numLeaves && c < cEnd; ++l, ++c) { 434 for (; c < leafLocal[l] && c < cEnd; ++c) { 435 ierr = DMLabelSetValue(vtkLabel, c, 1);CHKERRQ(ierr); 436 } 437 if (leafLocal[l] >= cEnd) break; 438 if (leafRemote[l].rank == rank) { 439 ierr = DMLabelSetValue(vtkLabel, c, 1);CHKERRQ(ierr); 440 } else { 441 ierr = DMLabelSetValue(ghostLabel, c, 2);CHKERRQ(ierr); 442 } 443 } 444 for (; c < cEnd; ++c) { 445 ierr = DMLabelSetValue(vtkLabel, c, 1);CHKERRQ(ierr); 446 } 447 if (0) { 448 ierr = DMLabelView(vtkLabel, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 449 } 450 ierr = DMPlexGetHeightStratum(dmNew, 1, &fStart, &fEnd);CHKERRQ(ierr); 451 for (f = fStart; f < fEnd; ++f) { 452 PetscInt numCells; 453 454 ierr = DMPlexGetSupportSize(dmNew, f, &numCells);CHKERRQ(ierr); 455 if (numCells < 2) { 456 ierr = DMLabelSetValue(ghostLabel, f, 1);CHKERRQ(ierr); 457 } else { 458 const PetscInt *cells = NULL; 459 PetscInt vA, vB; 460 461 ierr = DMPlexGetSupport(dmNew, f, &cells);CHKERRQ(ierr); 462 ierr = DMLabelGetValue(vtkLabel, cells[0], &vA);CHKERRQ(ierr); 463 ierr = DMLabelGetValue(vtkLabel, cells[1], &vB);CHKERRQ(ierr); 464 if (!vA && !vB) {ierr = DMLabelSetValue(ghostLabel, f, 1);CHKERRQ(ierr);} 465 } 466 } 467 if (0) { 468 ierr = DMLabelView(ghostLabel, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 469 } 470 PetscFunctionReturn(0); 471 } 472 473 #undef __FUNCT__ 474 #define __FUNCT__ "DMPlexConstructGhostCells_Internal" 475 static PetscErrorCode DMPlexConstructGhostCells_Internal(DM dm, DMLabel label, PetscInt *numGhostCells, DM gdm) 476 { 477 PetscSF sf; 478 IS valueIS; 479 const PetscInt *values, *leaves; 480 PetscInt *depthShift; 481 PetscInt depth = 0, nleaves, loc, Ng, numFS, fs, fStart, fEnd, ghostCell, cEnd, c; 482 PetscErrorCode ierr; 483 484 PetscFunctionBegin; 485 ierr = DMGetPointSF(dm, &sf);CHKERRQ(ierr); 486 ierr = PetscSFGetGraph(sf, NULL, &nleaves, &leaves, NULL);CHKERRQ(ierr); 487 nleaves = PetscMax(0, nleaves); 488 ierr = DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);CHKERRQ(ierr); 489 /* Count ghost cells */ 490 ierr = DMLabelGetValueIS(label, &valueIS);CHKERRQ(ierr); 491 ierr = ISGetLocalSize(valueIS, &numFS);CHKERRQ(ierr); 492 ierr = ISGetIndices(valueIS, &values);CHKERRQ(ierr); 493 Ng = 0; 494 for (fs = 0; fs < numFS; ++fs) { 495 IS faceIS; 496 const PetscInt *faces; 497 PetscInt numFaces, f, numBdFaces = 0; 498 499 ierr = DMLabelGetStratumIS(label, values[fs], &faceIS);CHKERRQ(ierr); 500 ierr = ISGetLocalSize(faceIS, &numFaces);CHKERRQ(ierr); 501 ierr = ISGetIndices(faceIS, &faces);CHKERRQ(ierr); 502 for (f = 0; f < numFaces; ++f) { 503 ierr = PetscFindInt(faces[f], nleaves, leaves, &loc);CHKERRQ(ierr); 504 if (loc >= 0) continue; 505 if ((faces[f] >= fStart) && (faces[f] < fEnd)) ++numBdFaces; 506 } 507 Ng += numBdFaces; 508 ierr = ISDestroy(&faceIS);CHKERRQ(ierr); 509 } 510 ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr); 511 ierr = PetscMalloc1(depth+1, &depthShift);CHKERRQ(ierr); 512 ierr = PetscMemzero(depthShift, (depth+1) * sizeof(PetscInt));CHKERRQ(ierr); 513 if (depth >= 0) depthShift[depth] = Ng; 514 ierr = DMPlexShiftSizes_Internal(dm, depthShift, gdm);CHKERRQ(ierr); 515 /* Step 3: Set cone/support sizes for new points */ 516 ierr = DMPlexGetHeightStratum(dm, 0, NULL, &cEnd);CHKERRQ(ierr); 517 ierr = DMPlexSetHybridBounds(gdm, cEnd, PETSC_DETERMINE, PETSC_DETERMINE, PETSC_DETERMINE);CHKERRQ(ierr); 518 for (c = cEnd; c < cEnd + Ng; ++c) { 519 ierr = DMPlexSetConeSize(gdm, c, 1);CHKERRQ(ierr); 520 } 521 for (fs = 0; fs < numFS; ++fs) { 522 IS faceIS; 523 const PetscInt *faces; 524 PetscInt numFaces, f; 525 526 ierr = DMLabelGetStratumIS(label, values[fs], &faceIS);CHKERRQ(ierr); 527 ierr = ISGetLocalSize(faceIS, &numFaces);CHKERRQ(ierr); 528 ierr = ISGetIndices(faceIS, &faces);CHKERRQ(ierr); 529 for (f = 0; f < numFaces; ++f) { 530 PetscInt size; 531 532 ierr = PetscFindInt(faces[f], nleaves, leaves, &loc);CHKERRQ(ierr); 533 if (loc >= 0) continue; 534 if ((faces[f] < fStart) || (faces[f] >= fEnd)) continue; 535 ierr = DMPlexGetSupportSize(dm, faces[f], &size);CHKERRQ(ierr); 536 if (size != 1) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "DM has boundary face %d with %d support cells", faces[f], size); 537 ierr = DMPlexSetSupportSize(gdm, faces[f] + Ng, 2);CHKERRQ(ierr); 538 } 539 ierr = ISRestoreIndices(faceIS, &faces);CHKERRQ(ierr); 540 ierr = ISDestroy(&faceIS);CHKERRQ(ierr); 541 } 542 /* Step 4: Setup ghosted DM */ 543 ierr = DMSetUp(gdm);CHKERRQ(ierr); 544 ierr = DMPlexShiftPoints_Internal(dm, depthShift, gdm);CHKERRQ(ierr); 545 /* Step 6: Set cones and supports for new points */ 546 ghostCell = cEnd; 547 for (fs = 0; fs < numFS; ++fs) { 548 IS faceIS; 549 const PetscInt *faces; 550 PetscInt numFaces, f; 551 552 ierr = DMLabelGetStratumIS(label, values[fs], &faceIS);CHKERRQ(ierr); 553 ierr = ISGetLocalSize(faceIS, &numFaces);CHKERRQ(ierr); 554 ierr = ISGetIndices(faceIS, &faces);CHKERRQ(ierr); 555 for (f = 0; f < numFaces; ++f) { 556 PetscInt newFace = faces[f] + Ng; 557 558 ierr = PetscFindInt(faces[f], nleaves, leaves, &loc);CHKERRQ(ierr); 559 if (loc >= 0) continue; 560 if ((faces[f] < fStart) || (faces[f] >= fEnd)) continue; 561 ierr = DMPlexSetCone(gdm, ghostCell, &newFace);CHKERRQ(ierr); 562 ierr = DMPlexInsertSupport(gdm, newFace, 1, ghostCell);CHKERRQ(ierr); 563 ++ghostCell; 564 } 565 ierr = ISRestoreIndices(faceIS, &faces);CHKERRQ(ierr); 566 ierr = ISDestroy(&faceIS);CHKERRQ(ierr); 567 } 568 ierr = ISRestoreIndices(valueIS, &values);CHKERRQ(ierr); 569 ierr = ISDestroy(&valueIS);CHKERRQ(ierr); 570 /* Step 7: Stratify */ 571 ierr = DMPlexStratify(gdm);CHKERRQ(ierr); 572 ierr = DMPlexShiftCoordinates_Internal(dm, depthShift, gdm);CHKERRQ(ierr); 573 ierr = DMPlexShiftSF_Internal(dm, depthShift, gdm);CHKERRQ(ierr); 574 ierr = DMPlexShiftLabels_Internal(dm, depthShift, gdm);CHKERRQ(ierr); 575 ierr = PetscFree(depthShift);CHKERRQ(ierr); 576 /* Step 7: Periodicity */ 577 if (dm->maxCell) { 578 const PetscReal *maxCell, *L; 579 const DMBoundaryType *bd; 580 ierr = DMGetPeriodicity(dm, &maxCell, &L, &bd);CHKERRQ(ierr); 581 ierr = DMSetPeriodicity(gdm, maxCell, L, bd);CHKERRQ(ierr); 582 } 583 if (numGhostCells) *numGhostCells = Ng; 584 PetscFunctionReturn(0); 585 } 586 587 #undef __FUNCT__ 588 #define __FUNCT__ "DMPlexConstructGhostCells" 589 /*@C 590 DMPlexConstructGhostCells - Construct ghost cells which connect to every boundary face 591 592 Collective on dm 593 594 Input Parameters: 595 + dm - The original DM 596 - labelName - The label specifying the boundary faces, or "Face Sets" if this is NULL 597 598 Output Parameters: 599 + numGhostCells - The number of ghost cells added to the DM 600 - dmGhosted - The new DM 601 602 Note: If no label exists of that name, one will be created marking all boundary faces 603 604 Level: developer 605 606 .seealso: DMCreate() 607 @*/ 608 PetscErrorCode DMPlexConstructGhostCells(DM dm, const char labelName[], PetscInt *numGhostCells, DM *dmGhosted) 609 { 610 DM gdm; 611 DMLabel label; 612 const char *name = labelName ? labelName : "Face Sets"; 613 PetscInt dim; 614 PetscBool flag; 615 PetscErrorCode ierr; 616 617 PetscFunctionBegin; 618 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 619 if (numGhostCells) PetscValidPointer(numGhostCells, 3); 620 PetscValidPointer(dmGhosted, 4); 621 ierr = DMCreate(PetscObjectComm((PetscObject)dm), &gdm);CHKERRQ(ierr); 622 ierr = DMSetType(gdm, DMPLEX);CHKERRQ(ierr); 623 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 624 ierr = DMSetDimension(gdm, dim);CHKERRQ(ierr); 625 ierr = DMPlexGetAdjacencyUseCone(dm, &flag);CHKERRQ(ierr); 626 ierr = DMPlexSetAdjacencyUseCone(gdm, flag);CHKERRQ(ierr); 627 ierr = DMPlexGetAdjacencyUseClosure(dm, &flag);CHKERRQ(ierr); 628 ierr = DMPlexSetAdjacencyUseClosure(gdm, flag);CHKERRQ(ierr); 629 ierr = DMPlexGetLabel(dm, name, &label);CHKERRQ(ierr); 630 if (!label) { 631 /* Get label for boundary faces */ 632 ierr = DMPlexCreateLabel(dm, name);CHKERRQ(ierr); 633 ierr = DMPlexGetLabel(dm, name, &label);CHKERRQ(ierr); 634 ierr = DMPlexMarkBoundaryFaces(dm, label);CHKERRQ(ierr); 635 } 636 ierr = DMPlexConstructGhostCells_Internal(dm, label, numGhostCells, gdm);CHKERRQ(ierr); 637 ierr = DMPlexCopyBoundary(dm, gdm);CHKERRQ(ierr); 638 *dmGhosted = gdm; 639 PetscFunctionReturn(0); 640 } 641 642 #undef __FUNCT__ 643 #define __FUNCT__ "DMPlexConstructCohesiveCells_Internal" 644 /* 645 We are adding three kinds of points here: 646 Replicated: Copies of points which exist in the mesh, such as vertices identified across a fault 647 Non-replicated: Points which exist on the fault, but are not replicated 648 Hybrid: Entirely new points, such as cohesive cells 649 650 When creating subsequent cohesive cells, we shift the old hybrid cells to the end of the numbering at 651 each depth so that the new split/hybrid points can be inserted as a block. 652 */ 653 static PetscErrorCode DMPlexConstructCohesiveCells_Internal(DM dm, DMLabel label, DM sdm) 654 { 655 MPI_Comm comm; 656 IS valueIS; 657 PetscInt numSP = 0; /* The number of depths for which we have replicated points */ 658 const PetscInt *values; /* List of depths for which we have replicated points */ 659 IS *splitIS; 660 IS *unsplitIS; 661 PetscInt *numSplitPoints; /* The number of replicated points at each depth */ 662 PetscInt *numUnsplitPoints; /* The number of non-replicated points at each depth which still give rise to hybrid points */ 663 PetscInt *numHybridPoints; /* The number of new hybrid points at each depth */ 664 PetscInt *numHybridPointsOld; /* The number of existing hybrid points at each depth */ 665 const PetscInt **splitPoints; /* Replicated points for each depth */ 666 const PetscInt **unsplitPoints; /* Non-replicated points for each depth */ 667 PetscSection coordSection; 668 Vec coordinates; 669 PetscScalar *coords; 670 PetscInt depths[4]; /* Depths in the order that plex points are numbered */ 671 PetscInt *depthMax; /* The first hybrid point at each depth in the original mesh */ 672 PetscInt *depthEnd; /* The point limit at each depth in the original mesh */ 673 PetscInt *depthShift; /* Number of replicated+hybrid points at each depth */ 674 PetscInt *depthOffset; /* Prefix sums of depthShift */ 675 PetscInt *pMaxNew; /* The first replicated point at each depth in the new mesh, hybrids come after this */ 676 PetscInt *coneNew, *coneONew, *supportNew; 677 PetscInt shift = 100, shift2 = 200, depth = 0, dep, dim, d, sp, maxConeSize, maxSupportSize, maxConeSizeNew, maxSupportSizeNew, numLabels, vStart, vEnd, pEnd, p, v; 678 PetscErrorCode ierr; 679 680 PetscFunctionBegin; 681 ierr = PetscObjectGetComm((PetscObject)dm,&comm);CHKERRQ(ierr); 682 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 683 ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr); 684 ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr); 685 depths[0] = depth; 686 depths[1] = 0; 687 depths[2] = depth-1; 688 depths[3] = 1; 689 /* Count split points and add cohesive cells */ 690 ierr = DMPlexGetMaxSizes(dm, &maxConeSize, &maxSupportSize);CHKERRQ(ierr); 691 ierr = PetscMalloc6(depth+1,&depthMax,depth+1,&depthEnd,depth+1,&depthShift,depth+1,&depthOffset,depth+1,&pMaxNew,depth+1,&numHybridPointsOld);CHKERRQ(ierr); 692 ierr = PetscMalloc7(depth+1,&splitIS,depth+1,&unsplitIS,depth+1,&numSplitPoints,depth+1,&numUnsplitPoints,depth+1,&numHybridPoints,depth+1,&splitPoints,depth+1,&unsplitPoints);CHKERRQ(ierr); 693 ierr = PetscMemzero(depthShift, (depth+1) * sizeof(PetscInt));CHKERRQ(ierr); 694 ierr = PetscMemzero(depthOffset, (depth+1) * sizeof(PetscInt));CHKERRQ(ierr); 695 ierr = DMPlexGetHybridBounds(dm, depth >= 0 ? &depthMax[depth] : NULL, depth>1 ? &depthMax[depth-1] : NULL, depth>2 ? &depthMax[1] : NULL, depth >= 0 ? &depthMax[0] : NULL);CHKERRQ(ierr); 696 for (d = 0; d <= depth; ++d) { 697 ierr = DMPlexGetDepthStratum(dm, d, NULL, &pMaxNew[d]);CHKERRQ(ierr); 698 depthEnd[d] = pMaxNew[d]; 699 depthMax[d] = depthMax[d] < 0 ? depthEnd[d] : depthMax[d]; 700 numSplitPoints[d] = 0; 701 numUnsplitPoints[d] = 0; 702 numHybridPoints[d] = 0; 703 numHybridPointsOld[d] = depthMax[d] < 0 ? 0 : depthEnd[d] - depthMax[d]; 704 splitPoints[d] = NULL; 705 unsplitPoints[d] = NULL; 706 splitIS[d] = NULL; 707 unsplitIS[d] = NULL; 708 } 709 if (label) { 710 ierr = DMLabelGetValueIS(label, &valueIS);CHKERRQ(ierr); 711 ierr = ISGetLocalSize(valueIS, &numSP);CHKERRQ(ierr); 712 ierr = ISGetIndices(valueIS, &values);CHKERRQ(ierr); 713 } 714 for (sp = 0; sp < numSP; ++sp) { 715 const PetscInt dep = values[sp]; 716 717 if ((dep < 0) || (dep > depth)) continue; 718 ierr = DMLabelGetStratumIS(label, dep, &splitIS[dep]);CHKERRQ(ierr); 719 if (splitIS[dep]) { 720 ierr = ISGetLocalSize(splitIS[dep], &numSplitPoints[dep]);CHKERRQ(ierr); 721 ierr = ISGetIndices(splitIS[dep], &splitPoints[dep]);CHKERRQ(ierr); 722 } 723 ierr = DMLabelGetStratumIS(label, shift2+dep, &unsplitIS[dep]);CHKERRQ(ierr); 724 if (unsplitIS[dep]) { 725 ierr = ISGetLocalSize(unsplitIS[dep], &numUnsplitPoints[dep]);CHKERRQ(ierr); 726 ierr = ISGetIndices(unsplitIS[dep], &unsplitPoints[dep]);CHKERRQ(ierr); 727 } 728 } 729 /* Calculate number of hybrid points */ 730 for (d = 1; d <= depth; ++d) numHybridPoints[d] = numSplitPoints[d-1] + numUnsplitPoints[d-1]; /* There is a hybrid cell/face/edge for every split face/edge/vertex */ 731 for (d = 0; d <= depth; ++d) depthShift[d] = numSplitPoints[d] + numHybridPoints[d]; 732 for (d = 1; d <= depth; ++d) depthOffset[depths[d]] = depthOffset[depths[d-1]] + depthShift[depths[d-1]]; 733 for (d = 0; d <= depth; ++d) pMaxNew[d] += depthOffset[d] - numHybridPointsOld[d]; 734 ierr = DMPlexShiftSizes_Internal(dm, depthShift, sdm);CHKERRQ(ierr); 735 /* Step 3: Set cone/support sizes for new points */ 736 for (dep = 0; dep <= depth; ++dep) { 737 for (p = 0; p < numSplitPoints[dep]; ++p) { 738 const PetscInt oldp = splitPoints[dep][p]; 739 const PetscInt newp = DMPlexShiftPoint_Internal(oldp, depth, depthMax, depthEnd, depthShift) /*oldp + depthOffset[dep]*/; 740 const PetscInt splitp = p + pMaxNew[dep]; 741 const PetscInt *support; 742 PetscInt coneSize, supportSize, qf, qn, qp, e; 743 744 ierr = DMPlexGetConeSize(dm, oldp, &coneSize);CHKERRQ(ierr); 745 ierr = DMPlexSetConeSize(sdm, splitp, coneSize);CHKERRQ(ierr); 746 ierr = DMPlexGetSupportSize(dm, oldp, &supportSize);CHKERRQ(ierr); 747 ierr = DMPlexSetSupportSize(sdm, splitp, supportSize);CHKERRQ(ierr); 748 if (dep == depth-1) { 749 const PetscInt hybcell = p + pMaxNew[dep+1] + numSplitPoints[dep+1]; 750 751 /* Add cohesive cells, they are prisms */ 752 ierr = DMPlexSetConeSize(sdm, hybcell, 2 + coneSize);CHKERRQ(ierr); 753 } else if (dep == 0) { 754 const PetscInt hybedge = p + pMaxNew[dep+1] + numSplitPoints[dep+1]; 755 756 ierr = DMPlexGetSupport(dm, oldp, &support);CHKERRQ(ierr); 757 for (e = 0, qn = 0, qp = 0, qf = 0; e < supportSize; ++e) { 758 PetscInt val; 759 760 ierr = DMLabelGetValue(label, support[e], &val);CHKERRQ(ierr); 761 if (val == 1) ++qf; 762 if ((val == 1) || (val == (shift + 1))) ++qn; 763 if ((val == 1) || (val == -(shift + 1))) ++qp; 764 } 765 /* Split old vertex: Edges into original vertex and new cohesive edge */ 766 ierr = DMPlexSetSupportSize(sdm, newp, qn+1);CHKERRQ(ierr); 767 /* Split new vertex: Edges into split vertex and new cohesive edge */ 768 ierr = DMPlexSetSupportSize(sdm, splitp, qp+1);CHKERRQ(ierr); 769 /* Add hybrid edge */ 770 ierr = DMPlexSetConeSize(sdm, hybedge, 2);CHKERRQ(ierr); 771 ierr = DMPlexSetSupportSize(sdm, hybedge, qf);CHKERRQ(ierr); 772 } else if (dep == dim-2) { 773 const PetscInt hybface = p + pMaxNew[dep+1] + numSplitPoints[dep+1]; 774 775 ierr = DMPlexGetSupport(dm, oldp, &support);CHKERRQ(ierr); 776 for (e = 0, qn = 0, qp = 0, qf = 0; e < supportSize; ++e) { 777 PetscInt val; 778 779 ierr = DMLabelGetValue(label, support[e], &val);CHKERRQ(ierr); 780 if (val == dim-1) ++qf; 781 if ((val == dim-1) || (val == (shift + dim-1))) ++qn; 782 if ((val == dim-1) || (val == -(shift + dim-1))) ++qp; 783 } 784 /* Split old edge: Faces into original edge and cohesive face (positive side?) */ 785 ierr = DMPlexSetSupportSize(sdm, newp, qn+1);CHKERRQ(ierr); 786 /* Split new edge: Faces into split edge and cohesive face (negative side?) */ 787 ierr = DMPlexSetSupportSize(sdm, splitp, qp+1);CHKERRQ(ierr); 788 /* Add hybrid face */ 789 ierr = DMPlexSetConeSize(sdm, hybface, 4);CHKERRQ(ierr); 790 ierr = DMPlexSetSupportSize(sdm, hybface, qf);CHKERRQ(ierr); 791 } 792 } 793 } 794 for (dep = 0; dep <= depth; ++dep) { 795 for (p = 0; p < numUnsplitPoints[dep]; ++p) { 796 const PetscInt oldp = unsplitPoints[dep][p]; 797 const PetscInt newp = DMPlexShiftPoint_Internal(oldp, depth, depthMax, depthEnd, depthShift) /*oldp + depthOffset[dep]*/; 798 const PetscInt *support; 799 PetscInt coneSize, supportSize, qf, e, s; 800 801 ierr = DMPlexGetConeSize(dm, oldp, &coneSize);CHKERRQ(ierr); 802 ierr = DMPlexGetSupportSize(dm, oldp, &supportSize);CHKERRQ(ierr); 803 ierr = DMPlexGetSupport(dm, oldp, &support);CHKERRQ(ierr); 804 if (dep == 0) { 805 const PetscInt hybedge = p + pMaxNew[dep+1] + numSplitPoints[dep+1] + numSplitPoints[dep]; 806 807 /* Unsplit vertex: Edges into original vertex, split edges, and new cohesive edge twice */ 808 for (s = 0, qf = 0; s < supportSize; ++s, ++qf) { 809 ierr = PetscFindInt(support[s], numSplitPoints[dep+1], splitPoints[dep+1], &e);CHKERRQ(ierr); 810 if (e >= 0) ++qf; 811 } 812 ierr = DMPlexSetSupportSize(sdm, newp, qf+2);CHKERRQ(ierr); 813 /* Add hybrid edge */ 814 ierr = DMPlexSetConeSize(sdm, hybedge, 2);CHKERRQ(ierr); 815 for (e = 0, qf = 0; e < supportSize; ++e) { 816 PetscInt val; 817 818 ierr = DMLabelGetValue(label, support[e], &val);CHKERRQ(ierr); 819 /* Split and unsplit edges produce hybrid faces */ 820 if (val == 1) ++qf; 821 if (val == (shift2 + 1)) ++qf; 822 } 823 ierr = DMPlexSetSupportSize(sdm, hybedge, qf);CHKERRQ(ierr); 824 } else if (dep == dim-2) { 825 const PetscInt hybface = p + pMaxNew[dep+1] + numSplitPoints[dep+1] + numSplitPoints[dep]; 826 PetscInt val; 827 828 for (e = 0, qf = 0; e < supportSize; ++e) { 829 ierr = DMLabelGetValue(label, support[e], &val);CHKERRQ(ierr); 830 if (val == dim-1) qf += 2; 831 else ++qf; 832 } 833 /* Unsplit edge: Faces into original edge, split face, and cohesive face twice */ 834 ierr = DMPlexSetSupportSize(sdm, newp, qf+2);CHKERRQ(ierr); 835 /* Add hybrid face */ 836 for (e = 0, qf = 0; e < supportSize; ++e) { 837 ierr = DMLabelGetValue(label, support[e], &val);CHKERRQ(ierr); 838 if (val == dim-1) ++qf; 839 } 840 ierr = DMPlexSetConeSize(sdm, hybface, 4);CHKERRQ(ierr); 841 ierr = DMPlexSetSupportSize(sdm, hybface, qf);CHKERRQ(ierr); 842 } 843 } 844 } 845 /* Step 4: Setup split DM */ 846 ierr = DMSetUp(sdm);CHKERRQ(ierr); 847 ierr = DMPlexShiftPoints_Internal(dm, depthShift, sdm);CHKERRQ(ierr); 848 ierr = DMPlexGetMaxSizes(sdm, &maxConeSizeNew, &maxSupportSizeNew);CHKERRQ(ierr); 849 ierr = PetscMalloc3(PetscMax(maxConeSize, maxConeSizeNew)*3,&coneNew,PetscMax(maxConeSize, maxConeSizeNew)*3,&coneONew,PetscMax(maxSupportSize, maxSupportSizeNew),&supportNew);CHKERRQ(ierr); 850 /* Step 6: Set cones and supports for new points */ 851 for (dep = 0; dep <= depth; ++dep) { 852 for (p = 0; p < numSplitPoints[dep]; ++p) { 853 const PetscInt oldp = splitPoints[dep][p]; 854 const PetscInt newp = DMPlexShiftPoint_Internal(oldp, depth, depthMax, depthEnd, depthShift) /*oldp + depthOffset[dep]*/; 855 const PetscInt splitp = p + pMaxNew[dep]; 856 const PetscInt *cone, *support, *ornt; 857 PetscInt coneSize, supportSize, q, qf, qn, qp, v, e, s; 858 859 ierr = DMPlexGetConeSize(dm, oldp, &coneSize);CHKERRQ(ierr); 860 ierr = DMPlexGetCone(dm, oldp, &cone);CHKERRQ(ierr); 861 ierr = DMPlexGetConeOrientation(dm, oldp, &ornt);CHKERRQ(ierr); 862 ierr = DMPlexGetSupportSize(dm, oldp, &supportSize);CHKERRQ(ierr); 863 ierr = DMPlexGetSupport(dm, oldp, &support);CHKERRQ(ierr); 864 if (dep == depth-1) { 865 PetscBool hasUnsplit = PETSC_FALSE; 866 const PetscInt hybcell = p + pMaxNew[dep+1] + numSplitPoints[dep+1]; 867 const PetscInt *supportF; 868 869 /* Split face: copy in old face to new face to start */ 870 ierr = DMPlexGetSupport(sdm, newp, &supportF);CHKERRQ(ierr); 871 ierr = DMPlexSetSupport(sdm, splitp, supportF);CHKERRQ(ierr); 872 /* Split old face: old vertices/edges in cone so no change */ 873 /* Split new face: new vertices/edges in cone */ 874 for (q = 0; q < coneSize; ++q) { 875 ierr = PetscFindInt(cone[q], numSplitPoints[dep-1], splitPoints[dep-1], &v);CHKERRQ(ierr); 876 if (v < 0) { 877 ierr = PetscFindInt(cone[q], numUnsplitPoints[dep-1], unsplitPoints[dep-1], &v);CHKERRQ(ierr); 878 if (v < 0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Could not locate point %d in split or unsplit points of depth %d", cone[q], dep-1); 879 coneNew[2+q] = DMPlexShiftPoint_Internal(cone[q], depth, depthMax, depthEnd, depthShift) /*cone[q] + depthOffset[dep-1]*/; 880 hasUnsplit = PETSC_TRUE; 881 } else { 882 coneNew[2+q] = v + pMaxNew[dep-1]; 883 if (dep > 1) { 884 const PetscInt *econe; 885 PetscInt econeSize, r, vs, vu; 886 887 ierr = DMPlexGetConeSize(dm, cone[q], &econeSize);CHKERRQ(ierr); 888 ierr = DMPlexGetCone(dm, cone[q], &econe);CHKERRQ(ierr); 889 for (r = 0; r < econeSize; ++r) { 890 ierr = PetscFindInt(econe[r], numSplitPoints[dep-2], splitPoints[dep-2], &vs);CHKERRQ(ierr); 891 ierr = PetscFindInt(econe[r], numUnsplitPoints[dep-2], unsplitPoints[dep-2], &vu);CHKERRQ(ierr); 892 if (vs >= 0) continue; 893 if (vu < 0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Could not locate point %d in split or unsplit points of depth %d", econe[r], dep-2); 894 hasUnsplit = PETSC_TRUE; 895 } 896 } 897 } 898 } 899 ierr = DMPlexSetCone(sdm, splitp, &coneNew[2]);CHKERRQ(ierr); 900 ierr = DMPlexSetConeOrientation(sdm, splitp, ornt);CHKERRQ(ierr); 901 /* Face support */ 902 for (s = 0; s < supportSize; ++s) { 903 PetscInt val; 904 905 ierr = DMLabelGetValue(label, support[s], &val);CHKERRQ(ierr); 906 if (val < 0) { 907 /* Split old face: Replace negative side cell with cohesive cell */ 908 ierr = DMPlexInsertSupport(sdm, newp, s, hybcell);CHKERRQ(ierr); 909 } else { 910 /* Split new face: Replace positive side cell with cohesive cell */ 911 ierr = DMPlexInsertSupport(sdm, splitp, s, hybcell);CHKERRQ(ierr); 912 /* Get orientation for cohesive face */ 913 { 914 const PetscInt *ncone, *nconeO; 915 PetscInt nconeSize, nc; 916 917 ierr = DMPlexGetConeSize(dm, support[s], &nconeSize);CHKERRQ(ierr); 918 ierr = DMPlexGetCone(dm, support[s], &ncone);CHKERRQ(ierr); 919 ierr = DMPlexGetConeOrientation(dm, support[s], &nconeO);CHKERRQ(ierr); 920 for (nc = 0; nc < nconeSize; ++nc) { 921 if (ncone[nc] == oldp) { 922 coneONew[0] = nconeO[nc]; 923 break; 924 } 925 } 926 if (nc >= nconeSize) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Could not locate face %d in neighboring cell %d", oldp, support[s]); 927 } 928 } 929 } 930 /* Cohesive cell: Old and new split face, then new cohesive faces */ 931 coneNew[0] = newp; /* Extracted negative side orientation above */ 932 coneNew[1] = splitp; 933 coneONew[1] = coneONew[0]; 934 for (q = 0; q < coneSize; ++q) { 935 ierr = PetscFindInt(cone[q], numSplitPoints[dep-1], splitPoints[dep-1], &v);CHKERRQ(ierr); 936 if (v < 0) { 937 ierr = PetscFindInt(cone[q], numUnsplitPoints[dep-1], unsplitPoints[dep-1], &v);CHKERRQ(ierr); 938 coneNew[2+q] = v + pMaxNew[dep] + numSplitPoints[dep] + numSplitPoints[dep-1]; 939 coneONew[2+q] = 0; 940 } else { 941 coneNew[2+q] = v + pMaxNew[dep] + numSplitPoints[dep]; 942 } 943 coneONew[2+q] = 0; 944 } 945 ierr = DMPlexSetCone(sdm, hybcell, coneNew);CHKERRQ(ierr); 946 ierr = DMPlexSetConeOrientation(sdm, hybcell, coneONew);CHKERRQ(ierr); 947 /* Label the hybrid cells on the boundary of the split */ 948 if (hasUnsplit) {ierr = DMLabelSetValue(label, -hybcell, dim);CHKERRQ(ierr);} 949 } else if (dep == 0) { 950 const PetscInt hybedge = p + pMaxNew[dep+1] + numSplitPoints[dep+1]; 951 952 /* Split old vertex: Edges in old split faces and new cohesive edge */ 953 for (e = 0, qn = 0; e < supportSize; ++e) { 954 PetscInt val; 955 956 ierr = DMLabelGetValue(label, support[e], &val);CHKERRQ(ierr); 957 if ((val == 1) || (val == (shift + 1))) { 958 supportNew[qn++] = DMPlexShiftPoint_Internal(support[e], depth, depthMax, depthEnd, depthShift) /*support[e] + depthOffset[dep+1]*/; 959 } 960 } 961 supportNew[qn] = hybedge; 962 ierr = DMPlexSetSupport(sdm, newp, supportNew);CHKERRQ(ierr); 963 /* Split new vertex: Edges in new split faces and new cohesive edge */ 964 for (e = 0, qp = 0; e < supportSize; ++e) { 965 PetscInt val, edge; 966 967 ierr = DMLabelGetValue(label, support[e], &val);CHKERRQ(ierr); 968 if (val == 1) { 969 ierr = PetscFindInt(support[e], numSplitPoints[dep+1], splitPoints[dep+1], &edge);CHKERRQ(ierr); 970 if (edge < 0) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Edge %d is not a split edge", support[e]); 971 supportNew[qp++] = edge + pMaxNew[dep+1]; 972 } else if (val == -(shift + 1)) { 973 supportNew[qp++] = DMPlexShiftPoint_Internal(support[e], depth, depthMax, depthEnd, depthShift) /*support[e] + depthOffset[dep+1]*/; 974 } 975 } 976 supportNew[qp] = hybedge; 977 ierr = DMPlexSetSupport(sdm, splitp, supportNew);CHKERRQ(ierr); 978 /* Hybrid edge: Old and new split vertex */ 979 coneNew[0] = newp; 980 coneNew[1] = splitp; 981 ierr = DMPlexSetCone(sdm, hybedge, coneNew);CHKERRQ(ierr); 982 for (e = 0, qf = 0; e < supportSize; ++e) { 983 PetscInt val, edge; 984 985 ierr = DMLabelGetValue(label, support[e], &val);CHKERRQ(ierr); 986 if (val == 1) { 987 ierr = PetscFindInt(support[e], numSplitPoints[dep+1], splitPoints[dep+1], &edge);CHKERRQ(ierr); 988 if (edge < 0) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Edge %d is not a split edge", support[e]); 989 supportNew[qf++] = edge + pMaxNew[dep+2] + numSplitPoints[dep+2]; 990 } 991 } 992 ierr = DMPlexSetSupport(sdm, hybedge, supportNew);CHKERRQ(ierr); 993 } else if (dep == dim-2) { 994 const PetscInt hybface = p + pMaxNew[dep+1] + numSplitPoints[dep+1]; 995 996 /* Split old edge: old vertices in cone so no change */ 997 /* Split new edge: new vertices in cone */ 998 for (q = 0; q < coneSize; ++q) { 999 ierr = PetscFindInt(cone[q], numSplitPoints[dep-1], splitPoints[dep-1], &v);CHKERRQ(ierr); 1000 if (v < 0) { 1001 ierr = PetscFindInt(cone[q], numUnsplitPoints[dep-1], unsplitPoints[dep-1], &v);CHKERRQ(ierr); 1002 if (v < 0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Could not locate point %d in split or unsplit points of depth %d", cone[q], dep-1); 1003 coneNew[q] = DMPlexShiftPoint_Internal(cone[q], depth, depthMax, depthEnd, depthShift) /*cone[q] + depthOffset[dep-1]*/; 1004 } else { 1005 coneNew[q] = v + pMaxNew[dep-1]; 1006 } 1007 } 1008 ierr = DMPlexSetCone(sdm, splitp, coneNew);CHKERRQ(ierr); 1009 /* Split old edge: Faces in positive side cells and old split faces */ 1010 for (e = 0, q = 0; e < supportSize; ++e) { 1011 PetscInt val; 1012 1013 ierr = DMLabelGetValue(label, support[e], &val);CHKERRQ(ierr); 1014 if (val == dim-1) { 1015 supportNew[q++] = DMPlexShiftPoint_Internal(support[e], depth, depthMax, depthEnd, depthShift) /*support[e] + depthOffset[dep+1]*/; 1016 } else if (val == (shift + dim-1)) { 1017 supportNew[q++] = DMPlexShiftPoint_Internal(support[e], depth, depthMax, depthEnd, depthShift) /*support[e] + depthOffset[dep+1]*/; 1018 } 1019 } 1020 supportNew[q++] = p + pMaxNew[dep+1] + numSplitPoints[dep+1]; 1021 ierr = DMPlexSetSupport(sdm, newp, supportNew);CHKERRQ(ierr); 1022 /* Split new edge: Faces in negative side cells and new split faces */ 1023 for (e = 0, q = 0; e < supportSize; ++e) { 1024 PetscInt val, face; 1025 1026 ierr = DMLabelGetValue(label, support[e], &val);CHKERRQ(ierr); 1027 if (val == dim-1) { 1028 ierr = PetscFindInt(support[e], numSplitPoints[dep+1], splitPoints[dep+1], &face);CHKERRQ(ierr); 1029 if (face < 0) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Face %d is not a split face", support[e]); 1030 supportNew[q++] = face + pMaxNew[dep+1]; 1031 } else if (val == -(shift + dim-1)) { 1032 supportNew[q++] = DMPlexShiftPoint_Internal(support[e], depth, depthMax, depthEnd, depthShift) /*support[e] + depthOffset[dep+1]*/; 1033 } 1034 } 1035 supportNew[q++] = p + pMaxNew[dep+1] + numSplitPoints[dep+1]; 1036 ierr = DMPlexSetSupport(sdm, splitp, supportNew);CHKERRQ(ierr); 1037 /* Hybrid face */ 1038 coneNew[0] = newp; 1039 coneNew[1] = splitp; 1040 for (v = 0; v < coneSize; ++v) { 1041 PetscInt vertex; 1042 ierr = PetscFindInt(cone[v], numSplitPoints[dep-1], splitPoints[dep-1], &vertex);CHKERRQ(ierr); 1043 if (vertex < 0) { 1044 ierr = PetscFindInt(cone[v], numUnsplitPoints[dep-1], unsplitPoints[dep-1], &vertex);CHKERRQ(ierr); 1045 if (vertex < 0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Could not locate point %d in split or unsplit points of depth %d", cone[v], dep-1); 1046 coneNew[2+v] = vertex + pMaxNew[dep] + numSplitPoints[dep] + numSplitPoints[dep-1]; 1047 } else { 1048 coneNew[2+v] = vertex + pMaxNew[dep] + numSplitPoints[dep]; 1049 } 1050 } 1051 ierr = DMPlexSetCone(sdm, hybface, coneNew);CHKERRQ(ierr); 1052 for (e = 0, qf = 0; e < supportSize; ++e) { 1053 PetscInt val, face; 1054 1055 ierr = DMLabelGetValue(label, support[e], &val);CHKERRQ(ierr); 1056 if (val == dim-1) { 1057 ierr = PetscFindInt(support[e], numSplitPoints[dep+1], splitPoints[dep+1], &face);CHKERRQ(ierr); 1058 if (face < 0) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Face %d is not a split face", support[e]); 1059 supportNew[qf++] = face + pMaxNew[dep+2] + numSplitPoints[dep+2]; 1060 } 1061 } 1062 ierr = DMPlexSetSupport(sdm, hybface, supportNew);CHKERRQ(ierr); 1063 } 1064 } 1065 } 1066 for (dep = 0; dep <= depth; ++dep) { 1067 for (p = 0; p < numUnsplitPoints[dep]; ++p) { 1068 const PetscInt oldp = unsplitPoints[dep][p]; 1069 const PetscInt newp = DMPlexShiftPoint_Internal(oldp, depth, depthMax, depthEnd, depthShift) /*oldp + depthOffset[dep]*/; 1070 const PetscInt *cone, *support, *ornt; 1071 PetscInt coneSize, supportSize, supportSizeNew, q, qf, e, f, s; 1072 1073 ierr = DMPlexGetConeSize(dm, oldp, &coneSize);CHKERRQ(ierr); 1074 ierr = DMPlexGetCone(dm, oldp, &cone);CHKERRQ(ierr); 1075 ierr = DMPlexGetConeOrientation(dm, oldp, &ornt);CHKERRQ(ierr); 1076 ierr = DMPlexGetSupportSize(dm, oldp, &supportSize);CHKERRQ(ierr); 1077 ierr = DMPlexGetSupport(dm, oldp, &support);CHKERRQ(ierr); 1078 if (dep == 0) { 1079 const PetscInt hybedge = p + pMaxNew[dep+1] + numSplitPoints[dep+1] + numSplitPoints[dep]; 1080 1081 /* Unsplit vertex */ 1082 ierr = DMPlexGetSupportSize(sdm, newp, &supportSizeNew);CHKERRQ(ierr); 1083 for (s = 0, q = 0; s < supportSize; ++s) { 1084 supportNew[q++] = DMPlexShiftPoint_Internal(support[s], depth, depthMax, depthEnd, depthShift) /*support[s] + depthOffset[dep+1]*/; 1085 ierr = PetscFindInt(support[s], numSplitPoints[dep+1], splitPoints[dep+1], &e);CHKERRQ(ierr); 1086 if (e >= 0) { 1087 supportNew[q++] = e + pMaxNew[dep+1]; 1088 } 1089 } 1090 supportNew[q++] = hybedge; 1091 supportNew[q++] = hybedge; 1092 if (q != supportSizeNew) SETERRQ3(comm, PETSC_ERR_ARG_WRONG, "Support size %d != %d for vertex %d", q, supportSizeNew, newp); 1093 ierr = DMPlexSetSupport(sdm, newp, supportNew);CHKERRQ(ierr); 1094 /* Hybrid edge */ 1095 coneNew[0] = newp; 1096 coneNew[1] = newp; 1097 ierr = DMPlexSetCone(sdm, hybedge, coneNew);CHKERRQ(ierr); 1098 for (e = 0, qf = 0; e < supportSize; ++e) { 1099 PetscInt val, edge; 1100 1101 ierr = DMLabelGetValue(label, support[e], &val);CHKERRQ(ierr); 1102 if (val == 1) { 1103 ierr = PetscFindInt(support[e], numSplitPoints[dep+1], splitPoints[dep+1], &edge);CHKERRQ(ierr); 1104 if (edge < 0) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Edge %d is not a split edge", support[e]); 1105 supportNew[qf++] = edge + pMaxNew[dep+2] + numSplitPoints[dep+2]; 1106 } else if (val == (shift2 + 1)) { 1107 ierr = PetscFindInt(support[e], numUnsplitPoints[dep+1], unsplitPoints[dep+1], &edge);CHKERRQ(ierr); 1108 if (edge < 0) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Edge %d is not a unsplit edge", support[e]); 1109 supportNew[qf++] = edge + pMaxNew[dep+2] + numSplitPoints[dep+2] + numSplitPoints[dep+1]; 1110 } 1111 } 1112 ierr = DMPlexSetSupport(sdm, hybedge, supportNew);CHKERRQ(ierr); 1113 } else if (dep == dim-2) { 1114 const PetscInt hybface = p + pMaxNew[dep+1] + numSplitPoints[dep+1] + numSplitPoints[dep]; 1115 1116 /* Unsplit edge: Faces into original edge, split face, and hybrid face twice */ 1117 for (f = 0, qf = 0; f < supportSize; ++f) { 1118 PetscInt val, face; 1119 1120 ierr = DMLabelGetValue(label, support[f], &val);CHKERRQ(ierr); 1121 if (val == dim-1) { 1122 ierr = PetscFindInt(support[f], numSplitPoints[dep+1], splitPoints[dep+1], &face);CHKERRQ(ierr); 1123 if (face < 0) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Face %d is not a split face", support[f]); 1124 supportNew[qf++] = DMPlexShiftPoint_Internal(support[f], depth, depthMax, depthEnd, depthShift) /*support[f] + depthOffset[dep+1]*/; 1125 supportNew[qf++] = face + pMaxNew[dep+1]; 1126 } else { 1127 supportNew[qf++] = DMPlexShiftPoint_Internal(support[f], depth, depthMax, depthEnd, depthShift) /*support[f] + depthOffset[dep+1]*/; 1128 } 1129 } 1130 supportNew[qf++] = hybface; 1131 supportNew[qf++] = hybface; 1132 ierr = DMPlexGetSupportSize(sdm, newp, &supportSizeNew);CHKERRQ(ierr); 1133 if (qf != supportSizeNew) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Support size for unsplit edge %d is %d != %d\n", newp, qf, supportSizeNew); 1134 ierr = DMPlexSetSupport(sdm, newp, supportNew);CHKERRQ(ierr); 1135 /* Add hybrid face */ 1136 coneNew[0] = newp; 1137 coneNew[1] = newp; 1138 ierr = PetscFindInt(cone[0], numUnsplitPoints[dep-1], unsplitPoints[dep-1], &v);CHKERRQ(ierr); 1139 if (v < 0) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Vertex %d is not an unsplit vertex", cone[0]); 1140 coneNew[2] = v + pMaxNew[dep] + numSplitPoints[dep] + numSplitPoints[dep-1]; 1141 ierr = PetscFindInt(cone[1], numUnsplitPoints[dep-1], unsplitPoints[dep-1], &v);CHKERRQ(ierr); 1142 if (v < 0) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Vertex %d is not an unsplit vertex", cone[1]); 1143 coneNew[3] = v + pMaxNew[dep] + numSplitPoints[dep] + numSplitPoints[dep-1]; 1144 ierr = DMPlexSetCone(sdm, hybface, coneNew);CHKERRQ(ierr); 1145 for (f = 0, qf = 0; f < supportSize; ++f) { 1146 PetscInt val, face; 1147 1148 ierr = DMLabelGetValue(label, support[f], &val);CHKERRQ(ierr); 1149 if (val == dim-1) { 1150 ierr = PetscFindInt(support[f], numSplitPoints[dep+1], splitPoints[dep+1], &face);CHKERRQ(ierr); 1151 supportNew[qf++] = face + pMaxNew[dep+2] + numSplitPoints[dep+2]; 1152 } 1153 } 1154 ierr = DMPlexGetSupportSize(sdm, hybface, &supportSizeNew);CHKERRQ(ierr); 1155 if (qf != supportSizeNew) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Support size for hybrid face %d is %d != %d\n", hybface, qf, supportSizeNew); 1156 ierr = DMPlexSetSupport(sdm, hybface, supportNew);CHKERRQ(ierr); 1157 } 1158 } 1159 } 1160 /* Step 6b: Replace split points in negative side cones */ 1161 for (sp = 0; sp < numSP; ++sp) { 1162 PetscInt dep = values[sp]; 1163 IS pIS; 1164 PetscInt numPoints; 1165 const PetscInt *points; 1166 1167 if (dep >= 0) continue; 1168 ierr = DMLabelGetStratumIS(label, dep, &pIS);CHKERRQ(ierr); 1169 if (!pIS) continue; 1170 dep = -dep - shift; 1171 ierr = ISGetLocalSize(pIS, &numPoints);CHKERRQ(ierr); 1172 ierr = ISGetIndices(pIS, &points);CHKERRQ(ierr); 1173 for (p = 0; p < numPoints; ++p) { 1174 const PetscInt oldp = points[p]; 1175 const PetscInt newp = DMPlexShiftPoint_Internal(oldp, depth, depthMax, depthEnd, depthShift) /*depthOffset[dep] + oldp*/; 1176 const PetscInt *cone; 1177 PetscInt coneSize, c; 1178 /* PetscBool replaced = PETSC_FALSE; */ 1179 1180 /* Negative edge: replace split vertex */ 1181 /* Negative cell: replace split face */ 1182 ierr = DMPlexGetConeSize(sdm, newp, &coneSize);CHKERRQ(ierr); 1183 ierr = DMPlexGetCone(sdm, newp, &cone);CHKERRQ(ierr); 1184 for (c = 0; c < coneSize; ++c) { 1185 const PetscInt coldp = cone[c] - depthOffset[dep-1]; 1186 PetscInt csplitp, cp, val; 1187 1188 ierr = DMLabelGetValue(label, coldp, &val);CHKERRQ(ierr); 1189 if (val == dep-1) { 1190 ierr = PetscFindInt(coldp, numSplitPoints[dep-1], splitPoints[dep-1], &cp);CHKERRQ(ierr); 1191 if (cp < 0) SETERRQ2(comm, PETSC_ERR_ARG_WRONG, "Point %d is not a split point of dimension %d", oldp, dep-1); 1192 csplitp = pMaxNew[dep-1] + cp; 1193 ierr = DMPlexInsertCone(sdm, newp, c, csplitp);CHKERRQ(ierr); 1194 /* replaced = PETSC_TRUE; */ 1195 } 1196 } 1197 /* Cells with only a vertex or edge on the submesh have no replacement */ 1198 /* if (!replaced) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "The cone of point %d does not contain split points", oldp); */ 1199 } 1200 ierr = ISRestoreIndices(pIS, &points);CHKERRQ(ierr); 1201 ierr = ISDestroy(&pIS);CHKERRQ(ierr); 1202 } 1203 /* Step 7: Stratify */ 1204 ierr = DMPlexStratify(sdm);CHKERRQ(ierr); 1205 /* Step 8: Coordinates */ 1206 ierr = DMPlexShiftCoordinates_Internal(dm, depthShift, sdm);CHKERRQ(ierr); 1207 ierr = DMGetCoordinateSection(sdm, &coordSection);CHKERRQ(ierr); 1208 ierr = DMGetCoordinatesLocal(sdm, &coordinates);CHKERRQ(ierr); 1209 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 1210 for (v = 0; v < (numSplitPoints ? numSplitPoints[0] : 0); ++v) { 1211 const PetscInt newp = DMPlexShiftPoint_Internal(splitPoints[0][v], depth, depthMax, depthEnd, depthShift) /*depthOffset[0] + splitPoints[0][v]*/; 1212 const PetscInt splitp = pMaxNew[0] + v; 1213 PetscInt dof, off, soff, d; 1214 1215 ierr = PetscSectionGetDof(coordSection, newp, &dof);CHKERRQ(ierr); 1216 ierr = PetscSectionGetOffset(coordSection, newp, &off);CHKERRQ(ierr); 1217 ierr = PetscSectionGetOffset(coordSection, splitp, &soff);CHKERRQ(ierr); 1218 for (d = 0; d < dof; ++d) coords[soff+d] = coords[off+d]; 1219 } 1220 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 1221 /* Step 9: SF, if I can figure this out we can split the mesh in parallel */ 1222 ierr = DMPlexShiftSF_Internal(dm, depthShift, sdm);CHKERRQ(ierr); 1223 /* Step 10: Labels */ 1224 ierr = DMPlexShiftLabels_Internal(dm, depthShift, sdm);CHKERRQ(ierr); 1225 ierr = DMPlexGetNumLabels(sdm, &numLabels);CHKERRQ(ierr); 1226 for (dep = 0; dep <= depth; ++dep) { 1227 for (p = 0; p < numSplitPoints[dep]; ++p) { 1228 const PetscInt newp = DMPlexShiftPoint_Internal(splitPoints[dep][p], depth, depthMax, depthEnd, depthShift) /*depthOffset[dep] + splitPoints[dep][p]*/; 1229 const PetscInt splitp = pMaxNew[dep] + p; 1230 PetscInt l; 1231 1232 for (l = 0; l < numLabels; ++l) { 1233 DMLabel mlabel; 1234 const char *lname; 1235 PetscInt val; 1236 PetscBool isDepth; 1237 1238 ierr = DMPlexGetLabelName(sdm, l, &lname);CHKERRQ(ierr); 1239 ierr = PetscStrcmp(lname, "depth", &isDepth);CHKERRQ(ierr); 1240 if (isDepth) continue; 1241 ierr = DMPlexGetLabel(sdm, lname, &mlabel);CHKERRQ(ierr); 1242 ierr = DMLabelGetValue(mlabel, newp, &val);CHKERRQ(ierr); 1243 if (val >= 0) { 1244 ierr = DMLabelSetValue(mlabel, splitp, val);CHKERRQ(ierr); 1245 #if 0 1246 /* Do not put cohesive edges into the label */ 1247 if (dep == 0) { 1248 const PetscInt cedge = p + pMaxNew[dep+1] + numSplitPoints[dep+1]; 1249 ierr = DMLabelSetValue(mlabel, cedge, val);CHKERRQ(ierr); 1250 } else if (dep == dim-2) { 1251 const PetscInt cface = p + pMaxNew[dep+1] + numSplitPoints[dep+1]; 1252 ierr = DMLabelSetValue(mlabel, cface, val);CHKERRQ(ierr); 1253 } 1254 /* Do not put cohesive faces into the label */ 1255 #endif 1256 } 1257 } 1258 } 1259 } 1260 for (sp = 0; sp < numSP; ++sp) { 1261 const PetscInt dep = values[sp]; 1262 1263 if ((dep < 0) || (dep > depth)) continue; 1264 if (splitIS[dep]) {ierr = ISRestoreIndices(splitIS[dep], &splitPoints[dep]);CHKERRQ(ierr);} 1265 ierr = ISDestroy(&splitIS[dep]);CHKERRQ(ierr); 1266 if (unsplitIS[dep]) {ierr = ISRestoreIndices(unsplitIS[dep], &unsplitPoints[dep]);CHKERRQ(ierr);} 1267 ierr = ISDestroy(&unsplitIS[dep]);CHKERRQ(ierr); 1268 } 1269 if (label) { 1270 ierr = ISRestoreIndices(valueIS, &values);CHKERRQ(ierr); 1271 ierr = ISDestroy(&valueIS);CHKERRQ(ierr); 1272 } 1273 for (d = 0; d <= depth; ++d) { 1274 ierr = DMPlexGetDepthStratum(sdm, d, NULL, &pEnd);CHKERRQ(ierr); 1275 pMaxNew[d] = pEnd - numHybridPoints[d] - numHybridPointsOld[d]; 1276 } 1277 ierr = DMPlexSetHybridBounds(sdm, depth >= 0 ? pMaxNew[depth] : PETSC_DETERMINE, depth>1 ? pMaxNew[depth-1] : PETSC_DETERMINE, depth>2 ? pMaxNew[1] : PETSC_DETERMINE, depth >= 0 ? pMaxNew[0] : PETSC_DETERMINE);CHKERRQ(ierr); 1278 ierr = PetscFree3(coneNew, coneONew, supportNew);CHKERRQ(ierr); 1279 ierr = PetscFree6(depthMax, depthEnd, depthShift, depthOffset, pMaxNew, numHybridPointsOld);CHKERRQ(ierr); 1280 ierr = PetscFree7(splitIS, unsplitIS, numSplitPoints, numUnsplitPoints, numHybridPoints, splitPoints, unsplitPoints);CHKERRQ(ierr); 1281 PetscFunctionReturn(0); 1282 } 1283 1284 #undef __FUNCT__ 1285 #define __FUNCT__ "DMPlexConstructCohesiveCells" 1286 /*@C 1287 DMPlexConstructCohesiveCells - Construct cohesive cells which split the face along an internal interface 1288 1289 Collective on dm 1290 1291 Input Parameters: 1292 + dm - The original DM 1293 - label - The label specifying the boundary faces (this could be auto-generated) 1294 1295 Output Parameters: 1296 - dmSplit - The new DM 1297 1298 Level: developer 1299 1300 .seealso: DMCreate(), DMPlexLabelCohesiveComplete() 1301 @*/ 1302 PetscErrorCode DMPlexConstructCohesiveCells(DM dm, DMLabel label, DM *dmSplit) 1303 { 1304 DM sdm; 1305 PetscInt dim; 1306 PetscErrorCode ierr; 1307 1308 PetscFunctionBegin; 1309 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 1310 PetscValidPointer(dmSplit, 3); 1311 ierr = DMCreate(PetscObjectComm((PetscObject)dm), &sdm);CHKERRQ(ierr); 1312 ierr = DMSetType(sdm, DMPLEX);CHKERRQ(ierr); 1313 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 1314 ierr = DMSetDimension(sdm, dim);CHKERRQ(ierr); 1315 switch (dim) { 1316 case 2: 1317 case 3: 1318 ierr = DMPlexConstructCohesiveCells_Internal(dm, label, sdm);CHKERRQ(ierr); 1319 break; 1320 default: 1321 SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Cannot construct cohesive cells for dimension %d", dim); 1322 } 1323 *dmSplit = sdm; 1324 PetscFunctionReturn(0); 1325 } 1326 1327 #undef __FUNCT__ 1328 #define __FUNCT__ "GetSurfaceSide_Static" 1329 /* Returns the side of the surface for a given cell with a face on the surface */ 1330 static PetscErrorCode GetSurfaceSide_Static(DM dm, DM subdm, PetscInt numSubpoints, const PetscInt *subpoints, PetscInt cell, PetscInt face, PetscBool *pos) 1331 { 1332 const PetscInt *cone, *ornt; 1333 PetscInt dim, coneSize, c; 1334 PetscErrorCode ierr; 1335 1336 PetscFunctionBegin; 1337 *pos = PETSC_TRUE; 1338 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 1339 ierr = DMPlexGetConeSize(dm, cell, &coneSize);CHKERRQ(ierr); 1340 ierr = DMPlexGetCone(dm, cell, &cone);CHKERRQ(ierr); 1341 ierr = DMPlexGetConeOrientation(dm, cell, &ornt);CHKERRQ(ierr); 1342 for (c = 0; c < coneSize; ++c) { 1343 if (cone[c] == face) { 1344 PetscInt o = ornt[c]; 1345 1346 if (subdm) { 1347 const PetscInt *subcone, *subornt; 1348 PetscInt subpoint, subface, subconeSize, sc; 1349 1350 ierr = PetscFindInt(cell, numSubpoints, subpoints, &subpoint);CHKERRQ(ierr); 1351 ierr = PetscFindInt(face, numSubpoints, subpoints, &subface);CHKERRQ(ierr); 1352 ierr = DMPlexGetConeSize(subdm, subpoint, &subconeSize);CHKERRQ(ierr); 1353 ierr = DMPlexGetCone(subdm, subpoint, &subcone);CHKERRQ(ierr); 1354 ierr = DMPlexGetConeOrientation(subdm, subpoint, &subornt);CHKERRQ(ierr); 1355 for (sc = 0; sc < subconeSize; ++sc) { 1356 if (subcone[sc] == subface) { 1357 o = subornt[0]; 1358 break; 1359 } 1360 } 1361 if (sc >= subconeSize) SETERRQ4(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Could not find subpoint %d (%d) in cone for subpoint %d (%d)", subface, face, subpoint, cell); 1362 } 1363 if (o >= 0) *pos = PETSC_TRUE; 1364 else *pos = PETSC_FALSE; 1365 break; 1366 } 1367 } 1368 if (c == coneSize) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Cell %d in split face %d support does not have it in the cone", cell, face); 1369 PetscFunctionReturn(0); 1370 } 1371 1372 #undef __FUNCT__ 1373 #define __FUNCT__ "DMPlexLabelCohesiveComplete" 1374 /*@ 1375 DMPlexLabelCohesiveComplete - Starting with a label marking points on an internal surface, we add all other mesh pieces 1376 to complete the surface 1377 1378 Input Parameters: 1379 + dm - The DM 1380 . label - A DMLabel marking the surface 1381 . blabel - A DMLabel marking the vertices on the boundary which will not be duplicated, or NULL to find them automatically 1382 . flip - Flag to flip the submesh normal and replace points on the other side 1383 - subdm - The subDM associated with the label, or NULL 1384 1385 Output Parameter: 1386 . label - A DMLabel marking all surface points 1387 1388 Note: The vertices in blabel are called "unsplit" in the terminology from hybrid cell creation. 1389 1390 Level: developer 1391 1392 .seealso: DMPlexConstructCohesiveCells(), DMPlexLabelComplete() 1393 @*/ 1394 PetscErrorCode DMPlexLabelCohesiveComplete(DM dm, DMLabel label, DMLabel blabel, PetscBool flip, DM subdm) 1395 { 1396 DMLabel depthLabel; 1397 IS dimIS, subpointIS, facePosIS, faceNegIS, crossEdgeIS = NULL; 1398 const PetscInt *points, *subpoints; 1399 const PetscInt rev = flip ? -1 : 1; 1400 PetscInt *pMax; 1401 PetscInt shift = 100, shift2 = 200, dim, depth, pSize, dep, cStart, cEnd, cMax, fStart, fEnd, vStart, vEnd, numPoints, numSubpoints, p, val; 1402 PetscErrorCode ierr; 1403 1404 PetscFunctionBegin; 1405 ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr); 1406 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 1407 pSize = PetscMax(depth, dim) + 1; 1408 ierr = PetscMalloc1(pSize,&pMax);CHKERRQ(ierr); 1409 ierr = DMPlexGetHybridBounds(dm, depth >= 0 ? &pMax[depth] : NULL, depth>1 ? &pMax[depth-1] : NULL, depth>2 ? &pMax[1] : NULL, &pMax[0]);CHKERRQ(ierr); 1410 ierr = DMPlexGetDepthLabel(dm, &depthLabel);CHKERRQ(ierr); 1411 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 1412 if (subdm) { 1413 ierr = DMPlexCreateSubpointIS(subdm, &subpointIS);CHKERRQ(ierr); 1414 if (subpointIS) { 1415 ierr = ISGetLocalSize(subpointIS, &numSubpoints);CHKERRQ(ierr); 1416 ierr = ISGetIndices(subpointIS, &subpoints);CHKERRQ(ierr); 1417 } 1418 } 1419 /* Mark cell on the fault, and its faces which touch the fault: cell orientation for face gives the side of the fault */ 1420 ierr = DMLabelGetStratumIS(label, dim-1, &dimIS);CHKERRQ(ierr); 1421 if (!dimIS) { 1422 ierr = PetscFree(pMax);CHKERRQ(ierr); 1423 PetscFunctionReturn(0); 1424 } 1425 ierr = ISGetLocalSize(dimIS, &numPoints);CHKERRQ(ierr); 1426 ierr = ISGetIndices(dimIS, &points);CHKERRQ(ierr); 1427 for (p = 0; p < numPoints; ++p) { /* Loop over fault faces */ 1428 const PetscInt *support; 1429 PetscInt supportSize, s; 1430 1431 ierr = DMPlexGetSupportSize(dm, points[p], &supportSize);CHKERRQ(ierr); 1432 if (supportSize != 2) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Split face %d has %d != 2 supports", points[p], supportSize); 1433 ierr = DMPlexGetSupport(dm, points[p], &support);CHKERRQ(ierr); 1434 for (s = 0; s < supportSize; ++s) { 1435 const PetscInt *cone; 1436 PetscInt coneSize, c; 1437 PetscBool pos; 1438 1439 ierr = GetSurfaceSide_Static(dm, subdm, numSubpoints, subpoints, support[s], points[p], &pos);CHKERRQ(ierr); 1440 if (pos) {ierr = DMLabelSetValue(label, support[s], rev*(shift+dim));CHKERRQ(ierr);} 1441 else {ierr = DMLabelSetValue(label, support[s], -rev*(shift+dim));CHKERRQ(ierr);} 1442 if (rev < 0) pos = !pos ? PETSC_TRUE : PETSC_FALSE; 1443 /* Put faces touching the fault in the label */ 1444 ierr = DMPlexGetConeSize(dm, support[s], &coneSize);CHKERRQ(ierr); 1445 ierr = DMPlexGetCone(dm, support[s], &cone);CHKERRQ(ierr); 1446 for (c = 0; c < coneSize; ++c) { 1447 const PetscInt point = cone[c]; 1448 1449 ierr = DMLabelGetValue(label, point, &val);CHKERRQ(ierr); 1450 if (val == -1) { 1451 PetscInt *closure = NULL; 1452 PetscInt closureSize, cl; 1453 1454 ierr = DMPlexGetTransitiveClosure(dm, point, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 1455 for (cl = 0; cl < closureSize*2; cl += 2) { 1456 const PetscInt clp = closure[cl]; 1457 PetscInt bval = -1; 1458 1459 ierr = DMLabelGetValue(label, clp, &val);CHKERRQ(ierr); 1460 if (blabel) {ierr = DMLabelGetValue(blabel, clp, &bval);CHKERRQ(ierr);} 1461 if ((val >= 0) && (val < dim-1) && (bval < 0)) { 1462 ierr = DMLabelSetValue(label, point, pos == PETSC_TRUE ? shift+dim-1 : -(shift+dim-1));CHKERRQ(ierr); 1463 break; 1464 } 1465 } 1466 ierr = DMPlexRestoreTransitiveClosure(dm, point, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 1467 } 1468 } 1469 } 1470 } 1471 ierr = ISRestoreIndices(dimIS, &points);CHKERRQ(ierr); 1472 ierr = ISDestroy(&dimIS);CHKERRQ(ierr); 1473 if (subdm) { 1474 if (subpointIS) {ierr = ISRestoreIndices(subpointIS, &subpoints);CHKERRQ(ierr);} 1475 ierr = ISDestroy(&subpointIS);CHKERRQ(ierr); 1476 } 1477 /* Mark boundary points as unsplit */ 1478 if (blabel) { 1479 ierr = DMLabelGetStratumIS(blabel, 1, &dimIS);CHKERRQ(ierr); 1480 ierr = ISGetLocalSize(dimIS, &numPoints);CHKERRQ(ierr); 1481 ierr = ISGetIndices(dimIS, &points);CHKERRQ(ierr); 1482 for (p = 0; p < numPoints; ++p) { 1483 const PetscInt point = points[p]; 1484 PetscInt val, bval; 1485 1486 ierr = DMLabelGetValue(blabel, point, &bval);CHKERRQ(ierr); 1487 if (bval >= 0) { 1488 ierr = DMLabelGetValue(label, point, &val);CHKERRQ(ierr); 1489 if ((val < 0) || (val > dim)) { 1490 /* This could be a point added from splitting a vertex on an adjacent fault, otherwise its just wrong */ 1491 ierr = DMLabelClearValue(blabel, point, bval);CHKERRQ(ierr); 1492 } 1493 } 1494 } 1495 for (p = 0; p < numPoints; ++p) { 1496 const PetscInt point = points[p]; 1497 PetscInt val, bval; 1498 1499 ierr = DMLabelGetValue(blabel, point, &bval);CHKERRQ(ierr); 1500 if (bval >= 0) { 1501 const PetscInt *cone, *support; 1502 PetscInt coneSize, supportSize, s, valA, valB, valE; 1503 1504 /* Mark as unsplit */ 1505 ierr = DMLabelGetValue(label, point, &val);CHKERRQ(ierr); 1506 if ((val < 0) || (val > dim)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Point %d has label value %d, should be part of the fault", point, val); 1507 ierr = DMLabelClearValue(label, point, val);CHKERRQ(ierr); 1508 ierr = DMLabelSetValue(label, point, shift2+val);CHKERRQ(ierr); 1509 /* Check for cross-edge 1510 A cross-edge has endpoints which are both on the boundary of the surface, but the edge itself is not. */ 1511 if (val != 0) continue; 1512 ierr = DMPlexGetSupport(dm, point, &support);CHKERRQ(ierr); 1513 ierr = DMPlexGetSupportSize(dm, point, &supportSize);CHKERRQ(ierr); 1514 for (s = 0; s < supportSize; ++s) { 1515 ierr = DMPlexGetCone(dm, support[s], &cone);CHKERRQ(ierr); 1516 ierr = DMPlexGetConeSize(dm, support[s], &coneSize);CHKERRQ(ierr); 1517 if (coneSize != 2) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Edge %D has %D vertices != 2", support[s], coneSize); 1518 ierr = DMLabelGetValue(blabel, cone[0], &valA);CHKERRQ(ierr); 1519 ierr = DMLabelGetValue(blabel, cone[1], &valB);CHKERRQ(ierr); 1520 ierr = DMLabelGetValue(blabel, support[s], &valE);CHKERRQ(ierr); 1521 if ((valE < 0) && (valA >= 0) && (valB >= 0) && (cone[0] != cone[1])) {ierr = DMLabelSetValue(blabel, support[s], 2);CHKERRQ(ierr);} 1522 } 1523 } 1524 } 1525 ierr = ISRestoreIndices(dimIS, &points);CHKERRQ(ierr); 1526 ierr = ISDestroy(&dimIS);CHKERRQ(ierr); 1527 } 1528 /* Search for other cells/faces/edges connected to the fault by a vertex */ 1529 ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr); 1530 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 1531 ierr = DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);CHKERRQ(ierr); 1532 cMax = cMax < 0 ? cEnd : cMax; 1533 ierr = DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);CHKERRQ(ierr); 1534 ierr = DMLabelGetStratumIS(label, 0, &dimIS);CHKERRQ(ierr); 1535 if (blabel) {ierr = DMLabelGetStratumIS(blabel, 2, &crossEdgeIS);CHKERRQ(ierr);} 1536 if (dimIS && crossEdgeIS) { 1537 IS vertIS = dimIS; 1538 1539 ierr = ISExpand(vertIS, crossEdgeIS, &dimIS);CHKERRQ(ierr); 1540 ierr = ISDestroy(&crossEdgeIS);CHKERRQ(ierr); 1541 ierr = ISDestroy(&vertIS);CHKERRQ(ierr); 1542 } 1543 if (!dimIS) { 1544 ierr = PetscFree(pMax);CHKERRQ(ierr); 1545 PetscFunctionReturn(0); 1546 } 1547 ierr = ISGetLocalSize(dimIS, &numPoints);CHKERRQ(ierr); 1548 ierr = ISGetIndices(dimIS, &points);CHKERRQ(ierr); 1549 for (p = 0; p < numPoints; ++p) { /* Loop over fault vertices */ 1550 PetscInt *star = NULL; 1551 PetscInt starSize, s; 1552 PetscInt again = 1; /* 0: Finished 1: Keep iterating after a change 2: No change */ 1553 1554 /* All points connected to the fault are inside a cell, so at the top level we will only check cells */ 1555 ierr = DMPlexGetTransitiveClosure(dm, points[p], PETSC_FALSE, &starSize, &star);CHKERRQ(ierr); 1556 while (again) { 1557 if (again > 1) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Could not classify all cells connected to the fault"); 1558 again = 0; 1559 for (s = 0; s < starSize*2; s += 2) { 1560 const PetscInt point = star[s]; 1561 const PetscInt *cone; 1562 PetscInt coneSize, c; 1563 1564 if ((point < cStart) || (point >= cMax)) continue; 1565 ierr = DMLabelGetValue(label, point, &val);CHKERRQ(ierr); 1566 if (val != -1) continue; 1567 again = again == 1 ? 1 : 2; 1568 ierr = DMPlexGetConeSize(dm, point, &coneSize);CHKERRQ(ierr); 1569 ierr = DMPlexGetCone(dm, point, &cone);CHKERRQ(ierr); 1570 for (c = 0; c < coneSize; ++c) { 1571 ierr = DMLabelGetValue(label, cone[c], &val);CHKERRQ(ierr); 1572 if (val != -1) { 1573 const PetscInt *ccone; 1574 PetscInt cconeSize, cc, side; 1575 1576 if (PetscAbs(val) < shift) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Face %d on cell %d has an invalid label %d", cone[c], point, val); 1577 if (val > 0) side = 1; 1578 else side = -1; 1579 ierr = DMLabelSetValue(label, point, side*(shift+dim));CHKERRQ(ierr); 1580 /* Mark cell faces which touch the fault */ 1581 ierr = DMPlexGetConeSize(dm, point, &cconeSize);CHKERRQ(ierr); 1582 ierr = DMPlexGetCone(dm, point, &ccone);CHKERRQ(ierr); 1583 for (cc = 0; cc < cconeSize; ++cc) { 1584 PetscInt *closure = NULL; 1585 PetscInt closureSize, cl; 1586 1587 ierr = DMLabelGetValue(label, ccone[cc], &val);CHKERRQ(ierr); 1588 if (val != -1) continue; 1589 ierr = DMPlexGetTransitiveClosure(dm, ccone[cc], PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 1590 for (cl = 0; cl < closureSize*2; cl += 2) { 1591 const PetscInt clp = closure[cl]; 1592 1593 ierr = DMLabelGetValue(label, clp, &val);CHKERRQ(ierr); 1594 if (val == -1) continue; 1595 ierr = DMLabelSetValue(label, ccone[cc], side*(shift+dim-1));CHKERRQ(ierr); 1596 break; 1597 } 1598 ierr = DMPlexRestoreTransitiveClosure(dm, ccone[cc], PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 1599 } 1600 again = 1; 1601 break; 1602 } 1603 } 1604 } 1605 } 1606 /* Classify the rest by cell membership */ 1607 for (s = 0; s < starSize*2; s += 2) { 1608 const PetscInt point = star[s]; 1609 1610 ierr = DMLabelGetValue(label, point, &val);CHKERRQ(ierr); 1611 if (val == -1) { 1612 PetscInt *sstar = NULL; 1613 PetscInt sstarSize, ss; 1614 PetscBool marked = PETSC_FALSE; 1615 1616 ierr = DMPlexGetTransitiveClosure(dm, point, PETSC_FALSE, &sstarSize, &sstar);CHKERRQ(ierr); 1617 for (ss = 0; ss < sstarSize*2; ss += 2) { 1618 const PetscInt spoint = sstar[ss]; 1619 1620 if ((spoint < cStart) || (spoint >= cMax)) continue; 1621 ierr = DMLabelGetValue(label, spoint, &val);CHKERRQ(ierr); 1622 if (val == -1) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Cell %d in star of %d does not have a valid label", spoint, point); 1623 ierr = DMLabelGetValue(depthLabel, point, &dep);CHKERRQ(ierr); 1624 if (val > 0) { 1625 ierr = DMLabelSetValue(label, point, shift+dep);CHKERRQ(ierr); 1626 } else { 1627 ierr = DMLabelSetValue(label, point, -(shift+dep));CHKERRQ(ierr); 1628 } 1629 marked = PETSC_TRUE; 1630 break; 1631 } 1632 ierr = DMPlexRestoreTransitiveClosure(dm, point, PETSC_FALSE, &sstarSize, &sstar);CHKERRQ(ierr); 1633 ierr = DMLabelGetValue(depthLabel, point, &dep);CHKERRQ(ierr); 1634 if (point < pMax[dep] && !marked) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d could not be classified", point); 1635 } 1636 } 1637 ierr = DMPlexRestoreTransitiveClosure(dm, points[p], PETSC_FALSE, &starSize, &star);CHKERRQ(ierr); 1638 } 1639 ierr = ISRestoreIndices(dimIS, &points);CHKERRQ(ierr); 1640 ierr = ISDestroy(&dimIS);CHKERRQ(ierr); 1641 /* If any faces touching the fault divide cells on either side, split them */ 1642 ierr = DMLabelGetStratumIS(label, shift+dim-1, &facePosIS);CHKERRQ(ierr); 1643 ierr = DMLabelGetStratumIS(label, -(shift+dim-1), &faceNegIS);CHKERRQ(ierr); 1644 ierr = ISExpand(facePosIS, faceNegIS, &dimIS);CHKERRQ(ierr); 1645 ierr = ISDestroy(&facePosIS);CHKERRQ(ierr); 1646 ierr = ISDestroy(&faceNegIS);CHKERRQ(ierr); 1647 ierr = ISGetLocalSize(dimIS, &numPoints);CHKERRQ(ierr); 1648 ierr = ISGetIndices(dimIS, &points);CHKERRQ(ierr); 1649 for (p = 0; p < numPoints; ++p) { 1650 const PetscInt point = points[p]; 1651 const PetscInt *support; 1652 PetscInt supportSize, valA, valB; 1653 1654 ierr = DMPlexGetSupportSize(dm, point, &supportSize);CHKERRQ(ierr); 1655 if (supportSize != 2) continue; 1656 ierr = DMPlexGetSupport(dm, point, &support);CHKERRQ(ierr); 1657 ierr = DMLabelGetValue(label, support[0], &valA);CHKERRQ(ierr); 1658 ierr = DMLabelGetValue(label, support[1], &valB);CHKERRQ(ierr); 1659 if ((valA == -1) || (valB == -1)) continue; 1660 if (valA*valB > 0) continue; 1661 /* Split the face */ 1662 ierr = DMLabelGetValue(label, point, &valA);CHKERRQ(ierr); 1663 ierr = DMLabelClearValue(label, point, valA);CHKERRQ(ierr); 1664 ierr = DMLabelSetValue(label, point, dim-1);CHKERRQ(ierr); 1665 /* Label its closure: 1666 unmarked: label as unsplit 1667 incident: relabel as split 1668 split: do nothing 1669 */ 1670 { 1671 PetscInt *closure = NULL; 1672 PetscInt closureSize, cl; 1673 1674 ierr = DMPlexGetTransitiveClosure(dm, point, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 1675 for (cl = 0; cl < closureSize*2; cl += 2) { 1676 ierr = DMLabelGetValue(label, closure[cl], &valA);CHKERRQ(ierr); 1677 if (valA == -1) { /* Mark as unsplit */ 1678 ierr = DMLabelGetValue(depthLabel, closure[cl], &dep);CHKERRQ(ierr); 1679 ierr = DMLabelSetValue(label, closure[cl], shift2+dep);CHKERRQ(ierr); 1680 } else if (((valA >= shift) && (valA < shift2)) || ((valA <= -shift) && (valA > -shift2))) { 1681 ierr = DMLabelGetValue(depthLabel, closure[cl], &dep);CHKERRQ(ierr); 1682 ierr = DMLabelClearValue(label, closure[cl], valA);CHKERRQ(ierr); 1683 ierr = DMLabelSetValue(label, closure[cl], dep);CHKERRQ(ierr); 1684 } 1685 } 1686 ierr = DMPlexRestoreTransitiveClosure(dm, point, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 1687 } 1688 } 1689 ierr = ISRestoreIndices(dimIS, &points);CHKERRQ(ierr); 1690 ierr = ISDestroy(&dimIS);CHKERRQ(ierr); 1691 ierr = PetscFree(pMax);CHKERRQ(ierr); 1692 PetscFunctionReturn(0); 1693 } 1694 1695 #undef __FUNCT__ 1696 #define __FUNCT__ "DMPlexCreateHybridMesh" 1697 /*@C 1698 DMPlexCreateHybridMesh - Create a mesh with hybrid cells along an internal interface 1699 1700 Collective on dm 1701 1702 Input Parameters: 1703 + dm - The original DM 1704 - labelName - The label specifying the interface vertices 1705 1706 Output Parameters: 1707 + hybridLabel - The label fully marking the interface 1708 - dmHybrid - The new DM 1709 1710 Level: developer 1711 1712 .seealso: DMPlexConstructCohesiveCells(), DMPlexLabelCohesiveComplete(), DMCreate() 1713 @*/ 1714 PetscErrorCode DMPlexCreateHybridMesh(DM dm, DMLabel label, DMLabel *hybridLabel, DM *dmHybrid) 1715 { 1716 DM idm; 1717 DMLabel subpointMap, hlabel; 1718 PetscInt dim; 1719 PetscErrorCode ierr; 1720 1721 PetscFunctionBegin; 1722 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 1723 if (hybridLabel) PetscValidPointer(hybridLabel, 3); 1724 PetscValidPointer(dmHybrid, 4); 1725 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 1726 ierr = DMPlexCreateSubmesh(dm, label, 1, &idm);CHKERRQ(ierr); 1727 ierr = DMPlexOrient(idm);CHKERRQ(ierr); 1728 ierr = DMPlexGetSubpointMap(idm, &subpointMap);CHKERRQ(ierr); 1729 ierr = DMLabelDuplicate(subpointMap, &hlabel);CHKERRQ(ierr); 1730 ierr = DMLabelClearStratum(hlabel, dim);CHKERRQ(ierr); 1731 ierr = DMPlexLabelCohesiveComplete(dm, hlabel, NULL, PETSC_FALSE, idm);CHKERRQ(ierr); 1732 ierr = DMDestroy(&idm);CHKERRQ(ierr); 1733 ierr = DMPlexConstructCohesiveCells(dm, hlabel, dmHybrid);CHKERRQ(ierr); 1734 if (hybridLabel) *hybridLabel = hlabel; 1735 else {ierr = DMLabelDestroy(&hlabel);CHKERRQ(ierr);} 1736 PetscFunctionReturn(0); 1737 } 1738 1739 #undef __FUNCT__ 1740 #define __FUNCT__ "DMPlexMarkSubmesh_Uninterpolated" 1741 /* Here we need the explicit assumption that: 1742 1743 For any marked cell, the marked vertices constitute a single face 1744 */ 1745 static PetscErrorCode DMPlexMarkSubmesh_Uninterpolated(DM dm, DMLabel vertexLabel, PetscInt value, DMLabel subpointMap, PetscInt *numFaces, PetscInt *nFV, DM subdm) 1746 { 1747 IS subvertexIS = NULL; 1748 const PetscInt *subvertices; 1749 PetscInt *pStart, *pEnd, *pMax, pSize; 1750 PetscInt depth, dim, d, numSubVerticesInitial = 0, v; 1751 PetscErrorCode ierr; 1752 1753 PetscFunctionBegin; 1754 *numFaces = 0; 1755 *nFV = 0; 1756 ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr); 1757 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 1758 pSize = PetscMax(depth, dim) + 1; 1759 ierr = PetscMalloc3(pSize,&pStart,pSize,&pEnd,pSize,&pMax);CHKERRQ(ierr); 1760 ierr = DMPlexGetHybridBounds(dm, depth >= 0 ? &pMax[depth] : NULL, depth>1 ? &pMax[depth-1] : NULL, depth>2 ? &pMax[1] : NULL, &pMax[0]);CHKERRQ(ierr); 1761 for (d = 0; d <= depth; ++d) { 1762 ierr = DMPlexGetDepthStratum(dm, d, &pStart[d], &pEnd[d]);CHKERRQ(ierr); 1763 if (pMax[d] >= 0) pEnd[d] = PetscMin(pEnd[d], pMax[d]); 1764 } 1765 /* Loop over initial vertices and mark all faces in the collective star() */ 1766 if (vertexLabel) {ierr = DMLabelGetStratumIS(vertexLabel, value, &subvertexIS);CHKERRQ(ierr);} 1767 if (subvertexIS) { 1768 ierr = ISGetSize(subvertexIS, &numSubVerticesInitial);CHKERRQ(ierr); 1769 ierr = ISGetIndices(subvertexIS, &subvertices);CHKERRQ(ierr); 1770 } 1771 for (v = 0; v < numSubVerticesInitial; ++v) { 1772 const PetscInt vertex = subvertices[v]; 1773 PetscInt *star = NULL; 1774 PetscInt starSize, s, numCells = 0, c; 1775 1776 ierr = DMPlexGetTransitiveClosure(dm, vertex, PETSC_FALSE, &starSize, &star);CHKERRQ(ierr); 1777 for (s = 0; s < starSize*2; s += 2) { 1778 const PetscInt point = star[s]; 1779 if ((point >= pStart[depth]) && (point < pEnd[depth])) star[numCells++] = point; 1780 } 1781 for (c = 0; c < numCells; ++c) { 1782 const PetscInt cell = star[c]; 1783 PetscInt *closure = NULL; 1784 PetscInt closureSize, cl; 1785 PetscInt cellLoc, numCorners = 0, faceSize = 0; 1786 1787 ierr = DMLabelGetValue(subpointMap, cell, &cellLoc);CHKERRQ(ierr); 1788 if (cellLoc == 2) continue; 1789 if (cellLoc >= 0) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Cell %d has dimension %d in the surface label", cell, cellLoc); 1790 ierr = DMPlexGetTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 1791 for (cl = 0; cl < closureSize*2; cl += 2) { 1792 const PetscInt point = closure[cl]; 1793 PetscInt vertexLoc; 1794 1795 if ((point >= pStart[0]) && (point < pEnd[0])) { 1796 ++numCorners; 1797 ierr = DMLabelGetValue(vertexLabel, point, &vertexLoc);CHKERRQ(ierr); 1798 if (vertexLoc == value) closure[faceSize++] = point; 1799 } 1800 } 1801 if (!(*nFV)) {ierr = DMPlexGetNumFaceVertices(dm, dim, numCorners, nFV);CHKERRQ(ierr);} 1802 if (faceSize > *nFV) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Invalid submesh: Too many vertices %d of an element on the surface", faceSize); 1803 if (faceSize == *nFV) { 1804 const PetscInt *cells = NULL; 1805 PetscInt numCells, nc; 1806 1807 ++(*numFaces); 1808 for (cl = 0; cl < faceSize; ++cl) { 1809 ierr = DMLabelSetValue(subpointMap, closure[cl], 0);CHKERRQ(ierr); 1810 } 1811 ierr = DMPlexGetJoin(dm, faceSize, closure, &numCells, &cells);CHKERRQ(ierr); 1812 for (nc = 0; nc < numCells; ++nc) { 1813 ierr = DMLabelSetValue(subpointMap, cells[nc], 2);CHKERRQ(ierr); 1814 } 1815 ierr = DMPlexRestoreJoin(dm, faceSize, closure, &numCells, &cells);CHKERRQ(ierr); 1816 } 1817 ierr = DMPlexRestoreTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 1818 } 1819 ierr = DMPlexRestoreTransitiveClosure(dm, vertex, PETSC_FALSE, &starSize, &star);CHKERRQ(ierr); 1820 } 1821 if (subvertexIS) { 1822 ierr = ISRestoreIndices(subvertexIS, &subvertices);CHKERRQ(ierr); 1823 } 1824 ierr = ISDestroy(&subvertexIS);CHKERRQ(ierr); 1825 ierr = PetscFree3(pStart,pEnd,pMax);CHKERRQ(ierr); 1826 PetscFunctionReturn(0); 1827 } 1828 1829 #undef __FUNCT__ 1830 #define __FUNCT__ "DMPlexMarkSubmesh_Interpolated" 1831 static PetscErrorCode DMPlexMarkSubmesh_Interpolated(DM dm, DMLabel vertexLabel, PetscInt value, DMLabel subpointMap, DM subdm) 1832 { 1833 IS subvertexIS = NULL; 1834 const PetscInt *subvertices; 1835 PetscInt *pStart, *pEnd, *pMax; 1836 PetscInt dim, d, numSubVerticesInitial = 0, v; 1837 PetscErrorCode ierr; 1838 1839 PetscFunctionBegin; 1840 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 1841 ierr = PetscMalloc3(dim+1,&pStart,dim+1,&pEnd,dim+1,&pMax);CHKERRQ(ierr); 1842 ierr = DMPlexGetHybridBounds(dm, &pMax[dim], dim>1 ? &pMax[dim-1] : NULL, dim > 2 ? &pMax[1] : NULL, &pMax[0]);CHKERRQ(ierr); 1843 for (d = 0; d <= dim; ++d) { 1844 ierr = DMPlexGetDepthStratum(dm, d, &pStart[d], &pEnd[d]);CHKERRQ(ierr); 1845 if (pMax[d] >= 0) pEnd[d] = PetscMin(pEnd[d], pMax[d]); 1846 } 1847 /* Loop over initial vertices and mark all faces in the collective star() */ 1848 if (vertexLabel) { 1849 ierr = DMLabelGetStratumIS(vertexLabel, value, &subvertexIS);CHKERRQ(ierr); 1850 if (subvertexIS) { 1851 ierr = ISGetSize(subvertexIS, &numSubVerticesInitial);CHKERRQ(ierr); 1852 ierr = ISGetIndices(subvertexIS, &subvertices);CHKERRQ(ierr); 1853 } 1854 } 1855 for (v = 0; v < numSubVerticesInitial; ++v) { 1856 const PetscInt vertex = subvertices[v]; 1857 PetscInt *star = NULL; 1858 PetscInt starSize, s, numFaces = 0, f; 1859 1860 ierr = DMPlexGetTransitiveClosure(dm, vertex, PETSC_FALSE, &starSize, &star);CHKERRQ(ierr); 1861 for (s = 0; s < starSize*2; s += 2) { 1862 const PetscInt point = star[s]; 1863 if ((point >= pStart[dim-1]) && (point < pEnd[dim-1])) star[numFaces++] = point; 1864 } 1865 for (f = 0; f < numFaces; ++f) { 1866 const PetscInt face = star[f]; 1867 PetscInt *closure = NULL; 1868 PetscInt closureSize, c; 1869 PetscInt faceLoc; 1870 1871 ierr = DMLabelGetValue(subpointMap, face, &faceLoc);CHKERRQ(ierr); 1872 if (faceLoc == dim-1) continue; 1873 if (faceLoc >= 0) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Face %d has dimension %d in the surface label", face, faceLoc); 1874 ierr = DMPlexGetTransitiveClosure(dm, face, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 1875 for (c = 0; c < closureSize*2; c += 2) { 1876 const PetscInt point = closure[c]; 1877 PetscInt vertexLoc; 1878 1879 if ((point >= pStart[0]) && (point < pEnd[0])) { 1880 ierr = DMLabelGetValue(vertexLabel, point, &vertexLoc);CHKERRQ(ierr); 1881 if (vertexLoc != value) break; 1882 } 1883 } 1884 if (c == closureSize*2) { 1885 const PetscInt *support; 1886 PetscInt supportSize, s; 1887 1888 for (c = 0; c < closureSize*2; c += 2) { 1889 const PetscInt point = closure[c]; 1890 1891 for (d = 0; d < dim; ++d) { 1892 if ((point >= pStart[d]) && (point < pEnd[d])) { 1893 ierr = DMLabelSetValue(subpointMap, point, d);CHKERRQ(ierr); 1894 break; 1895 } 1896 } 1897 } 1898 ierr = DMPlexGetSupportSize(dm, face, &supportSize);CHKERRQ(ierr); 1899 ierr = DMPlexGetSupport(dm, face, &support);CHKERRQ(ierr); 1900 for (s = 0; s < supportSize; ++s) { 1901 ierr = DMLabelSetValue(subpointMap, support[s], dim);CHKERRQ(ierr); 1902 } 1903 } 1904 ierr = DMPlexRestoreTransitiveClosure(dm, face, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 1905 } 1906 ierr = DMPlexRestoreTransitiveClosure(dm, vertex, PETSC_FALSE, &starSize, &star);CHKERRQ(ierr); 1907 } 1908 if (subvertexIS) {ierr = ISRestoreIndices(subvertexIS, &subvertices);CHKERRQ(ierr);} 1909 ierr = ISDestroy(&subvertexIS);CHKERRQ(ierr); 1910 ierr = PetscFree3(pStart,pEnd,pMax);CHKERRQ(ierr); 1911 PetscFunctionReturn(0); 1912 } 1913 1914 #undef __FUNCT__ 1915 #define __FUNCT__ "DMPlexMarkCohesiveSubmesh_Uninterpolated" 1916 static PetscErrorCode DMPlexMarkCohesiveSubmesh_Uninterpolated(DM dm, PetscBool hasLagrange, const char labelname[], PetscInt value, DMLabel subpointMap, PetscInt *numFaces, PetscInt *nFV, PetscInt *subCells[], DM subdm) 1917 { 1918 DMLabel label = NULL; 1919 const PetscInt *cone; 1920 PetscInt dim, cMax, cEnd, c, subc = 0, p, coneSize; 1921 PetscErrorCode ierr; 1922 1923 PetscFunctionBegin; 1924 *numFaces = 0; 1925 *nFV = 0; 1926 if (labelname) {ierr = DMPlexGetLabel(dm, labelname, &label);CHKERRQ(ierr);} 1927 *subCells = NULL; 1928 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 1929 ierr = DMPlexGetHeightStratum(dm, 0, NULL, &cEnd);CHKERRQ(ierr); 1930 ierr = DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);CHKERRQ(ierr); 1931 if (cMax < 0) PetscFunctionReturn(0); 1932 if (label) { 1933 for (c = cMax; c < cEnd; ++c) { 1934 PetscInt val; 1935 1936 ierr = DMLabelGetValue(label, c, &val);CHKERRQ(ierr); 1937 if (val == value) { 1938 ++(*numFaces); 1939 ierr = DMPlexGetConeSize(dm, c, &coneSize);CHKERRQ(ierr); 1940 } 1941 } 1942 } else { 1943 *numFaces = cEnd - cMax; 1944 ierr = DMPlexGetConeSize(dm, cMax, &coneSize);CHKERRQ(ierr); 1945 } 1946 *nFV = hasLagrange ? coneSize/3 : coneSize/2; 1947 ierr = PetscMalloc1(*numFaces *2, subCells);CHKERRQ(ierr); 1948 for (c = cMax; c < cEnd; ++c) { 1949 const PetscInt *cells; 1950 PetscInt numCells; 1951 1952 if (label) { 1953 PetscInt val; 1954 1955 ierr = DMLabelGetValue(label, c, &val);CHKERRQ(ierr); 1956 if (val != value) continue; 1957 } 1958 ierr = DMPlexGetCone(dm, c, &cone);CHKERRQ(ierr); 1959 for (p = 0; p < *nFV; ++p) { 1960 ierr = DMLabelSetValue(subpointMap, cone[p], 0);CHKERRQ(ierr); 1961 } 1962 /* Negative face */ 1963 ierr = DMPlexGetJoin(dm, *nFV, cone, &numCells, &cells);CHKERRQ(ierr); 1964 /* Not true in parallel 1965 if (numCells != 2) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cohesive cells should separate two cells"); */ 1966 for (p = 0; p < numCells; ++p) { 1967 ierr = DMLabelSetValue(subpointMap, cells[p], 2);CHKERRQ(ierr); 1968 (*subCells)[subc++] = cells[p]; 1969 } 1970 ierr = DMPlexRestoreJoin(dm, *nFV, cone, &numCells, &cells);CHKERRQ(ierr); 1971 /* Positive face is not included */ 1972 } 1973 PetscFunctionReturn(0); 1974 } 1975 1976 #undef __FUNCT__ 1977 #define __FUNCT__ "DMPlexMarkCohesiveSubmesh_Interpolated" 1978 static PetscErrorCode DMPlexMarkCohesiveSubmesh_Interpolated(DM dm, DMLabel label, PetscInt value, DMLabel subpointMap, DM subdm) 1979 { 1980 PetscInt *pStart, *pEnd; 1981 PetscInt dim, cMax, cEnd, c, d; 1982 PetscErrorCode ierr; 1983 1984 PetscFunctionBegin; 1985 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 1986 ierr = DMPlexGetHeightStratum(dm, 0, NULL, &cEnd);CHKERRQ(ierr); 1987 ierr = DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);CHKERRQ(ierr); 1988 if (cMax < 0) PetscFunctionReturn(0); 1989 ierr = PetscMalloc2(dim+1,&pStart,dim+1,&pEnd);CHKERRQ(ierr); 1990 for (d = 0; d <= dim; ++d) {ierr = DMPlexGetDepthStratum(dm, d, &pStart[d], &pEnd[d]);CHKERRQ(ierr);} 1991 for (c = cMax; c < cEnd; ++c) { 1992 const PetscInt *cone; 1993 PetscInt *closure = NULL; 1994 PetscInt fconeSize, coneSize, closureSize, cl, val; 1995 1996 if (label) { 1997 ierr = DMLabelGetValue(label, c, &val);CHKERRQ(ierr); 1998 if (val != value) continue; 1999 } 2000 ierr = DMPlexGetConeSize(dm, c, &coneSize);CHKERRQ(ierr); 2001 ierr = DMPlexGetCone(dm, c, &cone);CHKERRQ(ierr); 2002 ierr = DMPlexGetConeSize(dm, cone[0], &fconeSize);CHKERRQ(ierr); 2003 if (coneSize != (fconeSize ? fconeSize : 1) + 2) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cohesive cells should separate two cells"); 2004 /* Negative face */ 2005 ierr = DMPlexGetTransitiveClosure(dm, cone[0], PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 2006 for (cl = 0; cl < closureSize*2; cl += 2) { 2007 const PetscInt point = closure[cl]; 2008 2009 for (d = 0; d <= dim; ++d) { 2010 if ((point >= pStart[d]) && (point < pEnd[d])) { 2011 ierr = DMLabelSetValue(subpointMap, point, d);CHKERRQ(ierr); 2012 break; 2013 } 2014 } 2015 } 2016 ierr = DMPlexRestoreTransitiveClosure(dm, cone[0], PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 2017 /* Cells -- positive face is not included */ 2018 for (cl = 0; cl < 1; ++cl) { 2019 const PetscInt *support; 2020 PetscInt supportSize, s; 2021 2022 ierr = DMPlexGetSupportSize(dm, cone[cl], &supportSize);CHKERRQ(ierr); 2023 /* if (supportSize != 2) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cohesive faces should separate two cells"); */ 2024 ierr = DMPlexGetSupport(dm, cone[cl], &support);CHKERRQ(ierr); 2025 for (s = 0; s < supportSize; ++s) { 2026 ierr = DMLabelSetValue(subpointMap, support[s], dim);CHKERRQ(ierr); 2027 } 2028 } 2029 } 2030 ierr = PetscFree2(pStart, pEnd);CHKERRQ(ierr); 2031 PetscFunctionReturn(0); 2032 } 2033 2034 #undef __FUNCT__ 2035 #define __FUNCT__ "DMPlexGetFaceOrientation" 2036 PetscErrorCode DMPlexGetFaceOrientation(DM dm, PetscInt cell, PetscInt numCorners, PetscInt indices[], PetscInt oppositeVertex, PetscInt origVertices[], PetscInt faceVertices[], PetscBool *posOriented) 2037 { 2038 MPI_Comm comm; 2039 PetscBool posOrient = PETSC_FALSE; 2040 const PetscInt debug = 0; 2041 PetscInt cellDim, faceSize, f; 2042 PetscErrorCode ierr; 2043 2044 PetscFunctionBegin; 2045 ierr = PetscObjectGetComm((PetscObject)dm,&comm);CHKERRQ(ierr); 2046 ierr = DMGetDimension(dm, &cellDim);CHKERRQ(ierr); 2047 if (debug) {PetscPrintf(comm, "cellDim: %d numCorners: %d\n", cellDim, numCorners);CHKERRQ(ierr);} 2048 2049 if (cellDim == 1 && numCorners == 2) { 2050 /* Triangle */ 2051 faceSize = numCorners-1; 2052 posOrient = !(oppositeVertex%2) ? PETSC_TRUE : PETSC_FALSE; 2053 } else if (cellDim == 2 && numCorners == 3) { 2054 /* Triangle */ 2055 faceSize = numCorners-1; 2056 posOrient = !(oppositeVertex%2) ? PETSC_TRUE : PETSC_FALSE; 2057 } else if (cellDim == 3 && numCorners == 4) { 2058 /* Tetrahedron */ 2059 faceSize = numCorners-1; 2060 posOrient = (oppositeVertex%2) ? PETSC_TRUE : PETSC_FALSE; 2061 } else if (cellDim == 1 && numCorners == 3) { 2062 /* Quadratic line */ 2063 faceSize = 1; 2064 posOrient = PETSC_TRUE; 2065 } else if (cellDim == 2 && numCorners == 4) { 2066 /* Quads */ 2067 faceSize = 2; 2068 if ((indices[1] > indices[0]) && (indices[1] - indices[0] == 1)) { 2069 posOrient = PETSC_TRUE; 2070 } else if ((indices[0] == 3) && (indices[1] == 0)) { 2071 posOrient = PETSC_TRUE; 2072 } else { 2073 if (((indices[0] > indices[1]) && (indices[0] - indices[1] == 1)) || ((indices[0] == 0) && (indices[1] == 3))) { 2074 posOrient = PETSC_FALSE; 2075 } else SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid quad crossedge"); 2076 } 2077 } else if (cellDim == 2 && numCorners == 6) { 2078 /* Quadratic triangle (I hate this) */ 2079 /* Edges are determined by the first 2 vertices (corners of edges) */ 2080 const PetscInt faceSizeTri = 3; 2081 PetscInt sortedIndices[3], i, iFace; 2082 PetscBool found = PETSC_FALSE; 2083 PetscInt faceVerticesTriSorted[9] = { 2084 0, 3, 4, /* bottom */ 2085 1, 4, 5, /* right */ 2086 2, 3, 5, /* left */ 2087 }; 2088 PetscInt faceVerticesTri[9] = { 2089 0, 3, 4, /* bottom */ 2090 1, 4, 5, /* right */ 2091 2, 5, 3, /* left */ 2092 }; 2093 2094 faceSize = faceSizeTri; 2095 for (i = 0; i < faceSizeTri; ++i) sortedIndices[i] = indices[i]; 2096 ierr = PetscSortInt(faceSizeTri, sortedIndices);CHKERRQ(ierr); 2097 for (iFace = 0; iFace < 3; ++iFace) { 2098 const PetscInt ii = iFace*faceSizeTri; 2099 PetscInt fVertex, cVertex; 2100 2101 if ((sortedIndices[0] == faceVerticesTriSorted[ii+0]) && 2102 (sortedIndices[1] == faceVerticesTriSorted[ii+1])) { 2103 for (fVertex = 0; fVertex < faceSizeTri; ++fVertex) { 2104 for (cVertex = 0; cVertex < faceSizeTri; ++cVertex) { 2105 if (indices[cVertex] == faceVerticesTri[ii+fVertex]) { 2106 faceVertices[fVertex] = origVertices[cVertex]; 2107 break; 2108 } 2109 } 2110 } 2111 found = PETSC_TRUE; 2112 break; 2113 } 2114 } 2115 if (!found) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid tri crossface"); 2116 if (posOriented) *posOriented = PETSC_TRUE; 2117 PetscFunctionReturn(0); 2118 } else if (cellDim == 2 && numCorners == 9) { 2119 /* Quadratic quad (I hate this) */ 2120 /* Edges are determined by the first 2 vertices (corners of edges) */ 2121 const PetscInt faceSizeQuad = 3; 2122 PetscInt sortedIndices[3], i, iFace; 2123 PetscBool found = PETSC_FALSE; 2124 PetscInt faceVerticesQuadSorted[12] = { 2125 0, 1, 4, /* bottom */ 2126 1, 2, 5, /* right */ 2127 2, 3, 6, /* top */ 2128 0, 3, 7, /* left */ 2129 }; 2130 PetscInt faceVerticesQuad[12] = { 2131 0, 1, 4, /* bottom */ 2132 1, 2, 5, /* right */ 2133 2, 3, 6, /* top */ 2134 3, 0, 7, /* left */ 2135 }; 2136 2137 faceSize = faceSizeQuad; 2138 for (i = 0; i < faceSizeQuad; ++i) sortedIndices[i] = indices[i]; 2139 ierr = PetscSortInt(faceSizeQuad, sortedIndices);CHKERRQ(ierr); 2140 for (iFace = 0; iFace < 4; ++iFace) { 2141 const PetscInt ii = iFace*faceSizeQuad; 2142 PetscInt fVertex, cVertex; 2143 2144 if ((sortedIndices[0] == faceVerticesQuadSorted[ii+0]) && 2145 (sortedIndices[1] == faceVerticesQuadSorted[ii+1])) { 2146 for (fVertex = 0; fVertex < faceSizeQuad; ++fVertex) { 2147 for (cVertex = 0; cVertex < faceSizeQuad; ++cVertex) { 2148 if (indices[cVertex] == faceVerticesQuad[ii+fVertex]) { 2149 faceVertices[fVertex] = origVertices[cVertex]; 2150 break; 2151 } 2152 } 2153 } 2154 found = PETSC_TRUE; 2155 break; 2156 } 2157 } 2158 if (!found) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid quad crossface"); 2159 if (posOriented) *posOriented = PETSC_TRUE; 2160 PetscFunctionReturn(0); 2161 } else if (cellDim == 3 && numCorners == 8) { 2162 /* Hexes 2163 A hex is two oriented quads with the normal of the first 2164 pointing up at the second. 2165 2166 7---6 2167 /| /| 2168 4---5 | 2169 | 1-|-2 2170 |/ |/ 2171 0---3 2172 2173 Faces are determined by the first 4 vertices (corners of faces) */ 2174 const PetscInt faceSizeHex = 4; 2175 PetscInt sortedIndices[4], i, iFace; 2176 PetscBool found = PETSC_FALSE; 2177 PetscInt faceVerticesHexSorted[24] = { 2178 0, 1, 2, 3, /* bottom */ 2179 4, 5, 6, 7, /* top */ 2180 0, 3, 4, 5, /* front */ 2181 2, 3, 5, 6, /* right */ 2182 1, 2, 6, 7, /* back */ 2183 0, 1, 4, 7, /* left */ 2184 }; 2185 PetscInt faceVerticesHex[24] = { 2186 1, 2, 3, 0, /* bottom */ 2187 4, 5, 6, 7, /* top */ 2188 0, 3, 5, 4, /* front */ 2189 3, 2, 6, 5, /* right */ 2190 2, 1, 7, 6, /* back */ 2191 1, 0, 4, 7, /* left */ 2192 }; 2193 2194 faceSize = faceSizeHex; 2195 for (i = 0; i < faceSizeHex; ++i) sortedIndices[i] = indices[i]; 2196 ierr = PetscSortInt(faceSizeHex, sortedIndices);CHKERRQ(ierr); 2197 for (iFace = 0; iFace < 6; ++iFace) { 2198 const PetscInt ii = iFace*faceSizeHex; 2199 PetscInt fVertex, cVertex; 2200 2201 if ((sortedIndices[0] == faceVerticesHexSorted[ii+0]) && 2202 (sortedIndices[1] == faceVerticesHexSorted[ii+1]) && 2203 (sortedIndices[2] == faceVerticesHexSorted[ii+2]) && 2204 (sortedIndices[3] == faceVerticesHexSorted[ii+3])) { 2205 for (fVertex = 0; fVertex < faceSizeHex; ++fVertex) { 2206 for (cVertex = 0; cVertex < faceSizeHex; ++cVertex) { 2207 if (indices[cVertex] == faceVerticesHex[ii+fVertex]) { 2208 faceVertices[fVertex] = origVertices[cVertex]; 2209 break; 2210 } 2211 } 2212 } 2213 found = PETSC_TRUE; 2214 break; 2215 } 2216 } 2217 if (!found) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid hex crossface"); 2218 if (posOriented) *posOriented = PETSC_TRUE; 2219 PetscFunctionReturn(0); 2220 } else if (cellDim == 3 && numCorners == 10) { 2221 /* Quadratic tet */ 2222 /* Faces are determined by the first 3 vertices (corners of faces) */ 2223 const PetscInt faceSizeTet = 6; 2224 PetscInt sortedIndices[6], i, iFace; 2225 PetscBool found = PETSC_FALSE; 2226 PetscInt faceVerticesTetSorted[24] = { 2227 0, 1, 2, 6, 7, 8, /* bottom */ 2228 0, 3, 4, 6, 7, 9, /* front */ 2229 1, 4, 5, 7, 8, 9, /* right */ 2230 2, 3, 5, 6, 8, 9, /* left */ 2231 }; 2232 PetscInt faceVerticesTet[24] = { 2233 0, 1, 2, 6, 7, 8, /* bottom */ 2234 0, 4, 3, 6, 7, 9, /* front */ 2235 1, 5, 4, 7, 8, 9, /* right */ 2236 2, 3, 5, 8, 6, 9, /* left */ 2237 }; 2238 2239 faceSize = faceSizeTet; 2240 for (i = 0; i < faceSizeTet; ++i) sortedIndices[i] = indices[i]; 2241 ierr = PetscSortInt(faceSizeTet, sortedIndices);CHKERRQ(ierr); 2242 for (iFace=0; iFace < 4; ++iFace) { 2243 const PetscInt ii = iFace*faceSizeTet; 2244 PetscInt fVertex, cVertex; 2245 2246 if ((sortedIndices[0] == faceVerticesTetSorted[ii+0]) && 2247 (sortedIndices[1] == faceVerticesTetSorted[ii+1]) && 2248 (sortedIndices[2] == faceVerticesTetSorted[ii+2]) && 2249 (sortedIndices[3] == faceVerticesTetSorted[ii+3])) { 2250 for (fVertex = 0; fVertex < faceSizeTet; ++fVertex) { 2251 for (cVertex = 0; cVertex < faceSizeTet; ++cVertex) { 2252 if (indices[cVertex] == faceVerticesTet[ii+fVertex]) { 2253 faceVertices[fVertex] = origVertices[cVertex]; 2254 break; 2255 } 2256 } 2257 } 2258 found = PETSC_TRUE; 2259 break; 2260 } 2261 } 2262 if (!found) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid tet crossface"); 2263 if (posOriented) *posOriented = PETSC_TRUE; 2264 PetscFunctionReturn(0); 2265 } else if (cellDim == 3 && numCorners == 27) { 2266 /* Quadratic hexes (I hate this) 2267 A hex is two oriented quads with the normal of the first 2268 pointing up at the second. 2269 2270 7---6 2271 /| /| 2272 4---5 | 2273 | 3-|-2 2274 |/ |/ 2275 0---1 2276 2277 Faces are determined by the first 4 vertices (corners of faces) */ 2278 const PetscInt faceSizeQuadHex = 9; 2279 PetscInt sortedIndices[9], i, iFace; 2280 PetscBool found = PETSC_FALSE; 2281 PetscInt faceVerticesQuadHexSorted[54] = { 2282 0, 1, 2, 3, 8, 9, 10, 11, 24, /* bottom */ 2283 4, 5, 6, 7, 12, 13, 14, 15, 25, /* top */ 2284 0, 1, 4, 5, 8, 12, 16, 17, 22, /* front */ 2285 1, 2, 5, 6, 9, 13, 17, 18, 21, /* right */ 2286 2, 3, 6, 7, 10, 14, 18, 19, 23, /* back */ 2287 0, 3, 4, 7, 11, 15, 16, 19, 20, /* left */ 2288 }; 2289 PetscInt faceVerticesQuadHex[54] = { 2290 3, 2, 1, 0, 10, 9, 8, 11, 24, /* bottom */ 2291 4, 5, 6, 7, 12, 13, 14, 15, 25, /* top */ 2292 0, 1, 5, 4, 8, 17, 12, 16, 22, /* front */ 2293 1, 2, 6, 5, 9, 18, 13, 17, 21, /* right */ 2294 2, 3, 7, 6, 10, 19, 14, 18, 23, /* back */ 2295 3, 0, 4, 7, 11, 16, 15, 19, 20 /* left */ 2296 }; 2297 2298 faceSize = faceSizeQuadHex; 2299 for (i = 0; i < faceSizeQuadHex; ++i) sortedIndices[i] = indices[i]; 2300 ierr = PetscSortInt(faceSizeQuadHex, sortedIndices);CHKERRQ(ierr); 2301 for (iFace = 0; iFace < 6; ++iFace) { 2302 const PetscInt ii = iFace*faceSizeQuadHex; 2303 PetscInt fVertex, cVertex; 2304 2305 if ((sortedIndices[0] == faceVerticesQuadHexSorted[ii+0]) && 2306 (sortedIndices[1] == faceVerticesQuadHexSorted[ii+1]) && 2307 (sortedIndices[2] == faceVerticesQuadHexSorted[ii+2]) && 2308 (sortedIndices[3] == faceVerticesQuadHexSorted[ii+3])) { 2309 for (fVertex = 0; fVertex < faceSizeQuadHex; ++fVertex) { 2310 for (cVertex = 0; cVertex < faceSizeQuadHex; ++cVertex) { 2311 if (indices[cVertex] == faceVerticesQuadHex[ii+fVertex]) { 2312 faceVertices[fVertex] = origVertices[cVertex]; 2313 break; 2314 } 2315 } 2316 } 2317 found = PETSC_TRUE; 2318 break; 2319 } 2320 } 2321 if (!found) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid hex crossface"); 2322 if (posOriented) *posOriented = PETSC_TRUE; 2323 PetscFunctionReturn(0); 2324 } else SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Unknown cell type for faceOrientation()."); 2325 if (!posOrient) { 2326 if (debug) {ierr = PetscPrintf(comm, " Reversing initial face orientation\n");CHKERRQ(ierr);} 2327 for (f = 0; f < faceSize; ++f) faceVertices[f] = origVertices[faceSize-1 - f]; 2328 } else { 2329 if (debug) {ierr = PetscPrintf(comm, " Keeping initial face orientation\n");CHKERRQ(ierr);} 2330 for (f = 0; f < faceSize; ++f) faceVertices[f] = origVertices[f]; 2331 } 2332 if (posOriented) *posOriented = posOrient; 2333 PetscFunctionReturn(0); 2334 } 2335 2336 #undef __FUNCT__ 2337 #define __FUNCT__ "DMPlexGetOrientedFace" 2338 /* 2339 Given a cell and a face, as a set of vertices, 2340 return the oriented face, as a set of vertices, in faceVertices 2341 The orientation is such that the face normal points out of the cell 2342 */ 2343 PetscErrorCode DMPlexGetOrientedFace(DM dm, PetscInt cell, PetscInt faceSize, const PetscInt face[], PetscInt numCorners, PetscInt indices[], PetscInt origVertices[], PetscInt faceVertices[], PetscBool *posOriented) 2344 { 2345 const PetscInt *cone = NULL; 2346 PetscInt coneSize, v, f, v2; 2347 PetscInt oppositeVertex = -1; 2348 PetscErrorCode ierr; 2349 2350 PetscFunctionBegin; 2351 ierr = DMPlexGetConeSize(dm, cell, &coneSize);CHKERRQ(ierr); 2352 ierr = DMPlexGetCone(dm, cell, &cone);CHKERRQ(ierr); 2353 for (v = 0, v2 = 0; v < coneSize; ++v) { 2354 PetscBool found = PETSC_FALSE; 2355 2356 for (f = 0; f < faceSize; ++f) { 2357 if (face[f] == cone[v]) { 2358 found = PETSC_TRUE; break; 2359 } 2360 } 2361 if (found) { 2362 indices[v2] = v; 2363 origVertices[v2] = cone[v]; 2364 ++v2; 2365 } else { 2366 oppositeVertex = v; 2367 } 2368 } 2369 ierr = DMPlexGetFaceOrientation(dm, cell, numCorners, indices, oppositeVertex, origVertices, faceVertices, posOriented);CHKERRQ(ierr); 2370 PetscFunctionReturn(0); 2371 } 2372 2373 #undef __FUNCT__ 2374 #define __FUNCT__ "DMPlexInsertFace_Internal" 2375 /* 2376 DMPlexInsertFace_Internal - Puts a face into the mesh 2377 2378 Not collective 2379 2380 Input Parameters: 2381 + dm - The DMPlex 2382 . numFaceVertex - The number of vertices in the face 2383 . faceVertices - The vertices in the face for dm 2384 . subfaceVertices - The vertices in the face for subdm 2385 . numCorners - The number of vertices in the cell 2386 . cell - A cell in dm containing the face 2387 . subcell - A cell in subdm containing the face 2388 . firstFace - First face in the mesh 2389 - newFacePoint - Next face in the mesh 2390 2391 Output Parameters: 2392 . newFacePoint - Contains next face point number on input, updated on output 2393 2394 Level: developer 2395 */ 2396 static PetscErrorCode DMPlexInsertFace_Internal(DM dm, DM subdm, PetscInt numFaceVertices, const PetscInt faceVertices[], const PetscInt subfaceVertices[], PetscInt numCorners, PetscInt cell, PetscInt subcell, PetscInt firstFace, PetscInt *newFacePoint) 2397 { 2398 MPI_Comm comm; 2399 DM_Plex *submesh = (DM_Plex*) subdm->data; 2400 const PetscInt *faces; 2401 PetscInt numFaces, coneSize; 2402 PetscErrorCode ierr; 2403 2404 PetscFunctionBegin; 2405 ierr = PetscObjectGetComm((PetscObject)dm,&comm);CHKERRQ(ierr); 2406 ierr = DMPlexGetConeSize(subdm, subcell, &coneSize);CHKERRQ(ierr); 2407 if (coneSize != 1) SETERRQ2(comm, PETSC_ERR_ARG_OUTOFRANGE, "Cone size of cell %d is %d != 1", cell, coneSize); 2408 #if 0 2409 /* Cannot use this because support() has not been constructed yet */ 2410 ierr = DMPlexGetJoin(subdm, numFaceVertices, subfaceVertices, &numFaces, &faces);CHKERRQ(ierr); 2411 #else 2412 { 2413 PetscInt f; 2414 2415 numFaces = 0; 2416 ierr = DMGetWorkArray(subdm, 1, PETSC_INT, (void **) &faces);CHKERRQ(ierr); 2417 for (f = firstFace; f < *newFacePoint; ++f) { 2418 PetscInt dof, off, d; 2419 2420 ierr = PetscSectionGetDof(submesh->coneSection, f, &dof);CHKERRQ(ierr); 2421 ierr = PetscSectionGetOffset(submesh->coneSection, f, &off);CHKERRQ(ierr); 2422 /* Yes, I know this is quadratic, but I expect the sizes to be <5 */ 2423 for (d = 0; d < dof; ++d) { 2424 const PetscInt p = submesh->cones[off+d]; 2425 PetscInt v; 2426 2427 for (v = 0; v < numFaceVertices; ++v) { 2428 if (subfaceVertices[v] == p) break; 2429 } 2430 if (v == numFaceVertices) break; 2431 } 2432 if (d == dof) { 2433 numFaces = 1; 2434 ((PetscInt*) faces)[0] = f; 2435 } 2436 } 2437 } 2438 #endif 2439 if (numFaces > 1) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Vertex set had %d faces, not one", numFaces); 2440 else if (numFaces == 1) { 2441 /* Add the other cell neighbor for this face */ 2442 ierr = DMPlexSetCone(subdm, subcell, faces);CHKERRQ(ierr); 2443 } else { 2444 PetscInt *indices, *origVertices, *orientedVertices, *orientedSubVertices, v, ov; 2445 PetscBool posOriented; 2446 2447 ierr = DMGetWorkArray(subdm, 4*numFaceVertices * sizeof(PetscInt), PETSC_INT, &orientedVertices);CHKERRQ(ierr); 2448 origVertices = &orientedVertices[numFaceVertices]; 2449 indices = &orientedVertices[numFaceVertices*2]; 2450 orientedSubVertices = &orientedVertices[numFaceVertices*3]; 2451 ierr = DMPlexGetOrientedFace(dm, cell, numFaceVertices, faceVertices, numCorners, indices, origVertices, orientedVertices, &posOriented);CHKERRQ(ierr); 2452 /* TODO: I know that routine should return a permutation, not the indices */ 2453 for (v = 0; v < numFaceVertices; ++v) { 2454 const PetscInt vertex = faceVertices[v], subvertex = subfaceVertices[v]; 2455 for (ov = 0; ov < numFaceVertices; ++ov) { 2456 if (orientedVertices[ov] == vertex) { 2457 orientedSubVertices[ov] = subvertex; 2458 break; 2459 } 2460 } 2461 if (ov == numFaceVertices) SETERRQ1(comm, PETSC_ERR_PLIB, "Could not find face vertex %d in orientated set", vertex); 2462 } 2463 ierr = DMPlexSetCone(subdm, *newFacePoint, orientedSubVertices);CHKERRQ(ierr); 2464 ierr = DMPlexSetCone(subdm, subcell, newFacePoint);CHKERRQ(ierr); 2465 ierr = DMRestoreWorkArray(subdm, 4*numFaceVertices * sizeof(PetscInt), PETSC_INT, &orientedVertices);CHKERRQ(ierr); 2466 ++(*newFacePoint); 2467 } 2468 #if 0 2469 ierr = DMPlexRestoreJoin(subdm, numFaceVertices, subfaceVertices, &numFaces, &faces);CHKERRQ(ierr); 2470 #else 2471 ierr = DMRestoreWorkArray(subdm, 1, PETSC_INT, (void **) &faces);CHKERRQ(ierr); 2472 #endif 2473 PetscFunctionReturn(0); 2474 } 2475 2476 #undef __FUNCT__ 2477 #define __FUNCT__ "DMPlexCreateSubmesh_Uninterpolated" 2478 static PetscErrorCode DMPlexCreateSubmesh_Uninterpolated(DM dm, DMLabel vertexLabel, PetscInt value, DM subdm) 2479 { 2480 MPI_Comm comm; 2481 DMLabel subpointMap; 2482 IS subvertexIS, subcellIS; 2483 const PetscInt *subVertices, *subCells; 2484 PetscInt numSubVertices, firstSubVertex, numSubCells; 2485 PetscInt *subface, maxConeSize, numSubFaces = 0, firstSubFace, newFacePoint, nFV = 0; 2486 PetscInt vStart, vEnd, c, f; 2487 PetscErrorCode ierr; 2488 2489 PetscFunctionBegin; 2490 ierr = PetscObjectGetComm((PetscObject)dm,&comm);CHKERRQ(ierr); 2491 /* Create subpointMap which marks the submesh */ 2492 ierr = DMLabelCreate("subpoint_map", &subpointMap);CHKERRQ(ierr); 2493 ierr = DMPlexSetSubpointMap(subdm, subpointMap);CHKERRQ(ierr); 2494 ierr = DMLabelDestroy(&subpointMap);CHKERRQ(ierr); 2495 if (vertexLabel) {ierr = DMPlexMarkSubmesh_Uninterpolated(dm, vertexLabel, value, subpointMap, &numSubFaces, &nFV, subdm);CHKERRQ(ierr);} 2496 /* Setup chart */ 2497 ierr = DMLabelGetStratumSize(subpointMap, 0, &numSubVertices);CHKERRQ(ierr); 2498 ierr = DMLabelGetStratumSize(subpointMap, 2, &numSubCells);CHKERRQ(ierr); 2499 ierr = DMPlexSetChart(subdm, 0, numSubCells+numSubFaces+numSubVertices);CHKERRQ(ierr); 2500 ierr = DMPlexSetVTKCellHeight(subdm, 1);CHKERRQ(ierr); 2501 /* Set cone sizes */ 2502 firstSubVertex = numSubCells; 2503 firstSubFace = numSubCells+numSubVertices; 2504 newFacePoint = firstSubFace; 2505 ierr = DMLabelGetStratumIS(subpointMap, 0, &subvertexIS);CHKERRQ(ierr); 2506 if (subvertexIS) {ierr = ISGetIndices(subvertexIS, &subVertices);CHKERRQ(ierr);} 2507 ierr = DMLabelGetStratumIS(subpointMap, 2, &subcellIS);CHKERRQ(ierr); 2508 if (subcellIS) {ierr = ISGetIndices(subcellIS, &subCells);CHKERRQ(ierr);} 2509 for (c = 0; c < numSubCells; ++c) { 2510 ierr = DMPlexSetConeSize(subdm, c, 1);CHKERRQ(ierr); 2511 } 2512 for (f = firstSubFace; f < firstSubFace+numSubFaces; ++f) { 2513 ierr = DMPlexSetConeSize(subdm, f, nFV);CHKERRQ(ierr); 2514 } 2515 ierr = DMSetUp(subdm);CHKERRQ(ierr); 2516 /* Create face cones */ 2517 ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr); 2518 ierr = DMPlexGetMaxSizes(dm, &maxConeSize, NULL);CHKERRQ(ierr); 2519 ierr = DMGetWorkArray(subdm, maxConeSize, PETSC_INT, (void**) &subface);CHKERRQ(ierr); 2520 for (c = 0; c < numSubCells; ++c) { 2521 const PetscInt cell = subCells[c]; 2522 const PetscInt subcell = c; 2523 PetscInt *closure = NULL; 2524 PetscInt closureSize, cl, numCorners = 0, faceSize = 0; 2525 2526 ierr = DMPlexGetTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 2527 for (cl = 0; cl < closureSize*2; cl += 2) { 2528 const PetscInt point = closure[cl]; 2529 PetscInt subVertex; 2530 2531 if ((point >= vStart) && (point < vEnd)) { 2532 ++numCorners; 2533 ierr = PetscFindInt(point, numSubVertices, subVertices, &subVertex);CHKERRQ(ierr); 2534 if (subVertex >= 0) { 2535 closure[faceSize] = point; 2536 subface[faceSize] = firstSubVertex+subVertex; 2537 ++faceSize; 2538 } 2539 } 2540 } 2541 if (faceSize > nFV) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Invalid submesh: Too many vertices %d of an element on the surface", faceSize); 2542 if (faceSize == nFV) { 2543 ierr = DMPlexInsertFace_Internal(dm, subdm, faceSize, closure, subface, numCorners, cell, subcell, firstSubFace, &newFacePoint);CHKERRQ(ierr); 2544 } 2545 ierr = DMPlexRestoreTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 2546 } 2547 ierr = DMRestoreWorkArray(subdm, maxConeSize, PETSC_INT, (void**) &subface);CHKERRQ(ierr); 2548 ierr = DMPlexSymmetrize(subdm);CHKERRQ(ierr); 2549 ierr = DMPlexStratify(subdm);CHKERRQ(ierr); 2550 /* Build coordinates */ 2551 { 2552 PetscSection coordSection, subCoordSection; 2553 Vec coordinates, subCoordinates; 2554 PetscScalar *coords, *subCoords; 2555 PetscInt numComp, coordSize, v; 2556 2557 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 2558 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr); 2559 ierr = DMGetCoordinateSection(subdm, &subCoordSection);CHKERRQ(ierr); 2560 ierr = PetscSectionSetNumFields(subCoordSection, 1);CHKERRQ(ierr); 2561 ierr = PetscSectionGetFieldComponents(coordSection, 0, &numComp);CHKERRQ(ierr); 2562 ierr = PetscSectionSetFieldComponents(subCoordSection, 0, numComp);CHKERRQ(ierr); 2563 ierr = PetscSectionSetChart(subCoordSection, firstSubVertex, firstSubVertex+numSubVertices);CHKERRQ(ierr); 2564 for (v = 0; v < numSubVertices; ++v) { 2565 const PetscInt vertex = subVertices[v]; 2566 const PetscInt subvertex = firstSubVertex+v; 2567 PetscInt dof; 2568 2569 ierr = PetscSectionGetDof(coordSection, vertex, &dof);CHKERRQ(ierr); 2570 ierr = PetscSectionSetDof(subCoordSection, subvertex, dof);CHKERRQ(ierr); 2571 ierr = PetscSectionSetFieldDof(subCoordSection, subvertex, 0, dof);CHKERRQ(ierr); 2572 } 2573 ierr = PetscSectionSetUp(subCoordSection);CHKERRQ(ierr); 2574 ierr = PetscSectionGetStorageSize(subCoordSection, &coordSize);CHKERRQ(ierr); 2575 ierr = VecCreate(comm, &subCoordinates);CHKERRQ(ierr); 2576 ierr = VecSetSizes(subCoordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 2577 ierr = VecSetType(subCoordinates,VECSTANDARD);CHKERRQ(ierr); 2578 if (coordSize) { 2579 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 2580 ierr = VecGetArray(subCoordinates, &subCoords);CHKERRQ(ierr); 2581 for (v = 0; v < numSubVertices; ++v) { 2582 const PetscInt vertex = subVertices[v]; 2583 const PetscInt subvertex = firstSubVertex+v; 2584 PetscInt dof, off, sdof, soff, d; 2585 2586 ierr = PetscSectionGetDof(coordSection, vertex, &dof);CHKERRQ(ierr); 2587 ierr = PetscSectionGetOffset(coordSection, vertex, &off);CHKERRQ(ierr); 2588 ierr = PetscSectionGetDof(subCoordSection, subvertex, &sdof);CHKERRQ(ierr); 2589 ierr = PetscSectionGetOffset(subCoordSection, subvertex, &soff);CHKERRQ(ierr); 2590 if (dof != sdof) SETERRQ4(comm, PETSC_ERR_PLIB, "Coordinate dimension %d on subvertex %d, vertex %d should be %d", sdof, subvertex, vertex, dof); 2591 for (d = 0; d < dof; ++d) subCoords[soff+d] = coords[off+d]; 2592 } 2593 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 2594 ierr = VecRestoreArray(subCoordinates, &subCoords);CHKERRQ(ierr); 2595 } 2596 ierr = DMSetCoordinatesLocal(subdm, subCoordinates);CHKERRQ(ierr); 2597 ierr = VecDestroy(&subCoordinates);CHKERRQ(ierr); 2598 } 2599 /* Cleanup */ 2600 if (subvertexIS) {ierr = ISRestoreIndices(subvertexIS, &subVertices);CHKERRQ(ierr);} 2601 ierr = ISDestroy(&subvertexIS);CHKERRQ(ierr); 2602 if (subcellIS) {ierr = ISRestoreIndices(subcellIS, &subCells);CHKERRQ(ierr);} 2603 ierr = ISDestroy(&subcellIS);CHKERRQ(ierr); 2604 PetscFunctionReturn(0); 2605 } 2606 2607 #undef __FUNCT__ 2608 #define __FUNCT__ "DMPlexFilterPoint_Internal" 2609 PETSC_STATIC_INLINE PetscInt DMPlexFilterPoint_Internal(PetscInt point, PetscInt firstSubPoint, PetscInt numSubPoints, const PetscInt subPoints[]) 2610 { 2611 PetscInt subPoint; 2612 PetscErrorCode ierr; 2613 2614 ierr = PetscFindInt(point, numSubPoints, subPoints, &subPoint); if (ierr < 0) return ierr; 2615 return subPoint < 0 ? subPoint : firstSubPoint+subPoint; 2616 } 2617 2618 #undef __FUNCT__ 2619 #define __FUNCT__ "DMPlexCreateSubmeshGeneric_Interpolated" 2620 static PetscErrorCode DMPlexCreateSubmeshGeneric_Interpolated(DM dm, DMLabel label, PetscInt value, PetscBool isCohesive, PetscInt cellHeight, DM subdm) 2621 { 2622 MPI_Comm comm; 2623 DMLabel subpointMap; 2624 IS *subpointIS; 2625 const PetscInt **subpoints; 2626 PetscInt *numSubPoints, *firstSubPoint, *coneNew, *orntNew; 2627 PetscInt totSubPoints = 0, maxConeSize, dim, p, d, v; 2628 PetscErrorCode ierr; 2629 2630 PetscFunctionBegin; 2631 ierr = PetscObjectGetComm((PetscObject)dm,&comm);CHKERRQ(ierr); 2632 /* Create subpointMap which marks the submesh */ 2633 ierr = DMLabelCreate("subpoint_map", &subpointMap);CHKERRQ(ierr); 2634 ierr = DMPlexSetSubpointMap(subdm, subpointMap);CHKERRQ(ierr); 2635 if (cellHeight) { 2636 if (isCohesive) {ierr = DMPlexMarkCohesiveSubmesh_Interpolated(dm, label, value, subpointMap, subdm);CHKERRQ(ierr);} 2637 else {ierr = DMPlexMarkSubmesh_Interpolated(dm, label, value, subpointMap, subdm);CHKERRQ(ierr);} 2638 } else { 2639 DMLabel depth; 2640 IS pointIS; 2641 const PetscInt *points; 2642 PetscInt numPoints; 2643 2644 ierr = DMPlexGetDepthLabel(dm, &depth);CHKERRQ(ierr); 2645 ierr = DMLabelGetStratumSize(label, value, &numPoints);CHKERRQ(ierr); 2646 ierr = DMLabelGetStratumIS(label, value, &pointIS);CHKERRQ(ierr); 2647 ierr = ISGetIndices(pointIS, &points);CHKERRQ(ierr); 2648 for (p = 0; p < numPoints; ++p) { 2649 PetscInt *closure = NULL; 2650 PetscInt closureSize, c, pdim; 2651 2652 ierr = DMPlexGetTransitiveClosure(dm, points[p], PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 2653 for (c = 0; c < closureSize*2; c += 2) { 2654 ierr = DMLabelGetValue(depth, closure[c], &pdim);CHKERRQ(ierr); 2655 ierr = DMLabelSetValue(subpointMap, closure[c], pdim);CHKERRQ(ierr); 2656 } 2657 ierr = DMPlexRestoreTransitiveClosure(dm, points[p], PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 2658 } 2659 ierr = ISRestoreIndices(pointIS, &points);CHKERRQ(ierr); 2660 ierr = ISDestroy(&pointIS);CHKERRQ(ierr); 2661 } 2662 ierr = DMLabelDestroy(&subpointMap);CHKERRQ(ierr); 2663 /* Setup chart */ 2664 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 2665 ierr = PetscMalloc4(dim+1,&numSubPoints,dim+1,&firstSubPoint,dim+1,&subpointIS,dim+1,&subpoints);CHKERRQ(ierr); 2666 for (d = 0; d <= dim; ++d) { 2667 ierr = DMLabelGetStratumSize(subpointMap, d, &numSubPoints[d]);CHKERRQ(ierr); 2668 totSubPoints += numSubPoints[d]; 2669 } 2670 ierr = DMPlexSetChart(subdm, 0, totSubPoints);CHKERRQ(ierr); 2671 ierr = DMPlexSetVTKCellHeight(subdm, cellHeight);CHKERRQ(ierr); 2672 /* Set cone sizes */ 2673 firstSubPoint[dim] = 0; 2674 firstSubPoint[0] = firstSubPoint[dim] + numSubPoints[dim]; 2675 if (dim > 1) {firstSubPoint[dim-1] = firstSubPoint[0] + numSubPoints[0];} 2676 if (dim > 2) {firstSubPoint[dim-2] = firstSubPoint[dim-1] + numSubPoints[dim-1];} 2677 for (d = 0; d <= dim; ++d) { 2678 ierr = DMLabelGetStratumIS(subpointMap, d, &subpointIS[d]);CHKERRQ(ierr); 2679 if (subpointIS[d]) {ierr = ISGetIndices(subpointIS[d], &subpoints[d]);CHKERRQ(ierr);} 2680 } 2681 for (d = 0; d <= dim; ++d) { 2682 for (p = 0; p < numSubPoints[d]; ++p) { 2683 const PetscInt point = subpoints[d][p]; 2684 const PetscInt subpoint = firstSubPoint[d] + p; 2685 const PetscInt *cone; 2686 PetscInt coneSize, coneSizeNew, c, val; 2687 2688 ierr = DMPlexGetConeSize(dm, point, &coneSize);CHKERRQ(ierr); 2689 ierr = DMPlexSetConeSize(subdm, subpoint, coneSize);CHKERRQ(ierr); 2690 if (cellHeight && (d == dim)) { 2691 ierr = DMPlexGetCone(dm, point, &cone);CHKERRQ(ierr); 2692 for (c = 0, coneSizeNew = 0; c < coneSize; ++c) { 2693 ierr = DMLabelGetValue(subpointMap, cone[c], &val);CHKERRQ(ierr); 2694 if (val >= 0) coneSizeNew++; 2695 } 2696 ierr = DMPlexSetConeSize(subdm, subpoint, coneSizeNew);CHKERRQ(ierr); 2697 } 2698 } 2699 } 2700 ierr = DMSetUp(subdm);CHKERRQ(ierr); 2701 /* Set cones */ 2702 ierr = DMPlexGetMaxSizes(dm, &maxConeSize, NULL);CHKERRQ(ierr); 2703 ierr = PetscMalloc2(maxConeSize,&coneNew,maxConeSize,&orntNew);CHKERRQ(ierr); 2704 for (d = 0; d <= dim; ++d) { 2705 for (p = 0; p < numSubPoints[d]; ++p) { 2706 const PetscInt point = subpoints[d][p]; 2707 const PetscInt subpoint = firstSubPoint[d] + p; 2708 const PetscInt *cone, *ornt; 2709 PetscInt coneSize, subconeSize, coneSizeNew, c, subc; 2710 2711 ierr = DMPlexGetConeSize(dm, point, &coneSize);CHKERRQ(ierr); 2712 ierr = DMPlexGetConeSize(subdm, subpoint, &subconeSize);CHKERRQ(ierr); 2713 ierr = DMPlexGetCone(dm, point, &cone);CHKERRQ(ierr); 2714 ierr = DMPlexGetConeOrientation(dm, point, &ornt);CHKERRQ(ierr); 2715 for (c = 0, coneSizeNew = 0; c < coneSize; ++c) { 2716 ierr = PetscFindInt(cone[c], numSubPoints[d-1], subpoints[d-1], &subc);CHKERRQ(ierr); 2717 if (subc >= 0) { 2718 coneNew[coneSizeNew] = firstSubPoint[d-1] + subc; 2719 orntNew[coneSizeNew] = ornt[c]; 2720 ++coneSizeNew; 2721 } 2722 } 2723 if (coneSizeNew != subconeSize) SETERRQ2(comm, PETSC_ERR_PLIB, "Number of cone points located %d does not match subcone size %d", coneSizeNew, subconeSize); 2724 ierr = DMPlexSetCone(subdm, subpoint, coneNew);CHKERRQ(ierr); 2725 ierr = DMPlexSetConeOrientation(subdm, subpoint, orntNew);CHKERRQ(ierr); 2726 } 2727 } 2728 ierr = PetscFree2(coneNew,orntNew);CHKERRQ(ierr); 2729 ierr = DMPlexSymmetrize(subdm);CHKERRQ(ierr); 2730 ierr = DMPlexStratify(subdm);CHKERRQ(ierr); 2731 /* Build coordinates */ 2732 { 2733 PetscSection coordSection, subCoordSection; 2734 Vec coordinates, subCoordinates; 2735 PetscScalar *coords, *subCoords; 2736 PetscInt numComp, coordSize; 2737 2738 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 2739 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr); 2740 ierr = DMGetCoordinateSection(subdm, &subCoordSection);CHKERRQ(ierr); 2741 ierr = PetscSectionSetNumFields(subCoordSection, 1);CHKERRQ(ierr); 2742 ierr = PetscSectionGetFieldComponents(coordSection, 0, &numComp);CHKERRQ(ierr); 2743 ierr = PetscSectionSetFieldComponents(subCoordSection, 0, numComp);CHKERRQ(ierr); 2744 ierr = PetscSectionSetChart(subCoordSection, firstSubPoint[0], firstSubPoint[0]+numSubPoints[0]);CHKERRQ(ierr); 2745 for (v = 0; v < numSubPoints[0]; ++v) { 2746 const PetscInt vertex = subpoints[0][v]; 2747 const PetscInt subvertex = firstSubPoint[0]+v; 2748 PetscInt dof; 2749 2750 ierr = PetscSectionGetDof(coordSection, vertex, &dof);CHKERRQ(ierr); 2751 ierr = PetscSectionSetDof(subCoordSection, subvertex, dof);CHKERRQ(ierr); 2752 ierr = PetscSectionSetFieldDof(subCoordSection, subvertex, 0, dof);CHKERRQ(ierr); 2753 } 2754 ierr = PetscSectionSetUp(subCoordSection);CHKERRQ(ierr); 2755 ierr = PetscSectionGetStorageSize(subCoordSection, &coordSize);CHKERRQ(ierr); 2756 ierr = VecCreate(comm, &subCoordinates);CHKERRQ(ierr); 2757 ierr = VecSetSizes(subCoordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 2758 ierr = VecSetType(subCoordinates,VECSTANDARD);CHKERRQ(ierr); 2759 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 2760 ierr = VecGetArray(subCoordinates, &subCoords);CHKERRQ(ierr); 2761 for (v = 0; v < numSubPoints[0]; ++v) { 2762 const PetscInt vertex = subpoints[0][v]; 2763 const PetscInt subvertex = firstSubPoint[0]+v; 2764 PetscInt dof, off, sdof, soff, d; 2765 2766 ierr = PetscSectionGetDof(coordSection, vertex, &dof);CHKERRQ(ierr); 2767 ierr = PetscSectionGetOffset(coordSection, vertex, &off);CHKERRQ(ierr); 2768 ierr = PetscSectionGetDof(subCoordSection, subvertex, &sdof);CHKERRQ(ierr); 2769 ierr = PetscSectionGetOffset(subCoordSection, subvertex, &soff);CHKERRQ(ierr); 2770 if (dof != sdof) SETERRQ4(comm, PETSC_ERR_PLIB, "Coordinate dimension %d on subvertex %d, vertex %d should be %d", sdof, subvertex, vertex, dof); 2771 for (d = 0; d < dof; ++d) subCoords[soff+d] = coords[off+d]; 2772 } 2773 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 2774 ierr = VecRestoreArray(subCoordinates, &subCoords);CHKERRQ(ierr); 2775 ierr = DMSetCoordinatesLocal(subdm, subCoordinates);CHKERRQ(ierr); 2776 ierr = VecDestroy(&subCoordinates);CHKERRQ(ierr); 2777 } 2778 /* Build SF: We need this complexity because subpoints might not be selected on the owning process */ 2779 { 2780 PetscSF sfPoint, sfPointSub; 2781 IS subpIS; 2782 const PetscSFNode *remotePoints; 2783 PetscSFNode *sremotePoints, *newLocalPoints, *newOwners; 2784 const PetscInt *localPoints, *subpoints; 2785 PetscInt *slocalPoints; 2786 PetscInt numRoots, numLeaves, numSubpoints = 0, numSubroots, numSubleaves = 0, l, sl, ll, pStart, pEnd, p; 2787 PetscMPIInt rank; 2788 2789 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank);CHKERRQ(ierr); 2790 ierr = DMGetPointSF(dm, &sfPoint);CHKERRQ(ierr); 2791 ierr = DMGetPointSF(subdm, &sfPointSub);CHKERRQ(ierr); 2792 ierr = DMPlexGetChart(dm, &pStart, &pEnd);CHKERRQ(ierr); 2793 ierr = DMPlexGetChart(subdm, NULL, &numSubroots);CHKERRQ(ierr); 2794 ierr = DMPlexCreateSubpointIS(subdm, &subpIS);CHKERRQ(ierr); 2795 if (subpIS) { 2796 ierr = ISGetIndices(subpIS, &subpoints);CHKERRQ(ierr); 2797 ierr = ISGetLocalSize(subpIS, &numSubpoints);CHKERRQ(ierr); 2798 } 2799 ierr = PetscSFGetGraph(sfPoint, &numRoots, &numLeaves, &localPoints, &remotePoints);CHKERRQ(ierr); 2800 if (numRoots >= 0) { 2801 ierr = PetscMalloc2(pEnd-pStart,&newLocalPoints,numRoots,&newOwners);CHKERRQ(ierr); 2802 for (p = 0; p < pEnd-pStart; ++p) { 2803 newLocalPoints[p].rank = -2; 2804 newLocalPoints[p].index = -2; 2805 } 2806 /* Set subleaves */ 2807 for (l = 0; l < numLeaves; ++l) { 2808 const PetscInt point = localPoints[l]; 2809 const PetscInt subpoint = DMPlexFilterPoint_Internal(point, 0, numSubpoints, subpoints); 2810 2811 if (subpoint < 0) continue; 2812 newLocalPoints[point-pStart].rank = rank; 2813 newLocalPoints[point-pStart].index = subpoint; 2814 ++numSubleaves; 2815 } 2816 /* Must put in owned subpoints */ 2817 for (p = pStart; p < pEnd; ++p) { 2818 const PetscInt subpoint = DMPlexFilterPoint_Internal(p, 0, numSubpoints, subpoints); 2819 2820 if (subpoint < 0) { 2821 newOwners[p-pStart].rank = -3; 2822 newOwners[p-pStart].index = -3; 2823 } else { 2824 newOwners[p-pStart].rank = rank; 2825 newOwners[p-pStart].index = subpoint; 2826 } 2827 } 2828 ierr = PetscSFReduceBegin(sfPoint, MPIU_2INT, newLocalPoints, newOwners, MPI_MAXLOC);CHKERRQ(ierr); 2829 ierr = PetscSFReduceEnd(sfPoint, MPIU_2INT, newLocalPoints, newOwners, MPI_MAXLOC);CHKERRQ(ierr); 2830 ierr = PetscSFBcastBegin(sfPoint, MPIU_2INT, newOwners, newLocalPoints);CHKERRQ(ierr); 2831 ierr = PetscSFBcastEnd(sfPoint, MPIU_2INT, newOwners, newLocalPoints);CHKERRQ(ierr); 2832 ierr = PetscMalloc1(numSubleaves, &slocalPoints);CHKERRQ(ierr); 2833 ierr = PetscMalloc1(numSubleaves, &sremotePoints);CHKERRQ(ierr); 2834 for (l = 0, sl = 0, ll = 0; l < numLeaves; ++l) { 2835 const PetscInt point = localPoints[l]; 2836 const PetscInt subpoint = DMPlexFilterPoint_Internal(point, 0, numSubpoints, subpoints); 2837 2838 if (subpoint < 0) continue; 2839 if (newLocalPoints[point].rank == rank) {++ll; continue;} 2840 slocalPoints[sl] = subpoint; 2841 sremotePoints[sl].rank = newLocalPoints[point].rank; 2842 sremotePoints[sl].index = newLocalPoints[point].index; 2843 if (sremotePoints[sl].rank < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid remote rank for local point %d", point); 2844 if (sremotePoints[sl].index < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid remote subpoint for local point %d", point); 2845 ++sl; 2846 } 2847 if (sl + ll != numSubleaves) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Mismatch in number of subleaves %d + %d != %d", sl, ll, numSubleaves); 2848 ierr = PetscFree2(newLocalPoints,newOwners);CHKERRQ(ierr); 2849 ierr = PetscSFSetGraph(sfPointSub, numSubroots, sl, slocalPoints, PETSC_OWN_POINTER, sremotePoints, PETSC_OWN_POINTER);CHKERRQ(ierr); 2850 } 2851 if (subpIS) { 2852 ierr = ISRestoreIndices(subpIS, &subpoints);CHKERRQ(ierr); 2853 ierr = ISDestroy(&subpIS);CHKERRQ(ierr); 2854 } 2855 } 2856 /* Cleanup */ 2857 for (d = 0; d <= dim; ++d) { 2858 if (subpointIS[d]) {ierr = ISRestoreIndices(subpointIS[d], &subpoints[d]);CHKERRQ(ierr);} 2859 ierr = ISDestroy(&subpointIS[d]);CHKERRQ(ierr); 2860 } 2861 ierr = PetscFree4(numSubPoints,firstSubPoint,subpointIS,subpoints);CHKERRQ(ierr); 2862 PetscFunctionReturn(0); 2863 } 2864 2865 #undef __FUNCT__ 2866 #define __FUNCT__ "DMPlexCreateSubmesh_Interpolated" 2867 static PetscErrorCode DMPlexCreateSubmesh_Interpolated(DM dm, DMLabel vertexLabel, PetscInt value, DM subdm) 2868 { 2869 PetscErrorCode ierr; 2870 2871 PetscFunctionBegin; 2872 ierr = DMPlexCreateSubmeshGeneric_Interpolated(dm, vertexLabel, value, PETSC_FALSE, 1, subdm);CHKERRQ(ierr); 2873 PetscFunctionReturn(0); 2874 } 2875 2876 #undef __FUNCT__ 2877 #define __FUNCT__ "DMPlexCreateSubmesh" 2878 /*@ 2879 DMPlexCreateSubmesh - Extract a hypersurface from the mesh using vertices defined by a label 2880 2881 Input Parameters: 2882 + dm - The original mesh 2883 . vertexLabel - The DMLabel marking vertices contained in the surface 2884 - value - The label value to use 2885 2886 Output Parameter: 2887 . subdm - The surface mesh 2888 2889 Note: This function produces a DMLabel mapping original points in the submesh to their depth. This can be obtained using DMPlexGetSubpointMap(). 2890 2891 Level: developer 2892 2893 .seealso: DMPlexGetSubpointMap(), DMPlexGetLabel(), DMLabelSetValue() 2894 @*/ 2895 PetscErrorCode DMPlexCreateSubmesh(DM dm, DMLabel vertexLabel, PetscInt value, DM *subdm) 2896 { 2897 PetscInt dim, depth; 2898 PetscErrorCode ierr; 2899 2900 PetscFunctionBegin; 2901 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 2902 PetscValidPointer(subdm, 3); 2903 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 2904 ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr); 2905 ierr = DMCreate(PetscObjectComm((PetscObject)dm), subdm);CHKERRQ(ierr); 2906 ierr = DMSetType(*subdm, DMPLEX);CHKERRQ(ierr); 2907 ierr = DMSetDimension(*subdm, dim-1);CHKERRQ(ierr); 2908 if (depth == dim) { 2909 ierr = DMPlexCreateSubmesh_Interpolated(dm, vertexLabel, value, *subdm);CHKERRQ(ierr); 2910 } else { 2911 ierr = DMPlexCreateSubmesh_Uninterpolated(dm, vertexLabel, value, *subdm);CHKERRQ(ierr); 2912 } 2913 PetscFunctionReturn(0); 2914 } 2915 2916 #undef __FUNCT__ 2917 #define __FUNCT__ "DMPlexCreateCohesiveSubmesh_Uninterpolated" 2918 static PetscErrorCode DMPlexCreateCohesiveSubmesh_Uninterpolated(DM dm, PetscBool hasLagrange, const char label[], PetscInt value, DM subdm) 2919 { 2920 MPI_Comm comm; 2921 DMLabel subpointMap; 2922 IS subvertexIS; 2923 const PetscInt *subVertices; 2924 PetscInt numSubVertices, firstSubVertex, numSubCells, *subCells = NULL; 2925 PetscInt *subface, maxConeSize, numSubFaces, firstSubFace, newFacePoint, nFV; 2926 PetscInt cMax, c, f; 2927 PetscErrorCode ierr; 2928 2929 PetscFunctionBegin; 2930 ierr = PetscObjectGetComm((PetscObject)dm, &comm);CHKERRQ(ierr); 2931 /* Create subpointMap which marks the submesh */ 2932 ierr = DMLabelCreate("subpoint_map", &subpointMap);CHKERRQ(ierr); 2933 ierr = DMPlexSetSubpointMap(subdm, subpointMap);CHKERRQ(ierr); 2934 ierr = DMLabelDestroy(&subpointMap);CHKERRQ(ierr); 2935 ierr = DMPlexMarkCohesiveSubmesh_Uninterpolated(dm, hasLagrange, label, value, subpointMap, &numSubFaces, &nFV, &subCells, subdm);CHKERRQ(ierr); 2936 /* Setup chart */ 2937 ierr = DMLabelGetStratumSize(subpointMap, 0, &numSubVertices);CHKERRQ(ierr); 2938 ierr = DMLabelGetStratumSize(subpointMap, 2, &numSubCells);CHKERRQ(ierr); 2939 ierr = DMPlexSetChart(subdm, 0, numSubCells+numSubFaces+numSubVertices);CHKERRQ(ierr); 2940 ierr = DMPlexSetVTKCellHeight(subdm, 1);CHKERRQ(ierr); 2941 /* Set cone sizes */ 2942 firstSubVertex = numSubCells; 2943 firstSubFace = numSubCells+numSubVertices; 2944 newFacePoint = firstSubFace; 2945 ierr = DMLabelGetStratumIS(subpointMap, 0, &subvertexIS);CHKERRQ(ierr); 2946 if (subvertexIS) {ierr = ISGetIndices(subvertexIS, &subVertices);CHKERRQ(ierr);} 2947 for (c = 0; c < numSubCells; ++c) { 2948 ierr = DMPlexSetConeSize(subdm, c, 1);CHKERRQ(ierr); 2949 } 2950 for (f = firstSubFace; f < firstSubFace+numSubFaces; ++f) { 2951 ierr = DMPlexSetConeSize(subdm, f, nFV);CHKERRQ(ierr); 2952 } 2953 ierr = DMSetUp(subdm);CHKERRQ(ierr); 2954 /* Create face cones */ 2955 ierr = DMPlexGetMaxSizes(dm, &maxConeSize, NULL);CHKERRQ(ierr); 2956 ierr = DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);CHKERRQ(ierr); 2957 ierr = DMGetWorkArray(subdm, maxConeSize, PETSC_INT, (void**) &subface);CHKERRQ(ierr); 2958 for (c = 0; c < numSubCells; ++c) { 2959 const PetscInt cell = subCells[c]; 2960 const PetscInt subcell = c; 2961 const PetscInt *cone, *cells; 2962 PetscInt numCells, subVertex, p, v; 2963 2964 if (cell < cMax) continue; 2965 ierr = DMPlexGetCone(dm, cell, &cone);CHKERRQ(ierr); 2966 for (v = 0; v < nFV; ++v) { 2967 ierr = PetscFindInt(cone[v], numSubVertices, subVertices, &subVertex);CHKERRQ(ierr); 2968 subface[v] = firstSubVertex+subVertex; 2969 } 2970 ierr = DMPlexSetCone(subdm, newFacePoint, subface);CHKERRQ(ierr); 2971 ierr = DMPlexSetCone(subdm, subcell, &newFacePoint);CHKERRQ(ierr); 2972 ierr = DMPlexGetJoin(dm, nFV, cone, &numCells, &cells);CHKERRQ(ierr); 2973 /* Not true in parallel 2974 if (numCells != 2) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cohesive cells should separate two cells"); */ 2975 for (p = 0; p < numCells; ++p) { 2976 PetscInt negsubcell; 2977 2978 if (cells[p] >= cMax) continue; 2979 /* I know this is a crap search */ 2980 for (negsubcell = 0; negsubcell < numSubCells; ++negsubcell) { 2981 if (subCells[negsubcell] == cells[p]) break; 2982 } 2983 if (negsubcell == numSubCells) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Could not find negative face neighbor for cohesive cell %d", cell); 2984 ierr = DMPlexSetCone(subdm, negsubcell, &newFacePoint);CHKERRQ(ierr); 2985 } 2986 ierr = DMPlexRestoreJoin(dm, nFV, cone, &numCells, &cells);CHKERRQ(ierr); 2987 ++newFacePoint; 2988 } 2989 ierr = DMRestoreWorkArray(subdm, maxConeSize, PETSC_INT, (void**) &subface);CHKERRQ(ierr); 2990 ierr = DMPlexSymmetrize(subdm);CHKERRQ(ierr); 2991 ierr = DMPlexStratify(subdm);CHKERRQ(ierr); 2992 /* Build coordinates */ 2993 { 2994 PetscSection coordSection, subCoordSection; 2995 Vec coordinates, subCoordinates; 2996 PetscScalar *coords, *subCoords; 2997 PetscInt numComp, coordSize, v; 2998 2999 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 3000 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr); 3001 ierr = DMGetCoordinateSection(subdm, &subCoordSection);CHKERRQ(ierr); 3002 ierr = PetscSectionSetNumFields(subCoordSection, 1);CHKERRQ(ierr); 3003 ierr = PetscSectionGetFieldComponents(coordSection, 0, &numComp);CHKERRQ(ierr); 3004 ierr = PetscSectionSetFieldComponents(subCoordSection, 0, numComp);CHKERRQ(ierr); 3005 ierr = PetscSectionSetChart(subCoordSection, firstSubVertex, firstSubVertex+numSubVertices);CHKERRQ(ierr); 3006 for (v = 0; v < numSubVertices; ++v) { 3007 const PetscInt vertex = subVertices[v]; 3008 const PetscInt subvertex = firstSubVertex+v; 3009 PetscInt dof; 3010 3011 ierr = PetscSectionGetDof(coordSection, vertex, &dof);CHKERRQ(ierr); 3012 ierr = PetscSectionSetDof(subCoordSection, subvertex, dof);CHKERRQ(ierr); 3013 ierr = PetscSectionSetFieldDof(subCoordSection, subvertex, 0, dof);CHKERRQ(ierr); 3014 } 3015 ierr = PetscSectionSetUp(subCoordSection);CHKERRQ(ierr); 3016 ierr = PetscSectionGetStorageSize(subCoordSection, &coordSize);CHKERRQ(ierr); 3017 ierr = VecCreate(comm, &subCoordinates);CHKERRQ(ierr); 3018 ierr = VecSetSizes(subCoordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 3019 ierr = VecSetType(subCoordinates,VECSTANDARD);CHKERRQ(ierr); 3020 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 3021 ierr = VecGetArray(subCoordinates, &subCoords);CHKERRQ(ierr); 3022 for (v = 0; v < numSubVertices; ++v) { 3023 const PetscInt vertex = subVertices[v]; 3024 const PetscInt subvertex = firstSubVertex+v; 3025 PetscInt dof, off, sdof, soff, d; 3026 3027 ierr = PetscSectionGetDof(coordSection, vertex, &dof);CHKERRQ(ierr); 3028 ierr = PetscSectionGetOffset(coordSection, vertex, &off);CHKERRQ(ierr); 3029 ierr = PetscSectionGetDof(subCoordSection, subvertex, &sdof);CHKERRQ(ierr); 3030 ierr = PetscSectionGetOffset(subCoordSection, subvertex, &soff);CHKERRQ(ierr); 3031 if (dof != sdof) SETERRQ4(comm, PETSC_ERR_PLIB, "Coordinate dimension %d on subvertex %d, vertex %d should be %d", sdof, subvertex, vertex, dof); 3032 for (d = 0; d < dof; ++d) subCoords[soff+d] = coords[off+d]; 3033 } 3034 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 3035 ierr = VecRestoreArray(subCoordinates, &subCoords);CHKERRQ(ierr); 3036 ierr = DMSetCoordinatesLocal(subdm, subCoordinates);CHKERRQ(ierr); 3037 ierr = VecDestroy(&subCoordinates);CHKERRQ(ierr); 3038 } 3039 /* Build SF */ 3040 CHKMEMQ; 3041 { 3042 PetscSF sfPoint, sfPointSub; 3043 const PetscSFNode *remotePoints; 3044 PetscSFNode *sremotePoints, *newLocalPoints, *newOwners; 3045 const PetscInt *localPoints; 3046 PetscInt *slocalPoints; 3047 PetscInt numRoots, numLeaves, numSubRoots = numSubCells+numSubFaces+numSubVertices, numSubLeaves = 0, l, sl, ll, pStart, pEnd, p, vStart, vEnd; 3048 PetscMPIInt rank; 3049 3050 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank);CHKERRQ(ierr); 3051 ierr = DMGetPointSF(dm, &sfPoint);CHKERRQ(ierr); 3052 ierr = DMGetPointSF(subdm, &sfPointSub);CHKERRQ(ierr); 3053 ierr = DMPlexGetChart(dm, &pStart, &pEnd);CHKERRQ(ierr); 3054 ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr); 3055 ierr = PetscSFGetGraph(sfPoint, &numRoots, &numLeaves, &localPoints, &remotePoints);CHKERRQ(ierr); 3056 if (numRoots >= 0) { 3057 /* Only vertices should be shared */ 3058 ierr = PetscMalloc2(pEnd-pStart,&newLocalPoints,numRoots,&newOwners);CHKERRQ(ierr); 3059 for (p = 0; p < pEnd-pStart; ++p) { 3060 newLocalPoints[p].rank = -2; 3061 newLocalPoints[p].index = -2; 3062 } 3063 /* Set subleaves */ 3064 for (l = 0; l < numLeaves; ++l) { 3065 const PetscInt point = localPoints[l]; 3066 const PetscInt subPoint = DMPlexFilterPoint_Internal(point, firstSubVertex, numSubVertices, subVertices); 3067 3068 if ((point < vStart) && (point >= vEnd)) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Should not be mapping anything but vertices, %d", point); 3069 if (subPoint < 0) continue; 3070 newLocalPoints[point-pStart].rank = rank; 3071 newLocalPoints[point-pStart].index = subPoint; 3072 ++numSubLeaves; 3073 } 3074 /* Must put in owned subpoints */ 3075 for (p = pStart; p < pEnd; ++p) { 3076 const PetscInt subPoint = DMPlexFilterPoint_Internal(p, firstSubVertex, numSubVertices, subVertices); 3077 3078 if (subPoint < 0) { 3079 newOwners[p-pStart].rank = -3; 3080 newOwners[p-pStart].index = -3; 3081 } else { 3082 newOwners[p-pStart].rank = rank; 3083 newOwners[p-pStart].index = subPoint; 3084 } 3085 } 3086 ierr = PetscSFReduceBegin(sfPoint, MPIU_2INT, newLocalPoints, newOwners, MPI_MAXLOC);CHKERRQ(ierr); 3087 ierr = PetscSFReduceEnd(sfPoint, MPIU_2INT, newLocalPoints, newOwners, MPI_MAXLOC);CHKERRQ(ierr); 3088 ierr = PetscSFBcastBegin(sfPoint, MPIU_2INT, newOwners, newLocalPoints);CHKERRQ(ierr); 3089 ierr = PetscSFBcastEnd(sfPoint, MPIU_2INT, newOwners, newLocalPoints);CHKERRQ(ierr); 3090 ierr = PetscMalloc1(numSubLeaves, &slocalPoints);CHKERRQ(ierr); 3091 ierr = PetscMalloc1(numSubLeaves, &sremotePoints);CHKERRQ(ierr); 3092 for (l = 0, sl = 0, ll = 0; l < numLeaves; ++l) { 3093 const PetscInt point = localPoints[l]; 3094 const PetscInt subPoint = DMPlexFilterPoint_Internal(point, firstSubVertex, numSubVertices, subVertices); 3095 3096 if (subPoint < 0) continue; 3097 if (newLocalPoints[point].rank == rank) {++ll; continue;} 3098 slocalPoints[sl] = subPoint; 3099 sremotePoints[sl].rank = newLocalPoints[point].rank; 3100 sremotePoints[sl].index = newLocalPoints[point].index; 3101 if (sremotePoints[sl].rank < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid remote rank for local point %d", point); 3102 if (sremotePoints[sl].index < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid remote subpoint for local point %d", point); 3103 ++sl; 3104 } 3105 ierr = PetscFree2(newLocalPoints,newOwners);CHKERRQ(ierr); 3106 if (sl + ll != numSubLeaves) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Mismatch in number of subleaves %d + %d != %d", sl, ll, numSubLeaves); 3107 ierr = PetscSFSetGraph(sfPointSub, numSubRoots, sl, slocalPoints, PETSC_OWN_POINTER, sremotePoints, PETSC_OWN_POINTER);CHKERRQ(ierr); 3108 } 3109 } 3110 CHKMEMQ; 3111 /* Cleanup */ 3112 if (subvertexIS) {ierr = ISRestoreIndices(subvertexIS, &subVertices);CHKERRQ(ierr);} 3113 ierr = ISDestroy(&subvertexIS);CHKERRQ(ierr); 3114 ierr = PetscFree(subCells);CHKERRQ(ierr); 3115 PetscFunctionReturn(0); 3116 } 3117 3118 #undef __FUNCT__ 3119 #define __FUNCT__ "DMPlexCreateCohesiveSubmesh_Interpolated" 3120 static PetscErrorCode DMPlexCreateCohesiveSubmesh_Interpolated(DM dm, const char labelname[], PetscInt value, DM subdm) 3121 { 3122 DMLabel label = NULL; 3123 PetscErrorCode ierr; 3124 3125 PetscFunctionBegin; 3126 if (labelname) {ierr = DMPlexGetLabel(dm, labelname, &label);CHKERRQ(ierr);} 3127 ierr = DMPlexCreateSubmeshGeneric_Interpolated(dm, label, value, PETSC_TRUE, 1, subdm);CHKERRQ(ierr); 3128 PetscFunctionReturn(0); 3129 } 3130 3131 #undef __FUNCT__ 3132 #define __FUNCT__ "DMPlexCreateCohesiveSubmesh" 3133 /* 3134 DMPlexCreateCohesiveSubmesh - Extract from a mesh with cohesive cells the hypersurface defined by one face of the cells. Optionally, a Label an be given to restrict the cells. 3135 3136 Input Parameters: 3137 + dm - The original mesh 3138 . hasLagrange - The mesh has Lagrange unknowns in the cohesive cells 3139 . label - A label name, or NULL 3140 - value - A label value 3141 3142 Output Parameter: 3143 . subdm - The surface mesh 3144 3145 Note: This function produces a DMLabel mapping original points in the submesh to their depth. This can be obtained using DMPlexGetSubpointMap(). 3146 3147 Level: developer 3148 3149 .seealso: DMPlexGetSubpointMap(), DMPlexCreateSubmesh() 3150 */ 3151 PetscErrorCode DMPlexCreateCohesiveSubmesh(DM dm, PetscBool hasLagrange, const char label[], PetscInt value, DM *subdm) 3152 { 3153 PetscInt dim, depth; 3154 PetscErrorCode ierr; 3155 3156 PetscFunctionBegin; 3157 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 3158 PetscValidPointer(subdm, 5); 3159 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 3160 ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr); 3161 ierr = DMCreate(PetscObjectComm((PetscObject)dm), subdm);CHKERRQ(ierr); 3162 ierr = DMSetType(*subdm, DMPLEX);CHKERRQ(ierr); 3163 ierr = DMSetDimension(*subdm, dim-1);CHKERRQ(ierr); 3164 if (depth == dim) { 3165 ierr = DMPlexCreateCohesiveSubmesh_Interpolated(dm, label, value, *subdm);CHKERRQ(ierr); 3166 } else { 3167 ierr = DMPlexCreateCohesiveSubmesh_Uninterpolated(dm, hasLagrange, label, value, *subdm);CHKERRQ(ierr); 3168 } 3169 PetscFunctionReturn(0); 3170 } 3171 3172 #undef __FUNCT__ 3173 #define __FUNCT__ "DMPlexFilter" 3174 /*@ 3175 DMPlexFilter - Extract a subset of mesh cells defined by a label as a separate mesh 3176 3177 Input Parameters: 3178 + dm - The original mesh 3179 . cellLabel - The DMLabel marking cells contained in the new mesh 3180 - value - The label value to use 3181 3182 Output Parameter: 3183 . subdm - The new mesh 3184 3185 Note: This function produces a DMLabel mapping original points in the submesh to their depth. This can be obtained using DMPlexGetSubpointMap(). 3186 3187 Level: developer 3188 3189 .seealso: DMPlexGetSubpointMap(), DMPlexGetLabel(), DMLabelSetValue() 3190 @*/ 3191 PetscErrorCode DMPlexFilter(DM dm, DMLabel cellLabel, PetscInt value, DM *subdm) 3192 { 3193 PetscInt dim; 3194 PetscErrorCode ierr; 3195 3196 PetscFunctionBegin; 3197 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 3198 PetscValidPointer(subdm, 3); 3199 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 3200 ierr = DMCreate(PetscObjectComm((PetscObject) dm), subdm);CHKERRQ(ierr); 3201 ierr = DMSetType(*subdm, DMPLEX);CHKERRQ(ierr); 3202 ierr = DMSetDimension(*subdm, dim);CHKERRQ(ierr); 3203 /* Extract submesh in place, could be empty on some procs, could have inconsistency if procs do not both extract a shared cell */ 3204 ierr = DMPlexCreateSubmeshGeneric_Interpolated(dm, cellLabel, value, PETSC_FALSE, 0, *subdm);CHKERRQ(ierr); 3205 PetscFunctionReturn(0); 3206 } 3207 3208 #undef __FUNCT__ 3209 #define __FUNCT__ "DMPlexGetSubpointMap" 3210 /*@ 3211 DMPlexGetSubpointMap - Returns a DMLabel with point dimension as values 3212 3213 Input Parameter: 3214 . dm - The submesh DM 3215 3216 Output Parameter: 3217 . subpointMap - The DMLabel of all the points from the original mesh in this submesh, or NULL if this is not a submesh 3218 3219 Level: developer 3220 3221 .seealso: DMPlexCreateSubmesh(), DMPlexCreateSubpointIS() 3222 @*/ 3223 PetscErrorCode DMPlexGetSubpointMap(DM dm, DMLabel *subpointMap) 3224 { 3225 PetscFunctionBegin; 3226 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 3227 PetscValidPointer(subpointMap, 2); 3228 *subpointMap = ((DM_Plex*) dm->data)->subpointMap; 3229 PetscFunctionReturn(0); 3230 } 3231 3232 #undef __FUNCT__ 3233 #define __FUNCT__ "DMPlexSetSubpointMap" 3234 /* Note: Should normally not be called by the user, since it is set in DMPlexCreateSubmesh() */ 3235 PetscErrorCode DMPlexSetSubpointMap(DM dm, DMLabel subpointMap) 3236 { 3237 DM_Plex *mesh = (DM_Plex *) dm->data; 3238 DMLabel tmp; 3239 PetscErrorCode ierr; 3240 3241 PetscFunctionBegin; 3242 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 3243 tmp = mesh->subpointMap; 3244 mesh->subpointMap = subpointMap; 3245 ++mesh->subpointMap->refct; 3246 ierr = DMLabelDestroy(&tmp);CHKERRQ(ierr); 3247 PetscFunctionReturn(0); 3248 } 3249 3250 #undef __FUNCT__ 3251 #define __FUNCT__ "DMPlexCreateSubpointIS" 3252 /*@ 3253 DMPlexCreateSubpointIS - Creates an IS covering the entire subdm chart with the original points as data 3254 3255 Input Parameter: 3256 . dm - The submesh DM 3257 3258 Output Parameter: 3259 . subpointIS - The IS of all the points from the original mesh in this submesh, or NULL if this is not a submesh 3260 3261 Note: This IS is guaranteed to be sorted by the construction of the submesh 3262 3263 Level: developer 3264 3265 .seealso: DMPlexCreateSubmesh(), DMPlexGetSubpointMap() 3266 @*/ 3267 PetscErrorCode DMPlexCreateSubpointIS(DM dm, IS *subpointIS) 3268 { 3269 MPI_Comm comm; 3270 DMLabel subpointMap; 3271 IS is; 3272 const PetscInt *opoints; 3273 PetscInt *points, *depths; 3274 PetscInt depth, depStart, depEnd, d, pStart, pEnd, p, n, off; 3275 PetscErrorCode ierr; 3276 3277 PetscFunctionBegin; 3278 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 3279 PetscValidPointer(subpointIS, 2); 3280 ierr = PetscObjectGetComm((PetscObject)dm,&comm);CHKERRQ(ierr); 3281 *subpointIS = NULL; 3282 ierr = DMPlexGetSubpointMap(dm, &subpointMap);CHKERRQ(ierr); 3283 ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr); 3284 if (subpointMap && depth >= 0) { 3285 ierr = DMPlexGetChart(dm, &pStart, &pEnd);CHKERRQ(ierr); 3286 if (pStart) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Submeshes must start the point numbering at 0, not %d", pStart); 3287 ierr = DMGetWorkArray(dm, depth+1, PETSC_INT, &depths);CHKERRQ(ierr); 3288 depths[0] = depth; 3289 depths[1] = 0; 3290 for(d = 2; d <= depth; ++d) {depths[d] = depth+1 - d;} 3291 ierr = PetscMalloc1(pEnd, &points);CHKERRQ(ierr); 3292 for(d = 0, off = 0; d <= depth; ++d) { 3293 const PetscInt dep = depths[d]; 3294 3295 ierr = DMPlexGetDepthStratum(dm, dep, &depStart, &depEnd);CHKERRQ(ierr); 3296 ierr = DMLabelGetStratumSize(subpointMap, dep, &n);CHKERRQ(ierr); 3297 if (((d < 2) && (depth > 1)) || (d == 1)) { /* Only check vertices and cells for now since the map is broken for others */ 3298 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); 3299 } else { 3300 if (!n) { 3301 if (d == 0) { 3302 /* Missing cells */ 3303 for(p = 0; p < depEnd-depStart; ++p, ++off) points[off] = -1; 3304 } else { 3305 /* Missing faces */ 3306 for(p = 0; p < depEnd-depStart; ++p, ++off) points[off] = PETSC_MAX_INT; 3307 } 3308 } 3309 } 3310 if (n) { 3311 ierr = DMLabelGetStratumIS(subpointMap, dep, &is);CHKERRQ(ierr); 3312 ierr = ISGetIndices(is, &opoints);CHKERRQ(ierr); 3313 for(p = 0; p < n; ++p, ++off) points[off] = opoints[p]; 3314 ierr = ISRestoreIndices(is, &opoints);CHKERRQ(ierr); 3315 ierr = ISDestroy(&is);CHKERRQ(ierr); 3316 } 3317 } 3318 ierr = DMRestoreWorkArray(dm, depth+1, PETSC_INT, &depths);CHKERRQ(ierr); 3319 if (off != pEnd) SETERRQ2(comm, PETSC_ERR_ARG_WRONG, "The number of mapped submesh points %d should be %d", off, pEnd); 3320 ierr = ISCreateGeneral(PETSC_COMM_SELF, pEnd, points, PETSC_OWN_POINTER, subpointIS);CHKERRQ(ierr); 3321 } 3322 PetscFunctionReturn(0); 3323 } 3324