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