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