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;
15: DMDASetUniformCoordinates(da, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0);
16: DMGetCoordinateDM(da, &cda);
17: DMGetCoordinates(da, &global);
18: DMGetCoordinatesLocal(da, &local);
19: DMDAVecGetArray(cda, global, &coors);
20: DMDAVecGetArrayRead(cda, local, &coorslocal);
21: 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: DMDAVecRestoreArray(cda, global, &coors);
26: DMDAVecRestoreArrayRead(cda, local, &coorslocal);
27: DMGlobalToLocalBegin(cda, global, INSERT_VALUES, local);
28: DMGlobalToLocalEnd(cda, global, INSERT_VALUES, local);
29: return 0;
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;
40: DMDASetUniformCoordinates(da, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0);
41: DMGetCoordinateDM(da, &cda);
42: DMGetCoordinates(da, &global);
43: DMGetCoordinatesLocal(da, &local);
44: DMDAVecGetArray(cda, global, &coors);
45: DMDAVecGetArrayRead(cda, local, &coorslocal);
46: 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: DMDAVecRestoreArray(cda, global, &coors);
54: DMDAVecRestoreArrayRead(cda, local, &coorslocal);
56: DMGlobalToLocalBegin(cda, global, INSERT_VALUES, local);
57: DMGlobalToLocalEnd(cda, global, INSERT_VALUES, local);
58: return 0;
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;
69: DMDASetUniformCoordinates(da, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0);
70: DMGetCoordinateDM(da, &cda);
71: DMGetCoordinates(da, &global);
72: DMGetCoordinatesLocal(da, &local);
73: DMDAVecGetArray(cda, global, &coors);
74: DMDAVecGetArrayRead(cda, local, &coorslocal);
75: 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: DMDAVecRestoreArray(cda, global, &coors);
86: DMDAVecRestoreArrayRead(cda, local, &coorslocal);
87: DMGlobalToLocalBegin(cda, global, INSERT_VALUES, local);
88: DMGlobalToLocalEnd(cda, global, INSERT_VALUES, local);
89: return 0;
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;
101: PetscInitialize(&argc, &argv, (char *)0, help);
102: PetscOptionsGetInt(NULL, NULL, "-M", &M, NULL);
103: PetscOptionsGetInt(NULL, NULL, "-N", &N, NULL);
104: PetscOptionsGetInt(NULL, NULL, "-P", &P, NULL);
105: PetscOptionsGetInt(NULL, NULL, "-m", &m, NULL);
106: PetscOptionsGetInt(NULL, NULL, "-n", &n, NULL);
107: PetscOptionsGetInt(NULL, NULL, "-p", &p, NULL);
108: PetscOptionsGetInt(NULL, NULL, "-dim", &dim, NULL);
110: /* Create distributed array and get vectors */
111: if (dim == 1) {
112: DMDACreate1d(PETSC_COMM_WORLD, bx, M, 1, 1, NULL, &dac);
113: } else if (dim == 2) {
114: DMDACreate2d(PETSC_COMM_WORLD, bx, by, stype, M, N, PETSC_DECIDE, PETSC_DECIDE, 1, 1, NULL, NULL, &dac);
115: } else if (dim == 3) {
116: 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: DMSetFromOptions(dac);
119: DMSetUp(dac);
121: DMRefine(dac, PETSC_COMM_WORLD, &daf);
123: DMDASetUniformCoordinates(dac, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0);
124: if (dim == 1) {
125: SetCoordinates1d(daf);
126: } else if (dim == 2) {
127: SetCoordinates2d(daf);
128: } else if (dim == 3) {
129: SetCoordinates3d(daf);
130: }
131: DMCreateInterpolation(dac, daf, &A, 0);
133: /* Free memory */
134: DMDestroy(&dac);
135: DMDestroy(&daf);
136: MatDestroy(&A);
137: 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*/