Actual source code: ex26.c

  1: static char help[] = "Test FEM layout with DM and ExodusII storage\n\n";

  3: /*
  4:   In order to see the vectors which are being tested, use

  6:      -ua_vec_view -s_vec_view
  7: */

  9: #include <petsc.h>
 10: #include <exodusII.h>

 12: #include <petsc/private/dmpleximpl.h>

 14: int main(int argc, char **argv)
 15: {
 16:   DM              dm, pdm, dmU, dmA, dmS, dmUA, dmUA2, *dmList;
 17:   Vec             X, U, A, S, UA, UA2;
 18:   IS              isU, isA, isS, isUA;
 19:   PetscSection    section;
 20:   const PetscInt  fieldU     = 0;
 21:   const PetscInt  fieldA     = 2;
 22:   const PetscInt  fieldS     = 1;
 23:   const PetscInt  fieldUA[2] = {0, 2};
 24:   char            ifilename[PETSC_MAX_PATH_LEN], ofilename[PETSC_MAX_PATH_LEN];
 25:   int             exoid = -1;
 26:   IS              csIS;
 27:   const PetscInt *csID;
 28:   PetscInt       *pStartDepth, *pEndDepth;
 29:   PetscInt        order = 1;
 30:   PetscInt        sdim, d, pStart, pEnd, p, numCS, set;
 31:   PetscMPIInt     rank, size;
 32:   PetscViewer     viewer;

 34:   PetscFunctionBeginUser;
 35:   PetscCall(PetscInitialize(&argc, &argv, NULL, help));
 36:   PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
 37:   PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
 38:   PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "FEM Layout Options", "ex26");
 39:   PetscCall(PetscOptionsString("-i", "Filename to read", "ex26", ifilename, ifilename, sizeof(ifilename), NULL));
 40:   PetscCall(PetscOptionsString("-o", "Filename to write", "ex26", ofilename, ofilename, sizeof(ofilename), NULL));
 41:   PetscCall(PetscOptionsBoundedInt("-order", "FEM polynomial order", "ex26", order, &order, NULL, 1));
 42:   PetscOptionsEnd();
 43:   PetscCheck((order >= 1) && (order <= 2), PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Unsupported polynomial order %" PetscInt_FMT " not in [1, 2]", order);

 45:   /* Read the mesh from a file in any supported format */
 46:   PetscCall(DMPlexCreateFromFile(PETSC_COMM_WORLD, ifilename, NULL, PETSC_TRUE, &dm));
 47:   PetscCall(DMPlexDistributeSetDefault(dm, PETSC_FALSE));
 48:   PetscCall(DMSetFromOptions(dm));
 49:   PetscCall(DMViewFromOptions(dm, NULL, "-dm_view"));
 50:   PetscCall(DMGetDimension(dm, &sdim));

 52:   /* Create the exodus result file */
 53:   {
 54:     PetscInt numstep = 3, step;
 55:     char    *nodalVarName[4];
 56:     char    *zonalVarName[6];
 57:     int     *truthtable;
 58:     PetscInt numNodalVar, numZonalVar, i;

 60:     /* enable exodus debugging information */
 61:     ex_opts(EX_VERBOSE | EX_DEBUG);
 62:     /* Create the exodus file */
 63:     PetscCall(PetscViewerExodusIIOpen(PETSC_COMM_WORLD, ofilename, FILE_MODE_WRITE, &viewer));
 64:     /* The long way would be */
 65:     /*
 66:       PetscCall(PetscViewerCreate(PETSC_COMM_WORLD,&viewer));
 67:       PetscCall(PetscViewerSetType(viewer,PETSCVIEWEREXODUSII));
 68:       PetscCall(PetscViewerFileSetMode(viewer,FILE_MODE_APPEND));
 69:       PetscCall(PetscViewerFileSetName(viewer,ofilename));
 70:     */
 71:     /* set the mesh order */
 72:     PetscCall(PetscViewerExodusIISetOrder(viewer, order));
 73:     PetscCall(PetscViewerView(viewer, PETSC_VIEWER_STDOUT_WORLD));
 74:     /*
 75:       Notice how the exodus file is actually NOT open at this point (exoid is -1)
 76:       Since we are overwriting the file (mode is FILE_MODE_WRITE), we are going to have to
 77:       write the geometry (the DM), which can only be done on a brand new file.
 78:     */

 80:     /* Save the geometry to the file, erasing all previous content */
 81:     PetscCall(DMView(dm, viewer));
 82:     PetscCall(PetscViewerView(viewer, PETSC_VIEWER_STDOUT_WORLD));
 83:     /*
 84:       Note how the exodus file is now open
 85:     */

 87:     /* "Format" the exodus result file, i.e. allocate space for nodal and zonal variables */
 88:     switch (sdim) {
 89:     case 2:
 90:       numNodalVar     = 3;
 91:       nodalVarName[0] = (char *)"U_x";
 92:       nodalVarName[1] = (char *)"U_y";
 93:       nodalVarName[2] = (char *)"Alpha";
 94:       numZonalVar     = 3;
 95:       zonalVarName[0] = (char *)"Sigma_11";
 96:       zonalVarName[1] = (char *)"Sigma_22";
 97:       zonalVarName[2] = (char *)"Sigma_12";
 98:       break;
 99:     case 3:
100:       numNodalVar     = 4;
101:       nodalVarName[0] = (char *)"U_x";
102:       nodalVarName[1] = (char *)"U_y";
103:       nodalVarName[2] = (char *)"U_z";
104:       nodalVarName[3] = (char *)"Alpha";
105:       numZonalVar     = 6;
106:       zonalVarName[0] = (char *)"Sigma_11";
107:       zonalVarName[1] = (char *)"Sigma_22";
108:       zonalVarName[2] = (char *)"Sigma_33";
109:       zonalVarName[3] = (char *)"Sigma_23";
110:       zonalVarName[4] = (char *)"Sigma_13";
111:       zonalVarName[5] = (char *)"Sigma_12";
112:       break;
113:     default:
114:       SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "No layout for dimension %" PetscInt_FMT, sdim);
115:     }
116:     PetscCall(PetscViewerExodusIIGetId(viewer, &exoid));
117:     PetscCallExternal(ex_put_variable_param, exoid, EX_ELEM_BLOCK, numZonalVar);
118:     PetscCallExternal(ex_put_variable_names, exoid, EX_ELEM_BLOCK, numZonalVar, zonalVarName);
119:     PetscCallExternal(ex_put_variable_param, exoid, EX_NODAL, numNodalVar);
120:     PetscCallExternal(ex_put_variable_names, exoid, EX_NODAL, numNodalVar, nodalVarName);
121:     numCS = ex_inquire_int(exoid, EX_INQ_ELEM_BLK);

123:     /*
124:       An exodusII truth table specifies which fields are saved at which time step
125:       It speeds up I/O but reserving space for fieldsin the file ahead of time.
126:     */
127:     PetscCall(PetscMalloc1(numZonalVar * numCS, &truthtable));
128:     for (i = 0; i < numZonalVar * numCS; ++i) truthtable[i] = 1;
129:     PetscCallExternal(ex_put_truth_table, exoid, EX_ELEM_BLOCK, numCS, numZonalVar, truthtable);
130:     PetscCall(PetscFree(truthtable));

132:     /* Writing time step information in the file. Note that this is currently broken in the exodus library for netcdf4 (HDF5-based) files */
133:     for (step = 0; step < numstep; ++step) {
134:       PetscReal time = step;
135:       PetscCallExternal(ex_put_time, exoid, step + 1, &time);
136:     }
137:   }

139:   /* Create the main section containing all fields */
140:   PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)dm), &section));
141:   PetscCall(PetscSectionSetNumFields(section, 3));
142:   PetscCall(PetscSectionSetFieldName(section, fieldU, "U"));
143:   PetscCall(PetscSectionSetFieldName(section, fieldA, "Alpha"));
144:   PetscCall(PetscSectionSetFieldName(section, fieldS, "Sigma"));
145:   PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
146:   PetscCall(PetscSectionSetChart(section, pStart, pEnd));
147:   PetscCall(PetscMalloc2(sdim + 1, &pStartDepth, sdim + 1, &pEndDepth));
148:   for (d = 0; d <= sdim; ++d) PetscCall(DMPlexGetDepthStratum(dm, d, &pStartDepth[d], &pEndDepth[d]));
149:   /* Vector field U, Scalar field Alpha, Tensor field Sigma */
150:   PetscCall(PetscSectionSetFieldComponents(section, fieldU, sdim));
151:   PetscCall(PetscSectionSetFieldComponents(section, fieldA, 1));
152:   PetscCall(PetscSectionSetFieldComponents(section, fieldS, sdim * (sdim + 1) / 2));

154:   /* Going through cell sets then cells, and setting up storage for the sections */
155:   PetscCall(DMGetLabelSize(dm, "Cell Sets", &numCS));
156:   PetscCall(DMGetLabelIdIS(dm, "Cell Sets", &csIS));
157:   if (csIS) PetscCall(ISGetIndices(csIS, &csID));
158:   for (set = 0; set < numCS; set++) {
159:     IS              cellIS;
160:     const PetscInt *cellID;
161:     PetscInt        numCells, cell, closureSize, *closureA = NULL;

163:     PetscCall(DMGetStratumSize(dm, "Cell Sets", csID[set], &numCells));
164:     PetscCall(DMGetStratumIS(dm, "Cell Sets", csID[set], &cellIS));
165:     if (numCells > 0) {
166:       /* dof layout ordered by increasing height in the DAG: cell, face, edge, vertex */
167:       PetscInt  dofUP1Tri[]  = {2, 0, 0};
168:       PetscInt  dofAP1Tri[]  = {1, 0, 0};
169:       PetscInt  dofUP2Tri[]  = {2, 2, 0};
170:       PetscInt  dofAP2Tri[]  = {1, 1, 0};
171:       PetscInt  dofUP1Quad[] = {2, 0, 0};
172:       PetscInt  dofAP1Quad[] = {1, 0, 0};
173:       PetscInt  dofUP2Quad[] = {2, 2, 2};
174:       PetscInt  dofAP2Quad[] = {1, 1, 1};
175:       PetscInt  dofS2D[]     = {0, 0, 3};
176:       PetscInt  dofUP1Tet[]  = {3, 0, 0, 0};
177:       PetscInt  dofAP1Tet[]  = {1, 0, 0, 0};
178:       PetscInt  dofUP2Tet[]  = {3, 3, 0, 0};
179:       PetscInt  dofAP2Tet[]  = {1, 1, 0, 0};
180:       PetscInt  dofUP1Hex[]  = {3, 0, 0, 0};
181:       PetscInt  dofAP1Hex[]  = {1, 0, 0, 0};
182:       PetscInt  dofUP2Hex[]  = {3, 3, 3, 3};
183:       PetscInt  dofAP2Hex[]  = {1, 1, 1, 1};
184:       PetscInt  dofS3D[]     = {0, 0, 0, 6};
185:       PetscInt *dofU, *dofA, *dofS;

187:       switch (sdim) {
188:       case 2:
189:         dofS = dofS2D;
190:         break;
191:       case 3:
192:         dofS = dofS3D;
193:         break;
194:       default:
195:         SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "No layout for dimension %" PetscInt_FMT, sdim);
196:       }

198:       /* Identify cell type based on closure size only. This works for Tri/Tet/Quad/Hex meshes
199:          It will not be enough to identify more exotic elements like pyramid or prisms...  */
200:       PetscCall(ISGetIndices(cellIS, &cellID));
201:       PetscCall(DMPlexGetTransitiveClosure(dm, cellID[0], PETSC_TRUE, &closureSize, &closureA));
202:       switch (closureSize) {
203:       case 7: /* Tri */
204:         if (order == 1) {
205:           dofU = dofUP1Tri;
206:           dofA = dofAP1Tri;
207:         } else {
208:           dofU = dofUP2Tri;
209:           dofA = dofAP2Tri;
210:         }
211:         break;
212:       case 9: /* Quad */
213:         if (order == 1) {
214:           dofU = dofUP1Quad;
215:           dofA = dofAP1Quad;
216:         } else {
217:           dofU = dofUP2Quad;
218:           dofA = dofAP2Quad;
219:         }
220:         break;
221:       case 15: /* Tet */
222:         if (order == 1) {
223:           dofU = dofUP1Tet;
224:           dofA = dofAP1Tet;
225:         } else {
226:           dofU = dofUP2Tet;
227:           dofA = dofAP2Tet;
228:         }
229:         break;
230:       case 27: /* Hex */
231:         if (order == 1) {
232:           dofU = dofUP1Hex;
233:           dofA = dofAP1Hex;
234:         } else {
235:           dofU = dofUP2Hex;
236:           dofA = dofAP2Hex;
237:         }
238:         break;
239:       default:
240:         SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Unknown element with closure size %" PetscInt_FMT, closureSize);
241:       }
242:       PetscCall(DMPlexRestoreTransitiveClosure(dm, cellID[0], PETSC_TRUE, &closureSize, &closureA));

244:       for (cell = 0; cell < numCells; cell++) {
245:         PetscInt *closure = NULL;

247:         PetscCall(DMPlexGetTransitiveClosure(dm, cellID[cell], PETSC_TRUE, &closureSize, &closure));
248:         for (p = 0; p < closureSize; ++p) {
249:           /* Find depth of p */
250:           for (d = 0; d <= sdim; ++d) {
251:             if ((closure[2 * p] >= pStartDepth[d]) && (closure[2 * p] < pEndDepth[d])) {
252:               PetscCall(PetscSectionSetDof(section, closure[2 * p], dofU[d] + dofA[d] + dofS[d]));
253:               PetscCall(PetscSectionSetFieldDof(section, closure[2 * p], fieldU, dofU[d]));
254:               PetscCall(PetscSectionSetFieldDof(section, closure[2 * p], fieldA, dofA[d]));
255:               PetscCall(PetscSectionSetFieldDof(section, closure[2 * p], fieldS, dofS[d]));
256:             }
257:           }
258:         }
259:         PetscCall(DMPlexRestoreTransitiveClosure(dm, cellID[cell], PETSC_TRUE, &closureSize, &closure));
260:       }
261:       PetscCall(ISRestoreIndices(cellIS, &cellID));
262:       PetscCall(ISDestroy(&cellIS));
263:     }
264:   }
265:   if (csIS) PetscCall(ISRestoreIndices(csIS, &csID));
266:   PetscCall(ISDestroy(&csIS));
267:   PetscCall(PetscSectionSetUp(section));
268:   PetscCall(DMSetLocalSection(dm, section));
269:   PetscCall(PetscObjectViewFromOptions((PetscObject)section, NULL, "-dm_section_view"));
270:   PetscCall(PetscSectionDestroy(&section));

272:   {
273:     PetscSF          migrationSF;
274:     PetscInt         ovlp = 0;
275:     PetscPartitioner part;

277:     PetscCall(DMSetUseNatural(dm, PETSC_TRUE));
278:     PetscCall(DMPlexGetPartitioner(dm, &part));
279:     PetscCall(PetscPartitionerSetFromOptions(part));
280:     PetscCall(DMPlexDistribute(dm, ovlp, &migrationSF, &pdm));
281:     if (!pdm) pdm = dm;
282:     /* Set the migrationSF is mandatory since useNatural = PETSC_TRUE */
283:     if (migrationSF) {
284:       PetscCall(DMPlexSetMigrationSF(pdm, migrationSF));
285:       PetscCall(PetscSFDestroy(&migrationSF));
286:     }
287:     PetscCall(DMViewFromOptions(pdm, NULL, "-dm_view"));
288:   }

290:   /* Get DM and IS for each field of dm */
291:   PetscCall(DMCreateSubDM(pdm, 1, &fieldU, &isU, &dmU));
292:   PetscCall(DMCreateSubDM(pdm, 1, &fieldA, &isA, &dmA));
293:   PetscCall(DMCreateSubDM(pdm, 1, &fieldS, &isS, &dmS));
294:   PetscCall(DMCreateSubDM(pdm, 2, fieldUA, &isUA, &dmUA));

296:   PetscCall(PetscMalloc1(2, &dmList));
297:   dmList[0] = dmU;
298:   dmList[1] = dmA;

300:   PetscCall(DMCreateSuperDM(dmList, 2, NULL, &dmUA2));
301:   PetscCall(PetscFree(dmList));

303:   PetscCall(DMGetGlobalVector(pdm, &X));
304:   PetscCall(DMGetGlobalVector(dmU, &U));
305:   PetscCall(DMGetGlobalVector(dmA, &A));
306:   PetscCall(DMGetGlobalVector(dmS, &S));
307:   PetscCall(DMGetGlobalVector(dmUA, &UA));
308:   PetscCall(DMGetGlobalVector(dmUA2, &UA2));

310:   PetscCall(PetscObjectSetName((PetscObject)U, "U"));
311:   PetscCall(PetscObjectSetName((PetscObject)A, "Alpha"));
312:   PetscCall(PetscObjectSetName((PetscObject)S, "Sigma"));
313:   PetscCall(PetscObjectSetName((PetscObject)UA, "UAlpha"));
314:   PetscCall(PetscObjectSetName((PetscObject)UA2, "UAlpha2"));
315:   PetscCall(VecSet(X, -111.));

317:   /* Setting u to [x,y,z]  and alpha to x^2+y^2+z^2 by writing in UAlpha then restricting to U and Alpha */
318:   {
319:     PetscSection sectionUA;
320:     Vec          UALoc;
321:     PetscSection coordSection;
322:     Vec          coord;
323:     PetscScalar *cval, *xyz;
324:     PetscInt     clSize, i, j;

326:     PetscCall(DMGetLocalSection(dmUA, &sectionUA));
327:     PetscCall(DMGetLocalVector(dmUA, &UALoc));
328:     PetscCall(VecGetArray(UALoc, &cval));
329:     PetscCall(DMGetCoordinateSection(dmUA, &coordSection));
330:     PetscCall(DMGetCoordinatesLocal(dmUA, &coord));
331:     PetscCall(DMPlexGetChart(dmUA, &pStart, &pEnd));

333:     for (p = pStart; p < pEnd; ++p) {
334:       PetscInt dofUA, offUA;

336:       PetscCall(PetscSectionGetDof(sectionUA, p, &dofUA));
337:       if (dofUA > 0) {
338:         xyz = NULL;
339:         PetscCall(PetscSectionGetOffset(sectionUA, p, &offUA));
340:         PetscCall(DMPlexVecGetClosure(dmUA, coordSection, coord, p, &clSize, &xyz));
341:         cval[offUA + sdim] = 0.;
342:         for (i = 0; i < sdim; ++i) {
343:           cval[offUA + i] = 0;
344:           for (j = 0; j < clSize / sdim; ++j) cval[offUA + i] += xyz[j * sdim + i];
345:           cval[offUA + i] = cval[offUA + i] * sdim / clSize;
346:           cval[offUA + sdim] += PetscSqr(cval[offUA + i]);
347:         }
348:         PetscCall(DMPlexVecRestoreClosure(dmUA, coordSection, coord, p, &clSize, &xyz));
349:       }
350:     }
351:     PetscCall(VecRestoreArray(UALoc, &cval));
352:     PetscCall(DMLocalToGlobalBegin(dmUA, UALoc, INSERT_VALUES, UA));
353:     PetscCall(DMLocalToGlobalEnd(dmUA, UALoc, INSERT_VALUES, UA));
354:     PetscCall(DMRestoreLocalVector(dmUA, &UALoc));

356:     /* Update X */
357:     PetscCall(VecISCopy(X, isUA, SCATTER_FORWARD, UA));
358:     PetscCall(VecViewFromOptions(UA, NULL, "-ua_vec_view"));

360:     /* Restrict to U and Alpha */
361:     PetscCall(VecISCopy(X, isU, SCATTER_REVERSE, U));
362:     PetscCall(VecISCopy(X, isA, SCATTER_REVERSE, A));

364:     /* restrict to UA2 */
365:     PetscCall(VecISCopy(X, isUA, SCATTER_REVERSE, UA2));
366:     PetscCall(VecViewFromOptions(UA2, NULL, "-ua2_vec_view"));
367:   }

369:   {
370:     Vec          tmpVec;
371:     PetscSection coordSection;
372:     Vec          coord;
373:     PetscReal    norm;
374:     PetscReal    time = 1.234;

376:     /* Writing nodal variables to ExodusII file */
377:     PetscCall(DMSetOutputSequenceNumber(dmU, 0, time));
378:     PetscCall(DMSetOutputSequenceNumber(dmA, 0, time));

380:     PetscCall(VecView(U, viewer));
381:     PetscCall(VecView(A, viewer));

383:     /* Saving U and Alpha in one shot.
384:        For this, we need to cheat and change the Vec's name
385:        Note that in the end we write variables one component at a time,
386:        so that there is no real values in doing this
387:     */

389:     PetscCall(DMSetOutputSequenceNumber(dmUA, 1, time));
390:     PetscCall(DMGetGlobalVector(dmUA, &tmpVec));
391:     PetscCall(VecCopy(UA, tmpVec));
392:     PetscCall(PetscObjectSetName((PetscObject)tmpVec, "U"));
393:     PetscCall(VecView(tmpVec, viewer));
394:     /* Reading nodal variables in Exodus file */
395:     PetscCall(VecSet(tmpVec, -1000.0));
396:     PetscCall(VecLoad(tmpVec, viewer));
397:     PetscCall(VecAXPY(UA, -1.0, tmpVec));
398:     PetscCall(VecNorm(UA, NORM_INFINITY, &norm));
399:     PetscCheck(norm <= PETSC_SQRT_MACHINE_EPSILON, PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "UAlpha ||Vin - Vout|| = %g", (double)norm);
400:     PetscCall(DMRestoreGlobalVector(dmUA, &tmpVec));

402:     /* same thing with the UA2 Vec obtained from the superDM */
403:     PetscCall(DMGetGlobalVector(dmUA2, &tmpVec));
404:     PetscCall(VecCopy(UA2, tmpVec));
405:     PetscCall(PetscObjectSetName((PetscObject)tmpVec, "U"));
406:     PetscCall(DMSetOutputSequenceNumber(dmUA2, 2, time));
407:     PetscCall(VecView(tmpVec, viewer));
408:     /* Reading nodal variables in Exodus file */
409:     PetscCall(VecSet(tmpVec, -1000.0));
410:     PetscCall(VecLoad(tmpVec, viewer));
411:     PetscCall(VecAXPY(UA2, -1.0, tmpVec));
412:     PetscCall(VecNorm(UA2, NORM_INFINITY, &norm));
413:     PetscCheck(norm <= PETSC_SQRT_MACHINE_EPSILON, PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "UAlpha2 ||Vin - Vout|| = %g", (double)norm);
414:     PetscCall(DMRestoreGlobalVector(dmUA2, &tmpVec));

416:     /* Building and saving Sigma
417:        We set sigma_0 = rank (to see partitioning)
418:               sigma_1 = cell set ID
419:               sigma_2 = x_coordinate of the cell center of mass
420:     */
421:     PetscCall(DMGetCoordinateSection(dmS, &coordSection));
422:     PetscCall(DMGetCoordinatesLocal(dmS, &coord));
423:     PetscCall(DMGetLabelIdIS(dmS, "Cell Sets", &csIS));
424:     PetscCall(DMGetLabelSize(dmS, "Cell Sets", &numCS));
425:     PetscCall(ISGetIndices(csIS, &csID));
426:     for (set = 0; set < numCS; ++set) {
427:       /* We know that all cells in a cell set have the same type, so we can dimension cval and xyz once for each cell set */
428:       IS              cellIS;
429:       const PetscInt *cellID;
430:       PetscInt        numCells, cell;
431:       PetscScalar    *cval = NULL, *xyz = NULL;
432:       PetscInt        clSize, cdimCoord, c;

434:       PetscCall(DMGetStratumIS(dmS, "Cell Sets", csID[set], &cellIS));
435:       PetscCall(ISGetIndices(cellIS, &cellID));
436:       PetscCall(ISGetSize(cellIS, &numCells));
437:       for (cell = 0; cell < numCells; cell++) {
438:         PetscCall(DMPlexVecGetClosure(dmS, NULL, S, cellID[cell], &clSize, &cval));
439:         PetscCall(DMPlexVecGetClosure(dmS, coordSection, coord, cellID[cell], &cdimCoord, &xyz));
440:         cval[0] = rank;
441:         cval[1] = csID[set];
442:         cval[2] = 0.;
443:         for (c = 0; c < cdimCoord / sdim; c++) cval[2] += xyz[c * sdim];
444:         cval[2] = cval[2] * sdim / cdimCoord;
445:         PetscCall(DMPlexVecSetClosure(dmS, NULL, S, cellID[cell], cval, INSERT_ALL_VALUES));
446:       }
447:       PetscCall(DMPlexVecRestoreClosure(dmS, NULL, S, cellID[0], &clSize, &cval));
448:       PetscCall(DMPlexVecRestoreClosure(dmS, coordSection, coord, cellID[0], NULL, &xyz));
449:       PetscCall(ISRestoreIndices(cellIS, &cellID));
450:       PetscCall(ISDestroy(&cellIS));
451:     }
452:     PetscCall(ISRestoreIndices(csIS, &csID));
453:     PetscCall(ISDestroy(&csIS));
454:     PetscCall(VecViewFromOptions(S, NULL, "-s_vec_view"));

456:     /* Writing zonal variables in Exodus file */
457:     PetscCall(DMSetOutputSequenceNumber(dmS, 0, time));
458:     PetscCall(VecView(S, viewer));

460:     /* Reading zonal variables in Exodus file */
461:     PetscCall(DMGetGlobalVector(dmS, &tmpVec));
462:     PetscCall(VecSet(tmpVec, -1000.0));
463:     PetscCall(PetscObjectSetName((PetscObject)tmpVec, "Sigma"));
464:     PetscCall(VecLoad(tmpVec, viewer));
465:     PetscCall(VecAXPY(S, -1.0, tmpVec));
466:     PetscCall(VecNorm(S, NORM_INFINITY, &norm));
467:     PetscCheck(norm <= PETSC_SQRT_MACHINE_EPSILON, PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Sigma ||Vin - Vout|| = %g", (double)norm);
468:     PetscCall(DMRestoreGlobalVector(dmS, &tmpVec));
469:   }
470:   PetscCall(PetscViewerDestroy(&viewer));

472:   PetscCall(DMRestoreGlobalVector(dmUA2, &UA2));
473:   PetscCall(DMRestoreGlobalVector(dmUA, &UA));
474:   PetscCall(DMRestoreGlobalVector(dmS, &S));
475:   PetscCall(DMRestoreGlobalVector(dmA, &A));
476:   PetscCall(DMRestoreGlobalVector(dmU, &U));
477:   PetscCall(DMRestoreGlobalVector(pdm, &X));
478:   PetscCall(DMDestroy(&dmU));
479:   PetscCall(ISDestroy(&isU));
480:   PetscCall(DMDestroy(&dmA));
481:   PetscCall(ISDestroy(&isA));
482:   PetscCall(DMDestroy(&dmS));
483:   PetscCall(ISDestroy(&isS));
484:   PetscCall(DMDestroy(&dmUA));
485:   PetscCall(ISDestroy(&isUA));
486:   PetscCall(DMDestroy(&dmUA2));
487:   PetscCall(DMDestroy(&pdm));
488:   if (size > 1) PetscCall(DMDestroy(&dm));
489:   PetscCall(PetscFree2(pStartDepth, pEndDepth));
490:   PetscCall(PetscFinalize());
491:   return 0;
492: }

494: /*TEST

496:   build:
497:     requires: exodusii pnetcdf !complex
498:   # 2D seq
499:   test:
500:     suffix: 0
501:     args: -i ${wPETSC_DIR}/share/petsc/datafiles/meshes/FourSquareT-large.exo -o FourSquareT-large_out.exo -dm_view -petscpartitioner_type simple -order 1
502:     #TODO: bug in call to NetCDF failed to complete invalid type definition in file id 65536 NetCDF: One or more variable sizes violate format constraints
503:   test:
504:     suffix: 1
505:     args: -i ${wPETSC_DIR}/share/petsc/datafiles/meshes/FourSquareQ-large.exo -o FourSquareQ-large_out.exo -dm_view -petscpartitioner_type simple -order 1
506:   test:
507:     suffix: 2
508:     args: -i ${wPETSC_DIR}/share/petsc/datafiles/meshes/FourSquareH-large.exo -o FourSquareH-large_out.exo -dm_view -petscpartitioner_type simple -order 1
509:     #TODO: bug in call to NetCDF failed to complete invalid type definition in file id 65536 NetCDF: One or more variable sizes violate format constraints
510:   test:
511:     suffix: 3
512:     args: -i ${wPETSC_DIR}/share/petsc/datafiles/meshes/FourSquareT-large.exo -o FourSquareT-large_out.exo -dm_view -petscpartitioner_type simple -order 2
513:   test:
514:     suffix: 4
515:     args: -i ${wPETSC_DIR}/share/petsc/datafiles/meshes/FourSquareQ-large.exo -o FourSquareQ-large_out.exo -dm_view -petscpartitioner_type simple -order 2
516:   test:
517:     suffix: 5
518:     args: -i ${wPETSC_DIR}/share/petsc/datafiles/meshes/FourSquareH-large.exo -o FourSquareH-large_out.exo -dm_view -petscpartitioner_type simple -order 2

520:   # 2D par
521:   test:
522:     suffix: 6
523:     nsize: 2
524:     args: -i ${wPETSC_DIR}/share/petsc/datafiles/meshes/FourSquareT-large.exo -o FourSquareT-large_out.exo -dm_view -petscpartitioner_type simple -order 1
525:     #TODO: bug in call to NetCDF failed to complete invalid type definition in file id 65536 NetCDF: One or more variable sizes violate format constraints
526:   test:
527:     suffix: 7
528:     nsize: 2
529:     args: -i ${wPETSC_DIR}/share/petsc/datafiles/meshes/FourSquareQ-large.exo -o FourSquareQ-large_out.exo -dm_view -petscpartitioner_type simple -order 1
530:   test:
531:     suffix: 8
532:     nsize: 2
533:     args: -i ${wPETSC_DIR}/share/petsc/datafiles/meshes/FourSquareH-large.exo -o FourSquareH-large_out.exo -dm_view -petscpartitioner_type simple -order 1
534:     #TODO: bug in call to NetCDF failed to complete invalid type definition in file id 65536 NetCDF: invalid dimension ID or name
535:   test:
536:     suffix: 9
537:     nsize: 2
538:     args: -i ${wPETSC_DIR}/share/petsc/datafiles/meshes/FourSquareT-large.exo -o FourSquareT-large_out.exo -dm_view -petscpartitioner_type simple -order 2
539:   test:
540:     suffix: 10
541:     nsize: 2
542:     args: -i ${wPETSC_DIR}/share/petsc/datafiles/meshes/FourSquareQ-large.exo -o FourSquareQ-large_out.exo -dm_view -petscpartitioner_type simple -order 2
543:   test:
544:     # Something is now broken with parallel read/write for whatever shape H is
545:     TODO: broken
546:     suffix: 11
547:     nsize: 2
548:     args: -i ${wPETSC_DIR}/share/petsc/datafiles/meshes/FourSquareH-large.exo -o FourSquareH-large_out.exo -dm_view -petscpartitioner_type simple -order 2

550:   #3d seq
551:   test:
552:     suffix: 12
553:     args: -i ${wPETSC_DIR}/share/petsc/datafiles/meshes/FourBrickHex-large.exo -o FourBrickHex-large_out.exo -dm_view -petscpartitioner_type simple -order 1
554:   test:
555:     suffix: 13
556:     args: -i ${wPETSC_DIR}/share/petsc/datafiles/meshes/FourBrickTet-large.exo -o FourBrickTet-large_out.exo -dm_view -petscpartitioner_type simple -order 1
557:     #TODO: bug in call to NetCDF failed to complete invalid type definition in file id 65536 NetCDF: One or more variable sizes violate format constraints
558:   test:
559:     suffix: 14
560:     args: -i ${wPETSC_DIR}/share/petsc/datafiles/meshes/FourBrickHex-large.exo -o FourBrickHex-large_out.exo -dm_view -petscpartitioner_type simple -order 2
561:   test:
562:     suffix: 15
563:     args: -i ${wPETSC_DIR}/share/petsc/datafiles/meshes/FourBrickTet-large.exo -o FourBrickTet-large_out.exo -dm_view -petscpartitioner_type simple -order 2
564:     #TODO: bug in call to NetCDF failed to complete invalid type definition in file id 65536 NetCDF: One or more variable sizes violate format constraints
565:   #3d par
566:   test:
567:     suffix: 16
568:     nsize: 2
569:     args: -i ${wPETSC_DIR}/share/petsc/datafiles/meshes/FourBrickHex-large.exo -o FourBrickHex-large_out.exo -dm_view -petscpartitioner_type simple -order 1
570:   test:
571:     suffix: 17
572:     nsize: 2
573:     args: -i ${wPETSC_DIR}/share/petsc/datafiles/meshes/FourBrickTet-large.exo -o FourBrickTet-large_out.exo -dm_view -petscpartitioner_type simple -order 1
574:     #TODO: bug in call to NetCDF failed to complete invalid type definition in file id 65536 NetCDF: One or more variable sizes violate format constraints
575:   test:
576:     suffix: 18
577:     nsize: 2
578:     args: -i ${wPETSC_DIR}/share/petsc/datafiles/meshes/FourBrickHex-large.exo -o FourBrickHex-large_out.exo -dm_view -petscpartitioner_type simple -order 2
579:   test:
580:     suffix: 19
581:     nsize: 2
582:     args: -i ${wPETSC_DIR}/share/petsc/datafiles/meshes/FourBrickTet-large.exo -o FourBrickTet-large_out.exo -dm_view -petscpartitioner_type simple -order 2
583:     #TODO: bug in call to NetCDF failed to complete invalid type definition in file id 65536 NetCDF: One or more variable sizes violate format constraints

585: TEST*/