Actual source code: ex1.c

      1: static char help[] = "Demonstration of creating and viewing DMFields objects.\n\n";

      3: #include <petscdmfield.h>
      4: #include <petscdmplex.h>
      5: #include <petscdmda.h>

      7: static PetscErrorCode ViewResults(PetscViewer viewer, PetscInt N, PetscInt dim, PetscScalar *B, PetscScalar *D, PetscScalar *H, PetscReal *rB, PetscReal *rD, PetscReal *rH)
      8: {

     12:   PetscViewerASCIIPrintf(viewer,"B:\n");
     13:   PetscScalarView(N,B,viewer);
     14:   PetscViewerASCIIPrintf(viewer,"D:\n");
     15:   PetscScalarView(N*dim,D,viewer);
     16:   PetscViewerASCIIPrintf(viewer,"H:\n");
     17:   PetscScalarView(N*dim*dim,H,viewer);

     19:   PetscViewerASCIIPrintf(viewer,"rB:\n");
     20:   PetscRealView(N,rB,viewer);
     21:   PetscViewerASCIIPrintf(viewer,"rD:\n");
     22:   PetscRealView(N*dim,rD,viewer);
     23:   PetscViewerASCIIPrintf(viewer,"rH:\n");
     24:   PetscRealView(N*dim*dim,rH,viewer);
     25:   return(0);
     26: }

     28: static PetscErrorCode TestEvaluate(DMField field, PetscInt n, PetscRandom rand)
     29: {
     30:   DM             dm;
     31:   PetscInt       dim, i, nc;
     32:   PetscScalar    *B, *D, *H;
     33:   PetscReal      *rB, *rD, *rH;
     34:   Vec            points;
     35:   PetscScalar    *array;
     36:   PetscViewer    viewer;
     37:   MPI_Comm       comm;

     41:   comm = PetscObjectComm((PetscObject)field);
     42:   DMFieldGetNumComponents(field,&nc);
     43:   DMFieldGetDM(field,&dm);
     44:   DMGetDimension(dm,&dim);
     45:   VecCreateMPI(PetscObjectComm((PetscObject)field),n * dim,PETSC_DETERMINE,&points);
     46:   VecSetBlockSize(points,dim);
     47:   VecGetArray(points,&array);
     48:   for (i = 0; i < n * dim; i++) {PetscRandomGetValue(rand,&array[i]);}
     49:   VecRestoreArray(points,&array);
     50:   PetscMalloc6(n*nc,&B,n*nc,&rB,n*nc*dim,&D,n*nc*dim,&rD,n*nc*dim*dim,&H,n*nc*dim*dim,&rH);
     51:   DMFieldEvaluate(field,points,PETSC_SCALAR,B,D,H);
     52:   DMFieldEvaluate(field,points,PETSC_REAL,rB,rD,rH);
     53:   viewer = PETSC_VIEWER_STDOUT_(comm);

     55:   PetscObjectSetName((PetscObject)points,"Test Points");
     56:   VecView(points,viewer);
     57:   ViewResults(viewer,n*nc,dim,B,D,H,rB,rD,rH);

     59:   PetscFree6(B,rB,D,rD,H,rH);
     60:   VecDestroy(&points);
     61:   return(0);
     62: }

     64: static PetscErrorCode TestEvaluateFE(DMField field, PetscInt n, PetscInt cStart, PetscInt cEnd, PetscQuadrature quad, PetscRandom rand)
     65: {
     66:   DM             dm;
     67:   PetscInt       dim, i, nc, nq;
     68:   PetscInt       N;
     69:   PetscScalar    *B, *D, *H;
     70:   PetscReal      *rB, *rD, *rH;
     71:   PetscInt       *cells;
     72:   IS             cellIS;
     73:   PetscViewer    viewer;
     74:   MPI_Comm       comm;

     78:   comm = PetscObjectComm((PetscObject)field);
     79:   DMFieldGetNumComponents(field,&nc);
     80:   DMFieldGetDM(field,&dm);
     81:   DMGetDimension(dm,&dim);
     82:   PetscRandomSetInterval(rand,(PetscScalar) cStart, (PetscScalar) cEnd);
     83:   PetscMalloc1(n,&cells);
     84:   for (i = 0; i < n; i++) {
     85:     PetscReal rc;

     87:     PetscRandomGetValueReal(rand,&rc);
     88:     cells[i] = PetscFloorReal(rc);
     89:   }
     90:   ISCreateGeneral(comm,n,cells,PETSC_OWN_POINTER,&cellIS);
     91:   PetscObjectSetName((PetscObject)cellIS,"FE Test Cells");
     92:   PetscQuadratureGetData(quad,NULL,NULL,&nq,NULL,NULL);
     93:   N    = n * nq * nc;
     94:   PetscMalloc6(N,&B,N,&rB,N*dim,&D,N*dim,&rD,N*dim*dim,&H,N*dim*dim,&rH);
     95:   DMFieldEvaluateFE(field,cellIS,quad,PETSC_SCALAR,B,D,H);
     96:   DMFieldEvaluateFE(field,cellIS,quad,PETSC_REAL,rB,rD,rH);
     97:   viewer = PETSC_VIEWER_STDOUT_(comm);

     99:   PetscObjectSetName((PetscObject)quad,"Test quadrature");
    100:   PetscQuadratureView(quad,viewer);
    101:   ISView(cellIS,viewer);
    102:   ViewResults(viewer,N,dim,B,D,H,rB,rD,rH);

    104:   PetscFree6(B,rB,D,rD,H,rH);
    105:   ISDestroy(&cellIS);
    106:   return(0);
    107: }

    109: static PetscErrorCode TestEvaluateFV(DMField field, PetscInt n, PetscInt cStart, PetscInt cEnd, PetscRandom rand)
    110: {
    111:   DM             dm;
    112:   PetscInt       dim, i, nc;
    113:   PetscInt       N;
    114:   PetscScalar    *B, *D, *H;
    115:   PetscReal      *rB, *rD, *rH;
    116:   PetscInt       *cells;
    117:   IS             cellIS;
    118:   PetscViewer    viewer;
    119:   MPI_Comm       comm;

    123:   comm = PetscObjectComm((PetscObject)field);
    124:   DMFieldGetNumComponents(field,&nc);
    125:   DMFieldGetDM(field,&dm);
    126:   DMGetDimension(dm,&dim);
    127:   PetscRandomSetInterval(rand,(PetscScalar) cStart, (PetscScalar) cEnd);
    128:   PetscMalloc1(n,&cells);
    129:   for (i = 0; i < n; i++) {
    130:     PetscReal rc;

    132:     PetscRandomGetValueReal(rand,&rc);
    133:     cells[i] = PetscFloorReal(rc);
    134:   }
    135:   ISCreateGeneral(comm,n,cells,PETSC_OWN_POINTER,&cellIS);
    136:   PetscObjectSetName((PetscObject)cellIS,"FV Test Cells");
    137:   N    = n * nc;
    138:   PetscMalloc6(N,&B,N,&rB,N*dim,&D,N*dim,&rD,N*dim*dim,&H,N*dim*dim,&rH);
    139:   DMFieldEvaluateFV(field,cellIS,PETSC_SCALAR,B,D,H);
    140:   DMFieldEvaluateFV(field,cellIS,PETSC_REAL,rB,rD,rH);
    141:   viewer = PETSC_VIEWER_STDOUT_(comm);

    143:   ISView(cellIS,viewer);
    144:   ViewResults(viewer,N,dim,B,D,H,rB,rD,rH);

    146:   PetscFree6(B,rB,D,rD,H,rH);
    147:   ISDestroy(&cellIS);
    148:   return(0);
    149: }

    151: static PetscErrorCode radiusSquared(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar u[], void *ctx)
    152: {
    153:   PetscInt       i;
    154:   PetscReal      r2 = 0.;

    157:   for (i = 0; i < dim; i++) {r2 += PetscSqr(x[i]);}
    158:   for (i = 0; i < Nf; i++) {
    159:     u[i] = (i + 1) * r2;
    160:   }
    161:   return(0);
    162: }

    164: static PetscErrorCode TestShellEvaluate(DMField field, Vec points, PetscDataType type, void *B, void *D, void *H)
    165: {
    166:   Vec                ctxVec = NULL;
    167:   const PetscScalar *mult;
    168:   PetscInt           dim;
    169:   const PetscScalar *x;
    170:   PetscInt           Nc, n, i, j, k, l;
    171:   PetscErrorCode     ierr;

    174:   DMFieldGetNumComponents(field, &Nc);
    175:   DMFieldShellGetContext(field, (void *) &ctxVec);
    176:   VecGetBlockSize(points, &dim);
    177:   VecGetLocalSize(points, &n);
    178:   n /= Nc;
    179:   VecGetArrayRead(ctxVec, &mult);
    180:   VecGetArrayRead(points, &x);
    181:   for (i = 0; i < n; i++) {
    182:     PetscReal r2 = 0.;

    184:     for (j = 0; j < dim; j++) {r2 += PetscSqr(PetscRealPart(x[i * dim + j]));}
    185:     for (j = 0; j < Nc; j++) {
    186:       PetscReal m = PetscRealPart(mult[j]);
    187:       if (B) {
    188:         if (type == PETSC_SCALAR) {
    189:           ((PetscScalar *)B)[i * Nc + j] = m * r2;
    190:         } else {
    191:           ((PetscReal *)B)[i * Nc + j] = m * r2;
    192:         }
    193:       }
    194:       if (D) {
    195:         if (type == PETSC_SCALAR) {
    196:           for (k = 0; k < dim; k++) ((PetscScalar *)D)[(i * Nc + j) * dim + k] = 2. * m * x[i * dim + k];
    197:         } else {
    198:           for (k = 0; k < dim; k++) ((PetscReal   *)D)[(i * Nc + j) * dim + k] = 2. * m * PetscRealPart(x[i * dim + k]);
    199:         }
    200:       }
    201:       if (H) {
    202:         if (type == PETSC_SCALAR) {
    203:           for (k = 0; k < dim; k++) for (l = 0; l < dim; l++) ((PetscScalar *)H)[((i * Nc + j) * dim + k) * dim + l] = (k == l) ? 2. * m : 0.;
    204:         } else {
    205:           for (k = 0; k < dim; k++) for (l = 0; l < dim; l++) ((PetscReal   *)H)[((i * Nc + j) * dim + k) * dim + l] = (k == l) ? 2. * m : 0.;
    206:         }
    207:       }
    208:     }
    209:   }
    210:   VecRestoreArrayRead(points, &x);
    211:   VecRestoreArrayRead(ctxVec, &mult);
    212:   return(0);
    213: }

    215: static PetscErrorCode TestShellDestroy(DMField field)
    216: {
    217:   Vec                ctxVec = NULL;
    218:   PetscErrorCode     ierr;

    221:   DMFieldShellGetContext(field, (void *) &ctxVec);
    222:   VecDestroy(&ctxVec);
    223:   return(0);
    224: }

    226: int main(int argc, char **argv)
    227: {
    228:   DM              dm = NULL;
    229:   MPI_Comm        comm;
    230:   char            type[256] = DMPLEX;
    231:   PetscBool       isda, isplex;
    232:   PetscInt        dim = 2;
    233:   DMField         field = NULL;
    234:   PetscInt        nc = 1;
    235:   PetscInt        cStart = -1, cEnd = -1;
    236:   PetscRandom     rand;
    237:   PetscQuadrature quad = NULL;
    238:   PetscInt        pointsPerEdge = 2;
    239:   PetscInt        numPoint = 0, numFE = 0, numFV = 0;
    240:   PetscBool       testShell = PETSC_FALSE;
    241:   PetscErrorCode  ierr;

    243:   PetscInitialize(&argc, &argv, NULL, help);if (ierr) return ierr;
    244:   comm = PETSC_COMM_WORLD;
    245:   PetscOptionsBegin(comm, "", "DMField Tutorial Options", "DM");
    246:   PetscOptionsFList("-dm_type","DM implementation on which to define field","ex1.c",DMList,type,type,256,NULL);
    247:   PetscOptionsRangeInt("-dim","DM intrinsic dimension", "ex1.c", dim, &dim, NULL,1,3);
    248:   PetscOptionsBoundedInt("-num_components","Number of components in field", "ex1.c", nc, &nc, NULL,1);
    249:   PetscOptionsBoundedInt("-num_quad_points","Number of quadrature points per dimension", "ex1.c", pointsPerEdge, &pointsPerEdge, NULL,1);
    250:   PetscOptionsBoundedInt("-num_point_tests", "Number of test points for DMFieldEvaluate()", "ex1.c", numPoint, &numPoint, NULL,0);
    251:   PetscOptionsBoundedInt("-num_fe_tests", "Number of test cells for DMFieldEvaluateFE()", "ex1.c", numFE, &numFE, NULL,0);
    252:   PetscOptionsBoundedInt("-num_fv_tests", "Number of test cells for DMFieldEvaluateFV()", "ex1.c", numFV, &numFV, NULL,0);
    253:   PetscOptionsBool("-test_shell", "Test the DMFIELDSHELL implementation of DMField", "ex1.c", testShell, &testShell, NULL);
    254:   PetscOptionsEnd();

    256:   if (dim > 3) SETERRQ1(comm,PETSC_ERR_ARG_OUTOFRANGE,"This examples works for dim <= 3, not %D",dim);
    257:   PetscStrncmp(type,DMPLEX,256,&isplex);
    258:   PetscStrncmp(type,DMDA,256,&isda);

    260:   PetscRandomCreate(PETSC_COMM_SELF,&rand);
    261:   PetscRandomSetFromOptions(rand);
    262:   if (isplex) {
    263:     PetscBool simplex = PETSC_TRUE;
    264:     PetscInt  overlap = 0;
    265:     Vec       fieldvec;
    266:     PetscBool interpolate = PETSC_TRUE;
    267:     PetscInt       inttypes[3];
    268:     DMBoundaryType types[3] = {DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE};
    269:     PetscInt       cells[3] = {3,3,3};
    270:     PetscBool      flags[3] = {PETSC_FALSE, PETSC_FALSE, PETSC_FALSE};
    271:     PetscInt       n = 3, i;
    272:     PetscFE   fe;

    274:     PetscOptionsBegin(comm, "", "DMField DMPlex Options", "DM");
    275:     PetscOptionsBool("-simplex","Create a simplicial DMPlex","ex1.c",simplex,&simplex,NULL);
    276:     PetscOptionsBoundedInt("-overlap","DMPlex parallel overlap","ex1.c",overlap,&overlap,NULL,0);
    277:     PetscOptionsBool("-interpolate","Interpolate the DMPlex","ex1.c",interpolate,&interpolate,NULL);
    278:     PetscOptionsEList("-boundary_x","type of boundary in x direction","ex1.c",DMBoundaryTypes,5,DMBoundaryTypes[types[0]],&inttypes[0],&flags[0]);
    279:     PetscOptionsEList("-boundary_y","type of boundary in y direction","ex1.c",DMBoundaryTypes,5,DMBoundaryTypes[types[1]],&inttypes[1],&flags[1]);
    280:     PetscOptionsEList("-boundary_z","type of boundary in z direction","ex1.c",DMBoundaryTypes,5,DMBoundaryTypes[types[2]],&inttypes[2],&flags[2]);
    281:     PetscOptionsIntArray("-cells","Cells per dimension","ex1.c",cells,&n,NULL);
    282:     PetscOptionsEnd();
    283:     for (i = 0; i < 3; i++) {
    284:       if (flags[i]) {types[i] = (DMBoundaryType) inttypes[i];}
    285:     }
    286:     DMPlexCreateBoxMesh(comm,dim,simplex,cells,NULL,NULL,types,interpolate,&dm);
    287:     if (simplex) {
    288:       PetscDTStroudConicalQuadrature(dim, 1, pointsPerEdge, -1.0, 1.0, &quad);
    289:     } else {
    290:       PetscDTGaussTensorQuadrature(dim, 1, pointsPerEdge, -1.0, 1.0, &quad);
    291:     }
    292:     {
    293:       PetscPartitioner part;
    294:       DM               dmDist;

    296:       DMPlexGetPartitioner(dm,&part);
    297:       PetscPartitionerSetType(part,PETSCPARTITIONERSIMPLE);
    298:       DMPlexDistribute(dm,overlap,NULL,&dmDist);
    299:       if (dmDist) {
    300:         DMDestroy(&dm);
    301:         dm = dmDist;
    302:       }
    303:     }
    304:     DMPlexGetHeightStratum(dm,0,&cStart,&cEnd);
    305:     if (testShell) {
    306:       Vec ctxVec;
    307:       PetscInt i;
    308:       PetscScalar *array;

    310:       VecCreateSeq(PETSC_COMM_SELF, nc, &ctxVec);
    311:       VecSetUp(ctxVec);
    312:       VecGetArray(ctxVec,&array);
    313:       for (i = 0; i < nc; i++) array[i] = i + 1.;
    314:       VecRestoreArray(ctxVec,&array);
    315:       DMFieldCreateShell(dm, nc, DMFIELD_VERTEX, (void *) ctxVec, &field);
    316:       DMFieldShellSetEvaluate(field, TestShellEvaluate);
    317:       DMFieldShellSetDestroy(field, TestShellDestroy);
    318:     } else {
    319:       PetscFECreateDefault(comm,dim,nc,simplex,NULL,PETSC_DEFAULT,&fe);
    320:       PetscFESetName(fe,"MyPetscFE");
    321:       DMSetField(dm,0,NULL,(PetscObject)fe);
    322:       DMCreateDS(dm);
    323:       PetscFEDestroy(&fe);
    324:       DMCreateLocalVector(dm,&fieldvec);
    325:       {
    326:         PetscErrorCode (*func[1]) (PetscInt,PetscReal,const PetscReal [],PetscInt, PetscScalar *,void *);
    327:         void            *ctxs[1];

    329:         func[0] = radiusSquared;
    330:         ctxs[0] = NULL;

    332:         DMProjectFunctionLocal(dm,0.0,func,ctxs,INSERT_ALL_VALUES,fieldvec);
    333:       }
    334:       DMFieldCreateDS(dm,0,fieldvec,&field);
    335:       VecDestroy(&fieldvec);
    336:     }
    337:   } else if (isda) {
    338:     PetscInt       i;
    339:     PetscScalar    *cv;

    341:     switch (dim) {
    342:     case 1:
    343:       DMDACreate1d(comm, DM_BOUNDARY_NONE, 3, 1, 1, NULL, &dm);
    344:       break;
    345:     case 2:
    346:       DMDACreate2d(comm, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_BOX, 3, 3, PETSC_DETERMINE, PETSC_DETERMINE, 1, 1, NULL, NULL, &dm);
    347:       break;
    348:     default:
    349:       DMDACreate3d(comm, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_BOX, 3, 3, 3, PETSC_DETERMINE, PETSC_DETERMINE, PETSC_DETERMINE, 1, 1, NULL, NULL, NULL, &dm);
    350:       break;
    351:     }
    352:     DMSetUp(dm);
    353:     DMDAGetHeightStratum(dm,0,&cStart,&cEnd);
    354:     PetscMalloc1(nc * (1 << dim),&cv);
    355:     for (i = 0; i < nc * (1 << dim); i++) {
    356:       PetscReal rv;

    358:       PetscRandomGetValueReal(rand,&rv);
    359:       cv[i] = rv;
    360:     }
    361:     DMFieldCreateDA(dm,nc,cv,&field);
    362:     PetscFree(cv);
    363:     PetscDTGaussTensorQuadrature(dim, 1, pointsPerEdge, -1.0, 1.0, &quad);
    364:   } else SETERRQ1(comm,PETSC_ERR_SUP,"This test does not run for DM type %s",type);

    366:   PetscObjectSetName((PetscObject)dm,"mesh");
    367:   DMViewFromOptions(dm,NULL,"-dm_view");
    368:   DMDestroy(&dm);
    369:   PetscObjectSetName((PetscObject)field,"field");
    370:   PetscObjectViewFromOptions((PetscObject)field,NULL,"-dmfield_view");
    371:   if (numPoint) {TestEvaluate(field,numPoint,rand);}
    372:   if (numFE) {TestEvaluateFE(field,numFE,cStart,cEnd,quad,rand);}
    373:   if (numFV) {TestEvaluateFV(field,numFV,cStart,cEnd,rand);}
    374:   PetscQuadratureDestroy(&quad);
    375:   DMFieldDestroy(&field);
    376:   PetscRandomDestroy(&rand);
    377:   PetscFinalize();
    378:   return ierr;
    379: }

    381: /*TEST

    383:   test:
    384:     suffix: da
    385:     requires: !complex
    386:     args: -dm_type da -dim 2 -num_components 2 -num_point_tests 2 -num_fe_tests 2 -num_fv_tests 2 -dmfield_view

    388:   test:
    389:     suffix: da_1
    390:     requires: !complex
    391:     args: -dm_type da -dim 1  -num_fe_tests 2

    393:   test:
    394:     suffix: da_2
    395:     requires: !complex
    396:     args: -dm_type da -dim 2  -num_fe_tests 2

    398:   test:
    399:     suffix: da_3
    400:     requires: !complex
    401:     args: -dm_type da -dim 3  -num_fe_tests 2

    403:   test:
    404:     suffix: ds
    405:     requires: !complex triangle
    406:     args: -dm_type plex -dim 2 -num_components 2 -num_point_tests 2 -num_fe_tests 2 -num_fv_tests 2 -dmfield_view -petscspace_degree 2 -num_quad_points 1

    408:   test:
    409:     suffix: ds_simplex_0
    410:     requires: !complex triangle
    411:     args: -dm_type plex -dim 2  -num_fe_tests 2  -petscspace_degree 0

    413:   test:
    414:     suffix: ds_simplex_1
    415:     requires: !complex triangle
    416:     args: -dm_type plex -dim 2  -num_fe_tests 2  -petscspace_degree 1

    418:   test:
    419:     suffix: ds_simplex_2
    420:     requires: !complex triangle
    421:     args: -dm_type plex -dim 2  -num_fe_tests 2  -petscspace_degree 2

    423:   test:
    424:     suffix: ds_tensor_2_0
    425:     requires: !complex
    426:     args: -dm_type plex -dim 2  -num_fe_tests 2  -petscspace_poly_tensor 1 -petscspace_degree 0 -simplex 0

    428:   test:
    429:     suffix: ds_tensor_2_1
    430:     requires: !complex
    431:     args: -dm_type plex -dim 2  -num_fe_tests 2  -petscspace_poly_tensor 1 -petscspace_degree 1 -simplex 0

    433:   test:
    434:     suffix: ds_tensor_2_2
    435:     requires: !complex
    436:     args: -dm_type plex -dim 2  -num_fe_tests 2  -petscspace_poly_tensor 1 -petscspace_degree 2 -simplex 0

    438:   test:
    439:     suffix: ds_tensor_3_0
    440:     requires: !complex
    441:     args: -dm_type plex -dim 3  -num_fe_tests 2  -petscspace_poly_tensor 1 -petscspace_degree 0 -simplex 0

    443:   test:
    444:     suffix: ds_tensor_3_1
    445:     requires: !complex
    446:     args: -dm_type plex -dim 3  -num_fe_tests 2  -petscspace_poly_tensor 1 -petscspace_degree 1 -simplex 0

    448:   test:
    449:     suffix: ds_tensor_3_2
    450:     requires: !complex
    451:     args: -dm_type plex -dim 3  -num_fe_tests 2  -petscspace_poly_tensor 1 -petscspace_degree 2 -simplex 0

    453:   test:
    454:     suffix: shell
    455:     requires: !complex triangle
    456:     args: -dm_type plex -dim 2 -num_components 2 -num_point_tests 2 -num_fe_tests 2 -num_fv_tests 2 -dmfield_view -num_quad_points 1 -test_shell

    458: TEST*/