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
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
00115
00116 if (nProc() > 1)
00117 {
00118 RCP<std::ifstream> parStream
00119 = openFile(parFilename_, "parallel info");
00120
00121
00122
00123
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
00135
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
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
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
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
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
00231
00232 RCP<std::ifstream> nodeStream = openFile(nodeFilename_, "node info");
00233
00234
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
00260
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
00278 mesh = createMesh(dimension);
00279
00280
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
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
00305
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
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
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
00382
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
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
00466
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;
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 }