Actual source code: plexexodusii.c

  1: #define PETSCDM_DLL
  2: #include <petsc/private/dmpleximpl.h>

  4: #if defined(PETSC_HAVE_EXODUSII)
  5:   #include <netcdf.h>
  6:   #include <exodusII.h>
  7: #endif

  9: #include <petsc/private/viewerimpl.h>
 10: #include <petsc/private/viewerexodusiiimpl.h>
 11: #if defined(PETSC_HAVE_EXODUSII)
 12: /*@C
 13:   PETSC_VIEWER_EXODUSII_ - Creates an `PETSCVIEWEREXODUSII` `PetscViewer` shared by all processors in a communicator.

 15:   Collective

 17:   Input Parameter:
 18: . comm - the MPI communicator to share the `PETSCVIEWEREXODUSII` `PetscViewer`

 20:   Level: intermediate

 22:   Note:
 23:   Unlike almost all other PETSc routines, PETSC_VIEWER_EXODUSII_ does not return
 24:   an error code.  The GLVIS PetscViewer is usually used in the form
 25: $       XXXView(XXX object, PETSC_VIEWER_EXODUSII_(comm));

 27:   Fortran Note:
 28:   No support in Fortran

 30: .seealso: `PETSCVIEWEREXODUSII`, `PetscViewer`, `PetscViewer`, `PetscViewerExodusIIOpen()`, `PetscViewerType`, `PetscViewerCreate()`, `PetscViewerDestroy()`
 31: @*/
 32: PetscViewer PETSC_VIEWER_EXODUSII_(MPI_Comm comm)
 33: {
 34:   PetscViewer    viewer;

 37:   PetscViewerExodusIIOpen(comm, "mesh.exo", FILE_MODE_WRITE, &viewer);
 38:   if (ierr) {
 39:     PetscError(PETSC_COMM_SELF, __LINE__, "PETSC_VIEWER_EXODUSII_", __FILE__, PETSC_ERR_PLIB, PETSC_ERROR_INITIAL, " ");
 40:     return NULL;
 41:   }
 42:   PetscObjectRegisterDestroy((PetscObject)viewer);
 43:   if (ierr) {
 44:     PetscError(PETSC_COMM_SELF, __LINE__, "PETSC_VIEWER_EXODUSII_", __FILE__, PETSC_ERR_PLIB, PETSC_ERROR_INITIAL, " ");
 45:     return NULL;
 46:   }
 47:   return viewer;
 48: }

 50: static PetscErrorCode PetscViewerView_ExodusII(PetscViewer v, PetscViewer viewer)
 51: {
 52:   PetscViewer_ExodusII *exo = (PetscViewer_ExodusII *)v->data;

 54:   if (exo->filename) PetscViewerASCIIPrintf(viewer, "Filename:    %s\n", exo->filename);
 55:   if (exo->exoid) PetscViewerASCIIPrintf(viewer, "exoid:       %d\n", exo->exoid);
 56:   if (exo->btype) PetscViewerASCIIPrintf(viewer, "IO Mode:     %d\n", exo->btype);
 57:   if (exo->order) PetscViewerASCIIPrintf(viewer, "Mesh order:  %" PetscInt_FMT "\n", exo->order);
 58:   return 0;
 59: }

 61: static PetscErrorCode PetscViewerSetFromOptions_ExodusII(PetscViewer v, PetscOptionItems *PetscOptionsObject)
 62: {
 63:   PetscOptionsHeadBegin(PetscOptionsObject, "ExodusII PetscViewer Options");
 64:   PetscOptionsHeadEnd();
 65:   return 0;
 66: }

 68: static PetscErrorCode PetscViewerSetUp_ExodusII(PetscViewer viewer)
 69: {
 70:   return 0;
 71: }

 73: static PetscErrorCode PetscViewerDestroy_ExodusII(PetscViewer viewer)
 74: {
 75:   PetscViewer_ExodusII *exo = (PetscViewer_ExodusII *)viewer->data;

 77:   if (exo->exoid >= 0) PetscCallExternal(ex_close, exo->exoid);
 78:   PetscFree(exo->filename);
 79:   PetscFree(exo);
 80:   PetscObjectComposeFunction((PetscObject)viewer, "PetscViewerFileSetName_C", NULL);
 81:   PetscObjectComposeFunction((PetscObject)viewer, "PetscViewerFileGetName_C", NULL);
 82:   PetscObjectComposeFunction((PetscObject)viewer, "PetscViewerFileSetMode_C", NULL);
 83:   PetscObjectComposeFunction((PetscObject)viewer, "PetscViewerFileGetMode_C", NULL);
 84:   PetscObjectComposeFunction((PetscObject)viewer, "PetscViewerGetId_C", NULL);
 85:   PetscObjectComposeFunction((PetscObject)viewer, "PetscViewerGetOrder_C", NULL);
 86:   PetscObjectComposeFunction((PetscObject)viewer, "PetscViewerSetOrder_C", NULL);
 87:   return 0;
 88: }

 90: static PetscErrorCode PetscViewerFileSetName_ExodusII(PetscViewer viewer, const char name[])
 91: {
 92:   PetscViewer_ExodusII *exo = (PetscViewer_ExodusII *)viewer->data;
 93:   PetscMPIInt           rank;
 94:   int                   CPU_word_size, IO_word_size, EXO_mode;
 95:   MPI_Info              mpi_info = MPI_INFO_NULL;
 96:   float                 EXO_version;

 98:   MPI_Comm_rank(PetscObjectComm((PetscObject)viewer), &rank);
 99:   CPU_word_size = sizeof(PetscReal);
100:   IO_word_size  = sizeof(PetscReal);

102:   if (exo->exoid >= 0) {
103:     PetscCallExternal(ex_close, exo->exoid);
104:     exo->exoid = -1;
105:   }
106:   if (exo->filename) PetscFree(exo->filename);
107:   PetscStrallocpy(name, &exo->filename);
108:   switch (exo->btype) {
109:   case FILE_MODE_READ:
110:     EXO_mode = EX_READ;
111:     break;
112:   case FILE_MODE_APPEND:
113:   case FILE_MODE_UPDATE:
114:   case FILE_MODE_APPEND_UPDATE:
115:     /* Will fail if the file does not already exist */
116:     EXO_mode = EX_WRITE;
117:     break;
118:   case FILE_MODE_WRITE:
119:     /*
120:       exodus only allows writing geometry upon file creation, so we will let DMView create the file.
121:     */
122:     return 0;
123:     break;
124:   default:
125:     SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ORDER, "Must call PetscViewerFileSetMode() before PetscViewerFileSetName()");
126:   }
127:   #if defined(PETSC_USE_64BIT_INDICES)
128:   EXO_mode += EX_ALL_INT64_API;
129:   #endif
130:   exo->exoid = ex_open_par(name, EXO_mode, &CPU_word_size, &IO_word_size, &EXO_version, PETSC_COMM_WORLD, mpi_info);
132:   return 0;
133: }

135: static PetscErrorCode PetscViewerFileGetName_ExodusII(PetscViewer viewer, const char **name)
136: {
137:   PetscViewer_ExodusII *exo = (PetscViewer_ExodusII *)viewer->data;

139:   *name = exo->filename;
140:   return 0;
141: }

143: static PetscErrorCode PetscViewerFileSetMode_ExodusII(PetscViewer viewer, PetscFileMode type)
144: {
145:   PetscViewer_ExodusII *exo = (PetscViewer_ExodusII *)viewer->data;

147:   exo->btype = type;
148:   return 0;
149: }

151: static PetscErrorCode PetscViewerFileGetMode_ExodusII(PetscViewer viewer, PetscFileMode *type)
152: {
153:   PetscViewer_ExodusII *exo = (PetscViewer_ExodusII *)viewer->data;

155:   *type = exo->btype;
156:   return 0;
157: }

159: static PetscErrorCode PetscViewerExodusIIGetId_ExodusII(PetscViewer viewer, int *exoid)
160: {
161:   PetscViewer_ExodusII *exo = (PetscViewer_ExodusII *)viewer->data;

163:   *exoid = exo->exoid;
164:   return 0;
165: }

167: static PetscErrorCode PetscViewerExodusIIGetOrder_ExodusII(PetscViewer viewer, PetscInt *order)
168: {
169:   PetscViewer_ExodusII *exo = (PetscViewer_ExodusII *)viewer->data;

171:   *order = exo->order;
172:   return 0;
173: }

175: static PetscErrorCode PetscViewerExodusIISetOrder_ExodusII(PetscViewer viewer, PetscInt order)
176: {
177:   PetscViewer_ExodusII *exo = (PetscViewer_ExodusII *)viewer->data;

179:   exo->order = order;
180:   return 0;
181: }

183: /*MC
184:    PETSCVIEWEREXODUSII - A viewer that writes to an Exodus II file

186: .seealso: `PetscViewerExodusIIOpen()`, `PetscViewerCreate()`, `PETSCVIEWERBINARY`, `PETSCVIEWERHDF5`, `DMView()`,
187:           `PetscViewerFileSetName()`, `PetscViewerFileSetMode()`, `PetscViewerFormat`, `PetscViewerType`, `PetscViewerSetType()`

189:   Level: beginner
190: M*/

192: PETSC_EXTERN PetscErrorCode PetscViewerCreate_ExodusII(PetscViewer v)
193: {
194:   PetscViewer_ExodusII *exo;

196:   PetscNew(&exo);

198:   v->data                = (void *)exo;
199:   v->ops->destroy        = PetscViewerDestroy_ExodusII;
200:   v->ops->setfromoptions = PetscViewerSetFromOptions_ExodusII;
201:   v->ops->setup          = PetscViewerSetUp_ExodusII;
202:   v->ops->view           = PetscViewerView_ExodusII;
203:   v->ops->flush          = 0;
204:   exo->btype             = (PetscFileMode)-1;
205:   exo->filename          = 0;
206:   exo->exoid             = -1;

208:   PetscObjectComposeFunction((PetscObject)v, "PetscViewerFileSetName_C", PetscViewerFileSetName_ExodusII);
209:   PetscObjectComposeFunction((PetscObject)v, "PetscViewerFileGetName_C", PetscViewerFileGetName_ExodusII);
210:   PetscObjectComposeFunction((PetscObject)v, "PetscViewerFileSetMode_C", PetscViewerFileSetMode_ExodusII);
211:   PetscObjectComposeFunction((PetscObject)v, "PetscViewerFileGetMode_C", PetscViewerFileGetMode_ExodusII);
212:   PetscObjectComposeFunction((PetscObject)v, "PetscViewerGetId_C", PetscViewerExodusIIGetId_ExodusII);
213:   PetscObjectComposeFunction((PetscObject)v, "PetscViewerSetOrder_C", PetscViewerExodusIISetOrder_ExodusII);
214:   PetscObjectComposeFunction((PetscObject)v, "PetscViewerGetOrder_C", PetscViewerExodusIIGetOrder_ExodusII);
215:   return 0;
216: }

218: /*
219:   EXOGetVarIndex - Locate a result in an exodus file based on its name

221:   Collective

223:   Input Parameters:
224: + exoid    - the exodus id of a file (obtained from ex_open or ex_create for instance)
225: . obj_type - the type of entity for instance EX_NODAL, EX_ELEM_BLOCK
226: - name     - the name of the result

228:   Output Parameters:
229: . varIndex - the location in the exodus file of the result

231:   Notes:
232:   The exodus variable index is obtained by comparing name and the
233:   names of zonal variables declared in the exodus file. For instance if name is "V"
234:   the location in the exodus file will be the first match of "V", "V_X", "V_XX", "V_1", or "V_11"
235:   amongst all variables of type obj_type.

237:   Level: beginner

239: .seealso: `EXOGetVarIndex()`, `DMPlexView_ExodusII_Internal()`, `VecViewPlex_ExodusII_Nodal_Internal()`, `VecLoadNodal_PlexEXO()`, `VecLoadZonal_PlexEXO()`
240: */
241: PetscErrorCode EXOGetVarIndex_Internal(int exoid, ex_entity_type obj_type, const char name[], int *varIndex)
242: {
243:   int       num_vars, i, j;
244:   char      ext_name[MAX_STR_LENGTH + 1], var_name[MAX_STR_LENGTH + 1];
245:   const int num_suffix = 5;
246:   char     *suffix[5];
247:   PetscBool flg;

249:   suffix[0] = (char *)"";
250:   suffix[1] = (char *)"_X";
251:   suffix[2] = (char *)"_XX";
252:   suffix[3] = (char *)"_1";
253:   suffix[4] = (char *)"_11";

255:   *varIndex = -1;
256:   PetscCallExternal(ex_get_variable_param, exoid, obj_type, &num_vars);
257:   for (i = 0; i < num_vars; ++i) {
258:     PetscCallExternal(ex_get_variable_name, exoid, obj_type, i + 1, var_name);
259:     for (j = 0; j < num_suffix; ++j) {
260:       PetscStrncpy(ext_name, name, MAX_STR_LENGTH);
261:       PetscStrlcat(ext_name, suffix[j], MAX_STR_LENGTH);
262:       PetscStrcasecmp(ext_name, var_name, &flg);
263:       if (flg) *varIndex = i + 1;
264:     }
265:   }
266:   return 0;
267: }

269: /*
270:   DMView_PlexExodusII - Write a DM to disk in exodus format

272:   Collective on dm

274:   Input Parameters:
275: + dm  - The dm to be written
276: . viewer - an exodusII viewer

278:   Notes:
279:   Not all DM can be written to disk this way. For instance, exodus assume that element blocks (mapped to "Cell sets" labels)
280:   consists of sequentially numbered cells. If this is not the case, the exodus file will be corrupted.

282:   If the dm has been distributed, only the part of the DM on MPI rank 0 (including "ghost" cells and vertices)
283:   will be written.

285:   DMPlex only represents geometry while most post-processing software expect that a mesh also provides information
286:   on the discretization space. This function assumes that the file represents Lagrange finite elements of order 1 or 2.
287:   The order of the mesh shall be set using PetscViewerExodusIISetOrder
288:   It should be extended to use PetscFE objects.

290:   This function will only handle TRI, TET, QUAD, and HEX cells.
291:   Level: beginner

293: .seealso:
294: */
295: PetscErrorCode DMView_PlexExodusII(DM dm, PetscViewer viewer)
296: {
297:   enum ElemType {
298:     TRI,
299:     QUAD,
300:     TET,
301:     HEX
302:   };
303:   MPI_Comm comm;
304:   PetscInt degree; /* the order of the mesh */
305:   /* Connectivity Variables */
306:   PetscInt cellsNotInConnectivity;
307:   /* Cell Sets */
308:   DMLabel         csLabel;
309:   IS              csIS;
310:   const PetscInt *csIdx;
311:   PetscInt        num_cs, cs;
312:   enum ElemType  *type;
313:   PetscBool       hasLabel;
314:   /* Coordinate Variables */
315:   DM           cdm;
316:   PetscSection coordSection;
317:   Vec          coord;
318:   PetscInt   **nodes;
319:   PetscInt     depth, d, dim, skipCells = 0;
320:   PetscInt     pStart, pEnd, p, cStart, cEnd, numCells, vStart, vEnd, numVertices, eStart, eEnd, numEdges, fStart, fEnd, numFaces, numNodes;
321:   PetscInt     num_vs, num_fs;
322:   PetscMPIInt  rank, size;
323:   const char  *dmName;
324:   PetscInt     nodesTriP1[4]  = {3, 0, 0, 0};
325:   PetscInt     nodesTriP2[4]  = {3, 3, 0, 0};
326:   PetscInt     nodesQuadP1[4] = {4, 0, 0, 0};
327:   PetscInt     nodesQuadP2[4] = {4, 4, 0, 1};
328:   PetscInt     nodesTetP1[4]  = {4, 0, 0, 0};
329:   PetscInt     nodesTetP2[4]  = {4, 6, 0, 0};
330:   PetscInt     nodesHexP1[4]  = {8, 0, 0, 0};
331:   PetscInt     nodesHexP2[4]  = {8, 12, 6, 1};
332:   int          CPU_word_size, IO_word_size, EXO_mode;
333:   float        EXO_version;

335:   PetscViewer_ExodusII *exo = (PetscViewer_ExodusII *)viewer->data;

337:   PetscObjectGetComm((PetscObject)dm, &comm);
338:   MPI_Comm_rank(comm, &rank);
339:   MPI_Comm_size(comm, &size);

341:   /*
342:     Creating coordSection is a collective operation so we do it somewhat out of sequence
343:   */
344:   PetscSectionCreate(comm, &coordSection);
345:   DMGetCoordinatesLocalSetUp(dm);
346:   /*
347:     Check that all points are on rank 0 since we don't know how to save distributed DM in exodus format
348:   */
349:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
350:   DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);
351:   DMPlexGetDepthStratum(dm, 1, &eStart, &eEnd);
352:   DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
353:   numCells    = cEnd - cStart;
354:   numEdges    = eEnd - eStart;
355:   numVertices = vEnd - vStart;
357:   if (rank == 0) {
358:     switch (exo->btype) {
359:     case FILE_MODE_READ:
360:     case FILE_MODE_APPEND:
361:     case FILE_MODE_UPDATE:
362:     case FILE_MODE_APPEND_UPDATE:
363:       /* exodusII does not allow writing geometry to an existing file */
364:       SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "cannot add geometry to existing file %s", exo->filename);
365:     case FILE_MODE_WRITE:
366:       /* Create an empty file if one already exists*/
367:       EXO_mode = EX_CLOBBER;
368:   #if defined(PETSC_USE_64BIT_INDICES)
369:       EXO_mode += EX_ALL_INT64_API;
370:   #endif
371:       CPU_word_size = sizeof(PetscReal);
372:       IO_word_size  = sizeof(PetscReal);
373:       exo->exoid    = ex_create(exo->filename, EXO_mode, &CPU_word_size, &IO_word_size);

376:       break;
377:     default:
378:       SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ORDER, "Must call PetscViewerFileSetMode() before PetscViewerFileSetName()");
379:     }

381:     /* --- Get DM info --- */
382:     PetscObjectGetName((PetscObject)dm, &dmName);
383:     DMPlexGetDepth(dm, &depth);
384:     DMGetDimension(dm, &dim);
385:     DMPlexGetChart(dm, &pStart, &pEnd);
386:     if (depth == 3) {
387:       numFaces = fEnd - fStart;
388:     } else {
389:       numFaces = 0;
390:     }
391:     DMGetLabelSize(dm, "Cell Sets", &num_cs);
392:     DMGetLabelSize(dm, "Vertex Sets", &num_vs);
393:     DMGetLabelSize(dm, "Face Sets", &num_fs);
394:     DMGetCoordinatesLocal(dm, &coord);
395:     DMGetCoordinateDM(dm, &cdm);
396:     if (num_cs > 0) {
397:       DMGetLabel(dm, "Cell Sets", &csLabel);
398:       DMLabelGetValueIS(csLabel, &csIS);
399:       ISGetIndices(csIS, &csIdx);
400:     }
401:     PetscMalloc1(num_cs, &nodes);
402:     /* Set element type for each block and compute total number of nodes */
403:     PetscMalloc1(num_cs, &type);
404:     numNodes = numVertices;

406:     PetscViewerExodusIIGetOrder(viewer, &degree);
407:     if (degree == 2) numNodes += numEdges;
408:     cellsNotInConnectivity = numCells;
409:     for (cs = 0; cs < num_cs; ++cs) {
410:       IS              stratumIS;
411:       const PetscInt *cells;
412:       PetscScalar    *xyz = NULL;
413:       PetscInt        csSize, closureSize;

415:       DMLabelGetStratumIS(csLabel, csIdx[cs], &stratumIS);
416:       ISGetIndices(stratumIS, &cells);
417:       ISGetSize(stratumIS, &csSize);
418:       DMPlexVecGetClosure(cdm, NULL, coord, cells[0], &closureSize, &xyz);
419:       switch (dim) {
420:       case 2:
421:         if (closureSize == 3 * dim) {
422:           type[cs] = TRI;
423:         } else if (closureSize == 4 * dim) {
424:           type[cs] = QUAD;
425:         } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Number of vertices %" PetscInt_FMT " in dimension %" PetscInt_FMT " has no ExodusII type", closureSize / dim, dim);
426:         break;
427:       case 3:
428:         if (closureSize == 4 * dim) {
429:           type[cs] = TET;
430:         } else if (closureSize == 8 * dim) {
431:           type[cs] = HEX;
432:         } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Number of vertices %" PetscInt_FMT " in dimension %" PetscInt_FMT " has no ExodusII type", closureSize / dim, dim);
433:         break;
434:       default:
435:         SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Dimension %" PetscInt_FMT " not handled by ExodusII viewer", dim);
436:       }
437:       if ((degree == 2) && (type[cs] == QUAD)) numNodes += csSize;
438:       if ((degree == 2) && (type[cs] == HEX)) {
439:         numNodes += csSize;
440:         numNodes += numFaces;
441:       }
442:       DMPlexVecRestoreClosure(cdm, NULL, coord, cells[0], &closureSize, &xyz);
443:       /* Set nodes and Element type */
444:       if (type[cs] == TRI) {
445:         if (degree == 1) nodes[cs] = nodesTriP1;
446:         else if (degree == 2) nodes[cs] = nodesTriP2;
447:       } else if (type[cs] == QUAD) {
448:         if (degree == 1) nodes[cs] = nodesQuadP1;
449:         else if (degree == 2) nodes[cs] = nodesQuadP2;
450:       } else if (type[cs] == TET) {
451:         if (degree == 1) nodes[cs] = nodesTetP1;
452:         else if (degree == 2) nodes[cs] = nodesTetP2;
453:       } else if (type[cs] == HEX) {
454:         if (degree == 1) nodes[cs] = nodesHexP1;
455:         else if (degree == 2) nodes[cs] = nodesHexP2;
456:       }
457:       /* Compute the number of cells not in the connectivity table */
458:       cellsNotInConnectivity -= nodes[cs][3] * csSize;

460:       ISRestoreIndices(stratumIS, &cells);
461:       ISDestroy(&stratumIS);
462:     }
463:     if (num_cs) PetscCallExternal(ex_put_init, exo->exoid, dmName, dim, numNodes, numCells, num_cs, num_vs, num_fs);
464:     /* --- Connectivity --- */
465:     for (cs = 0; cs < num_cs; ++cs) {
466:       IS              stratumIS;
467:       const PetscInt *cells;
468:       PetscInt       *connect, off = 0;
469:       PetscInt        edgesInClosure = 0, facesInClosure = 0, verticesInClosure = 0;
470:       PetscInt        csSize, c, connectSize, closureSize;
471:       char           *elem_type        = NULL;
472:       char            elem_type_tri3[] = "TRI3", elem_type_quad4[] = "QUAD4";
473:       char            elem_type_tri6[] = "TRI6", elem_type_quad9[] = "QUAD9";
474:       char            elem_type_tet4[] = "TET4", elem_type_hex8[] = "HEX8";
475:       char            elem_type_tet10[] = "TET10", elem_type_hex27[] = "HEX27";

477:       DMLabelGetStratumIS(csLabel, csIdx[cs], &stratumIS);
478:       ISGetIndices(stratumIS, &cells);
479:       ISGetSize(stratumIS, &csSize);
480:       /* Set Element type */
481:       if (type[cs] == TRI) {
482:         if (degree == 1) elem_type = elem_type_tri3;
483:         else if (degree == 2) elem_type = elem_type_tri6;
484:       } else if (type[cs] == QUAD) {
485:         if (degree == 1) elem_type = elem_type_quad4;
486:         else if (degree == 2) elem_type = elem_type_quad9;
487:       } else if (type[cs] == TET) {
488:         if (degree == 1) elem_type = elem_type_tet4;
489:         else if (degree == 2) elem_type = elem_type_tet10;
490:       } else if (type[cs] == HEX) {
491:         if (degree == 1) elem_type = elem_type_hex8;
492:         else if (degree == 2) elem_type = elem_type_hex27;
493:       }
494:       connectSize = nodes[cs][0] + nodes[cs][1] + nodes[cs][2] + nodes[cs][3];
495:       PetscMalloc1(PetscMax(27, connectSize) * csSize, &connect);
496:       PetscCallExternal(ex_put_block, exo->exoid, EX_ELEM_BLOCK, csIdx[cs], elem_type, csSize, connectSize, 0, 0, 1);
497:       /* Find number of vertices, edges, and faces in the closure */
498:       verticesInClosure = nodes[cs][0];
499:       if (depth > 1) {
500:         if (dim == 2) {
501:           DMPlexGetConeSize(dm, cells[0], &edgesInClosure);
502:         } else if (dim == 3) {
503:           PetscInt *closure = NULL;

505:           DMPlexGetConeSize(dm, cells[0], &facesInClosure);
506:           DMPlexGetTransitiveClosure(dm, cells[0], PETSC_TRUE, &closureSize, &closure);
507:           edgesInClosure = closureSize - facesInClosure - 1 - verticesInClosure;
508:           DMPlexRestoreTransitiveClosure(dm, cells[0], PETSC_TRUE, &closureSize, &closure);
509:         }
510:       }
511:       /* Get connectivity for each cell */
512:       for (c = 0; c < csSize; ++c) {
513:         PetscInt *closure = NULL;
514:         PetscInt  temp, i;

516:         DMPlexGetTransitiveClosure(dm, cells[c], PETSC_TRUE, &closureSize, &closure);
517:         for (i = 0; i < connectSize; ++i) {
518:           if (i < nodes[cs][0]) { /* Vertices */
519:             connect[i + off] = closure[(i + edgesInClosure + facesInClosure + 1) * 2] + 1;
520:             connect[i + off] -= cellsNotInConnectivity;
521:           } else if (i < nodes[cs][0] + nodes[cs][1]) { /* Edges */
522:             connect[i + off] = closure[(i - verticesInClosure + facesInClosure + 1) * 2] + 1;
523:             if (nodes[cs][2] == 0) connect[i + off] -= numFaces;
524:             connect[i + off] -= cellsNotInConnectivity;
525:           } else if (i < nodes[cs][0] + nodes[cs][1] + nodes[cs][3]) { /* Cells */
526:             connect[i + off] = closure[0] + 1;
527:             connect[i + off] -= skipCells;
528:           } else if (i < nodes[cs][0] + nodes[cs][1] + nodes[cs][3] + nodes[cs][2]) { /* Faces */
529:             connect[i + off] = closure[(i - edgesInClosure - verticesInClosure) * 2] + 1;
530:             connect[i + off] -= cellsNotInConnectivity;
531:           } else {
532:             connect[i + off] = -1;
533:           }
534:         }
535:         /* Tetrahedra are inverted */
536:         if (type[cs] == TET) {
537:           temp             = connect[0 + off];
538:           connect[0 + off] = connect[1 + off];
539:           connect[1 + off] = temp;
540:           if (degree == 2) {
541:             temp             = connect[5 + off];
542:             connect[5 + off] = connect[6 + off];
543:             connect[6 + off] = temp;
544:             temp             = connect[7 + off];
545:             connect[7 + off] = connect[8 + off];
546:             connect[8 + off] = temp;
547:           }
548:         }
549:         /* Hexahedra are inverted */
550:         if (type[cs] == HEX) {
551:           temp             = connect[1 + off];
552:           connect[1 + off] = connect[3 + off];
553:           connect[3 + off] = temp;
554:           if (degree == 2) {
555:             temp              = connect[8 + off];
556:             connect[8 + off]  = connect[11 + off];
557:             connect[11 + off] = temp;
558:             temp              = connect[9 + off];
559:             connect[9 + off]  = connect[10 + off];
560:             connect[10 + off] = temp;
561:             temp              = connect[16 + off];
562:             connect[16 + off] = connect[17 + off];
563:             connect[17 + off] = temp;
564:             temp              = connect[18 + off];
565:             connect[18 + off] = connect[19 + off];
566:             connect[19 + off] = temp;

568:             temp              = connect[12 + off];
569:             connect[12 + off] = connect[16 + off];
570:             connect[16 + off] = temp;
571:             temp              = connect[13 + off];
572:             connect[13 + off] = connect[17 + off];
573:             connect[17 + off] = temp;
574:             temp              = connect[14 + off];
575:             connect[14 + off] = connect[18 + off];
576:             connect[18 + off] = temp;
577:             temp              = connect[15 + off];
578:             connect[15 + off] = connect[19 + off];
579:             connect[19 + off] = temp;

581:             temp              = connect[23 + off];
582:             connect[23 + off] = connect[26 + off];
583:             connect[26 + off] = temp;
584:             temp              = connect[24 + off];
585:             connect[24 + off] = connect[25 + off];
586:             connect[25 + off] = temp;
587:             temp              = connect[25 + off];
588:             connect[25 + off] = connect[26 + off];
589:             connect[26 + off] = temp;
590:           }
591:         }
592:         off += connectSize;
593:         DMPlexRestoreTransitiveClosure(dm, cells[c], PETSC_TRUE, &closureSize, &closure);
594:       }
595:       PetscCallExternal(ex_put_conn, exo->exoid, EX_ELEM_BLOCK, csIdx[cs], connect, 0, 0);
596:       skipCells += (nodes[cs][3] == 0) * csSize;
597:       PetscFree(connect);
598:       ISRestoreIndices(stratumIS, &cells);
599:       ISDestroy(&stratumIS);
600:     }
601:     PetscFree(type);
602:     /* --- Coordinates --- */
603:     PetscSectionSetChart(coordSection, pStart, pEnd);
604:     if (num_cs) {
605:       for (d = 0; d < depth; ++d) {
606:         DMPlexGetDepthStratum(dm, d, &pStart, &pEnd);
607:         for (p = pStart; p < pEnd; ++p) PetscSectionSetDof(coordSection, p, nodes[0][d] > 0);
608:       }
609:     }
610:     for (cs = 0; cs < num_cs; ++cs) {
611:       IS              stratumIS;
612:       const PetscInt *cells;
613:       PetscInt        csSize, c;

615:       DMLabelGetStratumIS(csLabel, csIdx[cs], &stratumIS);
616:       ISGetIndices(stratumIS, &cells);
617:       ISGetSize(stratumIS, &csSize);
618:       for (c = 0; c < csSize; ++c) PetscSectionSetDof(coordSection, cells[c], nodes[cs][3] > 0);
619:       ISRestoreIndices(stratumIS, &cells);
620:       ISDestroy(&stratumIS);
621:     }
622:     if (num_cs) {
623:       ISRestoreIndices(csIS, &csIdx);
624:       ISDestroy(&csIS);
625:     }
626:     PetscFree(nodes);
627:     PetscSectionSetUp(coordSection);
628:     if (numNodes) {
629:       const char  *coordNames[3] = {"x", "y", "z"};
630:       PetscScalar *closure, *cval;
631:       PetscReal   *coords;
632:       PetscInt     hasDof, n = 0;

634:       /* There can't be more than 24 values in the closure of a point for the coord coordSection */
635:       PetscCalloc3(numNodes * 3, &coords, dim, &cval, 24, &closure);
636:       DMGetCoordinatesLocalNoncollective(dm, &coord);
637:       DMPlexGetChart(dm, &pStart, &pEnd);
638:       for (p = pStart; p < pEnd; ++p) {
639:         PetscSectionGetDof(coordSection, p, &hasDof);
640:         if (hasDof) {
641:           PetscInt closureSize = 24, j;

643:           DMPlexVecGetClosure(cdm, NULL, coord, p, &closureSize, &closure);
644:           for (d = 0; d < dim; ++d) {
645:             cval[d] = 0.0;
646:             for (j = 0; j < closureSize / dim; j++) cval[d] += closure[j * dim + d];
647:             coords[d * numNodes + n] = PetscRealPart(cval[d]) * dim / closureSize;
648:           }
649:           ++n;
650:         }
651:       }
652:       PetscCallExternal(ex_put_coord, exo->exoid, &coords[0 * numNodes], &coords[1 * numNodes], &coords[2 * numNodes]);
653:       PetscFree3(coords, cval, closure);
654:       PetscCallExternal(ex_put_coord_names, exo->exoid, (char **)coordNames);
655:     }

657:     /* --- Node Sets/Vertex Sets --- */
658:     DMHasLabel(dm, "Vertex Sets", &hasLabel);
659:     if (hasLabel) {
660:       PetscInt        i, vs, vsSize;
661:       const PetscInt *vsIdx, *vertices;
662:       PetscInt       *nodeList;
663:       IS              vsIS, stratumIS;
664:       DMLabel         vsLabel;
665:       DMGetLabel(dm, "Vertex Sets", &vsLabel);
666:       DMLabelGetValueIS(vsLabel, &vsIS);
667:       ISGetIndices(vsIS, &vsIdx);
668:       for (vs = 0; vs < num_vs; ++vs) {
669:         DMLabelGetStratumIS(vsLabel, vsIdx[vs], &stratumIS);
670:         ISGetIndices(stratumIS, &vertices);
671:         ISGetSize(stratumIS, &vsSize);
672:         PetscMalloc1(vsSize, &nodeList);
673:         for (i = 0; i < vsSize; ++i) nodeList[i] = vertices[i] - skipCells + 1;
674:         PetscCallExternal(ex_put_set_param, exo->exoid, EX_NODE_SET, vsIdx[vs], vsSize, 0);
675:         PetscCallExternal(ex_put_set, exo->exoid, EX_NODE_SET, vsIdx[vs], nodeList, NULL);
676:         ISRestoreIndices(stratumIS, &vertices);
677:         ISDestroy(&stratumIS);
678:         PetscFree(nodeList);
679:       }
680:       ISRestoreIndices(vsIS, &vsIdx);
681:       ISDestroy(&vsIS);
682:     }
683:     /* --- Side Sets/Face Sets --- */
684:     DMHasLabel(dm, "Face Sets", &hasLabel);
685:     if (hasLabel) {
686:       PetscInt        i, j, fs, fsSize;
687:       const PetscInt *fsIdx, *faces;
688:       IS              fsIS, stratumIS;
689:       DMLabel         fsLabel;
690:       PetscInt        numPoints, *points;
691:       PetscInt        elem_list_size = 0;
692:       PetscInt       *elem_list, *elem_ind, *side_list;

694:       DMGetLabel(dm, "Face Sets", &fsLabel);
695:       /* Compute size of Node List and Element List */
696:       DMLabelGetValueIS(fsLabel, &fsIS);
697:       ISGetIndices(fsIS, &fsIdx);
698:       for (fs = 0; fs < num_fs; ++fs) {
699:         DMLabelGetStratumIS(fsLabel, fsIdx[fs], &stratumIS);
700:         ISGetSize(stratumIS, &fsSize);
701:         elem_list_size += fsSize;
702:         ISDestroy(&stratumIS);
703:       }
704:       if (num_fs) {
705:         PetscMalloc3(num_fs, &elem_ind, elem_list_size, &elem_list, elem_list_size, &side_list);
706:         elem_ind[0] = 0;
707:         for (fs = 0; fs < num_fs; ++fs) {
708:           DMLabelGetStratumIS(fsLabel, fsIdx[fs], &stratumIS);
709:           ISGetIndices(stratumIS, &faces);
710:           ISGetSize(stratumIS, &fsSize);
711:           /* Set Parameters */
712:           PetscCallExternal(ex_put_set_param, exo->exoid, EX_SIDE_SET, fsIdx[fs], fsSize, 0);
713:           /* Indices */
714:           if (fs < num_fs - 1) elem_ind[fs + 1] = elem_ind[fs] + fsSize;

716:           for (i = 0; i < fsSize; ++i) {
717:             /* Element List */
718:             points = NULL;
719:             DMPlexGetTransitiveClosure(dm, faces[i], PETSC_FALSE, &numPoints, &points);
720:             elem_list[elem_ind[fs] + i] = points[2] + 1;
721:             DMPlexRestoreTransitiveClosure(dm, faces[i], PETSC_FALSE, &numPoints, &points);

723:             /* Side List */
724:             points = NULL;
725:             DMPlexGetTransitiveClosure(dm, elem_list[elem_ind[fs] + i] - 1, PETSC_TRUE, &numPoints, &points);
726:             for (j = 1; j < numPoints; ++j) {
727:               if (points[j * 2] == faces[i]) break;
728:             }
729:             /* Convert HEX sides */
730:             if (numPoints == 27) {
731:               if (j == 1) {
732:                 j = 5;
733:               } else if (j == 2) {
734:                 j = 6;
735:               } else if (j == 3) {
736:                 j = 1;
737:               } else if (j == 4) {
738:                 j = 3;
739:               } else if (j == 5) {
740:                 j = 2;
741:               } else if (j == 6) {
742:                 j = 4;
743:               }
744:             }
745:             /* Convert TET sides */
746:             if (numPoints == 15) {
747:               --j;
748:               if (j == 0) j = 4;
749:             }
750:             side_list[elem_ind[fs] + i] = j;
751:             DMPlexRestoreTransitiveClosure(dm, elem_list[elem_ind[fs] + i] - 1, PETSC_TRUE, &numPoints, &points);
752:           }
753:           ISRestoreIndices(stratumIS, &faces);
754:           ISDestroy(&stratumIS);
755:         }
756:         ISRestoreIndices(fsIS, &fsIdx);
757:         ISDestroy(&fsIS);

759:         /* Put side sets */
760:         for (fs = 0; fs < num_fs; ++fs) PetscCallExternal(ex_put_set, exo->exoid, EX_SIDE_SET, fsIdx[fs], &elem_list[elem_ind[fs]], &side_list[elem_ind[fs]]);
761:         PetscFree3(elem_ind, elem_list, side_list);
762:       }
763:     }
764:     /*
765:       close the exodus file
766:     */
767:     ex_close(exo->exoid);
768:     exo->exoid = -1;
769:   }
770:   PetscSectionDestroy(&coordSection);

772:   /*
773:     reopen the file in parallel
774:   */
775:   EXO_mode = EX_WRITE;
776:   #if defined(PETSC_USE_64BIT_INDICES)
777:   EXO_mode += EX_ALL_INT64_API;
778:   #endif
779:   CPU_word_size = sizeof(PetscReal);
780:   IO_word_size  = sizeof(PetscReal);
781:   exo->exoid    = ex_open_par(exo->filename, EXO_mode, &CPU_word_size, &IO_word_size, &EXO_version, comm, MPI_INFO_NULL);
783:   return 0;
784: }

786: /*
787:   VecView_PlexExodusII_Internal - Write a Vec corresponding in an exodus file

789:   Collective on v

791:   Input Parameters:
792: + v  - The vector to be written
793: - viewer - The PetscViewerExodusII viewer associate to an exodus file

795:   Notes:
796:   The exodus result variable index is obtained by comparing the Vec name and the
797:   names of variables declared in the exodus file. For instance for a Vec named "V"
798:   the location in the exodus file will be the first match of "V", "V_X", "V_XX", "V_1", or "V_11"
799:   amongst all variables.
800:   In the event where a nodal and zonal variable both match, the function will return an error instead of
801:   possibly corrupting the file

803:   Level: beginner

805: .seealso: `EXOGetVarIndex_Internal()`, `DMPlexView_ExodusII()`, `VecView_PlexExodusII()`
806: @*/
807: PetscErrorCode VecView_PlexExodusII_Internal(Vec v, PetscViewer viewer)
808: {
809:   DM          dm;
810:   MPI_Comm    comm;
811:   PetscMPIInt rank;

813:   int         exoid, offsetN = 0, offsetZ = 0;
814:   const char *vecname;
815:   PetscInt    step;

817:   PetscObjectGetComm((PetscObject)v, &comm);
818:   MPI_Comm_rank(comm, &rank);
819:   PetscViewerExodusIIGetId(viewer, &exoid);
820:   VecGetDM(v, &dm);
821:   PetscObjectGetName((PetscObject)v, &vecname);

823:   DMGetOutputSequenceNumber(dm, &step, NULL);
824:   EXOGetVarIndex_Internal(exoid, EX_NODAL, vecname, &offsetN);
825:   EXOGetVarIndex_Internal(exoid, EX_ELEM_BLOCK, vecname, &offsetZ);
827:   if (offsetN > 0) {
828:     VecViewPlex_ExodusII_Nodal_Internal(v, exoid, (int)step + 1, offsetN);
829:   } else if (offsetZ > 0) {
830:     VecViewPlex_ExodusII_Zonal_Internal(v, exoid, (int)step + 1, offsetZ);
831:   } else SETERRQ(comm, PETSC_ERR_FILE_UNEXPECTED, "Could not find nodal or zonal variable %s in exodus file. ", vecname);
832:   return 0;
833: }

835: /*
836:   VecLoad_PlexExodusII_Internal - Write a Vec corresponding in an exodus file

838:   Collective on v

840:   Input Parameters:
841: + v  - The vector to be written
842: - viewer - The PetscViewerExodusII viewer associate to an exodus file

844:   Notes:
845:   The exodus result variable index is obtained by comparing the Vec name and the
846:   names of variables declared in the exodus file. For instance for a Vec named "V"
847:   the location in the exodus file will be the first match of "V", "V_X", "V_XX", "V_1", or "V_11"
848:   amongst all variables.
849:   In the event where a nodal and zonal variable both match, the function will return an error instead of
850:   possibly corrupting the file

852:   Level: beginner

854: .seealso: `EXOGetVarIndex_Internal()`, `DMPlexView_ExodusII()`, `VecView_PlexExodusII()`
855: @*/
856: PetscErrorCode VecLoad_PlexExodusII_Internal(Vec v, PetscViewer viewer)
857: {
858:   DM          dm;
859:   MPI_Comm    comm;
860:   PetscMPIInt rank;

862:   int         exoid, offsetN = 0, offsetZ = 0;
863:   const char *vecname;
864:   PetscInt    step;

866:   PetscObjectGetComm((PetscObject)v, &comm);
867:   MPI_Comm_rank(comm, &rank);
868:   PetscViewerExodusIIGetId(viewer, &exoid);
869:   VecGetDM(v, &dm);
870:   PetscObjectGetName((PetscObject)v, &vecname);

872:   DMGetOutputSequenceNumber(dm, &step, NULL);
873:   EXOGetVarIndex_Internal(exoid, EX_NODAL, vecname, &offsetN);
874:   EXOGetVarIndex_Internal(exoid, EX_ELEM_BLOCK, vecname, &offsetZ);
876:   if (offsetN > 0) VecLoadPlex_ExodusII_Nodal_Internal(v, exoid, (int)step + 1, offsetN);
877:   else if (offsetZ > 0) VecLoadPlex_ExodusII_Zonal_Internal(v, exoid, (int)step + 1, offsetZ);
878:   else SETERRQ(comm, PETSC_ERR_FILE_UNEXPECTED, "Could not find nodal or zonal variable %s in exodus file. ", vecname);
879:   return 0;
880: }

882: /*
883:   VecViewPlex_ExodusII_Nodal_Internal - Write a Vec corresponding to a nodal field to an exodus file

885:   Collective on v

887:   Input Parameters:
888: + v  - The vector to be written
889: . exoid - the exodus id of a file (obtained from ex_open or ex_create for instance)
890: . step - the time step to write at (exodus steps are numbered starting from 1)
891: - offset - the location of the variable in the file

893:   Notes:
894:   The exodus result nodal variable index is obtained by comparing the Vec name and the
895:   names of nodal variables declared in the exodus file. For instance for a Vec named "V"
896:   the location in the exodus file will be the first match of "V", "V_X", "V_XX", "V_1", or "V_11"
897:   amongst all nodal variables.

899:   Level: beginner

901: .seealso: `EXOGetVarIndex_Internal()`, `DMPlexView_ExodusII_Internal()`, `VecLoadNodal_PlexEXO()`, `VecViewZonal_PlexEXO()`, `VecLoadZonal_PlexEXO()`
902: @*/
903: PetscErrorCode VecViewPlex_ExodusII_Nodal_Internal(Vec v, int exoid, int step, int offset)
904: {
905:   MPI_Comm           comm;
906:   PetscMPIInt        size;
907:   DM                 dm;
908:   Vec                vNatural, vComp;
909:   const PetscScalar *varray;
910:   PetscInt           xs, xe, bs;
911:   PetscBool          useNatural;

913:   PetscObjectGetComm((PetscObject)v, &comm);
914:   MPI_Comm_size(comm, &size);
915:   VecGetDM(v, &dm);
916:   DMGetUseNatural(dm, &useNatural);
917:   useNatural = useNatural && size > 1 ? PETSC_TRUE : PETSC_FALSE;
918:   if (useNatural) {
919:     DMPlexCreateNaturalVector(dm, &vNatural);
920:     DMPlexGlobalToNaturalBegin(dm, v, vNatural);
921:     DMPlexGlobalToNaturalEnd(dm, v, vNatural);
922:   } else {
923:     vNatural = v;
924:   }

926:   /* Write local chunk of the result in the exodus file
927:      exodus stores each component of a vector-valued field as a separate variable.
928:      We assume that they are stored sequentially */
929:   VecGetOwnershipRange(vNatural, &xs, &xe);
930:   VecGetBlockSize(vNatural, &bs);
931:   if (bs == 1) {
932:     VecGetArrayRead(vNatural, &varray);
933:     PetscCallExternal(ex_put_partial_var, exoid, step, EX_NODAL, offset, 1, xs + 1, xe - xs, varray);
934:     VecRestoreArrayRead(vNatural, &varray);
935:   } else {
936:     IS       compIS;
937:     PetscInt c;

939:     ISCreateStride(comm, (xe - xs) / bs, xs, bs, &compIS);
940:     for (c = 0; c < bs; ++c) {
941:       ISStrideSetStride(compIS, (xe - xs) / bs, xs + c, bs);
942:       VecGetSubVector(vNatural, compIS, &vComp);
943:       VecGetArrayRead(vComp, &varray);
944:       PetscCallExternal(ex_put_partial_var, exoid, step, EX_NODAL, offset + c, 1, xs / bs + 1, (xe - xs) / bs, varray);
945:       VecRestoreArrayRead(vComp, &varray);
946:       VecRestoreSubVector(vNatural, compIS, &vComp);
947:     }
948:     ISDestroy(&compIS);
949:   }
950:   if (useNatural) VecDestroy(&vNatural);
951:   return 0;
952: }

954: /*
955:   VecLoadPlex_ExodusII_Nodal_Internal - Read a Vec corresponding to a nodal field from an exodus file

957:   Collective on v

959:   Input Parameters:
960: + v  - The vector to be written
961: . exoid - the exodus id of a file (obtained from ex_open or ex_create for instance)
962: . step - the time step to read at (exodus steps are numbered starting from 1)
963: - offset - the location of the variable in the file

965:   Notes:
966:   The exodus result nodal variable index is obtained by comparing the Vec name and the
967:   names of nodal variables declared in the exodus file. For instance for a Vec named "V"
968:   the location in the exodus file will be the first match of "V", "V_X", "V_XX", "V_1", or "V_11"
969:   amongst all nodal variables.

971:   Level: beginner

973: .seealso: `EXOGetVarIndex_Internal()`, `DMPlexView_ExodusII_Internal()`, `VecViewPlex_ExodusII_Nodal_Internal()`, `VecViewZonal_PlexEXO()`, `VecLoadZonal_PlexEXO()`
974: */
975: PetscErrorCode VecLoadPlex_ExodusII_Nodal_Internal(Vec v, int exoid, int step, int offset)
976: {
977:   MPI_Comm     comm;
978:   PetscMPIInt  size;
979:   DM           dm;
980:   Vec          vNatural, vComp;
981:   PetscScalar *varray;
982:   PetscInt     xs, xe, bs;
983:   PetscBool    useNatural;

985:   PetscObjectGetComm((PetscObject)v, &comm);
986:   MPI_Comm_size(comm, &size);
987:   VecGetDM(v, &dm);
988:   DMGetUseNatural(dm, &useNatural);
989:   useNatural = useNatural && size > 1 ? PETSC_TRUE : PETSC_FALSE;
990:   if (useNatural) DMPlexCreateNaturalVector(dm, &vNatural);
991:   else vNatural = v;

993:   /* Read local chunk from the file */
994:   VecGetOwnershipRange(vNatural, &xs, &xe);
995:   VecGetBlockSize(vNatural, &bs);
996:   if (bs == 1) {
997:     VecGetArray(vNatural, &varray);
998:     PetscCallExternal(ex_get_partial_var, exoid, step, EX_NODAL, offset, 1, xs + 1, xe - xs, varray);
999:     VecRestoreArray(vNatural, &varray);
1000:   } else {
1001:     IS       compIS;
1002:     PetscInt c;

1004:     ISCreateStride(comm, (xe - xs) / bs, xs, bs, &compIS);
1005:     for (c = 0; c < bs; ++c) {
1006:       ISStrideSetStride(compIS, (xe - xs) / bs, xs + c, bs);
1007:       VecGetSubVector(vNatural, compIS, &vComp);
1008:       VecGetArray(vComp, &varray);
1009:       PetscCallExternal(ex_get_partial_var, exoid, step, EX_NODAL, offset + c, 1, xs / bs + 1, (xe - xs) / bs, varray);
1010:       VecRestoreArray(vComp, &varray);
1011:       VecRestoreSubVector(vNatural, compIS, &vComp);
1012:     }
1013:     ISDestroy(&compIS);
1014:   }
1015:   if (useNatural) {
1016:     DMPlexNaturalToGlobalBegin(dm, vNatural, v);
1017:     DMPlexNaturalToGlobalEnd(dm, vNatural, v);
1018:     VecDestroy(&vNatural);
1019:   }
1020:   return 0;
1021: }

1023: /*
1024:   VecViewPlex_ExodusII_Zonal_Internal - Write a Vec corresponding to a zonal (cell based) field to an exodus file

1026:   Collective on v

1028:   Input Parameters:
1029: + v  - The vector to be written
1030: . exoid - the exodus id of a file (obtained from ex_open or ex_create for instance)
1031: . step - the time step to write at (exodus steps are numbered starting from 1)
1032: - offset - the location of the variable in the file

1034:   Notes:
1035:   The exodus result zonal variable index is obtained by comparing the Vec name and the
1036:   names of zonal variables declared in the exodus file. For instance for a Vec named "V"
1037:   the location in the exodus file will be the first match of "V", "V_X", "V_XX", "V_1", or "V_11"
1038:   amongst all zonal variables.

1040:   Level: beginner

1042: .seealso: `EXOGetVarIndex_Internal()`, `DMPlexView_ExodusII_Internal()`, `VecViewPlex_ExodusII_Nodal_Internal()`, `VecLoadPlex_ExodusII_Nodal_Internal()`, `VecLoadPlex_ExodusII_Zonal_Internal()`
1043: */
1044: PetscErrorCode VecViewPlex_ExodusII_Zonal_Internal(Vec v, int exoid, int step, int offset)
1045: {
1046:   MPI_Comm           comm;
1047:   PetscMPIInt        size;
1048:   DM                 dm;
1049:   Vec                vNatural, vComp;
1050:   const PetscScalar *varray;
1051:   PetscInt           xs, xe, bs;
1052:   PetscBool          useNatural;
1053:   IS                 compIS;
1054:   PetscInt          *csSize, *csID;
1055:   PetscInt           numCS, set, csxs = 0;

1057:   PetscObjectGetComm((PetscObject)v, &comm);
1058:   MPI_Comm_size(comm, &size);
1059:   VecGetDM(v, &dm);
1060:   DMGetUseNatural(dm, &useNatural);
1061:   useNatural = useNatural && size > 1 ? PETSC_TRUE : PETSC_FALSE;
1062:   if (useNatural) {
1063:     DMPlexCreateNaturalVector(dm, &vNatural);
1064:     DMPlexGlobalToNaturalBegin(dm, v, vNatural);
1065:     DMPlexGlobalToNaturalEnd(dm, v, vNatural);
1066:   } else {
1067:     vNatural = v;
1068:   }

1070:   /* Write local chunk of the result in the exodus file
1071:      exodus stores each component of a vector-valued field as a separate variable.
1072:      We assume that they are stored sequentially
1073:      Zonal variables are accessed one element block at a time, so we loop through the cell sets,
1074:      but once the vector has been reordered to natural size, we cannot use the label information
1075:      to figure out what to save where. */
1076:   numCS = ex_inquire_int(exoid, EX_INQ_ELEM_BLK);
1077:   PetscMalloc2(numCS, &csID, numCS, &csSize);
1078:   PetscCallExternal(ex_get_ids, exoid, EX_ELEM_BLOCK, csID);
1079:   for (set = 0; set < numCS; ++set) {
1080:     ex_block block;

1082:     block.id   = csID[set];
1083:     block.type = EX_ELEM_BLOCK;
1084:     PetscCallExternal(ex_get_block_param, exoid, &block);
1085:     csSize[set] = block.num_entry;
1086:   }
1087:   VecGetOwnershipRange(vNatural, &xs, &xe);
1088:   VecGetBlockSize(vNatural, &bs);
1089:   if (bs > 1) ISCreateStride(comm, (xe - xs) / bs, xs, bs, &compIS);
1090:   for (set = 0; set < numCS; set++) {
1091:     PetscInt csLocalSize, c;

1093:     /* range of indices for set setID[set]: csxs:csxs + csSize[set]-1
1094:        local slice of zonal values:         xs/bs,xm/bs-1
1095:        intersection:                        max(xs/bs,csxs),min(xm/bs-1,csxs + csSize[set]-1) */
1096:     csLocalSize = PetscMax(0, PetscMin(xe / bs, csxs + csSize[set]) - PetscMax(xs / bs, csxs));
1097:     if (bs == 1) {
1098:       VecGetArrayRead(vNatural, &varray);
1099:       PetscCallExternal(ex_put_partial_var, exoid, step, EX_ELEM_BLOCK, offset, csID[set], PetscMax(xs - csxs, 0) + 1, csLocalSize, &varray[PetscMax(0, csxs - xs)]);
1100:       VecRestoreArrayRead(vNatural, &varray);
1101:     } else {
1102:       for (c = 0; c < bs; ++c) {
1103:         ISStrideSetStride(compIS, (xe - xs) / bs, xs + c, bs);
1104:         VecGetSubVector(vNatural, compIS, &vComp);
1105:         VecGetArrayRead(vComp, &varray);
1106:         PetscCallExternal(ex_put_partial_var, exoid, step, EX_ELEM_BLOCK, offset + c, csID[set], PetscMax(xs / bs - csxs, 0) + 1, csLocalSize, &varray[PetscMax(0, csxs - xs / bs)]);
1107:         VecRestoreArrayRead(vComp, &varray);
1108:         VecRestoreSubVector(vNatural, compIS, &vComp);
1109:       }
1110:     }
1111:     csxs += csSize[set];
1112:   }
1113:   PetscFree2(csID, csSize);
1114:   if (bs > 1) ISDestroy(&compIS);
1115:   if (useNatural) VecDestroy(&vNatural);
1116:   return 0;
1117: }

1119: /*
1120:   VecLoadPlex_ExodusII_Zonal_Internal - Read a Vec corresponding to a zonal (cell based) field from an exodus file

1122:   Collective on v

1124:   Input Parameters:
1125: + v  - The vector to be written
1126: . exoid - the exodus id of a file (obtained from ex_open or ex_create for instance)
1127: . step - the time step to read at (exodus steps are numbered starting from 1)
1128: - offset - the location of the variable in the file

1130:   Notes:
1131:   The exodus result zonal variable index is obtained by comparing the Vec name and the
1132:   names of zonal variables declared in the exodus file. For instance for a Vec named "V"
1133:   the location in the exodus file will be the first match of "V", "V_X", "V_XX", "V_1", or "V_11"
1134:   amongst all zonal variables.

1136:   Level: beginner

1138: .seealso: `EXOGetVarIndex_Internal()`, `DMPlexView_ExodusII_Internal()`, `VecViewPlex_ExodusII_Nodal_Internal()`, `VecLoadPlex_ExodusII_Nodal_Internal()`, `VecLoadPlex_ExodusII_Zonal_Internal()`
1139: */
1140: PetscErrorCode VecLoadPlex_ExodusII_Zonal_Internal(Vec v, int exoid, int step, int offset)
1141: {
1142:   MPI_Comm     comm;
1143:   PetscMPIInt  size;
1144:   DM           dm;
1145:   Vec          vNatural, vComp;
1146:   PetscScalar *varray;
1147:   PetscInt     xs, xe, bs;
1148:   PetscBool    useNatural;
1149:   IS           compIS;
1150:   PetscInt    *csSize, *csID;
1151:   PetscInt     numCS, set, csxs = 0;

1153:   PetscObjectGetComm((PetscObject)v, &comm);
1154:   MPI_Comm_size(comm, &size);
1155:   VecGetDM(v, &dm);
1156:   DMGetUseNatural(dm, &useNatural);
1157:   useNatural = useNatural && size > 1 ? PETSC_TRUE : PETSC_FALSE;
1158:   if (useNatural) DMPlexCreateNaturalVector(dm, &vNatural);
1159:   else vNatural = v;

1161:   /* Read local chunk of the result in the exodus file
1162:      exodus stores each component of a vector-valued field as a separate variable.
1163:      We assume that they are stored sequentially
1164:      Zonal variables are accessed one element block at a time, so we loop through the cell sets,
1165:      but once the vector has been reordered to natural size, we cannot use the label information
1166:      to figure out what to save where. */
1167:   numCS = ex_inquire_int(exoid, EX_INQ_ELEM_BLK);
1168:   PetscMalloc2(numCS, &csID, numCS, &csSize);
1169:   PetscCallExternal(ex_get_ids, exoid, EX_ELEM_BLOCK, csID);
1170:   for (set = 0; set < numCS; ++set) {
1171:     ex_block block;

1173:     block.id   = csID[set];
1174:     block.type = EX_ELEM_BLOCK;
1175:     PetscCallExternal(ex_get_block_param, exoid, &block);
1176:     csSize[set] = block.num_entry;
1177:   }
1178:   VecGetOwnershipRange(vNatural, &xs, &xe);
1179:   VecGetBlockSize(vNatural, &bs);
1180:   if (bs > 1) ISCreateStride(comm, (xe - xs) / bs, xs, bs, &compIS);
1181:   for (set = 0; set < numCS; ++set) {
1182:     PetscInt csLocalSize, c;

1184:     /* range of indices for set setID[set]: csxs:csxs + csSize[set]-1
1185:        local slice of zonal values:         xs/bs,xm/bs-1
1186:        intersection:                        max(xs/bs,csxs),min(xm/bs-1,csxs + csSize[set]-1) */
1187:     csLocalSize = PetscMax(0, PetscMin(xe / bs, csxs + csSize[set]) - PetscMax(xs / bs, csxs));
1188:     if (bs == 1) {
1189:       VecGetArray(vNatural, &varray);
1190:       PetscCallExternal(ex_get_partial_var, exoid, step, EX_ELEM_BLOCK, offset, csID[set], PetscMax(xs - csxs, 0) + 1, csLocalSize, &varray[PetscMax(0, csxs - xs)]);
1191:       VecRestoreArray(vNatural, &varray);
1192:     } else {
1193:       for (c = 0; c < bs; ++c) {
1194:         ISStrideSetStride(compIS, (xe - xs) / bs, xs + c, bs);
1195:         VecGetSubVector(vNatural, compIS, &vComp);
1196:         VecGetArray(vComp, &varray);
1197:         PetscCallExternal(ex_get_partial_var, exoid, step, EX_ELEM_BLOCK, offset + c, csID[set], PetscMax(xs / bs - csxs, 0) + 1, csLocalSize, &varray[PetscMax(0, csxs - xs / bs)]);
1198:         VecRestoreArray(vComp, &varray);
1199:         VecRestoreSubVector(vNatural, compIS, &vComp);
1200:       }
1201:     }
1202:     csxs += csSize[set];
1203:   }
1204:   PetscFree2(csID, csSize);
1205:   if (bs > 1) ISDestroy(&compIS);
1206:   if (useNatural) {
1207:     DMPlexNaturalToGlobalBegin(dm, vNatural, v);
1208:     DMPlexNaturalToGlobalEnd(dm, vNatural, v);
1209:     VecDestroy(&vNatural);
1210:   }
1211:   return 0;
1212: }
1213: #endif

1215: /*@
1216:   PetscViewerExodusIIGetId - Get the file id of the `PETSCVIEWEREXODUSII` file

1218:   Logically Collective on viewer

1220:   Input Parameter:
1221: .  viewer - the `PetscViewer`

1223:   Output Parameter:
1224: .  exoid - The ExodusII file id

1226:   Level: intermediate

1228: .seealso: `PETSCVIEWEREXODUSII`, `PetscViewer`, `PetscViewerFileSetMode()`, `PetscViewerCreate()`, `PetscViewerSetType()`, `PetscViewerBinaryOpen()`
1229: @*/
1230: PetscErrorCode PetscViewerExodusIIGetId(PetscViewer viewer, int *exoid)
1231: {
1233:   PetscTryMethod(viewer, "PetscViewerGetId_C", (PetscViewer, int *), (viewer, exoid));
1234:   return 0;
1235: }

1237: /*@
1238:    PetscViewerExodusIISetOrder - Set the elements order in the exodusII file.

1240:    Collective

1242:    Input Parameters:
1243: +  viewer - the `PETSCVIEWEREXODUSII` viewer
1244: -  order - elements order

1246:    Output Parameter:

1248:    Level: beginner

1250: .seealso: `PETSCVIEWEREXODUSII`, `PetscViewer`, `PetscViewerExodusIIGetId()`, `PetscViewerExodusIIGetOrder()`, `PetscViewerExodusIISetOrder()`
1251: @*/
1252: PetscErrorCode PetscViewerExodusIISetOrder(PetscViewer viewer, PetscInt order)
1253: {
1255:   PetscTryMethod(viewer, "PetscViewerSetOrder_C", (PetscViewer, PetscInt), (viewer, order));
1256:   return 0;
1257: }

1259: /*@
1260:    PetscViewerExodusIIGetOrder - Get the elements order in the exodusII file.

1262:    Collective

1264:    Input Parameters:
1265: +  viewer - the `PETSCVIEWEREXODUSII` viewer
1266: -  order - elements order

1268:    Output Parameter:

1270:    Level: beginner

1272: .seealso: `PETSCVIEWEREXODUSII`, `PetscViewer`, `PetscViewerExodusIIGetId()`, `PetscViewerExodusIIGetOrder()`, `PetscViewerExodusIISetOrder()`
1273: @*/
1274: PetscErrorCode PetscViewerExodusIIGetOrder(PetscViewer viewer, PetscInt *order)
1275: {
1277:   PetscTryMethod(viewer, "PetscViewerGetOrder_C", (PetscViewer, PetscInt *), (viewer, order));
1278:   return 0;
1279: }

1281: /*@C
1282:    PetscViewerExodusIIOpen - Opens a file for ExodusII input/output.

1284:    Collective

1286:    Input Parameters:
1287: +  comm - MPI communicator
1288: .  name - name of file
1289: -  type - type of file
1290: .vb
1291:     FILE_MODE_WRITE - create new file for binary output
1292:     FILE_MODE_READ - open existing file for binary input
1293:     FILE_MODE_APPEND - open existing file for binary output
1294: .ve

1296:    Output Parameter:
1297: .  exo - `PETSCVIEWEREXODUSII` `PetscViewer` for Exodus II input/output to use with the specified file

1299:    Level: beginner

1301: .seealso: `PETSCVIEWEREXODUSII`, `PetscViewer`, `PetscViewerPushFormat()`, `PetscViewerDestroy()`,
1302:           `DMLoad()`, `PetscFileMode`, `PetscViewer`, `PetscViewerSetType()`, `PetscViewerFileSetMode()`, `PetscViewerFileSetName()`
1303: @*/
1304: PetscErrorCode PetscViewerExodusIIOpen(MPI_Comm comm, const char name[], PetscFileMode type, PetscViewer *exo)
1305: {
1306:   PetscViewerCreate(comm, exo);
1307:   PetscViewerSetType(*exo, PETSCVIEWEREXODUSII);
1308:   PetscViewerFileSetMode(*exo, type);
1309:   PetscViewerFileSetName(*exo, name);
1310:   PetscViewerSetFromOptions(*exo);
1311:   return 0;
1312: }

1314: /*@C
1315:   DMPlexCreateExodusFromFile - Create a `DMPLEX` mesh from an ExodusII file.

1317:   Collective

1319:   Input Parameters:
1320: + comm  - The MPI communicator
1321: . filename - The name of the ExodusII file
1322: - interpolate - Create faces and edges in the mesh

1324:   Output Parameter:
1325: . dm  - The `DM` object representing the mesh

1327:   Level: beginner

1329: .seealso: [](chapter_unstructured), `DM`, `PETSCVIEWEREXODUSII`, `DMPLEX`, `DMCreate()`, `DMPlexCreateExodus()`
1330: @*/
1331: PetscErrorCode DMPlexCreateExodusFromFile(MPI_Comm comm, const char filename[], PetscBool interpolate, DM *dm)
1332: {
1333:   PetscMPIInt rank;
1334: #if defined(PETSC_HAVE_EXODUSII)
1335:   int   CPU_word_size = sizeof(PetscReal), IO_word_size = 0, exoid = -1;
1336:   float version;
1337: #endif

1340:   MPI_Comm_rank(comm, &rank);
1341: #if defined(PETSC_HAVE_EXODUSII)
1342:   if (rank == 0) {
1343:     exoid = ex_open(filename, EX_READ, &CPU_word_size, &IO_word_size, &version);
1345:   }
1346:   DMPlexCreateExodus(comm, exoid, interpolate, dm);
1347:   if (rank == 0) PetscCallExternal(ex_close, exoid);
1348:   return 0;
1349: #else
1350:   SETERRQ(comm, PETSC_ERR_SUP, "This method requires ExodusII support. Reconfigure using --download-exodusii");
1351: #endif
1352: }

1354: #if defined(PETSC_HAVE_EXODUSII)
1355: static PetscErrorCode ExodusGetCellType_Internal(const char *elem_type, DMPolytopeType *ct)
1356: {
1357:   PetscBool flg;

1359:   *ct = DM_POLYTOPE_UNKNOWN;
1360:   PetscStrcmp(elem_type, "TRI", &flg);
1361:   if (flg) {
1362:     *ct = DM_POLYTOPE_TRIANGLE;
1363:     goto done;
1364:   }
1365:   PetscStrcmp(elem_type, "TRI3", &flg);
1366:   if (flg) {
1367:     *ct = DM_POLYTOPE_TRIANGLE;
1368:     goto done;
1369:   }
1370:   PetscStrcmp(elem_type, "QUAD", &flg);
1371:   if (flg) {
1372:     *ct = DM_POLYTOPE_QUADRILATERAL;
1373:     goto done;
1374:   }
1375:   PetscStrcmp(elem_type, "QUAD4", &flg);
1376:   if (flg) {
1377:     *ct = DM_POLYTOPE_QUADRILATERAL;
1378:     goto done;
1379:   }
1380:   PetscStrcmp(elem_type, "SHELL4", &flg);
1381:   if (flg) {
1382:     *ct = DM_POLYTOPE_QUADRILATERAL;
1383:     goto done;
1384:   }
1385:   PetscStrcmp(elem_type, "TETRA", &flg);
1386:   if (flg) {
1387:     *ct = DM_POLYTOPE_TETRAHEDRON;
1388:     goto done;
1389:   }
1390:   PetscStrcmp(elem_type, "TET4", &flg);
1391:   if (flg) {
1392:     *ct = DM_POLYTOPE_TETRAHEDRON;
1393:     goto done;
1394:   }
1395:   PetscStrcmp(elem_type, "WEDGE", &flg);
1396:   if (flg) {
1397:     *ct = DM_POLYTOPE_TRI_PRISM;
1398:     goto done;
1399:   }
1400:   PetscStrcmp(elem_type, "HEX", &flg);
1401:   if (flg) {
1402:     *ct = DM_POLYTOPE_HEXAHEDRON;
1403:     goto done;
1404:   }
1405:   PetscStrcmp(elem_type, "HEX8", &flg);
1406:   if (flg) {
1407:     *ct = DM_POLYTOPE_HEXAHEDRON;
1408:     goto done;
1409:   }
1410:   PetscStrcmp(elem_type, "HEXAHEDRON", &flg);
1411:   if (flg) {
1412:     *ct = DM_POLYTOPE_HEXAHEDRON;
1413:     goto done;
1414:   }
1415:   SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unrecognized element type %s", elem_type);
1416: done:
1417:   return 0;
1418: }
1419: #endif

1421: /*@
1422:   DMPlexCreateExodus - Create a `DMPLEX` mesh from an ExodusII file ID.

1424:   Collective

1426:   Input Parameters:
1427: + comm  - The MPI communicator
1428: . exoid - The ExodusII id associated with a exodus file and obtained using ex_open
1429: - interpolate - Create faces and edges in the mesh

1431:   Output Parameter:
1432: . dm  - The `DM` object representing the mesh

1434:   Level: beginner

1436: .seealso: [](chapter_unstructured), `DM`, `PETSCVIEWEREXODUSII`, `DMPLEX`, `DMPLEX`, `DMCreate()`
1437: @*/
1438: PetscErrorCode DMPlexCreateExodus(MPI_Comm comm, PetscInt exoid, PetscBool interpolate, DM *dm)
1439: {
1440: #if defined(PETSC_HAVE_EXODUSII)
1441:   PetscMPIInt  num_proc, rank;
1442:   DMLabel      cellSets = NULL, faceSets = NULL, vertSets = NULL;
1443:   PetscSection coordSection;
1444:   Vec          coordinates;
1445:   PetscScalar *coords;
1446:   PetscInt     coordSize, v;
1447:   /* Read from ex_get_init() */
1448:   char title[PETSC_MAX_PATH_LEN + 1];
1449:   int  dim = 0, dimEmbed = 0, numVertices = 0, numCells = 0;
1450:   int  num_cs = 0, num_vs = 0, num_fs = 0;
1451: #endif

1453: #if defined(PETSC_HAVE_EXODUSII)
1454:   MPI_Comm_rank(comm, &rank);
1455:   MPI_Comm_size(comm, &num_proc);
1456:   DMCreate(comm, dm);
1457:   DMSetType(*dm, DMPLEX);
1458:   /* Open EXODUS II file and read basic information on rank 0, then broadcast to all processors */
1459:   if (rank == 0) {
1460:     PetscMemzero(title, PETSC_MAX_PATH_LEN + 1);
1461:     PetscCallExternal(ex_get_init, exoid, title, &dimEmbed, &numVertices, &numCells, &num_cs, &num_vs, &num_fs);
1463:   }
1464:   MPI_Bcast(title, PETSC_MAX_PATH_LEN + 1, MPI_CHAR, 0, comm);
1465:   MPI_Bcast(&dim, 1, MPI_INT, 0, comm);
1466:   PetscObjectSetName((PetscObject)*dm, title);
1467:   DMPlexSetChart(*dm, 0, numCells + numVertices);
1468:   /*   We do not want this label automatically computed, instead we compute it here */
1469:   DMCreateLabel(*dm, "celltype");

1471:   /* Read cell sets information */
1472:   if (rank == 0) {
1473:     PetscInt *cone;
1474:     int       c, cs, ncs, c_loc, v, v_loc;
1475:     /* Read from ex_get_elem_blk_ids() */
1476:     int *cs_id, *cs_order;
1477:     /* Read from ex_get_elem_block() */
1478:     char buffer[PETSC_MAX_PATH_LEN + 1];
1479:     int  num_cell_in_set, num_vertex_per_cell, num_hybrid, num_attr;
1480:     /* Read from ex_get_elem_conn() */
1481:     int *cs_connect;

1483:     /* Get cell sets IDs */
1484:     PetscMalloc2(num_cs, &cs_id, num_cs, &cs_order);
1485:     PetscCallExternal(ex_get_ids, exoid, EX_ELEM_BLOCK, cs_id);
1486:     /* Read the cell set connectivity table and build mesh topology
1487:        EXO standard requires that cells in cell sets be numbered sequentially and be pairwise disjoint. */
1488:     /* Check for a hybrid mesh */
1489:     for (cs = 0, num_hybrid = 0; cs < num_cs; ++cs) {
1490:       DMPolytopeType ct;
1491:       char           elem_type[PETSC_MAX_PATH_LEN];

1493:       PetscArrayzero(elem_type, sizeof(elem_type));
1494:       PetscCallExternal(ex_get_elem_type, exoid, cs_id[cs], elem_type);
1495:       ExodusGetCellType_Internal(elem_type, &ct);
1496:       dim = PetscMax(dim, DMPolytopeTypeGetDim(ct));
1497:       PetscCallExternal(ex_get_block, exoid, EX_ELEM_BLOCK, cs_id[cs], buffer, &num_cell_in_set, &num_vertex_per_cell, 0, 0, &num_attr);
1498:       switch (ct) {
1499:       case DM_POLYTOPE_TRI_PRISM:
1500:         cs_order[cs] = cs;
1501:         ++num_hybrid;
1502:         break;
1503:       default:
1504:         for (c = cs; c > cs - num_hybrid; --c) cs_order[c] = cs_order[c - 1];
1505:         cs_order[cs - num_hybrid] = cs;
1506:       }
1507:     }
1508:     /* First set sizes */
1509:     for (ncs = 0, c = 0; ncs < num_cs; ++ncs) {
1510:       DMPolytopeType ct;
1511:       char           elem_type[PETSC_MAX_PATH_LEN];
1512:       const PetscInt cs = cs_order[ncs];

1514:       PetscArrayzero(elem_type, sizeof(elem_type));
1515:       PetscCallExternal(ex_get_elem_type, exoid, cs_id[cs], elem_type);
1516:       ExodusGetCellType_Internal(elem_type, &ct);
1517:       PetscCallExternal(ex_get_block, exoid, EX_ELEM_BLOCK, cs_id[cs], buffer, &num_cell_in_set, &num_vertex_per_cell, 0, 0, &num_attr);
1518:       for (c_loc = 0; c_loc < num_cell_in_set; ++c_loc, ++c) {
1519:         DMPlexSetConeSize(*dm, c, num_vertex_per_cell);
1520:         DMPlexSetCellType(*dm, c, ct);
1521:       }
1522:     }
1523:     for (v = numCells; v < numCells + numVertices; ++v) DMPlexSetCellType(*dm, v, DM_POLYTOPE_POINT);
1524:     DMSetUp(*dm);
1525:     for (ncs = 0, c = 0; ncs < num_cs; ++ncs) {
1526:       const PetscInt cs = cs_order[ncs];
1527:       PetscCallExternal(ex_get_block, exoid, EX_ELEM_BLOCK, cs_id[cs], buffer, &num_cell_in_set, &num_vertex_per_cell, 0, 0, &num_attr);
1528:       PetscMalloc2(num_vertex_per_cell * num_cell_in_set, &cs_connect, num_vertex_per_cell, &cone);
1529:       PetscCallExternal(ex_get_conn, exoid, EX_ELEM_BLOCK, cs_id[cs], cs_connect, NULL, NULL);
1530:       /* EXO uses Fortran-based indexing, DMPlex uses C-style and numbers cell first then vertices. */
1531:       for (c_loc = 0, v = 0; c_loc < num_cell_in_set; ++c_loc, ++c) {
1532:         DMPolytopeType ct;

1534:         for (v_loc = 0; v_loc < num_vertex_per_cell; ++v_loc, ++v) cone[v_loc] = cs_connect[v] + numCells - 1;
1535:         DMPlexGetCellType(*dm, c, &ct);
1536:         DMPlexInvertCell(ct, cone);
1537:         DMPlexSetCone(*dm, c, cone);
1538:         DMSetLabelValue_Fast(*dm, &cellSets, "Cell Sets", c, cs_id[cs]);
1539:       }
1540:       PetscFree2(cs_connect, cone);
1541:     }
1542:     PetscFree2(cs_id, cs_order);
1543:   }
1544:   {
1545:     PetscInt ints[] = {dim, dimEmbed};

1547:     MPI_Bcast(ints, 2, MPIU_INT, 0, comm);
1548:     DMSetDimension(*dm, ints[0]);
1549:     DMSetCoordinateDim(*dm, ints[1]);
1550:     dim      = ints[0];
1551:     dimEmbed = ints[1];
1552:   }
1553:   DMPlexSymmetrize(*dm);
1554:   DMPlexStratify(*dm);
1555:   if (interpolate) {
1556:     DM idm;

1558:     DMPlexInterpolate(*dm, &idm);
1559:     DMDestroy(dm);
1560:     *dm = idm;
1561:   }

1563:   /* Create vertex set label */
1564:   if (rank == 0 && (num_vs > 0)) {
1565:     int vs, v;
1566:     /* Read from ex_get_node_set_ids() */
1567:     int *vs_id;
1568:     /* Read from ex_get_node_set_param() */
1569:     int num_vertex_in_set;
1570:     /* Read from ex_get_node_set() */
1571:     int *vs_vertex_list;

1573:     /* Get vertex set ids */
1574:     PetscMalloc1(num_vs, &vs_id);
1575:     PetscCallExternal(ex_get_ids, exoid, EX_NODE_SET, vs_id);
1576:     for (vs = 0; vs < num_vs; ++vs) {
1577:       PetscCallExternal(ex_get_set_param, exoid, EX_NODE_SET, vs_id[vs], &num_vertex_in_set, NULL);
1578:       PetscMalloc1(num_vertex_in_set, &vs_vertex_list);
1579:       PetscCallExternal(ex_get_set, exoid, EX_NODE_SET, vs_id[vs], vs_vertex_list, NULL);
1580:       for (v = 0; v < num_vertex_in_set; ++v) DMSetLabelValue_Fast(*dm, &vertSets, "Vertex Sets", vs_vertex_list[v] + numCells - 1, vs_id[vs]);
1581:       PetscFree(vs_vertex_list);
1582:     }
1583:     PetscFree(vs_id);
1584:   }
1585:   /* Read coordinates */
1586:   DMGetCoordinateSection(*dm, &coordSection);
1587:   PetscSectionSetNumFields(coordSection, 1);
1588:   PetscSectionSetFieldComponents(coordSection, 0, dimEmbed);
1589:   PetscSectionSetChart(coordSection, numCells, numCells + numVertices);
1590:   for (v = numCells; v < numCells + numVertices; ++v) {
1591:     PetscSectionSetDof(coordSection, v, dimEmbed);
1592:     PetscSectionSetFieldDof(coordSection, v, 0, dimEmbed);
1593:   }
1594:   PetscSectionSetUp(coordSection);
1595:   PetscSectionGetStorageSize(coordSection, &coordSize);
1596:   VecCreate(PETSC_COMM_SELF, &coordinates);
1597:   PetscObjectSetName((PetscObject)coordinates, "coordinates");
1598:   VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);
1599:   VecSetBlockSize(coordinates, dimEmbed);
1600:   VecSetType(coordinates, VECSTANDARD);
1601:   VecGetArray(coordinates, &coords);
1602:   if (rank == 0) {
1603:     PetscReal *x, *y, *z;

1605:     PetscMalloc3(numVertices, &x, numVertices, &y, numVertices, &z);
1606:     PetscCallExternal(ex_get_coord, exoid, x, y, z);
1607:     if (dimEmbed > 0) {
1608:       for (v = 0; v < numVertices; ++v) coords[v * dimEmbed + 0] = x[v];
1609:     }
1610:     if (dimEmbed > 1) {
1611:       for (v = 0; v < numVertices; ++v) coords[v * dimEmbed + 1] = y[v];
1612:     }
1613:     if (dimEmbed > 2) {
1614:       for (v = 0; v < numVertices; ++v) coords[v * dimEmbed + 2] = z[v];
1615:     }
1616:     PetscFree3(x, y, z);
1617:   }
1618:   VecRestoreArray(coordinates, &coords);
1619:   DMSetCoordinatesLocal(*dm, coordinates);
1620:   VecDestroy(&coordinates);

1622:   /* Create side set label */
1623:   if (rank == 0 && interpolate && (num_fs > 0)) {
1624:     int fs, f, voff;
1625:     /* Read from ex_get_side_set_ids() */
1626:     int *fs_id;
1627:     /* Read from ex_get_side_set_param() */
1628:     int num_side_in_set;
1629:     /* Read from ex_get_side_set_node_list() */
1630:     int *fs_vertex_count_list, *fs_vertex_list;
1631:     /* Read side set labels */
1632:     char   fs_name[MAX_STR_LENGTH + 1];
1633:     size_t fs_name_len;

1635:     /* Get side set ids */
1636:     PetscMalloc1(num_fs, &fs_id);
1637:     PetscCallExternal(ex_get_ids, exoid, EX_SIDE_SET, fs_id);
1638:     for (fs = 0; fs < num_fs; ++fs) {
1639:       PetscCallExternal(ex_get_set_param, exoid, EX_SIDE_SET, fs_id[fs], &num_side_in_set, NULL);
1640:       PetscMalloc2(num_side_in_set, &fs_vertex_count_list, num_side_in_set * 4, &fs_vertex_list);
1641:       PetscCallExternal(ex_get_side_set_node_list, exoid, fs_id[fs], fs_vertex_count_list, fs_vertex_list);
1642:       /* Get the specific name associated with this side set ID. */
1643:       int fs_name_err = ex_get_name(exoid, EX_SIDE_SET, fs_id[fs], fs_name);
1644:       if (!fs_name_err) {
1645:         PetscStrlen(fs_name, &fs_name_len);
1646:         if (fs_name_len == 0) PetscStrncpy(fs_name, "Face Sets", MAX_STR_LENGTH);
1647:       }
1648:       for (f = 0, voff = 0; f < num_side_in_set; ++f) {
1649:         const PetscInt *faces    = NULL;
1650:         PetscInt        faceSize = fs_vertex_count_list[f], numFaces;
1651:         PetscInt        faceVertices[4], v;

1654:         for (v = 0; v < faceSize; ++v, ++voff) faceVertices[v] = fs_vertex_list[voff] + numCells - 1;
1655:         DMPlexGetFullJoin(*dm, faceSize, faceVertices, &numFaces, &faces);
1657:         DMSetLabelValue_Fast(*dm, &faceSets, "Face Sets", faces[0], fs_id[fs]);
1658:         /* Only add the label if one has been detected for this side set. */
1659:         if (!fs_name_err) DMSetLabelValue(*dm, fs_name, faces[0], fs_id[fs]);
1660:         DMPlexRestoreJoin(*dm, faceSize, faceVertices, &numFaces, &faces);
1661:       }
1662:       PetscFree2(fs_vertex_count_list, fs_vertex_list);
1663:     }
1664:     PetscFree(fs_id);
1665:   }

1667:   { /* Create Cell/Face/Vertex Sets labels at all processes */
1668:     enum {
1669:       n = 3
1670:     };
1671:     PetscBool flag[n];

1673:     flag[0] = cellSets ? PETSC_TRUE : PETSC_FALSE;
1674:     flag[1] = faceSets ? PETSC_TRUE : PETSC_FALSE;
1675:     flag[2] = vertSets ? PETSC_TRUE : PETSC_FALSE;
1676:     MPI_Bcast(flag, n, MPIU_BOOL, 0, comm);
1677:     if (flag[0]) DMCreateLabel(*dm, "Cell Sets");
1678:     if (flag[1]) DMCreateLabel(*dm, "Face Sets");
1679:     if (flag[2]) DMCreateLabel(*dm, "Vertex Sets");
1680:   }
1681:   return 0;
1682: #else
1683:   SETERRQ(comm, PETSC_ERR_SUP, "This method requires ExodusII support. Reconfigure using --download-exodusii");
1684: #endif
1685: }