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