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;
17: options->Nf = 0;
18: options->Nc = NULL;
19: options->k = NULL;
21: PetscOptionsBegin(comm, "", "SEM Problem Options", "DMPLEX");
22: PetscOptionsBoundedInt("-num_fields", "The number of fields", "ex6.c", options->Nf, &options->Nf, NULL, 0);
23: if (options->Nf) {
24: len = options->Nf;
25: PetscMalloc1(len, &options->Nc);
26: PetscOptionsIntArray("-num_components", "The number of components per field", "ex6.c", options->Nc, &len, &flg);
28: len = options->Nf;
29: PetscMalloc1(len, &options->k);
30: PetscOptionsIntArray("-order", "The spectral order per field", "ex6.c", options->k, &len, &flg);
32: }
33: PetscOptionsEnd();
34: return 0;
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;
43: PetscMalloc1(clSize, &closure);
44: for (j = 0; j < Nj; ++j) {
45: for (i = 0; i < Ni; ++i) {
46: PetscInt ki, kj, o = 0;
47: 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: DMPlexVecSetClosure(dm, NULL, u, j * Ni + i, closure, INSERT_VALUES);
59: }
60: }
61: PetscFree(closure);
62: return 0;
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;
71: 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: 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: DMPlexVecSetClosure(dm, NULL, u, (k * Nj + j) * Ni + i, closure, INSERT_VALUES);
90: }
91: }
92: }
93: PetscFree(closure);
94: return 0;
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;
105: DMGetLocalSection(dm, &s);
106: VecGetArrayRead(u, &array);
107: DMPlexPointLocalRead(dm, point, array, &a);
108: PetscSectionGetDof(s, point, &dof);
109: PetscPrintf(PETSC_COMM_SELF, "Point %" PetscInt_FMT ": ", point);
110: for (d = 0; d < dof; ++d) {
111: if (d > 0) PetscPrintf(PETSC_COMM_SELF, ", ");
112: PetscPrintf(PETSC_COMM_SELF, "%2.0f", (double)PetscRealPart(a[d]));
113: }
114: PetscPrintf(PETSC_COMM_SELF, "\n");
115: VecRestoreArrayRead(u, &array);
116: return 0;
117: }
119: static PetscErrorCode ReadData2D(DM dm, Vec u, AppCtx *user)
120: {
121: PetscInt cStart, cEnd, cell;
124: 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: DMPlexVecGetClosure(dm, NULL, u, cell, &closureSize, &closure);
130: PetscPrintf(PETSC_COMM_SELF, "Cell %" PetscInt_FMT "\n", cell);
131: for (f = 0; f < user->Nf; ++f) {
132: 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) PetscPrintf(PETSC_COMM_SELF, " ");
136: for (c = 0; c < user->Nc[f]; ++c) {
137: if (c > 0) PetscPrintf(PETSC_COMM_SELF, ",");
138: PetscPrintf(PETSC_COMM_SELF, "%2.0f", (double)PetscRealPart(closure[(kj * (user->k[f] + 1) + ki) * user->Nc[f] + c + foff]));
139: }
140: }
141: PetscPrintf(PETSC_COMM_SELF, "\n");
142: }
143: PetscPrintf(PETSC_COMM_SELF, "\n\n");
144: foff += PetscSqr(user->k[f] + 1);
145: }
146: DMPlexVecRestoreClosure(dm, NULL, u, cell, &closureSize, &closure);
147: PetscPrintf(PETSC_COMM_SELF, "\n\n");
148: }
149: return 0;
150: }
152: static PetscErrorCode ReadData3D(DM dm, Vec u, AppCtx *user)
153: {
154: PetscInt cStart, cEnd, cell;
157: 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: DMPlexVecGetClosure(dm, NULL, u, cell, &closureSize, &closure);
163: PetscPrintf(PETSC_COMM_SELF, "Cell %" PetscInt_FMT "\n", cell);
164: for (f = 0; f < user->Nf; ++f) {
165: 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) PetscPrintf(PETSC_COMM_SELF, " ");
170: for (c = 0; c < user->Nc[f]; ++c) {
171: if (c > 0) PetscPrintf(PETSC_COMM_SELF, ",");
172: 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: PetscPrintf(PETSC_COMM_SELF, "\n");
176: }
177: PetscPrintf(PETSC_COMM_SELF, "\n");
178: }
179: PetscPrintf(PETSC_COMM_SELF, "\n\n");
180: foff += PetscSqr(user->k[f] + 1);
181: }
182: DMPlexVecRestoreClosure(dm, NULL, u, cell, &closureSize, &closure);
183: PetscPrintf(PETSC_COMM_SELF, "\n\n");
184: }
185: return 0;
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: DMGetDimension(dm, &dim);
194: DMGetLabel(dm, "depth", &depthLabel);
195: for (f = 0; f < user->Nf; f++) {
196: PetscSectionSym sym;
198: if (user->k[f] < 3) continue; /* No symmetries needed for order < 3, because no cell, facet, edge or vertex has more than one node */
199: PetscSectionSymCreateLabel(PetscObjectComm((PetscObject)s), depthLabel, &sym);
201: for (d = 0; d <= dim; d++) {
202: if (d == 1) {
203: PetscInt numDof = user->k[f] - 1;
204: PetscInt numComp = user->Nc[f];
205: PetscInt minOrnt = -1;
206: PetscInt maxOrnt = 1;
207: PetscInt **perms;
209: PetscCalloc1(maxOrnt - minOrnt, &perms);
210: for (o = minOrnt; o < maxOrnt; o++) {
211: PetscInt *perm;
213: if (!o) { /* identity */
214: perms[o - minOrnt] = NULL;
215: } else {
216: PetscMalloc1(numDof * numComp, &perm);
217: for (i = numDof - 1, k = 0; i >= 0; i--) {
218: for (j = 0; j < numComp; j++, k++) perm[k] = i * numComp + j;
219: }
220: perms[o - minOrnt] = perm;
221: }
222: }
223: PetscSectionSymLabelSetStratum(sym, d, numDof * numComp, minOrnt, maxOrnt, PETSC_OWN_POINTER, (const PetscInt **)perms, NULL);
224: } else if (d == 2) {
225: PetscInt perEdge = user->k[f] - 1;
226: PetscInt numDof = perEdge * perEdge;
227: PetscInt numComp = user->Nc[f];
228: PetscInt minOrnt = -4;
229: PetscInt maxOrnt = 4;
230: PetscInt **perms;
232: PetscCalloc1(maxOrnt - minOrnt, &perms);
233: for (o = minOrnt; o < maxOrnt; o++) {
234: PetscInt *perm;
236: if (!o) continue; /* identity */
237: PetscMalloc1(numDof * numComp, &perm);
238: /* We want to perm[k] to list which *localArray* position the *sectionArray* position k should go to for the given orientation*/
239: switch (o) {
240: case 0:
241: break; /* identity */
242: 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 */
243: for (i = 0, k = 0; i < perEdge; i++) {
244: for (j = 0; j < perEdge; j++, k++) {
245: for (c = 0; c < numComp; c++) perm[k * numComp + c] = (perEdge * j + i) * numComp + c;
246: }
247: }
248: break;
249: case -1: /* flip along (-1, 0)--( 1, 0), which swaps edges 0 and 2. This reverses the i variable */
250: for (i = 0, k = 0; i < perEdge; i++) {
251: for (j = 0; j < perEdge; j++, k++) {
252: for (c = 0; c < numComp; c++) perm[k * numComp + c] = (perEdge * (perEdge - 1 - i) + j) * numComp + c;
253: }
254: }
255: break;
256: 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 */
257: for (i = 0, k = 0; i < perEdge; i++) {
258: for (j = 0; j < perEdge; j++, k++) {
259: for (c = 0; c < numComp; c++) perm[k * numComp + c] = (perEdge * (perEdge - 1 - j) + (perEdge - 1 - i)) * numComp + c;
260: }
261: }
262: break;
263: case -3: /* flip along ( 0,-1)--( 0, 1), which swaps edges 3 and 1. This reverses the j variable */
264: for (i = 0, k = 0; i < perEdge; i++) {
265: for (j = 0; j < perEdge; j++, k++) {
266: for (c = 0; c < numComp; c++) perm[k * numComp + c] = (perEdge * i + (perEdge - 1 - j)) * numComp + c;
267: }
268: }
269: break;
270: case 1: /* rotate section edge 1 to local edge 0. This swaps the i and j variables and then reverses the j variable */
271: for (i = 0, k = 0; i < perEdge; i++) {
272: for (j = 0; j < perEdge; j++, k++) {
273: for (c = 0; c < numComp; c++) perm[k * numComp + c] = (perEdge * (perEdge - 1 - j) + i) * numComp + c;
274: }
275: }
276: break;
277: case 2: /* rotate section edge 2 to local edge 0. This reverse both i and j variables */
278: for (i = 0, k = 0; i < perEdge; i++) {
279: for (j = 0; j < perEdge; j++, k++) {
280: for (c = 0; c < numComp; c++) perm[k * numComp + c] = (perEdge * (perEdge - 1 - i) + (perEdge - 1 - j)) * numComp + c;
281: }
282: }
283: break;
284: case 3: /* rotate section edge 3 to local edge 0. This swaps the i and j variables and then reverses the i variable */
285: for (i = 0, k = 0; i < perEdge; i++) {
286: for (j = 0; j < perEdge; j++, k++) {
287: for (c = 0; c < numComp; c++) perm[k * numComp + c] = (perEdge * j + (perEdge - 1 - i)) * numComp + c;
288: }
289: }
290: break;
291: default:
292: break;
293: }
294: perms[o - minOrnt] = perm;
295: }
296: PetscSectionSymLabelSetStratum(sym, d, numDof * numComp, minOrnt, maxOrnt, PETSC_OWN_POINTER, (const PetscInt **)perms, NULL);
297: }
298: }
299: PetscSectionSetFieldSym(s, f, sym);
300: PetscSectionSymDestroy(&sym);
301: }
302: PetscSectionViewFromOptions(s, NULL, "-section_with_sym_view");
303: return 0;
304: }
306: int main(int argc, char **argv)
307: {
308: DM dm;
309: PetscSection s;
310: Vec u;
311: AppCtx user;
312: PetscInt dim, size = 0, f;
315: PetscInitialize(&argc, &argv, NULL, help);
316: ProcessOptions(PETSC_COMM_WORLD, &user);
317: DMCreate(PETSC_COMM_WORLD, &dm);
318: DMSetType(dm, DMPLEX);
319: DMSetFromOptions(dm);
320: DMViewFromOptions(dm, NULL, "-dm_view");
321: DMGetDimension(dm, &dim);
322: /* Create a section for SEM order k */
323: {
324: PetscInt *numDof, d;
326: PetscMalloc1(user.Nf * (dim + 1), &numDof);
327: for (f = 0; f < user.Nf; ++f) {
328: for (d = 0; d <= dim; ++d) numDof[f * (dim + 1) + d] = PetscPowInt(user.k[f] - 1, d) * user.Nc[f];
329: size += PetscPowInt(user.k[f] + 1, d) * user.Nc[f];
330: }
331: DMSetNumFields(dm, user.Nf);
332: DMPlexCreateSection(dm, NULL, user.Nc, numDof, 0, NULL, NULL, NULL, NULL, &s);
333: SetSymmetries(dm, s, &user);
334: PetscFree(numDof);
335: }
336: DMSetLocalSection(dm, s);
337: /* Create spectral ordering and load in data */
338: DMPlexSetClosurePermutationTensor(dm, PETSC_DETERMINE, NULL);
339: DMGetLocalVector(dm, &u);
340: switch (dim) {
341: case 2:
342: LoadData2D(dm, 2, 2, size, u, &user);
343: break;
344: case 3:
345: LoadData3D(dm, 2, 2, 2, size, u, &user);
346: break;
347: }
348: /* Remove ordering and check some values */
349: PetscSectionSetClosurePermutation(s, (PetscObject)dm, dim, NULL);
350: switch (dim) {
351: case 2:
352: CheckPoint(dm, u, 0, &user);
353: CheckPoint(dm, u, 13, &user);
354: CheckPoint(dm, u, 15, &user);
355: CheckPoint(dm, u, 19, &user);
356: break;
357: case 3:
358: CheckPoint(dm, u, 0, &user);
359: CheckPoint(dm, u, 13, &user);
360: CheckPoint(dm, u, 15, &user);
361: CheckPoint(dm, u, 19, &user);
362: break;
363: }
364: /* Recreate spectral ordering and read out data */
365: DMPlexSetClosurePermutationTensor(dm, PETSC_DETERMINE, s);
366: switch (dim) {
367: case 2:
368: ReadData2D(dm, u, &user);
369: break;
370: case 3:
371: ReadData3D(dm, u, &user);
372: break;
373: }
374: DMRestoreLocalVector(dm, &u);
375: PetscSectionDestroy(&s);
376: DMDestroy(&dm);
377: PetscFree(user.Nc);
378: PetscFree(user.k);
379: PetscFinalize();
380: return 0;
381: }
383: /*TEST
385: # Spectral ordering 2D 0-5
386: testset:
387: args: -dm_plex_simplex 0 -dm_plex_box_faces 2,2
389: test:
390: suffix: 0
391: args: -num_fields 1 -num_components 1 -order 2
392: test:
393: suffix: 1
394: args: -num_fields 1 -num_components 1 -order 3
395: test:
396: suffix: 2
397: args: -num_fields 1 -num_components 1 -order 5
398: test:
399: suffix: 3
400: args: -num_fields 1 -num_components 2 -order 2
401: test:
402: suffix: 4
403: args: -num_fields 2 -num_components 1,1 -order 2,2
404: test:
405: suffix: 5
406: args: -num_fields 2 -num_components 1,2 -order 2,3
408: # Spectral ordering 3D 6-11
409: testset:
410: args: -dm_plex_dim 3 -dm_plex_simplex 0 -dm_plex_box_faces 2,2,2
412: test:
413: suffix: 6
414: args: -num_fields 1 -num_components 1 -order 2
415: test:
416: suffix: 7
417: args: -num_fields 1 -num_components 1 -order 3
418: test:
419: suffix: 8
420: args: -num_fields 1 -num_components 1 -order 5
421: test:
422: suffix: 9
423: args: -num_fields 1 -num_components 2 -order 2
424: test:
425: suffix: 10
426: args: -num_fields 2 -num_components 1,1 -order 2,2
427: test:
428: suffix: 11
429: args: -num_fields 2 -num_components 1,2 -order 2,3
431: TEST*/