Actual source code: ex21.c

petsc-3.9.4 2018-09-11
Report Typos and Errors
  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>

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

 18:   bx = DM_BOUNDARY_NONE;
 19:   by = DM_BOUNDARY_NONE;
 20:   bz = DM_BOUNDARY_NONE;

 22:   periodicity = 0;

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

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

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

 40:   DMDASetUniformCoordinates(dac, -1.0,1.0, -1.0,1.0, -1.0,1.0);
 41:   DMDASetUniformCoordinates(daf, -1.0,1.0, -1.0,1.0, -1.0,1.0);

 43:   {
 44:     DM         cdaf,cdac;
 45:     Vec        coordsc,coordsf,coordsf2;
 46:     Mat        inject;
 47:     VecScatter vscat;
 48:     Mat        interp;
 49:     PetscReal  norm;

 51:     DMGetCoordinateDM(dac,&cdac);
 52:     DMGetCoordinateDM(daf,&cdaf);

 54:     DMGetCoordinates(dac,&coordsc);
 55:     DMGetCoordinates(daf,&coordsf);

 57:     DMCreateInjection(cdac,cdaf,&inject);
 58:     MatScatterGetVecScatter(inject,&vscat);
 59:     VecScatterBegin(vscat,coordsf,coordsc,INSERT_VALUES,SCATTER_FORWARD);
 60:     VecScatterEnd(vscat  ,coordsf,coordsc,INSERT_VALUES,SCATTER_FORWARD);
 61:     MatDestroy(&inject);

 63:     DMCreateInterpolation(cdac,cdaf,&interp,NULL);
 64:     VecDuplicate(coordsf,&coordsf2);
 65:     MatInterpolate(interp,coordsc,coordsf2);
 66:     VecAXPY(coordsf2,-1.0,coordsf);
 67:     VecNorm(coordsf2,NORM_MAX,&norm);
 68:     /* The fine coordinates are only reproduced in certain cases */
 69:     if (!bx && !by && !bz && norm > PETSC_SQRT_MACHINE_EPSILON) {PetscPrintf(PETSC_COMM_WORLD,"Norm %g\n",(double)norm);}
 70:     VecDestroy(&coordsf2);
 71:     MatDestroy(&interp);
 72:   }

 74:   if (0) {
 75:     DMCreateGlobalVector(dac,&ac);
 76:     VecZeroEntries(ac);

 78:     DMCreateGlobalVector(daf,&af);
 79:     VecZeroEntries(af);

 81:     PetscViewerASCIIOpen(PETSC_COMM_WORLD, "dac_7.vtk", &vv);
 82:     PetscViewerPushFormat(vv, PETSC_VIEWER_ASCII_VTK);
 83:     DMView(dac, vv);
 84:     VecView(ac, vv);
 85:     PetscViewerPopFormat(vv);
 86:     PetscViewerDestroy(&vv);

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

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

107:   PetscInitialize(&argc,&argv,0,help);if (ierr) return ierr;
108:   mx   = 2;
109:   my   = 2;
110:   mz   = 2;
111:   PetscOptionsGetInt(NULL,NULL,"-mx", &mx, 0);
112:   PetscOptionsGetInt(NULL,NULL,"-my", &my, 0);
113:   PetscOptionsGetInt(NULL,NULL,"-mz", &mz, 0);
114:   test1_DAInjection3d(mx,my,mz);
115:   PetscFinalize();
116:   return ierr;
117: }


120: /*TEST

122:       test:
123:          nsize: 10
124:          args: -mx 30 -my 30 -mz 30 -periodic 0 -da_processors_x 10

126:       test:
127:          suffix: 2
128:          nsize: 10
129:          args: -mx 29 -my 30 -mz 30 -periodic 1 -da_processors_x 10

131:       test:
132:          suffix: 3
133:          nsize: 10
134:          args: -mx 30 -my 29 -mz 30 -periodic 2 -da_processors_x 10

136:       test:
137:          suffix: 4
138:          nsize: 10
139:          args: -mx 30 -my 30 -mz 29 -periodic 3 -da_processors_x 10

141:       test:
142:          suffix: 5
143:          nsize: 10
144:          args: -mx 30 -my 30 -mz 30 -periodic 0 -da_processors_y 10

146:       test:
147:          suffix: 6
148:          nsize: 10
149:          args: -mx 29 -my 30 -mz 30 -periodic 1 -da_processors_y 10

151:       test:
152:          suffix: 7
153:          nsize: 10
154:          args: -mx 30 -my 29 -mz 30 -periodic 2 -da_processors_y 10

156:       test:
157:          suffix: 8
158:          nsize: 10
159:          args: -mx 30 -my 30 -mz 29 -periodic 3 -da_processors_y 10

161:       test:
162:          suffix: 9
163:          nsize: 10
164:          args: -mx 30 -my 30 -mz 30 -periodic 0 -da_processors_z 10

166:       test:
167:          suffix: 10
168:          nsize: 10
169:          args: -mx 29 -my 30 -mz 30 -periodic 1 -da_processors_z 10

171:       test:
172:          suffix: 11
173:          nsize: 10
174:          args: -mx 30 -my 29 -mz 30 -periodic 2 -da_processors_z 10

176:       test:
177:          suffix: 12
178:          nsize: 10
179:          args: -mx 30 -my 30 -mz 29 -periodic 3 -da_processors_z 10

181:       test:
182:          suffix: 13
183:          nsize: 10
184:          args: -mx 30 -my 30 -mz 30 -periodic 0

186:       test:
187:          suffix: 14
188:          nsize: 10
189:          args: -mx 29 -my 30 -mz 30 -periodic 1

191:       test:
192:          suffix: 15
193:          nsize: 10
194:          args: -mx 30 -my 29 -mz 30 -periodic 2

196:       test:
197:          suffix: 16
198:          nsize: 10
199:          args: -mx 30 -my 30 -mz 29 -periodic 3

201: TEST*/