Actual source code: plexexodusii.c

petsc-3.13.6 2020-09-29
Report Typos and Errors
  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: /*
 13:   PETSC_VIEWER_EXODUSII_ - Creates an ExodusII PetscViewer shared by all processors in a communicator.

 15:   Collective

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

 20:   Level: intermediate

 22:   Notes:
 23:     misses Fortran bindings

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

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

 38:   PetscViewerExodusIIOpen(comm, "mesh.exo", FILE_MODE_WRITE, &viewer);
 39:   if (ierr) {PetscError(PETSC_COMM_SELF,__LINE__,"PETSC_VIEWER_EXODUSII_",__FILE__,PETSC_ERR_PLIB,PETSC_ERROR_INITIAL," ");return(0);}
 40:   PetscObjectRegisterDestroy((PetscObject) viewer);
 41:   if (ierr) {PetscError(PETSC_COMM_SELF,__LINE__,"PETSC_VIEWER_EXODUSII_",__FILE__,PETSC_ERR_PLIB,PETSC_ERROR_INITIAL," ");return(0);}
 42:   PetscFunctionReturn(viewer);
 43: }

 45: static PetscErrorCode PetscViewerView_ExodusII(PetscViewer v, PetscViewer viewer)
 46: {
 47:   PetscViewer_ExodusII *exo = (PetscViewer_ExodusII *) viewer->data;
 48:   PetscErrorCode        ierr;

 51:   if (exo->filename) {PetscViewerASCIIPrintf(viewer, "Filename: %s\n", exo->filename);}
 52:   return(0);
 53: }

 55: static PetscErrorCode PetscViewerSetFromOptions_ExodusII(PetscOptionItems *PetscOptionsObject, PetscViewer v)
 56: {

 60:   PetscOptionsHead(PetscOptionsObject, "ExodusII PetscViewer Options");
 61:   PetscOptionsTail();
 62:   return(0);
 63: }

 65: static PetscErrorCode PetscViewerSetUp_ExodusII(PetscViewer viewer)
 66: {
 68:   return(0);
 69: }

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

 77:   if (exo->exoid >= 0) {ex_close(exo->exoid);}
 78:   PetscFree(exo);
 79:   PetscObjectComposeFunction((PetscObject)viewer,"PetscViewerFileSetName_C",NULL);
 80:   PetscObjectComposeFunction((PetscObject)viewer,"PetscViewerFileGetName_C",NULL);
 81:   PetscObjectComposeFunction((PetscObject)viewer,"PetscViewerFileSetMode_C",NULL);
 82:   PetscObjectComposeFunction((PetscObject)viewer,"PetscViewerExodusIIGetId",NULL);
 83:   return(0);
 84: }

 86: static PetscErrorCode PetscViewerFileSetName_ExodusII(PetscViewer viewer, const char name[])
 87: {
 88:   PetscViewer_ExodusII *exo = (PetscViewer_ExodusII *) viewer->data;
 89:   PetscMPIInt           rank;
 90:   int                   CPU_word_size, IO_word_size, EXO_mode;
 91:   PetscErrorCode        ierr;


 94:   MPI_Comm_rank(PetscObjectComm((PetscObject) viewer), &rank);
 95:   CPU_word_size = sizeof(PetscReal);
 96:   IO_word_size  = sizeof(PetscReal);

 99:   if (exo->exoid >= 0) {ex_close(exo->exoid); exo->exoid = -1;}
100:   if (exo->filename) {PetscFree(exo->filename);}
101:   PetscStrallocpy(name, &exo->filename);
102:   /* Create or open the file collectively */
103:   switch (exo->btype) {
104:   case FILE_MODE_READ:
105:     EXO_mode = EX_CLOBBER;
106:     break;
107:   case FILE_MODE_APPEND:
108:     EXO_mode = EX_CLOBBER;
109:     break;
110:   case FILE_MODE_WRITE:
111:     EXO_mode = EX_CLOBBER;
112:     break;
113:   default:
114:     SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER, "Must call PetscViewerFileSetMode() before PetscViewerFileSetName()");
115:   }
116:   #if defined(PETSC_USE_64BIT_INDICES)
117:   EXO_mode += EX_ALL_INT64_API;
118:   #endif
119:   exo->exoid = ex_create(name, EXO_mode, &CPU_word_size, &IO_word_size);
120:   if (exo->exoid < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_LIB, "ex_create failed for %s", name);
121:   return(0);
122: }

124: static PetscErrorCode PetscViewerFileGetName_ExodusII(PetscViewer viewer, const char **name)
125: {
126:   PetscViewer_ExodusII *exo = (PetscViewer_ExodusII *) viewer->data;

129:   *name = exo->filename;
130:   return(0);
131: }

133: static PetscErrorCode PetscViewerFileSetMode_ExodusII(PetscViewer viewer, PetscFileMode type)
134: {
135:   PetscViewer_ExodusII *exo = (PetscViewer_ExodusII *) viewer->data;

138:   exo->btype = type;
139:   return(0);
140: }

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

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

151: static PetscErrorCode PetscViewerExodusIIGetId_ExodusII(PetscViewer viewer, int *exoid)
152: {
153:   PetscViewer_ExodusII *exo = (PetscViewer_ExodusII *) viewer->data;

156:   *exoid = exo->exoid;
157:   return(0);
158: }

160: /*@C
161:    PetscViewerExodusIIOpen - Opens a file for ExodusII input/output.

163:    Collective

165:    Input Parameters:
166: +  comm - MPI communicator
167: .  name - name of file
168: -  type - type of file
169: $    FILE_MODE_WRITE - create new file for binary output
170: $    FILE_MODE_READ - open existing file for binary input
171: $    FILE_MODE_APPEND - open existing file for binary output

173:    Output Parameter:
174: .  exo - PetscViewer for Exodus II input/output to use with the specified file

176:    Level: beginner

178:    Note:
179:    This PetscViewer should be destroyed with PetscViewerDestroy().


182: .seealso: PetscViewerPushFormat(), PetscViewerDestroy(),
183:           DMLoad(), PetscFileMode, PetscViewer, PetscViewerSetType(), PetscViewerFileSetMode(), PetscViewerFileSetName()
184: @*/
185: PetscErrorCode PetscViewerExodusIIOpen(MPI_Comm comm, const char name[], PetscFileMode type, PetscViewer *exo)
186: {

190:   PetscViewerCreate(comm, exo);
191:   PetscViewerSetType(*exo, PETSCVIEWEREXODUSII);
192:   PetscViewerFileSetMode(*exo, type);
193:   PetscViewerFileSetName(*exo, name);
194:   PetscViewerSetFromOptions(*exo);
195:   return(0);
196: }

198: /*MC
199:    PETSCVIEWEREXODUSII - A viewer that writes to an Exodus II file


202: .seealso:  PetscViewerExodusIIOpen(), PetscViewerCreate(), PETSCVIEWERBINARY, PETSCVIEWERHDF5, DMView(),
203:            PetscViewerFileSetName(), PetscViewerFileSetMode(), PetscViewerFormat, PetscViewerType, PetscViewerSetType()

205:   Level: beginner
206: M*/

208: PETSC_EXTERN PetscErrorCode PetscViewerCreate_ExodusII(PetscViewer v)
209: {
210:   PetscViewer_ExodusII *exo;
211:   PetscErrorCode        ierr;

214:   PetscNewLog(v,&exo);

216:   v->data                = (void*) exo;
217:   v->ops->destroy        = PetscViewerDestroy_ExodusII;
218:   v->ops->setfromoptions = PetscViewerSetFromOptions_ExodusII;
219:   v->ops->setup          = PetscViewerSetUp_ExodusII;
220:   v->ops->view           = PetscViewerView_ExodusII;
221:   v->ops->flush          = 0;
222:   exo->btype             = (PetscFileMode) -1;
223:   exo->filename          = 0;
224:   exo->exoid             = -1;

226:   PetscObjectComposeFunction((PetscObject)v,"PetscViewerFileSetName_C",PetscViewerFileSetName_ExodusII);
227:   PetscObjectComposeFunction((PetscObject)v,"PetscViewerFileGetName_C",PetscViewerFileGetName_ExodusII);
228:   PetscObjectComposeFunction((PetscObject)v,"PetscViewerFileSetMode_C",PetscViewerFileSetMode_ExodusII);
229:   PetscObjectComposeFunction((PetscObject)v,"PetscViewerFileGetMode_C",PetscViewerFileGetMode_ExodusII);
230:   PetscObjectComposeFunction((PetscObject)v,"PetscViewerExodusIIGetId_C",PetscViewerExodusIIGetId_ExodusII);
231:   return(0);
232: }

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

237:   Not collective

239:   Input Parameters:
240: + exoid    - the exodus id of a file (obtained from ex_open or ex_create for instance)
241: . obj_type - the type of entity for instance EX_NODAL, EX_ELEM_BLOCK
242: - name     - the name of the result

244:   Output Parameters:
245: . varIndex - the location in the exodus file of the result

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

253:   Level: beginner

255: .seealso: EXOGetVarIndex(),DMPlexView_ExodusII_Internal(),VecViewPlex_ExodusII_Nodal_Internal(),VecLoadNodal_PlexEXO(),VecLoadZonal_PlexEXO()
256: */
257: static PetscErrorCode EXOGetVarIndex_Private(int exoid, ex_entity_type obj_type, const char name[], int *varIndex)
258: {
259:   int            num_vars, i, j;
260:   char           ext_name[MAX_STR_LENGTH+1], var_name[MAX_STR_LENGTH+1];
261:   const int      num_suffix = 5;
262:   char          *suffix[5];
263:   PetscBool      flg;

267:   suffix[0] = (char *) "";
268:   suffix[1] = (char *) "_X";
269:   suffix[2] = (char *) "_XX";
270:   suffix[3] = (char *) "_1";
271:   suffix[4] = (char *) "_11";

273:   *varIndex = 0;
274:   PetscStackCallStandard(ex_get_variable_param,(exoid, obj_type, &num_vars));
275:   for (i = 0; i < num_vars; ++i) {
276:     PetscStackCallStandard(ex_get_variable_name,(exoid, obj_type, i+1, var_name));
277:     for (j = 0; j < num_suffix; ++j){
278:       PetscStrncpy(ext_name, name, MAX_STR_LENGTH);
279:       PetscStrlcat(ext_name, suffix[j], MAX_STR_LENGTH);
280:       PetscStrcasecmp(ext_name, var_name, &flg);
281:       if (flg) {
282:         *varIndex = i+1;
283:         return(0);
284:       }
285:     }
286:   }
287:   SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unable to locate variable %s in ExodusII file.", name);
288:  PetscFunctionReturn(-1);
289: }

291: /*
292:   DMPlexView_ExodusII_Internal - Write a DM to disk in exodus format

294:   Collective on dm

296:   Input Parameters:
297: + dm  - The dm to be written
298: . exoid - the exodus id of a file (obtained from ex_open or ex_create for instance)
299: - degree - the degree of the interpolation space

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

305:   If the dm has been distributed and exoid points to different files on each MPI rank, only the "local" part of the DM
306:   (including "ghost" cells and vertices) will be written. If all exoid points to the same file, the resulting file will
307:   probably be corrupted.

309:   DMPlex only represents geometry while most post-processing software expect that a mesh also provides information
310:   on the discretization space. This function assumes that the file represents Lagrange finite elements of order 1 or 2.
311:   It should be extended to use PetscFE objects.

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

316: .seealso: EXOGetVarIndex_Private(), VecViewPlex_ExodusII_Nodal_Internal(), VecLoadNodal_PlexEXO(), VecViewZonal_PlexEXO(), VecLoadZonal_PlexEXO()
317: */
318: PetscErrorCode DMPlexView_ExodusII_Internal(DM dm, int exoid, PetscInt degree)
319: {
320:   enum ElemType {TRI, QUAD, TET, HEX};
321:   MPI_Comm        comm;
322:   /* Connectivity Variables */
323:   PetscInt        cellsNotInConnectivity;
324:   /* Cell Sets */
325:   DMLabel         csLabel;
326:   IS              csIS;
327:   const PetscInt *csIdx;
328:   PetscInt        num_cs, cs;
329:   enum ElemType  *type;
330:   PetscBool       hasLabel;
331:   /* Coordinate Variables */
332:   DM              cdm;
333:   PetscSection    section;
334:   Vec             coord;
335:   PetscInt      **nodes;
336:   PetscInt        depth, d, dim, skipCells = 0;
337:   PetscInt        pStart, pEnd, p, cStart, cEnd, numCells, vStart, vEnd, numVertices, eStart, eEnd, numEdges, fStart, fEnd, numFaces, numNodes;
338:   PetscInt        num_vs, num_fs;
339:   PetscMPIInt     rank, size;
340:   const char     *dmName;
341:   PetscInt        nodesTriP1[4]  = {3,0,0,0};
342:   PetscInt        nodesTriP2[4]  = {3,3,0,0};
343:   PetscInt        nodesQuadP1[4] = {4,0,0,0};
344:   PetscInt        nodesQuadP2[4] = {4,4,0,1};
345:   PetscInt        nodesTetP1[4]  = {4,0,0,0};
346:   PetscInt        nodesTetP2[4]  = {4,6,0,0};
347:   PetscInt        nodesHexP1[4]  = {8,0,0,0};
348:   PetscInt        nodesHexP2[4]  = {8,12,6,1};
349:   PetscErrorCode  ierr;

352:   PetscObjectGetComm((PetscObject)dm, &comm);
353:   MPI_Comm_rank(comm, &rank);
354:   MPI_Comm_size(comm, &size);
355:   /* --- Get DM info --- */
356:   PetscObjectGetName((PetscObject) dm, &dmName);
357:   DMPlexGetDepth(dm, &depth);
358:   DMGetDimension(dm, &dim);
359:   DMPlexGetChart(dm, &pStart, &pEnd);
360:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
361:   DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);
362:   DMPlexGetDepthStratum(dm, 1, &eStart, &eEnd);
363:   DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
364:   numCells    = cEnd - cStart;
365:   numEdges    = eEnd - eStart;
366:   numVertices = vEnd - vStart;
367:   if (depth == 3) {numFaces = fEnd - fStart;}
368:   else            {numFaces = 0;}
369:   DMGetLabelSize(dm, "Cell Sets", &num_cs);
370:   DMGetLabelSize(dm, "Vertex Sets", &num_vs);
371:   DMGetLabelSize(dm, "Face Sets", &num_fs);
372:   DMGetCoordinatesLocal(dm, &coord);
373:   DMGetCoordinateDM(dm, &cdm);
374:   if (num_cs > 0) {
375:     DMGetLabel(dm, "Cell Sets", &csLabel);
376:     DMLabelGetValueIS(csLabel, &csIS);
377:     ISGetIndices(csIS, &csIdx);
378:   }
379:   PetscMalloc1(num_cs, &nodes);
380:   /* Set element type for each block and compute total number of nodes */
381:   PetscMalloc1(num_cs, &type);
382:   numNodes = numVertices;
383:   if (degree == 2) {numNodes += numEdges;}
384:   cellsNotInConnectivity = numCells;
385:   for (cs = 0; cs < num_cs; ++cs) {
386:     IS              stratumIS;
387:     const PetscInt *cells;
388:     PetscScalar    *xyz = NULL;
389:     PetscInt        csSize, closureSize;

391:     DMLabelGetStratumIS(csLabel, csIdx[cs], &stratumIS);
392:     ISGetIndices(stratumIS, &cells);
393:     ISGetSize(stratumIS, &csSize);
394:     DMPlexVecGetClosure(cdm, NULL, coord, cells[0], &closureSize, &xyz);
395:     switch (dim) {
396:       case 2:
397:         if      (closureSize == 3*dim) {type[cs] = TRI;}
398:         else if (closureSize == 4*dim) {type[cs] = QUAD;}
399:         else SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Number of vertices %D in dimension %D has no ExodusII type", closureSize/dim, dim);
400:         break;
401:       case 3:
402:         if      (closureSize == 4*dim) {type[cs] = TET;}
403:         else if (closureSize == 8*dim) {type[cs] = HEX;}
404:         else SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Number of vertices %D in dimension %D has no ExodusII type", closureSize/dim, dim);
405:         break;
406:       default: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Dimension %D not handled by ExodusII viewer", dim);
407:     }
408:     if ((degree == 2) && (type[cs] == QUAD)) {numNodes += csSize;}
409:     if ((degree == 2) && (type[cs] == HEX))  {numNodes += csSize; numNodes += numFaces;}
410:     DMPlexVecRestoreClosure(cdm, NULL, coord, cells[0], &closureSize, &xyz);
411:     /* Set nodes and Element type */
412:     if (type[cs] == TRI) {
413:       if      (degree == 1) nodes[cs] = nodesTriP1;
414:       else if (degree == 2) nodes[cs] = nodesTriP2;
415:     } else if (type[cs] == QUAD) {
416:       if      (degree == 1) nodes[cs] = nodesQuadP1;
417:       else if (degree == 2) nodes[cs] = nodesQuadP2;
418:     } else if (type[cs] == TET) {
419:       if      (degree == 1) nodes[cs] = nodesTetP1;
420:       else if (degree == 2) nodes[cs] = nodesTetP2;
421:     } else if (type[cs] == HEX) {
422:       if      (degree == 1) nodes[cs] = nodesHexP1;
423:       else if (degree == 2) nodes[cs] = nodesHexP2;
424:     }
425:     /* Compute the number of cells not in the connectivity table */
426:     cellsNotInConnectivity -= nodes[cs][3]*csSize;

428:     ISRestoreIndices(stratumIS, &cells);
429:     ISDestroy(&stratumIS);
430:   }
431:   if (num_cs > 0) {PetscStackCallStandard(ex_put_init,(exoid, dmName, dim, numNodes, numCells, num_cs, num_vs, num_fs));}
432:   /* --- Connectivity --- */
433:   for (cs = 0; cs < num_cs; ++cs) {
434:     IS              stratumIS;
435:     const PetscInt *cells;
436:     PetscInt       *connect, off = 0;
437:     PetscInt        edgesInClosure = 0, facesInClosure = 0, verticesInClosure = 0;
438:     PetscInt        csSize, c, connectSize, closureSize;
439:     char           *elem_type = NULL;
440:     char            elem_type_tri3[]  = "TRI3",  elem_type_quad4[] = "QUAD4";
441:     char            elem_type_tri6[]  = "TRI6",  elem_type_quad9[] = "QUAD9";
442:     char            elem_type_tet4[]  = "TET4",  elem_type_hex8[]  = "HEX8";
443:     char            elem_type_tet10[] = "TET10", elem_type_hex27[] = "HEX27";

445:     DMLabelGetStratumIS(csLabel, csIdx[cs], &stratumIS);
446:     ISGetIndices(stratumIS, &cells);
447:     ISGetSize(stratumIS, &csSize);
448:     /* Set Element type */
449:     if (type[cs] == TRI) {
450:       if      (degree == 1) elem_type = elem_type_tri3;
451:       else if (degree == 2) elem_type = elem_type_tri6;
452:     } else if (type[cs] == QUAD) {
453:       if      (degree == 1) elem_type = elem_type_quad4;
454:       else if (degree == 2) elem_type = elem_type_quad9;
455:     } else if (type[cs] == TET) {
456:       if      (degree == 1) elem_type = elem_type_tet4;
457:       else if (degree == 2) elem_type = elem_type_tet10;
458:     } else if (type[cs] == HEX) {
459:       if      (degree == 1) elem_type = elem_type_hex8;
460:       else if (degree == 2) elem_type = elem_type_hex27;
461:     }
462:     connectSize = nodes[cs][0] + nodes[cs][1] + nodes[cs][2] + nodes[cs][3];
463:     PetscMalloc1(PetscMax(27,connectSize)*csSize, &connect);
464:     PetscStackCallStandard(ex_put_block,(exoid, EX_ELEM_BLOCK, csIdx[cs], elem_type, csSize, connectSize, 0, 0, 1));
465:     /* Find number of vertices, edges, and faces in the closure */
466:     verticesInClosure = nodes[cs][0];
467:     if (depth > 1) {
468:       if (dim == 2) {
469:         DMPlexGetConeSize(dm, cells[0], &edgesInClosure);
470:       } else if (dim == 3) {
471:         PetscInt *closure = NULL;

473:         DMPlexGetConeSize(dm, cells[0], &facesInClosure);
474:         DMPlexGetTransitiveClosure(dm, cells[0], PETSC_TRUE, &closureSize, &closure);
475:         edgesInClosure = closureSize - facesInClosure - 1 - verticesInClosure;
476:         DMPlexRestoreTransitiveClosure(dm, cells[0], PETSC_TRUE, &closureSize, &closure);
477:       }
478:     }
479:     /* Get connectivity for each cell */
480:     for (c = 0; c < csSize; ++c) {
481:       PetscInt *closure = NULL;
482:       PetscInt  temp, i;

484:       DMPlexGetTransitiveClosure(dm, cells[c], PETSC_TRUE, &closureSize, &closure);
485:       for (i = 0; i < connectSize; ++i) {
486:         if (i < nodes[cs][0]) {/* Vertices */
487:           connect[i+off] = closure[(i+edgesInClosure+facesInClosure+1)*2] + 1;
488:           connect[i+off] -= cellsNotInConnectivity;
489:         } else if (i < nodes[cs][0]+nodes[cs][1]) { /* Edges */
490:           connect[i+off] = closure[(i-verticesInClosure+facesInClosure+1)*2] + 1;
491:           if (nodes[cs][2] == 0) connect[i+off] -= numFaces;
492:           connect[i+off] -= cellsNotInConnectivity;
493:         } else if (i < nodes[cs][0]+nodes[cs][1]+nodes[cs][3]){ /* Cells */
494:           connect[i+off] = closure[0] + 1;
495:           connect[i+off] -= skipCells;
496:         } else if (i < nodes[cs][0]+nodes[cs][1]+nodes[cs][3]+nodes[cs][2]){ /* Faces */
497:           connect[i+off] = closure[(i-edgesInClosure-verticesInClosure)*2] + 1;
498:           connect[i+off] -= cellsNotInConnectivity;
499:         } else {
500:           connect[i+off] = -1;
501:         }
502:       }
503:       /* Tetrahedra are inverted */
504:       if (type[cs] == TET) {
505:         temp = connect[0+off]; connect[0+off] = connect[1+off]; connect[1+off] = temp;
506:         if (degree == 2) {
507:           temp = connect[5+off]; connect[5+off] = connect[6+off]; connect[6+off] = temp;
508:           temp = connect[7+off]; connect[7+off] = connect[8+off]; connect[8+off] = temp;
509:         }
510:       }
511:       /* Hexahedra are inverted */
512:       if (type[cs] == HEX) {
513:         temp = connect[1+off]; connect[1+off] = connect[3+off]; connect[3+off] = temp;
514:         if (degree == 2) {
515:           temp = connect[8+off];  connect[8+off]  = connect[11+off]; connect[11+off] = temp;
516:           temp = connect[9+off];  connect[9+off]  = connect[10+off]; connect[10+off] = temp;
517:           temp = connect[16+off]; connect[16+off] = connect[17+off]; connect[17+off] = temp;
518:           temp = connect[18+off]; connect[18+off] = connect[19+off]; connect[19+off] = temp;

520:           temp = connect[12+off]; connect[12+off] = connect[16+off]; connect[16+off] = temp;
521:           temp = connect[13+off]; connect[13+off] = connect[17+off]; connect[17+off] = temp;
522:           temp = connect[14+off]; connect[14+off] = connect[18+off]; connect[18+off] = temp;
523:           temp = connect[15+off]; connect[15+off] = connect[19+off]; connect[19+off] = temp;

525:           temp = connect[23+off]; connect[23+off] = connect[26+off]; connect[26+off] = temp;
526:           temp = connect[24+off]; connect[24+off] = connect[25+off]; connect[25+off] = temp;
527:           temp = connect[25+off]; connect[25+off] = connect[26+off]; connect[26+off] = temp;
528:         }
529:       }
530:       off += connectSize;
531:       DMPlexRestoreTransitiveClosure(dm, cells[c], PETSC_TRUE, &closureSize, &closure);
532:     }
533:     PetscStackCallStandard(ex_put_conn,(exoid, EX_ELEM_BLOCK, csIdx[cs], connect, 0, 0));
534:     skipCells += (nodes[cs][3] == 0)*csSize;
535:     PetscFree(connect);
536:     ISRestoreIndices(stratumIS, &cells);
537:     ISDestroy(&stratumIS);
538:   }
539:   PetscFree(type);
540:   /* --- Coordinates --- */
541:   PetscSectionCreate(comm, &section);
542:   PetscSectionSetChart(section, pStart, pEnd);
543:   if (num_cs) {
544:     for (d = 0; d < depth; ++d) {
545:       DMPlexGetDepthStratum(dm, d, &pStart, &pEnd);
546:       for (p = pStart; p < pEnd; ++p) {
547:         PetscSectionSetDof(section, p, nodes[0][d] > 0);
548:       }
549:     }
550:   }
551:   for (cs = 0; cs < num_cs; ++cs) {
552:     IS              stratumIS;
553:     const PetscInt *cells;
554:     PetscInt        csSize, c;

556:     DMLabelGetStratumIS(csLabel, csIdx[cs], &stratumIS);
557:     ISGetIndices(stratumIS, &cells);
558:     ISGetSize(stratumIS, &csSize);
559:     for (c = 0; c < csSize; ++c) {
560:       PetscSectionSetDof(section, cells[c], nodes[cs][3] > 0);
561:     }
562:     ISRestoreIndices(stratumIS, &cells);
563:     ISDestroy(&stratumIS);
564:   }
565:   if (num_cs > 0) {
566:     ISRestoreIndices(csIS, &csIdx);
567:     ISDestroy(&csIS);
568:   }
569:   PetscFree(nodes);
570:   PetscSectionSetUp(section);
571:   if (numNodes > 0) {
572:     const char  *coordNames[3] = {"x", "y", "z"};
573:     PetscScalar *coords, *closure;
574:     PetscReal   *cval;
575:     PetscInt     hasDof, n = 0;

577:     /* There can't be more than 24 values in the closure of a point for the coord section */
578:     PetscCalloc3(numNodes*3, &coords, dim, &cval, 24, &closure);
579:     DMGetCoordinatesLocal(dm, &coord);
580:     DMPlexGetChart(dm, &pStart, &pEnd);
581:     for (p = pStart; p < pEnd; ++p) {
582:       PetscSectionGetDof(section, p, &hasDof);
583:       if (hasDof) {
584:         PetscInt closureSize = 24, j;

586:         DMPlexVecGetClosure(cdm, NULL, coord, p, &closureSize, &closure);
587:         for (d = 0; d < dim; ++d) {
588:           cval[d] = 0.0;
589:           for (j = 0; j < closureSize/dim; j++) cval[d] += closure[j*dim+d];
590:           coords[d*numNodes+n] = cval[d] * dim / closureSize;
591:         }
592:         ++n;
593:       }
594:     }
595:     PetscStackCallStandard(ex_put_coord,(exoid, &coords[0*numNodes], &coords[1*numNodes], &coords[2*numNodes]));
596:     PetscFree3(coords, cval, closure);
597:     PetscStackCallStandard(ex_put_coord_names,(exoid, (char **) coordNames));
598:   }
599:   PetscSectionDestroy(&section);
600:   /* --- Node Sets/Vertex Sets --- */
601:   DMHasLabel(dm, "Vertex Sets", &hasLabel);
602:   if (hasLabel) {
603:     PetscInt        i, vs, vsSize;
604:     const PetscInt *vsIdx, *vertices;
605:     PetscInt       *nodeList;
606:     IS              vsIS, stratumIS;
607:     DMLabel         vsLabel;
608:     DMGetLabel(dm, "Vertex Sets", &vsLabel);
609:     DMLabelGetValueIS(vsLabel, &vsIS);
610:     ISGetIndices(vsIS, &vsIdx);
611:     for (vs=0; vs<num_vs; ++vs) {
612:       DMLabelGetStratumIS(vsLabel, vsIdx[vs], &stratumIS);
613:       ISGetIndices(stratumIS, &vertices);
614:       ISGetSize(stratumIS, &vsSize);
615:       PetscMalloc1(vsSize, &nodeList);
616:       for (i=0; i<vsSize; ++i) {
617:         nodeList[i] = vertices[i] - skipCells + 1;
618:       }
619:       PetscStackCallStandard(ex_put_set_param,(exoid, EX_NODE_SET, vsIdx[vs], vsSize, 0));
620:       PetscStackCallStandard(ex_put_set,(exoid, EX_NODE_SET, vsIdx[vs], nodeList, NULL));
621:       ISRestoreIndices(stratumIS, &vertices);
622:       ISDestroy(&stratumIS);
623:       PetscFree(nodeList);
624:     }
625:     ISRestoreIndices(vsIS, &vsIdx);
626:     ISDestroy(&vsIS);
627:   }
628:   /* --- Side Sets/Face Sets --- */
629:   DMHasLabel(dm, "Face Sets", &hasLabel);
630:   if (hasLabel) {
631:     PetscInt        i, j, fs, fsSize;
632:     const PetscInt *fsIdx, *faces;
633:     IS              fsIS, stratumIS;
634:     DMLabel         fsLabel;
635:     PetscInt        numPoints, *points;
636:     PetscInt        elem_list_size = 0;
637:     PetscInt       *elem_list, *elem_ind, *side_list;

639:     DMGetLabel(dm, "Face Sets", &fsLabel);
640:     /* Compute size of Node List and Element List */
641:     DMLabelGetValueIS(fsLabel, &fsIS);
642:     ISGetIndices(fsIS, &fsIdx);
643:     for (fs=0; fs<num_fs; ++fs) {
644:       DMLabelGetStratumIS(fsLabel, fsIdx[fs], &stratumIS);
645:       ISGetSize(stratumIS, &fsSize);
646:       elem_list_size += fsSize;
647:       ISDestroy(&stratumIS);
648:     }
649:     PetscMalloc3(num_fs, &elem_ind,
650:                         elem_list_size, &elem_list,
651:                         elem_list_size, &side_list);
652:     elem_ind[0] = 0;
653:     for (fs=0; fs<num_fs; ++fs) {
654:       DMLabelGetStratumIS(fsLabel, fsIdx[fs], &stratumIS);
655:       ISGetIndices(stratumIS, &faces);
656:       ISGetSize(stratumIS, &fsSize);
657:       /* Set Parameters */
658:       PetscStackCallStandard(ex_put_set_param,(exoid, EX_SIDE_SET, fsIdx[fs], fsSize, 0));
659:       /* Indices */
660:       if (fs<num_fs-1) {
661:         elem_ind[fs+1] = elem_ind[fs] + fsSize;
662:       }

664:       for (i=0; i<fsSize; ++i) {
665:         /* Element List */
666:         points = NULL;
667:         DMPlexGetTransitiveClosure(dm, faces[i], PETSC_FALSE, &numPoints, &points);
668:         elem_list[elem_ind[fs] + i] = points[2] +1;
669:         DMPlexRestoreTransitiveClosure(dm, faces[i], PETSC_FALSE, &numPoints, &points);

671:         /* Side List */
672:         points = NULL;
673:         DMPlexGetTransitiveClosure(dm, elem_list[elem_ind[fs] + i]-1, PETSC_TRUE, &numPoints, &points);
674:         for (j=1; j<numPoints; ++j) {
675:           if (points[j*2]==faces[i]) {break;}
676:         }
677:         /* Convert HEX sides */
678:         if (numPoints == 27) {
679:           if      (j == 1) {j = 5;}
680:           else if (j == 2) {j = 6;}
681:           else if (j == 3) {j = 1;}
682:           else if (j == 4) {j = 3;}
683:           else if (j == 5) {j = 2;}
684:           else if (j == 6) {j = 4;}
685:         }
686:         /* Convert TET sides */
687:         if (numPoints == 15) {
688:           --j;
689:           if (j == 0) {j = 4;}
690:         }
691:         side_list[elem_ind[fs] + i] = j;
692:         DMPlexRestoreTransitiveClosure(dm, elem_list[elem_ind[fs] + i]-1, PETSC_TRUE, &numPoints, &points);

694:       }
695:       ISRestoreIndices(stratumIS, &faces);
696:       ISDestroy(&stratumIS);
697:     }
698:     ISRestoreIndices(fsIS, &fsIdx);
699:     ISDestroy(&fsIS);

701:     /* Put side sets */
702:     for (fs=0; fs<num_fs; ++fs) {
703:       PetscStackCallStandard(ex_put_set,(exoid, EX_SIDE_SET, fsIdx[fs], &elem_list[elem_ind[fs]], &side_list[elem_ind[fs]]));
704:     }
705:     PetscFree3(elem_ind,elem_list,side_list);
706:   }
707:   return(0);
708: }

710: /*
711:   VecViewPlex_ExodusII_Nodal_Internal - Write a Vec corresponding to a nodal field to an exodus file

713:   Collective on v

715:   Input Parameters:
716: + v  - The vector to be written
717: . exoid - the exodus id of a file (obtained from ex_open or ex_create for instance)
718: - step - the time step to write at (exodus steps are numbered starting from 1)

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

726:   Level: beginner

728: .seealso: EXOGetVarIndex_Private(),DMPlexView_ExodusII_Internal(),VecLoadNodal_PlexEXO(),VecViewZonal_PlexEXO(),VecLoadZonal_PlexEXO()
729: @*/
730: PetscErrorCode VecViewPlex_ExodusII_Nodal_Internal(Vec v, int exoid, int step)
731: {
732:   MPI_Comm         comm;
733:   PetscMPIInt      size;
734:   DM               dm;
735:   Vec              vNatural, vComp;
736:   const PetscScalar *varray;
737:   const char      *vecname;
738:   PetscInt         xs, xe, bs;
739:   PetscBool        useNatural;
740:   int              offset;
741:   PetscErrorCode   ierr;

744:   PetscObjectGetComm((PetscObject) v, &comm);
745:   MPI_Comm_size(comm, &size);
746:   VecGetDM(v, &dm);
747:   DMGetUseNatural(dm, &useNatural);
748:   useNatural = useNatural && size > 1 ? PETSC_TRUE : PETSC_FALSE;
749:   if (useNatural) {
750:     DMGetGlobalVector(dm, &vNatural);
751:     DMPlexGlobalToNaturalBegin(dm, v, vNatural);
752:     DMPlexGlobalToNaturalEnd(dm, v, vNatural);
753:   } else {
754:     vNatural = v;
755:   }
756:   /* Get the location of the variable in the exodus file */
757:   PetscObjectGetName((PetscObject) v, &vecname);
758:   EXOGetVarIndex_Private(exoid, EX_NODAL, vecname, &offset);
759:   if (!offset) SETERRQ1(comm, PETSC_ERR_FILE_UNEXPECTED, "Cannot locate nodal variable %s in exodus file\n", vecname);
760:   /* Write local chunk of the result in the exodus file
761:      exodus stores each component of a vector-valued field as a separate variable.
762:      We assume that they are stored sequentially */
763:   VecGetOwnershipRange(vNatural, &xs, &xe);
764:   VecGetBlockSize(vNatural, &bs);
765:   if (bs == 1) {
766:     VecGetArrayRead(vNatural, &varray);
767:     PetscStackCallStandard(ex_put_partial_var,(exoid, step, EX_NODAL, offset, 1, xs+1, xe-xs, varray));
768:     VecRestoreArrayRead(vNatural, &varray);
769:   } else {
770:     IS       compIS;
771:     PetscInt c;

773:     ISCreateStride(PETSC_COMM_SELF, (xe-xs)/bs, xs, bs, &compIS);
774:     for (c = 0; c < bs; ++c) {
775:       ISStrideSetStride(compIS, (xe-xs)/bs, xs+c, bs);
776:       VecGetSubVector(vNatural, compIS, &vComp);
777:       VecGetArrayRead(vComp, &varray);
778:       PetscStackCallStandard(ex_put_partial_var,(exoid, step, EX_NODAL, offset+c, 1, xs/bs+1, (xe-xs)/bs, varray));
779:       VecRestoreArrayRead(vComp, &varray);
780:       VecRestoreSubVector(vNatural, compIS, &vComp);
781:     }
782:     ISDestroy(&compIS);
783:   }
784:   if (useNatural) {DMRestoreGlobalVector(dm, &vNatural);}
785:   return(0);
786: }

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

791:   Collective on v

793:   Input Parameters:
794: + v  - The vector to be written
795: . exoid - the exodus id of a file (obtained from ex_open or ex_create for instance)
796: - step - the time step to read at (exodus steps are numbered starting from 1)

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

804:   Level: beginner

806: .seealso: EXOGetVarIndex_Private(), DMPlexView_ExodusII_Internal(), VecViewPlex_ExodusII_Nodal_Internal(), VecViewZonal_PlexEXO(), VecLoadZonal_PlexEXO()
807: */
808: PetscErrorCode VecLoadPlex_ExodusII_Nodal_Internal(Vec v, int exoid, int step)
809: {
810:   MPI_Comm       comm;
811:   PetscMPIInt    size;
812:   DM             dm;
813:   Vec            vNatural, vComp;
814:   PetscScalar   *varray;
815:   const char    *vecname;
816:   PetscInt       xs, xe, bs;
817:   PetscBool      useNatural;
818:   int            offset;

822:   PetscObjectGetComm((PetscObject) v, &comm);
823:   MPI_Comm_size(comm, &size);
824:   VecGetDM(v,&dm);
825:   DMGetUseNatural(dm, &useNatural);
826:   useNatural = useNatural && size > 1 ? PETSC_TRUE : PETSC_FALSE;
827:   if (useNatural) {DMGetGlobalVector(dm,&vNatural);}
828:   else            {vNatural = v;}
829:   /* Get the location of the variable in the exodus file */
830:   PetscObjectGetName((PetscObject) v, &vecname);
831:   EXOGetVarIndex_Private(exoid, EX_NODAL, vecname, &offset);
832:   if (!offset) SETERRQ1(comm, PETSC_ERR_FILE_UNEXPECTED, "Cannot locate nodal variable %s in exodus file", vecname);
833:   /* Read local chunk from the file */
834:   VecGetOwnershipRange(vNatural, &xs, &xe);
835:   VecGetBlockSize(vNatural, &bs);
836:   if (bs == 1) {
837:     VecGetArray(vNatural, &varray);
838:     PetscStackCallStandard(ex_get_partial_var,(exoid, step, EX_NODAL, offset, 1, xs+1, xe-xs, varray));
839:     VecRestoreArray(vNatural, &varray);
840:   } else {
841:     IS       compIS;
842:     PetscInt c;

844:     ISCreateStride(PETSC_COMM_SELF, (xe-xs)/bs, xs, bs, &compIS);
845:     for (c = 0; c < bs; ++c) {
846:       ISStrideSetStride(compIS, (xe-xs)/bs, xs+c, bs);
847:       VecGetSubVector(vNatural, compIS, &vComp);
848:       VecGetArray(vComp, &varray);
849:       PetscStackCallStandard(ex_get_partial_var,(exoid, step, EX_NODAL, offset+c, 1, xs/bs+1, (xe-xs)/bs, varray));
850:       VecRestoreArray(vComp, &varray);
851:       VecRestoreSubVector(vNatural, compIS, &vComp);
852:     }
853:     ISDestroy(&compIS);
854:   }
855:   if (useNatural) {
856:     DMPlexNaturalToGlobalBegin(dm, vNatural, v);
857:     DMPlexNaturalToGlobalEnd(dm, vNatural, v);
858:     DMRestoreGlobalVector(dm, &vNatural);
859:   }
860:   return(0);
861: }

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

866:   Collective on v

868:   Input Parameters:
869: + v  - The vector to be written
870: . exoid - the exodus id of a file (obtained from ex_open or ex_create for instance)
871: - step - the time step to write at (exodus steps are numbered starting from 1)

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

879:   Level: beginner

881: .seealso: EXOGetVarIndex_Private(),DMPlexView_ExodusII_Internal(),VecViewPlex_ExodusII_Nodal_Internal(),VecLoadPlex_ExodusII_Nodal_Internal(),VecLoadPlex_ExodusII_Zonal_Internal()
882: */
883: PetscErrorCode VecViewPlex_ExodusII_Zonal_Internal(Vec v, int exoid, int step)
884: {
885:   MPI_Comm          comm;
886:   PetscMPIInt       size;
887:   DM                dm;
888:   Vec               vNatural, vComp;
889:   const PetscScalar *varray;
890:   const char       *vecname;
891:   PetscInt          xs, xe, bs;
892:   PetscBool         useNatural;
893:   int               offset;
894:   IS                compIS;
895:   PetscInt         *csSize, *csID;
896:   PetscInt          numCS, set, csxs = 0;
897:   PetscErrorCode    ierr;

900:   PetscObjectGetComm((PetscObject)v, &comm);
901:   MPI_Comm_size(comm, &size);
902:   VecGetDM(v, &dm);
903:   DMGetUseNatural(dm, &useNatural);
904:   useNatural = useNatural && size > 1 ? PETSC_TRUE : PETSC_FALSE;
905:   if (useNatural) {
906:     DMGetGlobalVector(dm, &vNatural);
907:     DMPlexGlobalToNaturalBegin(dm, v, vNatural);
908:     DMPlexGlobalToNaturalEnd(dm, v, vNatural);
909:   } else {
910:     vNatural = v;
911:   }
912:   /* Get the location of the variable in the exodus file */
913:   PetscObjectGetName((PetscObject) v, &vecname);
914:   EXOGetVarIndex_Private(exoid, EX_ELEM_BLOCK, vecname, &offset);
915:   if (!offset) SETERRQ1(comm, PETSC_ERR_FILE_UNEXPECTED, "Cannot locate zonal variable %s in exodus file\n", vecname);
916:   /* Write local chunk of the result in the exodus file
917:      exodus stores each component of a vector-valued field as a separate variable.
918:      We assume that they are stored sequentially
919:      Zonal variables are accessed one element block at a time, so we loop through the cell sets,
920:      but once the vector has been reordered to natural size, we cannot use the label informations
921:      to figure out what to save where. */
922:   numCS = ex_inquire_int(exoid, EX_INQ_ELEM_BLK);
923:   PetscMalloc2(numCS, &csID, numCS, &csSize);
924:   PetscStackCallStandard(ex_get_ids,(exoid, EX_ELEM_BLOCK, csID));
925:   for (set = 0; set < numCS; ++set) {
926:     ex_block block;

928:     block.id   = csID[set];
929:     block.type = EX_ELEM_BLOCK;
930:     PetscStackCallStandard(ex_get_block_param,(exoid, &block));
931:     csSize[set] = block.num_entry;
932:   }
933:   VecGetOwnershipRange(vNatural, &xs, &xe);
934:   VecGetBlockSize(vNatural, &bs);
935:   if (bs > 1) {ISCreateStride(comm, (xe-xs)/bs, xs, bs, &compIS);}
936:   for (set = 0; set < numCS; set++) {
937:     PetscInt csLocalSize, c;

939:     /* range of indices for set setID[set]: csxs:csxs + csSize[set]-1
940:        local slice of zonal values:         xs/bs,xm/bs-1
941:        intersection:                        max(xs/bs,csxs),min(xm/bs-1,csxs + csSize[set]-1) */
942:     csLocalSize = PetscMax(0, PetscMin(xe/bs, csxs+csSize[set]) - PetscMax(xs/bs, csxs));
943:     if (bs == 1) {
944:       VecGetArrayRead(vNatural, &varray);
945:       PetscStackCallStandard(ex_put_partial_var,(exoid, step, EX_ELEM_BLOCK, offset, csID[set], PetscMax(xs-csxs, 0)+1, csLocalSize, &varray[PetscMax(0, csxs-xs)]));
946:       VecRestoreArrayRead(vNatural, &varray);
947:     } else {
948:       for (c = 0; c < bs; ++c) {
949:         ISStrideSetStride(compIS, (xe-xs)/bs, xs+c, bs);
950:         VecGetSubVector(vNatural, compIS, &vComp);
951:         VecGetArrayRead(vComp, &varray);
952:         PetscStackCallStandard(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)]));
953:         VecRestoreArrayRead(vComp, &varray);
954:         VecRestoreSubVector(vNatural, compIS, &vComp);
955:       }
956:     }
957:     csxs += csSize[set];
958:   }
959:   PetscFree2(csID, csSize);
960:   if (bs > 1) {ISDestroy(&compIS);}
961:   if (useNatural) {DMRestoreGlobalVector(dm,&vNatural);}
962:   return(0);
963: }

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

968:   Collective on v

970:   Input Parameters:
971: + v  - The vector to be written
972: . exoid - the exodus id of a file (obtained from ex_open or ex_create for instance)
973: - step - the time step to read at (exodus steps are numbered starting from 1)

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

981:   Level: beginner

983: .seealso: EXOGetVarIndex_Private(), DMPlexView_ExodusII_Internal(), VecViewPlex_ExodusII_Nodal_Internal(), VecLoadPlex_ExodusII_Nodal_Internal(), VecLoadPlex_ExodusII_Zonal_Internal()
984: */
985: PetscErrorCode VecLoadPlex_ExodusII_Zonal_Internal(Vec v, int exoid, int step)
986: {
987:   MPI_Comm          comm;
988:   PetscMPIInt       size;
989:   DM                dm;
990:   Vec               vNatural, vComp;
991:   PetscScalar      *varray;
992:   const char       *vecname;
993:   PetscInt          xs, xe, bs;
994:   PetscBool         useNatural;
995:   int               offset;
996:   IS                compIS;
997:   PetscInt         *csSize, *csID;
998:   PetscInt          numCS, set, csxs = 0;
999:   PetscErrorCode    ierr;

1002:   PetscObjectGetComm((PetscObject)v,&comm);
1003:   MPI_Comm_size(comm, &size);
1004:   VecGetDM(v, &dm);
1005:   DMGetUseNatural(dm, &useNatural);
1006:   useNatural = useNatural && size > 1 ? PETSC_TRUE : PETSC_FALSE;
1007:   if (useNatural) {DMGetGlobalVector(dm,&vNatural);}
1008:   else            {vNatural = v;}
1009:   /* Get the location of the variable in the exodus file */
1010:   PetscObjectGetName((PetscObject) v, &vecname);
1011:   EXOGetVarIndex_Private(exoid, EX_ELEM_BLOCK, vecname, &offset);
1012:   if (!offset) SETERRQ1(comm, PETSC_ERR_FILE_UNEXPECTED, "Cannot locate zonal variable %s in exodus file\n", vecname);
1013:   /* Read local chunk of the result in the exodus file
1014:      exodus stores each component of a vector-valued field as a separate variable.
1015:      We assume that they are stored sequentially
1016:      Zonal variables are accessed one element block at a time, so we loop through the cell sets,
1017:      but once the vector has been reordered to natural size, we cannot use the label informations
1018:      to figure out what to save where. */
1019:   numCS = ex_inquire_int(exoid, EX_INQ_ELEM_BLK);
1020:   PetscMalloc2(numCS, &csID, numCS, &csSize);
1021:   PetscStackCallStandard(ex_get_ids,(exoid, EX_ELEM_BLOCK, csID));
1022:   for (set = 0; set < numCS; ++set) {
1023:     ex_block block;

1025:     block.id   = csID[set];
1026:     block.type = EX_ELEM_BLOCK;
1027:     PetscStackCallStandard(ex_get_block_param,(exoid, &block));
1028:     csSize[set] = block.num_entry;
1029:   }
1030:   VecGetOwnershipRange(vNatural, &xs, &xe);
1031:   VecGetBlockSize(vNatural, &bs);
1032:   if (bs > 1) {ISCreateStride(comm, (xe-xs)/bs, xs, bs, &compIS);}
1033:   for (set = 0; set < numCS; ++set) {
1034:     PetscInt csLocalSize, c;

1036:     /* range of indices for set setID[set]: csxs:csxs + csSize[set]-1
1037:        local slice of zonal values:         xs/bs,xm/bs-1
1038:        intersection:                        max(xs/bs,csxs),min(xm/bs-1,csxs + csSize[set]-1) */
1039:     csLocalSize = PetscMax(0, PetscMin(xe/bs, csxs+csSize[set]) - PetscMax(xs/bs, csxs));
1040:     if (bs == 1) {
1041:       VecGetArray(vNatural, &varray);
1042:       PetscStackCallStandard(ex_get_partial_var,(exoid, step, EX_ELEM_BLOCK, offset, csID[set], PetscMax(xs-csxs, 0)+1, csLocalSize, &varray[PetscMax(0, csxs-xs)]));
1043:       VecRestoreArray(vNatural, &varray);
1044:     } else {
1045:       for (c = 0; c < bs; ++c) {
1046:         ISStrideSetStride(compIS, (xe-xs)/bs, xs+c, bs);
1047:         VecGetSubVector(vNatural, compIS, &vComp);
1048:         VecGetArray(vComp, &varray);
1049:         PetscStackCallStandard(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)]));
1050:         VecRestoreArray(vComp, &varray);
1051:         VecRestoreSubVector(vNatural, compIS, &vComp);
1052:       }
1053:     }
1054:     csxs += csSize[set];
1055:   }
1056:   PetscFree2(csID, csSize);
1057:   if (bs > 1) {ISDestroy(&compIS);}
1058:   if (useNatural) {
1059:     DMPlexNaturalToGlobalBegin(dm, vNatural, v);
1060:     DMPlexNaturalToGlobalEnd(dm, vNatural, v);
1061:     DMRestoreGlobalVector(dm, &vNatural);
1062:   }
1063:   return(0);
1064: }
1065: #endif

1067: /*@
1068:   PetscViewerExodusIIGetId - Get the file id of the ExodusII file

1070:   Logically Collective on PetscViewer

1072:   Input Parameter:
1073: .  viewer - the PetscViewer

1075:   Output Parameter:
1076: -  exoid - The ExodusII file id

1078:   Level: intermediate

1080: .seealso: PetscViewerFileSetMode(), PetscViewerCreate(), PetscViewerSetType(), PetscViewerBinaryOpen()
1081: @*/
1082: PetscErrorCode PetscViewerExodusIIGetId(PetscViewer viewer, int *exoid)
1083: {

1088:   PetscTryMethod(viewer, "PetscViewerExodusIIGetId_C",(PetscViewer,int*),(viewer,exoid));
1089:   return(0);
1090: }

1092: /*@C
1093:   DMPlexCreateExodusFromFile - Create a DMPlex mesh from an ExodusII file.

1095:   Collective

1097:   Input Parameters:
1098: + comm  - The MPI communicator
1099: . filename - The name of the ExodusII file
1100: - interpolate - Create faces and edges in the mesh

1102:   Output Parameter:
1103: . dm  - The DM object representing the mesh

1105:   Level: beginner

1107: .seealso: DMPLEX, DMCreate(), DMPlexCreateExodus()
1108: @*/
1109: PetscErrorCode DMPlexCreateExodusFromFile(MPI_Comm comm, const char filename[], PetscBool interpolate, DM *dm)
1110: {
1111:   PetscMPIInt    rank;
1113: #if defined(PETSC_HAVE_EXODUSII)
1114:   int   CPU_word_size = sizeof(PetscReal), IO_word_size = 0, exoid = -1;
1115:   float version;
1116: #endif

1120:   MPI_Comm_rank(comm, &rank);
1121: #if defined(PETSC_HAVE_EXODUSII)
1122:   if (!rank) {
1123:     exoid = ex_open(filename, EX_READ, &CPU_word_size, &IO_word_size, &version);
1124:     if (exoid <= 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_LIB, "ex_open(\"%s\",...) did not return a valid file ID", filename);
1125:   }
1126:   DMPlexCreateExodus(comm, exoid, interpolate, dm);
1127:   if (!rank) {PetscStackCallStandard(ex_close,(exoid));}
1128: #else
1129:   SETERRQ(comm, PETSC_ERR_SUP, "This method requires ExodusII support. Reconfigure using --download-exodusii");
1130: #endif
1131:   return(0);
1132: }

1134: /*@
1135:   DMPlexCreateExodus - Create a DMPlex mesh from an ExodusII file ID.

1137:   Collective

1139:   Input Parameters:
1140: + comm  - The MPI communicator
1141: . exoid - The ExodusII id associated with a exodus file and obtained using ex_open
1142: - interpolate - Create faces and edges in the mesh

1144:   Output Parameter:
1145: . dm  - The DM object representing the mesh

1147:   Level: beginner

1149: .seealso: DMPLEX, DMCreate()
1150: @*/
1151: PetscErrorCode DMPlexCreateExodus(MPI_Comm comm, PetscInt exoid, PetscBool interpolate, DM *dm)
1152: {
1153: #if defined(PETSC_HAVE_EXODUSII)
1154:   PetscMPIInt    num_proc, rank;
1155:   PetscSection   coordSection;
1156:   Vec            coordinates;
1157:   PetscScalar    *coords;
1158:   PetscInt       coordSize, v;
1160:   /* Read from ex_get_init() */
1161:   char title[PETSC_MAX_PATH_LEN+1];
1162:   int  dim    = 0, numVertices = 0, numCells = 0, numHybridCells = 0;
1163:   int  num_cs = 0, num_vs = 0, num_fs = 0;
1164: #endif

1167: #if defined(PETSC_HAVE_EXODUSII)
1168:   MPI_Comm_rank(comm, &rank);
1169:   MPI_Comm_size(comm, &num_proc);
1170:   DMCreate(comm, dm);
1171:   DMSetType(*dm, DMPLEX);
1172:   /* Open EXODUS II file and read basic informations on rank 0, then broadcast to all processors */
1173:   if (!rank) {
1174:     PetscMemzero(title,PETSC_MAX_PATH_LEN+1);
1175:     PetscStackCallStandard(ex_get_init,(exoid, title, &dim, &numVertices, &numCells, &num_cs, &num_vs, &num_fs));
1176:     if (!num_cs) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Exodus file does not contain any cell set\n");
1177:   }
1178:   MPI_Bcast(title, PETSC_MAX_PATH_LEN+1, MPI_CHAR, 0, comm);
1179:   MPI_Bcast(&dim, 1, MPI_INT, 0, comm);
1180:   PetscObjectSetName((PetscObject) *dm, title);
1181:   DMSetDimension(*dm, dim);
1182:   DMPlexSetChart(*dm, 0, numCells+numVertices);
1183:   /*   We do not want this label automatically computed, instead we compute it here */
1184:   DMCreateLabel(*dm, "celltype");

1186:   /* Read cell sets information */
1187:   if (!rank) {
1188:     PetscInt *cone;
1189:     int      c, cs, ncs, c_loc, v, v_loc;
1190:     /* Read from ex_get_elem_blk_ids() */
1191:     int *cs_id, *cs_order;
1192:     /* Read from ex_get_elem_block() */
1193:     char buffer[PETSC_MAX_PATH_LEN+1];
1194:     int  num_cell_in_set, num_vertex_per_cell, num_hybrid, num_attr;
1195:     /* Read from ex_get_elem_conn() */
1196:     int *cs_connect;

1198:     /* Get cell sets IDs */
1199:     PetscMalloc2(num_cs, &cs_id, num_cs, &cs_order);
1200:     PetscStackCallStandard(ex_get_ids,(exoid, EX_ELEM_BLOCK, cs_id));
1201:     /* Read the cell set connectivity table and build mesh topology
1202:        EXO standard requires that cells in cell sets be numbered sequentially and be pairwise disjoint. */
1203:     /* Check for a hybrid mesh */
1204:     for (cs = 0, num_hybrid = 0; cs < num_cs; ++cs) {
1205:       PetscStackCallStandard(ex_get_block,(exoid, EX_ELEM_BLOCK, cs_id[cs], buffer, &num_cell_in_set,&num_vertex_per_cell, 0, 0, &num_attr));
1206:       switch (dim) {
1207:         case 3:
1208:         switch (num_vertex_per_cell) {
1209:           case 6:
1210:             cs_order[cs] = cs;
1211:             numHybridCells += num_cell_in_set;
1212:             ++num_hybrid;
1213:           break;
1214:           default:
1215:             for (c = cs; c > cs-num_hybrid; --c) cs_order[c] = cs_order[c-1];
1216:             cs_order[cs-num_hybrid] = cs;
1217:         }
1218:         break;
1219:       default:
1220:         for (c = cs; c > cs-num_hybrid; --c) cs_order[c] = cs_order[c-1];
1221:         cs_order[cs-num_hybrid] = cs;
1222:       }
1223:     }
1224:     /* First set sizes */
1225:     for (ncs = 0, c = 0; ncs < num_cs; ++ncs) {
1226:       const PetscInt cs = cs_order[ncs];
1227:       PetscStackCallStandard(ex_get_block,(exoid, EX_ELEM_BLOCK, cs_id[cs], buffer, &num_cell_in_set,&num_vertex_per_cell, 0, 0, &num_attr));
1228:       for (c_loc = 0; c_loc < num_cell_in_set; ++c_loc, ++c) {
1229:         DMPlexSetConeSize(*dm, c, num_vertex_per_cell);
1230:         if (c >= numCells-numHybridCells) {
1231:           DMPlexSetCellType(*dm, c, DM_POLYTOPE_TRI_PRISM);
1232:         } else {
1233:           DMPolytopeType ct;

1235:           DMPlexComputeCellType_Internal(*dm, c, 1, &ct);
1236:           DMPlexSetCellType(*dm, c, ct);
1237:         }
1238:       }
1239:     }
1240:     for (v = numCells; v < numCells+numVertices; ++v) {DMPlexSetCellType(*dm, v, DM_POLYTOPE_POINT);}
1241:     DMSetUp(*dm);
1242:     for (ncs = 0, c = 0; ncs < num_cs; ++ncs) {
1243:       const PetscInt cs = cs_order[ncs];
1244:       PetscStackCallStandard(ex_get_block,(exoid, EX_ELEM_BLOCK, cs_id[cs], buffer, &num_cell_in_set, &num_vertex_per_cell, 0, 0, &num_attr));
1245:       PetscMalloc2(num_vertex_per_cell*num_cell_in_set,&cs_connect,num_vertex_per_cell,&cone);
1246:       PetscStackCallStandard(ex_get_conn,(exoid, EX_ELEM_BLOCK, cs_id[cs], cs_connect,NULL,NULL));
1247:       /* EXO uses Fortran-based indexing, DMPlex uses C-style and numbers cell first then vertices. */
1248:       for (c_loc = 0, v = 0; c_loc < num_cell_in_set; ++c_loc, ++c) {
1249:         for (v_loc = 0; v_loc < num_vertex_per_cell; ++v_loc, ++v) {
1250:           cone[v_loc] = cs_connect[v]+numCells-1;
1251:         }
1252:         if (dim == 3) {
1253:           /* Tetrahedra are inverted */
1254:           if (num_vertex_per_cell == 4) {
1255:             PetscInt tmp = cone[0];
1256:             cone[0] = cone[1];
1257:             cone[1] = tmp;
1258:           }
1259:           /* Hexahedra are inverted */
1260:           if (num_vertex_per_cell == 8) {
1261:             PetscInt tmp = cone[1];
1262:             cone[1] = cone[3];
1263:             cone[3] = tmp;
1264:           }
1265:           /* Triangular prisms are inverted */
1266:           if (num_vertex_per_cell == 6) {
1267:             PetscInt tmp = cone[1];
1268:             cone[1] = cone[2];
1269:             cone[2] = tmp;
1270:           }
1271:         }
1272:         DMPlexSetCone(*dm, c, cone);
1273:         DMSetLabelValue(*dm, "Cell Sets", c, cs_id[cs]);
1274:       }
1275:       PetscFree2(cs_connect,cone);
1276:     }
1277:     PetscFree2(cs_id, cs_order);
1278:   }
1279:   DMPlexSymmetrize(*dm);
1280:   DMPlexStratify(*dm);
1281:   if (interpolate) {
1282:     DM idm;

1284:     DMPlexInterpolate(*dm, &idm);
1285:     DMDestroy(dm);
1286:     *dm  = idm;
1287:   }

1289:   /* Create vertex set label */
1290:   if (!rank && (num_vs > 0)) {
1291:     int vs, v;
1292:     /* Read from ex_get_node_set_ids() */
1293:     int *vs_id;
1294:     /* Read from ex_get_node_set_param() */
1295:     int num_vertex_in_set;
1296:     /* Read from ex_get_node_set() */
1297:     int *vs_vertex_list;

1299:     /* Get vertex set ids */
1300:     PetscMalloc1(num_vs, &vs_id);
1301:     PetscStackCallStandard(ex_get_ids,(exoid, EX_NODE_SET, vs_id));
1302:     for (vs = 0; vs < num_vs; ++vs) {
1303:       PetscStackCallStandard(ex_get_set_param,(exoid, EX_NODE_SET, vs_id[vs], &num_vertex_in_set, NULL));
1304:       PetscMalloc1(num_vertex_in_set, &vs_vertex_list);
1305:       PetscStackCallStandard(ex_get_set,(exoid, EX_NODE_SET, vs_id[vs], vs_vertex_list, NULL));
1306:       for (v = 0; v < num_vertex_in_set; ++v) {
1307:         DMSetLabelValue(*dm, "Vertex Sets", vs_vertex_list[v]+numCells-1, vs_id[vs]);
1308:       }
1309:       PetscFree(vs_vertex_list);
1310:     }
1311:     PetscFree(vs_id);
1312:   }
1313:   /* Read coordinates */
1314:   DMGetCoordinateSection(*dm, &coordSection);
1315:   PetscSectionSetNumFields(coordSection, 1);
1316:   PetscSectionSetFieldComponents(coordSection, 0, dim);
1317:   PetscSectionSetChart(coordSection, numCells, numCells + numVertices);
1318:   for (v = numCells; v < numCells+numVertices; ++v) {
1319:     PetscSectionSetDof(coordSection, v, dim);
1320:     PetscSectionSetFieldDof(coordSection, v, 0, dim);
1321:   }
1322:   PetscSectionSetUp(coordSection);
1323:   PetscSectionGetStorageSize(coordSection, &coordSize);
1324:   VecCreate(PETSC_COMM_SELF, &coordinates);
1325:   PetscObjectSetName((PetscObject) coordinates, "coordinates");
1326:   VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);
1327:   VecSetBlockSize(coordinates, dim);
1328:   VecSetType(coordinates,VECSTANDARD);
1329:   VecGetArray(coordinates, &coords);
1330:   if (!rank) {
1331:     PetscReal *x, *y, *z;

1333:     PetscMalloc3(numVertices,&x,numVertices,&y,numVertices,&z);
1334:     PetscStackCallStandard(ex_get_coord,(exoid, x, y, z));
1335:     if (dim > 0) {
1336:       for (v = 0; v < numVertices; ++v) coords[v*dim+0] = x[v];
1337:     }
1338:     if (dim > 1) {
1339:       for (v = 0; v < numVertices; ++v) coords[v*dim+1] = y[v];
1340:     }
1341:     if (dim > 2) {
1342:       for (v = 0; v < numVertices; ++v) coords[v*dim+2] = z[v];
1343:     }
1344:     PetscFree3(x,y,z);
1345:   }
1346:   VecRestoreArray(coordinates, &coords);
1347:   DMSetCoordinatesLocal(*dm, coordinates);
1348:   VecDestroy(&coordinates);

1350:   /* Create side set label */
1351:   if (!rank && interpolate && (num_fs > 0)) {
1352:     int fs, f, voff;
1353:     /* Read from ex_get_side_set_ids() */
1354:     int *fs_id;
1355:     /* Read from ex_get_side_set_param() */
1356:     int num_side_in_set;
1357:     /* Read from ex_get_side_set_node_list() */
1358:     int *fs_vertex_count_list, *fs_vertex_list;
1359:     /* Read side set labels */
1360:     char fs_name[MAX_STR_LENGTH+1];

1362:     /* Get side set ids */
1363:     PetscMalloc1(num_fs, &fs_id);
1364:     PetscStackCallStandard(ex_get_ids,(exoid, EX_SIDE_SET, fs_id));
1365:     for (fs = 0; fs < num_fs; ++fs) {
1366:       PetscStackCallStandard(ex_get_set_param,(exoid, EX_SIDE_SET, fs_id[fs], &num_side_in_set, NULL));
1367:       PetscMalloc2(num_side_in_set,&fs_vertex_count_list,num_side_in_set*4,&fs_vertex_list);
1368:       PetscStackCallStandard(ex_get_side_set_node_list,(exoid, fs_id[fs], fs_vertex_count_list, fs_vertex_list));
1369:       /* Get the specific name associated with this side set ID. */
1370:       int fs_name_err = ex_get_name(exoid, EX_SIDE_SET, fs_id[fs], fs_name);
1371:       for (f = 0, voff = 0; f < num_side_in_set; ++f) {
1372:         const PetscInt *faces   = NULL;
1373:         PetscInt       faceSize = fs_vertex_count_list[f], numFaces;
1374:         PetscInt       faceVertices[4], v;

1376:         if (faceSize > 4) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "ExodusII side cannot have %d > 4 vertices", faceSize);
1377:         for (v = 0; v < faceSize; ++v, ++voff) {
1378:           faceVertices[v] = fs_vertex_list[voff]+numCells-1;
1379:         }
1380:         DMPlexGetFullJoin(*dm, faceSize, faceVertices, &numFaces, &faces);
1381:         if (numFaces != 1) SETERRQ3(comm, PETSC_ERR_ARG_WRONG, "Invalid ExodusII side %d in set %d maps to %d faces", f, fs, numFaces);
1382:         DMSetLabelValue(*dm, "Face Sets", faces[0], fs_id[fs]);
1383:         /* Only add the label if one has been detected for this side set. */
1384:         if (!fs_name_err) {
1385:           DMSetLabelValue(*dm, fs_name, faces[0], fs_id[fs]);
1386:         }
1387:         DMPlexRestoreJoin(*dm, faceSize, faceVertices, &numFaces, &faces);
1388:       }
1389:       PetscFree2(fs_vertex_count_list,fs_vertex_list);
1390:     }
1391:     PetscFree(fs_id);
1392:   }
1393: #else
1394:   SETERRQ(comm, PETSC_ERR_SUP, "This method requires ExodusII support. Reconfigure using --download-exodusii");
1395: #endif
1396:   return(0);
1397: }