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

  6: #include <petsc-private/daimpl.h>    /*I   "petscdmda.h"   I*/

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

 13: PetscErrorCode DMView_DA_Matlab(DM da,PetscViewer viewer)
 14: {
 15:   PetscErrorCode   ierr;
 16:   PetscMPIInt      rank;
 17:   PetscInt         dim,m,n,p,dof,swidth;
 18:   DMDAStencilType  stencil;
 19:   DMDABoundaryType bx,by,bz;
 20:   mxArray          *mx;
 21:   const char       *fnames[] = {"dimension","m","n","p","dof","stencil_width","bx","by","bz","stencil_type"};

 24:   MPI_Comm_rank(((PetscObject)da)->comm,&rank);
 25:   if (!rank) {
 26:     DMDAGetInfo(da,&dim,&m,&n,&p,0,0,0,&dof,&swidth,&bx,&by,&bz,&stencil);
 27:     mx = mxCreateStructMatrix(1,1,8,(const char **)fnames);
 28:     if (!mx) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Unable to generate MATLAB struct array to hold DMDA informations");
 29:     mxSetFieldByNumber(mx,0,0,mxCreateDoubleScalar((double)dim));
 30:     mxSetFieldByNumber(mx,0,1,mxCreateDoubleScalar((double)m));
 31:     mxSetFieldByNumber(mx,0,2,mxCreateDoubleScalar((double)n));
 32:     mxSetFieldByNumber(mx,0,3,mxCreateDoubleScalar((double)p));
 33:     mxSetFieldByNumber(mx,0,4,mxCreateDoubleScalar((double)dof));
 34:     mxSetFieldByNumber(mx,0,5,mxCreateDoubleScalar((double)swidth));
 35:     mxSetFieldByNumber(mx,0,6,mxCreateDoubleScalar((double)bx));
 36:     mxSetFieldByNumber(mx,0,7,mxCreateDoubleScalar((double)by));
 37:     mxSetFieldByNumber(mx,0,8,mxCreateDoubleScalar((double)bz));
 38:     mxSetFieldByNumber(mx,0,9,mxCreateDoubleScalar((double)stencil));
 39:     PetscObjectName((PetscObject)da);
 40:     PetscViewerMatlabPutVariable(viewer,((PetscObject)da)->name,mx);
 41:   }
 42:   return(0);
 43: }
 44: #endif

 48: PetscErrorCode DMView_DA_Binary(DM da,PetscViewer viewer)
 49: {
 50:   PetscErrorCode   ierr;
 51:   PetscMPIInt      rank;
 52:   PetscInt         dim,m,n,p,dof,swidth,M,N,P;
 53:   DMDAStencilType  stencil;
 54:   DMDABoundaryType bx,by,bz;
 55:   MPI_Comm         comm;
 56:   DM_DA            *dd = (DM_DA*)da->data;
 57:   PetscInt         classid = DM_FILE_CLASSID,subclassid = DMDA_FILE_CLASSID ;
 58:   PetscBool        coors = PETSC_FALSE;

 61:   PetscObjectGetComm((PetscObject)da,&comm);

 63:   DMDAGetInfo(da,&dim,&m,&n,&p,&M,&N,&P,&dof,&swidth,&bx,&by,&bz,&stencil);
 64:   MPI_Comm_rank(comm,&rank);
 65:   if (!rank) {

 67:     PetscViewerBinaryWrite(viewer,&classid,1,PETSC_INT,PETSC_FALSE);
 68:     PetscViewerBinaryWrite(viewer,&subclassid,1,PETSC_INT,PETSC_FALSE);
 69:     PetscViewerBinaryWrite(viewer,&dim,1,PETSC_INT,PETSC_FALSE);
 70:     PetscViewerBinaryWrite(viewer,&m,1,PETSC_INT,PETSC_FALSE);
 71:     PetscViewerBinaryWrite(viewer,&n,1,PETSC_INT,PETSC_FALSE);
 72:     PetscViewerBinaryWrite(viewer,&p,1,PETSC_INT,PETSC_FALSE);
 73:     PetscViewerBinaryWrite(viewer,&dof,1,PETSC_INT,PETSC_FALSE);
 74:     PetscViewerBinaryWrite(viewer,&swidth,1,PETSC_INT,PETSC_FALSE);
 75:     PetscViewerBinaryWrite(viewer,&bx,1,PETSC_ENUM,PETSC_FALSE);
 76:     PetscViewerBinaryWrite(viewer,&by,1,PETSC_ENUM,PETSC_FALSE);
 77:     PetscViewerBinaryWrite(viewer,&bz,1,PETSC_ENUM,PETSC_FALSE);
 78:     PetscViewerBinaryWrite(viewer,&stencil,1,PETSC_ENUM,PETSC_FALSE);
 79:     if (dd->coordinates) coors = PETSC_TRUE;
 80:     PetscViewerBinaryWrite(viewer,&coors,1,PETSC_BOOL,PETSC_FALSE);
 81:   }

 83:   /* save the coordinates if they exist to disk (in the natural ordering) */
 84:   if (dd->coordinates) {
 85:     VecView(dd->coordinates,viewer);
 86:   }
 87:   return(0);
 88: }

 92: PetscErrorCode DMView_DA_VTK(DM da, PetscViewer viewer)
 93: {
 94:   PetscInt       dim, dof, M = 0, N = 0, P = 0;
 96:   DM_DA          *dd = (DM_DA*)da->data;

 99:   DMDAGetInfo(da, &dim, &M, &N, &P, PETSC_NULL, PETSC_NULL, PETSC_NULL, &dof, PETSC_NULL, PETSC_NULL, PETSC_NULL, PETSC_NULL, PETSC_NULL);
100:   /* if (dim != 3) {SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP, "VTK output only works for three dimensional DMDAs.");} */
101:   if (!dd->coordinates) SETERRQ(((PetscObject)da)->comm,PETSC_ERR_SUP, "VTK output requires DMDA coordinates.");
102:   /* Write Header */
103:   PetscViewerASCIIPrintf(viewer,"# vtk DataFile Version 2.0\n");
104:   PetscViewerASCIIPrintf(viewer,"Structured Mesh Example\n");
105:   PetscViewerASCIIPrintf(viewer,"ASCII\n");
106:   PetscViewerASCIIPrintf(viewer,"DATASET STRUCTURED_GRID\n");
107:   PetscViewerASCIIPrintf(viewer,"DIMENSIONS %d %d %d\n", M, N, P);
108:   PetscViewerASCIIPrintf(viewer,"POINTS %d double\n", M*N*P);
109:   if (dd->coordinates) {
110:     DM  dac;
111:     Vec natural;

113:     DMDAGetCoordinateDA(da, &dac);
114:     DMDACreateNaturalVector(dac, &natural);
115:     PetscObjectSetOptionsPrefix((PetscObject) natural, "coor_");
116:     DMDAGlobalToNaturalBegin(dac, dd->coordinates, INSERT_VALUES, natural);
117:     DMDAGlobalToNaturalEnd(dac, dd->coordinates, INSERT_VALUES, natural);
118:     PetscViewerPushFormat(viewer, PETSC_VIEWER_ASCII_VTK_COORDS);
119:     VecView(natural, viewer);
120:     PetscViewerPopFormat(viewer);
121:     VecDestroy(&natural);
122:   }
123:   return(0);
124: }

128: /*@C
129:    DMDAGetInfo - Gets information about a given distributed array.

131:    Not Collective

133:    Input Parameter:
134: .  da - the distributed array

136:    Output Parameters:
137: +  dim      - dimension of the distributed array (1, 2, or 3)
138: .  M, N, P  - global dimension in each direction of the array
139: .  m, n, p  - corresponding number of procs in each dimension
140: .  dof      - number of degrees of freedom per node
141: .  s        - stencil width
142: .  bx,by,bz - type of ghost nodes at boundary, one of DMDA_BOUNDARY_NONE, DMDA_BOUNDARY_GHOSTED, 
143:               DMDA_BOUNDARY_MIRROR, DMDA_BOUNDARY_PERIODIC
144: -  st       - stencil type, either DMDA_STENCIL_STAR or DMDA_STENCIL_BOX

146:    Level: beginner
147:   
148:    Note:
149:    Use PETSC_NULL (PETSC_NULL_INTEGER in Fortran) in place of any output parameter that is not of interest.

151: .keywords: distributed array, get, information

153: .seealso: DMView(), DMDAGetCorners(), DMDAGetLocalInfo()
154: @*/
155: PetscErrorCode  DMDAGetInfo(DM da,PetscInt *dim,PetscInt *M,PetscInt *N,PetscInt *P,PetscInt *m,PetscInt *n,PetscInt *p,PetscInt *dof,PetscInt *s,DMDABoundaryType *bx,DMDABoundaryType *by,DMDABoundaryType *bz,DMDAStencilType *st)
156: {
157:   DM_DA *dd = (DM_DA*)da->data;

161:   if (dim)  *dim  = dd->dim;
162:   if (M)    *M    = dd->M;
163:   if (N)    *N    = dd->N;
164:   if (P)    *P    = dd->P;
165:   if (m)    *m    = dd->m;
166:   if (n)    *n    = dd->n;
167:   if (p)    *p    = dd->p;
168:   if (dof)  *dof  = dd->w;
169:   if (s)    *s    = dd->s;
170:   if (bx) *bx = dd->bx;
171:   if (by) *by = dd->by;
172:   if (bz) *bz = dd->bz;
173:   if (st)   *st   = dd->stencil_type;
174:   return(0);
175: }

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

182:    Not Collective

184:    Input Parameter:
185: .  da - the distributed array

187:    Output Parameters:
188: .  dainfo - structure containing the information

190:    Level: beginner
191:   
192: .keywords: distributed array, get, information

194: .seealso: DMDAGetInfo(), DMDAGetCorners()
195: @*/
196: PetscErrorCode  DMDAGetLocalInfo(DM da,DMDALocalInfo *info)
197: {
198:   PetscInt w;
199:   DM_DA    *dd = (DM_DA*)da->data;

204:   info->da   = da;
205:   info->dim  = dd->dim;
206:   info->mx   = dd->M;
207:   info->my   = dd->N;
208:   info->mz   = dd->P;
209:   info->dof  = dd->w;
210:   info->sw   = dd->s;
211:   info->bx   = dd->bx;
212:   info->by   = dd->by;
213:   info->bz   = dd->bz;
214:   info->st   = dd->stencil_type;

216:   /* since the xs, xe ... have all been multiplied by the number of degrees 
217:      of freedom per cell, w = dd->w, we divide that out before returning.*/
218:   w = dd->w;
219:   info->xs = dd->xs/w;
220:   info->xm = (dd->xe - dd->xs)/w;
221:   /* the y and z have NOT been multiplied by w */
222:   info->ys = dd->ys;
223:   info->ym = (dd->ye - dd->ys);
224:   info->zs = dd->zs;
225:   info->zm = (dd->ze - dd->zs);

227:   info->gxs = dd->Xs/w;
228:   info->gxm = (dd->Xe - dd->Xs)/w;
229:   /* the y and z have NOT been multiplied by w */
230:   info->gys = dd->Ys;
231:   info->gym = (dd->Ye - dd->Ys);
232:   info->gzs = dd->Zs;
233:   info->gzm = (dd->Ze - dd->Zs);
234:   return(0);
235: }