Actual source code: plexexodusii.c
petsc-3.12.5 2020-03-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: #if defined(PETSC_HAVE_EXODUSII)
10: /*
11: EXOGetVarIndex - Locate a result in an exodus file based on its name
13: Not collective
15: Input Parameters:
16: + exoid - the exodus id of a file (obtained from ex_open or ex_create for instance)
17: . obj_type - the type of entity for instance EX_NODAL, EX_ELEM_BLOCK
18: - name - the name of the result
20: Output Parameters:
21: . varIndex - the location in the exodus file of the result
23: Notes:
24: The exodus variable index is obtained by comparing name and the
25: names of zonal variables declared in the exodus file. For instance if name is "V"
26: the location in the exodus file will be the first match of "V", "V_X", "V_XX", "V_1", or "V_11"
27: amongst all variables of type obj_type.
29: Level: beginner
31: .seealso: EXOGetVarIndex(),DMPlexView_ExodusII_Internal(),VecViewPlex_ExodusII_Nodal_Internal(),VecLoadNodal_PlexEXO(),VecLoadZonal_PlexEXO()
32: */
33: static PetscErrorCode EXOGetVarIndex_Private(int exoid, ex_entity_type obj_type, const char name[], int *varIndex)
34: {
35: int num_vars, i, j;
36: char ext_name[MAX_STR_LENGTH+1], var_name[MAX_STR_LENGTH+1];
37: const int num_suffix = 5;
38: char *suffix[5];
39: PetscBool flg;
43: suffix[0] = (char *) "";
44: suffix[1] = (char *) "_X";
45: suffix[2] = (char *) "_XX";
46: suffix[3] = (char *) "_1";
47: suffix[4] = (char *) "_11";
49: *varIndex = 0;
50: PetscStackCallStandard(ex_get_variable_param,(exoid, obj_type, &num_vars));
51: for (i = 0; i < num_vars; ++i) {
52: PetscStackCallStandard(ex_get_variable_name,(exoid, obj_type, i+1, var_name));
53: for (j = 0; j < num_suffix; ++j){
54: PetscStrncpy(ext_name, name, MAX_STR_LENGTH);
55: PetscStrlcat(ext_name, suffix[j], MAX_STR_LENGTH);
56: PetscStrcasecmp(ext_name, var_name, &flg);
57: if (flg) {
58: *varIndex = i+1;
59: return(0);
60: }
61: }
62: }
63: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unable to locate variable %s in ExodusII file.", name);
64: PetscFunctionReturn(-1);
65: }
67: /*
68: DMPlexView_ExodusII_Internal - Write a DM to disk in exodus format
70: Collective on dm
72: Input Parameters:
73: + dm - The dm to be written
74: . exoid - the exodus id of a file (obtained from ex_open or ex_create for instance)
75: - degree - the degree of the interpolation space
77: Notes:
78: Not all DM can be written to disk this way. For instance, exodus assume that element blocks (mapped to "Cell sets" labels)
79: consists of sequentially numbered cells. If this is not the case, the exodus file will be corrupted.
81: If the dm has been distributed and exoid points to different files on each MPI rank, only the "local" part of the DM
82: (including "ghost" cells and vertices) will be written. If all exoid points to the same file, the resulting file will
83: probably be corrupted.
85: DMPlex only represents geometry while most post-processing software expect that a mesh also provides information
86: on the discretization space. This function assumes that the file represents Lagrange finite elements of order 1 or 2.
87: It should be extended to use PetscFE objects.
89: This function will only handle TRI, TET, QUAD and HEX cells.
90: Level: beginner
92: .seealso: EXOGetVarIndex_Private(), VecViewPlex_ExodusII_Nodal_Internal(), VecLoadNodal_PlexEXO(), VecViewZonal_PlexEXO(), VecLoadZonal_PlexEXO()
93: */
94: PetscErrorCode DMPlexView_ExodusII_Internal(DM dm, int exoid, PetscInt degree)
95: {
96: enum ElemType {TRI, QUAD, TET, HEX};
97: MPI_Comm comm;
98: /* Connectivity Variables */
99: PetscInt cellsNotInConnectivity;
100: /* Cell Sets */
101: DMLabel csLabel;
102: IS csIS;
103: const PetscInt *csIdx;
104: PetscInt num_cs, cs;
105: enum ElemType *type;
106: PetscBool hasLabel;
107: /* Coordinate Variables */
108: DM cdm;
109: PetscSection section;
110: Vec coord;
111: PetscInt **nodes;
112: PetscInt depth, d, dim, skipCells = 0;
113: PetscInt pStart, pEnd, p, cStart, cEnd, numCells, vStart, vEnd, numVertices, eStart, eEnd, numEdges, fStart, fEnd, numFaces, numNodes;
114: PetscInt num_vs, num_fs;
115: PetscMPIInt rank, size;
116: const char *dmName;
117: PetscInt nodesTriP1[4] = {3,0,0,0};
118: PetscInt nodesTriP2[4] = {3,3,0,0};
119: PetscInt nodesQuadP1[4] = {4,0,0,0};
120: PetscInt nodesQuadP2[4] = {4,4,0,1};
121: PetscInt nodesTetP1[4] = {4,0,0,0};
122: PetscInt nodesTetP2[4] = {4,6,0,0};
123: PetscInt nodesHexP1[4] = {8,0,0,0};
124: PetscInt nodesHexP2[4] = {8,12,6,1};
125: PetscErrorCode ierr;
128: PetscObjectGetComm((PetscObject)dm, &comm);
129: MPI_Comm_rank(comm, &rank);
130: MPI_Comm_size(comm, &size);
131: /* --- Get DM info --- */
132: PetscObjectGetName((PetscObject) dm, &dmName);
133: DMPlexGetDepth(dm, &depth);
134: DMGetDimension(dm, &dim);
135: DMPlexGetChart(dm, &pStart, &pEnd);
136: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
137: DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);
138: DMPlexGetDepthStratum(dm, 1, &eStart, &eEnd);
139: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
140: numCells = cEnd - cStart;
141: numEdges = eEnd - eStart;
142: numVertices = vEnd - vStart;
143: if (depth == 3) {numFaces = fEnd - fStart;}
144: else {numFaces = 0;}
145: DMGetLabelSize(dm, "Cell Sets", &num_cs);
146: DMGetLabelSize(dm, "Vertex Sets", &num_vs);
147: DMGetLabelSize(dm, "Face Sets", &num_fs);
148: DMGetCoordinatesLocal(dm, &coord);
149: DMGetCoordinateDM(dm, &cdm);
150: if (num_cs > 0) {
151: DMGetLabel(dm, "Cell Sets", &csLabel);
152: DMLabelGetValueIS(csLabel, &csIS);
153: ISGetIndices(csIS, &csIdx);
154: }
155: PetscMalloc1(num_cs, &nodes);
156: /* Set element type for each block and compute total number of nodes */
157: PetscMalloc1(num_cs, &type);
158: numNodes = numVertices;
159: if (degree == 2) {numNodes += numEdges;}
160: cellsNotInConnectivity = numCells;
161: for (cs = 0; cs < num_cs; ++cs) {
162: IS stratumIS;
163: const PetscInt *cells;
164: PetscScalar *xyz = NULL;
165: PetscInt csSize, closureSize;
167: DMLabelGetStratumIS(csLabel, csIdx[cs], &stratumIS);
168: ISGetIndices(stratumIS, &cells);
169: ISGetSize(stratumIS, &csSize);
170: DMPlexVecGetClosure(cdm, NULL, coord, cells[0], &closureSize, &xyz);
171: switch (dim) {
172: case 2:
173: if (closureSize == 3*dim) {type[cs] = TRI;}
174: else if (closureSize == 4*dim) {type[cs] = QUAD;}
175: else SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Number of vertices %D in dimension %D has no ExodusII type", closureSize/dim, dim);
176: break;
177: case 3:
178: if (closureSize == 4*dim) {type[cs] = TET;}
179: else if (closureSize == 8*dim) {type[cs] = HEX;}
180: else SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Number of vertices %D in dimension %D has no ExodusII type", closureSize/dim, dim);
181: break;
182: default: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Dimension %D not handled by ExodusII viewer", dim);
183: }
184: if ((degree == 2) && (type[cs] == QUAD)) {numNodes += csSize;}
185: if ((degree == 2) && (type[cs] == HEX)) {numNodes += csSize; numNodes += numFaces;}
186: DMPlexVecRestoreClosure(cdm, NULL, coord, cells[0], &closureSize, &xyz);
187: /* Set nodes and Element type */
188: if (type[cs] == TRI) {
189: if (degree == 1) nodes[cs] = nodesTriP1;
190: else if (degree == 2) nodes[cs] = nodesTriP2;
191: } else if (type[cs] == QUAD) {
192: if (degree == 1) nodes[cs] = nodesQuadP1;
193: else if (degree == 2) nodes[cs] = nodesQuadP2;
194: } else if (type[cs] == TET) {
195: if (degree == 1) nodes[cs] = nodesTetP1;
196: else if (degree == 2) nodes[cs] = nodesTetP2;
197: } else if (type[cs] == HEX) {
198: if (degree == 1) nodes[cs] = nodesHexP1;
199: else if (degree == 2) nodes[cs] = nodesHexP2;
200: }
201: /* Compute the number of cells not in the connectivity table */
202: cellsNotInConnectivity -= nodes[cs][3]*csSize;
204: ISRestoreIndices(stratumIS, &cells);
205: ISDestroy(&stratumIS);
206: }
207: if (num_cs > 0) {PetscStackCallStandard(ex_put_init,(exoid, dmName, dim, numNodes, numCells, num_cs, num_vs, num_fs));}
208: /* --- Connectivity --- */
209: for (cs = 0; cs < num_cs; ++cs) {
210: IS stratumIS;
211: const PetscInt *cells;
212: PetscInt *connect, off = 0;
213: PetscInt edgesInClosure = 0, facesInClosure = 0, verticesInClosure = 0;
214: PetscInt csSize, c, connectSize, closureSize;
215: char *elem_type = NULL;
216: char elem_type_tri3[] = "TRI3", elem_type_quad4[] = "QUAD4";
217: char elem_type_tri6[] = "TRI6", elem_type_quad9[] = "QUAD9";
218: char elem_type_tet4[] = "TET4", elem_type_hex8[] = "HEX8";
219: char elem_type_tet10[] = "TET10", elem_type_hex27[] = "HEX27";
221: DMLabelGetStratumIS(csLabel, csIdx[cs], &stratumIS);
222: ISGetIndices(stratumIS, &cells);
223: ISGetSize(stratumIS, &csSize);
224: /* Set Element type */
225: if (type[cs] == TRI) {
226: if (degree == 1) elem_type = elem_type_tri3;
227: else if (degree == 2) elem_type = elem_type_tri6;
228: } else if (type[cs] == QUAD) {
229: if (degree == 1) elem_type = elem_type_quad4;
230: else if (degree == 2) elem_type = elem_type_quad9;
231: } else if (type[cs] == TET) {
232: if (degree == 1) elem_type = elem_type_tet4;
233: else if (degree == 2) elem_type = elem_type_tet10;
234: } else if (type[cs] == HEX) {
235: if (degree == 1) elem_type = elem_type_hex8;
236: else if (degree == 2) elem_type = elem_type_hex27;
237: }
238: connectSize = nodes[cs][0] + nodes[cs][1] + nodes[cs][2] + nodes[cs][3];
239: PetscMalloc1(PetscMax(27,connectSize)*csSize, &connect);
240: PetscStackCallStandard(ex_put_block,(exoid, EX_ELEM_BLOCK, csIdx[cs], elem_type, csSize, connectSize, 0, 0, 1));
241: /* Find number of vertices, edges, and faces in the closure */
242: verticesInClosure = nodes[cs][0];
243: if (depth > 1) {
244: if (dim == 2) {
245: DMPlexGetConeSize(dm, cells[0], &edgesInClosure);
246: } else if (dim == 3) {
247: PetscInt *closure = NULL;
249: DMPlexGetConeSize(dm, cells[0], &facesInClosure);
250: DMPlexGetTransitiveClosure(dm, cells[0], PETSC_TRUE, &closureSize, &closure);
251: edgesInClosure = closureSize - facesInClosure - 1 - verticesInClosure;
252: DMPlexRestoreTransitiveClosure(dm, cells[0], PETSC_TRUE, &closureSize, &closure);
253: }
254: }
255: /* Get connectivity for each cell */
256: for (c = 0; c < csSize; ++c) {
257: PetscInt *closure = NULL;
258: PetscInt temp, i;
260: DMPlexGetTransitiveClosure(dm, cells[c], PETSC_TRUE, &closureSize, &closure);
261: for (i = 0; i < connectSize; ++i) {
262: if (i < nodes[cs][0]) {/* Vertices */
263: connect[i+off] = closure[(i+edgesInClosure+facesInClosure+1)*2] + 1;
264: connect[i+off] -= cellsNotInConnectivity;
265: } else if (i < nodes[cs][0]+nodes[cs][1]) { /* Edges */
266: connect[i+off] = closure[(i-verticesInClosure+facesInClosure+1)*2] + 1;
267: if (nodes[cs][2] == 0) connect[i+off] -= numFaces;
268: connect[i+off] -= cellsNotInConnectivity;
269: } else if (i < nodes[cs][0]+nodes[cs][1]+nodes[cs][3]){ /* Cells */
270: connect[i+off] = closure[0] + 1;
271: connect[i+off] -= skipCells;
272: } else if (i < nodes[cs][0]+nodes[cs][1]+nodes[cs][3]+nodes[cs][2]){ /* Faces */
273: connect[i+off] = closure[(i-edgesInClosure-verticesInClosure)*2] + 1;
274: connect[i+off] -= cellsNotInConnectivity;
275: } else {
276: connect[i+off] = -1;
277: }
278: }
279: /* Tetrahedra are inverted */
280: if (type[cs] == TET) {
281: temp = connect[0+off]; connect[0+off] = connect[1+off]; connect[1+off] = temp;
282: if (degree == 2) {
283: temp = connect[5+off]; connect[5+off] = connect[6+off]; connect[6+off] = temp;
284: temp = connect[7+off]; connect[7+off] = connect[8+off]; connect[8+off] = temp;
285: }
286: }
287: /* Hexahedra are inverted */
288: if (type[cs] == HEX) {
289: temp = connect[1+off]; connect[1+off] = connect[3+off]; connect[3+off] = temp;
290: if (degree == 2) {
291: temp = connect[8+off]; connect[8+off] = connect[11+off]; connect[11+off] = temp;
292: temp = connect[9+off]; connect[9+off] = connect[10+off]; connect[10+off] = temp;
293: temp = connect[16+off]; connect[16+off] = connect[17+off]; connect[17+off] = temp;
294: temp = connect[18+off]; connect[18+off] = connect[19+off]; connect[19+off] = temp;
296: temp = connect[12+off]; connect[12+off] = connect[16+off]; connect[16+off] = temp;
297: temp = connect[13+off]; connect[13+off] = connect[17+off]; connect[17+off] = temp;
298: temp = connect[14+off]; connect[14+off] = connect[18+off]; connect[18+off] = temp;
299: temp = connect[15+off]; connect[15+off] = connect[19+off]; connect[19+off] = temp;
301: temp = connect[23+off]; connect[23+off] = connect[26+off]; connect[26+off] = temp;
302: temp = connect[24+off]; connect[24+off] = connect[25+off]; connect[25+off] = temp;
303: temp = connect[25+off]; connect[25+off] = connect[26+off]; connect[26+off] = temp;
304: }
305: }
306: off += connectSize;
307: DMPlexRestoreTransitiveClosure(dm, cells[c], PETSC_TRUE, &closureSize, &closure);
308: }
309: PetscStackCallStandard(ex_put_conn,(exoid, EX_ELEM_BLOCK, csIdx[cs], connect, 0, 0));
310: skipCells += (nodes[cs][3] == 0)*csSize;
311: PetscFree(connect);
312: ISRestoreIndices(stratumIS, &cells);
313: ISDestroy(&stratumIS);
314: }
315: PetscFree(type);
316: /* --- Coordinates --- */
317: PetscSectionCreate(comm, §ion);
318: PetscSectionSetChart(section, pStart, pEnd);
319: for (d = 0; d < depth; ++d) {
320: DMPlexGetDepthStratum(dm, d, &pStart, &pEnd);
321: for (p = pStart; p < pEnd; ++p) {
322: PetscSectionSetDof(section, p, nodes[0][d] > 0);
323: }
324: }
325: for (cs = 0; cs < num_cs; ++cs) {
326: IS stratumIS;
327: const PetscInt *cells;
328: PetscInt csSize, c;
330: DMLabelGetStratumIS(csLabel, csIdx[cs], &stratumIS);
331: ISGetIndices(stratumIS, &cells);
332: ISGetSize(stratumIS, &csSize);
333: for (c = 0; c < csSize; ++c) {
334: PetscSectionSetDof(section, cells[c], nodes[cs][3] > 0);
335: }
336: ISRestoreIndices(stratumIS, &cells);
337: ISDestroy(&stratumIS);
338: }
339: if (num_cs > 0) {
340: ISRestoreIndices(csIS, &csIdx);
341: ISDestroy(&csIS);
342: }
343: PetscFree(nodes);
344: PetscSectionSetUp(section);
345: if (numNodes > 0) {
346: const char *coordNames[3] = {"x", "y", "z"};
347: PetscScalar *coords, *closure;
348: PetscReal *cval;
349: PetscInt hasDof, n = 0;
351: /* There can't be more than 24 values in the closure of a point for the coord section */
352: PetscCalloc3(numNodes*3, &coords, dim, &cval, 24, &closure);
353: DMGetCoordinatesLocal(dm, &coord);
354: DMPlexGetChart(dm, &pStart, &pEnd);
355: for (p = pStart; p < pEnd; ++p) {
356: PetscSectionGetDof(section, p, &hasDof);
357: if (hasDof) {
358: PetscInt closureSize = 24, j;
360: DMPlexVecGetClosure(cdm, NULL, coord, p, &closureSize, &closure);
361: for (d = 0; d < dim; ++d) {
362: cval[d] = 0.0;
363: for (j = 0; j < closureSize/dim; j++) cval[d] += closure[j*dim+d];
364: coords[d*numNodes+n] = cval[d] * dim / closureSize;
365: }
366: ++n;
367: }
368: }
369: PetscStackCallStandard(ex_put_coord,(exoid, &coords[0*numNodes], &coords[1*numNodes], &coords[2*numNodes]));
370: PetscFree3(coords, cval, closure);
371: PetscStackCallStandard(ex_put_coord_names,(exoid, (char **) coordNames));
372: }
373: PetscSectionDestroy(§ion);
374: /* --- Node Sets/Vertex Sets --- */
375: DMHasLabel(dm, "Vertex Sets", &hasLabel);
376: if (hasLabel) {
377: PetscInt i, vs, vsSize;
378: const PetscInt *vsIdx, *vertices;
379: PetscInt *nodeList;
380: IS vsIS, stratumIS;
381: DMLabel vsLabel;
382: DMGetLabel(dm, "Vertex Sets", &vsLabel);
383: DMLabelGetValueIS(vsLabel, &vsIS);
384: ISGetIndices(vsIS, &vsIdx);
385: for (vs=0; vs<num_vs; ++vs) {
386: DMLabelGetStratumIS(vsLabel, vsIdx[vs], &stratumIS);
387: ISGetIndices(stratumIS, &vertices);
388: ISGetSize(stratumIS, &vsSize);
389: PetscMalloc1(vsSize, &nodeList);
390: for (i=0; i<vsSize; ++i) {
391: nodeList[i] = vertices[i] - skipCells + 1;
392: }
393: PetscStackCallStandard(ex_put_set_param,(exoid, EX_NODE_SET, vsIdx[vs], vsSize, 0));
394: PetscStackCallStandard(ex_put_set,(exoid, EX_NODE_SET, vsIdx[vs], nodeList, NULL));
395: ISRestoreIndices(stratumIS, &vertices);
396: ISDestroy(&stratumIS);
397: PetscFree(nodeList);
398: }
399: ISRestoreIndices(vsIS, &vsIdx);
400: ISDestroy(&vsIS);
401: }
402: /* --- Side Sets/Face Sets --- */
403: DMHasLabel(dm, "Face Sets", &hasLabel);
404: if (hasLabel) {
405: PetscInt i, j, fs, fsSize;
406: const PetscInt *fsIdx, *faces;
407: IS fsIS, stratumIS;
408: DMLabel fsLabel;
409: PetscInt numPoints, *points;
410: PetscInt elem_list_size = 0;
411: PetscInt *elem_list, *elem_ind, *side_list;
413: DMGetLabel(dm, "Face Sets", &fsLabel);
414: /* Compute size of Node List and Element List */
415: DMLabelGetValueIS(fsLabel, &fsIS);
416: ISGetIndices(fsIS, &fsIdx);
417: for (fs=0; fs<num_fs; ++fs) {
418: DMLabelGetStratumIS(fsLabel, fsIdx[fs], &stratumIS);
419: ISGetSize(stratumIS, &fsSize);
420: elem_list_size += fsSize;
421: ISDestroy(&stratumIS);
422: }
423: PetscMalloc3(num_fs, &elem_ind,
424: elem_list_size, &elem_list,
425: elem_list_size, &side_list);
426: elem_ind[0] = 0;
427: for (fs=0; fs<num_fs; ++fs) {
428: DMLabelGetStratumIS(fsLabel, fsIdx[fs], &stratumIS);
429: ISGetIndices(stratumIS, &faces);
430: ISGetSize(stratumIS, &fsSize);
431: /* Set Parameters */
432: PetscStackCallStandard(ex_put_set_param,(exoid, EX_SIDE_SET, fsIdx[fs], fsSize, 0));
433: /* Indices */
434: if (fs<num_fs-1) {
435: elem_ind[fs+1] = elem_ind[fs] + fsSize;
436: }
438: for (i=0; i<fsSize; ++i) {
439: /* Element List */
440: points = NULL;
441: DMPlexGetTransitiveClosure(dm, faces[i], PETSC_FALSE, &numPoints, &points);
442: elem_list[elem_ind[fs] + i] = points[2] +1;
443: DMPlexRestoreTransitiveClosure(dm, faces[i], PETSC_FALSE, &numPoints, &points);
445: /* Side List */
446: points = NULL;
447: DMPlexGetTransitiveClosure(dm, elem_list[elem_ind[fs] + i]-1, PETSC_TRUE, &numPoints, &points);
448: for (j=1; j<numPoints; ++j) {
449: if (points[j*2]==faces[i]) {break;}
450: }
451: /* Convert HEX sides */
452: if (numPoints == 27) {
453: if (j == 1) {j = 5;}
454: else if (j == 2) {j = 6;}
455: else if (j == 3) {j = 1;}
456: else if (j == 4) {j = 3;}
457: else if (j == 5) {j = 2;}
458: else if (j == 6) {j = 4;}
459: }
460: /* Convert TET sides */
461: if (numPoints == 15) {
462: --j;
463: if (j == 0) {j = 4;}
464: }
465: side_list[elem_ind[fs] + i] = j;
466: DMPlexRestoreTransitiveClosure(dm, elem_list[elem_ind[fs] + i]-1, PETSC_TRUE, &numPoints, &points);
468: }
469: ISRestoreIndices(stratumIS, &faces);
470: ISDestroy(&stratumIS);
471: }
472: ISRestoreIndices(fsIS, &fsIdx);
473: ISDestroy(&fsIS);
475: /* Put side sets */
476: for (fs=0; fs<num_fs; ++fs) {
477: PetscStackCallStandard(ex_put_set,(exoid, EX_SIDE_SET, fsIdx[fs], &elem_list[elem_ind[fs]], &side_list[elem_ind[fs]]));
478: }
479: PetscFree3(elem_ind,elem_list,side_list);
480: }
481: return(0);
482: }
484: /*
485: VecViewPlex_ExodusII_Nodal_Internal - Write a Vec corresponding to a nodal field to an exodus file
487: Collective on v
489: Input Parameters:
490: + v - The vector to be written
491: . exoid - the exodus id of a file (obtained from ex_open or ex_create for instance)
492: - step - the time step to write at (exodus steps are numbered starting from 1)
494: Notes:
495: The exodus result nodal variable index is obtained by comparing the Vec name and the
496: names of nodal variables declared in the exodus file. For instance for a Vec named "V"
497: the location in the exodus file will be the first match of "V", "V_X", "V_XX", "V_1", or "V_11"
498: amongst all nodal variables.
500: Level: beginner
502: .seealso: EXOGetVarIndex_Private(),DMPlexView_ExodusII_Internal(),VecLoadNodal_PlexEXO(),VecViewZonal_PlexEXO(),VecLoadZonal_PlexEXO()
503: @*/
504: PetscErrorCode VecViewPlex_ExodusII_Nodal_Internal(Vec v, int exoid, int step)
505: {
506: MPI_Comm comm;
507: PetscMPIInt size;
508: DM dm;
509: Vec vNatural, vComp;
510: const PetscScalar *varray;
511: const char *vecname;
512: PetscInt xs, xe, bs;
513: PetscBool useNatural;
514: int offset;
515: PetscErrorCode ierr;
518: PetscObjectGetComm((PetscObject) v, &comm);
519: MPI_Comm_size(comm, &size);
520: VecGetDM(v, &dm);
521: DMGetUseNatural(dm, &useNatural);
522: useNatural = useNatural && size > 1 ? PETSC_TRUE : PETSC_FALSE;
523: if (useNatural) {
524: DMGetGlobalVector(dm, &vNatural);
525: DMPlexGlobalToNaturalBegin(dm, v, vNatural);
526: DMPlexGlobalToNaturalEnd(dm, v, vNatural);
527: } else {
528: vNatural = v;
529: }
530: /* Get the location of the variable in the exodus file */
531: PetscObjectGetName((PetscObject) v, &vecname);
532: EXOGetVarIndex_Private(exoid, EX_NODAL, vecname, &offset);
533: if (!offset) SETERRQ1(comm, PETSC_ERR_FILE_UNEXPECTED, "Cannot locate nodal variable %s in exodus file\n", vecname);
534: /* Write local chunk of the result in the exodus file
535: exodus stores each component of a vector-valued field as a separate variable.
536: We assume that they are stored sequentially */
537: VecGetOwnershipRange(vNatural, &xs, &xe);
538: VecGetBlockSize(vNatural, &bs);
539: if (bs == 1) {
540: VecGetArrayRead(vNatural, &varray);
541: PetscStackCallStandard(ex_put_partial_var,(exoid, step, EX_NODAL, offset, 1, xs+1, xe-xs, varray));
542: VecRestoreArrayRead(vNatural, &varray);
543: } else {
544: IS compIS;
545: PetscInt c;
547: ISCreateStride(PETSC_COMM_SELF, (xe-xs)/bs, xs, bs, &compIS);
548: for (c = 0; c < bs; ++c) {
549: ISStrideSetStride(compIS, (xe-xs)/bs, xs+c, bs);
550: VecGetSubVector(vNatural, compIS, &vComp);
551: VecGetArrayRead(vComp, &varray);
552: PetscStackCallStandard(ex_put_partial_var,(exoid, step, EX_NODAL, offset+c, 1, xs/bs+1, (xe-xs)/bs, varray));
553: VecRestoreArrayRead(vComp, &varray);
554: VecRestoreSubVector(vNatural, compIS, &vComp);
555: }
556: ISDestroy(&compIS);
557: }
558: if (useNatural) {DMRestoreGlobalVector(dm, &vNatural);}
559: return(0);
560: }
562: /*
563: VecLoadPlex_ExodusII_Nodal_Internal - Read a Vec corresponding to a nodal field from an exodus file
565: Collective on v
567: Input Parameters:
568: + v - The vector to be written
569: . exoid - the exodus id of a file (obtained from ex_open or ex_create for instance)
570: - step - the time step to read at (exodus steps are numbered starting from 1)
572: Notes:
573: The exodus result nodal variable index is obtained by comparing the Vec name and the
574: names of nodal variables declared in the exodus file. For instance for a Vec named "V"
575: the location in the exodus file will be the first match of "V", "V_X", "V_XX", "V_1", or "V_11"
576: amongst all nodal variables.
578: Level: beginner
580: .seealso: EXOGetVarIndex_Private(), DMPlexView_ExodusII_Internal(), VecViewPlex_ExodusII_Nodal_Internal(), VecViewZonal_PlexEXO(), VecLoadZonal_PlexEXO()
581: */
582: PetscErrorCode VecLoadPlex_ExodusII_Nodal_Internal(Vec v, int exoid, int step)
583: {
584: MPI_Comm comm;
585: PetscMPIInt size;
586: DM dm;
587: Vec vNatural, vComp;
588: PetscScalar *varray;
589: const char *vecname;
590: PetscInt xs, xe, bs;
591: PetscBool useNatural;
592: int offset;
596: PetscObjectGetComm((PetscObject) v, &comm);
597: MPI_Comm_size(comm, &size);
598: VecGetDM(v,&dm);
599: DMGetUseNatural(dm, &useNatural);
600: useNatural = useNatural && size > 1 ? PETSC_TRUE : PETSC_FALSE;
601: if (useNatural) {DMGetGlobalVector(dm,&vNatural);}
602: else {vNatural = v;}
603: /* Get the location of the variable in the exodus file */
604: PetscObjectGetName((PetscObject) v, &vecname);
605: EXOGetVarIndex_Private(exoid, EX_NODAL, vecname, &offset);
606: if (!offset) SETERRQ1(comm, PETSC_ERR_FILE_UNEXPECTED, "Cannot locate nodal variable %s in exodus file", vecname);
607: /* Read local chunk from the file */
608: VecGetOwnershipRange(vNatural, &xs, &xe);
609: VecGetBlockSize(vNatural, &bs);
610: if (bs == 1) {
611: VecGetArray(vNatural, &varray);
612: PetscStackCallStandard(ex_get_partial_var,(exoid, step, EX_NODAL, offset, 1, xs+1, xe-xs, varray));
613: VecRestoreArray(vNatural, &varray);
614: } else {
615: IS compIS;
616: PetscInt c;
618: ISCreateStride(PETSC_COMM_SELF, (xe-xs)/bs, xs, bs, &compIS);
619: for (c = 0; c < bs; ++c) {
620: ISStrideSetStride(compIS, (xe-xs)/bs, xs+c, bs);
621: VecGetSubVector(vNatural, compIS, &vComp);
622: VecGetArray(vComp, &varray);
623: PetscStackCallStandard(ex_get_partial_var,(exoid, step, EX_NODAL, offset+c, 1, xs/bs+1, (xe-xs)/bs, varray));
624: VecRestoreArray(vComp, &varray);
625: VecRestoreSubVector(vNatural, compIS, &vComp);
626: }
627: ISDestroy(&compIS);
628: }
629: if (useNatural) {
630: DMPlexNaturalToGlobalBegin(dm, vNatural, v);
631: DMPlexNaturalToGlobalEnd(dm, vNatural, v);
632: DMRestoreGlobalVector(dm, &vNatural);
633: }
634: return(0);
635: }
637: /*
638: VecViewPlex_ExodusII_Zonal_Internal - Write a Vec corresponding to a zonal (cell based) field to an exodus file
640: Collective on v
642: Input Parameters:
643: + v - The vector to be written
644: . exoid - the exodus id of a file (obtained from ex_open or ex_create for instance)
645: - step - the time step to write at (exodus steps are numbered starting from 1)
647: Notes:
648: The exodus result zonal variable index is obtained by comparing the Vec name and the
649: names of zonal variables declared in the exodus file. For instance for a Vec named "V"
650: the location in the exodus file will be the first match of "V", "V_X", "V_XX", "V_1", or "V_11"
651: amongst all zonal variables.
653: Level: beginner
655: .seealso: EXOGetVarIndex_Private(),DMPlexView_ExodusII_Internal(),VecViewPlex_ExodusII_Nodal_Internal(),VecLoadPlex_ExodusII_Nodal_Internal(),VecLoadPlex_ExodusII_Zonal_Internal()
656: */
657: PetscErrorCode VecViewPlex_ExodusII_Zonal_Internal(Vec v, int exoid, int step)
658: {
659: MPI_Comm comm;
660: PetscMPIInt size;
661: DM dm;
662: Vec vNatural, vComp;
663: const PetscScalar *varray;
664: const char *vecname;
665: PetscInt xs, xe, bs;
666: PetscBool useNatural;
667: int offset;
668: IS compIS;
669: PetscInt *csSize, *csID;
670: PetscInt numCS, set, csxs = 0;
671: PetscErrorCode ierr;
674: PetscObjectGetComm((PetscObject)v, &comm);
675: MPI_Comm_size(comm, &size);
676: VecGetDM(v, &dm);
677: DMGetUseNatural(dm, &useNatural);
678: useNatural = useNatural && size > 1 ? PETSC_TRUE : PETSC_FALSE;
679: if (useNatural) {
680: DMGetGlobalVector(dm, &vNatural);
681: DMPlexGlobalToNaturalBegin(dm, v, vNatural);
682: DMPlexGlobalToNaturalEnd(dm, v, vNatural);
683: } else {
684: vNatural = v;
685: }
686: /* Get the location of the variable in the exodus file */
687: PetscObjectGetName((PetscObject) v, &vecname);
688: EXOGetVarIndex_Private(exoid, EX_ELEM_BLOCK, vecname, &offset);
689: if (!offset) SETERRQ1(comm, PETSC_ERR_FILE_UNEXPECTED, "Cannot locate zonal variable %s in exodus file\n", vecname);
690: /* Write local chunk of the result in the exodus file
691: exodus stores each component of a vector-valued field as a separate variable.
692: We assume that they are stored sequentially
693: Zonal variables are accessed one element block at a time, so we loop through the cell sets,
694: but once the vector has been reordered to natural size, we cannot use the label informations
695: to figure out what to save where. */
696: numCS = ex_inquire_int(exoid, EX_INQ_ELEM_BLK);
697: PetscMalloc2(numCS, &csID, numCS, &csSize);
698: PetscStackCallStandard(ex_get_ids,(exoid, EX_ELEM_BLOCK, csID));
699: for (set = 0; set < numCS; ++set) {
700: ex_block block;
702: block.id = csID[set];
703: block.type = EX_ELEM_BLOCK;
704: PetscStackCallStandard(ex_get_block_param,(exoid, &block));
705: csSize[set] = block.num_entry;
706: }
707: VecGetOwnershipRange(vNatural, &xs, &xe);
708: VecGetBlockSize(vNatural, &bs);
709: if (bs > 1) {ISCreateStride(comm, (xe-xs)/bs, xs, bs, &compIS);}
710: for (set = 0; set < numCS; set++) {
711: PetscInt csLocalSize, c;
713: /* range of indices for set setID[set]: csxs:csxs + csSize[set]-1
714: local slice of zonal values: xs/bs,xm/bs-1
715: intersection: max(xs/bs,csxs),min(xm/bs-1,csxs + csSize[set]-1) */
716: csLocalSize = PetscMax(0, PetscMin(xe/bs, csxs+csSize[set]) - PetscMax(xs/bs, csxs));
717: if (bs == 1) {
718: VecGetArrayRead(vNatural, &varray);
719: PetscStackCallStandard(ex_put_partial_var,(exoid, step, EX_ELEM_BLOCK, offset, csID[set], PetscMax(xs-csxs, 0)+1, csLocalSize, &varray[PetscMax(0, csxs-xs)]));
720: VecRestoreArrayRead(vNatural, &varray);
721: } else {
722: for (c = 0; c < bs; ++c) {
723: ISStrideSetStride(compIS, (xe-xs)/bs, xs+c, bs);
724: VecGetSubVector(vNatural, compIS, &vComp);
725: VecGetArrayRead(vComp, &varray);
726: 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)]));
727: VecRestoreArrayRead(vComp, &varray);
728: VecRestoreSubVector(vNatural, compIS, &vComp);
729: }
730: }
731: csxs += csSize[set];
732: }
733: PetscFree2(csID, csSize);
734: if (bs > 1) {ISDestroy(&compIS);}
735: if (useNatural) {DMRestoreGlobalVector(dm,&vNatural);}
736: return(0);
737: }
739: /*
740: VecLoadPlex_ExodusII_Zonal_Internal - Read a Vec corresponding to a zonal (cell based) field from an exodus file
742: Collective on v
744: Input Parameters:
745: + v - The vector to be written
746: . exoid - the exodus id of a file (obtained from ex_open or ex_create for instance)
747: - step - the time step to read at (exodus steps are numbered starting from 1)
749: Notes:
750: The exodus result zonal variable index is obtained by comparing the Vec name and the
751: names of zonal variables declared in the exodus file. For instance for a Vec named "V"
752: the location in the exodus file will be the first match of "V", "V_X", "V_XX", "V_1", or "V_11"
753: amongst all zonal variables.
755: Level: beginner
757: .seealso: EXOGetVarIndex_Private(), DMPlexView_ExodusII_Internal(), VecViewPlex_ExodusII_Nodal_Internal(), VecLoadPlex_ExodusII_Nodal_Internal(), VecLoadPlex_ExodusII_Zonal_Internal()
758: */
759: PetscErrorCode VecLoadPlex_ExodusII_Zonal_Internal(Vec v, int exoid, int step)
760: {
761: MPI_Comm comm;
762: PetscMPIInt size;
763: DM dm;
764: Vec vNatural, vComp;
765: PetscScalar *varray;
766: const char *vecname;
767: PetscInt xs, xe, bs;
768: PetscBool useNatural;
769: int offset;
770: IS compIS;
771: PetscInt *csSize, *csID;
772: PetscInt numCS, set, csxs = 0;
773: PetscErrorCode ierr;
776: PetscObjectGetComm((PetscObject)v,&comm);
777: MPI_Comm_size(comm, &size);
778: VecGetDM(v, &dm);
779: DMGetUseNatural(dm, &useNatural);
780: useNatural = useNatural && size > 1 ? PETSC_TRUE : PETSC_FALSE;
781: if (useNatural) {DMGetGlobalVector(dm,&vNatural);}
782: else {vNatural = v;}
783: /* Get the location of the variable in the exodus file */
784: PetscObjectGetName((PetscObject) v, &vecname);
785: EXOGetVarIndex_Private(exoid, EX_ELEM_BLOCK, vecname, &offset);
786: if (!offset) SETERRQ1(comm, PETSC_ERR_FILE_UNEXPECTED, "Cannot locate zonal variable %s in exodus file\n", vecname);
787: /* Read local chunk of the result in the exodus file
788: exodus stores each component of a vector-valued field as a separate variable.
789: We assume that they are stored sequentially
790: Zonal variables are accessed one element block at a time, so we loop through the cell sets,
791: but once the vector has been reordered to natural size, we cannot use the label informations
792: to figure out what to save where. */
793: numCS = ex_inquire_int(exoid, EX_INQ_ELEM_BLK);
794: PetscMalloc2(numCS, &csID, numCS, &csSize);
795: PetscStackCallStandard(ex_get_ids,(exoid, EX_ELEM_BLOCK, csID));
796: for (set = 0; set < numCS; ++set) {
797: ex_block block;
799: block.id = csID[set];
800: block.type = EX_ELEM_BLOCK;
801: PetscStackCallStandard(ex_get_block_param,(exoid, &block));
802: csSize[set] = block.num_entry;
803: }
804: VecGetOwnershipRange(vNatural, &xs, &xe);
805: VecGetBlockSize(vNatural, &bs);
806: if (bs > 1) {ISCreateStride(comm, (xe-xs)/bs, xs, bs, &compIS);}
807: for (set = 0; set < numCS; ++set) {
808: PetscInt csLocalSize, c;
810: /* range of indices for set setID[set]: csxs:csxs + csSize[set]-1
811: local slice of zonal values: xs/bs,xm/bs-1
812: intersection: max(xs/bs,csxs),min(xm/bs-1,csxs + csSize[set]-1) */
813: csLocalSize = PetscMax(0, PetscMin(xe/bs, csxs+csSize[set]) - PetscMax(xs/bs, csxs));
814: if (bs == 1) {
815: VecGetArray(vNatural, &varray);
816: PetscStackCallStandard(ex_get_partial_var,(exoid, step, EX_ELEM_BLOCK, offset, csID[set], PetscMax(xs-csxs, 0)+1, csLocalSize, &varray[PetscMax(0, csxs-xs)]));
817: VecRestoreArray(vNatural, &varray);
818: } else {
819: for (c = 0; c < bs; ++c) {
820: ISStrideSetStride(compIS, (xe-xs)/bs, xs+c, bs);
821: VecGetSubVector(vNatural, compIS, &vComp);
822: VecGetArray(vComp, &varray);
823: 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)]));
824: VecRestoreArray(vComp, &varray);
825: VecRestoreSubVector(vNatural, compIS, &vComp);
826: }
827: }
828: csxs += csSize[set];
829: }
830: PetscFree2(csID, csSize);
831: if (bs > 1) {ISDestroy(&compIS);}
832: if (useNatural) {
833: DMPlexNaturalToGlobalBegin(dm, vNatural, v);
834: DMPlexNaturalToGlobalEnd(dm, vNatural, v);
835: DMRestoreGlobalVector(dm, &vNatural);
836: }
837: return(0);
838: }
839: #endif
841: /*@C
842: DMPlexCreateExodusFromFile - Create a DMPlex mesh from an ExodusII file.
844: Collective
846: Input Parameters:
847: + comm - The MPI communicator
848: . filename - The name of the ExodusII file
849: - interpolate - Create faces and edges in the mesh
851: Output Parameter:
852: . dm - The DM object representing the mesh
854: Level: beginner
856: .seealso: DMPLEX, DMCreate(), DMPlexCreateExodus()
857: @*/
858: PetscErrorCode DMPlexCreateExodusFromFile(MPI_Comm comm, const char filename[], PetscBool interpolate, DM *dm)
859: {
860: PetscMPIInt rank;
862: #if defined(PETSC_HAVE_EXODUSII)
863: int CPU_word_size = sizeof(PetscReal), IO_word_size = 0, exoid = -1;
864: float version;
865: #endif
869: MPI_Comm_rank(comm, &rank);
870: #if defined(PETSC_HAVE_EXODUSII)
871: if (!rank) {
872: exoid = ex_open(filename, EX_READ, &CPU_word_size, &IO_word_size, &version);
873: if (exoid <= 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_LIB, "ex_open(\"%s\",...) did not return a valid file ID", filename);
874: }
875: DMPlexCreateExodus(comm, exoid, interpolate, dm);
876: if (!rank) {PetscStackCallStandard(ex_close,(exoid));}
877: #else
878: SETERRQ(comm, PETSC_ERR_SUP, "This method requires ExodusII support. Reconfigure using --download-exodusii");
879: #endif
880: return(0);
881: }
883: /*@
884: DMPlexCreateExodus - Create a DMPlex mesh from an ExodusII file ID.
886: Collective
888: Input Parameters:
889: + comm - The MPI communicator
890: . exoid - The ExodusII id associated with a exodus file and obtained using ex_open
891: - interpolate - Create faces and edges in the mesh
893: Output Parameter:
894: . dm - The DM object representing the mesh
896: Level: beginner
898: .seealso: DMPLEX, DMCreate()
899: @*/
900: PetscErrorCode DMPlexCreateExodus(MPI_Comm comm, PetscInt exoid, PetscBool interpolate, DM *dm)
901: {
902: #if defined(PETSC_HAVE_EXODUSII)
903: PetscMPIInt num_proc, rank;
904: PetscSection coordSection;
905: Vec coordinates;
906: PetscScalar *coords;
907: PetscInt coordSize, v;
909: /* Read from ex_get_init() */
910: char title[PETSC_MAX_PATH_LEN+1];
911: int dim = 0, numVertices = 0, numCells = 0, numHybridCells = 0;
912: int num_cs = 0, num_vs = 0, num_fs = 0;
913: #endif
916: #if defined(PETSC_HAVE_EXODUSII)
917: MPI_Comm_rank(comm, &rank);
918: MPI_Comm_size(comm, &num_proc);
919: DMCreate(comm, dm);
920: DMSetType(*dm, DMPLEX);
921: /* Open EXODUS II file and read basic informations on rank 0, then broadcast to all processors */
922: if (!rank) {
923: PetscMemzero(title,PETSC_MAX_PATH_LEN+1);
924: PetscStackCallStandard(ex_get_init,(exoid, title, &dim, &numVertices, &numCells, &num_cs, &num_vs, &num_fs));
925: if (!num_cs) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Exodus file does not contain any cell set\n");
926: }
927: MPI_Bcast(title, PETSC_MAX_PATH_LEN+1, MPI_CHAR, 0, comm);
928: MPI_Bcast(&dim, 1, MPI_INT, 0, comm);
929: PetscObjectSetName((PetscObject) *dm, title);
930: DMSetDimension(*dm, dim);
931: DMPlexSetChart(*dm, 0, numCells+numVertices);
933: /* Read cell sets information */
934: if (!rank) {
935: PetscInt *cone;
936: int c, cs, ncs, c_loc, v, v_loc;
937: /* Read from ex_get_elem_blk_ids() */
938: int *cs_id, *cs_order;
939: /* Read from ex_get_elem_block() */
940: char buffer[PETSC_MAX_PATH_LEN+1];
941: int num_cell_in_set, num_vertex_per_cell, num_hybrid, num_attr;
942: /* Read from ex_get_elem_conn() */
943: int *cs_connect;
945: /* Get cell sets IDs */
946: PetscMalloc2(num_cs, &cs_id, num_cs, &cs_order);
947: PetscStackCallStandard(ex_get_ids,(exoid, EX_ELEM_BLOCK, cs_id));
948: /* Read the cell set connectivity table and build mesh topology
949: EXO standard requires that cells in cell sets be numbered sequentially and be pairwise disjoint. */
950: /* Check for a hybrid mesh */
951: for (cs = 0, num_hybrid = 0; cs < num_cs; ++cs) {
952: PetscStackCallStandard(ex_get_block,(exoid, EX_ELEM_BLOCK, cs_id[cs], buffer, &num_cell_in_set,&num_vertex_per_cell, 0, 0, &num_attr));
953: switch (dim) {
954: case 3:
955: switch (num_vertex_per_cell) {
956: case 6:
957: cs_order[cs] = cs;
958: numHybridCells += num_cell_in_set;
959: ++num_hybrid;
960: break;
961: default:
962: for (c = cs; c > cs-num_hybrid; --c) cs_order[c] = cs_order[c-1];
963: cs_order[cs-num_hybrid] = cs;
964: }
965: break;
966: default:
967: for (c = cs; c > cs-num_hybrid; --c) cs_order[c] = cs_order[c-1];
968: cs_order[cs-num_hybrid] = cs;
969: }
970: }
971: if (num_hybrid) {DMPlexSetHybridBounds(*dm, numCells-numHybridCells, PETSC_DETERMINE, PETSC_DETERMINE, PETSC_DETERMINE);}
972: /* First set sizes */
973: for (ncs = 0, c = 0; ncs < num_cs; ++ncs) {
974: const PetscInt cs = cs_order[ncs];
975: PetscStackCallStandard(ex_get_block,(exoid, EX_ELEM_BLOCK, cs_id[cs], buffer, &num_cell_in_set,&num_vertex_per_cell, 0, 0, &num_attr));
976: for (c_loc = 0; c_loc < num_cell_in_set; ++c_loc, ++c) {
977: DMPlexSetConeSize(*dm, c, num_vertex_per_cell);
978: }
979: }
980: DMSetUp(*dm);
981: for (ncs = 0, c = 0; ncs < num_cs; ++ncs) {
982: const PetscInt cs = cs_order[ncs];
983: PetscStackCallStandard(ex_get_block,(exoid, EX_ELEM_BLOCK, cs_id[cs], buffer, &num_cell_in_set, &num_vertex_per_cell, 0, 0, &num_attr));
984: PetscMalloc2(num_vertex_per_cell*num_cell_in_set,&cs_connect,num_vertex_per_cell,&cone);
985: PetscStackCallStandard(ex_get_conn,(exoid, EX_ELEM_BLOCK, cs_id[cs], cs_connect,NULL,NULL));
986: /* EXO uses Fortran-based indexing, DMPlex uses C-style and numbers cell first then vertices. */
987: for (c_loc = 0, v = 0; c_loc < num_cell_in_set; ++c_loc, ++c) {
988: for (v_loc = 0; v_loc < num_vertex_per_cell; ++v_loc, ++v) {
989: cone[v_loc] = cs_connect[v]+numCells-1;
990: }
991: if (dim == 3) {
992: /* Tetrahedra are inverted */
993: if (num_vertex_per_cell == 4) {
994: PetscInt tmp = cone[0];
995: cone[0] = cone[1];
996: cone[1] = tmp;
997: }
998: /* Hexahedra are inverted */
999: if (num_vertex_per_cell == 8) {
1000: PetscInt tmp = cone[1];
1001: cone[1] = cone[3];
1002: cone[3] = tmp;
1003: }
1004: }
1005: DMPlexSetCone(*dm, c, cone);
1006: DMSetLabelValue(*dm, "Cell Sets", c, cs_id[cs]);
1007: }
1008: PetscFree2(cs_connect,cone);
1009: }
1010: PetscFree2(cs_id, cs_order);
1011: }
1012: DMPlexSymmetrize(*dm);
1013: DMPlexStratify(*dm);
1014: if (interpolate) {
1015: DM idm;
1017: DMPlexInterpolate(*dm, &idm);
1018: DMDestroy(dm);
1019: *dm = idm;
1020: }
1022: /* Create vertex set label */
1023: if (!rank && (num_vs > 0)) {
1024: int vs, v;
1025: /* Read from ex_get_node_set_ids() */
1026: int *vs_id;
1027: /* Read from ex_get_node_set_param() */
1028: int num_vertex_in_set;
1029: /* Read from ex_get_node_set() */
1030: int *vs_vertex_list;
1032: /* Get vertex set ids */
1033: PetscMalloc1(num_vs, &vs_id);
1034: PetscStackCallStandard(ex_get_ids,(exoid, EX_NODE_SET, vs_id));
1035: for (vs = 0; vs < num_vs; ++vs) {
1036: PetscStackCallStandard(ex_get_set_param,(exoid, EX_NODE_SET, vs_id[vs], &num_vertex_in_set, NULL));
1037: PetscMalloc1(num_vertex_in_set, &vs_vertex_list);
1038: PetscStackCallStandard(ex_get_set,(exoid, EX_NODE_SET, vs_id[vs], vs_vertex_list, NULL));
1039: for (v = 0; v < num_vertex_in_set; ++v) {
1040: DMSetLabelValue(*dm, "Vertex Sets", vs_vertex_list[v]+numCells-1, vs_id[vs]);
1041: }
1042: PetscFree(vs_vertex_list);
1043: }
1044: PetscFree(vs_id);
1045: }
1046: /* Read coordinates */
1047: DMGetCoordinateSection(*dm, &coordSection);
1048: PetscSectionSetNumFields(coordSection, 1);
1049: PetscSectionSetFieldComponents(coordSection, 0, dim);
1050: PetscSectionSetChart(coordSection, numCells, numCells + numVertices);
1051: for (v = numCells; v < numCells+numVertices; ++v) {
1052: PetscSectionSetDof(coordSection, v, dim);
1053: PetscSectionSetFieldDof(coordSection, v, 0, dim);
1054: }
1055: PetscSectionSetUp(coordSection);
1056: PetscSectionGetStorageSize(coordSection, &coordSize);
1057: VecCreate(PETSC_COMM_SELF, &coordinates);
1058: PetscObjectSetName((PetscObject) coordinates, "coordinates");
1059: VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);
1060: VecSetBlockSize(coordinates, dim);
1061: VecSetType(coordinates,VECSTANDARD);
1062: VecGetArray(coordinates, &coords);
1063: if (!rank) {
1064: PetscReal *x, *y, *z;
1066: PetscMalloc3(numVertices,&x,numVertices,&y,numVertices,&z);
1067: PetscStackCallStandard(ex_get_coord,(exoid, x, y, z));
1068: if (dim > 0) {
1069: for (v = 0; v < numVertices; ++v) coords[v*dim+0] = x[v];
1070: }
1071: if (dim > 1) {
1072: for (v = 0; v < numVertices; ++v) coords[v*dim+1] = y[v];
1073: }
1074: if (dim > 2) {
1075: for (v = 0; v < numVertices; ++v) coords[v*dim+2] = z[v];
1076: }
1077: PetscFree3(x,y,z);
1078: }
1079: VecRestoreArray(coordinates, &coords);
1080: DMSetCoordinatesLocal(*dm, coordinates);
1081: VecDestroy(&coordinates);
1083: /* Create side set label */
1084: if (!rank && interpolate && (num_fs > 0)) {
1085: int fs, f, voff;
1086: /* Read from ex_get_side_set_ids() */
1087: int *fs_id;
1088: /* Read from ex_get_side_set_param() */
1089: int num_side_in_set;
1090: /* Read from ex_get_side_set_node_list() */
1091: int *fs_vertex_count_list, *fs_vertex_list;
1092: /* Read side set labels */
1093: char fs_name[MAX_STR_LENGTH+1];
1095: /* Get side set ids */
1096: PetscMalloc1(num_fs, &fs_id);
1097: PetscStackCallStandard(ex_get_ids,(exoid, EX_SIDE_SET, fs_id));
1098: for (fs = 0; fs < num_fs; ++fs) {
1099: PetscStackCallStandard(ex_get_set_param,(exoid, EX_SIDE_SET, fs_id[fs], &num_side_in_set, NULL));
1100: PetscMalloc2(num_side_in_set,&fs_vertex_count_list,num_side_in_set*4,&fs_vertex_list);
1101: PetscStackCallStandard(ex_get_side_set_node_list,(exoid, fs_id[fs], fs_vertex_count_list, fs_vertex_list));
1102: /* Get the specific name associated with this side set ID. */
1103: int fs_name_err = ex_get_name(exoid, EX_SIDE_SET, fs_id[fs], fs_name);
1104: for (f = 0, voff = 0; f < num_side_in_set; ++f) {
1105: const PetscInt *faces = NULL;
1106: PetscInt faceSize = fs_vertex_count_list[f], numFaces;
1107: PetscInt faceVertices[4], v;
1109: if (faceSize > 4) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "ExodusII side cannot have %d > 4 vertices", faceSize);
1110: for (v = 0; v < faceSize; ++v, ++voff) {
1111: faceVertices[v] = fs_vertex_list[voff]+numCells-1;
1112: }
1113: DMPlexGetFullJoin(*dm, faceSize, faceVertices, &numFaces, &faces);
1114: if (numFaces != 1) SETERRQ3(comm, PETSC_ERR_ARG_WRONG, "Invalid ExodusII side %d in set %d maps to %d faces", f, fs, numFaces);
1115: DMSetLabelValue(*dm, "Face Sets", faces[0], fs_id[fs]);
1116: /* Only add the label if one has been detected for this side set. */
1117: if (!fs_name_err) {
1118: DMSetLabelValue(*dm, fs_name, faces[0], fs_id[fs]);
1119: }
1120: DMPlexRestoreJoin(*dm, faceSize, faceVertices, &numFaces, &faces);
1121: }
1122: PetscFree2(fs_vertex_count_list,fs_vertex_list);
1123: }
1124: PetscFree(fs_id);
1125: }
1126: #else
1127: SETERRQ(comm, PETSC_ERR_SUP, "This method requires ExodusII support. Reconfigure using --download-exodusii");
1128: #endif
1129: return(0);
1130: }