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