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: {
10: PetscInt i,start,m;
11: Vec local,global;
12: PetscScalar *coors,*coorslocal;
13: DM cda;
16: DMDASetUniformCoordinates(da,0.0,1.0,0.0,1.0,0.0,1.0);
17: DMGetCoordinateDM(da,&cda);
18: DMGetCoordinates(da,&global);
19: DMGetCoordinatesLocal(da,&local);
20: DMDAVecGetArray(cda,global,&coors);
21: DMDAVecGetArrayRead(cda,local,&coorslocal);
22: DMDAGetCorners(cda,&start,0,0,&m,0,0);
23: for (i=start; i<start+m; i++) {
24: if (i % 2) {
25: coors[i] = coorslocal[i-1] + .1*(coorslocal[i+1] - coorslocal[i-1]);
26: }
27: }
28: DMDAVecRestoreArray(cda,global,&coors);
29: DMDAVecRestoreArrayRead(cda,local,&coorslocal);
30: DMGlobalToLocalBegin(cda,global,INSERT_VALUES,local);
31: DMGlobalToLocalEnd(cda,global,INSERT_VALUES,local);
32: return(0);
33: }
35: PetscErrorCode SetCoordinates2d(DM da)
36: {
38: PetscInt i,j,mstart,m,nstart,n;
39: Vec local,global;
40: DMDACoor2d **coors,**coorslocal;
41: DM cda;
44: DMDASetUniformCoordinates(da,0.0,1.0,0.0,1.0,0.0,1.0);
45: DMGetCoordinateDM(da,&cda);
46: DMGetCoordinates(da,&global);
47: DMGetCoordinatesLocal(da,&local);
48: DMDAVecGetArray(cda,global,&coors);
49: DMDAVecGetArrayRead(cda,local,&coorslocal);
50: DMDAGetCorners(cda,&mstart,&nstart,0,&m,&n,0);
51: for (i=mstart; i<mstart+m; i++) {
52: for (j=nstart; j<nstart+n; j++) {
53: if (i % 2) {
54: coors[j][i].x = coorslocal[j][i-1].x + .1*(coorslocal[j][i+1].x - coorslocal[j][i-1].x);
55: }
56: if (j % 2) {
57: coors[j][i].y = coorslocal[j-1][i].y + .3*(coorslocal[j+1][i].y - coorslocal[j-1][i].y);
58: }
59: }
60: }
61: DMDAVecRestoreArray(cda,global,&coors);
62: DMDAVecRestoreArrayRead(cda,local,&coorslocal);
64: DMGlobalToLocalBegin(cda,global,INSERT_VALUES,local);
65: DMGlobalToLocalEnd(cda,global,INSERT_VALUES,local);
66: return(0);
67: }
69: PetscErrorCode SetCoordinates3d(DM da)
70: {
72: PetscInt i,j,mstart,m,nstart,n,pstart,p,k;
73: Vec local,global;
74: DMDACoor3d ***coors,***coorslocal;
75: DM cda;
78: DMDASetUniformCoordinates(da,0.0,1.0,0.0,1.0,0.0,1.0);
79: DMGetCoordinateDM(da,&cda);
80: DMGetCoordinates(da,&global);
81: DMGetCoordinatesLocal(da,&local);
82: DMDAVecGetArray(cda,global,&coors);
83: DMDAVecGetArrayRead(cda,local,&coorslocal);
84: DMDAGetCorners(cda,&mstart,&nstart,&pstart,&m,&n,&p);
85: for (i=mstart; i<mstart+m; i++) {
86: for (j=nstart; j<nstart+n; j++) {
87: for (k=pstart; k<pstart+p; k++) {
88: if (i % 2) {
89: coors[k][j][i].x = coorslocal[k][j][i-1].x + .1*(coorslocal[k][j][i+1].x - coorslocal[k][j][i-1].x);
90: }
91: if (j % 2) {
92: coors[k][j][i].y = coorslocal[k][j-1][i].y + .3*(coorslocal[k][j+1][i].y - coorslocal[k][j-1][i].y);
93: }
94: if (k % 2) {
95: coors[k][j][i].z = coorslocal[k-1][j][i].z + .4*(coorslocal[k+1][j][i].z - coorslocal[k-1][j][i].z);
96: }
97: }
98: }
99: }
100: DMDAVecRestoreArray(cda,global,&coors);
101: DMDAVecRestoreArrayRead(cda,local,&coorslocal);
102: DMGlobalToLocalBegin(cda,global,INSERT_VALUES,local);
103: DMGlobalToLocalEnd(cda,global,INSERT_VALUES,local);
104: return(0);
105: }
107: int main(int argc,char **argv)
108: {
109: PetscInt M = 5,N = 4,P = 3, m = PETSC_DECIDE,n = PETSC_DECIDE,p = PETSC_DECIDE,dim = 1;
110: PetscErrorCode ierr;
111: DM dac,daf;
112: DMBoundaryType bx = DM_BOUNDARY_NONE,by=DM_BOUNDARY_NONE,bz=DM_BOUNDARY_NONE;
113: DMDAStencilType stype = DMDA_STENCIL_BOX;
114: Mat A;
116: PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
117: PetscOptionsGetInt(NULL,NULL,"-M",&M,NULL);
118: PetscOptionsGetInt(NULL,NULL,"-N",&N,NULL);
119: PetscOptionsGetInt(NULL,NULL,"-P",&P,NULL);
120: PetscOptionsGetInt(NULL,NULL,"-m",&m,NULL);
121: PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL);
122: PetscOptionsGetInt(NULL,NULL,"-p",&p,NULL);
123: PetscOptionsGetInt(NULL,NULL,"-dim",&dim,NULL);
125: /* Create distributed array and get vectors */
126: if (dim == 1) {
127: DMDACreate1d(PETSC_COMM_WORLD,bx,M,1,1,NULL,&dac);
128: } else if (dim == 2) {
129: DMDACreate2d(PETSC_COMM_WORLD,bx,by,stype,M,N,PETSC_DECIDE,PETSC_DECIDE,1,1,NULL,NULL,&dac);
130: } else if (dim == 3) {
131: DMDACreate3d(PETSC_COMM_WORLD,bx,by,bz,stype,M,N,P,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,1,1,NULL,NULL,NULL,&dac);
132: } else SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"dim must be 1,2, or 3");
133: DMSetFromOptions(dac);
134: DMSetUp(dac);
136: DMRefine(dac,PETSC_COMM_WORLD,&daf);
138: DMDASetUniformCoordinates(dac,0.0,1.0,0.0,1.0,0.0,1.0);
139: if (dim == 1) {
140: SetCoordinates1d(daf);
141: } else if (dim == 2) {
142: SetCoordinates2d(daf);
143: } else if (dim == 3) {
144: SetCoordinates3d(daf);
145: }
146: DMCreateInterpolation(dac,daf,&A,0);
148: /* Free memory */
149: DMDestroy(&dac);
150: DMDestroy(&daf);
151: MatDestroy(&A);
152: PetscFinalize();
153: return ierr;
154: }
158: /*TEST
160: test:
161: nsize: 3
162: args: -mat_view
164: test:
165: suffix: 2
166: nsize: 3
167: args: -mat_view -dim 2
169: test:
170: suffix: 3
171: nsize: 3
172: args: -mat_view -dim 3
174: TEST*/