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