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