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