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