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