Actual source code: ex48.c

petsc-3.11.3 2019-06-26
Report Typos and Errors
  1: static char help[] = "Test VTK structured (.vts)  and rectilinear (.vtr) viewer support with multi-dof DMDAs\n\n";

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

  6: /*
  7:   Write 3D DMDA vector with coordinates in VTK format
  8: */
  9: PetscErrorCode test_3d(const char filename[],PetscInt dof)
 10: {
 11:   MPI_Comm          comm = MPI_COMM_WORLD;
 12:   const PetscInt    M=10,N=15,P=30,sw=1;
 13:   const PetscScalar Lx=1.0,Ly=1.0,Lz=1.0;
 14:   DM                da;
 15:   Vec               v;
 16:   PetscViewer       view;
 17:   DMDALocalInfo     info;
 18:   PetscScalar       ****va;
 19:   PetscInt          i,j,k,c;
 20:   PetscErrorCode    ierr;

 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:   DMDAVecGetArrayDOF(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:         const PetscScalar x = (Lx*i)/M;
 34:         const PetscScalar y = (Ly*j)/N;
 35:         const PetscScalar z = (Lz*k)/P;
 36:         for (c=0; c<dof; ++ c) {
 37:         va[k][j][i][c] = PetscPowScalarInt(x-0.5*Lx,2)+PetscPowScalarInt(y-0.5*Ly,2)+PetscPowScalarInt(z-0.5*Lz,2) + 10.0*c;
 38:         }
 39:       }
 40:     }
 41:   }
 42:   DMDAVecRestoreArrayDOF(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: }


 52: /*
 53:   Write 2D DMDA vector with coordinates in VTK format
 54: */
 55: PetscErrorCode test_2d(const char filename[],PetscInt dof)
 56: {
 57:   MPI_Comm          comm = MPI_COMM_WORLD;
 58:   const PetscInt    M=10,N=20,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,c;
 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:   DMDAVecGetArrayDOF(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:       const PetscScalar x = (Lx*i)/M;
 78:       const PetscScalar y = (Ly*j)/N;
 79:       for (c=0; c<dof; ++c) {
 80:         va[j][i][c] = PetscPowScalarInt(x-0.5*Lx,2)+PetscPowScalarInt(y-0.5*Ly,2) + 10.0*c;
 81:       }
 82:     }
 83:   }
 84:   DMDAVecRestoreArrayDOF(da,v,&va);
 85:   PetscViewerVTKOpen(comm,filename,FILE_MODE_WRITE,&view);
 86:   VecView(v,view);
 87:   PetscViewerDestroy(&view);
 88:   VecDestroy(&v);
 89:   DMDestroy(&da);
 90:   return 0;
 91: }

 93: /*
 94:   Write a scalar and a vector field from two compatible 3d DMDAs
 95: */
 96: PetscErrorCode test_3d_compat(const char filename[],PetscInt dof)
 97: {
 98:   MPI_Comm          comm = MPI_COMM_WORLD;
 99:   const PetscInt    M=10,N=15,P=30,sw=1;
100:   const PetscScalar Lx=1.0,Ly=1.0,Lz=1.0;
101:   DM                da,daVector;
102:   Vec               v,vVector;
103:   PetscViewer       view;
104:   DMDALocalInfo     info;
105:   PetscScalar       ***va,****vVectora;
106:   PetscInt          i,j,k,c;
107:   PetscErrorCode    ierr;

109:   DMDACreate3d(comm,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DMDA_STENCIL_STAR, M,N,P,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,/* dof:*/1,sw,NULL,NULL,NULL,&da);
110:   DMSetFromOptions(da);
111:   DMSetUp(da);

113:   DMDASetUniformCoordinates(da,0.0,Lx,0.0,Ly,0.0,Lz);
114:   DMDAGetLocalInfo(da,&info);
115:   DMDACreateCompatibleDMDA(da,dof,&daVector);
116:   DMCreateGlobalVector(da,&v);
117:   DMCreateGlobalVector(daVector,&vVector);
118:   DMDAVecGetArray(da,v,&va);
119:   DMDAVecGetArrayDOF(daVector,vVector,&vVectora);
120:   for (k=info.zs; k<info.zs+info.zm; k++) {
121:     for (j=info.ys; j<info.ys+info.ym; j++) {
122:       for (i=info.xs; i<info.xs+info.xm; i++) {
123:         const PetscScalar x = (Lx*i)/M;
124:         const PetscScalar y = (Ly*j)/N;
125:         const PetscScalar z = (Lz*k)/P;
126:         va[k][j][i] = PetscPowScalarInt(x-0.5*Lx,2)+PetscPowScalarInt(y-0.5*Ly,2)+PetscPowScalarInt(z-0.5*Lz,2);
127:         for (c=0; c<dof; ++c) {
128:           vVectora[k][j][i][c] = PetscPowScalarInt(x-0.5*Lx,2)+PetscPowScalarInt(y-0.5*Ly,2)+PetscPowScalarInt(z-0.5*Lz,2) + 10.0*c;
129:         }
130:       }
131:     }
132:   }
133:   DMDAVecRestoreArray(da,v,&va);
134:   DMDAVecRestoreArrayDOF(da,v,&vVectora);
135:   PetscViewerVTKOpen(comm,filename,FILE_MODE_WRITE,&view);
136:   VecView(v,view);
137:   VecView(vVector,view);
138:   PetscViewerDestroy(&view);
139:   VecDestroy(&v);
140:   VecDestroy(&vVector);
141:   DMDestroy(&da);
142:   DMDestroy(&daVector);
143:   return 0;
144: }

146: /*
147:   Write a scalar and a vector field from two compatible 2d DMDAs
148: */
149: PetscErrorCode test_2d_compat(const char filename[],PetscInt dof)
150: {
151:   MPI_Comm          comm = MPI_COMM_WORLD;
152:   const PetscInt    M=10,N=20,sw=1;
153:   const PetscScalar Lx=1.0,Ly=1.0,Lz=1.0;
154:   DM                da, daVector;
155:   Vec               v,vVector;
156:   PetscViewer       view;
157:   DMDALocalInfo     info;
158:   PetscScalar       **va,***vVectora;
159:   PetscInt          i,j,c;
160:   PetscErrorCode    ierr;

162:   DMDACreate2d(comm,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DMDA_STENCIL_STAR, M,N,PETSC_DECIDE,PETSC_DECIDE,/* dof:*/ 1,sw,NULL,NULL,&da);
163:   DMSetFromOptions(da);
164:   DMSetUp(da);
165:   DMDASetUniformCoordinates(da,0.0,Lx,0.0,Ly,0.0,Lz);
166:   DMDACreateCompatibleDMDA(da,dof,&daVector);
167:   DMDAGetLocalInfo(da,&info);
168:   DMCreateGlobalVector(da,&v);
169:   DMCreateGlobalVector(daVector,&vVector);
170:   DMDAVecGetArray(da,v,&va);
171:   DMDAVecGetArrayDOF(daVector,vVector,&vVectora);
172:   for (j=info.ys; j<info.ys+info.ym; j++) {
173:     for (i=info.xs; i<info.xs+info.xm; i++) {
174:       const PetscScalar x = (Lx*i)/M;
175:       const PetscScalar y = (Ly*j)/N;
176:       va[j][i] = PetscPowScalarInt(x-0.5*Lx,2)+PetscPowScalarInt(y-0.5*Ly,2);
177:       for (c=0; c<dof; ++c) {
178:         vVectora[j][i][c] = PetscPowScalarInt(x-0.5*Lx,2)+PetscPowScalarInt(y-0.5*Ly,2) + 10.0*c;
179:       }
180:     }
181:   }
182:   DMDAVecRestoreArray(da,v,&va);
183:   DMDAVecRestoreArrayDOF(daVector,vVector,&vVectora);
184:   PetscViewerVTKOpen(comm,filename,FILE_MODE_WRITE,&view);
185:   VecView(v,view);
186:   VecView(vVector,view);
187:   PetscViewerDestroy(&view);
188:   VecDestroy(&v);
189:   VecDestroy(&vVector);
190:   DMDestroy(&da);
191:   DMDestroy(&daVector);
192:   return 0;
193: }

195: int main(int argc, char *argv[])
196: {
198:   PetscInt       dof;

200:   PetscInitialize(&argc,&argv,0,help);if (ierr) return ierr;
201:   dof = 2;
202:   PetscOptionsGetInt(NULL,NULL,"-dof",&dof,NULL);
203:   test_3d("3d.vtr",dof);
204:   test_2d("2d.vtr",dof);
205:   test_3d_compat("3d_compat.vtr",dof);
206:   test_2d_compat("2d_compat.vtr",dof);
207:   test_3d("3d.vts",dof);
208:   test_2d("2d.vts",dof);
209:   test_3d_compat("3d_compat.vts",dof);
210:   test_2d_compat("2d_compat.vts",dof);
211:   PetscFinalize();
212:   return ierr;
213: }

215: /*TEST

217:    build:
218:       requires: !complex

220:    test:
221:       suffix: 1
222:       nsize: 2
223:       args: -dof 2

225:    test:
226:       suffix: 2
227:       nsize: 2
228:       args: -dof 2

230:    test:
231:       suffix: 3
232:       nsize: 2
233:       args: -dof 3

235: TEST*/