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