Sundance 2.1

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     /* If we're running in parallel, we'd better have consistent numbers
00137      * of points in the .exo and .pex files. */
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   /* now we can build the mesh */
00147   mesh = createMesh(dim);
00148 
00149 
00150 
00151   /* Read the points */
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   /* add the points to the mesh */
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   /* Set up the element numbering */
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     /* If we're running in parallel, we'd better have consistent numbers
00206      * of elements in the .exo and .pex files. */
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   /* Read the elements for each block */
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 //  Sundance::Map<int,Array<int> > elemVerts;
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     /* get element connectivity */
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 //        elemVerts.put(elemGID[count], orderedVerts); 
00284       }
00285 
00286       int elemLID 
00287         = mesh.addElement(elemGID[count], orderedVerts, 
00288           elemOwner[count], bid);
00289       if (blockIsSimplicial) 
00290       {
00291 //        elemVerts.put(elemLID, orderedVerts); 
00292         permKey[elemLID]=key;
00293       }
00294     }
00295   }
00296 
00297   /* The topology of the mesh is now fully defined and we can call
00298    * topological functions */
00299   mesh.freezeTopology();
00300  
00301 
00302   /* Read the node sets */
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   /* Read the side sets */
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   /* Read the nodal attributes */
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   /* Read the element attributes */
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     /* if we're running in parallel, read the info on processor 
00446      * distribution */
00447     if (nProc() > 1)
00448     {
00449       RCP<std::ifstream> parStream 
00450         = openFile(parFilename_, "parallel info");
00451      
00452       /* read the number of processors and the processor rank in 
00453        * the file. These must be consistent with the current number of
00454        * processors and the current rank */
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       /* check consistency with the current number of
00466        * processors and the current rank */
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       /* read the number of points */
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       /* read the global ID and the owner PID for each point */
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       /* Read the number of elements */
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       /* read the global ID and the owner PID for each element */
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 }
SundanceExodusMeshReader.cpp
