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



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


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

 25:   PetscInitialize(&argc,&argv, (char*)0, help);
 26: #if defined(PETSC_USE_COMPLEX)
 27:   SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"Not for complex numbers");
 28: #else
 29:   DMDACreate2d(PETSC_COMM_WORLD,DM_BOUNDARY_PERIODIC,DM_BOUNDARY_PERIODIC,DMDA_STENCIL_BOX,1024,1024,PETSC_DECIDE,PETSC_DECIDE, 1, 1,NULL,NULL,&daf);
 30:   DMCreateGlobalVector(daf,&x);
 31:   VecGetArray(x,&values);

 33:   DMCoarsen(daf,PETSC_COMM_WORLD,&dac1);
 34:   DMCoarsen(dac1,PETSC_COMM_WORLD,&dac2);
 35:   DMCoarsen(dac2,PETSC_COMM_WORLD,&dac3);
 36:   DMCoarsen(dac3,PETSC_COMM_WORLD,&dac4);
 37:   DMRefine(daf,PETSC_COMM_WORLD,&daf1);
 38:   DMRefine(daf1,PETSC_COMM_WORLD,&daf2);
 39:   DMRefine(daf2,PETSC_COMM_WORLD,&daf3);
 40:   DMRefine(daf3,PETSC_COMM_WORLD,&daf4);

 42:   DMCreateGlobalVector(dac1,&yp1);
 43:   DMCreateGlobalVector(dac2,&yp2);
 44:   DMCreateGlobalVector(dac3,&yp3);
 45:   DMCreateGlobalVector(dac4,&yp4);
 46:   DMCreateGlobalVector(daf1,&ym1);
 47:   DMCreateGlobalVector(daf2,&ym2);
 48:   DMCreateGlobalVector(daf3,&ym3);
 49:   DMCreateGlobalVector(daf4,&ym4);

 51:   DMCreateInterpolation(dac1,daf,&interp_p1,&scaling_p1);
 52:   DMCreateInterpolation(dac2,dac1,&interp_p2,&scaling_p2);
 53:   DMCreateInterpolation(dac3,dac2,&interp_p3,&scaling_p3);
 54:   DMCreateInterpolation(dac4,dac3,&interp_p4,&scaling_p4);
 55:   DMCreateInterpolation(daf,daf1,&interp_m1,NULL);
 56:   DMCreateInterpolation(daf1,daf2,&interp_m2,NULL);
 57:   DMCreateInterpolation(daf2,daf3,&interp_m3,NULL);
 58:   DMCreateInterpolation(daf3,daf4,&interp_m4,NULL);

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

 75:   PetscViewerBinaryOpen(PETSC_COMM_WORLD,"phi1",FILE_MODE_WRITE,&viewer_outp1);
 76:   PetscViewerBinaryOpen(PETSC_COMM_WORLD,"phi2",FILE_MODE_WRITE,&viewer_outp2);
 77:   PetscViewerBinaryOpen(PETSC_COMM_WORLD,"phi3",FILE_MODE_WRITE,&viewer_outp3);
 78:   PetscViewerBinaryOpen(PETSC_COMM_WORLD,"phi4",FILE_MODE_WRITE,&viewer_outp4);
 79:   PetscViewerBinaryOpen(PETSC_COMM_WORLD,"phim1",FILE_MODE_WRITE,&viewer_outm1);
 80:   PetscViewerBinaryOpen(PETSC_COMM_WORLD,"phim2",FILE_MODE_WRITE,&viewer_outm2);
 81:   PetscViewerBinaryOpen(PETSC_COMM_WORLD,"phim3",FILE_MODE_WRITE,&viewer_outm3);
 82:   PetscViewerBinaryOpen(PETSC_COMM_WORLD,"phim4",FILE_MODE_WRITE,&viewer_outm4);

 84:   VecView(yp1,viewer_outp1);
 85:   VecView(x,viewer_outp1);
 86:   VecView(yp2,viewer_outp2);
 87:   VecView(yp3,viewer_outp3);
 88:   VecView(yp4,viewer_outp4);
 89:   VecView(ym1,viewer_outm1);
 90:   VecView(ym2,viewer_outm2);
 91:   VecView(ym3,viewer_outm3);
 92:   VecView(ym4,viewer_outm4);

 94:   PetscViewerDestroy(&viewer_in);
 95:   PetscViewerDestroy(&viewer_outp1);
 96:   PetscViewerDestroy(&viewer_outp2);
 97:   PetscViewerDestroy(&viewer_outp3);
 98:   PetscViewerDestroy(&viewer_outp4);

100:   PetscViewerDestroy(&viewer_outm1);
101:   PetscViewerDestroy(&viewer_outm2);
102:   PetscViewerDestroy(&viewer_outm3);
103:   PetscViewerDestroy(&viewer_outm4);

105:   VecDestroy(&x);
106:   VecDestroy(&yp1);
107:   VecDestroy(&yp2);
108:   VecDestroy(&yp3);
109:   VecDestroy(&yp4);
110:   VecDestroy(&ym1);
111:   VecDestroy(&ym2);
112:   VecDestroy(&ym3);
113:   VecDestroy(&ym4);
114: #endif
115:   PetscFinalize();
116:   return 0;
117: }