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 #include <petsc/private/viewerimpl.h> 10 #include <petsc/private/viewerexodusiiimpl.h> 11 #if defined(PETSC_HAVE_EXODUSII) 12 /* 13 PETSC_VIEWER_EXODUSII_ - Creates an ExodusII PetscViewer shared by all processors in a communicator. 14 15 Collective 16 17 Input Parameter: 18 . comm - the MPI communicator to share the ExodusII PetscViewer 19 20 Level: intermediate 21 22 Notes: 23 misses Fortran bindings 24 25 Notes: 26 Unlike almost all other PETSc routines, PETSC_VIEWER_EXODUSII_ does not return 27 an error code. The GLVIS PetscViewer is usually used in the form 28 $ XXXView(XXX object, PETSC_VIEWER_EXODUSII_(comm)); 29 30 .seealso: PetscViewerExodusIIOpen(), PetscViewerType, PetscViewerCreate(), PetscViewerDestroy() 31 */ 32 PetscViewer PETSC_VIEWER_EXODUSII_(MPI_Comm comm) 33 { 34 PetscViewer viewer; 35 PetscErrorCode ierr; 36 37 PetscFunctionBegin; 38 ierr = PetscViewerExodusIIOpen(comm, "mesh.exo", FILE_MODE_WRITE, &viewer); 39 if (ierr) {PetscError(PETSC_COMM_SELF,__LINE__,"PETSC_VIEWER_EXODUSII_",__FILE__,PETSC_ERR_PLIB,PETSC_ERROR_INITIAL," ");PetscFunctionReturn(0);} 40 ierr = PetscObjectRegisterDestroy((PetscObject) viewer); 41 if (ierr) {PetscError(PETSC_COMM_SELF,__LINE__,"PETSC_VIEWER_EXODUSII_",__FILE__,PETSC_ERR_PLIB,PETSC_ERROR_INITIAL," ");PetscFunctionReturn(0);} 42 PetscFunctionReturn(viewer); 43 } 44 45 static PetscErrorCode PetscViewerView_ExodusII(PetscViewer v, PetscViewer viewer) 46 { 47 PetscViewer_ExodusII *exo = (PetscViewer_ExodusII *) viewer->data; 48 PetscErrorCode ierr; 49 50 PetscFunctionBegin; 51 if (exo->filename) {ierr = PetscViewerASCIIPrintf(viewer, "Filename: %s\n", exo->filename);CHKERRQ(ierr);} 52 PetscFunctionReturn(0); 53 } 54 55 static PetscErrorCode PetscViewerSetFromOptions_ExodusII(PetscOptionItems *PetscOptionsObject, PetscViewer v) 56 { 57 PetscErrorCode ierr; 58 59 PetscFunctionBegin; 60 ierr = PetscOptionsHead(PetscOptionsObject, "ExodusII PetscViewer Options");CHKERRQ(ierr); 61 ierr = PetscOptionsTail();CHKERRQ(ierr); 62 PetscFunctionReturn(0); 63 } 64 65 static PetscErrorCode PetscViewerSetUp_ExodusII(PetscViewer viewer) 66 { 67 PetscFunctionBegin; 68 PetscFunctionReturn(0); 69 } 70 71 static PetscErrorCode PetscViewerDestroy_ExodusII(PetscViewer viewer) 72 { 73 PetscViewer_ExodusII *exo = (PetscViewer_ExodusII *) viewer->data; 74 PetscErrorCode ierr; 75 76 PetscFunctionBegin; 77 if (exo->exoid >= 0) {ierr = ex_close(exo->exoid);CHKERRQ(ierr);} 78 ierr = PetscFree(exo);CHKERRQ(ierr); 79 ierr = PetscObjectComposeFunction((PetscObject)viewer,"PetscViewerFileSetName_C",NULL);CHKERRQ(ierr); 80 ierr = PetscObjectComposeFunction((PetscObject)viewer,"PetscViewerFileGetName_C",NULL);CHKERRQ(ierr); 81 ierr = PetscObjectComposeFunction((PetscObject)viewer,"PetscViewerFileSetMode_C",NULL);CHKERRQ(ierr); 82 ierr = PetscObjectComposeFunction((PetscObject)viewer,"PetscViewerExodusIIGetId",NULL);CHKERRQ(ierr); 83 PetscFunctionReturn(0); 84 } 85 86 static PetscErrorCode PetscViewerFileSetName_ExodusII(PetscViewer viewer, const char name[]) 87 { 88 PetscViewer_ExodusII *exo = (PetscViewer_ExodusII *) viewer->data; 89 PetscMPIInt rank; 90 int CPU_word_size, IO_word_size, EXO_mode; 91 PetscErrorCode ierr; 92 93 94 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject) viewer), &rank);CHKERRQ(ierr); 95 CPU_word_size = sizeof(PetscReal); 96 IO_word_size = sizeof(PetscReal); 97 98 PetscFunctionBegin; 99 if (exo->exoid >= 0) {ex_close(exo->exoid); exo->exoid = -1;} 100 if (exo->filename) {ierr = PetscFree(exo->filename);CHKERRQ(ierr);} 101 ierr = PetscStrallocpy(name, &exo->filename);CHKERRQ(ierr); 102 /* Create or open the file collectively */ 103 switch (exo->btype) { 104 case FILE_MODE_READ: 105 EXO_mode = EX_CLOBBER; 106 break; 107 case FILE_MODE_APPEND: 108 EXO_mode = EX_CLOBBER; 109 break; 110 case FILE_MODE_WRITE: 111 EXO_mode = EX_CLOBBER; 112 break; 113 default: 114 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER, "Must call PetscViewerFileSetMode() before PetscViewerFileSetName()"); 115 } 116 #if defined(PETSC_USE_64BIT_INDICES) 117 EXO_mode += EX_ALL_INT64_API; 118 #endif 119 exo->exoid = ex_create(name, EXO_mode, &CPU_word_size, &IO_word_size); 120 if (exo->exoid < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_LIB, "ex_create failed for %s", name); 121 PetscFunctionReturn(0); 122 } 123 124 static PetscErrorCode PetscViewerFileGetName_ExodusII(PetscViewer viewer, const char **name) 125 { 126 PetscViewer_ExodusII *exo = (PetscViewer_ExodusII *) viewer->data; 127 128 PetscFunctionBegin; 129 *name = exo->filename; 130 PetscFunctionReturn(0); 131 } 132 133 static PetscErrorCode PetscViewerFileSetMode_ExodusII(PetscViewer viewer, PetscFileMode type) 134 { 135 PetscViewer_ExodusII *exo = (PetscViewer_ExodusII *) viewer->data; 136 137 PetscFunctionBegin; 138 exo->btype = type; 139 PetscFunctionReturn(0); 140 } 141 142 static PetscErrorCode PetscViewerFileGetMode_ExodusII(PetscViewer viewer, PetscFileMode *type) 143 { 144 PetscViewer_ExodusII *exo = (PetscViewer_ExodusII *) viewer->data; 145 146 PetscFunctionBegin; 147 *type = exo->btype; 148 PetscFunctionReturn(0); 149 } 150 151 static PetscErrorCode PetscViewerExodusIIGetId_ExodusII(PetscViewer viewer, int *exoid) 152 { 153 PetscViewer_ExodusII *exo = (PetscViewer_ExodusII *) viewer->data; 154 155 PetscFunctionBegin; 156 *exoid = exo->exoid; 157 PetscFunctionReturn(0); 158 } 159 160 /*@C 161 PetscViewerExodusIIOpen - Opens a file for ExodusII input/output. 162 163 Collective 164 165 Input Parameters: 166 + comm - MPI communicator 167 . name - name of file 168 - type - type of file 169 $ FILE_MODE_WRITE - create new file for binary output 170 $ FILE_MODE_READ - open existing file for binary input 171 $ FILE_MODE_APPEND - open existing file for binary output 172 173 Output Parameter: 174 . exo - PetscViewer for Exodus II input/output to use with the specified file 175 176 Level: beginner 177 178 Note: 179 This PetscViewer should be destroyed with PetscViewerDestroy(). 180 181 182 .seealso: PetscViewerPushFormat(), PetscViewerDestroy(), 183 DMLoad(), PetscFileMode, PetscViewer, PetscViewerSetType(), PetscViewerFileSetMode(), PetscViewerFileSetName() 184 @*/ 185 PetscErrorCode PetscViewerExodusIIOpen(MPI_Comm comm, const char name[], PetscFileMode type, PetscViewer *exo) 186 { 187 PetscErrorCode ierr; 188 189 PetscFunctionBegin; 190 ierr = PetscViewerCreate(comm, exo);CHKERRQ(ierr); 191 ierr = PetscViewerSetType(*exo, PETSCVIEWEREXODUSII);CHKERRQ(ierr); 192 ierr = PetscViewerFileSetMode(*exo, type);CHKERRQ(ierr); 193 ierr = PetscViewerFileSetName(*exo, name);CHKERRQ(ierr); 194 ierr = PetscViewerSetFromOptions(*exo);CHKERRQ(ierr); 195 PetscFunctionReturn(0); 196 } 197 198 /*MC 199 PETSCVIEWEREXODUSII - A viewer that writes to an Exodus II file 200 201 202 .seealso: PetscViewerExodusIIOpen(), PetscViewerCreate(), PETSCVIEWERBINARY, PETSCVIEWERHDF5, DMView(), 203 PetscViewerFileSetName(), PetscViewerFileSetMode(), PetscViewerFormat, PetscViewerType, PetscViewerSetType() 204 205 Level: beginner 206 M*/ 207 208 PETSC_EXTERN PetscErrorCode PetscViewerCreate_ExodusII(PetscViewer v) 209 { 210 PetscViewer_ExodusII *exo; 211 PetscErrorCode ierr; 212 213 PetscFunctionBegin; 214 ierr = PetscNewLog(v,&exo);CHKERRQ(ierr); 215 216 v->data = (void*) exo; 217 v->ops->destroy = PetscViewerDestroy_ExodusII; 218 v->ops->setfromoptions = PetscViewerSetFromOptions_ExodusII; 219 v->ops->setup = PetscViewerSetUp_ExodusII; 220 v->ops->view = PetscViewerView_ExodusII; 221 v->ops->flush = 0; 222 exo->btype = (PetscFileMode) -1; 223 exo->filename = 0; 224 exo->exoid = -1; 225 226 ierr = PetscObjectComposeFunction((PetscObject)v,"PetscViewerFileSetName_C",PetscViewerFileSetName_ExodusII);CHKERRQ(ierr); 227 ierr = PetscObjectComposeFunction((PetscObject)v,"PetscViewerFileGetName_C",PetscViewerFileGetName_ExodusII);CHKERRQ(ierr); 228 ierr = PetscObjectComposeFunction((PetscObject)v,"PetscViewerFileSetMode_C",PetscViewerFileSetMode_ExodusII);CHKERRQ(ierr); 229 ierr = PetscObjectComposeFunction((PetscObject)v,"PetscViewerFileGetMode_C",PetscViewerFileGetMode_ExodusII);CHKERRQ(ierr); 230 ierr = PetscObjectComposeFunction((PetscObject)v,"PetscViewerExodusIIGetId_C",PetscViewerExodusIIGetId_ExodusII);CHKERRQ(ierr); 231 PetscFunctionReturn(0); 232 } 233 234 /* 235 EXOGetVarIndex - Locate a result in an exodus file based on its name 236 237 Not collective 238 239 Input Parameters: 240 + exoid - the exodus id of a file (obtained from ex_open or ex_create for instance) 241 . obj_type - the type of entity for instance EX_NODAL, EX_ELEM_BLOCK 242 - name - the name of the result 243 244 Output Parameters: 245 . varIndex - the location in the exodus file of the result 246 247 Notes: 248 The exodus variable index is obtained by comparing name and the 249 names of zonal variables declared in the exodus file. For instance if name is "V" 250 the location in the exodus file will be the first match of "V", "V_X", "V_XX", "V_1", or "V_11" 251 amongst all variables of type obj_type. 252 253 Level: beginner 254 255 .seealso: EXOGetVarIndex(),DMPlexView_ExodusII_Internal(),VecViewPlex_ExodusII_Nodal_Internal(),VecLoadNodal_PlexEXO(),VecLoadZonal_PlexEXO() 256 */ 257 static PetscErrorCode EXOGetVarIndex_Private(int exoid, ex_entity_type obj_type, const char name[], int *varIndex) 258 { 259 int num_vars, i, j; 260 char ext_name[MAX_STR_LENGTH+1], var_name[MAX_STR_LENGTH+1]; 261 const int num_suffix = 5; 262 char *suffix[5]; 263 PetscBool flg; 264 PetscErrorCode ierr; 265 266 PetscFunctionBegin; 267 suffix[0] = (char *) ""; 268 suffix[1] = (char *) "_X"; 269 suffix[2] = (char *) "_XX"; 270 suffix[3] = (char *) "_1"; 271 suffix[4] = (char *) "_11"; 272 273 *varIndex = 0; 274 PetscStackCallStandard(ex_get_variable_param,(exoid, obj_type, &num_vars)); 275 for (i = 0; i < num_vars; ++i) { 276 PetscStackCallStandard(ex_get_variable_name,(exoid, obj_type, i+1, var_name)); 277 for (j = 0; j < num_suffix; ++j){ 278 ierr = PetscStrncpy(ext_name, name, MAX_STR_LENGTH);CHKERRQ(ierr); 279 ierr = PetscStrlcat(ext_name, suffix[j], MAX_STR_LENGTH);CHKERRQ(ierr); 280 ierr = PetscStrcasecmp(ext_name, var_name, &flg); 281 if (flg) { 282 *varIndex = i+1; 283 PetscFunctionReturn(0); 284 } 285 } 286 } 287 SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unable to locate variable %s in ExodusII file.", name); 288 PetscFunctionReturn(-1); 289 } 290 291 /* 292 DMPlexView_ExodusII_Internal - Write a DM to disk in exodus format 293 294 Collective on dm 295 296 Input Parameters: 297 + dm - The dm to be written 298 . exoid - the exodus id of a file (obtained from ex_open or ex_create for instance) 299 - degree - the degree of the interpolation space 300 301 Notes: 302 Not all DM can be written to disk this way. For instance, exodus assume that element blocks (mapped to "Cell sets" labels) 303 consists of sequentially numbered cells. If this is not the case, the exodus file will be corrupted. 304 305 If the dm has been distributed and exoid points to different files on each MPI rank, only the "local" part of the DM 306 (including "ghost" cells and vertices) will be written. If all exoid points to the same file, the resulting file will 307 probably be corrupted. 308 309 DMPlex only represents geometry while most post-processing software expect that a mesh also provides information 310 on the discretization space. This function assumes that the file represents Lagrange finite elements of order 1 or 2. 311 It should be extended to use PetscFE objects. 312 313 This function will only handle TRI, TET, QUAD and HEX cells. 314 Level: beginner 315 316 .seealso: EXOGetVarIndex_Private(), VecViewPlex_ExodusII_Nodal_Internal(), VecLoadNodal_PlexEXO(), VecViewZonal_PlexEXO(), VecLoadZonal_PlexEXO() 317 */ 318 PetscErrorCode DMPlexView_ExodusII_Internal(DM dm, int exoid, PetscInt degree) 319 { 320 enum ElemType {TRI, QUAD, TET, HEX}; 321 MPI_Comm comm; 322 /* Connectivity Variables */ 323 PetscInt cellsNotInConnectivity; 324 /* Cell Sets */ 325 DMLabel csLabel; 326 IS csIS; 327 const PetscInt *csIdx; 328 PetscInt num_cs, cs; 329 enum ElemType *type; 330 PetscBool hasLabel; 331 /* Coordinate Variables */ 332 DM cdm; 333 PetscSection section; 334 Vec coord; 335 PetscInt **nodes; 336 PetscInt depth, d, dim, skipCells = 0; 337 PetscInt pStart, pEnd, p, cStart, cEnd, numCells, vStart, vEnd, numVertices, eStart, eEnd, numEdges, fStart, fEnd, numFaces, numNodes; 338 PetscInt num_vs, num_fs; 339 PetscMPIInt rank, size; 340 const char *dmName; 341 PetscInt nodesTriP1[4] = {3,0,0,0}; 342 PetscInt nodesTriP2[4] = {3,3,0,0}; 343 PetscInt nodesQuadP1[4] = {4,0,0,0}; 344 PetscInt nodesQuadP2[4] = {4,4,0,1}; 345 PetscInt nodesTetP1[4] = {4,0,0,0}; 346 PetscInt nodesTetP2[4] = {4,6,0,0}; 347 PetscInt nodesHexP1[4] = {8,0,0,0}; 348 PetscInt nodesHexP2[4] = {8,12,6,1}; 349 PetscErrorCode ierr; 350 351 PetscFunctionBegin; 352 ierr = PetscObjectGetComm((PetscObject)dm, &comm);CHKERRQ(ierr); 353 ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr); 354 ierr = MPI_Comm_size(comm, &size);CHKERRQ(ierr); 355 /* --- Get DM info --- */ 356 ierr = PetscObjectGetName((PetscObject) dm, &dmName);CHKERRQ(ierr); 357 ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr); 358 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 359 ierr = DMPlexGetChart(dm, &pStart, &pEnd);CHKERRQ(ierr); 360 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 361 ierr = DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);CHKERRQ(ierr); 362 ierr = DMPlexGetDepthStratum(dm, 1, &eStart, &eEnd);CHKERRQ(ierr); 363 ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr); 364 numCells = cEnd - cStart; 365 numEdges = eEnd - eStart; 366 numVertices = vEnd - vStart; 367 if (depth == 3) {numFaces = fEnd - fStart;} 368 else {numFaces = 0;} 369 ierr = DMGetLabelSize(dm, "Cell Sets", &num_cs);CHKERRQ(ierr); 370 ierr = DMGetLabelSize(dm, "Vertex Sets", &num_vs);CHKERRQ(ierr); 371 ierr = DMGetLabelSize(dm, "Face Sets", &num_fs);CHKERRQ(ierr); 372 ierr = DMGetCoordinatesLocal(dm, &coord);CHKERRQ(ierr); 373 ierr = DMGetCoordinateDM(dm, &cdm);CHKERRQ(ierr); 374 if (num_cs > 0) { 375 ierr = DMGetLabel(dm, "Cell Sets", &csLabel);CHKERRQ(ierr); 376 ierr = DMLabelGetValueIS(csLabel, &csIS);CHKERRQ(ierr); 377 ierr = ISGetIndices(csIS, &csIdx);CHKERRQ(ierr); 378 } 379 ierr = PetscMalloc1(num_cs, &nodes);CHKERRQ(ierr); 380 /* Set element type for each block and compute total number of nodes */ 381 ierr = PetscMalloc1(num_cs, &type);CHKERRQ(ierr); 382 numNodes = numVertices; 383 if (degree == 2) {numNodes += numEdges;} 384 cellsNotInConnectivity = numCells; 385 for (cs = 0; cs < num_cs; ++cs) { 386 IS stratumIS; 387 const PetscInt *cells; 388 PetscScalar *xyz = NULL; 389 PetscInt csSize, closureSize; 390 391 ierr = DMLabelGetStratumIS(csLabel, csIdx[cs], &stratumIS);CHKERRQ(ierr); 392 ierr = ISGetIndices(stratumIS, &cells);CHKERRQ(ierr); 393 ierr = ISGetSize(stratumIS, &csSize);CHKERRQ(ierr); 394 ierr = DMPlexVecGetClosure(cdm, NULL, coord, cells[0], &closureSize, &xyz);CHKERRQ(ierr); 395 switch (dim) { 396 case 2: 397 if (closureSize == 3*dim) {type[cs] = TRI;} 398 else if (closureSize == 4*dim) {type[cs] = QUAD;} 399 else SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Number of vertices %D in dimension %D has no ExodusII type", closureSize/dim, dim); 400 break; 401 case 3: 402 if (closureSize == 4*dim) {type[cs] = TET;} 403 else if (closureSize == 8*dim) {type[cs] = HEX;} 404 else SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Number of vertices %D in dimension %D has no ExodusII type", closureSize/dim, dim); 405 break; 406 default: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Dimension %D not handled by ExodusII viewer", dim); 407 } 408 if ((degree == 2) && (type[cs] == QUAD)) {numNodes += csSize;} 409 if ((degree == 2) && (type[cs] == HEX)) {numNodes += csSize; numNodes += numFaces;} 410 ierr = DMPlexVecRestoreClosure(cdm, NULL, coord, cells[0], &closureSize, &xyz);CHKERRQ(ierr); 411 /* Set nodes and Element type */ 412 if (type[cs] == TRI) { 413 if (degree == 1) nodes[cs] = nodesTriP1; 414 else if (degree == 2) nodes[cs] = nodesTriP2; 415 } else if (type[cs] == QUAD) { 416 if (degree == 1) nodes[cs] = nodesQuadP1; 417 else if (degree == 2) nodes[cs] = nodesQuadP2; 418 } else if (type[cs] == TET) { 419 if (degree == 1) nodes[cs] = nodesTetP1; 420 else if (degree == 2) nodes[cs] = nodesTetP2; 421 } else if (type[cs] == HEX) { 422 if (degree == 1) nodes[cs] = nodesHexP1; 423 else if (degree == 2) nodes[cs] = nodesHexP2; 424 } 425 /* Compute the number of cells not in the connectivity table */ 426 cellsNotInConnectivity -= nodes[cs][3]*csSize; 427 428 ierr = ISRestoreIndices(stratumIS, &cells);CHKERRQ(ierr); 429 ierr = ISDestroy(&stratumIS);CHKERRQ(ierr); 430 } 431 if (num_cs > 0) {PetscStackCallStandard(ex_put_init,(exoid, dmName, dim, numNodes, numCells, num_cs, num_vs, num_fs));} 432 /* --- Connectivity --- */ 433 for (cs = 0; cs < num_cs; ++cs) { 434 IS stratumIS; 435 const PetscInt *cells; 436 PetscInt *connect, off = 0; 437 PetscInt edgesInClosure = 0, facesInClosure = 0, verticesInClosure = 0; 438 PetscInt csSize, c, connectSize, closureSize; 439 char *elem_type = NULL; 440 char elem_type_tri3[] = "TRI3", elem_type_quad4[] = "QUAD4"; 441 char elem_type_tri6[] = "TRI6", elem_type_quad9[] = "QUAD9"; 442 char elem_type_tet4[] = "TET4", elem_type_hex8[] = "HEX8"; 443 char elem_type_tet10[] = "TET10", elem_type_hex27[] = "HEX27"; 444 445 ierr = DMLabelGetStratumIS(csLabel, csIdx[cs], &stratumIS);CHKERRQ(ierr); 446 ierr = ISGetIndices(stratumIS, &cells);CHKERRQ(ierr); 447 ierr = ISGetSize(stratumIS, &csSize);CHKERRQ(ierr); 448 /* Set Element type */ 449 if (type[cs] == TRI) { 450 if (degree == 1) elem_type = elem_type_tri3; 451 else if (degree == 2) elem_type = elem_type_tri6; 452 } else if (type[cs] == QUAD) { 453 if (degree == 1) elem_type = elem_type_quad4; 454 else if (degree == 2) elem_type = elem_type_quad9; 455 } else if (type[cs] == TET) { 456 if (degree == 1) elem_type = elem_type_tet4; 457 else if (degree == 2) elem_type = elem_type_tet10; 458 } else if (type[cs] == HEX) { 459 if (degree == 1) elem_type = elem_type_hex8; 460 else if (degree == 2) elem_type = elem_type_hex27; 461 } 462 connectSize = nodes[cs][0] + nodes[cs][1] + nodes[cs][2] + nodes[cs][3]; 463 ierr = PetscMalloc1(PetscMax(27,connectSize)*csSize, &connect);CHKERRQ(ierr); 464 PetscStackCallStandard(ex_put_block,(exoid, EX_ELEM_BLOCK, csIdx[cs], elem_type, csSize, connectSize, 0, 0, 1)); 465 /* Find number of vertices, edges, and faces in the closure */ 466 verticesInClosure = nodes[cs][0]; 467 if (depth > 1) { 468 if (dim == 2) { 469 ierr = DMPlexGetConeSize(dm, cells[0], &edgesInClosure);CHKERRQ(ierr); 470 } else if (dim == 3) { 471 PetscInt *closure = NULL; 472 473 ierr = DMPlexGetConeSize(dm, cells[0], &facesInClosure);CHKERRQ(ierr); 474 ierr = DMPlexGetTransitiveClosure(dm, cells[0], PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 475 edgesInClosure = closureSize - facesInClosure - 1 - verticesInClosure; 476 ierr = DMPlexRestoreTransitiveClosure(dm, cells[0], PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 477 } 478 } 479 /* Get connectivity for each cell */ 480 for (c = 0; c < csSize; ++c) { 481 PetscInt *closure = NULL; 482 PetscInt temp, i; 483 484 ierr = DMPlexGetTransitiveClosure(dm, cells[c], PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 485 for (i = 0; i < connectSize; ++i) { 486 if (i < nodes[cs][0]) {/* Vertices */ 487 connect[i+off] = closure[(i+edgesInClosure+facesInClosure+1)*2] + 1; 488 connect[i+off] -= cellsNotInConnectivity; 489 } else if (i < nodes[cs][0]+nodes[cs][1]) { /* Edges */ 490 connect[i+off] = closure[(i-verticesInClosure+facesInClosure+1)*2] + 1; 491 if (nodes[cs][2] == 0) connect[i+off] -= numFaces; 492 connect[i+off] -= cellsNotInConnectivity; 493 } else if (i < nodes[cs][0]+nodes[cs][1]+nodes[cs][3]){ /* Cells */ 494 connect[i+off] = closure[0] + 1; 495 connect[i+off] -= skipCells; 496 } else if (i < nodes[cs][0]+nodes[cs][1]+nodes[cs][3]+nodes[cs][2]){ /* Faces */ 497 connect[i+off] = closure[(i-edgesInClosure-verticesInClosure)*2] + 1; 498 connect[i+off] -= cellsNotInConnectivity; 499 } else { 500 connect[i+off] = -1; 501 } 502 } 503 /* Tetrahedra are inverted */ 504 if (type[cs] == TET) { 505 temp = connect[0+off]; connect[0+off] = connect[1+off]; connect[1+off] = temp; 506 if (degree == 2) { 507 temp = connect[5+off]; connect[5+off] = connect[6+off]; connect[6+off] = temp; 508 temp = connect[7+off]; connect[7+off] = connect[8+off]; connect[8+off] = temp; 509 } 510 } 511 /* Hexahedra are inverted */ 512 if (type[cs] == HEX) { 513 temp = connect[1+off]; connect[1+off] = connect[3+off]; connect[3+off] = temp; 514 if (degree == 2) { 515 temp = connect[8+off]; connect[8+off] = connect[11+off]; connect[11+off] = temp; 516 temp = connect[9+off]; connect[9+off] = connect[10+off]; connect[10+off] = temp; 517 temp = connect[16+off]; connect[16+off] = connect[17+off]; connect[17+off] = temp; 518 temp = connect[18+off]; connect[18+off] = connect[19+off]; connect[19+off] = temp; 519 520 temp = connect[12+off]; connect[12+off] = connect[16+off]; connect[16+off] = temp; 521 temp = connect[13+off]; connect[13+off] = connect[17+off]; connect[17+off] = temp; 522 temp = connect[14+off]; connect[14+off] = connect[18+off]; connect[18+off] = temp; 523 temp = connect[15+off]; connect[15+off] = connect[19+off]; connect[19+off] = temp; 524 525 temp = connect[23+off]; connect[23+off] = connect[26+off]; connect[26+off] = temp; 526 temp = connect[24+off]; connect[24+off] = connect[25+off]; connect[25+off] = temp; 527 temp = connect[25+off]; connect[25+off] = connect[26+off]; connect[26+off] = temp; 528 } 529 } 530 off += connectSize; 531 ierr = DMPlexRestoreTransitiveClosure(dm, cells[c], PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 532 } 533 PetscStackCallStandard(ex_put_conn,(exoid, EX_ELEM_BLOCK, csIdx[cs], connect, 0, 0)); 534 skipCells += (nodes[cs][3] == 0)*csSize; 535 ierr = PetscFree(connect);CHKERRQ(ierr); 536 ierr = ISRestoreIndices(stratumIS, &cells);CHKERRQ(ierr); 537 ierr = ISDestroy(&stratumIS);CHKERRQ(ierr); 538 } 539 ierr = PetscFree(type);CHKERRQ(ierr); 540 /* --- Coordinates --- */ 541 ierr = PetscSectionCreate(comm, §ion);CHKERRQ(ierr); 542 ierr = PetscSectionSetChart(section, pStart, pEnd);CHKERRQ(ierr); 543 if (num_cs) { 544 for (d = 0; d < depth; ++d) { 545 ierr = DMPlexGetDepthStratum(dm, d, &pStart, &pEnd);CHKERRQ(ierr); 546 for (p = pStart; p < pEnd; ++p) { 547 ierr = PetscSectionSetDof(section, p, nodes[0][d] > 0);CHKERRQ(ierr); 548 } 549 } 550 } 551 for (cs = 0; cs < num_cs; ++cs) { 552 IS stratumIS; 553 const PetscInt *cells; 554 PetscInt csSize, c; 555 556 ierr = DMLabelGetStratumIS(csLabel, csIdx[cs], &stratumIS);CHKERRQ(ierr); 557 ierr = ISGetIndices(stratumIS, &cells);CHKERRQ(ierr); 558 ierr = ISGetSize(stratumIS, &csSize);CHKERRQ(ierr); 559 for (c = 0; c < csSize; ++c) { 560 ierr = PetscSectionSetDof(section, cells[c], nodes[cs][3] > 0);CHKERRQ(ierr); 561 } 562 ierr = ISRestoreIndices(stratumIS, &cells);CHKERRQ(ierr); 563 ierr = ISDestroy(&stratumIS);CHKERRQ(ierr); 564 } 565 if (num_cs > 0) { 566 ierr = ISRestoreIndices(csIS, &csIdx);CHKERRQ(ierr); 567 ierr = ISDestroy(&csIS);CHKERRQ(ierr); 568 } 569 ierr = PetscFree(nodes);CHKERRQ(ierr); 570 ierr = PetscSectionSetUp(section);CHKERRQ(ierr); 571 if (numNodes > 0) { 572 const char *coordNames[3] = {"x", "y", "z"}; 573 PetscScalar *coords, *closure; 574 PetscReal *cval; 575 PetscInt hasDof, n = 0; 576 577 /* There can't be more than 24 values in the closure of a point for the coord section */ 578 ierr = PetscCalloc3(numNodes*3, &coords, dim, &cval, 24, &closure);CHKERRQ(ierr); 579 ierr = DMGetCoordinatesLocal(dm, &coord);CHKERRQ(ierr); 580 ierr = DMPlexGetChart(dm, &pStart, &pEnd);CHKERRQ(ierr); 581 for (p = pStart; p < pEnd; ++p) { 582 ierr = PetscSectionGetDof(section, p, &hasDof); 583 if (hasDof) { 584 PetscInt closureSize = 24, j; 585 586 ierr = DMPlexVecGetClosure(cdm, NULL, coord, p, &closureSize, &closure);CHKERRQ(ierr); 587 for (d = 0; d < dim; ++d) { 588 cval[d] = 0.0; 589 for (j = 0; j < closureSize/dim; j++) cval[d] += closure[j*dim+d]; 590 coords[d*numNodes+n] = cval[d] * dim / closureSize; 591 } 592 ++n; 593 } 594 } 595 PetscStackCallStandard(ex_put_coord,(exoid, &coords[0*numNodes], &coords[1*numNodes], &coords[2*numNodes])); 596 ierr = PetscFree3(coords, cval, closure);CHKERRQ(ierr); 597 PetscStackCallStandard(ex_put_coord_names,(exoid, (char **) coordNames)); 598 } 599 ierr = PetscSectionDestroy(§ion);CHKERRQ(ierr); 600 /* --- Node Sets/Vertex Sets --- */ 601 DMHasLabel(dm, "Vertex Sets", &hasLabel); 602 if (hasLabel) { 603 PetscInt i, vs, vsSize; 604 const PetscInt *vsIdx, *vertices; 605 PetscInt *nodeList; 606 IS vsIS, stratumIS; 607 DMLabel vsLabel; 608 ierr = DMGetLabel(dm, "Vertex Sets", &vsLabel);CHKERRQ(ierr); 609 ierr = DMLabelGetValueIS(vsLabel, &vsIS);CHKERRQ(ierr); 610 ierr = ISGetIndices(vsIS, &vsIdx);CHKERRQ(ierr); 611 for (vs=0; vs<num_vs; ++vs) { 612 ierr = DMLabelGetStratumIS(vsLabel, vsIdx[vs], &stratumIS);CHKERRQ(ierr); 613 ierr = ISGetIndices(stratumIS, &vertices);CHKERRQ(ierr); 614 ierr = ISGetSize(stratumIS, &vsSize);CHKERRQ(ierr); 615 ierr = PetscMalloc1(vsSize, &nodeList); 616 for (i=0; i<vsSize; ++i) { 617 nodeList[i] = vertices[i] - skipCells + 1; 618 } 619 PetscStackCallStandard(ex_put_set_param,(exoid, EX_NODE_SET, vsIdx[vs], vsSize, 0)); 620 PetscStackCallStandard(ex_put_set,(exoid, EX_NODE_SET, vsIdx[vs], nodeList, NULL)); 621 ierr = ISRestoreIndices(stratumIS, &vertices);CHKERRQ(ierr); 622 ierr = ISDestroy(&stratumIS);CHKERRQ(ierr); 623 ierr = PetscFree(nodeList);CHKERRQ(ierr); 624 } 625 ierr = ISRestoreIndices(vsIS, &vsIdx);CHKERRQ(ierr); 626 ierr = ISDestroy(&vsIS);CHKERRQ(ierr); 627 } 628 /* --- Side Sets/Face Sets --- */ 629 ierr = DMHasLabel(dm, "Face Sets", &hasLabel);CHKERRQ(ierr); 630 if (hasLabel) { 631 PetscInt i, j, fs, fsSize; 632 const PetscInt *fsIdx, *faces; 633 IS fsIS, stratumIS; 634 DMLabel fsLabel; 635 PetscInt numPoints, *points; 636 PetscInt elem_list_size = 0; 637 PetscInt *elem_list, *elem_ind, *side_list; 638 639 ierr = DMGetLabel(dm, "Face Sets", &fsLabel);CHKERRQ(ierr); 640 /* Compute size of Node List and Element List */ 641 ierr = DMLabelGetValueIS(fsLabel, &fsIS);CHKERRQ(ierr); 642 ierr = ISGetIndices(fsIS, &fsIdx);CHKERRQ(ierr); 643 for (fs=0; fs<num_fs; ++fs) { 644 ierr = DMLabelGetStratumIS(fsLabel, fsIdx[fs], &stratumIS);CHKERRQ(ierr); 645 ierr = ISGetSize(stratumIS, &fsSize);CHKERRQ(ierr); 646 elem_list_size += fsSize; 647 ierr = ISDestroy(&stratumIS);CHKERRQ(ierr); 648 } 649 ierr = PetscMalloc3(num_fs, &elem_ind, 650 elem_list_size, &elem_list, 651 elem_list_size, &side_list);CHKERRQ(ierr); 652 elem_ind[0] = 0; 653 for (fs=0; fs<num_fs; ++fs) { 654 ierr = DMLabelGetStratumIS(fsLabel, fsIdx[fs], &stratumIS);CHKERRQ(ierr); 655 ierr = ISGetIndices(stratumIS, &faces);CHKERRQ(ierr); 656 ierr = ISGetSize(stratumIS, &fsSize);CHKERRQ(ierr); 657 /* Set Parameters */ 658 PetscStackCallStandard(ex_put_set_param,(exoid, EX_SIDE_SET, fsIdx[fs], fsSize, 0)); 659 /* Indices */ 660 if (fs<num_fs-1) { 661 elem_ind[fs+1] = elem_ind[fs] + fsSize; 662 } 663 664 for (i=0; i<fsSize; ++i) { 665 /* Element List */ 666 points = NULL; 667 ierr = DMPlexGetTransitiveClosure(dm, faces[i], PETSC_FALSE, &numPoints, &points);CHKERRQ(ierr); 668 elem_list[elem_ind[fs] + i] = points[2] +1; 669 ierr = DMPlexRestoreTransitiveClosure(dm, faces[i], PETSC_FALSE, &numPoints, &points);CHKERRQ(ierr); 670 671 /* Side List */ 672 points = NULL; 673 ierr = DMPlexGetTransitiveClosure(dm, elem_list[elem_ind[fs] + i]-1, PETSC_TRUE, &numPoints, &points);CHKERRQ(ierr); 674 for (j=1; j<numPoints; ++j) { 675 if (points[j*2]==faces[i]) {break;} 676 } 677 /* Convert HEX sides */ 678 if (numPoints == 27) { 679 if (j == 1) {j = 5;} 680 else if (j == 2) {j = 6;} 681 else if (j == 3) {j = 1;} 682 else if (j == 4) {j = 3;} 683 else if (j == 5) {j = 2;} 684 else if (j == 6) {j = 4;} 685 } 686 /* Convert TET sides */ 687 if (numPoints == 15) { 688 --j; 689 if (j == 0) {j = 4;} 690 } 691 side_list[elem_ind[fs] + i] = j; 692 ierr = DMPlexRestoreTransitiveClosure(dm, elem_list[elem_ind[fs] + i]-1, PETSC_TRUE, &numPoints, &points);CHKERRQ(ierr); 693 694 } 695 ierr = ISRestoreIndices(stratumIS, &faces);CHKERRQ(ierr); 696 ierr = ISDestroy(&stratumIS);CHKERRQ(ierr); 697 } 698 ierr = ISRestoreIndices(fsIS, &fsIdx);CHKERRQ(ierr); 699 ierr = ISDestroy(&fsIS);CHKERRQ(ierr); 700 701 /* Put side sets */ 702 for (fs=0; fs<num_fs; ++fs) { 703 PetscStackCallStandard(ex_put_set,(exoid, EX_SIDE_SET, fsIdx[fs], &elem_list[elem_ind[fs]], &side_list[elem_ind[fs]])); 704 } 705 ierr = PetscFree3(elem_ind,elem_list,side_list);CHKERRQ(ierr); 706 } 707 PetscFunctionReturn(0); 708 } 709 710 /* 711 VecViewPlex_ExodusII_Nodal_Internal - Write a Vec corresponding to a nodal field to an exodus file 712 713 Collective on v 714 715 Input Parameters: 716 + v - The vector to be written 717 . exoid - the exodus id of a file (obtained from ex_open or ex_create for instance) 718 - step - the time step to write at (exodus steps are numbered starting from 1) 719 720 Notes: 721 The exodus result nodal variable index is obtained by comparing the Vec name and the 722 names of nodal variables declared in the exodus file. For instance for a Vec named "V" 723 the location in the exodus file will be the first match of "V", "V_X", "V_XX", "V_1", or "V_11" 724 amongst all nodal variables. 725 726 Level: beginner 727 728 .seealso: EXOGetVarIndex_Private(),DMPlexView_ExodusII_Internal(),VecLoadNodal_PlexEXO(),VecViewZonal_PlexEXO(),VecLoadZonal_PlexEXO() 729 @*/ 730 PetscErrorCode VecViewPlex_ExodusII_Nodal_Internal(Vec v, int exoid, int step) 731 { 732 MPI_Comm comm; 733 PetscMPIInt size; 734 DM dm; 735 Vec vNatural, vComp; 736 const PetscScalar *varray; 737 const char *vecname; 738 PetscInt xs, xe, bs; 739 PetscBool useNatural; 740 int offset; 741 PetscErrorCode ierr; 742 743 PetscFunctionBegin; 744 ierr = PetscObjectGetComm((PetscObject) v, &comm);CHKERRQ(ierr); 745 ierr = MPI_Comm_size(comm, &size);CHKERRQ(ierr); 746 ierr = VecGetDM(v, &dm);CHKERRQ(ierr); 747 ierr = DMGetUseNatural(dm, &useNatural);CHKERRQ(ierr); 748 useNatural = useNatural && size > 1 ? PETSC_TRUE : PETSC_FALSE; 749 if (useNatural) { 750 ierr = DMGetGlobalVector(dm, &vNatural);CHKERRQ(ierr); 751 ierr = DMPlexGlobalToNaturalBegin(dm, v, vNatural);CHKERRQ(ierr); 752 ierr = DMPlexGlobalToNaturalEnd(dm, v, vNatural);CHKERRQ(ierr); 753 } else { 754 vNatural = v; 755 } 756 /* Get the location of the variable in the exodus file */ 757 ierr = PetscObjectGetName((PetscObject) v, &vecname);CHKERRQ(ierr); 758 ierr = EXOGetVarIndex_Private(exoid, EX_NODAL, vecname, &offset);CHKERRQ(ierr); 759 if (!offset) SETERRQ1(comm, PETSC_ERR_FILE_UNEXPECTED, "Cannot locate nodal variable %s in exodus file\n", vecname); 760 /* Write local chunk of the result in the exodus file 761 exodus stores each component of a vector-valued field as a separate variable. 762 We assume that they are stored sequentially */ 763 ierr = VecGetOwnershipRange(vNatural, &xs, &xe);CHKERRQ(ierr); 764 ierr = VecGetBlockSize(vNatural, &bs);CHKERRQ(ierr); 765 if (bs == 1) { 766 ierr = VecGetArrayRead(vNatural, &varray);CHKERRQ(ierr); 767 PetscStackCallStandard(ex_put_partial_var,(exoid, step, EX_NODAL, offset, 1, xs+1, xe-xs, varray)); 768 ierr = VecRestoreArrayRead(vNatural, &varray);CHKERRQ(ierr); 769 } else { 770 IS compIS; 771 PetscInt c; 772 773 ierr = ISCreateStride(comm, (xe-xs)/bs, xs, bs, &compIS);CHKERRQ(ierr); 774 for (c = 0; c < bs; ++c) { 775 ierr = ISStrideSetStride(compIS, (xe-xs)/bs, xs+c, bs);CHKERRQ(ierr); 776 ierr = VecGetSubVector(vNatural, compIS, &vComp);CHKERRQ(ierr); 777 ierr = VecGetArrayRead(vComp, &varray);CHKERRQ(ierr); 778 PetscStackCallStandard(ex_put_partial_var,(exoid, step, EX_NODAL, offset+c, 1, xs/bs+1, (xe-xs)/bs, varray)); 779 ierr = VecRestoreArrayRead(vComp, &varray);CHKERRQ(ierr); 780 ierr = VecRestoreSubVector(vNatural, compIS, &vComp);CHKERRQ(ierr); 781 } 782 ierr = ISDestroy(&compIS);CHKERRQ(ierr); 783 } 784 if (useNatural) {ierr = DMRestoreGlobalVector(dm, &vNatural);CHKERRQ(ierr);} 785 PetscFunctionReturn(0); 786 } 787 788 /* 789 VecLoadPlex_ExodusII_Nodal_Internal - Read a Vec corresponding to a nodal field from an exodus file 790 791 Collective on v 792 793 Input Parameters: 794 + v - The vector to be written 795 . exoid - the exodus id of a file (obtained from ex_open or ex_create for instance) 796 - step - the time step to read at (exodus steps are numbered starting from 1) 797 798 Notes: 799 The exodus result nodal variable index is obtained by comparing the Vec name and the 800 names of nodal variables declared in the exodus file. For instance for a Vec named "V" 801 the location in the exodus file will be the first match of "V", "V_X", "V_XX", "V_1", or "V_11" 802 amongst all nodal variables. 803 804 Level: beginner 805 806 .seealso: EXOGetVarIndex_Private(), DMPlexView_ExodusII_Internal(), VecViewPlex_ExodusII_Nodal_Internal(), VecViewZonal_PlexEXO(), VecLoadZonal_PlexEXO() 807 */ 808 PetscErrorCode VecLoadPlex_ExodusII_Nodal_Internal(Vec v, int exoid, int step) 809 { 810 MPI_Comm comm; 811 PetscMPIInt size; 812 DM dm; 813 Vec vNatural, vComp; 814 PetscScalar *varray; 815 const char *vecname; 816 PetscInt xs, xe, bs; 817 PetscBool useNatural; 818 int offset; 819 PetscErrorCode ierr; 820 821 PetscFunctionBegin; 822 ierr = PetscObjectGetComm((PetscObject) v, &comm);CHKERRQ(ierr); 823 ierr = MPI_Comm_size(comm, &size);CHKERRQ(ierr); 824 ierr = VecGetDM(v,&dm);CHKERRQ(ierr); 825 ierr = DMGetUseNatural(dm, &useNatural);CHKERRQ(ierr); 826 useNatural = useNatural && size > 1 ? PETSC_TRUE : PETSC_FALSE; 827 if (useNatural) {ierr = DMGetGlobalVector(dm,&vNatural);CHKERRQ(ierr);} 828 else {vNatural = v;} 829 /* Get the location of the variable in the exodus file */ 830 ierr = PetscObjectGetName((PetscObject) v, &vecname);CHKERRQ(ierr); 831 ierr = EXOGetVarIndex_Private(exoid, EX_NODAL, vecname, &offset);CHKERRQ(ierr); 832 if (!offset) SETERRQ1(comm, PETSC_ERR_FILE_UNEXPECTED, "Cannot locate nodal variable %s in exodus file", vecname); 833 /* Read local chunk from the file */ 834 ierr = VecGetOwnershipRange(vNatural, &xs, &xe);CHKERRQ(ierr); 835 ierr = VecGetBlockSize(vNatural, &bs);CHKERRQ(ierr); 836 if (bs == 1) { 837 ierr = VecGetArray(vNatural, &varray);CHKERRQ(ierr); 838 PetscStackCallStandard(ex_get_partial_var,(exoid, step, EX_NODAL, offset, 1, xs+1, xe-xs, varray)); 839 ierr = VecRestoreArray(vNatural, &varray);CHKERRQ(ierr); 840 } else { 841 IS compIS; 842 PetscInt c; 843 844 ierr = ISCreateStride(comm, (xe-xs)/bs, xs, bs, &compIS);CHKERRQ(ierr); 845 for (c = 0; c < bs; ++c) { 846 ierr = ISStrideSetStride(compIS, (xe-xs)/bs, xs+c, bs);CHKERRQ(ierr); 847 ierr = VecGetSubVector(vNatural, compIS, &vComp);CHKERRQ(ierr); 848 ierr = VecGetArray(vComp, &varray);CHKERRQ(ierr); 849 PetscStackCallStandard(ex_get_partial_var,(exoid, step, EX_NODAL, offset+c, 1, xs/bs+1, (xe-xs)/bs, varray)); 850 ierr = VecRestoreArray(vComp, &varray);CHKERRQ(ierr); 851 ierr = VecRestoreSubVector(vNatural, compIS, &vComp);CHKERRQ(ierr); 852 } 853 ierr = ISDestroy(&compIS);CHKERRQ(ierr); 854 } 855 if (useNatural) { 856 ierr = DMPlexNaturalToGlobalBegin(dm, vNatural, v);CHKERRQ(ierr); 857 ierr = DMPlexNaturalToGlobalEnd(dm, vNatural, v);CHKERRQ(ierr); 858 ierr = DMRestoreGlobalVector(dm, &vNatural);CHKERRQ(ierr); 859 } 860 PetscFunctionReturn(0); 861 } 862 863 /* 864 VecViewPlex_ExodusII_Zonal_Internal - Write a Vec corresponding to a zonal (cell based) field to an exodus file 865 866 Collective on v 867 868 Input Parameters: 869 + v - The vector to be written 870 . exoid - the exodus id of a file (obtained from ex_open or ex_create for instance) 871 - step - the time step to write at (exodus steps are numbered starting from 1) 872 873 Notes: 874 The exodus result zonal variable index is obtained by comparing the Vec name and the 875 names of zonal variables declared in the exodus file. For instance for a Vec named "V" 876 the location in the exodus file will be the first match of "V", "V_X", "V_XX", "V_1", or "V_11" 877 amongst all zonal variables. 878 879 Level: beginner 880 881 .seealso: EXOGetVarIndex_Private(),DMPlexView_ExodusII_Internal(),VecViewPlex_ExodusII_Nodal_Internal(),VecLoadPlex_ExodusII_Nodal_Internal(),VecLoadPlex_ExodusII_Zonal_Internal() 882 */ 883 PetscErrorCode VecViewPlex_ExodusII_Zonal_Internal(Vec v, int exoid, int step) 884 { 885 MPI_Comm comm; 886 PetscMPIInt size; 887 DM dm; 888 Vec vNatural, vComp; 889 const PetscScalar *varray; 890 const char *vecname; 891 PetscInt xs, xe, bs; 892 PetscBool useNatural; 893 int offset; 894 IS compIS; 895 PetscInt *csSize, *csID; 896 PetscInt numCS, set, csxs = 0; 897 PetscErrorCode ierr; 898 899 PetscFunctionBegin; 900 ierr = PetscObjectGetComm((PetscObject)v, &comm);CHKERRQ(ierr); 901 ierr = MPI_Comm_size(comm, &size);CHKERRQ(ierr); 902 ierr = VecGetDM(v, &dm);CHKERRQ(ierr); 903 ierr = DMGetUseNatural(dm, &useNatural);CHKERRQ(ierr); 904 useNatural = useNatural && size > 1 ? PETSC_TRUE : PETSC_FALSE; 905 if (useNatural) { 906 ierr = DMGetGlobalVector(dm, &vNatural);CHKERRQ(ierr); 907 ierr = DMPlexGlobalToNaturalBegin(dm, v, vNatural);CHKERRQ(ierr); 908 ierr = DMPlexGlobalToNaturalEnd(dm, v, vNatural);CHKERRQ(ierr); 909 } else { 910 vNatural = v; 911 } 912 /* Get the location of the variable in the exodus file */ 913 ierr = PetscObjectGetName((PetscObject) v, &vecname);CHKERRQ(ierr); 914 ierr = EXOGetVarIndex_Private(exoid, EX_ELEM_BLOCK, vecname, &offset);CHKERRQ(ierr); 915 if (!offset) SETERRQ1(comm, PETSC_ERR_FILE_UNEXPECTED, "Cannot locate zonal variable %s in exodus file\n", vecname); 916 /* Write local chunk of the result in the exodus file 917 exodus stores each component of a vector-valued field as a separate variable. 918 We assume that they are stored sequentially 919 Zonal variables are accessed one element block at a time, so we loop through the cell sets, 920 but once the vector has been reordered to natural size, we cannot use the label informations 921 to figure out what to save where. */ 922 numCS = ex_inquire_int(exoid, EX_INQ_ELEM_BLK); 923 ierr = PetscMalloc2(numCS, &csID, numCS, &csSize);CHKERRQ(ierr); 924 PetscStackCallStandard(ex_get_ids,(exoid, EX_ELEM_BLOCK, csID)); 925 for (set = 0; set < numCS; ++set) { 926 ex_block block; 927 928 block.id = csID[set]; 929 block.type = EX_ELEM_BLOCK; 930 PetscStackCallStandard(ex_get_block_param,(exoid, &block)); 931 csSize[set] = block.num_entry; 932 } 933 ierr = VecGetOwnershipRange(vNatural, &xs, &xe);CHKERRQ(ierr); 934 ierr = VecGetBlockSize(vNatural, &bs);CHKERRQ(ierr); 935 if (bs > 1) {ierr = ISCreateStride(comm, (xe-xs)/bs, xs, bs, &compIS);CHKERRQ(ierr);} 936 for (set = 0; set < numCS; set++) { 937 PetscInt csLocalSize, c; 938 939 /* range of indices for set setID[set]: csxs:csxs + csSize[set]-1 940 local slice of zonal values: xs/bs,xm/bs-1 941 intersection: max(xs/bs,csxs),min(xm/bs-1,csxs + csSize[set]-1) */ 942 csLocalSize = PetscMax(0, PetscMin(xe/bs, csxs+csSize[set]) - PetscMax(xs/bs, csxs)); 943 if (bs == 1) { 944 ierr = VecGetArrayRead(vNatural, &varray);CHKERRQ(ierr); 945 PetscStackCallStandard(ex_put_partial_var,(exoid, step, EX_ELEM_BLOCK, offset, csID[set], PetscMax(xs-csxs, 0)+1, csLocalSize, &varray[PetscMax(0, csxs-xs)])); 946 ierr = VecRestoreArrayRead(vNatural, &varray);CHKERRQ(ierr); 947 } else { 948 for (c = 0; c < bs; ++c) { 949 ierr = ISStrideSetStride(compIS, (xe-xs)/bs, xs+c, bs);CHKERRQ(ierr); 950 ierr = VecGetSubVector(vNatural, compIS, &vComp);CHKERRQ(ierr); 951 ierr = VecGetArrayRead(vComp, &varray);CHKERRQ(ierr); 952 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)])); 953 ierr = VecRestoreArrayRead(vComp, &varray);CHKERRQ(ierr); 954 ierr = VecRestoreSubVector(vNatural, compIS, &vComp);CHKERRQ(ierr); 955 } 956 } 957 csxs += csSize[set]; 958 } 959 ierr = PetscFree2(csID, csSize);CHKERRQ(ierr); 960 if (bs > 1) {ierr = ISDestroy(&compIS);CHKERRQ(ierr);} 961 if (useNatural) {ierr = DMRestoreGlobalVector(dm,&vNatural);CHKERRQ(ierr);} 962 PetscFunctionReturn(0); 963 } 964 965 /* 966 VecLoadPlex_ExodusII_Zonal_Internal - Read a Vec corresponding to a zonal (cell based) field from an exodus file 967 968 Collective on v 969 970 Input Parameters: 971 + v - The vector to be written 972 . exoid - the exodus id of a file (obtained from ex_open or ex_create for instance) 973 - step - the time step to read at (exodus steps are numbered starting from 1) 974 975 Notes: 976 The exodus result zonal variable index is obtained by comparing the Vec name and the 977 names of zonal variables declared in the exodus file. For instance for a Vec named "V" 978 the location in the exodus file will be the first match of "V", "V_X", "V_XX", "V_1", or "V_11" 979 amongst all zonal variables. 980 981 Level: beginner 982 983 .seealso: EXOGetVarIndex_Private(), DMPlexView_ExodusII_Internal(), VecViewPlex_ExodusII_Nodal_Internal(), VecLoadPlex_ExodusII_Nodal_Internal(), VecLoadPlex_ExodusII_Zonal_Internal() 984 */ 985 PetscErrorCode VecLoadPlex_ExodusII_Zonal_Internal(Vec v, int exoid, int step) 986 { 987 MPI_Comm comm; 988 PetscMPIInt size; 989 DM dm; 990 Vec vNatural, vComp; 991 PetscScalar *varray; 992 const char *vecname; 993 PetscInt xs, xe, bs; 994 PetscBool useNatural; 995 int offset; 996 IS compIS; 997 PetscInt *csSize, *csID; 998 PetscInt numCS, set, csxs = 0; 999 PetscErrorCode ierr; 1000 1001 PetscFunctionBegin; 1002 ierr = PetscObjectGetComm((PetscObject)v,&comm);CHKERRQ(ierr); 1003 ierr = MPI_Comm_size(comm, &size);CHKERRQ(ierr); 1004 ierr = VecGetDM(v, &dm);CHKERRQ(ierr); 1005 ierr = DMGetUseNatural(dm, &useNatural);CHKERRQ(ierr); 1006 useNatural = useNatural && size > 1 ? PETSC_TRUE : PETSC_FALSE; 1007 if (useNatural) {ierr = DMGetGlobalVector(dm,&vNatural);CHKERRQ(ierr);} 1008 else {vNatural = v;} 1009 /* Get the location of the variable in the exodus file */ 1010 ierr = PetscObjectGetName((PetscObject) v, &vecname);CHKERRQ(ierr); 1011 ierr = EXOGetVarIndex_Private(exoid, EX_ELEM_BLOCK, vecname, &offset);CHKERRQ(ierr); 1012 if (!offset) SETERRQ1(comm, PETSC_ERR_FILE_UNEXPECTED, "Cannot locate zonal variable %s in exodus file\n", vecname); 1013 /* Read local chunk of the result in the exodus file 1014 exodus stores each component of a vector-valued field as a separate variable. 1015 We assume that they are stored sequentially 1016 Zonal variables are accessed one element block at a time, so we loop through the cell sets, 1017 but once the vector has been reordered to natural size, we cannot use the label informations 1018 to figure out what to save where. */ 1019 numCS = ex_inquire_int(exoid, EX_INQ_ELEM_BLK); 1020 ierr = PetscMalloc2(numCS, &csID, numCS, &csSize);CHKERRQ(ierr); 1021 PetscStackCallStandard(ex_get_ids,(exoid, EX_ELEM_BLOCK, csID)); 1022 for (set = 0; set < numCS; ++set) { 1023 ex_block block; 1024 1025 block.id = csID[set]; 1026 block.type = EX_ELEM_BLOCK; 1027 PetscStackCallStandard(ex_get_block_param,(exoid, &block)); 1028 csSize[set] = block.num_entry; 1029 } 1030 ierr = VecGetOwnershipRange(vNatural, &xs, &xe);CHKERRQ(ierr); 1031 ierr = VecGetBlockSize(vNatural, &bs);CHKERRQ(ierr); 1032 if (bs > 1) {ierr = ISCreateStride(comm, (xe-xs)/bs, xs, bs, &compIS);CHKERRQ(ierr);} 1033 for (set = 0; set < numCS; ++set) { 1034 PetscInt csLocalSize, c; 1035 1036 /* range of indices for set setID[set]: csxs:csxs + csSize[set]-1 1037 local slice of zonal values: xs/bs,xm/bs-1 1038 intersection: max(xs/bs,csxs),min(xm/bs-1,csxs + csSize[set]-1) */ 1039 csLocalSize = PetscMax(0, PetscMin(xe/bs, csxs+csSize[set]) - PetscMax(xs/bs, csxs)); 1040 if (bs == 1) { 1041 ierr = VecGetArray(vNatural, &varray);CHKERRQ(ierr); 1042 PetscStackCallStandard(ex_get_partial_var,(exoid, step, EX_ELEM_BLOCK, offset, csID[set], PetscMax(xs-csxs, 0)+1, csLocalSize, &varray[PetscMax(0, csxs-xs)])); 1043 ierr = VecRestoreArray(vNatural, &varray);CHKERRQ(ierr); 1044 } else { 1045 for (c = 0; c < bs; ++c) { 1046 ierr = ISStrideSetStride(compIS, (xe-xs)/bs, xs+c, bs);CHKERRQ(ierr); 1047 ierr = VecGetSubVector(vNatural, compIS, &vComp);CHKERRQ(ierr); 1048 ierr = VecGetArray(vComp, &varray);CHKERRQ(ierr); 1049 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)])); 1050 ierr = VecRestoreArray(vComp, &varray);CHKERRQ(ierr); 1051 ierr = VecRestoreSubVector(vNatural, compIS, &vComp);CHKERRQ(ierr); 1052 } 1053 } 1054 csxs += csSize[set]; 1055 } 1056 ierr = PetscFree2(csID, csSize);CHKERRQ(ierr); 1057 if (bs > 1) {ierr = ISDestroy(&compIS);CHKERRQ(ierr);} 1058 if (useNatural) { 1059 ierr = DMPlexNaturalToGlobalBegin(dm, vNatural, v);CHKERRQ(ierr); 1060 ierr = DMPlexNaturalToGlobalEnd(dm, vNatural, v);CHKERRQ(ierr); 1061 ierr = DMRestoreGlobalVector(dm, &vNatural);CHKERRQ(ierr); 1062 } 1063 PetscFunctionReturn(0); 1064 } 1065 #endif 1066 1067 /*@ 1068 PetscViewerExodusIIGetId - Get the file id of the ExodusII file 1069 1070 Logically Collective on PetscViewer 1071 1072 Input Parameter: 1073 . viewer - the PetscViewer 1074 1075 Output Parameter: 1076 - exoid - The ExodusII file id 1077 1078 Level: intermediate 1079 1080 .seealso: PetscViewerFileSetMode(), PetscViewerCreate(), PetscViewerSetType(), PetscViewerBinaryOpen() 1081 @*/ 1082 PetscErrorCode PetscViewerExodusIIGetId(PetscViewer viewer, int *exoid) 1083 { 1084 PetscErrorCode ierr; 1085 1086 PetscFunctionBegin; 1087 PetscValidHeaderSpecific(viewer, PETSC_VIEWER_CLASSID, 1); 1088 ierr = PetscTryMethod(viewer, "PetscViewerExodusIIGetId_C",(PetscViewer,int*),(viewer,exoid));CHKERRQ(ierr); 1089 PetscFunctionReturn(0); 1090 } 1091 1092 /*@C 1093 DMPlexCreateExodusFromFile - Create a DMPlex mesh from an ExodusII file. 1094 1095 Collective 1096 1097 Input Parameters: 1098 + comm - The MPI communicator 1099 . filename - The name of the ExodusII file 1100 - interpolate - Create faces and edges in the mesh 1101 1102 Output Parameter: 1103 . dm - The DM object representing the mesh 1104 1105 Level: beginner 1106 1107 .seealso: DMPLEX, DMCreate(), DMPlexCreateExodus() 1108 @*/ 1109 PetscErrorCode DMPlexCreateExodusFromFile(MPI_Comm comm, const char filename[], PetscBool interpolate, DM *dm) 1110 { 1111 PetscMPIInt rank; 1112 PetscErrorCode ierr; 1113 #if defined(PETSC_HAVE_EXODUSII) 1114 int CPU_word_size = sizeof(PetscReal), IO_word_size = 0, exoid = -1; 1115 float version; 1116 #endif 1117 1118 PetscFunctionBegin; 1119 PetscValidCharPointer(filename, 2); 1120 ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr); 1121 #if defined(PETSC_HAVE_EXODUSII) 1122 if (!rank) { 1123 exoid = ex_open(filename, EX_READ, &CPU_word_size, &IO_word_size, &version); 1124 if (exoid <= 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_LIB, "ex_open(\"%s\",...) did not return a valid file ID", filename); 1125 } 1126 ierr = DMPlexCreateExodus(comm, exoid, interpolate, dm);CHKERRQ(ierr); 1127 if (!rank) {PetscStackCallStandard(ex_close,(exoid));} 1128 PetscFunctionReturn(0); 1129 #else 1130 SETERRQ(comm, PETSC_ERR_SUP, "This method requires ExodusII support. Reconfigure using --download-exodusii"); 1131 #endif 1132 } 1133 1134 #if defined(PETSC_HAVE_EXODUSII) 1135 static PetscErrorCode ExodusGetCellType_Private(const char *elem_type, DMPolytopeType *ct) 1136 { 1137 PetscBool flg; 1138 PetscErrorCode ierr; 1139 1140 PetscFunctionBegin; 1141 *ct = DM_POLYTOPE_UNKNOWN; 1142 ierr = PetscStrcmp(elem_type, "TRI", &flg);CHKERRQ(ierr); 1143 if (flg) {*ct = DM_POLYTOPE_TRIANGLE; goto done;} 1144 ierr = PetscStrcmp(elem_type, "TRI3", &flg);CHKERRQ(ierr); 1145 if (flg) {*ct = DM_POLYTOPE_TRIANGLE; goto done;} 1146 ierr = PetscStrcmp(elem_type, "QUAD", &flg);CHKERRQ(ierr); 1147 if (flg) {*ct = DM_POLYTOPE_QUADRILATERAL; goto done;} 1148 ierr = PetscStrcmp(elem_type, "QUAD4", &flg);CHKERRQ(ierr); 1149 if (flg) {*ct = DM_POLYTOPE_QUADRILATERAL; goto done;} 1150 ierr = PetscStrcmp(elem_type, "SHELL4", &flg);CHKERRQ(ierr); 1151 if (flg) {*ct = DM_POLYTOPE_QUADRILATERAL; goto done;} 1152 ierr = PetscStrcmp(elem_type, "TETRA", &flg);CHKERRQ(ierr); 1153 if (flg) {*ct = DM_POLYTOPE_TETRAHEDRON; goto done;} 1154 ierr = PetscStrcmp(elem_type, "TET4", &flg);CHKERRQ(ierr); 1155 if (flg) {*ct = DM_POLYTOPE_TETRAHEDRON; goto done;} 1156 ierr = PetscStrcmp(elem_type, "WEDGE", &flg);CHKERRQ(ierr); 1157 if (flg) {*ct = DM_POLYTOPE_TRI_PRISM; goto done;} 1158 ierr = PetscStrcmp(elem_type, "HEX", &flg);CHKERRQ(ierr); 1159 if (flg) {*ct = DM_POLYTOPE_HEXAHEDRON; goto done;} 1160 ierr = PetscStrcmp(elem_type, "HEX8", &flg);CHKERRQ(ierr); 1161 if (flg) {*ct = DM_POLYTOPE_HEXAHEDRON; goto done;} 1162 ierr = PetscStrcmp(elem_type, "HEXAHEDRON", &flg);CHKERRQ(ierr); 1163 if (flg) {*ct = DM_POLYTOPE_HEXAHEDRON; goto done;} 1164 SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unrecognized element type %s", elem_type); 1165 done: 1166 PetscFunctionReturn(0); 1167 } 1168 #endif 1169 1170 /*@ 1171 DMPlexCreateExodus - Create a DMPlex mesh from an ExodusII file ID. 1172 1173 Collective 1174 1175 Input Parameters: 1176 + comm - The MPI communicator 1177 . exoid - The ExodusII id associated with a exodus file and obtained using ex_open 1178 - interpolate - Create faces and edges in the mesh 1179 1180 Output Parameter: 1181 . dm - The DM object representing the mesh 1182 1183 Level: beginner 1184 1185 .seealso: DMPLEX, DMCreate() 1186 @*/ 1187 PetscErrorCode DMPlexCreateExodus(MPI_Comm comm, PetscInt exoid, PetscBool interpolate, DM *dm) 1188 { 1189 #if defined(PETSC_HAVE_EXODUSII) 1190 PetscMPIInt num_proc, rank; 1191 PetscSection coordSection; 1192 Vec coordinates; 1193 PetscScalar *coords; 1194 PetscInt coordSize, v; 1195 PetscErrorCode ierr; 1196 /* Read from ex_get_init() */ 1197 char title[PETSC_MAX_PATH_LEN+1]; 1198 int dim = 0, dimEmbed = 0, numVertices = 0, numCells = 0, numHybridCells = 0; 1199 int num_cs = 0, num_vs = 0, num_fs = 0; 1200 #endif 1201 1202 PetscFunctionBegin; 1203 #if defined(PETSC_HAVE_EXODUSII) 1204 ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr); 1205 ierr = MPI_Comm_size(comm, &num_proc);CHKERRQ(ierr); 1206 ierr = DMCreate(comm, dm);CHKERRQ(ierr); 1207 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr); 1208 /* Open EXODUS II file and read basic informations on rank 0, then broadcast to all processors */ 1209 if (!rank) { 1210 ierr = PetscMemzero(title,PETSC_MAX_PATH_LEN+1);CHKERRQ(ierr); 1211 PetscStackCallStandard(ex_get_init,(exoid, title, &dimEmbed, &numVertices, &numCells, &num_cs, &num_vs, &num_fs)); 1212 if (!num_cs) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Exodus file does not contain any cell set\n"); 1213 } 1214 ierr = MPI_Bcast(title, PETSC_MAX_PATH_LEN+1, MPI_CHAR, 0, comm);CHKERRQ(ierr); 1215 ierr = MPI_Bcast(&dim, 1, MPI_INT, 0, comm);CHKERRQ(ierr); 1216 ierr = PetscObjectSetName((PetscObject) *dm, title);CHKERRQ(ierr); 1217 ierr = DMPlexSetChart(*dm, 0, numCells+numVertices);CHKERRQ(ierr); 1218 /* We do not want this label automatically computed, instead we compute it here */ 1219 ierr = DMCreateLabel(*dm, "celltype");CHKERRQ(ierr); 1220 1221 /* Read cell sets information */ 1222 if (!rank) { 1223 PetscInt *cone; 1224 int c, cs, ncs, c_loc, v, v_loc; 1225 /* Read from ex_get_elem_blk_ids() */ 1226 int *cs_id, *cs_order; 1227 /* Read from ex_get_elem_block() */ 1228 char buffer[PETSC_MAX_PATH_LEN+1]; 1229 int num_cell_in_set, num_vertex_per_cell, num_hybrid, num_attr; 1230 /* Read from ex_get_elem_conn() */ 1231 int *cs_connect; 1232 1233 /* Get cell sets IDs */ 1234 ierr = PetscMalloc2(num_cs, &cs_id, num_cs, &cs_order);CHKERRQ(ierr); 1235 PetscStackCallStandard(ex_get_ids,(exoid, EX_ELEM_BLOCK, cs_id)); 1236 /* Read the cell set connectivity table and build mesh topology 1237 EXO standard requires that cells in cell sets be numbered sequentially and be pairwise disjoint. */ 1238 /* Check for a hybrid mesh */ 1239 for (cs = 0, num_hybrid = 0; cs < num_cs; ++cs) { 1240 DMPolytopeType ct; 1241 char elem_type[PETSC_MAX_PATH_LEN]; 1242 1243 ierr = PetscArrayzero(elem_type, sizeof(elem_type));CHKERRQ(ierr); 1244 PetscStackCallStandard(ex_get_elem_type,(exoid, cs_id[cs], elem_type)); 1245 ierr = ExodusGetCellType_Private(elem_type, &ct);CHKERRQ(ierr); 1246 dim = PetscMax(dim, DMPolytopeTypeGetDim(ct)); 1247 PetscStackCallStandard(ex_get_block,(exoid, EX_ELEM_BLOCK, cs_id[cs], buffer, &num_cell_in_set,&num_vertex_per_cell, 0, 0, &num_attr)); 1248 switch (ct) { 1249 case DM_POLYTOPE_TRI_PRISM: 1250 cs_order[cs] = cs; 1251 numHybridCells += num_cell_in_set; 1252 ++num_hybrid; 1253 break; 1254 default: 1255 for (c = cs; c > cs-num_hybrid; --c) cs_order[c] = cs_order[c-1]; 1256 cs_order[cs-num_hybrid] = cs; 1257 } 1258 } 1259 /* First set sizes */ 1260 for (ncs = 0, c = 0; ncs < num_cs; ++ncs) { 1261 DMPolytopeType ct; 1262 char elem_type[PETSC_MAX_PATH_LEN]; 1263 const PetscInt cs = cs_order[ncs]; 1264 1265 ierr = PetscArrayzero(elem_type, sizeof(elem_type));CHKERRQ(ierr); 1266 PetscStackCallStandard(ex_get_elem_type,(exoid, cs_id[cs], elem_type)); 1267 ierr = ExodusGetCellType_Private(elem_type, &ct);CHKERRQ(ierr); 1268 PetscStackCallStandard(ex_get_block,(exoid, EX_ELEM_BLOCK, cs_id[cs], buffer, &num_cell_in_set,&num_vertex_per_cell, 0, 0, &num_attr)); 1269 for (c_loc = 0; c_loc < num_cell_in_set; ++c_loc, ++c) { 1270 ierr = DMPlexSetConeSize(*dm, c, num_vertex_per_cell);CHKERRQ(ierr); 1271 ierr = DMPlexSetCellType(*dm, c, ct);CHKERRQ(ierr); 1272 } 1273 } 1274 for (v = numCells; v < numCells+numVertices; ++v) {ierr = DMPlexSetCellType(*dm, v, DM_POLYTOPE_POINT);CHKERRQ(ierr);} 1275 ierr = DMSetUp(*dm);CHKERRQ(ierr); 1276 for (ncs = 0, c = 0; ncs < num_cs; ++ncs) { 1277 const PetscInt cs = cs_order[ncs]; 1278 PetscStackCallStandard(ex_get_block,(exoid, EX_ELEM_BLOCK, cs_id[cs], buffer, &num_cell_in_set, &num_vertex_per_cell, 0, 0, &num_attr)); 1279 ierr = PetscMalloc2(num_vertex_per_cell*num_cell_in_set,&cs_connect,num_vertex_per_cell,&cone);CHKERRQ(ierr); 1280 PetscStackCallStandard(ex_get_conn,(exoid, EX_ELEM_BLOCK, cs_id[cs], cs_connect,NULL,NULL)); 1281 /* EXO uses Fortran-based indexing, DMPlex uses C-style and numbers cell first then vertices. */ 1282 for (c_loc = 0, v = 0; c_loc < num_cell_in_set; ++c_loc, ++c) { 1283 DMPolytopeType ct; 1284 1285 for (v_loc = 0; v_loc < num_vertex_per_cell; ++v_loc, ++v) { 1286 cone[v_loc] = cs_connect[v]+numCells-1; 1287 } 1288 ierr = DMPlexGetCellType(*dm, c, &ct);CHKERRQ(ierr); 1289 ierr = DMPlexInvertCell(ct, cone);CHKERRQ(ierr); 1290 ierr = DMPlexSetCone(*dm, c, cone);CHKERRQ(ierr); 1291 ierr = DMSetLabelValue(*dm, "Cell Sets", c, cs_id[cs]);CHKERRQ(ierr); 1292 } 1293 ierr = PetscFree2(cs_connect,cone);CHKERRQ(ierr); 1294 } 1295 ierr = PetscFree2(cs_id, cs_order);CHKERRQ(ierr); 1296 } 1297 { 1298 PetscInt ints[] = {dim, dimEmbed}; 1299 1300 ierr = MPI_Bcast(ints, 2, MPIU_INT, 0, comm);CHKERRQ(ierr); 1301 ierr = DMSetDimension(*dm, ints[0]);CHKERRQ(ierr); 1302 ierr = DMSetCoordinateDim(*dm, ints[1]);CHKERRQ(ierr); 1303 dim = ints[0]; 1304 dimEmbed = ints[1]; 1305 } 1306 ierr = DMPlexSymmetrize(*dm);CHKERRQ(ierr); 1307 ierr = DMPlexStratify(*dm);CHKERRQ(ierr); 1308 if (interpolate) { 1309 DM idm; 1310 1311 ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr); 1312 ierr = DMDestroy(dm);CHKERRQ(ierr); 1313 *dm = idm; 1314 } 1315 1316 /* Create vertex set label */ 1317 if (!rank && (num_vs > 0)) { 1318 int vs, v; 1319 /* Read from ex_get_node_set_ids() */ 1320 int *vs_id; 1321 /* Read from ex_get_node_set_param() */ 1322 int num_vertex_in_set; 1323 /* Read from ex_get_node_set() */ 1324 int *vs_vertex_list; 1325 1326 /* Get vertex set ids */ 1327 ierr = PetscMalloc1(num_vs, &vs_id);CHKERRQ(ierr); 1328 PetscStackCallStandard(ex_get_ids,(exoid, EX_NODE_SET, vs_id)); 1329 for (vs = 0; vs < num_vs; ++vs) { 1330 PetscStackCallStandard(ex_get_set_param,(exoid, EX_NODE_SET, vs_id[vs], &num_vertex_in_set, NULL)); 1331 ierr = PetscMalloc1(num_vertex_in_set, &vs_vertex_list);CHKERRQ(ierr); 1332 PetscStackCallStandard(ex_get_set,(exoid, EX_NODE_SET, vs_id[vs], vs_vertex_list, NULL)); 1333 for (v = 0; v < num_vertex_in_set; ++v) { 1334 ierr = DMSetLabelValue(*dm, "Vertex Sets", vs_vertex_list[v]+numCells-1, vs_id[vs]);CHKERRQ(ierr); 1335 } 1336 ierr = PetscFree(vs_vertex_list);CHKERRQ(ierr); 1337 } 1338 ierr = PetscFree(vs_id);CHKERRQ(ierr); 1339 } 1340 /* Read coordinates */ 1341 ierr = DMGetCoordinateSection(*dm, &coordSection);CHKERRQ(ierr); 1342 ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr); 1343 ierr = PetscSectionSetFieldComponents(coordSection, 0, dimEmbed);CHKERRQ(ierr); 1344 ierr = PetscSectionSetChart(coordSection, numCells, numCells + numVertices);CHKERRQ(ierr); 1345 for (v = numCells; v < numCells+numVertices; ++v) { 1346 ierr = PetscSectionSetDof(coordSection, v, dimEmbed);CHKERRQ(ierr); 1347 ierr = PetscSectionSetFieldDof(coordSection, v, 0, dimEmbed);CHKERRQ(ierr); 1348 } 1349 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 1350 ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr); 1351 ierr = VecCreate(PETSC_COMM_SELF, &coordinates);CHKERRQ(ierr); 1352 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr); 1353 ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 1354 ierr = VecSetBlockSize(coordinates, dimEmbed);CHKERRQ(ierr); 1355 ierr = VecSetType(coordinates,VECSTANDARD);CHKERRQ(ierr); 1356 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 1357 if (!rank) { 1358 PetscReal *x, *y, *z; 1359 1360 ierr = PetscMalloc3(numVertices,&x,numVertices,&y,numVertices,&z);CHKERRQ(ierr); 1361 PetscStackCallStandard(ex_get_coord,(exoid, x, y, z)); 1362 if (dimEmbed > 0) { 1363 for (v = 0; v < numVertices; ++v) coords[v*dimEmbed+0] = x[v]; 1364 } 1365 if (dimEmbed > 1) { 1366 for (v = 0; v < numVertices; ++v) coords[v*dimEmbed+1] = y[v]; 1367 } 1368 if (dimEmbed > 2) { 1369 for (v = 0; v < numVertices; ++v) coords[v*dimEmbed+2] = z[v]; 1370 } 1371 ierr = PetscFree3(x,y,z);CHKERRQ(ierr); 1372 } 1373 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 1374 ierr = DMSetCoordinatesLocal(*dm, coordinates);CHKERRQ(ierr); 1375 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 1376 1377 /* Create side set label */ 1378 if (!rank && interpolate && (num_fs > 0)) { 1379 int fs, f, voff; 1380 /* Read from ex_get_side_set_ids() */ 1381 int *fs_id; 1382 /* Read from ex_get_side_set_param() */ 1383 int num_side_in_set; 1384 /* Read from ex_get_side_set_node_list() */ 1385 int *fs_vertex_count_list, *fs_vertex_list; 1386 /* Read side set labels */ 1387 char fs_name[MAX_STR_LENGTH+1]; 1388 1389 /* Get side set ids */ 1390 ierr = PetscMalloc1(num_fs, &fs_id);CHKERRQ(ierr); 1391 PetscStackCallStandard(ex_get_ids,(exoid, EX_SIDE_SET, fs_id)); 1392 for (fs = 0; fs < num_fs; ++fs) { 1393 PetscStackCallStandard(ex_get_set_param,(exoid, EX_SIDE_SET, fs_id[fs], &num_side_in_set, NULL)); 1394 ierr = PetscMalloc2(num_side_in_set,&fs_vertex_count_list,num_side_in_set*4,&fs_vertex_list);CHKERRQ(ierr); 1395 PetscStackCallStandard(ex_get_side_set_node_list,(exoid, fs_id[fs], fs_vertex_count_list, fs_vertex_list)); 1396 /* Get the specific name associated with this side set ID. */ 1397 int fs_name_err = ex_get_name(exoid, EX_SIDE_SET, fs_id[fs], fs_name); 1398 for (f = 0, voff = 0; f < num_side_in_set; ++f) { 1399 const PetscInt *faces = NULL; 1400 PetscInt faceSize = fs_vertex_count_list[f], numFaces; 1401 PetscInt faceVertices[4], v; 1402 1403 if (faceSize > 4) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "ExodusII side cannot have %d > 4 vertices", faceSize); 1404 for (v = 0; v < faceSize; ++v, ++voff) { 1405 faceVertices[v] = fs_vertex_list[voff]+numCells-1; 1406 } 1407 ierr = DMPlexGetFullJoin(*dm, faceSize, faceVertices, &numFaces, &faces);CHKERRQ(ierr); 1408 if (numFaces != 1) SETERRQ3(comm, PETSC_ERR_ARG_WRONG, "Invalid ExodusII side %d in set %d maps to %d faces", f, fs, numFaces); 1409 ierr = DMSetLabelValue(*dm, "Face Sets", faces[0], fs_id[fs]);CHKERRQ(ierr); 1410 /* Only add the label if one has been detected for this side set. */ 1411 if (!fs_name_err) { 1412 ierr = DMSetLabelValue(*dm, fs_name, faces[0], fs_id[fs]);CHKERRQ(ierr); 1413 } 1414 ierr = DMPlexRestoreJoin(*dm, faceSize, faceVertices, &numFaces, &faces);CHKERRQ(ierr); 1415 } 1416 ierr = PetscFree2(fs_vertex_count_list,fs_vertex_list);CHKERRQ(ierr); 1417 } 1418 ierr = PetscFree(fs_id);CHKERRQ(ierr); 1419 } 1420 PetscFunctionReturn(0); 1421 #else 1422 SETERRQ(comm, PETSC_ERR_SUP, "This method requires ExodusII support. Reconfigure using --download-exodusii"); 1423 #endif 1424 } 1425