1 #include <petsc-private/dmpleximpl.h> /*I "petscdmplex.h" I*/ 2 #include <petscsf.h> 3 4 #undef __FUNCT__ 5 #define __FUNCT__ "DMPlexMarkBoundaryFaces_Internal" 6 PetscErrorCode DMPlexMarkBoundaryFaces_Internal(DM dm, PetscInt cellHeight, DMLabel label) 7 { 8 PetscInt fStart, fEnd, f; 9 PetscErrorCode ierr; 10 11 PetscFunctionBegin; 12 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 13 ierr = DMPlexGetHeightStratum(dm, cellHeight+1, &fStart, &fEnd);CHKERRQ(ierr); 14 for (f = fStart; f < fEnd; ++f) { 15 PetscInt supportSize; 16 17 ierr = DMPlexGetSupportSize(dm, f, &supportSize);CHKERRQ(ierr); 18 if (supportSize == 1) {ierr = DMLabelSetValue(label, f, 1);CHKERRQ(ierr);} 19 } 20 PetscFunctionReturn(0); 21 } 22 23 #undef __FUNCT__ 24 #define __FUNCT__ "DMPlexMarkBoundaryFaces" 25 /*@ 26 DMPlexMarkBoundaryFaces - Mark all faces on the boundary 27 28 Not Collective 29 30 Input Parameter: 31 . dm - The original DM 32 33 Output Parameter: 34 . label - The DMLabel marking boundary faces with value 1 35 36 Level: developer 37 38 .seealso: DMLabelCreate(), DMPlexCreateLabel() 39 @*/ 40 PetscErrorCode DMPlexMarkBoundaryFaces(DM dm, DMLabel label) 41 { 42 PetscErrorCode ierr; 43 44 PetscFunctionBegin; 45 ierr = DMPlexMarkBoundaryFaces_Internal(dm, 0, label);CHKERRQ(ierr); 46 PetscFunctionReturn(0); 47 } 48 49 #undef __FUNCT__ 50 #define __FUNCT__ "DMPlexLabelComplete" 51 /*@ 52 DMPlexLabelComplete - Starting with a label marking points on a surface, we add the transitive closure to the surface 53 54 Input Parameters: 55 + dm - The DM 56 - label - A DMLabel marking the surface points 57 58 Output Parameter: 59 . label - A DMLabel marking all surface points in the transitive closure 60 61 Level: developer 62 63 .seealso: DMPlexLabelCohesiveComplete() 64 @*/ 65 PetscErrorCode DMPlexLabelComplete(DM dm, DMLabel label) 66 { 67 IS valueIS; 68 const PetscInt *values; 69 PetscInt numValues, v; 70 PetscErrorCode ierr; 71 72 PetscFunctionBegin; 73 ierr = DMLabelGetNumValues(label, &numValues);CHKERRQ(ierr); 74 ierr = DMLabelGetValueIS(label, &valueIS);CHKERRQ(ierr); 75 ierr = ISGetIndices(valueIS, &values);CHKERRQ(ierr); 76 for (v = 0; v < numValues; ++v) { 77 IS pointIS; 78 const PetscInt *points; 79 PetscInt numPoints, p; 80 81 ierr = DMLabelGetStratumSize(label, values[v], &numPoints);CHKERRQ(ierr); 82 ierr = DMLabelGetStratumIS(label, values[v], &pointIS);CHKERRQ(ierr); 83 ierr = ISGetIndices(pointIS, &points);CHKERRQ(ierr); 84 for (p = 0; p < numPoints; ++p) { 85 PetscInt *closure = NULL; 86 PetscInt closureSize, c; 87 88 ierr = DMPlexGetTransitiveClosure(dm, points[p], PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 89 for (c = 0; c < closureSize*2; c += 2) { 90 ierr = DMLabelSetValue(label, closure[c], values[v]);CHKERRQ(ierr); 91 } 92 ierr = DMPlexRestoreTransitiveClosure(dm, points[p], PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 93 } 94 ierr = ISRestoreIndices(pointIS, &points);CHKERRQ(ierr); 95 ierr = ISDestroy(&pointIS);CHKERRQ(ierr); 96 } 97 ierr = ISRestoreIndices(valueIS, &values);CHKERRQ(ierr); 98 ierr = ISDestroy(&valueIS);CHKERRQ(ierr); 99 PetscFunctionReturn(0); 100 } 101 102 #undef __FUNCT__ 103 #define __FUNCT__ "DMPlexLabelAddCells" 104 /*@ 105 DMPlexLabelAddCells - Starting with a label marking faces on a surface, we add a cell for each face 106 107 Input Parameters: 108 + dm - The DM 109 - label - A DMLabel marking the surface points 110 111 Output Parameter: 112 . label - A DMLabel incorporating cells 113 114 Level: developer 115 116 Note: The cells allow FEM boundary conditions to be applied using the cell geometry 117 118 .seealso: DMPlexLabelComplete(), DMPlexLabelCohesiveComplete() 119 @*/ 120 PetscErrorCode DMPlexLabelAddCells(DM dm, DMLabel label) 121 { 122 IS valueIS; 123 const PetscInt *values; 124 PetscInt numValues, v, cStart, cEnd; 125 PetscErrorCode ierr; 126 127 PetscFunctionBegin; 128 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 129 ierr = DMLabelGetNumValues(label, &numValues);CHKERRQ(ierr); 130 ierr = DMLabelGetValueIS(label, &valueIS);CHKERRQ(ierr); 131 ierr = ISGetIndices(valueIS, &values);CHKERRQ(ierr); 132 for (v = 0; v < numValues; ++v) { 133 IS pointIS; 134 const PetscInt *points; 135 PetscInt numPoints, p; 136 137 ierr = DMLabelGetStratumSize(label, values[v], &numPoints);CHKERRQ(ierr); 138 ierr = DMLabelGetStratumIS(label, values[v], &pointIS);CHKERRQ(ierr); 139 ierr = ISGetIndices(pointIS, &points);CHKERRQ(ierr); 140 for (p = 0; p < numPoints; ++p) { 141 PetscInt *closure = NULL; 142 PetscInt closureSize, point; 143 144 ierr = DMPlexGetTransitiveClosure(dm, points[p], PETSC_FALSE, &closureSize, &closure);CHKERRQ(ierr); 145 point = closure[(closureSize-1)*2]; 146 if ((point >= cStart) && (point < cEnd)) {ierr = DMLabelSetValue(label, point, values[v]);CHKERRQ(ierr);} 147 ierr = DMPlexRestoreTransitiveClosure(dm, points[p], PETSC_FALSE, &closureSize, &closure);CHKERRQ(ierr); 148 } 149 ierr = ISRestoreIndices(pointIS, &points);CHKERRQ(ierr); 150 ierr = ISDestroy(&pointIS);CHKERRQ(ierr); 151 } 152 ierr = ISRestoreIndices(valueIS, &values);CHKERRQ(ierr); 153 ierr = ISDestroy(&valueIS);CHKERRQ(ierr); 154 PetscFunctionReturn(0); 155 } 156 157 #undef __FUNCT__ 158 #define __FUNCT__ "DMPlexShiftPoint_Internal" 159 PETSC_STATIC_INLINE PetscInt DMPlexShiftPoint_Internal(PetscInt p, PetscInt depth, PetscInt depthMax[], PetscInt depthEnd[], PetscInt depthShift[]) 160 { 161 if (depth < 0) return p; 162 /* Normal Cells */ if (p < depthMax[depth]) return p; 163 /* Hybrid Cells+Normal Vertices */ if (p < depthMax[0]) return p + depthShift[depth]; 164 /* Hybrid Vertices+Normal Faces */ if (depth < 2 || p < depthMax[depth-1]) return p + depthShift[depth] + depthShift[0]; 165 /* Hybrid Faces+Normal Edges */ if (depth < 3 || p < depthMax[depth-2]) return p + depthShift[depth] + depthShift[0] + depthShift[depth-1]; 166 /* Hybrid Edges */ return p + depthShift[depth] + depthShift[0] + depthShift[depth-1] + depthShift[depth-2]; 167 } 168 169 #undef __FUNCT__ 170 #define __FUNCT__ "DMPlexShiftSizes_Internal" 171 static PetscErrorCode DMPlexShiftSizes_Internal(DM dm, PetscInt depthShift[], DM dmNew) 172 { 173 PetscInt *depthMax, *depthEnd; 174 PetscInt depth = 0, d, pStart, pEnd, p; 175 PetscErrorCode ierr; 176 177 PetscFunctionBegin; 178 ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr); 179 if (depth < 0) PetscFunctionReturn(0); 180 ierr = PetscMalloc2(depth+1,&depthMax,depth+1,&depthEnd);CHKERRQ(ierr); 181 /* Step 1: Expand chart */ 182 ierr = DMPlexGetChart(dm, &pStart, &pEnd);CHKERRQ(ierr); 183 ierr = DMPlexGetHybridBounds(dm, &depthMax[depth], depth > 0 ? &depthMax[depth-1] : NULL, &depthMax[1], &depthMax[0]);CHKERRQ(ierr); 184 for (d = 0; d <= depth; ++d) { 185 pEnd += depthShift[d]; 186 ierr = DMPlexGetDepthStratum(dm, d, NULL, &depthEnd[d]);CHKERRQ(ierr); 187 depthMax[d] = depthMax[d] < 0 ? depthEnd[d] : depthMax[d]; 188 } 189 ierr = DMPlexSetChart(dmNew, pStart, pEnd);CHKERRQ(ierr); 190 /* Step 2: Set cone and support sizes */ 191 for (d = 0; d <= depth; ++d) { 192 ierr = DMPlexGetDepthStratum(dm, d, &pStart, &pEnd);CHKERRQ(ierr); 193 for (p = pStart; p < pEnd; ++p) { 194 PetscInt newp = DMPlexShiftPoint_Internal(p, depth, depthMax, depthEnd, depthShift); 195 PetscInt size; 196 197 ierr = DMPlexGetConeSize(dm, p, &size);CHKERRQ(ierr); 198 ierr = DMPlexSetConeSize(dmNew, newp, size);CHKERRQ(ierr); 199 ierr = DMPlexGetSupportSize(dm, p, &size);CHKERRQ(ierr); 200 ierr = DMPlexSetSupportSize(dmNew, newp, size);CHKERRQ(ierr); 201 } 202 } 203 ierr = PetscFree2(depthMax,depthEnd);CHKERRQ(ierr); 204 PetscFunctionReturn(0); 205 } 206 207 #undef __FUNCT__ 208 #define __FUNCT__ "DMPlexShiftPoints_Internal" 209 static PetscErrorCode DMPlexShiftPoints_Internal(DM dm, PetscInt depthShift[], DM dmNew) 210 { 211 PetscInt *depthEnd, *depthMax, *newpoints; 212 PetscInt depth = 0, d, maxConeSize, maxSupportSize, maxConeSizeNew, maxSupportSizeNew, pStart, pEnd, p; 213 PetscErrorCode ierr; 214 215 PetscFunctionBegin; 216 ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr); 217 if (depth < 0) PetscFunctionReturn(0); 218 ierr = DMPlexGetMaxSizes(dm, &maxConeSize, &maxSupportSize);CHKERRQ(ierr); 219 ierr = DMPlexGetMaxSizes(dmNew, &maxConeSizeNew, &maxSupportSizeNew);CHKERRQ(ierr); 220 ierr = PetscMalloc3(depth+1,&depthEnd,depth+1,&depthMax,PetscMax(PetscMax(maxConeSize, maxSupportSize), PetscMax(maxConeSizeNew, maxSupportSizeNew)),&newpoints);CHKERRQ(ierr); 221 ierr = DMPlexGetHybridBounds(dm, &depthMax[depth], depth > 0 ? &depthMax[depth-1] : NULL, &depthMax[1], &depthMax[0]);CHKERRQ(ierr); 222 for (d = 0; d <= depth; ++d) { 223 ierr = DMPlexGetDepthStratum(dm, d, NULL, &depthEnd[d]);CHKERRQ(ierr); 224 depthMax[d] = depthMax[d] < 0 ? depthEnd[d] : depthMax[d]; 225 } 226 /* Step 5: Set cones and supports */ 227 ierr = DMPlexGetChart(dm, &pStart, &pEnd);CHKERRQ(ierr); 228 for (p = pStart; p < pEnd; ++p) { 229 const PetscInt *points = NULL, *orientations = NULL; 230 PetscInt size,sizeNew, i, newp = DMPlexShiftPoint_Internal(p, depth, depthMax, depthEnd, depthShift); 231 232 ierr = DMPlexGetConeSize(dm, p, &size);CHKERRQ(ierr); 233 ierr = DMPlexGetCone(dm, p, &points);CHKERRQ(ierr); 234 ierr = DMPlexGetConeOrientation(dm, p, &orientations);CHKERRQ(ierr); 235 for (i = 0; i < size; ++i) { 236 newpoints[i] = DMPlexShiftPoint_Internal(points[i], depth, depthMax, depthEnd, depthShift); 237 } 238 ierr = DMPlexSetCone(dmNew, newp, newpoints);CHKERRQ(ierr); 239 ierr = DMPlexSetConeOrientation(dmNew, newp, orientations);CHKERRQ(ierr); 240 ierr = DMPlexGetSupportSize(dm, p, &size);CHKERRQ(ierr); 241 ierr = DMPlexGetSupportSize(dmNew, newp, &sizeNew);CHKERRQ(ierr); 242 ierr = DMPlexGetSupport(dm, p, &points);CHKERRQ(ierr); 243 for (i = 0; i < size; ++i) { 244 newpoints[i] = DMPlexShiftPoint_Internal(points[i], depth, depthMax, depthEnd, depthShift); 245 } 246 for (i = size; i < sizeNew; ++i) newpoints[i] = 0; 247 ierr = DMPlexSetSupport(dmNew, newp, newpoints);CHKERRQ(ierr); 248 } 249 ierr = PetscFree3(depthEnd,depthMax,newpoints);CHKERRQ(ierr); 250 PetscFunctionReturn(0); 251 } 252 253 #undef __FUNCT__ 254 #define __FUNCT__ "DMPlexShiftCoordinates_Internal" 255 static PetscErrorCode DMPlexShiftCoordinates_Internal(DM dm, PetscInt depthShift[], DM dmNew) 256 { 257 PetscSection coordSection, newCoordSection; 258 Vec coordinates, newCoordinates; 259 PetscScalar *coords, *newCoords; 260 PetscInt *depthEnd, coordSize; 261 PetscInt dim, depth = 0, d, vStart, vEnd, vStartNew, vEndNew, v; 262 PetscErrorCode ierr; 263 264 PetscFunctionBegin; 265 ierr = 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, crossEdgeIS = NULL; 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 PetscInt bval = -1; 1428 1429 ierr = DMLabelGetValue(label, clp, &val);CHKERRQ(ierr); 1430 if (blabel) {ierr = DMLabelGetValue(blabel, clp, &bval);CHKERRQ(ierr);} 1431 if ((val >= 0) && (val < dim-1) && (bval < 0)) { 1432 ierr = DMLabelSetValue(label, point, pos == PETSC_TRUE ? shift+dim-1 : -(shift+dim-1));CHKERRQ(ierr); 1433 break; 1434 } 1435 } 1436 ierr = DMPlexRestoreTransitiveClosure(dm, point, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 1437 } 1438 } 1439 } 1440 } 1441 ierr = ISRestoreIndices(dimIS, &points);CHKERRQ(ierr); 1442 ierr = ISDestroy(&dimIS);CHKERRQ(ierr); 1443 if (subdm) { 1444 if (subpointIS) {ierr = ISRestoreIndices(subpointIS, &subpoints);CHKERRQ(ierr);} 1445 ierr = ISDestroy(&subpointIS);CHKERRQ(ierr); 1446 } 1447 /* Mark boundary points as unsplit */ 1448 if (blabel) { 1449 ierr = DMLabelGetStratumIS(blabel, 1, &dimIS);CHKERRQ(ierr); 1450 ierr = ISGetLocalSize(dimIS, &numPoints);CHKERRQ(ierr); 1451 ierr = ISGetIndices(dimIS, &points);CHKERRQ(ierr); 1452 for (p = 0; p < numPoints; ++p) { 1453 const PetscInt point = points[p]; 1454 PetscInt val, bval; 1455 1456 ierr = DMLabelGetValue(blabel, point, &bval);CHKERRQ(ierr); 1457 if (bval >= 0) { 1458 const PetscInt *cone, *support; 1459 PetscInt coneSize, supportSize, s, valA, valB, valE; 1460 1461 /* Mark as unsplit */ 1462 ierr = DMLabelGetValue(label, point, &val);CHKERRQ(ierr); 1463 if ((val < 0) || (val > dim)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Point %d has label value %d, should be part of the fault", point, val); 1464 ierr = DMLabelClearValue(label, point, val);CHKERRQ(ierr); 1465 ierr = DMLabelSetValue(label, point, shift2+val);CHKERRQ(ierr); 1466 /* Check for cross-edge */ 1467 if (val != 0) continue; 1468 ierr = DMPlexGetSupport(dm, point, &support);CHKERRQ(ierr); 1469 ierr = DMPlexGetSupportSize(dm, point, &supportSize);CHKERRQ(ierr); 1470 for (s = 0; s < supportSize; ++s) { 1471 ierr = DMPlexGetCone(dm, support[s], &cone);CHKERRQ(ierr); 1472 ierr = DMPlexGetConeSize(dm, support[s], &coneSize);CHKERRQ(ierr); 1473 if (coneSize != 2) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Edge %D has %D vertices != 2", support[s], coneSize); 1474 ierr = DMLabelGetValue(blabel, cone[0], &valA);CHKERRQ(ierr); 1475 ierr = DMLabelGetValue(blabel, cone[1], &valB);CHKERRQ(ierr); 1476 ierr = DMLabelGetValue(blabel, support[s], &valE);CHKERRQ(ierr); 1477 if ((valE < 0) && (valA >= 0) && (valB >= 0)) {ierr = DMLabelSetValue(blabel, support[s], 2);CHKERRQ(ierr);} 1478 } 1479 } 1480 } 1481 ierr = ISRestoreIndices(dimIS, &points);CHKERRQ(ierr); 1482 ierr = ISDestroy(&dimIS);CHKERRQ(ierr); 1483 } 1484 /* Search for other cells/faces/edges connected to the fault by a vertex */ 1485 ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr); 1486 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 1487 ierr = DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);CHKERRQ(ierr); 1488 ierr = DMLabelGetStratumIS(label, 0, &dimIS);CHKERRQ(ierr); 1489 if (blabel) {ierr = DMLabelGetStratumIS(blabel, 2, &crossEdgeIS);CHKERRQ(ierr);} 1490 if (dimIS && crossEdgeIS) { 1491 IS vertIS = dimIS; 1492 1493 ierr = ISExpand(vertIS, crossEdgeIS, &dimIS);CHKERRQ(ierr); 1494 ierr = ISDestroy(&crossEdgeIS);CHKERRQ(ierr); 1495 ierr = ISDestroy(&vertIS);CHKERRQ(ierr); 1496 } 1497 if (!dimIS) PetscFunctionReturn(0); 1498 ierr = ISGetLocalSize(dimIS, &numPoints);CHKERRQ(ierr); 1499 ierr = ISGetIndices(dimIS, &points);CHKERRQ(ierr); 1500 for (p = 0; p < numPoints; ++p) { /* Loop over fault vertices */ 1501 PetscInt *star = NULL; 1502 PetscInt starSize, s; 1503 PetscInt again = 1; /* 0: Finished 1: Keep iterating after a change 2: No change */ 1504 1505 /* All points connected to the fault are inside a cell, so at the top level we will only check cells */ 1506 ierr = DMPlexGetTransitiveClosure(dm, points[p], PETSC_FALSE, &starSize, &star);CHKERRQ(ierr); 1507 while (again) { 1508 if (again > 1) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Could not classify all cells connected to the fault"); 1509 again = 0; 1510 for (s = 0; s < starSize*2; s += 2) { 1511 const PetscInt point = star[s]; 1512 const PetscInt *cone; 1513 PetscInt coneSize, c; 1514 1515 if ((point < cStart) || (point >= cEnd)) continue; 1516 ierr = DMLabelGetValue(label, point, &val);CHKERRQ(ierr); 1517 if (val != -1) continue; 1518 again = again == 1 ? 1 : 2; 1519 ierr = DMPlexGetConeSize(dm, point, &coneSize);CHKERRQ(ierr); 1520 ierr = DMPlexGetCone(dm, point, &cone);CHKERRQ(ierr); 1521 for (c = 0; c < coneSize; ++c) { 1522 ierr = DMLabelGetValue(label, cone[c], &val);CHKERRQ(ierr); 1523 if (val != -1) { 1524 const PetscInt *ccone; 1525 PetscInt cconeSize, cc, side; 1526 1527 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); 1528 if (val > 0) side = 1; 1529 else side = -1; 1530 ierr = DMLabelSetValue(label, point, side*(shift+dim));CHKERRQ(ierr); 1531 /* Mark cell faces which touch the fault */ 1532 ierr = DMPlexGetConeSize(dm, point, &cconeSize);CHKERRQ(ierr); 1533 ierr = DMPlexGetCone(dm, point, &ccone);CHKERRQ(ierr); 1534 for (cc = 0; cc < cconeSize; ++cc) { 1535 PetscInt *closure = NULL; 1536 PetscInt closureSize, cl; 1537 1538 ierr = DMLabelGetValue(label, ccone[cc], &val);CHKERRQ(ierr); 1539 if (val != -1) continue; 1540 ierr = DMPlexGetTransitiveClosure(dm, ccone[cc], PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 1541 for (cl = 0; cl < closureSize*2; cl += 2) { 1542 const PetscInt clp = closure[cl]; 1543 1544 ierr = DMLabelGetValue(label, clp, &val);CHKERRQ(ierr); 1545 if (val == -1) continue; 1546 ierr = DMLabelSetValue(label, ccone[cc], side*(shift+dim-1));CHKERRQ(ierr); 1547 break; 1548 } 1549 ierr = DMPlexRestoreTransitiveClosure(dm, ccone[cc], PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 1550 } 1551 again = 1; 1552 break; 1553 } 1554 } 1555 } 1556 } 1557 /* Classify the rest by cell membership */ 1558 for (s = 0; s < starSize*2; s += 2) { 1559 const PetscInt point = star[s]; 1560 1561 ierr = DMLabelGetValue(label, point, &val);CHKERRQ(ierr); 1562 if (val == -1) { 1563 PetscInt *sstar = NULL; 1564 PetscInt sstarSize, ss; 1565 PetscBool marked = PETSC_FALSE; 1566 1567 ierr = DMPlexGetTransitiveClosure(dm, point, PETSC_FALSE, &sstarSize, &sstar);CHKERRQ(ierr); 1568 for (ss = 0; ss < sstarSize*2; ss += 2) { 1569 const PetscInt spoint = sstar[ss]; 1570 1571 if ((spoint < cStart) || (spoint >= cEnd)) continue; 1572 ierr = DMLabelGetValue(label, spoint, &val);CHKERRQ(ierr); 1573 if (val == -1) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Cell %d in star of %d does not have a valid label", spoint, point); 1574 ierr = DMLabelGetValue(depthLabel, point, &dep);CHKERRQ(ierr); 1575 if (val > 0) { 1576 ierr = DMLabelSetValue(label, point, shift+dep);CHKERRQ(ierr); 1577 } else { 1578 ierr = DMLabelSetValue(label, point, -(shift+dep));CHKERRQ(ierr); 1579 } 1580 marked = PETSC_TRUE; 1581 break; 1582 } 1583 ierr = DMPlexRestoreTransitiveClosure(dm, point, PETSC_FALSE, &sstarSize, &sstar);CHKERRQ(ierr); 1584 if (!marked) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d could not be classified", point); 1585 } 1586 } 1587 ierr = DMPlexRestoreTransitiveClosure(dm, points[p], PETSC_FALSE, &starSize, &star);CHKERRQ(ierr); 1588 } 1589 ierr = ISRestoreIndices(dimIS, &points);CHKERRQ(ierr); 1590 ierr = ISDestroy(&dimIS);CHKERRQ(ierr); 1591 /* If any faces touching the fault divide cells on either side, split them */ 1592 ierr = DMLabelGetStratumIS(label, shift+dim-1, &facePosIS);CHKERRQ(ierr); 1593 ierr = DMLabelGetStratumIS(label, -(shift+dim-1), &faceNegIS);CHKERRQ(ierr); 1594 ierr = ISExpand(facePosIS, faceNegIS, &dimIS);CHKERRQ(ierr); 1595 ierr = ISDestroy(&facePosIS);CHKERRQ(ierr); 1596 ierr = ISDestroy(&faceNegIS);CHKERRQ(ierr); 1597 ierr = ISGetLocalSize(dimIS, &numPoints);CHKERRQ(ierr); 1598 ierr = ISGetIndices(dimIS, &points);CHKERRQ(ierr); 1599 for (p = 0; p < numPoints; ++p) { 1600 const PetscInt point = points[p]; 1601 const PetscInt *support; 1602 PetscInt supportSize, valA, valB; 1603 1604 ierr = DMPlexGetSupportSize(dm, point, &supportSize);CHKERRQ(ierr); 1605 if (supportSize != 2) continue; 1606 ierr = DMPlexGetSupport(dm, point, &support);CHKERRQ(ierr); 1607 ierr = DMLabelGetValue(label, support[0], &valA);CHKERRQ(ierr); 1608 ierr = DMLabelGetValue(label, support[1], &valB);CHKERRQ(ierr); 1609 if ((valA == -1) || (valB == -1)) continue; 1610 if (valA*valB > 0) continue; 1611 /* Split the face */ 1612 ierr = DMLabelGetValue(label, point, &valA);CHKERRQ(ierr); 1613 ierr = DMLabelClearValue(label, point, valA);CHKERRQ(ierr); 1614 ierr = DMLabelSetValue(label, point, dim-1);CHKERRQ(ierr); 1615 /* Label its closure: 1616 unmarked: label as unsplit 1617 incident: relabel as split 1618 split: do nothing 1619 */ 1620 { 1621 PetscInt *closure = NULL; 1622 PetscInt closureSize, cl; 1623 1624 ierr = DMPlexGetTransitiveClosure(dm, point, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 1625 for (cl = 0; cl < closureSize*2; cl += 2) { 1626 ierr = DMLabelGetValue(label, closure[cl], &valA);CHKERRQ(ierr); 1627 if (valA == -1) { /* Mark as unsplit */ 1628 ierr = DMLabelGetValue(depthLabel, closure[cl], &dep);CHKERRQ(ierr); 1629 ierr = DMLabelSetValue(label, closure[cl], shift2+dep);CHKERRQ(ierr); 1630 } else if (((valA >= shift) && (valA < shift2)) || ((valA <= -shift) && (valA > -shift2))) { 1631 ierr = DMLabelGetValue(depthLabel, closure[cl], &dep);CHKERRQ(ierr); 1632 ierr = DMLabelClearValue(label, closure[cl], valA);CHKERRQ(ierr); 1633 ierr = DMLabelSetValue(label, closure[cl], dep);CHKERRQ(ierr); 1634 } 1635 } 1636 ierr = DMPlexRestoreTransitiveClosure(dm, point, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 1637 } 1638 } 1639 ierr = ISRestoreIndices(dimIS, &points);CHKERRQ(ierr); 1640 ierr = ISDestroy(&dimIS);CHKERRQ(ierr); 1641 PetscFunctionReturn(0); 1642 } 1643 1644 #undef __FUNCT__ 1645 #define __FUNCT__ "DMPlexCreateHybridMesh" 1646 /*@C 1647 DMPlexCreateHybridMesh - Create a mesh with hybrid cells along an internal interface 1648 1649 Collective on dm 1650 1651 Input Parameters: 1652 + dm - The original DM 1653 - labelName - The label specifying the interface vertices 1654 1655 Output Parameters: 1656 + hybridLabel - The label fully marking the interface 1657 - dmHybrid - The new DM 1658 1659 Level: developer 1660 1661 .seealso: DMPlexConstructCohesiveCells(), DMPlexLabelCohesiveComplete(), DMCreate() 1662 @*/ 1663 PetscErrorCode DMPlexCreateHybridMesh(DM dm, DMLabel label, DMLabel *hybridLabel, DM *dmHybrid) 1664 { 1665 DM idm; 1666 DMLabel subpointMap, hlabel; 1667 PetscInt dim; 1668 PetscErrorCode ierr; 1669 1670 PetscFunctionBegin; 1671 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 1672 if (hybridLabel) PetscValidPointer(hybridLabel, 3); 1673 PetscValidPointer(dmHybrid, 4); 1674 ierr = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr); 1675 ierr = DMPlexCreateSubmesh(dm, label, 1, &idm);CHKERRQ(ierr); 1676 ierr = DMPlexOrient(idm);CHKERRQ(ierr); 1677 ierr = DMPlexGetSubpointMap(idm, &subpointMap);CHKERRQ(ierr); 1678 ierr = DMLabelDuplicate(subpointMap, &hlabel);CHKERRQ(ierr); 1679 ierr = DMLabelClearStratum(hlabel, dim);CHKERRQ(ierr); 1680 ierr = DMPlexLabelCohesiveComplete(dm, hlabel, NULL, PETSC_FALSE, idm);CHKERRQ(ierr); 1681 ierr = DMDestroy(&idm);CHKERRQ(ierr); 1682 ierr = DMPlexConstructCohesiveCells(dm, hlabel, dmHybrid);CHKERRQ(ierr); 1683 if (hybridLabel) *hybridLabel = hlabel; 1684 else {ierr = DMLabelDestroy(&hlabel);CHKERRQ(ierr);} 1685 PetscFunctionReturn(0); 1686 } 1687 1688 #undef __FUNCT__ 1689 #define __FUNCT__ "DMPlexMarkSubmesh_Uninterpolated" 1690 /* Here we need the explicit assumption that: 1691 1692 For any marked cell, the marked vertices constitute a single face 1693 */ 1694 static PetscErrorCode DMPlexMarkSubmesh_Uninterpolated(DM dm, DMLabel vertexLabel, PetscInt value, DMLabel subpointMap, PetscInt *numFaces, PetscInt *nFV, DM subdm) 1695 { 1696 IS subvertexIS = NULL; 1697 const PetscInt *subvertices; 1698 PetscInt *pStart, *pEnd, *pMax, pSize; 1699 PetscInt depth, dim, d, numSubVerticesInitial = 0, v; 1700 PetscErrorCode ierr; 1701 1702 PetscFunctionBegin; 1703 *numFaces = 0; 1704 *nFV = 0; 1705 ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr); 1706 ierr = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr); 1707 pSize = PetscMax(depth, dim) + 1; 1708 ierr = PetscMalloc3(pSize,&pStart,pSize,&pEnd,pSize,&pMax);CHKERRQ(ierr); 1709 ierr = DMPlexGetHybridBounds(dm, depth >= 0 ? &pMax[depth] : NULL, depth>1 ? &pMax[depth-1] : NULL, depth>2 ? &pMax[1] : NULL, &pMax[0]);CHKERRQ(ierr); 1710 for (d = 0; d <= depth; ++d) { 1711 ierr = DMPlexGetDepthStratum(dm, d, &pStart[d], &pEnd[d]);CHKERRQ(ierr); 1712 if (pMax[d] >= 0) pEnd[d] = PetscMin(pEnd[d], pMax[d]); 1713 } 1714 /* Loop over initial vertices and mark all faces in the collective star() */ 1715 if (vertexLabel) {ierr = DMLabelGetStratumIS(vertexLabel, value, &subvertexIS);CHKERRQ(ierr);} 1716 if (subvertexIS) { 1717 ierr = ISGetSize(subvertexIS, &numSubVerticesInitial);CHKERRQ(ierr); 1718 ierr = ISGetIndices(subvertexIS, &subvertices);CHKERRQ(ierr); 1719 } 1720 for (v = 0; v < numSubVerticesInitial; ++v) { 1721 const PetscInt vertex = subvertices[v]; 1722 PetscInt *star = NULL; 1723 PetscInt starSize, s, numCells = 0, c; 1724 1725 ierr = DMPlexGetTransitiveClosure(dm, vertex, PETSC_FALSE, &starSize, &star);CHKERRQ(ierr); 1726 for (s = 0; s < starSize*2; s += 2) { 1727 const PetscInt point = star[s]; 1728 if ((point >= pStart[depth]) && (point < pEnd[depth])) star[numCells++] = point; 1729 } 1730 for (c = 0; c < numCells; ++c) { 1731 const PetscInt cell = star[c]; 1732 PetscInt *closure = NULL; 1733 PetscInt closureSize, cl; 1734 PetscInt cellLoc, numCorners = 0, faceSize = 0; 1735 1736 ierr = DMLabelGetValue(subpointMap, cell, &cellLoc);CHKERRQ(ierr); 1737 if (cellLoc == 2) continue; 1738 if (cellLoc >= 0) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Cell %d has dimension %d in the surface label", cell, cellLoc); 1739 ierr = DMPlexGetTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 1740 for (cl = 0; cl < closureSize*2; cl += 2) { 1741 const PetscInt point = closure[cl]; 1742 PetscInt vertexLoc; 1743 1744 if ((point >= pStart[0]) && (point < pEnd[0])) { 1745 ++numCorners; 1746 ierr = DMLabelGetValue(vertexLabel, point, &vertexLoc);CHKERRQ(ierr); 1747 if (vertexLoc == value) closure[faceSize++] = point; 1748 } 1749 } 1750 if (!(*nFV)) {ierr = DMPlexGetNumFaceVertices(dm, dim, numCorners, nFV);CHKERRQ(ierr);} 1751 if (faceSize > *nFV) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Invalid submesh: Too many vertices %d of an element on the surface", faceSize); 1752 if (faceSize == *nFV) { 1753 const PetscInt *cells = NULL; 1754 PetscInt numCells, nc; 1755 1756 ++(*numFaces); 1757 for (cl = 0; cl < faceSize; ++cl) { 1758 ierr = DMLabelSetValue(subpointMap, closure[cl], 0);CHKERRQ(ierr); 1759 } 1760 ierr = DMPlexGetJoin(dm, faceSize, closure, &numCells, &cells);CHKERRQ(ierr); 1761 for (nc = 0; nc < numCells; ++nc) { 1762 ierr = DMLabelSetValue(subpointMap, cells[nc], 2);CHKERRQ(ierr); 1763 } 1764 ierr = DMPlexRestoreJoin(dm, faceSize, closure, &numCells, &cells);CHKERRQ(ierr); 1765 } 1766 ierr = DMPlexRestoreTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 1767 } 1768 ierr = DMPlexRestoreTransitiveClosure(dm, vertex, PETSC_FALSE, &starSize, &star);CHKERRQ(ierr); 1769 } 1770 if (subvertexIS) { 1771 ierr = ISRestoreIndices(subvertexIS, &subvertices);CHKERRQ(ierr); 1772 } 1773 ierr = ISDestroy(&subvertexIS);CHKERRQ(ierr); 1774 ierr = PetscFree3(pStart,pEnd,pMax);CHKERRQ(ierr); 1775 PetscFunctionReturn(0); 1776 } 1777 1778 #undef __FUNCT__ 1779 #define __FUNCT__ "DMPlexMarkSubmesh_Interpolated" 1780 static PetscErrorCode DMPlexMarkSubmesh_Interpolated(DM dm, DMLabel vertexLabel, PetscInt value, DMLabel subpointMap, DM subdm) 1781 { 1782 IS subvertexIS; 1783 const PetscInt *subvertices; 1784 PetscInt *pStart, *pEnd, *pMax; 1785 PetscInt dim, d, numSubVerticesInitial = 0, v; 1786 PetscErrorCode ierr; 1787 1788 PetscFunctionBegin; 1789 ierr = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr); 1790 ierr = PetscMalloc3(dim+1,&pStart,dim+1,&pEnd,dim+1,&pMax);CHKERRQ(ierr); 1791 ierr = DMPlexGetHybridBounds(dm, &pMax[dim], dim>1 ? &pMax[dim-1] : NULL, dim > 2 ? &pMax[1] : NULL, &pMax[0]);CHKERRQ(ierr); 1792 for (d = 0; d <= dim; ++d) { 1793 ierr = DMPlexGetDepthStratum(dm, d, &pStart[d], &pEnd[d]);CHKERRQ(ierr); 1794 if (pMax[d] >= 0) pEnd[d] = PetscMin(pEnd[d], pMax[d]); 1795 } 1796 /* Loop over initial vertices and mark all faces in the collective star() */ 1797 ierr = DMLabelGetStratumIS(vertexLabel, value, &subvertexIS);CHKERRQ(ierr); 1798 if (subvertexIS) { 1799 ierr = ISGetSize(subvertexIS, &numSubVerticesInitial);CHKERRQ(ierr); 1800 ierr = ISGetIndices(subvertexIS, &subvertices);CHKERRQ(ierr); 1801 } 1802 for (v = 0; v < numSubVerticesInitial; ++v) { 1803 const PetscInt vertex = subvertices[v]; 1804 PetscInt *star = NULL; 1805 PetscInt starSize, s, numFaces = 0, f; 1806 1807 ierr = DMPlexGetTransitiveClosure(dm, vertex, PETSC_FALSE, &starSize, &star);CHKERRQ(ierr); 1808 for (s = 0; s < starSize*2; s += 2) { 1809 const PetscInt point = star[s]; 1810 if ((point >= pStart[dim-1]) && (point < pEnd[dim-1])) star[numFaces++] = point; 1811 } 1812 for (f = 0; f < numFaces; ++f) { 1813 const PetscInt face = star[f]; 1814 PetscInt *closure = NULL; 1815 PetscInt closureSize, c; 1816 PetscInt faceLoc; 1817 1818 ierr = DMLabelGetValue(subpointMap, face, &faceLoc);CHKERRQ(ierr); 1819 if (faceLoc == dim-1) continue; 1820 if (faceLoc >= 0) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Face %d has dimension %d in the surface label", face, faceLoc); 1821 ierr = DMPlexGetTransitiveClosure(dm, face, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 1822 for (c = 0; c < closureSize*2; c += 2) { 1823 const PetscInt point = closure[c]; 1824 PetscInt vertexLoc; 1825 1826 if ((point >= pStart[0]) && (point < pEnd[0])) { 1827 ierr = DMLabelGetValue(vertexLabel, point, &vertexLoc);CHKERRQ(ierr); 1828 if (vertexLoc != value) break; 1829 } 1830 } 1831 if (c == closureSize*2) { 1832 const PetscInt *support; 1833 PetscInt supportSize, s; 1834 1835 for (c = 0; c < closureSize*2; c += 2) { 1836 const PetscInt point = closure[c]; 1837 1838 for (d = 0; d < dim; ++d) { 1839 if ((point >= pStart[d]) && (point < pEnd[d])) { 1840 ierr = DMLabelSetValue(subpointMap, point, d);CHKERRQ(ierr); 1841 break; 1842 } 1843 } 1844 } 1845 ierr = DMPlexGetSupportSize(dm, face, &supportSize);CHKERRQ(ierr); 1846 ierr = DMPlexGetSupport(dm, face, &support);CHKERRQ(ierr); 1847 for (s = 0; s < supportSize; ++s) { 1848 ierr = DMLabelSetValue(subpointMap, support[s], dim);CHKERRQ(ierr); 1849 } 1850 } 1851 ierr = DMPlexRestoreTransitiveClosure(dm, face, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 1852 } 1853 ierr = DMPlexRestoreTransitiveClosure(dm, vertex, PETSC_FALSE, &starSize, &star);CHKERRQ(ierr); 1854 } 1855 if (subvertexIS) { 1856 ierr = ISRestoreIndices(subvertexIS, &subvertices);CHKERRQ(ierr); 1857 } 1858 ierr = ISDestroy(&subvertexIS);CHKERRQ(ierr); 1859 ierr = PetscFree3(pStart,pEnd,pMax);CHKERRQ(ierr); 1860 PetscFunctionReturn(0); 1861 } 1862 1863 #undef __FUNCT__ 1864 #define __FUNCT__ "DMPlexMarkCohesiveSubmesh_Uninterpolated" 1865 static PetscErrorCode DMPlexMarkCohesiveSubmesh_Uninterpolated(DM dm, PetscBool hasLagrange, const char labelname[], PetscInt value, DMLabel subpointMap, PetscInt *numFaces, PetscInt *nFV, PetscInt *subCells[], DM subdm) 1866 { 1867 DMLabel label = NULL; 1868 const PetscInt *cone; 1869 PetscInt dim, cMax, cEnd, c, subc = 0, p, coneSize; 1870 PetscErrorCode ierr; 1871 1872 PetscFunctionBegin; 1873 *numFaces = 0; 1874 *nFV = 0; 1875 if (labelname) {ierr = DMPlexGetLabel(dm, labelname, &label);CHKERRQ(ierr);} 1876 *subCells = NULL; 1877 ierr = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr); 1878 ierr = DMPlexGetHeightStratum(dm, 0, NULL, &cEnd);CHKERRQ(ierr); 1879 ierr = DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);CHKERRQ(ierr); 1880 if (cMax < 0) PetscFunctionReturn(0); 1881 if (label) { 1882 for (c = cMax; c < cEnd; ++c) { 1883 PetscInt val; 1884 1885 ierr = DMLabelGetValue(label, c, &val);CHKERRQ(ierr); 1886 if (val == value) { 1887 ++(*numFaces); 1888 ierr = DMPlexGetConeSize(dm, c, &coneSize);CHKERRQ(ierr); 1889 } 1890 } 1891 } else { 1892 *numFaces = cEnd - cMax; 1893 ierr = DMPlexGetConeSize(dm, cMax, &coneSize);CHKERRQ(ierr); 1894 } 1895 *nFV = hasLagrange ? coneSize/3 : coneSize/2; 1896 ierr = PetscMalloc1(*numFaces *2, subCells);CHKERRQ(ierr); 1897 for (c = cMax; c < cEnd; ++c) { 1898 const PetscInt *cells; 1899 PetscInt numCells; 1900 1901 if (label) { 1902 PetscInt val; 1903 1904 ierr = DMLabelGetValue(label, c, &val);CHKERRQ(ierr); 1905 if (val != value) continue; 1906 } 1907 ierr = DMPlexGetCone(dm, c, &cone);CHKERRQ(ierr); 1908 for (p = 0; p < *nFV; ++p) { 1909 ierr = DMLabelSetValue(subpointMap, cone[p], 0);CHKERRQ(ierr); 1910 } 1911 /* Negative face */ 1912 ierr = DMPlexGetJoin(dm, *nFV, cone, &numCells, &cells);CHKERRQ(ierr); 1913 /* Not true in parallel 1914 if (numCells != 2) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cohesive cells should separate two cells"); */ 1915 for (p = 0; p < numCells; ++p) { 1916 ierr = DMLabelSetValue(subpointMap, cells[p], 2);CHKERRQ(ierr); 1917 (*subCells)[subc++] = cells[p]; 1918 } 1919 ierr = DMPlexRestoreJoin(dm, *nFV, cone, &numCells, &cells);CHKERRQ(ierr); 1920 /* Positive face is not included */ 1921 } 1922 PetscFunctionReturn(0); 1923 } 1924 1925 #undef __FUNCT__ 1926 #define __FUNCT__ "DMPlexMarkCohesiveSubmesh_Interpolated" 1927 static PetscErrorCode DMPlexMarkCohesiveSubmesh_Interpolated(DM dm, const char labelname[], PetscInt value, DMLabel subpointMap, DM subdm) 1928 { 1929 DMLabel label = NULL; 1930 PetscInt *pStart, *pEnd; 1931 PetscInt dim, cMax, cEnd, c, d; 1932 PetscErrorCode ierr; 1933 1934 PetscFunctionBegin; 1935 if (labelname) {ierr = DMPlexGetLabel(dm, labelname, &label);CHKERRQ(ierr);} 1936 ierr = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr); 1937 ierr = DMPlexGetHeightStratum(dm, 0, NULL, &cEnd);CHKERRQ(ierr); 1938 ierr = DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);CHKERRQ(ierr); 1939 if (cMax < 0) PetscFunctionReturn(0); 1940 ierr = PetscMalloc2(dim+1,&pStart,dim+1,&pEnd);CHKERRQ(ierr); 1941 for (d = 0; d <= dim; ++d) {ierr = DMPlexGetDepthStratum(dm, d, &pStart[d], &pEnd[d]);CHKERRQ(ierr);} 1942 for (c = cMax; c < cEnd; ++c) { 1943 const PetscInt *cone; 1944 PetscInt *closure = NULL; 1945 PetscInt fconeSize, coneSize, closureSize, cl, val; 1946 1947 if (label) { 1948 ierr = DMLabelGetValue(label, c, &val);CHKERRQ(ierr); 1949 if (val != value) continue; 1950 } 1951 ierr = DMPlexGetConeSize(dm, c, &coneSize);CHKERRQ(ierr); 1952 ierr = DMPlexGetCone(dm, c, &cone);CHKERRQ(ierr); 1953 ierr = DMPlexGetConeSize(dm, cone[0], &fconeSize);CHKERRQ(ierr); 1954 if (coneSize != (fconeSize ? fconeSize : 1) + 2) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cohesive cells should separate two cells"); 1955 /* Negative face */ 1956 ierr = DMPlexGetTransitiveClosure(dm, cone[0], PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 1957 for (cl = 0; cl < closureSize*2; cl += 2) { 1958 const PetscInt point = closure[cl]; 1959 1960 for (d = 0; d <= dim; ++d) { 1961 if ((point >= pStart[d]) && (point < pEnd[d])) { 1962 ierr = DMLabelSetValue(subpointMap, point, d);CHKERRQ(ierr); 1963 break; 1964 } 1965 } 1966 } 1967 ierr = DMPlexRestoreTransitiveClosure(dm, cone[0], PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 1968 /* Cells -- positive face is not included */ 1969 for (cl = 0; cl < 1; ++cl) { 1970 const PetscInt *support; 1971 PetscInt supportSize, s; 1972 1973 ierr = DMPlexGetSupportSize(dm, cone[cl], &supportSize);CHKERRQ(ierr); 1974 /* if (supportSize != 2) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cohesive faces should separate two cells"); */ 1975 ierr = DMPlexGetSupport(dm, cone[cl], &support);CHKERRQ(ierr); 1976 for (s = 0; s < supportSize; ++s) { 1977 ierr = DMLabelSetValue(subpointMap, support[s], dim);CHKERRQ(ierr); 1978 } 1979 } 1980 } 1981 ierr = PetscFree2(pStart, pEnd);CHKERRQ(ierr); 1982 PetscFunctionReturn(0); 1983 } 1984 1985 #undef __FUNCT__ 1986 #define __FUNCT__ "DMPlexGetFaceOrientation" 1987 PetscErrorCode DMPlexGetFaceOrientation(DM dm, PetscInt cell, PetscInt numCorners, PetscInt indices[], PetscInt oppositeVertex, PetscInt origVertices[], PetscInt faceVertices[], PetscBool *posOriented) 1988 { 1989 MPI_Comm comm; 1990 PetscBool posOrient = PETSC_FALSE; 1991 const PetscInt debug = 0; 1992 PetscInt cellDim, faceSize, f; 1993 PetscErrorCode ierr; 1994 1995 PetscFunctionBegin; 1996 ierr = PetscObjectGetComm((PetscObject)dm,&comm);CHKERRQ(ierr); 1997 ierr = DMPlexGetDimension(dm, &cellDim);CHKERRQ(ierr); 1998 if (debug) {PetscPrintf(comm, "cellDim: %d numCorners: %d\n", cellDim, numCorners);CHKERRQ(ierr);} 1999 2000 if (cellDim == 1 && numCorners == 2) { 2001 /* Triangle */ 2002 faceSize = numCorners-1; 2003 posOrient = !(oppositeVertex%2) ? PETSC_TRUE : PETSC_FALSE; 2004 } else if (cellDim == 2 && numCorners == 3) { 2005 /* Triangle */ 2006 faceSize = numCorners-1; 2007 posOrient = !(oppositeVertex%2) ? PETSC_TRUE : PETSC_FALSE; 2008 } else if (cellDim == 3 && numCorners == 4) { 2009 /* Tetrahedron */ 2010 faceSize = numCorners-1; 2011 posOrient = (oppositeVertex%2) ? PETSC_TRUE : PETSC_FALSE; 2012 } else if (cellDim == 1 && numCorners == 3) { 2013 /* Quadratic line */ 2014 faceSize = 1; 2015 posOrient = PETSC_TRUE; 2016 } else if (cellDim == 2 && numCorners == 4) { 2017 /* Quads */ 2018 faceSize = 2; 2019 if ((indices[1] > indices[0]) && (indices[1] - indices[0] == 1)) { 2020 posOrient = PETSC_TRUE; 2021 } else if ((indices[0] == 3) && (indices[1] == 0)) { 2022 posOrient = PETSC_TRUE; 2023 } else { 2024 if (((indices[0] > indices[1]) && (indices[0] - indices[1] == 1)) || ((indices[0] == 0) && (indices[1] == 3))) { 2025 posOrient = PETSC_FALSE; 2026 } else SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid quad crossedge"); 2027 } 2028 } else if (cellDim == 2 && numCorners == 6) { 2029 /* Quadratic triangle (I hate this) */ 2030 /* Edges are determined by the first 2 vertices (corners of edges) */ 2031 const PetscInt faceSizeTri = 3; 2032 PetscInt sortedIndices[3], i, iFace; 2033 PetscBool found = PETSC_FALSE; 2034 PetscInt faceVerticesTriSorted[9] = { 2035 0, 3, 4, /* bottom */ 2036 1, 4, 5, /* right */ 2037 2, 3, 5, /* left */ 2038 }; 2039 PetscInt faceVerticesTri[9] = { 2040 0, 3, 4, /* bottom */ 2041 1, 4, 5, /* right */ 2042 2, 5, 3, /* left */ 2043 }; 2044 2045 faceSize = faceSizeTri; 2046 for (i = 0; i < faceSizeTri; ++i) sortedIndices[i] = indices[i]; 2047 ierr = PetscSortInt(faceSizeTri, sortedIndices);CHKERRQ(ierr); 2048 for (iFace = 0; iFace < 3; ++iFace) { 2049 const PetscInt ii = iFace*faceSizeTri; 2050 PetscInt fVertex, cVertex; 2051 2052 if ((sortedIndices[0] == faceVerticesTriSorted[ii+0]) && 2053 (sortedIndices[1] == faceVerticesTriSorted[ii+1])) { 2054 for (fVertex = 0; fVertex < faceSizeTri; ++fVertex) { 2055 for (cVertex = 0; cVertex < faceSizeTri; ++cVertex) { 2056 if (indices[cVertex] == faceVerticesTri[ii+fVertex]) { 2057 faceVertices[fVertex] = origVertices[cVertex]; 2058 break; 2059 } 2060 } 2061 } 2062 found = PETSC_TRUE; 2063 break; 2064 } 2065 } 2066 if (!found) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid tri crossface"); 2067 if (posOriented) *posOriented = PETSC_TRUE; 2068 PetscFunctionReturn(0); 2069 } else if (cellDim == 2 && numCorners == 9) { 2070 /* Quadratic quad (I hate this) */ 2071 /* Edges are determined by the first 2 vertices (corners of edges) */ 2072 const PetscInt faceSizeQuad = 3; 2073 PetscInt sortedIndices[3], i, iFace; 2074 PetscBool found = PETSC_FALSE; 2075 PetscInt faceVerticesQuadSorted[12] = { 2076 0, 1, 4, /* bottom */ 2077 1, 2, 5, /* right */ 2078 2, 3, 6, /* top */ 2079 0, 3, 7, /* left */ 2080 }; 2081 PetscInt faceVerticesQuad[12] = { 2082 0, 1, 4, /* bottom */ 2083 1, 2, 5, /* right */ 2084 2, 3, 6, /* top */ 2085 3, 0, 7, /* left */ 2086 }; 2087 2088 faceSize = faceSizeQuad; 2089 for (i = 0; i < faceSizeQuad; ++i) sortedIndices[i] = indices[i]; 2090 ierr = PetscSortInt(faceSizeQuad, sortedIndices);CHKERRQ(ierr); 2091 for (iFace = 0; iFace < 4; ++iFace) { 2092 const PetscInt ii = iFace*faceSizeQuad; 2093 PetscInt fVertex, cVertex; 2094 2095 if ((sortedIndices[0] == faceVerticesQuadSorted[ii+0]) && 2096 (sortedIndices[1] == faceVerticesQuadSorted[ii+1])) { 2097 for (fVertex = 0; fVertex < faceSizeQuad; ++fVertex) { 2098 for (cVertex = 0; cVertex < faceSizeQuad; ++cVertex) { 2099 if (indices[cVertex] == faceVerticesQuad[ii+fVertex]) { 2100 faceVertices[fVertex] = origVertices[cVertex]; 2101 break; 2102 } 2103 } 2104 } 2105 found = PETSC_TRUE; 2106 break; 2107 } 2108 } 2109 if (!found) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid quad crossface"); 2110 if (posOriented) *posOriented = PETSC_TRUE; 2111 PetscFunctionReturn(0); 2112 } else if (cellDim == 3 && numCorners == 8) { 2113 /* Hexes 2114 A hex is two oriented quads with the normal of the first 2115 pointing up at the second. 2116 2117 7---6 2118 /| /| 2119 4---5 | 2120 | 1-|-2 2121 |/ |/ 2122 0---3 2123 2124 Faces are determined by the first 4 vertices (corners of faces) */ 2125 const PetscInt faceSizeHex = 4; 2126 PetscInt sortedIndices[4], i, iFace; 2127 PetscBool found = PETSC_FALSE; 2128 PetscInt faceVerticesHexSorted[24] = { 2129 0, 1, 2, 3, /* bottom */ 2130 4, 5, 6, 7, /* top */ 2131 0, 3, 4, 5, /* front */ 2132 2, 3, 5, 6, /* right */ 2133 1, 2, 6, 7, /* back */ 2134 0, 1, 4, 7, /* left */ 2135 }; 2136 PetscInt faceVerticesHex[24] = { 2137 1, 2, 3, 0, /* bottom */ 2138 4, 5, 6, 7, /* top */ 2139 0, 3, 5, 4, /* front */ 2140 3, 2, 6, 5, /* right */ 2141 2, 1, 7, 6, /* back */ 2142 1, 0, 4, 7, /* left */ 2143 }; 2144 2145 faceSize = faceSizeHex; 2146 for (i = 0; i < faceSizeHex; ++i) sortedIndices[i] = indices[i]; 2147 ierr = PetscSortInt(faceSizeHex, sortedIndices);CHKERRQ(ierr); 2148 for (iFace = 0; iFace < 6; ++iFace) { 2149 const PetscInt ii = iFace*faceSizeHex; 2150 PetscInt fVertex, cVertex; 2151 2152 if ((sortedIndices[0] == faceVerticesHexSorted[ii+0]) && 2153 (sortedIndices[1] == faceVerticesHexSorted[ii+1]) && 2154 (sortedIndices[2] == faceVerticesHexSorted[ii+2]) && 2155 (sortedIndices[3] == faceVerticesHexSorted[ii+3])) { 2156 for (fVertex = 0; fVertex < faceSizeHex; ++fVertex) { 2157 for (cVertex = 0; cVertex < faceSizeHex; ++cVertex) { 2158 if (indices[cVertex] == faceVerticesHex[ii+fVertex]) { 2159 faceVertices[fVertex] = origVertices[cVertex]; 2160 break; 2161 } 2162 } 2163 } 2164 found = PETSC_TRUE; 2165 break; 2166 } 2167 } 2168 if (!found) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid hex crossface"); 2169 if (posOriented) *posOriented = PETSC_TRUE; 2170 PetscFunctionReturn(0); 2171 } else if (cellDim == 3 && numCorners == 10) { 2172 /* Quadratic tet */ 2173 /* Faces are determined by the first 3 vertices (corners of faces) */ 2174 const PetscInt faceSizeTet = 6; 2175 PetscInt sortedIndices[6], i, iFace; 2176 PetscBool found = PETSC_FALSE; 2177 PetscInt faceVerticesTetSorted[24] = { 2178 0, 1, 2, 6, 7, 8, /* bottom */ 2179 0, 3, 4, 6, 7, 9, /* front */ 2180 1, 4, 5, 7, 8, 9, /* right */ 2181 2, 3, 5, 6, 8, 9, /* left */ 2182 }; 2183 PetscInt faceVerticesTet[24] = { 2184 0, 1, 2, 6, 7, 8, /* bottom */ 2185 0, 4, 3, 6, 7, 9, /* front */ 2186 1, 5, 4, 7, 8, 9, /* right */ 2187 2, 3, 5, 8, 6, 9, /* left */ 2188 }; 2189 2190 faceSize = faceSizeTet; 2191 for (i = 0; i < faceSizeTet; ++i) sortedIndices[i] = indices[i]; 2192 ierr = PetscSortInt(faceSizeTet, sortedIndices);CHKERRQ(ierr); 2193 for (iFace=0; iFace < 4; ++iFace) { 2194 const PetscInt ii = iFace*faceSizeTet; 2195 PetscInt fVertex, cVertex; 2196 2197 if ((sortedIndices[0] == faceVerticesTetSorted[ii+0]) && 2198 (sortedIndices[1] == faceVerticesTetSorted[ii+1]) && 2199 (sortedIndices[2] == faceVerticesTetSorted[ii+2]) && 2200 (sortedIndices[3] == faceVerticesTetSorted[ii+3])) { 2201 for (fVertex = 0; fVertex < faceSizeTet; ++fVertex) { 2202 for (cVertex = 0; cVertex < faceSizeTet; ++cVertex) { 2203 if (indices[cVertex] == faceVerticesTet[ii+fVertex]) { 2204 faceVertices[fVertex] = origVertices[cVertex]; 2205 break; 2206 } 2207 } 2208 } 2209 found = PETSC_TRUE; 2210 break; 2211 } 2212 } 2213 if (!found) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid tet crossface"); 2214 if (posOriented) *posOriented = PETSC_TRUE; 2215 PetscFunctionReturn(0); 2216 } else if (cellDim == 3 && numCorners == 27) { 2217 /* Quadratic hexes (I hate this) 2218 A hex is two oriented quads with the normal of the first 2219 pointing up at the second. 2220 2221 7---6 2222 /| /| 2223 4---5 | 2224 | 3-|-2 2225 |/ |/ 2226 0---1 2227 2228 Faces are determined by the first 4 vertices (corners of faces) */ 2229 const PetscInt faceSizeQuadHex = 9; 2230 PetscInt sortedIndices[9], i, iFace; 2231 PetscBool found = PETSC_FALSE; 2232 PetscInt faceVerticesQuadHexSorted[54] = { 2233 0, 1, 2, 3, 8, 9, 10, 11, 24, /* bottom */ 2234 4, 5, 6, 7, 12, 13, 14, 15, 25, /* top */ 2235 0, 1, 4, 5, 8, 12, 16, 17, 22, /* front */ 2236 1, 2, 5, 6, 9, 13, 17, 18, 21, /* right */ 2237 2, 3, 6, 7, 10, 14, 18, 19, 23, /* back */ 2238 0, 3, 4, 7, 11, 15, 16, 19, 20, /* left */ 2239 }; 2240 PetscInt faceVerticesQuadHex[54] = { 2241 3, 2, 1, 0, 10, 9, 8, 11, 24, /* bottom */ 2242 4, 5, 6, 7, 12, 13, 14, 15, 25, /* top */ 2243 0, 1, 5, 4, 8, 17, 12, 16, 22, /* front */ 2244 1, 2, 6, 5, 9, 18, 13, 17, 21, /* right */ 2245 2, 3, 7, 6, 10, 19, 14, 18, 23, /* back */ 2246 3, 0, 4, 7, 11, 16, 15, 19, 20 /* left */ 2247 }; 2248 2249 faceSize = faceSizeQuadHex; 2250 for (i = 0; i < faceSizeQuadHex; ++i) sortedIndices[i] = indices[i]; 2251 ierr = PetscSortInt(faceSizeQuadHex, sortedIndices);CHKERRQ(ierr); 2252 for (iFace = 0; iFace < 6; ++iFace) { 2253 const PetscInt ii = iFace*faceSizeQuadHex; 2254 PetscInt fVertex, cVertex; 2255 2256 if ((sortedIndices[0] == faceVerticesQuadHexSorted[ii+0]) && 2257 (sortedIndices[1] == faceVerticesQuadHexSorted[ii+1]) && 2258 (sortedIndices[2] == faceVerticesQuadHexSorted[ii+2]) && 2259 (sortedIndices[3] == faceVerticesQuadHexSorted[ii+3])) { 2260 for (fVertex = 0; fVertex < faceSizeQuadHex; ++fVertex) { 2261 for (cVertex = 0; cVertex < faceSizeQuadHex; ++cVertex) { 2262 if (indices[cVertex] == faceVerticesQuadHex[ii+fVertex]) { 2263 faceVertices[fVertex] = origVertices[cVertex]; 2264 break; 2265 } 2266 } 2267 } 2268 found = PETSC_TRUE; 2269 break; 2270 } 2271 } 2272 if (!found) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid hex crossface"); 2273 if (posOriented) *posOriented = PETSC_TRUE; 2274 PetscFunctionReturn(0); 2275 } else SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Unknown cell type for faceOrientation()."); 2276 if (!posOrient) { 2277 if (debug) {ierr = PetscPrintf(comm, " Reversing initial face orientation\n");CHKERRQ(ierr);} 2278 for (f = 0; f < faceSize; ++f) faceVertices[f] = origVertices[faceSize-1 - f]; 2279 } else { 2280 if (debug) {ierr = PetscPrintf(comm, " Keeping initial face orientation\n");CHKERRQ(ierr);} 2281 for (f = 0; f < faceSize; ++f) faceVertices[f] = origVertices[f]; 2282 } 2283 if (posOriented) *posOriented = posOrient; 2284 PetscFunctionReturn(0); 2285 } 2286 2287 #undef __FUNCT__ 2288 #define __FUNCT__ "DMPlexGetOrientedFace" 2289 /* 2290 Given a cell and a face, as a set of vertices, 2291 return the oriented face, as a set of vertices, in faceVertices 2292 The orientation is such that the face normal points out of the cell 2293 */ 2294 PetscErrorCode DMPlexGetOrientedFace(DM dm, PetscInt cell, PetscInt faceSize, const PetscInt face[], PetscInt numCorners, PetscInt indices[], PetscInt origVertices[], PetscInt faceVertices[], PetscBool *posOriented) 2295 { 2296 const PetscInt *cone = NULL; 2297 PetscInt coneSize, v, f, v2; 2298 PetscInt oppositeVertex = -1; 2299 PetscErrorCode ierr; 2300 2301 PetscFunctionBegin; 2302 ierr = DMPlexGetConeSize(dm, cell, &coneSize);CHKERRQ(ierr); 2303 ierr = DMPlexGetCone(dm, cell, &cone);CHKERRQ(ierr); 2304 for (v = 0, v2 = 0; v < coneSize; ++v) { 2305 PetscBool found = PETSC_FALSE; 2306 2307 for (f = 0; f < faceSize; ++f) { 2308 if (face[f] == cone[v]) { 2309 found = PETSC_TRUE; break; 2310 } 2311 } 2312 if (found) { 2313 indices[v2] = v; 2314 origVertices[v2] = cone[v]; 2315 ++v2; 2316 } else { 2317 oppositeVertex = v; 2318 } 2319 } 2320 ierr = DMPlexGetFaceOrientation(dm, cell, numCorners, indices, oppositeVertex, origVertices, faceVertices, posOriented);CHKERRQ(ierr); 2321 PetscFunctionReturn(0); 2322 } 2323 2324 #undef __FUNCT__ 2325 #define __FUNCT__ "DMPlexInsertFace_Internal" 2326 /* 2327 DMPlexInsertFace_Internal - Puts a face into the mesh 2328 2329 Not collective 2330 2331 Input Parameters: 2332 + dm - The DMPlex 2333 . numFaceVertex - The number of vertices in the face 2334 . faceVertices - The vertices in the face for dm 2335 . subfaceVertices - The vertices in the face for subdm 2336 . numCorners - The number of vertices in the cell 2337 . cell - A cell in dm containing the face 2338 . subcell - A cell in subdm containing the face 2339 . firstFace - First face in the mesh 2340 - newFacePoint - Next face in the mesh 2341 2342 Output Parameters: 2343 . newFacePoint - Contains next face point number on input, updated on output 2344 2345 Level: developer 2346 */ 2347 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) 2348 { 2349 MPI_Comm comm; 2350 DM_Plex *submesh = (DM_Plex*) subdm->data; 2351 const PetscInt *faces; 2352 PetscInt numFaces, coneSize; 2353 PetscErrorCode ierr; 2354 2355 PetscFunctionBegin; 2356 ierr = PetscObjectGetComm((PetscObject)dm,&comm);CHKERRQ(ierr); 2357 ierr = DMPlexGetConeSize(subdm, subcell, &coneSize);CHKERRQ(ierr); 2358 if (coneSize != 1) SETERRQ2(comm, PETSC_ERR_ARG_OUTOFRANGE, "Cone size of cell %d is %d != 1", cell, coneSize); 2359 #if 0 2360 /* Cannot use this because support() has not been constructed yet */ 2361 ierr = DMPlexGetJoin(subdm, numFaceVertices, subfaceVertices, &numFaces, &faces);CHKERRQ(ierr); 2362 #else 2363 { 2364 PetscInt f; 2365 2366 numFaces = 0; 2367 ierr = DMGetWorkArray(subdm, 1, PETSC_INT, (void **) &faces);CHKERRQ(ierr); 2368 for (f = firstFace; f < *newFacePoint; ++f) { 2369 PetscInt dof, off, d; 2370 2371 ierr = PetscSectionGetDof(submesh->coneSection, f, &dof);CHKERRQ(ierr); 2372 ierr = PetscSectionGetOffset(submesh->coneSection, f, &off);CHKERRQ(ierr); 2373 /* Yes, I know this is quadratic, but I expect the sizes to be <5 */ 2374 for (d = 0; d < dof; ++d) { 2375 const PetscInt p = submesh->cones[off+d]; 2376 PetscInt v; 2377 2378 for (v = 0; v < numFaceVertices; ++v) { 2379 if (subfaceVertices[v] == p) break; 2380 } 2381 if (v == numFaceVertices) break; 2382 } 2383 if (d == dof) { 2384 numFaces = 1; 2385 ((PetscInt*) faces)[0] = f; 2386 } 2387 } 2388 } 2389 #endif 2390 if (numFaces > 1) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Vertex set had %d faces, not one", numFaces); 2391 else if (numFaces == 1) { 2392 /* Add the other cell neighbor for this face */ 2393 ierr = DMPlexSetCone(subdm, subcell, faces);CHKERRQ(ierr); 2394 } else { 2395 PetscInt *indices, *origVertices, *orientedVertices, *orientedSubVertices, v, ov; 2396 PetscBool posOriented; 2397 2398 ierr = DMGetWorkArray(subdm, 4*numFaceVertices * sizeof(PetscInt), PETSC_INT, &orientedVertices);CHKERRQ(ierr); 2399 origVertices = &orientedVertices[numFaceVertices]; 2400 indices = &orientedVertices[numFaceVertices*2]; 2401 orientedSubVertices = &orientedVertices[numFaceVertices*3]; 2402 ierr = DMPlexGetOrientedFace(dm, cell, numFaceVertices, faceVertices, numCorners, indices, origVertices, orientedVertices, &posOriented);CHKERRQ(ierr); 2403 /* TODO: I know that routine should return a permutation, not the indices */ 2404 for (v = 0; v < numFaceVertices; ++v) { 2405 const PetscInt vertex = faceVertices[v], subvertex = subfaceVertices[v]; 2406 for (ov = 0; ov < numFaceVertices; ++ov) { 2407 if (orientedVertices[ov] == vertex) { 2408 orientedSubVertices[ov] = subvertex; 2409 break; 2410 } 2411 } 2412 if (ov == numFaceVertices) SETERRQ1(comm, PETSC_ERR_PLIB, "Could not find face vertex %d in orientated set", vertex); 2413 } 2414 ierr = DMPlexSetCone(subdm, *newFacePoint, orientedSubVertices);CHKERRQ(ierr); 2415 ierr = DMPlexSetCone(subdm, subcell, newFacePoint);CHKERRQ(ierr); 2416 ierr = DMRestoreWorkArray(subdm, 4*numFaceVertices * sizeof(PetscInt), PETSC_INT, &orientedVertices);CHKERRQ(ierr); 2417 ++(*newFacePoint); 2418 } 2419 #if 0 2420 ierr = DMPlexRestoreJoin(subdm, numFaceVertices, subfaceVertices, &numFaces, &faces);CHKERRQ(ierr); 2421 #else 2422 ierr = DMRestoreWorkArray(subdm, 1, PETSC_INT, (void **) &faces);CHKERRQ(ierr); 2423 #endif 2424 PetscFunctionReturn(0); 2425 } 2426 2427 #undef __FUNCT__ 2428 #define __FUNCT__ "DMPlexCreateSubmesh_Uninterpolated" 2429 static PetscErrorCode DMPlexCreateSubmesh_Uninterpolated(DM dm, DMLabel vertexLabel, PetscInt value, DM subdm) 2430 { 2431 MPI_Comm comm; 2432 DMLabel subpointMap; 2433 IS subvertexIS, subcellIS; 2434 const PetscInt *subVertices, *subCells; 2435 PetscInt numSubVertices, firstSubVertex, numSubCells; 2436 PetscInt *subface, maxConeSize, numSubFaces = 0, firstSubFace, newFacePoint, nFV = 0; 2437 PetscInt vStart, vEnd, c, f; 2438 PetscErrorCode ierr; 2439 2440 PetscFunctionBegin; 2441 ierr = PetscObjectGetComm((PetscObject)dm,&comm);CHKERRQ(ierr); 2442 /* Create subpointMap which marks the submesh */ 2443 ierr = DMLabelCreate("subpoint_map", &subpointMap);CHKERRQ(ierr); 2444 ierr = DMPlexSetSubpointMap(subdm, subpointMap);CHKERRQ(ierr); 2445 ierr = DMLabelDestroy(&subpointMap);CHKERRQ(ierr); 2446 if (vertexLabel) {ierr = DMPlexMarkSubmesh_Uninterpolated(dm, vertexLabel, value, subpointMap, &numSubFaces, &nFV, subdm);CHKERRQ(ierr);} 2447 /* Setup chart */ 2448 ierr = DMLabelGetStratumSize(subpointMap, 0, &numSubVertices);CHKERRQ(ierr); 2449 ierr = DMLabelGetStratumSize(subpointMap, 2, &numSubCells);CHKERRQ(ierr); 2450 ierr = DMPlexSetChart(subdm, 0, numSubCells+numSubFaces+numSubVertices);CHKERRQ(ierr); 2451 ierr = DMPlexSetVTKCellHeight(subdm, 1);CHKERRQ(ierr); 2452 /* Set cone sizes */ 2453 firstSubVertex = numSubCells; 2454 firstSubFace = numSubCells+numSubVertices; 2455 newFacePoint = firstSubFace; 2456 ierr = DMLabelGetStratumIS(subpointMap, 0, &subvertexIS);CHKERRQ(ierr); 2457 if (subvertexIS) {ierr = ISGetIndices(subvertexIS, &subVertices);CHKERRQ(ierr);} 2458 ierr = DMLabelGetStratumIS(subpointMap, 2, &subcellIS);CHKERRQ(ierr); 2459 if (subcellIS) {ierr = ISGetIndices(subcellIS, &subCells);CHKERRQ(ierr);} 2460 for (c = 0; c < numSubCells; ++c) { 2461 ierr = DMPlexSetConeSize(subdm, c, 1);CHKERRQ(ierr); 2462 } 2463 for (f = firstSubFace; f < firstSubFace+numSubFaces; ++f) { 2464 ierr = DMPlexSetConeSize(subdm, f, nFV);CHKERRQ(ierr); 2465 } 2466 ierr = DMSetUp(subdm);CHKERRQ(ierr); 2467 /* Create face cones */ 2468 ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr); 2469 ierr = DMPlexGetMaxSizes(dm, &maxConeSize, NULL);CHKERRQ(ierr); 2470 ierr = DMGetWorkArray(subdm, maxConeSize, PETSC_INT, (void**) &subface);CHKERRQ(ierr); 2471 for (c = 0; c < numSubCells; ++c) { 2472 const PetscInt cell = subCells[c]; 2473 const PetscInt subcell = c; 2474 PetscInt *closure = NULL; 2475 PetscInt closureSize, cl, numCorners = 0, faceSize = 0; 2476 2477 ierr = DMPlexGetTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 2478 for (cl = 0; cl < closureSize*2; cl += 2) { 2479 const PetscInt point = closure[cl]; 2480 PetscInt subVertex; 2481 2482 if ((point >= vStart) && (point < vEnd)) { 2483 ++numCorners; 2484 ierr = PetscFindInt(point, numSubVertices, subVertices, &subVertex);CHKERRQ(ierr); 2485 if (subVertex >= 0) { 2486 closure[faceSize] = point; 2487 subface[faceSize] = firstSubVertex+subVertex; 2488 ++faceSize; 2489 } 2490 } 2491 } 2492 if (faceSize > nFV) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Invalid submesh: Too many vertices %d of an element on the surface", faceSize); 2493 if (faceSize == nFV) { 2494 ierr = DMPlexInsertFace_Internal(dm, subdm, faceSize, closure, subface, numCorners, cell, subcell, firstSubFace, &newFacePoint);CHKERRQ(ierr); 2495 } 2496 ierr = DMPlexRestoreTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 2497 } 2498 ierr = DMRestoreWorkArray(subdm, maxConeSize, PETSC_INT, (void**) &subface);CHKERRQ(ierr); 2499 ierr = DMPlexSymmetrize(subdm);CHKERRQ(ierr); 2500 ierr = DMPlexStratify(subdm);CHKERRQ(ierr); 2501 /* Build coordinates */ 2502 { 2503 PetscSection coordSection, subCoordSection; 2504 Vec coordinates, subCoordinates; 2505 PetscScalar *coords, *subCoords; 2506 PetscInt numComp, coordSize, v; 2507 2508 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 2509 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr); 2510 ierr = DMGetCoordinateSection(subdm, &subCoordSection);CHKERRQ(ierr); 2511 ierr = PetscSectionSetNumFields(subCoordSection, 1);CHKERRQ(ierr); 2512 ierr = PetscSectionGetFieldComponents(coordSection, 0, &numComp);CHKERRQ(ierr); 2513 ierr = PetscSectionSetFieldComponents(subCoordSection, 0, numComp);CHKERRQ(ierr); 2514 ierr = PetscSectionSetChart(subCoordSection, firstSubVertex, firstSubVertex+numSubVertices);CHKERRQ(ierr); 2515 for (v = 0; v < numSubVertices; ++v) { 2516 const PetscInt vertex = subVertices[v]; 2517 const PetscInt subvertex = firstSubVertex+v; 2518 PetscInt dof; 2519 2520 ierr = PetscSectionGetDof(coordSection, vertex, &dof);CHKERRQ(ierr); 2521 ierr = PetscSectionSetDof(subCoordSection, subvertex, dof);CHKERRQ(ierr); 2522 ierr = PetscSectionSetFieldDof(subCoordSection, subvertex, 0, dof);CHKERRQ(ierr); 2523 } 2524 ierr = PetscSectionSetUp(subCoordSection);CHKERRQ(ierr); 2525 ierr = PetscSectionGetStorageSize(subCoordSection, &coordSize);CHKERRQ(ierr); 2526 ierr = VecCreate(comm, &subCoordinates);CHKERRQ(ierr); 2527 ierr = VecSetSizes(subCoordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 2528 ierr = VecSetType(subCoordinates,VECSTANDARD);CHKERRQ(ierr); 2529 if (coordSize) { 2530 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 2531 ierr = VecGetArray(subCoordinates, &subCoords);CHKERRQ(ierr); 2532 for (v = 0; v < numSubVertices; ++v) { 2533 const PetscInt vertex = subVertices[v]; 2534 const PetscInt subvertex = firstSubVertex+v; 2535 PetscInt dof, off, sdof, soff, d; 2536 2537 ierr = PetscSectionGetDof(coordSection, vertex, &dof);CHKERRQ(ierr); 2538 ierr = PetscSectionGetOffset(coordSection, vertex, &off);CHKERRQ(ierr); 2539 ierr = PetscSectionGetDof(subCoordSection, subvertex, &sdof);CHKERRQ(ierr); 2540 ierr = PetscSectionGetOffset(subCoordSection, subvertex, &soff);CHKERRQ(ierr); 2541 if (dof != sdof) SETERRQ4(comm, PETSC_ERR_PLIB, "Coordinate dimension %d on subvertex %d, vertex %d should be %d", sdof, subvertex, vertex, dof); 2542 for (d = 0; d < dof; ++d) subCoords[soff+d] = coords[off+d]; 2543 } 2544 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 2545 ierr = VecRestoreArray(subCoordinates, &subCoords);CHKERRQ(ierr); 2546 } 2547 ierr = DMSetCoordinatesLocal(subdm, subCoordinates);CHKERRQ(ierr); 2548 ierr = VecDestroy(&subCoordinates);CHKERRQ(ierr); 2549 } 2550 /* Cleanup */ 2551 if (subvertexIS) {ierr = ISRestoreIndices(subvertexIS, &subVertices);CHKERRQ(ierr);} 2552 ierr = ISDestroy(&subvertexIS);CHKERRQ(ierr); 2553 if (subcellIS) {ierr = ISRestoreIndices(subcellIS, &subCells);CHKERRQ(ierr);} 2554 ierr = ISDestroy(&subcellIS);CHKERRQ(ierr); 2555 PetscFunctionReturn(0); 2556 } 2557 2558 #undef __FUNCT__ 2559 #define __FUNCT__ "DMPlexCreateSubmesh_Interpolated" 2560 static PetscErrorCode DMPlexCreateSubmesh_Interpolated(DM dm, DMLabel vertexLabel, PetscInt value, DM subdm) 2561 { 2562 MPI_Comm comm; 2563 DMLabel subpointMap; 2564 IS *subpointIS; 2565 const PetscInt **subpoints; 2566 PetscInt *numSubPoints, *firstSubPoint, *coneNew, *coneONew; 2567 PetscInt totSubPoints = 0, maxConeSize, dim, p, d, v; 2568 PetscErrorCode ierr; 2569 2570 PetscFunctionBegin; 2571 ierr = PetscObjectGetComm((PetscObject)dm,&comm);CHKERRQ(ierr); 2572 /* Create subpointMap which marks the submesh */ 2573 ierr = DMLabelCreate("subpoint_map", &subpointMap);CHKERRQ(ierr); 2574 ierr = DMPlexSetSubpointMap(subdm, subpointMap);CHKERRQ(ierr); 2575 ierr = DMLabelDestroy(&subpointMap);CHKERRQ(ierr); 2576 if (vertexLabel) {ierr = DMPlexMarkSubmesh_Interpolated(dm, vertexLabel, value, subpointMap, subdm);CHKERRQ(ierr);} 2577 /* Setup chart */ 2578 ierr = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr); 2579 ierr = PetscMalloc4(dim+1,&numSubPoints,dim+1,&firstSubPoint,dim+1,&subpointIS,dim+1,&subpoints);CHKERRQ(ierr); 2580 for (d = 0; d <= dim; ++d) { 2581 ierr = DMLabelGetStratumSize(subpointMap, d, &numSubPoints[d]);CHKERRQ(ierr); 2582 totSubPoints += numSubPoints[d]; 2583 } 2584 ierr = DMPlexSetChart(subdm, 0, totSubPoints);CHKERRQ(ierr); 2585 ierr = DMPlexSetVTKCellHeight(subdm, 1);CHKERRQ(ierr); 2586 /* Set cone sizes */ 2587 firstSubPoint[dim] = 0; 2588 firstSubPoint[0] = firstSubPoint[dim] + numSubPoints[dim]; 2589 if (dim > 1) {firstSubPoint[dim-1] = firstSubPoint[0] + numSubPoints[0];} 2590 if (dim > 2) {firstSubPoint[dim-2] = firstSubPoint[dim-1] + numSubPoints[dim-1];} 2591 for (d = 0; d <= dim; ++d) { 2592 ierr = DMLabelGetStratumIS(subpointMap, d, &subpointIS[d]);CHKERRQ(ierr); 2593 if (subpointIS[d]) {ierr = ISGetIndices(subpointIS[d], &subpoints[d]);CHKERRQ(ierr);} 2594 } 2595 for (d = 0; d <= dim; ++d) { 2596 for (p = 0; p < numSubPoints[d]; ++p) { 2597 const PetscInt point = subpoints[d][p]; 2598 const PetscInt subpoint = firstSubPoint[d] + p; 2599 const PetscInt *cone; 2600 PetscInt coneSize, coneSizeNew, c, val; 2601 2602 ierr = DMPlexGetConeSize(dm, point, &coneSize);CHKERRQ(ierr); 2603 ierr = DMPlexSetConeSize(subdm, subpoint, coneSize);CHKERRQ(ierr); 2604 if (d == dim) { 2605 ierr = DMPlexGetCone(dm, point, &cone);CHKERRQ(ierr); 2606 for (c = 0, coneSizeNew = 0; c < coneSize; ++c) { 2607 ierr = DMLabelGetValue(subpointMap, cone[c], &val);CHKERRQ(ierr); 2608 if (val >= 0) coneSizeNew++; 2609 } 2610 ierr = DMPlexSetConeSize(subdm, subpoint, coneSizeNew);CHKERRQ(ierr); 2611 } 2612 } 2613 } 2614 ierr = DMSetUp(subdm);CHKERRQ(ierr); 2615 /* Set cones */ 2616 ierr = DMPlexGetMaxSizes(dm, &maxConeSize, NULL);CHKERRQ(ierr); 2617 ierr = PetscMalloc2(maxConeSize,&coneNew,maxConeSize,&coneONew);CHKERRQ(ierr); 2618 for (d = 0; d <= dim; ++d) { 2619 for (p = 0; p < numSubPoints[d]; ++p) { 2620 const PetscInt point = subpoints[d][p]; 2621 const PetscInt subpoint = firstSubPoint[d] + p; 2622 const PetscInt *cone, *ornt; 2623 PetscInt coneSize, subconeSize, coneSizeNew, c, subc; 2624 2625 ierr = DMPlexGetConeSize(dm, point, &coneSize);CHKERRQ(ierr); 2626 ierr = DMPlexGetConeSize(subdm, subpoint, &subconeSize);CHKERRQ(ierr); 2627 ierr = DMPlexGetCone(dm, point, &cone);CHKERRQ(ierr); 2628 ierr = DMPlexGetConeOrientation(dm, point, &ornt);CHKERRQ(ierr); 2629 for (c = 0, coneSizeNew = 0; c < coneSize; ++c) { 2630 ierr = PetscFindInt(cone[c], numSubPoints[d-1], subpoints[d-1], &subc);CHKERRQ(ierr); 2631 if (subc >= 0) { 2632 coneNew[coneSizeNew] = firstSubPoint[d-1] + subc; 2633 coneONew[coneSizeNew] = ornt[c]; 2634 ++coneSizeNew; 2635 } 2636 } 2637 if (coneSizeNew != subconeSize) SETERRQ2(comm, PETSC_ERR_PLIB, "Number of cone points located %d does not match subcone size %d", coneSizeNew, subconeSize); 2638 ierr = DMPlexSetCone(subdm, subpoint, coneNew);CHKERRQ(ierr); 2639 ierr = DMPlexSetConeOrientation(subdm, subpoint, coneONew);CHKERRQ(ierr); 2640 } 2641 } 2642 ierr = PetscFree2(coneNew,coneONew);CHKERRQ(ierr); 2643 ierr = DMPlexSymmetrize(subdm);CHKERRQ(ierr); 2644 ierr = DMPlexStratify(subdm);CHKERRQ(ierr); 2645 /* Build coordinates */ 2646 { 2647 PetscSection coordSection, subCoordSection; 2648 Vec coordinates, subCoordinates; 2649 PetscScalar *coords, *subCoords; 2650 PetscInt numComp, coordSize; 2651 2652 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 2653 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr); 2654 ierr = DMGetCoordinateSection(subdm, &subCoordSection);CHKERRQ(ierr); 2655 ierr = PetscSectionSetNumFields(subCoordSection, 1);CHKERRQ(ierr); 2656 ierr = PetscSectionGetFieldComponents(coordSection, 0, &numComp);CHKERRQ(ierr); 2657 ierr = PetscSectionSetFieldComponents(subCoordSection, 0, numComp);CHKERRQ(ierr); 2658 ierr = PetscSectionSetChart(subCoordSection, firstSubPoint[0], firstSubPoint[0]+numSubPoints[0]);CHKERRQ(ierr); 2659 for (v = 0; v < numSubPoints[0]; ++v) { 2660 const PetscInt vertex = subpoints[0][v]; 2661 const PetscInt subvertex = firstSubPoint[0]+v; 2662 PetscInt dof; 2663 2664 ierr = PetscSectionGetDof(coordSection, vertex, &dof);CHKERRQ(ierr); 2665 ierr = PetscSectionSetDof(subCoordSection, subvertex, dof);CHKERRQ(ierr); 2666 ierr = PetscSectionSetFieldDof(subCoordSection, subvertex, 0, dof);CHKERRQ(ierr); 2667 } 2668 ierr = PetscSectionSetUp(subCoordSection);CHKERRQ(ierr); 2669 ierr = PetscSectionGetStorageSize(subCoordSection, &coordSize);CHKERRQ(ierr); 2670 ierr = VecCreate(comm, &subCoordinates);CHKERRQ(ierr); 2671 ierr = VecSetSizes(subCoordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 2672 ierr = VecSetType(subCoordinates,VECSTANDARD);CHKERRQ(ierr); 2673 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 2674 ierr = VecGetArray(subCoordinates, &subCoords);CHKERRQ(ierr); 2675 for (v = 0; v < numSubPoints[0]; ++v) { 2676 const PetscInt vertex = subpoints[0][v]; 2677 const PetscInt subvertex = firstSubPoint[0]+v; 2678 PetscInt dof, off, sdof, soff, d; 2679 2680 ierr = PetscSectionGetDof(coordSection, vertex, &dof);CHKERRQ(ierr); 2681 ierr = PetscSectionGetOffset(coordSection, vertex, &off);CHKERRQ(ierr); 2682 ierr = PetscSectionGetDof(subCoordSection, subvertex, &sdof);CHKERRQ(ierr); 2683 ierr = PetscSectionGetOffset(subCoordSection, subvertex, &soff);CHKERRQ(ierr); 2684 if (dof != sdof) SETERRQ4(comm, PETSC_ERR_PLIB, "Coordinate dimension %d on subvertex %d, vertex %d should be %d", sdof, subvertex, vertex, dof); 2685 for (d = 0; d < dof; ++d) subCoords[soff+d] = coords[off+d]; 2686 } 2687 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 2688 ierr = VecRestoreArray(subCoordinates, &subCoords);CHKERRQ(ierr); 2689 ierr = DMSetCoordinatesLocal(subdm, subCoordinates);CHKERRQ(ierr); 2690 ierr = VecDestroy(&subCoordinates);CHKERRQ(ierr); 2691 } 2692 /* Cleanup */ 2693 for (d = 0; d <= dim; ++d) { 2694 if (subpointIS[d]) {ierr = ISRestoreIndices(subpointIS[d], &subpoints[d]);CHKERRQ(ierr);} 2695 ierr = ISDestroy(&subpointIS[d]);CHKERRQ(ierr); 2696 } 2697 ierr = PetscFree4(numSubPoints,firstSubPoint,subpointIS,subpoints);CHKERRQ(ierr); 2698 PetscFunctionReturn(0); 2699 } 2700 2701 #undef __FUNCT__ 2702 #define __FUNCT__ "DMPlexCreateSubmesh" 2703 /*@ 2704 DMPlexCreateSubmesh - Extract a hypersurface from the mesh using vertices defined by a label 2705 2706 Input Parameters: 2707 + dm - The original mesh 2708 . vertexLabel - The DMLabel marking vertices contained in the surface 2709 - value - The label value to use 2710 2711 Output Parameter: 2712 . subdm - The surface mesh 2713 2714 Note: This function produces a DMLabel mapping original points in the submesh to their depth. This can be obtained using DMPlexGetSubpointMap(). 2715 2716 Level: developer 2717 2718 .seealso: DMPlexGetSubpointMap(), DMPlexGetLabel(), DMLabelSetValue() 2719 @*/ 2720 PetscErrorCode DMPlexCreateSubmesh(DM dm, DMLabel vertexLabel, PetscInt value, DM *subdm) 2721 { 2722 PetscInt dim, depth; 2723 PetscErrorCode ierr; 2724 2725 PetscFunctionBegin; 2726 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 2727 PetscValidPointer(subdm, 3); 2728 ierr = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr); 2729 ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr); 2730 ierr = DMCreate(PetscObjectComm((PetscObject)dm), subdm);CHKERRQ(ierr); 2731 ierr = DMSetType(*subdm, DMPLEX);CHKERRQ(ierr); 2732 ierr = DMPlexSetDimension(*subdm, dim-1);CHKERRQ(ierr); 2733 if (depth == dim) { 2734 ierr = DMPlexCreateSubmesh_Interpolated(dm, vertexLabel, value, *subdm);CHKERRQ(ierr); 2735 } else { 2736 ierr = DMPlexCreateSubmesh_Uninterpolated(dm, vertexLabel, value, *subdm);CHKERRQ(ierr); 2737 } 2738 PetscFunctionReturn(0); 2739 } 2740 2741 #undef __FUNCT__ 2742 #define __FUNCT__ "DMPlexFilterPoint_Internal" 2743 PETSC_STATIC_INLINE PetscInt DMPlexFilterPoint_Internal(PetscInt point, PetscInt firstSubPoint, PetscInt numSubPoints, const PetscInt subPoints[]) 2744 { 2745 PetscInt subPoint; 2746 PetscErrorCode ierr; 2747 2748 ierr = PetscFindInt(point, numSubPoints, subPoints, &subPoint); if (ierr < 0) return ierr; 2749 return subPoint < 0 ? subPoint : firstSubPoint+subPoint; 2750 } 2751 2752 #undef __FUNCT__ 2753 #define __FUNCT__ "DMPlexCreateCohesiveSubmesh_Uninterpolated" 2754 static PetscErrorCode DMPlexCreateCohesiveSubmesh_Uninterpolated(DM dm, PetscBool hasLagrange, const char label[], PetscInt value, DM subdm) 2755 { 2756 MPI_Comm comm; 2757 DMLabel subpointMap; 2758 IS subvertexIS; 2759 const PetscInt *subVertices; 2760 PetscInt numSubVertices, firstSubVertex, numSubCells, *subCells = NULL; 2761 PetscInt *subface, maxConeSize, numSubFaces, firstSubFace, newFacePoint, nFV; 2762 PetscInt cMax, c, f; 2763 PetscErrorCode ierr; 2764 2765 PetscFunctionBegin; 2766 ierr = PetscObjectGetComm((PetscObject)dm, &comm);CHKERRQ(ierr); 2767 /* Create subpointMap which marks the submesh */ 2768 ierr = DMLabelCreate("subpoint_map", &subpointMap);CHKERRQ(ierr); 2769 ierr = DMPlexSetSubpointMap(subdm, subpointMap);CHKERRQ(ierr); 2770 ierr = DMLabelDestroy(&subpointMap);CHKERRQ(ierr); 2771 ierr = DMPlexMarkCohesiveSubmesh_Uninterpolated(dm, hasLagrange, label, value, subpointMap, &numSubFaces, &nFV, &subCells, subdm);CHKERRQ(ierr); 2772 /* Setup chart */ 2773 ierr = DMLabelGetStratumSize(subpointMap, 0, &numSubVertices);CHKERRQ(ierr); 2774 ierr = DMLabelGetStratumSize(subpointMap, 2, &numSubCells);CHKERRQ(ierr); 2775 ierr = DMPlexSetChart(subdm, 0, numSubCells+numSubFaces+numSubVertices);CHKERRQ(ierr); 2776 ierr = DMPlexSetVTKCellHeight(subdm, 1);CHKERRQ(ierr); 2777 /* Set cone sizes */ 2778 firstSubVertex = numSubCells; 2779 firstSubFace = numSubCells+numSubVertices; 2780 newFacePoint = firstSubFace; 2781 ierr = DMLabelGetStratumIS(subpointMap, 0, &subvertexIS);CHKERRQ(ierr); 2782 if (subvertexIS) {ierr = ISGetIndices(subvertexIS, &subVertices);CHKERRQ(ierr);} 2783 for (c = 0; c < numSubCells; ++c) { 2784 ierr = DMPlexSetConeSize(subdm, c, 1);CHKERRQ(ierr); 2785 } 2786 for (f = firstSubFace; f < firstSubFace+numSubFaces; ++f) { 2787 ierr = DMPlexSetConeSize(subdm, f, nFV);CHKERRQ(ierr); 2788 } 2789 ierr = DMSetUp(subdm);CHKERRQ(ierr); 2790 /* Create face cones */ 2791 ierr = DMPlexGetMaxSizes(dm, &maxConeSize, NULL);CHKERRQ(ierr); 2792 ierr = DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);CHKERRQ(ierr); 2793 ierr = DMGetWorkArray(subdm, maxConeSize, PETSC_INT, (void**) &subface);CHKERRQ(ierr); 2794 for (c = 0; c < numSubCells; ++c) { 2795 const PetscInt cell = subCells[c]; 2796 const PetscInt subcell = c; 2797 const PetscInt *cone, *cells; 2798 PetscInt numCells, subVertex, p, v; 2799 2800 if (cell < cMax) continue; 2801 ierr = DMPlexGetCone(dm, cell, &cone);CHKERRQ(ierr); 2802 for (v = 0; v < nFV; ++v) { 2803 ierr = PetscFindInt(cone[v], numSubVertices, subVertices, &subVertex);CHKERRQ(ierr); 2804 subface[v] = firstSubVertex+subVertex; 2805 } 2806 ierr = DMPlexSetCone(subdm, newFacePoint, subface);CHKERRQ(ierr); 2807 ierr = DMPlexSetCone(subdm, subcell, &newFacePoint);CHKERRQ(ierr); 2808 ierr = DMPlexGetJoin(dm, nFV, cone, &numCells, &cells);CHKERRQ(ierr); 2809 /* Not true in parallel 2810 if (numCells != 2) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cohesive cells should separate two cells"); */ 2811 for (p = 0; p < numCells; ++p) { 2812 PetscInt negsubcell; 2813 2814 if (cells[p] >= cMax) continue; 2815 /* I know this is a crap search */ 2816 for (negsubcell = 0; negsubcell < numSubCells; ++negsubcell) { 2817 if (subCells[negsubcell] == cells[p]) break; 2818 } 2819 if (negsubcell == numSubCells) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Could not find negative face neighbor for cohesive cell %d", cell); 2820 ierr = DMPlexSetCone(subdm, negsubcell, &newFacePoint);CHKERRQ(ierr); 2821 } 2822 ierr = DMPlexRestoreJoin(dm, nFV, cone, &numCells, &cells);CHKERRQ(ierr); 2823 ++newFacePoint; 2824 } 2825 ierr = DMRestoreWorkArray(subdm, maxConeSize, PETSC_INT, (void**) &subface);CHKERRQ(ierr); 2826 ierr = DMPlexSymmetrize(subdm);CHKERRQ(ierr); 2827 ierr = DMPlexStratify(subdm);CHKERRQ(ierr); 2828 /* Build coordinates */ 2829 { 2830 PetscSection coordSection, subCoordSection; 2831 Vec coordinates, subCoordinates; 2832 PetscScalar *coords, *subCoords; 2833 PetscInt numComp, coordSize, v; 2834 2835 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 2836 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr); 2837 ierr = DMGetCoordinateSection(subdm, &subCoordSection);CHKERRQ(ierr); 2838 ierr = PetscSectionSetNumFields(subCoordSection, 1);CHKERRQ(ierr); 2839 ierr = PetscSectionGetFieldComponents(coordSection, 0, &numComp);CHKERRQ(ierr); 2840 ierr = PetscSectionSetFieldComponents(subCoordSection, 0, numComp);CHKERRQ(ierr); 2841 ierr = PetscSectionSetChart(subCoordSection, firstSubVertex, firstSubVertex+numSubVertices);CHKERRQ(ierr); 2842 for (v = 0; v < numSubVertices; ++v) { 2843 const PetscInt vertex = subVertices[v]; 2844 const PetscInt subvertex = firstSubVertex+v; 2845 PetscInt dof; 2846 2847 ierr = PetscSectionGetDof(coordSection, vertex, &dof);CHKERRQ(ierr); 2848 ierr = PetscSectionSetDof(subCoordSection, subvertex, dof);CHKERRQ(ierr); 2849 ierr = PetscSectionSetFieldDof(subCoordSection, subvertex, 0, dof);CHKERRQ(ierr); 2850 } 2851 ierr = PetscSectionSetUp(subCoordSection);CHKERRQ(ierr); 2852 ierr = PetscSectionGetStorageSize(subCoordSection, &coordSize);CHKERRQ(ierr); 2853 ierr = VecCreate(comm, &subCoordinates);CHKERRQ(ierr); 2854 ierr = VecSetSizes(subCoordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 2855 ierr = VecSetType(subCoordinates,VECSTANDARD);CHKERRQ(ierr); 2856 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 2857 ierr = VecGetArray(subCoordinates, &subCoords);CHKERRQ(ierr); 2858 for (v = 0; v < numSubVertices; ++v) { 2859 const PetscInt vertex = subVertices[v]; 2860 const PetscInt subvertex = firstSubVertex+v; 2861 PetscInt dof, off, sdof, soff, d; 2862 2863 ierr = PetscSectionGetDof(coordSection, vertex, &dof);CHKERRQ(ierr); 2864 ierr = PetscSectionGetOffset(coordSection, vertex, &off);CHKERRQ(ierr); 2865 ierr = PetscSectionGetDof(subCoordSection, subvertex, &sdof);CHKERRQ(ierr); 2866 ierr = PetscSectionGetOffset(subCoordSection, subvertex, &soff);CHKERRQ(ierr); 2867 if (dof != sdof) SETERRQ4(comm, PETSC_ERR_PLIB, "Coordinate dimension %d on subvertex %d, vertex %d should be %d", sdof, subvertex, vertex, dof); 2868 for (d = 0; d < dof; ++d) subCoords[soff+d] = coords[off+d]; 2869 } 2870 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 2871 ierr = VecRestoreArray(subCoordinates, &subCoords);CHKERRQ(ierr); 2872 ierr = DMSetCoordinatesLocal(subdm, subCoordinates);CHKERRQ(ierr); 2873 ierr = VecDestroy(&subCoordinates);CHKERRQ(ierr); 2874 } 2875 /* Build SF */ 2876 CHKMEMQ; 2877 { 2878 PetscSF sfPoint, sfPointSub; 2879 const PetscSFNode *remotePoints; 2880 PetscSFNode *sremotePoints, *newLocalPoints, *newOwners; 2881 const PetscInt *localPoints; 2882 PetscInt *slocalPoints; 2883 PetscInt numRoots, numLeaves, numSubRoots = numSubCells+numSubFaces+numSubVertices, numSubLeaves = 0, l, sl, ll, pStart, pEnd, p, vStart, vEnd; 2884 PetscMPIInt rank; 2885 2886 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank);CHKERRQ(ierr); 2887 ierr = DMGetPointSF(dm, &sfPoint);CHKERRQ(ierr); 2888 ierr = DMGetPointSF(subdm, &sfPointSub);CHKERRQ(ierr); 2889 ierr = DMPlexGetChart(dm, &pStart, &pEnd);CHKERRQ(ierr); 2890 ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr); 2891 ierr = PetscSFGetGraph(sfPoint, &numRoots, &numLeaves, &localPoints, &remotePoints);CHKERRQ(ierr); 2892 if (numRoots >= 0) { 2893 /* Only vertices should be shared */ 2894 ierr = PetscMalloc2(pEnd-pStart,&newLocalPoints,numRoots,&newOwners);CHKERRQ(ierr); 2895 for (p = 0; p < pEnd-pStart; ++p) { 2896 newLocalPoints[p].rank = -2; 2897 newLocalPoints[p].index = -2; 2898 } 2899 /* Set subleaves */ 2900 for (l = 0; l < numLeaves; ++l) { 2901 const PetscInt point = localPoints[l]; 2902 const PetscInt subPoint = DMPlexFilterPoint_Internal(point, firstSubVertex, numSubVertices, subVertices); 2903 2904 if ((point < vStart) && (point >= vEnd)) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Should not be mapping anything but vertices, %d", point); 2905 if (subPoint < 0) continue; 2906 newLocalPoints[point-pStart].rank = rank; 2907 newLocalPoints[point-pStart].index = subPoint; 2908 ++numSubLeaves; 2909 } 2910 /* Must put in owned subpoints */ 2911 for (p = pStart; p < pEnd; ++p) { 2912 const PetscInt subPoint = DMPlexFilterPoint_Internal(p, firstSubVertex, numSubVertices, subVertices); 2913 2914 if (subPoint < 0) { 2915 newOwners[p-pStart].rank = -3; 2916 newOwners[p-pStart].index = -3; 2917 } else { 2918 newOwners[p-pStart].rank = rank; 2919 newOwners[p-pStart].index = subPoint; 2920 } 2921 } 2922 ierr = PetscSFReduceBegin(sfPoint, MPIU_2INT, newLocalPoints, newOwners, MPI_MAXLOC);CHKERRQ(ierr); 2923 ierr = PetscSFReduceEnd(sfPoint, MPIU_2INT, newLocalPoints, newOwners, MPI_MAXLOC);CHKERRQ(ierr); 2924 ierr = PetscSFBcastBegin(sfPoint, MPIU_2INT, newOwners, newLocalPoints);CHKERRQ(ierr); 2925 ierr = PetscSFBcastEnd(sfPoint, MPIU_2INT, newOwners, newLocalPoints);CHKERRQ(ierr); 2926 ierr = PetscMalloc1(numSubLeaves, &slocalPoints);CHKERRQ(ierr); 2927 ierr = PetscMalloc1(numSubLeaves, &sremotePoints);CHKERRQ(ierr); 2928 for (l = 0, sl = 0, ll = 0; l < numLeaves; ++l) { 2929 const PetscInt point = localPoints[l]; 2930 const PetscInt subPoint = DMPlexFilterPoint_Internal(point, firstSubVertex, numSubVertices, subVertices); 2931 2932 if (subPoint < 0) continue; 2933 if (newLocalPoints[point].rank == rank) {++ll; continue;} 2934 slocalPoints[sl] = subPoint; 2935 sremotePoints[sl].rank = newLocalPoints[point].rank; 2936 sremotePoints[sl].index = newLocalPoints[point].index; 2937 if (sremotePoints[sl].rank < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid remote rank for local point %d", point); 2938 if (sremotePoints[sl].index < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid remote subpoint for local point %d", point); 2939 ++sl; 2940 } 2941 ierr = PetscFree2(newLocalPoints,newOwners);CHKERRQ(ierr); 2942 if (sl + ll != numSubLeaves) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Mismatch in number of subleaves %d + %d != %d", sl, ll, numSubLeaves); 2943 ierr = PetscSFSetGraph(sfPointSub, numSubRoots, sl, slocalPoints, PETSC_OWN_POINTER, sremotePoints, PETSC_OWN_POINTER);CHKERRQ(ierr); 2944 } 2945 } 2946 CHKMEMQ; 2947 /* Cleanup */ 2948 if (subvertexIS) {ierr = ISRestoreIndices(subvertexIS, &subVertices);CHKERRQ(ierr);} 2949 ierr = ISDestroy(&subvertexIS);CHKERRQ(ierr); 2950 ierr = PetscFree(subCells);CHKERRQ(ierr); 2951 PetscFunctionReturn(0); 2952 } 2953 2954 #undef __FUNCT__ 2955 #define __FUNCT__ "DMPlexCreateCohesiveSubmesh_Interpolated" 2956 static PetscErrorCode DMPlexCreateCohesiveSubmesh_Interpolated(DM dm, const char label[], PetscInt value, DM subdm) 2957 { 2958 MPI_Comm comm; 2959 DMLabel subpointMap; 2960 IS *subpointIS; 2961 const PetscInt **subpoints; 2962 PetscInt *numSubPoints, *firstSubPoint, *coneNew, *orntNew; 2963 PetscInt totSubPoints = 0, maxConeSize, dim, p, d, v; 2964 PetscErrorCode ierr; 2965 2966 PetscFunctionBegin; 2967 ierr = PetscObjectGetComm((PetscObject)dm,&comm);CHKERRQ(ierr); 2968 /* Create subpointMap which marks the submesh */ 2969 ierr = DMLabelCreate("subpoint_map", &subpointMap);CHKERRQ(ierr); 2970 ierr = DMPlexSetSubpointMap(subdm, subpointMap);CHKERRQ(ierr); 2971 ierr = DMLabelDestroy(&subpointMap);CHKERRQ(ierr); 2972 ierr = DMPlexMarkCohesiveSubmesh_Interpolated(dm, label, value, subpointMap, subdm);CHKERRQ(ierr); 2973 /* Setup chart */ 2974 ierr = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr); 2975 ierr = PetscMalloc4(dim+1,&numSubPoints,dim+1,&firstSubPoint,dim+1,&subpointIS,dim+1,&subpoints);CHKERRQ(ierr); 2976 for (d = 0; d <= dim; ++d) { 2977 ierr = DMLabelGetStratumSize(subpointMap, d, &numSubPoints[d]);CHKERRQ(ierr); 2978 totSubPoints += numSubPoints[d]; 2979 } 2980 ierr = DMPlexSetChart(subdm, 0, totSubPoints);CHKERRQ(ierr); 2981 ierr = DMPlexSetVTKCellHeight(subdm, 1);CHKERRQ(ierr); 2982 /* Set cone sizes */ 2983 firstSubPoint[dim] = 0; 2984 firstSubPoint[0] = firstSubPoint[dim] + numSubPoints[dim]; 2985 if (dim > 1) {firstSubPoint[dim-1] = firstSubPoint[0] + numSubPoints[0];} 2986 if (dim > 2) {firstSubPoint[dim-2] = firstSubPoint[dim-1] + numSubPoints[dim-1];} 2987 for (d = 0; d <= dim; ++d) { 2988 ierr = DMLabelGetStratumIS(subpointMap, d, &subpointIS[d]);CHKERRQ(ierr); 2989 if (subpointIS[d]) {ierr = ISGetIndices(subpointIS[d], &subpoints[d]);CHKERRQ(ierr);} 2990 } 2991 for (d = 0; d <= dim; ++d) { 2992 for (p = 0; p < numSubPoints[d]; ++p) { 2993 const PetscInt point = subpoints[d][p]; 2994 const PetscInt subpoint = firstSubPoint[d] + p; 2995 const PetscInt *cone; 2996 PetscInt coneSize, coneSizeNew, c, val; 2997 2998 ierr = DMPlexGetConeSize(dm, point, &coneSize);CHKERRQ(ierr); 2999 ierr = DMPlexSetConeSize(subdm, subpoint, coneSize);CHKERRQ(ierr); 3000 if (d == dim) { 3001 ierr = DMPlexGetCone(dm, point, &cone);CHKERRQ(ierr); 3002 for (c = 0, coneSizeNew = 0; c < coneSize; ++c) { 3003 ierr = DMLabelGetValue(subpointMap, cone[c], &val);CHKERRQ(ierr); 3004 if (val >= 0) coneSizeNew++; 3005 } 3006 ierr = DMPlexSetConeSize(subdm, subpoint, coneSizeNew);CHKERRQ(ierr); 3007 } 3008 } 3009 } 3010 ierr = DMSetUp(subdm);CHKERRQ(ierr); 3011 /* Set cones */ 3012 ierr = DMPlexGetMaxSizes(dm, &maxConeSize, NULL);CHKERRQ(ierr); 3013 ierr = PetscMalloc2(maxConeSize,&coneNew,maxConeSize,&orntNew);CHKERRQ(ierr); 3014 for (d = 0; d <= dim; ++d) { 3015 for (p = 0; p < numSubPoints[d]; ++p) { 3016 const PetscInt point = subpoints[d][p]; 3017 const PetscInt subpoint = firstSubPoint[d] + p; 3018 const PetscInt *cone, *ornt; 3019 PetscInt coneSize, subconeSize, coneSizeNew, c, subc; 3020 3021 ierr = DMPlexGetConeSize(dm, point, &coneSize);CHKERRQ(ierr); 3022 ierr = DMPlexGetConeSize(subdm, subpoint, &subconeSize);CHKERRQ(ierr); 3023 ierr = DMPlexGetCone(dm, point, &cone);CHKERRQ(ierr); 3024 ierr = DMPlexGetConeOrientation(dm, point, &ornt);CHKERRQ(ierr); 3025 for (c = 0, coneSizeNew = 0; c < coneSize; ++c) { 3026 ierr = PetscFindInt(cone[c], numSubPoints[d-1], subpoints[d-1], &subc);CHKERRQ(ierr); 3027 if (subc >= 0) { 3028 coneNew[coneSizeNew] = firstSubPoint[d-1] + subc; 3029 orntNew[coneSizeNew++] = ornt[c]; 3030 } 3031 } 3032 if (coneSizeNew != subconeSize) SETERRQ2(comm, PETSC_ERR_PLIB, "Number of cone points located %d does not match subcone size %d", coneSizeNew, subconeSize); 3033 ierr = DMPlexSetCone(subdm, subpoint, coneNew);CHKERRQ(ierr); 3034 ierr = DMPlexSetConeOrientation(subdm, subpoint, orntNew);CHKERRQ(ierr); 3035 } 3036 } 3037 ierr = PetscFree2(coneNew,orntNew);CHKERRQ(ierr); 3038 ierr = DMPlexSymmetrize(subdm);CHKERRQ(ierr); 3039 ierr = DMPlexStratify(subdm);CHKERRQ(ierr); 3040 /* Build coordinates */ 3041 { 3042 PetscSection coordSection, subCoordSection; 3043 Vec coordinates, subCoordinates; 3044 PetscScalar *coords, *subCoords; 3045 PetscInt numComp, coordSize; 3046 3047 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 3048 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr); 3049 ierr = DMGetCoordinateSection(subdm, &subCoordSection);CHKERRQ(ierr); 3050 ierr = PetscSectionSetNumFields(subCoordSection, 1);CHKERRQ(ierr); 3051 ierr = PetscSectionGetFieldComponents(coordSection, 0, &numComp);CHKERRQ(ierr); 3052 ierr = PetscSectionSetFieldComponents(subCoordSection, 0, numComp);CHKERRQ(ierr); 3053 ierr = PetscSectionSetChart(subCoordSection, firstSubPoint[0], firstSubPoint[0]+numSubPoints[0]);CHKERRQ(ierr); 3054 for (v = 0; v < numSubPoints[0]; ++v) { 3055 const PetscInt vertex = subpoints[0][v]; 3056 const PetscInt subvertex = firstSubPoint[0]+v; 3057 PetscInt dof; 3058 3059 ierr = PetscSectionGetDof(coordSection, vertex, &dof);CHKERRQ(ierr); 3060 ierr = PetscSectionSetDof(subCoordSection, subvertex, dof);CHKERRQ(ierr); 3061 ierr = PetscSectionSetFieldDof(subCoordSection, subvertex, 0, dof);CHKERRQ(ierr); 3062 } 3063 ierr = PetscSectionSetUp(subCoordSection);CHKERRQ(ierr); 3064 ierr = PetscSectionGetStorageSize(subCoordSection, &coordSize);CHKERRQ(ierr); 3065 ierr = VecCreate(comm, &subCoordinates);CHKERRQ(ierr); 3066 ierr = VecSetSizes(subCoordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 3067 ierr = VecSetType(subCoordinates,VECSTANDARD);CHKERRQ(ierr); 3068 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 3069 ierr = VecGetArray(subCoordinates, &subCoords);CHKERRQ(ierr); 3070 for (v = 0; v < numSubPoints[0]; ++v) { 3071 const PetscInt vertex = subpoints[0][v]; 3072 const PetscInt subvertex = firstSubPoint[0]+v; 3073 PetscInt dof, off, sdof, soff, d; 3074 3075 ierr = PetscSectionGetDof(coordSection, vertex, &dof);CHKERRQ(ierr); 3076 ierr = PetscSectionGetOffset(coordSection, vertex, &off);CHKERRQ(ierr); 3077 ierr = PetscSectionGetDof(subCoordSection, subvertex, &sdof);CHKERRQ(ierr); 3078 ierr = PetscSectionGetOffset(subCoordSection, subvertex, &soff);CHKERRQ(ierr); 3079 if (dof != sdof) SETERRQ4(comm, PETSC_ERR_PLIB, "Coordinate dimension %d on subvertex %d, vertex %d should be %d", sdof, subvertex, vertex, dof); 3080 for (d = 0; d < dof; ++d) subCoords[soff+d] = coords[off+d]; 3081 } 3082 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 3083 ierr = VecRestoreArray(subCoordinates, &subCoords);CHKERRQ(ierr); 3084 ierr = DMSetCoordinatesLocal(subdm, subCoordinates);CHKERRQ(ierr); 3085 ierr = VecDestroy(&subCoordinates);CHKERRQ(ierr); 3086 } 3087 /* Build SF: We need this complexity because subpoints might not be selected on the owning process */ 3088 { 3089 PetscSF sfPoint, sfPointSub; 3090 IS subpIS; 3091 const PetscSFNode *remotePoints; 3092 PetscSFNode *sremotePoints, *newLocalPoints, *newOwners; 3093 const PetscInt *localPoints, *subpoints; 3094 PetscInt *slocalPoints; 3095 PetscInt numRoots, numLeaves, numSubpoints = 0, numSubroots, numSubleaves = 0, l, sl, ll, pStart, pEnd, p; 3096 PetscMPIInt rank; 3097 3098 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank);CHKERRQ(ierr); 3099 ierr = DMGetPointSF(dm, &sfPoint);CHKERRQ(ierr); 3100 ierr = DMGetPointSF(subdm, &sfPointSub);CHKERRQ(ierr); 3101 ierr = DMPlexGetChart(dm, &pStart, &pEnd);CHKERRQ(ierr); 3102 ierr = DMPlexGetChart(subdm, NULL, &numSubroots);CHKERRQ(ierr); 3103 ierr = DMPlexCreateSubpointIS(subdm, &subpIS);CHKERRQ(ierr); 3104 if (subpIS) { 3105 ierr = ISGetIndices(subpIS, &subpoints);CHKERRQ(ierr); 3106 ierr = ISGetLocalSize(subpIS, &numSubpoints);CHKERRQ(ierr); 3107 } 3108 ierr = PetscSFGetGraph(sfPoint, &numRoots, &numLeaves, &localPoints, &remotePoints);CHKERRQ(ierr); 3109 if (numRoots >= 0) { 3110 ierr = PetscMalloc2(pEnd-pStart,&newLocalPoints,numRoots,&newOwners);CHKERRQ(ierr); 3111 for (p = 0; p < pEnd-pStart; ++p) { 3112 newLocalPoints[p].rank = -2; 3113 newLocalPoints[p].index = -2; 3114 } 3115 /* Set subleaves */ 3116 for (l = 0; l < numLeaves; ++l) { 3117 const PetscInt point = localPoints[l]; 3118 const PetscInt subpoint = DMPlexFilterPoint_Internal(point, 0, numSubpoints, subpoints); 3119 3120 if (subpoint < 0) continue; 3121 newLocalPoints[point-pStart].rank = rank; 3122 newLocalPoints[point-pStart].index = subpoint; 3123 ++numSubleaves; 3124 } 3125 /* Must put in owned subpoints */ 3126 for (p = pStart; p < pEnd; ++p) { 3127 const PetscInt subpoint = DMPlexFilterPoint_Internal(p, 0, numSubpoints, subpoints); 3128 3129 if (subpoint < 0) { 3130 newOwners[p-pStart].rank = -3; 3131 newOwners[p-pStart].index = -3; 3132 } else { 3133 newOwners[p-pStart].rank = rank; 3134 newOwners[p-pStart].index = subpoint; 3135 } 3136 } 3137 ierr = PetscSFReduceBegin(sfPoint, MPIU_2INT, newLocalPoints, newOwners, MPI_MAXLOC);CHKERRQ(ierr); 3138 ierr = PetscSFReduceEnd(sfPoint, MPIU_2INT, newLocalPoints, newOwners, MPI_MAXLOC);CHKERRQ(ierr); 3139 ierr = PetscSFBcastBegin(sfPoint, MPIU_2INT, newOwners, newLocalPoints);CHKERRQ(ierr); 3140 ierr = PetscSFBcastEnd(sfPoint, MPIU_2INT, newOwners, newLocalPoints);CHKERRQ(ierr); 3141 ierr = PetscMalloc(numSubleaves * sizeof(PetscInt), &slocalPoints);CHKERRQ(ierr); 3142 ierr = PetscMalloc(numSubleaves * sizeof(PetscSFNode), &sremotePoints);CHKERRQ(ierr); 3143 for (l = 0, sl = 0, ll = 0; l < numLeaves; ++l) { 3144 const PetscInt point = localPoints[l]; 3145 const PetscInt subpoint = DMPlexFilterPoint_Internal(point, 0, numSubpoints, subpoints); 3146 3147 if (subpoint < 0) continue; 3148 if (newLocalPoints[point].rank == rank) {++ll; continue;} 3149 slocalPoints[sl] = subpoint; 3150 sremotePoints[sl].rank = newLocalPoints[point].rank; 3151 sremotePoints[sl].index = newLocalPoints[point].index; 3152 if (sremotePoints[sl].rank < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid remote rank for local point %d", point); 3153 if (sremotePoints[sl].index < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid remote subpoint for local point %d", point); 3154 ++sl; 3155 } 3156 if (sl + ll != numSubleaves) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Mismatch in number of subleaves %d + %d != %d", sl, ll, numSubleaves); 3157 ierr = PetscFree2(newLocalPoints,newOwners);CHKERRQ(ierr); 3158 ierr = PetscSFSetGraph(sfPointSub, numSubroots, sl, slocalPoints, PETSC_OWN_POINTER, sremotePoints, PETSC_OWN_POINTER);CHKERRQ(ierr); 3159 } 3160 if (subpIS) { 3161 ierr = ISRestoreIndices(subpIS, &subpoints);CHKERRQ(ierr); 3162 ierr = ISDestroy(&subpIS);CHKERRQ(ierr); 3163 } 3164 } 3165 /* Cleanup */ 3166 for (d = 0; d <= dim; ++d) { 3167 if (subpointIS[d]) {ierr = ISRestoreIndices(subpointIS[d], &subpoints[d]);CHKERRQ(ierr);} 3168 ierr = ISDestroy(&subpointIS[d]);CHKERRQ(ierr); 3169 } 3170 ierr = PetscFree4(numSubPoints,firstSubPoint,subpointIS,subpoints);CHKERRQ(ierr); 3171 PetscFunctionReturn(0); 3172 } 3173 3174 #undef __FUNCT__ 3175 #define __FUNCT__ "DMPlexCreateCohesiveSubmesh" 3176 /* 3177 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. 3178 3179 Input Parameters: 3180 + dm - The original mesh 3181 . hasLagrange - The mesh has Lagrange unknowns in the cohesive cells 3182 . label - A label name, or NULL 3183 - value - A label value 3184 3185 Output Parameter: 3186 . subdm - The surface mesh 3187 3188 Note: This function produces a DMLabel mapping original points in the submesh to their depth. This can be obtained using DMPlexGetSubpointMap(). 3189 3190 Level: developer 3191 3192 .seealso: DMPlexGetSubpointMap(), DMPlexCreateSubmesh() 3193 */ 3194 PetscErrorCode DMPlexCreateCohesiveSubmesh(DM dm, PetscBool hasLagrange, const char label[], PetscInt value, DM *subdm) 3195 { 3196 PetscInt dim, depth; 3197 PetscErrorCode ierr; 3198 3199 PetscFunctionBegin; 3200 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 3201 PetscValidPointer(subdm, 5); 3202 ierr = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr); 3203 ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr); 3204 ierr = DMCreate(PetscObjectComm((PetscObject)dm), subdm);CHKERRQ(ierr); 3205 ierr = DMSetType(*subdm, DMPLEX);CHKERRQ(ierr); 3206 ierr = DMPlexSetDimension(*subdm, dim-1);CHKERRQ(ierr); 3207 if (depth == dim) { 3208 ierr = DMPlexCreateCohesiveSubmesh_Interpolated(dm, label, value, *subdm);CHKERRQ(ierr); 3209 } else { 3210 ierr = DMPlexCreateCohesiveSubmesh_Uninterpolated(dm, hasLagrange, label, value, *subdm);CHKERRQ(ierr); 3211 } 3212 PetscFunctionReturn(0); 3213 } 3214 3215 #undef __FUNCT__ 3216 #define __FUNCT__ "DMPlexGetSubpointMap" 3217 /*@ 3218 DMPlexGetSubpointMap - Returns a DMLabel with point dimension as values 3219 3220 Input Parameter: 3221 . dm - The submesh DM 3222 3223 Output Parameter: 3224 . subpointMap - The DMLabel of all the points from the original mesh in this submesh, or NULL if this is not a submesh 3225 3226 Level: developer 3227 3228 .seealso: DMPlexCreateSubmesh(), DMPlexCreateSubpointIS() 3229 @*/ 3230 PetscErrorCode DMPlexGetSubpointMap(DM dm, DMLabel *subpointMap) 3231 { 3232 DM_Plex *mesh = (DM_Plex*) dm->data; 3233 3234 PetscFunctionBegin; 3235 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 3236 PetscValidPointer(subpointMap, 2); 3237 *subpointMap = mesh->subpointMap; 3238 PetscFunctionReturn(0); 3239 } 3240 3241 #undef __FUNCT__ 3242 #define __FUNCT__ "DMPlexSetSubpointMap" 3243 /* Note: Should normally not be called by the user, since it is set in DMPlexCreateSubmesh() */ 3244 PetscErrorCode DMPlexSetSubpointMap(DM dm, DMLabel subpointMap) 3245 { 3246 DM_Plex *mesh = (DM_Plex *) dm->data; 3247 DMLabel tmp; 3248 PetscErrorCode ierr; 3249 3250 PetscFunctionBegin; 3251 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 3252 tmp = mesh->subpointMap; 3253 mesh->subpointMap = subpointMap; 3254 ++mesh->subpointMap->refct; 3255 ierr = DMLabelDestroy(&tmp);CHKERRQ(ierr); 3256 PetscFunctionReturn(0); 3257 } 3258 3259 #undef __FUNCT__ 3260 #define __FUNCT__ "DMPlexCreateSubpointIS" 3261 /*@ 3262 DMPlexCreateSubpointIS - Creates an IS covering the entire subdm chart with the original points as data 3263 3264 Input Parameter: 3265 . dm - The submesh DM 3266 3267 Output Parameter: 3268 . subpointIS - The IS of all the points from the original mesh in this submesh, or NULL if this is not a submesh 3269 3270 Note: This is IS is guaranteed to be sorted by the construction of the submesh 3271 3272 Level: developer 3273 3274 .seealso: DMPlexCreateSubmesh(), DMPlexGetSubpointMap() 3275 @*/ 3276 PetscErrorCode DMPlexCreateSubpointIS(DM dm, IS *subpointIS) 3277 { 3278 MPI_Comm comm; 3279 DMLabel subpointMap; 3280 IS is; 3281 const PetscInt *opoints; 3282 PetscInt *points, *depths; 3283 PetscInt depth, depStart, depEnd, d, pStart, pEnd, p, n, off; 3284 PetscErrorCode ierr; 3285 3286 PetscFunctionBegin; 3287 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 3288 PetscValidPointer(subpointIS, 2); 3289 ierr = PetscObjectGetComm((PetscObject)dm,&comm);CHKERRQ(ierr); 3290 *subpointIS = NULL; 3291 ierr = DMPlexGetSubpointMap(dm, &subpointMap);CHKERRQ(ierr); 3292 ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr); 3293 if (subpointMap && depth >= 0) { 3294 ierr = DMPlexGetChart(dm, &pStart, &pEnd);CHKERRQ(ierr); 3295 if (pStart) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Submeshes must start the point numbering at 0, not %d", pStart); 3296 ierr = DMGetWorkArray(dm, depth+1, PETSC_INT, &depths);CHKERRQ(ierr); 3297 depths[0] = depth; 3298 depths[1] = 0; 3299 for(d = 2; d <= depth; ++d) {depths[d] = depth+1 - d;} 3300 ierr = PetscMalloc1(pEnd, &points);CHKERRQ(ierr); 3301 for(d = 0, off = 0; d <= depth; ++d) { 3302 const PetscInt dep = depths[d]; 3303 3304 ierr = DMPlexGetDepthStratum(dm, dep, &depStart, &depEnd);CHKERRQ(ierr); 3305 ierr = DMLabelGetStratumSize(subpointMap, dep, &n);CHKERRQ(ierr); 3306 if (((d < 2) && (depth > 1)) || (d == 1)) { /* Only check vertices and cells for now since the map is broken for others */ 3307 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); 3308 } else { 3309 if (!n) { 3310 if (d == 0) { 3311 /* Missing cells */ 3312 for(p = 0; p < depEnd-depStart; ++p, ++off) points[off] = -1; 3313 } else { 3314 /* Missing faces */ 3315 for(p = 0; p < depEnd-depStart; ++p, ++off) points[off] = PETSC_MAX_INT; 3316 } 3317 } 3318 } 3319 if (n) { 3320 ierr = DMLabelGetStratumIS(subpointMap, dep, &is);CHKERRQ(ierr); 3321 ierr = ISGetIndices(is, &opoints);CHKERRQ(ierr); 3322 for(p = 0; p < n; ++p, ++off) points[off] = opoints[p]; 3323 ierr = ISRestoreIndices(is, &opoints);CHKERRQ(ierr); 3324 ierr = ISDestroy(&is);CHKERRQ(ierr); 3325 } 3326 } 3327 ierr = DMRestoreWorkArray(dm, depth+1, PETSC_INT, &depths);CHKERRQ(ierr); 3328 if (off != pEnd) SETERRQ2(comm, PETSC_ERR_ARG_WRONG, "The number of mapped submesh points %d should be %d", off, pEnd); 3329 ierr = ISCreateGeneral(PETSC_COMM_SELF, pEnd, points, PETSC_OWN_POINTER, subpointIS);CHKERRQ(ierr); 3330 } 3331 PetscFunctionReturn(0); 3332 } 3333