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