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