Actual source code: ex3.c


  2: static char help[] = "Tests DMCreateInterpolation() for nonuniform DMDA coordinates.\n\n";

  4: #include <petscdm.h>
  5: #include <petscdmda.h>

  7: PetscErrorCode SetCoordinates1d(DM da)
  8: {
  9:   PetscInt     i, start, m;
 10:   Vec          local, global;
 11:   PetscScalar *coors, *coorslocal;
 12:   DM           cda;

 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: }

 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;

 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));

 56:   PetscCall(DMGlobalToLocalBegin(cda, global, INSERT_VALUES, local));
 57:   PetscCall(DMGlobalToLocalEnd(cda, global, INSERT_VALUES, local));
 58:   PetscFunctionReturn(PETSC_SUCCESS);
 59: }

 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;

 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: }

 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;

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));

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));

121:   PetscCall(DMRefine(dac, PETSC_COMM_WORLD, &daf));

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));

133:   /* Free memory */
134:   PetscCall(DMDestroy(&dac));
135:   PetscCall(DMDestroy(&daf));
136:   PetscCall(MatDestroy(&A));
137:   PetscCall(PetscFinalize());
138:   return 0;
139: }

141: /*TEST

143:    test:
144:       nsize: 3
145:       args: -mat_view

147:    test:
148:       suffix: 2
149:       nsize: 3
150:       args: -mat_view -dim 2

152:    test:
153:       suffix: 3
154:       nsize: 3
155:       args: -mat_view -dim 3

157: TEST*/