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