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