1 static char help[] = "Tests DMCreateInterpolation() for nonuniform DMDA coordinates.\n\n";
2
3 #include <petscdm.h>
4 #include <petscdmda.h>
5
SetCoordinates1d(DM da)6 PetscErrorCode SetCoordinates1d(DM da)
7 {
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(PETSC_SUCCESS);
29 }
30
SetCoordinates2d(DM da)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;
37
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));
54
55 PetscCall(DMGlobalToLocalBegin(cda, global, INSERT_VALUES, local));
56 PetscCall(DMGlobalToLocalEnd(cda, global, INSERT_VALUES, local));
57 PetscFunctionReturn(PETSC_SUCCESS);
58 }
59
SetCoordinates3d(DM da)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;
66
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 }
90
main(int argc,char ** argv)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;
98
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));
108
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));
119
120 PetscCall(DMRefine(dac, PETSC_COMM_WORLD, &daf));
121
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));
131
132 /* Free memory */
133 PetscCall(DMDestroy(&dac));
134 PetscCall(DMDestroy(&daf));
135 PetscCall(MatDestroy(&A));
136 PetscCall(PetscFinalize());
137 return 0;
138 }
139
140 /*TEST
141
142 test:
143 nsize: 3
144 args: -mat_view
145
146 test:
147 suffix: 2
148 nsize: 3
149 args: -mat_view -dim 2
150
151 test:
152 suffix: 3
153 nsize: 3
154 args: -mat_view -dim 3
155
156 TEST*/
157