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 PetscInt i,start,m; 10 Vec local,global; 11 PetscScalar *coors,*coorslocal; 12 DM cda; 13 14 PetscFunctionBeginUser; 15 PetscCall(DMDASetUniformCoordinates(da,0.0,1.0,0.0,1.0,0.0,1.0)); 16 PetscCall(DMGetCoordinateDM(da,&cda)); 17 PetscCall(DMGetCoordinates(da,&global)); 18 PetscCall(DMGetCoordinatesLocal(da,&local)); 19 PetscCall(DMDAVecGetArray(cda,global,&coors)); 20 PetscCall(DMDAVecGetArrayRead(cda,local,&coorslocal)); 21 PetscCall(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 PetscCall(DMDAVecRestoreArray(cda,global,&coors)); 28 PetscCall(DMDAVecRestoreArrayRead(cda,local,&coorslocal)); 29 PetscCall(DMGlobalToLocalBegin(cda,global,INSERT_VALUES,local)); 30 PetscCall(DMGlobalToLocalEnd(cda,global,INSERT_VALUES,local)); 31 PetscFunctionReturn(0); 32 } 33 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; 40 41 PetscFunctionBeginUser; 42 PetscCall(DMDASetUniformCoordinates(da,0.0,1.0,0.0,1.0,0.0,1.0)); 43 PetscCall(DMGetCoordinateDM(da,&cda)); 44 PetscCall(DMGetCoordinates(da,&global)); 45 PetscCall(DMGetCoordinatesLocal(da,&local)); 46 PetscCall(DMDAVecGetArray(cda,global,&coors)); 47 PetscCall(DMDAVecGetArrayRead(cda,local,&coorslocal)); 48 PetscCall(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 PetscCall(DMDAVecRestoreArray(cda,global,&coors)); 60 PetscCall(DMDAVecRestoreArrayRead(cda,local,&coorslocal)); 61 62 PetscCall(DMGlobalToLocalBegin(cda,global,INSERT_VALUES,local)); 63 PetscCall(DMGlobalToLocalEnd(cda,global,INSERT_VALUES,local)); 64 PetscFunctionReturn(0); 65 } 66 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; 73 74 PetscFunctionBeginUser; 75 PetscCall(DMDASetUniformCoordinates(da,0.0,1.0,0.0,1.0,0.0,1.0)); 76 PetscCall(DMGetCoordinateDM(da,&cda)); 77 PetscCall(DMGetCoordinates(da,&global)); 78 PetscCall(DMGetCoordinatesLocal(da,&local)); 79 PetscCall(DMDAVecGetArray(cda,global,&coors)); 80 PetscCall(DMDAVecGetArrayRead(cda,local,&coorslocal)); 81 PetscCall(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 PetscCall(DMDAVecRestoreArray(cda,global,&coors)); 98 PetscCall(DMDAVecRestoreArrayRead(cda,local,&coorslocal)); 99 PetscCall(DMGlobalToLocalBegin(cda,global,INSERT_VALUES,local)); 100 PetscCall(DMGlobalToLocalEnd(cda,global,INSERT_VALUES,local)); 101 PetscFunctionReturn(0); 102 } 103 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; 111 112 PetscFunctionBeginUser; 113 PetscCall(PetscInitialize(&argc,&argv,(char*)0,help)); 114 PetscCall(PetscOptionsGetInt(NULL,NULL,"-M",&M,NULL)); 115 PetscCall(PetscOptionsGetInt(NULL,NULL,"-N",&N,NULL)); 116 PetscCall(PetscOptionsGetInt(NULL,NULL,"-P",&P,NULL)); 117 PetscCall(PetscOptionsGetInt(NULL,NULL,"-m",&m,NULL)); 118 PetscCall(PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL)); 119 PetscCall(PetscOptionsGetInt(NULL,NULL,"-p",&p,NULL)); 120 PetscCall(PetscOptionsGetInt(NULL,NULL,"-dim",&dim,NULL)); 121 122 /* Create distributed array and get vectors */ 123 if (dim == 1) { 124 PetscCall(DMDACreate1d(PETSC_COMM_WORLD,bx,M,1,1,NULL,&dac)); 125 } else if (dim == 2) { 126 PetscCall(DMDACreate2d(PETSC_COMM_WORLD,bx,by,stype,M,N,PETSC_DECIDE,PETSC_DECIDE,1,1,NULL,NULL,&dac)); 127 } else if (dim == 3) { 128 PetscCall(DMDACreate3d(PETSC_COMM_WORLD,bx,by,bz,stype,M,N,P,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,1,1,NULL,NULL,NULL,&dac)); 129 } else SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"dim must be 1,2, or 3"); 130 PetscCall(DMSetFromOptions(dac)); 131 PetscCall(DMSetUp(dac)); 132 133 PetscCall(DMRefine(dac,PETSC_COMM_WORLD,&daf)); 134 135 PetscCall(DMDASetUniformCoordinates(dac,0.0,1.0,0.0,1.0,0.0,1.0)); 136 if (dim == 1) { 137 PetscCall(SetCoordinates1d(daf)); 138 } else if (dim == 2) { 139 PetscCall(SetCoordinates2d(daf)); 140 } else if (dim == 3) { 141 PetscCall(SetCoordinates3d(daf)); 142 } 143 PetscCall(DMCreateInterpolation(dac,daf,&A,0)); 144 145 /* Free memory */ 146 PetscCall(DMDestroy(&dac)); 147 PetscCall(DMDestroy(&daf)); 148 PetscCall(MatDestroy(&A)); 149 PetscCall(PetscFinalize()); 150 return 0; 151 } 152 153 /*TEST 154 155 test: 156 nsize: 3 157 args: -mat_view 158 159 test: 160 suffix: 2 161 nsize: 3 162 args: -mat_view -dim 2 163 164 test: 165 suffix: 3 166 nsize: 3 167 args: -mat_view -dim 3 168 169 TEST*/ 170