1 static char help[] = "Takes a patch of a large DMDA vector to one process.\n\n"; 2 3 #include <petscdm.h> 4 #include <petscdmda.h> 5 #include <petscdmpatch.h> 6 #include <petscsf.h> 7 8 typedef struct { 9 PetscScalar x, y; 10 } Field; 11 12 int main(int argc, char **argv) 13 { 14 Vec xy, sxy; 15 DM da, sda = NULL; 16 PetscSF sf; 17 PetscInt m = 10, n = 10, dof = 2; 18 MatStencil lower = {0, 3, 2, 0}, upper = {0, 7, 8, 0}; /* These are in the order of the z, y, x, logical coordinates, the fourth entry is ignored */ 19 MPI_Comm comm; 20 PetscMPIInt rank; 21 22 PetscFunctionBeginUser; 23 PetscCall(PetscInitialize(&argc, &argv, (char *)0, help)); 24 25 /* create the large DMDA and set coordinates (which we will copy down to the small DA). */ 26 PetscCall(DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_BOX, m, n, PETSC_DECIDE, PETSC_DECIDE, dof, 1, 0, 0, &da)); 27 PetscCall(DMSetFromOptions(da)); 28 PetscCall(DMSetUp(da)); 29 PetscCall(DMDASetUniformCoordinates(da, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0)); 30 /* Just as a simple example we use the coordinates as the variables in the vectors we wish to examine. */ 31 PetscCall(DMGetCoordinates(da, &xy)); 32 /* The vector entries are displayed in the "natural" ordering on the two dimensional grid; interlaced x and y with the x variable increasing more rapidly than the y */ 33 PetscCall(VecView(xy, 0)); 34 35 PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank)); 36 if (rank == 0) comm = MPI_COMM_SELF; 37 else comm = MPI_COMM_NULL; 38 39 PetscCall(DMPatchZoom(da, lower, upper, comm, &sda, NULL, &sf)); 40 if (rank == 0) { 41 PetscCall(DMCreateGlobalVector(sda, &sxy)); 42 } else { 43 PetscCall(VecCreateSeq(PETSC_COMM_SELF, 0, &sxy)); 44 } 45 /* A PetscSF can also be used as a VecScatter context */ 46 PetscCall(VecScatterBegin(sf, xy, sxy, INSERT_VALUES, SCATTER_FORWARD)); 47 PetscCall(VecScatterEnd(sf, xy, sxy, INSERT_VALUES, SCATTER_FORWARD)); 48 /* Only rank == 0 has the entries of the patch, so run code only at that rank */ 49 if (rank == 0) { 50 Field **vars; 51 DMDALocalInfo info; 52 PetscInt i, j; 53 PetscScalar sum = 0; 54 55 /* The vector entries of the patch are displayed in the "natural" ordering on the two grid; interlaced x and y with the x variable increasing more rapidly */ 56 PetscCall(VecView(sxy, PETSC_VIEWER_STDOUT_SELF)); 57 /* Compute some trivial statistic of the coordinates */ 58 PetscCall(DMDAGetLocalInfo(sda, &info)); 59 PetscCall(DMDAVecGetArray(sda, sxy, &vars)); 60 /* Loop over the patch of the entire domain */ 61 for (j = info.ys; j < info.ys + info.ym; j++) { 62 for (i = info.xs; i < info.xs + info.xm; i++) sum += vars[j][i].x; 63 } 64 PetscCall(PetscPrintf(PETSC_COMM_SELF, "The sum of the x coordinates is %g\n", (double)PetscRealPart(sum))); 65 PetscCall(DMDAVecRestoreArray(sda, sxy, &vars)); 66 } 67 68 PetscCall(VecDestroy(&sxy)); 69 PetscCall(PetscSFDestroy(&sf)); 70 PetscCall(DMDestroy(&sda)); 71 PetscCall(DMDestroy(&da)); 72 PetscCall(PetscFinalize()); 73 return 0; 74 } 75 76 /*TEST 77 78 test: 79 80 test: 81 nsize: 4 82 suffix: 2 83 84 TEST*/ 85