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            gc,global;
 14:   PetscScalar    *coors;
 15:   DM             cda;

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

 37: PetscErrorCode SetCoordinates2d(DM da)
 38: {
 40:   PetscInt       i,j,mstart,m,nstart,n;
 41:   Vec            gc,global;
 42:   DMDACoor2d     **coors;
 43:   DM             cda;

 46:   DMDASetUniformCoordinates(da,0.0,1.0,0.0,1.0,0.0,1.0);
 47:   DMGetCoordinateDM(da,&cda);
 48:   DMGetCoordinatesLocal(da,&gc);
 49:   DMDAVecGetArray(cda,gc,&coors);
 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 = coors[j][i-1].x + .1*(coors[j][i+1].x - coors[j][i-1].x);
 55:       }
 56:       if (j % 2) {
 57:         coors[j][i].y = coors[j-1][i].y + .3*(coors[j+1][i].y - coors[j-1][i].y);
 58:       }
 59:     }
 60:   }
 61:   DMDAVecRestoreArray(cda,gc,&coors);
 62:   DMGetCoordinates(da,&global);
 63:   DMLocalToGlobalBegin(cda,gc,INSERT_VALUES,global);
 64:   DMLocalToGlobalEnd(cda,gc,INSERT_VALUES,global);
 65:   return(0);
 66: }

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

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

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

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

119:   /* Read options */
120:   PetscOptionsGetInt(NULL,"-M",&M,NULL);
121:   PetscOptionsGetInt(NULL,"-N",&N,NULL);
122:   PetscOptionsGetInt(NULL,"-P",&P,NULL);
123:   PetscOptionsGetInt(NULL,"-m",&m,NULL);
124:   PetscOptionsGetInt(NULL,"-n",&n,NULL);
125:   PetscOptionsGetInt(NULL,"-p",&p,NULL);
126:   PetscOptionsGetInt(NULL,"-dim",&dim,NULL);

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

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

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

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