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