Actual source code: plexexodusii.c
petsc-3.13.6 2020-09-29
1: #define PETSCDM_DLL
2: #include <petsc/private/dmpleximpl.h>
4: #if defined(PETSC_HAVE_EXODUSII)
5: #include <netcdf.h>
6: #include <exodusII.h>
7: #endif
9: #include <petsc/private/viewerimpl.h>
10: #include <petsc/private/viewerexodusiiimpl.h>
11: #if defined(PETSC_HAVE_EXODUSII)
12: /*
13: PETSC_VIEWER_EXODUSII_ - Creates an ExodusII PetscViewer shared by all processors in a communicator.
15: Collective
17: Input Parameter:
18: . comm - the MPI communicator to share the ExodusII PetscViewer
20: Level: intermediate
22: Notes:
23: misses Fortran bindings
25: Notes:
26: Unlike almost all other PETSc routines, PETSC_VIEWER_EXODUSII_ does not return
27: an error code. The GLVIS PetscViewer is usually used in the form
28: $ XXXView(XXX object, PETSC_VIEWER_EXODUSII_(comm));
30: .seealso: PetscViewerExodusIIOpen(), PetscViewerType, PetscViewerCreate(), PetscViewerDestroy()
31: */
32: PetscViewer PETSC_VIEWER_EXODUSII_(MPI_Comm comm)
33: {
34: PetscViewer viewer;
38: PetscViewerExodusIIOpen(comm, "mesh.exo", FILE_MODE_WRITE, &viewer);
39: if (ierr) {PetscError(PETSC_COMM_SELF,__LINE__,"PETSC_VIEWER_EXODUSII_",__FILE__,PETSC_ERR_PLIB,PETSC_ERROR_INITIAL," ");return(0);}
40: PetscObjectRegisterDestroy((PetscObject) viewer);
41: if (ierr) {PetscError(PETSC_COMM_SELF,__LINE__,"PETSC_VIEWER_EXODUSII_",__FILE__,PETSC_ERR_PLIB,PETSC_ERROR_INITIAL," ");return(0);}
42: PetscFunctionReturn(viewer);
43: }
45: static PetscErrorCode PetscViewerView_ExodusII(PetscViewer v, PetscViewer viewer)
46: {
47: PetscViewer_ExodusII *exo = (PetscViewer_ExodusII *) viewer->data;
48: PetscErrorCode ierr;
51: if (exo->filename) {PetscViewerASCIIPrintf(viewer, "Filename: %s\n", exo->filename);}
52: return(0);
53: }
55: static PetscErrorCode PetscViewerSetFromOptions_ExodusII(PetscOptionItems *PetscOptionsObject, PetscViewer v)
56: {
60: PetscOptionsHead(PetscOptionsObject, "ExodusII PetscViewer Options");
61: PetscOptionsTail();
62: return(0);
63: }
65: static PetscErrorCode PetscViewerSetUp_ExodusII(PetscViewer viewer)
66: {
68: return(0);
69: }
71: static PetscErrorCode PetscViewerDestroy_ExodusII(PetscViewer viewer)
72: {
73: PetscViewer_ExodusII *exo = (PetscViewer_ExodusII *) viewer->data;
74: PetscErrorCode ierr;
77: if (exo->exoid >= 0) {ex_close(exo->exoid);}
78: PetscFree(exo);
79: PetscObjectComposeFunction((PetscObject)viewer,"PetscViewerFileSetName_C",NULL);
80: PetscObjectComposeFunction((PetscObject)viewer,"PetscViewerFileGetName_C",NULL);
81: PetscObjectComposeFunction((PetscObject)viewer,"PetscViewerFileSetMode_C",NULL);
82: PetscObjectComposeFunction((PetscObject)viewer,"PetscViewerExodusIIGetId",NULL);
83: return(0);
84: }
86: static PetscErrorCode PetscViewerFileSetName_ExodusII(PetscViewer viewer, const char name[])
87: {
88: PetscViewer_ExodusII *exo = (PetscViewer_ExodusII *) viewer->data;
89: PetscMPIInt rank;
90: int CPU_word_size, IO_word_size, EXO_mode;
91: PetscErrorCode ierr;
94: MPI_Comm_rank(PetscObjectComm((PetscObject) viewer), &rank);
95: CPU_word_size = sizeof(PetscReal);
96: IO_word_size = sizeof(PetscReal);
99: if (exo->exoid >= 0) {ex_close(exo->exoid); exo->exoid = -1;}
100: if (exo->filename) {PetscFree(exo->filename);}
101: PetscStrallocpy(name, &exo->filename);
102: /* Create or open the file collectively */
103: switch (exo->btype) {
104: case FILE_MODE_READ:
105: EXO_mode = EX_CLOBBER;
106: break;
107: case FILE_MODE_APPEND:
108: EXO_mode = EX_CLOBBER;
109: break;
110: case FILE_MODE_WRITE:
111: EXO_mode = EX_CLOBBER;
112: break;
113: default:
114: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER, "Must call PetscViewerFileSetMode() before PetscViewerFileSetName()");
115: }
116: #if defined(PETSC_USE_64BIT_INDICES)
117: EXO_mode += EX_ALL_INT64_API;
118: #endif
119: exo->exoid = ex_create(name, EXO_mode, &CPU_word_size, &IO_word_size);
120: if (exo->exoid < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_LIB, "ex_create failed for %s", name);
121: return(0);
122: }
124: static PetscErrorCode PetscViewerFileGetName_ExodusII(PetscViewer viewer, const char **name)
125: {
126: PetscViewer_ExodusII *exo = (PetscViewer_ExodusII *) viewer->data;
129: *name = exo->filename;
130: return(0);
131: }
133: static PetscErrorCode PetscViewerFileSetMode_ExodusII(PetscViewer viewer, PetscFileMode type)
134: {
135: PetscViewer_ExodusII *exo = (PetscViewer_ExodusII *) viewer->data;
138: exo->btype = type;
139: return(0);
140: }
142: static PetscErrorCode PetscViewerFileGetMode_ExodusII(PetscViewer viewer, PetscFileMode *type)
143: {
144: PetscViewer_ExodusII *exo = (PetscViewer_ExodusII *) viewer->data;
147: *type = exo->btype;
148: return(0);
149: }
151: static PetscErrorCode PetscViewerExodusIIGetId_ExodusII(PetscViewer viewer, int *exoid)
152: {
153: PetscViewer_ExodusII *exo = (PetscViewer_ExodusII *) viewer->data;
156: *exoid = exo->exoid;
157: return(0);
158: }
160: /*@C
161: PetscViewerExodusIIOpen - Opens a file for ExodusII input/output.
163: Collective
165: Input Parameters:
166: + comm - MPI communicator
167: . name - name of file
168: - type - type of file
169: $ FILE_MODE_WRITE - create new file for binary output
170: $ FILE_MODE_READ - open existing file for binary input
171: $ FILE_MODE_APPEND - open existing file for binary output
173: Output Parameter:
174: . exo - PetscViewer for Exodus II input/output to use with the specified file
176: Level: beginner
178: Note:
179: This PetscViewer should be destroyed with PetscViewerDestroy().
182: .seealso: PetscViewerPushFormat(), PetscViewerDestroy(),
183: DMLoad(), PetscFileMode, PetscViewer, PetscViewerSetType(), PetscViewerFileSetMode(), PetscViewerFileSetName()
184: @*/
185: PetscErrorCode PetscViewerExodusIIOpen(MPI_Comm comm, const char name[], PetscFileMode type, PetscViewer *exo)
186: {
190: PetscViewerCreate(comm, exo);
191: PetscViewerSetType(*exo, PETSCVIEWEREXODUSII);
192: PetscViewerFileSetMode(*exo, type);
193: PetscViewerFileSetName(*exo, name);
194: PetscViewerSetFromOptions(*exo);
195: return(0);
196: }
198: /*MC
199: PETSCVIEWEREXODUSII - A viewer that writes to an Exodus II file
202: .seealso: PetscViewerExodusIIOpen(), PetscViewerCreate(), PETSCVIEWERBINARY, PETSCVIEWERHDF5, DMView(),
203: PetscViewerFileSetName(), PetscViewerFileSetMode(), PetscViewerFormat, PetscViewerType, PetscViewerSetType()
205: Level: beginner
206: M*/
208: PETSC_EXTERN PetscErrorCode PetscViewerCreate_ExodusII(PetscViewer v)
209: {
210: PetscViewer_ExodusII *exo;
211: PetscErrorCode ierr;
214: PetscNewLog(v,&exo);
216: v->data = (void*) exo;
217: v->ops->destroy = PetscViewerDestroy_ExodusII;
218: v->ops->setfromoptions = PetscViewerSetFromOptions_ExodusII;
219: v->ops->setup = PetscViewerSetUp_ExodusII;
220: v->ops->view = PetscViewerView_ExodusII;
221: v->ops->flush = 0;
222: exo->btype = (PetscFileMode) -1;
223: exo->filename = 0;
224: exo->exoid = -1;
226: PetscObjectComposeFunction((PetscObject)v,"PetscViewerFileSetName_C",PetscViewerFileSetName_ExodusII);
227: PetscObjectComposeFunction((PetscObject)v,"PetscViewerFileGetName_C",PetscViewerFileGetName_ExodusII);
228: PetscObjectComposeFunction((PetscObject)v,"PetscViewerFileSetMode_C",PetscViewerFileSetMode_ExodusII);
229: PetscObjectComposeFunction((PetscObject)v,"PetscViewerFileGetMode_C",PetscViewerFileGetMode_ExodusII);
230: PetscObjectComposeFunction((PetscObject)v,"PetscViewerExodusIIGetId_C",PetscViewerExodusIIGetId_ExodusII);
231: return(0);
232: }
234: /*
235: EXOGetVarIndex - Locate a result in an exodus file based on its name
237: Not collective
239: Input Parameters:
240: + exoid - the exodus id of a file (obtained from ex_open or ex_create for instance)
241: . obj_type - the type of entity for instance EX_NODAL, EX_ELEM_BLOCK
242: - name - the name of the result
244: Output Parameters:
245: . varIndex - the location in the exodus file of the result
247: Notes:
248: The exodus variable index is obtained by comparing name and the
249: names of zonal variables declared in the exodus file. For instance if name is "V"
250: the location in the exodus file will be the first match of "V", "V_X", "V_XX", "V_1", or "V_11"
251: amongst all variables of type obj_type.
253: Level: beginner
255: .seealso: EXOGetVarIndex(),DMPlexView_ExodusII_Internal(),VecViewPlex_ExodusII_Nodal_Internal(),VecLoadNodal_PlexEXO(),VecLoadZonal_PlexEXO()
256: */
257: static PetscErrorCode EXOGetVarIndex_Private(int exoid, ex_entity_type obj_type, const char name[], int *varIndex)
258: {
259: int num_vars, i, j;
260: char ext_name[MAX_STR_LENGTH+1], var_name[MAX_STR_LENGTH+1];
261: const int num_suffix = 5;
262: char *suffix[5];
263: PetscBool flg;
267: suffix[0] = (char *) "";
268: suffix[1] = (char *) "_X";
269: suffix[2] = (char *) "_XX";
270: suffix[3] = (char *) "_1";
271: suffix[4] = (char *) "_11";
273: *varIndex = 0;
274: PetscStackCallStandard(ex_get_variable_param,(exoid, obj_type, &num_vars));
275: for (i = 0; i < num_vars; ++i) {
276: PetscStackCallStandard(ex_get_variable_name,(exoid, obj_type, i+1, var_name));
277: for (j = 0; j < num_suffix; ++j){
278: PetscStrncpy(ext_name, name, MAX_STR_LENGTH);
279: PetscStrlcat(ext_name, suffix[j], MAX_STR_LENGTH);
280: PetscStrcasecmp(ext_name, var_name, &flg);
281: if (flg) {
282: *varIndex = i+1;
283: return(0);
284: }
285: }
286: }
287: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unable to locate variable %s in ExodusII file.", name);
288: PetscFunctionReturn(-1);
289: }
291: /*
292: DMPlexView_ExodusII_Internal - Write a DM to disk in exodus format
294: Collective on dm
296: Input Parameters:
297: + dm - The dm to be written
298: . exoid - the exodus id of a file (obtained from ex_open or ex_create for instance)
299: - degree - the degree of the interpolation space
301: Notes:
302: Not all DM can be written to disk this way. For instance, exodus assume that element blocks (mapped to "Cell sets" labels)
303: consists of sequentially numbered cells. If this is not the case, the exodus file will be corrupted.
305: If the dm has been distributed and exoid points to different files on each MPI rank, only the "local" part of the DM
306: (including "ghost" cells and vertices) will be written. If all exoid points to the same file, the resulting file will
307: probably be corrupted.
309: DMPlex only represents geometry while most post-processing software expect that a mesh also provides information
310: on the discretization space. This function assumes that the file represents Lagrange finite elements of order 1 or 2.
311: It should be extended to use PetscFE objects.
313: This function will only handle TRI, TET, QUAD and HEX cells.
314: Level: beginner
316: .seealso: EXOGetVarIndex_Private(), VecViewPlex_ExodusII_Nodal_Internal(), VecLoadNodal_PlexEXO(), VecViewZonal_PlexEXO(), VecLoadZonal_PlexEXO()
317: */
318: PetscErrorCode DMPlexView_ExodusII_Internal(DM dm, int exoid, PetscInt degree)
319: {
320: enum ElemType {TRI, QUAD, TET, HEX};
321: MPI_Comm comm;
322: /* Connectivity Variables */
323: PetscInt cellsNotInConnectivity;
324: /* Cell Sets */
325: DMLabel csLabel;
326: IS csIS;
327: const PetscInt *csIdx;
328: PetscInt num_cs, cs;
329: enum ElemType *type;
330: PetscBool hasLabel;
331: /* Coordinate Variables */
332: DM cdm;
333: PetscSection section;
334: Vec coord;
335: PetscInt **nodes;
336: PetscInt depth, d, dim, skipCells = 0;
337: PetscInt pStart, pEnd, p, cStart, cEnd, numCells, vStart, vEnd, numVertices, eStart, eEnd, numEdges, fStart, fEnd, numFaces, numNodes;
338: PetscInt num_vs, num_fs;
339: PetscMPIInt rank, size;
340: const char *dmName;
341: PetscInt nodesTriP1[4] = {3,0,0,0};
342: PetscInt nodesTriP2[4] = {3,3,0,0};
343: PetscInt nodesQuadP1[4] = {4,0,0,0};
344: PetscInt nodesQuadP2[4] = {4,4,0,1};
345: PetscInt nodesTetP1[4] = {4,0,0,0};
346: PetscInt nodesTetP2[4] = {4,6,0,0};
347: PetscInt nodesHexP1[4] = {8,0,0,0};
348: PetscInt nodesHexP2[4] = {8,12,6,1};
349: PetscErrorCode ierr;
352: PetscObjectGetComm((PetscObject)dm, &comm);
353: MPI_Comm_rank(comm, &rank);
354: MPI_Comm_size(comm, &size);
355: /* --- Get DM info --- */
356: PetscObjectGetName((PetscObject) dm, &dmName);
357: DMPlexGetDepth(dm, &depth);
358: DMGetDimension(dm, &dim);
359: DMPlexGetChart(dm, &pStart, &pEnd);
360: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
361: DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);
362: DMPlexGetDepthStratum(dm, 1, &eStart, &eEnd);
363: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
364: numCells = cEnd - cStart;
365: numEdges = eEnd - eStart;
366: numVertices = vEnd - vStart;
367: if (depth == 3) {numFaces = fEnd - fStart;}
368: else {numFaces = 0;}
369: DMGetLabelSize(dm, "Cell Sets", &num_cs);
370: DMGetLabelSize(dm, "Vertex Sets", &num_vs);
371: DMGetLabelSize(dm, "Face Sets", &num_fs);
372: DMGetCoordinatesLocal(dm, &coord);
373: DMGetCoordinateDM(dm, &cdm);
374: if (num_cs > 0) {
375: DMGetLabel(dm, "Cell Sets", &csLabel);
376: DMLabelGetValueIS(csLabel, &csIS);
377: ISGetIndices(csIS, &csIdx);
378: }
379: PetscMalloc1(num_cs, &nodes);
380: /* Set element type for each block and compute total number of nodes */
381: PetscMalloc1(num_cs, &type);
382: numNodes = numVertices;
383: if (degree == 2) {numNodes += numEdges;}
384: cellsNotInConnectivity = numCells;
385: for (cs = 0; cs < num_cs; ++cs) {
386: IS stratumIS;
387: const PetscInt *cells;
388: PetscScalar *xyz = NULL;
389: PetscInt csSize, closureSize;
391: DMLabelGetStratumIS(csLabel, csIdx[cs], &stratumIS);
392: ISGetIndices(stratumIS, &cells);
393: ISGetSize(stratumIS, &csSize);
394: DMPlexVecGetClosure(cdm, NULL, coord, cells[0], &closureSize, &xyz);
395: switch (dim) {
396: case 2:
397: if (closureSize == 3*dim) {type[cs] = TRI;}
398: else if (closureSize == 4*dim) {type[cs] = QUAD;}
399: else SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Number of vertices %D in dimension %D has no ExodusII type", closureSize/dim, dim);
400: break;
401: case 3:
402: if (closureSize == 4*dim) {type[cs] = TET;}
403: else if (closureSize == 8*dim) {type[cs] = HEX;}
404: else SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Number of vertices %D in dimension %D has no ExodusII type", closureSize/dim, dim);
405: break;
406: default: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Dimension %D not handled by ExodusII viewer", dim);
407: }
408: if ((degree == 2) && (type[cs] == QUAD)) {numNodes += csSize;}
409: if ((degree == 2) && (type[cs] == HEX)) {numNodes += csSize; numNodes += numFaces;}
410: DMPlexVecRestoreClosure(cdm, NULL, coord, cells[0], &closureSize, &xyz);
411: /* Set nodes and Element type */
412: if (type[cs] == TRI) {
413: if (degree == 1) nodes[cs] = nodesTriP1;
414: else if (degree == 2) nodes[cs] = nodesTriP2;
415: } else if (type[cs] == QUAD) {
416: if (degree == 1) nodes[cs] = nodesQuadP1;
417: else if (degree == 2) nodes[cs] = nodesQuadP2;
418: } else if (type[cs] == TET) {
419: if (degree == 1) nodes[cs] = nodesTetP1;
420: else if (degree == 2) nodes[cs] = nodesTetP2;
421: } else if (type[cs] == HEX) {
422: if (degree == 1) nodes[cs] = nodesHexP1;
423: else if (degree == 2) nodes[cs] = nodesHexP2;
424: }
425: /* Compute the number of cells not in the connectivity table */
426: cellsNotInConnectivity -= nodes[cs][3]*csSize;
428: ISRestoreIndices(stratumIS, &cells);
429: ISDestroy(&stratumIS);
430: }
431: if (num_cs > 0) {PetscStackCallStandard(ex_put_init,(exoid, dmName, dim, numNodes, numCells, num_cs, num_vs, num_fs));}
432: /* --- Connectivity --- */
433: for (cs = 0; cs < num_cs; ++cs) {
434: IS stratumIS;
435: const PetscInt *cells;
436: PetscInt *connect, off = 0;
437: PetscInt edgesInClosure = 0, facesInClosure = 0, verticesInClosure = 0;
438: PetscInt csSize, c, connectSize, closureSize;
439: char *elem_type = NULL;
440: char elem_type_tri3[] = "TRI3", elem_type_quad4[] = "QUAD4";
441: char elem_type_tri6[] = "TRI6", elem_type_quad9[] = "QUAD9";
442: char elem_type_tet4[] = "TET4", elem_type_hex8[] = "HEX8";
443: char elem_type_tet10[] = "TET10", elem_type_hex27[] = "HEX27";
445: DMLabelGetStratumIS(csLabel, csIdx[cs], &stratumIS);
446: ISGetIndices(stratumIS, &cells);
447: ISGetSize(stratumIS, &csSize);
448: /* Set Element type */
449: if (type[cs] == TRI) {
450: if (degree == 1) elem_type = elem_type_tri3;
451: else if (degree == 2) elem_type = elem_type_tri6;
452: } else if (type[cs] == QUAD) {
453: if (degree == 1) elem_type = elem_type_quad4;
454: else if (degree == 2) elem_type = elem_type_quad9;
455: } else if (type[cs] == TET) {
456: if (degree == 1) elem_type = elem_type_tet4;
457: else if (degree == 2) elem_type = elem_type_tet10;
458: } else if (type[cs] == HEX) {
459: if (degree == 1) elem_type = elem_type_hex8;
460: else if (degree == 2) elem_type = elem_type_hex27;
461: }
462: connectSize = nodes[cs][0] + nodes[cs][1] + nodes[cs][2] + nodes[cs][3];
463: PetscMalloc1(PetscMax(27,connectSize)*csSize, &connect);
464: PetscStackCallStandard(ex_put_block,(exoid, EX_ELEM_BLOCK, csIdx[cs], elem_type, csSize, connectSize, 0, 0, 1));
465: /* Find number of vertices, edges, and faces in the closure */
466: verticesInClosure = nodes[cs][0];
467: if (depth > 1) {
468: if (dim == 2) {
469: DMPlexGetConeSize(dm, cells[0], &edgesInClosure);
470: } else if (dim == 3) {
471: PetscInt *closure = NULL;
473: DMPlexGetConeSize(dm, cells[0], &facesInClosure);
474: DMPlexGetTransitiveClosure(dm, cells[0], PETSC_TRUE, &closureSize, &closure);
475: edgesInClosure = closureSize - facesInClosure - 1 - verticesInClosure;
476: DMPlexRestoreTransitiveClosure(dm, cells[0], PETSC_TRUE, &closureSize, &closure);
477: }
478: }
479: /* Get connectivity for each cell */
480: for (c = 0; c < csSize; ++c) {
481: PetscInt *closure = NULL;
482: PetscInt temp, i;
484: DMPlexGetTransitiveClosure(dm, cells[c], PETSC_TRUE, &closureSize, &closure);
485: for (i = 0; i < connectSize; ++i) {
486: if (i < nodes[cs][0]) {/* Vertices */
487: connect[i+off] = closure[(i+edgesInClosure+facesInClosure+1)*2] + 1;
488: connect[i+off] -= cellsNotInConnectivity;
489: } else if (i < nodes[cs][0]+nodes[cs][1]) { /* Edges */
490: connect[i+off] = closure[(i-verticesInClosure+facesInClosure+1)*2] + 1;
491: if (nodes[cs][2] == 0) connect[i+off] -= numFaces;
492: connect[i+off] -= cellsNotInConnectivity;
493: } else if (i < nodes[cs][0]+nodes[cs][1]+nodes[cs][3]){ /* Cells */
494: connect[i+off] = closure[0] + 1;
495: connect[i+off] -= skipCells;
496: } else if (i < nodes[cs][0]+nodes[cs][1]+nodes[cs][3]+nodes[cs][2]){ /* Faces */
497: connect[i+off] = closure[(i-edgesInClosure-verticesInClosure)*2] + 1;
498: connect[i+off] -= cellsNotInConnectivity;
499: } else {
500: connect[i+off] = -1;
501: }
502: }
503: /* Tetrahedra are inverted */
504: if (type[cs] == TET) {
505: temp = connect[0+off]; connect[0+off] = connect[1+off]; connect[1+off] = temp;
506: if (degree == 2) {
507: temp = connect[5+off]; connect[5+off] = connect[6+off]; connect[6+off] = temp;
508: temp = connect[7+off]; connect[7+off] = connect[8+off]; connect[8+off] = temp;
509: }
510: }
511: /* Hexahedra are inverted */
512: if (type[cs] == HEX) {
513: temp = connect[1+off]; connect[1+off] = connect[3+off]; connect[3+off] = temp;
514: if (degree == 2) {
515: temp = connect[8+off]; connect[8+off] = connect[11+off]; connect[11+off] = temp;
516: temp = connect[9+off]; connect[9+off] = connect[10+off]; connect[10+off] = temp;
517: temp = connect[16+off]; connect[16+off] = connect[17+off]; connect[17+off] = temp;
518: temp = connect[18+off]; connect[18+off] = connect[19+off]; connect[19+off] = temp;
520: temp = connect[12+off]; connect[12+off] = connect[16+off]; connect[16+off] = temp;
521: temp = connect[13+off]; connect[13+off] = connect[17+off]; connect[17+off] = temp;
522: temp = connect[14+off]; connect[14+off] = connect[18+off]; connect[18+off] = temp;
523: temp = connect[15+off]; connect[15+off] = connect[19+off]; connect[19+off] = temp;
525: temp = connect[23+off]; connect[23+off] = connect[26+off]; connect[26+off] = temp;
526: temp = connect[24+off]; connect[24+off] = connect[25+off]; connect[25+off] = temp;
527: temp = connect[25+off]; connect[25+off] = connect[26+off]; connect[26+off] = temp;
528: }
529: }
530: off += connectSize;
531: DMPlexRestoreTransitiveClosure(dm, cells[c], PETSC_TRUE, &closureSize, &closure);
532: }
533: PetscStackCallStandard(ex_put_conn,(exoid, EX_ELEM_BLOCK, csIdx[cs], connect, 0, 0));
534: skipCells += (nodes[cs][3] == 0)*csSize;
535: PetscFree(connect);
536: ISRestoreIndices(stratumIS, &cells);
537: ISDestroy(&stratumIS);
538: }
539: PetscFree(type);
540: /* --- Coordinates --- */
541: PetscSectionCreate(comm, §ion);
542: PetscSectionSetChart(section, pStart, pEnd);
543: if (num_cs) {
544: for (d = 0; d < depth; ++d) {
545: DMPlexGetDepthStratum(dm, d, &pStart, &pEnd);
546: for (p = pStart; p < pEnd; ++p) {
547: PetscSectionSetDof(section, p, nodes[0][d] > 0);
548: }
549: }
550: }
551: for (cs = 0; cs < num_cs; ++cs) {
552: IS stratumIS;
553: const PetscInt *cells;
554: PetscInt csSize, c;
556: DMLabelGetStratumIS(csLabel, csIdx[cs], &stratumIS);
557: ISGetIndices(stratumIS, &cells);
558: ISGetSize(stratumIS, &csSize);
559: for (c = 0; c < csSize; ++c) {
560: PetscSectionSetDof(section, cells[c], nodes[cs][3] > 0);
561: }
562: ISRestoreIndices(stratumIS, &cells);
563: ISDestroy(&stratumIS);
564: }
565: if (num_cs > 0) {
566: ISRestoreIndices(csIS, &csIdx);
567: ISDestroy(&csIS);
568: }
569: PetscFree(nodes);
570: PetscSectionSetUp(section);
571: if (numNodes > 0) {
572: const char *coordNames[3] = {"x", "y", "z"};
573: PetscScalar *coords, *closure;
574: PetscReal *cval;
575: PetscInt hasDof, n = 0;
577: /* There can't be more than 24 values in the closure of a point for the coord section */
578: PetscCalloc3(numNodes*3, &coords, dim, &cval, 24, &closure);
579: DMGetCoordinatesLocal(dm, &coord);
580: DMPlexGetChart(dm, &pStart, &pEnd);
581: for (p = pStart; p < pEnd; ++p) {
582: PetscSectionGetDof(section, p, &hasDof);
583: if (hasDof) {
584: PetscInt closureSize = 24, j;
586: DMPlexVecGetClosure(cdm, NULL, coord, p, &closureSize, &closure);
587: for (d = 0; d < dim; ++d) {
588: cval[d] = 0.0;
589: for (j = 0; j < closureSize/dim; j++) cval[d] += closure[j*dim+d];
590: coords[d*numNodes+n] = cval[d] * dim / closureSize;
591: }
592: ++n;
593: }
594: }
595: PetscStackCallStandard(ex_put_coord,(exoid, &coords[0*numNodes], &coords[1*numNodes], &coords[2*numNodes]));
596: PetscFree3(coords, cval, closure);
597: PetscStackCallStandard(ex_put_coord_names,(exoid, (char **) coordNames));
598: }
599: PetscSectionDestroy(§ion);
600: /* --- Node Sets/Vertex Sets --- */
601: DMHasLabel(dm, "Vertex Sets", &hasLabel);
602: if (hasLabel) {
603: PetscInt i, vs, vsSize;
604: const PetscInt *vsIdx, *vertices;
605: PetscInt *nodeList;
606: IS vsIS, stratumIS;
607: DMLabel vsLabel;
608: DMGetLabel(dm, "Vertex Sets", &vsLabel);
609: DMLabelGetValueIS(vsLabel, &vsIS);
610: ISGetIndices(vsIS, &vsIdx);
611: for (vs=0; vs<num_vs; ++vs) {
612: DMLabelGetStratumIS(vsLabel, vsIdx[vs], &stratumIS);
613: ISGetIndices(stratumIS, &vertices);
614: ISGetSize(stratumIS, &vsSize);
615: PetscMalloc1(vsSize, &nodeList);
616: for (i=0; i<vsSize; ++i) {
617: nodeList[i] = vertices[i] - skipCells + 1;
618: }
619: PetscStackCallStandard(ex_put_set_param,(exoid, EX_NODE_SET, vsIdx[vs], vsSize, 0));
620: PetscStackCallStandard(ex_put_set,(exoid, EX_NODE_SET, vsIdx[vs], nodeList, NULL));
621: ISRestoreIndices(stratumIS, &vertices);
622: ISDestroy(&stratumIS);
623: PetscFree(nodeList);
624: }
625: ISRestoreIndices(vsIS, &vsIdx);
626: ISDestroy(&vsIS);
627: }
628: /* --- Side Sets/Face Sets --- */
629: DMHasLabel(dm, "Face Sets", &hasLabel);
630: if (hasLabel) {
631: PetscInt i, j, fs, fsSize;
632: const PetscInt *fsIdx, *faces;
633: IS fsIS, stratumIS;
634: DMLabel fsLabel;
635: PetscInt numPoints, *points;
636: PetscInt elem_list_size = 0;
637: PetscInt *elem_list, *elem_ind, *side_list;
639: DMGetLabel(dm, "Face Sets", &fsLabel);
640: /* Compute size of Node List and Element List */
641: DMLabelGetValueIS(fsLabel, &fsIS);
642: ISGetIndices(fsIS, &fsIdx);
643: for (fs=0; fs<num_fs; ++fs) {
644: DMLabelGetStratumIS(fsLabel, fsIdx[fs], &stratumIS);
645: ISGetSize(stratumIS, &fsSize);
646: elem_list_size += fsSize;
647: ISDestroy(&stratumIS);
648: }
649: PetscMalloc3(num_fs, &elem_ind,
650: elem_list_size, &elem_list,
651: elem_list_size, &side_list);
652: elem_ind[0] = 0;
653: for (fs=0; fs<num_fs; ++fs) {
654: DMLabelGetStratumIS(fsLabel, fsIdx[fs], &stratumIS);
655: ISGetIndices(stratumIS, &faces);
656: ISGetSize(stratumIS, &fsSize);
657: /* Set Parameters */
658: PetscStackCallStandard(ex_put_set_param,(exoid, EX_SIDE_SET, fsIdx[fs], fsSize, 0));
659: /* Indices */
660: if (fs<num_fs-1) {
661: elem_ind[fs+1] = elem_ind[fs] + fsSize;
662: }
664: for (i=0; i<fsSize; ++i) {
665: /* Element List */
666: points = NULL;
667: DMPlexGetTransitiveClosure(dm, faces[i], PETSC_FALSE, &numPoints, &points);
668: elem_list[elem_ind[fs] + i] = points[2] +1;
669: DMPlexRestoreTransitiveClosure(dm, faces[i], PETSC_FALSE, &numPoints, &points);
671: /* Side List */
672: points = NULL;
673: DMPlexGetTransitiveClosure(dm, elem_list[elem_ind[fs] + i]-1, PETSC_TRUE, &numPoints, &points);
674: for (j=1; j<numPoints; ++j) {
675: if (points[j*2]==faces[i]) {break;}
676: }
677: /* Convert HEX sides */
678: if (numPoints == 27) {
679: if (j == 1) {j = 5;}
680: else if (j == 2) {j = 6;}
681: else if (j == 3) {j = 1;}
682: else if (j == 4) {j = 3;}
683: else if (j == 5) {j = 2;}
684: else if (j == 6) {j = 4;}
685: }
686: /* Convert TET sides */
687: if (numPoints == 15) {
688: --j;
689: if (j == 0) {j = 4;}
690: }
691: side_list[elem_ind[fs] + i] = j;
692: DMPlexRestoreTransitiveClosure(dm, elem_list[elem_ind[fs] + i]-1, PETSC_TRUE, &numPoints, &points);
694: }
695: ISRestoreIndices(stratumIS, &faces);
696: ISDestroy(&stratumIS);
697: }
698: ISRestoreIndices(fsIS, &fsIdx);
699: ISDestroy(&fsIS);
701: /* Put side sets */
702: for (fs=0; fs<num_fs; ++fs) {
703: PetscStackCallStandard(ex_put_set,(exoid, EX_SIDE_SET, fsIdx[fs], &elem_list[elem_ind[fs]], &side_list[elem_ind[fs]]));
704: }
705: PetscFree3(elem_ind,elem_list,side_list);
706: }
707: return(0);
708: }
710: /*
711: VecViewPlex_ExodusII_Nodal_Internal - Write a Vec corresponding to a nodal field to an exodus file
713: Collective on v
715: Input Parameters:
716: + v - The vector to be written
717: . exoid - the exodus id of a file (obtained from ex_open or ex_create for instance)
718: - step - the time step to write at (exodus steps are numbered starting from 1)
720: Notes:
721: The exodus result nodal variable index is obtained by comparing the Vec name and the
722: names of nodal variables declared in the exodus file. For instance for a Vec named "V"
723: the location in the exodus file will be the first match of "V", "V_X", "V_XX", "V_1", or "V_11"
724: amongst all nodal variables.
726: Level: beginner
728: .seealso: EXOGetVarIndex_Private(),DMPlexView_ExodusII_Internal(),VecLoadNodal_PlexEXO(),VecViewZonal_PlexEXO(),VecLoadZonal_PlexEXO()
729: @*/
730: PetscErrorCode VecViewPlex_ExodusII_Nodal_Internal(Vec v, int exoid, int step)
731: {
732: MPI_Comm comm;
733: PetscMPIInt size;
734: DM dm;
735: Vec vNatural, vComp;
736: const PetscScalar *varray;
737: const char *vecname;
738: PetscInt xs, xe, bs;
739: PetscBool useNatural;
740: int offset;
741: PetscErrorCode ierr;
744: PetscObjectGetComm((PetscObject) v, &comm);
745: MPI_Comm_size(comm, &size);
746: VecGetDM(v, &dm);
747: DMGetUseNatural(dm, &useNatural);
748: useNatural = useNatural && size > 1 ? PETSC_TRUE : PETSC_FALSE;
749: if (useNatural) {
750: DMGetGlobalVector(dm, &vNatural);
751: DMPlexGlobalToNaturalBegin(dm, v, vNatural);
752: DMPlexGlobalToNaturalEnd(dm, v, vNatural);
753: } else {
754: vNatural = v;
755: }
756: /* Get the location of the variable in the exodus file */
757: PetscObjectGetName((PetscObject) v, &vecname);
758: EXOGetVarIndex_Private(exoid, EX_NODAL, vecname, &offset);
759: if (!offset) SETERRQ1(comm, PETSC_ERR_FILE_UNEXPECTED, "Cannot locate nodal variable %s in exodus file\n", vecname);
760: /* Write local chunk of the result in the exodus file
761: exodus stores each component of a vector-valued field as a separate variable.
762: We assume that they are stored sequentially */
763: VecGetOwnershipRange(vNatural, &xs, &xe);
764: VecGetBlockSize(vNatural, &bs);
765: if (bs == 1) {
766: VecGetArrayRead(vNatural, &varray);
767: PetscStackCallStandard(ex_put_partial_var,(exoid, step, EX_NODAL, offset, 1, xs+1, xe-xs, varray));
768: VecRestoreArrayRead(vNatural, &varray);
769: } else {
770: IS compIS;
771: PetscInt c;
773: ISCreateStride(PETSC_COMM_SELF, (xe-xs)/bs, xs, bs, &compIS);
774: for (c = 0; c < bs; ++c) {
775: ISStrideSetStride(compIS, (xe-xs)/bs, xs+c, bs);
776: VecGetSubVector(vNatural, compIS, &vComp);
777: VecGetArrayRead(vComp, &varray);
778: PetscStackCallStandard(ex_put_partial_var,(exoid, step, EX_NODAL, offset+c, 1, xs/bs+1, (xe-xs)/bs, varray));
779: VecRestoreArrayRead(vComp, &varray);
780: VecRestoreSubVector(vNatural, compIS, &vComp);
781: }
782: ISDestroy(&compIS);
783: }
784: if (useNatural) {DMRestoreGlobalVector(dm, &vNatural);}
785: return(0);
786: }
788: /*
789: VecLoadPlex_ExodusII_Nodal_Internal - Read a Vec corresponding to a nodal field from an exodus file
791: Collective on v
793: Input Parameters:
794: + v - The vector to be written
795: . exoid - the exodus id of a file (obtained from ex_open or ex_create for instance)
796: - step - the time step to read at (exodus steps are numbered starting from 1)
798: Notes:
799: The exodus result nodal variable index is obtained by comparing the Vec name and the
800: names of nodal variables declared in the exodus file. For instance for a Vec named "V"
801: the location in the exodus file will be the first match of "V", "V_X", "V_XX", "V_1", or "V_11"
802: amongst all nodal variables.
804: Level: beginner
806: .seealso: EXOGetVarIndex_Private(), DMPlexView_ExodusII_Internal(), VecViewPlex_ExodusII_Nodal_Internal(), VecViewZonal_PlexEXO(), VecLoadZonal_PlexEXO()
807: */
808: PetscErrorCode VecLoadPlex_ExodusII_Nodal_Internal(Vec v, int exoid, int step)
809: {
810: MPI_Comm comm;
811: PetscMPIInt size;
812: DM dm;
813: Vec vNatural, vComp;
814: PetscScalar *varray;
815: const char *vecname;
816: PetscInt xs, xe, bs;
817: PetscBool useNatural;
818: int offset;
822: PetscObjectGetComm((PetscObject) v, &comm);
823: MPI_Comm_size(comm, &size);
824: VecGetDM(v,&dm);
825: DMGetUseNatural(dm, &useNatural);
826: useNatural = useNatural && size > 1 ? PETSC_TRUE : PETSC_FALSE;
827: if (useNatural) {DMGetGlobalVector(dm,&vNatural);}
828: else {vNatural = v;}
829: /* Get the location of the variable in the exodus file */
830: PetscObjectGetName((PetscObject) v, &vecname);
831: EXOGetVarIndex_Private(exoid, EX_NODAL, vecname, &offset);
832: if (!offset) SETERRQ1(comm, PETSC_ERR_FILE_UNEXPECTED, "Cannot locate nodal variable %s in exodus file", vecname);
833: /* Read local chunk from the file */
834: VecGetOwnershipRange(vNatural, &xs, &xe);
835: VecGetBlockSize(vNatural, &bs);
836: if (bs == 1) {
837: VecGetArray(vNatural, &varray);
838: PetscStackCallStandard(ex_get_partial_var,(exoid, step, EX_NODAL, offset, 1, xs+1, xe-xs, varray));
839: VecRestoreArray(vNatural, &varray);
840: } else {
841: IS compIS;
842: PetscInt c;
844: ISCreateStride(PETSC_COMM_SELF, (xe-xs)/bs, xs, bs, &compIS);
845: for (c = 0; c < bs; ++c) {
846: ISStrideSetStride(compIS, (xe-xs)/bs, xs+c, bs);
847: VecGetSubVector(vNatural, compIS, &vComp);
848: VecGetArray(vComp, &varray);
849: PetscStackCallStandard(ex_get_partial_var,(exoid, step, EX_NODAL, offset+c, 1, xs/bs+1, (xe-xs)/bs, varray));
850: VecRestoreArray(vComp, &varray);
851: VecRestoreSubVector(vNatural, compIS, &vComp);
852: }
853: ISDestroy(&compIS);
854: }
855: if (useNatural) {
856: DMPlexNaturalToGlobalBegin(dm, vNatural, v);
857: DMPlexNaturalToGlobalEnd(dm, vNatural, v);
858: DMRestoreGlobalVector(dm, &vNatural);
859: }
860: return(0);
861: }
863: /*
864: VecViewPlex_ExodusII_Zonal_Internal - Write a Vec corresponding to a zonal (cell based) field to an exodus file
866: Collective on v
868: Input Parameters:
869: + v - The vector to be written
870: . exoid - the exodus id of a file (obtained from ex_open or ex_create for instance)
871: - step - the time step to write at (exodus steps are numbered starting from 1)
873: Notes:
874: The exodus result zonal variable index is obtained by comparing the Vec name and the
875: names of zonal variables declared in the exodus file. For instance for a Vec named "V"
876: the location in the exodus file will be the first match of "V", "V_X", "V_XX", "V_1", or "V_11"
877: amongst all zonal variables.
879: Level: beginner
881: .seealso: EXOGetVarIndex_Private(),DMPlexView_ExodusII_Internal(),VecViewPlex_ExodusII_Nodal_Internal(),VecLoadPlex_ExodusII_Nodal_Internal(),VecLoadPlex_ExodusII_Zonal_Internal()
882: */
883: PetscErrorCode VecViewPlex_ExodusII_Zonal_Internal(Vec v, int exoid, int step)
884: {
885: MPI_Comm comm;
886: PetscMPIInt size;
887: DM dm;
888: Vec vNatural, vComp;
889: const PetscScalar *varray;
890: const char *vecname;
891: PetscInt xs, xe, bs;
892: PetscBool useNatural;
893: int offset;
894: IS compIS;
895: PetscInt *csSize, *csID;
896: PetscInt numCS, set, csxs = 0;
897: PetscErrorCode ierr;
900: PetscObjectGetComm((PetscObject)v, &comm);
901: MPI_Comm_size(comm, &size);
902: VecGetDM(v, &dm);
903: DMGetUseNatural(dm, &useNatural);
904: useNatural = useNatural && size > 1 ? PETSC_TRUE : PETSC_FALSE;
905: if (useNatural) {
906: DMGetGlobalVector(dm, &vNatural);
907: DMPlexGlobalToNaturalBegin(dm, v, vNatural);
908: DMPlexGlobalToNaturalEnd(dm, v, vNatural);
909: } else {
910: vNatural = v;
911: }
912: /* Get the location of the variable in the exodus file */
913: PetscObjectGetName((PetscObject) v, &vecname);
914: EXOGetVarIndex_Private(exoid, EX_ELEM_BLOCK, vecname, &offset);
915: if (!offset) SETERRQ1(comm, PETSC_ERR_FILE_UNEXPECTED, "Cannot locate zonal variable %s in exodus file\n", vecname);
916: /* Write local chunk of the result in the exodus file
917: exodus stores each component of a vector-valued field as a separate variable.
918: We assume that they are stored sequentially
919: Zonal variables are accessed one element block at a time, so we loop through the cell sets,
920: but once the vector has been reordered to natural size, we cannot use the label informations
921: to figure out what to save where. */
922: numCS = ex_inquire_int(exoid, EX_INQ_ELEM_BLK);
923: PetscMalloc2(numCS, &csID, numCS, &csSize);
924: PetscStackCallStandard(ex_get_ids,(exoid, EX_ELEM_BLOCK, csID));
925: for (set = 0; set < numCS; ++set) {
926: ex_block block;
928: block.id = csID[set];
929: block.type = EX_ELEM_BLOCK;
930: PetscStackCallStandard(ex_get_block_param,(exoid, &block));
931: csSize[set] = block.num_entry;
932: }
933: VecGetOwnershipRange(vNatural, &xs, &xe);
934: VecGetBlockSize(vNatural, &bs);
935: if (bs > 1) {ISCreateStride(comm, (xe-xs)/bs, xs, bs, &compIS);}
936: for (set = 0; set < numCS; set++) {
937: PetscInt csLocalSize, c;
939: /* range of indices for set setID[set]: csxs:csxs + csSize[set]-1
940: local slice of zonal values: xs/bs,xm/bs-1
941: intersection: max(xs/bs,csxs),min(xm/bs-1,csxs + csSize[set]-1) */
942: csLocalSize = PetscMax(0, PetscMin(xe/bs, csxs+csSize[set]) - PetscMax(xs/bs, csxs));
943: if (bs == 1) {
944: VecGetArrayRead(vNatural, &varray);
945: PetscStackCallStandard(ex_put_partial_var,(exoid, step, EX_ELEM_BLOCK, offset, csID[set], PetscMax(xs-csxs, 0)+1, csLocalSize, &varray[PetscMax(0, csxs-xs)]));
946: VecRestoreArrayRead(vNatural, &varray);
947: } else {
948: for (c = 0; c < bs; ++c) {
949: ISStrideSetStride(compIS, (xe-xs)/bs, xs+c, bs);
950: VecGetSubVector(vNatural, compIS, &vComp);
951: VecGetArrayRead(vComp, &varray);
952: PetscStackCallStandard(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)]));
953: VecRestoreArrayRead(vComp, &varray);
954: VecRestoreSubVector(vNatural, compIS, &vComp);
955: }
956: }
957: csxs += csSize[set];
958: }
959: PetscFree2(csID, csSize);
960: if (bs > 1) {ISDestroy(&compIS);}
961: if (useNatural) {DMRestoreGlobalVector(dm,&vNatural);}
962: return(0);
963: }
965: /*
966: VecLoadPlex_ExodusII_Zonal_Internal - Read a Vec corresponding to a zonal (cell based) field from an exodus file
968: Collective on v
970: Input Parameters:
971: + v - The vector to be written
972: . exoid - the exodus id of a file (obtained from ex_open or ex_create for instance)
973: - step - the time step to read at (exodus steps are numbered starting from 1)
975: Notes:
976: The exodus result zonal variable index is obtained by comparing the Vec name and the
977: names of zonal variables declared in the exodus file. For instance for a Vec named "V"
978: the location in the exodus file will be the first match of "V", "V_X", "V_XX", "V_1", or "V_11"
979: amongst all zonal variables.
981: Level: beginner
983: .seealso: EXOGetVarIndex_Private(), DMPlexView_ExodusII_Internal(), VecViewPlex_ExodusII_Nodal_Internal(), VecLoadPlex_ExodusII_Nodal_Internal(), VecLoadPlex_ExodusII_Zonal_Internal()
984: */
985: PetscErrorCode VecLoadPlex_ExodusII_Zonal_Internal(Vec v, int exoid, int step)
986: {
987: MPI_Comm comm;
988: PetscMPIInt size;
989: DM dm;
990: Vec vNatural, vComp;
991: PetscScalar *varray;
992: const char *vecname;
993: PetscInt xs, xe, bs;
994: PetscBool useNatural;
995: int offset;
996: IS compIS;
997: PetscInt *csSize, *csID;
998: PetscInt numCS, set, csxs = 0;
999: PetscErrorCode ierr;
1002: PetscObjectGetComm((PetscObject)v,&comm);
1003: MPI_Comm_size(comm, &size);
1004: VecGetDM(v, &dm);
1005: DMGetUseNatural(dm, &useNatural);
1006: useNatural = useNatural && size > 1 ? PETSC_TRUE : PETSC_FALSE;
1007: if (useNatural) {DMGetGlobalVector(dm,&vNatural);}
1008: else {vNatural = v;}
1009: /* Get the location of the variable in the exodus file */
1010: PetscObjectGetName((PetscObject) v, &vecname);
1011: EXOGetVarIndex_Private(exoid, EX_ELEM_BLOCK, vecname, &offset);
1012: if (!offset) SETERRQ1(comm, PETSC_ERR_FILE_UNEXPECTED, "Cannot locate zonal variable %s in exodus file\n", vecname);
1013: /* Read local chunk of the result in the exodus file
1014: exodus stores each component of a vector-valued field as a separate variable.
1015: We assume that they are stored sequentially
1016: Zonal variables are accessed one element block at a time, so we loop through the cell sets,
1017: but once the vector has been reordered to natural size, we cannot use the label informations
1018: to figure out what to save where. */
1019: numCS = ex_inquire_int(exoid, EX_INQ_ELEM_BLK);
1020: PetscMalloc2(numCS, &csID, numCS, &csSize);
1021: PetscStackCallStandard(ex_get_ids,(exoid, EX_ELEM_BLOCK, csID));
1022: for (set = 0; set < numCS; ++set) {
1023: ex_block block;
1025: block.id = csID[set];
1026: block.type = EX_ELEM_BLOCK;
1027: PetscStackCallStandard(ex_get_block_param,(exoid, &block));
1028: csSize[set] = block.num_entry;
1029: }
1030: VecGetOwnershipRange(vNatural, &xs, &xe);
1031: VecGetBlockSize(vNatural, &bs);
1032: if (bs > 1) {ISCreateStride(comm, (xe-xs)/bs, xs, bs, &compIS);}
1033: for (set = 0; set < numCS; ++set) {
1034: PetscInt csLocalSize, c;
1036: /* range of indices for set setID[set]: csxs:csxs + csSize[set]-1
1037: local slice of zonal values: xs/bs,xm/bs-1
1038: intersection: max(xs/bs,csxs),min(xm/bs-1,csxs + csSize[set]-1) */
1039: csLocalSize = PetscMax(0, PetscMin(xe/bs, csxs+csSize[set]) - PetscMax(xs/bs, csxs));
1040: if (bs == 1) {
1041: VecGetArray(vNatural, &varray);
1042: PetscStackCallStandard(ex_get_partial_var,(exoid, step, EX_ELEM_BLOCK, offset, csID[set], PetscMax(xs-csxs, 0)+1, csLocalSize, &varray[PetscMax(0, csxs-xs)]));
1043: VecRestoreArray(vNatural, &varray);
1044: } else {
1045: for (c = 0; c < bs; ++c) {
1046: ISStrideSetStride(compIS, (xe-xs)/bs, xs+c, bs);
1047: VecGetSubVector(vNatural, compIS, &vComp);
1048: VecGetArray(vComp, &varray);
1049: PetscStackCallStandard(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)]));
1050: VecRestoreArray(vComp, &varray);
1051: VecRestoreSubVector(vNatural, compIS, &vComp);
1052: }
1053: }
1054: csxs += csSize[set];
1055: }
1056: PetscFree2(csID, csSize);
1057: if (bs > 1) {ISDestroy(&compIS);}
1058: if (useNatural) {
1059: DMPlexNaturalToGlobalBegin(dm, vNatural, v);
1060: DMPlexNaturalToGlobalEnd(dm, vNatural, v);
1061: DMRestoreGlobalVector(dm, &vNatural);
1062: }
1063: return(0);
1064: }
1065: #endif
1067: /*@
1068: PetscViewerExodusIIGetId - Get the file id of the ExodusII file
1070: Logically Collective on PetscViewer
1072: Input Parameter:
1073: . viewer - the PetscViewer
1075: Output Parameter:
1076: - exoid - The ExodusII file id
1078: Level: intermediate
1080: .seealso: PetscViewerFileSetMode(), PetscViewerCreate(), PetscViewerSetType(), PetscViewerBinaryOpen()
1081: @*/
1082: PetscErrorCode PetscViewerExodusIIGetId(PetscViewer viewer, int *exoid)
1083: {
1088: PetscTryMethod(viewer, "PetscViewerExodusIIGetId_C",(PetscViewer,int*),(viewer,exoid));
1089: return(0);
1090: }
1092: /*@C
1093: DMPlexCreateExodusFromFile - Create a DMPlex mesh from an ExodusII file.
1095: Collective
1097: Input Parameters:
1098: + comm - The MPI communicator
1099: . filename - The name of the ExodusII file
1100: - interpolate - Create faces and edges in the mesh
1102: Output Parameter:
1103: . dm - The DM object representing the mesh
1105: Level: beginner
1107: .seealso: DMPLEX, DMCreate(), DMPlexCreateExodus()
1108: @*/
1109: PetscErrorCode DMPlexCreateExodusFromFile(MPI_Comm comm, const char filename[], PetscBool interpolate, DM *dm)
1110: {
1111: PetscMPIInt rank;
1113: #if defined(PETSC_HAVE_EXODUSII)
1114: int CPU_word_size = sizeof(PetscReal), IO_word_size = 0, exoid = -1;
1115: float version;
1116: #endif
1120: MPI_Comm_rank(comm, &rank);
1121: #if defined(PETSC_HAVE_EXODUSII)
1122: if (!rank) {
1123: exoid = ex_open(filename, EX_READ, &CPU_word_size, &IO_word_size, &version);
1124: if (exoid <= 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_LIB, "ex_open(\"%s\",...) did not return a valid file ID", filename);
1125: }
1126: DMPlexCreateExodus(comm, exoid, interpolate, dm);
1127: if (!rank) {PetscStackCallStandard(ex_close,(exoid));}
1128: #else
1129: SETERRQ(comm, PETSC_ERR_SUP, "This method requires ExodusII support. Reconfigure using --download-exodusii");
1130: #endif
1131: return(0);
1132: }
1134: /*@
1135: DMPlexCreateExodus - Create a DMPlex mesh from an ExodusII file ID.
1137: Collective
1139: Input Parameters:
1140: + comm - The MPI communicator
1141: . exoid - The ExodusII id associated with a exodus file and obtained using ex_open
1142: - interpolate - Create faces and edges in the mesh
1144: Output Parameter:
1145: . dm - The DM object representing the mesh
1147: Level: beginner
1149: .seealso: DMPLEX, DMCreate()
1150: @*/
1151: PetscErrorCode DMPlexCreateExodus(MPI_Comm comm, PetscInt exoid, PetscBool interpolate, DM *dm)
1152: {
1153: #if defined(PETSC_HAVE_EXODUSII)
1154: PetscMPIInt num_proc, rank;
1155: PetscSection coordSection;
1156: Vec coordinates;
1157: PetscScalar *coords;
1158: PetscInt coordSize, v;
1160: /* Read from ex_get_init() */
1161: char title[PETSC_MAX_PATH_LEN+1];
1162: int dim = 0, numVertices = 0, numCells = 0, numHybridCells = 0;
1163: int num_cs = 0, num_vs = 0, num_fs = 0;
1164: #endif
1167: #if defined(PETSC_HAVE_EXODUSII)
1168: MPI_Comm_rank(comm, &rank);
1169: MPI_Comm_size(comm, &num_proc);
1170: DMCreate(comm, dm);
1171: DMSetType(*dm, DMPLEX);
1172: /* Open EXODUS II file and read basic informations on rank 0, then broadcast to all processors */
1173: if (!rank) {
1174: PetscMemzero(title,PETSC_MAX_PATH_LEN+1);
1175: PetscStackCallStandard(ex_get_init,(exoid, title, &dim, &numVertices, &numCells, &num_cs, &num_vs, &num_fs));
1176: if (!num_cs) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Exodus file does not contain any cell set\n");
1177: }
1178: MPI_Bcast(title, PETSC_MAX_PATH_LEN+1, MPI_CHAR, 0, comm);
1179: MPI_Bcast(&dim, 1, MPI_INT, 0, comm);
1180: PetscObjectSetName((PetscObject) *dm, title);
1181: DMSetDimension(*dm, dim);
1182: DMPlexSetChart(*dm, 0, numCells+numVertices);
1183: /* We do not want this label automatically computed, instead we compute it here */
1184: DMCreateLabel(*dm, "celltype");
1186: /* Read cell sets information */
1187: if (!rank) {
1188: PetscInt *cone;
1189: int c, cs, ncs, c_loc, v, v_loc;
1190: /* Read from ex_get_elem_blk_ids() */
1191: int *cs_id, *cs_order;
1192: /* Read from ex_get_elem_block() */
1193: char buffer[PETSC_MAX_PATH_LEN+1];
1194: int num_cell_in_set, num_vertex_per_cell, num_hybrid, num_attr;
1195: /* Read from ex_get_elem_conn() */
1196: int *cs_connect;
1198: /* Get cell sets IDs */
1199: PetscMalloc2(num_cs, &cs_id, num_cs, &cs_order);
1200: PetscStackCallStandard(ex_get_ids,(exoid, EX_ELEM_BLOCK, cs_id));
1201: /* Read the cell set connectivity table and build mesh topology
1202: EXO standard requires that cells in cell sets be numbered sequentially and be pairwise disjoint. */
1203: /* Check for a hybrid mesh */
1204: for (cs = 0, num_hybrid = 0; cs < num_cs; ++cs) {
1205: PetscStackCallStandard(ex_get_block,(exoid, EX_ELEM_BLOCK, cs_id[cs], buffer, &num_cell_in_set,&num_vertex_per_cell, 0, 0, &num_attr));
1206: switch (dim) {
1207: case 3:
1208: switch (num_vertex_per_cell) {
1209: case 6:
1210: cs_order[cs] = cs;
1211: numHybridCells += num_cell_in_set;
1212: ++num_hybrid;
1213: break;
1214: default:
1215: for (c = cs; c > cs-num_hybrid; --c) cs_order[c] = cs_order[c-1];
1216: cs_order[cs-num_hybrid] = cs;
1217: }
1218: break;
1219: default:
1220: for (c = cs; c > cs-num_hybrid; --c) cs_order[c] = cs_order[c-1];
1221: cs_order[cs-num_hybrid] = cs;
1222: }
1223: }
1224: /* First set sizes */
1225: for (ncs = 0, c = 0; ncs < num_cs; ++ncs) {
1226: const PetscInt cs = cs_order[ncs];
1227: PetscStackCallStandard(ex_get_block,(exoid, EX_ELEM_BLOCK, cs_id[cs], buffer, &num_cell_in_set,&num_vertex_per_cell, 0, 0, &num_attr));
1228: for (c_loc = 0; c_loc < num_cell_in_set; ++c_loc, ++c) {
1229: DMPlexSetConeSize(*dm, c, num_vertex_per_cell);
1230: if (c >= numCells-numHybridCells) {
1231: DMPlexSetCellType(*dm, c, DM_POLYTOPE_TRI_PRISM);
1232: } else {
1233: DMPolytopeType ct;
1235: DMPlexComputeCellType_Internal(*dm, c, 1, &ct);
1236: DMPlexSetCellType(*dm, c, ct);
1237: }
1238: }
1239: }
1240: for (v = numCells; v < numCells+numVertices; ++v) {DMPlexSetCellType(*dm, v, DM_POLYTOPE_POINT);}
1241: DMSetUp(*dm);
1242: for (ncs = 0, c = 0; ncs < num_cs; ++ncs) {
1243: const PetscInt cs = cs_order[ncs];
1244: PetscStackCallStandard(ex_get_block,(exoid, EX_ELEM_BLOCK, cs_id[cs], buffer, &num_cell_in_set, &num_vertex_per_cell, 0, 0, &num_attr));
1245: PetscMalloc2(num_vertex_per_cell*num_cell_in_set,&cs_connect,num_vertex_per_cell,&cone);
1246: PetscStackCallStandard(ex_get_conn,(exoid, EX_ELEM_BLOCK, cs_id[cs], cs_connect,NULL,NULL));
1247: /* EXO uses Fortran-based indexing, DMPlex uses C-style and numbers cell first then vertices. */
1248: for (c_loc = 0, v = 0; c_loc < num_cell_in_set; ++c_loc, ++c) {
1249: for (v_loc = 0; v_loc < num_vertex_per_cell; ++v_loc, ++v) {
1250: cone[v_loc] = cs_connect[v]+numCells-1;
1251: }
1252: if (dim == 3) {
1253: /* Tetrahedra are inverted */
1254: if (num_vertex_per_cell == 4) {
1255: PetscInt tmp = cone[0];
1256: cone[0] = cone[1];
1257: cone[1] = tmp;
1258: }
1259: /* Hexahedra are inverted */
1260: if (num_vertex_per_cell == 8) {
1261: PetscInt tmp = cone[1];
1262: cone[1] = cone[3];
1263: cone[3] = tmp;
1264: }
1265: /* Triangular prisms are inverted */
1266: if (num_vertex_per_cell == 6) {
1267: PetscInt tmp = cone[1];
1268: cone[1] = cone[2];
1269: cone[2] = tmp;
1270: }
1271: }
1272: DMPlexSetCone(*dm, c, cone);
1273: DMSetLabelValue(*dm, "Cell Sets", c, cs_id[cs]);
1274: }
1275: PetscFree2(cs_connect,cone);
1276: }
1277: PetscFree2(cs_id, cs_order);
1278: }
1279: DMPlexSymmetrize(*dm);
1280: DMPlexStratify(*dm);
1281: if (interpolate) {
1282: DM idm;
1284: DMPlexInterpolate(*dm, &idm);
1285: DMDestroy(dm);
1286: *dm = idm;
1287: }
1289: /* Create vertex set label */
1290: if (!rank && (num_vs > 0)) {
1291: int vs, v;
1292: /* Read from ex_get_node_set_ids() */
1293: int *vs_id;
1294: /* Read from ex_get_node_set_param() */
1295: int num_vertex_in_set;
1296: /* Read from ex_get_node_set() */
1297: int *vs_vertex_list;
1299: /* Get vertex set ids */
1300: PetscMalloc1(num_vs, &vs_id);
1301: PetscStackCallStandard(ex_get_ids,(exoid, EX_NODE_SET, vs_id));
1302: for (vs = 0; vs < num_vs; ++vs) {
1303: PetscStackCallStandard(ex_get_set_param,(exoid, EX_NODE_SET, vs_id[vs], &num_vertex_in_set, NULL));
1304: PetscMalloc1(num_vertex_in_set, &vs_vertex_list);
1305: PetscStackCallStandard(ex_get_set,(exoid, EX_NODE_SET, vs_id[vs], vs_vertex_list, NULL));
1306: for (v = 0; v < num_vertex_in_set; ++v) {
1307: DMSetLabelValue(*dm, "Vertex Sets", vs_vertex_list[v]+numCells-1, vs_id[vs]);
1308: }
1309: PetscFree(vs_vertex_list);
1310: }
1311: PetscFree(vs_id);
1312: }
1313: /* Read coordinates */
1314: DMGetCoordinateSection(*dm, &coordSection);
1315: PetscSectionSetNumFields(coordSection, 1);
1316: PetscSectionSetFieldComponents(coordSection, 0, dim);
1317: PetscSectionSetChart(coordSection, numCells, numCells + numVertices);
1318: for (v = numCells; v < numCells+numVertices; ++v) {
1319: PetscSectionSetDof(coordSection, v, dim);
1320: PetscSectionSetFieldDof(coordSection, v, 0, dim);
1321: }
1322: PetscSectionSetUp(coordSection);
1323: PetscSectionGetStorageSize(coordSection, &coordSize);
1324: VecCreate(PETSC_COMM_SELF, &coordinates);
1325: PetscObjectSetName((PetscObject) coordinates, "coordinates");
1326: VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);
1327: VecSetBlockSize(coordinates, dim);
1328: VecSetType(coordinates,VECSTANDARD);
1329: VecGetArray(coordinates, &coords);
1330: if (!rank) {
1331: PetscReal *x, *y, *z;
1333: PetscMalloc3(numVertices,&x,numVertices,&y,numVertices,&z);
1334: PetscStackCallStandard(ex_get_coord,(exoid, x, y, z));
1335: if (dim > 0) {
1336: for (v = 0; v < numVertices; ++v) coords[v*dim+0] = x[v];
1337: }
1338: if (dim > 1) {
1339: for (v = 0; v < numVertices; ++v) coords[v*dim+1] = y[v];
1340: }
1341: if (dim > 2) {
1342: for (v = 0; v < numVertices; ++v) coords[v*dim+2] = z[v];
1343: }
1344: PetscFree3(x,y,z);
1345: }
1346: VecRestoreArray(coordinates, &coords);
1347: DMSetCoordinatesLocal(*dm, coordinates);
1348: VecDestroy(&coordinates);
1350: /* Create side set label */
1351: if (!rank && interpolate && (num_fs > 0)) {
1352: int fs, f, voff;
1353: /* Read from ex_get_side_set_ids() */
1354: int *fs_id;
1355: /* Read from ex_get_side_set_param() */
1356: int num_side_in_set;
1357: /* Read from ex_get_side_set_node_list() */
1358: int *fs_vertex_count_list, *fs_vertex_list;
1359: /* Read side set labels */
1360: char fs_name[MAX_STR_LENGTH+1];
1362: /* Get side set ids */
1363: PetscMalloc1(num_fs, &fs_id);
1364: PetscStackCallStandard(ex_get_ids,(exoid, EX_SIDE_SET, fs_id));
1365: for (fs = 0; fs < num_fs; ++fs) {
1366: PetscStackCallStandard(ex_get_set_param,(exoid, EX_SIDE_SET, fs_id[fs], &num_side_in_set, NULL));
1367: PetscMalloc2(num_side_in_set,&fs_vertex_count_list,num_side_in_set*4,&fs_vertex_list);
1368: PetscStackCallStandard(ex_get_side_set_node_list,(exoid, fs_id[fs], fs_vertex_count_list, fs_vertex_list));
1369: /* Get the specific name associated with this side set ID. */
1370: int fs_name_err = ex_get_name(exoid, EX_SIDE_SET, fs_id[fs], fs_name);
1371: for (f = 0, voff = 0; f < num_side_in_set; ++f) {
1372: const PetscInt *faces = NULL;
1373: PetscInt faceSize = fs_vertex_count_list[f], numFaces;
1374: PetscInt faceVertices[4], v;
1376: if (faceSize > 4) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "ExodusII side cannot have %d > 4 vertices", faceSize);
1377: for (v = 0; v < faceSize; ++v, ++voff) {
1378: faceVertices[v] = fs_vertex_list[voff]+numCells-1;
1379: }
1380: DMPlexGetFullJoin(*dm, faceSize, faceVertices, &numFaces, &faces);
1381: if (numFaces != 1) SETERRQ3(comm, PETSC_ERR_ARG_WRONG, "Invalid ExodusII side %d in set %d maps to %d faces", f, fs, numFaces);
1382: DMSetLabelValue(*dm, "Face Sets", faces[0], fs_id[fs]);
1383: /* Only add the label if one has been detected for this side set. */
1384: if (!fs_name_err) {
1385: DMSetLabelValue(*dm, fs_name, faces[0], fs_id[fs]);
1386: }
1387: DMPlexRestoreJoin(*dm, faceSize, faceVertices, &numFaces, &faces);
1388: }
1389: PetscFree2(fs_vertex_count_list,fs_vertex_list);
1390: }
1391: PetscFree(fs_id);
1392: }
1393: #else
1394: SETERRQ(comm, PETSC_ERR_SUP, "This method requires ExodusII support. Reconfigure using --download-exodusii");
1395: #endif
1396: return(0);
1397: }