Actual source code: ex47.c

  1: static char help[] = "Test VTK structured grid (.vts) viewer support\n\n";

  3: #include <petscdm.h>
  4: #include <petscdmda.h>

  6: /*
  7:   Write 3D DMDA vector with coordinates in VTK .vts format

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

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

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

 49: /*
 50:   Write 2D DMDA vector with coordinates in VTK .vts format

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

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

 88: /*
 89:   Write 2D DMDA vector without coordinates in VTK .vts format

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

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

126: /*
127:   Write 3D DMDA vector without coordinates in VTK .vts format

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

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

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

168: int main(int argc, char *argv[])
169: {

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

180: /*TEST

182:    build:
183:       requires: !complex

185:    test:
186:       nsize: 2

188: TEST*/