Actual source code: ex42.c
petsc-3.8.4 2018-03-24
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,DMDA_STENCIL_STAR, M,N,P,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,dof,sw,NULL,NULL,NULL,&da);
27: DMSetFromOptions(da);
28: DMSetUp(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,DMDA_STENCIL_STAR, M,N,PETSC_DECIDE,PETSC_DECIDE,dof,sw,NULL,NULL,&da);
72: DMSetFromOptions(da);
73: DMSetUp(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,DMDA_STENCIL_STAR,M,N,PETSC_DECIDE,PETSC_DECIDE,dof,sw,NULL,NULL,&da);
113: DMSetFromOptions(da);
114: DMSetUp(da);
115: DMDAGetLocalInfo(da,&info);
116: DMCreateGlobalVector(da,&v);
117: DMDAVecGetArray(da,v,&va);
118: for (j=info.ys; j<info.ys+info.ym; j++) {
119: for (i=info.xs; i<info.xs+info.xm; i++) {
120: PetscScalar x = (Lx*i)/M;
121: PetscScalar y = (Ly*j)/N;
122: va[j][i] = PetscPowScalarInt(x-0.5*Lx,2)+PetscPowScalarInt(y-0.5*Ly,2);
123: }
124: }
125: DMDAVecRestoreArray(da,v,&va);
126: PetscViewerVTKOpen(comm,filename,FILE_MODE_WRITE,&view);
127: VecView(v,view);
128: PetscViewerDestroy(&view);
129: VecDestroy(&v);
130: DMDestroy(&da);
131: return 0;
132: }
135: /*
136: Write 3D DMDA vector without coordinates in VTK VTR format
138: */
139: PetscErrorCode test_3d_nocoord(const char filename[])
140: {
141: MPI_Comm comm = MPI_COMM_WORLD;
142: const PetscInt M=10, N=20, P=30, dof=1, sw=1;
143: const PetscScalar Lx=1.0, Ly=1.0, Lz=1.0;
144: DM da;
145: Vec v;
146: PetscViewer view;
147: DMDALocalInfo info;
148: PetscScalar ***va;
149: PetscInt i,j,k;
150: PetscErrorCode ierr;
152: 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);
153: DMSetFromOptions(da);
154: DMSetUp(da);
156: DMDAGetLocalInfo(da,&info);
157: DMCreateGlobalVector(da,&v);
158: DMDAVecGetArray(da,v,&va);
159: for (k=info.zs; k<info.zs+info.zm; k++) {
160: for (j=info.ys; j<info.ys+info.ym; j++) {
161: for (i=info.xs; i<info.xs+info.xm; i++) {
162: PetscScalar x = (Lx*i)/M;
163: PetscScalar y = (Ly*j)/N;
164: PetscScalar z = (Lz*k)/P;
165: va[k][j][i] = PetscPowScalarInt(x-0.5*Lx,2)+PetscPowScalarInt(y-0.5*Ly,2)+PetscPowScalarInt(z-0.5*Lz,2);
166: }
167: }
168: }
169: DMDAVecRestoreArray(da,v,&va);
170: PetscViewerVTKOpen(comm,filename,FILE_MODE_WRITE,&view);
171: VecView(v,view);
172: PetscViewerDestroy(&view);
173: VecDestroy(&v);
174: DMDestroy(&da);
175: return 0;
176: }
178: int main(int argc, char *argv[])
179: {
181:
182: PetscInitialize(&argc,&argv,0,help);if (ierr) return ierr;
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 ierr;
189: }