Actual source code: ex48.c

  1: static char help[] = "Test VTK structured (.vts)  and rectilinear (.vtr) viewer support with multi-dof DMDAs.\n\
  2:                       Supply the -namefields flag to test with field names.\n\n";

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

  7: /* Helper function to name DMDA fields */
  8: PetscErrorCode NameFields(DM da,PetscInt dof)
  9: {
 10:   PetscInt       c;

 13:   for (c=0; c<dof; ++c) {
 14:     char fieldname[256];
 15:     PetscSNPrintf(fieldname,sizeof(fieldname),"field_%D",c);
 16:     DMDASetFieldName(da,c,fieldname);
 17:   }
 18:   return 0;
 19: }

 21: /*
 22:   Write 3D DMDA vector with coordinates in VTK format
 23: */
 24: PetscErrorCode test_3d(const char filename[],PetscInt dof,PetscBool namefields)
 25: {
 26:   MPI_Comm          comm = MPI_COMM_WORLD;
 27:   const PetscInt    M=10,N=15,P=30,sw=1;
 28:   const PetscScalar Lx=1.0,Ly=1.0,Lz=1.0;
 29:   DM                da;
 30:   Vec               v;
 31:   PetscViewer       view;
 32:   DMDALocalInfo     info;
 33:   PetscScalar       ****va;
 34:   PetscInt          i,j,k,c;

 36:   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);
 37:   DMSetFromOptions(da);
 38:   DMSetUp(da);
 39:   if (namefields) NameFields(da,dof);

 41:   DMDASetUniformCoordinates(da,0.0,Lx,0.0,Ly,0.0,Lz);
 42:   DMDAGetLocalInfo(da,&info);
 43:   DMCreateGlobalVector(da,&v);
 44:   DMDAVecGetArrayDOF(da,v,&va);
 45:   for (k=info.zs; k<info.zs+info.zm; k++) {
 46:     for (j=info.ys; j<info.ys+info.ym; j++) {
 47:       for (i=info.xs; i<info.xs+info.xm; i++) {
 48:         const PetscScalar x = (Lx*i)/M;
 49:         const PetscScalar y = (Ly*j)/N;
 50:         const PetscScalar z = (Lz*k)/P;
 51:         for (c=0; c<dof; ++ c) {
 52:         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;
 53:         }
 54:       }
 55:     }
 56:   }
 57:   DMDAVecRestoreArrayDOF(da,v,&va);
 58:   PetscViewerVTKOpen(comm,filename,FILE_MODE_WRITE,&view);
 59:   VecView(v,view);
 60:   PetscViewerDestroy(&view);
 61:   VecDestroy(&v);
 62:   DMDestroy(&da);
 63:   return 0;
 64: }

 66: /*
 67:   Write 2D DMDA vector with coordinates in VTK format
 68: */
 69: PetscErrorCode test_2d(const char filename[],PetscInt dof,PetscBool namefields)
 70: {
 71:   MPI_Comm          comm = MPI_COMM_WORLD;
 72:   const PetscInt    M=10,N=20,sw=1;
 73:   const PetscScalar Lx=1.0,Ly=1.0,Lz=1.0;
 74:   DM                da;
 75:   Vec               v;
 76:   PetscViewer       view;
 77:   DMDALocalInfo     info;
 78:   PetscScalar       ***va;
 79:   PetscInt          i,j,c;

 81:   DMDACreate2d(comm,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DMDA_STENCIL_STAR, M,N,PETSC_DECIDE,PETSC_DECIDE,dof,sw,NULL,NULL,&da);
 82:   DMSetFromOptions(da);
 83:   DMSetUp(da);
 84:   if (namefields) NameFields(da,dof);
 85:   DMDASetUniformCoordinates(da,0.0,Lx,0.0,Ly,0.0,Lz);
 86:   DMDAGetLocalInfo(da,&info);
 87:   DMCreateGlobalVector(da,&v);
 88:   DMDAVecGetArrayDOF(da,v,&va);
 89:   for (j=info.ys; j<info.ys+info.ym; j++) {
 90:     for (i=info.xs; i<info.xs+info.xm; i++) {
 91:       const PetscScalar x = (Lx*i)/M;
 92:       const PetscScalar y = (Ly*j)/N;
 93:       for (c=0; c<dof; ++c) {
 94:         va[j][i][c] = PetscPowScalarInt(x-0.5*Lx,2)+PetscPowScalarInt(y-0.5*Ly,2) + 10.0*c;
 95:       }
 96:     }
 97:   }
 98:   DMDAVecRestoreArrayDOF(da,v,&va);
 99:   PetscViewerVTKOpen(comm,filename,FILE_MODE_WRITE,&view);
100:   VecView(v,view);
101:   PetscViewerDestroy(&view);
102:   VecDestroy(&v);
103:   DMDestroy(&da);
104:   return 0;
105: }

107: /*
108:   Write a scalar and a vector field from two compatible 3d DMDAs
109: */
110: PetscErrorCode test_3d_compat(const char filename[],PetscInt dof,PetscBool namefields)
111: {
112:   MPI_Comm          comm = MPI_COMM_WORLD;
113:   const PetscInt    M=10,N=15,P=30,sw=1;
114:   const PetscScalar Lx=1.0,Ly=1.0,Lz=1.0;
115:   DM                da,daVector;
116:   Vec               v,vVector;
117:   PetscViewer       view;
118:   DMDALocalInfo     info;
119:   PetscScalar       ***va,****vVectora;
120:   PetscInt          i,j,k,c;

122:   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);
123:   DMSetFromOptions(da);
124:   DMSetUp(da);
125:   if (namefields) NameFields(da,1);

127:   DMDASetUniformCoordinates(da,0.0,Lx,0.0,Ly,0.0,Lz);
128:   DMDAGetLocalInfo(da,&info);
129:   DMDACreateCompatibleDMDA(da,dof,&daVector);
130:   if (namefields) NameFields(daVector,dof);
131:   DMCreateGlobalVector(da,&v);
132:   DMCreateGlobalVector(daVector,&vVector);
133:   DMDAVecGetArray(da,v,&va);
134:   DMDAVecGetArrayDOF(daVector,vVector,&vVectora);
135:   for (k=info.zs; k<info.zs+info.zm; k++) {
136:     for (j=info.ys; j<info.ys+info.ym; j++) {
137:       for (i=info.xs; i<info.xs+info.xm; i++) {
138:         const PetscScalar x = (Lx*i)/M;
139:         const PetscScalar y = (Ly*j)/N;
140:         const PetscScalar z = (Lz*k)/P;
141:         va[k][j][i] = PetscPowScalarInt(x-0.5*Lx,2)+PetscPowScalarInt(y-0.5*Ly,2)+PetscPowScalarInt(z-0.5*Lz,2);
142:         for (c=0; c<dof; ++c) {
143:           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;
144:         }
145:       }
146:     }
147:   }
148:   DMDAVecRestoreArray(da,v,&va);
149:   DMDAVecRestoreArrayDOF(da,v,&vVectora);
150:   PetscViewerVTKOpen(comm,filename,FILE_MODE_WRITE,&view);
151:   VecView(v,view);
152:   VecView(vVector,view);
153:   PetscViewerDestroy(&view);
154:   VecDestroy(&v);
155:   VecDestroy(&vVector);
156:   DMDestroy(&da);
157:   DMDestroy(&daVector);
158:   return 0;
159: }

161: /*
162:   Write a scalar and a vector field from two compatible 2d DMDAs
163: */
164: PetscErrorCode test_2d_compat(const char filename[],PetscInt dof,PetscBool namefields)
165: {
166:   MPI_Comm          comm = MPI_COMM_WORLD;
167:   const PetscInt    M=10,N=20,sw=1;
168:   const PetscScalar Lx=1.0,Ly=1.0,Lz=1.0;
169:   DM                da, daVector;
170:   Vec               v,vVector;
171:   PetscViewer       view;
172:   DMDALocalInfo     info;
173:   PetscScalar       **va,***vVectora;
174:   PetscInt          i,j,c;

176:   DMDACreate2d(comm,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DMDA_STENCIL_STAR, M,N,PETSC_DECIDE,PETSC_DECIDE,/* dof:*/ 1,sw,NULL,NULL,&da);
177:   DMSetFromOptions(da);
178:   DMSetUp(da);
179:   if (namefields) NameFields(da,1);
180:   DMDASetUniformCoordinates(da,0.0,Lx,0.0,Ly,0.0,Lz);
181:   DMDACreateCompatibleDMDA(da,dof,&daVector);
182:   if (namefields) NameFields(daVector,dof);
183:   DMDAGetLocalInfo(da,&info);
184:   DMCreateGlobalVector(da,&v);
185:   DMCreateGlobalVector(daVector,&vVector);
186:   DMDAVecGetArray(da,v,&va);
187:   DMDAVecGetArrayDOF(daVector,vVector,&vVectora);
188:   for (j=info.ys; j<info.ys+info.ym; j++) {
189:     for (i=info.xs; i<info.xs+info.xm; i++) {
190:       const PetscScalar x = (Lx*i)/M;
191:       const PetscScalar y = (Ly*j)/N;
192:       va[j][i] = PetscPowScalarInt(x-0.5*Lx,2)+PetscPowScalarInt(y-0.5*Ly,2);
193:       for (c=0; c<dof; ++c) {
194:         vVectora[j][i][c] = PetscPowScalarInt(x-0.5*Lx,2)+PetscPowScalarInt(y-0.5*Ly,2) + 10.0*c;
195:       }
196:     }
197:   }
198:   DMDAVecRestoreArray(da,v,&va);
199:   DMDAVecRestoreArrayDOF(daVector,vVector,&vVectora);
200:   PetscViewerVTKOpen(comm,filename,FILE_MODE_WRITE,&view);
201:   VecView(v,view);
202:   VecView(vVector,view);
203:   PetscViewerDestroy(&view);
204:   VecDestroy(&v);
205:   VecDestroy(&vVector);
206:   DMDestroy(&da);
207:   DMDestroy(&daVector);
208:   return 0;
209: }

211: int main(int argc, char *argv[])
212: {
213:   PetscInt       dof;
214:   PetscBool      namefields;

216:   PetscInitialize(&argc,&argv,0,help);
217:   dof = 2;
218:   PetscOptionsGetInt(NULL,NULL,"-dof",&dof,NULL);
219:   namefields = PETSC_FALSE;
220:   PetscOptionsGetBool(NULL,NULL,"-namefields",&namefields,NULL);
221:   test_3d("3d.vtr",dof,namefields);
222:   test_2d("2d.vtr",dof,namefields);
223:   test_3d_compat("3d_compat.vtr",dof,namefields);
224:   test_2d_compat("2d_compat.vtr",dof,namefields);
225:   test_3d("3d.vts",dof,namefields);
226:   test_2d("2d.vts",dof,namefields);
227:   test_3d_compat("3d_compat.vts",dof,namefields);
228:   test_2d_compat("2d_compat.vts",dof,namefields);
229:   PetscFinalize();
230:   return 0;
231: }

233: /*TEST

235:    build:
236:       requires: !complex

238:    test:
239:       suffix: 1
240:       nsize: 2
241:       args: -dof 2

243:    test:
244:       suffix: 2
245:       nsize: 2
246:       args: -dof 2

248:    test:
249:       suffix: 3
250:       nsize: 2
251:       args: -dof 3

253:    test:
254:       suffix: 4
255:       nsize: 1
256:       args: -dof 2 -namefields

258:    test:
259:       suffix: 5
260:       nsize: 2
261:       args: -dof 4 -namefields

263: TEST*/