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*/