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