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++) PetscCall(PetscFree(exo->zonalVariableNames[i])); 83 PetscCall(PetscFree(exo->zonalVariableNames)); 84 for (PetscInt i = 0; i < exo->numNodalVariables; i++) PetscCall(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 break; 893 case 2: 894 if (closureSize == 3 * dim) { 895 type[cs] = TRI; 896 } else if (closureSize == 4 * dim) { 897 type[cs] = QUAD; 898 } 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); 899 break; 900 case 3: 901 if (closureSize == 4 * dim) { 902 type[cs] = TET; 903 } else if (closureSize == 8 * dim) { 904 type[cs] = HEX; 905 } 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); 906 break; 907 default: 908 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Dimension %" PetscInt_FMT " not handled by ExodusII viewer", dim); 909 } 910 if ((degree == 2) && (type[cs] == SEGMENT)) numNodes += csSize; 911 if ((degree == 2) && (type[cs] == QUAD)) numNodes += csSize; 912 if ((degree == 2) && (type[cs] == HEX)) { 913 numNodes += csSize; 914 numNodes += numFaces; 915 } 916 PetscCall(DMPlexVecRestoreClosure(cdm, NULL, coord, cells[0], &closureSize, &xyz)); 917 /* Set nodes and Element type */ 918 if (type[cs] == SEGMENT) { 919 if (degree == 1) nodes[cs] = nodesLineP1; 920 else if (degree == 2) nodes[cs] = nodesLineP2; 921 } else if (type[cs] == TRI) { 922 if (degree == 1) nodes[cs] = nodesTriP1; 923 else if (degree == 2) nodes[cs] = nodesTriP2; 924 } else if (type[cs] == QUAD) { 925 if (degree == 1) nodes[cs] = nodesQuadP1; 926 else if (degree == 2) nodes[cs] = nodesQuadP2; 927 } else if (type[cs] == TET) { 928 if (degree == 1) nodes[cs] = nodesTetP1; 929 else if (degree == 2) nodes[cs] = nodesTetP2; 930 } else if (type[cs] == HEX) { 931 if (degree == 1) nodes[cs] = nodesHexP1; 932 else if (degree == 2) nodes[cs] = nodesHexP2; 933 } 934 /* Compute the number of cells not in the connectivity table */ 935 cellsNotInConnectivity -= nodes[cs][3] * csSize; 936 937 PetscCall(ISRestoreIndices(stratumIS, &cells)); 938 PetscCall(ISDestroy(&stratumIS)); 939 } 940 if (num_cs) PetscCallExternal(ex_put_init, exo->exoid, dmName, dim, numNodes, numCells, num_cs, num_vs, num_fs); 941 /* --- Connectivity --- */ 942 for (cs = 0; cs < num_cs; ++cs) { 943 IS stratumIS; 944 const PetscInt *cells; 945 PetscInt *connect, off = 0; 946 PetscInt edgesInClosure = 0, facesInClosure = 0, verticesInClosure = 0; 947 PetscInt csSize, c, connectSize, closureSize; 948 char *elem_type = NULL; 949 char elem_type_bar2[] = "BAR2", elem_type_bar3[] = "BAR3"; 950 char elem_type_tri3[] = "TRI3", elem_type_quad4[] = "QUAD4"; 951 char elem_type_tri6[] = "TRI6", elem_type_quad9[] = "QUAD9"; 952 char elem_type_tet4[] = "TET4", elem_type_hex8[] = "HEX8"; 953 char elem_type_tet10[] = "TET10", elem_type_hex27[] = "HEX27"; 954 955 PetscCall(DMLabelGetStratumIS(csLabel, csIdx[cs], &stratumIS)); 956 PetscCall(ISGetIndices(stratumIS, &cells)); 957 PetscCall(ISGetSize(stratumIS, &csSize)); 958 /* Set Element type */ 959 if (type[cs] == SEGMENT) { 960 if (degree == 1) elem_type = elem_type_bar2; 961 else if (degree == 2) elem_type = elem_type_bar3; 962 } else if (type[cs] == TRI) { 963 if (degree == 1) elem_type = elem_type_tri3; 964 else if (degree == 2) elem_type = elem_type_tri6; 965 } else if (type[cs] == QUAD) { 966 if (degree == 1) elem_type = elem_type_quad4; 967 else if (degree == 2) elem_type = elem_type_quad9; 968 } else if (type[cs] == TET) { 969 if (degree == 1) elem_type = elem_type_tet4; 970 else if (degree == 2) elem_type = elem_type_tet10; 971 } else if (type[cs] == HEX) { 972 if (degree == 1) elem_type = elem_type_hex8; 973 else if (degree == 2) elem_type = elem_type_hex27; 974 } 975 connectSize = nodes[cs][0] + nodes[cs][1] + nodes[cs][2] + nodes[cs][3]; 976 PetscCall(PetscMalloc1(PetscMax(27, connectSize) * csSize, &connect)); 977 PetscCallExternal(ex_put_block, exo->exoid, EX_ELEM_BLOCK, csIdx[cs], elem_type, csSize, connectSize, 0, 0, 1); 978 /* Find number of vertices, edges, and faces in the closure */ 979 verticesInClosure = nodes[cs][0]; 980 if (depth > 1) { 981 if (dim == 2) { 982 PetscCall(DMPlexGetConeSize(dm, cells[0], &edgesInClosure)); 983 } else if (dim == 3) { 984 PetscInt *closure = NULL; 985 986 PetscCall(DMPlexGetConeSize(dm, cells[0], &facesInClosure)); 987 PetscCall(DMPlexGetTransitiveClosure(dm, cells[0], PETSC_TRUE, &closureSize, &closure)); 988 edgesInClosure = closureSize - facesInClosure - 1 - verticesInClosure; 989 PetscCall(DMPlexRestoreTransitiveClosure(dm, cells[0], PETSC_TRUE, &closureSize, &closure)); 990 } 991 } 992 /* Get connectivity for each cell */ 993 for (c = 0; c < csSize; ++c) { 994 PetscInt *closure = NULL; 995 PetscInt temp, i; 996 997 PetscCall(DMPlexGetTransitiveClosure(dm, cells[c], PETSC_TRUE, &closureSize, &closure)); 998 for (i = 0; i < connectSize; ++i) { 999 if (i < nodes[cs][0]) { /* Vertices */ 1000 connect[i + off] = closure[(i + edgesInClosure + facesInClosure + 1) * 2] + 1; 1001 connect[i + off] -= cellsNotInConnectivity; 1002 } else if (i < nodes[cs][0] + nodes[cs][1]) { /* Edges */ 1003 connect[i + off] = closure[(i - verticesInClosure + facesInClosure + 1) * 2] + 1; 1004 if (nodes[cs][2] == 0) connect[i + off] -= numFaces; 1005 connect[i + off] -= cellsNotInConnectivity; 1006 } else if (i < nodes[cs][0] + nodes[cs][1] + nodes[cs][3]) { /* Cells */ 1007 connect[i + off] = closure[0] + 1; 1008 connect[i + off] -= skipCells; 1009 } else if (i < nodes[cs][0] + nodes[cs][1] + nodes[cs][3] + nodes[cs][2]) { /* Faces */ 1010 connect[i + off] = closure[(i - edgesInClosure - verticesInClosure) * 2] + 1; 1011 connect[i + off] -= cellsNotInConnectivity; 1012 } else { 1013 connect[i + off] = -1; 1014 } 1015 } 1016 /* Tetrahedra are inverted */ 1017 if (type[cs] == TET) { 1018 temp = connect[0 + off]; 1019 connect[0 + off] = connect[1 + off]; 1020 connect[1 + off] = temp; 1021 if (degree == 2) { 1022 temp = connect[5 + off]; 1023 connect[5 + off] = connect[6 + off]; 1024 connect[6 + off] = temp; 1025 temp = connect[7 + off]; 1026 connect[7 + off] = connect[8 + off]; 1027 connect[8 + off] = temp; 1028 } 1029 } 1030 /* Hexahedra are inverted */ 1031 if (type[cs] == HEX) { 1032 temp = connect[1 + off]; 1033 connect[1 + off] = connect[3 + off]; 1034 connect[3 + off] = temp; 1035 if (degree == 2) { 1036 temp = connect[8 + off]; 1037 connect[8 + off] = connect[11 + off]; 1038 connect[11 + off] = temp; 1039 temp = connect[9 + off]; 1040 connect[9 + off] = connect[10 + off]; 1041 connect[10 + off] = temp; 1042 temp = connect[16 + off]; 1043 connect[16 + off] = connect[17 + off]; 1044 connect[17 + off] = temp; 1045 temp = connect[18 + off]; 1046 connect[18 + off] = connect[19 + off]; 1047 connect[19 + off] = temp; 1048 1049 temp = connect[12 + off]; 1050 connect[12 + off] = connect[16 + off]; 1051 connect[16 + off] = temp; 1052 temp = connect[13 + off]; 1053 connect[13 + off] = connect[17 + off]; 1054 connect[17 + off] = temp; 1055 temp = connect[14 + off]; 1056 connect[14 + off] = connect[18 + off]; 1057 connect[18 + off] = temp; 1058 temp = connect[15 + off]; 1059 connect[15 + off] = connect[19 + off]; 1060 connect[19 + off] = temp; 1061 1062 temp = connect[23 + off]; 1063 connect[23 + off] = connect[26 + off]; 1064 connect[26 + off] = temp; 1065 temp = connect[24 + off]; 1066 connect[24 + off] = connect[25 + off]; 1067 connect[25 + off] = temp; 1068 temp = connect[25 + off]; 1069 connect[25 + off] = connect[26 + off]; 1070 connect[26 + off] = temp; 1071 } 1072 } 1073 off += connectSize; 1074 PetscCall(DMPlexRestoreTransitiveClosure(dm, cells[c], PETSC_TRUE, &closureSize, &closure)); 1075 } 1076 PetscCallExternal(ex_put_conn, exo->exoid, EX_ELEM_BLOCK, csIdx[cs], connect, 0, 0); 1077 skipCells += (nodes[cs][3] == 0) * csSize; 1078 PetscCall(PetscFree(connect)); 1079 PetscCall(ISRestoreIndices(stratumIS, &cells)); 1080 PetscCall(ISDestroy(&stratumIS)); 1081 } 1082 PetscCall(PetscFree(type)); 1083 /* --- Coordinates --- */ 1084 PetscCall(PetscSectionSetChart(coordSection, pStart, pEnd)); 1085 if (num_cs) { 1086 for (d = 0; d < depth; ++d) { 1087 PetscCall(DMPlexGetDepthStratum(dm, d, &pStart, &pEnd)); 1088 for (p = pStart; p < pEnd; ++p) PetscCall(PetscSectionSetDof(coordSection, p, nodes[0][d] > 0)); 1089 } 1090 } 1091 for (cs = 0; cs < num_cs; ++cs) { 1092 IS stratumIS; 1093 const PetscInt *cells; 1094 PetscInt csSize, c; 1095 1096 PetscCall(DMLabelGetStratumIS(csLabel, csIdx[cs], &stratumIS)); 1097 PetscCall(ISGetIndices(stratumIS, &cells)); 1098 PetscCall(ISGetSize(stratumIS, &csSize)); 1099 for (c = 0; c < csSize; ++c) PetscCall(PetscSectionSetDof(coordSection, cells[c], nodes[cs][3] > 0)); 1100 PetscCall(ISRestoreIndices(stratumIS, &cells)); 1101 PetscCall(ISDestroy(&stratumIS)); 1102 } 1103 if (num_cs) { 1104 PetscCall(ISRestoreIndices(csIS, &csIdx)); 1105 PetscCall(ISDestroy(&csIS)); 1106 } 1107 PetscCall(PetscFree(nodes)); 1108 PetscCall(PetscSectionSetUp(coordSection)); 1109 if (numNodes) { 1110 const char *coordNames[3] = {"x", "y", "z"}; 1111 PetscScalar *closure, *cval; 1112 PetscReal *coords; 1113 PetscInt hasDof, n = 0; 1114 1115 /* There can't be more than 24 values in the closure of a point for the coord coordSection */ 1116 PetscCall(PetscCalloc3(numNodes * 3, &coords, dim, &cval, 24, &closure)); 1117 PetscCall(DMGetCoordinatesLocalNoncollective(dm, &coord)); 1118 PetscCall(DMPlexGetChart(dm, &pStart, &pEnd)); 1119 for (p = pStart; p < pEnd; ++p) { 1120 PetscCall(PetscSectionGetDof(coordSection, p, &hasDof)); 1121 if (hasDof) { 1122 PetscInt closureSize = 24, j; 1123 1124 PetscCall(DMPlexVecGetClosure(cdm, NULL, coord, p, &closureSize, &closure)); 1125 for (d = 0; d < dim; ++d) { 1126 cval[d] = 0.0; 1127 for (j = 0; j < closureSize / dim; j++) cval[d] += closure[j * dim + d]; 1128 coords[d * numNodes + n] = PetscRealPart(cval[d]) * dim / closureSize; 1129 } 1130 ++n; 1131 } 1132 } 1133 PetscCallExternal(ex_put_coord, exo->exoid, &coords[0 * numNodes], &coords[1 * numNodes], &coords[2 * numNodes]); 1134 PetscCall(PetscFree3(coords, cval, closure)); 1135 PetscCallExternal(ex_put_coord_names, exo->exoid, (char **)coordNames); 1136 } 1137 1138 /* --- Node Sets/Vertex Sets --- */ 1139 PetscCall(DMHasLabel(dm, "Vertex Sets", &hasLabel)); 1140 if (hasLabel) { 1141 PetscInt i, vs, vsSize; 1142 const PetscInt *vsIdx, *vertices; 1143 PetscInt *nodeList; 1144 IS vsIS, stratumIS; 1145 DMLabel vsLabel; 1146 PetscCall(DMGetLabel(dm, "Vertex Sets", &vsLabel)); 1147 PetscCall(DMLabelGetValueIS(vsLabel, &vsIS)); 1148 PetscCall(ISGetIndices(vsIS, &vsIdx)); 1149 for (vs = 0; vs < num_vs; ++vs) { 1150 PetscCall(DMLabelGetStratumIS(vsLabel, vsIdx[vs], &stratumIS)); 1151 PetscCall(ISGetIndices(stratumIS, &vertices)); 1152 PetscCall(ISGetSize(stratumIS, &vsSize)); 1153 PetscCall(PetscMalloc1(vsSize, &nodeList)); 1154 for (i = 0; i < vsSize; ++i) nodeList[i] = vertices[i] - skipCells + 1; 1155 PetscCallExternal(ex_put_set_param, exo->exoid, EX_NODE_SET, vsIdx[vs], vsSize, 0); 1156 PetscCallExternal(ex_put_set, exo->exoid, EX_NODE_SET, vsIdx[vs], nodeList, NULL); 1157 PetscCall(ISRestoreIndices(stratumIS, &vertices)); 1158 PetscCall(ISDestroy(&stratumIS)); 1159 PetscCall(PetscFree(nodeList)); 1160 } 1161 PetscCall(ISRestoreIndices(vsIS, &vsIdx)); 1162 PetscCall(ISDestroy(&vsIS)); 1163 } 1164 /* --- Side Sets/Face Sets --- */ 1165 PetscCall(DMHasLabel(dm, "Face Sets", &hasLabel)); 1166 if (hasLabel) { 1167 PetscInt i, j, fs, fsSize; 1168 const PetscInt *fsIdx, *faces; 1169 IS fsIS, stratumIS; 1170 DMLabel fsLabel; 1171 PetscInt numPoints, *points; 1172 PetscInt elem_list_size = 0; 1173 PetscInt *elem_list, *elem_ind, *side_list; 1174 1175 PetscCall(DMGetLabel(dm, "Face Sets", &fsLabel)); 1176 /* Compute size of Node List and Element List */ 1177 PetscCall(DMLabelGetValueIS(fsLabel, &fsIS)); 1178 PetscCall(ISGetIndices(fsIS, &fsIdx)); 1179 for (fs = 0; fs < num_fs; ++fs) { 1180 PetscCall(DMLabelGetStratumIS(fsLabel, fsIdx[fs], &stratumIS)); 1181 PetscCall(ISGetSize(stratumIS, &fsSize)); 1182 elem_list_size += fsSize; 1183 PetscCall(ISDestroy(&stratumIS)); 1184 } 1185 if (num_fs) { 1186 PetscCall(PetscMalloc3(num_fs, &elem_ind, elem_list_size, &elem_list, elem_list_size, &side_list)); 1187 elem_ind[0] = 0; 1188 for (fs = 0; fs < num_fs; ++fs) { 1189 PetscCall(DMLabelGetStratumIS(fsLabel, fsIdx[fs], &stratumIS)); 1190 PetscCall(ISGetIndices(stratumIS, &faces)); 1191 PetscCall(ISGetSize(stratumIS, &fsSize)); 1192 /* Set Parameters */ 1193 PetscCallExternal(ex_put_set_param, exo->exoid, EX_SIDE_SET, fsIdx[fs], fsSize, 0); 1194 /* Indices */ 1195 if (fs < num_fs - 1) elem_ind[fs + 1] = elem_ind[fs] + fsSize; 1196 1197 for (i = 0; i < fsSize; ++i) { 1198 /* Element List */ 1199 points = NULL; 1200 PetscCall(DMPlexGetTransitiveClosure(dm, faces[i], PETSC_FALSE, &numPoints, &points)); 1201 elem_list[elem_ind[fs] + i] = points[2] + 1; 1202 PetscCall(DMPlexRestoreTransitiveClosure(dm, faces[i], PETSC_FALSE, &numPoints, &points)); 1203 1204 /* Side List */ 1205 points = NULL; 1206 PetscCall(DMPlexGetTransitiveClosure(dm, elem_list[elem_ind[fs] + i] - 1, PETSC_TRUE, &numPoints, &points)); 1207 for (j = 1; j < numPoints; ++j) { 1208 if (points[j * 2] == faces[i]) break; 1209 } 1210 /* Convert HEX sides */ 1211 if (numPoints == 27) { 1212 if (j == 1) { 1213 j = 5; 1214 } else if (j == 2) { 1215 j = 6; 1216 } else if (j == 3) { 1217 j = 1; 1218 } else if (j == 4) { 1219 j = 3; 1220 } else if (j == 5) { 1221 j = 2; 1222 } else if (j == 6) { 1223 j = 4; 1224 } 1225 } 1226 /* Convert TET sides */ 1227 if (numPoints == 15) { 1228 --j; 1229 if (j == 0) j = 4; 1230 } 1231 side_list[elem_ind[fs] + i] = j; 1232 PetscCall(DMPlexRestoreTransitiveClosure(dm, elem_list[elem_ind[fs] + i] - 1, PETSC_TRUE, &numPoints, &points)); 1233 } 1234 PetscCall(ISRestoreIndices(stratumIS, &faces)); 1235 PetscCall(ISDestroy(&stratumIS)); 1236 } 1237 PetscCall(ISRestoreIndices(fsIS, &fsIdx)); 1238 PetscCall(ISDestroy(&fsIS)); 1239 1240 /* Put side sets */ 1241 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]]); 1242 PetscCall(PetscFree3(elem_ind, elem_list, side_list)); 1243 } 1244 } 1245 /* 1246 close the exodus file 1247 */ 1248 ex_close(exo->exoid); 1249 exo->exoid = -1; 1250 } 1251 PetscCall(PetscSectionDestroy(&coordSection)); 1252 1253 /* 1254 reopen the file in parallel 1255 */ 1256 EXO_mode = EX_WRITE; 1257 #if defined(PETSC_USE_64BIT_INDICES) 1258 EXO_mode += EX_ALL_INT64_API; 1259 #endif 1260 CPU_word_size = sizeof(PetscReal); 1261 IO_word_size = sizeof(PetscReal); 1262 exo->exoid = ex_open_par(exo->filename, EXO_mode, &CPU_word_size, &IO_word_size, &EXO_version, comm, MPI_INFO_NULL); 1263 PetscCheck(exo->exoid >= 0, PETSC_COMM_SELF, PETSC_ERR_LIB, "ex_open_par failed for %s", exo->filename); 1264 PetscFunctionReturn(PETSC_SUCCESS); 1265 } 1266 1267 static PetscErrorCode VecViewPlex_ExodusII_Nodal_Internal(Vec v, PetscExodusIIInt exoid, PetscExodusIIInt step, PetscExodusIIInt offset); 1268 static PetscErrorCode VecViewPlex_ExodusII_Zonal_Internal(Vec v, PetscExodusIIInt exoid, PetscExodusIIInt step, PetscExodusIIInt offset); 1269 static PetscErrorCode VecLoadPlex_ExodusII_Nodal_Internal(Vec v, PetscExodusIIInt exoid, PetscExodusIIInt step, PetscExodusIIInt offset); 1270 static PetscErrorCode VecLoadPlex_ExodusII_Zonal_Internal(Vec v, PetscExodusIIInt exoid, PetscExodusIIInt step, PetscExodusIIInt offset); 1271 1272 PetscErrorCode VecView_PlexExodusII_Internal(Vec v, PetscViewer viewer) 1273 { 1274 DM dm; 1275 MPI_Comm comm; 1276 PetscMPIInt rank; 1277 1278 PetscExodusIIInt exoid, offsetN = -1, offsetZ = -1; 1279 const char *vecname; 1280 PetscInt step; 1281 1282 PetscFunctionBegin; 1283 PetscCall(PetscObjectGetComm((PetscObject)v, &comm)); 1284 PetscCallMPI(MPI_Comm_rank(comm, &rank)); 1285 PetscCall(PetscViewerExodusIIGetId(viewer, &exoid)); 1286 PetscCall(VecGetDM(v, &dm)); 1287 PetscCall(PetscObjectGetName((PetscObject)v, &vecname)); 1288 1289 PetscCall(DMGetOutputSequenceNumber(dm, &step, NULL)); 1290 PetscCall(PetscViewerExodusIIGetNodalVariableIndex(viewer, vecname, &offsetN)); 1291 PetscCall(PetscViewerExodusIIGetZonalVariableIndex(viewer, vecname, &offsetZ)); 1292 PetscCheck(!(offsetN >= 0 && offsetZ >= 0), comm, PETSC_ERR_FILE_UNEXPECTED, "Found both nodal and zonal variable %s in exodus file. ", vecname); 1293 if (offsetN >= 0) { 1294 PetscCall(VecViewPlex_ExodusII_Nodal_Internal(v, exoid, (PetscExodusIIInt)step + 1, offsetN + 1)); 1295 } else if (offsetZ >= 0) { 1296 PetscCall(VecViewPlex_ExodusII_Zonal_Internal(v, exoid, (PetscExodusIIInt)step + 1, offsetZ + 1)); 1297 } else SETERRQ(comm, PETSC_ERR_FILE_UNEXPECTED, "Could not find nodal or zonal variable %s in exodus file. ", vecname); 1298 PetscFunctionReturn(PETSC_SUCCESS); 1299 } 1300 1301 PetscErrorCode VecLoad_PlexExodusII_Internal(Vec v, PetscViewer viewer) 1302 { 1303 DM dm; 1304 MPI_Comm comm; 1305 PetscMPIInt rank; 1306 1307 PetscExodusIIInt exoid, offsetN = 0, offsetZ = 0; 1308 const char *vecname; 1309 PetscInt step; 1310 1311 PetscFunctionBegin; 1312 PetscCall(PetscObjectGetComm((PetscObject)v, &comm)); 1313 PetscCallMPI(MPI_Comm_rank(comm, &rank)); 1314 PetscCall(PetscViewerExodusIIGetId(viewer, &exoid)); 1315 PetscCall(VecGetDM(v, &dm)); 1316 PetscCall(PetscObjectGetName((PetscObject)v, &vecname)); 1317 1318 PetscCall(DMGetOutputSequenceNumber(dm, &step, NULL)); 1319 PetscCall(PetscViewerExodusIIGetNodalVariableIndex(viewer, vecname, &offsetN)); 1320 PetscCall(PetscViewerExodusIIGetZonalVariableIndex(viewer, vecname, &offsetZ)); 1321 PetscCheck(!(offsetN >= 0 && offsetZ >= 0), comm, PETSC_ERR_FILE_UNEXPECTED, "Found both nodal and zonal variable %s in exodus file. ", vecname); 1322 if (offsetN >= 0) PetscCall(VecLoadPlex_ExodusII_Nodal_Internal(v, exoid, (PetscExodusIIInt)step + 1, offsetN + 1)); 1323 else if (offsetZ >= 0) PetscCall(VecLoadPlex_ExodusII_Zonal_Internal(v, exoid, (PetscExodusIIInt)step + 1, offsetZ + 1)); 1324 else SETERRQ(comm, PETSC_ERR_FILE_UNEXPECTED, "Could not find nodal or zonal variable %s in exodus file. ", vecname); 1325 PetscFunctionReturn(PETSC_SUCCESS); 1326 } 1327 1328 static PetscErrorCode VecViewPlex_ExodusII_Nodal_Internal(Vec v, PetscExodusIIInt exoid, PetscExodusIIInt step, PetscExodusIIInt offset) 1329 { 1330 MPI_Comm comm; 1331 PetscMPIInt size; 1332 DM dm; 1333 Vec vNatural, vComp; 1334 const PetscScalar *varray; 1335 PetscInt xs, xe, bs; 1336 PetscBool useNatural; 1337 1338 PetscFunctionBegin; 1339 PetscCall(PetscObjectGetComm((PetscObject)v, &comm)); 1340 PetscCallMPI(MPI_Comm_size(comm, &size)); 1341 PetscCall(VecGetDM(v, &dm)); 1342 PetscCall(DMGetUseNatural(dm, &useNatural)); 1343 useNatural = useNatural && size > 1 ? PETSC_TRUE : PETSC_FALSE; 1344 if (useNatural) { 1345 PetscCall(DMPlexCreateNaturalVector(dm, &vNatural)); 1346 PetscCall(DMPlexGlobalToNaturalBegin(dm, v, vNatural)); 1347 PetscCall(DMPlexGlobalToNaturalEnd(dm, v, vNatural)); 1348 } else { 1349 vNatural = v; 1350 } 1351 1352 /* Write local chunk of the result in the exodus file 1353 exodus stores each component of a vector-valued field as a separate variable. 1354 We assume that they are stored sequentially */ 1355 PetscCall(VecGetOwnershipRange(vNatural, &xs, &xe)); 1356 PetscCall(VecGetBlockSize(vNatural, &bs)); 1357 if (bs == 1) { 1358 PetscCall(VecGetArrayRead(vNatural, &varray)); 1359 PetscCallExternal(ex_put_partial_var, exoid, step, EX_NODAL, offset, 1, xs + 1, xe - xs, varray); 1360 PetscCall(VecRestoreArrayRead(vNatural, &varray)); 1361 } else { 1362 IS compIS; 1363 PetscInt c; 1364 1365 PetscCall(ISCreateStride(comm, (xe - xs) / bs, xs, bs, &compIS)); 1366 for (c = 0; c < bs; ++c) { 1367 PetscCall(ISStrideSetStride(compIS, (xe - xs) / bs, xs + c, bs)); 1368 PetscCall(VecGetSubVector(vNatural, compIS, &vComp)); 1369 PetscCall(VecGetArrayRead(vComp, &varray)); 1370 PetscCallExternal(ex_put_partial_var, exoid, step, EX_NODAL, offset + c, 1, xs / bs + 1, (xe - xs) / bs, varray); 1371 PetscCall(VecRestoreArrayRead(vComp, &varray)); 1372 PetscCall(VecRestoreSubVector(vNatural, compIS, &vComp)); 1373 } 1374 PetscCall(ISDestroy(&compIS)); 1375 } 1376 if (useNatural) PetscCall(VecDestroy(&vNatural)); 1377 PetscFunctionReturn(PETSC_SUCCESS); 1378 } 1379 1380 static PetscErrorCode VecLoadPlex_ExodusII_Nodal_Internal(Vec v, PetscExodusIIInt exoid, PetscExodusIIInt step, PetscExodusIIInt offset) 1381 { 1382 MPI_Comm comm; 1383 PetscMPIInt size; 1384 DM dm; 1385 Vec vNatural, vComp; 1386 PetscScalar *varray; 1387 PetscInt xs, xe, bs; 1388 PetscBool useNatural; 1389 1390 PetscFunctionBegin; 1391 PetscCall(PetscObjectGetComm((PetscObject)v, &comm)); 1392 PetscCallMPI(MPI_Comm_size(comm, &size)); 1393 PetscCall(VecGetDM(v, &dm)); 1394 PetscCall(DMGetUseNatural(dm, &useNatural)); 1395 useNatural = useNatural && size > 1 ? PETSC_TRUE : PETSC_FALSE; 1396 if (useNatural) PetscCall(DMPlexCreateNaturalVector(dm, &vNatural)); 1397 else vNatural = v; 1398 1399 /* Read local chunk from the file */ 1400 PetscCall(VecGetOwnershipRange(vNatural, &xs, &xe)); 1401 PetscCall(VecGetBlockSize(vNatural, &bs)); 1402 if (bs == 1) { 1403 PetscCall(VecGetArray(vNatural, &varray)); 1404 PetscCallExternal(ex_get_partial_var, exoid, step, EX_NODAL, offset, 1, xs + 1, xe - xs, varray); 1405 PetscCall(VecRestoreArray(vNatural, &varray)); 1406 } else { 1407 IS compIS; 1408 PetscInt c; 1409 1410 PetscCall(ISCreateStride(comm, (xe - xs) / bs, xs, bs, &compIS)); 1411 for (c = 0; c < bs; ++c) { 1412 PetscCall(ISStrideSetStride(compIS, (xe - xs) / bs, xs + c, bs)); 1413 PetscCall(VecGetSubVector(vNatural, compIS, &vComp)); 1414 PetscCall(VecGetArray(vComp, &varray)); 1415 PetscCallExternal(ex_get_partial_var, exoid, step, EX_NODAL, offset + c, 1, xs / bs + 1, (xe - xs) / bs, varray); 1416 PetscCall(VecRestoreArray(vComp, &varray)); 1417 PetscCall(VecRestoreSubVector(vNatural, compIS, &vComp)); 1418 } 1419 PetscCall(ISDestroy(&compIS)); 1420 } 1421 if (useNatural) { 1422 PetscCall(DMPlexNaturalToGlobalBegin(dm, vNatural, v)); 1423 PetscCall(DMPlexNaturalToGlobalEnd(dm, vNatural, v)); 1424 PetscCall(VecDestroy(&vNatural)); 1425 } 1426 PetscFunctionReturn(PETSC_SUCCESS); 1427 } 1428 1429 static PetscErrorCode VecViewPlex_ExodusII_Zonal_Internal(Vec v, PetscExodusIIInt exoid, PetscExodusIIInt step, PetscExodusIIInt offset) 1430 { 1431 MPI_Comm comm; 1432 PetscMPIInt size; 1433 DM dm; 1434 Vec vNatural, vComp; 1435 const PetscScalar *varray; 1436 PetscInt xs, xe, bs; 1437 PetscBool useNatural; 1438 IS compIS; 1439 PetscInt *csSize, *csID; 1440 PetscExodusIIInt numCS, set, csxs = 0; 1441 1442 PetscFunctionBegin; 1443 PetscCall(PetscObjectGetComm((PetscObject)v, &comm)); 1444 PetscCallMPI(MPI_Comm_size(comm, &size)); 1445 PetscCall(VecGetDM(v, &dm)); 1446 PetscCall(DMGetUseNatural(dm, &useNatural)); 1447 useNatural = useNatural && size > 1 ? PETSC_TRUE : PETSC_FALSE; 1448 if (useNatural) { 1449 PetscCall(DMPlexCreateNaturalVector(dm, &vNatural)); 1450 PetscCall(DMPlexGlobalToNaturalBegin(dm, v, vNatural)); 1451 PetscCall(DMPlexGlobalToNaturalEnd(dm, v, vNatural)); 1452 } else { 1453 vNatural = v; 1454 } 1455 1456 /* Write local chunk of the result in the exodus file 1457 exodus stores each component of a vector-valued field as a separate variable. 1458 We assume that they are stored sequentially 1459 Zonal variables are accessed one element block at a time, so we loop through the cell sets, 1460 but once the vector has been reordered to natural size, we cannot use the label information 1461 to figure out what to save where. */ 1462 numCS = ex_inquire_int(exoid, EX_INQ_ELEM_BLK); 1463 PetscCall(PetscMalloc2(numCS, &csID, numCS, &csSize)); 1464 PetscCallExternal(ex_get_ids, exoid, EX_ELEM_BLOCK, csID); 1465 for (set = 0; set < numCS; ++set) { 1466 ex_block block; 1467 1468 block.id = csID[set]; 1469 block.type = EX_ELEM_BLOCK; 1470 PetscCallExternal(ex_get_block_param, exoid, &block); 1471 csSize[set] = block.num_entry; 1472 } 1473 PetscCall(VecGetOwnershipRange(vNatural, &xs, &xe)); 1474 PetscCall(VecGetBlockSize(vNatural, &bs)); 1475 if (bs > 1) PetscCall(ISCreateStride(comm, (xe - xs) / bs, xs, bs, &compIS)); 1476 for (set = 0; set < numCS; set++) { 1477 PetscInt csLocalSize, c; 1478 1479 /* range of indices for set setID[set]: csxs:csxs + csSize[set]-1 1480 local slice of zonal values: xs/bs,xm/bs-1 1481 intersection: max(xs/bs,csxs),min(xm/bs-1,csxs + csSize[set]-1) */ 1482 csLocalSize = PetscMax(0, PetscMin(xe / bs, csxs + csSize[set]) - PetscMax(xs / bs, csxs)); 1483 if (bs == 1) { 1484 PetscCall(VecGetArrayRead(vNatural, &varray)); 1485 PetscCallExternal(ex_put_partial_var, exoid, step, EX_ELEM_BLOCK, offset, csID[set], PetscMax(xs - csxs, 0) + 1, csLocalSize, &varray[PetscMax(0, csxs - xs)]); 1486 PetscCall(VecRestoreArrayRead(vNatural, &varray)); 1487 } else { 1488 for (c = 0; c < bs; ++c) { 1489 PetscCall(ISStrideSetStride(compIS, (xe - xs) / bs, xs + c, bs)); 1490 PetscCall(VecGetSubVector(vNatural, compIS, &vComp)); 1491 PetscCall(VecGetArrayRead(vComp, &varray)); 1492 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)]); 1493 PetscCall(VecRestoreArrayRead(vComp, &varray)); 1494 PetscCall(VecRestoreSubVector(vNatural, compIS, &vComp)); 1495 } 1496 } 1497 csxs += csSize[set]; 1498 } 1499 PetscCall(PetscFree2(csID, csSize)); 1500 if (bs > 1) PetscCall(ISDestroy(&compIS)); 1501 if (useNatural) PetscCall(VecDestroy(&vNatural)); 1502 PetscFunctionReturn(PETSC_SUCCESS); 1503 } 1504 1505 static PetscErrorCode VecLoadPlex_ExodusII_Zonal_Internal(Vec v, PetscExodusIIInt exoid, PetscExodusIIInt step, PetscExodusIIInt offset) 1506 { 1507 MPI_Comm comm; 1508 PetscMPIInt size; 1509 DM dm; 1510 Vec vNatural, vComp; 1511 PetscScalar *varray; 1512 PetscInt xs, xe, bs; 1513 PetscBool useNatural; 1514 IS compIS; 1515 PetscInt *csSize, *csID; 1516 PetscExodusIIInt numCS, set, csxs = 0; 1517 1518 PetscFunctionBegin; 1519 PetscCall(PetscObjectGetComm((PetscObject)v, &comm)); 1520 PetscCallMPI(MPI_Comm_size(comm, &size)); 1521 PetscCall(VecGetDM(v, &dm)); 1522 PetscCall(DMGetUseNatural(dm, &useNatural)); 1523 useNatural = useNatural && size > 1 ? PETSC_TRUE : PETSC_FALSE; 1524 if (useNatural) PetscCall(DMPlexCreateNaturalVector(dm, &vNatural)); 1525 else vNatural = v; 1526 1527 /* Read local chunk of the result in the exodus file 1528 exodus stores each component of a vector-valued field as a separate variable. 1529 We assume that they are stored sequentially 1530 Zonal variables are accessed one element block at a time, so we loop through the cell sets, 1531 but once the vector has been reordered to natural size, we cannot use the label information 1532 to figure out what to save where. */ 1533 numCS = ex_inquire_int(exoid, EX_INQ_ELEM_BLK); 1534 PetscCall(PetscMalloc2(numCS, &csID, numCS, &csSize)); 1535 PetscCallExternal(ex_get_ids, exoid, EX_ELEM_BLOCK, csID); 1536 for (set = 0; set < numCS; ++set) { 1537 ex_block block; 1538 1539 block.id = csID[set]; 1540 block.type = EX_ELEM_BLOCK; 1541 PetscCallExternal(ex_get_block_param, exoid, &block); 1542 csSize[set] = block.num_entry; 1543 } 1544 PetscCall(VecGetOwnershipRange(vNatural, &xs, &xe)); 1545 PetscCall(VecGetBlockSize(vNatural, &bs)); 1546 if (bs > 1) PetscCall(ISCreateStride(comm, (xe - xs) / bs, xs, bs, &compIS)); 1547 for (set = 0; set < numCS; ++set) { 1548 PetscInt csLocalSize, c; 1549 1550 /* range of indices for set setID[set]: csxs:csxs + csSize[set]-1 1551 local slice of zonal values: xs/bs,xm/bs-1 1552 intersection: max(xs/bs,csxs),min(xm/bs-1,csxs + csSize[set]-1) */ 1553 csLocalSize = PetscMax(0, PetscMin(xe / bs, csxs + csSize[set]) - PetscMax(xs / bs, csxs)); 1554 if (bs == 1) { 1555 PetscCall(VecGetArray(vNatural, &varray)); 1556 PetscCallExternal(ex_get_partial_var, exoid, step, EX_ELEM_BLOCK, offset, csID[set], PetscMax(xs - csxs, 0) + 1, csLocalSize, &varray[PetscMax(0, csxs - xs)]); 1557 PetscCall(VecRestoreArray(vNatural, &varray)); 1558 } else { 1559 for (c = 0; c < bs; ++c) { 1560 PetscCall(ISStrideSetStride(compIS, (xe - xs) / bs, xs + c, bs)); 1561 PetscCall(VecGetSubVector(vNatural, compIS, &vComp)); 1562 PetscCall(VecGetArray(vComp, &varray)); 1563 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)]); 1564 PetscCall(VecRestoreArray(vComp, &varray)); 1565 PetscCall(VecRestoreSubVector(vNatural, compIS, &vComp)); 1566 } 1567 } 1568 csxs += csSize[set]; 1569 } 1570 PetscCall(PetscFree2(csID, csSize)); 1571 if (bs > 1) PetscCall(ISDestroy(&compIS)); 1572 if (useNatural) { 1573 PetscCall(DMPlexNaturalToGlobalBegin(dm, vNatural, v)); 1574 PetscCall(DMPlexNaturalToGlobalEnd(dm, vNatural, v)); 1575 PetscCall(VecDestroy(&vNatural)); 1576 } 1577 PetscFunctionReturn(PETSC_SUCCESS); 1578 } 1579 1580 /*@ 1581 PetscViewerExodusIIGetId - Get the file id of the `PETSCVIEWEREXODUSII` file 1582 1583 Logically Collective 1584 1585 Input Parameter: 1586 . viewer - the `PetscViewer` 1587 1588 Output Parameter: 1589 . exoid - The ExodusII file id 1590 1591 Level: intermediate 1592 1593 .seealso: `PETSCVIEWEREXODUSII`, `PetscViewer`, `PetscViewerFileSetMode()`, `PetscViewerCreate()`, `PetscViewerSetType()`, `PetscViewerBinaryOpen()` 1594 @*/ 1595 PetscErrorCode PetscViewerExodusIIGetId(PetscViewer viewer, PetscExodusIIInt *exoid) 1596 { 1597 PetscFunctionBegin; 1598 PetscValidHeaderSpecific(viewer, PETSC_VIEWER_CLASSID, 1); 1599 PetscTryMethod(viewer, "PetscViewerGetId_C", (PetscViewer, PetscExodusIIInt *), (viewer, exoid)); 1600 PetscFunctionReturn(PETSC_SUCCESS); 1601 } 1602 1603 /*@ 1604 PetscViewerExodusIISetOrder - Set the elements order in the exodusII file. 1605 1606 Collective 1607 1608 Input Parameters: 1609 + viewer - the `PETSCVIEWEREXODUSII` viewer 1610 - order - elements order 1611 1612 Output Parameter: 1613 1614 Level: beginner 1615 1616 .seealso: `PETSCVIEWEREXODUSII`, `PetscViewer`, `PetscViewerExodusIIGetId()`, `PetscViewerExodusIIGetOrder()` 1617 @*/ 1618 PetscErrorCode PetscViewerExodusIISetOrder(PetscViewer viewer, PetscInt order) 1619 { 1620 PetscFunctionBegin; 1621 PetscValidHeaderSpecific(viewer, PETSC_VIEWER_CLASSID, 1); 1622 PetscTryMethod(viewer, "PetscViewerSetOrder_C", (PetscViewer, PetscInt), (viewer, order)); 1623 PetscFunctionReturn(PETSC_SUCCESS); 1624 } 1625 1626 /*@ 1627 PetscViewerExodusIIGetOrder - Get the elements order in the exodusII file. 1628 1629 Collective 1630 1631 Input Parameters: 1632 + viewer - the `PETSCVIEWEREXODUSII` viewer 1633 - order - elements order 1634 1635 Output Parameter: 1636 1637 Level: beginner 1638 1639 .seealso: `PETSCVIEWEREXODUSII`, `PetscViewer`, `PetscViewerExodusIIGetId()`, `PetscViewerExodusIISetOrder()` 1640 @*/ 1641 PetscErrorCode PetscViewerExodusIIGetOrder(PetscViewer viewer, PetscInt *order) 1642 { 1643 PetscFunctionBegin; 1644 PetscValidHeaderSpecific(viewer, PETSC_VIEWER_CLASSID, 1); 1645 PetscTryMethod(viewer, "PetscViewerGetOrder_C", (PetscViewer, PetscInt *), (viewer, order)); 1646 PetscFunctionReturn(PETSC_SUCCESS); 1647 } 1648 1649 /*@ 1650 PetscViewerExodusIIOpen - Opens a file for ExodusII input/output. 1651 1652 Collective 1653 1654 Input Parameters: 1655 + comm - MPI communicator 1656 . name - name of file 1657 - mode - access mode 1658 .vb 1659 FILE_MODE_WRITE - create new file for binary output 1660 FILE_MODE_READ - open existing file for binary input 1661 FILE_MODE_APPEND - open existing file for binary output 1662 .ve 1663 1664 Output Parameter: 1665 . exo - `PETSCVIEWEREXODUSII` `PetscViewer` for Exodus II input/output to use with the specified file 1666 1667 Level: beginner 1668 1669 .seealso: `PETSCVIEWEREXODUSII`, `PetscViewer`, `PetscViewerPushFormat()`, `PetscViewerDestroy()`, 1670 `DMLoad()`, `PetscFileMode`, `PetscViewerSetType()`, `PetscViewerFileSetMode()`, `PetscViewerFileSetName()` 1671 @*/ 1672 PetscErrorCode PetscViewerExodusIIOpen(MPI_Comm comm, const char name[], PetscFileMode mode, PetscViewer *exo) 1673 { 1674 PetscFunctionBegin; 1675 PetscCall(PetscViewerCreate(comm, exo)); 1676 PetscCall(PetscViewerSetType(*exo, PETSCVIEWEREXODUSII)); 1677 PetscCall(PetscViewerFileSetMode(*exo, mode)); 1678 PetscCall(PetscViewerFileSetName(*exo, name)); 1679 PetscCall(PetscViewerSetFromOptions(*exo)); 1680 PetscFunctionReturn(PETSC_SUCCESS); 1681 } 1682 1683 static PetscErrorCode ExodusGetCellType_Internal(const char *elem_type, DMPolytopeType *ct) 1684 { 1685 PetscBool flg; 1686 1687 PetscFunctionBegin; 1688 *ct = DM_POLYTOPE_UNKNOWN; 1689 PetscCall(PetscStrcmp(elem_type, "BAR2", &flg)); 1690 if (flg) { 1691 *ct = DM_POLYTOPE_SEGMENT; 1692 goto done; 1693 } 1694 PetscCall(PetscStrcmp(elem_type, "BAR3", &flg)); 1695 if (flg) { 1696 *ct = DM_POLYTOPE_SEGMENT; 1697 goto done; 1698 } 1699 PetscCall(PetscStrcmp(elem_type, "TRI", &flg)); 1700 if (flg) { 1701 *ct = DM_POLYTOPE_TRIANGLE; 1702 goto done; 1703 } 1704 PetscCall(PetscStrcmp(elem_type, "TRI3", &flg)); 1705 if (flg) { 1706 *ct = DM_POLYTOPE_TRIANGLE; 1707 goto done; 1708 } 1709 PetscCall(PetscStrcmp(elem_type, "QUAD", &flg)); 1710 if (flg) { 1711 *ct = DM_POLYTOPE_QUADRILATERAL; 1712 goto done; 1713 } 1714 PetscCall(PetscStrcmp(elem_type, "QUAD4", &flg)); 1715 if (flg) { 1716 *ct = DM_POLYTOPE_QUADRILATERAL; 1717 goto done; 1718 } 1719 PetscCall(PetscStrcmp(elem_type, "SHELL4", &flg)); 1720 if (flg) { 1721 *ct = DM_POLYTOPE_QUADRILATERAL; 1722 goto done; 1723 } 1724 PetscCall(PetscStrcmp(elem_type, "TETRA", &flg)); 1725 if (flg) { 1726 *ct = DM_POLYTOPE_TETRAHEDRON; 1727 goto done; 1728 } 1729 PetscCall(PetscStrcmp(elem_type, "TET4", &flg)); 1730 if (flg) { 1731 *ct = DM_POLYTOPE_TETRAHEDRON; 1732 goto done; 1733 } 1734 PetscCall(PetscStrcmp(elem_type, "WEDGE", &flg)); 1735 if (flg) { 1736 *ct = DM_POLYTOPE_TRI_PRISM; 1737 goto done; 1738 } 1739 PetscCall(PetscStrcmp(elem_type, "HEX", &flg)); 1740 if (flg) { 1741 *ct = DM_POLYTOPE_HEXAHEDRON; 1742 goto done; 1743 } 1744 PetscCall(PetscStrcmp(elem_type, "HEX8", &flg)); 1745 if (flg) { 1746 *ct = DM_POLYTOPE_HEXAHEDRON; 1747 goto done; 1748 } 1749 PetscCall(PetscStrcmp(elem_type, "HEXAHEDRON", &flg)); 1750 if (flg) { 1751 *ct = DM_POLYTOPE_HEXAHEDRON; 1752 goto done; 1753 } 1754 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unrecognized element type %s", elem_type); 1755 done: 1756 PetscFunctionReturn(PETSC_SUCCESS); 1757 } 1758 1759 /*@ 1760 DMPlexCreateExodus - Create a `DMPLEX` mesh from an ExodusII file ID. 1761 1762 Collective 1763 1764 Input Parameters: 1765 + comm - The MPI communicator 1766 . exoid - The ExodusII id associated with a exodus file and obtained using ex_open 1767 - interpolate - Create faces and edges in the mesh 1768 1769 Output Parameter: 1770 . dm - The `DM` object representing the mesh 1771 1772 Level: beginner 1773 1774 .seealso: [](ch_unstructured), `DM`, `PETSCVIEWEREXODUSII`, `DMPLEX`, `DMCreate()` 1775 @*/ 1776 PetscErrorCode DMPlexCreateExodus(MPI_Comm comm, PetscExodusIIInt exoid, PetscBool interpolate, DM *dm) 1777 { 1778 PetscMPIInt num_proc, rank; 1779 DMLabel cellSets = NULL, faceSets = NULL, vertSets = NULL; 1780 PetscSection coordSection; 1781 Vec coordinates; 1782 PetscScalar *coords; 1783 PetscInt coordSize, v; 1784 /* Read from ex_get_init() */ 1785 char title[PETSC_MAX_PATH_LEN + 1]; 1786 int dim = 0, dimEmbed = 0, numVertices = 0, numCells = 0; 1787 int num_cs = 0, num_vs = 0, num_fs = 0; 1788 1789 PetscFunctionBegin; 1790 PetscCallMPI(MPI_Comm_rank(comm, &rank)); 1791 PetscCallMPI(MPI_Comm_size(comm, &num_proc)); 1792 PetscCall(DMCreate(comm, dm)); 1793 PetscCall(DMSetType(*dm, DMPLEX)); 1794 /* Open EXODUS II file and read basic information on rank 0, then broadcast to all processors */ 1795 if (rank == 0) { 1796 PetscCall(PetscMemzero(title, PETSC_MAX_PATH_LEN + 1)); 1797 PetscCallExternal(ex_get_init, exoid, title, &dimEmbed, &numVertices, &numCells, &num_cs, &num_vs, &num_fs); 1798 PetscCheck(num_cs, PETSC_COMM_SELF, PETSC_ERR_SUP, "Exodus file does not contain any cell set"); 1799 } 1800 PetscCallMPI(MPI_Bcast(title, PETSC_MAX_PATH_LEN + 1, MPI_CHAR, 0, comm)); 1801 PetscCallMPI(MPI_Bcast(&dim, 1, MPI_INT, 0, comm)); 1802 PetscCall(PetscObjectSetName((PetscObject)*dm, title)); 1803 PetscCall(DMPlexSetChart(*dm, 0, numCells + numVertices)); 1804 /* We do not want this label automatically computed, instead we compute it here */ 1805 PetscCall(DMCreateLabel(*dm, "celltype")); 1806 1807 /* Read cell sets information */ 1808 if (rank == 0) { 1809 PetscInt *cone; 1810 int c, cs, ncs, c_loc, v, v_loc; 1811 /* Read from ex_get_elem_blk_ids() */ 1812 int *cs_id, *cs_order; 1813 /* Read from ex_get_elem_block() */ 1814 char buffer[PETSC_MAX_PATH_LEN + 1]; 1815 int num_cell_in_set, num_vertex_per_cell, num_hybrid, num_attr; 1816 /* Read from ex_get_elem_conn() */ 1817 int *cs_connect; 1818 1819 /* Get cell sets IDs */ 1820 PetscCall(PetscMalloc2(num_cs, &cs_id, num_cs, &cs_order)); 1821 PetscCallExternal(ex_get_ids, exoid, EX_ELEM_BLOCK, cs_id); 1822 /* Read the cell set connectivity table and build mesh topology 1823 EXO standard requires that cells in cell sets be numbered sequentially and be pairwise disjoint. */ 1824 /* Check for a hybrid mesh */ 1825 for (cs = 0, num_hybrid = 0; cs < num_cs; ++cs) { 1826 DMPolytopeType ct; 1827 char elem_type[PETSC_MAX_PATH_LEN]; 1828 1829 PetscCall(PetscArrayzero(elem_type, sizeof(elem_type))); 1830 PetscCallExternal(ex_get_elem_type, exoid, cs_id[cs], elem_type); 1831 PetscCall(ExodusGetCellType_Internal(elem_type, &ct)); 1832 dim = PetscMax(dim, DMPolytopeTypeGetDim(ct)); 1833 PetscCallExternal(ex_get_block, exoid, EX_ELEM_BLOCK, cs_id[cs], buffer, &num_cell_in_set, &num_vertex_per_cell, 0, 0, &num_attr); 1834 switch (ct) { 1835 case DM_POLYTOPE_TRI_PRISM: 1836 cs_order[cs] = cs; 1837 ++num_hybrid; 1838 break; 1839 default: 1840 for (c = cs; c > cs - num_hybrid; --c) cs_order[c] = cs_order[c - 1]; 1841 cs_order[cs - num_hybrid] = cs; 1842 } 1843 } 1844 /* First set sizes */ 1845 for (ncs = 0, c = 0; ncs < num_cs; ++ncs) { 1846 DMPolytopeType ct; 1847 char elem_type[PETSC_MAX_PATH_LEN]; 1848 const PetscInt cs = cs_order[ncs]; 1849 1850 PetscCall(PetscArrayzero(elem_type, sizeof(elem_type))); 1851 PetscCallExternal(ex_get_elem_type, exoid, cs_id[cs], elem_type); 1852 PetscCall(ExodusGetCellType_Internal(elem_type, &ct)); 1853 PetscCallExternal(ex_get_block, exoid, EX_ELEM_BLOCK, cs_id[cs], buffer, &num_cell_in_set, &num_vertex_per_cell, 0, 0, &num_attr); 1854 for (c_loc = 0; c_loc < num_cell_in_set; ++c_loc, ++c) { 1855 PetscCall(DMPlexSetConeSize(*dm, c, num_vertex_per_cell)); 1856 PetscCall(DMPlexSetCellType(*dm, c, ct)); 1857 } 1858 } 1859 for (v = numCells; v < numCells + numVertices; ++v) PetscCall(DMPlexSetCellType(*dm, v, DM_POLYTOPE_POINT)); 1860 PetscCall(DMSetUp(*dm)); 1861 for (ncs = 0, c = 0; ncs < num_cs; ++ncs) { 1862 const PetscInt cs = cs_order[ncs]; 1863 PetscCallExternal(ex_get_block, exoid, EX_ELEM_BLOCK, cs_id[cs], buffer, &num_cell_in_set, &num_vertex_per_cell, 0, 0, &num_attr); 1864 PetscCall(PetscMalloc2(num_vertex_per_cell * num_cell_in_set, &cs_connect, num_vertex_per_cell, &cone)); 1865 PetscCallExternal(ex_get_conn, exoid, EX_ELEM_BLOCK, cs_id[cs], cs_connect, NULL, NULL); 1866 /* EXO uses Fortran-based indexing, DMPlex uses C-style and numbers cell first then vertices. */ 1867 for (c_loc = 0, v = 0; c_loc < num_cell_in_set; ++c_loc, ++c) { 1868 DMPolytopeType ct; 1869 1870 for (v_loc = 0; v_loc < num_vertex_per_cell; ++v_loc, ++v) cone[v_loc] = cs_connect[v] + numCells - 1; 1871 PetscCall(DMPlexGetCellType(*dm, c, &ct)); 1872 PetscCall(DMPlexInvertCell(ct, cone)); 1873 PetscCall(DMPlexSetCone(*dm, c, cone)); 1874 PetscCall(DMSetLabelValue_Fast(*dm, &cellSets, "Cell Sets", c, cs_id[cs])); 1875 } 1876 PetscCall(PetscFree2(cs_connect, cone)); 1877 } 1878 PetscCall(PetscFree2(cs_id, cs_order)); 1879 } 1880 { 1881 PetscInt ints[] = {dim, dimEmbed}; 1882 1883 PetscCallMPI(MPI_Bcast(ints, 2, MPIU_INT, 0, comm)); 1884 PetscCall(DMSetDimension(*dm, ints[0])); 1885 PetscCall(DMSetCoordinateDim(*dm, ints[1])); 1886 dim = ints[0]; 1887 dimEmbed = ints[1]; 1888 } 1889 PetscCall(DMPlexSymmetrize(*dm)); 1890 PetscCall(DMPlexStratify(*dm)); 1891 if (interpolate) { 1892 DM idm; 1893 1894 PetscCall(DMPlexInterpolate(*dm, &idm)); 1895 PetscCall(DMDestroy(dm)); 1896 *dm = idm; 1897 } 1898 1899 /* Create vertex set label */ 1900 if (rank == 0 && (num_vs > 0)) { 1901 int vs, v; 1902 /* Read from ex_get_node_set_ids() */ 1903 int *vs_id; 1904 /* Read from ex_get_node_set_param() */ 1905 int num_vertex_in_set; 1906 /* Read from ex_get_node_set() */ 1907 int *vs_vertex_list; 1908 1909 /* Get vertex set ids */ 1910 PetscCall(PetscMalloc1(num_vs, &vs_id)); 1911 PetscCallExternal(ex_get_ids, exoid, EX_NODE_SET, vs_id); 1912 for (vs = 0; vs < num_vs; ++vs) { 1913 PetscCallExternal(ex_get_set_param, exoid, EX_NODE_SET, vs_id[vs], &num_vertex_in_set, NULL); 1914 PetscCall(PetscMalloc1(num_vertex_in_set, &vs_vertex_list)); 1915 PetscCallExternal(ex_get_set, exoid, EX_NODE_SET, vs_id[vs], vs_vertex_list, NULL); 1916 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])); 1917 PetscCall(PetscFree(vs_vertex_list)); 1918 } 1919 PetscCall(PetscFree(vs_id)); 1920 } 1921 /* Read coordinates */ 1922 PetscCall(DMGetCoordinateSection(*dm, &coordSection)); 1923 PetscCall(PetscSectionSetNumFields(coordSection, 1)); 1924 PetscCall(PetscSectionSetFieldComponents(coordSection, 0, dimEmbed)); 1925 PetscCall(PetscSectionSetChart(coordSection, numCells, numCells + numVertices)); 1926 for (v = numCells; v < numCells + numVertices; ++v) { 1927 PetscCall(PetscSectionSetDof(coordSection, v, dimEmbed)); 1928 PetscCall(PetscSectionSetFieldDof(coordSection, v, 0, dimEmbed)); 1929 } 1930 PetscCall(PetscSectionSetUp(coordSection)); 1931 PetscCall(PetscSectionGetStorageSize(coordSection, &coordSize)); 1932 PetscCall(VecCreate(PETSC_COMM_SELF, &coordinates)); 1933 PetscCall(PetscObjectSetName((PetscObject)coordinates, "coordinates")); 1934 PetscCall(VecSetSizes(coordinates, coordSize, PETSC_DETERMINE)); 1935 PetscCall(VecSetBlockSize(coordinates, dimEmbed)); 1936 PetscCall(VecSetType(coordinates, VECSTANDARD)); 1937 PetscCall(VecGetArray(coordinates, &coords)); 1938 if (rank == 0) { 1939 PetscReal *x, *y, *z; 1940 1941 PetscCall(PetscMalloc3(numVertices, &x, numVertices, &y, numVertices, &z)); 1942 PetscCallExternal(ex_get_coord, exoid, x, y, z); 1943 if (dimEmbed > 0) { 1944 for (v = 0; v < numVertices; ++v) coords[v * dimEmbed + 0] = x[v]; 1945 } 1946 if (dimEmbed > 1) { 1947 for (v = 0; v < numVertices; ++v) coords[v * dimEmbed + 1] = y[v]; 1948 } 1949 if (dimEmbed > 2) { 1950 for (v = 0; v < numVertices; ++v) coords[v * dimEmbed + 2] = z[v]; 1951 } 1952 PetscCall(PetscFree3(x, y, z)); 1953 } 1954 PetscCall(VecRestoreArray(coordinates, &coords)); 1955 PetscCall(DMSetCoordinatesLocal(*dm, coordinates)); 1956 PetscCall(VecDestroy(&coordinates)); 1957 1958 /* Create side set label */ 1959 if (rank == 0 && interpolate && (num_fs > 0)) { 1960 int fs, f, voff; 1961 /* Read from ex_get_side_set_ids() */ 1962 int *fs_id; 1963 /* Read from ex_get_side_set_param() */ 1964 int num_side_in_set; 1965 /* Read from ex_get_side_set_node_list() */ 1966 int *fs_vertex_count_list, *fs_vertex_list, *fs_side_list; 1967 /* Read side set labels */ 1968 char fs_name[MAX_STR_LENGTH + 1]; 1969 size_t fs_name_len; 1970 1971 /* Get side set ids */ 1972 PetscCall(PetscMalloc1(num_fs, &fs_id)); 1973 PetscCallExternal(ex_get_ids, exoid, EX_SIDE_SET, fs_id); 1974 // Ids 1 and 2 are reserved by ExodusII for indicating things in 3D 1975 for (fs = 0; fs < num_fs; ++fs) { 1976 PetscCallExternal(ex_get_set_param, exoid, EX_SIDE_SET, fs_id[fs], &num_side_in_set, NULL); 1977 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)); 1978 PetscCallExternal(ex_get_side_set_node_list, exoid, fs_id[fs], fs_vertex_count_list, fs_vertex_list); 1979 PetscCallExternal(ex_get_set, exoid, EX_SIDE_SET, fs_id[fs], NULL, fs_side_list); 1980 1981 /* Get the specific name associated with this side set ID. */ 1982 int fs_name_err = ex_get_name(exoid, EX_SIDE_SET, fs_id[fs], fs_name); 1983 if (!fs_name_err) { 1984 PetscCall(PetscStrlen(fs_name, &fs_name_len)); 1985 if (fs_name_len == 0) PetscCall(PetscStrncpy(fs_name, "Face Sets", MAX_STR_LENGTH)); 1986 } 1987 for (f = 0, voff = 0; f < num_side_in_set; ++f) { 1988 const PetscInt *faces = NULL; 1989 PetscInt faceSize = fs_vertex_count_list[f], numFaces; 1990 PetscInt faceVertices[4], v; 1991 1992 PetscCheck(faceSize <= 4, comm, PETSC_ERR_ARG_WRONG, "ExodusII side cannot have %" PetscInt_FMT " > 4 vertices", faceSize); 1993 for (v = 0; v < faceSize; ++v, ++voff) faceVertices[v] = fs_vertex_list[voff] + numCells - 1; 1994 PetscCall(DMPlexGetFullJoin(*dm, faceSize, faceVertices, &numFaces, &faces)); 1995 PetscCheck(numFaces == 1, comm, PETSC_ERR_ARG_WRONG, "Invalid ExodusII side %d in set %d maps to %" PetscInt_FMT " faces", f, fs, numFaces); 1996 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]); 1997 PetscCall(DMSetLabelValue_Fast(*dm, &faceSets, "Face Sets", faces[0], fs_id[fs])); 1998 /* Only add the label if one has been detected for this side set. */ 1999 if (!fs_name_err) PetscCall(DMSetLabelValue(*dm, fs_name, faces[0], fs_id[fs])); 2000 PetscCall(DMPlexRestoreJoin(*dm, faceSize, faceVertices, &numFaces, &faces)); 2001 } 2002 PetscCall(PetscFree3(fs_vertex_count_list, fs_vertex_list, fs_side_list)); 2003 } 2004 PetscCall(PetscFree(fs_id)); 2005 } 2006 2007 { /* Create Cell/Face/Vertex Sets labels at all processes */ 2008 enum { 2009 n = 3 2010 }; 2011 PetscBool flag[n]; 2012 2013 flag[0] = cellSets ? PETSC_TRUE : PETSC_FALSE; 2014 flag[1] = faceSets ? PETSC_TRUE : PETSC_FALSE; 2015 flag[2] = vertSets ? PETSC_TRUE : PETSC_FALSE; 2016 PetscCallMPI(MPI_Bcast(flag, n, MPIU_BOOL, 0, comm)); 2017 if (flag[0]) PetscCall(DMCreateLabel(*dm, "Cell Sets")); 2018 if (flag[1]) PetscCall(DMCreateLabel(*dm, "Face Sets")); 2019 if (flag[2]) PetscCall(DMCreateLabel(*dm, "Vertex Sets")); 2020 } 2021 PetscFunctionReturn(PETSC_SUCCESS); 2022 } 2023