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