Actual source code: ex47.c
1: static char help[] = "Test VTK structured grid (.vts) viewer support\n\n";
3: #include <petscdm.h>
4: #include <petscdmda.h>
6: /*
7: Write 3D DMDA vector with coordinates in VTK .vts format
9: */
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;
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));
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: }
49: /*
50: Write 2D DMDA vector with coordinates in VTK .vts format
52: */
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;
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: }
88: /*
89: Write 2D DMDA vector without coordinates in VTK .vts format
91: */
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;
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: }
126: /*
127: Write 3D DMDA vector without coordinates in VTK .vts format
129: */
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;
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));
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: }
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: }
180: /*TEST
182: build:
183: requires: !complex
185: test:
186: nsize: 2
188: TEST*/