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 // ***** Overview of ZCode ******* 21 // The SFC uses integer indexing for each dimension and encodes them into a single integer by interleaving the bits of each index. 22 // This is known as Morton encoding, and is refered to as ZCode in this code. 23 // So for index i having bits [i2,i1,i0], and similar for indexes j and k, the ZCode (Morton number) would be: 24 // [k2,j2,i2,k1,j1,i1,k0,j0,i0] 25 // This encoding allows for easier traversal of the SFC structure (see https://en.wikipedia.org/wiki/Z-order_curve and `ZStepOct()`). 26 // `ZEncode()` is used to go from indices to ZCode, while `ZCodeSplit()` goes from ZCode back to indices. 27 28 // Decodes the leading interleaved index from a ZCode 29 // e.g. [k2,j2,i2,k1,j1,i1,k0,j0,i0] -> [i2,i1,i0] 30 // Magic numbers taken from https://stackoverflow.com/a/18528775/7564988 (translated to octal) 31 static unsigned ZCodeSplit1(ZCode z) 32 { 33 z &= 0111111111111111111111; 34 z = (z | z >> 2) & 0103030303030303030303; 35 z = (z | z >> 4) & 0100170017001700170017; 36 z = (z | z >> 8) & 0370000037700000377; 37 z = (z | z >> 16) & 0370000000000177777; 38 z = (z | z >> 32) & 07777777; 39 return (unsigned)z; 40 } 41 42 // Encodes the leading interleaved index from a ZCode 43 // e.g. [i2,i1,i0] -> [0,0,i2,0,0,i1,0,0,i0] 44 static ZCode ZEncode1(unsigned t) 45 { 46 ZCode z = t; 47 z &= 07777777; 48 z = (z | z << 32) & 0370000000000177777; 49 z = (z | z << 16) & 0370000037700000377; 50 z = (z | z << 8) & 0100170017001700170017; 51 z = (z | z << 4) & 0103030303030303030303; 52 z = (z | z << 2) & 0111111111111111111111; 53 return z; 54 } 55 56 // Decodes i j k indices from a ZCode. 57 // Uses `ZCodeSplit1()` by shifting ZCode so that the leading index is the desired one to decode 58 static Ijk ZCodeSplit(ZCode z) 59 { 60 Ijk c; 61 c.i = ZCodeSplit1(z >> 2); 62 c.j = ZCodeSplit1(z >> 1); 63 c.k = ZCodeSplit1(z >> 0); 64 return c; 65 } 66 67 // Encodes i j k indices to a ZCode. 68 // Uses `ZCodeEncode1()` by shifting resulting ZCode to the appropriate bit placement 69 static ZCode ZEncode(Ijk c) 70 { 71 ZCode z = (ZEncode1((unsigned int)c.i) << 2) | (ZEncode1((unsigned int)c.j) << 1) | ZEncode1((unsigned int)c.k); 72 return z; 73 } 74 75 static PetscBool IjkActive(Ijk extent, Ijk l) 76 { 77 if (l.i < extent.i && l.j < extent.j && l.k < extent.k) return PETSC_TRUE; 78 return PETSC_FALSE; 79 } 80 81 // If z is not the base of an octet (last 3 bits 0), return 0. 82 // 83 // If z is the base of an octet, we recursively grow to the biggest structured octet. This is typically useful when a z 84 // is outside the domain and we wish to skip a (possibly recursively large) octet to find our next interesting point. 85 static ZCode ZStepOct(ZCode z) 86 { 87 if (PetscUnlikely(z == 0)) return 0; // Infinite loop below if z == 0 88 ZCode step = 07; 89 for (; (z & step) == 0; step = (step << 3) | 07) { } 90 return step >> 3; 91 } 92 93 // Since element/vertex box extents are typically not equal powers of 2, Z codes that lie within the domain are not contiguous. 94 static PetscErrorCode ZLayoutCreate(PetscMPIInt size, const PetscInt eextent[3], const PetscInt vextent[3], ZLayout *layout) 95 { 96 PetscFunctionBegin; 97 layout->eextent.i = eextent[0]; 98 layout->eextent.j = eextent[1]; 99 layout->eextent.k = eextent[2]; 100 layout->vextent.i = vextent[0]; 101 layout->vextent.j = vextent[1]; 102 layout->vextent.k = vextent[2]; 103 layout->comm_size = size; 104 layout->zstarts = NULL; 105 PetscCall(PetscMalloc1(size + 1, &layout->zstarts)); 106 107 PetscInt total_elems = eextent[0] * eextent[1] * eextent[2]; 108 ZCode z = 0; 109 layout->zstarts[0] = 0; 110 // This loop traverses all vertices in the global domain, so is worth making fast. We use ZStepBound 111 for (PetscMPIInt r = 0; r < size; r++) { 112 PetscInt elems_needed = (total_elems / size) + (total_elems % size > r), count; 113 for (count = 0; count < elems_needed; z++) { 114 ZCode skip = ZStepOct(z); // optimistically attempt a longer step 115 for (ZCode s = skip;; s >>= 3) { 116 Ijk trial = ZCodeSplit(z + s); 117 if (IjkActive(layout->eextent, trial)) { 118 while (count + s + 1 > (ZCode)elems_needed) s >>= 3; // Shrink the octet 119 count += s + 1; 120 z += s; 121 break; 122 } 123 if (s == 0) { // the whole skip octet is inactive 124 z += skip; 125 break; 126 } 127 } 128 } 129 // Pick up any extra vertices in the Z ordering before the next rank's first owned element. 130 // 131 // This leads to poorly balanced vertices when eextent is a power of 2, since all the fringe vertices end up 132 // on the last rank. A possible solution is to balance the Z-order vertices independently from the cells, which will 133 // result in a lot of element closures being remote. We could finish marking boundary conditions, then do a round of 134 // vertex ownership smoothing (which would reorder and redistribute vertices without touching element distribution). 135 // Another would be to have an analytic ownership criteria for vertices in the fringe veextent - eextent. This would 136 // complicate the job of identifying an owner and its offset. 137 // 138 // The current recommended approach is to let `-dm_distribute 1` (default) resolve vertex ownership. This is 139 // *mandatory* with isoperiodicity (except in special cases) to remove standed vertices from local spaces. Here's 140 // the issue: 141 // 142 // Consider this partition on rank 0 (left) and rank 1. 143 // 144 // 4 -------- 5 -- 14 --10 -- 21 --11 145 // | | | 146 // 7 -- 16 -- 8 | | | 147 // | | 3 ------- 7 ------- 9 148 // | | | | 149 // 4 -------- 6 ------ 10 | | 150 // | | | 6 -- 16 -- 8 151 // | | | 152 // 3 ---11--- 5 --18-- 9 153 // 154 // The periodic face SF looks like 155 // [0] Number of roots=21, leaves=1, remote ranks=1 156 // [0] 16 <- (0,11) 157 // [1] Number of roots=22, leaves=2, remote ranks=2 158 // [1] 14 <- (0,18) 159 // [1] 21 <- (1,16) 160 // 161 // 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 162 // the point SF links to (1,4) and (1,5). Rank 1 learns about the periodic mapping of (1,5) while handling face 163 // (1,14), but never learns that vertex (1,4) has been mapped to (0,3) by face (0,16). 164 // 165 // We can relatively easily inform vertex (1,4) of this mapping, but it stays in rank 1's local space despite not 166 // being in the closure and thus not being contributed to. This would be mostly harmless except that some viewer 167 // routines expect all local points to be somehow significant. It is not easy to analytically remove the (1,4) 168 // vertex because the point SF and isoperiodic face SF would need to be updated to account for removal of the 169 // stranded vertices. 170 for (; z <= ZEncode(layout->vextent); z++) { 171 Ijk loc = ZCodeSplit(z); 172 if (IjkActive(layout->eextent, loc)) break; 173 z += ZStepOct(z); 174 } 175 layout->zstarts[r + 1] = z; 176 } 177 layout->zstarts[size] = ZEncode(layout->vextent); 178 PetscFunctionReturn(PETSC_SUCCESS); 179 } 180 181 static PetscInt ZLayoutElementsOnRank(const ZLayout *layout, PetscMPIInt rank) 182 { 183 PetscInt remote_elem = 0; 184 for (ZCode rz = layout->zstarts[rank]; rz < layout->zstarts[rank + 1]; rz++) { 185 Ijk loc = ZCodeSplit(rz); 186 if (IjkActive(layout->eextent, loc)) remote_elem++; 187 else rz += ZStepOct(rz); 188 } 189 return remote_elem; 190 } 191 192 static PetscInt ZCodeFind(ZCode key, PetscInt n, const ZCode X[]) 193 { 194 PetscInt lo = 0, hi = n; 195 196 if (n == 0) return -1; 197 while (hi - lo > 1) { 198 PetscInt mid = lo + (hi - lo) / 2; 199 if (key < X[mid]) hi = mid; 200 else lo = mid; 201 } 202 return key == X[lo] ? lo : -(lo + (key > X[lo]) + 1); 203 } 204 205 static inline PetscBool IsPointInsideStratum(PetscInt point, PetscInt pStart, PetscInt pEnd) 206 { 207 return (point >= pStart && point < pEnd) ? PETSC_TRUE : PETSC_FALSE; 208 } 209 210 static PetscErrorCode DMPlexCreateBoxMesh_Tensor_SFC_Periodicity_Private(DM dm, const ZLayout *layout, const ZCode *vert_z, PetscSegBuffer per_faces[3], const PetscReal *lower, const PetscReal *upper, const DMBoundaryType *periodicity, PetscSegBuffer donor_face_closure[3], PetscSegBuffer my_donor_faces[3]) 211 { 212 MPI_Comm comm; 213 PetscInt dim, vStart, vEnd; 214 PetscMPIInt size; 215 PetscSF face_sfs[3]; 216 PetscScalar transforms[3][4][4] = {{{0}}}; 217 218 PetscFunctionBegin; 219 PetscCall(PetscObjectGetComm((PetscObject)dm, &comm)); 220 PetscCallMPI(MPI_Comm_size(comm, &size)); 221 PetscCall(DMGetDimension(dm, &dim)); 222 const PetscInt csize = PetscPowInt(2, dim - 1); 223 PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd)); 224 225 PetscInt num_directions = 0; 226 for (PetscInt direction = 0; direction < dim; direction++) { 227 PetscCount num_faces; 228 PetscInt *faces; 229 ZCode *donor_verts, *donor_minz; 230 PetscSFNode *leaf; 231 PetscCount num_multiroots = 0; 232 PetscInt pStart, pEnd; 233 PetscBool sorted; 234 PetscInt inum_faces; 235 236 if (periodicity[direction] != DM_BOUNDARY_PERIODIC) continue; 237 PetscCall(PetscSegBufferGetSize(per_faces[direction], &num_faces)); 238 PetscCall(PetscSegBufferExtractInPlace(per_faces[direction], &faces)); 239 PetscCall(PetscSegBufferExtractInPlace(donor_face_closure[direction], &donor_verts)); 240 PetscCall(PetscMalloc1(num_faces, &donor_minz)); 241 PetscCall(PetscMalloc1(num_faces, &leaf)); 242 for (PetscCount i = 0; i < num_faces; i++) { 243 ZCode minz = donor_verts[i * csize]; 244 245 for (PetscInt j = 1; j < csize; j++) minz = PetscMin(minz, donor_verts[i * csize + j]); 246 donor_minz[i] = minz; 247 } 248 PetscCall(PetscIntCast(num_faces, &inum_faces)); 249 PetscCall(PetscSortedInt64(inum_faces, (const PetscInt64 *)donor_minz, &sorted)); 250 // If a donor vertex were chosen to broker multiple faces, we would have a logic error. 251 // Checking for sorting is a cheap check that there are no duplicates. 252 PetscCheck(sorted, PETSC_COMM_SELF, PETSC_ERR_PLIB, "minz not sorted; possible duplicates not checked"); 253 for (PetscCount i = 0; i < num_faces;) { 254 ZCode z = donor_minz[i]; 255 PetscMPIInt remote_rank, remote_count = 0; 256 257 PetscCall(PetscMPIIntCast(ZCodeFind(z, size + 1, layout->zstarts), &remote_rank)); 258 if (remote_rank < 0) remote_rank = -(remote_rank + 1) - 1; 259 // Process all the vertices on this rank 260 for (ZCode rz = layout->zstarts[remote_rank]; rz < layout->zstarts[remote_rank + 1]; rz++) { 261 Ijk loc = ZCodeSplit(rz); 262 263 if (rz == z) { 264 leaf[i].rank = remote_rank; 265 leaf[i].index = remote_count; 266 i++; 267 if (i == num_faces) break; 268 z = donor_minz[i]; 269 } 270 if (IjkActive(layout->vextent, loc)) remote_count++; 271 } 272 } 273 PetscCall(PetscFree(donor_minz)); 274 PetscCall(PetscSFCreate(PetscObjectComm((PetscObject)dm), &face_sfs[num_directions])); 275 PetscCall(PetscSFSetGraph(face_sfs[num_directions], vEnd - vStart, inum_faces, NULL, PETSC_USE_POINTER, leaf, PETSC_USE_POINTER)); 276 const PetscInt *my_donor_degree; 277 PetscCall(PetscSFComputeDegreeBegin(face_sfs[num_directions], &my_donor_degree)); 278 PetscCall(PetscSFComputeDegreeEnd(face_sfs[num_directions], &my_donor_degree)); 279 280 for (PetscInt i = 0; i < vEnd - vStart; i++) { 281 num_multiroots += my_donor_degree[i]; 282 if (my_donor_degree[i] == 0) continue; 283 PetscAssert(my_donor_degree[i] == 1, PETSC_COMM_SELF, PETSC_ERR_SUP, "Local vertex has multiple faces"); 284 } 285 PetscInt *my_donors, *donor_indices, *my_donor_indices; 286 PetscCount num_my_donors; 287 288 PetscCall(PetscSegBufferGetSize(my_donor_faces[direction], &num_my_donors)); 289 PetscCheck(num_my_donors == num_multiroots, PETSC_COMM_SELF, PETSC_ERR_SUP, "Donor request (%" PetscCount_FMT ") does not match expected donors (%" PetscCount_FMT ")", num_my_donors, num_multiroots); 290 PetscCall(PetscSegBufferExtractInPlace(my_donor_faces[direction], &my_donors)); 291 PetscCall(PetscMalloc1(vEnd - vStart, &my_donor_indices)); 292 for (PetscCount i = 0; i < num_my_donors; i++) { 293 PetscInt f = my_donors[i]; 294 PetscInt num_points, *points = NULL, minv = PETSC_INT_MAX; 295 296 PetscCall(DMPlexGetTransitiveClosure(dm, f, PETSC_TRUE, &num_points, &points)); 297 for (PetscInt j = 0; j < num_points; j++) { 298 PetscInt p = points[2 * j]; 299 if (!IsPointInsideStratum(p, vStart, vEnd)) continue; 300 minv = PetscMin(minv, p); 301 } 302 PetscCall(DMPlexRestoreTransitiveClosure(dm, f, PETSC_TRUE, &num_points, &points)); 303 PetscAssert(my_donor_degree[minv - vStart] == 1, PETSC_COMM_SELF, PETSC_ERR_SUP, "Local vertex not requested"); 304 my_donor_indices[minv - vStart] = f; 305 } 306 PetscCall(PetscMalloc1(num_faces, &donor_indices)); 307 PetscCall(PetscSFBcastBegin(face_sfs[num_directions], MPIU_INT, my_donor_indices, donor_indices, MPI_REPLACE)); 308 PetscCall(PetscSFBcastEnd(face_sfs[num_directions], MPIU_INT, my_donor_indices, donor_indices, MPI_REPLACE)); 309 PetscCall(PetscFree(my_donor_indices)); 310 // Modify our leafs so they point to donor faces instead of donor minz. Additionally, give them indices as faces. 311 for (PetscCount i = 0; i < num_faces; i++) leaf[i].index = donor_indices[i]; 312 PetscCall(PetscFree(donor_indices)); 313 PetscCall(DMPlexGetChart(dm, &pStart, &pEnd)); 314 PetscCall(PetscSFSetGraph(face_sfs[num_directions], pEnd - pStart, inum_faces, faces, PETSC_COPY_VALUES, leaf, PETSC_OWN_POINTER)); 315 { 316 char face_sf_name[PETSC_MAX_PATH_LEN]; 317 PetscCall(PetscSNPrintf(face_sf_name, sizeof face_sf_name, "Z-order Isoperiodic Faces #%" PetscInt_FMT, num_directions)); 318 PetscCall(PetscObjectSetName((PetscObject)face_sfs[num_directions], face_sf_name)); 319 } 320 321 transforms[num_directions][0][0] = 1; 322 transforms[num_directions][1][1] = 1; 323 transforms[num_directions][2][2] = 1; 324 transforms[num_directions][3][3] = 1; 325 transforms[num_directions][direction][3] = upper[direction] - lower[direction]; 326 num_directions++; 327 } 328 329 PetscCall(DMPlexSetIsoperiodicFaceSF(dm, num_directions, face_sfs)); 330 PetscCall(DMPlexSetIsoperiodicFaceTransform(dm, num_directions, (PetscScalar *)transforms)); 331 332 for (PetscInt i = 0; i < num_directions; i++) PetscCall(PetscSFDestroy(&face_sfs[i])); 333 PetscFunctionReturn(PETSC_SUCCESS); 334 } 335 336 // This is a DMGlobalToLocalHook that applies the affine offsets. When extended for rotated periodicity, it'll need to 337 // apply a rotatonal transform and similar operations will be needed for fields (e.g., to rotate a velocity vector). 338 // We use this crude approach here so we don't have to write new GPU kernels yet. 339 static PetscErrorCode DMCoordAddPeriodicOffsets_Private(DM dm, Vec g, InsertMode mode, Vec l, void *ctx) 340 { 341 PetscFunctionBegin; 342 // These `VecScatter`s should be merged to improve efficiency; the scatters cannot be overlapped. 343 for (PetscInt i = 0; i < dm->periodic.num_affines; i++) { 344 PetscCall(VecScatterBegin(dm->periodic.affine_to_local[i], dm->periodic.affine[i], l, ADD_VALUES, SCATTER_FORWARD)); 345 PetscCall(VecScatterEnd(dm->periodic.affine_to_local[i], dm->periodic.affine[i], l, ADD_VALUES, SCATTER_FORWARD)); 346 } 347 PetscFunctionReturn(PETSC_SUCCESS); 348 } 349 350 // Modify Vec based on the transformation of `point` for the given section and field 351 static PetscErrorCode DMPlexOrientFieldPointVec(DM dm, PetscSection section, PetscInt field, Vec V, PetscInt point, PetscInt orientation) 352 { 353 PetscScalar *copy, *V_arr; 354 PetscInt dof, off, point_ornt[2] = {point, orientation}; 355 const PetscInt **perms; 356 const PetscScalar **rots; 357 358 PetscFunctionBeginUser; 359 PetscCall(PetscSectionGetDof(section, point, &dof)); 360 PetscCall(PetscSectionGetOffset(section, point, &off)); 361 PetscCall(VecGetArray(V, &V_arr)); 362 PetscCall(DMGetWorkArray(dm, dof, MPIU_SCALAR, ©)); 363 PetscArraycpy(copy, &V_arr[off], dof); 364 365 PetscCall(PetscSectionGetFieldPointSyms(section, field, 1, point_ornt, &perms, &rots)); 366 for (PetscInt i = 0; i < dof; i++) { 367 if (perms[0]) V_arr[off + perms[0][i]] = copy[i]; 368 if (rots[0]) V_arr[off + perms[0][i]] *= rots[0][i]; 369 } 370 371 PetscCall(PetscSectionRestoreFieldPointSyms(section, field, 1, point_ornt, &perms, &rots)); 372 PetscCall(DMRestoreWorkArray(dm, dof, MPIU_SCALAR, ©)); 373 PetscCall(VecRestoreArray(V, &V_arr)); 374 PetscFunctionReturn(PETSC_SUCCESS); 375 } 376 377 // Reorient the point in the DMPlex while also applying necessary corrections to other structures (e.g. coordinates) 378 static PetscErrorCode DMPlexOrientPointWithCorrections(DM dm, PetscInt point, PetscInt ornt) 379 { 380 // TODO: Potential speed up if we early exit for ornt == 0 (i.e. if ornt is identity, we don't need to do anything) 381 PetscFunctionBeginUser; 382 PetscCall(DMPlexOrientPoint(dm, point, ornt)); 383 384 { // Correct coordinates based on new cone ordering 385 DM cdm; 386 PetscSection csection; 387 Vec coordinates; 388 PetscInt pStart, pEnd; 389 390 PetscCall(DMGetCoordinatesLocal(dm, &coordinates)); 391 PetscCall(DMGetCoordinateDM(dm, &cdm)); 392 PetscCall(DMGetLocalSection(cdm, &csection)); 393 PetscCall(PetscSectionGetChart(csection, &pStart, &pEnd)); 394 if (IsPointInsideStratum(point, pStart, pEnd)) PetscCall(DMPlexOrientFieldPointVec(cdm, csection, 0, coordinates, point, ornt)); 395 } 396 // TODO: Correct sfNatural 397 PetscFunctionReturn(PETSC_SUCCESS); 398 } 399 400 // Creates SF to communicate data from donor to periodic faces. The data can be different sizes per donor-periodic pair and is given in `point_sizes[]` 401 static PetscErrorCode CreateDonorToPeriodicSF(DM dm, PetscSF face_sf, PetscInt pStart, PetscInt pEnd, const PetscInt point_sizes[], PetscInt *rootbuffersize, PetscInt *leafbuffersize, PetscBT *rootbt, PetscSF *sf_closure) 402 { 403 MPI_Comm comm; 404 PetscMPIInt rank; 405 PetscInt nroots, nleaves; 406 PetscInt *rootdata, *leafdata; 407 const PetscInt *filocal; 408 const PetscSFNode *firemote; 409 410 PetscFunctionBeginUser; 411 PetscCall(PetscObjectGetComm((PetscObject)dm, &comm)); 412 PetscCallMPI(MPI_Comm_rank(comm, &rank)); 413 414 PetscCall(PetscSFGetGraph(face_sf, &nroots, &nleaves, &filocal, &firemote)); 415 PetscCall(PetscCalloc2(2 * nroots, &rootdata, 2 * nroots, &leafdata)); 416 for (PetscInt i = 0; i < nleaves; i++) { 417 PetscInt point = filocal[i]; 418 PetscCheck(IsPointInsideStratum(point, pStart, pEnd), PETSC_COMM_SELF, PETSC_ERR_PLIB, "Point %" PetscInt_FMT " in leaves exists outside of stratum [%" PetscInt_FMT ", %" PetscInt_FMT ")", point, pStart, pEnd); 419 leafdata[point] = point_sizes[point - pStart]; 420 } 421 PetscCall(PetscSFReduceBegin(face_sf, MPIU_INT, leafdata, rootdata + nroots, MPIU_SUM)); 422 PetscCall(PetscSFReduceEnd(face_sf, MPIU_INT, leafdata, rootdata + nroots, MPIU_SUM)); 423 424 PetscInt root_offset = 0; 425 PetscCall(PetscBTCreate(nroots, rootbt)); 426 for (PetscInt p = 0; p < nroots; p++) { 427 const PetscInt *donor_dof = rootdata + nroots; 428 if (donor_dof[p] == 0) { 429 rootdata[2 * p] = -1; 430 rootdata[2 * p + 1] = -1; 431 continue; 432 } 433 PetscCall(PetscBTSet(*rootbt, p)); 434 PetscCheck(IsPointInsideStratum(p, pStart, pEnd), PETSC_COMM_SELF, PETSC_ERR_PLIB, "Point %" PetscInt_FMT " in roots exists outside of stratum [%" PetscInt_FMT ", %" PetscInt_FMT ")", p, pStart, pEnd); 435 PetscInt p_size = point_sizes[p - pStart]; 436 PetscCheck(donor_dof[p] == p_size, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Reduced leaf data size (%" PetscInt_FMT ") does not match root data size (%" PetscInt_FMT ")", donor_dof[p], p_size); 437 rootdata[2 * p] = root_offset; 438 rootdata[2 * p + 1] = p_size; 439 root_offset += p_size; 440 } 441 PetscCall(PetscSFBcastBegin(face_sf, MPIU_2INT, rootdata, leafdata, MPI_REPLACE)); 442 PetscCall(PetscSFBcastEnd(face_sf, MPIU_2INT, rootdata, leafdata, MPI_REPLACE)); 443 // Count how many leaves we need to communicate the closures 444 PetscInt leaf_offset = 0; 445 for (PetscInt i = 0; i < nleaves; i++) { 446 PetscInt point = filocal[i]; 447 if (leafdata[2 * point + 1] < 0) continue; 448 leaf_offset += leafdata[2 * point + 1]; 449 } 450 451 PetscSFNode *closure_leaf; 452 PetscCall(PetscMalloc1(leaf_offset, &closure_leaf)); 453 leaf_offset = 0; 454 for (PetscInt i = 0; i < nleaves; i++) { 455 PetscInt point = filocal[i]; 456 PetscInt cl_size = leafdata[2 * point + 1]; 457 if (cl_size < 0) continue; 458 for (PetscInt j = 0; j < cl_size; j++) { 459 closure_leaf[leaf_offset].rank = firemote[i].rank; 460 closure_leaf[leaf_offset].index = leafdata[2 * point] + j; 461 leaf_offset++; 462 } 463 } 464 465 PetscCall(PetscSFCreate(comm, sf_closure)); 466 PetscCall(PetscSFSetGraph(*sf_closure, root_offset, leaf_offset, NULL, PETSC_USE_POINTER, closure_leaf, PETSC_OWN_POINTER)); 467 *rootbuffersize = root_offset; 468 *leafbuffersize = leaf_offset; 469 PetscCall(PetscFree2(rootdata, leafdata)); 470 PetscFunctionReturn(PETSC_SUCCESS); 471 } 472 473 // Determine if `key` is in `array`. `array` does NOT need to be sorted 474 static inline PetscBool SearchIntArray(PetscInt key, PetscInt array_size, const PetscInt array[]) 475 { 476 for (PetscInt i = 0; i < array_size; i++) 477 if (array[i] == key) return PETSC_TRUE; 478 return PETSC_FALSE; 479 } 480 481 // Translate a cone in periodic points to the cone in donor points based on the `periodic2donor` array 482 static inline PetscErrorCode TranslateConeP2D(const PetscInt periodic_cone[], PetscInt cone_size, const PetscInt periodic2donor[], PetscInt p2d_count, PetscInt p2d_cone[]) 483 { 484 PetscFunctionBeginUser; 485 for (PetscInt p = 0; p < cone_size; p++) { 486 PetscInt p2d_index = -1; 487 for (PetscInt p2d = 0; p2d < p2d_count; p2d++) { 488 if (periodic2donor[p2d * 2] == periodic_cone[p]) p2d_index = p2d; 489 } 490 PetscCheck(p2d_index >= 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Could not find periodic point in periodic-to-donor array"); 491 p2d_cone[p] = periodic2donor[2 * p2d_index + 1]; 492 } 493 PetscFunctionReturn(PETSC_SUCCESS); 494 } 495 496 // Corrects the cone order of periodic faces (and their transitive closure's cones) to match their donor face pair. 497 // 498 // This is done by: 499 // 1. Communicating the donor's vertex coordinates and recursive cones (i.e. its own cone and those of it's constituent edges) to it's periodic pairs 500 // - The donor vertices have the isoperiodic transform applied to them such that they should match exactly 501 // 2. Translating the periodic vertices into the donor vertices point IDs 502 // 3. Translating the cone of each periodic point into the donor point IDs 503 // 4. Comparing the periodic-to-donor cone to the donor cone for each point 504 // 5. Apply the necessary transformation to the periodic cone to make it match the donor cone 505 static PetscErrorCode DMPlexCorrectOrientationForIsoperiodic(DM dm) 506 { 507 MPI_Comm comm; 508 DM_Plex *plex = (DM_Plex *)dm->data; 509 PetscInt nroots, nleaves; 510 const PetscInt *filocal; 511 DM cdm; 512 PetscSection csection; 513 Vec coordinates; 514 PetscInt coords_field_id = 0; 515 PetscBool debug_printing = PETSC_FALSE; 516 517 PetscFunctionBeginUser; 518 PetscCall(PetscObjectGetComm((PetscObject)dm, &comm)); 519 PetscCall(DMGetCoordinatesLocal(dm, &coordinates)); 520 PetscCheck(coordinates, comm, PETSC_ERR_ARG_WRONGSTATE, "DM must have coordinates to setup isoperiodic"); 521 PetscCall(DMGetCoordinateDM(dm, &cdm)); 522 PetscCall(DMGetLocalSection(cdm, &csection)); 523 524 for (PetscInt f = 0; f < plex->periodic.num_face_sfs; f++) { 525 PetscSF face_sf = plex->periodic.face_sfs[f]; 526 const PetscScalar(*transform)[4] = (const PetscScalar(*)[4])plex->periodic.transform[f]; 527 PetscInt *face_vertices_size, *face_cones_size; 528 PetscInt fStart, fEnd, vStart, vEnd, rootnumvert, leafnumvert, rootconesize, leafconesize, dim; 529 PetscSF sf_vert_coords, sf_face_cones; 530 PetscBT rootbt; 531 532 PetscCall(DMGetCoordinateDim(dm, &dim)); 533 PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd)); 534 PetscCall(DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd)); 535 PetscCall(PetscCalloc2(fEnd - fStart, &face_vertices_size, fEnd - fStart, &face_cones_size)); 536 537 // Create SFs to communicate donor vertices and donor cones to periodic faces 538 for (PetscInt f = fStart, index = 0; f < fEnd; f++, index++) { 539 PetscInt cl_size, *closure = NULL, num_vertices = 0; 540 PetscCall(DMPlexGetTransitiveClosure(dm, f, PETSC_TRUE, &cl_size, &closure)); 541 for (PetscInt p = 0; p < cl_size; p++) { 542 PetscInt cl_point = closure[2 * p]; 543 if (IsPointInsideStratum(cl_point, vStart, vEnd)) num_vertices++; 544 else { 545 PetscInt cone_size; 546 PetscCall(DMPlexGetConeSize(dm, cl_point, &cone_size)); 547 face_cones_size[index] += cone_size + 2; 548 } 549 } 550 face_vertices_size[index] = num_vertices; 551 face_cones_size[index] += num_vertices; 552 PetscCall(DMPlexRestoreTransitiveClosure(dm, f, PETSC_TRUE, &cl_size, &closure)); 553 } 554 PetscCall(CreateDonorToPeriodicSF(dm, face_sf, fStart, fEnd, face_vertices_size, &rootnumvert, &leafnumvert, &rootbt, &sf_vert_coords)); 555 PetscCall(PetscBTDestroy(&rootbt)); 556 PetscCall(CreateDonorToPeriodicSF(dm, face_sf, fStart, fEnd, face_cones_size, &rootconesize, &leafconesize, &rootbt, &sf_face_cones)); 557 558 PetscCall(PetscSFGetGraph(face_sf, &nroots, &nleaves, &filocal, NULL)); 559 560 PetscReal *leaf_donor_coords; 561 PetscInt *leaf_donor_cones; 562 563 { // Communicate donor coords and cones to the periodic faces 564 PetscReal *mydonor_vertices; 565 PetscInt *mydonor_cones; 566 const PetscScalar *coords_arr; 567 568 PetscCall(PetscCalloc2(rootnumvert * dim, &mydonor_vertices, rootconesize, &mydonor_cones)); 569 PetscCall(VecGetArrayRead(coordinates, &coords_arr)); 570 for (PetscInt donor_face = 0, donor_vert_offset = 0, donor_cone_offset = 0; donor_face < nroots; donor_face++) { 571 if (!PetscBTLookup(rootbt, donor_face)) continue; 572 PetscInt cl_size, *closure = NULL; 573 574 PetscCall(DMPlexGetTransitiveClosure(dm, donor_face, PETSC_TRUE, &cl_size, &closure)); 575 // Pack vertex coordinates 576 for (PetscInt p = 0; p < cl_size; p++) { 577 PetscInt cl_point = closure[2 * p], dof, offset; 578 if (!IsPointInsideStratum(cl_point, vStart, vEnd)) continue; 579 mydonor_cones[donor_cone_offset++] = cl_point; 580 PetscCall(PetscSectionGetFieldDof(csection, cl_point, coords_field_id, &dof)); 581 PetscCall(PetscSectionGetFieldOffset(csection, cl_point, coords_field_id, &offset)); 582 PetscAssert(dof == dim, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Point %" PetscInt_FMT " has dof size %" PetscInt_FMT ", but should match dimension size %" PetscInt_FMT, cl_point, dof, dim); 583 // Apply isoperiodic transform to donor vertices such that corresponding periodic vertices should match exactly 584 for (PetscInt d = 0; d < dof; d++) mydonor_vertices[donor_vert_offset * dim + d] = PetscRealPart(coords_arr[offset + d]) + PetscRealPart(transform[d][3]); 585 donor_vert_offset++; 586 } 587 // Pack cones of face points (including face itself) 588 for (PetscInt p = 0; p < cl_size; p++) { 589 PetscInt cl_point = closure[2 * p], cone_size, depth; 590 const PetscInt *cone; 591 592 PetscCall(DMPlexGetConeSize(dm, cl_point, &cone_size)); 593 PetscCall(DMPlexGetCone(dm, cl_point, &cone)); 594 PetscCall(DMPlexGetPointDepth(dm, cl_point, &depth)); 595 if (depth == 0) continue; // don't include vertex depth 596 mydonor_cones[donor_cone_offset++] = cone_size; 597 mydonor_cones[donor_cone_offset++] = cl_point; 598 PetscArraycpy(&mydonor_cones[donor_cone_offset], cone, cone_size); 599 donor_cone_offset += cone_size; 600 } 601 PetscCall(DMPlexRestoreTransitiveClosure(dm, donor_face, PETSC_TRUE, &cl_size, &closure)); 602 } 603 PetscCall(VecRestoreArrayRead(coordinates, &coords_arr)); 604 PetscCall(PetscBTDestroy(&rootbt)); 605 606 MPI_Datatype vertex_unit; 607 PetscMPIInt n; 608 PetscCall(PetscMPIIntCast(dim, &n)); 609 PetscCallMPI(MPI_Type_contiguous(n, MPIU_REAL, &vertex_unit)); 610 PetscCallMPI(MPI_Type_commit(&vertex_unit)); 611 PetscCall(PetscMalloc2(leafnumvert * 3, &leaf_donor_coords, leafconesize, &leaf_donor_cones)); 612 PetscCall(PetscSFBcastBegin(sf_vert_coords, vertex_unit, mydonor_vertices, leaf_donor_coords, MPI_REPLACE)); 613 PetscCall(PetscSFBcastBegin(sf_face_cones, MPIU_INT, mydonor_cones, leaf_donor_cones, MPI_REPLACE)); 614 PetscCall(PetscSFBcastEnd(sf_vert_coords, vertex_unit, mydonor_vertices, leaf_donor_coords, MPI_REPLACE)); 615 PetscCall(PetscSFBcastEnd(sf_face_cones, MPIU_INT, mydonor_cones, leaf_donor_cones, MPI_REPLACE)); 616 PetscCall(PetscSFDestroy(&sf_vert_coords)); 617 PetscCall(PetscSFDestroy(&sf_face_cones)); 618 PetscCallMPI(MPI_Type_free(&vertex_unit)); 619 PetscCall(PetscFree2(mydonor_vertices, mydonor_cones)); 620 } 621 622 { // Determine periodic orientation w/rt donor vertices and reorient 623 PetscReal tol = PetscSqr(PETSC_MACHINE_EPSILON * 1e3); 624 PetscInt *periodic2donor, dm_depth, maxConeSize; 625 PetscInt coords_offset = 0, cones_offset = 0; 626 627 PetscCall(DMPlexGetDepth(dm, &dm_depth)); 628 PetscCall(DMPlexGetMaxSizes(dm, &maxConeSize, NULL)); 629 PetscCall(DMGetWorkArray(dm, 2 * PetscPowInt(maxConeSize, dm_depth - 1), MPIU_INT, &periodic2donor)); 630 631 // Translate the periodic face vertices into the donor vertices 632 // Translation stored in periodic2donor 633 for (PetscInt i = 0; i < nleaves; i++) { 634 PetscInt periodic_face = filocal[i], cl_size, num_verts = face_vertices_size[periodic_face - fStart]; 635 PetscInt cones_size = face_cones_size[periodic_face - fStart], p2d_count = 0; 636 PetscInt *closure = NULL; 637 638 PetscCall(DMPlexGetTransitiveClosure(dm, periodic_face, PETSC_TRUE, &cl_size, &closure)); 639 for (PetscInt p = 0; p < cl_size; p++) { 640 PetscInt cl_point = closure[2 * p], coords_size, donor_vertex = -1; 641 PetscScalar *coords = NULL; 642 643 if (!IsPointInsideStratum(cl_point, vStart, vEnd)) continue; 644 PetscCall(DMPlexVecGetClosure(dm, csection, coordinates, cl_point, &coords_size, &coords)); 645 PetscAssert(coords_size == dim, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Point %" PetscInt_FMT " has dof size %" PetscInt_FMT ", but should match dimension size %" PetscInt_FMT, cl_point, coords_size, dim); 646 647 for (PetscInt v = 0; v < num_verts; v++) { 648 PetscReal dist_sqr = 0; 649 for (PetscInt d = 0; d < coords_size; d++) dist_sqr += PetscSqr(PetscRealPart(coords[d]) - leaf_donor_coords[(v + coords_offset) * dim + d]); 650 if (dist_sqr < tol) { 651 donor_vertex = leaf_donor_cones[cones_offset + v]; 652 break; 653 } 654 } 655 PetscCheck(donor_vertex >= 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Periodic face %" PetscInt_FMT " could not find matching donor vertex for vertex %" PetscInt_FMT, periodic_face, cl_point); 656 if (PetscDefined(USE_DEBUG)) { 657 for (PetscInt c = 0; c < p2d_count; c++) PetscCheck(periodic2donor[2 * c + 1] != donor_vertex, comm, PETSC_ERR_PLIB, "Found repeated cone_point in periodic_ordering"); 658 } 659 660 periodic2donor[2 * p2d_count + 0] = cl_point; 661 periodic2donor[2 * p2d_count + 1] = donor_vertex; 662 p2d_count++; 663 PetscCall(DMPlexVecRestoreClosure(dm, csection, coordinates, cl_point, &coords_size, &coords)); 664 } 665 coords_offset += num_verts; 666 PetscCall(DMPlexRestoreTransitiveClosure(dm, periodic_face, PETSC_TRUE, &cl_size, &closure)); 667 668 { // Determine periodic orientation w/rt donor vertices and reorient 669 PetscInt depth, *p2d_cone, face_is_array[1] = {periodic_face}; 670 IS *is_arrays, face_is; 671 PetscSection *section_arrays; 672 PetscInt *donor_cone_array = &leaf_donor_cones[cones_offset + num_verts]; 673 674 PetscCall(ISCreateGeneral(PETSC_COMM_SELF, 1, face_is_array, PETSC_USE_POINTER, &face_is)); 675 PetscCall(DMPlexGetConeRecursive(dm, face_is, &depth, &is_arrays, §ion_arrays)); 676 PetscCall(ISDestroy(&face_is)); 677 PetscCall(DMGetWorkArray(dm, maxConeSize, MPIU_INT, &p2d_cone)); 678 for (PetscInt d = 0; d < depth - 1; d++) { 679 PetscInt pStart, pEnd; 680 PetscSection section = section_arrays[d]; 681 const PetscInt *periodic_cone_arrays, *periodic_point_arrays; 682 683 PetscCall(ISGetIndices(is_arrays[d], &periodic_cone_arrays)); 684 PetscCall(ISGetIndices(is_arrays[d + 1], &periodic_point_arrays)); // Points at d+1 correspond to the cones at d 685 PetscCall(PetscSectionGetChart(section_arrays[d], &pStart, &pEnd)); 686 for (PetscInt p = pStart; p < pEnd; p++) { 687 PetscInt periodic_cone_size, periodic_cone_offset, periodic_point = periodic_point_arrays[p]; 688 689 PetscCall(PetscSectionGetDof(section, p, &periodic_cone_size)); 690 PetscCall(PetscSectionGetOffset(section, p, &periodic_cone_offset)); 691 const PetscInt *periodic_cone = &periodic_cone_arrays[periodic_cone_offset]; 692 PetscCall(TranslateConeP2D(periodic_cone, periodic_cone_size, periodic2donor, p2d_count, p2d_cone)); 693 694 // Find the donor cone that matches the periodic point's cone 695 PetscInt donor_cone_offset = 0, donor_point = -1, *donor_cone = NULL; 696 PetscBool cone_matches = PETSC_FALSE; 697 while (donor_cone_offset < cones_size - num_verts) { 698 PetscInt donor_cone_size = donor_cone_array[donor_cone_offset]; 699 donor_point = donor_cone_array[donor_cone_offset + 1]; 700 donor_cone = &donor_cone_array[donor_cone_offset + 2]; 701 702 if (donor_cone_size != periodic_cone_size) goto next_cone; 703 for (PetscInt c = 0; c < periodic_cone_size; c++) { 704 cone_matches = SearchIntArray(donor_cone[c], periodic_cone_size, p2d_cone); 705 if (!cone_matches) goto next_cone; 706 } 707 // Save the found donor cone's point to the translation array. These will be used for higher depth points. 708 // i.e. we save the edge translations for when we look for face cones 709 periodic2donor[2 * p2d_count + 0] = periodic_point; 710 periodic2donor[2 * p2d_count + 1] = donor_point; 711 p2d_count++; 712 break; 713 714 next_cone: 715 donor_cone_offset += donor_cone_size + 2; 716 } 717 PetscCheck(donor_cone_offset < cones_size - num_verts, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Could not find donor cone equivalent to cone of periodic point %" PetscInt_FMT, periodic_point); 718 719 { // Compare the donor cone with the translated periodic cone and reorient 720 PetscInt ornt; 721 DMPolytopeType cell_type; 722 PetscBool found; 723 PetscCall(DMPlexGetCellType(dm, periodic_point, &cell_type)); 724 PetscCall(DMPolytopeMatchOrientation(cell_type, donor_cone, p2d_cone, &ornt, &found)); 725 PetscCheck(found, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Could not find transformation between donor (%" PetscInt_FMT ") and periodic (%" PetscInt_FMT ") cone's", periodic_point, donor_point); 726 if (debug_printing) PetscCall(PetscPrintf(PETSC_COMM_SELF, "Reorienting point %" PetscInt_FMT " by %" PetscInt_FMT "\n", periodic_point, ornt)); 727 PetscCall(DMPlexOrientPointWithCorrections(dm, periodic_point, ornt)); 728 } 729 } 730 PetscCall(ISRestoreIndices(is_arrays[d], &periodic_cone_arrays)); 731 PetscCall(ISRestoreIndices(is_arrays[d + 1], &periodic_point_arrays)); 732 } 733 PetscCall(DMRestoreWorkArray(dm, maxConeSize, MPIU_INT, &p2d_cone)); 734 PetscCall(DMPlexRestoreConeRecursive(dm, face_is, &depth, &is_arrays, §ion_arrays)); 735 } 736 737 PetscCall(DMPlexRestoreTransitiveClosure(dm, periodic_face, PETSC_TRUE, &cl_size, &closure)); 738 cones_offset += cones_size; 739 } 740 PetscCall(DMRestoreWorkArray(dm, 2 * PetscPowInt(maxConeSize, dm_depth - 1), MPIU_INT, &periodic2donor)); 741 } 742 743 PetscCall(PetscFree2(leaf_donor_coords, leaf_donor_cones)); 744 PetscCall(PetscFree2(face_vertices_size, face_cones_size)); 745 } 746 PetscFunctionReturn(PETSC_SUCCESS); 747 } 748 749 // Start with an SF for a positive depth (e.g., faces) and create a new SF for matched closure. 750 // 751 // Output Arguments: 752 // 753 // + closure_sf - augmented point SF (see `DMGetPointSF()`) that includes the faces and all points in its closure. This 754 // can be used to create a global section and section SF. 755 // - is_points - array of index sets for just the points in the closure of `face_sf`. These may be used to apply an affine 756 // transformation to periodic dofs; see DMPeriodicCoordinateSetUp_Internal(). 757 // 758 static PetscErrorCode DMPlexCreateIsoperiodicPointSF_Private(DM dm, PetscInt num_face_sfs, PetscSF *face_sfs, PetscSF *closure_sf, IS **is_points) 759 { 760 MPI_Comm comm; 761 PetscMPIInt rank; 762 PetscSF point_sf; 763 PetscInt nroots, nleaves; 764 const PetscInt *filocal; 765 const PetscSFNode *firemote; 766 767 PetscFunctionBegin; 768 PetscCall(PetscObjectGetComm((PetscObject)dm, &comm)); 769 PetscCallMPI(MPI_Comm_rank(comm, &rank)); 770 PetscCall(DMGetPointSF(dm, &point_sf)); // Point SF has remote points 771 PetscCall(PetscMalloc1(num_face_sfs, is_points)); 772 773 PetscCall(DMPlexCorrectOrientationForIsoperiodic(dm)); 774 775 for (PetscInt f = 0; f < num_face_sfs; f++) { 776 PetscSF face_sf = face_sfs[f]; 777 PetscInt *cl_sizes; 778 PetscInt fStart, fEnd, rootbuffersize, leafbuffersize; 779 PetscSF sf_closure; 780 PetscBT rootbt; 781 782 PetscCall(DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd)); 783 PetscCall(PetscMalloc1(fEnd - fStart, &cl_sizes)); 784 for (PetscInt f = fStart, index = 0; f < fEnd; f++, index++) { 785 PetscInt cl_size, *closure = NULL; 786 PetscCall(DMPlexGetTransitiveClosure(dm, f, PETSC_TRUE, &cl_size, &closure)); 787 cl_sizes[index] = cl_size - 1; 788 PetscCall(DMPlexRestoreTransitiveClosure(dm, f, PETSC_TRUE, &cl_size, &closure)); 789 } 790 791 PetscCall(CreateDonorToPeriodicSF(dm, face_sf, fStart, fEnd, cl_sizes, &rootbuffersize, &leafbuffersize, &rootbt, &sf_closure)); 792 PetscCall(PetscFree(cl_sizes)); 793 PetscCall(PetscSFGetGraph(face_sf, &nroots, &nleaves, &filocal, &firemote)); 794 795 PetscSFNode *leaf_donor_closure; 796 { // Pack root buffer with owner for every point in the root cones 797 PetscSFNode *donor_closure; 798 const PetscInt *pilocal; 799 const PetscSFNode *piremote; 800 PetscInt npoints; 801 802 PetscCall(PetscSFGetGraph(point_sf, NULL, &npoints, &pilocal, &piremote)); 803 PetscCall(PetscCalloc1(rootbuffersize, &donor_closure)); 804 for (PetscInt p = 0, root_offset = 0; p < nroots; p++) { 805 if (!PetscBTLookup(rootbt, p)) continue; 806 PetscInt cl_size, *closure = NULL; 807 PetscCall(DMPlexGetTransitiveClosure(dm, p, PETSC_TRUE, &cl_size, &closure)); 808 for (PetscInt j = 1; j < cl_size; j++) { 809 PetscInt c = closure[2 * j]; 810 if (pilocal) { 811 PetscInt found = -1; 812 if (npoints > 0) PetscCall(PetscFindInt(c, npoints, pilocal, &found)); 813 if (found >= 0) { 814 donor_closure[root_offset++] = piremote[found]; 815 continue; 816 } 817 } 818 // we own c 819 donor_closure[root_offset].rank = rank; 820 donor_closure[root_offset].index = c; 821 root_offset++; 822 } 823 PetscCall(DMPlexRestoreTransitiveClosure(dm, p, PETSC_TRUE, &cl_size, &closure)); 824 } 825 826 PetscCall(PetscMalloc1(leafbuffersize, &leaf_donor_closure)); 827 PetscCall(PetscSFBcastBegin(sf_closure, MPIU_SF_NODE, donor_closure, leaf_donor_closure, MPI_REPLACE)); 828 PetscCall(PetscSFBcastEnd(sf_closure, MPIU_SF_NODE, donor_closure, leaf_donor_closure, MPI_REPLACE)); 829 PetscCall(PetscSFDestroy(&sf_closure)); 830 PetscCall(PetscFree(donor_closure)); 831 } 832 833 PetscSFNode *new_iremote; 834 PetscCall(PetscCalloc1(nroots, &new_iremote)); 835 for (PetscInt i = 0; i < nroots; i++) new_iremote[i].rank = -1; 836 // Walk leaves and match vertices 837 for (PetscInt i = 0, leaf_offset = 0; i < nleaves; i++) { 838 PetscInt point = filocal[i], cl_size; 839 PetscInt *closure = NULL; 840 PetscCall(DMPlexGetTransitiveClosure(dm, point, PETSC_TRUE, &cl_size, &closure)); 841 for (PetscInt j = 1; j < cl_size; j++) { 842 PetscInt c = closure[2 * j]; 843 PetscSFNode lc = leaf_donor_closure[leaf_offset]; 844 // printf("[%d] face %d.%d: %d ?-- (%d,%d)\n", rank, point, j, c, lc.rank, lc.index); 845 if (new_iremote[c].rank == -1) { 846 new_iremote[c] = lc; 847 } 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"); 848 leaf_offset++; 849 } 850 PetscCall(DMPlexRestoreTransitiveClosure(dm, point, PETSC_TRUE, &cl_size, &closure)); 851 } 852 PetscCall(PetscFree(leaf_donor_closure)); 853 854 // Include face points in closure SF 855 for (PetscInt i = 0; i < nleaves; i++) new_iremote[filocal[i]] = firemote[i]; 856 // consolidate leaves 857 PetscInt *leafdata; 858 PetscCall(PetscMalloc1(nroots, &leafdata)); 859 PetscInt num_new_leaves = 0; 860 for (PetscInt i = 0; i < nroots; i++) { 861 if (new_iremote[i].rank == -1) continue; 862 new_iremote[num_new_leaves] = new_iremote[i]; 863 leafdata[num_new_leaves] = i; 864 num_new_leaves++; 865 } 866 PetscCall(ISCreateGeneral(PETSC_COMM_SELF, num_new_leaves, leafdata, PETSC_COPY_VALUES, &(*is_points)[f])); 867 868 PetscSF csf; 869 PetscCall(PetscSFCreate(comm, &csf)); 870 PetscCall(PetscSFSetGraph(csf, nroots, num_new_leaves, leafdata, PETSC_COPY_VALUES, new_iremote, PETSC_COPY_VALUES)); 871 PetscCall(PetscFree(new_iremote)); // copy and delete because new_iremote is longer than it needs to be 872 PetscCall(PetscFree(leafdata)); 873 PetscCall(PetscBTDestroy(&rootbt)); 874 875 PetscInt npoints; 876 PetscCall(PetscSFGetGraph(point_sf, NULL, &npoints, NULL, NULL)); 877 if (npoints < 0) { // empty point_sf 878 *closure_sf = csf; 879 } else { 880 PetscCall(PetscSFMerge(point_sf, csf, closure_sf)); 881 PetscCall(PetscSFDestroy(&csf)); 882 } 883 if (f > 0) PetscCall(PetscSFDestroy(&point_sf)); // Only destroy if point_sf is from previous calls to PetscSFMerge 884 point_sf = *closure_sf; // Use combined point + isoperiodic SF to define point ownership for further face_sf 885 } 886 PetscCall(PetscObjectSetName((PetscObject)*closure_sf, "Composed Periodic Points")); 887 PetscFunctionReturn(PETSC_SUCCESS); 888 } 889 890 static PetscErrorCode DMGetIsoperiodicPointSF_Plex(DM dm, PetscSF *sf) 891 { 892 DM_Plex *plex = (DM_Plex *)dm->data; 893 894 PetscFunctionBegin; 895 if (!plex->periodic.composed_sf) PetscCall(DMPlexCreateIsoperiodicPointSF_Private(dm, plex->periodic.num_face_sfs, plex->periodic.face_sfs, &plex->periodic.composed_sf, &plex->periodic.periodic_points)); 896 if (sf) *sf = plex->periodic.composed_sf; 897 PetscFunctionReturn(PETSC_SUCCESS); 898 } 899 900 PetscErrorCode DMPlexMigrateIsoperiodicFaceSF_Internal(DM old_dm, DM dm, PetscSF sf_migration) 901 { 902 DM_Plex *plex = (DM_Plex *)old_dm->data; 903 PetscSF sf_point, *new_face_sfs; 904 PetscMPIInt rank; 905 906 PetscFunctionBegin; 907 if (!plex->periodic.face_sfs) PetscFunctionReturn(PETSC_SUCCESS); 908 PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank)); 909 PetscCall(DMGetPointSF(dm, &sf_point)); 910 PetscCall(PetscMalloc1(plex->periodic.num_face_sfs, &new_face_sfs)); 911 912 for (PetscInt f = 0; f < plex->periodic.num_face_sfs; f++) { 913 PetscInt old_npoints, new_npoints, old_nleaf, new_nleaf, point_nleaf; 914 PetscSFNode *new_leafdata, *rootdata, *leafdata; 915 const PetscInt *old_local, *point_local; 916 const PetscSFNode *old_remote, *point_remote; 917 918 PetscCall(PetscSFGetGraph(plex->periodic.face_sfs[f], &old_npoints, &old_nleaf, &old_local, &old_remote)); 919 PetscCall(PetscSFGetGraph(sf_migration, NULL, &new_nleaf, NULL, NULL)); 920 PetscCall(PetscSFGetGraph(sf_point, &new_npoints, &point_nleaf, &point_local, &point_remote)); 921 PetscAssert(new_nleaf == new_npoints, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Expected migration leaf space to match new point root space"); 922 PetscCall(PetscMalloc3(old_npoints, &rootdata, old_npoints, &leafdata, new_npoints, &new_leafdata)); 923 924 // Fill new_leafdata with new owners of all points 925 for (PetscInt i = 0; i < new_npoints; i++) { 926 new_leafdata[i].rank = rank; 927 new_leafdata[i].index = i; 928 } 929 for (PetscInt i = 0; i < point_nleaf; i++) { 930 PetscInt j = point_local[i]; 931 new_leafdata[j] = point_remote[i]; 932 } 933 // REPLACE is okay because every leaf agrees about the new owners 934 PetscCall(PetscSFReduceBegin(sf_migration, MPIU_SF_NODE, new_leafdata, rootdata, MPI_REPLACE)); 935 PetscCall(PetscSFReduceEnd(sf_migration, MPIU_SF_NODE, new_leafdata, rootdata, MPI_REPLACE)); 936 // rootdata now contains the new owners 937 938 // Send to leaves of old space 939 for (PetscInt i = 0; i < old_npoints; i++) { 940 leafdata[i].rank = -1; 941 leafdata[i].index = -1; 942 } 943 PetscCall(PetscSFBcastBegin(plex->periodic.face_sfs[f], MPIU_SF_NODE, rootdata, leafdata, MPI_REPLACE)); 944 PetscCall(PetscSFBcastEnd(plex->periodic.face_sfs[f], MPIU_SF_NODE, rootdata, leafdata, MPI_REPLACE)); 945 946 // Send to new leaf space 947 PetscCall(PetscSFBcastBegin(sf_migration, MPIU_SF_NODE, leafdata, new_leafdata, MPI_REPLACE)); 948 PetscCall(PetscSFBcastEnd(sf_migration, MPIU_SF_NODE, leafdata, new_leafdata, MPI_REPLACE)); 949 950 PetscInt nface = 0, *new_local; 951 PetscSFNode *new_remote; 952 for (PetscInt i = 0; i < new_npoints; i++) nface += (new_leafdata[i].rank >= 0); 953 PetscCall(PetscMalloc1(nface, &new_local)); 954 PetscCall(PetscMalloc1(nface, &new_remote)); 955 nface = 0; 956 for (PetscInt i = 0; i < new_npoints; i++) { 957 if (new_leafdata[i].rank == -1) continue; 958 new_local[nface] = i; 959 new_remote[nface] = new_leafdata[i]; 960 nface++; 961 } 962 PetscCall(PetscFree3(rootdata, leafdata, new_leafdata)); 963 PetscCall(PetscSFCreate(PetscObjectComm((PetscObject)dm), &new_face_sfs[f])); 964 PetscCall(PetscSFSetGraph(new_face_sfs[f], new_npoints, nface, new_local, PETSC_OWN_POINTER, new_remote, PETSC_OWN_POINTER)); 965 { 966 char new_face_sf_name[PETSC_MAX_PATH_LEN]; 967 PetscCall(PetscSNPrintf(new_face_sf_name, sizeof new_face_sf_name, "Migrated Isoperiodic Faces #%" PetscInt_FMT, f)); 968 PetscCall(PetscObjectSetName((PetscObject)new_face_sfs[f], new_face_sf_name)); 969 } 970 } 971 972 PetscCall(DMPlexSetIsoperiodicFaceSF(dm, plex->periodic.num_face_sfs, new_face_sfs)); 973 PetscCall(DMPlexSetIsoperiodicFaceTransform(dm, plex->periodic.num_face_sfs, (PetscScalar *)plex->periodic.transform)); 974 for (PetscInt f = 0; f < plex->periodic.num_face_sfs; f++) PetscCall(PetscSFDestroy(&new_face_sfs[f])); 975 PetscCall(PetscFree(new_face_sfs)); 976 PetscFunctionReturn(PETSC_SUCCESS); 977 } 978 979 PetscErrorCode DMPeriodicCoordinateSetUp_Internal(DM dm) 980 { 981 DM_Plex *plex = (DM_Plex *)dm->data; 982 PetscCount count; 983 IS isdof; 984 PetscInt dim; 985 986 PetscFunctionBegin; 987 if (!plex->periodic.face_sfs) PetscFunctionReturn(PETSC_SUCCESS); 988 PetscCheck(plex->periodic.periodic_points, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Isoperiodic PointSF must be created before this function is called"); 989 990 PetscCall(DMGetCoordinateDim(dm, &dim)); 991 dm->periodic.num_affines = plex->periodic.num_face_sfs; 992 PetscCall(PetscFree2(dm->periodic.affine_to_local, dm->periodic.affine)); 993 PetscCall(PetscMalloc2(dm->periodic.num_affines, &dm->periodic.affine_to_local, dm->periodic.num_affines, &dm->periodic.affine)); 994 995 for (PetscInt f = 0; f < plex->periodic.num_face_sfs; f++) { 996 PetscInt npoints, vsize, isize; 997 const PetscInt *points; 998 IS is = plex->periodic.periodic_points[f]; 999 PetscSegBuffer seg; 1000 PetscSection section; 1001 PetscInt *ind; 1002 Vec L, P; 1003 VecType vec_type; 1004 VecScatter scatter; 1005 PetscScalar *x; 1006 1007 PetscCall(DMGetLocalSection(dm, §ion)); 1008 PetscCall(PetscSegBufferCreate(sizeof(PetscInt), 32, &seg)); 1009 PetscCall(ISGetSize(is, &npoints)); 1010 PetscCall(ISGetIndices(is, &points)); 1011 for (PetscInt i = 0; i < npoints; i++) { 1012 PetscInt point = points[i], off, dof; 1013 1014 PetscCall(PetscSectionGetOffset(section, point, &off)); 1015 PetscCall(PetscSectionGetDof(section, point, &dof)); 1016 PetscAssert(dof % dim == 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Unexpected dof %" PetscInt_FMT " not divisible by dimension %" PetscInt_FMT, dof, dim); 1017 for (PetscInt j = 0; j < dof / dim; j++) { 1018 PetscInt *slot; 1019 1020 PetscCall(PetscSegBufferGetInts(seg, 1, &slot)); 1021 *slot = off / dim + j; 1022 } 1023 } 1024 PetscCall(PetscSegBufferGetSize(seg, &count)); 1025 PetscCall(PetscSegBufferExtractAlloc(seg, &ind)); 1026 PetscCall(PetscSegBufferDestroy(&seg)); 1027 PetscCall(PetscIntCast(count, &isize)); 1028 PetscCall(ISCreateBlock(PETSC_COMM_SELF, dim, isize, ind, PETSC_OWN_POINTER, &isdof)); 1029 1030 PetscCall(PetscIntCast(count * dim, &vsize)); 1031 PetscCall(DMGetLocalVector(dm, &L)); 1032 PetscCall(VecCreate(PETSC_COMM_SELF, &P)); 1033 PetscCall(VecSetSizes(P, vsize, vsize)); 1034 PetscCall(VecGetType(L, &vec_type)); 1035 PetscCall(VecSetType(P, vec_type)); 1036 PetscCall(VecScatterCreate(P, NULL, L, isdof, &scatter)); 1037 PetscCall(DMRestoreLocalVector(dm, &L)); 1038 PetscCall(ISDestroy(&isdof)); 1039 1040 PetscCall(VecGetArrayWrite(P, &x)); 1041 for (PetscCount i = 0; i < count; i++) { 1042 for (PetscInt j = 0; j < dim; j++) x[i * dim + j] = plex->periodic.transform[f][j][3]; 1043 } 1044 PetscCall(VecRestoreArrayWrite(P, &x)); 1045 1046 dm->periodic.affine_to_local[f] = scatter; 1047 dm->periodic.affine[f] = P; 1048 } 1049 PetscCall(DMGlobalToLocalHookAdd(dm, NULL, DMCoordAddPeriodicOffsets_Private, NULL)); 1050 PetscFunctionReturn(PETSC_SUCCESS); 1051 } 1052 1053 PetscErrorCode DMPlexCreateBoxMesh_Tensor_SFC_Internal(DM dm, PetscInt dim, const PetscInt faces[], const PetscReal lower[], const PetscReal upper[], const DMBoundaryType periodicity[], PetscBool interpolate) 1054 { 1055 PetscInt eextent[3] = {1, 1, 1}, vextent[3] = {1, 1, 1}; 1056 const Ijk closure_1[] = { 1057 {0, 0, 0}, 1058 {1, 0, 0}, 1059 }; 1060 const Ijk closure_2[] = { 1061 {0, 0, 0}, 1062 {1, 0, 0}, 1063 {1, 1, 0}, 1064 {0, 1, 0}, 1065 }; 1066 const Ijk closure_3[] = { 1067 {0, 0, 0}, 1068 {0, 1, 0}, 1069 {1, 1, 0}, 1070 {1, 0, 0}, 1071 {0, 0, 1}, 1072 {1, 0, 1}, 1073 {1, 1, 1}, 1074 {0, 1, 1}, 1075 }; 1076 const Ijk *const closure_dim[] = {NULL, closure_1, closure_2, closure_3}; 1077 // This must be kept consistent with DMPlexCreateCubeMesh_Internal 1078 const PetscInt face_marker_1[] = {1, 2}; 1079 const PetscInt face_marker_2[] = {4, 2, 1, 3}; 1080 const PetscInt face_marker_3[] = {6, 5, 3, 4, 1, 2}; 1081 const PetscInt *const face_marker_dim[] = {NULL, face_marker_1, face_marker_2, face_marker_3}; 1082 // Orient faces so the normal is in the positive axis and the first vertex is the one closest to zero. 1083 // These orientations can be determined by examining cones of a reference quad and hex element. 1084 const PetscInt face_orient_1[] = {0, 0}; 1085 const PetscInt face_orient_2[] = {-1, 0, 0, -1}; 1086 const PetscInt face_orient_3[] = {-2, 0, -2, 1, -2, 0}; 1087 const PetscInt *const face_orient_dim[] = {NULL, face_orient_1, face_orient_2, face_orient_3}; 1088 1089 PetscFunctionBegin; 1090 PetscCall(PetscLogEventBegin(DMPLEX_CreateBoxSFC, dm, 0, 0, 0)); 1091 PetscAssertPointer(dm, 1); 1092 PetscValidLogicalCollectiveInt(dm, dim, 2); 1093 PetscCall(DMSetDimension(dm, dim)); 1094 PetscMPIInt rank, size; 1095 PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)dm), &size)); 1096 PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank)); 1097 for (PetscInt i = 0; i < dim; i++) { 1098 eextent[i] = faces[i]; 1099 vextent[i] = faces[i] + 1; 1100 } 1101 ZLayout layout; 1102 PetscCall(ZLayoutCreate(size, eextent, vextent, &layout)); 1103 PetscZSet vset; // set of all vertices in the closure of the owned elements 1104 PetscCall(PetscZSetCreate(&vset)); 1105 PetscInt local_elems = 0; 1106 for (ZCode z = layout.zstarts[rank]; z < layout.zstarts[rank + 1]; z++) { 1107 Ijk loc = ZCodeSplit(z); 1108 if (IjkActive(layout.vextent, loc)) PetscCall(PetscZSetAdd(vset, z)); 1109 else { 1110 z += ZStepOct(z); 1111 continue; 1112 } 1113 if (IjkActive(layout.eextent, loc)) { 1114 local_elems++; 1115 // Add all neighboring vertices to set 1116 for (PetscInt n = 0; n < PetscPowInt(2, dim); n++) { 1117 Ijk inc = closure_dim[dim][n]; 1118 Ijk nloc = {loc.i + inc.i, loc.j + inc.j, loc.k + inc.k}; 1119 ZCode v = ZEncode(nloc); 1120 PetscCall(PetscZSetAdd(vset, v)); 1121 } 1122 } 1123 } 1124 PetscInt local_verts, off = 0; 1125 ZCode *vert_z; 1126 PetscCall(PetscZSetGetSize(vset, &local_verts)); 1127 PetscCall(PetscMalloc1(local_verts, &vert_z)); 1128 PetscCall(PetscZSetGetElems(vset, &off, vert_z)); 1129 PetscCall(PetscZSetDestroy(&vset)); 1130 // ZCode is unsigned for bitwise convenience, but highest bit should never be set, so can interpret as signed 1131 PetscCall(PetscSortInt64(local_verts, (PetscInt64 *)vert_z)); 1132 1133 PetscCall(DMPlexSetChart(dm, 0, local_elems + local_verts)); 1134 for (PetscInt e = 0; e < local_elems; e++) PetscCall(DMPlexSetConeSize(dm, e, PetscPowInt(2, dim))); 1135 PetscCall(DMSetUp(dm)); 1136 { 1137 PetscInt e = 0; 1138 for (ZCode z = layout.zstarts[rank]; z < layout.zstarts[rank + 1]; z++) { 1139 Ijk loc = ZCodeSplit(z); 1140 if (!IjkActive(layout.eextent, loc)) { 1141 z += ZStepOct(z); 1142 continue; 1143 } 1144 PetscInt cone[8], orient[8] = {0}; 1145 for (PetscInt n = 0; n < PetscPowInt(2, dim); n++) { 1146 Ijk inc = closure_dim[dim][n]; 1147 Ijk nloc = {loc.i + inc.i, loc.j + inc.j, loc.k + inc.k}; 1148 ZCode v = ZEncode(nloc); 1149 PetscInt ci = ZCodeFind(v, local_verts, vert_z); 1150 PetscAssert(ci >= 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Did not find neighbor vertex in set"); 1151 cone[n] = local_elems + ci; 1152 } 1153 PetscCall(DMPlexSetCone(dm, e, cone)); 1154 PetscCall(DMPlexSetConeOrientation(dm, e, orient)); 1155 e++; 1156 } 1157 } 1158 1159 PetscCall(DMPlexSymmetrize(dm)); 1160 PetscCall(DMPlexStratify(dm)); 1161 1162 { // Create point SF 1163 PetscSF sf; 1164 PetscCall(PetscSFCreate(PetscObjectComm((PetscObject)dm), &sf)); 1165 PetscInt owned_verts = ZCodeFind(layout.zstarts[rank + 1], local_verts, vert_z); 1166 if (owned_verts < 0) owned_verts = -(owned_verts + 1); // We don't care whether the key was found 1167 PetscInt num_ghosts = local_verts - owned_verts; // Due to sorting, owned vertices always come first 1168 PetscInt *local_ghosts; 1169 PetscSFNode *ghosts; 1170 PetscCall(PetscMalloc1(num_ghosts, &local_ghosts)); 1171 PetscCall(PetscMalloc1(num_ghosts, &ghosts)); 1172 for (PetscInt i = 0; i < num_ghosts;) { 1173 ZCode z = vert_z[owned_verts + i]; 1174 PetscMPIInt remote_rank, remote_count = 0; 1175 1176 PetscCall(PetscMPIIntCast(ZCodeFind(z, size + 1, layout.zstarts), &remote_rank)); 1177 if (remote_rank < 0) remote_rank = -(remote_rank + 1) - 1; 1178 // We have a new remote rank; find all the ghost indices (which are contiguous in vert_z) 1179 1180 // Count the elements on remote_rank 1181 PetscInt remote_elem = ZLayoutElementsOnRank(&layout, remote_rank); 1182 1183 // Traverse vertices and make ghost links 1184 for (ZCode rz = layout.zstarts[remote_rank]; rz < layout.zstarts[remote_rank + 1]; rz++) { 1185 Ijk loc = ZCodeSplit(rz); 1186 if (rz == z) { 1187 local_ghosts[i] = local_elems + owned_verts + i; 1188 ghosts[i].rank = remote_rank; 1189 ghosts[i].index = remote_elem + remote_count; 1190 i++; 1191 if (i == num_ghosts) break; 1192 z = vert_z[owned_verts + i]; 1193 } 1194 if (IjkActive(layout.vextent, loc)) remote_count++; 1195 else rz += ZStepOct(rz); 1196 } 1197 } 1198 PetscCall(PetscSFSetGraph(sf, local_elems + local_verts, num_ghosts, local_ghosts, PETSC_OWN_POINTER, ghosts, PETSC_OWN_POINTER)); 1199 PetscCall(PetscObjectSetName((PetscObject)sf, "SFC Point SF")); 1200 PetscCall(DMSetPointSF(dm, sf)); 1201 PetscCall(PetscSFDestroy(&sf)); 1202 } 1203 { 1204 Vec coordinates; 1205 PetscScalar *coords; 1206 PetscSection coord_section; 1207 PetscInt coord_size; 1208 PetscCall(DMGetCoordinateSection(dm, &coord_section)); 1209 PetscCall(PetscSectionSetNumFields(coord_section, 1)); 1210 PetscCall(PetscSectionSetFieldComponents(coord_section, 0, dim)); 1211 PetscCall(PetscSectionSetChart(coord_section, local_elems, local_elems + local_verts)); 1212 for (PetscInt v = 0; v < local_verts; v++) { 1213 PetscInt point = local_elems + v; 1214 PetscCall(PetscSectionSetDof(coord_section, point, dim)); 1215 PetscCall(PetscSectionSetFieldDof(coord_section, point, 0, dim)); 1216 } 1217 PetscCall(PetscSectionSetUp(coord_section)); 1218 PetscCall(PetscSectionGetStorageSize(coord_section, &coord_size)); 1219 PetscCall(VecCreate(PETSC_COMM_SELF, &coordinates)); 1220 PetscCall(PetscObjectSetName((PetscObject)coordinates, "coordinates")); 1221 PetscCall(VecSetSizes(coordinates, coord_size, PETSC_DETERMINE)); 1222 PetscCall(VecSetBlockSize(coordinates, dim)); 1223 PetscCall(VecSetType(coordinates, VECSTANDARD)); 1224 PetscCall(VecGetArray(coordinates, &coords)); 1225 for (PetscInt v = 0; v < local_verts; v++) { 1226 Ijk loc = ZCodeSplit(vert_z[v]); 1227 coords[v * dim + 0] = lower[0] + loc.i * (upper[0] - lower[0]) / layout.eextent.i; 1228 if (dim > 1) coords[v * dim + 1] = lower[1] + loc.j * (upper[1] - lower[1]) / layout.eextent.j; 1229 if (dim > 2) coords[v * dim + 2] = lower[2] + loc.k * (upper[2] - lower[2]) / layout.eextent.k; 1230 } 1231 PetscCall(VecRestoreArray(coordinates, &coords)); 1232 PetscCall(DMSetCoordinatesLocal(dm, coordinates)); 1233 PetscCall(VecDestroy(&coordinates)); 1234 } 1235 if (interpolate) { 1236 PetscCall(DMPlexInterpolateInPlace_Internal(dm)); 1237 1238 DMLabel label; 1239 PetscCall(DMCreateLabel(dm, "Face Sets")); 1240 PetscCall(DMGetLabel(dm, "Face Sets", &label)); 1241 PetscSegBuffer per_faces[3], donor_face_closure[3], my_donor_faces[3]; 1242 for (PetscInt i = 0; i < 3; i++) { 1243 PetscCall(PetscSegBufferCreate(sizeof(PetscInt), 64, &per_faces[i])); 1244 PetscCall(PetscSegBufferCreate(sizeof(PetscInt), 64, &my_donor_faces[i])); 1245 PetscCall(PetscSegBufferCreate(sizeof(ZCode), 64 * PetscPowInt(2, dim), &donor_face_closure[i])); 1246 } 1247 PetscInt fStart, fEnd, vStart, vEnd; 1248 PetscCall(DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd)); 1249 PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd)); 1250 for (PetscInt f = fStart; f < fEnd; f++) { 1251 PetscInt npoints, *points = NULL, num_fverts = 0, fverts[8]; 1252 PetscCall(DMPlexGetTransitiveClosure(dm, f, PETSC_TRUE, &npoints, &points)); 1253 PetscInt bc_count[6] = {0}; 1254 for (PetscInt i = 0; i < npoints; i++) { 1255 PetscInt p = points[2 * i]; 1256 if (!IsPointInsideStratum(p, vStart, vEnd)) continue; 1257 fverts[num_fverts++] = p; 1258 Ijk loc = ZCodeSplit(vert_z[p - vStart]); 1259 // Convention here matches DMPlexCreateCubeMesh_Internal 1260 bc_count[0] += loc.i == 0; 1261 bc_count[1] += loc.i == layout.vextent.i - 1; 1262 bc_count[2] += loc.j == 0; 1263 bc_count[3] += loc.j == layout.vextent.j - 1; 1264 bc_count[4] += loc.k == 0; 1265 bc_count[5] += loc.k == layout.vextent.k - 1; 1266 } 1267 PetscCall(DMPlexRestoreTransitiveClosure(dm, f, PETSC_TRUE, &npoints, &points)); 1268 for (PetscInt bc = 0, bc_match = 0; bc < 2 * dim; bc++) { 1269 if (bc_count[bc] == PetscPowInt(2, dim - 1)) { 1270 PetscCall(DMPlexOrientPoint(dm, f, face_orient_dim[dim][bc])); 1271 if (periodicity[bc / 2] == DM_BOUNDARY_PERIODIC) { 1272 PetscInt *put; 1273 if (bc % 2 == 0) { // donor face; no label 1274 PetscCall(PetscSegBufferGet(my_donor_faces[bc / 2], 1, &put)); 1275 *put = f; 1276 } else { // periodic face 1277 PetscCall(PetscSegBufferGet(per_faces[bc / 2], 1, &put)); 1278 *put = f; 1279 ZCode *zput; 1280 PetscCall(PetscSegBufferGet(donor_face_closure[bc / 2], num_fverts, &zput)); 1281 for (PetscInt i = 0; i < num_fverts; i++) { 1282 Ijk loc = ZCodeSplit(vert_z[fverts[i] - vStart]); 1283 switch (bc / 2) { 1284 case 0: 1285 loc.i = 0; 1286 break; 1287 case 1: 1288 loc.j = 0; 1289 break; 1290 case 2: 1291 loc.k = 0; 1292 break; 1293 } 1294 *zput++ = ZEncode(loc); 1295 } 1296 } 1297 continue; 1298 } 1299 PetscAssert(bc_match == 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Face matches multiple face sets"); 1300 PetscCall(DMLabelSetValue(label, f, face_marker_dim[dim][bc])); 1301 bc_match++; 1302 } 1303 } 1304 } 1305 // Ensure that the Coordinate DM has our new boundary labels 1306 DM cdm; 1307 PetscCall(DMGetCoordinateDM(dm, &cdm)); 1308 PetscCall(DMCopyLabels(dm, cdm, PETSC_COPY_VALUES, PETSC_FALSE, DM_COPY_LABELS_FAIL)); 1309 if (periodicity[0] == DM_BOUNDARY_PERIODIC || (dim > 1 && periodicity[1] == DM_BOUNDARY_PERIODIC) || (dim > 2 && periodicity[2] == DM_BOUNDARY_PERIODIC)) { 1310 PetscCall(DMPlexCreateBoxMesh_Tensor_SFC_Periodicity_Private(dm, &layout, vert_z, per_faces, lower, upper, periodicity, donor_face_closure, my_donor_faces)); 1311 } 1312 for (PetscInt i = 0; i < 3; i++) { 1313 PetscCall(PetscSegBufferDestroy(&per_faces[i])); 1314 PetscCall(PetscSegBufferDestroy(&donor_face_closure[i])); 1315 PetscCall(PetscSegBufferDestroy(&my_donor_faces[i])); 1316 } 1317 } 1318 PetscCall(PetscFree(layout.zstarts)); 1319 PetscCall(PetscFree(vert_z)); 1320 PetscCall(PetscLogEventEnd(DMPLEX_CreateBoxSFC, dm, 0, 0, 0)); 1321 PetscFunctionReturn(PETSC_SUCCESS); 1322 } 1323 1324 /*@ 1325 DMPlexSetIsoperiodicFaceSF - Express periodicity from an existing mesh 1326 1327 Logically Collective 1328 1329 Input Parameters: 1330 + dm - The `DMPLEX` on which to set periodicity 1331 . num_face_sfs - Number of `PetscSF`s in `face_sfs` 1332 - face_sfs - Array of `PetscSF` in which roots are (owned) donor faces and leaves are faces that must be matched to a (possibly remote) donor face. 1333 1334 Level: advanced 1335 1336 Note: 1337 One can use `-dm_plex_shape zbox` to use this mode of periodicity, wherein the periodic points are distinct both globally 1338 and locally, but are paired when creating a global dof space. 1339 1340 .seealso: [](ch_unstructured), `DMPLEX`, `DMGetGlobalSection()`, `DMPlexGetIsoperiodicFaceSF()` 1341 @*/ 1342 PetscErrorCode DMPlexSetIsoperiodicFaceSF(DM dm, PetscInt num_face_sfs, PetscSF *face_sfs) 1343 { 1344 DM_Plex *plex = (DM_Plex *)dm->data; 1345 1346 PetscFunctionBegin; 1347 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 1348 if (num_face_sfs) PetscCall(PetscObjectComposeFunction((PetscObject)dm, "DMGetIsoperiodicPointSF_C", DMGetIsoperiodicPointSF_Plex)); 1349 else PetscCall(PetscObjectComposeFunction((PetscObject)dm, "DMGetIsoperiodicPointSF_C", NULL)); 1350 if (num_face_sfs == plex->periodic.num_face_sfs && (num_face_sfs == 0 || face_sfs == plex->periodic.face_sfs)) PetscFunctionReturn(PETSC_SUCCESS); 1351 PetscCall(DMSetGlobalSection(dm, NULL)); 1352 1353 for (PetscInt i = 0; i < num_face_sfs; i++) PetscCall(PetscObjectReference((PetscObject)face_sfs[i])); 1354 1355 if (plex->periodic.num_face_sfs > 0) { 1356 for (PetscInt i = 0; i < plex->periodic.num_face_sfs; i++) PetscCall(PetscSFDestroy(&plex->periodic.face_sfs[i])); 1357 PetscCall(PetscFree(plex->periodic.face_sfs)); 1358 } 1359 1360 plex->periodic.num_face_sfs = num_face_sfs; 1361 PetscCall(PetscCalloc1(num_face_sfs, &plex->periodic.face_sfs)); 1362 for (PetscInt i = 0; i < num_face_sfs; i++) plex->periodic.face_sfs[i] = face_sfs[i]; 1363 1364 DM cdm = dm->coordinates[0].dm; // Can't DMGetCoordinateDM because it automatically creates one 1365 if (cdm) { 1366 PetscCall(DMPlexSetIsoperiodicFaceSF(cdm, num_face_sfs, face_sfs)); 1367 if (face_sfs) cdm->periodic.setup = DMPeriodicCoordinateSetUp_Internal; 1368 } 1369 PetscFunctionReturn(PETSC_SUCCESS); 1370 } 1371 1372 /*@C 1373 DMPlexGetIsoperiodicFaceSF - Obtain periodicity for a mesh 1374 1375 Logically Collective 1376 1377 Input Parameter: 1378 . dm - The `DMPLEX` for which to obtain periodic relation 1379 1380 Output Parameters: 1381 + num_face_sfs - Number of `PetscSF`s in the array 1382 - face_sfs - Array of `PetscSF` in which roots are (owned) donor faces and leaves are faces that must be matched to a (possibly remote) donor face. 1383 1384 Level: advanced 1385 1386 .seealso: [](ch_unstructured), `DMPLEX`, `DMGetGlobalSection()`, `DMPlexSetIsoperiodicFaceSF()` 1387 @*/ 1388 PetscErrorCode DMPlexGetIsoperiodicFaceSF(DM dm, PetscInt *num_face_sfs, const PetscSF **face_sfs) 1389 { 1390 DM_Plex *plex = (DM_Plex *)dm->data; 1391 1392 PetscFunctionBegin; 1393 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 1394 *face_sfs = plex->periodic.face_sfs; 1395 *num_face_sfs = plex->periodic.num_face_sfs; 1396 PetscFunctionReturn(PETSC_SUCCESS); 1397 } 1398 1399 /*@C 1400 DMPlexSetIsoperiodicFaceTransform - set geometric transform from donor to periodic points 1401 1402 Logically Collective 1403 1404 Input Parameters: 1405 + dm - `DMPLEX` that has been configured with `DMPlexSetIsoperiodicFaceSF()` 1406 . n - Number of transforms in array 1407 - t - Array of 4x4 affine transformation basis. 1408 1409 Level: advanced 1410 1411 Notes: 1412 Affine transforms are 4x4 matrices in which the leading 3x3 block expresses a rotation (or identity for no rotation), 1413 the last column contains a translation vector, and the bottom row is all zero except the last entry, which must always 1414 be 1. This representation is common in geometric modeling and allows affine transformations to be composed using 1415 simple matrix multiplication. 1416 1417 Although the interface accepts a general affine transform, only affine translation is supported at present. 1418 1419 Developer Notes: 1420 This interface should be replaced by making BasisTransform public, expanding it to support affine representations, and 1421 adding GPU implementations to apply the G2L/L2G transforms. 1422 1423 .seealso: [](ch_unstructured), `DMPLEX`, `DMGetGlobalSection()`, `DMPlexSetIsoperiodicFaceSF()` 1424 @*/ 1425 PetscErrorCode DMPlexSetIsoperiodicFaceTransform(DM dm, PetscInt n, const PetscScalar t[]) 1426 { 1427 DM_Plex *plex = (DM_Plex *)dm->data; 1428 1429 PetscFunctionBegin; 1430 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 1431 PetscCheck(n == plex->periodic.num_face_sfs, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Number of transforms (%" PetscInt_FMT ") must equal number of isoperiodc face SFs (%" PetscInt_FMT ")", n, plex->periodic.num_face_sfs); 1432 1433 PetscCall(PetscFree(plex->periodic.transform)); 1434 PetscCall(PetscMalloc1(n, &plex->periodic.transform)); 1435 for (PetscInt i = 0; i < n; i++) { 1436 for (PetscInt j = 0; j < 4; j++) { 1437 for (PetscInt k = 0; k < 4; k++) { 1438 PetscCheck(j != k || t[i * 16 + j * 4 + k] == 1., PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Rotated transforms not supported"); 1439 plex->periodic.transform[i][j][k] = t[i * 16 + j * 4 + k]; 1440 } 1441 } 1442 } 1443 PetscFunctionReturn(PETSC_SUCCESS); 1444 } 1445