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) PetscCall(VecView(da->coordinates,viewer)); 74 PetscFunctionReturn(0); 75 } 76 77 PetscErrorCode DMView_DA_VTK(DM da, PetscViewer viewer) 78 { 79 PetscInt dim, dof, M = 0, N = 0, P = 0; 80 81 PetscFunctionBegin; 82 PetscCall(DMDAGetInfo(da, &dim, &M, &N, &P, NULL, NULL, NULL, &dof, NULL, NULL, NULL, NULL, NULL)); 83 PetscCheck(da->coordinates,PetscObjectComm((PetscObject)da),PETSC_ERR_SUP, "VTK output requires DMDA coordinates."); 84 /* Write Header */ 85 PetscCall(PetscViewerASCIIPrintf(viewer,"# vtk DataFile Version 2.0\n")); 86 PetscCall(PetscViewerASCIIPrintf(viewer,"Structured Mesh Example\n")); 87 PetscCall(PetscViewerASCIIPrintf(viewer,"ASCII\n")); 88 PetscCall(PetscViewerASCIIPrintf(viewer,"DATASET STRUCTURED_GRID\n")); 89 PetscCall(PetscViewerASCIIPrintf(viewer,"DIMENSIONS %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT "\n", M, N, P)); 90 PetscCall(PetscViewerASCIIPrintf(viewer,"POINTS %" PetscInt_FMT " double\n", M*N*P)); 91 if (da->coordinates) { 92 DM dac; 93 Vec natural; 94 95 PetscCall(DMGetCoordinateDM(da, &dac)); 96 PetscCall(DMDACreateNaturalVector(dac, &natural)); 97 PetscCall(PetscObjectSetOptionsPrefix((PetscObject) natural, "coor_")); 98 PetscCall(DMDAGlobalToNaturalBegin(dac, da->coordinates, INSERT_VALUES, natural)); 99 PetscCall(DMDAGlobalToNaturalEnd(dac, da->coordinates, INSERT_VALUES, natural)); 100 PetscCall(PetscViewerPushFormat(viewer, PETSC_VIEWER_ASCII_VTK_COORDS_DEPRECATED)); 101 PetscCall(VecView(natural, viewer)); 102 PetscCall(PetscViewerPopFormat(viewer)); 103 PetscCall(VecDestroy(&natural)); 104 } 105 PetscFunctionReturn(0); 106 } 107 108 /*@C 109 DMDAGetInfo - Gets information about a given distributed array. 110 111 Not Collective 112 113 Input Parameter: 114 . da - the distributed array 115 116 Output Parameters: 117 + dim - dimension of the distributed array (1, 2, or 3) 118 . M - global dimension in first direction of the array 119 . N - global dimension in second direction of the array 120 . P - global dimension in third direction of the array 121 . m - corresponding number of procs in first dimension 122 . n - corresponding number of procs in second dimension 123 . p - corresponding number of procs in third dimension 124 . dof - number of degrees of freedom per node 125 . s - stencil width 126 . bx - type of ghost nodes at boundary in first dimension 127 . by - type of ghost nodes at boundary in second dimension 128 . bz - type of ghost nodes at boundary in third dimension 129 - st - stencil type, either DMDA_STENCIL_STAR or DMDA_STENCIL_BOX 130 131 Level: beginner 132 133 Note: 134 Use NULL (NULL_INTEGER in Fortran) in place of any output parameter that is not of interest. 135 136 .seealso: `DMView()`, `DMDAGetCorners()`, `DMDAGetLocalInfo()` 137 @*/ 138 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) 139 { 140 DM_DA *dd = (DM_DA*)da->data; 141 142 PetscFunctionBegin; 143 PetscValidHeaderSpecificType(da,DM_CLASSID,1,DMDA); 144 if (dim) *dim = da->dim; 145 if (M) { 146 if (dd->Mo < 0) *M = dd->M; 147 else *M = dd->Mo; 148 } 149 if (N) { 150 if (dd->No < 0) *N = dd->N; 151 else *N = dd->No; 152 } 153 if (P) { 154 if (dd->Po < 0) *P = dd->P; 155 else *P = dd->Po; 156 } 157 if (m) *m = dd->m; 158 if (n) *n = dd->n; 159 if (p) *p = dd->p; 160 if (dof) *dof = dd->w; 161 if (s) *s = dd->s; 162 if (bx) *bx = dd->bx; 163 if (by) *by = dd->by; 164 if (bz) *bz = dd->bz; 165 if (st) *st = dd->stencil_type; 166 PetscFunctionReturn(0); 167 } 168 169 /*@C 170 DMDAGetLocalInfo - Gets information about a given distributed array and this processors location in it 171 172 Not Collective 173 174 Input Parameter: 175 . da - the distributed array 176 177 Output Parameters: 178 . dainfo - structure containing the information 179 180 Level: beginner 181 182 Notes: 183 See DMDALocalInfo for the information that is returned 184 185 .seealso: `DMDAGetInfo()`, `DMDAGetCorners()`, `DMDALocalInfo` 186 @*/ 187 PetscErrorCode DMDAGetLocalInfo(DM da,DMDALocalInfo *info) 188 { 189 PetscInt w; 190 DM_DA *dd = (DM_DA*)da->data; 191 192 PetscFunctionBegin; 193 PetscValidHeaderSpecificType(da,DM_CLASSID,1,DMDA); 194 PetscValidPointer(info,2); 195 info->da = da; 196 info->dim = da->dim; 197 if (dd->Mo < 0) info->mx = dd->M; 198 else info->mx = dd->Mo; 199 if (dd->No < 0) info->my = dd->N; 200 else info->my = dd->No; 201 if (dd->Po < 0) info->mz = dd->P; 202 else info->mz = dd->Po; 203 info->dof = dd->w; 204 info->sw = dd->s; 205 info->bx = dd->bx; 206 info->by = dd->by; 207 info->bz = dd->bz; 208 info->st = dd->stencil_type; 209 210 /* since the xs, xe ... have all been multiplied by the number of degrees 211 of freedom per cell, w = dd->w, we divide that out before returning.*/ 212 w = dd->w; 213 info->xs = dd->xs/w + dd->xo; 214 info->xm = (dd->xe - dd->xs)/w; 215 /* the y and z have NOT been multiplied by w */ 216 info->ys = dd->ys + dd->yo; 217 info->ym = (dd->ye - dd->ys); 218 info->zs = dd->zs + dd->zo; 219 info->zm = (dd->ze - dd->zs); 220 221 info->gxs = dd->Xs/w + dd->xo; 222 info->gxm = (dd->Xe - dd->Xs)/w; 223 /* the y and z have NOT been multiplied by w */ 224 info->gys = dd->Ys + dd->yo; 225 info->gym = (dd->Ye - dd->Ys); 226 info->gzs = dd->Zs + dd->zo; 227 info->gzm = (dd->Ze - dd->Zs); 228 PetscFunctionReturn(0); 229 } 230