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 #if defined(PETSC_HAVE_EXODUSII) 10 /* 11 EXOGetVarIndex - Locate a result in an exodus file based on its name 12 13 Not collective 14 15 Input Parameters: 16 + exoid - the exodus id of a file (obtained from ex_open or ex_create for instance) 17 . obj_type - the type of entity for instance EX_NODAL, EX_ELEM_BLOCK 18 - name - the name of the result 19 20 Output Parameters: 21 . varIndex - the location in the exodus file of the result 22 23 Notes: 24 The exodus variable index is obtained by comparing name and the 25 names of zonal variables declared in the exodus file. For instance if name is "V" 26 the location in the exodus file will be the first match of "V", "V_X", "V_XX", "V_1", or "V_11" 27 amongst all variables of type obj_type. 28 29 Level: beginner 30 31 .keywords: mesh,ExodusII 32 .seealso: EXOGetVarIndex(),DMPlexView_ExodusII_Internal(),VecViewPlex_ExodusII_Nodal_Internal(),VecLoadNodal_PlexEXO(),VecLoadZonal_PlexEXO() 33 */ 34 static PetscErrorCode EXOGetVarIndex_Private(int exoid, ex_entity_type obj_type, const char name[], int *varIndex) 35 { 36 int num_vars, i, j; 37 char ext_name[MAX_STR_LENGTH+1], var_name[MAX_STR_LENGTH+1]; 38 const int num_suffix = 5; 39 char *suffix[5]; 40 PetscBool flg; 41 PetscErrorCode ierr; 42 43 PetscFunctionBegin; 44 suffix[0] = (char *) ""; 45 suffix[1] = (char *) "_X"; 46 suffix[2] = (char *) "_XX"; 47 suffix[3] = (char *) "_1"; 48 suffix[4] = (char *) "_11"; 49 50 *varIndex = 0; 51 PetscStackCallStandard(ex_get_variable_param,(exoid, obj_type, &num_vars)); 52 for (i = 0; i < num_vars; ++i) { 53 PetscStackCallStandard(ex_get_variable_name,(exoid, obj_type, i+1, var_name)); 54 for (j = 0; j < num_suffix; ++j){ 55 ierr = PetscStrncpy(ext_name, name, MAX_STR_LENGTH);CHKERRQ(ierr); 56 ierr = PetscStrlcat(ext_name, suffix[j], MAX_STR_LENGTH);CHKERRQ(ierr); 57 ierr = PetscStrcasecmp(ext_name, var_name, &flg); 58 if (flg) { 59 *varIndex = i+1; 60 PetscFunctionReturn(0); 61 } 62 } 63 } 64 SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unable to locate variable %s in ExodusII file.", name); 65 PetscFunctionReturn(-1); 66 } 67 68 /* 69 DMPlexView_ExodusII_Internal - Write a DM to disk in exodus format 70 71 Collective on dm 72 73 Input Parameters: 74 + dm - The dm to be written 75 . exoid - the exodus id of a file (obtained from ex_open or ex_create for instance) 76 - degree - the degree of the interpolation space 77 78 Notes: 79 Not all DM can be written to disk this way. For instance, exodus assume that element blocks (mapped to "Cell sets" labels) 80 consists of sequentially numbered cells. If this is not the case, the exodus file will be corrupted. 81 82 If the dm has been distributed and exoid points to different files on each MPI rank, only the "local" part of the DM 83 (including "ghost" cells and vertices) will be written. If all exoid points to the same file, the resulting file will 84 probably be corrupted. 85 86 DMPlex only represents geometry while most post-processing software expect that a mesh also provides information 87 on the discretization space. This function assumes that the file represents Lagrange finite elements of order 1 or 2. 88 It should be extended to use PetscFE objects. 89 90 This function will only handle TRI, TET, QUAD and HEX cells. 91 Level: beginner 92 93 .keywords: mesh, ExodusII 94 .seealso: EXOGetVarIndex_Private(), VecViewPlex_ExodusII_Nodal_Internal(), VecLoadNodal_PlexEXO(), VecViewZonal_PlexEXO(), VecLoadZonal_PlexEXO() 95 */ 96 PetscErrorCode DMPlexView_ExodusII_Internal(DM dm, int exoid, PetscInt degree) 97 { 98 enum ElemType {TRI, QUAD, TET, HEX}; 99 MPI_Comm comm; 100 /* Connectivity Variables */ 101 PetscInt cellsNotInConnectivity; 102 /* Cell Sets */ 103 DMLabel csLabel; 104 IS csIS; 105 const PetscInt *csIdx; 106 PetscInt num_cs, cs; 107 enum ElemType *type; 108 PetscBool hasLabel; 109 /* Coordinate Variables */ 110 DM cdm; 111 PetscSection section; 112 Vec coord; 113 PetscInt **nodes; 114 PetscInt depth, d, dim, skipCells = 0; 115 PetscInt pStart, pEnd, p, cStart, cEnd, numCells, vStart, vEnd, numVertices, eStart, eEnd, numEdges, fStart, fEnd, numFaces, numNodes; 116 PetscInt num_vs, num_fs; 117 PetscMPIInt rank, size; 118 const char *dmName; 119 PetscInt nodesTriP1[4] = {3,0,0,0}; 120 PetscInt nodesTriP2[4] = {3,3,0,0}; 121 PetscInt nodesQuadP1[4] = {4,0,0,0}; 122 PetscInt nodesQuadP2[4] = {4,4,0,1}; 123 PetscInt nodesTetP1[4] = {4,0,0,0}; 124 PetscInt nodesTetP2[4] = {4,6,0,0}; 125 PetscInt nodesHexP1[4] = {8,0,0,0}; 126 PetscInt nodesHexP2[4] = {8,12,6,1}; 127 PetscErrorCode ierr; 128 129 PetscFunctionBegin; 130 ierr = PetscObjectGetComm((PetscObject)dm, &comm);CHKERRQ(ierr); 131 ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr); 132 ierr = MPI_Comm_size(comm, &size);CHKERRQ(ierr); 133 /* --- Get DM info --- */ 134 ierr = PetscObjectGetName((PetscObject) dm, &dmName);CHKERRQ(ierr); 135 ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr); 136 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 137 ierr = DMPlexGetChart(dm, &pStart, &pEnd);CHKERRQ(ierr); 138 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 139 ierr = DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);CHKERRQ(ierr); 140 ierr = DMPlexGetDepthStratum(dm, 1, &eStart, &eEnd);CHKERRQ(ierr); 141 ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr); 142 numCells = cEnd - cStart; 143 numEdges = eEnd - eStart; 144 numVertices = vEnd - vStart; 145 if (depth == 3) {numFaces = fEnd - fStart;} 146 else {numFaces = 0;} 147 ierr = DMGetLabelSize(dm, "Cell Sets", &num_cs);CHKERRQ(ierr); 148 ierr = DMGetLabelSize(dm, "Vertex Sets", &num_vs);CHKERRQ(ierr); 149 ierr = DMGetLabelSize(dm, "Face Sets", &num_fs);CHKERRQ(ierr); 150 ierr = DMGetCoordinatesLocal(dm, &coord);CHKERRQ(ierr); 151 ierr = DMGetCoordinateDM(dm, &cdm);CHKERRQ(ierr); 152 if (num_cs > 0) { 153 ierr = DMGetLabel(dm, "Cell Sets", &csLabel);CHKERRQ(ierr); 154 ierr = DMLabelGetValueIS(csLabel, &csIS);CHKERRQ(ierr); 155 ierr = ISGetIndices(csIS, &csIdx);CHKERRQ(ierr); 156 } 157 ierr = PetscMalloc1(num_cs, &nodes);CHKERRQ(ierr); 158 /* Set element type for each block and compute total number of nodes */ 159 ierr = PetscMalloc1(num_cs, &type);CHKERRQ(ierr); 160 numNodes = numVertices; 161 if (degree == 2) {numNodes += numEdges;} 162 cellsNotInConnectivity = numCells; 163 for (cs = 0; cs < num_cs; ++cs) { 164 IS stratumIS; 165 const PetscInt *cells; 166 PetscScalar *xyz = NULL; 167 PetscInt csSize, closureSize; 168 169 ierr = DMLabelGetStratumIS(csLabel, csIdx[cs], &stratumIS);CHKERRQ(ierr); 170 ierr = ISGetIndices(stratumIS, &cells);CHKERRQ(ierr); 171 ierr = ISGetSize(stratumIS, &csSize);CHKERRQ(ierr); 172 ierr = DMPlexVecGetClosure(cdm, NULL, coord, cells[0], &closureSize, &xyz);CHKERRQ(ierr); 173 switch (dim) { 174 case 2: 175 if (closureSize == 3*dim) {type[cs] = TRI;} 176 else if (closureSize == 4*dim) {type[cs] = QUAD;} 177 else SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Number of vertices %D in dimension %D has no ExodusII type", closureSize/dim, dim); 178 break; 179 case 3: 180 if (closureSize == 4*dim) {type[cs] = TET;} 181 else if (closureSize == 8*dim) {type[cs] = HEX;} 182 else SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Number of vertices %D in dimension %D has no ExodusII type", closureSize/dim, dim); 183 break; 184 default: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Dimension %D not handled by ExodusII viewer", dim); 185 } 186 if ((degree == 2) && (type[cs] == QUAD)) {numNodes += csSize;} 187 if ((degree == 2) && (type[cs] == HEX)) {numNodes += csSize; numNodes += numFaces;} 188 ierr = DMPlexVecRestoreClosure(cdm, NULL, coord, cells[0], &closureSize, &xyz);CHKERRQ(ierr); 189 /* Set nodes and Element type */ 190 if (type[cs] == TRI) { 191 if (degree == 1) nodes[cs] = nodesTriP1; 192 else if (degree == 2) nodes[cs] = nodesTriP2; 193 } else if (type[cs] == QUAD) { 194 if (degree == 1) nodes[cs] = nodesQuadP1; 195 else if (degree == 2) nodes[cs] = nodesQuadP2; 196 } else if (type[cs] == TET) { 197 if (degree == 1) nodes[cs] = nodesTetP1; 198 else if (degree == 2) nodes[cs] = nodesTetP2; 199 } else if (type[cs] == HEX) { 200 if (degree == 1) nodes[cs] = nodesHexP1; 201 else if (degree == 2) nodes[cs] = nodesHexP2; 202 } 203 /* Compute the number of cells not in the connectivity table */ 204 cellsNotInConnectivity -= nodes[cs][3]*csSize; 205 206 ierr = ISRestoreIndices(stratumIS, &cells);CHKERRQ(ierr); 207 ierr = ISDestroy(&stratumIS);CHKERRQ(ierr); 208 } 209 if (num_cs > 0) {PetscStackCallStandard(ex_put_init,(exoid, dmName, dim, numNodes, numCells, num_cs, num_vs, num_fs));} 210 /* --- Connectivity --- */ 211 for (cs = 0; cs < num_cs; ++cs) { 212 IS stratumIS; 213 const PetscInt *cells; 214 PetscInt *connect, off = 0; 215 PetscInt edgesInClosure = 0, facesInClosure = 0, verticesInClosure = 0; 216 PetscInt csSize, c, connectSize, closureSize; 217 char *elem_type = NULL; 218 char elem_type_tri3[] = "TRI3", elem_type_quad4[] = "QUAD4"; 219 char elem_type_tri6[] = "TRI6", elem_type_quad9[] = "QUAD9"; 220 char elem_type_tet4[] = "TET4", elem_type_hex8[] = "HEX8"; 221 char elem_type_tet10[] = "TET10", elem_type_hex27[] = "HEX27"; 222 223 ierr = DMLabelGetStratumIS(csLabel, csIdx[cs], &stratumIS);CHKERRQ(ierr); 224 ierr = ISGetIndices(stratumIS, &cells);CHKERRQ(ierr); 225 ierr = ISGetSize(stratumIS, &csSize);CHKERRQ(ierr); 226 /* Set Element type */ 227 if (type[cs] == TRI) { 228 if (degree == 1) elem_type = elem_type_tri3; 229 else if (degree == 2) elem_type = elem_type_tri6; 230 } else if (type[cs] == QUAD) { 231 if (degree == 1) elem_type = elem_type_quad4; 232 else if (degree == 2) elem_type = elem_type_quad9; 233 } else if (type[cs] == TET) { 234 if (degree == 1) elem_type = elem_type_tet4; 235 else if (degree == 2) elem_type = elem_type_tet10; 236 } else if (type[cs] == HEX) { 237 if (degree == 1) elem_type = elem_type_hex8; 238 else if (degree == 2) elem_type = elem_type_hex27; 239 } 240 connectSize = nodes[cs][0] + nodes[cs][1] + nodes[cs][2] + nodes[cs][3]; 241 ierr = PetscMalloc1(PetscMax(27,connectSize)*csSize, &connect);CHKERRQ(ierr); 242 PetscStackCallStandard(ex_put_block,(exoid, EX_ELEM_BLOCK, csIdx[cs], elem_type, csSize, connectSize, 0, 0, 1)); 243 /* Find number of vertices, edges, and faces in the closure */ 244 verticesInClosure = nodes[cs][0]; 245 if (depth > 1) { 246 if (dim == 2) { 247 ierr = DMPlexGetConeSize(dm, cells[0], &edgesInClosure);CHKERRQ(ierr); 248 } else if (dim == 3) { 249 PetscInt *closure = NULL; 250 251 ierr = DMPlexGetConeSize(dm, cells[0], &facesInClosure);CHKERRQ(ierr); 252 ierr = DMPlexGetTransitiveClosure(dm, cells[0], PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 253 edgesInClosure = closureSize - facesInClosure - 1 - verticesInClosure; 254 ierr = DMPlexRestoreTransitiveClosure(dm, cells[0], PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 255 } 256 } 257 /* Get connectivity for each cell */ 258 for (c = 0; c < csSize; ++c) { 259 PetscInt *closure = NULL; 260 PetscInt temp, i; 261 262 ierr = DMPlexGetTransitiveClosure(dm, cells[c], PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 263 for (i = 0; i < connectSize; ++i) { 264 if (i < nodes[cs][0]) {/* Vertices */ 265 connect[i+off] = closure[(i+edgesInClosure+facesInClosure+1)*2] + 1; 266 connect[i+off] -= cellsNotInConnectivity; 267 } else if (i < nodes[cs][0]+nodes[cs][1]) { /* Edges */ 268 connect[i+off] = closure[(i-verticesInClosure+facesInClosure+1)*2] + 1; 269 if (nodes[cs][2] == 0) connect[i+off] -= numFaces; 270 connect[i+off] -= cellsNotInConnectivity; 271 } else if (i < nodes[cs][0]+nodes[cs][1]+nodes[cs][3]){ /* Cells */ 272 connect[i+off] = closure[0] + 1; 273 connect[i+off] -= skipCells; 274 } else if (i < nodes[cs][0]+nodes[cs][1]+nodes[cs][3]+nodes[cs][2]){ /* Faces */ 275 connect[i+off] = closure[(i-edgesInClosure-verticesInClosure)*2] + 1; 276 connect[i+off] -= cellsNotInConnectivity; 277 } else { 278 connect[i+off] = -1; 279 } 280 } 281 /* Tetrahedra are inverted */ 282 if (type[cs] == TET) { 283 temp = connect[0+off]; connect[0+off] = connect[1+off]; connect[1+off] = temp; 284 if (degree == 2) { 285 temp = connect[5+off]; connect[5+off] = connect[6+off]; connect[6+off] = temp; 286 temp = connect[7+off]; connect[7+off] = connect[8+off]; connect[8+off] = temp; 287 } 288 } 289 /* Hexahedra are inverted */ 290 if (type[cs] == HEX) { 291 temp = connect[1+off]; connect[1+off] = connect[3+off]; connect[3+off] = temp; 292 if (degree == 2) { 293 temp = connect[8+off]; connect[8+off] = connect[11+off]; connect[11+off] = temp; 294 temp = connect[9+off]; connect[9+off] = connect[10+off]; connect[10+off] = temp; 295 temp = connect[16+off]; connect[16+off] = connect[17+off]; connect[17+off] = temp; 296 temp = connect[18+off]; connect[18+off] = connect[19+off]; connect[19+off] = temp; 297 298 temp = connect[12+off]; connect[12+off] = connect[16+off]; connect[16+off] = temp; 299 temp = connect[13+off]; connect[13+off] = connect[17+off]; connect[17+off] = temp; 300 temp = connect[14+off]; connect[14+off] = connect[18+off]; connect[18+off] = temp; 301 temp = connect[15+off]; connect[15+off] = connect[19+off]; connect[19+off] = temp; 302 303 temp = connect[23+off]; connect[23+off] = connect[26+off]; connect[26+off] = temp; 304 temp = connect[24+off]; connect[24+off] = connect[25+off]; connect[25+off] = temp; 305 temp = connect[25+off]; connect[25+off] = connect[26+off]; connect[26+off] = temp; 306 } 307 } 308 off += connectSize; 309 ierr = DMPlexRestoreTransitiveClosure(dm, cells[c], PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 310 } 311 PetscStackCallStandard(ex_put_conn,(exoid, EX_ELEM_BLOCK, csIdx[cs], connect, 0, 0)); 312 skipCells += (nodes[cs][3] == 0)*csSize; 313 ierr = PetscFree(connect);CHKERRQ(ierr); 314 ierr = ISRestoreIndices(stratumIS, &cells);CHKERRQ(ierr); 315 ierr = ISDestroy(&stratumIS);CHKERRQ(ierr); 316 } 317 ierr = PetscFree(type);CHKERRQ(ierr); 318 /* --- Coordinates --- */ 319 ierr = PetscSectionCreate(comm, §ion);CHKERRQ(ierr); 320 ierr = PetscSectionSetChart(section, pStart, pEnd);CHKERRQ(ierr); 321 for (d = 0; d < depth; ++d) { 322 ierr = DMPlexGetDepthStratum(dm, d, &pStart, &pEnd);CHKERRQ(ierr); 323 for (p = pStart; p < pEnd; ++p) { 324 ierr = PetscSectionSetDof(section, p, nodes[0][d] > 0);CHKERRQ(ierr); 325 } 326 } 327 for (cs = 0; cs < num_cs; ++cs) { 328 IS stratumIS; 329 const PetscInt *cells; 330 PetscInt csSize, c; 331 332 ierr = DMLabelGetStratumIS(csLabel, csIdx[cs], &stratumIS);CHKERRQ(ierr); 333 ierr = ISGetIndices(stratumIS, &cells);CHKERRQ(ierr); 334 ierr = ISGetSize(stratumIS, &csSize);CHKERRQ(ierr); 335 for (c = 0; c < csSize; ++c) { 336 ierr = PetscSectionSetDof(section, cells[c], nodes[cs][3] > 0);CHKERRQ(ierr); 337 } 338 ierr = ISRestoreIndices(stratumIS, &cells);CHKERRQ(ierr); 339 ierr = ISDestroy(&stratumIS);CHKERRQ(ierr); 340 } 341 if (num_cs > 0) { 342 ierr = ISRestoreIndices(csIS, &csIdx);CHKERRQ(ierr); 343 ierr = ISDestroy(&csIS);CHKERRQ(ierr); 344 } 345 ierr = PetscFree(nodes);CHKERRQ(ierr); 346 ierr = PetscSectionSetUp(section);CHKERRQ(ierr); 347 if (numNodes > 0) { 348 const char *coordNames[3] = {"x", "y", "z"}; 349 PetscScalar *coords, *closure; 350 PetscReal *cval; 351 PetscInt hasDof, n = 0; 352 353 /* There can't be more than 24 values in the closure of a point for the coord section */ 354 ierr = PetscCalloc3(numNodes*3, &coords, dim, &cval, 24, &closure);CHKERRQ(ierr); 355 ierr = DMGetCoordinatesLocal(dm, &coord);CHKERRQ(ierr); 356 ierr = DMPlexGetChart(dm, &pStart, &pEnd);CHKERRQ(ierr); 357 for (p = pStart; p < pEnd; ++p) { 358 ierr = PetscSectionGetDof(section, p, &hasDof); 359 if (hasDof) { 360 PetscInt closureSize = 24, j; 361 362 ierr = DMPlexVecGetClosure(cdm, NULL, coord, p, &closureSize, &closure);CHKERRQ(ierr); 363 for (d = 0; d < dim; ++d) { 364 cval[d] = 0.0; 365 for (j = 0; j < closureSize/dim; j++) cval[d] += closure[j*dim+d]; 366 coords[d*numNodes+n] = cval[d] * dim / closureSize; 367 } 368 ++n; 369 } 370 } 371 PetscStackCallStandard(ex_put_coord,(exoid, &coords[0*numNodes], &coords[1*numNodes], &coords[2*numNodes])); 372 ierr = PetscFree3(coords, cval, closure);CHKERRQ(ierr); 373 PetscStackCallStandard(ex_put_coord_names,(exoid, (char **) coordNames)); 374 } 375 ierr = PetscSectionDestroy(§ion);CHKERRQ(ierr); 376 /* --- Node Sets/Vertex Sets --- */ 377 DMHasLabel(dm, "Vertex Sets", &hasLabel); 378 if (hasLabel) { 379 PetscInt i, vs, vsSize; 380 const PetscInt *vsIdx, *vertices; 381 PetscInt *nodeList; 382 IS vsIS, stratumIS; 383 DMLabel vsLabel; 384 ierr = DMGetLabel(dm, "Vertex Sets", &vsLabel);CHKERRQ(ierr); 385 ierr = DMLabelGetValueIS(vsLabel, &vsIS);CHKERRQ(ierr); 386 ierr = ISGetIndices(vsIS, &vsIdx);CHKERRQ(ierr); 387 for (vs=0; vs<num_vs; ++vs) { 388 ierr = DMLabelGetStratumIS(vsLabel, vsIdx[vs], &stratumIS);CHKERRQ(ierr); 389 ierr = ISGetIndices(stratumIS, &vertices);CHKERRQ(ierr); 390 ierr = ISGetSize(stratumIS, &vsSize);CHKERRQ(ierr); 391 ierr = PetscMalloc1(vsSize, &nodeList); 392 for (i=0; i<vsSize; ++i) { 393 nodeList[i] = vertices[i] - skipCells + 1; 394 } 395 PetscStackCallStandard(ex_put_set_param,(exoid, EX_NODE_SET, vsIdx[vs], vsSize, 0)); 396 PetscStackCallStandard(ex_put_set,(exoid, EX_NODE_SET, vsIdx[vs], nodeList, NULL)); 397 ierr = ISRestoreIndices(stratumIS, &vertices);CHKERRQ(ierr); 398 ierr = ISDestroy(&stratumIS);CHKERRQ(ierr); 399 ierr = PetscFree(nodeList);CHKERRQ(ierr); 400 } 401 ierr = ISRestoreIndices(vsIS, &vsIdx);CHKERRQ(ierr); 402 ierr = ISDestroy(&vsIS);CHKERRQ(ierr); 403 } 404 /* --- Side Sets/Face Sets --- */ 405 ierr = DMHasLabel(dm, "Face Sets", &hasLabel);CHKERRQ(ierr); 406 if (hasLabel) { 407 PetscInt i, j, fs, fsSize; 408 const PetscInt *fsIdx, *faces; 409 IS fsIS, stratumIS; 410 DMLabel fsLabel; 411 PetscInt numPoints, *points; 412 PetscInt elem_list_size = 0; 413 PetscInt *elem_list, *elem_ind, *side_list; 414 415 ierr = DMGetLabel(dm, "Face Sets", &fsLabel);CHKERRQ(ierr); 416 /* Compute size of Node List and Element List */ 417 ierr = DMLabelGetValueIS(fsLabel, &fsIS);CHKERRQ(ierr); 418 ierr = ISGetIndices(fsIS, &fsIdx);CHKERRQ(ierr); 419 for (fs=0; fs<num_fs; ++fs) { 420 ierr = DMLabelGetStratumIS(fsLabel, fsIdx[fs], &stratumIS);CHKERRQ(ierr);; 421 ierr = ISGetSize(stratumIS, &fsSize);CHKERRQ(ierr); 422 elem_list_size += fsSize; 423 ierr = ISDestroy(&stratumIS);CHKERRQ(ierr); 424 } 425 ierr = PetscMalloc3(num_fs, &elem_ind, 426 elem_list_size, &elem_list, 427 elem_list_size, &side_list);CHKERRQ(ierr); 428 elem_ind[0] = 0; 429 for (fs=0; fs<num_fs; ++fs) { 430 ierr = DMLabelGetStratumIS(fsLabel, fsIdx[fs], &stratumIS);CHKERRQ(ierr); 431 ierr = ISGetIndices(stratumIS, &faces);CHKERRQ(ierr); 432 ierr = ISGetSize(stratumIS, &fsSize);CHKERRQ(ierr); 433 /* Set Parameters */ 434 PetscStackCallStandard(ex_put_set_param,(exoid, EX_SIDE_SET, fsIdx[fs], fsSize, 0)); 435 /* Indices */ 436 if (fs<num_fs-1) { 437 elem_ind[fs+1] = elem_ind[fs] + fsSize; 438 } 439 440 for (i=0; i<fsSize; ++i) { 441 /* Element List */ 442 points = NULL; 443 ierr = DMPlexGetTransitiveClosure(dm, faces[i], PETSC_FALSE, &numPoints, &points);CHKERRQ(ierr); 444 elem_list[elem_ind[fs] + i] = points[2] +1; 445 ierr = DMPlexRestoreTransitiveClosure(dm, faces[i], PETSC_FALSE, &numPoints, &points);CHKERRQ(ierr); 446 447 /* Side List */ 448 points = NULL; 449 ierr = DMPlexGetTransitiveClosure(dm, elem_list[elem_ind[fs] + i]-1, PETSC_TRUE, &numPoints, &points);CHKERRQ(ierr); 450 for (j=1; j<numPoints; ++j) { 451 if (points[j*2]==faces[i]) {break;} 452 } 453 /* Convert HEX sides */ 454 if (numPoints == 27) { 455 if (j == 1) {j = 5;} 456 else if (j == 2) {j = 6;} 457 else if (j == 3) {j = 1;} 458 else if (j == 4) {j = 3;} 459 else if (j == 5) {j = 2;} 460 else if (j == 6) {j = 4;} 461 } 462 /* Convert TET sides */ 463 if (numPoints == 15) { 464 --j; 465 if (j == 0) {j = 4;} 466 } 467 side_list[elem_ind[fs] + i] = j; 468 ierr = DMPlexRestoreTransitiveClosure(dm, elem_list[elem_ind[fs] + i]-1, PETSC_TRUE, &numPoints, &points);CHKERRQ(ierr); 469 470 } 471 ierr = ISRestoreIndices(stratumIS, &faces);CHKERRQ(ierr); 472 ierr = ISDestroy(&stratumIS);CHKERRQ(ierr); 473 } 474 ierr = ISRestoreIndices(fsIS, &fsIdx);CHKERRQ(ierr); 475 ierr = ISDestroy(&fsIS);CHKERRQ(ierr); 476 477 /* Put side sets */ 478 for (fs=0; fs<num_fs; ++fs) { 479 PetscStackCallStandard(ex_put_set,(exoid, EX_SIDE_SET, fsIdx[fs], &elem_list[elem_ind[fs]], &side_list[elem_ind[fs]])); 480 } 481 ierr = PetscFree3(elem_ind,elem_list,side_list);CHKERRQ(ierr); 482 } 483 PetscFunctionReturn(0); 484 } 485 486 /* 487 VecViewPlex_ExodusII_Nodal_Internal - Write a Vec corresponding to a nodal field to an exodus file 488 489 Collective on v 490 491 Input Parameters: 492 + v - The vector to be written 493 . exoid - the exodus id of a file (obtained from ex_open or ex_create for instance) 494 - step - the time step to write at (exodus steps are numbered starting from 1) 495 496 Notes: 497 The exodus result nodal variable index is obtained by comparing the Vec name and the 498 names of nodal variables declared in the exodus file. For instance for a Vec named "V" 499 the location in the exodus file will be the first match of "V", "V_X", "V_XX", "V_1", or "V_11" 500 amongst all nodal variables. 501 502 Level: beginner 503 504 .keywords: mesh,ExodusII 505 .seealso: EXOGetVarIndex_Private(),DMPlexView_ExodusII_Internal(),VecLoadNodal_PlexEXO(),VecViewZonal_PlexEXO(),VecLoadZonal_PlexEXO() 506 @*/ 507 PetscErrorCode VecViewPlex_ExodusII_Nodal_Internal(Vec v, int exoid, int step) 508 { 509 MPI_Comm comm; 510 PetscMPIInt size; 511 DM dm; 512 Vec vNatural, vComp; 513 const PetscReal *varray; 514 const char *vecname; 515 PetscInt xs, xe, bs; 516 PetscBool useNatural; 517 int offset; 518 PetscErrorCode ierr; 519 520 PetscFunctionBegin; 521 ierr = PetscObjectGetComm((PetscObject) v, &comm);CHKERRQ(ierr); 522 ierr = MPI_Comm_size(comm, &size);CHKERRQ(ierr); 523 ierr = VecGetDM(v, &dm);CHKERRQ(ierr); 524 ierr = DMGetUseNatural(dm, &useNatural);CHKERRQ(ierr); 525 useNatural = useNatural && size > 1 ? PETSC_TRUE : PETSC_FALSE; 526 if (useNatural) { 527 ierr = DMGetGlobalVector(dm, &vNatural);CHKERRQ(ierr); 528 ierr = DMPlexGlobalToNaturalBegin(dm, v, vNatural);CHKERRQ(ierr); 529 ierr = DMPlexGlobalToNaturalEnd(dm, v, vNatural);CHKERRQ(ierr); 530 } else { 531 vNatural = v; 532 } 533 /* Get the location of the variable in the exodus file */ 534 ierr = PetscObjectGetName((PetscObject) v, &vecname);CHKERRQ(ierr); 535 ierr = EXOGetVarIndex_Private(exoid, EX_NODAL, vecname, &offset);CHKERRQ(ierr); 536 if (!offset) SETERRQ1(comm, PETSC_ERR_FILE_UNEXPECTED, "Cannot locate nodal variable %s in exodus file\n", vecname); 537 /* Write local chunk of the result in the exodus file 538 exodus stores each component of a vector-valued field as a separate variable. 539 We assume that they are stored sequentially */ 540 ierr = VecGetOwnershipRange(vNatural, &xs, &xe);CHKERRQ(ierr); 541 ierr = VecGetBlockSize(vNatural, &bs);CHKERRQ(ierr); 542 if (bs == 1) { 543 ierr = VecGetArrayRead(vNatural, &varray);CHKERRQ(ierr); 544 PetscStackCallStandard(ex_put_partial_var,(exoid, step, EX_NODAL, offset, 1, xs+1, xe-xs, varray)); 545 ierr = VecRestoreArrayRead(vNatural, &varray);CHKERRQ(ierr); 546 } else { 547 IS compIS; 548 PetscInt c; 549 550 ierr = ISCreateStride(PETSC_COMM_SELF, (xe-xs)/bs, xs, bs, &compIS);CHKERRQ(ierr); 551 for (c = 0; c < bs; ++c) { 552 ierr = ISStrideSetStride(compIS, (xe-xs)/bs, xs+c, bs);CHKERRQ(ierr); 553 ierr = VecGetSubVector(vNatural, compIS, &vComp);CHKERRQ(ierr); 554 ierr = VecGetArrayRead(vComp, &varray);CHKERRQ(ierr); 555 PetscStackCallStandard(ex_put_partial_var,(exoid, step, EX_NODAL, offset+c, 1, xs/bs+1, (xe-xs)/bs, varray)); 556 ierr = VecRestoreArrayRead(vComp, &varray);CHKERRQ(ierr); 557 ierr = VecRestoreSubVector(vNatural, compIS, &vComp);CHKERRQ(ierr); 558 } 559 ierr = ISDestroy(&compIS);CHKERRQ(ierr); 560 } 561 if (useNatural) {ierr = DMRestoreGlobalVector(dm, &vNatural);CHKERRQ(ierr);} 562 PetscFunctionReturn(0); 563 } 564 565 /* 566 VecLoadPlex_ExodusII_Nodal_Internal - Read a Vec corresponding to a nodal field from an exodus file 567 568 Collective on v 569 570 Input Parameters: 571 + v - The vector to be written 572 . exoid - the exodus id of a file (obtained from ex_open or ex_create for instance) 573 - step - the time step to read at (exodus steps are numbered starting from 1) 574 575 Notes: 576 The exodus result nodal variable index is obtained by comparing the Vec name and the 577 names of nodal variables declared in the exodus file. For instance for a Vec named "V" 578 the location in the exodus file will be the first match of "V", "V_X", "V_XX", "V_1", or "V_11" 579 amongst all nodal variables. 580 581 Level: beginner 582 583 .keywords: mesh, ExodusII 584 .seealso: EXOGetVarIndex_Private(), DMPlexView_ExodusII_Internal(), VecViewPlex_ExodusII_Nodal_Internal(), VecViewZonal_PlexEXO(), VecLoadZonal_PlexEXO() 585 */ 586 PetscErrorCode VecLoadPlex_ExodusII_Nodal_Internal(Vec v, int exoid, int step) 587 { 588 MPI_Comm comm; 589 PetscMPIInt size; 590 DM dm; 591 Vec vNatural, vComp; 592 PetscReal *varray; 593 const char *vecname; 594 PetscInt xs, xe, bs; 595 PetscBool useNatural; 596 int offset; 597 PetscErrorCode ierr; 598 599 PetscFunctionBegin; 600 ierr = PetscObjectGetComm((PetscObject) v, &comm);CHKERRQ(ierr); 601 ierr = MPI_Comm_size(comm, &size);CHKERRQ(ierr); 602 ierr = VecGetDM(v,&dm);CHKERRQ(ierr); 603 ierr = DMGetUseNatural(dm, &useNatural);CHKERRQ(ierr); 604 useNatural = useNatural && size > 1 ? PETSC_TRUE : PETSC_FALSE; 605 if (useNatural) {ierr = DMGetGlobalVector(dm,&vNatural);CHKERRQ(ierr);} 606 else {vNatural = v;} 607 /* Get the location of the variable in the exodus file */ 608 ierr = PetscObjectGetName((PetscObject) v, &vecname);CHKERRQ(ierr); 609 ierr = EXOGetVarIndex_Private(exoid, EX_NODAL, vecname, &offset);CHKERRQ(ierr); 610 if (!offset) SETERRQ1(comm, PETSC_ERR_FILE_UNEXPECTED, "Cannot locate nodal variable %s in exodus file", vecname); 611 /* Read local chunk from the file */ 612 ierr = VecGetOwnershipRange(vNatural, &xs, &xe);CHKERRQ(ierr); 613 ierr = VecGetBlockSize(vNatural, &bs);CHKERRQ(ierr); 614 if (bs == 1) { 615 ierr = VecGetArray(vNatural, &varray);CHKERRQ(ierr); 616 PetscStackCallStandard(ex_get_partial_var,(exoid, step, EX_NODAL, offset, 1, xs+1, xe-xs, varray)); 617 ierr = VecRestoreArray(vNatural, &varray);CHKERRQ(ierr); 618 } else { 619 IS compIS; 620 PetscInt c; 621 622 ierr = ISCreateStride(PETSC_COMM_SELF, (xe-xs)/bs, xs, bs, &compIS);CHKERRQ(ierr); 623 for (c = 0; c < bs; ++c) { 624 ierr = ISStrideSetStride(compIS, (xe-xs)/bs, xs+c, bs);CHKERRQ(ierr); 625 ierr = VecGetSubVector(vNatural, compIS, &vComp);CHKERRQ(ierr); 626 ierr = VecGetArray(vComp, &varray);CHKERRQ(ierr); 627 PetscStackCallStandard(ex_get_partial_var,(exoid, step, EX_NODAL, offset+c, 1, xs/bs+1, (xe-xs)/bs, varray)); 628 ierr = VecRestoreArray(vComp, &varray);CHKERRQ(ierr); 629 ierr = VecRestoreSubVector(vNatural, compIS, &vComp);CHKERRQ(ierr); 630 } 631 ierr = ISDestroy(&compIS);CHKERRQ(ierr); 632 } 633 if (useNatural) { 634 ierr = DMPlexNaturalToGlobalBegin(dm, vNatural, v);CHKERRQ(ierr); 635 ierr = DMPlexNaturalToGlobalEnd(dm, vNatural, v);CHKERRQ(ierr); 636 ierr = DMRestoreGlobalVector(dm, &vNatural);CHKERRQ(ierr); 637 } 638 PetscFunctionReturn(0); 639 } 640 641 /* 642 VecViewPlex_ExodusII_Zonal_Internal - Write a Vec corresponding to a zonal (cell based) field to an exodus file 643 644 Collective on v 645 646 Input Parameters: 647 + v - The vector to be written 648 . exoid - the exodus id of a file (obtained from ex_open or ex_create for instance) 649 - step - the time step to write at (exodus steps are numbered starting from 1) 650 651 Notes: 652 The exodus result zonal variable index is obtained by comparing the Vec name and the 653 names of zonal variables declared in the exodus file. For instance for a Vec named "V" 654 the location in the exodus file will be the first match of "V", "V_X", "V_XX", "V_1", or "V_11" 655 amongst all zonal variables. 656 657 Level: beginner 658 659 .keywords: mesh,ExodusII 660 .seealso: EXOGetVarIndex_Private(),DMPlexView_ExodusII_Internal(),VecViewPlex_ExodusII_Nodal_Internal(),VecLoadPlex_ExodusII_Nodal_Internal(),VecLoadPlex_ExodusII_Zonal_Internal() 661 */ 662 PetscErrorCode VecViewPlex_ExodusII_Zonal_Internal(Vec v, int exoid, int step) 663 { 664 MPI_Comm comm; 665 PetscMPIInt size; 666 DM dm; 667 Vec vNatural, vComp; 668 const PetscReal *varray; 669 const char *vecname; 670 PetscInt xs, xe, bs; 671 PetscBool useNatural; 672 int offset; 673 IS compIS; 674 PetscInt *csSize, *csID; 675 PetscInt numCS, set, csxs = 0; 676 PetscErrorCode ierr; 677 678 PetscFunctionBegin; 679 ierr = PetscObjectGetComm((PetscObject)v, &comm);CHKERRQ(ierr); 680 ierr = MPI_Comm_size(comm, &size);CHKERRQ(ierr); 681 ierr = VecGetDM(v, &dm);CHKERRQ(ierr); 682 ierr = DMGetUseNatural(dm, &useNatural);CHKERRQ(ierr); 683 useNatural = useNatural && size > 1 ? PETSC_TRUE : PETSC_FALSE; 684 if (useNatural) { 685 ierr = DMGetGlobalVector(dm, &vNatural);CHKERRQ(ierr); 686 ierr = DMPlexGlobalToNaturalBegin(dm, v, vNatural);CHKERRQ(ierr); 687 ierr = DMPlexGlobalToNaturalEnd(dm, v, vNatural);CHKERRQ(ierr); 688 } else { 689 vNatural = v; 690 } 691 /* Get the location of the variable in the exodus file */ 692 ierr = PetscObjectGetName((PetscObject) v, &vecname);CHKERRQ(ierr); 693 ierr = EXOGetVarIndex_Private(exoid, EX_ELEM_BLOCK, vecname, &offset);CHKERRQ(ierr); 694 if (!offset) SETERRQ1(comm, PETSC_ERR_FILE_UNEXPECTED, "Cannot locate zonal variable %s in exodus file\n", vecname); 695 /* Write local chunk of the result in the exodus file 696 exodus stores each component of a vector-valued field as a separate variable. 697 We assume that they are stored sequentially 698 Zonal variables are accessed one element block at a time, so we loop through the cell sets, 699 but once the vector has been reordered to natural size, we cannot use the label informations 700 to figure out what to save where. */ 701 numCS = ex_inquire_int(exoid, EX_INQ_ELEM_BLK); 702 ierr = PetscMalloc2(numCS, &csID, numCS, &csSize);CHKERRQ(ierr); 703 PetscStackCallStandard(ex_get_ids,(exoid, EX_ELEM_BLOCK, csID)); 704 for (set = 0; set < numCS; ++set) { 705 ex_block block; 706 707 block.id = csID[set]; 708 block.type = EX_ELEM_BLOCK; 709 PetscStackCallStandard(ex_get_block_param,(exoid, &block)); 710 csSize[set] = block.num_entry; 711 } 712 ierr = VecGetOwnershipRange(vNatural, &xs, &xe);CHKERRQ(ierr); 713 ierr = VecGetBlockSize(vNatural, &bs);CHKERRQ(ierr); 714 if (bs > 1) {ierr = ISCreateStride(comm, (xe-xs)/bs, xs, bs, &compIS);CHKERRQ(ierr);} 715 for (set = 0; set < numCS; set++) { 716 PetscInt csLocalSize, c; 717 718 /* range of indices for set setID[set]: csxs:csxs + csSize[set]-1 719 local slice of zonal values: xs/bs,xm/bs-1 720 intersection: max(xs/bs,csxs),min(xm/bs-1,csxs + csSize[set]-1) */ 721 csLocalSize = PetscMax(0, PetscMin(xe/bs, csxs+csSize[set]) - PetscMax(xs/bs, csxs)); 722 if (bs == 1) { 723 ierr = VecGetArrayRead(vNatural, &varray);CHKERRQ(ierr); 724 PetscStackCallStandard(ex_put_partial_var,(exoid, step, EX_ELEM_BLOCK, offset, csID[set], PetscMax(xs-csxs, 0)+1, csLocalSize, &varray[PetscMax(0, csxs-xs)])); 725 ierr = VecRestoreArrayRead(vNatural, &varray);CHKERRQ(ierr); 726 } else { 727 for (c = 0; c < bs; ++c) { 728 ierr = ISStrideSetStride(compIS, (xe-xs)/bs, xs+c, bs);CHKERRQ(ierr); 729 ierr = VecGetSubVector(vNatural, compIS, &vComp);CHKERRQ(ierr); 730 ierr = VecGetArrayRead(vComp, &varray);CHKERRQ(ierr); 731 PetscStackCallStandard(ex_put_partial_var,(exoid, step, EX_ELEM_BLOCK, offset+c, csID[set], PetscMax(xs/bs-csxs, 0)+1, csLocalSize, &varray[PetscMax(0, csxs-xs/bs)])); 732 ierr = VecRestoreArrayRead(vComp, &varray);CHKERRQ(ierr); 733 ierr = VecRestoreSubVector(vNatural, compIS, &vComp);CHKERRQ(ierr); 734 } 735 } 736 csxs += csSize[set]; 737 } 738 ierr = PetscFree2(csID, csSize);CHKERRQ(ierr); 739 if (bs > 1) {ierr = ISDestroy(&compIS);CHKERRQ(ierr);} 740 if (useNatural) {ierr = DMRestoreGlobalVector(dm,&vNatural);CHKERRQ(ierr);} 741 PetscFunctionReturn(0); 742 } 743 744 /* 745 VecLoadPlex_ExodusII_Zonal_Internal - Read a Vec corresponding to a zonal (cell based) field from an exodus file 746 747 Collective on v 748 749 Input Parameters: 750 + v - The vector to be written 751 . exoid - the exodus id of a file (obtained from ex_open or ex_create for instance) 752 - step - the time step to read at (exodus steps are numbered starting from 1) 753 754 Notes: 755 The exodus result zonal variable index is obtained by comparing the Vec name and the 756 names of zonal variables declared in the exodus file. For instance for a Vec named "V" 757 the location in the exodus file will be the first match of "V", "V_X", "V_XX", "V_1", or "V_11" 758 amongst all zonal variables. 759 760 Level: beginner 761 762 .keywords: mesh,ExodusII 763 .seealso: EXOGetVarIndex_Private(), DMPlexView_ExodusII_Internal(), VecViewPlex_ExodusII_Nodal_Internal(), VecLoadPlex_ExodusII_Nodal_Internal(), VecLoadPlex_ExodusII_Zonal_Internal() 764 */ 765 PetscErrorCode VecLoadPlex_ExodusII_Zonal_Internal(Vec v, int exoid, int step) 766 { 767 MPI_Comm comm; 768 PetscMPIInt size; 769 DM dm; 770 Vec vNatural, vComp; 771 PetscReal *varray; 772 const char *vecname; 773 PetscInt xs, xe, bs; 774 PetscBool useNatural; 775 int offset; 776 IS compIS; 777 PetscInt *csSize, *csID; 778 PetscInt numCS, set, csxs = 0; 779 PetscErrorCode ierr; 780 781 PetscFunctionBegin; 782 ierr = PetscObjectGetComm((PetscObject)v,&comm);CHKERRQ(ierr); 783 ierr = MPI_Comm_size(comm, &size);CHKERRQ(ierr); 784 ierr = VecGetDM(v, &dm);CHKERRQ(ierr); 785 ierr = DMGetUseNatural(dm, &useNatural);CHKERRQ(ierr); 786 useNatural = useNatural && size > 1 ? PETSC_TRUE : PETSC_FALSE; 787 if (useNatural) {ierr = DMGetGlobalVector(dm,&vNatural);CHKERRQ(ierr);} 788 else {vNatural = v;} 789 /* Get the location of the variable in the exodus file */ 790 ierr = PetscObjectGetName((PetscObject) v, &vecname);CHKERRQ(ierr); 791 ierr = EXOGetVarIndex_Private(exoid, EX_ELEM_BLOCK, vecname, &offset);CHKERRQ(ierr); 792 if (!offset) SETERRQ1(comm, PETSC_ERR_FILE_UNEXPECTED, "Cannot locate zonal variable %s in exodus file\n", vecname); 793 /* Read local chunk of the result in the exodus file 794 exodus stores each component of a vector-valued field as a separate variable. 795 We assume that they are stored sequentially 796 Zonal variables are accessed one element block at a time, so we loop through the cell sets, 797 but once the vector has been reordered to natural size, we cannot use the label informations 798 to figure out what to save where. */ 799 numCS = ex_inquire_int(exoid, EX_INQ_ELEM_BLK); 800 ierr = PetscMalloc2(numCS, &csID, numCS, &csSize);CHKERRQ(ierr); 801 PetscStackCallStandard(ex_get_ids,(exoid, EX_ELEM_BLOCK, csID)); 802 for (set = 0; set < numCS; ++set) { 803 ex_block block; 804 805 block.id = csID[set]; 806 block.type = EX_ELEM_BLOCK; 807 PetscStackCallStandard(ex_get_block_param,(exoid, &block)); 808 csSize[set] = block.num_entry; 809 } 810 ierr = VecGetOwnershipRange(vNatural, &xs, &xe);CHKERRQ(ierr); 811 ierr = VecGetBlockSize(vNatural, &bs);CHKERRQ(ierr); 812 if (bs > 1) {ierr = ISCreateStride(comm, (xe-xs)/bs, xs, bs, &compIS);CHKERRQ(ierr);} 813 for (set = 0; set < numCS; ++set) { 814 PetscInt csLocalSize, c; 815 816 /* range of indices for set setID[set]: csxs:csxs + csSize[set]-1 817 local slice of zonal values: xs/bs,xm/bs-1 818 intersection: max(xs/bs,csxs),min(xm/bs-1,csxs + csSize[set]-1) */ 819 csLocalSize = PetscMax(0, PetscMin(xe/bs, csxs+csSize[set]) - PetscMax(xs/bs, csxs)); 820 if (bs == 1) { 821 ierr = VecGetArray(vNatural, &varray);CHKERRQ(ierr); 822 PetscStackCallStandard(ex_get_partial_var,(exoid, step, EX_ELEM_BLOCK, offset, csID[set], PetscMax(xs-csxs, 0)+1, csLocalSize, &varray[PetscMax(0, csxs-xs)])); 823 ierr = VecRestoreArray(vNatural, &varray);CHKERRQ(ierr); 824 } else { 825 for (c = 0; c < bs; ++c) { 826 ierr = ISStrideSetStride(compIS, (xe-xs)/bs, xs+c, bs);CHKERRQ(ierr); 827 ierr = VecGetSubVector(vNatural, compIS, &vComp);CHKERRQ(ierr); 828 ierr = VecGetArray(vComp, &varray);CHKERRQ(ierr); 829 PetscStackCallStandard(ex_get_partial_var,(exoid, step, EX_ELEM_BLOCK, offset+c, csID[set], PetscMax(xs/bs-csxs, 0)+1, csLocalSize, &varray[PetscMax(0, csxs-xs/bs)])); 830 ierr = VecRestoreArray(vComp, &varray);CHKERRQ(ierr); 831 ierr = VecRestoreSubVector(vNatural, compIS, &vComp);CHKERRQ(ierr); 832 } 833 } 834 csxs += csSize[set]; 835 } 836 ierr = PetscFree2(csID, csSize);CHKERRQ(ierr); 837 if (bs > 1) {ierr = ISDestroy(&compIS);CHKERRQ(ierr);} 838 if (useNatural) { 839 ierr = DMPlexNaturalToGlobalBegin(dm, vNatural, v);CHKERRQ(ierr); 840 ierr = DMPlexNaturalToGlobalEnd(dm, vNatural, v);CHKERRQ(ierr); 841 ierr = DMRestoreGlobalVector(dm, &vNatural);CHKERRQ(ierr); 842 } 843 PetscFunctionReturn(0); 844 } 845 #endif 846 847 /*@C 848 DMPlexCreateExodusFromFile - Create a DMPlex mesh from an ExodusII file. 849 850 Collective on comm 851 852 Input Parameters: 853 + comm - The MPI communicator 854 . filename - The name of the ExodusII file 855 - interpolate - Create faces and edges in the mesh 856 857 Output Parameter: 858 . dm - The DM object representing the mesh 859 860 Level: beginner 861 862 .keywords: mesh,ExodusII 863 .seealso: DMPLEX, DMCreate(), DMPlexCreateExodus() 864 @*/ 865 PetscErrorCode DMPlexCreateExodusFromFile(MPI_Comm comm, const char filename[], PetscBool interpolate, DM *dm) 866 { 867 PetscMPIInt rank; 868 PetscErrorCode ierr; 869 #if defined(PETSC_HAVE_EXODUSII) 870 int CPU_word_size = sizeof(PetscReal), IO_word_size = 0, exoid = -1; 871 float version; 872 #endif 873 874 PetscFunctionBegin; 875 PetscValidCharPointer(filename, 2); 876 ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr); 877 #if defined(PETSC_HAVE_EXODUSII) 878 if (!rank) { 879 exoid = ex_open(filename, EX_READ, &CPU_word_size, &IO_word_size, &version); 880 if (exoid <= 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_LIB, "ex_open(\"%s\",...) did not return a valid file ID", filename); 881 } 882 ierr = DMPlexCreateExodus(comm, exoid, interpolate, dm);CHKERRQ(ierr); 883 if (!rank) {PetscStackCallStandard(ex_close,(exoid));} 884 #else 885 SETERRQ(comm, PETSC_ERR_SUP, "This method requires ExodusII support. Reconfigure using --download-exodusii"); 886 #endif 887 PetscFunctionReturn(0); 888 } 889 890 /*@ 891 DMPlexCreateExodus - Create a DMPlex mesh from an ExodusII file ID. 892 893 Collective on comm 894 895 Input Parameters: 896 + comm - The MPI communicator 897 . exoid - The ExodusII id associated with a exodus file and obtained using ex_open 898 - interpolate - Create faces and edges in the mesh 899 900 Output Parameter: 901 . dm - The DM object representing the mesh 902 903 Level: beginner 904 905 .keywords: mesh,ExodusII 906 .seealso: DMPLEX, DMCreate() 907 @*/ 908 PetscErrorCode DMPlexCreateExodus(MPI_Comm comm, PetscInt exoid, PetscBool interpolate, DM *dm) 909 { 910 #if defined(PETSC_HAVE_EXODUSII) 911 PetscMPIInt num_proc, rank; 912 PetscSection coordSection; 913 Vec coordinates; 914 PetscScalar *coords; 915 PetscInt coordSize, v; 916 PetscErrorCode ierr; 917 /* Read from ex_get_init() */ 918 char title[PETSC_MAX_PATH_LEN+1]; 919 int dim = 0, numVertices = 0, numCells = 0; 920 int num_cs = 0, num_vs = 0, num_fs = 0; 921 #endif 922 923 PetscFunctionBegin; 924 #if defined(PETSC_HAVE_EXODUSII) 925 ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr); 926 ierr = MPI_Comm_size(comm, &num_proc);CHKERRQ(ierr); 927 ierr = DMCreate(comm, dm);CHKERRQ(ierr); 928 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr); 929 /* Open EXODUS II file and read basic informations on rank 0, then broadcast to all processors */ 930 if (!rank) { 931 ierr = PetscMemzero(title,(PETSC_MAX_PATH_LEN+1)*sizeof(char));CHKERRQ(ierr); 932 PetscStackCallStandard(ex_get_init,(exoid, title, &dim, &numVertices, &numCells, &num_cs, &num_vs, &num_fs)); 933 if (!num_cs) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Exodus file does not contain any cell set\n"); 934 } 935 ierr = MPI_Bcast(title, PETSC_MAX_PATH_LEN+1, MPI_CHAR, 0, comm);CHKERRQ(ierr); 936 ierr = MPI_Bcast(&dim, 1, MPI_INT, 0, comm);CHKERRQ(ierr); 937 ierr = PetscObjectSetName((PetscObject) *dm, title);CHKERRQ(ierr); 938 ierr = DMSetDimension(*dm, dim);CHKERRQ(ierr); 939 ierr = DMPlexSetChart(*dm, 0, numCells+numVertices);CHKERRQ(ierr); 940 941 /* Read cell sets information */ 942 if (!rank) { 943 PetscInt *cone; 944 int c, cs, c_loc, v, v_loc; 945 /* Read from ex_get_elem_blk_ids() */ 946 int *cs_id; 947 /* Read from ex_get_elem_block() */ 948 char buffer[PETSC_MAX_PATH_LEN+1]; 949 int num_cell_in_set, num_vertex_per_cell, num_attr; 950 /* Read from ex_get_elem_conn() */ 951 int *cs_connect; 952 953 /* Get cell sets IDs */ 954 ierr = PetscMalloc1(num_cs, &cs_id);CHKERRQ(ierr); 955 PetscStackCallStandard(ex_get_ids,(exoid, EX_ELEM_BLOCK, cs_id)); 956 /* Read the cell set connectivity table and build mesh topology 957 EXO standard requires that cells in cell sets be numbered sequentially and be pairwise disjoint. */ 958 /* First set sizes */ 959 for (cs = 0, c = 0; cs < num_cs; cs++) { 960 PetscStackCallStandard(ex_get_block,(exoid, EX_ELEM_BLOCK, cs_id[cs], buffer, &num_cell_in_set,&num_vertex_per_cell, 0, 0, &num_attr)); 961 for (c_loc = 0; c_loc < num_cell_in_set; ++c_loc, ++c) { 962 ierr = DMPlexSetConeSize(*dm, c, num_vertex_per_cell);CHKERRQ(ierr); 963 } 964 } 965 ierr = DMSetUp(*dm);CHKERRQ(ierr); 966 for (cs = 0, c = 0; cs < num_cs; cs++) { 967 PetscStackCallStandard(ex_get_block,(exoid, EX_ELEM_BLOCK, cs_id[cs], buffer, &num_cell_in_set, &num_vertex_per_cell, 0, 0, &num_attr)); 968 ierr = PetscMalloc2(num_vertex_per_cell*num_cell_in_set,&cs_connect,num_vertex_per_cell,&cone);CHKERRQ(ierr); 969 PetscStackCallStandard(ex_get_conn,(exoid, EX_ELEM_BLOCK, cs_id[cs], cs_connect,NULL,NULL)); 970 /* EXO uses Fortran-based indexing, DMPlex uses C-style and numbers cell first then vertices. */ 971 for (c_loc = 0, v = 0; c_loc < num_cell_in_set; ++c_loc, ++c) { 972 for (v_loc = 0; v_loc < num_vertex_per_cell; ++v_loc, ++v) { 973 cone[v_loc] = cs_connect[v]+numCells-1; 974 } 975 if (dim == 3) { 976 /* Tetrahedra are inverted */ 977 if (num_vertex_per_cell == 4) { 978 PetscInt tmp = cone[0]; 979 cone[0] = cone[1]; 980 cone[1] = tmp; 981 } 982 /* Hexahedra are inverted */ 983 if (num_vertex_per_cell == 8) { 984 PetscInt tmp = cone[1]; 985 cone[1] = cone[3]; 986 cone[3] = tmp; 987 } 988 } 989 ierr = DMPlexSetCone(*dm, c, cone);CHKERRQ(ierr); 990 ierr = DMSetLabelValue(*dm, "Cell Sets", c, cs_id[cs]);CHKERRQ(ierr); 991 } 992 ierr = PetscFree2(cs_connect,cone);CHKERRQ(ierr); 993 } 994 ierr = PetscFree(cs_id);CHKERRQ(ierr); 995 } 996 ierr = DMPlexSymmetrize(*dm);CHKERRQ(ierr); 997 ierr = DMPlexStratify(*dm);CHKERRQ(ierr); 998 if (interpolate) { 999 DM idm; 1000 1001 ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr); 1002 ierr = DMDestroy(dm);CHKERRQ(ierr); 1003 *dm = idm; 1004 } 1005 1006 /* Create vertex set label */ 1007 if (!rank && (num_vs > 0)) { 1008 int vs, v; 1009 /* Read from ex_get_node_set_ids() */ 1010 int *vs_id; 1011 /* Read from ex_get_node_set_param() */ 1012 int num_vertex_in_set; 1013 /* Read from ex_get_node_set() */ 1014 int *vs_vertex_list; 1015 1016 /* Get vertex set ids */ 1017 ierr = PetscMalloc1(num_vs, &vs_id);CHKERRQ(ierr); 1018 PetscStackCallStandard(ex_get_ids,(exoid, EX_NODE_SET, vs_id)); 1019 for (vs = 0; vs < num_vs; ++vs) { 1020 PetscStackCallStandard(ex_get_set_param,(exoid, EX_NODE_SET, vs_id[vs], &num_vertex_in_set, NULL)); 1021 ierr = PetscMalloc1(num_vertex_in_set, &vs_vertex_list);CHKERRQ(ierr); 1022 PetscStackCallStandard(ex_get_set,(exoid, EX_NODE_SET, vs_id[vs], vs_vertex_list, NULL)); 1023 for (v = 0; v < num_vertex_in_set; ++v) { 1024 ierr = DMSetLabelValue(*dm, "Vertex Sets", vs_vertex_list[v]+numCells-1, vs_id[vs]);CHKERRQ(ierr); 1025 } 1026 ierr = PetscFree(vs_vertex_list);CHKERRQ(ierr); 1027 } 1028 ierr = PetscFree(vs_id);CHKERRQ(ierr); 1029 } 1030 /* Read coordinates */ 1031 ierr = DMGetCoordinateSection(*dm, &coordSection);CHKERRQ(ierr); 1032 ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr); 1033 ierr = PetscSectionSetFieldComponents(coordSection, 0, dim);CHKERRQ(ierr); 1034 ierr = PetscSectionSetChart(coordSection, numCells, numCells + numVertices);CHKERRQ(ierr); 1035 for (v = numCells; v < numCells+numVertices; ++v) { 1036 ierr = PetscSectionSetDof(coordSection, v, dim);CHKERRQ(ierr); 1037 ierr = PetscSectionSetFieldDof(coordSection, v, 0, dim);CHKERRQ(ierr); 1038 } 1039 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 1040 ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr); 1041 ierr = VecCreate(PETSC_COMM_SELF, &coordinates);CHKERRQ(ierr); 1042 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr); 1043 ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 1044 ierr = VecSetBlockSize(coordinates, dim);CHKERRQ(ierr); 1045 ierr = VecSetType(coordinates,VECSTANDARD);CHKERRQ(ierr); 1046 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 1047 if (!rank) { 1048 PetscReal *x, *y, *z; 1049 1050 ierr = PetscMalloc3(numVertices,&x,numVertices,&y,numVertices,&z);CHKERRQ(ierr); 1051 PetscStackCallStandard(ex_get_coord,(exoid, x, y, z)); 1052 if (dim > 0) { 1053 for (v = 0; v < numVertices; ++v) coords[v*dim+0] = x[v]; 1054 } 1055 if (dim > 1) { 1056 for (v = 0; v < numVertices; ++v) coords[v*dim+1] = y[v]; 1057 } 1058 if (dim > 2) { 1059 for (v = 0; v < numVertices; ++v) coords[v*dim+2] = z[v]; 1060 } 1061 ierr = PetscFree3(x,y,z);CHKERRQ(ierr); 1062 } 1063 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 1064 ierr = DMSetCoordinatesLocal(*dm, coordinates);CHKERRQ(ierr); 1065 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 1066 1067 /* Create side set label */ 1068 if (!rank && interpolate && (num_fs > 0)) { 1069 int fs, f, voff; 1070 /* Read from ex_get_side_set_ids() */ 1071 int *fs_id; 1072 /* Read from ex_get_side_set_param() */ 1073 int num_side_in_set; 1074 /* Read from ex_get_side_set_node_list() */ 1075 int *fs_vertex_count_list, *fs_vertex_list; 1076 /* Read side set labels */ 1077 char fs_name[MAX_STR_LENGTH+1]; 1078 1079 /* Get side set ids */ 1080 ierr = PetscMalloc1(num_fs, &fs_id);CHKERRQ(ierr); 1081 PetscStackCallStandard(ex_get_ids,(exoid, EX_SIDE_SET, fs_id)); 1082 for (fs = 0; fs < num_fs; ++fs) { 1083 PetscStackCallStandard(ex_get_set_param,(exoid, EX_SIDE_SET, fs_id[fs], &num_side_in_set, NULL)); 1084 ierr = PetscMalloc2(num_side_in_set,&fs_vertex_count_list,num_side_in_set*4,&fs_vertex_list);CHKERRQ(ierr); 1085 PetscStackCallStandard(ex_get_side_set_node_list,(exoid, fs_id[fs], fs_vertex_count_list, fs_vertex_list)); 1086 /* Get the specific name associated with this side set ID. */ 1087 int fs_name_err = ex_get_name(exoid, EX_SIDE_SET, fs_id[fs], fs_name); 1088 for (f = 0, voff = 0; f < num_side_in_set; ++f) { 1089 const PetscInt *faces = NULL; 1090 PetscInt faceSize = fs_vertex_count_list[f], numFaces; 1091 PetscInt faceVertices[4], v; 1092 1093 if (faceSize > 4) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "ExodusII side cannot have %d > 4 vertices", faceSize); 1094 for (v = 0; v < faceSize; ++v, ++voff) { 1095 faceVertices[v] = fs_vertex_list[voff]+numCells-1; 1096 } 1097 ierr = DMPlexGetFullJoin(*dm, faceSize, faceVertices, &numFaces, &faces);CHKERRQ(ierr); 1098 if (numFaces != 1) SETERRQ3(comm, PETSC_ERR_ARG_WRONG, "Invalid ExodusII side %d in set %d maps to %d faces", f, fs, numFaces); 1099 ierr = DMSetLabelValue(*dm, "Face Sets", faces[0], fs_id[fs]);CHKERRQ(ierr); 1100 /* Only add the label if one has been detected for this side set. */ 1101 if (!fs_name_err) { 1102 ierr = DMSetLabelValue(*dm, fs_name, faces[0], fs_id[fs]);CHKERRQ(ierr); 1103 } 1104 ierr = DMPlexRestoreJoin(*dm, faceSize, faceVertices, &numFaces, &faces);CHKERRQ(ierr); 1105 } 1106 ierr = PetscFree2(fs_vertex_count_list,fs_vertex_list);CHKERRQ(ierr); 1107 } 1108 ierr = PetscFree(fs_id);CHKERRQ(ierr); 1109 } 1110 #else 1111 SETERRQ(comm, PETSC_ERR_SUP, "This method requires ExodusII support. Reconfigure using --download-exodusii"); 1112 #endif 1113 PetscFunctionReturn(0); 1114 } 1115