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