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