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