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