Actual source code: ex47.c

petsc-3.12.5 2020-03-29
Report Typos and Errors
  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;
 21:   PetscErrorCode    ierr;

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

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


 51: /*
 52:   Write 2D DMDA vector with coordinates in VTK .vts 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;
 66:   PetscErrorCode    ierr;

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


 92: /*
 93:   Write 2D DMDA vector without coordinates in VTK .vts 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: }


132: /*
133:   Write 3D DMDA vector without coordinates in VTK .vts format

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

149:   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);
150:   DMSetFromOptions(da);
151:   DMSetUp(da);

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

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

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


189: /*TEST

191:    build:
192:       requires: !complex

194:    test:
195:       nsize: 2

197: TEST*/