1 2 /* 3 Code for manipulating distributed regular arrays in parallel. 4 */ 5 6 #include <petsc/private/dmdaimpl.h> /*I "petscdmda.h" I*/ 7 8 #if defined(PETSC_HAVE_MATLAB_ENGINE) 9 #include <mat.h> /* MATLAB include file */ 10 11 PetscErrorCode DMView_DA_Matlab(DM da,PetscViewer viewer) 12 { 13 PetscErrorCode ierr; 14 PetscMPIInt rank; 15 PetscInt dim,m,n,p,dof,swidth; 16 DMDAStencilType stencil; 17 DMBoundaryType bx,by,bz; 18 mxArray *mx; 19 const char *fnames[] = {"dimension","m","n","p","dof","stencil_width","bx","by","bz","stencil_type"}; 20 21 PetscFunctionBegin; 22 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)da),&rank);CHKERRMPI(ierr); 23 if (rank == 0) { 24 ierr = DMDAGetInfo(da,&dim,&m,&n,&p,0,0,0,&dof,&swidth,&bx,&by,&bz,&stencil);CHKERRQ(ierr); 25 mx = mxCreateStructMatrix(1,1,8,(const char**)fnames); 26 if (!mx) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Unable to generate MATLAB struct array to hold DMDA information"); 27 mxSetFieldByNumber(mx,0,0,mxCreateDoubleScalar((double)dim)); 28 mxSetFieldByNumber(mx,0,1,mxCreateDoubleScalar((double)m)); 29 mxSetFieldByNumber(mx,0,2,mxCreateDoubleScalar((double)n)); 30 mxSetFieldByNumber(mx,0,3,mxCreateDoubleScalar((double)p)); 31 mxSetFieldByNumber(mx,0,4,mxCreateDoubleScalar((double)dof)); 32 mxSetFieldByNumber(mx,0,5,mxCreateDoubleScalar((double)swidth)); 33 mxSetFieldByNumber(mx,0,6,mxCreateDoubleScalar((double)bx)); 34 mxSetFieldByNumber(mx,0,7,mxCreateDoubleScalar((double)by)); 35 mxSetFieldByNumber(mx,0,8,mxCreateDoubleScalar((double)bz)); 36 mxSetFieldByNumber(mx,0,9,mxCreateDoubleScalar((double)stencil)); 37 ierr = PetscObjectName((PetscObject)da);CHKERRQ(ierr); 38 ierr = PetscViewerMatlabPutVariable(viewer,((PetscObject)da)->name,mx);CHKERRQ(ierr); 39 } 40 PetscFunctionReturn(0); 41 } 42 #endif 43 44 PetscErrorCode DMView_DA_Binary(DM da,PetscViewer viewer) 45 { 46 PetscErrorCode ierr; 47 PetscMPIInt rank; 48 PetscInt dim,m,n,p,dof,swidth,M,N,P; 49 DMDAStencilType stencil; 50 DMBoundaryType bx,by,bz; 51 MPI_Comm comm; 52 PetscBool coors = PETSC_FALSE; 53 54 PetscFunctionBegin; 55 ierr = PetscObjectGetComm((PetscObject)da,&comm);CHKERRQ(ierr); 56 57 ierr = DMDAGetInfo(da,&dim,&m,&n,&p,&M,&N,&P,&dof,&swidth,&bx,&by,&bz,&stencil);CHKERRQ(ierr); 58 ierr = MPI_Comm_rank(comm,&rank);CHKERRMPI(ierr); 59 if (rank == 0) { 60 ierr = PetscViewerBinaryWrite(viewer,&dim,1,PETSC_INT);CHKERRQ(ierr); 61 ierr = PetscViewerBinaryWrite(viewer,&m,1,PETSC_INT);CHKERRQ(ierr); 62 ierr = PetscViewerBinaryWrite(viewer,&n,1,PETSC_INT);CHKERRQ(ierr); 63 ierr = PetscViewerBinaryWrite(viewer,&p,1,PETSC_INT);CHKERRQ(ierr); 64 ierr = PetscViewerBinaryWrite(viewer,&dof,1,PETSC_INT);CHKERRQ(ierr); 65 ierr = PetscViewerBinaryWrite(viewer,&swidth,1,PETSC_INT);CHKERRQ(ierr); 66 ierr = PetscViewerBinaryWrite(viewer,&bx,1,PETSC_ENUM);CHKERRQ(ierr); 67 ierr = PetscViewerBinaryWrite(viewer,&by,1,PETSC_ENUM);CHKERRQ(ierr); 68 ierr = PetscViewerBinaryWrite(viewer,&bz,1,PETSC_ENUM);CHKERRQ(ierr); 69 ierr = PetscViewerBinaryWrite(viewer,&stencil,1,PETSC_ENUM);CHKERRQ(ierr); 70 if (da->coordinates) coors = PETSC_TRUE; 71 ierr = PetscViewerBinaryWrite(viewer,&coors,1,PETSC_BOOL);CHKERRQ(ierr); 72 } 73 74 /* save the coordinates if they exist to disk (in the natural ordering) */ 75 if (da->coordinates) { 76 ierr = VecView(da->coordinates,viewer);CHKERRQ(ierr); 77 } 78 PetscFunctionReturn(0); 79 } 80 81 PetscErrorCode DMView_DA_VTK(DM da, PetscViewer viewer) 82 { 83 PetscInt dim, dof, M = 0, N = 0, P = 0; 84 PetscErrorCode ierr; 85 86 PetscFunctionBegin; 87 ierr = DMDAGetInfo(da, &dim, &M, &N, &P, NULL, NULL, NULL, &dof, NULL, NULL, NULL, NULL, NULL);CHKERRQ(ierr); 88 if (!da->coordinates) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP, "VTK output requires DMDA coordinates."); 89 /* Write Header */ 90 ierr = PetscViewerASCIIPrintf(viewer,"# vtk DataFile Version 2.0\n");CHKERRQ(ierr); 91 ierr = PetscViewerASCIIPrintf(viewer,"Structured Mesh Example\n");CHKERRQ(ierr); 92 ierr = PetscViewerASCIIPrintf(viewer,"ASCII\n");CHKERRQ(ierr); 93 ierr = PetscViewerASCIIPrintf(viewer,"DATASET STRUCTURED_GRID\n");CHKERRQ(ierr); 94 ierr = PetscViewerASCIIPrintf(viewer,"DIMENSIONS %d %d %d\n", M, N, P);CHKERRQ(ierr); 95 ierr = PetscViewerASCIIPrintf(viewer,"POINTS %d double\n", M*N*P);CHKERRQ(ierr); 96 if (da->coordinates) { 97 DM dac; 98 Vec natural; 99 100 ierr = DMGetCoordinateDM(da, &dac);CHKERRQ(ierr); 101 ierr = DMDACreateNaturalVector(dac, &natural);CHKERRQ(ierr); 102 ierr = PetscObjectSetOptionsPrefix((PetscObject) natural, "coor_");CHKERRQ(ierr); 103 ierr = DMDAGlobalToNaturalBegin(dac, da->coordinates, INSERT_VALUES, natural);CHKERRQ(ierr); 104 ierr = DMDAGlobalToNaturalEnd(dac, da->coordinates, INSERT_VALUES, natural);CHKERRQ(ierr); 105 ierr = PetscViewerPushFormat(viewer, PETSC_VIEWER_ASCII_VTK_COORDS_DEPRECATED);CHKERRQ(ierr); 106 ierr = VecView(natural, viewer);CHKERRQ(ierr); 107 ierr = PetscViewerPopFormat(viewer);CHKERRQ(ierr); 108 ierr = VecDestroy(&natural);CHKERRQ(ierr); 109 } 110 PetscFunctionReturn(0); 111 } 112 113 /*@C 114 DMDAGetInfo - Gets information about a given distributed array. 115 116 Not Collective 117 118 Input Parameter: 119 . da - the distributed array 120 121 Output Parameters: 122 + dim - dimension of the distributed array (1, 2, or 3) 123 . M - global dimension in first direction of the array 124 . N - global dimension in second direction of the array 125 . P - global dimension in third direction of the array 126 . m - corresponding number of procs in first dimension 127 . n - corresponding number of procs in second dimension 128 . p - corresponding number of procs in third dimension 129 . dof - number of degrees of freedom per node 130 . s - stencil width 131 . bx - type of ghost nodes at boundary in first dimension 132 . by - type of ghost nodes at boundary in second dimension 133 . bz - type of ghost nodes at boundary in third dimension 134 - st - stencil type, either DMDA_STENCIL_STAR or DMDA_STENCIL_BOX 135 136 Level: beginner 137 138 Note: 139 Use NULL (NULL_INTEGER in Fortran) in place of any output parameter that is not of interest. 140 141 .seealso: DMView(), DMDAGetCorners(), DMDAGetLocalInfo() 142 @*/ 143 PetscErrorCode DMDAGetInfo(DM da,PetscInt *dim,PetscInt *M,PetscInt *N,PetscInt *P,PetscInt *m,PetscInt *n,PetscInt *p,PetscInt *dof,PetscInt *s,DMBoundaryType *bx,DMBoundaryType *by,DMBoundaryType *bz,DMDAStencilType *st) 144 { 145 DM_DA *dd = (DM_DA*)da->data; 146 147 PetscFunctionBegin; 148 PetscValidHeaderSpecificType(da,DM_CLASSID,1,DMDA); 149 if (dim) *dim = da->dim; 150 if (M) { 151 if (dd->Mo < 0) *M = dd->M; 152 else *M = dd->Mo; 153 } 154 if (N) { 155 if (dd->No < 0) *N = dd->N; 156 else *N = dd->No; 157 } 158 if (P) { 159 if (dd->Po < 0) *P = dd->P; 160 else *P = dd->Po; 161 } 162 if (m) *m = dd->m; 163 if (n) *n = dd->n; 164 if (p) *p = dd->p; 165 if (dof) *dof = dd->w; 166 if (s) *s = dd->s; 167 if (bx) *bx = dd->bx; 168 if (by) *by = dd->by; 169 if (bz) *bz = dd->bz; 170 if (st) *st = dd->stencil_type; 171 PetscFunctionReturn(0); 172 } 173 174 /*@C 175 DMDAGetLocalInfo - Gets information about a given distributed array and this processors location in it 176 177 Not Collective 178 179 Input Parameter: 180 . da - the distributed array 181 182 Output Parameters: 183 . dainfo - structure containing the information 184 185 Level: beginner 186 187 Notes: 188 See DMDALocalInfo for the information that is returned 189 190 .seealso: DMDAGetInfo(), DMDAGetCorners(), DMDALocalInfo 191 @*/ 192 PetscErrorCode DMDAGetLocalInfo(DM da,DMDALocalInfo *info) 193 { 194 PetscInt w; 195 DM_DA *dd = (DM_DA*)da->data; 196 197 PetscFunctionBegin; 198 PetscValidHeaderSpecificType(da,DM_CLASSID,1,DMDA); 199 PetscValidPointer(info,2); 200 info->da = da; 201 info->dim = da->dim; 202 if (dd->Mo < 0) info->mx = dd->M; 203 else info->mx = dd->Mo; 204 if (dd->No < 0) info->my = dd->N; 205 else info->my = dd->No; 206 if (dd->Po < 0) info->mz = dd->P; 207 else info->mz = dd->Po; 208 info->dof = dd->w; 209 info->sw = dd->s; 210 info->bx = dd->bx; 211 info->by = dd->by; 212 info->bz = dd->bz; 213 info->st = dd->stencil_type; 214 215 /* since the xs, xe ... have all been multiplied by the number of degrees 216 of freedom per cell, w = dd->w, we divide that out before returning.*/ 217 w = dd->w; 218 info->xs = dd->xs/w + dd->xo; 219 info->xm = (dd->xe - dd->xs)/w; 220 /* the y and z have NOT been multiplied by w */ 221 info->ys = dd->ys + dd->yo; 222 info->ym = (dd->ye - dd->ys); 223 info->zs = dd->zs + dd->zo; 224 info->zm = (dd->ze - dd->zs); 225 226 info->gxs = dd->Xs/w + dd->xo; 227 info->gxm = (dd->Xe - dd->Xs)/w; 228 /* the y and z have NOT been multiplied by w */ 229 info->gys = dd->Ys + dd->yo; 230 info->gym = (dd->Ye - dd->Ys); 231 info->gzs = dd->Zs + dd->zo; 232 info->gzm = (dd->Ze - dd->Zs); 233 PetscFunctionReturn(0); 234 } 235