Actual source code: plexexodusii.c
petsc-3.8.4 2018-03-24
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: }