00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031 #include "SundanceVTKWriter.hpp"
00032 #include "PlayaExceptions.hpp"
00033 #include "SundanceOut.hpp"
00034 #include "PlayaTabs.hpp"
00035 #include "Teuchos_XMLObject.hpp"
00036
00037
00038
00039 using namespace Sundance;
00040 using namespace Sundance;
00041 using namespace Sundance;
00042 using namespace Teuchos;
00043
00044
00045
00046
00047 void VTKWriter::write() const
00048 {
00049 lowLevelWrite(filename(), false);
00050 if (nProc() > 1 && myRank()==0) lowLevelWrite(filename(), true);
00051 }
00052
00053 void VTKWriter::lowLevelWrite(const std::string& filename, bool isPHeader) const
00054 {
00055 std::string PHeader = "";
00056 if (isPHeader) PHeader="P";
00057
00058 std::string f = filename;
00059
00060 if (isPHeader) f = f + ".pvtu";
00061 else if (nProc() > 1)
00062 {
00063 f = f + Teuchos::toString(myRank()) + ".vtu";
00064 }
00065 else
00066 {
00067 f = f + ".vtu";
00068 }
00069
00070 SUNDANCE_VERB_MEDIUM("writing VTK file " << f);
00071
00072 std::ofstream os(f.c_str());
00073
00074 XMLObject head("VTKFile");
00075 head.addAttribute("type", PHeader + "UnstructuredGrid");
00076 head.addAttribute("version", "0.1");
00077
00078 os << head.header() << std::endl;
00079
00080 for (int i=0; i<comments().length(); i++)
00081 {
00082 os << "<!-- " << comments()[i] << " -->" << std::endl;
00083 }
00084
00085 XMLObject ug(PHeader + "UnstructuredGrid");
00086 os << ug.header() << std::endl;
00087
00088 if (isPHeader)
00089 {
00090 writePoints(os, isPHeader);
00091 writePointData(os, isPHeader);
00092 writeCellData(os, isPHeader);
00093 for (int p=0; p<nProc(); p++)
00094 {
00095 XMLObject pc("Piece");
00096 std::string pfile = filename + Teuchos::toString(p) + ".vtu";
00097 size_t pos = pfile.find_last_of("/");
00098 if (pos != std::string::npos)
00099 {
00100 pfile = pfile.substr(pos+1);
00101 }
00102 pc.addAttribute("Source", pfile);
00103 os << pc << std::endl;
00104 }
00105 }
00106 else
00107 {
00108 XMLObject pc("Piece");
00109 pc.addAttribute("NumberOfPoints", Teuchos::toString(mesh().numCells(0)));
00110 pc.addAttribute("NumberOfCells", Teuchos::toString(mesh().numCells(mesh().spatialDim())));
00111
00112 os << pc.header() << std::endl;
00113
00114 writePoints(os, false);
00115 writeCells(os);
00116 writePointData(os, false);
00117 writeCellData(os, false);
00118
00119 os << pc.footer() << std::endl;
00120 }
00121
00122 os << ug.footer() << std::endl;
00123 os << head.footer() << std::endl;
00124 }
00125
00126 void VTKWriter::writePoints(std::ostream& os, bool isPHeader) const
00127 {
00128 std::string PHeader = "";
00129 if (isPHeader) PHeader="P";
00130 XMLObject pts(PHeader + "Points");
00131
00132 os << pts.header() << std::endl;
00133
00134 XMLObject xml(PHeader + "DataArray");
00135 xml.addAttribute("NumberOfComponents", "3");
00136 xml.addAttribute("type", "Float32");
00137 xml.addAttribute("format", "ascii");
00138
00139 os << xml.header() << std::endl;
00140
00141
00142 if (!isPHeader)
00143 {
00144 int np = mesh().numCells(0);
00145 int dim = mesh().spatialDim();
00146
00147 for (int i=0; i<np; i++)
00148 {
00149 const Point& x = mesh().nodePosition(i);
00150
00151 for (int d=0; d<dim; d++)
00152 {
00153 os << x[d] << " ";
00154 }
00155 for (int d=dim; d<3; d++)
00156 {
00157 os << "0.0 ";
00158 }
00159 os << std::endl;
00160 }
00161 }
00162
00163 os << xml.footer() << std::endl;
00164
00165 os << pts.footer() << std::endl;
00166 }
00167
00168
00169 void VTKWriter::writeCells(std::ostream& os) const
00170 {
00171 XMLObject cells("Cells");
00172 os << cells.header() << std::endl;
00173
00174 XMLObject conn("DataArray");
00175 conn.addAttribute("type", "Int32");
00176 conn.addAttribute("Name", "connectivity");
00177 conn.addAttribute("format", "ascii");
00178
00179 int dim = mesh().spatialDim();
00180 int nc = mesh().numCells(dim);
00181 int dummySign;
00182
00183 os << conn.header() << std::endl;
00184 CellType cellType = mesh().cellType(dim);
00185
00186 for (int c=0; c<nc; c++)
00187 {
00188 int nNodes = mesh().numFacets( dim , c , 0 );
00189
00190 switch(cellType)
00191 {
00192 case TriangleCell: case LineCell: case TetCell:
00193 for (int i=0; i<nNodes; i++)
00194 {
00195 os << " " << mesh().facetLID(dim,c,0,i,dummySign);
00196 }
00197 os << std::endl;
00198 break;
00199 case QuadCell:
00200 os << " " << mesh().facetLID(dim,c,0,0,dummySign);
00201 os << " " << mesh().facetLID(dim,c,0,1,dummySign);
00202 os << " " << mesh().facetLID(dim,c,0,3,dummySign);
00203 os << " " << mesh().facetLID(dim,c,0,2,dummySign);
00204 os << std::endl;
00205 break;
00206 case BrickCell:
00207 os << " " << mesh().facetLID(dim,c,0,0,dummySign);
00208 os << " " << mesh().facetLID(dim,c,0,1,dummySign);
00209 os << " " << mesh().facetLID(dim,c,0,2,dummySign);
00210 os << " " << mesh().facetLID(dim,c,0,3,dummySign);
00211 os << " " << mesh().facetLID(dim,c,0,4,dummySign);
00212 os << " " << mesh().facetLID(dim,c,0,5,dummySign);
00213 os << " " << mesh().facetLID(dim,c,0,6,dummySign);
00214 os << " " << mesh().facetLID(dim,c,0,7,dummySign);
00215 os << std::endl;
00216 break;
00217 default:
00218 TEUCHOS_TEST_FOR_EXCEPTION(true, std::runtime_error, "call type " << cellType <<
00219 " not handled in VTKWriter::writeCells()");
00220 }
00221 }
00222
00223 os << conn.footer() << std::endl;
00224
00225
00226 XMLObject offsets("DataArray");
00227 offsets.addAttribute("type", "Int32");
00228 offsets.addAttribute("Name", "offsets");
00229 offsets.addAttribute("format", "ascii");
00230
00231 os << offsets.header() << std::endl;
00232
00233 int count = 0;
00234 for (int c=0; c<nc; c++)
00235 {
00236 count += mesh().numFacets(dim, c, 0);
00237 os << count << std::endl;
00238 }
00239
00240 os << offsets.footer() << std::endl;
00241
00242 XMLObject types("DataArray");
00243 types.addAttribute("type", "UInt8");
00244 types.addAttribute("Name", "types");
00245 types.addAttribute("format", "ascii");
00246
00247 os << types.header() << std::endl;
00248
00249 for (int c=0; c<nc; c++)
00250 {
00251
00252 int vtkCode = 0;
00253 switch(cellType)
00254 {
00255 case TriangleCell:
00256 vtkCode = 5;
00257 break;
00258 case QuadCell:
00259 vtkCode = 9;
00260 break;
00261 case TetCell:
00262 vtkCode = 10;
00263 break;
00264 case BrickCell:
00265 vtkCode = 11;
00266 break;
00267 default:
00268 TEUCHOS_TEST_FOR_EXCEPTION(true, std::runtime_error,
00269 "call type " << cellType << " not handled "
00270 "in VTKWriter::writeCells()");
00271 }
00272 os << vtkCode << std::endl;
00273 }
00274
00275 os << types.footer() << std::endl;
00276
00277 os << cells.footer() << std::endl;
00278 }
00279
00280 void VTKWriter::writePointData(std::ostream& os, bool isPHeader) const
00281 {
00282 std::string PHeader = "";
00283 if (isPHeader) PHeader="P";
00284
00285 XMLObject xml(PHeader + "PointData");
00286
00287 if (pointVectorNames().length() > 0) xml.addAttribute("Vectors", pointVectorNames()[0]);
00288 if (pointScalarNames().length() > 0) xml.addAttribute("Scalars", pointScalarNames()[0]);
00289
00290 os << xml.header() << std::endl;
00291
00292 for (int i=0; i<pointScalarNames().length(); i++)
00293 {
00294 writeDataArray(os, pointScalarNames()[i], pointScalarFields()[i], isPHeader, true);
00295 }
00296
00297 for (int i=0; i<pointVectorNames().length(); i++)
00298 {
00299 writeDataArray(os, pointVectorNames()[i], pointVectorFields()[i], isPHeader, true);
00300 }
00301
00302 os << xml.footer() << std::endl;
00303 }
00304
00305 void VTKWriter::writeCellData(std::ostream& os, bool isPHeader) const
00306 {
00307 std::string PHeader = "";
00308 if (isPHeader) PHeader="P";
00309
00310 XMLObject xml(PHeader + "CellData");
00311
00312 if (cellVectorNames().length() > 0) xml.addAttribute("Vectors", cellVectorNames()[0]);
00313 if (cellScalarNames().length() > 0) xml.addAttribute("Scalars", cellScalarNames()[0]);
00314
00315 os << xml.header() << std::endl;
00316
00317 for (int i=0; i<cellScalarNames().length(); i++)
00318 {
00319 writeDataArray(os, cellScalarNames()[i], cellScalarFields()[i], isPHeader, false);
00320 }
00321
00322
00323 for (int i=0; i<cellVectorNames().length(); i++)
00324 {
00325 writeDataArray(os, cellVectorNames()[i], cellVectorFields()[i], isPHeader, false);
00326 }
00327
00328 os << xml.footer() << std::endl;
00329 }
00330
00331
00332 void VTKWriter::writeDataArray(std::ostream& os, const std::string& name,
00333 const RCP<FieldBase>& expr, bool isPHeader, bool isPointData) const
00334 {
00335 std::string PHeader = "";
00336 if (isPHeader) PHeader="P";
00337
00338 XMLObject xml(PHeader + "DataArray");
00339 xml.addAttribute("type", "Float32");
00340 xml.addAttribute("Name", name);
00341 xml.addAttribute("format", "ascii");
00342
00343 if (expr->numElems() > 1)
00344 {
00345
00346 xml.addAttribute("NumberOfComponents", "3" );
00347 }
00348
00349 os << xml.header() << std::endl;
00350
00351
00352 if (!isPHeader)
00353 {
00354
00355 if (isPointData)
00356 {
00357 int numNodes = mesh().numCells(0);
00358
00359 for (int i=0; i<numNodes; i++)
00360 {
00361 for (int j=0; j < expr->numElems(); j++)
00362 {
00363 if (expr->isDefined(0,i,j)){
00364 double val = expr->getData(0, i, j);
00365 val = (fabs(val) > 1e-16) ? val : 0.0;
00366 os << (float) val << std::endl;
00367 }else
00368 os << undefinedValue() << std::endl;
00369 }
00370
00371 if (expr->numElems() > 1)
00372 for (int j= expr->numElems(); j < 3 ; j++)
00373 os << "0.0 " << std::endl;
00374 }
00375 }
00376 else
00377 {
00378 int dim = mesh().spatialDim();
00379 int nc = mesh().numCells(dim);
00380
00381 for (int c=0; c<nc; c++)
00382 {
00383 for (int j=0; j < expr->numElems(); j++)
00384 {
00385 if (expr->isDefined(dim,c,j)){
00386 double val = expr->getData(dim, c, j);
00387 val = (fabs(val) > 1e-16) ? val : 0.0;
00388 os << (float) val << std::endl;
00389 }else
00390 os << undefinedValue() << std::endl;
00391 }
00392
00393 if (expr->numElems() > 1)
00394 for (int j= expr->numElems(); j < 3 ; j++)
00395 os << "0.0 " << std::endl;
00396 }
00397 }
00398 }
00399
00400 os << xml.footer() << std::endl;
00401 }
00402