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*/