1 #include <petsc/private/dmpleximpl.h> /*I "petscdmplex.h" I*/ 2 #include <petscsf.h> 3 4 /*@ 5 DMPlexOrientPoint - Act with the given orientation on the cone points of this mesh point, and update its use in the mesh. 6 7 Not Collective 8 9 Input Parameters: 10 + dm - The `DM` 11 . p - The mesh point 12 - o - The orientation 13 14 Level: intermediate 15 16 .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexOrient()`, `DMPlexGetCone()`, `DMPlexGetConeOrientation()`, `DMPlexInterpolate()`, `DMPlexGetChart()` 17 @*/ 18 PetscErrorCode DMPlexOrientPoint(DM dm, PetscInt p, PetscInt o) 19 { 20 DMPolytopeType ct; 21 const PetscInt *arr, *cone, *ornt, *support; 22 PetscInt *newcone, *newornt; 23 PetscInt coneSize, c, supportSize, s; 24 25 PetscFunctionBegin; 26 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 27 PetscCall(DMPlexGetCellType(dm, p, &ct)); 28 arr = DMPolytopeTypeGetArrangement(ct, o); 29 if (!arr) PetscFunctionReturn(PETSC_SUCCESS); 30 PetscCall(DMPlexGetConeSize(dm, p, &coneSize)); 31 PetscCall(DMPlexGetCone(dm, p, &cone)); 32 PetscCall(DMPlexGetConeOrientation(dm, p, &ornt)); 33 PetscCall(DMGetWorkArray(dm, coneSize, MPIU_INT, &newcone)); 34 PetscCall(DMGetWorkArray(dm, coneSize, MPIU_INT, &newornt)); 35 for (c = 0; c < coneSize; ++c) { 36 DMPolytopeType ft; 37 PetscInt nO; 38 39 PetscCall(DMPlexGetCellType(dm, cone[c], &ft)); 40 nO = DMPolytopeTypeGetNumArrangements(ft) / 2; 41 newcone[c] = cone[arr[c * 2 + 0]]; 42 newornt[c] = DMPolytopeTypeComposeOrientation(ft, arr[c * 2 + 1], ornt[arr[c * 2 + 0]]); 43 PetscCheck(!newornt[c] || !(newornt[c] >= nO || newornt[c] < -nO), PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid orientation %" PetscInt_FMT " not in [%" PetscInt_FMT ",%" PetscInt_FMT ") for %s %" PetscInt_FMT, newornt[c], -nO, nO, DMPolytopeTypes[ft], cone[c]); 44 } 45 PetscCall(DMPlexSetCone(dm, p, newcone)); 46 PetscCall(DMPlexSetConeOrientation(dm, p, newornt)); 47 PetscCall(DMRestoreWorkArray(dm, coneSize, MPIU_INT, &newcone)); 48 PetscCall(DMRestoreWorkArray(dm, coneSize, MPIU_INT, &newornt)); 49 /* Update orientation of this point in the support points */ 50 PetscCall(DMPlexGetSupportSize(dm, p, &supportSize)); 51 PetscCall(DMPlexGetSupport(dm, p, &support)); 52 for (s = 0; s < supportSize; ++s) { 53 PetscCall(DMPlexGetConeSize(dm, support[s], &coneSize)); 54 PetscCall(DMPlexGetCone(dm, support[s], &cone)); 55 PetscCall(DMPlexGetConeOrientation(dm, support[s], &ornt)); 56 for (c = 0; c < coneSize; ++c) { 57 PetscInt po; 58 59 if (cone[c] != p) continue; 60 /* ornt[c] * 0 = target = po * o so that po = ornt[c] * o^{-1} */ 61 po = DMPolytopeTypeComposeOrientationInv(ct, ornt[c], o); 62 PetscCall(DMPlexInsertConeOrientation(dm, support[s], c, po)); 63 } 64 } 65 PetscFunctionReturn(PETSC_SUCCESS); 66 } 67 68 static PetscInt GetPointIndex(PetscInt point, PetscInt pStart, PetscInt pEnd, const PetscInt points[]) 69 { 70 if (points) { 71 PetscInt loc; 72 73 PetscCallAbort(PETSC_COMM_SELF, PetscFindInt(point, pEnd - pStart, points, &loc)); 74 if (loc >= 0) return loc; 75 } else { 76 if (point >= pStart && point < pEnd) return point - pStart; 77 } 78 return -1; 79 } 80 81 /* 82 - Checks face match 83 - Flips non-matching 84 - Inserts faces of support cells in FIFO 85 */ 86 static PetscErrorCode DMPlexCheckFace_Internal(DM dm, PetscInt *faceFIFO, PetscInt *fTop, PetscInt *fBottom, IS cellIS, IS faceIS, PetscBT seenCells, PetscBT flippedCells, PetscBT seenFaces) 87 { 88 const PetscInt *supp, *coneA, *coneB, *coneOA, *coneOB; 89 PetscInt suppSize, Ns = 0, coneSizeA, coneSizeB, posA = -1, posB = -1; 90 PetscInt face, dim, indC[3], indS[3], seenA, flippedA, seenB, flippedB, mismatch; 91 const PetscInt *cells, *faces; 92 PetscInt cStart, cEnd, fStart, fEnd; 93 94 PetscFunctionBegin; 95 face = faceFIFO[(*fTop)++]; 96 PetscCall(ISGetPointRange(cellIS, &cStart, &cEnd, &cells)); 97 PetscCall(ISGetPointRange(faceIS, &fStart, &fEnd, &faces)); 98 PetscCall(DMPlexGetPointDepth(dm, cells ? cells[cStart] : cStart, &dim)); 99 PetscCall(DMPlexGetSupportSize(dm, face, &suppSize)); 100 PetscCall(DMPlexGetSupport(dm, face, &supp)); 101 // Filter the support 102 for (PetscInt s = 0; s < suppSize; ++s) { 103 // Filter support 104 indC[Ns] = GetPointIndex(supp[s], cStart, cEnd, cells); 105 indS[Ns] = s; 106 if (indC[Ns] >= 0) ++Ns; 107 } 108 if (Ns < 2) PetscFunctionReturn(PETSC_SUCCESS); 109 PetscCheck(Ns == 2, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Faces should separate only two cells, not %" PetscInt_FMT, Ns); 110 PetscCheck(indC[0] >= 0 && indC[1] >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Support cells %" PetscInt_FMT " (%" PetscInt_FMT ") and %" PetscInt_FMT " (%" PetscInt_FMT ") are not both valid", supp[0], indC[0], supp[1], indC[1]); 111 seenA = PetscBTLookup(seenCells, indC[0]); 112 flippedA = PetscBTLookup(flippedCells, indC[0]) ? 1 : 0; 113 seenB = PetscBTLookup(seenCells, indC[1]); 114 flippedB = PetscBTLookup(flippedCells, indC[1]) ? 1 : 0; 115 116 PetscCall(DMPlexGetConeSize(dm, supp[indS[0]], &coneSizeA)); 117 PetscCall(DMPlexGetConeSize(dm, supp[indS[1]], &coneSizeB)); 118 PetscCall(DMPlexGetCone(dm, supp[indS[0]], &coneA)); 119 PetscCall(DMPlexGetCone(dm, supp[indS[1]], &coneB)); 120 PetscCall(DMPlexGetConeOrientation(dm, supp[indS[0]], &coneOA)); 121 PetscCall(DMPlexGetConeOrientation(dm, supp[indS[1]], &coneOB)); 122 for (PetscInt c = 0; c < coneSizeA; ++c) { 123 const PetscInt indF = GetPointIndex(coneA[c], fStart, fEnd, faces); 124 125 // Filter cone 126 if (indF < 0) continue; 127 if (!PetscBTLookup(seenFaces, indF)) { 128 faceFIFO[(*fBottom)++] = coneA[c]; 129 PetscCall(PetscBTSet(seenFaces, indF)); 130 } 131 if (coneA[c] == face) posA = c; 132 PetscCheck(*fBottom <= fEnd - fStart, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Face %" PetscInt_FMT " was pushed exceeding capacity %" PetscInt_FMT " > %" PetscInt_FMT, coneA[c], *fBottom, fEnd - fStart); 133 } 134 PetscCheck(posA >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Face %" PetscInt_FMT " could not be located in cell %" PetscInt_FMT, face, supp[indS[0]]); 135 for (PetscInt c = 0; c < coneSizeB; ++c) { 136 const PetscInt indF = GetPointIndex(coneB[c], fStart, fEnd, faces); 137 138 // Filter cone 139 if (indF < 0) continue; 140 if (!PetscBTLookup(seenFaces, indF)) { 141 faceFIFO[(*fBottom)++] = coneB[c]; 142 PetscCall(PetscBTSet(seenFaces, indF)); 143 } 144 if (coneB[c] == face) posB = c; 145 PetscCheck(*fBottom <= fEnd - fStart, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Face %" PetscInt_FMT " was pushed exceeding capacity %" PetscInt_FMT " > %" PetscInt_FMT, coneA[c], *fBottom, fEnd - fStart); 146 } 147 PetscCheck(posB >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Face %" PetscInt_FMT " could not be located in cell %" PetscInt_FMT, face, supp[indS[1]]); 148 149 if (dim == 1) { 150 mismatch = posA == posB; 151 } else { 152 mismatch = coneOA[posA] == coneOB[posB]; 153 } 154 155 if (mismatch ^ (flippedA ^ flippedB)) { 156 PetscCheck(!seenA || !seenB, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Previously seen cells %" PetscInt_FMT " and %" PetscInt_FMT " do not match: Fault mesh is non-orientable", supp[indS[0]], supp[indS[1]]); 157 if (!seenA && !flippedA) { 158 PetscCall(PetscBTSet(flippedCells, indC[0])); 159 } else if (!seenB && !flippedB) { 160 PetscCall(PetscBTSet(flippedCells, indC[1])); 161 } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Inconsistent mesh orientation: Fault mesh is non-orientable"); 162 } else PetscCheck(!mismatch || !flippedA || !flippedB, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Attempt to flip already flipped cell: Fault mesh is non-orientable"); 163 PetscCall(PetscBTSet(seenCells, indC[0])); 164 PetscCall(PetscBTSet(seenCells, indC[1])); 165 PetscFunctionReturn(PETSC_SUCCESS); 166 } 167 168 static PetscErrorCode DMPlexCheckFace_Old_Internal(DM dm, PetscInt *faceFIFO, PetscInt *fTop, PetscInt *fBottom, PetscInt cStart, PetscInt fStart, PetscInt fEnd, PetscBT seenCells, PetscBT flippedCells, PetscBT seenFaces) 169 { 170 const PetscInt *support, *coneA, *coneB, *coneOA, *coneOB; 171 PetscInt supportSize, coneSizeA, coneSizeB, posA = -1, posB = -1; 172 PetscInt face, dim, seenA, flippedA, seenB, flippedB, mismatch, c; 173 174 PetscFunctionBegin; 175 face = faceFIFO[(*fTop)++]; 176 PetscCall(DMGetDimension(dm, &dim)); 177 PetscCall(DMPlexGetSupportSize(dm, face, &supportSize)); 178 PetscCall(DMPlexGetSupport(dm, face, &support)); 179 if (supportSize < 2) PetscFunctionReturn(PETSC_SUCCESS); 180 PetscCheck(supportSize == 2, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Faces should separate only two cells, not %" PetscInt_FMT, supportSize); 181 seenA = PetscBTLookup(seenCells, support[0] - cStart); 182 flippedA = PetscBTLookup(flippedCells, support[0] - cStart) ? 1 : 0; 183 seenB = PetscBTLookup(seenCells, support[1] - cStart); 184 flippedB = PetscBTLookup(flippedCells, support[1] - cStart) ? 1 : 0; 185 186 PetscCall(DMPlexGetConeSize(dm, support[0], &coneSizeA)); 187 PetscCall(DMPlexGetConeSize(dm, support[1], &coneSizeB)); 188 PetscCall(DMPlexGetCone(dm, support[0], &coneA)); 189 PetscCall(DMPlexGetCone(dm, support[1], &coneB)); 190 PetscCall(DMPlexGetConeOrientation(dm, support[0], &coneOA)); 191 PetscCall(DMPlexGetConeOrientation(dm, support[1], &coneOB)); 192 for (c = 0; c < coneSizeA; ++c) { 193 if (!PetscBTLookup(seenFaces, coneA[c] - fStart)) { 194 faceFIFO[(*fBottom)++] = coneA[c]; 195 PetscCall(PetscBTSet(seenFaces, coneA[c] - fStart)); 196 } 197 if (coneA[c] == face) posA = c; 198 PetscCheck(*fBottom <= fEnd - fStart, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Face %" PetscInt_FMT " was pushed exceeding capacity %" PetscInt_FMT " > %" PetscInt_FMT, coneA[c], *fBottom, fEnd - fStart); 199 } 200 PetscCheck(posA >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Face %" PetscInt_FMT " could not be located in cell %" PetscInt_FMT, face, support[0]); 201 for (c = 0; c < coneSizeB; ++c) { 202 if (!PetscBTLookup(seenFaces, coneB[c] - fStart)) { 203 faceFIFO[(*fBottom)++] = coneB[c]; 204 PetscCall(PetscBTSet(seenFaces, coneB[c] - fStart)); 205 } 206 if (coneB[c] == face) posB = c; 207 PetscCheck(*fBottom <= fEnd - fStart, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Face %" PetscInt_FMT " was pushed exceeding capacity %" PetscInt_FMT " > %" PetscInt_FMT, coneA[c], *fBottom, fEnd - fStart); 208 } 209 PetscCheck(posB >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Face %" PetscInt_FMT " could not be located in cell %" PetscInt_FMT, face, support[1]); 210 211 if (dim == 1) { 212 mismatch = posA == posB; 213 } else { 214 mismatch = coneOA[posA] == coneOB[posB]; 215 } 216 217 if (mismatch ^ (flippedA ^ flippedB)) { 218 PetscCheck(!seenA || !seenB, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Previously seen cells %" PetscInt_FMT " and %" PetscInt_FMT " do not match: Fault mesh is non-orientable", support[0], support[1]); 219 if (!seenA && !flippedA) { 220 PetscCall(PetscBTSet(flippedCells, support[0] - cStart)); 221 } else if (!seenB && !flippedB) { 222 PetscCall(PetscBTSet(flippedCells, support[1] - cStart)); 223 } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Inconsistent mesh orientation: Fault mesh is non-orientable"); 224 } else PetscCheck(!mismatch || !flippedA || !flippedB, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Attempt to flip already flipped cell: Fault mesh is non-orientable"); 225 PetscCall(PetscBTSet(seenCells, support[0] - cStart)); 226 PetscCall(PetscBTSet(seenCells, support[1] - cStart)); 227 PetscFunctionReturn(PETSC_SUCCESS); 228 } 229 230 /* 231 DMPlexOrient_Serial - Compute valid orientation for local connected components 232 233 Not collective 234 235 Input Parameters: 236 + dm - The `DM` 237 - cellHeight - The height of k-cells to be oriented 238 239 Output Parameters: 240 + Ncomp - The number of connected component 241 . cellComp - The connected component for each local cell 242 . faceComp - The connected component for each local face 243 - flippedCells - Marked cells should be inverted 244 245 Level: developer 246 247 .seealso: `DMPlexOrient()` 248 */ 249 static PetscErrorCode DMPlexOrient_Serial(DM dm, IS cellIS, IS faceIS, PetscInt *Ncomp, PetscInt cellComp[], PetscInt faceComp[], PetscBT flippedCells) 250 { 251 PetscBT seenCells, seenFaces; 252 PetscInt *faceFIFO; 253 const PetscInt *cells = NULL, *faces = NULL; 254 PetscInt cStart = 0, cEnd = 0, fStart = 0, fEnd = 0; 255 256 PetscFunctionBegin; 257 if (cellIS) PetscCall(ISGetPointRange(cellIS, &cStart, &cEnd, &cells)); 258 if (faceIS) PetscCall(ISGetPointRange(faceIS, &fStart, &fEnd, &faces)); 259 PetscCall(PetscBTCreate(cEnd - cStart, &seenCells)); 260 PetscCall(PetscBTMemzero(cEnd - cStart, seenCells)); 261 PetscCall(PetscBTCreate(fEnd - fStart, &seenFaces)); 262 PetscCall(PetscBTMemzero(fEnd - fStart, seenFaces)); 263 PetscCall(PetscMalloc1(fEnd - fStart, &faceFIFO)); 264 *Ncomp = 0; 265 for (PetscInt c = 0; c < cEnd - cStart; ++c) cellComp[c] = -1; 266 do { 267 PetscInt cc, fTop, fBottom; 268 269 // Look for first unmarked cell 270 for (cc = cStart; cc < cEnd; ++cc) 271 if (cellComp[cc - cStart] < 0) break; 272 if (cc >= cEnd) break; 273 // Initialize FIFO with first cell in component 274 { 275 const PetscInt cell = cells ? cells[cc] : cc; 276 const PetscInt *cone; 277 PetscInt coneSize; 278 279 fTop = fBottom = 0; 280 PetscCall(DMPlexGetConeSize(dm, cell, &coneSize)); 281 PetscCall(DMPlexGetCone(dm, cell, &cone)); 282 for (PetscInt c = 0; c < coneSize; ++c) { 283 const PetscInt idx = GetPointIndex(cone[c], fStart, fEnd, faces); 284 285 // Cell faces are guaranteed to be in the face set 286 PetscCheck(idx >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Face %" PetscInt_FMT " of cell %" PetscInt_FMT " is not present in the label", cone[c], cell); 287 faceFIFO[fBottom++] = cone[c]; 288 PetscCall(PetscBTSet(seenFaces, idx)); 289 } 290 PetscCall(PetscBTSet(seenCells, cc - cStart)); 291 } 292 // Consider each face in FIFO 293 while (fTop < fBottom) PetscCall(DMPlexCheckFace_Internal(dm, faceFIFO, &fTop, &fBottom, cellIS, faceIS, seenCells, flippedCells, seenFaces)); 294 // Set component for cells and faces 295 for (PetscInt c = 0; c < cEnd - cStart; ++c) { 296 if (PetscBTLookup(seenCells, c)) cellComp[c] = *Ncomp; 297 } 298 for (PetscInt f = 0; f < fEnd - fStart; ++f) { 299 if (PetscBTLookup(seenFaces, f)) faceComp[f] = *Ncomp; 300 } 301 // Wipe seenCells and seenFaces for next component 302 PetscCall(PetscBTMemzero(fEnd - fStart, seenFaces)); 303 PetscCall(PetscBTMemzero(cEnd - cStart, seenCells)); 304 ++(*Ncomp); 305 } while (1); 306 PetscCall(PetscBTDestroy(&seenCells)); 307 PetscCall(PetscBTDestroy(&seenFaces)); 308 PetscCall(PetscFree(faceFIFO)); 309 PetscFunctionReturn(PETSC_SUCCESS); 310 } 311 312 /*@ 313 DMPlexOrient - Give a consistent orientation to the input mesh 314 315 Input Parameter: 316 . dm - The `DM` 317 318 Note: 319 The orientation data for the `DM` are change in-place. 320 321 This routine will fail for non-orientable surfaces, such as the Moebius strip. 322 323 Level: advanced 324 325 .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMCreate()` 326 @*/ 327 PetscErrorCode DMPlexOrient(DM dm) 328 { 329 #if 0 330 IS cellIS, faceIS; 331 332 PetscFunctionBegin; 333 PetscCall(DMPlexGetAllCells_Internal(dm, &cellIS)); 334 PetscCall(DMPlexGetAllFaces_Internal(dm, &faceIS)); 335 PetscCall(DMPlexOrientCells_Internal(dm, cellIS, faceIS)); 336 PetscCall(ISDestroy(&cellIS)); 337 PetscCall(ISDestroy(&faceIS)); 338 PetscFunctionReturn(PETSC_SUCCESS); 339 #else 340 MPI_Comm comm; 341 PetscSF sf; 342 const PetscInt *lpoints; 343 const PetscSFNode *rpoints; 344 PetscSFNode *rorntComp = NULL, *lorntComp = NULL; 345 PetscInt *numNeighbors, **neighbors, *locSupport = NULL; 346 PetscSFNode *nrankComp; 347 PetscBool *match, *flipped; 348 PetscBT seenCells, flippedCells, seenFaces; 349 PetscInt *faceFIFO, fTop, fBottom, *cellComp, *faceComp; 350 PetscInt numLeaves, numRoots, dim, h, cStart, cEnd, c, cell, fStart, fEnd, face, off, totNeighbors = 0; 351 PetscMPIInt rank, size, numComponents, comp = 0; 352 PetscBool flg, flg2; 353 PetscViewer viewer = NULL, selfviewer = NULL; 354 355 PetscFunctionBegin; 356 PetscCall(PetscObjectGetComm((PetscObject)dm, &comm)); 357 PetscCallMPI(MPI_Comm_rank(comm, &rank)); 358 PetscCallMPI(MPI_Comm_size(comm, &size)); 359 PetscCall(PetscOptionsHasName(((PetscObject)dm)->options, ((PetscObject)dm)->prefix, "-orientation_view", &flg)); 360 PetscCall(PetscOptionsHasName(((PetscObject)dm)->options, ((PetscObject)dm)->prefix, "-orientation_view_synchronized", &flg2)); 361 PetscCall(DMGetPointSF(dm, &sf)); 362 PetscCall(PetscSFGetGraph(sf, &numRoots, &numLeaves, &lpoints, &rpoints)); 363 /* Truth Table 364 mismatch flips do action mismatch flipA ^ flipB action 365 F 0 flips no F F F 366 F 1 flip yes F T T 367 F 2 flips no T F T 368 T 0 flips yes T T F 369 T 1 flip no 370 T 2 flips yes 371 */ 372 PetscCall(DMGetDimension(dm, &dim)); 373 PetscCall(DMPlexGetVTKCellHeight(dm, &h)); 374 PetscCall(DMPlexGetHeightStratum(dm, h, &cStart, &cEnd)); 375 PetscCall(DMPlexGetHeightStratum(dm, h + 1, &fStart, &fEnd)); 376 PetscCall(PetscBTCreate(cEnd - cStart, &seenCells)); 377 PetscCall(PetscBTMemzero(cEnd - cStart, seenCells)); 378 PetscCall(PetscBTCreate(cEnd - cStart, &flippedCells)); 379 PetscCall(PetscBTMemzero(cEnd - cStart, flippedCells)); 380 PetscCall(PetscBTCreate(fEnd - fStart, &seenFaces)); 381 PetscCall(PetscBTMemzero(fEnd - fStart, seenFaces)); 382 PetscCall(PetscCalloc3(fEnd - fStart, &faceFIFO, cEnd - cStart, &cellComp, fEnd - fStart, &faceComp)); 383 /* 384 OLD STYLE 385 - Add an integer array over cells and faces (component) for connected component number 386 Foreach component 387 - Mark the initial cell as seen 388 - Process component as usual 389 - Set component for all seenCells 390 - Wipe seenCells and seenFaces (flippedCells can stay) 391 - Generate parallel adjacency for component using SF and seenFaces 392 - Collect numComponents adj data from each proc to 0 393 - Build same serial graph 394 - Use same solver 395 - Use Scatterv to send back flipped flags for each component 396 - Negate flippedCells by component 397 398 NEW STYLE 399 - Create the adj on each process 400 - Bootstrap to complete graph on proc 0 401 */ 402 /* Loop over components */ 403 for (cell = cStart; cell < cEnd; ++cell) cellComp[cell - cStart] = -1; 404 do { 405 /* Look for first unmarked cell */ 406 for (cell = cStart; cell < cEnd; ++cell) 407 if (cellComp[cell - cStart] < 0) break; 408 if (cell >= cEnd) break; 409 /* Initialize FIFO with first cell in component */ 410 { 411 const PetscInt *cone; 412 PetscInt coneSize; 413 414 fTop = fBottom = 0; 415 PetscCall(DMPlexGetConeSize(dm, cell, &coneSize)); 416 PetscCall(DMPlexGetCone(dm, cell, &cone)); 417 for (c = 0; c < coneSize; ++c) { 418 faceFIFO[fBottom++] = cone[c]; 419 PetscCall(PetscBTSet(seenFaces, cone[c] - fStart)); 420 } 421 PetscCall(PetscBTSet(seenCells, cell - cStart)); 422 } 423 /* Consider each face in FIFO */ 424 while (fTop < fBottom) PetscCall(DMPlexCheckFace_Old_Internal(dm, faceFIFO, &fTop, &fBottom, cStart, fStart, fEnd, seenCells, flippedCells, seenFaces)); 425 /* Set component for cells and faces */ 426 for (cell = 0; cell < cEnd - cStart; ++cell) { 427 if (PetscBTLookup(seenCells, cell)) cellComp[cell] = comp; 428 } 429 for (face = 0; face < fEnd - fStart; ++face) { 430 if (PetscBTLookup(seenFaces, face)) faceComp[face] = comp; 431 } 432 /* Wipe seenCells and seenFaces for next component */ 433 PetscCall(PetscBTMemzero(fEnd - fStart, seenFaces)); 434 PetscCall(PetscBTMemzero(cEnd - cStart, seenCells)); 435 ++comp; 436 } while (1); 437 numComponents = comp; 438 if (flg) { 439 PetscViewer v; 440 441 PetscCall(PetscViewerASCIIGetStdout(comm, &v)); 442 PetscCall(PetscViewerASCIIPushSynchronized(v)); 443 PetscCall(PetscViewerASCIISynchronizedPrintf(v, "[%d]BT for serial flipped cells:\n", rank)); 444 PetscCall(PetscBTView(cEnd - cStart, flippedCells, v)); 445 PetscCall(PetscViewerFlush(v)); 446 PetscCall(PetscViewerASCIIPopSynchronized(v)); 447 } 448 /* Now all subdomains are oriented, but we need a consistent parallel orientation */ 449 if (numLeaves >= 0) { 450 PetscInt maxSupportSize, neighbor; 451 452 /* Store orientations of boundary faces*/ 453 PetscCall(DMPlexGetMaxSizes(dm, NULL, &maxSupportSize)); 454 PetscCall(PetscCalloc3(numRoots, &rorntComp, numRoots, &lorntComp, maxSupportSize, &locSupport)); 455 for (face = fStart; face < fEnd; ++face) { 456 const PetscInt *cone, *support, *ornt; 457 PetscInt coneSize, supportSize, Ns = 0, s, l; 458 459 PetscCall(DMPlexGetSupportSize(dm, face, &supportSize)); 460 /* Ignore overlapping cells */ 461 PetscCall(DMPlexGetSupport(dm, face, &support)); 462 for (s = 0; s < supportSize; ++s) { 463 if (lpoints) PetscCall(PetscFindInt(support[s], numLeaves, lpoints, &l)); 464 else { 465 if (support[s] >= 0 && support[s] < numLeaves) l = support[s]; 466 else l = -1; 467 } 468 if (l >= 0) continue; 469 locSupport[Ns++] = support[s]; 470 } 471 if (Ns != 1) continue; 472 neighbor = locSupport[0]; 473 PetscCall(DMPlexGetCone(dm, neighbor, &cone)); 474 PetscCall(DMPlexGetConeSize(dm, neighbor, &coneSize)); 475 PetscCall(DMPlexGetConeOrientation(dm, neighbor, &ornt)); 476 for (c = 0; c < coneSize; ++c) 477 if (cone[c] == face) break; 478 if (dim == 1) { 479 /* Use cone position instead, shifted to -1 or 1 */ 480 if (PetscBTLookup(flippedCells, neighbor - cStart)) rorntComp[face].rank = 1 - c * 2; 481 else rorntComp[face].rank = c * 2 - 1; 482 } else { 483 if (PetscBTLookup(flippedCells, neighbor - cStart)) rorntComp[face].rank = ornt[c] < 0 ? -1 : 1; 484 else rorntComp[face].rank = ornt[c] < 0 ? 1 : -1; 485 } 486 rorntComp[face].index = faceComp[face - fStart]; 487 } 488 /* Communicate boundary edge orientations */ 489 PetscCall(PetscSFBcastBegin(sf, MPIU_SF_NODE, rorntComp, lorntComp, MPI_REPLACE)); 490 PetscCall(PetscSFBcastEnd(sf, MPIU_SF_NODE, rorntComp, lorntComp, MPI_REPLACE)); 491 } 492 /* Get process adjacency */ 493 PetscCall(PetscMalloc2(numComponents, &numNeighbors, numComponents, &neighbors)); 494 viewer = PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)dm)); 495 if (flg2) PetscCall(PetscViewerASCIIPushSynchronized(viewer)); 496 PetscCall(PetscViewerGetSubViewer(viewer, PETSC_COMM_SELF, &selfviewer)); 497 for (comp = 0; comp < numComponents; ++comp) { 498 PetscInt l, n; 499 500 numNeighbors[comp] = 0; 501 PetscCall(PetscMalloc1(PetscMax(numLeaves, 0), &neighbors[comp])); 502 /* I know this is p^2 time in general, but for bounded degree its alright */ 503 for (l = 0; l < numLeaves; ++l) { 504 const PetscInt face = lpoints ? lpoints[l] : l; 505 506 /* Find a representative face (edge) separating pairs of procs */ 507 if ((face >= fStart) && (face < fEnd) && (faceComp[face - fStart] == comp) && rorntComp[face].rank) { 508 const PetscInt rrank = rpoints[l].rank; 509 const PetscInt rcomp = lorntComp[face].index; 510 511 for (n = 0; n < numNeighbors[comp]; ++n) 512 if ((rrank == rpoints[neighbors[comp][n]].rank) && (rcomp == lorntComp[lpoints[neighbors[comp][n]]].index)) break; 513 if (n >= numNeighbors[comp]) { 514 PetscInt supportSize; 515 516 PetscCall(DMPlexGetSupportSize(dm, face, &supportSize)); 517 // We can have internal faces in the SF if we have cells in the SF 518 if (supportSize > 1) continue; 519 if (flg) 520 PetscCall(PetscViewerASCIIPrintf(selfviewer, "[%d]: component %d, Found representative leaf %" PetscInt_FMT " (face %" PetscInt_FMT ") connecting to face %" PetscInt_FMT " on (%" PetscInt_FMT ", %" PetscInt_FMT ") with orientation %" PetscInt_FMT "\n", rank, comp, l, face, 521 rpoints[l].index, rrank, rcomp, lorntComp[face].rank)); 522 neighbors[comp][numNeighbors[comp]++] = l; 523 } 524 } 525 } 526 totNeighbors += numNeighbors[comp]; 527 } 528 PetscCall(PetscViewerRestoreSubViewer(viewer, PETSC_COMM_SELF, &selfviewer)); 529 if (flg2) PetscCall(PetscViewerASCIIPopSynchronized(viewer)); 530 PetscCall(PetscMalloc2(totNeighbors, &nrankComp, totNeighbors, &match)); 531 for (comp = 0, off = 0; comp < numComponents; ++comp) { 532 PetscInt n; 533 534 for (n = 0; n < numNeighbors[comp]; ++n, ++off) { 535 const PetscInt face = lpoints ? lpoints[neighbors[comp][n]] : neighbors[comp][n]; 536 const PetscInt o = rorntComp[face].rank * lorntComp[face].rank; 537 538 if (o < 0) match[off] = PETSC_TRUE; 539 else if (o > 0) match[off] = PETSC_FALSE; 540 else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid face %" PetscInt_FMT " (%" PetscInt_FMT ", %" PetscInt_FMT ") neighbor: %" PetscInt_FMT " comp: %d", face, rorntComp[face].rank, lorntComp[face].rank, neighbors[comp][n], comp); 541 nrankComp[off].rank = rpoints[neighbors[comp][n]].rank; 542 nrankComp[off].index = lorntComp[lpoints ? lpoints[neighbors[comp][n]] : neighbors[comp][n]].index; 543 } 544 PetscCall(PetscFree(neighbors[comp])); 545 } 546 /* Collect the graph on 0 */ 547 if (numLeaves >= 0) { 548 Mat G; 549 PetscBT seenProcs, flippedProcs; 550 PetscInt *procFIFO, pTop, pBottom; 551 PetscInt *N = NULL, *Noff; 552 PetscSFNode *adj = NULL; 553 PetscBool *val = NULL; 554 PetscMPIInt *recvcounts = NULL, *displs = NULL, *Nc, p, o, itotNeighbors; 555 PetscMPIInt size = 0; 556 557 PetscCall(PetscCalloc1(numComponents, &flipped)); 558 if (rank == 0) PetscCallMPI(MPI_Comm_size(comm, &size)); 559 PetscCall(PetscCalloc4(size, &recvcounts, size + 1, &displs, size, &Nc, size + 1, &Noff)); 560 PetscCallMPI(MPI_Gather(&numComponents, 1, MPI_INT, Nc, 1, MPI_INT, 0, comm)); 561 for (p = 0; p < size; ++p) displs[p + 1] = displs[p] + Nc[p]; 562 if (rank == 0) PetscCall(PetscMalloc1(displs[size], &N)); 563 PetscCallMPI(MPI_Gatherv(numNeighbors, numComponents, MPIU_INT, N, Nc, displs, MPIU_INT, 0, comm)); 564 for (p = 0, o = 0; p < size; ++p) { 565 recvcounts[p] = 0; 566 for (c = 0; c < Nc[p]; ++c, ++o) recvcounts[p] += N[o]; 567 displs[p + 1] = displs[p] + recvcounts[p]; 568 } 569 if (rank == 0) PetscCall(PetscMalloc2(displs[size], &adj, displs[size], &val)); 570 PetscCall(PetscMPIIntCast(totNeighbors, &itotNeighbors)); 571 PetscCallMPI(MPI_Gatherv(nrankComp, itotNeighbors, MPIU_SF_NODE, adj, recvcounts, displs, MPIU_SF_NODE, 0, comm)); 572 PetscCallMPI(MPI_Gatherv(match, itotNeighbors, MPIU_BOOL, val, recvcounts, displs, MPIU_BOOL, 0, comm)); 573 PetscCall(PetscFree2(numNeighbors, neighbors)); 574 if (rank == 0) { 575 for (p = 1; p <= size; ++p) Noff[p] = Noff[p - 1] + Nc[p - 1]; 576 if (flg) { 577 PetscInt n; 578 579 for (p = 0, off = 0; p < size; ++p) { 580 for (c = 0; c < Nc[p]; ++c) { 581 PetscCall(PetscPrintf(PETSC_COMM_SELF, "Proc %d Comp %" PetscInt_FMT ":\n", p, c)); 582 for (n = 0; n < N[Noff[p] + c]; ++n, ++off) PetscCall(PetscPrintf(PETSC_COMM_SELF, " edge (%" PetscInt_FMT ", %" PetscInt_FMT ") (%s):\n", adj[off].rank, adj[off].index, PetscBools[val[off]])); 583 } 584 } 585 } 586 /* Symmetrize the graph */ 587 PetscCall(MatCreate(PETSC_COMM_SELF, &G)); 588 PetscCall(MatSetSizes(G, Noff[size], Noff[size], Noff[size], Noff[size])); 589 PetscCall(MatSetUp(G)); 590 for (p = 0, off = 0; p < size; ++p) { 591 for (c = 0; c < Nc[p]; ++c) { 592 const PetscInt r = Noff[p] + c; 593 PetscInt n; 594 595 for (n = 0; n < N[r]; ++n, ++off) { 596 const PetscInt q = Noff[adj[off].rank] + adj[off].index; 597 const PetscScalar o = val[off] ? 1.0 : 0.0; 598 599 PetscCall(MatSetValues(G, 1, &r, 1, &q, &o, INSERT_VALUES)); 600 PetscCall(MatSetValues(G, 1, &q, 1, &r, &o, INSERT_VALUES)); 601 } 602 } 603 } 604 PetscCall(MatAssemblyBegin(G, MAT_FINAL_ASSEMBLY)); 605 PetscCall(MatAssemblyEnd(G, MAT_FINAL_ASSEMBLY)); 606 607 PetscCall(PetscBTCreate(Noff[size], &seenProcs)); 608 PetscCall(PetscBTMemzero(Noff[size], seenProcs)); 609 PetscCall(PetscBTCreate(Noff[size], &flippedProcs)); 610 PetscCall(PetscBTMemzero(Noff[size], flippedProcs)); 611 PetscCall(PetscMalloc1(Noff[size], &procFIFO)); 612 pTop = pBottom = 0; 613 for (p = 0; p < Noff[size]; ++p) { 614 if (PetscBTLookup(seenProcs, p)) continue; 615 /* Initialize FIFO with next proc */ 616 procFIFO[pBottom++] = p; 617 PetscCall(PetscBTSet(seenProcs, p)); 618 /* Consider each proc in FIFO */ 619 while (pTop < pBottom) { 620 const PetscScalar *ornt; 621 const PetscInt *neighbors; 622 PetscInt proc, nproc, seen, flippedA, flippedB, mismatch, numNeighbors, n; 623 624 proc = procFIFO[pTop++]; 625 flippedA = PetscBTLookup(flippedProcs, proc) ? 1 : 0; 626 PetscCall(MatGetRow(G, proc, &numNeighbors, &neighbors, &ornt)); 627 /* Loop over neighboring procs */ 628 for (n = 0; n < numNeighbors; ++n) { 629 nproc = neighbors[n]; 630 mismatch = PetscRealPart(ornt[n]) > 0.5 ? 0 : 1; 631 seen = PetscBTLookup(seenProcs, nproc); 632 flippedB = PetscBTLookup(flippedProcs, nproc) ? 1 : 0; 633 634 if (mismatch ^ (flippedA ^ flippedB)) { 635 PetscCheck(!seen, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Previously seen procs %" PetscInt_FMT " and %" PetscInt_FMT " do not match: Fault mesh is non-orientable", proc, nproc); 636 if (!flippedB) { 637 PetscCall(PetscBTSet(flippedProcs, nproc)); 638 } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Inconsistent mesh orientation: Fault mesh is non-orientable"); 639 } else PetscCheck(!mismatch || !flippedA || !flippedB, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Attempt to flip already flipped cell: Fault mesh is non-orientable"); 640 if (!seen) { 641 procFIFO[pBottom++] = nproc; 642 PetscCall(PetscBTSet(seenProcs, nproc)); 643 } 644 } 645 } 646 } 647 PetscCall(PetscFree(procFIFO)); 648 PetscCall(MatDestroy(&G)); 649 PetscCall(PetscFree2(adj, val)); 650 PetscCall(PetscBTDestroy(&seenProcs)); 651 } 652 /* Scatter flip flags */ 653 { 654 PetscBool *flips = NULL; 655 656 if (rank == 0) { 657 PetscCall(PetscMalloc1(Noff[size], &flips)); 658 for (p = 0; p < Noff[size]; ++p) { 659 flips[p] = PetscBTLookup(flippedProcs, p) ? PETSC_TRUE : PETSC_FALSE; 660 if (flg && flips[p]) PetscCall(PetscPrintf(comm, "Flipping Proc+Comp %d:\n", p)); 661 } 662 for (p = 0; p < size; ++p) displs[p + 1] = displs[p] + Nc[p]; 663 } 664 PetscCallMPI(MPI_Scatterv(flips, Nc, displs, MPIU_BOOL, flipped, numComponents, MPIU_BOOL, 0, comm)); 665 PetscCall(PetscFree(flips)); 666 } 667 if (rank == 0) PetscCall(PetscBTDestroy(&flippedProcs)); 668 PetscCall(PetscFree(N)); 669 PetscCall(PetscFree4(recvcounts, displs, Nc, Noff)); 670 PetscCall(PetscFree2(nrankComp, match)); 671 672 /* Decide whether to flip cells in each component */ 673 for (c = 0; c < cEnd - cStart; ++c) { 674 if (flipped[cellComp[c]]) PetscCall(PetscBTNegate(flippedCells, c)); 675 } 676 PetscCall(PetscFree(flipped)); 677 } 678 if (flg) { 679 PetscViewer v; 680 681 PetscCall(PetscViewerASCIIGetStdout(comm, &v)); 682 PetscCall(PetscViewerASCIIPushSynchronized(v)); 683 PetscCall(PetscViewerASCIISynchronizedPrintf(v, "[%d]BT for parallel flipped cells:\n", rank)); 684 PetscCall(PetscBTView(cEnd - cStart, flippedCells, v)); 685 PetscCall(PetscViewerFlush(v)); 686 PetscCall(PetscViewerASCIIPopSynchronized(v)); 687 } 688 /* Reverse flipped cells in the mesh */ 689 for (c = cStart; c < cEnd; ++c) { 690 if (PetscBTLookup(flippedCells, c - cStart)) PetscCall(DMPlexOrientPoint(dm, c, -1)); 691 } 692 PetscCall(PetscBTDestroy(&seenCells)); 693 PetscCall(PetscBTDestroy(&flippedCells)); 694 PetscCall(PetscBTDestroy(&seenFaces)); 695 PetscCall(PetscFree2(numNeighbors, neighbors)); 696 PetscCall(PetscFree3(rorntComp, lorntComp, locSupport)); 697 PetscCall(PetscFree3(faceFIFO, cellComp, faceComp)); 698 PetscFunctionReturn(PETSC_SUCCESS); 699 #endif 700 } 701 702 static PetscErrorCode CreateCellAndFaceIS_Private(DM dm, DMLabel label, IS *cellIS, IS *faceIS) 703 { 704 IS valueIS; 705 const PetscInt *values; 706 PetscInt Nv, depth = 0; 707 708 PetscFunctionBegin; 709 PetscCall(DMLabelGetValueIS(label, &valueIS)); 710 PetscCall(ISGetLocalSize(valueIS, &Nv)); 711 PetscCall(ISGetIndices(valueIS, &values)); 712 for (PetscInt v = 0; v < Nv; ++v) { 713 const PetscInt val = values[v] < 0 || values[v] >= 100 ? 0 : values[v]; 714 PetscInt n; 715 716 PetscCall(DMLabelGetStratumSize(label, val, &n)); 717 if (!n) continue; 718 depth = PetscMax(val, depth); 719 } 720 PetscCall(ISDestroy(&valueIS)); 721 PetscCheck(depth >= 1 || !Nv, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Depth for interface must be at least 1, not %" PetscInt_FMT, depth); 722 PetscCall(DMLabelGetStratumIS(label, depth, cellIS)); 723 PetscCall(DMLabelGetStratumIS(label, depth - 1, faceIS)); 724 if (!*cellIS) PetscCall(ISCreateStride(PETSC_COMM_SELF, 0, 0, 1, cellIS)); 725 if (!*faceIS) PetscCall(ISCreateStride(PETSC_COMM_SELF, 0, 0, 1, faceIS)); 726 PetscFunctionReturn(PETSC_SUCCESS); 727 } 728 729 PetscErrorCode DMPlexOrientLabel(DM dm, DMLabel label) 730 { 731 IS cellIS, faceIS; 732 733 PetscFunctionBegin; 734 PetscCall(CreateCellAndFaceIS_Private(dm, label, &cellIS, &faceIS)); 735 PetscCall(DMPlexOrientCells_Internal(dm, cellIS, faceIS)); 736 PetscCall(ISDestroy(&cellIS)); 737 PetscCall(ISDestroy(&faceIS)); 738 PetscFunctionReturn(PETSC_SUCCESS); 739 } 740 741 PetscErrorCode DMPlexOrientCells_Internal(DM dm, IS cellIS, IS faceIS) 742 { 743 MPI_Comm comm; 744 PetscSF sf; 745 const PetscInt *lpoints; 746 const PetscSFNode *rpoints; 747 PetscSFNode *rorntComp = NULL, *lorntComp = NULL; 748 PetscInt *numNeighbors, **neighbors, *locSupp = NULL; 749 PetscSFNode *nrankComp; 750 PetscBool *match, *flipped; 751 PetscBT flippedCells; 752 PetscInt *cellComp, *faceComp; 753 const PetscInt *cells = NULL, *faces = NULL; 754 PetscInt cStart = 0, cEnd = 0, fStart = 0, fEnd = 0; 755 PetscInt numLeaves, numRoots, dim, Ncomp, totNeighbors = 0; 756 PetscMPIInt rank, size; 757 PetscBool view, viewSync; 758 PetscViewer viewer = NULL, selfviewer = NULL; 759 760 PetscFunctionBegin; 761 PetscCall(PetscObjectGetComm((PetscObject)dm, &comm)); 762 PetscCallMPI(MPI_Comm_rank(comm, &rank)); 763 PetscCallMPI(MPI_Comm_size(comm, &size)); 764 PetscCall(PetscOptionsHasName(((PetscObject)dm)->options, ((PetscObject)dm)->prefix, "-orientation_view", &view)); 765 PetscCall(PetscOptionsHasName(((PetscObject)dm)->options, ((PetscObject)dm)->prefix, "-orientation_view_synchronized", &viewSync)); 766 767 if (cellIS) PetscCall(ISGetPointRange(cellIS, &cStart, &cEnd, &cells)); 768 if (faceIS) PetscCall(ISGetPointRange(faceIS, &fStart, &fEnd, &faces)); 769 PetscCall(DMGetPointSF(dm, &sf)); 770 PetscCall(PetscSFGetGraph(sf, &numRoots, &numLeaves, &lpoints, &rpoints)); 771 /* Truth Table 772 mismatch flips do action mismatch flipA ^ flipB action 773 F 0 flips no F F F 774 F 1 flip yes F T T 775 F 2 flips no T F T 776 T 0 flips yes T T F 777 T 1 flip no 778 T 2 flips yes 779 */ 780 PetscCall(DMGetDimension(dm, &dim)); 781 PetscCall(PetscBTCreate(cEnd - cStart, &flippedCells)); 782 PetscCall(PetscBTMemzero(cEnd - cStart, flippedCells)); 783 PetscCall(PetscCalloc2(cEnd - cStart, &cellComp, fEnd - fStart, &faceComp)); 784 /* 785 OLD STYLE 786 - Add an integer array over cells and faces (component) for connected component number 787 Foreach component 788 - Mark the initial cell as seen 789 - Process component as usual 790 - Set component for all seenCells 791 - Wipe seenCells and seenFaces (flippedCells can stay) 792 - Generate parallel adjacency for component using SF and seenFaces 793 - Collect Ncomp adj data from each proc to 0 794 - Build same serial graph 795 - Use same solver 796 - Use Scatterv to send back flipped flags for each component 797 - Negate flippedCells by component 798 799 NEW STYLE 800 - Create the adj on each process 801 - Bootstrap to complete graph on proc 0 802 */ 803 PetscCall(DMPlexOrient_Serial(dm, cellIS, faceIS, &Ncomp, cellComp, faceComp, flippedCells)); 804 if (view) { 805 PetscViewer v; 806 PetscInt cdepth = -1; 807 808 PetscCall(PetscViewerASCIIGetStdout(comm, &v)); 809 PetscCall(PetscViewerASCIIPushSynchronized(v)); 810 if (cEnd > cStart) PetscCall(DMPlexGetPointDepth(dm, cells ? cells[cStart] : cStart, &cdepth)); 811 PetscCall(PetscViewerASCIISynchronizedPrintf(v, "[%d]New Orientation %" PetscInt_FMT " cells (depth %" PetscInt_FMT ") and %" PetscInt_FMT " faces\n", rank, cEnd - cStart, cdepth, fEnd - fStart)); 812 PetscCall(PetscViewerASCIISynchronizedPrintf(v, "[%d]BT for serial flipped cells:\n", rank)); 813 PetscCall(PetscBTView(cEnd - cStart, flippedCells, v)); 814 PetscCall(PetscViewerFlush(v)); 815 PetscCall(PetscViewerASCIIPopSynchronized(v)); 816 } 817 /* Now all subdomains are oriented, but we need a consistent parallel orientation */ 818 // TODO: This all has to be rewritten to filter cones/supports to the ISes 819 if (numLeaves >= 0) { 820 PetscInt maxSuppSize, neighbor; 821 822 // Store orientations of boundary faces 823 PetscCall(DMPlexGetMaxSizes(dm, NULL, &maxSuppSize)); 824 PetscCall(PetscCalloc3(numRoots, &rorntComp, numRoots, &lorntComp, maxSuppSize, &locSupp)); 825 for (PetscInt f = fStart; f < fEnd; ++f) { 826 const PetscInt face = faces ? faces[f] : f; 827 const PetscInt *cone, *supp, *ornt; 828 PetscInt coneSize, suppSize, nind, c, Ns = 0; 829 830 PetscCall(DMPlexGetSupportSize(dm, face, &suppSize)); 831 PetscCall(DMPlexGetSupport(dm, face, &supp)); 832 for (PetscInt s = 0; s < suppSize; ++s) { 833 PetscInt ind, l; 834 835 // Filter support 836 ind = GetPointIndex(supp[s], cStart, cEnd, cells); 837 if (ind < 0) continue; 838 // Ignore overlapping cells 839 PetscCall(PetscFindInt(supp[s], numLeaves, lpoints, &l)); 840 if (l >= 0) continue; 841 locSupp[Ns++] = supp[s]; 842 } 843 PetscCheck(Ns < maxSuppSize, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Index %" PetscInt_FMT " exceeds array size %" PetscInt_FMT, Ns, maxSuppSize); 844 if (Ns != 1) continue; 845 neighbor = locSupp[0]; 846 nind = GetPointIndex(neighbor, cStart, cEnd, cells); 847 PetscCall(DMPlexGetCone(dm, neighbor, &cone)); 848 PetscCall(DMPlexGetConeSize(dm, neighbor, &coneSize)); 849 PetscCall(DMPlexGetConeOrientation(dm, neighbor, &ornt)); 850 for (c = 0; c < coneSize; ++c) 851 if (cone[c] == face) break; 852 if (dim == 1) { 853 /* Use cone position instead, shifted to -1 or 1 */ 854 if (PetscBTLookup(flippedCells, nind)) rorntComp[face].rank = 1 - c * 2; 855 else rorntComp[face].rank = c * 2 - 1; 856 } else { 857 if (PetscBTLookup(flippedCells, nind)) rorntComp[face].rank = ornt[c] < 0 ? -1 : 1; 858 else rorntComp[face].rank = ornt[c] < 0 ? 1 : -1; 859 } 860 rorntComp[face].index = faceComp[GetPointIndex(face, fStart, fEnd, faces)]; 861 } 862 // Communicate boundary edge orientations 863 PetscCall(PetscSFBcastBegin(sf, MPIU_SF_NODE, rorntComp, lorntComp, MPI_REPLACE)); 864 PetscCall(PetscSFBcastEnd(sf, MPIU_SF_NODE, rorntComp, lorntComp, MPI_REPLACE)); 865 } 866 /* Get process adjacency */ 867 PetscCall(PetscMalloc2(Ncomp, &numNeighbors, Ncomp, &neighbors)); 868 viewer = PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)dm)); 869 if (viewSync) PetscCall(PetscViewerASCIIPushSynchronized(viewer)); 870 PetscCall(PetscViewerGetSubViewer(viewer, PETSC_COMM_SELF, &selfviewer)); 871 for (PetscInt comp = 0; comp < Ncomp; ++comp) { 872 PetscInt n; 873 874 numNeighbors[comp] = 0; 875 PetscCall(PetscMalloc1(PetscMax(numLeaves, 0), &neighbors[comp])); 876 /* I know this is p^2 time in general, but for bounded degree its alright */ 877 for (PetscInt l = 0; l < numLeaves; ++l) { 878 const PetscInt face = lpoints[l]; 879 PetscInt find; 880 881 /* Find a representative face (edge) separating pairs of procs */ 882 find = GetPointIndex(face, fStart, fEnd, faces); 883 if ((find >= 0) && (faceComp[find] == comp) && rorntComp[face].rank) { 884 const PetscInt rrank = rpoints[l].rank; 885 const PetscInt rcomp = lorntComp[face].index; 886 887 for (n = 0; n < numNeighbors[comp]; ++n) 888 if ((rrank == rpoints[neighbors[comp][n]].rank) && (rcomp == lorntComp[lpoints[neighbors[comp][n]]].index)) break; 889 if (n >= numNeighbors[comp]) { 890 const PetscInt *supp; 891 PetscInt suppSize, Ns = 0; 892 893 PetscCall(DMPlexGetSupport(dm, face, &supp)); 894 PetscCall(DMPlexGetSupportSize(dm, face, &suppSize)); 895 for (PetscInt s = 0; s < suppSize; ++s) { 896 // Filter support 897 if (GetPointIndex(supp[s], cStart, cEnd, cells) >= 0) ++Ns; 898 } 899 PetscCheck(Ns == 1, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Boundary face %" PetscInt_FMT " should see one cell, not %" PetscInt_FMT, face, Ns); 900 if (view) 901 PetscCall(PetscViewerASCIIPrintf(selfviewer, "[%d]: component %" PetscInt_FMT ", Found representative leaf %" PetscInt_FMT " (face %" PetscInt_FMT ") connecting to face %" PetscInt_FMT " on (%" PetscInt_FMT ", %" PetscInt_FMT ") with orientation %" PetscInt_FMT "\n", rank, comp, l, face, 902 rpoints[l].index, rrank, rcomp, lorntComp[face].rank)); 903 neighbors[comp][numNeighbors[comp]++] = l; 904 } 905 } 906 } 907 totNeighbors += numNeighbors[comp]; 908 } 909 PetscCall(PetscViewerRestoreSubViewer(viewer, PETSC_COMM_SELF, &selfviewer)); 910 if (viewSync) PetscCall(PetscViewerASCIIPopSynchronized(viewer)); 911 PetscCall(PetscMalloc2(totNeighbors, &nrankComp, totNeighbors, &match)); 912 for (PetscInt comp = 0, off = 0; comp < Ncomp; ++comp) { 913 for (PetscInt n = 0; n < numNeighbors[comp]; ++n, ++off) { 914 const PetscInt face = lpoints[neighbors[comp][n]]; 915 const PetscInt o = rorntComp[face].rank * lorntComp[face].rank; 916 917 if (o < 0) match[off] = PETSC_TRUE; 918 else if (o > 0) match[off] = PETSC_FALSE; 919 else 920 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid face %" PetscInt_FMT " (%" PetscInt_FMT ", %" PetscInt_FMT ") neighbor: %" PetscInt_FMT " comp: %" PetscInt_FMT, face, rorntComp[face].rank, lorntComp[face].rank, neighbors[comp][n], comp); 921 nrankComp[off].rank = rpoints[neighbors[comp][n]].rank; 922 nrankComp[off].index = lorntComp[lpoints[neighbors[comp][n]]].index; 923 } 924 PetscCall(PetscFree(neighbors[comp])); 925 } 926 /* Collect the graph on 0 */ 927 if (numLeaves >= 0) { 928 Mat G; 929 PetscBT seenProcs, flippedProcs; 930 PetscInt *procFIFO, pTop, pBottom; 931 PetscInt *N = NULL, *Noff; 932 PetscSFNode *adj = NULL; 933 PetscBool *val = NULL; 934 PetscMPIInt *recvcounts = NULL, *displs = NULL, *Nc; 935 PetscMPIInt size = 0, iNcomp, itotNeighbors; 936 937 PetscCall(PetscCalloc1(Ncomp, &flipped)); 938 if (rank == 0) PetscCallMPI(MPI_Comm_size(comm, &size)); 939 PetscCall(PetscCalloc4(size, &recvcounts, size + 1, &displs, size, &Nc, size + 1, &Noff)); 940 PetscCallMPI(MPI_Gather(&Ncomp, 1, MPI_INT, Nc, 1, MPI_INT, 0, comm)); 941 for (PetscInt p = 0; p < size; ++p) displs[p + 1] = displs[p] + Nc[p]; 942 if (rank == 0) PetscCall(PetscMalloc1(displs[size], &N)); 943 PetscCall(PetscMPIIntCast(Ncomp, &iNcomp)); 944 PetscCallMPI(MPI_Gatherv(numNeighbors, iNcomp, MPIU_INT, N, Nc, displs, MPIU_INT, 0, comm)); 945 for (PetscInt p = 0, o = 0; p < size; ++p) { 946 recvcounts[p] = 0; 947 for (PetscInt c = 0; c < Nc[p]; ++c, ++o) recvcounts[p] += N[o]; 948 displs[p + 1] = displs[p] + recvcounts[p]; 949 } 950 if (rank == 0) PetscCall(PetscMalloc2(displs[size], &adj, displs[size], &val)); 951 PetscCall(PetscMPIIntCast(totNeighbors, &itotNeighbors)); 952 PetscCallMPI(MPI_Gatherv(nrankComp, itotNeighbors, MPIU_SF_NODE, adj, recvcounts, displs, MPIU_SF_NODE, 0, comm)); 953 PetscCallMPI(MPI_Gatherv(match, itotNeighbors, MPIU_BOOL, val, recvcounts, displs, MPIU_BOOL, 0, comm)); 954 PetscCall(PetscFree2(numNeighbors, neighbors)); 955 if (rank == 0) { 956 for (PetscInt p = 1; p <= size; ++p) Noff[p] = Noff[p - 1] + Nc[p - 1]; 957 if (view) { 958 for (PetscInt p = 0, off = 0; p < size; ++p) { 959 for (PetscInt c = 0; c < Nc[p]; ++c) { 960 PetscCall(PetscPrintf(PETSC_COMM_SELF, "Proc %" PetscInt_FMT " Comp %" PetscInt_FMT ":\n", p, c)); 961 for (PetscInt n = 0; n < N[Noff[p] + c]; ++n, ++off) PetscCall(PetscPrintf(PETSC_COMM_SELF, " edge (%" PetscInt_FMT ", %" PetscInt_FMT ") (%s):\n", adj[off].rank, adj[off].index, PetscBools[val[off]])); 962 } 963 } 964 } 965 /* Symmetrize the graph */ 966 PetscCall(MatCreate(PETSC_COMM_SELF, &G)); 967 PetscCall(MatSetSizes(G, Noff[size], Noff[size], Noff[size], Noff[size])); 968 PetscCall(MatSetUp(G)); 969 for (PetscInt p = 0, off = 0; p < size; ++p) { 970 for (PetscInt c = 0; c < Nc[p]; ++c) { 971 const PetscInt r = Noff[p] + c; 972 973 for (PetscInt n = 0; n < N[r]; ++n, ++off) { 974 const PetscInt q = Noff[adj[off].rank] + adj[off].index; 975 const PetscScalar o = val[off] ? 1.0 : 0.0; 976 977 PetscCall(MatSetValues(G, 1, &r, 1, &q, &o, INSERT_VALUES)); 978 PetscCall(MatSetValues(G, 1, &q, 1, &r, &o, INSERT_VALUES)); 979 } 980 } 981 } 982 PetscCall(MatAssemblyBegin(G, MAT_FINAL_ASSEMBLY)); 983 PetscCall(MatAssemblyEnd(G, MAT_FINAL_ASSEMBLY)); 984 985 PetscCall(PetscBTCreate(Noff[size], &seenProcs)); 986 PetscCall(PetscBTMemzero(Noff[size], seenProcs)); 987 PetscCall(PetscBTCreate(Noff[size], &flippedProcs)); 988 PetscCall(PetscBTMemzero(Noff[size], flippedProcs)); 989 PetscCall(PetscMalloc1(Noff[size], &procFIFO)); 990 pTop = pBottom = 0; 991 for (PetscInt p = 0; p < Noff[size]; ++p) { 992 if (PetscBTLookup(seenProcs, p)) continue; 993 /* Initialize FIFO with next proc */ 994 procFIFO[pBottom++] = p; 995 PetscCall(PetscBTSet(seenProcs, p)); 996 /* Consider each proc in FIFO */ 997 while (pTop < pBottom) { 998 const PetscScalar *ornt; 999 const PetscInt *neighbors; 1000 PetscInt proc, nproc, seen, flippedA, flippedB, mismatch, numNeighbors; 1001 1002 proc = procFIFO[pTop++]; 1003 flippedA = PetscBTLookup(flippedProcs, proc) ? 1 : 0; 1004 PetscCall(MatGetRow(G, proc, &numNeighbors, &neighbors, &ornt)); 1005 /* Loop over neighboring procs */ 1006 for (PetscInt n = 0; n < numNeighbors; ++n) { 1007 nproc = neighbors[n]; 1008 mismatch = PetscRealPart(ornt[n]) > 0.5 ? 0 : 1; 1009 seen = PetscBTLookup(seenProcs, nproc); 1010 flippedB = PetscBTLookup(flippedProcs, nproc) ? 1 : 0; 1011 1012 if (mismatch ^ (flippedA ^ flippedB)) { 1013 PetscCheck(!seen, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Previously seen procs %" PetscInt_FMT " and %" PetscInt_FMT " do not match: Fault mesh is non-orientable", proc, nproc); 1014 if (!flippedB) { 1015 PetscCall(PetscBTSet(flippedProcs, nproc)); 1016 } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Inconsistent mesh orientation: Fault mesh is non-orientable"); 1017 } else PetscCheck(!mismatch || !flippedA || !flippedB, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Attempt to flip already flipped cell: Fault mesh is non-orientable"); 1018 if (!seen) { 1019 procFIFO[pBottom++] = nproc; 1020 PetscCall(PetscBTSet(seenProcs, nproc)); 1021 } 1022 } 1023 } 1024 } 1025 PetscCall(PetscFree(procFIFO)); 1026 PetscCall(MatDestroy(&G)); 1027 PetscCall(PetscFree2(adj, val)); 1028 PetscCall(PetscBTDestroy(&seenProcs)); 1029 } 1030 /* Scatter flip flags */ 1031 { 1032 PetscBool *flips = NULL; 1033 1034 if (rank == 0) { 1035 PetscCall(PetscMalloc1(Noff[size], &flips)); 1036 for (PetscInt p = 0; p < Noff[size]; ++p) { 1037 flips[p] = PetscBTLookup(flippedProcs, p) ? PETSC_TRUE : PETSC_FALSE; 1038 if (view && flips[p]) PetscCall(PetscPrintf(comm, "Flipping Proc+Comp %" PetscInt_FMT ":\n", p)); 1039 } 1040 for (PetscInt p = 0; p < size; ++p) displs[p + 1] = displs[p] + Nc[p]; 1041 } 1042 PetscCall(PetscMPIIntCast(Ncomp, &iNcomp)); 1043 PetscCallMPI(MPI_Scatterv(flips, Nc, displs, MPIU_BOOL, flipped, iNcomp, MPIU_BOOL, 0, comm)); 1044 PetscCall(PetscFree(flips)); 1045 } 1046 if (rank == 0) PetscCall(PetscBTDestroy(&flippedProcs)); 1047 PetscCall(PetscFree(N)); 1048 PetscCall(PetscFree4(recvcounts, displs, Nc, Noff)); 1049 PetscCall(PetscFree2(nrankComp, match)); 1050 1051 /* Decide whether to flip cells in each component */ 1052 for (PetscInt c = 0; c < cEnd - cStart; ++c) { 1053 if (flipped[cellComp[c]]) PetscCall(PetscBTNegate(flippedCells, c)); 1054 } 1055 PetscCall(PetscFree(flipped)); 1056 } 1057 if (view) { 1058 PetscViewer v; 1059 1060 PetscCall(PetscViewerASCIIGetStdout(comm, &v)); 1061 PetscCall(PetscViewerASCIIPushSynchronized(v)); 1062 PetscCall(PetscViewerASCIISynchronizedPrintf(v, "[%d]BT for parallel flipped cells:\n", rank)); 1063 PetscCall(PetscBTView(cEnd - cStart, flippedCells, v)); 1064 PetscCall(PetscViewerFlush(v)); 1065 PetscCall(PetscViewerASCIIPopSynchronized(v)); 1066 } 1067 // Reverse flipped cells in the mesh 1068 PetscViewer v; 1069 const PetscInt *degree = NULL; 1070 PetscInt *points; 1071 PetscInt pStart, pEnd; 1072 1073 if (view) { 1074 PetscCall(PetscViewerASCIIGetStdout(comm, &v)); 1075 PetscCall(PetscViewerASCIIPushSynchronized(v)); 1076 } 1077 PetscCall(DMPlexGetChart(dm, &pStart, &pEnd)); 1078 if (numRoots >= 0) { 1079 PetscCall(PetscSFComputeDegreeBegin(sf, °ree)); 1080 PetscCall(PetscSFComputeDegreeEnd(sf, °ree)); 1081 } 1082 PetscCall(PetscCalloc1(pEnd - pStart, &points)); 1083 for (PetscInt c = cStart; c < cEnd; ++c) { 1084 if (PetscBTLookup(flippedCells, c - cStart)) { 1085 const PetscInt cell = cells ? cells[c] : c; 1086 1087 PetscCall(DMPlexOrientPoint(dm, cell, -1)); 1088 if (degree && degree[cell]) points[cell] = 1; 1089 if (view) PetscCall(PetscViewerASCIISynchronizedPrintf(v, "[%d]Flipping cell %" PetscInt_FMT "%s\n", rank, cell, degree && degree[cell] ? " and sending to overlap" : "")); 1090 } 1091 } 1092 // Must propagate flips for cells in the overlap 1093 if (numRoots >= 0) { 1094 PetscCall(PetscSFBcastBegin(sf, MPIU_INT, points, points, MPI_SUM)); 1095 PetscCall(PetscSFBcastEnd(sf, MPIU_INT, points, points, MPI_SUM)); 1096 } 1097 for (PetscInt c = cStart; c < cEnd; ++c) { 1098 const PetscInt cell = cells ? cells[c] : c; 1099 1100 if (points[cell] && !PetscBTLookup(flippedCells, c - cStart)) { 1101 PetscCall(DMPlexOrientPoint(dm, cell, -1)); 1102 if (view) PetscCall(PetscViewerASCIISynchronizedPrintf(v, "[%d]Flipping cell %" PetscInt_FMT " through overlap\n", rank, cell)); 1103 } 1104 } 1105 if (view) { 1106 PetscCall(PetscViewerFlush(v)); 1107 PetscCall(PetscViewerASCIIPopSynchronized(v)); 1108 } 1109 PetscCall(PetscFree(points)); 1110 PetscCall(PetscBTDestroy(&flippedCells)); 1111 PetscCall(PetscFree2(numNeighbors, neighbors)); 1112 PetscCall(PetscFree3(rorntComp, lorntComp, locSupp)); 1113 PetscCall(PetscFree2(cellComp, faceComp)); 1114 PetscFunctionReturn(PETSC_SUCCESS); 1115 } 1116