00001 #include "SundanceExodusMeshReader.hpp"
00002 #include "SundanceVertexSort.hpp"
00003 #include "SundanceOut.hpp"
00004 #include "PlayaExceptions.hpp"
00005 #include "SundancePathUtils.hpp"
00006 #include "Teuchos_Time.hpp"
00007 #include "Teuchos_TimeMonitor.hpp"
00008
00009 #ifdef HAVE_SUNDANCE_EXODUS
00010 #include "exodusII.h"
00011 #endif
00012
00013 using namespace Teuchos;
00014 using namespace Sundance;
00015 using namespace std;
00016
00017 static Time& getExoTimer()
00018 {
00019 static RCP<Time> rtn
00020 = TimeMonitor::getNewTimer("raw exodus reader functions");
00021 return *rtn;
00022 }
00023
00024
00025 static Time& getExoFillTimer()
00026 {
00027 static RCP<Time> rtn
00028 = TimeMonitor::getNewTimer("exodus reader fillMesh");
00029 return *rtn;
00030 }
00031
00032
00033 ExodusMeshReader::ExodusMeshReader(const std::string& fname,
00034 const MeshType& meshType,
00035 int verbosity,
00036 const MPIComm& comm)
00037 : MeshReaderBase(fname, meshType, verbosity, comm),
00038 exoFilename_(fname),
00039 parFilename_(fname)
00040 {
00041 Tabs tab0(0);
00042
00043 if (nProc() > 1)
00044 {
00045 std::string suffix = "-" + Teuchos::toString(nProc())
00046 + "-" + Teuchos::toString(myRank());
00047 exoFilename_ = exoFilename_ + suffix;
00048 parFilename_ = parFilename_ + suffix;
00049 }
00050 exoFilename_ = exoFilename_ + ".exo";
00051 parFilename_ = parFilename_ + ".pxo";
00052
00053 PLAYA_MSG1(verb(), "exodus filename = " << exoFilename_);
00054
00055 if (nProc() > 1)
00056 {
00057 SUNDANCE_OUT(this->verb() > 1,
00058 "par filename = " << parFilename_);
00059 }
00060 }
00061
00062
00063
00064
00065
00066 Mesh ExodusMeshReader::fillMesh() const
00067 {
00068 TimeMonitor fillTimer(getExoFillTimer());
00069 Tabs tab0(0);
00070 Tabs tab1;
00071 Mesh mesh;
00072 #ifndef HAVE_SUNDANCE_EXODUS
00073 TEUCHOS_TEST_FOR_EXCEPTION(true, std::runtime_error,
00074 "ExodusMeshReader called for a build without ExodusII");
00075 #else
00076 PLAYA_ROOT_MSG1(verb(), tab0 << "in ExodusMeshReader::fillMesh()");
00077
00078 int CPU_word_size = 8;
00079 int IO_word_size = 0;
00080 float version;
00081
00082 Array<int> ptGID;
00083 Array<int> ptOwner;
00084 Array<int> elemGID;
00085 Array<int> elemOwner;
00086
00087 readParallelInfo(ptGID, ptOwner, elemGID, elemOwner);
00088
00089 if (verb() > 2) ex_opts(EX_DEBUG | EX_VERBOSE);
00090
00091 PLAYA_MSG3(verb(), tab1 << "opening file");
00092 std::string resolvedName = searchForFile(exoFilename_);
00093 int exoID = ex_open(resolvedName.c_str(), EX_READ,
00094 &CPU_word_size, &IO_word_size, &version);
00095
00096 TEUCHOS_TEST_FOR_EXCEPTION(exoID < 0, std::runtime_error, "ExodusMeshReader unable to "
00097 "open file: " << exoFilename_);
00098
00099 PLAYA_MSG3(verb(), tab1 << "exoID=" << exoID);
00100
00101 TEUCHOS_TEST_FOR_EXCEPT(IO_word_size != 8 || CPU_word_size != 8);
00102
00103 char title[MAX_LINE_LENGTH+1];
00104
00105 int dim = 0;
00106 int numNodes = 0;
00107 int numElems = 0;
00108 int numElemBlocks = 0;
00109 int numNodeSets = 0;
00110 int numSideSets = 0;
00111
00112 int ierr = ex_get_init(exoID, title, &dim, &numNodes, &numElems,
00113 &numElemBlocks, &numNodeSets, &numSideSets);
00114
00115 TEUCHOS_TEST_FOR_EXCEPTION(numNodes <= 0, std::runtime_error, "invalid numNodes="
00116 << numNodes);
00117 TEUCHOS_TEST_FOR_EXCEPTION(numElems <= 0, std::runtime_error, "invalid numElems="
00118 << numElems);
00119
00120 PLAYA_MSG2(verb(), tab1 << "numNodes=" << numNodes);
00121 PLAYA_MSG2(verb(), tab1 << "numElems=" << numElems);
00122
00123
00124 if (nProc()==1)
00125 {
00126 ptGID.resize(numNodes);
00127 ptOwner.resize(numNodes);
00128 for (int i=0; i<numNodes; i++)
00129 {
00130 ptGID[i] = i;
00131 ptOwner[i] = 0;
00132 }
00133 }
00134 else
00135 {
00136
00137
00138 TEUCHOS_TEST_FOR_EXCEPTION((int)ptGID.size() != numNodes, std::runtime_error,
00139 "ExodusMeshReader::getMesh() found inconsistent "
00140 "numbers of points in .exo file and .pex files. Exo "
00141 "file " << exoFilename_ << " had nPoints="
00142 << numNodes << " but .pex file "
00143 << parFilename_ << " had nPoints=" << ptGID.size());
00144 }
00145
00146
00147 mesh = createMesh(dim);
00148
00149
00150
00151
00152 PLAYA_MSG2(verb(), tab1 << "reading vertices");
00153 Array<double> x(numNodes);
00154 Array<double> y(numNodes);
00155 Array<double> z(numNodes * (dim > 2));
00156
00157 if (dim == 2)
00158 {
00159 ierr = ex_get_coord(exoID, &(x[0]), &(y[0]), (void*) 0);
00160 TEUCHOS_TEST_FOR_EXCEPT(ierr < 0);
00161 TEUCHOS_TEST_FOR_EXCEPT(ierr > 0);
00162 }
00163 else if (dim==3)
00164 {
00165 ierr = ex_get_coord(exoID, (void*) &(x[0]), (void*) &(y[0]), (void*) &(z[0]));
00166 TEUCHOS_TEST_FOR_EXCEPT(ierr < 0);
00167 }
00168 else
00169 {
00170 TEUCHOS_TEST_FOR_EXCEPTION(dim < 2 || dim > 3, std::runtime_error,
00171 "invalid dimension=" << dim << " in ExodusMeshReader");
00172 }
00173
00174
00175
00176 PLAYA_MSG2(verb(), tab1 << "adding vertices to mesh");
00177 for (int n=0; n<numNodes; n++)
00178 {
00179 Point p;
00180 if (dim==2)
00181 {
00182 p = Point(x[n], y[n]);
00183 }
00184 else
00185 {
00186 p = Point(x[n], y[n], z[n]);
00187 }
00188 mesh.addVertex(ptGID[n], p, ptOwner[n], 0);
00189 }
00190
00191
00192
00193 if (nProc()==1)
00194 {
00195 elemGID.resize(numElems);
00196 elemOwner.resize(numElems);
00197 for (int i=0; i<numElems; i++)
00198 {
00199 elemGID[i] = i;
00200 elemOwner[i] = 0;
00201 }
00202 }
00203 else
00204 {
00205
00206
00207 TEUCHOS_TEST_FOR_EXCEPTION((int)elemGID.size() != numElems, std::runtime_error,
00208 "ExodusMeshReader::readElems() found inconsistent "
00209 "numbers of elements in .exo file and .pex files. Exodus "
00210 "file " << exoFilename_ << " had nElems="
00211 << numElems << " but .pex file "
00212 << parFilename_ << " had nElems=" << elemGID.size());
00213 }
00214
00215
00216
00217 PLAYA_MSG2(verb(), tab1 << "reading elements");
00218 Array<int> blockIDs(numElemBlocks);
00219 if (numElemBlocks > 0)
00220 {
00221 TimeMonitor timer(getExoTimer());
00222 ierr = ex_get_elem_blk_ids(exoID, &(blockIDs[0]));
00223 }
00224 TEUCHOS_TEST_FOR_EXCEPT(ierr < 0);
00225 PLAYA_MSG2(verb(), tab1 << "file contains block IDs " << blockIDs);
00226 int count = 0;
00227 Array<int> permKey;
00228 permKey.resize(numElems);
00229
00230 Array<int> blockIsSimplicial(numElemBlocks);
00231 bool allBlocksAreSimplicial = true;
00232
00233 for (int b=0; b<numElemBlocks; b++)
00234 {
00235 Tabs tab2;
00236 char elemType[MAX_LINE_LENGTH+1];
00237 int elsInBlock=0;
00238 int nodesPerEl=0;
00239 int numAttrs=0;
00240 int bid = blockIDs[b];
00241 PLAYA_MSG2(verb(), tab2 << "reading elems for block ID=" << bid);
00242 {
00243 TimeMonitor timer(getExoTimer());
00244 ierr = ex_get_elem_block(exoID, bid, elemType, &elsInBlock,
00245 &nodesPerEl, &numAttrs);
00246 }
00247 TEUCHOS_TEST_FOR_EXCEPTION(ierr < 0, std::runtime_error,
00248 "ex_get_elem_block returned error code ierr=" << ierr);
00249
00250 bool blockIsSimplicial = true;
00251 if (nodesPerEl != dim+1)
00252 {
00253 blockIsSimplicial=false;
00254 allBlocksAreSimplicial=false;
00255 }
00256
00257
00258
00259 PLAYA_MSG2(verb(), tab2 << "reading connectivity for block ID=" << bid);
00260 Array<int> connect(elsInBlock * nodesPerEl);
00261
00262 {
00263 TimeMonitor timer(getExoTimer());
00264 ierr = ex_get_elem_conn(exoID, bid, &(connect[0]));
00265 }
00266 TEUCHOS_TEST_FOR_EXCEPT(ierr < 0);
00267 int n=0;
00268 Array<int> orderedVerts(nodesPerEl);
00269 Array<int> exVerts(nodesPerEl);
00270
00271 for (int e=0; e<elsInBlock; e++, n+=nodesPerEl, count++)
00272 {
00273
00274 for (int v=0; v<nodesPerEl; v++)
00275 {
00276 orderedVerts[v] = ptGID[connect[n+v]-1];
00277 }
00278 exVerts = orderedVerts;
00279 int key = -1;
00280 if (blockIsSimplicial)
00281 {
00282 vertexSort(orderedVerts, &key);
00283
00284 }
00285
00286 int elemLID
00287 = mesh.addElement(elemGID[count], orderedVerts,
00288 elemOwner[count], bid);
00289 if (blockIsSimplicial)
00290 {
00291
00292 permKey[elemLID]=key;
00293 }
00294 }
00295 }
00296
00297
00298
00299 mesh.freezeTopology();
00300
00301
00302
00303 Array<int> nsIDs(numNodeSets);
00304 if (numNodeSets > 0)
00305 {
00306 TimeMonitor timer(getExoTimer());
00307 ierr = ex_get_node_set_ids(exoID, &(nsIDs[0]));
00308 }
00309 TEUCHOS_TEST_FOR_EXCEPT(ierr < 0);
00310 for (int ns=0; ns<numNodeSets; ns++)
00311 {
00312 int nNodes;
00313 int nDist;
00314 int nsID = nsIDs[ns];
00315 {
00316 TimeMonitor timer(getExoTimer());
00317 ierr = ex_get_node_set_param(exoID, nsID, &nNodes, &nDist);
00318 }
00319 TEUCHOS_TEST_FOR_EXCEPT(ierr < 0);
00320 Array<int> nodes(nNodes);
00321 {
00322 TimeMonitor timer(getExoTimer());
00323 ierr = ex_get_node_set(exoID, nsID, &(nodes[0]));
00324 }
00325 TEUCHOS_TEST_FOR_EXCEPT(ierr < 0);
00326 for (int n=0; n<nNodes; n++)
00327 {
00328 mesh.setLabel(0, nodes[n]-1, nsID);
00329 }
00330 }
00331
00332
00333
00334 Array<int> ssIDs(numSideSets);
00335 if (numSideSets > 0)
00336 {
00337 TimeMonitor timer(getExoTimer());
00338 ierr = ex_get_side_set_ids(exoID, &(ssIDs[0]));
00339 }
00340 TEUCHOS_TEST_FOR_EXCEPT(ierr < 0);
00341 for (int ss=0; ss<numSideSets; ss++)
00342 {
00343 int nSides;
00344 int nDist;
00345 int ssID = ssIDs[ss];
00346 {
00347 TimeMonitor timer(getExoTimer());
00348 ierr = ex_get_side_set_param(exoID, ssID, &nSides, &nDist);
00349 }
00350 TEUCHOS_TEST_FOR_EXCEPT(ierr < 0);
00351 Array<int> sides(nSides);
00352 Array<int> elems(nSides);
00353 {
00354 TimeMonitor timer(getExoTimer());
00355 ierr = ex_get_side_set(exoID, ssID, &(elems[0]), &(sides[0]));
00356 }
00357 TEUCHOS_TEST_FOR_EXCEPT(ierr < 0);
00358 for (int n=0; n<nSides; n++)
00359 {
00360 int elemID = elems[n]-1;
00361 int facetNum = sides[n]-1;
00362 int fsign;
00363 if (allBlocksAreSimplicial)
00364 {
00365 int key = permKey[elemID];
00366 facetNum = exFacetIndexToUFCFacetIndex(dim, key, facetNum);
00367 }
00368 int sideLID = mesh.facetLID(dim, elemID, dim-1,
00369 facetNum, fsign);
00370 mesh.setLabel(dim-1, sideLID, ssID);
00371 }
00372 }
00373
00374
00375
00376 int nNodalVars = 0;
00377 {
00378 TimeMonitor timer(getExoTimer());
00379 ierr = ex_get_var_param(exoID, "N", &nNodalVars);
00380 }
00381 TEUCHOS_TEST_FOR_EXCEPT(ierr < 0);
00382
00383 Array<Array<double> >& funcVals = *nodeAttributes();
00384 funcVals.resize(nNodalVars);
00385
00386 for (int i=0; i<nNodalVars; i++)
00387 {
00388 int t = 1;
00389 funcVals[i].resize(mesh.numCells(0));
00390 {
00391 TimeMonitor timer(getExoTimer());
00392 ierr = ex_get_nodal_var(exoID, t, i+1, mesh.numCells(0), &(funcVals[i][0]));
00393 }
00394 TEUCHOS_TEST_FOR_EXCEPT(ierr < 0);
00395 }
00396
00397
00398 int nElemVars = 0;
00399 {
00400 TimeMonitor timer(getExoTimer());
00401 ierr = ex_get_var_param(exoID, "E", &nElemVars);
00402 }
00403 TEUCHOS_TEST_FOR_EXCEPT(ierr < 0);
00404
00405 Array<Array<double> >& eFuncVals = *elemAttributes();
00406 eFuncVals.resize(nElemVars);
00407
00408 for (int i=0; i<nElemVars; i++)
00409 {
00410 int t = 1;
00411 eFuncVals[i].resize(mesh.numCells(mesh.spatialDim()));
00412 {
00413 TimeMonitor timer(getExoTimer());
00414 ierr = ex_get_elem_var(exoID, t, i+1, 1, mesh.numCells(mesh.spatialDim()), &(eFuncVals[i][0]));
00415 }
00416 TEUCHOS_TEST_FOR_EXCEPT(ierr < 0);
00417 }
00418
00419 ierr = ex_close(exoID);
00420 TEUCHOS_TEST_FOR_EXCEPT(ierr < 0);
00421
00422 #endif
00423 return mesh;
00424 }
00425
00426
00427
00428
00429 void ExodusMeshReader::readParallelInfo(Array<int>& ptGID,
00430 Array<int>& ptOwner,
00431 Array<int>& cellGID,
00432 Array<int>& cellOwner) const
00433 {
00434 int nPoints;
00435 int nElems;
00436 std::string line;
00437 Array<string> tokens;
00438 try
00439 {
00440 ptGID.resize(0);
00441 ptOwner.resize(0);
00442 cellGID.resize(0);
00443 cellOwner.resize(0);
00444
00445
00446
00447 if (nProc() > 1)
00448 {
00449 RCP<std::ifstream> parStream
00450 = openFile(parFilename_, "parallel info");
00451
00452
00453
00454
00455 getNextLine(*parStream, line, tokens, '#');
00456
00457 TEUCHOS_TEST_FOR_EXCEPTION(tokens.length() != 2, std::runtime_error,
00458 "ExodusMeshReader::getMesh() expects 2 entries "
00459 "on the first line of .par file. In "
00460 << parFilename_ << " I found \n[" << line << "]\n");
00461
00462 int np = atoi(tokens[1]);
00463 int pid = atoi(tokens[0]);
00464
00465
00466
00467
00468 TEUCHOS_TEST_FOR_EXCEPTION(np != nProc(), std::runtime_error,
00469 "ExodusMeshReader::getMesh() found "
00470 "a mismatch between the current number of processors="
00471 << nProc() <<
00472 "and the number of processors=" << np
00473 << "in the file " << parFilename_);
00474
00475 TEUCHOS_TEST_FOR_EXCEPTION(pid != myRank(), std::runtime_error,
00476 "ExodusMeshReader::getMesh() found "
00477 "a mismatch between the current processor rank="
00478 << myRank() << "and the processor rank="
00479 << pid << " in the file " << parFilename_);
00480
00481
00482 getNextLine(*parStream, line, tokens, '#');
00483
00484 TEUCHOS_TEST_FOR_EXCEPTION(tokens.length() != 1, std::runtime_error,
00485 "ExodusMeshReader::getMesh() requires 1 entry "
00486 "on the second line of .pxo file. Found line \n["
00487 << line << "]\n in file " << parFilename_);
00488
00489 nPoints = StrUtils::atoi(tokens[0]);
00490
00491
00492 ptGID.resize(nPoints);
00493 ptOwner.resize(nPoints);
00494
00495 for (int i=0; i<nPoints; i++)
00496 {
00497 getNextLine(*parStream, line, tokens, '#');
00498
00499 TEUCHOS_TEST_FOR_EXCEPTION(tokens.length() != 3, std::runtime_error,
00500 "ExodusMeshReader::getMesh() requires 3 "
00501 "entries on each line of the point section in "
00502 "the .pxo file. Found line \n[" << line
00503 << "]\n in file " << parFilename_);
00504
00505 ptGID[i] = StrUtils::atoi(tokens[1]);
00506 ptOwner[i] = StrUtils::atoi(tokens[2]);
00507 }
00508
00509
00510
00511
00512 getNextLine(*parStream, line, tokens, '#');
00513
00514 TEUCHOS_TEST_FOR_EXCEPTION(tokens.length() != 1, std::runtime_error,
00515 "ExodusMeshReader::getMesh() requires 1 entry "
00516 "on the cell count line of .pxo file. Found line \n["
00517 << line << "]\n in file " << parFilename_);
00518
00519 nElems = StrUtils::atoi(tokens[0]);
00520
00521 SUNDANCE_OUT(this->verb() > 1,
00522 "read nElems = " << nElems);
00523
00524
00525
00526
00527 cellGID.resize(nElems);
00528 cellOwner.resize(nElems);
00529 for (int i=0; i<nElems; i++)
00530 {
00531 getNextLine(*parStream, line, tokens, '#');
00532
00533 TEUCHOS_TEST_FOR_EXCEPTION(tokens.length() != 3, std::runtime_error,
00534 "ExodusMeshReader::getMesh() requires 3 "
00535 "entries on each line of the element section in "
00536 "the .pxo file. Found line \n[" << line
00537 << "]\n in file " << parFilename_);
00538
00539 cellGID[i] = StrUtils::atoi(tokens[1]);
00540 cellOwner[i] = StrUtils::atoi(tokens[2]);
00541 }
00542 }
00543
00544 nPoints = ptGID.length();
00545 nElems = cellGID.length();
00546 }
00547 catch(std::exception& e)
00548 {
00549 SUNDANCE_TRACE(e);
00550 }
00551 }