Actual source code: ex42.c
1: /* -*- Mode: C++; c-basic-offset:2 ; indent-tabs-mode:nil ; -*- */
3: static char help[] = "Test VTK Rectilinear grid (.vtr) viewer support\n\n";
5: #include <petscdm.h>
6: #include <petscdmda.h>
8: /*
9: Write 3D DMDA vector with coordinates in VTK VTR format
11: */
12: PetscErrorCode test_3d(const char filename[])
13: {
14: MPI_Comm comm = MPI_COMM_WORLD;
15: const PetscInt M=10, N=15, P=30, dof=1, sw=1;
16: const PetscScalar Lx=1.0, Ly=1.0, Lz=1.0;
17: DM da;
18: Vec v;
19: PetscViewer view;
20: DMDALocalInfo info;
21: PetscScalar ***va;
22: PetscInt i,j,k;
24: 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);
25: DMSetFromOptions(da);
26: DMSetUp(da);
28: DMDASetUniformCoordinates(da,0.0,Lx,0.0,Ly,0.0,Lz);
29: DMDAGetLocalInfo(da,&info);
30: DMCreateGlobalVector(da,&v);
31: DMDAVecGetArray(da,v,&va);
32: for (k=info.zs; k<info.zs+info.zm; k++) {
33: for (j=info.ys; j<info.ys+info.ym; j++) {
34: for (i=info.xs; i<info.xs+info.xm; i++) {
35: PetscScalar x = (Lx*i)/M;
36: PetscScalar y = (Ly*j)/N;
37: PetscScalar z = (Lz*k)/P;
38: va[k][j][i] = PetscPowScalarInt(x-0.5*Lx,2)+PetscPowScalarInt(y-0.5*Ly,2)+PetscPowScalarInt(z-0.5*Lz,2);
39: }
40: }
41: }
42: DMDAVecRestoreArray(da,v,&va);
43: PetscViewerVTKOpen(comm,filename,FILE_MODE_WRITE,&view);
44: VecView(v,view);
45: PetscViewerDestroy(&view);
46: VecDestroy(&v);
47: DMDestroy(&da);
48: return 0;
49: }
51: /*
52: Write 2D DMDA vector with coordinates in VTK VTR format
54: */
55: PetscErrorCode test_2d(const char filename[])
56: {
57: MPI_Comm comm = MPI_COMM_WORLD;
58: const PetscInt M=10, N=20, dof=1, sw=1;
59: const PetscScalar Lx=1.0, Ly=1.0, Lz=1.0;
60: DM da;
61: Vec v;
62: PetscViewer view;
63: DMDALocalInfo info;
64: PetscScalar **va;
65: PetscInt i,j;
67: DMDACreate2d(comm,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DMDA_STENCIL_STAR, M,N,PETSC_DECIDE,PETSC_DECIDE,dof,sw,NULL,NULL,&da);
68: DMSetFromOptions(da);
69: DMSetUp(da);
70: DMDASetUniformCoordinates(da,0.0,Lx,0.0,Ly,0.0,Lz);
71: DMDAGetLocalInfo(da,&info);
72: DMCreateGlobalVector(da,&v);
73: DMDAVecGetArray(da,v,&va);
74: for (j=info.ys; j<info.ys+info.ym; j++) {
75: for (i=info.xs; i<info.xs+info.xm; i++) {
76: PetscScalar x = (Lx*i)/M;
77: PetscScalar y = (Ly*j)/N;
78: va[j][i] = PetscPowScalarInt(x-0.5*Lx,2)+PetscPowScalarInt(y-0.5*Ly,2);
79: }
80: }
81: DMDAVecRestoreArray(da,v,&va);
82: PetscViewerVTKOpen(comm,filename,FILE_MODE_WRITE,&view);
83: VecView(v,view);
84: PetscViewerDestroy(&view);
85: VecDestroy(&v);
86: DMDestroy(&da);
87: return 0;
88: }
90: /*
91: Write 2D DMDA vector without coordinates in VTK VTR format
93: */
94: PetscErrorCode test_2d_nocoord(const char filename[])
95: {
96: MPI_Comm comm = MPI_COMM_WORLD;
97: const PetscInt M=10, N=20, dof=1, sw=1;
98: const PetscScalar Lx=1.0, Ly=1.0;
99: DM da;
100: Vec v;
101: PetscViewer view;
102: DMDALocalInfo info;
103: PetscScalar **va;
104: PetscInt i,j;
106: DMDACreate2d(comm,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DMDA_STENCIL_STAR,M,N,PETSC_DECIDE,PETSC_DECIDE,dof,sw,NULL,NULL,&da);
107: DMSetFromOptions(da);
108: DMSetUp(da);
109: DMDAGetLocalInfo(da,&info);
110: DMCreateGlobalVector(da,&v);
111: DMDAVecGetArray(da,v,&va);
112: for (j=info.ys; j<info.ys+info.ym; j++) {
113: for (i=info.xs; i<info.xs+info.xm; i++) {
114: PetscScalar x = (Lx*i)/M;
115: PetscScalar y = (Ly*j)/N;
116: va[j][i] = PetscPowScalarInt(x-0.5*Lx,2)+PetscPowScalarInt(y-0.5*Ly,2);
117: }
118: }
119: DMDAVecRestoreArray(da,v,&va);
120: PetscViewerVTKOpen(comm,filename,FILE_MODE_WRITE,&view);
121: VecView(v,view);
122: PetscViewerDestroy(&view);
123: VecDestroy(&v);
124: DMDestroy(&da);
125: return 0;
126: }
128: /*
129: Write 3D DMDA vector without coordinates in VTK VTR format
131: */
132: PetscErrorCode test_3d_nocoord(const char filename[])
133: {
134: MPI_Comm comm = MPI_COMM_WORLD;
135: const PetscInt M=10, N=20, P=30, dof=1, sw=1;
136: const PetscScalar Lx=1.0, Ly=1.0, Lz=1.0;
137: DM da;
138: Vec v;
139: PetscViewer view;
140: DMDALocalInfo info;
141: PetscScalar ***va;
142: PetscInt i,j,k;
144: 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);
145: DMSetFromOptions(da);
146: DMSetUp(da);
148: DMDAGetLocalInfo(da,&info);
149: DMCreateGlobalVector(da,&v);
150: DMDAVecGetArray(da,v,&va);
151: for (k=info.zs; k<info.zs+info.zm; k++) {
152: for (j=info.ys; j<info.ys+info.ym; j++) {
153: for (i=info.xs; i<info.xs+info.xm; i++) {
154: PetscScalar x = (Lx*i)/M;
155: PetscScalar y = (Ly*j)/N;
156: PetscScalar z = (Lz*k)/P;
157: va[k][j][i] = PetscPowScalarInt(x-0.5*Lx,2)+PetscPowScalarInt(y-0.5*Ly,2)+PetscPowScalarInt(z-0.5*Lz,2);
158: }
159: }
160: }
161: DMDAVecRestoreArray(da,v,&va);
162: PetscViewerVTKOpen(comm,filename,FILE_MODE_WRITE,&view);
163: VecView(v,view);
164: PetscViewerDestroy(&view);
165: VecDestroy(&v);
166: DMDestroy(&da);
167: return 0;
168: }
170: int main(int argc, char *argv[])
171: {
173: PetscInitialize(&argc,&argv,0,help);
174: test_3d("3d.vtr");
175: test_2d("2d.vtr");
176: test_2d_nocoord("2d_nocoord.vtr");
177: test_3d_nocoord("3d_nocoord.vtr");
178: PetscFinalize();
179: return 0;
180: }
182: /*TEST
184: build:
185: requires: !complex
187: test:
188: nsize: 2
190: TEST*/