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