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