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