Actual source code: ex6.c
1: static char help[] = "Spectral element access patterns with Plex\n\n";
3: #include <petscdmplex.h>
5: typedef struct {
6: PetscInt Nf; /* Number of fields */
7: PetscInt *Nc; /* Number of components per field */
8: PetscInt *k; /* Spectral order per field */
9: } AppCtx;
11: static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
12: {
13: PetscInt len;
14: PetscBool flg;
16: PetscFunctionBeginUser;
17: options->Nf = 0;
18: options->Nc = NULL;
19: options->k = NULL;
21: PetscOptionsBegin(comm, "", "SEM Problem Options", "DMPLEX");
22: PetscCall(PetscOptionsBoundedInt("-num_fields", "The number of fields", "ex6.c", options->Nf, &options->Nf, NULL, 0));
23: if (options->Nf) {
24: len = options->Nf;
25: PetscCall(PetscMalloc1(len, &options->Nc));
26: PetscCall(PetscOptionsIntArray("-num_components", "The number of components per field", "ex6.c", options->Nc, &len, &flg));
27: PetscCheck(!flg || !(len != options->Nf), PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Length of components array is %" PetscInt_FMT " should be %" PetscInt_FMT, len, options->Nf);
28: len = options->Nf;
29: PetscCall(PetscMalloc1(len, &options->k));
30: PetscCall(PetscOptionsIntArray("-order", "The spectral order per field", "ex6.c", options->k, &len, &flg));
31: PetscCheck(!flg || !(len != options->Nf), PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Length of order array is %" PetscInt_FMT " should be %" PetscInt_FMT, len, options->Nf);
32: }
33: PetscOptionsEnd();
34: PetscFunctionReturn(PETSC_SUCCESS);
35: }
37: static PetscErrorCode LoadData2D(DM dm, PetscInt Ni, PetscInt Nj, PetscInt clSize, Vec u, AppCtx *user)
38: {
39: PetscInt i, j, f, c;
40: PetscScalar *closure;
42: PetscFunctionBeginUser;
43: PetscCall(PetscMalloc1(clSize, &closure));
44: for (j = 0; j < Nj; ++j) {
45: for (i = 0; i < Ni; ++i) {
46: PetscInt ki, kj, o = 0;
47: PetscCall(PetscArrayzero(closure, clSize));
49: for (f = 0; f < user->Nf; ++f) {
50: PetscInt ioff = i * user->k[f], joff = j * user->k[f];
52: for (kj = 0; kj <= user->k[f]; ++kj) {
53: for (ki = 0; ki <= user->k[f]; ++ki) {
54: for (c = 0; c < user->Nc[f]; ++c) closure[o++] = ((kj + joff) * (Ni * user->k[f] + 1) + ki + ioff) * user->Nc[f] + c;
55: }
56: }
57: }
58: PetscCall(DMPlexVecSetClosure(dm, NULL, u, j * Ni + i, closure, INSERT_VALUES));
59: }
60: }
61: PetscCall(PetscFree(closure));
62: PetscFunctionReturn(PETSC_SUCCESS);
63: }
65: static PetscErrorCode LoadData3D(DM dm, PetscInt Ni, PetscInt Nj, PetscInt Nk, PetscInt clSize, Vec u, AppCtx *user)
66: {
67: PetscInt i, j, k, f, c;
68: PetscScalar *closure;
70: PetscFunctionBeginUser;
71: PetscCall(PetscMalloc1(clSize, &closure));
72: for (k = 0; k < Nk; ++k) {
73: for (j = 0; j < Nj; ++j) {
74: for (i = 0; i < Ni; ++i) {
75: PetscInt ki, kj, kk, o = 0;
76: PetscCall(PetscArrayzero(closure, clSize));
78: for (f = 0; f < user->Nf; ++f) {
79: PetscInt ioff = i * user->k[f], joff = j * user->k[f], koff = k * user->k[f];
81: for (kk = 0; kk <= user->k[f]; ++kk) {
82: for (kj = 0; kj <= user->k[f]; ++kj) {
83: for (ki = 0; ki <= user->k[f]; ++ki) {
84: for (c = 0; c < user->Nc[f]; ++c) closure[o++] = (((kk + koff) * (Nj * user->k[f] + 1) + kj + joff) * (Ni * user->k[f] + 1) + ki + ioff) * user->Nc[f] + c;
85: }
86: }
87: }
88: }
89: PetscCall(DMPlexVecSetClosure(dm, NULL, u, (k * Nj + j) * Ni + i, closure, INSERT_VALUES));
90: }
91: }
92: }
93: PetscCall(PetscFree(closure));
94: PetscFunctionReturn(PETSC_SUCCESS);
95: }
97: static PetscErrorCode CheckPoint(DM dm, Vec u, PetscInt point, AppCtx *user)
98: {
99: PetscSection s;
100: PetscScalar *a;
101: const PetscScalar *array;
102: PetscInt dof, d;
104: PetscFunctionBeginUser;
105: PetscCall(DMGetLocalSection(dm, &s));
106: PetscCall(VecGetArrayRead(u, &array));
107: PetscCall(DMPlexPointLocalRead(dm, point, array, &a));
108: PetscCall(PetscSectionGetDof(s, point, &dof));
109: PetscCall(PetscPrintf(PETSC_COMM_SELF, "Point %" PetscInt_FMT ": ", point));
110: for (d = 0; d < dof; ++d) {
111: if (d > 0) PetscCall(PetscPrintf(PETSC_COMM_SELF, ", "));
112: PetscCall(PetscPrintf(PETSC_COMM_SELF, "%2.0f", (double)PetscRealPart(a[d])));
113: }
114: PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n"));
115: PetscCall(VecRestoreArrayRead(u, &array));
116: PetscFunctionReturn(PETSC_SUCCESS);
117: }
119: static PetscErrorCode ReadData2D(DM dm, Vec u, AppCtx *user)
120: {
121: PetscInt cStart, cEnd, cell;
123: PetscFunctionBeginUser;
124: PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
125: for (cell = cStart; cell < cEnd; ++cell) {
126: PetscScalar *closure = NULL;
127: PetscInt closureSize, ki, kj, f, c, foff = 0;
129: PetscCall(DMPlexVecGetClosure(dm, NULL, u, cell, &closureSize, &closure));
130: PetscCall(PetscPrintf(PETSC_COMM_SELF, "Cell %" PetscInt_FMT "\n", cell));
131: for (f = 0; f < user->Nf; ++f) {
132: PetscCall(PetscPrintf(PETSC_COMM_SELF, " Field %" PetscInt_FMT "\n", f));
133: for (kj = user->k[f]; kj >= 0; --kj) {
134: for (ki = 0; ki <= user->k[f]; ++ki) {
135: if (ki > 0) PetscCall(PetscPrintf(PETSC_COMM_SELF, " "));
136: for (c = 0; c < user->Nc[f]; ++c) {
137: if (c > 0) PetscCall(PetscPrintf(PETSC_COMM_SELF, ","));
138: PetscCall(PetscPrintf(PETSC_COMM_SELF, "%2.0f", (double)PetscRealPart(closure[(kj * (user->k[f] + 1) + ki) * user->Nc[f] + c + foff])));
139: }
140: }
141: PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n"));
142: }
143: PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n\n"));
144: foff += PetscSqr(user->k[f] + 1);
145: }
146: PetscCall(DMPlexVecRestoreClosure(dm, NULL, u, cell, &closureSize, &closure));
147: PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n\n"));
148: }
149: PetscFunctionReturn(PETSC_SUCCESS);
150: }
152: static PetscErrorCode ReadData3D(DM dm, Vec u, AppCtx *user)
153: {
154: PetscInt cStart, cEnd, cell;
156: PetscFunctionBeginUser;
157: PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
158: for (cell = cStart; cell < cEnd; ++cell) {
159: PetscScalar *closure = NULL;
160: PetscInt closureSize, ki, kj, kk, f, c, foff = 0;
162: PetscCall(DMPlexVecGetClosure(dm, NULL, u, cell, &closureSize, &closure));
163: PetscCall(PetscPrintf(PETSC_COMM_SELF, "Cell %" PetscInt_FMT "\n", cell));
164: for (f = 0; f < user->Nf; ++f) {
165: PetscCall(PetscPrintf(PETSC_COMM_SELF, " Field %" PetscInt_FMT "\n", f));
166: for (kk = user->k[f]; kk >= 0; --kk) {
167: for (kj = user->k[f]; kj >= 0; --kj) {
168: for (ki = 0; ki <= user->k[f]; ++ki) {
169: if (ki > 0) PetscCall(PetscPrintf(PETSC_COMM_SELF, " "));
170: for (c = 0; c < user->Nc[f]; ++c) {
171: if (c > 0) PetscCall(PetscPrintf(PETSC_COMM_SELF, ","));
172: PetscCall(PetscPrintf(PETSC_COMM_SELF, "%2.0f", (double)PetscRealPart(closure[((kk * (user->k[f] + 1) + kj) * (user->k[f] + 1) + ki) * user->Nc[f] + c + foff])));
173: }
174: }
175: PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n"));
176: }
177: PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n"));
178: }
179: PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n\n"));
180: foff += PetscSqr(user->k[f] + 1);
181: }
182: PetscCall(DMPlexVecRestoreClosure(dm, NULL, u, cell, &closureSize, &closure));
183: PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n\n"));
184: }
185: PetscFunctionReturn(PETSC_SUCCESS);
186: }
188: static PetscErrorCode SetSymmetries(DM dm, PetscSection s, AppCtx *user)
189: {
190: PetscInt dim, f, o, i, j, k, c, d;
191: DMLabel depthLabel;
193: PetscFunctionBegin;
194: PetscCall(DMGetDimension(dm, &dim));
195: PetscCall(DMGetLabel(dm, "depth", &depthLabel));
196: for (f = 0; f < user->Nf; f++) {
197: PetscSectionSym sym;
199: if (user->k[f] < 3) continue; /* No symmetries needed for order < 3, because no cell, facet, edge or vertex has more than one node */
200: PetscCall(PetscSectionSymCreateLabel(PetscObjectComm((PetscObject)s), depthLabel, &sym));
202: for (d = 0; d <= dim; d++) {
203: if (d == 1) {
204: PetscInt numDof = user->k[f] - 1;
205: PetscInt numComp = user->Nc[f];
206: PetscInt minOrnt = -1;
207: PetscInt maxOrnt = 1;
208: PetscInt **perms;
210: PetscCall(PetscCalloc1(maxOrnt - minOrnt, &perms));
211: for (o = minOrnt; o < maxOrnt; o++) {
212: PetscInt *perm;
214: if (!o) { /* identity */
215: perms[o - minOrnt] = NULL;
216: } else {
217: PetscCall(PetscMalloc1(numDof * numComp, &perm));
218: for (i = numDof - 1, k = 0; i >= 0; i--) {
219: for (j = 0; j < numComp; j++, k++) perm[k] = i * numComp + j;
220: }
221: perms[o - minOrnt] = perm;
222: }
223: }
224: PetscCall(PetscSectionSymLabelSetStratum(sym, d, numDof * numComp, minOrnt, maxOrnt, PETSC_OWN_POINTER, (const PetscInt **)perms, NULL));
225: } else if (d == 2) {
226: PetscInt perEdge = user->k[f] - 1;
227: PetscInt numDof = perEdge * perEdge;
228: PetscInt numComp = user->Nc[f];
229: PetscInt minOrnt = -4;
230: PetscInt maxOrnt = 4;
231: PetscInt **perms;
233: PetscCall(PetscCalloc1(maxOrnt - minOrnt, &perms));
234: for (o = minOrnt; o < maxOrnt; o++) {
235: PetscInt *perm;
237: if (!o) continue; /* identity */
238: PetscCall(PetscMalloc1(numDof * numComp, &perm));
239: /* We want to perm[k] to list which *localArray* position the *sectionArray* position k should go to for the given orientation*/
240: switch (o) {
241: case 0:
242: break; /* identity */
243: case -2: /* flip along (-1,-1)--( 1, 1), which swaps edges 0 and 3 and edges 1 and 2. This swaps the i and j variables */
244: for (i = 0, k = 0; i < perEdge; i++) {
245: for (j = 0; j < perEdge; j++, k++) {
246: for (c = 0; c < numComp; c++) perm[k * numComp + c] = (perEdge * j + i) * numComp + c;
247: }
248: }
249: break;
250: case -1: /* flip along (-1, 0)--( 1, 0), which swaps edges 0 and 2. This reverses the i variable */
251: for (i = 0, k = 0; i < perEdge; i++) {
252: for (j = 0; j < perEdge; j++, k++) {
253: for (c = 0; c < numComp; c++) perm[k * numComp + c] = (perEdge * (perEdge - 1 - i) + j) * numComp + c;
254: }
255: }
256: break;
257: case -4: /* flip along ( 1,-1)--(-1, 1), which swaps edges 0 and 1 and edges 2 and 3. This swaps the i and j variables and reverse both */
258: for (i = 0, k = 0; i < perEdge; i++) {
259: for (j = 0; j < perEdge; j++, k++) {
260: for (c = 0; c < numComp; c++) perm[k * numComp + c] = (perEdge * (perEdge - 1 - j) + (perEdge - 1 - i)) * numComp + c;
261: }
262: }
263: break;
264: case -3: /* flip along ( 0,-1)--( 0, 1), which swaps edges 3 and 1. This reverses the j variable */
265: for (i = 0, k = 0; i < perEdge; i++) {
266: for (j = 0; j < perEdge; j++, k++) {
267: for (c = 0; c < numComp; c++) perm[k * numComp + c] = (perEdge * i + (perEdge - 1 - j)) * numComp + c;
268: }
269: }
270: break;
271: case 1: /* rotate section edge 1 to local edge 0. This swaps the i and j variables and then reverses the j variable */
272: for (i = 0, k = 0; i < perEdge; i++) {
273: for (j = 0; j < perEdge; j++, k++) {
274: for (c = 0; c < numComp; c++) perm[k * numComp + c] = (perEdge * (perEdge - 1 - j) + i) * numComp + c;
275: }
276: }
277: break;
278: case 2: /* rotate section edge 2 to local edge 0. This reverse both i and j variables */
279: for (i = 0, k = 0; i < perEdge; i++) {
280: for (j = 0; j < perEdge; j++, k++) {
281: for (c = 0; c < numComp; c++) perm[k * numComp + c] = (perEdge * (perEdge - 1 - i) + (perEdge - 1 - j)) * numComp + c;
282: }
283: }
284: break;
285: case 3: /* rotate section edge 3 to local edge 0. This swaps the i and j variables and then reverses the i variable */
286: for (i = 0, k = 0; i < perEdge; i++) {
287: for (j = 0; j < perEdge; j++, k++) {
288: for (c = 0; c < numComp; c++) perm[k * numComp + c] = (perEdge * j + (perEdge - 1 - i)) * numComp + c;
289: }
290: }
291: break;
292: default:
293: break;
294: }
295: perms[o - minOrnt] = perm;
296: }
297: PetscCall(PetscSectionSymLabelSetStratum(sym, d, numDof * numComp, minOrnt, maxOrnt, PETSC_OWN_POINTER, (const PetscInt **)perms, NULL));
298: }
299: }
300: PetscCall(PetscSectionSetFieldSym(s, f, sym));
301: PetscCall(PetscSectionSymDestroy(&sym));
302: }
303: PetscCall(PetscSectionViewFromOptions(s, NULL, "-section_with_sym_view"));
304: PetscFunctionReturn(PETSC_SUCCESS);
305: }
307: int main(int argc, char **argv)
308: {
309: DM dm;
310: PetscSection s;
311: Vec u;
312: AppCtx user;
313: PetscInt dim, size = 0, f;
315: PetscFunctionBeginUser;
316: PetscCall(PetscInitialize(&argc, &argv, NULL, help));
317: PetscCall(ProcessOptions(PETSC_COMM_WORLD, &user));
318: PetscCall(DMCreate(PETSC_COMM_WORLD, &dm));
319: PetscCall(DMSetType(dm, DMPLEX));
320: PetscCall(DMSetFromOptions(dm));
321: PetscCall(DMViewFromOptions(dm, NULL, "-dm_view"));
322: PetscCall(DMGetDimension(dm, &dim));
323: /* Create a section for SEM order k */
324: {
325: PetscInt *numDof, d;
327: PetscCall(PetscMalloc1(user.Nf * (dim + 1), &numDof));
328: for (f = 0; f < user.Nf; ++f) {
329: for (d = 0; d <= dim; ++d) numDof[f * (dim + 1) + d] = PetscPowInt(user.k[f] - 1, d) * user.Nc[f];
330: size += PetscPowInt(user.k[f] + 1, d) * user.Nc[f];
331: }
332: PetscCall(DMSetNumFields(dm, user.Nf));
333: PetscCall(DMPlexCreateSection(dm, NULL, user.Nc, numDof, 0, NULL, NULL, NULL, NULL, &s));
334: PetscCall(SetSymmetries(dm, s, &user));
335: PetscCall(PetscFree(numDof));
336: }
337: PetscCall(DMSetLocalSection(dm, s));
338: /* Create spectral ordering and load in data */
339: PetscCall(DMPlexSetClosurePermutationTensor(dm, PETSC_DETERMINE, NULL));
340: PetscCall(DMGetLocalVector(dm, &u));
341: switch (dim) {
342: case 2:
343: PetscCall(LoadData2D(dm, 2, 2, size, u, &user));
344: break;
345: case 3:
346: PetscCall(LoadData3D(dm, 2, 2, 2, size, u, &user));
347: break;
348: }
349: /* Remove ordering and check some values */
350: PetscCall(PetscSectionSetClosurePermutation(s, (PetscObject)dm, dim, NULL));
351: switch (dim) {
352: case 2:
353: PetscCall(CheckPoint(dm, u, 0, &user));
354: PetscCall(CheckPoint(dm, u, 13, &user));
355: PetscCall(CheckPoint(dm, u, 15, &user));
356: PetscCall(CheckPoint(dm, u, 19, &user));
357: break;
358: case 3:
359: PetscCall(CheckPoint(dm, u, 0, &user));
360: PetscCall(CheckPoint(dm, u, 13, &user));
361: PetscCall(CheckPoint(dm, u, 15, &user));
362: PetscCall(CheckPoint(dm, u, 19, &user));
363: break;
364: }
365: /* Recreate spectral ordering and read out data */
366: PetscCall(DMPlexSetClosurePermutationTensor(dm, PETSC_DETERMINE, s));
367: switch (dim) {
368: case 2:
369: PetscCall(ReadData2D(dm, u, &user));
370: break;
371: case 3:
372: PetscCall(ReadData3D(dm, u, &user));
373: break;
374: }
375: PetscCall(DMRestoreLocalVector(dm, &u));
376: PetscCall(PetscSectionDestroy(&s));
377: PetscCall(DMDestroy(&dm));
378: PetscCall(PetscFree(user.Nc));
379: PetscCall(PetscFree(user.k));
380: PetscCall(PetscFinalize());
381: return 0;
382: }
384: /*TEST
386: # Spectral ordering 2D 0-5
387: testset:
388: args: -dm_plex_simplex 0 -dm_plex_box_faces 2,2
390: test:
391: suffix: 0
392: args: -num_fields 1 -num_components 1 -order 2
393: test:
394: suffix: 1
395: args: -num_fields 1 -num_components 1 -order 3
396: test:
397: suffix: 2
398: args: -num_fields 1 -num_components 1 -order 5
399: test:
400: suffix: 3
401: args: -num_fields 1 -num_components 2 -order 2
402: test:
403: suffix: 4
404: args: -num_fields 2 -num_components 1,1 -order 2,2
405: test:
406: suffix: 5
407: args: -num_fields 2 -num_components 1,2 -order 2,3
409: # Spectral ordering 3D 6-11
410: testset:
411: args: -dm_plex_dim 3 -dm_plex_simplex 0 -dm_plex_box_faces 2,2,2
413: test:
414: suffix: 6
415: args: -num_fields 1 -num_components 1 -order 2
416: test:
417: suffix: 7
418: args: -num_fields 1 -num_components 1 -order 3
419: test:
420: suffix: 8
421: args: -num_fields 1 -num_components 1 -order 5
422: test:
423: suffix: 9
424: args: -num_fields 1 -num_components 2 -order 2
425: test:
426: suffix: 10
427: args: -num_fields 2 -num_components 1,1 -order 2,2
428: test:
429: suffix: 11
430: args: -num_fields 2 -num_components 1,2 -order 2,3
432: TEST*/