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: {
 11:   PetscInt       c;

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

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

 38:   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);
 39:   DMSetFromOptions(da);
 40:   DMSetUp(da);
 41:   if (namefields) {NameFields(da,dof);}

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

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

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

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

126:   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);
127:   DMSetFromOptions(da);
128:   DMSetUp(da);
129:   if (namefields) {NameFields(da,1);}

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

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

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

216: int main(int argc, char *argv[])
217: {
219:   PetscInt       dof;
220:   PetscBool      namefields;

222:   PetscInitialize(&argc,&argv,0,help);if (ierr) return ierr;
223:   dof = 2;
224:   PetscOptionsGetInt(NULL,NULL,"-dof",&dof,NULL);
225:   namefields = PETSC_FALSE;
226:   PetscOptionsGetBool(NULL,NULL,"-namefields",&namefields,NULL);
227:   test_3d("3d.vtr",dof,namefields);
228:   test_2d("2d.vtr",dof,namefields);
229:   test_3d_compat("3d_compat.vtr",dof,namefields);
230:   test_2d_compat("2d_compat.vtr",dof,namefields);
231:   test_3d("3d.vts",dof,namefields);
232:   test_2d("2d.vts",dof,namefields);
233:   test_3d_compat("3d_compat.vts",dof,namefields);
234:   test_2d_compat("2d_compat.vts",dof,namefields);
235:   PetscFinalize();
236:   return ierr;
237: }

239: /*TEST

241:    build:
242:       requires: !complex

244:    test:
245:       suffix: 1
246:       nsize: 2
247:       args: -dof 2

249:    test:
250:       suffix: 2
251:       nsize: 2
252:       args: -dof 2

254:    test:
255:       suffix: 3
256:       nsize: 2
257:       args: -dof 3

259:    test:
260:       suffix: 4
261:       nsize: 1
262:       args: -dof 2 -namefields

264:    test:
265:       suffix: 5
266:       nsize: 2
267:       args: -dof 4 -namefields

269: TEST*/