Actual source code: gr1.c

      1: #define PETSCDM_DLL

      3: /*
      4:    Plots vectors obtained with DACreate1d()
      5: */

     7:  #include petscda.h

     11: /*@
     12:     DASetUniformCoordinates - Sets a DA coordinates to be a uniform grid

     14:   Collective on DA

     16:   Input Parameters:
     17: +  da - the distributed array object
     18: .  xmin,xmax - extremes in the x direction
     19: .  ymin,ymax - extremes in the y direction (use PETSC_NULL for 1 dimensional problems)
     20: -  zmin,zmax - extremes in the z direction (use PETSC_NULL for 1 or 2 dimensional problems)

     22:   Level: beginner

     24: .seealso: DASetCoordinates(), DAGetCoordinates(), DACreate1d(), DACreate2d(), DACreate3d()

     26: @*/
     27: PetscErrorCode PETSCDM_DLLEXPORT DASetUniformCoordinates(DA da,PetscReal xmin,PetscReal xmax,PetscReal ymin,PetscReal ymax,PetscReal zmin,PetscReal zmax)
     28: {
     30:   PetscInt            i,j,k,M,N,P,istart,isize,jstart,jsize,kstart,ksize,dim,cnt;
     31:   PetscReal      hx,hy,hz_;
     32:   Vec            xcoor;
     33:   DAPeriodicType periodic;
     34:   PetscScalar    *coors;

     37:   if (xmax <= xmin) SETERRQ2(PETSC_ERR_ARG_INCOMP,"xmax must be larger than xmin %g %g",xmin,xmax);

     39:   DAGetInfo(da,&dim,&M,&N,&P,0,0,0,0,0,&periodic,0);
     40:   DAGetCorners(da,&istart,&jstart,&kstart,&isize,&jsize,&ksize);

     42:   if (dim == 1) {
     43:     VecCreateMPI(PETSC_COMM_WORLD,isize,PETSC_DETERMINE,&xcoor);
     44:     if (periodic == DA_NONPERIODIC) hx = (xmax-xmin)/(M-1);
     45:     else                            hx = (xmax-xmin)/M;
     46:     VecGetArray(xcoor,&coors);
     47:     for (i=0; i<isize; i++) {
     48:       coors[i] = xmin + hx*(i+istart);
     49:     }
     50:     VecRestoreArray(xcoor,&coors);
     51:   } else if (dim == 2) {
     52:     if (ymax <= ymin) SETERRQ2(PETSC_ERR_ARG_INCOMP,"ymax must be larger than ymin %g %g",ymin,ymax);
     53:     VecCreateMPI(PETSC_COMM_WORLD,2*isize*jsize,PETSC_DETERMINE,&xcoor);
     54:     VecSetBlockSize(xcoor,2);
     55:     if (DAXPeriodic(periodic)) hx = (xmax-xmin)/(M);
     56:     else                       hx = (xmax-xmin)/(M-1);
     57:     if (DAYPeriodic(periodic)) hy = (ymax-ymin)/(N);
     58:     else                       hy = (ymax-ymin)/(N-1);
     59:     VecGetArray(xcoor,&coors);
     60:     cnt  = 0;
     61:     for (j=0; j<jsize; j++) {
     62:       for (i=0; i<isize; i++) {
     63:         coors[cnt++] = xmin + hx*(i+istart);
     64:         coors[cnt++] = ymin + hy*(j+jstart);
     65:       }
     66:     }
     67:     VecRestoreArray(xcoor,&coors);
     68:   } else if (dim == 3) {
     69:     if (ymax <= ymin) SETERRQ2(PETSC_ERR_ARG_INCOMP,"ymax must be larger than ymin %g %g",ymin,ymax);
     70:     if (zmax <= zmin) SETERRQ2(PETSC_ERR_ARG_INCOMP,"zmax must be larger than zmin %g %g",zmin,zmax);
     71:     VecCreateMPI(PETSC_COMM_WORLD,3*isize*jsize*ksize,PETSC_DETERMINE,&xcoor);
     72:     VecSetBlockSize(xcoor,3);
     73:     if (DAXPeriodic(periodic)) hx = (xmax-xmin)/(M);
     74:     else                       hx = (xmax-xmin)/(M-1);
     75:     if (DAYPeriodic(periodic)) hy = (ymax-ymin)/(N);
     76:     else                       hy = (ymax-ymin)/(N-1);
     77:     if (DAZPeriodic(periodic)) hz_ = (zmax-zmin)/(P);
     78:     else                       hz_ = (zmax-zmin)/(P-1);
     79:     VecGetArray(xcoor,&coors);
     80:     cnt  = 0;
     81:     for (k=0; k<ksize; k++) {
     82:       for (j=0; j<jsize; j++) {
     83:         for (i=0; i<isize; i++) {
     84:           coors[cnt++] = xmin + hx*(i+istart);
     85:           coors[cnt++] = ymin + hy*(j+jstart);
     86:           coors[cnt++] = zmin + hz_*(k+kstart);
     87:         }
     88:       }
     89:     }
     90:     VecRestoreArray(xcoor,&coors);
     91:   } else {
     92:     SETERRQ1(PETSC_ERR_SUP,"Cannot create uniform coordinates for this dimension %D\n",dim);
     93:   }
     94:   DASetCoordinates(da,xcoor);
     95:   PetscLogObjectParent(da,xcoor);

     97:   return(0);
     98: }

    102: PetscErrorCode VecView_MPI_Draw_DA1d(Vec xin,PetscViewer v)
    103: {
    104:   DA             da;
    106:   PetscMPIInt    rank,size,tag1,tag2;
    107:   PetscInt       i,n,N,step,istart,isize,j;
    108:   MPI_Status     status;
    109:   PetscReal      coors[4],ymin,ymax,min,max,xmin,xmax,tmp,xgtmp;
    110:   PetscScalar    *array,*xg;
    111:   PetscDraw      draw;
    112:   PetscTruth     isnull,showpoints;
    113:   MPI_Comm       comm;
    114:   PetscDrawAxis  axis;
    115:   Vec            xcoor;
    116:   DAPeriodicType periodic;

    119:   PetscViewerDrawGetDraw(v,0,&draw);
    120:   PetscDrawIsNull(draw,&isnull); if (isnull) return(0);

    122:   PetscObjectQuery((PetscObject)xin,"DA",(PetscObject*)&da);
    123:   if (!da) SETERRQ(PETSC_ERR_ARG_WRONG,"Vector not generated from a DA");

    125:   PetscOptionsHasName(PETSC_NULL,"-draw_vec_mark_points",&showpoints);

    127:   DAGetInfo(da,0,&N,0,0,0,0,0,&step,0,&periodic,0);
    128:   DAGetCorners(da,&istart,0,0,&isize,0,0);
    129:   VecGetArray(xin,&array);
    130:   VecGetLocalSize(xin,&n);
    131:   n    = n/step;

    133:   /* get coordinates of nodes */
    134:   DAGetCoordinates(da,&xcoor);
    135:   if (!xcoor) {
    136:     DASetUniformCoordinates(da,0.0,1.0,0.0,0.0,0.0,0.0);
    137:     DAGetCoordinates(da,&xcoor);
    138:   }
    139:   VecGetArray(xcoor,&xg);

    141:   PetscObjectGetComm((PetscObject)xin,&comm);
    142:   MPI_Comm_size(comm,&size);
    143:   MPI_Comm_rank(comm,&rank);

    145:   /*
    146:       Determine the min and max x coordinate in plot
    147:   */
    148:   if (!rank) {
    149:     xmin = PetscRealPart(xg[0]);
    150:   }
    151:   if (rank == size-1) {
    152:     xmax = PetscRealPart(xg[n-1]);
    153:   }
    154:   MPI_Bcast(&xmin,1,MPIU_REAL,0,comm);
    155:   MPI_Bcast(&xmax,1,MPIU_REAL,size-1,comm);

    157:   for (j=0; j<step; j++) {
    158:     PetscViewerDrawGetDraw(v,j,&draw);
    159:     PetscDrawCheckResizedWindow(draw);

    161:     /*
    162:         Determine the min and max y coordinate in plot
    163:     */
    164:     min = 1.e20; max = -1.e20;
    165:     for (i=0; i<n; i++) {
    166: #if defined(PETSC_USE_COMPLEX)
    167:       if (PetscRealPart(array[j+i*step]) < min) min = PetscRealPart(array[j+i*step]);
    168:       if (PetscRealPart(array[j+i*step]) > max) max = PetscRealPart(array[j+i*step]);
    169: #else
    170:       if (array[j+i*step] < min) min = array[j+i*step];
    171:       if (array[j+i*step] > max) max = array[j+i*step];
    172: #endif
    173:     }
    174:     if (min + 1.e-10 > max) {
    175:       min -= 1.e-5;
    176:       max += 1.e-5;
    177:     }
    178:     MPI_Reduce(&min,&ymin,1,MPIU_REAL,MPI_MIN,0,comm);
    179:     MPI_Reduce(&max,&ymax,1,MPIU_REAL,MPI_MAX,0,comm);

    181:     PetscDrawSynchronizedClear(draw);
    182:     PetscViewerDrawGetDrawAxis(v,j,&axis);
    183:     PetscLogObjectParent(draw,axis);
    184:     if (!rank) {
    185:       char *title;

    187:       PetscDrawAxisSetLimits(axis,xmin,xmax,ymin,ymax);
    188:       PetscDrawAxisDraw(axis);
    189:       PetscDrawGetCoordinates(draw,coors,coors+1,coors+2,coors+3);
    190:       DAGetFieldName(da,j,&title);
    191:       if (title) {PetscDrawSetTitle(draw,title);}
    192:     }
    193:     MPI_Bcast(coors,4,MPIU_REAL,0,comm);
    194:     if (rank) {
    195:       PetscDrawSetCoordinates(draw,coors[0],coors[1],coors[2],coors[3]);
    196:     }

    198:     /* draw local part of vector */
    199:     PetscObjectGetNewTag((PetscObject)xin,&tag1);
    200:     PetscObjectGetNewTag((PetscObject)xin,&tag2);
    201:     if (rank < size-1) { /*send value to right */
    202:       MPI_Send(&array[j+(n-1)*step],1,MPIU_REAL,rank+1,tag1,comm);
    203:       MPI_Send(&xg[n-1],1,MPIU_REAL,rank+1,tag1,comm);
    204:     }
    205:     if (!rank && periodic) { /* first processor sends first value to last */
    206:       MPI_Send(&array[j],1,MPIU_REAL,size-1,tag2,comm);
    207:     }

    209:     for (i=1; i<n; i++) {
    210: #if !defined(PETSC_USE_COMPLEX)
    211:       PetscDrawLine(draw,xg[i-1],array[j+step*(i-1)],xg[i],array[j+step*i],
    212:                       PETSC_DRAW_RED);
    213: #else
    214:       PetscDrawLine(draw,PetscRealPart(xg[i-1]),PetscRealPart(array[j+step*(i-1)]),
    215:                       PetscRealPart(xg[i]),PetscRealPart(array[j+step*i]),PETSC_DRAW_RED);
    216: #endif
    217:       if (showpoints) {
    218:         PetscDrawPoint(draw,PetscRealPart(xg[i-1]),PetscRealPart(array[j+step*(i-1)]),PETSC_DRAW_BLACK);
    219:       }
    220:     }
    221:     if (rank) { /* receive value from left */
    222:       MPI_Recv(&tmp,1,MPIU_REAL,rank-1,tag1,comm,&status);
    223:       MPI_Recv(&xgtmp,1,MPIU_REAL,rank-1,tag1,comm,&status);
    224: #if !defined(PETSC_USE_COMPLEX)
    225:       PetscDrawLine(draw,xgtmp,tmp,xg[0],array[j],PETSC_DRAW_RED);
    226: #else
    227:       PetscDrawLine(draw,xgtmp,tmp,PetscRealPart(xg[0]),PetscRealPart(array[j]),
    228:                       PETSC_DRAW_RED);
    229: #endif
    230:       if (showpoints) {
    231:         PetscDrawPoint(draw,xgtmp,tmp,PETSC_DRAW_BLACK);
    232:       }
    233:     }
    234:     if (rank == size-1 && periodic) {
    235:       MPI_Recv(&tmp,1,MPIU_REAL,0,tag2,comm,&status);
    236: #if !defined(PETSC_USE_COMPLEX)
    237:       PetscDrawLine(draw,xg[n-2],array[j+step*(n-1)],xg[n-1],tmp,PETSC_DRAW_RED);
    238: #else
    239:       PetscDrawLine(draw,PetscRealPart(xg[n-2]),PetscRealPart(array[j+step*(n-1)]),
    240:                       PetscRealPart(xg[n-1]),tmp,PETSC_DRAW_RED);
    241: #endif
    242:       if (showpoints) {
    243:         PetscDrawPoint(draw,PetscRealPart(xg[n-2]),PetscRealPart(array[j+step*(n-1)]),PETSC_DRAW_BLACK);
    244:       }
    245:     }
    246:     PetscDrawSynchronizedFlush(draw);
    247:     PetscDrawPause(draw);
    248:   }
    249:   VecRestoreArray(xcoor,&xg);
    250:   VecRestoreArray(xin,&array);
    251:   return(0);
    252: }