Actual source code: ex42.c
1: /* -*- Mode: C++; c-basic-offset:2 ; indent-tabs-mode:nil ; -*- */
3: static char help[] = "Test VTK Rectilinear grid (.vtr) viewer support\n\n";
5: #include <petscdm.h>
6: #include <petscdmda.h>
8: /*
9: Write 3D DMDA vector with coordinates in VTK VTR format
11: */
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;
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));
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: }
51: /*
52: Write 2D DMDA vector with coordinates in VTK VTR format
54: */
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;
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: }
90: /*
91: Write 2D DMDA vector without coordinates in VTK VTR format
93: */
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;
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: }
128: /*
129: Write 3D DMDA vector without coordinates in VTK VTR format
131: */
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;
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));
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: }
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: }
182: /*TEST
184: build:
185: requires: !complex
187: test:
188: nsize: 2
190: TEST*/