1 static char help[] = "Test VTK structured (.vts) and rectilinear (.vtr) viewer support with multi-dof DMDAs.\n\
2 Supply the -namefields flag to test with field names.\n\n";
3
4 #include <petscdm.h>
5 #include <petscdmda.h>
6
7 /* Helper function to name DMDA fields */
NameFields(DM da,PetscInt dof)8 PetscErrorCode NameFields(DM da, PetscInt dof)
9 {
10 PetscInt c;
11
12 PetscFunctionBeginUser;
13 for (c = 0; c < dof; ++c) {
14 char fieldname[256];
15 PetscCall(PetscSNPrintf(fieldname, sizeof(fieldname), "field_%" PetscInt_FMT, c));
16 PetscCall(DMDASetFieldName(da, c, fieldname));
17 }
18 PetscFunctionReturn(PETSC_SUCCESS);
19 }
20
21 /*
22 Write 3D DMDA vector with coordinates in VTK format
23 */
test_3d(const char filename[],PetscInt dof,PetscBool namefields)24 PetscErrorCode test_3d(const char filename[], PetscInt dof, PetscBool namefields)
25 {
26 MPI_Comm comm = MPI_COMM_WORLD;
27 const PetscInt M = 10, N = 15, P = 30, sw = 1;
28 const PetscScalar Lx = 1.0, Ly = 1.0, Lz = 1.0;
29 DM da;
30 Vec v;
31 PetscViewer view;
32 DMDALocalInfo info;
33 PetscScalar ****va;
34 PetscInt i, j, k, c;
35
36 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));
37 PetscCall(DMSetFromOptions(da));
38 PetscCall(DMSetUp(da));
39 if (namefields) PetscCall(NameFields(da, dof));
40
41 PetscCall(DMDASetUniformCoordinates(da, 0.0, Lx, 0.0, Ly, 0.0, Lz));
42 PetscCall(DMDAGetLocalInfo(da, &info));
43 PetscCall(DMCreateGlobalVector(da, &v));
44 PetscCall(DMDAVecGetArrayDOF(da, v, &va));
45 for (k = info.zs; k < info.zs + info.zm; k++) {
46 for (j = info.ys; j < info.ys + info.ym; j++) {
47 for (i = info.xs; i < info.xs + info.xm; i++) {
48 const PetscScalar x = (Lx * i) / M;
49 const PetscScalar y = (Ly * j) / N;
50 const PetscScalar z = (Lz * k) / P;
51 for (c = 0; c < dof; ++c) va[k][j][i][c] = PetscPowScalarInt(x - 0.5 * Lx, 2) + PetscPowScalarInt(y - 0.5 * Ly, 2) + PetscPowScalarInt(z - 0.5 * Lz, 2) + 10.0 * c;
52 }
53 }
54 }
55 PetscCall(DMDAVecRestoreArrayDOF(da, v, &va));
56 PetscCall(PetscViewerVTKOpen(comm, filename, FILE_MODE_WRITE, &view));
57 PetscCall(VecView(v, view));
58 PetscCall(PetscViewerDestroy(&view));
59 PetscCall(VecDestroy(&v));
60 PetscCall(DMDestroy(&da));
61 return PETSC_SUCCESS;
62 }
63
64 /*
65 Write 2D DMDA vector with coordinates in VTK format
66 */
test_2d(const char filename[],PetscInt dof,PetscBool namefields)67 PetscErrorCode test_2d(const char filename[], PetscInt dof, PetscBool namefields)
68 {
69 MPI_Comm comm = MPI_COMM_WORLD;
70 const PetscInt M = 10, N = 20, sw = 1;
71 const PetscScalar Lx = 1.0, Ly = 1.0, Lz = 1.0;
72 DM da;
73 Vec v;
74 PetscViewer view;
75 DMDALocalInfo info;
76 PetscScalar ***va;
77 PetscInt i, j, c;
78
79 PetscCall(DMDACreate2d(comm, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_STAR, M, N, PETSC_DECIDE, PETSC_DECIDE, dof, sw, NULL, NULL, &da));
80 PetscCall(DMSetFromOptions(da));
81 PetscCall(DMSetUp(da));
82 if (namefields) PetscCall(NameFields(da, dof));
83 PetscCall(DMDASetUniformCoordinates(da, 0.0, Lx, 0.0, Ly, 0.0, Lz));
84 PetscCall(DMDAGetLocalInfo(da, &info));
85 PetscCall(DMCreateGlobalVector(da, &v));
86 PetscCall(DMDAVecGetArrayDOF(da, v, &va));
87 for (j = info.ys; j < info.ys + info.ym; j++) {
88 for (i = info.xs; i < info.xs + info.xm; i++) {
89 const PetscScalar x = (Lx * i) / M;
90 const PetscScalar y = (Ly * j) / N;
91 for (c = 0; c < dof; ++c) va[j][i][c] = PetscPowScalarInt(x - 0.5 * Lx, 2) + PetscPowScalarInt(y - 0.5 * Ly, 2) + 10.0 * c;
92 }
93 }
94 PetscCall(DMDAVecRestoreArrayDOF(da, v, &va));
95 PetscCall(PetscViewerVTKOpen(comm, filename, FILE_MODE_WRITE, &view));
96 PetscCall(VecView(v, view));
97 PetscCall(PetscViewerDestroy(&view));
98 PetscCall(VecDestroy(&v));
99 PetscCall(DMDestroy(&da));
100 return PETSC_SUCCESS;
101 }
102
103 /*
104 Write a scalar and a vector field from two compatible 3d DMDAs
105 */
test_3d_compat(const char filename[],PetscInt dof,PetscBool namefields)106 PetscErrorCode test_3d_compat(const char filename[], PetscInt dof, PetscBool namefields)
107 {
108 MPI_Comm comm = MPI_COMM_WORLD;
109 const PetscInt M = 10, N = 15, P = 30, sw = 1;
110 const PetscScalar Lx = 1.0, Ly = 1.0, Lz = 1.0;
111 DM da, daVector;
112 Vec v, vVector;
113 PetscViewer view;
114 DMDALocalInfo info;
115 PetscScalar ***va, ****vVectora;
116 PetscInt i, j, k, c;
117
118 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:*/ 1, sw, NULL, NULL, NULL, &da));
119 PetscCall(DMSetFromOptions(da));
120 PetscCall(DMSetUp(da));
121 if (namefields) PetscCall(NameFields(da, 1));
122
123 PetscCall(DMDASetUniformCoordinates(da, 0.0, Lx, 0.0, Ly, 0.0, Lz));
124 PetscCall(DMDAGetLocalInfo(da, &info));
125 PetscCall(DMDACreateCompatibleDMDA(da, dof, &daVector));
126 if (namefields) PetscCall(NameFields(daVector, dof));
127 PetscCall(DMCreateGlobalVector(da, &v));
128 PetscCall(DMCreateGlobalVector(daVector, &vVector));
129 PetscCall(DMDAVecGetArray(da, v, &va));
130 PetscCall(DMDAVecGetArrayDOF(daVector, vVector, &vVectora));
131 for (k = info.zs; k < info.zs + info.zm; k++) {
132 for (j = info.ys; j < info.ys + info.ym; j++) {
133 for (i = info.xs; i < info.xs + info.xm; i++) {
134 const PetscScalar x = (Lx * i) / M;
135 const PetscScalar y = (Ly * j) / N;
136 const PetscScalar z = (Lz * k) / P;
137 va[k][j][i] = PetscPowScalarInt(x - 0.5 * Lx, 2) + PetscPowScalarInt(y - 0.5 * Ly, 2) + PetscPowScalarInt(z - 0.5 * Lz, 2);
138 for (c = 0; c < dof; ++c) vVectora[k][j][i][c] = PetscPowScalarInt(x - 0.5 * Lx, 2) + PetscPowScalarInt(y - 0.5 * Ly, 2) + PetscPowScalarInt(z - 0.5 * Lz, 2) + 10.0 * c;
139 }
140 }
141 }
142 PetscCall(DMDAVecRestoreArray(da, v, &va));
143 PetscCall(DMDAVecRestoreArrayDOF(da, v, &vVectora));
144 PetscCall(PetscViewerVTKOpen(comm, filename, FILE_MODE_WRITE, &view));
145 PetscCall(VecView(v, view));
146 PetscCall(VecView(vVector, view));
147 PetscCall(PetscViewerDestroy(&view));
148 PetscCall(VecDestroy(&v));
149 PetscCall(VecDestroy(&vVector));
150 PetscCall(DMDestroy(&da));
151 PetscCall(DMDestroy(&daVector));
152 return PETSC_SUCCESS;
153 }
154
155 /*
156 Write a scalar and a vector field from two compatible 2d DMDAs
157 */
test_2d_compat(const char filename[],PetscInt dof,PetscBool namefields)158 PetscErrorCode test_2d_compat(const char filename[], PetscInt dof, PetscBool namefields)
159 {
160 MPI_Comm comm = MPI_COMM_WORLD;
161 const PetscInt M = 10, N = 20, sw = 1;
162 const PetscScalar Lx = 1.0, Ly = 1.0, Lz = 1.0;
163 DM da, daVector;
164 Vec v, vVector;
165 PetscViewer view;
166 DMDALocalInfo info;
167 PetscScalar **va, ***vVectora;
168 PetscInt i, j, c;
169
170 PetscCall(DMDACreate2d(comm, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_STAR, M, N, PETSC_DECIDE, PETSC_DECIDE, /* dof:*/ 1, sw, NULL, NULL, &da));
171 PetscCall(DMSetFromOptions(da));
172 PetscCall(DMSetUp(da));
173 if (namefields) PetscCall(NameFields(da, 1));
174 PetscCall(DMDASetUniformCoordinates(da, 0.0, Lx, 0.0, Ly, 0.0, Lz));
175 PetscCall(DMDACreateCompatibleDMDA(da, dof, &daVector));
176 if (namefields) PetscCall(NameFields(daVector, dof));
177 PetscCall(DMDAGetLocalInfo(da, &info));
178 PetscCall(DMCreateGlobalVector(da, &v));
179 PetscCall(DMCreateGlobalVector(daVector, &vVector));
180 PetscCall(DMDAVecGetArray(da, v, &va));
181 PetscCall(DMDAVecGetArrayDOF(daVector, vVector, &vVectora));
182 for (j = info.ys; j < info.ys + info.ym; j++) {
183 for (i = info.xs; i < info.xs + info.xm; i++) {
184 const PetscScalar x = (Lx * i) / M;
185 const PetscScalar y = (Ly * j) / N;
186 va[j][i] = PetscPowScalarInt(x - 0.5 * Lx, 2) + PetscPowScalarInt(y - 0.5 * Ly, 2);
187 for (c = 0; c < dof; ++c) vVectora[j][i][c] = PetscPowScalarInt(x - 0.5 * Lx, 2) + PetscPowScalarInt(y - 0.5 * Ly, 2) + 10.0 * c;
188 }
189 }
190 PetscCall(DMDAVecRestoreArray(da, v, &va));
191 PetscCall(DMDAVecRestoreArrayDOF(daVector, vVector, &vVectora));
192 PetscCall(PetscViewerVTKOpen(comm, filename, FILE_MODE_WRITE, &view));
193 PetscCall(VecView(v, view));
194 PetscCall(VecView(vVector, view));
195 PetscCall(PetscViewerDestroy(&view));
196 PetscCall(VecDestroy(&v));
197 PetscCall(VecDestroy(&vVector));
198 PetscCall(DMDestroy(&da));
199 PetscCall(DMDestroy(&daVector));
200 return PETSC_SUCCESS;
201 }
202
main(int argc,char * argv[])203 int main(int argc, char *argv[])
204 {
205 PetscInt dof;
206 PetscBool namefields;
207
208 PetscFunctionBeginUser;
209 PetscCall(PetscInitialize(&argc, &argv, 0, help));
210 dof = 2;
211 PetscCall(PetscOptionsGetInt(NULL, NULL, "-dof", &dof, NULL));
212 namefields = PETSC_FALSE;
213 PetscCall(PetscOptionsGetBool(NULL, NULL, "-namefields", &namefields, NULL));
214 PetscCall(test_3d("3d.vtr", dof, namefields));
215 PetscCall(test_2d("2d.vtr", dof, namefields));
216 PetscCall(test_3d_compat("3d_compat.vtr", dof, namefields));
217 PetscCall(test_2d_compat("2d_compat.vtr", dof, namefields));
218 PetscCall(test_3d("3d.vts", dof, namefields));
219 PetscCall(test_2d("2d.vts", dof, namefields));
220 PetscCall(test_3d_compat("3d_compat.vts", dof, namefields));
221 PetscCall(test_2d_compat("2d_compat.vts", dof, namefields));
222 PetscCall(PetscFinalize());
223 return 0;
224 }
225
226 /*TEST
227
228 build:
229 requires: !complex
230
231 test:
232 suffix: 1
233 nsize: 2
234 args: -dof 2
235 output_file: output/empty.out
236
237 test:
238 suffix: 2
239 nsize: 2
240 args: -dof 2
241 output_file: output/empty.out
242
243 test:
244 suffix: 3
245 nsize: 2
246 args: -dof 3
247 output_file: output/empty.out
248
249 test:
250 suffix: 4
251 nsize: 1
252 args: -dof 2 -namefields
253 output_file: output/empty.out
254
255 test:
256 suffix: 5
257 nsize: 2
258 args: -dof 4 -namefields
259 output_file: output/empty.out
260
261 TEST*/
262