Sundance 2.1

00001 #include "SundanceTriangleMeshReader.hpp"
00002 #include "SundanceOut.hpp"
00003 #include "PlayaExceptions.hpp"
00004 
00005 using namespace Sundance;
00006 using namespace Sundance;
00007 
00008 using namespace Teuchos;
00009 using namespace Sundance;
00010 
00011 
00012 TriangleMeshReader::TriangleMeshReader(const std::string& fname,
00013   const MeshType& meshType,
00014   int verbosity,
00015   const MPIComm& comm)
00016   : MeshReaderBase(fname, meshType, verbosity, comm),
00017     nodeFilename_(filename()),
00018     elemFilename_(filename()),
00019     parFilename_(filename()),
00020     sideFilename_(filename()),
00021     offset_(0)
00022 {
00023   
00024 
00025   if (nProc() > 1)
00026     {
00027       std::string suffix =  "." + Teuchos::toString(nProc()) 
00028         + "." + Teuchos::toString(myRank());
00029       nodeFilename_ = nodeFilename_ + suffix;
00030       parFilename_ = parFilename_ + suffix;
00031       elemFilename_ = elemFilename_ + suffix;
00032       sideFilename_ = sideFilename_ + suffix;
00033     }
00034   nodeFilename_ = nodeFilename_ + ".node";
00035   elemFilename_ = elemFilename_ + ".ele";
00036   parFilename_ = parFilename_ + ".par";
00037   sideFilename_ = sideFilename_ + ".side";
00038   
00039   //  verbosity() = 5;
00040   SUNDANCE_OUT(this->verb() > 1,
00041                "node filename = " << nodeFilename_);
00042   
00043   SUNDANCE_OUT(this->verb() > 1,
00044                "elem filename = " << elemFilename_);
00045   
00046 }
00047 TriangleMeshReader::TriangleMeshReader(const ParameterList& params)
00048   : MeshReaderBase(params),
00049     nodeFilename_(filename()),
00050     elemFilename_(filename()),
00051     parFilename_(filename())
00052 {
00053   
00054 
00055   if (nProc() > 1)
00056     {
00057       std::string suffix =  "." + Teuchos::toString(nProc()) 
00058         + "." + Teuchos::toString(myRank());
00059       nodeFilename_ = nodeFilename_ + suffix;
00060       parFilename_ = parFilename_ + suffix;
00061       elemFilename_ = elemFilename_ + suffix;
00062     }
00063   nodeFilename_ = nodeFilename_ + ".node";
00064   elemFilename_ = elemFilename_ + ".ele";
00065   parFilename_ = parFilename_ + ".par";
00066   
00067   SUNDANCE_OUT(this->verb() > 1,
00068                "node filename = " << nodeFilename_);
00069   
00070   SUNDANCE_OUT(this->verb() > 1,
00071                "elem filename = " << elemFilename_);
00072   
00073 }
00074 
00075 
00076 Mesh TriangleMeshReader::fillMesh() const 
00077 {
00078   Mesh mesh;
00079 
00080   Array<int> ptGID;
00081   Array<int> ptOwner;
00082   Array<int> cellGID;
00083   Array<int> cellOwner;
00084 
00085   readParallelInfo(ptGID, ptOwner, cellGID, cellOwner);
00086 
00087   mesh = readNodes(ptGID, ptOwner);
00088 
00089   readElems(mesh, ptGID, cellGID, cellOwner);
00090 
00091   mesh.freezeTopology();
00092 
00093   readSides(mesh);
00094 
00095   return mesh;
00096 }
00097 
00098 void TriangleMeshReader::readParallelInfo(Array<int>& ptGID, 
00099                                           Array<int>& ptOwner,
00100                                           Array<int>& cellGID, 
00101                                           Array<int>& cellOwner) const
00102 {
00103   int nPoints;
00104   int nElems;
00105   std::string line;
00106   Array<string> tokens;
00107   try
00108     {
00109       ptGID.resize(0);
00110       ptOwner.resize(0);
00111       cellGID.resize(0);
00112       cellOwner.resize(0);
00113 
00114       /* if we're running in parallel, read the info on processor 
00115        * distribution */
00116       if (nProc() > 1)
00117         {
00118           RCP<std::ifstream> parStream 
00119             = openFile(parFilename_, "parallel info");
00120      
00121           /* read the number of processors and the processor rank in 
00122            * the file. These must be consistent with the current number of
00123            * processors and the current rank */
00124           getNextLine(*parStream, line, tokens, '#');
00125       
00126           TEUCHOS_TEST_FOR_EXCEPTION(tokens.length() != 2, std::runtime_error,
00127                              "TriangleMeshReader::getMesh() expects 2 entries "
00128                              "on the first line of .par file. In " 
00129                              << parFilename_ << " I found \n[" << line << "]\n");
00130 
00131           int np = atoi(tokens[1]);
00132           int pid = atoi(tokens[0]);
00133 
00134           /* check consistency with the current number of
00135            * processors and the current rank */
00136       
00137           TEUCHOS_TEST_FOR_EXCEPTION(np != nProc(), std::runtime_error,
00138                              "TriangleMeshReader::getMesh() found "
00139                              "a mismatch between the current number of processors="
00140                              << nProc() << 
00141                              "and the number of processors=" << np
00142                              << "in the file " << parFilename_);
00143 
00144           TEUCHOS_TEST_FOR_EXCEPTION(pid != myRank(), std::runtime_error,
00145                              "TriangleMeshReader::getMesh() found "
00146                              "a mismatch between the current processor rank="
00147                              << myRank() << "and the processor rank="
00148                              << pid << " in the file " << parFilename_);
00149 
00150           /* read the number of points */
00151           getNextLine(*parStream, line, tokens, '#');
00152 
00153           TEUCHOS_TEST_FOR_EXCEPTION(tokens.length() != 1, std::runtime_error,
00154                              "TriangleMeshReader::getMesh() requires 1 entry "
00155                              "on the second line of .par file. Found line \n[" 
00156                              << line << "]\n in file " << parFilename_);
00157       
00158           nPoints = StrUtils::atoi(tokens[0]);
00159 
00160           /* read the global ID and the owner PID for each point */
00161           ptGID.resize(nPoints);
00162           ptOwner.resize(nPoints);
00163 
00164           for (int i=0; i<nPoints; i++)
00165             {
00166               getNextLine(*parStream, line, tokens, '#');
00167 
00168               TEUCHOS_TEST_FOR_EXCEPTION(tokens.length() != 3, std::runtime_error,
00169                                  "TriangleMeshReader::getMesh() requires 3 "
00170                                  "entries on each line of the point section in "
00171                                  "the .par file. Found line \n[" << line
00172                                  << "]\n in file " << parFilename_);
00173 
00174               ptGID[i] = StrUtils::atoi(tokens[1]);
00175               ptOwner[i] = StrUtils::atoi(tokens[2]);
00176             }
00177 
00178 
00179           /* Read the number of elements */
00180 
00181           getNextLine(*parStream, line, tokens, '#');
00182 
00183           TEUCHOS_TEST_FOR_EXCEPTION(tokens.length() != 1, std::runtime_error,
00184                              "TriangleMeshReader::getMesh() requires 1 entry "
00185                              "on the cell count line of .par file. Found line \n[" 
00186                              << line << "]\n in file " << parFilename_);
00187 
00188           nElems = StrUtils::atoi(tokens[0]);
00189 
00190           SUNDANCE_OUT(this->verb() > 1,
00191                        "read nElems = " << nElems);
00192 
00193 
00194           /* read the global ID and the owner PID for each element */
00195 
00196           cellGID.resize(nElems);
00197           cellOwner.resize(nElems);
00198           for (int i=0; i<nElems; i++)
00199             {
00200               getNextLine(*parStream, line, tokens, '#');
00201 
00202               TEUCHOS_TEST_FOR_EXCEPTION(tokens.length() != 3, std::runtime_error,
00203                                  "TriangleMeshReader::getMesh() requires 3 "
00204                                  "entries on each line of the element section in "
00205                                  "the .par file. Found line \n[" << line
00206                                  << "]\n in file " << parFilename_);
00207 
00208               cellGID[i] = StrUtils::atoi(tokens[1]);
00209               cellOwner[i] = StrUtils::atoi(tokens[2]);
00210             }
00211         }
00212 
00213       nPoints = ptGID.length();
00214       nElems = cellGID.length();
00215     }
00216   catch(std::exception& e)
00217     {
00218       SUNDANCE_TRACE(e);
00219     }
00220 }
00221 
00222 Mesh TriangleMeshReader::readNodes(Array<int>& ptGID,
00223                                    Array<int>& ptOwner) const 
00224 {
00225   Mesh mesh;
00226   std::string line;
00227   Array<string> tokens;
00228   int nPoints = -1;
00229 
00230   /* Open the node file so we can read in the nodes */
00231   
00232   RCP<std::ifstream> nodeStream = openFile(nodeFilename_, "node info");
00233   
00234   /* read the header line */
00235   getNextLine(*nodeStream, line, tokens, '#');
00236   TEUCHOS_TEST_FOR_EXCEPTION(tokens.length() != 4, std::runtime_error,
00237                      "TriangleMeshReader::getMesh() requires 4 "
00238                      "entries on the header line in "
00239                      "the .node file. Found line \n[" << line
00240                      << "]\n in file " << nodeFilename_);
00241   string headerLine = line;
00242   SUNDANCE_OUT(this->verb() > 2,
00243                "read point header " << line);
00244   
00245   
00246   if (nProc()==1)
00247     {
00248       nPoints = atoi(tokens[0]);
00249       ptGID.resize(nPoints);
00250       ptOwner.resize(nPoints);
00251       for (int i=0; i<nPoints; i++)
00252         {
00253           ptGID[i] = i;
00254           ptOwner[i] = 0;
00255         }
00256     }
00257   else
00258     {
00259       /* If we're running in parallel, we'd better have consistent numbers
00260        * of points in the .node and .par file. */
00261       nPoints = ptGID.length();
00262       TEUCHOS_TEST_FOR_EXCEPTION(atoi(tokens[0]) != nPoints, std::runtime_error,
00263                          "TriangleMeshReader::getMesh() found inconsistent "
00264                          "numbers of points in .node file and par file. Node "
00265                          "file " << nodeFilename_ << " had nPoints=" 
00266                          << atoi(tokens[0]) << " but .par file " 
00267                          << parFilename_ << " had nPoints=" << nPoints);
00268     }
00269 
00270   SUNDANCE_OUT(this->verb() > 3,
00271                "expecting to read " << nPoints << " points");
00272   
00273   int dimension = atoi(tokens[1]);
00274   int nAttributes = atoi(tokens[2]);
00275   int nBdryMarkers = atoi(tokens[3]);
00276 
00277   /* now that we know the dimension, we can build the mesh object */
00278   mesh = createMesh(dimension);
00279 
00280   /* size point-related arrays */
00281   Array<int> ptIndices(nPoints);
00282   Array<bool> usedPoint(nPoints);
00283   nodeAttributes()->resize(nAttributes);
00284   for (int i=0; i<nAttributes; i++)
00285   {
00286     (*nodeAttributes())[i].resize(nPoints);
00287   }
00288   offset_=-1;
00289   bool first = true;
00290 
00291   /* read all the points */
00292   for (int count=0; count<nPoints; count++)
00293     {
00294       getNextLine(*nodeStream, line, tokens, '#');
00295       
00296       TEUCHOS_TEST_FOR_EXCEPTION(tokens.length() 
00297                          != (1 + dimension + nAttributes + nBdryMarkers),
00298                          std::runtime_error,
00299                          "TriangleMeshReader::getMesh() found bad node input "
00300                          "line. Expected " 
00301                          << (1 + dimension + nAttributes + nBdryMarkers)
00302                          << " entries but found line \n[" << 
00303                          line << "]\n in file " << nodeFilename_);
00304       /* Triangle files can use either 0-offset or 1-offset numbering. We'll
00305        * inspect the first node line to decide which numbering to use. */
00306       if (first)
00307         {
00308           offset_ = atoi(tokens[0]);
00309           TEUCHOS_TEST_FOR_EXCEPTION(offset_ < 0 || offset_ > 1, std::runtime_error,
00310                              "TriangleMeshReader::getMesh() expected "
00311                              "either 0-offset or 1-offset numbering. Found an "
00312                              "initial offset of " << offset_ << " in line \n["
00313                              << line << "]\n of file " << nodeFilename_);
00314           first = false;
00315         }
00316       
00317       /* now we can add the point to the mesh */
00318       double x = atof(tokens[1]); 
00319       double y = atof(tokens[2]);
00320       double z = 0.0;
00321       Point pt;
00322       int ptLabel = 0;
00323 
00324       if (dimension==3)
00325         {
00326           z = atof(tokens[3]);
00327           pt = Point(x,y,z);
00328         }
00329       else 
00330         {
00331           pt = Point(x,y);
00332         }
00333 
00334       ptIndices[count] 
00335         = mesh.addVertex(ptGID[count], pt, ptOwner[count], ptLabel);
00336 
00337       for (int i=0; i<nAttributes; i++)
00338         {
00339           (*nodeAttributes())[i][count] = atof(tokens[dimension+1+i]);
00340         }
00341     }
00342   return mesh;
00343 }
00344 
00345 void TriangleMeshReader::readElems(Mesh& mesh,
00346                                    const Array<int>& ptGID,
00347                                    Array<int>& elemGID,
00348                                    Array<int>& elemOwner) const 
00349 {
00350   try
00351     {
00352       std::string line;  
00353       Array<string> tokens;
00354       /* Open the element file */
00355   
00356       RCP<std::ifstream> elemStream = openFile(elemFilename_, "element info");
00357 
00358       getNextLine(*elemStream, line, tokens, '#');
00359 
00360       TEUCHOS_TEST_FOR_EXCEPTION(tokens.length() != 3, std::runtime_error,
00361                          "TriangleMeshReader::getMesh() requires 3 "
00362                          "entries on the header line in "
00363                          "the .ele file. Found line \n[" << line
00364                          << "]\n in file " << elemFilename_);
00365                    
00366       int nElems = -1;
00367 
00368       if (nProc()==1)
00369         {
00370           nElems = atoi(tokens[0]);
00371           elemGID.resize(nElems);
00372           elemOwner.resize(nElems);
00373           for (int i=0; i<nElems; i++)
00374             {
00375               elemGID[i] = i;
00376               elemOwner[i] = 0;
00377             }
00378         }
00379       else
00380         {
00381           /* If we're running in parallel, we'd better have consistent numbers
00382            * of points in the .node and .par file. */
00383           nElems = elemGID.length();
00384           TEUCHOS_TEST_FOR_EXCEPTION(atoi(tokens[0]) != nElems, std::runtime_error,
00385                              "TriangleMeshReader::readElems() found inconsistent "
00386                              "numbers of elements in .ele file and par file. Elem "
00387                              "file " << elemFilename_ << " had nElems=" 
00388                              << atoi(tokens[0]) << " but .par file " 
00389                              << parFilename_ << " had nElems=" << nElems);
00390         }
00391 
00392       int ptsPerElem = atoi(tokens[1]);
00393 
00394       TEUCHOS_TEST_FOR_EXCEPTION(ptsPerElem != mesh.spatialDim()+1, std::runtime_error,
00395                          "TriangleMeshReader::readElems() found inconsistency "
00396                          "between number of points per element=" << ptsPerElem 
00397                          << " and dimension=" << mesh.spatialDim() << ". Number of pts "
00398                          "per element should be dimension + 1");
00399 
00400       int nAttributes = atoi(tokens[2]);
00401       elemAttributes()->resize(nElems);
00402 
00403       int dim = mesh.spatialDim();
00404       Array<int> nodes(dim+1);
00405 
00406       for (int count=0; count<nElems; count++)
00407         {
00408           getNextLine(*elemStream, line, tokens, '#');
00409       
00410           TEUCHOS_TEST_FOR_EXCEPTION(tokens.length() 
00411                              != (1 + ptsPerElem + nAttributes),
00412                              std::runtime_error,
00413                              "TriangleMeshReader::readElems() found bad elem "
00414                              "input line. Expected " 
00415                              << (1 + ptsPerElem + nAttributes)
00416                              << " entries but found line \n[" << 
00417                              line << "]\n in file " << elemFilename_);
00418 
00419           for (int d=0; d<=dim; d++)
00420             {
00421               nodes[d] = ptGID[atoi(tokens[d+1])-offset_];
00422             }
00423 
00424           int elemLabel = 0;
00425           try
00426             {
00427               mesh.addElement(elemGID[count], nodes, elemOwner[count], elemLabel);
00428             }
00429           catch(std::exception& ex1)
00430             {
00431               SUNDANCE_TRACE(ex1);
00432             }
00433       
00434           (*elemAttributes())[count].resize(nAttributes);
00435           for (int i=0; i<nAttributes; i++)
00436             {
00437               (*elemAttributes())[count][i] = atof(tokens[1+ptsPerElem+i]);
00438             }
00439         }
00440     }
00441   catch(std::exception& ex)
00442     {
00443       SUNDANCE_TRACE(ex);
00444     }
00445 }
00446 
00447 
00448 void TriangleMeshReader::readSides(Mesh& mesh) const 
00449 {
00450   try
00451     {
00452       std::string line;  
00453       Array<string> tokens;
00454       /* Open the side file */
00455       RCP<std::ifstream> sideStream;
00456       bool fileOK = false;
00457 
00458       try
00459         {
00460           sideStream = openFile(sideFilename_, "side info");
00461           fileOK = true;
00462         }
00463       catch(std::exception& e) {;}
00464 
00465       /* Not all meshes will have sides files.
00466        * If the sides file doesn't exist, return. */
00467       if (!fileOK)
00468         {
00469           SUNDANCE_VERB_LOW("side file [" << sideFilename_ << "] not found");
00470           return;
00471         }
00472 
00473       getNextLine(*sideStream, line, tokens, '#');
00474 
00475       TEUCHOS_TEST_FOR_EXCEPTION(tokens.length() != 1, std::runtime_error,
00476                          "TriangleMeshReader::readSides() requires 1 "
00477                          "entry on the header line in "
00478                          "the .side file. Found line \n[" << line
00479                          << "]\n in file " << sideFilename_);
00480 
00481       int nSides = atoi(tokens[0]);
00482 
00483       int elemDim = mesh.spatialDim();
00484       int sideDim = elemDim - 1;
00485 
00486       for (int i=0; i<nSides; i++)
00487         {
00488           getNextLine(*sideStream, line, tokens, '#');
00489       
00490           TEUCHOS_TEST_FOR_EXCEPTION(tokens.length() != 4,
00491                              std::runtime_error,
00492                              "TriangleMeshReader::readSides() found bad side "
00493                              "input line. Expected 4 entries but found line \n[" 
00494                              << line << "]\n in file " << sideFilename_);
00495 
00496           int elemGID = atoi(tokens[1]);
00497           int elemFacet = atoi(tokens[2]);
00498           TEUCHOS_TEST_FOR_EXCEPTION(!mesh.hasGID(elemDim, elemGID), std::runtime_error,
00499                              "element GID " << elemGID << " not found");
00500           int elemLID = mesh.mapGIDToLID(elemDim, elemGID);
00501           int o=0; // dummy orientation variable; not needed here
00502           int sideLID = mesh.facetLID(elemDim, elemLID, sideDim, 
00503             elemFacet, o);
00504           
00505           int sideLabel = atoi(tokens[3]);
00506 
00507           mesh.setLabel(sideDim, sideLID, sideLabel);
00508         }
00509     }
00510   catch(std::exception& ex)
00511     {
00512       SUNDANCE_TRACE(ex);
00513     }
00514 }
SundanceTriangleMeshReader.cpp
