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