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