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