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  dim; /* Topological problem dimension */
      7:   PetscInt  Nf;  /* Number of fields */
      8:   PetscInt *Nc;  /* Number of components per field */
      9:   PetscInt *k;   /* Spectral order per field */
     10: } AppCtx;

     12: static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
     13: {
     14:   PetscInt       len;
     15:   PetscBool      flg;

     19:   options->dim = 2;
     20:   options->Nf  = 0;
     21:   options->Nc  = NULL;
     22:   options->k   = NULL;

     24:   PetscOptionsBegin(comm, "", "SEM Problem Options", "DMPLEX");
     25:   PetscOptionsRangeInt("-dim", "Problem dimension", "ex6.c", options->dim, &options->dim, NULL,1,3);
     26:   PetscOptionsBoundedInt("-num_fields", "The number of fields", "ex6.c", options->Nf, &options->Nf, NULL,0);
     27:   if (options->Nf) {
     28:     len  = options->Nf;
     29:     PetscMalloc1(len, &options->Nc);
     30:     PetscOptionsIntArray("-num_components", "The number of components per field", "ex6.c", options->Nc, &len, &flg);
     31:     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);
     32:     len  = options->Nf;
     33:     PetscMalloc1(len, &options->k);
     34:     PetscOptionsIntArray("-order", "The spectral order per field", "ex6.c", options->k, &len, &flg);
     35:     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);
     36:   }
     37:   PetscOptionsEnd();
     38:   return(0);
     39: }

     41: static PetscErrorCode LoadData2D(DM dm, PetscInt Ni, PetscInt Nj, PetscInt clSize, Vec u, AppCtx *user)
     42: {
     43:   PetscInt       i, j, f, c;
     45:   PetscScalar *closure;

     48:   PetscMalloc1(clSize,&closure);
     49:   for (j = 0; j < Nj; ++j) {
     50:     for (i = 0; i < Ni; ++i) {
     51:       PetscInt    ki, kj, o = 0;
     52:       PetscArrayzero(closure,clSize);

     54:       for (f = 0; f < user->Nf; ++f) {
     55:         PetscInt ioff = i*user->k[f], joff = j*user->k[f];

     57:         for (kj = 0; kj <= user->k[f]; ++kj) {
     58:           for (ki = 0; ki <= user->k[f]; ++ki) {
     59:             for (c = 0; c < user->Nc[f]; ++c) {
     60:               closure[o++] = ((kj + joff)*(Ni*user->k[f]+1) + ki + ioff)*user->Nc[f]+c;
     61:             }
     62:           }
     63:         }
     64:       }
     65:       DMPlexVecSetClosure(dm, NULL, u, j*Ni+i, closure, INSERT_VALUES);
     66:     }
     67:   }
     68:   PetscFree(closure);
     69:   return(0);
     70: }

     72: static PetscErrorCode LoadData3D(DM dm, PetscInt Ni, PetscInt Nj, PetscInt Nk, PetscInt clSize, Vec u, AppCtx *user)
     73: {
     74:   PetscInt       i, j, k, f, c;
     76:   PetscScalar *closure;

     79:   PetscMalloc1(clSize,&closure);
     80:   for (k = 0; k < Nk; ++k) {
     81:     for (j = 0; j < Nj; ++j) {
     82:       for (i = 0; i < Ni; ++i) {
     83:         PetscInt    ki, kj, kk, o = 0;
     84:         PetscArrayzero(closure,clSize);

     86:         for (f = 0; f < user->Nf; ++f) {
     87:           PetscInt ioff = i*user->k[f], joff = j*user->k[f], koff = k*user->k[f];

     89:           for (kk = 0; kk <= user->k[f]; ++kk) {
     90:             for (kj = 0; kj <= user->k[f]; ++kj) {
     91:               for (ki = 0; ki <= user->k[f]; ++ki) {
     92:                 for (c = 0; c < user->Nc[f]; ++c) {
     93:                   closure[o++] = (((kk + koff)*(Nj*user->k[f]+1) + kj + joff)*(Ni*user->k[f]+1) + ki + ioff)*user->Nc[f]+c;
     94:                 }
     95:               }
     96:             }
     97:           }
     98:         }
     99:         DMPlexVecSetClosure(dm, NULL, u, (k*Nj+j)*Ni+i, closure, INSERT_VALUES);
    100:       }
    101:     }
    102:   }
    103:   PetscFree(closure);
    104:   return(0);
    105: }

    107: static PetscErrorCode CheckPoint(DM dm, Vec u, PetscInt point, AppCtx *user)
    108: {
    109:   PetscSection       s;
    110:   PetscScalar        *a;
    111:   const PetscScalar  *array;
    112:   PetscInt           dof, d;
    113:   PetscErrorCode     ierr;

    116:   DMGetLocalSection(dm, &s);
    117:   VecGetArrayRead(u, &array);
    118:   DMPlexPointLocalRead(dm, point, array, &a);
    119:   PetscSectionGetDof(s, point, &dof);
    120:   PetscPrintf(PETSC_COMM_SELF, "Point %D: ", point);
    121:   for (d = 0; d < dof; ++d) {
    122:     if (d > 0) {PetscPrintf(PETSC_COMM_SELF, ", ");}
    123:     PetscPrintf(PETSC_COMM_SELF, "%2.0f", (double) PetscRealPart(a[d]));
    124:   }
    125:   PetscPrintf(PETSC_COMM_SELF, "\n");
    126:   VecRestoreArrayRead(u, &array);
    127:   return(0);
    128: }

    130: static PetscErrorCode ReadData2D(DM dm, Vec u, AppCtx *user)
    131: {
    132:   PetscInt       cStart, cEnd, cell;

    136:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
    137:   for (cell = cStart; cell < cEnd; ++cell) {
    138:     PetscScalar *closure = NULL;
    139:     PetscInt     closureSize, ki, kj, f, c, foff = 0;

    141:     DMPlexVecGetClosure(dm, NULL, u, cell, &closureSize, &closure);
    142:     PetscPrintf(PETSC_COMM_SELF, "Cell %D\n", cell);
    143:     for (f = 0; f < user->Nf; ++f) {
    144:       PetscPrintf(PETSC_COMM_SELF, "  Field %D\n", f);
    145:       for (kj = user->k[f]; kj >= 0; --kj) {
    146:         for (ki = 0; ki <= user->k[f]; ++ki) {
    147:           if (ki > 0) {PetscPrintf(PETSC_COMM_SELF, "  ");}
    148:           for (c = 0; c < user->Nc[f]; ++c) {
    149:             if (c > 0) PetscPrintf(PETSC_COMM_SELF, ",");
    150:             PetscPrintf(PETSC_COMM_SELF, "%2.0f", (double) PetscRealPart(closure[(kj*(user->k[f]+1) + ki)*user->Nc[f]+c + foff]));
    151:           }
    152:         }
    153:         PetscPrintf(PETSC_COMM_SELF, "\n");
    154:       }
    155:       PetscPrintf(PETSC_COMM_SELF, "\n\n");
    156:       foff += PetscSqr(user->k[f]+1);
    157:     }
    158:     DMPlexVecRestoreClosure(dm, NULL, u, cell, &closureSize, &closure);
    159:     PetscPrintf(PETSC_COMM_SELF, "\n\n");
    160:   }
    161:   return(0);
    162: }

    164: static PetscErrorCode ReadData3D(DM dm, Vec u, AppCtx *user)
    165: {
    166:   PetscInt       cStart, cEnd, cell;

    170:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
    171:   for (cell = cStart; cell < cEnd; ++cell) {
    172:     PetscScalar *closure = NULL;
    173:     PetscInt     closureSize, ki, kj, kk, f, c, foff = 0;

    175:     DMPlexVecGetClosure(dm, NULL, u, cell, &closureSize, &closure);
    176:     PetscPrintf(PETSC_COMM_SELF, "Cell %D\n", cell);
    177:     for (f = 0; f < user->Nf; ++f) {
    178:       PetscPrintf(PETSC_COMM_SELF, "  Field %D\n", f);
    179:       for (kk = user->k[f]; kk >= 0; --kk) {
    180:         for (kj = user->k[f]; kj >= 0; --kj) {
    181:           for (ki = 0; ki <= user->k[f]; ++ki) {
    182:             if (ki > 0) {PetscPrintf(PETSC_COMM_SELF, "  ");}
    183:             for (c = 0; c < user->Nc[f]; ++c) {
    184:               if (c > 0) PetscPrintf(PETSC_COMM_SELF, ",");
    185:               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]));
    186:             }
    187:           }
    188:           PetscPrintf(PETSC_COMM_SELF, "\n");
    189:         }
    190:         PetscPrintf(PETSC_COMM_SELF, "\n");
    191:       }
    192:       PetscPrintf(PETSC_COMM_SELF, "\n\n");
    193:       foff += PetscSqr(user->k[f]+1);
    194:     }
    195:     DMPlexVecRestoreClosure(dm, NULL, u, cell, &closureSize, &closure);
    196:     PetscPrintf(PETSC_COMM_SELF, "\n\n");
    197:   }
    198:   return(0);
    199: }

    201: static PetscErrorCode SetSymmetries(DM dm, PetscSection s, AppCtx *user)
    202: {
    203:   PetscInt       f, o, i, j, k, c, d;
    204:   DMLabel        depthLabel;

    208:   DMGetLabel(dm,"depth",&depthLabel);
    209:   for (f = 0; f < user->Nf; f++) {
    210:     PetscSectionSym sym;

    212:     if (user->k[f] < 3) continue; /* No symmetries needed for order < 3, because no cell, facet, edge or vertex has more than one node */
    213:     PetscSectionSymCreateLabel(PetscObjectComm((PetscObject)s),depthLabel,&sym);

    215:     for (d = 0; d <= user->dim; d++) {
    216:       if (d == 1) {
    217:         PetscInt        numDof  = user->k[f] - 1;
    218:         PetscInt        numComp = user->Nc[f];
    219:         PetscInt        minOrnt = -2;
    220:         PetscInt        maxOrnt = 2;
    221:         PetscInt        **perms;

    223:         PetscCalloc1(maxOrnt - minOrnt,&perms);
    224:         for (o = minOrnt; o < maxOrnt; o++) {
    225:           PetscInt *perm;

    227:           if (o == -1 || !o) { /* identity */
    228:             perms[o - minOrnt] = NULL;
    229:           } else {
    230:             PetscMalloc1(numDof * numComp, &perm);
    231:             for (i = numDof - 1, k = 0; i >= 0; i--) {
    232:               for (j = 0; j < numComp; j++, k++) perm[k] = i * numComp + j;
    233:             }
    234:             perms[o - minOrnt] = perm;
    235:           }
    236:         }
    237:         PetscSectionSymLabelSetStratum(sym,d,numDof*numComp,minOrnt,maxOrnt,PETSC_OWN_POINTER,(const PetscInt **) perms,NULL);
    238:       } else if (d == 2) {
    239:         PetscInt        perEdge = user->k[f] - 1;
    240:         PetscInt        numDof  = perEdge * perEdge;
    241:         PetscInt        numComp = user->Nc[f];
    242:         PetscInt        minOrnt = -4;
    243:         PetscInt        maxOrnt = 4;
    244:         PetscInt        **perms;

    246:         PetscCalloc1(maxOrnt-minOrnt,&perms);
    247:         for (o = minOrnt; o < maxOrnt; o++) {
    248:           PetscInt *perm;

    250:           if (!o) continue; /* identity */
    251:           PetscMalloc1(numDof * numComp, &perm);
    252:           /* We want to perm[k] to list which *localArray* position the *sectionArray* position k should go to for the given orientation*/
    253:           switch (o) {
    254:           case 0:
    255:             break; /* identity */
    256:           case -4: /* flip along (-1,-1)--( 1, 1), which swaps edges 0 and 3 and edges 1 and 2.  This swaps the i and j variables */
    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 * j + i) * numComp + c;
    261:                 }
    262:               }
    263:             }
    264:             break;
    265:           case -3: /* flip along (-1, 0)--( 1, 0), which swaps edges 0 and 2.  This reverses the i variable */
    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 - i) + j) * numComp + c;
    270:                 }
    271:               }
    272:             }
    273:             break;
    274:           case -2: /* 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 */
    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 * (perEdge - 1 - j) + (perEdge - 1 - i)) * numComp + c;
    279:                 }
    280:               }
    281:             }
    282:             break;
    283:           case -1: /* flip along ( 0,-1)--( 0, 1), which swaps edges 3 and 1.  This 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 * i + (perEdge - 1 - j)) * numComp + c;
    288:                 }
    289:               }
    290:             }
    291:             break;
    292:           case  1: /* rotate section edge 1 to local edge 0.  This swaps the i and j variables and then reverses the j variable */
    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 - j) + i) * numComp + c;
    297:                 }
    298:               }
    299:             }
    300:             break;
    301:           case  2: /* rotate section edge 2 to local edge 0.  This reverse both i and j variables */
    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 * (perEdge - 1 - i) + (perEdge - 1 - j)) * numComp + c;
    306:                 }
    307:               }
    308:             }
    309:             break;
    310:           case  3: /* rotate section edge 3 to local edge 0.  This swaps the i and j variables and then reverses the i variable */
    311:             for (i = 0, k = 0; i < perEdge; i++) {
    312:               for (j = 0; j < perEdge; j++, k++) {
    313:                 for (c = 0; c < numComp; c++) {
    314:                   perm[k * numComp + c] = (perEdge * j + (perEdge - 1 - i)) * numComp + c;
    315:                 }
    316:               }
    317:             }
    318:             break;
    319:           default:
    320:             break;
    321:           }
    322:           perms[o - minOrnt] = perm;
    323:         }
    324:         PetscSectionSymLabelSetStratum(sym,d,numDof*numComp,minOrnt,maxOrnt,PETSC_OWN_POINTER,(const PetscInt **) perms,NULL);
    325:       }
    326:     }
    327:     PetscSectionSetFieldSym(s,f,sym);
    328:     PetscSectionSymDestroy(&sym);
    329:   }
    330:   PetscSectionViewFromOptions(s,NULL,"-section_with_sym_view");
    331:   return(0);
    332: }

    334: int main(int argc, char **argv)
    335: {
    336:   DM             dm;
    337:   PetscSection   s;
    338:   Vec            u;
    339:   AppCtx         user;
    340:   PetscInt       cells[3] = {2, 2, 2};
    341:   PetscInt       size = 0, f;

    344:   PetscInitialize(&argc, &argv, NULL,help);if (ierr) return ierr;
    345:   ProcessOptions(PETSC_COMM_WORLD, &user);
    346:   DMPlexCreateBoxMesh(PETSC_COMM_WORLD, user.dim, PETSC_FALSE, cells, NULL, NULL, NULL, PETSC_TRUE, &dm);
    347:   DMSetFromOptions(dm);
    348:   DMViewFromOptions(dm, NULL, "-dm_view");
    349:   /* Create a section for SEM order k */
    350:   {
    351:     PetscInt *numDof, d;

    353:     PetscMalloc1(user.Nf*(user.dim+1), &numDof);
    354:     for (f = 0; f < user.Nf; ++f) {
    355:       for (d = 0; d <= user.dim; ++d) numDof[f*(user.dim+1)+d] = PetscPowInt(user.k[f]-1, d)*user.Nc[f];
    356:       size += PetscPowInt(user.k[f]+1, d)*user.Nc[f];
    357:     }
    358:     DMSetNumFields(dm, user.Nf);
    359:     DMPlexCreateSection(dm, NULL, user.Nc, numDof, 0, NULL, NULL, NULL, NULL, &s);
    360:     SetSymmetries(dm, s, &user);
    361:     PetscFree(numDof);
    362:   }
    363:   DMSetLocalSection(dm, s);
    364:   /* Create spectral ordering and load in data */
    365:   DMPlexSetClosurePermutationTensor(dm, PETSC_DETERMINE, NULL);
    366:   DMGetLocalVector(dm, &u);
    367:   switch (user.dim) {
    368:   case 2: LoadData2D(dm, 2, 2, size, u, &user);break;
    369:   case 3: LoadData3D(dm, 2, 2, 2, size, u, &user);break;
    370:   }
    371:   /* Remove ordering and check some values */
    372:   PetscSectionSetClosurePermutation(s, (PetscObject) dm, user.dim, NULL);
    373:   switch (user.dim) {
    374:   case 2:
    375:     CheckPoint(dm, u,  0, &user);
    376:     CheckPoint(dm, u, 13, &user);
    377:     CheckPoint(dm, u, 15, &user);
    378:     CheckPoint(dm, u, 19, &user);
    379:     break;
    380:   case 3:
    381:     CheckPoint(dm, u,  0, &user);
    382:     CheckPoint(dm, u, 13, &user);
    383:     CheckPoint(dm, u, 15, &user);
    384:     CheckPoint(dm, u, 19, &user);
    385:     break;
    386:   }
    387:   /* Recreate spectral ordering and read out data */
    388:   DMPlexSetClosurePermutationTensor(dm, PETSC_DETERMINE, s);
    389:   switch (user.dim) {
    390:   case 2: ReadData2D(dm, u, &user);break;
    391:   case 3: ReadData3D(dm, u, &user);break;
    392:   }
    393:   DMRestoreLocalVector(dm, &u);
    394:   PetscSectionDestroy(&s);
    395:   DMDestroy(&dm);
    396:   PetscFree(user.Nc);
    397:   PetscFree(user.k);
    398:   PetscFinalize();
    399:   return ierr;
    400: }

    402: /*TEST

    404:   # Spectral ordering 2D 0-5
    405:   test:
    406:     suffix: 0
    407:     args: -dim 2 -num_fields 1 -num_components 1 -order 2
    408:   test:
    409:     suffix: 1
    410:     args: -dim 2 -num_fields 1 -num_components 1 -order 3
    411:   test:
    412:     suffix: 2
    413:     args: -dim 2 -num_fields 1 -num_components 1 -order 5
    414:   test:
    415:     suffix: 3
    416:     args: -dim 2 -num_fields 1 -num_components 2 -order 2
    417:   test:
    418:     suffix: 4
    419:     args: -dim 2 -num_fields 2 -num_components 1,1 -order 2,2
    420:   test:
    421:     suffix: 5
    422:     args: -dim 2 -num_fields 2 -num_components 1,2 -order 2,3
    423:   # Spectral ordering 3D 6-11
    424:   test:
    425:     suffix: 6
    426:     args: -dim 3 -num_fields 1 -num_components 1 -order 2
    427:   test:
    428:     suffix: 7
    429:     args: -dim 3 -num_fields 1 -num_components 1 -order 3
    430:   test:
    431:     suffix: 8
    432:     args: -dim 3 -num_fields 1 -num_components 1 -order 5
    433:   test:
    434:     suffix: 9
    435:     args: -dim 3 -num_fields 1 -num_components 2 -order 2
    436:   test:
    437:     suffix: 10
    438:     args: -dim 3 -num_fields 2 -num_components 1,1 -order 2,2
    439:   test:
    440:     suffix: 11
    441:     args: -dim 3 -num_fields 2 -num_components 1,2 -order 2,3

    443: TEST*/