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*/