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

  4: #include <petscdmda.h>

  8: PetscErrorCode SetCoordinates1d(DM da)
  9: {
 11:   PetscInt       i,start,m;
 12:   Vec            gc,global;
 13:   PetscScalar    *coors;
 14:   DM             cda;

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

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

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

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

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

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:   DMDABoundaryType bx    = DMDA_BOUNDARY_NONE,by=DMDA_BOUNDARY_NONE,bz=DMDA_BOUNDARY_NONE;
113:   DMDAStencilType  stype = DMDA_STENCIL_BOX;
114:   Mat              A;

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

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

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

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);
147:   MatViewFromOptions(A,NULL,"-mat_view");


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