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