1 #define PETSCDM_DLL 2 3 /* 4 Code for manipulating distributed regular arrays in parallel. 5 */ 6 7 #include "private/daimpl.h" /*I "petscdm.h" I*/ 8 9 #if defined(PETSC_HAVE_MATLAB_ENGINE) 10 #include "mat.h" /* MATLAB include file */ 11 12 #undef __FUNCT__ 13 #define __FUNCT__ "DMView_DA_Matlab" 14 PetscErrorCode DMView_DA_Matlab(DM da,PetscViewer viewer) 15 { 16 PetscErrorCode ierr; 17 PetscMPIInt rank; 18 PetscInt dim,m,n,p,dof,swidth; 19 DMDAStencilType stencil; 20 DMDABoundaryType bx,by,bz; 21 mxArray *mx; 22 const char *fnames[] = {"dimension","m","n","p","dof","stencil_width","bx","by","bz","stencil_type"}; 23 24 PetscFunctionBegin; 25 ierr = MPI_Comm_rank(((PetscObject)da)->comm,&rank);CHKERRQ(ierr); 26 if (!rank) { 27 ierr = DMDAGetInfo(da,&dim,&m,&n,&p,0,0,0,&dof,&swidth,&bx,&by,&bz,&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 DMDA 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)bx)); 37 mxSetFieldByNumber(mx,0,7,mxCreateDoubleScalar((double)by)); 38 mxSetFieldByNumber(mx,0,8,mxCreateDoubleScalar((double)bz)); 39 mxSetFieldByNumber(mx,0,9,mxCreateDoubleScalar((double)stencil)); 40 ierr = PetscObjectName((PetscObject)da);CHKERRQ(ierr); 41 ierr = PetscViewerMatlabPutVariable(viewer,((PetscObject)da)->name,mx);CHKERRQ(ierr); 42 } 43 PetscFunctionReturn(0); 44 } 45 #endif 46 47 #undef __FUNCT__ 48 #define __FUNCT__ "DMView_DA_Binary" 49 PetscErrorCode DMView_DA_Binary(DM da,PetscViewer viewer) 50 { 51 PetscErrorCode ierr; 52 PetscMPIInt rank; 53 PetscInt i,dim,m,n,p,dof,swidth,M,N,P; 54 size_t j,len; 55 DMDAStencilType stencil; 56 DMDABoundaryType bx,by,bz; 57 MPI_Comm comm; 58 DM_DA *dd = (DM_DA*)da->data; 59 60 PetscFunctionBegin; 61 ierr = PetscObjectGetComm((PetscObject)da,&comm);CHKERRQ(ierr); 62 63 ierr = DMDAGetInfo(da,&dim,&m,&n,&p,&M,&N,&P,&dof,&swidth,&bx,&by,&bz,&stencil);CHKERRQ(ierr); 64 ierr = MPI_Comm_rank(comm,&rank);CHKERRQ(ierr); 65 if (!rank) { 66 FILE *file; 67 68 ierr = PetscViewerBinaryGetInfoPointer(viewer,&file);CHKERRQ(ierr); 69 if (file) { 70 char fieldname[PETSC_MAX_PATH_LEN]; 71 72 ierr = PetscFPrintf(PETSC_COMM_SELF,file,"-daload_info %D,%D,%D,%D,%D,%D,%D,%D,%D,%D\n",dim,m,n,p,dof,swidth,stencil,bx,by,bz);CHKERRQ(ierr); 73 for (i=0; i<dof; i++) { 74 if (dd->fieldname[i]) { 75 ierr = PetscStrncpy(fieldname,dd->fieldname[i],PETSC_MAX_PATH_LEN);CHKERRQ(ierr); 76 ierr = PetscStrlen(fieldname,&len);CHKERRQ(ierr); 77 len = PetscMin(PETSC_MAX_PATH_LEN,len);CHKERRQ(ierr); 78 for (j=0; j<len; j++) { 79 if (fieldname[j] == ' ') fieldname[j] = '_'; 80 } 81 ierr = PetscFPrintf(PETSC_COMM_SELF,file,"-daload_fieldname_%D %s\n",i,fieldname);CHKERRQ(ierr); 82 } 83 } 84 if (dd->coordinates) { /* save the DMDA's coordinates */ 85 ierr = PetscFPrintf(PETSC_COMM_SELF,file,"-daload_coordinates\n");CHKERRQ(ierr); 86 } 87 } 88 } 89 90 /* save the coordinates if they exist to disk (in the natural ordering) */ 91 if (dd->coordinates) { 92 DM dac; 93 const PetscInt *lx,*ly,*lz; 94 Vec natural; 95 96 /* create the appropriate DMDA to map to natural ordering */ 97 ierr = DMDAGetOwnershipRanges(da,&lx,&ly,&lz);CHKERRQ(ierr); 98 if (dim == 1) { 99 ierr = DMDACreate1d(comm,DMDA_BOUNDARY_NONE,m,dim,0,lx,&dac);CHKERRQ(ierr); 100 } else if (dim == 2) { 101 ierr = DMDACreate2d(comm,DMDA_BOUNDARY_NONE,DMDA_BOUNDARY_NONE,DMDA_STENCIL_BOX,m,n,M,N,dim,0,lx,ly,&dac);CHKERRQ(ierr); 102 } else if (dim == 3) { 103 ierr = DMDACreate3d(comm,DMDA_BOUNDARY_NONE,DMDA_BOUNDARY_NONE,DMDA_BOUNDARY_NONE,DMDA_STENCIL_BOX,m,n,p,M,N,P,dim,0,lx,ly,lz,&dac);CHKERRQ(ierr); 104 } else SETERRQ1(comm,PETSC_ERR_ARG_CORRUPT,"Dimension is not 1 2 or 3: %D\n",dim); 105 ierr = DMDACreateNaturalVector(dac,&natural);CHKERRQ(ierr); 106 ierr = PetscObjectSetOptionsPrefix((PetscObject)natural,"coor_");CHKERRQ(ierr); 107 ierr = DMDAGlobalToNaturalBegin(dac,dd->coordinates,INSERT_VALUES,natural);CHKERRQ(ierr); 108 ierr = DMDAGlobalToNaturalEnd(dac,dd->coordinates,INSERT_VALUES,natural);CHKERRQ(ierr); 109 ierr = VecView(natural,viewer);CHKERRQ(ierr); 110 ierr = VecDestroy(natural);CHKERRQ(ierr); 111 ierr = DMDestroy(dac);CHKERRQ(ierr); 112 } 113 PetscFunctionReturn(0); 114 } 115 116 #undef __FUNCT__ 117 #define __FUNCT__ "DMView_DA_VTK" 118 PetscErrorCode DMView_DA_VTK(DM da, PetscViewer viewer) 119 { 120 PetscInt dim, dof, M = 0, N = 0, P = 0; 121 PetscErrorCode ierr; 122 DM_DA *dd = (DM_DA*)da->data; 123 124 PetscFunctionBegin; 125 ierr = DMDAGetInfo(da, &dim, &M, &N, &P, PETSC_NULL, PETSC_NULL, PETSC_NULL, &dof, PETSC_NULL, PETSC_NULL, PETSC_NULL, PETSC_NULL, PETSC_NULL);CHKERRQ(ierr); 126 /* if (dim != 3) {SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP, "VTK output only works for three dimensional DMDAs.");} */ 127 if (!dd->coordinates) SETERRQ(((PetscObject)da)->comm,PETSC_ERR_SUP, "VTK output requires DMDA coordinates."); 128 /* Write Header */ 129 ierr = PetscViewerASCIIPrintf(viewer,"# vtk DataFile Version 2.0\n");CHKERRQ(ierr); 130 ierr = PetscViewerASCIIPrintf(viewer,"Structured Mesh Example\n");CHKERRQ(ierr); 131 ierr = PetscViewerASCIIPrintf(viewer,"ASCII\n");CHKERRQ(ierr); 132 ierr = PetscViewerASCIIPrintf(viewer,"DATASET STRUCTURED_GRID\n");CHKERRQ(ierr); 133 ierr = PetscViewerASCIIPrintf(viewer,"DIMENSIONS %d %d %d\n", M, N, P);CHKERRQ(ierr); 134 ierr = PetscViewerASCIIPrintf(viewer,"POINTS %d double\n", M*N*P);CHKERRQ(ierr); 135 if (dd->coordinates) { 136 DM dac; 137 Vec natural; 138 139 ierr = DMDAGetCoordinateDA(da, &dac);CHKERRQ(ierr); 140 ierr = DMDACreateNaturalVector(dac, &natural);CHKERRQ(ierr); 141 ierr = PetscObjectSetOptionsPrefix((PetscObject) natural, "coor_");CHKERRQ(ierr); 142 ierr = DMDAGlobalToNaturalBegin(dac, dd->coordinates, INSERT_VALUES, natural);CHKERRQ(ierr); 143 ierr = DMDAGlobalToNaturalEnd(dac, dd->coordinates, INSERT_VALUES, natural);CHKERRQ(ierr); 144 ierr = PetscViewerPushFormat(viewer, PETSC_VIEWER_ASCII_VTK_COORDS);CHKERRQ(ierr); 145 ierr = VecView(natural, viewer);CHKERRQ(ierr); 146 ierr = PetscViewerPopFormat(viewer);CHKERRQ(ierr); 147 ierr = VecDestroy(natural);CHKERRQ(ierr); 148 } 149 PetscFunctionReturn(0); 150 } 151 152 #undef __FUNCT__ 153 #define __FUNCT__ "DMDAGetInfo" 154 /*@C 155 DMDAGetInfo - Gets information about a given distributed array. 156 157 Not Collective 158 159 Input Parameter: 160 . da - the distributed array 161 162 Output Parameters: 163 + dim - dimension of the distributed array (1, 2, or 3) 164 . M, N, P - global dimension in each direction of the array 165 . m, n, p - corresponding number of procs in each dimension 166 . dof - number of degrees of freedom per node 167 . s - stencil width 168 . bx,by,bz - type of ghost nodes at boundary, one of DMDA_BOUNDARY_NONE, DMDA_BOUNDARY_GHOSTED, 169 DMDA_BOUNDARY_MIRROR, DMDA_BOUNDARY_PERIODIC 170 - st - stencil type, either DMDA_STENCIL_STAR or DMDA_STENCIL_BOX 171 172 Level: beginner 173 174 Note: 175 Use PETSC_NULL (PETSC_NULL_INTEGER in Fortran) in place of any output parameter that is not of interest. 176 177 .keywords: distributed array, get, information 178 179 .seealso: DMView(), DMDAGetCorners(), DMDAGetLocalInfo() 180 @*/ 181 PetscErrorCode DMDAGetInfo(DM da,PetscInt *dim,PetscInt *M,PetscInt *N,PetscInt *P,PetscInt *m,PetscInt *n,PetscInt *p,PetscInt *dof,PetscInt *s,DMDABoundaryType *bx,DMDABoundaryType *by,DMDABoundaryType *bz,DMDAStencilType *st) 182 { 183 DM_DA *dd = (DM_DA*)da->data; 184 185 PetscFunctionBegin; 186 PetscValidHeaderSpecific(da,DM_CLASSID,1); 187 if (dim) *dim = dd->dim; 188 if (M) *M = dd->M; 189 if (N) *N = dd->N; 190 if (P) *P = dd->P; 191 if (m) *m = dd->m; 192 if (n) *n = dd->n; 193 if (p) *p = dd->p; 194 if (dof) *dof = dd->w; 195 if (s) *s = dd->s; 196 if (bx) *bx = dd->bx; 197 if (by) *by = dd->by; 198 if (bz) *bz = dd->bz; 199 if (st) *st = dd->stencil_type; 200 PetscFunctionReturn(0); 201 } 202 203 #undef __FUNCT__ 204 #define __FUNCT__ "DMDAGetLocalInfo" 205 /*@C 206 DMDAGetLocalInfo - Gets information about a given distributed array and this processors location in it 207 208 Not Collective 209 210 Input Parameter: 211 . da - the distributed array 212 213 Output Parameters: 214 . dainfo - structure containing the information 215 216 Level: beginner 217 218 .keywords: distributed array, get, information 219 220 .seealso: DMDAGetInfo(), DMDAGetCorners() 221 @*/ 222 PetscErrorCode DMDAGetLocalInfo(DM da,DMDALocalInfo *info) 223 { 224 PetscInt w; 225 DM_DA *dd = (DM_DA*)da->data; 226 227 PetscFunctionBegin; 228 PetscValidHeaderSpecific(da,DM_CLASSID,1); 229 PetscValidPointer(info,2); 230 info->da = da; 231 info->dim = dd->dim; 232 info->mx = dd->M; 233 info->my = dd->N; 234 info->mz = dd->P; 235 info->dof = dd->w; 236 info->sw = dd->s; 237 info->bx = dd->bx; 238 info->by = dd->by; 239 info->bz = dd->bz; 240 info->st = dd->stencil_type; 241 242 /* since the xs, xe ... have all been multiplied by the number of degrees 243 of freedom per cell, w = dd->w, we divide that out before returning.*/ 244 w = dd->w; 245 info->xs = dd->xs/w; 246 info->xm = (dd->xe - dd->xs)/w; 247 /* the y and z have NOT been multiplied by w */ 248 info->ys = dd->ys; 249 info->ym = (dd->ye - dd->ys); 250 info->zs = dd->zs; 251 info->zm = (dd->ze - dd->zs); 252 253 info->gxs = dd->Xs/w; 254 info->gxm = (dd->Xe - dd->Xs)/w; 255 /* the y and z have NOT been multiplied by w */ 256 info->gys = dd->Ys; 257 info->gym = (dd->Ye - dd->Ys); 258 info->gzs = dd->Zs; 259 info->gzm = (dd->Ze - dd->Zs); 260 PetscFunctionReturn(0); 261 } 262 263