1 #define PETSCDM_DLL 2 #include <petsc/private/dmpleximpl.h> /*I "petscdmplex.h" I*/ 3 4 #if defined(PETSC_HAVE_EXODUSII) 5 #include <netcdf.h> 6 #include <exodusII.h> 7 #endif 8 9 /*@C 10 DMPlexCreateExodus - Create a DMPlex mesh from an ExodusII file. 11 12 Collective on comm 13 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 18 19 Output Parameter: 20 . dm - The DM object representing the mesh 21 22 Level: beginner 23 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; 30 PetscErrorCode ierr; 31 #if defined(PETSC_HAVE_EXODUSII) 32 int CPU_word_size = 0, IO_word_size = 0, exoid = -1; 33 float version; 34 #endif 35 36 PetscFunctionBegin; 37 PetscValidCharPointer(filename, 2); 38 ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr); 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 ierr = DMPlexCreateExodus(comm, exoid, interpolate, dm);CHKERRQ(ierr); 45 if (!rank) {ierr = ex_close(exoid);CHKERRQ(ierr);} 46 #else 47 SETERRQ(comm, PETSC_ERR_SUP, "This method requires ExodusII support. Reconfigure using --download-exodusii"); 48 #endif 49 PetscFunctionReturn(0); 50 } 51 52 /*@ 53 DMPlexCreateExodus - Create a DMPlex mesh from an ExodusII file ID. 54 55 Collective on comm 56 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 61 62 Output Parameter: 63 . dm - The DM object representing the mesh 64 65 Level: beginner 66 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; 78 PetscErrorCode ierr; 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 84 85 PetscFunctionBegin; 86 #if defined(PETSC_HAVE_EXODUSII) 87 ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr); 88 ierr = MPI_Comm_size(comm, &num_proc);CHKERRQ(ierr); 89 ierr = DMCreate(comm, dm);CHKERRQ(ierr); 90 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr); 91 /* Open EXODUS II file and read basic informations on rank 0, then broadcast to all processors */ 92 if (!rank) { 93 ierr = PetscMemzero(title,(PETSC_MAX_PATH_LEN+1)*sizeof(char));CHKERRQ(ierr); 94 ierr = 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 ierr = MPI_Bcast(title, PETSC_MAX_PATH_LEN+1, MPI_CHAR, 0, comm);CHKERRQ(ierr); 99 ierr = MPI_Bcast(&dim, 1, MPI_INT, 0, comm);CHKERRQ(ierr); 100 ierr = PetscObjectSetName((PetscObject) *dm, title);CHKERRQ(ierr); 101 ierr = DMSetDimension(*dm, dim);CHKERRQ(ierr); 102 ierr = DMPlexSetChart(*dm, 0, numCells+numVertices);CHKERRQ(ierr); 103 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; 115 116 /* Get cell sets IDs */ 117 ierr = PetscMalloc1(num_cs, &cs_id);CHKERRQ(ierr); 118 ierr = ex_get_elem_blk_ids(exoid, cs_id);CHKERRQ(ierr); 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 ierr = ex_get_elem_block(exoid, cs_id[cs], buffer, &num_cell_in_set, &num_vertex_per_cell, &num_attr);CHKERRQ(ierr); 124 for (c_loc = 0; c_loc < num_cell_in_set; ++c_loc, ++c) { 125 ierr = DMPlexSetConeSize(*dm, c, num_vertex_per_cell);CHKERRQ(ierr); 126 } 127 } 128 ierr = DMSetUp(*dm);CHKERRQ(ierr); 129 for (cs = 0, c = 0; cs < num_cs; cs++) { 130 ierr = ex_get_elem_block(exoid, cs_id[cs], buffer, &num_cell_in_set, &num_vertex_per_cell, &num_attr);CHKERRQ(ierr); 131 ierr = PetscMalloc2(num_vertex_per_cell*num_cell_in_set,&cs_connect,num_vertex_per_cell,&cone);CHKERRQ(ierr); 132 ierr = ex_get_elem_conn(exoid, cs_id[cs], cs_connect);CHKERRQ(ierr); 133 /* EXO uses Fortran-based indexing, sieve 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 ierr = DMPlexSetCone(*dm, c, cone);CHKERRQ(ierr); 153 ierr = DMSetLabelValue(*dm, "Cell Sets", c, cs_id[cs]);CHKERRQ(ierr); 154 } 155 ierr = PetscFree2(cs_connect,cone);CHKERRQ(ierr); 156 } 157 ierr = PetscFree(cs_id);CHKERRQ(ierr); 158 } 159 ierr = DMPlexSymmetrize(*dm);CHKERRQ(ierr); 160 ierr = DMPlexStratify(*dm);CHKERRQ(ierr); 161 if (interpolate) { 162 DM idm = NULL; 163 164 ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr); 165 /* Maintain Cell Sets label */ 166 { 167 DMLabel label; 168 169 ierr = DMRemoveLabel(*dm, "Cell Sets", &label);CHKERRQ(ierr); 170 if (label) {ierr = DMAddLabel(idm, label);CHKERRQ(ierr);} 171 } 172 ierr = DMDestroy(dm);CHKERRQ(ierr); 173 *dm = idm; 174 } 175 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, num_attr; 183 /* Read from ex_get_node_set() */ 184 int *vs_vertex_list; 185 186 /* Get vertex set ids */ 187 ierr = PetscMalloc1(num_vs, &vs_id);CHKERRQ(ierr); 188 ierr = ex_get_node_set_ids(exoid, vs_id);CHKERRQ(ierr); 189 for (vs = 0; vs < num_vs; ++vs) { 190 ierr = ex_get_node_set_param(exoid, vs_id[vs], &num_vertex_in_set, &num_attr);CHKERRQ(ierr); 191 ierr = PetscMalloc1(num_vertex_in_set, &vs_vertex_list);CHKERRQ(ierr); 192 ierr = ex_get_node_set(exoid, vs_id[vs], vs_vertex_list);CHKERRQ(ierr); 193 for (v = 0; v < num_vertex_in_set; ++v) { 194 ierr = DMSetLabelValue(*dm, "Vertex Sets", vs_vertex_list[v]+numCells-1, vs_id[vs]);CHKERRQ(ierr); 195 } 196 ierr = PetscFree(vs_vertex_list);CHKERRQ(ierr); 197 } 198 ierr = PetscFree(vs_id);CHKERRQ(ierr); 199 } 200 /* Read coordinates */ 201 ierr = DMGetCoordinateSection(*dm, &coordSection);CHKERRQ(ierr); 202 ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr); 203 ierr = PetscSectionSetFieldComponents(coordSection, 0, dim);CHKERRQ(ierr); 204 ierr = PetscSectionSetChart(coordSection, numCells, numCells + numVertices);CHKERRQ(ierr); 205 for (v = numCells; v < numCells+numVertices; ++v) { 206 ierr = PetscSectionSetDof(coordSection, v, dim);CHKERRQ(ierr); 207 ierr = PetscSectionSetFieldDof(coordSection, v, 0, dim);CHKERRQ(ierr); 208 } 209 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 210 ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr); 211 ierr = VecCreate(PETSC_COMM_SELF, &coordinates);CHKERRQ(ierr); 212 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr); 213 ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 214 ierr = VecSetBlockSize(coordinates, dim);CHKERRQ(ierr); 215 ierr = VecSetType(coordinates,VECSTANDARD);CHKERRQ(ierr); 216 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 217 if (!rank) { 218 float *x, *y, *z; 219 220 ierr = PetscMalloc3(numVertices,&x,numVertices,&y,numVertices,&z);CHKERRQ(ierr); 221 ierr = ex_get_coord(exoid, x, y, z);CHKERRQ(ierr); 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 ierr = PetscFree3(x,y,z);CHKERRQ(ierr); 232 } 233 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 234 ierr = DMSetCoordinatesLocal(*dm, coordinates);CHKERRQ(ierr); 235 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 236 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, num_dist_fact_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]; 248 249 /* Get side set ids */ 250 ierr = PetscMalloc1(num_fs, &fs_id);CHKERRQ(ierr); 251 ierr = ex_get_side_set_ids(exoid, fs_id);CHKERRQ(ierr); 252 for (fs = 0; fs < num_fs; ++fs) { 253 ierr = ex_get_side_set_param(exoid, fs_id[fs], &num_side_in_set, &num_dist_fact_in_set);CHKERRQ(ierr); 254 ierr = PetscMalloc2(num_side_in_set,&fs_vertex_count_list,num_side_in_set*4,&fs_vertex_list);CHKERRQ(ierr); 255 ierr = ex_get_side_set_node_list(exoid, fs_id[fs], fs_vertex_count_list, fs_vertex_list);CHKERRQ(ierr); 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; 262 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 ierr = DMPlexGetFullJoin(*dm, faceSize, faceVertices, &numFaces, &faces);CHKERRQ(ierr); 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 ierr = DMSetLabelValue(*dm, "Face Sets", faces[0], fs_id[fs]);CHKERRQ(ierr); 270 /* Only add the label if one has been detected for this side set. */ 271 if (!fs_name_err) { 272 ierr = DMSetLabelValue(*dm, fs_name, faces[0], fs_id[fs]);CHKERRQ(ierr); 273 } 274 ierr = DMPlexRestoreJoin(*dm, faceSize, faceVertices, &numFaces, &faces);CHKERRQ(ierr); 275 } 276 ierr = PetscFree2(fs_vertex_count_list,fs_vertex_list);CHKERRQ(ierr); 277 } 278 ierr = PetscFree(fs_id);CHKERRQ(ierr); 279 } 280 #else 281 SETERRQ(comm, PETSC_ERR_SUP, "This method requires ExodusII support. Reconfigure using --download-exodusii"); 282 #endif 283 PetscFunctionReturn(0); 284 } 285