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