Actual source code: plexexodusii.c
petsc-3.6.1 2015-08-06
1: #define PETSCDM_DLL
2: #include <petsc/private/dmpleximpl.h> /*I "petscdmplex.h" I*/
4: #if defined(PETSC_HAVE_EXODUSII)
5: #include <netcdf.h>
6: #include <exodusII.h>
7: #endif
11: /*@C
12: DMPlexCreateExodus - Create a DMPlex mesh from an ExodusII file.
14: Collective on comm
16: Input Parameters:
17: + comm - The MPI communicator
18: . filename - The name of the ExodusII file
19: - interpolate - Create faces and edges in the mesh
21: Output Parameter:
22: . dm - The DM object representing the mesh
24: Level: beginner
26: .keywords: mesh,ExodusII
27: .seealso: DMPLEX, DMCreate(), DMPlexCreateExodus()
28: @*/
29: PetscErrorCode DMPlexCreateExodusFromFile(MPI_Comm comm, const char filename[], PetscBool interpolate, DM *dm)
30: {
31: PetscMPIInt rank;
33: #if defined(PETSC_HAVE_EXODUSII)
34: int CPU_word_size = 0, IO_word_size = 0, exoid = -1;
35: float version;
36: #endif
40: MPI_Comm_rank(comm, &rank);
41: #if defined(PETSC_HAVE_EXODUSII)
42: if (!rank) {
43: exoid = ex_open(filename, EX_READ, &CPU_word_size, &IO_word_size, &version);
44: if (exoid <= 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_LIB, "ex_open(\"%s\",...) did not return a valid file ID", filename);
45: }
46: DMPlexCreateExodus(comm, exoid, interpolate, dm);
47: if (!rank) {ex_close(exoid);}
48: #else
49: SETERRQ(comm, PETSC_ERR_SUP, "This method requires ExodusII support. Reconfigure using --download-exodusii");
50: #endif
51: return(0);
52: }
56: /*@
57: DMPlexCreateExodus - Create a DMPlex mesh from an ExodusII file ID.
59: Collective on comm
61: Input Parameters:
62: + comm - The MPI communicator
63: . exoid - The ExodusII id associated with a exodus file and obtained using ex_open
64: - interpolate - Create faces and edges in the mesh
66: Output Parameter:
67: . dm - The DM object representing the mesh
69: Level: beginner
71: .keywords: mesh,ExodusII
72: .seealso: DMPLEX, DMCreate()
73: @*/
74: PetscErrorCode DMPlexCreateExodus(MPI_Comm comm, PetscInt exoid, PetscBool interpolate, DM *dm)
75: {
76: #if defined(PETSC_HAVE_EXODUSII)
77: PetscMPIInt num_proc, rank;
78: PetscSection coordSection;
79: Vec coordinates;
80: PetscScalar *coords;
81: PetscInt coordSize, v;
83: /* Read from ex_get_init() */
84: char title[PETSC_MAX_PATH_LEN+1];
85: int dim = 0, numVertices = 0, numCells = 0;
86: int num_cs = 0, num_vs = 0, num_fs = 0;
87: #endif
90: #if defined(PETSC_HAVE_EXODUSII)
91: MPI_Comm_rank(comm, &rank);
92: MPI_Comm_size(comm, &num_proc);
93: DMCreate(comm, dm);
94: DMSetType(*dm, DMPLEX);
95: /* Open EXODUS II file and read basic informations on rank 0, then broadcast to all processors */
96: if (!rank) {
97: PetscMemzero(title,(PETSC_MAX_PATH_LEN+1)*sizeof(char));
98: ex_get_init(exoid, title, &dim, &numVertices, &numCells, &num_cs, &num_vs, &num_fs);
99: if (ierr) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"ExodusII ex_get_init() failed with error code %D\n",ierr);
100: if (!num_cs) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Exodus file does not contain any cell set\n");
101: }
102: MPI_Bcast(title, PETSC_MAX_PATH_LEN+1, MPI_CHAR, 0, comm);
103: MPI_Bcast(&dim, 1, MPI_INT, 0, comm);
104: PetscObjectSetName((PetscObject) *dm, title);
105: DMSetDimension(*dm, dim);
106: DMPlexSetChart(*dm, 0, numCells+numVertices);
108: /* Read cell sets information */
109: if (!rank) {
110: PetscInt *cone;
111: int c, cs, c_loc, v, v_loc;
112: /* Read from ex_get_elem_blk_ids() */
113: int *cs_id;
114: /* Read from ex_get_elem_block() */
115: char buffer[PETSC_MAX_PATH_LEN+1];
116: int num_cell_in_set, num_vertex_per_cell, num_attr;
117: /* Read from ex_get_elem_conn() */
118: int *cs_connect;
120: /* Get cell sets IDs */
121: PetscMalloc1(num_cs, &cs_id);
122: ex_get_elem_blk_ids(exoid, cs_id);
123: /* Read the cell set connectivity table and build mesh topology
124: EXO standard requires that cells in cell sets be numbered sequentially and be pairwise disjoint. */
125: /* First set sizes */
126: for (cs = 0, c = 0; cs < num_cs; cs++) {
127: ex_get_elem_block(exoid, cs_id[cs], buffer, &num_cell_in_set, &num_vertex_per_cell, &num_attr);
128: for (c_loc = 0; c_loc < num_cell_in_set; ++c_loc, ++c) {
129: DMPlexSetConeSize(*dm, c, num_vertex_per_cell);
130: }
131: }
132: DMSetUp(*dm);
133: for (cs = 0, c = 0; cs < num_cs; cs++) {
134: ex_get_elem_block(exoid, cs_id[cs], buffer, &num_cell_in_set, &num_vertex_per_cell, &num_attr);
135: PetscMalloc2(num_vertex_per_cell*num_cell_in_set,&cs_connect,num_vertex_per_cell,&cone);
136: ex_get_elem_conn(exoid, cs_id[cs], cs_connect);
137: /* EXO uses Fortran-based indexing, sieve uses C-style and numbers cell first then vertices. */
138: for (c_loc = 0, v = 0; c_loc < num_cell_in_set; ++c_loc, ++c) {
139: for (v_loc = 0; v_loc < num_vertex_per_cell; ++v_loc, ++v) {
140: cone[v_loc] = cs_connect[v]+numCells-1;
141: }
142: if (dim == 3) {
143: /* Tetrahedra are inverted */
144: if (num_vertex_per_cell == 4) {
145: PetscInt tmp = cone[0];
146: cone[0] = cone[1];
147: cone[1] = tmp;
148: }
149: /* Hexahedra are inverted */
150: if (num_vertex_per_cell == 8) {
151: PetscInt tmp = cone[1];
152: cone[1] = cone[3];
153: cone[3] = tmp;
154: }
155: }
156: DMPlexSetCone(*dm, c, cone);
157: DMPlexSetLabelValue(*dm, "Cell Sets", c, cs_id[cs]);
158: }
159: PetscFree2(cs_connect,cone);
160: }
161: PetscFree(cs_id);
162: }
163: DMPlexSymmetrize(*dm);
164: DMPlexStratify(*dm);
165: if (interpolate) {
166: DM idm = NULL;
168: DMPlexInterpolate(*dm, &idm);
169: /* Maintain Cell Sets label */
170: {
171: DMLabel label;
173: DMPlexRemoveLabel(*dm, "Cell Sets", &label);
174: if (label) {DMPlexAddLabel(idm, label);}
175: }
176: DMDestroy(dm);
177: *dm = idm;
178: }
180: /* Create vertex set label */
181: if (!rank && (num_vs > 0)) {
182: int vs, v;
183: /* Read from ex_get_node_set_ids() */
184: int *vs_id;
185: /* Read from ex_get_node_set_param() */
186: int num_vertex_in_set, num_attr;
187: /* Read from ex_get_node_set() */
188: int *vs_vertex_list;
190: /* Get vertex set ids */
191: PetscMalloc1(num_vs, &vs_id);
192: ex_get_node_set_ids(exoid, vs_id);
193: for (vs = 0; vs < num_vs; ++vs) {
194: ex_get_node_set_param(exoid, vs_id[vs], &num_vertex_in_set, &num_attr);
195: PetscMalloc1(num_vertex_in_set, &vs_vertex_list);
196: ex_get_node_set(exoid, vs_id[vs], vs_vertex_list);
197: for (v = 0; v < num_vertex_in_set; ++v) {
198: DMPlexSetLabelValue(*dm, "Vertex Sets", vs_vertex_list[v]+numCells-1, vs_id[vs]);
199: }
200: PetscFree(vs_vertex_list);
201: }
202: PetscFree(vs_id);
203: }
204: /* Read coordinates */
205: DMGetCoordinateSection(*dm, &coordSection);
206: PetscSectionSetNumFields(coordSection, 1);
207: PetscSectionSetFieldComponents(coordSection, 0, dim);
208: PetscSectionSetChart(coordSection, numCells, numCells + numVertices);
209: for (v = numCells; v < numCells+numVertices; ++v) {
210: PetscSectionSetDof(coordSection, v, dim);
211: PetscSectionSetFieldDof(coordSection, v, 0, dim);
212: }
213: PetscSectionSetUp(coordSection);
214: PetscSectionGetStorageSize(coordSection, &coordSize);
215: VecCreate(comm, &coordinates);
216: PetscObjectSetName((PetscObject) coordinates, "coordinates");
217: VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);
218: VecSetBlockSize(coordinates, dim);
219: VecSetType(coordinates,VECSTANDARD);
220: VecGetArray(coordinates, &coords);
221: if (!rank) {
222: float *x, *y, *z;
224: PetscMalloc3(numVertices,&x,numVertices,&y,numVertices,&z);
225: ex_get_coord(exoid, x, y, z);
226: if (dim > 0) {
227: for (v = 0; v < numVertices; ++v) coords[v*dim+0] = x[v];
228: }
229: if (dim > 1) {
230: for (v = 0; v < numVertices; ++v) coords[v*dim+1] = y[v];
231: }
232: if (dim > 2) {
233: for (v = 0; v < numVertices; ++v) coords[v*dim+2] = z[v];
234: }
235: PetscFree3(x,y,z);
236: }
237: VecRestoreArray(coordinates, &coords);
238: DMSetCoordinatesLocal(*dm, coordinates);
239: VecDestroy(&coordinates);
241: /* Create side set label */
242: if (!rank && interpolate && (num_fs > 0)) {
243: int fs, f, voff;
244: /* Read from ex_get_side_set_ids() */
245: int *fs_id;
246: /* Read from ex_get_side_set_param() */
247: int num_side_in_set, num_dist_fact_in_set;
248: /* Read from ex_get_side_set_node_list() */
249: int *fs_vertex_count_list, *fs_vertex_list;
251: /* Get side set ids */
252: PetscMalloc1(num_fs, &fs_id);
253: ex_get_side_set_ids(exoid, fs_id);
254: for (fs = 0; fs < num_fs; ++fs) {
255: ex_get_side_set_param(exoid, fs_id[fs], &num_side_in_set, &num_dist_fact_in_set);
256: PetscMalloc2(num_side_in_set,&fs_vertex_count_list,num_side_in_set*4,&fs_vertex_list);
257: ex_get_side_set_node_list(exoid, fs_id[fs], fs_vertex_count_list, fs_vertex_list);
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: DMPlexSetLabelValue(*dm, "Face Sets", faces[0], fs_id[fs]);
270: DMPlexRestoreJoin(*dm, faceSize, faceVertices, &numFaces, &faces);
271: }
272: PetscFree2(fs_vertex_count_list,fs_vertex_list);
273: }
274: PetscFree(fs_id);
275: }
276: #else
277: SETERRQ(comm, PETSC_ERR_SUP, "This method requires ExodusII support. Reconfigure using --download-exodusii");
278: #endif
279: return(0);
280: }