1 #define PETSCDM_DLL 2 #include <petsc-private/pleximpl.h> /*I "petscdmplex.h" I*/ 3 4 #if defined(PETSC_HAVE_EXODUSII) 5 #include <netcdf.h> 6 #include <exodusII.h> 7 #endif 8 9 PETSC_EXTERN PetscErrorCode DMPlexInterpolate_2D(DM, DM*); 10 PETSC_EXTERN PetscErrorCode DMPlexInterpolate_3D(DM, DM*); 11 12 #undef __FUNCT__ 13 #define __FUNCT__ "DMPlexCreateExodus" 14 /*@ 15 DMPlexCreateExodus - Create a DMPlex mesh from an ExodusII file. 16 17 Collective on comm 18 19 Input Parameters: 20 + comm - The MPI communicator 21 . exoid - The ExodusII id associated with a exodus file and obtained using ex_open 22 - interpolate - Create faces and edges in the mesh 23 24 Output Parameter: 25 . dm - The DM object representing the mesh 26 27 Level: beginner 28 29 .keywords: mesh,ExodusII 30 .seealso: MeshCreate(), MeshCreateExodus() 31 @*/ 32 PetscErrorCode DMPlexCreateExodus(MPI_Comm comm, PetscInt exoid, PetscBool interpolate, DM *dm) 33 { 34 #if defined(PETSC_HAVE_EXODUSII) 35 PetscMPIInt num_proc, rank; 36 PetscSection coordSection; 37 Vec coordinates; 38 PetscScalar *coords; 39 PetscInt coordSize, v; 40 PetscErrorCode ierr; 41 /* Read from ex_get_init() */ 42 char title[PETSC_MAX_PATH_LEN+1]; 43 int dim = 0, numVertices = 0, numCells = 0; 44 int num_cs = 0, num_vs = 0, num_fs = 0; 45 #endif 46 47 PetscFunctionBegin; 48 #if defined(PETSC_HAVE_EXODUSII) 49 ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr); 50 ierr = MPI_Comm_size(comm, &num_proc);CHKERRQ(ierr); 51 ierr = DMCreate(comm, dm);CHKERRQ(ierr); 52 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr); 53 /* Open EXODUS II file and read basic informations on rank 0, then broadcast to all processors */ 54 if (!rank) { 55 ierr = ex_get_init(exoid, title, &dim, &numVertices, &numCells, &num_cs, &num_vs, &num_fs); 56 if (ierr) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"ExodusII ex_get_init() failed with error code %D\n",ierr); 57 if (!num_cs) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Exodus file does not contain any cell set\n"); 58 } 59 ierr = MPI_Bcast(title, PETSC_MAX_PATH_LEN+1, MPI_CHAR, 0, comm);CHKERRQ(ierr); 60 ierr = MPI_Bcast(&dim, 1, MPI_INT, 0, comm);CHKERRQ(ierr); 61 ierr = PetscObjectSetName((PetscObject) *dm, title);CHKERRQ(ierr); 62 ierr = DMPlexSetDimension(*dm, dim);CHKERRQ(ierr); 63 ierr = DMPlexSetChart(*dm, 0, numCells+numVertices);CHKERRQ(ierr); 64 65 /* Read cell sets information */ 66 if (!rank) { 67 PetscInt *cone; 68 int c, cs, c_loc, v, v_loc; 69 /* Read from ex_get_elem_blk_ids() */ 70 int *cs_id; 71 /* Read from ex_get_elem_block() */ 72 char buffer[PETSC_MAX_PATH_LEN+1]; 73 int num_cell_in_set, num_vertex_per_cell, num_attr; 74 /* Read from ex_get_elem_conn() */ 75 int *cs_connect; 76 77 /* Get cell sets IDs */ 78 ierr = PetscMalloc(num_cs * sizeof(int), &cs_id);CHKERRQ(ierr); 79 ierr = ex_get_elem_blk_ids(exoid, cs_id);CHKERRQ(ierr); 80 /* Read the cell set connectivity table and build mesh topology 81 EXO standard requires that cells in cell sets be numbered sequentially and be pairwise disjoint. */ 82 /* First set sizes */ 83 for (cs = 0, c = 0; cs < num_cs; cs++) { 84 ierr = ex_get_elem_block(exoid, cs_id[cs], buffer, &num_cell_in_set, &num_vertex_per_cell, &num_attr);CHKERRQ(ierr); 85 for (c_loc = 0; c_loc < num_cell_in_set; ++c_loc, ++c) { 86 ierr = DMPlexSetConeSize(*dm, c, num_vertex_per_cell);CHKERRQ(ierr); 87 } 88 } 89 ierr = DMSetUp(*dm);CHKERRQ(ierr); 90 for (cs = 0, c = 0; cs < num_cs; cs++) { 91 ierr = ex_get_elem_block(exoid, cs_id[cs], buffer, &num_cell_in_set, &num_vertex_per_cell, &num_attr);CHKERRQ(ierr); 92 ierr = PetscMalloc2(num_vertex_per_cell*num_cell_in_set,int,&cs_connect,num_vertex_per_cell,PetscInt,&cone);CHKERRQ(ierr); 93 ierr = ex_get_elem_conn(exoid, cs_id[cs], cs_connect);CHKERRQ(ierr); 94 /* EXO uses Fortran-based indexing, sieve uses C-style and numbers cell first then vertices. */ 95 for (c_loc = 0, v = 0; c_loc < num_cell_in_set; ++c_loc, ++c) { 96 for (v_loc = 0; v_loc < num_vertex_per_cell; ++v_loc, ++v) { 97 cone[v_loc] = cs_connect[v]+numCells-1; 98 } 99 ierr = DMPlexSetCone(*dm, c, cone);CHKERRQ(ierr); 100 ierr = DMPlexSetLabelValue(*dm, "Cell Sets", c, cs_id[cs]);CHKERRQ(ierr); 101 } 102 ierr = PetscFree2(cs_connect,cone);CHKERRQ(ierr); 103 } 104 ierr = PetscFree(cs_id);CHKERRQ(ierr); 105 } 106 ierr = DMPlexSymmetrize(*dm);CHKERRQ(ierr); 107 ierr = DMPlexStratify(*dm);CHKERRQ(ierr); 108 if (interpolate) { 109 DM idm; 110 111 switch (dim) { 112 case 2: 113 ierr = DMPlexInterpolate_2D(*dm, &idm);CHKERRQ(ierr);break; 114 case 3: 115 ierr = DMPlexInterpolate_3D(*dm, &idm);CHKERRQ(ierr);break; 116 default: 117 SETERRQ1(comm, PETSC_ERR_ARG_OUTOFRANGE, "No mesh interpolation support for dimension %D", dim); 118 } 119 /* Maintain Cell Sets label */ 120 { 121 DMLabel label; 122 123 ierr = DMPlexRemoveLabel(*dm, "Cell Sets", &label);CHKERRQ(ierr); 124 if (label) {ierr = DMPlexAddLabel(idm, label);CHKERRQ(ierr);} 125 } 126 ierr = DMDestroy(dm);CHKERRQ(ierr); 127 *dm = idm; 128 } 129 130 /* Create vertex set label */ 131 if (!rank && (num_vs > 0)) { 132 int vs, v; 133 /* Read from ex_get_node_set_ids() */ 134 int *vs_id; 135 /* Read from ex_get_node_set_param() */ 136 int num_vertex_in_set, num_attr; 137 /* Read from ex_get_node_set() */ 138 int *vs_vertex_list; 139 140 /* Get vertex set ids */ 141 ierr = PetscMalloc(num_vs * sizeof(int), &vs_id);CHKERRQ(ierr); 142 ierr = ex_get_node_set_ids(exoid, vs_id);CHKERRQ(ierr); 143 for (vs = 0; vs < num_vs; ++vs) { 144 ierr = ex_get_node_set_param(exoid, vs_id[vs], &num_vertex_in_set, &num_attr);CHKERRQ(ierr); 145 ierr = PetscMalloc(num_vertex_in_set * sizeof(int), &vs_vertex_list);CHKERRQ(ierr); 146 ierr = ex_get_node_set(exoid, vs_id[vs], vs_vertex_list);CHKERRQ(ierr); 147 for (v = 0; v < num_vertex_in_set; ++v) { 148 ierr = DMPlexSetLabelValue(*dm, "Vertex Sets", vs_vertex_list[v]+numCells-1, vs_id[vs]);CHKERRQ(ierr); 149 } 150 ierr = PetscFree(vs_vertex_list);CHKERRQ(ierr); 151 } 152 ierr = PetscFree(vs_id);CHKERRQ(ierr); 153 } 154 /* Read coordinates */ 155 ierr = DMPlexGetCoordinateSection(*dm, &coordSection);CHKERRQ(ierr); 156 ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr); 157 ierr = PetscSectionSetFieldComponents(coordSection, 0, dim);CHKERRQ(ierr); 158 ierr = PetscSectionSetChart(coordSection, numCells, numCells + numVertices);CHKERRQ(ierr); 159 for (v = numCells; v < numCells+numVertices; ++v) { 160 ierr = PetscSectionSetDof(coordSection, v, dim);CHKERRQ(ierr); 161 ierr = PetscSectionSetFieldDof(coordSection, v, 0, dim);CHKERRQ(ierr); 162 } 163 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 164 ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr); 165 ierr = VecCreate(comm, &coordinates);CHKERRQ(ierr); 166 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr); 167 ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 168 ierr = VecSetFromOptions(coordinates);CHKERRQ(ierr); 169 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 170 if (!rank) { 171 float *x, *y, *z; 172 173 ierr = PetscMalloc3(numVertices,float,&x,numVertices,float,&y,numVertices,float,&z);CHKERRQ(ierr); 174 ierr = ex_get_coord(exoid, x, y, z);CHKERRQ(ierr); 175 if (dim > 0) { 176 for (v = 0; v < numVertices; ++v) coords[v*dim+0] = x[v]; 177 } 178 if (dim > 1) { 179 for (v = 0; v < numVertices; ++v) coords[v*dim+1] = y[v]; 180 } 181 if (dim > 2) { 182 for (v = 0; v < numVertices; ++v) coords[v*dim+2] = z[v]; 183 } 184 ierr = PetscFree3(x,y,z);CHKERRQ(ierr); 185 } 186 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 187 ierr = DMSetCoordinatesLocal(*dm, coordinates);CHKERRQ(ierr); 188 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 189 190 /* Create side set label */ 191 if (!rank && interpolate && (num_fs > 0)) { 192 int fs, f, voff; 193 /* Read from ex_get_side_set_ids() */ 194 int *fs_id; 195 /* Read from ex_get_side_set_param() */ 196 int num_side_in_set, num_dist_fact_in_set; 197 /* Read from ex_get_side_set_node_list() */ 198 int *fs_vertex_count_list, *fs_vertex_list; 199 200 /* Get side set ids */ 201 ierr = PetscMalloc(num_fs * sizeof(int), &fs_id);CHKERRQ(ierr); 202 ierr = ex_get_side_set_ids(exoid, fs_id);CHKERRQ(ierr); 203 for (fs = 0; fs < num_fs; ++fs) { 204 ierr = ex_get_side_set_param(exoid, fs_id[fs], &num_side_in_set, &num_dist_fact_in_set);CHKERRQ(ierr); 205 ierr = PetscMalloc2(num_side_in_set,int,&fs_vertex_count_list,num_side_in_set*4,int,&fs_vertex_list);CHKERRQ(ierr); 206 ierr = ex_get_side_set_node_list(exoid, fs_id[fs], fs_vertex_count_list, fs_vertex_list);CHKERRQ(ierr); 207 for (f = 0, voff = 0; f < num_side_in_set; ++f) { 208 const PetscInt *faces = PETSC_NULL; 209 PetscInt faceSize = fs_vertex_count_list[f], numFaces; 210 PetscInt faceVertices[4], v; 211 212 if (faceSize > 4) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "ExodusII side cannot have %d > 4 vertices", faceSize); 213 for (v = 0; v < faceSize; ++v, ++voff) { 214 faceVertices[v] = fs_vertex_list[voff]+numCells-1; 215 } 216 ierr = DMPlexGetFullJoin(*dm, faceSize, faceVertices, &numFaces, &faces);CHKERRQ(ierr); 217 if (numFaces != 1) SETERRQ3(comm, PETSC_ERR_ARG_WRONG, "Invalid ExodusII side %d in set %d maps to %d faces", f, fs, numFaces); 218 ierr = DMPlexSetLabelValue(*dm, "Face Sets", faces[0], fs_id[fs]);CHKERRQ(ierr); 219 ierr = DMPlexRestoreJoin(*dm, faceSize, faceVertices, &numFaces, &faces);CHKERRQ(ierr); 220 } 221 ierr = PetscFree2(fs_vertex_count_list,fs_vertex_list);CHKERRQ(ierr); 222 } 223 ierr = PetscFree(fs_id);CHKERRQ(ierr); 224 } 225 #else 226 SETERRQ(comm, PETSC_ERR_SUP, "This method requires ExodusII support. Reconfigure using --download-exodusii"); 227 #endif 228 PetscFunctionReturn(0); 229 } 230