Actual source code: ex3.c
2: static char help[] = "Tests DMCreateInterpolation() for nonuniform DMDA coordinates.\n\n";
4: #include <petscdm.h>
5: #include <petscdmda.h>
7: PetscErrorCode SetCoordinates1d(DM da)
8: {
9: PetscInt i,start,m;
10: Vec local,global;
11: PetscScalar *coors,*coorslocal;
12: DM cda;
15: DMDASetUniformCoordinates(da,0.0,1.0,0.0,1.0,0.0,1.0);
16: DMGetCoordinateDM(da,&cda);
17: DMGetCoordinates(da,&global);
18: DMGetCoordinatesLocal(da,&local);
19: DMDAVecGetArray(cda,global,&coors);
20: DMDAVecGetArrayRead(cda,local,&coorslocal);
21: DMDAGetCorners(cda,&start,0,0,&m,0,0);
22: for (i=start; i<start+m; i++) {
23: if (i % 2) {
24: coors[i] = coorslocal[i-1] + .1*(coorslocal[i+1] - coorslocal[i-1]);
25: }
26: }
27: DMDAVecRestoreArray(cda,global,&coors);
28: DMDAVecRestoreArrayRead(cda,local,&coorslocal);
29: DMGlobalToLocalBegin(cda,global,INSERT_VALUES,local);
30: DMGlobalToLocalEnd(cda,global,INSERT_VALUES,local);
31: return 0;
32: }
34: PetscErrorCode SetCoordinates2d(DM da)
35: {
36: PetscInt i,j,mstart,m,nstart,n;
37: Vec local,global;
38: DMDACoor2d **coors,**coorslocal;
39: DM cda;
42: DMDASetUniformCoordinates(da,0.0,1.0,0.0,1.0,0.0,1.0);
43: DMGetCoordinateDM(da,&cda);
44: DMGetCoordinates(da,&global);
45: DMGetCoordinatesLocal(da,&local);
46: DMDAVecGetArray(cda,global,&coors);
47: DMDAVecGetArrayRead(cda,local,&coorslocal);
48: DMDAGetCorners(cda,&mstart,&nstart,0,&m,&n,0);
49: for (i=mstart; i<mstart+m; i++) {
50: for (j=nstart; j<nstart+n; j++) {
51: if (i % 2) {
52: coors[j][i].x = coorslocal[j][i-1].x + .1*(coorslocal[j][i+1].x - coorslocal[j][i-1].x);
53: }
54: if (j % 2) {
55: coors[j][i].y = coorslocal[j-1][i].y + .3*(coorslocal[j+1][i].y - coorslocal[j-1][i].y);
56: }
57: }
58: }
59: DMDAVecRestoreArray(cda,global,&coors);
60: DMDAVecRestoreArrayRead(cda,local,&coorslocal);
62: DMGlobalToLocalBegin(cda,global,INSERT_VALUES,local);
63: DMGlobalToLocalEnd(cda,global,INSERT_VALUES,local);
64: return 0;
65: }
67: PetscErrorCode SetCoordinates3d(DM da)
68: {
69: PetscInt i,j,mstart,m,nstart,n,pstart,p,k;
70: Vec local,global;
71: DMDACoor3d ***coors,***coorslocal;
72: DM cda;
75: DMDASetUniformCoordinates(da,0.0,1.0,0.0,1.0,0.0,1.0);
76: DMGetCoordinateDM(da,&cda);
77: DMGetCoordinates(da,&global);
78: DMGetCoordinatesLocal(da,&local);
79: DMDAVecGetArray(cda,global,&coors);
80: DMDAVecGetArrayRead(cda,local,&coorslocal);
81: DMDAGetCorners(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 = coorslocal[k][j][i-1].x + .1*(coorslocal[k][j][i+1].x - coorslocal[k][j][i-1].x);
87: }
88: if (j % 2) {
89: coors[k][j][i].y = coorslocal[k][j-1][i].y + .3*(coorslocal[k][j+1][i].y - coorslocal[k][j-1][i].y);
90: }
91: if (k % 2) {
92: coors[k][j][i].z = coorslocal[k-1][j][i].z + .4*(coorslocal[k+1][j][i].z - coorslocal[k-1][j][i].z);
93: }
94: }
95: }
96: }
97: DMDAVecRestoreArray(cda,global,&coors);
98: DMDAVecRestoreArrayRead(cda,local,&coorslocal);
99: DMGlobalToLocalBegin(cda,global,INSERT_VALUES,local);
100: DMGlobalToLocalEnd(cda,global,INSERT_VALUES,local);
101: return 0;
102: }
104: int main(int argc,char **argv)
105: {
106: PetscInt M = 5,N = 4,P = 3, m = PETSC_DECIDE,n = PETSC_DECIDE,p = PETSC_DECIDE,dim = 1;
107: DM dac,daf;
108: DMBoundaryType bx = DM_BOUNDARY_NONE,by=DM_BOUNDARY_NONE,bz=DM_BOUNDARY_NONE;
109: DMDAStencilType stype = DMDA_STENCIL_BOX;
110: Mat A;
112: PetscInitialize(&argc,&argv,(char*)0,help);
113: PetscOptionsGetInt(NULL,NULL,"-M",&M,NULL);
114: PetscOptionsGetInt(NULL,NULL,"-N",&N,NULL);
115: PetscOptionsGetInt(NULL,NULL,"-P",&P,NULL);
116: PetscOptionsGetInt(NULL,NULL,"-m",&m,NULL);
117: PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL);
118: PetscOptionsGetInt(NULL,NULL,"-p",&p,NULL);
119: PetscOptionsGetInt(NULL,NULL,"-dim",&dim,NULL);
121: /* Create distributed array and get vectors */
122: if (dim == 1) {
123: DMDACreate1d(PETSC_COMM_WORLD,bx,M,1,1,NULL,&dac);
124: } else if (dim == 2) {
125: DMDACreate2d(PETSC_COMM_WORLD,bx,by,stype,M,N,PETSC_DECIDE,PETSC_DECIDE,1,1,NULL,NULL,&dac);
126: } else if (dim == 3) {
127: DMDACreate3d(PETSC_COMM_WORLD,bx,by,bz,stype,M,N,P,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,1,1,NULL,NULL,NULL,&dac);
128: } else SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"dim must be 1,2, or 3");
129: DMSetFromOptions(dac);
130: DMSetUp(dac);
132: DMRefine(dac,PETSC_COMM_WORLD,&daf);
134: DMDASetUniformCoordinates(dac,0.0,1.0,0.0,1.0,0.0,1.0);
135: if (dim == 1) {
136: SetCoordinates1d(daf);
137: } else if (dim == 2) {
138: SetCoordinates2d(daf);
139: } else if (dim == 3) {
140: SetCoordinates3d(daf);
141: }
142: DMCreateInterpolation(dac,daf,&A,0);
144: /* Free memory */
145: DMDestroy(&dac);
146: DMDestroy(&daf);
147: MatDestroy(&A);
148: PetscFinalize();
149: return 0;
150: }
152: /*TEST
154: test:
155: nsize: 3
156: args: -mat_view
158: test:
159: suffix: 2
160: nsize: 3
161: args: -mat_view -dim 2
163: test:
164: suffix: 3
165: nsize: 3
166: args: -mat_view -dim 3
168: TEST*/