1 #define PETSCDM_DLL 2 3 /* 4 Code for manipulating distributed regular arrays in parallel. 5 */ 6 7 #include "private/daimpl.h" /*I "petscda.h" I*/ 8 9 #if defined(PETSC_HAVE_MATLAB_ENGINE) 10 #include "mat.h" /* Matlab include file */ 11 12 #undef __FUNCT__ 13 #define __FUNCT__ "DAView_Matlab" 14 PetscErrorCode DAView_Matlab(DA da,PetscViewer viewer) 15 { 16 PetscErrorCode ierr; 17 PetscMPIInt rank; 18 PetscInt dim,m,n,p,dof,swidth; 19 DAStencilType stencil; 20 DAPeriodicType periodic; 21 mxArray *mx; 22 const char *fnames[] = {"dimension","m","n","p","dof","stencil_width","periodicity","stencil_type"}; 23 24 PetscFunctionBegin; 25 ierr = MPI_Comm_rank(((PetscObject)da)->comm,&rank);CHKERRQ(ierr); 26 if (!rank) { 27 ierr = DAGetInfo(da,&dim,&m,&n,&p,0,0,0,&dof,&swidth,&periodic,&stencil);CHKERRQ(ierr); 28 mx = mxCreateStructMatrix(1,1,8,(const char **)fnames); 29 if (!mx) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Unable to generate Matlab struct array to hold DA informations"); 30 mxSetFieldByNumber(mx,0,0,mxCreateDoubleScalar((double)dim)); 31 mxSetFieldByNumber(mx,0,1,mxCreateDoubleScalar((double)m)); 32 mxSetFieldByNumber(mx,0,2,mxCreateDoubleScalar((double)n)); 33 mxSetFieldByNumber(mx,0,3,mxCreateDoubleScalar((double)p)); 34 mxSetFieldByNumber(mx,0,4,mxCreateDoubleScalar((double)dof)); 35 mxSetFieldByNumber(mx,0,5,mxCreateDoubleScalar((double)swidth)); 36 mxSetFieldByNumber(mx,0,6,mxCreateDoubleScalar((double)periodic)); 37 mxSetFieldByNumber(mx,0,7,mxCreateDoubleScalar((double)stencil)); 38 ierr = PetscObjectName((PetscObject)da);CHKERRQ(ierr); 39 ierr = PetscViewerMatlabPutVariable(viewer,((PetscObject)da)->name,mx);CHKERRQ(ierr); 40 } 41 PetscFunctionReturn(0); 42 } 43 #endif 44 45 #undef __FUNCT__ 46 #define __FUNCT__ "DAView_Binary" 47 PetscErrorCode DAView_Binary(DA da,PetscViewer viewer) 48 { 49 PetscErrorCode ierr; 50 PetscMPIInt rank; 51 PetscInt i,dim,m,n,p,dof,swidth,M,N,P; 52 size_t j,len; 53 DAStencilType stencil; 54 DAPeriodicType periodic; 55 MPI_Comm comm; 56 DM_DA *dd = (DM_DA*)da->data; 57 58 PetscFunctionBegin; 59 ierr = PetscObjectGetComm((PetscObject)da,&comm);CHKERRQ(ierr); 60 61 ierr = DAGetInfo(da,&dim,&m,&n,&p,&M,&N,&P,&dof,&swidth,&periodic,&stencil);CHKERRQ(ierr); 62 ierr = MPI_Comm_rank(comm,&rank);CHKERRQ(ierr); 63 if (!rank) { 64 FILE *file; 65 66 ierr = PetscViewerBinaryGetInfoPointer(viewer,&file);CHKERRQ(ierr); 67 if (file) { 68 char fieldname[PETSC_MAX_PATH_LEN]; 69 70 ierr = PetscFPrintf(PETSC_COMM_SELF,file,"-daload_info %D,%D,%D,%D,%D,%D,%D,%D\n",dim,m,n,p,dof,swidth,stencil,periodic);CHKERRQ(ierr); 71 for (i=0; i<dof; i++) { 72 if (dd->fieldname[i]) { 73 ierr = PetscStrncpy(fieldname,dd->fieldname[i],PETSC_MAX_PATH_LEN);CHKERRQ(ierr); 74 ierr = PetscStrlen(fieldname,&len);CHKERRQ(ierr); 75 len = PetscMin(PETSC_MAX_PATH_LEN,len);CHKERRQ(ierr); 76 for (j=0; j<len; j++) { 77 if (fieldname[j] == ' ') fieldname[j] = '_'; 78 } 79 ierr = PetscFPrintf(PETSC_COMM_SELF,file,"-daload_fieldname_%D %s\n",i,fieldname);CHKERRQ(ierr); 80 } 81 } 82 if (dd->coordinates) { /* save the DA's coordinates */ 83 ierr = PetscFPrintf(PETSC_COMM_SELF,file,"-daload_coordinates\n");CHKERRQ(ierr); 84 } 85 } 86 } 87 88 /* save the coordinates if they exist to disk (in the natural ordering) */ 89 if (dd->coordinates) { 90 DA dac; 91 const PetscInt *lx,*ly,*lz; 92 Vec natural; 93 94 /* create the appropriate DA to map to natural ordering */ 95 ierr = DAGetOwnershipRanges(da,&lx,&ly,&lz);CHKERRQ(ierr); 96 if (dim == 1) { 97 ierr = DACreate1d(comm,DA_NONPERIODIC,m,dim,0,lx,&dac);CHKERRQ(ierr); 98 } else if (dim == 2) { 99 ierr = DACreate2d(comm,DA_NONPERIODIC,DA_STENCIL_BOX,m,n,M,N,dim,0,lx,ly,&dac);CHKERRQ(ierr); 100 } else if (dim == 3) { 101 ierr = DACreate3d(comm,DA_NONPERIODIC,DA_STENCIL_BOX,m,n,p,M,N,P,dim,0,lx,ly,lz,&dac);CHKERRQ(ierr); 102 } else { 103 SETERRQ1(comm,PETSC_ERR_ARG_CORRUPT,"Dimension is not 1 2 or 3: %D\n",dim); 104 } 105 ierr = DACreateNaturalVector(dac,&natural);CHKERRQ(ierr); 106 ierr = PetscObjectSetOptionsPrefix((PetscObject)natural,"coor_");CHKERRQ(ierr); 107 ierr = DAGlobalToNaturalBegin(dac,dd->coordinates,INSERT_VALUES,natural);CHKERRQ(ierr); 108 ierr = DAGlobalToNaturalEnd(dac,dd->coordinates,INSERT_VALUES,natural);CHKERRQ(ierr); 109 ierr = VecView(natural,viewer);CHKERRQ(ierr); 110 ierr = VecDestroy(natural);CHKERRQ(ierr); 111 ierr = DADestroy(dac);CHKERRQ(ierr); 112 } 113 114 PetscFunctionReturn(0); 115 } 116 117 #undef __FUNCT__ 118 #define __FUNCT__ "DAView_VTK" 119 PetscErrorCode DAView_VTK(DA da, PetscViewer viewer) 120 { 121 PetscInt dim, dof, M = 0, N = 0, P = 0; 122 PetscErrorCode ierr; 123 DM_DA *dd = (DM_DA*)da->data; 124 125 PetscFunctionBegin; 126 ierr = DAGetInfo(da, &dim, &M, &N, &P, PETSC_NULL, PETSC_NULL, PETSC_NULL, &dof, PETSC_NULL, PETSC_NULL, PETSC_NULL);CHKERRQ(ierr); 127 /* if (dim != 3) {SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP, "VTK output only works for three dimensional DAs.");} */ 128 if (!dd->coordinates) SETERRQ(((PetscObject)da)->comm,PETSC_ERR_SUP, "VTK output requires DA coordinates."); 129 /* Write Header */ 130 ierr = PetscViewerASCIIPrintf(viewer,"# vtk DataFile Version 2.0\n");CHKERRQ(ierr); 131 ierr = PetscViewerASCIIPrintf(viewer,"Structured Mesh Example\n");CHKERRQ(ierr); 132 ierr = PetscViewerASCIIPrintf(viewer,"ASCII\n");CHKERRQ(ierr); 133 ierr = PetscViewerASCIIPrintf(viewer,"DATASET STRUCTURED_GRID\n");CHKERRQ(ierr); 134 ierr = PetscViewerASCIIPrintf(viewer,"DIMENSIONS %d %d %d\n", M, N, P);CHKERRQ(ierr); 135 ierr = PetscViewerASCIIPrintf(viewer,"POINTS %d double\n", M*N*P);CHKERRQ(ierr); 136 if (dd->coordinates) { 137 DA dac; 138 Vec natural; 139 140 ierr = DAGetCoordinateDA(da, &dac);CHKERRQ(ierr); 141 ierr = DACreateNaturalVector(dac, &natural);CHKERRQ(ierr); 142 ierr = PetscObjectSetOptionsPrefix((PetscObject) natural, "coor_");CHKERRQ(ierr); 143 ierr = DAGlobalToNaturalBegin(dac, dd->coordinates, INSERT_VALUES, natural);CHKERRQ(ierr); 144 ierr = DAGlobalToNaturalEnd(dac, dd->coordinates, INSERT_VALUES, natural);CHKERRQ(ierr); 145 ierr = PetscViewerPushFormat(viewer, PETSC_VIEWER_ASCII_VTK_COORDS);CHKERRQ(ierr); 146 ierr = VecView(natural, viewer);CHKERRQ(ierr); 147 ierr = PetscViewerPopFormat(viewer);CHKERRQ(ierr); 148 ierr = VecDestroy(natural);CHKERRQ(ierr); 149 } 150 PetscFunctionReturn(0); 151 } 152 153 #undef __FUNCT__ 154 #define __FUNCT__ "DAView" 155 /*@C 156 DAView - Visualizes a distributed array object. 157 158 Collective on DA 159 160 Input Parameters: 161 + da - the distributed array 162 - ptr - an optional visualization context 163 164 Notes: 165 The available visualization contexts include 166 + PETSC_VIEWER_STDOUT_SELF - standard output (default) 167 . PETSC_VIEWER_STDOUT_WORLD - synchronized standard 168 output where only the first processor opens 169 the file. All other processors send their 170 data to the first processor to print. 171 - PETSC_VIEWER_DRAW_WORLD - to default window 172 173 The user can open alternative visualization contexts with 174 + PetscViewerASCIIOpen() - Outputs vector to a specified file 175 - PetscViewerDrawOpen() - Outputs vector to an X window display 176 177 Default Output Format: 178 (for 3d arrays) 179 .vb 180 Processor [proc] M N P m n p w s 181 X range: xs xe, Y range: ys, ye, Z range: zs, ze 182 183 where 184 M,N,P - global dimension in each direction of the array 185 m,n,p - corresponding number of procs in each dimension 186 w - number of degrees of freedom per node 187 s - stencil width 188 xs, xe - internal local starting/ending grid points 189 in x-direction, (augmented to handle multiple 190 degrees of freedom per node) 191 ys, ye - local starting/ending grid points in y-direction 192 zs, ze - local starting/ending grid points in z-direction 193 .ve 194 195 Options Database Key: 196 . -da_view - Calls DAView() at the conclusion of DACreate1d(), 197 DACreate2d(), and DACreate3d() 198 199 Level: beginner 200 201 Notes: 202 Use DAGetCorners() and DAGetGhostCorners() to get the starting 203 and ending grid points (ghost points) in each direction. 204 205 When drawing the DA grid it only draws the logical grid and does not 206 respect the grid coordinates set with DASetCoordinates() 207 208 .keywords: distributed array, view, visualize 209 210 .seealso: PetscViewerASCIIOpen(), PetscViewerDrawOpen(), DAGetInfo(), DAGetCorners(), 211 DAGetGhostCorners(), DAGetOwnershipRanges(), DACreate(), DACreate1d(), DACreate2d(), DACreate3d() 212 @*/ 213 PetscErrorCode PETSCDM_DLLEXPORT DAView(DA da,PetscViewer viewer) 214 { 215 PetscErrorCode ierr; 216 DM_DA *dd = (DM_DA*)da->data; 217 PetscInt i,dof = dd->w; 218 PetscBool iascii,fieldsnamed = PETSC_FALSE,isbinary; 219 #if defined(PETSC_HAVE_MATLAB_ENGINE) 220 PetscBool ismatlab; 221 #endif 222 223 PetscFunctionBegin; 224 PetscValidHeaderSpecific(da,DM_CLASSID,1); 225 if (!viewer) { 226 ierr = PetscViewerASCIIGetStdout(((PetscObject)da)->comm,&viewer);CHKERRQ(ierr); 227 } 228 PetscValidHeaderSpecific(viewer,PETSC_VIEWER_CLASSID,2); 229 230 ierr = PetscTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 231 ierr = PetscTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);CHKERRQ(ierr); 232 #if defined(PETSC_HAVE_MATLAB_ENGINE) 233 ierr = PetscTypeCompare((PetscObject)viewer,PETSCVIEWERMATLAB,&ismatlab);CHKERRQ(ierr); 234 #endif 235 if (iascii) { 236 PetscViewerFormat format; 237 238 ierr = PetscViewerGetFormat(viewer, &format);CHKERRQ(ierr); 239 if (format == PETSC_VIEWER_ASCII_VTK || format == PETSC_VIEWER_ASCII_VTK_CELL) { 240 ierr = DAView_VTK(da, viewer);CHKERRQ(ierr); 241 } else { 242 ierr = PetscObjectPrintClassNamePrefixType((PetscObject)da,viewer,"DA Object");CHKERRQ(ierr); 243 for (i=0; i<dof; i++) { 244 if (dd->fieldname[i]) { 245 fieldsnamed = PETSC_TRUE; 246 break; 247 } 248 } 249 if (fieldsnamed) { 250 ierr = PetscViewerASCIIPrintf(viewer,"FieldNames: ");CHKERRQ(ierr); 251 for (i=0; i<dof; i++) { 252 if (dd->fieldname[i]) { 253 ierr = PetscViewerASCIIPrintf(viewer,"%s ",dd->fieldname[i]);CHKERRQ(ierr); 254 } else { 255 ierr = PetscViewerASCIIPrintf(viewer,"(not named) ");CHKERRQ(ierr); 256 } 257 } 258 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 259 } 260 } 261 } 262 if (isbinary){ 263 ierr = DAView_Binary(da,viewer);CHKERRQ(ierr); 264 #if defined(PETSC_HAVE_MATLAB_ENGINE) 265 } else if (ismatlab) { 266 ierr = DAView_Matlab(da,viewer);CHKERRQ(ierr); 267 #endif 268 } else { 269 ierr = (*da->ops->view)(da,viewer);CHKERRQ(ierr); 270 } 271 PetscFunctionReturn(0); 272 } 273 274 #undef __FUNCT__ 275 #define __FUNCT__ "DAGetInfo" 276 /*@C 277 DAGetInfo - Gets information about a given distributed array. 278 279 Not Collective 280 281 Input Parameter: 282 . da - the distributed array 283 284 Output Parameters: 285 + dim - dimension of the distributed array (1, 2, or 3) 286 . M, N, P - global dimension in each direction of the array 287 . m, n, p - corresponding number of procs in each dimension 288 . dof - number of degrees of freedom per node 289 . s - stencil width 290 . wrap - type of periodicity, one of DA_NONPERIODIC, DA_XPERIODIC, DA_YPERIODIC, 291 DA_XYPERIODIC, DA_XYZPERIODIC, DA_XZPERIODIC, DA_YZPERIODIC,DA_ZPERIODIC 292 - st - stencil type, either DA_STENCIL_STAR or DA_STENCIL_BOX 293 294 Level: beginner 295 296 Note: 297 Use PETSC_NULL (PETSC_NULL_INTEGER in Fortran) in place of any output parameter that is not of interest. 298 299 .keywords: distributed array, get, information 300 301 .seealso: DAView(), DAGetCorners(), DAGetLocalInfo() 302 @*/ 303 PetscErrorCode PETSCDM_DLLEXPORT DAGetInfo(DA da,PetscInt *dim,PetscInt *M,PetscInt *N,PetscInt *P,PetscInt *m,PetscInt *n,PetscInt *p,PetscInt *dof,PetscInt *s,DAPeriodicType *wrap,DAStencilType *st) 304 { 305 DM_DA *dd = (DM_DA*)da->data; 306 307 PetscFunctionBegin; 308 PetscValidHeaderSpecific(da,DM_CLASSID,1); 309 if (dim) *dim = dd->dim; 310 if (M) *M = dd->M; 311 if (N) *N = dd->N; 312 if (P) *P = dd->P; 313 if (m) *m = dd->m; 314 if (n) *n = dd->n; 315 if (p) *p = dd->p; 316 if (dof) *dof = dd->w; 317 if (s) *s = dd->s; 318 if (wrap) *wrap = dd->wrap; 319 if (st) *st = dd->stencil_type; 320 PetscFunctionReturn(0); 321 } 322 323 #undef __FUNCT__ 324 #define __FUNCT__ "DAGetLocalInfo" 325 /*@C 326 DAGetLocalInfo - Gets information about a given distributed array and this processors location in it 327 328 Not Collective 329 330 Input Parameter: 331 . da - the distributed array 332 333 Output Parameters: 334 . dainfo - structure containing the information 335 336 Level: beginner 337 338 .keywords: distributed array, get, information 339 340 .seealso: DAGetInfo(), DAGetCorners() 341 @*/ 342 PetscErrorCode PETSCDM_DLLEXPORT DAGetLocalInfo(DA da,DALocalInfo *info) 343 { 344 PetscInt w; 345 DM_DA *dd = (DM_DA*)da->data; 346 347 PetscFunctionBegin; 348 PetscValidHeaderSpecific(da,DM_CLASSID,1); 349 PetscValidPointer(info,2); 350 info->da = da; 351 info->dim = dd->dim; 352 info->mx = dd->M; 353 info->my = dd->N; 354 info->mz = dd->P; 355 info->dof = dd->w; 356 info->sw = dd->s; 357 info->pt = dd->wrap; 358 info->st = dd->stencil_type; 359 360 /* since the xs, xe ... have all been multiplied by the number of degrees 361 of freedom per cell, w = dd->w, we divide that out before returning.*/ 362 w = dd->w; 363 info->xs = dd->xs/w; 364 info->xm = (dd->xe - dd->xs)/w; 365 /* the y and z have NOT been multiplied by w */ 366 info->ys = dd->ys; 367 info->ym = (dd->ye - dd->ys); 368 info->zs = dd->zs; 369 info->zm = (dd->ze - dd->zs); 370 371 info->gxs = dd->Xs/w; 372 info->gxm = (dd->Xe - dd->Xs)/w; 373 /* the y and z have NOT been multiplied by w */ 374 info->gys = dd->Ys; 375 info->gym = (dd->Ye - dd->Ys); 376 info->gzs = dd->Zs; 377 info->gzm = (dd->Ze - dd->Zs); 378 PetscFunctionReturn(0); 379 } 380 381