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: 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);
884: type[cs] = SEGMENT;
885: break;
886: case 2:
887: if (closureSize == 3 * dim) {
888: type[cs] = TRI;
889: } else if (closureSize == 4 * dim) {
890: type[cs] = QUAD;
891: } 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);
892: break;
893: case 3:
894: if (closureSize == 4 * dim) {
895: type[cs] = TET;
896: } else if (closureSize == 8 * dim) {
897: type[cs] = HEX;
898: } 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);
899: break;
900: default:
901: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Dimension %" PetscInt_FMT " not handled by ExodusII viewer", dim);
902: }
903: if ((degree == 2) && (type[cs] == SEGMENT)) numNodes += csSize;
904: if ((degree == 2) && (type[cs] == QUAD)) numNodes += csSize;
905: if ((degree == 2) && (type[cs] == HEX)) {
906: numNodes += csSize;
907: numNodes += numFaces;
908: }
909: PetscCall(DMPlexVecRestoreClosure(cdm, NULL, coord, cells[0], &closureSize, &xyz));
910: /* Set nodes and Element type */
911: if (type[cs] == SEGMENT) {
912: if (degree == 1) nodes[cs] = nodesLineP1;
913: else if (degree == 2) nodes[cs] = nodesLineP2;
914: } else if (type[cs] == TRI) {
915: if (degree == 1) nodes[cs] = nodesTriP1;
916: else if (degree == 2) nodes[cs] = nodesTriP2;
917: } else if (type[cs] == QUAD) {
918: if (degree == 1) nodes[cs] = nodesQuadP1;
919: else if (degree == 2) nodes[cs] = nodesQuadP2;
920: } else if (type[cs] == TET) {
921: if (degree == 1) nodes[cs] = nodesTetP1;
922: else if (degree == 2) nodes[cs] = nodesTetP2;
923: } else if (type[cs] == HEX) {
924: if (degree == 1) nodes[cs] = nodesHexP1;
925: else if (degree == 2) nodes[cs] = nodesHexP2;
926: }
927: /* Compute the number of cells not in the connectivity table */
928: cellsNotInConnectivity -= nodes[cs][3] * csSize;
930: PetscCall(ISRestoreIndices(stratumIS, &cells));
931: PetscCall(ISDestroy(&stratumIS));
932: }
933: if (num_cs) PetscCallExternal(ex_put_init, exo->exoid, dmName, dim, numNodes, numCells, num_cs, num_vs, num_fs);
934: /* --- Connectivity --- */
935: for (cs = 0; cs < num_cs; ++cs) {
936: IS stratumIS;
937: const PetscInt *cells;
938: PetscInt *connect, off = 0;
939: PetscInt edgesInClosure = 0, facesInClosure = 0, verticesInClosure = 0;
940: PetscInt csSize, c, connectSize, closureSize;
941: char *elem_type = NULL;
942: char elem_type_bar2[] = "BAR2", elem_type_bar3[] = "BAR3";
943: char elem_type_tri3[] = "TRI3", elem_type_quad4[] = "QUAD4";
944: char elem_type_tri6[] = "TRI6", elem_type_quad9[] = "QUAD9";
945: char elem_type_tet4[] = "TET4", elem_type_hex8[] = "HEX8";
946: char elem_type_tet10[] = "TET10", elem_type_hex27[] = "HEX27";
948: PetscCall(DMLabelGetStratumIS(csLabel, csIdx[cs], &stratumIS));
949: PetscCall(ISGetIndices(stratumIS, &cells));
950: PetscCall(ISGetSize(stratumIS, &csSize));
951: /* Set Element type */
952: if (type[cs] == SEGMENT) {
953: if (degree == 1) elem_type = elem_type_bar2;
954: else if (degree == 2) elem_type = elem_type_bar3;
955: } else if (type[cs] == TRI) {
956: if (degree == 1) elem_type = elem_type_tri3;
957: else if (degree == 2) elem_type = elem_type_tri6;
958: } else if (type[cs] == QUAD) {
959: if (degree == 1) elem_type = elem_type_quad4;
960: else if (degree == 2) elem_type = elem_type_quad9;
961: } else if (type[cs] == TET) {
962: if (degree == 1) elem_type = elem_type_tet4;
963: else if (degree == 2) elem_type = elem_type_tet10;
964: } else if (type[cs] == HEX) {
965: if (degree == 1) elem_type = elem_type_hex8;
966: else if (degree == 2) elem_type = elem_type_hex27;
967: }
968: connectSize = nodes[cs][0] + nodes[cs][1] + nodes[cs][2] + nodes[cs][3];
969: PetscCall(PetscMalloc1(PetscMax(27, connectSize) * csSize, &connect));
970: PetscCallExternal(ex_put_block, exo->exoid, EX_ELEM_BLOCK, csIdx[cs], elem_type, csSize, connectSize, 0, 0, 1);
971: /* Find number of vertices, edges, and faces in the closure */
972: verticesInClosure = nodes[cs][0];
973: if (depth > 1) {
974: if (dim == 2) {
975: PetscCall(DMPlexGetConeSize(dm, cells[0], &edgesInClosure));
976: } else if (dim == 3) {
977: PetscInt *closure = NULL;
979: PetscCall(DMPlexGetConeSize(dm, cells[0], &facesInClosure));
980: PetscCall(DMPlexGetTransitiveClosure(dm, cells[0], PETSC_TRUE, &closureSize, &closure));
981: edgesInClosure = closureSize - facesInClosure - 1 - verticesInClosure;
982: PetscCall(DMPlexRestoreTransitiveClosure(dm, cells[0], PETSC_TRUE, &closureSize, &closure));
983: }
984: }
985: /* Get connectivity for each cell */
986: for (c = 0; c < csSize; ++c) {
987: PetscInt *closure = NULL;
988: PetscInt temp, i;
990: PetscCall(DMPlexGetTransitiveClosure(dm, cells[c], PETSC_TRUE, &closureSize, &closure));
991: for (i = 0; i < connectSize; ++i) {
992: if (i < nodes[cs][0]) { /* Vertices */
993: connect[i + off] = closure[(i + edgesInClosure + facesInClosure + 1) * 2] + 1;
994: connect[i + off] -= cellsNotInConnectivity;
995: } else if (i < nodes[cs][0] + nodes[cs][1]) { /* Edges */
996: connect[i + off] = closure[(i - verticesInClosure + facesInClosure + 1) * 2] + 1;
997: if (nodes[cs][2] == 0) connect[i + off] -= numFaces;
998: connect[i + off] -= cellsNotInConnectivity;
999: } else if (i < nodes[cs][0] + nodes[cs][1] + nodes[cs][3]) { /* Cells */
1000: connect[i + off] = closure[0] + 1;
1001: connect[i + off] -= skipCells;
1002: } else if (i < nodes[cs][0] + nodes[cs][1] + nodes[cs][3] + nodes[cs][2]) { /* Faces */
1003: connect[i + off] = closure[(i - edgesInClosure - verticesInClosure) * 2] + 1;
1004: connect[i + off] -= cellsNotInConnectivity;
1005: } else {
1006: connect[i + off] = -1;
1007: }
1008: }
1009: /* Tetrahedra are inverted */
1010: if (type[cs] == TET) {
1011: temp = connect[0 + off];
1012: connect[0 + off] = connect[1 + off];
1013: connect[1 + off] = temp;
1014: if (degree == 2) {
1015: temp = connect[5 + off];
1016: connect[5 + off] = connect[6 + off];
1017: connect[6 + off] = temp;
1018: temp = connect[7 + off];
1019: connect[7 + off] = connect[8 + off];
1020: connect[8 + off] = temp;
1021: }
1022: }
1023: /* Hexahedra are inverted */
1024: if (type[cs] == HEX) {
1025: temp = connect[1 + off];
1026: connect[1 + off] = connect[3 + off];
1027: connect[3 + off] = temp;
1028: if (degree == 2) {
1029: temp = connect[8 + off];
1030: connect[8 + off] = connect[11 + off];
1031: connect[11 + off] = temp;
1032: temp = connect[9 + off];
1033: connect[9 + off] = connect[10 + off];
1034: connect[10 + off] = temp;
1035: temp = connect[16 + off];
1036: connect[16 + off] = connect[17 + off];
1037: connect[17 + off] = temp;
1038: temp = connect[18 + off];
1039: connect[18 + off] = connect[19 + off];
1040: connect[19 + off] = temp;
1042: temp = connect[12 + off];
1043: connect[12 + off] = connect[16 + off];
1044: connect[16 + off] = temp;
1045: temp = connect[13 + off];
1046: connect[13 + off] = connect[17 + off];
1047: connect[17 + off] = temp;
1048: temp = connect[14 + off];
1049: connect[14 + off] = connect[18 + off];
1050: connect[18 + off] = temp;
1051: temp = connect[15 + off];
1052: connect[15 + off] = connect[19 + off];
1053: connect[19 + off] = temp;
1055: temp = connect[23 + off];
1056: connect[23 + off] = connect[26 + off];
1057: connect[26 + off] = temp;
1058: temp = connect[24 + off];
1059: connect[24 + off] = connect[25 + off];
1060: connect[25 + off] = temp;
1061: temp = connect[25 + off];
1062: connect[25 + off] = connect[26 + off];
1063: connect[26 + off] = temp;
1064: }
1065: }
1066: off += connectSize;
1067: PetscCall(DMPlexRestoreTransitiveClosure(dm, cells[c], PETSC_TRUE, &closureSize, &closure));
1068: }
1069: PetscCallExternal(ex_put_conn, exo->exoid, EX_ELEM_BLOCK, csIdx[cs], connect, 0, 0);
1070: skipCells += (nodes[cs][3] == 0) * csSize;
1071: PetscCall(PetscFree(connect));
1072: PetscCall(ISRestoreIndices(stratumIS, &cells));
1073: PetscCall(ISDestroy(&stratumIS));
1074: }
1075: PetscCall(PetscFree(type));
1076: /* --- Coordinates --- */
1077: PetscCall(PetscSectionSetChart(coordSection, pStart, pEnd));
1078: if (num_cs) {
1079: for (d = 0; d < depth; ++d) {
1080: PetscCall(DMPlexGetDepthStratum(dm, d, &pStart, &pEnd));
1081: for (p = pStart; p < pEnd; ++p) PetscCall(PetscSectionSetDof(coordSection, p, nodes[0][d] > 0));
1082: }
1083: }
1084: for (cs = 0; cs < num_cs; ++cs) {
1085: IS stratumIS;
1086: const PetscInt *cells;
1087: PetscInt csSize, c;
1089: PetscCall(DMLabelGetStratumIS(csLabel, csIdx[cs], &stratumIS));
1090: PetscCall(ISGetIndices(stratumIS, &cells));
1091: PetscCall(ISGetSize(stratumIS, &csSize));
1092: for (c = 0; c < csSize; ++c) PetscCall(PetscSectionSetDof(coordSection, cells[c], nodes[cs][3] > 0));
1093: PetscCall(ISRestoreIndices(stratumIS, &cells));
1094: PetscCall(ISDestroy(&stratumIS));
1095: }
1096: if (num_cs) {
1097: PetscCall(ISRestoreIndices(csIS, &csIdx));
1098: PetscCall(ISDestroy(&csIS));
1099: }
1100: PetscCall(PetscFree(nodes));
1101: PetscCall(PetscSectionSetUp(coordSection));
1102: if (numNodes) {
1103: const char *coordNames[3] = {"x", "y", "z"};
1104: PetscScalar *closure, *cval;
1105: PetscReal *coords;
1106: PetscInt hasDof, n = 0;
1108: /* There can't be more than 24 values in the closure of a point for the coord coordSection */
1109: PetscCall(PetscCalloc3(numNodes * 3, &coords, dim, &cval, 24, &closure));
1110: PetscCall(DMGetCoordinatesLocalNoncollective(dm, &coord));
1111: PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
1112: for (p = pStart; p < pEnd; ++p) {
1113: PetscCall(PetscSectionGetDof(coordSection, p, &hasDof));
1114: if (hasDof) {
1115: PetscInt closureSize = 24, j;
1117: PetscCall(DMPlexVecGetClosure(cdm, NULL, coord, p, &closureSize, &closure));
1118: for (d = 0; d < dim; ++d) {
1119: cval[d] = 0.0;
1120: for (j = 0; j < closureSize / dim; j++) cval[d] += closure[j * dim + d];
1121: coords[d * numNodes + n] = PetscRealPart(cval[d]) * dim / closureSize;
1122: }
1123: ++n;
1124: }
1125: }
1126: PetscCallExternal(ex_put_coord, exo->exoid, &coords[0 * numNodes], &coords[1 * numNodes], &coords[2 * numNodes]);
1127: PetscCall(PetscFree3(coords, cval, closure));
1128: PetscCallExternal(ex_put_coord_names, exo->exoid, (char **)coordNames);
1129: }
1131: /* --- Node Sets/Vertex Sets --- */
1132: PetscCall(DMHasLabel(dm, "Vertex Sets", &hasLabel));
1133: if (hasLabel) {
1134: PetscInt i, vs, vsSize;
1135: const PetscInt *vsIdx, *vertices;
1136: PetscInt *nodeList;
1137: IS vsIS, stratumIS;
1138: DMLabel vsLabel;
1139: PetscCall(DMGetLabel(dm, "Vertex Sets", &vsLabel));
1140: PetscCall(DMLabelGetValueIS(vsLabel, &vsIS));
1141: PetscCall(ISGetIndices(vsIS, &vsIdx));
1142: for (vs = 0; vs < num_vs; ++vs) {
1143: PetscCall(DMLabelGetStratumIS(vsLabel, vsIdx[vs], &stratumIS));
1144: PetscCall(ISGetIndices(stratumIS, &vertices));
1145: PetscCall(ISGetSize(stratumIS, &vsSize));
1146: PetscCall(PetscMalloc1(vsSize, &nodeList));
1147: for (i = 0; i < vsSize; ++i) nodeList[i] = vertices[i] - skipCells + 1;
1148: PetscCallExternal(ex_put_set_param, exo->exoid, EX_NODE_SET, vsIdx[vs], vsSize, 0);
1149: PetscCallExternal(ex_put_set, exo->exoid, EX_NODE_SET, vsIdx[vs], nodeList, NULL);
1150: PetscCall(ISRestoreIndices(stratumIS, &vertices));
1151: PetscCall(ISDestroy(&stratumIS));
1152: PetscCall(PetscFree(nodeList));
1153: }
1154: PetscCall(ISRestoreIndices(vsIS, &vsIdx));
1155: PetscCall(ISDestroy(&vsIS));
1156: }
1157: /* --- Side Sets/Face Sets --- */
1158: PetscCall(DMHasLabel(dm, "Face Sets", &hasLabel));
1159: if (hasLabel) {
1160: PetscInt i, j, fs, fsSize;
1161: const PetscInt *fsIdx, *faces;
1162: IS fsIS, stratumIS;
1163: DMLabel fsLabel;
1164: PetscInt numPoints, *points;
1165: PetscInt elem_list_size = 0;
1166: PetscInt *elem_list, *elem_ind, *side_list;
1168: PetscCall(DMGetLabel(dm, "Face Sets", &fsLabel));
1169: /* Compute size of Node List and Element List */
1170: PetscCall(DMLabelGetValueIS(fsLabel, &fsIS));
1171: PetscCall(ISGetIndices(fsIS, &fsIdx));
1172: for (fs = 0; fs < num_fs; ++fs) {
1173: PetscCall(DMLabelGetStratumIS(fsLabel, fsIdx[fs], &stratumIS));
1174: PetscCall(ISGetSize(stratumIS, &fsSize));
1175: elem_list_size += fsSize;
1176: PetscCall(ISDestroy(&stratumIS));
1177: }
1178: if (num_fs) {
1179: PetscCall(PetscMalloc3(num_fs, &elem_ind, elem_list_size, &elem_list, elem_list_size, &side_list));
1180: elem_ind[0] = 0;
1181: for (fs = 0; fs < num_fs; ++fs) {
1182: PetscCall(DMLabelGetStratumIS(fsLabel, fsIdx[fs], &stratumIS));
1183: PetscCall(ISGetIndices(stratumIS, &faces));
1184: PetscCall(ISGetSize(stratumIS, &fsSize));
1185: /* Set Parameters */
1186: PetscCallExternal(ex_put_set_param, exo->exoid, EX_SIDE_SET, fsIdx[fs], fsSize, 0);
1187: /* Indices */
1188: if (fs < num_fs - 1) elem_ind[fs + 1] = elem_ind[fs] + fsSize;
1190: for (i = 0; i < fsSize; ++i) {
1191: /* Element List */
1192: points = NULL;
1193: PetscCall(DMPlexGetTransitiveClosure(dm, faces[i], PETSC_FALSE, &numPoints, &points));
1194: elem_list[elem_ind[fs] + i] = points[2] + 1;
1195: PetscCall(DMPlexRestoreTransitiveClosure(dm, faces[i], PETSC_FALSE, &numPoints, &points));
1197: /* Side List */
1198: points = NULL;
1199: PetscCall(DMPlexGetTransitiveClosure(dm, elem_list[elem_ind[fs] + i] - 1, PETSC_TRUE, &numPoints, &points));
1200: for (j = 1; j < numPoints; ++j) {
1201: if (points[j * 2] == faces[i]) break;
1202: }
1203: /* Convert HEX sides */
1204: if (numPoints == 27) {
1205: if (j == 1) {
1206: j = 5;
1207: } else if (j == 2) {
1208: j = 6;
1209: } else if (j == 3) {
1210: j = 1;
1211: } else if (j == 4) {
1212: j = 3;
1213: } else if (j == 5) {
1214: j = 2;
1215: } else if (j == 6) {
1216: j = 4;
1217: }
1218: }
1219: /* Convert TET sides */
1220: if (numPoints == 15) {
1221: --j;
1222: if (j == 0) j = 4;
1223: }
1224: side_list[elem_ind[fs] + i] = j;
1225: PetscCall(DMPlexRestoreTransitiveClosure(dm, elem_list[elem_ind[fs] + i] - 1, PETSC_TRUE, &numPoints, &points));
1226: }
1227: PetscCall(ISRestoreIndices(stratumIS, &faces));
1228: PetscCall(ISDestroy(&stratumIS));
1229: }
1230: PetscCall(ISRestoreIndices(fsIS, &fsIdx));
1231: PetscCall(ISDestroy(&fsIS));
1233: /* Put side sets */
1234: 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]]);
1235: PetscCall(PetscFree3(elem_ind, elem_list, side_list));
1236: }
1237: }
1238: /*
1239: close the exodus file
1240: */
1241: ex_close(exo->exoid);
1242: exo->exoid = -1;
1243: }
1244: PetscCall(PetscSectionDestroy(&coordSection));
1246: /*
1247: reopen the file in parallel
1248: */
1249: EXO_mode = EX_WRITE;
1250: #if defined(PETSC_USE_64BIT_INDICES)
1251: EXO_mode += EX_ALL_INT64_API;
1252: #endif
1253: CPU_word_size = sizeof(PetscReal);
1254: IO_word_size = sizeof(PetscReal);
1255: exo->exoid = ex_open_par(exo->filename, EXO_mode, &CPU_word_size, &IO_word_size, &EXO_version, comm, MPI_INFO_NULL);
1256: PetscCheck(exo->exoid >= 0, PETSC_COMM_SELF, PETSC_ERR_LIB, "ex_open_par failed for %s", exo->filename);
1257: PetscFunctionReturn(PETSC_SUCCESS);
1258: }
1260: static PetscErrorCode VecViewPlex_ExodusII_Nodal_Internal(Vec v, PetscExodusIIInt exoid, PetscExodusIIInt step, PetscExodusIIInt offset);
1261: static PetscErrorCode VecViewPlex_ExodusII_Zonal_Internal(Vec v, PetscExodusIIInt exoid, PetscExodusIIInt step, PetscExodusIIInt offset);
1262: static PetscErrorCode VecLoadPlex_ExodusII_Nodal_Internal(Vec v, PetscExodusIIInt exoid, PetscExodusIIInt step, PetscExodusIIInt offset);
1263: static PetscErrorCode VecLoadPlex_ExodusII_Zonal_Internal(Vec v, PetscExodusIIInt exoid, PetscExodusIIInt step, PetscExodusIIInt offset);
1265: PetscErrorCode VecView_PlexExodusII_Internal(Vec v, PetscViewer viewer)
1266: {
1267: DM dm;
1268: MPI_Comm comm;
1269: PetscMPIInt rank;
1271: PetscExodusIIInt exoid, offsetN = -1, offsetZ = -1;
1272: const char *vecname;
1273: PetscInt step;
1275: PetscFunctionBegin;
1276: PetscCall(PetscObjectGetComm((PetscObject)v, &comm));
1277: PetscCallMPI(MPI_Comm_rank(comm, &rank));
1278: PetscCall(PetscViewerExodusIIGetId(viewer, &exoid));
1279: PetscCall(VecGetDM(v, &dm));
1280: PetscCall(PetscObjectGetName((PetscObject)v, &vecname));
1282: PetscCall(DMGetOutputSequenceNumber(dm, &step, NULL));
1283: PetscCall(PetscViewerExodusIIGetNodalVariableIndex(viewer, vecname, &offsetN));
1284: PetscCall(PetscViewerExodusIIGetZonalVariableIndex(viewer, vecname, &offsetZ));
1285: PetscCheck(!(offsetN >= 0 && offsetZ >= 0), comm, PETSC_ERR_FILE_UNEXPECTED, "Found both nodal and zonal variable %s in exodus file. ", vecname);
1286: if (offsetN >= 0) {
1287: PetscCall(VecViewPlex_ExodusII_Nodal_Internal(v, exoid, (PetscExodusIIInt)step + 1, offsetN + 1));
1288: } else if (offsetZ >= 0) {
1289: PetscCall(VecViewPlex_ExodusII_Zonal_Internal(v, exoid, (PetscExodusIIInt)step + 1, offsetZ + 1));
1290: } else SETERRQ(comm, PETSC_ERR_FILE_UNEXPECTED, "Could not find nodal or zonal variable %s in exodus file. ", vecname);
1291: PetscFunctionReturn(PETSC_SUCCESS);
1292: }
1294: PetscErrorCode VecLoad_PlexExodusII_Internal(Vec v, PetscViewer viewer)
1295: {
1296: DM dm;
1297: MPI_Comm comm;
1298: PetscMPIInt rank;
1300: PetscExodusIIInt exoid, offsetN = 0, offsetZ = 0;
1301: const char *vecname;
1302: PetscInt step;
1304: PetscFunctionBegin;
1305: PetscCall(PetscObjectGetComm((PetscObject)v, &comm));
1306: PetscCallMPI(MPI_Comm_rank(comm, &rank));
1307: PetscCall(PetscViewerExodusIIGetId(viewer, &exoid));
1308: PetscCall(VecGetDM(v, &dm));
1309: PetscCall(PetscObjectGetName((PetscObject)v, &vecname));
1311: PetscCall(DMGetOutputSequenceNumber(dm, &step, NULL));
1312: PetscCall(PetscViewerExodusIIGetNodalVariableIndex(viewer, vecname, &offsetN));
1313: PetscCall(PetscViewerExodusIIGetZonalVariableIndex(viewer, vecname, &offsetZ));
1314: PetscCheck(!(offsetN >= 0 && offsetZ >= 0), comm, PETSC_ERR_FILE_UNEXPECTED, "Found both nodal and zonal variable %s in exodus file. ", vecname);
1315: if (offsetN >= 0) PetscCall(VecLoadPlex_ExodusII_Nodal_Internal(v, exoid, (PetscExodusIIInt)step + 1, offsetN + 1));
1316: else if (offsetZ >= 0) PetscCall(VecLoadPlex_ExodusII_Zonal_Internal(v, exoid, (PetscExodusIIInt)step + 1, offsetZ + 1));
1317: else SETERRQ(comm, PETSC_ERR_FILE_UNEXPECTED, "Could not find nodal or zonal variable %s in exodus file. ", vecname);
1318: PetscFunctionReturn(PETSC_SUCCESS);
1319: }
1321: static PetscErrorCode VecViewPlex_ExodusII_Nodal_Internal(Vec v, PetscExodusIIInt exoid, PetscExodusIIInt step, PetscExodusIIInt offset)
1322: {
1323: MPI_Comm comm;
1324: PetscMPIInt size;
1325: DM dm;
1326: Vec vNatural, vComp;
1327: const PetscScalar *varray;
1328: PetscInt xs, xe, bs;
1329: PetscBool useNatural;
1331: PetscFunctionBegin;
1332: PetscCall(PetscObjectGetComm((PetscObject)v, &comm));
1333: PetscCallMPI(MPI_Comm_size(comm, &size));
1334: PetscCall(VecGetDM(v, &dm));
1335: PetscCall(DMGetUseNatural(dm, &useNatural));
1336: useNatural = useNatural && size > 1 ? PETSC_TRUE : PETSC_FALSE;
1337: if (useNatural) {
1338: PetscCall(DMPlexCreateNaturalVector(dm, &vNatural));
1339: PetscCall(DMPlexGlobalToNaturalBegin(dm, v, vNatural));
1340: PetscCall(DMPlexGlobalToNaturalEnd(dm, v, vNatural));
1341: } else {
1342: vNatural = v;
1343: }
1345: /* Write local chunk of the result in the exodus file
1346: exodus stores each component of a vector-valued field as a separate variable.
1347: We assume that they are stored sequentially */
1348: PetscCall(VecGetOwnershipRange(vNatural, &xs, &xe));
1349: PetscCall(VecGetBlockSize(vNatural, &bs));
1350: if (bs == 1) {
1351: PetscCall(VecGetArrayRead(vNatural, &varray));
1352: PetscCallExternal(ex_put_partial_var, exoid, step, EX_NODAL, offset, 1, xs + 1, xe - xs, varray);
1353: PetscCall(VecRestoreArrayRead(vNatural, &varray));
1354: } else {
1355: IS compIS;
1356: PetscInt c;
1358: PetscCall(ISCreateStride(comm, (xe - xs) / bs, xs, bs, &compIS));
1359: for (c = 0; c < bs; ++c) {
1360: PetscCall(ISStrideSetStride(compIS, (xe - xs) / bs, xs + c, bs));
1361: PetscCall(VecGetSubVector(vNatural, compIS, &vComp));
1362: PetscCall(VecGetArrayRead(vComp, &varray));
1363: PetscCallExternal(ex_put_partial_var, exoid, step, EX_NODAL, offset + c, 1, xs / bs + 1, (xe - xs) / bs, varray);
1364: PetscCall(VecRestoreArrayRead(vComp, &varray));
1365: PetscCall(VecRestoreSubVector(vNatural, compIS, &vComp));
1366: }
1367: PetscCall(ISDestroy(&compIS));
1368: }
1369: if (useNatural) PetscCall(VecDestroy(&vNatural));
1370: PetscFunctionReturn(PETSC_SUCCESS);
1371: }
1373: static PetscErrorCode VecLoadPlex_ExodusII_Nodal_Internal(Vec v, PetscExodusIIInt exoid, PetscExodusIIInt step, PetscExodusIIInt offset)
1374: {
1375: MPI_Comm comm;
1376: PetscMPIInt size;
1377: DM dm;
1378: Vec vNatural, vComp;
1379: PetscScalar *varray;
1380: PetscInt xs, xe, bs;
1381: PetscBool useNatural;
1383: PetscFunctionBegin;
1384: PetscCall(PetscObjectGetComm((PetscObject)v, &comm));
1385: PetscCallMPI(MPI_Comm_size(comm, &size));
1386: PetscCall(VecGetDM(v, &dm));
1387: PetscCall(DMGetUseNatural(dm, &useNatural));
1388: useNatural = useNatural && size > 1 ? PETSC_TRUE : PETSC_FALSE;
1389: if (useNatural) PetscCall(DMPlexCreateNaturalVector(dm, &vNatural));
1390: else vNatural = v;
1392: /* Read local chunk from the file */
1393: PetscCall(VecGetOwnershipRange(vNatural, &xs, &xe));
1394: PetscCall(VecGetBlockSize(vNatural, &bs));
1395: if (bs == 1) {
1396: PetscCall(VecGetArray(vNatural, &varray));
1397: PetscCallExternal(ex_get_partial_var, exoid, step, EX_NODAL, offset, 1, xs + 1, xe - xs, varray);
1398: PetscCall(VecRestoreArray(vNatural, &varray));
1399: } else {
1400: IS compIS;
1401: PetscInt c;
1403: PetscCall(ISCreateStride(comm, (xe - xs) / bs, xs, bs, &compIS));
1404: for (c = 0; c < bs; ++c) {
1405: PetscCall(ISStrideSetStride(compIS, (xe - xs) / bs, xs + c, bs));
1406: PetscCall(VecGetSubVector(vNatural, compIS, &vComp));
1407: PetscCall(VecGetArray(vComp, &varray));
1408: PetscCallExternal(ex_get_partial_var, exoid, step, EX_NODAL, offset + c, 1, xs / bs + 1, (xe - xs) / bs, varray);
1409: PetscCall(VecRestoreArray(vComp, &varray));
1410: PetscCall(VecRestoreSubVector(vNatural, compIS, &vComp));
1411: }
1412: PetscCall(ISDestroy(&compIS));
1413: }
1414: if (useNatural) {
1415: PetscCall(DMPlexNaturalToGlobalBegin(dm, vNatural, v));
1416: PetscCall(DMPlexNaturalToGlobalEnd(dm, vNatural, v));
1417: PetscCall(VecDestroy(&vNatural));
1418: }
1419: PetscFunctionReturn(PETSC_SUCCESS);
1420: }
1422: static PetscErrorCode VecViewPlex_ExodusII_Zonal_Internal(Vec v, PetscExodusIIInt exoid, PetscExodusIIInt step, PetscExodusIIInt offset)
1423: {
1424: MPI_Comm comm;
1425: PetscMPIInt size;
1426: DM dm;
1427: Vec vNatural, vComp;
1428: const PetscScalar *varray;
1429: PetscInt xs, xe, bs;
1430: PetscBool useNatural;
1431: IS compIS;
1432: PetscInt *csSize, *csID;
1433: PetscExodusIIInt numCS, set, csxs = 0;
1435: PetscFunctionBegin;
1436: PetscCall(PetscObjectGetComm((PetscObject)v, &comm));
1437: PetscCallMPI(MPI_Comm_size(comm, &size));
1438: PetscCall(VecGetDM(v, &dm));
1439: PetscCall(DMGetUseNatural(dm, &useNatural));
1440: useNatural = useNatural && size > 1 ? PETSC_TRUE : PETSC_FALSE;
1441: if (useNatural) {
1442: PetscCall(DMPlexCreateNaturalVector(dm, &vNatural));
1443: PetscCall(DMPlexGlobalToNaturalBegin(dm, v, vNatural));
1444: PetscCall(DMPlexGlobalToNaturalEnd(dm, v, vNatural));
1445: } else {
1446: vNatural = v;
1447: }
1449: /* Write local chunk of the result in the exodus file
1450: exodus stores each component of a vector-valued field as a separate variable.
1451: We assume that they are stored sequentially
1452: Zonal variables are accessed one element block at a time, so we loop through the cell sets,
1453: but once the vector has been reordered to natural size, we cannot use the label information
1454: to figure out what to save where. */
1455: numCS = (PetscExodusIIInt)ex_inquire_int(exoid, EX_INQ_ELEM_BLK); // This is an int64_t
1456: PetscCall(PetscMalloc2(numCS, &csID, numCS, &csSize));
1457: PetscCallExternal(ex_get_ids, exoid, EX_ELEM_BLOCK, csID);
1458: for (set = 0; set < numCS; ++set) {
1459: ex_block block;
1461: block.id = csID[set];
1462: block.type = EX_ELEM_BLOCK;
1463: PetscCallExternal(ex_get_block_param, exoid, &block);
1464: PetscCall(PetscIntCast(block.num_entry, &csSize[set])); // This is an int64_t
1465: }
1466: PetscCall(VecGetOwnershipRange(vNatural, &xs, &xe));
1467: PetscCall(VecGetBlockSize(vNatural, &bs));
1468: if (bs > 1) PetscCall(ISCreateStride(comm, (xe - xs) / bs, xs, bs, &compIS));
1469: for (set = 0; set < numCS; set++) {
1470: PetscInt csLocalSize, c;
1472: /* range of indices for set setID[set]: csxs:csxs + csSize[set]-1
1473: local slice of zonal values: xs/bs,xm/bs-1
1474: intersection: max(xs/bs,csxs),min(xm/bs-1,csxs + csSize[set]-1) */
1475: csLocalSize = PetscMax(0, PetscMin(xe / bs, csxs + csSize[set]) - PetscMax(xs / bs, csxs));
1476: if (bs == 1) {
1477: PetscCall(VecGetArrayRead(vNatural, &varray));
1478: PetscCallExternal(ex_put_partial_var, exoid, step, EX_ELEM_BLOCK, offset, csID[set], PetscMax(xs - csxs, 0) + 1, csLocalSize, &varray[PetscMax(0, csxs - xs)]);
1479: PetscCall(VecRestoreArrayRead(vNatural, &varray));
1480: } else {
1481: for (c = 0; c < bs; ++c) {
1482: PetscCall(ISStrideSetStride(compIS, (xe - xs) / bs, xs + c, bs));
1483: PetscCall(VecGetSubVector(vNatural, compIS, &vComp));
1484: PetscCall(VecGetArrayRead(vComp, &varray));
1485: 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)]);
1486: PetscCall(VecRestoreArrayRead(vComp, &varray));
1487: PetscCall(VecRestoreSubVector(vNatural, compIS, &vComp));
1488: }
1489: }
1490: csxs += csSize[set];
1491: }
1492: PetscCall(PetscFree2(csID, csSize));
1493: if (bs > 1) PetscCall(ISDestroy(&compIS));
1494: if (useNatural) PetscCall(VecDestroy(&vNatural));
1495: PetscFunctionReturn(PETSC_SUCCESS);
1496: }
1498: static PetscErrorCode VecLoadPlex_ExodusII_Zonal_Internal(Vec v, PetscExodusIIInt exoid, PetscExodusIIInt step, PetscExodusIIInt offset)
1499: {
1500: MPI_Comm comm;
1501: PetscMPIInt size;
1502: DM dm;
1503: Vec vNatural, vComp;
1504: PetscScalar *varray;
1505: PetscInt xs, xe, bs;
1506: PetscBool useNatural;
1507: IS compIS;
1508: PetscInt *csSize, *csID;
1509: PetscExodusIIInt numCS, set, csxs = 0;
1511: PetscFunctionBegin;
1512: PetscCall(PetscObjectGetComm((PetscObject)v, &comm));
1513: PetscCallMPI(MPI_Comm_size(comm, &size));
1514: PetscCall(VecGetDM(v, &dm));
1515: PetscCall(DMGetUseNatural(dm, &useNatural));
1516: useNatural = useNatural && size > 1 ? PETSC_TRUE : PETSC_FALSE;
1517: if (useNatural) PetscCall(DMPlexCreateNaturalVector(dm, &vNatural));
1518: else vNatural = v;
1520: /* Read local chunk of the result in the exodus file
1521: exodus stores each component of a vector-valued field as a separate variable.
1522: We assume that they are stored sequentially
1523: Zonal variables are accessed one element block at a time, so we loop through the cell sets,
1524: but once the vector has been reordered to natural size, we cannot use the label information
1525: to figure out what to save where. */
1526: numCS = (PetscExodusIIInt)ex_inquire_int(exoid, EX_INQ_ELEM_BLK); // This is an int64_t
1527: PetscCall(PetscMalloc2(numCS, &csID, numCS, &csSize));
1528: PetscCallExternal(ex_get_ids, exoid, EX_ELEM_BLOCK, csID);
1529: for (set = 0; set < numCS; ++set) {
1530: ex_block block;
1532: block.id = csID[set];
1533: block.type = EX_ELEM_BLOCK;
1534: PetscCallExternal(ex_get_block_param, exoid, &block);
1535: PetscCall(PetscIntCast(block.num_entry, &csSize[set])); // This is an int64_t
1536: }
1537: PetscCall(VecGetOwnershipRange(vNatural, &xs, &xe));
1538: PetscCall(VecGetBlockSize(vNatural, &bs));
1539: if (bs > 1) PetscCall(ISCreateStride(comm, (xe - xs) / bs, xs, bs, &compIS));
1540: for (set = 0; set < numCS; ++set) {
1541: PetscInt csLocalSize, c;
1543: /* range of indices for set setID[set]: csxs:csxs + csSize[set]-1
1544: local slice of zonal values: xs/bs,xm/bs-1
1545: intersection: max(xs/bs,csxs),min(xm/bs-1,csxs + csSize[set]-1) */
1546: csLocalSize = PetscMax(0, PetscMin(xe / bs, csxs + csSize[set]) - PetscMax(xs / bs, csxs));
1547: if (bs == 1) {
1548: PetscCall(VecGetArray(vNatural, &varray));
1549: PetscCallExternal(ex_get_partial_var, exoid, step, EX_ELEM_BLOCK, offset, csID[set], PetscMax(xs - csxs, 0) + 1, csLocalSize, &varray[PetscMax(0, csxs - xs)]);
1550: PetscCall(VecRestoreArray(vNatural, &varray));
1551: } else {
1552: for (c = 0; c < bs; ++c) {
1553: PetscCall(ISStrideSetStride(compIS, (xe - xs) / bs, xs + c, bs));
1554: PetscCall(VecGetSubVector(vNatural, compIS, &vComp));
1555: PetscCall(VecGetArray(vComp, &varray));
1556: 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)]);
1557: PetscCall(VecRestoreArray(vComp, &varray));
1558: PetscCall(VecRestoreSubVector(vNatural, compIS, &vComp));
1559: }
1560: }
1561: csxs += csSize[set];
1562: }
1563: PetscCall(PetscFree2(csID, csSize));
1564: if (bs > 1) PetscCall(ISDestroy(&compIS));
1565: if (useNatural) {
1566: PetscCall(DMPlexNaturalToGlobalBegin(dm, vNatural, v));
1567: PetscCall(DMPlexNaturalToGlobalEnd(dm, vNatural, v));
1568: PetscCall(VecDestroy(&vNatural));
1569: }
1570: PetscFunctionReturn(PETSC_SUCCESS);
1571: }
1573: static PetscErrorCode ExodusGetCellType_Internal(const char *elem_type, DMPolytopeType *ct)
1574: {
1575: PetscBool flg;
1577: PetscFunctionBegin;
1578: *ct = DM_POLYTOPE_UNKNOWN;
1579: PetscCall(PetscStrcmp(elem_type, "BAR2", &flg));
1580: if (flg) {
1581: *ct = DM_POLYTOPE_SEGMENT;
1582: goto done;
1583: }
1584: PetscCall(PetscStrcmp(elem_type, "BAR3", &flg));
1585: if (flg) {
1586: *ct = DM_POLYTOPE_SEGMENT;
1587: goto done;
1588: }
1589: PetscCall(PetscStrcmp(elem_type, "TRI", &flg));
1590: if (flg) {
1591: *ct = DM_POLYTOPE_TRIANGLE;
1592: goto done;
1593: }
1594: PetscCall(PetscStrcmp(elem_type, "TRI3", &flg));
1595: if (flg) {
1596: *ct = DM_POLYTOPE_TRIANGLE;
1597: goto done;
1598: }
1599: PetscCall(PetscStrcmp(elem_type, "QUAD", &flg));
1600: if (flg) {
1601: *ct = DM_POLYTOPE_QUADRILATERAL;
1602: goto done;
1603: }
1604: PetscCall(PetscStrcmp(elem_type, "QUAD4", &flg));
1605: if (flg) {
1606: *ct = DM_POLYTOPE_QUADRILATERAL;
1607: goto done;
1608: }
1609: PetscCall(PetscStrcmp(elem_type, "SHELL4", &flg));
1610: if (flg) {
1611: *ct = DM_POLYTOPE_QUADRILATERAL;
1612: goto done;
1613: }
1614: PetscCall(PetscStrcmp(elem_type, "TETRA", &flg));
1615: if (flg) {
1616: *ct = DM_POLYTOPE_TETRAHEDRON;
1617: goto done;
1618: }
1619: PetscCall(PetscStrcmp(elem_type, "TET4", &flg));
1620: if (flg) {
1621: *ct = DM_POLYTOPE_TETRAHEDRON;
1622: goto done;
1623: }
1624: PetscCall(PetscStrcmp(elem_type, "WEDGE", &flg));
1625: if (flg) {
1626: *ct = DM_POLYTOPE_TRI_PRISM;
1627: goto done;
1628: }
1629: PetscCall(PetscStrcmp(elem_type, "HEX", &flg));
1630: if (flg) {
1631: *ct = DM_POLYTOPE_HEXAHEDRON;
1632: goto done;
1633: }
1634: PetscCall(PetscStrcmp(elem_type, "HEX8", &flg));
1635: if (flg) {
1636: *ct = DM_POLYTOPE_HEXAHEDRON;
1637: goto done;
1638: }
1639: PetscCall(PetscStrcmp(elem_type, "HEXAHEDRON", &flg));
1640: if (flg) {
1641: *ct = DM_POLYTOPE_HEXAHEDRON;
1642: goto done;
1643: }
1644: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unrecognized element type %s", elem_type);
1645: done:
1646: PetscFunctionReturn(PETSC_SUCCESS);
1647: }
1649: /*@
1650: DMPlexCreateExodus - Create a `DMPLEX` mesh from an ExodusII file ID.
1652: Collective
1654: Input Parameters:
1655: + comm - The MPI communicator
1656: . exoid - The ExodusII id associated with a exodus file and obtained using ex_open
1657: - interpolate - Create faces and edges in the mesh
1659: Output Parameter:
1660: . dm - The `DM` object representing the mesh
1662: Level: beginner
1664: .seealso: [](ch_unstructured), `DM`, `PETSCVIEWEREXODUSII`, `DMPLEX`, `DMCreate()`
1665: @*/
1666: PetscErrorCode DMPlexCreateExodus(MPI_Comm comm, PetscExodusIIInt exoid, PetscBool interpolate, DM *dm)
1667: {
1668: PetscMPIInt num_proc, rank;
1669: DMLabel cellSets = NULL, faceSets = NULL, vertSets = NULL;
1670: PetscSection coordSection;
1671: Vec coordinates;
1672: PetscScalar *coords;
1673: PetscInt coordSize, v;
1674: /* Read from ex_get_init() */
1675: char title[PETSC_MAX_PATH_LEN + 1];
1676: int dim = 0, dimEmbed = 0, numVertices = 0, numCells = 0;
1677: int num_cs = 0, num_vs = 0, num_fs = 0;
1679: PetscFunctionBegin;
1680: PetscCallMPI(MPI_Comm_rank(comm, &rank));
1681: PetscCallMPI(MPI_Comm_size(comm, &num_proc));
1682: PetscCall(DMCreate(comm, dm));
1683: PetscCall(DMSetType(*dm, DMPLEX));
1684: /* Open EXODUS II file and read basic information on rank 0, then broadcast to all processors */
1685: if (rank == 0) {
1686: PetscCall(PetscMemzero(title, PETSC_MAX_PATH_LEN + 1));
1687: PetscCallExternal(ex_get_init, exoid, title, &dimEmbed, &numVertices, &numCells, &num_cs, &num_vs, &num_fs);
1688: PetscCheck(num_cs, PETSC_COMM_SELF, PETSC_ERR_SUP, "Exodus file does not contain any cell set");
1689: }
1690: PetscCallMPI(MPI_Bcast(title, PETSC_MAX_PATH_LEN + 1, MPI_CHAR, 0, comm));
1691: PetscCallMPI(MPI_Bcast(&dim, 1, MPI_INT, 0, comm));
1692: PetscCall(PetscObjectSetName((PetscObject)*dm, title));
1693: PetscCall(DMPlexSetChart(*dm, 0, numCells + numVertices));
1694: /* We do not want this label automatically computed, instead we compute it here */
1695: PetscCall(DMCreateLabel(*dm, "celltype"));
1697: /* Read cell sets information */
1698: if (rank == 0) {
1699: PetscInt *cone;
1700: int c, cs, ncs, c_loc, v, v_loc;
1701: /* Read from ex_get_elem_blk_ids() */
1702: int *cs_id, *cs_order;
1703: /* Read from ex_get_elem_block() */
1704: char buffer[PETSC_MAX_PATH_LEN + 1];
1705: int num_cell_in_set, num_vertex_per_cell, num_hybrid, num_attr;
1706: /* Read from ex_get_elem_conn() */
1707: int *cs_connect;
1709: /* Get cell sets IDs */
1710: PetscCall(PetscMalloc2(num_cs, &cs_id, num_cs, &cs_order));
1711: PetscCallExternal(ex_get_ids, exoid, EX_ELEM_BLOCK, cs_id);
1712: /* Read the cell set connectivity table and build mesh topology
1713: EXO standard requires that cells in cell sets be numbered sequentially and be pairwise disjoint. */
1714: /* Check for a hybrid mesh */
1715: for (cs = 0, num_hybrid = 0; cs < num_cs; ++cs) {
1716: DMPolytopeType ct;
1717: char elem_type[PETSC_MAX_PATH_LEN];
1719: PetscCall(PetscArrayzero(elem_type, sizeof(elem_type)));
1720: PetscCallExternal(ex_get_elem_type, exoid, cs_id[cs], elem_type);
1721: PetscCall(ExodusGetCellType_Internal(elem_type, &ct));
1722: dim = PetscMax(dim, DMPolytopeTypeGetDim(ct));
1723: PetscCallExternal(ex_get_block, exoid, EX_ELEM_BLOCK, cs_id[cs], buffer, &num_cell_in_set, &num_vertex_per_cell, 0, 0, &num_attr);
1724: switch (ct) {
1725: case DM_POLYTOPE_TRI_PRISM:
1726: cs_order[cs] = cs;
1727: ++num_hybrid;
1728: break;
1729: default:
1730: for (c = cs; c > cs - num_hybrid; --c) cs_order[c] = cs_order[c - 1];
1731: cs_order[cs - num_hybrid] = cs;
1732: }
1733: }
1734: /* First set sizes */
1735: for (ncs = 0, c = 0; ncs < num_cs; ++ncs) {
1736: DMPolytopeType ct;
1737: char elem_type[PETSC_MAX_PATH_LEN];
1738: const PetscInt cs = cs_order[ncs];
1740: PetscCall(PetscArrayzero(elem_type, sizeof(elem_type)));
1741: PetscCallExternal(ex_get_elem_type, exoid, cs_id[cs], elem_type);
1742: PetscCall(ExodusGetCellType_Internal(elem_type, &ct));
1743: PetscCallExternal(ex_get_block, exoid, EX_ELEM_BLOCK, cs_id[cs], buffer, &num_cell_in_set, &num_vertex_per_cell, 0, 0, &num_attr);
1744: for (c_loc = 0; c_loc < num_cell_in_set; ++c_loc, ++c) {
1745: PetscCall(DMPlexSetConeSize(*dm, c, num_vertex_per_cell));
1746: PetscCall(DMPlexSetCellType(*dm, c, ct));
1747: }
1748: }
1749: for (v = numCells; v < numCells + numVertices; ++v) PetscCall(DMPlexSetCellType(*dm, v, DM_POLYTOPE_POINT));
1750: PetscCall(DMSetUp(*dm));
1751: for (ncs = 0, c = 0; ncs < num_cs; ++ncs) {
1752: const PetscInt cs = cs_order[ncs];
1753: PetscCallExternal(ex_get_block, exoid, EX_ELEM_BLOCK, cs_id[cs], buffer, &num_cell_in_set, &num_vertex_per_cell, 0, 0, &num_attr);
1754: PetscCall(PetscMalloc2(num_vertex_per_cell * num_cell_in_set, &cs_connect, num_vertex_per_cell, &cone));
1755: PetscCallExternal(ex_get_conn, exoid, EX_ELEM_BLOCK, cs_id[cs], cs_connect, NULL, NULL);
1756: /* EXO uses Fortran-based indexing, DMPlex uses C-style and numbers cell first then vertices. */
1757: for (c_loc = 0, v = 0; c_loc < num_cell_in_set; ++c_loc, ++c) {
1758: DMPolytopeType ct;
1760: for (v_loc = 0; v_loc < num_vertex_per_cell; ++v_loc, ++v) cone[v_loc] = cs_connect[v] + numCells - 1;
1761: PetscCall(DMPlexGetCellType(*dm, c, &ct));
1762: PetscCall(DMPlexInvertCell(ct, cone));
1763: PetscCall(DMPlexSetCone(*dm, c, cone));
1764: PetscCall(DMSetLabelValue_Fast(*dm, &cellSets, "Cell Sets", c, cs_id[cs]));
1765: }
1766: PetscCall(PetscFree2(cs_connect, cone));
1767: }
1768: PetscCall(PetscFree2(cs_id, cs_order));
1769: }
1770: {
1771: PetscInt ints[] = {dim, dimEmbed};
1773: PetscCallMPI(MPI_Bcast(ints, 2, MPIU_INT, 0, comm));
1774: PetscCall(DMSetDimension(*dm, ints[0]));
1775: PetscCall(DMSetCoordinateDim(*dm, ints[1]));
1776: dim = ints[0];
1777: dimEmbed = ints[1];
1778: }
1779: PetscCall(DMPlexSymmetrize(*dm));
1780: PetscCall(DMPlexStratify(*dm));
1781: if (interpolate) {
1782: DM idm;
1784: PetscCall(DMPlexInterpolate(*dm, &idm));
1785: PetscCall(DMDestroy(dm));
1786: *dm = idm;
1787: }
1789: /* Create vertex set label */
1790: if (rank == 0 && (num_vs > 0)) {
1791: int vs, v;
1792: /* Read from ex_get_node_set_ids() */
1793: int *vs_id;
1794: /* Read from ex_get_node_set_param() */
1795: int num_vertex_in_set;
1796: /* Read from ex_get_node_set() */
1797: int *vs_vertex_list;
1799: /* Get vertex set ids */
1800: PetscCall(PetscMalloc1(num_vs, &vs_id));
1801: PetscCallExternal(ex_get_ids, exoid, EX_NODE_SET, vs_id);
1802: for (vs = 0; vs < num_vs; ++vs) {
1803: PetscCallExternal(ex_get_set_param, exoid, EX_NODE_SET, vs_id[vs], &num_vertex_in_set, NULL);
1804: PetscCall(PetscMalloc1(num_vertex_in_set, &vs_vertex_list));
1805: PetscCallExternal(ex_get_set, exoid, EX_NODE_SET, vs_id[vs], vs_vertex_list, NULL);
1806: 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]));
1807: PetscCall(PetscFree(vs_vertex_list));
1808: }
1809: PetscCall(PetscFree(vs_id));
1810: }
1811: /* Read coordinates */
1812: PetscCall(DMGetCoordinateSection(*dm, &coordSection));
1813: PetscCall(PetscSectionSetNumFields(coordSection, 1));
1814: PetscCall(PetscSectionSetFieldComponents(coordSection, 0, dimEmbed));
1815: PetscCall(PetscSectionSetChart(coordSection, numCells, numCells + numVertices));
1816: for (v = numCells; v < numCells + numVertices; ++v) {
1817: PetscCall(PetscSectionSetDof(coordSection, v, dimEmbed));
1818: PetscCall(PetscSectionSetFieldDof(coordSection, v, 0, dimEmbed));
1819: }
1820: PetscCall(PetscSectionSetUp(coordSection));
1821: PetscCall(PetscSectionGetStorageSize(coordSection, &coordSize));
1822: PetscCall(VecCreate(PETSC_COMM_SELF, &coordinates));
1823: PetscCall(PetscObjectSetName((PetscObject)coordinates, "coordinates"));
1824: PetscCall(VecSetSizes(coordinates, coordSize, PETSC_DETERMINE));
1825: PetscCall(VecSetBlockSize(coordinates, dimEmbed));
1826: PetscCall(VecSetType(coordinates, VECSTANDARD));
1827: PetscCall(VecGetArray(coordinates, &coords));
1828: if (rank == 0) {
1829: PetscReal *x, *y, *z;
1831: PetscCall(PetscMalloc3(numVertices, &x, numVertices, &y, numVertices, &z));
1832: PetscCallExternal(ex_get_coord, exoid, x, y, z);
1833: if (dimEmbed > 0) {
1834: for (v = 0; v < numVertices; ++v) coords[v * dimEmbed + 0] = x[v];
1835: }
1836: if (dimEmbed > 1) {
1837: for (v = 0; v < numVertices; ++v) coords[v * dimEmbed + 1] = y[v];
1838: }
1839: if (dimEmbed > 2) {
1840: for (v = 0; v < numVertices; ++v) coords[v * dimEmbed + 2] = z[v];
1841: }
1842: PetscCall(PetscFree3(x, y, z));
1843: }
1844: PetscCall(VecRestoreArray(coordinates, &coords));
1845: PetscCall(DMSetCoordinatesLocal(*dm, coordinates));
1846: PetscCall(VecDestroy(&coordinates));
1848: /* Create side set label */
1849: if (rank == 0 && interpolate && (num_fs > 0)) {
1850: int fs, f, voff;
1851: /* Read from ex_get_side_set_ids() */
1852: int *fs_id;
1853: /* Read from ex_get_side_set_param() */
1854: int num_side_in_set;
1855: /* Read from ex_get_side_set_node_list() */
1856: int *fs_vertex_count_list, *fs_vertex_list, *fs_side_list;
1857: /* Read side set labels */
1858: char fs_name[MAX_STR_LENGTH + 1];
1859: size_t fs_name_len;
1861: /* Get side set ids */
1862: PetscCall(PetscMalloc1(num_fs, &fs_id));
1863: PetscCallExternal(ex_get_ids, exoid, EX_SIDE_SET, fs_id);
1864: // Ids 1 and 2 are reserved by ExodusII for indicating things in 3D
1865: for (fs = 0; fs < num_fs; ++fs) {
1866: PetscCallExternal(ex_get_set_param, exoid, EX_SIDE_SET, fs_id[fs], &num_side_in_set, NULL);
1867: 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));
1868: PetscCallExternal(ex_get_side_set_node_list, exoid, fs_id[fs], fs_vertex_count_list, fs_vertex_list);
1869: PetscCallExternal(ex_get_set, exoid, EX_SIDE_SET, fs_id[fs], NULL, fs_side_list);
1871: /* Get the specific name associated with this side set ID. */
1872: int fs_name_err = ex_get_name(exoid, EX_SIDE_SET, fs_id[fs], fs_name);
1873: if (!fs_name_err) {
1874: PetscCall(PetscStrlen(fs_name, &fs_name_len));
1875: if (fs_name_len == 0) PetscCall(PetscStrncpy(fs_name, "Face Sets", MAX_STR_LENGTH));
1876: }
1877: for (f = 0, voff = 0; f < num_side_in_set; ++f) {
1878: const PetscInt *faces = NULL;
1879: PetscInt faceSize = fs_vertex_count_list[f], numFaces;
1880: PetscInt faceVertices[4], v;
1882: PetscCheck(faceSize <= 4, comm, PETSC_ERR_ARG_WRONG, "ExodusII side cannot have %" PetscInt_FMT " > 4 vertices", faceSize);
1883: for (v = 0; v < faceSize; ++v, ++voff) faceVertices[v] = fs_vertex_list[voff] + numCells - 1;
1884: PetscCall(DMPlexGetFullJoin(*dm, faceSize, faceVertices, &numFaces, &faces));
1885: PetscCheck(numFaces == 1, comm, PETSC_ERR_ARG_WRONG, "Invalid ExodusII side %d in set %d maps to %" PetscInt_FMT " faces", f, fs, numFaces);
1886: 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]);
1887: PetscCall(DMSetLabelValue_Fast(*dm, &faceSets, "Face Sets", faces[0], fs_id[fs]));
1888: /* Only add the label if one has been detected for this side set. */
1889: if (!fs_name_err) PetscCall(DMSetLabelValue(*dm, fs_name, faces[0], fs_id[fs]));
1890: PetscCall(DMPlexRestoreJoin(*dm, faceSize, faceVertices, &numFaces, &faces));
1891: }
1892: PetscCall(PetscFree3(fs_vertex_count_list, fs_vertex_list, fs_side_list));
1893: }
1894: PetscCall(PetscFree(fs_id));
1895: }
1897: { /* Create Cell/Face/Vertex Sets labels at all processes */
1898: enum {
1899: n = 3
1900: };
1901: PetscBool flag[n];
1903: flag[0] = cellSets ? PETSC_TRUE : PETSC_FALSE;
1904: flag[1] = faceSets ? PETSC_TRUE : PETSC_FALSE;
1905: flag[2] = vertSets ? PETSC_TRUE : PETSC_FALSE;
1906: PetscCallMPI(MPI_Bcast(flag, n, MPIU_BOOL, 0, comm));
1907: if (flag[0]) PetscCall(DMCreateLabel(*dm, "Cell Sets"));
1908: if (flag[1]) PetscCall(DMCreateLabel(*dm, "Face Sets"));
1909: if (flag[2]) PetscCall(DMCreateLabel(*dm, "Vertex Sets"));
1910: }
1911: PetscFunctionReturn(PETSC_SUCCESS);
1912: }