Actual source code: ex23.c
1: static char help[] = "Test modifying DMStag coordinates, when represented as a product of 1d coordinate arrays\n\n";
3: #include <petscdm.h>
4: #include <petscdmstag.h>
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;
12: PetscFunctionBeginUser;
13: PetscCall(PetscInitialize(&argc, &argv, NULL, help));
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));
20: PetscCall(DMStagGetCorners(dm, &start[0], &start[1], &start[2], &n[0], &n[1], &n[2], &nExtra[0], &nExtra[1], &nExtra[2]));
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));
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;
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: }
80: PetscCall(DMDestroy(&dm));
81: PetscCall(PetscFinalize());
82: return 0;
83: }
85: /*TEST
87: test:
88: suffix: 1
89: nsize: 1
91: test:
92: suffix: 2
93: nsize: 2
95: TEST*/