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