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