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