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:   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:   DMSetFromOptions(da);
 24:   DMSetUp(da);

 26:   DMDASetUniformCoordinates(da, 0.0, Lx, 0.0, Ly, 0.0, Lz);
 27:   DMDAGetLocalInfo(da, &info);
 28:   DMCreateGlobalVector(da, &v);
 29:   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:   DMDAVecRestoreArray(da, v, &va);
 41:   PetscViewerVTKOpen(comm, filename, FILE_MODE_WRITE, &view);
 42:   VecView(v, view);
 43:   PetscViewerDestroy(&view);
 44:   VecDestroy(&v);
 45:   DMDestroy(&da);
 46:   return 0;
 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:   DMDACreate2d(comm, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_STAR, M, N, PETSC_DECIDE, PETSC_DECIDE, dof, sw, NULL, NULL, &da);
 66:   DMSetFromOptions(da);
 67:   DMSetUp(da);
 68:   DMDASetUniformCoordinates(da, 0.0, Lx, 0.0, Ly, 0.0, Lz);
 69:   DMDAGetLocalInfo(da, &info);
 70:   DMCreateGlobalVector(da, &v);
 71:   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:   DMDAVecRestoreArray(da, v, &va);
 80:   PetscViewerVTKOpen(comm, filename, FILE_MODE_WRITE, &view);
 81:   VecView(v, view);
 82:   PetscViewerDestroy(&view);
 83:   VecDestroy(&v);
 84:   DMDestroy(&da);
 85:   return 0;
 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:   DMDACreate2d(comm, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_STAR, M, N, PETSC_DECIDE, PETSC_DECIDE, dof, sw, NULL, NULL, &da);
105:   DMSetFromOptions(da);
106:   DMSetUp(da);
107:   DMDAGetLocalInfo(da, &info);
108:   DMCreateGlobalVector(da, &v);
109:   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:   DMDAVecRestoreArray(da, v, &va);
118:   PetscViewerVTKOpen(comm, filename, FILE_MODE_WRITE, &view);
119:   VecView(v, view);
120:   PetscViewerDestroy(&view);
121:   VecDestroy(&v);
122:   DMDestroy(&da);
123:   return 0;
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:   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:   DMSetFromOptions(da);
144:   DMSetUp(da);

146:   DMDAGetLocalInfo(da, &info);
147:   DMCreateGlobalVector(da, &v);
148:   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:   DMDAVecRestoreArray(da, v, &va);
160:   PetscViewerVTKOpen(comm, filename, FILE_MODE_WRITE, &view);
161:   VecView(v, view);
162:   PetscViewerDestroy(&view);
163:   VecDestroy(&v);
164:   DMDestroy(&da);
165:   return 0;
166: }

168: int main(int argc, char *argv[])
169: {
171:   PetscInitialize(&argc, &argv, 0, help);
172:   test_3d("3d.vts");
173:   test_2d("2d.vts");
174:   test_2d_nocoord("2d_nocoord.vts");
175:   test_3d_nocoord("3d_nocoord.vts");
176:   PetscFinalize();
177:   return 0;
178: }

180: /*TEST

182:    build:
183:       requires: !complex

185:    test:
186:       nsize: 2

188: TEST*/