Actual source code: ex19.c
1: static char help[] = "(Partially) test DMStag default interpolation, 2d faces-only.\n\n";
3: #include <petscdm.h>
4: #include <petscdmstag.h>
5: #include <petscksp.h>
7: PetscErrorCode CreateSystem(DM dm, Mat *A, Vec *b);
9: int main(int argc, char **argv)
10: {
11: DM dm, dmCoarse;
12: Mat Ai;
14: PetscFunctionBeginUser;
15: PetscCall(PetscInitialize(&argc, &argv, NULL, help));
16: PetscCall(DMStagCreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, 2, 4, PETSC_DECIDE, PETSC_DECIDE, 0, 1, 1, DMSTAG_STENCIL_BOX, 1, NULL, NULL, &dm));
17: PetscCall(DMSetFromOptions(dm));
18: PetscCall(DMSetUp(dm));
19: PetscCall(DMCoarsen(dm, MPI_COMM_NULL, &dmCoarse));
20: PetscCall(DMCreateInterpolation(dmCoarse, dm, &Ai, NULL));
22: /* See what happens to a constant value on each sub-grid */
23: {
24: Vec localCoarse, globalCoarse, globalFine, localFine;
25: PetscCall(DMGetGlobalVector(dm, &globalFine));
26: PetscCall(DMGetGlobalVector(dmCoarse, &globalCoarse));
27: PetscCall(DMGetLocalVector(dmCoarse, &localCoarse));
28: PetscCall(DMGetLocalVector(dm, &localFine));
29: PetscCall(VecSet(localCoarse, -1.0));
30: PetscCall(VecSet(localFine, -1.0));
31: {
32: PetscInt i, j, startx, starty, nx, ny, extrax, extray;
33: PetscInt p, vx, vy;
34: PetscScalar ***arr;
35: PetscCall(DMStagGetCorners(dmCoarse, &startx, &starty, NULL, &nx, &ny, NULL, &extrax, &extray, NULL));
36: PetscCall(DMStagVecGetArray(dmCoarse, localCoarse, &arr));
37: PetscCall(DMStagGetLocationSlot(dmCoarse, DMSTAG_LEFT, 0, &vx));
38: PetscCall(DMStagGetLocationSlot(dmCoarse, DMSTAG_DOWN, 0, &vy));
39: PetscCall(DMStagGetLocationSlot(dmCoarse, DMSTAG_ELEMENT, 0, &p));
40: for (j = starty; j < starty + ny + extray; ++j) {
41: for (i = startx; i < startx + nx + extrax; ++i) {
42: arr[j][i][vy] = (i < startx + nx) ? 10.0 : -1;
43: arr[j][i][vx] = (j < starty + ny) ? 20.0 : -1;
44: arr[j][i][p] = (i < startx + nx) && (j < starty + ny) ? 30.0 : -1;
45: }
46: }
47: PetscCall(DMStagVecRestoreArray(dmCoarse, localCoarse, &arr));
48: }
49: PetscCall(DMLocalToGlobal(dmCoarse, localCoarse, INSERT_VALUES, globalCoarse));
50: PetscCall(MatInterpolate(Ai, globalCoarse, globalFine));
51: PetscCall(DMGlobalToLocal(dm, globalFine, INSERT_VALUES, localFine));
52: {
53: PetscInt i, j, startx, starty, nx, ny, extrax, extray;
54: PetscInt p, vx, vy;
55: PetscScalar ***arr;
56: PetscCall(DMStagGetCorners(dm, &startx, &starty, NULL, &nx, &ny, NULL, &extrax, &extray, NULL));
57: PetscCall(DMStagVecGetArrayRead(dm, localFine, &arr));
58: PetscCall(DMStagGetLocationSlot(dm, DMSTAG_LEFT, 0, &vx));
59: PetscCall(DMStagGetLocationSlot(dm, DMSTAG_DOWN, 0, &vy));
60: PetscCall(DMStagGetLocationSlot(dm, DMSTAG_ELEMENT, 0, &p));
61: for (j = starty; j < starty + ny + extray; ++j) {
62: for (i = startx; i < startx + nx + extrax; ++i) {
63: const PetscScalar expected_vy = (i < startx + nx) ? 10.0 : -1;
64: const PetscScalar expected_vx = (j < starty + ny) ? 20.0 : -1;
65: const PetscScalar expected_p = (i < startx + nx) && (j < starty + ny) ? 30.0 : -1;
66: if (arr[j][i][vy] != expected_vy) PetscCall(PetscPrintf(PETSC_COMM_SELF, "wrong %" PetscInt_FMT " %" PetscInt_FMT "\n", i, j));
67: if (arr[j][i][vx] != expected_vx) PetscCall(PetscPrintf(PETSC_COMM_SELF, "wrong %" PetscInt_FMT " %" PetscInt_FMT "\n", i, j));
68: if (arr[j][i][p] != expected_p) PetscCall(PetscPrintf(PETSC_COMM_SELF, "wrong %" PetscInt_FMT " %" PetscInt_FMT "\n", i, j));
69: }
70: }
71: PetscCall(DMStagVecRestoreArrayRead(dm, localFine, &arr));
72: }
73: PetscCall(DMRestoreLocalVector(dmCoarse, &localCoarse));
74: PetscCall(DMRestoreLocalVector(dm, &localFine));
75: PetscCall(DMRestoreGlobalVector(dmCoarse, &globalCoarse));
76: PetscCall(DMRestoreGlobalVector(dm, &globalFine));
77: }
79: PetscCall(MatDestroy(&Ai));
80: PetscCall(DMDestroy(&dm));
81: PetscCall(DMDestroy(&dmCoarse));
82: PetscCall(PetscFinalize());
83: return 0;
84: }
86: /*TEST
88: test:
89: suffix: 1
90: nsize: 1
91: args:
93: test:
94: suffix: 2
95: nsize: 4
96: args: -stag_grid_x 8 -stag_grid_y 4
98: TEST*/