2: static char help[] = "Tests DAGetInterpolation for nonuniform DA coordinates.\n\n";

 4:  #include petscda.h
 5:  #include petscsys.h

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

 18:   DASetUniformCoordinates(da,0.0,1.0,0.0,1.0,0.0,1.0);
 19:   DAGetCoordinateDA(da,&cda);
 20:   DAGetGhostedCoordinates(da,&gc);
 21:   DAVecGetArray(cda,gc,&coors);
 22:   DAGetCorners(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:   DAVecRestoreArray(cda,gc,&coors);
 29:   DAGetCoordinates(da,&global);
 30:   DALocalToGlobal(cda,gc,INSERT_VALUES,global);
 31:   return(0);
 32: }

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

 45:   DASetUniformCoordinates(da,0.0,1.0,0.0,1.0,0.0,1.0);
 46:   DAGetCoordinateDA(da,&cda);
 47:   DAGetGhostedCoordinates(da,&gc);
 48:   DAVecGetArray(cda,gc,&coors);
 49:   DAGetCorners(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:   DAVecRestoreArray(cda,gc,&coors);
 61:   DAGetCoordinates(da,&global);
 62:   DALocalToGlobal(cda,gc,INSERT_VALUES,global);
 63:   return(0);
 64: }

 68: PetscErrorCode SetCoordinates3d(DA da)
 69: {
 71:   PetscInt       i,j,mstart,m,nstart,n,pstart,p,k;
 72:   Vec            gc,global;
 73:   DACoor3d       ***coors;
 74:   DA             cda;

 77:   DASetUniformCoordinates(da,0.0,1.0,0.0,1.0,0.0,1.0);
 78:   DAGetCoordinateDA(da,&cda);
 79:   DAGetGhostedCoordinates(da,&gc);
 80:   DAVecGetArray(cda,gc,&coors);
 81:   DAGetCorners(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 = coors[k][j][i-1].x + .1*(coors[k][j][i+1].x - coors[k][j][i-1].x);
 87:         }
 88:         if (j % 2) {
 89:           coors[k][j][i].y = coors[k][j-1][i].y + .3*(coors[k][j+1][i].y - coors[k][j-1][i].y);
 90:         }
 91:         if (k % 2) {
 92:           coors[k][j][i].z = coors[k-1][j][i].z + .4*(coors[k+1][j][i].z - coors[k-1][j][i].z);
 93:         }
 94:       }
 95:     }
 96:   }
 97:   DAVecRestoreArray(cda,gc,&coors);
 98:   DAGetCoordinates(da,&global);
 99:   DALocalToGlobal(cda,gc,INSERT_VALUES,global);
100:   return(0);
101: }

105: int main(int argc,char **argv)
106: {
107:   PetscInt       M = 5,N = 4,P = 3, m = PETSC_DECIDE,n = PETSC_DECIDE,p = PETSC_DECIDE,dim = 1;
109:   DA             dac,daf;
110:   DAPeriodicType ptype = DA_NONPERIODIC;
111:   DAStencilType  stype = DA_STENCIL_BOX;
112:   Mat            A;

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

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

125:   /* Create distributed array and get vectors */
126:   if (dim == 1) {
127:     DACreate1d(PETSC_COMM_WORLD,ptype,M,1,1,PETSC_NULL,&dac);
128:   } else if (dim == 2) {
129:     DACreate2d(PETSC_COMM_WORLD,ptype,stype,M,N,PETSC_DECIDE,PETSC_DECIDE,1,1,PETSC_NULL,PETSC_NULL,&dac);
130:   } else if (dim == 3) {
131:     DACreate3d(PETSC_COMM_WORLD,ptype,stype,M,N,P,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,1,1,PETSC_NULL,PETSC_NULL,PETSC_NULL,&dac);
132:   }

134:   DARefine(dac,PETSC_COMM_WORLD,&daf);

136:   DASetUniformCoordinates(dac,0.0,1.0,0.0,1.0,0.0,1.0);
137:   if (dim == 1) {
138:     SetCoordinates1d(daf);
139:   } else if (dim == 2) {
140:     SetCoordinates2d(daf);
141:   } else if (dim == 3) {
142:     SetCoordinates3d(daf);
143:   }
144:   DAGetInterpolation(dac,daf,&A,0);


147:   /* Free memory */
148:   DADestroy(dac);
149:   DADestroy(daf);
150:   MatDestroy(A);
151:   PetscFinalize();
152:   return 0;
153: }
154: