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

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

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

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

 39: PetscErrorCode SetCoordinates2d(DM da)
 40: {
 42:   PetscInt       i,j,mstart,m,nstart,n;
 43:   Vec            local,global;
 44:   DMDACoor2d     **coors,**coorslocal;
 45:   DM             cda;

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

 68:   DMGlobalToLocalBegin(cda,global,INSERT_VALUES,local);
 69:   DMGlobalToLocalEnd(cda,global,INSERT_VALUES,local);
 70:   return(0);
 71: }

 75: PetscErrorCode SetCoordinates3d(DM da)
 76: {
 78:   PetscInt       i,j,mstart,m,nstart,n,pstart,p,k;
 79:   Vec            local,global;
 80:   DMDACoor3d     ***coors,***coorslocal;
 81:   DM             cda;

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

115: int main(int argc,char **argv)
116: {
117:   PetscInt         M = 5,N = 4,P = 3, m = PETSC_DECIDE,n = PETSC_DECIDE,p = PETSC_DECIDE,dim = 1;
118:   PetscErrorCode   ierr;
119:   DM               dac,daf;
120:   DMBoundaryType   bx    = DM_BOUNDARY_NONE,by=DM_BOUNDARY_NONE,bz=DM_BOUNDARY_NONE;
121:   DMDAStencilType  stype = DMDA_STENCIL_BOX;
122:   Mat              A;

124:   PetscInitialize(&argc,&argv,(char*)0,help);

126:   /* Read options */
127:   PetscOptionsGetInt(NULL,NULL,"-M",&M,NULL);
128:   PetscOptionsGetInt(NULL,NULL,"-N",&N,NULL);
129:   PetscOptionsGetInt(NULL,NULL,"-P",&P,NULL);
130:   PetscOptionsGetInt(NULL,NULL,"-m",&m,NULL);
131:   PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL);
132:   PetscOptionsGetInt(NULL,NULL,"-p",&p,NULL);
133:   PetscOptionsGetInt(NULL,NULL,"-dim",&dim,NULL);

135:   /* Create distributed array and get vectors */
136:   if (dim == 1) {
137:     DMDACreate1d(PETSC_COMM_WORLD,bx,M,1,1,NULL,&dac);
138:   } else if (dim == 2) {
139:     DMDACreate2d(PETSC_COMM_WORLD,bx,by,stype,M,N,PETSC_DECIDE,PETSC_DECIDE,1,1,NULL,NULL,&dac);
140:   } else if (dim == 3) {
141:     DMDACreate3d(PETSC_COMM_WORLD,bx,by,bz,stype,M,N,P,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,1,1,NULL,NULL,NULL,&dac);
142:   } else SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"dim must be 1,2, or 3");

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

146:   DMDASetUniformCoordinates(dac,0.0,1.0,0.0,1.0,0.0,1.0);
147:   if (dim == 1) {
148:     SetCoordinates1d(daf);
149:   } else if (dim == 2) {
150:     SetCoordinates2d(daf);
151:   } else if (dim == 3) {
152:     SetCoordinates3d(daf);
153:   }
154:   DMCreateInterpolation(dac,daf,&A,0);

156:   /* Free memory */
157:   DMDestroy(&dac);
158:   DMDestroy(&daf);
159:   MatDestroy(&A);
160:   PetscFinalize();
161:   return 0;
162: }