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;

 13:   for (c = 0; c < dof; ++c) {
 14:     char fieldname[256];
 15:     PetscSNPrintf(fieldname, sizeof(fieldname), "field_%" PetscInt_FMT, c);
 16:     DMDASetFieldName(da, c, fieldname);
 17:   }
 18:   return 0;
 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:   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:   DMSetFromOptions(da);
 38:   DMSetUp(da);
 39:   if (namefields) NameFields(da, dof);

 41:   DMDASetUniformCoordinates(da, 0.0, Lx, 0.0, Ly, 0.0, Lz);
 42:   DMDAGetLocalInfo(da, &info);
 43:   DMCreateGlobalVector(da, &v);
 44:   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:   DMDAVecRestoreArrayDOF(da, v, &va);
 56:   PetscViewerVTKOpen(comm, filename, FILE_MODE_WRITE, &view);
 57:   VecView(v, view);
 58:   PetscViewerDestroy(&view);
 59:   VecDestroy(&v);
 60:   DMDestroy(&da);
 61:   return 0;
 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:   DMDACreate2d(comm, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_STAR, M, N, PETSC_DECIDE, PETSC_DECIDE, dof, sw, NULL, NULL, &da);
 80:   DMSetFromOptions(da);
 81:   DMSetUp(da);
 82:   if (namefields) NameFields(da, dof);
 83:   DMDASetUniformCoordinates(da, 0.0, Lx, 0.0, Ly, 0.0, Lz);
 84:   DMDAGetLocalInfo(da, &info);
 85:   DMCreateGlobalVector(da, &v);
 86:   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:   DMDAVecRestoreArrayDOF(da, v, &va);
 95:   PetscViewerVTKOpen(comm, filename, FILE_MODE_WRITE, &view);
 96:   VecView(v, view);
 97:   PetscViewerDestroy(&view);
 98:   VecDestroy(&v);
 99:   DMDestroy(&da);
100:   return 0;
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:   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:   DMSetFromOptions(da);
120:   DMSetUp(da);
121:   if (namefields) NameFields(da, 1);

123:   DMDASetUniformCoordinates(da, 0.0, Lx, 0.0, Ly, 0.0, Lz);
124:   DMDAGetLocalInfo(da, &info);
125:   DMDACreateCompatibleDMDA(da, dof, &daVector);
126:   if (namefields) NameFields(daVector, dof);
127:   DMCreateGlobalVector(da, &v);
128:   DMCreateGlobalVector(daVector, &vVector);
129:   DMDAVecGetArray(da, v, &va);
130:   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:   DMDAVecRestoreArray(da, v, &va);
143:   DMDAVecRestoreArrayDOF(da, v, &vVectora);
144:   PetscViewerVTKOpen(comm, filename, FILE_MODE_WRITE, &view);
145:   VecView(v, view);
146:   VecView(vVector, view);
147:   PetscViewerDestroy(&view);
148:   VecDestroy(&v);
149:   VecDestroy(&vVector);
150:   DMDestroy(&da);
151:   DMDestroy(&daVector);
152:   return 0;
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:   DMDACreate2d(comm, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_STAR, M, N, PETSC_DECIDE, PETSC_DECIDE, /* dof:*/ 1, sw, NULL, NULL, &da);
171:   DMSetFromOptions(da);
172:   DMSetUp(da);
173:   if (namefields) NameFields(da, 1);
174:   DMDASetUniformCoordinates(da, 0.0, Lx, 0.0, Ly, 0.0, Lz);
175:   DMDACreateCompatibleDMDA(da, dof, &daVector);
176:   if (namefields) NameFields(daVector, dof);
177:   DMDAGetLocalInfo(da, &info);
178:   DMCreateGlobalVector(da, &v);
179:   DMCreateGlobalVector(daVector, &vVector);
180:   DMDAVecGetArray(da, v, &va);
181:   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:   DMDAVecRestoreArray(da, v, &va);
191:   DMDAVecRestoreArrayDOF(daVector, vVector, &vVectora);
192:   PetscViewerVTKOpen(comm, filename, FILE_MODE_WRITE, &view);
193:   VecView(v, view);
194:   VecView(vVector, view);
195:   PetscViewerDestroy(&view);
196:   VecDestroy(&v);
197:   VecDestroy(&vVector);
198:   DMDestroy(&da);
199:   DMDestroy(&daVector);
200:   return 0;
201: }

203: int main(int argc, char *argv[])
204: {
205:   PetscInt  dof;
206:   PetscBool namefields;

209:   PetscInitialize(&argc, &argv, 0, help);
210:   dof = 2;
211:   PetscOptionsGetInt(NULL, NULL, "-dof", &dof, NULL);
212:   namefields = PETSC_FALSE;
213:   PetscOptionsGetBool(NULL, NULL, "-namefields", &namefields, NULL);
214:   test_3d("3d.vtr", dof, namefields);
215:   test_2d("2d.vtr", dof, namefields);
216:   test_3d_compat("3d_compat.vtr", dof, namefields);
217:   test_2d_compat("2d_compat.vtr", dof, namefields);
218:   test_3d("3d.vts", dof, namefields);
219:   test_2d("2d.vts", dof, namefields);
220:   test_3d_compat("3d_compat.vts", dof, namefields);
221:   test_2d_compat("2d_compat.vts", dof, namefields);
222:   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*/