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