1 #define PETSCDM_DLL 2 #include <petsc/private/dmpleximpl.h> /*I "petscdmplex.h" I*/ 3 4 #include <netcdf.h> 5 #include <exodusII.h> 6 7 #include <petsc/private/viewerimpl.h> 8 #include <petsc/private/viewerexodusiiimpl.h> 9 /*@C 10 PETSC_VIEWER_EXODUSII_ - Creates an `PETSCVIEWEREXODUSII` `PetscViewer` shared by all processors in a communicator. 11 12 Collective; No Fortran Support 13 14 Input Parameter: 15 . comm - the MPI communicator to share the `PETSCVIEWEREXODUSII` `PetscViewer` 16 17 Level: intermediate 18 19 Note: 20 Unlike almost all other PETSc routines, `PETSC_VIEWER_EXODUSII_()` does not return 21 an error code. The GLVIS PetscViewer is usually used in the form 22 $ XXXView(XXX object, PETSC_VIEWER_EXODUSII_(comm)); 23 24 .seealso: `PETSCVIEWEREXODUSII`, `PetscViewer`, `PetscViewer`, `PetscViewerExodusIIOpen()`, `PetscViewerType`, `PetscViewerCreate()`, `PetscViewerDestroy()` 25 @*/ 26 PetscViewer PETSC_VIEWER_EXODUSII_(MPI_Comm comm) 27 { 28 PetscViewer viewer; 29 30 PetscFunctionBegin; 31 PetscCallNull(PetscViewerExodusIIOpen(comm, "mesh.exo", FILE_MODE_WRITE, &viewer)); 32 PetscCallNull(PetscObjectRegisterDestroy((PetscObject)viewer)); 33 PetscFunctionReturn(viewer); 34 } 35 36 static PetscErrorCode PetscViewerView_ExodusII(PetscViewer v, PetscViewer viewer) 37 { 38 PetscViewer_ExodusII *exo = (PetscViewer_ExodusII *)v->data; 39 40 PetscFunctionBegin; 41 if (exo->filename) PetscCall(PetscViewerASCIIPrintf(viewer, "Filename: %s\n", exo->filename)); 42 if (exo->exoid) PetscCall(PetscViewerASCIIPrintf(viewer, "exoid: %" PetscExodusIIInt_FMT "\n", exo->exoid)); 43 if (exo->btype) PetscCall(PetscViewerASCIIPrintf(viewer, "IO Mode: %d\n", exo->btype)); 44 if (exo->order) PetscCall(PetscViewerASCIIPrintf(viewer, "Mesh order: %" PetscInt_FMT "\n", exo->order)); 45 PetscCall(PetscViewerASCIIPrintf(viewer, "Number of nodal variables: %" PetscExodusIIInt_FMT "\n", exo->numNodalVariables)); 46 for (int i = 0; i < exo->numNodalVariables; i++) PetscCall(PetscViewerASCIIPrintf(viewer, " %d: %s\n", i, exo->nodalVariableNames[i])); 47 PetscCall(PetscViewerASCIIPrintf(viewer, "Number of zonal variables: %" PetscExodusIIInt_FMT "\n", exo->numZonalVariables)); 48 for (int i = 0; i < exo->numZonalVariables; i++) PetscCall(PetscViewerASCIIPrintf(viewer, " %d: %s\n", i, exo->zonalVariableNames[i])); 49 PetscFunctionReturn(PETSC_SUCCESS); 50 } 51 52 static PetscErrorCode PetscViewerFlush_ExodusII(PetscViewer v) 53 { 54 PetscViewer_ExodusII *exo = (PetscViewer_ExodusII *)v->data; 55 56 PetscFunctionBegin; 57 if (exo->exoid >= 0) PetscCallExternal(ex_update, exo->exoid); 58 PetscFunctionReturn(PETSC_SUCCESS); 59 } 60 61 static PetscErrorCode PetscViewerSetFromOptions_ExodusII(PetscViewer v, PetscOptionItems *PetscOptionsObject) 62 { 63 PetscFunctionBegin; 64 PetscOptionsHeadBegin(PetscOptionsObject, "ExodusII PetscViewer Options"); 65 PetscOptionsHeadEnd(); 66 PetscFunctionReturn(PETSC_SUCCESS); 67 } 68 69 static PetscErrorCode PetscViewerSetUp_ExodusII(PetscViewer viewer) 70 { 71 PetscFunctionBegin; 72 PetscFunctionReturn(PETSC_SUCCESS); 73 } 74 75 static PetscErrorCode PetscViewerDestroy_ExodusII(PetscViewer viewer) 76 { 77 PetscViewer_ExodusII *exo = (PetscViewer_ExodusII *)viewer->data; 78 79 PetscFunctionBegin; 80 if (exo->exoid >= 0) PetscCallExternal(ex_close, exo->exoid); 81 for (PetscInt i = 0; i < exo->numZonalVariables; i++) PetscCall(PetscFree(exo->zonalVariableNames[i])); 82 PetscCall(PetscFree(exo->zonalVariableNames)); 83 for (PetscInt i = 0; i < exo->numNodalVariables; i++) PetscCall(PetscFree(exo->nodalVariableNames[i])); 84 PetscCall(PetscFree(exo->nodalVariableNames)); 85 PetscCall(PetscFree(exo->filename)); 86 PetscCall(PetscFree(exo)); 87 PetscCall(PetscObjectComposeFunction((PetscObject)viewer, "PetscViewerFileSetName_C", NULL)); 88 PetscCall(PetscObjectComposeFunction((PetscObject)viewer, "PetscViewerFileGetName_C", NULL)); 89 PetscCall(PetscObjectComposeFunction((PetscObject)viewer, "PetscViewerFileSetMode_C", NULL)); 90 PetscCall(PetscObjectComposeFunction((PetscObject)viewer, "PetscViewerFileGetMode_C", NULL)); 91 PetscCall(PetscObjectComposeFunction((PetscObject)viewer, "PetscViewerGetId_C", NULL)); 92 PetscCall(PetscObjectComposeFunction((PetscObject)viewer, "PetscViewerGetOrder_C", NULL)); 93 PetscCall(PetscObjectComposeFunction((PetscObject)viewer, "PetscViewerSetOrder_C", NULL)); 94 PetscFunctionReturn(PETSC_SUCCESS); 95 } 96 97 static PetscErrorCode PetscViewerFileSetName_ExodusII(PetscViewer viewer, const char name[]) 98 { 99 PetscViewer_ExodusII *exo = (PetscViewer_ExodusII *)viewer->data; 100 PetscMPIInt rank; 101 PetscExodusIIInt CPU_word_size, IO_word_size, EXO_mode; 102 MPI_Info mpi_info = MPI_INFO_NULL; 103 PetscExodusIIFloat EXO_version; 104 105 PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)viewer), &rank)); 106 CPU_word_size = sizeof(PetscReal); 107 IO_word_size = sizeof(PetscReal); 108 109 PetscFunctionBegin; 110 if (exo->exoid >= 0) { 111 PetscCallExternal(ex_close, exo->exoid); 112 exo->exoid = -1; 113 } 114 if (exo->filename) PetscCall(PetscFree(exo->filename)); 115 PetscCall(PetscStrallocpy(name, &exo->filename)); 116 switch (exo->btype) { 117 case FILE_MODE_READ: 118 EXO_mode = EX_READ; 119 break; 120 case FILE_MODE_APPEND: 121 case FILE_MODE_UPDATE: 122 case FILE_MODE_APPEND_UPDATE: 123 /* Will fail if the file does not already exist */ 124 EXO_mode = EX_WRITE; 125 break; 126 case FILE_MODE_WRITE: 127 /* 128 exodus only allows writing geometry upon file creation, so we will let DMView create the file. 129 */ 130 PetscFunctionReturn(PETSC_SUCCESS); 131 break; 132 default: 133 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ORDER, "Must call PetscViewerFileSetMode() before PetscViewerFileSetName()"); 134 } 135 #if defined(PETSC_USE_64BIT_INDICES) 136 EXO_mode += EX_ALL_INT64_API; 137 #endif 138 exo->exoid = ex_open_par(name, EXO_mode, &CPU_word_size, &IO_word_size, &EXO_version, PETSC_COMM_WORLD, mpi_info); 139 PetscCheck(exo->exoid >= 0, PETSC_COMM_SELF, PETSC_ERR_LIB, "ex_open_par failed for %s", name); 140 PetscFunctionReturn(PETSC_SUCCESS); 141 } 142 143 static PetscErrorCode PetscViewerFileGetName_ExodusII(PetscViewer viewer, const char *name[]) 144 { 145 PetscViewer_ExodusII *exo = (PetscViewer_ExodusII *)viewer->data; 146 147 PetscFunctionBegin; 148 *name = exo->filename; 149 PetscFunctionReturn(PETSC_SUCCESS); 150 } 151 152 static PetscErrorCode PetscViewerFileSetMode_ExodusII(PetscViewer viewer, PetscFileMode type) 153 { 154 PetscViewer_ExodusII *exo = (PetscViewer_ExodusII *)viewer->data; 155 156 PetscFunctionBegin; 157 exo->btype = type; 158 PetscFunctionReturn(PETSC_SUCCESS); 159 } 160 161 static PetscErrorCode PetscViewerFileGetMode_ExodusII(PetscViewer viewer, PetscFileMode *type) 162 { 163 PetscViewer_ExodusII *exo = (PetscViewer_ExodusII *)viewer->data; 164 165 PetscFunctionBegin; 166 *type = exo->btype; 167 PetscFunctionReturn(PETSC_SUCCESS); 168 } 169 170 static PetscErrorCode PetscViewerExodusIIGetId_ExodusII(PetscViewer viewer, PetscExodusIIInt *exoid) 171 { 172 PetscViewer_ExodusII *exo = (PetscViewer_ExodusII *)viewer->data; 173 174 PetscFunctionBegin; 175 *exoid = exo->exoid; 176 PetscFunctionReturn(PETSC_SUCCESS); 177 } 178 179 static PetscErrorCode PetscViewerExodusIIGetOrder_ExodusII(PetscViewer viewer, PetscInt *order) 180 { 181 PetscViewer_ExodusII *exo = (PetscViewer_ExodusII *)viewer->data; 182 183 PetscFunctionBegin; 184 *order = exo->order; 185 PetscFunctionReturn(PETSC_SUCCESS); 186 } 187 188 static PetscErrorCode PetscViewerExodusIISetOrder_ExodusII(PetscViewer viewer, PetscInt order) 189 { 190 PetscViewer_ExodusII *exo = (PetscViewer_ExodusII *)viewer->data; 191 192 PetscFunctionBegin; 193 exo->order = order; 194 PetscFunctionReturn(PETSC_SUCCESS); 195 } 196 197 /*@ 198 PetscViewerExodusIISetZonalVariable - Sets the number of zonal variables in an exodusII file 199 200 Collective; 201 202 Input Parameters: 203 + viewer - a `PetscViewer` of type `PETSCVIEWEREXODUSII` 204 - num - the number of zonal variables in the exodusII file 205 206 Level: intermediate 207 208 Notes: 209 The exodusII API does not allow changing the number of variables in a file so this function will return an error 210 if called twice, called on a read-only file, or called on file for which the number of variables has already been specified 211 212 .seealso: `PETSCVIEWEREXODUSII`, `PetscViewer`, `PetscViewerCreate()`, `PetscViewerDestroy()`, `PetscViewerExodusIIOpen()`, `PetscViewerSetType()`, `PetscViewerType`, `PetscViewerExodusIIGetZonalVariable()` 213 @*/ 214 PetscErrorCode PetscViewerExodusIISetZonalVariable(PetscViewer viewer, PetscExodusIIInt num) 215 { 216 PetscViewer_ExodusII *exo = (PetscViewer_ExodusII *)viewer->data; 217 MPI_Comm comm; 218 PetscExodusIIInt exoid = -1; 219 220 PetscFunctionBegin; 221 PetscCall(PetscObjectGetComm((PetscObject)viewer, &comm)); 222 PetscCheck(exo->numZonalVariables == -1, comm, PETSC_ERR_SUP, "The number of zonal variables has already been set to %d and cannot be overwritten", exo->numZonalVariables); 223 PetscCheck((exo->btype != FILE_MODE_READ) && (exo->btype != FILE_MODE_UNDEFINED), comm, PETSC_ERR_FILE_WRITE, "Cannot set the number of variables because the file is not writable"); 224 225 exo->numZonalVariables = num; 226 PetscCall(PetscMalloc1(num, &exo->zonalVariableNames)); 227 for (int i = 0; i < num; i++) { exo->zonalVariableNames[i] = NULL; } 228 PetscCall(PetscViewerExodusIIGetId(viewer, &exoid)); 229 PetscCallExternal(ex_put_variable_param, exoid, EX_ELEM_BLOCK, num); 230 PetscFunctionReturn(PETSC_SUCCESS); 231 } 232 233 /*@ 234 PetscViewerExodusIISetNodalVariable - Sets the number of nodal variables in an exodusII file 235 236 Collective; 237 238 Input Parameters: 239 + viewer - a `PetscViewer` of type `PETSCVIEWEREXODUSII` 240 - num - the number of nodal variables in the exodusII file 241 242 Level: intermediate 243 244 Notes: 245 The exodusII API does not allow changing the number of variables in a file so this function will return an error 246 if called twice, called on a read-only file, or called on file for which the number of variables has already been specified 247 248 .seealso: `PETSCVIEWEREXODUSII`, `PetscViewer`, `PetscViewerCreate()`, `PetscViewerDestroy()`, `PetscViewerExodusIIOpen()`, `PetscViewerSetType()`, `PetscViewerType`, `PetscViewerExodusIIGetNodalVariable()` 249 @*/ 250 PetscErrorCode PetscViewerExodusIISetNodalVariable(PetscViewer viewer, PetscExodusIIInt num) 251 { 252 PetscViewer_ExodusII *exo = (PetscViewer_ExodusII *)viewer->data; 253 MPI_Comm comm; 254 PetscExodusIIInt exoid = -1; 255 256 PetscFunctionBegin; 257 PetscCall(PetscObjectGetComm((PetscObject)viewer, &comm)); 258 PetscCheck(exo->numNodalVariables == -1, comm, PETSC_ERR_SUP, "The number of nodal variables has already been set to %d and cannot be overwritten", exo->numNodalVariables); 259 PetscCheck((exo->btype != FILE_MODE_READ) && (exo->btype != FILE_MODE_UNDEFINED), comm, PETSC_ERR_FILE_WRITE, "Cannot set the number of variables because the file is not writable"); 260 261 exo->numNodalVariables = num; 262 PetscCall(PetscMalloc1(num, &exo->nodalVariableNames)); 263 for (int i = 0; i < num; i++) { exo->nodalVariableNames[i] = NULL; } 264 PetscCall(PetscViewerExodusIIGetId(viewer, &exoid)); 265 PetscCallExternal(ex_put_variable_param, exoid, EX_NODAL, num); 266 PetscFunctionReturn(PETSC_SUCCESS); 267 } 268 269 /*@ 270 PetscViewerExodusIIGetZonalVariable - Gets the number of zonal variables in an exodusII file 271 272 Collective 273 274 Input Parameters: 275 . viewer - a `PetscViewer` of type `PETSCVIEWEREXODUSII` 276 277 Output Parameters: 278 . num - the number variables in the exodusII file 279 280 Level: intermediate 281 282 Notes: 283 The number of variables in the exodusII file is cached in the viewer 284 285 .seealso: `PETSCVIEWEREXODUSII`, `PetscViewer`, `PetscViewerCreate()`, `PetscViewerDestroy()`, `PetscViewerExodusIIOpen()`, `PetscViewerSetType()`, `PetscViewerType`, `PetscViewerExodusIIsetZonalVariable()` 286 @*/ 287 PetscErrorCode PetscViewerExodusIIGetZonalVariable(PetscViewer viewer, PetscExodusIIInt *num) 288 { 289 PetscViewer_ExodusII *exo = (PetscViewer_ExodusII *)viewer->data; 290 MPI_Comm comm; 291 PetscExodusIIInt exoid = -1; 292 293 PetscFunctionBegin; 294 if (exo->numZonalVariables > -1) { 295 *num = exo->numZonalVariables; 296 } else { 297 PetscCall(PetscViewerExodusIIGetId(viewer, &exoid)); 298 PetscCall(PetscObjectGetComm((PetscObject)viewer, &comm)); 299 PetscCheck(exoid > 0, comm, PETSC_ERR_FILE_OPEN, "Exodus file is not open"); 300 PetscCallExternal(ex_get_variable_param, exoid, EX_ELEM_BLOCK, num); 301 exo->numZonalVariables = *num; 302 PetscCall(PetscMalloc1(*num, &exo->zonalVariableNames)); 303 for (int i = 0; i < *num; i++) { exo->zonalVariableNames[i] = NULL; } 304 } 305 PetscFunctionReturn(PETSC_SUCCESS); 306 } 307 308 /*@ 309 PetscViewerExodusIIGetNodalVariable - Gets the number of nodal variables in an exodusII file 310 311 Collective 312 313 Input Parameters: 314 . viewer - a `PetscViewer` of type `PETSCVIEWEREXODUSII` 315 316 Output Parameters: 317 . num - the number variables in the exodusII file 318 319 Level: intermediate 320 321 Notes: 322 This function gets the number of nodal variables and saves it in the address of num. 323 324 .seealso: `PETSCVIEWEREXODUSII`, `PetscViewer`, `PetscViewerCreate()`, `PetscViewerDestroy()`, `PetscViewerExodusIIOpen()`, `PetscViewerSetType()`, `PetscViewerType`, `PetscViewerExodusIISetNodalVariable()` 325 @*/ 326 PetscErrorCode PetscViewerExodusIIGetNodalVariable(PetscViewer viewer, PetscExodusIIInt *num) 327 { 328 PetscViewer_ExodusII *exo = (PetscViewer_ExodusII *)viewer->data; 329 MPI_Comm comm; 330 PetscExodusIIInt exoid = -1; 331 332 PetscFunctionBegin; 333 if (exo->numNodalVariables > -1) { 334 *num = exo->numNodalVariables; 335 } else { 336 PetscCall(PetscViewerExodusIIGetId(viewer, &exoid)); 337 PetscCall(PetscObjectGetComm((PetscObject)viewer, &comm)); 338 PetscCheck(exoid > 0, comm, PETSC_ERR_FILE_OPEN, "Exodus file is not open"); 339 PetscCallExternal(ex_get_variable_param, exoid, EX_NODAL, num); 340 exo->numNodalVariables = *num; 341 PetscCall(PetscMalloc1(*num, &exo->nodalVariableNames)); 342 for (int i = 0; i < *num; i++) { exo->nodalVariableNames[i] = NULL; } 343 } 344 PetscFunctionReturn(PETSC_SUCCESS); 345 } 346 347 /*@ 348 PetscViewerExodusIISetZonalVariableName - Sets the name of a zonal variable. 349 350 Collective; 351 352 Input Parameters: 353 + viewer - a `PetscViewer` of type `PETSCVIEWEREXODUSII` 354 . idx - the index for which you want to save the name 355 - name - string containing the name characters 356 357 Level: intermediate 358 359 .seealso: `PETSCVIEWEREXODUSII`, `PetscViewer`, `PetscViewerCreate()`, `PetscViewerDestroy()`, `PetscViewerExodusIIOpen()`, `PetscViewerSetType()`, `PetscViewerType`, `PetscViewerExodusIIGetZonalVariableName()` 360 @*/ 361 PetscErrorCode PetscViewerExodusIISetZonalVariableName(PetscViewer viewer, PetscExodusIIInt idx, const char name[]) 362 { 363 PetscViewer_ExodusII *exo = (PetscViewer_ExodusII *)viewer->data; 364 365 PetscFunctionBegin; 366 PetscCheck((idx >= 0) && (idx < exo->numZonalVariables), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Variable index out of range. Was PetscViewerExodusIISetZonalVariable called?"); 367 PetscCall(PetscStrallocpy(name, (char **)&exo->zonalVariableNames[idx])); 368 PetscCallExternal(ex_put_variable_name, exo->exoid, EX_ELEM_BLOCK, idx + 1, name); 369 PetscFunctionReturn(PETSC_SUCCESS); 370 } 371 372 /*@ 373 PetscViewerExodusIISetNodalVariableName - Sets the name of a nodal variable. 374 375 Collective; 376 377 Input Parameters: 378 + viewer - a `PetscViewer` of type `PETSCVIEWEREXODUSII` 379 . idx - the index for which you want to save the name 380 - name - string containing the name characters 381 382 Level: intermediate 383 384 .seealso: `PETSCVIEWEREXODUSII`, `PetscViewer`, `PetscViewerCreate()`, `PetscViewerDestroy()`, `PetscViewerExodusIIOpen()`, `PetscViewerSetType()`, `PetscViewerType`, `PetscViewerExodusIIGetNodalVariableName()` 385 @*/ 386 PetscErrorCode PetscViewerExodusIISetNodalVariableName(PetscViewer viewer, PetscExodusIIInt idx, const char name[]) 387 { 388 PetscViewer_ExodusII *exo = (PetscViewer_ExodusII *)viewer->data; 389 390 PetscFunctionBegin; 391 PetscCheck((idx >= 0) && (idx < exo->numNodalVariables), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Variable index out of range. Was PetscViewerExodusIISetNodalVariable called?"); 392 PetscCall(PetscStrallocpy(name, (char **)&exo->nodalVariableNames[idx])); 393 PetscCallExternal(ex_put_variable_name, exo->exoid, EX_NODAL, idx + 1, name); 394 PetscFunctionReturn(PETSC_SUCCESS); 395 } 396 397 /*@ 398 PetscViewerExodusIIGetZonalVariableName - Gets the name of a zonal variable. 399 400 Collective; 401 402 Input Parameters: 403 + viewer - a `PetscViewer` of type `PETSCVIEWEREXODUSII` 404 - idx - the index for which you want to get the name 405 406 Output Parameter: 407 . name - pointer to the string containing the name characters 408 409 Level: intermediate 410 411 .seealso: `PETSCVIEWEREXODUSII`, `PetscViewer`, `PetscViewerCreate()`, `PetscViewerDestroy()`, `PetscViewerExodusIIOpen()`, `PetscViewerSetType()`, `PetscViewerType`, `PetscViewerExodusIISetZonalVariableName()` 412 @*/ 413 PetscErrorCode PetscViewerExodusIIGetZonalVariableName(PetscViewer viewer, PetscExodusIIInt idx, const char *name[]) 414 { 415 PetscViewer_ExodusII *exo = (PetscViewer_ExodusII *)viewer->data; 416 PetscExodusIIInt exoid = -1; 417 char tmpName[MAX_NAME_LENGTH + 1]; 418 419 PetscFunctionBegin; 420 PetscCheck(idx >= 0 && idx < exo->numZonalVariables, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Variable index out of range. Was PetscViewerExodusIISetZonalVariable called?"); 421 if (!exo->zonalVariableNames[idx]) { 422 PetscCall(PetscViewerExodusIIGetId(viewer, &exoid)); 423 PetscCallExternal(ex_get_variable_name, exoid, EX_ELEM_BLOCK, idx + 1, tmpName); 424 PetscCall(PetscStrallocpy(tmpName, (char **)&exo->zonalVariableNames[idx])); 425 } 426 *name = exo->zonalVariableNames[idx]; 427 PetscFunctionReturn(PETSC_SUCCESS); 428 } 429 430 /*@ 431 PetscViewerExodusIIGetNodalVariableName - Gets the name of a nodal variable. 432 433 Collective; 434 435 Input Parameters: 436 + viewer - a `PetscViewer` of type `PETSCVIEWEREXODUSII` 437 - idx - the index for which you want to save the name 438 439 Output Parameter: 440 . name - string array containing name characters 441 442 Level: intermediate 443 444 .seealso: `PETSCVIEWEREXODUSII`, `PetscViewer`, `PetscViewerCreate()`, `PetscViewerDestroy()`, `PetscViewerExodusIIOpen()`, `PetscViewerSetType()`, `PetscViewerType`, `PetscViewerExodusIISetNodalVariableName()` 445 @*/ 446 PetscErrorCode PetscViewerExodusIIGetNodalVariableName(PetscViewer viewer, PetscExodusIIInt idx, const char *name[]) 447 { 448 PetscViewer_ExodusII *exo = (PetscViewer_ExodusII *)viewer->data; 449 PetscExodusIIInt exoid = -1; 450 char tmpName[MAX_NAME_LENGTH + 1]; 451 452 PetscFunctionBegin; 453 PetscCheck((idx >= 0) && (idx < exo->numNodalVariables), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Variable index out of range. Was PetscViewerExodusIISetNodalVariable called?"); 454 if (!exo->nodalVariableNames[idx]) { 455 PetscCall(PetscViewerExodusIIGetId(viewer, &exoid)); 456 PetscCallExternal(ex_get_variable_name, exoid, EX_NODAL, idx + 1, tmpName); 457 PetscCall(PetscStrallocpy(tmpName, (char **)&exo->nodalVariableNames[idx])); 458 } 459 *name = exo->nodalVariableNames[idx]; 460 PetscFunctionReturn(PETSC_SUCCESS); 461 } 462 463 /*@C 464 PetscViewerExodusIISetZonalVariableNames - Sets the names of all nodal variables 465 466 Collective; No Fortran Support 467 468 Input Parameters: 469 + viewer - a `PetscViewer` of type `PETSCVIEWEREXODUSII` 470 - names - 2D array containing the array of string names to be set 471 472 Level: intermediate 473 474 Notes: 475 This function allows users to set multiple zonal variable names at a time. 476 477 .seealso: `PETSCVIEWEREXODUSII`, `PetscViewer`, `PetscViewerCreate()`, `PetscViewerDestroy()`, `PetscViewerExodusIIOpen()`, `PetscViewerSetType()`, `PetscViewerType`, `PetscViewerExodusIIGetZonalVariableNames()` 478 @*/ 479 PetscErrorCode PetscViewerExodusIISetZonalVariableNames(PetscViewer viewer, const char *names[]) 480 { 481 PetscExodusIIInt numNames; 482 PetscExodusIIInt exoid = -1; 483 PetscViewer_ExodusII *exo = (PetscViewer_ExodusII *)viewer->data; 484 485 PetscFunctionBegin; 486 PetscCall(PetscViewerExodusIIGetZonalVariable(viewer, &numNames)); 487 PetscCheck(numNames >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Number of zonal variables not set. Was PetscViewerExodusIISetZonalVariable called?"); 488 489 PetscCall(PetscViewerExodusIIGetId(viewer, &exoid)); 490 for (int i = 0; i < numNames; i++) { 491 PetscCall(PetscStrallocpy(names[i], &exo->zonalVariableNames[i])); 492 PetscCallExternal(ex_put_variable_name, exoid, EX_ELEM_BLOCK, i + 1, names[i]); 493 } 494 PetscFunctionReturn(PETSC_SUCCESS); 495 } 496 497 /*@C 498 PetscViewerExodusIISetNodalVariableNames - Sets the names of all nodal variables. 499 500 Collective; No Fortran Support 501 502 Input Parameters: 503 + viewer - a `PetscViewer` of type `PETSCVIEWEREXODUSII` 504 - names - 2D array containing the array of string names to be set 505 506 Level: intermediate 507 508 Notes: 509 This function allows users to set multiple nodal variable names at a time. 510 511 .seealso: `PETSCVIEWEREXODUSII`, `PetscViewer`, `PetscViewerCreate()`, `PetscViewerDestroy()`, `PetscViewerExodusIIOpen()`, `PetscViewerSetType()`, `PetscViewerType`, `PetscViewerExodusIIGetNodalVariableNames()` 512 @*/ 513 PetscErrorCode PetscViewerExodusIISetNodalVariableNames(PetscViewer viewer, const char *names[]) 514 { 515 PetscExodusIIInt numNames; 516 PetscExodusIIInt exoid = -1; 517 PetscViewer_ExodusII *exo = (PetscViewer_ExodusII *)viewer->data; 518 519 PetscFunctionBegin; 520 PetscCall(PetscViewerExodusIIGetNodalVariable(viewer, &numNames)); 521 PetscCheck(numNames >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Number of nodal variables not set. Was PetscViewerExodusIISetNodalVariable called?"); 522 523 PetscCall(PetscViewerExodusIIGetId(viewer, &exoid)); 524 for (int i = 0; i < numNames; i++) { 525 PetscCall(PetscStrallocpy(names[i], &exo->nodalVariableNames[i])); 526 PetscCallExternal(ex_put_variable_name, exoid, EX_NODAL, i + 1, names[i]); 527 } 528 PetscFunctionReturn(PETSC_SUCCESS); 529 } 530 531 /*@C 532 PetscViewerExodusIIGetZonalVariableNames - Gets the names of all zonal variables. 533 534 Collective; No Fortran Support 535 536 Input Parameters: 537 + viewer - a `PetscViewer` of type `PETSCVIEWEREXODUSII` 538 - numVars - the number of zonal variable names to retrieve 539 540 Output Parameters: 541 . varNames - pointer to a 2D array where the zonal variable names will be saved 542 543 Level: intermediate 544 545 Notes: 546 This function returns a borrowed pointer which should not be freed. 547 548 .seealso: `PETSCVIEWEREXODUSII`, `PetscViewer`, `PetscViewerCreate()`, `PetscViewerDestroy()`, `PetscViewerExodusIIOpen()`, `PetscViewerSetType()`, `PetscViewerType`, `PetscViewerExodusIISetZonalVariableNames()` 549 @*/ 550 PetscErrorCode PetscViewerExodusIIGetZonalVariableNames(PetscViewer viewer, PetscExodusIIInt *numVars, char ***varNames) 551 { 552 PetscViewer_ExodusII *exo = (PetscViewer_ExodusII *)viewer->data; 553 PetscExodusIIInt idx; 554 char tmpName[MAX_NAME_LENGTH + 1]; 555 PetscExodusIIInt exoid = -1; 556 557 PetscFunctionBegin; 558 PetscCall(PetscViewerExodusIIGetZonalVariable(viewer, numVars)); 559 /* 560 Cache variable names if necessary 561 */ 562 for (idx = 0; idx < *numVars; idx++) { 563 if (!exo->zonalVariableNames[idx]) { 564 PetscCall(PetscViewerExodusIIGetId(viewer, &exoid)); 565 PetscCallExternal(ex_get_variable_name, exoid, EX_ELEM_BLOCK, idx + 1, tmpName); 566 PetscCall(PetscStrallocpy(tmpName, (char **)&exo->zonalVariableNames[idx])); 567 } 568 } 569 *varNames = (char **)exo->zonalVariableNames; 570 PetscFunctionReturn(PETSC_SUCCESS); 571 } 572 573 /*@C 574 PetscViewerExodusIIGetNodalVariableNames - Gets the names of all nodal variables. 575 576 Collective; No Fortran Support 577 578 Input Parameters: 579 + viewer - a `PetscViewer` of type `PETSCVIEWEREXODUSII` 580 - numVars - the number of nodal variable names to retrieve 581 582 Output Parameters: 583 . varNames - 2D array where the nodal variable names will be saved 584 585 Level: intermediate 586 587 Notes: 588 This function returns a borrowed pointer which should not be freed. 589 590 .seealso: `PETSCVIEWEREXODUSII`, `PetscViewer`, `PetscViewerCreate()`, `PetscViewerDestroy()`, `PetscViewerExodusIIOpen()`, `PetscViewerSetType()`, `PetscViewerType`, `PetscViewerExodusIISetNodalVariableNames()` 591 @*/ 592 PetscErrorCode PetscViewerExodusIIGetNodalVariableNames(PetscViewer viewer, PetscExodusIIInt *numVars, char ***varNames) 593 { 594 PetscViewer_ExodusII *exo = (PetscViewer_ExodusII *)viewer->data; 595 PetscExodusIIInt idx; 596 char tmpName[MAX_NAME_LENGTH + 1]; 597 PetscExodusIIInt exoid = -1; 598 599 PetscFunctionBegin; 600 PetscCall(PetscViewerExodusIIGetNodalVariable(viewer, numVars)); 601 /* 602 Cache variable names if necessary 603 */ 604 for (idx = 0; idx < *numVars; idx++) { 605 if (!exo->nodalVariableNames[idx]) { 606 PetscCall(PetscViewerExodusIIGetId(viewer, &exoid)); 607 PetscCallExternal(ex_get_variable_name, exoid, EX_NODAL, idx + 1, tmpName); 608 PetscCall(PetscStrallocpy(tmpName, (char **)&exo->nodalVariableNames[idx])); 609 } 610 } 611 *varNames = (char **)exo->nodalVariableNames; 612 PetscFunctionReturn(PETSC_SUCCESS); 613 } 614 615 /*MC 616 PETSCVIEWEREXODUSII - A viewer that writes to an Exodus II file 617 618 Level: beginner 619 620 .seealso: `PetscViewerExodusIIOpen()`, `PetscViewerCreate()`, `PETSCVIEWERBINARY`, `PETSCVIEWERHDF5`, `DMView()`, 621 `PetscViewerFileSetName()`, `PetscViewerFileSetMode()`, `PetscViewerFormat`, `PetscViewerType`, `PetscViewerSetType()` 622 M*/ 623 PETSC_EXTERN PetscErrorCode PetscViewerCreate_ExodusII(PetscViewer v) 624 { 625 PetscViewer_ExodusII *exo; 626 627 PetscFunctionBegin; 628 PetscCall(PetscNew(&exo)); 629 630 v->data = (void *)exo; 631 v->ops->destroy = PetscViewerDestroy_ExodusII; 632 v->ops->setfromoptions = PetscViewerSetFromOptions_ExodusII; 633 v->ops->setup = PetscViewerSetUp_ExodusII; 634 v->ops->view = PetscViewerView_ExodusII; 635 v->ops->flush = PetscViewerFlush_ExodusII; 636 exo->btype = FILE_MODE_UNDEFINED; 637 exo->filename = 0; 638 exo->exoid = -1; 639 exo->numNodalVariables = -1; 640 exo->numZonalVariables = -1; 641 exo->nodalVariableNames = NULL; 642 exo->zonalVariableNames = NULL; 643 644 PetscCall(PetscObjectComposeFunction((PetscObject)v, "PetscViewerFileSetName_C", PetscViewerFileSetName_ExodusII)); 645 PetscCall(PetscObjectComposeFunction((PetscObject)v, "PetscViewerFileGetName_C", PetscViewerFileGetName_ExodusII)); 646 PetscCall(PetscObjectComposeFunction((PetscObject)v, "PetscViewerFileSetMode_C", PetscViewerFileSetMode_ExodusII)); 647 PetscCall(PetscObjectComposeFunction((PetscObject)v, "PetscViewerFileGetMode_C", PetscViewerFileGetMode_ExodusII)); 648 PetscCall(PetscObjectComposeFunction((PetscObject)v, "PetscViewerGetId_C", PetscViewerExodusIIGetId_ExodusII)); 649 PetscCall(PetscObjectComposeFunction((PetscObject)v, "PetscViewerSetOrder_C", PetscViewerExodusIISetOrder_ExodusII)); 650 PetscCall(PetscObjectComposeFunction((PetscObject)v, "PetscViewerGetOrder_C", PetscViewerExodusIIGetOrder_ExodusII)); 651 PetscFunctionReturn(PETSC_SUCCESS); 652 } 653 654 /*@ 655 PetscViewerExodusIIGetNodalVariableIndex - return the location of a nodal variable in an exodusII file given its name 656 657 Collective 658 659 Input Parameters: 660 + viewer - a `PetscViewer` of type `PETSCVIEWEREXODUSII` 661 - name - the name of the result 662 663 Output Parameter: 664 . varIndex - the location of the variable in the exodus file or -1 if the variable is not found 665 666 Level: beginner 667 668 Notes: 669 The exodus variable index is obtained by comparing the name argument to the 670 names of zonal variables declared in the exodus file. For instance if name is "V" 671 the location in the exodus file will be the first match of "V", "V_X", "V_XX", "V_1", or "V_11" 672 amongst all variables of type obj_type. 673 674 .seealso: `PetscViewerExodusIISetNodalVariable()`, `PetscViewerExodusIIGetNodalVariable()`, `PetscViewerExodusIISetNodalVariableName()`, `PetscViewerExodusIISetNodalVariableNames()`, `PetscViewerExodusIIGetNodalVariableName()`, `PetscViewerExodusIIGetNodalVariableNames()` 675 @*/ 676 PetscErrorCode PetscViewerExodusIIGetNodalVariableIndex(PetscViewer viewer, const char name[], PetscExodusIIInt *varIndex) 677 { 678 PetscExodusIIInt num_vars = 0, i, j; 679 char ext_name[MAX_STR_LENGTH + 1]; 680 char **var_names; 681 const int num_suffix = 5; 682 char *suffix[5]; 683 PetscBool flg; 684 685 PetscFunctionBegin; 686 suffix[0] = (char *)""; 687 suffix[1] = (char *)"_X"; 688 suffix[2] = (char *)"_XX"; 689 suffix[3] = (char *)"_1"; 690 suffix[4] = (char *)"_11"; 691 *varIndex = -1; 692 693 PetscCall(PetscViewerExodusIIGetNodalVariableNames(viewer, &num_vars, &var_names)); 694 for (i = 0; i < num_vars; ++i) { 695 for (j = 0; j < num_suffix; ++j) { 696 PetscCall(PetscStrncpy(ext_name, name, MAX_STR_LENGTH)); 697 PetscCall(PetscStrlcat(ext_name, suffix[j], MAX_STR_LENGTH)); 698 PetscCall(PetscStrcasecmp(ext_name, var_names[i], &flg)); 699 if (flg) *varIndex = i; 700 } 701 if (flg) break; 702 } 703 PetscFunctionReturn(PETSC_SUCCESS); 704 } 705 706 /*@ 707 PetscViewerExodusIIGetZonalVariableIndex - return the location of a zonal variable in an exodusII file given its name 708 709 Collective 710 711 Input Parameters: 712 + viewer - a `PetscViewer` of type `PETSCVIEWEREXODUSII` 713 - name - the name of the result 714 715 Output Parameter: 716 . varIndex - the location of the variable in the exodus file or -1 if the variable is not found 717 718 Level: beginner 719 720 Notes: 721 The exodus variable index is obtained by comparing the name argument to the 722 names of zonal variables declared in the exodus file. For instance if name is "V" 723 the location in the exodus file will be the first match of "V", "V_X", "V_XX", "V_1", or "V_11" 724 amongst all variables of type obj_type. 725 726 .seealso: `PetscViewerExodusIISetNodalVariable()`, `PetscViewerExodusIIGetNodalVariable()`, `PetscViewerExodusIISetNodalVariableName()`, `PetscViewerExodusIISetNodalVariableNames()`, `PetscViewerExodusIIGetNodalVariableName()`, `PetscViewerExodusIIGetNodalVariableNames()` 727 @*/ 728 PetscErrorCode PetscViewerExodusIIGetZonalVariableIndex(PetscViewer viewer, const char name[], int *varIndex) 729 { 730 PetscExodusIIInt num_vars = 0, i, j; 731 char ext_name[MAX_STR_LENGTH + 1]; 732 char **var_names; 733 const int num_suffix = 5; 734 char *suffix[5]; 735 PetscBool flg; 736 737 PetscFunctionBegin; 738 suffix[0] = (char *)""; 739 suffix[1] = (char *)"_X"; 740 suffix[2] = (char *)"_XX"; 741 suffix[3] = (char *)"_1"; 742 suffix[4] = (char *)"_11"; 743 *varIndex = -1; 744 745 PetscCall(PetscViewerExodusIIGetZonalVariableNames(viewer, &num_vars, &var_names)); 746 for (i = 0; i < num_vars; ++i) { 747 for (j = 0; j < num_suffix; ++j) { 748 PetscCall(PetscStrncpy(ext_name, name, MAX_STR_LENGTH)); 749 PetscCall(PetscStrlcat(ext_name, suffix[j], MAX_STR_LENGTH)); 750 PetscCall(PetscStrcasecmp(ext_name, var_names[i], &flg)); 751 if (flg) *varIndex = i; 752 } 753 if (flg) break; 754 } 755 PetscFunctionReturn(PETSC_SUCCESS); 756 } 757 758 PetscErrorCode DMView_PlexExodusII(DM dm, PetscViewer viewer) 759 { 760 enum ElemType { 761 SEGMENT, 762 TRI, 763 QUAD, 764 TET, 765 HEX 766 }; 767 MPI_Comm comm; 768 PetscInt degree; /* the order of the mesh */ 769 /* Connectivity Variables */ 770 PetscInt cellsNotInConnectivity; 771 /* Cell Sets */ 772 DMLabel csLabel; 773 IS csIS; 774 const PetscInt *csIdx; 775 PetscInt num_cs, cs; 776 enum ElemType *type; 777 PetscBool hasLabel; 778 /* Coordinate Variables */ 779 DM cdm; 780 PetscSection coordSection; 781 Vec coord; 782 PetscInt **nodes; 783 PetscInt depth, d, dim, skipCells = 0; 784 PetscInt pStart, pEnd, p, cStart, cEnd, numCells, vStart, vEnd, numVertices, eStart, eEnd, numEdges, fStart, fEnd, numFaces, numNodes; 785 PetscInt num_vs, num_fs; 786 PetscMPIInt rank, size; 787 const char *dmName; 788 PetscInt nodesLineP1[4] = {2, 0, 0, 0}; 789 PetscInt nodesLineP2[4] = {2, 0, 0, 1}; 790 PetscInt nodesTriP1[4] = {3, 0, 0, 0}; 791 PetscInt nodesTriP2[4] = {3, 3, 0, 0}; 792 PetscInt nodesQuadP1[4] = {4, 0, 0, 0}; 793 PetscInt nodesQuadP2[4] = {4, 4, 0, 1}; 794 PetscInt nodesTetP1[4] = {4, 0, 0, 0}; 795 PetscInt nodesTetP2[4] = {4, 6, 0, 0}; 796 PetscInt nodesHexP1[4] = {8, 0, 0, 0}; 797 PetscInt nodesHexP2[4] = {8, 12, 6, 1}; 798 PetscExodusIIInt CPU_word_size, IO_word_size, EXO_mode; 799 PetscExodusIIFloat EXO_version; 800 801 PetscViewer_ExodusII *exo = (PetscViewer_ExodusII *)viewer->data; 802 803 PetscFunctionBegin; 804 PetscCall(PetscObjectGetComm((PetscObject)dm, &comm)); 805 PetscCallMPI(MPI_Comm_rank(comm, &rank)); 806 PetscCallMPI(MPI_Comm_size(comm, &size)); 807 808 /* 809 Creating coordSection is a collective operation so we do it somewhat out of sequence 810 */ 811 PetscCall(PetscSectionCreate(comm, &coordSection)); 812 PetscCall(DMGetCoordinatesLocalSetUp(dm)); 813 /* 814 Check that all points are on rank 0 since we don't know how to save distributed DM in exodus format 815 */ 816 PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd)); 817 PetscCall(DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd)); 818 PetscCall(DMPlexGetDepthStratum(dm, 1, &eStart, &eEnd)); 819 PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd)); 820 numCells = cEnd - cStart; 821 numEdges = eEnd - eStart; 822 numVertices = vEnd - vStart; 823 PetscCheck(!(rank && (numCells || numEdges || numVertices)), PETSC_COMM_SELF, PETSC_ERR_SUP, "Writing distributed DM in exodusII format not supported"); 824 if (rank == 0) { 825 switch (exo->btype) { 826 case FILE_MODE_READ: 827 case FILE_MODE_APPEND: 828 case FILE_MODE_UPDATE: 829 case FILE_MODE_APPEND_UPDATE: 830 /* exodusII does not allow writing geometry to an existing file */ 831 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "cannot add geometry to existing file %s", exo->filename); 832 case FILE_MODE_WRITE: 833 /* Create an empty file if one already exists*/ 834 EXO_mode = EX_CLOBBER; 835 #if defined(PETSC_USE_64BIT_INDICES) 836 EXO_mode += EX_ALL_INT64_API; 837 #endif 838 CPU_word_size = sizeof(PetscReal); 839 IO_word_size = sizeof(PetscReal); 840 exo->exoid = ex_create(exo->filename, EXO_mode, &CPU_word_size, &IO_word_size); 841 PetscCheck(exo->exoid >= 0, PETSC_COMM_SELF, PETSC_ERR_LIB, "ex_create failed for %s", exo->filename); 842 843 break; 844 default: 845 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ORDER, "Must call PetscViewerFileSetMode() before PetscViewerFileSetName()"); 846 } 847 848 /* --- Get DM info --- */ 849 PetscCall(PetscObjectGetName((PetscObject)dm, &dmName)); 850 PetscCall(DMPlexGetDepth(dm, &depth)); 851 PetscCall(DMGetDimension(dm, &dim)); 852 PetscCall(DMPlexGetChart(dm, &pStart, &pEnd)); 853 if (depth == 3) { 854 numFaces = fEnd - fStart; 855 } else { 856 numFaces = 0; 857 } 858 PetscCall(DMGetLabelSize(dm, "Cell Sets", &num_cs)); 859 PetscCall(DMGetLabelSize(dm, "Vertex Sets", &num_vs)); 860 PetscCall(DMGetLabelSize(dm, "Face Sets", &num_fs)); 861 PetscCall(DMGetCoordinatesLocal(dm, &coord)); 862 PetscCall(DMGetCoordinateDM(dm, &cdm)); 863 if (num_cs > 0) { 864 PetscCall(DMGetLabel(dm, "Cell Sets", &csLabel)); 865 PetscCall(DMLabelGetValueIS(csLabel, &csIS)); 866 PetscCall(ISGetIndices(csIS, &csIdx)); 867 } 868 PetscCall(PetscMalloc1(num_cs, &nodes)); 869 /* Set element type for each block and compute total number of nodes */ 870 PetscCall(PetscMalloc1(num_cs, &type)); 871 numNodes = numVertices; 872 873 PetscCall(PetscViewerExodusIIGetOrder(viewer, °ree)); 874 if (degree == 2) numNodes += numEdges; 875 cellsNotInConnectivity = numCells; 876 for (cs = 0; cs < num_cs; ++cs) { 877 IS stratumIS; 878 const PetscInt *cells; 879 PetscScalar *xyz = NULL; 880 PetscInt csSize, closureSize; 881 882 PetscCall(DMLabelGetStratumIS(csLabel, csIdx[cs], &stratumIS)); 883 PetscCall(ISGetIndices(stratumIS, &cells)); 884 PetscCall(ISGetSize(stratumIS, &csSize)); 885 PetscCall(DMPlexVecGetClosure(cdm, NULL, coord, cells[0], &closureSize, &xyz)); 886 switch (dim) { 887 case 1: 888 if (closureSize == 2 * dim) { 889 type[cs] = SEGMENT; 890 } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Number of vertices %" PetscInt_FMT " in dimension %" PetscInt_FMT " has no ExodusII type", closureSize / dim, dim); 891 break; 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 = (PetscExodusIIInt)ex_inquire_int(exoid, EX_INQ_ELEM_BLK); // This is an int64_t 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 PetscCall(PetscIntCast(block.num_entry, &csSize[set])); // This is an int64_t 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 = (PetscExodusIIInt)ex_inquire_int(exoid, EX_INQ_ELEM_BLK); // This is an int64_t 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 PetscCall(PetscIntCast(block.num_entry, &csSize[set])); // This is an int64_t 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