Actual source code: plexexodusii2.c
1: #include "petscsystypes.h"
2: #define PETSCDM_DLL
3: #include <petsc/private/dmpleximpl.h>
5: #include <netcdf.h>
6: #include <exodusII.h>
8: #include <petsc/private/viewerimpl.h>
9: #include <petsc/private/viewerexodusiiimpl.h>
10: /*@C
11: PETSC_VIEWER_EXODUSII_ - Creates an `PETSCVIEWEREXODUSII` `PetscViewer` shared by all processors in a communicator.
13: Collective; No Fortran Support
15: Input Parameter:
16: . comm - the MPI communicator to share the `PETSCVIEWEREXODUSII` `PetscViewer`
18: Level: intermediate
20: Note:
21: Unlike almost all other PETSc routines, `PETSC_VIEWER_EXODUSII_()` does not return
22: an error code. The GLVIS PetscViewer is usually used in the form
23: $ XXXView(XXX object, PETSC_VIEWER_EXODUSII_(comm));
25: .seealso: `PETSCVIEWEREXODUSII`, `PetscViewer`, `PetscViewer`, `PetscViewerExodusIIOpen()`, `PetscViewerType`, `PetscViewerCreate()`, `PetscViewerDestroy()`
26: @*/
27: PetscViewer PETSC_VIEWER_EXODUSII_(MPI_Comm comm)
28: {
29: PetscViewer viewer;
31: PetscFunctionBegin;
32: PetscCallNull(PetscViewerExodusIIOpen(comm, "mesh.exo", FILE_MODE_WRITE, &viewer));
33: PetscCallNull(PetscObjectRegisterDestroy((PetscObject)viewer));
34: PetscFunctionReturn(viewer);
35: }
37: static PetscErrorCode PetscViewerView_ExodusII(PetscViewer v, PetscViewer viewer)
38: {
39: PetscViewer_ExodusII *exo = (PetscViewer_ExodusII *)v->data;
41: PetscFunctionBegin;
42: if (exo->filename) PetscCall(PetscViewerASCIIPrintf(viewer, "Filename: %s\n", exo->filename));
43: if (exo->exoid) PetscCall(PetscViewerASCIIPrintf(viewer, "exoid: %" PetscExodusIIInt_FMT "\n", exo->exoid));
44: if (exo->btype) PetscCall(PetscViewerASCIIPrintf(viewer, "IO Mode: %d\n", exo->btype));
45: if (exo->order) PetscCall(PetscViewerASCIIPrintf(viewer, "Mesh order: %" PetscInt_FMT "\n", exo->order));
46: PetscCall(PetscViewerASCIIPrintf(viewer, "Number of nodal variables: %" PetscExodusIIInt_FMT "\n", exo->numNodalVariables));
47: for (int i = 0; i < exo->numNodalVariables; i++) PetscCall(PetscViewerASCIIPrintf(viewer, " %d: %s\n", i, exo->nodalVariableNames[i]));
48: PetscCall(PetscViewerASCIIPrintf(viewer, "Number of zonal variables: %" PetscExodusIIInt_FMT "\n", exo->numZonalVariables));
49: for (int i = 0; i < exo->numZonalVariables; i++) PetscCall(PetscViewerASCIIPrintf(viewer, " %d: %s\n", i, exo->zonalVariableNames[i]));
50: PetscFunctionReturn(PETSC_SUCCESS);
51: }
53: static PetscErrorCode PetscViewerFlush_ExodusII(PetscViewer v)
54: {
55: PetscViewer_ExodusII *exo = (PetscViewer_ExodusII *)v->data;
57: PetscFunctionBegin;
58: if (exo->exoid >= 0) PetscCallExternal(ex_update, exo->exoid);
59: PetscFunctionReturn(PETSC_SUCCESS);
60: }
62: static PetscErrorCode PetscViewerSetFromOptions_ExodusII(PetscViewer v, PetscOptionItems *PetscOptionsObject)
63: {
64: PetscFunctionBegin;
65: PetscOptionsHeadBegin(PetscOptionsObject, "ExodusII PetscViewer Options");
66: PetscOptionsHeadEnd();
67: PetscFunctionReturn(PETSC_SUCCESS);
68: }
70: static PetscErrorCode PetscViewerSetUp_ExodusII(PetscViewer viewer)
71: {
72: PetscFunctionBegin;
73: PetscFunctionReturn(PETSC_SUCCESS);
74: }
76: static PetscErrorCode PetscViewerDestroy_ExodusII(PetscViewer viewer)
77: {
78: PetscViewer_ExodusII *exo = (PetscViewer_ExodusII *)viewer->data;
80: PetscFunctionBegin;
81: if (exo->exoid >= 0) PetscCallExternal(ex_close, exo->exoid);
82: for (PetscInt i = 0; i < exo->numZonalVariables; i++) PetscCall(PetscFree(exo->zonalVariableNames[i]));
83: PetscCall(PetscFree(exo->zonalVariableNames));
84: for (PetscInt i = 0; i < exo->numNodalVariables; i++) PetscCall(PetscFree(exo->nodalVariableNames[i]));
85: PetscCall(PetscFree(exo->nodalVariableNames));
86: PetscCall(PetscFree(exo->filename));
87: PetscCall(PetscFree(exo));
88: PetscCall(PetscObjectComposeFunction((PetscObject)viewer, "PetscViewerFileSetName_C", NULL));
89: PetscCall(PetscObjectComposeFunction((PetscObject)viewer, "PetscViewerFileGetName_C", NULL));
90: PetscCall(PetscObjectComposeFunction((PetscObject)viewer, "PetscViewerFileSetMode_C", NULL));
91: PetscCall(PetscObjectComposeFunction((PetscObject)viewer, "PetscViewerFileGetMode_C", NULL));
92: PetscCall(PetscObjectComposeFunction((PetscObject)viewer, "PetscViewerGetId_C", NULL));
93: PetscCall(PetscObjectComposeFunction((PetscObject)viewer, "PetscViewerGetOrder_C", NULL));
94: PetscCall(PetscObjectComposeFunction((PetscObject)viewer, "PetscViewerSetOrder_C", NULL));
95: PetscFunctionReturn(PETSC_SUCCESS);
96: }
98: static PetscErrorCode PetscViewerFileSetName_ExodusII(PetscViewer viewer, const char name[])
99: {
100: PetscViewer_ExodusII *exo = (PetscViewer_ExodusII *)viewer->data;
101: PetscMPIInt rank;
102: PetscExodusIIInt CPU_word_size, IO_word_size, EXO_mode;
103: MPI_Info mpi_info = MPI_INFO_NULL;
104: PetscExodusIIFloat EXO_version;
106: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)viewer), &rank));
107: CPU_word_size = sizeof(PetscReal);
108: IO_word_size = sizeof(PetscReal);
110: PetscFunctionBegin;
111: if (exo->exoid >= 0) {
112: PetscCallExternal(ex_close, exo->exoid);
113: exo->exoid = -1;
114: }
115: if (exo->filename) PetscCall(PetscFree(exo->filename));
116: PetscCall(PetscStrallocpy(name, &exo->filename));
117: switch (exo->btype) {
118: case FILE_MODE_READ:
119: EXO_mode = EX_READ;
120: break;
121: case FILE_MODE_APPEND:
122: case FILE_MODE_UPDATE:
123: case FILE_MODE_APPEND_UPDATE:
124: /* Will fail if the file does not already exist */
125: EXO_mode = EX_WRITE;
126: break;
127: case FILE_MODE_WRITE:
128: /*
129: exodus only allows writing geometry upon file creation, so we will let DMView create the file.
130: */
131: PetscFunctionReturn(PETSC_SUCCESS);
132: break;
133: default:
134: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ORDER, "Must call PetscViewerFileSetMode() before PetscViewerFileSetName()");
135: }
136: #if defined(PETSC_USE_64BIT_INDICES)
137: EXO_mode += EX_ALL_INT64_API;
138: #endif
139: exo->exoid = ex_open_par(name, EXO_mode, &CPU_word_size, &IO_word_size, &EXO_version, PETSC_COMM_WORLD, mpi_info);
140: PetscCheck(exo->exoid >= 0, PETSC_COMM_SELF, PETSC_ERR_LIB, "ex_open_par failed for %s", name);
141: PetscFunctionReturn(PETSC_SUCCESS);
142: }
144: static PetscErrorCode PetscViewerFileGetName_ExodusII(PetscViewer viewer, const char *name[])
145: {
146: PetscViewer_ExodusII *exo = (PetscViewer_ExodusII *)viewer->data;
148: PetscFunctionBegin;
149: *name = exo->filename;
150: PetscFunctionReturn(PETSC_SUCCESS);
151: }
153: static PetscErrorCode PetscViewerFileSetMode_ExodusII(PetscViewer viewer, PetscFileMode type)
154: {
155: PetscViewer_ExodusII *exo = (PetscViewer_ExodusII *)viewer->data;
157: PetscFunctionBegin;
158: exo->btype = type;
159: PetscFunctionReturn(PETSC_SUCCESS);
160: }
162: static PetscErrorCode PetscViewerFileGetMode_ExodusII(PetscViewer viewer, PetscFileMode *type)
163: {
164: PetscViewer_ExodusII *exo = (PetscViewer_ExodusII *)viewer->data;
166: PetscFunctionBegin;
167: *type = exo->btype;
168: PetscFunctionReturn(PETSC_SUCCESS);
169: }
171: static PetscErrorCode PetscViewerExodusIIGetId_ExodusII(PetscViewer viewer, PetscExodusIIInt *exoid)
172: {
173: PetscViewer_ExodusII *exo = (PetscViewer_ExodusII *)viewer->data;
175: PetscFunctionBegin;
176: *exoid = exo->exoid;
177: PetscFunctionReturn(PETSC_SUCCESS);
178: }
180: static PetscErrorCode PetscViewerExodusIIGetOrder_ExodusII(PetscViewer viewer, PetscInt *order)
181: {
182: PetscViewer_ExodusII *exo = (PetscViewer_ExodusII *)viewer->data;
184: PetscFunctionBegin;
185: *order = exo->order;
186: PetscFunctionReturn(PETSC_SUCCESS);
187: }
189: static PetscErrorCode PetscViewerExodusIISetOrder_ExodusII(PetscViewer viewer, PetscInt order)
190: {
191: PetscViewer_ExodusII *exo = (PetscViewer_ExodusII *)viewer->data;
193: PetscFunctionBegin;
194: exo->order = order;
195: PetscFunctionReturn(PETSC_SUCCESS);
196: }
198: /*@
199: PetscViewerExodusIISetZonalVariable - Sets the number of zonal variables in an exodusII file
201: Collective;
203: Input Parameters:
204: + viewer - a `PetscViewer` of type `PETSCVIEWEREXODUSII`
205: - num - the number of zonal variables in the exodusII file
207: Level: intermediate
209: Notes:
210: The exodusII API does not allow changing the number of variables in a file so this function will return an error
211: if called twice, called on a read-only file, or called on file for which the number of variables has already been specified
213: .seealso: `PETSCVIEWEREXODUSII`, `PetscViewer`, `PetscViewerCreate()`, `PetscViewerDestroy()`, `PetscViewerExodusIIOpen()`, `PetscViewerSetType()`, `PetscViewerType`, `PetscViewerExodusIIGetZonalVariable()`
214: @*/
215: PetscErrorCode PetscViewerExodusIISetZonalVariable(PetscViewer viewer, PetscExodusIIInt num)
216: {
217: PetscViewer_ExodusII *exo = (PetscViewer_ExodusII *)viewer->data;
218: MPI_Comm comm;
219: PetscExodusIIInt exoid = -1;
221: PetscFunctionBegin;
222: PetscCall(PetscObjectGetComm((PetscObject)viewer, &comm));
223: PetscCheck(exo->numZonalVariables == -1, comm, PETSC_ERR_SUP, "The number of zonal variables has already been set to %d and cannot be overwritten", exo->numZonalVariables);
224: PetscCheck((exo->btype != FILE_MODE_READ) && (exo->btype != FILE_MODE_UNDEFINED), comm, PETSC_ERR_FILE_WRITE, "Cannot set the number of variables because the file is not writable");
226: exo->numZonalVariables = num;
227: PetscCall(PetscMalloc1(num, &exo->zonalVariableNames));
228: for (int i = 0; i < num; i++) { exo->zonalVariableNames[i] = NULL; }
229: PetscCall(PetscViewerExodusIIGetId(viewer, &exoid));
230: PetscCallExternal(ex_put_variable_param, exoid, EX_ELEM_BLOCK, num);
231: PetscFunctionReturn(PETSC_SUCCESS);
232: }
234: /*@
235: PetscViewerExodusIISetNodalVariable - Sets the number of nodal variables in an exodusII file
237: Collective;
239: Input Parameters:
240: + viewer - a `PetscViewer` of type `PETSCVIEWEREXODUSII`
241: - num - the number of nodal variables in the exodusII file
243: Level: intermediate
245: Notes:
246: The exodusII API does not allow changing the number of variables in a file so this function will return an error
247: if called twice, called on a read-only file, or called on file for which the number of variables has already been specified
249: .seealso: `PETSCVIEWEREXODUSII`, `PetscViewer`, `PetscViewerCreate()`, `PetscViewerDestroy()`, `PetscViewerExodusIIOpen()`, `PetscViewerSetType()`, `PetscViewerType`, `PetscViewerExodusIIGetNodalVariable()`
250: @*/
251: PetscErrorCode PetscViewerExodusIISetNodalVariable(PetscViewer viewer, PetscExodusIIInt num)
252: {
253: PetscViewer_ExodusII *exo = (PetscViewer_ExodusII *)viewer->data;
254: MPI_Comm comm;
255: PetscExodusIIInt exoid = -1;
257: PetscFunctionBegin;
258: PetscCall(PetscObjectGetComm((PetscObject)viewer, &comm));
259: PetscCheck(exo->numNodalVariables == -1, comm, PETSC_ERR_SUP, "The number of nodal variables has already been set to %d and cannot be overwritten", exo->numNodalVariables);
260: PetscCheck((exo->btype != FILE_MODE_READ) && (exo->btype != FILE_MODE_UNDEFINED), comm, PETSC_ERR_FILE_WRITE, "Cannot set the number of variables because the file is not writable");
262: exo->numNodalVariables = num;
263: PetscCall(PetscMalloc1(num, &exo->nodalVariableNames));
264: for (int i = 0; i < num; i++) { exo->nodalVariableNames[i] = NULL; }
265: PetscCall(PetscViewerExodusIIGetId(viewer, &exoid));
266: PetscCallExternal(ex_put_variable_param, exoid, EX_NODAL, num);
267: PetscFunctionReturn(PETSC_SUCCESS);
268: }
270: /*@
271: PetscViewerExodusIIGetZonalVariable - Gets the number of zonal variables in an exodusII file
273: Collective
275: Input Parameters:
276: . viewer - a `PetscViewer` of type `PETSCVIEWEREXODUSII`
278: Output Parameters:
279: . num - the number variables in the exodusII file
281: Level: intermediate
283: Notes:
284: The number of variables in the exodusII file is cached in the viewer
286: .seealso: `PETSCVIEWEREXODUSII`, `PetscViewer`, `PetscViewerCreate()`, `PetscViewerDestroy()`, `PetscViewerExodusIIOpen()`, `PetscViewerSetType()`, `PetscViewerType`, `PetscViewerExodusIIsetZonalVariable()`
287: @*/
288: PetscErrorCode PetscViewerExodusIIGetZonalVariable(PetscViewer viewer, PetscExodusIIInt *num)
289: {
290: PetscViewer_ExodusII *exo = (PetscViewer_ExodusII *)viewer->data;
291: MPI_Comm comm;
292: PetscExodusIIInt exoid = -1;
294: PetscFunctionBegin;
295: if (exo->numZonalVariables > -1) {
296: *num = exo->numZonalVariables;
297: } else {
298: PetscCall(PetscViewerExodusIIGetId(viewer, &exoid));
299: PetscCall(PetscObjectGetComm((PetscObject)viewer, &comm));
300: PetscCheck(exoid > 0, comm, PETSC_ERR_FILE_OPEN, "Exodus file is not open");
301: PetscCallExternal(ex_get_variable_param, exoid, EX_ELEM_BLOCK, num);
302: exo->numZonalVariables = *num;
303: PetscCall(PetscMalloc1(*num, &exo->zonalVariableNames));
304: for (int i = 0; i < *num; i++) { exo->zonalVariableNames[i] = NULL; }
305: }
306: PetscFunctionReturn(PETSC_SUCCESS);
307: }
309: /*@
310: PetscViewerExodusIIGetNodalVariable - Gets the number of nodal variables in an exodusII file
312: Collective
314: Input Parameters:
315: . viewer - a `PetscViewer` of type `PETSCVIEWEREXODUSII`
317: Output Parameters:
318: . num - the number variables in the exodusII file
320: Level: intermediate
322: Notes:
323: This function gets the number of nodal variables and saves it in the address of num.
325: .seealso: `PETSCVIEWEREXODUSII`, `PetscViewer`, `PetscViewerCreate()`, `PetscViewerDestroy()`, `PetscViewerExodusIIOpen()`, `PetscViewerSetType()`, `PetscViewerType`, `PetscViewerExodusIISetNodalVariable()`
326: @*/
327: PetscErrorCode PetscViewerExodusIIGetNodalVariable(PetscViewer viewer, PetscExodusIIInt *num)
328: {
329: PetscViewer_ExodusII *exo = (PetscViewer_ExodusII *)viewer->data;
330: MPI_Comm comm;
331: PetscExodusIIInt exoid = -1;
333: PetscFunctionBegin;
334: if (exo->numNodalVariables > -1) {
335: *num = exo->numNodalVariables;
336: } else {
337: PetscCall(PetscViewerExodusIIGetId(viewer, &exoid));
338: PetscCall(PetscObjectGetComm((PetscObject)viewer, &comm));
339: PetscCheck(exoid > 0, comm, PETSC_ERR_FILE_OPEN, "Exodus file is not open");
340: PetscCallExternal(ex_get_variable_param, exoid, EX_NODAL, num);
341: exo->numNodalVariables = *num;
342: PetscCall(PetscMalloc1(*num, &exo->nodalVariableNames));
343: for (int i = 0; i < *num; i++) { exo->nodalVariableNames[i] = NULL; }
344: }
345: PetscFunctionReturn(PETSC_SUCCESS);
346: }
348: /*@
349: PetscViewerExodusIISetZonalVariableName - Sets the name of a zonal variable.
351: Collective;
353: Input Parameters:
354: + viewer - a `PetscViewer` of type `PETSCVIEWEREXODUSII`
355: . idx - the index for which you want to save the name
356: - name - string containing the name characters
358: Level: intermediate
360: .seealso: `PETSCVIEWEREXODUSII`, `PetscViewer`, `PetscViewerCreate()`, `PetscViewerDestroy()`, `PetscViewerExodusIIOpen()`, `PetscViewerSetType()`, `PetscViewerType`, `PetscViewerExodusIIGetZonalVariableName()`
361: @*/
362: PetscErrorCode PetscViewerExodusIISetZonalVariableName(PetscViewer viewer, PetscExodusIIInt idx, const char name[])
363: {
364: PetscViewer_ExodusII *exo = (PetscViewer_ExodusII *)viewer->data;
366: PetscFunctionBegin;
367: PetscCheck((idx >= 0) && (idx < exo->numZonalVariables), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Variable index out of range. Was PetscViewerExodusIISetZonalVariable called?");
368: PetscCall(PetscStrallocpy(name, (char **)&exo->zonalVariableNames[idx]));
369: PetscCallExternal(ex_put_variable_name, exo->exoid, EX_ELEM_BLOCK, idx + 1, name);
370: PetscFunctionReturn(PETSC_SUCCESS);
371: }
373: /*@
374: PetscViewerExodusIISetNodalVariableName - Sets the name of a nodal variable.
376: Collective;
378: Input Parameters:
379: + viewer - a `PetscViewer` of type `PETSCVIEWEREXODUSII`
380: . idx - the index for which you want to save the name
381: - name - string containing the name characters
383: Level: intermediate
385: .seealso: `PETSCVIEWEREXODUSII`, `PetscViewer`, `PetscViewerCreate()`, `PetscViewerDestroy()`, `PetscViewerExodusIIOpen()`, `PetscViewerSetType()`, `PetscViewerType`, `PetscViewerExodusIIGetNodalVariableName()`
386: @*/
387: PetscErrorCode PetscViewerExodusIISetNodalVariableName(PetscViewer viewer, PetscExodusIIInt idx, const char name[])
388: {
389: PetscViewer_ExodusII *exo = (PetscViewer_ExodusII *)viewer->data;
391: PetscFunctionBegin;
392: PetscCheck((idx >= 0) && (idx < exo->numNodalVariables), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Variable index out of range. Was PetscViewerExodusIISetNodalVariable called?");
393: PetscCall(PetscStrallocpy(name, (char **)&exo->nodalVariableNames[idx]));
394: PetscCallExternal(ex_put_variable_name, exo->exoid, EX_NODAL, idx + 1, name);
395: PetscFunctionReturn(PETSC_SUCCESS);
396: }
398: /*@
399: PetscViewerExodusIIGetZonalVariableName - Gets the name of a zonal variable.
401: Collective;
403: Input Parameters:
404: + viewer - a `PetscViewer` of type `PETSCVIEWEREXODUSII`
405: - idx - the index for which you want to get the name
407: Output Parameter:
408: . name - pointer to the string containing the name characters
410: Level: intermediate
412: .seealso: `PETSCVIEWEREXODUSII`, `PetscViewer`, `PetscViewerCreate()`, `PetscViewerDestroy()`, `PetscViewerExodusIIOpen()`, `PetscViewerSetType()`, `PetscViewerType`, `PetscViewerExodusIISetZonalVariableName()`
413: @*/
414: PetscErrorCode PetscViewerExodusIIGetZonalVariableName(PetscViewer viewer, PetscExodusIIInt idx, const char *name[])
415: {
416: PetscViewer_ExodusII *exo = (PetscViewer_ExodusII *)viewer->data;
417: PetscExodusIIInt exoid = -1;
418: char tmpName[MAX_NAME_LENGTH + 1];
420: PetscFunctionBegin;
421: PetscCheck(idx >= 0 && idx < exo->numZonalVariables, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Variable index out of range. Was PetscViewerExodusIISetZonalVariable called?");
422: if (!exo->zonalVariableNames[idx]) {
423: PetscCall(PetscViewerExodusIIGetId(viewer, &exoid));
424: PetscCallExternal(ex_get_variable_name, exoid, EX_ELEM_BLOCK, idx + 1, tmpName);
425: PetscCall(PetscStrallocpy(tmpName, (char **)&exo->zonalVariableNames[idx]));
426: }
427: *name = exo->zonalVariableNames[idx];
428: PetscFunctionReturn(PETSC_SUCCESS);
429: }
431: /*@
432: PetscViewerExodusIIGetNodalVariableName - Gets the name of a nodal variable.
434: Collective;
436: Input Parameters:
437: + viewer - a `PetscViewer` of type `PETSCVIEWEREXODUSII`
438: - idx - the index for which you want to save the name
440: Output Parameter:
441: . name - string array containing name characters
443: Level: intermediate
445: .seealso: `PETSCVIEWEREXODUSII`, `PetscViewer`, `PetscViewerCreate()`, `PetscViewerDestroy()`, `PetscViewerExodusIIOpen()`, `PetscViewerSetType()`, `PetscViewerType`, `PetscViewerExodusIISetNodalVariableName()`
446: @*/
447: PetscErrorCode PetscViewerExodusIIGetNodalVariableName(PetscViewer viewer, PetscExodusIIInt idx, const char *name[])
448: {
449: PetscViewer_ExodusII *exo = (PetscViewer_ExodusII *)viewer->data;
450: PetscExodusIIInt exoid = -1;
451: char tmpName[MAX_NAME_LENGTH + 1];
453: PetscFunctionBegin;
454: PetscCheck((idx >= 0) && (idx < exo->numNodalVariables), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Variable index out of range. Was PetscViewerExodusIISetNodalVariable called?");
455: if (!exo->nodalVariableNames[idx]) {
456: PetscCall(PetscViewerExodusIIGetId(viewer, &exoid));
457: PetscCallExternal(ex_get_variable_name, exoid, EX_NODAL, idx + 1, tmpName);
458: PetscCall(PetscStrallocpy(tmpName, (char **)&exo->nodalVariableNames[idx]));
459: }
460: *name = exo->nodalVariableNames[idx];
461: PetscFunctionReturn(PETSC_SUCCESS);
462: }
464: /*@C
465: PetscViewerExodusIISetZonalVariableNames - Sets the names of all nodal variables
467: Collective; No Fortran Support
469: Input Parameters:
470: + viewer - a `PetscViewer` of type `PETSCVIEWEREXODUSII`
471: - names - 2D array containing the array of string names to be set
473: Level: intermediate
475: Notes:
476: This function allows users to set multiple zonal variable names at a time.
478: .seealso: `PETSCVIEWEREXODUSII`, `PetscViewer`, `PetscViewerCreate()`, `PetscViewerDestroy()`, `PetscViewerExodusIIOpen()`, `PetscViewerSetType()`, `PetscViewerType`, `PetscViewerExodusIIGetZonalVariableNames()`
479: @*/
480: PetscErrorCode PetscViewerExodusIISetZonalVariableNames(PetscViewer viewer, const char *names[])
481: {
482: PetscExodusIIInt numNames;
483: PetscExodusIIInt exoid = -1;
484: PetscViewer_ExodusII *exo = (PetscViewer_ExodusII *)viewer->data;
486: PetscFunctionBegin;
487: PetscCall(PetscViewerExodusIIGetZonalVariable(viewer, &numNames));
488: PetscCheck(numNames >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Number of zonal variables not set. Was PetscViewerExodusIISetZonalVariable called?");
490: PetscCall(PetscViewerExodusIIGetId(viewer, &exoid));
491: for (int i = 0; i < numNames; i++) {
492: PetscCall(PetscStrallocpy(names[i], &exo->zonalVariableNames[i]));
493: PetscCallExternal(ex_put_variable_name, exoid, EX_ELEM_BLOCK, i + 1, names[i]);
494: }
495: PetscFunctionReturn(PETSC_SUCCESS);
496: }
498: /*@C
499: PetscViewerExodusIISetNodalVariableNames - Sets the names of all nodal variables.
501: Collective; No Fortran Support
503: Input Parameters:
504: + viewer - a `PetscViewer` of type `PETSCVIEWEREXODUSII`
505: - names - 2D array containing the array of string names to be set
507: Level: intermediate
509: Notes:
510: This function allows users to set multiple nodal variable names at a time.
512: .seealso: `PETSCVIEWEREXODUSII`, `PetscViewer`, `PetscViewerCreate()`, `PetscViewerDestroy()`, `PetscViewerExodusIIOpen()`, `PetscViewerSetType()`, `PetscViewerType`, `PetscViewerExodusIIGetNodalVariableNames()`
513: @*/
514: PetscErrorCode PetscViewerExodusIISetNodalVariableNames(PetscViewer viewer, const char *names[])
515: {
516: PetscExodusIIInt numNames;
517: PetscExodusIIInt exoid = -1;
518: PetscViewer_ExodusII *exo = (PetscViewer_ExodusII *)viewer->data;
520: PetscFunctionBegin;
521: PetscCall(PetscViewerExodusIIGetNodalVariable(viewer, &numNames));
522: PetscCheck(numNames >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Number of nodal variables not set. Was PetscViewerExodusIISetNodalVariable called?");
524: PetscCall(PetscViewerExodusIIGetId(viewer, &exoid));
525: for (int i = 0; i < numNames; i++) {
526: PetscCall(PetscStrallocpy(names[i], &exo->nodalVariableNames[i]));
527: PetscCallExternal(ex_put_variable_name, exoid, EX_NODAL, i + 1, names[i]);
528: }
529: PetscFunctionReturn(PETSC_SUCCESS);
530: }
532: /*@C
533: PetscViewerExodusIIGetZonalVariableNames - Gets the names of all zonal variables.
535: Collective; No Fortran Support
537: Input Parameters:
538: + viewer - a `PetscViewer` of type `PETSCVIEWEREXODUSII`
539: - numVars - the number of zonal variable names to retrieve
541: Output Parameters:
542: . varNames - pointer to a 2D array where the zonal variable names will be saved
544: Level: intermediate
546: Notes:
547: This function returns a borrowed pointer which should not be freed.
549: .seealso: `PETSCVIEWEREXODUSII`, `PetscViewer`, `PetscViewerCreate()`, `PetscViewerDestroy()`, `PetscViewerExodusIIOpen()`, `PetscViewerSetType()`, `PetscViewerType`, `PetscViewerExodusIISetZonalVariableNames()`
550: @*/
551: PetscErrorCode PetscViewerExodusIIGetZonalVariableNames(PetscViewer viewer, PetscExodusIIInt *numVars, char ***varNames)
552: {
553: PetscViewer_ExodusII *exo = (PetscViewer_ExodusII *)viewer->data;
554: PetscExodusIIInt idx;
555: char tmpName[MAX_NAME_LENGTH + 1];
556: PetscExodusIIInt exoid = -1;
558: PetscFunctionBegin;
559: PetscCall(PetscViewerExodusIIGetZonalVariable(viewer, numVars));
560: /*
561: Cache variable names if necessary
562: */
563: for (idx = 0; idx < *numVars; idx++) {
564: if (!exo->zonalVariableNames[idx]) {
565: PetscCall(PetscViewerExodusIIGetId(viewer, &exoid));
566: PetscCallExternal(ex_get_variable_name, exoid, EX_ELEM_BLOCK, idx + 1, tmpName);
567: PetscCall(PetscStrallocpy(tmpName, (char **)&exo->zonalVariableNames[idx]));
568: }
569: }
570: *varNames = (char **)exo->zonalVariableNames;
571: PetscFunctionReturn(PETSC_SUCCESS);
572: }
574: /*@C
575: PetscViewerExodusIIGetNodalVariableNames - Gets the names of all nodal variables.
577: Collective; No Fortran Support
579: Input Parameters:
580: + viewer - a `PetscViewer` of type `PETSCVIEWEREXODUSII`
581: - numVars - the number of nodal variable names to retrieve
583: Output Parameters:
584: . varNames - 2D array where the nodal variable names will be saved
586: Level: intermediate
588: Notes:
589: This function returns a borrowed pointer which should not be freed.
591: .seealso: `PETSCVIEWEREXODUSII`, `PetscViewer`, `PetscViewerCreate()`, `PetscViewerDestroy()`, `PetscViewerExodusIIOpen()`, `PetscViewerSetType()`, `PetscViewerType`, `PetscViewerExodusIISetNodalVariableNames()`
592: @*/
593: PetscErrorCode PetscViewerExodusIIGetNodalVariableNames(PetscViewer viewer, PetscExodusIIInt *numVars, char ***varNames)
594: {
595: PetscViewer_ExodusII *exo = (PetscViewer_ExodusII *)viewer->data;
596: PetscExodusIIInt idx;
597: char tmpName[MAX_NAME_LENGTH + 1];
598: PetscExodusIIInt exoid = -1;
600: PetscFunctionBegin;
601: PetscCall(PetscViewerExodusIIGetNodalVariable(viewer, numVars));
602: /*
603: Cache variable names if necessary
604: */
605: for (idx = 0; idx < *numVars; idx++) {
606: if (!exo->nodalVariableNames[idx]) {
607: PetscCall(PetscViewerExodusIIGetId(viewer, &exoid));
608: PetscCallExternal(ex_get_variable_name, exoid, EX_NODAL, idx + 1, tmpName);
609: PetscCall(PetscStrallocpy(tmpName, (char **)&exo->nodalVariableNames[idx]));
610: }
611: }
612: *varNames = (char **)exo->nodalVariableNames;
613: PetscFunctionReturn(PETSC_SUCCESS);
614: }
616: /*MC
617: PETSCVIEWEREXODUSII - A viewer that writes to an Exodus II file
619: Level: beginner
621: .seealso: `PetscViewerExodusIIOpen()`, `PetscViewerCreate()`, `PETSCVIEWERBINARY`, `PETSCVIEWERHDF5`, `DMView()`,
622: `PetscViewerFileSetName()`, `PetscViewerFileSetMode()`, `PetscViewerFormat`, `PetscViewerType`, `PetscViewerSetType()`
623: M*/
624: PETSC_EXTERN PetscErrorCode PetscViewerCreate_ExodusII(PetscViewer v)
625: {
626: PetscViewer_ExodusII *exo;
628: PetscFunctionBegin;
629: PetscCall(PetscNew(&exo));
631: v->data = (void *)exo;
632: v->ops->destroy = PetscViewerDestroy_ExodusII;
633: v->ops->setfromoptions = PetscViewerSetFromOptions_ExodusII;
634: v->ops->setup = PetscViewerSetUp_ExodusII;
635: v->ops->view = PetscViewerView_ExodusII;
636: v->ops->flush = PetscViewerFlush_ExodusII;
637: exo->btype = FILE_MODE_UNDEFINED;
638: exo->filename = 0;
639: exo->exoid = -1;
640: exo->numNodalVariables = -1;
641: exo->numZonalVariables = -1;
642: exo->nodalVariableNames = NULL;
643: exo->zonalVariableNames = NULL;
645: PetscCall(PetscObjectComposeFunction((PetscObject)v, "PetscViewerFileSetName_C", PetscViewerFileSetName_ExodusII));
646: PetscCall(PetscObjectComposeFunction((PetscObject)v, "PetscViewerFileGetName_C", PetscViewerFileGetName_ExodusII));
647: PetscCall(PetscObjectComposeFunction((PetscObject)v, "PetscViewerFileSetMode_C", PetscViewerFileSetMode_ExodusII));
648: PetscCall(PetscObjectComposeFunction((PetscObject)v, "PetscViewerFileGetMode_C", PetscViewerFileGetMode_ExodusII));
649: PetscCall(PetscObjectComposeFunction((PetscObject)v, "PetscViewerGetId_C", PetscViewerExodusIIGetId_ExodusII));
650: PetscCall(PetscObjectComposeFunction((PetscObject)v, "PetscViewerSetOrder_C", PetscViewerExodusIISetOrder_ExodusII));
651: PetscCall(PetscObjectComposeFunction((PetscObject)v, "PetscViewerGetOrder_C", PetscViewerExodusIIGetOrder_ExodusII));
652: PetscFunctionReturn(PETSC_SUCCESS);
653: }
655: /*@
656: PetscViewerExodusIIGetNodalVariableIndex - return the location of a nodal variable in an exodusII file given its name
658: Collective
660: Input Parameters:
661: + viewer - a `PetscViewer` of type `PETSCVIEWEREXODUSII`
662: - name - the name of the result
664: Output Parameter:
665: . varIndex - the location of the variable in the exodus file or -1 if the variable is not found
667: Level: beginner
669: Notes:
670: The exodus variable index is obtained by comparing the name argument to the
671: names of zonal variables declared in the exodus file. For instance if name is "V"
672: the location in the exodus file will be the first match of "V", "V_X", "V_XX", "V_1", or "V_11"
673: amongst all variables of type obj_type.
675: .seealso: `PetscViewerExodusIISetNodalVariable()`, `PetscViewerExodusIIGetNodalVariable()`, `PetscViewerExodusIISetNodalVariableName()`, `PetscViewerExodusIISetNodalVariableNames()`, `PetscViewerExodusIIGetNodalVariableName()`, `PetscViewerExodusIIGetNodalVariableNames()`
676: @*/
677: PetscErrorCode PetscViewerExodusIIGetNodalVariableIndex(PetscViewer viewer, const char name[], PetscExodusIIInt *varIndex)
678: {
679: PetscExodusIIInt num_vars = 0, i, j;
680: char ext_name[MAX_STR_LENGTH + 1];
681: char **var_names;
682: const int num_suffix = 5;
683: char *suffix[5];
684: PetscBool flg;
686: PetscFunctionBegin;
687: suffix[0] = (char *)"";
688: suffix[1] = (char *)"_X";
689: suffix[2] = (char *)"_XX";
690: suffix[3] = (char *)"_1";
691: suffix[4] = (char *)"_11";
692: *varIndex = -1;
694: PetscCall(PetscViewerExodusIIGetNodalVariableNames(viewer, &num_vars, &var_names));
695: for (i = 0; i < num_vars; ++i) {
696: for (j = 0; j < num_suffix; ++j) {
697: PetscCall(PetscStrncpy(ext_name, name, MAX_STR_LENGTH));
698: PetscCall(PetscStrlcat(ext_name, suffix[j], MAX_STR_LENGTH));
699: PetscCall(PetscStrcasecmp(ext_name, var_names[i], &flg));
700: if (flg) *varIndex = i;
701: }
702: if (flg) break;
703: }
704: PetscFunctionReturn(PETSC_SUCCESS);
705: }
707: /*@
708: PetscViewerExodusIIGetZonalVariableIndex - return the location of a zonal variable in an exodusII file given its name
710: Collective
712: Input Parameters:
713: + viewer - a `PetscViewer` of type `PETSCVIEWEREXODUSII`
714: - name - the name of the result
716: Output Parameter:
717: . varIndex - the location of the variable in the exodus file or -1 if the variable is not found
719: Level: beginner
721: Notes:
722: The exodus variable index is obtained by comparing the name argument to the
723: names of zonal variables declared in the exodus file. For instance if name is "V"
724: the location in the exodus file will be the first match of "V", "V_X", "V_XX", "V_1", or "V_11"
725: amongst all variables of type obj_type.
727: .seealso: `PetscViewerExodusIISetNodalVariable()`, `PetscViewerExodusIIGetNodalVariable()`, `PetscViewerExodusIISetNodalVariableName()`, `PetscViewerExodusIISetNodalVariableNames()`, `PetscViewerExodusIIGetNodalVariableName()`, `PetscViewerExodusIIGetNodalVariableNames()`
728: @*/
729: PetscErrorCode PetscViewerExodusIIGetZonalVariableIndex(PetscViewer viewer, const char name[], int *varIndex)
730: {
731: PetscExodusIIInt num_vars = 0, i, j;
732: char ext_name[MAX_STR_LENGTH + 1];
733: char **var_names;
734: const int num_suffix = 5;
735: char *suffix[5];
736: PetscBool flg;
738: PetscFunctionBegin;
739: suffix[0] = (char *)"";
740: suffix[1] = (char *)"_X";
741: suffix[2] = (char *)"_XX";
742: suffix[3] = (char *)"_1";
743: suffix[4] = (char *)"_11";
744: *varIndex = -1;
746: PetscCall(PetscViewerExodusIIGetZonalVariableNames(viewer, &num_vars, &var_names));
747: for (i = 0; i < num_vars; ++i) {
748: for (j = 0; j < num_suffix; ++j) {
749: PetscCall(PetscStrncpy(ext_name, name, MAX_STR_LENGTH));
750: PetscCall(PetscStrlcat(ext_name, suffix[j], MAX_STR_LENGTH));
751: PetscCall(PetscStrcasecmp(ext_name, var_names[i], &flg));
752: if (flg) *varIndex = i;
753: }
754: if (flg) break;
755: }
756: PetscFunctionReturn(PETSC_SUCCESS);
757: }
759: PetscErrorCode DMView_PlexExodusII(DM dm, PetscViewer viewer)
760: {
761: enum ElemType {
762: SEGMENT,
763: TRI,
764: QUAD,
765: TET,
766: HEX
767: };
768: MPI_Comm comm;
769: PetscInt degree; /* the order of the mesh */
770: /* Connectivity Variables */
771: PetscInt cellsNotInConnectivity;
772: /* Cell Sets */
773: DMLabel csLabel;
774: IS csIS;
775: const PetscInt *csIdx;
776: PetscInt num_cs, cs;
777: enum ElemType *type;
778: PetscBool hasLabel;
779: /* Coordinate Variables */
780: DM cdm;
781: PetscSection coordSection;
782: Vec coord;
783: PetscInt **nodes;
784: PetscInt depth, d, dim, skipCells = 0;
785: PetscInt pStart, pEnd, p, cStart, cEnd, numCells, vStart, vEnd, numVertices, eStart, eEnd, numEdges, fStart, fEnd, numFaces, numNodes;
786: PetscInt num_vs, num_fs;
787: PetscMPIInt rank, size;
788: const char *dmName;
789: PetscInt nodesLineP1[4] = {2, 0, 0, 0};
790: PetscInt nodesLineP2[4] = {2, 0, 0, 1};
791: PetscInt nodesTriP1[4] = {3, 0, 0, 0};
792: PetscInt nodesTriP2[4] = {3, 3, 0, 0};
793: PetscInt nodesQuadP1[4] = {4, 0, 0, 0};
794: PetscInt nodesQuadP2[4] = {4, 4, 0, 1};
795: PetscInt nodesTetP1[4] = {4, 0, 0, 0};
796: PetscInt nodesTetP2[4] = {4, 6, 0, 0};
797: PetscInt nodesHexP1[4] = {8, 0, 0, 0};
798: PetscInt nodesHexP2[4] = {8, 12, 6, 1};
799: PetscExodusIIInt CPU_word_size, IO_word_size, EXO_mode;
800: PetscExodusIIFloat EXO_version;
802: PetscViewer_ExodusII *exo = (PetscViewer_ExodusII *)viewer->data;
804: PetscFunctionBegin;
805: PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
806: PetscCallMPI(MPI_Comm_rank(comm, &rank));
807: PetscCallMPI(MPI_Comm_size(comm, &size));
809: /*
810: Creating coordSection is a collective operation so we do it somewhat out of sequence
811: */
812: PetscCall(PetscSectionCreate(comm, &coordSection));
813: PetscCall(DMGetCoordinatesLocalSetUp(dm));
814: /*
815: Check that all points are on rank 0 since we don't know how to save distributed DM in exodus format
816: */
817: PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
818: PetscCall(DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd));
819: PetscCall(DMPlexGetDepthStratum(dm, 1, &eStart, &eEnd));
820: PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd));
821: numCells = cEnd - cStart;
822: numEdges = eEnd - eStart;
823: numVertices = vEnd - vStart;
824: PetscCheck(!(rank && (numCells || numEdges || numVertices)), PETSC_COMM_SELF, PETSC_ERR_SUP, "Writing distributed DM in exodusII format not supported");
825: if (rank == 0) {
826: switch (exo->btype) {
827: case FILE_MODE_READ:
828: case FILE_MODE_APPEND:
829: case FILE_MODE_UPDATE:
830: case FILE_MODE_APPEND_UPDATE:
831: /* exodusII does not allow writing geometry to an existing file */
832: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "cannot add geometry to existing file %s", exo->filename);
833: case FILE_MODE_WRITE:
834: /* Create an empty file if one already exists*/
835: EXO_mode = EX_CLOBBER;
836: #if defined(PETSC_USE_64BIT_INDICES)
837: EXO_mode += EX_ALL_INT64_API;
838: #endif
839: CPU_word_size = sizeof(PetscReal);
840: IO_word_size = sizeof(PetscReal);
841: exo->exoid = ex_create(exo->filename, EXO_mode, &CPU_word_size, &IO_word_size);
842: PetscCheck(exo->exoid >= 0, PETSC_COMM_SELF, PETSC_ERR_LIB, "ex_create failed for %s", exo->filename);
844: break;
845: default:
846: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ORDER, "Must call PetscViewerFileSetMode() before PetscViewerFileSetName()");
847: }
849: /* --- Get DM info --- */
850: PetscCall(PetscObjectGetName((PetscObject)dm, &dmName));
851: PetscCall(DMPlexGetDepth(dm, &depth));
852: PetscCall(DMGetDimension(dm, &dim));
853: PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
854: if (depth == 3) {
855: numFaces = fEnd - fStart;
856: } else {
857: numFaces = 0;
858: }
859: PetscCall(DMGetLabelSize(dm, "Cell Sets", &num_cs));
860: PetscCall(DMGetLabelSize(dm, "Vertex Sets", &num_vs));
861: PetscCall(DMGetLabelSize(dm, "Face Sets", &num_fs));
862: PetscCall(DMGetCoordinatesLocal(dm, &coord));
863: PetscCall(DMGetCoordinateDM(dm, &cdm));
864: if (num_cs > 0) {
865: PetscCall(DMGetLabel(dm, "Cell Sets", &csLabel));
866: PetscCall(DMLabelGetValueIS(csLabel, &csIS));
867: PetscCall(ISGetIndices(csIS, &csIdx));
868: }
869: PetscCall(PetscMalloc1(num_cs, &nodes));
870: /* Set element type for each block and compute total number of nodes */
871: PetscCall(PetscMalloc1(num_cs, &type));
872: numNodes = numVertices;
874: PetscCall(PetscViewerExodusIIGetOrder(viewer, °ree));
875: if (degree == 2) numNodes += numEdges;
876: cellsNotInConnectivity = numCells;
877: for (cs = 0; cs < num_cs; ++cs) {
878: IS stratumIS;
879: const PetscInt *cells;
880: PetscScalar *xyz = NULL;
881: PetscInt csSize, closureSize;
883: PetscCall(DMLabelGetStratumIS(csLabel, csIdx[cs], &stratumIS));
884: PetscCall(ISGetIndices(stratumIS, &cells));
885: PetscCall(ISGetSize(stratumIS, &csSize));
886: PetscCall(DMPlexVecGetClosure(cdm, NULL, coord, cells[0], &closureSize, &xyz));
887: switch (dim) {
888: case 1:
889: if (closureSize == 2 * dim) {
890: type[cs] = SEGMENT;
891: } 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);
892: break;
893: case 2:
894: if (closureSize == 3 * dim) {
895: type[cs] = TRI;
896: } else if (closureSize == 4 * dim) {
897: type[cs] = QUAD;
898: } 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);
899: break;
900: case 3:
901: if (closureSize == 4 * dim) {
902: type[cs] = TET;
903: } else if (closureSize == 8 * dim) {
904: type[cs] = HEX;
905: } 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);
906: break;
907: default:
908: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Dimension %" PetscInt_FMT " not handled by ExodusII viewer", dim);
909: }
910: if ((degree == 2) && (type[cs] == SEGMENT)) numNodes += csSize;
911: if ((degree == 2) && (type[cs] == QUAD)) numNodes += csSize;
912: if ((degree == 2) && (type[cs] == HEX)) {
913: numNodes += csSize;
914: numNodes += numFaces;
915: }
916: PetscCall(DMPlexVecRestoreClosure(cdm, NULL, coord, cells[0], &closureSize, &xyz));
917: /* Set nodes and Element type */
918: if (type[cs] == SEGMENT) {
919: if (degree == 1) nodes[cs] = nodesLineP1;
920: else if (degree == 2) nodes[cs] = nodesLineP2;
921: } else if (type[cs] == TRI) {
922: if (degree == 1) nodes[cs] = nodesTriP1;
923: else if (degree == 2) nodes[cs] = nodesTriP2;
924: } else if (type[cs] == QUAD) {
925: if (degree == 1) nodes[cs] = nodesQuadP1;
926: else if (degree == 2) nodes[cs] = nodesQuadP2;
927: } else if (type[cs] == TET) {
928: if (degree == 1) nodes[cs] = nodesTetP1;
929: else if (degree == 2) nodes[cs] = nodesTetP2;
930: } else if (type[cs] == HEX) {
931: if (degree == 1) nodes[cs] = nodesHexP1;
932: else if (degree == 2) nodes[cs] = nodesHexP2;
933: }
934: /* Compute the number of cells not in the connectivity table */
935: cellsNotInConnectivity -= nodes[cs][3] * csSize;
937: PetscCall(ISRestoreIndices(stratumIS, &cells));
938: PetscCall(ISDestroy(&stratumIS));
939: }
940: if (num_cs) PetscCallExternal(ex_put_init, exo->exoid, dmName, dim, numNodes, numCells, num_cs, num_vs, num_fs);
941: /* --- Connectivity --- */
942: for (cs = 0; cs < num_cs; ++cs) {
943: IS stratumIS;
944: const PetscInt *cells;
945: PetscInt *connect, off = 0;
946: PetscInt edgesInClosure = 0, facesInClosure = 0, verticesInClosure = 0;
947: PetscInt csSize, c, connectSize, closureSize;
948: char *elem_type = NULL;
949: char elem_type_bar2[] = "BAR2", elem_type_bar3[] = "BAR3";
950: char elem_type_tri3[] = "TRI3", elem_type_quad4[] = "QUAD4";
951: char elem_type_tri6[] = "TRI6", elem_type_quad9[] = "QUAD9";
952: char elem_type_tet4[] = "TET4", elem_type_hex8[] = "HEX8";
953: char elem_type_tet10[] = "TET10", elem_type_hex27[] = "HEX27";
955: PetscCall(DMLabelGetStratumIS(csLabel, csIdx[cs], &stratumIS));
956: PetscCall(ISGetIndices(stratumIS, &cells));
957: PetscCall(ISGetSize(stratumIS, &csSize));
958: /* Set Element type */
959: if (type[cs] == SEGMENT) {
960: if (degree == 1) elem_type = elem_type_bar2;
961: else if (degree == 2) elem_type = elem_type_bar3;
962: } else if (type[cs] == TRI) {
963: if (degree == 1) elem_type = elem_type_tri3;
964: else if (degree == 2) elem_type = elem_type_tri6;
965: } else if (type[cs] == QUAD) {
966: if (degree == 1) elem_type = elem_type_quad4;
967: else if (degree == 2) elem_type = elem_type_quad9;
968: } else if (type[cs] == TET) {
969: if (degree == 1) elem_type = elem_type_tet4;
970: else if (degree == 2) elem_type = elem_type_tet10;
971: } else if (type[cs] == HEX) {
972: if (degree == 1) elem_type = elem_type_hex8;
973: else if (degree == 2) elem_type = elem_type_hex27;
974: }
975: connectSize = nodes[cs][0] + nodes[cs][1] + nodes[cs][2] + nodes[cs][3];
976: PetscCall(PetscMalloc1(PetscMax(27, connectSize) * csSize, &connect));
977: PetscCallExternal(ex_put_block, exo->exoid, EX_ELEM_BLOCK, csIdx[cs], elem_type, csSize, connectSize, 0, 0, 1);
978: /* Find number of vertices, edges, and faces in the closure */
979: verticesInClosure = nodes[cs][0];
980: if (depth > 1) {
981: if (dim == 2) {
982: PetscCall(DMPlexGetConeSize(dm, cells[0], &edgesInClosure));
983: } else if (dim == 3) {
984: PetscInt *closure = NULL;
986: PetscCall(DMPlexGetConeSize(dm, cells[0], &facesInClosure));
987: PetscCall(DMPlexGetTransitiveClosure(dm, cells[0], PETSC_TRUE, &closureSize, &closure));
988: edgesInClosure = closureSize - facesInClosure - 1 - verticesInClosure;
989: PetscCall(DMPlexRestoreTransitiveClosure(dm, cells[0], PETSC_TRUE, &closureSize, &closure));
990: }
991: }
992: /* Get connectivity for each cell */
993: for (c = 0; c < csSize; ++c) {
994: PetscInt *closure = NULL;
995: PetscInt temp, i;
997: PetscCall(DMPlexGetTransitiveClosure(dm, cells[c], PETSC_TRUE, &closureSize, &closure));
998: for (i = 0; i < connectSize; ++i) {
999: if (i < nodes[cs][0]) { /* Vertices */
1000: connect[i + off] = closure[(i + edgesInClosure + facesInClosure + 1) * 2] + 1;
1001: connect[i + off] -= cellsNotInConnectivity;
1002: } else if (i < nodes[cs][0] + nodes[cs][1]) { /* Edges */
1003: connect[i + off] = closure[(i - verticesInClosure + facesInClosure + 1) * 2] + 1;
1004: if (nodes[cs][2] == 0) connect[i + off] -= numFaces;
1005: connect[i + off] -= cellsNotInConnectivity;
1006: } else if (i < nodes[cs][0] + nodes[cs][1] + nodes[cs][3]) { /* Cells */
1007: connect[i + off] = closure[0] + 1;
1008: connect[i + off] -= skipCells;
1009: } else if (i < nodes[cs][0] + nodes[cs][1] + nodes[cs][3] + nodes[cs][2]) { /* Faces */
1010: connect[i + off] = closure[(i - edgesInClosure - verticesInClosure) * 2] + 1;
1011: connect[i + off] -= cellsNotInConnectivity;
1012: } else {
1013: connect[i + off] = -1;
1014: }
1015: }
1016: /* Tetrahedra are inverted */
1017: if (type[cs] == TET) {
1018: temp = connect[0 + off];
1019: connect[0 + off] = connect[1 + off];
1020: connect[1 + off] = temp;
1021: if (degree == 2) {
1022: temp = connect[5 + off];
1023: connect[5 + off] = connect[6 + off];
1024: connect[6 + off] = temp;
1025: temp = connect[7 + off];
1026: connect[7 + off] = connect[8 + off];
1027: connect[8 + off] = temp;
1028: }
1029: }
1030: /* Hexahedra are inverted */
1031: if (type[cs] == HEX) {
1032: temp = connect[1 + off];
1033: connect[1 + off] = connect[3 + off];
1034: connect[3 + off] = temp;
1035: if (degree == 2) {
1036: temp = connect[8 + off];
1037: connect[8 + off] = connect[11 + off];
1038: connect[11 + off] = temp;
1039: temp = connect[9 + off];
1040: connect[9 + off] = connect[10 + off];
1041: connect[10 + off] = temp;
1042: temp = connect[16 + off];
1043: connect[16 + off] = connect[17 + off];
1044: connect[17 + off] = temp;
1045: temp = connect[18 + off];
1046: connect[18 + off] = connect[19 + off];
1047: connect[19 + off] = temp;
1049: temp = connect[12 + off];
1050: connect[12 + off] = connect[16 + off];
1051: connect[16 + off] = temp;
1052: temp = connect[13 + off];
1053: connect[13 + off] = connect[17 + off];
1054: connect[17 + off] = temp;
1055: temp = connect[14 + off];
1056: connect[14 + off] = connect[18 + off];
1057: connect[18 + off] = temp;
1058: temp = connect[15 + off];
1059: connect[15 + off] = connect[19 + off];
1060: connect[19 + off] = temp;
1062: temp = connect[23 + off];
1063: connect[23 + off] = connect[26 + off];
1064: connect[26 + off] = temp;
1065: temp = connect[24 + off];
1066: connect[24 + off] = connect[25 + off];
1067: connect[25 + off] = temp;
1068: temp = connect[25 + off];
1069: connect[25 + off] = connect[26 + off];
1070: connect[26 + off] = temp;
1071: }
1072: }
1073: off += connectSize;
1074: PetscCall(DMPlexRestoreTransitiveClosure(dm, cells[c], PETSC_TRUE, &closureSize, &closure));
1075: }
1076: PetscCallExternal(ex_put_conn, exo->exoid, EX_ELEM_BLOCK, csIdx[cs], connect, 0, 0);
1077: skipCells += (nodes[cs][3] == 0) * csSize;
1078: PetscCall(PetscFree(connect));
1079: PetscCall(ISRestoreIndices(stratumIS, &cells));
1080: PetscCall(ISDestroy(&stratumIS));
1081: }
1082: PetscCall(PetscFree(type));
1083: /* --- Coordinates --- */
1084: PetscCall(PetscSectionSetChart(coordSection, pStart, pEnd));
1085: if (num_cs) {
1086: for (d = 0; d < depth; ++d) {
1087: PetscCall(DMPlexGetDepthStratum(dm, d, &pStart, &pEnd));
1088: for (p = pStart; p < pEnd; ++p) PetscCall(PetscSectionSetDof(coordSection, p, nodes[0][d] > 0));
1089: }
1090: }
1091: for (cs = 0; cs < num_cs; ++cs) {
1092: IS stratumIS;
1093: const PetscInt *cells;
1094: PetscInt csSize, c;
1096: PetscCall(DMLabelGetStratumIS(csLabel, csIdx[cs], &stratumIS));
1097: PetscCall(ISGetIndices(stratumIS, &cells));
1098: PetscCall(ISGetSize(stratumIS, &csSize));
1099: for (c = 0; c < csSize; ++c) PetscCall(PetscSectionSetDof(coordSection, cells[c], nodes[cs][3] > 0));
1100: PetscCall(ISRestoreIndices(stratumIS, &cells));
1101: PetscCall(ISDestroy(&stratumIS));
1102: }
1103: if (num_cs) {
1104: PetscCall(ISRestoreIndices(csIS, &csIdx));
1105: PetscCall(ISDestroy(&csIS));
1106: }
1107: PetscCall(PetscFree(nodes));
1108: PetscCall(PetscSectionSetUp(coordSection));
1109: if (numNodes) {
1110: const char *coordNames[3] = {"x", "y", "z"};
1111: PetscScalar *closure, *cval;
1112: PetscReal *coords;
1113: PetscInt hasDof, n = 0;
1115: /* There can't be more than 24 values in the closure of a point for the coord coordSection */
1116: PetscCall(PetscCalloc3(numNodes * 3, &coords, dim, &cval, 24, &closure));
1117: PetscCall(DMGetCoordinatesLocalNoncollective(dm, &coord));
1118: PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
1119: for (p = pStart; p < pEnd; ++p) {
1120: PetscCall(PetscSectionGetDof(coordSection, p, &hasDof));
1121: if (hasDof) {
1122: PetscInt closureSize = 24, j;
1124: PetscCall(DMPlexVecGetClosure(cdm, NULL, coord, p, &closureSize, &closure));
1125: for (d = 0; d < dim; ++d) {
1126: cval[d] = 0.0;
1127: for (j = 0; j < closureSize / dim; j++) cval[d] += closure[j * dim + d];
1128: coords[d * numNodes + n] = PetscRealPart(cval[d]) * dim / closureSize;
1129: }
1130: ++n;
1131: }
1132: }
1133: PetscCallExternal(ex_put_coord, exo->exoid, &coords[0 * numNodes], &coords[1 * numNodes], &coords[2 * numNodes]);
1134: PetscCall(PetscFree3(coords, cval, closure));
1135: PetscCallExternal(ex_put_coord_names, exo->exoid, (char **)coordNames);
1136: }
1138: /* --- Node Sets/Vertex Sets --- */
1139: PetscCall(DMHasLabel(dm, "Vertex Sets", &hasLabel));
1140: if (hasLabel) {
1141: PetscInt i, vs, vsSize;
1142: const PetscInt *vsIdx, *vertices;
1143: PetscInt *nodeList;
1144: IS vsIS, stratumIS;
1145: DMLabel vsLabel;
1146: PetscCall(DMGetLabel(dm, "Vertex Sets", &vsLabel));
1147: PetscCall(DMLabelGetValueIS(vsLabel, &vsIS));
1148: PetscCall(ISGetIndices(vsIS, &vsIdx));
1149: for (vs = 0; vs < num_vs; ++vs) {
1150: PetscCall(DMLabelGetStratumIS(vsLabel, vsIdx[vs], &stratumIS));
1151: PetscCall(ISGetIndices(stratumIS, &vertices));
1152: PetscCall(ISGetSize(stratumIS, &vsSize));
1153: PetscCall(PetscMalloc1(vsSize, &nodeList));
1154: for (i = 0; i < vsSize; ++i) nodeList[i] = vertices[i] - skipCells + 1;
1155: PetscCallExternal(ex_put_set_param, exo->exoid, EX_NODE_SET, vsIdx[vs], vsSize, 0);
1156: PetscCallExternal(ex_put_set, exo->exoid, EX_NODE_SET, vsIdx[vs], nodeList, NULL);
1157: PetscCall(ISRestoreIndices(stratumIS, &vertices));
1158: PetscCall(ISDestroy(&stratumIS));
1159: PetscCall(PetscFree(nodeList));
1160: }
1161: PetscCall(ISRestoreIndices(vsIS, &vsIdx));
1162: PetscCall(ISDestroy(&vsIS));
1163: }
1164: /* --- Side Sets/Face Sets --- */
1165: PetscCall(DMHasLabel(dm, "Face Sets", &hasLabel));
1166: if (hasLabel) {
1167: PetscInt i, j, fs, fsSize;
1168: const PetscInt *fsIdx, *faces;
1169: IS fsIS, stratumIS;
1170: DMLabel fsLabel;
1171: PetscInt numPoints, *points;
1172: PetscInt elem_list_size = 0;
1173: PetscInt *elem_list, *elem_ind, *side_list;
1175: PetscCall(DMGetLabel(dm, "Face Sets", &fsLabel));
1176: /* Compute size of Node List and Element List */
1177: PetscCall(DMLabelGetValueIS(fsLabel, &fsIS));
1178: PetscCall(ISGetIndices(fsIS, &fsIdx));
1179: for (fs = 0; fs < num_fs; ++fs) {
1180: PetscCall(DMLabelGetStratumIS(fsLabel, fsIdx[fs], &stratumIS));
1181: PetscCall(ISGetSize(stratumIS, &fsSize));
1182: elem_list_size += fsSize;
1183: PetscCall(ISDestroy(&stratumIS));
1184: }
1185: if (num_fs) {
1186: PetscCall(PetscMalloc3(num_fs, &elem_ind, elem_list_size, &elem_list, elem_list_size, &side_list));
1187: elem_ind[0] = 0;
1188: for (fs = 0; fs < num_fs; ++fs) {
1189: PetscCall(DMLabelGetStratumIS(fsLabel, fsIdx[fs], &stratumIS));
1190: PetscCall(ISGetIndices(stratumIS, &faces));
1191: PetscCall(ISGetSize(stratumIS, &fsSize));
1192: /* Set Parameters */
1193: PetscCallExternal(ex_put_set_param, exo->exoid, EX_SIDE_SET, fsIdx[fs], fsSize, 0);
1194: /* Indices */
1195: if (fs < num_fs - 1) elem_ind[fs + 1] = elem_ind[fs] + fsSize;
1197: for (i = 0; i < fsSize; ++i) {
1198: /* Element List */
1199: points = NULL;
1200: PetscCall(DMPlexGetTransitiveClosure(dm, faces[i], PETSC_FALSE, &numPoints, &points));
1201: elem_list[elem_ind[fs] + i] = points[2] + 1;
1202: PetscCall(DMPlexRestoreTransitiveClosure(dm, faces[i], PETSC_FALSE, &numPoints, &points));
1204: /* Side List */
1205: points = NULL;
1206: PetscCall(DMPlexGetTransitiveClosure(dm, elem_list[elem_ind[fs] + i] - 1, PETSC_TRUE, &numPoints, &points));
1207: for (j = 1; j < numPoints; ++j) {
1208: if (points[j * 2] == faces[i]) break;
1209: }
1210: /* Convert HEX sides */
1211: if (numPoints == 27) {
1212: if (j == 1) {
1213: j = 5;
1214: } else if (j == 2) {
1215: j = 6;
1216: } else if (j == 3) {
1217: j = 1;
1218: } else if (j == 4) {
1219: j = 3;
1220: } else if (j == 5) {
1221: j = 2;
1222: } else if (j == 6) {
1223: j = 4;
1224: }
1225: }
1226: /* Convert TET sides */
1227: if (numPoints == 15) {
1228: --j;
1229: if (j == 0) j = 4;
1230: }
1231: side_list[elem_ind[fs] + i] = j;
1232: PetscCall(DMPlexRestoreTransitiveClosure(dm, elem_list[elem_ind[fs] + i] - 1, PETSC_TRUE, &numPoints, &points));
1233: }
1234: PetscCall(ISRestoreIndices(stratumIS, &faces));
1235: PetscCall(ISDestroy(&stratumIS));
1236: }
1237: PetscCall(ISRestoreIndices(fsIS, &fsIdx));
1238: PetscCall(ISDestroy(&fsIS));
1240: /* Put side sets */
1241: 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]]);
1242: PetscCall(PetscFree3(elem_ind, elem_list, side_list));
1243: }
1244: }
1245: /*
1246: close the exodus file
1247: */
1248: ex_close(exo->exoid);
1249: exo->exoid = -1;
1250: }
1251: PetscCall(PetscSectionDestroy(&coordSection));
1253: /*
1254: reopen the file in parallel
1255: */
1256: EXO_mode = EX_WRITE;
1257: #if defined(PETSC_USE_64BIT_INDICES)
1258: EXO_mode += EX_ALL_INT64_API;
1259: #endif
1260: CPU_word_size = sizeof(PetscReal);
1261: IO_word_size = sizeof(PetscReal);
1262: exo->exoid = ex_open_par(exo->filename, EXO_mode, &CPU_word_size, &IO_word_size, &EXO_version, comm, MPI_INFO_NULL);
1263: PetscCheck(exo->exoid >= 0, PETSC_COMM_SELF, PETSC_ERR_LIB, "ex_open_par failed for %s", exo->filename);
1264: PetscFunctionReturn(PETSC_SUCCESS);
1265: }
1267: static PetscErrorCode VecViewPlex_ExodusII_Nodal_Internal(Vec v, PetscExodusIIInt exoid, PetscExodusIIInt step, PetscExodusIIInt offset);
1268: static PetscErrorCode VecViewPlex_ExodusII_Zonal_Internal(Vec v, PetscExodusIIInt exoid, PetscExodusIIInt step, PetscExodusIIInt offset);
1269: static PetscErrorCode VecLoadPlex_ExodusII_Nodal_Internal(Vec v, PetscExodusIIInt exoid, PetscExodusIIInt step, PetscExodusIIInt offset);
1270: static PetscErrorCode VecLoadPlex_ExodusII_Zonal_Internal(Vec v, PetscExodusIIInt exoid, PetscExodusIIInt step, PetscExodusIIInt offset);
1272: PetscErrorCode VecView_PlexExodusII_Internal(Vec v, PetscViewer viewer)
1273: {
1274: DM dm;
1275: MPI_Comm comm;
1276: PetscMPIInt rank;
1278: PetscExodusIIInt exoid, offsetN = -1, offsetZ = -1;
1279: const char *vecname;
1280: PetscInt step;
1282: PetscFunctionBegin;
1283: PetscCall(PetscObjectGetComm((PetscObject)v, &comm));
1284: PetscCallMPI(MPI_Comm_rank(comm, &rank));
1285: PetscCall(PetscViewerExodusIIGetId(viewer, &exoid));
1286: PetscCall(VecGetDM(v, &dm));
1287: PetscCall(PetscObjectGetName((PetscObject)v, &vecname));
1289: PetscCall(DMGetOutputSequenceNumber(dm, &step, NULL));
1290: PetscCall(PetscViewerExodusIIGetNodalVariableIndex(viewer, vecname, &offsetN));
1291: PetscCall(PetscViewerExodusIIGetZonalVariableIndex(viewer, vecname, &offsetZ));
1292: PetscCheck(!(offsetN >= 0 && offsetZ >= 0), comm, PETSC_ERR_FILE_UNEXPECTED, "Found both nodal and zonal variable %s in exodus file. ", vecname);
1293: if (offsetN >= 0) {
1294: PetscCall(VecViewPlex_ExodusII_Nodal_Internal(v, exoid, (PetscExodusIIInt)step + 1, offsetN + 1));
1295: } else if (offsetZ >= 0) {
1296: PetscCall(VecViewPlex_ExodusII_Zonal_Internal(v, exoid, (PetscExodusIIInt)step + 1, offsetZ + 1));
1297: } else SETERRQ(comm, PETSC_ERR_FILE_UNEXPECTED, "Could not find nodal or zonal variable %s in exodus file. ", vecname);
1298: PetscFunctionReturn(PETSC_SUCCESS);
1299: }
1301: PetscErrorCode VecLoad_PlexExodusII_Internal(Vec v, PetscViewer viewer)
1302: {
1303: DM dm;
1304: MPI_Comm comm;
1305: PetscMPIInt rank;
1307: PetscExodusIIInt exoid, offsetN = 0, offsetZ = 0;
1308: const char *vecname;
1309: PetscInt step;
1311: PetscFunctionBegin;
1312: PetscCall(PetscObjectGetComm((PetscObject)v, &comm));
1313: PetscCallMPI(MPI_Comm_rank(comm, &rank));
1314: PetscCall(PetscViewerExodusIIGetId(viewer, &exoid));
1315: PetscCall(VecGetDM(v, &dm));
1316: PetscCall(PetscObjectGetName((PetscObject)v, &vecname));
1318: PetscCall(DMGetOutputSequenceNumber(dm, &step, NULL));
1319: PetscCall(PetscViewerExodusIIGetNodalVariableIndex(viewer, vecname, &offsetN));
1320: PetscCall(PetscViewerExodusIIGetZonalVariableIndex(viewer, vecname, &offsetZ));
1321: PetscCheck(!(offsetN >= 0 && offsetZ >= 0), comm, PETSC_ERR_FILE_UNEXPECTED, "Found both nodal and zonal variable %s in exodus file. ", vecname);
1322: if (offsetN >= 0) PetscCall(VecLoadPlex_ExodusII_Nodal_Internal(v, exoid, (PetscExodusIIInt)step + 1, offsetN + 1));
1323: else if (offsetZ >= 0) PetscCall(VecLoadPlex_ExodusII_Zonal_Internal(v, exoid, (PetscExodusIIInt)step + 1, offsetZ + 1));
1324: else SETERRQ(comm, PETSC_ERR_FILE_UNEXPECTED, "Could not find nodal or zonal variable %s in exodus file. ", vecname);
1325: PetscFunctionReturn(PETSC_SUCCESS);
1326: }
1328: static PetscErrorCode VecViewPlex_ExodusII_Nodal_Internal(Vec v, PetscExodusIIInt exoid, PetscExodusIIInt step, PetscExodusIIInt offset)
1329: {
1330: MPI_Comm comm;
1331: PetscMPIInt size;
1332: DM dm;
1333: Vec vNatural, vComp;
1334: const PetscScalar *varray;
1335: PetscInt xs, xe, bs;
1336: PetscBool useNatural;
1338: PetscFunctionBegin;
1339: PetscCall(PetscObjectGetComm((PetscObject)v, &comm));
1340: PetscCallMPI(MPI_Comm_size(comm, &size));
1341: PetscCall(VecGetDM(v, &dm));
1342: PetscCall(DMGetUseNatural(dm, &useNatural));
1343: useNatural = useNatural && size > 1 ? PETSC_TRUE : PETSC_FALSE;
1344: if (useNatural) {
1345: PetscCall(DMPlexCreateNaturalVector(dm, &vNatural));
1346: PetscCall(DMPlexGlobalToNaturalBegin(dm, v, vNatural));
1347: PetscCall(DMPlexGlobalToNaturalEnd(dm, v, vNatural));
1348: } else {
1349: vNatural = v;
1350: }
1352: /* Write local chunk of the result in the exodus file
1353: exodus stores each component of a vector-valued field as a separate variable.
1354: We assume that they are stored sequentially */
1355: PetscCall(VecGetOwnershipRange(vNatural, &xs, &xe));
1356: PetscCall(VecGetBlockSize(vNatural, &bs));
1357: if (bs == 1) {
1358: PetscCall(VecGetArrayRead(vNatural, &varray));
1359: PetscCallExternal(ex_put_partial_var, exoid, step, EX_NODAL, offset, 1, xs + 1, xe - xs, varray);
1360: PetscCall(VecRestoreArrayRead(vNatural, &varray));
1361: } else {
1362: IS compIS;
1363: PetscInt c;
1365: PetscCall(ISCreateStride(comm, (xe - xs) / bs, xs, bs, &compIS));
1366: for (c = 0; c < bs; ++c) {
1367: PetscCall(ISStrideSetStride(compIS, (xe - xs) / bs, xs + c, bs));
1368: PetscCall(VecGetSubVector(vNatural, compIS, &vComp));
1369: PetscCall(VecGetArrayRead(vComp, &varray));
1370: PetscCallExternal(ex_put_partial_var, exoid, step, EX_NODAL, offset + c, 1, xs / bs + 1, (xe - xs) / bs, varray);
1371: PetscCall(VecRestoreArrayRead(vComp, &varray));
1372: PetscCall(VecRestoreSubVector(vNatural, compIS, &vComp));
1373: }
1374: PetscCall(ISDestroy(&compIS));
1375: }
1376: if (useNatural) PetscCall(VecDestroy(&vNatural));
1377: PetscFunctionReturn(PETSC_SUCCESS);
1378: }
1380: static PetscErrorCode VecLoadPlex_ExodusII_Nodal_Internal(Vec v, PetscExodusIIInt exoid, PetscExodusIIInt step, PetscExodusIIInt offset)
1381: {
1382: MPI_Comm comm;
1383: PetscMPIInt size;
1384: DM dm;
1385: Vec vNatural, vComp;
1386: PetscScalar *varray;
1387: PetscInt xs, xe, bs;
1388: PetscBool useNatural;
1390: PetscFunctionBegin;
1391: PetscCall(PetscObjectGetComm((PetscObject)v, &comm));
1392: PetscCallMPI(MPI_Comm_size(comm, &size));
1393: PetscCall(VecGetDM(v, &dm));
1394: PetscCall(DMGetUseNatural(dm, &useNatural));
1395: useNatural = useNatural && size > 1 ? PETSC_TRUE : PETSC_FALSE;
1396: if (useNatural) PetscCall(DMPlexCreateNaturalVector(dm, &vNatural));
1397: else vNatural = v;
1399: /* Read local chunk from the file */
1400: PetscCall(VecGetOwnershipRange(vNatural, &xs, &xe));
1401: PetscCall(VecGetBlockSize(vNatural, &bs));
1402: if (bs == 1) {
1403: PetscCall(VecGetArray(vNatural, &varray));
1404: PetscCallExternal(ex_get_partial_var, exoid, step, EX_NODAL, offset, 1, xs + 1, xe - xs, varray);
1405: PetscCall(VecRestoreArray(vNatural, &varray));
1406: } else {
1407: IS compIS;
1408: PetscInt c;
1410: PetscCall(ISCreateStride(comm, (xe - xs) / bs, xs, bs, &compIS));
1411: for (c = 0; c < bs; ++c) {
1412: PetscCall(ISStrideSetStride(compIS, (xe - xs) / bs, xs + c, bs));
1413: PetscCall(VecGetSubVector(vNatural, compIS, &vComp));
1414: PetscCall(VecGetArray(vComp, &varray));
1415: PetscCallExternal(ex_get_partial_var, exoid, step, EX_NODAL, offset + c, 1, xs / bs + 1, (xe - xs) / bs, varray);
1416: PetscCall(VecRestoreArray(vComp, &varray));
1417: PetscCall(VecRestoreSubVector(vNatural, compIS, &vComp));
1418: }
1419: PetscCall(ISDestroy(&compIS));
1420: }
1421: if (useNatural) {
1422: PetscCall(DMPlexNaturalToGlobalBegin(dm, vNatural, v));
1423: PetscCall(DMPlexNaturalToGlobalEnd(dm, vNatural, v));
1424: PetscCall(VecDestroy(&vNatural));
1425: }
1426: PetscFunctionReturn(PETSC_SUCCESS);
1427: }
1429: static PetscErrorCode VecViewPlex_ExodusII_Zonal_Internal(Vec v, PetscExodusIIInt exoid, PetscExodusIIInt step, PetscExodusIIInt offset)
1430: {
1431: MPI_Comm comm;
1432: PetscMPIInt size;
1433: DM dm;
1434: Vec vNatural, vComp;
1435: const PetscScalar *varray;
1436: PetscInt xs, xe, bs;
1437: PetscBool useNatural;
1438: IS compIS;
1439: PetscInt *csSize, *csID;
1440: PetscExodusIIInt numCS, set, csxs = 0;
1442: PetscFunctionBegin;
1443: PetscCall(PetscObjectGetComm((PetscObject)v, &comm));
1444: PetscCallMPI(MPI_Comm_size(comm, &size));
1445: PetscCall(VecGetDM(v, &dm));
1446: PetscCall(DMGetUseNatural(dm, &useNatural));
1447: useNatural = useNatural && size > 1 ? PETSC_TRUE : PETSC_FALSE;
1448: if (useNatural) {
1449: PetscCall(DMPlexCreateNaturalVector(dm, &vNatural));
1450: PetscCall(DMPlexGlobalToNaturalBegin(dm, v, vNatural));
1451: PetscCall(DMPlexGlobalToNaturalEnd(dm, v, vNatural));
1452: } else {
1453: vNatural = v;
1454: }
1456: /* Write local chunk of the result in the exodus file
1457: exodus stores each component of a vector-valued field as a separate variable.
1458: We assume that they are stored sequentially
1459: Zonal variables are accessed one element block at a time, so we loop through the cell sets,
1460: but once the vector has been reordered to natural size, we cannot use the label information
1461: to figure out what to save where. */
1462: numCS = ex_inquire_int(exoid, EX_INQ_ELEM_BLK);
1463: PetscCall(PetscMalloc2(numCS, &csID, numCS, &csSize));
1464: PetscCallExternal(ex_get_ids, exoid, EX_ELEM_BLOCK, csID);
1465: for (set = 0; set < numCS; ++set) {
1466: ex_block block;
1468: block.id = csID[set];
1469: block.type = EX_ELEM_BLOCK;
1470: PetscCallExternal(ex_get_block_param, exoid, &block);
1471: csSize[set] = block.num_entry;
1472: }
1473: PetscCall(VecGetOwnershipRange(vNatural, &xs, &xe));
1474: PetscCall(VecGetBlockSize(vNatural, &bs));
1475: if (bs > 1) PetscCall(ISCreateStride(comm, (xe - xs) / bs, xs, bs, &compIS));
1476: for (set = 0; set < numCS; set++) {
1477: PetscInt csLocalSize, c;
1479: /* range of indices for set setID[set]: csxs:csxs + csSize[set]-1
1480: local slice of zonal values: xs/bs,xm/bs-1
1481: intersection: max(xs/bs,csxs),min(xm/bs-1,csxs + csSize[set]-1) */
1482: csLocalSize = PetscMax(0, PetscMin(xe / bs, csxs + csSize[set]) - PetscMax(xs / bs, csxs));
1483: if (bs == 1) {
1484: PetscCall(VecGetArrayRead(vNatural, &varray));
1485: PetscCallExternal(ex_put_partial_var, exoid, step, EX_ELEM_BLOCK, offset, csID[set], PetscMax(xs - csxs, 0) + 1, csLocalSize, &varray[PetscMax(0, csxs - xs)]);
1486: PetscCall(VecRestoreArrayRead(vNatural, &varray));
1487: } else {
1488: for (c = 0; c < bs; ++c) {
1489: PetscCall(ISStrideSetStride(compIS, (xe - xs) / bs, xs + c, bs));
1490: PetscCall(VecGetSubVector(vNatural, compIS, &vComp));
1491: PetscCall(VecGetArrayRead(vComp, &varray));
1492: 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)]);
1493: PetscCall(VecRestoreArrayRead(vComp, &varray));
1494: PetscCall(VecRestoreSubVector(vNatural, compIS, &vComp));
1495: }
1496: }
1497: csxs += csSize[set];
1498: }
1499: PetscCall(PetscFree2(csID, csSize));
1500: if (bs > 1) PetscCall(ISDestroy(&compIS));
1501: if (useNatural) PetscCall(VecDestroy(&vNatural));
1502: PetscFunctionReturn(PETSC_SUCCESS);
1503: }
1505: static PetscErrorCode VecLoadPlex_ExodusII_Zonal_Internal(Vec v, PetscExodusIIInt exoid, PetscExodusIIInt step, PetscExodusIIInt offset)
1506: {
1507: MPI_Comm comm;
1508: PetscMPIInt size;
1509: DM dm;
1510: Vec vNatural, vComp;
1511: PetscScalar *varray;
1512: PetscInt xs, xe, bs;
1513: PetscBool useNatural;
1514: IS compIS;
1515: PetscInt *csSize, *csID;
1516: PetscExodusIIInt numCS, set, csxs = 0;
1518: PetscFunctionBegin;
1519: PetscCall(PetscObjectGetComm((PetscObject)v, &comm));
1520: PetscCallMPI(MPI_Comm_size(comm, &size));
1521: PetscCall(VecGetDM(v, &dm));
1522: PetscCall(DMGetUseNatural(dm, &useNatural));
1523: useNatural = useNatural && size > 1 ? PETSC_TRUE : PETSC_FALSE;
1524: if (useNatural) PetscCall(DMPlexCreateNaturalVector(dm, &vNatural));
1525: else vNatural = v;
1527: /* Read local chunk of the result in the exodus file
1528: exodus stores each component of a vector-valued field as a separate variable.
1529: We assume that they are stored sequentially
1530: Zonal variables are accessed one element block at a time, so we loop through the cell sets,
1531: but once the vector has been reordered to natural size, we cannot use the label information
1532: to figure out what to save where. */
1533: numCS = ex_inquire_int(exoid, EX_INQ_ELEM_BLK);
1534: PetscCall(PetscMalloc2(numCS, &csID, numCS, &csSize));
1535: PetscCallExternal(ex_get_ids, exoid, EX_ELEM_BLOCK, csID);
1536: for (set = 0; set < numCS; ++set) {
1537: ex_block block;
1539: block.id = csID[set];
1540: block.type = EX_ELEM_BLOCK;
1541: PetscCallExternal(ex_get_block_param, exoid, &block);
1542: csSize[set] = block.num_entry;
1543: }
1544: PetscCall(VecGetOwnershipRange(vNatural, &xs, &xe));
1545: PetscCall(VecGetBlockSize(vNatural, &bs));
1546: if (bs > 1) PetscCall(ISCreateStride(comm, (xe - xs) / bs, xs, bs, &compIS));
1547: for (set = 0; set < numCS; ++set) {
1548: PetscInt csLocalSize, c;
1550: /* range of indices for set setID[set]: csxs:csxs + csSize[set]-1
1551: local slice of zonal values: xs/bs,xm/bs-1
1552: intersection: max(xs/bs,csxs),min(xm/bs-1,csxs + csSize[set]-1) */
1553: csLocalSize = PetscMax(0, PetscMin(xe / bs, csxs + csSize[set]) - PetscMax(xs / bs, csxs));
1554: if (bs == 1) {
1555: PetscCall(VecGetArray(vNatural, &varray));
1556: PetscCallExternal(ex_get_partial_var, exoid, step, EX_ELEM_BLOCK, offset, csID[set], PetscMax(xs - csxs, 0) + 1, csLocalSize, &varray[PetscMax(0, csxs - xs)]);
1557: PetscCall(VecRestoreArray(vNatural, &varray));
1558: } else {
1559: for (c = 0; c < bs; ++c) {
1560: PetscCall(ISStrideSetStride(compIS, (xe - xs) / bs, xs + c, bs));
1561: PetscCall(VecGetSubVector(vNatural, compIS, &vComp));
1562: PetscCall(VecGetArray(vComp, &varray));
1563: 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)]);
1564: PetscCall(VecRestoreArray(vComp, &varray));
1565: PetscCall(VecRestoreSubVector(vNatural, compIS, &vComp));
1566: }
1567: }
1568: csxs += csSize[set];
1569: }
1570: PetscCall(PetscFree2(csID, csSize));
1571: if (bs > 1) PetscCall(ISDestroy(&compIS));
1572: if (useNatural) {
1573: PetscCall(DMPlexNaturalToGlobalBegin(dm, vNatural, v));
1574: PetscCall(DMPlexNaturalToGlobalEnd(dm, vNatural, v));
1575: PetscCall(VecDestroy(&vNatural));
1576: }
1577: PetscFunctionReturn(PETSC_SUCCESS);
1578: }
1580: /*@
1581: PetscViewerExodusIIGetId - Get the file id of the `PETSCVIEWEREXODUSII` file
1583: Logically Collective
1585: Input Parameter:
1586: . viewer - the `PetscViewer`
1588: Output Parameter:
1589: . exoid - The ExodusII file id
1591: Level: intermediate
1593: .seealso: `PETSCVIEWEREXODUSII`, `PetscViewer`, `PetscViewerFileSetMode()`, `PetscViewerCreate()`, `PetscViewerSetType()`, `PetscViewerBinaryOpen()`
1594: @*/
1595: PetscErrorCode PetscViewerExodusIIGetId(PetscViewer viewer, PetscExodusIIInt *exoid)
1596: {
1597: PetscFunctionBegin;
1599: PetscTryMethod(viewer, "PetscViewerGetId_C", (PetscViewer, PetscExodusIIInt *), (viewer, exoid));
1600: PetscFunctionReturn(PETSC_SUCCESS);
1601: }
1603: /*@
1604: PetscViewerExodusIISetOrder - Set the elements order in the exodusII file.
1606: Collective
1608: Input Parameters:
1609: + viewer - the `PETSCVIEWEREXODUSII` viewer
1610: - order - elements order
1612: Output Parameter:
1614: Level: beginner
1616: .seealso: `PETSCVIEWEREXODUSII`, `PetscViewer`, `PetscViewerExodusIIGetId()`, `PetscViewerExodusIIGetOrder()`
1617: @*/
1618: PetscErrorCode PetscViewerExodusIISetOrder(PetscViewer viewer, PetscInt order)
1619: {
1620: PetscFunctionBegin;
1622: PetscTryMethod(viewer, "PetscViewerSetOrder_C", (PetscViewer, PetscInt), (viewer, order));
1623: PetscFunctionReturn(PETSC_SUCCESS);
1624: }
1626: /*@
1627: PetscViewerExodusIIGetOrder - Get the elements order in the exodusII file.
1629: Collective
1631: Input Parameters:
1632: + viewer - the `PETSCVIEWEREXODUSII` viewer
1633: - order - elements order
1635: Output Parameter:
1637: Level: beginner
1639: .seealso: `PETSCVIEWEREXODUSII`, `PetscViewer`, `PetscViewerExodusIIGetId()`, `PetscViewerExodusIISetOrder()`
1640: @*/
1641: PetscErrorCode PetscViewerExodusIIGetOrder(PetscViewer viewer, PetscInt *order)
1642: {
1643: PetscFunctionBegin;
1645: PetscTryMethod(viewer, "PetscViewerGetOrder_C", (PetscViewer, PetscInt *), (viewer, order));
1646: PetscFunctionReturn(PETSC_SUCCESS);
1647: }
1649: /*@
1650: PetscViewerExodusIIOpen - Opens a file for ExodusII input/output.
1652: Collective
1654: Input Parameters:
1655: + comm - MPI communicator
1656: . name - name of file
1657: - mode - access mode
1658: .vb
1659: FILE_MODE_WRITE - create new file for binary output
1660: FILE_MODE_READ - open existing file for binary input
1661: FILE_MODE_APPEND - open existing file for binary output
1662: .ve
1664: Output Parameter:
1665: . exo - `PETSCVIEWEREXODUSII` `PetscViewer` for Exodus II input/output to use with the specified file
1667: Level: beginner
1669: .seealso: `PETSCVIEWEREXODUSII`, `PetscViewer`, `PetscViewerPushFormat()`, `PetscViewerDestroy()`,
1670: `DMLoad()`, `PetscFileMode`, `PetscViewerSetType()`, `PetscViewerFileSetMode()`, `PetscViewerFileSetName()`
1671: @*/
1672: PetscErrorCode PetscViewerExodusIIOpen(MPI_Comm comm, const char name[], PetscFileMode mode, PetscViewer *exo)
1673: {
1674: PetscFunctionBegin;
1675: PetscCall(PetscViewerCreate(comm, exo));
1676: PetscCall(PetscViewerSetType(*exo, PETSCVIEWEREXODUSII));
1677: PetscCall(PetscViewerFileSetMode(*exo, mode));
1678: PetscCall(PetscViewerFileSetName(*exo, name));
1679: PetscCall(PetscViewerSetFromOptions(*exo));
1680: PetscFunctionReturn(PETSC_SUCCESS);
1681: }
1683: static PetscErrorCode ExodusGetCellType_Internal(const char *elem_type, DMPolytopeType *ct)
1684: {
1685: PetscBool flg;
1687: PetscFunctionBegin;
1688: *ct = DM_POLYTOPE_UNKNOWN;
1689: PetscCall(PetscStrcmp(elem_type, "BAR2", &flg));
1690: if (flg) {
1691: *ct = DM_POLYTOPE_SEGMENT;
1692: goto done;
1693: }
1694: PetscCall(PetscStrcmp(elem_type, "BAR3", &flg));
1695: if (flg) {
1696: *ct = DM_POLYTOPE_SEGMENT;
1697: goto done;
1698: }
1699: PetscCall(PetscStrcmp(elem_type, "TRI", &flg));
1700: if (flg) {
1701: *ct = DM_POLYTOPE_TRIANGLE;
1702: goto done;
1703: }
1704: PetscCall(PetscStrcmp(elem_type, "TRI3", &flg));
1705: if (flg) {
1706: *ct = DM_POLYTOPE_TRIANGLE;
1707: goto done;
1708: }
1709: PetscCall(PetscStrcmp(elem_type, "QUAD", &flg));
1710: if (flg) {
1711: *ct = DM_POLYTOPE_QUADRILATERAL;
1712: goto done;
1713: }
1714: PetscCall(PetscStrcmp(elem_type, "QUAD4", &flg));
1715: if (flg) {
1716: *ct = DM_POLYTOPE_QUADRILATERAL;
1717: goto done;
1718: }
1719: PetscCall(PetscStrcmp(elem_type, "SHELL4", &flg));
1720: if (flg) {
1721: *ct = DM_POLYTOPE_QUADRILATERAL;
1722: goto done;
1723: }
1724: PetscCall(PetscStrcmp(elem_type, "TETRA", &flg));
1725: if (flg) {
1726: *ct = DM_POLYTOPE_TETRAHEDRON;
1727: goto done;
1728: }
1729: PetscCall(PetscStrcmp(elem_type, "TET4", &flg));
1730: if (flg) {
1731: *ct = DM_POLYTOPE_TETRAHEDRON;
1732: goto done;
1733: }
1734: PetscCall(PetscStrcmp(elem_type, "WEDGE", &flg));
1735: if (flg) {
1736: *ct = DM_POLYTOPE_TRI_PRISM;
1737: goto done;
1738: }
1739: PetscCall(PetscStrcmp(elem_type, "HEX", &flg));
1740: if (flg) {
1741: *ct = DM_POLYTOPE_HEXAHEDRON;
1742: goto done;
1743: }
1744: PetscCall(PetscStrcmp(elem_type, "HEX8", &flg));
1745: if (flg) {
1746: *ct = DM_POLYTOPE_HEXAHEDRON;
1747: goto done;
1748: }
1749: PetscCall(PetscStrcmp(elem_type, "HEXAHEDRON", &flg));
1750: if (flg) {
1751: *ct = DM_POLYTOPE_HEXAHEDRON;
1752: goto done;
1753: }
1754: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unrecognized element type %s", elem_type);
1755: done:
1756: PetscFunctionReturn(PETSC_SUCCESS);
1757: }
1759: /*@
1760: DMPlexCreateExodus - Create a `DMPLEX` mesh from an ExodusII file ID.
1762: Collective
1764: Input Parameters:
1765: + comm - The MPI communicator
1766: . exoid - The ExodusII id associated with a exodus file and obtained using ex_open
1767: - interpolate - Create faces and edges in the mesh
1769: Output Parameter:
1770: . dm - The `DM` object representing the mesh
1772: Level: beginner
1774: .seealso: [](ch_unstructured), `DM`, `PETSCVIEWEREXODUSII`, `DMPLEX`, `DMCreate()`
1775: @*/
1776: PetscErrorCode DMPlexCreateExodus(MPI_Comm comm, PetscExodusIIInt exoid, PetscBool interpolate, DM *dm)
1777: {
1778: PetscMPIInt num_proc, rank;
1779: DMLabel cellSets = NULL, faceSets = NULL, vertSets = NULL;
1780: PetscSection coordSection;
1781: Vec coordinates;
1782: PetscScalar *coords;
1783: PetscInt coordSize, v;
1784: /* Read from ex_get_init() */
1785: char title[PETSC_MAX_PATH_LEN + 1];
1786: int dim = 0, dimEmbed = 0, numVertices = 0, numCells = 0;
1787: int num_cs = 0, num_vs = 0, num_fs = 0;
1789: PetscFunctionBegin;
1790: PetscCallMPI(MPI_Comm_rank(comm, &rank));
1791: PetscCallMPI(MPI_Comm_size(comm, &num_proc));
1792: PetscCall(DMCreate(comm, dm));
1793: PetscCall(DMSetType(*dm, DMPLEX));
1794: /* Open EXODUS II file and read basic information on rank 0, then broadcast to all processors */
1795: if (rank == 0) {
1796: PetscCall(PetscMemzero(title, PETSC_MAX_PATH_LEN + 1));
1797: PetscCallExternal(ex_get_init, exoid, title, &dimEmbed, &numVertices, &numCells, &num_cs, &num_vs, &num_fs);
1798: PetscCheck(num_cs, PETSC_COMM_SELF, PETSC_ERR_SUP, "Exodus file does not contain any cell set");
1799: }
1800: PetscCallMPI(MPI_Bcast(title, PETSC_MAX_PATH_LEN + 1, MPI_CHAR, 0, comm));
1801: PetscCallMPI(MPI_Bcast(&dim, 1, MPI_INT, 0, comm));
1802: PetscCall(PetscObjectSetName((PetscObject)*dm, title));
1803: PetscCall(DMPlexSetChart(*dm, 0, numCells + numVertices));
1804: /* We do not want this label automatically computed, instead we compute it here */
1805: PetscCall(DMCreateLabel(*dm, "celltype"));
1807: /* Read cell sets information */
1808: if (rank == 0) {
1809: PetscInt *cone;
1810: int c, cs, ncs, c_loc, v, v_loc;
1811: /* Read from ex_get_elem_blk_ids() */
1812: int *cs_id, *cs_order;
1813: /* Read from ex_get_elem_block() */
1814: char buffer[PETSC_MAX_PATH_LEN + 1];
1815: int num_cell_in_set, num_vertex_per_cell, num_hybrid, num_attr;
1816: /* Read from ex_get_elem_conn() */
1817: int *cs_connect;
1819: /* Get cell sets IDs */
1820: PetscCall(PetscMalloc2(num_cs, &cs_id, num_cs, &cs_order));
1821: PetscCallExternal(ex_get_ids, exoid, EX_ELEM_BLOCK, cs_id);
1822: /* Read the cell set connectivity table and build mesh topology
1823: EXO standard requires that cells in cell sets be numbered sequentially and be pairwise disjoint. */
1824: /* Check for a hybrid mesh */
1825: for (cs = 0, num_hybrid = 0; cs < num_cs; ++cs) {
1826: DMPolytopeType ct;
1827: char elem_type[PETSC_MAX_PATH_LEN];
1829: PetscCall(PetscArrayzero(elem_type, sizeof(elem_type)));
1830: PetscCallExternal(ex_get_elem_type, exoid, cs_id[cs], elem_type);
1831: PetscCall(ExodusGetCellType_Internal(elem_type, &ct));
1832: dim = PetscMax(dim, DMPolytopeTypeGetDim(ct));
1833: PetscCallExternal(ex_get_block, exoid, EX_ELEM_BLOCK, cs_id[cs], buffer, &num_cell_in_set, &num_vertex_per_cell, 0, 0, &num_attr);
1834: switch (ct) {
1835: case DM_POLYTOPE_TRI_PRISM:
1836: cs_order[cs] = cs;
1837: ++num_hybrid;
1838: break;
1839: default:
1840: for (c = cs; c > cs - num_hybrid; --c) cs_order[c] = cs_order[c - 1];
1841: cs_order[cs - num_hybrid] = cs;
1842: }
1843: }
1844: /* First set sizes */
1845: for (ncs = 0, c = 0; ncs < num_cs; ++ncs) {
1846: DMPolytopeType ct;
1847: char elem_type[PETSC_MAX_PATH_LEN];
1848: const PetscInt cs = cs_order[ncs];
1850: PetscCall(PetscArrayzero(elem_type, sizeof(elem_type)));
1851: PetscCallExternal(ex_get_elem_type, exoid, cs_id[cs], elem_type);
1852: PetscCall(ExodusGetCellType_Internal(elem_type, &ct));
1853: PetscCallExternal(ex_get_block, exoid, EX_ELEM_BLOCK, cs_id[cs], buffer, &num_cell_in_set, &num_vertex_per_cell, 0, 0, &num_attr);
1854: for (c_loc = 0; c_loc < num_cell_in_set; ++c_loc, ++c) {
1855: PetscCall(DMPlexSetConeSize(*dm, c, num_vertex_per_cell));
1856: PetscCall(DMPlexSetCellType(*dm, c, ct));
1857: }
1858: }
1859: for (v = numCells; v < numCells + numVertices; ++v) PetscCall(DMPlexSetCellType(*dm, v, DM_POLYTOPE_POINT));
1860: PetscCall(DMSetUp(*dm));
1861: for (ncs = 0, c = 0; ncs < num_cs; ++ncs) {
1862: const PetscInt cs = cs_order[ncs];
1863: PetscCallExternal(ex_get_block, exoid, EX_ELEM_BLOCK, cs_id[cs], buffer, &num_cell_in_set, &num_vertex_per_cell, 0, 0, &num_attr);
1864: PetscCall(PetscMalloc2(num_vertex_per_cell * num_cell_in_set, &cs_connect, num_vertex_per_cell, &cone));
1865: PetscCallExternal(ex_get_conn, exoid, EX_ELEM_BLOCK, cs_id[cs], cs_connect, NULL, NULL);
1866: /* EXO uses Fortran-based indexing, DMPlex uses C-style and numbers cell first then vertices. */
1867: for (c_loc = 0, v = 0; c_loc < num_cell_in_set; ++c_loc, ++c) {
1868: DMPolytopeType ct;
1870: for (v_loc = 0; v_loc < num_vertex_per_cell; ++v_loc, ++v) cone[v_loc] = cs_connect[v] + numCells - 1;
1871: PetscCall(DMPlexGetCellType(*dm, c, &ct));
1872: PetscCall(DMPlexInvertCell(ct, cone));
1873: PetscCall(DMPlexSetCone(*dm, c, cone));
1874: PetscCall(DMSetLabelValue_Fast(*dm, &cellSets, "Cell Sets", c, cs_id[cs]));
1875: }
1876: PetscCall(PetscFree2(cs_connect, cone));
1877: }
1878: PetscCall(PetscFree2(cs_id, cs_order));
1879: }
1880: {
1881: PetscInt ints[] = {dim, dimEmbed};
1883: PetscCallMPI(MPI_Bcast(ints, 2, MPIU_INT, 0, comm));
1884: PetscCall(DMSetDimension(*dm, ints[0]));
1885: PetscCall(DMSetCoordinateDim(*dm, ints[1]));
1886: dim = ints[0];
1887: dimEmbed = ints[1];
1888: }
1889: PetscCall(DMPlexSymmetrize(*dm));
1890: PetscCall(DMPlexStratify(*dm));
1891: if (interpolate) {
1892: DM idm;
1894: PetscCall(DMPlexInterpolate(*dm, &idm));
1895: PetscCall(DMDestroy(dm));
1896: *dm = idm;
1897: }
1899: /* Create vertex set label */
1900: if (rank == 0 && (num_vs > 0)) {
1901: int vs, v;
1902: /* Read from ex_get_node_set_ids() */
1903: int *vs_id;
1904: /* Read from ex_get_node_set_param() */
1905: int num_vertex_in_set;
1906: /* Read from ex_get_node_set() */
1907: int *vs_vertex_list;
1909: /* Get vertex set ids */
1910: PetscCall(PetscMalloc1(num_vs, &vs_id));
1911: PetscCallExternal(ex_get_ids, exoid, EX_NODE_SET, vs_id);
1912: for (vs = 0; vs < num_vs; ++vs) {
1913: PetscCallExternal(ex_get_set_param, exoid, EX_NODE_SET, vs_id[vs], &num_vertex_in_set, NULL);
1914: PetscCall(PetscMalloc1(num_vertex_in_set, &vs_vertex_list));
1915: PetscCallExternal(ex_get_set, exoid, EX_NODE_SET, vs_id[vs], vs_vertex_list, NULL);
1916: for (v = 0; v < num_vertex_in_set; ++v) PetscCall(DMSetLabelValue_Fast(*dm, &vertSets, "Vertex Sets", vs_vertex_list[v] + numCells - 1, vs_id[vs]));
1917: PetscCall(PetscFree(vs_vertex_list));
1918: }
1919: PetscCall(PetscFree(vs_id));
1920: }
1921: /* Read coordinates */
1922: PetscCall(DMGetCoordinateSection(*dm, &coordSection));
1923: PetscCall(PetscSectionSetNumFields(coordSection, 1));
1924: PetscCall(PetscSectionSetFieldComponents(coordSection, 0, dimEmbed));
1925: PetscCall(PetscSectionSetChart(coordSection, numCells, numCells + numVertices));
1926: for (v = numCells; v < numCells + numVertices; ++v) {
1927: PetscCall(PetscSectionSetDof(coordSection, v, dimEmbed));
1928: PetscCall(PetscSectionSetFieldDof(coordSection, v, 0, dimEmbed));
1929: }
1930: PetscCall(PetscSectionSetUp(coordSection));
1931: PetscCall(PetscSectionGetStorageSize(coordSection, &coordSize));
1932: PetscCall(VecCreate(PETSC_COMM_SELF, &coordinates));
1933: PetscCall(PetscObjectSetName((PetscObject)coordinates, "coordinates"));
1934: PetscCall(VecSetSizes(coordinates, coordSize, PETSC_DETERMINE));
1935: PetscCall(VecSetBlockSize(coordinates, dimEmbed));
1936: PetscCall(VecSetType(coordinates, VECSTANDARD));
1937: PetscCall(VecGetArray(coordinates, &coords));
1938: if (rank == 0) {
1939: PetscReal *x, *y, *z;
1941: PetscCall(PetscMalloc3(numVertices, &x, numVertices, &y, numVertices, &z));
1942: PetscCallExternal(ex_get_coord, exoid, x, y, z);
1943: if (dimEmbed > 0) {
1944: for (v = 0; v < numVertices; ++v) coords[v * dimEmbed + 0] = x[v];
1945: }
1946: if (dimEmbed > 1) {
1947: for (v = 0; v < numVertices; ++v) coords[v * dimEmbed + 1] = y[v];
1948: }
1949: if (dimEmbed > 2) {
1950: for (v = 0; v < numVertices; ++v) coords[v * dimEmbed + 2] = z[v];
1951: }
1952: PetscCall(PetscFree3(x, y, z));
1953: }
1954: PetscCall(VecRestoreArray(coordinates, &coords));
1955: PetscCall(DMSetCoordinatesLocal(*dm, coordinates));
1956: PetscCall(VecDestroy(&coordinates));
1958: /* Create side set label */
1959: if (rank == 0 && interpolate && (num_fs > 0)) {
1960: int fs, f, voff;
1961: /* Read from ex_get_side_set_ids() */
1962: int *fs_id;
1963: /* Read from ex_get_side_set_param() */
1964: int num_side_in_set;
1965: /* Read from ex_get_side_set_node_list() */
1966: int *fs_vertex_count_list, *fs_vertex_list, *fs_side_list;
1967: /* Read side set labels */
1968: char fs_name[MAX_STR_LENGTH + 1];
1969: size_t fs_name_len;
1971: /* Get side set ids */
1972: PetscCall(PetscMalloc1(num_fs, &fs_id));
1973: PetscCallExternal(ex_get_ids, exoid, EX_SIDE_SET, fs_id);
1974: // Ids 1 and 2 are reserved by ExodusII for indicating things in 3D
1975: for (fs = 0; fs < num_fs; ++fs) {
1976: PetscCallExternal(ex_get_set_param, exoid, EX_SIDE_SET, fs_id[fs], &num_side_in_set, NULL);
1977: PetscCall(PetscMalloc3(num_side_in_set, &fs_vertex_count_list, num_side_in_set * 4, &fs_vertex_list, num_side_in_set, &fs_side_list));
1978: PetscCallExternal(ex_get_side_set_node_list, exoid, fs_id[fs], fs_vertex_count_list, fs_vertex_list);
1979: PetscCallExternal(ex_get_set, exoid, EX_SIDE_SET, fs_id[fs], NULL, fs_side_list);
1981: /* Get the specific name associated with this side set ID. */
1982: int fs_name_err = ex_get_name(exoid, EX_SIDE_SET, fs_id[fs], fs_name);
1983: if (!fs_name_err) {
1984: PetscCall(PetscStrlen(fs_name, &fs_name_len));
1985: if (fs_name_len == 0) PetscCall(PetscStrncpy(fs_name, "Face Sets", MAX_STR_LENGTH));
1986: }
1987: for (f = 0, voff = 0; f < num_side_in_set; ++f) {
1988: const PetscInt *faces = NULL;
1989: PetscInt faceSize = fs_vertex_count_list[f], numFaces;
1990: PetscInt faceVertices[4], v;
1992: PetscCheck(faceSize <= 4, comm, PETSC_ERR_ARG_WRONG, "ExodusII side cannot have %" PetscInt_FMT " > 4 vertices", faceSize);
1993: for (v = 0; v < faceSize; ++v, ++voff) faceVertices[v] = fs_vertex_list[voff] + numCells - 1;
1994: PetscCall(DMPlexGetFullJoin(*dm, faceSize, faceVertices, &numFaces, &faces));
1995: PetscCheck(numFaces == 1, comm, PETSC_ERR_ARG_WRONG, "Invalid ExodusII side %d in set %d maps to %" PetscInt_FMT " faces", f, fs, numFaces);
1996: PetscCheck(dim == 1 || faces[0] >= numCells + numVertices, comm, PETSC_ERR_ARG_WRONG, "Invalid ExodusII side %d in set %d maps to point %" PetscInt_FMT " which is not a face", f, fs, faces[0]);
1997: PetscCall(DMSetLabelValue_Fast(*dm, &faceSets, "Face Sets", faces[0], fs_id[fs]));
1998: /* Only add the label if one has been detected for this side set. */
1999: if (!fs_name_err) PetscCall(DMSetLabelValue(*dm, fs_name, faces[0], fs_id[fs]));
2000: PetscCall(DMPlexRestoreJoin(*dm, faceSize, faceVertices, &numFaces, &faces));
2001: }
2002: PetscCall(PetscFree3(fs_vertex_count_list, fs_vertex_list, fs_side_list));
2003: }
2004: PetscCall(PetscFree(fs_id));
2005: }
2007: { /* Create Cell/Face/Vertex Sets labels at all processes */
2008: enum {
2009: n = 3
2010: };
2011: PetscBool flag[n];
2013: flag[0] = cellSets ? PETSC_TRUE : PETSC_FALSE;
2014: flag[1] = faceSets ? PETSC_TRUE : PETSC_FALSE;
2015: flag[2] = vertSets ? PETSC_TRUE : PETSC_FALSE;
2016: PetscCallMPI(MPI_Bcast(flag, n, MPIU_BOOL, 0, comm));
2017: if (flag[0]) PetscCall(DMCreateLabel(*dm, "Cell Sets"));
2018: if (flag[1]) PetscCall(DMCreateLabel(*dm, "Face Sets"));
2019: if (flag[2]) PetscCall(DMCreateLabel(*dm, "Vertex Sets"));
2020: }
2021: PetscFunctionReturn(PETSC_SUCCESS);
2022: }