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;
23: PetscErrorCode ierr;
25: 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);
26: DMSetFromOptions(da);
27: DMSetUp(da);
29: DMDASetUniformCoordinates(da,0.0,Lx,0.0,Ly,0.0,Lz);
30: DMDAGetLocalInfo(da,&info);
31: DMCreateGlobalVector(da,&v);
32: DMDAVecGetArray(da,v,&va);
33: for (k=info.zs; k<info.zs+info.zm; k++) {
34: for (j=info.ys; j<info.ys+info.ym; j++) {
35: for (i=info.xs; i<info.xs+info.xm; i++) {
36: PetscScalar x = (Lx*i)/M;
37: PetscScalar y = (Ly*j)/N;
38: PetscScalar z = (Lz*k)/P;
39: va[k][j][i] = PetscPowScalarInt(x-0.5*Lx,2)+PetscPowScalarInt(y-0.5*Ly,2)+PetscPowScalarInt(z-0.5*Lz,2);
40: }
41: }
42: }
43: DMDAVecRestoreArray(da,v,&va);
44: PetscViewerVTKOpen(comm,filename,FILE_MODE_WRITE,&view);
45: VecView(v,view);
46: PetscViewerDestroy(&view);
47: VecDestroy(&v);
48: DMDestroy(&da);
49: return 0;
50: }
52: /*
53: Write 2D DMDA vector with coordinates in VTK VTR format
55: */
56: PetscErrorCode test_2d(const char filename[])
57: {
58: MPI_Comm comm = MPI_COMM_WORLD;
59: const PetscInt M=10, N=20, dof=1, sw=1;
60: const PetscScalar Lx=1.0, Ly=1.0, Lz=1.0;
61: DM da;
62: Vec v;
63: PetscViewer view;
64: DMDALocalInfo info;
65: PetscScalar **va;
66: PetscInt i,j;
67: PetscErrorCode ierr;
69: DMDACreate2d(comm,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DMDA_STENCIL_STAR, M,N,PETSC_DECIDE,PETSC_DECIDE,dof,sw,NULL,NULL,&da);
70: DMSetFromOptions(da);
71: DMSetUp(da);
72: DMDASetUniformCoordinates(da,0.0,Lx,0.0,Ly,0.0,Lz);
73: DMDAGetLocalInfo(da,&info);
74: DMCreateGlobalVector(da,&v);
75: DMDAVecGetArray(da,v,&va);
76: for (j=info.ys; j<info.ys+info.ym; j++) {
77: for (i=info.xs; i<info.xs+info.xm; i++) {
78: PetscScalar x = (Lx*i)/M;
79: PetscScalar y = (Ly*j)/N;
80: va[j][i] = PetscPowScalarInt(x-0.5*Lx,2)+PetscPowScalarInt(y-0.5*Ly,2);
81: }
82: }
83: DMDAVecRestoreArray(da,v,&va);
84: PetscViewerVTKOpen(comm,filename,FILE_MODE_WRITE,&view);
85: VecView(v,view);
86: PetscViewerDestroy(&view);
87: VecDestroy(&v);
88: DMDestroy(&da);
89: return 0;
90: }
92: /*
93: Write 2D DMDA vector without coordinates in VTK VTR format
95: */
96: PetscErrorCode test_2d_nocoord(const char filename[])
97: {
98: MPI_Comm comm = MPI_COMM_WORLD;
99: const PetscInt M=10, N=20, dof=1, sw=1;
100: const PetscScalar Lx=1.0, Ly=1.0;
101: DM da;
102: Vec v;
103: PetscViewer view;
104: DMDALocalInfo info;
105: PetscScalar **va;
106: PetscInt i,j;
107: PetscErrorCode ierr;
109: DMDACreate2d(comm,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DMDA_STENCIL_STAR,M,N,PETSC_DECIDE,PETSC_DECIDE,dof,sw,NULL,NULL,&da);
110: DMSetFromOptions(da);
111: DMSetUp(da);
112: DMDAGetLocalInfo(da,&info);
113: DMCreateGlobalVector(da,&v);
114: DMDAVecGetArray(da,v,&va);
115: for (j=info.ys; j<info.ys+info.ym; j++) {
116: for (i=info.xs; i<info.xs+info.xm; i++) {
117: PetscScalar x = (Lx*i)/M;
118: PetscScalar y = (Ly*j)/N;
119: va[j][i] = PetscPowScalarInt(x-0.5*Lx,2)+PetscPowScalarInt(y-0.5*Ly,2);
120: }
121: }
122: DMDAVecRestoreArray(da,v,&va);
123: PetscViewerVTKOpen(comm,filename,FILE_MODE_WRITE,&view);
124: VecView(v,view);
125: PetscViewerDestroy(&view);
126: VecDestroy(&v);
127: DMDestroy(&da);
128: return 0;
129: }
131: /*
132: Write 3D DMDA vector without coordinates in VTK VTR format
134: */
135: PetscErrorCode test_3d_nocoord(const char filename[])
136: {
137: MPI_Comm comm = MPI_COMM_WORLD;
138: const PetscInt M=10, N=20, P=30, dof=1, sw=1;
139: const PetscScalar Lx=1.0, Ly=1.0, Lz=1.0;
140: DM da;
141: Vec v;
142: PetscViewer view;
143: DMDALocalInfo info;
144: PetscScalar ***va;
145: PetscInt i,j,k;
146: PetscErrorCode ierr;
148: 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);
149: DMSetFromOptions(da);
150: DMSetUp(da);
152: DMDAGetLocalInfo(da,&info);
153: DMCreateGlobalVector(da,&v);
154: DMDAVecGetArray(da,v,&va);
155: for (k=info.zs; k<info.zs+info.zm; k++) {
156: for (j=info.ys; j<info.ys+info.ym; j++) {
157: for (i=info.xs; i<info.xs+info.xm; i++) {
158: PetscScalar x = (Lx*i)/M;
159: PetscScalar y = (Ly*j)/N;
160: PetscScalar z = (Lz*k)/P;
161: va[k][j][i] = PetscPowScalarInt(x-0.5*Lx,2)+PetscPowScalarInt(y-0.5*Ly,2)+PetscPowScalarInt(z-0.5*Lz,2);
162: }
163: }
164: }
165: DMDAVecRestoreArray(da,v,&va);
166: PetscViewerVTKOpen(comm,filename,FILE_MODE_WRITE,&view);
167: VecView(v,view);
168: PetscViewerDestroy(&view);
169: VecDestroy(&v);
170: DMDestroy(&da);
171: return 0;
172: }
174: int main(int argc, char *argv[])
175: {
178: PetscInitialize(&argc,&argv,0,help);if (ierr) return ierr;
179: test_3d("3d.vtr");
180: test_2d("2d.vtr");
181: test_2d_nocoord("2d_nocoord.vtr");
182: test_3d_nocoord("3d_nocoord.vtr");
183: PetscFinalize();
184: return ierr;
185: }
187: /*TEST
189: build:
190: requires: !complex
192: test:
193: nsize: 2
195: TEST*/