Actual source code: ex3.c

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

 4:  #include petscda.h

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

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

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

 47:   DASetUniformCoordinates(da,0.0,1.0,0.0,1.0,0.0,1.0);
 48:   DAGetCoordinateDA(da,&cda);
 49:   DAGetGhostedCoordinates(da,&gc);
 50:   DAVecGetArray(cda,gc,&coors);
 51:   DAGetCorners(cda,&mstart,&nstart,0,&m,&n,0);
 52:   for (i=mstart; i<mstart+m; i++) {
 53:     for (j=nstart; j<nstart+n; j++) {
 54:       if (i % 2) {
 55:         coors[j][i].x = coors[j][i-1].x + .1*(coors[j][i+1].x - coors[j][i-1].x);
 56:       }
 57:       if (j % 2) {
 58:         coors[j][i].y = coors[j-1][i].y + .3*(coors[j+1][i].y - coors[j-1][i].y);
 59:       }
 60:     }
 61:   }
 62:   DAVecRestoreArray(cda,gc,&coors);
 63:   DAGetCoordinates(da,&global);
 64:   DALocalToGlobal(cda,gc,INSERT_VALUES,global);
 65:   VecDestroy(gc);
 66:   VecDestroy(global);
 67:   DADestroy(cda);
 68:   return(0);
 69: }

 73: PetscErrorCode SetCoordinates3d(DA da)
 74: {
 76:   PetscInt       i,j,mstart,m,nstart,n,pstart,p,k;
 77:   Vec            gc,global;
 78:   DACoor3d       ***coors;
 79:   DA             cda;

 82:   DASetUniformCoordinates(da,0.0,1.0,0.0,1.0,0.0,1.0);
 83:   DAGetCoordinateDA(da,&cda);
 84:   DAGetGhostedCoordinates(da,&gc);
 85:   DAVecGetArray(cda,gc,&coors);
 86:   DAGetCorners(cda,&mstart,&nstart,&pstart,&m,&n,&p);
 87:   for (i=mstart; i<mstart+m; i++) {
 88:     for (j=nstart; j<nstart+n; j++) {
 89:       for (k=pstart; k<pstart+p; k++) {
 90:         if (i % 2) {
 91:           coors[k][j][i].x = coors[k][j][i-1].x + .1*(coors[k][j][i+1].x - coors[k][j][i-1].x);
 92:         }
 93:         if (j % 2) {
 94:           coors[k][j][i].y = coors[k][j-1][i].y + .3*(coors[k][j+1][i].y - coors[k][j-1][i].y);
 95:         }
 96:         if (k % 2) {
 97:           coors[k][j][i].z = coors[k-1][j][i].z + .4*(coors[k+1][j][i].z - coors[k-1][j][i].z);
 98:         }
 99:       }
100:     }
101:   }
102:   DAVecRestoreArray(cda,gc,&coors);
103:   DAGetCoordinates(da,&global);
104:   DALocalToGlobal(cda,gc,INSERT_VALUES,global);
105:   VecDestroy(gc);
106:   VecDestroy(global);
107:   DADestroy(cda);
108:   return(0);
109: }

113: int main(int argc,char **argv)
114: {
115:   PetscInt       M = 5,N = 4,P = 3, m = PETSC_DECIDE,n = PETSC_DECIDE,p = PETSC_DECIDE,dim = 1;
117:   DA             dac,daf;
118:   DAPeriodicType ptype = DA_NONPERIODIC;
119:   DAStencilType  stype = DA_STENCIL_BOX;
120:   Mat            A;

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

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

133:   /* Create distributed array and get vectors */
134:   if (dim == 1) {
135:     DACreate1d(PETSC_COMM_WORLD,ptype,M,1,1,PETSC_NULL,&dac);
136:   } else if (dim == 2) {
137:     DACreate2d(PETSC_COMM_WORLD,ptype,stype,M,N,PETSC_DECIDE,PETSC_DECIDE,1,1,PETSC_NULL,PETSC_NULL,&dac);
138:   } else if (dim == 3) {
139:     DACreate3d(PETSC_COMM_WORLD,ptype,stype,M,N,P,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,1,1,PETSC_NULL,PETSC_NULL,PETSC_NULL,&dac);
140:   }

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

144:   DASetUniformCoordinates(dac,0.0,1.0,0.0,1.0,0.0,1.0);
145:   if (dim == 1) {
146:     SetCoordinates1d(daf);
147:   } else if (dim == 2) {
148:     SetCoordinates2d(daf);
149:   } else if (dim == 3) {
150:     SetCoordinates3d(daf);
151:   }
152:   DAGetInterpolation(dac,daf,&A,0);


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