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