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