00001 #include "SundanceExodusNetCDFMeshReader.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 ExodusNetCDFMeshReader::ExodusNetCDFMeshReader(const std::string& fname,
00013 const MeshType& meshType,
00014 int verbosity,
00015 const MPIComm& comm)
00016 : MeshReaderBase(fname, meshType, verbosity, comm)
00017 {
00018 TEUCHOS_TEST_FOR_EXCEPTION(nProc() > 1, std::runtime_error,
00019 "ExodusNetCDFMeshReader not useable with parallel meshes");
00020 }
00021
00022 ExodusNetCDFMeshReader::ExodusNetCDFMeshReader(const ParameterList& params)
00023 : MeshReaderBase(params)
00024 {
00025 TEUCHOS_TEST_FOR_EXCEPTION(nProc() > 1, std::runtime_error,
00026 "ExodusNetCDFMeshReader not useable with parallel meshes");
00027 }
00028
00029
00030
00031
00032 Mesh ExodusNetCDFMeshReader::fillMesh() const
00033 {
00034 Mesh mesh;
00035
00036 RCP<std::ifstream> is = openFile(filename(), "NetCDF");
00037
00038 std::string line;
00039 Array<string> tokens;
00040
00041
00042 getNextLine(*is, line, tokens, '#');
00043
00044 TEUCHOS_TEST_FOR_EXCEPTION(tokens[0] != "netcdf", std::runtime_error,
00045 "ExodusNetCDF reader expected to find [netcdf] as first "
00046 "token, found " << tokens[0]);
00047
00048
00049 getNextLine(*is, line, tokens, '#');
00050
00051
00052 TEUCHOS_TEST_FOR_EXCEPTION(tokens[0] != "dimensions:", std::runtime_error,
00053 "ExodusNetCDF reader expected to find [dimension:] as first "
00054 "token, found " << tokens[0]);
00055
00056 int nElem = 0;
00057 int nNodes = 0;
00058 int nElemBlocks = 0;
00059 int nSideSets = 0;
00060 int nNodeSets = 0;
00061 int dimension = 0 ;
00062 Array<int> blockSizes;
00063 Array<int> sideSetSizes;
00064 Array<int> nodeSetSizes;
00065 Array<int> nodesPerElem;
00066 while (true)
00067 {
00068 getNextLine(*is, line, tokens, '#');
00069
00070 if (tokens[0] == "variables:") break;
00071 std::string keyword = tokens[0];
00072 std::string equals = tokens[1];
00073 TEUCHOS_TEST_FOR_EXCEPTION(equals!="=", std::runtime_error, "ExodusNetCDF reader "
00074 "expected [=] as second token, found "
00075 << equals);
00076
00077 TEUCHOS_TEST_FOR_EXCEPTION(tokens.size() < 4, std::runtime_error,
00078 "ExodusNetCDF reader found a dimension line with "
00079 "fewer than 4 tokens");
00080
00081 int val = atoi(tokens[2]);
00082 if (keyword=="num_dim")
00083 {
00084 dimension = val;
00085 }
00086 else if (keyword=="num_nodes")
00087 {
00088 nNodes = val;
00089 }
00090 else if (keyword=="num_elem")
00091 {
00092 nElem = val;
00093 }
00094 else if (keyword=="num_el_blk")
00095 {
00096 nElemBlocks = val;
00097 blockSizes.resize(nElemBlocks);
00098 nodesPerElem.resize(nElemBlocks);
00099 }
00100 else if (keyword=="num_side_sets")
00101 {
00102 nSideSets = val;
00103 sideSetSizes.resize(nSideSets);
00104 std::cerr << "num side sets = " << nSideSets << std::endl;
00105 }
00106 else if (keyword=="num_node_sets")
00107 {
00108 nNodeSets = val;
00109 nodeSetSizes.resize(nNodeSets);
00110 std::cerr << "num node sets = " << nNodeSets << std::endl;
00111 }
00112 else
00113 {
00114 for (int i=0; i<nElemBlocks; i++)
00115 {
00116 if (keyword=="num_el_in_blk" + Teuchos::toString(i+1))
00117 {
00118 blockSizes[i] = val;
00119 break;
00120 }
00121 if (keyword=="num_nod_per_el" + Teuchos::toString(i+1))
00122 {
00123 nodesPerElem[i] = val;
00124 break;
00125 }
00126 }
00127 for (int i=0; i<nSideSets; i++)
00128 {
00129 if (keyword=="num_side_ss" + Teuchos::toString(i+1))
00130 {
00131 sideSetSizes[i] = val;
00132 break;
00133 }
00134 }
00135 for (int i=0; i<nNodeSets; i++)
00136 {
00137 if (keyword=="num_nod_ns" + Teuchos::toString(i+1))
00138 {
00139 nodeSetSizes[i] = val;
00140 break;
00141 }
00142 }
00143 }
00144 }
00145
00146
00147
00148
00149 while (true)
00150 {
00151 getNextLine(*is, line, tokens, '#');
00152
00153 if (tokens[0] == "data:") break;
00154 }
00155
00156
00157
00158
00159 Array<double> coords;
00160 coords.reserve(nNodes*dimension);
00161 Array<Array<int> > connect(nElemBlocks);
00162 Array<Array<int> > sideSetElems(nSideSets);
00163 Array<Array<int> > sideSetFacets(nSideSets);
00164 Array<Array<int> > nodeSetNodes(nNodeSets);
00165
00166 bool doneWithData = false;
00167 while(!doneWithData)
00168 {
00169 if (!getNextLine(*is, line, tokens, '#'))
00170 {
00171 doneWithData = true;
00172 break;
00173 }
00174
00175 if (tokens.size()==0)
00176 {
00177
00178 doneWithData=true;
00179 break;
00180 }
00181
00182 if (tokens[0]=="coord")
00183 {
00184 bool done = false;
00185 for (int i=1; i<tokens.size(); i++)
00186 {
00187 if (tokens[i] == "=") continue;
00188 if (tokens[i] == ";")
00189 {
00190 done = true;
00191 break;
00192 }
00193 coords.append(atof(StrUtils::before(tokens[i], ",")));
00194 }
00195 while (!done)
00196 {
00197 if (!getNextLine(*is, line, tokens, '#'))
00198 {
00199 doneWithData = true;
00200 break;
00201 }
00202
00203 for (int i=0; i<tokens.size(); i++)
00204 {
00205 if (tokens[i] == "=") continue;
00206 if (tokens[i] == ";")
00207 {
00208 done = true;
00209 break;
00210 }
00211 coords.append(atof(StrUtils::before(tokens[i], ",")));
00212 }
00213 }
00214 continue;
00215 }
00216
00217 for (int b=0; b<nElemBlocks; b++)
00218 {
00219 if (tokens[0] == "connect" + Teuchos::toString(b+1))
00220 {
00221 connect[b].reserve(blockSizes[b]);
00222 bool done = false;
00223 for (int i=1; i<tokens.size(); i++)
00224 {
00225 if (tokens[i] == "=") continue;
00226 if (tokens[i] == ";")
00227 {
00228 done = true;
00229 break;
00230 }
00231 connect[b].append(atoi(StrUtils::before(tokens[i], ",")));
00232 }
00233 while (!done)
00234 {
00235 if (!getNextLine(*is, line, tokens, '#'))
00236 {
00237 doneWithData = true;
00238 break;
00239 }
00240
00241 for (int i=0; i<tokens.size(); i++)
00242 {
00243 if (tokens[i] == "=") continue;
00244 if (tokens[i] == ";")
00245 {
00246 done = true;
00247 break;
00248 }
00249 connect[b].append(atoi(StrUtils::before(tokens[i], ",")));
00250 }
00251 }
00252 continue;
00253 }
00254 }
00255
00256
00257 for (int s=0; s<nSideSets; s++)
00258 {
00259 if (tokens[0] == "elem_ss" + Teuchos::toString(s+1))
00260 {
00261 sideSetElems[s].reserve(sideSetSizes[s]);
00262 bool done = false;
00263 for (int i=1; i<tokens.size(); i++)
00264 {
00265 if (tokens[i] == "=") continue;
00266 if (tokens[i] == ";")
00267 {
00268 done = true;
00269 break;
00270 }
00271 sideSetElems[s].append(atoi(StrUtils::before(tokens[i], ",")));
00272 }
00273 while (!done)
00274 {
00275 if (!getNextLine(*is, line, tokens, '#'))
00276 {
00277 doneWithData = true;
00278 break;
00279 }
00280
00281 for (int i=0; i<tokens.size(); i++)
00282 {
00283 if (tokens[i] == "=") continue;
00284 if (tokens[i] == ";")
00285 {
00286 done = true;
00287 break;
00288 }
00289 sideSetElems[s].append(atoi(StrUtils::before(tokens[i], ",")));
00290 }
00291 }
00292 continue;
00293 }
00294 }
00295
00296
00297 for (int s=0; s<nSideSets; s++)
00298 {
00299 if (tokens[0] == "side_ss" + Teuchos::toString(s+1))
00300 {
00301 sideSetFacets[s].reserve(sideSetSizes[s]);
00302 bool done = false;
00303 for (int i=1; i<tokens.size(); i++)
00304 {
00305 if (tokens[i] == "=") continue;
00306 if (tokens[i] == ";")
00307 {
00308 done = true;
00309 break;
00310 }
00311 sideSetFacets[s].append(atoi(StrUtils::before(tokens[i], ",")));
00312 }
00313 while (!done)
00314 {
00315 if (!getNextLine(*is, line, tokens, '#'))
00316 {
00317 doneWithData = true;
00318 break;
00319 }
00320 for (int i=0; i<tokens.size(); i++)
00321 {
00322 if (tokens[i] == "=") continue;
00323 if (tokens[i] == ";")
00324 {
00325 done = true;
00326 break;
00327 }
00328 sideSetFacets[s].append(atoi(StrUtils::before(tokens[i], ",")));
00329 }
00330 }
00331 continue;
00332 }
00333 }
00334
00335
00336 for (int s=0; s<nNodeSets; s++)
00337 {
00338 if (tokens[0] == "node_ns" + Teuchos::toString(s+1))
00339 {
00340 nodeSetNodes[s].reserve(nodeSetSizes[s]);
00341 bool done = false;
00342 for (int i=1; i<tokens.size(); i++)
00343 {
00344 if (tokens[i] == "=") continue;
00345 if (tokens[i] == ";")
00346 {
00347 done = true;
00348 break;
00349 }
00350 nodeSetNodes[s].append(atoi(StrUtils::before(tokens[i], ",")));
00351 }
00352 while (!done)
00353 {
00354 if (!getNextLine(*is, line, tokens, '#'))
00355 {
00356 doneWithData = true;
00357 break;
00358 }
00359 for (int i=0; i<tokens.size(); i++)
00360 {
00361 if (tokens[i] == "=") continue;
00362 if (tokens[i] == ";")
00363 {
00364 done = true;
00365 break;
00366 }
00367 nodeSetNodes[s].append(atoi(StrUtils::before(tokens[i], ",")));
00368 }
00369 }
00370 continue;
00371 }
00372 }
00373 }
00374
00375
00376
00377
00378 mesh = createMesh(dimension);
00379
00380
00381 TEUCHOS_TEST_FOR_EXCEPTION(dimension * nNodes != coords.size(), std::runtime_error,
00382 "bad coordinate array in exodus reader");
00383
00384 for (int b=0; b<nElemBlocks; b++)
00385 {
00386 TEUCHOS_TEST_FOR_EXCEPTION(blockSizes[b]*nodesPerElem[b] != connect[b].size(), std::runtime_error,
00387 "bad connectivity array for block " << b << " in exodus reader");
00388 }
00389
00390 for (int s=0; s<nSideSets; s++)
00391 {
00392 TEUCHOS_TEST_FOR_EXCEPTION(sideSetElems[s].size() != sideSetFacets[s].size(), std::runtime_error,
00393 "inconsistent side set specification for ss=" << s
00394 << " in exodus reader ");
00395 }
00396
00397
00398 for (int n=0; n<nNodes; n++)
00399 {
00400 Point x;
00401 if (dimension==2)
00402 {
00403 x = Point(coords[n], coords[nNodes+n]);
00404 }
00405 else
00406 {
00407 x = Point(coords[n], coords[nNodes+n], coords[2*nNodes+n]);
00408 }
00409 mesh.addVertex(n, x, 0, 0);
00410 }
00411
00412
00413
00414 int lid=0;
00415 for (int b=0; b<nElemBlocks; b++)
00416 {
00417 int n = 0;
00418 for (int e=0; e<blockSizes[b]; e++, n+=nodesPerElem[b], lid++)
00419 {
00420 if (dimension==2)
00421 {
00422 mesh.addElement(lid, tuple(connect[b][n]-1, connect[b][n+1]-1, connect[b][n+2]-1), 0, b+1);
00423 }
00424 else
00425 {
00426 mesh.addElement(lid,
00427 tuple(connect[b][n]-1, connect[b][n+1]-1, connect[b][n+2]-1, connect[b][n+3]-1),
00428 0, b+1);
00429 }
00430 }
00431 }
00432
00433 mesh.freezeTopology();
00434
00435
00436
00437 for (int s=0; s<nSideSets; s++)
00438 {
00439 for (int n=0; n<sideSetSizes[s]; n++)
00440 {
00441 int elemID = sideSetElems[s][n];
00442 int facetNum = sideSetFacets[s][n];
00443 int fsign;
00444 int sideLID = mesh.facetLID(dimension, elemID-1, dimension-1,
00445 facetNum-1,fsign);
00446 mesh.setLabel(dimension-1, sideLID, s+1);
00447 }
00448 }
00449
00450
00451 for (int s=0; s<nNodeSets; s++)
00452 {
00453 for (int n=0; n<nodeSetSizes[s]; n++)
00454 {
00455 int nodeNum = nodeSetNodes[s][n]-1;
00456 mesh.setLabel(0, nodeNum, s+1);
00457 }
00458 }
00459
00460 return mesh;
00461 }
00462