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