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