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