1 /* -*- Mode: C++; c-basic-offset:2 ; indent-tabs-mode:nil ; -*- */
2
3 static char help[] = "Test VTK Rectilinear grid (.vtr) viewer support\n\n";
4
5 #include <petscdm.h>
6 #include <petscdmda.h>
7
8 /*
9 Write 3D DMDA vector with coordinates in VTK VTR format
10
11 */
test_3d(const char filename[])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
24 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));
25 PetscCall(DMSetFromOptions(da));
26 PetscCall(DMSetUp(da));
27
28 PetscCall(DMDASetUniformCoordinates(da, 0.0, Lx, 0.0, Ly, 0.0, Lz));
29 PetscCall(DMDAGetLocalInfo(da, &info));
30 PetscCall(DMCreateGlobalVector(da, &v));
31 PetscCall(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 PetscCall(DMDAVecRestoreArray(da, v, &va));
43 PetscCall(PetscViewerVTKOpen(comm, filename, FILE_MODE_WRITE, &view));
44 PetscCall(VecView(v, view));
45 PetscCall(PetscViewerDestroy(&view));
46 PetscCall(VecDestroy(&v));
47 PetscCall(DMDestroy(&da));
48 return PETSC_SUCCESS;
49 }
50
51 /*
52 Write 2D DMDA vector with coordinates in VTK VTR format
53
54 */
test_2d(const char filename[])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;
66
67 PetscCall(DMDACreate2d(comm, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_STAR, M, N, PETSC_DECIDE, PETSC_DECIDE, dof, sw, NULL, NULL, &da));
68 PetscCall(DMSetFromOptions(da));
69 PetscCall(DMSetUp(da));
70 PetscCall(DMDASetUniformCoordinates(da, 0.0, Lx, 0.0, Ly, 0.0, Lz));
71 PetscCall(DMDAGetLocalInfo(da, &info));
72 PetscCall(DMCreateGlobalVector(da, &v));
73 PetscCall(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 PetscCall(DMDAVecRestoreArray(da, v, &va));
82 PetscCall(PetscViewerVTKOpen(comm, filename, FILE_MODE_WRITE, &view));
83 PetscCall(VecView(v, view));
84 PetscCall(PetscViewerDestroy(&view));
85 PetscCall(VecDestroy(&v));
86 PetscCall(DMDestroy(&da));
87 return PETSC_SUCCESS;
88 }
89
90 /*
91 Write 2D DMDA vector without coordinates in VTK VTR format
92
93 */
test_2d_nocoord(const char filename[])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;
105
106 PetscCall(DMDACreate2d(comm, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_STAR, M, N, PETSC_DECIDE, PETSC_DECIDE, dof, sw, NULL, NULL, &da));
107 PetscCall(DMSetFromOptions(da));
108 PetscCall(DMSetUp(da));
109 PetscCall(DMDAGetLocalInfo(da, &info));
110 PetscCall(DMCreateGlobalVector(da, &v));
111 PetscCall(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 PetscCall(DMDAVecRestoreArray(da, v, &va));
120 PetscCall(PetscViewerVTKOpen(comm, filename, FILE_MODE_WRITE, &view));
121 PetscCall(VecView(v, view));
122 PetscCall(PetscViewerDestroy(&view));
123 PetscCall(VecDestroy(&v));
124 PetscCall(DMDestroy(&da));
125 return PETSC_SUCCESS;
126 }
127
128 /*
129 Write 3D DMDA vector without coordinates in VTK VTR format
130
131 */
test_3d_nocoord(const char filename[])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;
143
144 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));
145 PetscCall(DMSetFromOptions(da));
146 PetscCall(DMSetUp(da));
147
148 PetscCall(DMDAGetLocalInfo(da, &info));
149 PetscCall(DMCreateGlobalVector(da, &v));
150 PetscCall(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 PetscCall(DMDAVecRestoreArray(da, v, &va));
162 PetscCall(PetscViewerVTKOpen(comm, filename, FILE_MODE_WRITE, &view));
163 PetscCall(VecView(v, view));
164 PetscCall(PetscViewerDestroy(&view));
165 PetscCall(VecDestroy(&v));
166 PetscCall(DMDestroy(&da));
167 return PETSC_SUCCESS;
168 }
169
main(int argc,char * argv[])170 int main(int argc, char *argv[])
171 {
172 PetscFunctionBeginUser;
173 PetscCall(PetscInitialize(&argc, &argv, 0, help));
174 PetscCall(test_3d("3d.vtr"));
175 PetscCall(test_2d("2d.vtr"));
176 PetscCall(test_2d_nocoord("2d_nocoord.vtr"));
177 PetscCall(test_3d_nocoord("3d_nocoord.vtr"));
178 PetscCall(PetscFinalize());
179 return 0;
180 }
181
182 /*TEST
183
184 build:
185 requires: !complex
186
187 test:
188 nsize: 2
189 output_file: output/empty.out
190
191 TEST*/
192