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: