Actual source code: ex21.c

  1: static const char help[] = "Test DMGetInjection 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: {
 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;
 22:   periodicity = 0;
 23:   PetscOptionsGetInt(PETSC_NULL,"-periodic", &periodicity, PETSC_NULL);
 24:   if (periodicity==1) {
 25:     bx = DMDA_BOUNDARY_PERIODIC;
 26:   } else if (periodicity==2) {
 27:     by = DMDA_BOUNDARY_PERIODIC;
 28:   } else if (periodicity==3) {
 29:     bz = DMDA_BOUNDARY_PERIODIC;
 30:   }

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

 40:   DMSetFromOptions(daf);

 42:   DMCoarsen(daf,PETSC_NULL,&dac);

 44:   DMDASetUniformCoordinates(dac, -1.0,1.0, -1.0,1.0, -1.0,1.0 );
 45:   DMDASetUniformCoordinates(daf, -1.0,1.0, -1.0,1.0, -1.0,1.0 );

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

 54:     DMDAGetCoordinateDA(dac,&cdac);
 55:     DMDAGetCoordinateDA(daf,&cdaf);

 57:     DMDAGetCoordinates(dac,&coordsc);
 58:     DMDAGetCoordinates(daf,&coordsf);

 60:     DMGetInjection(cdac,cdaf,&inject);

 62:     VecScatterBegin(inject,coordsf,coordsc,INSERT_VALUES,SCATTER_FORWARD);
 63:     VecScatterEnd(inject  ,coordsf,coordsc,INSERT_VALUES,SCATTER_FORWARD);
 64:     VecScatterDestroy(&inject);

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

 77:   if (0) {
 78:     DMCreateGlobalVector(dac,&ac);
 79:     VecZeroEntries(ac);

 81:     DMCreateGlobalVector(daf,&af);
 82:     VecZeroEntries(af);

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

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

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

110:   PetscInitialize(&argc,&argv,0,0);
111:   mx = 2;
112:   my = 2;
113:   mz = 2;
114:   PetscOptionsGetInt(PETSC_NULL,"-mx", &mx, 0 );
115:   PetscOptionsGetInt(PETSC_NULL,"-my", &my, 0 );
116:   PetscOptionsGetInt(PETSC_NULL,"-mz", &mz, 0 );

118:   test1_DAInjection3d(mx,my,mz);

120:   PetscFinalize();
121:   return 0;
122: }