Actual source code: ex42.c

  1: /* -*- Mode: C++; c-basic-offset:2 ; indent-tabs-mode:nil ; -*- */

  3: static char help[] = "Test VTK Rectilinear grid (.vtr) viewer support\n\n";

  5: #include <petscdm.h>
  6: #include <petscdmda.h>

  8: /*
  9:   Write 3D DMDA vector with coordinates in VTK VTR format

 11: */
 12: PetscErrorCode test_3d(const char filename[])
 13: {
 14:   MPI_Comm          comm = MPI_COMM_WORLD;
 15:   const PetscInt    M=10, N=15, P=30, dof=1, sw=1;
 16:   const PetscScalar Lx=1.0, Ly=1.0, Lz=1.0;
 17:   DM                da;
 18:   Vec               v;
 19:   PetscViewer       view;
 20:   DMDALocalInfo     info;
 21:   PetscScalar       ***va;
 22:   PetscInt          i,j,k;
 23:   PetscErrorCode    ierr;

 25:   DMDACreate3d(comm,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DMDA_STENCIL_STAR, M,N,P,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,dof,sw,NULL,NULL,NULL,&da);
 26:   DMSetFromOptions(da);
 27:   DMSetUp(da);

 29:   DMDASetUniformCoordinates(da,0.0,Lx,0.0,Ly,0.0,Lz);
 30:   DMDAGetLocalInfo(da,&info);
 31:   DMCreateGlobalVector(da,&v);
 32:   DMDAVecGetArray(da,v,&va);
 33:   for (k=info.zs; k<info.zs+info.zm; k++) {
 34:     for (j=info.ys; j<info.ys+info.ym; j++) {
 35:       for (i=info.xs; i<info.xs+info.xm; i++) {
 36:         PetscScalar x = (Lx*i)/M;
 37:         PetscScalar y = (Ly*j)/N;
 38:         PetscScalar z = (Lz*k)/P;
 39:         va[k][j][i] = PetscPowScalarInt(x-0.5*Lx,2)+PetscPowScalarInt(y-0.5*Ly,2)+PetscPowScalarInt(z-0.5*Lz,2);
 40:       }
 41:     }
 42:   }
 43:   DMDAVecRestoreArray(da,v,&va);
 44:   PetscViewerVTKOpen(comm,filename,FILE_MODE_WRITE,&view);
 45:   VecView(v,view);
 46:   PetscViewerDestroy(&view);
 47:   VecDestroy(&v);
 48:   DMDestroy(&da);
 49:   return 0;
 50: }

 52: /*
 53:   Write 2D DMDA vector with coordinates in VTK VTR format

 55: */
 56: PetscErrorCode test_2d(const char filename[])
 57: {
 58:   MPI_Comm          comm = MPI_COMM_WORLD;
 59:   const PetscInt    M=10, N=20, dof=1, sw=1;
 60:   const PetscScalar Lx=1.0, Ly=1.0, Lz=1.0;
 61:   DM                da;
 62:   Vec               v;
 63:   PetscViewer       view;
 64:   DMDALocalInfo     info;
 65:   PetscScalar       **va;
 66:   PetscInt          i,j;
 67:   PetscErrorCode    ierr;

 69:   DMDACreate2d(comm,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DMDA_STENCIL_STAR, M,N,PETSC_DECIDE,PETSC_DECIDE,dof,sw,NULL,NULL,&da);
 70:   DMSetFromOptions(da);
 71:   DMSetUp(da);
 72:   DMDASetUniformCoordinates(da,0.0,Lx,0.0,Ly,0.0,Lz);
 73:   DMDAGetLocalInfo(da,&info);
 74:   DMCreateGlobalVector(da,&v);
 75:   DMDAVecGetArray(da,v,&va);
 76:   for (j=info.ys; j<info.ys+info.ym; j++) {
 77:     for (i=info.xs; i<info.xs+info.xm; i++) {
 78:       PetscScalar x = (Lx*i)/M;
 79:       PetscScalar y = (Ly*j)/N;
 80:       va[j][i] = PetscPowScalarInt(x-0.5*Lx,2)+PetscPowScalarInt(y-0.5*Ly,2);
 81:     }
 82:   }
 83:   DMDAVecRestoreArray(da,v,&va);
 84:   PetscViewerVTKOpen(comm,filename,FILE_MODE_WRITE,&view);
 85:   VecView(v,view);
 86:   PetscViewerDestroy(&view);
 87:   VecDestroy(&v);
 88:   DMDestroy(&da);
 89:   return 0;
 90: }

 92: /*
 93:   Write 2D DMDA vector without coordinates in VTK VTR format

 95: */
 96: PetscErrorCode test_2d_nocoord(const char filename[])
 97: {
 98:   MPI_Comm          comm = MPI_COMM_WORLD;
 99:   const PetscInt    M=10, N=20, dof=1, sw=1;
100:   const PetscScalar Lx=1.0, Ly=1.0;
101:   DM                da;
102:   Vec               v;
103:   PetscViewer       view;
104:   DMDALocalInfo     info;
105:   PetscScalar       **va;
106:   PetscInt          i,j;
107:   PetscErrorCode    ierr;

109:   DMDACreate2d(comm,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DMDA_STENCIL_STAR,M,N,PETSC_DECIDE,PETSC_DECIDE,dof,sw,NULL,NULL,&da);
110:   DMSetFromOptions(da);
111:   DMSetUp(da);
112:   DMDAGetLocalInfo(da,&info);
113:   DMCreateGlobalVector(da,&v);
114:   DMDAVecGetArray(da,v,&va);
115:   for (j=info.ys; j<info.ys+info.ym; j++) {
116:     for (i=info.xs; i<info.xs+info.xm; i++) {
117:       PetscScalar x = (Lx*i)/M;
118:       PetscScalar y = (Ly*j)/N;
119:       va[j][i] = PetscPowScalarInt(x-0.5*Lx,2)+PetscPowScalarInt(y-0.5*Ly,2);
120:     }
121:   }
122:   DMDAVecRestoreArray(da,v,&va);
123:   PetscViewerVTKOpen(comm,filename,FILE_MODE_WRITE,&view);
124:   VecView(v,view);
125:   PetscViewerDestroy(&view);
126:   VecDestroy(&v);
127:   DMDestroy(&da);
128:   return 0;
129: }

131: /*
132:   Write 3D DMDA vector without coordinates in VTK VTR format

134: */
135: PetscErrorCode test_3d_nocoord(const char filename[])
136: {
137:   MPI_Comm          comm = MPI_COMM_WORLD;
138:   const PetscInt    M=10, N=20, P=30, dof=1, sw=1;
139:   const PetscScalar Lx=1.0, Ly=1.0, Lz=1.0;
140:   DM                da;
141:   Vec               v;
142:   PetscViewer       view;
143:   DMDALocalInfo     info;
144:   PetscScalar       ***va;
145:   PetscInt          i,j,k;
146:   PetscErrorCode    ierr;

148:   DMDACreate3d(comm,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DMDA_STENCIL_STAR, M,N,P,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,dof,sw,NULL,NULL,NULL,&da);
149:   DMSetFromOptions(da);
150:   DMSetUp(da);

152:   DMDAGetLocalInfo(da,&info);
153:   DMCreateGlobalVector(da,&v);
154:   DMDAVecGetArray(da,v,&va);
155:   for (k=info.zs; k<info.zs+info.zm; k++) {
156:     for (j=info.ys; j<info.ys+info.ym; j++) {
157:       for (i=info.xs; i<info.xs+info.xm; i++) {
158:         PetscScalar x = (Lx*i)/M;
159:         PetscScalar y = (Ly*j)/N;
160:         PetscScalar z = (Lz*k)/P;
161:         va[k][j][i] = PetscPowScalarInt(x-0.5*Lx,2)+PetscPowScalarInt(y-0.5*Ly,2)+PetscPowScalarInt(z-0.5*Lz,2);
162:       }
163:     }
164:   }
165:   DMDAVecRestoreArray(da,v,&va);
166:   PetscViewerVTKOpen(comm,filename,FILE_MODE_WRITE,&view);
167:   VecView(v,view);
168:   PetscViewerDestroy(&view);
169:   VecDestroy(&v);
170:   DMDestroy(&da);
171:   return 0;
172: }

174: int main(int argc, char *argv[])
175: {

178:   PetscInitialize(&argc,&argv,0,help);if (ierr) return ierr;
179:   test_3d("3d.vtr");
180:   test_2d("2d.vtr");
181:   test_2d_nocoord("2d_nocoord.vtr");
182:   test_3d_nocoord("3d_nocoord.vtr");
183:   PetscFinalize();
184:   return ierr;
185: }

187: /*TEST

189:    build:
190:       requires: !complex

192:    test:
193:       nsize: 2

195: TEST*/