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, °ree));
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;
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(PetscMalloc2(num_side_in_set, &fs_vertex_count_list, num_side_in_set * 4, &fs_vertex_list));
1698: PetscCallExternal(ex_get_side_set_node_list, exoid, fs_id[fs], fs_vertex_count_list, fs_vertex_list);
1699: /* Get the specific name associated with this side set ID. */
1700: int fs_name_err = ex_get_name(exoid, EX_SIDE_SET, fs_id[fs], fs_name);
1701: if (!fs_name_err) {
1702: PetscCall(PetscStrlen(fs_name, &fs_name_len));
1703: if (fs_name_len == 0) PetscCall(PetscStrncpy(fs_name, "Face Sets", MAX_STR_LENGTH));
1704: }
1705: PetscCheck(fs_id[fs] != 1 && fs_id[fs] != 2, comm, PETSC_ERR_ARG_WRONG, "Side set %s marker cannot be %d since this is reserved by ExodusII", fs_name, fs_id[fs]);
1706: for (f = 0, voff = 0; f < num_side_in_set; ++f) {
1707: const PetscInt *faces = NULL;
1708: PetscInt faceSize = fs_vertex_count_list[f], numFaces;
1709: PetscInt faceVertices[4], v;
1711: PetscCheck(faceSize <= 4, comm, PETSC_ERR_ARG_WRONG, "ExodusII side cannot have %" PetscInt_FMT " > 4 vertices", faceSize);
1712: for (v = 0; v < faceSize; ++v, ++voff) faceVertices[v] = fs_vertex_list[voff] + numCells - 1;
1713: PetscCall(DMPlexGetFullJoin(*dm, faceSize, faceVertices, &numFaces, &faces));
1714: PetscCheck(numFaces == 1, comm, PETSC_ERR_ARG_WRONG, "Invalid ExodusII side %d in set %d maps to %" PetscInt_FMT " faces", f, fs, numFaces);
1715: PetscCall(DMSetLabelValue_Fast(*dm, &faceSets, "Face Sets", faces[0], fs_id[fs]));
1716: /* Only add the label if one has been detected for this side set. */
1717: if (!fs_name_err) PetscCall(DMSetLabelValue(*dm, fs_name, faces[0], fs_id[fs]));
1718: PetscCall(DMPlexRestoreJoin(*dm, faceSize, faceVertices, &numFaces, &faces));
1719: }
1720: PetscCall(PetscFree2(fs_vertex_count_list, fs_vertex_list));
1721: }
1722: PetscCall(PetscFree(fs_id));
1723: }
1725: { /* Create Cell/Face/Vertex Sets labels at all processes */
1726: enum {
1727: n = 3
1728: };
1729: PetscBool flag[n];
1731: flag[0] = cellSets ? PETSC_TRUE : PETSC_FALSE;
1732: flag[1] = faceSets ? PETSC_TRUE : PETSC_FALSE;
1733: flag[2] = vertSets ? PETSC_TRUE : PETSC_FALSE;
1734: PetscCallMPI(MPI_Bcast(flag, n, MPIU_BOOL, 0, comm));
1735: if (flag[0]) PetscCall(DMCreateLabel(*dm, "Cell Sets"));
1736: if (flag[1]) PetscCall(DMCreateLabel(*dm, "Face Sets"));
1737: if (flag[2]) PetscCall(DMCreateLabel(*dm, "Vertex Sets"));
1738: }
1739: PetscFunctionReturn(PETSC_SUCCESS);
1740: #else
1741: SETERRQ(comm, PETSC_ERR_SUP, "This method requires ExodusII support. Reconfigure using --download-exodusii");
1742: #endif
1743: }