1 static char help[] = "Test modifying DMStag coordinates, when represented as a product of 1d coordinate arrays\n\n";
2
3 #include <petscdm.h>
4 #include <petscdmstag.h>
5
main(int argc,char ** argv)6 int main(int argc, char **argv)
7 {
8 DM dm, cdm;
9 PetscInt ex, ey, ez, n[3], start[3], nExtra[3], iNext, iPrev, iCenter, d, round;
10 PetscScalar **cArrX, **cArrY, **cArrZ;
11
12 PetscFunctionBeginUser;
13 PetscCall(PetscInitialize(&argc, &argv, NULL, help));
14
15 PetscCall(DMStagCreate3d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_GHOSTED, DM_BOUNDARY_PERIODIC, 4, 3, 2, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE, 1, 1, 1, 1, DMSTAG_STENCIL_BOX, 2, NULL, NULL, NULL, &dm));
16 PetscCall(DMSetFromOptions(dm));
17 PetscCall(DMSetUp(dm));
18 PetscCall(DMStagSetUniformCoordinatesProduct(dm, -1.0, 0.0, -2.0, 0.0, -3.0, 0.0));
19
20 PetscCall(DMStagGetCorners(dm, &start[0], &start[1], &start[2], &n[0], &n[1], &n[2], &nExtra[0], &nExtra[1], &nExtra[2]));
21
22 for (round = 1; round <= 2; ++round) {
23 PetscCall(DMStagGetProductCoordinateArrays(dm, &cArrX, &cArrY, &cArrZ));
24 PetscCall(DMStagGetProductCoordinateLocationSlot(dm, DMSTAG_LEFT, &iPrev));
25 PetscCall(DMStagGetProductCoordinateLocationSlot(dm, DMSTAG_RIGHT, &iNext));
26 PetscCall(DMStagGetProductCoordinateLocationSlot(dm, DMSTAG_ELEMENT, &iCenter));
27 if (round == 1) {
28 /* On first round, do a stretching operation */
29 for (ex = start[0]; ex < start[0] + n[0]; ++ex) {
30 cArrX[ex][iPrev] *= 1.1;
31 cArrX[ex][iNext] = cArrX[ex][iPrev] + 0.1;
32 cArrX[ex][iCenter] = 0.5 * (cArrX[ex][iPrev] + cArrX[ex][iNext]);
33 }
34 for (ey = start[1]; ey < start[1] + n[1]; ++ey) {
35 cArrY[ey][iPrev] *= 1.1;
36 cArrY[ey][iNext] = cArrY[ey][iPrev] + 0.1;
37 cArrY[ey][iCenter] = 0.5 * (cArrY[ey][iPrev] + cArrY[ey][iNext]);
38 }
39 for (ez = start[2]; ez < start[2] + n[2]; ++ez) {
40 cArrZ[ez][iPrev] *= 1.1;
41 cArrZ[ez][iNext] = cArrZ[ez][iPrev] + 0.1;
42 cArrZ[ez][iCenter] = 0.5 * (cArrZ[ez][iPrev] + cArrZ[ez][iNext]);
43 }
44 } else {
45 /* On second round, set everything to 2.0 */
46 for (ex = start[0]; ex < start[0] + n[0]; ++ex) {
47 cArrX[ex][iPrev] = 2.0;
48 cArrX[ex][iNext] = 2.0;
49 cArrX[ex][iCenter] = 2.0;
50 }
51 for (ey = start[1]; ey < start[1] + n[1]; ++ey) {
52 cArrY[ey][iPrev] = 2.0;
53 cArrY[ey][iNext] = 2.0;
54 cArrY[ey][iCenter] = 2.0;
55 }
56 for (ez = start[2]; ez < start[2] + n[2]; ++ez) {
57 cArrZ[ez][iPrev] = 2.0;
58 cArrZ[ez][iNext] = 2.0;
59 cArrZ[ez][iCenter] = 2.0;
60 }
61 }
62 PetscCall(DMStagRestoreProductCoordinateArrays(dm, &cArrX, &cArrY, &cArrZ));
63
64 /* View the global coordinates, after explicitly calling a local-global scatter */
65 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "####### Round %" PetscInt_FMT " #######\n", round));
66 PetscCall(DMGetCoordinateDM(dm, &cdm));
67 for (d = 0; d < 3; ++d) {
68 DM subdm;
69 Vec coor, coor_local;
70
71 PetscCall(DMProductGetDM(cdm, d, &subdm));
72 PetscCall(DMGetCoordinates(subdm, &coor));
73 PetscCall(DMGetCoordinatesLocal(subdm, &coor_local));
74 PetscCall(DMLocalToGlobal(subdm, coor_local, INSERT_VALUES, coor));
75 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Coordinates dim %" PetscInt_FMT ":\n", d));
76 PetscCall(VecView(coor, PETSC_VIEWER_STDOUT_WORLD));
77 }
78 }
79
80 PetscCall(DMDestroy(&dm));
81 PetscCall(PetscFinalize());
82 return 0;
83 }
84
85 /*TEST
86
87 test:
88 suffix: 1
89 nsize: 1
90
91 test:
92 suffix: 2
93 nsize: 2
94
95 TEST*/
96