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