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