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*/