Actual source code: ex21.c
1: static const char help[] = "Test DMCreateInjection() for mapping coordinates in 3D";
3: #include <petscvec.h>
4: #include <petscmat.h>
5: #include <petscdm.h>
6: #include <petscdmda.h>
8: PetscErrorCode test1_DAInjection3d(PetscInt mx, PetscInt my, PetscInt mz)
9: {
10: DM dac, daf;
11: PetscViewer vv;
12: Vec ac, af;
13: PetscInt periodicity;
14: DMBoundaryType bx, by, bz;
16: PetscFunctionBeginUser;
17: bx = DM_BOUNDARY_NONE;
18: by = DM_BOUNDARY_NONE;
19: bz = DM_BOUNDARY_NONE;
21: periodicity = 0;
23: PetscCall(PetscOptionsGetInt(NULL, NULL, "-periodic", &periodicity, NULL));
24: if (periodicity == 1) {
25: bx = DM_BOUNDARY_PERIODIC;
26: } else if (periodicity == 2) {
27: by = DM_BOUNDARY_PERIODIC;
28: } else if (periodicity == 3) {
29: bz = DM_BOUNDARY_PERIODIC;
30: }
32: PetscCall(DMDACreate3d(PETSC_COMM_WORLD, bx, by, bz, DMDA_STENCIL_BOX, mx + 1, my + 1, mz + 1, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE, 1, /* 1 dof */
33: 1, /* stencil = 1 */ NULL, NULL, NULL, &daf));
34: PetscCall(DMSetFromOptions(daf));
35: PetscCall(DMSetUp(daf));
37: PetscCall(DMCoarsen(daf, MPI_COMM_NULL, &dac));
39: PetscCall(DMDASetUniformCoordinates(dac, -1.0, 1.0, -1.0, 1.0, -1.0, 1.0));
40: PetscCall(DMDASetUniformCoordinates(daf, -1.0, 1.0, -1.0, 1.0, -1.0, 1.0));
42: {
43: DM cdaf, cdac;
44: Vec coordsc, coordsf, coordsf2;
45: Mat inject;
46: VecScatter vscat;
47: Mat interp;
48: PetscReal norm;
50: PetscCall(DMGetCoordinateDM(dac, &cdac));
51: PetscCall(DMGetCoordinateDM(daf, &cdaf));
53: PetscCall(DMGetCoordinates(dac, &coordsc));
54: PetscCall(DMGetCoordinates(daf, &coordsf));
56: PetscCall(DMCreateInjection(cdac, cdaf, &inject));
57: PetscCall(MatScatterGetVecScatter(inject, &vscat));
58: PetscCall(VecScatterBegin(vscat, coordsf, coordsc, INSERT_VALUES, SCATTER_FORWARD));
59: PetscCall(VecScatterEnd(vscat, coordsf, coordsc, INSERT_VALUES, SCATTER_FORWARD));
60: PetscCall(MatDestroy(&inject));
62: PetscCall(DMCreateInterpolation(cdac, cdaf, &interp, NULL));
63: PetscCall(VecDuplicate(coordsf, &coordsf2));
64: PetscCall(MatInterpolate(interp, coordsc, coordsf2));
65: PetscCall(VecAXPY(coordsf2, -1.0, coordsf));
66: PetscCall(VecNorm(coordsf2, NORM_MAX, &norm));
67: /* The fine coordinates are only reproduced in certain cases */
68: if (!bx && !by && !bz && norm > PETSC_SQRT_MACHINE_EPSILON) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Norm %g\n", (double)norm));
69: PetscCall(VecDestroy(&coordsf2));
70: PetscCall(MatDestroy(&interp));
71: }
73: if (0) {
74: PetscCall(DMCreateGlobalVector(dac, &ac));
75: PetscCall(VecZeroEntries(ac));
77: PetscCall(DMCreateGlobalVector(daf, &af));
78: PetscCall(VecZeroEntries(af));
80: PetscCall(PetscViewerASCIIOpen(PETSC_COMM_WORLD, "dac_7.vtu", &vv));
81: PetscCall(VecView(ac, vv));
82: PetscCall(PetscViewerDestroy(&vv));
84: PetscCall(PetscViewerASCIIOpen(PETSC_COMM_WORLD, "daf_7.vtu", &vv));
85: PetscCall(VecView(af, vv));
86: PetscCall(PetscViewerDestroy(&vv));
87: PetscCall(VecDestroy(&ac));
88: PetscCall(VecDestroy(&af));
89: }
90: PetscCall(DMDestroy(&dac));
91: PetscCall(DMDestroy(&daf));
92: PetscFunctionReturn(PETSC_SUCCESS);
93: }
95: int main(int argc, char **argv)
96: {
97: PetscInt mx, my, mz;
99: PetscFunctionBeginUser;
100: PetscCall(PetscInitialize(&argc, &argv, 0, help));
101: mx = 2;
102: my = 2;
103: mz = 2;
104: PetscCall(PetscOptionsGetInt(NULL, NULL, "-mx", &mx, 0));
105: PetscCall(PetscOptionsGetInt(NULL, NULL, "-my", &my, 0));
106: PetscCall(PetscOptionsGetInt(NULL, NULL, "-mz", &mz, 0));
107: PetscCall(test1_DAInjection3d(mx, my, mz));
108: PetscCall(PetscFinalize());
109: return 0;
110: }
112: /*TEST
114: test:
115: nsize: 5
116: args: -mx 30 -my 30 -mz 30 -periodic 0 -da_processors_x 5
118: test:
119: suffix: 2
120: nsize: 5
121: args: -mx 29 -my 30 -mz 30 -periodic 1 -da_processors_x 5
123: test:
124: suffix: 3
125: nsize: 5
126: args: -mx 30 -my 29 -mz 30 -periodic 2 -da_processors_x 5
128: test:
129: suffix: 4
130: nsize: 5
131: args: -mx 30 -my 30 -mz 29 -periodic 3 -da_processors_x 5
133: test:
134: suffix: 5
135: nsize: 5
136: args: -mx 30 -my 30 -mz 30 -periodic 0 -da_processors_y 5
138: test:
139: suffix: 6
140: nsize: 5
141: args: -mx 29 -my 30 -mz 30 -periodic 1 -da_processors_y 5
143: test:
144: suffix: 7
145: nsize: 5
146: args: -mx 30 -my 29 -mz 30 -periodic 2 -da_processors_y 5
148: test:
149: suffix: 8
150: nsize: 5
151: args: -mx 30 -my 30 -mz 29 -periodic 3 -da_processors_y 5
153: test:
154: suffix: 9
155: nsize: 5
156: args: -mx 30 -my 30 -mz 30 -periodic 0 -da_processors_z 5
158: test:
159: suffix: 10
160: nsize: 5
161: args: -mx 29 -my 30 -mz 30 -periodic 1 -da_processors_z 5
163: test:
164: suffix: 11
165: nsize: 5
166: args: -mx 30 -my 29 -mz 30 -periodic 2 -da_processors_z 5
168: test:
169: suffix: 12
170: nsize: 5
171: args: -mx 30 -my 30 -mz 29 -periodic 3 -da_processors_z 5
173: test:
174: suffix: 13
175: nsize: 5
176: args: -mx 30 -my 30 -mz 30 -periodic 0
178: test:
179: suffix: 14
180: nsize: 5
181: args: -mx 29 -my 30 -mz 30 -periodic 1
183: test:
184: suffix: 15
185: nsize: 5
186: args: -mx 30 -my 29 -mz 30 -periodic 2
188: test:
189: suffix: 16
190: nsize: 5
191: args: -mx 30 -my 30 -mz 29 -periodic 3
193: TEST*/