Actual source code: ex3.c
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: VecDestroy(gc);
30: DAGetCoordinates(da,&global);
31: DALocalToGlobal(cda,gc,INSERT_VALUES,global);
32: VecDestroy(global);
33: DADestroy(cda);
34: return(0);
35: }
39: PetscErrorCode SetCoordinates2d(DA da)
40: {
42: PetscInt i,j,mstart,m,nstart,n;
43: Vec gc,global;
44: DACoor2d **coors;
45: DA cda;
48: DASetUniformCoordinates(da,0.0,1.0,0.0,1.0,0.0,1.0);
49: DAGetCoordinateDA(da,&cda);
50: DAGetGhostedCoordinates(da,&gc);
51: DAVecGetArray(cda,gc,&coors);
52: DAGetCorners(cda,&mstart,&nstart,0,&m,&n,0);
53: for (i=mstart; i<mstart+m; i++) {
54: for (j=nstart; j<nstart+n; j++) {
55: if (i % 2) {
56: coors[j][i].x = coors[j][i-1].x + .1*(coors[j][i+1].x - coors[j][i-1].x);
57: }
58: if (j % 2) {
59: coors[j][i].y = coors[j-1][i].y + .3*(coors[j+1][i].y - coors[j-1][i].y);
60: }
61: }
62: }
63: DAVecRestoreArray(cda,gc,&coors);
64: VecDestroy(gc);
65: DAGetCoordinates(da,&global);
66: DALocalToGlobal(cda,gc,INSERT_VALUES,global);
67: VecDestroy(global);
68: DADestroy(cda);
69: return(0);
70: }
74: PetscErrorCode SetCoordinates3d(DA da)
75: {
77: PetscInt i,j,mstart,m,nstart,n,pstart,p,k;
78: Vec gc,global;
79: DACoor3d ***coors;
80: DA cda;
83: DASetUniformCoordinates(da,0.0,1.0,0.0,1.0,0.0,1.0);
84: DAGetCoordinateDA(da,&cda);
85: DAGetGhostedCoordinates(da,&gc);
86: DAVecGetArray(cda,gc,&coors);
87: DAGetCorners(cda,&mstart,&nstart,&pstart,&m,&n,&p);
88: for (i=mstart; i<mstart+m; i++) {
89: for (j=nstart; j<nstart+n; j++) {
90: for (k=pstart; k<pstart+p; k++) {
91: if (i % 2) {
92: coors[k][j][i].x = coors[k][j][i-1].x + .1*(coors[k][j][i+1].x - coors[k][j][i-1].x);
93: }
94: if (j % 2) {
95: coors[k][j][i].y = coors[k][j-1][i].y + .3*(coors[k][j+1][i].y - coors[k][j-1][i].y);
96: }
97: if (k % 2) {
98: coors[k][j][i].z = coors[k-1][j][i].z + .4*(coors[k+1][j][i].z - coors[k-1][j][i].z);
99: }
100: }
101: }
102: }
103: DAVecRestoreArray(cda,gc,&coors);
104: VecDestroy(gc);
105: DAGetCoordinates(da,&global);
106: DALocalToGlobal(cda,gc,INSERT_VALUES,global);
107: VecDestroy(global);
108: DADestroy(cda);
109: return(0);
110: }
114: int main(int argc,char **argv)
115: {
116: PetscInt M = 5,N = 4,P = 3, m = PETSC_DECIDE,n = PETSC_DECIDE,p = PETSC_DECIDE,dim = 1;
118: DA dac,daf;
119: DAPeriodicType ptype = DA_NONPERIODIC;
120: DAStencilType stype = DA_STENCIL_BOX;
121: Mat A;
123: PetscInitialize(&argc,&argv,(char*)0,help);
125: /* Read options */
126: PetscOptionsGetInt(PETSC_NULL,"-M",&M,PETSC_NULL);
127: PetscOptionsGetInt(PETSC_NULL,"-N",&N,PETSC_NULL);
128: PetscOptionsGetInt(PETSC_NULL,"-P",&P,PETSC_NULL);
129: PetscOptionsGetInt(PETSC_NULL,"-m",&m,PETSC_NULL);
130: PetscOptionsGetInt(PETSC_NULL,"-n",&n,PETSC_NULL);
131: PetscOptionsGetInt(PETSC_NULL,"-p",&p,PETSC_NULL);
132: PetscOptionsGetInt(PETSC_NULL,"-dim",&dim,PETSC_NULL);
134: /* Create distributed array and get vectors */
135: if (dim == 1) {
136: DACreate1d(PETSC_COMM_WORLD,ptype,M,1,1,PETSC_NULL,&dac);
137: } else if (dim == 2) {
138: DACreate2d(PETSC_COMM_WORLD,ptype,stype,M,N,PETSC_DECIDE,PETSC_DECIDE,1,1,PETSC_NULL,PETSC_NULL,&dac);
139: } else if (dim == 3) {
140: DACreate3d(PETSC_COMM_WORLD,ptype,stype,M,N,P,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,1,1,PETSC_NULL,PETSC_NULL,PETSC_NULL,&dac);
141: }
143: DARefine(dac,PETSC_COMM_WORLD,&daf);
145: DASetUniformCoordinates(dac,0.0,1.0,0.0,1.0,0.0,1.0);
146: if (dim == 1) {
147: SetCoordinates1d(daf);
148: } else if (dim == 2) {
149: SetCoordinates2d(daf);
150: } else if (dim == 3) {
151: SetCoordinates3d(daf);
152: }
153: DAGetInterpolation(dac,daf,&A,0);
156: /* Free memory */
157: DADestroy(dac);
158: DADestroy(daf);
159: MatDestroy(A);
160: PetscFinalize();
161: return 0;
162: }
163: