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 // Re-set local coordinates (i.e. destroy global coordinates if they were modified 743 PetscCall(DMSetCoordinatesLocal(dm, coordinates)); 744 745 PetscCall(PetscFree2(leaf_donor_coords, leaf_donor_cones)); 746 PetscCall(PetscFree2(face_vertices_size, face_cones_size)); 747 } 748 PetscFunctionReturn(PETSC_SUCCESS); 749 } 750 751 // Start with an SF for a positive depth (e.g., faces) and create a new SF for matched closure. 752 // 753 // Output Arguments: 754 // 755 // + closure_sf - augmented point SF (see `DMGetPointSF()`) that includes the faces and all points in its closure. This 756 // can be used to create a global section and section SF. 757 // - is_points - array of index sets for just the points in the closure of `face_sf`. These may be used to apply an affine 758 // transformation to periodic dofs; see DMPeriodicCoordinateSetUp_Internal(). 759 // 760 static PetscErrorCode DMPlexCreateIsoperiodicPointSF_Private(DM dm, PetscInt num_face_sfs, PetscSF *face_sfs, PetscSF *closure_sf, IS **is_points) 761 { 762 MPI_Comm comm; 763 PetscMPIInt rank; 764 PetscSF point_sf; 765 PetscInt nroots, nleaves; 766 const PetscInt *filocal; 767 const PetscSFNode *firemote; 768 769 PetscFunctionBegin; 770 PetscCall(PetscObjectGetComm((PetscObject)dm, &comm)); 771 PetscCallMPI(MPI_Comm_rank(comm, &rank)); 772 PetscCall(DMGetPointSF(dm, &point_sf)); // Point SF has remote points 773 PetscCall(PetscMalloc1(num_face_sfs, is_points)); 774 775 PetscCall(DMPlexCorrectOrientationForIsoperiodic(dm)); 776 777 for (PetscInt f = 0; f < num_face_sfs; f++) { 778 PetscSF face_sf = face_sfs[f]; 779 PetscInt *cl_sizes; 780 PetscInt fStart, fEnd, rootbuffersize, leafbuffersize; 781 PetscSF sf_closure; 782 PetscBT rootbt; 783 784 PetscCall(DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd)); 785 PetscCall(PetscMalloc1(fEnd - fStart, &cl_sizes)); 786 for (PetscInt f = fStart, index = 0; f < fEnd; f++, index++) { 787 PetscInt cl_size, *closure = NULL; 788 PetscCall(DMPlexGetTransitiveClosure(dm, f, PETSC_TRUE, &cl_size, &closure)); 789 cl_sizes[index] = cl_size - 1; 790 PetscCall(DMPlexRestoreTransitiveClosure(dm, f, PETSC_TRUE, &cl_size, &closure)); 791 } 792 793 PetscCall(CreateDonorToPeriodicSF(dm, face_sf, fStart, fEnd, cl_sizes, &rootbuffersize, &leafbuffersize, &rootbt, &sf_closure)); 794 PetscCall(PetscFree(cl_sizes)); 795 PetscCall(PetscSFGetGraph(face_sf, &nroots, &nleaves, &filocal, &firemote)); 796 797 PetscSFNode *leaf_donor_closure; 798 { // Pack root buffer with owner for every point in the root cones 799 PetscSFNode *donor_closure; 800 const PetscInt *pilocal; 801 const PetscSFNode *piremote; 802 PetscInt npoints; 803 804 PetscCall(PetscSFGetGraph(point_sf, NULL, &npoints, &pilocal, &piremote)); 805 PetscCall(PetscCalloc1(rootbuffersize, &donor_closure)); 806 for (PetscInt p = 0, root_offset = 0; p < nroots; p++) { 807 if (!PetscBTLookup(rootbt, p)) continue; 808 PetscInt cl_size, *closure = NULL; 809 PetscCall(DMPlexGetTransitiveClosure(dm, p, PETSC_TRUE, &cl_size, &closure)); 810 for (PetscInt j = 1; j < cl_size; j++) { 811 PetscInt c = closure[2 * j]; 812 if (pilocal) { 813 PetscInt found = -1; 814 if (npoints > 0) PetscCall(PetscFindInt(c, npoints, pilocal, &found)); 815 if (found >= 0) { 816 donor_closure[root_offset++] = piremote[found]; 817 continue; 818 } 819 } 820 // we own c 821 donor_closure[root_offset].rank = rank; 822 donor_closure[root_offset].index = c; 823 root_offset++; 824 } 825 PetscCall(DMPlexRestoreTransitiveClosure(dm, p, PETSC_TRUE, &cl_size, &closure)); 826 } 827 828 PetscCall(PetscMalloc1(leafbuffersize, &leaf_donor_closure)); 829 PetscCall(PetscSFBcastBegin(sf_closure, MPIU_SF_NODE, donor_closure, leaf_donor_closure, MPI_REPLACE)); 830 PetscCall(PetscSFBcastEnd(sf_closure, MPIU_SF_NODE, donor_closure, leaf_donor_closure, MPI_REPLACE)); 831 PetscCall(PetscSFDestroy(&sf_closure)); 832 PetscCall(PetscFree(donor_closure)); 833 } 834 835 PetscSFNode *new_iremote; 836 PetscCall(PetscCalloc1(nroots, &new_iremote)); 837 for (PetscInt i = 0; i < nroots; i++) new_iremote[i].rank = -1; 838 // Walk leaves and match vertices 839 for (PetscInt i = 0, leaf_offset = 0; i < nleaves; i++) { 840 PetscInt point = filocal[i], cl_size; 841 PetscInt *closure = NULL; 842 PetscCall(DMPlexGetTransitiveClosure(dm, point, PETSC_TRUE, &cl_size, &closure)); 843 for (PetscInt j = 1; j < cl_size; j++) { 844 PetscInt c = closure[2 * j]; 845 PetscSFNode lc = leaf_donor_closure[leaf_offset]; 846 // printf("[%d] face %d.%d: %d ?-- (%d,%d)\n", rank, point, j, c, lc.rank, lc.index); 847 if (new_iremote[c].rank == -1) { 848 new_iremote[c] = lc; 849 } 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"); 850 leaf_offset++; 851 } 852 PetscCall(DMPlexRestoreTransitiveClosure(dm, point, PETSC_TRUE, &cl_size, &closure)); 853 } 854 PetscCall(PetscFree(leaf_donor_closure)); 855 856 // Include face points in closure SF 857 for (PetscInt i = 0; i < nleaves; i++) new_iremote[filocal[i]] = firemote[i]; 858 // consolidate leaves 859 PetscInt *leafdata; 860 PetscCall(PetscMalloc1(nroots, &leafdata)); 861 PetscInt num_new_leaves = 0; 862 for (PetscInt i = 0; i < nroots; i++) { 863 if (new_iremote[i].rank == -1) continue; 864 new_iremote[num_new_leaves] = new_iremote[i]; 865 leafdata[num_new_leaves] = i; 866 num_new_leaves++; 867 } 868 PetscCall(ISCreateGeneral(PETSC_COMM_SELF, num_new_leaves, leafdata, PETSC_COPY_VALUES, &(*is_points)[f])); 869 870 PetscSF csf; 871 PetscCall(PetscSFCreate(comm, &csf)); 872 PetscCall(PetscSFSetGraph(csf, nroots, num_new_leaves, leafdata, PETSC_COPY_VALUES, new_iremote, PETSC_COPY_VALUES)); 873 PetscCall(PetscFree(new_iremote)); // copy and delete because new_iremote is longer than it needs to be 874 PetscCall(PetscFree(leafdata)); 875 PetscCall(PetscBTDestroy(&rootbt)); 876 877 PetscInt npoints; 878 PetscCall(PetscSFGetGraph(point_sf, NULL, &npoints, NULL, NULL)); 879 if (npoints < 0) { // empty point_sf 880 *closure_sf = csf; 881 } else { 882 PetscCall(PetscSFMerge(point_sf, csf, closure_sf)); 883 PetscCall(PetscSFDestroy(&csf)); 884 } 885 if (f > 0) PetscCall(PetscSFDestroy(&point_sf)); // Only destroy if point_sf is from previous calls to PetscSFMerge 886 point_sf = *closure_sf; // Use combined point + isoperiodic SF to define point ownership for further face_sf 887 } 888 PetscCall(PetscObjectSetName((PetscObject)*closure_sf, "Composed Periodic Points")); 889 PetscFunctionReturn(PETSC_SUCCESS); 890 } 891 892 static PetscErrorCode DMGetIsoperiodicPointSF_Plex(DM dm, PetscSF *sf) 893 { 894 DM_Plex *plex = (DM_Plex *)dm->data; 895 896 PetscFunctionBegin; 897 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)); 898 if (sf) *sf = plex->periodic.composed_sf; 899 PetscFunctionReturn(PETSC_SUCCESS); 900 } 901 902 PetscErrorCode DMPlexMigrateIsoperiodicFaceSF_Internal(DM old_dm, DM dm, PetscSF sf_migration) 903 { 904 DM_Plex *plex = (DM_Plex *)old_dm->data; 905 PetscSF sf_point, *new_face_sfs; 906 PetscMPIInt rank; 907 908 PetscFunctionBegin; 909 if (!plex->periodic.face_sfs) PetscFunctionReturn(PETSC_SUCCESS); 910 PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank)); 911 PetscCall(DMGetPointSF(dm, &sf_point)); 912 PetscCall(PetscMalloc1(plex->periodic.num_face_sfs, &new_face_sfs)); 913 914 for (PetscInt f = 0; f < plex->periodic.num_face_sfs; f++) { 915 PetscInt old_npoints, new_npoints, old_nleaf, new_nleaf, point_nleaf; 916 PetscSFNode *new_leafdata, *rootdata, *leafdata; 917 const PetscInt *old_local, *point_local; 918 const PetscSFNode *old_remote, *point_remote; 919 920 PetscCall(PetscSFGetGraph(plex->periodic.face_sfs[f], &old_npoints, &old_nleaf, &old_local, &old_remote)); 921 PetscCall(PetscSFGetGraph(sf_migration, NULL, &new_nleaf, NULL, NULL)); 922 PetscCall(PetscSFGetGraph(sf_point, &new_npoints, &point_nleaf, &point_local, &point_remote)); 923 PetscAssert(new_nleaf == new_npoints, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Expected migration leaf space to match new point root space"); 924 PetscCall(PetscMalloc3(old_npoints, &rootdata, old_npoints, &leafdata, new_npoints, &new_leafdata)); 925 926 // Fill new_leafdata with new owners of all points 927 for (PetscInt i = 0; i < new_npoints; i++) { 928 new_leafdata[i].rank = rank; 929 new_leafdata[i].index = i; 930 } 931 for (PetscInt i = 0; i < point_nleaf; i++) { 932 PetscInt j = point_local[i]; 933 new_leafdata[j] = point_remote[i]; 934 } 935 // REPLACE is okay because every leaf agrees about the new owners 936 PetscCall(PetscSFReduceBegin(sf_migration, MPIU_SF_NODE, new_leafdata, rootdata, MPI_REPLACE)); 937 PetscCall(PetscSFReduceEnd(sf_migration, MPIU_SF_NODE, new_leafdata, rootdata, MPI_REPLACE)); 938 // rootdata now contains the new owners 939 940 // Send to leaves of old space 941 for (PetscInt i = 0; i < old_npoints; i++) { 942 leafdata[i].rank = -1; 943 leafdata[i].index = -1; 944 } 945 PetscCall(PetscSFBcastBegin(plex->periodic.face_sfs[f], MPIU_SF_NODE, rootdata, leafdata, MPI_REPLACE)); 946 PetscCall(PetscSFBcastEnd(plex->periodic.face_sfs[f], MPIU_SF_NODE, rootdata, leafdata, MPI_REPLACE)); 947 948 // Send to new leaf space 949 PetscCall(PetscSFBcastBegin(sf_migration, MPIU_SF_NODE, leafdata, new_leafdata, MPI_REPLACE)); 950 PetscCall(PetscSFBcastEnd(sf_migration, MPIU_SF_NODE, leafdata, new_leafdata, MPI_REPLACE)); 951 952 PetscInt nface = 0, *new_local; 953 PetscSFNode *new_remote; 954 for (PetscInt i = 0; i < new_npoints; i++) nface += (new_leafdata[i].rank >= 0); 955 PetscCall(PetscMalloc1(nface, &new_local)); 956 PetscCall(PetscMalloc1(nface, &new_remote)); 957 nface = 0; 958 for (PetscInt i = 0; i < new_npoints; i++) { 959 if (new_leafdata[i].rank == -1) continue; 960 new_local[nface] = i; 961 new_remote[nface] = new_leafdata[i]; 962 nface++; 963 } 964 PetscCall(PetscFree3(rootdata, leafdata, new_leafdata)); 965 PetscCall(PetscSFCreate(PetscObjectComm((PetscObject)dm), &new_face_sfs[f])); 966 PetscCall(PetscSFSetGraph(new_face_sfs[f], new_npoints, nface, new_local, PETSC_OWN_POINTER, new_remote, PETSC_OWN_POINTER)); 967 { 968 char new_face_sf_name[PETSC_MAX_PATH_LEN]; 969 PetscCall(PetscSNPrintf(new_face_sf_name, sizeof new_face_sf_name, "Migrated Isoperiodic Faces #%" PetscInt_FMT, f)); 970 PetscCall(PetscObjectSetName((PetscObject)new_face_sfs[f], new_face_sf_name)); 971 } 972 } 973 974 PetscCall(DMPlexSetIsoperiodicFaceSF(dm, plex->periodic.num_face_sfs, new_face_sfs)); 975 PetscCall(DMPlexSetIsoperiodicFaceTransform(dm, plex->periodic.num_face_sfs, (PetscScalar *)plex->periodic.transform)); 976 for (PetscInt f = 0; f < plex->periodic.num_face_sfs; f++) PetscCall(PetscSFDestroy(&new_face_sfs[f])); 977 PetscCall(PetscFree(new_face_sfs)); 978 PetscFunctionReturn(PETSC_SUCCESS); 979 } 980 981 PetscErrorCode DMPeriodicCoordinateSetUp_Internal(DM dm) 982 { 983 DM_Plex *plex = (DM_Plex *)dm->data; 984 PetscCount count; 985 IS isdof; 986 PetscInt dim; 987 988 PetscFunctionBegin; 989 if (!plex->periodic.face_sfs) PetscFunctionReturn(PETSC_SUCCESS); 990 PetscCheck(plex->periodic.periodic_points, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Isoperiodic PointSF must be created before this function is called"); 991 992 PetscCall(DMGetCoordinateDim(dm, &dim)); 993 dm->periodic.num_affines = plex->periodic.num_face_sfs; 994 PetscCall(PetscFree2(dm->periodic.affine_to_local, dm->periodic.affine)); 995 PetscCall(PetscMalloc2(dm->periodic.num_affines, &dm->periodic.affine_to_local, dm->periodic.num_affines, &dm->periodic.affine)); 996 997 for (PetscInt f = 0; f < plex->periodic.num_face_sfs; f++) { 998 PetscInt npoints, vsize, isize; 999 const PetscInt *points; 1000 IS is = plex->periodic.periodic_points[f]; 1001 PetscSegBuffer seg; 1002 PetscSection section; 1003 PetscInt *ind; 1004 Vec L, P; 1005 VecType vec_type; 1006 VecScatter scatter; 1007 PetscScalar *x; 1008 1009 PetscCall(DMGetLocalSection(dm, §ion)); 1010 PetscCall(PetscSegBufferCreate(sizeof(PetscInt), 32, &seg)); 1011 PetscCall(ISGetSize(is, &npoints)); 1012 PetscCall(ISGetIndices(is, &points)); 1013 for (PetscInt i = 0; i < npoints; i++) { 1014 PetscInt point = points[i], off, dof; 1015 1016 PetscCall(PetscSectionGetOffset(section, point, &off)); 1017 PetscCall(PetscSectionGetDof(section, point, &dof)); 1018 PetscAssert(dof % dim == 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Unexpected dof %" PetscInt_FMT " not divisible by dimension %" PetscInt_FMT, dof, dim); 1019 for (PetscInt j = 0; j < dof / dim; j++) { 1020 PetscInt *slot; 1021 1022 PetscCall(PetscSegBufferGetInts(seg, 1, &slot)); 1023 *slot = off / dim + j; 1024 } 1025 } 1026 PetscCall(PetscSegBufferGetSize(seg, &count)); 1027 PetscCall(PetscSegBufferExtractAlloc(seg, &ind)); 1028 PetscCall(PetscSegBufferDestroy(&seg)); 1029 PetscCall(PetscIntCast(count, &isize)); 1030 PetscCall(ISCreateBlock(PETSC_COMM_SELF, dim, isize, ind, PETSC_OWN_POINTER, &isdof)); 1031 1032 PetscCall(PetscIntCast(count * dim, &vsize)); 1033 PetscCall(DMGetLocalVector(dm, &L)); 1034 PetscCall(VecCreate(PETSC_COMM_SELF, &P)); 1035 PetscCall(VecSetSizes(P, vsize, vsize)); 1036 PetscCall(VecGetType(L, &vec_type)); 1037 PetscCall(VecSetType(P, vec_type)); 1038 PetscCall(VecScatterCreate(P, NULL, L, isdof, &scatter)); 1039 PetscCall(DMRestoreLocalVector(dm, &L)); 1040 PetscCall(ISDestroy(&isdof)); 1041 1042 PetscCall(VecGetArrayWrite(P, &x)); 1043 for (PetscCount i = 0; i < count; i++) { 1044 for (PetscInt j = 0; j < dim; j++) x[i * dim + j] = plex->periodic.transform[f][j][3]; 1045 } 1046 PetscCall(VecRestoreArrayWrite(P, &x)); 1047 1048 dm->periodic.affine_to_local[f] = scatter; 1049 dm->periodic.affine[f] = P; 1050 } 1051 PetscCall(DMGlobalToLocalHookAdd(dm, NULL, DMCoordAddPeriodicOffsets_Private, NULL)); 1052 PetscFunctionReturn(PETSC_SUCCESS); 1053 } 1054 1055 PetscErrorCode DMPlexCreateBoxMesh_Tensor_SFC_Internal(DM dm, PetscInt dim, const PetscInt faces[], const PetscReal lower[], const PetscReal upper[], const DMBoundaryType periodicity[], PetscBool interpolate) 1056 { 1057 PetscInt eextent[3] = {1, 1, 1}, vextent[3] = {1, 1, 1}; 1058 const Ijk closure_1[] = { 1059 {0, 0, 0}, 1060 {1, 0, 0}, 1061 }; 1062 const Ijk closure_2[] = { 1063 {0, 0, 0}, 1064 {1, 0, 0}, 1065 {1, 1, 0}, 1066 {0, 1, 0}, 1067 }; 1068 const Ijk closure_3[] = { 1069 {0, 0, 0}, 1070 {0, 1, 0}, 1071 {1, 1, 0}, 1072 {1, 0, 0}, 1073 {0, 0, 1}, 1074 {1, 0, 1}, 1075 {1, 1, 1}, 1076 {0, 1, 1}, 1077 }; 1078 const Ijk *const closure_dim[] = {NULL, closure_1, closure_2, closure_3}; 1079 // This must be kept consistent with DMPlexCreateCubeMesh_Internal 1080 const PetscInt face_marker_1[] = {1, 2}; 1081 const PetscInt face_marker_2[] = {4, 2, 1, 3}; 1082 const PetscInt face_marker_3[] = {6, 5, 3, 4, 1, 2}; 1083 const PetscInt *const face_marker_dim[] = {NULL, face_marker_1, face_marker_2, face_marker_3}; 1084 // Orient faces so the normal is in the positive axis and the first vertex is the one closest to zero. 1085 // These orientations can be determined by examining cones of a reference quad and hex element. 1086 const PetscInt face_orient_1[] = {0, 0}; 1087 const PetscInt face_orient_2[] = {-1, 0, 0, -1}; 1088 const PetscInt face_orient_3[] = {-2, 0, -2, 1, -2, 0}; 1089 const PetscInt *const face_orient_dim[] = {NULL, face_orient_1, face_orient_2, face_orient_3}; 1090 1091 PetscFunctionBegin; 1092 PetscCall(PetscLogEventBegin(DMPLEX_CreateBoxSFC, dm, 0, 0, 0)); 1093 PetscAssertPointer(dm, 1); 1094 PetscValidLogicalCollectiveInt(dm, dim, 2); 1095 PetscCall(DMSetDimension(dm, dim)); 1096 PetscMPIInt rank, size; 1097 PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)dm), &size)); 1098 PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank)); 1099 for (PetscInt i = 0; i < dim; i++) { 1100 eextent[i] = faces[i]; 1101 vextent[i] = faces[i] + 1; 1102 } 1103 ZLayout layout; 1104 PetscCall(ZLayoutCreate(size, eextent, vextent, &layout)); 1105 PetscZSet vset; // set of all vertices in the closure of the owned elements 1106 PetscCall(PetscZSetCreate(&vset)); 1107 PetscInt local_elems = 0; 1108 for (ZCode z = layout.zstarts[rank]; z < layout.zstarts[rank + 1]; z++) { 1109 Ijk loc = ZCodeSplit(z); 1110 if (IjkActive(layout.vextent, loc)) PetscCall(PetscZSetAdd(vset, z)); 1111 else { 1112 z += ZStepOct(z); 1113 continue; 1114 } 1115 if (IjkActive(layout.eextent, loc)) { 1116 local_elems++; 1117 // Add all neighboring vertices to set 1118 for (PetscInt n = 0; n < PetscPowInt(2, dim); n++) { 1119 Ijk inc = closure_dim[dim][n]; 1120 Ijk nloc = {loc.i + inc.i, loc.j + inc.j, loc.k + inc.k}; 1121 ZCode v = ZEncode(nloc); 1122 PetscCall(PetscZSetAdd(vset, v)); 1123 } 1124 } 1125 } 1126 PetscInt local_verts, off = 0; 1127 ZCode *vert_z; 1128 PetscCall(PetscZSetGetSize(vset, &local_verts)); 1129 PetscCall(PetscMalloc1(local_verts, &vert_z)); 1130 PetscCall(PetscZSetGetElems(vset, &off, vert_z)); 1131 PetscCall(PetscZSetDestroy(&vset)); 1132 // ZCode is unsigned for bitwise convenience, but highest bit should never be set, so can interpret as signed 1133 PetscCall(PetscSortInt64(local_verts, (PetscInt64 *)vert_z)); 1134 1135 PetscCall(DMPlexSetChart(dm, 0, local_elems + local_verts)); 1136 for (PetscInt e = 0; e < local_elems; e++) PetscCall(DMPlexSetConeSize(dm, e, PetscPowInt(2, dim))); 1137 PetscCall(DMSetUp(dm)); 1138 { 1139 PetscInt e = 0; 1140 for (ZCode z = layout.zstarts[rank]; z < layout.zstarts[rank + 1]; z++) { 1141 Ijk loc = ZCodeSplit(z); 1142 if (!IjkActive(layout.eextent, loc)) { 1143 z += ZStepOct(z); 1144 continue; 1145 } 1146 PetscInt cone[8], orient[8] = {0}; 1147 for (PetscInt n = 0; n < PetscPowInt(2, dim); n++) { 1148 Ijk inc = closure_dim[dim][n]; 1149 Ijk nloc = {loc.i + inc.i, loc.j + inc.j, loc.k + inc.k}; 1150 ZCode v = ZEncode(nloc); 1151 PetscInt ci = ZCodeFind(v, local_verts, vert_z); 1152 PetscAssert(ci >= 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Did not find neighbor vertex in set"); 1153 cone[n] = local_elems + ci; 1154 } 1155 PetscCall(DMPlexSetCone(dm, e, cone)); 1156 PetscCall(DMPlexSetConeOrientation(dm, e, orient)); 1157 e++; 1158 } 1159 } 1160 1161 PetscCall(DMPlexSymmetrize(dm)); 1162 PetscCall(DMPlexStratify(dm)); 1163 1164 { // Create point SF 1165 PetscSF sf; 1166 PetscCall(PetscSFCreate(PetscObjectComm((PetscObject)dm), &sf)); 1167 PetscInt owned_verts = ZCodeFind(layout.zstarts[rank + 1], local_verts, vert_z); 1168 if (owned_verts < 0) owned_verts = -(owned_verts + 1); // We don't care whether the key was found 1169 PetscInt num_ghosts = local_verts - owned_verts; // Due to sorting, owned vertices always come first 1170 PetscInt *local_ghosts; 1171 PetscSFNode *ghosts; 1172 PetscCall(PetscMalloc1(num_ghosts, &local_ghosts)); 1173 PetscCall(PetscMalloc1(num_ghosts, &ghosts)); 1174 for (PetscInt i = 0; i < num_ghosts;) { 1175 ZCode z = vert_z[owned_verts + i]; 1176 PetscMPIInt remote_rank, remote_count = 0; 1177 1178 PetscCall(PetscMPIIntCast(ZCodeFind(z, size + 1, layout.zstarts), &remote_rank)); 1179 if (remote_rank < 0) remote_rank = -(remote_rank + 1) - 1; 1180 // We have a new remote rank; find all the ghost indices (which are contiguous in vert_z) 1181 1182 // Count the elements on remote_rank 1183 PetscInt remote_elem = ZLayoutElementsOnRank(&layout, remote_rank); 1184 1185 // Traverse vertices and make ghost links 1186 for (ZCode rz = layout.zstarts[remote_rank]; rz < layout.zstarts[remote_rank + 1]; rz++) { 1187 Ijk loc = ZCodeSplit(rz); 1188 if (rz == z) { 1189 local_ghosts[i] = local_elems + owned_verts + i; 1190 ghosts[i].rank = remote_rank; 1191 ghosts[i].index = remote_elem + remote_count; 1192 i++; 1193 if (i == num_ghosts) break; 1194 z = vert_z[owned_verts + i]; 1195 } 1196 if (IjkActive(layout.vextent, loc)) remote_count++; 1197 else rz += ZStepOct(rz); 1198 } 1199 } 1200 PetscCall(PetscSFSetGraph(sf, local_elems + local_verts, num_ghosts, local_ghosts, PETSC_OWN_POINTER, ghosts, PETSC_OWN_POINTER)); 1201 PetscCall(PetscObjectSetName((PetscObject)sf, "SFC Point SF")); 1202 PetscCall(DMSetPointSF(dm, sf)); 1203 PetscCall(PetscSFDestroy(&sf)); 1204 } 1205 { 1206 Vec coordinates; 1207 PetscScalar *coords; 1208 PetscSection coord_section; 1209 PetscInt coord_size; 1210 PetscCall(DMGetCoordinateSection(dm, &coord_section)); 1211 PetscCall(PetscSectionSetNumFields(coord_section, 1)); 1212 PetscCall(PetscSectionSetFieldComponents(coord_section, 0, dim)); 1213 PetscCall(PetscSectionSetChart(coord_section, local_elems, local_elems + local_verts)); 1214 for (PetscInt v = 0; v < local_verts; v++) { 1215 PetscInt point = local_elems + v; 1216 PetscCall(PetscSectionSetDof(coord_section, point, dim)); 1217 PetscCall(PetscSectionSetFieldDof(coord_section, point, 0, dim)); 1218 } 1219 PetscCall(PetscSectionSetUp(coord_section)); 1220 PetscCall(PetscSectionGetStorageSize(coord_section, &coord_size)); 1221 PetscCall(VecCreate(PETSC_COMM_SELF, &coordinates)); 1222 PetscCall(PetscObjectSetName((PetscObject)coordinates, "coordinates")); 1223 PetscCall(VecSetSizes(coordinates, coord_size, PETSC_DETERMINE)); 1224 PetscCall(VecSetBlockSize(coordinates, dim)); 1225 PetscCall(VecSetType(coordinates, VECSTANDARD)); 1226 PetscCall(VecGetArray(coordinates, &coords)); 1227 for (PetscInt v = 0; v < local_verts; v++) { 1228 Ijk loc = ZCodeSplit(vert_z[v]); 1229 coords[v * dim + 0] = lower[0] + loc.i * (upper[0] - lower[0]) / layout.eextent.i; 1230 if (dim > 1) coords[v * dim + 1] = lower[1] + loc.j * (upper[1] - lower[1]) / layout.eextent.j; 1231 if (dim > 2) coords[v * dim + 2] = lower[2] + loc.k * (upper[2] - lower[2]) / layout.eextent.k; 1232 } 1233 PetscCall(VecRestoreArray(coordinates, &coords)); 1234 PetscCall(DMSetCoordinatesLocal(dm, coordinates)); 1235 PetscCall(VecDestroy(&coordinates)); 1236 } 1237 if (interpolate) { 1238 PetscCall(DMPlexInterpolateInPlace_Internal(dm)); 1239 1240 DMLabel label; 1241 PetscCall(DMCreateLabel(dm, "Face Sets")); 1242 PetscCall(DMGetLabel(dm, "Face Sets", &label)); 1243 PetscSegBuffer per_faces[3], donor_face_closure[3], my_donor_faces[3]; 1244 for (PetscInt i = 0; i < 3; i++) { 1245 PetscCall(PetscSegBufferCreate(sizeof(PetscInt), 64, &per_faces[i])); 1246 PetscCall(PetscSegBufferCreate(sizeof(PetscInt), 64, &my_donor_faces[i])); 1247 PetscCall(PetscSegBufferCreate(sizeof(ZCode), 64 * PetscPowInt(2, dim), &donor_face_closure[i])); 1248 } 1249 PetscInt fStart, fEnd, vStart, vEnd; 1250 PetscCall(DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd)); 1251 PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd)); 1252 for (PetscInt f = fStart; f < fEnd; f++) { 1253 PetscInt npoints, *points = NULL, num_fverts = 0, fverts[8]; 1254 PetscCall(DMPlexGetTransitiveClosure(dm, f, PETSC_TRUE, &npoints, &points)); 1255 PetscInt bc_count[6] = {0}; 1256 for (PetscInt i = 0; i < npoints; i++) { 1257 PetscInt p = points[2 * i]; 1258 if (!IsPointInsideStratum(p, vStart, vEnd)) continue; 1259 fverts[num_fverts++] = p; 1260 Ijk loc = ZCodeSplit(vert_z[p - vStart]); 1261 // Convention here matches DMPlexCreateCubeMesh_Internal 1262 bc_count[0] += loc.i == 0; 1263 bc_count[1] += loc.i == layout.vextent.i - 1; 1264 bc_count[2] += loc.j == 0; 1265 bc_count[3] += loc.j == layout.vextent.j - 1; 1266 bc_count[4] += loc.k == 0; 1267 bc_count[5] += loc.k == layout.vextent.k - 1; 1268 } 1269 PetscCall(DMPlexRestoreTransitiveClosure(dm, f, PETSC_TRUE, &npoints, &points)); 1270 for (PetscInt bc = 0, bc_match = 0; bc < 2 * dim; bc++) { 1271 if (bc_count[bc] == PetscPowInt(2, dim - 1)) { 1272 PetscCall(DMPlexOrientPoint(dm, f, face_orient_dim[dim][bc])); 1273 if (periodicity[bc / 2] == DM_BOUNDARY_PERIODIC) { 1274 PetscInt *put; 1275 if (bc % 2 == 0) { // donor face; no label 1276 PetscCall(PetscSegBufferGet(my_donor_faces[bc / 2], 1, &put)); 1277 *put = f; 1278 } else { // periodic face 1279 PetscCall(PetscSegBufferGet(per_faces[bc / 2], 1, &put)); 1280 *put = f; 1281 ZCode *zput; 1282 PetscCall(PetscSegBufferGet(donor_face_closure[bc / 2], num_fverts, &zput)); 1283 for (PetscInt i = 0; i < num_fverts; i++) { 1284 Ijk loc = ZCodeSplit(vert_z[fverts[i] - vStart]); 1285 switch (bc / 2) { 1286 case 0: 1287 loc.i = 0; 1288 break; 1289 case 1: 1290 loc.j = 0; 1291 break; 1292 case 2: 1293 loc.k = 0; 1294 break; 1295 } 1296 *zput++ = ZEncode(loc); 1297 } 1298 } 1299 continue; 1300 } 1301 PetscAssert(bc_match == 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Face matches multiple face sets"); 1302 PetscCall(DMLabelSetValue(label, f, face_marker_dim[dim][bc])); 1303 bc_match++; 1304 } 1305 } 1306 } 1307 // Ensure that the Coordinate DM has our new boundary labels 1308 DM cdm; 1309 PetscCall(DMGetCoordinateDM(dm, &cdm)); 1310 PetscCall(DMCopyLabels(dm, cdm, PETSC_COPY_VALUES, PETSC_FALSE, DM_COPY_LABELS_FAIL)); 1311 if (periodicity[0] == DM_BOUNDARY_PERIODIC || (dim > 1 && periodicity[1] == DM_BOUNDARY_PERIODIC) || (dim > 2 && periodicity[2] == DM_BOUNDARY_PERIODIC)) { 1312 PetscCall(DMPlexCreateBoxMesh_Tensor_SFC_Periodicity_Private(dm, &layout, vert_z, per_faces, lower, upper, periodicity, donor_face_closure, my_donor_faces)); 1313 } 1314 for (PetscInt i = 0; i < 3; i++) { 1315 PetscCall(PetscSegBufferDestroy(&per_faces[i])); 1316 PetscCall(PetscSegBufferDestroy(&donor_face_closure[i])); 1317 PetscCall(PetscSegBufferDestroy(&my_donor_faces[i])); 1318 } 1319 } 1320 PetscCall(PetscFree(layout.zstarts)); 1321 PetscCall(PetscFree(vert_z)); 1322 PetscCall(PetscLogEventEnd(DMPLEX_CreateBoxSFC, dm, 0, 0, 0)); 1323 PetscFunctionReturn(PETSC_SUCCESS); 1324 } 1325 1326 /*@ 1327 DMPlexSetIsoperiodicFaceSF - Express periodicity from an existing mesh 1328 1329 Logically Collective 1330 1331 Input Parameters: 1332 + dm - The `DMPLEX` on which to set periodicity 1333 . num_face_sfs - Number of `PetscSF`s in `face_sfs` 1334 - 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. 1335 1336 Level: advanced 1337 1338 Note: 1339 One can use `-dm_plex_shape zbox` to use this mode of periodicity, wherein the periodic points are distinct both globally 1340 and locally, but are paired when creating a global dof space. 1341 1342 .seealso: [](ch_unstructured), `DMPLEX`, `DMGetGlobalSection()`, `DMPlexGetIsoperiodicFaceSF()` 1343 @*/ 1344 PetscErrorCode DMPlexSetIsoperiodicFaceSF(DM dm, PetscInt num_face_sfs, PetscSF *face_sfs) 1345 { 1346 DM_Plex *plex = (DM_Plex *)dm->data; 1347 1348 PetscFunctionBegin; 1349 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 1350 if (num_face_sfs) PetscCall(PetscObjectComposeFunction((PetscObject)dm, "DMGetIsoperiodicPointSF_C", DMGetIsoperiodicPointSF_Plex)); 1351 else PetscCall(PetscObjectComposeFunction((PetscObject)dm, "DMGetIsoperiodicPointSF_C", NULL)); 1352 if (num_face_sfs == plex->periodic.num_face_sfs && (num_face_sfs == 0 || face_sfs == plex->periodic.face_sfs)) PetscFunctionReturn(PETSC_SUCCESS); 1353 PetscCall(DMSetGlobalSection(dm, NULL)); 1354 1355 for (PetscInt i = 0; i < num_face_sfs; i++) PetscCall(PetscObjectReference((PetscObject)face_sfs[i])); 1356 1357 if (plex->periodic.num_face_sfs > 0) { 1358 for (PetscInt i = 0; i < plex->periodic.num_face_sfs; i++) PetscCall(PetscSFDestroy(&plex->periodic.face_sfs[i])); 1359 PetscCall(PetscFree(plex->periodic.face_sfs)); 1360 } 1361 1362 plex->periodic.num_face_sfs = num_face_sfs; 1363 PetscCall(PetscCalloc1(num_face_sfs, &plex->periodic.face_sfs)); 1364 for (PetscInt i = 0; i < num_face_sfs; i++) plex->periodic.face_sfs[i] = face_sfs[i]; 1365 1366 DM cdm = dm->coordinates[0].dm; // Can't DMGetCoordinateDM because it automatically creates one 1367 if (cdm) { 1368 PetscCall(DMPlexSetIsoperiodicFaceSF(cdm, num_face_sfs, face_sfs)); 1369 if (face_sfs) cdm->periodic.setup = DMPeriodicCoordinateSetUp_Internal; 1370 } 1371 PetscFunctionReturn(PETSC_SUCCESS); 1372 } 1373 1374 /*@C 1375 DMPlexGetIsoperiodicFaceSF - Obtain periodicity for a mesh 1376 1377 Logically Collective 1378 1379 Input Parameter: 1380 . dm - The `DMPLEX` for which to obtain periodic relation 1381 1382 Output Parameters: 1383 + num_face_sfs - Number of `PetscSF`s in the array 1384 - 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. 1385 1386 Level: advanced 1387 1388 .seealso: [](ch_unstructured), `DMPLEX`, `DMGetGlobalSection()`, `DMPlexSetIsoperiodicFaceSF()` 1389 @*/ 1390 PetscErrorCode DMPlexGetIsoperiodicFaceSF(DM dm, PetscInt *num_face_sfs, const PetscSF **face_sfs) 1391 { 1392 DM_Plex *plex = (DM_Plex *)dm->data; 1393 1394 PetscFunctionBegin; 1395 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 1396 *face_sfs = plex->periodic.face_sfs; 1397 *num_face_sfs = plex->periodic.num_face_sfs; 1398 PetscFunctionReturn(PETSC_SUCCESS); 1399 } 1400 1401 /*@C 1402 DMPlexSetIsoperiodicFaceTransform - set geometric transform from donor to periodic points 1403 1404 Logically Collective 1405 1406 Input Parameters: 1407 + dm - `DMPLEX` that has been configured with `DMPlexSetIsoperiodicFaceSF()` 1408 . n - Number of transforms in array 1409 - t - Array of 4x4 affine transformation basis. 1410 1411 Level: advanced 1412 1413 Notes: 1414 Affine transforms are 4x4 matrices in which the leading 3x3 block expresses a rotation (or identity for no rotation), 1415 the last column contains a translation vector, and the bottom row is all zero except the last entry, which must always 1416 be 1. This representation is common in geometric modeling and allows affine transformations to be composed using 1417 simple matrix multiplication. 1418 1419 Although the interface accepts a general affine transform, only affine translation is supported at present. 1420 1421 Developer Notes: 1422 This interface should be replaced by making BasisTransform public, expanding it to support affine representations, and 1423 adding GPU implementations to apply the G2L/L2G transforms. 1424 1425 .seealso: [](ch_unstructured), `DMPLEX`, `DMGetGlobalSection()`, `DMPlexSetIsoperiodicFaceSF()` 1426 @*/ 1427 PetscErrorCode DMPlexSetIsoperiodicFaceTransform(DM dm, PetscInt n, const PetscScalar t[]) 1428 { 1429 DM_Plex *plex = (DM_Plex *)dm->data; 1430 1431 PetscFunctionBegin; 1432 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 1433 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); 1434 1435 PetscCall(PetscFree(plex->periodic.transform)); 1436 PetscCall(PetscMalloc1(n, &plex->periodic.transform)); 1437 for (PetscInt i = 0; i < n; i++) { 1438 for (PetscInt j = 0; j < 4; j++) { 1439 for (PetscInt k = 0; k < 4; k++) { 1440 PetscCheck(j != k || t[i * 16 + j * 4 + k] == 1., PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Rotated transforms not supported"); 1441 plex->periodic.transform[i][j][k] = t[i * 16 + j * 4 + k]; 1442 } 1443 } 1444 } 1445 PetscFunctionReturn(PETSC_SUCCESS); 1446 } 1447