Actual source code: plexexodusii.c

petsc-3.9.0 2018-04-07
Report Typos and Errors
  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(),DMView_Exodus(),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, &section);
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(&section);
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: }