1 #define PETSCDM_DLL 2 3 /* 4 Plots vectors obtained with DMDACreate1d() 5 */ 6 7 #include "petscdm.h" /*I "petscdm.h" I*/ 8 9 #undef __FUNCT__ 10 #define __FUNCT__ "DMDASetUniformCoordinates" 11 /*@ 12 DMDASetUniformCoordinates - Sets a DMDA coordinates to be a uniform grid 13 14 Collective on DMDA 15 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) 21 22 Level: beginner 23 24 .seealso: DMDASetCoordinates(), DMDAGetCoordinates(), DMDACreate1d(), DMDACreate2d(), DMDACreate3d() 25 26 @*/ 27 PetscErrorCode PETSCDM_DLLEXPORT DMDASetUniformCoordinates(DM da,PetscReal xmin,PetscReal xmax,PetscReal ymin,PetscReal ymax,PetscReal zmin,PetscReal zmax) 28 { 29 MPI_Comm comm; 30 DM cda; 31 DMDAPeriodicType periodic; 32 Vec xcoor; 33 PetscScalar *coors; 34 PetscReal hx,hy,hz_; 35 PetscInt i,j,k,M,N,P,istart,isize,jstart,jsize,kstart,ksize,dim,cnt; 36 PetscErrorCode ierr; 37 38 PetscFunctionBegin; 39 if (xmax <= xmin) SETERRQ2(((PetscObject)da)->comm,PETSC_ERR_ARG_INCOMP,"xmax must be larger than xmin %G %G",xmin,xmax); 40 41 ierr = PetscObjectGetComm((PetscObject)da,&comm);CHKERRQ(ierr); 42 ierr = DMDAGetInfo(da,&dim,&M,&N,&P,0,0,0,0,0,&periodic,0);CHKERRQ(ierr); 43 ierr = DMDAGetCorners(da,&istart,&jstart,&kstart,&isize,&jsize,&ksize);CHKERRQ(ierr); 44 ierr = DMDAGetCoordinateDA(da, &cda);CHKERRQ(ierr); 45 ierr = DMCreateGlobalVector(cda, &xcoor);CHKERRQ(ierr); 46 if (dim == 1) { 47 if (periodic == DMDA_NONPERIODIC) hx = (xmax-xmin)/(M-1); 48 else hx = (xmax-xmin)/M; 49 ierr = VecGetArray(xcoor,&coors);CHKERRQ(ierr); 50 for (i=0; i<isize; i++) { 51 coors[i] = xmin + hx*(i+istart); 52 } 53 ierr = VecRestoreArray(xcoor,&coors);CHKERRQ(ierr); 54 } else if (dim == 2) { 55 if (ymax <= ymin) SETERRQ2(((PetscObject)da)->comm,PETSC_ERR_ARG_INCOMP,"ymax must be larger than ymin %G %G",ymin,ymax); 56 if (DMDAXPeriodic(periodic)) hx = (xmax-xmin)/(M); 57 else hx = (xmax-xmin)/(M-1); 58 if (DMDAYPeriodic(periodic)) hy = (ymax-ymin)/(N); 59 else hy = (ymax-ymin)/(N-1); 60 ierr = VecGetArray(xcoor,&coors);CHKERRQ(ierr); 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 ierr = VecRestoreArray(xcoor,&coors);CHKERRQ(ierr); 69 } else if (dim == 3) { 70 if (ymax <= ymin) SETERRQ2(((PetscObject)da)->comm,PETSC_ERR_ARG_INCOMP,"ymax must be larger than ymin %G %G",ymin,ymax); 71 if (zmax <= zmin) SETERRQ2(((PetscObject)da)->comm,PETSC_ERR_ARG_INCOMP,"zmax must be larger than zmin %G %G",zmin,zmax); 72 if (DMDAXPeriodic(periodic)) hx = (xmax-xmin)/(M); 73 else hx = (xmax-xmin)/(M-1); 74 if (DMDAYPeriodic(periodic)) hy = (ymax-ymin)/(N); 75 else hy = (ymax-ymin)/(N-1); 76 if (DMDAZPeriodic(periodic)) hz_ = (zmax-zmin)/(P); 77 else hz_ = (zmax-zmin)/(P-1); 78 ierr = VecGetArray(xcoor,&coors);CHKERRQ(ierr); 79 cnt = 0; 80 for (k=0; k<ksize; k++) { 81 for (j=0; j<jsize; j++) { 82 for (i=0; i<isize; i++) { 83 coors[cnt++] = xmin + hx*(i+istart); 84 coors[cnt++] = ymin + hy*(j+jstart); 85 coors[cnt++] = zmin + hz_*(k+kstart); 86 } 87 } 88 } 89 ierr = VecRestoreArray(xcoor,&coors);CHKERRQ(ierr); 90 } else { 91 SETERRQ1(((PetscObject)da)->comm,PETSC_ERR_SUP,"Cannot create uniform coordinates for this dimension %D\n",dim); 92 } 93 ierr = DMDASetCoordinates(da,xcoor);CHKERRQ(ierr); 94 ierr = PetscLogObjectParent(da,xcoor);CHKERRQ(ierr); 95 ierr = VecDestroy(xcoor);CHKERRQ(ierr); 96 PetscFunctionReturn(0); 97 } 98 99 #undef __FUNCT__ 100 #define __FUNCT__ "VecView_MPI_Draw_DA1d" 101 PetscErrorCode VecView_MPI_Draw_DA1d(Vec xin,PetscViewer v) 102 { 103 DM da; 104 PetscErrorCode ierr; 105 PetscMPIInt rank,size,tag1,tag2; 106 PetscInt i,n,N,step,istart,isize,j; 107 MPI_Status status; 108 PetscReal coors[4],ymin,ymax,min,max,xmin,xmax,tmp,xgtmp; 109 const PetscScalar *array,*xg; 110 PetscDraw draw; 111 PetscBool isnull,showpoints = PETSC_FALSE; 112 MPI_Comm comm; 113 PetscDrawAxis axis; 114 Vec xcoor; 115 DMDAPeriodicType periodic; 116 117 PetscFunctionBegin; 118 ierr = PetscViewerDrawGetDraw(v,0,&draw);CHKERRQ(ierr); 119 ierr = PetscDrawIsNull(draw,&isnull);CHKERRQ(ierr); if (isnull) PetscFunctionReturn(0); 120 121 ierr = PetscObjectQuery((PetscObject)xin,"DMDA",(PetscObject*)&da);CHKERRQ(ierr); 122 if (!da) SETERRQ(((PetscObject)xin)->comm,PETSC_ERR_ARG_WRONG,"Vector not generated from a DMDA"); 123 124 ierr = PetscOptionsGetBool(PETSC_NULL,"-draw_vec_mark_points",&showpoints,PETSC_NULL);CHKERRQ(ierr); 125 126 ierr = DMDAGetInfo(da,0,&N,0,0,0,0,0,&step,0,&periodic,0);CHKERRQ(ierr); 127 ierr = DMDAGetCorners(da,&istart,0,0,&isize,0,0);CHKERRQ(ierr); 128 ierr = VecGetArrayRead(xin,&array);CHKERRQ(ierr); 129 ierr = VecGetLocalSize(xin,&n);CHKERRQ(ierr); 130 n = n/step; 131 132 /* get coordinates of nodes */ 133 ierr = DMDAGetCoordinates(da,&xcoor);CHKERRQ(ierr); 134 if (!xcoor) { 135 ierr = DMDASetUniformCoordinates(da,0.0,1.0,0.0,0.0,0.0,0.0);CHKERRQ(ierr); 136 ierr = DMDAGetCoordinates(da,&xcoor);CHKERRQ(ierr); 137 } 138 ierr = VecGetArrayRead(xcoor,&xg);CHKERRQ(ierr); 139 140 ierr = PetscObjectGetComm((PetscObject)xin,&comm);CHKERRQ(ierr); 141 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 142 ierr = MPI_Comm_rank(comm,&rank);CHKERRQ(ierr); 143 144 /* 145 Determine the min and max x coordinate in plot 146 */ 147 if (!rank) { 148 xmin = PetscRealPart(xg[0]); 149 } 150 if (rank == size-1) { 151 xmax = PetscRealPart(xg[n-1]); 152 } 153 ierr = MPI_Bcast(&xmin,1,MPIU_REAL,0,comm);CHKERRQ(ierr); 154 ierr = MPI_Bcast(&xmax,1,MPIU_REAL,size-1,comm);CHKERRQ(ierr); 155 156 for (j=0; j<step; j++) { 157 ierr = PetscViewerDrawGetDraw(v,j,&draw);CHKERRQ(ierr); 158 ierr = PetscDrawCheckResizedWindow(draw);CHKERRQ(ierr); 159 160 /* 161 Determine the min and max y coordinate in plot 162 */ 163 min = 1.e20; max = -1.e20; 164 for (i=0; i<n; i++) { 165 #if defined(PETSC_USE_COMPLEX) 166 if (PetscRealPart(array[j+i*step]) < min) min = PetscRealPart(array[j+i*step]); 167 if (PetscRealPart(array[j+i*step]) > max) max = PetscRealPart(array[j+i*step]); 168 #else 169 if (array[j+i*step] < min) min = array[j+i*step]; 170 if (array[j+i*step] > max) max = array[j+i*step]; 171 #endif 172 } 173 if (min + 1.e-10 > max) { 174 min -= 1.e-5; 175 max += 1.e-5; 176 } 177 ierr = MPI_Reduce(&min,&ymin,1,MPIU_REAL,MPI_MIN,0,comm);CHKERRQ(ierr); 178 ierr = MPI_Reduce(&max,&ymax,1,MPIU_REAL,MPI_MAX,0,comm);CHKERRQ(ierr); 179 180 ierr = PetscDrawSynchronizedClear(draw);CHKERRQ(ierr); 181 ierr = PetscViewerDrawGetDrawAxis(v,j,&axis);CHKERRQ(ierr); 182 ierr = PetscLogObjectParent(draw,axis);CHKERRQ(ierr); 183 if (!rank) { 184 const char *title; 185 186 ierr = PetscDrawAxisSetLimits(axis,xmin,xmax,ymin,ymax);CHKERRQ(ierr); 187 ierr = PetscDrawAxisDraw(axis);CHKERRQ(ierr); 188 ierr = PetscDrawGetCoordinates(draw,coors,coors+1,coors+2,coors+3);CHKERRQ(ierr); 189 ierr = DMDAGetFieldName(da,j,&title);CHKERRQ(ierr); 190 if (title) {ierr = PetscDrawSetTitle(draw,title);CHKERRQ(ierr);} 191 } 192 ierr = MPI_Bcast(coors,4,MPIU_REAL,0,comm);CHKERRQ(ierr); 193 if (rank) { 194 ierr = PetscDrawSetCoordinates(draw,coors[0],coors[1],coors[2],coors[3]);CHKERRQ(ierr); 195 } 196 197 /* draw local part of vector */ 198 PetscObjectGetNewTag((PetscObject)xin,&tag1);CHKERRQ(ierr); 199 PetscObjectGetNewTag((PetscObject)xin,&tag2);CHKERRQ(ierr); 200 if (rank < size-1) { /*send value to right */ 201 ierr = MPI_Send((void*)&array[j+(n-1)*step],1,MPIU_REAL,rank+1,tag1,comm);CHKERRQ(ierr); 202 ierr = MPI_Send((void*)&xg[n-1],1,MPIU_REAL,rank+1,tag1,comm);CHKERRQ(ierr); 203 } 204 if (!rank && periodic && size > 1) { /* first processor sends first value to last */ 205 ierr = MPI_Send((void*)&array[j],1,MPIU_REAL,size-1,tag2,comm);CHKERRQ(ierr); 206 } 207 208 for (i=1; i<n; i++) { 209 #if !defined(PETSC_USE_COMPLEX) 210 ierr = PetscDrawLine(draw,xg[i-1],array[j+step*(i-1)],xg[i],array[j+step*i],PETSC_DRAW_RED);CHKERRQ(ierr); 211 #else 212 ierr = PetscDrawLine(draw,PetscRealPart(xg[i-1]),PetscRealPart(array[j+step*(i-1)]),PetscRealPart(xg[i]),PetscRealPart(array[j+step*i]),PETSC_DRAW_RED);CHKERRQ(ierr); 213 #endif 214 if (showpoints) { 215 ierr = PetscDrawPoint(draw,PetscRealPart(xg[i-1]),PetscRealPart(array[j+step*(i-1)]),PETSC_DRAW_BLACK);CHKERRQ(ierr); 216 } 217 } 218 if (rank) { /* receive value from left */ 219 ierr = MPI_Recv(&tmp,1,MPIU_REAL,rank-1,tag1,comm,&status);CHKERRQ(ierr); 220 ierr = MPI_Recv(&xgtmp,1,MPIU_REAL,rank-1,tag1,comm,&status);CHKERRQ(ierr); 221 #if !defined(PETSC_USE_COMPLEX) 222 ierr = PetscDrawLine(draw,xgtmp,tmp,xg[0],array[j],PETSC_DRAW_RED);CHKERRQ(ierr); 223 #else 224 ierr = PetscDrawLine(draw,xgtmp,tmp,PetscRealPart(xg[0]),PetscRealPart(array[j]),PETSC_DRAW_RED);CHKERRQ(ierr); 225 #endif 226 if (showpoints) { 227 ierr = PetscDrawPoint(draw,xgtmp,tmp,PETSC_DRAW_BLACK);CHKERRQ(ierr); 228 } 229 } 230 if (rank == size-1 && periodic && size > 1) { 231 ierr = MPI_Recv(&tmp,1,MPIU_REAL,0,tag2,comm,&status);CHKERRQ(ierr); 232 #if !defined(PETSC_USE_COMPLEX) 233 ierr = PetscDrawLine(draw,xg[n-2],array[j+step*(n-1)],xg[n-1],tmp,PETSC_DRAW_RED);CHKERRQ(ierr); 234 #else 235 ierr = PetscDrawLine(draw,PetscRealPart(xg[n-2]),PetscRealPart(array[j+step*(n-1)]), 236 PetscRealPart(xg[n-1]),tmp,PETSC_DRAW_RED);CHKERRQ(ierr); 237 #endif 238 if (showpoints) { 239 ierr = PetscDrawPoint(draw,PetscRealPart(xg[n-2]),PetscRealPart(array[j+step*(n-1)]),PETSC_DRAW_BLACK);CHKERRQ(ierr); 240 } 241 } 242 ierr = PetscDrawSynchronizedFlush(draw);CHKERRQ(ierr); 243 ierr = PetscDrawPause(draw);CHKERRQ(ierr); 244 } 245 ierr = VecRestoreArrayRead(xcoor,&xg);CHKERRQ(ierr); 246 ierr = VecRestoreArrayRead(xin,&array);CHKERRQ(ierr); 247 PetscFunctionReturn(0); 248 } 249 250