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 Parameter: 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 Parameter: 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 - an array of string names to be set, the strings are copied into the `PetscViewer` 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 *const 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 - an array of string names to be set, the strings are copied into the `PetscViewer` 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 *const 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 Parameter: 535 . varNames - returns an array of char pointers where the zonal variable names are 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, const char *const *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 = (const char *const *)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 Parameter: 577 . varNames - returns an array of char pointers where the nodal variable names are 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, const char *const *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 = (const char *const *)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 const char *const *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 const char *const *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 static PetscErrorCode ExodusGetCellType_Internal(const char *elem_type, DMPolytopeType *ct) 1573 { 1574 PetscBool flg; 1575 1576 PetscFunctionBegin; 1577 *ct = DM_POLYTOPE_UNKNOWN; 1578 PetscCall(PetscStrcmp(elem_type, "BAR2", &flg)); 1579 if (flg) { 1580 *ct = DM_POLYTOPE_SEGMENT; 1581 goto done; 1582 } 1583 PetscCall(PetscStrcmp(elem_type, "BAR3", &flg)); 1584 if (flg) { 1585 *ct = DM_POLYTOPE_SEGMENT; 1586 goto done; 1587 } 1588 PetscCall(PetscStrcmp(elem_type, "TRI", &flg)); 1589 if (flg) { 1590 *ct = DM_POLYTOPE_TRIANGLE; 1591 goto done; 1592 } 1593 PetscCall(PetscStrcmp(elem_type, "TRI3", &flg)); 1594 if (flg) { 1595 *ct = DM_POLYTOPE_TRIANGLE; 1596 goto done; 1597 } 1598 PetscCall(PetscStrcmp(elem_type, "QUAD", &flg)); 1599 if (flg) { 1600 *ct = DM_POLYTOPE_QUADRILATERAL; 1601 goto done; 1602 } 1603 PetscCall(PetscStrcmp(elem_type, "QUAD4", &flg)); 1604 if (flg) { 1605 *ct = DM_POLYTOPE_QUADRILATERAL; 1606 goto done; 1607 } 1608 PetscCall(PetscStrcmp(elem_type, "SHELL4", &flg)); 1609 if (flg) { 1610 *ct = DM_POLYTOPE_QUADRILATERAL; 1611 goto done; 1612 } 1613 PetscCall(PetscStrcmp(elem_type, "TETRA", &flg)); 1614 if (flg) { 1615 *ct = DM_POLYTOPE_TETRAHEDRON; 1616 goto done; 1617 } 1618 PetscCall(PetscStrcmp(elem_type, "TET4", &flg)); 1619 if (flg) { 1620 *ct = DM_POLYTOPE_TETRAHEDRON; 1621 goto done; 1622 } 1623 PetscCall(PetscStrcmp(elem_type, "WEDGE", &flg)); 1624 if (flg) { 1625 *ct = DM_POLYTOPE_TRI_PRISM; 1626 goto done; 1627 } 1628 PetscCall(PetscStrcmp(elem_type, "HEX", &flg)); 1629 if (flg) { 1630 *ct = DM_POLYTOPE_HEXAHEDRON; 1631 goto done; 1632 } 1633 PetscCall(PetscStrcmp(elem_type, "HEX8", &flg)); 1634 if (flg) { 1635 *ct = DM_POLYTOPE_HEXAHEDRON; 1636 goto done; 1637 } 1638 PetscCall(PetscStrcmp(elem_type, "HEXAHEDRON", &flg)); 1639 if (flg) { 1640 *ct = DM_POLYTOPE_HEXAHEDRON; 1641 goto done; 1642 } 1643 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unrecognized element type %s", elem_type); 1644 done: 1645 PetscFunctionReturn(PETSC_SUCCESS); 1646 } 1647 1648 /*@ 1649 DMPlexCreateExodus - Create a `DMPLEX` mesh from an ExodusII file ID. 1650 1651 Collective 1652 1653 Input Parameters: 1654 + comm - The MPI communicator 1655 . exoid - The ExodusII id associated with a exodus file and obtained using ex_open 1656 - interpolate - Create faces and edges in the mesh 1657 1658 Output Parameter: 1659 . dm - The `DM` object representing the mesh 1660 1661 Level: beginner 1662 1663 .seealso: [](ch_unstructured), `DM`, `PETSCVIEWEREXODUSII`, `DMPLEX`, `DMCreate()` 1664 @*/ 1665 PetscErrorCode DMPlexCreateExodus(MPI_Comm comm, PetscExodusIIInt exoid, PetscBool interpolate, DM *dm) 1666 { 1667 PetscMPIInt num_proc, rank; 1668 DMLabel cellSets = NULL, faceSets = NULL, vertSets = NULL; 1669 PetscSection coordSection; 1670 Vec coordinates; 1671 PetscScalar *coords; 1672 PetscInt coordSize, v; 1673 /* Read from ex_get_init() */ 1674 char title[PETSC_MAX_PATH_LEN + 1]; 1675 int dim = 0, dimEmbed = 0, numVertices = 0, numCells = 0; 1676 int num_cs = 0, num_vs = 0, num_fs = 0; 1677 1678 PetscFunctionBegin; 1679 PetscCallMPI(MPI_Comm_rank(comm, &rank)); 1680 PetscCallMPI(MPI_Comm_size(comm, &num_proc)); 1681 PetscCall(DMCreate(comm, dm)); 1682 PetscCall(DMSetType(*dm, DMPLEX)); 1683 /* Open EXODUS II file and read basic information on rank 0, then broadcast to all processors */ 1684 if (rank == 0) { 1685 PetscCall(PetscMemzero(title, PETSC_MAX_PATH_LEN + 1)); 1686 PetscCallExternal(ex_get_init, exoid, title, &dimEmbed, &numVertices, &numCells, &num_cs, &num_vs, &num_fs); 1687 PetscCheck(num_cs, PETSC_COMM_SELF, PETSC_ERR_SUP, "Exodus file does not contain any cell set"); 1688 } 1689 PetscCallMPI(MPI_Bcast(title, PETSC_MAX_PATH_LEN + 1, MPI_CHAR, 0, comm)); 1690 PetscCallMPI(MPI_Bcast(&dim, 1, MPI_INT, 0, comm)); 1691 PetscCall(PetscObjectSetName((PetscObject)*dm, title)); 1692 PetscCall(DMPlexSetChart(*dm, 0, numCells + numVertices)); 1693 /* We do not want this label automatically computed, instead we compute it here */ 1694 PetscCall(DMCreateLabel(*dm, "celltype")); 1695 1696 /* Read cell sets information */ 1697 if (rank == 0) { 1698 PetscInt *cone; 1699 int c, cs, ncs, c_loc, v, v_loc; 1700 /* Read from ex_get_elem_blk_ids() */ 1701 int *cs_id, *cs_order; 1702 /* Read from ex_get_elem_block() */ 1703 char buffer[PETSC_MAX_PATH_LEN + 1]; 1704 int num_cell_in_set, num_vertex_per_cell, num_hybrid, num_attr; 1705 /* Read from ex_get_elem_conn() */ 1706 int *cs_connect; 1707 1708 /* Get cell sets IDs */ 1709 PetscCall(PetscMalloc2(num_cs, &cs_id, num_cs, &cs_order)); 1710 PetscCallExternal(ex_get_ids, exoid, EX_ELEM_BLOCK, cs_id); 1711 /* Read the cell set connectivity table and build mesh topology 1712 EXO standard requires that cells in cell sets be numbered sequentially and be pairwise disjoint. */ 1713 /* Check for a hybrid mesh */ 1714 for (cs = 0, num_hybrid = 0; cs < num_cs; ++cs) { 1715 DMPolytopeType ct; 1716 char elem_type[PETSC_MAX_PATH_LEN]; 1717 1718 PetscCall(PetscArrayzero(elem_type, sizeof(elem_type))); 1719 PetscCallExternal(ex_get_elem_type, exoid, cs_id[cs], elem_type); 1720 PetscCall(ExodusGetCellType_Internal(elem_type, &ct)); 1721 dim = PetscMax(dim, DMPolytopeTypeGetDim(ct)); 1722 PetscCallExternal(ex_get_block, exoid, EX_ELEM_BLOCK, cs_id[cs], buffer, &num_cell_in_set, &num_vertex_per_cell, 0, 0, &num_attr); 1723 switch (ct) { 1724 case DM_POLYTOPE_TRI_PRISM: 1725 cs_order[cs] = cs; 1726 ++num_hybrid; 1727 break; 1728 default: 1729 for (c = cs; c > cs - num_hybrid; --c) cs_order[c] = cs_order[c - 1]; 1730 cs_order[cs - num_hybrid] = cs; 1731 } 1732 } 1733 /* First set sizes */ 1734 for (ncs = 0, c = 0; ncs < num_cs; ++ncs) { 1735 DMPolytopeType ct; 1736 char elem_type[PETSC_MAX_PATH_LEN]; 1737 const PetscInt cs = cs_order[ncs]; 1738 1739 PetscCall(PetscArrayzero(elem_type, sizeof(elem_type))); 1740 PetscCallExternal(ex_get_elem_type, exoid, cs_id[cs], elem_type); 1741 PetscCall(ExodusGetCellType_Internal(elem_type, &ct)); 1742 PetscCallExternal(ex_get_block, exoid, EX_ELEM_BLOCK, cs_id[cs], buffer, &num_cell_in_set, &num_vertex_per_cell, 0, 0, &num_attr); 1743 for (c_loc = 0; c_loc < num_cell_in_set; ++c_loc, ++c) { 1744 PetscCall(DMPlexSetConeSize(*dm, c, num_vertex_per_cell)); 1745 PetscCall(DMPlexSetCellType(*dm, c, ct)); 1746 } 1747 } 1748 for (v = numCells; v < numCells + numVertices; ++v) PetscCall(DMPlexSetCellType(*dm, v, DM_POLYTOPE_POINT)); 1749 PetscCall(DMSetUp(*dm)); 1750 for (ncs = 0, c = 0; ncs < num_cs; ++ncs) { 1751 const PetscInt cs = cs_order[ncs]; 1752 PetscCallExternal(ex_get_block, exoid, EX_ELEM_BLOCK, cs_id[cs], buffer, &num_cell_in_set, &num_vertex_per_cell, 0, 0, &num_attr); 1753 PetscCall(PetscMalloc2(num_vertex_per_cell * num_cell_in_set, &cs_connect, num_vertex_per_cell, &cone)); 1754 PetscCallExternal(ex_get_conn, exoid, EX_ELEM_BLOCK, cs_id[cs], cs_connect, NULL, NULL); 1755 /* EXO uses Fortran-based indexing, DMPlex uses C-style and numbers cell first then vertices. */ 1756 for (c_loc = 0, v = 0; c_loc < num_cell_in_set; ++c_loc, ++c) { 1757 DMPolytopeType ct; 1758 1759 for (v_loc = 0; v_loc < num_vertex_per_cell; ++v_loc, ++v) cone[v_loc] = cs_connect[v] + numCells - 1; 1760 PetscCall(DMPlexGetCellType(*dm, c, &ct)); 1761 PetscCall(DMPlexInvertCell(ct, cone)); 1762 PetscCall(DMPlexSetCone(*dm, c, cone)); 1763 PetscCall(DMSetLabelValue_Fast(*dm, &cellSets, "Cell Sets", c, cs_id[cs])); 1764 } 1765 PetscCall(PetscFree2(cs_connect, cone)); 1766 } 1767 PetscCall(PetscFree2(cs_id, cs_order)); 1768 } 1769 { 1770 PetscInt ints[] = {dim, dimEmbed}; 1771 1772 PetscCallMPI(MPI_Bcast(ints, 2, MPIU_INT, 0, comm)); 1773 PetscCall(DMSetDimension(*dm, ints[0])); 1774 PetscCall(DMSetCoordinateDim(*dm, ints[1])); 1775 dim = ints[0]; 1776 dimEmbed = ints[1]; 1777 } 1778 PetscCall(DMPlexSymmetrize(*dm)); 1779 PetscCall(DMPlexStratify(*dm)); 1780 if (interpolate) { 1781 DM idm; 1782 1783 PetscCall(DMPlexInterpolate(*dm, &idm)); 1784 PetscCall(DMDestroy(dm)); 1785 *dm = idm; 1786 } 1787 1788 /* Create vertex set label */ 1789 if (rank == 0 && (num_vs > 0)) { 1790 int vs, v; 1791 /* Read from ex_get_node_set_ids() */ 1792 int *vs_id; 1793 /* Read from ex_get_node_set_param() */ 1794 int num_vertex_in_set; 1795 /* Read from ex_get_node_set() */ 1796 int *vs_vertex_list; 1797 1798 /* Get vertex set ids */ 1799 PetscCall(PetscMalloc1(num_vs, &vs_id)); 1800 PetscCallExternal(ex_get_ids, exoid, EX_NODE_SET, vs_id); 1801 for (vs = 0; vs < num_vs; ++vs) { 1802 PetscCallExternal(ex_get_set_param, exoid, EX_NODE_SET, vs_id[vs], &num_vertex_in_set, NULL); 1803 PetscCall(PetscMalloc1(num_vertex_in_set, &vs_vertex_list)); 1804 PetscCallExternal(ex_get_set, exoid, EX_NODE_SET, vs_id[vs], vs_vertex_list, NULL); 1805 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])); 1806 PetscCall(PetscFree(vs_vertex_list)); 1807 } 1808 PetscCall(PetscFree(vs_id)); 1809 } 1810 /* Read coordinates */ 1811 PetscCall(DMGetCoordinateSection(*dm, &coordSection)); 1812 PetscCall(PetscSectionSetNumFields(coordSection, 1)); 1813 PetscCall(PetscSectionSetFieldComponents(coordSection, 0, dimEmbed)); 1814 PetscCall(PetscSectionSetChart(coordSection, numCells, numCells + numVertices)); 1815 for (v = numCells; v < numCells + numVertices; ++v) { 1816 PetscCall(PetscSectionSetDof(coordSection, v, dimEmbed)); 1817 PetscCall(PetscSectionSetFieldDof(coordSection, v, 0, dimEmbed)); 1818 } 1819 PetscCall(PetscSectionSetUp(coordSection)); 1820 PetscCall(PetscSectionGetStorageSize(coordSection, &coordSize)); 1821 PetscCall(VecCreate(PETSC_COMM_SELF, &coordinates)); 1822 PetscCall(PetscObjectSetName((PetscObject)coordinates, "coordinates")); 1823 PetscCall(VecSetSizes(coordinates, coordSize, PETSC_DETERMINE)); 1824 PetscCall(VecSetBlockSize(coordinates, dimEmbed)); 1825 PetscCall(VecSetType(coordinates, VECSTANDARD)); 1826 PetscCall(VecGetArray(coordinates, &coords)); 1827 if (rank == 0) { 1828 PetscReal *x, *y, *z; 1829 1830 PetscCall(PetscMalloc3(numVertices, &x, numVertices, &y, numVertices, &z)); 1831 PetscCallExternal(ex_get_coord, exoid, x, y, z); 1832 if (dimEmbed > 0) { 1833 for (v = 0; v < numVertices; ++v) coords[v * dimEmbed + 0] = x[v]; 1834 } 1835 if (dimEmbed > 1) { 1836 for (v = 0; v < numVertices; ++v) coords[v * dimEmbed + 1] = y[v]; 1837 } 1838 if (dimEmbed > 2) { 1839 for (v = 0; v < numVertices; ++v) coords[v * dimEmbed + 2] = z[v]; 1840 } 1841 PetscCall(PetscFree3(x, y, z)); 1842 } 1843 PetscCall(VecRestoreArray(coordinates, &coords)); 1844 PetscCall(DMSetCoordinatesLocal(*dm, coordinates)); 1845 PetscCall(VecDestroy(&coordinates)); 1846 1847 /* Create side set label */ 1848 if (rank == 0 && interpolate && (num_fs > 0)) { 1849 int fs, f, voff; 1850 /* Read from ex_get_side_set_ids() */ 1851 int *fs_id; 1852 /* Read from ex_get_side_set_param() */ 1853 int num_side_in_set; 1854 /* Read from ex_get_side_set_node_list() */ 1855 int *fs_vertex_count_list, *fs_vertex_list, *fs_side_list; 1856 /* Read side set labels */ 1857 char fs_name[MAX_STR_LENGTH + 1]; 1858 size_t fs_name_len; 1859 1860 /* Get side set ids */ 1861 PetscCall(PetscMalloc1(num_fs, &fs_id)); 1862 PetscCallExternal(ex_get_ids, exoid, EX_SIDE_SET, fs_id); 1863 // Ids 1 and 2 are reserved by ExodusII for indicating things in 3D 1864 for (fs = 0; fs < num_fs; ++fs) { 1865 PetscCallExternal(ex_get_set_param, exoid, EX_SIDE_SET, fs_id[fs], &num_side_in_set, NULL); 1866 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)); 1867 PetscCallExternal(ex_get_side_set_node_list, exoid, fs_id[fs], fs_vertex_count_list, fs_vertex_list); 1868 PetscCallExternal(ex_get_set, exoid, EX_SIDE_SET, fs_id[fs], NULL, fs_side_list); 1869 1870 /* Get the specific name associated with this side set ID. */ 1871 int fs_name_err = ex_get_name(exoid, EX_SIDE_SET, fs_id[fs], fs_name); 1872 if (!fs_name_err) { 1873 PetscCall(PetscStrlen(fs_name, &fs_name_len)); 1874 if (fs_name_len == 0) PetscCall(PetscStrncpy(fs_name, "Face Sets", MAX_STR_LENGTH)); 1875 } 1876 for (f = 0, voff = 0; f < num_side_in_set; ++f) { 1877 const PetscInt *faces = NULL; 1878 PetscInt faceSize = fs_vertex_count_list[f], numFaces; 1879 PetscInt faceVertices[4], v; 1880 1881 PetscCheck(faceSize <= 4, comm, PETSC_ERR_ARG_WRONG, "ExodusII side cannot have %" PetscInt_FMT " > 4 vertices", faceSize); 1882 for (v = 0; v < faceSize; ++v, ++voff) faceVertices[v] = fs_vertex_list[voff] + numCells - 1; 1883 PetscCall(DMPlexGetFullJoin(*dm, faceSize, faceVertices, &numFaces, &faces)); 1884 PetscCheck(numFaces == 1, comm, PETSC_ERR_ARG_WRONG, "Invalid ExodusII side %d in set %d maps to %" PetscInt_FMT " faces", f, fs, numFaces); 1885 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]); 1886 PetscCall(DMSetLabelValue_Fast(*dm, &faceSets, "Face Sets", faces[0], fs_id[fs])); 1887 /* Only add the label if one has been detected for this side set. */ 1888 if (!fs_name_err) PetscCall(DMSetLabelValue(*dm, fs_name, faces[0], fs_id[fs])); 1889 PetscCall(DMPlexRestoreJoin(*dm, faceSize, faceVertices, &numFaces, &faces)); 1890 } 1891 PetscCall(PetscFree3(fs_vertex_count_list, fs_vertex_list, fs_side_list)); 1892 } 1893 PetscCall(PetscFree(fs_id)); 1894 } 1895 1896 { /* Create Cell/Face/Vertex Sets labels at all processes */ 1897 enum { 1898 n = 3 1899 }; 1900 PetscBool flag[n]; 1901 1902 flag[0] = cellSets ? PETSC_TRUE : PETSC_FALSE; 1903 flag[1] = faceSets ? PETSC_TRUE : PETSC_FALSE; 1904 flag[2] = vertSets ? PETSC_TRUE : PETSC_FALSE; 1905 PetscCallMPI(MPI_Bcast(flag, n, MPIU_BOOL, 0, comm)); 1906 if (flag[0]) PetscCall(DMCreateLabel(*dm, "Cell Sets")); 1907 if (flag[1]) PetscCall(DMCreateLabel(*dm, "Face Sets")); 1908 if (flag[2]) PetscCall(DMCreateLabel(*dm, "Vertex Sets")); 1909 } 1910 PetscFunctionReturn(PETSC_SUCCESS); 1911 } 1912