xref: /petsc/src/dm/tests/ex48.c (revision c4762a1b19cd2af06abeed90e8f9d34fb975dd94) 
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 /*
70   Write 2D DMDA vector with coordinates in VTK format
71 */
72 PetscErrorCode test_2d(const char filename[],PetscInt dof,PetscBool namefields)
73 {
74   MPI_Comm          comm = MPI_COMM_WORLD;
75   const PetscInt    M=10,N=20,sw=1;
76   const PetscScalar Lx=1.0,Ly=1.0,Lz=1.0;
77   DM                da;
78   Vec               v;
79   PetscViewer       view;
80   DMDALocalInfo     info;
81   PetscScalar       ***va;
82   PetscInt          i,j,c;
83   PetscErrorCode    ierr;
84 
85   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);
86   ierr = DMSetFromOptions(da);CHKERRQ(ierr);
87   ierr = DMSetUp(da);CHKERRQ(ierr);
88   if (namefields) {ierr = NameFields(da,dof);CHKERRQ(ierr);}
89   ierr = DMDASetUniformCoordinates(da,0.0,Lx,0.0,Ly,0.0,Lz);CHKERRQ(ierr);
90   ierr = DMDAGetLocalInfo(da,&info);CHKERRQ(ierr);
91   ierr = DMCreateGlobalVector(da,&v);CHKERRQ(ierr);
92   ierr = DMDAVecGetArrayDOF(da,v,&va);CHKERRQ(ierr);
93   for (j=info.ys; j<info.ys+info.ym; j++) {
94     for (i=info.xs; i<info.xs+info.xm; i++) {
95       const PetscScalar x = (Lx*i)/M;
96       const PetscScalar y = (Ly*j)/N;
97       for (c=0; c<dof; ++c) {
98         va[j][i][c] = PetscPowScalarInt(x-0.5*Lx,2)+PetscPowScalarInt(y-0.5*Ly,2) + 10.0*c;
99       }
100     }
101   }
102   ierr = DMDAVecRestoreArrayDOF(da,v,&va);CHKERRQ(ierr);
103   ierr = PetscViewerVTKOpen(comm,filename,FILE_MODE_WRITE,&view);CHKERRQ(ierr);
104   ierr = VecView(v,view);CHKERRQ(ierr);
105   ierr = PetscViewerDestroy(&view);CHKERRQ(ierr);
106   ierr = VecDestroy(&v);CHKERRQ(ierr);
107   ierr = DMDestroy(&da);CHKERRQ(ierr);
108   return 0;
109 }
110 
111 /*
112   Write a scalar and a vector field from two compatible 3d DMDAs
113 */
114 PetscErrorCode test_3d_compat(const char filename[],PetscInt dof,PetscBool namefields)
115 {
116   MPI_Comm          comm = MPI_COMM_WORLD;
117   const PetscInt    M=10,N=15,P=30,sw=1;
118   const PetscScalar Lx=1.0,Ly=1.0,Lz=1.0;
119   DM                da,daVector;
120   Vec               v,vVector;
121   PetscViewer       view;
122   DMDALocalInfo     info;
123   PetscScalar       ***va,****vVectora;
124   PetscInt          i,j,k,c;
125   PetscErrorCode    ierr;
126 
127   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);
128   ierr = DMSetFromOptions(da);CHKERRQ(ierr);
129   ierr = DMSetUp(da);CHKERRQ(ierr);
130   if (namefields) {ierr = NameFields(da,1);CHKERRQ(ierr);}
131 
132   ierr = DMDASetUniformCoordinates(da,0.0,Lx,0.0,Ly,0.0,Lz);CHKERRQ(ierr);
133   ierr = DMDAGetLocalInfo(da,&info);CHKERRQ(ierr);
134   ierr = DMDACreateCompatibleDMDA(da,dof,&daVector);CHKERRQ(ierr);
135   if (namefields) {ierr = NameFields(daVector,dof);CHKERRQ(ierr);}
136   ierr = DMCreateGlobalVector(da,&v);CHKERRQ(ierr);
137   ierr = DMCreateGlobalVector(daVector,&vVector);CHKERRQ(ierr);
138   ierr = DMDAVecGetArray(da,v,&va);CHKERRQ(ierr);
139   ierr = DMDAVecGetArrayDOF(daVector,vVector,&vVectora);CHKERRQ(ierr);
140   for (k=info.zs; k<info.zs+info.zm; k++) {
141     for (j=info.ys; j<info.ys+info.ym; j++) {
142       for (i=info.xs; i<info.xs+info.xm; i++) {
143         const PetscScalar x = (Lx*i)/M;
144         const PetscScalar y = (Ly*j)/N;
145         const PetscScalar z = (Lz*k)/P;
146         va[k][j][i] = PetscPowScalarInt(x-0.5*Lx,2)+PetscPowScalarInt(y-0.5*Ly,2)+PetscPowScalarInt(z-0.5*Lz,2);
147         for (c=0; c<dof; ++c) {
148           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;
149         }
150       }
151     }
152   }
153   ierr = DMDAVecRestoreArray(da,v,&va);CHKERRQ(ierr);
154   ierr = DMDAVecRestoreArrayDOF(da,v,&vVectora);CHKERRQ(ierr);
155   ierr = PetscViewerVTKOpen(comm,filename,FILE_MODE_WRITE,&view);CHKERRQ(ierr);
156   ierr = VecView(v,view);CHKERRQ(ierr);
157   ierr = VecView(vVector,view);CHKERRQ(ierr);
158   ierr = PetscViewerDestroy(&view);CHKERRQ(ierr);
159   ierr = VecDestroy(&v);CHKERRQ(ierr);
160   ierr = VecDestroy(&vVector);CHKERRQ(ierr);
161   ierr = DMDestroy(&da);CHKERRQ(ierr);
162   ierr = DMDestroy(&daVector);CHKERRQ(ierr);
163   return 0;
164 }
165 
166 /*
167   Write a scalar and a vector field from two compatible 2d DMDAs
168 */
169 PetscErrorCode test_2d_compat(const char filename[],PetscInt dof,PetscBool namefields)
170 {
171   MPI_Comm          comm = MPI_COMM_WORLD;
172   const PetscInt    M=10,N=20,sw=1;
173   const PetscScalar Lx=1.0,Ly=1.0,Lz=1.0;
174   DM                da, daVector;
175   Vec               v,vVector;
176   PetscViewer       view;
177   DMDALocalInfo     info;
178   PetscScalar       **va,***vVectora;
179   PetscInt          i,j,c;
180   PetscErrorCode    ierr;
181 
182   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);
183   ierr = DMSetFromOptions(da);CHKERRQ(ierr);
184   ierr = DMSetUp(da);CHKERRQ(ierr);
185   if (namefields) {ierr = NameFields(da,1);CHKERRQ(ierr);}
186   ierr = DMDASetUniformCoordinates(da,0.0,Lx,0.0,Ly,0.0,Lz);CHKERRQ(ierr);
187   ierr = DMDACreateCompatibleDMDA(da,dof,&daVector);CHKERRQ(ierr);
188   if (namefields) {ierr = NameFields(daVector,dof);CHKERRQ(ierr);}
189   ierr = DMDAGetLocalInfo(da,&info);CHKERRQ(ierr);
190   ierr = DMCreateGlobalVector(da,&v);CHKERRQ(ierr);
191   ierr = DMCreateGlobalVector(daVector,&vVector);CHKERRQ(ierr);
192   ierr = DMDAVecGetArray(da,v,&va);CHKERRQ(ierr);
193   ierr = DMDAVecGetArrayDOF(daVector,vVector,&vVectora);CHKERRQ(ierr);
194   for (j=info.ys; j<info.ys+info.ym; j++) {
195     for (i=info.xs; i<info.xs+info.xm; i++) {
196       const PetscScalar x = (Lx*i)/M;
197       const PetscScalar y = (Ly*j)/N;
198       va[j][i] = PetscPowScalarInt(x-0.5*Lx,2)+PetscPowScalarInt(y-0.5*Ly,2);
199       for (c=0; c<dof; ++c) {
200         vVectora[j][i][c] = PetscPowScalarInt(x-0.5*Lx,2)+PetscPowScalarInt(y-0.5*Ly,2) + 10.0*c;
201       }
202     }
203   }
204   ierr = DMDAVecRestoreArray(da,v,&va);CHKERRQ(ierr);
205   ierr = DMDAVecRestoreArrayDOF(daVector,vVector,&vVectora);CHKERRQ(ierr);
206   ierr = PetscViewerVTKOpen(comm,filename,FILE_MODE_WRITE,&view);CHKERRQ(ierr);
207   ierr = VecView(v,view);CHKERRQ(ierr);
208   ierr = VecView(vVector,view);CHKERRQ(ierr);
209   ierr = PetscViewerDestroy(&view);CHKERRQ(ierr);
210   ierr = VecDestroy(&v);CHKERRQ(ierr);
211   ierr = VecDestroy(&vVector);CHKERRQ(ierr);
212   ierr = DMDestroy(&da);CHKERRQ(ierr);
213   ierr = DMDestroy(&daVector);CHKERRQ(ierr);
214   return 0;
215 }
216 
217 int main(int argc, char *argv[])
218 {
219   PetscErrorCode ierr;
220   PetscInt       dof;
221   PetscBool      namefields;
222 
223   ierr = PetscInitialize(&argc,&argv,0,help);if (ierr) return ierr;
224   dof = 2;
225   ierr = PetscOptionsGetInt(NULL,NULL,"-dof",&dof,NULL);CHKERRQ(ierr);
226   namefields = PETSC_FALSE;
227   ierr = PetscOptionsGetBool(NULL,NULL,"-namefields",&namefields,NULL);CHKERRQ(ierr);
228   ierr = test_3d("3d.vtr",dof,namefields);CHKERRQ(ierr);
229   ierr = test_2d("2d.vtr",dof,namefields);CHKERRQ(ierr);
230   ierr = test_3d_compat("3d_compat.vtr",dof,namefields);CHKERRQ(ierr);
231   ierr = test_2d_compat("2d_compat.vtr",dof,namefields);CHKERRQ(ierr);
232   ierr = test_3d("3d.vts",dof,namefields);CHKERRQ(ierr);
233   ierr = test_2d("2d.vts",dof,namefields);CHKERRQ(ierr);
234   ierr = test_3d_compat("3d_compat.vts",dof,namefields);CHKERRQ(ierr);
235   ierr = test_2d_compat("2d_compat.vts",dof,namefields);CHKERRQ(ierr);
236   ierr = PetscFinalize();
237   return ierr;
238 }
239 
240 /*TEST
241 
242    build:
243       requires: !complex
244 
245    test:
246       suffix: 1
247       nsize: 2
248       args: -dof 2
249 
250    test:
251       suffix: 2
252       nsize: 2
253       args: -dof 2
254 
255    test:
256       suffix: 3
257       nsize: 2
258       args: -dof 3
259 
260    test:
261       suffix: 4
262       nsize: 1
263       args: -dof 2 -namefields
264 
265    test:
266       suffix: 5
267       nsize: 2
268       args: -dof 4 -namefields
269 
270 TEST*/
271