Actual source code: gr1.c

      2: /*
      3:    Plots vectors obtained with DMDACreate1d()
      4: */

      6: #include <petsc/private/dmdaimpl.h>

      8: /*@
      9:     DMDASetUniformCoordinates - Sets a DMDA coordinates to be a uniform grid

     11:   Collective on da

     13:   Input Parameters:
     14: +  da - the distributed array object
     15: .  xmin,xmax - extremes in the x direction
     16: .  ymin,ymax - extremes in the y direction (value ignored for 1 dimensional problems)
     17: -  zmin,zmax - extremes in the z direction (value ignored for 1 or 2 dimensional problems)

     19:   Level: beginner

     21: .seealso: DMSetCoordinates(), DMGetCoordinates(), DMDACreate1d(), DMDACreate2d(), DMDACreate3d(), DMStagSetUniformCoordinates()

     23: @*/
     24: PetscErrorCode  DMDASetUniformCoordinates(DM da,PetscReal xmin,PetscReal xmax,PetscReal ymin,PetscReal ymax,PetscReal zmin,PetscReal zmax)
     25: {
     26:   MPI_Comm         comm;
     27:   DM               cda;
     28:   DM_DA            *dd = (DM_DA*)da->data;
     29:   DMBoundaryType   bx,by,bz;
     30:   Vec              xcoor;
     31:   PetscScalar      *coors;
     32:   PetscReal        hx,hy,hz_;
     33:   PetscInt         i,j,k,M,N,P,istart,isize,jstart,jsize,kstart,ksize,dim,cnt;
     34:   PetscErrorCode   ierr;

     38:   if (!dd->gtol) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_WRONGSTATE,"Cannot set coordinates until after DMDA has been setup");
     39:   DMDAGetInfo(da,&dim,&M,&N,&P,NULL,NULL,NULL,NULL,NULL,&bx,&by,&bz,NULL);
     40:   if (xmax < xmin) SETERRQ2(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_INCOMP,"xmax must be larger than xmin %g %g",(double)xmin,(double)xmax);
     41:   if ((dim > 1) && (ymax < ymin)) SETERRQ2(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_INCOMP,"ymax must be larger than ymin %g %g",(double)ymin,(double)ymax);
     42:   if ((dim > 2) && (zmax < zmin)) SETERRQ2(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_INCOMP,"zmax must be larger than zmin %g %g",(double)zmin,(double)zmax);
     43:   PetscObjectGetComm((PetscObject)da,&comm);
     44:   DMDAGetCorners(da,&istart,&jstart,&kstart,&isize,&jsize,&ksize);
     45:   DMGetCoordinateDM(da, &cda);
     46:   DMCreateGlobalVector(cda, &xcoor);
     47:   if (dim == 1) {
     48:     if (bx == DM_BOUNDARY_PERIODIC) hx = (xmax-xmin)/M;
     49:     else hx = (xmax-xmin)/(M-1);
     50:     VecGetArray(xcoor,&coors);
     51:     for (i=0; i<isize; i++) {
     52:       coors[i] = xmin + hx*(i+istart);
     53:     }
     54:     VecRestoreArray(xcoor,&coors);
     55:   } else if (dim == 2) {
     56:     if (bx == DM_BOUNDARY_PERIODIC) hx = (xmax-xmin)/(M);
     57:     else hx = (xmax-xmin)/(M-1);
     58:     if (by == DM_BOUNDARY_PERIODIC) hy = (ymax-ymin)/(N);
     59:     else hy = (ymax-ymin)/(N-1);
     60:     VecGetArray(xcoor,&coors);
     61:     cnt  = 0;
     62:     for (j=0; j<jsize; j++) {
     63:       for (i=0; i<isize; i++) {
     64:         coors[cnt++] = xmin + hx*(i+istart);
     65:         coors[cnt++] = ymin + hy*(j+jstart);
     66:       }
     67:     }
     68:     VecRestoreArray(xcoor,&coors);
     69:   } else if (dim == 3) {
     70:     if (bx == DM_BOUNDARY_PERIODIC) hx = (xmax-xmin)/(M);
     71:     else hx = (xmax-xmin)/(M-1);
     72:     if (by == DM_BOUNDARY_PERIODIC) hy = (ymax-ymin)/(N);
     73:     else hy = (ymax-ymin)/(N-1);
     74:     if (bz == DM_BOUNDARY_PERIODIC) hz_ = (zmax-zmin)/(P);
     75:     else hz_ = (zmax-zmin)/(P-1);
     76:     VecGetArray(xcoor,&coors);
     77:     cnt  = 0;
     78:     for (k=0; k<ksize; k++) {
     79:       for (j=0; j<jsize; j++) {
     80:         for (i=0; i<isize; i++) {
     81:           coors[cnt++] = xmin + hx*(i+istart);
     82:           coors[cnt++] = ymin + hy*(j+jstart);
     83:           coors[cnt++] = zmin + hz_*(k+kstart);
     84:         }
     85:       }
     86:     }
     87:     VecRestoreArray(xcoor,&coors);
     88:   } else SETERRQ1(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"Cannot create uniform coordinates for this dimension %D\n",dim);
     89:   DMSetCoordinates(da,xcoor);
     90:   PetscLogObjectParent((PetscObject)da,(PetscObject)xcoor);
     91:   VecDestroy(&xcoor);
     92:   return(0);
     93: }

     95: /*
     96:     Allows a user to select a subset of the fields to be drawn by VecView() when the vector comes from a DMDA
     97: */
     98: PetscErrorCode DMDASelectFields(DM da,PetscInt *outfields,PetscInt **fields)
     99: {
    101:   PetscInt       step,ndisplayfields,*displayfields,k,j;
    102:   PetscBool      flg;

    105:   DMDAGetInfo(da,NULL,NULL,NULL,NULL,NULL,NULL,NULL,&step,NULL,NULL,NULL,NULL,NULL);
    106:   PetscMalloc1(step,&displayfields);
    107:   for (k=0; k<step; k++) displayfields[k] = k;
    108:   ndisplayfields = step;
    109:   PetscOptionsGetIntArray(NULL,NULL,"-draw_fields",displayfields,&ndisplayfields,&flg);
    110:   if (!ndisplayfields) ndisplayfields = step;
    111:   if (!flg) {
    112:     char       **fields;
    113:     const char *fieldname;
    114:     PetscInt   nfields = step;
    115:     PetscMalloc1(step,&fields);
    116:     PetscOptionsGetStringArray(NULL,NULL,"-draw_fields_by_name",fields,&nfields,&flg);
    117:     if (flg) {
    118:       ndisplayfields = 0;
    119:       for (k=0; k<nfields;k++) {
    120:         for (j=0; j<step; j++) {
    121:           DMDAGetFieldName(da,j,&fieldname);
    122:           PetscStrcmp(fieldname,fields[k],&flg);
    123:           if (flg) {
    124:             goto found;
    125:           }
    126:         }
    127:         SETERRQ1(PetscObjectComm((PetscObject)da),PETSC_ERR_USER,"Unknown fieldname %s",fields[k]);
    128: found:  displayfields[ndisplayfields++] = j;
    129:       }
    130:     }
    131:     for (k=0; k<nfields; k++) {
    132:       PetscFree(fields[k]);
    133:     }
    134:     PetscFree(fields);
    135:   }
    136:   *fields    = displayfields;
    137:   *outfields = ndisplayfields;
    138:   return(0);
    139: }

    141: #include <petscdraw.h>

    143: PetscErrorCode VecView_MPI_Draw_DA1d(Vec xin,PetscViewer v)
    144: {
    145:   DM                da;
    146:   PetscErrorCode    ierr;
    147:   PetscMPIInt       rank,size,tag;
    148:   PetscInt          i,n,N,dof,istart,isize,j,nbounds;
    149:   MPI_Status        status;
    150:   PetscReal         min,max,xmin = 0.0,xmax = 0.0,tmp = 0.0,xgtmp = 0.0;
    151:   const PetscScalar *array,*xg;
    152:   PetscDraw         draw;
    153:   PetscBool         isnull,useports = PETSC_FALSE,showmarkers = PETSC_FALSE;
    154:   MPI_Comm          comm;
    155:   PetscDrawAxis     axis;
    156:   Vec               xcoor;
    157:   DMBoundaryType    bx;
    158:   const char        *tlabel = NULL,*xlabel = NULL;
    159:   const PetscReal   *bounds;
    160:   PetscInt          *displayfields;
    161:   PetscInt          k,ndisplayfields;
    162:   PetscBool         hold;
    163:   PetscDrawViewPorts *ports = NULL;
    164:   PetscViewerFormat  format;

    167:   PetscViewerDrawGetDraw(v,0,&draw);
    168:   PetscDrawIsNull(draw,&isnull);
    169:   if (isnull) return(0);
    170:   PetscViewerDrawGetBounds(v,&nbounds,&bounds);

    172:   VecGetDM(xin,&da);
    173:   if (!da) SETERRQ(PetscObjectComm((PetscObject)xin),PETSC_ERR_ARG_WRONG,"Vector not generated from a DMDA");
    174:   PetscObjectGetComm((PetscObject)xin,&comm);
    175:   MPI_Comm_size(comm,&size);
    176:   MPI_Comm_rank(comm,&rank);

    178:   PetscOptionsGetBool(NULL,NULL,"-draw_vec_use_markers",&showmarkers,NULL);

    180:   DMDAGetInfo(da,NULL,&N,NULL,NULL,NULL,NULL,NULL,&dof,NULL,&bx,NULL,NULL,NULL);
    181:   DMDAGetCorners(da,&istart,NULL,NULL,&isize,NULL,NULL);
    182:   VecGetArrayRead(xin,&array);
    183:   VecGetLocalSize(xin,&n);
    184:   n    = n/dof;

    186:   /* Get coordinates of nodes */
    187:   DMGetCoordinates(da,&xcoor);
    188:   if (!xcoor) {
    189:     DMDASetUniformCoordinates(da,0.0,1.0,0.0,0.0,0.0,0.0);
    190:     DMGetCoordinates(da,&xcoor);
    191:   }
    192:   VecGetArrayRead(xcoor,&xg);
    193:   DMDAGetCoordinateName(da,0,&xlabel);

    195:   /* Determine the min and max coordinate in plot */
    196:   if (!rank) xmin = PetscRealPart(xg[0]);
    197:   if (rank == size-1) xmax = PetscRealPart(xg[n-1]);
    198:   MPI_Bcast(&xmin,1,MPIU_REAL,0,comm);
    199:   MPI_Bcast(&xmax,1,MPIU_REAL,size-1,comm);

    201:   DMDASelectFields(da,&ndisplayfields,&displayfields);
    202:   PetscViewerGetFormat(v,&format);
    203:   PetscOptionsGetBool(NULL,NULL,"-draw_ports",&useports,NULL);
    204:   if (format == PETSC_VIEWER_DRAW_PORTS) useports = PETSC_TRUE;
    205:   if (useports) {
    206:     PetscViewerDrawGetDraw(v,0,&draw);
    207:     PetscViewerDrawGetDrawAxis(v,0,&axis);
    208:     PetscDrawCheckResizedWindow(draw);
    209:     PetscDrawClear(draw);
    210:     PetscDrawViewPortsCreate(draw,ndisplayfields,&ports);
    211:   }

    213:   /* Loop over each field; drawing each in a different window */
    214:   for (k=0; k<ndisplayfields; k++) {
    215:     j = displayfields[k];

    217:     /* determine the min and max value in plot */
    218:     VecStrideMin(xin,j,NULL,&min);
    219:     VecStrideMax(xin,j,NULL,&max);
    220:     if (j < nbounds) {
    221:       min = PetscMin(min,bounds[2*j]);
    222:       max = PetscMax(max,bounds[2*j+1]);
    223:     }
    224:     if (min == max) {
    225:       min -= 1.e-5;
    226:       max += 1.e-5;
    227:     }

    229:     if (useports) {
    230:       PetscDrawViewPortsSet(ports,k);
    231:       DMDAGetFieldName(da,j,&tlabel);
    232:     } else {
    233:       const char *title;
    234:       PetscViewerDrawGetHold(v,&hold);
    235:       PetscViewerDrawGetDraw(v,k,&draw);
    236:       PetscViewerDrawGetDrawAxis(v,k,&axis);
    237:       DMDAGetFieldName(da,j,&title);
    238:       if (title) {PetscDrawSetTitle(draw,title);}
    239:       PetscDrawCheckResizedWindow(draw);
    240:       if (!hold) {PetscDrawClear(draw);}
    241:     }
    242:     PetscDrawAxisSetLabels(axis,tlabel,xlabel,NULL);
    243:     PetscDrawAxisSetLimits(axis,xmin,xmax,min,max);
    244:     PetscDrawAxisDraw(axis);

    246:     /* draw local part of vector */
    247:     PetscObjectGetNewTag((PetscObject)xin,&tag);
    248:     if (rank < size-1) { /*send value to right */
    249:       MPI_Send((void*)&xg[n-1],1,MPIU_REAL,rank+1,tag,comm);
    250:       MPI_Send((void*)&array[j+(n-1)*dof],1,MPIU_REAL,rank+1,tag,comm);
    251:     }
    252:     if (rank) { /* receive value from left */
    253:       MPI_Recv(&xgtmp,1,MPIU_REAL,rank-1,tag,comm,&status);
    254:       MPI_Recv(&tmp,1,MPIU_REAL,rank-1,tag,comm,&status);
    255:     }
    256:     PetscDrawCollectiveBegin(draw);
    257:     if (rank) {
    258:       PetscDrawLine(draw,xgtmp,tmp,PetscRealPart(xg[0]),PetscRealPart(array[j]),PETSC_DRAW_RED);
    259:       if (showmarkers) {PetscDrawPoint(draw,xgtmp,tmp,PETSC_DRAW_BLACK);}
    260:     }
    261:     for (i=1; i<n; i++) {
    262:       PetscDrawLine(draw,PetscRealPart(xg[i-1]),PetscRealPart(array[j+dof*(i-1)]),PetscRealPart(xg[i]),PetscRealPart(array[j+dof*i]),PETSC_DRAW_RED);
    263:       if (showmarkers) {PetscDrawMarker(draw,PetscRealPart(xg[i-1]),PetscRealPart(array[j+dof*(i-1)]),PETSC_DRAW_BLACK);}
    264:     }
    265:     if (rank == size-1) {
    266:       if (showmarkers) {PetscDrawMarker(draw,PetscRealPart(xg[n-1]),PetscRealPart(array[j+dof*(n-1)]),PETSC_DRAW_BLACK);}
    267:     }
    268:     PetscDrawCollectiveEnd(draw);
    269:     PetscDrawFlush(draw);
    270:     PetscDrawPause(draw);
    271:     if (!useports) {PetscDrawSave(draw);}
    272:   }
    273:   if (useports) {
    274:     PetscViewerDrawGetDraw(v,0,&draw);
    275:     PetscDrawSave(draw);
    276:   }

    278:   PetscDrawViewPortsDestroy(ports);
    279:   PetscFree(displayfields);
    280:   VecRestoreArrayRead(xcoor,&xg);
    281:   VecRestoreArrayRead(xin,&array);
    282:   return(0);
    283: }