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