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