Actual source code: plexexodusii.c
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;
37: PetscViewerExodusIIOpen(comm, "mesh.exo", FILE_MODE_WRITE, &viewer);
38: if (ierr) {PetscError(PETSC_COMM_SELF,__LINE__,"PETSC_VIEWER_EXODUSII_",__FILE__,PETSC_ERR_PLIB,PETSC_ERROR_INITIAL," ");return 0;}
39: PetscObjectRegisterDestroy((PetscObject) viewer);
40: if (ierr) {PetscError(PETSC_COMM_SELF,__LINE__,"PETSC_VIEWER_EXODUSII_",__FILE__,PETSC_ERR_PLIB,PETSC_ERROR_INITIAL," ");return 0;}
41: return viewer;
42: }
44: static PetscErrorCode PetscViewerView_ExodusII(PetscViewer v, PetscViewer viewer)
45: {
46: PetscViewer_ExodusII *exo = (PetscViewer_ExodusII *) v->data;
48: if (exo->filename) PetscViewerASCIIPrintf(viewer, "Filename: %s\n", exo->filename);
49: if (exo->exoid) PetscViewerASCIIPrintf(viewer, "exoid: %d\n", exo->exoid);
50: if (exo->btype) PetscViewerASCIIPrintf(viewer, "IO Mode: %d\n", exo->btype);
51: if (exo->order) PetscViewerASCIIPrintf(viewer, "Mesh order: %d\n", exo->order);
52: return 0;
53: }
55: static PetscErrorCode PetscViewerSetFromOptions_ExodusII(PetscOptionItems *PetscOptionsObject, PetscViewer v)
56: {
57: PetscOptionsHead(PetscOptionsObject, "ExodusII PetscViewer Options");
58: PetscOptionsTail();
59: return 0;
60: }
62: static PetscErrorCode PetscViewerSetUp_ExodusII(PetscViewer viewer)
63: {
64: return 0;
65: }
67: static PetscErrorCode PetscViewerDestroy_ExodusII(PetscViewer viewer)
68: {
69: PetscViewer_ExodusII *exo = (PetscViewer_ExodusII *) viewer->data;
71: if (exo->exoid >= 0) {PetscStackCallStandard(ex_close,exo->exoid);}
72: PetscFree(exo->filename);
73: PetscFree(exo);
74: PetscObjectComposeFunction((PetscObject)viewer,"PetscViewerFileSetName_C",NULL);
75: PetscObjectComposeFunction((PetscObject)viewer,"PetscViewerFileGetName_C",NULL);
76: PetscObjectComposeFunction((PetscObject)viewer,"PetscViewerFileSetMode_C",NULL);
77: PetscObjectComposeFunction((PetscObject)viewer,"PetscViewerGetId_C",NULL);
78: PetscObjectComposeFunction((PetscObject)viewer,"PetscViewerGetOrder_C",NULL);
79: PetscObjectComposeFunction((PetscObject)viewer,"PetscViewerSetOrder_C",NULL);
80: return 0;
81: }
83: static PetscErrorCode PetscViewerFileSetName_ExodusII(PetscViewer viewer, const char name[])
84: {
85: PetscViewer_ExodusII *exo = (PetscViewer_ExodusII *) viewer->data;
86: PetscMPIInt rank;
87: int CPU_word_size, IO_word_size, EXO_mode;
88: MPI_Info mpi_info = MPI_INFO_NULL;
89: float EXO_version;
91: MPI_Comm_rank(PetscObjectComm((PetscObject) viewer), &rank);
92: CPU_word_size = sizeof(PetscReal);
93: IO_word_size = sizeof(PetscReal);
95: if (exo->exoid >= 0) {
96: PetscStackCallStandard(ex_close,exo->exoid);
97: exo->exoid = -1;
98: }
99: if (exo->filename) PetscFree(exo->filename);
100: PetscStrallocpy(name, &exo->filename);
101: switch (exo->btype) {
102: case FILE_MODE_READ:
103: EXO_mode = EX_READ;
104: break;
105: case FILE_MODE_APPEND:
106: case FILE_MODE_UPDATE:
107: case FILE_MODE_APPEND_UPDATE:
108: /* Will fail if the file does not already exist */
109: EXO_mode = EX_WRITE;
110: break;
111: case FILE_MODE_WRITE:
112: /*
113: exodus only allows writing geometry upon file creation, so we will let DMView create the file.
114: */
115: return 0;
116: break;
117: default:
118: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER, "Must call PetscViewerFileSetMode() before PetscViewerFileSetName()");
119: }
120: #if defined(PETSC_USE_64BIT_INDICES)
121: EXO_mode += EX_ALL_INT64_API;
122: #endif
123: exo->exoid = ex_open_par(name,EXO_mode,&CPU_word_size,&IO_word_size,&EXO_version,PETSC_COMM_WORLD,mpi_info);
125: return 0;
126: }
128: static PetscErrorCode PetscViewerFileGetName_ExodusII(PetscViewer viewer, const char **name)
129: {
130: PetscViewer_ExodusII *exo = (PetscViewer_ExodusII *) viewer->data;
132: *name = exo->filename;
133: return 0;
134: }
136: static PetscErrorCode PetscViewerFileSetMode_ExodusII(PetscViewer viewer, PetscFileMode type)
137: {
138: PetscViewer_ExodusII *exo = (PetscViewer_ExodusII *) viewer->data;
140: exo->btype = type;
141: return 0;
142: }
144: static PetscErrorCode PetscViewerFileGetMode_ExodusII(PetscViewer viewer, PetscFileMode *type)
145: {
146: PetscViewer_ExodusII *exo = (PetscViewer_ExodusII *) viewer->data;
148: *type = exo->btype;
149: return 0;
150: }
152: static PetscErrorCode PetscViewerExodusIIGetId_ExodusII(PetscViewer viewer, int *exoid)
153: {
154: PetscViewer_ExodusII *exo = (PetscViewer_ExodusII *) viewer->data;
156: *exoid = exo->exoid;
157: return 0;
158: }
160: static PetscErrorCode PetscViewerExodusIIGetOrder_ExodusII(PetscViewer viewer, PetscInt *order)
161: {
162: PetscViewer_ExodusII *exo = (PetscViewer_ExodusII *) viewer->data;
164: *order = exo->order;
165: return 0;
166: }
168: static PetscErrorCode PetscViewerExodusIISetOrder_ExodusII(PetscViewer viewer, PetscInt order)
169: {
170: PetscViewer_ExodusII *exo = (PetscViewer_ExodusII *) viewer->data;
172: exo->order = order;
173: return 0;
174: }
176: /*MC
177: PETSCVIEWEREXODUSII - A viewer that writes to an Exodus II file
179: .seealso: PetscViewerExodusIIOpen(), PetscViewerCreate(), PETSCVIEWERBINARY, PETSCVIEWERHDF5, DMView(),
180: PetscViewerFileSetName(), PetscViewerFileSetMode(), PetscViewerFormat, PetscViewerType, PetscViewerSetType()
182: Level: beginner
183: M*/
185: PETSC_EXTERN PetscErrorCode PetscViewerCreate_ExodusII(PetscViewer v)
186: {
187: PetscViewer_ExodusII *exo;
189: PetscNewLog(v,&exo);
191: v->data = (void*) exo;
192: v->ops->destroy = PetscViewerDestroy_ExodusII;
193: v->ops->setfromoptions = PetscViewerSetFromOptions_ExodusII;
194: v->ops->setup = PetscViewerSetUp_ExodusII;
195: v->ops->view = PetscViewerView_ExodusII;
196: v->ops->flush = 0;
197: exo->btype = (PetscFileMode) -1;
198: exo->filename = 0;
199: exo->exoid = -1;
201: PetscObjectComposeFunction((PetscObject)v,"PetscViewerFileSetName_C",PetscViewerFileSetName_ExodusII);
202: PetscObjectComposeFunction((PetscObject)v,"PetscViewerFileGetName_C",PetscViewerFileGetName_ExodusII);
203: PetscObjectComposeFunction((PetscObject)v,"PetscViewerFileSetMode_C",PetscViewerFileSetMode_ExodusII);
204: PetscObjectComposeFunction((PetscObject)v,"PetscViewerFileGetMode_C",PetscViewerFileGetMode_ExodusII);
205: PetscObjectComposeFunction((PetscObject)v,"PetscViewerGetId_C",PetscViewerExodusIIGetId_ExodusII);
206: PetscObjectComposeFunction((PetscObject)v,"PetscViewerSetOrder_C",PetscViewerExodusIISetOrder_ExodusII);
207: PetscObjectComposeFunction((PetscObject)v,"PetscViewerGetOrder_C",PetscViewerExodusIIGetOrder_ExodusII);
208: return 0;
209: }
211: /*
212: EXOGetVarIndex - Locate a result in an exodus file based on its name
214: Collective
216: Input Parameters:
217: + exoid - the exodus id of a file (obtained from ex_open or ex_create for instance)
218: . obj_type - the type of entity for instance EX_NODAL, EX_ELEM_BLOCK
219: - name - the name of the result
221: Output Parameters:
222: . varIndex - the location in the exodus file of the result
224: Notes:
225: The exodus variable index is obtained by comparing name and the
226: names of zonal variables declared in the exodus file. For instance if name is "V"
227: the location in the exodus file will be the first match of "V", "V_X", "V_XX", "V_1", or "V_11"
228: amongst all variables of type obj_type.
230: Level: beginner
232: .seealso: EXOGetVarIndex(),DMPlexView_ExodusII_Internal(),VecViewPlex_ExodusII_Nodal_Internal(),VecLoadNodal_PlexEXO(),VecLoadZonal_PlexEXO()
233: */
234: PetscErrorCode EXOGetVarIndex_Internal(int exoid, ex_entity_type obj_type, const char name[], int *varIndex)
235: {
236: int num_vars, i, j;
237: char ext_name[MAX_STR_LENGTH+1], var_name[MAX_STR_LENGTH+1];
238: const int num_suffix = 5;
239: char *suffix[5];
240: PetscBool flg;
242: suffix[0] = (char *) "";
243: suffix[1] = (char *) "_X";
244: suffix[2] = (char *) "_XX";
245: suffix[3] = (char *) "_1";
246: suffix[4] = (char *) "_11";
248: *varIndex = -1;
249: PetscStackCallStandard(ex_get_variable_param,exoid, obj_type, &num_vars);
250: for (i = 0; i < num_vars; ++i) {
251: PetscStackCallStandard(ex_get_variable_name,exoid, obj_type, i+1, var_name);
252: for (j = 0; j < num_suffix; ++j) {
253: PetscStrncpy(ext_name, name, MAX_STR_LENGTH);
254: PetscStrlcat(ext_name, suffix[j], MAX_STR_LENGTH);
255: PetscStrcasecmp(ext_name, var_name, &flg);
256: if (flg) {
257: *varIndex = i+1;
258: }
259: }
260: }
261: return 0;
262: }
264: /*
265: DMView_PlexExodusII - Write a DM to disk in exodus format
267: Collective on dm
269: Input Parameters:
270: + dm - The dm to be written
271: . viewer - an exodusII viewer
273: Notes:
274: Not all DM can be written to disk this way. For instance, exodus assume that element blocks (mapped to "Cell sets" labels)
275: consists of sequentially numbered cells. If this is not the case, the exodus file will be corrupted.
277: If the dm has been distributed, only the part of the DM on MPI rank 0 (including "ghost" cells and vertices)
278: will be written.
280: DMPlex only represents geometry while most post-processing software expect that a mesh also provides information
281: on the discretization space. This function assumes that the file represents Lagrange finite elements of order 1 or 2.
282: The order of the mesh shall be set using PetscViewerExodusIISetOrder
283: It should be extended to use PetscFE objects.
285: This function will only handle TRI, TET, QUAD, and HEX cells.
286: Level: beginner
288: .seealso:
289: */
290: PetscErrorCode DMView_PlexExodusII(DM dm, PetscViewer viewer)
291: {
292: enum ElemType {TRI, QUAD, TET, HEX};
293: MPI_Comm comm;
294: PetscInt degree; /* the order of the mesh */
295: /* Connectivity Variables */
296: PetscInt cellsNotInConnectivity;
297: /* Cell Sets */
298: DMLabel csLabel;
299: IS csIS;
300: const PetscInt *csIdx;
301: PetscInt num_cs, cs;
302: enum ElemType *type;
303: PetscBool hasLabel;
304: /* Coordinate Variables */
305: DM cdm;
306: PetscSection coordSection;
307: Vec coord;
308: PetscInt **nodes;
309: PetscInt depth, d, dim, skipCells = 0;
310: PetscInt pStart, pEnd, p, cStart, cEnd, numCells, vStart, vEnd, numVertices, eStart, eEnd, numEdges, fStart, fEnd, numFaces, numNodes;
311: PetscInt num_vs, num_fs;
312: PetscMPIInt rank, size;
313: const char *dmName;
314: PetscInt nodesTriP1[4] = {3,0,0,0};
315: PetscInt nodesTriP2[4] = {3,3,0,0};
316: PetscInt nodesQuadP1[4] = {4,0,0,0};
317: PetscInt nodesQuadP2[4] = {4,4,0,1};
318: PetscInt nodesTetP1[4] = {4,0,0,0};
319: PetscInt nodesTetP2[4] = {4,6,0,0};
320: PetscInt nodesHexP1[4] = {8,0,0,0};
321: PetscInt nodesHexP2[4] = {8,12,6,1};
322: PetscErrorCode ierr;
324: int CPU_word_size, IO_word_size, EXO_mode;
325: float EXO_version;
327: PetscViewer_ExodusII *exo = (PetscViewer_ExodusII *) viewer->data;
329: PetscObjectGetComm((PetscObject)dm, &comm);
330: MPI_Comm_rank(comm, &rank);
331: MPI_Comm_size(comm, &size);
333: /*
334: Creating coordSection is a collective operation so we do it somewhat out of sequence
335: */
336: PetscSectionCreate(comm, &coordSection);
337: DMGetCoordinatesLocalSetUp(dm);
338: if (!rank) {
339: switch (exo->btype) {
340: case FILE_MODE_READ:
341: case FILE_MODE_APPEND:
342: case FILE_MODE_UPDATE:
343: case FILE_MODE_APPEND_UPDATE:
344: /* exodusII does not allow writing geometry to an existing file */
345: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "cannot add geometry to existing file %s", exo->filename);
346: case FILE_MODE_WRITE:
347: /* Create an empty file if one already exists*/
348: EXO_mode = EX_CLOBBER;
349: #if defined(PETSC_USE_64BIT_INDICES)
350: EXO_mode += EX_ALL_INT64_API;
351: #endif
352: CPU_word_size = sizeof(PetscReal);
353: IO_word_size = sizeof(PetscReal);
354: exo->exoid = ex_create(exo->filename, EXO_mode, &CPU_word_size, &IO_word_size);
357: break;
358: default:
359: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER, "Must call PetscViewerFileSetMode() before PetscViewerFileSetName()");
360: }
362: /* --- Get DM info --- */
363: PetscObjectGetName((PetscObject) dm, &dmName);
364: DMPlexGetDepth(dm, &depth);
365: DMGetDimension(dm, &dim);
366: DMPlexGetChart(dm, &pStart, &pEnd);
367: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
368: DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);
369: DMPlexGetDepthStratum(dm, 1, &eStart, &eEnd);
370: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
371: numCells = cEnd - cStart;
372: numEdges = eEnd - eStart;
373: numVertices = vEnd - vStart;
374: if (depth == 3) {numFaces = fEnd - fStart;}
375: else {numFaces = 0;}
376: DMGetLabelSize(dm, "Cell Sets", &num_cs);
377: DMGetLabelSize(dm, "Vertex Sets", &num_vs);
378: DMGetLabelSize(dm, "Face Sets", &num_fs);
379: DMGetCoordinatesLocal(dm, &coord);
380: DMGetCoordinateDM(dm, &cdm);
381: if (num_cs > 0) {
382: DMGetLabel(dm, "Cell Sets", &csLabel);
383: DMLabelGetValueIS(csLabel, &csIS);
384: ISGetIndices(csIS, &csIdx);
385: }
386: PetscMalloc1(num_cs, &nodes);
387: /* Set element type for each block and compute total number of nodes */
388: PetscMalloc1(num_cs, &type);
389: numNodes = numVertices;
391: PetscViewerExodusIIGetOrder(viewer, °ree);
392: if (degree == 2) {numNodes += numEdges;}
393: cellsNotInConnectivity = numCells;
394: for (cs = 0; cs < num_cs; ++cs) {
395: IS stratumIS;
396: const PetscInt *cells;
397: PetscScalar *xyz = NULL;
398: PetscInt csSize, closureSize;
400: DMLabelGetStratumIS(csLabel, csIdx[cs], &stratumIS);
401: ISGetIndices(stratumIS, &cells);
402: ISGetSize(stratumIS, &csSize);
403: DMPlexVecGetClosure(cdm, NULL, coord, cells[0], &closureSize, &xyz);
404: switch (dim) {
405: case 2:
406: if (closureSize == 3*dim) {type[cs] = TRI;}
407: else if (closureSize == 4*dim) {type[cs] = QUAD;}
408: else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Number of vertices %D in dimension %D has no ExodusII type", closureSize/dim, dim);
409: break;
410: case 3:
411: if (closureSize == 4*dim) {type[cs] = TET;}
412: else if (closureSize == 8*dim) {type[cs] = HEX;}
413: else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Number of vertices %D in dimension %D has no ExodusII type", closureSize/dim, dim);
414: break;
415: default: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Dimension %D not handled by ExodusII viewer", dim);
416: }
417: if ((degree == 2) && (type[cs] == QUAD)) {numNodes += csSize;}
418: if ((degree == 2) && (type[cs] == HEX)) {numNodes += csSize; numNodes += numFaces;}
419: DMPlexVecRestoreClosure(cdm, NULL, coord, cells[0], &closureSize, &xyz);
420: /* Set nodes and Element type */
421: if (type[cs] == TRI) {
422: if (degree == 1) nodes[cs] = nodesTriP1;
423: else if (degree == 2) nodes[cs] = nodesTriP2;
424: } else if (type[cs] == QUAD) {
425: if (degree == 1) nodes[cs] = nodesQuadP1;
426: else if (degree == 2) nodes[cs] = nodesQuadP2;
427: } else if (type[cs] == TET) {
428: if (degree == 1) nodes[cs] = nodesTetP1;
429: else if (degree == 2) nodes[cs] = nodesTetP2;
430: } else if (type[cs] == HEX) {
431: if (degree == 1) nodes[cs] = nodesHexP1;
432: else if (degree == 2) nodes[cs] = nodesHexP2;
433: }
434: /* Compute the number of cells not in the connectivity table */
435: cellsNotInConnectivity -= nodes[cs][3]*csSize;
437: ISRestoreIndices(stratumIS, &cells);
438: ISDestroy(&stratumIS);
439: }
440: if (num_cs > 0) {PetscStackCallStandard(ex_put_init,exo->exoid, dmName, dim, numNodes, numCells, num_cs, num_vs, num_fs);}
441: /* --- Connectivity --- */
442: for (cs = 0; cs < num_cs; ++cs) {
443: IS stratumIS;
444: const PetscInt *cells;
445: PetscInt *connect, off = 0;
446: PetscInt edgesInClosure = 0, facesInClosure = 0, verticesInClosure = 0;
447: PetscInt csSize, c, connectSize, closureSize;
448: char *elem_type = NULL;
449: char elem_type_tri3[] = "TRI3", elem_type_quad4[] = "QUAD4";
450: char elem_type_tri6[] = "TRI6", elem_type_quad9[] = "QUAD9";
451: char elem_type_tet4[] = "TET4", elem_type_hex8[] = "HEX8";
452: char elem_type_tet10[] = "TET10", elem_type_hex27[] = "HEX27";
454: DMLabelGetStratumIS(csLabel, csIdx[cs], &stratumIS);
455: ISGetIndices(stratumIS, &cells);
456: ISGetSize(stratumIS, &csSize);
457: /* Set Element type */
458: if (type[cs] == TRI) {
459: if (degree == 1) elem_type = elem_type_tri3;
460: else if (degree == 2) elem_type = elem_type_tri6;
461: } else if (type[cs] == QUAD) {
462: if (degree == 1) elem_type = elem_type_quad4;
463: else if (degree == 2) elem_type = elem_type_quad9;
464: } else if (type[cs] == TET) {
465: if (degree == 1) elem_type = elem_type_tet4;
466: else if (degree == 2) elem_type = elem_type_tet10;
467: } else if (type[cs] == HEX) {
468: if (degree == 1) elem_type = elem_type_hex8;
469: else if (degree == 2) elem_type = elem_type_hex27;
470: }
471: connectSize = nodes[cs][0] + nodes[cs][1] + nodes[cs][2] + nodes[cs][3];
472: PetscMalloc1(PetscMax(27,connectSize)*csSize, &connect);
473: PetscStackCallStandard(ex_put_block,exo->exoid, EX_ELEM_BLOCK, csIdx[cs], elem_type, csSize, connectSize, 0, 0, 1);
474: /* Find number of vertices, edges, and faces in the closure */
475: verticesInClosure = nodes[cs][0];
476: if (depth > 1) {
477: if (dim == 2) {
478: DMPlexGetConeSize(dm, cells[0], &edgesInClosure);
479: } else if (dim == 3) {
480: PetscInt *closure = NULL;
482: DMPlexGetConeSize(dm, cells[0], &facesInClosure);
483: DMPlexGetTransitiveClosure(dm, cells[0], PETSC_TRUE, &closureSize, &closure);
484: edgesInClosure = closureSize - facesInClosure - 1 - verticesInClosure;
485: DMPlexRestoreTransitiveClosure(dm, cells[0], PETSC_TRUE, &closureSize, &closure);
486: }
487: }
488: /* Get connectivity for each cell */
489: for (c = 0; c < csSize; ++c) {
490: PetscInt *closure = NULL;
491: PetscInt temp, i;
493: DMPlexGetTransitiveClosure(dm, cells[c], PETSC_TRUE, &closureSize, &closure);
494: for (i = 0; i < connectSize; ++i) {
495: if (i < nodes[cs][0]) {/* Vertices */
496: connect[i+off] = closure[(i+edgesInClosure+facesInClosure+1)*2] + 1;
497: connect[i+off] -= cellsNotInConnectivity;
498: } else if (i < nodes[cs][0]+nodes[cs][1]) { /* Edges */
499: connect[i+off] = closure[(i-verticesInClosure+facesInClosure+1)*2] + 1;
500: if (nodes[cs][2] == 0) connect[i+off] -= numFaces;
501: connect[i+off] -= cellsNotInConnectivity;
502: } else if (i < nodes[cs][0]+nodes[cs][1]+nodes[cs][3]) { /* Cells */
503: connect[i+off] = closure[0] + 1;
504: connect[i+off] -= skipCells;
505: } else if (i < nodes[cs][0]+nodes[cs][1]+nodes[cs][3]+nodes[cs][2]) { /* Faces */
506: connect[i+off] = closure[(i-edgesInClosure-verticesInClosure)*2] + 1;
507: connect[i+off] -= cellsNotInConnectivity;
508: } else {
509: connect[i+off] = -1;
510: }
511: }
512: /* Tetrahedra are inverted */
513: if (type[cs] == TET) {
514: temp = connect[0+off]; connect[0+off] = connect[1+off]; connect[1+off] = temp;
515: if (degree == 2) {
516: temp = connect[5+off]; connect[5+off] = connect[6+off]; connect[6+off] = temp;
517: temp = connect[7+off]; connect[7+off] = connect[8+off]; connect[8+off] = temp;
518: }
519: }
520: /* Hexahedra are inverted */
521: if (type[cs] == HEX) {
522: temp = connect[1+off]; connect[1+off] = connect[3+off]; connect[3+off] = temp;
523: if (degree == 2) {
524: temp = connect[8+off]; connect[8+off] = connect[11+off]; connect[11+off] = temp;
525: temp = connect[9+off]; connect[9+off] = connect[10+off]; connect[10+off] = temp;
526: temp = connect[16+off]; connect[16+off] = connect[17+off]; connect[17+off] = temp;
527: temp = connect[18+off]; connect[18+off] = connect[19+off]; connect[19+off] = temp;
529: temp = connect[12+off]; connect[12+off] = connect[16+off]; connect[16+off] = temp;
530: temp = connect[13+off]; connect[13+off] = connect[17+off]; connect[17+off] = temp;
531: temp = connect[14+off]; connect[14+off] = connect[18+off]; connect[18+off] = temp;
532: temp = connect[15+off]; connect[15+off] = connect[19+off]; connect[19+off] = temp;
534: temp = connect[23+off]; connect[23+off] = connect[26+off]; connect[26+off] = temp;
535: temp = connect[24+off]; connect[24+off] = connect[25+off]; connect[25+off] = temp;
536: temp = connect[25+off]; connect[25+off] = connect[26+off]; connect[26+off] = temp;
537: }
538: }
539: off += connectSize;
540: DMPlexRestoreTransitiveClosure(dm, cells[c], PETSC_TRUE, &closureSize, &closure);
541: }
542: PetscStackCallStandard(ex_put_conn,exo->exoid, EX_ELEM_BLOCK, csIdx[cs], connect, 0, 0);
543: skipCells += (nodes[cs][3] == 0)*csSize;
544: PetscFree(connect);
545: ISRestoreIndices(stratumIS, &cells);
546: ISDestroy(&stratumIS);
547: }
548: PetscFree(type);
549: /* --- Coordinates --- */
550: PetscSectionSetChart(coordSection, pStart, pEnd);
551: if (num_cs) {
552: for (d = 0; d < depth; ++d) {
553: DMPlexGetDepthStratum(dm, d, &pStart, &pEnd);
554: for (p = pStart; p < pEnd; ++p) {
555: PetscSectionSetDof(coordSection, p, nodes[0][d] > 0);
556: }
557: }
558: }
559: for (cs = 0; cs < num_cs; ++cs) {
560: IS stratumIS;
561: const PetscInt *cells;
562: PetscInt csSize, c;
564: DMLabelGetStratumIS(csLabel, csIdx[cs], &stratumIS);
565: ISGetIndices(stratumIS, &cells);
566: ISGetSize(stratumIS, &csSize);
567: for (c = 0; c < csSize; ++c) {
568: PetscSectionSetDof(coordSection, cells[c], nodes[cs][3] > 0);
569: }
570: ISRestoreIndices(stratumIS, &cells);
571: ISDestroy(&stratumIS);
572: }
573: if (num_cs > 0) {
574: ISRestoreIndices(csIS, &csIdx);
575: ISDestroy(&csIS);
576: }
577: PetscFree(nodes);
578: PetscSectionSetUp(coordSection);
579: if (numNodes > 0) {
580: const char *coordNames[3] = {"x", "y", "z"};
581: PetscScalar *closure, *cval;
582: PetscReal *coords;
583: PetscInt hasDof, n = 0;
585: /* There can't be more than 24 values in the closure of a point for the coord coordSection */
586: PetscCalloc3(numNodes*3, &coords, dim, &cval, 24, &closure);
587: DMGetCoordinatesLocalNoncollective(dm, &coord);
588: DMPlexGetChart(dm, &pStart, &pEnd);
589: for (p = pStart; p < pEnd; ++p) {
590: PetscSectionGetDof(coordSection, p, &hasDof);
591: if (hasDof) {
592: PetscInt closureSize = 24, j;
594: DMPlexVecGetClosure(cdm, NULL, coord, p, &closureSize, &closure);
595: for (d = 0; d < dim; ++d) {
596: cval[d] = 0.0;
597: for (j = 0; j < closureSize/dim; j++) cval[d] += closure[j*dim+d];
598: coords[d*numNodes+n] = PetscRealPart(cval[d]) * dim / closureSize;
599: }
600: ++n;
601: }
602: }
603: PetscStackCallStandard(ex_put_coord,exo->exoid, &coords[0*numNodes], &coords[1*numNodes], &coords[2*numNodes]);
604: PetscFree3(coords, cval, closure);
605: PetscStackCallStandard(ex_put_coord_names,exo->exoid, (char **) coordNames);
606: }
608: /* --- Node Sets/Vertex Sets --- */
609: DMHasLabel(dm, "Vertex Sets", &hasLabel);
610: if (hasLabel) {
611: PetscInt i, vs, vsSize;
612: const PetscInt *vsIdx, *vertices;
613: PetscInt *nodeList;
614: IS vsIS, stratumIS;
615: DMLabel vsLabel;
616: DMGetLabel(dm, "Vertex Sets", &vsLabel);
617: DMLabelGetValueIS(vsLabel, &vsIS);
618: ISGetIndices(vsIS, &vsIdx);
619: for (vs=0; vs<num_vs; ++vs) {
620: DMLabelGetStratumIS(vsLabel, vsIdx[vs], &stratumIS);
621: ISGetIndices(stratumIS, &vertices);
622: ISGetSize(stratumIS, &vsSize);
623: PetscMalloc1(vsSize, &nodeList);
624: for (i=0; i<vsSize; ++i) {
625: nodeList[i] = vertices[i] - skipCells + 1;
626: }
627: PetscStackCallStandard(ex_put_set_param,exo->exoid, EX_NODE_SET, vsIdx[vs], vsSize, 0);
628: PetscStackCallStandard(ex_put_set,exo->exoid, EX_NODE_SET, vsIdx[vs], nodeList, NULL);
629: ISRestoreIndices(stratumIS, &vertices);
630: ISDestroy(&stratumIS);
631: PetscFree(nodeList);
632: }
633: ISRestoreIndices(vsIS, &vsIdx);
634: ISDestroy(&vsIS);
635: }
636: /* --- Side Sets/Face Sets --- */
637: DMHasLabel(dm, "Face Sets", &hasLabel);
638: if (hasLabel) {
639: PetscInt i, j, fs, fsSize;
640: const PetscInt *fsIdx, *faces;
641: IS fsIS, stratumIS;
642: DMLabel fsLabel;
643: PetscInt numPoints, *points;
644: PetscInt elem_list_size = 0;
645: PetscInt *elem_list, *elem_ind, *side_list;
647: DMGetLabel(dm, "Face Sets", &fsLabel);
648: /* Compute size of Node List and Element List */
649: DMLabelGetValueIS(fsLabel, &fsIS);
650: ISGetIndices(fsIS, &fsIdx);
651: for (fs=0; fs<num_fs; ++fs) {
652: DMLabelGetStratumIS(fsLabel, fsIdx[fs], &stratumIS);
653: ISGetSize(stratumIS, &fsSize);
654: elem_list_size += fsSize;
655: ISDestroy(&stratumIS);
656: }
657: PetscMalloc3(num_fs, &elem_ind,
658: elem_list_size, &elem_list,
659: elem_list_size, &side_list);
660: elem_ind[0] = 0;
661: for (fs=0; fs<num_fs; ++fs) {
662: DMLabelGetStratumIS(fsLabel, fsIdx[fs], &stratumIS);
663: ISGetIndices(stratumIS, &faces);
664: ISGetSize(stratumIS, &fsSize);
665: /* Set Parameters */
666: PetscStackCallStandard(ex_put_set_param,exo->exoid, EX_SIDE_SET, fsIdx[fs], fsSize, 0);
667: /* Indices */
668: if (fs<num_fs-1) {
669: elem_ind[fs+1] = elem_ind[fs] + fsSize;
670: }
672: for (i=0; i<fsSize; ++i) {
673: /* Element List */
674: points = NULL;
675: DMPlexGetTransitiveClosure(dm, faces[i], PETSC_FALSE, &numPoints, &points);
676: elem_list[elem_ind[fs] + i] = points[2] +1;
677: DMPlexRestoreTransitiveClosure(dm, faces[i], PETSC_FALSE, &numPoints, &points);
679: /* Side List */
680: points = NULL;
681: DMPlexGetTransitiveClosure(dm, elem_list[elem_ind[fs] + i]-1, PETSC_TRUE, &numPoints, &points);
682: for (j=1; j<numPoints; ++j) {
683: if (points[j*2]==faces[i]) {break;}
684: }
685: /* Convert HEX sides */
686: if (numPoints == 27) {
687: if (j == 1) {j = 5;}
688: else if (j == 2) {j = 6;}
689: else if (j == 3) {j = 1;}
690: else if (j == 4) {j = 3;}
691: else if (j == 5) {j = 2;}
692: else if (j == 6) {j = 4;}
693: }
694: /* Convert TET sides */
695: if (numPoints == 15) {
696: --j;
697: if (j == 0) {j = 4;}
698: }
699: side_list[elem_ind[fs] + i] = j;
700: DMPlexRestoreTransitiveClosure(dm, elem_list[elem_ind[fs] + i]-1, PETSC_TRUE, &numPoints, &points);
702: }
703: ISRestoreIndices(stratumIS, &faces);
704: ISDestroy(&stratumIS);
705: }
706: ISRestoreIndices(fsIS, &fsIdx);
707: ISDestroy(&fsIS);
709: /* Put side sets */
710: for (fs=0; fs<num_fs; ++fs) {
711: PetscStackCallStandard(ex_put_set,exo->exoid, EX_SIDE_SET, fsIdx[fs], &elem_list[elem_ind[fs]], &side_list[elem_ind[fs]]);
712: }
713: PetscFree3(elem_ind,elem_list,side_list);
714: }
715: /*
716: close the exodus file
717: */
718: ex_close(exo->exoid);
719: exo->exoid = -1;
720: }
721: PetscSectionDestroy(&coordSection);
723: /*
724: reopen the file in parallel
725: */
726: EXO_mode = EX_WRITE;
727: #if defined(PETSC_USE_64BIT_INDICES)
728: EXO_mode += EX_ALL_INT64_API;
729: #endif
730: CPU_word_size = sizeof(PetscReal);
731: IO_word_size = sizeof(PetscReal);
732: exo->exoid = ex_open_par(exo->filename,EXO_mode,&CPU_word_size,&IO_word_size,&EXO_version,comm,MPI_INFO_NULL);
734: return 0;
735: }
737: /*
738: VecView_PlexExodusII_Internal - Write a Vec corresponding in an exodus file
740: Collective on v
742: Input Parameters:
743: + v - The vector to be written
744: - viewer - The PetscViewerExodusII viewer associate to an exodus file
746: Notes:
747: The exodus result variable index is obtained by comparing the Vec name and the
748: names of variables declared in the exodus file. For instance for a Vec named "V"
749: the location in the exodus file will be the first match of "V", "V_X", "V_XX", "V_1", or "V_11"
750: amongst all variables.
751: In the event where a nodal and zonal variable both match, the function will return an error instead of
752: possibly corrupting the file
754: Level: beginner
756: .seealso: EXOGetVarIndex_Internal(),DMPlexView_ExodusII(),VecView_PlexExodusII()
757: @*/
758: PetscErrorCode VecView_PlexExodusII_Internal(Vec v, PetscViewer viewer)
759: {
760: DM dm;
761: MPI_Comm comm;
762: PetscMPIInt rank;
764: int exoid,offsetN = 0, offsetZ = 0;
765: const char *vecname;
766: PetscInt step;
768: PetscObjectGetComm((PetscObject) v, &comm);
769: MPI_Comm_rank(comm, &rank);
770: PetscViewerExodusIIGetId(viewer,&exoid);
771: VecGetDM(v, &dm);
772: PetscObjectGetName((PetscObject) v, &vecname);
774: DMGetOutputSequenceNumber(dm,&step,NULL);
775: EXOGetVarIndex_Internal(exoid,EX_NODAL,vecname,&offsetN);
776: EXOGetVarIndex_Internal(exoid,EX_ELEM_BLOCK,vecname,&offsetZ);
778: if (offsetN > 0) {
779: VecViewPlex_ExodusII_Nodal_Internal(v,exoid,(int) step+1,offsetN);
780: } else if (offsetZ > 0) {
781: VecViewPlex_ExodusII_Zonal_Internal(v,exoid,(int) step+1,offsetZ);
782: } else SETERRQ(comm, PETSC_ERR_FILE_UNEXPECTED, "Could not find nodal or zonal variable %s in exodus file. ", vecname);
783: return 0;
784: }
786: /*
787: VecLoad_PlexExodusII_Internal - Write a Vec corresponding in an exodus file
789: Collective on v
791: Input Parameters:
792: + v - The vector to be written
793: - viewer - The PetscViewerExodusII viewer associate to an exodus file
795: Notes:
796: The exodus result variable index is obtained by comparing the Vec name and the
797: names of variables declared in the exodus file. For instance for a Vec named "V"
798: the location in the exodus file will be the first match of "V", "V_X", "V_XX", "V_1", or "V_11"
799: amongst all variables.
800: In the event where a nodal and zonal variable both match, the function will return an error instead of
801: possibly corrupting the file
803: Level: beginner
805: .seealso: EXOGetVarIndex_Internal(),DMPlexView_ExodusII(),VecView_PlexExodusII()
806: @*/
807: PetscErrorCode VecLoad_PlexExodusII_Internal(Vec v, PetscViewer viewer)
808: {
809: DM dm;
810: MPI_Comm comm;
811: PetscMPIInt rank;
813: int exoid,offsetN = 0, offsetZ = 0;
814: const char *vecname;
815: PetscInt step;
817: PetscObjectGetComm((PetscObject) v, &comm);
818: MPI_Comm_rank(comm, &rank);
819: PetscViewerExodusIIGetId(viewer,&exoid);
820: VecGetDM(v, &dm);
821: PetscObjectGetName((PetscObject) v, &vecname);
823: DMGetOutputSequenceNumber(dm,&step,NULL);
824: EXOGetVarIndex_Internal(exoid,EX_NODAL,vecname,&offsetN);
825: EXOGetVarIndex_Internal(exoid,EX_ELEM_BLOCK,vecname,&offsetZ);
827: if (offsetN > 0) {
828: VecLoadPlex_ExodusII_Nodal_Internal(v,exoid,(int) step+1,offsetN);
829: } else if (offsetZ > 0) {
830: VecLoadPlex_ExodusII_Zonal_Internal(v,exoid,(int) step+1,offsetZ);
831: } else SETERRQ(comm, PETSC_ERR_FILE_UNEXPECTED, "Could not find nodal or zonal variable %s in exodus file. ", vecname);
832: return 0;
833: }
835: /*
836: VecViewPlex_ExodusII_Nodal_Internal - Write a Vec corresponding to a nodal field to an exodus file
838: Collective on v
840: Input Parameters:
841: + v - The vector to be written
842: . exoid - the exodus id of a file (obtained from ex_open or ex_create for instance)
843: . step - the time step to write at (exodus steps are numbered starting from 1)
844: - offset - the location of the variable in the file
846: Notes:
847: The exodus result nodal variable index is obtained by comparing the Vec name and the
848: names of nodal variables declared in the exodus file. For instance for a Vec named "V"
849: the location in the exodus file will be the first match of "V", "V_X", "V_XX", "V_1", or "V_11"
850: amongst all nodal variables.
852: Level: beginner
854: .seealso: EXOGetVarIndex_Internal(),DMPlexView_ExodusII_Internal(),VecLoadNodal_PlexEXO(),VecViewZonal_PlexEXO(),VecLoadZonal_PlexEXO()
855: @*/
856: PetscErrorCode VecViewPlex_ExodusII_Nodal_Internal(Vec v, int exoid, int step, int offset)
857: {
858: MPI_Comm comm;
859: PetscMPIInt size;
860: DM dm;
861: Vec vNatural, vComp;
862: const PetscScalar *varray;
863: PetscInt xs, xe, bs;
864: PetscBool useNatural;
866: PetscObjectGetComm((PetscObject) v, &comm);
867: MPI_Comm_size(comm, &size);
868: VecGetDM(v, &dm);
869: DMGetUseNatural(dm, &useNatural);
870: useNatural = useNatural && size > 1 ? PETSC_TRUE : PETSC_FALSE;
871: if (useNatural) {
872: DMGetGlobalVector(dm, &vNatural);
873: DMPlexGlobalToNaturalBegin(dm, v, vNatural);
874: DMPlexGlobalToNaturalEnd(dm, v, vNatural);
875: } else {
876: vNatural = v;
877: }
879: /* Write local chunk of the result in the exodus file
880: exodus stores each component of a vector-valued field as a separate variable.
881: We assume that they are stored sequentially */
882: VecGetOwnershipRange(vNatural, &xs, &xe);
883: VecGetBlockSize(vNatural, &bs);
884: if (bs == 1) {
885: VecGetArrayRead(vNatural, &varray);
886: PetscStackCallStandard(ex_put_partial_var,exoid, step, EX_NODAL, offset, 1, xs+1, xe-xs, varray);
887: VecRestoreArrayRead(vNatural, &varray);
888: } else {
889: IS compIS;
890: PetscInt c;
892: ISCreateStride(comm, (xe-xs)/bs, xs, bs, &compIS);
893: for (c = 0; c < bs; ++c) {
894: ISStrideSetStride(compIS, (xe-xs)/bs, xs+c, bs);
895: VecGetSubVector(vNatural, compIS, &vComp);
896: VecGetArrayRead(vComp, &varray);
897: PetscStackCallStandard(ex_put_partial_var,exoid, step, EX_NODAL, offset+c, 1, xs/bs+1, (xe-xs)/bs, varray);
898: VecRestoreArrayRead(vComp, &varray);
899: VecRestoreSubVector(vNatural, compIS, &vComp);
900: }
901: ISDestroy(&compIS);
902: }
903: if (useNatural) DMRestoreGlobalVector(dm, &vNatural);
904: return 0;
905: }
907: /*
908: VecLoadPlex_ExodusII_Nodal_Internal - Read a Vec corresponding to a nodal field from an exodus file
910: Collective on v
912: Input Parameters:
913: + v - The vector to be written
914: . exoid - the exodus id of a file (obtained from ex_open or ex_create for instance)
915: . step - the time step to read at (exodus steps are numbered starting from 1)
916: - offset - the location of the variable in the file
918: Notes:
919: The exodus result nodal variable index is obtained by comparing the Vec name and the
920: names of nodal variables declared in the exodus file. For instance for a Vec named "V"
921: the location in the exodus file will be the first match of "V", "V_X", "V_XX", "V_1", or "V_11"
922: amongst all nodal variables.
924: Level: beginner
926: .seealso: EXOGetVarIndex_Internal(), DMPlexView_ExodusII_Internal(), VecViewPlex_ExodusII_Nodal_Internal(), VecViewZonal_PlexEXO(), VecLoadZonal_PlexEXO()
927: */
928: PetscErrorCode VecLoadPlex_ExodusII_Nodal_Internal(Vec v, int exoid, int step, int offset)
929: {
930: MPI_Comm comm;
931: PetscMPIInt size;
932: DM dm;
933: Vec vNatural, vComp;
934: PetscScalar *varray;
935: PetscInt xs, xe, bs;
936: PetscBool useNatural;
938: PetscObjectGetComm((PetscObject) v, &comm);
939: MPI_Comm_size(comm, &size);
940: VecGetDM(v,&dm);
941: DMGetUseNatural(dm, &useNatural);
942: useNatural = useNatural && size > 1 ? PETSC_TRUE : PETSC_FALSE;
943: if (useNatural) DMGetGlobalVector(dm,&vNatural);
944: else {vNatural = v;}
946: /* Read local chunk from the file */
947: VecGetOwnershipRange(vNatural, &xs, &xe);
948: VecGetBlockSize(vNatural, &bs);
949: if (bs == 1) {
950: VecGetArray(vNatural, &varray);
951: PetscStackCallStandard(ex_get_partial_var,exoid, step, EX_NODAL, offset, 1, xs+1, xe-xs, varray);
952: VecRestoreArray(vNatural, &varray);
953: } else {
954: IS compIS;
955: PetscInt c;
957: ISCreateStride(comm, (xe-xs)/bs, xs, bs, &compIS);
958: for (c = 0; c < bs; ++c) {
959: ISStrideSetStride(compIS, (xe-xs)/bs, xs+c, bs);
960: VecGetSubVector(vNatural, compIS, &vComp);
961: VecGetArray(vComp, &varray);
962: PetscStackCallStandard(ex_get_partial_var,exoid, step, EX_NODAL, offset+c, 1, xs/bs+1, (xe-xs)/bs, varray);
963: VecRestoreArray(vComp, &varray);
964: VecRestoreSubVector(vNatural, compIS, &vComp);
965: }
966: ISDestroy(&compIS);
967: }
968: if (useNatural) {
969: DMPlexNaturalToGlobalBegin(dm, vNatural, v);
970: DMPlexNaturalToGlobalEnd(dm, vNatural, v);
971: DMRestoreGlobalVector(dm, &vNatural);
972: }
973: return 0;
974: }
976: /*
977: VecViewPlex_ExodusII_Zonal_Internal - Write a Vec corresponding to a zonal (cell based) field to an exodus file
979: Collective on v
981: Input Parameters:
982: + v - The vector to be written
983: . exoid - the exodus id of a file (obtained from ex_open or ex_create for instance)
984: . step - the time step to write at (exodus steps are numbered starting from 1)
985: - offset - the location of the variable in the file
987: Notes:
988: The exodus result zonal variable index is obtained by comparing the Vec name and the
989: names of zonal variables declared in the exodus file. For instance for a Vec named "V"
990: the location in the exodus file will be the first match of "V", "V_X", "V_XX", "V_1", or "V_11"
991: amongst all zonal variables.
993: Level: beginner
995: .seealso: EXOGetVarIndex_Internal(),DMPlexView_ExodusII_Internal(),VecViewPlex_ExodusII_Nodal_Internal(),VecLoadPlex_ExodusII_Nodal_Internal(),VecLoadPlex_ExodusII_Zonal_Internal()
996: */
997: PetscErrorCode VecViewPlex_ExodusII_Zonal_Internal(Vec v, int exoid, int step, int offset)
998: {
999: MPI_Comm comm;
1000: PetscMPIInt size;
1001: DM dm;
1002: Vec vNatural, vComp;
1003: const PetscScalar *varray;
1004: PetscInt xs, xe, bs;
1005: PetscBool useNatural;
1006: IS compIS;
1007: PetscInt *csSize, *csID;
1008: PetscInt numCS, set, csxs = 0;
1010: PetscObjectGetComm((PetscObject)v, &comm);
1011: MPI_Comm_size(comm, &size);
1012: VecGetDM(v, &dm);
1013: DMGetUseNatural(dm, &useNatural);
1014: useNatural = useNatural && size > 1 ? PETSC_TRUE : PETSC_FALSE;
1015: if (useNatural) {
1016: DMGetGlobalVector(dm, &vNatural);
1017: DMPlexGlobalToNaturalBegin(dm, v, vNatural);
1018: DMPlexGlobalToNaturalEnd(dm, v, vNatural);
1019: } else {
1020: vNatural = v;
1021: }
1023: /* Write local chunk of the result in the exodus file
1024: exodus stores each component of a vector-valued field as a separate variable.
1025: We assume that they are stored sequentially
1026: Zonal variables are accessed one element block at a time, so we loop through the cell sets,
1027: but once the vector has been reordered to natural size, we cannot use the label information
1028: to figure out what to save where. */
1029: numCS = ex_inquire_int(exoid, EX_INQ_ELEM_BLK);
1030: PetscMalloc2(numCS, &csID, numCS, &csSize);
1031: PetscStackCallStandard(ex_get_ids,exoid, EX_ELEM_BLOCK, csID);
1032: for (set = 0; set < numCS; ++set) {
1033: ex_block block;
1035: block.id = csID[set];
1036: block.type = EX_ELEM_BLOCK;
1037: PetscStackCallStandard(ex_get_block_param,exoid, &block);
1038: csSize[set] = block.num_entry;
1039: }
1040: VecGetOwnershipRange(vNatural, &xs, &xe);
1041: VecGetBlockSize(vNatural, &bs);
1042: if (bs > 1) ISCreateStride(comm, (xe-xs)/bs, xs, bs, &compIS);
1043: for (set = 0; set < numCS; set++) {
1044: PetscInt csLocalSize, c;
1046: /* range of indices for set setID[set]: csxs:csxs + csSize[set]-1
1047: local slice of zonal values: xs/bs,xm/bs-1
1048: intersection: max(xs/bs,csxs),min(xm/bs-1,csxs + csSize[set]-1) */
1049: csLocalSize = PetscMax(0, PetscMin(xe/bs, csxs+csSize[set]) - PetscMax(xs/bs, csxs));
1050: if (bs == 1) {
1051: VecGetArrayRead(vNatural, &varray);
1052: PetscStackCallStandard(ex_put_partial_var,exoid, step, EX_ELEM_BLOCK, offset, csID[set], PetscMax(xs-csxs, 0)+1, csLocalSize, &varray[PetscMax(0, csxs-xs)]);
1053: VecRestoreArrayRead(vNatural, &varray);
1054: } else {
1055: for (c = 0; c < bs; ++c) {
1056: ISStrideSetStride(compIS, (xe-xs)/bs, xs+c, bs);
1057: VecGetSubVector(vNatural, compIS, &vComp);
1058: VecGetArrayRead(vComp, &varray);
1059: 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)]);
1060: VecRestoreArrayRead(vComp, &varray);
1061: VecRestoreSubVector(vNatural, compIS, &vComp);
1062: }
1063: }
1064: csxs += csSize[set];
1065: }
1066: PetscFree2(csID, csSize);
1067: if (bs > 1) ISDestroy(&compIS);
1068: if (useNatural) DMRestoreGlobalVector(dm,&vNatural);
1069: return 0;
1070: }
1072: /*
1073: VecLoadPlex_ExodusII_Zonal_Internal - Read a Vec corresponding to a zonal (cell based) field from an exodus file
1075: Collective on v
1077: Input Parameters:
1078: + v - The vector to be written
1079: . exoid - the exodus id of a file (obtained from ex_open or ex_create for instance)
1080: . step - the time step to read at (exodus steps are numbered starting from 1)
1081: - offset - the location of the variable in the file
1083: Notes:
1084: The exodus result zonal variable index is obtained by comparing the Vec name and the
1085: names of zonal variables declared in the exodus file. For instance for a Vec named "V"
1086: the location in the exodus file will be the first match of "V", "V_X", "V_XX", "V_1", or "V_11"
1087: amongst all zonal variables.
1089: Level: beginner
1091: .seealso: EXOGetVarIndex_Internal(), DMPlexView_ExodusII_Internal(), VecViewPlex_ExodusII_Nodal_Internal(), VecLoadPlex_ExodusII_Nodal_Internal(), VecLoadPlex_ExodusII_Zonal_Internal()
1092: */
1093: PetscErrorCode VecLoadPlex_ExodusII_Zonal_Internal(Vec v, int exoid, int step, int offset)
1094: {
1095: MPI_Comm comm;
1096: PetscMPIInt size;
1097: DM dm;
1098: Vec vNatural, vComp;
1099: PetscScalar *varray;
1100: PetscInt xs, xe, bs;
1101: PetscBool useNatural;
1102: IS compIS;
1103: PetscInt *csSize, *csID;
1104: PetscInt numCS, set, csxs = 0;
1106: PetscObjectGetComm((PetscObject)v,&comm);
1107: MPI_Comm_size(comm, &size);
1108: VecGetDM(v, &dm);
1109: DMGetUseNatural(dm, &useNatural);
1110: useNatural = useNatural && size > 1 ? PETSC_TRUE : PETSC_FALSE;
1111: if (useNatural) DMGetGlobalVector(dm,&vNatural);
1112: else {vNatural = v;}
1114: /* Read local chunk of the result in the exodus file
1115: exodus stores each component of a vector-valued field as a separate variable.
1116: We assume that they are stored sequentially
1117: Zonal variables are accessed one element block at a time, so we loop through the cell sets,
1118: but once the vector has been reordered to natural size, we cannot use the label information
1119: to figure out what to save where. */
1120: numCS = ex_inquire_int(exoid, EX_INQ_ELEM_BLK);
1121: PetscMalloc2(numCS, &csID, numCS, &csSize);
1122: PetscStackCallStandard(ex_get_ids,exoid, EX_ELEM_BLOCK, csID);
1123: for (set = 0; set < numCS; ++set) {
1124: ex_block block;
1126: block.id = csID[set];
1127: block.type = EX_ELEM_BLOCK;
1128: PetscStackCallStandard(ex_get_block_param,exoid, &block);
1129: csSize[set] = block.num_entry;
1130: }
1131: VecGetOwnershipRange(vNatural, &xs, &xe);
1132: VecGetBlockSize(vNatural, &bs);
1133: if (bs > 1) ISCreateStride(comm, (xe-xs)/bs, xs, bs, &compIS);
1134: for (set = 0; set < numCS; ++set) {
1135: PetscInt csLocalSize, c;
1137: /* range of indices for set setID[set]: csxs:csxs + csSize[set]-1
1138: local slice of zonal values: xs/bs,xm/bs-1
1139: intersection: max(xs/bs,csxs),min(xm/bs-1,csxs + csSize[set]-1) */
1140: csLocalSize = PetscMax(0, PetscMin(xe/bs, csxs+csSize[set]) - PetscMax(xs/bs, csxs));
1141: if (bs == 1) {
1142: VecGetArray(vNatural, &varray);
1143: PetscStackCallStandard(ex_get_partial_var,exoid, step, EX_ELEM_BLOCK, offset, csID[set], PetscMax(xs-csxs, 0)+1, csLocalSize, &varray[PetscMax(0, csxs-xs)]);
1144: VecRestoreArray(vNatural, &varray);
1145: } else {
1146: for (c = 0; c < bs; ++c) {
1147: ISStrideSetStride(compIS, (xe-xs)/bs, xs+c, bs);
1148: VecGetSubVector(vNatural, compIS, &vComp);
1149: VecGetArray(vComp, &varray);
1150: 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)]);
1151: VecRestoreArray(vComp, &varray);
1152: VecRestoreSubVector(vNatural, compIS, &vComp);
1153: }
1154: }
1155: csxs += csSize[set];
1156: }
1157: PetscFree2(csID, csSize);
1158: if (bs > 1) ISDestroy(&compIS);
1159: if (useNatural) {
1160: DMPlexNaturalToGlobalBegin(dm, vNatural, v);
1161: DMPlexNaturalToGlobalEnd(dm, vNatural, v);
1162: DMRestoreGlobalVector(dm, &vNatural);
1163: }
1164: return 0;
1165: }
1166: #endif
1168: /*@
1169: PetscViewerExodusIIGetId - Get the file id of the ExodusII file
1171: Logically Collective on PetscViewer
1173: Input Parameter:
1174: . viewer - the PetscViewer
1176: Output Parameter:
1177: . exoid - The ExodusII file id
1179: Level: intermediate
1181: .seealso: PetscViewerFileSetMode(), PetscViewerCreate(), PetscViewerSetType(), PetscViewerBinaryOpen()
1182: @*/
1183: PetscErrorCode PetscViewerExodusIIGetId(PetscViewer viewer, int *exoid)
1184: {
1186: PetscTryMethod(viewer, "PetscViewerGetId_C",(PetscViewer,int*),(viewer,exoid));
1187: return 0;
1188: }
1190: /*@
1191: PetscViewerExodusIISetOrder - Set the elements order in the exodusII file.
1193: Collective
1195: Input Parameters:
1196: + viewer - the viewer
1197: - order - elements order
1199: Output Parameter:
1201: Level: beginner
1203: Note:
1205: .seealso: PetscViewerExodusIIGetId(), PetscViewerExodusIIGetOrder(), PetscViewerExodusIISetOrder()
1206: @*/
1207: PetscErrorCode PetscViewerExodusIISetOrder(PetscViewer viewer, PetscInt order)
1208: {
1210: PetscTryMethod(viewer, "PetscViewerSetOrder_C",(PetscViewer,PetscInt),(viewer,order));
1211: return 0;
1212: }
1214: /*@
1215: PetscViewerExodusIIGetOrder - Get the elements order in the exodusII file.
1217: Collective
1219: Input Parameters:
1220: + viewer - the viewer
1221: - order - elements order
1223: Output Parameter:
1225: Level: beginner
1227: Note:
1229: .seealso: PetscViewerExodusIIGetId(), PetscViewerExodusIIGetOrder(), PetscViewerExodusIISetOrder()
1230: @*/
1231: PetscErrorCode PetscViewerExodusIIGetOrder(PetscViewer viewer, PetscInt *order)
1232: {
1234: PetscTryMethod(viewer, "PetscViewerGetOrder_C",(PetscViewer,PetscInt*),(viewer,order));
1235: return 0;
1236: }
1238: /*@C
1239: PetscViewerExodusIIOpen - Opens a file for ExodusII input/output.
1241: Collective
1243: Input Parameters:
1244: + comm - MPI communicator
1245: . name - name of file
1246: - type - type of file
1247: $ FILE_MODE_WRITE - create new file for binary output
1248: $ FILE_MODE_READ - open existing file for binary input
1249: $ FILE_MODE_APPEND - open existing file for binary output
1251: Output Parameter:
1252: . exo - PetscViewer for Exodus II input/output to use with the specified file
1254: Level: beginner
1256: Note:
1257: This PetscViewer should be destroyed with PetscViewerDestroy().
1259: .seealso: PetscViewerPushFormat(), PetscViewerDestroy(),
1260: DMLoad(), PetscFileMode, PetscViewer, PetscViewerSetType(), PetscViewerFileSetMode(), PetscViewerFileSetName()
1261: @*/
1262: PetscErrorCode PetscViewerExodusIIOpen(MPI_Comm comm, const char name[], PetscFileMode type, PetscViewer *exo)
1263: {
1264: PetscViewerCreate(comm, exo);
1265: PetscViewerSetType(*exo, PETSCVIEWEREXODUSII);
1266: PetscViewerFileSetMode(*exo, type);
1267: PetscViewerFileSetName(*exo, name);
1268: PetscViewerSetFromOptions(*exo);
1269: return 0;
1270: }
1272: /*@C
1273: DMPlexCreateExodusFromFile - Create a DMPlex mesh from an ExodusII file.
1275: Collective
1277: Input Parameters:
1278: + comm - The MPI communicator
1279: . filename - The name of the ExodusII file
1280: - interpolate - Create faces and edges in the mesh
1282: Output Parameter:
1283: . dm - The DM object representing the mesh
1285: Level: beginner
1287: .seealso: DMPLEX, DMCreate(), DMPlexCreateExodus()
1288: @*/
1289: PetscErrorCode DMPlexCreateExodusFromFile(MPI_Comm comm, const char filename[], PetscBool interpolate, DM *dm)
1290: {
1291: PetscMPIInt rank;
1292: #if defined(PETSC_HAVE_EXODUSII)
1293: int CPU_word_size = sizeof(PetscReal), IO_word_size = 0, exoid = -1;
1294: float version;
1295: #endif
1298: MPI_Comm_rank(comm, &rank);
1299: #if defined(PETSC_HAVE_EXODUSII)
1300: if (rank == 0) {
1301: exoid = ex_open(filename, EX_READ, &CPU_word_size, &IO_word_size, &version);
1303: }
1304: DMPlexCreateExodus(comm, exoid, interpolate, dm);
1305: if (rank == 0) {PetscStackCallStandard(ex_close,exoid);}
1306: return 0;
1307: #else
1308: SETERRQ(comm, PETSC_ERR_SUP, "This method requires ExodusII support. Reconfigure using --download-exodusii");
1309: #endif
1310: }
1312: #if defined(PETSC_HAVE_EXODUSII)
1313: static PetscErrorCode ExodusGetCellType_Internal(const char *elem_type, DMPolytopeType *ct)
1314: {
1315: PetscBool flg;
1317: *ct = DM_POLYTOPE_UNKNOWN;
1318: PetscStrcmp(elem_type, "TRI", &flg);
1319: if (flg) {*ct = DM_POLYTOPE_TRIANGLE; goto done;}
1320: PetscStrcmp(elem_type, "TRI3", &flg);
1321: if (flg) {*ct = DM_POLYTOPE_TRIANGLE; goto done;}
1322: PetscStrcmp(elem_type, "QUAD", &flg);
1323: if (flg) {*ct = DM_POLYTOPE_QUADRILATERAL; goto done;}
1324: PetscStrcmp(elem_type, "QUAD4", &flg);
1325: if (flg) {*ct = DM_POLYTOPE_QUADRILATERAL; goto done;}
1326: PetscStrcmp(elem_type, "SHELL4", &flg);
1327: if (flg) {*ct = DM_POLYTOPE_QUADRILATERAL; goto done;}
1328: PetscStrcmp(elem_type, "TETRA", &flg);
1329: if (flg) {*ct = DM_POLYTOPE_TETRAHEDRON; goto done;}
1330: PetscStrcmp(elem_type, "TET4", &flg);
1331: if (flg) {*ct = DM_POLYTOPE_TETRAHEDRON; goto done;}
1332: PetscStrcmp(elem_type, "WEDGE", &flg);
1333: if (flg) {*ct = DM_POLYTOPE_TRI_PRISM; goto done;}
1334: PetscStrcmp(elem_type, "HEX", &flg);
1335: if (flg) {*ct = DM_POLYTOPE_HEXAHEDRON; goto done;}
1336: PetscStrcmp(elem_type, "HEX8", &flg);
1337: if (flg) {*ct = DM_POLYTOPE_HEXAHEDRON; goto done;}
1338: PetscStrcmp(elem_type, "HEXAHEDRON", &flg);
1339: if (flg) {*ct = DM_POLYTOPE_HEXAHEDRON; goto done;}
1340: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unrecognized element type %s", elem_type);
1341: done:
1342: return 0;
1343: }
1344: #endif
1346: /*@
1347: DMPlexCreateExodus - Create a DMPlex mesh from an ExodusII file ID.
1349: Collective
1351: Input Parameters:
1352: + comm - The MPI communicator
1353: . exoid - The ExodusII id associated with a exodus file and obtained using ex_open
1354: - interpolate - Create faces and edges in the mesh
1356: Output Parameter:
1357: . dm - The DM object representing the mesh
1359: Level: beginner
1361: .seealso: DMPLEX, DMCreate()
1362: @*/
1363: PetscErrorCode DMPlexCreateExodus(MPI_Comm comm, PetscInt exoid, PetscBool interpolate, DM *dm)
1364: {
1365: #if defined(PETSC_HAVE_EXODUSII)
1366: PetscMPIInt num_proc, rank;
1367: DMLabel cellSets = NULL, faceSets = NULL, vertSets = NULL;
1368: PetscSection coordSection;
1369: Vec coordinates;
1370: PetscScalar *coords;
1371: PetscInt coordSize, v;
1372: /* Read from ex_get_init() */
1373: char title[PETSC_MAX_PATH_LEN+1];
1374: int dim = 0, dimEmbed = 0, numVertices = 0, numCells = 0;
1375: int num_cs = 0, num_vs = 0, num_fs = 0;
1376: #endif
1378: #if defined(PETSC_HAVE_EXODUSII)
1379: MPI_Comm_rank(comm, &rank);
1380: MPI_Comm_size(comm, &num_proc);
1381: DMCreate(comm, dm);
1382: DMSetType(*dm, DMPLEX);
1383: /* Open EXODUS II file and read basic information on rank 0, then broadcast to all processors */
1384: if (rank == 0) {
1385: PetscMemzero(title,PETSC_MAX_PATH_LEN+1);
1386: PetscStackCallStandard(ex_get_init,exoid, title, &dimEmbed, &numVertices, &numCells, &num_cs, &num_vs, &num_fs);
1388: }
1389: MPI_Bcast(title, PETSC_MAX_PATH_LEN+1, MPI_CHAR, 0, comm);
1390: MPI_Bcast(&dim, 1, MPI_INT, 0, comm);
1391: PetscObjectSetName((PetscObject) *dm, title);
1392: DMPlexSetChart(*dm, 0, numCells+numVertices);
1393: /* We do not want this label automatically computed, instead we compute it here */
1394: DMCreateLabel(*dm, "celltype");
1396: /* Read cell sets information */
1397: if (rank == 0) {
1398: PetscInt *cone;
1399: int c, cs, ncs, c_loc, v, v_loc;
1400: /* Read from ex_get_elem_blk_ids() */
1401: int *cs_id, *cs_order;
1402: /* Read from ex_get_elem_block() */
1403: char buffer[PETSC_MAX_PATH_LEN+1];
1404: int num_cell_in_set, num_vertex_per_cell, num_hybrid, num_attr;
1405: /* Read from ex_get_elem_conn() */
1406: int *cs_connect;
1408: /* Get cell sets IDs */
1409: PetscMalloc2(num_cs, &cs_id, num_cs, &cs_order);
1410: PetscStackCallStandard(ex_get_ids,exoid, EX_ELEM_BLOCK, cs_id);
1411: /* Read the cell set connectivity table and build mesh topology
1412: EXO standard requires that cells in cell sets be numbered sequentially and be pairwise disjoint. */
1413: /* Check for a hybrid mesh */
1414: for (cs = 0, num_hybrid = 0; cs < num_cs; ++cs) {
1415: DMPolytopeType ct;
1416: char elem_type[PETSC_MAX_PATH_LEN];
1418: PetscArrayzero(elem_type, sizeof(elem_type));
1419: PetscStackCallStandard(ex_get_elem_type,exoid, cs_id[cs], elem_type);
1420: ExodusGetCellType_Internal(elem_type, &ct);
1421: dim = PetscMax(dim, DMPolytopeTypeGetDim(ct));
1422: PetscStackCallStandard(ex_get_block,exoid, EX_ELEM_BLOCK, cs_id[cs], buffer, &num_cell_in_set,&num_vertex_per_cell, 0, 0, &num_attr);
1423: switch (ct) {
1424: case DM_POLYTOPE_TRI_PRISM:
1425: cs_order[cs] = cs;
1426: ++num_hybrid;
1427: break;
1428: default:
1429: for (c = cs; c > cs-num_hybrid; --c) cs_order[c] = cs_order[c-1];
1430: cs_order[cs-num_hybrid] = cs;
1431: }
1432: }
1433: /* First set sizes */
1434: for (ncs = 0, c = 0; ncs < num_cs; ++ncs) {
1435: DMPolytopeType ct;
1436: char elem_type[PETSC_MAX_PATH_LEN];
1437: const PetscInt cs = cs_order[ncs];
1439: PetscArrayzero(elem_type, sizeof(elem_type));
1440: PetscStackCallStandard(ex_get_elem_type,exoid, cs_id[cs], elem_type);
1441: ExodusGetCellType_Internal(elem_type, &ct);
1442: PetscStackCallStandard(ex_get_block,exoid, EX_ELEM_BLOCK, cs_id[cs], buffer, &num_cell_in_set,&num_vertex_per_cell, 0, 0, &num_attr);
1443: for (c_loc = 0; c_loc < num_cell_in_set; ++c_loc, ++c) {
1444: DMPlexSetConeSize(*dm, c, num_vertex_per_cell);
1445: DMPlexSetCellType(*dm, c, ct);
1446: }
1447: }
1448: for (v = numCells; v < numCells+numVertices; ++v) DMPlexSetCellType(*dm, v, DM_POLYTOPE_POINT);
1449: DMSetUp(*dm);
1450: for (ncs = 0, c = 0; ncs < num_cs; ++ncs) {
1451: const PetscInt cs = cs_order[ncs];
1452: PetscStackCallStandard(ex_get_block,exoid, EX_ELEM_BLOCK, cs_id[cs], buffer, &num_cell_in_set, &num_vertex_per_cell, 0, 0, &num_attr);
1453: PetscMalloc2(num_vertex_per_cell*num_cell_in_set,&cs_connect,num_vertex_per_cell,&cone);
1454: PetscStackCallStandard(ex_get_conn,exoid, EX_ELEM_BLOCK, cs_id[cs], cs_connect,NULL,NULL);
1455: /* EXO uses Fortran-based indexing, DMPlex uses C-style and numbers cell first then vertices. */
1456: for (c_loc = 0, v = 0; c_loc < num_cell_in_set; ++c_loc, ++c) {
1457: DMPolytopeType ct;
1459: for (v_loc = 0; v_loc < num_vertex_per_cell; ++v_loc, ++v) {
1460: cone[v_loc] = cs_connect[v]+numCells-1;
1461: }
1462: DMPlexGetCellType(*dm, c, &ct);
1463: DMPlexInvertCell(ct, cone);
1464: DMPlexSetCone(*dm, c, cone);
1465: DMSetLabelValue_Fast(*dm, &cellSets, "Cell Sets", c, cs_id[cs]);
1466: }
1467: PetscFree2(cs_connect,cone);
1468: }
1469: PetscFree2(cs_id, cs_order);
1470: }
1471: {
1472: PetscInt ints[] = {dim, dimEmbed};
1474: MPI_Bcast(ints, 2, MPIU_INT, 0, comm);
1475: DMSetDimension(*dm, ints[0]);
1476: DMSetCoordinateDim(*dm, ints[1]);
1477: dim = ints[0];
1478: dimEmbed = ints[1];
1479: }
1480: DMPlexSymmetrize(*dm);
1481: DMPlexStratify(*dm);
1482: if (interpolate) {
1483: DM idm;
1485: DMPlexInterpolate(*dm, &idm);
1486: DMDestroy(dm);
1487: *dm = idm;
1488: }
1490: /* Create vertex set label */
1491: if (rank == 0 && (num_vs > 0)) {
1492: int vs, v;
1493: /* Read from ex_get_node_set_ids() */
1494: int *vs_id;
1495: /* Read from ex_get_node_set_param() */
1496: int num_vertex_in_set;
1497: /* Read from ex_get_node_set() */
1498: int *vs_vertex_list;
1500: /* Get vertex set ids */
1501: PetscMalloc1(num_vs, &vs_id);
1502: PetscStackCallStandard(ex_get_ids,exoid, EX_NODE_SET, vs_id);
1503: for (vs = 0; vs < num_vs; ++vs) {
1504: PetscStackCallStandard(ex_get_set_param,exoid, EX_NODE_SET, vs_id[vs], &num_vertex_in_set, NULL);
1505: PetscMalloc1(num_vertex_in_set, &vs_vertex_list);
1506: PetscStackCallStandard(ex_get_set,exoid, EX_NODE_SET, vs_id[vs], vs_vertex_list, NULL);
1507: for (v = 0; v < num_vertex_in_set; ++v) {
1508: DMSetLabelValue_Fast(*dm, &vertSets, "Vertex Sets", vs_vertex_list[v]+numCells-1, vs_id[vs]);
1509: }
1510: PetscFree(vs_vertex_list);
1511: }
1512: PetscFree(vs_id);
1513: }
1514: /* Read coordinates */
1515: DMGetCoordinateSection(*dm, &coordSection);
1516: PetscSectionSetNumFields(coordSection, 1);
1517: PetscSectionSetFieldComponents(coordSection, 0, dimEmbed);
1518: PetscSectionSetChart(coordSection, numCells, numCells + numVertices);
1519: for (v = numCells; v < numCells+numVertices; ++v) {
1520: PetscSectionSetDof(coordSection, v, dimEmbed);
1521: PetscSectionSetFieldDof(coordSection, v, 0, dimEmbed);
1522: }
1523: PetscSectionSetUp(coordSection);
1524: PetscSectionGetStorageSize(coordSection, &coordSize);
1525: VecCreate(PETSC_COMM_SELF, &coordinates);
1526: PetscObjectSetName((PetscObject) coordinates, "coordinates");
1527: VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);
1528: VecSetBlockSize(coordinates, dimEmbed);
1529: VecSetType(coordinates,VECSTANDARD);
1530: VecGetArray(coordinates, &coords);
1531: if (rank == 0) {
1532: PetscReal *x, *y, *z;
1534: PetscMalloc3(numVertices,&x,numVertices,&y,numVertices,&z);
1535: PetscStackCallStandard(ex_get_coord,exoid, x, y, z);
1536: if (dimEmbed > 0) {
1537: for (v = 0; v < numVertices; ++v) coords[v*dimEmbed+0] = x[v];
1538: }
1539: if (dimEmbed > 1) {
1540: for (v = 0; v < numVertices; ++v) coords[v*dimEmbed+1] = y[v];
1541: }
1542: if (dimEmbed > 2) {
1543: for (v = 0; v < numVertices; ++v) coords[v*dimEmbed+2] = z[v];
1544: }
1545: PetscFree3(x,y,z);
1546: }
1547: VecRestoreArray(coordinates, &coords);
1548: DMSetCoordinatesLocal(*dm, coordinates);
1549: VecDestroy(&coordinates);
1551: /* Create side set label */
1552: if (rank == 0 && interpolate && (num_fs > 0)) {
1553: int fs, f, voff;
1554: /* Read from ex_get_side_set_ids() */
1555: int *fs_id;
1556: /* Read from ex_get_side_set_param() */
1557: int num_side_in_set;
1558: /* Read from ex_get_side_set_node_list() */
1559: int *fs_vertex_count_list, *fs_vertex_list;
1560: /* Read side set labels */
1561: char fs_name[MAX_STR_LENGTH+1];
1563: /* Get side set ids */
1564: PetscMalloc1(num_fs, &fs_id);
1565: PetscStackCallStandard(ex_get_ids,exoid, EX_SIDE_SET, fs_id);
1566: for (fs = 0; fs < num_fs; ++fs) {
1567: PetscStackCallStandard(ex_get_set_param,exoid, EX_SIDE_SET, fs_id[fs], &num_side_in_set, NULL);
1568: PetscMalloc2(num_side_in_set,&fs_vertex_count_list,num_side_in_set*4,&fs_vertex_list);
1569: PetscStackCallStandard(ex_get_side_set_node_list,exoid, fs_id[fs], fs_vertex_count_list, fs_vertex_list);
1570: /* Get the specific name associated with this side set ID. */
1571: int fs_name_err = ex_get_name(exoid, EX_SIDE_SET, fs_id[fs], fs_name);
1572: for (f = 0, voff = 0; f < num_side_in_set; ++f) {
1573: const PetscInt *faces = NULL;
1574: PetscInt faceSize = fs_vertex_count_list[f], numFaces;
1575: PetscInt faceVertices[4], v;
1578: for (v = 0; v < faceSize; ++v, ++voff) {
1579: faceVertices[v] = fs_vertex_list[voff]+numCells-1;
1580: }
1581: DMPlexGetFullJoin(*dm, faceSize, faceVertices, &numFaces, &faces);
1583: DMSetLabelValue_Fast(*dm, &faceSets, "Face Sets", faces[0], fs_id[fs]);
1584: /* Only add the label if one has been detected for this side set. */
1585: if (!fs_name_err) {
1586: DMSetLabelValue(*dm, fs_name, faces[0], fs_id[fs]);
1587: }
1588: DMPlexRestoreJoin(*dm, faceSize, faceVertices, &numFaces, &faces);
1589: }
1590: PetscFree2(fs_vertex_count_list,fs_vertex_list);
1591: }
1592: PetscFree(fs_id);
1593: }
1595: { /* Create Cell/Face/Vertex Sets labels at all processes */
1596: enum {n = 3};
1597: PetscBool flag[n];
1599: flag[0] = cellSets ? PETSC_TRUE : PETSC_FALSE;
1600: flag[1] = faceSets ? PETSC_TRUE : PETSC_FALSE;
1601: flag[2] = vertSets ? PETSC_TRUE : PETSC_FALSE;
1602: MPI_Bcast(flag, n, MPIU_BOOL, 0, comm);
1603: if (flag[0]) DMCreateLabel(*dm, "Cell Sets");
1604: if (flag[1]) DMCreateLabel(*dm, "Face Sets");
1605: if (flag[2]) DMCreateLabel(*dm, "Vertex Sets");
1606: }
1607: return 0;
1608: #else
1609: SETERRQ(comm, PETSC_ERR_SUP, "This method requires ExodusII support. Reconfigure using --download-exodusii");
1610: #endif
1611: }