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