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