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