1 #include <petsc/private/dmpleximpl.h> /*I "petscdmplex.h" I*/ 2 #include <petscsf.h> 3 #include <petsc/private/hashset.h> 4 5 typedef uint64_t ZCode; 6 7 PETSC_HASH_SET(ZSet, ZCode, PetscHash_UInt64, PetscHashEqual) 8 9 typedef struct { 10 PetscInt i, j, k; 11 } Ijk; 12 13 typedef struct { 14 Ijk eextent; 15 Ijk vextent; 16 PetscMPIInt comm_size; 17 ZCode *zstarts; 18 } ZLayout; 19 20 static unsigned ZCodeSplit1(ZCode z) 21 { 22 z = ((z & 01001001001001001) | ((z >> 2) & 02002002002002002) | ((z >> 4) & 04004004004004004)); 23 z = (z | (z >> 6) | (z >> 12)) & 0000000777000000777; 24 z = (z | (z >> 18)) & 0777777; 25 return (unsigned)z; 26 } 27 28 static ZCode ZEncode1(unsigned t) 29 { 30 ZCode z = t; 31 z = (z | (z << 18)) & 0777000000777; 32 z = (z | (z << 6) | (z << 12)) & 07007007007007007; 33 z = (z | (z << 2) | (z << 4)) & 0111111111111111111; 34 return z; 35 } 36 37 static Ijk ZCodeSplit(ZCode z) 38 { 39 Ijk c; 40 c.i = ZCodeSplit1(z >> 2); 41 c.j = ZCodeSplit1(z >> 1); 42 c.k = ZCodeSplit1(z >> 0); 43 return c; 44 } 45 46 static ZCode ZEncode(Ijk c) 47 { 48 ZCode z = (ZEncode1(c.i) << 2) | (ZEncode1(c.j) << 1) | ZEncode1(c.k); 49 return z; 50 } 51 52 static PetscBool IjkActive(Ijk extent, Ijk l) 53 { 54 if (l.i < extent.i && l.j < extent.j && l.k < extent.k) return PETSC_TRUE; 55 return PETSC_FALSE; 56 } 57 58 // Since element/vertex box extents are typically not equal powers of 2, Z codes that lie within the domain are not contiguous. 59 static PetscErrorCode ZLayoutCreate(PetscMPIInt size, const PetscInt eextent[3], const PetscInt vextent[3], ZLayout *zlayout) 60 { 61 ZLayout layout; 62 63 PetscFunctionBegin; 64 layout.eextent.i = eextent[0]; 65 layout.eextent.j = eextent[1]; 66 layout.eextent.k = eextent[2]; 67 layout.vextent.i = vextent[0]; 68 layout.vextent.j = vextent[1]; 69 layout.vextent.k = vextent[2]; 70 layout.comm_size = size; 71 PetscCall(PetscMalloc1(size + 1, &layout.zstarts)); 72 73 PetscInt total_elems = eextent[0] * eextent[1] * eextent[2]; 74 ZCode z = 0; 75 layout.zstarts[0] = 0; 76 for (PetscMPIInt r = 0; r < size; r++) { 77 PetscInt elems_needed = (total_elems / size) + (total_elems % size > r), count; 78 for (count = 0; count < elems_needed; z++) { 79 Ijk loc = ZCodeSplit(z); 80 if (IjkActive(layout.eextent, loc)) count++; 81 } 82 // Pick up any extra vertices in the Z ordering before the next rank's first owned element. 83 // 84 // TODO: This leads to poorly balanced vertices when eextent is a power of 2, since all the fringe vertices end up 85 // on the last rank. A possible solution is to balance the Z-order vertices independently from the cells, which will 86 // result in a lot of element closures being remote. We could finish marking boundary conditions, then do a round of 87 // vertex ownership smoothing (which would reorder and redistribute vertices without touching element distribution). 88 // Another would be to have an analytic ownership criteria for vertices in the fringe veextent - eextent. This would 89 // complicate the job of identifying an owner and its offset. 90 for (; z <= ZEncode(layout.vextent); z++) { 91 Ijk loc = ZCodeSplit(z); 92 if (IjkActive(layout.eextent, loc)) break; 93 } 94 layout.zstarts[r + 1] = z; 95 } 96 *zlayout = layout; 97 PetscFunctionReturn(0); 98 } 99 100 static PetscInt ZLayoutElementsOnRank(const ZLayout *layout, PetscMPIInt rank) 101 { 102 PetscInt remote_elem = 0; 103 for (ZCode rz = layout->zstarts[rank]; rz < layout->zstarts[rank + 1]; rz++) { 104 Ijk loc = ZCodeSplit(rz); 105 if (IjkActive(layout->eextent, loc)) remote_elem++; 106 } 107 return remote_elem; 108 } 109 110 PetscInt ZCodeFind(ZCode key, PetscInt n, const ZCode X[]) 111 { 112 PetscInt lo = 0, hi = n; 113 114 if (n == 0) return -1; 115 while (hi - lo > 1) { 116 PetscInt mid = lo + (hi - lo) / 2; 117 if (key < X[mid]) hi = mid; 118 else lo = mid; 119 } 120 return key == X[lo] ? lo : -(lo + (key > X[lo]) + 1); 121 } 122 123 static PetscErrorCode DMPlexCreateBoxMesh_Tensor_SFC_Periodicity_Private(DM dm, const ZLayout *layout, const ZCode *vert_z, PetscSegBuffer per_faces, PetscSegBuffer donor_face_closure, PetscSegBuffer my_donor_faces) 124 { 125 MPI_Comm comm; 126 size_t num_faces; 127 PetscInt dim, *faces, vStart, vEnd; 128 PetscMPIInt size; 129 ZCode *donor_verts, *donor_minz; 130 PetscSFNode *leaf; 131 132 PetscFunctionBegin; 133 PetscCall(PetscObjectGetComm((PetscObject)dm, &comm)); 134 PetscCallMPI(MPI_Comm_size(comm, &size)); 135 PetscCall(DMGetDimension(dm, &dim)); 136 const PetscInt csize = PetscPowInt(2, dim - 1); 137 PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd)); 138 PetscCall(PetscSegBufferGetSize(per_faces, &num_faces)); 139 PetscCall(PetscSegBufferExtractInPlace(per_faces, &faces)); 140 PetscCall(PetscSegBufferExtractInPlace(donor_face_closure, &donor_verts)); 141 PetscCall(PetscMalloc1(num_faces, &donor_minz)); 142 PetscCall(PetscMalloc1(num_faces, &leaf)); 143 for (PetscInt i = 0; i < (PetscInt)num_faces; i++) { 144 ZCode minz = donor_verts[i * csize]; 145 for (PetscInt j = 1; j < csize; j++) minz = PetscMin(minz, donor_verts[i * csize + j]); 146 donor_minz[i] = minz; 147 } 148 { 149 PetscBool sorted; 150 PetscCall(PetscSortedInt64(num_faces, (const PetscInt64 *)donor_minz, &sorted)); 151 PetscCheck(sorted, comm, PETSC_ERR_PLIB, "minz not sorted; periodicity in multiple dimensions not yet supported"); 152 } 153 for (PetscInt i = 0; i < (PetscInt)num_faces;) { 154 ZCode z = donor_minz[i]; 155 PetscInt remote_rank = ZCodeFind(z, size + 1, layout->zstarts), remote_count = 0; 156 if (remote_rank < 0) remote_rank = -(remote_rank + 1) - 1; 157 // Process all the vertices on this rank 158 for (ZCode rz = layout->zstarts[remote_rank]; rz < layout->zstarts[remote_rank + 1]; rz++) { 159 Ijk loc = ZCodeSplit(rz); 160 if (rz == z) { 161 leaf[i].rank = remote_rank; 162 leaf[i].index = remote_count; 163 i++; 164 if (i == (PetscInt)num_faces) break; 165 z = donor_minz[i]; 166 } 167 if (IjkActive(layout->vextent, loc)) remote_count++; 168 } 169 } 170 PetscCall(PetscFree(donor_minz)); 171 PetscSF sfper; 172 PetscCall(PetscSFCreate(PetscObjectComm((PetscObject)dm), &sfper)); 173 PetscCall(PetscSFSetGraph(sfper, vEnd - vStart, num_faces, PETSC_NULL, PETSC_USE_POINTER, leaf, PETSC_USE_POINTER)); 174 const PetscInt *my_donor_degree; 175 PetscCall(PetscSFComputeDegreeBegin(sfper, &my_donor_degree)); 176 PetscCall(PetscSFComputeDegreeEnd(sfper, &my_donor_degree)); 177 PetscInt num_multiroots = 0; 178 for (PetscInt i = 0; i < vEnd - vStart; i++) { 179 num_multiroots += my_donor_degree[i]; 180 if (my_donor_degree[i] == 0) continue; 181 PetscAssert(my_donor_degree[i] == 1, comm, PETSC_ERR_SUP, "Local vertex has multiple faces"); 182 } 183 PetscInt *my_donors, *donor_indices, *my_donor_indices; 184 size_t num_my_donors; 185 PetscCall(PetscSegBufferGetSize(my_donor_faces, &num_my_donors)); 186 PetscCheck((PetscInt)num_my_donors == num_multiroots, PETSC_COMM_SELF, PETSC_ERR_SUP, "Donor request does not match expected donors"); 187 PetscCall(PetscSegBufferExtractInPlace(my_donor_faces, &my_donors)); 188 PetscCall(PetscMalloc1(vEnd - vStart, &my_donor_indices)); 189 for (PetscInt i = 0; i < (PetscInt)num_my_donors; i++) { 190 PetscInt f = my_donors[i]; 191 PetscInt num_points, *points = NULL, minv = PETSC_MAX_INT; 192 PetscCall(DMPlexGetTransitiveClosure(dm, f, PETSC_TRUE, &num_points, &points)); 193 for (PetscInt j = 0; j < num_points; j++) { 194 PetscInt p = points[2 * j]; 195 if (p < vStart || vEnd <= p) continue; 196 minv = PetscMin(minv, p); 197 } 198 PetscCall(DMPlexRestoreTransitiveClosure(dm, f, PETSC_TRUE, &num_points, &points)); 199 PetscAssert(my_donor_degree[minv - vStart] == 1, comm, PETSC_ERR_SUP, "Local vertex not requested"); 200 my_donor_indices[minv - vStart] = f; 201 } 202 PetscCall(PetscMalloc1(num_faces, &donor_indices)); 203 PetscCall(PetscSFBcastBegin(sfper, MPIU_INT, my_donor_indices, donor_indices, MPI_REPLACE)); 204 PetscCall(PetscSFBcastEnd(sfper, MPIU_INT, my_donor_indices, donor_indices, MPI_REPLACE)); 205 PetscCall(PetscFree(my_donor_indices)); 206 // Modify our leafs so they point to donor faces instead of donor minz. Additionally, give them indices as faces. 207 for (PetscInt i = 0; i < (PetscInt)num_faces; i++) leaf[i].index = donor_indices[i]; 208 PetscCall(PetscFree(donor_indices)); 209 PetscInt pStart, pEnd; 210 PetscCall(DMPlexGetChart(dm, &pStart, &pEnd)); 211 PetscCall(PetscSFSetGraph(sfper, pEnd - pStart, num_faces, faces, PETSC_COPY_VALUES, leaf, PETSC_OWN_POINTER)); 212 PetscCall(PetscObjectSetName((PetscObject)sfper, "Periodic Faces")); 213 PetscCall(PetscSFViewFromOptions(sfper, NULL, "-sfper_view")); 214 215 PetscCall(PetscSegBufferDestroy(&per_faces)); 216 PetscCall(PetscSegBufferDestroy(&donor_face_closure)); 217 PetscCall(PetscSegBufferDestroy(&my_donor_faces)); 218 PetscCall(DMPlexSetPeriodicFaceSF(dm, sfper)); 219 PetscCall(PetscSFDestroy(&sfper)); 220 PetscFunctionReturn(0); 221 } 222 223 PetscErrorCode DMPlexCreateBoxMesh_Tensor_SFC_Internal(DM dm, PetscInt dim, const PetscInt faces[], const PetscReal lower[], const PetscReal upper[], const DMBoundaryType periodicity[], PetscBool interpolate) 224 { 225 PetscInt eextent[3] = {1, 1, 1}, vextent[3] = {1, 1, 1}; 226 const Ijk closure_1[] = { 227 {0, 0, 0}, 228 {1, 0, 0}, 229 }; 230 const Ijk closure_2[] = { 231 {0, 0, 0}, 232 {1, 0, 0}, 233 {1, 1, 0}, 234 {0, 1, 0}, 235 }; 236 const Ijk closure_3[] = { 237 {0, 0, 0}, 238 {0, 1, 0}, 239 {1, 1, 0}, 240 {1, 0, 0}, 241 {0, 0, 1}, 242 {1, 0, 1}, 243 {1, 1, 1}, 244 {0, 1, 1}, 245 }; 246 const Ijk *const closure_dim[] = {NULL, closure_1, closure_2, closure_3}; 247 // This must be kept consistent with DMPlexCreateCubeMesh_Internal 248 const PetscInt face_marker_1[] = {1, 2}; 249 const PetscInt face_marker_2[] = {4, 2, 1, 3}; 250 const PetscInt face_marker_3[] = {6, 5, 3, 4, 1, 2}; 251 const PetscInt *const face_marker_dim[] = {NULL, face_marker_1, face_marker_2, face_marker_3}; 252 253 PetscFunctionBegin; 254 PetscValidPointer(dm, 1); 255 PetscValidLogicalCollectiveInt(dm, dim, 2); 256 PetscCall(DMSetDimension(dm, dim)); 257 PetscMPIInt rank, size; 258 PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)dm), &size)); 259 PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank)); 260 for (PetscInt i = 0; i < dim; i++) { 261 eextent[i] = faces[i]; 262 vextent[i] = faces[i] + 1; 263 } 264 ZLayout layout; 265 PetscCall(ZLayoutCreate(size, eextent, vextent, &layout)); 266 PetscZSet vset; // set of all vertices in the closure of the owned elements 267 PetscCall(PetscZSetCreate(&vset)); 268 PetscInt local_elems = 0; 269 for (ZCode z = layout.zstarts[rank]; z < layout.zstarts[rank + 1]; z++) { 270 Ijk loc = ZCodeSplit(z); 271 if (IjkActive(layout.vextent, loc)) PetscZSetAdd(vset, z); 272 if (IjkActive(layout.eextent, loc)) { 273 local_elems++; 274 // Add all neighboring vertices to set 275 for (PetscInt n = 0; n < PetscPowInt(2, dim); n++) { 276 Ijk inc = closure_dim[dim][n]; 277 Ijk nloc = {loc.i + inc.i, loc.j + inc.j, loc.k + inc.k}; 278 ZCode v = ZEncode(nloc); 279 PetscZSetAdd(vset, v); 280 } 281 } 282 } 283 PetscInt local_verts, off = 0; 284 ZCode *vert_z; 285 PetscCall(PetscZSetGetSize(vset, &local_verts)); 286 PetscCall(PetscMalloc1(local_verts, &vert_z)); 287 PetscCall(PetscZSetGetElems(vset, &off, vert_z)); 288 PetscCall(PetscZSetDestroy(&vset)); 289 // ZCode is unsigned for bitwise convenience, but highest bit should never be set, so can interpret as signed 290 PetscCall(PetscSortInt64(local_verts, (PetscInt64 *)vert_z)); 291 292 PetscCall(DMPlexSetChart(dm, 0, local_elems + local_verts)); 293 for (PetscInt e = 0; e < local_elems; e++) PetscCall(DMPlexSetConeSize(dm, e, PetscPowInt(2, dim))); 294 PetscCall(DMSetUp(dm)); 295 { 296 PetscInt e = 0; 297 for (ZCode z = layout.zstarts[rank]; z < layout.zstarts[rank + 1]; z++) { 298 Ijk loc = ZCodeSplit(z); 299 if (!IjkActive(layout.eextent, loc)) continue; 300 PetscInt cone[8], orient[8] = {0}; 301 for (PetscInt n = 0; n < PetscPowInt(2, dim); n++) { 302 Ijk inc = closure_dim[dim][n]; 303 Ijk nloc = {loc.i + inc.i, loc.j + inc.j, loc.k + inc.k}; 304 ZCode v = ZEncode(nloc); 305 PetscInt ci = ZCodeFind(v, local_verts, vert_z); 306 PetscAssert(ci >= 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Did not find neighbor vertex in set"); 307 cone[n] = local_elems + ci; 308 } 309 PetscCall(DMPlexSetCone(dm, e, cone)); 310 PetscCall(DMPlexSetConeOrientation(dm, e, orient)); 311 e++; 312 } 313 } 314 315 if (0) { 316 DMLabel depth; 317 PetscCall(DMCreateLabel(dm, "depth")); 318 PetscCall(DMPlexGetDepthLabel(dm, &depth)); 319 PetscCall(DMLabelSetStratumBounds(depth, 0, local_elems, local_elems + local_verts)); 320 PetscCall(DMLabelSetStratumBounds(depth, 1, 0, local_elems)); 321 } else { 322 PetscCall(DMPlexSymmetrize(dm)); 323 PetscCall(DMPlexStratify(dm)); 324 } 325 { // Create point SF 326 PetscSF sf; 327 PetscSFCreate(PetscObjectComm((PetscObject)dm), &sf); 328 PetscInt owned_verts = ZCodeFind(layout.zstarts[rank + 1], local_verts, vert_z); 329 if (owned_verts < 0) owned_verts = -(owned_verts + 1); // We don't care whether the key was found 330 PetscInt num_ghosts = local_verts - owned_verts; // Due to sorting, owned vertices always come first 331 PetscInt *local_ghosts; 332 PetscSFNode *ghosts; 333 PetscCall(PetscMalloc1(num_ghosts, &local_ghosts)); 334 PetscCall(PetscMalloc1(num_ghosts, &ghosts)); 335 for (PetscInt i = 0; i < num_ghosts;) { 336 ZCode z = vert_z[owned_verts + i]; 337 PetscInt remote_rank = ZCodeFind(z, size + 1, layout.zstarts), remote_count = 0; 338 if (remote_rank < 0) remote_rank = -(remote_rank + 1) - 1; 339 // We have a new remote rank; find all the ghost indices (which are contiguous in vert_z) 340 341 // Count the elements on remote_rank 342 PetscInt remote_elem = ZLayoutElementsOnRank(&layout, remote_rank); 343 344 // Traverse vertices and make ghost links 345 for (ZCode rz = layout.zstarts[remote_rank]; rz < layout.zstarts[remote_rank + 1]; rz++) { 346 Ijk loc = ZCodeSplit(rz); 347 if (rz == z) { 348 local_ghosts[i] = local_elems + owned_verts + i; 349 ghosts[i].rank = remote_rank; 350 ghosts[i].index = remote_elem + remote_count; 351 i++; 352 if (i == num_ghosts) break; 353 z = vert_z[owned_verts + i]; 354 } 355 if (IjkActive(layout.vextent, loc)) remote_count++; 356 } 357 } 358 PetscCall(PetscSFSetGraph(sf, local_elems + local_verts, num_ghosts, local_ghosts, PETSC_OWN_POINTER, ghosts, PETSC_OWN_POINTER)); 359 PetscCall(PetscObjectSetName((PetscObject)sf, "SFC Point SF")); 360 PetscCall(DMSetPointSF(dm, sf)); 361 PetscCall(PetscSFDestroy(&sf)); 362 } 363 { 364 Vec coordinates; 365 PetscScalar *coords; 366 PetscSection coord_section; 367 PetscInt coord_size; 368 PetscCall(DMGetCoordinateSection(dm, &coord_section)); 369 PetscCall(PetscSectionSetNumFields(coord_section, 1)); 370 PetscCall(PetscSectionSetFieldComponents(coord_section, 0, dim)); 371 PetscCall(PetscSectionSetChart(coord_section, local_elems, local_elems + local_verts)); 372 for (PetscInt v = 0; v < local_verts; v++) { 373 PetscInt point = local_elems + v; 374 PetscCall(PetscSectionSetDof(coord_section, point, dim)); 375 PetscCall(PetscSectionSetFieldDof(coord_section, point, 0, dim)); 376 } 377 PetscCall(PetscSectionSetUp(coord_section)); 378 PetscCall(PetscSectionGetStorageSize(coord_section, &coord_size)); 379 PetscCall(VecCreate(PETSC_COMM_SELF, &coordinates)); 380 PetscCall(PetscObjectSetName((PetscObject)coordinates, "coordinates")); 381 PetscCall(VecSetSizes(coordinates, coord_size, PETSC_DETERMINE)); 382 PetscCall(VecSetBlockSize(coordinates, dim)); 383 PetscCall(VecSetType(coordinates, VECSTANDARD)); 384 PetscCall(VecGetArray(coordinates, &coords)); 385 for (PetscInt v = 0; v < local_verts; v++) { 386 Ijk loc = ZCodeSplit(vert_z[v]); 387 coords[v * dim + 0] = lower[0] + loc.i * (upper[0] - lower[0]) / layout.eextent.i; 388 if (dim > 1) coords[v * dim + 1] = lower[1] + loc.j * (upper[1] - lower[1]) / layout.eextent.j; 389 if (dim > 2) coords[v * dim + 2] = lower[2] + loc.k * (upper[2] - lower[2]) / layout.eextent.k; 390 } 391 PetscCall(VecRestoreArray(coordinates, &coords)); 392 PetscCall(DMSetCoordinatesLocal(dm, coordinates)); 393 PetscCall(VecDestroy(&coordinates)); 394 } 395 if (interpolate) { 396 PetscCall(DMPlexInterpolateInPlace_Internal(dm)); 397 398 DMLabel label; 399 PetscCall(DMCreateLabel(dm, "Face Sets")); 400 PetscCall(DMGetLabel(dm, "Face Sets", &label)); 401 PetscSegBuffer per_faces, donor_face_closure, my_donor_faces; 402 PetscCall(PetscSegBufferCreate(sizeof(PetscInt), 64, &per_faces)); 403 PetscCall(PetscSegBufferCreate(sizeof(PetscInt), 64, &my_donor_faces)); 404 PetscCall(PetscSegBufferCreate(sizeof(ZCode), 64 * PetscPowInt(2, dim), &donor_face_closure)); 405 PetscInt fStart, fEnd, vStart, vEnd; 406 PetscCall(DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd)); 407 PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd)); 408 for (PetscInt f = fStart; f < fEnd; f++) { 409 PetscInt npoints, *points = NULL, num_fverts = 0, fverts[8]; 410 PetscCall(DMPlexGetTransitiveClosure(dm, f, PETSC_TRUE, &npoints, &points)); 411 PetscInt bc_count[6] = {0}; 412 for (PetscInt i = 0; i < npoints; i++) { 413 PetscInt p = points[2 * i]; 414 if (p < vStart || vEnd <= p) continue; 415 fverts[num_fverts++] = p; 416 Ijk loc = ZCodeSplit(vert_z[p - vStart]); 417 // Convention here matches DMPlexCreateCubeMesh_Internal 418 bc_count[0] += loc.i == 0; 419 bc_count[1] += loc.i == layout.vextent.i - 1; 420 bc_count[2] += loc.j == 0; 421 bc_count[3] += loc.j == layout.vextent.j - 1; 422 bc_count[4] += loc.k == 0; 423 bc_count[5] += loc.k == layout.vextent.k - 1; 424 } 425 PetscCall(DMPlexRestoreTransitiveClosure(dm, f, PETSC_TRUE, &npoints, &points)); 426 for (PetscInt bc = 0, bc_match = 0; bc < 2 * dim; bc++) { 427 if (bc_count[bc] == PetscPowInt(2, dim - 1)) { 428 if (periodicity[bc / 2] == DM_BOUNDARY_PERIODIC) { 429 PetscInt *put; 430 if (bc % 2 == 0) { // donor face; no label 431 PetscCall(PetscSegBufferGet(my_donor_faces, 1, &put)); 432 *put = f; 433 } else { // periodic face 434 PetscCall(PetscSegBufferGet(per_faces, 1, &put)); 435 *put = f; 436 ZCode *zput; 437 PetscCall(PetscSegBufferGet(donor_face_closure, num_fverts, &zput)); 438 for (PetscInt i = 0; i < num_fverts; i++) { 439 Ijk loc = ZCodeSplit(vert_z[fverts[i] - vStart]); 440 switch (bc / 2) { 441 case 0: 442 loc.i = 0; 443 break; 444 case 1: 445 loc.j = 0; 446 break; 447 case 2: 448 loc.k = 0; 449 break; 450 } 451 *zput++ = ZEncode(loc); 452 } 453 } 454 continue; 455 } 456 PetscAssert(bc_match == 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Face matches multiple face sets"); 457 PetscCall(DMLabelSetValue(label, f, face_marker_dim[dim][bc])); 458 bc_match++; 459 } 460 } 461 } 462 PetscCall(DMPlexCreateBoxMesh_Tensor_SFC_Periodicity_Private(dm, &layout, vert_z, per_faces, donor_face_closure, my_donor_faces)); 463 // Ensure that the Coordinate DM has our new boundary labels 464 DM cdm; 465 PetscCall(DMGetCoordinateDM(dm, &cdm)); 466 PetscCall(DMCopyLabels(dm, cdm, PETSC_COPY_VALUES, PETSC_FALSE, DM_COPY_LABELS_FAIL)); 467 PetscSF sfper; 468 PetscCall(DMPlexGetPeriodicFaceSF(dm, &sfper)); 469 PetscCall(DMPlexSetPeriodicFaceSF(cdm, sfper)); 470 } 471 PetscCall(PetscFree(layout.zstarts)); 472 PetscCall(PetscFree(vert_z)); 473 PetscFunctionReturn(0); 474 } 475 476 /*@ 477 DMPlexSetPeriodicFaceSF - Express periodicity from an existing mesh 478 479 Logically collective 480 481 Input Parameters: 482 + dm - The `DMPLEX` on which to set periodicity 483 - face_sf - SF in which roots are (owned) donor faces and leaves are faces that must be matched to a (possibly remote) donor face. 484 485 Level: advanced 486 487 Notes: 488 489 One can use `-dm_plex_box_sfc` to use this mode of periodicity, wherein the periodic points are distinct both globally 490 and locally, but are paired when creating a global dof space. 491 492 .seealso: [](chapter_unstructured), `DMPLEX`, `DMGetGlobalSection()`, `DMPlexGetPeriodicFaceSF()` 493 @*/ 494 PetscErrorCode DMPlexSetPeriodicFaceSF(DM dm, PetscSF face_sf) 495 { 496 DM_Plex *plex = (DM_Plex *)dm->data; 497 PetscFunctionBegin; 498 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 499 PetscCall(PetscObjectReference((PetscObject)face_sf)); 500 PetscCall(PetscSFDestroy(&plex->periodic.face_sf)); 501 plex->periodic.face_sf = face_sf; 502 PetscFunctionReturn(0); 503 } 504 505 /*@ 506 DMPlexGetPeriodicFaceSF - Obtain periodicity for a mesh 507 508 Logically collective 509 510 Input Parameters: 511 . dm - The `DMPLEX` for which to obtain periodic relation 512 513 Output Parameters: 514 . face_sf - SF in which roots are (owned) donor faces and leaves are faces that must be matched to a (possibly remote) donor face. 515 516 Level: advanced 517 518 .seealso: [](chapter_unstructured), `DMPLEX`, `DMGetGlobalSection()`, `DMPlexSetPeriodicFaceSF()` 519 @*/ 520 PetscErrorCode DMPlexGetPeriodicFaceSF(DM dm, PetscSF *face_sf) 521 { 522 DM_Plex *plex = (DM_Plex *)dm->data; 523 PetscFunctionBegin; 524 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 525 *face_sf = plex->periodic.face_sf; 526 PetscFunctionReturn(0); 527 } 528