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: /*
347: Check that all points are on rank 0 since we don't know how to save distributed DM in exodus format
348: */
349: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
350: DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);
351: DMPlexGetDepthStratum(dm, 1, &eStart, &eEnd);
352: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
353: numCells = cEnd - cStart;
354: numEdges = eEnd - eStart;
355: numVertices = vEnd - vStart;
357: if (rank == 0) {
358: switch (exo->btype) {
359: case FILE_MODE_READ:
360: case FILE_MODE_APPEND:
361: case FILE_MODE_UPDATE:
362: case FILE_MODE_APPEND_UPDATE:
363: /* exodusII does not allow writing geometry to an existing file */
364: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "cannot add geometry to existing file %s", exo->filename);
365: case FILE_MODE_WRITE:
366: /* Create an empty file if one already exists*/
367: EXO_mode = EX_CLOBBER;
368: #if defined(PETSC_USE_64BIT_INDICES)
369: EXO_mode += EX_ALL_INT64_API;
370: #endif
371: CPU_word_size = sizeof(PetscReal);
372: IO_word_size = sizeof(PetscReal);
373: exo->exoid = ex_create(exo->filename, EXO_mode, &CPU_word_size, &IO_word_size);
376: break;
377: default:
378: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ORDER, "Must call PetscViewerFileSetMode() before PetscViewerFileSetName()");
379: }
381: /* --- Get DM info --- */
382: PetscObjectGetName((PetscObject)dm, &dmName);
383: DMPlexGetDepth(dm, &depth);
384: DMGetDimension(dm, &dim);
385: DMPlexGetChart(dm, &pStart, &pEnd);
386: if (depth == 3) {
387: numFaces = fEnd - fStart;
388: } else {
389: numFaces = 0;
390: }
391: DMGetLabelSize(dm, "Cell Sets", &num_cs);
392: DMGetLabelSize(dm, "Vertex Sets", &num_vs);
393: DMGetLabelSize(dm, "Face Sets", &num_fs);
394: DMGetCoordinatesLocal(dm, &coord);
395: DMGetCoordinateDM(dm, &cdm);
396: if (num_cs > 0) {
397: DMGetLabel(dm, "Cell Sets", &csLabel);
398: DMLabelGetValueIS(csLabel, &csIS);
399: ISGetIndices(csIS, &csIdx);
400: }
401: PetscMalloc1(num_cs, &nodes);
402: /* Set element type for each block and compute total number of nodes */
403: PetscMalloc1(num_cs, &type);
404: numNodes = numVertices;
406: PetscViewerExodusIIGetOrder(viewer, °ree);
407: if (degree == 2) numNodes += numEdges;
408: cellsNotInConnectivity = numCells;
409: for (cs = 0; cs < num_cs; ++cs) {
410: IS stratumIS;
411: const PetscInt *cells;
412: PetscScalar *xyz = NULL;
413: PetscInt csSize, closureSize;
415: DMLabelGetStratumIS(csLabel, csIdx[cs], &stratumIS);
416: ISGetIndices(stratumIS, &cells);
417: ISGetSize(stratumIS, &csSize);
418: DMPlexVecGetClosure(cdm, NULL, coord, cells[0], &closureSize, &xyz);
419: switch (dim) {
420: case 2:
421: if (closureSize == 3 * dim) {
422: type[cs] = TRI;
423: } else if (closureSize == 4 * dim) {
424: type[cs] = QUAD;
425: } 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);
426: break;
427: case 3:
428: if (closureSize == 4 * dim) {
429: type[cs] = TET;
430: } else if (closureSize == 8 * dim) {
431: type[cs] = HEX;
432: } 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);
433: break;
434: default:
435: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Dimension %" PetscInt_FMT " not handled by ExodusII viewer", dim);
436: }
437: if ((degree == 2) && (type[cs] == QUAD)) numNodes += csSize;
438: if ((degree == 2) && (type[cs] == HEX)) {
439: numNodes += csSize;
440: numNodes += numFaces;
441: }
442: DMPlexVecRestoreClosure(cdm, NULL, coord, cells[0], &closureSize, &xyz);
443: /* Set nodes and Element type */
444: if (type[cs] == TRI) {
445: if (degree == 1) nodes[cs] = nodesTriP1;
446: else if (degree == 2) nodes[cs] = nodesTriP2;
447: } else if (type[cs] == QUAD) {
448: if (degree == 1) nodes[cs] = nodesQuadP1;
449: else if (degree == 2) nodes[cs] = nodesQuadP2;
450: } else if (type[cs] == TET) {
451: if (degree == 1) nodes[cs] = nodesTetP1;
452: else if (degree == 2) nodes[cs] = nodesTetP2;
453: } else if (type[cs] == HEX) {
454: if (degree == 1) nodes[cs] = nodesHexP1;
455: else if (degree == 2) nodes[cs] = nodesHexP2;
456: }
457: /* Compute the number of cells not in the connectivity table */
458: cellsNotInConnectivity -= nodes[cs][3] * csSize;
460: ISRestoreIndices(stratumIS, &cells);
461: ISDestroy(&stratumIS);
462: }
463: if (num_cs) PetscCallExternal(ex_put_init, exo->exoid, dmName, dim, numNodes, numCells, num_cs, num_vs, num_fs);
464: /* --- Connectivity --- */
465: for (cs = 0; cs < num_cs; ++cs) {
466: IS stratumIS;
467: const PetscInt *cells;
468: PetscInt *connect, off = 0;
469: PetscInt edgesInClosure = 0, facesInClosure = 0, verticesInClosure = 0;
470: PetscInt csSize, c, connectSize, closureSize;
471: char *elem_type = NULL;
472: char elem_type_tri3[] = "TRI3", elem_type_quad4[] = "QUAD4";
473: char elem_type_tri6[] = "TRI6", elem_type_quad9[] = "QUAD9";
474: char elem_type_tet4[] = "TET4", elem_type_hex8[] = "HEX8";
475: char elem_type_tet10[] = "TET10", elem_type_hex27[] = "HEX27";
477: DMLabelGetStratumIS(csLabel, csIdx[cs], &stratumIS);
478: ISGetIndices(stratumIS, &cells);
479: ISGetSize(stratumIS, &csSize);
480: /* Set Element type */
481: if (type[cs] == TRI) {
482: if (degree == 1) elem_type = elem_type_tri3;
483: else if (degree == 2) elem_type = elem_type_tri6;
484: } else if (type[cs] == QUAD) {
485: if (degree == 1) elem_type = elem_type_quad4;
486: else if (degree == 2) elem_type = elem_type_quad9;
487: } else if (type[cs] == TET) {
488: if (degree == 1) elem_type = elem_type_tet4;
489: else if (degree == 2) elem_type = elem_type_tet10;
490: } else if (type[cs] == HEX) {
491: if (degree == 1) elem_type = elem_type_hex8;
492: else if (degree == 2) elem_type = elem_type_hex27;
493: }
494: connectSize = nodes[cs][0] + nodes[cs][1] + nodes[cs][2] + nodes[cs][3];
495: PetscMalloc1(PetscMax(27, connectSize) * csSize, &connect);
496: PetscCallExternal(ex_put_block, exo->exoid, EX_ELEM_BLOCK, csIdx[cs], elem_type, csSize, connectSize, 0, 0, 1);
497: /* Find number of vertices, edges, and faces in the closure */
498: verticesInClosure = nodes[cs][0];
499: if (depth > 1) {
500: if (dim == 2) {
501: DMPlexGetConeSize(dm, cells[0], &edgesInClosure);
502: } else if (dim == 3) {
503: PetscInt *closure = NULL;
505: DMPlexGetConeSize(dm, cells[0], &facesInClosure);
506: DMPlexGetTransitiveClosure(dm, cells[0], PETSC_TRUE, &closureSize, &closure);
507: edgesInClosure = closureSize - facesInClosure - 1 - verticesInClosure;
508: DMPlexRestoreTransitiveClosure(dm, cells[0], PETSC_TRUE, &closureSize, &closure);
509: }
510: }
511: /* Get connectivity for each cell */
512: for (c = 0; c < csSize; ++c) {
513: PetscInt *closure = NULL;
514: PetscInt temp, i;
516: DMPlexGetTransitiveClosure(dm, cells[c], PETSC_TRUE, &closureSize, &closure);
517: for (i = 0; i < connectSize; ++i) {
518: if (i < nodes[cs][0]) { /* Vertices */
519: connect[i + off] = closure[(i + edgesInClosure + facesInClosure + 1) * 2] + 1;
520: connect[i + off] -= cellsNotInConnectivity;
521: } else if (i < nodes[cs][0] + nodes[cs][1]) { /* Edges */
522: connect[i + off] = closure[(i - verticesInClosure + facesInClosure + 1) * 2] + 1;
523: if (nodes[cs][2] == 0) connect[i + off] -= numFaces;
524: connect[i + off] -= cellsNotInConnectivity;
525: } else if (i < nodes[cs][0] + nodes[cs][1] + nodes[cs][3]) { /* Cells */
526: connect[i + off] = closure[0] + 1;
527: connect[i + off] -= skipCells;
528: } else if (i < nodes[cs][0] + nodes[cs][1] + nodes[cs][3] + nodes[cs][2]) { /* Faces */
529: connect[i + off] = closure[(i - edgesInClosure - verticesInClosure) * 2] + 1;
530: connect[i + off] -= cellsNotInConnectivity;
531: } else {
532: connect[i + off] = -1;
533: }
534: }
535: /* Tetrahedra are inverted */
536: if (type[cs] == TET) {
537: temp = connect[0 + off];
538: connect[0 + off] = connect[1 + off];
539: connect[1 + off] = temp;
540: if (degree == 2) {
541: temp = connect[5 + off];
542: connect[5 + off] = connect[6 + off];
543: connect[6 + off] = temp;
544: temp = connect[7 + off];
545: connect[7 + off] = connect[8 + off];
546: connect[8 + off] = temp;
547: }
548: }
549: /* Hexahedra are inverted */
550: if (type[cs] == HEX) {
551: temp = connect[1 + off];
552: connect[1 + off] = connect[3 + off];
553: connect[3 + off] = temp;
554: if (degree == 2) {
555: temp = connect[8 + off];
556: connect[8 + off] = connect[11 + off];
557: connect[11 + off] = temp;
558: temp = connect[9 + off];
559: connect[9 + off] = connect[10 + off];
560: connect[10 + off] = temp;
561: temp = connect[16 + off];
562: connect[16 + off] = connect[17 + off];
563: connect[17 + off] = temp;
564: temp = connect[18 + off];
565: connect[18 + off] = connect[19 + off];
566: connect[19 + off] = temp;
568: temp = connect[12 + off];
569: connect[12 + off] = connect[16 + off];
570: connect[16 + off] = temp;
571: temp = connect[13 + off];
572: connect[13 + off] = connect[17 + off];
573: connect[17 + off] = temp;
574: temp = connect[14 + off];
575: connect[14 + off] = connect[18 + off];
576: connect[18 + off] = temp;
577: temp = connect[15 + off];
578: connect[15 + off] = connect[19 + off];
579: connect[19 + off] = temp;
581: temp = connect[23 + off];
582: connect[23 + off] = connect[26 + off];
583: connect[26 + off] = temp;
584: temp = connect[24 + off];
585: connect[24 + off] = connect[25 + off];
586: connect[25 + off] = temp;
587: temp = connect[25 + off];
588: connect[25 + off] = connect[26 + off];
589: connect[26 + off] = temp;
590: }
591: }
592: off += connectSize;
593: DMPlexRestoreTransitiveClosure(dm, cells[c], PETSC_TRUE, &closureSize, &closure);
594: }
595: PetscCallExternal(ex_put_conn, exo->exoid, EX_ELEM_BLOCK, csIdx[cs], connect, 0, 0);
596: skipCells += (nodes[cs][3] == 0) * csSize;
597: PetscFree(connect);
598: ISRestoreIndices(stratumIS, &cells);
599: ISDestroy(&stratumIS);
600: }
601: PetscFree(type);
602: /* --- Coordinates --- */
603: PetscSectionSetChart(coordSection, pStart, pEnd);
604: if (num_cs) {
605: for (d = 0; d < depth; ++d) {
606: DMPlexGetDepthStratum(dm, d, &pStart, &pEnd);
607: for (p = pStart; p < pEnd; ++p) PetscSectionSetDof(coordSection, p, nodes[0][d] > 0);
608: }
609: }
610: for (cs = 0; cs < num_cs; ++cs) {
611: IS stratumIS;
612: const PetscInt *cells;
613: PetscInt csSize, c;
615: DMLabelGetStratumIS(csLabel, csIdx[cs], &stratumIS);
616: ISGetIndices(stratumIS, &cells);
617: ISGetSize(stratumIS, &csSize);
618: for (c = 0; c < csSize; ++c) PetscSectionSetDof(coordSection, cells[c], nodes[cs][3] > 0);
619: ISRestoreIndices(stratumIS, &cells);
620: ISDestroy(&stratumIS);
621: }
622: if (num_cs) {
623: ISRestoreIndices(csIS, &csIdx);
624: ISDestroy(&csIS);
625: }
626: PetscFree(nodes);
627: PetscSectionSetUp(coordSection);
628: if (numNodes) {
629: const char *coordNames[3] = {"x", "y", "z"};
630: PetscScalar *closure, *cval;
631: PetscReal *coords;
632: PetscInt hasDof, n = 0;
634: /* There can't be more than 24 values in the closure of a point for the coord coordSection */
635: PetscCalloc3(numNodes * 3, &coords, dim, &cval, 24, &closure);
636: DMGetCoordinatesLocalNoncollective(dm, &coord);
637: DMPlexGetChart(dm, &pStart, &pEnd);
638: for (p = pStart; p < pEnd; ++p) {
639: PetscSectionGetDof(coordSection, p, &hasDof);
640: if (hasDof) {
641: PetscInt closureSize = 24, j;
643: DMPlexVecGetClosure(cdm, NULL, coord, p, &closureSize, &closure);
644: for (d = 0; d < dim; ++d) {
645: cval[d] = 0.0;
646: for (j = 0; j < closureSize / dim; j++) cval[d] += closure[j * dim + d];
647: coords[d * numNodes + n] = PetscRealPart(cval[d]) * dim / closureSize;
648: }
649: ++n;
650: }
651: }
652: PetscCallExternal(ex_put_coord, exo->exoid, &coords[0 * numNodes], &coords[1 * numNodes], &coords[2 * numNodes]);
653: PetscFree3(coords, cval, closure);
654: PetscCallExternal(ex_put_coord_names, exo->exoid, (char **)coordNames);
655: }
657: /* --- Node Sets/Vertex Sets --- */
658: DMHasLabel(dm, "Vertex Sets", &hasLabel);
659: if (hasLabel) {
660: PetscInt i, vs, vsSize;
661: const PetscInt *vsIdx, *vertices;
662: PetscInt *nodeList;
663: IS vsIS, stratumIS;
664: DMLabel vsLabel;
665: DMGetLabel(dm, "Vertex Sets", &vsLabel);
666: DMLabelGetValueIS(vsLabel, &vsIS);
667: ISGetIndices(vsIS, &vsIdx);
668: for (vs = 0; vs < num_vs; ++vs) {
669: DMLabelGetStratumIS(vsLabel, vsIdx[vs], &stratumIS);
670: ISGetIndices(stratumIS, &vertices);
671: ISGetSize(stratumIS, &vsSize);
672: PetscMalloc1(vsSize, &nodeList);
673: for (i = 0; i < vsSize; ++i) nodeList[i] = vertices[i] - skipCells + 1;
674: PetscCallExternal(ex_put_set_param, exo->exoid, EX_NODE_SET, vsIdx[vs], vsSize, 0);
675: PetscCallExternal(ex_put_set, exo->exoid, EX_NODE_SET, vsIdx[vs], nodeList, NULL);
676: ISRestoreIndices(stratumIS, &vertices);
677: ISDestroy(&stratumIS);
678: PetscFree(nodeList);
679: }
680: ISRestoreIndices(vsIS, &vsIdx);
681: ISDestroy(&vsIS);
682: }
683: /* --- Side Sets/Face Sets --- */
684: DMHasLabel(dm, "Face Sets", &hasLabel);
685: if (hasLabel) {
686: PetscInt i, j, fs, fsSize;
687: const PetscInt *fsIdx, *faces;
688: IS fsIS, stratumIS;
689: DMLabel fsLabel;
690: PetscInt numPoints, *points;
691: PetscInt elem_list_size = 0;
692: PetscInt *elem_list, *elem_ind, *side_list;
694: DMGetLabel(dm, "Face Sets", &fsLabel);
695: /* Compute size of Node List and Element List */
696: DMLabelGetValueIS(fsLabel, &fsIS);
697: ISGetIndices(fsIS, &fsIdx);
698: for (fs = 0; fs < num_fs; ++fs) {
699: DMLabelGetStratumIS(fsLabel, fsIdx[fs], &stratumIS);
700: ISGetSize(stratumIS, &fsSize);
701: elem_list_size += fsSize;
702: ISDestroy(&stratumIS);
703: }
704: if (num_fs) {
705: PetscMalloc3(num_fs, &elem_ind, elem_list_size, &elem_list, elem_list_size, &side_list);
706: elem_ind[0] = 0;
707: for (fs = 0; fs < num_fs; ++fs) {
708: DMLabelGetStratumIS(fsLabel, fsIdx[fs], &stratumIS);
709: ISGetIndices(stratumIS, &faces);
710: ISGetSize(stratumIS, &fsSize);
711: /* Set Parameters */
712: PetscCallExternal(ex_put_set_param, exo->exoid, EX_SIDE_SET, fsIdx[fs], fsSize, 0);
713: /* Indices */
714: if (fs < num_fs - 1) elem_ind[fs + 1] = elem_ind[fs] + fsSize;
716: for (i = 0; i < fsSize; ++i) {
717: /* Element List */
718: points = NULL;
719: DMPlexGetTransitiveClosure(dm, faces[i], PETSC_FALSE, &numPoints, &points);
720: elem_list[elem_ind[fs] + i] = points[2] + 1;
721: DMPlexRestoreTransitiveClosure(dm, faces[i], PETSC_FALSE, &numPoints, &points);
723: /* Side List */
724: points = NULL;
725: DMPlexGetTransitiveClosure(dm, elem_list[elem_ind[fs] + i] - 1, PETSC_TRUE, &numPoints, &points);
726: for (j = 1; j < numPoints; ++j) {
727: if (points[j * 2] == faces[i]) break;
728: }
729: /* Convert HEX sides */
730: if (numPoints == 27) {
731: if (j == 1) {
732: j = 5;
733: } else if (j == 2) {
734: j = 6;
735: } else if (j == 3) {
736: j = 1;
737: } else if (j == 4) {
738: j = 3;
739: } else if (j == 5) {
740: j = 2;
741: } else if (j == 6) {
742: j = 4;
743: }
744: }
745: /* Convert TET sides */
746: if (numPoints == 15) {
747: --j;
748: if (j == 0) j = 4;
749: }
750: side_list[elem_ind[fs] + i] = j;
751: DMPlexRestoreTransitiveClosure(dm, elem_list[elem_ind[fs] + i] - 1, PETSC_TRUE, &numPoints, &points);
752: }
753: ISRestoreIndices(stratumIS, &faces);
754: ISDestroy(&stratumIS);
755: }
756: ISRestoreIndices(fsIS, &fsIdx);
757: ISDestroy(&fsIS);
759: /* Put side sets */
760: 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]]);
761: PetscFree3(elem_ind, elem_list, side_list);
762: }
763: }
764: /*
765: close the exodus file
766: */
767: ex_close(exo->exoid);
768: exo->exoid = -1;
769: }
770: PetscSectionDestroy(&coordSection);
772: /*
773: reopen the file in parallel
774: */
775: EXO_mode = EX_WRITE;
776: #if defined(PETSC_USE_64BIT_INDICES)
777: EXO_mode += EX_ALL_INT64_API;
778: #endif
779: CPU_word_size = sizeof(PetscReal);
780: IO_word_size = sizeof(PetscReal);
781: exo->exoid = ex_open_par(exo->filename, EXO_mode, &CPU_word_size, &IO_word_size, &EXO_version, comm, MPI_INFO_NULL);
783: return 0;
784: }
786: /*
787: VecView_PlexExodusII_Internal - Write a Vec corresponding in an exodus file
789: Collective on v
791: Input Parameters:
792: + v - The vector to be written
793: - viewer - The PetscViewerExodusII viewer associate to an exodus file
795: Notes:
796: The exodus result variable index is obtained by comparing the Vec name and the
797: names of variables declared in the exodus file. For instance for a Vec named "V"
798: the location in the exodus file will be the first match of "V", "V_X", "V_XX", "V_1", or "V_11"
799: amongst all variables.
800: In the event where a nodal and zonal variable both match, the function will return an error instead of
801: possibly corrupting the file
803: Level: beginner
805: .seealso: `EXOGetVarIndex_Internal()`, `DMPlexView_ExodusII()`, `VecView_PlexExodusII()`
806: @*/
807: PetscErrorCode VecView_PlexExodusII_Internal(Vec v, PetscViewer viewer)
808: {
809: DM dm;
810: MPI_Comm comm;
811: PetscMPIInt rank;
813: int exoid, offsetN = 0, offsetZ = 0;
814: const char *vecname;
815: PetscInt step;
817: PetscObjectGetComm((PetscObject)v, &comm);
818: MPI_Comm_rank(comm, &rank);
819: PetscViewerExodusIIGetId(viewer, &exoid);
820: VecGetDM(v, &dm);
821: PetscObjectGetName((PetscObject)v, &vecname);
823: DMGetOutputSequenceNumber(dm, &step, NULL);
824: EXOGetVarIndex_Internal(exoid, EX_NODAL, vecname, &offsetN);
825: EXOGetVarIndex_Internal(exoid, EX_ELEM_BLOCK, vecname, &offsetZ);
827: if (offsetN > 0) {
828: VecViewPlex_ExodusII_Nodal_Internal(v, exoid, (int)step + 1, offsetN);
829: } else if (offsetZ > 0) {
830: VecViewPlex_ExodusII_Zonal_Internal(v, exoid, (int)step + 1, offsetZ);
831: } else SETERRQ(comm, PETSC_ERR_FILE_UNEXPECTED, "Could not find nodal or zonal variable %s in exodus file. ", vecname);
832: return 0;
833: }
835: /*
836: VecLoad_PlexExodusII_Internal - Write a Vec corresponding in an exodus file
838: Collective on v
840: Input Parameters:
841: + v - The vector to be written
842: - viewer - The PetscViewerExodusII viewer associate to an exodus file
844: Notes:
845: The exodus result variable index is obtained by comparing the Vec name and the
846: names of variables declared in the exodus file. For instance for a Vec named "V"
847: the location in the exodus file will be the first match of "V", "V_X", "V_XX", "V_1", or "V_11"
848: amongst all variables.
849: In the event where a nodal and zonal variable both match, the function will return an error instead of
850: possibly corrupting the file
852: Level: beginner
854: .seealso: `EXOGetVarIndex_Internal()`, `DMPlexView_ExodusII()`, `VecView_PlexExodusII()`
855: @*/
856: PetscErrorCode VecLoad_PlexExodusII_Internal(Vec v, PetscViewer viewer)
857: {
858: DM dm;
859: MPI_Comm comm;
860: PetscMPIInt rank;
862: int exoid, offsetN = 0, offsetZ = 0;
863: const char *vecname;
864: PetscInt step;
866: PetscObjectGetComm((PetscObject)v, &comm);
867: MPI_Comm_rank(comm, &rank);
868: PetscViewerExodusIIGetId(viewer, &exoid);
869: VecGetDM(v, &dm);
870: PetscObjectGetName((PetscObject)v, &vecname);
872: DMGetOutputSequenceNumber(dm, &step, NULL);
873: EXOGetVarIndex_Internal(exoid, EX_NODAL, vecname, &offsetN);
874: EXOGetVarIndex_Internal(exoid, EX_ELEM_BLOCK, vecname, &offsetZ);
876: if (offsetN > 0) VecLoadPlex_ExodusII_Nodal_Internal(v, exoid, (int)step + 1, offsetN);
877: else if (offsetZ > 0) VecLoadPlex_ExodusII_Zonal_Internal(v, exoid, (int)step + 1, offsetZ);
878: else SETERRQ(comm, PETSC_ERR_FILE_UNEXPECTED, "Could not find nodal or zonal variable %s in exodus file. ", vecname);
879: return 0;
880: }
882: /*
883: VecViewPlex_ExodusII_Nodal_Internal - Write a Vec corresponding to a nodal field to an exodus file
885: Collective on v
887: Input Parameters:
888: + v - The vector to be written
889: . exoid - the exodus id of a file (obtained from ex_open or ex_create for instance)
890: . step - the time step to write at (exodus steps are numbered starting from 1)
891: - offset - the location of the variable in the file
893: Notes:
894: The exodus result nodal variable index is obtained by comparing the Vec name and the
895: names of nodal variables declared in the exodus file. For instance for a Vec named "V"
896: the location in the exodus file will be the first match of "V", "V_X", "V_XX", "V_1", or "V_11"
897: amongst all nodal variables.
899: Level: beginner
901: .seealso: `EXOGetVarIndex_Internal()`, `DMPlexView_ExodusII_Internal()`, `VecLoadNodal_PlexEXO()`, `VecViewZonal_PlexEXO()`, `VecLoadZonal_PlexEXO()`
902: @*/
903: PetscErrorCode VecViewPlex_ExodusII_Nodal_Internal(Vec v, int exoid, int step, int offset)
904: {
905: MPI_Comm comm;
906: PetscMPIInt size;
907: DM dm;
908: Vec vNatural, vComp;
909: const PetscScalar *varray;
910: PetscInt xs, xe, bs;
911: PetscBool useNatural;
913: PetscObjectGetComm((PetscObject)v, &comm);
914: MPI_Comm_size(comm, &size);
915: VecGetDM(v, &dm);
916: DMGetUseNatural(dm, &useNatural);
917: useNatural = useNatural && size > 1 ? PETSC_TRUE : PETSC_FALSE;
918: if (useNatural) {
919: DMPlexCreateNaturalVector(dm, &vNatural);
920: DMPlexGlobalToNaturalBegin(dm, v, vNatural);
921: DMPlexGlobalToNaturalEnd(dm, v, vNatural);
922: } else {
923: vNatural = v;
924: }
926: /* Write local chunk of the result in the exodus file
927: exodus stores each component of a vector-valued field as a separate variable.
928: We assume that they are stored sequentially */
929: VecGetOwnershipRange(vNatural, &xs, &xe);
930: VecGetBlockSize(vNatural, &bs);
931: if (bs == 1) {
932: VecGetArrayRead(vNatural, &varray);
933: PetscCallExternal(ex_put_partial_var, exoid, step, EX_NODAL, offset, 1, xs + 1, xe - xs, varray);
934: VecRestoreArrayRead(vNatural, &varray);
935: } else {
936: IS compIS;
937: PetscInt c;
939: ISCreateStride(comm, (xe - xs) / bs, xs, bs, &compIS);
940: for (c = 0; c < bs; ++c) {
941: ISStrideSetStride(compIS, (xe - xs) / bs, xs + c, bs);
942: VecGetSubVector(vNatural, compIS, &vComp);
943: VecGetArrayRead(vComp, &varray);
944: PetscCallExternal(ex_put_partial_var, exoid, step, EX_NODAL, offset + c, 1, xs / bs + 1, (xe - xs) / bs, varray);
945: VecRestoreArrayRead(vComp, &varray);
946: VecRestoreSubVector(vNatural, compIS, &vComp);
947: }
948: ISDestroy(&compIS);
949: }
950: if (useNatural) VecDestroy(&vNatural);
951: return 0;
952: }
954: /*
955: VecLoadPlex_ExodusII_Nodal_Internal - Read a Vec corresponding to a nodal field from an exodus file
957: Collective on v
959: Input Parameters:
960: + v - The vector to be written
961: . exoid - the exodus id of a file (obtained from ex_open or ex_create for instance)
962: . step - the time step to read at (exodus steps are numbered starting from 1)
963: - offset - the location of the variable in the file
965: Notes:
966: The exodus result nodal variable index is obtained by comparing the Vec name and the
967: names of nodal variables declared in the exodus file. For instance for a Vec named "V"
968: the location in the exodus file will be the first match of "V", "V_X", "V_XX", "V_1", or "V_11"
969: amongst all nodal variables.
971: Level: beginner
973: .seealso: `EXOGetVarIndex_Internal()`, `DMPlexView_ExodusII_Internal()`, `VecViewPlex_ExodusII_Nodal_Internal()`, `VecViewZonal_PlexEXO()`, `VecLoadZonal_PlexEXO()`
974: */
975: PetscErrorCode VecLoadPlex_ExodusII_Nodal_Internal(Vec v, int exoid, int step, int offset)
976: {
977: MPI_Comm comm;
978: PetscMPIInt size;
979: DM dm;
980: Vec vNatural, vComp;
981: PetscScalar *varray;
982: PetscInt xs, xe, bs;
983: PetscBool useNatural;
985: PetscObjectGetComm((PetscObject)v, &comm);
986: MPI_Comm_size(comm, &size);
987: VecGetDM(v, &dm);
988: DMGetUseNatural(dm, &useNatural);
989: useNatural = useNatural && size > 1 ? PETSC_TRUE : PETSC_FALSE;
990: if (useNatural) DMPlexCreateNaturalVector(dm, &vNatural);
991: else vNatural = v;
993: /* Read local chunk from the file */
994: VecGetOwnershipRange(vNatural, &xs, &xe);
995: VecGetBlockSize(vNatural, &bs);
996: if (bs == 1) {
997: VecGetArray(vNatural, &varray);
998: PetscCallExternal(ex_get_partial_var, exoid, step, EX_NODAL, offset, 1, xs + 1, xe - xs, varray);
999: VecRestoreArray(vNatural, &varray);
1000: } else {
1001: IS compIS;
1002: PetscInt c;
1004: ISCreateStride(comm, (xe - xs) / bs, xs, bs, &compIS);
1005: for (c = 0; c < bs; ++c) {
1006: ISStrideSetStride(compIS, (xe - xs) / bs, xs + c, bs);
1007: VecGetSubVector(vNatural, compIS, &vComp);
1008: VecGetArray(vComp, &varray);
1009: PetscCallExternal(ex_get_partial_var, exoid, step, EX_NODAL, offset + c, 1, xs / bs + 1, (xe - xs) / bs, varray);
1010: VecRestoreArray(vComp, &varray);
1011: VecRestoreSubVector(vNatural, compIS, &vComp);
1012: }
1013: ISDestroy(&compIS);
1014: }
1015: if (useNatural) {
1016: DMPlexNaturalToGlobalBegin(dm, vNatural, v);
1017: DMPlexNaturalToGlobalEnd(dm, vNatural, v);
1018: VecDestroy(&vNatural);
1019: }
1020: return 0;
1021: }
1023: /*
1024: VecViewPlex_ExodusII_Zonal_Internal - Write a Vec corresponding to a zonal (cell based) field to an exodus file
1026: Collective on v
1028: Input Parameters:
1029: + v - The vector to be written
1030: . exoid - the exodus id of a file (obtained from ex_open or ex_create for instance)
1031: . step - the time step to write at (exodus steps are numbered starting from 1)
1032: - offset - the location of the variable in the file
1034: Notes:
1035: The exodus result zonal variable index is obtained by comparing the Vec name and the
1036: names of zonal variables declared in the exodus file. For instance for a Vec named "V"
1037: the location in the exodus file will be the first match of "V", "V_X", "V_XX", "V_1", or "V_11"
1038: amongst all zonal variables.
1040: Level: beginner
1042: .seealso: `EXOGetVarIndex_Internal()`, `DMPlexView_ExodusII_Internal()`, `VecViewPlex_ExodusII_Nodal_Internal()`, `VecLoadPlex_ExodusII_Nodal_Internal()`, `VecLoadPlex_ExodusII_Zonal_Internal()`
1043: */
1044: PetscErrorCode VecViewPlex_ExodusII_Zonal_Internal(Vec v, int exoid, int step, int offset)
1045: {
1046: MPI_Comm comm;
1047: PetscMPIInt size;
1048: DM dm;
1049: Vec vNatural, vComp;
1050: const PetscScalar *varray;
1051: PetscInt xs, xe, bs;
1052: PetscBool useNatural;
1053: IS compIS;
1054: PetscInt *csSize, *csID;
1055: PetscInt numCS, set, csxs = 0;
1057: PetscObjectGetComm((PetscObject)v, &comm);
1058: MPI_Comm_size(comm, &size);
1059: VecGetDM(v, &dm);
1060: DMGetUseNatural(dm, &useNatural);
1061: useNatural = useNatural && size > 1 ? PETSC_TRUE : PETSC_FALSE;
1062: if (useNatural) {
1063: DMPlexCreateNaturalVector(dm, &vNatural);
1064: DMPlexGlobalToNaturalBegin(dm, v, vNatural);
1065: DMPlexGlobalToNaturalEnd(dm, v, vNatural);
1066: } else {
1067: vNatural = v;
1068: }
1070: /* Write local chunk of the result in the exodus file
1071: exodus stores each component of a vector-valued field as a separate variable.
1072: We assume that they are stored sequentially
1073: Zonal variables are accessed one element block at a time, so we loop through the cell sets,
1074: but once the vector has been reordered to natural size, we cannot use the label information
1075: to figure out what to save where. */
1076: numCS = ex_inquire_int(exoid, EX_INQ_ELEM_BLK);
1077: PetscMalloc2(numCS, &csID, numCS, &csSize);
1078: PetscCallExternal(ex_get_ids, exoid, EX_ELEM_BLOCK, csID);
1079: for (set = 0; set < numCS; ++set) {
1080: ex_block block;
1082: block.id = csID[set];
1083: block.type = EX_ELEM_BLOCK;
1084: PetscCallExternal(ex_get_block_param, exoid, &block);
1085: csSize[set] = block.num_entry;
1086: }
1087: VecGetOwnershipRange(vNatural, &xs, &xe);
1088: VecGetBlockSize(vNatural, &bs);
1089: if (bs > 1) ISCreateStride(comm, (xe - xs) / bs, xs, bs, &compIS);
1090: for (set = 0; set < numCS; set++) {
1091: PetscInt csLocalSize, c;
1093: /* range of indices for set setID[set]: csxs:csxs + csSize[set]-1
1094: local slice of zonal values: xs/bs,xm/bs-1
1095: intersection: max(xs/bs,csxs),min(xm/bs-1,csxs + csSize[set]-1) */
1096: csLocalSize = PetscMax(0, PetscMin(xe / bs, csxs + csSize[set]) - PetscMax(xs / bs, csxs));
1097: if (bs == 1) {
1098: VecGetArrayRead(vNatural, &varray);
1099: PetscCallExternal(ex_put_partial_var, exoid, step, EX_ELEM_BLOCK, offset, csID[set], PetscMax(xs - csxs, 0) + 1, csLocalSize, &varray[PetscMax(0, csxs - xs)]);
1100: VecRestoreArrayRead(vNatural, &varray);
1101: } else {
1102: for (c = 0; c < bs; ++c) {
1103: ISStrideSetStride(compIS, (xe - xs) / bs, xs + c, bs);
1104: VecGetSubVector(vNatural, compIS, &vComp);
1105: VecGetArrayRead(vComp, &varray);
1106: 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)]);
1107: VecRestoreArrayRead(vComp, &varray);
1108: VecRestoreSubVector(vNatural, compIS, &vComp);
1109: }
1110: }
1111: csxs += csSize[set];
1112: }
1113: PetscFree2(csID, csSize);
1114: if (bs > 1) ISDestroy(&compIS);
1115: if (useNatural) VecDestroy(&vNatural);
1116: return 0;
1117: }
1119: /*
1120: VecLoadPlex_ExodusII_Zonal_Internal - Read a Vec corresponding to a zonal (cell based) field from an exodus file
1122: Collective on v
1124: Input Parameters:
1125: + v - The vector to be written
1126: . exoid - the exodus id of a file (obtained from ex_open or ex_create for instance)
1127: . step - the time step to read at (exodus steps are numbered starting from 1)
1128: - offset - the location of the variable in the file
1130: Notes:
1131: The exodus result zonal variable index is obtained by comparing the Vec name and the
1132: names of zonal variables declared in the exodus file. For instance for a Vec named "V"
1133: the location in the exodus file will be the first match of "V", "V_X", "V_XX", "V_1", or "V_11"
1134: amongst all zonal variables.
1136: Level: beginner
1138: .seealso: `EXOGetVarIndex_Internal()`, `DMPlexView_ExodusII_Internal()`, `VecViewPlex_ExodusII_Nodal_Internal()`, `VecLoadPlex_ExodusII_Nodal_Internal()`, `VecLoadPlex_ExodusII_Zonal_Internal()`
1139: */
1140: PetscErrorCode VecLoadPlex_ExodusII_Zonal_Internal(Vec v, int exoid, int step, int offset)
1141: {
1142: MPI_Comm comm;
1143: PetscMPIInt size;
1144: DM dm;
1145: Vec vNatural, vComp;
1146: PetscScalar *varray;
1147: PetscInt xs, xe, bs;
1148: PetscBool useNatural;
1149: IS compIS;
1150: PetscInt *csSize, *csID;
1151: PetscInt numCS, set, csxs = 0;
1153: PetscObjectGetComm((PetscObject)v, &comm);
1154: MPI_Comm_size(comm, &size);
1155: VecGetDM(v, &dm);
1156: DMGetUseNatural(dm, &useNatural);
1157: useNatural = useNatural && size > 1 ? PETSC_TRUE : PETSC_FALSE;
1158: if (useNatural) DMPlexCreateNaturalVector(dm, &vNatural);
1159: else vNatural = v;
1161: /* Read local chunk of the result in the exodus file
1162: exodus stores each component of a vector-valued field as a separate variable.
1163: We assume that they are stored sequentially
1164: Zonal variables are accessed one element block at a time, so we loop through the cell sets,
1165: but once the vector has been reordered to natural size, we cannot use the label information
1166: to figure out what to save where. */
1167: numCS = ex_inquire_int(exoid, EX_INQ_ELEM_BLK);
1168: PetscMalloc2(numCS, &csID, numCS, &csSize);
1169: PetscCallExternal(ex_get_ids, exoid, EX_ELEM_BLOCK, csID);
1170: for (set = 0; set < numCS; ++set) {
1171: ex_block block;
1173: block.id = csID[set];
1174: block.type = EX_ELEM_BLOCK;
1175: PetscCallExternal(ex_get_block_param, exoid, &block);
1176: csSize[set] = block.num_entry;
1177: }
1178: VecGetOwnershipRange(vNatural, &xs, &xe);
1179: VecGetBlockSize(vNatural, &bs);
1180: if (bs > 1) ISCreateStride(comm, (xe - xs) / bs, xs, bs, &compIS);
1181: for (set = 0; set < numCS; ++set) {
1182: PetscInt csLocalSize, c;
1184: /* range of indices for set setID[set]: csxs:csxs + csSize[set]-1
1185: local slice of zonal values: xs/bs,xm/bs-1
1186: intersection: max(xs/bs,csxs),min(xm/bs-1,csxs + csSize[set]-1) */
1187: csLocalSize = PetscMax(0, PetscMin(xe / bs, csxs + csSize[set]) - PetscMax(xs / bs, csxs));
1188: if (bs == 1) {
1189: VecGetArray(vNatural, &varray);
1190: PetscCallExternal(ex_get_partial_var, exoid, step, EX_ELEM_BLOCK, offset, csID[set], PetscMax(xs - csxs, 0) + 1, csLocalSize, &varray[PetscMax(0, csxs - xs)]);
1191: VecRestoreArray(vNatural, &varray);
1192: } else {
1193: for (c = 0; c < bs; ++c) {
1194: ISStrideSetStride(compIS, (xe - xs) / bs, xs + c, bs);
1195: VecGetSubVector(vNatural, compIS, &vComp);
1196: VecGetArray(vComp, &varray);
1197: 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)]);
1198: VecRestoreArray(vComp, &varray);
1199: VecRestoreSubVector(vNatural, compIS, &vComp);
1200: }
1201: }
1202: csxs += csSize[set];
1203: }
1204: PetscFree2(csID, csSize);
1205: if (bs > 1) ISDestroy(&compIS);
1206: if (useNatural) {
1207: DMPlexNaturalToGlobalBegin(dm, vNatural, v);
1208: DMPlexNaturalToGlobalEnd(dm, vNatural, v);
1209: VecDestroy(&vNatural);
1210: }
1211: return 0;
1212: }
1213: #endif
1215: /*@
1216: PetscViewerExodusIIGetId - Get the file id of the `PETSCVIEWEREXODUSII` file
1218: Logically Collective on viewer
1220: Input Parameter:
1221: . viewer - the `PetscViewer`
1223: Output Parameter:
1224: . exoid - The ExodusII file id
1226: Level: intermediate
1228: .seealso: `PETSCVIEWEREXODUSII`, `PetscViewer`, `PetscViewerFileSetMode()`, `PetscViewerCreate()`, `PetscViewerSetType()`, `PetscViewerBinaryOpen()`
1229: @*/
1230: PetscErrorCode PetscViewerExodusIIGetId(PetscViewer viewer, int *exoid)
1231: {
1233: PetscTryMethod(viewer, "PetscViewerGetId_C", (PetscViewer, int *), (viewer, exoid));
1234: return 0;
1235: }
1237: /*@
1238: PetscViewerExodusIISetOrder - Set the elements order in the exodusII file.
1240: Collective
1242: Input Parameters:
1243: + viewer - the `PETSCVIEWEREXODUSII` viewer
1244: - order - elements order
1246: Output Parameter:
1248: Level: beginner
1250: .seealso: `PETSCVIEWEREXODUSII`, `PetscViewer`, `PetscViewerExodusIIGetId()`, `PetscViewerExodusIIGetOrder()`, `PetscViewerExodusIISetOrder()`
1251: @*/
1252: PetscErrorCode PetscViewerExodusIISetOrder(PetscViewer viewer, PetscInt order)
1253: {
1255: PetscTryMethod(viewer, "PetscViewerSetOrder_C", (PetscViewer, PetscInt), (viewer, order));
1256: return 0;
1257: }
1259: /*@
1260: PetscViewerExodusIIGetOrder - Get the elements order in the exodusII file.
1262: Collective
1264: Input Parameters:
1265: + viewer - the `PETSCVIEWEREXODUSII` viewer
1266: - order - elements order
1268: Output Parameter:
1270: Level: beginner
1272: .seealso: `PETSCVIEWEREXODUSII`, `PetscViewer`, `PetscViewerExodusIIGetId()`, `PetscViewerExodusIIGetOrder()`, `PetscViewerExodusIISetOrder()`
1273: @*/
1274: PetscErrorCode PetscViewerExodusIIGetOrder(PetscViewer viewer, PetscInt *order)
1275: {
1277: PetscTryMethod(viewer, "PetscViewerGetOrder_C", (PetscViewer, PetscInt *), (viewer, order));
1278: return 0;
1279: }
1281: /*@C
1282: PetscViewerExodusIIOpen - Opens a file for ExodusII input/output.
1284: Collective
1286: Input Parameters:
1287: + comm - MPI communicator
1288: . name - name of file
1289: - type - type of file
1290: .vb
1291: FILE_MODE_WRITE - create new file for binary output
1292: FILE_MODE_READ - open existing file for binary input
1293: FILE_MODE_APPEND - open existing file for binary output
1294: .ve
1296: Output Parameter:
1297: . exo - `PETSCVIEWEREXODUSII` `PetscViewer` for Exodus II input/output to use with the specified file
1299: Level: beginner
1301: .seealso: `PETSCVIEWEREXODUSII`, `PetscViewer`, `PetscViewerPushFormat()`, `PetscViewerDestroy()`,
1302: `DMLoad()`, `PetscFileMode`, `PetscViewer`, `PetscViewerSetType()`, `PetscViewerFileSetMode()`, `PetscViewerFileSetName()`
1303: @*/
1304: PetscErrorCode PetscViewerExodusIIOpen(MPI_Comm comm, const char name[], PetscFileMode type, PetscViewer *exo)
1305: {
1306: PetscViewerCreate(comm, exo);
1307: PetscViewerSetType(*exo, PETSCVIEWEREXODUSII);
1308: PetscViewerFileSetMode(*exo, type);
1309: PetscViewerFileSetName(*exo, name);
1310: PetscViewerSetFromOptions(*exo);
1311: return 0;
1312: }
1314: /*@C
1315: DMPlexCreateExodusFromFile - Create a `DMPLEX` mesh from an ExodusII file.
1317: Collective
1319: Input Parameters:
1320: + comm - The MPI communicator
1321: . filename - The name of the ExodusII file
1322: - interpolate - Create faces and edges in the mesh
1324: Output Parameter:
1325: . dm - The `DM` object representing the mesh
1327: Level: beginner
1329: .seealso: [](chapter_unstructured), `DM`, `PETSCVIEWEREXODUSII`, `DMPLEX`, `DMCreate()`, `DMPlexCreateExodus()`
1330: @*/
1331: PetscErrorCode DMPlexCreateExodusFromFile(MPI_Comm comm, const char filename[], PetscBool interpolate, DM *dm)
1332: {
1333: PetscMPIInt rank;
1334: #if defined(PETSC_HAVE_EXODUSII)
1335: int CPU_word_size = sizeof(PetscReal), IO_word_size = 0, exoid = -1;
1336: float version;
1337: #endif
1340: MPI_Comm_rank(comm, &rank);
1341: #if defined(PETSC_HAVE_EXODUSII)
1342: if (rank == 0) {
1343: exoid = ex_open(filename, EX_READ, &CPU_word_size, &IO_word_size, &version);
1345: }
1346: DMPlexCreateExodus(comm, exoid, interpolate, dm);
1347: if (rank == 0) PetscCallExternal(ex_close, exoid);
1348: return 0;
1349: #else
1350: SETERRQ(comm, PETSC_ERR_SUP, "This method requires ExodusII support. Reconfigure using --download-exodusii");
1351: #endif
1352: }
1354: #if defined(PETSC_HAVE_EXODUSII)
1355: static PetscErrorCode ExodusGetCellType_Internal(const char *elem_type, DMPolytopeType *ct)
1356: {
1357: PetscBool flg;
1359: *ct = DM_POLYTOPE_UNKNOWN;
1360: PetscStrcmp(elem_type, "TRI", &flg);
1361: if (flg) {
1362: *ct = DM_POLYTOPE_TRIANGLE;
1363: goto done;
1364: }
1365: PetscStrcmp(elem_type, "TRI3", &flg);
1366: if (flg) {
1367: *ct = DM_POLYTOPE_TRIANGLE;
1368: goto done;
1369: }
1370: PetscStrcmp(elem_type, "QUAD", &flg);
1371: if (flg) {
1372: *ct = DM_POLYTOPE_QUADRILATERAL;
1373: goto done;
1374: }
1375: PetscStrcmp(elem_type, "QUAD4", &flg);
1376: if (flg) {
1377: *ct = DM_POLYTOPE_QUADRILATERAL;
1378: goto done;
1379: }
1380: PetscStrcmp(elem_type, "SHELL4", &flg);
1381: if (flg) {
1382: *ct = DM_POLYTOPE_QUADRILATERAL;
1383: goto done;
1384: }
1385: PetscStrcmp(elem_type, "TETRA", &flg);
1386: if (flg) {
1387: *ct = DM_POLYTOPE_TETRAHEDRON;
1388: goto done;
1389: }
1390: PetscStrcmp(elem_type, "TET4", &flg);
1391: if (flg) {
1392: *ct = DM_POLYTOPE_TETRAHEDRON;
1393: goto done;
1394: }
1395: PetscStrcmp(elem_type, "WEDGE", &flg);
1396: if (flg) {
1397: *ct = DM_POLYTOPE_TRI_PRISM;
1398: goto done;
1399: }
1400: PetscStrcmp(elem_type, "HEX", &flg);
1401: if (flg) {
1402: *ct = DM_POLYTOPE_HEXAHEDRON;
1403: goto done;
1404: }
1405: PetscStrcmp(elem_type, "HEX8", &flg);
1406: if (flg) {
1407: *ct = DM_POLYTOPE_HEXAHEDRON;
1408: goto done;
1409: }
1410: PetscStrcmp(elem_type, "HEXAHEDRON", &flg);
1411: if (flg) {
1412: *ct = DM_POLYTOPE_HEXAHEDRON;
1413: goto done;
1414: }
1415: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unrecognized element type %s", elem_type);
1416: done:
1417: return 0;
1418: }
1419: #endif
1421: /*@
1422: DMPlexCreateExodus - Create a `DMPLEX` mesh from an ExodusII file ID.
1424: Collective
1426: Input Parameters:
1427: + comm - The MPI communicator
1428: . exoid - The ExodusII id associated with a exodus file and obtained using ex_open
1429: - interpolate - Create faces and edges in the mesh
1431: Output Parameter:
1432: . dm - The `DM` object representing the mesh
1434: Level: beginner
1436: .seealso: [](chapter_unstructured), `DM`, `PETSCVIEWEREXODUSII`, `DMPLEX`, `DMPLEX`, `DMCreate()`
1437: @*/
1438: PetscErrorCode DMPlexCreateExodus(MPI_Comm comm, PetscInt exoid, PetscBool interpolate, DM *dm)
1439: {
1440: #if defined(PETSC_HAVE_EXODUSII)
1441: PetscMPIInt num_proc, rank;
1442: DMLabel cellSets = NULL, faceSets = NULL, vertSets = NULL;
1443: PetscSection coordSection;
1444: Vec coordinates;
1445: PetscScalar *coords;
1446: PetscInt coordSize, v;
1447: /* Read from ex_get_init() */
1448: char title[PETSC_MAX_PATH_LEN + 1];
1449: int dim = 0, dimEmbed = 0, numVertices = 0, numCells = 0;
1450: int num_cs = 0, num_vs = 0, num_fs = 0;
1451: #endif
1453: #if defined(PETSC_HAVE_EXODUSII)
1454: MPI_Comm_rank(comm, &rank);
1455: MPI_Comm_size(comm, &num_proc);
1456: DMCreate(comm, dm);
1457: DMSetType(*dm, DMPLEX);
1458: /* Open EXODUS II file and read basic information on rank 0, then broadcast to all processors */
1459: if (rank == 0) {
1460: PetscMemzero(title, PETSC_MAX_PATH_LEN + 1);
1461: PetscCallExternal(ex_get_init, exoid, title, &dimEmbed, &numVertices, &numCells, &num_cs, &num_vs, &num_fs);
1463: }
1464: MPI_Bcast(title, PETSC_MAX_PATH_LEN + 1, MPI_CHAR, 0, comm);
1465: MPI_Bcast(&dim, 1, MPI_INT, 0, comm);
1466: PetscObjectSetName((PetscObject)*dm, title);
1467: DMPlexSetChart(*dm, 0, numCells + numVertices);
1468: /* We do not want this label automatically computed, instead we compute it here */
1469: DMCreateLabel(*dm, "celltype");
1471: /* Read cell sets information */
1472: if (rank == 0) {
1473: PetscInt *cone;
1474: int c, cs, ncs, c_loc, v, v_loc;
1475: /* Read from ex_get_elem_blk_ids() */
1476: int *cs_id, *cs_order;
1477: /* Read from ex_get_elem_block() */
1478: char buffer[PETSC_MAX_PATH_LEN + 1];
1479: int num_cell_in_set, num_vertex_per_cell, num_hybrid, num_attr;
1480: /* Read from ex_get_elem_conn() */
1481: int *cs_connect;
1483: /* Get cell sets IDs */
1484: PetscMalloc2(num_cs, &cs_id, num_cs, &cs_order);
1485: PetscCallExternal(ex_get_ids, exoid, EX_ELEM_BLOCK, cs_id);
1486: /* Read the cell set connectivity table and build mesh topology
1487: EXO standard requires that cells in cell sets be numbered sequentially and be pairwise disjoint. */
1488: /* Check for a hybrid mesh */
1489: for (cs = 0, num_hybrid = 0; cs < num_cs; ++cs) {
1490: DMPolytopeType ct;
1491: char elem_type[PETSC_MAX_PATH_LEN];
1493: PetscArrayzero(elem_type, sizeof(elem_type));
1494: PetscCallExternal(ex_get_elem_type, exoid, cs_id[cs], elem_type);
1495: ExodusGetCellType_Internal(elem_type, &ct);
1496: dim = PetscMax(dim, DMPolytopeTypeGetDim(ct));
1497: PetscCallExternal(ex_get_block, exoid, EX_ELEM_BLOCK, cs_id[cs], buffer, &num_cell_in_set, &num_vertex_per_cell, 0, 0, &num_attr);
1498: switch (ct) {
1499: case DM_POLYTOPE_TRI_PRISM:
1500: cs_order[cs] = cs;
1501: ++num_hybrid;
1502: break;
1503: default:
1504: for (c = cs; c > cs - num_hybrid; --c) cs_order[c] = cs_order[c - 1];
1505: cs_order[cs - num_hybrid] = cs;
1506: }
1507: }
1508: /* First set sizes */
1509: for (ncs = 0, c = 0; ncs < num_cs; ++ncs) {
1510: DMPolytopeType ct;
1511: char elem_type[PETSC_MAX_PATH_LEN];
1512: const PetscInt cs = cs_order[ncs];
1514: PetscArrayzero(elem_type, sizeof(elem_type));
1515: PetscCallExternal(ex_get_elem_type, exoid, cs_id[cs], elem_type);
1516: ExodusGetCellType_Internal(elem_type, &ct);
1517: PetscCallExternal(ex_get_block, exoid, EX_ELEM_BLOCK, cs_id[cs], buffer, &num_cell_in_set, &num_vertex_per_cell, 0, 0, &num_attr);
1518: for (c_loc = 0; c_loc < num_cell_in_set; ++c_loc, ++c) {
1519: DMPlexSetConeSize(*dm, c, num_vertex_per_cell);
1520: DMPlexSetCellType(*dm, c, ct);
1521: }
1522: }
1523: for (v = numCells; v < numCells + numVertices; ++v) DMPlexSetCellType(*dm, v, DM_POLYTOPE_POINT);
1524: DMSetUp(*dm);
1525: for (ncs = 0, c = 0; ncs < num_cs; ++ncs) {
1526: const PetscInt cs = cs_order[ncs];
1527: PetscCallExternal(ex_get_block, exoid, EX_ELEM_BLOCK, cs_id[cs], buffer, &num_cell_in_set, &num_vertex_per_cell, 0, 0, &num_attr);
1528: PetscMalloc2(num_vertex_per_cell * num_cell_in_set, &cs_connect, num_vertex_per_cell, &cone);
1529: PetscCallExternal(ex_get_conn, exoid, EX_ELEM_BLOCK, cs_id[cs], cs_connect, NULL, NULL);
1530: /* EXO uses Fortran-based indexing, DMPlex uses C-style and numbers cell first then vertices. */
1531: for (c_loc = 0, v = 0; c_loc < num_cell_in_set; ++c_loc, ++c) {
1532: DMPolytopeType ct;
1534: for (v_loc = 0; v_loc < num_vertex_per_cell; ++v_loc, ++v) cone[v_loc] = cs_connect[v] + numCells - 1;
1535: DMPlexGetCellType(*dm, c, &ct);
1536: DMPlexInvertCell(ct, cone);
1537: DMPlexSetCone(*dm, c, cone);
1538: DMSetLabelValue_Fast(*dm, &cellSets, "Cell Sets", c, cs_id[cs]);
1539: }
1540: PetscFree2(cs_connect, cone);
1541: }
1542: PetscFree2(cs_id, cs_order);
1543: }
1544: {
1545: PetscInt ints[] = {dim, dimEmbed};
1547: MPI_Bcast(ints, 2, MPIU_INT, 0, comm);
1548: DMSetDimension(*dm, ints[0]);
1549: DMSetCoordinateDim(*dm, ints[1]);
1550: dim = ints[0];
1551: dimEmbed = ints[1];
1552: }
1553: DMPlexSymmetrize(*dm);
1554: DMPlexStratify(*dm);
1555: if (interpolate) {
1556: DM idm;
1558: DMPlexInterpolate(*dm, &idm);
1559: DMDestroy(dm);
1560: *dm = idm;
1561: }
1563: /* Create vertex set label */
1564: if (rank == 0 && (num_vs > 0)) {
1565: int vs, v;
1566: /* Read from ex_get_node_set_ids() */
1567: int *vs_id;
1568: /* Read from ex_get_node_set_param() */
1569: int num_vertex_in_set;
1570: /* Read from ex_get_node_set() */
1571: int *vs_vertex_list;
1573: /* Get vertex set ids */
1574: PetscMalloc1(num_vs, &vs_id);
1575: PetscCallExternal(ex_get_ids, exoid, EX_NODE_SET, vs_id);
1576: for (vs = 0; vs < num_vs; ++vs) {
1577: PetscCallExternal(ex_get_set_param, exoid, EX_NODE_SET, vs_id[vs], &num_vertex_in_set, NULL);
1578: PetscMalloc1(num_vertex_in_set, &vs_vertex_list);
1579: PetscCallExternal(ex_get_set, exoid, EX_NODE_SET, vs_id[vs], vs_vertex_list, NULL);
1580: for (v = 0; v < num_vertex_in_set; ++v) DMSetLabelValue_Fast(*dm, &vertSets, "Vertex Sets", vs_vertex_list[v] + numCells - 1, vs_id[vs]);
1581: PetscFree(vs_vertex_list);
1582: }
1583: PetscFree(vs_id);
1584: }
1585: /* Read coordinates */
1586: DMGetCoordinateSection(*dm, &coordSection);
1587: PetscSectionSetNumFields(coordSection, 1);
1588: PetscSectionSetFieldComponents(coordSection, 0, dimEmbed);
1589: PetscSectionSetChart(coordSection, numCells, numCells + numVertices);
1590: for (v = numCells; v < numCells + numVertices; ++v) {
1591: PetscSectionSetDof(coordSection, v, dimEmbed);
1592: PetscSectionSetFieldDof(coordSection, v, 0, dimEmbed);
1593: }
1594: PetscSectionSetUp(coordSection);
1595: PetscSectionGetStorageSize(coordSection, &coordSize);
1596: VecCreate(PETSC_COMM_SELF, &coordinates);
1597: PetscObjectSetName((PetscObject)coordinates, "coordinates");
1598: VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);
1599: VecSetBlockSize(coordinates, dimEmbed);
1600: VecSetType(coordinates, VECSTANDARD);
1601: VecGetArray(coordinates, &coords);
1602: if (rank == 0) {
1603: PetscReal *x, *y, *z;
1605: PetscMalloc3(numVertices, &x, numVertices, &y, numVertices, &z);
1606: PetscCallExternal(ex_get_coord, exoid, x, y, z);
1607: if (dimEmbed > 0) {
1608: for (v = 0; v < numVertices; ++v) coords[v * dimEmbed + 0] = x[v];
1609: }
1610: if (dimEmbed > 1) {
1611: for (v = 0; v < numVertices; ++v) coords[v * dimEmbed + 1] = y[v];
1612: }
1613: if (dimEmbed > 2) {
1614: for (v = 0; v < numVertices; ++v) coords[v * dimEmbed + 2] = z[v];
1615: }
1616: PetscFree3(x, y, z);
1617: }
1618: VecRestoreArray(coordinates, &coords);
1619: DMSetCoordinatesLocal(*dm, coordinates);
1620: VecDestroy(&coordinates);
1622: /* Create side set label */
1623: if (rank == 0 && interpolate && (num_fs > 0)) {
1624: int fs, f, voff;
1625: /* Read from ex_get_side_set_ids() */
1626: int *fs_id;
1627: /* Read from ex_get_side_set_param() */
1628: int num_side_in_set;
1629: /* Read from ex_get_side_set_node_list() */
1630: int *fs_vertex_count_list, *fs_vertex_list;
1631: /* Read side set labels */
1632: char fs_name[MAX_STR_LENGTH + 1];
1633: size_t fs_name_len;
1635: /* Get side set ids */
1636: PetscMalloc1(num_fs, &fs_id);
1637: PetscCallExternal(ex_get_ids, exoid, EX_SIDE_SET, fs_id);
1638: for (fs = 0; fs < num_fs; ++fs) {
1639: PetscCallExternal(ex_get_set_param, exoid, EX_SIDE_SET, fs_id[fs], &num_side_in_set, NULL);
1640: PetscMalloc2(num_side_in_set, &fs_vertex_count_list, num_side_in_set * 4, &fs_vertex_list);
1641: PetscCallExternal(ex_get_side_set_node_list, exoid, fs_id[fs], fs_vertex_count_list, fs_vertex_list);
1642: /* Get the specific name associated with this side set ID. */
1643: int fs_name_err = ex_get_name(exoid, EX_SIDE_SET, fs_id[fs], fs_name);
1644: if (!fs_name_err) {
1645: PetscStrlen(fs_name, &fs_name_len);
1646: if (fs_name_len == 0) PetscStrncpy(fs_name, "Face Sets", MAX_STR_LENGTH);
1647: }
1648: for (f = 0, voff = 0; f < num_side_in_set; ++f) {
1649: const PetscInt *faces = NULL;
1650: PetscInt faceSize = fs_vertex_count_list[f], numFaces;
1651: PetscInt faceVertices[4], v;
1654: for (v = 0; v < faceSize; ++v, ++voff) faceVertices[v] = fs_vertex_list[voff] + numCells - 1;
1655: DMPlexGetFullJoin(*dm, faceSize, faceVertices, &numFaces, &faces);
1657: DMSetLabelValue_Fast(*dm, &faceSets, "Face Sets", faces[0], fs_id[fs]);
1658: /* Only add the label if one has been detected for this side set. */
1659: if (!fs_name_err) DMSetLabelValue(*dm, fs_name, faces[0], fs_id[fs]);
1660: DMPlexRestoreJoin(*dm, faceSize, faceVertices, &numFaces, &faces);
1661: }
1662: PetscFree2(fs_vertex_count_list, fs_vertex_list);
1663: }
1664: PetscFree(fs_id);
1665: }
1667: { /* Create Cell/Face/Vertex Sets labels at all processes */
1668: enum {
1669: n = 3
1670: };
1671: PetscBool flag[n];
1673: flag[0] = cellSets ? PETSC_TRUE : PETSC_FALSE;
1674: flag[1] = faceSets ? PETSC_TRUE : PETSC_FALSE;
1675: flag[2] = vertSets ? PETSC_TRUE : PETSC_FALSE;
1676: MPI_Bcast(flag, n, MPIU_BOOL, 0, comm);
1677: if (flag[0]) DMCreateLabel(*dm, "Cell Sets");
1678: if (flag[1]) DMCreateLabel(*dm, "Face Sets");
1679: if (flag[2]) DMCreateLabel(*dm, "Vertex Sets");
1680: }
1681: return 0;
1682: #else
1683: SETERRQ(comm, PETSC_ERR_SUP, "This method requires ExodusII support. Reconfigure using --download-exodusii");
1684: #endif
1685: }