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