00001 #include "SundanceBamgMeshReader.hpp"
00002 #include "SundanceOut.hpp"
00003 #include "PlayaExceptions.hpp"
00004 #include "Teuchos_StrUtils.hpp"
00005
00006 using namespace Teuchos;
00007 using namespace Sundance;
00008
00009
00010 BamgMeshReader::BamgMeshReader(const std::string& fname,
00011 const MeshType& meshType, bool bbAttr,
00012 int verbosity,
00013 const MPIComm& comm)
00014 : MeshReaderBase(fname, meshType, verbosity, comm),
00015 nodeFilename_(filename()),
00016 elemFilename_(filename()),
00017 parFilename_(filename()),
00018 meshFilename_(filename()),
00019 bbFilename_(filename()),
00020 bbAttr_(bbAttr)
00021
00022 {
00023
00024
00025 if (nProc() > 1)
00026 {
00027 nodeFilename_ = nodeFilename_ + Teuchos::toString(myRank());
00028 parFilename_ = parFilename_ + Teuchos::toString(myRank());
00029 elemFilename_ = elemFilename_ + Teuchos::toString(myRank());
00030 }
00031
00032
00033 nodeFilename_ = nodeFilename_ + ".node";
00034 elemFilename_ = elemFilename_ + ".ele";
00035 parFilename_ = parFilename_ + ".par";
00036 meshFilename_ = meshFilename_ + ".mesh";
00037 bbFilename_ = bbFilename_ + ".bb";
00038
00039 SUNDANCE_OUT(this->verb() > 1,
00040 "node filename = " << nodeFilename_);
00041
00042 SUNDANCE_OUT(this->verb() > 1,
00043 "elem filename = " << elemFilename_);
00044
00045 SUNDANCE_OUT(this->verb() > 1,
00046 "mesh filename = " << meshFilename_);
00047
00048 SUNDANCE_OUT(this->verb() > 1,
00049 "bb filename = " << bbFilename_);
00050
00051 }
00052 BamgMeshReader::BamgMeshReader(const ParameterList& params)
00053 : MeshReaderBase(params),
00054 nodeFilename_(filename()),
00055 elemFilename_(filename()),
00056 parFilename_(filename()),
00057 meshFilename_(filename())
00058 {
00059
00060
00061 if (nProc() > 1)
00062 {
00063 nodeFilename_ = nodeFilename_ + Teuchos::toString(myRank());
00064 parFilename_ = parFilename_ + Teuchos::toString(myRank());
00065 elemFilename_ = elemFilename_ + Teuchos::toString(myRank());
00066 }
00067
00068
00069 nodeFilename_ = nodeFilename_ + ".node";
00070 elemFilename_ = elemFilename_ + ".ele";
00071 parFilename_ = parFilename_ + ".par";
00072 meshFilename_ = meshFilename_ + ".mesh";
00073
00074 SUNDANCE_OUT(this->verb() > 1,
00075 "node filename = " << nodeFilename_);
00076
00077 SUNDANCE_OUT(this->verb() > 1,
00078 "elem filename = " << elemFilename_);
00079
00080 SUNDANCE_OUT(this->verb() > 1,
00081 "mesh filename = " << meshFilename_);
00082
00083 SUNDANCE_OUT(this->verb() > 1,
00084 "bb filename = " << bbFilename_);
00085
00086 }
00087
00088
00089 Mesh BamgMeshReader::fillMesh() const
00090 {
00091 Mesh mesh;
00092
00093 Array<int> ptGID;
00094 Array<int> ptOwner;
00095 Array<int> cellGID;
00096 Array<int> cellOwner;
00097
00098 readParallelInfo(ptGID, ptOwner, cellGID, cellOwner);
00099
00100
00101
00102
00103
00104 mesh = readMesh(ptGID, ptOwner);
00105
00106
00107
00108 return mesh;
00109 }
00110
00111 void BamgMeshReader::readParallelInfo(Array<int>& ptGID,
00112 Array<int>& ptOwner,
00113 Array<int>& cellGID,
00114 Array<int>& cellOwner) const
00115 {
00116 int nPoints;
00117 int nElems;
00118 std::string line;
00119 Array<string> tokens;
00120 try
00121 {
00122 ptGID.resize(0);
00123 ptOwner.resize(0);
00124 cellGID.resize(0);
00125 cellOwner.resize(0);
00126
00127
00128
00129
00130 if (nProc() > 1)
00131 {
00132 RCP<std::ifstream> parStream
00133 = openFile(parFilename_, "parallel info");
00134
00135
00136
00137
00138 getNextLine(*parStream, line, tokens, '#');
00139
00140 TEUCHOS_TEST_FOR_EXCEPTION(tokens.length() != 2, std::runtime_error,
00141 "TriangleMeshReader::getMesh() expects 2 entries "
00142 "on the first line of .par file. In "
00143 << parFilename_ << " I found \n[" << line << "]\n");
00144
00145 int np = atoi(tokens[1]);
00146 int pid = atoi(tokens[0]);
00147
00148
00149
00150
00151 TEUCHOS_TEST_FOR_EXCEPTION(np != nProc(), std::runtime_error,
00152 "TriangleMeshReader::getMesh() found "
00153 "a mismatch between the current number of processors="
00154 << nProc() <<
00155 "and the number of processors=" << np
00156 << "in the file " << parFilename_);
00157
00158 TEUCHOS_TEST_FOR_EXCEPTION(pid != myRank(), std::runtime_error,
00159 "TriangleMeshReader::getMesh() found "
00160 "a mismatch between the current processor rank="
00161 << myRank() << "and the processor rank="
00162 << pid << " in the file " << parFilename_);
00163
00164
00165 getNextLine(*parStream, line, tokens, '#');
00166
00167 TEUCHOS_TEST_FOR_EXCEPTION(tokens.length() != 1, std::runtime_error,
00168 "TriangleMeshReader::getMesh() requires 1 entry "
00169 "on the second line of .par file. Found line \n["
00170 << line << "]\n in file " << parFilename_);
00171
00172 nPoints = StrUtils::atoi(tokens[0]);
00173
00174
00175 ptGID.resize(nPoints);
00176 ptOwner.resize(nPoints);
00177
00178 for (int i=0; i<nPoints; i++)
00179 {
00180 getNextLine(*parStream, line, tokens, '#');
00181
00182 TEUCHOS_TEST_FOR_EXCEPTION(tokens.length() != 3, std::runtime_error,
00183 "TriangleMeshReader::getMesh() requires 3 "
00184 "entries on each line of the point section in "
00185 "the .par file. Found line \n[" << line
00186 << "]\n in file " << parFilename_);
00187
00188 ptGID[i] = StrUtils::atoi(tokens[1]);
00189 ptOwner[i] = StrUtils::atoi(tokens[2]);
00190 }
00191
00192
00193
00194
00195 getNextLine(*parStream, line, tokens, '#');
00196
00197 TEUCHOS_TEST_FOR_EXCEPTION(tokens.length() != 1, std::runtime_error,
00198 "TriangleMeshReader::getMesh() requires 1 entry "
00199 "on the cell count line of .par file. Found line \n["
00200 << line << "]\n in file " << parFilename_);
00201
00202 nElems = StrUtils::atoi(tokens[0]);
00203
00204 SUNDANCE_OUT(this->verb() > 1,
00205 "read nElems = " << nElems);
00206
00207
00208
00209
00210 cellGID.resize(nElems);
00211 cellOwner.resize(nElems);
00212 for (int i=0; i<nElems; i++)
00213 {
00214 getNextLine(*parStream, line, tokens, '#');
00215
00216 TEUCHOS_TEST_FOR_EXCEPTION(tokens.length() != 3, std::runtime_error,
00217 "TriangleMeshReader::getMesh() requires 3 "
00218 "entries on each line of the element section in "
00219 "the .par file. Found line \n[" << line
00220 << "]\n in file " << parFilename_);
00221
00222 cellGID[i] = StrUtils::atoi(tokens[1]);
00223 cellOwner[i] = StrUtils::atoi(tokens[2]);
00224 }
00225 }
00226
00227
00228 nPoints = ptGID.length();
00229 nElems = cellGID.length();
00230 }
00231 catch(std::exception& e)
00232 {
00233 SUNDANCE_TRACE(e);
00234 }
00235 }
00236
00237
00238
00239
00240 Mesh BamgMeshReader::readMesh(Array<int>& ptGID,
00241 Array<int>& ptOwner) const
00242 {
00243 Mesh mesh;
00244 std::string line;
00245 Array<string> tokens;
00246 int nPoints=0;
00247
00248
00249
00250
00251
00252
00253
00254
00255
00256
00257
00258
00259
00260
00261
00262
00263
00264
00265
00266
00267
00268
00269
00270
00271
00272
00273
00274
00275
00276
00277
00278
00279
00280
00281
00282
00283
00284
00285
00286
00287
00288
00289
00290
00291
00292
00293
00294
00295
00296
00297
00298
00299
00300
00301
00302
00303
00304
00305
00306
00307
00308
00309
00310
00311
00312
00313
00314
00315
00316
00317
00318
00319
00320
00321
00322
00323
00324
00325
00326
00327
00328
00329
00330
00331
00332
00333
00334
00335
00336
00337
00338
00339
00340 std::cerr << "starting to read meshFile" << std::endl;
00341
00342 RCP<std::ifstream> meshStream = openFile(meshFilename_,
00343 "node & elem info");
00344
00345 Array<string> liners = StrUtils::readFile(*meshStream, '#');
00346
00347 std::cerr << "defined liners" << std::endl;
00348
00349
00350 SUNDANCE_OUT(this->verb() > 3,
00351 "reading nodes from " + meshFilename_);
00352
00353
00354 int dimension=0;
00355 int dimensionIndex = 0;
00356 int lineIndex = 0;
00357 int verticesIndex = 0;
00358 int edgesIndex = 0;
00359 int triangleIndex = 0;
00360 int nCells=0;
00361 std::cerr << "DimensionIndex = " << dimensionIndex << std::endl;
00362
00363
00364 int linerssize = liners.length();
00365 std::cerr << "number of lines in liners = " << linerssize << std::endl;
00366
00367 for (int i = lineIndex; i < linerssize; i++)
00368 {
00369 tokens = StrUtils::stringTokenizer(liners[i]);
00370 std::cerr << "test: i = " << i << "tokens[0] = " << tokens[0] << std::endl;
00371 if (i < 20)
00372 {
00373 std::cerr << "i = " << i << ": number of tokens = " << tokens.length()
00374 << std::endl;
00375 for(int j = 0; j < tokens.length();j++)
00376 {
00377 std::cerr << " token " << j << " = " << tokens[j] << std::endl;
00378 std::cerr << " length of token[" << j << "] = "
00379 << tokens[j].length() << std::endl;
00380 }
00381 }
00382 if (tokens[0] == "Dimension")
00383 {
00384 dimensionIndex = i;
00385 std::cerr << "token[0] = Dimension" << std::endl;
00386 break;
00387 }
00388 }
00389 std::cerr << "DimensionIndex = " << dimensionIndex << std::endl;
00390 if (dimensionIndex > 0)
00391 {
00392 tokens = StrUtils::stringTokenizer(liners[dimensionIndex + 1]);
00393 dimension = atoi(tokens[0]);
00394 std::cerr << "dimension = " << dimension << std::endl;
00395 lineIndex = dimensionIndex + 2;
00396 }
00397
00398 for (int i = lineIndex; i < liners.size(); i++)
00399 {
00400 tokens = StrUtils::stringTokenizer(liners[i]);
00401 if (tokens[0] == "Vertices") verticesIndex = i;
00402 if (tokens[0] == "Vertices") break;
00403 }
00404 if (verticesIndex > 0)
00405 {
00406 tokens = StrUtils::stringTokenizer(liners[verticesIndex + 1]);
00407 nPoints = atoi(tokens[0]);
00408 std::cerr << "nPoints = " << nPoints << std::endl;
00409 ptGID.resize(nPoints);
00410 ptOwner.resize(nPoints);
00411 for (int i=0; i<nPoints; i++)
00412 {
00413 ptGID[i] = i;
00414 ptOwner[i] = 0;
00415 }
00416 lineIndex = verticesIndex + 2;
00417
00418 }
00419
00420 Array<int> ptIndices(nPoints);
00421 Array<int> rawIndices(nPoints);
00422
00423 Array<double> velVector(2*nPoints);
00424 int numbbAttr = 0;
00425
00426
00427
00428 if (bbAttr_)
00429 {
00430 RCP<std::ifstream> bbStream = openFile(bbFilename_,
00431 "velocity info");
00432 Array<string> bbliners = StrUtils::readFile(*bbStream, '#');
00433
00434
00435
00436
00437 int bbdimension;
00438 int bbnumSolns;
00439 int bbnumPoints=0;
00440 int bbsolnType;
00441 int bbdimensionIndex = 0;
00442 int bblineIndex = 0;
00443
00444 std::cerr << "bbDimensionIndex = " << bbdimensionIndex << std::endl;
00445
00446 int bblinerssize = bbliners.size();
00447 std::cerr << "number of lines in bbliners = " << bblinerssize << std::endl;
00448
00449 for (int i = bblineIndex; i < bblinerssize; i++)
00450 {
00451 Array<string> bbtokens = StrUtils::stringTokenizer(bbliners[i]);
00452 if(bbtokens.length() > 0)
00453 {
00454 if(bbtokens.length() != 4)
00455 {
00456 std::cerr <<
00457 "warning: .bb header line should have 4 tokens, read "
00458 << bbtokens.length() << std::endl;
00459 }
00460 std::cerr << "bbline " << i << ": read bbtokens " << bbtokens[0]
00461 << " " << bbtokens[1] << " " << bbtokens[2] << " "
00462 << bbtokens[3] << std::endl;
00463 bbdimensionIndex = i;
00464 break;
00465 }
00466 }
00467 bblineIndex = bbdimensionIndex + 1;
00468 std::cerr << "bblineIndex = " << bblineIndex << std::endl;
00469 if (bblineIndex > 0)
00470 {
00471 Array<string> bbtokens =
00472 StrUtils::stringTokenizer(bbliners[bbdimensionIndex]);
00473 bbdimension = StrUtils::atoi(bbtokens[0]);
00474 std::cerr << "bbdimension = " << bbdimension << std::endl;
00475 if (bbdimension != 2)
00476 std::cerr << "Error! bbdimension should be 2" << std::endl;
00477 bbnumSolns = StrUtils::atoi(bbtokens[1]);
00478 numbbAttr = bbnumSolns;
00479 std::cerr << "number of solutions per vertex = " << bbnumSolns << std::endl;
00480
00481
00482 bbnumPoints = StrUtils::atoi(bbtokens[2]);
00483 std::cerr << "number of vertices = " << bbnumPoints << std::endl;
00484 if (bbnumPoints != nPoints)
00485 std::cerr << "Error!! number of bb points != nPoints" << std::endl;
00486 bbsolnType = StrUtils::atoi(bbtokens[3]);
00487 std::cerr << "bbsolution type = " << bbsolnType << std::endl;
00488 }
00489
00490
00491
00492 for (int i = bblineIndex; i < bbnumPoints + bblineIndex; i++)
00493 {
00494 Array<string> bbtokens = StrUtils::stringTokenizer(bbliners[i]);
00495 if (bbtokens.length() == 0)
00496 std::cerr << "warning: encountered a blank soln line" << std::endl;
00497 int ii = i - bblineIndex;
00498
00499 velVector[ii] = StrUtils::atof(bbtokens[0]);
00500 velVector[ii+bbnumPoints] = StrUtils::atof(bbtokens[1]);
00501
00502
00503
00504 if( i < 5 || i > bbnumPoints - 5)
00505 {
00506
00507
00508
00509
00510 std::cerr << "vel0[" << ii << "] = " << velVector[ii] << std::endl;
00511 std::cerr << "vel1[" << ii << "] = " << velVector[ii+bbnumPoints]
00512 << std::endl;
00513 }
00514
00515 }
00516 }
00517
00518
00519
00520 int nAttributes = 0;
00521 if (bbAttr_) nAttributes = numbbAttr;
00522 std::cerr << "nAttributes = " << nAttributes << std::endl;
00523
00524
00525
00526
00527 mesh = createMesh(dimension);
00528 int count=0;
00529 int offset=0;
00530
00531 Array<bool> usedPoint(nPoints);
00532
00533
00534 nodeAttributes()->resize(nPoints);
00535
00536
00537
00538
00539 for (int i = lineIndex; i < liners.size(); i++)
00540
00541 {
00542 tokens = StrUtils::stringTokenizer(liners[i]);
00543 if (tokens.length() > 1)
00544 {
00545 double x = atof(tokens[0]);
00546 double y = atof(tokens[1]);
00547 count = i - lineIndex;
00548 rawIndices[count] = count;
00549 if ((i == lineIndex) || (i == lineIndex + 1) ||
00550 (i > lineIndex + nPoints - 2))
00551 {
00552 std::cerr << "i = " << i << "; node = (" << x << "," << y << ")"
00553 << std::endl;
00554 std::cerr << "count = " << count << std::endl;
00555 }
00556
00557 double z = 0.0;
00558 Point pt;
00559 int ptLabel = 0;
00560
00561 if (dimension==3)
00562 {
00563 z = atof(tokens[3]);
00564 pt = Point(x,y,z);
00565 }
00566 else
00567 {
00568 pt = Point(x,y);
00569 }
00570
00571 ptIndices[count]
00572 = mesh.addVertex(ptGID[count], pt, ptOwner[count], ptLabel);
00573
00574 (*nodeAttributes())[count].resize(nAttributes);
00575 for (int i=0; i<nAttributes; i++)
00576 {
00577
00578 if(i == 0) (*nodeAttributes())[count][i] = velVector[count];
00579 if(i == 1) (*nodeAttributes())[count][i] =
00580 velVector[count + nPoints];
00581 }
00582 }
00583
00584 if (tokens[0] == "Edges")
00585 {
00586 edgesIndex = i;
00587 break;
00588 }
00589 }
00590 std::cerr << "edgesIndex = " << edgesIndex << std::endl;
00591 std::cerr << "count = " << count << "; npoints - 1 = " << nPoints - 1 << std::endl;
00592
00593 if (count != nPoints - 1) std::cout << "error: # of nodes != # of vertices"
00594 << std::endl;
00595 else std::cerr << "successfully read node data" << std::endl;
00596 lineIndex = edgesIndex + 1;
00597
00598
00599
00600 std::cerr << "lineIndex = " << lineIndex << "; triangleIndex = "
00601 << triangleIndex << std::endl;
00602 for (int i = lineIndex; i < liners.size(); i++)
00603 {
00604 tokens = StrUtils::stringTokenizer(liners[i]);
00605 if (tokens[0] == "Triangles") {triangleIndex = i; break;}
00606 if (tokens[0] == "CrackedEdges") {triangleIndex = i+2; break;}
00607
00608 }
00609 std::cerr << "lineIndex = " << lineIndex << "; triangleIndex = "
00610 << triangleIndex << std::endl;
00611 std::cerr << "tokens[0] = " << tokens[0] << std::endl;
00612
00613 lineIndex = triangleIndex + 1;
00614
00615 Array<int> elemGID;
00616 Array<int> elemOwner;
00617
00618 if (triangleIndex > 0)
00619 {
00620 tokens = StrUtils::stringTokenizer(liners[lineIndex]);
00621 std::cerr << "lineIndex = " << lineIndex << "; tokens[0] = " << tokens[0]
00622 << std::endl;
00623 nCells = atoi(tokens[0]);
00624 elemGID.resize(nCells);
00625 elemOwner.resize(nCells);
00626 std::cerr << "nCells = " << nCells << std::endl;
00627 for (int i=0; i<nCells; i++)
00628 {
00629 elemGID[i] = i;
00630 elemOwner[i] = 0;
00631 }
00632 lineIndex = lineIndex + 1;
00633 }
00634
00635
00636
00637
00638
00639
00640
00641
00642
00643
00644
00645
00646
00647
00648
00649
00650
00651
00652
00653
00654
00655
00656
00657
00658
00659
00660
00661
00662
00663
00664
00665
00666
00667
00668
00669
00670
00671
00672
00673
00674
00675
00676
00677
00678
00679
00680
00681
00682
00683
00684
00685
00686
00687
00688
00689
00690
00691
00692
00693
00694
00695
00696
00697
00698
00699
00700
00701
00702
00703
00704
00705
00706
00707
00708
00709
00710
00711
00712
00713
00714
00715
00716
00717
00718
00719
00720
00721
00722
00723
00724
00725
00726
00727
00728
00729
00730
00731
00732
00733 std::cerr << "# of triangles = " << nCells
00734 << "; starting to read triangle data" << std::endl;
00735 SUNDANCE_OUT(this->verb() > 3,
00736 "done reading nodes, ready to read elements from " + meshFilename_);
00737
00738 int nElems = nCells;
00739 int ptsPerElem = dimension + 1;
00740 elemGID.resize(nElems);
00741 elemOwner.resize(nElems);
00742 offset = 1;
00743
00744 for (int i=0; i<nElems; i++)
00745 {
00746 elemGID[i] = i;
00747 elemOwner[i] = 0;
00748 }
00749 nAttributes = 0;
00750
00751 elemAttributes()->resize(nElems);
00752 count = 0;
00753
00754 int dim = mesh.spatialDim();
00755 Array<int> nodes(dim+1);
00756 if (dim != dimension) std::cerr << "ERROR: dim = " << dim << "!= dimension = "
00757 << dimension << std::endl;
00758
00759 std::cerr << "lineIndex = " << lineIndex << std::endl;
00760 std::cerr << "size of liners = " << liners.size() << std::endl;
00761 std::cerr << "lineIndex+nCells-1 = " << lineIndex+nCells-1 << std::endl;
00762 for (int i = lineIndex; i < liners.size();i++)
00763
00764 {
00765 tokens = StrUtils::stringTokenizer(liners[i]);
00766 if (tokens.length() > 1)
00767 {
00768 if (i < lineIndex + 5) std::cerr << "i = " << i << " first triangles: "
00769 << tokens[0] << " " << tokens[1]
00770 << " " << tokens[2] << std::endl;
00771 if (i == lineIndex+nCells-1) std::cerr << "i = " << i
00772 << " last triangle: " << tokens[0]
00773 << " " << tokens[1] << " "
00774 << tokens[2] << std::endl;
00775
00776 for (int d=0; d<=dim; d++)
00777 {
00778
00779
00780
00781 nodes[d] = ptGID[atoi(tokens[d])-offset];
00782
00783
00784 }
00785
00786 count = i - lineIndex;
00787 int elemLabel = 0;
00788 mesh.addElement(elemGID[count], nodes, elemOwner[count], elemLabel);
00789
00790 (*elemAttributes())[count].resize(nAttributes);
00791 for (int i=0; i<nAttributes; i++)
00792 {
00793
00794
00795
00796 (*elemAttributes())[count][i] = atof(tokens[ptsPerElem+i]);
00797
00798
00799 }
00800 }
00801 if (tokens[0] == "SubDomainFromMesh") break;
00802 }
00803
00804 std::cerr << "# of cells read = " << count + 1 << std::endl;
00805 if (count != nCells - 1) std::cerr << "error: # of cells read != nCells" << std::endl;
00806
00807
00808
00809 return mesh;
00810 }
00811
00812
00813
00814
00815
00816
00817
00818
00819
00820
00821
00822
00823
00824
00825
00826
00827
00828
00829
00830
00831
00832
00833
00834
00835
00836
00837
00838
00839
00840
00841
00842
00843
00844
00845
00846
00847
00848
00849
00850
00851
00852
00853
00854
00855
00856
00857
00858
00859
00860
00861
00862
00863
00864
00865
00866
00867
00868
00869
00870
00871
00872
00873
00874
00875
00876
00877
00878
00879
00880
00881
00882
00883
00884
00885
00886
00887
00888
00889
00890
00891
00892
00893
00894
00895
00896
00897
00898
00899
00900
00901
00902
00903
00904
00905
00906
00907