1: static char help[] = "Test DMStag 2d star stencil\n\n";
  2: #include <petscdm.h>
  3: #include <petscdmstag.h>

  5: int main(int argc, char **argv)
  6: {
  7:   DM             dm;
  8:   Vec            vec, vecLocal1, vecLocal2;
  9:   PetscScalar   *a, ***a1, ***a2, expected, sum;
 10:   PetscInt       startx, starty, nx, ny, i, j, d, is, js, dof0, dof1, dof2, dofTotal, stencilWidth, ngx, ngy;
 11:   DMBoundaryType boundaryTypex, boundaryTypey;
 12:   PetscMPIInt    rank;

 14:   PetscFunctionBeginUser;
 15:   PetscCall(PetscInitialize(&argc, &argv, NULL, help));
 16:   PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
 17:   dof0         = 1;
 18:   dof1         = 1;
 19:   dof2         = 1;
 20:   stencilWidth = 2;
 21:   PetscCall(DMStagCreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_PERIODIC, DM_BOUNDARY_PERIODIC, 4, 4, PETSC_DECIDE, PETSC_DECIDE, dof0, dof1, dof2, DMSTAG_STENCIL_STAR, stencilWidth, NULL, NULL, &dm));
 22:   PetscCall(DMSetFromOptions(dm));
 23:   PetscCall(DMSetUp(dm));
 24:   PetscCall(DMStagGetDOF(dm, &dof0, &dof1, &dof2, NULL));
 25:   dofTotal = dof0 + 2 * dof1 + dof2;
 26:   PetscCall(DMStagGetStencilWidth(dm, &stencilWidth));

 28:   PetscCall(DMCreateLocalVector(dm, &vecLocal1));
 29:   PetscCall(VecDuplicate(vecLocal1, &vecLocal2));

 31:   PetscCall(DMCreateGlobalVector(dm, &vec));
 32:   PetscCall(VecSet(vec, 1.0));
 33:   PetscCall(DMGlobalToLocalBegin(dm, vec, INSERT_VALUES, vecLocal1));
 34:   PetscCall(DMGlobalToLocalEnd(dm, vec, INSERT_VALUES, vecLocal1));

 36:   PetscCall(DMStagGetCorners(dm, &startx, &starty, NULL, &nx, &ny, NULL, NULL, NULL, NULL));
 37:   PetscCall(DMStagVecGetArrayRead(dm, vecLocal1, &a1));
 38:   PetscCall(DMStagVecGetArray(dm, vecLocal2, &a2));
 39:   for (j = starty; j < starty + ny; ++j) {
 40:     for (i = startx; i < startx + nx; ++i) {
 41:       for (d = 0; d < dofTotal; ++d) {
 42:         if (a1[j][i][d] != 1.0) PetscCall(PetscPrintf(PETSC_COMM_SELF, "[%d] Unexpected value %g (expecting %g)\n", rank, (double)PetscRealPart(a1[j][i][d]), 1.0));
 43:         a2[j][i][d] = 0.0;
 44:         for (js = -stencilWidth; js <= stencilWidth; ++js) a2[j][i][d] += a1[j + js][i][d];
 45:         for (is = -stencilWidth; is <= stencilWidth; ++is) a2[j][i][d] += a1[j][i + is][d];
 46:         a2[j][i][d] -= a1[j][i][d];
 47:       }
 48:     }
 49:   }
 50:   PetscCall(DMStagVecRestoreArrayRead(dm, vecLocal1, &a1));
 51:   PetscCall(DMStagVecRestoreArray(dm, vecLocal2, &a2));

 53:   PetscCall(DMLocalToGlobalBegin(dm, vecLocal2, INSERT_VALUES, vec));
 54:   PetscCall(DMLocalToGlobalEnd(dm, vecLocal2, INSERT_VALUES, vec));

 56:   /* For the all-periodic case, some additional checks */
 57:   PetscCall(DMStagGetBoundaryTypes(dm, &boundaryTypex, &boundaryTypey, NULL));
 58:   if (boundaryTypex == DM_BOUNDARY_PERIODIC && boundaryTypey == DM_BOUNDARY_PERIODIC) {
 59:     PetscCall(DMStagGetGhostCorners(dm, NULL, NULL, NULL, &ngx, &ngy, NULL));
 60:     expected = (ngx * ngy - 4 * stencilWidth * stencilWidth) * dofTotal;
 61:     PetscCall(VecSum(vecLocal1, &sum));
 62:     if (sum != expected) PetscCall(PetscPrintf(PETSC_COMM_SELF, "[%d] Unexpected sum of local entries %g (expected %g)\n", rank, (double)PetscRealPart(sum), (double)PetscRealPart(expected)));

 64:     PetscCall(VecGetArray(vec, &a));
 65:     expected = 1 + 4 * stencilWidth;
 66:     for (i = 0; i < ny * nx * dofTotal; ++i) {
 67:       if (a[i] != expected) PetscCall(PetscPrintf(PETSC_COMM_SELF, "[%d] Unexpected value %g (expecting %g)\n", rank, (double)PetscRealPart(a[i]), (double)PetscRealPart(expected)));
 68:     }
 69:     PetscCall(VecRestoreArray(vec, &a));
 70:   }

 72:   PetscCall(VecDestroy(&vec));
 73:   PetscCall(VecDestroy(&vecLocal1));
 74:   PetscCall(VecDestroy(&vecLocal2));
 75:   PetscCall(DMDestroy(&dm));
 76:   PetscCall(PetscFinalize());
 77:   return 0;
 78: }

 80: /*TEST

 82:    test:
 83:       suffix: 1
 84:       nsize: 4
 85:       args: -stag_ranks_x 2 -stag_ranks_y 2 -stag_stencil_width 1
 86:       output_file: output/empty.out

 88:    test:
 89:       suffix: 2
 90:       nsize: 6
 91:       args: -stag_ranks_x 3 -stag_ranks_y 2 -stag_dof_0 2 -stag_grid_x 6
 92:       output_file: output/empty.out

 94:    test:
 95:       suffix: 3
 96:       nsize: 4
 97:       args: -stag_dof_0 3 -stag_dof_1 2 -stag_dof_2 4 2 -stag_stencil_width 3 -stag_grid_x 6 -stag_grid_y 6
 98:       output_file: output/empty.out

100:    test:
101:       suffix: 4
102:       nsize: 4
103:       args: -stag_stencil_width 1 -stag_grid_x 2 -stag_grid_y 2 -stag_boundary_type_x ghosted
104:       output_file: output/empty.out

106:    test:
107:       suffix: 5
108:       nsize: 4
109:       args: -stag_stencil_width 1 -stag_grid_x 2 -stag_grid_y 2 -stag_boundary_type_y ghosted
110:       output_file: output/empty.out

112:    test:
113:       suffix: 6
114:       nsize: 4
115:       args: -stag_stencil_width 1 -stag_grid_x 3 -stag_grid_y 2 -stag_boundary_type_x ghosted -stag_boundary_type_y ghosted
116:       output_file: output/empty.out

118:    test:
119:       suffix: 7
120:       nsize: 4
121:       args: -stag_stencil_width 1 -stag_grid_x 2 -stag_grid_y 2 -stag_boundary_type_y ghosted
122:       output_file: output/empty.out

124:    test:
125:       suffix: 8
126:       nsize: 6
127:       args: -stag_stencil_width 1 -stag_grid_y 2 -stag_grid_x 19 -stag_boundary_type_y ghosted -stag_ranks_x 6
128:       output_file: output/empty.out
129: TEST*/