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