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