Actual source code: ex3.c
1: static char help[] = "Tests DMCreateInterpolation() for nonuniform DMDA coordinates.\n\n";
3: #include <petscdm.h>
4: #include <petscdmda.h>
6: PetscErrorCode SetCoordinates1d(DM da)
7: {
8: PetscInt i, start, m;
9: Vec local, global;
10: PetscScalar *coors, *coorslocal;
11: DM cda;
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: }
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;
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));
55: PetscCall(DMGlobalToLocalBegin(cda, global, INSERT_VALUES, local));
56: PetscCall(DMGlobalToLocalEnd(cda, global, INSERT_VALUES, local));
57: PetscFunctionReturn(PETSC_SUCCESS);
58: }
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;
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: }
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;
99: PetscFunctionBeginUser;
100: PetscCall(PetscInitialize(&argc, &argv, NULL, 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));
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));
120: PetscCall(DMRefine(dac, PETSC_COMM_WORLD, &daf));
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));
132: /* Free memory */
133: PetscCall(DMDestroy(&dac));
134: PetscCall(DMDestroy(&daf));
135: PetscCall(MatDestroy(&A));
136: PetscCall(PetscFinalize());
137: return 0;
138: }
140: /*TEST
142: test:
143: nsize: 3
144: args: -mat_view
146: test:
147: suffix: 2
148: nsize: 3
149: args: -mat_view -dim 2
151: test:
152: suffix: 3
153: nsize: 3
154: args: -mat_view -dim 3
156: TEST*/