1 #define PETSCDM_DLL 2 #include <petsc/private/dmpleximpl.h> /*I "petscdmplex.h" I*/ 3 4 #if defined(PETSC_HAVE_EXODUSII) 5 #include <netcdf.h> 6 #include <exodusII.h> 7 #endif 8 9 #include <petsc/private/viewerimpl.h> 10 #include <petsc/private/viewerexodusiiimpl.h> 11 #if defined(PETSC_HAVE_EXODUSII) 12 /* 13 PETSC_VIEWER_EXODUSII_ - Creates an ExodusII PetscViewer shared by all processors in a communicator. 14 15 Collective 16 17 Input Parameter: 18 . comm - the MPI communicator to share the ExodusII PetscViewer 19 20 Level: intermediate 21 22 Notes: 23 misses Fortran bindings 24 25 Notes: 26 Unlike almost all other PETSc routines, PETSC_VIEWER_EXODUSII_ does not return 27 an error code. The GLVIS PetscViewer is usually used in the form 28 $ XXXView(XXX object, PETSC_VIEWER_EXODUSII_(comm)); 29 30 .seealso: PetscViewerExodusIIOpen(), PetscViewerType, PetscViewerCreate(), PetscViewerDestroy() 31 */ 32 PetscViewer PETSC_VIEWER_EXODUSII_(MPI_Comm comm) 33 { 34 PetscViewer viewer; 35 PetscErrorCode ierr; 36 37 PetscFunctionBegin; 38 ierr = PetscViewerExodusIIOpen(comm, "mesh.exo", FILE_MODE_WRITE, &viewer); 39 if (ierr) {PetscError(PETSC_COMM_SELF,__LINE__,"PETSC_VIEWER_EXODUSII_",__FILE__,PETSC_ERR_PLIB,PETSC_ERROR_INITIAL," ");PetscFunctionReturn(0);} 40 ierr = PetscObjectRegisterDestroy((PetscObject) viewer); 41 if (ierr) {PetscError(PETSC_COMM_SELF,__LINE__,"PETSC_VIEWER_EXODUSII_",__FILE__,PETSC_ERR_PLIB,PETSC_ERROR_INITIAL," ");PetscFunctionReturn(0);} 42 PetscFunctionReturn(viewer); 43 } 44 45 static PetscErrorCode PetscViewerView_ExodusII(PetscViewer v, PetscViewer viewer) 46 { 47 PetscViewer_ExodusII *exo = (PetscViewer_ExodusII *) v->data; 48 PetscErrorCode ierr; 49 50 PetscFunctionBegin; 51 if (exo->filename) {ierr = PetscViewerASCIIPrintf(viewer, "Filename: %s\n", exo->filename);CHKERRQ(ierr);} 52 if (exo->exoid) {ierr = PetscViewerASCIIPrintf(viewer, "exoid: %d\n", exo->exoid) ;CHKERRQ(ierr);} 53 if (exo->btype) {ierr = PetscViewerASCIIPrintf(viewer, "IO Mode: %d\n", exo->btype) ;CHKERRQ(ierr);} 54 if (exo->order) {ierr = PetscViewerASCIIPrintf(viewer, "Mesh order: %d\n", exo->order) ;CHKERRQ(ierr);} 55 PetscFunctionReturn(0); 56 } 57 58 static PetscErrorCode PetscViewerSetFromOptions_ExodusII(PetscOptionItems *PetscOptionsObject, PetscViewer v) 59 { 60 PetscErrorCode ierr; 61 62 PetscFunctionBegin; 63 ierr = PetscOptionsHead(PetscOptionsObject, "ExodusII PetscViewer Options");CHKERRQ(ierr); 64 ierr = PetscOptionsTail();CHKERRQ(ierr); 65 PetscFunctionReturn(0); 66 } 67 68 static PetscErrorCode PetscViewerSetUp_ExodusII(PetscViewer viewer) 69 { 70 PetscFunctionBegin; 71 PetscFunctionReturn(0); 72 } 73 74 static PetscErrorCode PetscViewerDestroy_ExodusII(PetscViewer viewer) 75 { 76 PetscViewer_ExodusII *exo = (PetscViewer_ExodusII *) viewer->data; 77 PetscErrorCode ierr; 78 79 PetscFunctionBegin; 80 if (exo->exoid >= 0) {PetscStackCallStandard(ex_close,(exo->exoid));} 81 ierr = PetscFree(exo->filename);CHKERRQ(ierr); 82 ierr = PetscFree(exo);CHKERRQ(ierr); 83 ierr = PetscObjectComposeFunction((PetscObject)viewer,"PetscViewerFileSetName_C",NULL);CHKERRQ(ierr); 84 ierr = PetscObjectComposeFunction((PetscObject)viewer,"PetscViewerFileGetName_C",NULL);CHKERRQ(ierr); 85 ierr = PetscObjectComposeFunction((PetscObject)viewer,"PetscViewerFileSetMode_C",NULL);CHKERRQ(ierr); 86 ierr = PetscObjectComposeFunction((PetscObject)viewer,"PetscViewerGetId_C",NULL);CHKERRQ(ierr); 87 ierr = PetscObjectComposeFunction((PetscObject)viewer,"PetscViewerGetOrder_C",NULL);CHKERRQ(ierr); 88 ierr = PetscObjectComposeFunction((PetscObject)viewer,"PetscViewerSetOrder_C",NULL);CHKERRQ(ierr); 89 PetscFunctionReturn(0); 90 } 91 92 static PetscErrorCode PetscViewerFileSetName_ExodusII(PetscViewer viewer, const char name[]) 93 { 94 PetscViewer_ExodusII *exo = (PetscViewer_ExodusII *) viewer->data; 95 PetscMPIInt rank; 96 int CPU_word_size, IO_word_size, EXO_mode; 97 PetscErrorCode ierr; 98 MPI_Info mpi_info = MPI_INFO_NULL; 99 float EXO_version; 100 101 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject) viewer), &rank);CHKERRMPI(ierr); 102 CPU_word_size = sizeof(PetscReal); 103 IO_word_size = sizeof(PetscReal); 104 105 PetscFunctionBegin; 106 if (exo->exoid >= 0) { 107 PetscStackCallStandard(ex_close,(exo->exoid)); 108 exo->exoid = -1; 109 } 110 if (exo->filename) {ierr = PetscFree(exo->filename);CHKERRQ(ierr);} 111 ierr = PetscStrallocpy(name, &exo->filename);CHKERRQ(ierr); 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(0); 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 if (exo->exoid < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_LIB, "ex_open_par failed for %s", name); 136 PetscFunctionReturn(0); 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(0); 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(0); 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(0); 164 } 165 166 static PetscErrorCode PetscViewerExodusIIGetId_ExodusII(PetscViewer viewer, int *exoid) 167 { 168 PetscViewer_ExodusII *exo = (PetscViewer_ExodusII *) viewer->data; 169 170 PetscFunctionBegin; 171 *exoid = exo->exoid; 172 PetscFunctionReturn(0); 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(0); 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(0); 191 } 192 193 /*MC 194 PETSCVIEWEREXODUSII - A viewer that writes to an Exodus II file 195 196 197 .seealso: PetscViewerExodusIIOpen(), PetscViewerCreate(), PETSCVIEWERBINARY, PETSCVIEWERHDF5, DMView(), 198 PetscViewerFileSetName(), PetscViewerFileSetMode(), PetscViewerFormat, PetscViewerType, PetscViewerSetType() 199 200 Level: beginner 201 M*/ 202 203 PETSC_EXTERN PetscErrorCode PetscViewerCreate_ExodusII(PetscViewer v) 204 { 205 PetscViewer_ExodusII *exo; 206 PetscErrorCode ierr; 207 208 PetscFunctionBegin; 209 ierr = PetscNewLog(v,&exo);CHKERRQ(ierr); 210 211 v->data = (void*) exo; 212 v->ops->destroy = PetscViewerDestroy_ExodusII; 213 v->ops->setfromoptions = PetscViewerSetFromOptions_ExodusII; 214 v->ops->setup = PetscViewerSetUp_ExodusII; 215 v->ops->view = PetscViewerView_ExodusII; 216 v->ops->flush = 0; 217 exo->btype = (PetscFileMode) -1; 218 exo->filename = 0; 219 exo->exoid = -1; 220 221 ierr = PetscObjectComposeFunction((PetscObject)v,"PetscViewerFileSetName_C",PetscViewerFileSetName_ExodusII);CHKERRQ(ierr); 222 ierr = PetscObjectComposeFunction((PetscObject)v,"PetscViewerFileGetName_C",PetscViewerFileGetName_ExodusII);CHKERRQ(ierr); 223 ierr = PetscObjectComposeFunction((PetscObject)v,"PetscViewerFileSetMode_C",PetscViewerFileSetMode_ExodusII);CHKERRQ(ierr); 224 ierr = PetscObjectComposeFunction((PetscObject)v,"PetscViewerFileGetMode_C",PetscViewerFileGetMode_ExodusII);CHKERRQ(ierr); 225 ierr = PetscObjectComposeFunction((PetscObject)v,"PetscViewerGetId_C",PetscViewerExodusIIGetId_ExodusII);CHKERRQ(ierr); 226 ierr = PetscObjectComposeFunction((PetscObject)v,"PetscViewerSetOrder_C",PetscViewerExodusIISetOrder_ExodusII);CHKERRQ(ierr); 227 ierr = PetscObjectComposeFunction((PetscObject)v,"PetscViewerGetOrder_C",PetscViewerExodusIIGetOrder_ExodusII);CHKERRQ(ierr); 228 PetscFunctionReturn(0); 229 } 230 231 /* 232 EXOGetVarIndex - Locate a result in an exodus file based on its name 233 234 Collective 235 236 Input Parameters: 237 + exoid - the exodus id of a file (obtained from ex_open or ex_create for instance) 238 . obj_type - the type of entity for instance EX_NODAL, EX_ELEM_BLOCK 239 - name - the name of the result 240 241 Output Parameters: 242 . varIndex - the location in the exodus file of the result 243 244 Notes: 245 The exodus variable index is obtained by comparing name and the 246 names of zonal variables declared in the exodus file. For instance if name is "V" 247 the location in the exodus file will be the first match of "V", "V_X", "V_XX", "V_1", or "V_11" 248 amongst all variables of type obj_type. 249 250 Level: beginner 251 252 .seealso: EXOGetVarIndex(),DMPlexView_ExodusII_Internal(),VecViewPlex_ExodusII_Nodal_Internal(),VecLoadNodal_PlexEXO(),VecLoadZonal_PlexEXO() 253 */ 254 PetscErrorCode EXOGetVarIndex_Internal(int exoid, ex_entity_type obj_type, const char name[], int *varIndex) 255 { 256 int num_vars, i, j; 257 char ext_name[MAX_STR_LENGTH+1], var_name[MAX_STR_LENGTH+1]; 258 const int num_suffix = 5; 259 char *suffix[5]; 260 PetscBool flg; 261 PetscErrorCode ierr; 262 263 PetscFunctionBegin; 264 suffix[0] = (char *) ""; 265 suffix[1] = (char *) "_X"; 266 suffix[2] = (char *) "_XX"; 267 suffix[3] = (char *) "_1"; 268 suffix[4] = (char *) "_11"; 269 270 *varIndex = -1; 271 PetscStackCallStandard(ex_get_variable_param,(exoid, obj_type, &num_vars)); 272 for (i = 0; i < num_vars; ++i) { 273 PetscStackCallStandard(ex_get_variable_name,(exoid, obj_type, i+1, var_name)); 274 for (j = 0; j < num_suffix; ++j){ 275 ierr = PetscStrncpy(ext_name, name, MAX_STR_LENGTH);CHKERRQ(ierr); 276 ierr = PetscStrlcat(ext_name, suffix[j], MAX_STR_LENGTH);CHKERRQ(ierr); 277 ierr = PetscStrcasecmp(ext_name, var_name, &flg); 278 if (flg) { 279 *varIndex = i+1; 280 } 281 } 282 } 283 PetscFunctionReturn(0); 284 } 285 286 /* 287 DMView_PlexExodusII - Write a DM to disk in exodus format 288 289 Collective on dm 290 291 Input Parameters: 292 + dm - The dm to be written 293 . viewer - an exodusII viewer 294 295 Notes: 296 Not all DM can be written to disk this way. For instance, exodus assume that element blocks (mapped to "Cell sets" labels) 297 consists of sequentially numbered cells. If this is not the case, the exodus file will be corrupted. 298 299 If the dm has been distributed, only the part of the DM on MPI rank 0 (including "ghost" cells and vertices) 300 will be written. 301 302 DMPlex only represents geometry while most post-processing software expect that a mesh also provides information 303 on the discretization space. This function assumes that the file represents Lagrange finite elements of order 1 or 2. 304 The order of the mesh shall be set using PetscViewerExodusIISetOrder 305 It should be extended to use PetscFE objects. 306 307 This function will only handle TRI, TET, QUAD, and HEX cells. 308 Level: beginner 309 310 .seealso: 311 */ 312 PetscErrorCode DMView_PlexExodusII(DM dm, PetscViewer viewer) 313 { 314 enum ElemType {TRI, QUAD, TET, HEX}; 315 MPI_Comm comm; 316 PetscInt degree; /* the order of the mesh */ 317 /* Connectivity Variables */ 318 PetscInt cellsNotInConnectivity; 319 /* Cell Sets */ 320 DMLabel csLabel; 321 IS csIS; 322 const PetscInt *csIdx; 323 PetscInt num_cs, cs; 324 enum ElemType *type; 325 PetscBool hasLabel; 326 /* Coordinate Variables */ 327 DM cdm; 328 PetscSection coordSection; 329 Vec coord; 330 PetscInt **nodes; 331 PetscInt depth, d, dim, skipCells = 0; 332 PetscInt pStart, pEnd, p, cStart, cEnd, numCells, vStart, vEnd, numVertices, eStart, eEnd, numEdges, fStart, fEnd, numFaces, numNodes; 333 PetscInt num_vs, num_fs; 334 PetscMPIInt rank, size; 335 const char *dmName; 336 PetscInt nodesTriP1[4] = {3,0,0,0}; 337 PetscInt nodesTriP2[4] = {3,3,0,0}; 338 PetscInt nodesQuadP1[4] = {4,0,0,0}; 339 PetscInt nodesQuadP2[4] = {4,4,0,1}; 340 PetscInt nodesTetP1[4] = {4,0,0,0}; 341 PetscInt nodesTetP2[4] = {4,6,0,0}; 342 PetscInt nodesHexP1[4] = {8,0,0,0}; 343 PetscInt nodesHexP2[4] = {8,12,6,1}; 344 PetscErrorCode ierr; 345 346 int CPU_word_size, IO_word_size, EXO_mode; 347 float EXO_version; 348 349 PetscViewer_ExodusII *exo = (PetscViewer_ExodusII *) viewer->data; 350 351 PetscFunctionBegin; 352 ierr = PetscObjectGetComm((PetscObject)dm, &comm);CHKERRQ(ierr); 353 ierr = MPI_Comm_rank(comm, &rank);CHKERRMPI(ierr); 354 ierr = MPI_Comm_size(comm, &size);CHKERRMPI(ierr); 355 356 /* 357 Creating coordSection is a collective operation so we do it somewhat out of sequence 358 */ 359 ierr = PetscSectionCreate(comm, &coordSection);CHKERRQ(ierr); 360 ierr = DMGetCoordinatesLocalSetUp(dm);CHKERRQ(ierr); 361 if (!rank) { 362 switch (exo->btype) { 363 case FILE_MODE_READ: 364 case FILE_MODE_APPEND: 365 case FILE_MODE_UPDATE: 366 case FILE_MODE_APPEND_UPDATE: 367 /* exodusII does not allow writing geometry to an existing file */ 368 SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_LIB, "cannot add geometry to existing file %s", exo->filename); 369 case FILE_MODE_WRITE: 370 /* Create an empty file if one already exists*/ 371 EXO_mode = EX_CLOBBER; 372 #if defined(PETSC_USE_64BIT_INDICES) 373 EXO_mode += EX_ALL_INT64_API; 374 #endif 375 CPU_word_size = sizeof(PetscReal); 376 IO_word_size = sizeof(PetscReal); 377 exo->exoid = ex_create(exo->filename, EXO_mode, &CPU_word_size, &IO_word_size); 378 if (exo->exoid < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_LIB, "ex_create failed for %s", exo->filename); 379 380 break; 381 default: 382 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER, "Must call PetscViewerFileSetMode() before PetscViewerFileSetName()"); 383 } 384 385 /* --- Get DM info --- */ 386 ierr = PetscObjectGetName((PetscObject) dm, &dmName);CHKERRQ(ierr); 387 ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr); 388 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 389 ierr = DMPlexGetChart(dm, &pStart, &pEnd);CHKERRQ(ierr); 390 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 391 ierr = DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);CHKERRQ(ierr); 392 ierr = DMPlexGetDepthStratum(dm, 1, &eStart, &eEnd);CHKERRQ(ierr); 393 ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr); 394 numCells = cEnd - cStart; 395 numEdges = eEnd - eStart; 396 numVertices = vEnd - vStart; 397 if (depth == 3) {numFaces = fEnd - fStart;} 398 else {numFaces = 0;} 399 ierr = DMGetLabelSize(dm, "Cell Sets", &num_cs);CHKERRQ(ierr); 400 ierr = DMGetLabelSize(dm, "Vertex Sets", &num_vs);CHKERRQ(ierr); 401 ierr = DMGetLabelSize(dm, "Face Sets", &num_fs);CHKERRQ(ierr); 402 ierr = DMGetCoordinatesLocal(dm, &coord);CHKERRQ(ierr); 403 ierr = DMGetCoordinateDM(dm, &cdm);CHKERRQ(ierr); 404 if (num_cs > 0) { 405 ierr = DMGetLabel(dm, "Cell Sets", &csLabel);CHKERRQ(ierr); 406 ierr = DMLabelGetValueIS(csLabel, &csIS);CHKERRQ(ierr); 407 ierr = ISGetIndices(csIS, &csIdx);CHKERRQ(ierr); 408 } 409 ierr = PetscMalloc1(num_cs, &nodes);CHKERRQ(ierr); 410 /* Set element type for each block and compute total number of nodes */ 411 ierr = PetscMalloc1(num_cs, &type);CHKERRQ(ierr); 412 numNodes = numVertices; 413 414 ierr = PetscViewerExodusIIGetOrder(viewer, °ree);CHKERRQ(ierr); 415 if (degree == 2) {numNodes += numEdges;} 416 cellsNotInConnectivity = numCells; 417 for (cs = 0; cs < num_cs; ++cs) { 418 IS stratumIS; 419 const PetscInt *cells; 420 PetscScalar *xyz = NULL; 421 PetscInt csSize, closureSize; 422 423 ierr = DMLabelGetStratumIS(csLabel, csIdx[cs], &stratumIS);CHKERRQ(ierr); 424 ierr = ISGetIndices(stratumIS, &cells);CHKERRQ(ierr); 425 ierr = ISGetSize(stratumIS, &csSize);CHKERRQ(ierr); 426 ierr = DMPlexVecGetClosure(cdm, NULL, coord, cells[0], &closureSize, &xyz);CHKERRQ(ierr); 427 switch (dim) { 428 case 2: 429 if (closureSize == 3*dim) {type[cs] = TRI;} 430 else if (closureSize == 4*dim) {type[cs] = QUAD;} 431 else SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Number of vertices %D in dimension %D has no ExodusII type", closureSize/dim, dim); 432 break; 433 case 3: 434 if (closureSize == 4*dim) {type[cs] = TET;} 435 else if (closureSize == 8*dim) {type[cs] = HEX;} 436 else SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Number of vertices %D in dimension %D has no ExodusII type", closureSize/dim, dim); 437 break; 438 default: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Dimension %D not handled by ExodusII viewer", dim); 439 } 440 if ((degree == 2) && (type[cs] == QUAD)) {numNodes += csSize;} 441 if ((degree == 2) && (type[cs] == HEX)) {numNodes += csSize; numNodes += numFaces;} 442 ierr = DMPlexVecRestoreClosure(cdm, NULL, coord, cells[0], &closureSize, &xyz);CHKERRQ(ierr); 443 /* Set nodes and Element type */ 444 if (type[cs] == TRI) { 445 if (degree == 1) nodes[cs] = nodesTriP1; 446 else if (degree == 2) nodes[cs] = nodesTriP2; 447 } else if (type[cs] == QUAD) { 448 if (degree == 1) nodes[cs] = nodesQuadP1; 449 else if (degree == 2) nodes[cs] = nodesQuadP2; 450 } else if (type[cs] == TET) { 451 if (degree == 1) nodes[cs] = nodesTetP1; 452 else if (degree == 2) nodes[cs] = nodesTetP2; 453 } else if (type[cs] == HEX) { 454 if (degree == 1) nodes[cs] = nodesHexP1; 455 else if (degree == 2) nodes[cs] = nodesHexP2; 456 } 457 /* Compute the number of cells not in the connectivity table */ 458 cellsNotInConnectivity -= nodes[cs][3]*csSize; 459 460 ierr = ISRestoreIndices(stratumIS, &cells);CHKERRQ(ierr); 461 ierr = ISDestroy(&stratumIS);CHKERRQ(ierr); 462 } 463 if (num_cs > 0) {PetscStackCallStandard(ex_put_init,(exo->exoid, dmName, dim, numNodes, numCells, num_cs, num_vs, num_fs));} 464 /* --- Connectivity --- */ 465 for (cs = 0; cs < num_cs; ++cs) { 466 IS stratumIS; 467 const PetscInt *cells; 468 PetscInt *connect, off = 0; 469 PetscInt edgesInClosure = 0, facesInClosure = 0, verticesInClosure = 0; 470 PetscInt csSize, c, connectSize, closureSize; 471 char *elem_type = NULL; 472 char elem_type_tri3[] = "TRI3", elem_type_quad4[] = "QUAD4"; 473 char elem_type_tri6[] = "TRI6", elem_type_quad9[] = "QUAD9"; 474 char elem_type_tet4[] = "TET4", elem_type_hex8[] = "HEX8"; 475 char elem_type_tet10[] = "TET10", elem_type_hex27[] = "HEX27"; 476 477 ierr = DMLabelGetStratumIS(csLabel, csIdx[cs], &stratumIS);CHKERRQ(ierr); 478 ierr = ISGetIndices(stratumIS, &cells);CHKERRQ(ierr); 479 ierr = ISGetSize(stratumIS, &csSize);CHKERRQ(ierr); 480 /* Set Element type */ 481 if (type[cs] == TRI) { 482 if (degree == 1) elem_type = elem_type_tri3; 483 else if (degree == 2) elem_type = elem_type_tri6; 484 } else if (type[cs] == QUAD) { 485 if (degree == 1) elem_type = elem_type_quad4; 486 else if (degree == 2) elem_type = elem_type_quad9; 487 } else if (type[cs] == TET) { 488 if (degree == 1) elem_type = elem_type_tet4; 489 else if (degree == 2) elem_type = elem_type_tet10; 490 } else if (type[cs] == HEX) { 491 if (degree == 1) elem_type = elem_type_hex8; 492 else if (degree == 2) elem_type = elem_type_hex27; 493 } 494 connectSize = nodes[cs][0] + nodes[cs][1] + nodes[cs][2] + nodes[cs][3]; 495 ierr = PetscMalloc1(PetscMax(27,connectSize)*csSize, &connect);CHKERRQ(ierr); 496 PetscStackCallStandard(ex_put_block,(exo->exoid, EX_ELEM_BLOCK, csIdx[cs], elem_type, csSize, connectSize, 0, 0, 1)); 497 /* Find number of vertices, edges, and faces in the closure */ 498 verticesInClosure = nodes[cs][0]; 499 if (depth > 1) { 500 if (dim == 2) { 501 ierr = DMPlexGetConeSize(dm, cells[0], &edgesInClosure);CHKERRQ(ierr); 502 } else if (dim == 3) { 503 PetscInt *closure = NULL; 504 505 ierr = DMPlexGetConeSize(dm, cells[0], &facesInClosure);CHKERRQ(ierr); 506 ierr = DMPlexGetTransitiveClosure(dm, cells[0], PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 507 edgesInClosure = closureSize - facesInClosure - 1 - verticesInClosure; 508 ierr = DMPlexRestoreTransitiveClosure(dm, cells[0], PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 509 } 510 } 511 /* Get connectivity for each cell */ 512 for (c = 0; c < csSize; ++c) { 513 PetscInt *closure = NULL; 514 PetscInt temp, i; 515 516 ierr = DMPlexGetTransitiveClosure(dm, cells[c], PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 517 for (i = 0; i < connectSize; ++i) { 518 if (i < nodes[cs][0]) {/* Vertices */ 519 connect[i+off] = closure[(i+edgesInClosure+facesInClosure+1)*2] + 1; 520 connect[i+off] -= cellsNotInConnectivity; 521 } else if (i < nodes[cs][0]+nodes[cs][1]) { /* Edges */ 522 connect[i+off] = closure[(i-verticesInClosure+facesInClosure+1)*2] + 1; 523 if (nodes[cs][2] == 0) connect[i+off] -= numFaces; 524 connect[i+off] -= cellsNotInConnectivity; 525 } else if (i < nodes[cs][0]+nodes[cs][1]+nodes[cs][3]){ /* Cells */ 526 connect[i+off] = closure[0] + 1; 527 connect[i+off] -= skipCells; 528 } else if (i < nodes[cs][0]+nodes[cs][1]+nodes[cs][3]+nodes[cs][2]){ /* Faces */ 529 connect[i+off] = closure[(i-edgesInClosure-verticesInClosure)*2] + 1; 530 connect[i+off] -= cellsNotInConnectivity; 531 } else { 532 connect[i+off] = -1; 533 } 534 } 535 /* Tetrahedra are inverted */ 536 if (type[cs] == TET) { 537 temp = connect[0+off]; connect[0+off] = connect[1+off]; connect[1+off] = temp; 538 if (degree == 2) { 539 temp = connect[5+off]; connect[5+off] = connect[6+off]; connect[6+off] = temp; 540 temp = connect[7+off]; connect[7+off] = connect[8+off]; connect[8+off] = temp; 541 } 542 } 543 /* Hexahedra are inverted */ 544 if (type[cs] == HEX) { 545 temp = connect[1+off]; connect[1+off] = connect[3+off]; connect[3+off] = temp; 546 if (degree == 2) { 547 temp = connect[8+off]; connect[8+off] = connect[11+off]; connect[11+off] = temp; 548 temp = connect[9+off]; connect[9+off] = connect[10+off]; connect[10+off] = temp; 549 temp = connect[16+off]; connect[16+off] = connect[17+off]; connect[17+off] = temp; 550 temp = connect[18+off]; connect[18+off] = connect[19+off]; connect[19+off] = temp; 551 552 temp = connect[12+off]; connect[12+off] = connect[16+off]; connect[16+off] = temp; 553 temp = connect[13+off]; connect[13+off] = connect[17+off]; connect[17+off] = temp; 554 temp = connect[14+off]; connect[14+off] = connect[18+off]; connect[18+off] = temp; 555 temp = connect[15+off]; connect[15+off] = connect[19+off]; connect[19+off] = temp; 556 557 temp = connect[23+off]; connect[23+off] = connect[26+off]; connect[26+off] = temp; 558 temp = connect[24+off]; connect[24+off] = connect[25+off]; connect[25+off] = temp; 559 temp = connect[25+off]; connect[25+off] = connect[26+off]; connect[26+off] = temp; 560 } 561 } 562 off += connectSize; 563 ierr = DMPlexRestoreTransitiveClosure(dm, cells[c], PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 564 } 565 PetscStackCallStandard(ex_put_conn,(exo->exoid, EX_ELEM_BLOCK, csIdx[cs], connect, 0, 0)); 566 skipCells += (nodes[cs][3] == 0)*csSize; 567 ierr = PetscFree(connect);CHKERRQ(ierr); 568 ierr = ISRestoreIndices(stratumIS, &cells);CHKERRQ(ierr); 569 ierr = ISDestroy(&stratumIS);CHKERRQ(ierr); 570 } 571 ierr = PetscFree(type);CHKERRQ(ierr); 572 /* --- Coordinates --- */ 573 ierr = PetscSectionSetChart(coordSection, pStart, pEnd);CHKERRQ(ierr); 574 if (num_cs) { 575 for (d = 0; d < depth; ++d) { 576 ierr = DMPlexGetDepthStratum(dm, d, &pStart, &pEnd);CHKERRQ(ierr); 577 for (p = pStart; p < pEnd; ++p) { 578 ierr = PetscSectionSetDof(coordSection, p, nodes[0][d] > 0);CHKERRQ(ierr); 579 } 580 } 581 } 582 for (cs = 0; cs < num_cs; ++cs) { 583 IS stratumIS; 584 const PetscInt *cells; 585 PetscInt csSize, c; 586 587 ierr = DMLabelGetStratumIS(csLabel, csIdx[cs], &stratumIS);CHKERRQ(ierr); 588 ierr = ISGetIndices(stratumIS, &cells);CHKERRQ(ierr); 589 ierr = ISGetSize(stratumIS, &csSize);CHKERRQ(ierr); 590 for (c = 0; c < csSize; ++c) { 591 ierr = PetscSectionSetDof(coordSection, cells[c], nodes[cs][3] > 0);CHKERRQ(ierr); 592 } 593 ierr = ISRestoreIndices(stratumIS, &cells);CHKERRQ(ierr); 594 ierr = ISDestroy(&stratumIS);CHKERRQ(ierr); 595 } 596 if (num_cs > 0) { 597 ierr = ISRestoreIndices(csIS, &csIdx);CHKERRQ(ierr); 598 ierr = ISDestroy(&csIS);CHKERRQ(ierr); 599 } 600 ierr = PetscFree(nodes);CHKERRQ(ierr); 601 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 602 if (numNodes > 0) { 603 const char *coordNames[3] = {"x", "y", "z"}; 604 PetscScalar *coords, *closure; 605 PetscReal *cval; 606 PetscInt hasDof, n = 0; 607 608 /* There can't be more than 24 values in the closure of a point for the coord coordSection */ 609 ierr = PetscCalloc3(numNodes*3, &coords, dim, &cval, 24, &closure);CHKERRQ(ierr); 610 ierr = DMGetCoordinatesLocalNoncollective(dm, &coord);CHKERRQ(ierr); 611 ierr = DMPlexGetChart(dm, &pStart, &pEnd);CHKERRQ(ierr); 612 for (p = pStart; p < pEnd; ++p) { 613 ierr = PetscSectionGetDof(coordSection, p, &hasDof); 614 if (hasDof) { 615 PetscInt closureSize = 24, j; 616 617 ierr = DMPlexVecGetClosure(cdm, NULL, coord, p, &closureSize, &closure);CHKERRQ(ierr); 618 for (d = 0; d < dim; ++d) { 619 cval[d] = 0.0; 620 for (j = 0; j < closureSize/dim; j++) cval[d] += closure[j*dim+d]; 621 coords[d*numNodes+n] = cval[d] * dim / closureSize; 622 } 623 ++n; 624 } 625 } 626 PetscStackCallStandard(ex_put_coord,(exo->exoid, &coords[0*numNodes], &coords[1*numNodes], &coords[2*numNodes])); 627 ierr = PetscFree3(coords, cval, closure);CHKERRQ(ierr); 628 PetscStackCallStandard(ex_put_coord_names,(exo->exoid, (char **) coordNames)); 629 } 630 631 /* --- Node Sets/Vertex Sets --- */ 632 DMHasLabel(dm, "Vertex Sets", &hasLabel); 633 if (hasLabel) { 634 PetscInt i, vs, vsSize; 635 const PetscInt *vsIdx, *vertices; 636 PetscInt *nodeList; 637 IS vsIS, stratumIS; 638 DMLabel vsLabel; 639 ierr = DMGetLabel(dm, "Vertex Sets", &vsLabel);CHKERRQ(ierr); 640 ierr = DMLabelGetValueIS(vsLabel, &vsIS);CHKERRQ(ierr); 641 ierr = ISGetIndices(vsIS, &vsIdx);CHKERRQ(ierr); 642 for (vs=0; vs<num_vs; ++vs) { 643 ierr = DMLabelGetStratumIS(vsLabel, vsIdx[vs], &stratumIS);CHKERRQ(ierr); 644 ierr = ISGetIndices(stratumIS, &vertices);CHKERRQ(ierr); 645 ierr = ISGetSize(stratumIS, &vsSize);CHKERRQ(ierr); 646 ierr = PetscMalloc1(vsSize, &nodeList); 647 for (i=0; i<vsSize; ++i) { 648 nodeList[i] = vertices[i] - skipCells + 1; 649 } 650 PetscStackCallStandard(ex_put_set_param,(exo->exoid, EX_NODE_SET, vsIdx[vs], vsSize, 0)); 651 PetscStackCallStandard(ex_put_set,(exo->exoid, EX_NODE_SET, vsIdx[vs], nodeList, NULL)); 652 ierr = ISRestoreIndices(stratumIS, &vertices);CHKERRQ(ierr); 653 ierr = ISDestroy(&stratumIS);CHKERRQ(ierr); 654 ierr = PetscFree(nodeList);CHKERRQ(ierr); 655 } 656 ierr = ISRestoreIndices(vsIS, &vsIdx);CHKERRQ(ierr); 657 ierr = ISDestroy(&vsIS);CHKERRQ(ierr); 658 } 659 /* --- Side Sets/Face Sets --- */ 660 ierr = DMHasLabel(dm, "Face Sets", &hasLabel);CHKERRQ(ierr); 661 if (hasLabel) { 662 PetscInt i, j, fs, fsSize; 663 const PetscInt *fsIdx, *faces; 664 IS fsIS, stratumIS; 665 DMLabel fsLabel; 666 PetscInt numPoints, *points; 667 PetscInt elem_list_size = 0; 668 PetscInt *elem_list, *elem_ind, *side_list; 669 670 ierr = DMGetLabel(dm, "Face Sets", &fsLabel);CHKERRQ(ierr); 671 /* Compute size of Node List and Element List */ 672 ierr = DMLabelGetValueIS(fsLabel, &fsIS);CHKERRQ(ierr); 673 ierr = ISGetIndices(fsIS, &fsIdx);CHKERRQ(ierr); 674 for (fs=0; fs<num_fs; ++fs) { 675 ierr = DMLabelGetStratumIS(fsLabel, fsIdx[fs], &stratumIS);CHKERRQ(ierr); 676 ierr = ISGetSize(stratumIS, &fsSize);CHKERRQ(ierr); 677 elem_list_size += fsSize; 678 ierr = ISDestroy(&stratumIS);CHKERRQ(ierr); 679 } 680 ierr = PetscMalloc3(num_fs, &elem_ind, 681 elem_list_size, &elem_list, 682 elem_list_size, &side_list);CHKERRQ(ierr); 683 elem_ind[0] = 0; 684 for (fs=0; fs<num_fs; ++fs) { 685 ierr = DMLabelGetStratumIS(fsLabel, fsIdx[fs], &stratumIS);CHKERRQ(ierr); 686 ierr = ISGetIndices(stratumIS, &faces);CHKERRQ(ierr); 687 ierr = ISGetSize(stratumIS, &fsSize);CHKERRQ(ierr); 688 /* Set Parameters */ 689 PetscStackCallStandard(ex_put_set_param,(exo->exoid, EX_SIDE_SET, fsIdx[fs], fsSize, 0)); 690 /* Indices */ 691 if (fs<num_fs-1) { 692 elem_ind[fs+1] = elem_ind[fs] + fsSize; 693 } 694 695 for (i=0; i<fsSize; ++i) { 696 /* Element List */ 697 points = NULL; 698 ierr = DMPlexGetTransitiveClosure(dm, faces[i], PETSC_FALSE, &numPoints, &points);CHKERRQ(ierr); 699 elem_list[elem_ind[fs] + i] = points[2] +1; 700 ierr = DMPlexRestoreTransitiveClosure(dm, faces[i], PETSC_FALSE, &numPoints, &points);CHKERRQ(ierr); 701 702 /* Side List */ 703 points = NULL; 704 ierr = DMPlexGetTransitiveClosure(dm, elem_list[elem_ind[fs] + i]-1, PETSC_TRUE, &numPoints, &points);CHKERRQ(ierr); 705 for (j=1; j<numPoints; ++j) { 706 if (points[j*2]==faces[i]) {break;} 707 } 708 /* Convert HEX sides */ 709 if (numPoints == 27) { 710 if (j == 1) {j = 5;} 711 else if (j == 2) {j = 6;} 712 else if (j == 3) {j = 1;} 713 else if (j == 4) {j = 3;} 714 else if (j == 5) {j = 2;} 715 else if (j == 6) {j = 4;} 716 } 717 /* Convert TET sides */ 718 if (numPoints == 15) { 719 --j; 720 if (j == 0) {j = 4;} 721 } 722 side_list[elem_ind[fs] + i] = j; 723 ierr = DMPlexRestoreTransitiveClosure(dm, elem_list[elem_ind[fs] + i]-1, PETSC_TRUE, &numPoints, &points);CHKERRQ(ierr); 724 725 } 726 ierr = ISRestoreIndices(stratumIS, &faces);CHKERRQ(ierr); 727 ierr = ISDestroy(&stratumIS);CHKERRQ(ierr); 728 } 729 ierr = ISRestoreIndices(fsIS, &fsIdx);CHKERRQ(ierr); 730 ierr = ISDestroy(&fsIS);CHKERRQ(ierr); 731 732 /* Put side sets */ 733 for (fs=0; fs<num_fs; ++fs) { 734 PetscStackCallStandard(ex_put_set,(exo->exoid, EX_SIDE_SET, fsIdx[fs], &elem_list[elem_ind[fs]], &side_list[elem_ind[fs]])); 735 } 736 ierr = PetscFree3(elem_ind,elem_list,side_list);CHKERRQ(ierr); 737 } 738 /* 739 close the exodus file 740 */ 741 ex_close(exo->exoid); 742 exo->exoid = -1; 743 } 744 ierr = PetscSectionDestroy(&coordSection);CHKERRQ(ierr); 745 746 /* 747 reopen the file in parallel 748 */ 749 EXO_mode = EX_WRITE; 750 #if defined(PETSC_USE_64BIT_INDICES) 751 EXO_mode += EX_ALL_INT64_API; 752 #endif 753 CPU_word_size = sizeof(PetscReal); 754 IO_word_size = sizeof(PetscReal); 755 exo->exoid = ex_open_par(exo->filename,EXO_mode,&CPU_word_size,&IO_word_size,&EXO_version,comm,MPI_INFO_NULL); 756 if (exo->exoid < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_LIB, "ex_open_par failed for %s", exo->filename); 757 PetscFunctionReturn(0); 758 } 759 760 /* 761 VecView_PlexExodusII_Internal - Write a Vec corresponding in an exodus file 762 763 Collective on v 764 765 Input Parameters: 766 + v - The vector to be written 767 - viewer - The PetscViewerExodusII viewer associate to an exodus file 768 769 Notes: 770 The exodus result variable index is obtained by comparing the Vec name and the 771 names of variables declared in the exodus file. For instance for a Vec named "V" 772 the location in the exodus file will be the first match of "V", "V_X", "V_XX", "V_1", or "V_11" 773 amongst all variables. 774 In the event where a nodal and zonal variable both match, the function will return an error instead of 775 possibly corrupting the file 776 777 Level: beginner 778 779 .seealso: EXOGetVarIndex_Internal(),DMPlexView_ExodusII(),VecView_PlexExodusII() 780 @*/ 781 PetscErrorCode VecView_PlexExodusII_Internal(Vec v, PetscViewer viewer) 782 { 783 DM dm; 784 MPI_Comm comm; 785 PetscMPIInt rank; 786 787 788 int exoid,offsetN = 0, offsetZ = 0; 789 const char *vecname; 790 PetscInt step; 791 PetscErrorCode ierr; 792 793 PetscFunctionBegin; 794 ierr = PetscObjectGetComm((PetscObject) v, &comm);CHKERRQ(ierr); 795 ierr = MPI_Comm_rank(comm, &rank);CHKERRMPI(ierr); 796 ierr = PetscViewerExodusIIGetId(viewer,&exoid);CHKERRQ(ierr); 797 ierr = VecGetDM(v, &dm);CHKERRQ(ierr); 798 ierr = PetscObjectGetName((PetscObject) v, &vecname);CHKERRQ(ierr); 799 800 ierr = DMGetOutputSequenceNumber(dm,&step,NULL);CHKERRQ(ierr); 801 ierr = EXOGetVarIndex_Internal(exoid,EX_NODAL,vecname,&offsetN);CHKERRQ(ierr); 802 ierr = EXOGetVarIndex_Internal(exoid,EX_ELEM_BLOCK,vecname,&offsetZ);CHKERRQ(ierr); 803 if (offsetN <= 0 && offsetZ <= 0) SETERRQ1(comm, PETSC_ERR_FILE_UNEXPECTED, "Found both nodal and zonal variable %s in exodus file. \n", vecname); 804 if (offsetN > 0) { 805 ierr = VecViewPlex_ExodusII_Nodal_Internal(v,exoid,(int) step+1,offsetN);CHKERRQ(ierr); 806 } else if (offsetZ > 0) { 807 ierr = VecViewPlex_ExodusII_Zonal_Internal(v,exoid,(int) step+1,offsetZ);CHKERRQ(ierr); 808 } else SETERRQ1(comm, PETSC_ERR_FILE_UNEXPECTED, "Could not find nodal or zonal variable %s in exodus file. \n", vecname); 809 PetscFunctionReturn(0); 810 } 811 812 /* 813 VecLoad_PlexExodusII_Internal - Write a Vec corresponding in an exodus file 814 815 Collective on v 816 817 Input Parameters: 818 + v - The vector to be written 819 - viewer - The PetscViewerExodusII viewer associate to an exodus file 820 821 Notes: 822 The exodus result variable index is obtained by comparing the Vec name and the 823 names of variables declared in the exodus file. For instance for a Vec named "V" 824 the location in the exodus file will be the first match of "V", "V_X", "V_XX", "V_1", or "V_11" 825 amongst all variables. 826 In the event where a nodal and zonal variable both match, the function will return an error instead of 827 possibly corrupting the file 828 829 Level: beginner 830 831 .seealso: EXOGetVarIndex_Internal(),DMPlexView_ExodusII(),VecView_PlexExodusII() 832 @*/ 833 PetscErrorCode VecLoad_PlexExodusII_Internal(Vec v, PetscViewer viewer) 834 { 835 DM dm; 836 MPI_Comm comm; 837 PetscMPIInt rank; 838 839 840 int exoid,offsetN = 0, offsetZ = 0; 841 const char *vecname; 842 PetscInt step; 843 PetscErrorCode ierr; 844 845 PetscFunctionBegin; 846 ierr = PetscObjectGetComm((PetscObject) v, &comm);CHKERRQ(ierr); 847 ierr = MPI_Comm_rank(comm, &rank);CHKERRMPI(ierr); 848 ierr = PetscViewerExodusIIGetId(viewer,&exoid);CHKERRQ(ierr); 849 ierr = VecGetDM(v, &dm);CHKERRQ(ierr); 850 ierr = PetscObjectGetName((PetscObject) v, &vecname);CHKERRQ(ierr); 851 852 ierr = DMGetOutputSequenceNumber(dm,&step,NULL);CHKERRQ(ierr); 853 ierr = EXOGetVarIndex_Internal(exoid,EX_NODAL,vecname,&offsetN);CHKERRQ(ierr); 854 ierr = EXOGetVarIndex_Internal(exoid,EX_ELEM_BLOCK,vecname,&offsetZ);CHKERRQ(ierr); 855 if (offsetN <= 0 && offsetZ <= 0) SETERRQ1(comm, PETSC_ERR_FILE_UNEXPECTED, "Found both nodal and zonal variable %s in exodus file. \n", vecname); 856 if (offsetN > 0) { 857 ierr = VecLoadPlex_ExodusII_Nodal_Internal(v,exoid,(int) step+1,offsetN);CHKERRQ(ierr); 858 } else if (offsetZ > 0) { 859 ierr = VecLoadPlex_ExodusII_Zonal_Internal(v,exoid,(int) step+1,offsetZ);CHKERRQ(ierr); 860 } else SETERRQ1(comm, PETSC_ERR_FILE_UNEXPECTED, "Could not find nodal or zonal variable %s in exodus file. \n", vecname); 861 PetscFunctionReturn(0); 862 } 863 864 /* 865 VecViewPlex_ExodusII_Nodal_Internal - Write a Vec corresponding to a nodal field to an exodus file 866 867 Collective on v 868 869 Input Parameters: 870 + v - The vector to be written 871 . exoid - the exodus id of a file (obtained from ex_open or ex_create for instance) 872 . step - the time step to write at (exodus steps are numbered starting from 1) 873 - offset - the location of the variable in the file 874 875 Notes: 876 The exodus result nodal variable index is obtained by comparing the Vec name and the 877 names of nodal variables declared in the exodus file. For instance for a Vec named "V" 878 the location in the exodus file will be the first match of "V", "V_X", "V_XX", "V_1", or "V_11" 879 amongst all nodal variables. 880 881 Level: beginner 882 883 .seealso: EXOGetVarIndex_Internal(),DMPlexView_ExodusII_Internal(),VecLoadNodal_PlexEXO(),VecViewZonal_PlexEXO(),VecLoadZonal_PlexEXO() 884 @*/ 885 PetscErrorCode VecViewPlex_ExodusII_Nodal_Internal(Vec v, int exoid, int step, int offset) 886 { 887 MPI_Comm comm; 888 PetscMPIInt size; 889 DM dm; 890 Vec vNatural, vComp; 891 const PetscScalar *varray; 892 PetscInt xs, xe, bs; 893 PetscBool useNatural; 894 PetscErrorCode ierr; 895 896 PetscFunctionBegin; 897 ierr = PetscObjectGetComm((PetscObject) v, &comm);CHKERRQ(ierr); 898 ierr = MPI_Comm_size(comm, &size);CHKERRMPI(ierr); 899 ierr = VecGetDM(v, &dm);CHKERRQ(ierr); 900 ierr = DMGetUseNatural(dm, &useNatural);CHKERRQ(ierr); 901 useNatural = useNatural && size > 1 ? PETSC_TRUE : PETSC_FALSE; 902 if (useNatural) { 903 ierr = DMGetGlobalVector(dm, &vNatural);CHKERRQ(ierr); 904 ierr = DMPlexGlobalToNaturalBegin(dm, v, vNatural);CHKERRQ(ierr); 905 ierr = DMPlexGlobalToNaturalEnd(dm, v, vNatural);CHKERRQ(ierr); 906 } else { 907 vNatural = v; 908 } 909 910 /* Write local chunk of the result in the exodus file 911 exodus stores each component of a vector-valued field as a separate variable. 912 We assume that they are stored sequentially */ 913 ierr = VecGetOwnershipRange(vNatural, &xs, &xe);CHKERRQ(ierr); 914 ierr = VecGetBlockSize(vNatural, &bs);CHKERRQ(ierr); 915 if (bs == 1) { 916 ierr = VecGetArrayRead(vNatural, &varray);CHKERRQ(ierr); 917 PetscStackCallStandard(ex_put_partial_var,(exoid, step, EX_NODAL, offset, 1, xs+1, xe-xs, varray)); 918 ierr = VecRestoreArrayRead(vNatural, &varray);CHKERRQ(ierr); 919 } else { 920 IS compIS; 921 PetscInt c; 922 923 ierr = ISCreateStride(comm, (xe-xs)/bs, xs, bs, &compIS);CHKERRQ(ierr); 924 for (c = 0; c < bs; ++c) { 925 ierr = ISStrideSetStride(compIS, (xe-xs)/bs, xs+c, bs);CHKERRQ(ierr); 926 ierr = VecGetSubVector(vNatural, compIS, &vComp);CHKERRQ(ierr); 927 ierr = VecGetArrayRead(vComp, &varray);CHKERRQ(ierr); 928 PetscStackCallStandard(ex_put_partial_var,(exoid, step, EX_NODAL, offset+c, 1, xs/bs+1, (xe-xs)/bs, varray)); 929 ierr = VecRestoreArrayRead(vComp, &varray);CHKERRQ(ierr); 930 ierr = VecRestoreSubVector(vNatural, compIS, &vComp);CHKERRQ(ierr); 931 } 932 ierr = ISDestroy(&compIS);CHKERRQ(ierr); 933 } 934 if (useNatural) {ierr = DMRestoreGlobalVector(dm, &vNatural);CHKERRQ(ierr);} 935 PetscFunctionReturn(0); 936 } 937 938 /* 939 VecLoadPlex_ExodusII_Nodal_Internal - Read a Vec corresponding to a nodal field from an exodus file 940 941 Collective on v 942 943 Input Parameters: 944 + v - The vector to be written 945 . exoid - the exodus id of a file (obtained from ex_open or ex_create for instance) 946 . step - the time step to read at (exodus steps are numbered starting from 1) 947 - offset - the location of the variable in the file 948 949 Notes: 950 The exodus result nodal variable index is obtained by comparing the Vec name and the 951 names of nodal variables declared in the exodus file. For instance for a Vec named "V" 952 the location in the exodus file will be the first match of "V", "V_X", "V_XX", "V_1", or "V_11" 953 amongst all nodal variables. 954 955 Level: beginner 956 957 .seealso: EXOGetVarIndex_Internal(), DMPlexView_ExodusII_Internal(), VecViewPlex_ExodusII_Nodal_Internal(), VecViewZonal_PlexEXO(), VecLoadZonal_PlexEXO() 958 */ 959 PetscErrorCode VecLoadPlex_ExodusII_Nodal_Internal(Vec v, int exoid, int step, int offset) 960 { 961 MPI_Comm comm; 962 PetscMPIInt size; 963 DM dm; 964 Vec vNatural, vComp; 965 PetscScalar *varray; 966 PetscInt xs, xe, bs; 967 PetscBool useNatural; 968 PetscErrorCode ierr; 969 970 PetscFunctionBegin; 971 ierr = PetscObjectGetComm((PetscObject) v, &comm);CHKERRQ(ierr); 972 ierr = MPI_Comm_size(comm, &size);CHKERRMPI(ierr); 973 ierr = VecGetDM(v,&dm);CHKERRQ(ierr); 974 ierr = DMGetUseNatural(dm, &useNatural);CHKERRQ(ierr); 975 useNatural = useNatural && size > 1 ? PETSC_TRUE : PETSC_FALSE; 976 if (useNatural) {ierr = DMGetGlobalVector(dm,&vNatural);CHKERRQ(ierr);} 977 else {vNatural = v;} 978 979 /* Read local chunk from the file */ 980 ierr = VecGetOwnershipRange(vNatural, &xs, &xe);CHKERRQ(ierr); 981 ierr = VecGetBlockSize(vNatural, &bs);CHKERRQ(ierr); 982 if (bs == 1) { 983 ierr = VecGetArray(vNatural, &varray);CHKERRQ(ierr); 984 PetscStackCallStandard(ex_get_partial_var,(exoid, step, EX_NODAL, offset, 1, xs+1, xe-xs, varray)); 985 ierr = VecRestoreArray(vNatural, &varray);CHKERRQ(ierr); 986 } else { 987 IS compIS; 988 PetscInt c; 989 990 ierr = ISCreateStride(comm, (xe-xs)/bs, xs, bs, &compIS);CHKERRQ(ierr); 991 for (c = 0; c < bs; ++c) { 992 ierr = ISStrideSetStride(compIS, (xe-xs)/bs, xs+c, bs);CHKERRQ(ierr); 993 ierr = VecGetSubVector(vNatural, compIS, &vComp);CHKERRQ(ierr); 994 ierr = VecGetArray(vComp, &varray);CHKERRQ(ierr); 995 PetscStackCallStandard(ex_get_partial_var,(exoid, step, EX_NODAL, offset+c, 1, xs/bs+1, (xe-xs)/bs, varray)); 996 ierr = VecRestoreArray(vComp, &varray);CHKERRQ(ierr); 997 ierr = VecRestoreSubVector(vNatural, compIS, &vComp);CHKERRQ(ierr); 998 } 999 ierr = ISDestroy(&compIS);CHKERRQ(ierr); 1000 } 1001 if (useNatural) { 1002 ierr = DMPlexNaturalToGlobalBegin(dm, vNatural, v);CHKERRQ(ierr); 1003 ierr = DMPlexNaturalToGlobalEnd(dm, vNatural, v);CHKERRQ(ierr); 1004 ierr = DMRestoreGlobalVector(dm, &vNatural);CHKERRQ(ierr); 1005 } 1006 PetscFunctionReturn(0); 1007 } 1008 1009 /* 1010 VecViewPlex_ExodusII_Zonal_Internal - Write a Vec corresponding to a zonal (cell based) field to an exodus file 1011 1012 Collective on v 1013 1014 Input Parameters: 1015 + v - The vector to be written 1016 . exoid - the exodus id of a file (obtained from ex_open or ex_create for instance) 1017 . step - the time step to write at (exodus steps are numbered starting from 1) 1018 - offset - the location of the variable in the file 1019 1020 Notes: 1021 The exodus result zonal variable index is obtained by comparing the Vec name and the 1022 names of zonal variables declared in the exodus file. For instance for a Vec named "V" 1023 the location in the exodus file will be the first match of "V", "V_X", "V_XX", "V_1", or "V_11" 1024 amongst all zonal variables. 1025 1026 Level: beginner 1027 1028 .seealso: EXOGetVarIndex_Internal(),DMPlexView_ExodusII_Internal(),VecViewPlex_ExodusII_Nodal_Internal(),VecLoadPlex_ExodusII_Nodal_Internal(),VecLoadPlex_ExodusII_Zonal_Internal() 1029 */ 1030 PetscErrorCode VecViewPlex_ExodusII_Zonal_Internal(Vec v, int exoid, int step, int offset) 1031 { 1032 MPI_Comm comm; 1033 PetscMPIInt size; 1034 DM dm; 1035 Vec vNatural, vComp; 1036 const PetscScalar *varray; 1037 PetscInt xs, xe, bs; 1038 PetscBool useNatural; 1039 IS compIS; 1040 PetscInt *csSize, *csID; 1041 PetscInt numCS, set, csxs = 0; 1042 PetscErrorCode ierr; 1043 1044 PetscFunctionBegin; 1045 ierr = PetscObjectGetComm((PetscObject)v, &comm);CHKERRQ(ierr); 1046 ierr = MPI_Comm_size(comm, &size);CHKERRMPI(ierr); 1047 ierr = VecGetDM(v, &dm);CHKERRQ(ierr); 1048 ierr = DMGetUseNatural(dm, &useNatural);CHKERRQ(ierr); 1049 useNatural = useNatural && size > 1 ? PETSC_TRUE : PETSC_FALSE; 1050 if (useNatural) { 1051 ierr = DMGetGlobalVector(dm, &vNatural);CHKERRQ(ierr); 1052 ierr = DMPlexGlobalToNaturalBegin(dm, v, vNatural);CHKERRQ(ierr); 1053 ierr = DMPlexGlobalToNaturalEnd(dm, v, vNatural);CHKERRQ(ierr); 1054 } else { 1055 vNatural = v; 1056 } 1057 1058 /* Write local chunk of the result in the exodus file 1059 exodus stores each component of a vector-valued field as a separate variable. 1060 We assume that they are stored sequentially 1061 Zonal variables are accessed one element block at a time, so we loop through the cell sets, 1062 but once the vector has been reordered to natural size, we cannot use the label informations 1063 to figure out what to save where. */ 1064 numCS = ex_inquire_int(exoid, EX_INQ_ELEM_BLK); 1065 ierr = PetscMalloc2(numCS, &csID, numCS, &csSize);CHKERRQ(ierr); 1066 PetscStackCallStandard(ex_get_ids,(exoid, EX_ELEM_BLOCK, csID)); 1067 for (set = 0; set < numCS; ++set) { 1068 ex_block block; 1069 1070 block.id = csID[set]; 1071 block.type = EX_ELEM_BLOCK; 1072 PetscStackCallStandard(ex_get_block_param,(exoid, &block)); 1073 csSize[set] = block.num_entry; 1074 } 1075 ierr = VecGetOwnershipRange(vNatural, &xs, &xe);CHKERRQ(ierr); 1076 ierr = VecGetBlockSize(vNatural, &bs);CHKERRQ(ierr); 1077 if (bs > 1) {ierr = ISCreateStride(comm, (xe-xs)/bs, xs, bs, &compIS);CHKERRQ(ierr);} 1078 for (set = 0; set < numCS; set++) { 1079 PetscInt csLocalSize, c; 1080 1081 /* range of indices for set setID[set]: csxs:csxs + csSize[set]-1 1082 local slice of zonal values: xs/bs,xm/bs-1 1083 intersection: max(xs/bs,csxs),min(xm/bs-1,csxs + csSize[set]-1) */ 1084 csLocalSize = PetscMax(0, PetscMin(xe/bs, csxs+csSize[set]) - PetscMax(xs/bs, csxs)); 1085 if (bs == 1) { 1086 ierr = VecGetArrayRead(vNatural, &varray);CHKERRQ(ierr); 1087 PetscStackCallStandard(ex_put_partial_var,(exoid, step, EX_ELEM_BLOCK, offset, csID[set], PetscMax(xs-csxs, 0)+1, csLocalSize, &varray[PetscMax(0, csxs-xs)])); 1088 ierr = VecRestoreArrayRead(vNatural, &varray);CHKERRQ(ierr); 1089 } else { 1090 for (c = 0; c < bs; ++c) { 1091 ierr = ISStrideSetStride(compIS, (xe-xs)/bs, xs+c, bs);CHKERRQ(ierr); 1092 ierr = VecGetSubVector(vNatural, compIS, &vComp);CHKERRQ(ierr); 1093 ierr = VecGetArrayRead(vComp, &varray);CHKERRQ(ierr); 1094 PetscStackCallStandard(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)])); 1095 ierr = VecRestoreArrayRead(vComp, &varray);CHKERRQ(ierr); 1096 ierr = VecRestoreSubVector(vNatural, compIS, &vComp);CHKERRQ(ierr); 1097 } 1098 } 1099 csxs += csSize[set]; 1100 } 1101 ierr = PetscFree2(csID, csSize);CHKERRQ(ierr); 1102 if (bs > 1) {ierr = ISDestroy(&compIS);CHKERRQ(ierr);} 1103 if (useNatural) {ierr = DMRestoreGlobalVector(dm,&vNatural);CHKERRQ(ierr);} 1104 PetscFunctionReturn(0); 1105 } 1106 1107 /* 1108 VecLoadPlex_ExodusII_Zonal_Internal - Read a Vec corresponding to a zonal (cell based) field from an exodus file 1109 1110 Collective on v 1111 1112 Input Parameters: 1113 + v - The vector to be written 1114 . exoid - the exodus id of a file (obtained from ex_open or ex_create for instance) 1115 . step - the time step to read at (exodus steps are numbered starting from 1) 1116 - offset - the location of the variable in the file 1117 1118 Notes: 1119 The exodus result zonal variable index is obtained by comparing the Vec name and the 1120 names of zonal variables declared in the exodus file. For instance for a Vec named "V" 1121 the location in the exodus file will be the first match of "V", "V_X", "V_XX", "V_1", or "V_11" 1122 amongst all zonal variables. 1123 1124 Level: beginner 1125 1126 .seealso: EXOGetVarIndex_Internal(), DMPlexView_ExodusII_Internal(), VecViewPlex_ExodusII_Nodal_Internal(), VecLoadPlex_ExodusII_Nodal_Internal(), VecLoadPlex_ExodusII_Zonal_Internal() 1127 */ 1128 PetscErrorCode VecLoadPlex_ExodusII_Zonal_Internal(Vec v, int exoid, int step, int offset) 1129 { 1130 MPI_Comm comm; 1131 PetscMPIInt size; 1132 DM dm; 1133 Vec vNatural, vComp; 1134 PetscScalar *varray; 1135 PetscInt xs, xe, bs; 1136 PetscBool useNatural; 1137 IS compIS; 1138 PetscInt *csSize, *csID; 1139 PetscInt numCS, set, csxs = 0; 1140 PetscErrorCode ierr; 1141 1142 PetscFunctionBegin; 1143 ierr = PetscObjectGetComm((PetscObject)v,&comm);CHKERRQ(ierr); 1144 ierr = MPI_Comm_size(comm, &size);CHKERRMPI(ierr); 1145 ierr = VecGetDM(v, &dm);CHKERRQ(ierr); 1146 ierr = DMGetUseNatural(dm, &useNatural);CHKERRQ(ierr); 1147 useNatural = useNatural && size > 1 ? PETSC_TRUE : PETSC_FALSE; 1148 if (useNatural) {ierr = DMGetGlobalVector(dm,&vNatural);CHKERRQ(ierr);} 1149 else {vNatural = v;} 1150 1151 /* Read local chunk of the result in the exodus file 1152 exodus stores each component of a vector-valued field as a separate variable. 1153 We assume that they are stored sequentially 1154 Zonal variables are accessed one element block at a time, so we loop through the cell sets, 1155 but once the vector has been reordered to natural size, we cannot use the label informations 1156 to figure out what to save where. */ 1157 numCS = ex_inquire_int(exoid, EX_INQ_ELEM_BLK); 1158 ierr = PetscMalloc2(numCS, &csID, numCS, &csSize);CHKERRQ(ierr); 1159 PetscStackCallStandard(ex_get_ids,(exoid, EX_ELEM_BLOCK, csID)); 1160 for (set = 0; set < numCS; ++set) { 1161 ex_block block; 1162 1163 block.id = csID[set]; 1164 block.type = EX_ELEM_BLOCK; 1165 PetscStackCallStandard(ex_get_block_param,(exoid, &block)); 1166 csSize[set] = block.num_entry; 1167 } 1168 ierr = VecGetOwnershipRange(vNatural, &xs, &xe);CHKERRQ(ierr); 1169 ierr = VecGetBlockSize(vNatural, &bs);CHKERRQ(ierr); 1170 if (bs > 1) {ierr = ISCreateStride(comm, (xe-xs)/bs, xs, bs, &compIS);CHKERRQ(ierr);} 1171 for (set = 0; set < numCS; ++set) { 1172 PetscInt csLocalSize, c; 1173 1174 /* range of indices for set setID[set]: csxs:csxs + csSize[set]-1 1175 local slice of zonal values: xs/bs,xm/bs-1 1176 intersection: max(xs/bs,csxs),min(xm/bs-1,csxs + csSize[set]-1) */ 1177 csLocalSize = PetscMax(0, PetscMin(xe/bs, csxs+csSize[set]) - PetscMax(xs/bs, csxs)); 1178 if (bs == 1) { 1179 ierr = VecGetArray(vNatural, &varray);CHKERRQ(ierr); 1180 PetscStackCallStandard(ex_get_partial_var,(exoid, step, EX_ELEM_BLOCK, offset, csID[set], PetscMax(xs-csxs, 0)+1, csLocalSize, &varray[PetscMax(0, csxs-xs)])); 1181 ierr = VecRestoreArray(vNatural, &varray);CHKERRQ(ierr); 1182 } else { 1183 for (c = 0; c < bs; ++c) { 1184 ierr = ISStrideSetStride(compIS, (xe-xs)/bs, xs+c, bs);CHKERRQ(ierr); 1185 ierr = VecGetSubVector(vNatural, compIS, &vComp);CHKERRQ(ierr); 1186 ierr = VecGetArray(vComp, &varray);CHKERRQ(ierr); 1187 PetscStackCallStandard(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)])); 1188 ierr = VecRestoreArray(vComp, &varray);CHKERRQ(ierr); 1189 ierr = VecRestoreSubVector(vNatural, compIS, &vComp);CHKERRQ(ierr); 1190 } 1191 } 1192 csxs += csSize[set]; 1193 } 1194 ierr = PetscFree2(csID, csSize);CHKERRQ(ierr); 1195 if (bs > 1) {ierr = ISDestroy(&compIS);CHKERRQ(ierr);} 1196 if (useNatural) { 1197 ierr = DMPlexNaturalToGlobalBegin(dm, vNatural, v);CHKERRQ(ierr); 1198 ierr = DMPlexNaturalToGlobalEnd(dm, vNatural, v);CHKERRQ(ierr); 1199 ierr = DMRestoreGlobalVector(dm, &vNatural);CHKERRQ(ierr); 1200 } 1201 PetscFunctionReturn(0); 1202 } 1203 #endif 1204 1205 /*@ 1206 PetscViewerExodusIIGetId - Get the file id of the ExodusII file 1207 1208 Logically Collective on PetscViewer 1209 1210 Input Parameter: 1211 . viewer - the PetscViewer 1212 1213 Output Parameter: 1214 - exoid - The ExodusII file id 1215 1216 Level: intermediate 1217 1218 .seealso: PetscViewerFileSetMode(), PetscViewerCreate(), PetscViewerSetType(), PetscViewerBinaryOpen() 1219 @*/ 1220 PetscErrorCode PetscViewerExodusIIGetId(PetscViewer viewer, int *exoid) 1221 { 1222 PetscErrorCode ierr; 1223 1224 PetscFunctionBegin; 1225 PetscValidHeaderSpecific(viewer, PETSC_VIEWER_CLASSID, 1); 1226 ierr = PetscTryMethod(viewer, "PetscViewerGetId_C",(PetscViewer,int*),(viewer,exoid));CHKERRQ(ierr); 1227 PetscFunctionReturn(0); 1228 } 1229 1230 1231 1232 /*@ 1233 PetscViewerExodusIISetOrder - Set the elements order in the exodusII file. 1234 1235 Collective 1236 1237 Input Parameters: 1238 + viewer - the viewer 1239 - order - elements order 1240 1241 Output Parameter: 1242 1243 Level: beginner 1244 1245 Note: 1246 1247 1248 .seealso: PetscViewerExodusIIGetId(), PetscViewerExodusIIGetOrder(), PetscViewerExodusIISetOrder() 1249 @*/ 1250 PetscErrorCode PetscViewerExodusIISetOrder(PetscViewer viewer, PetscInt order) 1251 { 1252 PetscErrorCode ierr; 1253 1254 PetscFunctionBegin; 1255 PetscValidHeaderSpecific(viewer, PETSC_VIEWER_CLASSID, 1); 1256 ierr = PetscTryMethod(viewer, "PetscViewerSetOrder_C",(PetscViewer,PetscInt),(viewer,order));CHKERRQ(ierr); 1257 PetscFunctionReturn(0); 1258 } 1259 1260 /*@ 1261 PetscViewerExodusIIGetOrder - Get the elements order in the exodusII file. 1262 1263 Collective 1264 1265 Input Parameters: 1266 + viewer - the viewer 1267 - order - elements order 1268 1269 Output Parameter: 1270 1271 Level: beginner 1272 1273 Note: 1274 1275 1276 .seealso: PetscViewerExodusIIGetId(), PetscViewerExodusIIGetOrder(), PetscViewerExodusIISetOrder() 1277 @*/ 1278 PetscErrorCode PetscViewerExodusIIGetOrder(PetscViewer viewer, PetscInt *order) 1279 { 1280 PetscErrorCode ierr; 1281 1282 PetscFunctionBegin; 1283 PetscValidHeaderSpecific(viewer, PETSC_VIEWER_CLASSID, 1); 1284 ierr = PetscTryMethod(viewer, "PetscViewerGetOrder_C",(PetscViewer,PetscInt*),(viewer,order));CHKERRQ(ierr); 1285 PetscFunctionReturn(0); 1286 } 1287 1288 /*@C 1289 PetscViewerExodusIIOpen - Opens a file for ExodusII input/output. 1290 1291 Collective 1292 1293 Input Parameters: 1294 + comm - MPI communicator 1295 . name - name of file 1296 - type - type of file 1297 $ FILE_MODE_WRITE - create new file for binary output 1298 $ FILE_MODE_READ - open existing file for binary input 1299 $ FILE_MODE_APPEND - open existing file for binary output 1300 1301 Output Parameter: 1302 . exo - PetscViewer for Exodus II input/output to use with the specified file 1303 1304 Level: beginner 1305 1306 Note: 1307 This PetscViewer should be destroyed with PetscViewerDestroy(). 1308 1309 1310 .seealso: PetscViewerPushFormat(), PetscViewerDestroy(), 1311 DMLoad(), PetscFileMode, PetscViewer, PetscViewerSetType(), PetscViewerFileSetMode(), PetscViewerFileSetName() 1312 @*/ 1313 PetscErrorCode PetscViewerExodusIIOpen(MPI_Comm comm, const char name[], PetscFileMode type, PetscViewer *exo) 1314 { 1315 PetscErrorCode ierr; 1316 1317 PetscFunctionBegin; 1318 ierr = PetscViewerCreate(comm, exo);CHKERRQ(ierr); 1319 ierr = PetscViewerSetType(*exo, PETSCVIEWEREXODUSII);CHKERRQ(ierr); 1320 ierr = PetscViewerFileSetMode(*exo, type);CHKERRQ(ierr); 1321 ierr = PetscViewerFileSetName(*exo, name);CHKERRQ(ierr); 1322 ierr = PetscViewerSetFromOptions(*exo);CHKERRQ(ierr); 1323 PetscFunctionReturn(0); 1324 } 1325 1326 /*@C 1327 DMPlexCreateExodusFromFile - Create a DMPlex mesh from an ExodusII file. 1328 1329 Collective 1330 1331 Input Parameters: 1332 + comm - The MPI communicator 1333 . filename - The name of the ExodusII file 1334 - interpolate - Create faces and edges in the mesh 1335 1336 Output Parameter: 1337 . dm - The DM object representing the mesh 1338 1339 Level: beginner 1340 1341 .seealso: DMPLEX, DMCreate(), DMPlexCreateExodus() 1342 @*/ 1343 PetscErrorCode DMPlexCreateExodusFromFile(MPI_Comm comm, const char filename[], PetscBool interpolate, DM *dm) 1344 { 1345 PetscMPIInt rank; 1346 PetscErrorCode ierr; 1347 #if defined(PETSC_HAVE_EXODUSII) 1348 int CPU_word_size = sizeof(PetscReal), IO_word_size = 0, exoid = -1; 1349 float version; 1350 #endif 1351 1352 PetscFunctionBegin; 1353 PetscValidCharPointer(filename, 2); 1354 ierr = MPI_Comm_rank(comm, &rank);CHKERRMPI(ierr); 1355 #if defined(PETSC_HAVE_EXODUSII) 1356 if (!rank) { 1357 exoid = ex_open(filename, EX_READ, &CPU_word_size, &IO_word_size, &version); 1358 if (exoid <= 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_LIB, "ex_open(\"%s\",...) did not return a valid file ID", filename); 1359 } 1360 ierr = DMPlexCreateExodus(comm, exoid, interpolate, dm);CHKERRQ(ierr); 1361 if (!rank) {PetscStackCallStandard(ex_close,(exoid));} 1362 PetscFunctionReturn(0); 1363 #else 1364 SETERRQ(comm, PETSC_ERR_SUP, "This method requires ExodusII support. Reconfigure using --download-exodusii"); 1365 #endif 1366 } 1367 1368 #if defined(PETSC_HAVE_EXODUSII) 1369 static PetscErrorCode ExodusGetCellType_Internal(const char *elem_type, DMPolytopeType *ct) 1370 { 1371 PetscBool flg; 1372 PetscErrorCode ierr; 1373 1374 PetscFunctionBegin; 1375 *ct = DM_POLYTOPE_UNKNOWN; 1376 ierr = PetscStrcmp(elem_type, "TRI", &flg);CHKERRQ(ierr); 1377 if (flg) {*ct = DM_POLYTOPE_TRIANGLE; goto done;} 1378 ierr = PetscStrcmp(elem_type, "TRI3", &flg);CHKERRQ(ierr); 1379 if (flg) {*ct = DM_POLYTOPE_TRIANGLE; goto done;} 1380 ierr = PetscStrcmp(elem_type, "QUAD", &flg);CHKERRQ(ierr); 1381 if (flg) {*ct = DM_POLYTOPE_QUADRILATERAL; goto done;} 1382 ierr = PetscStrcmp(elem_type, "QUAD4", &flg);CHKERRQ(ierr); 1383 if (flg) {*ct = DM_POLYTOPE_QUADRILATERAL; goto done;} 1384 ierr = PetscStrcmp(elem_type, "SHELL4", &flg);CHKERRQ(ierr); 1385 if (flg) {*ct = DM_POLYTOPE_QUADRILATERAL; goto done;} 1386 ierr = PetscStrcmp(elem_type, "TETRA", &flg);CHKERRQ(ierr); 1387 if (flg) {*ct = DM_POLYTOPE_TETRAHEDRON; goto done;} 1388 ierr = PetscStrcmp(elem_type, "TET4", &flg);CHKERRQ(ierr); 1389 if (flg) {*ct = DM_POLYTOPE_TETRAHEDRON; goto done;} 1390 ierr = PetscStrcmp(elem_type, "WEDGE", &flg);CHKERRQ(ierr); 1391 if (flg) {*ct = DM_POLYTOPE_TRI_PRISM; goto done;} 1392 ierr = PetscStrcmp(elem_type, "HEX", &flg);CHKERRQ(ierr); 1393 if (flg) {*ct = DM_POLYTOPE_HEXAHEDRON; goto done;} 1394 ierr = PetscStrcmp(elem_type, "HEX8", &flg);CHKERRQ(ierr); 1395 if (flg) {*ct = DM_POLYTOPE_HEXAHEDRON; goto done;} 1396 ierr = PetscStrcmp(elem_type, "HEXAHEDRON", &flg);CHKERRQ(ierr); 1397 if (flg) {*ct = DM_POLYTOPE_HEXAHEDRON; goto done;} 1398 SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unrecognized element type %s", elem_type); 1399 done: 1400 PetscFunctionReturn(0); 1401 } 1402 #endif 1403 1404 /*@ 1405 DMPlexCreateExodus - Create a DMPlex mesh from an ExodusII file ID. 1406 1407 Collective 1408 1409 Input Parameters: 1410 + comm - The MPI communicator 1411 . exoid - The ExodusII id associated with a exodus file and obtained using ex_open 1412 - interpolate - Create faces and edges in the mesh 1413 1414 Output Parameter: 1415 . dm - The DM object representing the mesh 1416 1417 Level: beginner 1418 1419 .seealso: DMPLEX, DMCreate() 1420 @*/ 1421 PetscErrorCode DMPlexCreateExodus(MPI_Comm comm, PetscInt exoid, PetscBool interpolate, DM *dm) 1422 { 1423 #if defined(PETSC_HAVE_EXODUSII) 1424 PetscMPIInt num_proc, rank; 1425 PetscSection coordSection; 1426 Vec coordinates; 1427 PetscScalar *coords; 1428 PetscInt coordSize, v; 1429 PetscErrorCode ierr; 1430 /* Read from ex_get_init() */ 1431 char title[PETSC_MAX_PATH_LEN+1]; 1432 int dim = 0, dimEmbed = 0, numVertices = 0, numCells = 0, numHybridCells = 0; 1433 int num_cs = 0, num_vs = 0, num_fs = 0; 1434 #endif 1435 1436 PetscFunctionBegin; 1437 #if defined(PETSC_HAVE_EXODUSII) 1438 ierr = MPI_Comm_rank(comm, &rank);CHKERRMPI(ierr); 1439 ierr = MPI_Comm_size(comm, &num_proc);CHKERRMPI(ierr); 1440 ierr = DMCreate(comm, dm);CHKERRQ(ierr); 1441 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr); 1442 /* Open EXODUS II file and read basic informations on rank 0, then broadcast to all processors */ 1443 if (!rank) { 1444 ierr = PetscMemzero(title,PETSC_MAX_PATH_LEN+1);CHKERRQ(ierr); 1445 PetscStackCallStandard(ex_get_init,(exoid, title, &dimEmbed, &numVertices, &numCells, &num_cs, &num_vs, &num_fs)); 1446 if (!num_cs) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Exodus file does not contain any cell set\n"); 1447 } 1448 ierr = MPI_Bcast(title, PETSC_MAX_PATH_LEN+1, MPI_CHAR, 0, comm);CHKERRMPI(ierr); 1449 ierr = MPI_Bcast(&dim, 1, MPI_INT, 0, comm);CHKERRMPI(ierr); 1450 ierr = PetscObjectSetName((PetscObject) *dm, title);CHKERRQ(ierr); 1451 ierr = DMPlexSetChart(*dm, 0, numCells+numVertices);CHKERRQ(ierr); 1452 /* We do not want this label automatically computed, instead we compute it here */ 1453 ierr = DMCreateLabel(*dm, "celltype");CHKERRQ(ierr); 1454 1455 /* Read cell sets information */ 1456 if (!rank) { 1457 PetscInt *cone; 1458 int c, cs, ncs, c_loc, v, v_loc; 1459 /* Read from ex_get_elem_blk_ids() */ 1460 int *cs_id, *cs_order; 1461 /* Read from ex_get_elem_block() */ 1462 char buffer[PETSC_MAX_PATH_LEN+1]; 1463 int num_cell_in_set, num_vertex_per_cell, num_hybrid, num_attr; 1464 /* Read from ex_get_elem_conn() */ 1465 int *cs_connect; 1466 1467 /* Get cell sets IDs */ 1468 ierr = PetscMalloc2(num_cs, &cs_id, num_cs, &cs_order);CHKERRQ(ierr); 1469 PetscStackCallStandard(ex_get_ids,(exoid, EX_ELEM_BLOCK, cs_id)); 1470 /* Read the cell set connectivity table and build mesh topology 1471 EXO standard requires that cells in cell sets be numbered sequentially and be pairwise disjoint. */ 1472 /* Check for a hybrid mesh */ 1473 for (cs = 0, num_hybrid = 0; cs < num_cs; ++cs) { 1474 DMPolytopeType ct; 1475 char elem_type[PETSC_MAX_PATH_LEN]; 1476 1477 ierr = PetscArrayzero(elem_type, sizeof(elem_type));CHKERRQ(ierr); 1478 PetscStackCallStandard(ex_get_elem_type,(exoid, cs_id[cs], elem_type)); 1479 ierr = ExodusGetCellType_Internal(elem_type, &ct);CHKERRQ(ierr); 1480 dim = PetscMax(dim, DMPolytopeTypeGetDim(ct)); 1481 PetscStackCallStandard(ex_get_block,(exoid, EX_ELEM_BLOCK, cs_id[cs], buffer, &num_cell_in_set,&num_vertex_per_cell, 0, 0, &num_attr)); 1482 switch (ct) { 1483 case DM_POLYTOPE_TRI_PRISM: 1484 cs_order[cs] = cs; 1485 numHybridCells += num_cell_in_set; 1486 ++num_hybrid; 1487 break; 1488 default: 1489 for (c = cs; c > cs-num_hybrid; --c) cs_order[c] = cs_order[c-1]; 1490 cs_order[cs-num_hybrid] = cs; 1491 } 1492 } 1493 /* First set sizes */ 1494 for (ncs = 0, c = 0; ncs < num_cs; ++ncs) { 1495 DMPolytopeType ct; 1496 char elem_type[PETSC_MAX_PATH_LEN]; 1497 const PetscInt cs = cs_order[ncs]; 1498 1499 ierr = PetscArrayzero(elem_type, sizeof(elem_type));CHKERRQ(ierr); 1500 PetscStackCallStandard(ex_get_elem_type,(exoid, cs_id[cs], elem_type)); 1501 ierr = ExodusGetCellType_Internal(elem_type, &ct);CHKERRQ(ierr); 1502 PetscStackCallStandard(ex_get_block,(exoid, EX_ELEM_BLOCK, cs_id[cs], buffer, &num_cell_in_set,&num_vertex_per_cell, 0, 0, &num_attr)); 1503 for (c_loc = 0; c_loc < num_cell_in_set; ++c_loc, ++c) { 1504 ierr = DMPlexSetConeSize(*dm, c, num_vertex_per_cell);CHKERRQ(ierr); 1505 ierr = DMPlexSetCellType(*dm, c, ct);CHKERRQ(ierr); 1506 } 1507 } 1508 for (v = numCells; v < numCells+numVertices; ++v) {ierr = DMPlexSetCellType(*dm, v, DM_POLYTOPE_POINT);CHKERRQ(ierr);} 1509 ierr = DMSetUp(*dm);CHKERRQ(ierr); 1510 for (ncs = 0, c = 0; ncs < num_cs; ++ncs) { 1511 const PetscInt cs = cs_order[ncs]; 1512 PetscStackCallStandard(ex_get_block,(exoid, EX_ELEM_BLOCK, cs_id[cs], buffer, &num_cell_in_set, &num_vertex_per_cell, 0, 0, &num_attr)); 1513 ierr = PetscMalloc2(num_vertex_per_cell*num_cell_in_set,&cs_connect,num_vertex_per_cell,&cone);CHKERRQ(ierr); 1514 PetscStackCallStandard(ex_get_conn,(exoid, EX_ELEM_BLOCK, cs_id[cs], cs_connect,NULL,NULL)); 1515 /* EXO uses Fortran-based indexing, DMPlex uses C-style and numbers cell first then vertices. */ 1516 for (c_loc = 0, v = 0; c_loc < num_cell_in_set; ++c_loc, ++c) { 1517 DMPolytopeType ct; 1518 1519 for (v_loc = 0; v_loc < num_vertex_per_cell; ++v_loc, ++v) { 1520 cone[v_loc] = cs_connect[v]+numCells-1; 1521 } 1522 ierr = DMPlexGetCellType(*dm, c, &ct);CHKERRQ(ierr); 1523 ierr = DMPlexInvertCell(ct, cone);CHKERRQ(ierr); 1524 ierr = DMPlexSetCone(*dm, c, cone);CHKERRQ(ierr); 1525 ierr = DMSetLabelValue(*dm, "Cell Sets", c, cs_id[cs]);CHKERRQ(ierr); 1526 } 1527 ierr = PetscFree2(cs_connect,cone);CHKERRQ(ierr); 1528 } 1529 ierr = PetscFree2(cs_id, cs_order);CHKERRQ(ierr); 1530 } 1531 { 1532 PetscInt ints[] = {dim, dimEmbed}; 1533 1534 ierr = MPI_Bcast(ints, 2, MPIU_INT, 0, comm);CHKERRMPI(ierr); 1535 ierr = DMSetDimension(*dm, ints[0]);CHKERRQ(ierr); 1536 ierr = DMSetCoordinateDim(*dm, ints[1]);CHKERRQ(ierr); 1537 dim = ints[0]; 1538 dimEmbed = ints[1]; 1539 } 1540 ierr = DMPlexSymmetrize(*dm);CHKERRQ(ierr); 1541 ierr = DMPlexStratify(*dm);CHKERRQ(ierr); 1542 if (interpolate) { 1543 DM idm; 1544 1545 ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr); 1546 ierr = DMDestroy(dm);CHKERRQ(ierr); 1547 *dm = idm; 1548 } 1549 1550 /* Create vertex set label */ 1551 if (!rank && (num_vs > 0)) { 1552 int vs, v; 1553 /* Read from ex_get_node_set_ids() */ 1554 int *vs_id; 1555 /* Read from ex_get_node_set_param() */ 1556 int num_vertex_in_set; 1557 /* Read from ex_get_node_set() */ 1558 int *vs_vertex_list; 1559 1560 /* Get vertex set ids */ 1561 ierr = PetscMalloc1(num_vs, &vs_id);CHKERRQ(ierr); 1562 PetscStackCallStandard(ex_get_ids,(exoid, EX_NODE_SET, vs_id)); 1563 for (vs = 0; vs < num_vs; ++vs) { 1564 PetscStackCallStandard(ex_get_set_param,(exoid, EX_NODE_SET, vs_id[vs], &num_vertex_in_set, NULL)); 1565 ierr = PetscMalloc1(num_vertex_in_set, &vs_vertex_list);CHKERRQ(ierr); 1566 PetscStackCallStandard(ex_get_set,(exoid, EX_NODE_SET, vs_id[vs], vs_vertex_list, NULL)); 1567 for (v = 0; v < num_vertex_in_set; ++v) { 1568 ierr = DMSetLabelValue(*dm, "Vertex Sets", vs_vertex_list[v]+numCells-1, vs_id[vs]);CHKERRQ(ierr); 1569 } 1570 ierr = PetscFree(vs_vertex_list);CHKERRQ(ierr); 1571 } 1572 ierr = PetscFree(vs_id);CHKERRQ(ierr); 1573 } 1574 /* Read coordinates */ 1575 ierr = DMGetCoordinateSection(*dm, &coordSection);CHKERRQ(ierr); 1576 ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr); 1577 ierr = PetscSectionSetFieldComponents(coordSection, 0, dimEmbed);CHKERRQ(ierr); 1578 ierr = PetscSectionSetChart(coordSection, numCells, numCells + numVertices);CHKERRQ(ierr); 1579 for (v = numCells; v < numCells+numVertices; ++v) { 1580 ierr = PetscSectionSetDof(coordSection, v, dimEmbed);CHKERRQ(ierr); 1581 ierr = PetscSectionSetFieldDof(coordSection, v, 0, dimEmbed);CHKERRQ(ierr); 1582 } 1583 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 1584 ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr); 1585 ierr = VecCreate(PETSC_COMM_SELF, &coordinates);CHKERRQ(ierr); 1586 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr); 1587 ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 1588 ierr = VecSetBlockSize(coordinates, dimEmbed);CHKERRQ(ierr); 1589 ierr = VecSetType(coordinates,VECSTANDARD);CHKERRQ(ierr); 1590 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 1591 if (!rank) { 1592 PetscReal *x, *y, *z; 1593 1594 ierr = PetscMalloc3(numVertices,&x,numVertices,&y,numVertices,&z);CHKERRQ(ierr); 1595 PetscStackCallStandard(ex_get_coord,(exoid, x, y, z)); 1596 if (dimEmbed > 0) { 1597 for (v = 0; v < numVertices; ++v) coords[v*dimEmbed+0] = x[v]; 1598 } 1599 if (dimEmbed > 1) { 1600 for (v = 0; v < numVertices; ++v) coords[v*dimEmbed+1] = y[v]; 1601 } 1602 if (dimEmbed > 2) { 1603 for (v = 0; v < numVertices; ++v) coords[v*dimEmbed+2] = z[v]; 1604 } 1605 ierr = PetscFree3(x,y,z);CHKERRQ(ierr); 1606 } 1607 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 1608 ierr = DMSetCoordinatesLocal(*dm, coordinates);CHKERRQ(ierr); 1609 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 1610 1611 /* Create side set label */ 1612 if (!rank && interpolate && (num_fs > 0)) { 1613 int fs, f, voff; 1614 /* Read from ex_get_side_set_ids() */ 1615 int *fs_id; 1616 /* Read from ex_get_side_set_param() */ 1617 int num_side_in_set; 1618 /* Read from ex_get_side_set_node_list() */ 1619 int *fs_vertex_count_list, *fs_vertex_list; 1620 /* Read side set labels */ 1621 char fs_name[MAX_STR_LENGTH+1]; 1622 1623 /* Get side set ids */ 1624 ierr = PetscMalloc1(num_fs, &fs_id);CHKERRQ(ierr); 1625 PetscStackCallStandard(ex_get_ids,(exoid, EX_SIDE_SET, fs_id)); 1626 for (fs = 0; fs < num_fs; ++fs) { 1627 PetscStackCallStandard(ex_get_set_param,(exoid, EX_SIDE_SET, fs_id[fs], &num_side_in_set, NULL)); 1628 ierr = PetscMalloc2(num_side_in_set,&fs_vertex_count_list,num_side_in_set*4,&fs_vertex_list);CHKERRQ(ierr); 1629 PetscStackCallStandard(ex_get_side_set_node_list,(exoid, fs_id[fs], fs_vertex_count_list, fs_vertex_list)); 1630 /* Get the specific name associated with this side set ID. */ 1631 int fs_name_err = ex_get_name(exoid, EX_SIDE_SET, fs_id[fs], fs_name); 1632 for (f = 0, voff = 0; f < num_side_in_set; ++f) { 1633 const PetscInt *faces = NULL; 1634 PetscInt faceSize = fs_vertex_count_list[f], numFaces; 1635 PetscInt faceVertices[4], v; 1636 1637 if (faceSize > 4) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "ExodusII side cannot have %d > 4 vertices", faceSize); 1638 for (v = 0; v < faceSize; ++v, ++voff) { 1639 faceVertices[v] = fs_vertex_list[voff]+numCells-1; 1640 } 1641 ierr = DMPlexGetFullJoin(*dm, faceSize, faceVertices, &numFaces, &faces);CHKERRQ(ierr); 1642 if (numFaces != 1) SETERRQ3(comm, PETSC_ERR_ARG_WRONG, "Invalid ExodusII side %d in set %d maps to %d faces", f, fs, numFaces); 1643 ierr = DMSetLabelValue(*dm, "Face Sets", faces[0], fs_id[fs]);CHKERRQ(ierr); 1644 /* Only add the label if one has been detected for this side set. */ 1645 if (!fs_name_err) { 1646 ierr = DMSetLabelValue(*dm, fs_name, faces[0], fs_id[fs]);CHKERRQ(ierr); 1647 } 1648 ierr = DMPlexRestoreJoin(*dm, faceSize, faceVertices, &numFaces, &faces);CHKERRQ(ierr); 1649 } 1650 ierr = PetscFree2(fs_vertex_count_list,fs_vertex_list);CHKERRQ(ierr); 1651 } 1652 ierr = PetscFree(fs_id);CHKERRQ(ierr); 1653 } 1654 PetscFunctionReturn(0); 1655 #else 1656 SETERRQ(comm, PETSC_ERR_SUP, "This method requires ExodusII support. Reconfigure using --download-exodusii"); 1657 #endif 1658 } 1659