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