Actual source code: ex21.c
petsc-3.4.2 2013-07-02

      1: static const char help[] = "Test DMCreateInjection() for mapping coordinates in 3D";

      3: #include <petscvec.h>
      4: #include <petscmat.h>
      5: #include <petscdmda.h>

      9: PetscErrorCode test1_DAInjection3d(PetscInt mx, PetscInt my, PetscInt mz)
     10: {
     11:   PetscErrorCode   ierr;
     12:   DM               dac,daf;
     13:   PetscViewer      vv;
     14:   Vec              ac,af;
     15:   PetscInt         periodicity;
     16:   DMDABoundaryType bx,by,bz;

     19:   bx = DMDA_BOUNDARY_NONE;
     20:   by = DMDA_BOUNDARY_NONE;
     21:   bz = DMDA_BOUNDARY_NONE;

     23:   periodicity = 0;

     25:   PetscOptionsGetInt(NULL,"-periodic", &periodicity, NULL);
     26:   if (periodicity==1) {
     27:     bx = DMDA_BOUNDARY_PERIODIC;
     28:   } else if (periodicity==2) {
     29:     by = DMDA_BOUNDARY_PERIODIC;
     30:   } else if (periodicity==3) {
     31:     bz = DMDA_BOUNDARY_PERIODIC;
     32:   }

     34:   DMDACreate3d(PETSC_COMM_WORLD, bx,by,bz, DMDA_STENCIL_BOX,
     35:                       mx+1, my+1,mz+1,
     36:                       PETSC_DECIDE, PETSC_DECIDE,PETSC_DECIDE,
     37:                       1, /* 1 dof */
     38:                       1, /* stencil = 1 */
     39:                       NULL,NULL,NULL,
     40:                       &daf);

     42:   DMSetFromOptions(daf);

     44:   DMCoarsen(daf,MPI_COMM_NULL,&dac);

     46:   DMDASetUniformCoordinates(dac, -1.0,1.0, -1.0,1.0, -1.0,1.0);
     47:   DMDASetUniformCoordinates(daf, -1.0,1.0, -1.0,1.0, -1.0,1.0);

     49:   {
     50:     DM         cdaf,cdac;
     51:     Vec        coordsc,coordsf,coordsf2;
     52:     VecScatter inject;
     53:     Mat        interp;
     54:     PetscReal  norm;

     56:     DMGetCoordinateDM(dac,&cdac);
     57:     DMGetCoordinateDM(daf,&cdaf);

     59:     DMGetCoordinates(dac,&coordsc);
     60:     DMGetCoordinates(daf,&coordsf);

     62:     DMCreateInjection(cdac,cdaf,&inject);

     64:     VecScatterBegin(inject,coordsf,coordsc,INSERT_VALUES,SCATTER_FORWARD);
     65:     VecScatterEnd(inject  ,coordsf,coordsc,INSERT_VALUES,SCATTER_FORWARD);
     66:     VecScatterDestroy(&inject);

     68:     DMCreateInterpolation(cdac,cdaf,&interp,NULL);
     69:     VecDuplicate(coordsf,&coordsf2);
     70:     MatInterpolate(interp,coordsc,coordsf2);
     71:     VecAXPY(coordsf2,-1.0,coordsf);
     72:     VecNorm(coordsf2,NORM_MAX,&norm);
     73:     /* The fine coordinates are only reproduced in certain cases */
     74:     if (!bx && !by && !bz && norm > 1.e-10) {PetscPrintf(PETSC_COMM_WORLD,"Norm %G\n",norm);}
     75:     VecDestroy(&coordsf2);
     76:     MatDestroy(&interp);
     77:   }

     79:   if (0) {
     80:     DMCreateGlobalVector(dac,&ac);
     81:     VecZeroEntries(ac);

     83:     DMCreateGlobalVector(daf,&af);
     84:     VecZeroEntries(af);

     86:     PetscViewerASCIIOpen(PETSC_COMM_WORLD, "dac_7.vtk", &vv);
     87:     PetscViewerSetFormat(vv, PETSC_VIEWER_ASCII_VTK);
     88:     DMView(dac, vv);
     89:     VecView(ac, vv);
     90:     PetscViewerDestroy(&vv);

     92:     PetscViewerASCIIOpen(PETSC_COMM_WORLD, "daf_7.vtk", &vv);
     93:     PetscViewerSetFormat(vv, PETSC_VIEWER_ASCII_VTK);
     94:     DMView(daf, vv);
     95:     VecView(af, vv);
     96:     PetscViewerDestroy(&vv);
     97:     VecDestroy(&ac);
     98:     VecDestroy(&af);
     99:   }
    100:   DMDestroy(&dac);
    101:   DMDestroy(&daf);
    102:   return(0);
    103: }

    107: int main(int argc,char **argv)
    108: {
    110:   PetscInt       mx,my,mz;

    112:   PetscInitialize(&argc,&argv,0,help);
    113:   mx   = 2;
    114:   my   = 2;
    115:   mz   = 2;
    116:   PetscOptionsGetInt(NULL,"-mx", &mx, 0);
    117:   PetscOptionsGetInt(NULL,"-my", &my, 0);
    118:   PetscOptionsGetInt(NULL,"-mz", &mz, 0);

    120:   test1_DAInjection3d(mx,my,mz);

    122:   PetscFinalize();
    123:   return 0;
    124: }