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