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;
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: }
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: 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));
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: }
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: 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: }
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: 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));
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: }
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: 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: }
203: int main(int argc, char *argv[])
204: {
205: PetscInt dof;
206: PetscBool namefields;
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: }
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*/