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