Actual source code: plexexodusii.c

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

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

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

 15:   Collective; No Fortran Support

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

 20:   Level: intermediate

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

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

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

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

 52:   PetscFunctionBegin;
 53:   if (exo->filename) PetscCall(PetscViewerASCIIPrintf(viewer, "Filename:    %s\n", exo->filename));
 54:   if (exo->exoid) PetscCall(PetscViewerASCIIPrintf(viewer, "exoid:       %d\n", exo->exoid));
 55:   if (exo->btype) PetscCall(PetscViewerASCIIPrintf(viewer, "IO Mode:     %d\n", exo->btype));
 56:   if (exo->order) PetscCall(PetscViewerASCIIPrintf(viewer, "Mesh order:  %" PetscInt_FMT "\n", exo->order));
 57:   PetscFunctionReturn(PETSC_SUCCESS);
 58: }

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

 68: static PetscErrorCode PetscViewerSetUp_ExodusII(PetscViewer viewer)
 69: {
 70:   PetscFunctionBegin;
 71:   PetscFunctionReturn(PETSC_SUCCESS);
 72: }

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

192: /*MC
193:    PETSCVIEWEREXODUSII - A viewer that writes to an Exodus II file

195:   Level: beginner

197: .seealso: `PetscViewerExodusIIOpen()`, `PetscViewerCreate()`, `PETSCVIEWERBINARY`, `PETSCVIEWERHDF5`, `DMView()`,
198:           `PetscViewerFileSetName()`, `PetscViewerFileSetMode()`, `PetscViewerFormat`, `PetscViewerType`, `PetscViewerSetType()`
199: M*/

201: PETSC_EXTERN PetscErrorCode PetscViewerCreate_ExodusII(PetscViewer v)
202: {
203:   PetscViewer_ExodusII *exo;

205:   PetscFunctionBegin;
206:   PetscCall(PetscNew(&exo));

208:   v->data                = (void *)exo;
209:   v->ops->destroy        = PetscViewerDestroy_ExodusII;
210:   v->ops->setfromoptions = PetscViewerSetFromOptions_ExodusII;
211:   v->ops->setup          = PetscViewerSetUp_ExodusII;
212:   v->ops->view           = PetscViewerView_ExodusII;
213:   v->ops->flush          = 0;
214:   exo->btype             = (PetscFileMode)-1;
215:   exo->filename          = 0;
216:   exo->exoid             = -1;

218:   PetscCall(PetscObjectComposeFunction((PetscObject)v, "PetscViewerFileSetName_C", PetscViewerFileSetName_ExodusII));
219:   PetscCall(PetscObjectComposeFunction((PetscObject)v, "PetscViewerFileGetName_C", PetscViewerFileGetName_ExodusII));
220:   PetscCall(PetscObjectComposeFunction((PetscObject)v, "PetscViewerFileSetMode_C", PetscViewerFileSetMode_ExodusII));
221:   PetscCall(PetscObjectComposeFunction((PetscObject)v, "PetscViewerFileGetMode_C", PetscViewerFileGetMode_ExodusII));
222:   PetscCall(PetscObjectComposeFunction((PetscObject)v, "PetscViewerGetId_C", PetscViewerExodusIIGetId_ExodusII));
223:   PetscCall(PetscObjectComposeFunction((PetscObject)v, "PetscViewerSetOrder_C", PetscViewerExodusIISetOrder_ExodusII));
224:   PetscCall(PetscObjectComposeFunction((PetscObject)v, "PetscViewerGetOrder_C", PetscViewerExodusIIGetOrder_ExodusII));
225:   PetscFunctionReturn(PETSC_SUCCESS);
226: }

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

231:   Collective

233:   Input Parameters:
234: + exoid    - the exodus id of a file (obtained from ex_open or ex_create for instance)
235: . obj_type - the type of entity for instance EX_NODAL, EX_ELEM_BLOCK
236: - name     - the name of the result

238:   Output Parameter:
239: . varIndex - the location in the exodus file of the result

241:   Level: beginner

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

249: .seealso: `EXOGetVarIndex()`, `DMPlexView_ExodusII_Internal()`, `VecViewPlex_ExodusII_Nodal_Internal()`, `VecLoadNodal_PlexEXO()`, `VecLoadZonal_PlexEXO()`
250: */
251: static PetscErrorCode EXOGetVarIndex_Internal(int exoid, ex_entity_type obj_type, const char name[], int *varIndex)
252: {
253:   int       num_vars, i, j;
254:   char      ext_name[MAX_STR_LENGTH + 1], var_name[MAX_STR_LENGTH + 1];
255:   const int num_suffix = 5;
256:   char     *suffix[5];
257:   PetscBool flg;

259:   PetscFunctionBegin;
260:   suffix[0] = (char *)"";
261:   suffix[1] = (char *)"_X";
262:   suffix[2] = (char *)"_XX";
263:   suffix[3] = (char *)"_1";
264:   suffix[4] = (char *)"_11";

266:   *varIndex = -1;
267:   PetscCallExternal(ex_get_variable_param, exoid, obj_type, &num_vars);
268:   for (i = 0; i < num_vars; ++i) {
269:     PetscCallExternal(ex_get_variable_name, exoid, obj_type, i + 1, var_name);
270:     for (j = 0; j < num_suffix; ++j) {
271:       PetscCall(PetscStrncpy(ext_name, name, MAX_STR_LENGTH));
272:       PetscCall(PetscStrlcat(ext_name, suffix[j], MAX_STR_LENGTH));
273:       PetscCall(PetscStrcasecmp(ext_name, var_name, &flg));
274:       if (flg) *varIndex = i + 1;
275:     }
276:   }
277:   PetscFunctionReturn(PETSC_SUCCESS);
278: }

280: /*
281:   DMView_PlexExodusII - Write a `DM` to disk in exodus format

283:   Collective

285:   Input Parameters:
286: + dm     - The dm to be written
287: - viewer - an exodusII viewer

289:   Level: beginner

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

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

298:   `DMPLEX` only represents geometry while most post-processing software expect that a mesh also provides information
299:   on the discretization space. This function assumes that the file represents Lagrange finite elements of order 1 or 2.
300:   The order of the mesh shall be set using `PetscViewerExodusIISetOrder()`
301:   It should be extended to use `PetscFE` objects.

303:   This function will only handle TRI, TET, QUAD, and HEX cells.

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

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

352:   PetscFunctionBegin;
353:   PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
354:   PetscCallMPI(MPI_Comm_rank(comm, &rank));
355:   PetscCallMPI(MPI_Comm_size(comm, &size));

357:   /*
358:     Creating coordSection is a collective operation so we do it somewhat out of sequence
359:   */
360:   PetscCall(PetscSectionCreate(comm, &coordSection));
361:   PetscCall(DMGetCoordinatesLocalSetUp(dm));
362:   /*
363:     Check that all points are on rank 0 since we don't know how to save distributed DM in exodus format
364:   */
365:   PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
366:   PetscCall(DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd));
367:   PetscCall(DMPlexGetDepthStratum(dm, 1, &eStart, &eEnd));
368:   PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd));
369:   numCells    = cEnd - cStart;
370:   numEdges    = eEnd - eStart;
371:   numVertices = vEnd - vStart;
372:   PetscCheck(!(rank && (numCells || numEdges || numVertices)), PETSC_COMM_SELF, PETSC_ERR_SUP, "Writing distributed DM in exodusII format not supported");
373:   if (rank == 0) {
374:     switch (exo->btype) {
375:     case FILE_MODE_READ:
376:     case FILE_MODE_APPEND:
377:     case FILE_MODE_UPDATE:
378:     case FILE_MODE_APPEND_UPDATE:
379:       /* exodusII does not allow writing geometry to an existing file */
380:       SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "cannot add geometry to existing file %s", exo->filename);
381:     case FILE_MODE_WRITE:
382:       /* Create an empty file if one already exists*/
383:       EXO_mode = EX_CLOBBER;
384:   #if defined(PETSC_USE_64BIT_INDICES)
385:       EXO_mode += EX_ALL_INT64_API;
386:   #endif
387:       CPU_word_size = sizeof(PetscReal);
388:       IO_word_size  = sizeof(PetscReal);
389:       exo->exoid    = ex_create(exo->filename, EXO_mode, &CPU_word_size, &IO_word_size);
390:       PetscCheck(exo->exoid >= 0, PETSC_COMM_SELF, PETSC_ERR_LIB, "ex_create failed for %s", exo->filename);

392:       break;
393:     default:
394:       SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ORDER, "Must call PetscViewerFileSetMode() before PetscViewerFileSetName()");
395:     }

397:     /* --- Get DM info --- */
398:     PetscCall(PetscObjectGetName((PetscObject)dm, &dmName));
399:     PetscCall(DMPlexGetDepth(dm, &depth));
400:     PetscCall(DMGetDimension(dm, &dim));
401:     PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
402:     if (depth == 3) {
403:       numFaces = fEnd - fStart;
404:     } else {
405:       numFaces = 0;
406:     }
407:     PetscCall(DMGetLabelSize(dm, "Cell Sets", &num_cs));
408:     PetscCall(DMGetLabelSize(dm, "Vertex Sets", &num_vs));
409:     PetscCall(DMGetLabelSize(dm, "Face Sets", &num_fs));
410:     PetscCall(DMGetCoordinatesLocal(dm, &coord));
411:     PetscCall(DMGetCoordinateDM(dm, &cdm));
412:     if (num_cs > 0) {
413:       PetscCall(DMGetLabel(dm, "Cell Sets", &csLabel));
414:       PetscCall(DMLabelGetValueIS(csLabel, &csIS));
415:       PetscCall(ISGetIndices(csIS, &csIdx));
416:     }
417:     PetscCall(PetscMalloc1(num_cs, &nodes));
418:     /* Set element type for each block and compute total number of nodes */
419:     PetscCall(PetscMalloc1(num_cs, &type));
420:     numNodes = numVertices;

422:     PetscCall(PetscViewerExodusIIGetOrder(viewer, &degree));
423:     if (degree == 2) numNodes += numEdges;
424:     cellsNotInConnectivity = numCells;
425:     for (cs = 0; cs < num_cs; ++cs) {
426:       IS              stratumIS;
427:       const PetscInt *cells;
428:       PetscScalar    *xyz = NULL;
429:       PetscInt        csSize, closureSize;

431:       PetscCall(DMLabelGetStratumIS(csLabel, csIdx[cs], &stratumIS));
432:       PetscCall(ISGetIndices(stratumIS, &cells));
433:       PetscCall(ISGetSize(stratumIS, &csSize));
434:       PetscCall(DMPlexVecGetClosure(cdm, NULL, coord, cells[0], &closureSize, &xyz));
435:       switch (dim) {
436:       case 1:
437:         if (closureSize == 2 * dim) {
438:           type[cs] = SEGMENT;
439:         } 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);
440:       case 2:
441:         if (closureSize == 3 * dim) {
442:           type[cs] = TRI;
443:         } else if (closureSize == 4 * dim) {
444:           type[cs] = QUAD;
445:         } 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);
446:         break;
447:       case 3:
448:         if (closureSize == 4 * dim) {
449:           type[cs] = TET;
450:         } else if (closureSize == 8 * dim) {
451:           type[cs] = HEX;
452:         } 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);
453:         break;
454:       default:
455:         SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Dimension %" PetscInt_FMT " not handled by ExodusII viewer", dim);
456:       }
457:       if ((degree == 2) && (type[cs] == SEGMENT)) numNodes += csSize;
458:       if ((degree == 2) && (type[cs] == QUAD)) numNodes += csSize;
459:       if ((degree == 2) && (type[cs] == HEX)) {
460:         numNodes += csSize;
461:         numNodes += numFaces;
462:       }
463:       PetscCall(DMPlexVecRestoreClosure(cdm, NULL, coord, cells[0], &closureSize, &xyz));
464:       /* Set nodes and Element type */
465:       if (type[cs] == SEGMENT) {
466:         if (degree == 1) nodes[cs] = nodesLineP1;
467:         else if (degree == 2) nodes[cs] = nodesLineP2;
468:       } else if (type[cs] == TRI) {
469:         if (degree == 1) nodes[cs] = nodesTriP1;
470:         else if (degree == 2) nodes[cs] = nodesTriP2;
471:       } else if (type[cs] == QUAD) {
472:         if (degree == 1) nodes[cs] = nodesQuadP1;
473:         else if (degree == 2) nodes[cs] = nodesQuadP2;
474:       } else if (type[cs] == TET) {
475:         if (degree == 1) nodes[cs] = nodesTetP1;
476:         else if (degree == 2) nodes[cs] = nodesTetP2;
477:       } else if (type[cs] == HEX) {
478:         if (degree == 1) nodes[cs] = nodesHexP1;
479:         else if (degree == 2) nodes[cs] = nodesHexP2;
480:       }
481:       /* Compute the number of cells not in the connectivity table */
482:       cellsNotInConnectivity -= nodes[cs][3] * csSize;

484:       PetscCall(ISRestoreIndices(stratumIS, &cells));
485:       PetscCall(ISDestroy(&stratumIS));
486:     }
487:     if (num_cs) PetscCallExternal(ex_put_init, exo->exoid, dmName, dim, numNodes, numCells, num_cs, num_vs, num_fs);
488:     /* --- Connectivity --- */
489:     for (cs = 0; cs < num_cs; ++cs) {
490:       IS              stratumIS;
491:       const PetscInt *cells;
492:       PetscInt       *connect, off = 0;
493:       PetscInt        edgesInClosure = 0, facesInClosure = 0, verticesInClosure = 0;
494:       PetscInt        csSize, c, connectSize, closureSize;
495:       char           *elem_type        = NULL;
496:       char            elem_type_bar2[] = "BAR2", elem_type_bar3[] = "BAR3";
497:       char            elem_type_tri3[] = "TRI3", elem_type_quad4[] = "QUAD4";
498:       char            elem_type_tri6[] = "TRI6", elem_type_quad9[] = "QUAD9";
499:       char            elem_type_tet4[] = "TET4", elem_type_hex8[] = "HEX8";
500:       char            elem_type_tet10[] = "TET10", elem_type_hex27[] = "HEX27";

502:       PetscCall(DMLabelGetStratumIS(csLabel, csIdx[cs], &stratumIS));
503:       PetscCall(ISGetIndices(stratumIS, &cells));
504:       PetscCall(ISGetSize(stratumIS, &csSize));
505:       /* Set Element type */
506:       if (type[cs] == SEGMENT) {
507:         if (degree == 1) elem_type = elem_type_bar2;
508:         else if (degree == 2) elem_type = elem_type_bar3;
509:       } else if (type[cs] == TRI) {
510:         if (degree == 1) elem_type = elem_type_tri3;
511:         else if (degree == 2) elem_type = elem_type_tri6;
512:       } else if (type[cs] == QUAD) {
513:         if (degree == 1) elem_type = elem_type_quad4;
514:         else if (degree == 2) elem_type = elem_type_quad9;
515:       } else if (type[cs] == TET) {
516:         if (degree == 1) elem_type = elem_type_tet4;
517:         else if (degree == 2) elem_type = elem_type_tet10;
518:       } else if (type[cs] == HEX) {
519:         if (degree == 1) elem_type = elem_type_hex8;
520:         else if (degree == 2) elem_type = elem_type_hex27;
521:       }
522:       connectSize = nodes[cs][0] + nodes[cs][1] + nodes[cs][2] + nodes[cs][3];
523:       PetscCall(PetscMalloc1(PetscMax(27, connectSize) * csSize, &connect));
524:       PetscCallExternal(ex_put_block, exo->exoid, EX_ELEM_BLOCK, csIdx[cs], elem_type, csSize, connectSize, 0, 0, 1);
525:       /* Find number of vertices, edges, and faces in the closure */
526:       verticesInClosure = nodes[cs][0];
527:       if (depth > 1) {
528:         if (dim == 2) {
529:           PetscCall(DMPlexGetConeSize(dm, cells[0], &edgesInClosure));
530:         } else if (dim == 3) {
531:           PetscInt *closure = NULL;

533:           PetscCall(DMPlexGetConeSize(dm, cells[0], &facesInClosure));
534:           PetscCall(DMPlexGetTransitiveClosure(dm, cells[0], PETSC_TRUE, &closureSize, &closure));
535:           edgesInClosure = closureSize - facesInClosure - 1 - verticesInClosure;
536:           PetscCall(DMPlexRestoreTransitiveClosure(dm, cells[0], PETSC_TRUE, &closureSize, &closure));
537:         }
538:       }
539:       /* Get connectivity for each cell */
540:       for (c = 0; c < csSize; ++c) {
541:         PetscInt *closure = NULL;
542:         PetscInt  temp, i;

544:         PetscCall(DMPlexGetTransitiveClosure(dm, cells[c], PETSC_TRUE, &closureSize, &closure));
545:         for (i = 0; i < connectSize; ++i) {
546:           if (i < nodes[cs][0]) { /* Vertices */
547:             connect[i + off] = closure[(i + edgesInClosure + facesInClosure + 1) * 2] + 1;
548:             connect[i + off] -= cellsNotInConnectivity;
549:           } else if (i < nodes[cs][0] + nodes[cs][1]) { /* Edges */
550:             connect[i + off] = closure[(i - verticesInClosure + facesInClosure + 1) * 2] + 1;
551:             if (nodes[cs][2] == 0) connect[i + off] -= numFaces;
552:             connect[i + off] -= cellsNotInConnectivity;
553:           } else if (i < nodes[cs][0] + nodes[cs][1] + nodes[cs][3]) { /* Cells */
554:             connect[i + off] = closure[0] + 1;
555:             connect[i + off] -= skipCells;
556:           } else if (i < nodes[cs][0] + nodes[cs][1] + nodes[cs][3] + nodes[cs][2]) { /* Faces */
557:             connect[i + off] = closure[(i - edgesInClosure - verticesInClosure) * 2] + 1;
558:             connect[i + off] -= cellsNotInConnectivity;
559:           } else {
560:             connect[i + off] = -1;
561:           }
562:         }
563:         /* Tetrahedra are inverted */
564:         if (type[cs] == TET) {
565:           temp             = connect[0 + off];
566:           connect[0 + off] = connect[1 + off];
567:           connect[1 + off] = temp;
568:           if (degree == 2) {
569:             temp             = connect[5 + off];
570:             connect[5 + off] = connect[6 + off];
571:             connect[6 + off] = temp;
572:             temp             = connect[7 + off];
573:             connect[7 + off] = connect[8 + off];
574:             connect[8 + off] = temp;
575:           }
576:         }
577:         /* Hexahedra are inverted */
578:         if (type[cs] == HEX) {
579:           temp             = connect[1 + off];
580:           connect[1 + off] = connect[3 + off];
581:           connect[3 + off] = temp;
582:           if (degree == 2) {
583:             temp              = connect[8 + off];
584:             connect[8 + off]  = connect[11 + off];
585:             connect[11 + off] = temp;
586:             temp              = connect[9 + off];
587:             connect[9 + off]  = connect[10 + off];
588:             connect[10 + off] = temp;
589:             temp              = connect[16 + off];
590:             connect[16 + off] = connect[17 + off];
591:             connect[17 + off] = temp;
592:             temp              = connect[18 + off];
593:             connect[18 + off] = connect[19 + off];
594:             connect[19 + off] = temp;

596:             temp              = connect[12 + off];
597:             connect[12 + off] = connect[16 + off];
598:             connect[16 + off] = temp;
599:             temp              = connect[13 + off];
600:             connect[13 + off] = connect[17 + off];
601:             connect[17 + off] = temp;
602:             temp              = connect[14 + off];
603:             connect[14 + off] = connect[18 + off];
604:             connect[18 + off] = temp;
605:             temp              = connect[15 + off];
606:             connect[15 + off] = connect[19 + off];
607:             connect[19 + off] = temp;

609:             temp              = connect[23 + off];
610:             connect[23 + off] = connect[26 + off];
611:             connect[26 + off] = temp;
612:             temp              = connect[24 + off];
613:             connect[24 + off] = connect[25 + off];
614:             connect[25 + off] = temp;
615:             temp              = connect[25 + off];
616:             connect[25 + off] = connect[26 + off];
617:             connect[26 + off] = temp;
618:           }
619:         }
620:         off += connectSize;
621:         PetscCall(DMPlexRestoreTransitiveClosure(dm, cells[c], PETSC_TRUE, &closureSize, &closure));
622:       }
623:       PetscCallExternal(ex_put_conn, exo->exoid, EX_ELEM_BLOCK, csIdx[cs], connect, 0, 0);
624:       skipCells += (nodes[cs][3] == 0) * csSize;
625:       PetscCall(PetscFree(connect));
626:       PetscCall(ISRestoreIndices(stratumIS, &cells));
627:       PetscCall(ISDestroy(&stratumIS));
628:     }
629:     PetscCall(PetscFree(type));
630:     /* --- Coordinates --- */
631:     PetscCall(PetscSectionSetChart(coordSection, pStart, pEnd));
632:     if (num_cs) {
633:       for (d = 0; d < depth; ++d) {
634:         PetscCall(DMPlexGetDepthStratum(dm, d, &pStart, &pEnd));
635:         for (p = pStart; p < pEnd; ++p) PetscCall(PetscSectionSetDof(coordSection, p, nodes[0][d] > 0));
636:       }
637:     }
638:     for (cs = 0; cs < num_cs; ++cs) {
639:       IS              stratumIS;
640:       const PetscInt *cells;
641:       PetscInt        csSize, c;

643:       PetscCall(DMLabelGetStratumIS(csLabel, csIdx[cs], &stratumIS));
644:       PetscCall(ISGetIndices(stratumIS, &cells));
645:       PetscCall(ISGetSize(stratumIS, &csSize));
646:       for (c = 0; c < csSize; ++c) PetscCall(PetscSectionSetDof(coordSection, cells[c], nodes[cs][3] > 0));
647:       PetscCall(ISRestoreIndices(stratumIS, &cells));
648:       PetscCall(ISDestroy(&stratumIS));
649:     }
650:     if (num_cs) {
651:       PetscCall(ISRestoreIndices(csIS, &csIdx));
652:       PetscCall(ISDestroy(&csIS));
653:     }
654:     PetscCall(PetscFree(nodes));
655:     PetscCall(PetscSectionSetUp(coordSection));
656:     if (numNodes) {
657:       const char  *coordNames[3] = {"x", "y", "z"};
658:       PetscScalar *closure, *cval;
659:       PetscReal   *coords;
660:       PetscInt     hasDof, n = 0;

662:       /* There can't be more than 24 values in the closure of a point for the coord coordSection */
663:       PetscCall(PetscCalloc3(numNodes * 3, &coords, dim, &cval, 24, &closure));
664:       PetscCall(DMGetCoordinatesLocalNoncollective(dm, &coord));
665:       PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
666:       for (p = pStart; p < pEnd; ++p) {
667:         PetscCall(PetscSectionGetDof(coordSection, p, &hasDof));
668:         if (hasDof) {
669:           PetscInt closureSize = 24, j;

671:           PetscCall(DMPlexVecGetClosure(cdm, NULL, coord, p, &closureSize, &closure));
672:           for (d = 0; d < dim; ++d) {
673:             cval[d] = 0.0;
674:             for (j = 0; j < closureSize / dim; j++) cval[d] += closure[j * dim + d];
675:             coords[d * numNodes + n] = PetscRealPart(cval[d]) * dim / closureSize;
676:           }
677:           ++n;
678:         }
679:       }
680:       PetscCallExternal(ex_put_coord, exo->exoid, &coords[0 * numNodes], &coords[1 * numNodes], &coords[2 * numNodes]);
681:       PetscCall(PetscFree3(coords, cval, closure));
682:       PetscCallExternal(ex_put_coord_names, exo->exoid, (char **)coordNames);
683:     }

685:     /* --- Node Sets/Vertex Sets --- */
686:     PetscCall(DMHasLabel(dm, "Vertex Sets", &hasLabel));
687:     if (hasLabel) {
688:       PetscInt        i, vs, vsSize;
689:       const PetscInt *vsIdx, *vertices;
690:       PetscInt       *nodeList;
691:       IS              vsIS, stratumIS;
692:       DMLabel         vsLabel;
693:       PetscCall(DMGetLabel(dm, "Vertex Sets", &vsLabel));
694:       PetscCall(DMLabelGetValueIS(vsLabel, &vsIS));
695:       PetscCall(ISGetIndices(vsIS, &vsIdx));
696:       for (vs = 0; vs < num_vs; ++vs) {
697:         PetscCall(DMLabelGetStratumIS(vsLabel, vsIdx[vs], &stratumIS));
698:         PetscCall(ISGetIndices(stratumIS, &vertices));
699:         PetscCall(ISGetSize(stratumIS, &vsSize));
700:         PetscCall(PetscMalloc1(vsSize, &nodeList));
701:         for (i = 0; i < vsSize; ++i) nodeList[i] = vertices[i] - skipCells + 1;
702:         PetscCallExternal(ex_put_set_param, exo->exoid, EX_NODE_SET, vsIdx[vs], vsSize, 0);
703:         PetscCallExternal(ex_put_set, exo->exoid, EX_NODE_SET, vsIdx[vs], nodeList, NULL);
704:         PetscCall(ISRestoreIndices(stratumIS, &vertices));
705:         PetscCall(ISDestroy(&stratumIS));
706:         PetscCall(PetscFree(nodeList));
707:       }
708:       PetscCall(ISRestoreIndices(vsIS, &vsIdx));
709:       PetscCall(ISDestroy(&vsIS));
710:     }
711:     /* --- Side Sets/Face Sets --- */
712:     PetscCall(DMHasLabel(dm, "Face Sets", &hasLabel));
713:     if (hasLabel) {
714:       PetscInt        i, j, fs, fsSize;
715:       const PetscInt *fsIdx, *faces;
716:       IS              fsIS, stratumIS;
717:       DMLabel         fsLabel;
718:       PetscInt        numPoints, *points;
719:       PetscInt        elem_list_size = 0;
720:       PetscInt       *elem_list, *elem_ind, *side_list;

722:       PetscCall(DMGetLabel(dm, "Face Sets", &fsLabel));
723:       /* Compute size of Node List and Element List */
724:       PetscCall(DMLabelGetValueIS(fsLabel, &fsIS));
725:       PetscCall(ISGetIndices(fsIS, &fsIdx));
726:       for (fs = 0; fs < num_fs; ++fs) {
727:         PetscCall(DMLabelGetStratumIS(fsLabel, fsIdx[fs], &stratumIS));
728:         PetscCall(ISGetSize(stratumIS, &fsSize));
729:         elem_list_size += fsSize;
730:         PetscCall(ISDestroy(&stratumIS));
731:       }
732:       if (num_fs) {
733:         PetscCall(PetscMalloc3(num_fs, &elem_ind, elem_list_size, &elem_list, elem_list_size, &side_list));
734:         elem_ind[0] = 0;
735:         for (fs = 0; fs < num_fs; ++fs) {
736:           PetscCall(DMLabelGetStratumIS(fsLabel, fsIdx[fs], &stratumIS));
737:           PetscCall(ISGetIndices(stratumIS, &faces));
738:           PetscCall(ISGetSize(stratumIS, &fsSize));
739:           /* Set Parameters */
740:           PetscCallExternal(ex_put_set_param, exo->exoid, EX_SIDE_SET, fsIdx[fs], fsSize, 0);
741:           /* Indices */
742:           if (fs < num_fs - 1) elem_ind[fs + 1] = elem_ind[fs] + fsSize;

744:           for (i = 0; i < fsSize; ++i) {
745:             /* Element List */
746:             points = NULL;
747:             PetscCall(DMPlexGetTransitiveClosure(dm, faces[i], PETSC_FALSE, &numPoints, &points));
748:             elem_list[elem_ind[fs] + i] = points[2] + 1;
749:             PetscCall(DMPlexRestoreTransitiveClosure(dm, faces[i], PETSC_FALSE, &numPoints, &points));

751:             /* Side List */
752:             points = NULL;
753:             PetscCall(DMPlexGetTransitiveClosure(dm, elem_list[elem_ind[fs] + i] - 1, PETSC_TRUE, &numPoints, &points));
754:             for (j = 1; j < numPoints; ++j) {
755:               if (points[j * 2] == faces[i]) break;
756:             }
757:             /* Convert HEX sides */
758:             if (numPoints == 27) {
759:               if (j == 1) {
760:                 j = 5;
761:               } else if (j == 2) {
762:                 j = 6;
763:               } else if (j == 3) {
764:                 j = 1;
765:               } else if (j == 4) {
766:                 j = 3;
767:               } else if (j == 5) {
768:                 j = 2;
769:               } else if (j == 6) {
770:                 j = 4;
771:               }
772:             }
773:             /* Convert TET sides */
774:             if (numPoints == 15) {
775:               --j;
776:               if (j == 0) j = 4;
777:             }
778:             side_list[elem_ind[fs] + i] = j;
779:             PetscCall(DMPlexRestoreTransitiveClosure(dm, elem_list[elem_ind[fs] + i] - 1, PETSC_TRUE, &numPoints, &points));
780:           }
781:           PetscCall(ISRestoreIndices(stratumIS, &faces));
782:           PetscCall(ISDestroy(&stratumIS));
783:         }
784:         PetscCall(ISRestoreIndices(fsIS, &fsIdx));
785:         PetscCall(ISDestroy(&fsIS));

787:         /* Put side sets */
788:         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]]);
789:         PetscCall(PetscFree3(elem_ind, elem_list, side_list));
790:       }
791:     }
792:     /*
793:       close the exodus file
794:     */
795:     ex_close(exo->exoid);
796:     exo->exoid = -1;
797:   }
798:   PetscCall(PetscSectionDestroy(&coordSection));

800:   /*
801:     reopen the file in parallel
802:   */
803:   EXO_mode = EX_WRITE;
804:   #if defined(PETSC_USE_64BIT_INDICES)
805:   EXO_mode += EX_ALL_INT64_API;
806:   #endif
807:   CPU_word_size = sizeof(PetscReal);
808:   IO_word_size  = sizeof(PetscReal);
809:   exo->exoid    = ex_open_par(exo->filename, EXO_mode, &CPU_word_size, &IO_word_size, &EXO_version, comm, MPI_INFO_NULL);
810:   PetscCheck(exo->exoid >= 0, PETSC_COMM_SELF, PETSC_ERR_LIB, "ex_open_par failed for %s", exo->filename);
811:   PetscFunctionReturn(PETSC_SUCCESS);
812: }

814: static PetscErrorCode VecViewPlex_ExodusII_Nodal_Internal(Vec v, int exoid, int step, int offset);
815: static PetscErrorCode VecViewPlex_ExodusII_Zonal_Internal(Vec v, int exoid, int step, int offset);
816: static PetscErrorCode VecLoadPlex_ExodusII_Nodal_Internal(Vec v, int exoid, int step, int offset);
817: static PetscErrorCode VecLoadPlex_ExodusII_Zonal_Internal(Vec v, int exoid, int step, int offset);

819: /*
820:   VecView_PlexExodusII_Internal - Write a Vec corresponding in an exodus file

822:   Collective

824:   Input Parameters:
825: + v  - The vector to be written
826: - viewer - The PetscViewerExodusII viewer associate to an exodus file

828:   Level: beginner

830:   Notes:
831:   The exodus result variable index is obtained by comparing the Vec name and the
832:   names of variables declared in the exodus file. For instance for a Vec named "V"
833:   the location in the exodus file will be the first match of "V", "V_X", "V_XX", "V_1", or "V_11"
834:   amongst all variables.
835:   In the event where a nodal and zonal variable both match, the function will return an error instead of
836:   possibly corrupting the file

838: .seealso: `EXOGetVarIndex_Internal()`, `DMPlexView_ExodusII()`, `VecView_PlexExodusII()`
839: @*/
840: PetscErrorCode VecView_PlexExodusII_Internal(Vec v, PetscViewer viewer)
841: {
842:   DM          dm;
843:   MPI_Comm    comm;
844:   PetscMPIInt rank;

846:   int         exoid, offsetN = 0, offsetZ = 0;
847:   const char *vecname;
848:   PetscInt    step;

850:   PetscFunctionBegin;
851:   PetscCall(PetscObjectGetComm((PetscObject)v, &comm));
852:   PetscCallMPI(MPI_Comm_rank(comm, &rank));
853:   PetscCall(PetscViewerExodusIIGetId(viewer, &exoid));
854:   PetscCall(VecGetDM(v, &dm));
855:   PetscCall(PetscObjectGetName((PetscObject)v, &vecname));

857:   PetscCall(DMGetOutputSequenceNumber(dm, &step, NULL));
858:   PetscCall(EXOGetVarIndex_Internal(exoid, EX_NODAL, vecname, &offsetN));
859:   PetscCall(EXOGetVarIndex_Internal(exoid, EX_ELEM_BLOCK, vecname, &offsetZ));
860:   PetscCheck(offsetN > 0 || offsetZ > 0, comm, PETSC_ERR_FILE_UNEXPECTED, "Found both nodal and zonal variable %s in exodus file. ", vecname);
861:   if (offsetN > 0) {
862:     PetscCall(VecViewPlex_ExodusII_Nodal_Internal(v, exoid, (int)step + 1, offsetN));
863:   } else if (offsetZ > 0) {
864:     PetscCall(VecViewPlex_ExodusII_Zonal_Internal(v, exoid, (int)step + 1, offsetZ));
865:   } else SETERRQ(comm, PETSC_ERR_FILE_UNEXPECTED, "Could not find nodal or zonal variable %s in exodus file. ", vecname);
866:   PetscFunctionReturn(PETSC_SUCCESS);
867: }

869: /*
870:   VecLoad_PlexExodusII_Internal - Write a Vec corresponding in an exodus file

872:   Collective

874:   Input Parameters:
875: + v  - The vector to be written
876: - viewer - The PetscViewerExodusII viewer associate to an exodus file

878:   Level: beginner

880:   Notes:
881:   The exodus result variable index is obtained by comparing the Vec name and the
882:   names of variables declared in the exodus file. For instance for a Vec named "V"
883:   the location in the exodus file will be the first match of "V", "V_X", "V_XX", "V_1", or "V_11"
884:   amongst all variables.
885:   In the event where a nodal and zonal variable both match, the function will return an error instead of
886:   possibly corrupting the file

888: .seealso: `EXOGetVarIndex_Internal()`, `DMPlexView_ExodusII()`, `VecView_PlexExodusII()`
889: @*/
890: PetscErrorCode VecLoad_PlexExodusII_Internal(Vec v, PetscViewer viewer)
891: {
892:   DM          dm;
893:   MPI_Comm    comm;
894:   PetscMPIInt rank;

896:   int         exoid, offsetN = 0, offsetZ = 0;
897:   const char *vecname;
898:   PetscInt    step;

900:   PetscFunctionBegin;
901:   PetscCall(PetscObjectGetComm((PetscObject)v, &comm));
902:   PetscCallMPI(MPI_Comm_rank(comm, &rank));
903:   PetscCall(PetscViewerExodusIIGetId(viewer, &exoid));
904:   PetscCall(VecGetDM(v, &dm));
905:   PetscCall(PetscObjectGetName((PetscObject)v, &vecname));

907:   PetscCall(DMGetOutputSequenceNumber(dm, &step, NULL));
908:   PetscCall(EXOGetVarIndex_Internal(exoid, EX_NODAL, vecname, &offsetN));
909:   PetscCall(EXOGetVarIndex_Internal(exoid, EX_ELEM_BLOCK, vecname, &offsetZ));
910:   PetscCheck(offsetN > 0 || offsetZ > 0, comm, PETSC_ERR_FILE_UNEXPECTED, "Found both nodal and zonal variable %s in exodus file. ", vecname);
911:   if (offsetN > 0) PetscCall(VecLoadPlex_ExodusII_Nodal_Internal(v, exoid, (int)step + 1, offsetN));
912:   else if (offsetZ > 0) PetscCall(VecLoadPlex_ExodusII_Zonal_Internal(v, exoid, (int)step + 1, offsetZ));
913:   else SETERRQ(comm, PETSC_ERR_FILE_UNEXPECTED, "Could not find nodal or zonal variable %s in exodus file. ", vecname);
914:   PetscFunctionReturn(PETSC_SUCCESS);
915: }

917: /*
918:   VecViewPlex_ExodusII_Nodal_Internal - Write a Vec corresponding to a nodal field to an exodus file

920:   Collective

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

928:   Level: beginner

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

936: .seealso: `EXOGetVarIndex_Internal()`, `DMPlexView_ExodusII_Internal()`, `VecLoadNodal_PlexEXO()`, `VecViewZonal_PlexEXO()`, `VecLoadZonal_PlexEXO()`
937: @*/
938: static PetscErrorCode VecViewPlex_ExodusII_Nodal_Internal(Vec v, int exoid, int step, int offset)
939: {
940:   MPI_Comm           comm;
941:   PetscMPIInt        size;
942:   DM                 dm;
943:   Vec                vNatural, vComp;
944:   const PetscScalar *varray;
945:   PetscInt           xs, xe, bs;
946:   PetscBool          useNatural;

948:   PetscFunctionBegin;
949:   PetscCall(PetscObjectGetComm((PetscObject)v, &comm));
950:   PetscCallMPI(MPI_Comm_size(comm, &size));
951:   PetscCall(VecGetDM(v, &dm));
952:   PetscCall(DMGetUseNatural(dm, &useNatural));
953:   useNatural = useNatural && size > 1 ? PETSC_TRUE : PETSC_FALSE;
954:   if (useNatural) {
955:     PetscCall(DMPlexCreateNaturalVector(dm, &vNatural));
956:     PetscCall(DMPlexGlobalToNaturalBegin(dm, v, vNatural));
957:     PetscCall(DMPlexGlobalToNaturalEnd(dm, v, vNatural));
958:   } else {
959:     vNatural = v;
960:   }

962:   /* Write local chunk of the result in the exodus file
963:      exodus stores each component of a vector-valued field as a separate variable.
964:      We assume that they are stored sequentially */
965:   PetscCall(VecGetOwnershipRange(vNatural, &xs, &xe));
966:   PetscCall(VecGetBlockSize(vNatural, &bs));
967:   if (bs == 1) {
968:     PetscCall(VecGetArrayRead(vNatural, &varray));
969:     PetscCallExternal(ex_put_partial_var, exoid, step, EX_NODAL, offset, 1, xs + 1, xe - xs, varray);
970:     PetscCall(VecRestoreArrayRead(vNatural, &varray));
971:   } else {
972:     IS       compIS;
973:     PetscInt c;

975:     PetscCall(ISCreateStride(comm, (xe - xs) / bs, xs, bs, &compIS));
976:     for (c = 0; c < bs; ++c) {
977:       PetscCall(ISStrideSetStride(compIS, (xe - xs) / bs, xs + c, bs));
978:       PetscCall(VecGetSubVector(vNatural, compIS, &vComp));
979:       PetscCall(VecGetArrayRead(vComp, &varray));
980:       PetscCallExternal(ex_put_partial_var, exoid, step, EX_NODAL, offset + c, 1, xs / bs + 1, (xe - xs) / bs, varray);
981:       PetscCall(VecRestoreArrayRead(vComp, &varray));
982:       PetscCall(VecRestoreSubVector(vNatural, compIS, &vComp));
983:     }
984:     PetscCall(ISDestroy(&compIS));
985:   }
986:   if (useNatural) PetscCall(VecDestroy(&vNatural));
987:   PetscFunctionReturn(PETSC_SUCCESS);
988: }

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

993:   Collective

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

1001:   Level: beginner

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

1009: .seealso: `EXOGetVarIndex_Internal()`, `DMPlexView_ExodusII_Internal()`, `VecViewPlex_ExodusII_Nodal_Internal()`, `VecViewZonal_PlexEXO()`, `VecLoadZonal_PlexEXO()`
1010: */
1011: static PetscErrorCode VecLoadPlex_ExodusII_Nodal_Internal(Vec v, int exoid, int step, int offset)
1012: {
1013:   MPI_Comm     comm;
1014:   PetscMPIInt  size;
1015:   DM           dm;
1016:   Vec          vNatural, vComp;
1017:   PetscScalar *varray;
1018:   PetscInt     xs, xe, bs;
1019:   PetscBool    useNatural;

1021:   PetscFunctionBegin;
1022:   PetscCall(PetscObjectGetComm((PetscObject)v, &comm));
1023:   PetscCallMPI(MPI_Comm_size(comm, &size));
1024:   PetscCall(VecGetDM(v, &dm));
1025:   PetscCall(DMGetUseNatural(dm, &useNatural));
1026:   useNatural = useNatural && size > 1 ? PETSC_TRUE : PETSC_FALSE;
1027:   if (useNatural) PetscCall(DMPlexCreateNaturalVector(dm, &vNatural));
1028:   else vNatural = v;

1030:   /* Read local chunk from the file */
1031:   PetscCall(VecGetOwnershipRange(vNatural, &xs, &xe));
1032:   PetscCall(VecGetBlockSize(vNatural, &bs));
1033:   if (bs == 1) {
1034:     PetscCall(VecGetArray(vNatural, &varray));
1035:     PetscCallExternal(ex_get_partial_var, exoid, step, EX_NODAL, offset, 1, xs + 1, xe - xs, varray);
1036:     PetscCall(VecRestoreArray(vNatural, &varray));
1037:   } else {
1038:     IS       compIS;
1039:     PetscInt c;

1041:     PetscCall(ISCreateStride(comm, (xe - xs) / bs, xs, bs, &compIS));
1042:     for (c = 0; c < bs; ++c) {
1043:       PetscCall(ISStrideSetStride(compIS, (xe - xs) / bs, xs + c, bs));
1044:       PetscCall(VecGetSubVector(vNatural, compIS, &vComp));
1045:       PetscCall(VecGetArray(vComp, &varray));
1046:       PetscCallExternal(ex_get_partial_var, exoid, step, EX_NODAL, offset + c, 1, xs / bs + 1, (xe - xs) / bs, varray);
1047:       PetscCall(VecRestoreArray(vComp, &varray));
1048:       PetscCall(VecRestoreSubVector(vNatural, compIS, &vComp));
1049:     }
1050:     PetscCall(ISDestroy(&compIS));
1051:   }
1052:   if (useNatural) {
1053:     PetscCall(DMPlexNaturalToGlobalBegin(dm, vNatural, v));
1054:     PetscCall(DMPlexNaturalToGlobalEnd(dm, vNatural, v));
1055:     PetscCall(VecDestroy(&vNatural));
1056:   }
1057:   PetscFunctionReturn(PETSC_SUCCESS);
1058: }

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

1063:   Collective

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

1071:   Level: beginner

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

1079: .seealso: `EXOGetVarIndex_Internal()`, `DMPlexView_ExodusII_Internal()`, `VecViewPlex_ExodusII_Nodal_Internal()`, `VecLoadPlex_ExodusII_Nodal_Internal()`, `VecLoadPlex_ExodusII_Zonal_Internal()`
1080: */
1081: static PetscErrorCode VecViewPlex_ExodusII_Zonal_Internal(Vec v, int exoid, int step, int offset)
1082: {
1083:   MPI_Comm           comm;
1084:   PetscMPIInt        size;
1085:   DM                 dm;
1086:   Vec                vNatural, vComp;
1087:   const PetscScalar *varray;
1088:   PetscInt           xs, xe, bs;
1089:   PetscBool          useNatural;
1090:   IS                 compIS;
1091:   PetscInt          *csSize, *csID;
1092:   PetscInt           numCS, set, csxs = 0;

1094:   PetscFunctionBegin;
1095:   PetscCall(PetscObjectGetComm((PetscObject)v, &comm));
1096:   PetscCallMPI(MPI_Comm_size(comm, &size));
1097:   PetscCall(VecGetDM(v, &dm));
1098:   PetscCall(DMGetUseNatural(dm, &useNatural));
1099:   useNatural = useNatural && size > 1 ? PETSC_TRUE : PETSC_FALSE;
1100:   if (useNatural) {
1101:     PetscCall(DMPlexCreateNaturalVector(dm, &vNatural));
1102:     PetscCall(DMPlexGlobalToNaturalBegin(dm, v, vNatural));
1103:     PetscCall(DMPlexGlobalToNaturalEnd(dm, v, vNatural));
1104:   } else {
1105:     vNatural = v;
1106:   }

1108:   /* Write local chunk of the result in the exodus file
1109:      exodus stores each component of a vector-valued field as a separate variable.
1110:      We assume that they are stored sequentially
1111:      Zonal variables are accessed one element block at a time, so we loop through the cell sets,
1112:      but once the vector has been reordered to natural size, we cannot use the label information
1113:      to figure out what to save where. */
1114:   numCS = ex_inquire_int(exoid, EX_INQ_ELEM_BLK);
1115:   PetscCall(PetscMalloc2(numCS, &csID, numCS, &csSize));
1116:   PetscCallExternal(ex_get_ids, exoid, EX_ELEM_BLOCK, csID);
1117:   for (set = 0; set < numCS; ++set) {
1118:     ex_block block;

1120:     block.id   = csID[set];
1121:     block.type = EX_ELEM_BLOCK;
1122:     PetscCallExternal(ex_get_block_param, exoid, &block);
1123:     csSize[set] = block.num_entry;
1124:   }
1125:   PetscCall(VecGetOwnershipRange(vNatural, &xs, &xe));
1126:   PetscCall(VecGetBlockSize(vNatural, &bs));
1127:   if (bs > 1) PetscCall(ISCreateStride(comm, (xe - xs) / bs, xs, bs, &compIS));
1128:   for (set = 0; set < numCS; set++) {
1129:     PetscInt csLocalSize, c;

1131:     /* range of indices for set setID[set]: csxs:csxs + csSize[set]-1
1132:        local slice of zonal values:         xs/bs,xm/bs-1
1133:        intersection:                        max(xs/bs,csxs),min(xm/bs-1,csxs + csSize[set]-1) */
1134:     csLocalSize = PetscMax(0, PetscMin(xe / bs, csxs + csSize[set]) - PetscMax(xs / bs, csxs));
1135:     if (bs == 1) {
1136:       PetscCall(VecGetArrayRead(vNatural, &varray));
1137:       PetscCallExternal(ex_put_partial_var, exoid, step, EX_ELEM_BLOCK, offset, csID[set], PetscMax(xs - csxs, 0) + 1, csLocalSize, &varray[PetscMax(0, csxs - xs)]);
1138:       PetscCall(VecRestoreArrayRead(vNatural, &varray));
1139:     } else {
1140:       for (c = 0; c < bs; ++c) {
1141:         PetscCall(ISStrideSetStride(compIS, (xe - xs) / bs, xs + c, bs));
1142:         PetscCall(VecGetSubVector(vNatural, compIS, &vComp));
1143:         PetscCall(VecGetArrayRead(vComp, &varray));
1144:         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)]);
1145:         PetscCall(VecRestoreArrayRead(vComp, &varray));
1146:         PetscCall(VecRestoreSubVector(vNatural, compIS, &vComp));
1147:       }
1148:     }
1149:     csxs += csSize[set];
1150:   }
1151:   PetscCall(PetscFree2(csID, csSize));
1152:   if (bs > 1) PetscCall(ISDestroy(&compIS));
1153:   if (useNatural) PetscCall(VecDestroy(&vNatural));
1154:   PetscFunctionReturn(PETSC_SUCCESS);
1155: }

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

1160:   Collective

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

1168:   Level: beginner

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

1176: .seealso: `EXOGetVarIndex_Internal()`, `DMPlexView_ExodusII_Internal()`, `VecViewPlex_ExodusII_Nodal_Internal()`, `VecLoadPlex_ExodusII_Nodal_Internal()`, `VecLoadPlex_ExodusII_Zonal_Internal()`
1177: */
1178: static PetscErrorCode VecLoadPlex_ExodusII_Zonal_Internal(Vec v, int exoid, int step, int offset)
1179: {
1180:   MPI_Comm     comm;
1181:   PetscMPIInt  size;
1182:   DM           dm;
1183:   Vec          vNatural, vComp;
1184:   PetscScalar *varray;
1185:   PetscInt     xs, xe, bs;
1186:   PetscBool    useNatural;
1187:   IS           compIS;
1188:   PetscInt    *csSize, *csID;
1189:   PetscInt     numCS, set, csxs = 0;

1191:   PetscFunctionBegin;
1192:   PetscCall(PetscObjectGetComm((PetscObject)v, &comm));
1193:   PetscCallMPI(MPI_Comm_size(comm, &size));
1194:   PetscCall(VecGetDM(v, &dm));
1195:   PetscCall(DMGetUseNatural(dm, &useNatural));
1196:   useNatural = useNatural && size > 1 ? PETSC_TRUE : PETSC_FALSE;
1197:   if (useNatural) PetscCall(DMPlexCreateNaturalVector(dm, &vNatural));
1198:   else vNatural = v;

1200:   /* Read local chunk of the result in the exodus file
1201:      exodus stores each component of a vector-valued field as a separate variable.
1202:      We assume that they are stored sequentially
1203:      Zonal variables are accessed one element block at a time, so we loop through the cell sets,
1204:      but once the vector has been reordered to natural size, we cannot use the label information
1205:      to figure out what to save where. */
1206:   numCS = ex_inquire_int(exoid, EX_INQ_ELEM_BLK);
1207:   PetscCall(PetscMalloc2(numCS, &csID, numCS, &csSize));
1208:   PetscCallExternal(ex_get_ids, exoid, EX_ELEM_BLOCK, csID);
1209:   for (set = 0; set < numCS; ++set) {
1210:     ex_block block;

1212:     block.id   = csID[set];
1213:     block.type = EX_ELEM_BLOCK;
1214:     PetscCallExternal(ex_get_block_param, exoid, &block);
1215:     csSize[set] = block.num_entry;
1216:   }
1217:   PetscCall(VecGetOwnershipRange(vNatural, &xs, &xe));
1218:   PetscCall(VecGetBlockSize(vNatural, &bs));
1219:   if (bs > 1) PetscCall(ISCreateStride(comm, (xe - xs) / bs, xs, bs, &compIS));
1220:   for (set = 0; set < numCS; ++set) {
1221:     PetscInt csLocalSize, c;

1223:     /* range of indices for set setID[set]: csxs:csxs + csSize[set]-1
1224:        local slice of zonal values:         xs/bs,xm/bs-1
1225:        intersection:                        max(xs/bs,csxs),min(xm/bs-1,csxs + csSize[set]-1) */
1226:     csLocalSize = PetscMax(0, PetscMin(xe / bs, csxs + csSize[set]) - PetscMax(xs / bs, csxs));
1227:     if (bs == 1) {
1228:       PetscCall(VecGetArray(vNatural, &varray));
1229:       PetscCallExternal(ex_get_partial_var, exoid, step, EX_ELEM_BLOCK, offset, csID[set], PetscMax(xs - csxs, 0) + 1, csLocalSize, &varray[PetscMax(0, csxs - xs)]);
1230:       PetscCall(VecRestoreArray(vNatural, &varray));
1231:     } else {
1232:       for (c = 0; c < bs; ++c) {
1233:         PetscCall(ISStrideSetStride(compIS, (xe - xs) / bs, xs + c, bs));
1234:         PetscCall(VecGetSubVector(vNatural, compIS, &vComp));
1235:         PetscCall(VecGetArray(vComp, &varray));
1236:         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)]);
1237:         PetscCall(VecRestoreArray(vComp, &varray));
1238:         PetscCall(VecRestoreSubVector(vNatural, compIS, &vComp));
1239:       }
1240:     }
1241:     csxs += csSize[set];
1242:   }
1243:   PetscCall(PetscFree2(csID, csSize));
1244:   if (bs > 1) PetscCall(ISDestroy(&compIS));
1245:   if (useNatural) {
1246:     PetscCall(DMPlexNaturalToGlobalBegin(dm, vNatural, v));
1247:     PetscCall(DMPlexNaturalToGlobalEnd(dm, vNatural, v));
1248:     PetscCall(VecDestroy(&vNatural));
1249:   }
1250:   PetscFunctionReturn(PETSC_SUCCESS);
1251: }
1252: #endif

1254: /*@
1255:   PetscViewerExodusIIGetId - Get the file id of the `PETSCVIEWEREXODUSII` file

1257:   Logically Collective

1259:   Input Parameter:
1260: . viewer - the `PetscViewer`

1262:   Output Parameter:
1263: . exoid - The ExodusII file id

1265:   Level: intermediate

1267: .seealso: `PETSCVIEWEREXODUSII`, `PetscViewer`, `PetscViewerFileSetMode()`, `PetscViewerCreate()`, `PetscViewerSetType()`, `PetscViewerBinaryOpen()`
1268: @*/
1269: PetscErrorCode PetscViewerExodusIIGetId(PetscViewer viewer, int *exoid)
1270: {
1271:   PetscFunctionBegin;
1273:   PetscTryMethod(viewer, "PetscViewerGetId_C", (PetscViewer, int *), (viewer, exoid));
1274:   PetscFunctionReturn(PETSC_SUCCESS);
1275: }

1277: /*@
1278:   PetscViewerExodusIISetOrder - Set the elements order in the exodusII file.

1280:   Collective

1282:   Input Parameters:
1283: + viewer - the `PETSCVIEWEREXODUSII` viewer
1284: - order  - elements order

1286:   Output Parameter:

1288:   Level: beginner

1290: .seealso: `PETSCVIEWEREXODUSII`, `PetscViewer`, `PetscViewerExodusIIGetId()`, `PetscViewerExodusIIGetOrder()`
1291: @*/
1292: PetscErrorCode PetscViewerExodusIISetOrder(PetscViewer viewer, PetscInt order)
1293: {
1294:   PetscFunctionBegin;
1296:   PetscTryMethod(viewer, "PetscViewerSetOrder_C", (PetscViewer, PetscInt), (viewer, order));
1297:   PetscFunctionReturn(PETSC_SUCCESS);
1298: }

1300: /*@
1301:   PetscViewerExodusIIGetOrder - Get the elements order in the exodusII file.

1303:   Collective

1305:   Input Parameters:
1306: + viewer - the `PETSCVIEWEREXODUSII` viewer
1307: - order  - elements order

1309:   Output Parameter:

1311:   Level: beginner

1313: .seealso: `PETSCVIEWEREXODUSII`, `PetscViewer`, `PetscViewerExodusIIGetId()`, `PetscViewerExodusIISetOrder()`
1314: @*/
1315: PetscErrorCode PetscViewerExodusIIGetOrder(PetscViewer viewer, PetscInt *order)
1316: {
1317:   PetscFunctionBegin;
1319:   PetscTryMethod(viewer, "PetscViewerGetOrder_C", (PetscViewer, PetscInt *), (viewer, order));
1320:   PetscFunctionReturn(PETSC_SUCCESS);
1321: }

1323: /*@C
1324:   PetscViewerExodusIIOpen - Opens a file for ExodusII input/output.

1326:   Collective

1328:   Input Parameters:
1329: + comm - MPI communicator
1330: . name - name of file
1331: - type - type of file
1332: .vb
1333:     FILE_MODE_WRITE - create new file for binary output
1334:     FILE_MODE_READ - open existing file for binary input
1335:     FILE_MODE_APPEND - open existing file for binary output
1336: .ve

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

1341:   Level: beginner

1343: .seealso: `PETSCVIEWEREXODUSII`, `PetscViewer`, `PetscViewerPushFormat()`, `PetscViewerDestroy()`,
1344:           `DMLoad()`, `PetscFileMode`, `PetscViewerSetType()`, `PetscViewerFileSetMode()`, `PetscViewerFileSetName()`
1345: @*/
1346: PetscErrorCode PetscViewerExodusIIOpen(MPI_Comm comm, const char name[], PetscFileMode type, PetscViewer *exo)
1347: {
1348:   PetscFunctionBegin;
1349:   PetscCall(PetscViewerCreate(comm, exo));
1350:   PetscCall(PetscViewerSetType(*exo, PETSCVIEWEREXODUSII));
1351:   PetscCall(PetscViewerFileSetMode(*exo, type));
1352:   PetscCall(PetscViewerFileSetName(*exo, name));
1353:   PetscCall(PetscViewerSetFromOptions(*exo));
1354:   PetscFunctionReturn(PETSC_SUCCESS);
1355: }

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

1360:   Collective

1362:   Input Parameters:
1363: + comm        - The MPI communicator
1364: . filename    - The name of the ExodusII file
1365: - interpolate - Create faces and edges in the mesh

1367:   Output Parameter:
1368: . dm - The `DM` object representing the mesh

1370:   Level: beginner

1372: .seealso: [](ch_unstructured), `DM`, `PETSCVIEWEREXODUSII`, `DMPLEX`, `DMCreate()`, `DMPlexCreateExodus()`
1373: @*/
1374: PetscErrorCode DMPlexCreateExodusFromFile(MPI_Comm comm, const char filename[], PetscBool interpolate, DM *dm)
1375: {
1376:   PetscMPIInt rank;
1377: #if defined(PETSC_HAVE_EXODUSII)
1378:   int   CPU_word_size = sizeof(PetscReal), IO_word_size = 0, exoid = -1;
1379:   float version;
1380: #endif

1382:   PetscFunctionBegin;
1383:   PetscAssertPointer(filename, 2);
1384:   PetscCallMPI(MPI_Comm_rank(comm, &rank));
1385: #if defined(PETSC_HAVE_EXODUSII)
1386:   if (rank == 0) {
1387:     exoid = ex_open(filename, EX_READ, &CPU_word_size, &IO_word_size, &version);
1388:     PetscCheck(exoid > 0, PETSC_COMM_SELF, PETSC_ERR_LIB, "ex_open(\"%s\",...) did not return a valid file ID", filename);
1389:   }
1390:   PetscCall(DMPlexCreateExodus(comm, exoid, interpolate, dm));
1391:   if (rank == 0) PetscCallExternal(ex_close, exoid);
1392:   PetscFunctionReturn(PETSC_SUCCESS);
1393: #else
1394:   SETERRQ(comm, PETSC_ERR_SUP, "This method requires ExodusII support. Reconfigure using --download-exodusii");
1395: #endif
1396: }

1398: #if defined(PETSC_HAVE_EXODUSII)
1399: static PetscErrorCode ExodusGetCellType_Internal(const char *elem_type, DMPolytopeType *ct)
1400: {
1401:   PetscBool flg;

1403:   PetscFunctionBegin;
1404:   *ct = DM_POLYTOPE_UNKNOWN;
1405:   PetscCall(PetscStrcmp(elem_type, "BAR2", &flg));
1406:   if (flg) {
1407:     *ct = DM_POLYTOPE_SEGMENT;
1408:     goto done;
1409:   }
1410:   PetscCall(PetscStrcmp(elem_type, "BAR3", &flg));
1411:   if (flg) {
1412:     *ct = DM_POLYTOPE_SEGMENT;
1413:     goto done;
1414:   }
1415:   PetscCall(PetscStrcmp(elem_type, "TRI", &flg));
1416:   if (flg) {
1417:     *ct = DM_POLYTOPE_TRIANGLE;
1418:     goto done;
1419:   }
1420:   PetscCall(PetscStrcmp(elem_type, "TRI3", &flg));
1421:   if (flg) {
1422:     *ct = DM_POLYTOPE_TRIANGLE;
1423:     goto done;
1424:   }
1425:   PetscCall(PetscStrcmp(elem_type, "QUAD", &flg));
1426:   if (flg) {
1427:     *ct = DM_POLYTOPE_QUADRILATERAL;
1428:     goto done;
1429:   }
1430:   PetscCall(PetscStrcmp(elem_type, "QUAD4", &flg));
1431:   if (flg) {
1432:     *ct = DM_POLYTOPE_QUADRILATERAL;
1433:     goto done;
1434:   }
1435:   PetscCall(PetscStrcmp(elem_type, "SHELL4", &flg));
1436:   if (flg) {
1437:     *ct = DM_POLYTOPE_QUADRILATERAL;
1438:     goto done;
1439:   }
1440:   PetscCall(PetscStrcmp(elem_type, "TETRA", &flg));
1441:   if (flg) {
1442:     *ct = DM_POLYTOPE_TETRAHEDRON;
1443:     goto done;
1444:   }
1445:   PetscCall(PetscStrcmp(elem_type, "TET4", &flg));
1446:   if (flg) {
1447:     *ct = DM_POLYTOPE_TETRAHEDRON;
1448:     goto done;
1449:   }
1450:   PetscCall(PetscStrcmp(elem_type, "WEDGE", &flg));
1451:   if (flg) {
1452:     *ct = DM_POLYTOPE_TRI_PRISM;
1453:     goto done;
1454:   }
1455:   PetscCall(PetscStrcmp(elem_type, "HEX", &flg));
1456:   if (flg) {
1457:     *ct = DM_POLYTOPE_HEXAHEDRON;
1458:     goto done;
1459:   }
1460:   PetscCall(PetscStrcmp(elem_type, "HEX8", &flg));
1461:   if (flg) {
1462:     *ct = DM_POLYTOPE_HEXAHEDRON;
1463:     goto done;
1464:   }
1465:   PetscCall(PetscStrcmp(elem_type, "HEXAHEDRON", &flg));
1466:   if (flg) {
1467:     *ct = DM_POLYTOPE_HEXAHEDRON;
1468:     goto done;
1469:   }
1470:   SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unrecognized element type %s", elem_type);
1471: done:
1472:   PetscFunctionReturn(PETSC_SUCCESS);
1473: }
1474: #endif

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

1479:   Collective

1481:   Input Parameters:
1482: + comm        - The MPI communicator
1483: . exoid       - The ExodusII id associated with a exodus file and obtained using ex_open
1484: - interpolate - Create faces and edges in the mesh

1486:   Output Parameter:
1487: . dm - The `DM` object representing the mesh

1489:   Level: beginner

1491: .seealso: [](ch_unstructured), `DM`, `PETSCVIEWEREXODUSII`, `DMPLEX`, `DMCreate()`
1492: @*/
1493: PetscErrorCode DMPlexCreateExodus(MPI_Comm comm, PetscInt exoid, PetscBool interpolate, DM *dm)
1494: {
1495: #if defined(PETSC_HAVE_EXODUSII)
1496:   PetscMPIInt  num_proc, rank;
1497:   DMLabel      cellSets = NULL, faceSets = NULL, vertSets = NULL;
1498:   PetscSection coordSection;
1499:   Vec          coordinates;
1500:   PetscScalar *coords;
1501:   PetscInt     coordSize, v;
1502:   /* Read from ex_get_init() */
1503:   char title[PETSC_MAX_PATH_LEN + 1];
1504:   int  dim = 0, dimEmbed = 0, numVertices = 0, numCells = 0;
1505:   int  num_cs = 0, num_vs = 0, num_fs = 0;
1506: #endif

1508:   PetscFunctionBegin;
1509: #if defined(PETSC_HAVE_EXODUSII)
1510:   PetscCallMPI(MPI_Comm_rank(comm, &rank));
1511:   PetscCallMPI(MPI_Comm_size(comm, &num_proc));
1512:   PetscCall(DMCreate(comm, dm));
1513:   PetscCall(DMSetType(*dm, DMPLEX));
1514:   /* Open EXODUS II file and read basic information on rank 0, then broadcast to all processors */
1515:   if (rank == 0) {
1516:     PetscCall(PetscMemzero(title, PETSC_MAX_PATH_LEN + 1));
1517:     PetscCallExternal(ex_get_init, exoid, title, &dimEmbed, &numVertices, &numCells, &num_cs, &num_vs, &num_fs);
1518:     PetscCheck(num_cs, PETSC_COMM_SELF, PETSC_ERR_SUP, "Exodus file does not contain any cell set");
1519:   }
1520:   PetscCallMPI(MPI_Bcast(title, PETSC_MAX_PATH_LEN + 1, MPI_CHAR, 0, comm));
1521:   PetscCallMPI(MPI_Bcast(&dim, 1, MPI_INT, 0, comm));
1522:   PetscCall(PetscObjectSetName((PetscObject)*dm, title));
1523:   PetscCall(DMPlexSetChart(*dm, 0, numCells + numVertices));
1524:   /*   We do not want this label automatically computed, instead we compute it here */
1525:   PetscCall(DMCreateLabel(*dm, "celltype"));

1527:   /* Read cell sets information */
1528:   if (rank == 0) {
1529:     PetscInt *cone;
1530:     int       c, cs, ncs, c_loc, v, v_loc;
1531:     /* Read from ex_get_elem_blk_ids() */
1532:     int *cs_id, *cs_order;
1533:     /* Read from ex_get_elem_block() */
1534:     char buffer[PETSC_MAX_PATH_LEN + 1];
1535:     int  num_cell_in_set, num_vertex_per_cell, num_hybrid, num_attr;
1536:     /* Read from ex_get_elem_conn() */
1537:     int *cs_connect;

1539:     /* Get cell sets IDs */
1540:     PetscCall(PetscMalloc2(num_cs, &cs_id, num_cs, &cs_order));
1541:     PetscCallExternal(ex_get_ids, exoid, EX_ELEM_BLOCK, cs_id);
1542:     /* Read the cell set connectivity table and build mesh topology
1543:        EXO standard requires that cells in cell sets be numbered sequentially and be pairwise disjoint. */
1544:     /* Check for a hybrid mesh */
1545:     for (cs = 0, num_hybrid = 0; cs < num_cs; ++cs) {
1546:       DMPolytopeType ct;
1547:       char           elem_type[PETSC_MAX_PATH_LEN];

1549:       PetscCall(PetscArrayzero(elem_type, sizeof(elem_type)));
1550:       PetscCallExternal(ex_get_elem_type, exoid, cs_id[cs], elem_type);
1551:       PetscCall(ExodusGetCellType_Internal(elem_type, &ct));
1552:       dim = PetscMax(dim, DMPolytopeTypeGetDim(ct));
1553:       PetscCallExternal(ex_get_block, exoid, EX_ELEM_BLOCK, cs_id[cs], buffer, &num_cell_in_set, &num_vertex_per_cell, 0, 0, &num_attr);
1554:       switch (ct) {
1555:       case DM_POLYTOPE_TRI_PRISM:
1556:         cs_order[cs] = cs;
1557:         ++num_hybrid;
1558:         break;
1559:       default:
1560:         for (c = cs; c > cs - num_hybrid; --c) cs_order[c] = cs_order[c - 1];
1561:         cs_order[cs - num_hybrid] = cs;
1562:       }
1563:     }
1564:     /* First set sizes */
1565:     for (ncs = 0, c = 0; ncs < num_cs; ++ncs) {
1566:       DMPolytopeType ct;
1567:       char           elem_type[PETSC_MAX_PATH_LEN];
1568:       const PetscInt cs = cs_order[ncs];

1570:       PetscCall(PetscArrayzero(elem_type, sizeof(elem_type)));
1571:       PetscCallExternal(ex_get_elem_type, exoid, cs_id[cs], elem_type);
1572:       PetscCall(ExodusGetCellType_Internal(elem_type, &ct));
1573:       PetscCallExternal(ex_get_block, exoid, EX_ELEM_BLOCK, cs_id[cs], buffer, &num_cell_in_set, &num_vertex_per_cell, 0, 0, &num_attr);
1574:       for (c_loc = 0; c_loc < num_cell_in_set; ++c_loc, ++c) {
1575:         PetscCall(DMPlexSetConeSize(*dm, c, num_vertex_per_cell));
1576:         PetscCall(DMPlexSetCellType(*dm, c, ct));
1577:       }
1578:     }
1579:     for (v = numCells; v < numCells + numVertices; ++v) PetscCall(DMPlexSetCellType(*dm, v, DM_POLYTOPE_POINT));
1580:     PetscCall(DMSetUp(*dm));
1581:     for (ncs = 0, c = 0; ncs < num_cs; ++ncs) {
1582:       const PetscInt cs = cs_order[ncs];
1583:       PetscCallExternal(ex_get_block, exoid, EX_ELEM_BLOCK, cs_id[cs], buffer, &num_cell_in_set, &num_vertex_per_cell, 0, 0, &num_attr);
1584:       PetscCall(PetscMalloc2(num_vertex_per_cell * num_cell_in_set, &cs_connect, num_vertex_per_cell, &cone));
1585:       PetscCallExternal(ex_get_conn, exoid, EX_ELEM_BLOCK, cs_id[cs], cs_connect, NULL, NULL);
1586:       /* EXO uses Fortran-based indexing, DMPlex uses C-style and numbers cell first then vertices. */
1587:       for (c_loc = 0, v = 0; c_loc < num_cell_in_set; ++c_loc, ++c) {
1588:         DMPolytopeType ct;

1590:         for (v_loc = 0; v_loc < num_vertex_per_cell; ++v_loc, ++v) cone[v_loc] = cs_connect[v] + numCells - 1;
1591:         PetscCall(DMPlexGetCellType(*dm, c, &ct));
1592:         PetscCall(DMPlexInvertCell(ct, cone));
1593:         PetscCall(DMPlexSetCone(*dm, c, cone));
1594:         PetscCall(DMSetLabelValue_Fast(*dm, &cellSets, "Cell Sets", c, cs_id[cs]));
1595:       }
1596:       PetscCall(PetscFree2(cs_connect, cone));
1597:     }
1598:     PetscCall(PetscFree2(cs_id, cs_order));
1599:   }
1600:   {
1601:     PetscInt ints[] = {dim, dimEmbed};

1603:     PetscCallMPI(MPI_Bcast(ints, 2, MPIU_INT, 0, comm));
1604:     PetscCall(DMSetDimension(*dm, ints[0]));
1605:     PetscCall(DMSetCoordinateDim(*dm, ints[1]));
1606:     dim      = ints[0];
1607:     dimEmbed = ints[1];
1608:   }
1609:   PetscCall(DMPlexSymmetrize(*dm));
1610:   PetscCall(DMPlexStratify(*dm));
1611:   if (interpolate) {
1612:     DM idm;

1614:     PetscCall(DMPlexInterpolate(*dm, &idm));
1615:     PetscCall(DMDestroy(dm));
1616:     *dm = idm;
1617:   }

1619:   /* Create vertex set label */
1620:   if (rank == 0 && (num_vs > 0)) {
1621:     int vs, v;
1622:     /* Read from ex_get_node_set_ids() */
1623:     int *vs_id;
1624:     /* Read from ex_get_node_set_param() */
1625:     int num_vertex_in_set;
1626:     /* Read from ex_get_node_set() */
1627:     int *vs_vertex_list;

1629:     /* Get vertex set ids */
1630:     PetscCall(PetscMalloc1(num_vs, &vs_id));
1631:     PetscCallExternal(ex_get_ids, exoid, EX_NODE_SET, vs_id);
1632:     for (vs = 0; vs < num_vs; ++vs) {
1633:       PetscCallExternal(ex_get_set_param, exoid, EX_NODE_SET, vs_id[vs], &num_vertex_in_set, NULL);
1634:       PetscCall(PetscMalloc1(num_vertex_in_set, &vs_vertex_list));
1635:       PetscCallExternal(ex_get_set, exoid, EX_NODE_SET, vs_id[vs], vs_vertex_list, NULL);
1636:       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]));
1637:       PetscCall(PetscFree(vs_vertex_list));
1638:     }
1639:     PetscCall(PetscFree(vs_id));
1640:   }
1641:   /* Read coordinates */
1642:   PetscCall(DMGetCoordinateSection(*dm, &coordSection));
1643:   PetscCall(PetscSectionSetNumFields(coordSection, 1));
1644:   PetscCall(PetscSectionSetFieldComponents(coordSection, 0, dimEmbed));
1645:   PetscCall(PetscSectionSetChart(coordSection, numCells, numCells + numVertices));
1646:   for (v = numCells; v < numCells + numVertices; ++v) {
1647:     PetscCall(PetscSectionSetDof(coordSection, v, dimEmbed));
1648:     PetscCall(PetscSectionSetFieldDof(coordSection, v, 0, dimEmbed));
1649:   }
1650:   PetscCall(PetscSectionSetUp(coordSection));
1651:   PetscCall(PetscSectionGetStorageSize(coordSection, &coordSize));
1652:   PetscCall(VecCreate(PETSC_COMM_SELF, &coordinates));
1653:   PetscCall(PetscObjectSetName((PetscObject)coordinates, "coordinates"));
1654:   PetscCall(VecSetSizes(coordinates, coordSize, PETSC_DETERMINE));
1655:   PetscCall(VecSetBlockSize(coordinates, dimEmbed));
1656:   PetscCall(VecSetType(coordinates, VECSTANDARD));
1657:   PetscCall(VecGetArray(coordinates, &coords));
1658:   if (rank == 0) {
1659:     PetscReal *x, *y, *z;

1661:     PetscCall(PetscMalloc3(numVertices, &x, numVertices, &y, numVertices, &z));
1662:     PetscCallExternal(ex_get_coord, exoid, x, y, z);
1663:     if (dimEmbed > 0) {
1664:       for (v = 0; v < numVertices; ++v) coords[v * dimEmbed + 0] = x[v];
1665:     }
1666:     if (dimEmbed > 1) {
1667:       for (v = 0; v < numVertices; ++v) coords[v * dimEmbed + 1] = y[v];
1668:     }
1669:     if (dimEmbed > 2) {
1670:       for (v = 0; v < numVertices; ++v) coords[v * dimEmbed + 2] = z[v];
1671:     }
1672:     PetscCall(PetscFree3(x, y, z));
1673:   }
1674:   PetscCall(VecRestoreArray(coordinates, &coords));
1675:   PetscCall(DMSetCoordinatesLocal(*dm, coordinates));
1676:   PetscCall(VecDestroy(&coordinates));

1678:   /* Create side set label */
1679:   if (rank == 0 && interpolate && (num_fs > 0)) {
1680:     int fs, f, voff;
1681:     /* Read from ex_get_side_set_ids() */
1682:     int *fs_id;
1683:     /* Read from ex_get_side_set_param() */
1684:     int num_side_in_set;
1685:     /* Read from ex_get_side_set_node_list() */
1686:     int *fs_vertex_count_list, *fs_vertex_list, *fs_side_list;
1687:     /* Read side set labels */
1688:     char   fs_name[MAX_STR_LENGTH + 1];
1689:     size_t fs_name_len;

1691:     /* Get side set ids */
1692:     PetscCall(PetscMalloc1(num_fs, &fs_id));
1693:     PetscCallExternal(ex_get_ids, exoid, EX_SIDE_SET, fs_id);
1694:     // Ids 1 and 2 are reserved by ExodusII for indicating things in 3D
1695:     for (fs = 0; fs < num_fs; ++fs) {
1696:       PetscCallExternal(ex_get_set_param, exoid, EX_SIDE_SET, fs_id[fs], &num_side_in_set, NULL);
1697:       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));
1698:       PetscCallExternal(ex_get_side_set_node_list, exoid, fs_id[fs], fs_vertex_count_list, fs_vertex_list);
1699:       PetscCallExternal(ex_get_set, exoid, EX_SIDE_SET, fs_id[fs], NULL, fs_side_list);

1701:       /* Get the specific name associated with this side set ID. */
1702:       int fs_name_err = ex_get_name(exoid, EX_SIDE_SET, fs_id[fs], fs_name);
1703:       if (!fs_name_err) {
1704:         PetscCall(PetscStrlen(fs_name, &fs_name_len));
1705:         if (fs_name_len == 0) PetscCall(PetscStrncpy(fs_name, "Face Sets", MAX_STR_LENGTH));
1706:       }
1707:       for (f = 0, voff = 0; f < num_side_in_set; ++f) {
1708:         const PetscInt *faces    = NULL;
1709:         PetscInt        faceSize = fs_vertex_count_list[f], numFaces;
1710:         PetscInt        faceVertices[4], v;

1712:         PetscCheck(faceSize <= 4, comm, PETSC_ERR_ARG_WRONG, "ExodusII side cannot have %" PetscInt_FMT " > 4 vertices", faceSize);
1713:         for (v = 0; v < faceSize; ++v, ++voff) faceVertices[v] = fs_vertex_list[voff] + numCells - 1;
1714:         PetscCall(DMPlexGetFullJoin(*dm, faceSize, faceVertices, &numFaces, &faces));
1715:         PetscCheck(numFaces == 1, comm, PETSC_ERR_ARG_WRONG, "Invalid ExodusII side %d in set %d maps to %" PetscInt_FMT " faces", f, fs, numFaces);
1716:         PetscCall(DMSetLabelValue_Fast(*dm, &faceSets, "Face Sets", faces[0], fs_id[fs]));
1717:         /* Only add the label if one has been detected for this side set. */
1718:         if (!fs_name_err) PetscCall(DMSetLabelValue(*dm, fs_name, faces[0], fs_id[fs]));
1719:         PetscCall(DMPlexRestoreJoin(*dm, faceSize, faceVertices, &numFaces, &faces));
1720:       }
1721:       PetscCall(PetscFree3(fs_vertex_count_list, fs_vertex_list, fs_side_list));
1722:     }
1723:     PetscCall(PetscFree(fs_id));
1724:   }

1726:   { /* Create Cell/Face/Vertex Sets labels at all processes */
1727:     enum {
1728:       n = 3
1729:     };
1730:     PetscBool flag[n];

1732:     flag[0] = cellSets ? PETSC_TRUE : PETSC_FALSE;
1733:     flag[1] = faceSets ? PETSC_TRUE : PETSC_FALSE;
1734:     flag[2] = vertSets ? PETSC_TRUE : PETSC_FALSE;
1735:     PetscCallMPI(MPI_Bcast(flag, n, MPIU_BOOL, 0, comm));
1736:     if (flag[0]) PetscCall(DMCreateLabel(*dm, "Cell Sets"));
1737:     if (flag[1]) PetscCall(DMCreateLabel(*dm, "Face Sets"));
1738:     if (flag[2]) PetscCall(DMCreateLabel(*dm, "Vertex Sets"));
1739:   }
1740:   PetscFunctionReturn(PETSC_SUCCESS);
1741: #else
1742:   SETERRQ(comm, PETSC_ERR_SUP, "This method requires ExodusII support. Reconfigure using --download-exodusii");
1743: #endif
1744: }