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: PetscMPIInt rank;
96: PetscExodusIIInt CPU_word_size, IO_word_size, EXO_mode;
97: MPI_Info mpi_info = MPI_INFO_NULL;
98: PetscExodusIIFloat EXO_version;
100: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)viewer), &rank));
101: CPU_word_size = sizeof(PetscReal);
102: IO_word_size = sizeof(PetscReal);
104: PetscFunctionBegin;
105: if (exo->exoid >= 0) {
106: PetscCallExternal(ex_close, exo->exoid);
107: exo->exoid = -1;
108: }
109: if (exo->filename) PetscCall(PetscFree(exo->filename));
110: PetscCall(PetscStrallocpy(name, &exo->filename));
111: switch (exo->btype) {
112: case FILE_MODE_READ:
113: EXO_mode = EX_READ;
114: break;
115: case FILE_MODE_APPEND:
116: case FILE_MODE_UPDATE:
117: case FILE_MODE_APPEND_UPDATE:
118: /* Will fail if the file does not already exist */
119: EXO_mode = EX_WRITE;
120: break;
121: case FILE_MODE_WRITE:
122: /*
123: exodus only allows writing geometry upon file creation, so we will let DMView create the file.
124: */
125: PetscFunctionReturn(PETSC_SUCCESS);
126: break;
127: default:
128: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ORDER, "Must call PetscViewerFileSetMode() before PetscViewerFileSetName()");
129: }
130: #if defined(PETSC_USE_64BIT_INDICES)
131: EXO_mode += EX_ALL_INT64_API;
132: #endif
133: exo->exoid = ex_open_par(name, EXO_mode, &CPU_word_size, &IO_word_size, &EXO_version, PETSC_COMM_WORLD, mpi_info);
134: PetscCheck(exo->exoid >= 0, PETSC_COMM_SELF, PETSC_ERR_LIB, "ex_open_par failed for %s", name);
135: PetscFunctionReturn(PETSC_SUCCESS);
136: }
138: static PetscErrorCode PetscViewerFileGetName_ExodusII(PetscViewer viewer, const char *name[])
139: {
140: PetscViewer_ExodusII *exo = (PetscViewer_ExodusII *)viewer->data;
142: PetscFunctionBegin;
143: *name = exo->filename;
144: PetscFunctionReturn(PETSC_SUCCESS);
145: }
147: static PetscErrorCode PetscViewerFileSetMode_ExodusII(PetscViewer viewer, PetscFileMode type)
148: {
149: PetscViewer_ExodusII *exo = (PetscViewer_ExodusII *)viewer->data;
151: PetscFunctionBegin;
152: exo->btype = type;
153: PetscFunctionReturn(PETSC_SUCCESS);
154: }
156: static PetscErrorCode PetscViewerFileGetMode_ExodusII(PetscViewer viewer, PetscFileMode *type)
157: {
158: PetscViewer_ExodusII *exo = (PetscViewer_ExodusII *)viewer->data;
160: PetscFunctionBegin;
161: *type = exo->btype;
162: PetscFunctionReturn(PETSC_SUCCESS);
163: }
165: static PetscErrorCode PetscViewerExodusIIGetId_ExodusII(PetscViewer viewer, PetscExodusIIInt *exoid)
166: {
167: PetscViewer_ExodusII *exo = (PetscViewer_ExodusII *)viewer->data;
169: PetscFunctionBegin;
170: *exoid = exo->exoid;
171: PetscFunctionReturn(PETSC_SUCCESS);
172: }
174: static PetscErrorCode PetscViewerExodusIIGetOrder_ExodusII(PetscViewer viewer, PetscInt *order)
175: {
176: PetscViewer_ExodusII *exo = (PetscViewer_ExodusII *)viewer->data;
178: PetscFunctionBegin;
179: *order = exo->order;
180: PetscFunctionReturn(PETSC_SUCCESS);
181: }
183: static PetscErrorCode PetscViewerExodusIISetOrder_ExodusII(PetscViewer viewer, PetscInt order)
184: {
185: PetscViewer_ExodusII *exo = (PetscViewer_ExodusII *)viewer->data;
187: PetscFunctionBegin;
188: exo->order = order;
189: PetscFunctionReturn(PETSC_SUCCESS);
190: }
192: /*@
193: PetscViewerExodusIISetZonalVariable - Sets the number of zonal variables in an ExodusII file
195: Collective;
197: Input Parameters:
198: + viewer - a `PetscViewer` of type `PETSCVIEWEREXODUSII`
199: - num - the number of zonal variables in the ExodusII file
201: Level: intermediate
203: Notes:
204: The ExodusII API does not allow changing the number of variables in a file so this function will return an error
205: if called twice, called on a read-only file, or called on file for which the number of variables has already been specified
207: .seealso: `PETSCVIEWEREXODUSII`, `PetscViewer`, `PetscViewerCreate()`, `PetscViewerDestroy()`, `PetscViewerExodusIIOpen()`, `PetscViewerSetType()`, `PetscViewerType`, `PetscViewerExodusIIGetZonalVariable()`
208: @*/
209: PetscErrorCode PetscViewerExodusIISetZonalVariable(PetscViewer viewer, PetscExodusIIInt num)
210: {
211: PetscViewer_ExodusII *exo = (PetscViewer_ExodusII *)viewer->data;
212: MPI_Comm comm;
213: PetscExodusIIInt exoid = -1;
215: PetscFunctionBegin;
216: PetscCall(PetscObjectGetComm((PetscObject)viewer, &comm));
217: 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);
218: 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");
220: exo->numZonalVariables = num;
221: PetscCall(PetscMalloc1(num, &exo->zonalVariableNames));
222: for (int i = 0; i < num; i++) exo->zonalVariableNames[i] = NULL;
223: PetscCall(PetscViewerExodusIIGetId(viewer, &exoid));
224: PetscCallExternal(ex_put_variable_param, exoid, EX_ELEM_BLOCK, num);
225: PetscFunctionReturn(PETSC_SUCCESS);
226: }
228: /*@
229: PetscViewerExodusIISetNodalVariable - Sets the number of nodal variables in an ExodusII file
231: Collective;
233: Input Parameters:
234: + viewer - a `PetscViewer` of type `PETSCVIEWEREXODUSII`
235: - num - the number of nodal variables in the ExodusII file
237: Level: intermediate
239: Notes:
240: The ExodusII API does not allow changing the number of variables in a file so this function will return an error
241: if called twice, called on a read-only file, or called on file for which the number of variables has already been specified
243: .seealso: `PETSCVIEWEREXODUSII`, `PetscViewer`, `PetscViewerCreate()`, `PetscViewerDestroy()`, `PetscViewerExodusIIOpen()`, `PetscViewerSetType()`, `PetscViewerType`, `PetscViewerExodusIIGetNodalVariable()`
244: @*/
245: PetscErrorCode PetscViewerExodusIISetNodalVariable(PetscViewer viewer, PetscExodusIIInt num)
246: {
247: PetscViewer_ExodusII *exo = (PetscViewer_ExodusII *)viewer->data;
248: MPI_Comm comm;
249: PetscExodusIIInt exoid = -1;
251: PetscFunctionBegin;
252: PetscCall(PetscObjectGetComm((PetscObject)viewer, &comm));
253: 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);
254: 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");
256: exo->numNodalVariables = num;
257: PetscCall(PetscMalloc1(num, &exo->nodalVariableNames));
258: for (int i = 0; i < num; i++) exo->nodalVariableNames[i] = NULL;
259: PetscCall(PetscViewerExodusIIGetId(viewer, &exoid));
260: PetscCallExternal(ex_put_variable_param, exoid, EX_NODAL, num);
261: PetscFunctionReturn(PETSC_SUCCESS);
262: }
264: /*@
265: PetscViewerExodusIIGetZonalVariable - Gets the number of zonal variables in an ExodusII file
267: Collective
269: Input Parameters:
270: . viewer - a `PetscViewer` of type `PETSCVIEWEREXODUSII`
272: Output Parameter:
273: . num - the number variables in the ExodusII file
275: Level: intermediate
277: Notes:
278: The number of variables in the ExodusII file is cached in the viewer
280: .seealso: `PETSCVIEWEREXODUSII`, `PetscViewer`, `PetscViewerCreate()`, `PetscViewerDestroy()`, `PetscViewerExodusIIOpen()`, `PetscViewerSetType()`, `PetscViewerType`, `PetscViewerExodusIIsetZonalVariable()`
281: @*/
282: PetscErrorCode PetscViewerExodusIIGetZonalVariable(PetscViewer viewer, PetscExodusIIInt *num)
283: {
284: PetscViewer_ExodusII *exo = (PetscViewer_ExodusII *)viewer->data;
285: MPI_Comm comm;
286: PetscExodusIIInt exoid = -1;
288: PetscFunctionBegin;
289: if (exo->numZonalVariables > -1) {
290: *num = exo->numZonalVariables;
291: } else {
292: PetscCall(PetscViewerExodusIIGetId(viewer, &exoid));
293: PetscCall(PetscObjectGetComm((PetscObject)viewer, &comm));
294: PetscCheck(exoid > 0, comm, PETSC_ERR_FILE_OPEN, "Exodus file is not open");
295: PetscCallExternal(ex_get_variable_param, exoid, EX_ELEM_BLOCK, num);
296: exo->numZonalVariables = *num;
297: PetscCall(PetscMalloc1(*num, &exo->zonalVariableNames));
298: for (int i = 0; i < *num; i++) exo->zonalVariableNames[i] = NULL;
299: }
300: PetscFunctionReturn(PETSC_SUCCESS);
301: }
303: /*@
304: PetscViewerExodusIIGetNodalVariable - Gets the number of nodal variables in an ExodusII file
306: Collective
308: Input Parameters:
309: . viewer - a `PetscViewer` of type `PETSCVIEWEREXODUSII`
311: Output Parameter:
312: . num - the number variables in the ExodusII file
314: Level: intermediate
316: Notes:
317: This function gets the number of nodal variables and saves it in the address of num.
319: .seealso: `PETSCVIEWEREXODUSII`, `PetscViewer`, `PetscViewerCreate()`, `PetscViewerDestroy()`, `PetscViewerExodusIIOpen()`, `PetscViewerSetType()`, `PetscViewerType`, `PetscViewerExodusIISetNodalVariable()`
320: @*/
321: PetscErrorCode PetscViewerExodusIIGetNodalVariable(PetscViewer viewer, PetscExodusIIInt *num)
322: {
323: PetscViewer_ExodusII *exo = (PetscViewer_ExodusII *)viewer->data;
324: MPI_Comm comm;
325: PetscExodusIIInt exoid = -1;
327: PetscFunctionBegin;
328: if (exo->numNodalVariables > -1) {
329: *num = exo->numNodalVariables;
330: } else {
331: PetscCall(PetscViewerExodusIIGetId(viewer, &exoid));
332: PetscCall(PetscObjectGetComm((PetscObject)viewer, &comm));
333: PetscCheck(exoid > 0, comm, PETSC_ERR_FILE_OPEN, "Exodus file is not open");
334: PetscCallExternal(ex_get_variable_param, exoid, EX_NODAL, num);
335: exo->numNodalVariables = *num;
336: PetscCall(PetscMalloc1(*num, &exo->nodalVariableNames));
337: for (int i = 0; i < *num; i++) exo->nodalVariableNames[i] = NULL;
338: }
339: PetscFunctionReturn(PETSC_SUCCESS);
340: }
342: /*@
343: PetscViewerExodusIISetZonalVariableName - Sets the name of a zonal variable.
345: Collective;
347: Input Parameters:
348: + viewer - a `PetscViewer` of type `PETSCVIEWEREXODUSII`
349: . idx - the index for which you want to save the name
350: - name - string containing the name characters
352: Level: intermediate
354: .seealso: `PETSCVIEWEREXODUSII`, `PetscViewer`, `PetscViewerCreate()`, `PetscViewerDestroy()`, `PetscViewerExodusIIOpen()`, `PetscViewerSetType()`, `PetscViewerType`, `PetscViewerExodusIIGetZonalVariableName()`
355: @*/
356: PetscErrorCode PetscViewerExodusIISetZonalVariableName(PetscViewer viewer, PetscExodusIIInt idx, const char name[])
357: {
358: PetscViewer_ExodusII *exo = (PetscViewer_ExodusII *)viewer->data;
360: PetscFunctionBegin;
361: PetscCheck((idx >= 0) && (idx < exo->numZonalVariables), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Variable index out of range. Was PetscViewerExodusIISetZonalVariable called?");
362: PetscCall(PetscStrallocpy(name, (char **)&exo->zonalVariableNames[idx]));
363: PetscCallExternal(ex_put_variable_name, exo->exoid, EX_ELEM_BLOCK, idx + 1, name);
364: PetscFunctionReturn(PETSC_SUCCESS);
365: }
367: /*@
368: PetscViewerExodusIISetNodalVariableName - Sets the name of a nodal variable.
370: Collective;
372: Input Parameters:
373: + viewer - a `PetscViewer` of type `PETSCVIEWEREXODUSII`
374: . idx - the index for which you want to save the name
375: - name - string containing the name characters
377: Level: intermediate
379: .seealso: `PETSCVIEWEREXODUSII`, `PetscViewer`, `PetscViewerCreate()`, `PetscViewerDestroy()`, `PetscViewerExodusIIOpen()`, `PetscViewerSetType()`, `PetscViewerType`, `PetscViewerExodusIIGetNodalVariableName()`
380: @*/
381: PetscErrorCode PetscViewerExodusIISetNodalVariableName(PetscViewer viewer, PetscExodusIIInt idx, const char name[])
382: {
383: PetscViewer_ExodusII *exo = (PetscViewer_ExodusII *)viewer->data;
385: PetscFunctionBegin;
386: PetscCheck((idx >= 0) && (idx < exo->numNodalVariables), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Variable index out of range. Was PetscViewerExodusIISetNodalVariable called?");
387: PetscCall(PetscStrallocpy(name, (char **)&exo->nodalVariableNames[idx]));
388: PetscCallExternal(ex_put_variable_name, exo->exoid, EX_NODAL, idx + 1, name);
389: PetscFunctionReturn(PETSC_SUCCESS);
390: }
392: /*@
393: PetscViewerExodusIIGetZonalVariableName - Gets the name of a zonal variable.
395: Collective;
397: Input Parameters:
398: + viewer - a `PetscViewer` of type `PETSCVIEWEREXODUSII`
399: - idx - the index for which you want to get the name
401: Output Parameter:
402: . name - pointer to the string containing the name characters
404: Level: intermediate
406: .seealso: `PETSCVIEWEREXODUSII`, `PetscViewer`, `PetscViewerCreate()`, `PetscViewerDestroy()`, `PetscViewerExodusIIOpen()`, `PetscViewerSetType()`, `PetscViewerType`, `PetscViewerExodusIISetZonalVariableName()`
407: @*/
408: PetscErrorCode PetscViewerExodusIIGetZonalVariableName(PetscViewer viewer, PetscExodusIIInt idx, const char *name[])
409: {
410: PetscViewer_ExodusII *exo = (PetscViewer_ExodusII *)viewer->data;
411: PetscExodusIIInt exoid = -1;
412: char tmpName[MAX_NAME_LENGTH + 1];
414: PetscFunctionBegin;
415: PetscCheck(idx >= 0 && idx < exo->numZonalVariables, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Variable index out of range. Was PetscViewerExodusIISetZonalVariable called?");
416: if (!exo->zonalVariableNames[idx]) {
417: PetscCall(PetscViewerExodusIIGetId(viewer, &exoid));
418: PetscCallExternal(ex_get_variable_name, exoid, EX_ELEM_BLOCK, idx + 1, tmpName);
419: PetscCall(PetscStrallocpy(tmpName, (char **)&exo->zonalVariableNames[idx]));
420: }
421: *name = exo->zonalVariableNames[idx];
422: PetscFunctionReturn(PETSC_SUCCESS);
423: }
425: /*@
426: PetscViewerExodusIIGetNodalVariableName - Gets the name of a nodal variable.
428: Collective;
430: Input Parameters:
431: + viewer - a `PetscViewer` of type `PETSCVIEWEREXODUSII`
432: - idx - the index for which you want to save the name
434: Output Parameter:
435: . name - string array containing name characters
437: Level: intermediate
439: .seealso: `PETSCVIEWEREXODUSII`, `PetscViewer`, `PetscViewerCreate()`, `PetscViewerDestroy()`, `PetscViewerExodusIIOpen()`, `PetscViewerSetType()`, `PetscViewerType`, `PetscViewerExodusIISetNodalVariableName()`
440: @*/
441: PetscErrorCode PetscViewerExodusIIGetNodalVariableName(PetscViewer viewer, PetscExodusIIInt idx, const char *name[])
442: {
443: PetscViewer_ExodusII *exo = (PetscViewer_ExodusII *)viewer->data;
444: PetscExodusIIInt exoid = -1;
445: char tmpName[MAX_NAME_LENGTH + 1];
447: PetscFunctionBegin;
448: PetscCheck((idx >= 0) && (idx < exo->numNodalVariables), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Variable index out of range. Was PetscViewerExodusIISetNodalVariable called?");
449: if (!exo->nodalVariableNames[idx]) {
450: PetscCall(PetscViewerExodusIIGetId(viewer, &exoid));
451: PetscCallExternal(ex_get_variable_name, exoid, EX_NODAL, idx + 1, tmpName);
452: PetscCall(PetscStrallocpy(tmpName, (char **)&exo->nodalVariableNames[idx]));
453: }
454: *name = exo->nodalVariableNames[idx];
455: PetscFunctionReturn(PETSC_SUCCESS);
456: }
458: /*@C
459: PetscViewerExodusIISetZonalVariableNames - Sets the names of all nodal variables
461: Collective; No Fortran Support
463: Input Parameters:
464: + viewer - a `PetscViewer` of type `PETSCVIEWEREXODUSII`
465: - names - an array of string names to be set, the strings are copied into the `PetscViewer`
467: Level: intermediate
469: Notes:
470: This function allows users to set multiple zonal variable names at a time.
472: .seealso: `PETSCVIEWEREXODUSII`, `PetscViewer`, `PetscViewerCreate()`, `PetscViewerDestroy()`, `PetscViewerExodusIIOpen()`, `PetscViewerSetType()`, `PetscViewerType`, `PetscViewerExodusIIGetZonalVariableNames()`
473: @*/
474: PetscErrorCode PetscViewerExodusIISetZonalVariableNames(PetscViewer viewer, const char *const names[])
475: {
476: PetscExodusIIInt numNames;
477: PetscExodusIIInt exoid = -1;
478: PetscViewer_ExodusII *exo = (PetscViewer_ExodusII *)viewer->data;
480: PetscFunctionBegin;
481: PetscCall(PetscViewerExodusIIGetZonalVariable(viewer, &numNames));
482: PetscCheck(numNames >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Number of zonal variables not set. Was PetscViewerExodusIISetZonalVariable called?");
484: PetscCall(PetscViewerExodusIIGetId(viewer, &exoid));
485: for (int i = 0; i < numNames; i++) {
486: PetscCall(PetscStrallocpy(names[i], &exo->zonalVariableNames[i]));
487: PetscCallExternal(ex_put_variable_name, exoid, EX_ELEM_BLOCK, i + 1, names[i]);
488: }
489: PetscFunctionReturn(PETSC_SUCCESS);
490: }
492: /*@C
493: PetscViewerExodusIISetNodalVariableNames - Sets the names of all nodal variables.
495: Collective; No Fortran Support
497: Input Parameters:
498: + viewer - a `PetscViewer` of type `PETSCVIEWEREXODUSII`
499: - names - an array of string names to be set, the strings are copied into the `PetscViewer`
501: Level: intermediate
503: Notes:
504: This function allows users to set multiple nodal variable names at a time.
506: .seealso: `PETSCVIEWEREXODUSII`, `PetscViewer`, `PetscViewerCreate()`, `PetscViewerDestroy()`, `PetscViewerExodusIIOpen()`, `PetscViewerSetType()`, `PetscViewerType`, `PetscViewerExodusIIGetNodalVariableNames()`
507: @*/
508: PetscErrorCode PetscViewerExodusIISetNodalVariableNames(PetscViewer viewer, const char *const names[])
509: {
510: PetscExodusIIInt numNames;
511: PetscExodusIIInt exoid = -1;
512: PetscViewer_ExodusII *exo = (PetscViewer_ExodusII *)viewer->data;
514: PetscFunctionBegin;
515: PetscCall(PetscViewerExodusIIGetNodalVariable(viewer, &numNames));
516: PetscCheck(numNames >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Number of nodal variables not set. Was PetscViewerExodusIISetNodalVariable called?");
518: PetscCall(PetscViewerExodusIIGetId(viewer, &exoid));
519: for (int i = 0; i < numNames; i++) {
520: PetscCall(PetscStrallocpy(names[i], &exo->nodalVariableNames[i]));
521: PetscCallExternal(ex_put_variable_name, exoid, EX_NODAL, i + 1, names[i]);
522: }
523: PetscFunctionReturn(PETSC_SUCCESS);
524: }
526: /*@C
527: PetscViewerExodusIIGetZonalVariableNames - Gets the names of all zonal variables.
529: Collective; No Fortran Support
531: Input Parameters:
532: + viewer - a `PetscViewer` of type `PETSCVIEWEREXODUSII`
533: - numVars - the number of zonal variable names to retrieve
535: Output Parameter:
536: . varNames - returns an array of char pointers where the zonal variable names are
538: Level: intermediate
540: Notes:
541: This function returns a borrowed pointer which should not be freed.
543: .seealso: `PETSCVIEWEREXODUSII`, `PetscViewer`, `PetscViewerCreate()`, `PetscViewerDestroy()`, `PetscViewerExodusIIOpen()`, `PetscViewerSetType()`, `PetscViewerType`, `PetscViewerExodusIISetZonalVariableNames()`
544: @*/
545: PetscErrorCode PetscViewerExodusIIGetZonalVariableNames(PetscViewer viewer, PetscExodusIIInt *numVars, const char *const *varNames[])
546: {
547: PetscViewer_ExodusII *exo = (PetscViewer_ExodusII *)viewer->data;
548: PetscExodusIIInt idx;
549: char tmpName[MAX_NAME_LENGTH + 1];
550: PetscExodusIIInt exoid = -1;
552: PetscFunctionBegin;
553: PetscCall(PetscViewerExodusIIGetZonalVariable(viewer, numVars));
554: /*
555: Cache variable names if necessary
556: */
557: for (idx = 0; idx < *numVars; idx++) {
558: if (!exo->zonalVariableNames[idx]) {
559: PetscCall(PetscViewerExodusIIGetId(viewer, &exoid));
560: PetscCallExternal(ex_get_variable_name, exoid, EX_ELEM_BLOCK, idx + 1, tmpName);
561: PetscCall(PetscStrallocpy(tmpName, (char **)&exo->zonalVariableNames[idx]));
562: }
563: }
564: *varNames = (const char *const *)exo->zonalVariableNames;
565: PetscFunctionReturn(PETSC_SUCCESS);
566: }
568: /*@C
569: PetscViewerExodusIIGetNodalVariableNames - Gets the names of all nodal variables.
571: Collective; No Fortran Support
573: Input Parameters:
574: + viewer - a `PetscViewer` of type `PETSCVIEWEREXODUSII`
575: - numVars - the number of nodal variable names to retrieve
577: Output Parameter:
578: . varNames - returns an array of char pointers where the nodal variable names are
580: Level: intermediate
582: Notes:
583: This function returns a borrowed pointer which should not be freed.
585: .seealso: `PETSCVIEWEREXODUSII`, `PetscViewer`, `PetscViewerCreate()`, `PetscViewerDestroy()`, `PetscViewerExodusIIOpen()`, `PetscViewerSetType()`, `PetscViewerType`, `PetscViewerExodusIISetNodalVariableNames()`
586: @*/
587: PetscErrorCode PetscViewerExodusIIGetNodalVariableNames(PetscViewer viewer, PetscExodusIIInt *numVars, const char *const *varNames[])
588: {
589: PetscViewer_ExodusII *exo = (PetscViewer_ExodusII *)viewer->data;
590: PetscExodusIIInt idx;
591: char tmpName[MAX_NAME_LENGTH + 1];
592: PetscExodusIIInt exoid = -1;
594: PetscFunctionBegin;
595: PetscCall(PetscViewerExodusIIGetNodalVariable(viewer, numVars));
596: /*
597: Cache variable names if necessary
598: */
599: for (idx = 0; idx < *numVars; idx++) {
600: if (!exo->nodalVariableNames[idx]) {
601: PetscCall(PetscViewerExodusIIGetId(viewer, &exoid));
602: PetscCallExternal(ex_get_variable_name, exoid, EX_NODAL, idx + 1, tmpName);
603: PetscCall(PetscStrallocpy(tmpName, (char **)&exo->nodalVariableNames[idx]));
604: }
605: }
606: *varNames = (const char *const *)exo->nodalVariableNames;
607: PetscFunctionReturn(PETSC_SUCCESS);
608: }
610: /*MC
611: PETSCVIEWEREXODUSII - A viewer that writes to an Exodus II file
613: Level: beginner
615: .seealso: `PetscViewerExodusIIOpen()`, `PetscViewerCreate()`, `PETSCVIEWERBINARY`, `PETSCVIEWERHDF5`, `DMView()`,
616: `PetscViewerFileSetName()`, `PetscViewerFileSetMode()`, `PetscViewerFormat`, `PetscViewerType`, `PetscViewerSetType()`
617: M*/
618: PETSC_EXTERN PetscErrorCode PetscViewerCreate_ExodusII(PetscViewer v)
619: {
620: PetscViewer_ExodusII *exo;
622: PetscFunctionBegin;
623: PetscCall(PetscNew(&exo));
625: v->data = (void *)exo;
626: v->ops->destroy = PetscViewerDestroy_ExodusII;
627: v->ops->setfromoptions = PetscViewerSetFromOptions_ExodusII;
628: v->ops->view = PetscViewerView_ExodusII;
629: v->ops->flush = PetscViewerFlush_ExodusII;
630: exo->btype = FILE_MODE_UNDEFINED;
631: exo->filename = 0;
632: exo->exoid = -1;
633: exo->numNodalVariables = -1;
634: exo->numZonalVariables = -1;
635: exo->nodalVariableNames = NULL;
636: exo->zonalVariableNames = NULL;
638: PetscCall(PetscObjectComposeFunction((PetscObject)v, "PetscViewerFileSetName_C", PetscViewerFileSetName_ExodusII));
639: PetscCall(PetscObjectComposeFunction((PetscObject)v, "PetscViewerFileGetName_C", PetscViewerFileGetName_ExodusII));
640: PetscCall(PetscObjectComposeFunction((PetscObject)v, "PetscViewerFileSetMode_C", PetscViewerFileSetMode_ExodusII));
641: PetscCall(PetscObjectComposeFunction((PetscObject)v, "PetscViewerFileGetMode_C", PetscViewerFileGetMode_ExodusII));
642: PetscCall(PetscObjectComposeFunction((PetscObject)v, "PetscViewerGetId_C", PetscViewerExodusIIGetId_ExodusII));
643: PetscCall(PetscObjectComposeFunction((PetscObject)v, "PetscViewerSetOrder_C", PetscViewerExodusIISetOrder_ExodusII));
644: PetscCall(PetscObjectComposeFunction((PetscObject)v, "PetscViewerGetOrder_C", PetscViewerExodusIIGetOrder_ExodusII));
645: PetscFunctionReturn(PETSC_SUCCESS);
646: }
648: /*@
649: PetscViewerExodusIIGetNodalVariableIndex - return the location of a nodal variable in an ExodusII file given its name
651: Collective
653: Input Parameters:
654: + viewer - a `PetscViewer` of type `PETSCVIEWEREXODUSII`
655: - name - the name of the result
657: Output Parameter:
658: . varIndex - the location of the variable in the exodus file or -1 if the variable is not found
660: Level: beginner
662: Notes:
663: The exodus variable index is obtained by comparing the name argument to the
664: names of zonal variables declared in the exodus file. For instance if name is "V"
665: the location in the exodus file will be the first match of "V", "V_X", "V_XX", "V_1", or "V_11"
666: amongst all variables of type obj_type.
668: .seealso: `PetscViewerExodusIISetNodalVariable()`, `PetscViewerExodusIIGetNodalVariable()`, `PetscViewerExodusIISetNodalVariableName()`, `PetscViewerExodusIISetNodalVariableNames()`, `PetscViewerExodusIIGetNodalVariableName()`, `PetscViewerExodusIIGetNodalVariableNames()`
669: @*/
670: PetscErrorCode PetscViewerExodusIIGetNodalVariableIndex(PetscViewer viewer, const char name[], PetscExodusIIInt *varIndex)
671: {
672: PetscExodusIIInt num_vars = 0, i, j;
673: char ext_name[MAX_STR_LENGTH + 1];
674: const char *const *var_names;
675: const int num_suffix = 5;
676: char *suffix[5];
677: PetscBool flg;
679: PetscFunctionBegin;
680: suffix[0] = (char *)"";
681: suffix[1] = (char *)"_X";
682: suffix[2] = (char *)"_XX";
683: suffix[3] = (char *)"_1";
684: suffix[4] = (char *)"_11";
685: *varIndex = -1;
687: PetscCall(PetscViewerExodusIIGetNodalVariableNames(viewer, &num_vars, &var_names));
688: for (i = 0; i < num_vars; ++i) {
689: for (j = 0; j < num_suffix; ++j) {
690: PetscCall(PetscStrncpy(ext_name, name, MAX_STR_LENGTH));
691: PetscCall(PetscStrlcat(ext_name, suffix[j], MAX_STR_LENGTH));
692: PetscCall(PetscStrcasecmp(ext_name, var_names[i], &flg));
693: if (flg) *varIndex = i;
694: }
695: if (flg) break;
696: }
697: PetscFunctionReturn(PETSC_SUCCESS);
698: }
700: /*@
701: PetscViewerExodusIIGetZonalVariableIndex - return the location of a zonal variable in an ExodusII file given its name
703: Collective
705: Input Parameters:
706: + viewer - a `PetscViewer` of type `PETSCVIEWEREXODUSII`
707: - name - the name of the result
709: Output Parameter:
710: . varIndex - the location of the variable in the exodus file or -1 if the variable is not found
712: Level: beginner
714: Notes:
715: The exodus variable index is obtained by comparing the name argument to the
716: names of zonal variables declared in the exodus file. For instance if name is "V"
717: the location in the exodus file will be the first match of "V", "V_X", "V_XX", "V_1", or "V_11"
718: amongst all variables of type obj_type.
720: .seealso: `PetscViewerExodusIISetNodalVariable()`, `PetscViewerExodusIIGetNodalVariable()`, `PetscViewerExodusIISetNodalVariableName()`, `PetscViewerExodusIISetNodalVariableNames()`, `PetscViewerExodusIIGetNodalVariableName()`, `PetscViewerExodusIIGetNodalVariableNames()`
721: @*/
722: PetscErrorCode PetscViewerExodusIIGetZonalVariableIndex(PetscViewer viewer, const char name[], int *varIndex)
723: {
724: PetscExodusIIInt num_vars = 0, i, j;
725: char ext_name[MAX_STR_LENGTH + 1];
726: const char *const *var_names;
727: const int num_suffix = 5;
728: char *suffix[5];
729: PetscBool flg;
731: PetscFunctionBegin;
732: suffix[0] = (char *)"";
733: suffix[1] = (char *)"_X";
734: suffix[2] = (char *)"_XX";
735: suffix[3] = (char *)"_1";
736: suffix[4] = (char *)"_11";
737: *varIndex = -1;
739: PetscCall(PetscViewerExodusIIGetZonalVariableNames(viewer, &num_vars, &var_names));
740: for (i = 0; i < num_vars; ++i) {
741: for (j = 0; j < num_suffix; ++j) {
742: PetscCall(PetscStrncpy(ext_name, name, MAX_STR_LENGTH));
743: PetscCall(PetscStrlcat(ext_name, suffix[j], MAX_STR_LENGTH));
744: PetscCall(PetscStrcasecmp(ext_name, var_names[i], &flg));
745: if (flg) *varIndex = i;
746: }
747: if (flg) break;
748: }
749: PetscFunctionReturn(PETSC_SUCCESS);
750: }
752: PetscErrorCode DMView_PlexExodusII(DM dm, PetscViewer viewer)
753: {
754: enum ElemType {
755: SEGMENT,
756: TRI,
757: QUAD,
758: TET,
759: HEX
760: };
761: MPI_Comm comm;
762: PetscInt degree; /* the order of the mesh */
763: /* Connectivity Variables */
764: PetscInt cellsNotInConnectivity;
765: /* Cell Sets */
766: DMLabel csLabel;
767: IS csIS;
768: const PetscInt *csIdx;
769: PetscInt num_cs, cs;
770: enum ElemType *type;
771: PetscBool hasLabel;
772: /* Coordinate Variables */
773: DM cdm;
774: PetscSection coordSection;
775: Vec coord;
776: PetscInt **nodes;
777: PetscInt depth, d, dim, skipCells = 0;
778: PetscInt pStart, pEnd, p, cStart, cEnd, numCells, vStart, vEnd, numVertices, eStart, eEnd, numEdges, fStart, fEnd, numFaces, numNodes;
779: PetscInt num_vs, num_fs;
780: PetscMPIInt rank, size;
781: const char *dmName;
782: PetscInt nodesLineP1[4] = {2, 0, 0, 0};
783: PetscInt nodesLineP2[4] = {2, 0, 0, 1};
784: PetscInt nodesTriP1[4] = {3, 0, 0, 0};
785: PetscInt nodesTriP2[4] = {3, 3, 0, 0};
786: PetscInt nodesQuadP1[4] = {4, 0, 0, 0};
787: PetscInt nodesQuadP2[4] = {4, 4, 0, 1};
788: PetscInt nodesTetP1[4] = {4, 0, 0, 0};
789: PetscInt nodesTetP2[4] = {4, 6, 0, 0};
790: PetscInt nodesHexP1[4] = {8, 0, 0, 0};
791: PetscInt nodesHexP2[4] = {8, 12, 6, 1};
792: PetscExodusIIInt CPU_word_size, IO_word_size, EXO_mode;
793: PetscExodusIIFloat EXO_version;
795: PetscViewer_ExodusII *exo = (PetscViewer_ExodusII *)viewer->data;
797: PetscFunctionBegin;
798: PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
799: PetscCallMPI(MPI_Comm_rank(comm, &rank));
800: PetscCallMPI(MPI_Comm_size(comm, &size));
802: /*
803: Creating coordSection is a collective operation so we do it somewhat out of sequence
804: */
805: PetscCall(PetscSectionCreate(comm, &coordSection));
806: PetscCall(DMGetCoordinatesLocalSetUp(dm));
807: /*
808: Check that all points are on rank 0 since we don't know how to save distributed DM in exodus format
809: */
810: PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
811: PetscCall(DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd));
812: PetscCall(DMPlexGetDepthStratum(dm, 1, &eStart, &eEnd));
813: PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd));
814: numCells = cEnd - cStart;
815: numEdges = eEnd - eStart;
816: numVertices = vEnd - vStart;
817: PetscCheck(!(rank && (numCells || numEdges || numVertices)), PETSC_COMM_SELF, PETSC_ERR_SUP, "Writing distributed DM in ExodusII format not supported");
818: if (rank == 0) {
819: switch (exo->btype) {
820: case FILE_MODE_READ:
821: case FILE_MODE_APPEND:
822: case FILE_MODE_UPDATE:
823: case FILE_MODE_APPEND_UPDATE:
824: /* ExodusII does not allow writing geometry to an existing file */
825: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "cannot add geometry to existing file %s", exo->filename);
826: case FILE_MODE_WRITE:
827: /* Create an empty file if one already exists*/
828: EXO_mode = EX_CLOBBER;
829: #if defined(PETSC_USE_64BIT_INDICES)
830: EXO_mode += EX_ALL_INT64_API;
831: #endif
832: CPU_word_size = sizeof(PetscReal);
833: IO_word_size = sizeof(PetscReal);
834: exo->exoid = ex_create(exo->filename, EXO_mode, &CPU_word_size, &IO_word_size);
835: PetscCheck(exo->exoid >= 0, PETSC_COMM_SELF, PETSC_ERR_LIB, "ex_create failed for %s", exo->filename);
837: break;
838: default:
839: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ORDER, "Must call PetscViewerFileSetMode() before PetscViewerFileSetName()");
840: }
842: /* --- Get DM info --- */
843: PetscCall(PetscObjectGetName((PetscObject)dm, &dmName));
844: PetscCall(DMPlexGetDepth(dm, &depth));
845: PetscCall(DMGetDimension(dm, &dim));
846: PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
847: if (depth == 3) {
848: numFaces = fEnd - fStart;
849: } else {
850: numFaces = 0;
851: }
852: PetscCall(DMGetLabelSize(dm, "Cell Sets", &num_cs));
853: PetscCall(DMGetLabelSize(dm, "Vertex Sets", &num_vs));
854: PetscCall(DMGetLabelSize(dm, "Face Sets", &num_fs));
855: PetscCall(DMGetCoordinatesLocal(dm, &coord));
856: PetscCall(DMGetCoordinateDM(dm, &cdm));
857: if (num_cs > 0) {
858: PetscCall(DMGetLabel(dm, "Cell Sets", &csLabel));
859: PetscCall(DMLabelGetValueIS(csLabel, &csIS));
860: PetscCall(ISGetIndices(csIS, &csIdx));
861: }
862: PetscCall(PetscMalloc1(num_cs, &nodes));
863: /* Set element type for each block and compute total number of nodes */
864: PetscCall(PetscMalloc1(num_cs, &type));
865: numNodes = numVertices;
867: PetscCall(PetscViewerExodusIIGetOrder(viewer, °ree));
868: if (degree == 2) numNodes += numEdges;
869: cellsNotInConnectivity = numCells;
870: for (cs = 0; cs < num_cs; ++cs) {
871: IS stratumIS;
872: const PetscInt *cells;
873: PetscScalar *xyz = NULL;
874: PetscInt csSize, closureSize;
876: PetscCall(DMLabelGetStratumIS(csLabel, csIdx[cs], &stratumIS));
877: PetscCall(ISGetIndices(stratumIS, &cells));
878: PetscCall(ISGetSize(stratumIS, &csSize));
879: PetscCall(DMPlexVecGetClosure(cdm, NULL, coord, cells[0], &closureSize, &xyz));
880: switch (dim) {
881: case 1:
882: 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);
883: type[cs] = SEGMENT;
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, MPI_C_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: }