1: static char help[] = "Tests coarsening with DM.\n";



  5:  #include <petscsys.h>
  6:  #include <petscvec.h>
  7:  #include <petscdmda.h>


 10: int main(int argc, char **argv)
 11: {
 13:   Vec            x,yp1,yp2,yp3,yp4,ym1,ym2,ym3,ym4;
 14:   PetscReal      *values;
 15:   PetscViewer    viewer_in,viewer_outp1,viewer_outp2,viewer_outp3,viewer_outp4;
 16:   PetscViewer    viewer_outm1,viewer_outm2,viewer_outm3,viewer_outm4;
 17:   DM             daf,dac1,dac2,dac3,dac4,daf1,daf2,daf3,daf4;
 18:   Vec            scaling_p1,scaling_p2,scaling_p3,scaling_p4;
 19:   Mat            interp_p1,interp_p2,interp_p3,interp_p4,interp_m1,interp_m2,interp_m3,interp_m4;

 21:   PetscInitialize(&argc,&argv, (char*)0, help);if (ierr) return ierr;
 22:   DMDACreate2d(PETSC_COMM_WORLD,DM_BOUNDARY_PERIODIC,DM_BOUNDARY_PERIODIC,DMDA_STENCIL_BOX,1024,1024,PETSC_DECIDE,PETSC_DECIDE,1,1,NULL,NULL,&daf);
 23:   DMSetFromOptions(daf);
 24:   DMSetUp(daf);
 25:   DMCreateGlobalVector(daf,&x);
 26:   VecGetArray(x,&values);

 28:   DMCoarsen(daf,PETSC_COMM_WORLD,&dac1);
 29:   DMCoarsen(dac1,PETSC_COMM_WORLD,&dac2);
 30:   DMCoarsen(dac2,PETSC_COMM_WORLD,&dac3);
 31:   DMCoarsen(dac3,PETSC_COMM_WORLD,&dac4);
 32:   DMRefine(daf,PETSC_COMM_WORLD,&daf1);
 33:   DMRefine(daf1,PETSC_COMM_WORLD,&daf2);
 34:   DMRefine(daf2,PETSC_COMM_WORLD,&daf3);
 35:   DMRefine(daf3,PETSC_COMM_WORLD,&daf4);

 37:   DMCreateGlobalVector(dac1,&yp1);
 38:   DMCreateGlobalVector(dac2,&yp2);
 39:   DMCreateGlobalVector(dac3,&yp3);
 40:   DMCreateGlobalVector(dac4,&yp4);
 41:   DMCreateGlobalVector(daf1,&ym1);
 42:   DMCreateGlobalVector(daf2,&ym2);
 43:   DMCreateGlobalVector(daf3,&ym3);
 44:   DMCreateGlobalVector(daf4,&ym4);

 46:   DMCreateInterpolation(dac1,daf,&interp_p1,&scaling_p1);
 47:   DMCreateInterpolation(dac2,dac1,&interp_p2,&scaling_p2);
 48:   DMCreateInterpolation(dac3,dac2,&interp_p3,&scaling_p3);
 49:   DMCreateInterpolation(dac4,dac3,&interp_p4,&scaling_p4);
 50:   DMCreateInterpolation(daf,daf1,&interp_m1,NULL);
 51:   DMCreateInterpolation(daf1,daf2,&interp_m2,NULL);
 52:   DMCreateInterpolation(daf2,daf3,&interp_m3,NULL);
 53:   DMCreateInterpolation(daf3,daf4,&interp_m4,NULL);

 55:   PetscViewerBinaryOpen(PETSC_COMM_WORLD,"phi",FILE_MODE_READ,&viewer_in);
 56:   PetscViewerBinaryRead(viewer_in,values,1048576,NULL,PETSC_DOUBLE);
 57:   MatRestrict(interp_p1,x,yp1);
 58:   VecPointwiseMult(yp1,yp1,scaling_p1);
 59:   MatRestrict(interp_p2,yp1,yp2);
 60:   VecPointwiseMult(yp2,yp2,scaling_p2);
 61:   MatRestrict(interp_p3,yp2,yp3);
 62:   VecPointwiseMult(yp3,yp3,scaling_p3);
 63:   MatRestrict(interp_p4,yp3,yp4);
 64:   VecPointwiseMult(yp4,yp4,scaling_p4);
 65:   MatRestrict(interp_m1,x,ym1);
 66:   MatRestrict(interp_m2,ym1,ym2);
 67:   MatRestrict(interp_m3,ym2,ym3);
 68:   MatRestrict(interp_m4,ym3,ym4);

 70:   PetscViewerBinaryOpen(PETSC_COMM_WORLD,"phi1",FILE_MODE_WRITE,&viewer_outp1);
 71:   PetscViewerBinaryOpen(PETSC_COMM_WORLD,"phi2",FILE_MODE_WRITE,&viewer_outp2);
 72:   PetscViewerBinaryOpen(PETSC_COMM_WORLD,"phi3",FILE_MODE_WRITE,&viewer_outp3);
 73:   PetscViewerBinaryOpen(PETSC_COMM_WORLD,"phi4",FILE_MODE_WRITE,&viewer_outp4);
 74:   PetscViewerBinaryOpen(PETSC_COMM_WORLD,"phim1",FILE_MODE_WRITE,&viewer_outm1);
 75:   PetscViewerBinaryOpen(PETSC_COMM_WORLD,"phim2",FILE_MODE_WRITE,&viewer_outm2);
 76:   PetscViewerBinaryOpen(PETSC_COMM_WORLD,"phim3",FILE_MODE_WRITE,&viewer_outm3);
 77:   PetscViewerBinaryOpen(PETSC_COMM_WORLD,"phim4",FILE_MODE_WRITE,&viewer_outm4);

 79:   VecView(yp1,viewer_outp1);
 80:   VecView(x,viewer_outp1);
 81:   VecView(yp2,viewer_outp2);
 82:   VecView(yp3,viewer_outp3);
 83:   VecView(yp4,viewer_outp4);
 84:   VecView(ym1,viewer_outm1);
 85:   VecView(ym2,viewer_outm2);
 86:   VecView(ym3,viewer_outm3);
 87:   VecView(ym4,viewer_outm4);

 89:   PetscViewerDestroy(&viewer_in);
 90:   PetscViewerDestroy(&viewer_outp1);
 91:   PetscViewerDestroy(&viewer_outp2);
 92:   PetscViewerDestroy(&viewer_outp3);
 93:   PetscViewerDestroy(&viewer_outp4);

 95:   PetscViewerDestroy(&viewer_outm1);
 96:   PetscViewerDestroy(&viewer_outm2);
 97:   PetscViewerDestroy(&viewer_outm3);
 98:   PetscViewerDestroy(&viewer_outm4);

100:   VecDestroy(&x);
101:   VecDestroy(&yp1);
102:   VecDestroy(&yp2);
103:   VecDestroy(&yp3);
104:   VecDestroy(&yp4);
105:   VecDestroy(&ym1);
106:   VecDestroy(&ym2);
107:   VecDestroy(&ym3);
108:   VecDestroy(&ym4);
109:   PetscFinalize();
110:   return ierr;
111: }