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