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