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 *layout) 60 { 61 PetscFunctionBegin; 62 layout->eextent.i = eextent[0]; 63 layout->eextent.j = eextent[1]; 64 layout->eextent.k = eextent[2]; 65 layout->vextent.i = vextent[0]; 66 layout->vextent.j = vextent[1]; 67 layout->vextent.k = vextent[2]; 68 layout->comm_size = size; 69 layout->zstarts = NULL; 70 PetscCall(PetscMalloc1(size + 1, &layout->zstarts)); 71 72 PetscInt total_elems = eextent[0] * eextent[1] * eextent[2]; 73 ZCode z = 0; 74 layout->zstarts[0] = 0; 75 for (PetscMPIInt r = 0; r < size; r++) { 76 PetscInt elems_needed = (total_elems / size) + (total_elems % size > r), count; 77 for (count = 0; count < elems_needed; z++) { 78 Ijk loc = ZCodeSplit(z); 79 if (IjkActive(layout->eextent, loc)) count++; 80 } 81 // Pick up any extra vertices in the Z ordering before the next rank's first owned element. 82 // 83 // This leads to poorly balanced vertices when eextent is a power of 2, since all the fringe vertices end up 84 // on the last rank. A possible solution is to balance the Z-order vertices independently from the cells, which will 85 // result in a lot of element closures being remote. We could finish marking boundary conditions, then do a round of 86 // vertex ownership smoothing (which would reorder and redistribute vertices without touching element distribution). 87 // Another would be to have an analytic ownership criteria for vertices in the fringe veextent - eextent. This would 88 // complicate the job of identifying an owner and its offset. 89 // 90 // The current recommended approach is to let `-dm_distribute 1` (default) resolve vertex ownership. This is 91 // *mandatory* with isoperiodicity (except in special cases) to remove standed vertices from local spaces. Here's 92 // the issue: 93 // 94 // Consider this partition on rank 0 (left) and rank 1. 95 // 96 // 4 -------- 5 -- 14 --10 -- 21 --11 97 // | | | 98 // 7 -- 16 -- 8 | | | 99 // | | 3 ------- 7 ------- 9 100 // | | | | 101 // 4 -------- 6 ------ 10 | | 102 // | | | 6 -- 16 -- 8 103 // | | | 104 // 3 ---11--- 5 --18-- 9 105 // 106 // The periodic face SF looks like 107 // [0] Number of roots=21, leaves=1, remote ranks=1 108 // [0] 16 <- (0,11) 109 // [1] Number of roots=22, leaves=2, remote ranks=2 110 // [1] 14 <- (0,18) 111 // [1] 21 <- (1,16) 112 // 113 // In handling face (0,16), rank 0 learns that (0,7) and (0,8) map to (0,3) and (0,5) respectively, thus we won't use 114 // the point SF links to (1,4) and (1,5). Rank 1 learns about the periodic mapping of (1,5) while handling face 115 // (1,14), but never learns that vertex (1,4) has been mapped to (0,3) by face (0,16). 116 // 117 // We can relatively easily inform vertex (1,4) of this mapping, but it stays in rank 1's local space despite not 118 // being in the closure and thus not being contributed to. This would be mostly harmless except that some viewer 119 // routines expect all local points to be somehow significant. It is not easy to analytically remove the (1,4) 120 // vertex because the point SF and isoperiodic face SF would need to be updated to account for removal of the 121 // stranded vertices. 122 for (; z <= ZEncode(layout->vextent); z++) { 123 Ijk loc = ZCodeSplit(z); 124 if (IjkActive(layout->eextent, loc)) break; 125 } 126 layout->zstarts[r + 1] = z; 127 } 128 PetscFunctionReturn(0); 129 } 130 131 static PetscInt ZLayoutElementsOnRank(const ZLayout *layout, PetscMPIInt rank) 132 { 133 PetscInt remote_elem = 0; 134 for (ZCode rz = layout->zstarts[rank]; rz < layout->zstarts[rank + 1]; rz++) { 135 Ijk loc = ZCodeSplit(rz); 136 if (IjkActive(layout->eextent, loc)) remote_elem++; 137 } 138 return remote_elem; 139 } 140 141 PetscInt ZCodeFind(ZCode key, PetscInt n, const ZCode X[]) 142 { 143 PetscInt lo = 0, hi = n; 144 145 if (n == 0) return -1; 146 while (hi - lo > 1) { 147 PetscInt mid = lo + (hi - lo) / 2; 148 if (key < X[mid]) hi = mid; 149 else lo = mid; 150 } 151 return key == X[lo] ? lo : -(lo + (key > X[lo]) + 1); 152 } 153 154 static PetscErrorCode DMPlexCreateBoxMesh_Tensor_SFC_Periodicity_Private(DM dm, const ZLayout *layout, const ZCode *vert_z, PetscSegBuffer per_faces, const PetscReal *lower, const PetscReal *upper, const DMBoundaryType *periodicity, PetscSegBuffer donor_face_closure, PetscSegBuffer my_donor_faces) 155 { 156 MPI_Comm comm; 157 size_t num_faces; 158 PetscInt dim, *faces, vStart, vEnd; 159 PetscMPIInt size; 160 ZCode *donor_verts, *donor_minz; 161 PetscSFNode *leaf; 162 163 PetscFunctionBegin; 164 PetscCall(PetscObjectGetComm((PetscObject)dm, &comm)); 165 PetscCallMPI(MPI_Comm_size(comm, &size)); 166 PetscCall(DMGetDimension(dm, &dim)); 167 const PetscInt csize = PetscPowInt(2, dim - 1); 168 PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd)); 169 PetscCall(PetscSegBufferGetSize(per_faces, &num_faces)); 170 PetscCall(PetscSegBufferExtractInPlace(per_faces, &faces)); 171 PetscCall(PetscSegBufferExtractInPlace(donor_face_closure, &donor_verts)); 172 PetscCall(PetscMalloc1(num_faces, &donor_minz)); 173 PetscCall(PetscMalloc1(num_faces, &leaf)); 174 for (PetscInt i = 0; i < (PetscInt)num_faces; i++) { 175 ZCode minz = donor_verts[i * csize]; 176 for (PetscInt j = 1; j < csize; j++) minz = PetscMin(minz, donor_verts[i * csize + j]); 177 donor_minz[i] = minz; 178 } 179 { 180 PetscBool sorted; 181 PetscCall(PetscSortedInt64(num_faces, (const PetscInt64 *)donor_minz, &sorted)); 182 PetscCheck(sorted, PETSC_COMM_SELF, PETSC_ERR_PLIB, "minz not sorted; periodicity in multiple dimensions not yet supported"); 183 } 184 for (PetscInt i = 0; i < (PetscInt)num_faces;) { 185 ZCode z = donor_minz[i]; 186 PetscInt remote_rank = ZCodeFind(z, size + 1, layout->zstarts), remote_count = 0; 187 if (remote_rank < 0) remote_rank = -(remote_rank + 1) - 1; 188 // Process all the vertices on this rank 189 for (ZCode rz = layout->zstarts[remote_rank]; rz < layout->zstarts[remote_rank + 1]; rz++) { 190 Ijk loc = ZCodeSplit(rz); 191 if (rz == z) { 192 leaf[i].rank = remote_rank; 193 leaf[i].index = remote_count; 194 i++; 195 if (i == (PetscInt)num_faces) break; 196 z = donor_minz[i]; 197 } 198 if (IjkActive(layout->vextent, loc)) remote_count++; 199 } 200 } 201 PetscCall(PetscFree(donor_minz)); 202 PetscSF sfper; 203 PetscCall(PetscSFCreate(PetscObjectComm((PetscObject)dm), &sfper)); 204 PetscCall(PetscSFSetGraph(sfper, vEnd - vStart, num_faces, PETSC_NULL, PETSC_USE_POINTER, leaf, PETSC_USE_POINTER)); 205 const PetscInt *my_donor_degree; 206 PetscCall(PetscSFComputeDegreeBegin(sfper, &my_donor_degree)); 207 PetscCall(PetscSFComputeDegreeEnd(sfper, &my_donor_degree)); 208 PetscInt num_multiroots = 0; 209 for (PetscInt i = 0; i < vEnd - vStart; i++) { 210 num_multiroots += my_donor_degree[i]; 211 if (my_donor_degree[i] == 0) continue; 212 PetscAssert(my_donor_degree[i] == 1, PETSC_COMM_SELF, PETSC_ERR_SUP, "Local vertex has multiple faces"); 213 } 214 PetscInt *my_donors, *donor_indices, *my_donor_indices; 215 size_t num_my_donors; 216 PetscCall(PetscSegBufferGetSize(my_donor_faces, &num_my_donors)); 217 PetscCheck((PetscInt)num_my_donors == num_multiroots, PETSC_COMM_SELF, PETSC_ERR_SUP, "Donor request does not match expected donors"); 218 PetscCall(PetscSegBufferExtractInPlace(my_donor_faces, &my_donors)); 219 PetscCall(PetscMalloc1(vEnd - vStart, &my_donor_indices)); 220 for (PetscInt i = 0; i < (PetscInt)num_my_donors; i++) { 221 PetscInt f = my_donors[i]; 222 PetscInt num_points, *points = NULL, minv = PETSC_MAX_INT; 223 PetscCall(DMPlexGetTransitiveClosure(dm, f, PETSC_TRUE, &num_points, &points)); 224 for (PetscInt j = 0; j < num_points; j++) { 225 PetscInt p = points[2 * j]; 226 if (p < vStart || vEnd <= p) continue; 227 minv = PetscMin(minv, p); 228 } 229 PetscCall(DMPlexRestoreTransitiveClosure(dm, f, PETSC_TRUE, &num_points, &points)); 230 PetscAssert(my_donor_degree[minv - vStart] == 1, PETSC_COMM_SELF, PETSC_ERR_SUP, "Local vertex not requested"); 231 my_donor_indices[minv - vStart] = f; 232 } 233 PetscCall(PetscMalloc1(num_faces, &donor_indices)); 234 PetscCall(PetscSFBcastBegin(sfper, MPIU_INT, my_donor_indices, donor_indices, MPI_REPLACE)); 235 PetscCall(PetscSFBcastEnd(sfper, MPIU_INT, my_donor_indices, donor_indices, MPI_REPLACE)); 236 PetscCall(PetscFree(my_donor_indices)); 237 // Modify our leafs so they point to donor faces instead of donor minz. Additionally, give them indices as faces. 238 for (PetscInt i = 0; i < (PetscInt)num_faces; i++) leaf[i].index = donor_indices[i]; 239 PetscCall(PetscFree(donor_indices)); 240 PetscInt pStart, pEnd; 241 PetscCall(DMPlexGetChart(dm, &pStart, &pEnd)); 242 PetscCall(PetscSFSetGraph(sfper, pEnd - pStart, num_faces, faces, PETSC_COPY_VALUES, leaf, PETSC_OWN_POINTER)); 243 PetscCall(PetscObjectSetName((PetscObject)sfper, "Z-order Isoperiodic Faces")); 244 245 PetscCall(DMPlexSetIsoperiodicFaceSF(dm, sfper)); 246 247 PetscScalar t[4][4] = {{0}}; 248 t[0][0] = 1; 249 t[1][1] = 1; 250 t[2][2] = 1; 251 t[3][3] = 1; 252 for (PetscInt i = 0; i < dim; i++) 253 if (periodicity[i] == DM_BOUNDARY_PERIODIC) t[i][3] = upper[i] - lower[i]; 254 PetscCall(DMPlexSetIsoperiodicFaceTransform(dm, &t[0][0])); 255 PetscCall(PetscSFDestroy(&sfper)); 256 PetscFunctionReturn(0); 257 } 258 259 // This is a DMGlobalToLocalHook that applies the affine offsets. When extended for rotated periodicity, it'll need to 260 // apply a rotatonal transform and similar operations will be needed for fields (e.g., to rotate a velocity vector). 261 // We use this crude approach here so we don't have to write new GPU kernels yet. 262 static PetscErrorCode DMCoordAddPeriodicOffsets_Private(DM dm, Vec g, InsertMode mode, Vec l, void *ctx) 263 { 264 PetscFunctionBegin; 265 PetscCall(VecScatterBegin(dm->periodic.affine_to_local, dm->periodic.affine, l, ADD_VALUES, SCATTER_FORWARD)); 266 PetscCall(VecScatterEnd(dm->periodic.affine_to_local, dm->periodic.affine, l, ADD_VALUES, SCATTER_FORWARD)); 267 PetscFunctionReturn(0); 268 } 269 270 // Start with an SF for a positive depth (e.g., faces) and create a new SF for matched closure. The caller must ensure 271 // that both the donor (root) face and the periodic (leaf) face have consistent orientation, meaning that their closures 272 // are isomorphic. It may be useful/necessary for this restriction to be loosened. 273 // 274 // Output Arguments: 275 // 276 // + closure_sf - augmented point SF (see `DMGetPointSF()`) that includes the faces and all points in its closure. This 277 // can be used to create a global section and section SF. 278 // - is_points - index set for just the points in the closure of `face_sf`. These may be used to apply an affine 279 // transformation to periodic dofs; see DMPeriodicCoordinateSetUp_Internal(). 280 // 281 static PetscErrorCode DMPlexCreateIsoperiodicPointSF_Private(DM dm, PetscSF face_sf, PetscSF *closure_sf, IS *is_points) 282 { 283 MPI_Comm comm; 284 PetscInt nroots, nleaves, npoints; 285 const PetscInt *filocal, *pilocal; 286 const PetscSFNode *firemote, *piremote; 287 PetscMPIInt rank; 288 PetscSF point_sf; 289 290 PetscFunctionBegin; 291 PetscCall(PetscObjectGetComm((PetscObject)dm, &comm)); 292 PetscCallMPI(MPI_Comm_rank(comm, &rank)); 293 PetscCall(PetscSFGetGraph(face_sf, &nroots, &nleaves, &filocal, &firemote)); 294 PetscCall(DMGetPointSF(dm, &point_sf)); // Point SF has remote points 295 PetscCall(PetscSFGetGraph(point_sf, NULL, &npoints, &pilocal, &piremote)); 296 PetscInt *rootdata, *leafdata; 297 PetscCall(PetscCalloc2(2 * nroots, &rootdata, 2 * nroots, &leafdata)); 298 for (PetscInt i = 0; i < nleaves; i++) { 299 PetscInt point = filocal[i], cl_size, *closure = NULL; 300 PetscCall(DMPlexGetTransitiveClosure(dm, point, PETSC_TRUE, &cl_size, &closure)); 301 leafdata[point] = cl_size - 1; 302 PetscCall(DMPlexRestoreTransitiveClosure(dm, point, PETSC_TRUE, &cl_size, &closure)); 303 } 304 PetscCall(PetscSFReduceBegin(face_sf, MPIU_INT, leafdata, rootdata + nroots, MPIU_SUM)); 305 PetscCall(PetscSFReduceEnd(face_sf, MPIU_INT, leafdata, rootdata + nroots, MPIU_SUM)); 306 307 PetscInt root_offset = 0; 308 for (PetscInt p = 0; p < nroots; p++) { 309 const PetscInt *donor_dof = rootdata + nroots; 310 if (donor_dof[p] == 0) { 311 rootdata[2 * p] = -1; 312 rootdata[2 * p + 1] = -1; 313 continue; 314 } 315 PetscInt cl_size; 316 PetscInt *closure = NULL; 317 PetscCall(DMPlexGetTransitiveClosure(dm, p, PETSC_TRUE, &cl_size, &closure)); 318 // cl_size - 1 = points not including self 319 PetscAssert(donor_dof[p] == cl_size - 1, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Reduced leaf cone sizes do not match root cone sizes"); 320 rootdata[2 * p] = root_offset; 321 rootdata[2 * p + 1] = cl_size - 1; 322 root_offset += cl_size - 1; 323 PetscCall(DMPlexRestoreTransitiveClosure(dm, p, PETSC_TRUE, &cl_size, &closure)); 324 } 325 PetscCall(PetscSFBcastBegin(face_sf, MPIU_2INT, rootdata, leafdata, MPI_REPLACE)); 326 PetscCall(PetscSFBcastEnd(face_sf, MPIU_2INT, rootdata, leafdata, MPI_REPLACE)); 327 // Count how many leaves we need to communicate the closures 328 PetscInt leaf_offset = 0; 329 for (PetscInt i = 0; i < nleaves; i++) { 330 PetscInt point = filocal[i]; 331 if (leafdata[2 * point + 1] < 0) continue; 332 leaf_offset += leafdata[2 * point + 1]; 333 } 334 335 PetscSFNode *closure_leaf; 336 PetscCall(PetscMalloc1(leaf_offset, &closure_leaf)); 337 leaf_offset = 0; 338 for (PetscInt i = 0; i < nleaves; i++) { 339 PetscInt point = filocal[i]; 340 PetscInt cl_size = leafdata[2 * point + 1]; 341 if (cl_size < 0) continue; 342 for (PetscInt j = 0; j < cl_size; j++) { 343 closure_leaf[leaf_offset].rank = firemote[i].rank; 344 closure_leaf[leaf_offset].index = leafdata[2 * point] + j; 345 leaf_offset++; 346 } 347 } 348 349 PetscSF sf_closure; 350 PetscCall(PetscSFCreate(comm, &sf_closure)); 351 PetscCall(PetscSFSetGraph(sf_closure, root_offset, leaf_offset, NULL, PETSC_USE_POINTER, closure_leaf, PETSC_OWN_POINTER)); 352 353 // Pack root buffer with owner for every point in the root cones 354 PetscSFNode *donor_closure; 355 PetscCall(PetscCalloc1(root_offset, &donor_closure)); 356 root_offset = 0; 357 for (PetscInt p = 0; p < nroots; p++) { 358 if (rootdata[2 * p] < 0) continue; 359 PetscInt cl_size; 360 PetscInt *closure = NULL; 361 PetscCall(DMPlexGetTransitiveClosure(dm, p, PETSC_TRUE, &cl_size, &closure)); 362 for (PetscInt j = 1; j < cl_size; j++) { 363 PetscInt c = closure[2 * j]; 364 if (pilocal) { 365 PetscInt found = -1; 366 if (npoints > 0) PetscCall(PetscFindInt(c, npoints, pilocal, &found)); 367 if (found >= 0) { 368 donor_closure[root_offset++] = piremote[found]; 369 continue; 370 } 371 } 372 // we own c 373 donor_closure[root_offset].rank = rank; 374 donor_closure[root_offset].index = c; 375 root_offset++; 376 } 377 PetscCall(DMPlexRestoreTransitiveClosure(dm, p, PETSC_TRUE, &cl_size, &closure)); 378 } 379 380 PetscSFNode *leaf_donor_closure; 381 PetscCall(PetscMalloc1(leaf_offset, &leaf_donor_closure)); 382 PetscCall(PetscSFBcastBegin(sf_closure, MPIU_2INT, donor_closure, leaf_donor_closure, MPI_REPLACE)); 383 PetscCall(PetscSFBcastEnd(sf_closure, MPIU_2INT, donor_closure, leaf_donor_closure, MPI_REPLACE)); 384 PetscCall(PetscSFDestroy(&sf_closure)); 385 PetscCall(PetscFree(donor_closure)); 386 387 PetscSFNode *new_iremote; 388 PetscCall(PetscCalloc1(nroots, &new_iremote)); 389 for (PetscInt i = 0; i < nroots; i++) new_iremote[i].rank = -1; 390 // Walk leaves and match vertices 391 leaf_offset = 0; 392 for (PetscInt i = 0; i < nleaves; i++) { 393 PetscInt point = filocal[i], cl_size; 394 PetscInt *closure = NULL; 395 PetscCall(DMPlexGetTransitiveClosure(dm, point, PETSC_TRUE, &cl_size, &closure)); 396 for (PetscInt j = 1; j < cl_size; j++) { // TODO: should we send donor edge orientations so we can flip for consistency? 397 PetscInt c = closure[2 * j]; 398 PetscSFNode lc = leaf_donor_closure[leaf_offset]; 399 // printf("[%d] face %d.%d: %d ?-- (%d,%d)\n", rank, point, j, c, lc.rank, lc.index); 400 if (new_iremote[c].rank == -1) { 401 new_iremote[c] = lc; 402 } else PetscCheck(new_iremote[c].rank == lc.rank && new_iremote[c].index == lc.index, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Mismatched cone ordering between faces"); 403 leaf_offset++; 404 } 405 PetscCall(DMPlexRestoreTransitiveClosure(dm, point, PETSC_TRUE, &cl_size, &closure)); 406 } 407 PetscCall(PetscFree(leaf_donor_closure)); 408 409 // Include face points in closure SF 410 for (PetscInt i = 0; i < nleaves; i++) new_iremote[filocal[i]] = firemote[i]; 411 // consolidate leaves 412 PetscInt num_new_leaves = 0; 413 for (PetscInt i = 0; i < nroots; i++) { 414 if (new_iremote[i].rank == -1) continue; 415 new_iremote[num_new_leaves] = new_iremote[i]; 416 leafdata[num_new_leaves] = i; 417 num_new_leaves++; 418 } 419 PetscCall(ISCreateGeneral(PETSC_COMM_SELF, num_new_leaves, leafdata, PETSC_COPY_VALUES, is_points)); 420 421 PetscSF csf; 422 PetscCall(PetscSFCreate(comm, &csf)); 423 PetscCall(PetscSFSetGraph(csf, nroots, num_new_leaves, leafdata, PETSC_COPY_VALUES, new_iremote, PETSC_COPY_VALUES)); 424 PetscCall(PetscFree(new_iremote)); // copy and delete because new_iremote is longer than it needs to be 425 PetscCall(PetscFree2(rootdata, leafdata)); 426 427 if (npoints < 0) { // empty point_sf 428 *closure_sf = csf; 429 } else { 430 PetscCall(PetscSFMerge(point_sf, csf, closure_sf)); 431 PetscCall(PetscSFDestroy(&csf)); 432 } 433 PetscCall(PetscObjectSetName((PetscObject)*closure_sf, "Composed Periodic Points")); 434 PetscFunctionReturn(0); 435 } 436 437 static PetscErrorCode DMGetIsoperiodicPointSF_Plex(DM dm, PetscSF *sf) 438 { 439 DM_Plex *plex = (DM_Plex *)dm->data; 440 441 PetscFunctionBegin; 442 if (!plex->periodic.composed_sf) { 443 PetscSF face_sf = plex->periodic.face_sf; 444 445 PetscCall(DMPlexCreateIsoperiodicPointSF_Private(dm, face_sf, &plex->periodic.composed_sf, &plex->periodic.periodic_points)); 446 } 447 if (sf) *sf = plex->periodic.composed_sf; 448 PetscFunctionReturn(0); 449 } 450 451 PetscErrorCode DMPlexMigrateIsoperiodicFaceSF_Internal(DM old_dm, DM dm, PetscSF sf_migration) 452 { 453 DM_Plex *plex = (DM_Plex *)old_dm->data; 454 PetscSF sf_point; 455 PetscMPIInt rank; 456 457 PetscFunctionBegin; 458 if (!plex->periodic.face_sf) PetscFunctionReturn(0); 459 PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank)); 460 PetscCall(DMGetPointSF(dm, &sf_point)); 461 PetscInt old_npoints, new_npoints, old_nleaf, new_nleaf, point_nleaf; 462 PetscSFNode *new_leafdata, *rootdata, *leafdata; 463 const PetscInt *old_local, *point_local; 464 const PetscSFNode *old_remote, *point_remote; 465 PetscCall(PetscSFGetGraph(plex->periodic.face_sf, &old_npoints, &old_nleaf, &old_local, &old_remote)); 466 PetscCall(PetscSFGetGraph(sf_migration, NULL, &new_nleaf, NULL, NULL)); 467 PetscCall(PetscSFGetGraph(sf_point, &new_npoints, &point_nleaf, &point_local, &point_remote)); 468 PetscAssert(new_nleaf == new_npoints, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Expected migration leaf space to match new point root space"); 469 PetscCall(PetscMalloc3(old_npoints, &rootdata, old_npoints, &leafdata, new_npoints, &new_leafdata)); 470 471 // Fill new_leafdata with new owners of all points 472 for (PetscInt i = 0; i < new_npoints; i++) { 473 new_leafdata[i].rank = rank; 474 new_leafdata[i].index = i; 475 } 476 for (PetscInt i = 0; i < point_nleaf; i++) { 477 PetscInt j = point_local[i]; 478 new_leafdata[j] = point_remote[i]; 479 } 480 // REPLACE is okay because every leaf agrees about the new owners 481 PetscCall(PetscSFReduceBegin(sf_migration, MPIU_2INT, new_leafdata, rootdata, MPI_REPLACE)); 482 PetscCall(PetscSFReduceEnd(sf_migration, MPIU_2INT, new_leafdata, rootdata, MPI_REPLACE)); 483 // rootdata now contains the new owners 484 485 // Send to leaves of old space 486 for (PetscInt i = 0; i < old_npoints; i++) { 487 leafdata[i].rank = -1; 488 leafdata[i].index = -1; 489 } 490 PetscCall(PetscSFBcastBegin(plex->periodic.face_sf, MPIU_2INT, rootdata, leafdata, MPI_REPLACE)); 491 PetscCall(PetscSFBcastEnd(plex->periodic.face_sf, MPIU_2INT, rootdata, leafdata, MPI_REPLACE)); 492 493 // Send to new leaf space 494 PetscCall(PetscSFBcastBegin(sf_migration, MPIU_2INT, leafdata, new_leafdata, MPI_REPLACE)); 495 PetscCall(PetscSFBcastEnd(sf_migration, MPIU_2INT, leafdata, new_leafdata, MPI_REPLACE)); 496 497 PetscInt nface = 0, *new_local; 498 PetscSFNode *new_remote; 499 for (PetscInt i = 0; i < new_npoints; i++) nface += (new_leafdata[i].rank >= 0); 500 PetscCall(PetscMalloc1(nface, &new_local)); 501 PetscCall(PetscMalloc1(nface, &new_remote)); 502 nface = 0; 503 for (PetscInt i = 0; i < new_npoints; i++) { 504 if (new_leafdata[i].rank == -1) continue; 505 new_local[nface] = i; 506 new_remote[nface] = new_leafdata[i]; 507 nface++; 508 } 509 PetscCall(PetscFree3(rootdata, leafdata, new_leafdata)); 510 PetscSF sf_face; 511 PetscCall(PetscSFCreate(PetscObjectComm((PetscObject)dm), &sf_face)); 512 PetscCall(PetscSFSetGraph(sf_face, new_npoints, nface, new_local, PETSC_OWN_POINTER, new_remote, PETSC_OWN_POINTER)); 513 PetscCall(PetscObjectSetName((PetscObject)sf_face, "Migrated Isoperiodic Faces")); 514 PetscCall(DMPlexSetIsoperiodicFaceSF(dm, sf_face)); 515 PetscCall(DMPlexSetIsoperiodicFaceTransform(dm, &plex->periodic.transform[0][0])); 516 PetscCall(PetscSFDestroy(&sf_face)); 517 PetscFunctionReturn(0); 518 } 519 520 PetscErrorCode DMPeriodicCoordinateSetUp_Internal(DM dm) 521 { 522 DM_Plex *plex = (DM_Plex *)dm->data; 523 PetscFunctionBegin; 524 if (!plex->periodic.face_sf) PetscFunctionReturn(0); 525 PetscCall(DMGetIsoperiodicPointSF_Plex(dm, NULL)); 526 PetscCall(PetscObjectComposeFunction((PetscObject)dm, "DMGetIsoperiodicPointSF_C", DMGetIsoperiodicPointSF_Plex)); 527 528 PetscInt dim; 529 PetscCall(DMGetDimension(dm, &dim)); 530 size_t count; 531 IS isdof; 532 { 533 PetscInt npoints; 534 const PetscInt *points; 535 IS is = plex->periodic.periodic_points; 536 PetscSegBuffer seg; 537 PetscSection section; 538 PetscCall(DMGetLocalSection(dm, §ion)); 539 PetscCall(PetscSegBufferCreate(sizeof(PetscInt), 32, &seg)); 540 PetscCall(ISGetSize(is, &npoints)); 541 PetscCall(ISGetIndices(is, &points)); 542 for (PetscInt i = 0; i < npoints; i++) { 543 PetscInt point = points[i], off, dof; 544 PetscCall(PetscSectionGetOffset(section, point, &off)); 545 PetscCall(PetscSectionGetDof(section, point, &dof)); 546 PetscAssert(dof % dim == 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Unexpected dof %" PetscInt_FMT " not divisible by dimension %" PetscInt_FMT, dof, dim); 547 for (PetscInt j = 0; j < dof / dim; j++) { 548 PetscInt *slot; 549 PetscCall(PetscSegBufferGetInts(seg, 1, &slot)); 550 *slot = off / dim; 551 } 552 } 553 PetscInt *ind; 554 PetscCall(PetscSegBufferGetSize(seg, &count)); 555 PetscCall(PetscSegBufferExtractAlloc(seg, &ind)); 556 PetscCall(PetscSegBufferDestroy(&seg)); 557 PetscCall(ISCreateBlock(PETSC_COMM_SELF, dim, count, ind, PETSC_OWN_POINTER, &isdof)); 558 } 559 Vec L, P; 560 VecType vec_type; 561 VecScatter scatter; 562 PetscCall(DMGetLocalVector(dm, &L)); 563 PetscCall(VecCreate(PETSC_COMM_SELF, &P)); 564 PetscCall(VecSetSizes(P, count * dim, count * dim)); 565 PetscCall(VecGetType(L, &vec_type)); 566 PetscCall(VecSetType(P, vec_type)); 567 PetscCall(VecScatterCreate(P, NULL, L, isdof, &scatter)); 568 PetscCall(DMRestoreLocalVector(dm, &L)); 569 PetscCall(ISDestroy(&isdof)); 570 571 { 572 PetscScalar *x; 573 PetscCall(VecGetArrayWrite(P, &x)); 574 for (PetscInt i = 0; i < (PetscInt)count; i++) { 575 for (PetscInt j = 0; j < dim; j++) x[i * dim + j] = plex->periodic.transform[j][3]; 576 } 577 PetscCall(VecRestoreArrayWrite(P, &x)); 578 } 579 580 dm->periodic.affine_to_local = scatter; 581 dm->periodic.affine = P; 582 PetscCall(DMGlobalToLocalHookAdd(dm, NULL, DMCoordAddPeriodicOffsets_Private, NULL)); 583 PetscFunctionReturn(0); 584 } 585 586 // We'll just orient all the edges, though only periodic boundary edges need orientation 587 static PetscErrorCode DMPlexOrientPositiveEdges_Private(DM dm) 588 { 589 PetscInt dim, eStart, eEnd; 590 PetscFunctionBegin; 591 PetscCall(DMGetDimension(dm, &dim)); 592 if (dim < 3) PetscFunctionReturn(0); // not necessary 593 PetscCall(DMPlexGetDepthStratum(dm, 1, &eStart, &eEnd)); 594 for (PetscInt e = eStart; e < eEnd; e++) { 595 const PetscInt *cone; 596 PetscCall(DMPlexGetCone(dm, e, &cone)); 597 if (cone[0] > cone[1]) PetscCall(DMPlexOrientPoint(dm, e, -1)); 598 } 599 PetscFunctionReturn(0); 600 } 601 602 PetscErrorCode DMPlexCreateBoxMesh_Tensor_SFC_Internal(DM dm, PetscInt dim, const PetscInt faces[], const PetscReal lower[], const PetscReal upper[], const DMBoundaryType periodicity[], PetscBool interpolate) 603 { 604 PetscInt eextent[3] = {1, 1, 1}, vextent[3] = {1, 1, 1}; 605 const Ijk closure_1[] = { 606 {0, 0, 0}, 607 {1, 0, 0}, 608 }; 609 const Ijk closure_2[] = { 610 {0, 0, 0}, 611 {1, 0, 0}, 612 {1, 1, 0}, 613 {0, 1, 0}, 614 }; 615 const Ijk closure_3[] = { 616 {0, 0, 0}, 617 {0, 1, 0}, 618 {1, 1, 0}, 619 {1, 0, 0}, 620 {0, 0, 1}, 621 {1, 0, 1}, 622 {1, 1, 1}, 623 {0, 1, 1}, 624 }; 625 const Ijk *const closure_dim[] = {NULL, closure_1, closure_2, closure_3}; 626 // This must be kept consistent with DMPlexCreateCubeMesh_Internal 627 const PetscInt face_marker_1[] = {1, 2}; 628 const PetscInt face_marker_2[] = {4, 2, 1, 3}; 629 const PetscInt face_marker_3[] = {6, 5, 3, 4, 1, 2}; 630 const PetscInt *const face_marker_dim[] = {NULL, face_marker_1, face_marker_2, face_marker_3}; 631 // Orient faces so the normal is in the positive axis and the first vertex is the one closest to zero. 632 // These orientations can be determined by examining cones of a reference quad and hex element. 633 const PetscInt face_orient_1[] = {0, 0}; 634 const PetscInt face_orient_2[] = {-1, 0, 0, -1}; 635 const PetscInt face_orient_3[] = {-2, 0, -2, 1, -2, 0}; 636 const PetscInt *const face_orient_dim[] = {NULL, face_orient_1, face_orient_2, face_orient_3}; 637 638 PetscFunctionBegin; 639 PetscValidPointer(dm, 1); 640 PetscValidLogicalCollectiveInt(dm, dim, 2); 641 PetscCall(DMSetDimension(dm, dim)); 642 PetscMPIInt rank, size; 643 PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)dm), &size)); 644 PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank)); 645 for (PetscInt i = 0; i < dim; i++) { 646 eextent[i] = faces[i]; 647 vextent[i] = faces[i] + 1; 648 } 649 ZLayout layout; 650 PetscCall(ZLayoutCreate(size, eextent, vextent, &layout)); 651 PetscZSet vset; // set of all vertices in the closure of the owned elements 652 PetscCall(PetscZSetCreate(&vset)); 653 PetscInt local_elems = 0; 654 for (ZCode z = layout.zstarts[rank]; z < layout.zstarts[rank + 1]; z++) { 655 Ijk loc = ZCodeSplit(z); 656 if (IjkActive(layout.vextent, loc)) PetscZSetAdd(vset, z); 657 if (IjkActive(layout.eextent, loc)) { 658 local_elems++; 659 // Add all neighboring vertices to set 660 for (PetscInt n = 0; n < PetscPowInt(2, dim); n++) { 661 Ijk inc = closure_dim[dim][n]; 662 Ijk nloc = {loc.i + inc.i, loc.j + inc.j, loc.k + inc.k}; 663 ZCode v = ZEncode(nloc); 664 PetscZSetAdd(vset, v); 665 } 666 } 667 } 668 PetscInt local_verts, off = 0; 669 ZCode *vert_z; 670 PetscCall(PetscZSetGetSize(vset, &local_verts)); 671 PetscCall(PetscMalloc1(local_verts, &vert_z)); 672 PetscCall(PetscZSetGetElems(vset, &off, vert_z)); 673 PetscCall(PetscZSetDestroy(&vset)); 674 // ZCode is unsigned for bitwise convenience, but highest bit should never be set, so can interpret as signed 675 PetscCall(PetscSortInt64(local_verts, (PetscInt64 *)vert_z)); 676 677 PetscCall(DMPlexSetChart(dm, 0, local_elems + local_verts)); 678 for (PetscInt e = 0; e < local_elems; e++) PetscCall(DMPlexSetConeSize(dm, e, PetscPowInt(2, dim))); 679 PetscCall(DMSetUp(dm)); 680 { 681 PetscInt e = 0; 682 for (ZCode z = layout.zstarts[rank]; z < layout.zstarts[rank + 1]; z++) { 683 Ijk loc = ZCodeSplit(z); 684 if (!IjkActive(layout.eextent, loc)) continue; 685 PetscInt cone[8], orient[8] = {0}; 686 for (PetscInt n = 0; n < PetscPowInt(2, dim); n++) { 687 Ijk inc = closure_dim[dim][n]; 688 Ijk nloc = {loc.i + inc.i, loc.j + inc.j, loc.k + inc.k}; 689 ZCode v = ZEncode(nloc); 690 PetscInt ci = ZCodeFind(v, local_verts, vert_z); 691 PetscAssert(ci >= 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Did not find neighbor vertex in set"); 692 cone[n] = local_elems + ci; 693 } 694 PetscCall(DMPlexSetCone(dm, e, cone)); 695 PetscCall(DMPlexSetConeOrientation(dm, e, orient)); 696 e++; 697 } 698 } 699 700 PetscCall(DMPlexSymmetrize(dm)); 701 PetscCall(DMPlexStratify(dm)); 702 703 { // Create point SF 704 PetscSF sf; 705 PetscSFCreate(PetscObjectComm((PetscObject)dm), &sf); 706 PetscInt owned_verts = ZCodeFind(layout.zstarts[rank + 1], local_verts, vert_z); 707 if (owned_verts < 0) owned_verts = -(owned_verts + 1); // We don't care whether the key was found 708 PetscInt num_ghosts = local_verts - owned_verts; // Due to sorting, owned vertices always come first 709 PetscInt *local_ghosts; 710 PetscSFNode *ghosts; 711 PetscCall(PetscMalloc1(num_ghosts, &local_ghosts)); 712 PetscCall(PetscMalloc1(num_ghosts, &ghosts)); 713 for (PetscInt i = 0; i < num_ghosts;) { 714 ZCode z = vert_z[owned_verts + i]; 715 PetscInt remote_rank = ZCodeFind(z, size + 1, layout.zstarts), remote_count = 0; 716 if (remote_rank < 0) remote_rank = -(remote_rank + 1) - 1; 717 // We have a new remote rank; find all the ghost indices (which are contiguous in vert_z) 718 719 // Count the elements on remote_rank 720 PetscInt remote_elem = ZLayoutElementsOnRank(&layout, remote_rank); 721 722 // Traverse vertices and make ghost links 723 for (ZCode rz = layout.zstarts[remote_rank]; rz < layout.zstarts[remote_rank + 1]; rz++) { 724 Ijk loc = ZCodeSplit(rz); 725 if (rz == z) { 726 local_ghosts[i] = local_elems + owned_verts + i; 727 ghosts[i].rank = remote_rank; 728 ghosts[i].index = remote_elem + remote_count; 729 i++; 730 if (i == num_ghosts) break; 731 z = vert_z[owned_verts + i]; 732 } 733 if (IjkActive(layout.vextent, loc)) remote_count++; 734 } 735 } 736 PetscCall(PetscSFSetGraph(sf, local_elems + local_verts, num_ghosts, local_ghosts, PETSC_OWN_POINTER, ghosts, PETSC_OWN_POINTER)); 737 PetscCall(PetscObjectSetName((PetscObject)sf, "SFC Point SF")); 738 PetscCall(DMSetPointSF(dm, sf)); 739 PetscCall(PetscSFDestroy(&sf)); 740 } 741 { 742 Vec coordinates; 743 PetscScalar *coords; 744 PetscSection coord_section; 745 PetscInt coord_size; 746 PetscCall(DMGetCoordinateSection(dm, &coord_section)); 747 PetscCall(PetscSectionSetNumFields(coord_section, 1)); 748 PetscCall(PetscSectionSetFieldComponents(coord_section, 0, dim)); 749 PetscCall(PetscSectionSetChart(coord_section, local_elems, local_elems + local_verts)); 750 for (PetscInt v = 0; v < local_verts; v++) { 751 PetscInt point = local_elems + v; 752 PetscCall(PetscSectionSetDof(coord_section, point, dim)); 753 PetscCall(PetscSectionSetFieldDof(coord_section, point, 0, dim)); 754 } 755 PetscCall(PetscSectionSetUp(coord_section)); 756 PetscCall(PetscSectionGetStorageSize(coord_section, &coord_size)); 757 PetscCall(VecCreate(PETSC_COMM_SELF, &coordinates)); 758 PetscCall(PetscObjectSetName((PetscObject)coordinates, "coordinates")); 759 PetscCall(VecSetSizes(coordinates, coord_size, PETSC_DETERMINE)); 760 PetscCall(VecSetBlockSize(coordinates, dim)); 761 PetscCall(VecSetType(coordinates, VECSTANDARD)); 762 PetscCall(VecGetArray(coordinates, &coords)); 763 for (PetscInt v = 0; v < local_verts; v++) { 764 Ijk loc = ZCodeSplit(vert_z[v]); 765 coords[v * dim + 0] = lower[0] + loc.i * (upper[0] - lower[0]) / layout.eextent.i; 766 if (dim > 1) coords[v * dim + 1] = lower[1] + loc.j * (upper[1] - lower[1]) / layout.eextent.j; 767 if (dim > 2) coords[v * dim + 2] = lower[2] + loc.k * (upper[2] - lower[2]) / layout.eextent.k; 768 } 769 PetscCall(VecRestoreArray(coordinates, &coords)); 770 PetscCall(DMSetCoordinatesLocal(dm, coordinates)); 771 PetscCall(VecDestroy(&coordinates)); 772 } 773 if (interpolate) { 774 PetscCall(DMPlexInterpolateInPlace_Internal(dm)); 775 // It's currently necessary to orient the donor and periodic edges consistently. An easy way to ensure that is ot 776 // give all edges positive orientation. Since vertices are created in Z-order, all ranks will agree about the 777 // ordering cone[0] < cone[1]. This is overkill and it would be nice to remove this preparation and make 778 // DMPlexCreateIsoperiodicClosureSF_Private() more resilient, so it fixes any inconsistent orientations. That might 779 // be needed in a general CGNS reader, for example. 780 PetscCall(DMPlexOrientPositiveEdges_Private(dm)); 781 782 DMLabel label; 783 PetscCall(DMCreateLabel(dm, "Face Sets")); 784 PetscCall(DMGetLabel(dm, "Face Sets", &label)); 785 PetscSegBuffer per_faces, donor_face_closure, my_donor_faces; 786 PetscCall(PetscSegBufferCreate(sizeof(PetscInt), 64, &per_faces)); 787 PetscCall(PetscSegBufferCreate(sizeof(PetscInt), 64, &my_donor_faces)); 788 PetscCall(PetscSegBufferCreate(sizeof(ZCode), 64 * PetscPowInt(2, dim), &donor_face_closure)); 789 PetscInt fStart, fEnd, vStart, vEnd; 790 PetscCall(DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd)); 791 PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd)); 792 for (PetscInt f = fStart; f < fEnd; f++) { 793 PetscInt npoints, *points = NULL, num_fverts = 0, fverts[8]; 794 PetscCall(DMPlexGetTransitiveClosure(dm, f, PETSC_TRUE, &npoints, &points)); 795 PetscInt bc_count[6] = {0}; 796 for (PetscInt i = 0; i < npoints; i++) { 797 PetscInt p = points[2 * i]; 798 if (p < vStart || vEnd <= p) continue; 799 fverts[num_fverts++] = p; 800 Ijk loc = ZCodeSplit(vert_z[p - vStart]); 801 // Convention here matches DMPlexCreateCubeMesh_Internal 802 bc_count[0] += loc.i == 0; 803 bc_count[1] += loc.i == layout.vextent.i - 1; 804 bc_count[2] += loc.j == 0; 805 bc_count[3] += loc.j == layout.vextent.j - 1; 806 bc_count[4] += loc.k == 0; 807 bc_count[5] += loc.k == layout.vextent.k - 1; 808 } 809 PetscCall(DMPlexRestoreTransitiveClosure(dm, f, PETSC_TRUE, &npoints, &points)); 810 for (PetscInt bc = 0, bc_match = 0; bc < 2 * dim; bc++) { 811 if (bc_count[bc] == PetscPowInt(2, dim - 1)) { 812 PetscCall(DMPlexOrientPoint(dm, f, face_orient_dim[dim][bc])); 813 if (periodicity[bc / 2] == DM_BOUNDARY_PERIODIC) { 814 PetscInt *put; 815 if (bc % 2 == 0) { // donor face; no label 816 PetscCall(PetscSegBufferGet(my_donor_faces, 1, &put)); 817 *put = f; 818 } else { // periodic face 819 PetscCall(PetscSegBufferGet(per_faces, 1, &put)); 820 *put = f; 821 ZCode *zput; 822 PetscCall(PetscSegBufferGet(donor_face_closure, num_fverts, &zput)); 823 for (PetscInt i = 0; i < num_fverts; i++) { 824 Ijk loc = ZCodeSplit(vert_z[fverts[i] - vStart]); 825 switch (bc / 2) { 826 case 0: 827 loc.i = 0; 828 break; 829 case 1: 830 loc.j = 0; 831 break; 832 case 2: 833 loc.k = 0; 834 break; 835 } 836 *zput++ = ZEncode(loc); 837 } 838 } 839 continue; 840 } 841 PetscAssert(bc_match == 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Face matches multiple face sets"); 842 PetscCall(DMLabelSetValue(label, f, face_marker_dim[dim][bc])); 843 bc_match++; 844 } 845 } 846 } 847 // Ensure that the Coordinate DM has our new boundary labels 848 DM cdm; 849 PetscCall(DMGetCoordinateDM(dm, &cdm)); 850 PetscCall(DMCopyLabels(dm, cdm, PETSC_COPY_VALUES, PETSC_FALSE, DM_COPY_LABELS_FAIL)); 851 if (periodicity[0] == DM_BOUNDARY_PERIODIC || (dim > 1 && periodicity[1] == DM_BOUNDARY_PERIODIC) || (dim > 2 && periodicity[2] == DM_BOUNDARY_PERIODIC)) { 852 PetscCall(DMPlexCreateBoxMesh_Tensor_SFC_Periodicity_Private(dm, &layout, vert_z, per_faces, lower, upper, periodicity, donor_face_closure, my_donor_faces)); 853 } 854 PetscCall(PetscSegBufferDestroy(&per_faces)); 855 PetscCall(PetscSegBufferDestroy(&donor_face_closure)); 856 PetscCall(PetscSegBufferDestroy(&my_donor_faces)); 857 } 858 PetscCall(PetscFree(layout.zstarts)); 859 PetscCall(PetscFree(vert_z)); 860 PetscFunctionReturn(0); 861 } 862 863 /*@ 864 DMPlexSetIsoperiodicFaceSF - Express periodicity from an existing mesh 865 866 Logically collective 867 868 Input Parameters: 869 + dm - The `DMPLEX` on which to set periodicity 870 - face_sf - SF in which roots are (owned) donor faces and leaves are faces that must be matched to a (possibly remote) donor face. 871 872 Level: advanced 873 874 Notes: 875 876 One can use `-dm_plex_box_sfc` to use this mode of periodicity, wherein the periodic points are distinct both globally 877 and locally, but are paired when creating a global dof space. 878 879 .seealso: [](chapter_unstructured), `DMPLEX`, `DMGetGlobalSection()`, `DMPlexGetIsoperiodicFaceSF()` 880 @*/ 881 PetscErrorCode DMPlexSetIsoperiodicFaceSF(DM dm, PetscSF face_sf) 882 { 883 DM_Plex *plex = (DM_Plex *)dm->data; 884 PetscFunctionBegin; 885 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 886 PetscCall(PetscObjectReference((PetscObject)face_sf)); 887 PetscCall(PetscSFDestroy(&plex->periodic.face_sf)); 888 plex->periodic.face_sf = face_sf; 889 if (face_sf) PetscCall(PetscObjectComposeFunction((PetscObject)dm, "DMGetIsoperiodicPointSF_C", DMGetIsoperiodicPointSF_Plex)); 890 891 DM cdm = dm->coordinates[0].dm; // Can't DMGetCoordinateDM because it automatically creates one 892 if (cdm) { 893 PetscCall(DMPlexSetIsoperiodicFaceSF(cdm, face_sf)); 894 if (face_sf) cdm->periodic.setup = DMPeriodicCoordinateSetUp_Internal; 895 } 896 PetscFunctionReturn(0); 897 } 898 899 /*@ 900 DMPlexGetIsoperiodicFaceSF - Obtain periodicity for a mesh 901 902 Logically collective 903 904 Input Parameters: 905 . dm - The `DMPLEX` for which to obtain periodic relation 906 907 Output Parameters: 908 . face_sf - SF in which roots are (owned) donor faces and leaves are faces that must be matched to a (possibly remote) donor face. 909 910 Level: advanced 911 912 .seealso: [](chapter_unstructured), `DMPLEX`, `DMGetGlobalSection()`, `DMPlexSetIsoperiodicFaceSF()` 913 @*/ 914 PetscErrorCode DMPlexGetIsoperiodicFaceSF(DM dm, PetscSF *face_sf) 915 { 916 DM_Plex *plex = (DM_Plex *)dm->data; 917 PetscFunctionBegin; 918 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 919 *face_sf = plex->periodic.face_sf; 920 PetscFunctionReturn(0); 921 } 922 923 /*@C 924 DMPlexSetIsoperiodicFaceTransform - set geometric transform from donor to periodic points 925 926 Logically Collective 927 928 Input Arguments: 929 + dm - `DMPlex` that has been configured with `DMPlexSetIsoperiodicFaceSF()` 930 - t - 4x4 affine transformation basis. 931 932 Notes: 933 Affine transforms are 4x4 matrices in which the leading 3x3 block expresses a rotation (or identity for no rotation), 934 the last column contains a translation vector, and the bottom row is all zero except the last entry, which must always 935 be 1. This representation is common in geometric modeling and allows affine transformations to be composed using 936 simple matrix mulitplication. 937 938 Although the interface accepts a general affine transform, only affine translation is supported at present. 939 940 Developer Notes: 941 This interface should be replaced by making BasisTransform public, expanding it to support affine representations, and 942 adding GPU implementations to apply the G2L/L2G transforms. 943 @*/ 944 PetscErrorCode DMPlexSetIsoperiodicFaceTransform(DM dm, const PetscScalar t[]) 945 { 946 DM_Plex *plex = (DM_Plex *)dm->data; 947 PetscFunctionBegin; 948 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 949 for (PetscInt i = 0; i < 4; i++) { 950 for (PetscInt j = 0; j < 4; j++) { 951 PetscCheck(i != j || t[i * 4 + j] == 1., PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Rotated transforms not supported"); 952 plex->periodic.transform[i][j] = t[i * 4 + j]; 953 } 954 } 955 PetscFunctionReturn(0); 956 } 957