Actual source code: ex42.c
petsc-3.5.4 2015-05-23
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 <mpi.h>
6: #include <petscdm.h>
7: #include petscdmda.h
9: /*
10: Write 3D DMDA vector with coordinates in VTK VTR format
12: */
13: PetscErrorCode test_3d(const char filename[])
14: {
15: MPI_Comm comm = MPI_COMM_WORLD;
16: const PetscInt M=10, N=15, P=30, dof=1, sw=1;
17: const PetscScalar Lx=1.0, Ly=1.0, Lz=1.0;
18: DM da;
19: Vec v;
20: PetscViewer view;
21: DMDALocalInfo info;
22: PetscScalar ***va;
23: PetscInt i,j,k;
24: PetscErrorCode ierr;
26: DMDACreate3d(comm,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,
27: DMDA_STENCIL_STAR, M,N,P,
28: PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,dof,sw,NULL,NULL,NULL,&da);
30: DMDASetUniformCoordinates(da,0.0,Lx,0.0,Ly,0.0,Lz);
31: DMDAGetLocalInfo(da,&info);
32: DMCreateGlobalVector(da,&v);
33: DMDAVecGetArray(da,v,&va);
34: for (k=info.zs; k<info.zs+info.zm; k++) {
35: for (j=info.ys; j<info.ys+info.ym; j++) {
36: for (i=info.xs; i<info.xs+info.xm; i++) {
37: PetscScalar x = (Lx*i)/M;
38: PetscScalar y = (Ly*j)/N;
39: PetscScalar z = (Lz*k)/P;
40: va[k][j][i] = PetscPowScalarInt(x-0.5*Lx,2)+PetscPowScalarInt(y-0.5*Ly,2)+PetscPowScalarInt(z-0.5*Lz,2);
41: }
42: }
43: }
44: DMDAVecRestoreArray(da,v,&va);
45: PetscViewerVTKOpen(comm,filename,FILE_MODE_WRITE,&view);
46: VecView(v,view);
47: PetscViewerDestroy(&view);
48: VecDestroy(&v);
49: DMDestroy(&da);
50: return 0;
51: }
54: /*
55: Write 2D DMDA vector with coordinates in VTK VTR format
57: */
58: PetscErrorCode test_2d(const char filename[])
59: {
60: MPI_Comm comm = MPI_COMM_WORLD;
61: const PetscInt M=10, N=20, dof=1, sw=1;
62: const PetscScalar Lx=1.0, Ly=1.0, Lz=1.0;
63: DM da;
64: Vec v;
65: PetscViewer view;
66: DMDALocalInfo info;
67: PetscScalar **va;
68: PetscInt i,j;
69: PetscErrorCode ierr;
71: DMDACreate2d(comm,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,
72: DMDA_STENCIL_STAR, M,N,
73: PETSC_DECIDE,PETSC_DECIDE,dof,sw,NULL,NULL,&da);
74: DMDASetUniformCoordinates(da,0.0,Lx,0.0,Ly,0.0,Lz);
75: DMDAGetLocalInfo(da,&info);
76: DMCreateGlobalVector(da,&v);
77: DMDAVecGetArray(da,v,&va);
78: for (j=info.ys; j<info.ys+info.ym; j++) {
79: for (i=info.xs; i<info.xs+info.xm; i++) {
80: PetscScalar x = (Lx*i)/M;
81: PetscScalar y = (Ly*j)/N;
82: va[j][i] = PetscPowScalarInt(x-0.5*Lx,2)+PetscPowScalarInt(y-0.5*Ly,2);
83: }
84: }
85: DMDAVecRestoreArray(da,v,&va);
86: PetscViewerVTKOpen(comm,filename,FILE_MODE_WRITE,&view);
87: VecView(v,view);
88: PetscViewerDestroy(&view);
89: VecDestroy(&v);
90: DMDestroy(&da);
91: return 0;
92: }
95: /*
96: Write 2D DMDA vector without coordinates in VTK VTR format
98: */
99: PetscErrorCode test_2d_nocoord(const char filename[])
100: {
101: MPI_Comm comm = MPI_COMM_WORLD;
102: const PetscInt M=10, N=20, dof=1, sw=1;
103: const PetscScalar Lx=1.0, Ly=1.0;
104: DM da;
105: Vec v;
106: PetscViewer view;
107: DMDALocalInfo info;
108: PetscScalar **va;
109: PetscInt i,j;
110: PetscErrorCode ierr;
112: DMDACreate2d(comm,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,
113: DMDA_STENCIL_STAR, M,N,
114: PETSC_DECIDE,PETSC_DECIDE,dof,sw,NULL,NULL,&da);
116: DMDAGetLocalInfo(da,&info);
117: DMCreateGlobalVector(da,&v);
118: DMDAVecGetArray(da,v,&va);
119: for (j=info.ys; j<info.ys+info.ym; j++) {
120: for (i=info.xs; i<info.xs+info.xm; i++) {
121: PetscScalar x = (Lx*i)/M;
122: PetscScalar y = (Ly*j)/N;
123: va[j][i] = PetscPowScalarInt(x-0.5*Lx,2)+PetscPowScalarInt(y-0.5*Ly,2);
124: }
125: }
126: DMDAVecRestoreArray(da,v,&va);
127: PetscViewerVTKOpen(comm,filename,FILE_MODE_WRITE,&view);
128: VecView(v,view);
129: PetscViewerDestroy(&view);
130: VecDestroy(&v);
131: DMDestroy(&da);
132: return 0;
133: }
136: /*
137: Write 3D DMDA vector without coordinates in VTK VTR format
139: */
140: PetscErrorCode test_3d_nocoord(const char filename[])
141: {
142: MPI_Comm comm = MPI_COMM_WORLD;
143: const PetscInt M=10, N=20, P=30, dof=1, sw=1;
144: const PetscScalar Lx=1.0, Ly=1.0, Lz=1.0;
145: DM da;
146: Vec v;
147: PetscViewer view;
148: DMDALocalInfo info;
149: PetscScalar ***va;
150: PetscInt i,j,k;
151: PetscErrorCode ierr;
153: DMDACreate3d(comm,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,
154: DMDA_STENCIL_STAR, M,N,P,
155: PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,dof,sw,NULL,NULL,NULL,&da);
157: DMDAGetLocalInfo(da,&info);
158: DMCreateGlobalVector(da,&v);
159: DMDAVecGetArray(da,v,&va);
160: for (k=info.zs; k<info.zs+info.zm; k++) {
161: for (j=info.ys; j<info.ys+info.ym; j++) {
162: for (i=info.xs; i<info.xs+info.xm; i++) {
163: PetscScalar x = (Lx*i)/M;
164: PetscScalar y = (Ly*j)/N;
165: PetscScalar z = (Lz*k)/P;
166: va[k][j][i] = PetscPowScalarInt(x-0.5*Lx,2)+PetscPowScalarInt(y-0.5*Ly,2)+PetscPowScalarInt(z-0.5*Lz,2);
167: }
168: }
169: }
170: DMDAVecRestoreArray(da,v,&va);
171: PetscViewerVTKOpen(comm,filename,FILE_MODE_WRITE,&view);
172: VecView(v,view);
173: PetscViewerDestroy(&view);
174: VecDestroy(&v);
175: DMDestroy(&da);
176: return 0;
177: }
179: int main(int argc, char *argv[])
180: {
182: PetscInitialize(&argc,&argv,0,help);
183: test_3d("3d.vtr");
184: test_2d("2d.vtr");
185: test_2d_nocoord("2d_nocoord.vtr");
186: test_3d_nocoord("3d_nocoord.vtr");
187: PetscFinalize();
188: return 0;
189: }