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

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

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

 20:   bx = DM_BOUNDARY_NONE;
 21:   by = DM_BOUNDARY_NONE;
 22:   bz = DM_BOUNDARY_NONE;

 24:   periodicity = 0;

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

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

 43:   DMSetFromOptions(daf);

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

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

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

 58:     DMGetCoordinateDM(dac,&cdac);
 59:     DMGetCoordinateDM(daf,&cdaf);

 61:     DMGetCoordinates(dac,&coordsc);
 62:     DMGetCoordinates(daf,&coordsf);

 64:     DMCreateInjection(cdac,cdaf,&inject);
 65:     MatScatterGetVecScatter(inject,&vscat);
 66:     VecScatterBegin(vscat,coordsf,coordsc,INSERT_VALUES,SCATTER_FORWARD);
 67:     VecScatterEnd(vscat  ,coordsf,coordsc,INSERT_VALUES,SCATTER_FORWARD);
 68:     MatDestroy(&inject);

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

 81:   if (0) {
 82:     DMCreateGlobalVector(dac,&ac);
 83:     VecZeroEntries(ac);

 85:     DMCreateGlobalVector(daf,&af);
 86:     VecZeroEntries(af);

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

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

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

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

122:   test1_DAInjection3d(mx,my,mz);

124:   PetscFinalize();
125:   return 0;
126: }