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