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