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