Actual source code: ex3.c


  2: static char help[] = "Tests DMCreateInterpolation() for nonuniform DMDA coordinates.\n\n";

  4: #include <petscdm.h>
  5: #include <petscdmda.h>

  7: PetscErrorCode SetCoordinates1d(DM da)
  8: {
  9:   PetscInt       i,start,m;
 10:   Vec            local,global;
 11:   PetscScalar    *coors,*coorslocal;
 12:   DM             cda;

 15:   DMDASetUniformCoordinates(da,0.0,1.0,0.0,1.0,0.0,1.0);
 16:   DMGetCoordinateDM(da,&cda);
 17:   DMGetCoordinates(da,&global);
 18:   DMGetCoordinatesLocal(da,&local);
 19:   DMDAVecGetArray(cda,global,&coors);
 20:   DMDAVecGetArrayRead(cda,local,&coorslocal);
 21:   DMDAGetCorners(cda,&start,0,0,&m,0,0);
 22:   for (i=start; i<start+m; i++) {
 23:     if (i % 2) {
 24:       coors[i] = coorslocal[i-1] + .1*(coorslocal[i+1] - coorslocal[i-1]);
 25:     }
 26:   }
 27:   DMDAVecRestoreArray(cda,global,&coors);
 28:   DMDAVecRestoreArrayRead(cda,local,&coorslocal);
 29:   DMGlobalToLocalBegin(cda,global,INSERT_VALUES,local);
 30:   DMGlobalToLocalEnd(cda,global,INSERT_VALUES,local);
 31:   return 0;
 32: }

 34: PetscErrorCode SetCoordinates2d(DM da)
 35: {
 36:   PetscInt       i,j,mstart,m,nstart,n;
 37:   Vec            local,global;
 38:   DMDACoor2d     **coors,**coorslocal;
 39:   DM             cda;

 42:   DMDASetUniformCoordinates(da,0.0,1.0,0.0,1.0,0.0,1.0);
 43:   DMGetCoordinateDM(da,&cda);
 44:   DMGetCoordinates(da,&global);
 45:   DMGetCoordinatesLocal(da,&local);
 46:   DMDAVecGetArray(cda,global,&coors);
 47:   DMDAVecGetArrayRead(cda,local,&coorslocal);
 48:   DMDAGetCorners(cda,&mstart,&nstart,0,&m,&n,0);
 49:   for (i=mstart; i<mstart+m; i++) {
 50:     for (j=nstart; j<nstart+n; j++) {
 51:       if (i % 2) {
 52:         coors[j][i].x = coorslocal[j][i-1].x + .1*(coorslocal[j][i+1].x - coorslocal[j][i-1].x);
 53:       }
 54:       if (j % 2) {
 55:         coors[j][i].y = coorslocal[j-1][i].y + .3*(coorslocal[j+1][i].y - coorslocal[j-1][i].y);
 56:       }
 57:     }
 58:   }
 59:   DMDAVecRestoreArray(cda,global,&coors);
 60:   DMDAVecRestoreArrayRead(cda,local,&coorslocal);

 62:   DMGlobalToLocalBegin(cda,global,INSERT_VALUES,local);
 63:   DMGlobalToLocalEnd(cda,global,INSERT_VALUES,local);
 64:   return 0;
 65: }

 67: PetscErrorCode SetCoordinates3d(DM da)
 68: {
 69:   PetscInt       i,j,mstart,m,nstart,n,pstart,p,k;
 70:   Vec            local,global;
 71:   DMDACoor3d     ***coors,***coorslocal;
 72:   DM             cda;

 75:   DMDASetUniformCoordinates(da,0.0,1.0,0.0,1.0,0.0,1.0);
 76:   DMGetCoordinateDM(da,&cda);
 77:   DMGetCoordinates(da,&global);
 78:   DMGetCoordinatesLocal(da,&local);
 79:   DMDAVecGetArray(cda,global,&coors);
 80:   DMDAVecGetArrayRead(cda,local,&coorslocal);
 81:   DMDAGetCorners(cda,&mstart,&nstart,&pstart,&m,&n,&p);
 82:   for (i=mstart; i<mstart+m; i++) {
 83:     for (j=nstart; j<nstart+n; j++) {
 84:       for (k=pstart; k<pstart+p; k++) {
 85:         if (i % 2) {
 86:           coors[k][j][i].x = coorslocal[k][j][i-1].x + .1*(coorslocal[k][j][i+1].x - coorslocal[k][j][i-1].x);
 87:         }
 88:         if (j % 2) {
 89:           coors[k][j][i].y = coorslocal[k][j-1][i].y + .3*(coorslocal[k][j+1][i].y - coorslocal[k][j-1][i].y);
 90:         }
 91:         if (k % 2) {
 92:           coors[k][j][i].z = coorslocal[k-1][j][i].z + .4*(coorslocal[k+1][j][i].z - coorslocal[k-1][j][i].z);
 93:         }
 94:       }
 95:     }
 96:   }
 97:   DMDAVecRestoreArray(cda,global,&coors);
 98:   DMDAVecRestoreArrayRead(cda,local,&coorslocal);
 99:   DMGlobalToLocalBegin(cda,global,INSERT_VALUES,local);
100:   DMGlobalToLocalEnd(cda,global,INSERT_VALUES,local);
101:   return 0;
102: }

104: int main(int argc,char **argv)
105: {
106:   PetscInt         M = 5,N = 4,P = 3, m = PETSC_DECIDE,n = PETSC_DECIDE,p = PETSC_DECIDE,dim = 1;
107:   DM               dac,daf;
108:   DMBoundaryType   bx    = DM_BOUNDARY_NONE,by=DM_BOUNDARY_NONE,bz=DM_BOUNDARY_NONE;
109:   DMDAStencilType  stype = DMDA_STENCIL_BOX;
110:   Mat              A;

112:   PetscInitialize(&argc,&argv,(char*)0,help);
113:   PetscOptionsGetInt(NULL,NULL,"-M",&M,NULL);
114:   PetscOptionsGetInt(NULL,NULL,"-N",&N,NULL);
115:   PetscOptionsGetInt(NULL,NULL,"-P",&P,NULL);
116:   PetscOptionsGetInt(NULL,NULL,"-m",&m,NULL);
117:   PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL);
118:   PetscOptionsGetInt(NULL,NULL,"-p",&p,NULL);
119:   PetscOptionsGetInt(NULL,NULL,"-dim",&dim,NULL);

121:   /* Create distributed array and get vectors */
122:   if (dim == 1) {
123:     DMDACreate1d(PETSC_COMM_WORLD,bx,M,1,1,NULL,&dac);
124:   } else if (dim == 2) {
125:     DMDACreate2d(PETSC_COMM_WORLD,bx,by,stype,M,N,PETSC_DECIDE,PETSC_DECIDE,1,1,NULL,NULL,&dac);
126:   } else if (dim == 3) {
127:     DMDACreate3d(PETSC_COMM_WORLD,bx,by,bz,stype,M,N,P,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,1,1,NULL,NULL,NULL,&dac);
128:   } else SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"dim must be 1,2, or 3");
129:   DMSetFromOptions(dac);
130:   DMSetUp(dac);

132:   DMRefine(dac,PETSC_COMM_WORLD,&daf);

134:   DMDASetUniformCoordinates(dac,0.0,1.0,0.0,1.0,0.0,1.0);
135:   if (dim == 1) {
136:     SetCoordinates1d(daf);
137:   } else if (dim == 2) {
138:     SetCoordinates2d(daf);
139:   } else if (dim == 3) {
140:     SetCoordinates3d(daf);
141:   }
142:   DMCreateInterpolation(dac,daf,&A,0);

144:   /* Free memory */
145:   DMDestroy(&dac);
146:   DMDestroy(&daf);
147:   MatDestroy(&A);
148:   PetscFinalize();
149:   return 0;
150: }

152: /*TEST

154:    test:
155:       nsize: 3
156:       args: -mat_view

158:    test:
159:       suffix: 2
160:       nsize: 3
161:       args: -mat_view -dim 2

163:    test:
164:       suffix: 3
165:       nsize: 3
166:       args: -mat_view -dim 3

168: TEST*/