Actual source code: daview.c

petsc-master 2017-09-24
Report Typos and Errors

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

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

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

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

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

 44: PetscErrorCode DMView_DA_Binary(DM da,PetscViewer viewer)
 45: {
 46:   PetscErrorCode   ierr;
 47:   PetscMPIInt      rank;
 48:   PetscInt         dim,m,n,p,dof,swidth,M,N,P;
 49:   DMDAStencilType  stencil;
 50:   DMBoundaryType   bx,by,bz;
 51:   MPI_Comm         comm;
 52:   PetscBool        coors = PETSC_FALSE;

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

 57:   DMDAGetInfo(da,&dim,&m,&n,&p,&M,&N,&P,&dof,&swidth,&bx,&by,&bz,&stencil);
 58:   MPI_Comm_rank(comm,&rank);
 59:   if (!rank) {
 60:     PetscViewerBinaryWrite(viewer,&dim,1,PETSC_INT,PETSC_FALSE);
 61:     PetscViewerBinaryWrite(viewer,&m,1,PETSC_INT,PETSC_FALSE);
 62:     PetscViewerBinaryWrite(viewer,&n,1,PETSC_INT,PETSC_FALSE);
 63:     PetscViewerBinaryWrite(viewer,&p,1,PETSC_INT,PETSC_FALSE);
 64:     PetscViewerBinaryWrite(viewer,&dof,1,PETSC_INT,PETSC_FALSE);
 65:     PetscViewerBinaryWrite(viewer,&swidth,1,PETSC_INT,PETSC_FALSE);
 66:     PetscViewerBinaryWrite(viewer,&bx,1,PETSC_ENUM,PETSC_FALSE);
 67:     PetscViewerBinaryWrite(viewer,&by,1,PETSC_ENUM,PETSC_FALSE);
 68:     PetscViewerBinaryWrite(viewer,&bz,1,PETSC_ENUM,PETSC_FALSE);
 69:     PetscViewerBinaryWrite(viewer,&stencil,1,PETSC_ENUM,PETSC_FALSE);
 70:     if (da->coordinates) coors = PETSC_TRUE;
 71:     PetscViewerBinaryWrite(viewer,&coors,1,PETSC_BOOL,PETSC_FALSE);
 72:   }

 74:   /* save the coordinates if they exist to disk (in the natural ordering) */
 75:   if (da->coordinates) {
 76:     VecView(da->coordinates,viewer);
 77:   }
 78:   return(0);
 79: }

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

 87:   DMDAGetInfo(da, &dim, &M, &N, &P, NULL, NULL, NULL, &dof, NULL, NULL, NULL, NULL, NULL);
 88:   if (!da->coordinates) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP, "VTK output requires DMDA coordinates.");
 89:   /* Write Header */
 90:   PetscViewerASCIIPrintf(viewer,"# vtk DataFile Version 2.0\n");
 91:   PetscViewerASCIIPrintf(viewer,"Structured Mesh Example\n");
 92:   PetscViewerASCIIPrintf(viewer,"ASCII\n");
 93:   PetscViewerASCIIPrintf(viewer,"DATASET STRUCTURED_GRID\n");
 94:   PetscViewerASCIIPrintf(viewer,"DIMENSIONS %d %d %d\n", M, N, P);
 95:   PetscViewerASCIIPrintf(viewer,"POINTS %d double\n", M*N*P);
 96:   if (da->coordinates) {
 97:     DM  dac;
 98:     Vec natural;

100:     DMGetCoordinateDM(da, &dac);
101:     DMDACreateNaturalVector(dac, &natural);
102:     PetscObjectSetOptionsPrefix((PetscObject) natural, "coor_");
103:     DMDAGlobalToNaturalBegin(dac, da->coordinates, INSERT_VALUES, natural);
104:     DMDAGlobalToNaturalEnd(dac, da->coordinates, INSERT_VALUES, natural);
105:     PetscViewerPushFormat(viewer, PETSC_VIEWER_ASCII_VTK_COORDS);
106:     VecView(natural, viewer);
107:     PetscViewerPopFormat(viewer);
108:     VecDestroy(&natural);
109:   }
110:   return(0);
111: }

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

116:    Not Collective

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

121:    Output Parameters:
122: +  dim      - dimension of the distributed array (1, 2, or 3)
123: .  M, N, P  - global dimension in each direction of the array
124: .  m, n, p  - corresponding number of procs in each dimension
125: .  dof      - number of degrees of freedom per node
126: .  s        - stencil width
127: .  bx,by,bz - type of ghost nodes at boundary, one of DM_BOUNDARY_NONE, DM_BOUNDARY_GHOSTED,
128:               DM_BOUNDARY_MIRROR, DM_BOUNDARY_PERIODIC
129: -  st       - stencil type, either DMDA_STENCIL_STAR or DMDA_STENCIL_BOX

131:    Level: beginner

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

136: .keywords: distributed array, get, information

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;

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:   return(0);
169: }

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

174:    Not Collective

176:    Input Parameter:
177: .  da - the distributed array

179:    Output Parameters:
180: .  dainfo - structure containing the information

182:    Level: beginner

184:    Notes: See DMDALocalInfo for the information that is returned

186:    Fortran Notes: In Fortran the routine is DMDAGetLocalInfoF90(), see DMDALocalInfo for how to access the values

188:    Developer Notes: Not sure why the Fortran function has a F90() in the name since it does not utilize F90 constructs.

190: .keywords: distributed array, get, information

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

202:   info->da  = da;
203:   info->dim = da->dim;
204:   if (dd->Mo < 0) info->mx = dd->M;
205:   else info->mx = dd->Mo;
206:   if (dd->No < 0) info->my = dd->N;
207:   else info->my = dd->No;
208:   if (dd->Po < 0) info->mz = dd->P;
209:   else info->mz = dd->Po;
210:   info->dof = dd->w;
211:   info->sw  = dd->s;
212:   info->bx  = dd->bx;
213:   info->by  = dd->by;
214:   info->bz  = dd->bz;
215:   info->st  = dd->stencil_type;

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

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