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);CHKERRQ(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);CHKERRQ(ierr); 354 ierr = MPI_Comm_size(comm, &size);CHKERRQ(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);CHKERRQ(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(PETSC_COMM_SELF, (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);CHKERRQ(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(PETSC_COMM_SELF, (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);CHKERRQ(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);CHKERRQ(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);CHKERRQ(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 #else 1129 SETERRQ(comm, PETSC_ERR_SUP, "This method requires ExodusII support. Reconfigure using --download-exodusii"); 1130 #endif 1131 PetscFunctionReturn(0); 1132 } 1133 1134 /*@ 1135 DMPlexCreateExodus - Create a DMPlex mesh from an ExodusII file ID. 1136 1137 Collective 1138 1139 Input Parameters: 1140 + comm - The MPI communicator 1141 . exoid - The ExodusII id associated with a exodus file and obtained using ex_open 1142 - interpolate - Create faces and edges in the mesh 1143 1144 Output Parameter: 1145 . dm - The DM object representing the mesh 1146 1147 Level: beginner 1148 1149 .seealso: DMPLEX, DMCreate() 1150 @*/ 1151 PetscErrorCode DMPlexCreateExodus(MPI_Comm comm, PetscInt exoid, PetscBool interpolate, DM *dm) 1152 { 1153 #if defined(PETSC_HAVE_EXODUSII) 1154 PetscMPIInt num_proc, rank; 1155 PetscSection coordSection; 1156 Vec coordinates; 1157 PetscScalar *coords; 1158 PetscInt coordSize, v; 1159 PetscErrorCode ierr; 1160 /* Read from ex_get_init() */ 1161 char title[PETSC_MAX_PATH_LEN+1]; 1162 int dim = 0, numVertices = 0, numCells = 0, numHybridCells = 0; 1163 int num_cs = 0, num_vs = 0, num_fs = 0; 1164 #endif 1165 1166 PetscFunctionBegin; 1167 #if defined(PETSC_HAVE_EXODUSII) 1168 ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr); 1169 ierr = MPI_Comm_size(comm, &num_proc);CHKERRQ(ierr); 1170 ierr = DMCreate(comm, dm);CHKERRQ(ierr); 1171 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr); 1172 /* Open EXODUS II file and read basic informations on rank 0, then broadcast to all processors */ 1173 if (!rank) { 1174 ierr = PetscMemzero(title,PETSC_MAX_PATH_LEN+1);CHKERRQ(ierr); 1175 PetscStackCallStandard(ex_get_init,(exoid, title, &dim, &numVertices, &numCells, &num_cs, &num_vs, &num_fs)); 1176 if (!num_cs) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Exodus file does not contain any cell set\n"); 1177 } 1178 ierr = MPI_Bcast(title, PETSC_MAX_PATH_LEN+1, MPI_CHAR, 0, comm);CHKERRQ(ierr); 1179 ierr = MPI_Bcast(&dim, 1, MPI_INT, 0, comm);CHKERRQ(ierr); 1180 ierr = PetscObjectSetName((PetscObject) *dm, title);CHKERRQ(ierr); 1181 ierr = DMSetDimension(*dm, dim);CHKERRQ(ierr); 1182 ierr = DMPlexSetChart(*dm, 0, numCells+numVertices);CHKERRQ(ierr); 1183 /* We do not want this label automatically computed, instead we compute it here */ 1184 ierr = DMCreateLabel(*dm, "celltype");CHKERRQ(ierr); 1185 1186 /* Read cell sets information */ 1187 if (!rank) { 1188 PetscInt *cone; 1189 int c, cs, ncs, c_loc, v, v_loc; 1190 /* Read from ex_get_elem_blk_ids() */ 1191 int *cs_id, *cs_order; 1192 /* Read from ex_get_elem_block() */ 1193 char buffer[PETSC_MAX_PATH_LEN+1]; 1194 int num_cell_in_set, num_vertex_per_cell, num_hybrid, num_attr; 1195 /* Read from ex_get_elem_conn() */ 1196 int *cs_connect; 1197 1198 /* Get cell sets IDs */ 1199 ierr = PetscMalloc2(num_cs, &cs_id, num_cs, &cs_order);CHKERRQ(ierr); 1200 PetscStackCallStandard(ex_get_ids,(exoid, EX_ELEM_BLOCK, cs_id)); 1201 /* Read the cell set connectivity table and build mesh topology 1202 EXO standard requires that cells in cell sets be numbered sequentially and be pairwise disjoint. */ 1203 /* Check for a hybrid mesh */ 1204 for (cs = 0, num_hybrid = 0; cs < num_cs; ++cs) { 1205 PetscStackCallStandard(ex_get_block,(exoid, EX_ELEM_BLOCK, cs_id[cs], buffer, &num_cell_in_set,&num_vertex_per_cell, 0, 0, &num_attr)); 1206 switch (dim) { 1207 case 3: 1208 switch (num_vertex_per_cell) { 1209 case 6: 1210 cs_order[cs] = cs; 1211 numHybridCells += num_cell_in_set; 1212 ++num_hybrid; 1213 break; 1214 default: 1215 for (c = cs; c > cs-num_hybrid; --c) cs_order[c] = cs_order[c-1]; 1216 cs_order[cs-num_hybrid] = cs; 1217 } 1218 break; 1219 default: 1220 for (c = cs; c > cs-num_hybrid; --c) cs_order[c] = cs_order[c-1]; 1221 cs_order[cs-num_hybrid] = cs; 1222 } 1223 } 1224 /* First set sizes */ 1225 for (ncs = 0, c = 0; ncs < num_cs; ++ncs) { 1226 const PetscInt cs = cs_order[ncs]; 1227 PetscStackCallStandard(ex_get_block,(exoid, EX_ELEM_BLOCK, cs_id[cs], buffer, &num_cell_in_set,&num_vertex_per_cell, 0, 0, &num_attr)); 1228 for (c_loc = 0; c_loc < num_cell_in_set; ++c_loc, ++c) { 1229 ierr = DMPlexSetConeSize(*dm, c, num_vertex_per_cell);CHKERRQ(ierr); 1230 if (c >= numCells-numHybridCells) { 1231 ierr = DMPlexSetCellType(*dm, c, DM_POLYTOPE_TRI_PRISM);CHKERRQ(ierr); 1232 } else { 1233 DMPolytopeType ct; 1234 1235 ierr = DMPlexComputeCellType_Internal(*dm, c, 1, &ct);CHKERRQ(ierr); 1236 ierr = DMPlexSetCellType(*dm, c, ct);CHKERRQ(ierr); 1237 } 1238 } 1239 } 1240 for (v = numCells; v < numCells+numVertices; ++v) {ierr = DMPlexSetCellType(*dm, v, DM_POLYTOPE_POINT);CHKERRQ(ierr);} 1241 ierr = DMSetUp(*dm);CHKERRQ(ierr); 1242 for (ncs = 0, c = 0; ncs < num_cs; ++ncs) { 1243 const PetscInt cs = cs_order[ncs]; 1244 PetscStackCallStandard(ex_get_block,(exoid, EX_ELEM_BLOCK, cs_id[cs], buffer, &num_cell_in_set, &num_vertex_per_cell, 0, 0, &num_attr)); 1245 ierr = PetscMalloc2(num_vertex_per_cell*num_cell_in_set,&cs_connect,num_vertex_per_cell,&cone);CHKERRQ(ierr); 1246 PetscStackCallStandard(ex_get_conn,(exoid, EX_ELEM_BLOCK, cs_id[cs], cs_connect,NULL,NULL)); 1247 /* EXO uses Fortran-based indexing, DMPlex uses C-style and numbers cell first then vertices. */ 1248 for (c_loc = 0, v = 0; c_loc < num_cell_in_set; ++c_loc, ++c) { 1249 for (v_loc = 0; v_loc < num_vertex_per_cell; ++v_loc, ++v) { 1250 cone[v_loc] = cs_connect[v]+numCells-1; 1251 } 1252 if (dim == 3) { 1253 /* Tetrahedra are inverted */ 1254 if (num_vertex_per_cell == 4) { 1255 PetscInt tmp = cone[0]; 1256 cone[0] = cone[1]; 1257 cone[1] = tmp; 1258 } 1259 /* Hexahedra are inverted */ 1260 if (num_vertex_per_cell == 8) { 1261 PetscInt tmp = cone[1]; 1262 cone[1] = cone[3]; 1263 cone[3] = tmp; 1264 } 1265 /* Triangular prisms are inverted */ 1266 if (num_vertex_per_cell == 6) { 1267 PetscInt tmp = cone[1]; 1268 cone[1] = cone[2]; 1269 cone[2] = tmp; 1270 } 1271 } 1272 ierr = DMPlexSetCone(*dm, c, cone);CHKERRQ(ierr); 1273 ierr = DMSetLabelValue(*dm, "Cell Sets", c, cs_id[cs]);CHKERRQ(ierr); 1274 } 1275 ierr = PetscFree2(cs_connect,cone);CHKERRQ(ierr); 1276 } 1277 ierr = PetscFree2(cs_id, cs_order);CHKERRQ(ierr); 1278 } 1279 ierr = DMPlexSymmetrize(*dm);CHKERRQ(ierr); 1280 ierr = DMPlexStratify(*dm);CHKERRQ(ierr); 1281 if (interpolate) { 1282 DM idm; 1283 1284 ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr); 1285 ierr = DMDestroy(dm);CHKERRQ(ierr); 1286 *dm = idm; 1287 } 1288 1289 /* Create vertex set label */ 1290 if (!rank && (num_vs > 0)) { 1291 int vs, v; 1292 /* Read from ex_get_node_set_ids() */ 1293 int *vs_id; 1294 /* Read from ex_get_node_set_param() */ 1295 int num_vertex_in_set; 1296 /* Read from ex_get_node_set() */ 1297 int *vs_vertex_list; 1298 1299 /* Get vertex set ids */ 1300 ierr = PetscMalloc1(num_vs, &vs_id);CHKERRQ(ierr); 1301 PetscStackCallStandard(ex_get_ids,(exoid, EX_NODE_SET, vs_id)); 1302 for (vs = 0; vs < num_vs; ++vs) { 1303 PetscStackCallStandard(ex_get_set_param,(exoid, EX_NODE_SET, vs_id[vs], &num_vertex_in_set, NULL)); 1304 ierr = PetscMalloc1(num_vertex_in_set, &vs_vertex_list);CHKERRQ(ierr); 1305 PetscStackCallStandard(ex_get_set,(exoid, EX_NODE_SET, vs_id[vs], vs_vertex_list, NULL)); 1306 for (v = 0; v < num_vertex_in_set; ++v) { 1307 ierr = DMSetLabelValue(*dm, "Vertex Sets", vs_vertex_list[v]+numCells-1, vs_id[vs]);CHKERRQ(ierr); 1308 } 1309 ierr = PetscFree(vs_vertex_list);CHKERRQ(ierr); 1310 } 1311 ierr = PetscFree(vs_id);CHKERRQ(ierr); 1312 } 1313 /* Read coordinates */ 1314 ierr = DMGetCoordinateSection(*dm, &coordSection);CHKERRQ(ierr); 1315 ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr); 1316 ierr = PetscSectionSetFieldComponents(coordSection, 0, dim);CHKERRQ(ierr); 1317 ierr = PetscSectionSetChart(coordSection, numCells, numCells + numVertices);CHKERRQ(ierr); 1318 for (v = numCells; v < numCells+numVertices; ++v) { 1319 ierr = PetscSectionSetDof(coordSection, v, dim);CHKERRQ(ierr); 1320 ierr = PetscSectionSetFieldDof(coordSection, v, 0, dim);CHKERRQ(ierr); 1321 } 1322 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 1323 ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr); 1324 ierr = VecCreate(PETSC_COMM_SELF, &coordinates);CHKERRQ(ierr); 1325 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr); 1326 ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 1327 ierr = VecSetBlockSize(coordinates, dim);CHKERRQ(ierr); 1328 ierr = VecSetType(coordinates,VECSTANDARD);CHKERRQ(ierr); 1329 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 1330 if (!rank) { 1331 PetscReal *x, *y, *z; 1332 1333 ierr = PetscMalloc3(numVertices,&x,numVertices,&y,numVertices,&z);CHKERRQ(ierr); 1334 PetscStackCallStandard(ex_get_coord,(exoid, x, y, z)); 1335 if (dim > 0) { 1336 for (v = 0; v < numVertices; ++v) coords[v*dim+0] = x[v]; 1337 } 1338 if (dim > 1) { 1339 for (v = 0; v < numVertices; ++v) coords[v*dim+1] = y[v]; 1340 } 1341 if (dim > 2) { 1342 for (v = 0; v < numVertices; ++v) coords[v*dim+2] = z[v]; 1343 } 1344 ierr = PetscFree3(x,y,z);CHKERRQ(ierr); 1345 } 1346 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 1347 ierr = DMSetCoordinatesLocal(*dm, coordinates);CHKERRQ(ierr); 1348 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 1349 1350 /* Create side set label */ 1351 if (!rank && interpolate && (num_fs > 0)) { 1352 int fs, f, voff; 1353 /* Read from ex_get_side_set_ids() */ 1354 int *fs_id; 1355 /* Read from ex_get_side_set_param() */ 1356 int num_side_in_set; 1357 /* Read from ex_get_side_set_node_list() */ 1358 int *fs_vertex_count_list, *fs_vertex_list; 1359 /* Read side set labels */ 1360 char fs_name[MAX_STR_LENGTH+1]; 1361 1362 /* Get side set ids */ 1363 ierr = PetscMalloc1(num_fs, &fs_id);CHKERRQ(ierr); 1364 PetscStackCallStandard(ex_get_ids,(exoid, EX_SIDE_SET, fs_id)); 1365 for (fs = 0; fs < num_fs; ++fs) { 1366 PetscStackCallStandard(ex_get_set_param,(exoid, EX_SIDE_SET, fs_id[fs], &num_side_in_set, NULL)); 1367 ierr = PetscMalloc2(num_side_in_set,&fs_vertex_count_list,num_side_in_set*4,&fs_vertex_list);CHKERRQ(ierr); 1368 PetscStackCallStandard(ex_get_side_set_node_list,(exoid, fs_id[fs], fs_vertex_count_list, fs_vertex_list)); 1369 /* Get the specific name associated with this side set ID. */ 1370 int fs_name_err = ex_get_name(exoid, EX_SIDE_SET, fs_id[fs], fs_name); 1371 for (f = 0, voff = 0; f < num_side_in_set; ++f) { 1372 const PetscInt *faces = NULL; 1373 PetscInt faceSize = fs_vertex_count_list[f], numFaces; 1374 PetscInt faceVertices[4], v; 1375 1376 if (faceSize > 4) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "ExodusII side cannot have %d > 4 vertices", faceSize); 1377 for (v = 0; v < faceSize; ++v, ++voff) { 1378 faceVertices[v] = fs_vertex_list[voff]+numCells-1; 1379 } 1380 ierr = DMPlexGetFullJoin(*dm, faceSize, faceVertices, &numFaces, &faces);CHKERRQ(ierr); 1381 if (numFaces != 1) SETERRQ3(comm, PETSC_ERR_ARG_WRONG, "Invalid ExodusII side %d in set %d maps to %d faces", f, fs, numFaces); 1382 ierr = DMSetLabelValue(*dm, "Face Sets", faces[0], fs_id[fs]);CHKERRQ(ierr); 1383 /* Only add the label if one has been detected for this side set. */ 1384 if (!fs_name_err) { 1385 ierr = DMSetLabelValue(*dm, fs_name, faces[0], fs_id[fs]);CHKERRQ(ierr); 1386 } 1387 ierr = DMPlexRestoreJoin(*dm, faceSize, faceVertices, &numFaces, &faces);CHKERRQ(ierr); 1388 } 1389 ierr = PetscFree2(fs_vertex_count_list,fs_vertex_list);CHKERRQ(ierr); 1390 } 1391 ierr = PetscFree(fs_id);CHKERRQ(ierr); 1392 } 1393 #else 1394 SETERRQ(comm, PETSC_ERR_SUP, "This method requires ExodusII support. Reconfigure using --download-exodusii"); 1395 #endif 1396 PetscFunctionReturn(0); 1397 } 1398