Actual source code: daview.c


  2: /*
  3:   Code for manipulating distributed regular arrays in parallel.
  4: */

  6: #include <petsc/private/dmdaimpl.h>

  8: #if defined(PETSC_HAVE_MATLAB)
  9:   #include <mat.h> /* MATLAB include file */

 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"};

 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

 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;

 53:   PetscFunctionBegin;
 54:   PetscCall(PetscObjectGetComm((PetscObject)da, &comm));

 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:   }

 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: }

 79: PetscErrorCode DMView_DA_VTK(DM da, PetscViewer viewer)
 80: {
 81:   Vec      coordinates;
 82:   PetscInt dim, dof, M = 0, N = 0, P = 0;

 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;

 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: }

112: /*@C
113:    DMDAGetInfo - Gets information about a given distributed array.

115:    Not Collective

117:    Input Parameter:
118: .  da - the distributed array

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`

135:    Level: beginner

137:    Note:
138:    Use NULL (NULL_INTEGER in Fortran) in place of any output parameter that is not of interest.

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;

146:   PetscFunctionBegin;
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: }

173: /*@C
174:    DMDAGetLocalInfo - Gets information about a given distributed array and this processors location in it

176:    Not Collective

178:    Input Parameter:
179: .  da - the distributed array

181:    Output Parameters:
182: .  dainfo - structure containing the information

184:    Level: beginner

186:    Note:
187:     See `DMDALocalInfo` for the information that is returned

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;

196:   PetscFunctionBegin;
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;

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);

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: }