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