1 2 static char help[] = "Tests DMCreateInterpolation() for nonuniform DMDA coordinates.\n\n"; 3 4 #include <petscdm.h> 5 #include <petscdmda.h> 6 7 PetscErrorCode SetCoordinates1d(DM da) 8 { 9 PetscErrorCode ierr; 10 PetscInt i,start,m; 11 Vec local,global; 12 PetscScalar *coors,*coorslocal; 13 DM cda; 14 15 PetscFunctionBeginUser; 16 ierr = DMDASetUniformCoordinates(da,0.0,1.0,0.0,1.0,0.0,1.0);CHKERRQ(ierr); 17 ierr = DMGetCoordinateDM(da,&cda);CHKERRQ(ierr); 18 ierr = DMGetCoordinates(da,&global);CHKERRQ(ierr); 19 ierr = DMGetCoordinatesLocal(da,&local);CHKERRQ(ierr); 20 ierr = DMDAVecGetArray(cda,global,&coors);CHKERRQ(ierr); 21 ierr = DMDAVecGetArrayRead(cda,local,&coorslocal);CHKERRQ(ierr); 22 ierr = DMDAGetCorners(cda,&start,0,0,&m,0,0);CHKERRQ(ierr); 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 ierr = DMDAVecRestoreArray(cda,global,&coors);CHKERRQ(ierr); 29 ierr = DMDAVecRestoreArrayRead(cda,local,&coorslocal);CHKERRQ(ierr); 30 ierr = DMGlobalToLocalBegin(cda,global,INSERT_VALUES,local);CHKERRQ(ierr); 31 ierr = DMGlobalToLocalEnd(cda,global,INSERT_VALUES,local);CHKERRQ(ierr); 32 PetscFunctionReturn(0); 33 } 34 35 PetscErrorCode SetCoordinates2d(DM da) 36 { 37 PetscErrorCode ierr; 38 PetscInt i,j,mstart,m,nstart,n; 39 Vec local,global; 40 DMDACoor2d **coors,**coorslocal; 41 DM cda; 42 43 PetscFunctionBeginUser; 44 ierr = DMDASetUniformCoordinates(da,0.0,1.0,0.0,1.0,0.0,1.0);CHKERRQ(ierr); 45 ierr = DMGetCoordinateDM(da,&cda);CHKERRQ(ierr); 46 ierr = DMGetCoordinates(da,&global);CHKERRQ(ierr); 47 ierr = DMGetCoordinatesLocal(da,&local);CHKERRQ(ierr); 48 ierr = DMDAVecGetArray(cda,global,&coors);CHKERRQ(ierr); 49 ierr = DMDAVecGetArrayRead(cda,local,&coorslocal);CHKERRQ(ierr); 50 ierr = DMDAGetCorners(cda,&mstart,&nstart,0,&m,&n,0);CHKERRQ(ierr); 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 ierr = DMDAVecRestoreArray(cda,global,&coors);CHKERRQ(ierr); 62 ierr = DMDAVecRestoreArrayRead(cda,local,&coorslocal);CHKERRQ(ierr); 63 64 ierr = DMGlobalToLocalBegin(cda,global,INSERT_VALUES,local);CHKERRQ(ierr); 65 ierr = DMGlobalToLocalEnd(cda,global,INSERT_VALUES,local);CHKERRQ(ierr); 66 PetscFunctionReturn(0); 67 } 68 69 PetscErrorCode SetCoordinates3d(DM da) 70 { 71 PetscErrorCode ierr; 72 PetscInt i,j,mstart,m,nstart,n,pstart,p,k; 73 Vec local,global; 74 DMDACoor3d ***coors,***coorslocal; 75 DM cda; 76 77 PetscFunctionBeginUser; 78 ierr = DMDASetUniformCoordinates(da,0.0,1.0,0.0,1.0,0.0,1.0);CHKERRQ(ierr); 79 ierr = DMGetCoordinateDM(da,&cda);CHKERRQ(ierr); 80 ierr = DMGetCoordinates(da,&global);CHKERRQ(ierr); 81 ierr = DMGetCoordinatesLocal(da,&local);CHKERRQ(ierr); 82 ierr = DMDAVecGetArray(cda,global,&coors);CHKERRQ(ierr); 83 ierr = DMDAVecGetArrayRead(cda,local,&coorslocal);CHKERRQ(ierr); 84 ierr = DMDAGetCorners(cda,&mstart,&nstart,&pstart,&m,&n,&p);CHKERRQ(ierr); 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 ierr = DMDAVecRestoreArray(cda,global,&coors);CHKERRQ(ierr); 101 ierr = DMDAVecRestoreArrayRead(cda,local,&coorslocal);CHKERRQ(ierr); 102 ierr = DMGlobalToLocalBegin(cda,global,INSERT_VALUES,local);CHKERRQ(ierr); 103 ierr = DMGlobalToLocalEnd(cda,global,INSERT_VALUES,local);CHKERRQ(ierr); 104 PetscFunctionReturn(0); 105 } 106 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; 115 116 ierr = PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr; 117 ierr = PetscOptionsGetInt(NULL,NULL,"-M",&M,NULL);CHKERRQ(ierr); 118 ierr = PetscOptionsGetInt(NULL,NULL,"-N",&N,NULL);CHKERRQ(ierr); 119 ierr = PetscOptionsGetInt(NULL,NULL,"-P",&P,NULL);CHKERRQ(ierr); 120 ierr = PetscOptionsGetInt(NULL,NULL,"-m",&m,NULL);CHKERRQ(ierr); 121 ierr = PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL);CHKERRQ(ierr); 122 ierr = PetscOptionsGetInt(NULL,NULL,"-p",&p,NULL);CHKERRQ(ierr); 123 ierr = PetscOptionsGetInt(NULL,NULL,"-dim",&dim,NULL);CHKERRQ(ierr); 124 125 /* Create distributed array and get vectors */ 126 if (dim == 1) { 127 ierr = DMDACreate1d(PETSC_COMM_WORLD,bx,M,1,1,NULL,&dac);CHKERRQ(ierr); 128 } else if (dim == 2) { 129 ierr = DMDACreate2d(PETSC_COMM_WORLD,bx,by,stype,M,N,PETSC_DECIDE,PETSC_DECIDE,1,1,NULL,NULL,&dac);CHKERRQ(ierr); 130 } else if (dim == 3) { 131 ierr = DMDACreate3d(PETSC_COMM_WORLD,bx,by,bz,stype,M,N,P,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,1,1,NULL,NULL,NULL,&dac);CHKERRQ(ierr); 132 } else SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"dim must be 1,2, or 3"); 133 ierr = DMSetFromOptions(dac);CHKERRQ(ierr); 134 ierr = DMSetUp(dac);CHKERRQ(ierr); 135 136 ierr = DMRefine(dac,PETSC_COMM_WORLD,&daf);CHKERRQ(ierr); 137 138 ierr = DMDASetUniformCoordinates(dac,0.0,1.0,0.0,1.0,0.0,1.0);CHKERRQ(ierr); 139 if (dim == 1) { 140 ierr = SetCoordinates1d(daf);CHKERRQ(ierr); 141 } else if (dim == 2) { 142 ierr = SetCoordinates2d(daf);CHKERRQ(ierr); 143 } else if (dim == 3) { 144 ierr = SetCoordinates3d(daf);CHKERRQ(ierr); 145 } 146 ierr = DMCreateInterpolation(dac,daf,&A,0);CHKERRQ(ierr); 147 148 /* Free memory */ 149 ierr = DMDestroy(&dac);CHKERRQ(ierr); 150 ierr = DMDestroy(&daf);CHKERRQ(ierr); 151 ierr = MatDestroy(&A);CHKERRQ(ierr); 152 ierr = PetscFinalize(); 153 return ierr; 154 } 155 156 /*TEST 157 158 test: 159 nsize: 3 160 args: -mat_view 161 162 test: 163 suffix: 2 164 nsize: 3 165 args: -mat_view -dim 2 166 167 test: 168 suffix: 3 169 nsize: 3 170 args: -mat_view -dim 3 171 172 TEST*/ 173