xref: /petsc/src/dm/tests/ex48.c (revision 7d5fd1e4d9337468ad3f05b65b7facdcd2dfd2a4) 
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";
3 
4 #include <petscdm.h>
5 #include <petscdmda.h>
6 
7 /* Helper function to name DMDA fields */
8 PetscErrorCode NameFields(DM da,PetscInt dof)
9 {
10   PetscErrorCode ierr;
11   PetscInt       c;
12 
13   PetscFunctionBeginUser;
14   for (c=0; c<dof; ++c) {
15     char fieldname[256];
16     ierr = PetscSNPrintf(fieldname,sizeof(fieldname),"field_%D",c);CHKERRQ(ierr);
17     ierr = DMDASetFieldName(da,c,fieldname);CHKERRQ(ierr);
18   }
19   PetscFunctionReturn(0);
20 }
21 
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;
37 
38   ierr = 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);CHKERRQ(ierr);
39   ierr = DMSetFromOptions(da);CHKERRQ(ierr);
40   ierr = DMSetUp(da);CHKERRQ(ierr);
41   if (namefields) {ierr = NameFields(da,dof);CHKERRQ(ierr);}
42 
43   ierr = DMDASetUniformCoordinates(da,0.0,Lx,0.0,Ly,0.0,Lz);CHKERRQ(ierr);
44   ierr = DMDAGetLocalInfo(da,&info);CHKERRQ(ierr);
45   ierr = DMCreateGlobalVector(da,&v);CHKERRQ(ierr);
46   ierr = DMDAVecGetArrayDOF(da,v,&va);CHKERRQ(ierr);
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   ierr = DMDAVecRestoreArrayDOF(da,v,&va);CHKERRQ(ierr);
60   ierr = PetscViewerVTKOpen(comm,filename,FILE_MODE_WRITE,&view);CHKERRQ(ierr);
61   ierr = VecView(v,view);CHKERRQ(ierr);
62   ierr = PetscViewerDestroy(&view);CHKERRQ(ierr);
63   ierr = VecDestroy(&v);CHKERRQ(ierr);
64   ierr = DMDestroy(&da);CHKERRQ(ierr);
65   return 0;
66 }
67 
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;
83 
84   ierr = DMDACreate2d(comm,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DMDA_STENCIL_STAR, M,N,PETSC_DECIDE,PETSC_DECIDE,dof,sw,NULL,NULL,&da);CHKERRQ(ierr);
85   ierr = DMSetFromOptions(da);CHKERRQ(ierr);
86   ierr = DMSetUp(da);CHKERRQ(ierr);
87   if (namefields) {ierr = NameFields(da,dof);CHKERRQ(ierr);}
88   ierr = DMDASetUniformCoordinates(da,0.0,Lx,0.0,Ly,0.0,Lz);CHKERRQ(ierr);
89   ierr = DMDAGetLocalInfo(da,&info);CHKERRQ(ierr);
90   ierr = DMCreateGlobalVector(da,&v);CHKERRQ(ierr);
91   ierr = DMDAVecGetArrayDOF(da,v,&va);CHKERRQ(ierr);
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   ierr = DMDAVecRestoreArrayDOF(da,v,&va);CHKERRQ(ierr);
102   ierr = PetscViewerVTKOpen(comm,filename,FILE_MODE_WRITE,&view);CHKERRQ(ierr);
103   ierr = VecView(v,view);CHKERRQ(ierr);
104   ierr = PetscViewerDestroy(&view);CHKERRQ(ierr);
105   ierr = VecDestroy(&v);CHKERRQ(ierr);
106   ierr = DMDestroy(&da);CHKERRQ(ierr);
107   return 0;
108 }
109 
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;
125 
126   ierr = 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);CHKERRQ(ierr);
127   ierr = DMSetFromOptions(da);CHKERRQ(ierr);
128   ierr = DMSetUp(da);CHKERRQ(ierr);
129   if (namefields) {ierr = NameFields(da,1);CHKERRQ(ierr);}
130 
131   ierr = DMDASetUniformCoordinates(da,0.0,Lx,0.0,Ly,0.0,Lz);CHKERRQ(ierr);
132   ierr = DMDAGetLocalInfo(da,&info);CHKERRQ(ierr);
133   ierr = DMDACreateCompatibleDMDA(da,dof,&daVector);CHKERRQ(ierr);
134   if (namefields) {ierr = NameFields(daVector,dof);CHKERRQ(ierr);}
135   ierr = DMCreateGlobalVector(da,&v);CHKERRQ(ierr);
136   ierr = DMCreateGlobalVector(daVector,&vVector);CHKERRQ(ierr);
137   ierr = DMDAVecGetArray(da,v,&va);CHKERRQ(ierr);
138   ierr = DMDAVecGetArrayDOF(daVector,vVector,&vVectora);CHKERRQ(ierr);
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   ierr = DMDAVecRestoreArray(da,v,&va);CHKERRQ(ierr);
153   ierr = DMDAVecRestoreArrayDOF(da,v,&vVectora);CHKERRQ(ierr);
154   ierr = PetscViewerVTKOpen(comm,filename,FILE_MODE_WRITE,&view);CHKERRQ(ierr);
155   ierr = VecView(v,view);CHKERRQ(ierr);
156   ierr = VecView(vVector,view);CHKERRQ(ierr);
157   ierr = PetscViewerDestroy(&view);CHKERRQ(ierr);
158   ierr = VecDestroy(&v);CHKERRQ(ierr);
159   ierr = VecDestroy(&vVector);CHKERRQ(ierr);
160   ierr = DMDestroy(&da);CHKERRQ(ierr);
161   ierr = DMDestroy(&daVector);CHKERRQ(ierr);
162   return 0;
163 }
164 
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;
180 
181   ierr = DMDACreate2d(comm,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DMDA_STENCIL_STAR, M,N,PETSC_DECIDE,PETSC_DECIDE,/* dof:*/ 1,sw,NULL,NULL,&da);CHKERRQ(ierr);
182   ierr = DMSetFromOptions(da);CHKERRQ(ierr);
183   ierr = DMSetUp(da);CHKERRQ(ierr);
184   if (namefields) {ierr = NameFields(da,1);CHKERRQ(ierr);}
185   ierr = DMDASetUniformCoordinates(da,0.0,Lx,0.0,Ly,0.0,Lz);CHKERRQ(ierr);
186   ierr = DMDACreateCompatibleDMDA(da,dof,&daVector);CHKERRQ(ierr);
187   if (namefields) {ierr = NameFields(daVector,dof);CHKERRQ(ierr);}
188   ierr = DMDAGetLocalInfo(da,&info);CHKERRQ(ierr);
189   ierr = DMCreateGlobalVector(da,&v);CHKERRQ(ierr);
190   ierr = DMCreateGlobalVector(daVector,&vVector);CHKERRQ(ierr);
191   ierr = DMDAVecGetArray(da,v,&va);CHKERRQ(ierr);
192   ierr = DMDAVecGetArrayDOF(daVector,vVector,&vVectora);CHKERRQ(ierr);
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   ierr = DMDAVecRestoreArray(da,v,&va);CHKERRQ(ierr);
204   ierr = DMDAVecRestoreArrayDOF(daVector,vVector,&vVectora);CHKERRQ(ierr);
205   ierr = PetscViewerVTKOpen(comm,filename,FILE_MODE_WRITE,&view);CHKERRQ(ierr);
206   ierr = VecView(v,view);CHKERRQ(ierr);
207   ierr = VecView(vVector,view);CHKERRQ(ierr);
208   ierr = PetscViewerDestroy(&view);CHKERRQ(ierr);
209   ierr = VecDestroy(&v);CHKERRQ(ierr);
210   ierr = VecDestroy(&vVector);CHKERRQ(ierr);
211   ierr = DMDestroy(&da);CHKERRQ(ierr);
212   ierr = DMDestroy(&daVector);CHKERRQ(ierr);
213   return 0;
214 }
215 
216 int main(int argc, char *argv[])
217 {
218   PetscErrorCode ierr;
219   PetscInt       dof;
220   PetscBool      namefields;
221 
222   ierr = PetscInitialize(&argc,&argv,0,help);if (ierr) return ierr;
223   dof = 2;
224   ierr = PetscOptionsGetInt(NULL,NULL,"-dof",&dof,NULL);CHKERRQ(ierr);
225   namefields = PETSC_FALSE;
226   ierr = PetscOptionsGetBool(NULL,NULL,"-namefields",&namefields,NULL);CHKERRQ(ierr);
227   ierr = test_3d("3d.vtr",dof,namefields);CHKERRQ(ierr);
228   ierr = test_2d("2d.vtr",dof,namefields);CHKERRQ(ierr);
229   ierr = test_3d_compat("3d_compat.vtr",dof,namefields);CHKERRQ(ierr);
230   ierr = test_2d_compat("2d_compat.vtr",dof,namefields);CHKERRQ(ierr);
231   ierr = test_3d("3d.vts",dof,namefields);CHKERRQ(ierr);
232   ierr = test_2d("2d.vts",dof,namefields);CHKERRQ(ierr);
233   ierr = test_3d_compat("3d_compat.vts",dof,namefields);CHKERRQ(ierr);
234   ierr = test_2d_compat("2d_compat.vts",dof,namefields);CHKERRQ(ierr);
235   ierr = PetscFinalize();
236   return ierr;
237 }
238 
239 /*TEST
240 
241    build:
242       requires: !complex
243 
244    test:
245       suffix: 1
246       nsize: 2
247       args: -dof 2
248 
249    test:
250       suffix: 2
251       nsize: 2
252       args: -dof 2
253 
254    test:
255       suffix: 3
256       nsize: 2
257       args: -dof 3
258 
259    test:
260       suffix: 4
261       nsize: 1
262       args: -dof 2 -namefields
263 
264    test:
265       suffix: 5
266       nsize: 2
267       args: -dof 4 -namefields
268 
269 TEST*/
270