Actual source code: ex48.c
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";
4: #include <petscdm.h>
5: #include <petscdmda.h>
7: /* Helper function to name DMDA fields */
8: PetscErrorCode NameFields(DM da, PetscInt dof)
9: {
10: PetscInt c;
13: for (c = 0; c < dof; ++c) {
14: char fieldname[256];
15: PetscSNPrintf(fieldname, sizeof(fieldname), "field_%" PetscInt_FMT, c);
16: DMDASetFieldName(da, c, fieldname);
17: }
18: return 0;
19: }
21: /*
22: Write 3D DMDA vector with coordinates in VTK format
23: */
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;
36: 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: DMSetFromOptions(da);
38: DMSetUp(da);
39: if (namefields) NameFields(da, dof);
41: DMDASetUniformCoordinates(da, 0.0, Lx, 0.0, Ly, 0.0, Lz);
42: DMDAGetLocalInfo(da, &info);
43: DMCreateGlobalVector(da, &v);
44: 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: DMDAVecRestoreArrayDOF(da, v, &va);
56: PetscViewerVTKOpen(comm, filename, FILE_MODE_WRITE, &view);
57: VecView(v, view);
58: PetscViewerDestroy(&view);
59: VecDestroy(&v);
60: DMDestroy(&da);
61: return 0;
62: }
64: /*
65: Write 2D DMDA vector with coordinates in VTK format
66: */
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;
79: DMDACreate2d(comm, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_STAR, M, N, PETSC_DECIDE, PETSC_DECIDE, dof, sw, NULL, NULL, &da);
80: DMSetFromOptions(da);
81: DMSetUp(da);
82: if (namefields) NameFields(da, dof);
83: DMDASetUniformCoordinates(da, 0.0, Lx, 0.0, Ly, 0.0, Lz);
84: DMDAGetLocalInfo(da, &info);
85: DMCreateGlobalVector(da, &v);
86: 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: DMDAVecRestoreArrayDOF(da, v, &va);
95: PetscViewerVTKOpen(comm, filename, FILE_MODE_WRITE, &view);
96: VecView(v, view);
97: PetscViewerDestroy(&view);
98: VecDestroy(&v);
99: DMDestroy(&da);
100: return 0;
101: }
103: /*
104: Write a scalar and a vector field from two compatible 3d DMDAs
105: */
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;
118: 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: DMSetFromOptions(da);
120: DMSetUp(da);
121: if (namefields) NameFields(da, 1);
123: DMDASetUniformCoordinates(da, 0.0, Lx, 0.0, Ly, 0.0, Lz);
124: DMDAGetLocalInfo(da, &info);
125: DMDACreateCompatibleDMDA(da, dof, &daVector);
126: if (namefields) NameFields(daVector, dof);
127: DMCreateGlobalVector(da, &v);
128: DMCreateGlobalVector(daVector, &vVector);
129: DMDAVecGetArray(da, v, &va);
130: 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: DMDAVecRestoreArray(da, v, &va);
143: DMDAVecRestoreArrayDOF(da, v, &vVectora);
144: PetscViewerVTKOpen(comm, filename, FILE_MODE_WRITE, &view);
145: VecView(v, view);
146: VecView(vVector, view);
147: PetscViewerDestroy(&view);
148: VecDestroy(&v);
149: VecDestroy(&vVector);
150: DMDestroy(&da);
151: DMDestroy(&daVector);
152: return 0;
153: }
155: /*
156: Write a scalar and a vector field from two compatible 2d DMDAs
157: */
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;
170: DMDACreate2d(comm, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_STAR, M, N, PETSC_DECIDE, PETSC_DECIDE, /* dof:*/ 1, sw, NULL, NULL, &da);
171: DMSetFromOptions(da);
172: DMSetUp(da);
173: if (namefields) NameFields(da, 1);
174: DMDASetUniformCoordinates(da, 0.0, Lx, 0.0, Ly, 0.0, Lz);
175: DMDACreateCompatibleDMDA(da, dof, &daVector);
176: if (namefields) NameFields(daVector, dof);
177: DMDAGetLocalInfo(da, &info);
178: DMCreateGlobalVector(da, &v);
179: DMCreateGlobalVector(daVector, &vVector);
180: DMDAVecGetArray(da, v, &va);
181: 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: DMDAVecRestoreArray(da, v, &va);
191: DMDAVecRestoreArrayDOF(daVector, vVector, &vVectora);
192: PetscViewerVTKOpen(comm, filename, FILE_MODE_WRITE, &view);
193: VecView(v, view);
194: VecView(vVector, view);
195: PetscViewerDestroy(&view);
196: VecDestroy(&v);
197: VecDestroy(&vVector);
198: DMDestroy(&da);
199: DMDestroy(&daVector);
200: return 0;
201: }
203: int main(int argc, char *argv[])
204: {
205: PetscInt dof;
206: PetscBool namefields;
209: PetscInitialize(&argc, &argv, 0, help);
210: dof = 2;
211: PetscOptionsGetInt(NULL, NULL, "-dof", &dof, NULL);
212: namefields = PETSC_FALSE;
213: PetscOptionsGetBool(NULL, NULL, "-namefields", &namefields, NULL);
214: test_3d("3d.vtr", dof, namefields);
215: test_2d("2d.vtr", dof, namefields);
216: test_3d_compat("3d_compat.vtr", dof, namefields);
217: test_2d_compat("2d_compat.vtr", dof, namefields);
218: test_3d("3d.vts", dof, namefields);
219: test_2d("2d.vts", dof, namefields);
220: test_3d_compat("3d_compat.vts", dof, namefields);
221: test_2d_compat("2d_compat.vts", dof, namefields);
222: PetscFinalize();
223: return 0;
224: }
226: /*TEST
228: build:
229: requires: !complex
231: test:
232: suffix: 1
233: nsize: 2
234: args: -dof 2
236: test:
237: suffix: 2
238: nsize: 2
239: args: -dof 2
241: test:
242: suffix: 3
243: nsize: 2
244: args: -dof 3
246: test:
247: suffix: 4
248: nsize: 1
249: args: -dof 2 -namefields
251: test:
252: suffix: 5
253: nsize: 2
254: args: -dof 4 -namefields
256: TEST*/