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