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