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: {
 10:   PetscInt       i,start,m;
 11:   Vec            local,global;
 12:   PetscScalar    *coors,*coorslocal;
 13:   DM             cda;

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

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

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

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

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

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

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

116:   PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
117:   PetscOptionsGetInt(NULL,NULL,"-M",&M,NULL);
118:   PetscOptionsGetInt(NULL,NULL,"-N",&N,NULL);
119:   PetscOptionsGetInt(NULL,NULL,"-P",&P,NULL);
120:   PetscOptionsGetInt(NULL,NULL,"-m",&m,NULL);
121:   PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL);
122:   PetscOptionsGetInt(NULL,NULL,"-p",&p,NULL);
123:   PetscOptionsGetInt(NULL,NULL,"-dim",&dim,NULL);

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

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

138:   DMDASetUniformCoordinates(dac,0.0,1.0,0.0,1.0,0.0,1.0);
139:   if (dim == 1) {
140:     SetCoordinates1d(daf);
141:   } else if (dim == 2) {
142:     SetCoordinates2d(daf);
143:   } else if (dim == 3) {
144:     SetCoordinates3d(daf);
145:   }
146:   DMCreateInterpolation(dac,daf,&A,0);

148:   /* Free memory */
149:   DMDestroy(&dac);
150:   DMDestroy(&daf);
151:   MatDestroy(&A);
152:   PetscFinalize();
153:   return ierr;
154: }