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