Actual source code: plexexodusii.c

petsc-3.8.3 2017-12-09
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: /*@C
 10:   DMPlexCreateExodusFromFile - Create a DMPlex mesh from an ExodusII file.

 12:   Collective on comm

 14:   Input Parameters:
 15: + comm  - The MPI communicator
 16: . filename - The name of the ExodusII file
 17: - interpolate - Create faces and edges in the mesh

 19:   Output Parameter:
 20: . dm  - The DM object representing the mesh

 22:   Level: beginner

 24: .keywords: mesh,ExodusII
 25: .seealso: DMPLEX, DMCreate(), DMPlexCreateExodus()
 26: @*/
 27: PetscErrorCode DMPlexCreateExodusFromFile(MPI_Comm comm, const char filename[], PetscBool interpolate, DM *dm)
 28: {
 29:   PetscMPIInt    rank;
 31: #if defined(PETSC_HAVE_EXODUSII)
 32:   int   CPU_word_size = 0, IO_word_size = 0, exoid = -1;
 33:   float version;
 34: #endif

 38:   MPI_Comm_rank(comm, &rank);
 39: #if defined(PETSC_HAVE_EXODUSII)
 40:   if (!rank) {
 41:     exoid = ex_open(filename, EX_READ, &CPU_word_size, &IO_word_size, &version);
 42:     if (exoid <= 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_LIB, "ex_open(\"%s\",...) did not return a valid file ID", filename);
 43:   }
 44:   DMPlexCreateExodus(comm, exoid, interpolate, dm);
 45:   if (!rank) {ex_close(exoid);}
 46: #else
 47:   SETERRQ(comm, PETSC_ERR_SUP, "This method requires ExodusII support. Reconfigure using --download-exodusii");
 48: #endif
 49:   return(0);
 50: }

 52: /*@
 53:   DMPlexCreateExodus - Create a DMPlex mesh from an ExodusII file ID.

 55:   Collective on comm

 57:   Input Parameters:
 58: + comm  - The MPI communicator
 59: . exoid - The ExodusII id associated with a exodus file and obtained using ex_open
 60: - interpolate - Create faces and edges in the mesh

 62:   Output Parameter:
 63: . dm  - The DM object representing the mesh

 65:   Level: beginner

 67: .keywords: mesh,ExodusII
 68: .seealso: DMPLEX, DMCreate()
 69: @*/
 70: PetscErrorCode DMPlexCreateExodus(MPI_Comm comm, PetscInt exoid, PetscBool interpolate, DM *dm)
 71: {
 72: #if defined(PETSC_HAVE_EXODUSII)
 73:   PetscMPIInt    num_proc, rank;
 74:   PetscSection   coordSection;
 75:   Vec            coordinates;
 76:   PetscScalar    *coords;
 77:   PetscInt       coordSize, v;
 79:   /* Read from ex_get_init() */
 80:   char title[PETSC_MAX_PATH_LEN+1];
 81:   int  dim    = 0, numVertices = 0, numCells = 0;
 82:   int  num_cs = 0, num_vs = 0, num_fs = 0;
 83: #endif

 86: #if defined(PETSC_HAVE_EXODUSII)
 87:   MPI_Comm_rank(comm, &rank);
 88:   MPI_Comm_size(comm, &num_proc);
 89:   DMCreate(comm, dm);
 90:   DMSetType(*dm, DMPLEX);
 91:   /* Open EXODUS II file and read basic informations on rank 0, then broadcast to all processors */
 92:   if (!rank) {
 93:     PetscMemzero(title,(PETSC_MAX_PATH_LEN+1)*sizeof(char));
 94:     ex_get_init(exoid, title, &dim, &numVertices, &numCells, &num_cs, &num_vs, &num_fs);
 95:     if (ierr) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"ExodusII ex_get_init() failed with error code %D\n",ierr);
 96:     if (!num_cs) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Exodus file does not contain any cell set\n");
 97:   }
 98:   MPI_Bcast(title, PETSC_MAX_PATH_LEN+1, MPI_CHAR, 0, comm);
 99:   MPI_Bcast(&dim, 1, MPI_INT, 0, comm);
100:   PetscObjectSetName((PetscObject) *dm, title);
101:   DMSetDimension(*dm, dim);
102:   DMPlexSetChart(*dm, 0, numCells+numVertices);

104:   /* Read cell sets information */
105:   if (!rank) {
106:     PetscInt *cone;
107:     int      c, cs, c_loc, v, v_loc;
108:     /* Read from ex_get_elem_blk_ids() */
109:     int *cs_id;
110:     /* Read from ex_get_elem_block() */
111:     char buffer[PETSC_MAX_PATH_LEN+1];
112:     int  num_cell_in_set, num_vertex_per_cell, num_attr;
113:     /* Read from ex_get_elem_conn() */
114:     int *cs_connect;

116:     /* Get cell sets IDs */
117:     PetscMalloc1(num_cs, &cs_id);
118:     ex_get_ids (exoid, EX_ELEM_BLOCK, cs_id);
119:     /* Read the cell set connectivity table and build mesh topology
120:        EXO standard requires that cells in cell sets be numbered sequentially and be pairwise disjoint. */
121:     /* First set sizes */
122:     for (cs = 0, c = 0; cs < num_cs; cs++) {
123:       ex_get_block(exoid, EX_ELEM_BLOCK, cs_id[cs], buffer, &num_cell_in_set,&num_vertex_per_cell, 0, 0, &num_attr);
124:       for (c_loc = 0; c_loc < num_cell_in_set; ++c_loc, ++c) {
125:         DMPlexSetConeSize(*dm, c, num_vertex_per_cell);
126:       }
127:     }
128:     DMSetUp(*dm);
129:     for (cs = 0, c = 0; cs < num_cs; cs++) {
130:       ex_get_block(exoid, EX_ELEM_BLOCK, cs_id[cs], buffer, &num_cell_in_set, &num_vertex_per_cell, 0, 0, &num_attr);
131:       PetscMalloc2(num_vertex_per_cell*num_cell_in_set,&cs_connect,num_vertex_per_cell,&cone);
132:       ex_get_conn (exoid, EX_ELEM_BLOCK, cs_id[cs], cs_connect,NULL,NULL);
133:       /* EXO uses Fortran-based indexing, DMPlex uses C-style and numbers cell first then vertices. */
134:       for (c_loc = 0, v = 0; c_loc < num_cell_in_set; ++c_loc, ++c) {
135:         for (v_loc = 0; v_loc < num_vertex_per_cell; ++v_loc, ++v) {
136:           cone[v_loc] = cs_connect[v]+numCells-1;
137:         }
138:         if (dim == 3) {
139:           /* Tetrahedra are inverted */
140:           if (num_vertex_per_cell == 4) {
141:             PetscInt tmp = cone[0];
142:             cone[0] = cone[1];
143:             cone[1] = tmp;
144:           }
145:           /* Hexahedra are inverted */
146:           if (num_vertex_per_cell == 8) {
147:             PetscInt tmp = cone[1];
148:             cone[1] = cone[3];
149:             cone[3] = tmp;
150:           }
151:         }
152:         DMPlexSetCone(*dm, c, cone);
153:         DMSetLabelValue(*dm, "Cell Sets", c, cs_id[cs]);
154:       }
155:       PetscFree2(cs_connect,cone);
156:     }
157:     PetscFree(cs_id);
158:   }
159:   DMPlexSymmetrize(*dm);
160:   DMPlexStratify(*dm);
161:   if (interpolate) {
162:     DM idm = NULL;

164:     DMPlexInterpolate(*dm, &idm);
165:     /* Maintain Cell Sets label */
166:     {
167:       DMLabel label;

169:       DMRemoveLabel(*dm, "Cell Sets", &label);
170:       if (label) {DMAddLabel(idm, label);}
171:     }
172:     DMDestroy(dm);
173:     *dm  = idm;
174:   }

176:   /* Create vertex set label */
177:   if (!rank && (num_vs > 0)) {
178:     int vs, v;
179:     /* Read from ex_get_node_set_ids() */
180:     int *vs_id;
181:     /* Read from ex_get_node_set_param() */
182:     int num_vertex_in_set;
183:     /* Read from ex_get_node_set() */
184:     int *vs_vertex_list;

186:     /* Get vertex set ids */
187:     PetscMalloc1(num_vs, &vs_id);
188:     ex_get_ids(exoid, EX_NODE_SET, vs_id);
189:     for (vs = 0; vs < num_vs; ++vs) {
190:       ex_get_set_param(exoid, EX_NODE_SET, vs_id[vs], &num_vertex_in_set, NULL);
191:       PetscMalloc1(num_vertex_in_set, &vs_vertex_list);
192:       ex_get_set(exoid, EX_NODE_SET, vs_id[vs], vs_vertex_list, NULL);
193:       for (v = 0; v < num_vertex_in_set; ++v) {
194:         DMSetLabelValue(*dm, "Vertex Sets", vs_vertex_list[v]+numCells-1, vs_id[vs]);
195:       }
196:       PetscFree(vs_vertex_list);
197:     }
198:     PetscFree(vs_id);
199:   }
200:   /* Read coordinates */
201:   DMGetCoordinateSection(*dm, &coordSection);
202:   PetscSectionSetNumFields(coordSection, 1);
203:   PetscSectionSetFieldComponents(coordSection, 0, dim);
204:   PetscSectionSetChart(coordSection, numCells, numCells + numVertices);
205:   for (v = numCells; v < numCells+numVertices; ++v) {
206:     PetscSectionSetDof(coordSection, v, dim);
207:     PetscSectionSetFieldDof(coordSection, v, 0, dim);
208:   }
209:   PetscSectionSetUp(coordSection);
210:   PetscSectionGetStorageSize(coordSection, &coordSize);
211:   VecCreate(PETSC_COMM_SELF, &coordinates);
212:   PetscObjectSetName((PetscObject) coordinates, "coordinates");
213:   VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);
214:   VecSetBlockSize(coordinates, dim);
215:   VecSetType(coordinates,VECSTANDARD);
216:   VecGetArray(coordinates, &coords);
217:   if (!rank) {
218:     float *x, *y, *z;

220:     PetscMalloc3(numVertices,&x,numVertices,&y,numVertices,&z);
221:     ex_get_coord(exoid, x, y, z);
222:     if (dim > 0) {
223:       for (v = 0; v < numVertices; ++v) coords[v*dim+0] = x[v];
224:     }
225:     if (dim > 1) {
226:       for (v = 0; v < numVertices; ++v) coords[v*dim+1] = y[v];
227:     }
228:     if (dim > 2) {
229:       for (v = 0; v < numVertices; ++v) coords[v*dim+2] = z[v];
230:     }
231:     PetscFree3(x,y,z);
232:   }
233:   VecRestoreArray(coordinates, &coords);
234:   DMSetCoordinatesLocal(*dm, coordinates);
235:   VecDestroy(&coordinates);

237:   /* Create side set label */
238:   if (!rank && interpolate && (num_fs > 0)) {
239:     int fs, f, voff;
240:     /* Read from ex_get_side_set_ids() */
241:     int *fs_id;
242:     /* Read from ex_get_side_set_param() */
243:     int num_side_in_set;
244:     /* Read from ex_get_side_set_node_list() */
245:     int *fs_vertex_count_list, *fs_vertex_list;
246:     /* Read side set labels */
247:     char fs_name[MAX_STR_LENGTH+1];

249:     /* Get side set ids */
250:     PetscMalloc1(num_fs, &fs_id);
251:     ex_get_ids(exoid, EX_SIDE_SET, fs_id);
252:     for (fs = 0; fs < num_fs; ++fs) {
253:       ex_get_set_param(exoid, EX_SIDE_SET, fs_id[fs], &num_side_in_set, NULL);
254:       PetscMalloc2(num_side_in_set,&fs_vertex_count_list,num_side_in_set*4,&fs_vertex_list);
255:       ex_get_side_set_node_list(exoid, fs_id[fs], fs_vertex_count_list, fs_vertex_list);
256:       /* Get the specific name associated with this side set ID. */
257:       int fs_name_err = ex_get_name(exoid, EX_SIDE_SET, fs_id[fs], fs_name);
258:       for (f = 0, voff = 0; f < num_side_in_set; ++f) {
259:         const PetscInt *faces   = NULL;
260:         PetscInt       faceSize = fs_vertex_count_list[f], numFaces;
261:         PetscInt       faceVertices[4], v;

263:         if (faceSize > 4) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "ExodusII side cannot have %d > 4 vertices", faceSize);
264:         for (v = 0; v < faceSize; ++v, ++voff) {
265:           faceVertices[v] = fs_vertex_list[voff]+numCells-1;
266:         }
267:         DMPlexGetFullJoin(*dm, faceSize, faceVertices, &numFaces, &faces);
268:         if (numFaces != 1) SETERRQ3(comm, PETSC_ERR_ARG_WRONG, "Invalid ExodusII side %d in set %d maps to %d faces", f, fs, numFaces);
269:         DMSetLabelValue(*dm, "Face Sets", faces[0], fs_id[fs]);
270:         /* Only add the label if one has been detected for this side set. */
271:         if (!fs_name_err) {
272:           DMSetLabelValue(*dm, fs_name, faces[0], fs_id[fs]);
273:         }
274:         DMPlexRestoreJoin(*dm, faceSize, faceVertices, &numFaces, &faces);
275:       }
276:       PetscFree2(fs_vertex_count_list,fs_vertex_list);
277:     }
278:     PetscFree(fs_id);
279:   }
280: #else
281:   SETERRQ(comm, PETSC_ERR_SUP, "This method requires ExodusII support. Reconfigure using --download-exodusii");
282: #endif
283:   return(0);
284: }