Actual source code: plexexodusii.c
1: #define PETSCDM_DLL
2: #include <petsc/private/dmpleximpl.h>
4: #if defined(PETSC_HAVE_EXODUSII)
5: #include <netcdf.h>
6: #include <exodusII.h>
7: #endif
9: #include <petsc/private/viewerimpl.h>
10: #include <petsc/private/viewerexodusiiimpl.h>
11: #if defined(PETSC_HAVE_EXODUSII)
12: /*@C
13: PETSC_VIEWER_EXODUSII_ - Creates an `PETSCVIEWEREXODUSII` `PetscViewer` shared by all processors in a communicator.
15: Collective
17: Input Parameter:
18: . comm - the MPI communicator to share the `PETSCVIEWEREXODUSII` `PetscViewer`
20: Level: intermediate
22: Note:
23: Unlike almost all other PETSc routines, PETSC_VIEWER_EXODUSII_ does not return
24: an error code. The GLVIS PetscViewer is usually used in the form
25: $ XXXView(XXX object, PETSC_VIEWER_EXODUSII_(comm));
27: Fortran Note:
28: No support in Fortran
30: .seealso: `PETSCVIEWEREXODUSII`, `PetscViewer`, `PetscViewer`, `PetscViewerExodusIIOpen()`, `PetscViewerType`, `PetscViewerCreate()`, `PetscViewerDestroy()`
31: @*/
32: PetscViewer PETSC_VIEWER_EXODUSII_(MPI_Comm comm)
33: {
34: PetscViewer viewer;
37: PetscViewerExodusIIOpen(comm, "mesh.exo", FILE_MODE_WRITE, &viewer);
38: if (ierr) {
39: PetscError(PETSC_COMM_SELF, __LINE__, "PETSC_VIEWER_EXODUSII_", __FILE__, PETSC_ERR_PLIB, PETSC_ERROR_INITIAL, " ");
40: return NULL;
41: }
42: PetscObjectRegisterDestroy((PetscObject)viewer);
43: if (ierr) {
44: PetscError(PETSC_COMM_SELF, __LINE__, "PETSC_VIEWER_EXODUSII_", __FILE__, PETSC_ERR_PLIB, PETSC_ERROR_INITIAL, " ");
45: return NULL;
46: }
47: return viewer;
48: }
50: static PetscErrorCode PetscViewerView_ExodusII(PetscViewer v, PetscViewer viewer)
51: {
52: PetscViewer_ExodusII *exo = (PetscViewer_ExodusII *)v->data;
54: if (exo->filename) PetscViewerASCIIPrintf(viewer, "Filename: %s\n", exo->filename);
55: if (exo->exoid) PetscViewerASCIIPrintf(viewer, "exoid: %d\n", exo->exoid);
56: if (exo->btype) PetscViewerASCIIPrintf(viewer, "IO Mode: %d\n", exo->btype);
57: if (exo->order) PetscViewerASCIIPrintf(viewer, "Mesh order: %" PetscInt_FMT "\n", exo->order);
58: return 0;
59: }
61: static PetscErrorCode PetscViewerSetFromOptions_ExodusII(PetscViewer v, PetscOptionItems *PetscOptionsObject)
62: {
63: PetscOptionsHeadBegin(PetscOptionsObject, "ExodusII PetscViewer Options");
64: PetscOptionsHeadEnd();
65: return 0;
66: }
68: static PetscErrorCode PetscViewerSetUp_ExodusII(PetscViewer viewer)
69: {
70: return 0;
71: }
73: static PetscErrorCode PetscViewerDestroy_ExodusII(PetscViewer viewer)
74: {
75: PetscViewer_ExodusII *exo = (PetscViewer_ExodusII *)viewer->data;
77: if (exo->exoid >= 0) PetscCallExternal(ex_close, exo->exoid);
78: PetscFree(exo->filename);
79: PetscFree(exo);
80: PetscObjectComposeFunction((PetscObject)viewer, "PetscViewerFileSetName_C", NULL);
81: PetscObjectComposeFunction((PetscObject)viewer, "PetscViewerFileGetName_C", NULL);
82: PetscObjectComposeFunction((PetscObject)viewer, "PetscViewerFileSetMode_C", NULL);
83: PetscObjectComposeFunction((PetscObject)viewer, "PetscViewerFileGetMode_C", NULL);
84: PetscObjectComposeFunction((PetscObject)viewer, "PetscViewerGetId_C", NULL);
85: PetscObjectComposeFunction((PetscObject)viewer, "PetscViewerGetOrder_C", NULL);
86: PetscObjectComposeFunction((PetscObject)viewer, "PetscViewerSetOrder_C", NULL);
87: return 0;
88: }
90: static PetscErrorCode PetscViewerFileSetName_ExodusII(PetscViewer viewer, const char name[])
91: {
92: PetscViewer_ExodusII *exo = (PetscViewer_ExodusII *)viewer->data;
93: PetscMPIInt rank;
94: int CPU_word_size, IO_word_size, EXO_mode;
95: MPI_Info mpi_info = MPI_INFO_NULL;
96: float EXO_version;
98: MPI_Comm_rank(PetscObjectComm((PetscObject)viewer), &rank);
99: CPU_word_size = sizeof(PetscReal);
100: IO_word_size = sizeof(PetscReal);
102: if (exo->exoid >= 0) {
103: PetscCallExternal(ex_close, exo->exoid);
104: exo->exoid = -1;
105: }
106: if (exo->filename) PetscFree(exo->filename);
107: PetscStrallocpy(name, &exo->filename);
108: switch (exo->btype) {
109: case FILE_MODE_READ:
110: EXO_mode = EX_READ;
111: break;
112: case FILE_MODE_APPEND:
113: case FILE_MODE_UPDATE:
114: case FILE_MODE_APPEND_UPDATE:
115: /* Will fail if the file does not already exist */
116: EXO_mode = EX_WRITE;
117: break;
118: case FILE_MODE_WRITE:
119: /*
120: exodus only allows writing geometry upon file creation, so we will let DMView create the file.
121: */
122: return 0;
123: break;
124: default:
125: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ORDER, "Must call PetscViewerFileSetMode() before PetscViewerFileSetName()");
126: }
127: #if defined(PETSC_USE_64BIT_INDICES)
128: EXO_mode += EX_ALL_INT64_API;
129: #endif
130: exo->exoid = ex_open_par(name, EXO_mode, &CPU_word_size, &IO_word_size, &EXO_version, PETSC_COMM_WORLD, mpi_info);
132: return 0;
133: }
135: static PetscErrorCode PetscViewerFileGetName_ExodusII(PetscViewer viewer, const char **name)
136: {
137: PetscViewer_ExodusII *exo = (PetscViewer_ExodusII *)viewer->data;
139: *name = exo->filename;
140: return 0;
141: }
143: static PetscErrorCode PetscViewerFileSetMode_ExodusII(PetscViewer viewer, PetscFileMode type)
144: {
145: PetscViewer_ExodusII *exo = (PetscViewer_ExodusII *)viewer->data;
147: exo->btype = type;
148: return 0;
149: }
151: static PetscErrorCode PetscViewerFileGetMode_ExodusII(PetscViewer viewer, PetscFileMode *type)
152: {
153: PetscViewer_ExodusII *exo = (PetscViewer_ExodusII *)viewer->data;
155: *type = exo->btype;
156: return 0;
157: }
159: static PetscErrorCode PetscViewerExodusIIGetId_ExodusII(PetscViewer viewer, int *exoid)
160: {
161: PetscViewer_ExodusII *exo = (PetscViewer_ExodusII *)viewer->data;
163: *exoid = exo->exoid;
164: return 0;
165: }
167: static PetscErrorCode PetscViewerExodusIIGetOrder_ExodusII(PetscViewer viewer, PetscInt *order)
168: {
169: PetscViewer_ExodusII *exo = (PetscViewer_ExodusII *)viewer->data;
171: *order = exo->order;
172: return 0;
173: }
175: static PetscErrorCode PetscViewerExodusIISetOrder_ExodusII(PetscViewer viewer, PetscInt order)
176: {
177: PetscViewer_ExodusII *exo = (PetscViewer_ExodusII *)viewer->data;
179: exo->order = order;
180: return 0;
181: }
183: /*MC
184: PETSCVIEWEREXODUSII - A viewer that writes to an Exodus II file
186: .seealso: `PetscViewerExodusIIOpen()`, `PetscViewerCreate()`, `PETSCVIEWERBINARY`, `PETSCVIEWERHDF5`, `DMView()`,
187: `PetscViewerFileSetName()`, `PetscViewerFileSetMode()`, `PetscViewerFormat`, `PetscViewerType`, `PetscViewerSetType()`
189: Level: beginner
190: M*/
192: PETSC_EXTERN PetscErrorCode PetscViewerCreate_ExodusII(PetscViewer v)
193: {
194: PetscViewer_ExodusII *exo;
196: PetscNew(&exo);
198: v->data = (void *)exo;
199: v->ops->destroy = PetscViewerDestroy_ExodusII;
200: v->ops->setfromoptions = PetscViewerSetFromOptions_ExodusII;
201: v->ops->setup = PetscViewerSetUp_ExodusII;
202: v->ops->view = PetscViewerView_ExodusII;
203: v->ops->flush = 0;
204: exo->btype = (PetscFileMode)-1;
205: exo->filename = 0;
206: exo->exoid = -1;
208: PetscObjectComposeFunction((PetscObject)v, "PetscViewerFileSetName_C", PetscViewerFileSetName_ExodusII);
209: PetscObjectComposeFunction((PetscObject)v, "PetscViewerFileGetName_C", PetscViewerFileGetName_ExodusII);
210: PetscObjectComposeFunction((PetscObject)v, "PetscViewerFileSetMode_C", PetscViewerFileSetMode_ExodusII);
211: PetscObjectComposeFunction((PetscObject)v, "PetscViewerFileGetMode_C", PetscViewerFileGetMode_ExodusII);
212: PetscObjectComposeFunction((PetscObject)v, "PetscViewerGetId_C", PetscViewerExodusIIGetId_ExodusII);
213: PetscObjectComposeFunction((PetscObject)v, "PetscViewerSetOrder_C", PetscViewerExodusIISetOrder_ExodusII);
214: PetscObjectComposeFunction((PetscObject)v, "PetscViewerGetOrder_C", PetscViewerExodusIIGetOrder_ExodusII);
215: return 0;
216: }
218: /*
219: EXOGetVarIndex - Locate a result in an exodus file based on its name
221: Collective
223: Input Parameters:
224: + exoid - the exodus id of a file (obtained from ex_open or ex_create for instance)
225: . obj_type - the type of entity for instance EX_NODAL, EX_ELEM_BLOCK
226: - name - the name of the result
228: Output Parameters:
229: . varIndex - the location in the exodus file of the result
231: Notes:
232: The exodus variable index is obtained by comparing name and the
233: names of zonal variables declared in the exodus file. For instance if name is "V"
234: the location in the exodus file will be the first match of "V", "V_X", "V_XX", "V_1", or "V_11"
235: amongst all variables of type obj_type.
237: Level: beginner
239: .seealso: `EXOGetVarIndex()`, `DMPlexView_ExodusII_Internal()`, `VecViewPlex_ExodusII_Nodal_Internal()`, `VecLoadNodal_PlexEXO()`, `VecLoadZonal_PlexEXO()`
240: */
241: PetscErrorCode EXOGetVarIndex_Internal(int exoid, ex_entity_type obj_type, const char name[], int *varIndex)
242: {
243: int num_vars, i, j;
244: char ext_name[MAX_STR_LENGTH + 1], var_name[MAX_STR_LENGTH + 1];
245: const int num_suffix = 5;
246: char *suffix[5];
247: PetscBool flg;
249: suffix[0] = (char *)"";
250: suffix[1] = (char *)"_X";
251: suffix[2] = (char *)"_XX";
252: suffix[3] = (char *)"_1";
253: suffix[4] = (char *)"_11";
255: *varIndex = -1;
256: PetscCallExternal(ex_get_variable_param, exoid, obj_type, &num_vars);
257: for (i = 0; i < num_vars; ++i) {
258: PetscCallExternal(ex_get_variable_name, exoid, obj_type, i + 1, var_name);
259: for (j = 0; j < num_suffix; ++j) {
260: PetscStrncpy(ext_name, name, MAX_STR_LENGTH);
261: PetscStrlcat(ext_name, suffix[j], MAX_STR_LENGTH);
262: PetscStrcasecmp(ext_name, var_name, &flg);
263: if (flg) *varIndex = i + 1;
264: }
265: }
266: return 0;
267: }
269: /*
270: DMView_PlexExodusII - Write a DM to disk in exodus format
272: Collective on dm
274: Input Parameters:
275: + dm - The dm to be written
276: . viewer - an exodusII viewer
278: Notes:
279: Not all DM can be written to disk this way. For instance, exodus assume that element blocks (mapped to "Cell sets" labels)
280: consists of sequentially numbered cells. If this is not the case, the exodus file will be corrupted.
282: If the dm has been distributed, only the part of the DM on MPI rank 0 (including "ghost" cells and vertices)
283: will be written.
285: DMPlex only represents geometry while most post-processing software expect that a mesh also provides information
286: on the discretization space. This function assumes that the file represents Lagrange finite elements of order 1 or 2.
287: The order of the mesh shall be set using PetscViewerExodusIISetOrder
288: It should be extended to use PetscFE objects.
290: This function will only handle TRI, TET, QUAD, and HEX cells.
291: Level: beginner
293: .seealso:
294: */
295: PetscErrorCode DMView_PlexExodusII(DM dm, PetscViewer viewer)
296: {
297: enum ElemType {
298: TRI,
299: QUAD,
300: TET,
301: HEX
302: };
303: MPI_Comm comm;
304: PetscInt degree; /* the order of the mesh */
305: /* Connectivity Variables */
306: PetscInt cellsNotInConnectivity;
307: /* Cell Sets */
308: DMLabel csLabel;
309: IS csIS;
310: const PetscInt *csIdx;
311: PetscInt num_cs, cs;
312: enum ElemType *type;
313: PetscBool hasLabel;
314: /* Coordinate Variables */
315: DM cdm;
316: PetscSection coordSection;
317: Vec coord;
318: PetscInt **nodes;
319: PetscInt depth, d, dim, skipCells = 0;
320: PetscInt pStart, pEnd, p, cStart, cEnd, numCells, vStart, vEnd, numVertices, eStart, eEnd, numEdges, fStart, fEnd, numFaces, numNodes;
321: PetscInt num_vs, num_fs;
322: PetscMPIInt rank, size;
323: const char *dmName;
324: PetscInt nodesTriP1[4] = {3, 0, 0, 0};
325: PetscInt nodesTriP2[4] = {3, 3, 0, 0};
326: PetscInt nodesQuadP1[4] = {4, 0, 0, 0};
327: PetscInt nodesQuadP2[4] = {4, 4, 0, 1};
328: PetscInt nodesTetP1[4] = {4, 0, 0, 0};
329: PetscInt nodesTetP2[4] = {4, 6, 0, 0};
330: PetscInt nodesHexP1[4] = {8, 0, 0, 0};
331: PetscInt nodesHexP2[4] = {8, 12, 6, 1};
332: int CPU_word_size, IO_word_size, EXO_mode;
333: float EXO_version;
335: PetscViewer_ExodusII *exo = (PetscViewer_ExodusII *)viewer->data;
337: PetscObjectGetComm((PetscObject)dm, &comm);
338: MPI_Comm_rank(comm, &rank);
339: MPI_Comm_size(comm, &size);
341: /*
342: Creating coordSection is a collective operation so we do it somewhat out of sequence
343: */
344: PetscSectionCreate(comm, &coordSection);
345: DMGetCoordinatesLocalSetUp(dm);
346: if (rank == 0) {
347: switch (exo->btype) {
348: case FILE_MODE_READ:
349: case FILE_MODE_APPEND:
350: case FILE_MODE_UPDATE:
351: case FILE_MODE_APPEND_UPDATE:
352: /* exodusII does not allow writing geometry to an existing file */
353: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "cannot add geometry to existing file %s", exo->filename);
354: case FILE_MODE_WRITE:
355: /* Create an empty file if one already exists*/
356: EXO_mode = EX_CLOBBER;
357: #if defined(PETSC_USE_64BIT_INDICES)
358: EXO_mode += EX_ALL_INT64_API;
359: #endif
360: CPU_word_size = sizeof(PetscReal);
361: IO_word_size = sizeof(PetscReal);
362: exo->exoid = ex_create(exo->filename, EXO_mode, &CPU_word_size, &IO_word_size);
365: break;
366: default:
367: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ORDER, "Must call PetscViewerFileSetMode() before PetscViewerFileSetName()");
368: }
370: /* --- Get DM info --- */
371: PetscObjectGetName((PetscObject)dm, &dmName);
372: DMPlexGetDepth(dm, &depth);
373: DMGetDimension(dm, &dim);
374: DMPlexGetChart(dm, &pStart, &pEnd);
375: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
376: DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);
377: DMPlexGetDepthStratum(dm, 1, &eStart, &eEnd);
378: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
379: numCells = cEnd - cStart;
380: numEdges = eEnd - eStart;
381: numVertices = vEnd - vStart;
382: if (depth == 3) {
383: numFaces = fEnd - fStart;
384: } else {
385: numFaces = 0;
386: }
387: DMGetLabelSize(dm, "Cell Sets", &num_cs);
388: DMGetLabelSize(dm, "Vertex Sets", &num_vs);
389: DMGetLabelSize(dm, "Face Sets", &num_fs);
390: DMGetCoordinatesLocal(dm, &coord);
391: DMGetCoordinateDM(dm, &cdm);
392: if (num_cs > 0) {
393: DMGetLabel(dm, "Cell Sets", &csLabel);
394: DMLabelGetValueIS(csLabel, &csIS);
395: ISGetIndices(csIS, &csIdx);
396: }
397: PetscMalloc1(num_cs, &nodes);
398: /* Set element type for each block and compute total number of nodes */
399: PetscMalloc1(num_cs, &type);
400: numNodes = numVertices;
402: PetscViewerExodusIIGetOrder(viewer, °ree);
403: if (degree == 2) numNodes += numEdges;
404: cellsNotInConnectivity = numCells;
405: for (cs = 0; cs < num_cs; ++cs) {
406: IS stratumIS;
407: const PetscInt *cells;
408: PetscScalar *xyz = NULL;
409: PetscInt csSize, closureSize;
411: DMLabelGetStratumIS(csLabel, csIdx[cs], &stratumIS);
412: ISGetIndices(stratumIS, &cells);
413: ISGetSize(stratumIS, &csSize);
414: DMPlexVecGetClosure(cdm, NULL, coord, cells[0], &closureSize, &xyz);
415: switch (dim) {
416: case 2:
417: if (closureSize == 3 * dim) {
418: type[cs] = TRI;
419: } else if (closureSize == 4 * dim) {
420: type[cs] = QUAD;
421: } 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);
422: break;
423: case 3:
424: if (closureSize == 4 * dim) {
425: type[cs] = TET;
426: } else if (closureSize == 8 * dim) {
427: type[cs] = HEX;
428: } 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);
429: break;
430: default:
431: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Dimension %" PetscInt_FMT " not handled by ExodusII viewer", dim);
432: }
433: if ((degree == 2) && (type[cs] == QUAD)) numNodes += csSize;
434: if ((degree == 2) && (type[cs] == HEX)) {
435: numNodes += csSize;
436: numNodes += numFaces;
437: }
438: DMPlexVecRestoreClosure(cdm, NULL, coord, cells[0], &closureSize, &xyz);
439: /* Set nodes and Element type */
440: if (type[cs] == TRI) {
441: if (degree == 1) nodes[cs] = nodesTriP1;
442: else if (degree == 2) nodes[cs] = nodesTriP2;
443: } else if (type[cs] == QUAD) {
444: if (degree == 1) nodes[cs] = nodesQuadP1;
445: else if (degree == 2) nodes[cs] = nodesQuadP2;
446: } else if (type[cs] == TET) {
447: if (degree == 1) nodes[cs] = nodesTetP1;
448: else if (degree == 2) nodes[cs] = nodesTetP2;
449: } else if (type[cs] == HEX) {
450: if (degree == 1) nodes[cs] = nodesHexP1;
451: else if (degree == 2) nodes[cs] = nodesHexP2;
452: }
453: /* Compute the number of cells not in the connectivity table */
454: cellsNotInConnectivity -= nodes[cs][3] * csSize;
456: ISRestoreIndices(stratumIS, &cells);
457: ISDestroy(&stratumIS);
458: }
459: if (num_cs > 0) PetscCallExternal(ex_put_init, exo->exoid, dmName, dim, numNodes, numCells, num_cs, num_vs, num_fs);
460: /* --- Connectivity --- */
461: for (cs = 0; cs < num_cs; ++cs) {
462: IS stratumIS;
463: const PetscInt *cells;
464: PetscInt *connect, off = 0;
465: PetscInt edgesInClosure = 0, facesInClosure = 0, verticesInClosure = 0;
466: PetscInt csSize, c, connectSize, closureSize;
467: char *elem_type = NULL;
468: char elem_type_tri3[] = "TRI3", elem_type_quad4[] = "QUAD4";
469: char elem_type_tri6[] = "TRI6", elem_type_quad9[] = "QUAD9";
470: char elem_type_tet4[] = "TET4", elem_type_hex8[] = "HEX8";
471: char elem_type_tet10[] = "TET10", elem_type_hex27[] = "HEX27";
473: DMLabelGetStratumIS(csLabel, csIdx[cs], &stratumIS);
474: ISGetIndices(stratumIS, &cells);
475: ISGetSize(stratumIS, &csSize);
476: /* Set Element type */
477: if (type[cs] == TRI) {
478: if (degree == 1) elem_type = elem_type_tri3;
479: else if (degree == 2) elem_type = elem_type_tri6;
480: } else if (type[cs] == QUAD) {
481: if (degree == 1) elem_type = elem_type_quad4;
482: else if (degree == 2) elem_type = elem_type_quad9;
483: } else if (type[cs] == TET) {
484: if (degree == 1) elem_type = elem_type_tet4;
485: else if (degree == 2) elem_type = elem_type_tet10;
486: } else if (type[cs] == HEX) {
487: if (degree == 1) elem_type = elem_type_hex8;
488: else if (degree == 2) elem_type = elem_type_hex27;
489: }
490: connectSize = nodes[cs][0] + nodes[cs][1] + nodes[cs][2] + nodes[cs][3];
491: PetscMalloc1(PetscMax(27, connectSize) * csSize, &connect);
492: PetscCallExternal(ex_put_block, exo->exoid, EX_ELEM_BLOCK, csIdx[cs], elem_type, csSize, connectSize, 0, 0, 1);
493: /* Find number of vertices, edges, and faces in the closure */
494: verticesInClosure = nodes[cs][0];
495: if (depth > 1) {
496: if (dim == 2) {
497: DMPlexGetConeSize(dm, cells[0], &edgesInClosure);
498: } else if (dim == 3) {
499: PetscInt *closure = NULL;
501: DMPlexGetConeSize(dm, cells[0], &facesInClosure);
502: DMPlexGetTransitiveClosure(dm, cells[0], PETSC_TRUE, &closureSize, &closure);
503: edgesInClosure = closureSize - facesInClosure - 1 - verticesInClosure;
504: DMPlexRestoreTransitiveClosure(dm, cells[0], PETSC_TRUE, &closureSize, &closure);
505: }
506: }
507: /* Get connectivity for each cell */
508: for (c = 0; c < csSize; ++c) {
509: PetscInt *closure = NULL;
510: PetscInt temp, i;
512: DMPlexGetTransitiveClosure(dm, cells[c], PETSC_TRUE, &closureSize, &closure);
513: for (i = 0; i < connectSize; ++i) {
514: if (i < nodes[cs][0]) { /* Vertices */
515: connect[i + off] = closure[(i + edgesInClosure + facesInClosure + 1) * 2] + 1;
516: connect[i + off] -= cellsNotInConnectivity;
517: } else if (i < nodes[cs][0] + nodes[cs][1]) { /* Edges */
518: connect[i + off] = closure[(i - verticesInClosure + facesInClosure + 1) * 2] + 1;
519: if (nodes[cs][2] == 0) connect[i + off] -= numFaces;
520: connect[i + off] -= cellsNotInConnectivity;
521: } else if (i < nodes[cs][0] + nodes[cs][1] + nodes[cs][3]) { /* Cells */
522: connect[i + off] = closure[0] + 1;
523: connect[i + off] -= skipCells;
524: } else if (i < nodes[cs][0] + nodes[cs][1] + nodes[cs][3] + nodes[cs][2]) { /* Faces */
525: connect[i + off] = closure[(i - edgesInClosure - verticesInClosure) * 2] + 1;
526: connect[i + off] -= cellsNotInConnectivity;
527: } else {
528: connect[i + off] = -1;
529: }
530: }
531: /* Tetrahedra are inverted */
532: if (type[cs] == TET) {
533: temp = connect[0 + off];
534: connect[0 + off] = connect[1 + off];
535: connect[1 + off] = temp;
536: if (degree == 2) {
537: temp = connect[5 + off];
538: connect[5 + off] = connect[6 + off];
539: connect[6 + off] = temp;
540: temp = connect[7 + off];
541: connect[7 + off] = connect[8 + off];
542: connect[8 + off] = temp;
543: }
544: }
545: /* Hexahedra are inverted */
546: if (type[cs] == HEX) {
547: temp = connect[1 + off];
548: connect[1 + off] = connect[3 + off];
549: connect[3 + off] = temp;
550: if (degree == 2) {
551: temp = connect[8 + off];
552: connect[8 + off] = connect[11 + off];
553: connect[11 + off] = temp;
554: temp = connect[9 + off];
555: connect[9 + off] = connect[10 + off];
556: connect[10 + off] = temp;
557: temp = connect[16 + off];
558: connect[16 + off] = connect[17 + off];
559: connect[17 + off] = temp;
560: temp = connect[18 + off];
561: connect[18 + off] = connect[19 + off];
562: connect[19 + off] = temp;
564: temp = connect[12 + off];
565: connect[12 + off] = connect[16 + off];
566: connect[16 + off] = temp;
567: temp = connect[13 + off];
568: connect[13 + off] = connect[17 + off];
569: connect[17 + off] = temp;
570: temp = connect[14 + off];
571: connect[14 + off] = connect[18 + off];
572: connect[18 + off] = temp;
573: temp = connect[15 + off];
574: connect[15 + off] = connect[19 + off];
575: connect[19 + off] = temp;
577: temp = connect[23 + off];
578: connect[23 + off] = connect[26 + off];
579: connect[26 + off] = temp;
580: temp = connect[24 + off];
581: connect[24 + off] = connect[25 + off];
582: connect[25 + off] = temp;
583: temp = connect[25 + off];
584: connect[25 + off] = connect[26 + off];
585: connect[26 + off] = temp;
586: }
587: }
588: off += connectSize;
589: DMPlexRestoreTransitiveClosure(dm, cells[c], PETSC_TRUE, &closureSize, &closure);
590: }
591: PetscCallExternal(ex_put_conn, exo->exoid, EX_ELEM_BLOCK, csIdx[cs], connect, 0, 0);
592: skipCells += (nodes[cs][3] == 0) * csSize;
593: PetscFree(connect);
594: ISRestoreIndices(stratumIS, &cells);
595: ISDestroy(&stratumIS);
596: }
597: PetscFree(type);
598: /* --- Coordinates --- */
599: PetscSectionSetChart(coordSection, pStart, pEnd);
600: if (num_cs) {
601: for (d = 0; d < depth; ++d) {
602: DMPlexGetDepthStratum(dm, d, &pStart, &pEnd);
603: for (p = pStart; p < pEnd; ++p) PetscSectionSetDof(coordSection, p, nodes[0][d] > 0);
604: }
605: }
606: for (cs = 0; cs < num_cs; ++cs) {
607: IS stratumIS;
608: const PetscInt *cells;
609: PetscInt csSize, c;
611: DMLabelGetStratumIS(csLabel, csIdx[cs], &stratumIS);
612: ISGetIndices(stratumIS, &cells);
613: ISGetSize(stratumIS, &csSize);
614: for (c = 0; c < csSize; ++c) PetscSectionSetDof(coordSection, cells[c], nodes[cs][3] > 0);
615: ISRestoreIndices(stratumIS, &cells);
616: ISDestroy(&stratumIS);
617: }
618: if (num_cs > 0) {
619: ISRestoreIndices(csIS, &csIdx);
620: ISDestroy(&csIS);
621: }
622: PetscFree(nodes);
623: PetscSectionSetUp(coordSection);
624: if (numNodes > 0) {
625: const char *coordNames[3] = {"x", "y", "z"};
626: PetscScalar *closure, *cval;
627: PetscReal *coords;
628: PetscInt hasDof, n = 0;
630: /* There can't be more than 24 values in the closure of a point for the coord coordSection */
631: PetscCalloc3(numNodes * 3, &coords, dim, &cval, 24, &closure);
632: DMGetCoordinatesLocalNoncollective(dm, &coord);
633: DMPlexGetChart(dm, &pStart, &pEnd);
634: for (p = pStart; p < pEnd; ++p) {
635: PetscSectionGetDof(coordSection, p, &hasDof);
636: if (hasDof) {
637: PetscInt closureSize = 24, j;
639: DMPlexVecGetClosure(cdm, NULL, coord, p, &closureSize, &closure);
640: for (d = 0; d < dim; ++d) {
641: cval[d] = 0.0;
642: for (j = 0; j < closureSize / dim; j++) cval[d] += closure[j * dim + d];
643: coords[d * numNodes + n] = PetscRealPart(cval[d]) * dim / closureSize;
644: }
645: ++n;
646: }
647: }
648: PetscCallExternal(ex_put_coord, exo->exoid, &coords[0 * numNodes], &coords[1 * numNodes], &coords[2 * numNodes]);
649: PetscFree3(coords, cval, closure);
650: PetscCallExternal(ex_put_coord_names, exo->exoid, (char **)coordNames);
651: }
653: /* --- Node Sets/Vertex Sets --- */
654: DMHasLabel(dm, "Vertex Sets", &hasLabel);
655: if (hasLabel) {
656: PetscInt i, vs, vsSize;
657: const PetscInt *vsIdx, *vertices;
658: PetscInt *nodeList;
659: IS vsIS, stratumIS;
660: DMLabel vsLabel;
661: DMGetLabel(dm, "Vertex Sets", &vsLabel);
662: DMLabelGetValueIS(vsLabel, &vsIS);
663: ISGetIndices(vsIS, &vsIdx);
664: for (vs = 0; vs < num_vs; ++vs) {
665: DMLabelGetStratumIS(vsLabel, vsIdx[vs], &stratumIS);
666: ISGetIndices(stratumIS, &vertices);
667: ISGetSize(stratumIS, &vsSize);
668: PetscMalloc1(vsSize, &nodeList);
669: for (i = 0; i < vsSize; ++i) nodeList[i] = vertices[i] - skipCells + 1;
670: PetscCallExternal(ex_put_set_param, exo->exoid, EX_NODE_SET, vsIdx[vs], vsSize, 0);
671: PetscCallExternal(ex_put_set, exo->exoid, EX_NODE_SET, vsIdx[vs], nodeList, NULL);
672: ISRestoreIndices(stratumIS, &vertices);
673: ISDestroy(&stratumIS);
674: PetscFree(nodeList);
675: }
676: ISRestoreIndices(vsIS, &vsIdx);
677: ISDestroy(&vsIS);
678: }
679: /* --- Side Sets/Face Sets --- */
680: DMHasLabel(dm, "Face Sets", &hasLabel);
681: if (hasLabel) {
682: PetscInt i, j, fs, fsSize;
683: const PetscInt *fsIdx, *faces;
684: IS fsIS, stratumIS;
685: DMLabel fsLabel;
686: PetscInt numPoints, *points;
687: PetscInt elem_list_size = 0;
688: PetscInt *elem_list, *elem_ind, *side_list;
690: DMGetLabel(dm, "Face Sets", &fsLabel);
691: /* Compute size of Node List and Element List */
692: DMLabelGetValueIS(fsLabel, &fsIS);
693: ISGetIndices(fsIS, &fsIdx);
694: for (fs = 0; fs < num_fs; ++fs) {
695: DMLabelGetStratumIS(fsLabel, fsIdx[fs], &stratumIS);
696: ISGetSize(stratumIS, &fsSize);
697: elem_list_size += fsSize;
698: ISDestroy(&stratumIS);
699: }
700: PetscMalloc3(num_fs, &elem_ind, elem_list_size, &elem_list, elem_list_size, &side_list);
701: elem_ind[0] = 0;
702: for (fs = 0; fs < num_fs; ++fs) {
703: DMLabelGetStratumIS(fsLabel, fsIdx[fs], &stratumIS);
704: ISGetIndices(stratumIS, &faces);
705: ISGetSize(stratumIS, &fsSize);
706: /* Set Parameters */
707: PetscCallExternal(ex_put_set_param, exo->exoid, EX_SIDE_SET, fsIdx[fs], fsSize, 0);
708: /* Indices */
709: if (fs < num_fs - 1) elem_ind[fs + 1] = elem_ind[fs] + fsSize;
711: for (i = 0; i < fsSize; ++i) {
712: /* Element List */
713: points = NULL;
714: DMPlexGetTransitiveClosure(dm, faces[i], PETSC_FALSE, &numPoints, &points);
715: elem_list[elem_ind[fs] + i] = points[2] + 1;
716: DMPlexRestoreTransitiveClosure(dm, faces[i], PETSC_FALSE, &numPoints, &points);
718: /* Side List */
719: points = NULL;
720: DMPlexGetTransitiveClosure(dm, elem_list[elem_ind[fs] + i] - 1, PETSC_TRUE, &numPoints, &points);
721: for (j = 1; j < numPoints; ++j) {
722: if (points[j * 2] == faces[i]) break;
723: }
724: /* Convert HEX sides */
725: if (numPoints == 27) {
726: if (j == 1) {
727: j = 5;
728: } else if (j == 2) {
729: j = 6;
730: } else if (j == 3) {
731: j = 1;
732: } else if (j == 4) {
733: j = 3;
734: } else if (j == 5) {
735: j = 2;
736: } else if (j == 6) {
737: j = 4;
738: }
739: }
740: /* Convert TET sides */
741: if (numPoints == 15) {
742: --j;
743: if (j == 0) j = 4;
744: }
745: side_list[elem_ind[fs] + i] = j;
746: DMPlexRestoreTransitiveClosure(dm, elem_list[elem_ind[fs] + i] - 1, PETSC_TRUE, &numPoints, &points);
747: }
748: ISRestoreIndices(stratumIS, &faces);
749: ISDestroy(&stratumIS);
750: }
751: ISRestoreIndices(fsIS, &fsIdx);
752: ISDestroy(&fsIS);
754: /* Put side sets */
755: 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]]);
756: PetscFree3(elem_ind, elem_list, side_list);
757: }
758: /*
759: close the exodus file
760: */
761: ex_close(exo->exoid);
762: exo->exoid = -1;
763: }
764: PetscSectionDestroy(&coordSection);
766: /*
767: reopen the file in parallel
768: */
769: EXO_mode = EX_WRITE;
770: #if defined(PETSC_USE_64BIT_INDICES)
771: EXO_mode += EX_ALL_INT64_API;
772: #endif
773: CPU_word_size = sizeof(PetscReal);
774: IO_word_size = sizeof(PetscReal);
775: exo->exoid = ex_open_par(exo->filename, EXO_mode, &CPU_word_size, &IO_word_size, &EXO_version, comm, MPI_INFO_NULL);
777: return 0;
778: }
780: /*
781: VecView_PlexExodusII_Internal - Write a Vec corresponding in an exodus file
783: Collective on v
785: Input Parameters:
786: + v - The vector to be written
787: - viewer - The PetscViewerExodusII viewer associate to an exodus file
789: Notes:
790: The exodus result variable index is obtained by comparing the Vec name and the
791: names of variables declared in the exodus file. For instance for a Vec named "V"
792: the location in the exodus file will be the first match of "V", "V_X", "V_XX", "V_1", or "V_11"
793: amongst all variables.
794: In the event where a nodal and zonal variable both match, the function will return an error instead of
795: possibly corrupting the file
797: Level: beginner
799: .seealso: `EXOGetVarIndex_Internal()`, `DMPlexView_ExodusII()`, `VecView_PlexExodusII()`
800: @*/
801: PetscErrorCode VecView_PlexExodusII_Internal(Vec v, PetscViewer viewer)
802: {
803: DM dm;
804: MPI_Comm comm;
805: PetscMPIInt rank;
807: int exoid, offsetN = 0, offsetZ = 0;
808: const char *vecname;
809: PetscInt step;
811: PetscObjectGetComm((PetscObject)v, &comm);
812: MPI_Comm_rank(comm, &rank);
813: PetscViewerExodusIIGetId(viewer, &exoid);
814: VecGetDM(v, &dm);
815: PetscObjectGetName((PetscObject)v, &vecname);
817: DMGetOutputSequenceNumber(dm, &step, NULL);
818: EXOGetVarIndex_Internal(exoid, EX_NODAL, vecname, &offsetN);
819: EXOGetVarIndex_Internal(exoid, EX_ELEM_BLOCK, vecname, &offsetZ);
821: if (offsetN > 0) {
822: VecViewPlex_ExodusII_Nodal_Internal(v, exoid, (int)step + 1, offsetN);
823: } else if (offsetZ > 0) {
824: VecViewPlex_ExodusII_Zonal_Internal(v, exoid, (int)step + 1, offsetZ);
825: } else SETERRQ(comm, PETSC_ERR_FILE_UNEXPECTED, "Could not find nodal or zonal variable %s in exodus file. ", vecname);
826: return 0;
827: }
829: /*
830: VecLoad_PlexExodusII_Internal - Write a Vec corresponding in an exodus file
832: Collective on v
834: Input Parameters:
835: + v - The vector to be written
836: - viewer - The PetscViewerExodusII viewer associate to an exodus file
838: Notes:
839: The exodus result variable index is obtained by comparing the Vec name and the
840: names of variables declared in the exodus file. For instance for a Vec named "V"
841: the location in the exodus file will be the first match of "V", "V_X", "V_XX", "V_1", or "V_11"
842: amongst all variables.
843: In the event where a nodal and zonal variable both match, the function will return an error instead of
844: possibly corrupting the file
846: Level: beginner
848: .seealso: `EXOGetVarIndex_Internal()`, `DMPlexView_ExodusII()`, `VecView_PlexExodusII()`
849: @*/
850: PetscErrorCode VecLoad_PlexExodusII_Internal(Vec v, PetscViewer viewer)
851: {
852: DM dm;
853: MPI_Comm comm;
854: PetscMPIInt rank;
856: int exoid, offsetN = 0, offsetZ = 0;
857: const char *vecname;
858: PetscInt step;
860: PetscObjectGetComm((PetscObject)v, &comm);
861: MPI_Comm_rank(comm, &rank);
862: PetscViewerExodusIIGetId(viewer, &exoid);
863: VecGetDM(v, &dm);
864: PetscObjectGetName((PetscObject)v, &vecname);
866: DMGetOutputSequenceNumber(dm, &step, NULL);
867: EXOGetVarIndex_Internal(exoid, EX_NODAL, vecname, &offsetN);
868: EXOGetVarIndex_Internal(exoid, EX_ELEM_BLOCK, vecname, &offsetZ);
870: if (offsetN > 0) VecLoadPlex_ExodusII_Nodal_Internal(v, exoid, (int)step + 1, offsetN);
871: else if (offsetZ > 0) VecLoadPlex_ExodusII_Zonal_Internal(v, exoid, (int)step + 1, offsetZ);
872: else SETERRQ(comm, PETSC_ERR_FILE_UNEXPECTED, "Could not find nodal or zonal variable %s in exodus file. ", vecname);
873: return 0;
874: }
876: /*
877: VecViewPlex_ExodusII_Nodal_Internal - Write a Vec corresponding to a nodal field to an exodus file
879: Collective on v
881: Input Parameters:
882: + v - The vector to be written
883: . exoid - the exodus id of a file (obtained from ex_open or ex_create for instance)
884: . step - the time step to write at (exodus steps are numbered starting from 1)
885: - offset - the location of the variable in the file
887: Notes:
888: The exodus result nodal variable index is obtained by comparing the Vec name and the
889: names of nodal variables declared in the exodus file. For instance for a Vec named "V"
890: the location in the exodus file will be the first match of "V", "V_X", "V_XX", "V_1", or "V_11"
891: amongst all nodal variables.
893: Level: beginner
895: .seealso: `EXOGetVarIndex_Internal()`, `DMPlexView_ExodusII_Internal()`, `VecLoadNodal_PlexEXO()`, `VecViewZonal_PlexEXO()`, `VecLoadZonal_PlexEXO()`
896: @*/
897: PetscErrorCode VecViewPlex_ExodusII_Nodal_Internal(Vec v, int exoid, int step, int offset)
898: {
899: MPI_Comm comm;
900: PetscMPIInt size;
901: DM dm;
902: Vec vNatural, vComp;
903: const PetscScalar *varray;
904: PetscInt xs, xe, bs;
905: PetscBool useNatural;
907: PetscObjectGetComm((PetscObject)v, &comm);
908: MPI_Comm_size(comm, &size);
909: VecGetDM(v, &dm);
910: DMGetUseNatural(dm, &useNatural);
911: useNatural = useNatural && size > 1 ? PETSC_TRUE : PETSC_FALSE;
912: if (useNatural) {
913: DMPlexCreateNaturalVector(dm, &vNatural);
914: DMPlexGlobalToNaturalBegin(dm, v, vNatural);
915: DMPlexGlobalToNaturalEnd(dm, v, vNatural);
916: } else {
917: vNatural = v;
918: }
920: /* Write local chunk of the result in the exodus file
921: exodus stores each component of a vector-valued field as a separate variable.
922: We assume that they are stored sequentially */
923: VecGetOwnershipRange(vNatural, &xs, &xe);
924: VecGetBlockSize(vNatural, &bs);
925: if (bs == 1) {
926: VecGetArrayRead(vNatural, &varray);
927: PetscCallExternal(ex_put_partial_var, exoid, step, EX_NODAL, offset, 1, xs + 1, xe - xs, varray);
928: VecRestoreArrayRead(vNatural, &varray);
929: } else {
930: IS compIS;
931: PetscInt c;
933: ISCreateStride(comm, (xe - xs) / bs, xs, bs, &compIS);
934: for (c = 0; c < bs; ++c) {
935: ISStrideSetStride(compIS, (xe - xs) / bs, xs + c, bs);
936: VecGetSubVector(vNatural, compIS, &vComp);
937: VecGetArrayRead(vComp, &varray);
938: PetscCallExternal(ex_put_partial_var, exoid, step, EX_NODAL, offset + c, 1, xs / bs + 1, (xe - xs) / bs, varray);
939: VecRestoreArrayRead(vComp, &varray);
940: VecRestoreSubVector(vNatural, compIS, &vComp);
941: }
942: ISDestroy(&compIS);
943: }
944: if (useNatural) VecDestroy(&vNatural);
945: return 0;
946: }
948: /*
949: VecLoadPlex_ExodusII_Nodal_Internal - Read a Vec corresponding to a nodal field from an exodus file
951: Collective on v
953: Input Parameters:
954: + v - The vector to be written
955: . exoid - the exodus id of a file (obtained from ex_open or ex_create for instance)
956: . step - the time step to read at (exodus steps are numbered starting from 1)
957: - offset - the location of the variable in the file
959: Notes:
960: The exodus result nodal variable index is obtained by comparing the Vec name and the
961: names of nodal variables declared in the exodus file. For instance for a Vec named "V"
962: the location in the exodus file will be the first match of "V", "V_X", "V_XX", "V_1", or "V_11"
963: amongst all nodal variables.
965: Level: beginner
967: .seealso: `EXOGetVarIndex_Internal()`, `DMPlexView_ExodusII_Internal()`, `VecViewPlex_ExodusII_Nodal_Internal()`, `VecViewZonal_PlexEXO()`, `VecLoadZonal_PlexEXO()`
968: */
969: PetscErrorCode VecLoadPlex_ExodusII_Nodal_Internal(Vec v, int exoid, int step, int offset)
970: {
971: MPI_Comm comm;
972: PetscMPIInt size;
973: DM dm;
974: Vec vNatural, vComp;
975: PetscScalar *varray;
976: PetscInt xs, xe, bs;
977: PetscBool useNatural;
979: PetscObjectGetComm((PetscObject)v, &comm);
980: MPI_Comm_size(comm, &size);
981: VecGetDM(v, &dm);
982: DMGetUseNatural(dm, &useNatural);
983: useNatural = useNatural && size > 1 ? PETSC_TRUE : PETSC_FALSE;
984: if (useNatural) DMPlexCreateNaturalVector(dm, &vNatural);
985: else vNatural = v;
987: /* Read local chunk from the file */
988: VecGetOwnershipRange(vNatural, &xs, &xe);
989: VecGetBlockSize(vNatural, &bs);
990: if (bs == 1) {
991: VecGetArray(vNatural, &varray);
992: PetscCallExternal(ex_get_partial_var, exoid, step, EX_NODAL, offset, 1, xs + 1, xe - xs, varray);
993: VecRestoreArray(vNatural, &varray);
994: } else {
995: IS compIS;
996: PetscInt c;
998: ISCreateStride(comm, (xe - xs) / bs, xs, bs, &compIS);
999: for (c = 0; c < bs; ++c) {
1000: ISStrideSetStride(compIS, (xe - xs) / bs, xs + c, bs);
1001: VecGetSubVector(vNatural, compIS, &vComp);
1002: VecGetArray(vComp, &varray);
1003: PetscCallExternal(ex_get_partial_var, exoid, step, EX_NODAL, offset + c, 1, xs / bs + 1, (xe - xs) / bs, varray);
1004: VecRestoreArray(vComp, &varray);
1005: VecRestoreSubVector(vNatural, compIS, &vComp);
1006: }
1007: ISDestroy(&compIS);
1008: }
1009: if (useNatural) {
1010: DMPlexNaturalToGlobalBegin(dm, vNatural, v);
1011: DMPlexNaturalToGlobalEnd(dm, vNatural, v);
1012: VecDestroy(&vNatural);
1013: }
1014: return 0;
1015: }
1017: /*
1018: VecViewPlex_ExodusII_Zonal_Internal - Write a Vec corresponding to a zonal (cell based) field to an exodus file
1020: Collective on v
1022: Input Parameters:
1023: + v - The vector to be written
1024: . exoid - the exodus id of a file (obtained from ex_open or ex_create for instance)
1025: . step - the time step to write at (exodus steps are numbered starting from 1)
1026: - offset - the location of the variable in the file
1028: Notes:
1029: The exodus result zonal variable index is obtained by comparing the Vec name and the
1030: names of zonal variables declared in the exodus file. For instance for a Vec named "V"
1031: the location in the exodus file will be the first match of "V", "V_X", "V_XX", "V_1", or "V_11"
1032: amongst all zonal variables.
1034: Level: beginner
1036: .seealso: `EXOGetVarIndex_Internal()`, `DMPlexView_ExodusII_Internal()`, `VecViewPlex_ExodusII_Nodal_Internal()`, `VecLoadPlex_ExodusII_Nodal_Internal()`, `VecLoadPlex_ExodusII_Zonal_Internal()`
1037: */
1038: PetscErrorCode VecViewPlex_ExodusII_Zonal_Internal(Vec v, int exoid, int step, int offset)
1039: {
1040: MPI_Comm comm;
1041: PetscMPIInt size;
1042: DM dm;
1043: Vec vNatural, vComp;
1044: const PetscScalar *varray;
1045: PetscInt xs, xe, bs;
1046: PetscBool useNatural;
1047: IS compIS;
1048: PetscInt *csSize, *csID;
1049: PetscInt numCS, set, csxs = 0;
1051: PetscObjectGetComm((PetscObject)v, &comm);
1052: MPI_Comm_size(comm, &size);
1053: VecGetDM(v, &dm);
1054: DMGetUseNatural(dm, &useNatural);
1055: useNatural = useNatural && size > 1 ? PETSC_TRUE : PETSC_FALSE;
1056: if (useNatural) {
1057: DMPlexCreateNaturalVector(dm, &vNatural);
1058: DMPlexGlobalToNaturalBegin(dm, v, vNatural);
1059: DMPlexGlobalToNaturalEnd(dm, v, vNatural);
1060: } else {
1061: vNatural = v;
1062: }
1064: /* Write local chunk of the result in the exodus file
1065: exodus stores each component of a vector-valued field as a separate variable.
1066: We assume that they are stored sequentially
1067: Zonal variables are accessed one element block at a time, so we loop through the cell sets,
1068: but once the vector has been reordered to natural size, we cannot use the label information
1069: to figure out what to save where. */
1070: numCS = ex_inquire_int(exoid, EX_INQ_ELEM_BLK);
1071: PetscMalloc2(numCS, &csID, numCS, &csSize);
1072: PetscCallExternal(ex_get_ids, exoid, EX_ELEM_BLOCK, csID);
1073: for (set = 0; set < numCS; ++set) {
1074: ex_block block;
1076: block.id = csID[set];
1077: block.type = EX_ELEM_BLOCK;
1078: PetscCallExternal(ex_get_block_param, exoid, &block);
1079: csSize[set] = block.num_entry;
1080: }
1081: VecGetOwnershipRange(vNatural, &xs, &xe);
1082: VecGetBlockSize(vNatural, &bs);
1083: if (bs > 1) ISCreateStride(comm, (xe - xs) / bs, xs, bs, &compIS);
1084: for (set = 0; set < numCS; set++) {
1085: PetscInt csLocalSize, c;
1087: /* range of indices for set setID[set]: csxs:csxs + csSize[set]-1
1088: local slice of zonal values: xs/bs,xm/bs-1
1089: intersection: max(xs/bs,csxs),min(xm/bs-1,csxs + csSize[set]-1) */
1090: csLocalSize = PetscMax(0, PetscMin(xe / bs, csxs + csSize[set]) - PetscMax(xs / bs, csxs));
1091: if (bs == 1) {
1092: VecGetArrayRead(vNatural, &varray);
1093: PetscCallExternal(ex_put_partial_var, exoid, step, EX_ELEM_BLOCK, offset, csID[set], PetscMax(xs - csxs, 0) + 1, csLocalSize, &varray[PetscMax(0, csxs - xs)]);
1094: VecRestoreArrayRead(vNatural, &varray);
1095: } else {
1096: for (c = 0; c < bs; ++c) {
1097: ISStrideSetStride(compIS, (xe - xs) / bs, xs + c, bs);
1098: VecGetSubVector(vNatural, compIS, &vComp);
1099: VecGetArrayRead(vComp, &varray);
1100: 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)]);
1101: VecRestoreArrayRead(vComp, &varray);
1102: VecRestoreSubVector(vNatural, compIS, &vComp);
1103: }
1104: }
1105: csxs += csSize[set];
1106: }
1107: PetscFree2(csID, csSize);
1108: if (bs > 1) ISDestroy(&compIS);
1109: if (useNatural) VecDestroy(&vNatural);
1110: return 0;
1111: }
1113: /*
1114: VecLoadPlex_ExodusII_Zonal_Internal - Read a Vec corresponding to a zonal (cell based) field from an exodus file
1116: Collective on v
1118: Input Parameters:
1119: + v - The vector to be written
1120: . exoid - the exodus id of a file (obtained from ex_open or ex_create for instance)
1121: . step - the time step to read at (exodus steps are numbered starting from 1)
1122: - offset - the location of the variable in the file
1124: Notes:
1125: The exodus result zonal variable index is obtained by comparing the Vec name and the
1126: names of zonal variables declared in the exodus file. For instance for a Vec named "V"
1127: the location in the exodus file will be the first match of "V", "V_X", "V_XX", "V_1", or "V_11"
1128: amongst all zonal variables.
1130: Level: beginner
1132: .seealso: `EXOGetVarIndex_Internal()`, `DMPlexView_ExodusII_Internal()`, `VecViewPlex_ExodusII_Nodal_Internal()`, `VecLoadPlex_ExodusII_Nodal_Internal()`, `VecLoadPlex_ExodusII_Zonal_Internal()`
1133: */
1134: PetscErrorCode VecLoadPlex_ExodusII_Zonal_Internal(Vec v, int exoid, int step, int offset)
1135: {
1136: MPI_Comm comm;
1137: PetscMPIInt size;
1138: DM dm;
1139: Vec vNatural, vComp;
1140: PetscScalar *varray;
1141: PetscInt xs, xe, bs;
1142: PetscBool useNatural;
1143: IS compIS;
1144: PetscInt *csSize, *csID;
1145: PetscInt numCS, set, csxs = 0;
1147: PetscObjectGetComm((PetscObject)v, &comm);
1148: MPI_Comm_size(comm, &size);
1149: VecGetDM(v, &dm);
1150: DMGetUseNatural(dm, &useNatural);
1151: useNatural = useNatural && size > 1 ? PETSC_TRUE : PETSC_FALSE;
1152: if (useNatural) DMPlexCreateNaturalVector(dm, &vNatural);
1153: else vNatural = v;
1155: /* Read local chunk of the result in the exodus file
1156: exodus stores each component of a vector-valued field as a separate variable.
1157: We assume that they are stored sequentially
1158: Zonal variables are accessed one element block at a time, so we loop through the cell sets,
1159: but once the vector has been reordered to natural size, we cannot use the label information
1160: to figure out what to save where. */
1161: numCS = ex_inquire_int(exoid, EX_INQ_ELEM_BLK);
1162: PetscMalloc2(numCS, &csID, numCS, &csSize);
1163: PetscCallExternal(ex_get_ids, exoid, EX_ELEM_BLOCK, csID);
1164: for (set = 0; set < numCS; ++set) {
1165: ex_block block;
1167: block.id = csID[set];
1168: block.type = EX_ELEM_BLOCK;
1169: PetscCallExternal(ex_get_block_param, exoid, &block);
1170: csSize[set] = block.num_entry;
1171: }
1172: VecGetOwnershipRange(vNatural, &xs, &xe);
1173: VecGetBlockSize(vNatural, &bs);
1174: if (bs > 1) ISCreateStride(comm, (xe - xs) / bs, xs, bs, &compIS);
1175: for (set = 0; set < numCS; ++set) {
1176: PetscInt csLocalSize, c;
1178: /* range of indices for set setID[set]: csxs:csxs + csSize[set]-1
1179: local slice of zonal values: xs/bs,xm/bs-1
1180: intersection: max(xs/bs,csxs),min(xm/bs-1,csxs + csSize[set]-1) */
1181: csLocalSize = PetscMax(0, PetscMin(xe / bs, csxs + csSize[set]) - PetscMax(xs / bs, csxs));
1182: if (bs == 1) {
1183: VecGetArray(vNatural, &varray);
1184: PetscCallExternal(ex_get_partial_var, exoid, step, EX_ELEM_BLOCK, offset, csID[set], PetscMax(xs - csxs, 0) + 1, csLocalSize, &varray[PetscMax(0, csxs - xs)]);
1185: VecRestoreArray(vNatural, &varray);
1186: } else {
1187: for (c = 0; c < bs; ++c) {
1188: ISStrideSetStride(compIS, (xe - xs) / bs, xs + c, bs);
1189: VecGetSubVector(vNatural, compIS, &vComp);
1190: VecGetArray(vComp, &varray);
1191: 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)]);
1192: VecRestoreArray(vComp, &varray);
1193: VecRestoreSubVector(vNatural, compIS, &vComp);
1194: }
1195: }
1196: csxs += csSize[set];
1197: }
1198: PetscFree2(csID, csSize);
1199: if (bs > 1) ISDestroy(&compIS);
1200: if (useNatural) {
1201: DMPlexNaturalToGlobalBegin(dm, vNatural, v);
1202: DMPlexNaturalToGlobalEnd(dm, vNatural, v);
1203: VecDestroy(&vNatural);
1204: }
1205: return 0;
1206: }
1207: #endif
1209: /*@
1210: PetscViewerExodusIIGetId - Get the file id of the `PETSCVIEWEREXODUSII` file
1212: Logically Collective on viewer
1214: Input Parameter:
1215: . viewer - the `PetscViewer`
1217: Output Parameter:
1218: . exoid - The ExodusII file id
1220: Level: intermediate
1222: .seealso: `PETSCVIEWEREXODUSII`, `PetscViewer`, `PetscViewerFileSetMode()`, `PetscViewerCreate()`, `PetscViewerSetType()`, `PetscViewerBinaryOpen()`
1223: @*/
1224: PetscErrorCode PetscViewerExodusIIGetId(PetscViewer viewer, int *exoid)
1225: {
1227: PetscTryMethod(viewer, "PetscViewerGetId_C", (PetscViewer, int *), (viewer, exoid));
1228: return 0;
1229: }
1231: /*@
1232: PetscViewerExodusIISetOrder - Set the elements order in the exodusII file.
1234: Collective
1236: Input Parameters:
1237: + viewer - the `PETSCVIEWEREXODUSII` viewer
1238: - order - elements order
1240: Output Parameter:
1242: Level: beginner
1244: .seealso: `PETSCVIEWEREXODUSII`, `PetscViewer`, `PetscViewerExodusIIGetId()`, `PetscViewerExodusIIGetOrder()`, `PetscViewerExodusIISetOrder()`
1245: @*/
1246: PetscErrorCode PetscViewerExodusIISetOrder(PetscViewer viewer, PetscInt order)
1247: {
1249: PetscTryMethod(viewer, "PetscViewerSetOrder_C", (PetscViewer, PetscInt), (viewer, order));
1250: return 0;
1251: }
1253: /*@
1254: PetscViewerExodusIIGetOrder - Get the elements order in the exodusII file.
1256: Collective
1258: Input Parameters:
1259: + viewer - the `PETSCVIEWEREXODUSII` viewer
1260: - order - elements order
1262: Output Parameter:
1264: Level: beginner
1266: .seealso: `PETSCVIEWEREXODUSII`, `PetscViewer`, `PetscViewerExodusIIGetId()`, `PetscViewerExodusIIGetOrder()`, `PetscViewerExodusIISetOrder()`
1267: @*/
1268: PetscErrorCode PetscViewerExodusIIGetOrder(PetscViewer viewer, PetscInt *order)
1269: {
1271: PetscTryMethod(viewer, "PetscViewerGetOrder_C", (PetscViewer, PetscInt *), (viewer, order));
1272: return 0;
1273: }
1275: /*@C
1276: PetscViewerExodusIIOpen - Opens a file for ExodusII input/output.
1278: Collective
1280: Input Parameters:
1281: + comm - MPI communicator
1282: . name - name of file
1283: - type - type of file
1284: .vb
1285: FILE_MODE_WRITE - create new file for binary output
1286: FILE_MODE_READ - open existing file for binary input
1287: FILE_MODE_APPEND - open existing file for binary output
1288: .ve
1290: Output Parameter:
1291: . exo - `PETSCVIEWEREXODUSII` `PetscViewer` for Exodus II input/output to use with the specified file
1293: Level: beginner
1295: .seealso: `PETSCVIEWEREXODUSII`, `PetscViewer`, `PetscViewerPushFormat()`, `PetscViewerDestroy()`,
1296: `DMLoad()`, `PetscFileMode`, `PetscViewer`, `PetscViewerSetType()`, `PetscViewerFileSetMode()`, `PetscViewerFileSetName()`
1297: @*/
1298: PetscErrorCode PetscViewerExodusIIOpen(MPI_Comm comm, const char name[], PetscFileMode type, PetscViewer *exo)
1299: {
1300: PetscViewerCreate(comm, exo);
1301: PetscViewerSetType(*exo, PETSCVIEWEREXODUSII);
1302: PetscViewerFileSetMode(*exo, type);
1303: PetscViewerFileSetName(*exo, name);
1304: PetscViewerSetFromOptions(*exo);
1305: return 0;
1306: }
1308: /*@C
1309: DMPlexCreateExodusFromFile - Create a `DMPLEX` mesh from an ExodusII file.
1311: Collective
1313: Input Parameters:
1314: + comm - The MPI communicator
1315: . filename - The name of the ExodusII file
1316: - interpolate - Create faces and edges in the mesh
1318: Output Parameter:
1319: . dm - The `DM` object representing the mesh
1321: Level: beginner
1323: .seealso: [](chapter_unstructured), `DM`, `PETSCVIEWEREXODUSII`, `DMPLEX`, `DMCreate()`, `DMPlexCreateExodus()`
1324: @*/
1325: PetscErrorCode DMPlexCreateExodusFromFile(MPI_Comm comm, const char filename[], PetscBool interpolate, DM *dm)
1326: {
1327: PetscMPIInt rank;
1328: #if defined(PETSC_HAVE_EXODUSII)
1329: int CPU_word_size = sizeof(PetscReal), IO_word_size = 0, exoid = -1;
1330: float version;
1331: #endif
1334: MPI_Comm_rank(comm, &rank);
1335: #if defined(PETSC_HAVE_EXODUSII)
1336: if (rank == 0) {
1337: exoid = ex_open(filename, EX_READ, &CPU_word_size, &IO_word_size, &version);
1339: }
1340: DMPlexCreateExodus(comm, exoid, interpolate, dm);
1341: if (rank == 0) PetscCallExternal(ex_close, exoid);
1342: return 0;
1343: #else
1344: SETERRQ(comm, PETSC_ERR_SUP, "This method requires ExodusII support. Reconfigure using --download-exodusii");
1345: #endif
1346: }
1348: #if defined(PETSC_HAVE_EXODUSII)
1349: static PetscErrorCode ExodusGetCellType_Internal(const char *elem_type, DMPolytopeType *ct)
1350: {
1351: PetscBool flg;
1353: *ct = DM_POLYTOPE_UNKNOWN;
1354: PetscStrcmp(elem_type, "TRI", &flg);
1355: if (flg) {
1356: *ct = DM_POLYTOPE_TRIANGLE;
1357: goto done;
1358: }
1359: PetscStrcmp(elem_type, "TRI3", &flg);
1360: if (flg) {
1361: *ct = DM_POLYTOPE_TRIANGLE;
1362: goto done;
1363: }
1364: PetscStrcmp(elem_type, "QUAD", &flg);
1365: if (flg) {
1366: *ct = DM_POLYTOPE_QUADRILATERAL;
1367: goto done;
1368: }
1369: PetscStrcmp(elem_type, "QUAD4", &flg);
1370: if (flg) {
1371: *ct = DM_POLYTOPE_QUADRILATERAL;
1372: goto done;
1373: }
1374: PetscStrcmp(elem_type, "SHELL4", &flg);
1375: if (flg) {
1376: *ct = DM_POLYTOPE_QUADRILATERAL;
1377: goto done;
1378: }
1379: PetscStrcmp(elem_type, "TETRA", &flg);
1380: if (flg) {
1381: *ct = DM_POLYTOPE_TETRAHEDRON;
1382: goto done;
1383: }
1384: PetscStrcmp(elem_type, "TET4", &flg);
1385: if (flg) {
1386: *ct = DM_POLYTOPE_TETRAHEDRON;
1387: goto done;
1388: }
1389: PetscStrcmp(elem_type, "WEDGE", &flg);
1390: if (flg) {
1391: *ct = DM_POLYTOPE_TRI_PRISM;
1392: goto done;
1393: }
1394: PetscStrcmp(elem_type, "HEX", &flg);
1395: if (flg) {
1396: *ct = DM_POLYTOPE_HEXAHEDRON;
1397: goto done;
1398: }
1399: PetscStrcmp(elem_type, "HEX8", &flg);
1400: if (flg) {
1401: *ct = DM_POLYTOPE_HEXAHEDRON;
1402: goto done;
1403: }
1404: PetscStrcmp(elem_type, "HEXAHEDRON", &flg);
1405: if (flg) {
1406: *ct = DM_POLYTOPE_HEXAHEDRON;
1407: goto done;
1408: }
1409: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unrecognized element type %s", elem_type);
1410: done:
1411: return 0;
1412: }
1413: #endif
1415: /*@
1416: DMPlexCreateExodus - Create a `DMPLEX` mesh from an ExodusII file ID.
1418: Collective
1420: Input Parameters:
1421: + comm - The MPI communicator
1422: . exoid - The ExodusII id associated with a exodus file and obtained using ex_open
1423: - interpolate - Create faces and edges in the mesh
1425: Output Parameter:
1426: . dm - The `DM` object representing the mesh
1428: Level: beginner
1430: .seealso: [](chapter_unstructured), `DM`, `PETSCVIEWEREXODUSII`, `DMPLEX`, `DMPLEX`, `DMCreate()`
1431: @*/
1432: PetscErrorCode DMPlexCreateExodus(MPI_Comm comm, PetscInt exoid, PetscBool interpolate, DM *dm)
1433: {
1434: #if defined(PETSC_HAVE_EXODUSII)
1435: PetscMPIInt num_proc, rank;
1436: DMLabel cellSets = NULL, faceSets = NULL, vertSets = NULL;
1437: PetscSection coordSection;
1438: Vec coordinates;
1439: PetscScalar *coords;
1440: PetscInt coordSize, v;
1441: /* Read from ex_get_init() */
1442: char title[PETSC_MAX_PATH_LEN + 1];
1443: int dim = 0, dimEmbed = 0, numVertices = 0, numCells = 0;
1444: int num_cs = 0, num_vs = 0, num_fs = 0;
1445: #endif
1447: #if defined(PETSC_HAVE_EXODUSII)
1448: MPI_Comm_rank(comm, &rank);
1449: MPI_Comm_size(comm, &num_proc);
1450: DMCreate(comm, dm);
1451: DMSetType(*dm, DMPLEX);
1452: /* Open EXODUS II file and read basic information on rank 0, then broadcast to all processors */
1453: if (rank == 0) {
1454: PetscMemzero(title, PETSC_MAX_PATH_LEN + 1);
1455: PetscCallExternal(ex_get_init, exoid, title, &dimEmbed, &numVertices, &numCells, &num_cs, &num_vs, &num_fs);
1457: }
1458: MPI_Bcast(title, PETSC_MAX_PATH_LEN + 1, MPI_CHAR, 0, comm);
1459: MPI_Bcast(&dim, 1, MPI_INT, 0, comm);
1460: PetscObjectSetName((PetscObject)*dm, title);
1461: DMPlexSetChart(*dm, 0, numCells + numVertices);
1462: /* We do not want this label automatically computed, instead we compute it here */
1463: DMCreateLabel(*dm, "celltype");
1465: /* Read cell sets information */
1466: if (rank == 0) {
1467: PetscInt *cone;
1468: int c, cs, ncs, c_loc, v, v_loc;
1469: /* Read from ex_get_elem_blk_ids() */
1470: int *cs_id, *cs_order;
1471: /* Read from ex_get_elem_block() */
1472: char buffer[PETSC_MAX_PATH_LEN + 1];
1473: int num_cell_in_set, num_vertex_per_cell, num_hybrid, num_attr;
1474: /* Read from ex_get_elem_conn() */
1475: int *cs_connect;
1477: /* Get cell sets IDs */
1478: PetscMalloc2(num_cs, &cs_id, num_cs, &cs_order);
1479: PetscCallExternal(ex_get_ids, exoid, EX_ELEM_BLOCK, cs_id);
1480: /* Read the cell set connectivity table and build mesh topology
1481: EXO standard requires that cells in cell sets be numbered sequentially and be pairwise disjoint. */
1482: /* Check for a hybrid mesh */
1483: for (cs = 0, num_hybrid = 0; cs < num_cs; ++cs) {
1484: DMPolytopeType ct;
1485: char elem_type[PETSC_MAX_PATH_LEN];
1487: PetscArrayzero(elem_type, sizeof(elem_type));
1488: PetscCallExternal(ex_get_elem_type, exoid, cs_id[cs], elem_type);
1489: ExodusGetCellType_Internal(elem_type, &ct);
1490: dim = PetscMax(dim, DMPolytopeTypeGetDim(ct));
1491: PetscCallExternal(ex_get_block, exoid, EX_ELEM_BLOCK, cs_id[cs], buffer, &num_cell_in_set, &num_vertex_per_cell, 0, 0, &num_attr);
1492: switch (ct) {
1493: case DM_POLYTOPE_TRI_PRISM:
1494: cs_order[cs] = cs;
1495: ++num_hybrid;
1496: break;
1497: default:
1498: for (c = cs; c > cs - num_hybrid; --c) cs_order[c] = cs_order[c - 1];
1499: cs_order[cs - num_hybrid] = cs;
1500: }
1501: }
1502: /* First set sizes */
1503: for (ncs = 0, c = 0; ncs < num_cs; ++ncs) {
1504: DMPolytopeType ct;
1505: char elem_type[PETSC_MAX_PATH_LEN];
1506: const PetscInt cs = cs_order[ncs];
1508: PetscArrayzero(elem_type, sizeof(elem_type));
1509: PetscCallExternal(ex_get_elem_type, exoid, cs_id[cs], elem_type);
1510: ExodusGetCellType_Internal(elem_type, &ct);
1511: PetscCallExternal(ex_get_block, exoid, EX_ELEM_BLOCK, cs_id[cs], buffer, &num_cell_in_set, &num_vertex_per_cell, 0, 0, &num_attr);
1512: for (c_loc = 0; c_loc < num_cell_in_set; ++c_loc, ++c) {
1513: DMPlexSetConeSize(*dm, c, num_vertex_per_cell);
1514: DMPlexSetCellType(*dm, c, ct);
1515: }
1516: }
1517: for (v = numCells; v < numCells + numVertices; ++v) DMPlexSetCellType(*dm, v, DM_POLYTOPE_POINT);
1518: DMSetUp(*dm);
1519: for (ncs = 0, c = 0; ncs < num_cs; ++ncs) {
1520: const PetscInt cs = cs_order[ncs];
1521: PetscCallExternal(ex_get_block, exoid, EX_ELEM_BLOCK, cs_id[cs], buffer, &num_cell_in_set, &num_vertex_per_cell, 0, 0, &num_attr);
1522: PetscMalloc2(num_vertex_per_cell * num_cell_in_set, &cs_connect, num_vertex_per_cell, &cone);
1523: PetscCallExternal(ex_get_conn, exoid, EX_ELEM_BLOCK, cs_id[cs], cs_connect, NULL, NULL);
1524: /* EXO uses Fortran-based indexing, DMPlex uses C-style and numbers cell first then vertices. */
1525: for (c_loc = 0, v = 0; c_loc < num_cell_in_set; ++c_loc, ++c) {
1526: DMPolytopeType ct;
1528: for (v_loc = 0; v_loc < num_vertex_per_cell; ++v_loc, ++v) cone[v_loc] = cs_connect[v] + numCells - 1;
1529: DMPlexGetCellType(*dm, c, &ct);
1530: DMPlexInvertCell(ct, cone);
1531: DMPlexSetCone(*dm, c, cone);
1532: DMSetLabelValue_Fast(*dm, &cellSets, "Cell Sets", c, cs_id[cs]);
1533: }
1534: PetscFree2(cs_connect, cone);
1535: }
1536: PetscFree2(cs_id, cs_order);
1537: }
1538: {
1539: PetscInt ints[] = {dim, dimEmbed};
1541: MPI_Bcast(ints, 2, MPIU_INT, 0, comm);
1542: DMSetDimension(*dm, ints[0]);
1543: DMSetCoordinateDim(*dm, ints[1]);
1544: dim = ints[0];
1545: dimEmbed = ints[1];
1546: }
1547: DMPlexSymmetrize(*dm);
1548: DMPlexStratify(*dm);
1549: if (interpolate) {
1550: DM idm;
1552: DMPlexInterpolate(*dm, &idm);
1553: DMDestroy(dm);
1554: *dm = idm;
1555: }
1557: /* Create vertex set label */
1558: if (rank == 0 && (num_vs > 0)) {
1559: int vs, v;
1560: /* Read from ex_get_node_set_ids() */
1561: int *vs_id;
1562: /* Read from ex_get_node_set_param() */
1563: int num_vertex_in_set;
1564: /* Read from ex_get_node_set() */
1565: int *vs_vertex_list;
1567: /* Get vertex set ids */
1568: PetscMalloc1(num_vs, &vs_id);
1569: PetscCallExternal(ex_get_ids, exoid, EX_NODE_SET, vs_id);
1570: for (vs = 0; vs < num_vs; ++vs) {
1571: PetscCallExternal(ex_get_set_param, exoid, EX_NODE_SET, vs_id[vs], &num_vertex_in_set, NULL);
1572: PetscMalloc1(num_vertex_in_set, &vs_vertex_list);
1573: PetscCallExternal(ex_get_set, exoid, EX_NODE_SET, vs_id[vs], vs_vertex_list, NULL);
1574: for (v = 0; v < num_vertex_in_set; ++v) DMSetLabelValue_Fast(*dm, &vertSets, "Vertex Sets", vs_vertex_list[v] + numCells - 1, vs_id[vs]);
1575: PetscFree(vs_vertex_list);
1576: }
1577: PetscFree(vs_id);
1578: }
1579: /* Read coordinates */
1580: DMGetCoordinateSection(*dm, &coordSection);
1581: PetscSectionSetNumFields(coordSection, 1);
1582: PetscSectionSetFieldComponents(coordSection, 0, dimEmbed);
1583: PetscSectionSetChart(coordSection, numCells, numCells + numVertices);
1584: for (v = numCells; v < numCells + numVertices; ++v) {
1585: PetscSectionSetDof(coordSection, v, dimEmbed);
1586: PetscSectionSetFieldDof(coordSection, v, 0, dimEmbed);
1587: }
1588: PetscSectionSetUp(coordSection);
1589: PetscSectionGetStorageSize(coordSection, &coordSize);
1590: VecCreate(PETSC_COMM_SELF, &coordinates);
1591: PetscObjectSetName((PetscObject)coordinates, "coordinates");
1592: VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);
1593: VecSetBlockSize(coordinates, dimEmbed);
1594: VecSetType(coordinates, VECSTANDARD);
1595: VecGetArray(coordinates, &coords);
1596: if (rank == 0) {
1597: PetscReal *x, *y, *z;
1599: PetscMalloc3(numVertices, &x, numVertices, &y, numVertices, &z);
1600: PetscCallExternal(ex_get_coord, exoid, x, y, z);
1601: if (dimEmbed > 0) {
1602: for (v = 0; v < numVertices; ++v) coords[v * dimEmbed + 0] = x[v];
1603: }
1604: if (dimEmbed > 1) {
1605: for (v = 0; v < numVertices; ++v) coords[v * dimEmbed + 1] = y[v];
1606: }
1607: if (dimEmbed > 2) {
1608: for (v = 0; v < numVertices; ++v) coords[v * dimEmbed + 2] = z[v];
1609: }
1610: PetscFree3(x, y, z);
1611: }
1612: VecRestoreArray(coordinates, &coords);
1613: DMSetCoordinatesLocal(*dm, coordinates);
1614: VecDestroy(&coordinates);
1616: /* Create side set label */
1617: if (rank == 0 && interpolate && (num_fs > 0)) {
1618: int fs, f, voff;
1619: /* Read from ex_get_side_set_ids() */
1620: int *fs_id;
1621: /* Read from ex_get_side_set_param() */
1622: int num_side_in_set;
1623: /* Read from ex_get_side_set_node_list() */
1624: int *fs_vertex_count_list, *fs_vertex_list;
1625: /* Read side set labels */
1626: char fs_name[MAX_STR_LENGTH + 1];
1627: size_t fs_name_len;
1629: /* Get side set ids */
1630: PetscMalloc1(num_fs, &fs_id);
1631: PetscCallExternal(ex_get_ids, exoid, EX_SIDE_SET, fs_id);
1632: for (fs = 0; fs < num_fs; ++fs) {
1633: PetscCallExternal(ex_get_set_param, exoid, EX_SIDE_SET, fs_id[fs], &num_side_in_set, NULL);
1634: PetscMalloc2(num_side_in_set, &fs_vertex_count_list, num_side_in_set * 4, &fs_vertex_list);
1635: PetscCallExternal(ex_get_side_set_node_list, exoid, fs_id[fs], fs_vertex_count_list, fs_vertex_list);
1636: /* Get the specific name associated with this side set ID. */
1637: int fs_name_err = ex_get_name(exoid, EX_SIDE_SET, fs_id[fs], fs_name);
1638: if (!fs_name_err) {
1639: PetscStrlen(fs_name, &fs_name_len);
1640: if (fs_name_len == 0) PetscStrncpy(fs_name, "Face Sets", MAX_STR_LENGTH);
1641: }
1642: for (f = 0, voff = 0; f < num_side_in_set; ++f) {
1643: const PetscInt *faces = NULL;
1644: PetscInt faceSize = fs_vertex_count_list[f], numFaces;
1645: PetscInt faceVertices[4], v;
1648: for (v = 0; v < faceSize; ++v, ++voff) faceVertices[v] = fs_vertex_list[voff] + numCells - 1;
1649: DMPlexGetFullJoin(*dm, faceSize, faceVertices, &numFaces, &faces);
1651: DMSetLabelValue_Fast(*dm, &faceSets, "Face Sets", faces[0], fs_id[fs]);
1652: /* Only add the label if one has been detected for this side set. */
1653: if (!fs_name_err) DMSetLabelValue(*dm, fs_name, faces[0], fs_id[fs]);
1654: DMPlexRestoreJoin(*dm, faceSize, faceVertices, &numFaces, &faces);
1655: }
1656: PetscFree2(fs_vertex_count_list, fs_vertex_list);
1657: }
1658: PetscFree(fs_id);
1659: }
1661: { /* Create Cell/Face/Vertex Sets labels at all processes */
1662: enum {
1663: n = 3
1664: };
1665: PetscBool flag[n];
1667: flag[0] = cellSets ? PETSC_TRUE : PETSC_FALSE;
1668: flag[1] = faceSets ? PETSC_TRUE : PETSC_FALSE;
1669: flag[2] = vertSets ? PETSC_TRUE : PETSC_FALSE;
1670: MPI_Bcast(flag, n, MPIU_BOOL, 0, comm);
1671: if (flag[0]) DMCreateLabel(*dm, "Cell Sets");
1672: if (flag[1]) DMCreateLabel(*dm, "Face Sets");
1673: if (flag[2]) DMCreateLabel(*dm, "Vertex Sets");
1674: }
1675: return 0;
1676: #else
1677: SETERRQ(comm, PETSC_ERR_SUP, "This method requires ExodusII support. Reconfigure using --download-exodusii");
1678: #endif
1679: }