Actual source code: plexexodusii2.c

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

      4: #include <netcdf.h>
      5: #include <exodusII.h>

      7: #include <petsc/private/viewerimpl.h>
      8: #include <petsc/private/viewerexodusiiimpl.h>
      9: /*@C
     10:   PETSC_VIEWER_EXODUSII_ - Creates an `PETSCVIEWEREXODUSII` `PetscViewer` shared by all processors in a communicator.

     12:   Collective; No Fortran Support

     14:   Input Parameter:
     15: . comm - the MPI communicator to share the `PETSCVIEWEREXODUSII` `PetscViewer`

     17:   Level: intermediate

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

     26: .seealso: `PETSCVIEWEREXODUSII`, `PetscViewer`, `PetscViewer`, `PetscViewerExodusIIOpen()`, `PetscViewerType`, `PetscViewerCreate()`, `PetscViewerDestroy()`
     27: @*/
     28: PetscViewer PETSC_VIEWER_EXODUSII_(MPI_Comm comm)
     29: {
     30:   PetscViewer viewer;

     32:   PetscFunctionBegin;
     33:   PetscCallNull(PetscViewerExodusIIOpen(comm, "mesh.exo", FILE_MODE_WRITE, &viewer));
     34:   PetscCallNull(PetscObjectRegisterDestroy((PetscObject)viewer));
     35:   PetscFunctionReturn(viewer);
     36: }

     38: static PetscErrorCode PetscViewerView_ExodusII(PetscViewer v, PetscViewer viewer)
     39: {
     40:   PetscViewer_ExodusII *exo = (PetscViewer_ExodusII *)v->data;

     42:   PetscFunctionBegin;
     43:   if (exo->filename) PetscCall(PetscViewerASCIIPrintf(viewer, "Filename:    %s\n", exo->filename));
     44:   if (exo->exoid) PetscCall(PetscViewerASCIIPrintf(viewer, "exoid:       %" PetscExodusIIInt_FMT "\n", exo->exoid));
     45:   if (exo->btype) PetscCall(PetscViewerASCIIPrintf(viewer, "IO Mode:     %d\n", exo->btype));
     46:   if (exo->order) PetscCall(PetscViewerASCIIPrintf(viewer, "Mesh order:  %" PetscInt_FMT "\n", exo->order));
     47:   PetscCall(PetscViewerASCIIPrintf(viewer, "Number of nodal variables:  %" PetscExodusIIInt_FMT "\n", exo->numNodalVariables));
     48:   for (int i = 0; i < exo->numNodalVariables; i++) PetscCall(PetscViewerASCIIPrintf(viewer, "   %d: %s\n", i, exo->nodalVariableNames[i]));
     49:   PetscCall(PetscViewerASCIIPrintf(viewer, "Number of zonal variables:  %" PetscExodusIIInt_FMT "\n", exo->numZonalVariables));
     50:   for (int i = 0; i < exo->numZonalVariables; i++) PetscCall(PetscViewerASCIIPrintf(viewer, "   %d: %s\n", i, exo->zonalVariableNames[i]));
     51:   PetscFunctionReturn(PETSC_SUCCESS);
     52: }

     54: static PetscErrorCode PetscViewerFlush_ExodusII(PetscViewer v)
     55: {
     56:   PetscViewer_ExodusII *exo = (PetscViewer_ExodusII *)v->data;

     58:   PetscFunctionBegin;
     59:   if (exo->exoid >= 0) PetscCallExternal(ex_update, exo->exoid);
     60:   PetscFunctionReturn(PETSC_SUCCESS);
     61: }

     63: static PetscErrorCode PetscViewerSetFromOptions_ExodusII(PetscViewer v, PetscOptionItems PetscOptionsObject)
     64: {
     65:   PetscFunctionBegin;
     66:   PetscOptionsHeadBegin(PetscOptionsObject, "ExodusII PetscViewer Options");
     67:   PetscOptionsHeadEnd();
     68:   PetscFunctionReturn(PETSC_SUCCESS);
     69: }

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

     75:   PetscFunctionBegin;
     76:   if (exo->exoid >= 0) PetscCallExternal(ex_close, exo->exoid);
     77:   for (PetscInt i = 0; i < exo->numZonalVariables; i++) PetscCall(PetscFree(exo->zonalVariableNames[i]));
     78:   PetscCall(PetscFree(exo->zonalVariableNames));
     79:   for (PetscInt i = 0; i < exo->numNodalVariables; i++) PetscCall(PetscFree(exo->nodalVariableNames[i]));
     80:   PetscCall(PetscFree(exo->nodalVariableNames));
     81:   PetscCall(PetscFree(exo->filename));
     82:   PetscCall(PetscFree(exo));
     83:   PetscCall(PetscObjectComposeFunction((PetscObject)viewer, "PetscViewerFileSetName_C", NULL));
     84:   PetscCall(PetscObjectComposeFunction((PetscObject)viewer, "PetscViewerFileGetName_C", NULL));
     85:   PetscCall(PetscObjectComposeFunction((PetscObject)viewer, "PetscViewerFileSetMode_C", NULL));
     86:   PetscCall(PetscObjectComposeFunction((PetscObject)viewer, "PetscViewerFileGetMode_C", NULL));
     87:   PetscCall(PetscObjectComposeFunction((PetscObject)viewer, "PetscViewerGetId_C", NULL));
     88:   PetscCall(PetscObjectComposeFunction((PetscObject)viewer, "PetscViewerGetOrder_C", NULL));
     89:   PetscCall(PetscObjectComposeFunction((PetscObject)viewer, "PetscViewerSetOrder_C", NULL));
     90:   PetscFunctionReturn(PETSC_SUCCESS);
     91: }

     93: static PetscErrorCode PetscViewerFileSetName_ExodusII(PetscViewer viewer, const char name[])
     94: {
     95:   PetscViewer_ExodusII *exo = (PetscViewer_ExodusII *)viewer->data;
     96:   PetscMPIInt           rank;
     97:   PetscExodusIIInt      CPU_word_size, IO_word_size, EXO_mode;
     98:   MPI_Info              mpi_info = MPI_INFO_NULL;
     99:   PetscExodusIIFloat    EXO_version;

    101:   PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)viewer), &rank));
    102:   CPU_word_size = sizeof(PetscReal);
    103:   IO_word_size  = sizeof(PetscReal);

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

    139: static PetscErrorCode PetscViewerFileGetName_ExodusII(PetscViewer viewer, const char *name[])
    140: {
    141:   PetscViewer_ExodusII *exo = (PetscViewer_ExodusII *)viewer->data;

    143:   PetscFunctionBegin;
    144:   *name = exo->filename;
    145:   PetscFunctionReturn(PETSC_SUCCESS);
    146: }

    148: static PetscErrorCode PetscViewerFileSetMode_ExodusII(PetscViewer viewer, PetscFileMode type)
    149: {
    150:   PetscViewer_ExodusII *exo = (PetscViewer_ExodusII *)viewer->data;

    152:   PetscFunctionBegin;
    153:   exo->btype = type;
    154:   PetscFunctionReturn(PETSC_SUCCESS);
    155: }

    157: static PetscErrorCode PetscViewerFileGetMode_ExodusII(PetscViewer viewer, PetscFileMode *type)
    158: {
    159:   PetscViewer_ExodusII *exo = (PetscViewer_ExodusII *)viewer->data;

    161:   PetscFunctionBegin;
    162:   *type = exo->btype;
    163:   PetscFunctionReturn(PETSC_SUCCESS);
    164: }

    166: static PetscErrorCode PetscViewerExodusIIGetId_ExodusII(PetscViewer viewer, PetscExodusIIInt *exoid)
    167: {
    168:   PetscViewer_ExodusII *exo = (PetscViewer_ExodusII *)viewer->data;

    170:   PetscFunctionBegin;
    171:   *exoid = exo->exoid;
    172:   PetscFunctionReturn(PETSC_SUCCESS);
    173: }

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

    179:   PetscFunctionBegin;
    180:   *order = exo->order;
    181:   PetscFunctionReturn(PETSC_SUCCESS);
    182: }

    184: static PetscErrorCode PetscViewerExodusIISetOrder_ExodusII(PetscViewer viewer, PetscInt order)
    185: {
    186:   PetscViewer_ExodusII *exo = (PetscViewer_ExodusII *)viewer->data;

    188:   PetscFunctionBegin;
    189:   exo->order = order;
    190:   PetscFunctionReturn(PETSC_SUCCESS);
    191: }

    193: /*@
    194:   PetscViewerExodusIISetZonalVariable - Sets the number of zonal variables in an ExodusII file

    196:   Collective;

    198:   Input Parameters:
    199: + viewer - a `PetscViewer` of type `PETSCVIEWEREXODUSII`
    200: - num    - the number of zonal variables in the ExodusII file

    202:   Level: intermediate

    204:   Notes:
    205:   The ExodusII API does not allow changing the number of variables in a file so this function will return an error
    206:   if called twice, called on a read-only file, or called on file for which the number of variables has already been specified

    208: .seealso: `PETSCVIEWEREXODUSII`, `PetscViewer`, `PetscViewerCreate()`, `PetscViewerDestroy()`, `PetscViewerExodusIIOpen()`, `PetscViewerSetType()`, `PetscViewerType`, `PetscViewerExodusIIGetZonalVariable()`
    209: @*/
    210: PetscErrorCode PetscViewerExodusIISetZonalVariable(PetscViewer viewer, PetscExodusIIInt num)
    211: {
    212:   PetscViewer_ExodusII *exo = (PetscViewer_ExodusII *)viewer->data;
    213:   MPI_Comm              comm;
    214:   PetscExodusIIInt      exoid = -1;

    216:   PetscFunctionBegin;
    217:   PetscCall(PetscObjectGetComm((PetscObject)viewer, &comm));
    218:   PetscCheck(exo->numZonalVariables == -1, comm, PETSC_ERR_SUP, "The number of zonal variables has already been set to %d and cannot be overwritten", exo->numZonalVariables);
    219:   PetscCheck((exo->btype != FILE_MODE_READ) && (exo->btype != FILE_MODE_UNDEFINED), comm, PETSC_ERR_FILE_WRITE, "Cannot set the number of variables because the file is not writable");

    221:   exo->numZonalVariables = num;
    222:   PetscCall(PetscMalloc1(num, &exo->zonalVariableNames));
    223:   for (int i = 0; i < num; i++) exo->zonalVariableNames[i] = NULL;
    224:   PetscCall(PetscViewerExodusIIGetId(viewer, &exoid));
    225:   PetscCallExternal(ex_put_variable_param, exoid, EX_ELEM_BLOCK, num);
    226:   PetscFunctionReturn(PETSC_SUCCESS);
    227: }

    229: /*@
    230:   PetscViewerExodusIISetNodalVariable - Sets the number of nodal variables in an ExodusII file

    232:   Collective;

    234:   Input Parameters:
    235: + viewer - a `PetscViewer` of type `PETSCVIEWEREXODUSII`
    236: - num    - the number of nodal variables in the ExodusII file

    238:   Level: intermediate

    240:   Notes:
    241:   The ExodusII API does not allow changing the number of variables in a file so this function will return an error
    242:   if called twice, called on a read-only file, or called on file for which the number of variables has already been specified

    244: .seealso: `PETSCVIEWEREXODUSII`, `PetscViewer`, `PetscViewerCreate()`, `PetscViewerDestroy()`, `PetscViewerExodusIIOpen()`, `PetscViewerSetType()`, `PetscViewerType`, `PetscViewerExodusIIGetNodalVariable()`
    245: @*/
    246: PetscErrorCode PetscViewerExodusIISetNodalVariable(PetscViewer viewer, PetscExodusIIInt num)
    247: {
    248:   PetscViewer_ExodusII *exo = (PetscViewer_ExodusII *)viewer->data;
    249:   MPI_Comm              comm;
    250:   PetscExodusIIInt      exoid = -1;

    252:   PetscFunctionBegin;
    253:   PetscCall(PetscObjectGetComm((PetscObject)viewer, &comm));
    254:   PetscCheck(exo->numNodalVariables == -1, comm, PETSC_ERR_SUP, "The number of nodal variables has already been set to %d and cannot be overwritten", exo->numNodalVariables);
    255:   PetscCheck((exo->btype != FILE_MODE_READ) && (exo->btype != FILE_MODE_UNDEFINED), comm, PETSC_ERR_FILE_WRITE, "Cannot set the number of variables because the file is not writable");

    257:   exo->numNodalVariables = num;
    258:   PetscCall(PetscMalloc1(num, &exo->nodalVariableNames));
    259:   for (int i = 0; i < num; i++) exo->nodalVariableNames[i] = NULL;
    260:   PetscCall(PetscViewerExodusIIGetId(viewer, &exoid));
    261:   PetscCallExternal(ex_put_variable_param, exoid, EX_NODAL, num);
    262:   PetscFunctionReturn(PETSC_SUCCESS);
    263: }

    265: /*@
    266:   PetscViewerExodusIIGetZonalVariable - Gets the number of zonal variables in an ExodusII file

    268:   Collective

    270:   Input Parameters:
    271: . viewer - a `PetscViewer` of type `PETSCVIEWEREXODUSII`

    273:   Output Parameter:
    274: . num - the number variables in the ExodusII file

    276:   Level: intermediate

    278:   Notes:
    279:   The number of variables in the ExodusII file is cached in the viewer

    281: .seealso: `PETSCVIEWEREXODUSII`, `PetscViewer`, `PetscViewerCreate()`, `PetscViewerDestroy()`, `PetscViewerExodusIIOpen()`, `PetscViewerSetType()`, `PetscViewerType`, `PetscViewerExodusIIsetZonalVariable()`
    282: @*/
    283: PetscErrorCode PetscViewerExodusIIGetZonalVariable(PetscViewer viewer, PetscExodusIIInt *num)
    284: {
    285:   PetscViewer_ExodusII *exo = (PetscViewer_ExodusII *)viewer->data;
    286:   MPI_Comm              comm;
    287:   PetscExodusIIInt      exoid = -1;

    289:   PetscFunctionBegin;
    290:   if (exo->numZonalVariables > -1) {
    291:     *num = exo->numZonalVariables;
    292:   } else {
    293:     PetscCall(PetscViewerExodusIIGetId(viewer, &exoid));
    294:     PetscCall(PetscObjectGetComm((PetscObject)viewer, &comm));
    295:     PetscCheck(exoid > 0, comm, PETSC_ERR_FILE_OPEN, "Exodus file is not open");
    296:     PetscCallExternal(ex_get_variable_param, exoid, EX_ELEM_BLOCK, num);
    297:     exo->numZonalVariables = *num;
    298:     PetscCall(PetscMalloc1(*num, &exo->zonalVariableNames));
    299:     for (int i = 0; i < *num; i++) exo->zonalVariableNames[i] = NULL;
    300:   }
    301:   PetscFunctionReturn(PETSC_SUCCESS);
    302: }

    304: /*@
    305:   PetscViewerExodusIIGetNodalVariable - Gets the number of nodal variables in an ExodusII file

    307:   Collective

    309:   Input Parameters:
    310: . viewer - a `PetscViewer` of type `PETSCVIEWEREXODUSII`

    312:   Output Parameter:
    313: . num - the number variables in the ExodusII file

    315:   Level: intermediate

    317:   Notes:
    318:   This function gets the number of nodal variables and saves it in the address of num.

    320: .seealso: `PETSCVIEWEREXODUSII`, `PetscViewer`, `PetscViewerCreate()`, `PetscViewerDestroy()`, `PetscViewerExodusIIOpen()`, `PetscViewerSetType()`, `PetscViewerType`, `PetscViewerExodusIISetNodalVariable()`
    321: @*/
    322: PetscErrorCode PetscViewerExodusIIGetNodalVariable(PetscViewer viewer, PetscExodusIIInt *num)
    323: {
    324:   PetscViewer_ExodusII *exo = (PetscViewer_ExodusII *)viewer->data;
    325:   MPI_Comm              comm;
    326:   PetscExodusIIInt      exoid = -1;

    328:   PetscFunctionBegin;
    329:   if (exo->numNodalVariables > -1) {
    330:     *num = exo->numNodalVariables;
    331:   } else {
    332:     PetscCall(PetscViewerExodusIIGetId(viewer, &exoid));
    333:     PetscCall(PetscObjectGetComm((PetscObject)viewer, &comm));
    334:     PetscCheck(exoid > 0, comm, PETSC_ERR_FILE_OPEN, "Exodus file is not open");
    335:     PetscCallExternal(ex_get_variable_param, exoid, EX_NODAL, num);
    336:     exo->numNodalVariables = *num;
    337:     PetscCall(PetscMalloc1(*num, &exo->nodalVariableNames));
    338:     for (int i = 0; i < *num; i++) exo->nodalVariableNames[i] = NULL;
    339:   }
    340:   PetscFunctionReturn(PETSC_SUCCESS);
    341: }

    343: /*@
    344:   PetscViewerExodusIISetZonalVariableName - Sets the name of a zonal variable.

    346:   Collective;

    348:   Input Parameters:
    349: + viewer - a `PetscViewer` of type `PETSCVIEWEREXODUSII`
    350: . idx    - the index for which you want to save the name
    351: - name   - string containing the name characters

    353:   Level: intermediate

    355: .seealso: `PETSCVIEWEREXODUSII`, `PetscViewer`, `PetscViewerCreate()`, `PetscViewerDestroy()`, `PetscViewerExodusIIOpen()`, `PetscViewerSetType()`, `PetscViewerType`, `PetscViewerExodusIIGetZonalVariableName()`
    356: @*/
    357: PetscErrorCode PetscViewerExodusIISetZonalVariableName(PetscViewer viewer, PetscExodusIIInt idx, const char name[])
    358: {
    359:   PetscViewer_ExodusII *exo = (PetscViewer_ExodusII *)viewer->data;

    361:   PetscFunctionBegin;
    362:   PetscCheck((idx >= 0) && (idx < exo->numZonalVariables), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Variable index out of range. Was PetscViewerExodusIISetZonalVariable called?");
    363:   PetscCall(PetscStrallocpy(name, (char **)&exo->zonalVariableNames[idx]));
    364:   PetscCallExternal(ex_put_variable_name, exo->exoid, EX_ELEM_BLOCK, idx + 1, name);
    365:   PetscFunctionReturn(PETSC_SUCCESS);
    366: }

    368: /*@
    369:   PetscViewerExodusIISetNodalVariableName - Sets the name of a nodal variable.

    371:   Collective;

    373:   Input Parameters:
    374: + viewer - a `PetscViewer` of type `PETSCVIEWEREXODUSII`
    375: . idx    - the index for which you want to save the name
    376: - name   - string containing the name characters

    378:   Level: intermediate

    380: .seealso: `PETSCVIEWEREXODUSII`, `PetscViewer`, `PetscViewerCreate()`, `PetscViewerDestroy()`, `PetscViewerExodusIIOpen()`, `PetscViewerSetType()`, `PetscViewerType`, `PetscViewerExodusIIGetNodalVariableName()`
    381: @*/
    382: PetscErrorCode PetscViewerExodusIISetNodalVariableName(PetscViewer viewer, PetscExodusIIInt idx, const char name[])
    383: {
    384:   PetscViewer_ExodusII *exo = (PetscViewer_ExodusII *)viewer->data;

    386:   PetscFunctionBegin;
    387:   PetscCheck((idx >= 0) && (idx < exo->numNodalVariables), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Variable index out of range. Was PetscViewerExodusIISetNodalVariable called?");
    388:   PetscCall(PetscStrallocpy(name, (char **)&exo->nodalVariableNames[idx]));
    389:   PetscCallExternal(ex_put_variable_name, exo->exoid, EX_NODAL, idx + 1, name);
    390:   PetscFunctionReturn(PETSC_SUCCESS);
    391: }

    393: /*@
    394:   PetscViewerExodusIIGetZonalVariableName - Gets the name of a zonal variable.

    396:   Collective;

    398:   Input Parameters:
    399: + viewer - a `PetscViewer` of type `PETSCVIEWEREXODUSII`
    400: - idx    - the index for which you want to get the name

    402:   Output Parameter:
    403: . name - pointer to the string containing the name characters

    405:   Level: intermediate

    407: .seealso: `PETSCVIEWEREXODUSII`, `PetscViewer`, `PetscViewerCreate()`, `PetscViewerDestroy()`, `PetscViewerExodusIIOpen()`, `PetscViewerSetType()`, `PetscViewerType`, `PetscViewerExodusIISetZonalVariableName()`
    408: @*/
    409: PetscErrorCode PetscViewerExodusIIGetZonalVariableName(PetscViewer viewer, PetscExodusIIInt idx, const char *name[])
    410: {
    411:   PetscViewer_ExodusII *exo   = (PetscViewer_ExodusII *)viewer->data;
    412:   PetscExodusIIInt      exoid = -1;
    413:   char                  tmpName[MAX_NAME_LENGTH + 1];

    415:   PetscFunctionBegin;
    416:   PetscCheck(idx >= 0 && idx < exo->numZonalVariables, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Variable index out of range. Was PetscViewerExodusIISetZonalVariable called?");
    417:   if (!exo->zonalVariableNames[idx]) {
    418:     PetscCall(PetscViewerExodusIIGetId(viewer, &exoid));
    419:     PetscCallExternal(ex_get_variable_name, exoid, EX_ELEM_BLOCK, idx + 1, tmpName);
    420:     PetscCall(PetscStrallocpy(tmpName, (char **)&exo->zonalVariableNames[idx]));
    421:   }
    422:   *name = exo->zonalVariableNames[idx];
    423:   PetscFunctionReturn(PETSC_SUCCESS);
    424: }

    426: /*@
    427:   PetscViewerExodusIIGetNodalVariableName - Gets the name of a nodal variable.

    429:   Collective;

    431:   Input Parameters:
    432: + viewer - a `PetscViewer` of type `PETSCVIEWEREXODUSII`
    433: - idx    - the index for which you want to save the name

    435:   Output Parameter:
    436: . name - string array containing name characters

    438:   Level: intermediate

    440: .seealso: `PETSCVIEWEREXODUSII`, `PetscViewer`, `PetscViewerCreate()`, `PetscViewerDestroy()`, `PetscViewerExodusIIOpen()`, `PetscViewerSetType()`, `PetscViewerType`, `PetscViewerExodusIISetNodalVariableName()`
    441: @*/
    442: PetscErrorCode PetscViewerExodusIIGetNodalVariableName(PetscViewer viewer, PetscExodusIIInt idx, const char *name[])
    443: {
    444:   PetscViewer_ExodusII *exo   = (PetscViewer_ExodusII *)viewer->data;
    445:   PetscExodusIIInt      exoid = -1;
    446:   char                  tmpName[MAX_NAME_LENGTH + 1];

    448:   PetscFunctionBegin;
    449:   PetscCheck((idx >= 0) && (idx < exo->numNodalVariables), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Variable index out of range. Was PetscViewerExodusIISetNodalVariable called?");
    450:   if (!exo->nodalVariableNames[idx]) {
    451:     PetscCall(PetscViewerExodusIIGetId(viewer, &exoid));
    452:     PetscCallExternal(ex_get_variable_name, exoid, EX_NODAL, idx + 1, tmpName);
    453:     PetscCall(PetscStrallocpy(tmpName, (char **)&exo->nodalVariableNames[idx]));
    454:   }
    455:   *name = exo->nodalVariableNames[idx];
    456:   PetscFunctionReturn(PETSC_SUCCESS);
    457: }

    459: /*@C
    460:   PetscViewerExodusIISetZonalVariableNames - Sets the names of all nodal variables

    462:   Collective; No Fortran Support

    464:   Input Parameters:
    465: + viewer - a `PetscViewer` of type `PETSCVIEWEREXODUSII`
    466: - names  - an array of string names to be set, the strings are copied into the `PetscViewer`

    468:   Level: intermediate

    470:   Notes:
    471:   This function allows users to set multiple zonal variable names at a time.

    473: .seealso: `PETSCVIEWEREXODUSII`, `PetscViewer`, `PetscViewerCreate()`, `PetscViewerDestroy()`, `PetscViewerExodusIIOpen()`, `PetscViewerSetType()`, `PetscViewerType`, `PetscViewerExodusIIGetZonalVariableNames()`
    474: @*/
    475: PetscErrorCode PetscViewerExodusIISetZonalVariableNames(PetscViewer viewer, const char *const names[])
    476: {
    477:   PetscExodusIIInt      numNames;
    478:   PetscExodusIIInt      exoid = -1;
    479:   PetscViewer_ExodusII *exo   = (PetscViewer_ExodusII *)viewer->data;

    481:   PetscFunctionBegin;
    482:   PetscCall(PetscViewerExodusIIGetZonalVariable(viewer, &numNames));
    483:   PetscCheck(numNames >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Number of zonal variables not set. Was PetscViewerExodusIISetZonalVariable called?");

    485:   PetscCall(PetscViewerExodusIIGetId(viewer, &exoid));
    486:   for (int i = 0; i < numNames; i++) {
    487:     PetscCall(PetscStrallocpy(names[i], &exo->zonalVariableNames[i]));
    488:     PetscCallExternal(ex_put_variable_name, exoid, EX_ELEM_BLOCK, i + 1, names[i]);
    489:   }
    490:   PetscFunctionReturn(PETSC_SUCCESS);
    491: }

    493: /*@C
    494:   PetscViewerExodusIISetNodalVariableNames - Sets the names of all nodal variables.

    496:   Collective; No Fortran Support

    498:   Input Parameters:
    499: + viewer - a `PetscViewer` of type `PETSCVIEWEREXODUSII`
    500: - names  - an array of string names to be set, the strings are copied into the `PetscViewer`

    502:   Level: intermediate

    504:   Notes:
    505:   This function allows users to set multiple nodal variable names at a time.

    507: .seealso: `PETSCVIEWEREXODUSII`, `PetscViewer`, `PetscViewerCreate()`, `PetscViewerDestroy()`, `PetscViewerExodusIIOpen()`, `PetscViewerSetType()`, `PetscViewerType`, `PetscViewerExodusIIGetNodalVariableNames()`
    508: @*/
    509: PetscErrorCode PetscViewerExodusIISetNodalVariableNames(PetscViewer viewer, const char *const names[])
    510: {
    511:   PetscExodusIIInt      numNames;
    512:   PetscExodusIIInt      exoid = -1;
    513:   PetscViewer_ExodusII *exo   = (PetscViewer_ExodusII *)viewer->data;

    515:   PetscFunctionBegin;
    516:   PetscCall(PetscViewerExodusIIGetNodalVariable(viewer, &numNames));
    517:   PetscCheck(numNames >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Number of nodal variables not set. Was PetscViewerExodusIISetNodalVariable called?");

    519:   PetscCall(PetscViewerExodusIIGetId(viewer, &exoid));
    520:   for (int i = 0; i < numNames; i++) {
    521:     PetscCall(PetscStrallocpy(names[i], &exo->nodalVariableNames[i]));
    522:     PetscCallExternal(ex_put_variable_name, exoid, EX_NODAL, i + 1, names[i]);
    523:   }
    524:   PetscFunctionReturn(PETSC_SUCCESS);
    525: }

    527: /*@C
    528:   PetscViewerExodusIIGetZonalVariableNames - Gets the names of all zonal variables.

    530:   Collective; No Fortran Support

    532:   Input Parameters:
    533: + viewer  - a `PetscViewer` of type `PETSCVIEWEREXODUSII`
    534: - numVars - the number of zonal variable names to retrieve

    536:   Output Parameter:
    537: . varNames - returns an array of char pointers where the zonal variable names are

    539:   Level: intermediate

    541:   Notes:
    542:   This function returns a borrowed pointer which should not be freed.

    544: .seealso: `PETSCVIEWEREXODUSII`, `PetscViewer`, `PetscViewerCreate()`, `PetscViewerDestroy()`, `PetscViewerExodusIIOpen()`, `PetscViewerSetType()`, `PetscViewerType`, `PetscViewerExodusIISetZonalVariableNames()`
    545: @*/
    546: PetscErrorCode PetscViewerExodusIIGetZonalVariableNames(PetscViewer viewer, PetscExodusIIInt *numVars, const char *const *varNames[])
    547: {
    548:   PetscViewer_ExodusII *exo = (PetscViewer_ExodusII *)viewer->data;
    549:   PetscExodusIIInt      idx;
    550:   char                  tmpName[MAX_NAME_LENGTH + 1];
    551:   PetscExodusIIInt      exoid = -1;

    553:   PetscFunctionBegin;
    554:   PetscCall(PetscViewerExodusIIGetZonalVariable(viewer, numVars));
    555:   /*
    556:     Cache variable names if necessary
    557:   */
    558:   for (idx = 0; idx < *numVars; idx++) {
    559:     if (!exo->zonalVariableNames[idx]) {
    560:       PetscCall(PetscViewerExodusIIGetId(viewer, &exoid));
    561:       PetscCallExternal(ex_get_variable_name, exoid, EX_ELEM_BLOCK, idx + 1, tmpName);
    562:       PetscCall(PetscStrallocpy(tmpName, (char **)&exo->zonalVariableNames[idx]));
    563:     }
    564:   }
    565:   *varNames = (const char *const *)exo->zonalVariableNames;
    566:   PetscFunctionReturn(PETSC_SUCCESS);
    567: }

    569: /*@C
    570:   PetscViewerExodusIIGetNodalVariableNames - Gets the names of all nodal variables.

    572:   Collective; No Fortran Support

    574:   Input Parameters:
    575: + viewer  - a `PetscViewer` of type `PETSCVIEWEREXODUSII`
    576: - numVars - the number of nodal variable names to retrieve

    578:   Output Parameter:
    579: . varNames - returns an array of char pointers where the nodal variable names are

    581:   Level: intermediate

    583:   Notes:
    584:   This function returns a borrowed pointer which should not be freed.

    586: .seealso: `PETSCVIEWEREXODUSII`, `PetscViewer`, `PetscViewerCreate()`, `PetscViewerDestroy()`, `PetscViewerExodusIIOpen()`, `PetscViewerSetType()`, `PetscViewerType`, `PetscViewerExodusIISetNodalVariableNames()`
    587: @*/
    588: PetscErrorCode PetscViewerExodusIIGetNodalVariableNames(PetscViewer viewer, PetscExodusIIInt *numVars, const char *const *varNames[])
    589: {
    590:   PetscViewer_ExodusII *exo = (PetscViewer_ExodusII *)viewer->data;
    591:   PetscExodusIIInt      idx;
    592:   char                  tmpName[MAX_NAME_LENGTH + 1];
    593:   PetscExodusIIInt      exoid = -1;

    595:   PetscFunctionBegin;
    596:   PetscCall(PetscViewerExodusIIGetNodalVariable(viewer, numVars));
    597:   /*
    598:     Cache variable names if necessary
    599:   */
    600:   for (idx = 0; idx < *numVars; idx++) {
    601:     if (!exo->nodalVariableNames[idx]) {
    602:       PetscCall(PetscViewerExodusIIGetId(viewer, &exoid));
    603:       PetscCallExternal(ex_get_variable_name, exoid, EX_NODAL, idx + 1, tmpName);
    604:       PetscCall(PetscStrallocpy(tmpName, (char **)&exo->nodalVariableNames[idx]));
    605:     }
    606:   }
    607:   *varNames = (const char *const *)exo->nodalVariableNames;
    608:   PetscFunctionReturn(PETSC_SUCCESS);
    609: }

    611: /*MC
    612:    PETSCVIEWEREXODUSII - A viewer that writes to an Exodus II file

    614:   Level: beginner

    616: .seealso: `PetscViewerExodusIIOpen()`, `PetscViewerCreate()`, `PETSCVIEWERBINARY`, `PETSCVIEWERHDF5`, `DMView()`,
    617:           `PetscViewerFileSetName()`, `PetscViewerFileSetMode()`, `PetscViewerFormat`, `PetscViewerType`, `PetscViewerSetType()`
    618: M*/
    619: PETSC_EXTERN PetscErrorCode PetscViewerCreate_ExodusII(PetscViewer v)
    620: {
    621:   PetscViewer_ExodusII *exo;

    623:   PetscFunctionBegin;
    624:   PetscCall(PetscNew(&exo));

    626:   v->data                 = (void *)exo;
    627:   v->ops->destroy         = PetscViewerDestroy_ExodusII;
    628:   v->ops->setfromoptions  = PetscViewerSetFromOptions_ExodusII;
    629:   v->ops->view            = PetscViewerView_ExodusII;
    630:   v->ops->flush           = PetscViewerFlush_ExodusII;
    631:   exo->btype              = FILE_MODE_UNDEFINED;
    632:   exo->filename           = 0;
    633:   exo->exoid              = -1;
    634:   exo->numNodalVariables  = -1;
    635:   exo->numZonalVariables  = -1;
    636:   exo->nodalVariableNames = NULL;
    637:   exo->zonalVariableNames = NULL;

    639:   PetscCall(PetscObjectComposeFunction((PetscObject)v, "PetscViewerFileSetName_C", PetscViewerFileSetName_ExodusII));
    640:   PetscCall(PetscObjectComposeFunction((PetscObject)v, "PetscViewerFileGetName_C", PetscViewerFileGetName_ExodusII));
    641:   PetscCall(PetscObjectComposeFunction((PetscObject)v, "PetscViewerFileSetMode_C", PetscViewerFileSetMode_ExodusII));
    642:   PetscCall(PetscObjectComposeFunction((PetscObject)v, "PetscViewerFileGetMode_C", PetscViewerFileGetMode_ExodusII));
    643:   PetscCall(PetscObjectComposeFunction((PetscObject)v, "PetscViewerGetId_C", PetscViewerExodusIIGetId_ExodusII));
    644:   PetscCall(PetscObjectComposeFunction((PetscObject)v, "PetscViewerSetOrder_C", PetscViewerExodusIISetOrder_ExodusII));
    645:   PetscCall(PetscObjectComposeFunction((PetscObject)v, "PetscViewerGetOrder_C", PetscViewerExodusIIGetOrder_ExodusII));
    646:   PetscFunctionReturn(PETSC_SUCCESS);
    647: }

    649: /*@
    650:   PetscViewerExodusIIGetNodalVariableIndex - return the location of a nodal variable in an ExodusII file given its name

    652:   Collective

    654:   Input Parameters:
    655: + viewer - a `PetscViewer` of type `PETSCVIEWEREXODUSII`
    656: - name   - the name of the result

    658:   Output Parameter:
    659: . varIndex - the location of the variable in the exodus file or -1 if the variable is not found

    661:   Level: beginner

    663:   Notes:
    664:   The exodus variable index is obtained by comparing the name argument to the
    665:   names of zonal variables declared in the exodus file. For instance if name is "V"
    666:   the location in the exodus file will be the first match of "V", "V_X", "V_XX", "V_1", or "V_11"
    667:   amongst all variables of type obj_type.

    669: .seealso: `PetscViewerExodusIISetNodalVariable()`, `PetscViewerExodusIIGetNodalVariable()`, `PetscViewerExodusIISetNodalVariableName()`, `PetscViewerExodusIISetNodalVariableNames()`, `PetscViewerExodusIIGetNodalVariableName()`, `PetscViewerExodusIIGetNodalVariableNames()`
    670: @*/
    671: PetscErrorCode PetscViewerExodusIIGetNodalVariableIndex(PetscViewer viewer, const char name[], PetscExodusIIInt *varIndex)
    672: {
    673:   PetscExodusIIInt   num_vars = 0, i, j;
    674:   char               ext_name[MAX_STR_LENGTH + 1];
    675:   const char *const *var_names;
    676:   const int          num_suffix = 5;
    677:   char              *suffix[5];
    678:   PetscBool          flg;

    680:   PetscFunctionBegin;
    681:   suffix[0] = (char *)"";
    682:   suffix[1] = (char *)"_X";
    683:   suffix[2] = (char *)"_XX";
    684:   suffix[3] = (char *)"_1";
    685:   suffix[4] = (char *)"_11";
    686:   *varIndex = -1;

    688:   PetscCall(PetscViewerExodusIIGetNodalVariableNames(viewer, &num_vars, &var_names));
    689:   for (i = 0; i < num_vars; ++i) {
    690:     for (j = 0; j < num_suffix; ++j) {
    691:       PetscCall(PetscStrncpy(ext_name, name, MAX_STR_LENGTH));
    692:       PetscCall(PetscStrlcat(ext_name, suffix[j], MAX_STR_LENGTH));
    693:       PetscCall(PetscStrcasecmp(ext_name, var_names[i], &flg));
    694:       if (flg) *varIndex = i;
    695:     }
    696:     if (flg) break;
    697:   }
    698:   PetscFunctionReturn(PETSC_SUCCESS);
    699: }

    701: /*@
    702:   PetscViewerExodusIIGetZonalVariableIndex - return the location of a zonal variable in an ExodusII file given its name

    704:   Collective

    706:   Input Parameters:
    707: + viewer - a `PetscViewer` of type `PETSCVIEWEREXODUSII`
    708: - name   - the name of the result

    710:   Output Parameter:
    711: . varIndex - the location of the variable in the exodus file or -1 if the variable is not found

    713:   Level: beginner

    715:   Notes:
    716:   The exodus variable index is obtained by comparing the name argument to the
    717:   names of zonal variables declared in the exodus file. For instance if name is "V"
    718:   the location in the exodus file will be the first match of "V", "V_X", "V_XX", "V_1", or "V_11"
    719:   amongst all variables of type obj_type.

    721: .seealso: `PetscViewerExodusIISetNodalVariable()`, `PetscViewerExodusIIGetNodalVariable()`, `PetscViewerExodusIISetNodalVariableName()`, `PetscViewerExodusIISetNodalVariableNames()`, `PetscViewerExodusIIGetNodalVariableName()`, `PetscViewerExodusIIGetNodalVariableNames()`
    722: @*/
    723: PetscErrorCode PetscViewerExodusIIGetZonalVariableIndex(PetscViewer viewer, const char name[], int *varIndex)
    724: {
    725:   PetscExodusIIInt   num_vars = 0, i, j;
    726:   char               ext_name[MAX_STR_LENGTH + 1];
    727:   const char *const *var_names;
    728:   const int          num_suffix = 5;
    729:   char              *suffix[5];
    730:   PetscBool          flg;

    732:   PetscFunctionBegin;
    733:   suffix[0] = (char *)"";
    734:   suffix[1] = (char *)"_X";
    735:   suffix[2] = (char *)"_XX";
    736:   suffix[3] = (char *)"_1";
    737:   suffix[4] = (char *)"_11";
    738:   *varIndex = -1;

    740:   PetscCall(PetscViewerExodusIIGetZonalVariableNames(viewer, &num_vars, &var_names));
    741:   for (i = 0; i < num_vars; ++i) {
    742:     for (j = 0; j < num_suffix; ++j) {
    743:       PetscCall(PetscStrncpy(ext_name, name, MAX_STR_LENGTH));
    744:       PetscCall(PetscStrlcat(ext_name, suffix[j], MAX_STR_LENGTH));
    745:       PetscCall(PetscStrcasecmp(ext_name, var_names[i], &flg));
    746:       if (flg) *varIndex = i;
    747:     }
    748:     if (flg) break;
    749:   }
    750:   PetscFunctionReturn(PETSC_SUCCESS);
    751: }

    753: PetscErrorCode DMView_PlexExodusII(DM dm, PetscViewer viewer)
    754: {
    755:   enum ElemType {
    756:     SEGMENT,
    757:     TRI,
    758:     QUAD,
    759:     TET,
    760:     HEX
    761:   };
    762:   MPI_Comm comm;
    763:   PetscInt degree; /* the order of the mesh */
    764:   /* Connectivity Variables */
    765:   PetscInt cellsNotInConnectivity;
    766:   /* Cell Sets */
    767:   DMLabel         csLabel;
    768:   IS              csIS;
    769:   const PetscInt *csIdx;
    770:   PetscInt        num_cs, cs;
    771:   enum ElemType  *type;
    772:   PetscBool       hasLabel;
    773:   /* Coordinate Variables */
    774:   DM                 cdm;
    775:   PetscSection       coordSection;
    776:   Vec                coord;
    777:   PetscInt         **nodes;
    778:   PetscInt           depth, d, dim, skipCells = 0;
    779:   PetscInt           pStart, pEnd, p, cStart, cEnd, numCells, vStart, vEnd, numVertices, eStart, eEnd, numEdges, fStart, fEnd, numFaces, numNodes;
    780:   PetscInt           num_vs, num_fs;
    781:   PetscMPIInt        rank, size;
    782:   const char        *dmName;
    783:   PetscInt           nodesLineP1[4] = {2, 0, 0, 0};
    784:   PetscInt           nodesLineP2[4] = {2, 0, 0, 1};
    785:   PetscInt           nodesTriP1[4]  = {3, 0, 0, 0};
    786:   PetscInt           nodesTriP2[4]  = {3, 3, 0, 0};
    787:   PetscInt           nodesQuadP1[4] = {4, 0, 0, 0};
    788:   PetscInt           nodesQuadP2[4] = {4, 4, 0, 1};
    789:   PetscInt           nodesTetP1[4]  = {4, 0, 0, 0};
    790:   PetscInt           nodesTetP2[4]  = {4, 6, 0, 0};
    791:   PetscInt           nodesHexP1[4]  = {8, 0, 0, 0};
    792:   PetscInt           nodesHexP2[4]  = {8, 12, 6, 1};
    793:   PetscExodusIIInt   CPU_word_size, IO_word_size, EXO_mode;
    794:   PetscExodusIIFloat EXO_version;

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

    798:   PetscFunctionBegin;
    799:   PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
    800:   PetscCallMPI(MPI_Comm_rank(comm, &rank));
    801:   PetscCallMPI(MPI_Comm_size(comm, &size));

    803:   /*
    804:     Creating coordSection is a collective operation so we do it somewhat out of sequence
    805:   */
    806:   PetscCall(PetscSectionCreate(comm, &coordSection));
    807:   PetscCall(DMGetCoordinatesLocalSetUp(dm));
    808:   /*
    809:     Check that all points are on rank 0 since we don't know how to save distributed DM in exodus format
    810:   */
    811:   PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
    812:   PetscCall(DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd));
    813:   PetscCall(DMPlexGetDepthStratum(dm, 1, &eStart, &eEnd));
    814:   PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd));
    815:   numCells    = cEnd - cStart;
    816:   numEdges    = eEnd - eStart;
    817:   numVertices = vEnd - vStart;
    818:   PetscCheck(!(rank && (numCells || numEdges || numVertices)), PETSC_COMM_SELF, PETSC_ERR_SUP, "Writing distributed DM in ExodusII format not supported");
    819:   if (rank == 0) {
    820:     switch (exo->btype) {
    821:     case FILE_MODE_READ:
    822:     case FILE_MODE_APPEND:
    823:     case FILE_MODE_UPDATE:
    824:     case FILE_MODE_APPEND_UPDATE:
    825:       /* ExodusII does not allow writing geometry to an existing file */
    826:       SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "cannot add geometry to existing file %s", exo->filename);
    827:     case FILE_MODE_WRITE:
    828:       /* Create an empty file if one already exists*/
    829:       EXO_mode = EX_CLOBBER;
    830: #if defined(PETSC_USE_64BIT_INDICES)
    831:       EXO_mode += EX_ALL_INT64_API;
    832: #endif
    833:       CPU_word_size = sizeof(PetscReal);
    834:       IO_word_size  = sizeof(PetscReal);
    835:       exo->exoid    = ex_create(exo->filename, EXO_mode, &CPU_word_size, &IO_word_size);
    836:       PetscCheck(exo->exoid >= 0, PETSC_COMM_SELF, PETSC_ERR_LIB, "ex_create failed for %s", exo->filename);

    838:       break;
    839:     default:
    840:       SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ORDER, "Must call PetscViewerFileSetMode() before PetscViewerFileSetName()");
    841:     }

    843:     /* --- Get DM info --- */
    844:     PetscCall(PetscObjectGetName((PetscObject)dm, &dmName));
    845:     PetscCall(DMPlexGetDepth(dm, &depth));
    846:     PetscCall(DMGetDimension(dm, &dim));
    847:     PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
    848:     if (depth == 3) {
    849:       numFaces = fEnd - fStart;
    850:     } else {
    851:       numFaces = 0;
    852:     }
    853:     PetscCall(DMGetLabelSize(dm, "Cell Sets", &num_cs));
    854:     PetscCall(DMGetLabelSize(dm, "Vertex Sets", &num_vs));
    855:     PetscCall(DMGetLabelSize(dm, "Face Sets", &num_fs));
    856:     PetscCall(DMGetCoordinatesLocal(dm, &coord));
    857:     PetscCall(DMGetCoordinateDM(dm, &cdm));
    858:     if (num_cs > 0) {
    859:       PetscCall(DMGetLabel(dm, "Cell Sets", &csLabel));
    860:       PetscCall(DMLabelGetValueIS(csLabel, &csIS));
    861:       PetscCall(ISGetIndices(csIS, &csIdx));
    862:     }
    863:     PetscCall(PetscMalloc1(num_cs, &nodes));
    864:     /* Set element type for each block and compute total number of nodes */
    865:     PetscCall(PetscMalloc1(num_cs, &type));
    866:     numNodes = numVertices;

    868:     PetscCall(PetscViewerExodusIIGetOrder(viewer, &degree));
    869:     if (degree == 2) numNodes += numEdges;
    870:     cellsNotInConnectivity = numCells;
    871:     for (cs = 0; cs < num_cs; ++cs) {
    872:       IS              stratumIS;
    873:       const PetscInt *cells;
    874:       PetscScalar    *xyz = NULL;
    875:       PetscInt        csSize, closureSize;

    877:       PetscCall(DMLabelGetStratumIS(csLabel, csIdx[cs], &stratumIS));
    878:       PetscCall(ISGetIndices(stratumIS, &cells));
    879:       PetscCall(ISGetSize(stratumIS, &csSize));
    880:       PetscCall(DMPlexVecGetClosure(cdm, NULL, coord, cells[0], &closureSize, &xyz));
    881:       switch (dim) {
    882:       case 1:
    883:         PetscCheck(closureSize == 2 * dim, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Number of vertices %" PetscInt_FMT " in dimension %" PetscInt_FMT " has no ExodusII type", closureSize / dim, dim);
    884:         type[cs] = SEGMENT;
    885:         break;
    886:       case 2:
    887:         if (closureSize == 3 * dim) {
    888:           type[cs] = TRI;
    889:         } else if (closureSize == 4 * dim) {
    890:           type[cs] = QUAD;
    891:         } 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);
    892:         break;
    893:       case 3:
    894:         if (closureSize == 4 * dim) {
    895:           type[cs] = TET;
    896:         } else if (closureSize == 8 * dim) {
    897:           type[cs] = HEX;
    898:         } 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);
    899:         break;
    900:       default:
    901:         SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Dimension %" PetscInt_FMT " not handled by ExodusII viewer", dim);
    902:       }
    903:       if ((degree == 2) && (type[cs] == SEGMENT)) numNodes += csSize;
    904:       if ((degree == 2) && (type[cs] == QUAD)) numNodes += csSize;
    905:       if ((degree == 2) && (type[cs] == HEX)) {
    906:         numNodes += csSize;
    907:         numNodes += numFaces;
    908:       }
    909:       PetscCall(DMPlexVecRestoreClosure(cdm, NULL, coord, cells[0], &closureSize, &xyz));
    910:       /* Set nodes and Element type */
    911:       if (type[cs] == SEGMENT) {
    912:         if (degree == 1) nodes[cs] = nodesLineP1;
    913:         else if (degree == 2) nodes[cs] = nodesLineP2;
    914:       } else if (type[cs] == TRI) {
    915:         if (degree == 1) nodes[cs] = nodesTriP1;
    916:         else if (degree == 2) nodes[cs] = nodesTriP2;
    917:       } else if (type[cs] == QUAD) {
    918:         if (degree == 1) nodes[cs] = nodesQuadP1;
    919:         else if (degree == 2) nodes[cs] = nodesQuadP2;
    920:       } else if (type[cs] == TET) {
    921:         if (degree == 1) nodes[cs] = nodesTetP1;
    922:         else if (degree == 2) nodes[cs] = nodesTetP2;
    923:       } else if (type[cs] == HEX) {
    924:         if (degree == 1) nodes[cs] = nodesHexP1;
    925:         else if (degree == 2) nodes[cs] = nodesHexP2;
    926:       }
    927:       /* Compute the number of cells not in the connectivity table */
    928:       cellsNotInConnectivity -= nodes[cs][3] * csSize;

    930:       PetscCall(ISRestoreIndices(stratumIS, &cells));
    931:       PetscCall(ISDestroy(&stratumIS));
    932:     }
    933:     if (num_cs) PetscCallExternal(ex_put_init, exo->exoid, dmName, dim, numNodes, numCells, num_cs, num_vs, num_fs);
    934:     /* --- Connectivity --- */
    935:     for (cs = 0; cs < num_cs; ++cs) {
    936:       IS              stratumIS;
    937:       const PetscInt *cells;
    938:       PetscInt       *connect, off = 0;
    939:       PetscInt        edgesInClosure = 0, facesInClosure = 0, verticesInClosure = 0;
    940:       PetscInt        csSize, c, connectSize, closureSize;
    941:       char           *elem_type        = NULL;
    942:       char            elem_type_bar2[] = "BAR2", elem_type_bar3[] = "BAR3";
    943:       char            elem_type_tri3[] = "TRI3", elem_type_quad4[] = "QUAD4";
    944:       char            elem_type_tri6[] = "TRI6", elem_type_quad9[] = "QUAD9";
    945:       char            elem_type_tet4[] = "TET4", elem_type_hex8[] = "HEX8";
    946:       char            elem_type_tet10[] = "TET10", elem_type_hex27[] = "HEX27";

    948:       PetscCall(DMLabelGetStratumIS(csLabel, csIdx[cs], &stratumIS));
    949:       PetscCall(ISGetIndices(stratumIS, &cells));
    950:       PetscCall(ISGetSize(stratumIS, &csSize));
    951:       /* Set Element type */
    952:       if (type[cs] == SEGMENT) {
    953:         if (degree == 1) elem_type = elem_type_bar2;
    954:         else if (degree == 2) elem_type = elem_type_bar3;
    955:       } else if (type[cs] == TRI) {
    956:         if (degree == 1) elem_type = elem_type_tri3;
    957:         else if (degree == 2) elem_type = elem_type_tri6;
    958:       } else if (type[cs] == QUAD) {
    959:         if (degree == 1) elem_type = elem_type_quad4;
    960:         else if (degree == 2) elem_type = elem_type_quad9;
    961:       } else if (type[cs] == TET) {
    962:         if (degree == 1) elem_type = elem_type_tet4;
    963:         else if (degree == 2) elem_type = elem_type_tet10;
    964:       } else if (type[cs] == HEX) {
    965:         if (degree == 1) elem_type = elem_type_hex8;
    966:         else if (degree == 2) elem_type = elem_type_hex27;
    967:       }
    968:       connectSize = nodes[cs][0] + nodes[cs][1] + nodes[cs][2] + nodes[cs][3];
    969:       PetscCall(PetscMalloc1(PetscMax(27, connectSize) * csSize, &connect));
    970:       PetscCallExternal(ex_put_block, exo->exoid, EX_ELEM_BLOCK, csIdx[cs], elem_type, csSize, connectSize, 0, 0, 1);
    971:       /* Find number of vertices, edges, and faces in the closure */
    972:       verticesInClosure = nodes[cs][0];
    973:       if (depth > 1) {
    974:         if (dim == 2) {
    975:           PetscCall(DMPlexGetConeSize(dm, cells[0], &edgesInClosure));
    976:         } else if (dim == 3) {
    977:           PetscInt *closure = NULL;

    979:           PetscCall(DMPlexGetConeSize(dm, cells[0], &facesInClosure));
    980:           PetscCall(DMPlexGetTransitiveClosure(dm, cells[0], PETSC_TRUE, &closureSize, &closure));
    981:           edgesInClosure = closureSize - facesInClosure - 1 - verticesInClosure;
    982:           PetscCall(DMPlexRestoreTransitiveClosure(dm, cells[0], PETSC_TRUE, &closureSize, &closure));
    983:         }
    984:       }
    985:       /* Get connectivity for each cell */
    986:       for (c = 0; c < csSize; ++c) {
    987:         PetscInt *closure = NULL;
    988:         PetscInt  temp, i;

    990:         PetscCall(DMPlexGetTransitiveClosure(dm, cells[c], PETSC_TRUE, &closureSize, &closure));
    991:         for (i = 0; i < connectSize; ++i) {
    992:           if (i < nodes[cs][0]) { /* Vertices */
    993:             connect[i + off] = closure[(i + edgesInClosure + facesInClosure + 1) * 2] + 1;
    994:             connect[i + off] -= cellsNotInConnectivity;
    995:           } else if (i < nodes[cs][0] + nodes[cs][1]) { /* Edges */
    996:             connect[i + off] = closure[(i - verticesInClosure + facesInClosure + 1) * 2] + 1;
    997:             if (nodes[cs][2] == 0) connect[i + off] -= numFaces;
    998:             connect[i + off] -= cellsNotInConnectivity;
    999:           } else if (i < nodes[cs][0] + nodes[cs][1] + nodes[cs][3]) { /* Cells */
    1000:             connect[i + off] = closure[0] + 1;
    1001:             connect[i + off] -= skipCells;
    1002:           } else if (i < nodes[cs][0] + nodes[cs][1] + nodes[cs][3] + nodes[cs][2]) { /* Faces */
    1003:             connect[i + off] = closure[(i - edgesInClosure - verticesInClosure) * 2] + 1;
    1004:             connect[i + off] -= cellsNotInConnectivity;
    1005:           } else {
    1006:             connect[i + off] = -1;
    1007:           }
    1008:         }
    1009:         /* Tetrahedra are inverted */
    1010:         if (type[cs] == TET) {
    1011:           temp             = connect[0 + off];
    1012:           connect[0 + off] = connect[1 + off];
    1013:           connect[1 + off] = temp;
    1014:           if (degree == 2) {
    1015:             temp             = connect[5 + off];
    1016:             connect[5 + off] = connect[6 + off];
    1017:             connect[6 + off] = temp;
    1018:             temp             = connect[7 + off];
    1019:             connect[7 + off] = connect[8 + off];
    1020:             connect[8 + off] = temp;
    1021:           }
    1022:         }
    1023:         /* Hexahedra are inverted */
    1024:         if (type[cs] == HEX) {
    1025:           temp             = connect[1 + off];
    1026:           connect[1 + off] = connect[3 + off];
    1027:           connect[3 + off] = temp;
    1028:           if (degree == 2) {
    1029:             temp              = connect[8 + off];
    1030:             connect[8 + off]  = connect[11 + off];
    1031:             connect[11 + off] = temp;
    1032:             temp              = connect[9 + off];
    1033:             connect[9 + off]  = connect[10 + off];
    1034:             connect[10 + off] = temp;
    1035:             temp              = connect[16 + off];
    1036:             connect[16 + off] = connect[17 + off];
    1037:             connect[17 + off] = temp;
    1038:             temp              = connect[18 + off];
    1039:             connect[18 + off] = connect[19 + off];
    1040:             connect[19 + off] = temp;

    1042:             temp              = connect[12 + off];
    1043:             connect[12 + off] = connect[16 + off];
    1044:             connect[16 + off] = temp;
    1045:             temp              = connect[13 + off];
    1046:             connect[13 + off] = connect[17 + off];
    1047:             connect[17 + off] = temp;
    1048:             temp              = connect[14 + off];
    1049:             connect[14 + off] = connect[18 + off];
    1050:             connect[18 + off] = temp;
    1051:             temp              = connect[15 + off];
    1052:             connect[15 + off] = connect[19 + off];
    1053:             connect[19 + off] = temp;

    1055:             temp              = connect[23 + off];
    1056:             connect[23 + off] = connect[26 + off];
    1057:             connect[26 + off] = temp;
    1058:             temp              = connect[24 + off];
    1059:             connect[24 + off] = connect[25 + off];
    1060:             connect[25 + off] = temp;
    1061:             temp              = connect[25 + off];
    1062:             connect[25 + off] = connect[26 + off];
    1063:             connect[26 + off] = temp;
    1064:           }
    1065:         }
    1066:         off += connectSize;
    1067:         PetscCall(DMPlexRestoreTransitiveClosure(dm, cells[c], PETSC_TRUE, &closureSize, &closure));
    1068:       }
    1069:       PetscCallExternal(ex_put_conn, exo->exoid, EX_ELEM_BLOCK, csIdx[cs], connect, 0, 0);
    1070:       skipCells += (nodes[cs][3] == 0) * csSize;
    1071:       PetscCall(PetscFree(connect));
    1072:       PetscCall(ISRestoreIndices(stratumIS, &cells));
    1073:       PetscCall(ISDestroy(&stratumIS));
    1074:     }
    1075:     PetscCall(PetscFree(type));
    1076:     /* --- Coordinates --- */
    1077:     PetscCall(PetscSectionSetChart(coordSection, pStart, pEnd));
    1078:     if (num_cs) {
    1079:       for (d = 0; d < depth; ++d) {
    1080:         PetscCall(DMPlexGetDepthStratum(dm, d, &pStart, &pEnd));
    1081:         for (p = pStart; p < pEnd; ++p) PetscCall(PetscSectionSetDof(coordSection, p, nodes[0][d] > 0));
    1082:       }
    1083:     }
    1084:     for (cs = 0; cs < num_cs; ++cs) {
    1085:       IS              stratumIS;
    1086:       const PetscInt *cells;
    1087:       PetscInt        csSize, c;

    1089:       PetscCall(DMLabelGetStratumIS(csLabel, csIdx[cs], &stratumIS));
    1090:       PetscCall(ISGetIndices(stratumIS, &cells));
    1091:       PetscCall(ISGetSize(stratumIS, &csSize));
    1092:       for (c = 0; c < csSize; ++c) PetscCall(PetscSectionSetDof(coordSection, cells[c], nodes[cs][3] > 0));
    1093:       PetscCall(ISRestoreIndices(stratumIS, &cells));
    1094:       PetscCall(ISDestroy(&stratumIS));
    1095:     }
    1096:     if (num_cs) {
    1097:       PetscCall(ISRestoreIndices(csIS, &csIdx));
    1098:       PetscCall(ISDestroy(&csIS));
    1099:     }
    1100:     PetscCall(PetscFree(nodes));
    1101:     PetscCall(PetscSectionSetUp(coordSection));
    1102:     if (numNodes) {
    1103:       const char  *coordNames[3] = {"x", "y", "z"};
    1104:       PetscScalar *closure, *cval;
    1105:       PetscReal   *coords;
    1106:       PetscInt     hasDof, n = 0;

    1108:       /* There can't be more than 24 values in the closure of a point for the coord coordSection */
    1109:       PetscCall(PetscCalloc3(numNodes * 3, &coords, dim, &cval, 24, &closure));
    1110:       PetscCall(DMGetCoordinatesLocalNoncollective(dm, &coord));
    1111:       PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
    1112:       for (p = pStart; p < pEnd; ++p) {
    1113:         PetscCall(PetscSectionGetDof(coordSection, p, &hasDof));
    1114:         if (hasDof) {
    1115:           PetscInt closureSize = 24, j;

    1117:           PetscCall(DMPlexVecGetClosure(cdm, NULL, coord, p, &closureSize, &closure));
    1118:           for (d = 0; d < dim; ++d) {
    1119:             cval[d] = 0.0;
    1120:             for (j = 0; j < closureSize / dim; j++) cval[d] += closure[j * dim + d];
    1121:             coords[d * numNodes + n] = PetscRealPart(cval[d]) * dim / closureSize;
    1122:           }
    1123:           ++n;
    1124:         }
    1125:       }
    1126:       PetscCallExternal(ex_put_coord, exo->exoid, &coords[0 * numNodes], &coords[1 * numNodes], &coords[2 * numNodes]);
    1127:       PetscCall(PetscFree3(coords, cval, closure));
    1128:       PetscCallExternal(ex_put_coord_names, exo->exoid, (char **)coordNames);
    1129:     }

    1131:     /* --- Node Sets/Vertex Sets --- */
    1132:     PetscCall(DMHasLabel(dm, "Vertex Sets", &hasLabel));
    1133:     if (hasLabel) {
    1134:       PetscInt        i, vs, vsSize;
    1135:       const PetscInt *vsIdx, *vertices;
    1136:       PetscInt       *nodeList;
    1137:       IS              vsIS, stratumIS;
    1138:       DMLabel         vsLabel;
    1139:       PetscCall(DMGetLabel(dm, "Vertex Sets", &vsLabel));
    1140:       PetscCall(DMLabelGetValueIS(vsLabel, &vsIS));
    1141:       PetscCall(ISGetIndices(vsIS, &vsIdx));
    1142:       for (vs = 0; vs < num_vs; ++vs) {
    1143:         PetscCall(DMLabelGetStratumIS(vsLabel, vsIdx[vs], &stratumIS));
    1144:         PetscCall(ISGetIndices(stratumIS, &vertices));
    1145:         PetscCall(ISGetSize(stratumIS, &vsSize));
    1146:         PetscCall(PetscMalloc1(vsSize, &nodeList));
    1147:         for (i = 0; i < vsSize; ++i) nodeList[i] = vertices[i] - skipCells + 1;
    1148:         PetscCallExternal(ex_put_set_param, exo->exoid, EX_NODE_SET, vsIdx[vs], vsSize, 0);
    1149:         PetscCallExternal(ex_put_set, exo->exoid, EX_NODE_SET, vsIdx[vs], nodeList, NULL);
    1150:         PetscCall(ISRestoreIndices(stratumIS, &vertices));
    1151:         PetscCall(ISDestroy(&stratumIS));
    1152:         PetscCall(PetscFree(nodeList));
    1153:       }
    1154:       PetscCall(ISRestoreIndices(vsIS, &vsIdx));
    1155:       PetscCall(ISDestroy(&vsIS));
    1156:     }
    1157:     /* --- Side Sets/Face Sets --- */
    1158:     PetscCall(DMHasLabel(dm, "Face Sets", &hasLabel));
    1159:     if (hasLabel) {
    1160:       PetscInt        i, j, fs, fsSize;
    1161:       const PetscInt *fsIdx, *faces;
    1162:       IS              fsIS, stratumIS;
    1163:       DMLabel         fsLabel;
    1164:       PetscInt        numPoints, *points;
    1165:       PetscInt        elem_list_size = 0;
    1166:       PetscInt       *elem_list, *elem_ind, *side_list;

    1168:       PetscCall(DMGetLabel(dm, "Face Sets", &fsLabel));
    1169:       /* Compute size of Node List and Element List */
    1170:       PetscCall(DMLabelGetValueIS(fsLabel, &fsIS));
    1171:       PetscCall(ISGetIndices(fsIS, &fsIdx));
    1172:       for (fs = 0; fs < num_fs; ++fs) {
    1173:         PetscCall(DMLabelGetStratumIS(fsLabel, fsIdx[fs], &stratumIS));
    1174:         PetscCall(ISGetSize(stratumIS, &fsSize));
    1175:         elem_list_size += fsSize;
    1176:         PetscCall(ISDestroy(&stratumIS));
    1177:       }
    1178:       if (num_fs) {
    1179:         PetscCall(PetscMalloc3(num_fs, &elem_ind, elem_list_size, &elem_list, elem_list_size, &side_list));
    1180:         elem_ind[0] = 0;
    1181:         for (fs = 0; fs < num_fs; ++fs) {
    1182:           PetscCall(DMLabelGetStratumIS(fsLabel, fsIdx[fs], &stratumIS));
    1183:           PetscCall(ISGetIndices(stratumIS, &faces));
    1184:           PetscCall(ISGetSize(stratumIS, &fsSize));
    1185:           /* Set Parameters */
    1186:           PetscCallExternal(ex_put_set_param, exo->exoid, EX_SIDE_SET, fsIdx[fs], fsSize, 0);
    1187:           /* Indices */
    1188:           if (fs < num_fs - 1) elem_ind[fs + 1] = elem_ind[fs] + fsSize;

    1190:           for (i = 0; i < fsSize; ++i) {
    1191:             /* Element List */
    1192:             points = NULL;
    1193:             PetscCall(DMPlexGetTransitiveClosure(dm, faces[i], PETSC_FALSE, &numPoints, &points));
    1194:             elem_list[elem_ind[fs] + i] = points[2] + 1;
    1195:             PetscCall(DMPlexRestoreTransitiveClosure(dm, faces[i], PETSC_FALSE, &numPoints, &points));

    1197:             /* Side List */
    1198:             points = NULL;
    1199:             PetscCall(DMPlexGetTransitiveClosure(dm, elem_list[elem_ind[fs] + i] - 1, PETSC_TRUE, &numPoints, &points));
    1200:             for (j = 1; j < numPoints; ++j) {
    1201:               if (points[j * 2] == faces[i]) break;
    1202:             }
    1203:             /* Convert HEX sides */
    1204:             if (numPoints == 27) {
    1205:               if (j == 1) {
    1206:                 j = 5;
    1207:               } else if (j == 2) {
    1208:                 j = 6;
    1209:               } else if (j == 3) {
    1210:                 j = 1;
    1211:               } else if (j == 4) {
    1212:                 j = 3;
    1213:               } else if (j == 5) {
    1214:                 j = 2;
    1215:               } else if (j == 6) {
    1216:                 j = 4;
    1217:               }
    1218:             }
    1219:             /* Convert TET sides */
    1220:             if (numPoints == 15) {
    1221:               --j;
    1222:               if (j == 0) j = 4;
    1223:             }
    1224:             side_list[elem_ind[fs] + i] = j;
    1225:             PetscCall(DMPlexRestoreTransitiveClosure(dm, elem_list[elem_ind[fs] + i] - 1, PETSC_TRUE, &numPoints, &points));
    1226:           }
    1227:           PetscCall(ISRestoreIndices(stratumIS, &faces));
    1228:           PetscCall(ISDestroy(&stratumIS));
    1229:         }
    1230:         PetscCall(ISRestoreIndices(fsIS, &fsIdx));
    1231:         PetscCall(ISDestroy(&fsIS));

    1233:         /* Put side sets */
    1234:         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]]);
    1235:         PetscCall(PetscFree3(elem_ind, elem_list, side_list));
    1236:       }
    1237:     }
    1238:     /*
    1239:       close the exodus file
    1240:     */
    1241:     ex_close(exo->exoid);
    1242:     exo->exoid = -1;
    1243:   }
    1244:   PetscCall(PetscSectionDestroy(&coordSection));

    1246:   /*
    1247:     reopen the file in parallel
    1248:   */
    1249:   EXO_mode = EX_WRITE;
    1250: #if defined(PETSC_USE_64BIT_INDICES)
    1251:   EXO_mode += EX_ALL_INT64_API;
    1252: #endif
    1253:   CPU_word_size = sizeof(PetscReal);
    1254:   IO_word_size  = sizeof(PetscReal);
    1255:   exo->exoid    = ex_open_par(exo->filename, EXO_mode, &CPU_word_size, &IO_word_size, &EXO_version, comm, MPI_INFO_NULL);
    1256:   PetscCheck(exo->exoid >= 0, PETSC_COMM_SELF, PETSC_ERR_LIB, "ex_open_par failed for %s", exo->filename);
    1257:   PetscFunctionReturn(PETSC_SUCCESS);
    1258: }

    1260: static PetscErrorCode VecViewPlex_ExodusII_Nodal_Internal(Vec v, PetscExodusIIInt exoid, PetscExodusIIInt step, PetscExodusIIInt offset);
    1261: static PetscErrorCode VecViewPlex_ExodusII_Zonal_Internal(Vec v, PetscExodusIIInt exoid, PetscExodusIIInt step, PetscExodusIIInt offset);
    1262: static PetscErrorCode VecLoadPlex_ExodusII_Nodal_Internal(Vec v, PetscExodusIIInt exoid, PetscExodusIIInt step, PetscExodusIIInt offset);
    1263: static PetscErrorCode VecLoadPlex_ExodusII_Zonal_Internal(Vec v, PetscExodusIIInt exoid, PetscExodusIIInt step, PetscExodusIIInt offset);

    1265: PetscErrorCode VecView_PlexExodusII_Internal(Vec v, PetscViewer viewer)
    1266: {
    1267:   DM          dm;
    1268:   MPI_Comm    comm;
    1269:   PetscMPIInt rank;

    1271:   PetscExodusIIInt exoid, offsetN = -1, offsetZ = -1;
    1272:   const char      *vecname;
    1273:   PetscInt         step;

    1275:   PetscFunctionBegin;
    1276:   PetscCall(PetscObjectGetComm((PetscObject)v, &comm));
    1277:   PetscCallMPI(MPI_Comm_rank(comm, &rank));
    1278:   PetscCall(PetscViewerExodusIIGetId(viewer, &exoid));
    1279:   PetscCall(VecGetDM(v, &dm));
    1280:   PetscCall(PetscObjectGetName((PetscObject)v, &vecname));

    1282:   PetscCall(DMGetOutputSequenceNumber(dm, &step, NULL));
    1283:   PetscCall(PetscViewerExodusIIGetNodalVariableIndex(viewer, vecname, &offsetN));
    1284:   PetscCall(PetscViewerExodusIIGetZonalVariableIndex(viewer, vecname, &offsetZ));
    1285:   PetscCheck(!(offsetN >= 0 && offsetZ >= 0), comm, PETSC_ERR_FILE_UNEXPECTED, "Found both nodal and zonal variable %s in exodus file. ", vecname);
    1286:   if (offsetN >= 0) {
    1287:     PetscCall(VecViewPlex_ExodusII_Nodal_Internal(v, exoid, (PetscExodusIIInt)step + 1, offsetN + 1));
    1288:   } else if (offsetZ >= 0) {
    1289:     PetscCall(VecViewPlex_ExodusII_Zonal_Internal(v, exoid, (PetscExodusIIInt)step + 1, offsetZ + 1));
    1290:   } else SETERRQ(comm, PETSC_ERR_FILE_UNEXPECTED, "Could not find nodal or zonal variable %s in exodus file. ", vecname);
    1291:   PetscFunctionReturn(PETSC_SUCCESS);
    1292: }

    1294: PetscErrorCode VecLoad_PlexExodusII_Internal(Vec v, PetscViewer viewer)
    1295: {
    1296:   DM          dm;
    1297:   MPI_Comm    comm;
    1298:   PetscMPIInt rank;

    1300:   PetscExodusIIInt exoid, offsetN = 0, offsetZ = 0;
    1301:   const char      *vecname;
    1302:   PetscInt         step;

    1304:   PetscFunctionBegin;
    1305:   PetscCall(PetscObjectGetComm((PetscObject)v, &comm));
    1306:   PetscCallMPI(MPI_Comm_rank(comm, &rank));
    1307:   PetscCall(PetscViewerExodusIIGetId(viewer, &exoid));
    1308:   PetscCall(VecGetDM(v, &dm));
    1309:   PetscCall(PetscObjectGetName((PetscObject)v, &vecname));

    1311:   PetscCall(DMGetOutputSequenceNumber(dm, &step, NULL));
    1312:   PetscCall(PetscViewerExodusIIGetNodalVariableIndex(viewer, vecname, &offsetN));
    1313:   PetscCall(PetscViewerExodusIIGetZonalVariableIndex(viewer, vecname, &offsetZ));
    1314:   PetscCheck(!(offsetN >= 0 && offsetZ >= 0), comm, PETSC_ERR_FILE_UNEXPECTED, "Found both nodal and zonal variable %s in exodus file. ", vecname);
    1315:   if (offsetN >= 0) PetscCall(VecLoadPlex_ExodusII_Nodal_Internal(v, exoid, (PetscExodusIIInt)step + 1, offsetN + 1));
    1316:   else if (offsetZ >= 0) PetscCall(VecLoadPlex_ExodusII_Zonal_Internal(v, exoid, (PetscExodusIIInt)step + 1, offsetZ + 1));
    1317:   else SETERRQ(comm, PETSC_ERR_FILE_UNEXPECTED, "Could not find nodal or zonal variable %s in exodus file. ", vecname);
    1318:   PetscFunctionReturn(PETSC_SUCCESS);
    1319: }

    1321: static PetscErrorCode VecViewPlex_ExodusII_Nodal_Internal(Vec v, PetscExodusIIInt exoid, PetscExodusIIInt step, PetscExodusIIInt offset)
    1322: {
    1323:   MPI_Comm           comm;
    1324:   PetscMPIInt        size;
    1325:   DM                 dm;
    1326:   Vec                vNatural, vComp;
    1327:   const PetscScalar *varray;
    1328:   PetscInt           xs, xe, bs;
    1329:   PetscBool          useNatural;

    1331:   PetscFunctionBegin;
    1332:   PetscCall(PetscObjectGetComm((PetscObject)v, &comm));
    1333:   PetscCallMPI(MPI_Comm_size(comm, &size));
    1334:   PetscCall(VecGetDM(v, &dm));
    1335:   PetscCall(DMGetUseNatural(dm, &useNatural));
    1336:   useNatural = useNatural && size > 1 ? PETSC_TRUE : PETSC_FALSE;
    1337:   if (useNatural) {
    1338:     PetscCall(DMPlexCreateNaturalVector(dm, &vNatural));
    1339:     PetscCall(DMPlexGlobalToNaturalBegin(dm, v, vNatural));
    1340:     PetscCall(DMPlexGlobalToNaturalEnd(dm, v, vNatural));
    1341:   } else {
    1342:     vNatural = v;
    1343:   }

    1345:   /* Write local chunk of the result in the exodus file
    1346:      exodus stores each component of a vector-valued field as a separate variable.
    1347:      We assume that they are stored sequentially */
    1348:   PetscCall(VecGetOwnershipRange(vNatural, &xs, &xe));
    1349:   PetscCall(VecGetBlockSize(vNatural, &bs));
    1350:   if (bs == 1) {
    1351:     PetscCall(VecGetArrayRead(vNatural, &varray));
    1352:     PetscCallExternal(ex_put_partial_var, exoid, step, EX_NODAL, offset, 1, xs + 1, xe - xs, varray);
    1353:     PetscCall(VecRestoreArrayRead(vNatural, &varray));
    1354:   } else {
    1355:     IS       compIS;
    1356:     PetscInt c;

    1358:     PetscCall(ISCreateStride(comm, (xe - xs) / bs, xs, bs, &compIS));
    1359:     for (c = 0; c < bs; ++c) {
    1360:       PetscCall(ISStrideSetStride(compIS, (xe - xs) / bs, xs + c, bs));
    1361:       PetscCall(VecGetSubVector(vNatural, compIS, &vComp));
    1362:       PetscCall(VecGetArrayRead(vComp, &varray));
    1363:       PetscCallExternal(ex_put_partial_var, exoid, step, EX_NODAL, offset + c, 1, xs / bs + 1, (xe - xs) / bs, varray);
    1364:       PetscCall(VecRestoreArrayRead(vComp, &varray));
    1365:       PetscCall(VecRestoreSubVector(vNatural, compIS, &vComp));
    1366:     }
    1367:     PetscCall(ISDestroy(&compIS));
    1368:   }
    1369:   if (useNatural) PetscCall(VecDestroy(&vNatural));
    1370:   PetscFunctionReturn(PETSC_SUCCESS);
    1371: }

    1373: static PetscErrorCode VecLoadPlex_ExodusII_Nodal_Internal(Vec v, PetscExodusIIInt exoid, PetscExodusIIInt step, PetscExodusIIInt offset)
    1374: {
    1375:   MPI_Comm     comm;
    1376:   PetscMPIInt  size;
    1377:   DM           dm;
    1378:   Vec          vNatural, vComp;
    1379:   PetscScalar *varray;
    1380:   PetscInt     xs, xe, bs;
    1381:   PetscBool    useNatural;

    1383:   PetscFunctionBegin;
    1384:   PetscCall(PetscObjectGetComm((PetscObject)v, &comm));
    1385:   PetscCallMPI(MPI_Comm_size(comm, &size));
    1386:   PetscCall(VecGetDM(v, &dm));
    1387:   PetscCall(DMGetUseNatural(dm, &useNatural));
    1388:   useNatural = useNatural && size > 1 ? PETSC_TRUE : PETSC_FALSE;
    1389:   if (useNatural) PetscCall(DMPlexCreateNaturalVector(dm, &vNatural));
    1390:   else vNatural = v;

    1392:   /* Read local chunk from the file */
    1393:   PetscCall(VecGetOwnershipRange(vNatural, &xs, &xe));
    1394:   PetscCall(VecGetBlockSize(vNatural, &bs));
    1395:   if (bs == 1) {
    1396:     PetscCall(VecGetArray(vNatural, &varray));
    1397:     PetscCallExternal(ex_get_partial_var, exoid, step, EX_NODAL, offset, 1, xs + 1, xe - xs, varray);
    1398:     PetscCall(VecRestoreArray(vNatural, &varray));
    1399:   } else {
    1400:     IS       compIS;
    1401:     PetscInt c;

    1403:     PetscCall(ISCreateStride(comm, (xe - xs) / bs, xs, bs, &compIS));
    1404:     for (c = 0; c < bs; ++c) {
    1405:       PetscCall(ISStrideSetStride(compIS, (xe - xs) / bs, xs + c, bs));
    1406:       PetscCall(VecGetSubVector(vNatural, compIS, &vComp));
    1407:       PetscCall(VecGetArray(vComp, &varray));
    1408:       PetscCallExternal(ex_get_partial_var, exoid, step, EX_NODAL, offset + c, 1, xs / bs + 1, (xe - xs) / bs, varray);
    1409:       PetscCall(VecRestoreArray(vComp, &varray));
    1410:       PetscCall(VecRestoreSubVector(vNatural, compIS, &vComp));
    1411:     }
    1412:     PetscCall(ISDestroy(&compIS));
    1413:   }
    1414:   if (useNatural) {
    1415:     PetscCall(DMPlexNaturalToGlobalBegin(dm, vNatural, v));
    1416:     PetscCall(DMPlexNaturalToGlobalEnd(dm, vNatural, v));
    1417:     PetscCall(VecDestroy(&vNatural));
    1418:   }
    1419:   PetscFunctionReturn(PETSC_SUCCESS);
    1420: }

    1422: static PetscErrorCode VecViewPlex_ExodusII_Zonal_Internal(Vec v, PetscExodusIIInt exoid, PetscExodusIIInt step, PetscExodusIIInt offset)
    1423: {
    1424:   MPI_Comm           comm;
    1425:   PetscMPIInt        size;
    1426:   DM                 dm;
    1427:   Vec                vNatural, vComp;
    1428:   const PetscScalar *varray;
    1429:   PetscInt           xs, xe, bs;
    1430:   PetscBool          useNatural;
    1431:   IS                 compIS;
    1432:   PetscInt          *csSize, *csID;
    1433:   PetscExodusIIInt   numCS, set, csxs = 0;

    1435:   PetscFunctionBegin;
    1436:   PetscCall(PetscObjectGetComm((PetscObject)v, &comm));
    1437:   PetscCallMPI(MPI_Comm_size(comm, &size));
    1438:   PetscCall(VecGetDM(v, &dm));
    1439:   PetscCall(DMGetUseNatural(dm, &useNatural));
    1440:   useNatural = useNatural && size > 1 ? PETSC_TRUE : PETSC_FALSE;
    1441:   if (useNatural) {
    1442:     PetscCall(DMPlexCreateNaturalVector(dm, &vNatural));
    1443:     PetscCall(DMPlexGlobalToNaturalBegin(dm, v, vNatural));
    1444:     PetscCall(DMPlexGlobalToNaturalEnd(dm, v, vNatural));
    1445:   } else {
    1446:     vNatural = v;
    1447:   }

    1449:   /* Write local chunk of the result in the exodus file
    1450:      exodus stores each component of a vector-valued field as a separate variable.
    1451:      We assume that they are stored sequentially
    1452:      Zonal variables are accessed one element block at a time, so we loop through the cell sets,
    1453:      but once the vector has been reordered to natural size, we cannot use the label information
    1454:      to figure out what to save where. */
    1455:   numCS = (PetscExodusIIInt)ex_inquire_int(exoid, EX_INQ_ELEM_BLK); // This is an int64_t
    1456:   PetscCall(PetscMalloc2(numCS, &csID, numCS, &csSize));
    1457:   PetscCallExternal(ex_get_ids, exoid, EX_ELEM_BLOCK, csID);
    1458:   for (set = 0; set < numCS; ++set) {
    1459:     ex_block block;

    1461:     block.id   = csID[set];
    1462:     block.type = EX_ELEM_BLOCK;
    1463:     PetscCallExternal(ex_get_block_param, exoid, &block);
    1464:     PetscCall(PetscIntCast(block.num_entry, &csSize[set])); // This is an int64_t
    1465:   }
    1466:   PetscCall(VecGetOwnershipRange(vNatural, &xs, &xe));
    1467:   PetscCall(VecGetBlockSize(vNatural, &bs));
    1468:   if (bs > 1) PetscCall(ISCreateStride(comm, (xe - xs) / bs, xs, bs, &compIS));
    1469:   for (set = 0; set < numCS; set++) {
    1470:     PetscInt csLocalSize, c;

    1472:     /* range of indices for set setID[set]: csxs:csxs + csSize[set]-1
    1473:        local slice of zonal values:         xs/bs,xm/bs-1
    1474:        intersection:                        max(xs/bs,csxs),min(xm/bs-1,csxs + csSize[set]-1) */
    1475:     csLocalSize = PetscMax(0, PetscMin(xe / bs, csxs + csSize[set]) - PetscMax(xs / bs, csxs));
    1476:     if (bs == 1) {
    1477:       PetscCall(VecGetArrayRead(vNatural, &varray));
    1478:       PetscCallExternal(ex_put_partial_var, exoid, step, EX_ELEM_BLOCK, offset, csID[set], PetscMax(xs - csxs, 0) + 1, csLocalSize, &varray[PetscMax(0, csxs - xs)]);
    1479:       PetscCall(VecRestoreArrayRead(vNatural, &varray));
    1480:     } else {
    1481:       for (c = 0; c < bs; ++c) {
    1482:         PetscCall(ISStrideSetStride(compIS, (xe - xs) / bs, xs + c, bs));
    1483:         PetscCall(VecGetSubVector(vNatural, compIS, &vComp));
    1484:         PetscCall(VecGetArrayRead(vComp, &varray));
    1485:         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)]);
    1486:         PetscCall(VecRestoreArrayRead(vComp, &varray));
    1487:         PetscCall(VecRestoreSubVector(vNatural, compIS, &vComp));
    1488:       }
    1489:     }
    1490:     csxs += csSize[set];
    1491:   }
    1492:   PetscCall(PetscFree2(csID, csSize));
    1493:   if (bs > 1) PetscCall(ISDestroy(&compIS));
    1494:   if (useNatural) PetscCall(VecDestroy(&vNatural));
    1495:   PetscFunctionReturn(PETSC_SUCCESS);
    1496: }

    1498: static PetscErrorCode VecLoadPlex_ExodusII_Zonal_Internal(Vec v, PetscExodusIIInt exoid, PetscExodusIIInt step, PetscExodusIIInt offset)
    1499: {
    1500:   MPI_Comm         comm;
    1501:   PetscMPIInt      size;
    1502:   DM               dm;
    1503:   Vec              vNatural, vComp;
    1504:   PetscScalar     *varray;
    1505:   PetscInt         xs, xe, bs;
    1506:   PetscBool        useNatural;
    1507:   IS               compIS;
    1508:   PetscInt        *csSize, *csID;
    1509:   PetscExodusIIInt numCS, set, csxs = 0;

    1511:   PetscFunctionBegin;
    1512:   PetscCall(PetscObjectGetComm((PetscObject)v, &comm));
    1513:   PetscCallMPI(MPI_Comm_size(comm, &size));
    1514:   PetscCall(VecGetDM(v, &dm));
    1515:   PetscCall(DMGetUseNatural(dm, &useNatural));
    1516:   useNatural = useNatural && size > 1 ? PETSC_TRUE : PETSC_FALSE;
    1517:   if (useNatural) PetscCall(DMPlexCreateNaturalVector(dm, &vNatural));
    1518:   else vNatural = v;

    1520:   /* Read local chunk of the result in the exodus file
    1521:      exodus stores each component of a vector-valued field as a separate variable.
    1522:      We assume that they are stored sequentially
    1523:      Zonal variables are accessed one element block at a time, so we loop through the cell sets,
    1524:      but once the vector has been reordered to natural size, we cannot use the label information
    1525:      to figure out what to save where. */
    1526:   numCS = (PetscExodusIIInt)ex_inquire_int(exoid, EX_INQ_ELEM_BLK); // This is an int64_t
    1527:   PetscCall(PetscMalloc2(numCS, &csID, numCS, &csSize));
    1528:   PetscCallExternal(ex_get_ids, exoid, EX_ELEM_BLOCK, csID);
    1529:   for (set = 0; set < numCS; ++set) {
    1530:     ex_block block;

    1532:     block.id   = csID[set];
    1533:     block.type = EX_ELEM_BLOCK;
    1534:     PetscCallExternal(ex_get_block_param, exoid, &block);
    1535:     PetscCall(PetscIntCast(block.num_entry, &csSize[set])); // This is an int64_t
    1536:   }
    1537:   PetscCall(VecGetOwnershipRange(vNatural, &xs, &xe));
    1538:   PetscCall(VecGetBlockSize(vNatural, &bs));
    1539:   if (bs > 1) PetscCall(ISCreateStride(comm, (xe - xs) / bs, xs, bs, &compIS));
    1540:   for (set = 0; set < numCS; ++set) {
    1541:     PetscInt csLocalSize, c;

    1543:     /* range of indices for set setID[set]: csxs:csxs + csSize[set]-1
    1544:        local slice of zonal values:         xs/bs,xm/bs-1
    1545:        intersection:                        max(xs/bs,csxs),min(xm/bs-1,csxs + csSize[set]-1) */
    1546:     csLocalSize = PetscMax(0, PetscMin(xe / bs, csxs + csSize[set]) - PetscMax(xs / bs, csxs));
    1547:     if (bs == 1) {
    1548:       PetscCall(VecGetArray(vNatural, &varray));
    1549:       PetscCallExternal(ex_get_partial_var, exoid, step, EX_ELEM_BLOCK, offset, csID[set], PetscMax(xs - csxs, 0) + 1, csLocalSize, &varray[PetscMax(0, csxs - xs)]);
    1550:       PetscCall(VecRestoreArray(vNatural, &varray));
    1551:     } else {
    1552:       for (c = 0; c < bs; ++c) {
    1553:         PetscCall(ISStrideSetStride(compIS, (xe - xs) / bs, xs + c, bs));
    1554:         PetscCall(VecGetSubVector(vNatural, compIS, &vComp));
    1555:         PetscCall(VecGetArray(vComp, &varray));
    1556:         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)]);
    1557:         PetscCall(VecRestoreArray(vComp, &varray));
    1558:         PetscCall(VecRestoreSubVector(vNatural, compIS, &vComp));
    1559:       }
    1560:     }
    1561:     csxs += csSize[set];
    1562:   }
    1563:   PetscCall(PetscFree2(csID, csSize));
    1564:   if (bs > 1) PetscCall(ISDestroy(&compIS));
    1565:   if (useNatural) {
    1566:     PetscCall(DMPlexNaturalToGlobalBegin(dm, vNatural, v));
    1567:     PetscCall(DMPlexNaturalToGlobalEnd(dm, vNatural, v));
    1568:     PetscCall(VecDestroy(&vNatural));
    1569:   }
    1570:   PetscFunctionReturn(PETSC_SUCCESS);
    1571: }

    1573: static PetscErrorCode ExodusGetCellType_Internal(const char *elem_type, DMPolytopeType *ct)
    1574: {
    1575:   PetscBool flg;

    1577:   PetscFunctionBegin;
    1578:   *ct = DM_POLYTOPE_UNKNOWN;
    1579:   PetscCall(PetscStrcmp(elem_type, "BAR2", &flg));
    1580:   if (flg) {
    1581:     *ct = DM_POLYTOPE_SEGMENT;
    1582:     goto done;
    1583:   }
    1584:   PetscCall(PetscStrcmp(elem_type, "BAR3", &flg));
    1585:   if (flg) {
    1586:     *ct = DM_POLYTOPE_SEGMENT;
    1587:     goto done;
    1588:   }
    1589:   PetscCall(PetscStrcmp(elem_type, "TRI", &flg));
    1590:   if (flg) {
    1591:     *ct = DM_POLYTOPE_TRIANGLE;
    1592:     goto done;
    1593:   }
    1594:   PetscCall(PetscStrcmp(elem_type, "TRI3", &flg));
    1595:   if (flg) {
    1596:     *ct = DM_POLYTOPE_TRIANGLE;
    1597:     goto done;
    1598:   }
    1599:   PetscCall(PetscStrcmp(elem_type, "QUAD", &flg));
    1600:   if (flg) {
    1601:     *ct = DM_POLYTOPE_QUADRILATERAL;
    1602:     goto done;
    1603:   }
    1604:   PetscCall(PetscStrcmp(elem_type, "QUAD4", &flg));
    1605:   if (flg) {
    1606:     *ct = DM_POLYTOPE_QUADRILATERAL;
    1607:     goto done;
    1608:   }
    1609:   PetscCall(PetscStrcmp(elem_type, "SHELL4", &flg));
    1610:   if (flg) {
    1611:     *ct = DM_POLYTOPE_QUADRILATERAL;
    1612:     goto done;
    1613:   }
    1614:   PetscCall(PetscStrcmp(elem_type, "TETRA", &flg));
    1615:   if (flg) {
    1616:     *ct = DM_POLYTOPE_TETRAHEDRON;
    1617:     goto done;
    1618:   }
    1619:   PetscCall(PetscStrcmp(elem_type, "TET4", &flg));
    1620:   if (flg) {
    1621:     *ct = DM_POLYTOPE_TETRAHEDRON;
    1622:     goto done;
    1623:   }
    1624:   PetscCall(PetscStrcmp(elem_type, "WEDGE", &flg));
    1625:   if (flg) {
    1626:     *ct = DM_POLYTOPE_TRI_PRISM;
    1627:     goto done;
    1628:   }
    1629:   PetscCall(PetscStrcmp(elem_type, "HEX", &flg));
    1630:   if (flg) {
    1631:     *ct = DM_POLYTOPE_HEXAHEDRON;
    1632:     goto done;
    1633:   }
    1634:   PetscCall(PetscStrcmp(elem_type, "HEX8", &flg));
    1635:   if (flg) {
    1636:     *ct = DM_POLYTOPE_HEXAHEDRON;
    1637:     goto done;
    1638:   }
    1639:   PetscCall(PetscStrcmp(elem_type, "HEXAHEDRON", &flg));
    1640:   if (flg) {
    1641:     *ct = DM_POLYTOPE_HEXAHEDRON;
    1642:     goto done;
    1643:   }
    1644:   SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unrecognized element type %s", elem_type);
    1645: done:
    1646:   PetscFunctionReturn(PETSC_SUCCESS);
    1647: }

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

    1652:   Collective

    1654:   Input Parameters:
    1655: + comm        - The MPI communicator
    1656: . exoid       - The ExodusII id associated with a exodus file and obtained using ex_open
    1657: - interpolate - Create faces and edges in the mesh

    1659:   Output Parameter:
    1660: . dm - The `DM` object representing the mesh

    1662:   Level: beginner

    1664: .seealso: [](ch_unstructured), `DM`, `PETSCVIEWEREXODUSII`, `DMPLEX`, `DMCreate()`
    1665: @*/
    1666: PetscErrorCode DMPlexCreateExodus(MPI_Comm comm, PetscExodusIIInt exoid, PetscBool interpolate, DM *dm)
    1667: {
    1668:   PetscMPIInt  num_proc, rank;
    1669:   DMLabel      cellSets = NULL, faceSets = NULL, vertSets = NULL;
    1670:   PetscSection coordSection;
    1671:   Vec          coordinates;
    1672:   PetscScalar *coords;
    1673:   PetscInt     coordSize, v;
    1674:   /* Read from ex_get_init() */
    1675:   char title[PETSC_MAX_PATH_LEN + 1];
    1676:   int  dim = 0, dimEmbed = 0, numVertices = 0, numCells = 0;
    1677:   int  num_cs = 0, num_vs = 0, num_fs = 0;

    1679:   PetscFunctionBegin;
    1680:   PetscCallMPI(MPI_Comm_rank(comm, &rank));
    1681:   PetscCallMPI(MPI_Comm_size(comm, &num_proc));
    1682:   PetscCall(DMCreate(comm, dm));
    1683:   PetscCall(DMSetType(*dm, DMPLEX));
    1684:   /* Open EXODUS II file and read basic information on rank 0, then broadcast to all processors */
    1685:   if (rank == 0) {
    1686:     PetscCall(PetscMemzero(title, PETSC_MAX_PATH_LEN + 1));
    1687:     PetscCallExternal(ex_get_init, exoid, title, &dimEmbed, &numVertices, &numCells, &num_cs, &num_vs, &num_fs);
    1688:     PetscCheck(num_cs, PETSC_COMM_SELF, PETSC_ERR_SUP, "Exodus file does not contain any cell set");
    1689:   }
    1690:   PetscCallMPI(MPI_Bcast(title, PETSC_MAX_PATH_LEN + 1, MPI_CHAR, 0, comm));
    1691:   PetscCallMPI(MPI_Bcast(&dim, 1, MPI_INT, 0, comm));
    1692:   PetscCall(PetscObjectSetName((PetscObject)*dm, title));
    1693:   PetscCall(DMPlexSetChart(*dm, 0, numCells + numVertices));
    1694:   /*   We do not want this label automatically computed, instead we compute it here */
    1695:   PetscCall(DMCreateLabel(*dm, "celltype"));

    1697:   /* Read cell sets information */
    1698:   if (rank == 0) {
    1699:     PetscInt *cone;
    1700:     int       c, cs, ncs, c_loc, v, v_loc;
    1701:     /* Read from ex_get_elem_blk_ids() */
    1702:     int *cs_id, *cs_order;
    1703:     /* Read from ex_get_elem_block() */
    1704:     char buffer[PETSC_MAX_PATH_LEN + 1];
    1705:     int  num_cell_in_set, num_vertex_per_cell, num_hybrid, num_attr;
    1706:     /* Read from ex_get_elem_conn() */
    1707:     int *cs_connect;

    1709:     /* Get cell sets IDs */
    1710:     PetscCall(PetscMalloc2(num_cs, &cs_id, num_cs, &cs_order));
    1711:     PetscCallExternal(ex_get_ids, exoid, EX_ELEM_BLOCK, cs_id);
    1712:     /* Read the cell set connectivity table and build mesh topology
    1713:        EXO standard requires that cells in cell sets be numbered sequentially and be pairwise disjoint. */
    1714:     /* Check for a hybrid mesh */
    1715:     for (cs = 0, num_hybrid = 0; cs < num_cs; ++cs) {
    1716:       DMPolytopeType ct;
    1717:       char           elem_type[PETSC_MAX_PATH_LEN];

    1719:       PetscCall(PetscArrayzero(elem_type, sizeof(elem_type)));
    1720:       PetscCallExternal(ex_get_elem_type, exoid, cs_id[cs], elem_type);
    1721:       PetscCall(ExodusGetCellType_Internal(elem_type, &ct));
    1722:       dim = PetscMax(dim, DMPolytopeTypeGetDim(ct));
    1723:       PetscCallExternal(ex_get_block, exoid, EX_ELEM_BLOCK, cs_id[cs], buffer, &num_cell_in_set, &num_vertex_per_cell, 0, 0, &num_attr);
    1724:       switch (ct) {
    1725:       case DM_POLYTOPE_TRI_PRISM:
    1726:         cs_order[cs] = cs;
    1727:         ++num_hybrid;
    1728:         break;
    1729:       default:
    1730:         for (c = cs; c > cs - num_hybrid; --c) cs_order[c] = cs_order[c - 1];
    1731:         cs_order[cs - num_hybrid] = cs;
    1732:       }
    1733:     }
    1734:     /* First set sizes */
    1735:     for (ncs = 0, c = 0; ncs < num_cs; ++ncs) {
    1736:       DMPolytopeType ct;
    1737:       char           elem_type[PETSC_MAX_PATH_LEN];
    1738:       const PetscInt cs = cs_order[ncs];

    1740:       PetscCall(PetscArrayzero(elem_type, sizeof(elem_type)));
    1741:       PetscCallExternal(ex_get_elem_type, exoid, cs_id[cs], elem_type);
    1742:       PetscCall(ExodusGetCellType_Internal(elem_type, &ct));
    1743:       PetscCallExternal(ex_get_block, exoid, EX_ELEM_BLOCK, cs_id[cs], buffer, &num_cell_in_set, &num_vertex_per_cell, 0, 0, &num_attr);
    1744:       for (c_loc = 0; c_loc < num_cell_in_set; ++c_loc, ++c) {
    1745:         PetscCall(DMPlexSetConeSize(*dm, c, num_vertex_per_cell));
    1746:         PetscCall(DMPlexSetCellType(*dm, c, ct));
    1747:       }
    1748:     }
    1749:     for (v = numCells; v < numCells + numVertices; ++v) PetscCall(DMPlexSetCellType(*dm, v, DM_POLYTOPE_POINT));
    1750:     PetscCall(DMSetUp(*dm));
    1751:     for (ncs = 0, c = 0; ncs < num_cs; ++ncs) {
    1752:       const PetscInt cs = cs_order[ncs];
    1753:       PetscCallExternal(ex_get_block, exoid, EX_ELEM_BLOCK, cs_id[cs], buffer, &num_cell_in_set, &num_vertex_per_cell, 0, 0, &num_attr);
    1754:       PetscCall(PetscMalloc2(num_vertex_per_cell * num_cell_in_set, &cs_connect, num_vertex_per_cell, &cone));
    1755:       PetscCallExternal(ex_get_conn, exoid, EX_ELEM_BLOCK, cs_id[cs], cs_connect, NULL, NULL);
    1756:       /* EXO uses Fortran-based indexing, DMPlex uses C-style and numbers cell first then vertices. */
    1757:       for (c_loc = 0, v = 0; c_loc < num_cell_in_set; ++c_loc, ++c) {
    1758:         DMPolytopeType ct;

    1760:         for (v_loc = 0; v_loc < num_vertex_per_cell; ++v_loc, ++v) cone[v_loc] = cs_connect[v] + numCells - 1;
    1761:         PetscCall(DMPlexGetCellType(*dm, c, &ct));
    1762:         PetscCall(DMPlexInvertCell(ct, cone));
    1763:         PetscCall(DMPlexSetCone(*dm, c, cone));
    1764:         PetscCall(DMSetLabelValue_Fast(*dm, &cellSets, "Cell Sets", c, cs_id[cs]));
    1765:       }
    1766:       PetscCall(PetscFree2(cs_connect, cone));
    1767:     }
    1768:     PetscCall(PetscFree2(cs_id, cs_order));
    1769:   }
    1770:   {
    1771:     PetscInt ints[] = {dim, dimEmbed};

    1773:     PetscCallMPI(MPI_Bcast(ints, 2, MPIU_INT, 0, comm));
    1774:     PetscCall(DMSetDimension(*dm, ints[0]));
    1775:     PetscCall(DMSetCoordinateDim(*dm, ints[1]));
    1776:     dim      = ints[0];
    1777:     dimEmbed = ints[1];
    1778:   }
    1779:   PetscCall(DMPlexSymmetrize(*dm));
    1780:   PetscCall(DMPlexStratify(*dm));
    1781:   if (interpolate) {
    1782:     DM idm;

    1784:     PetscCall(DMPlexInterpolate(*dm, &idm));
    1785:     PetscCall(DMDestroy(dm));
    1786:     *dm = idm;
    1787:   }

    1789:   /* Create vertex set label */
    1790:   if (rank == 0 && (num_vs > 0)) {
    1791:     int vs, v;
    1792:     /* Read from ex_get_node_set_ids() */
    1793:     int *vs_id;
    1794:     /* Read from ex_get_node_set_param() */
    1795:     int num_vertex_in_set;
    1796:     /* Read from ex_get_node_set() */
    1797:     int *vs_vertex_list;

    1799:     /* Get vertex set ids */
    1800:     PetscCall(PetscMalloc1(num_vs, &vs_id));
    1801:     PetscCallExternal(ex_get_ids, exoid, EX_NODE_SET, vs_id);
    1802:     for (vs = 0; vs < num_vs; ++vs) {
    1803:       PetscCallExternal(ex_get_set_param, exoid, EX_NODE_SET, vs_id[vs], &num_vertex_in_set, NULL);
    1804:       PetscCall(PetscMalloc1(num_vertex_in_set, &vs_vertex_list));
    1805:       PetscCallExternal(ex_get_set, exoid, EX_NODE_SET, vs_id[vs], vs_vertex_list, NULL);
    1806:       for (v = 0; v < num_vertex_in_set; ++v) PetscCall(DMSetLabelValue_Fast(*dm, &vertSets, "Vertex Sets", vs_vertex_list[v] + numCells - 1, vs_id[vs]));
    1807:       PetscCall(PetscFree(vs_vertex_list));
    1808:     }
    1809:     PetscCall(PetscFree(vs_id));
    1810:   }
    1811:   /* Read coordinates */
    1812:   PetscCall(DMGetCoordinateSection(*dm, &coordSection));
    1813:   PetscCall(PetscSectionSetNumFields(coordSection, 1));
    1814:   PetscCall(PetscSectionSetFieldComponents(coordSection, 0, dimEmbed));
    1815:   PetscCall(PetscSectionSetChart(coordSection, numCells, numCells + numVertices));
    1816:   for (v = numCells; v < numCells + numVertices; ++v) {
    1817:     PetscCall(PetscSectionSetDof(coordSection, v, dimEmbed));
    1818:     PetscCall(PetscSectionSetFieldDof(coordSection, v, 0, dimEmbed));
    1819:   }
    1820:   PetscCall(PetscSectionSetUp(coordSection));
    1821:   PetscCall(PetscSectionGetStorageSize(coordSection, &coordSize));
    1822:   PetscCall(VecCreate(PETSC_COMM_SELF, &coordinates));
    1823:   PetscCall(PetscObjectSetName((PetscObject)coordinates, "coordinates"));
    1824:   PetscCall(VecSetSizes(coordinates, coordSize, PETSC_DETERMINE));
    1825:   PetscCall(VecSetBlockSize(coordinates, dimEmbed));
    1826:   PetscCall(VecSetType(coordinates, VECSTANDARD));
    1827:   PetscCall(VecGetArray(coordinates, &coords));
    1828:   if (rank == 0) {
    1829:     PetscReal *x, *y, *z;

    1831:     PetscCall(PetscMalloc3(numVertices, &x, numVertices, &y, numVertices, &z));
    1832:     PetscCallExternal(ex_get_coord, exoid, x, y, z);
    1833:     if (dimEmbed > 0) {
    1834:       for (v = 0; v < numVertices; ++v) coords[v * dimEmbed + 0] = x[v];
    1835:     }
    1836:     if (dimEmbed > 1) {
    1837:       for (v = 0; v < numVertices; ++v) coords[v * dimEmbed + 1] = y[v];
    1838:     }
    1839:     if (dimEmbed > 2) {
    1840:       for (v = 0; v < numVertices; ++v) coords[v * dimEmbed + 2] = z[v];
    1841:     }
    1842:     PetscCall(PetscFree3(x, y, z));
    1843:   }
    1844:   PetscCall(VecRestoreArray(coordinates, &coords));
    1845:   PetscCall(DMSetCoordinatesLocal(*dm, coordinates));
    1846:   PetscCall(VecDestroy(&coordinates));

    1848:   /* Create side set label */
    1849:   if (rank == 0 && interpolate && (num_fs > 0)) {
    1850:     int fs, f, voff;
    1851:     /* Read from ex_get_side_set_ids() */
    1852:     int *fs_id;
    1853:     /* Read from ex_get_side_set_param() */
    1854:     int num_side_in_set;
    1855:     /* Read from ex_get_side_set_node_list() */
    1856:     int *fs_vertex_count_list, *fs_vertex_list, *fs_side_list;
    1857:     /* Read side set labels */
    1858:     char   fs_name[MAX_STR_LENGTH + 1];
    1859:     size_t fs_name_len;

    1861:     /* Get side set ids */
    1862:     PetscCall(PetscMalloc1(num_fs, &fs_id));
    1863:     PetscCallExternal(ex_get_ids, exoid, EX_SIDE_SET, fs_id);
    1864:     // Ids 1 and 2 are reserved by ExodusII for indicating things in 3D
    1865:     for (fs = 0; fs < num_fs; ++fs) {
    1866:       PetscCallExternal(ex_get_set_param, exoid, EX_SIDE_SET, fs_id[fs], &num_side_in_set, NULL);
    1867:       PetscCall(PetscMalloc3(num_side_in_set, &fs_vertex_count_list, num_side_in_set * 4, &fs_vertex_list, num_side_in_set, &fs_side_list));
    1868:       PetscCallExternal(ex_get_side_set_node_list, exoid, fs_id[fs], fs_vertex_count_list, fs_vertex_list);
    1869:       PetscCallExternal(ex_get_set, exoid, EX_SIDE_SET, fs_id[fs], NULL, fs_side_list);

    1871:       /* Get the specific name associated with this side set ID. */
    1872:       int fs_name_err = ex_get_name(exoid, EX_SIDE_SET, fs_id[fs], fs_name);
    1873:       if (!fs_name_err) {
    1874:         PetscCall(PetscStrlen(fs_name, &fs_name_len));
    1875:         if (fs_name_len == 0) PetscCall(PetscStrncpy(fs_name, "Face Sets", MAX_STR_LENGTH));
    1876:       }
    1877:       for (f = 0, voff = 0; f < num_side_in_set; ++f) {
    1878:         const PetscInt *faces    = NULL;
    1879:         PetscInt        faceSize = fs_vertex_count_list[f], numFaces;
    1880:         PetscInt        faceVertices[4], v;

    1882:         PetscCheck(faceSize <= 4, comm, PETSC_ERR_ARG_WRONG, "ExodusII side cannot have %" PetscInt_FMT " > 4 vertices", faceSize);
    1883:         for (v = 0; v < faceSize; ++v, ++voff) faceVertices[v] = fs_vertex_list[voff] + numCells - 1;
    1884:         PetscCall(DMPlexGetFullJoin(*dm, faceSize, faceVertices, &numFaces, &faces));
    1885:         PetscCheck(numFaces == 1, comm, PETSC_ERR_ARG_WRONG, "Invalid ExodusII side %d in set %d maps to %" PetscInt_FMT " faces", f, fs, numFaces);
    1886:         PetscCheck(dim == 1 || faces[0] >= numCells + numVertices, comm, PETSC_ERR_ARG_WRONG, "Invalid ExodusII side %d in set %d maps to point %" PetscInt_FMT " which is not a face", f, fs, faces[0]);
    1887:         PetscCall(DMSetLabelValue_Fast(*dm, &faceSets, "Face Sets", faces[0], fs_id[fs]));
    1888:         /* Only add the label if one has been detected for this side set. */
    1889:         if (!fs_name_err) PetscCall(DMSetLabelValue(*dm, fs_name, faces[0], fs_id[fs]));
    1890:         PetscCall(DMPlexRestoreJoin(*dm, faceSize, faceVertices, &numFaces, &faces));
    1891:       }
    1892:       PetscCall(PetscFree3(fs_vertex_count_list, fs_vertex_list, fs_side_list));
    1893:     }
    1894:     PetscCall(PetscFree(fs_id));
    1895:   }

    1897:   { /* Create Cell/Face/Vertex Sets labels at all processes */
    1898:     enum {
    1899:       n = 3
    1900:     };
    1901:     PetscBool flag[n];

    1903:     flag[0] = cellSets ? PETSC_TRUE : PETSC_FALSE;
    1904:     flag[1] = faceSets ? PETSC_TRUE : PETSC_FALSE;
    1905:     flag[2] = vertSets ? PETSC_TRUE : PETSC_FALSE;
    1906:     PetscCallMPI(MPI_Bcast(flag, n, MPI_C_BOOL, 0, comm));
    1907:     if (flag[0]) PetscCall(DMCreateLabel(*dm, "Cell Sets"));
    1908:     if (flag[1]) PetscCall(DMCreateLabel(*dm, "Face Sets"));
    1909:     if (flag[2]) PetscCall(DMCreateLabel(*dm, "Vertex Sets"));
    1910:   }
    1911:   PetscFunctionReturn(PETSC_SUCCESS);
    1912: }