Actual source code: complexexodusii.c
petsc-3.3-p7 2013-05-11
1: #define PETSCDM_DLL
2: #include <petsc-private/compleximpl.h> /*I "petscdmcomplex.h" I*/
4: #ifdef PETSC_HAVE_EXODUSII
5: #include<netcdf.h>
6: #include<exodusII.h>
7: #endif
11: /*@
12: DMComplexCreateExodus - Create a DMComplex mesh from an ExodusII file.
14: Collective on comm
16: Input Parameters:
17: + comm - The MPI communicator
18: - exoid - The ExodusII id associated with a exodus file and obtained using ex_open
20: Output Parameter:
21: . dm - The DM object representing the mesh
23: ExodusII side sets are ignored
25: Interpolated meshes are not supported.
27: Level: beginner
29: .keywords: mesh,ExodusII
30: .seealso: MeshCreate(), MeshCreateExodus()
31: @*/
32: PetscErrorCode DMComplexCreateExodus(MPI_Comm comm, PetscInt exoid, DM *dm)
33: {
34: #if defined(PETSC_HAVE_EXODUSII)
35: DM_Complex *mesh;
36: PetscMPIInt num_proc, rank;
37: PetscScalar *coords;
38: PetscInt coordSize, v;
40: /* Read from ex_get_init() */
41: char title[PETSC_MAX_PATH_LEN+1];
42: int dim = 0, numVertices = 0, numCells = 0;
43: int num_cs = 0, num_vs = 0, num_fs = 0;
44: #endif
47: #if defined(PETSC_HAVE_EXODUSII)
48: MPI_Comm_rank(comm, &rank);
49: MPI_Comm_size(comm, &num_proc);
50: DMCreate(comm, dm);
51: DMSetType(*dm, DMCOMPLEX);
52: mesh = (DM_Complex *) (*dm)->data;
53: /* Open EXODUS II file and read basic informations on rank 0, then broadcast to all processors */
54: if (!rank) {
55: ex_get_init(exoid, title, &dim, &numVertices, &numCells, &num_cs, &num_vs, &num_fs);
56: if (!num_cs) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Exodus file does not contain any cell set\n");
57: }
58: MPI_Bcast(title, PETSC_MAX_PATH_LEN+1, MPI_CHAR, 0, comm);
59: MPI_Bcast(&dim, 1, MPI_INT, 0, comm);
60: PetscObjectSetName((PetscObject) *dm, title);
61: DMComplexSetDimension(*dm, dim);
62: DMComplexSetChart(*dm, 0, numCells+numVertices);
64: /* Read cell sets information */
65: if (!rank) {
66: PetscInt *cone;
67: int c, cs, c_loc, v, v_loc;
68: /* Read from ex_get_elem_blk_ids() */
69: int *cs_id;
70: /* Read from ex_get_elem_block() */
71: char buffer[PETSC_MAX_PATH_LEN+1];
72: int num_cell_in_set, num_vertex_per_cell, num_attr;
73: /* Read from ex_get_elem_conn() */
74: int *cs_connect;
76: /* Get cell sets IDs */
77: PetscMalloc(num_cs * sizeof(int), &cs_id);
78: ex_get_elem_blk_ids(exoid, cs_id);
79: /* Read the cell set connectivity table and build mesh topology
80: EXO standard requires that cells in cell sets be numbered sequentially and be pairwise disjoint. */
81: /* First set sizes */
82: for(cs = 0, c = 0; cs < num_cs; cs++) {
83: ex_get_elem_block(exoid, cs_id[cs], buffer, &num_cell_in_set, &num_vertex_per_cell, &num_attr);
84: for(c_loc = 0; c_loc < num_cell_in_set; ++c_loc, ++c) {
85: DMComplexSetConeSize(*dm, c, num_vertex_per_cell);
86: }
87: }
88: DMSetUp(*dm);
89: for(cs = 0, c = 0; cs < num_cs; cs++) {
90: ex_get_elem_block(exoid, cs_id[cs], buffer, &num_cell_in_set, &num_vertex_per_cell, &num_attr);
91: PetscMalloc2(num_vertex_per_cell*num_cell_in_set,int,&cs_connect,num_vertex_per_cell,PetscInt,&cone);
92: ex_get_elem_conn(exoid, cs_id[cs], cs_connect);
93: /* EXO uses Fortran-based indexing, sieve uses C-style and numbers cell first then vertices. */
94: for(c_loc = 0, v = 0; c_loc < num_cell_in_set; ++c_loc, ++c) {
95: for(v_loc = 0; v_loc < num_vertex_per_cell; ++v_loc, ++v) {
96: cone[v_loc] = cs_connect[v]+numCells-1;
97: }
98: DMComplexSetCone(*dm, c, cone);
99: DMComplexSetLabelValue(*dm, "Cell Sets", c, cs_id[cs]);
100: }
101: PetscFree2(cs_connect,cone);
102: }
103: PetscFree(cs_id);
104: }
105: DMComplexSymmetrize(*dm);
106: DMComplexStratify(*dm);
108: /* Create vertex set label */
109: if (!rank && (num_vs > 0)) {
110: int vs, v;
111: /* Read from ex_get_node_set_ids() */
112: int *vs_id;
113: /* Read from ex_get_node_set_param() */
114: int num_vertex_in_set, num_attr;
115: /* Read from ex_get_node_set() */
116: int *vs_vertex_list;
118: /* Get vertex set ids */
119: PetscMalloc(num_vs * sizeof(int), &vs_id);
120: ex_get_node_set_ids(exoid, vs_id);
121: for(vs = 0; vs < num_vs; vs++) {
122: ex_get_node_set_param(exoid, vs_id[vs], &num_vertex_in_set, &num_attr);
123: PetscMalloc(num_vertex_in_set * sizeof(int), &vs_vertex_list);
124: ex_get_node_set(exoid, vs_id[vs], vs_vertex_list);
125: for(v = 0; v < num_vertex_in_set; v++) {
126: DMComplexSetLabelValue(*dm, "Vertex Sets", vs_vertex_list[v]+numCells-1, vs_id[vs]);
127: }
128: PetscFree(vs_vertex_list);
129: }
130: PetscFree(vs_id);
131: }
132: /* Read coordinates */
133: PetscSectionSetChart(mesh->coordSection, numCells, numCells + numVertices);
134: for(v = numCells; v < numCells+numVertices; ++v) {
135: PetscSectionSetDof(mesh->coordSection, v, dim);
136: }
137: PetscSectionSetUp(mesh->coordSection);
138: PetscSectionGetStorageSize(mesh->coordSection, &coordSize);
139: VecSetSizes(mesh->coordinates, coordSize, PETSC_DETERMINE);
140: VecSetFromOptions(mesh->coordinates);
141: VecGetArray(mesh->coordinates, &coords);
142: if (rank == 0) {
143: PetscReal *x, *y, *z;
145: PetscMalloc3(numVertices,PetscReal,&x,numVertices,PetscReal,&y,numVertices,PetscReal,&z);
146: ex_get_coord(exoid, x, y, z);
147: if (dim > 0) {for (v = 0; v < numVertices; ++v) {coords[v*dim+0] = x[v];}}
148: if (dim > 1) {for (v = 0; v < numVertices; ++v) {coords[v*dim+1] = y[v];}}
149: if (dim > 2) {for (v = 0; v < numVertices; ++v) {coords[v*dim+2] = z[v];}}
150: PetscFree3(x,y,z);
151: }
152: VecRestoreArray(mesh->coordinates, &coords);
153: #else
154: SETERRQ(comm, PETSC_ERR_SUP, "This method requires ExodusII support. Reconfigure using --download-exodusii");
155: #endif
156: return(0);
157: }