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;

 24:   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);
 25:   DMSetFromOptions(da);
 26:   DMSetUp(da);

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

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

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

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

 90: /*
 91:   Write 2D DMDA vector without coordinates in VTK VTR format

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

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

128: /*
129:   Write 3D DMDA vector without coordinates in VTK VTR format

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

144:   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);
145:   DMSetFromOptions(da);
146:   DMSetUp(da);

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

170: int main(int argc, char *argv[])
171: {

173:   PetscInitialize(&argc,&argv,0,help);
174:   test_3d("3d.vtr");
175:   test_2d("2d.vtr");
176:   test_2d_nocoord("2d_nocoord.vtr");
177:   test_3d_nocoord("3d_nocoord.vtr");
178:   PetscFinalize();
179:   return 0;
180: }

182: /*TEST

184:    build:
185:       requires: !complex

187:    test:
188:       nsize: 2

190: TEST*/